ねこすたっと

ねこの気持ちと統計について悩む筆者の備忘録的ページ。

多重代入法(4)：補完後データを解析する（miceパッケージ）[R]

missing data handling tidyverse R neco

多重代入法を使ったデータ解析の流れのおさらい。

図1：多重補完を使ったデータ解析の流れ

基本的には、(1)mice( )で作成した補完後データに対して、(2)with( )を使って各々の補完後データに同じ解析を適応し、(3)pool( )を使って結果を統合するが、後述のtidyverseパッケージを使う方法もある。

基本：mice( )→with( )→pool( )

まずmice( )で欠測を補完する。ここではmiceパッケージのnhanes2データを使う。

library(mice)
imp <- mice(nhanes2, m=5, maxit=20, seed=1234)

次にwith( )で各々の補完後データに同じ解析を適応する。

fit <- with(imp, lm(chl ~ age + bmi + hyp))

glm()なら下のを参考に。

fit <- with(imp, glm(chl ~ age + bmi + hyp, family = gaussian(link="identity")))

最後にpool( )で結果を統合する。

est <- pool(fit)

結果は下のように表示される。

> est
Class: mipo    m = 5 
         term m   estimate       ubar           b           t dfcom        df       riv    lambda       fmi
1 (Intercept) 5 -25.989753 2970.56604 716.6785355 3830.580281    20 12.016685 0.2895119 0.2245128 0.3277962
2    age40-59 5  64.485671  324.47266 230.7211305  601.338020    20  6.473135 0.8532779 0.4604155 0.5743344
3    age60-99 5  94.013261  475.61448 433.3687166  995.656944    20  5.469221 1.0934117 0.5223109 0.6351168
4         bmi 5   6.855359    3.56238   0.7391608    4.449373    20 12.766118 0.2489889 0.1993524 0.3009180
5      hypyes 5 -24.262019  358.80462 199.4101443  598.096796    20  7.616070 0.6669150 0.4000894 0.5131087

結果の項目については?mipoで得られる。最後の3つの指標については資料を参照のこと。

  • estimate：統合した推定量（pooled estimate）
  • ubar：代入データ内分散（within-imputation variance of estimate）
  • b：代入データ間分散（between-imputation variance of estimate）
  • t：総分散（total variance of estimate）
  • dfcom：complete dataの自由度
  • df：t統計量の自由度
  • riv：relative increase in variance. RIV = \frac{V_B + \frac{V_B}{m}}{V_W}
  • lambda：proportion attributable to the missingness. Lambda = \frac{V_B + \frac{V_B}{m}}{V_T}
  • fmi：fraction of missing information. FMI=\frac{RIV+\frac{2}{df+3}}{1+RIV}
> summary(est)
         term   estimate std.error  statistic        df     p.value
1 (Intercept) -25.989753 61.891682 -0.4199232 12.016685 0.681952086
2    age40-59  64.485671 24.522194  2.6296860  6.473135 0.036419856
3    age60-99  94.013261 31.554032  2.9794373  5.469221 0.027593435
4         bmi   6.855359  2.109354  3.2499811 12.766118 0.006464121
5      hypyes -24.262019 24.456018 -0.9920674  7.616070 0.351632239

別の方法1：mice( )→complete( )→by( )

多重代入したデータを縦に繋げて、補完のセット毎にモデルを当てはめて統合する。パイプ演算子（%>%）が使えるのでコードがスッキリ見やすい。
by( )の中にある.$.impの、$前の.は%>%で受けるオブジェクト（前のステップで作成されたもの）を意味している。

library(tidyverse)
est <- nhanes2 %>%
  mice(seed = 1234, m=5, maxit=20, print = FALSE) %>%
  mice::complete("long", include = FALSE)  %>%
  by(as.factor(.$.imp), glm, formula = chl ~ age + bmi + hyp, family = gaussian("identity")) %>%
  pool()
summary(est)

ちなみに、complete("long")で取り出したデータは次のようになっている。include=FALSEとしてオリジナルのデータセットが含まれないように。

    .imp .id   age  bmi hyp chl
1      1   1 20-39 27.5  no 187
2      1   2 40-59 22.7  no 187
3      1   3 20-39 30.1  no 187
~~~~~（省略）~~~~~
123    5  23 20-39 27.5  no 131
124    5  24 60-99 24.9  no 284
125    5  25 40-59 27.4  no 186

一旦この状態で取り出して、変数を加工・追加してからby( )に渡してやるようにすれば、passive imputationを気にせずに解析ができそう。

別の方法2：mice( )→group_by( )→do( )

別法1と同様に%>%を使った書き方。ここでも前のステップで作成されたオブジェクトを受ける.が登場する。（do( )の中のdata=.とas.list( )の次の.[[2]]。.impの.は変数名の一部。）

est <- nhanes2 %>%
  mice(seed = 1234, m=5, maxit=20, print = FALSE) %>%
  mice::complete("long") %>%
  group_by(.imp) %>%
  do(model = lm(formula = chl ~ age + bmi + hyp, data = .)) %>%
  as.list() %>%
  .[[2]] %>%
  pool()
summary(est)

.[[2]]は前のステップで作成されたリストの2番目の要素を使う、という意味。1番目はimputation番号（ここでは1〜5）で、モデル当てはめの結果は2番目に格納されているため。

おわりに

  • 個人的には別法1が使いやすそうだと思った。
  • 我が家の猫は最近よく鳴いて呼ぶようになりました。

ランキング参加中
統計

参考資料