- パッケージとデータの準備
- mice( )で補完後データを作成する
- 補完値を決定するために使用する変数を指定する
- 補完手法を指定する
- 上手く補完できない場合（2021-12-08 追記）
- おわりに
- 参考資料
欠測に対して多重代入（multiple imputation, MI）する場合にはMARの仮定が満たされていることが必要があるが、これはデータから検証することが難しい。
多重代入法を使ったデータ解析の流れは下のとおり。mice( )で補完後データを複数作成し、with( )を使って各々の補完後データに同じ解析を適応し、pool( )を使って結果を統合する。
パッケージとデータの準備
ここではmiceパッケージのnhanes2データ（およびnhanesデータ）を使う（多分CDCのThe National Health and Nutrition Examination Survey（NHANES）のデータ？）
library(mice)
> head(nhanes2) age bmi hyp chl 1 20-39 NA <NA> NA 2 40-59 22.7 no 187 3 20-39 NA no 187 4 60-99 NA <NA> NA 5 20-39 20.4 no 113 6 60-99 NA <NA> 184
age：年齢のカテゴリーbmi：BMIhyp：高血圧の有無chl：血清コレステロール値（mg/dL）
これに対し、nhanesデータではageとhypが数値（numeric）として保持されている。
> head(nhanes, 3) age bmi hyp chl 1 1 NA NA NA 2 2 22.7 1 187 3 1 NA 1 187
後で例として使うので、nhanes2データにht（身長[m], 欠測なし）とwt（体重[kg], 欠測あり）を追加したorig_dデータを作成しておく。
library(tidyverse) orig_d <- nhanes2 %>% mutate(ht = round(runif(n=nrow(nhanes),min=1.5,max=1.8),2), wt = round(bmi*ht*ht,1))
> head(orig_d, 3) age bmi hyp chl ht wt 1 20-39 NA <NA> NA 1.56 NA 2 40-59 22.7 no 187 1.64 61.1 3 20-39 NA no 187 1.80 NA
mice( )で補完後データを作成する
mice( )は連鎖方程式を使って多重代入を行う（multivariate imputation by chained equations, MICE）。MICEでは下の図のように、補完したデータを別の変数の補完値のpredictorとして用いながら順次代入を行なっていく。
(1) 補完前データ
(2) 初期の補完値を適当に決める
(3), (4) 他の変数を使って予測される値を補完する
(5)-(8) 前のステップで代入された値を使いながら順次代入していく
(9) 適当に与えた初期値の影響が薄れるまで、何回か繰り返す（iteration）
mice( )で使いそうな引数は次のとおり。
m：作成する補完後データの個数（冒頭の図だとm=3）。デフォルトは5だが、安定した推定のためには100〜1000くらい必要。method：各変数の補完手法predictorMatrix：ある変数の補完のためにどの変数をpredictorに使うかを行列で示すvisitSequence：iterationの順序をどうするかroman：左→右arabic：右→左monotone：欠測が少ないもの→多いものrevmonotone：結束が多いもの→少ないもの
defaultMethod：数値データ、2値データ、多値データ（順序なし）、多値データ（順序あり）の順に、使用する補完手法を指定するmaxit：iterationの回数。通常15-20あれば十分。seed：乱数シードの指定
methodやpredictorMatrixを1から書くのは大変なので、下のコードのように補間をしないでmice()走らせて（dry runって言うみたい）、作成されたものを書き換えて利用する方が便利。
ini <- mice(orig_d, maxit=0)
> ini$method age bmi hyp chl ht wt "" "pmm" "logreg" "pmm" "" "pmm" > ini$predictorMatrix age bmi hyp chl ht wt age 0 1 1 1 1 1 bmi 1 0 1 1 1 1 hyp 1 1 0 1 1 1 chl 1 1 1 0 1 1 ht 1 1 1 1 0 1 wt 1 1 1 1 1 0 > ini$visitSequence [1] "age" "bmi" "hyp" "chl" "ht" "wt"
補完値を決定するために使用する変数を指定する
補完値の予測に使う変数は多い方が良いが、15-25個もあれば十分とのこと。
選択の基準は、
- 補完後の解析に用いる変数は全て含める
- 欠測の発生に関連している変数は含める
- 説明できる変動が大きい変数は含める
- 欠測が多すぎる変数は除外する
quickpred( )を使うと、欠測発生と他の変数の観測値との相関やproportion of usable cases（PUC）を参考に選んでくれる。
mincor：欠測発生と他の変数の観測値間との相関の下限。デフォルトは0.1minpuc：PUCの下限。デフォルトは0method：相関係数の計算方法。デフォルトは"pearson"になっている。
> quickpred(orig_d, mincor = 0.1, minpuc=0.25) age bmi hyp chl ht wt age 0 0 0 0 0 0 bmi 1 0 0 0 1 0 hyp 1 0 0 0 1 0 chl 1 1 1 0 1 1 ht 0 0 0 0 0 0 wt 1 0 0 0 1 0
このpredictorMatrixは、行が補完される変数（左端の変数名=target）で列が補完値の予測に使う変数（上端の変数名=predictor）を表している。
例えば、「hypの予測にはageとhtのみ用いる」という意味。
また各変数の予測には自分自身を用いないので対角線は0になっている。
補完手法を指定する
methodで指定できる補完手法（built-in univariate imputation methods）は?miceでヘルプを開くと参照できる。
デフォルトで使用される手法は次のとおり。
予測平均マッチング（predictive mean matching）は、線形回帰に基づく予測値が一番近い観測値を補完値として代入する方法（補完方法については野間先生の文献がとても分かりやすかった）。
代入の必要ない（つまり欠測がない）場合は""が指定される。
例えば、hypの補完方法をpmmにしたいときは次のようにmethodを書き換える。
ini <- mice(orig_d, maxit=0) meth <- ini$method meth["hyp"] <- "pmm" imp <- mice(orig_d, method=meth)
他の変数との関係から補完値に制約がある場合（passive imputation）
のように変数間の関係性が決まっている場合は、その関係性に従った補完値を代入する必要がある。
このためにはmethod, predictorMatrix, visitSequenceを適切に指定しなければダメ。
例示のためnhanesデータに変数を追加した。
library(tidyverse) orig_d <- nhanes2 %>% mutate(ht=round(runif(n=nrow(nhanes), min=1.5,max=1.8),2), wt=round(bmi*ht*ht,1))
> head(orig_d,3) age bmi hyp chl ht wt 1 20-39 NA <NA> NA 1.54 NA 2 40-59 22.7 no 187 1.69 64.8 3 20-39 NA no 187 1.51 NA
ht：身長（欠測なし）wt：体重（bmiと同じ箇所に欠測あり）
まずdry runでデフォルトの引数を取り出す。
ini <- mice(orig_d, maxit=0) meth <- ini$method pred <- ini$predictorMatrix seq <- ini$visitSequence
まずmethodを指定する。methodで特定の関係式を指定するためには、~（チルダ）とI()を使う。
meth["bmi"] <- "~I(wt/ht/ht)"
次にpredictorMatrixを指定する。予測の循環を回避するため、BMIの決定に身長を使ったら身長の決定にBMIを使わないようにする。
> pred[c("ht","wt"),"bmi"] <- 0 > pred age bmi hyp chl ht wt age 0 1 1 1 1 1 bmi 1 0 1 1 1 1 hyp 1 1 0 1 1 1 chl 1 1 1 0 1 1 ht 1 0 1 1 0 1 wt 1 0 1 1 1 0
bmiの予測にはhtとwtを使うが、htとwtの予測にはbmiを使わないように指定する。
最後にvisitSequenceを確認する。上記の例の場合、各症例において身長や体重を補完する前にBMIを補完してしまうと、前のiterationの身長・体重を使ってBMIを補完することになり、補完後データで の関係が成立しなくなってしまう。
> seq <- ini$visitSequence > seq [1] "age" "bmi" "hyp" "chl" "ht" "wt"
デフォルトではbmiがht, wtより先に補完されてしまうので、次のように順序を変えておく。
seq <- c("age","hyp","chl","ht","wt","bmi")
上記のように引数を修正して、次のコードで多重補完を行う。
imp <- mice(orig_d, maxit=20, method = meth, predictorMatrix = pred, visitSequence = seq)
作成した補完後データはcomplete( )で確認できる。作成したm個の補完後データのうち、何番目を表示したいかを数字で指定できる。
> complete(imp,1) age bmi hyp chl ht wt 1 20-39 30.31709 no 118 1.54 71.9 2 40-59 22.70000 no 187 1.69 64.8 3 20-39 28.41981 no 187 1.51 64.8 4 60-99 22.28395 yes 229 1.79 71.4 5 20-39 20.40000 no 113 1.74 61.8
上手く補完できない場合（2021-12-08 追記）
これまでの説明に沿ってやったら出来そうな気がしてたんですが、実際にやってみるとmice( )を実行した後も補完されないままNAが残ってしまうことがありました。
色々な人の苦労の跡を拾い読みしていくと、原因の1つとして変数間の共線性の問題があるようです。
補完の変数に完全に一致した（=多重共線性が強い）変数が含まれていると、以下のようなメッセージが出てきてうまく補完できません。
Error in edit.setup(data, setup, ...) : `mice` detected constant and/or collinear variables. No predictors were left after their removal.
多重共線性が問題になるような変数があるかどうかは、
mice:::find.collinear(orig_d)
で確認できます。
これでも上手くいかない場合、カテゴリー変数を文字列としてではなく数値として読み込んだら上手く補完されました。ちょっと面倒ですが、ダメなときは試してみるといいかもしれません。
おわりに
- Passive imputationは条件設定するポイントを抑えれば理解するのは難しくないが、実際に上手く作動させるには工夫が必要だったりする。
- Passive imputationを考えるより、imputation後に変数加工をする方がわかりやすい。
- 抱っこが嫌いなウチの猫も寒くなれば膝に乗ってくれるんではないかと期待中です。
参考資料
van Buuren, S., & Groothuis-Oudshoorn, K. (2011). mice: Multivariate Imputation by Chained Equations in R. Journal of Statistical Software, 45(3), 1–67.
実践的なRの使い方についてまとめられている良著
- Flexible Imputation of Missing Data, 2nd editionのオンラインバージョン： https://stefvanbuuren.name/fimd/
