7

I have a lmer model with three-way interaction and I want to set up a specific contrast testing for the significance of two-way interaction on each level of the third variable. I can do it by hand with a simple model, but I was hoping that there might be a more efficient way of doing this.

Here's an example:

library(lsmeans)
library(lme4)

#setting up data.frame
mydata<-data.frame(expand.grid(subjid=1:10, A=c('0','1'), B=c('X','Y'), C=c('A','B','C','D'))) 
mydata$dv<-rnorm(nrow(mydata))
# and here is the model
fit<-lmer(dv~A*B*C+(1|subjid), mydata)

The interaction I want to test can be written symbolically as (X.0-Y.0)-(X.1-Y.1)|C, that is, a test for significance of the differences by B between levels of A at each level of C.

I can do it by including an interaction term instead of A and B in the model and setting up the contrasts by hand:

mydata$BA<-interaction(mydata$B,mydata$A)
fit<-lmer(dv~BA*C+(1|subjid), mydata)

lsmf<-lsmeans(fit, c('BA','C'))

c_list <- list(c1 = c(0.5, -0.5, -0.5, 0.5, rep(0,12)),
               c2 = c(rep(0,4),0.5, -0.5, -0.5, 0.5, rep(0,8) ),
               c3 = c(rep(0,8),0.5, -0.5, -0.5, 0.5, rep(0,4) ),
               c4 = c(rep(0,12),0.5, -0.5, -0.5, 0.5 ))

summary(contrast(lsmf, c_list), adjust = "holm")

But this is a very clumsy way, especially if there are more than four levels of C or there are other factors in the model. Moreover, with manual coding of BA the model becomes slightly different I think. So is there a better way for setting up such contrasts?

Andrey Chetverikov
  • 522
  • 3
  • 11
  • Using the FIRST version of your 'fit' model, do `lsm = lsmeans(fit, ~A*B|C)` and then `contrast(lsm, list(c = c(1,0,0,-1))` – Russ Lenth Nov 24 '16 at 21:33
  • Oops, mis-read the question. The contrasts you want are obtained via `contrast(lsm, list(con = c(-1,1,-1,1))`. As documented, `lsm` remembers the `by` spec from the construction (or you can specify it explicitly if you like), and you only need to specify the contrast coefficients within a level of the `by` variable(s). – Russ Lenth Nov 25 '16 at 00:23
  • Thanks, @rvl! I did not know that `by` is remembered. The correct contrasts then seem to be c(0.5,-0.5,-0.5,0.5) as I am comparing average values. At least that gives me the same results as with manual interaction coding. Could you post your response as an answer so that I can accept it? – Andrey Chetverikov Nov 30 '16 at 12:40
  • Oops, I messed-up again. The really easy way to get THAT contrast is to specify `contrast(lsm, interaction = TRUE, "pairwise")` -- albeit it won't be divided by 2. – Russ Lenth Dec 01 '16 at 14:59
  • Thanks again! I'd be glad to close this question if you'd post your comment as an answer. – Andrey Chetverikov Dec 07 '16 at 14:02

1 Answers1

2

This can be done in the lsmeans package pretty simply:

lsm = lsmeans(fit, ~A*B|C)
contrast(lsm, interaction = "pairwise")

This code generates and tests the contrast with coefficients $(1,-1,-1,1)$ at each level of factor $C$. This contrast is generated by taking the product of coefficients $(1,-1,1,-1)$ (for factor $A$) and $(1,1,-1,-1)$ (for factor $B$).

Russ Lenth
  • 15,161
  • 20
  • 53
  • 1
    By the way, this question led me to consider the idea that it'd be useful to be able to see the contrast coefficients that were generated. So the *next* update of **lsmeans** (available in maybe a month or so) will include a `coef` method for obtaining these. – Russ Lenth Dec 07 '16 at 23:12