我可以使用 mi 包合并估算的随机效应模型估计吗？
Can I pool imputed random effect model estimates using the mi package?
似乎 mi 包在过去几年中的某个时候进行了相当大的重写。
"old" 做事的方式在以下教程中有很好的概述：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)
虽然函数名称发生了变化，但这实际上与 "old" 的处理方式非常相似。
最大的变化（从我的角度来看）是替换了以下 "old" 功能
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() 函数似乎读取在上述插补过程中分配给变量的族和 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" for "unordered categorical" 然后按照上面的方法进行。
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 可用的功能），否则我会假设添加一个随机公式中的 effect 应将 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
我的问题是：
- 是否可以使用
mi 轻松生成合并的随机效应估计值？和
- 如果是，怎么办？
您可以将 FUN 参数指定给 pool() 函数以更改估算器。在您的情况下，它将是 summary(pool(y ~ x1 + x2 + (1|as.numeric(as.character(group_var))), data = mydf2_imp, FUN = lmer))。这可能有效也可能无效，但它是合法的语法。如果失败，那么您可以使用 complete 函数创建完成的 data.frames，对每个函数调用 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))
只是为了提供一个替代方案，有一个包相当侧重于混合效应模型的 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.
似乎 mi 包在过去几年中的某个时候进行了相当大的重写。
"old" 做事的方式在以下教程中有很好的概述：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)
虽然函数名称发生了变化，但这实际上与 "old" 的处理方式非常相似。
最大的变化（从我的角度来看）是替换了以下 "old" 功能
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() 函数似乎读取在上述插补过程中分配给变量的族和 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" for "unordered categorical" 然后按照上面的方法进行。
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 可用的功能），否则我会假设添加一个随机公式中的 effect 应将 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
我的问题是：
- 是否可以使用
mi轻松生成合并的随机效应估计值？和 - 如果是，怎么办？
您可以将 FUN 参数指定给 pool() 函数以更改估算器。在您的情况下，它将是 summary(pool(y ~ x1 + x2 + (1|as.numeric(as.character(group_var))), data = mydf2_imp, FUN = lmer))。这可能有效也可能无效，但它是合法的语法。如果失败，那么您可以使用 complete 函数创建完成的 data.frames，对每个函数调用 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))
只是为了提供一个替代方案，有一个包相当侧重于混合效应模型的 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.