Solved – Random slopes for interactions without random slopes for main effects

glmmlme4-nlme

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.

Thanks!

Best Answer

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

set.seed(101)
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,
                 weights=rep(1,nrow(dd)),
                 newparams=list(beta=rep(1,4),
                                theta=rep(1,10)))[[1]]
library(lme4)

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.

Related Solutions

Solved – Testing significance of random effects in a linear mixed-effects model with interactions

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.

Solved – Maximal model for linear mixed-effects model for repeated mesaures design

The maximal structure would need to include also a random effect for the interaction between color and shape, that is:

Y ~ color * shape + (color + shape + color:shape | subject)

This will result in all your predictors (color, shape and their interaction) having a fixed effect (constant for all subjects), and a random effect (individual fluctuations around the estimated fixed effect). In this sense the model is the maximal one. Note that it might not be fully equivalent to a repeated-measures ANOVA as it doesn't make equally strict assumptions on the correlational structure (see Tom's answer).

If you don't include the interaction in the random effect part of the formula, individual variation in the interaction effect will not be considered as "random", and the model will not be equivalent to a repeated-measures ANOVA. Of course, the variance of the random deviates for the interaction (or any other random effect) might be so small that including it in the model do not improve much the fit. You can check this not only with the AIC, but with a likelihood ratio test, as model with vs without one random effect are nested one another. In principle if the likelihood ratio test is not significant, it means that you can safely remove that random effect. Simplifying the random effect structures by removing negligible components would be an example of what in the article you linked is called data-driven approach.

You can simplify the model in this way, and it would still be equivalent to a repeated-measures ANOVA:

Y ~ color*shape + (1|subject) + (0+color|subject) + (0+shape|subject) + (0+color:shape|subject)

This syntax tells lmer to not estimate the correlations of random deviates across subjects. The drawback here is that, for example, you won't be able to tell whether subjects that have a large effect of color tend to have also a larger effect of shape (or smaller effect, in case of negative correlation).

You can easily include a between-subjects predictor, the only difference is that you can't add a random effect for it. "gender" for example cannot have a random effect grouped according to subject, but it can interact with the other fixed effects, e.g.:

Y ~ color * shape * gender + (color + shape + color:shape | subject)

Best Answer

Related Solutions

Solved – Testing significance of random effects in a linear mixed-effects model with interactions

Solved – Maximal model for linear mixed-effects model for repeated mesaures design

Related Question