r - 我可以使用 mi 包汇集估算的随机效应模型估计值吗？

Question

看来 mi在过去几年的某个时刻，软件包进行了相当大的重写。

以下教程详细概述了“旧”的做事方式:http://thomasleeper.com/Rcourse/Tutorials/mi.html

"new"的做事方式(坚持使用 Leeper 的模拟演示)看起来像这样:

#load mi
library(mi)
#set seed
set.seed(10)
#simulate some data (with some observations missing)
x1 <- runif(100, 0, 5)
x2 <- rnorm(100)
y <- 2*x1 + 20*x2 + rnorm(100)
mydf <- cbind.data.frame(x1, x2, y)
mydf$x1[sample(1:nrow(mydf), 20, FALSE)] <- NA
mydf$x2[sample(1:nrow(mydf), 10, FALSE)] <- NA

# Convert to a missing_data.frame
mydf_mdf <- missing_data.frame(mydf)

# impute
mydf_imp <- mi(mydf_mdf)

尽管函数名称已更改，但这实际上与“旧”的处理方式非常相似。

最大的变化(从我的角度来看)是替换了以下“旧”功能

lm.mi(formula, mi.object, ...)

glm.mi(formula, mi.object, family = gaussian, ...)

bayesglm.mi(formula, mi.object, family = gaussian, ...)

polr.mi(formula, mi.object, ...)

bayespolr.mi(formula, mi.object, ...)

lmer.mi(formula, mi.object, rescale=FALSE, ...)

glmer.mi(formula, mi.object, family = gaussian, rescale=FALSE, ...) 。

以前，用户可以使用这些函数之一为每个估算数据集计算模型，然后使用 mi.pooled() 汇集结果。 (或者 coef.mi() 如果我们遵循 Leeper 示例)。

在mi的当前版本中(我安装了 v1.0)，最后这些步骤似乎已合并为一个函数，pool() 。 pool()函数似乎读取在上述插补过程中分配给变量的族和链接函数，然后使用 bayesglm 估计模型使用如下所示的指定公式。

# run models on imputed data and pool the results
summary(pool(y ~ x1 + x2, mydf_imp))

## 
## Call:
## pool(formula = y ~ x1 + x2, data = mydf_imp)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.98754  -0.40923   0.03393   0.46734   2.13848  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.34711    0.25979  -1.336    0.215    
## x1           2.07806    0.08738  23.783 1.46e-13 ***
## x2          19.90544    0.11068 179.844  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.7896688)
## 
##     Null deviance: 38594.916  on 99  degrees of freedom
## Residual deviance:    76.598  on 97  degrees of freedom
## AIC: 264.74
## 
## Number of Fisher Scoring iterations: 7

看起来我们即将恢复模拟 Beta 值(2 和 20)。换句话说，它的行为符合预期。

为了获得分组变量，我们采用一组稍微大一些的数据，并简单地模拟随机效应。

mydf2 <- data.frame(x1 = rep(runif(100, 0, 5), 20)
                   ,x2 = rep(rnorm(100, 0, 2.5), 20)
                   ,group_var = rep(1:20, each = 100)
                   ,noise = rep(rnorm(100), 20))

mydf2$y <- 2*mydf2$x1 + 20*mydf2$x2 + mydf2$noise

mydf2$x1[sample(1:nrow(mydf2), 200, FALSE)] <- NA
mydf2$x2[sample(1:nrow(mydf2), 100, FALSE)] <- NA

# Convert to a missing_data.frame
mydf2_mdf <- missing_data.frame(mydf2)

show(mydf2_mdf)

## Object of class missing_data.frame with 2000 observations on 5 variables
## 
## There are 4 missing data patterns
## 
## Append '@patterns' to this missing_data.frame to access the corresponding pattern for every observation or perhaps use table()
## 
##                 type missing method  model
## x1        continuous     200    ppd linear
## x2        continuous     100    ppd linear
## group_var continuous       0   <NA>   <NA>
## noise     continuous       0   <NA>   <NA>
## y         continuous       0   <NA>   <NA>
## 
##             family     link transformation
## x1        gaussian identity    standardize
## x2        gaussian identity    standardize
## group_var     <NA>     <NA>    standardize
## noise         <NA>     <NA>    standardize
## y             <NA>     <NA>    standardize

自 missing_data.frame()似乎正在解释group_var作为连续，我使用 change()来自 mi 的函数重新分配给"un"对于“无序分类”，然后按上述方式进行。

mydf2_mdf <- change(mydf2_mdf, y = "group_var", what = "type", to = "un"  )

# impute
mydf2_imp <- mi(mydf2_mdf)

现在，除非 mi 的 1.0 版本已删除以前版本的功能(即 lmer.mi 和 glmer.mi 可用的功能)，我假设在公式中添加随机效应应该指向 pool()到适当的lme4功能。然而，最初的错误消息表明情况并非如此。

# run models on imputed data and pool the results
summary(pool(y ~ x1 + x2 + (1|group_var), mydf2_imp))
## Warning in Ops.factor(1, group_var): '|' not meaningful for factors
## Warning in Ops.factor(1, group_var): '|' not meaningful for factors
## Error in if (prior.scale[j] < min.prior.scale) {: missing value where TRUE/FALSE needed

按照我的警告消息并从我的因子中提取整数确实可以得到一个估计，但结果表明 pool()仍在估计固定效应模型 bayesglm并保持我尝试的随机效应不变。

summary(pool(y ~ x1 + x2 + (1|as.numeric(as.character(group_var))), mydf2_imp))

## 
## Call:
## pool(formula = y ~ x1 + x2 + (1 | as.numeric(as.character(group_var))), 
##     data = mydf2_imp)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.93633  -0.69923   0.01073   0.56752   2.12167  
## 
## Coefficients:
##                                               Estimate Std. Error  t value
## (Intercept)                                  1.383e-01  2.596e+02    0.001
## x1                                           1.995e+00  1.463e-02  136.288
## x2                                           2.000e+01  8.004e-03 2499.077
## 1 | as.numeric(as.character(group_var))TRUE -3.105e-08  2.596e+02    0.000
##                                             Pr(>|t|)    
## (Intercept)                                        1    
## x1                                            <2e-16 ***
## x2                                            <2e-16 ***
## 1 | as.numeric(as.character(group_var))TRUE        1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.8586836)
## 
##     Null deviance: 5384205.2  on 1999  degrees of freedom
## Residual deviance:    1713.9  on 1996  degrees of freedom
## AIC: 5377
## 
## Number of Fisher Scoring iterations: 4

我的问题是:

1. 是否可以使用mi轻松生成汇总随机效应估计？，和
2. 如果是，怎么做？

Answer 1

只是为了提供一种替代方案，有一个包主要关注混合效应模型的 MI 以及汇集从中获得的结果(mitml，find it here)。

使用该包非常简单。它依赖于 pan 和 jomo 包进行插补，但它也可以处理来自其他 MI 包的输入 (?as.mitml.list) .

混合效应模型的池化估计大部分是自动化的，并包含在 testEstimates 函数中。

require(mitml)
require(lme4)

data(studentratings)

# impute example data using 'pan'
fml <- ReadDis + SES ~ ReadAchiev + (1|ID)
imp <- panImpute(studentratings, formula=fml, n.burn=1000, n.iter=100, m=5)

implist <- mitmlComplete(imp, print=1:5)

# fit model using lme4
fit.lmer <- with(implist, lmer(SES ~ (1|ID)))

# pool results using 'Rubin's rules'
testEstimates(fit.lmer, var.comp=TRUE)

输出:

# Call:

# testEstimates(model = fit.lmer, var.comp = TRUE)

# Final parameter estimates and inferences obtained from 5 imputed data sets.

#              Estimate Std.Error   t.value        df   p.value       RIV       FMI 
# (Intercept)    46.988     1.119    41.997   801.800     0.000     0.076     0.073 

#                         Estimate 
# Intercept~~Intercept|ID   38.272 
# Residual~~Residual       298.446 
# ICC|ID                     0.114 

# Unadjusted hypothesis test as appropriate in larger samples.