[R] coxph models with frailty

[Roberto Perrelli]

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