2

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:

  1. within a regression model that has an interaction between factor variables;

  2. 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.

gung - Reinstate Monica
  • 11,583
  • 7
  • 60
  • 79
Michael
  • 31
  • 3

1 Answers1

2

I think what you're looking for is

test(contrast(lsm, interaction = "consec"), joint=TRUE)

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 in test. Try it with and without the interaction keyword.

Russ Lenth
  • 5,922
  • 2
  • 13
  • 21
  • Thanks for the suggestion - but this only gets me one F-test results (oddly with 2, 12 df). I'm looking for one interaction F-test (each with 2,12 df) for each level of "a". – Michael Sep 06 '17 at 16:11
  • Did you try adding `by="a"` to the `test` call? – Russ Lenth Sep 06 '17 at 17:13