我可以使用mi包汇总imputed随机效应模型估计吗？

Question

似乎在过去几年中，`mi`包进行了相当大的重写。以下教程详细介绍了“旧”方法：http://thomasleeper.com/Rcourse/Tutorials/mi.html。采用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（）（或如果我们遵循Leeper示例，则使用coef.mi（））汇总结果。

在当前版本的mi中（我安装了v1.0），这些最后的步骤似乎已合并为一个单一的函数pool（）。 pool（）函数似乎会读取在上述填充过程中分配给变量的family和link函数，然后使用指定的公式估计一个带有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解释为连续变量，因此我使用mi中的change()函数来重新分配到"un"以表示"无序分类变量"，然后按照上述方式进行。

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

# impute
mydf2_imp <- mi(mydf2_mdf)

现在，除非1.0版本的mi已经去掉了之前版本的功能（即使用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. 如果可以，如何操作？

- joemienko

Answer 1

为了提供另一种选择，有一个专注于混合效应模型MI以及汇总结果的包(mitml, 在这里找到它)。

使用该包非常简单。它依赖于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. 

- SimonG

---

我正在使用mitml软件包。我也得到了一个p值为0.000。如何获得更具信息量的p值？ - Outsider

另外，您能否解释一下这里的Intercept~~Intercept|ID是什么意思？ - Outsider

为了获得更高精度的p值和其他结果，您可以保存汇总结果，然后查看其中存储的内部值：result <- testEstimates(fit)，接着是 result$estimates。Intercept~~Intercept|ID 的值是斜率方差的汇总估计值。 - SimonG

Answer 2

您可以在pool()函数中指定FUN参数来更改估计器。在您的情况下，它将是summary(pool(y ~ x1 + x2 + (1|as.numeric(as.character(group_var))), data = mydf2_imp, FUN = lmer))。这可能有效，也可能无效，但它是合法的语法。如果失败了，那么您可以使用complete函数创建完整的数据框，对每个数据框调用lmer，并自己平均结果，例如： dfs <- complete(mydf2_imp) estimates <- lapply(dfs, FUN = lme4, formula = y ~ x1 + x2 + (1|as.numeric(as.character(group_var)))) rowMeans(sapply(estimates, FUN = fixef))

- Ben Goodrich

---

谢谢！@ben-goodrich，如何组合随机效应有什么建议吗？在那里使用平均值是否合适？ - joemienko

当然：rowMeans(sapply(estimates, FUN = function(x) unlist(ranef(x)))) - Ben Goodrich