[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
More information about the R-help
mailing list