데이터마이너를 꿈꾸며

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.1생명표 방법 life table method

생명표 방법 개념 추가 설명

예] 협심증(angina pectoris)이 있는 2,418명의 남성들에 대한 생존자료이다. 이 생명표에서 경과기관에 따른 생존확률을 구해보자

1. 데이터 입력

> library(KMsurv)

> setwd('C:/Rwork')

> 협심증환자자료 = read.csv('협심증환자자료.csv')

> head(협심증환자자료)

  time censor freq

1  0.5      1  456

2  0.5      0    0

3  1.5      1  226

4  1.5      0   39

5  2.5      1  152

6  2.5      0   22

> attach(협심증환자자료)

2. lifetab 함수 연구

> lifetab

function (tis, ninit, nlost, nevent)

(이하 생략)

-> tis: 시간, 길이는 17, ninit = 전체 데이터, 2418명, nlost = 중도절단, nevent 사건발생 

3. nevent 자료 만들기

> which(censor==1) #사망 데이터

 [1]  1  3  5  7  9 11 13 15 17 19 21 23 25 27 29 31

> 집단.사망=협심증환자자료 [which(censor==1),]

> head(집단.사망)

   time censor freq

1   0.5      1  456

3   1.5      1  226

5   2.5      1  152

7   3.5      1  171

9   4.5      1  135

11  5.5      1  125

 4. nlost 자료 만들기

> which(censor==0) #절단 데이터

 [1]  2  4  6  8 10 12 14 16 18 20 22 24 26 28 30 32

> 집단.절단=협심증환자자료[which(censor==0),]

> head(집단.절단)

   time censor freq

2   0.5      0    0

4   1.5      0   39

6   2.5      0   22

8   3.5      0   23

10  4.5      0   24

12  5.5      0  107

5. ninit전체데이터 자료 만들기

> 사망자수=집단.사망[,3]

> 사망자수

 [1] 456 226 152 171 135 125  83  74  51  42  43  34  18   9   6   0

> 절단자수=집단.절단[,3]

> 절단자수

 [1]   0  39  22  23  24 107 133 102  68  64  45  53  33  27  23  30

> 합=사망자수+절단자수

> 합

 [1] 456 265 174 194 159 232 216 176 119 106  88  87  51  36  29  30

> sum(합) #2418

[1] 2418

6. tis 시간 변수 자료 만들기

> 년 = floor(집단.사망$time) #floor 내림값

> 년

 [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15

> lt=length(년)

> lt

[1] 16

> length(년)

[1] 17

> 년[lt+1] = NA

> 년[lt+1]

[1] NA

7. 생명표 만들기 

> 생명표=lifetab(년,sum(합),절단자수,사망자수)

> 생명표

      nsubs nlost  nrisk nevent      surv        pdf     hazard     se.surv

0-1    2418     0 2418.0    456 1.0000000 0.18858561 0.20821918 0.000000000

1-2    1962    39 1942.5    226 0.8114144 0.09440394 0.12353102 0.007955134

2-3    1697    22 1686.0    152 0.7170105 0.06464151 0.09440994 0.009179397

3-4    1523    23 1511.5    171 0.6523689 0.07380423 0.11991585 0.009734736

4-5    1329    24 1317.0    135 0.5785647 0.05930618 0.10804322 0.010138361

5-6    1170   107 1116.5    125 0.5192585 0.05813463 0.11859583 0.010304216

6-7     938   133  871.5     83 0.4611239 0.04391656 0.10000000 0.010379949

7-8     722   102  671.0     74 0.4172073 0.04601094 0.11671924 0.010450930

8-9     546    68  512.0     51 0.3711964 0.03697464 0.10483042 0.010578887

9-10    427    64  395.0     42 0.3342218 0.03553750 0.11229947 0.010717477

10-11   321    45  298.5     43 0.2986843 0.04302654 0.15523466 0.010890741

11-12   233    53  206.5     34 0.2556577 0.04209376 0.17941953 0.011124244

12-13   146    33  129.5     18 0.2135639 0.02968456 0.14937759 0.011396799

13-14    95    27   81.5      9 0.1838794 0.02030570 0.11688312 0.011765989

14-15    59    23   47.5      6 0.1635737 0.02066194 0.13483146 0.012259921

15-NA    30    30   15.0      0 0.1429117         NA         NA 0.013300258

           se.pdf   se.hazard

0-1   0.007955134 0.009697769

1-2   0.005975178 0.008201472

2-3   0.005069200 0.007649121

3-4   0.005428013 0.009153696

4-5   0.004945997 0.009285301

5-6   0.005033980 0.010588867

6-7   0.004690538 0.010962697

7-8   0.005175094 0.013545211

8-9   0.005024599 0.014659017

9-10  0.005307615 0.017300846

10-11 0.006269963 0.023601647

11-12 0.006847514 0.030646128

12-13 0.006682743 0.035110295

13-14 0.006514794 0.038894448

14-15 0.008035120 0.054919485

15-NA          NA          NA

-> 5번째 열 surv는 생존 함수, 7번째 열 hazard는 위험 함수. 협심증 환자의 약 19%가 1년 이내, 약 28%가 2년 이내에 사망하는 것으로 추정된다. 따라서 협심증 환자가 2년 이상 생존할 확률은 약 72%가 된다는 것을 알 수 있다. 또한 5년 이상 생존할 확률은 약 52%이다.

#전체 환자 중 생존확율이 70% 되는 지점

> names(생명표)

 [1] "nsubs"     "nlost"     "nrisk"     "nevent"    "surv"      "pdf"     

 [7] "hazard"    "se.surv"   "se.pdf"    "se.hazard"

> which.min(abs(생명표$surv-0.7))

[1] 3

> 생명표[3,]

    nsubs nlost nrisk nevent      surv        pdf     hazard     se.surv

2-3  1697    22  1686    152 0.7170105 0.06464151 0.09440994 0.009179397

       se.pdf   se.hazard

2-3 0.0050692 0.007649121

-> abs:절대값, which.min 최소값, 또는 생명표에서 3번째 행을 보면 survival이 70%에 가장 가까이 접근했음을 알 수 있다. 2년 이상 생존할 확률은 72%

#전체 환자 중 생존확율이 50% 되는 지점

> which.min(abs(생명표$surv-0.5))

[1] 6

> 생명표[6,]

    nsubs nlost  nrisk nevent      surv        pdf    hazard    se.surv

5-6  1170   107 1116.5    125 0.5192585 0.05813463 0.1185958 0.01030422

        se.pdf  se.hazard

5-6 0.00503398 0.01058887

-> 5년이 지나면 surv 즉 생존할 확률이 약 52%가 됨을 알 수 있다.

8. 데이터 시각화 – 생존함수 추정치 그래프

 > plot(년[1:lt], 생명표[,5],type='s',xlab='year',ylab='Survival function',ylim=c(0,1),lwd=2)

> abline(h=0.5)

> abline(v=5)

> plot(년[1:lt], 생명표[,5],type='o',xlab='year',ylab='Survival function',ylim=c(0,1),lwd=2)

> abline(h=0.5)

> abline(v=5)         

9. 데이터 시각화 – 위험함수 추정치 그래프

> names(생명표)

 [1] "nsubs"     "nlost"     "nrisk"     "nevent"    "surv"      "pdf"     

 [7] "hazard"    "se.surv"   "se.pdf"    "se.hazard"

> mean(생명표$hazard)

[1] NA

> 생명표$hazard

 [1] 0.20821918 0.12353102 0.09440994 0.11991585 0.10804322 0.11859583

 [7] 0.10000000 0.11671924 0.10483042 0.11229947 0.15523466 0.17941953

[13] 0.14937759 0.11688312 0.13483146         NA

> mean(생명표$hazard,na.rm=T) #NA 값 제외한 평균

[1] 0.1294874

> plot(년[1:lt], 생명표[,7],type='s',xlab='년',ylab='Hazard function',ylim=c(0,0.25),lwd=2)

> abline(h=0.1294874)

> plot(년[1:lt], 생명표[,7],type='o',xlab='년',ylab='Hazard function',ylim=c(0,0.25),lwd=2)>

> abline(h=0.1294874)

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