Should a grouping variable that is nested within a fixed effect be included as a random effect? 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?

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.
