D.3 Using flexsurv Package

R에 기본적으로 내장되어 있는 survival package에는 gompertz distribution을 지원하지 않는다. 따라서, flexsurv package의 flexsurvreg함수를 이용하였다.

require(flexsurv)

Loading required package: flexsurv

r1 = flexsurvreg(Surv(TIME, DV) ~ 1, data=d1, dist="gompertz") ; r1

Call:
flexsurvreg(formula = Surv(TIME, DV) ~ 1, data = d1, dist = "gompertz")

Estimates: 
       est      L95%     U95%     se     
shape  0.02786  0.01705  0.03866  0.00551
rate   0.00516  0.00272  0.00979  0.00169

N = 50,  Events: 48,  Censored: 2
Total time at risk: 2755
Log-likelihood = -230, df = 2
AIC = 463

위의 결과를 보면 우리가 구성한 목적함수(objective function, 여기서는 minus log likelihood)의 결과가 nlminb 함수를 이용하여 구한 것과 동일하다.

각 관찰시점(observation time point)에서의 생존확률(survival probability)은 다음과 같이 볼 수 있다.

summary(r1) 

 
   time    est    lcl   ucl
1     6 0.9669 0.9403 0.982
2    11 0.9358 0.8887 0.964
3    12 0.9291 0.8782 0.960
4    16 0.9012 0.8352 0.943
5    24 0.8385 0.7461 0.902
6    25 0.8300 0.7346 0.896
7    27 0.8125 0.7114 0.883
8    30 0.7852 0.6764 0.863
9    32 0.7662 0.6524 0.849
10   35 0.7366 0.6164 0.826
11   38 0.7058 0.5820 0.801
12   42 0.6627 0.5347 0.767
13   43 0.6517 0.5236 0.757
14   44 0.6405 0.5122 0.748
15   45 0.6292 0.5004 0.738
16   46 0.6177 0.4889 0.728
17   49 0.5829 0.4535 0.696
18   58 0.4741 0.3530 0.584
19   59 0.4618 0.3423 0.570
20   61 0.4371 0.3203 0.545
21   63 0.4125 0.2981 0.518
22   66 0.3757 0.2642 0.480
23   67 0.3636 0.2517 0.469
24   70 0.3275 0.2232 0.433
25   75 0.2697 0.1746 0.374
26   76 0.2585 0.1659 0.364
27   78 0.2367 0.1470 0.339
28   79 0.2261 0.1384 0.327
29   80 0.2156 0.1292 0.317
30   81 0.2054 0.1211 0.307
31   83 0.1856 0.1061 0.283
32   87 0.1489 0.0778 0.243
33   89 0.1321 0.0637 0.219
34   90 0.1241 0.0577 0.211
35   96 0.0821 0.0296 0.167
36   99 0.0650 0.0202 0.142
37  100 0.0598 0.0174 0.136