I'm attempting use lsmeans
and its contrast for an F-test on an interaction. Basically, I'd like to replicate what Stata does with its contrast
command. I'd like to do this for two reasons:
within a regression model that has an interaction between factor variables;
within an ANOVA to help decompose a three-way interaction.
For this question, I'll ask about the three-way interaction.
library(haven)
threeway <- read_spss("https://stats.idre.ucla.edu/wp-content/uploads/2016/02/threeway.sav")
threeway$ID <- row.names(threeway)
library(afex)
three_fit <- aov_ez("ID", "y", data = threeway, between = c("a", "b", "c"))
three_fit
> three_fit
Anova Table (Type 3 tests)
Response: y
Effect df MSE F ges p.value
1 a 1, 12 1.33 112.50 *** .90 <.0001
2 b 1, 12 1.33 0.50 .04 .49
3 c 2, 12 1.33 47.84 *** .89 <.0001
4 a:b 1, 12 1.33 120.13 *** .91 <.0001
5 a:c 2, 12 1.33 6.84 * .53 .01
6 b:c 2, 12 1.33 8.47 ** .59 .005
7 a:b:c 2, 12 1.33 6.97 ** .54 .010
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1
The three-way interaction is significant. Now, from the UCLA page on Stata, Stata can use the code:
contrast b#c@a
This will produce an F-test for the b*c interaction at the levels of a.
I'm trying to do the same thing with lsmeans
in R. But, I just can't get it. Here's what I've tried:
library(lsmeans)
lsm <- lsmeans(three_fit, c("b", "c"), by="a")
test(contrast(lsm, "consec"), joint=TRUE)
That gets me an F-test, but it is not correct (or at least it is not the one I want). Any help in replicating the Stata results would be appreciated. I'd really like to stay within lsmeans
to do this, but if something else works, I'll take it.