D.2 Manual Calculation

y1 = d1[d1$DV == 1, "TIME"] ; y1 # Death event times 

 [1]  81  11  43  16  42  46  24  43  46  49  96  12  75  89  80  38  87  83  67   6  76
[22]  45  96  79 100  38  58  61  75  90  63  78  16  30  59  30  70  63  99  35  66  44
[43]  27  25  24  32  12  30

y2 = d1[d1$DV == 0, "TIME"] ; y2 # Censored times

[1] 100 100

D.2.1 Survival function and hazard function

St = function(y, a, b) exp(b/a*(1 - exp(a*y))) # Gompertz survival function
ht = function(y, a, b) b*exp(a*y)              # Gompertz hazard function

위의 자세한 설명은 다른 자료를 참고한다.

아래는 목적함수(objective function)로 사용할 minus log likelihood function 이다.

mll = function(th) 
{
# th[1] = a = shape = theta/gamma
# th[2] = b = rate/scale = lambda
  r1 = sum(log(St(y1, th[1], th[2])) + log(ht(y1, th[1], th[2]))) # for event
  r2 = sum(log(St(y2, th[1], th[2]))) # for censored
  return(-(r1 + r2))
}

다음은 위의 목적함수(objecteive function)를 범용 최소화 함수(general purpose minimization function)인 nlminb를 사용하여 적합(fitting)한 것이다. 범용 최소화 함수로 optim도 있지만, 이 경우 nlminb가 정밀도가 더 좋아 nlminb함수를 사용하였다.

r3 = nlminb(c(0.01, 0.01), mll) ; r3

$par
[1] 0.02786 0.00516

$objective
[1] 230

$convergence
[1] 0

$iterations
[1] 17

$evaluations
function gradient 
      29       41 

$message
[1] "relative convergence (4)"

생존함수(Survival function)의 모수(parameter)인 a와 b의 점추정치(PE)는 약 0.02786과 0.00516로 추정되었다. 하지만, nlminb함수는 optim함수와 달리 hessian matrix를 output으로 주지 않는다. 그래서, 추정값의 표준오차(standard error)를 구하기 위해, 다음과 같이 numDeriv package의 hessian함수를 사용하여 hessian matrix를 구하고, 이를 이용하여 표준오차(standard error)를 구하였다.

require(numDeriv)
H0 = hessian(mll, r3$par)
SE0 = sqrt(diag(solve(H0))) ; SE0

[1] 0.00552 0.00169