6
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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)) $\endgroup$ – rvl Nov 24 '16 at 21:33
  • $\begingroup$ 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). $\endgroup$ – rvl Nov 25 '16 at 0:23
  • $\begingroup$ 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? $\endgroup$ – Andrey Chetverikov Nov 30 '16 at 12:40
  • $\begingroup$ 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. $\endgroup$ – rvl Dec 1 '16 at 14:59
  • $\begingroup$ Thanks again! I'd be glad to close this question if you'd post your comment as an answer. $\endgroup$ – Andrey Chetverikov Dec 7 '16 at 14:02
2
$\begingroup$

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$).

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – rvl Dec 7 '16 at 23:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.