coxme with frailty--variance of random effect?

Thank you very much for taking the time to answer. This pointed me in
the right direction.
For those interested, I found a useful explanation  of the derivation
of the estimated variance of the random effect in Hosmer & Lemeshow's
Applied Survival Analysis (1999), pp.321-26 (I love your book, Dr.
Therneay, but I needed something easier...). They proceed in 4 steps:

1. Obtain the cumulative hazard function for each subject.
2. Choose an arbitrary value for the variance parameter of the
frailties (call it theta).
3. Compute for each subject an estimate of the value of their
frailties, USING this variance parameter theta:
frailty_i= \frac(1+\theta \times c_i}{1+\theta \times H_i} (formula on
p. 321), where H is the cumulative hazard for the subject. So if theta
is 0 (no variance), then frailty=1 (i.e., no excess frailty). As theta
goes to infinity, the estimated frailty is simply the ratio
1/(cumulative risk so far) or 1/(cumulative risk so far), depending on
whether the subject is still alive or not.
4. fit the proportional hazard model with the same covariates as
before, but this time including the frailties in the hazard function.
5. calculate a log-likelihood for the model.

Repeat this for all possible values of the frailties (in practice,
sweep through them according to some algorithm), and pick the one with
the highest log-likelihood.

SO IF I understand correctly, the frailties are calculated GIVEN a
variance parameter of the frailties, and not the reverse. I.e., theta
is chosen first, and the random effects are estimated accordingly.
Which is why the variance of the estimated random effects is NOT the
same as theta. Did I get this right?

Thanks again,
Thomas


Best wishes,
Thomas Chadefaux

In answer to the several questions:

1. Variance of the random effect: ?Your example is not reproducable
since you didn't give the random number seed. ?Instead I'll use one of
the data sets that comes with the survival package.

library(coxme)
fit <- coxme(Surv(tstart, tstop, status) ~ treat + (1|center), cgd)
VarCorr(fit)$center

Intercept
0.1643779

var(ranef(fit)$center)

Intercept
[1] 0.06261608

?Yes, it is true that for a random effects model, the estimated variance
of the random effect is not equal to var(estimated per center effects).
Exactly the same is true of a linear mixed effects model: try the same
with lmer
? ? ? ?> library(lme4)
? ? ? ?> lfit(status ~ treat + (1|center), cgd)
? ? ? ?> VarCorr(lfit)$center
? ? ? ?> var(ranef(lfit)$center)

Why? It's a statistical insight that took me a while, so I don't think I
can explain it over email. ?Find someone familiar with mixed efffects
models and have a chat.

2. coxph(.... +frailty(center)) and coxme give different results.
?Read the documentation. ?One of them is fitting a gamma distribution
for the random effect, the other a Gaussian.

Terry Therneau