Plotting partial pooling in mixedeffects models
In this post, I demonstrate a few techniques for plotting information from a relatively simple mixedeffects model fit in R. These plots can help us develop intuitions about what these models are doing and what “partial pooling” means.
The sleepstudy
dataset
For these examples, I’m going to use the
sleepstudy
dataset from the lme4 package. The
outcome measure is reaction time, the predictor measure is days of sleep
deprivation, and these measurements are nested within participants—we
have 10 observations per participant. I am also going to add two fake
participants with incomplete data to illustrate partial pooling.
library(lme4)
#> Loading required package: Matrix
library(dplyr)
library(tibble)
# Convert to tibble for better printing. Convert factors to strings
sleepstudy < sleepstudy %>%
as_tibble() %>%
mutate(Subject = as.character(Subject))
# Add two fake participants
df_sleep < bind_rows(
sleepstudy,
tibble(Reaction = c(286, 288), Days = 0:1, Subject = "374"),
tibble(Reaction = 245, Days = 0, Subject = "373"))
df_sleep
#> # A tibble: 183 × 3
#> Reaction Days Subject
#> <dbl> <dbl> <chr>
#> 1 250. 0 308
#> 2 259. 1 308
#> 3 251. 2 308
#> 4 321. 3 308
#> 5 357. 4 308
#> 6 415. 5 308
#> 7 382. 6 308
#> 8 290. 7 308
#> 9 431. 8 308
#> 10 466. 9 308
#> # … with 173 more rows
We can visualize all the data in ggplot2 by using
facet_wrap()
to create subplots for each
participant and stat_smooth()
to create a
regression line in each subplot.
library(ggplot2)
xlab < "Days of sleep deprivation"
ylab < "Average reaction time (ms)"
ggplot(df_sleep) +
aes(x = Days, y = Reaction) +
stat_smooth(method = "lm", se = FALSE) +
# Put the points on top of lines
geom_point() +
facet_wrap("Subject") +
labs(x = xlab, y = ylab) +
# We also need to help the xaxis, so it doesn't
# create gridlines/ticks on 2.5 days
scale_x_continuous(breaks = 0:4 * 2)
#> `geom_smooth()` using formula 'y ~ x'
By the way, ggplot2 doesn’t draw the regression lines outside of the
range of the data unless we set fullrange = TRUE
. That’s a helpful
feature for 374!
Update: Douglas Bates did it first. Someone sent me a link to a
slide deck by Douglas
Bates,
lead author of the lme4 package, where he has some plots just like the
ones I demo in this post. He uses the sleepstudy
dataset too—it’s
his R package and his teaching dataset, after all—so the similarities
are uncanny but accidental. Origin of this post: I was asked on
twitter
how to make a facet plot of a mixed effects model, wrote up a quick
demo using the
convenient sleepstudy
dataset, and then fleshed that demo into a
tutorial. By using his teaching dataset to illustrate some partial
pooling concepts, I ended up recreating some of his work on accident.
😬 [Sept. 14, 2017]
Complete pooling and no pooling models
Each one of these panels plotted above shows an independently estimated regression line. This approach to fitting a separate line for each participant is sometimes called the no pooling model because none of the information from different participants is combined or pooled together.
We fit a separate line for each cluster of data, unaware that any of the
other participants exist. The lmList()
function in
lme4 automates this process.
df_no_pooling < lmList(Reaction ~ Days  Subject, df_sleep) %>%
coef() %>%
# Subject IDs are stored as rownames. Make them an explicit column
rownames_to_column("Subject") %>%
rename(Intercept = `(Intercept)`, Slope_Days = Days) %>%
add_column(Model = "No pooling") %>%
# Remove the participant who only had one datapoint
filter(Subject != "373")
head(df_no_pooling)
#> Subject Intercept Slope_Days Model
#> 1 308 244.1927 21.764702 No pooling
#> 2 309 205.0549 2.261785 No pooling
#> 3 310 203.4842 6.114899 No pooling
#> 4 330 289.6851 3.008073 No pooling
#> 5 331 285.7390 5.266019 No pooling
#> 6 332 264.2516 9.566768 No pooling
In contrast, we might consider a complete pooling model where all the information from the participants is combined together. We fit a single line for the combined data set, unaware that the data came from different participants.
# Fit a model on all the data pooled together
m_pooled < lm(Reaction ~ Days, df_sleep)
# Repeat the intercept and slope terms for each participant
df_pooled < tibble(
Model = "Complete pooling",
Subject = unique(df_sleep$Subject),
Intercept = coef(m_pooled)[1],
Slope_Days = coef(m_pooled)[2]
)
head(df_pooled)
#> # A tibble: 6 × 4
#> Model Subject Intercept Slope_Days
#> <chr> <chr> <dbl> <dbl>
#> 1 Complete pooling 308 252. 10.3
#> 2 Complete pooling 309 252. 10.3
#> 3 Complete pooling 310 252. 10.3
#> 4 Complete pooling 330 252. 10.3
#> 5 Complete pooling 331 252. 10.3
#> 6 Complete pooling 332 252. 10.3
We can compare these two approaches. Instead of calculating the
regression lines with stat_smooth()
, we can use
geom_abline()
to draw the lines from our dataframe
of intercept and slope parameters.
# Join the raw data so we can use plot the points and the lines.
df_models < bind_rows(df_pooled, df_no_pooling) %>%
left_join(df_sleep, by = "Subject")
p_model_comparison < ggplot(df_models) +
aes(x = Days, y = Reaction) +
# Set the color mapping in this layer so the points don't get a color
geom_abline(
aes(intercept = Intercept, slope = Slope_Days, color = Model),
size = .75
) +
geom_point() +
facet_wrap("Subject") +
labs(x = xlab, y = ylab) +
scale_x_continuous(breaks = 0:4 * 2) +
# Fix the color palette
scale_color_brewer(palette = "Dark2") +
theme(legend.position = "top", legend.justification = "left")
p_model_comparison
If we stare at this plot, a few things become apparent. The complete pooling model estimates a single line, and we see that same line drawn on every facet. One advantage is that the model can make a guess about the line for 373 who only has one observation. That model looks pretty terrible elsewhere—309, 310, etc.—because nobody is perfectly average. In contrast, the no pooling model can follow the data, fitting the sharp trend upwards in 308 and even capturing the negative slope in 335.
(Here’s a fun question: Which approach has the better guess for 374’s line?)
The no pooling model cannot make a guess about 373. In Statistical Rethinking, McElreath says these models have amnesia 😯:
Many statistical models also have anterograde amnesia. As the models move from one cluster—individual, group, location—in the data to another, estimating parameters for each cluster, they forget everything about the previous clusters. They behave this way, because the assumptions force them to. Any of the models from previous chapters that used dummy variables to handle categories are programmed for amnesia. These models implicitly assume that nothing learned about any one category informs estimates for the other categories—the parameters are independent of one another and learn from completely separate portions of the data. This would be like forgetting you had ever been in a café, each time you go to a new café. Cafés do differ, but they are also alike.
Once the no pooling model draws the line for 372, and it completely forgets everything it has seen and moves on to 373. It has to skip 373 because it cannot estimate a line from a single point, and it moves on.
Improving estimates with a mixedeffects model
We can do better with mixedeffects models. In these models, we pool information from all the lines together to improve our estimates of each individual line. This approach is sometimes called partial pooling. In particular, after seeing the 18 trend lines for the participants with complete data, we can make an informed guess about the trend lines for the two participants with incomplete data.
We can fit a classical mixedeffects model with the lme4 package:
m < lmer(Reaction ~ 1 + Days + (1 + Days  Subject), df_sleep)
arm::display(m)
#> lmer(formula = Reaction ~ 1 + Days + (1 + Days  Subject), data = df_sleep)
#> coef.est coef.se
#> (Intercept) 252.54 6.43
#> Days 10.45 1.54
#>
#> Error terms:
#> Groups Name Std.Dev. Corr
#> Subject (Intercept) 24.14
#> Days 5.92 0.07
#> Residual 25.48
#> 
#> number of obs: 183, groups: Subject, 20
#> AIC = 1783.4, DIC = 1787.8
#> deviance = 1779.6
The first two coef.est
items are the “fixed effects” estimates; they
reflect the average intercept and slope parameters. For this example,
the values are practically the same as the completepooling estimates.
This model assumes that each participant’s individual intercept and
slope parameters are deviations from this average, and these random
deviations drawn from a distribution of possible intercept and slope
parameters. These are “randomly varying” or “random effects”. The
information in the Error terms
area describes the distribution of the
effects. Because we have both fixed and random effects, we have a
“mixedeffects” model. Hence the name.
To visualize these estimates, we extract each participant’s intercept
and slope using coef()
.
# Make a dataframe with the fitted effects
df_partial_pooling < coef(m)[["Subject"]] %>%
rownames_to_column("Subject") %>%
as_tibble() %>%
rename(Intercept = `(Intercept)`, Slope_Days = Days) %>%
add_column(Model = "Partial pooling")
head(df_partial_pooling)
#> # A tibble: 6 × 4
#> Subject Intercept Slope_Days Model
#> <chr> <dbl> <dbl> <chr>
#> 1 308 254. 19.6 Partial pooling
#> 2 309 212. 1.73 Partial pooling
#> 3 310 213. 4.91 Partial pooling
#> 4 330 275. 5.64 Partial pooling
#> 5 331 274. 7.39 Partial pooling
#> 6 332 261. 10.2 Partial pooling
Update the previous plot with a dataframe of all three models’ estimates.
df_models < bind_rows(df_pooled, df_no_pooling, df_partial_pooling) %>%
left_join(df_sleep, by = "Subject")
# Replace the dataset of the last plot
p_model_comparison %+% df_models
Most of the time, the no pooling and partial pooling lines are on top of each other. But when the two differ, it’s because the partial pooling model’s line is pulled slightly towards the completepooling line.
We can appreciate the differences by zooming in on some participants.
df_zoom < df_models %>%
filter(Subject %in% c("335", "350", "373", "374"))
p_model_comparison %+% df_zoom
The negative line for 335 from the no pooling model gets a flatter slope in the partial pooling model. The model knows that negative trends are rather unlikely, so the it hedges its bets and pulls that line towards the group average. Something similar happens with 350 where a sharp slope is slightly attenuated. For the participants with incomplete data, the partial pooling model is much more like the complete pooling model. The complete pooling and the partial pooling lines are basically parallel—i.e, they have the same slope. That’s a reasonable guess given so little information.
It’s shrinkage
The partial pooling model pulls more extreme estimates towards an overall average. We can visualize this effect by plotting a scatterplot of intercept and slope parameters from each model and connecting estimates for the same participant.
# Also visualize the point for the fixed effects
df_fixef < tibble(
Model = "Partial pooling (average)",
Intercept = fixef(m)[1],
Slope_Days = fixef(m)[2]
)
# Complete pooling / fixed effects are center of gravity in the plot
df_gravity < df_pooled %>%
distinct(Model, Intercept, Slope_Days) %>%
bind_rows(df_fixef)
df_gravity
#> # A tibble: 2 × 3
#> Model Intercept Slope_Days
#> <chr> <dbl> <dbl>
#> 1 Complete pooling 252. 10.3
#> 2 Partial pooling (average) 253. 10.5
df_pulled < bind_rows(df_no_pooling, df_partial_pooling)
ggplot(df_pulled) +
aes(x = Intercept, y = Slope_Days, color = Model, shape = Model) +
geom_point(size = 2) +
geom_point(
data = df_gravity,
size = 5,
# Prevent size5 point from showing in legend keys
show.legend = FALSE
) +
# Draw an arrow connecting the observations between models
geom_path(
aes(group = Subject, color = NULL),
arrow = arrow(length = unit(.02, "npc")),
show.legend = FALSE
) +
# Use ggrepel to jitter the labels away from the points
ggrepel::geom_text_repel(
aes(label = Subject, color = NULL),
data = df_no_pooling,
show.legend = FALSE
) +
# Don't forget 373
ggrepel::geom_text_repel(
aes(label = Subject, color = NULL),
data = filter(df_partial_pooling, Subject == "373"),
show.legend = FALSE
) +
theme(
legend.position = "bottom",
legend.justification = "right"
) +
ggtitle("Pooling of regression parameters") +
xlab("Intercept estimate") +
ylab("Slope estimate") +
scale_shape_manual(values = c(15:18)) +
scale_color_brewer(palette = "Dark2")
The average intercept and slope act like a center of gravity, pulling values parameter estimates towards it. Hmm, maybe gravity is not quite the right analogy, because the pull is greater for more extreme values. The lines near that center point are very short; they get adjusted very little. The lines in general get longer as we move away from the complete pooling estimate. The fewer the observations in a cluster (here, participants), the more information is borrowed from other clusters, and the greater the pull towards the average estimate. Participant 373 had one observation, so their slope estimate is the average. Likewise, 374 had only two observations, so they get pulled the farthest and receive a slope estimate near the overall average.
This effect is sometimes called shrinkage, because more extreme values shrinkage are pulled towards a more reasonable, more average value. In the lme4 book, Douglas Bates provides an alternative to shrinkage:
The term “shrinkage” may have negative connotations. John Tukey preferred to refer to the process as the estimates for individual subjects “borrowing strength” from each other. This is a fundamental difference in the models underlying mixedeffects models versus strictly fixed effects models. In a mixedeffects model we assume that the levels of a grouping factor are a selection from a population and, as a result, can be expected to share characteristics to some degree. Consequently, the predictions from a mixedeffects model are attenuated relative to those from strictly fixedeffects models.
Shrinkage, borrowing strength 💪 … Another term would also be regularization if we think about how the model avoids overfitting by the taming extreme estimates.
This feature is why I use mixed effects models in my work. If I have a speech discrimination experiment and I want to describe a child’s speech perception ability, I am going to use the partially pooled, shrunken, strengthborrowing, regularized, modelderived estimate of their ability, because it uses more information. It’s that simple to me. Other disciplines might highlight other reasons to use these models, but for me, it’s partially pooling information that’s the most attractive feature.
A topographic map of parameters
For the next visualization, we are going to visualize the distribution of randomly varying effects. Honestly, I am partly including it just so that I have a working ggplot2 version of how to make this plot online. It’s not a routine visualization, but it reveals a little more about where estimates are being pulled towards.
I already remarked that the model estimates a distribution of intercept and slope effects. We know where the center of that distribution is: It’s the fixed effects estimate, the center of gravity in the last plot. What the model also needs to estimate is the variability/spread of values around that center. Also, intercepts and slopes might be correlated: Maybe the effect of an additional day on reaction time is diminished for participants who are slower to respond in general. So, the model also estimates the correlation of those effects too.
Imagine that the last plot is a landscape, and fixed effects point is the peak of a hill. What were are going to do is draw a topographic map with contour lines to show different elevation regions on that hill.
First, we need to extract the covariance matrix estimated by the model.
# Extract the matrix
cov_mat < VarCorr(m)[["Subject"]]
# Strip off some details so that just the useful part is printed
attr(cov_mat, "stddev") < NULL
attr(cov_mat, "correlation") < NULL
cov_mat
#> (Intercept) Days
#> (Intercept) 582.717345 9.897673
#> Days 9.897673 35.033088
The ellipse()
function takes a covariance matrix, a
center value, and quantile/confidence level and returns the points from
an oval around the center at the given confidence level. We create five
ellipses for different quantile levels.
library(ellipse)
#>
#> Attaching package: 'ellipse'
#> The following object is masked from 'package:graphics':
#>
#> pairs
# Helper function to make a dataframe of ellipse points that
# includes the level as a column
make_ellipse < function(cov_mat, center, level) {
ellipse(cov_mat, centre = center, level = level) %>%
as.data.frame() %>%
add_column(level = level) %>%
as_tibble()
}
center < fixef(m)
levels < c(.1, .3, .5, .7, .9)
# Create an ellipse dataframe for each of the levels defined
# above and combine them
df_ellipse < levels %>%
lapply(
function(x) make_ellipse(cov_mat, center, level = x)
) %>%
bind_rows() %>%
rename(Intercept = `(Intercept)`, Slope_Days = Days)
df_ellipse
#> # A tibble: 500 × 3
#> Intercept Slope_Days level
#> <dbl> <dbl> <dbl>
#> 1 261. 12.4 0.1
#> 2 260. 12.6 0.1
#> 3 260. 12.7 0.1
#> 4 259. 12.8 0.1
#> 5 258. 12.8 0.1
#> 6 258. 12.9 0.1
#> 7 257. 13.0 0.1
#> 8 257. 13.0 0.1
#> 9 256. 13.1 0.1
#> 10 255. 13.1 0.1
#> # … with 490 more rows
Then we add them onto our previous plot.
ggplot(df_pulled) +
aes(x = Intercept, y = Slope_Days, color = Model, shape = Model) +
# Draw contour lines from the distribution of effects
geom_path(
aes(group = level, color = NULL, shape = NULL),
data = df_ellipse,
linetype = "dashed",
color = "grey40"
) +
geom_point(
aes(shape = Model),
data = df_gravity,
size = 5,
show.legend = FALSE
) +
geom_point(size = 2) +
geom_path(
aes(group = Subject, color = NULL),
arrow = arrow(length = unit(.02, "npc")),
show.legend = FALSE
) +
theme(
legend.position = "bottom",
legend.justification = "right"
) +
ggtitle("Topographic map of regression parameters") +
xlab("Intercept estimate") +
ylab("Slope estimate") +
scale_color_brewer(palette = "Dark2") +
scale_shape_manual(values = c(15:18))
The ellipses provide a little more information about where the estimates are being pulled. Even if some of the points are not being pulled directly towards the center of gravity, nearly all of them are being pulled into a higher confidence region.
There are a few tweaks we might consider for this plot. I don’t think
the ovals need to be contained in the plot. The points are more
important, and the plotting boundaries should be set with respect to the
points. We can redefine the limits by using
coord_cartesian()
. (Your aesthetic preferences
may differ. That’s fine.)
last_plot() +
coord_cartesian(
xlim = range(df_pulled$Intercept),
ylim = range(df_pulled$Slope_Days),
expand = TRUE
)
To go all out 😎, let’s also label the contours with the confidence levels. I see that the lower left area is relatively free of points, so I can place the labels there. I filter down to just the ellipse points in the bottom 25% of x and y values. That will keep points in that lower left quadrant. Then I find the (x, y) point with the farthest distance from the center as the location for my label. I make it sound so easy but it took a lot of trial and error (including an an attempt to use cosines).
# Euclidean distance
contour_dist < function(xs, ys, center_x, center_y) {
x_diff < (center_x  xs) ^ 2
y_diff < (center_y  ys) ^ 2
sqrt(x_diff + y_diff)
}
# Find the point to label in each ellipse.
df_label_locations < df_ellipse %>%
group_by(level) %>%
filter(
Intercept < quantile(Intercept, .25),
Slope_Days < quantile(Slope_Days, .25)
) %>%
# Compute distance from center.
mutate(
dist = contour_dist(Intercept, Slope_Days, fixef(m)[1], fixef(m)[2])
) %>%
# Keep smallest values.
top_n(1, wt = dist) %>%
ungroup()
# Tweak the last plot one more time!
last_plot() +
geom_text(
aes(label = level, color = NULL, shape = NULL),
data = df_label_locations,
nudge_x = .5,
nudge_y = .8,
size = 3.5,
color = "grey40"
)
Are you feeling satisfied? I feel satisfied.
Bonus: Plotting lines from a Bayesian mixed effects model
This last part is more of a code demo than a walkthrough. I call myself a Bayesian. Visualizing uncertainty is one of my things here, so I would be remiss if I didn’t also demo how to do some plots using posterior samples.
Conceptually, the classical model above estimated a single set of partially pooled regression lines. With the Bayesian model, we can sample from a posterior distribution of partially pooled regression lines. Instead of one line for each participant, there’s an entire distribution of them for each participant. This distribution lets us quantify our uncertainty about each part of our model.
First, we fit the model in RStanARM with weakly informative priors.
library(rstanarm)
#> Loading required package: Rcpp
#> This is rstanarm version 2.21.3
#>  See https://mcstan.org/rstanarm/articles/priors for changes to default priors!
#>  Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.
#>  For execution on a local, multicore CPU with excess RAM we recommend calling
#> options(mc.cores = parallel::detectCores())
# Update 202102: Prior to mid2020 priors were autoscaled (so `autoscale =
# TRUE`) was implicity set. But now they are no longer autoscaled. The code has
# been updated to use the autoscaling.
b < stan_glmer(
Reaction ~ Days + (Days  Subject),
family = gaussian(),
data = df_sleep,
prior = normal(0, 2, autoscale = TRUE),
prior_intercept = normal(0, 5, autoscale = TRUE),
prior_covariance = decov(regularization = 2),
prior_aux = cauchy(0, 1, autoscale = TRUE),
# reproducible blogging
seed = 20211116
)
We get a similar overview as arm::display()
when
we print the model.
b
#> stan_glmer
#> family: gaussian [identity]
#> formula: Reaction ~ Days + (Days  Subject)
#> observations: 183
#> 
#> Median MAD_SD
#> (Intercept) 252.3 6.2
#> Days 10.4 1.7
#>
#> Auxiliary parameter(s):
#> Median MAD_SD
#> sigma 25.7 1.6
#>
#> Error terms:
#> Groups Name Std.Dev. Corr
#> Subject (Intercept) 24
#> Days 7 0.06
#> Residual 26
#> Num. levels: Subject 20
#>
#> 
#> * For help interpreting the printed output see ?print.stanreg
#> * For info on the priors used see ?prior_summary.stanreg
We have posterior distribution of values now! That means instead of one
“center of gravity” point, we have 4,000 plausible points for our
central value. The center of our former contour plot has its own contour
plot. That’s Bayes for you. We can plot that easily with
stat_density_2d()
. We set the coordinate limits
to be the same as the last plot, just so that we don’t exaggerate the
uncertainty around the central point by drawing a gigantic contour
surface.
# Get a dataframe: One row per posterior sample
df_posterior < b %>%
as.data.frame() %>%
as_tibble()
ggplot(df_posterior) +
aes(x = `(Intercept)`, y = `Days`) +
# Calculate the density
stat_density_2d(aes(fill = stat(level)), geom = "polygon") +
ggtitle("Where's the average intercept and slope?") +
xlab("Estimate for average intercept") +
ylab("Estimate for average slope") +
# Use the same coordinate limits as last plot
coord_cartesian(
xlim = range(df_pulled$Intercept),
ylim = range(df_pulled$Slope_Days),
expand = TRUE
) +
guides(fill = "none")
For each participant, we have 4,000 partiallypooled regression lines too, so we can visualize our uncertainty for each participant’s individual regression line.
Let’s finish by drawing a sample of those lines for a faceted plot. We have to do a bunch of data wrangling to get a dataframe with one row per subject per posterior sample.
# For each sample, add the average intercept and average slope values to each
# participant's deviation from that average. These yields the intercept and
# slope parameters for each participant.
df_effects < df_posterior %>%
mutate(
# Find all the columns with the pattern "b[(Intercept". Add the column
# `(Intercept)` to each of those columns.
across(
.cols = matches("b\\[\\(Intercept"),
.fns = ~ . + `(Intercept)`
),
# Again for slope
across(
.cols = matches("b\\[Day"),
.fns = ~ . + Days
)
)
# Convert to a long format
df_long_effects < df_effects %>%
select(matches("b\\[")) %>%
rowid_to_column("draw") %>%
tidyr::pivot_longer(
cols = c(draw),
# when we make new columns with pivot_ functions, the
# they get quotes
names_to = "Parameter",
values_to = "Value"
)
# Extract the effect type and subject number from each parameter name
df_long_effects < df_long_effects %>%
mutate(
Effect = Parameter %>%
stringr::str_detect("Intercept") %>%
ifelse(., "Intercept", "Slope_Day"),
Subject = Parameter %>%
stringr::str_extract("\\d\\d\\d")
) %>%
select(draw, Subject, Effect, Value)
# Finally!
df_long_effects
#> # A tibble: 160,000 × 4
#> draw Subject Effect Value
#> <int> <chr> <chr> <dbl>
#> 1 1 308 Intercept 256.
#> 2 1 308 Slope_Day 19.5
#> 3 1 309 Intercept 208.
#> 4 1 309 Slope_Day 2.49
#> 5 1 310 Intercept 197.
#> 6 1 310 Slope_Day 9.00
#> 7 1 330 Intercept 281.
#> 8 1 330 Slope_Day 5.21
#> 9 1 331 Intercept 307.
#> 10 1 331 Slope_Day 1.16
#> # … with 159,990 more rows
Now that we have the data in the right shape, we are going randomly choose 50 posterior samples and plot those lines alongside the observed data.
# For reproducibility
set.seed(20220330)
df_samples < df_long_effects %>%
filter(draw %in% sample(1:4000, size = 50)) %>%
tidyr::pivot_wider(names_from = Effect, values_from = Value)
df_samples
#> # A tibble: 1,000 × 4
#> draw Subject Intercept Slope_Day
#> <int> <chr> <dbl> <dbl>
#> 1 93 308 254. 17.6
#> 2 93 309 217. 2.87
#> 3 93 310 202. 5.61
#> 4 93 330 250. 11.7
#> 5 93 331 262. 8.65
#> 6 93 332 242. 12.0
#> 7 93 333 258. 15.8
#> 8 93 334 232. 12.2
#> 9 93 335 240. 0.869
#> 10 93 337 283. 18.9
#> # … with 990 more rows
ggplot(df_sleep) +
aes(x = Days, y = Reaction) +
geom_abline(
aes(intercept = Intercept, slope = Slope_Day),
data = df_samples,
color = "#3366FF",
alpha = .1
) +
geom_point() +
facet_wrap("Subject") +
scale_x_continuous(breaks = 0:4 * 2) +
labs(x = xlab, y = ylab)
For the participants with complete data, the lines pile up and form a narrow band, indicating a low degree of uncertainty. In the final two panels, however, we only have limited data, and the sample of lines fan out and cover many different plausible trajectories.
The uncertainty is more dramatic if we draw a contour plot for each participant—basically, drawing each participants’ mostly likely locations in the landscape of parameter values.
df_effects_all_wide < df_long_effects %>%
tidyr::pivot_wider(names_from = Effect, values_from = Value)
ggplot(df_effects_all_wide) +
aes(x = Intercept, y = Slope_Day) +
stat_density_2d(
aes(fill = stat(level)),
geom = "polygon",
# normalized density so all colors appear in each plot
contour_var = "ndensity"
) +
facet_wrap("Subject") +
xlab("Intercept estimate") +
ylab("Slope estimate") +
guides(fill = "none")
For 373 and 374, the contour regions/inksplats are very tall: A lot of slope values are plausible. The region for 374 is more off center and slightly narrow than that of 373: That extra data point matters.
Funnily enough, this post started as a quick writeup of a demo I wrote, but it kind of spiraled out of control. I hope this writeup helps students and users understand mixedeffects models at a more intuitive level.
I had formally learned about these models twice in graduate school. In psychology, we were told to use them if we wanted to make inferences about a larger population of subjects or stimulus items. In educational psychology, we were told to use them to capture the sources of variances in a nested dataset: Kids nested in classrooms nested in schools, etc. It wasn’t until I taught myself Bayesian stats that I learned about third reason to use them: They pool information across different units, providing regularized model estimates. I find this rationale most intuitive. The Gelman and Hill book and Statistical Rethinking both discuss the partial pooling description of these models. (Ooooh, as I added the Rethinking link, I just noticed that I created a ggplot2 version of the plot from the cover of that book. 😆)
Last knitted on 20220527. Source code on GitHub.^{1}

.session_info #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 4.2.0 (20220422 ucrt) #> os Windows 10 x64 (build 22000) #> system x86_64, mingw32 #> ui RTerm #> language (EN) #> collate English_United States.utf8 #> ctype English_United States.utf8 #> tz America/Chicago #> date 20220527 #> pandoc NA #> stan (rstan) 2.21.0 #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> ! package * version date (UTC) lib source #> abind 1.45 20160721 [1] CRAN (R 4.2.0) #> arm 1.122 20211015 [1] CRAN (R 4.2.0) #> assertthat 0.2.1 20190321 [1] CRAN (R 4.2.0) #> base64enc 0.13 20150728 [1] CRAN (R 4.2.0) #> bayesplot 1.9.0 20220310 [1] CRAN (R 4.2.0) #> boot 1.328 20210503 [2] CRAN (R 4.2.0) #> cachem 1.0.6 20210819 [1] CRAN (R 4.2.0) #> callr 3.7.0 20210420 [1] CRAN (R 4.2.0) #> cli 3.3.0 20220425 [1] CRAN (R 4.2.0) #> coda 0.194 20200930 [1] CRAN (R 4.2.0) #> codetools 0.218 20201104 [2] CRAN (R 4.2.0) #> colorspace 2.03 20220221 [1] CRAN (R 4.2.0) #> colourpicker 1.1.1 20211004 [1] CRAN (R 4.2.0) #> crayon 1.5.1 20220326 [1] CRAN (R 4.2.0) #> crosstalk 1.2.0 20211104 [1] CRAN (R 4.2.0) #> DBI 1.1.2 20211220 [1] CRAN (R 4.2.0) #> digest 0.6.29 20211201 [1] CRAN (R 4.2.0) #> downlit 0.4.0 20211029 [1] CRAN (R 4.2.0) #> dplyr * 1.0.9 20220428 [1] CRAN (R 4.2.0) #> DT 0.23 20220510 [1] CRAN (R 4.2.0) #> dygraphs 1.1.1.6 20180711 [1] CRAN (R 4.2.0) #> ellipse * 0.4.2 20200527 [1] CRAN (R 4.2.0) #> ellipsis 0.3.2 20210429 [1] CRAN (R 4.2.0) #> evaluate 0.15 20220218 [1] CRAN (R 4.2.0) #> fansi 1.0.3 20220324 [1] CRAN (R 4.2.0) #> farver 2.1.0 20210228 [1] CRAN (R 4.2.0) #> fastmap 1.1.0 20210125 [1] CRAN (R 4.2.0) #> generics 0.1.2 20220131 [1] CRAN (R 4.2.0) #> ggplot2 * 3.3.6 20220503 [1] CRAN (R 4.2.0) #> ggrepel 0.9.1 20210115 [1] CRAN (R 4.2.0) #> ggridges 0.5.3 20210108 [1] CRAN (R 4.2.0) #> git2r 0.30.1 20220316 [1] CRAN (R 4.2.0) #> glue 1.6.2 20220224 [1] CRAN (R 4.2.0) #> gridExtra 2.3 20170909 [1] CRAN (R 4.2.0) #> gtable 0.3.0 20190325 [1] CRAN (R 4.2.0) #> gtools 3.9.2.1 20220523 [1] CRAN (R 4.2.0) #> here 1.0.1 20201213 [1] CRAN (R 4.2.0) #> highr 0.9 20210416 [1] CRAN (R 4.2.0) #> htmltools 0.5.2 20210825 [1] CRAN (R 4.2.0) #> htmlwidgets 1.5.4 20210908 [1] CRAN (R 4.2.0) #> httpuv 1.6.5 20220105 [1] CRAN (R 4.2.0) #> igraph 1.3.1 20220420 [1] CRAN (R 4.2.0) #> inline 0.3.19 20210531 [1] CRAN (R 4.2.0) #> isoband 0.2.5 20210713 [1] CRAN (R 4.2.0) #> knitr * 1.39 20220426 [1] CRAN (R 4.2.0) #> labeling 0.4.2 20201020 [1] CRAN (R 4.2.0) #> later 1.3.0 20210818 [1] CRAN (R 4.2.0) #> lattice 0.2045 20210922 [2] CRAN (R 4.2.0) #> lifecycle 1.0.1 20210924 [1] CRAN (R 4.2.0) #> lme4 * 1.129 20220407 [1] CRAN (R 4.2.0) #> loo 2.5.1 20220324 [1] CRAN (R 4.2.0) #> magrittr 2.0.3 20220330 [1] CRAN (R 4.2.0) #> markdown 1.1 20190807 [1] CRAN (R 4.2.0) #> MASS 7.356 20220323 [2] CRAN (R 4.2.0) #> Matrix * 1.41 20220323 [2] CRAN (R 4.2.0) #> matrixStats 0.62.0 20220419 [1] CRAN (R 4.2.0) #> memoise 2.0.1 20211126 [1] CRAN (R 4.2.0) #> mgcv 1.840 20220329 [2] CRAN (R 4.2.0) #> mime 0.12 20210928 [1] CRAN (R 4.2.0) #> miniUI 0.1.1.1 20180518 [1] CRAN (R 4.2.0) #> minqa 1.2.4 20141009 [1] CRAN (R 4.2.0) #> munsell 0.5.0 20180612 [1] CRAN (R 4.2.0) #> nlme 3.1157 20220325 [2] CRAN (R 4.2.0) #> nloptr 2.0.2 20220519 [1] CRAN (R 4.2.0) #> pillar 1.7.0 20220201 [1] CRAN (R 4.2.0) #> pkgbuild 1.3.1 20211220 [1] CRAN (R 4.2.0) #> pkgconfig 2.0.3 20190922 [1] CRAN (R 4.2.0) #> plyr 1.8.7 20220324 [1] CRAN (R 4.2.0) #> prettyunits 1.1.1 20200124 [1] CRAN (R 4.2.0) #> processx 3.5.3 20220325 [1] CRAN (R 4.2.0) #> promises 1.2.0.1 20210211 [1] CRAN (R 4.2.0) #> ps 1.7.0 20220423 [1] CRAN (R 4.2.0) #> purrr 0.3.4 20200417 [1] CRAN (R 4.2.0) #> R6 2.5.1 20210819 [1] CRAN (R 4.2.0) #> ragg 1.2.2 20220221 [1] CRAN (R 4.2.0) #> RColorBrewer 1.13 20220403 [1] CRAN (R 4.2.0) #> Rcpp * 1.0.8.3 20220317 [1] CRAN (R 4.2.0) #> D RcppParallel 5.1.5 20220105 [1] CRAN (R 4.2.0) #> reshape2 1.4.4 20200409 [1] CRAN (R 4.2.0) #> rlang 1.0.2 20220304 [1] CRAN (R 4.2.0) #> rprojroot 2.0.3 20220402 [1] CRAN (R 4.2.0) #> rstan 2.21.5 20220411 [1] CRAN (R 4.2.0) #> rstanarm * 2.21.3 20220409 [1] CRAN (R 4.2.0) #> rstantools 2.2.0 20220408 [1] CRAN (R 4.2.0) #> rstudioapi 0.13 20201112 [1] CRAN (R 4.2.0) #> scales 1.2.0 20220413 [1] CRAN (R 4.2.0) #> sessioninfo 1.2.2 20211206 [1] CRAN (R 4.2.0) #> shiny 1.7.1 20211002 [1] CRAN (R 4.2.0) #> shinyjs 2.1.0 20211223 [1] CRAN (R 4.2.0) #> shinystan 2.6.0 20220303 [1] CRAN (R 4.2.0) #> shinythemes 1.2.0 20210125 [1] CRAN (R 4.2.0) #> StanHeaders 2.21.07 20201217 [1] CRAN (R 4.2.0) #> stringi 1.7.6 20211129 [1] CRAN (R 4.2.0) #> stringr 1.4.0 20190210 [1] CRAN (R 4.2.0) #> survival 3.31 20220303 [2] CRAN (R 4.2.0) #> systemfonts 1.0.4 20220211 [1] CRAN (R 4.2.0) #> textshaping 0.3.6 20211013 [1] CRAN (R 4.2.0) #> threejs 0.3.3 20200121 [1] CRAN (R 4.2.0) #> tibble * 3.1.7 20220503 [1] CRAN (R 4.2.0) #> tidyr 1.2.0 20220201 [1] CRAN (R 4.2.0) #> tidyselect 1.1.2 20220221 [1] CRAN (R 4.2.0) #> utf8 1.2.2 20210724 [1] CRAN (R 4.2.0) #> vctrs 0.4.1 20220413 [1] CRAN (R 4.2.0) #> withr 2.5.0 20220303 [1] CRAN (R 4.2.0) #> xfun 0.31 20220510 [1] CRAN (R 4.2.0) #> xtable 1.84 20190421 [1] CRAN (R 4.2.0) #> xts 0.12.1 20200909 [1] CRAN (R 4.2.0) #> zoo 1.810 20220415 [1] CRAN (R 4.2.0) #> #> [1] C:/Users/Tristan/AppData/Local/R/winlibrary/4.2 #> [2] C:/Program Files/R/R4.2.0/library #> #> D ── DLL MD5 mismatch, broken installation. #> #> ──────────────────────────────────────────────────────────────────────────────
Leave a comment