On Tue, 6 May 2003, Jerome Asselin wrote:
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).).
Looking at a wider context in the code
pfun <- function(coef, theta, ndeath) {
if (theta == 0)
list(recenter = 0, penalty = 0, flag = TRUE)
else {
recenter <- log(mean(exp(coef)))
coef <- coef - recenter
nu <- 1/theta
list(recenter = recenter, first = (exp(coef) - 1) *
nu, second = exp(coef) * nu, penalty = -sum(coef) *
nu, flag = FALSE)
}
}
so the recentering means the penalty is actually your penalty (2) -- not
surprising, as the code was written by Terry Therneau.
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.
The log density penalty doesn't give maximum likelihood (which you would get by integrating out the frailties). It gives the joint likelihood of the data and the random effects. For the gamma model, as you note, they are equivalent. I believe that the current state of knowledge is that the log density penalty generally gives consistent estimates but is not equivalent to maximum likelihood. However, I have not actually seen the arguments for consistency and it isn't obvious to me how the argument would go. -thomas