Should a grouping variable that is nested within a fixed effect be included as a random effect?

Question

I am wondering whether a grouping variable, which is nested within a fixed effect, should be included as a random effect in a linear mixed effects model.

I will exemplify this with simulated data in R. Say we have a fixed effect x that is an integer, ranging from 4 to 14. We have an outcome variable y that has a weak positive association with x. In addition, our data can be grouped by id, where observations with the same id have the same value for x. There can be multiple values of id which have the same value for x.

library(dplyr)
library(ggplot2)
library(lmerTest)

set.seed(261022)
data <- tibble(x = rep(c(4, 14, sample(x = 5:13, size = 21, replace = T)), 
                       each = 12),
               id = factor(rep(1:23, each = 12))) %>%
  group_by(id) %>%
  mutate(y = rnorm(n = n(), mean = (ifelse(x > 4 & x < 14,
                                           x - rbinom(1,1,0.5) * 0.75 * x,
                                           x)), 
                   sd = 1))

We can visualise this data:

ggplot(data, aes(x = x, y = y)) +
  geom_point(aes(colour = id)) +
  theme_bw()

I am wondering which of the following models would be more appropriate for this data:

lin_mod <- lm(y ~ x, data = data)
mixed_mod <- lmer(y ~ x + (1|id), data = data)

On the one hand, I believe a random effect should be included for id, since we have multiple observations per value of id. However, this results in varying intercepts across id, and by extension, parallel lines across id. Having a separate line for each id seems nonsensical, since an id can by definition only have one value of x.

Would either of these models be appropriate, or is there an alternative way to account for the grouping of the data by id?

Answer 1

With this data structure, the model

mixed_mod <- lmer(y ~ x + (1|id), data = data)

can be thought of as:

• An overall regression line defined by y ~ 1 + x
• Clusters within each id that deviate from this line. This variation is captured by the random intercept terms, (1|id)
• Data points within each cluster that deviate from the cluster mean. This variation is captured by the residual variance.

You're correct that there aren't really separate lines for each id: this would only be the case if x varied within id groups.

On the other hand, the model

lin_mod <- lm(y ~ x, data = data)

has only the overall line and the residual variance, but doesn't account for the clustering by id. This means that all deviations from the regression line are included in the residual variance.

Answer 2

For ease of writing, I assume here id refers to participants, i.e., each level of id = person.

This may not be the best formulated reply, but the way I would approach this: If I understand correctly, your example has a "level 1" dependent variable (i.e., you have several values of y for each participant), but only a "level 2" independent variable (you have one value of x for each participant). This is common in, say, experience sampling studies when you predict a momentary feeling (a level 1 variable, many values per participant) from, say, a personality trait (a level 2 variable, one value per participant).

In this case you cannot estimate a participant-specific slope between x and y (that you would get from (x|id) if you had a level 1 predictor), but you can estimate a participant-specific random intercept with (1|id).

When you use

mixed_mod <- lmer(y ~ x + (1|id), data = data)

and plot the results, you get, visually, a separate slope for each id (all slopes being parallel) - but this slope is in fact just the overall slope of y ~ x positioned at different levels for each participant as a function of their personal intercept + overall intercept. The usefulness of this is that you see how participants vary in terms of their average levels of outcome, given the grand average level of x. But you could say, as you do, that separate lines in this case are somewhat nonsensical. It could be said it would be more logical to plot participant-specific intercepts as, say, points and not lines and only plot the overall y ~ x effect as a line.

As for your questions: the mixed_mod code you post would be correct for modelling this type of data.

An alternative way would be to fit a single-level regression of y ~ x with clustered standard errors (r does not have a built-in way to do that, see here: https://www.r-bloggers.com/2021/05/clustered-standard-errors-with-r/).

Answer 3

To complement the answers by @Eoin and @Sointu (+1 to both), yes, there is an alternative way to account for the structure in the data due to having multiple observations from the same participant. If the outcome variable Y is normally distributed and only the population-level (not the subject-level) effects of the predictors are of interest, then you can also model the data with Generalized Least Squares (GLS).

Here is how to do it in R using nlme::gls. In the example, observations with the same id are equally correlated; this corresponds to a correlation structure with compound symmetry.

# Fit a linear mixed model with random intercepts
library("lme4")
m1 <- lmer(
  y ~ x + (1 | id),
  data = data
)

# Fit a linear model with compound symmetry structure
library("nlme")
m2 <- gls(
  y ~ x,
  correlation = corCompSymm(form = ~ 1 | id),
  data = data
)

Random intercepts correspond to compound symmetry correlation (see below); so the estimated population-level effects (i.e. the fixed effects) are the same in this case.

library("broom.mixed")
tidy(m1, "fixed")
#> # A tibble: 2 × 5
#>   effect term        estimate std.error statistic
#>   <chr>  <chr>          <dbl>     <dbl>     <dbl>
#> 1 fixed  (Intercept)    1.29      2.27      0.570
#> 2 fixed  x              0.453     0.252     1.80
tidy(m2)
#> # A tibble: 2 × 5
#>   term        estimate std.error statistic p.value
#>   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
#> 1 (Intercept)    1.29      2.27      0.570  0.569 
#> 2 x              0.453     0.252     1.80   0.0729

Appendix

Take the varying-intercepts model $Y_{ij} = \beta_0 + \beta_1x_{ij} + u_i + e_{ij}$. Then for any two observations $Y_{ij}$ and $Y_{ik}$ with the same id $i$, we have:

$$ \begin{aligned} \operatorname{Cov}\left(Y_{ij}, Y_{ik}\right) &= \operatorname{Cov}\left(\beta_0 + \beta_1x_{ij} + u_i + e_{ij}, \beta_0 + \beta_1x_{ik} + u_i + e_{ik}\right) \\ &= \operatorname{Cov}\left(u_i, u_i\right) + \operatorname{Cov}\left(e_{ij}, e_{ik}\right) \\ &= \operatorname{Var}\left(u_i\right) \end{aligned} $$