[R] coxme frailty model standard errors?

[Terry Therneau]

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