多変量Cox回帰を組もうとして変数を並べると、最後に数人から数十人が消えることがあります。追跡時間とイベントは全員分そろっているのに検査値が欠け、解析対象者数が少し減るだけなら気にしなくてよい、と扱われがちです。しかし、そこで残った人が元の研究集団と同じとは限らず、欠測をどう扱ったかはハザード比と95%信頼区間の両方を動かします。
生存時間解析の欠測では、共変量だけを補完すればよいと考えるのも注意が必要です。イベントの有無と追跡時間が補完に使われなければ、予後と関連する情報を捨てたまま値を作ることになります。ここではsurvivalパッケージのpbcデータセットを使用して、完全ケース解析がどの仮定で成り立つか、多重代入で何を混ぜるか、そして推定値の差をどう報告するかを確かめます。また、Cox回帰用の補完モデルへイベントと追跡時間を入れる理由、完全ケースと多重代入の差を恣意的に選ばず報告する仕方も、同じ例に沿って扱います。
最初に欠けるのは値ではなく対象集団
欠測がある行を除く完全ケース解析は、最も手軽です。ただし、これは解析対象者を変える操作でもあります。たとえば胆汁酸塩の高い患者ほど検査を受けやすい、あるいは重症で早くイベントが起きた患者ほど来院時の検査値が残っていない、という状況なら、欠測のない人だけのCox回帰は元の登録集団をそのまま代表しません。
欠測機構は、観測値とも未観測値とも無関係なMCAR、観測されている情報で説明できるMAR、未観測値そのものにも依存するMNARに分けて考えます。そして、完全ケース解析が係数推定で比較的扱いやすいのは、欠測の起こり方がアウトカムと独立という強い条件が満たされる場面です。欠測割合が小さいことはこの条件の代わりにならず、Sterneらは、欠測割合だけを見て方法を決めるのではなく、なぜ値が欠けたかを研究の情報から考える必要があると整理しています。
私なら、欠測率の表を最初に作り、次に欠測あり・なしでイベント率と追跡時間の分布を比べます。差が見えたからMNARと結論づけるわけではありませんが、完全ケースだけで済ませる根拠が弱いことは分かります。情報的打ち切りで扱ったように、観測されなくなる仕組みを無視して同じ集団とみなすと、推定対象がいつの間にか変わります。
完全ケース解析が使える条件
完全ケース解析は検査項目の欠測がごく少なく、欠測の理由がイベントや予後と関係しないと説明でき、外した人数とイベント数を明記できるなら、主解析としても読みやすい方法です。ただし「欠測が5%未満なので完全ケースにした」という説明だけは適切ではありません。5%という数字から、欠測機構は分からないためです。
Cox回帰では、共変量が1つ欠けただけでその患者の追跡時間全体が尤度から外れます。イベントが114件ある284人で推定したハザード比と、イベントを保持した312人で推定したハザード比は、同じモデル式でも別の情報に基づく値です。特に候補共変量を多く入れる予後モデルでは、各変数の小さな欠測が重なり、完全ケース数が予想より大きく減ります。予測モデルの検証で説明した性能評価も、完全ケースだけに限定すれば、欠測を含む実際の対象集団への一般化可能性を別に考える必要があります。
一方で、アウトカムのイベント判定そのものが欠ける場合は話が異なります。共変量欠測の多重代入を、そのままイベント欠測に当てはめることは適切ではありません。イベント確認の手順、打ち切りの定義、そして未確認の理由を先に分けて考える必要があり、エンドポイントの定義と整合させます。
多重代入で推定値を一つに戻す
多重代入では、観測された年齢、検査値、イベント、追跡時間などから欠測値のもっともらしい候補を作り、補完済みデータセットを複数用意します。各データセットに同じCox回帰を当てはめ、係数を平均するだけでは足りません。補完ごとに推定値が動く不確実性も分散へ加えます。
\[
\bar{Q}=\frac{1}{M}\sum_{m=1}^{M}Q_m,\qquad
T=\bar{U}+\left(1+\frac{1}{M}\right)B
\tag{1}
\]
$Q_m$は$m$番目の補完済みデータセットの係数、$\bar{Q}$はその平均、$\bar{U}$は各データセット内の分散の平均、$B$は補完同士の係数のばらつき、$T$は最終的な分散です。つまり、補完後の標準誤差は、通常のCox回帰の不確実性だけでなく、欠測値をどのように補ったかによる不確実性も含みます。
単一代入では、中央値を1回入れた後のデータを本物の観測値と同じように扱います。この方法は標準誤差を小さく見せやすいため、推定値が欲しいだけの探索的な確認を除き、最終報告には使いません。MICEのような連鎖方程式による多重代入は、連続値、二値、順序尺度が混在する研究データへ使いやすく、Whiteらが実務上の注意点をまとめています。
補完モデルに生存時間を入れる理由
予後と関係する共変量を補完するなら、イベント指標だけでなく追跡時間も補完モデルへ入れます。生存時間解析では、同じイベントなしでも30日で打ち切られた人と3,000日追跡された人は同じ情報量ではなく、したがってイベント指標だけを入れると、この違いを補完の予測に使えません。
Cox回帰の係数を主な関心とする場合、Nelson-Aalen累積ハザード推定量とイベント指標を補完モデルに含める方法がよく使われます。そして、追跡時間を生のまま入れるか、累積ハザードとして入れるかは、欠測変数の型や補完手法で変わります。注意点として、解析モデルに入る変数だけを補完するだけでは十分ではなく、欠測値をよく予測し、かつ欠測の起こり方と関係しそうな補助変数も加える方が、MARの仮定を置く際の説明が具体的になります。
ただし、補完モデルへ大量の変数を機械的に入れればよいわけでもありません。観測数が少なく欠測が多い変数や、明らかに不安定な予測式を混ぜると、補完値そのものが不自然になります。連続変数なら補完値の範囲、二値変数なら割合、補完前後の分布を確認します。Coxモデルの診断と変数選択と同じく、モデルを一度当てはめて数字だけを受け取る進め方には無理があります。
pbcデータで284人と312人を比べる
survivalパッケージのpbcデータセットから治療群が記録された312人を取り出し、死亡をイベントとして、対数変換したビリルビン、アルブミン、コレステロールをCox回帰に入れます。コレステロールには28人の欠測があり、3変数が全てそろう完全ケースは284人、死亡は114件でした。そして、予測平均マッチングでコレステロールを20回補完し、各回のモデルをRubinの規則で統合します。
set.seed(42)
library(survival)
library(mice)
data(pbc, package = "survival")
analysis_data <- subset(pbc, !is.na(trt))
analysis_data$death <- as.integer(analysis_data$status == 2)
analysis_data$log_bili <- log(analysis_data$bili)
variables <- c("time", "death", "log_bili", "albumin", "chol", "age", "edema")
complete_data <- analysis_data[complete.cases(analysis_data[, variables]), variables]
complete_fit <- coxph(Surv(time, death) ~ log_bili + albumin + chol,
data = complete_data)
imputation <- mice(analysis_data[, variables], m = 20, maxit = 10,
method = c("", "", "", "", "pmm", "", ""), printFlag = FALSE)
imputed_fits <- with(imputation,
coxph(Surv(time, death) ~ log_bili + albumin + chol))
pooled_fit <- pool(imputed_fits)
round(c(n_complete = nrow(complete_data), events_complete = sum(complete_data$death),
hr_complete = exp(coef(complete_fit)["log_bili"]),
hr_mi = exp(summary(pooled_fit, conf.int = TRUE)[1, "estimate"]),
ci_mi_low = exp(summary(pooled_fit, conf.int = TRUE)[1, "2.5 %"]),
ci_mi_high = exp(summary(pooled_fit, conf.int = TRUE)[1, "97.5 %"])), 3)#> n_complete events_complete hr_complete.log_bili
#> 284.000 114.000 2.898
#> hr_mi ci_mi_low ci_mi_high
#> 2.948 2.378 3.656対数ビリルビンのハザード比は完全ケースで2.898、多重代入後で2.948でした。今回の差は大きくありませんが、「近いから完全ケースでよい」と読むのは早計です。この一致はこの変数、欠測割合、補完モデルでの結果にすぎません。また、多重代入後の95%信頼区間2.378〜3.656には、補完間のばらつきも含まれます。
log_bili_hr <- sapply(imputed_fits$analyses,
function(fit) exp(coef(fit)["log_bili"]))
png("article29-mi-log-bili-hr.png", width = 1100, height = 688, res = 200)
plot(seq_along(log_bili_hr), log_bili_hr, pch = 16,
ylim = range(c(log_bili_hr, exp(coef(complete_fit)["log_bili"]))),
xlab = "Imputed data set", ylab = "Hazard ratio for log(bilirubin)")
abline(h = exp(coef(complete_fit)["log_bili"]), lty = 2, col = "gray40")
abline(h = exp(summary(pooled_fit)[1, "estimate"]), lty = 1, col = "#C23B22")
legend("topright", legend = c("complete-case", "pooled multiple imputation"),
lty = c(2, 1), col = c("gray40", "#C23B22"), bty = "n")
dev.off()
差の大きさではなく仮定を報告する
完全ケースと多重代入でハザード比が大きく異なる時、どちらか都合のよい値を採用する進め方は推奨しません。まず、完全ケースで誰が外れたか、補完モデルへ何を入れたか、MARを支持する観測情報が何かを確認し、そのうえで主解析を先に定め、もう一方を感度分析として併記します。推定値が近ければ頑健性を確認する材料になり、また、離れていれば欠測の仮定が結論に関与していると分かります。
MARは観測データだけから検定できまないため、欠測値が観測されない値にも依存するMNARの可能性が大きい研究では、MARの多重代入だけで結論を出さず、欠測値を一定量だけ悪化させるdelta調整や、パターン混合モデルによる感度分析も検討します。この部分は動的治療戦略のように、推定対象と仮定を先に言葉にしておくと、解析の選択が追いやすくなります。
論文には、変数ごとの欠測数、完全ケース数とイベント数、補完回数、補完法、補完モデルの変数、各推定値と95%信頼区間を書きます。「多重代入を行った」とだけ書いても、読者はMARをどの情報で支えたのか確認できません。ここまで読んで、自分の研究の欠測が完全ケースで済むのか迷われた場合は、Dr.データサイエンスで、欠測の理由から補完モデル、感度分析まで含めた解析方針をご相談いただけます。
参考文献
- 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.
- White IR, Royston P, Wood AM. Multiple imputation using chained equations: issues and guidance for practice. Statistics in Medicine. 2011;30(4):377-399.
- White IR, Royston P. Imputing missing covariate values for the Cox model. Statistics in Medicine. 2009;28(15):1982-1998.
- Carpenter JR, Kenward MG, White IR. Sensitivity analysis after multiple imputation under missing at random: a weighting approach. Statistical Methods in Medical Research. 2007;16(3):259-275.
- Moons KGM, Donders RART, Stijnen T, Harrell FE Jr. Using the outcome for imputation of missing predictor values was preferred. Journal of Clinical Epidemiology. 2006;59(10):1092-1101.
前回：情報的打ち切り ／ 次回：感度分析の設計と報告
