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.
I think what you're looking for is
That is, you want to test interaction contrasts, not comparisons among means. You can see what it is you're testing (I recommend that!) by running that
contrast
call without wrapping it intest
. Try it with and without theinteraction
keyword.