This weekend I posted two new blog articles about building reports that contained both data tables and plots on the same canvas (see HERE and HERE). As a follow up, James Baker asked if I could do some plotting of mixed model outputs. That got me thinking, I’ve done a few blog tutorials on mixed models (see HERE and HERE) and this got me thinking. Because he left it pretty wide open (“Do you have any guides on visualizing mixed models?”) I was trying to think about what aspects of the mixed models he’d like to visualize. R makes it relatively easy to plot random effects using the {lattice} package, but I figured we could go a little deeper and customize some of our own plots of the random effects as well as show how we might plot future predictions from a mixed model.
All of the code for this article is available on my GITHUB page.
Loading Packages & Data
As always we begin by loading some of the packages we require and the data. In this case, we will use the sleepstudy dataset, which is freely available from the {lme4} package.
## Load packages library(tidyverse) library(lme4) library(lattice) library(patchwork) theme_set(theme_bw()) ## load data dat <- sleepstudy dat %>% head()
Fit a mixed model
We will fit a mixed model that sets the dependent variable as Reaction time and the fixed effect as days of sleep deprivation. We will also allow both the intercept and slope to vary randomly by nesting the individual SubjectID within each Day of sleep deprivation.
## Fit mixed model fit_lmer <- lmer(Reaction ~ Days + (1 + Days|Subject), data = dat) summary(fit_lmer)
Inspect the random effects
We can see in the model output above that we have a random effect standard deviation for the Intercept (24.84) and for the slope, Days (5.92). We can extract out the random effect intercept and slope for each subject with the code below. This tells us how much each subject’s slope and intercept vary from the population fixed effects (251.4 and 10.5 for the intercept and slope, respectively).
# look at the random effects random_effects <- ranef(fit_lmer) %>% pluck(1) %>% rownames_to_column() %>% rename(Subject = rowname, Intercept = "(Intercept)") random_effects %>% knitr::kable()
Plotting the random effects
Aside from looking at a table of numbers, which can sometimes be difficult to draw conclusions from (especially if there are a large number of subjects) we can plot the data and make some observational inference.
The {lattice} package allows us to create waterfall plots of the random effects for each subject with the dotplot() function.
## plot random effects dotplot(ranef(fit_lmer))
That’s a pretty nice plot and easy to obtain with just a single line of code. But, we might want to create our own plot using {ggplot2} so that we have more control over the styling.
I’ll store the standard deviation of the random slope and intercept, from the model read out above, in their own element. Then, I’ll use the random effects table we made above, which contains the intercept and slope of each subject, to plot them and add the standard deviation to them as error bars.
## Make one in ggplot2 subject_intercept_sd <- 24.7 subject_days_sd <- 5.92 int_plt <- random_effects %>% mutate(Subject = as.factor(Subject)) %>% ggplot(aes(x = Intercept, y = reorder(Subject, Intercept))) + geom_errorbar(aes(xmin = Intercept - subject_intercept_sd, xmax = Intercept + subject_intercept_sd), width = 0, size = 1) + geom_point(size = 3, shape = 21, color = "black", fill = "white") + geom_vline(xintercept = 0, color = "red", size = 1, linetype = "dashed") + scale_x_continuous(breaks = seq(-60, 60, 20)) + labs(x = "Intercept", y = "Subject ID", title = "Random Intercepts") slope_plt <- random_effects %>% mutate(Subject = as.factor(Subject)) %>% ggplot(aes(x = Days, y = reorder(Subject, Days))) + geom_errorbar(aes(xmin = Days - subject_days_sd, xmax = Days + subject_days_sd), width = 0, size = 1) + geom_point(size = 3, shape = 21, color = "black", fill = "white") + geom_vline(xintercept = 0, color = "red", size = 1, linetype = "dashed") + xlim(-60, 60) + labs(x = "Slope", y = "Subject ID", title = "Random Slopes") slope_plt / int_plt
We get the same plot but now we have more control. We can color the dot specific subjects, or only choose to display specific subjects, or flip the x- and y-axes, etc.
Plotting the model residuals
We can also plot the model residuals. Using the residual() function we can get the residuals directly from our mixed model and the plot() function with automatically plot the Residual and Fitted values. These types of plots are useful for exploring assumptions such as normality of the residuals and homoscedasticity.
## Plot Residual plot(fit_lmer) hist(resid(fit_lmer))
As above, perhaps we want to have more control over the bottom plot, so that we can style it however we’d like. We can extract the fitted values and residuals and build our own plot using base R.
## Plotting our own residual ~ fitted
lmer_fitted <- predict(fit_lmer, newdata = dat, re.form = ~(1 + Days|Subject))
lmer_resid <- dat$Reaction - lmer_fitted
plot(x = lmer_fitted,
y = lmer_resid,
pch = 19,
main = "Resid ~ Fitted",
xlab = "Fitted",
ylab = "Residuals")
abline(h = 0,
col = "red",
lwd = 3,
lty = 2)
Plotting Predictions
The final plot I’ll build are the predictions of Reaction time as Days of sleep deprivation increase. This is time series data, so I’m going to extract the first 6 days of sleep deprivation for each subject and build the model using that data. Then, make predictions on the next 4 days of sleep deprivation for each subject and get both a predicted point estimate and 90% prediction interval. In this way, we can observe the next 4 days of sleep deprivation for each subject and see how far outside of what we would expect (from our mixed model predictions) those values fall.
### Plotting the time series on new data
# training set
dat_train <- dat %>%
group_by(Subject) %>%
slice(head(row_number(), 6)) %>%
ungroup()
# testing set
dat_test <- dat %>%
group_by(Subject) %>%
slice(tail(row_number(), 4)) %>%
ungroup()
## Fit mixed model
fit_lmer2 <- lmer(Reaction ~ Days + (1 + Days|Subject), data = dat_train)
summary(fit_lmer2)
# Predict on training set
train_preds <- merTools::predictInterval(fit_lmer2, newdata = dat_train, n.sims = 100, returnSims = TRUE, seed = 657, level = 0.9) %>%
as.data.frame()
dat_train <- dat_train %>% bind_cols(train_preds)
dat_train$group <- "train"
# Predict on test set with 90% prediction intervals
test_preds <- merTools::predictInterval(fit_lmer2, newdata = dat_test, n.sims = 100, returnSims = TRUE, seed = 657, level = 0.9) %>%
as.data.frame()
dat_test <- dat_test %>% bind_cols(test_preds)
dat_test$group <- "test"
## Combine the data together
combined_dat <- bind_rows(dat_train, dat_test) %>%
arrange(Subject, Days)
## Plot the time series of predictions and observed data
combined_dat %>%
mutate(group = factor(group, levels = c("train", "test"))) %>%
ggplot(aes(x = Days, y = Reaction)) +
geom_ribbon(aes(ymin = lwr,
ymax = upr),
fill = "light grey",
alpha = 0.8) +
geom_line(aes(y = fit),
col = "red",
size = 1) +
geom_point(aes(fill = group),
size = 3,
shape = 21) +
geom_line() +
facet_wrap(~Subject) +
theme(strip.background = element_rect(fill = "black"),
strip.text = element_text(face = "bold", color = "white"),
legend.position = "top") +
labs(x = "Days",
y = "Reaction Time",
title = "Reaction Time based on Days of Sleep Deprivation")
Wrapping Up
Above are a few different plot options we have with mixed model outputs. I’m not sure what James was after or what he had in mind because he left the question very wide open. Hopefully this article provides some useful ideas for your own mixed model plotting. If there are other things you are hoping to see or have other ideas of things to plot from the mixed model output, feel free to reach out!
The complete code for this article is available on my GITHUB page.