予後予測モデルを作る場面を考えます。候補となる臨床検査値が20個あって、まず単変量で当てはめてp値が0.20を下回ったものを拾い、それをステップワイズ法にかけて5変数のモデルに絞る。C統計量は0.68でした。査読者からは「Please report internally validated performance.」と返ってきます。
この0.68は、同じデータで作って同じデータで測った値です。ブートストラップで補正すると0.61になります。しかも選ばれた5変数のうち3つは、真の効果がゼロの変数でした。
ここではCoxモデルの診断に使う残差を確認したうえで、変数選択が何を壊すのか、そして選んだモデルの性能をどう見積もるかを見ていきます。
打ち切りがあるときの残差
線形回帰の残差は観測値から予測値を引いたものですが、生存時間解析では打ち切り例の真のイベント時刻が分かりません。代わりに使うのがマーチンゲール残差です。
$$M_i = \delta_i – \hat{H}(t_i \mid \mathbf{x}_i) \tag{1}$$
式(1)は、その個体に実際にイベントが起きたかどうか（0か1）から、モデルが予測した累積ハザードを引いたものです。モデルが正しければ期待値はゼロになります。値は1が上限で、下は負の方向に限りなく伸びます。イベントが起きた個体は正、打ち切り個体は負に出やすい非対称な量です。
この残差の使いどころは、連続共変量をどの形でモデルに入れるかの確認です。共変量を1つも入れないモデルの残差を、その共変量に対してプロットして平滑化曲線を重ねます。曲線の形が、log-ハザードの上でその変数がどう効いているかを教えてくれます。
library(survival)
liver <- na.omit(pbc[pbc$id <= 312, c("time", "status", "bili", "age", "albumin")])
liver$death <- as.numeric(liver$status == 2)
res <- resid(coxph(Surv(time, death) ~ 1, data = liver), type = "martingale")
plot(liver$bili, res, xlab = "Bilirubin", ylab = "Martingale residual")
lines(lowess(liver$bili, res, iter = 0), lwd = 2)
数値でも確かめられます。ビリルビンをそのまま入れたモデルのAICは1142.5、対数変換すると1105.8。C統計量は0.811から0.831に上がります。1変数の入れ方を変えただけで、AICが37縮みました。
臨床検査値には対数正規に近い分布のものが多く、ビリルビン、CRP、フェリチン、クレアチニンあたりは対数変換したほうが当てはまりがよくなることが珍しくありません。モデルを組む前にこのプロットを1枚描いておくと、あとで診断に戻る手間が減ります。
1例で結果が変わっていないか
係数がごく少数の症例に引っ張られていないかは、その症例を抜いたときに係数がどれだけ動くかで測ります。全例について1つずつ抜いて再推定するのは重いので、スコア残差と情報行列から近似します。
$$\mathrm{dfbeta}_i \approx I(\hat{\boldsymbol{\beta}})^{-1} L_i \tag{2}$$
式(2)は、個体$i$を除いたときの係数の変化量の近似です。$L_i$はその個体のスコア残差ベクトルで、部分尤度の推定方程式にどれだけ寄与しているかを表します。これを標準誤差で割って無次元にしたものがdfbetasで、変数間で大きさを比べられます。
fit <- coxph(Surv(time, death) ~ log(bili) + age + albumin, data = liver)
db <- residuals(fit, type = "dfbetas")
plot(db[, 1], xlab = "Observation", ylab = "dfbetas")
abline(h = c(-2, 2) / sqrt(nrow(liver)), lty = 2)
この1例を除いてモデルを組み直すと、log(ビリルビン)のハザード比が2.70から2.90に動きます。7%の変化です。
ここで注意が要るのは、目安を超えた症例をどうするかです。除外してはいけません。dfbetasが大きいことはデータが間違っている証拠ではなく、その症例が推定に効いているという事実を示すだけです。私なら、まずその症例のカルテを確認できるかを共著者に聞きます。検査値の入力ミスや単位の取り違えなら訂正します。臨床的に説明のつく経過なら残したうえで、その症例を除いた感度分析を補足資料に載せます。除外して本解析だけを報告すると、都合の悪い症例を落としたことになります。
この$2/\sqrt{n}$という目安自体、根拠の強いものではありません。この例で超えたのは312例中14例と17例、いずれも5%前後です。dfbetasがおおむね対称に散らばる以上、この目安はどんなデータでも数%を拾います。超えた数を数えることに意味はなく、飛び抜けて大きい少数を個別に見るほうが役に立ちます。
ステップワイズ法が作る幻
候補変数が多いとき、どれを残すかを自動で決めたくなります。前向き選択、後向き除去、その組み合わせ。Rならstep関数で1行です。
何が起きるかを確かめます。200例、イベント75件、候補変数20個。このうち真に効果があるのは2つだけで、残り18個は生存時間と完全に無関係な乱数です。後向き除去を100回のシミュレーションで走らせます。以下は説明のための架空データです。
set.seed(1)
n <- 200; p <- 20
covars <- matrix(rnorm(n * p), n, p)
colnames(covars) <- paste0("v", 1:p)
risk <- 0.5 * covars[, 1] - 0.4 * covars[, 2] # 効くのはv1とv2だけ
event_time <- rexp(n, rate = 0.05 * exp(risk))
cens_time <- runif(n, 0, 30)
sim <- data.frame(tm = pmin(event_time, cens_time),
st = as.numeric(event_time <= cens_time), covars)
chosen <- step(coxph(Surv(tm, st) ~ ., data = sim),
direction = "backward", trace = 0)
attr(terms(chosen), "term.labels")#> [1] "v1" "v2" "v5" "v8" "v15"
この動きは多重比較の問題そのものです。18個の無関係な変数を1つずつ検定すれば、そのうちいくつかは偶然p値が小さくなります。選択の過程で使ったp値を、選択後のモデルでもそのまま報告することはできません。
結果として、残った係数は過大に推定され、信頼区間は実際より狭くなります。データを少し変えれば選ばれる変数も変わるので、モデル自体が不安定です。別の施設のデータで再現しようとしても、同じ変数が選ばれません。
AICを基準にしても、この問題は消えません。AICは検定の多重性を調整する道具ではなく、当てはまりと複雑さのバランスを測る指標です。今回の100回はすべてAIC基準の後向き除去で、それでも18個の無効な変数が2割の頻度で残りました。
だから、変数選択でモデルを決めることは勧めません。臨床的に交絡が疑われる変数と、既知の予後因子を、データを見る前にプロトコルへ書く。これが基本です。候補が多すぎて絞らざるをえないなら、選択そのものをブートストラップの内側に入れて、選ばれる頻度と補正後の性能を報告します。
見かけの性能から楽観主義を引く
C統計量は、2人の患者を取り出したときに予後の悪いほうを正しく高リスクと判定できる割合です。0.5でコイン投げ、1.0で完全な識別になります。
問題は、モデルを作ったデータでC統計量を測ると、必ず高めに出ることです。この上振れを楽観主義と呼びます。ブートストラップで見積もれます。復元抽出した標本でモデルを作り直し、その標本での性能と、元データに当てはめたときの性能の差をとる。これを繰り返して平均します。
refit <- function(dat)
step(coxph(Surv(tm, st) ~ ., data = dat), direction = "backward", trace = 0)
optimism <- replicate(60, {
boot <- sim[sample(nrow(sim), replace = TRUE), ]
model <- refit(boot)
cindex(model, boot) - cindex(model, sim)
})
mean(optimism)#> 見かけのC統計量 0.676
#> 楽観主義 0.064
#> 補正後のC統計量 0.6120.676が0.612になりました。参考までに、真のv1とv2だけを入れたモデルの見かけのC統計量は0.652です。ステップワイズ法で選んだモデルは、補正すると正解のモデルより性能が低くなっています。
ここで重要なのは、変数選択をブートストラップの内側に入れることです。選択済みのモデルを固定してブートストラップすると、選択で生じた上振れが補正されず、楽観主義を小さく見積もります。上のコードでrefit関数を毎回呼んでいるのはこのためです。
報告には、見かけの値と補正後の値を両方書きます。TRIPOD声明も、予測モデルの論文には内部検証の方法と結果を記載することを求めています。補正後の値だけを書くと、どれだけ縮んだのかが読者に分かりません。
作る前に決めておくこと
予後モデルを組むなら、変数はデータを見る前に決めます。イベント数と候補変数の数の比を先に計算し、比が小さいなら候補を減らすか、罰則付きの推定を検討します。連続変数は図1のような形で入れ方を確かめ、切って階級にしない。カテゴリ化は情報を捨てるうえ、切る位置をデータから決めると変数選択と同じ問題が起きます。
組んだあとは、影響点の確認、比例ハザード性の検証、そして内部検証。比例ハザード仮定の検証はSchoenfeld残差で扱います。ここまでの推定の枠組みはCox比例ハザードモデルと部分尤度、時間とともに値が変わる共変量を入れる場合は時変共変量が続きです。
候補変数が多くイベント数が限られた予後モデルで、どこまで変数を入れるか。内部検証をどう組んで論文にどう書くか。このあたりはデータの規模と報告先の要求を見ないと決まりません。Dr.データサイエンスでは、こうしたモデル構築のご相談を承っています。
参考文献
Harrell FE, Lee KL, Mark DB. Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Statistics in Medicine. 1996;15(4):361-387.
Steyerberg EW, Harrell FE, Borsboom GJ, et al. Internal validation of predictive models: efficiency of some procedures for logistic regression analysis. Journal of Clinical Epidemiology. 2001;54(8):774-781.
Collins GS, Reitsma JB, Altman DG, Moons KGM. Transparent reporting of a multivariable prediction model for individual prognosis or diagnosis (TRIPOD). Annals of Internal Medicine. 2015;162(1):55-63.
Therneau TM, Grambsch PM, Fleming TR. Martingale-based residuals for survival models. Biometrika. 1990;77(1):147-160.
Riley RD, Snell KI, Ensor J, et al. Minimum sample size for developing a multivariable prediction model. Statistics in Medicine. 2019;38(7):1276-1296.
