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