데이터마이너를 꿈꾸며

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일임을 알 수 있다.

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

p.46 제2장 보건 정보 데이터의 기초 분석

2.3 자료의 기술 및 요약

> setwd('c:/Rwork/')

> 담즙과포화비율자료 = read.table('담즙과포화비율.txt',header=T)

> attach(담즙과포화비율자료)

> head(담즙과포화비율자료,2)

성별 담즙의과포화비율

1 남자               40

2 남자               88

> plot(담즙의과포화비율,type='p',xlab='자료',ylab='담즙과포화비율',main='담즙과포화비율')

> par(new=T)

> plot(담즙의과포화비율,type='h',xlab='자료',ylab='담즙과포화비율',main='담즙과포화비율')

> length(담즙의과포화비율)

[1] 60

> length(담즙의과포화비율) #길이

[1] 60

> sort(담즙의과포화비율,decreasing=T) #내림차순

[1] 146 142 137 128 127 123 123 120 118 116 112 111 110 110 107 106 106  98  91  90  89  88  88

[24]  88  87  87  86  86  86  84  84  82  80  80  80  79  78  77  76  76  75  74  73  73  69  67

[47]  66  66  65  65  58  58  57  56  55  52  52  47  40  35

> sum(담즙의과포화비율) #합

[1] 5185

> cumsum(담즙의과포화비율) #누적합

[1]   40  128  238  295  381  518  596  707  795  875  961 1041 1088 1194 1259 1333 1399 1478

[19] 1536 1659 1746 1834 1924 1980 2053 2165 2275 2393 2445 2551 2618 2683 2735 2819 2905 2940

[37] 3056 3132 3187 3260 3349 3476 3563 3705 3782 3858 3916 4007 4114 4212 4340 4424 4570 4645

[55] 4765 4845 4927 5050 5116 5185

> mean(담즙의과포화비율);median(담즙의과포화비율) #평균과 중앙값

[1] 86.41667

[1] 84

> mean(담즙의과포화비율,trim=1/10) #10% 절삭

[1] 85.4375

> var(담즙의과포화비율);sd(담즙의과포화비율) #분산과 표준편차

[1] 657.1624

[1] 25.63518

> fivenum(담즙의과포화비율) #다섯숫자요약

[1]  35.0  68.0  84.0 106.5 146.0

> quantile(담즙의과포화비율)

0%    25%    50%    75%   100%

35.00  68.50  84.00 106.25 146.00

> IQR(담즙의과포화비율) #3사분위수-1사분위수

[1] 37.75

> quantile(담즙의과포화비율)[4]-quantile(담즙의과포화비율)[2]

75%

37.75

> mad(담즙의과포화비율) #Median Absolut Deviation 각 데이터에서 중앙값을 뺀 후 절대값을 취한 값들의 중앙값

[1] 26.6868

> max(담즙의과포화비율);min(담즙의과포화비율)

[1] 146

[1] 35

> range(담즙의과포화비율)

[1]  35 146

> R=max(담즙의과포화비율)-min(담즙의과포화비율)

> R                          

[1] 111

2. 4 표와 그래프를 이용한 자료의 요약 

> head(담즙과포화비율자료,2)

성별 담즙의과포화비율

1 남자               40

2 남자               88

> 계급 = cut(담즙의과포화비율, breaks=c(20,40,60,80,100,120,140,160))

> head(계급)

[1] (20,40]   (80,100]  (100,120] (40,60]   (80,100]  (120,140]

Levels: (20,40] (40,60] (60,80] (80,100] (100,120] (120,140] (140,160]

> table(계급)

계급

(20,40]   (40,60]   (60,80]  (80,100] (100,120] (120,140] (140,160]

2         8        18        15        10         5         2

> par(mfrow=c(1,1))

> hist(담즙의과포화비율,breaks=c(20,40,60,80,100,120,140,160),main='히스토그램')

> rug(담즙의과포화비율)

> stem(담즙의과포화비율) #나무 줄기 그림

The decimal point is 1 digit(s) to the right of the | 

2 | 5

4 | 072256788

6 | 556679334566789

8 | 000244666778889018

10 | 667001268

12 | 033787

14 | 26

> stem(담즙의과포화비율,scale=2) #줄기의 마디 2배로 늘리기

The decimal point is 1 digit(s) to the right of the | 

3 | 5

4 | 07

5 | 2256788

6 | 556679

7 | 334566789

8 | 000244666778889

9 | 018

10 | 667

11 | 001268

12 | 03378

13 | 7

14 | 26

> library(vioplot)

> vioplot(담즙의과포화비율,names='담즙의과포화비율',col='yellow')

> n=length(담즙의과포화비율)

> plot(sort(담즙의과포화비율),(1:n)/n,type='s',ylim=c(0,1),main='Ogive of bile supersaturation',ylab='ECDF',xlab='담즙과포화비율') #ogive 오자이브 그래프

> rug(담즙의과포화비율)

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