testing contrast in two-way ANOVA using multcomp

Question

I am trying to compare two levels of a factor aggregated over the levels of another one using the multcomp package.

I use a "between-subjects-version" of the dataset obk.long contained in the afex package which I computed like this:

library(dplyr)
library(afex)
data(obk.long)

obk <- obk.long%>%
group_by(id)%>%
mutate(treatment,value = mean(value))%>%
distinct(id)

obk <- data.frame(obk)

It is now possible to test several contrasts by doing the following.

fit <- lm(value~treatment*gender,data=obk)
summary(fit)

# test for difference in A and B when gender = F
K <- matrix(c(0, 1,-1,0,0,0),1)
t <- glht(fit, linfct = K)
summary(t)

So far so good, but what if I want to compare A versus B aggregated over the levels of gender? I was thinking about this approach:

K <- matrix(c(0, 0.5,-0.5,0,0.5,-0.5),1)
t <- glht(fit, linfct = K)
summary(t)

Alternatively, I could fit another model without gender

fit2 <- lm(value~treatment,data=obk)
summary(fit2)

K <- matrix(c(0, 1,-1),1)
t <- glht(fit2, linfct = K)
summary(t)

but the results are not the same. This is probably due to the interaction?

What is the right (best) way to test this contrast using multcomp?

Answer 1

This can be solved by using the ingenious combination of afex with lsmeans (and also multcomp if one desires so, but this is usually not necessary). Furthermore, thanks to afex functionality to aggregate automatically, dplyr is not needed.

library(afex)
require(lsmeans)
require(multcomp)
data(obk.long)

# Step 1: set up the model using afex
# but use return = "aov" to obtain an object lsmeans can handle.
fit <- aov.car(value~treatment*gender + Error(id),data=obk.long, return = "aov")


# Step 2: set up reference grid on which we 
# we can test any type of tests using lsmeans functionality:
(ref1 <- lsmeans(fit, c("treatment", "gender")))
##  treatment gender   lsmean        SE df lower.CL upper.CL
##  control   F      4.333333 0.8718860 10 2.390650 6.276016
##  A         F      4.500000 0.8718860 10 2.557317 6.442683
##  B         F      5.833333 0.6165165 10 4.459649 7.207018
##  control   M      4.111111 0.7118919 10 2.524917 5.697305
##  A         M      8.000000 0.8718860 10 6.057317 9.942683
##  B         M      6.222222 0.7118919 10 4.636028 7.808416
## 
## Confidence level used: 0.95 

# we simply define the contrasts as a list on the reference grid:
c_list <- list(c1 = c(0, -1, 1, 0, 0, 0),
               c2 = c(0, -0.5, 0.5, 0, -0.5, 0.5))

# because we want to control for Type I errors we test this using
# the Bonferroni-Holm correction
summary(contrast(ref1, c_list), adjust = "holm")
##  contrast   estimate        SE df t.ratio p.value
##  c1        1.3333333 1.0678379 10   1.249  0.4805
##  c2       -0.2222222 0.7757662 10  -0.286  0.7804
## 
## P value adjustment: holm method for 2 tests 

# alternatively, we can pass it to multcomp for even cooler corrections:
summary(as.glht(contrast(ref1, c_list)), test = adjusted("free"))
## Note: df set to 10
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Linear Hypotheses:
##         Estimate Std. Error t value Pr(>|t|)
## c1 == 0   1.3333     1.0678   1.249    0.359
## c2 == 0  -0.2222     0.7758  -0.286    0.780
## (Adjusted p values reported -- free method)

Note that the first contrast (c1: A and B when gender = F) correspond to your first result while the latter one does not. You must have made an error.

The second contrast you wish can also be obtained easier (simply ignore the other two rows):

pairs(lsmeans(fit, "treatment"), adjust = "none")
## NOTE: Results may be misleading due to involvement in interactions
##  contrast      estimate        SE df t.ratio p.value
##  control - A -2.0277778 0.8347673 10  -2.429  0.0355
##  control - B -1.8055556 0.7338014 10  -2.461  0.0336
##  A - B        0.2222222 0.7757662 10   0.286  0.7804
## 
## Results are averaged over the levels of: gender