miパッケージを使用して帰属ランダム効果モデル推定値をプールすることはできますか？

Question

miパッケージの過去2〜3年のうちのいくつかの時点でかなり大きな書き換えが行われているようです。

物事の「古い」方法は、次のチュートリアルではよく概説である：（リーパーのシミュレーションデモにこだわっ）物事のhttp://thomasleeper.com/Rcourse/Tutorials/mi.html

「新しい」方法は、次のようになります。

​​

関数名は変更されていますが、これは実際には "古い"方法と似ています。

（私の見晴らしから）最大の変化は、以下の "古い" の機能の代替

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()（またはLeeerの例に従えばcoef.mi()）を使用してプールすることができました。

現在のバージョンmi（私はv1.0がインストールされています）では、これらの最後のステップが1つの関数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 

これは、模擬ベータ値（2と20）の回復に近づいているようです。つまり、期待通りに動作しています。

グループ化変数を取得するために、ちょっとしたシミュレーションをランダムに実行してみましょう。 missing_data.frame()ので

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 

私が「順不同カテゴリ」の"un"に再割り当てした後、上記のように進行するmiからchange()関数を使用し、連続としてgroup_varをintepretingているように見えます。バージョンmiの1.0以前のバージョン（lmer.miとglmer.miで利用可能、すなわち機能）の機能を削除していない限り

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

# impute 
mydf2_imp <- mi(mydf2_mdf) 

ここで、Iは、式中のランダム効果の添加が適切に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 

私の質問はありますか？A：はい場合

1. は、それが簡単にmiを使用してプールされたランダム効果推定値を生成することができる、と
2. 、どのように？

Answer 1

推定値を変更するpool()関数の引数にFUN引数を指定できます。あなたの場合は、summary(pool(y ~ x1 + x2 + (1|as.numeric(as.character(group_var))), data = mydf2_imp, FUN = lmer))になります。これは実際には機能しないかもしれませんが、法的な構文です。それが失敗した場合、あなたは、作成完了data.framesに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))

Answer 2

ようなものになるだろう結果を自分で、平均化することができ、そのパッケージがありますミックス・エフェクト・モデルのために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.