[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