[R] coxme with frailty--variance of random effect?

Thomas Chadefaux chadefaux at gmail.com
Fri Feb 3 12:25:33 CET 2012


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



More information about the R-help mailing list