I have some experimental data that I'm trying to analyze. I have 1 response variable and 3 explanatory variables (these are factor variables). The explanatory variables are the presence of a disease(positive and negative), a genetic profile (X and Y), and whether or not an MRI contrast agent was given (YES and NO).

Structure of data looks like this:

     measurement   profile  disease contrast
1    -1.76269      X        NEG       YES
2    -0.34492      X        NEG       NO
3     0.57455      X        POS       YES
4     2.16539      X        POS       NO
            .      .          .        .
            .      .          .        .
            .      .          .        .
77   -1.76269      Y        NEG       YES
78   -0.34492      Y        NEG       NO
79    0.57455      Y        POS       YES
80    2.16539      Y        POS       NO

I looked into using ANOVA for this analysis but the post hoc Tukey HSD looks at all possible combinations of the explanatory variables so it makes far more comparisons than I actually care about.

We have some specific hypotheses that, e.g.: X.NEG.NO will differ from Y.NEG.NO, X.NEG.NO will differ form X.NEG.YES, X.NEG.NO will differ from X.POS.NO, etc.

(notice that each compared group "consist" from interaction of all three variables)

How to get only some specific comparisons from TukeyHSD? Is appropriate to make Is there a better approach?

Reproducible example:

my.data <- data.frame(measurement = rnorm(80),
                      my.profile = rep(c("X","Y"), each = 40),
                      my.disease = rep(c("NEG","NEG","POS","POS"), times=20),
                      my.contrast = rep(c("NO","YES"), times = 40))
  • $\begingroup$ I don't have much background in DOE, but I'll try to answer the best I can. I don't think this could be considered an orthogonal design. The full set of comparisons I need to make is: X.NEG.NO vs. Y.NEG.NO, X.POS.NO vs Y.POS.NO; X.NEG.NO vs X.NEG.YES Y.NEG.NO vs Y.NEG.YES; X.POS.NO vs X.POS.YES, Y.POS.NO vs Y.POS.YES; X.NEG.NO vs X.POS.NO, Y.NEG.NO vs Y.POS.NO; X.NEG.YES vs X.POS.YES, Y.NEG.YES vs Y.POS.YES $\endgroup$ – NoMoreData Mar 6 '14 at 14:20
  • $\begingroup$ All other comparisons are meaningless.I also have 3 different regions that I want to consider all these comparisons in, but they are separate so it doesn't make sense to compare across them. $\endgroup$ – NoMoreData Mar 6 '14 at 14:28
  • $\begingroup$ For some comparisons they do but for others they don't (e.g. n=11 vs n=9). $\endgroup$ – NoMoreData Mar 6 '14 at 14:34
  • $\begingroup$ Thanks! Do I perform diagnostics with the summary command? Also, is there a name for this process so that I can cite it when I write up my results? And when you say I can only compare based on profile and contrast, that means I should perform a bivariate test ignoring the disease factor? (e.g. t.test(data$measurement[data$profile=="X"], data$measurement[data$profile == "Y"] $\endgroup$ – NoMoreData Mar 7 '14 at 13:35
  • $\begingroup$ No, please edit away if you think it makes it clearer. $\endgroup$ – NoMoreData Mar 7 '14 at 15:16

Suppose that we have one response variable and one explanatory variables (5 levels).

con.data <- data.frame(x = c(rnorm(75), c(rnorm(75)+5)),
                      category = rep(c("A","B1","B2","C1","C2"), each=30))

If we do ANOVA followed by classical TukeyHSD post-hoc...

m1 <- aov(con.data$x ~ con.data$category); summary(m1)

           diff        lwr      upr     p adj
B1-A  0.1877877 -0.9109837 1.286559 0.9897357
B2-A  2.2050543  1.1062829 3.303826 0.0000013
C1-A  4.5951737  3.4964022 5.693945 0.0000000
C2-A  4.7790488  3.6802773 5.877820 0.0000000
B2-B1 2.0172666  0.9184952 3.116038 0.0000117
C1-B1 4.4073860  3.3086145 5.506157 0.0000000
C2-B1 4.5912611  3.4924896 5.690033 0.0000000
C1-B2 2.3901193  1.2913479 3.488891 0.0000001
C2-B2 2.5739944  1.4752230 3.672766 0.0000000
C2-C1 0.1838751 -0.9148963 1.282647 0.9905238

...we obtain all possible combinations.

If you want to make only comparison that are interesting to you, use CONTRASTS.

In R there are several default contrast matrixes:

treatment, helmert, sum , polynomial...but you can create your own.

First - decide in which comparisons you are interested.

Than create an contrast matrix:

contrasts(con.data$category) <- cbind(c(1,-1/4,-1/4,-1/4,-1/4),

look on this table...

    [,1] [,2] [,3] [,4]
A   1.00  0.0  0.0  0.0
B1 -0.25 -0.5  0.0 -0.5
B2 -0.25 -0.5  0.0  0.5
C1 -0.25  0.5  0.5  0.0
C2 -0.25  0.5 -0.5  0.0

Take notice that sum of each column is equal to 0.

If you have 5 levels of a factor, there can be only 4 comparisons, due to degrees of freedom.

In first column you compare mean of "A" with mean of all others groups.

In second column you compare only "B" group with "C" group. (B1+B2 vs. C1+C2)

In third column you compare only within "C" category (C1 vs. C2)

In fourth column you compare only within "B" category (B1 vs. B2)

To see the results re-make the ANOVA with created contrast matrix.

summary(lm(con.data$x, con.data$category)) 

                    Estimate Std. Error t value Pr(>|t|)    
  (Intercept)         2.5910     0.1258  20.599  < 2e-16 ***
  con.data$category1  -2.3534     0.2516  -9.355  < 2e-16 ***
  con.data$category2   3.4907     0.2813  12.411  < 2e-16 ***
  con.data$category3  -0.1839     0.3978  -0.462    0.645    
  con.data$category4   2.0173     0.3978   5.072 1.19e-06 ***

...and each row in obtained table corresponds to each comparison (column in contrast matrix) made. (i.e. con.data$category1 is significant so there is significant difference between mean of "A" vs. mean of all others groups...etc.)

In short:

Try do make contrast matrix containing only comparisons which are interesting to you. With example above it should not be difficult.

However !!!

I would not use post-hoc (or contrasts) on data immediately. It is like to teach model to "run" before it can "walk". So as a first thing I would create a model containing all variables. Subsequently, remove all non-significant variables (or their interaction) according to marginality rule. This procedure (reduction) will determine if all your desired comparisons are necessary.

so try do make full model:

m1 <- glm(measurement ~ my.profile*my.contrast*my.disease, data = my.data)
anova(m1, test = "F")

For example if variable "disease" will not be significant (alone or in interaction - it should not by included in post-hoc.

Suppose the results of m1 model will looks like this:

my.profile                        0.00012 **
my.contrast                       0.00231 *
my.disease                        0.07690 .
my.profile:my.contrast            0.26159  
my.profile:my.disease             0.07709 .
my.contrast:my.disease            0.21256  
my.profile:my.contrast:my.disease 0.23319  

Use the rule of marginality:

update no. 1: you can see that triple interaction is non-significant - let's remove it

m2 <- update(m1, ~.-my.profile:my.contrast:my.disease)

Now, show the anova of updated model:

anova(m2, test = "F")

my.profile                        0.00012 **
my.contrast                       0.00231 *
my.disease                        0.07690 .
my.profile:my.contrast            0.26159  
my.profile:my.disease             0.07709 .
my.contrast:my.disease            0.21256  

update no. 2: you can see that double interaction are non-significant - let's remove it (start with double interaction with highest p- value)

m3 <- update(m2, ~.-my.profile:my.contrast)

anova(m3, test = "F")

my.profile                        0.00012 **
my.contrast                       0.00231 *
my.disease                        0.07690 .
my.profile:my.disease             0.07709 .
my.contrast:my.disease            0.21256  

update no. 3: remove another double interaction

m4 <- update(m3, ~.-my.contrast:my.disease)

anova(m4, test = "F")

my.profile                        0.00012 **
my.contrast                       0.00231 *
my.disease                        0.07690 .
my.profile:my.disease             0.07709 .

update no. 4: remove last double interaction

m5 <- update(m4, ~.-my.profile:my.disease)

anova(m5, test = "F")

my.profile                        0.00012 **
my.contrast                       0.00231 *
my.disease                        0.07690 .

update no. 5: remove non-significant variable

m6 <- update(m5, ~.-my.disease)

anova(m6, test = "F")

my.profile                        0.00012 **
my.contrast                       0.00231 *

Model m6 is our final model.

Unfortunately it is obvious that making comparisons (Y.NEG.NO, X.NEG.NO and others) make no sense, because the triple interaction and double also are non-significant. And it would be non appropriate to make them by selection of desired rows from TukeyHSD (even if such post-hoc will show some significant difference !!!). Believe me, such approach will be very hard to defend in peer-review process. So you can make only comparison in profile (X vs. Y) and contrast (NO vs. YES). Disease variable is non-significant. Do not be sad - even non-significant result is a result.

Marginality rule is the practical application of Occam's razor. The good model is always the simplest one - and explain the large amount of variability present in data.

You can publish this result in article as a "full model" (containing all the variables and all their interaction) or as a "minimal adequate model" (MAM - only significant variables. I would prefer to send both version into manuscript and let reviewer to decide which one he like more.

Point is to not fishing for p-values in post-hoc, when ANOVA results are non-significant.

  • $\begingroup$ Great, this is has really helped set me on the right course. But for multiple factors, do I need multiple matrices, a higher dimensional array, or can it be done with one matrix the way you describe? $\endgroup$ – NoMoreData Mar 6 '14 at 16:38
  • $\begingroup$ smallData <- data.frame(measurement = c(rnorm(5),rnorm(7),rnorm(23-(5+7))), profile = c(rep("X",7),rep("Y",23-7)), contrast = c(rep("YES",10),rep("NO",23-10)),disease = c(rep("POS",14),rep("NEG",23-14))) $\endgroup$ – NoMoreData Mar 6 '14 at 16:59
  • $\begingroup$ I'm looking to model measurement ~ profile + contrast + disease but I'm only interested in the comparisons listed above. $\endgroup$ – NoMoreData Mar 6 '14 at 17:02
  • $\begingroup$ There are 82 rows overall. $\endgroup$ – NoMoreData Mar 6 '14 at 18:35
  • $\begingroup$ I've added some example code to reproduce the format of my data. Thanks! $\endgroup$ – NoMoreData Mar 6 '14 at 19:18

Running step on your model will also work.

# Full model
m1.0 = glm (logICPMS ~ SITE*ISLAND*ORG*ORGL*CAGE*Metal,
            data = M.m, na.action = na.exclude)

# this removes each interacting factor to give you the MAM model

# MAM model
m1.1 = glm(formula = logICPMS ~ ISLAND + ORG + CAGE + Metal + ISLAND:ORG + 
                                ISLAND:CAGE + ISLAND:Metal + ORG:Metal, 
           data = M.m, na.action = na.exclude)
  • $\begingroup$ Stepwise selection methods are generally considered a bad idea (see here). $\endgroup$ – gung - Reinstate Monica May 20 '15 at 17:53
  • $\begingroup$ Since making post hoc comparisons is not the usual purpose of stepwise regression, could you please explain how and why this recommendation answers the question? $\endgroup$ – whuber May 20 '15 at 19:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.