# Mixed Effects

This vignette replicates some of the analyses in this excellent blog post by Solomon Kurz: Use `emmeans()` to include 95% CIs around your `lme4`-based fitted lines

Load libraries and fit two models of chick weights:

``````library(lme4)
library(tidyverse)
library(patchwork)
library(marginaleffects)

## unconditional linear growth model
fit1 <- lmer(
weight ~ 1 + Time + (1 + Time | Chick),
data = ChickWeight)

fit2 <- lmer(
weight ~ 1 + Time + I(Time^2) + Diet + Time:Diet + I(Time^2):Diet + (1 + Time + I(Time^2) | Chick),
data = ChickWeight)``````

### Unit-level predictions

Predict weight of each chick over time:

``````pred1 <- predictions(fit1,
newdata = datagrid(Chick = ChickWeight\$Chick,
Time = 0:21))

p1 <- ggplot(pred1, aes(Time, estimate, level = Chick)) +
geom_line() +
labs(y = "Predicted weight", x = "Time", title = "Linear growth model")

pred2 <- predictions(fit2,
newdata = datagrid(Chick = ChickWeight\$Chick,
Time = 0:21))

p2 <- ggplot(pred2, aes(Time, estimate, level = Chick)) +
geom_line() +
labs(y = "Predicted weight", x = "Time", title = "Quadratic growth model")

p1 + p2``````

Predictions for each chick, in the 4 counterfactual worlds with different values for the `Diet` variable:

``````pred <- predictions(fit2)

ggplot(pred, aes(Time, estimate, level = Chick)) +
geom_line() +
ylab("Predicted Weight") +
facet_wrap(~ Diet, labeller = label_both)``````

### Population-level predictions

To make population-level predictions, we set the `Chick` variable to `NA`, and set `re.form=NA`. This last argument is offered by the `lme4::predict` function which is used behind the scenes to compute predictions:

``````pred <- predictions(
fit2,
newdata = datagrid(Chick = NA,
Diet = 1:4,
Time = 0:21),
re.form = NA)

ggplot(pred, aes(x = Time, y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_ribbon(alpha = .1, fill = "red") +
geom_line() +
facet_wrap(~ Diet, labeller = label_both) +
labs(title = "Population-level trajectories")``````