[R] coxme with frailty

[Terry Therneau]

 A couple of clarifications for you.

1. I write mixed effects Cox models as exp(X beta + Z b), beta = fixed
effects coefficients and b = random effects coefficients.  I'm using
notation that is common in linear mixed effects models (on purpose).
  About 2/3 of the papers use exp(X beta)* c, i.e., pull the random
effects out of the exponent.  Does it make a difference?  Not much: b
will be Gaussian and c will be log-normal. However, it's much easier to
write down complicated (crossed, nested, ect) random effects in the
notation I chose.  The description you give from H&L is using the "c"
style.

 2. If var(b) is known, solution of beta and b is straightforward.  It's
a small modification of the ususal coxph internals for a fixed effects
model.  That is, a bit of c code that has been optimized and exercised
over a decade, it makes sense to re-use it.  So all the routines I know
of solve mixed effects in a nested fashion: an outer maximizer looks for
the parameters of var(b), and for each trial value calls the known
routine (which itself iterates) for the corresponding (beta, b) pair.
The H&L formula you quote only holds in a specific case.

Terry T

-------- begin included message ----------------------------------
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?