[R] Looking for Post-hoc tests (a la TukeyHSD) or interaction-level independent contrasts for survival analysis.
David Winsemius
dwinsemius at comcast.net
Wed Sep 16 22:10:30 CEST 2015
On Sep 15, 2015, at 12:09 PM, Huot, Matthieu wrote:
> Hi Tom
>
> I know the post is over 7-8 years old but I am having the same question. How to do a post-hoc test like TukeyHSD on coxph type output.
Create a new variable using the `interaction`-function, apply you contrasts to that object, and that should let you side-step all the errors of the person trying to follow Crawley's book.
--
David.
>
> Have you received any info in this matter?
> Thanks
> Matthieu
>
> Looking for Post-hoc tests (a la TukeyHSD) or interaction-level independent contrasts for survival analysis.
> Thomas Oliver toliver at stanford.edu <mailto:r-help%40r-project.org?Subject=Re%3A%20%5BR%5D%20Looking%20for%20Post-hoc%20tests%20%28a%20la%20TukeyHSD%29%20or%20interaction-level%0A%20independent%20contrasts%20for%20survival%20analysis.&In-Reply-To=%3C6.2.5.6.2.20080429131826.04deaeb0%40stanford.edu%3E>
> Tue Apr 29 23:12:33 CEST 2008
>
> * Previous message: [R] AIC extract and comparison <https://stat.ethz.ch/pipermail/r-help/2008-April/160916.html>
> * Next message: [R] for loop in nls function <https://stat.ethz.ch/pipermail/r-help/2008-April/160886.html>
> * Messages sorted by: [ date ]<https://stat.ethz.ch/pipermail/r-help/2008-April/date.html#160885> [ thread ]<https://stat.ethz.ch/pipermail/r-help/2008-April/thread.html#160885> [ subject ]<https://stat.ethz.ch/pipermail/r-help/2008-April/subject.html#160885> [ author ]<https://stat.ethz.ch/pipermail/r-help/2008-April/author.html#160885>
>
> ________________________________
>
> 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
>
>
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.
David Winsemius
Alameda, CA, USA
More information about the R-help
mailing list