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
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> fit1Call: 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