2
$\begingroup$

I would like to test the significance of all interactions in a 3 factor linear mixed-effects model. Factors A and B are fixed, and factor C is random. Using lmer, the full model is:

lmer(response ~ A*B + (1 + A*B|C), data.frame)

My question concerns getting p-values for the 2-way interactions and the main effect of C. I initially planned to use likelihood ratio tests to compare this full model to reduced models by removing one interaction at a time. But of course specifying the 3-way interaction automatically includes the 2-way interactions and main effect of C.

So I would like to know if there is a way to construct a reduced model that includes the 3-way interaction and omits only a single 2-way interaction (or the main effect of C).

But of course the LR-test has its limitations. I would appreciate any advice on other ways to obtain p-values for random effects in this model. I am aware of RLRsim, MCMC, bootMer. Are any of these preferable/more feasible for my situation?

$\endgroup$
2
  • $\begingroup$ stackoverflow.com/questions/13502815/… does this help $\endgroup$
    – jona
    Sep 2, 2015 at 12:14
  • $\begingroup$ Thanks, that's good to know. It doesn't seem to be working in my case, will have to investigate further. $\endgroup$
    – jjh
    Sep 3, 2015 at 12:30

1 Answer 1

3
$\begingroup$

tl;dr This is possible, but very tedious (especially for factors with more than 2 levels), and you might not want to do it after all. If you can decide exactly what you mean to test by testing the lower-level interactions, you can probably do it.

This is a little bit tricky partly for technical R-specific reasons, and partly for statistical/inferential reasons.

The statistically tricky part is that testing lower-level interactions in a model that also contains higher-level interactions is (depending on who you talk to) either (i) hard to do correctly or (ii) just plain silly (for the latter position, see part 5 of Bill Venables's "exegeses on linear models". The rubric for this is the principle of marginality. At the very least, the meaning of the lower-order terms depends sensitively on how contrasts in the model are coded (e.g. treatment vs. midpoint/sum-to-zero). My default rule is that if you're not sure you understand exactly why this might be a problem, you shouldn't violate the principle of marginality.

The technically tricky part is that R's formula language doesn't have a simple way to drop lower-order terms in the presence of higher-order interactions (in part because this stuff was all designed originally by researchers from the Nelder/Venables camp that doesn't think this is a sensible thing to do). I have done this in the past, in cases where I knew it was sensible, by constructing the model matrix, dropping the columns I didn't want, and using the remaining terms as explicit variables in the model. In particular, I had terms period (before vs. after) and ttt (treatment: control vs. removal) and wanted to measure the period:ttt interaction but setting the before-period difference between treatments to zero: for the default ("treatment") contrasts in R, this corresponds to setting $\beta_{\textrm{ttt}}$ to zero, or removing that column from the model matrix.

dd <- expand.grid(period=c("before","after"),
                  ttt=c("control","removal"))
## some trickery to get factor levels in the sensible order
dd[] <- lapply(dd,function(x) factor(x,levels=unique(x)))

You might think that ~period*ttt-period would be the way to specify the desired result, but it doesn't actually work:

colnames(model.matrix(~period*ttt-ttt,dd))
## [1] "(Intercept)"             "periodafter"            
## [3] "periodbefore:tttremoval" "periodafter:tttremoval" 

(neither does ~period+period:ttt, or any other sensible combination I could think of). Instead:

X <- model.matrix(~period*ttt,dd)
## also remove the intercept because this will get re-added
X <- X[,!colnames(X) %in% c("(Intercept)","tttremoval")]
colnames(X) <- c("period","period_ttt")  ## friendlier names

Now we can fit a model with y~period+period_ttt,data=X and get what we want.

How does this translate to your case?

the full model is:

lmer(response ~ A*B + (1 + A*B|C), data.frame)

My question concerns getting p-values for the 2-way interactions and the main effect of C.

If you want a p-value for the interaction of A:B you should probably compare the model above to

update(model, . ~ A + B + (1 + A + B | C))

I'm not sure what the "main effect of C" is -- is that the variation of the intercept among levels of C ? If so, you probably have to use the trick above to exclude the intercept from A*B, e.g.

update(model, . ~ A*B + (A2 + B2 + A2_B2 | C))

where A2, B2, and A2_B2 are dummy variables as above.

But it's up to you to make sure this makes sense.

$\endgroup$
1
  • $\begingroup$ Your answer is very helpful, thanks. Precisely the explanation I needed (and as far as I can tell is not available anywhere else). $\endgroup$
    – jjh
    Sep 3, 2015 at 12:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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