Hi
I would really appreciate some help with my code for coxme...
My data set
I'm interested in survival of animals after an experiment with 4
treatments, which was performed on males and females. I also have two
random factors:
Response variable: survival (death)
Factor 1: treatment (4 levels)
Factor 2: sex (male / female)
Random effects 1: person nested within day (2 people did the
experiment over 2 days)
Random effects 2: box nested within treatment (animals were kept in
boxes in groups of 6, and there were multiple boxes per treatment)
I've read the introductions to coxme by Terry Therneau, and something
like the following is what I think I should use:
model1<-coxme(Surv(death,censor)~treatment*sex+(1|day/person)+(1|
treatment/box))
Which gives me the following:
Cox mixed-effects model fit by maximum likelihood
events, n = 154, 291
Iterations= 57 305
NULL Integrated Penalized
Log-likelihood -823.276 -795.2354 -784.4807
Chisq df p
AIC BIC
Integrated loglik 56.08 11.00 4.9096e-08 34.08 0.67
Penalized loglik 77.59 17.91 2.0958e-09 41.78 -12.60
Model: Surv(death, censor) ~ treatment * sex + (1 | day/person) + (1
| treatment/box)
Fixed coefficients
coef
exp(coef) se(coef)
z p
teratmentb -0.0838877 0.9195345 0.3744511
-0.22 0.8200
treatmentb2 -0.4731922 0.6230103 0.3136199 -1.51
0.1300
treatmentn -1.0154149 0.3622521 0.4097754
-2.48 0.0130
sexmale -0.1838885 0.8320286 0.2602169
-0.71 0.4800
treatmentb:sexmale -0.3905856 0.6766605 0.2132936 -1.83
0.0670
treatmentb2:sexmale 0.6742202 1.9625020 0.3836907 1.76 0.0790
treatmentn:sexmale 1.2628977 3.5356520 0.4603589 2.74
0.0061
Random effects
Group Variable Std Dev Variance
day/person (Intercept) 0.32690104 0.10686429
day (Intercept) 0.49516113 0.24518455
treatment/box (Intercept) 0.26837158 0.07202330
treatment (Intercept) 0.29263637 0.08563604
My questions
(1) Does anyone know how I can test the significance of each of the
random effects in turn? For example, to find the significance of (1|
treatment/box) can I compare the results of the above model to one
without this term? (i.e. by subtracting the integrated loglikelihood
values of the model without (1|treatment/box) from the model with
that term).
(2) Can I include 'treatment' as a factor as well as including it as
part of a nested term? (incidentally I did wonder if I should include
it as (treatment|box), but an error message comes back that factors
cannot be used as a covariate within a random effect)
(3) Is it possible to test the proportionality assumption within
coxme. Previously I used >cox.zph(model1) with coxph, but that does
not work with coxme.
Very many thanks to anyone who can offer me some advice!
Sophie
[[alternative HTML version deleted]]
______________________________________________
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.
