I am confused on how the log-likelihood is calculated in a parametric
survival problem with frailty. I see a contradiction in the frailty() help
file vs. the source code of frailty.gamma(), frailty.gaussian() and
frailty.t().
The function frailty.gaussian() appears to calculate the penalty as the
negative log-density of independent Gaussian variables, as one would
expect:
frailty.gaussian
[...]
list([...], penalty = 0.5 * sum(coef^2/theta +
log(2 * pi * theta)), flag = FALSE)
[...]
Similarly, the frailty.t() appears to use the joint negative log-density
of Student t random variables. HOWEVER, frailty.gamma() uses:
frailty.gamma
[...]
list([...], penalty = -sum(coef) *
nu, flag = FALSE)
[...]
I would rather expect to see something like:
(1) penalty=sum(coef*nu-(nu-1)*log(coef)+lgamma(nu)-nu*log(nu))
which is the joint negative log-density of gamma variables. Alternately, I
could also expect to see something like this:
(2) penalty=sum(coef-exp(coef))*nu
which was shown to lead to the same EM solution as penalty (1) -- at least
in the case of a Cox proportional hazard model (Therneau and Grambsch,
2000. Modeling Survival Data, Extending the Cox Model. Springer, New York.
Page 254, Eq. (9.8).). Bare we me, I don't know whether this holds in the
case of a parametric model. I also have concerns about the validity of the
likelihood ratio tests obtained with the latter penalty function (2),
because this penalty is NOT equal to the negative log-likelihood (1).
Finally, it's not clear to me whether we gain significant convergence
speed and accuracy by using the penalty (2) as opposed to (1).
Furthermore, the help file for frailty() says,
"The penalised likelihood method is equivalent to maximum (partial)
likelihood for the gamma frailty but not for the others."
In the current state of the package, I would think that this should be the
other way around. That is,
"The penalised likelihood method is equivalent to maximum (partial)
likelihood for the gaussian and t frailty but not for the gamma."
However, my current comprehension of the problem leads me to recommend to
use the negative log-likelihood of gamma variables (2). Hence, both gamma,
Gaussian and t frailty would be equivalent to maximum (partial) likelihood.
Any comment on this issue would be much appreciated.
Sincerely,
Jerome Asselin
R 1.6.2 on Red Hat Linux 7.2
Package: survival
Version: 2.9-7