데이터마이너를 꿈꾸며

비모수적 방법 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일임을 알 수 있다.

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