これまでの内容を、実際のコードと出力で確かめてみます。predictorMatrixの組み方は前のアウトカムを補完してよいかで扱いましたので、今回使用するのはmiceパッケージに入っている25人の教材データで、この規模で医学的な結論を出すことは目的にせず、どの情報を使うと結果がどう変わるのかを、生成数と反復回数の意味を先に整理したうえで、補完方法の指定、統合した回帰係数、反復推移と分布の診断、そして、Rubinの統合則が何を足し合わせているのかまで見ていきます。
生成数と反復回数の決め方
mを増やすと、有限個の補完データセットを使うことによるモンテカルロ誤差を小さくできます。また、maxitを増やすと、各連鎖をより長く更新できます。
では、どちらをどれだけ増やせばよいでしょうか。初期値の影響が残ったままの連鎖を大量に作っても、反復不足への対処にはならず、反対に、1本の連鎖を長く回して1組だけ取り出しても、多重代入としてのデータセット間の不確実性を評価できません。これは、両者が解決する問題が異なるためです。
そこで、後で使う代入間分散を$B$、生成数を$m$とすると、条件付きでほぼ独立な補完を得ている場合、統合した点推定値のモンテカルロ標準誤差は次のように近似できます。
$$\operatorname{MCSE}(\bar Q) \approx \sqrt{\frac{B}{m}} \tag{1}$$
式(1)は、代入間の推定値のばらつきが大きいほど計算結果も乱数に左右されやすく、生成数を増やすほどその揺れが小さくなることを表しています。$\bar Q$は各補完データセットから得た推定値の平均です。
$m$を4倍にしても、この標準誤差はおおむね半分です。そして、これは統計的な標準誤差そのものを半分にするという意味ではありません。対象者を増やす場合とは異なり、生成数を増やしても元のデータの情報量は増えません。
そのため、実際の解析では、点推定値だけでなく、標準誤差や信頼区間の端点も乱数を変えた時にどれくらい動くかを見ます。欠測情報が多い係数や、閾値に近い判断を伴う解析では、より多くの生成数が必要になる場合があります。
また、欠測情報の割合は、データ表の欠測セルの割合とは異なります。欠測が同じ数だけあっても、欠けている変数が主要な説明変数なのか、予測にほとんど寄与しない補助変数なのかで、推定への影響は変わります。
25人のデータによる実行例
miceパッケージのnhanes2データセットを使用して、補完方法の指定と統合までを確かめます。年齢階級、BMI、高血圧の有無、総コレステロールを持つ25人の教材用データで、年齢階級と高血圧の有無は因子型です。
ここでは、アウトカムにした総コレステロールにも欠測があります。最初のコードは、その総コレステロールも補って25人全員を解析する通常の多重代入で、アウトカムを予測因子としてのみ使う解析やMIDではありません。
以下のコードはmiceパッケージ3.18.0で実行しています。
set.seed(42)
library(mice)
study <- nhanes2
colSums(is.na(study))
sum(complete.cases(study))
imputed <- mice(study, m = 30, maxit = 20, seed = 42, printFlag = FALSE)
imputed$method
imputed$loggedEvents
models <- with(imputed, lm(chl ~ age + bmi + hyp))
combined <- pool(models)
summary(combined, conf.int = TRUE)
coef(summary(lm(chl ~ age + bmi + hyp, data = study)))#> age bmi hyp chl
#> 0 9 8 10
#>
#> 完全ケース数
#> 13
#>
#> 補完方法
#> age bmi hyp chl
#> "" "pmm" "logreg" "pmm"
#>
#> loggedEvents
#> NULLBMIの欠測は9人で36.0%、高血圧の有無は8人で32.0%、総コレステロールは10人で40.0%です。そして、アウトカムが観測されている人は15人、全変数がそろった完全ケースは13人なので、完全ケース解析では25人のうち12人が外れます。
ここで、年齢階級の補完方法が空文字なのは、年齢階級に欠測がないためです。年齢階級を予測に使っていないという意味ではありません。また、BMIと総コレステロールにはPMM、高血圧の有無にはロジスティック回帰が選ばれています。
loggedEventsは今回NULLでした。ただし、記録対象のイベントがなかったことから、MARが成立した、モデルに問題がない、と判断することはできません。
統合した回帰係数
年齢階級と高血圧の有無を調整し、BMIと総コレステロールの関連を線形回帰で推定しました。出力のうち、切片を除く係数を小数第3位へ丸めると、次の値になります。
#> MICEによる統合結果
#> 項目 推定値 標準誤差 95%信頼区間
#> age40-59 48.085 21.180 2.242 ～ 93.928
#> age60-99 70.910 28.309 6.039 ～135.781
#> bmi 6.071 2.189 1.338 ～ 10.803
#> hypyes -4.337 25.289 -62.676 ～ 54.002
#>
#> 完全ケース解析のBMI係数
#> 推定値 7.160
#> 標準誤差 2.201このうちBMIの統合係数6.071は、この線形モデルで年齢階級と高血圧の有無を同じとした場合、BMIが1高いことに対応する総コレステロールの平均差です。
完全ケース解析のBMI係数は7.160で、多重代入では6.071となりました。ただし、この差には欠測共変量の補完だけでなく、アウトカム未観測者10人を補完アウトカムとともに解析へ含めたことも関わっています。共変量だけを補った効果として説明することはできず、値が小さくなることもMICEの一般的な性質ではありません。
また、総コレステロールが未観測だった10人の真値が判明したわけではなく、結果は補完モデルと欠測の仮定に依存しています。完全ケース解析と多重代入のどちらが真値に近いかは、この比較だけでは判断できません。
そして、今回は多重代入後の標準誤差がわずかに小さくなっています。不完全な対象者の情報を使う効果と、欠測に由来する不確実性を加える効果が同時に働くため、完全ケース解析との標準誤差の大小には一定の方向がありません。
同じ補完からMIDを計算する
通常の多重代入とMIDで何が変わるかを追うため、同じ30組の補完データセットに対して、元から総コレステロールが観測されていた15人だけで回帰モデルを推定します。以下では元データの観測・欠測の区別を保存し、補完後の値から対象者を選ばないようにしています。
observed_outcome <- !is.na(study$chl)
mid_models <- with(imputed, lm(chl ~ age + bmi + hyp,
subset = observed_outcome))
summary(pool(mid_models), conf.int = TRUE)BMI係数に絞ると、3通りの結果は表2のようになります。各補完データセット内の行順は元データと同じで、全30組のMIDで解析人数が15人になることも確認しています。
| 解析方法 | 回帰分析の人数 | 解析に含むアウトカム | BMI係数 | 標準誤差 | 95%信頼区間 |
|---|---|---|---|---|---|
| 完全ケース解析 | 13人 | 観測値13人分 | 7.160 | 2.201 | 2.085〜12.235 |
| 通常の多重代入 | 各組25人 | 観測値15人分と補完値10人分 | 6.071 | 2.189 | 1.338〜10.803 |
| MID | 各組15人 | 観測値15人分 | 6.880 | 2.312 | 1.512〜12.247 |
ここで、MIDでは、通常の多重代入に使った補完値をそのまま利用し、回帰分析へ渡す対象者だけを変えています。一方、完全ケース解析では共変量にも欠測のない13人に絞るため、アウトカム観測者に残る共変量の欠測も補うMIDとは人数が異なります。
ただし、MIDの係数が両者の間にあることから、妥当性の高い妥協案と読むことはできません。一つの小標本で比較した結果なので、MIDの選別に必要な仮定や、補助変数に伴う偏りが解決したとは判断できません。
なお、補完後に欠測判定をしても、総コレステロールの欠測は既に埋まっているため、元の未観測者を識別できません。実際の研究では、元の欠測フラグと対象者IDを保存し、並べ替えや結合でフラグと対象者の対応が変わらないようにします。
反復推移による診断
反復ごとに補完値の平均を描くと、初期値の影響が残っていないか、連鎖間で異なる領域にとどまっていないかを見られます。以下では、30本のうち最初の5本について、BMIの欠測9人へ入れた値の平均を描きます。
trace_mean <- imputed$chainMean["bmi", , 1:5]
png("mice-trace-200dpi.png", width = 1100, height = 688, res = 200, type = "cairo")
plot(1:20, y = trace_mean[, 1], type = "l", ylim = range(trace_mean),
xlab = "Iteration", ylab = "Mean of imputed BMI", col = 1)
for (chain in 2:5) lines(1:20, y = trace_mean[, chain], col = chain)
legend("top", inset = c(0, -0.28), legend = paste("Chain", 1:5),
col = 1:5, lty = 1, ncol = 5, cex = 0.6, bty = "n", xpd = NA)
dev.off()
ここで、多重代入では反復ごとに値が動くため、線が一直線へ落ち着くことを期待しません。むしろ、ある連鎖だけが長く高い領域に残る、反復が進むほど全体が上昇する、といったパターンがないかを見ます。
ただし、平均が安定していても、分散や裾の値、変数間の関連が安定しているとは限らず、図1は見やすさのため5本だけを描いているので、実際の診断では残りの連鎖と、BMI以外の補完対象も確認します。
miceには自己相関や潜在尺度縮小因子を計算する診断関数もありますが、これらは連鎖の挙動に関する情報です。欠測機構やモデル指定の正しさを判定する検定ではないことにご留意ください。
そこで、反復不足が疑われるなら反復を延ばし、その後も連鎖が分離するなら、観測数の少ないカテゴリ、共線性、極端な予測値、補完モデルの指定を調べます。maxitを増やす操作だけで済ませないことが、診断の次の段階です。
観測値と補完値の分布
補完BMIの範囲や集中する場所を確認するため、観測された16人のBMIと、9人分を30組生成した補完BMIを比較します。
observed <- study$bmi[!is.na(study$bmi)]
filled <- unlist(imputed$imp$bmi, use.names = FALSE)
png("mice-distribution-200dpi.png", width = 1100, height = 688, res = 200, type = "cairo")
plot(density(observed, from = 15, to = 40), xlim = c(15, 40), ylim = c(0, 0.15),
xlab = "BMI", ylab = "Density", main = "", lwd = 2)
lines(density(filled, from = 15, to = 40), col = 2, lwd = 2)
legend("topright", legend = c("Observed", "Imputed"), col = 1:2, lwd = 2)
dev.off()
図2では、補完値が集中する位置と観測値の分布が完全には一致していません。この違いを見ただけで補完失敗と判断するのは注意が必要です。欠測者の年齢や高血圧の構成が観測者と異なれば、その情報を反映した補完BMIの周辺分布が異なることはあり得ます。
そのうえ、PMMは観測値をドナーとして使うため、観測値と補完値の範囲が似ること自体は方法の仕組みから期待されます。範囲が似ていることを、欠測値が正しく復元された証拠にはできません。
さらに、密度曲線は平滑化の設定にも左右されます。特に観測者が16人と少ない今回の例では、曲線の小さな山や谷へ医学的な意味を付けることは避けます。曲線の裾が実測値の範囲を少し越えても、それは平滑化によって生じているので、PMMがその値を生成したことを意味しません。
そのため、実際の診断では、年齢階級別や群別の分布、補完に使った予測変数との散布図、まれなカテゴリの頻度も見ます。全体の分布が似ていても、群間差や変数間の関係が失われている場合があるためです。
Rubinの統合則
30組の補完データセットから30個のBMI係数を得た後、その平均だけを報告すると、欠測値の候補によって推定がどれくらい変わるかが抜け落ちます。
$j$番目のデータセットから得た推定値を$\hat Q_j$、その推定分散を$U_j$、生成数を$m$とすると、統合した推定値と平均的な代入内分散は次の形です。
$$\bar Q=\frac{1}{m}\sum_{j=1}^{m}\hat Q_j, \qquad \bar U=\frac{1}{m}\sum_{j=1}^{m}U_j \tag{2}$$
式(2)は、各データセットから得た推定値を平均し、同時に、それぞれの解析で計算された分散も平均することを表しています。$\bar U$は、補完データセットごとの解析に残る標本変動の大きさです。
次に、推定値が補完データセット間でどれくらい異なるかを計算します。
$$B=\frac{1}{m-1} \sum_{j=1}^{m}(\hat Q_j-\bar Q)^2 \tag{3}$$
式(3)は、補完の仕方によって推定値が動く程度を、代入間分散$B$として求めています。同じ観測値を使っていても、欠測部分の候補が変わることで生じる推定のばらつきです。
そして、統合分散と標準誤差は、代入内分散と代入間分散から求めます。
$$T=\bar U+\left(1+\frac{1}{m}\right)B, \qquad \operatorname{SE}(\bar Q)=\sqrt{T} \tag{4}$$
式(4)は、各解析の分散に、欠測部分の候補による推定のばらつきを加えることを表しています。$1+1/m$には、有限個の補完データセットを使うことによる補正が含まれます。
miceのpool()はこの統合を行い、信頼区間には小標本を意識した自由度の調整を使います。そのため、今回のような小さいデータでは、常に「推定値±1.96×標準誤差」と同じ区間になるわけではありません。
実行値との照合
今回の通常の多重代入によるBMI係数では、統合に使った量は次の値でした。
#> 統合推定値 6.070866
#> 平均代入内分散 3.491491
#> 代入間分散 1.256171
#> 統合分散 4.789534
#> 統合標準誤差 2.188500
#> 点推定値の近似MCSE 0.204627ここで、平均代入内分散3.491491だけで標準誤差を求めると、代入間分散1.256171に由来する部分が落ちます。今回はRubinの式から別に計算した統合推定値と統合分散が、pool()の返り値と一致することを確認しました。
さらに、近似MCSEは約0.205で、統合標準誤差約2.189の約9.3%に相当します。30組で得た点推定値にも、この程度の計算上の揺れが残り得るため、係数の小数第3位までを安定した研究上の差として読むことはできません。
欠測率と欠測情報割合の違い
変数別の欠測率は、測定対象のうち値が欠けている人の割合です。欠測情報割合、FMIはそれとは異なり、特定の係数などを推定するための情報が、欠測によってどれくらい失われているかを、補完モデルの下で評価します。同じ補完データでも、BMI係数と年齢階級の係数ではFMIが異なり得ます。
この違いを数える量として、大標本では$\lambda=(1+1/m)B/T$があります。これは、統合分散のうち、代入間のばらつきに関わる部分の割合で、miceの出力では小標本調整を含むFMIと区別されています。
実際、今回のBMI係数では、この$\lambda$は約0.271で、miceが小標本調整を加えて返すFMIは約0.363でした。BMIの欠測率36.0%、アウトカムの欠測率40.0%、完全ケース解析で外れる割合48.0%とは定義が異なり、数値が近い場合も同じ指標とは扱いません。また、$\lambda$と調整済みFMIも区別して記載します。
なお、補助変数を追加して欠測値を予測する情報が増えれば、欠測率を変えずに欠測情報割合が小さくなる場合があります。しかし、不十分なモデルが欠測部分の不確実性を小さく見積もることもあるため、FMIが低いという一点で補完の正しさを評価せず、モデルの根拠と診断を合わせて読みます。
データを平均しない理由
30組のデータセットについて、各人の補完値を平均して1組にまとめると、その人の欠測値にどの程度の幅があったかが失われます。また、30組を縦に結合して通常の回帰分析へ渡すと、同じ対象者を繰り返し独立な観測として扱ってしまいます。
多重代入では、各組で解析してから推定値と分散を統合します。
このとき、ロジスティック回帰なら対数オッズ比、Cox回帰なら対数ハザード比の尺度で係数と分散を統合し、その後で指数変換します。複数カテゴリを持つ変数全体の検定など、係数をまとめて評価したい場合には、係数間の共分散も使う多重代入用の検定が必要です。
MARという仮定そのものを動かす感度分析と、論文への書き方は、次のMNARの感度分析と論文への記載で扱います。
参考文献
van Buuren S, Groothuis-Oudshoorn K. mice: Multivariate Imputation by Chained Equations in R. Journal of Statistical Software. 2011;45(3):1-67.
White IR, Royston P, Wood AM. Multiple imputation using chained equations: Issues and guidance for practice. Statistics in Medicine. 2011;30(4):377-399.
Sterne JAC, White IR, Carlin JB, et al. Multiple imputation for missing data in epidemiological and clinical research: potential and pitfalls. BMJ. 2009;338:b2393.
