比例ハザード性の仮定が成立しているか検証する[R]

パッケージとデータの準備
補対数-対数プロットを使って検証する
Schoenfeld残差を使って検証する
- ggcoxdiagnostics( )を使って残差プロットを描く
時間依存性共変量を使って検証する
おわりに
参考資料

Cox比例ハザードモデル（Cox proportional hazard model, CPH model）では、
ハザード関数 $h(X,t)$ が時間依存性・共変量非依存性の部分 $h_0(t)$ と共変量依存性・時間非依存性の部分 $e ^{\Sigma \beta X}$ の積で構成される。

ベースラインハザード $h_0(t)$ が時間で変化しても共変量の層間での比は一定であるという「比例ハザード性」が満たされていることを仮定したモデルなので、実際のデータでもこの仮定が成り立っているか確認する必要がある。

検証には、

補対数-対数プロット
残差
時間依存性共変量

を用いた方法がある。

パッケージとデータの準備

library(survival)
library(survminer)
library(ggplot2)
library(tidyverse)
d <- veteran %>% 
  mutate(age_cat = if_else((age<60), "under60", "over60"),
         age_cat = factor(age_cat),
         trt = factor(trt, labels=c("standard", "test")),
         prior = factor(prior, labels=c("N"," Y")))

> head(d,3)
       trt celltype time status karno diagtime age prior age_cat
1 standard squamous   72      1    60        7  69     N  over60
2 standard squamous  411      1    70        5  64     Y  over60
3 standard squamous  228      1    60        3  38     N under60

補対数-対数プロットを使って検証する

生存関数を $S(t)$ , ハザード関数を $h(t)$ , 累積ハザード関数 $H(t)$ とすると、

$\displaystyle \begin{align} S(t) &= e^{-H(t)} \\ H(t) &= \int_0^t h(u) du \end{align}$

の関係が成り立つ。
生存関数の補対数対数（complementary log-log）を取ると、

$\displaystyle \begin{align} \log (-\log S(t)) &= \log(-\log e^{-H(t)}) \\ &= \log H(t) \\ &= \log \int_0^t h(u) du \\ \end{align}$

となる。

比例ハザードモデルでは $\displaystyle h(t) = h_0(t)e^{\Sigma \beta X}$ が成立しちている（と仮定している）ので、上記式はさらに次のように変形できる（ $S_0(t)$ はベースライン生存関数）。

$\displaystyle \begin{align} \log (-\log S(t)) &= \log \int_0^t h_0(u) e^{\Sigma \beta X} du \\ &= \log e^{\Sigma \beta X} \int_0^t h_0(u) du \\ &= \log e^{\Sigma \beta X} + \log \int_0^t h_0(u) du \\ &= \Sigma \beta X + \log (-\log S_0(t)) \end{align}$

時間tまたはlog(t)に対して $\log (-\log S(t))$ （= 累積ハザードの対数）をプロットすると、共変量で決まる分 $\Sigma \beta X$ だけ切片が異なる平行な線が描かれる。
逆に、層間で補対数対数プロットが並行でないときは比例ハザード性の仮定が成立していない。

例えば治療（trt）の層別に補対数対数プロットは、survfitオブジェクトをplot()に渡すときに引数fun="cloglog"を指定すれば描ける。

plot(survfit(Surv(time, status)~trt, data=d), fun="cloglog",
     col=c("red", "blue"), 
     lty=c("solid", "dotted"),
     lwd=2,
     xlab="log(Time)", ylab="log(-log S(t))")

Colletのテキストによると、横軸が時間tではなく対数時間になっている理由は

plotting the log-cumulative hazard functions against the logarithm of t, rather than t itself, is a useful diagnostic in parametric modeling and so this form of plot is generally used.

とのこと。生存関数が指数関数だったら補対数対数プロットが直線になる、といったことだろうか。
平行かどうかは主観的に判断する。明確な違反がなさそうなら比例ハザード性を仮定してもよい（図1はOK）。

層別変数が連続変数の場合はカテゴリー化が必要。
また共変量が複数ある場合は、それらの層を組み合わせたカテゴリーで補対数対数プロットを描くことになるが、各カテゴリーに含まれる数が少なくなると信頼性が低くなるので限界あり。

Schoenfeld残差を使って検証する

共変量が複数ある場合、補対数対数プロットでは個々の変数について比例ハザード性が成立しているかどうかは検証できない。そんなときはSchoenfeld（シェーンフェルド）残差を使って検証する。

Schoenfeld残差はイベントを発生した人の各共変量に対して定義される。「イベントを発生した人の共変量と、イベントを発生した時点におけるat risk集団の共変量を各人のハザードで重み付けした平均との差」らしいが、特徴だけ押さえておく。

特徴：

イベントを発生した人にだけ定義される
各共変量毎に定義される
総和=0

もし、ある変数についてSchoenfeld残差は時間と関連が認められなければ、その変数について比例ハザード性の仮定が成り立っていると考えられる。

cox.zph( )を使って検定をして、ggcoxzph( )でSchoenfeld残差プロットを作成できる。

> fit <- coxph(Surv(time, status)~trt + celltype + age_cat + prior, data=d)
> fit_test <- cox.zph(fit)
> fit_test
          chisq df     p
trt       0.906  1 0.341
celltype  9.109  3 0.028
age_cat   2.942  1 0.086
prior     4.000  1 0.046
GLOBAL   15.657  6 0.016

ggcoxzph(fit_test)

celltypeとpriorについてはP値が小さく、「比例ハザード性が成立する」をいう帰無仮説が棄却されているので比例ハザード性が成立しているかどうか疑わしい。

ggcoxdiagnostics( )を使って残差プロットを描く

typeで残差の種類、ox.scaleで横軸に取る変数を指定する。

typeとして指定できるもの（の一部）：

"martingale"：Martingale（マーチンゲール）残差。Cox-Snell（コックス-スネル）残差*1を平均が0になるように修正したもの。時間や線形予測子に対して特定の傾向がなければよい。
"deviance"：Deviance（デビアンス）残差。Martingale残差がもっと0に関して対称に分布するように修正したもの。時間や線形予測子に対して特定の傾向がなければよい。
"score"：Schoenfeld残差の修正版のようなもの。
"dfbeta"：ある観測値を抜いたときにβがどれくらい変化するか（Δβ）
"defbetas"：ΔβをSE(β)で割って標準化したもの。
"schoenfeld"：前述
"scaledsch"：係数βの共分散行列で重み付けしたSchoenfeld残差。オリジナルよりも鋭敏らしい。

ox.scaleとして指定できるもの：

"linear.prediction"：線形予測子 $\Sigma \beta X$ を横軸にとる。 $\Sigma \beta X$ はイベントの発生しやすさを表している（はず）なのでリスクスコアと呼ばれることもある。
"observation.id"：イベントを発生した順に並べたもの。
"time"：イベントを発生した時間。

指定の仕方は下のコードを参照。

ggcoxdiagnostics(fit, type="deviance", ox.scale="observation.id")

ggcoxdiagnostics(fit, type="dfbeta", ox.scale="time")

ggcoxdiagnostics( )を使ってSchoenfeld残差を時間に対してプロットしようとすると何故かエラーが出て上手くいかないので、ggcoxzph( )を使った方がよさそう。

時間依存性共変量を使って検証する

ある変数のハザードに対する影響に時間依存性があるかどうかを、変数と時間の交互作用項を使って検証する方法。
例えば、

$\displaystyle h(t) = h_0(t) e^{\beta_0 + \beta_1 X + \beta_2 Xt}$

のようにXとtの交互作用項を加え、その係数が $\beta_2$ が0に近ければ、変数Xのハザードに対する影響は時間非依存性と考えられる。

おわりに

比例ハザード性は重要な仮定なので必ずチェックする。生存曲線が交差するときは要注意。
キャットタワーを設置しましたが、なかなか登ってくれません。

※ 生存時間でよく登場する関数については別の記事でまとめ直しました（2022-07-17 追加）

necostat.hatenablog.jp

ランキング参加中

統計

参考資料

生存時間解析について自己学習できる詳細な解説付きテキスト

Survival Analysis: A Self-Learning Text, Third Edition (Statistics for Biology and Health) (English Edition)

作者:Kleinbaum, David G.
Springer

Amazon

数式に関する説明が詳細

Modelling Survival Data in Medical Research (Chapman & Hall/CRC Texts in Statistical Science Book 115) (English Edition)

作者:Collett, David
Chapman and Hall/CRC