I have a GLMM with the following form:

m1 <- glmer(Y ~ A*B + (A*B|Subject), data, family=binomial(logit))

where A and B are categorical factors with 2 levels, and Subject is my grouping factor, with a relatively small number of levels (e.g., 10). A likelihood ratio test indicates that the random slope for the interaction A:B makes a significant contribution to the fit of the model, so I keep it.

Let's say I have no prior hypothesis about the structure of the random effects, and I want to find the model that provides the best and most parsimonious description of the data. Does it makes sense to proceed testing random slopes for the main effects?

More specifically, is it appropriate to compare (with a likelihood ratio test) the model above (m1) with this reduced model where I removed the random component for one of the main effects?

m2 <- glmer(Y ~ A*B + (A + A:B|Subject), data, family=binomial(logit))

I would say that intuitively it makes sense, as m2 would be a model where the effect of B is pretty much constant for all subjects, while the magnitude of the interaction A:B varies from subject to subject. However I am wondering whether in practice one can distinguish between variability in the effect of B vs. variability in the interaction A:B. I know that, for fixed-effects predictors, it is generally wrong to fit a model that has an interaction term without including the main effects marginal to that interaction, so I guess that I am asking whether I should respect the marginality principle also for testing random slopes.



1 Answer 1


I'm not sure this makes sense.

First of all, it won't work the way you expect: A:B will expand to a model of the same dimensionality as A*B (and lme4 isn't quite as smart as it might be about handling the redundancies

dd <- expand.grid(A=factor(1:2),B=factor(1:2),Subject=1:40,rep=1:20)
dd$Y <- simulate(~A*B+(A*B|Subject),
                 newdata=dd, family=binomial,

Your original model:

m1 <- glmer(Y ~ A*B + (A*B|Subject), dd, family=binomial)

This gives convergence warnings, and ends up with a log-likelihood 0.003 less than the original (basically, numeric fuzz)

m2 <- update(m1, . ~ A*B+(A+A:B|Subject))

This is also equivalent:

m3 <- update(m1, . ~ A*B+(A:B-1|Subject))
all.equal(logLik(m1),logLik(m3),tolerance=1e-4) ## TRUE

If you really wanted to leave out the main effect of B you'd have to use dummy variables.

For the meaning of this: can you explain what the model would mean? In particular, the model that you intend (not what lme4 actually does) when you use (A+A:B|Subject) would have among-individual variation in the effects of A, no effects of B at the baseline value of A (assuming you are using the standard treatment contrasts), but variation in the effect of adding B to A. Another analogy would be to a model with fixed intercepts but variable slopes.

  • $\begingroup$ thanks! for the meaning, A indicates the experimental condition; B is a factor that has weak effects (with very small individual variation) in the condition corresponding to the baseline value of A, but larger effects (with larger individual variation) in the other condition. Would it make sense in this case to simplify the model by removing the random slope of B? I noticed that if I set manually the dummy variables (in your example I'd recode both factors in this way dd$A<-ifelse(dd$A=="2",1,0)) then lme4 seems to do what I want. Is that a correct approach? $\endgroup$
    – matteo
    Jun 15, 2016 at 7:08
  • $\begingroup$ You could, although I don't think you will gain very much by removing the random slope. Yes, coding the dummy variables in the way you suggest will work (I think); there's also a dummy() function in lme4 that may be helpful. $\endgroup$
    – Ben Bolker
    Jun 15, 2016 at 15:19

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.