1

I have been trying to compare a set of interaction contrasts using emmeans() and contrast(). However, I am having trouble applying a custom contrast and then compare it between groups. Here I added an example dataset so you can replicate what I am working with.

In summary, this dataset contains 3 factors: group, location, scenario. Location and scenario are within subject variables and group is a between subject variable. I would like to know if there is a higher measurement in location A compared to B, C and D (I set a contrast in the form c(1, -1/2, -1/2)) and then compare this contrast between groups (a contrast of contrast). So far I managed to get the contrast for each group (working with a single scenario at the moment):

library(afex)
library(emmeans)

data <- data.frame(
  id = rep(c(100:189), each = 9),
  group = rep(c("W", "X", "Y"), each = 90 * 3),
  location = rep(c("A", "B", "C"), each = 3, times = 90),
  scenario = rep(c("alpha", "beta", "gamma"), times = 270),
  measurement = c(rnorm(length(id)))
)

test_aov <- aov_car(measurement ~ group * scenario * location + 
                      Error(id/scenario * location),
                    data = data)

location_a_vs_all <- c(1, -1/2, -1/2)

test_effects <- emmeans(test_aov,
                        specs =  ~ location | scenario | group, 
                        at    = list(scenario = c("alpha"))) 
    
test_contrast <- contrast(test_effects,
                          method = list("Location A vs All" = location_a_vs_all),
                          adjust = "none")

But if I try to set the argument interaction to True in the contrast() function I get an error:

 test_contrast <- contrast(test_effects,
                           method = list("Location A vs All" = location_a_vs_all),
                           interaction = T,
                           adjust = "none")

Error in contrast.emmGrid(object, interaction[[i]], by = vars[-pos], name = nm, : 'method' must be a list, function, or the basename of an '.emmc' function

The same happens if I put my custom contrast directly in the interaction argument.

test_contrast <- contrast(test_effects,
                          interaction = list("Location A vs All" = location_a_vs_all),
                          adjust = "none")

I looked at the documentation on Interaction analysis in CRAN but I am lost on how to correctly implement what I would like to do, so I would appreciate any pointers from you. Thanks in advance!

Ramiro Reyes
  • 535
  • 2
  • 7
  • When you pipe one call into another, we have little idea what went wrong or why. If you do the steps separately and show the results or error messages for each step, I may be willing to help. – Russ Lenth Sep 20 '22 at 20:46
  • I added the error messages for clarity and removed the pipes. On a separate note, I may be reading this the wrong way but that "I may be willing to help" sounds very condescending, let's try to keep things civil here. – Ramiro Reyes Sep 22 '22 at 01:31
  • 1
    I'm sorry that came across as condescending. I was just trying to say I wasn't sure I'd be able to help. – Russ Lenth Sep 22 '22 at 21:51
  • 1
    It looks like you have an answer that works. But interaction = TRUE requires named contrast(s), and is not of any use here anyway since you have only one primary factor, location. The other factors are "by" factors. You need at least two primary factors for the interaction spec to produce interaction contrasts. – Russ Lenth Sep 22 '22 at 21:58
  • My bad, I definitely misinterpreted that. And now that makes sense, I was mistakenly considering the group as a primary factor, but yes, it is a "by" factor in this analysis. Thanks for your time! (and for your other answers, they were very helpful to figure out what I was doing wrong here). – Ramiro Reyes Sep 24 '22 at 11:35

1 Answers1

0

The solution is to split the analysis in three steps:

library(afex)
library(emmeans)

data <- data.frame(
  id = rep(c(100:189), each = 9),
  group = rep(c("W", "X", "Y"), each = 90 * 3),
  location = rep(c("A", "B", "C"), each = 3, times = 90),
  scenario = rep(c("alpha", "beta", "gamma"), times = 270),
  measurement = c(rnorm(length(id)))
)

test_aov <- aov_car(measurement ~ group * scenario * location + 
                      Error(id/scenario * location),
                    data = data)

First I got the effects for the specific scenario:

test_effects <- emmeans(test_aov,
                        specs =  ~ location | group, 
                        at    = list(scenario = c("alpha"))) 

Then I make the contrast comparing one of the locations to the other 2 (which I referred as my custom contrast in the question setup), using the "treatment vs control" method:

location_con <- contrast(test_effects, 
                         method = "trt.vs.ctrl", 
                         name = "A vs Others", 
                         by = "group", 
                         ref = 2:3)

Finally I apply a contrast of contrast to this, to compare the previous result between group W and the others, defining W as my reference level:

group_con <- contrast(location_con, 
                      method = "trt.vs.ctrl", 
                      by = "A vs Others", 
                      ref = "W")
Ramiro Reyes
  • 535
  • 2
  • 7