予後予測モデルを組む場面では、Coxモデルの代わりにRandom Survival Forestのような手法を使う選択肢があります。木を多数育てて平均する手法で、非線形な関係も変数間の交互作用も、あらかじめモデル式に書かなくても拾えます。
Coxモデルより当たる場面は確かにありますが、ただし条件つきです。リスクの形が線形に近ければCoxが上に来ますし、イベント数が足りなければどちらの手法も伸びません。どこで順位が入れ替わるのかを、乱数で生成したデータセットで確かめました。
木を育てて平均する仕組み
Random Survival Forestは、生存時間データ用に作られた決定木を多数育てて平均する手法です。Ishwaranらが2008年にAnnals of Applied Statisticsで提案しました。通常の回帰木と違うのは分割の基準で、ログランク検定の統計量を使います。各ノードで、候補となる変数と分割点の組み合わせのうち、2つの子ノードのログランク統計量が最大になるものを選びます。生存曲線が最もはっきり分かれる位置で切るということです。
終端ノードに落ちた患者の予測は、そのノードの症例から計算したNelson-Aalen推定量になります。1本の木では不安定なので、ブートストラップ標本ごとに木を育て、各分割では変数の一部だけを候補にして、最後に平均します。
$$\hat{H}_b(t \mid \mathbf{x}) = \sum_{t_l \le t} \frac{d_{l,b}}{n_{l,b}}, \qquad \hat{H}(t \mid \mathbf{x}) = \frac{1}{B} \sum_{b=1}^{B} \hat{H}_b(t \mid \mathbf{x}) \tag{1}$$
式(1)の$\hat{H}_b$は$b$番目の木が返す累積ハザードで、$d_{l,b}$と$n_{l,b}$は、その木で$\mathbf{x}$が落ちた終端ノードにおける時点$t_l$のイベント数とリスク集合の大きさです。$B$は木の本数で、森の予測はこの単純平均になります。なお、比例ハザード性も線形性も、この手続きのどこにも出てきません。
非線形成分が2割を超えると順位が逆転する
共変量を8個、いずれも標準正規分布から発生させ、対数ハザードを次の形にしました。
$$\log h(t \mid \mathbf{x}) = \log h_0(t) + 0.95\left(\sqrt{1-\alpha}\,\frac{\eta_{\mathrm{lin}}}{s_{\mathrm{lin}}} + \sqrt{\alpha}\,\frac{\eta_{\mathrm{non}}}{s_{\mathrm{non}}}\right) \tag{2}$$
式(2)の$\eta_{\mathrm{lin}}$は主効果だけの線形和（$0.60x_1 + 0.50x_2 + 0.40x_3 + 0.30x_4$）、$\eta_{\mathrm{non}}$は$x_1 x_2 + 0.65(x_3^2 – 1)$で、交互作用と2次項だけでできた部分です。$s_{\mathrm{lin}}$と$s_{\mathrm{non}}$はそれぞれの標準偏差で、これで割ることにより、$\alpha$を0から1まで動かしても線形予測子の標準偏差が0.95に保たれます。
訓練データは1,500例で、イベント数は約810件です。別に生成した4,000例で性能を測り、この手順を20回繰り返しました。真のリスクを使った場合のC統計量は、$\alpha$が0でも1でも0.71から0.73の間に収まっており、設計どおり難易度は変わっていません。
| 非線形成分の割合 | 主効果Cox | 森 | 森−Cox | 真のリスク |
|---|---|---|---|---|
| 0.0 | 0.725 | 0.712 | -0.013 | 0.726 |
| 0.2 | 0.705 | 0.709 | 0.004 | 0.725 |
| 0.4 | 0.677 | 0.702 | 0.025 | 0.720 |
| 0.6 | 0.644 | 0.693 | 0.049 | 0.716 |
| 0.8 | 0.601 | 0.689 | 0.089 | 0.714 |
| 1.0 | 0.506 | 0.674 | 0.168 | 0.709 |
森が主効果だけのCoxを上回ったのは、非線形成分が2割を超えたあたりからでした。割合0.2での差は0.004で、20回の反復での標準偏差が0.005ですから、ここは実質的に互角です。割合0.6で0.049、1.0で0.168まで開きます。
逆に、真のリスクが完全に線形の時は森が0.013下回りました。これは20回全てで同じ向きの差です。木は連続的な直線関係を階段で近似するので、その分だけ損をします。では、この0.013という差を実務でどう受け止めればよいでしょうか。臨床予測モデルのC統計量としては小さい差ですが、線形の関係しかないデータで森を使う理由にはなりません。
スプラインを入れたCoxが森に並んだ
ここまでの比較において、主効果だけのCoxは非線形性を書き込む機会を与えられていません。そこで、8変数それぞれに自由度3の自然スプラインを当て、さらに全28通りの2変数交互作用を入れたCoxモデルを、同じデータに当てはめました。パラメータ数は52になります。
結果は表の森の列を上回りました。非線形成分の割合0.4で0.712（森は0.702）、0.8で0.704（森は0.689）、1.0で0.699（森は0.674）です。真のリスクによる上限0.709との差は0.010しかありません。
この結果は、森が拾えた構造をCoxの枠内でも書けたことを意味します。式に書く手間を引き受ければ、セミパラメトリックなモデルのまま同じ精度に届いた、ということです。もちろん今回は真の構造が2変数の交互作用と2次項という、スプラインと交互作用項で表現しやすい形でした。真の構造が8変数の高次の組み合わせなら、この52パラメータのモデルでは追いつかないはずです。
私なら、森がCoxを明確に上回った時、それを最終モデルの選択とは受け取りません。データのどこかに非線形構造があるという情報が出たと考えて、部分依存プロットで形を見てから、Coxに項を足して組み直せないかを先に試します。係数と信頼区間が出るモデルを、性能を落とさずに保てるなら、その方が論文でも臨床でも扱いやすくなります。
イベント100件では真のリスクに届かない
非線形成分の割合を0.7に固定し、訓練データの大きさだけを200例から3,200例まで変えました。イベント数にすると108件から1,731件です。検証データは同じく4,000例です。
イベント108件の時点では、森は0.649、スプラインと交互作用を入れたCoxは0.650で、ほとんど差がありません。真のリスクによる上限は0.713ですから、どちらも0.06以上離れています。イベントが1,731件まで増えると森は0.700、スプライン入りCoxは0.711になり、上限との差はそれぞれ0.016と0.005まで縮みました。
柔軟な手法ほど多くのイベントを必要とする、という傾向は以前から報告されています。van der Ploegらは2014年のBMC Medical Research Methodologyで、ランダムフォレストやサポートベクターマシンが安定したAUCに達するには、ロジスティック回帰の10倍を超える変数あたりイベント数を要する場合があると書いています。二値アウトカムを対象にした報告ですが、今回の生存時間データでも同じ向きの結果になりました。
では、症例数を増やせばどの手法も上限に近づくのでしょうか。主効果だけのCoxは、イベント108件で0.610、1,731件でも0.627で、早い段階で頭打ちになっていました。データを増やしても、モデルが書けない構造は書けるようになりません。
順位が当たっても3年生存率は外れる
C統計量が高ければ予測モデルとして使える、と考えたくなりますが、これは正確ではありません。C統計量は誰が先にイベントを起こすかという順位付けだけを測っており、「3年生存率は72%です」という予測値が実際の72%に対応しているかは見ていないためです。この対応が較正で、詳しくは予測モデルの検証で扱っています。
非線形成分の割合0.7、訓練1,500例のモデルで、検証8,000例を予測3年生存率の十分位に分け、各群の平均予測値とKaplan-Meier推定値を比べました。
主効果だけのCoxでは、予測3年生存率が最も高い十分位で予測値0.874に対して実測0.787でした。8番目の群の実測0.816より低く、順位が逆転しています。信頼区間は重なるので、この逆転自体を確定した現象とまでは言えませんが、モデルが本当に予後の良い患者を選び分けられていない状態と整合します。森の方は10群全てで実測が単調に上がりました。
ただし森も較正が合っているわけではありません。予測値が最も低い十分位では、予測0.326に対して実測0.208で、0.118のズレがあります。予測値と実測値の差の平均は、Coxで0.052、森で0.037でした。順位付けが良い手法が、そのまま確率としても正しいとは限りません。
絶対リスクを臨床で使う予測モデルなら、C統計量だけで手法を決めるのは足りません。
1,500例のデータで両方を当ててみる
ここまでの数値を出した手順のうち、中心になる部分を載せます。survivalパッケージとrandomForestSRCパッケージ（バージョン3.5.1）を使います。まずデータを生成する関数です。
library(survival)
vn <- paste0("v", 1:8)
gen <- function(n) {
x <- matrix(rnorm(n * 8), n, 8); colnames(x) <- vn
lin <- 0.60*x[,1] + 0.50*x[,2] + 0.40*x[,3] + 0.30*x[,4] # 線形成分
non <- x[,1]*x[,2] + 0.65*(x[,3]^2 - 1) # 交互作用と2次項
lp <- 0.95 * (sqrt(0.3)*lin/sqrt(0.86) + sqrt(0.7)*non/sqrt(1.845))
ev <- rexp(n, rate = 0.10*exp(lp)); ce <- rexp(n, rate = 0.05)
data.frame(tm = pmin(ev, ce, 12), st = as.numeric(ev <= pmin(ce, 12)), lp, x)
}非線形成分の割合を0.7にした場合です。この関数で訓練1,500例と検証4,000例を作り、両方の手法を当てはめます。
library(randomForestSRC)
set.seed(42)
train <- gen(1500); test <- gen(4000)
cox_fit <- coxph(Surv(tm, st) ~ ., data = train[, c("tm", "st", vn)])
rsf_fit <- rfsrc(Surv(tm, st) ~ ., data = train[, c("tm", "st", vn)],
ntree = 300, nodesize = 15, mtry = 3, seed = -42)
score_cox <- as.vector(as.matrix(test[, vn]) %*% coef(cox_fit))
score_rsf <- predict(rsf_fit, newdata = test[, c("tm", "st", vn)])[["predicted"]]
cc <- function(s) 1 - concordance(Surv(tm, st) ~ s, data = test)[["concordance"]]
round(c(events = sum(train[["st"]]), cox = cc(score_cox),
rsf = cc(score_rsf), oracle = cc(test[["lp"]])), 3)#> events cox rsf oracle
#> 801.000 0.624 0.700 0.720訓練データのイベントは801件で、森の0.700がCoxの0.624を0.076上回りました。真のリスクによる上限0.720との差は0.020です。
読み違えやすいのは最後の行の意味です。ここでのoracleは真の線形予測子をそのまま順位付けに使った場合の値で、どんな手法も原理的にこれを超えられません。0.720という数字は、この設定では打ち切りと偶然の分だけ順位付けに限界があることを表しています。
もう一点、concordance()の返り値をそのまま使わずに1から引いています。Rのconcordanceは共変量が大きいほど生存時間が長いという向きで計算するため、リスクスコアを入れると0.5より小さい値が返ってきます。
森を使ったと書くだけでは再現できない
2024年にBMJで公表された声明は、回帰でも機械学習でも同じ項目を報告するよう求めています。手法の名前だけでなく、ハイパーパラメータとその決め方、訓練と検証の分け方、そして較正の結果です。
今回の設定なら、木の本数300、終端ノードの最小サイズ15、各分割での候補変数3個、という値を書くことになります。これらを書かずに「Random Survival Forestを用いた」とだけ書いた論文は、読者が再現できません。Coxモデルなら係数表が事実上の再現手段になりますが、森にはそれがないので、報告の負担はむしろ重くなります。
較正の報告も同じです。図3の左端で見たように、順位付けが良い手法でも予測確率が0.1以上ズレることがあります。C統計量だけを報告して較正を省いた論文は、絶対リスクを使う読者にとって判断材料が足りません。
報告の負担が重くなる分だけ、森を選ぶ理由もはっきりしている必要があります。
私が森を使うのは構造を探す時
私が相談を受けた時にまず確認するのはイベント数です。イベントが100件前後のデータであれば、図2の左端の状況にあたるので、森でもスプライン入りのCoxでも上限には遠く、変数を絞って主効果のCoxを丁寧に組む方が結果は安定します。正則化Coxモデルで扱った罰則も、この段階で考える選択肢になります。
イベントが数百件以上あり、非線形性や交互作用が想定される場合には、森を当てる価値があります。ただし主解析としてではなく、Cox比例ハザードモデルに足すべき項を探す道具としてです。効果量とその不確実性を報告する必要がある研究で、森を主解析にする理由は今回の結果からは出てきませんでした。効果の不確実性に外部の情報を足す方向の話はベイズ生存時間解析で扱っています。
森は比例ハザード性を仮定しない、という点は長所として語られます。裏返せば、仮定を確かめる枠組みも持たないということです。Coxモデルなら比例ハザード仮定の検証と診断で扱ったSchoenfeld残差で違反の形まで見えますが、森ではモデルが何を学んだのかを確かめる標準的な手続きがありません。部分依存プロットやSHAP値は振る舞いの説明であって、予後因子の同定に使える推定量ではないと考えます。
手元のデータで森を使う価値があるのか、Coxで十分なのか。イベント数、候補変数の数、想定される構造、そして報告先が求める内容まで見ないと決まりません。Dr.データサイエンスでは、こうした予後予測モデルの設計と検証についてのご相談を承っています。
参考文献
Ishwaran H, Kogalur UB, Blackstone EH, Lauer MS. Random survival forests. Annals of Applied Statistics. 2008;2(3):841-860.
van der Ploeg T, Austin PC, Steyerberg EW. Modern modelling techniques are data hungry: a simulation study for predicting dichotomous endpoints. BMC Medical Research Methodology. 2014;14:137.
Christodoulou E, Ma J, Collins GS, Steyerberg EW, Verbakel JY, Van Calster B. A systematic review shows no performance benefit of machine learning over logistic regression for clinical prediction models. Journal of Clinical Epidemiology. 2019;110:12-22.
Collins GS, Moons KGM, Dhiman P, et al. TRIPOD+AI statement: updated guidance for reporting clinical prediction models that use regression or machine learning methods. BMJ. 2024;385:e078378.
Steyerberg EW. Clinical Prediction Models: A Practical Approach to Development, Validation, and Updating. 2nd ed. Springer; 2019.
