【統計解析】生存時間解析　コックス回帰分析（コックス比例ハザードモデル）

コックス回帰分析（Cox Regression Analysis）のメモ。

【目次】

計算式等
Rでの実行
プログラムコード
参考

コックス回帰分析はハザード関数を目的変数とした回帰分析手法。コックス比例ハザードモデル、比例ハザードモデルなどとも呼ばれる。
任意の共変量（説明変数）を投入することができるセミパラメトリックモデル。

コックス回帰分析の回帰係数を指数変換したものはハザード比（Hazard Ratio、HR）。
ロジスティクス回帰の回帰係数を指数変換したものはオッズ比。

計算式等

■ ハザード関数

$\quad h\left(t\right)=Pr\left(T=t \ | \ T \geqq t \right)=\dfrac{f\left(t\right)}{S\left(t-1\right)}$

（イメージ図、離散値の場合）

* 縦軸が確率密度関数 f (T)、横軸が時間 T
* 死亡関数 F(t)＝時点ｔまでにイベント発生する確率
* 生存関数 S(t)＝時点ｔ以降も生存している（イベント発生しない）確率。
* ハザード関数 H(t)＝時点ｔの直前まで生存している集団（リスク集団、Number At Risk）のうち、時点ｔでイベント発生する確率

ハザード関数、生存関数、死亡関数、確率密度関数・・・などの関係性は別記事にまとめている。

■ コックス回帰分析モデル

$\quad h \left(x,\ t \right) = h_{0} \left( t \right) ^{ \ \beta ^{T} X } = h_{0} \left( t \right) \cdot \exp \{ \beta ^{T} X \}$

$h_{0} \left( t \right)$ は、ベースラインハザード関数（baseline hazard function）と呼ばれ、時点ｔのみに依存するノンパラメトリックな部分（共変量に依存しない潜在因子）。

$\exp { \beta ^{T} X }$ は、相対危険度関数（relative risk function）と呼ばれ、時点ｔに依存しないパラメトリックな部分。

ベースラインハザード関数は推定不要。
例えば、治療群（治療A＝１、治療B＝０）を共変量としたとき、ハザード比（基準：治療B）は次のとおりで、導出の過程でベースラインハザード関数は打ち消される。

$\verb|（治療A）| h_{A} \left( t \right) = h_{0} \left( t \right) \cdot \exp \{ \beta _{1} \times 1 \} = h_{0} \left( t \right) \cdot \exp \{ \beta_{1}\}$

$\verb|（治療B）| h_{B} \left( t \right) = h_{0} \left( t \right) \cdot \exp \{ \beta _{1} \times 0 \} = h_{0} \left( t \right)$

$\quad \dfrac{h_{A} \left( t \right)}{h_{B} \left( t \right)} = \dfrac{h_{0} \left( t \right) \cdot \exp \{ \beta_{1}\}}{h_{0} \left( t \right)}=\exp \{ \beta_{1} \}$

この例ではハザード比は exp{β1} となる。
（治療Aのイベント発生リスクは、治療Bと比べて exp{β1} 倍高まると解釈される。）

なお、コックス回帰分析は「ハザード比が時間に依存せず一定だ」という仮定（比例ハザード性）の上に成り立つ。
つまり、相対危険度関数は時点ｔに依存しないという強い仮定が必要。
上記の例では、治療群の影響は時間に依存せず一定という仮定となり、1年後に治療群の影響が逆転したら解釈が難しくなる。

■ 比例ハザード性の評価

① Schoenfeld残差プロットによる評価

　縦軸：Schoenfeld残差（シェーンフィールド残差）、横軸：時間のプロット。
　平滑化曲線に時間依存的な変化がなければ比例ハザード性を否定しない。
　Rの cox.zph では、Schoenfeld残差を用いた比例ハザード性の検定（Schoenfeld Individual Test）ができる。「帰無仮説：比例ハザード性を満たす」「対立仮説：比例ハザード性を満たさない」の仮説の下、p値が非有意で帰無仮説を棄却できなければ、比例ハザード性を否定しない。

② 二重対数プロット（log-logプロット）による評価

　縦軸：log(-log S(t))、横軸：時間のプロット。
　群間のプロットが平行で時間依存性がなければ、比例ハザード性を否定しない。
　累積ハザード関数と累積生存率の関係は以下。

$\qquad \log \widehat{H}\left( t\right) = \log\left(-\log\ \widehat{S}\left( t \right)\right)$

■ 残差診断
martingale残差（マルチンゲール残差）、Schoenfeld残差（シェーンフィールド残差）などの複数の方法が提案されている。

Rでの実行

以前のKaplan-meierと同じく、MASSライブラリーのgehanデータでコックス回帰を実施してみる。
【統計】生存時間解析カプランマイヤー法（Kaplan-meier method） - こちにぃるの日記

library(MASS)
library(survival)
library(survminer)
ADTTE <- data.frame(
   AVAL   = gehan$time
  ,EVENT  = gehan$cens
  ,CNSR   = ifelse(gehan$cens==0, 1, 0) #-- イベント=0, 打ち切り>0の形式
  ,TRT01A = factor(gehan$treat, levels=c("control","6-MP")) )

#-- Kaplan-meier plot
KM1 <- survfit(formula    = Surv(time=AVAL, event=EVENT) ~ TRT01A
               ,conf.type = "log-log"
               ,conf.int  = .95
               ,type      = "kaplan-meier"
               ,error     = "greenwood"
               ,data      = ADTTE)
 
ggsurvplot(KM1
           ,pval             = TRUE
           ,conf.int         = TRUE
           ,risk.table       = TRUE
           ,risk.table.col   = "strata"
           ,linetype         = "strata"
           ,surv.median.line = "hv"
           ,ggtheme          = theme_bw()
           ,palette          = c("blue", "darkgreen") )

死亡を目的変数、治療群（Control vs. 6-MP）を共変量としたコックス回帰を実行。
Control群を基準としたハザード比は 0.2076 で、信頼区間 (0.09251, 0.4659) は１を挟まず、治療群に有意差が見られる。

> COX1 <- coxph(Surv(time=AVAL, event=EVENT) ~ TRT01A, ADTTE, ties="efron")
> summary(COX1)
Call:
coxph(formula = Surv(time = AVAL, event = EVENT) ~ TRT01A, data = ADTTE, 
    ties = "efron")

  n= 42, number of events= 30 

              coef exp(coef) se(coef)      z Pr(>|z|)    
TRT01A6-MP -1.5721    0.2076   0.4124 -3.812 0.000138 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

           exp(coef) exp(-coef) lower .95 upper .95
TRT01A6-MP    0.2076      4.817   0.09251    0.4659

Concordance= 0.69  (se = 0.041 )
Likelihood ratio test= 16.35  on 1 df,   p=5e-05
Wald test            = 14.53  on 1 df,   p=1e-04
Score (logrank) test = 17.25  on 1 df,   p=3e-05

Schoenfeld残差プロットを作成し、比例ハザード性を評価する。
平滑化曲線に時間依存的な変化は見られず、かつ比例ハザード性の検定（Schoenfeld Individual Test）も非有意であるため、比例ハザード性は否定されない。

cox.zph1 <- cox.zph(COX1)
ggcoxzph(cox.zph1)

二重対数プロット（log-logプロット）でも比例ハザード性を評価してみる。
治療群間のプロットは平行で時間依存性は見られないため、本プロットでも比例ハザード性は否定されない。

ggsurvplot(KM1, fun="cloglog")

マルチンゲール残差とその標準化残差を用いて残差プロットを作成する。
０を中心にばらつき、明らかな偏りは見られない。（マルチンゲール残差は上限：１、下限：－∞）

scatter.smooth(residuals(COX1, type="martingale"))
scatter.smooth(residuals(COX1, type="deviance"))

プログラムコード

library(MASS)
library(survival)
library(survminer)
 
#-- Cox回帰
COX1 <- coxph(Surv(time=AVAL, event=EVENT) ~ covariate, ADTTE, ties="efron")
summary(COX1)

##-- Schoenfeld残差プロット＆比例ハザード性の検定（Schoenfeld Individual Test）
cox.zph1 <- cox.zph(COX1)
ggcoxzph(cox.zph1)

##-- 二重対数プロット（log-logプロット）
KM1 <- survfit(Surv(time=AVAL, event=EVENT) ~ TRT01A
               ,conf.type = "log-log"
               ,conf.int  = .95
               ,type      = "kaplan-meier"
               ,error     = "greenwood"
               ,data      = ADTTE)
ggsurvplot(KM1, fun="cloglog")

##-- 残差診断
scatter.smooth(residuals(COX1, type="martingale"))

タイの計算は "efron","breslow","exact" の３つ。デフォルトはefron。
残差診断は type で "martingale","deviance","schoenfeld" など指定。

proc phreg data=ADTTE plots(overlay)=(cumhaz survival) ;
    class covariate(ref="XXX") / param=glm ;
    model AVAL * CNSR(1) = covariate / rl ties=exact;
    hazardratio covariate / cl=wald ;
    baseline out=COX1 survival=all cumhaz=all ;
    output   out=COX2 atrisk=AtRisk survival=SURVIVAL
             resmart=RESMART resdev=RESDEV ressch=all ;
    id SUBJID ;
run;