rvlenth / emmeans

Estimated marginal means
https://rvlenth.github.io/emmeans/
341 stars 30 forks source link

Can we have "mmrm" package working with emmeans after MICE imputation using "list" rather than "mira"? #407

Closed Generalized closed 1 year ago

Generalized commented 1 year ago

Let's use this naive scenario (totally nonsensical, but just to illustrate)

nhanes <- structure(list(age = c(1, 2, 1, 3, 1, 3, 1, 1, 2, 2, 1, 2, 3, 
2, 1, 1, 3, 2, 1, 3, 1, 1, 1, 3, 2), bmi = c(NA, 22.7, NA, NA, 
20.4, NA, 22.5, 30.1, 22, NA, NA, NA, 21.7, 28.7, 29.6, NA, 27.2, 
26.3, 35.3, 25.5, NA, 33.2, 27.5, 24.9, 27.4), hyp = structure(c(NA, 
1L, 1L, NA, 1L, NA, 1L, 1L, 1L, NA, NA, NA, 1L, 2L, 1L, NA, 2L, 
2L, 1L, 2L, NA, 1L, 1L, 1L, 1L), levels = c("1", "2"), class = "factor"), 
    chl = c(NA, 187, 187, NA, 113, 184, 118, 187, 238, NA, NA, 
    NA, 206, 204, NA, NA, 284, 199, 218, NA, NA, 229, 131, NA, 
    186), id = structure(1:25, levels = c("1", "2", "3", "4", 
    "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", 
    "16", "17", "18", "19", "20", "21", "22", "23", "24", "25"
    ), class = "factor"), tim = structure(c(1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L), class = "factor", levels = "1")), row.names = c("1", 
"2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", 
"14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", 
"25"), class = "data.frame")

nhanes_miss <- mice(nhanes, m = 5, method = "pmm", seed = 123)

While "mice" works with the nlme::gls()

> imp_m_gls <- with(nhanes_miss, gls(bmi ~ hyp + chl))
> (gls_pooled <- pool(imp_m_gls))

Class: mipo    m = 5 
         term m   estimate         ubar            b            t dfcom        df       riv    lambda       fmi
1 (Intercept) 5 19.8365753 2.081502e+01 4.9135453308 2.671128e+01    22 13.230280 0.2832692 0.2207403 0.3167657
2        hyp2 5 -0.3542480 5.744374e+00 2.6254225741 8.894882e+00    22  9.270615 0.5484509 0.3541933 0.4594540
3         chl 5  0.0320654 6.100493e-04 0.0001111273 7.434021e-04    22 14.651663 0.2185934 0.1793817 0.2723609

> model_parameters(gls_pooled)
# Fixed Effects

Parameter   | Coefficient |   SE |         95% CI |     t |    df |     p
-------------------------------------------------------------------------
(Intercept) |       19.84 | 5.17 | [ 8.69, 30.98] |  3.84 | 13.23 | 0.002
hyp2        |       -0.35 | 2.98 | [-7.07,  6.36] | -0.12 |  9.27 | 0.908
chl         |        0.03 | 0.03 | [-0.03,  0.09] |  1.18 | 14.65 | 0.258

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald t-distribution approximation.
Warning messages:
1: Could not recover model data from environment. Please make sure your data is available in your workspace.
  Trying to retrieve data from the model frame now. 
2: Could not get model data. 

But it doesn't work with mmrm. I reported this to the authors ( https://github.com/openpharma/mmrm/issues/247#issuecomment-1458476424 )

> imp_m_mmrm <- with(nhanes_miss, mmrm(bmi ~ hyp + chl + us(tim|id)))
Error in attr(data, "dataname") <- toString(match.call()$data) : 
  cannot set attribute on a symbol

But still there's a workaround to make mmrm working with mice:

dat_mice <- complete(nhanes_miss, "all")

fit_mmrm <- function(dat) {
  mod <- mmrm(bmi ~ hyp + chl + us(tim|id),data=dat)
  return(mod)
}

imp_m_mmrm <- lapply(dat_mice, fit_mmrm)

> imp_m_mmrm
$`1`
mmrm fit

Formula:     bmi ~ hyp + chl + us(tim | id)
Data:        dat (used 25 observations from 25 subjects with maximum 1 timepoints)
Covariance:  unstructured (1 variance parameters)
Method:      REML
Deviance:    135.6773

Coefficients: 
(Intercept)        hyp2         chl 
18.27884063  0.90312733  0.03547781 

Model Inference Optimization:
Converged with code 0 and message: convergence: rel_reduction_of_f <= factr*epsmch

$`2`......
.....

Let's compare both results pooled:

> pool(imp_m_mmrm)
Class: mipo    m = 5 
         term m   estimate         ubar            b            t dfcom        df       riv    lambda       fmi
1 (Intercept) 5 19.8365753 2.081502e+01 4.9135453308 2.671128e+01    22 13.230279 0.2832692 0.2207403 0.3167657
2        hyp2 5 -0.3542480 5.744374e+00 2.6254225741 8.894881e+00    22  9.270615 0.5484509 0.3541933 0.4594540
3         chl 5  0.0320654 6.100493e-04 0.0001111273 7.434021e-04    22 14.651663 0.2185934 0.1793818 0.2723609

> pool(imp_m_gls)
Class: mipo    m = 5 
         term m   estimate         ubar            b            t dfcom        df       riv    lambda       fmi
1 (Intercept) 5 19.8365753 2.081502e+01 4.9135453308 2.671128e+01    22 13.230280 0.2832692 0.2207403 0.3167657
2        hyp2 5 -0.3542480 5.744374e+00 2.6254225741 8.894882e+00    22  9.270615 0.5484509 0.3541933 0.4594540
3         chl 5  0.0320654 6.100493e-04 0.0001111273 7.434021e-04    22 14.651663 0.2185934 0.1793817 0.2723609

> summary(pool(imp_m_mmrm))
         term   estimate std.error  statistic        df     p.value
1 (Intercept) 19.8365753 5.1682953  3.8381273 13.230279 0.001992168
2        hyp2 -0.3542480 2.9824287 -0.1187784  9.270615 0.907985019
3         chl  0.0320654 0.0272654  1.1760471 14.651663 0.258331811

> summary(pool(imp_m_gls))
         term   estimate std.error  statistic        df     p.value
1 (Intercept) 19.8365753 5.1682954  3.8381272 13.230280 0.001992168
2        hyp2 -0.3542480 2.9824288 -0.1187784  9.270615 0.907985021
3         chl  0.0320654 0.0272654  1.1760471 14.651663 0.258331829

But problem is that while gls gives proper "mira" object, mmrm - gives an ordinary list, which won't work with emmeans.

> emmeans(imp_m_mmrm, specs = ~hyp+chl)
Error in (function (object, at, cov.reduce = mean, cov.keep = get_emm_option("cov.keep"),  : 
  Can't handle an object of class  “list” 
 Use help("models", package = "emmeans") for information on supported models.

Would it be possible, to allow the same that can be done with mira, only using a list, as above?

rvlenth commented 1 year ago

No, this is not possible. My package is based on object classes and methods. While I suppose it's possible to define methods for class list, it would be inappropriate to do so because a list is a very general thing, whereas what you have here is specialized to multiple imputations.

Generalized commented 1 year ago

I understand the reason. Indeed, it would be dangerous... I will try to "wrap" the list into a "mira" class or wait until the "mmrm" will support the "mice". Thank you for checking this.