coxph models with frailty

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'))

These frailty models require multiple observations to share the same 
frailty value.  The model with independent observations and frailties is 
barely identifiable, and we aren't using maximum likelihood here, so it 
isn't surprising that it doesn't work.

Changing the simulation to
  v=rep(rgamma(n/2,shape=1,scale=1),2)     #gamma frailty
  id=rep(seq(1:(n/2)),2)

I get

fit1

Call:
coxph(formula = Surv(t, c) ~ x + frailty(id, dist = "gamma",
     sparse = TRUE, method = "em"))

                           coef  se(coef) se2   Chisq DF p
x                         -1.21 0.206    0.184  34.8  1 3.7e-09
frailty(id, dist = "gamma                      265.4 62 0.0e+00

Iterations: 6 outer, 56 Newton-Raphson
      Variance of random effect= 0.837   I-likelihood = -838.7
Degrees of freedom for terms=  0.8 62.0
Likelihood ratio test=230  on 62.8 df, p=0  n= 200

so the estimation of the variance of the random effect is reasonable.

 	-thomas