# How do you know the number of random effects in a mixed effects model?

I am trying to fit a random slope model in R and my code is as follows:

lmer(data=ds, Outcome ~ treatment + (0 + treatment|ID))

I get the following error message when I try running this code:

Error: number of observations (=2035) <= number of random effects (=2035) for term (0 + treatment | ID); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable

I have 407 unique IDs and 5 treatments, so the number of observations is 407*5 = 2035. However, I don't understand why I have 2035 random effects. Being a random slope only model, I would expect that I am estimating one random slope per unique ID and so I thought there would only be 407 random effects in my model.

I also tried fitting a random slope and intercept model as follows:

lmer(data=ds, Outcome ~ treatment + (1 + treatment|ID))

However, I still got the same error message as above, that I have 2035 observations and the number of random effects = 2035. I get that the number of observations should be more than the number of random effects being estimated, but I still don't understand why there are 2035 random effects in this model either. I would expect that since I am estimating a random slope and intercept for each unique ID (2 random effects for each ID), I would have 407*2=814 random effects. Clearly, there is something I am missing here about how to calculate the number of random effects for these two models. Any help understanding this would be highly appreciated. Thanks.

• The issue is that treatment is a categorical variable. This means that, due to dummy encoding, you will have several (in your case five) coefficients per ID. Your random effect can't include treatment because you lack the repeated measures. There are no degrees of freedom left for the residual variance. The correct random effect in your design should be (1 | ID). Commented Sep 21, 2020 at 10:40
• Yes, you are right @Roland. I don't have repeated measures by treatment. Thanks! Commented Sep 21, 2020 at 20:04

This is expected behavious whenever you try to fit a model with random slopes where the variable for the random slopes is categorical and there is only one observation per treatment/group combination.

This is because the levels of a categorical variable are represented by dummy variables - essentially they are treated as different variables. So in your case, when you fit random slopes only you are asking the software to estimate 5 random slopes for each group. When you fit random intercepts and random slopes there will be 407 random intercepts, but only 4 random slopes for each group (since one level will be treated as a reference group and included in the intercept), so either way you will have 5 x 407 random effects.

The only way to solve this is by either coding the variable as numeric, if that is plausible in your study/data, or not fitting random slopes, or having more than 1 observation per treatment per group.

It may be illustrative to see this with a toy dataset:

> set.seed(1)
> dt <- expand.grid(G = LETTERS[1:4], a = LETTERS[1:2])
> dt\$Y = rnorm(nrow(dt))
> dt
G a          Y
1 A A -0.6264538
2 B A  0.1836433
3 C A -0.8356286
4 D A  1.5952808
5 A B  0.3295078
6 B B -0.8204684
7 C B  0.4874291
8 D B  0.7383247

Now we fit the models, both of which will not run for the reasons explained above.

> lmer(Y ~ a + (0 + a | G), data = dt) %>% summary()
Error: number of observations (=8) <= number of random effects (=8) for term (0 + a | G); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable
> lmer(Y ~ a + (1 + a | G), data = dt) %>% summary()
Error: number of observations (=8) <= number of random effects (=8) for term (1 + a | G); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable

But now we add just 1 extra row to the dataset, and they run:

> (dt <- rbind(dt, dt[1, ]))
G a          Y
1 A A -0.6264538
2 B A  0.1836433
3 C A -0.8356286
4 D A  1.5952808
5 A B  0.3295078
6 B B -0.8204684
7 C B  0.4874291
8 D B  0.7383247
9 A A -0.6264538

> lmer(Y ~ a + (0 + a | G), data = dt) %>% summary()

Random effects:
Groups   Name Variance  Std.Dev.  Corr
G        aA   1.451e+00 1.205e+00
aB   3.224e-01 5.678e-01 -0.04
Residual      4.239e-15 6.511e-08

> lmer(Y ~ a + (1 + a | G), data = dt) %>% summary()

Random effects:
Groups   Name        Variance  Std.Dev.  Corr
G        (Intercept) 9.776e-01 9.887e-01
aB          1.222e+00 1.105e+00 -0.81
Residual             1.159e-14 1.077e-07
Number of obs: 9, groups:  G, 4

In the model with random slopes only we have 2 random slopes in 4 groups (8 random effects), and in the model with both random intercepts and random slopes we have 4 random intercepts and 4 random slopes.

• Thanks @Robert Long for the detailed explanation! It's definitely more clear now. Not fitting random slopes would probably be more ideal for this problem. Thanks. Commented Sep 21, 2020 at 20:13