I have a question regarding including my main predictors in the fixed as well as random effects term in a mixed model. I am using R and will provide a sample data set, but I am generally happy for theoretical input on this.

The data set is the following: It's a repeated measures dataset with two different interventions (food and training), each of them have two levels. All participants underwent each combination of conditions. Hormone samples were collected over the course of the day, but in this model I want to focus on the baseline differences. Age and BMI are included as covariates.

ID <- rep(letters[1:16],4)
training <- rep(rep(c("T1","T2"),each=16),2)
food <- c(rep(c("F1","F2"),each=16),rep(c("F2","F1"),each=16))
hormone <- rnorm(n = 64,mean=5,sd=2)
BMI <- as.numeric(rep(sample(18:35,size = 16,replace = T),4))
age <- as.numeric(rep(sample(40:60,size = 16,replace=T),4))
DF <- data.frame(ID,training,food,hormone,BMI,age)
indNA <- which(DF$hormone%in%sample(DF$hormone,2))
indNABMI <- which(DF$BMI%in%sample(DF$BMI,1))
DF$hormone[indNA] <- NA
DF$BMI <- scale(DF$BMI)
DF$age <- scale(DF$age)

In order to test for the interaction effect of both intervention types on baseline hormonal differences, I would have liked to setup the following model:

m1 <- lmer(hormone~training*food+age+BMI+(1+training*food|ID),DF)

and test it against this null-model (if you see any flaw in the model setup to begin with, I am of course happy to know):

m01 <- lmer(hormone~training+food+age+BMI+(1+training+food|ID),DF)

However, it is obvious I do not have enough observations to do so.

My questions now are:

  1. Would it be appropriate to reduce the random slopes term to this?

    m1red <- lmer(hormone~training*food+age+BMI+(1+training+food|ID),DF)
  2. As pointed out in this post here for example, Minimum number of levels for a random effects factor?, one should have at least 5-6 observations for a random effect to be included. Is this also true for the random slopes per subject (in my case, there are two observations per condition level per subject) and should my random slopes therefore be excluded entirely?

  • 1
    $\begingroup$ You are not obliged to have everything in the random part of the model which is also in the fixed. If I understand your last part you have as many random slopes as levels of ID. $\endgroup$
    – mdewey
    Jan 31 '18 at 13:47

(1) tl;dr yes, the reduced model (1+training+food|ID) is appropriate.

As you have realized, the (food $\times$ training $\times$ ID) interaction has one observation per combination: if you didn't know this already from thinking about the experimental design, you could cross-tabulate:

ctab <- with(DF,table(food,training,ID))
> all(ctab<=1)
[1] TRUE

This means that in a linear mixed model, which (in lme4 at least) always includes an observation-level residual term, the random effect for this interaction will be completely confounded with the residual variance term. (This is not true for GLMMs where the scale parameter is fixed to a theoretical value, such as Poisson or binomial models, where an observation-level random effect can be used to incorporate overdispersion in the model.) There's another example/description of this issue here.

(2) the issue described about minimal numbers of observations has to do with the number of groups, not with the numbers of observations per group, so you might be OK. In general one would want more groups (IDs) for a more complex random-effects model. That is, 5-6 groups might be OK (barely) for estimating a scalar random effect (1|ID); you would want more for a random effect with independent terms, i.e. (1+training+food||ID) (note double bar: 3 variance parameters estimated), and more still for a random effect correlated terms, (1+training+food|ID) (3*(3+1)/2) = 6 variance-covariance parameters estimated). I don't know of a well-supported/referenced rule of thumb. The most obvious consequences of overfitting will be singular fits, which are most easily detectable when you see variances of 0 or correlations of +/-1, see the GLMM FAQ for more details. (In this case your example gives a singular fit, but that's because it's constructed with all-null results; your real data might do better.)

  • $\begingroup$ +1 but (1+training+food||ID) when training and food are categorical (in this case binary) factors will not only estimate 3 variance parameters. That's the very confusing thing that happens with || syntax for categorical factors on the left side of the double bar. $\endgroup$
    – amoeba
    Jan 31 '18 at 16:24
  • $\begingroup$ oops, that's right. I think afex may have a solution for this. $\endgroup$
    – Ben Bolker
    Jan 31 '18 at 17:24
  • $\begingroup$ By the way, is it something on the lme4 to-do list, or is not going to be implemented there? I find it quite confusing to be honest. $\endgroup$
    – amoeba
    Jan 31 '18 at 18:41
  • 2
    $\begingroup$ it's sufficiently difficult/architectural that it's pretty far down the to-do list. I don't find myself with much time to do more than surface and/or critical bug fixes (luckily, it's pretty stable). With Henrik Singmann's permission (which I'm guessing he would give) I could try to pull expand_re from here and incorporate it into lme4 ... I'd be extremely happy to consider a pull request :-) $\endgroup$
    – Ben Bolker
    Jan 31 '18 at 19:56

This site is temporarily in read only mode and not accepting new answers.

Not the answer you're looking for? Browse other questions tagged .