- パッケージとデータの準備
- 欠測の発生と他の変数との相関を調べる
- md.pattern( ), md.pairs( )で欠測のパターンを示す
- aggr( )を使って欠測に関する統合的なグラフを描く
- marginplot( )で2変数間の観測値・欠測をプロットする
- matrixplot( )で全症例の欠測パターンを図示する
- おわりに
- 参考資料
欠測は無いに越したことはないが、実際のデータでは避けられない。欠測の割合が個々の変数においては小さかったとしても、変数が多くなると欠測が全くない症例はそれなりの数になる。ほとんどの解析手法は欠測がないことを前提にしているので、(1)欠測がある症例を除外するか、(2)欠測に適切な補完値を代入する（imputation）か、のどちらかが必要になる。
欠測は発生するメカニズムによって、次のように分類される。
- Missing complete at random（MCAR）：完全にランダムに欠測が発生している
- Missing at random（MAR）：ある変数の欠測の発生は他の変数の観測値と関連するが、欠測している変数の真の値とは関連しない（つまり「小さい値ほど観察されにくい」のようなことはない）
- Missing not at random（MNAR）：欠測している変数の真の値と欠測の発生が関連する
MARなら補完可能だが、データからMARといえるか検証することは困難。MNARだとどう頑張っても適切に対応することが難しい。
パッケージとデータの準備
ここではVIMパッケージ内のsleepデータを使う。sleepデータは哺乳類の睡眠時間をまとめたデータ。
library(mice) library(VIM)
> head(sleep) BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger 1 6654.000 5712.0 NA NA 3.3 38.6 645 3 5 3 2 1.000 6.6 6.3 2.0 8.3 4.5 42 3 1 3 3 3.385 44.5 NA NA 12.5 14.0 60 1 1 1 4 0.920 5.7 NA NA 16.5 NA 25 5 2 3 5 2547.000 4603.0 2.1 1.8 3.9 69.0 624 3 5 4 6 10.550 179.5 9.1 0.7 9.8 27.0 180 4 4 4
BodyWgt：体重（kg）BrainWgt：脳重量（g）NonD：ノンレム睡眠の長さDream：レム睡眠の長さSleep：睡眠の長さ（=NonD+Dream）Span：一生の長さ（年）Gest：妊娠期間（日）Pred,Exp,Danger：それぞれ「捕食されやすさ」、「睡眠中の曝露」、「総合的な危険度」を5段階で表したものらしい。R in Actionには次のような説明あり。
The ecological variables included degree to which species were preyed upon (Pred), degree of their exposure while sleeping (Exp), and overall danger (Danger) faced. The ecological variables were measured on 5-point rating scales that ranged from 1 (low) to 5 (high).
欠測の発生と他の変数との相関を調べる
基本関数のcor()を使う。見やすくなるようにround()を使って3桁に丸めてしまう。
> round(cor(x=!is.na(sleep), y=sleep, use="pairwise.complete.obs"),3) BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger BodyWgt NA NA NA NA NA NA NA NA NA NA BrainWgt NA NA NA NA NA NA NA NA NA NA NonD -0.227 -0.179 NA 0.189 0.080 -0.083 -0.202 -0.048 -0.245 -0.065 Dream -0.223 -0.163 NA NA 0.080 -0.060 -0.051 0.068 -0.127 0.067 Sleep -0.002 -0.008 NA 0.189 NA -0.005 -0.160 -0.202 -0.261 -0.209 Span 0.058 0.079 0.043 -0.117 -0.096 NA 0.175 -0.023 0.193 0.067 Gest 0.054 0.073 0.046 -0.228 -0.040 0.065 NA 0.201 0.193 0.204 Pred NA NA NA NA NA NA NA NA NA NA Exp NA NA NA NA NA NA NA NA NA NA Danger NA NA NA NA NA NA NA NA NA NA
行方向に並んだ変数が欠測することが、列方向の変数の値とどれくらい相関しているかを示している。
例えば、NonDはBodyWgt（r=-0.227）やExp（r=-0.245）の観測値との相関が他と比べると強い。また、BodyWgtは欠測がないので相関係数が計算されない。
md.pattern( ), md.pairs( )で欠測のパターンを示す
md.pattern( )を使う
miceパッケージのmd.pattern( )を使って欠測のパターンを見る。
> md.pattern(sleep, rotate.names=TRUE) BodyWgt BrainWgt Pred Exp Danger Sleep Span Gest Dream NonD 42 1 1 1 1 1 1 1 1 1 1 0 9 1 1 1 1 1 1 1 1 0 0 2 3 1 1 1 1 1 1 1 0 1 1 1 2 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 0 1 0 0 3 1 1 1 1 1 1 1 0 0 1 1 2 2 1 1 1 1 1 0 1 1 1 0 2 2 1 1 1 1 1 0 1 1 0 0 3 0 0 0 0 0 4 4 4 12 14 38
各行は欠測のパターンを示している。左端の数値はそのパターンに該当する症例数、右端の数値は欠測している変数の個数。
例えば2行目はDreamとNonDの2つの変数が欠測しているものが9症例あるということ。
引数のrotate.names=TRUEは次の図の変数名を垂直にして横と重ならないようにするためのもの。
VIMパッケージのaggr( )で同じようなグラフが描ける（後述）。
md.pairs( )を使う
次にmiceパッケージのmd.pairs関数を使って、2つの変数の組み合わせごとに、両方とも観察（rr）、行の変数は観察・列の変数は欠測（rm）、行は欠測・列は観察（mr）、両方とも欠測（mm）の例数を示す。
> md.pairs(sleep) $rr BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger BodyWgt 62 62 48 50 58 58 58 62 62 62 BrainWgt 62 62 48 50 58 58 58 62 62 62 NonD 48 48 48 48 48 45 44 48 48 48 ~~~（省略）~~~ $rm BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger BodyWgt 0 0 14 12 4 4 4 0 0 0 BrainWgt 0 0 14 12 4 4 4 0 0 0 NonD 0 0 0 0 0 3 4 0 0 0 ~~~（省略）~~~ $mr BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger BodyWgt 0 0 0 0 0 0 0 0 0 0 BrainWgt 0 0 0 0 0 0 0 0 0 0 NonD 14 14 0 2 10 13 14 14 14 14 ~~~（省略）~~~ $mm BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger BodyWgt 0 0 0 0 0 0 0 0 0 0 BrainWgt 0 0 0 0 0 0 0 0 0 0 NonD 0 0 14 12 4 1 0 0 0 0 ~~~（省略）~~~
次のようにして「ある変数について欠測がある場合に、他の変数でどれくらいの割合が観察されているか」（proportion of usable cases, PUC）を計算できる（ここでも見やすくするためにround( )を使って丸めている）。
> p <- md.pairs(sleep) > round(p$mr/(p$mr + p$mm),3) BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger BodyWgt NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN BrainWgt NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NonD 1 1 0.00 0.143 0.714 0.929 1.00 1 1 1 Dream 1 1 0.00 0.000 0.833 0.917 1.00 1 1 1 Sleep 1 1 0.00 0.500 0.000 1.000 1.00 1 1 1 Span 1 1 0.75 0.750 1.000 0.000 0.75 1 1 1 Gest 1 1 1.00 1.000 1.000 0.750 0.00 1 1 1 Pred NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN Exp NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN Danger NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
例えばNonDが欠測している場合、Sleepは約71%の症例で観察されているのに対し、Dreamは約14%の症例でしか観察されていない（残りの約86%は同時に欠測してしまっている）。
aggr( )を使って欠測に関する統合的なグラフを描く
VIMパッケージのaggr( )を使う。 左側には各変数の欠測割合（または欠測数）、右側にはmd.pattern( )で出力されるような欠測パターン別発生割合（または発生数）が示される。
次の引数を指定できる。
prop：割合で表示する場合は=TRUEとする（デフォルト）number：割合や数を表示するかどうか。デフォルトは=FALSE。
aggr(sleep)
次のようにしてmd.pattern( )のような結果を出力することもできる（出力省略）が、md.pattern( ) の結果の方が見やすいと思う。
a <- aggr() summary(a)
marginplot( )で2変数間の観測値・欠測をプロットする
VIMデータのmarginplot( )を使って、2つの変数の欠測情報を含めたプロットを作成する。下の例ではGestとDreamを指定した。
色やマーカー種類を変えたいときは、col=c("blue","red"), pch=c(10,20)のように観測・欠測の順に指定する。
そのままだと散布図のプロットが中抜きで見にくいので、pch=20などの方が見やすくなる。
marginplot(sleep[c("Gest", "Dream")], pch=20)
- 右上のプロット（上図では青）：両方とも観察されたデータの散布図
- 周辺のプロット（上図では赤）：もう一方の変数が欠測のものの値をプロットしたもの
- 周辺の箱ヒゲ図：もう一方の変数が観察されているもののデータは青で、欠測のものは赤で描かれている
marginmatrix(sleep)とすれば全ての変数の組み合わせについて上記のようなプロットを作成してくれるが、変数が多いと見にくい。
matrixplot( )で全症例の欠測パターンを図示する
必要があれば以下の引数を指定することができる。
sortby：ソートに使いたい変数を指定col：欠測を表す部分の色を指定（デフォルトは多分"red"）
matrixplot(sleep, sortby="Danger", col="blue")
上記のプロットでは、欠測は青、観測値は0〜1に再スケール化してグレースケールで示されている。
おわりに
- 欠測パターンのプロットにはVIMパッケージが便利。
参考資料
van Buuren, S., & Groothuis-Oudshoorn, K. (2011). mice: Multivariate Imputation by Chained Equations in R. Journal of Statistical Software, 45(3), 1–67.
実践的なRの使い方についてまとめられている良著
- Rの解析に役に立つ記事がたくさんあります。VIMパッケージの他の関数の使い方も説明されていました。 www.karada-good.net
