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