欠測に対して適切な値が代入されたかを評価する。まずはnhanes2データを使って補完データを準備。
library(mice) imp <- mice(nhanes2, m=5, maxit=20)
補完値の収束状況を確認する
mice( )が行う連鎖方程式を使って多重代入法(multivariate imputation by chained equations, MICE)では、補完したデータを別の変数の補完値のpredictorとして用いながら順次代入を行なっていく。
初期値の影響が取れた値が補完されているかを次のグラフで確認する。
plot(imp)
maxit=20
と指定したので、各欠測に対して20回(20周)補完されている- 色分けされた各線(鎖)は作成された各補完データセット(ここでは
m=5
)を表している - 連続変数のみ平均と標準偏差が示される。平均・標準偏差以外を表示したいときは、補完後データを抽出して自分でグラフ化することもできるみたい。
以下、収束プロットの例とその評価。
- 左上:初期値のまま固定されているのでダメ
- 右上:鎖が十分に混ざっていないのでダメ
- 左下:iterationが15回以降くらいになると収束しているので
maxit=20
なら十分 - 右下:鎖が十分に混ざっていてOK
complete( )を使って補完後のデータセットを抽出する
miceパッケージのcomplete( )で実際に作成された補完後データセットを抽出できる。
指定する主な引数は次のとおり。
action
- m以下の数字nを指定するとn番目の補完後データセットを抽出
"all"
:全てのデータセットのリスト"long"
:全てのデータを縦につなげたもの。何番目の補完後データかを示す.imp
という変数(オリジナルを含めるときは0)と、症例IDを示す.id
という変数が追加される。"stack"
:"long"と同じだが.imp
,.id
が追加されない"broad"
:全てのデータを横につなげたもの。補完後データをひとかたまりにつなげている。"repeated"
:上記を同じ変数が隣り合う様に並べ替えたもの
include
:補完前のオリジナルを含めたいときは=TRUE
を指定。
次のようにすれば意図しない補完がされていないかザッと調べられる(と思う)。
summary(complete(imp, action="long"))
補完データと元データの分布を比較する
densityplot( )を使って、代入された値の分布と元々の観察された値の分布を重ねてプロットすることで、補完のアルゴリズムが妥当かどうかチェックすることができる(densityplot( )は元々latticeパッケージの関数を利用)。
densityplot(imp)
欠測発生の傾向スコアを使う
もし多重代入に用いたモデルが適切であれば、欠測発生に対する傾向スコアが等しい観察値と補完値は似た様な値になるはず。
bmi_miss <- is.na(nhanes2$bmi) fit_bmi <- with(imp, glm(bmi_miss~age+hyp+chl, family=binomial(link="logit"))) ps <- rep(rowMeans(sapply(fit_bmi$analyses, fitted.values)) ,6) xyplot(imp, bmi~ps | .imp, pch=c(1,20))
コードの概要は以下のとおり。
- BMIが観察されていれば0, 欠測なら1となる2値変数
bmi_miss
を作成する - 作成した補完データのセット(オリジナルも含む)において、
bmi_miss
以外の変数を説明変数としたロジスティック回帰モデルを当てはめる - 2で当てはめたモデルから予測値
fitted.values
を取り出し、rowMeans()
を使って各症例毎の平均をとる。この後で、合計m+1個のデータのセット(オリジナルと補完データ)ごとにプロットを作るので、ここではrep()
で6回複製している。 - データセットごとにプロットを作成する。観察値は
pch=1
で、補完値はpch=20
で描いた。
青が観察値、赤が補完値。補完値が相対的に右側にあるのは当たり前。
傾向スコアの値が同じときに観察値と補完値の分布は近くなっていれば、上手く補完できているということ。観察値と補完値の分布が重なっていなければ、欠測値の予測に重要な変数が含まれていない可能性がある。
おわりに
- せめて分布がだいたい同じになっているかは確認した方がよい
参考資料
- van Buuren, S., & Groothuis-Oudshoorn, K. (2011). mice: Multivariate Imputation by Chained Equations in R. Journal of Statistical Software, 45(3), 1–67.