[R] coxph models with frailty
Roberto Perrelli
perrelli at uiuc.edu
Mon Nov 8 18:20:21 CET 2004
Dear R users:
I'm generating the following survival data:
set.seed(123)
n=200 #sample size
x=rbinom(n,size=1,prob=.5) #binomial treatment
v=rgamma(n,shape=1,scale=1) #gamma frailty
w=rweibull(n,shape=1,scale=1) #Weibull deviates
b=-log(2) #treatment's slope
t=exp( -x*b -log(v) + log(w) ) #failure times
c=rep(1,n) #uncensored indicator
id=seq(1:n) #individual frailty indicator
group=rep(1:(n/10), 10) #shared frailty indicator
Then I'm using the survival package in R 1.9.1 to estimate
the following Cox models with frailty:
fit1=coxph(Surv(t,c)~x+frailty
(id,dist='gamma',sparse=TRUE,method='em'))
fit2=coxph(Surv(t,c)~x+frailty
(group,dist='gamma',sparse=TRUE,method='em'))
fit3=coxph(Surv(t,c)~x+frailty
(id,dist='gamma',sparse=TRUE,method='aic',caic=TRUE))
fit4=coxph(Surv(t,c)~x+frailty
(group,dist='gamma',sparse=TRUE,method='aic',caic=TRUE))
In all cases, and after several replications, I am getting
estimates of the variance of the random effect that are
almost zero, whereas I thought that they should be around 1
(the variance of the gamma frailty in my data generating
process). Am I misunderstanding the procedures in some way,
or is this a known feature?
PS: Why the difference between the "penalty" (e.g. fit1
$history$frailty$theta) and the "variance of the random
effect" reported in the gamma frailty models above?
Cordially,
Roberto Perrelli
Department of Economics
University of Illinois
484 Wohlers Hall
1206 South Sixth Street
Champaign, Illinois 61820
USA
More information about the R-help
mailing list