Cox回帰

Rで生存時間解析をするためのパッケージ survival をあらかじめ読み込んでおきます：

library(survival)

このとき survival 中で必要なパッケージ splines も読み込まれます。

詳しくはパッケージ survival のマニュアル（PDF）を読めばいいのですが，ここではCox（コックス）回帰のための関数 coxph() について解説します。

Cox回帰（Cox比例ハザードモデル）は，ハザード（平たくいえば単位時間の死亡率）が時間 $t$ と例えば二つの変数 $x_1$，$x_2$ の関数であるとき，次のように時間の関数とそれ以外の変数の線形結合の指数関数に因数分解できると仮定します：

\[ h(t, x_1, x_2) = h_0(t) \exp(\beta_1 x_1 + \beta_2 x_2) \]

この左辺の関数をハザード関数（hazard function），右辺の時間だけに依存する因子 $h_0(t)$ を基準ハザード関数（baseline hazard function）と呼びます。

具体的に次のような簡単なデータ（マニュアルの coxph() の解説で使われています）を仮定しましょう：

timestatusxsex
4100
3120
1110
1010
2111
2101
3001

これを test1 というオブジェクトに読み込んでください。

このデータで，x は何らかの変数，sex は性別，time は時間で，status が1なら死亡を意味します。例えば最初のケースは4時間後に死亡しました。status が0なら生存を意味します。例えば最後のケースは3時間後にまだ生きていましたが時間切れでそれ以上観察しませんでした。

このような time と status を合わせたものをここではsurvivalオブジェクトと呼びます。Survivalオブジェクトを作るには Surv() 関数を使います。例として

Surv(test1$time, test1$status)

と打ち込んでみてください。画面に次のように現れるはずです：

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

この 3+ というのは少なくとも3時間生きていたことを表します。

さて，これを使ってCox回帰をしてみましょう：

result1 = coxph(Surv(time,status) ~ x + sex, data=test1)
summary(result1)

画面に次のように出ます：

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

Concordance= 0.571  (se = 0.206 )
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

上では基準ハザード関数は性別によらないと考えて計算しましたが，基準ハザード関数が性別によると仮定して，変数としては使わないのであれば，sex の代わりに strata(sex) とします。

Last modified: 2012-03-21 15:18:21