Cox比例ハザードモデル

R

Cox比例ハザードモデルについて勉強したくなったのでまとめ。パラメーター推計に関心があるので、そこを詳しくやる一方、他の点については省略。パラメーター推計以外にも興味がある場合は参考LINKを参照してください。

Cox比例ハザードモデルはハザードレート（単位時間当たりのデフォルト・死亡率的な量）をモデル化するために使うもんで、ハザードレートh(t)とすると、Cox比例ハザードモデルは
h(t|z_1,\cdots,z_n) = h_0(t) \exp(z_1 \omega_1 + \cdots + z_n \omega_n)
と記述できる。ここでz_i,i \in (1,\cdots,n)は共変量と呼ばれるハザードに関係した変数（医学だと性別・投薬の有無等等）で、\omega_iはその重み・回帰パラメーター、h_0(t)はベースラインハザードと呼ばれる（共変量に依存しない基準線）。やりたいことは実際のハザードデータから\omega_i,i \in (1,\cdots,n)を推計すること。


まずはマニュアルにあるサンプルの写経＆解釈から。サンプルデータを以下のように作成する。

> test1 <- list(time=c(4,3,1,1,2,2,3), 
              status=c(1,1,1,0,1,1,0), 
              x=c(0,2,1,1,1,0,0), 
              sex=c(0,0,0,0,1,1,1)) 
> test1
$time
[1] 4 3 1 1 2 2 3

$status
[1] 1 1 1 0 1 1 0

$x
[1] 0 2 1 1 1 0 0

$sex
[1] 0 0 0 0 1 1 1

各項目の意味合いは

  • time：時点
  • status：生きてるか死んでるか
  • x:ある共変量
  • sex:性別


このデータからある時点での生死を表すためデータオブジェクト「survival object」を作るためにはSurv関数を使う。第一引数に”時間”をぶち込んで、第二引数に生きてるのか死んでるのかを突っ込む。通常は「0：生きてる」・「1：死亡」として扱うらしい（マニュアルより）

> Surv(test1$time, test1$status)

[1] 4  3  1  1+ 2  2  3+

パラメータを推定するには関数coxphを以下のように使う。

> test1.coxph1 <- coxph(Surv(time, status) ~ x + sex, test1) 
> test1.coxph2 <- coxph(Surv(time, status) ~ x + strata(sex), test1)

test1.coxph1・test1.coxph2の２つの結果を算出した。ここではxとsexが共変量になっているに見えるが、test1.coxph2の方ではxのみが共変量となっている。実際に、summary関数で結果を見てみると

> summary(test1.coxph1)
Call:
coxph(formula = Surv(time, status) ~ x + sex, data = test1)

  n= 7, number of events= 5 

      coef exp(coef) se(coef)     z Pr(>|z|)
x   0.7812    2.1841   0.7976 0.979    0.327
sex 0.9338    2.5441   1.4081 0.663    0.507

    exp(coef) exp(-coef) lower .95 upper .95
x       2.184     0.4579    0.4575     10.43
sex     2.544     0.3931    0.1610     40.19

Rsquare= 0.15   (max possible= 0.822 )
Likelihood ratio test= 1.13  on 2 df,   p=0.567
Wald test            = 0.96  on 2 df,   p=0.619
Score (logrank) test = 1.05  on 2 df,   p=0.5927

> summary(test1.coxph2)
Call:
coxph(formula = Surv(time, status) ~ x + strata(sex), data = test1)

  n= 7, number of events= 5 

    coef exp(coef) se(coef)     z Pr(>|z|)
x 0.8023    2.2307   0.8224 0.976    0.329

  exp(coef) exp(-coef) lower .95 upper .95
x     2.231     0.4483    0.4451     11.18

Rsquare= 0.144   (max possible= 0.669 )
Likelihood ratio test= 1.09  on 1 df,   p=0.2971
Wald test            = 0.95  on 1 df,   p=0.3293
Score (logrank) test = 1.05  on 1 df,   p=0.3053

となっていて、test1.coxph1の方にはsexの計算結果も載っているがtest1.coxph2には載っていない。何が起きているのかというと、strata関数で指定した変数ごとにデータ集合が区切られて、パラメーター推計されるようになっている。これでパラメーター推計はおしまい。


生存確率の評価にはsurvfit関数を使う。デフォルトでは各共変量の平均値を使って生存確率の値が計算される。

> summary(survfit(test1.coxph1))
Call: survfit(formula = test1.coxph1)

 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    1      7       1   0.8821   0.115     6.84e-01            1
    2      5       2   0.6095   0.208     3.13e-01            1
    3      3       1   0.4455   0.217     1.71e-01            1
    4      1       1   0.0329   0.121     2.36e-05            1
> summary(survfit(test1.coxph2))
Call: survfit(formula = test1.coxph2)

                sex=0 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    1      4       1    0.844   0.164      0.57707            1
    3      2       1    0.627   0.296      0.24861            1
    4      1       1    0.106   0.207      0.00236            1

                sex=1 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      2.0000       3.0000       2.0000       0.3337       0.2765       0.0658       1.0000

生存確率 S(t)は累積ハザード確率 H(t)と H(t) = -\log(S(t))という関係で結ばれているので、累積ハザード確率を取得する関数basehazを使って

> exp(-basehaz(test1.coxph1)$hazard)
[1] 0.88211123 0.60951279 0.44545950 0.03285681

としても計算できる（test1.coxph1の計算結果survival列の値と一致）。




(追加更新予定）

参考LINK