I have the following coxph model:
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 predict() and 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 lme fit)?
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?
