R help - Feb 2012 - coxme with frailty--variance of random effect?

Dear all,

This probably stems from my lack of understanding of the model, but I
do not understand the variance of the random effect reported in coxme.
Consider the following toy example:

#------------------------------- BEGINNING OF CODE
------------------------------------------------
library(survival)
library(coxme)

#--- Generate toy data:
d <- data.frame(id = c(1:100),                   # create 100 individuals
                group=c(rep(1:10, each=10)),   #assign each to one of 10 groups
                surv=rnorm(100),                      #give them a survival time
                x=rnorm(100),                           # a predictor
                event=rep(1,100))                     #everyone dies

#--- run the cox model w/ random effects (frailty) and return the results
m1.coxme <- coxme(Surv(d$surv, d$event) ~ d$x + (1|d$group)); m1.coxme

#--- the estimated frailties for each of the 10 groups are reported under
m1.coxme$frail

#--- put these estimated frailty back in the data.frame
# (or should I just do var(m1.coxme$frail$d.group) ? Either way, the
results are different from the ones reported in coxme)
for(i in attributes(m1.coxme$frail$d.group)$names){
    d$frail[d$group==i] <- m1.coxme$frail$d.group[i]
    print(i)
}

#--- calculate variance of frailty
var(d$frail)   # result is 0.000001163377

#--- This doesn't match either coxph or coxme calculations:
m1.coxme  #Here I get var of the d.group random effect = 0.0003988903
m1.coxph <- coxph(Surv(d$surv, d$event) ~ d$x + frailty(d$group));
summary(m1.coxph)  #here it gives me 0.0000128

#--- note also that these two (coxme and coxph don't yield the same
estimate of the variance of the frailty (nor of the fixed effects
coefficients, for that matter).
#--------------------------------- END OF CODE
--------------------------------------------

I am quite sure I am missing something simple, and I apologize if this
is the case. I just cannot see what it is...

Thank you and best,
Thomas

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