3.3 시간 변수 (Time-to-event Variable)

R 기본 package인 survival package가 많이 사용되지만, 생명표(life table) 작성을 위해서는 KMsurv package의 lifetab이 좋다.

다음은 KMsurv의 lifetab 함수의 예제 자료이다. ?명령어로 자료와 설명을 볼 수 있다.

require(KMsurv)
?lifetab # copy and paste the example

다음은 위 함수 도움말에 나온 예제이다.

tis = c(0, 2, 3, 5, 7, 11, 17, 25, 37, 53, NA) # length 11
nsubs = c(927, 848, 774, 649, 565, 449, 296, 186, 112, 27) # length 10
nlost = c(2, 3, 6, 9, 7, 5, 3, rep(0, 3))
nevent = c(77, 71, 119, 75, 109, 148, 107, 74, 85, 27)
t1 = lifetab(tis, nsubs[1], nlost, nevent) ; t1
      nsubs nlost nrisk nevent    surv      pdf  hazard  se.surv    se.pdf se.hazard
0-2     927     2 926.0     77 1.00000 0.041577 0.04338 0.000000 0.0045368  0.004939
2-3     848     3 846.5     71 0.91685 0.076900 0.08755 0.009074 0.0087684  0.010380
3-5     774     6 771.0    119 0.83995 0.064821 0.08363 0.012058 0.0055430  0.007639
5-7     649     9 644.5     75 0.71030 0.041329 0.06178 0.014947 0.0045695  0.007120
7-11    565     7 561.5    109 0.62765 0.030460 0.05375 0.015967 0.0027313  0.005118
11-17   449     5 446.5    148 0.50581 0.027943 0.06622 0.016593 0.0020898  0.005335
17-25   296     3 294.5    107 0.33815 0.015357 0.05550 0.015812 0.0013853  0.005231
25-37   186     0 186.0     74 0.21529 0.007138 0.04139 0.013826 0.0007904  0.004660
37-53   112     0 112.0     85 0.12964 0.006149 0.07644 0.011358 0.0006305  0.006560
53-NA    27     0  27.0     27 0.03125       NA      NA 0.005912        NA        NA

 

tmed = (tis[1:10] + c(tis[2:10], 70))/2 # Last time is somewhat arbitrary one.

생존 확률 graph와 hazard rate graph는 다음과 같이 그릴 수 있다. 

oPar = par(mfrow=c(1, 2))
plot(tmed, t1$surv, type="o", xlab="Time (month)", ylab="Survival Probability",
     cex.axis=0.8, cex.lab=0.8)
plot(tmed, t1$hazard, type="o", xlab="Time (month)", ylab="Hazard Rate",
     cex.axis=0.8, cex.lab=0.8)
Survival and hazard plot

Figure 3.11: Survival and hazard plot

 

MASS package의 VA dataset으로 전체 자료의 survival은 다음과 같이 분석할 수 있다. 

require(MASS)     # for VA dataset
require(survival) # for survfit and Surv function
r3 = survfit(Surv(stime, status==1) ~ 1, VA) ; r3
Call: survfit(formula = Surv(stime, status == 1) ~ 1, data = VA)

       n events median 0.95LCL 0.95UCL
[1,] 137    128     80      52     105
# summary(r3) # Suppressed because the output is too long.

Graph는 다음과 같다. 중도 절단(censoring)부분을 +로 표시하기 위해서 mark.time=TRUE라는 option을 사용한다. 기본적으로 신뢰구간도 함께 그려준다.

plot(r3, mark.time=TRUE, xlab="Time (day)", ylab="Probability of Survival")
Survival curve of a group

Figure 3.12: Survival curve of a group