# [R] coxph models with frailty

Thomas Lumley tlumley at u.washington.edu
Mon Nov 8 22:30:58 CET 2004

```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

```