[R] Survival Analysis (Kaplan-Meier, Log-rank, Cox PH)

이전 포스팅에서 생존 분석의 개념, Kaplan-Meier 추정과 Log-rank test 그리고 Cox 비례위험모형까지 알아보았다. 이번 포스팅에서는 R에서 생존 분석을 수행하는 방법에 대해 알아보겠다.

 

분석 절차는 아래와 같다.

1. Kaplan-Meier 생존 곡선 

2. Log-rank Test

3. Cox 비례 위험 모형

 

 

1. Kaplan-Meier 생존 곡선 

먼저 'survival' 패키지의 survfit() 함수를 사용하여 카플란마이어 생존 곡선을 추정한다.

survfit의 종속변수는 Surv(시간, 발생여부) 형식으로 넣어주고, 독립변수로는 group을 넣어 주어 각 그룹의 시간에 따른 발생 확률의 변화를 알아본다.

  surv_obj <- Surv(time = df.surv$et, event = df.surv$event)
  fit <- survfit(surv_obj ~ group, data = df.surv)
  summary(fit)

 

'survminer' 패키지의 ggsurvplot() 함수를 사용하여 그래프를 그릴 수 있다.

ggsurvplot(
  fit,data = df.surv,size = 1,                
  palette = c("#2E9FDF","#E7B800"),
  censor = F, conf.int = TRUE, pval = F, risk.table = TRUE,  
  break.x.by = 2, xlab = "Years after Surgery", xlim = c(0,10),
  legend.labs = c('Achieved',"Failed"), risk.table.height = 0.25, 
  ggtheme = theme_bw(),            
  font.x = c(16),  # x축 제목 크기
  font.y = c(16),  # y축 제목 크기
  font.tickslab = c(14),  # 축 텍스트(눈금) 크기
  font.legend = c(12),  # 범례 텍스트 크기
  font.title = c(14)  # 범례 제목 크기     
)

 

2. Log-rank Test

 

두 그룹의 생존 곡선에 차이가 있는지 확인하기 위해 'survival' 패키지의 survdiff() 함수를 사용하여 Log-rank test를 수행한다.

위에서 ggsurvplot()으로 그래프를 그릴 때 pval=T로 하여 p-value를 확인할 수도 있다.

survdiff(surv_obj ~ group, data = df.surv)

> survdiff(surv_obj ~ group, data = df.surv)
Call:
survdiff(formula = surv_obj ~ group, data = df.surv)

                  N Observed Expected (O-E)^2/E (O-E)^2/V
df.surv$group=0 186       23     29.6      1.47       4.1
df.surv$group=1 107       24     17.4      2.50       4.1

 Chisq= 4.1  on 1 degrees of freedom, p= 0.04

 

 

3. Cox 비례위험모형

공변량을 포함하여 여러 설명 변수들이 생존 시간에 미치는 영향을 평가하기 위해 'survival' 패키지의 coxph() 함수를 사용하여 Cox 비례위험모형을 만든다. 

 # Cox PH 모형
fit.coxph <- coxph(surv_obj ~ group + age + sex, data = df.surv)
fit.coxph

> fit.coxph
Call:
coxph(formula = as.formula(surv_obj ~ group + age + sex), data = df.OS)

             coef exp(coef)  se(coef)      z       p
Group1    1.601895  4.962430  0.607102  2.639 0.00833
age1      1.635166  5.130310  0.779163  2.099 0.03585
sex1     -0.010500  0.989555  0.004639 -2.263 0.02362

Likelihood ratio test=21.36  on 4 df, p=0.0002684
n= 293, number of events= 14

 

Cox PH 모형에서는 비례위험가정이 중요하기 때문에 cox.zph() 함수로 가정을 확인하고, 가정을 만족하지 않는 공변량은 모형에서 제외한다.

# 비례위험가정 확인
cox.zph(fit.coxph)

>   cox.zph(fit.coxph, global = T) 
         chisq df    p
group   0.0152  1 0.90
age      0.5719  1 0.45
sex     1.6546  1 0.20
GLOBAL  2.2602  4 0.69

 

step() 함수로 stepwise 변수선택을 진행하여 최종 모델을 만든다.

# stepwise 변수 선택
step.coxph <- step(fit.coxph, direction = "backward")
summary(step.coxph)