[Rd] frailty models in survreg() -- survival package (PR#2933)
jerome at hivnet.ubc.ca
jerome at hivnet.ubc.ca
Wed May 7 01:37:02 MEST 2003
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
More information about the R-devel
mailing list