I have the following
Surv(Y) ~ genotype + treatment + age + frailty(batch) + genotype:frailty(batch)
batch is a random effect (two replicates of the experiment were done with different subjects). Events can only happen once per subject.
The summary output looks like this:
coef se(coef) se2 Chisq DF p genotypeTG -0.2504 0.44224 0.43550 0.32 1.00 5.7e-01 treatmentX -0.6695 0.16379 0.16370 16.71 1.00 4.4e-05 age 0.0135 0.00839 0.00829 2.59 1.00 1.1e-01 frailty(batch) 1.39 0.94 2.2e-01 genotypeTG:frailty(batch) -1.3479 0.40020 0.39353 11.34 1.00 7.6e-04 genotypeTG:frailty(batch) 0.00000 0.00000 1.00
My question is: given that there is a significant interaction between genotype and batch, is the genotype coefficient a summary of the two batches, or is it just for the batch the model is using as the reference group?
If the former, how do I extract the batch-specific estimates? If the latter, how do I construct the overall coefficient and its SE? Thanks.
Update: by playing around with
coef() I found:
> coef(myfit) genotypeTG treatmentX age gamma:batch1 gamma:batch2 genotypeTG:frailty(batch)1 genotypeTG:frailty(batch)2 -0.25044857 -0.66945067 0.01350297 -0.17817716 0.15117985 -1.34789198 NA > cbind(mynewdata,predict(myfit,mynewdata)) genotype treatment age batch predicted 1 WT Z 0 batch1 0.64416111 2 TG Z 0 batch1 -0.95417945 3 WT X 0 batch1 -0.02528956 4 TG X 0 batch1 -1.62363012 5 WT Z 0 batch2 0.97351812 6 TG Z 0 batch2 0.72306955 7 WT X 0 batch2 0.30406745 8 TG X 0 batch2 0.05361887
...assuming no bugs in
predict.coxph.penal(), this means that the first three coefficients are ONLY for
batch1. For groups that have a level which interacts with the frailty term, I add the interaction coefficient for that batch in addition to the coefficient/s for those variables. To get
batch2 estimates, I do the above but also SUBTRACT the
batch1 frailty intercept, add the
batch2 frailty intercept, and then add other coefficients as appropriate. The interaction term for
batch2 is NA but apparently it doesn't get added.
My question now reduces to, how do I get the overall coefficients and SEs (i.e. the equivalent of what
summary.lme() would report if this were an
Also, should I be concerned that one of the frailty interaction coefficients is NA and that the summary omits some of the frailty coefficients?
Or, am I creating unnecessary problems for myself by allowing an interaction between a frailty term and a fixed effect? My reasoning was that bidirectional
stepAIC() identified this interaction as improving the fit of the model; this implies that the response of transgenic organisms may have been systematically different from that of normal ones between the two replicate experiments. Therefore, ignoring the interaction would amount to ignoring this source of error by treating the experiments as though only the baseline hazard differs between them. Is there a flaw in this reasoning somewhere?