[R] Looking for Post-hoc tests (a la TukeyHSD) or interaction-level independent contrasts for survival analysis.
Thomas Oliver
toliver at stanford.edu
Tue Apr 29 23:12:33 CEST 2008
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
More information about the R-help
mailing list