# Linear mixed model with heteroscedastic residuals
mod <- galamm(
formula = y ~ x + (1 | id),
weights = ~ (1 | item),
data = hsced
)
# Extract fitted values and plot against x
plot(hsced$x, fitted(mod))
Run the code above in your browser using DataLab