Solved – Reporting linear regression with Post-hoc comparisons

anovapost-hocregression

I have a quantitative variable (Dim: which is actually factor scores) and 2 categorical variables (Category + Region) in my dataset. I am interested in knowing if the variance explained by each categorical variable (and their interaction) is significant. Additionally I also want to know how much overall variance is explained by the model (i.e. the R squared value). After consulting materials about this topic I came to know that following 2 commands in R can fulfill my requirements.

my_lm <- lm(Dim4 ~ Category*Region, data=df)
summary(my_lm)
Anova(my_lm, type = "III")

I'm also interested in reporting whether group comparisons within each category are significantly different from each other. For this purpose, TukeyHSD has been mentioned in the reference material.

TukeyHSD(aov(my_lm), ordered=TRUE)

This command outputs the whole range of group comparisons for example:

$Category
                            diff          lwr      upr     p adj
Tweets-Conversations   0.6843803 -0.279052730 1.647813 0.2961601
Comments-Conversations 1.6552203  0.665253423 2.645187 0.0000540
$Region
          diff       lwr      upr p adj
US-PK 1.146865 0.7036846 1.590046 5e-07

$`Category:Region`
                                        diff         lwr      upr     p adj
Conversations:US-Tweets:PK        0.04138269 -1.46974284 1.552508 1.0000000
Conversations:PK-Tweets:PK        0.30818926 -1.31628692 1.932665 0.9998612
Comments:PK-Tweets:PK             1.20670619 -0.35032424 2.763737 0.2931398
Tweets:US-Tweets:PK               1.78097509  0.24027735 3.321673 0.0097448

While I am interested in all comparisons of Category and Region individually, in the interaction subgroups I am only interested in certain comparisons, e.g. Tweets:US-Tweets-PK, Conversations:US-Conversations-PK.
I have 2 questions regarding reporting the results of these tests. Firstly,How should I report the results of linear regression model (e.g. in case of Anova there would be F-statistic, p value and degrees of freedom like this One-way ANOVA: F(9, 3085) = 70.23, p < .001;). But I am looking at my data in terms of linear regression and want to report r squared + significance of each variable and their interaction as well. Does that involve reporting degrees of freedom and F statistic?
Secondly, as mentioned above I am only interested in certain group comparisons when it comes to interaction effects, can I just report the selected group comparisons? Another relevant question is if I'm not interested in all comparisons, can I limit the output of post hoc test to my desired group comparisons only?

Best Answer

You could take a look at the emmeans (estimated marginal means) package in R. It provides a lot of flexibility for these kinds of purposes. Since it appears you do have an interaction of consequence, I recommend against doing main-effects comparisons of either factor, averaged over the other. Instead, get a good display of how the interaction plays out, e.g. an interaction plot:

require("emmeans")
emmip(my_lm, Category ~ Region)

and display the EMMs and comparisons of one factor's levels separately by the other:

emm = emmeans(my_lm, ~ Category * Region)
emm
pairs(emm, by = "Category")
pairs(emm, by = "Region")

If you do not want all pairwise comparisons even with a by variable, there are several other built-in contrast types, or you can specify your own. For example, to compare consecutive levels of Category for each Region, do

contrast(emm, "consec", by = "Region")

There are a lot of examples and discussion in the several vignettes in the package. Start with vignette("basics")

Related Solutions

Solved – How to get only desirable comparisons from post-hoc

Suppose that we have one response variable and one explanatory variable (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)
TukeyHSD(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 matrices (overview):

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

First - decide which comparisons you are interested in.

Then create a contrast matrix:

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

look at this table for reference:

contrasts(con.data$category)
    [,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 in a factor, there can be only 4 comparisons, due to degrees of freedom.

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

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

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

In the 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)) 


  Coefficients:
                    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 the 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 to make a contrast matrix containing only comparisons you are interested in. With the example above it should not be difficult.

However !!!

I would not use post-hoc (or contrasts) on data immediately. It is like to teach a 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 to make full model:

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

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

Suppose the results of m1 model will look 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 interactions are non-significant - let's remove them (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 the 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 factor

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) has no sense, because the triple and double interaction are non-significant as well. And it would not be correct to select the desired rows from TukeyHSD (even if such post-hoc will show 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 factor is non-significant. Do not be sad - even non-significant result is a result.

Marginality rule is the practical application of Occam's razor (See also in Crawley, M.J. Statistics: An Introduction Using R. 2nd ed. Wiley, 2014, Ch. 10 Multiple Regression, p. 195). A good model is always the simplest one - and explains the largest portion of variability in data.

You can publish this result in an article as a "full model" (containing all the factors and their interaction(s)) or as a "minimal adequate model" (MAM, containing only significant effects). I would prefer to include both versions into a manuscript and let reviewers to decide which one to prefer.

The point is not to fishing for p-values in post-hoc tests when ANOVA results are non-significant.

Solved – Interpreting a significant 2-way interaction with post-hoc tests

Pooling across the remaining group was correct. You don't need anymore tests. Take choice 1. Your interaction means the effect of difficulty (hard - easy) in controls is different from the effect of difficulty in clinical. Since what it means is exactly what you want to know a post hoc is completely unnecessary.

You might want to see Gelman and Stern (2006) for a lesson in why the path you were thinking of taking with post hocs was faulty logic.

Gelman, A., & Stern, H. (2006). The difference between “significant” and “not significant” is not itself statistically significant. The American Statistician, 60(4), 328-331. doi: 10.1198/000313006X152649

Some general information might help others in the same situation. There seems to be a broadly accepted mistaken belief that there must be a significant simple effect when you have a significant interaction. Not only is this false but there is no guarantee of any significant simple effect in an ANOVA with more than 2 levels (in which case it's just an alternative to a t-test). Interactions are about differences among differences. None of the individual differences needs to be significant, or even close to significant, for a significant interaction. Consider a situation where one effect is +5 and another is -5. The difference between those effects is 10, which is much larger than either simple effect. Therefore, it would be very easy to not have a significant simple effect of 5 but a significant interaction of 10.

With something as simple as a 2x2 designs one should never perform a post hoc on the interaction because the interaction told you all of the necessary information. Even with slightly more complex designs post hoc tests are usually unnecessary, especially when one variable only has 2 levels. People generally do far too many post hoc tests following ANOVAs. In school they tell you to do an ANOVA to avoid multiple comparison problems. Of what benefit is that when, after you perform the ANOVA, you perform a bunch of comparisons? Interpret the ANOVA first.

Best Answer

Related Solutions

Solved – How to get only desirable comparisons from post-hoc

Solved – Interpreting a significant 2-way interaction with post-hoc tests

Related Question