追跡期間に沿った生存確率を滑らかに推定したいのに、Weibull分布では曲線の形がデータに合わず、Kaplan–Meier推定量では段差が残ることがあり、こうした場面で候補になるのが柔軟なパラメトリック生存モデルです。ただし、「基準ハザードにスプラインを当てるモデル」とだけ覚えると何をモデル化しているのかが曖昧になるため、Royston–Parmarモデルが表す量、節点で変わる柔軟性、観察範囲外へ延ばすときの限界を分けて考えます。
この記事を読むと、比例ハザード尺度のRoyston–Parmarモデルが対数累積ハザードをどう表すかを説明でき、単純なパラメトリック分布やCox比例ハザードモデルとの違いを判断できます。さらに、スプラインの自由度をAICだけで決めず、推定曲線とデータが疎な期間を併せて見る理由が分かります。
累積ハザードを選ぶ理由
生存時間解析では、生存関数$S(t)$、ハザード関数$h(t)$、累積ハザード関数$H(t)$が同じ時間分布を異なる角度から表します。これらの関係に馴染みがなければ、先に生存関数とハザード関数を確認すると、以下のモデル式を追いやすくなります。
Royston–Parmarモデルの比例ハザード尺度では、時間$t$における累積ハザード$H(t)$の対数を、対数時間の制限付き三次スプライン$s$と共変量で表します。
$$\log H(t\mid\mathbf{x})=s\{\log(t);\boldsymbol{\gamma}\}+\mathbf{x}^{\top}\boldsymbol{\beta}.\tag{1}$$
式(1)は、$\mathbf{x}$を共変量ベクトル、$\boldsymbol{\beta}$をその係数、$\boldsymbol{\gamma}$をスプライン係数として、対数累積ハザードの時間変化を滑らかな曲線で記述するという意味です。
ここで大切なのは、基準ハザードそのものへ直接スプラインを当てるとは限らない点で、Royston and Parmar（2002）は、比例ハザード尺度では対数累積ハザードを、比例オッズ尺度では対数累積オッズをモデル化しました。本記事では前者だけを扱い、柔軟なパラメトリックAFTモデルなどへ範囲を広げません。
制限付き三次スプラインの役割
指数分布ではハザードが一定、Weibull分布では対数累積ハザードと対数時間の関係が直線になりますが、研究開始後にハザードが上がり、その後に下がるような時間分布では、この直線性が強すぎることがあります。Durrleman and Simon（1989）が扱った制限付き三次スプラインは、内部では曲率を持たせ、境界の外側では線形にする仕組みです。
節点は曲線を接続する位置であり、内部節点が増えるほど対数累積ハザードの形を細かく追えますが、イベントの少ない右端まで細かな曲率を許すと、少数の観測に形が左右されます。そのため、自由度は単なる調整つまみではなく、モデルが時間分布のどの程度の複雑さを許すかを決め、そこから得る生存確率やハザードの安定性にも影響します。
完全パラメトリックモデルの利点
Cox比例ハザードモデルは基準ハザードの形を特定せず、共変量効果を推定するセミパラメトリックモデルです。一方、Royston–Parmarモデルは時間軸を含む生存分布全体を有限個の係数で表すため、Jackson（2016）が説明するように完全尤度で推定でき、生存関数、累積ハザード関数、ハザード関数を滑らかな関数として得られます。
この性質は、特定時点の生存確率や時間とともに変わるハザードの形を報告したい研究で役立ち、Kaplan–Meier推定量は観測範囲内の分布を仮定せず要約するには向く一方、観測イベント時点で段階的に変わるため滑らかなハザードを直接返しません。両者の役割を同じにせず、モデルに依存しない曲線との一致も点検するのが自然です。
ただし、滑らかに描けることと観察範囲外でも信頼できることは別で、完全パラメトリックモデルの外挿は、端のスプライン形状とその先も同じ構造が続くという仮定に依存します。追跡期間の先を推定する目的なら、パラメトリック生存モデルで扱う外挿の判断も必要です。
AICだけでは決まらない自由度
今回の解析では、内部節点数$k$を0、1、2、3と変えた候補のAICが、それぞれ8,267.44、8,196.60、8,192.55、8,193.59で、候補内では$k=2$が最小、$k=3$との差は1.04でした。この数値は$k=2$が比較したモデルの中で最も小さい情報量規準を持つことを表しますが、曲線の臨床的な妥当性や外挿性能まで保証しません。
私なら、AIC最小という一行で自由度を確定せず、Kaplan–Meier曲線との対応、ハザードの形、イベントが疎になる右端、研究で結果を必要とする時間範囲を同時に見ます。差が小さい候補なら、わずかなAICの改善と曲線の複雑さのどちらを優先するかも説明に含めます。
実際、$k=0$のモデルと選択した$k=2$のモデルを比べると、365日の生存確率の差は$k=0$から$k=2$を引いて−0.0152でしたが、1,095日では0.0533、1,825日では0.0290でした。このデータでは差の向きが時点で変わっており、単純な形が生存確率を常に過大評価する、あるいは常に過小評価するとは読めません。
colonデータで形を確かめる
survivalパッケージのcolonデータセットから死亡までの記録929件を使い、死亡は452件、打ち切りは477件で、観測された追跡期間は23〜3,329日でした。このデータはMoertel et al.（1990）の結腸癌術後補助療法試験に由来しますが、ここでは全治療群をまとめ、治療効果や因果効果を推定しません。
set.seed(42)
library(flexsurv)
colon_death <- subset(survival::colon, etype == 2)
colon_death$event <- colon_death$status == 1
evaluation_times <- c(365, 1095, 1825)
fits <- lapply(0:3, function(k) {
flexsurvspline(Surv(time, event) ~ 1, data = colon_death,
k = k, scale = "hazard")
})
fit <- fits[[3]]
round(vapply(fits, AIC, numeric(1)), 2)#> [1] 8267.44 8196.60 8192.55 8193.59rp_s <- summary(fit, t = evaluation_times, type = "survival")[[1]]
rp_h <- summary(fit, t = evaluation_times, type = "hazard")[[1]]
rp_c <- summary(fit, t = evaluation_times, type = "cumhaz", ci = FALSE)[[1]]
km <- summary(survfit(Surv(time, event) ~ 1, data = colon_death),
times = evaluation_times, extend = TRUE)$surv
k0 <- summary(fits[[1]], t = evaluation_times,
type = "survival", ci = FALSE)[[1]]$est
round(data.frame(time = evaluation_times, survival = rp_s$est,
km = km, k0_minus_k2 = k0 - rp_s$est, hazard = rp_h$est,
hazard_lcl = rp_h$lcl, hazard_ucl = rp_h$ucl,
cumulative_hazard = rp_c$est), 6)#> time survival km k0_minus_k2 hazard hazard_lcl hazard_ucl cumulative_hazard
#> 1 365 0.914138 0.916039 -0.015160 0.000385 0.000326 0.000451 0.089773
#> 2 1095 0.673425 0.674667 0.053313 0.000354 0.000318 0.000392 0.395378
#> 3 1825 0.558537 0.564417 0.029039 0.000186 0.000165 0.000205 0.582434curve_times <- unique(round(exp(seq(log(min(colon_death$time)),
log(max(colon_death$time)), length.out = 80))))
curve_s <- summary(fit, t = curve_times, type = "survival", ci = FALSE)[[1]]
curve_h <- summary(fit, t = curve_times, type = "hazard", ci = FALSE)[[1]]
par(mar = c(4.3, 4.6, 1.2, 4.6))
plot(curve_times, curve_s$est, type = "l", lwd = 2,
xlab = "Follow-up time (days)", ylab = "Survival probability",
ylim = c(0, 1), col = "#176B87")
par(new = TRUE)
plot(curve_times, curve_h$est * 1000, type = "l", lwd = 2, lty = 2,
axes = FALSE, xlab = "", ylab = "", col = "#C44E52")
axis(4); mtext("Hazard per 1,000 person-days", 4, 3)
legend("bottomleft", c("Survival", "Hazard"), lwd = 2, lty = c(1, 2),
col = c("#176B87", "#C44E52"), bty = "n")
生存確率とハザードの読み分け
図1の青線は生存確率、赤の破線は1,000人日当たりハザードで、生存確率は時間とともに低下し続けますが、ハザードは約700日付近まで上昇した後に低下しています。生存確率の傾きと、その時点まで生存した人における瞬間的なイベント率を混同しないことが、曲線を読むための基本となります。
さらに、1,825日の推定ハザードは1日当たり0.000186、95%信頼区間は0.000165〜0.000205で、この値を個人の1,825日までの死亡確率とは解釈しません。同時点の累積ハザードは0.582であり、式$S(t)=\exp\{-H(t)\}$から得た生存確率0.559と一致します。
この一致は、ハザード、累積ハザード、生存確率が別々の結果ではなく、同じ当てはめから変換された量だと確認するために有用です。一方、治療群間の比較をしていないので、曲線の形から治療の有効性を論じることはできません。
公表研究で求められる説明
Royston and Parmar（2002）は、柔軟なパラメトリックモデルを予後モデリングと治療効果推定へ応用し、比例ハザードから非比例ハザードへ拡張できる枠組みを提示しました。Lambert and Royston（2009）は時間依存効果や予測量への発展を記述していますが、これらは自動的に備わる解釈ではなく、どの尺度とモデル項を選んだかに依存します。
論文で「柔軟なパラメトリックモデルを用いた」とだけ書くと再現に必要な情報が不足するため、私は、比例ハザードなどの尺度、スプラインの自由度または節点とその配置、共変量、時間依存効果の有無、モデル選択の根拠、報告した予測時点、使用ソフトウェアと版を記載することを勧めます。この推奨は、同じモデルを組み直して予測を確かめられるようにするためで、今回の解析はR 4.4.2、flexsurv 2.3.2、survival 3.8.3で実行しました。
AICだけで候補を選び、滑らかな図が得られたことを適合の証明にする進め方は適切ではありません。モデルに依存しない推定との対応と、研究目的に必要な期間での挙動が欠けると、柔軟性が判断材料を増やすのではなく、仮定を見えにくくするためです。
観察範囲外では仮定が前に出る
今回の最長追跡は3,329日ですが、右端まで同じ密度で情報があるわけではなく、Kaplan–Meier曲線の右端が少数例に左右される仕組みは、Kaplan–Meier推定量でも扱っています。スプラインで線が滑らかになっても、右端の情報量が増えるわけではありません。
観察範囲を越えると、制限付き三次スプラインは境界節点の外側で線形になる性質に強く支配されるため、長期予測が研究目的なら、自由度の異なるモデルや標準的なパラメトリック分布でも外挿を比べ、外部資料と矛盾しないかを見る必要があります。今回のAIC最小モデルだけを3,329日より先へ延ばして採用する根拠は、解析結果にはありません。
治癒を想定する研究では、長期の平坦部を通常のRoyston–Parmarモデルだけで表すかどうかも別の問題です。その目的には治癒率モデルという異なる構造があり、本記事の範囲には含めていません。
研究目的からモデルを選ぶ
柔軟なパラメトリック生存モデルを選ぶ中心的な理由は、複雑な手法を使うことではなく、時間分布全体を滑らかに表し、研究で必要な時点の生存確率やハザードを得ることです。群間の相対効果だけが目的ならCox比例ハザードモデルが候補に残り、観察範囲内の生存曲線を仮定少なく記述するならKaplan–Meier推定量が基準になります。
モデル選択では、何を推定して報告するか、どの時間範囲に判断が必要か、スプラインの形を支えるイベントがあるかを先に結び付けます。その上でAICと推定曲線を読み、自由度を変えたときに結論が動くなら、その動き自体を不確実性として扱う方が誠実です。
実際の研究データでは、節点配置と予測期間の判断が解析結果を左右します。自分の研究でどの尺度と自由度を選ぶべきか、観察期間を越えた予測をどこまで報告できるかを整理したい場合は、Dr.データサイエンスで解析方針をご相談いただけます。
参考文献
- Royston P, Parmar MKB. Flexible parametric proportional-hazards and proportional-odds models for censored survival data, with application to prognostic modelling and estimation of treatment effects. Statistics in Medicine. 2002;21(15):2175-2197.
- Lambert PC, Royston P. Further development of flexible parametric models for survival analysis. The Stata Journal. 2009;9(2):265-290.
- Jackson CH. flexsurv: A platform for parametric survival modeling in R. Journal of Statistical Software. 2016;70(8):1-33.
- Durrleman S, Simon R. Flexible regression models with cubic splines. Statistics in Medicine. 1989;8(5):551-561.
- Moertel CG, Fleming TR, Macdonald JS, et al. Levamisole and fluorouracil for adjuvant therapy of resected colon carcinoma. The New England Journal of Medicine. 1990;322(6):352-358.
