R help - Dec 2007 - coxme frailty model standard errors?

Hello,
I am running R 2.6.1 on windows xp 
I am trying to fit a cox proportional hazard model with a shared
Gaussian frailty term using coxme
My model is specified as:

nofit1<-coxme(Surv(Age,cen1new)~ Sex+bo2+bo3,random=~1|isl,data=mydat)

With x1-x3 being dummy variables, and isl being the community level
variable with 4 levels.

Does anyone know if there is a way to get the standard error for the
random effect, like in nofit1$var?  I would like to know if my random
effect is worth writing home about.

Any help would be most appreciated
Corey Sparks

I can get the following output
nofit1<-coxme(Surv(Age,cen1new)~ Sex+bo2+bo3,random=~1|isl, data=no1901)
nofit1
Cox mixed-effects model fit by maximum likelihood
  Data: no1901 
  n=959 (2313 observations deleted due to missingness)
  Iterations= 3 69 
                    NULL Integrated Penalized
Log-likelihood -600.0795  -581.1718 -577.9682

  Penalized loglik: chisq= 44.22 on 5.61 degrees of freedom, p= 4.3e-08 
 Integrated loglik: chisq= 37.82 on 4 degrees of freedom, p= 1.2e-07 

Fixed effects: Surv(Age, cen1new) ~ Sex + bo2 + bo3 
         coef exp(coef)  se(coef)    z      p
Sex 0.2269214  1.254731 0.2151837 1.05 0.2900
bo2 0.5046991  1.656487 0.2510523 2.01 0.0440
bo3 1.0606144  2.888145 0.2726000 3.89 0.0001

Random effects: ~1 | isl 
                isl
Variance: 0.3876189



Corey Sparks
Assistant Professor
Department of Demography and Organization Studies
University of Texas-San Antonio
One UTSA Circle
San Antonio TX 78249
Phone: 210 458 6858
corey.sparks at utsa.edu

Hello,
I am running R 2.6.1 on windows xp
I am trying to fit a cox proportional hazard model with a shared
Gaussian frailty term using coxme My model is specified as:

nofit1<-coxme(Surv(Age,cen1new)~ Sex+bo2+bo3,random=~1|isl,data=mydat)

With x1-x3 being dummy variables, and isl being the community level
variable with 4 levels.

Does anyone know if there is a way to get the standard error for the
random effect, like in nofit1$var?  I would like to know if my random
effect is worth writing home about.

Any help would be most appreciated
Corey Sparks

I can get the following output
nofit1<-coxme(Surv(Age,cen1new)~ Sex+bo2+bo3,random=~1|isl, data=no1901)
nofit1 Cox mixed-effects model fit by maximum likelihood
  Data: no1901
  n=959 (2313 observations deleted due to missingness)
  Iterations= 3 69 
                    NULL Integrated Penalized Log-likelihood -600.0795
-581.1718 -577.9682

  Penalized loglik: chisq= 44.22 on 5.61 degrees of freedom, p= 4.3e-08
Integrated loglik: chisq= 37.82 on 4 degrees of freedom, p= 1.2e-07 

Fixed effects: Surv(Age, cen1new) ~ Sex + bo2 + bo3 
         coef exp(coef)  se(coef)    z      p
Sex 0.2269214  1.254731 0.2151837 1.05 0.2900
bo2 0.5046991  1.656487 0.2510523 2.01 0.0440
bo3 1.0606144  2.888145 0.2726000 3.89 0.0001

Random effects: ~1 | isl 
                isl
Variance: 0.3876189

Corey Sparks
Assistant Professor
Department of Demography and Organization Studies
University of Texas-San Antonio
One UTSA Circle
San Antonio TX 78249
Phone: 210 458 6858
corey.sparks at utsa.edu

--- begin included message
I am running R 2.6.1 on windows xp
I am trying to fit a cox proportional hazard model with a shared
Gaussian frailty term using coxme My model is specified as:

nofit1<-coxme(Surv(Age,cen1new)~ Sex+bo2+bo3,random=~1|isl,data=mydat)

With x1-x3 being dummy variables, and isl being the community level
variable with 4 levels.

Does anyone know if there is a way to get the standard error for the
random effect, like in nofit1$var?  I would like to know if my random
effect is worth writing home about.

-- end included message

  Computation of the se of the random effect turns out to be very hard, and 
worse it isn't worth much when you have done so.  For both these reasons
coxme
doesn't even try (and it's not on the list of things to add).
  
   First, you can do a likelihood ratio test by comparing the fit to a coxph 
model that does not have the random effect.   Make sure that the null 
loglikelihood for the two fits are the same though!  (For instance, one obs was 
missing the "isl" variable above, so the two fits have differnt n).  

     fit1<-  coxme(Surv(Age,cen1new)~ Sex+bo2+bo3,random=~1|isl, data=mydat)
     fit2<- coxph(Surv(Age,cen1new)~ Sex+bo2+bo3, data=mydat)
     fit1$loglik[1:2] - fit2$loglik	
     
The first number printed should be 0, twice the second is distributed as a chisq
on "number of random effects" degrees of freedom.  (Not quite.  The
chisq test
is actually conservative since the random effect is constrained to be >=0. 
But
it is close, and I'll let someone else work out the details)

   Second, you can print a profile likelihood confidence interval.  You had a 
random effects variance of about .4, so make a guess that the confidence 
interval is somewhere between 0 and 2.
     vtemp <- seq(0.0001, 2, length=20) 
     ltemp <- 0*vtemp
     for (i in 1:length(vtemp)) {
     	tfit <- coxme(Surv(Age,cen1new)~ Sex+bo2+bo3,random=~1|isl, data=mydat,
     		    variance= vtemp[i])
        ltemp <- 2 * diff(tfit$loglik[1:2])
        }
     plot(vtemp, ltemp, xlab='Variance', ylab="LR test")
     abline(h= 2*diff(fit1$loglik[1:2]) - 3.84)
     
   The sequence of fits each have a fixed value for the variance of the random 
effect; the plot shows the profile likelihood.  The profile is often very 
asymmetric (which is why the se of the random effect isn't worth much).  The
intersection of the profile with a line 3.84 units down (chisq on 1df) is the 
profile based 95% confidence interval. 
   
   	Terry Therneau

While it may be true that for coxme models the "standard errors"
are not very good approximations, it is always useful to have them
to compare with other diagnostics such as likelihood ratios and
profile likelihoods.  It is interesting to hear that with the currently
used methodology

   "Computation of the se of the random effect turns out to be very
hard"

because if you simply use AD Model Builders Random Effects module to
formulate the model you will get the standard errors calculated for free
without any more effort.

    Cheers,

      Dave

-- 
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-rsch.com

Hi, but why we do the difference : ltemp <- 2 * diff(tfit$loglik[1:2]) ??

Where I can find information about Integrate Likelihooh and null like
lihood??


Thank you very much,
Roby

--
View this message in context:
http://r.789695.n4.nabble.com/coxme-frailty-model-standard-errors-tp842203p3695827.html
Sent from the R help mailing list archive at Nabble.com.