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

[tlumley at u.washington.edu]

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