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?

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/).
