Cox回帰（比例ハザードモデル）
R を使って、生存分析でもっともよく使われるCox回帰（比例ハザードモデル）をやってみる。基本的にSASのノート "Survival Analysis Using the Proportional Hazards Model" を参考にして、そのなかで使われているのと同じ、ヘロイン依存患者のメタドン治療についてのデータを使った（データはこの辺でも入手できる）。アウトカム（イベント）はクリニックでのメタドン治療からの脱落、共変量（独立変数）は治療を受けたクリニック（2ヶ所）、服役歴の有無、メタドンの用量（連続変数）。
まずはデータを読み込んで変数を名付けたら、イベント数の確認。
> meth <- read.table("methadone.txt")
> names(meth) <- c("id","clinic","status","time","prison","dose")
> meth[1:5,]
id clinic status time prison dose
1 1 1 1 428 0 50
2 2 1 1 275 1 55
3 3 1 1 262 0 55
4 4 1 1 183 0 30
5 5 1 1 259 1 65
> table(meth$status)
0 1
88 150イベント数（メタドン治療からの脱落者数）は150に対して、打ち切り例数（誤解しやすいけど、この場合は非脱落者数）は88。Cox回帰を使う1つの目安に「共変量1つにつき10のイベント数（あるいは非イベント数）」というのがあるけど、ここでは大丈夫みたい。実際にCox回帰を行なうには survival パッケージの coxph() を使う。
> library(survival)
> meth.cox <- coxph(Surv(time,status)~clinic+dose+prison,meth,method="efron")
> summary(meth.cox)
Call:
coxph(formula = Surv(time, status) ~ clinic + dose + prison,
data = meth, method = "efron")
n= 238
coef exp(coef) se(coef) z p
clinic -1.0099 0.364 0.21489 -4.70 2.6e-06
dose -0.0354 0.965 0.00638 -5.54 2.9e-08
prison 0.3266 1.386 0.16722 1.95 5.1e-02
exp(coef) exp(-coef) lower .95 upper .95
clinic 0.364 2.745 0.239 0.555
dose 0.965 1.036 0.953 0.977
prison 1.386 0.721 0.999 1.924
Rsquare= 0.238 (max possible= 0.997 )
Likelihood ratio test= 64.6 on 3 df, p=6.23e-14
Wald test = 54.1 on 3 df, p=1.06e-11
Score (logrank) test = 56.3 on 3 df, p=3.6e-12イベントが同時発生した場合（タイがある場合）の取り扱いについては method で Efron 法、Breslow 法、または Exact を指定できる。デフォルトは Efron で、上では明示的にこれを指定している。結果をみるとクリニックと用量で有意、服役歴については微妙なところ。また回帰係数の指数をとったハザード比も計算されている。クリニックを比較した場合、クリニック2では脱落のハザードがクリニック1の36.4%、またメタドンの用量が1単位増えるたびにハザードは3.5% 減っていくことになる。
coxph() の結果をプロットすると、各共変量が平均値で調整された生存曲線（とその信頼区間）が描かれる。
> plot(survfit(meth.cox),xlab="Days Spent in Clinic", + ylab="Survival Function Estimate",main="Adjusted Survival Function")
また例えばクリニックごとに調整された生存曲線をプロットするには、変数ごとにハザード関数に放り込む値を指定したデータフレームを作ってやって、survfit() の引数 newdata に指定してからプロットする。
> attach(meth)
> meth.adj <- data.frame(clinic=c(1,2),prison=rep(mean(prison),2),dose=rep(mean(dose),2))
> detach(meth)
> plot(survfit(meth.cox,newdata=meth.adj),lty=1:2,col=1:2,xlab="Days Spent in Clinic",
+ ylab="Survival Function Estimate",main="Adjusted Survival Function by Clinic")
> legend(650,1,lty=1:2,col=1:2,legend=c("Clinic 1","Clinic 2"))
> abline(h=0.5,lty=3,col=3)
同様にして、「クリニック2で治療を受け、服役歴があり、メタドン用量が90 mg/day」の患者の生存曲線をプロットしてみる。
> meth.pred <- data.frame(clinic=2,prison=1,dose=90) > plot(survfit(meth.cox,newdata=meth.pred),xlab="Days Spent in Clinic", + ylab="Survival Function Estimate") > abline(v=200,lty=3,col=3) > (data.frame(pred$time,pred$surv,pred$lower,pred$upper))[45:50,] pred.time pred.surv pred.lower pred.upper 45 181 0.9592266 0.9344855 0.9846226 46 183 0.9581918 0.9329043 0.9841647 47 190 0.9571393 0.9313015 0.9836939 48 192 0.9560832 0.9296940 0.9832215 49 193 0.9550238 0.9280821 0.9827475 50 204 0.9539589 0.9264626 0.9822711
この場合、その患者が200日間脱落せずにいられる確率は95%ほどで、その95%信頼区間は92-98%であることが分かる。
さて回帰診断だけど、Cox回帰のいちばんの前提は、ハザード比が時間によらず一定であること。この比例ハザードの仮定を吟味するのには cox.zph() の関数を使うのが簡単。どうやらこれはモデルに時間の関数を組み込んだときに、その関数にかかる係数が0であるかどうかを検定しているらしい。この辺参照。
> meth.zph <- cox.zph(meth.cox)
> meth.zph
rho chisq p
clinic -0.2578 11.185 0.000824
dose 0.0724 0.700 0.402749
prison -0.0382 0.220 0.639369
GLOBAL NA 12.616 0.005546結果を見ると、クリニックの変数で有意となり、比例ハザードの仮定が充たされていない様子。さらにそのアウトプットをそのままプロットすると、変数ごとの Scaled Schoenfeld residual plot が吐きだされる。
> windows(width=11,height=4) > par(mfrow=c(1,3)) > plot(meth.zph)
ただしこのプロットはやや見にくいので、自分でScaled/unscaled Schoenfeld 残差をプロットしてみよう。それぞれの残差は residual() で引数 type に scaledsch あるいは shoenfeld を指定して抽出できる。Schoenfeld 残差はイベント例にしか計算されないことに注意しよう。
> ##### Schoenfeld residuals against survival time
> windows(width=11,height=8)
> par(mfrow=c(2,3),cex=.8)
> scshoen <- resid(meth.cox,type="scaledsch")
> time.sort <- sort(meth$time[meth$status==1])
> xnames <- names(meth)[c(2,5,6)]
> for (i in 1:3){
+ plot(time.sort,scshoen [,i],xlab="Days",
+ ylab=paste("Scaled Schoenfeld Residuals for",xnames[i]))
+ abline(h=0,lty=2,col=3)
+ lines(lowess(time.sort,scshoen [,i],0.9),col=2)}
> schoen <- resid(meth.cox,type="schoenfeld")
> for (i in 1:3){
+ plot(time.sort,schoen[,i],xlab="Days",
+ ylab=paste("Schoenfeld Residuals for",xnames[i]))
+ abline(h=0,lty=2,col=3)
+ lines(lowess(time.sort,schoen[,i],0.9),col=2)}
Schoenfeld 残差プロットにパターンが見られなければ比例ハザードの仮定には問題はない。プロットをみると、クリニックの変数でやはり右下がりの傾向がみられて、比例ハザードが怪しいことが分かる。
Schoenfeld 残差プロットでもいいのだけど、比例ハザードを確かめるより一般的なやりかたは、生存関数の二重対数プロット（ログマイナスログ）をとってやること。これにはまず各共変量を Strata() で層別して、横軸に対数時間をとってプロットしてやる。ただし、用量のように共変量が連続変数の場合、適当にデータを区切って層別する必要がある。下のコードでは二重対数プロットの関数を定義してある。
> ##### Log-log survival functions
> meth.clinic <- coxph(Surv(time,status)~strata(clinic)+dose+prison,meth)
> meth.prison <- coxph(Surv(time,status)~clinic+dose+strata(prison),meth)
> dose.strata <- cut(meth$dose,c(0,50,60,70,120))
> meth.dose <- coxph(Surv(time,status)~clinic+strata(dose.strata)+prison,meth)
> loglog.plot <- function(coxph, title=NULL) {
+ group <- rep(1:length(survfit(coxph)$strata), survfit(coxph)$ntimes.strata)
+ plot(log(survfit(coxph)$time), log(-log(survfit(coxph)$surv)),type='n',
+ xlab="Log(time)",ylab="Log(-log(S(t)))")
+ for (i in 1:length(survfit(coxph)$strata)) {
+ lines(log(survfit(coxph)$time)[(group==i)],
+ log(-log(survfit(coxph)$surv))[(group==i)],type="l",lty=i,lwd=2,col=i)
+ }
+ if (!is.null(title)) {
+ title(title)
+ }
+ }
> windows(width=11,height=4)
> par(mfrow=c(1,3))
> loglog.plot(meth.clinic,"Log-log Survival Functions for Clinic")
> legend(5.3,-4,lty=1:2,col=1:2,legend=c("Clinic 1","Clinic 2"))
> loglog.plot(meth.prison,"Log-log Survival Functions for Prison")
> legend(4.3,-4,lty=1:2,col=1:2,legend=c("Prison Record: No","Prison Record: Yes"))
> loglog.plot(meth.dose,"Log-log Survival Functions for Dose")
> legend(5.4,-3,lty=1:4,col=1:4,legend=c("Dose 1","Dose 2","Dose 3","Dose 4"))
比例ハザードの仮定が充たされているならば、ほぼ平行な直線がプロットされるはず。プロットをみると、クリニックの変数で明らかに傾きが違っているけれど、用量と服役歴のグラフでははっきりしない。
また共変量が連続変数の場合、対数ハザードとその共変量との関係が線形であるという仮定がある。この線形の仮定を吟味するのに手っ取り早いのはMartingale 残差プロットとその成分プラス残差プロットをとってやること。プロットを lowess() で平滑化した直線がゼロの周辺、あるいは回帰直線上に乗っていれば線形の仮定に問題はない。
> ##### Martingale residual plot and Component-plus-residual plot > windows(width=7.5,height=4) > par(mfrow=c(1,2)) > martin <- resid(meth.cox,type='martingale') > plot(meth$dose,martin,xlab="Dose",ylab="Martingale Residual") > abline(h=0,lty=2,col=3) > lines(lowess(meth$dose,martin,iter=0),col=2) > plot(meth$dose,coef(meth.cox)[2]*meth$dose+martin,xlab="Dose",ylab="Component+Residual") > abline(lm(coef(meth.cox)[2]*meth$dose+martin~meth$dose),lty=2,col=3) > lines(lowess(meth$dose,coef(meth.cox)[2]*meth$dose+martin,iter=0),col=2)
どうやら問題はないみたい。線形の仮定をチェックするもう1つのやりかたは、共変量をいくつかに区切ってダミー変数を作り、Cox回帰を行なったときの係数がほぼ線形になっていることを確認する。ただしサンプルがあまり大きくない場合、ダミー変数の数に注意。「共変量1つにつき10のイベント数（あるいは非イベント数）」の目安を忘れずに。
> levels(dose.strata) <- 1:4
> x <- matrix(0,ncol=4,nrow=nrow(meth))
> for (i in 1:nrow(x)){
+ if (dose.strata[i]==1) {x[i,1] <- 1}
+ else if (dose.strata[i]==2) {x[i,2] <- 1}
+ else if (dose.strata[i]==3) {x[i,3] <- 1}
+ else if (dose.strata[i]==4) {x[i,4] <- 1}
+ }
> meth <- data.frame(meth,x)
> names(meth)[7:10] <- c("dose1","dose2","dose3","dose4")
> meth.dose <- coxph(Surv(time,status)~clinic+prison+dose1+dose2+dose3+dose4,meth)
> summary(meth.dose)
Call:
coxph(formula = Surv(time, status) ~ clinic + prison + dose1 +
dose2 + dose3 + dose4, data = meth)
n= 238
coef exp(coef) se(coef) z p
clinic -0.951 0.386 0.219 -4.35 1.4e-05
prison 0.303 1.354 0.167 1.81 7.0e-02
dose1 1.451 4.268 0.292 4.96 7.0e-07
dose2 1.209 3.350 0.285 4.24 2.2e-05
dose3 1.093 2.982 0.311 3.52 4.3e-04
dose4 NA NA 0.000 NA NA
exp(coef) exp(-coef) lower .95 upper .95
clinic 0.386 2.589 0.251 0.593
prison 1.354 0.739 0.976 1.879
dose1 4.268 0.234 2.406 7.570
dose2 3.350 0.299 1.916 5.857
dose3 2.982 0.335 1.623 5.481
dose4 NA NA NA NA
Rsquare= 0.242 (max possible= 0.997 )
Likelihood ratio test= 65.9 on 5 df, p=7.34e-13
Wald test = 50.8 on 5 df, p=9.68e-10
Score (logrank) test = 57.1 on 5 df, p=4.9e-11上の結果では dose4 が基準となって、係数はほぼ直線上に並んでいると言っていいだろう。
またはずれ値や影響の大きな観測値を見つけるには DfBeta をプロットするか、連続変数に対してはscore residual をプロットしてやるといい。DfBeta もscore residual も residual() で引数の type に指定してやることで抽出できる。
> ##### Detecting influential observations
> dfbeta <- resid(meth.cox,type="dfbeta")
> windows(width=11,height=4)
> par(mfrow=c(1,3),cex=1)
> for (i in 1:3){
+ plot(dfbeta[,i],ylab=paste("dfbeta for",xnames[i]))
+ abline(h=0,lty=2)
+ }
> score <- resid(meth.cox,type="score")
> plot(meth$dose,score[,2],xlab="Dose",ylab="Score Residual")
いずれの場合も、他から飛びぬけている観測値が影響の大きな観測値である。
