{"id":3056,"date":"2023-05-08T04:26:35","date_gmt":"2023-05-08T04:26:35","guid":{"rendered":"http:\/\/optimumsportsperformance.com\/blog\/?p=3056"},"modified":"2023-05-08T12:28:52","modified_gmt":"2023-05-08T12:28:52","slug":"plotting-mixed-model-outputs","status":"publish","type":"post","link":"https:\/\/optimumsportsperformance.com\/blog\/plotting-mixed-model-outputs\/","title":{"rendered":"Plotting Mixed Model Outputs"},"content":{"rendered":"

This weekend I posted two new blog articles about building reports that contained both data tables and plots on the same canvas (see HERE<\/span><\/strong> and HERE<\/a><\/span><\/strong>). As a follow up, James Baker<\/a><\/span><\/strong> 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<\/a><\/span><\/strong> and HERE<\/a><\/span><\/strong>) and this got me thinking. Because he left it pretty wide open (“Do you have any guides on visualizing mixed models?”<\/em>) 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<\/strong>} 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.<\/span><\/p>\n

All of the code for this article is available on my GITHUB page<\/a><\/span><\/strong>.<\/p>\n

Loading Packages & Data<\/strong><\/span><\/p>\n

As always we begin by loading some of the packages we require and the data. In this case, we will use the sleepstudy <\/strong><\/em>dataset, which is freely available from the {lme4<\/strong>} package.<\/p>\n

\r\n## Load packages\r\nlibrary(tidyverse)\r\nlibrary(lme4)\r\nlibrary(lattice)\r\nlibrary(patchwork)\r\n\r\ntheme_set(theme_bw())\r\n\r\n## load data\r\ndat <- sleepstudy dat %>%\r\n  head()\r\n<\/pre>\n
Fit a mixed model<\/strong><\/span><\/p>\n
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.<\/p>\n
\r\n## Fit mixed model\r\nfit_lmer <- lmer(Reaction ~ Days + (1 + Days|Subject), data = dat)\r\nsummary(fit_lmer)\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Inspect the random effects<\/strong><\/span><\/p>\n
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).<\/p>\n
\r\n# look at the random effects\r\nrandom_effects <- ranef(fit_lmer) %>%\r\n  pluck(1) %>%\r\n  rownames_to_column() %>%\r\n  rename(Subject = rowname, Intercept = "(Intercept)") \r\n\r\nrandom_effects %>%\r\n  knitr::kable()\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Plotting the random effects<\/strong><\/span><\/p>\n
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.<\/p>\n
The {lattice<\/strong>} package allows us to create waterfall plots of the random effects for each subject with the dotplot()<\/strong> function.<\/p>\n
\r\n## plot random effects\r\ndotplot(ranef(fit_lmer))\r\n<\/pre>\n
\"\"<\/a><\/p>\n
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<\/strong>} so that we have more control over the styling.<\/p>\n
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.<\/p>\n
\r\n## Make one in ggplot2\r\nsubject_intercept_sd <- 24.7\r\nsubject_days_sd <- 5.92\r\n\r\nint_plt <- random_effects %>%\r\nmutate(Subject = as.factor(Subject)) %>%\r\nggplot(aes(x = Intercept, y = reorder(Subject, Intercept))) +\r\ngeom_errorbar(aes(xmin = Intercept - subject_intercept_sd,\r\nxmax = Intercept + subject_intercept_sd),\r\nwidth = 0,\r\nsize = 1) +\r\ngeom_point(size = 3,\r\nshape = 21,\r\ncolor = "black",\r\nfill = "white") +\r\ngeom_vline(xintercept = 0,\r\ncolor = "red",\r\nsize = 1,\r\nlinetype = "dashed") +\r\nscale_x_continuous(breaks = seq(-60, 60, 20)) +\r\nlabs(x = "Intercept",\r\ny = "Subject ID",\r\ntitle = "Random Intercepts")\r\n\r\nslope_plt <- random_effects %>%\r\nmutate(Subject = as.factor(Subject)) %>%\r\nggplot(aes(x = Days, y = reorder(Subject, Days))) +\r\ngeom_errorbar(aes(xmin = Days - subject_days_sd,\r\nxmax = Days + subject_days_sd),\r\nwidth = 0,\r\nsize = 1) +\r\ngeom_point(size = 3,\r\nshape = 21,\r\ncolor = "black",\r\nfill = "white") +\r\ngeom_vline(xintercept = 0,\r\ncolor = "red",\r\nsize = 1,\r\nlinetype = "dashed") +\r\nxlim(-60, 60) +\r\nlabs(x = "Slope",\r\ny = "Subject ID",\r\ntitle = "Random Slopes")\r\n\r\nslope_plt \/ int_plt\r\n<\/pre>\n
\"\"<\/a><\/p>\n
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.<\/p>\n
Plotting the model residuals<\/strong><\/span><\/p>\n
We can also plot the model residuals. Using the residual() <\/strong>function we can get the residuals directly from our mixed model and the plot()<\/strong> 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.<\/p>\n
\r\n## Plot Residual\r\nplot(fit_lmer)\r\nhist(resid(fit_lmer))\r\n\r\n<\/pre>\n
\"\"<\/a> \"\"<\/a><\/p>\n
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.<\/p>\n
\r\n## Plotting our own residual ~ fitted\r\nlmer_fitted <- predict(fit_lmer, newdata = dat, re.form = ~(1 + Days|Subject))\r\nlmer_resid <- dat$Reaction - lmer_fitted\r\n\r\nplot(x = lmer_fitted,\r\n     y = lmer_resid,\r\n     pch = 19,\r\n     main = "Resid ~ Fitted",\r\n     xlab = "Fitted",\r\n     ylab = "Residuals")\r\nabline(h = 0,\r\n       col = "red",\r\n       lwd = 3,\r\n       lty = 2)\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Plotting Predictions<\/strong><\/span><\/p>\n
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.<\/p>\n
 <\/p>\n
\r\n### Plotting the time series on new data\r\n# training set\r\ndat_train <- dat %>%\r\n  group_by(Subject) %>%\r\n  slice(head(row_number(), 6)) %>%\r\n  ungroup()\r\n\r\n# testing set\r\ndat_test <- dat %>%\r\n  group_by(Subject) %>%\r\n  slice(tail(row_number(), 4)) %>%\r\n  ungroup()\r\n\r\n## Fit mixed model\r\nfit_lmer2 <- lmer(Reaction ~ Days + (1 + Days|Subject), data = dat_train)\r\nsummary(fit_lmer2)\r\n\r\n# Predict on training set\r\ntrain_preds  <- merTools::predictInterval(fit_lmer2, newdata = dat_train, n.sims = 100, returnSims = TRUE, seed = 657, level = 0.9) %>%\r\n  as.data.frame()\r\n\r\ndat_train <- dat_train %>% bind_cols(train_preds)\r\n\r\ndat_train$group <- "train"\r\n\r\n# Predict on test set with 90% prediction intervals\r\ntest_preds  <- merTools::predictInterval(fit_lmer2, newdata = dat_test, n.sims = 100, returnSims = TRUE, seed = 657, level = 0.9) %>%\r\n  as.data.frame()\r\n\r\ndat_test <- dat_test %>% bind_cols(test_preds)\r\n\r\ndat_test$group <- "test"\r\n\r\n## Combine the data together\r\ncombined_dat <- bind_rows(dat_train, dat_test) %>%\r\n  arrange(Subject, Days)\r\n\r\n## Plot the time series of predictions and observed data\r\ncombined_dat %>%\r\nmutate(group = factor(group, levels = c("train", "test"))) %>%\r\nggplot(aes(x = Days, y = Reaction)) +\r\n  geom_ribbon(aes(ymin = lwr,\r\n                  ymax = upr),\r\n              fill = "light grey",\r\n              alpha = 0.8) +\r\n  geom_line(aes(y = fit),\r\n            col = "red",\r\n            size = 1) +\r\n  geom_point(aes(fill = group),\r\n             size = 3,\r\n             shape = 21) +\r\n  geom_line() +\r\n  facet_wrap(~Subject) +\r\n  theme(strip.background = element_rect(fill = "black"),\r\n        strip.text = element_text(face = "bold", color = "white"),\r\n        legend.position = "top") +\r\n  labs(x = "Days",\r\n       y = "Reaction Time",\r\n       title = "Reaction Time based on Days of Sleep Deprivation")\r\n<\/pre>\n
\"\"<\/a><\/p>\n
Wrapping Up<\/strong><\/span><\/p>\n
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!<\/p>\n
The complete code for this article is available on my GITHUB page<\/a><\/span><\/strong>.<\/p>\n","protected":false},"excerpt":{"rendered":"
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 […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[47,43],"tags":[],"class_list":["post-3056","post","type-post","status-publish","format-standard","hentry","category-model-building-in-r","category-sports-analytics"],"_links":{"self":[{"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/posts\/3056","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/comments?post=3056"}],"version-history":[{"count":2,"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/posts\/3056\/revisions"}],"predecessor-version":[{"id":3067,"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/posts\/3056\/revisions\/3067"}],"wp:attachment":[{"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/media?parent=3056"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/categories?post=3056"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/optimumsportsperformance.com\/blog\/wp-json\/wp\/v2\/tags?post=3056"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}