[R] Looking for Post-hoc tests (a la TukeyHSD) or interaction-level independent contrasts for survival analysis.

[Thomas Oliver]

Hello all R-helpers,

    I've performed an experiment to test for differential effects of 
elevated temperatures on three different groups of corals.  I'm 
currently performing a cox proportional hazards regression with 
censoring on the survivorship (days to mortality) of each individual 
in the experiment with two factors: Temperature Treatment (2 levels: 
ambient and elevated) and experimental group (3 levels: say 1,2,3).

In my experiment, all three groups survived equally well in the 
ambient control treatment, but  two of three of the groups succumbed 
to heat stress in the elevated temperature treatment.  I can see that 
the third group had a small degree of mortality, but it appears to be 
significantly less than the other two and may be not significantly 
different from the ambient controls.

I would like to ask three questions:  1) Is group 3 different from 
controls? 2) Is group 3 different from group 1 and/or group 2 in the 
elevated treatment? and 3) are groups 1 and 2 different from each 
other in the elevated treatment?

   Because I'm testing for differential effects among the elevated 
temperature treatment group, and "I've seen the data" by now,  the 
analysis would be easiest for me if I performed a responsible 
multiple comparisons test like TukeyHSD to see how each of the six 
Treatment:Group subgroups compared to each other.  TukeyHSD does not 
appear to be defined for outputs from the function coxph -- (see 
survival library).

cph <- coxph(Surv(DayToMort, Censor) ~ Treatment*Group, data=subb)

--> Does anyone know of an implementation of TukeyHSD that would 
work, or of another post-hoc multiple comparison test?

I believe that another responsible tack would be to clearly define 
the contrasts I'd like to make within the interaction term. However 
this has yet to work as fully as I'd like it.

  I've successfully set the contrasts matrix for the three-level 
factor "Group" following Crawley's The R Book.

cmat<-cbind(c(-1,1,0),c(0,-1,1))
contrasts(subb$Group)<-cmat
contrasts(subb$Group)

By setting these contrasts and then looking at the interaction terms 
in the coxph model, this allows me to compare groups _within_ each 
separate treatment, and confirms both that #2) that groups 1 and 2 
are not sig. different in the elevated treatment, and #3) the group3 
corals survived significantly better than the other groups in the 
elevated treatment. BUT it does not allow me to say if the group3 
survival is or is not different from its own control.(#1 above).

To make this comparison, I've tried setting the contrast matrix for 
the Treatment:Group interaction term, with no success.  Whenever I 
attempt to do so, I run the code below:

cmat<-cbind(c(-1,0,0,1,0,0),c(0,-1,0,0,1,0),c(0,0,-1,0,0,1),c(0,0,0,-1,1,0),c(0,0,0,0,-1,1)) 
#Build a matrix
rownames(cmat)<-rownames(contrasts(subb$Treatment:subb$Group))  #give 
cmat the correct rownames
colnames(cmat)<-colnames(contrasts(subb$Treatment:subb$Group)) 
#give cmat the correct colnames
contrasts(subb$Treatment:subb$Group)<-cmat           #try to assign cmat

and I get this error message:
Error in contrasts(subb$Treatment:subb$Group, ) <- cmat :
   could not find function ":<-"

Alternatively I could run:
options(contrasts=c("contr.sum","contr.poly"))
where:
  contrasts(subb$Treatment:subb$Group)

             [,1] [,2] [,3] [,4] [,5]
Con:1   1    0    0    0    0
Con:2    0    1    0    0    0
Con:3    0    0    1    0    0
Exp:1    0    0    0    1    0
Exp:2    0    0    0    0    1
Exp:3   -1   -1   -1   -1   -1

But even that doesn't appear to affect the output of :

cph <- coxph(Surv(DayToMort, Censor) ~ Treatment*Group, data=subb)


--> Is what I'm trying to do statistically invalid and R is trying to 
quietly save me from statistical destruction, or it is just being a 
pain?  Is there a way around it?

--> Any other suggestions?

Many Thanks in Advance,

Tom Oliver