Adding true random effects to survreg is certainly on my list of useful
additions, but one with no start date in sight. That said, one can get an
alternate solution with
survreg(Surv(time, status) ~ x1 + x2 + frailty.gaussian(id, method='aic'))
Justification: one can view a random effects model as a penalized model, that
is, as
a. the addition of "factor(id)" - a coefficient b_i for each subject
b. a shrinkage penalty, -.5*k* sum(b_i^2), is added to the log-lik, and
we minimize the sum
c. the value of k is chosen to maximize an integrated likelihood, one
with the b's integrated out. 1/k is the variance of the random effect
The above code uses the AIC to choose k. You could also use a user-specified
degrees of freedom.
WARNING: Using "method='reml'" in the above won't work correctly. The fact
that no warning is given in this case is serious flaw in the survival package.
(In my defense, the 'penalty function' code for coxph and survreg was designed
to allow general user-written penalties; a side effect is that the penalty
functions can't easily tell which routine is calling them. Most, e.g.,
pspline() and ridge(), work for both coxph and survreg. But frailty with
either
an 'ml' or 'reml' argument computes the appropriate Cox model integral, not the
survreg one, and so gives nonsense answers when used with survreg.)
Terry Therneau
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
