8

次のような線形混合モデルで、各実験単位 (plot) のそれぞれの標準誤差で係数 (b0 と b1) を抽出するにはどうすればよいでしょうか。

線形モデルへのより適切な適合

これと同じデータセット (df)、および適合モデル (fitL1) の場合: どうすればこのようなデータ フレームを取得できますか...

   plot    b0      b0_se   b1    b1_se 
    1    2898.69   53.85   -7.5  4.3

   ...    ...       ...     ...   ...
4

2 に答える 2

11

最初のコメントは、これは実際には自明ではない理論的な問題であるということです。r-sig-mixed-modelsには、技術的な詳細の一部に入るかなり長いスレッドがあります。少し怖くても、絶対に見ておくべきです。基本的な問題は、各グループの推定係数値が、そのグループの固定効果パラメーターと BLUP/条件付きモードの合計であり、オブジェクトの異なるクラスであるということです (1 つはパラメーター、もう 1 つはオブジェクトの条件付き平均です)。これにより、技術的な問題が発生します。

2番目のポイントは、(残念ながら) でこれを行う簡単な方法がわからないlmeため、私の答えはlmer(lme4パッケージから) を使用することです。

最も簡単なことを行い、固定効果パラメーターと BLUP の間の (おそらく不明確な) 共分散を無視することに慣れている場合は、以下のコードを使用できます。

2 つの選択肢は、(1) モデルをベイジアン階層アプローチ (MCMCglmmパッケージなど) に適合させ、各レベルの事後予測の標準偏差を計算する (2) パラメトリック ブートストラップを使用して BLUP/条件付きモードを計算し、ブートストラップ分布の標準偏差。

いつものように、このアドバイスには保証がないことを覚えておいてください。

library(lme4)
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
cc <- coef(fm1)$Subject
## variances of fixed effects
fixed.vars <- diag(vcov(fm1))
## extract variances of conditional modes
r1 <- ranef(fm1,condVar=TRUE)
cmode.vars <- t(apply(cv <- attr(r1[[1]],"postVar"),3,diag))
seVals <- sqrt(sweep(cmode.vars,2,fixed.vars,"+"))
res <- cbind(cc,seVals)
res2 <- setNames(res[,c(1,3,2,4)],
                 c("int","int_se","slope","slope_se"))
##          int   int_se     slope slope_se
## 308 253.6637 13.86649 19.666258   2.7752
## 309 211.0065 13.86649  1.847583   2.7752
## 310 212.4449 13.86649  5.018406   2.7752
## 330 275.0956 13.86649  5.652955   2.7752
## 331 273.6653 13.86649  7.397391   2.7752
## 332 260.4446 13.86649 10.195115   2.7752
于 2014-10-05T19:43:57.873 に答える