signal non-convergence if theta_vcov is not positive definite (#468)

NEWS.md

@@ -15,6 +15,7 @@
 - `model.matrix` is updated to ensure that the `NA` values are dropped. Additionally, an argument `use_response` is added to decide whether records with `NA` values in the response should be discarded.
 - `model.frame` is updated to ensure that the `na.action` works.
 - `predict` is updated to allow duplicated subject IDs for unconditional prediction.
+- `mmrm` now checks on the positive definiteness of the covariance matrix `theta_vcov`. If it is not positive definite, non-convergence is messaged appropriately.
 
 # mmrm 0.3.12
 

R/tmb.R

@@ -334,6 +334,12 @@ h_mmrm_tmb_check_conv <- function(tmb_opt, mmrm_tmb) {
     warning("Model convergence problem: theta_vcov contains non-finite values.")
     return()
   }
+  eigen_vals <- eigen(theta_vcov, only.values = TRUE)$values
+  if (mode(eigen_vals) == "complex" || any(eigen_vals <= 0)) {
+    # Note: complex eigen values signal that the matrix is not symmetric, therefore not positive definite.
+    warning("Model convergence problem: theta_vcov is not positive definite.")
+    return()
+  }
   qr_rank <- qr(theta_vcov)$rank
   if (qr_rank < ncol(theta_vcov)) {
     warning("Model convergence problem: theta_vcov is numerically singular.")

tests/testthat/test-fit.R

@@ -419,16 +419,20 @@ test_that("mmrm works for specific small data example", {
 })
 
 test_that("mmrm works for custom optimizer", {
+  custom_optimizer <- optim
   fit <- mmrm(
     FEV1 ~ ARMCD + ar1(AVISIT | SEX / USUBJID),
     data = fev_data,
     reml = TRUE,
-    optimizer_fun = silly_optimizer,
-    optimizer_args = list(value_add = 2, message = "this is wrong"),
+    optimizer_fun = custom_optimizer,
+    optimizer_args = list(method = "Nelder-Mead", hessian = TRUE),
     start = c(1, 2, 3, 4),
     method = "Kenward-Roger"
   )
-  expect_identical(fit$theta_est, c(3, 4, 5, 6))
+  expect_equal(
+    fit$theta_est, c(2.2225, 0.4054, 2.2136, 0.4389),
+    tolerance = 1e-4
+  )
 })
 
 test_that("mmrm works for constructed control", {

tests/testthat/test-tmb.R

@@ -812,6 +812,98 @@ test_that("h_mmrm_tmb_check_conv warns if theta_vcov is singular", {
     list(theta_vcov = chol_m %*% t(chol_m)),
     class = "mmrm_tmb"
   )
+  expect_warning(
+    h_mmrm_tmb_check_conv(tmb_opt, mmrm_tmb),
+    "Model convergence problem: theta_vcov is not positive definite."
+  )
+})
+
+test_that("h_mmrm_tmb_check_conv warns if theta_vcov is not positive definite", {
+  tmb_opt <- list(
+    par = 1:5,
+    objective = 10,
+    convergence = 0,
+    message = NULL
+  )
+  not_positive_definite_matrix <- matrix(
+    c(
+      1, 2, 3, 4, 5,
+      2, 1, 4, 5, 6,
+      3, 4, 1, 6, 7,
+      4, 5, 6, 1, 8,
+      5, 6, 7, 8, -10
+    ),
+    nrow = 5,
+    byrow = TRUE
+  )
+  eigenvalues <- eigen(not_positive_definite_matrix)$values
+  assert_true(any(eigenvalues < 0))
+  mmrm_tmb <- structure(
+    list(theta_vcov = not_positive_definite_matrix),
+    class = "mmrm_tmb"
+  )
+  expect_warning(
+    h_mmrm_tmb_check_conv(tmb_opt, mmrm_tmb),
+    "Model convergence problem: theta_vcov is not positive definite."
+  )
+})
+
+test_that("h_mmrm_tmb_check_conv warns if theta_vcov is not symmetric", {
+  tmb_opt <- list(
+    par = 1:5,
+    objective = 10,
+    convergence = 0,
+    message = NULL
+  )
+  not_symmetric_matrix <- matrix(
+    c(
+      1,  1,  3,  4,  5,
+      2,  1,  4,  5,  6,
+      3,  4,  1,  6,  7,
+      4,  5,  6,  1,  8,
+      5,  6,  7,  8,  9
+    ),
+    nrow = 5,
+    byrow = TRUE
+  )
+  assert_false(all(not_symmetric_matrix == t(not_symmetric_matrix)))
+  mmrm_tmb <- structure(
+    list(theta_vcov = not_symmetric_matrix),
+    class = "mmrm_tmb"
+  )
+  expect_warning(
+    h_mmrm_tmb_check_conv(tmb_opt, mmrm_tmb),
+    "Model convergence problem: theta_vcov is not positive definite."
+  )
+})
+
+test_that("h_mmrm_tmb_check_conv warns if theta_vcov is numerically singular", {
+  tmb_opt <- list(
+    par = 1:5,
+    objective = 10,
+    convergence = 0,
+    message = NULL
+  )
+  singular_matrix <- matrix(
+    c(
+      1, 0.99, 0.98, 0.97, 0.96,
+      0.99, 1, 0.99, 0.98, 0.97,
+      0.98, 0.99, 1, 0.99, 0.98,
+      0.97, 0.98, 0.99, 1, 0.99,
+      0.96, 0.97, 0.98, 0.99, 1
+    ),
+    nrow = 5,
+    byrow = TRUE
+  )
+  # Slightly perturb the matrix to make it numerically singular
+  singular_matrix[5, ] <- singular_matrix[4, ] + 1e-10
+
+  assert_true(all(eigen(singular_matrix)$values > 0))
+  assert_true(qr(singular_matrix)$rank == 4)
+  mmrm_tmb <- structure(
+    list(theta_vcov = singular_matrix),
+    class = "mmrm_tmb"
+  )
   expect_warning(
     h_mmrm_tmb_check_conv(tmb_opt, mmrm_tmb),
     "Model convergence problem: theta_vcov is numerically singular."