【医療統計版】正則化Coxモデル：候補変数を増やすと性能はどうなるか

電子カルテやオミクスのデータでは、候補変数が数十から数千ある一方で、イベント数はそれほど多くないことがあります。この状態でCoxモデルを組むと、訓練データにはよく当てはまり、C統計量も高く見えることがあります。

問題は、その性能が別の患者でも保たれるかどうかです。Coxモデルの診断と変数選択で扱ったように、同じデータを使ってモデルを作り、そのまま同じデータで性能を評価すると、結果は楽観的になります。候補変数が多いほど、この問題は目立ちやすくなります。

では、候補変数を増やしたとき、見かけの性能と別データでの性能はどのくらい離れるのでしょうか。ランダムに生成したデータセットで確認すると、かなり分かりやすい形になりました。

80変数ではC統計量の差が0.285になった

訓練データは150例で、イベント数は平均88件です。真に生存時間と関連する変数は6個だけに固定し、それ以外はアウトカムと関係のないノイズ変数としました。候補変数を6個から80個まで増やし、別に生成した2,000例で性能を評価します。この手順を100回繰り返した平均が次の表です。

数字の動きが逆なのが分かります。候補変数を6個から80個に増やすと、訓練データのC統計量は0.750から0.888まで上がりました。ところが別データでは0.728から0.604まで下がっています。

ノイズ変数には本来、予測に使える情報はありません。それでも訓練データの中だけを見れば、偶然イベントと関連しているように見える変数が出てきます。変数が増えるほど、その偶然をモデルが拾える余地も増えます。

較正勾配も同じ方向へ崩れました。6変数では0.886でしたが、80変数では0.068です。患者間のリスク差を、実際よりはるかに大きく予測している状態です。

予後モデルでは「1変数あたり10イベント」という経験則が長く使われてきました。ただし、現在では10を絶対的な境界として扱う方法は勧められていません。必要な症例数は、候補パラメータ数、予測性能、アウトカムの頻度、許容する過学習の程度などを考えて決めます。Rileyらの方法は、この考え方を症例数設計として整理したものです。

それでも今回の表には、イベント数に比べて候補変数を増やしたときに何が起こるかがよく出ています。変数80個のモデルでは、同じデータだけを見ればC統計量0.888ですが、別データでは0.604です。

Rでは見かけの性能と別データを分ける

上のシミュレーションでは、候補変数を最大80個生成しておき、Coxモデルへ入れる数だけを変えています。

library(survival)
set.seed(42)
n_tr <- 150; n_te <- 2000; n_true <- 6
beta <- c(0.55, 0.45, 0.40, -0.40, -0.35, 0.30)

gen <- function(n, p = 80){
  X  <- matrix(rnorm(n*p), n, p); colnames(X) <- paste0("v", 1:p)
  lp <- as.vector(X[, 1:n_true] %*% beta)
  tt <- rexp(n, rate = 0.08*exp(lp)); cc <- rexp(n, rate = 0.05)
  data.frame(tm = pmin(tt, cc), st = as.numeric(tt <= cc), X)
}
cidx <- function(tm, st, lp){
  z <- concordance(Surv(tm, st) ~ lp)$concordance
  if (z < 0.5) 1 - z else z
}

tr <- gen(n_tr); te <- gen(n_te)
for (p in c(6, 10, 20, 40, 60, 80)){
  vn  <- paste0("v", 1:p)
  fit <- coxph(as.formula(paste("Surv(tm, st) ~", paste(vn, collapse = "+"))),
               data = tr)
  cf <- coef(fit); cf[is.na(cf)] <- 0
  lp_tr <- as.vector(as.matrix(tr[, vn]) %*% cf)
  lp_te <- as.vector(as.matrix(te[, vn]) %*% cf)
  cat(p, round(cidx(tr$tm, tr$st, lp_tr), 3),
         round(cidx(te$tm, te$st, lp_te), 3),
         round(coef(coxph(Surv(tm, st) ~ lp_te, data = te)), 3), "\n")
}

最後に計算している係数が較正勾配です。訓練データで作った線形予測子を、別データのCox回帰へ1つの共変量として入れています。

理想的には1です。1より小さければ、モデルが作った線形予測子の広がりが別データでは大きすぎることを表します。表の80変数で0.068になっているのは、かなり強い過学習が起きている状態です。

実際の研究では真の係数を知ることも、自由に2,000例の検証データを作ることもできません。そのため、内部検証にはブートストラップや交差検証を使います。ここで重要なのは、モデル構築の手続きを検証の内側に入れることです。

罰則は膨らんだ係数を0へ戻す

候補変数をそのまま残しながら過学習を抑える方法の一つが正則化です。部分尤度に罰則を加え、データだけに合わせて大きくなった係数を縮めます。

$$\ell_{\textrm{pen}}(\boldsymbol{\beta}) = \ell(\boldsymbol{\beta}) – \lambda \sum_{j} |\beta_j|^q \tag{1}$$

式(1)では、$\lambda$が罰則の強さを決めます。$q=2$ならリッジ、$q=1$ならラッソです。$\lambda$を大きくするほど、推定された係数は0へ強く縮められます。

リッジでは係数を連続的に縮めますが、通常はちょうど0にはしません。ラッソでは一部の係数が0になるため、推定と同時に変数を減らす働きがあります。

リッジとラッソの中間にあたるのがエラスティックネットです。相関の強い候補変数が多いデータでは、純粋なラッソより安定して変数群を扱える場合があります。ただし、「相関する変数を必ず一緒に残す」とまで考えるのは正確ではありません。

今回のシミュレーションでは、罰則なしの60変数モデルで、真の係数0.55に設定した変数が平均0.949まで膨らんでいました。イベント数に対して候補変数が多く、訓練データの偶然まで係数が拾っているためです。

正則化は、この膨らみを抑えます。

coxph(Surv(tm, st) ~ ridge(v1, v2, v3, theta = 50), data = tr)

ただし、係数を縮めれば元の情報量が増えるわけではありません。ここを取り違えると、「症例数が少なくても正則化すれば何とかなる」という使い方になります。

60変数のデータで罰則の強さを変えたところ、今回のシミュレーションでは別データのC統計量が0.639から最大0.660まで上がりました。改善幅は0.021です。

6変数だけを入れたモデルでは別データのC統計量が0.728でした。同じデータ生成条件で比べると、ノイズ変数を54個追加してから罰則で縮めても、最初から真の6変数だけを使った状態までは戻っていません。

もちろん、実データではどれが真の6変数か分かりません。ここで言いたいのは、罰則が「候補変数を無制限に増やしてよい理由」にはならないということです。

高次元データでは罰則の意味が変わる

遺伝子発現量のように、候補パラメータ数が症例数やイベント数に比べて非常に多い場合は状況がさらに厳しくなります。

罰則なしのCox回帰では、係数を一意かつ安定して推定するための情報が不足し、収束しなかったり、極端な係数が得られたりすることがあります。特に高次元データでは、正則化は性能を少し改善するための追加的な工夫というより、モデルを推定可能な形に制約する役割を持ちます。

このとき訓練データのC統計量は、ほとんど参考になりません。数千、数万個の候補変数からモデルを作れば、訓練データの偶然にかなり細かく適合できます。性能を評価するなら、モデル構築から完全に分離したデータが必要です。

遺伝子発現量では変数同士の相関も強くなります。ラッソでは、似た情報を持つ変数の中でどれが残るかがデータのわずかな違いで変わることがあります。

そのため、ラッソで残った遺伝子をそのまま「重要な予後遺伝子」と解釈するのは適切ではありません。予測変数として選ばれたことと、生物学的に特別な役割を持つことは別の主張だからです。

較正も見ておく必要があります。今回の60変数データでは、罰則なしの較正勾配は0.186でした。罰則を強めると0.751まで改善しましたが、さらに強くすると1を超え、罰則の強さ100付近では1.107、500では3.438になりました。

罰則を強くしすぎると、今度は予測値の広がりを縮めすぎます。C統計量は順位を見る指標なので、係数を全体に縮めても大きく変わらないことがあります。一方、較正勾配は予測値の広がりに反応します。

したがって、$\lambda$を交差検証で決めたあとも、較正まで別に確認します。臨床で絶対リスクを使う予測モデルなら、C統計量だけで正則化の強さを評価するのは足りません。

選択された変数より選択頻度を見る

変数選択には、もう一つ見えにくい問題があります。

変数20個、イベント約90件のデータからブートストラップ標本を200個作り、それぞれで後ろ向き選択を行いました。真に関連する6変数は89%から100%の標本で残りましたが、ノイズ14変数も平均40%の標本で選択されています。最も高いノイズ変数では57%でした。

1回の解析結果だけを見れば、その変数が真の予後因子なのか、たまたま今回のデータで残ったのかは分かりません。

こういう場合、最終的に選ばれた変数名だけを報告するのではなく、ブートストラップでの選択頻度を併記すれば、少なくとも選択結果がどの程度安定しているかは見えます。

ただし、選択頻度にも注意が必要です。相関の強い変数が3つあり、どれでも同程度にアウトカムを説明できる場合、3変数が交代で選ばれて個々の頻度は低くなることがあります。頻度が低いから情報を持たない、と単純には判断できません。

モデルを作る前に減らせる候補は、統計的な選択に回さず先に減らします。既往の研究、測定の再現性、欠測、実装時に取得できるかどうかなどから、研究目的に不要な候補を落とします。

同じ情報を持つ変数が多ければ、要約スコアや主成分を使う選択肢もあります。ただし主成分はアウトカムを見ずに作るため、予後に関係する方向と一致するとは限りません。自由度を減らせることと、予測性能が上がることは分けて考えます。

それでも候補変数が多い場合に正則化を使います。$\lambda$を交差検証で決めるなら、性能評価はさらに外側で行います。内側で$\lambda$を選び、外側でその一連のモデル構築手順を評価するネストされたクロスバリデーションの形です。

同じ交差検証結果を使って$\lambda$を選び、その値をそのまま最終的な性能として報告すると、モデル選択に使った情報が性能評価にも入ります。

論文では候補数と検証手順まで書く

正則化Coxモデルを報告するとき、最終的に係数が残った変数数だけを書いても、過学習の程度は分かりません。最初に何個の候補パラメータからモデル構築を始めたのかと、イベント数の両方を書きます。

「最終モデルは8変数だった」という情報と、「候補10変数から8変数になった」のか「候補2,000変数から8変数になった」のかでは、読者が受け取る意味がまったく違います。

性能は、モデル構築に使わなかった情報で評価した値を報告します。同じデータで計算したC統計量しかない場合、図1の上側の線を見ていることになります。

症例数が限られる予測モデルでは、単純に訓練群と検証群へ1回分割すると、モデル開発に使えるイベント数まで減ります。そのため、内部検証としてはブートストラップを使う方法が有力です。検証方法そのものは予測モデルの検証で扱いました。

罰則を使った場合は、リッジ、ラッソ、エラスティックネットのどれを使ったのか、$\lambda$をどう決めたのか、その選択を内部検証のどこに含めたのかを書きます。

較正も外せません。正則化の強さを変えてもC統計量がほとんど動かない一方で、較正勾配が大きく変わることは図3からも分かります。絶対リスクを出す予測モデルなら、識別能と較正を別々に評価します。

係数の読み方にも注意します。ラッソやリッジの係数は、罰則によって0方向へ縮めることを意図した推定値です。通常の最尤推定によるハザード比と同じ感覚で、「この変数の効果量はHR 1.4」と因果的・説明的に解釈する目的には向きません。

予測が目的なのか、個々の変数とアウトカムの関連を推定したいのか。ここを曖昧にしたまま正則化を使うと、モデル構築の目的と係数の解釈がずれます。

繰り返し測定された検査値を候補に含める場合には、別の問題もあります。観測値の測定誤差によって関連が薄く見えることがあるため、縦断データと生存データの結合モデルで扱った測定誤差も切り分けます。木を使った予測手法との違いは機械学習的生存時間解析で扱っています。

実際の解析では、まずイベント数と候補パラメータ数を見ます。そのうえで、事前に減らせる変数がないかを確認し、残った候補に正則化が必要かを判断します。罰則を入れることより、モデル構築から性能評価までを一つの手順として検証することの方が重要です。

候補変数をどこまで残すか、内部検証をどう組むか、外部検証用のデータを用意できるかによって解析方針は変わります。Dr.データサイエンスでは、こうした予後予測モデルの設計と検証についてのご相談を承っています。

参考文献

1. Tibshirani R. The lasso method for variable selection in the Cox model. Statistics in Medicine. 1997;16(4):385-395.
2. Simon N, Friedman J, Hastie T, Tibshirani R. Regularization paths for Cox’s proportional hazards model via coordinate descent. Journal of Statistical Software. 2011;39(5):1-13.
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. Steyerberg EW. Clinical Prediction Models. 2nd ed. Springer; 2019.
5. Harrell FE. Regression Modeling Strategies. 2nd ed. Springer; 2015.