fix emmeans on spatial covariance (#450)

* fix emmeans on spatial covariance * [skip style] [skip vbump] Restyle files * [skip roxygen] [skip vbump] Roxygen Man Pages Auto Update * fix lintr issue * update url * revert url * use doi url * Update NEWS.md Co-authored-by: Daniel Sabanes Bove <[email protected]> * use explicit NULL in if * update news --------- Co-authored-by: github-actions <41898282+github-actions[bot]@users.noreply.github.com> Co-authored-by: Daniel Sabanes Bove <[email protected]>
openpharma · Jul 15, 2024 · b0741f8 · b0741f8
1 parent b06da42
commit b0741f8
Show file tree

Hide file tree

Showing 5 changed files with 80 additions and 4 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -96,7 +96,7 @@ Language: en-US
 LazyData: true
 NeedsCompilation: yes
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 Collate:
     'between-within.R'
     'catch-routine-registration.R'

diff --git a/NEWS.md b/NEWS.md
@@ -1,5 +1,11 @@
 # mmrm 0.3.12.9000
 
+### Bug Fixes
+
+- Previously `emmeans` will return `NA` for spatial covariance structure. This is fixed now.
+
+# mmrm 0.3.12
+
 ### New Features
 
 - Add parameter `conditional` for `predict` method to control whether the prediction is conditional on the observation or not.

diff --git a/R/interop-emmeans.R b/R/interop-emmeans.R
@@ -30,8 +30,14 @@ recover_data.mmrm <- function(object, ...) { # nolint
   fun_call <- stats::getCall(object)
   # subject_var is excluded because it should not contain fixed effect.
   # visit_var is not excluded because emmeans can provide marginal mean
-  # by each visit.
-  model_frame <- stats::model.frame(object, include = c("visit_var", "response_var", "group_var"))
+  # by each visit if visit_var is not spatial.
+  model_frame <- stats::model.frame(
+    object,
+    include = c(
+      if (!object$formula_parts$is_spatial) "visit_var" else NULL,
+      "response_var", "group_var"
+    )
+  )
   model_terms <- stats::delete.response(stats::terms(model_frame))
   emmeans::recover_data(
     fun_call,
@@ -62,7 +68,13 @@ emm_basis.mmrm <- function(object, # nolint
     beta_hat[kept] <- component(object, "beta_est")
     orig_model_mat <- stats::model.matrix(
       trms,
-      stats::model.frame(object, include = c("visit_var", "response_var", "group_var")),
+      stats::model.frame(
+        object,
+        include = c(
+          if (!object$formula_parts$is_spatial) "visit_var" else NULL,
+          "response_var", "group_var"
+        )
+      ),
       contrasts.arg = contrasts
     )
     estimability::nonest.basis(orig_model_mat)

diff --git a/tests/testthat/_snaps/emmeans.md b/tests/testthat/_snaps/emmeans.md
@@ -0,0 +1,51 @@
+# emmeans works for spatial covariance
+
+    structure(list(contrast = c("PBO VIS1 - TRT VIS1", "PBO VIS1 - PBO VIS2", 
+    "PBO VIS1 - TRT VIS2", "PBO VIS1 - PBO VIS3", "PBO VIS1 - TRT VIS3", 
+    "PBO VIS1 - PBO VIS4", "PBO VIS1 - TRT VIS4", "TRT VIS1 - PBO VIS2", 
+    "TRT VIS1 - TRT VIS2", "TRT VIS1 - PBO VIS3", "TRT VIS1 - TRT VIS3", 
+    "TRT VIS1 - PBO VIS4", "TRT VIS1 - TRT VIS4", "PBO VIS2 - TRT VIS2", 
+    "PBO VIS2 - PBO VIS3", "PBO VIS2 - TRT VIS3", "PBO VIS2 - PBO VIS4", 
+    "PBO VIS2 - TRT VIS4", "TRT VIS2 - PBO VIS3", "TRT VIS2 - TRT VIS3", 
+    "TRT VIS2 - PBO VIS4", "TRT VIS2 - TRT VIS4", "PBO VIS3 - TRT VIS3", 
+    "PBO VIS3 - PBO VIS4", "PBO VIS3 - TRT VIS4", "TRT VIS3 - PBO VIS4", 
+    "TRT VIS3 - TRT VIS4", "PBO VIS4 - TRT VIS4"), estimate = c(-4.58256155613587, 
+    -5.05924050093087, -9.47063321820258, -10.3804344374154, -13.8559896148128, 
+    -15.5594383984546, -20.2882833923675, -0.476678944795, -4.88807166206672, 
+    -5.79787288127955, -9.27342805867693, -10.9768768423188, -15.7057218362316, 
+    -4.41139271727172, -5.32119393648455, -8.79674911388193, -10.5001978975238, 
+    -15.2290428914366, -0.909801219212832, -4.38535639661021, -6.08880518025206, 
+    -10.8176501741649, -3.47555517739738, -5.17900396103923, -9.90784895495207, 
+    -1.70344878364184, -6.43229377755469, -4.72884499391284), SE = c(1.17792215101562, 
+    1.05934929373609, 1.15256301754827, 1.10405218961219, 1.20738705061085, 
+    1.1510614228206, 1.17360022884585, 1.16893841060977, 1.03673949448934, 
+    1.15990389615065, 1.16564528295435, 1.17897264682249, 1.15768123768126, 
+    1.14311878568268, 1.03621794370689, 1.19809496664043, 1.11233916501112, 
+    1.16436979085874, 1.13409029249035, 1.06473988772671, 1.15466799006841, 
+    1.09802289103113, 1.18971053554843, 1.04231948270098, 1.15547240776813, 
+    1.21050168189854, 1.09557535929194, 1.17549724125026), df = c(515.783958214491, 
+    352.662219824974, 517.932040471532, 446.570856856199, 525.630607576163, 
+    514.176896437678, 513.679726240331, 522.480446555653, 309.76554043098, 
+    522.699920876888, 422.083464913346, 518.311431879116, 510.149541790452, 
+    524.018790363486, 316.433429268676, 527.789686140492, 432.248132345559, 
+    520.926106512673, 524.295212759084, 279.48183426502, 520.688672109129, 
+    417.532827009544, 527.86969791444, 332.816650858056, 521.188696397337, 
+    526.888205312726, 317.888086840647, 516.650695813446), t.ratio = c(-3.89037726490221, 
+    -4.77580013584382, -8.21701987137199, -9.40212295676138, -11.4760131043336, 
+    -13.5174701279861, -17.2872183335542, -0.407787904365588, -4.71485044029734, 
+    -4.99858040008388, -7.95561754016047, -9.31054411813805, -13.566533968961, 
+    -3.85908513841559, -5.13520728800442, -7.34228033571397, -9.4397448438479, 
+    -13.0792150492026, -0.802229968140369, -4.11871147794909, -5.27320860422512, 
+    -9.85193502114173, -2.92134521259428, -4.96872988272154, -8.57471704935794, 
+    -1.4072254579359, -5.87115593920607, -4.0228465265331), p.value = c(0.000113194389980527, 
+    2.62900632607059e-06, 1.68172460557064e-15, 2.80861841173368e-19, 
+    2.32945322197733e-27, 7.65844658326766e-36, 4.10753407600708e-53, 
+    0.68359643524018, 3.66225791754532e-06, 7.89314567399559e-07, 
+    1.65536288122504e-14, 3.59608207718733e-19, 5.10569464455113e-36, 
+    0.000128009414950824, 4.93439343271037e-07, 8.00101365627567e-13, 
+    2.3429997228639e-19, 5.28117117492745e-34, 0.422783412931216, 
+    5.02251675202135e-05, 1.96853134208825e-07, 9.97046669418354e-21, 
+    0.00363441315070957, 1.08001230717931e-06, 1.13734565979031e-16, 
+    0.159950166048537, 1.08898297285039e-08, 6.6089605299358e-05)), class = c("summary_emm", 
+    "data.frame"), row.names = c(NA, 28L), estName = "estimate", pri.vars = "contrast", adjust = "none", side = 0, delta = 0, type = "link")
+
diff --git a/tests/testthat/test-emmeans.R b/tests/testthat/test-emmeans.R
@@ -93,6 +93,13 @@ test_that("emmeans works as expected", {
   expect_equal(result_emmeans, expected_emmeans, tolerance = 1e-3)
 })
 
+test_that("emmeans works for spatial covariance", {
+  fit <- mmrm(FEV1 ~ FEV1_BL + ARMCD * AVISIT + sp_exp(VISITN, VISITN2 | USUBJID), data = fev_data)
+  lsmeans <- expect_silent(emmeans(fit, ~ ARMCD * AVISIT))
+  diffs <- expect_silent(summary(pairs(lsmeans, adjust = "none")))
+  expect_snapshot_tolerance(diffs)
+})
+
 test_that("emmeans works as expected for transformed variables", {
   skip_if_not_installed("emmeans", minimum_version = "1.6")