Can we have "mmrm" package working with emmeans after MICE imputation using "list" rather than "mira"?
Generalized opened this issue · 2 comments
Generalized commented
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 ( openpharma/mmrm#247 (comment) )
> 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
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
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.