--- 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
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
