[R] 생존분석(Survival analysis)

[R] 생존분석(Survival analysis) Start

BioinformaticsAndMe

#생존분석에 대한 개념적 설명은 아래 포스팅 참조

https://bioinformaticsandme.tistory.com/223

1. 생존분석 라이브러리 로딩

---

: 생존분석에 필요한 패키지(survival, survminer) 및 데이터프레임 관리 패키지(dplyr) 로딩

: 참고로, 패키지 설치 에러에는 여러 이유가 있겠으나, 일반적으로 R 최신 버전 상태에서 패키지를 설치하면 해결됨

*아래 튜토리얼은 'R version 3.6.1 (2019-07-05)'에서 수행됨

# 패키지 인스톨 후, 로딩
install.packages(c('survival', 'survminer', 'dplyr'))
library(survival)
library(survminer)
library(dplyr)

2. 예제 데이터셋 로딩

---

: 난소암(Ovarian cancer) 데이터를 생존분석의 예제로 사용

# glimpse 함수로 난소암 데이터의 내용/구조/타입을 파악
data(ovarian)
glimpse(ovarian)
Observations: 26
Variables: 6
$ futime   <dbl> 59, 115, 156, 421, 431, 448, 464, 475, 477, 563, 638,...  #생존 시간
$ fustat   <dbl> 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0,...  #Censored data 여부
$ resid.ds <dbl> 2, 2, 2, 2, 2, 1, 2, 2, 2, 1, 1, 1, 2, 2, 1, 1, 2, 1,...  #종양 퇴행
$ age      <dbl> 72.3315, 74.4932, 66.4658, 53.3644, 50.3397, 56.4301,...  #환자 나이
$ rx       <dbl> 1, 1, 1, 2, 1, 1, 2, 2, 1, 2, 1, 2, 2, 2, 1, 1, 1, 1,...  #치료 요법 
$ ecog.ps  <dbl> 1, 1, 2, 1, 1, 2, 2, 2, 1, 2, 2, 1, 2, 1, 1, 2, 2, 1,...  #환자 전신 활성도(ECOG)

# 데이터 Labeling
ovarian$rx <- factor(ovarian$rx, 
                     levels = c("1", "2"), 
                     labels = c("A", "B"))
ovarian$resid.ds <- factor(ovarian$resid.ds, 
                           levels = c("1", "2"), 
                           labels = c("no", "yes"))
ovarian$ecog.ps <- factor(ovarian$ecog.ps, 
                          levels = c("1", "2"), 
                          labels = c("good", "bad"))

# 히스토그램으로 나이 데이터 분포 확인
hist(ovarian$age) 

# 위 이항분포 히스토그램 결과에 근거하여, '50 years' 기준의 두 그룹으로 나눔
ovarian <- ovarian %>% mutate(age_group = ifelse(age >=50, "old", "young"))
ovarian$age_group <- factor(ovarian$age_group)

3. Kaplan-Meier(카플란-마이어) 생존분석

---

# Surv 함수로 Survival object 생성 (Survival object는 survfit 함수에 인식되도록 특정칼럼을 컴파일한 형태)
# surv_object 결과의 '+' 표시는 Censored data를 의미
surv_object <- Surv(time = ovarian$futime, event = ovarian$fustat)
surv_object 
[1]   59   115   156   421+  431   448+  464   475   477+  563   638 
[12]  744+  769+  770+  803+  855+ 1040+ 1106+ 1129+ 1206+ 1227+  268 
[23]  329   353   365   377+

# survfit 함수를 통해, Kaplan-Meier 곡선 Fitting
fit1 <- survfit(surv_object ~ rx, data = ovarian)
summary(fit1)
Call: survfit(formula = surv_object ~ rx, data = ovarian)

               rx=A 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
  59     13       1    0.923  0.0739        0.789        1.000
  115     12       1    0.846  0.1001        0.671        1.000
  156     11       1    0.769  0.1169        0.571        1.000
  268     10       1    0.692  0.1280        0.482        0.995
  329      9       1    0.615  0.1349        0.400        0.946
  431      8       1    0.538  0.1383        0.326        0.891
  638      5       1    0.431  0.1467        0.221        0.840

               rx=B 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
  353     13       1    0.923  0.0739        0.789        1.000
  365     12       1    0.846  0.1001        0.671        1.000
  464      9       1    0.752  0.1256        0.542        1.000
  475      8       1    0.658  0.1407        0.433        1.000
  563      7       1    0.564  0.1488        0.336        0.946

# ggsurvplot 함수로 생존분석 결과 시각화 (추가로 p value 표시)
# 그래프 내에 있는 수직바는 Censored data를 의미함
ggsurvplot(fit1, data = ovarian, pval = TRUE)

# resid.ds 기준으로 나뉜 두 그룹의 생존분석은 p-value가 약 0.05로 유의한 결과를 냄
fit2 <- survfit(surv_object ~ resid.ds, data = ovarian)
ggsurvplot(fit2, data = ovarian, pval = TRUE)

4. Cox proportional hazards model (콕스 비례위험모델)

---

: Hazard Ratio(HR)는 상대적 위험도이므로, 그 위험비는 항상 reference에 대한 특정값임

: HR>1(사망 위험 증가), HR<1(사망 위험 감소)

: 아래 그림에서 0.25는 Treatment B를 받은 환자들이 Treatment A를 받은 환자들에 비해 사망률이 0.25배 감소했다는 의미

*Treatment A는 reference로 사용됨

# coxph 및 ggforest 함수를 통해, Hazard ratio 결과 출력
fit.coxph <- coxph(surv_object ~ rx + resid.ds + age_group + ecog.ps, 
                   data = ovarian)
ggforest(fit.coxph, data = ovarian)

#생존분석에 대한 개념적 설명은 아래 포스팅 참조

https://bioinformaticsandme.tistory.com/223

#Reference

1) http://www.gums.ac.ir/Upload/Modules/Contents/asset68/Medical%20Statistic%20Made%20Easy.pdf

2) https://www.datacamp.com/community/tutorials/survival-analysis-R#fourth

3) https://ko.wikipedia.org/wiki/%EC%83%9D%EC%A1%B4%EB%B6%84%EC%84%9D

4) https://namu.wiki/w/%EC%83%9D%EC%A1%B4%20%EB%B6%84%EC%84%9D

5) https://www.partek.com/webinar/survival-analysis-with-partek-genomics-suite-software/

6) https://ko.wikipedia.org/wiki/%EC%A4%91%EB%8F%84%EC%A0%88%EB%8B%A8

7) https://ko.wikipedia.org/wiki/%EB%A1%9C%EA%B7%B8%EC%88%9C%EC%9C%84%EB%B2%95

8) https://www.youtube.com/watch?v=czQ3l0QXxnA

9) https://rexsoft.org/?page_id=485

10) http://www.e-urol-sci.com/viewimage.asp?img=UrolSci_2018_29_5_223_240363_t6.jpg

[R] 생존분석(Survival analysis) End

BioinformaticsAndMe