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?

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}
$$
