[Rd] Re: frailty models in survreg() -- survival package (PR#2934)
jerome at hivnet.ubc.ca
jerome at hivnet.ubc.ca
Wed May 7 02:54:17 MEST 2003
SEE ALSO ORIGINAL POSTING IN PR#2933
On May 6, 2003 03:58 pm, Thomas Lumley wrote:
>
> 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.
Ok, I forgot to account for the recentering. However, correct me if I am
wrong, but if Therneau and Grambsch (2000, Eq. (9.8)) are right, I think
it should be:
recenter <- mean(exp(coef))
instead of
recenter <- log(mean(exp(coef)))
>
> 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.
I'm confused... Do you mean that the log-likelihood of all model
parameters (which is maximised) is NOT the sum of the conditional
log-likelihood of the data and the log-likelihood of 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.
>
Sorry, I don't understand why you say that they are equivalent for the
gamma model. Therneau and Grambsch (p.254) specifically note,
"... the penalized loglikelihood and the observed-data loglikelihood in
equation (9.4) have the same solution, although these two equations are
NOT equal to one another."
For that reason, I think that the help file should read,
> > "The penalised likelihood method is equivalent to maximum
> > (partial) likelihood for the gaussian and t frailty but not for the
> > gamma."
Also, is it clear whether the current gamma penalty function remains valid
even in the case of parametric survival analysis? Is it computationally
better than (1)? Finally, what is the meaning of the "loglik" value in
gamma frailty? Can it still be used to calculate likelihood ratio tests?
Cheers,
Jerome
More information about the R-devel
mailing list