免疫チェックポイント阻害薬の臨床試験で、こんな生存曲線を見たことはないでしょうか。投与開始から数か月は2群の曲線がほとんど重なっていて、その後になって治療群が離れていく。あるいは逆に、外科的治療の比較で、術後しばらくは手術群の死亡が多く、時間が経つと逆転する。
曲線が交差しているわけではない。それでも、治療効果が追跡期間を通じて同じとは言いにくい。
ところが、この2群にCox回帰をあてはめてハザード比を1つ報告すると、その値は「前半と後半を混ぜた何か」になります。その値が何をどう反映しているのかは、イベントがいつ起きたかによって変わります。読者にはそれが見えません。
比例ハザード性の検証は、この事態が起きていないかを確かめる手続きです。以下では、Schoenfeld残差が何を見ているのか、検定が有意だったときに何をすればよいのか、そして有意でなかったときに何を主張してよいのかを順に見ていきます。
Coxモデルが置いている仮定
Cox回帰は、個体$i$のハザードを次の形に置きます。
$$ h_i(t) = h_0(t)\exp(\mathbf{x}_i^\top\boldsymbol{\beta}) \tag{1} $$
$h_0(t)$は基準ハザード、$\mathbf{x}_i$は共変量ベクトル、$\boldsymbol{\beta}$は回帰係数です。
式(1)で注目したいのは、右辺のうち時間$t$を含むのが$h_0(t)$だけだという点です。共変量の効果を表す$\exp(\mathbf{x}_i^\top\boldsymbol{\beta})$には$t$が入っていません。つまり効果の大きさは、時間によらず同じだと置かれています。
この構造から、2個体のハザード比は次のようになります。
$$ \frac{h_i(t)}{h_j(t)} = \exp\left((\mathbf{x}_i – \mathbf{x}_j)^\top\boldsymbol{\beta}\right) \tag{2} $$
式(2)には$t$が現れません。診断直後だろうと5年後だろうと、2群のハザードの比は同じ値をとる。これが比例ハザード性です。冒頭に挙げた免疫療法の例は、まさにこの仮定が成り立っていない状況にあたります。
仮定が破れている状況は、係数が時間の関数になっていると書けます。
$$ h(t\mid x) = h_0(t)\exp(\beta(t)\, x) \tag{3} $$
式(3)で$\beta(t)$が定数なら式(1)に戻ります。ですから確かめたいのは、$\beta(t)$が定数からどれだけ離れているか、ということになります。
Schoenfeld残差が見ているもの
イベントが起きた時刻$t_k$で、そのイベントを起こした個体の共変量と、その時点でまだリスク下にある集団の共変量の重み付き平均との差を取ります。これがSchoenfeld残差です。
$$ \hat{r}_k = \mathbf{x}_{(k)} – \hat{E}[\mathbf{X}\mid\mathcal{R}_k] \tag{4} $$
式(4)は、部分尤度のスコア関数をイベント時刻ごとに分解したものにあたります。
直感的に言えば、こういうことです。モデルが「共変量の値が高い個体ほどイベントが起きやすい」と言っているなら、実際にイベントを起こした個体の共変量は、その時点のリスク集合の平均より高いはずです。残差はその「予想からのずれ」を各時点で記録しています。効果の大きさが時間によらず一定なら、このずれは時間方向に系統的な傾向を持たないはずです。
GrambschとTherneauは1994年に、この残差を分散で尺度調整すると、その期待値が時間依存係数そのものに近似することを示しました。
$$ E[\widetilde{r}_k] \approx \beta(t_k) – \bar{\beta} \tag{5} $$
式(5)が言っているのは、尺度調整済み残差を時間に対してプロットすれば、その散布図が$\beta(t)$の形をおおよそ描く、ということです。傾きがあれば、係数が時間とともに動いています。検定はこの傾きに対する回帰として構成され、帰無仮説は「残差と時間に関連がない」、すなわち$\beta(t)$が定数であることになります。Rのcox.zphがこれを実装しています。
ここは少し注意が必要です。
「残差プロットに傾きがなければ比例ハザード性は大丈夫」と考えたくなりますが、そこまで単純ではありません。この方法が見ているのは、あくまで残差と時間の間に系統的な傾きがあるかどうかです。たとえば、効果が途中まで強くなり、その後で弱くなって元に戻るような変化があったとしましょう。前半の上向きと後半の下向きが打ち消し合って、全体としては傾きがほとんど見えなくなります。
ですから私は、検定の$p$値だけで判断せず、必ず残差プロットも見ることを勧めます。数字は1つしか返ってきませんが、絵には形が出ます。
veteranデータで確かめる
survivalパッケージのveteranデータを使います。肺癌患者137名、イベント128件の退役軍人向け臨床試験のデータです。Karnofsky performance status（karno）、治療群（trt）、年齢を共変量に入れます。
library(survival)
fit <- coxph(Surv(time, status) ~ trt + karno + age, data = veteran)
cox.zph(fit)
#> chisq df p
#> trt 0.284 1 0.59439
#> karno 12.003 1 0.00053
#> age 2.100 1 0.14735
#> GLOBAL 19.102 3 0.00026
karnoが$p=0.00053$で、比例ハザード性からの逸脱を示しています。GLOBALも有意です。
ここで、いったん立ち止まります。
「karnoのp値が0.00053だから比例ハザード性は破れている、したがってCox回帰は使えない」と考えたくなりますが、そこまで単純ではありません。問題は、どの程度破れているのかです。統計的に有意な逸脱があることと、その逸脱が研究結果の解釈を大きく変えることは、同じではありません。
trtとageについては、今回のデータからは明確な逸脱は確認されませんでした。ただ、ここで「この2つは問題ない」と言い切るのは避けた方がよいでしょう。理由は後で述べます。
このときkarnoの推定値はこうなっています。karnoの行だけを抜き出します。
#> coef exp(coef) se(coef) z Pr(>|z|)
#> karno -0.0344 0.9661 0.00523 -6.58 4.62e-11
1点あたりのハザード比が0.966、10点あたりに換算すると0.71です。$p$値は$10^{-11}$の桁ですから、この数字だけを見れば、きわめて頑健な結果に見えます。
しかし先ほどの検定は、この0.71が追跡期間を通じた一定の効果ではないと告げています。では、この値は何を表しているのでしょうか。
有意だったときに私がすること
私なら、ここですぐにモデルを組み替えることはしません。まず残差プロットを見て、どのあたりの時間帯でずれが生じているのかを確認します。検定が有意になったという事実だけでは、どういう形の非比例性が起きているのかが分からないからです。
plot(cox.zph(fit)[2])
そのうえで、期間を区切って係数がどう動いているかを見ます。プロットで見当をつけた区切りで、実際に推定してみるということです。
vet2 <- survSplit(Surv(time, status) ~ ., data = veteran,
cut = 90, episode = "period")
fit2 <- coxph(Surv(tstart, time, status) ~ karno:strata(period) + trt + age,
data = vet2)
summary(fit2)$coefficients
karnoの行だけを抜き出します。
#> coef exp(coef) se(coef) z Pr(>|z|)
#> karno:strata(period)period=1 -0.049837 0.9514 0.00641 -7.78 7.2e-15
#> karno:strata(period)period=2 0.000651 1.0007 0.00978 0.07 0.947
90日で区切ると、karnoの効果は前半で10点あたりハザード比0.61、後半では1.01でした。後半は$p=0.95$で、効果がほぼ消えています。
ここで、全期間の0.71が前半と後半の平均だ、と考えたくなります。ただ、そう単純ではありません。
全期間を1つの係数で表した0.71は、前半の効果と後半の効果を別々に示す値ではありません。どの時点でどれだけイベントが起き、そのときどのようなリスク集合が構成されていたかを反映した、全追跡期間を通じた要約値です。このデータでは90日を境に効果が大きく違っていましたから、0.71がどの期間の情報をどの程度反映しているかは、イベントの発生状況にも左右されます。
厄介なのは、この値が母集団のどの部分を指すのかが定まらない点です。追跡期間の設計が違う別の研究と、同じ意味の量として比べることもできません。Hernánが2010年に指摘したのは、この種の値が持つ解釈上の不安定さでした。
では、比例ハザード性が崩れていたらどうするのか。
まず考えられるのは、問題になっている変数を層別化する方法です。層別Coxでその共変量を層に移せば、比例ハザード性の仮定から外せます。次に、効果そのものが時間とともに変化することを、時変係数としてモデルに組み込む方法があります。そして、そもそもハザード比を主要な指標として使わず、制限付き平均生存時間（RMST）や特定時点の生存率差に切り替えるという選択もあります。
私は、最後の方法が適しているケースはかなり多いと考えています。
層別も時変係数も、非比例性を技術的に処理する方法としては正しいのですが、最終的に読者へ渡す数値は、依然としてハザード比か、時間ごとに変わる係数の束です。前者は解釈が定まらず、後者は臨床的な意味を読み取りにくい。これに対してRMSTは「追跡24か月のうち、治療群は平均して何か月長く生存したか」という形で答えが出るので、比例ハザード性を前提としませんし、臨床的な意味も伝わりやすい。Unoらが2014年に提案し、RoystonとParmarも同様の主張をしています。
ただし、すべての研究でRMSTに置き換えればよいわけではありません。
主要評価項目と解析方法を事前に規定している介入研究では、結果を見てから解析計画を変えることはできません。その場合は、事前規定どおりハザード比を主要指標として報告したうえで、比例ハザード性の検定結果と期間別の推定値を感度分析として併記する。これが現実的な落としどころだと思います。
逆に避けたいのは、検定結果に触れずに単一のハザード比だけを報告することです。これは検定を行わなかった場合より悪いと私は考えています。読者は、報告されたハザード比が追跡期間を通じて一定だと受け取るからです。書き手はそうでないと知っている。知っていて書かないのは、検定していないのとは意味が違います。
有意でなかったときに言えること
今度は逆の場合を見ておきます。lungデータで同じ検定をしてみます。
fit3 <- coxph(Surv(time, status) ~ sex + age, data = lung)
cox.zph(fit3)
#> chisq df p
#> sex 2.608 1 0.11
#> age 0.209 1 0.65
#> GLOBAL 2.771 2 0.25
plot(cox.zph(fit3)[1])
いずれも有意ではありません。ここで「比例ハザード性が成り立つことを確認した」と書きたくなります。実際、そう書かれた原稿をしばしば見かけます。
しかし、これは正確ではありません。
検定が有意でないことは、逸脱がないことの証明ではないからです。この検定の検出力はイベント数に依存します。イベントが数十件しかない研究では、実際にはかなりの非比例性があっても、それを検出するだけの力がありません。有意でなかったという結果は、この標本では逸脱を検出できなかった、という事実を述べているに過ぎません。
先ほどveteranのtrtとageについて「問題ない」と言い切るのを避けたのは、これが理由です。同じ論理が、有意でなかったすべての共変量にあてはまります。
StensrudとHernánは2020年に、そもそも検定の有意性を判断基準にすること自体を問題にしています。標本が大きければ、臨床的には無視してよい程度の逸脱でも有意になる。小さければ、重大な逸脱でも見逃される。検定の$p$値は、逸脱の大きさを表す量ではないからです。
では、どう判断すればよいでしょうか。
私は検定結果、残差プロットの形、そして期間によって推定値がどれくらい変わっているかを揃えて見ます。
その中でも、研究者として最終的に知りたいのは、その非比例性が研究結果の解釈を変えるほど大きいのか、ということです。$p$値そのものではありません。前半0.61と後半1.01ほど離れていれば、単一のハザード比だけで結果を説明するのは難しくなります。逆に、前半と後半でほとんど変わらないのであれば、検定が有意でも実際の解析への影響は小さいでしょう。
論文にどう書くか
観察研究であればSTROBEが、統計手法とその前提の検証方法を記載するよう求めています。比例ハザード性の検証は、この前提の検証にあたります。介入研究であればCONSORTが、事前規定された解析と事後解析を区別して報告することを要求します。
Methodsには、検証の手法と判断基準を書きます。「Schoenfeld残差にもとづくGrambsch-Therneau検定により比例ハザード性を評価し、あわせて尺度調整済み残差の時間に対するプロットを確認した」といった記述です。判断基準を$p<0.05$と書くのであれば、その閾値を超えたときに何をするかまで、あらかじめ決めておく必要があります。書いていない基準を結果が出てから持ち出すと、事後解析ではないかと疑われます。
Resultsには、検定統計量と$p$値を書きます。逸脱があった場合は、それをどう扱ったかまで書きます。期間別の推定値を出したなら、その結果も載せます。
査読でよく問われるのは、逸脱を認めながら対処を書いていない場合と、対処はしたけれども事前規定だったのか事後の判断だったのかが読み取れない場合です。この2点は先回りして明記しておくと、やり取りが1往復減ります。
Limitationsには、検定の検出力について触れておきます。有意でなかった場合に「比例ハザード性が成立することを確認した」と書くのは書き過ぎです。「本標本では比例ハザード性からの明確な逸脱は検出されなかった」が、言える範囲に合った書き方になります。
実務での判断
比例ハザード性の検定は、モデルの合否を判定する装置ではありません。これから報告しようとしている数値が何を指しているのかを、書き手が確かめるための手続きです。
有意だったなら、ハザード比は追跡期間の構成に依存する量に変わっています。期間を区切って、どれだけ動いているかを見てください。動きが大きければ、単一の値で報告することの妥当性を疑うべきです。有意でなかったなら、逸脱が検出されなかったと書き、成立を確認したとは書かない。
どちらの場合も、検定結果を論文に書かないという選択はありません。
ここまで読んで、では自分の研究ではどう判断すればよいのか、と思われた方もいるかもしれません。実際のデータでは、検定結果だけでなく、残差の形、イベントの発生時期、期間別の推定値、そして事前に定めた解析計画まで含めて判断することになります。Dr.データサイエンスでは、こうした解析方針のご相談を承っています。
参考文献
- Schoenfeld D. Partial residuals for the proportional hazards regression model. Biometrika. 1982;69(1):239-241.
- Grambsch PM, Therneau TM. Proportional hazards tests and diagnostics based on weighted residuals. Biometrika. 1994;81(3):515-526.
- Therneau TM, Grambsch PM. Modeling Survival Data: Extending the Cox Model. Springer; 2000.
- Hernán MA. The hazards of hazard ratios. Epidemiology. 2010;21(1):13-15.
- Uno H, Claggett B, Tian L, et al. Moving beyond the hazard ratio in quantifying the between-group difference in survival analysis. J Clin Oncol. 2014. doi:10.1200/JCO.2014.55.2208
- Stensrud MJ, Hernán MA. Why test for proportional hazards? JAMA. 2020;323(14):1401-1402. doi:10.1001/jama.2020.1267
前の記事：部分尤度：Coxモデルの推定理論
次の記事：Coxモデルの変数選択と診断：残差・影響点分析
