[Rd] Re: frailty models in survreg() -- survival package (PR#2934)

[jerome at hivnet.ubc.ca]

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