8. 생존분석
MK
2021 3 13
8.1 생존분석
- R에서 생존분석을 시행하는 방법은 다음과 같습니다. 먼저 survival 패키지를 불러드립니다.
## 'data.frame': 1858 obs. of 16 variables:
## $ id : num 1 1 2 2 3 3 4 4 5 5 ...
## $ study : num 1 1 1 1 1 1 1 1 1 1 ...
## $ rx : Factor w/ 3 levels "Obs","Lev","Lev+5FU": 3 3 3 3 1 1 3 3 1 1 ...
## $ sex : num 1 1 1 1 0 0 0 0 1 1 ...
## $ age : num 43 43 63 63 71 71 66 66 69 69 ...
## $ obstruct: num 0 0 0 0 0 0 1 1 0 0 ...
## $ perfor : num 0 0 0 0 0 0 0 0 0 0 ...
## $ adhere : num 0 0 0 0 1 1 0 0 0 0 ...
## $ nodes : num 5 5 1 1 7 7 6 6 22 22 ...
## $ status : num 1 1 0 0 1 1 1 1 1 1 ...
## $ differ : num 2 2 2 2 2 2 2 2 2 2 ...
## $ extent : num 3 3 3 3 2 2 3 3 3 3 ...
## $ surg : num 0 0 0 0 0 0 1 1 1 1 ...
## $ node4 : num 1 1 0 0 1 1 1 1 1 1 ...
## $ time : num 1521 968 3087 3087 963 ...
## $ etype : num 2 1 2 1 2 1 2 1 2 1 ...- colon 데이터에서 time은 생존기간을 뜻하며 status가 1은 사망한 경우이고 0은 현재 살아있는 경우 입니다.
8.2 Kaplan-Meier Curve
- Kaplan-Meier curve를 통해 추정된 생존함수는 survfit()을 사용하여 얻습니다. 전체환자의 데이터를 얻고자 할 때는 survfit 함수의 정규식에 “~1”을 붙여 사용합니다. summary로 전체 환자의 데이터를 얻을 수 있습니다.
## Call: survfit(formula = Surv(time, status == 1) ~ 1, data = colon)
##
## n events median 0.95LCL 0.95UCL
## 1858 920 2351 2018 2910
- 가운데 선이 KM curve이고 점선은 95% 신뢰구간입니다. 신뢰구간을 없애려면 conf.int 인수를 FALSE로 지정해줍니다.
- colon 데이터의 rx열은 환자를 치료군에 따라 세군으로 나눈 데이터 입니다. 치료에 따른 생존율의 차이를 보기 위해서는 survfit 함수를 사용할때 “~rx”를 넣어줍니다.
fit1 = survfit(Surv(time,status==1)~rx,data=colon)
plot(fit1,col=1:3, lty=1:3)
legend("topright", legend=levels(rx), col=1:3, lty=1:3)
- censored data의 mark를 보고싶다면 mark.time 인수를 TRUE로 줍니다.
- cumulative hazard 그래프를 보려면 fun=“cumhaz”로 줍니다.
plot(fit1,col=1:3, lty=1:3, fun="cumhaz", xlab="Days")
legend("topright", legend=levels(rx), col=1:3, lty=1:3)
8.3 Plot for survival curve
- Survival curve를 그리고 위한 몇가지 함수를 설명하겠습니다. GGally 패키지, ggkm() 함수, rms 패키지를 이용하여 생존곡선을 그리는 방법을 소개하겠습니다.
- GGally 패키지의 ggsurv() 함수를 이용하면 plot() 함수를 이용하는 것보다 보기 좋은 그래프를 얻을 수 있습니다.
require(GGally)
fit = survfit(Surv(time,status==1)~rx, data=colon)
ggsurv(fit, plot.cens=FALSE) # censored data mark 없음
- ggkm() 함수는 보기 좋은 그래프와 함께 number at risk table도 같이 plot에 넣을 수 있습니다.
require(ggplot2)
require(plyr)
# 제대로 작동하려면 인터넷 상에서 ggkm.R을 다운받고 실행합시다.
# https://gist.github.com/araastat/9927677
# source("ggkm.R")
ggkm <- function(sfit, returns = FALSE,
xlabs = "Time", ylabs = "survival probability",
ystratalabs = NULL, ystrataname = NULL,
timeby = 100, main = "Kaplan-Meier Plot",
pval = TRUE, ...) {
require(plyr)
require(ggplot2)
require(survival)
require(gridExtra)
if(is.null(ystratalabs)) {
ystratalabs <- as.character(levels(summary(sfit)$strata))
}
m <- max(nchar(ystratalabs))
if(is.null(ystrataname)) ystrataname <- "Strata"
times <- seq(0, max(sfit$time), by = timeby)
.df <- data.frame(time = sfit$time, n.risk = sfit$n.risk,
n.event = sfit$n.event, surv = sfit$surv,
strata = summary(sfit, censored =T)$strata,
upper = sfit$upper, lower = sfit$lower)
levels(.df$strata) <- ystratalabs
zeros <- data.frame(time = 0, surv = 1, strata = factor(ystratalabs, levels=levels(.df$strata)),
upper = 1, lower = 1)
.df <- rbind.fill(zeros, .df)
d <- length(levels(.df$strata))
p <- ggplot(.df, aes(time, surv, group = strata)) +
geom_step(aes(linetype = strata), size = 0.7) +
theme_bw() +
theme(axis.title.x = element_text(vjust = 0.5)) +
scale_x_continuous(xlabs, breaks = times, limits = c(0, max(sfit$time))) +
scale_y_continuous(ylabs, limits = c(0, 1)) +
theme(panel.grid.minor = element_blank()) +
theme(legend.position = c(ifelse(m < 10, .28, .35), ifelse(d < 4, .25, .35))) +
theme(legend.key = element_rect(colour = NA)) +
labs(linetype = ystrataname) +
theme(plot.margin = unit(c(0, 1, .5, ifelse(m < 10, 1.5, 2.5)), "lines")) +
ggtitle(main)
if(pval) {
sdiff <- survdiff(eval(sfit$call$formula), data = eval(sfit$call$data))
pval <- pchisq(sdiff$chisq, length(sdiff$n)-1, lower.tail = FALSE)
pvaltxt <- ifelse(pval < 0.0001, "p < 0.0001", paste("p =", signif(pval, 3)))
p <- p + annotate("text", x = 0.6 * max(sfit$time), y = 0.1, label = pvaltxt)
}
## Plotting the graphs
print(p)
if(returns) return(p)
}
- rms 패키지의 survplot()을 이용할 수도 있습니다.
require(rms)
S=Surv(time,status==1)
f=npsurv(S~rx)
survplot(f,xlab="Days", conf="none", n.risk=TRUE, y.n.risk=-0.3, cex.n.risk=0.6, col=3:1) # 95% 신뢰구간 표시 안함, number at risk 표시
8.4 Log-rank test
- 환자들의 생존확률이 치료에 따라 유의한 차이가 있는 지 검정하는 것이 log-rank test 입니다. survdiff() 함수를 통해 log-rank test를 할 수 있습니다.
## Call:
## survdiff(formula = Surv(time, status == 1) ~ rx, data = colon)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## rx=Obs 630 345 299 7.01 10.40
## rx=Lev 620 333 295 4.93 7.26
## rx=Lev+5FU 608 242 326 21.61 33.54
##
## Chisq= 33.6 on 2 degrees of freedom, p= 5e-088.5 Cox Regression
- Cox regression은 log-hazard를 회귀분석처럼 선형모형화하여 위험함수 h(t)를 추정합니다. coxph() 함수를 사용합시다.
## Call:
## coxph(formula = Surv(time, status == 1) ~ rx, data = colon)
##
## coef exp(coef) se(coef) z p
## rxLev -0.02090 0.97932 0.07683 -0.272 0.786
## rxLev+5FU -0.44101 0.64339 0.08391 -5.256 1.47e-07
##
## Likelihood ratio test=35.23 on 2 df, p=2.233e-08
## n= 1858, number of events= 920- Observation한 군을 기준으로 Lev를 사용한 군의 hazard ratio(HR) exp(coef)는 p-value가 의미 없지만 Lev+5-FU를 사용한 군은 HR 0.643 p-value가 의미있는 차이를 보입니다.
- HR의 95% 신뢰구간을 알고 싶으면 summary() 함수를 사용합니다. Observation한 군을 기준으로 Lev를 사용한 군의 HR의 95% 신뢰구간은 0.842~1.138이고 Lev+5-FU를 사용한 군의 HR의 95% 신뢰구간은 0.546~0.758입니다.
## Call:
## coxph(formula = Surv(time, status == 1) ~ rx, data = colon)
##
## n= 1858, number of events= 920
##
## coef exp(coef) se(coef) z Pr(>|z|)
## rxLev -0.02090 0.97932 0.07683 -0.272 0.786
## rxLev+5FU -0.44101 0.64339 0.08391 -5.256 1.47e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## rxLev 0.9793 1.021 0.8424 1.1385
## rxLev+5FU 0.6434 1.554 0.5458 0.7584
##
## Concordance= 0.545 (se = 0.009 )
## Likelihood ratio test= 35.23 on 2 df, p=2e-08
## Wald test = 33.11 on 2 df, p=6e-08
## Score (logrank) test = 33.63 on 2 df, p=5e-088.6 Cox 비례위험모형의 자동화
- 예후와 관련이 있어 보이는 변수들이 여러개 있습니다. Colon 데이터의 경우 17개의 변수가 존재하는 데 mycph() 함수를 통해 각각의 HR를 계산할 수 있습니다.
require(moonBook)
colon$TS = Surv(time,status==1)
out=mycph(TS~.-id-study-time-status-etype,data=colon)##
## mycph : perform coxph of individual expecting variables
##
## Call: TS ~ . - id - study - time - status - etype, data= colon## HR lcl ucl p
## rxLev 0.98 0.84 1.14 0.786
## rxLev+5FU 0.64 0.55 0.76 0.000
## sex 0.97 0.85 1.10 0.610
## age 1.00 0.99 1.00 0.382
## obstruct 1.27 1.09 1.49 0.003
## perfor 1.30 0.92 1.85 0.142
## adhere 1.37 1.16 1.62 0.000
## nodes 1.09 1.08 1.10 0.000
## differ 1.36 1.19 1.55 0.000
## extent 1.78 1.53 2.07 0.000
## surg 1.28 1.11 1.47 0.001
## node4 2.47 2.17 2.83 0.0002.그래프는 HRplot() 함수를 사용해서 그립니다.
- 너무 변수가 많을 경우유의한 것만 보기위해서는 sig.level을 지정해주면 됩니다.
8.7 rms 패키지를 이용한 Cox 비례위험모형의 시각화
- 이 패키지에 있는 함수를 이용하여 cph 모형을 시각화하는 방법을 소개하고자 합니다. 폐암데이터를 이용해 보겠습니다. 폐암 데이터 중 성별이 남자는 1, 여자는 2로 되어 있으니 그래프에서 알아보기 쉽도록 Male, Female로 변경하겠습니다.
require(survival)
require(rms)
data(lung)
lung$sex <- factor(lung$sex, levels = 1:2, labels=c("Male", "Female"))
TS <- with(lung, Surv(time,status)) # 객체 TS 생성
ddist <- datadist(lung) # rms분석준비
options(datadist="ddist")fit = coxph(TS~rcs(age,4)+sex, lung, x=T, y=T)
par(mfrow=c(1,2))
termplot(fit, se=T, rug=T, ylab=rep("Hazard Ratio", times=2),
main= rep("cph() plot", times=2),
col.se=rgb(.2,.2,1,.4), col.term = "black")
- termplot2() 함수는 조금 더 개선되었습니다.
par(mfrow=c(1,2))
termplot(fit, se=T, rug.type="density", rug=T, desity.proportion=.05,
se.type="polygon",
ylab=rep("Hazard Ratio", times=2),
main= rep("cph() plot", times=2),
col.se=rgb(.2,.2,1,.4), col.term = "black")
- rms 패키지의 cph() 함수 및 predict() 함수를 사용하여 age와 sex의 HR을 같이 그려볼 수 있습니다. cph() 함수는 coxph()함수와 비슷하다고 이해합시다.
## Cox Proportional Hazards Model
##
## cph(formula = TS ~ rcs(age, 4) + sex, data = lung, x = T, y = T)
##
## Model Tests Discrimination
## Indexes
## Obs 228 LR chi2 15.87 R2 0.067
## Events 165 d.f. 4 Dxy 0.184
## Center 1.7005 Pr(> chi2) 0.0032 g 0.359
## Score chi2 16.15 gr 1.431
## Pr(> chi2) 0.0028
##
## Coef S.E. Wald Z Pr(>|Z|)
## age 0.0351 0.0355 0.99 0.3220
## age' -0.0643 0.0705 -0.91 0.3616
## age'' 0.3937 0.3430 1.15 0.2511
## sex=Female -0.5055 0.1678 -3.01 0.0026
## 
8.8 Bladder 데이터를 이용한 생존분석
## id rx number size stop event enum
## 1 1 1 1 3 1 0 1
## 2 1 1 1 3 1 0 2
## 3 1 1 1 3 1 0 3
## 4 1 1 1 3 1 0 4
## 5 2 1 2 1 4 0 1
## 6 2 1 2 1 4 0 2## Call:
## coxph(formula = Surv(stop, event) ~ rx + number + size + enum,
## data = bladder)
##
## n= 340, number of events= 112
##
## coef exp(coef) se(coef) z Pr(>|z|)
## rx -0.59739 0.55024 0.20088 -2.974 0.00294 **
## number 0.21754 1.24301 0.04653 4.675 2.93e-06 ***
## size -0.05677 0.94481 0.07091 -0.801 0.42333
## enum -0.60385 0.54670 0.09401 -6.423 1.34e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## rx 0.5502 1.8174 0.3712 0.8157
## number 1.2430 0.8045 1.1347 1.3617
## size 0.9448 1.0584 0.8222 1.0857
## enum 0.5467 1.8291 0.4547 0.6573
##
## Concordance= 0.753 (se = 0.021 )
## Likelihood ratio test= 67.21 on 4 df, p=9e-14
## Wald test = 64.73 on 4 df, p=3e-13
## Score (logrank) test = 69.42 on 4 df, p=3e-14
## Call:
## coxph(formula = Surv(stop, event) ~ strata(rx) + number + size +
## enum, data = bladder)
##
## n= 340, number of events= 112
##
## coef exp(coef) se(coef) z Pr(>|z|)
## number 0.21371 1.23826 0.04648 4.598 4.27e-06 ***
## size -0.05485 0.94662 0.07097 -0.773 0.44
## enum -0.60695 0.54501 0.09408 -6.451 1.11e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## number 1.2383 0.8076 1.1304 1.3564
## size 0.9466 1.0564 0.8237 1.0879
## enum 0.5450 1.8348 0.4532 0.6554
##
## Concordance= 0.74 (se = 0.023 )
## Likelihood ratio test= 61.84 on 3 df, p=2e-13
## Wald test = 60.04 on 3 df, p=6e-13
## Score (logrank) test = 65.05 on 3 df, p=5e-14plot(survfit(cox2), xlab="Days", ylab="Survival Rate", conf.int=TRUE, col=1:2)
legend("bottomleft", legend=c(1,2), lty = 1, col = 1:2, text.col = 1:2, title = 'rx')
- 비례위험에 대한 가정 확인
## chisq df p
## rx 0.051 1 0.821
## number 3.235 1 0.072
## size 3.027 1 0.082
## enum 28.701 1 8.4e-08
## GLOBAL 32.868 4 1.3e-06## chisq df p
## number 3.24 1 0.072
## size 2.90 1 0.089
## enum 28.86 1 7.8e-08
## GLOBAL 32.92 3 3.3e-07- Anderson-Gill analysis는 생존분석과 유사하기는 하지만, 구간화된 데이터를 받고, 여러개의 이벤트를 다룰 수 있다는 점이 다릅니다.
ag1 <- coxph(Surv(start,stop,event)~rx+number+size+enum+cluster(id), data=bladder2)
ag2 <- coxph(Surv(start,stop,event)~strata(rx)+number+size+enum+cluster(id), data=bladder2)
plot(survfit(ag2), conf.int=TRUE, col=1:2)
legend("bottomleft", legend=c(1,2), lty = 1, col = 1:2, text.col = 1:2, title = 'rx')