데이터마이너를 꿈꾸며

6.3 비모수적 방법을 이용한 생존함수의 비교

[예] 흑색종(melanoma) 환자들에 대한 BCG와 CP(coryne-bacterium parvum)의 생존지속 효과를 비교하기 위한 연구에서 30명의 흑색종 환자 중 11명은 BCG 처리를 받고 나머지 19명은 CP처리를 받았다고 한다. 중도절단이 포함되어 있는 경우에 이 두 그룹의 생존분포를 비교하기 위한 방법을 알아보자.

1. 데이터 입력

> library(survival)

> setwd('c:/Rwork')

> melanoma = read.table('melanoma.txt',header=T)

> head(melanoma,3)

  癤퓍ime status   x

1    33.7      0 BCG

2     3.9      1 BCG

3    10.5      1 BCG

> colnames(melanoma)<-c('time','status','x')

> head(melanoma)

  time status   x

1 33.7      0 BCG

2  3.9      1 BCG

3 10.5      1 BCG

4  5.4      1 BCG

5 19.5      1 BCG

6 23.8      0 BCG

> attach(melanoma)

2. 누적한계추정치(Kaplan-Meier 추정치)

> fit2 = survfit(Surv(time,status)~x,data=melanoma)

> summary(fit2)

Call: survfit(formula = Surv(time, status) ~ x, data = melanoma)

                x=BCG

 time n.risk n.event survival std.err lower 95% CI upper 95% CI

  3.9     11       1    0.909  0.0867        0.754        1.000

  5.4     10       1    0.818  0.1163        0.619        1.000

  7.9      9       1    0.727  0.1343        0.506        1.000

 10.5      8       1    0.636  0.1450        0.407        0.995

 19.5      4       1    0.477  0.1755        0.232        0.981

                x=CP

 time n.risk n.event survival std.err lower 95% CI upper 95% CI

  6.9     19       1    0.947  0.0512        0.852        1.000

  7.7     18       1    0.895  0.0704        0.767        1.000

  8.0     16       1    0.839  0.0854        0.687        1.000

  8.3     13       1    0.774  0.1003        0.601        0.998

 24.4      3       1    0.516  0.2211        0.223        1.000

> fit2

Call: survfit(formula = Surv(time, status) ~ x, data = melanoma)

       n events median 0.95LCL 0.95UCL

x=BCG 11      5   19.5    10.5      NA

x=CP  19      5     NA    24.4      NA

3. 사망시점의 사분위수 추정치와 그의 신뢰구간

> quantile(fit2,probs=c(0.25,0.5,0.75),conf.int=T)

$quantile

        25   50 75

x=BCG  7.9 19.5 NA

x=CP  24.4   NA NA

$lower

       25   50   75

x=BCG 5.4 10.5 19.5

x=CP  8.0 24.4 24.4

$upper

      25 50 75

x=BCG NA NA NA

x=CP  NA NA NA

4. 데이터 시각화

> plot(fit2,xlab='time',ylab='survival function',lty=c(1,2),col=c(1,2))

> legend(5,0.2,c('cp 처리 그룹','BCG 처리 그룹'),lty=c(2,1),col=c(2,1))

> abline(h=0.5)

> abline(v=c(10.5,24.4))

5. 로그 순위 검정법(log-rank test)과 Gehan-Wilcoxon 검정법 비교

1) 로그 순위 검정법(log-rank test)

> survdiff(Surv(time,status)~x,data=melanoma)

Call:

survdiff(formula = Surv(time, status) ~ x, data = melanoma)

       N Observed Expected (O-E)^2/E (O-E)^2/V

x=BCG 11        5     3.68     0.469     0.747

x=CP  19        5     6.32     0.274     0.747

 Chisq= 0.7  on 1 degrees of freedom, p= 0.387

2) Gehan-Wilcoxon 검정법

> survdiff(Surv(time,status)~x,rho=1,data=melanoma)

Call:

survdiff(formula = Surv(time, status) ~ x, data = melanoma, rho = 1)

       N Observed Expected (O-E)^2/E (O-E)^2/V

x=BCG 11     4.31     3.07     0.500     0.929

x=CP  19     4.10     5.34     0.288     0.929

 Chisq= 0.9  on 1 degrees of freedom, p= 0.335  

두 검정법 모두 p value가 0.05보다 크기 떄문에 유의하지 않다. 즉 BCG, CP 두 그룹 간의 생존함수는 유의한 차이를 보이지 않는다.

이상 그래프에서 보듯이 두 그룹의 누적한계추정치의 그래프도 교차되지 않고 나란한 형태를 보이므로 로그-순위 검정법이 타당한 것이었음을 알 수 있다.

출처: 보건정보데이터 분석 (이태림, 이재원, 김주한, 장대흥 공저) 

비모수적 방법 non-parametric method

6.2.2 누적한계추정법 product-limit method

생명표 방법: 기간을 1년, 6개월 등 특정 단위로 나눠 구분

누적한계추정법: 매 사건이 발생할 때 마다 해당 시점을 표기

누적한계추정법(product-limit method)은 생존함수를 추정하는 대표적인 방법 중 하나로 연구자의 이름을 따서 Kaplan-Meier추정법이라고도 한다. 모든 환자의 생존시간 또는 중도절단 시간이 각각 관찰되었다고 하자. 모든 자료의 생존시간 도는 중도절단 x1,…,투을 순서대로 배열한 것을 t1<t2<…<tn이라 하고, δi = 0, 그렇지 않은 경우 δi =1로 정의한다. 즉 중도 절단 되지 않고 사건이 발생한 경우 status = 1, 중도절단된 경우 stats = 0 로 표기한다.

예] 신장이식수술을 받은 15명 환자들의 호전기간(remission duration)이 아래와 같다고 하자. 여기서 +로 표기된 것은 중도절단된 자료를 가르킨다. 각 시점에서 생존함수를 누적한계추정법을 이용하여 추정해보자

표. 신장이식 환자들의 호전기간 (단위 : 일)

1. 데이터 불러오기

> library(survival)

> setwd('c:/Rwork')

> kidney<-read.table('kidney.txt',header=T)

> kidney #status =1 사건, 0 = 절단, 절단 표시 매우 중요

time status

1   3.0      1

2   4.0      0

3   4.5      1

4   4.5      1

5   5.5      1

6   6.0      1

7   6.4      1

8   6.5      1

9   7.0      1

10  7.5      1

11  8.4      0

12 10.0      0

13 10.0      1

14 12.0      1

15 15.0      1

> attach(kidney)

2. 누적한계추정치(Kaplan-Meier 추정치)

> fit1 = survfit(Surv(time,status)~1, data=kidney) #~ 우측에는 공변량

> summary(fit1)

Call: survfit(formula = Surv(time, status) ~ 1, data = kidney)

time n.risk n.event survival std.err lower 95% CI upper 95% CI

3.0     15       1    0.933  0.0644       0.8153        1.000

4.5     13       2    0.790  0.1081       0.6039        1.000

5.5     11       1    0.718  0.1198       0.5177        0.996

6.0     10       1    0.646  0.1275       0.4389        0.951

6.4      9       1    0.574  0.1320       0.3660        0.901

6.5      8       1    0.503  0.1336       0.2984        0.846

7.0      7       1    0.431  0.1324       0.2358        0.787

7.5      6       1    0.359  0.1283       0.1781        0.723

10.0      4       1    0.269  0.1237       0.1094        0.663

12.0      2       1    0.135  0.1135       0.0258        0.703

15.0      1       1    0.000     NaN           NA           NA

> fit1

Call: survfit(formula = Surv(time, status) ~ 1, data = kidney)

n  events  median 0.95LCL 0.95UCL

15      12       7       6      NA

Error: invalid multibyte character in parser at line 1

-> fit1 = survfit(Surv(time,status)~1, data=kidney) #~ 우측에는 공변량

6.5일 에서의 생존확률은 50.3%, 7.0일 이 지나면 생존율이 43.1%로 50%이하가 된다.

3. 사망시점의 사분위수 추정치와 그의 신뢰구간, 가령 사망자가 50%되는 시점은 언제인가?

> quantile(fit1,probs=c(0.25,0.5,0.75),conf.int=T) #신뢰구간까지

$quantile

25   50   75

5.5  7.0 12.0

$lower

25  50  75

4.5 6.0 7.0

$upper

25  50  75

7.5  NA  NA

-> 생존확률이 75%인 경우, 즉 사망확률이 25%인 경우의 시점은 5.5일, 생존확률이 50%인 경우 7.0일, 생존확율이 25%인 경우의 시점은 12.0 일

4. 누적한계추정치의 95% 신뢰구간 그래프

> plot(fit1,xlab='time',ylab='Survival function',lwd=2)

> legend(0.2,0.3,c('KM estimate','95% CI'),lty=c(1,2))

-> 점선은 신뢰구간을 의미, 실선은 생존함수추정치.생존확율이 50%인 경우는 약 7.0일임을 알 수 있다.

출처: 보건정보 데이터분석 (이태림, 이재원, 김주한, 장대흥 공저)