[R] coxph models with frailty

[Thomas Lumley]

On Mon, 8 Nov 2004, Roberto Perrelli wrote:

> 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