費用対効果分析の担当者から、こう頼まれたとします。「試験データから生涯の平均生存期間を出してください」。手元にあるのは追跡9年の第III相試験のデータで、その時点でまだ半数以上が生存しています。Kaplan-Meier曲線は9年で途切れていて、その先は描けません。
ここで分布を仮定します。指数分布でも、ワイブル分布でも、対数正規分布でも、追跡期間内の当てはまりはほとんど変わりません。ところが30年まで伸ばすと、平均生存期間が11.1年から13.7年まで開きます。
パラメトリックモデルの使いどころは、この外挿にあります。同時に、いちばん危ないところでもあります。ここでは分布を仮定すると何ができるようになるのかを確認したうえで、選択をどう正当化するか、効果指標としてハザード比の代わりに何が使えるかを見ていきます。
分布を仮定すると何が決まるか
生存時間$T$に分布を仮定するというのは、$S(t)$の形をパラメータ数個で書き切るということです。いちばん単純な指数分布では、ハザードが時間によらず一定と置きます。
$$h(t) = \lambda, \qquad S(t) = \exp(-\lambda t) \tag{1}$$
式(1)は、いつの時点でも危険度が同じという仮定です。この仮定からは記憶なし性が出てきます。5年生きた患者のその後の生存分布が、診断直後の患者と同じになる。癌の術後経過にも術後死亡にも当てはまらないので、医学研究で指数分布を単独で使う場面はほとんどありません。
ワイブル分布は形状パラメータ$\gamma$を加えて、ハザードが単調に増えるか減るかを表せるようにしたものです。
$$h(t) = \frac{\gamma}{\lambda}\left(\frac{t}{\lambda}\right)^{\gamma-1} \tag{2}$$
式(2)は$\gamma > 1$なら時間とともに危険が増し、$\gamma < 1$なら減っていく形です。$\gamma = 1$のとき指数分布に戻ります。術後早期の危険が高くその後下がる経過なら$\gamma < 1$、加齢や再発の蓄積で危険が増すなら$\gamma > 1$になります。
対数正規分布と対数ロジスティック分布は、$\log T$に正規分布またはロジスティック分布を置きます。この2つは山型のハザードを表せるので、周術期のように危険がいったん上がってから下がる経過に使えます。ただし裾が厚く、長期の生存率を高めに見積もる方向に働きます。
共通しているのは、追跡が終わったあとの形まで仮定が決めてしまうことです。データが何も言っていない領域の生存率が、分布族の選択だけで決まります。
追跡期間の中では見分けがつかない
survivalパッケージのcolonは、大腸癌の術後補助化学療法の試験データです。Lev+5FU群304例、追跡は最長9年、死亡123件。この群に4つの分布を当てはめてみます。
library(survival)
adjuvant <- subset(colon, etype == 2 & rx == "Lev+5FU")
dists <- c("exponential", "weibull", "lognormal", "loglogistic")
fits <- lapply(dists, function(d)
survreg(Surv(time, status) ~ 1, data = adjuvant, dist = d))
sapply(fits, AIC)#> exponential weibull lognormal loglogistic
#> 2314.3 2315.4 2306.5 2310.1
AICは2306から2315の範囲に収まっています。最良の対数正規分布と最悪のワイブル分布で9の差ですが、指数分布と対数ロジスティック分布は4しか違いません。判断の根拠として強いとは言えない差です。
外挿すると答えが分かれる
同じ4つの当てはめを30年まで伸ばします。
費用対効果分析では、この平均生存期間がそのまま質調整生存年の計算に入ります。2.5年の違いは増分費用効果比を大きく動かし、償還の可否が変わることもあります。追跡期間内でAICが4しか違わない2つのモデルが、そこまで違う結論を出します。
だから、外挿を含む解析でAICだけを根拠に分布を選ぶことは勧めません。AICが測っているのは観測された範囲での当てはまりであって、その外側の形については何も言っていないからです。
では、どう決めるか。私なら、まず疾患の自然史から考えます。術後5年で再発しなければその後の死亡ハザードは一般人口に近づくのか、それとも遅発再発が続くのか。この知識は臨床の側にあるので、統計家だけでは決められません。共著者と話して、20年後の生存率がどのくらいなら医学的にあり得るかの範囲を先に決めておく。そのうえで、その範囲に収まる分布のなかから選びます。
選び切れないときは、複数の外挿を並べて報告します。英国のNICE Decision Support Unitの技術文書14は、この形を標準的な手順として示しています。「対数正規分布を採用したが、ワイブル分布ならICERはこうなる」と書いてあれば、読者は不確実性の幅を自分で判断できます。単一の分布を選んで根拠を書かない報告のほうが、査読では通りにくくなります。
当てはまりをどう診断するか
AICは相対的な指標なので、比べた分布がすべて外れていても、そのうちいちばんましなものを選んでしまいます。絶対的な当てはまりは別に見る必要があります。
ワイブル分布なら、式(2)から$\log(-\log S(t))$が$\log t$の1次式になります。
$$\log(-\log S(t)) = \gamma \log t – \gamma \log \lambda \tag{3}$$
式(3)は、Kaplan-Meier推定量をこの座標に置いたとき、ワイブル分布が正しければ点が直線に並ぶことを表しています。傾きが形状パラメータ$\gamma$です。
もう1つ、分布によらず使える診断としてCox-Snell残差があります。当てはめたモデルから$r_i = -\log \hat{S}(t_i)$を計算すると、モデルが正しければこの$r_i$は率1の指数分布に従います。$r_i$を生存時間だと思ってKaplan-Meier法にかけ、累積ハザードを描いたときに傾き1の直線になるかを見ます。
weib <- survreg(Surv(time, status) ~ 1, data = adjuvant, dist = "weibull")
resid <- -log(pweibull(adjuvant$time, shape = 1 / weib$scale,
scale = exp(coef(weib)), lower.tail = FALSE))
check <- survfit(Surv(resid, adjuvant$status) ~ 1)
plot(check$time, -log(check$surv), type = "s",
xlab = "Cox-Snell residual", ylab = "Cumulative hazard")
abline(0, 1, lty = 2)
図1では4本とも合っているように見えたのに、こちらの座標では外れが見えます。生存率の目盛りで見ると小さな差が、対数変換すると拡大されるためです。当てはまりの診断は、生存曲線の重ね書きだけで済ませないほうがよいでしょう。
ただし、この程度の外れで分布を却下すべきかは別の判断になります。図1の重ね書きが臨床的に十分な精度なら、追跡期間内の記述としては使えます。問題になるのは、その形のまま20年先まで伸ばしたときです。
効果を時間の比で表す
パラメトリックモデルにはもう1つの使い道があります。survregが当てはめているのは加速故障時間モデルで、共変量が生存時間そのものを何倍にするかという形をしています。
$$\log T = \beta_0 + \beta_1 x + \sigma W \tag{4}$$
式(4)は、共変量$x$が1増えると生存時間が$\exp(\beta_1)$倍になる、と読めます。この$\exp(\beta_1)$を時間比と呼びます。
aft <- survreg(Surv(time, status) ~ sex, data = lung, dist = "weibull")
exp(coef(aft)["sex"]) # 時間比
cox <- coxph(Surv(time, status) ~ sex, data = lung)
exp(coef(cox)) # ハザード比#> 時間比 1.485
#> ハザード比 0.588同じデータから、女性の生存時間は男性の1.49倍、女性の死亡ハザードは男性の0.588倍という2つの表現が出てきます。「生存期間が1.5倍」のほうが、共著者や患者への説明には通りやすい言い方でしょう。
ワイブル分布に限っては、比例ハザードと加速故障時間の両方の性質が成り立つので、この2つは$\mathrm{HR} = (\text{時間比})^{-\gamma}$で結ばれます。ほかの分布では比例ハザードが成り立たないため、時間比だけが定義されます。
ここが実務では効いてきます。比例ハザード性が崩れている状況では、Cox回帰の単一のハザード比が解釈しにくくなります。対数正規分布や対数ロジスティック分布の加速故障時間モデルなら、比例ハザードを仮定せずに1つの効果指標を出せます。分布の仮定を新たに負う代わりに、比例ハザードの仮定を外す形です。どちらの仮定のほうが自分のデータで納得できるかで選ぶことになります。
どちらの仮定を選ぶか
追跡期間内の記述と群間比較が目的なら、Kaplan-Meier法とCox回帰で足ります。分布を仮定する必要があるのは、追跡の先を推定するとき、比例ハザードが崩れていて別の効果指標が要るとき、そしてイベント数が少なく推定効率を上げたいときです。
比例ハザード性そのものの確認は比例ハザード仮定の検証、共変量調整の標準的な枠組みはCox比例ハザードモデルで扱います。累積ハザードの側から当てはまりを見る話はNelson-Aalen推定量、そもそも観測範囲の外を推定できない理由はKaplan-Meier推定量にあります。
費用対効果分析のために生涯の生存期間を出す、あるいは比例ハザードが崩れた試験で効果をどう報告するか。分布の選択とその正当化は、疾患の自然史と報告先の要求の両方を見ないと決まりません。Dr.データサイエンスでは、こうした解析方針のご相談を承っています。
参考文献
Collett D. Modelling Survival Data in Medical Research. 3rd ed. Chapman and Hall/CRC; 2015.
Latimer NR. Survival analysis for economic evaluations alongside clinical trials: extrapolation with patient-level data. Medical Decision Making. 2013;33(6):743-754.
Latimer NR. NICE DSU Technical Support Document 14: Survival analysis for economic evaluations alongside clinical trials. National Institute for Health and Care Excellence Decision Support Unit; 2011.
Royston P, Parmar MKB. Flexible parametric proportional-hazards and proportional-odds models for censored survival data. Statistics in Medicine. 2002;21(15):2175-2197.
Klein JP, Moeschberger ML. Survival Analysis: Techniques for Censored and Truncated Data. 2nd ed. Springer; 2003.
