데이터마이너를 꿈꾸며

비모수적 방법 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)

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