[R]R言語で生存時間分析
生存時間分析
ある時点から興味のあるイベントが起きるまでの時間とイベントとの間の関係に関する分析。
以下のような解析を行う
- 生存率の推定（Kaplan-Meier曲線）
- 2郡の生存率の比較 （Log-Rank検定）
- 生存率と共変量との関係の解明 (Cox回帰)
また、イベントの例としては、以下のようなものがある。
- 機械システムや製品の故障
- 疾患の病気の再発や死亡
データの準備
ここでは、survivalパッケージにあるデータセットcolonを用いる。
- colon : Stage B/Cの結腸癌患者を対象とした術後補助化学療法の比較臨床試験データ
> library(survival) 要求されたパッケージ splines をロード中です > colon.OS <- subset(colon, colon$etype==2) #死亡に関するデータのみを取り出す
Kaplan-Meier法による生存時間分布の推定
survfit関数を使って、生存時間分布を推定することが出来る
- 生存時間分布の推定
> result.KM <- survfit(Surv(time, status)~rx, data=colon.OS) > result.KM Call: survfit(formula = Surv(time, status) ~ rx, data = colon.OS) records n.max n.start events median 0.95LCL 0.95UCL rx=Obs 315 315 315 168 2083 1656 2789 rx=Lev 310 310 310 161 2152 1540 NA rx=Lev+5FU 304 304 304 123 NA 2725 NA
- Kaplan-Meier曲線
> plot(result.KM, lty=1:3, xlim=c(0, 3500), ylim=c(0,1), xlab="Time", ylab="Event Free Fate", main="Survival Curve") > abline(0.5, 0) > plot(result.KM, lty=1:3, xlim=c(0, 3500), ylim=c(0,1), xlab="Time", ylab="Event Free Fate", main="Survival Curve") > abline(0.5, 0)
- Log-Logプロット
> plot(result.KM, lty=1:3, fun="cloglog", xlab="Time", ylab="log(-log(Event Free Rate))", main="Log-log Curve") > legend(20, 0, legend=c("Obs", "Levamisole", "Levamisole+5-FU"), lty=1:3, bty="n")
Log-Rank検定
2郡以上の観測値が得られた場合、郡ごとに生存曲線の間の有意性について検定する必要がある。
最も広く用いられている検定法は、Log-Rank検定法で、survdiff関数で実行する事が出来る。
> result.LR <- survdiff(Surv(time, status)~rx, + subset=(rx %in% c("Obs", "Lev+5FU")), data=colon.OS) > result.LR Call: survdiff(formula = Surv(time, status) ~ rx, data = colon.OS, subset = (rx %in% c("Obs", "Lev+5FU"))) N Observed Expected (O-E)^2/E (O-E)^2/V rx=Obs 315 168 141 5.12 9.97 rx=Lev+5FU 304 123 150 4.82 9.97 Chisq= 10 on 1 degrees of freedom, p= 0.00159
Cox比例ハザードモデルによるハザード比の推定
Cox比例ハザードモデルは、イベントに影響を及ぼす複数の因子（共変量）の影響を解析する代表的なモデル。
Cox比例ハザードモデルのパラメータ推定するには、coxph関数を用いる。method引数で、推定法を指定する事が出来、ここでは、ブレスローの近似法を指定している。
> result.Cox <- coxph(Surv(time, status)~rx+node4+age+obstruct+perfor+adhere+nodes+differ+extent+surg+node4, method="breslow", data=na.omit(colon.OS)) > summary(result.Cox) Call: coxph(formula = Surv(time, status) ~ rx + node4 + age + obstruct + perfor + adhere + nodes + differ + extent + surg + node4, data = na.omit(colon.OS), method = "breslow") n= 888, number of events= 430 coef exp(coef) se(coef) z Pr(>|z|) rxLev -0.041393 0.959452 0.114366 -0.362 0.717404 rxLev+5FU -0.362772 0.695745 0.122066 -2.972 0.002959 ** node4 0.676559 1.967096 0.142361 4.752 2.01e-06 *** age 0.007537 1.007565 0.004181 1.803 0.071442 . obstruct 0.268871 1.308486 0.120310 2.235 0.025430 * perfor 0.017162 1.017310 0.269935 0.064 0.949307 adhere 0.170110 1.185435 0.131738 1.291 0.196611 nodes 0.043738 1.044709 0.015239 2.870 0.004102 ** differ 0.137941 1.147908 0.100818 1.368 0.171244 extent 0.446006 1.562061 0.118344 3.769 0.000164 *** surg 0.240548 1.271946 0.106118 2.267 0.023402 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 exp(coef) exp(-coef) lower .95 upper .95 rxLev 0.9595 1.0423 0.7668 1.2005 rxLev+5FU 0.6957 1.4373 0.5477 0.8838 node4 1.9671 0.5084 1.4882 2.6002 age 1.0076 0.9925 0.9993 1.0159 obstruct 1.3085 0.7642 1.0336 1.6564 perfor 1.0173 0.9830 0.5994 1.7267 adhere 1.1854 0.8436 0.9157 1.5347 nodes 1.0447 0.9572 1.0140 1.0764 differ 1.1479 0.8712 0.9421 1.3987 extent 1.5621 0.6402 1.2387 1.9698 surg 1.2719 0.7862 1.0331 1.5660 Rsquare= 0.145 (max possible= 0.998 ) Likelihood ratio test= 139.5 on 11 df, p=0 Wald test = 148.3 on 11 df, p=0 Score (logrank) test = 158.8 on 11 df, p=0
比例ハザード性の仮定が成り立っているかを確認
Cox比例ハザードモデルは、ハザード比が時間によらず一定であることを前提としているため、この仮定を吟味する必要がある。
これは、cox.zph関数を実行することで出来る。
> result.coxzph <- cox.zph(result.Cox) > result.coxzph rho chisq p rxLev -0.0930 3.8661 4.93e-02 rxLev+5FU -0.0711 2.2282 1.36e-01 node4 -0.1327 7.0797 7.80e-03 age -0.0626 1.8782 1.71e-01 obstruct -0.1496 10.1433 1.45e-03 perfor 0.0534 1.3010 2.54e-01 adhere 0.0326 0.4741 4.91e-01 nodes 0.1013 3.5367 6.00e-02 differ -0.1588 12.7542 3.55e-04 extent -0.0778 2.5097 1.13e-01 surg -0.0118 0.0597 8.07e-01 GLOBAL NA 37.5164 9.44e-05
また、plot関数で、比例性の診断プロットを出力することが出来る。
スプライン平滑か曲線の傾向に、時間に伴う明らかな変化パターンが見られなければ、比例ハザードの仮定には問題がないと言われている。
個別の変数に比例ハザードの仮定が充たされていない場合は、変数の交差項をモデルに導入するような試みが必要になる。
> op <- par(mfrow=c(2,2), mar=c(4.5, 4, 1, 1)) > plot(result.coxzph, df=2); par(op)
変数選択
回帰分析と同様に、Cox比例ハザードモデルでもsetAIC関数でモデルを構築出来る。
> library(MASS) > result.step <- stepAIC(result.Cox) ...（省略）... > summary(result.step) Call: coxph(formula = Surv(time, status) ~ rx + node4 + age + obstruct + nodes + differ + extent + surg, data = na.omit(colon.OS), method = "breslow") n= 888, number of events= 430 coef exp(coef) se(coef) z Pr(>|z|) rxLev -0.038729 0.962011 0.114243 -0.339 0.73461 rxLev+5FU -0.364097 0.694824 0.121994 -2.985 0.00284 ** node4 0.674685 1.963414 0.141758 4.759 1.94e-06 *** age 0.008063 1.008096 0.004159 1.939 0.05253 . obstruct 0.270757 1.310957 0.119359 2.268 0.02330 * nodes 0.043733 1.044703 0.015168 2.883 0.00394 ** differ 0.152120 1.164300 0.100346 1.516 0.12953 extent 0.461643 1.586679 0.118520 3.895 9.82e-05 *** surg 0.241874 1.273633 0.106109 2.279 0.02264 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 exp(coef) exp(-coef) lower .95 upper .95 rxLev 0.9620 1.0395 0.7690 1.2034 rxLev+5FU 0.6948 1.4392 0.5471 0.8825 node4 1.9634 0.5093 1.4871 2.5922 age 1.0081 0.9920 0.9999 1.0163 obstruct 1.3110 0.7628 1.0375 1.6565 nodes 1.0447 0.9572 1.0141 1.0762 differ 1.1643 0.8589 0.9564 1.4174 extent 1.5867 0.6302 1.2578 2.0016 surg 1.2736 0.7852 1.0345 1.5681 Rsquare= 0.144 (max possible= 0.998 ) Likelihood ratio test= 137.9 on 9 df, p=0 Wald test = 146.7 on 9 df, p=0 Score (logrank) test = 156.7 on 9 df, p=0
生存時間の推定
構築したモデルによる生存時間の当てはめは、survfit関数を用いて行うことが出来る。
また、plot関数で生存曲線および信頼区間を出力することが出来る。
> result.step.fit <- survfit(result.step) > plot(result.step.fit)
