# Tukey HSD Post Hoc test with the Anova() Command (Car-package)

I am trying to do a Tukey-Test with the Anova() command from the car-package. I am using Anova() because I have an unbalanced Dataset. A statistic professor from my Uni told me Anova() would solve this problem. However he didn't tell me, how one could do a Tukey Test with it afterwards.

Tukey.HSD from R; HST.Test from agricola and emmeans all don't seem to work. I already did some research and found the post below, which also doesn't provide an answer.

I am relatively new to R, so I might've missed something.

Below you see the two factorial code. I am trying to do a Tukey test for A and C. You can see how the significance changes depending on what factor comes first. The significance stays the same when using Anova from the car-package. The significance interaction A:C seems to stay the same with and with out an unbalanced Dataset. Only small difference is due to lm() oder aov(). So I just did the Tukey-test for it with an unbalanced dataset.

So is it possible to do a Tukey-test with Anova()? If so, how?

data_NIRS$$dv<-data_NIRS$$Sinapine
data_NIRS$$dv<-data_NIRS$$Sinapine
> m1 <-aov(data_NIRS$$dv~data_NIRS$$A*data_NIRS$$C) > m2 <-aov(data_NIRS$$dv~data_NIRS$$C*data_NIRS$$A)
> summary(m1)
Df Sum Sq Mean Sq F value   Pr(>F)
data_NIRS$$A 2 23.14 11.570 23.92 6.30e-10 *** data_NIRS$$C               2  11.04   5.522   11.42 2.16e-05 ***
data_NIRS$$A:data_NIRS$$C   4   3.87   0.967    2.00   0.0965 .
Residuals               178  86.10   0.484
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> summary(m2)
Df Sum Sq Mean Sq F value   Pr(>F)
data_NIRS$$C 2 17.69 8.846 18.29 5.98e-08 *** data_NIRS$$A               2  16.49   8.246   17.05 1.68e-07 ***
data_NIRS$$C:data_NIRS$$A   4   3.87   0.967    2.00   0.0965 .
Residuals               178  86.10   0.484
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> m3 <-lm(data_NIRS$$dv~data_NIRS$$A*data_NIRS$$C) > m4 <-lm(data_NIRS$$dv~data_NIRS$$C*data_NIRS$$A)
> anova(m3)
Analysis of Variance Table

Response: data_NIRS$$dv Df Sum Sq Mean Sq F value Pr(>F) data_NIRS$$A               2 23.141 11.5704 23.9196 6.299e-10 ***
data_NIRS$$C 2 11.045 5.5223 11.4163 2.163e-05 *** data_NIRS$$A:data_NIRS$C 4 3.870 0.9675 2.0001 0.09647 . Residuals 178 86.102 0.4837 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > anova(m4) Analysis of Variance Table Response: data_NIRS$$dv Df Sum Sq Mean Sq F value Pr(>F) data_NIRS$$C 2 17.693 8.8464 18.2883 5.978e-08 *** data_NIRS$$A 2 16.493 8.2463 17.0476 1.683e-07 *** data_NIRS$$C:data_NIRS$A   4  3.870  0.9675  2.0001   0.09647 .
Residuals               178 86.102  0.4837
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> Anova(m3)
Anova Table (Type II tests)

Response: data_NIRS$$dv Sum Sq Df F value Pr(>F) data_NIRS$$A             16.493   2 17.0476 1.683e-07 ***
data_NIRS$$C 11.045 2 11.4163 2.163e-05 *** data_NIRS$$A:data_NIRS$C 3.870 4 2.0001 0.09647 . Residuals 86.102 178 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > Anova(m4) Anova Table (Type II tests) Response: data_NIRS$$dv Sum Sq Df F value Pr(>F) data_NIRS$$C 11.045 2 11.4163 2.163e-05 *** data_NIRS$$A 16.493 2 17.0476 1.683e-07 *** data_NIRS$$C:data_NIRS$A  3.870   4  2.0001   0.09647 .
Residuals               86.102 178
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>


Relevant SO thread: How to do a Tukey HSD test with the Anova command (car package)

• In what way does emmeans not work ? Oct 18, 2021 at 19:58
• The statistic professor wrote this line of code befor we had the problem with the unbalanced dataset. lsm.A<-emmeans(m2, pairwise~A,adjust="Tukey", type="response",infer=TRUE);lsm.A est.A<-lsm.A$emmeans contr.A<-lsm.A$contrasts cld.A<-as.data.frame(cld(est.A, alpha=0.05, Letters=letters, adjust="tukey") ) Oct 19, 2021 at 9:06
• @Sal Mangiafico I get the following: Error in (function (object, at, cov.reduce = mean, cov.keep = get_emm_option("cov.keep"), : Can't handle an object of class “anova” Use help("models", package = "emmeans") for information on supported models. Oct 19, 2021 at 9:07
• Run emmeans on the model, not on the anova result. Oct 21, 2021 at 13:38

As to the use of emmeans here, this is really a programming question. The code just needs to be updated a little:

A = rep(c("I","II","III"), each=5)
Y = c(1,2,3,4,5,3,4,5,6,7,9,10,11,12,13)
m1 = lm(Y~A)

library(emmeans)

lsm.A<-emmeans(m1, ~A)
summary(lsm.A, type="response",infer=TRUE)
pairs(lsm.A)

library(multcomp)
cld(lsm.A,alpha=0.05,Letters=letters)


EDIT: model name changed from m2 to m1 to avoid confusion with the example in the comments.

• Hey, could you explain the function of your first 2 lines? Would I have to adjust them to my dataset?. Sorry... R is new to me. Oct 19, 2021 at 11:09
• Yes, the first two lines just make up some toy data so the example can be run as is. So, yes, instead you would use your actual data, loaded in whatever method you want (read.table, read.csv, etc.) Oct 19, 2021 at 11:16
• Awesome thanks! I will try it out. Oct 19, 2021 at 11:22
• I still get the same error. data_NIRS$dv<-data_NIRS$Sinapine m1 <-lm(data_NIRS$dv~data_NIRS$A*data_NIRS$C) m2<-Anova(m1) A=data_NIRS$A lsm.A<-emmeans(m2, ~A) summary(lsm.A, type="response",infer=TRUE) pairs(lsm.A) Oct 19, 2021 at 11:26
• I also dont know how to make the comments here more spaced out :D Oct 19, 2021 at 11:28