[Rd] frailty models in survreg() -- survival package (PR#2933)

[jerome at hivnet.ubc.ca]

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

-- 

Jerome Asselin (Jérôme), Statistical Analyst
British Columbia Centre for Excellence in HIV/AIDS
St. Paul's Hospital, 608 - 1081 Burrard Street
Vancouver, British Columbia, CANADA V6Z 1Y6
Email: jerome at hivnet.ubc.ca
Phone: 604 806-9112   Fax: 604 806-9044