Survival Analysis-II
Sung Woo Ko
2022년 1월 14일
R을 이용한 Cox 비례 위험 모형 적합
1. Cox 비례 위험 모형 (Cox Propotional Hazard Model): ’coxph’함수 사용
’Cox 비례 위험 모형’이란?
여러 독립변수들이 생존함수에 미치는 영향을 조사하기 위한 다변수 분석법(multivariable analysis)
library(survival) # survival package 불러옴
library(survminer)## 필요한 패키지를 로딩중입니다: ggplot2## 필요한 패키지를 로딩중입니다: ggpubrcolon<-get(data(colon))
str(colon) # 'str'은 데이터의 구조를 보는 명령어## '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 ...R 내장 데이터 ’colon’을 불러온다
colon 데이터의 구조를 확인한다
colon 데이터에서 rx(치료)는 ‘obs(관찰)’, ‘Lev’, ‘Lev+ 5FU’ 3군으로 이루어짐
result1<-coxph(Surv(time,status)~rx,data=colon) # 'rx(치료)'에 따른 위험비를 알아본다
result1## Call:
## coxph(formula = Surv(time, status) ~ 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
Cox 비례 위험 모형을 적합 하기 위해서 coxph 함수를 사용한다
coxph 결과에서 HR(Hazard ratio)를 알고 싶으면 ‘exp(coef)’를 보면된다
해석
obs에 비해서 ’Lev’치료군의 HR는 0.978 (p 0.786)로 유의한 차이는 없다
obs에 비해서 ’Lev+5FU’치료군의 HR는 0.643 (p 1.5x 10 -7)로 유의한 차이를 보인다
summary(result1)## Call:
## coxph(formula = Surv(time, status) ~ 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-08
HR의 95% 신뢰 구간은 summary 함수를 이용한다
<해석>
obs군과 비교하였을 때 Lev 군의 HR 95% 신뢰구간은 0.842-1.13
obs군과 비교하였을 때 Lev+5FU 군의 HR 95% 신뢰구간은 0.546-0.759
Lev+5FU의 HR 95% 신뢰구간이 1을 포함하지 않으므로 obs에 비해 유의한 차이가 있다
2. Forest plot 그리기 : ’ggforest’함수 사용
ggforest(result1,data=colon)
’coxph’함수를 통해 추정된 HR와 그 신뢰구간을 그림으로 표현할 수 있다
‘survminer’ 패키지에 내장된 ’ggforest’함수를 사용한다
3. Cox 비례 위험 가정의 확인
Cox 비례 위험 모형을 적합 하기 위해서는 어떤 시점에서라도 해당변수의 HR이 일정하다는 전제조건이 필요하다
즉, 연구 1년째 이든 5년째 이든, 해당 변수(ex. 치료법)에 따른 HR가 일정해야한다
이것을 ’Cox Propotional Hazard Assumption’이라고 하는데 ,이를 검증하는 3가지 방법을 간략히 소개한다
1) Goodness of fit test : ’cox.zph’함수 사용
cox.zph(result1,global = F)## chisq df p
## rx 0.648 2 0.72<해석>
‘cox.zph’ 에서 귀무가설은 해당 변수가 ’비례위험을 만족한다’이다
즉, p-value가 0.05 이상이면 ’비례위험 모형을 만족한다’라고 할 수 있으며, 위 결과에서 ’rx(치료법)’변수는 비례 위험 모형을 만족한다.
장점: p-value라는 객관적인 값으로 비례 위험 가정을 만족하는지 평가할 수 있다
단점: sample size가 작으면, p-value값이 작게 나오는것이 부정확할 수 있다
2) log-log plot
result2<-survfit(Surv(time,status)~rx,data=colon) # survfit 객체를 생성
plot (result2, fun='cloglog', # x,y축 값을 log로 변경하여 그래프 그림
lty = 1:3, # 그룹 마다 line type 다르게
col = 1:3) # 그룹 마다 line color 다르게
legend('bottomright',lty=1:3,col=1:3,legend = levels(colon$rx)) # 범례 작성 
log-log plot을 그리고, 그룹간에 평행하게 나오면 비례 위험 가정을 만족한다.
장점: 직관적
단점: ’평행’하다라는 객관적 기준이 없음
3) observed vs. expected plots
Kaplan-Meier 생존곡선과, Cox 비례위험 모형에서 나온 생존곡선의 차이를 본다.
result2<-survfit(Surv(time,status)~rx,data=colon) # survfit 객체를 생성
#observed
plot(result2, lty='dashed', col= 1:3, main="Observed Versus Expected Plots by Treatment", ylab="Survival probability", xlab="Time") # KM curve 생성
par(new=T)
#expected
result1<-coxph(Surv(time,status)~rx,data=colon) # coxph 객체를 생성
new_df <- data.frame(rx=c('Obs','Lev','Lev+5FU'))
expect <- survfit(result1, newdata=new_df)
plot(expect, lty="solid", col=1:3)
Kaplan Meier 생존곡선과, Cox 비례위험 모형에서 나온 생존곡선이 거의 일치한다
따라서, 비례 위험 가정을 만족한다
장점: 직관적
단점: 두 곡선이 겹친다라는 객관적 기준이 없음
4. Cox 비례위험 모형 만들기
Cox 비례 위험 모형을 만든다는 것은 수집한 여러 독립변수 중에 생존에 유의한 영향을 미치는 변수를 찾고,
그 변수들을 조합하여 최소한의 변수로 최대의 설명력을 가지는 모형을 찾는것이다.
1) 단변량 분석 : ’moonBook’패키지의 ’mycph’함수 활용
보통 Cox 비례위험 모형을 만들때는 단변량 분석으로 가능성 없는 변수 (보통 p>0.2)를 걸러내는 작업을 한다.
그러나 변수가 많다면, 그 변수의 갯수만큼 coxph 함수를 돌려야하는 번거로움이 있다
‘moonBook’ 패키지의 ‘mycph’함수를 활용하면 이런 번거로움을 해결할 수 있다
(출처:‘의학논문 작성을 위한 R통계와 그래프’ 성빈센트병원 순환기내과 문건웅 교수님)
attach(colon) #attach:데이터의 변수를 자유롭게 불러주는 기능
colon$TS<-Surv(time,status)
생존분석을 위한 Surv()함수의 결과를 colon$TS에 저장한다.
library(moonBook) # moonBook 패키지를 불러온다
result3<-mycph(TS~.-id-study-time-status-etype,data=colon)##
## mycph : perform coxph of individual expecting variables
##
## Call: TS ~ . - id - study - time - status - etype, data= colonresult3## 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.000‘moonBook’ 패키지의 ’mycph’함수로 단변량 분석을 한꺼번에 시행한다.
‘~.’은 모든 변수를 포함시키라는 명령어이고, 그 다음’-’는 해당변수를 제외하라는 명령어이다.
’id’는 환자 일련번호, ’study’는 모두1로 의미없는 변수이며, time, staus는 이미 사용한 변수이므로 제외하여야 한다.
’etype’은 event의 종류이니 제거하고, ’nodes’는 전이된 림프절의 수이고 ’node4’는 림프절이 4개 이상이나 이하이냐 하는 데이터이므로 둘 중하나는 제외해야여야 한다.
결과에서 p-value가 0.2보다 큰 sex, age를 제하고 나머지 변수들로 다변량 분석을 시행한다
2) 다변량 분석: ‘step’ function으로 최적의 모형을 적합한다
colon1<-na.omit(colon) # colon 데이터에서 결측치 제거
result4<-coxph(Surv(time,status)
~rx+obstruct+perfor+adhere+node4+differ+extent+surg,
data=colon1) # p value>0.2 변수를 제외하고 모든 변수를 포함시킴
finalmodel<-step(result4, direction = 'backward') # 후진제거법 적용 ## Start: AIC=12267.85
## Surv(time, status) ~ rx + obstruct + perfor + adhere + node4 +
## differ + extent + surg
##
## Df AIC
## - perfor 1 12266
## <none> 12268
## - adhere 1 12270
## - differ 1 12271
## - obstruct 1 12272
## - surg 1 12275
## - rx 2 12294
## - extent 1 12299
## - node4 1 12403
##
## Step: AIC=12266.18
## Surv(time, status) ~ rx + obstruct + adhere + node4 + differ +
## extent + surg
##
## Df AIC
## <none> 12266
## - adhere 1 12268
## - differ 1 12269
## - obstruct 1 12271
## - surg 1 12274
## - rx 2 12292
## - extent 1 12298
## - node4 1 12402다변량분석에서 모형적합에는 ’step’함수를 사용한다
’step’함수를 사용하기위해서는 ’na.omit(데이터)’으로 결측치 제거 과정이 필요하다.
‘step’함수 적용방식에는 ’forward’, ‘backward’, ‘both’ 3가지가 있다
모든 변수를 포함시킨 상태에서 유의하지 않은 변수를 차례로 제거해나가는 ’backward elimination’을 선택한다.
최종 모형에 선택된 변수는 ‘rx’, ‘obstruct’, ‘adhere’, ‘node4’,‘differ’, ‘extent’, ’surg’이다.
summary(finalmodel)## Call:
## coxph(formula = Surv(time, status) ~ rx + obstruct + adhere +
## node4 + differ + extent + surg, data = colon1)
##
## n= 1776, number of events= 876
##
## coef exp(coef) se(coef) z Pr(>|z|)
## rxLev -0.01632 0.98381 0.07929 -0.206 0.83693
## rxLev+5FU -0.41577 0.65983 0.08601 -4.834 1.34e-06 ***
## obstruct 0.21557 1.24057 0.08345 2.583 0.00979 **
## adhere 0.18892 1.20794 0.09069 2.083 0.03725 *
## node4 0.86615 2.37775 0.07071 12.250 < 2e-16 ***
## differ 0.15551 1.16826 0.06959 2.235 0.02543 *
## extent 0.47124 1.60198 0.08399 5.611 2.01e-08 ***
## surg 0.23131 1.26025 0.07408 3.122 0.00179 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## rxLev 0.9838 1.0165 0.8422 1.149
## rxLev+5FU 0.6598 1.5155 0.5575 0.781
## obstruct 1.2406 0.8061 1.0534 1.461
## adhere 1.2079 0.8279 1.0112 1.443
## node4 2.3777 0.4206 2.0700 2.731
## differ 1.1683 0.8560 1.0193 1.339
## extent 1.6020 0.6242 1.3588 1.889
## surg 1.2602 0.7935 1.0899 1.457
##
## Concordance= 0.659 (se = 0.009 )
## Likelihood ratio test= 250.1 on 8 df, p=<2e-16
## Wald test = 253.5 on 8 df, p=<2e-16
## Score (logrank) test = 266.1 on 8 df, p=<2e-16최종 모형에 선택된 변수 HR의 95% 신뢰구간은 ’summary’함수를 이용한다.
3) 최종 적합 모형을 Forest plot으로 표현: ’ggforest’함수 사용
ggforest(finalmodel,data=colon1)