予後予測モデルの論文に「C統計量0.75」と書かれていたとします。では、そのモデルを別の病院の患者に使ったときも0.75になるのでしょうか。
今回は、150例、イベント89件、候補変数15個という条件で乱数からデータセットを作り、同じモデルを3通りの方法で評価しました。検証方法を変えるだけで、数字がどのくらい動くのかを見てみます。
同じ数字が3つある
モデルを作ったデータで測ると0.749でした。一方、独立したデータでは0.631です。この差が、この条件で生じた過大評価です。
イベント89件に対して候補変数は15個なので、1変数あたりのイベント数は5.9件。正則化Coxモデルの記事で触れた目安の10を下回っており、過学習が起きても不思議ではない条件です。
ここで少し気になるのが、ブートストラップ補正後の0.703も、外部データで得られた0.631までは下がっていないことです。
ブートストラップで評価できるのは、基本的には「同じ集団のデータを使ってモデルを作ったことによる楽観度」で、病院が変わった、患者背景が変わった、測定方法が違う、といった集団間の違いまでは補正してくれません。
そのため、内部検証できちんと補正した数字であっても、別施設に持っていけばさらに性能が落ちることがあります。
楽観度をブートストラップで測る
ブートストラップでは、元データから復元抽出してB個の標本を作り、それぞれでモデルを最初から作り直します。
そして、そのモデルを「モデル作成に使ったブートストラップ標本」と「元データ」の両方で評価します。この2つの差が、その標本での楽観度です。
$$C_{\text{corrected}} = C_{\text{app}} – \frac{1}{B}\sum_{b=1}^{B}\left\{C\left(\hat{\beta}_b, D_b\right) – C\left(\hat{\beta}_b, D\right)\right\}$$
$C_{\text{app}}$はモデルを作ったデータで測ったC統計量、$\hat{\beta}_b$は$b$番目のブートストラップ標本$D_b$から推定した係数、$D$は元データです。
ここで大事なのは、「モデルを作り直す」という部分です。変数選択をしているなら、それもブートストラップのループの中で毎回やり直さなければいけません。
先に変数選択を済ませて、選ばれた変数だけを使ってブートストラップを回すと、変数選択によって生じた楽観度が補正されません。
今回の例では補正量は0.047でした。ただし、変数選択まで含めてモデル構築を毎回やり直せば、補正量はこれより大きくなる可能性があります。
1回の分割は当てにならない
データを訓練用と検証用に分ける方法は直感的で分かりやすく、実際によく使われます。ただし、症例数が多くない場合は、どの患者がどちらに入ったかで結果がかなり変わります。
そこで、同じ150例を70対30にランダム分割する操作を200回繰り返し、そのたびにC統計量を計算しました。
結果を見ると、かなり大きな差です。
ある分け方なら0.55、別の分け方なら0.81になり、どちらが論文に載ってもおかしくありません。1回だけの分割検証は、この分布の中からたまたま1点を選んで報告しているのに近いわけです。
もちろん、症例数が数千例もあれば振れ幅は小さくなりますが、問題になるのは、数百例程度、特にイベント数が3桁に届かないような研究です。
その規模なら、単純な1回の分割よりもブートストラップや交差検証を使った方が情報を無駄にしません。
交差検証も1回やれば十分、というわけではありません。たとえば10分割交差検証でも、どの症例をどのfoldに入れたかで結果は変わります。分割を変えて10回程度繰り返し、平均値とばらつきを示した方が結果を解釈しやすくなります。
順位が合っても確率が合わない
C統計量が評価しているのは、患者の順位付けです。
つまり、より早くイベントが起こる患者に高いリスクを付けられているかどうかを見ています。予測された生存確率そのものが正しいかどうかは、C統計量だけでは分かりません。
較正勾配は、外部データでモデルの線形予測子を1つの共変量としてCox回帰に入れたときの係数です。
理想は1です。1より小さい場合、モデルが予測するリスクの差が実際より大きすぎることを示します。
今回の較正勾配は0.433となり、かなり小さい値です。モデル上では患者間に大きなリスク差があるように見えていても、実際の差はそこまで大きくありません。
これは過学習でよく見られるパターンで、係数が必要以上に大きく推定されるため、高リスクの患者はさらに高く、低リスクの患者はさらに低く予測されます。
C統計量だけを見ていると、この問題には気づけません。
もう1つ、較正は評価する時点によって変わります。図3は15か月時点の較正です。1年時点ではよく合っているのに、5年時点では大きくずれる、といったこともあります。
比例ハザード仮定が成り立っていないモデルでは、時間とともに較正が崩れることがあります。そのため、臨床で実際に使う時点について個別に確認する必要があります。
係数を縮めても順位は変わらない
較正勾配が0.433なら、係数を小さくすれば改善するのでしょうか。
係数をどの程度縮めるべきかは、モデルを作ったデータから大まかに見積もる方法があります。モデルの尤度比カイ二乗を$\chi^2$、変数の数を$p$とすると、経験的な縮小係数は次のように書けます。
$$\hat{s} = \frac{\chi^2 – p}{\chi^2}$$
$\hat{s}$は各係数に掛ける倍率です。1に近ければほとんど縮小は不要で、小さくなるほど強い縮小が必要だと考えます。
今回のモデルではカイ二乗が67.6、変数が15個なので、$\hat{s}$は0.778でした。
ところが、外部データで実際に得られた較正勾配は0.433です。内部データから計算した0.778より、はるかに強い縮小が必要な状態でした。
考えてみれば当然で、縮小係数そのものもモデルを作ったデータから計算しています。そこにも楽観性が入り得ます。「縮小係数がそれほど小さくなかったから大丈夫」と判断するのは少々乱暴といえます。
試しに係数を0.778倍し、外部データでもう一度評価しました。
C統計量は0.631のままです。一方、較正勾配は0.433から0.556に改善しました。
これは、すべての係数を同じ割合で小さくしても、患者の順位は変わらないからです。順位が同じならC統計量も同じです。一方で、予測値の広がりは小さくなるので、較正は改善します。
ここからも、識別と較正は別の性能だと分かります。C統計量が上がらないからといって、縮小に意味がないわけではありません。
では、そもそもどのくらいのデータが必要だったのでしょうか。
縮小係数を0.9以上に保とうとすると、今回の条件ではカイ二乗が150程度必要になります。現在の67.6がイベント数にほぼ比例すると仮定すれば、必要なイベント数は約198件です。
しかし、実際には89件しかありません。つまり、今回の条件では15変数のモデルを安定して開発するには、最初からデータ量が不足していたと考えられます。
臨床的な便益では「全員治療」に負けることもある
識別と較正が分かったとしても、まだ1つ問題が残ります。
そのモデルを使うことで、実際の臨床判断が良くなるのかどうかです。
これを見る方法の1つが決定曲線分析です。横軸に治療を開始するリスク閾値、縦軸にその閾値でモデルを使った場合の正味の便益を取ります。
たとえば閾値0.20では、モデルの正味の便益が0.337、全員治療が0.354でした。
つまり、この閾値では予測モデルを使って患者を選別するより、単純に全員を治療する方が成績が良いという結果です。
一方、閾値0.50ではモデルが0.101、全員治療が−0.033となり、今度はモデルを使う方が大きく上回ります。
15か月時点のイベント割合は0.483です。そのため、かなり低い閾値を使う場合は、結局ほとんどの患者が治療対象になります。その状況では、わざわざモデルで患者を振り分けるメリットが小さくなります。
C統計量0.631という数字だけを見ても、この違いは分かりません。
新しいバイオマーカーを加えて「C統計量が0.02改善した」と報告されても、臨床で使う予定の閾値周辺で正味の便益が変わっていなければ、実際の意思決定はほとんど変わらないかもしれません。
モデルを使う閾値を先に考え、その範囲で本当に便益が増えるかを見る必要があります。
論文では何を報告すればよいか
図1では、同じモデルから0.749、0.703、0.631という3つの数字が出ました。
数字だけを書かれても、それがどの検証で得られたものか分からなければ解釈できません。
モデルを作ったデータで測ったのか、訓練・検証に分割したのか、ブートストラップで補正したのか、それとも完全に別のデータで外部検証したのか。まず、この区別を明記する必要があります。
分割検証を使った場合は、何回繰り返したかと、結果がどの程度ばらついたかも書いておきたいところです。1回しか分割していないのであれば、その点も明記した方がよいでしょう。
C統計量しか報告されていない論文では、分かるのは患者の順位付けの性能だけです。
そのモデルが出した20％や50％というリスクを、そのまま目の前の患者に当てはめてよいのかは判断できません。予測モデルとして評価するのであれば、較正勾配や較正プロットも必要です。
さらに、特定のリスク閾値を使って治療や検査を決めるのであれば、決定曲線まで示した方が臨床的な意味を判断しやすくなります。
モデル構築の手順も重要です。候補変数をどう決めたか、変数選択を行ったか、連続変数をどう扱ったか、欠測値をどう処理したか。そして、それらの処理のどこまでを内部検証の中で繰り返したか。
たとえば変数選択だけを先に済ませ、その後に内部検証している場合、報告された性能には選択時の楽観性が残っています。
TRIPOD声明では、こうした予測モデル研究で報告すべき項目が整理されています。投稿前のチェックリストとして使うと便利です。
Rで確認する
survivalパッケージのconcordance関数を使ってC統計量を計算し、ブートストラップで楽観度を補正します。
ci <- function(cf, d) { lp <- as.vector(as.matrix(d[, vs]) %*% cf) concordance(Surv(tm, st) ~ lp, data = d, reverse = TRUE)$concordance } optimism <- replicate(300, { i <- sample(nrow(train), replace = TRUE) fb <- coxph(fml, data = train[i, ]) ci(coef(fb), train[i, ]) - ci(coef(fb), train) }) c(apparent = ci(coef(fit), train), corrected = ci(coef(fit), train) - mean(optimism), external = ci(coef(fit), test))#> apparent corrected external #> 0.749 0.703 0.631次に、較正勾配と縮小係数を求めます。較正勾配は、外部データ上で線形予測子を共変量1つとしてCox回帰に入れれば計算できます。
lp_test <- as.vector(as.matrix(test[, vs]) %*% coef(fit)) slope <- coef(coxph(Surv(tm, st) ~ lp_test, data = test)) chisq <- 2 * (fit$loglik[2] - fit$loglik[1]) c(slope = slope, shrinkage = (chisq - length(vs)) / chisq)#> slope shrinkage #> 0.433 0.778この2つは似た数字に見えますが、推定しているものは違います。
較正勾配0.433は外部データで実際に観測した値です。一方、縮小係数0.778はモデルを開発したデータだけから計算した見積もりです。
0.778と0.433の差を見ると、内部データだけから過学習の程度を推測することには限界があると分かります。
もう1点、survivalパッケージのconcordance関数ではreverse引数にも注意が必要です。これを指定しないと向きが逆になり、今回なら0.631ではなく0.369が返ります。値だけを見ると解析ミスに気づきにくいので、確認しておいた方がよいところです。
次に確かめること
候補変数が多い場合の過学習の抑え方は正則化Coxモデル、変数選択そのものの不安定さはCoxモデルの診断と変数選択で扱いました。機械学習の手法を使った場合に較正がどう崩れるかは機械学習的生存時間解析にあります。競合リスクがある場合の予測はFine-Grayモデルが土台になります。
同じ施設のデータをどう分割しても、それは内部検証です。患者背景も測定方法も診療の流れも揃っているので、別の施設や別の時期のデータで確かめてはじめて、そのモデルが移植可能かどうかが分かります。外部検証の結果が悪かった場合も、モデルを捨てる必要はありません。切片や係数を新しい集団に合わせて再較正すれば使えることがあり、その手続きを報告すれば次の施設が同じことをできます。
検証の設計をどう組むか、外部データをどこに求めるか、較正が崩れたモデルをどう扱うか。このあたりは研究の規模と実装先を見ないと決まりません。Dr.データサイエンスでは、こうした予後モデルの設計についてのご相談を承っています。
参考文献
1. 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.
2. Collins GS, Reitsma JB, Altman DG, Moons KG. Transparent reporting of a multivariable prediction model for individual prognosis or diagnosis (TRIPOD). Annals of Internal Medicine. 2015;162(1):55-63.
3. Riley RD, Snell KI, Ensor J, et al. Minimum sample size for developing a multivariable prediction model: Part II – binary and time-to-event outcomes. Statistics in Medicine. 2019;38(7):1276-1296.
4. Van Calster B, McLernon DJ, van Smeden M, Wynants L, Steyerberg EW. Calibration: the Achilles heel of predictive analytics. BMC Medicine. 2019;17:230.
5. Vickers AJ, Elkin EB. Decision curve analysis: a novel method for evaluating prediction models. Medical Decision Making. 2006;26(6):565-574.
6. Steyerberg EW. Clinical Prediction Models. 2nd ed. Springer; 2019.
