카플란 메이어 (kaplan meier) 생존분석 - R

 

/*

 * http://sosal.kr/
 * made by so_Sal
 */

 

카플란메이어 생존 분포는 샘플에서 사망이 발생함에 따라 누적되는 형태의 계단형태 생존곡선을 나타낸다.

생존분석에서의 반응변수 (종속변수)는 결국 해당 시간에 생존률이 얼마가 되는지를 나타낸다.

 

 

1. Kaplan meier 기본 함수 및 사용방법

- survival 패키지 다운로드

- 환자의 생존상태 식별하기 (surv 함수)

- 시간에 따른 생존커브 구하기 (survfit)

- Cumulative hazard 그래프 그리기

- 다중 생존곡선 그리기

 

 

2. 공개데이터 colon을 이용한 kaplan meier 생존분석

- colon 전체데이터의 생존곡선 구하기

- 성별에 따른 다중 생존곡선 구하기

- rx 치료법에 따른 다중 생존곡선 구하기

- 나이에 따른 다중 생존곡선 구하기

 

 

 

1. Kaplan meier 기본 함수 및 사용방법

 

 

- survival 패키지 다운로드

> install.packages("survival")

> library(survival)

 

 

 

간단한 예로, 다음의 값들을 특정 질병에 걸린 8명의 환자들의 데이터라고 가정해보자

 

- Column 설명

Patient: 환자 목록

Days: 가장 마지막으로 연락된 기간, 사망한 경우 죽은 날을 나타낸다.

Status: 죽은 환자는 1, 살아있는 환자는 0

 

환자	생존기간	생존여부
환자1	3	1
환자2	5	0
환자3	7	1
환자4	11	0
환자5	12	1
환자6	13	1
환자7	17	1
환자8	18	1

 

 

> death_time = c(3,5,7,11,12,13,17,18)
> status=c(1,0,1, 0, 1, 1, 1, 1)

 

 

 

 

- 환자의 생존상태 식별하기

 

따라서 환자2, 환자4를 제외한 나머지 6명의 환자는 사망한 상태이다.

이것은 Surv 라는 함수로 표현이 가능하다.

 

> Surv(death_time, status)
[1]  3   5+  7  11+ 12  13  17  18

 

status 벡터 덕분에 환자2와 환자4는 +로 표시되어 생존해 있다는 것을 알 수 있다.

 

 

 

- 시간에 따른 생존커브 구하기 (coxph, kaplan meier 등)

 

> fit = survfit(Surv(death_time, status)~1)
> summary(fit)
Call: survfit(formula = Surv(death_time, status) ~ 1)

 

 time n.risk n.event   survival    std.err   lower 95%CI   upper 95%CI
     3      8           1       0.875       0.117          0.6734                    1
     7      6           1       0.729       0.165          0.4680                    1
   12     4          1       0.547       0.201          0.2665                    1
   13     3          1      0.365       0.200           0.1244                   1
   17     2         1       0.182       0.163          0.0315                    1
   18     1         1       0.000         NaN               NA                  NA

 

input으로 들어간 환자는 8명이지만, 환자2와 환자4의 경우에는 마지막 컨택시간만 있을 뿐, 그들의 생존상태는 알 수 없다.

따라서 그에따라 survival rate와 std.err 값도 보정된다.

 

survfit에서 '~1' 이라고 들어간 것은 사실 그룹을 나누는 인자다.

예를 들어, 백혈병 환자의 경우 ALL 백혈병과 AML 백혈병으로 나눌 수 있는데, 생존분석을 통해 두 백혈병을 비교하고자 한다면 ~group 으로 두가지의 survival curve를 한번에 그릴 수 있다.

 

 

survfit의 결과인 fit 변수를 plot으로 그린다면 Kaplan meier 곡선을 볼 수 있다.

> plot( fit, ylab="Survival", xlab="Days", mark.time=T)

 

 

 

 

점선은 해당 생존곡선의 신뢰 구간을 보여준다.

Lower 95% CI,  Upper 95% CI는 신뢰구간을 뜻하며 CI는 confidence interval를 뜻한다.

생존곡선 위 Days -5, 12에 그어진 짧은 선은 사망이 확인되진 않았지만, 마지막으로 연락된 환자2, 환자4를 뜻한다.

이를 없애고자 하면 plot에 conf.int=FALSE 라는 인자를 넣어주면 된다.

ex) plot( fit, ylab="Survival", xlab="Days", conf.int=FALSE)

 

 

 

- Cumulative hazard 그래프 그리기

survfit 함수로부터 리턴받은 fit을 그대로 사용하여 Cumulative hazard 그래프를 그려보자.

cumulative hazard

> plot( fit, ylab="Cumulative hazard H(ti)", xlab="Days", fun="cumhaz", mark.time=T)

 

Cumulative hazard plotting에 대한 장단점, 혹은 내용은 아래에서 확인

http://www.itl.nist.gov/div898/handbook/apr/section2/apr222.htm

 

 

 

- 다중 생존곡선 그리기

지금은 survfit 함수에서 '~1'로 하나의 환자군의 생존곡선 그래프를 그렸지만, '~' 틸데를 이용하여 그룹을 나누어 survfit 함수를 사용한다면, 다중 생존곡선을 그릴 수 있다.

 

1~50일 사이에 100명의 환자를 무작위로 나타내고 생존상태도 무작위로 넣을것이다.

100명의 환자중 50명으로 두 그룹을 나누어 survfit 함수를 사용해보자. 

 

> death_time = sample(1:50, 100, replace=TRUE)
> status = sample(1:0, 100, replace=TRUE)
> group = c(rep(0,50), rep(1, 50))

> group <- factor(group)                        # 0, 1 2개의 level로 factor화
> fit = survfit(Surv(death_time, status)~group)

 

Call: survfit(formula = Surv(death_time, status) ~ group)

 

                group=0
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    4     48       1    0.979  0.0206        0.940        1.000
    6     47       3    0.917  0.0399        0.842        0.998
    7     43       2    0.874  0.0481        0.785        0.974
                  - 중략 -
   46      9       1    0.400  0.0905        0.256        0.623
   48      7       1    0.343  0.0938        0.200        0.586
   49      4       1    0.257  0.1023        0.118        0.560

 

                group=1
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    5     49       1    0.980  0.0202        0.941        1.000
    6     47       1    0.959  0.0286        0.904        1.000
    7     46       1    0.938  0.0347        0.872        1.000

                       - 중략 -
   40      9       1    0.468  0.0990        0.309        0.708
   41      8       1    0.409  0.1025        0.251        0.668
   46      4       1    0.307  0.1173        0.145        0.649

 

각 그룹 0, 1에 따라서 서로 다른 생존곡선을 만든 것을 볼 수 있다.

직접 plot으로 그래프를 그려 확인해보자. 

 

> plot( fit, ylab="Survival", xlab="Days", col=c("red","blue"), lty=1:2, mark.time=T)

> legend("topright", legend=levels(group), col=c("red","blue"), lty=1:2)
  

 

그룹1과 그룹2의 명확한 구분이 없고 일양분포의 데이터이기 때문에 별 다른 차이를 볼 수 없습니다.

 

 

 

 

 

 

 

 

 

 

 

 

 

2. 공개데이터 colon을 이용한 kaplan meier 생존분석

 

데이터를 로드하고 컬럼변수를 확인해보자.

colon 데이터는 survival 패키지가 로드되어 있으면 바로 불러올 수 있다.

 

> data(colon)

> attach(colon) # colon 데이터의 변수들을 바로 사용할 수 있도록 해준다.

> str(colon)
'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 ...

 

 

https://stat.ethz.ch/R-manual/R-devel/library/survival/html/colon.html

위의 사이트에서 각 colomn의 자세한 정보를 알 수 있다.

 

colon cancer는 대장에 발생하는 악성 종양으로, 대장암 (결장직장암)이라고 불린다.

 

간단하게 데이터를 설명하면

rx:치료법으로 Lev+5FU, Obs, Lev 3가지 분류로 나눠진다. factors 라고 생각하면 되겠다.

sex: 0=female, 1=male

age: 나이

status: 1=death, 0은 생존시간이 더 길 수 있음을 말함.

time: 생존 기간, 혹은 마지막으로 연락이 닿은 기간.

정도가 되겠다.

 

 

 

- colon 전체데이터의 생존곡선 구하기

 

> fit = survfit(Surv(time, status)~1, data=colon)

> plot(fit, ylab="Survival", xlab="Days", mark.time=T)

 

 

 

- 성별에 따른 다중 생존곡선 구하기

> fit = survfit(Surv(time, status==1)~sex, data=colon)

> plot(fit, ylab="Survival", xlab="Days", col=c("red","blue"), lty=1:2, mark.time=T)

> legend("topright", legend=c("female", "male"), col=c("red","blue"), lty=1:2)

 

 

 

colon이란 질병은 성별에 따라 생존률이 다르진 않는 것 같다.

물론 남여에 따른 생존률이 다른가? 에 대한 것은 t.test나 다른 regression 방법으로 p-value를 직접 구하는게 맞다.

 

 

- rx 치료법에 따른 다중 생존곡선 구하기

> fit = survfit(Surv(time, status==1)~rx, data=colon)
> plot(fit, ylab="Survival", xlab="Days", col=1:3, lty=1:3, mark.time=T)

> legend("topright", legend=levels(rx), col=1:3, lty=1:3)

rx는 이미 factor화 된 데이터이기 때문에 levels()라는 함수로 legend에 넣으면 된다.

 

 

Lev+5FU 치료법으로 치료한 환자의 예후가 좋다는 것을 볼 수 있다.

확실히 Lev, Obs 치료법보다 좋다는걸 증명하기 위해선 통계적인 검정을 해야겠지만

일단 이 그래프를 통해 어느정도 도움이 될 것 같다고 볼 수 있을 것이다.

 

 

- 나이에 따른 다중 생존곡선 구하기

 

나이가 50세 이상, 혹은 미만으로 두 그룹을 나누어 생존곡선을 그려보자,

> age_group <- ifelse(age<50,0,1)

> fit = survfit(Surv(time, status==1)~age_group, data=colon)
> plot(fit, ylab="Survival", xlab="Days", col=1:2, lty=1:2, mark.time=T)

> legend("topright", legend=c("age<50","age>=50"), col=1:2, lty=1:2)

 

오히려 나이가 젊은사람이 약간 더 빨리 사망하는것을 볼 수 있습니다.

 

 

 

+

서로다른 그룹의 생존이 얼마나 다른지 통계적으로 확인하기 위해 LogRank test를 사용할 수 있습니다.

혹은 coxph survival을 통해서 어떤 변수가 생존에 영향을 주는지 알아볼 수 있습니다.