Kevin O'Neill
Kevin O'Neill PhD candidate in cognitive neuroscience at Duke University working with Dr. Felipe De Brigard and Dr. John Pearson on causal reasoning.

Intro to Bayesian Regression in R

Intro to Bayesian Regression in R



Welcome! This is an intro-level workshop about Bayesian mixed effects regression in R. We’ll cover the basics of Bayesian linear and logit models. You should have an intermediate-level understanding of R and Frequentist linear regression (using e.g. lm and lmer in R).


Acknowledgments: To make our analyses directly comparable to analyses we’ve already covered, this workshop is directly copied from Allie’s awesome workshop on Frequentist mixed-effect regression. That workshop was adapted from code provided by Gabriela K Hajduk, who in turn referenced a workshop developed by Liam Bailey. Parts of the tutorial are also adapted from a lesson on partial pooling by Tristan Mahr.

For further reading, please check out their tutorials and blogs here:
https://gkhajduk.github.io/2017-03-09-mixed-models/
https://www.tjmahr.com/plotting-partial-pooling-in-mixed-effects-models/




Setup


First, we’ll just get everything set up. We need to tweak some settings, load packages, and read our data.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#change some settings
options(contrasts = c("contr.sum","contr.poly")) 
#this tweaks makes sure that contrasts are interpretable as main effects

#time to load some packages!
library(lme4) #fit the models
library(lmerTest) #gives p-values and more info
library(car) #more settings for regression output
library(tidyr) #for data wrangling
library(dplyr) #for data wrangling
library(tibble) #for data wrangling
library(ggplot2) #plotting raw data
library(data.table) #for pretty HTML tables of model parameters

library(brms)       # bayesian regression!
library(emmeans)    # used to get predicted means per condition
library(modelr)     # used to get predicted means per condition
library(tidybayes)  # for accessing model posteriors 
library(bayestestR) # for testing over posteriors

#load the data
dragons <- read.csv("2020-10-21-dragon-data.csv")




Data


Let’s get familiar with our dataset. This is a fictional dataset about dragons. Each dragon has one row. We have information about each dragon’s body length and cognitive test score. Let’s say our first research question is whether the length of the dragon is related to its intelligence.

We also have some other information about each dragon: We know about the mountain range where it lives, its color, its diet, and whether or not it breathes fire.


Let’s take a look at the data and check the counts of our variables:

1
2
3
4
5
6
7
8
# take a peek at the header
head(dragons)

#check out counts for all our categorical variables
table(dragons$mountainRange)
table(dragons$diet)
table(dragons$color)
table(dragons$breathesFire)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
##    testScore bodyLength mountainRange color      diet breathesFire
## 1  0.0000000   175.5122      Bavarian  Blue Carnivore            1
## 2  0.7429138   190.6410      Bavarian  Blue Carnivore            1
## 3  2.5018247   169.7088      Bavarian  Blue Carnivore            1
## 4  3.3804301   188.8472      Bavarian  Blue Carnivore            1
## 5  4.5820954   174.2217      Bavarian  Blue Carnivore            0
## 6 12.4536350   183.0819      Bavarian  Blue Carnivore            1
## 
## Bavarian  Central Emmental   Julian Ligurian Maritime  Sarntal Southern 
##       60       60       60       60       60       60       60       60 
## 
##  Carnivore   Omnivore Vegetarian 
##        156        167        157 
## 
##   Blue    Red Yellow 
##    160    160    160 
## 
##   0   1 
## 229 251


Now let’s check distributions. Do test scores and body length measurements look approximately normal?

1
2
#check assumptions: do our continuous variables have approximately normal distributions?
hist(dragons$testScore)

1
hist(dragons$bodyLength)




Bayesian Linear Regression


Okay, let’s start fitting some lines! Key Question: Does body length predict test score?

One way to analyse this data would be to try fitting a linear model to all our data, ignoring all the other variables for now. To make sure that we can interpret our coefficients, we should mean-center our continuous measure of body length before using it in a model.

This is a “complete pooling” approach, where we “pool” together all the data and ignore the fact that some observations came from specific mountain ranges.

1
2
model <- lm(testScore ~ scale(bodyLength), data = dragons)
summary(model)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
## 
## Call:
## lm(formula = testScore ~ scale(bodyLength), data = dragons)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -56.962 -16.411  -0.783  15.193  55.200 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        50.3860     0.9676  52.072   <2e-16 ***
## scale(bodyLength)   8.9956     0.9686   9.287   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.2 on 478 degrees of freedom
## Multiple R-squared:  0.1529, Adjusted R-squared:  0.1511 
## F-statistic: 86.25 on 1 and 478 DF,  p-value: < 2.2e-16

Incredible! It’s super significant! We’re gonna publish in Nature Dragonology! How can we run an analogous Bayesian regression to get full posterior distributions of our coefficient? Since we standardized body length, it might be reasonable to set a normal(0, sd(testScore)) prior over the effect of body length, which says that we expect a unit increase in body length to yield on the order of somewhere between a unit decrease and a unit increase in test score. The code is mostly the same except that we use the brm function instead of lm, we specify a file to store our model in (so we don’t have to fit it multiple times), and we specify our normal prior for the model coefficients:

1
2
3
model.bayes <- brm(testScore ~ scale(bodyLength),
                   data=dragons, file='bodyLength',
                   prior=set_prior(paste0('normal(0, ', sd(dragons$testScore), ')')))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
## 
## SAMPLING FOR MODEL '49f78aecf53dc588372d55583a8e39cb' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 2.6e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.26 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.052493 seconds (Warm-up)
## Chain 1:                0.043537 seconds (Sampling)
## Chain 1:                0.09603 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '49f78aecf53dc588372d55583a8e39cb' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 1.3e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.064402 seconds (Warm-up)
## Chain 2:                0.035814 seconds (Sampling)
## Chain 2:                0.100216 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '49f78aecf53dc588372d55583a8e39cb' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 1.2e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.058688 seconds (Warm-up)
## Chain 3:                0.038065 seconds (Sampling)
## Chain 3:                0.096753 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '49f78aecf53dc588372d55583a8e39cb' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 1.3e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.052508 seconds (Warm-up)
## Chain 4:                0.041537 seconds (Sampling)
## Chain 4:                0.094045 seconds (Total)
## Chain 4:
1
summary(model.bayes, prior=TRUE)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: testScore ~ scale(bodyLength) 
##    Data: dragons (Number of observations: 480) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Priors: 
## b ~ normal(0, 23.0088071104418)
## Intercept ~ student_t(3, 51, 26.1)
## sigma ~ student_t(3, 0, 26.1)
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept          50.37      0.98    48.46    52.22 1.00     3608     2184
## scalebodyLength     8.98      0.95     7.14    10.86 1.00     3939     3218
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    21.25      0.67    19.96    22.59 1.00     3826     3030
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

After waiting a minute for the model to compile and fit, we can see that this summary gives us a little more information than lm did. First, it tells us some basics about our model: the noise distribution family is Gaussian with an identity link function, and so on. Next it tells us what our priors are. In this case, brms uses default Student t priors for the intercept and for the standard deviation, and our specified normal prior for the regression slope. Finally, brms tells us about our results. Here Estimate is the posterior mean (which is the most likely for unimodal/symmetric distributions), Est.Error is the standard deviation of the posterior (kind of like the standard error), and l-95% CI and u-95% CI are the credible intervals (values inside this range are the 95% most probable values). We also get some convergence information for each parameter (Rhat, Bulk_ESS, and Tail_ESS), which we’ll talk about later.

You also might be wondering what the heck is the deal with all those chain outputs. Those are just updates on the algorithm brms uses to compute the posterior. It’s called Markov-Chain Monte Carlo (MCMC) sampling and you can learn more about it here if you’re interested.

To perform significance testing on our Bayesian regression model, we can use the nifty describe_posterior function from the bayestestR package, which computes a bunch of different tests for us:

1
2
describe_posterior(model.bayes, ci=.95, rope_ci=.95,
                   test=c('pd', 'p_map', 'rope', 'bf'))
1
2
3
4
5
6
## # Description of Posterior Distributions
## 
## Parameter       | Median |           95% CI | p_MAP |      pd |        95% ROPE | % in ROPE |        BF |  Rhat |      ESS
## --------------------------------------------------------------------------------------------------------------------------
## Intercept       | 50.384 | [48.511, 52.268] |     0 | 100.00% | [-2.301, 2.301] |         0 | 6.641e+52 | 1.000 | 3539.342
## scalebodyLength |  8.971 | [ 7.218, 10.889] |     0 | 100.00% | [-2.301, 2.301] |         0 | 1.316e+06 | 1.000 | 3922.522

Like before, this gives us summaries of our posterior (in this case the median and 95% CI). But this time we also see two measures of effect existence, analogous to p-values on Frequentist approaches. p_MAP is the posterior density at 0 divided by the posterior density at the mode, and is on the same scale as a p-value. pd is the probability of direction, which is the percentage of the posterior that all has the same sign (whichever is greater). We also get two measures of effect significance: the % in ROPE tells us how much of the posterior is inside a null region (in this case a standardized effect size of < .1), and BF is a Bayes Factor, where BF > 1 indicates support for the alternative hypothesis and BF < 1 indicates support for the null hypothesis. There is a whole lot of discussion about which of these metrics are best to use, and we could take a whole meeting talking about this. But for now, we can take solace that all of the measures agree- just like we saw before, it looks like the effect of body length is significant!


We saw that the estimates and inferences from both models look similar, but let’s plot the data with ggplot2 to see how much they actually overlap.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
lm.emm <- emmeans(model, ~bodyLength,
                  at=list(bodyLength=seq(min(dragons$bodyLength),
                                         max(dragons$bodyLength)))) %>%
    as.data.frame
brm.emm <- emmeans(model.bayes, ~bodyLength,
                   at=list(bodyLength=seq(min(dragons$bodyLength),
                                          max(dragons$bodyLength)))) %>%
    as.data.frame

ggplot(dragons, aes(x=bodyLength, y=testScore)) +
    geom_point() +
    geom_line(aes(y=emmean, color='lm'), data=lm.emm) +
    geom_ribbon(aes(y=emmean, ymin=lower.CL, ymax=upper.CL, fill='lm'),
                data=lm.emm, alpha=0.4) +
    geom_line(aes(y=emmean, color='brm'), data=brm.emm) +
    geom_ribbon(aes(y=emmean, ymin=lower.HPD, ymax=upper.HPD, fill='brm'),
                data=brm.emm, alpha=0.4) +
    xlab("Body Length") +
    ylab("Test Score") + theme_minimal()

As you can see, the fitted lines and 95% CIs look almost exactly the same between the two different approaches. Then why bother waiting for a Bayesian model to fit? Because instead of just an estimate and a confidence interval, we get a full posterior distribution over our model coefficients that allows us to directly infer the most probable values:

1
plot(model.bayes)

From these plots, we can see that the average test score is most likely about 50, but could also be somewhere between 48 and 53. Similarly, we can see that a unit increase in body length is most likely to yield an increase in 9 points of test score, and that the standard deviation in test scores is probably around 21 test points. With our Frequentist regression, we can predict that these values make the observed data most likely, but we can’t directly make inferences about how probable the values are.

Then what are the squiggly things on the right side? Since most regression models don’t have easy analytic solutions, we have to sample from our posterior distribution instead of calculating it directly. The lines are called MCMC chains, the x-axis is the number of each sample and the y-axis is the value of each sample. For now, all you need to know is that if the chains look like nice fuzzy caterpillars (like they do here), then our model has converged. Otherwise, brms will give us some warnings that things went awry, and will give us advice for how to solve the problem.

Another benefit of brms is that we can use it to simulate test scores and see how well it covers the actual distribution of test scores. This is called a “posterior predictive check.” Here the dark line is the actual distribution of test scores, and each light line is the distribution of test scores predicted by a single sample of the posterior. Since the light lines mostly cover the solid line, it looks like our model fits the data fairly well:

1
pp_check(model.bayes)




But before we make any grand conclusions, we need to check that we met assumptions! We can use plot for the Frequentist lm model. To check the assumptions of the Bayesian model, we can use add_residual_draws from the tidybayes package to get the residual posterior for each data point.

1
2
3
4
5
6
7
draws <- full_join(add_fitted_draws(dragons, model.bayes),
                   add_residual_draws(dragons, model.bayes)) %>%
    group_by(.row) %>%
    median_hdi(.value, .residual)

## Let's plot the residuals from this model. Ideally, the red line should be flat.
plot(model, which = 1)  # not perfect, but looks alright

1
2
3
4
5
6
ggplot(draws, aes(x=.value, xmin=.value.lower, xmax=.value.upper,
                  y=.residual, ymin=.residual.lower, ymax=.residual.upper)) +
    geom_pointinterval() +
    geom_hline(yintercept=0, linetype='dashed') +
    stat_smooth(se=FALSE) +
    theme_classic()

1
2
## Have a quick look at the  qqplot too - point should ideally fall onto the diagonal dashed line
plot(model, which = 2)  # a bit off at the extremes, but that's often the case; again doesn't look too bad

1
2
3
4
ggplot(draws, aes(sample=.residual)) +
    geom_qq() +
    geom_qq_line() +
    theme_classic()




But linear models also assume that observations are INDEPENDENT. Uh oh.

We collected multiple samples from eight mountain ranges. It’s perfectly plausible that the data from within each mountain range are more similar to each other than the data from different mountain ranges - they are correlated. This could be a problem.

1
2
3
4
5
6
#Lets have a quick look at the data split by mountain range
#We use the facet_wrap to do that
ggplot(data = dragons, aes(x = bodyLength, y = testScore)) +
  geom_point() +
  facet_wrap(.~mountainRange) +
  xlab("length") + ylab("test score")

From the above plots it indeed looks like our mountain ranges vary both in the dragon body length and in their test scores. This confirms that our observations from within each of the ranges aren’t independent. We can’t ignore that.






Bayesian Multilevel Linear Regression


Mountain range clearly introduces a structured source of variance in our data. We need to control for that variation if we want to understand whether body length really predicts test scores.

Multilevel regression is a compromise: Partial pooling! We can let each mountain range have its own regression line, but make an informed guess about that line based on the group-level estimates. This is especially useful when some groups/participants have incomplete data.

Note that here we use the term “multilevel” or “hierarchical” instead of “mixed effect” regression. This is because in a Bayesian framework, all of your effects are modeled as “random” effects! In Frequentist regression, “random” effects are just normal effects that are modeled with an underlying standard normal distribution: in other words, they’re modeled with a standard normal prior. Since all effects in Bayesian regressions have priors, it’s easier and more precise to refer to population-level effects (akin to “fixed” effects) and group-level effects (akin to “random” effects).



Should Mountain Range be a Population-level or Group-level effect?

Since we want to estimate the mean effect of body-length over all mountain ranges, we want a population-level effect of body length. But since we also think that mountain ranges could have different mean test scores and different effects of body length, we also want separate group-level effects over mountain ranges.

Here’s how we did this with lmer:

1
2
3
4
5
#let's fit our first multilevel model!
multilevel_model <- lmer(testScore ~ scale(bodyLength) + (1+scale(bodyLength)|mountainRange), data = dragons)

#what's the verdict?
summary(multilevel_model)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## testScore ~ scale(bodyLength) + (1 + scale(bodyLength) | mountainRange)
##    Data: dragons
## 
## REML criterion at convergence: 3980.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5004 -0.6683  0.0207  0.6592  2.9449 
## 
## Random effects:
##  Groups        Name              Variance Std.Dev. Corr 
##  mountainRange (Intercept)       324.102  18.003        
##                scale(bodyLength)   9.905   3.147   -1.00
##  Residual                        221.578  14.885        
## Number of obs: 480, groups:  mountainRange, 8
## 
## Fixed effects:
##                   Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)       51.75302    6.40349  6.91037   8.082 9.16e-05 ***
## scale(bodyLength) -0.03326    1.68317 10.22337  -0.020    0.985    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(bdyLng) -0.674
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular

Again, the Bayesian version with brms is essentially the same to run:

1
2
3
4
5
6
#let's fit our first mixed model!
multilevel_model.bayes <- brm(testScore ~ scale(bodyLength) + (1+scale(bodyLength)|mountainRange),
                              data=dragons, file='bodyLength_multilevel',
                              prior=set_prior(paste0('normal(0, ', sd(dragons$testScore), ')')))

summary(multilevel_model.bayes, prior=TRUE)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: testScore ~ scale(bodyLength) + (1 + scale(bodyLength) | mountainRange) 
##    Data: dragons (Number of observations: 480) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Priors: 
## b ~ normal(0, 23.0088071104418)
## Intercept ~ student_t(3, 51, 26.1)
## L ~ lkj_corr_cholesky(1)
## sd ~ student_t(3, 0, 26.1)
## sigma ~ student_t(3, 0, 26.1)
## 
## Group-Level Effects: 
## ~mountainRange (Number of levels: 8) 
##                                Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                     20.63      5.92    12.19    35.13 1.00
## sd(scalebodyLength)                3.76      2.17     0.35     8.89 1.00
## cor(Intercept,scalebodyLength)    -0.55      0.39    -0.98     0.45 1.00
##                                Bulk_ESS Tail_ESS
## sd(Intercept)                      1362     1915
## sd(scalebodyLength)                1820     1699
## cor(Intercept,scalebodyLength)     3176     2620
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept          51.25      7.11    36.62    64.88 1.00     1040     1590
## scalebodyLength     0.12      1.97    -3.97     3.91 1.00     2340     1936
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    14.93      0.49    14.03    15.93 1.00     4904     2967
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
1
2
describe_posterior(multilevel_model.bayes, ci=.95, rope_ci=.95,
                   test=c('pd', 'p_map', 'rope', 'bf'))
1
2
3
4
5
6
## # Description of Posterior Distributions
## 
## Parameter       | Median |           95% CI | p_MAP |      pd |        95% ROPE | % in ROPE |        BF |  Rhat |      ESS
## --------------------------------------------------------------------------------------------------------------------------
## Intercept       | 51.279 | [36.597, 64.809] | 0.000 | 100.00% | [-2.301, 2.301] |     0.000 | 91530.079 | 1.000 | 1018.165
## scalebodyLength |  0.200 | [-3.993,  3.871] | 0.979 |  54.57% | [-2.301, 2.301] |    82.163 |     0.078 | 1.000 | 2210.472
1
pp_check(multilevel_model.bayes)

The summary output for brms is pretty much the same as before, except now we also have a section for group-level effects in addition to our earlier section for Population-level effects. Since we estimated separate intercepts and effects of body length on test score for each mountain range, our model tells us the standard deviation of each of those two effects, and also the correlation between them. The standard deviations look fairly similar to the lmer values, but the correlation looks much more reasonable than the lmer value (i.e., it’s no longer a perfect negative correlation).


Overall, it looks like when we account for the effect of mountain range, there is no relationship between body length and test scores. This is true for both the Frequentist and the Bayesian regressions. Well, so much for our Nature Dragonology paper!

Unless… What about our other variables? Let’s test whether diet is related to test scores instead.

1
2
3
4
5
model.diet <- brm(testScore ~ diet + (1+diet|mountainRange),
                        data=dragons, file='diet_multilevel',
                        prior=set_prior(paste0('normal(0, ', sd(dragons$testScore), ')')))

summary(model.diet, prior=TRUE)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: testScore ~ diet + (1 + diet | mountainRange) 
##    Data: dragons (Number of observations: 480) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Priors: 
## b ~ normal(0, 23.0088071104418)
## Intercept ~ student_t(3, 51, 26.1)
## L ~ lkj_corr_cholesky(1)
## sd ~ student_t(3, 0, 26.1)
## sigma ~ student_t(3, 0, 26.1)
## 
## Group-Level Effects: 
## ~mountainRange (Number of levels: 8) 
##                      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)           11.47      3.82     6.57    20.22 1.00     1311
## sd(diet1)                2.70      1.98     0.14     7.47 1.00      935
## sd(diet2)                4.92      2.01     1.85     9.86 1.00     1475
## cor(Intercept,diet1)     0.08      0.44    -0.79     0.83 1.00     3302
## cor(Intercept,diet2)     0.72      0.24     0.10     0.98 1.00     1950
## cor(diet1,diet2)        -0.29      0.43    -0.90     0.69 1.00     2261
##                      Tail_ESS
## sd(Intercept)            1923
## sd(diet1)                1867
## sd(diet2)                1961
## cor(Intercept,diet1)     2301
## cor(Intercept,diet2)     2866
## cor(diet1,diet2)         2887
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    49.45      4.16    41.08    57.96 1.00     1168     1650
## diet1       -15.02      1.61   -18.09   -11.41 1.00     2308     1715
## diet2        15.00      2.13    10.57    18.89 1.00     1653     1758
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    10.74      0.35    10.05    11.46 1.00     4776     2879
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
1
2
describe_posterior(model.diet, ci=.95, rope_ci=.95,
                   test=c('pd', 'p_map', 'rope', 'bf'))
1
2
3
4
5
6
7
## # Description of Posterior Distributions
## 
## Parameter |  Median |             95% CI | p_MAP |      pd |        95% ROPE | % in ROPE |        BF |  Rhat |      ESS
## -----------------------------------------------------------------------------------------------------------------------
## Intercept |  49.455 | [ 40.919,  57.793] |     0 | 100.00% | [-2.301, 2.301] |         0 | 9.108e+08 | 1.001 | 1161.725
## diet1     | -15.090 | [-18.170, -11.587] |     0 | 100.00% | [-2.301, 2.301] |         0 |  2381.239 | 1.001 | 2050.769
## diet2     |  15.072 | [ 10.533,  18.810] |     0 | 100.00% | [-2.301, 2.301] |         0 |  2135.790 | 1.001 | 1639.868
1
pp_check(model.diet)

Visualizing the effect of diet on test scores predicted by our model will help us better understand what we just found:

1
2
3
4
5
6
7
8
#Plot average test score by diet type
dragons %>%
    data_grid(diet) %>%
    add_fitted_draws(model.diet, re_formula=NA) %>%
    ggplot(aes(x=diet, y=.value)) +
    stat_halfeye(point_interval=median_hdi) + ylim(0, NA) +
    xlab("Diet") + ylab("Test Score") +
    theme_minimal()

1
2
3
4
5
6
7
8
9
#Let's also look at the effect across mountain ranges.
dragons %>%
    data_grid(diet, mountainRange) %>%
    add_fitted_draws(model.diet) %>%
    ggplot(aes(x=diet, y=.value)) +
    stat_halfeye(point_interval=median_hdi) + ylim(0, NA) +
    facet_wrap( ~ mountainRange) +
    xlab("Diet") + ylab("Test Score") +
    theme_minimal()

What are we looking at here? The black points represent the mode of the posterior, the thick black lines represent the 66% credible intervals, and the thinner black lines represnt the 95% credible intervals (also called highest density intervals). The curves in gray represent the full posterior distribution. A nice thing about using a Bayesian model is that instead of being stuck with ugly & often misleading bar charts, we can directly plot how likely each mean is along with the full posterior distribution.

In sum, these results look pretty consistent, but there’s clearly still variability among different mountains.




So far, so good. But what if we’re interested in testing multiple variables at the same time? We build a model with an interaction term!

1
2
3
4
5
model.diet.length <- brm(testScore ~ diet*scale(bodyLength) + (1 + diet*scale(bodyLength)|mountainRange),
                         data=dragons, file='diet_length',
                         prior=set_prior(paste0('normal(0, ', sd(dragons$testScore), ')')))

summary(model.diet.length, prior=TRUE)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: testScore ~ diet * scale(bodyLength) + (1 + diet * scale(bodyLength) | mountainRange) 
##    Data: dragons (Number of observations: 480) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Priors: 
## b ~ normal(0, 23.0088071104418)
## Intercept ~ student_t(3, 51, 26.1)
## L ~ lkj_corr_cholesky(1)
## sd ~ student_t(3, 0, 26.1)
## sigma ~ student_t(3, 0, 26.1)
## 
## Group-Level Effects: 
## ~mountainRange (Number of levels: 8) 
##                                                  Estimate Est.Error l-95% CI
## sd(Intercept)                                       11.52      4.04     6.30
## sd(diet1)                                            2.79      2.01     0.19
## sd(diet2)                                            3.33      2.07     0.31
## sd(scalebodyLength)                                  1.89      1.42     0.07
## sd(diet1:scalebodyLength)                            1.68      1.55     0.05
## sd(diet2:scalebodyLength)                            2.00      1.59     0.07
## cor(Intercept,diet1)                                 0.09      0.35    -0.60
## cor(Intercept,diet2)                                 0.37      0.34    -0.39
## cor(diet1,diet2)                                    -0.14      0.37    -0.77
## cor(Intercept,scalebodyLength)                      -0.14      0.36    -0.76
## cor(diet1,scalebodyLength)                           0.02      0.38    -0.71
## cor(diet2,scalebodyLength)                          -0.11      0.38    -0.75
## cor(Intercept,diet1:scalebodyLength)                -0.01      0.38    -0.74
## cor(diet1,diet1:scalebodyLength)                     0.05      0.38    -0.69
## cor(diet2,diet1:scalebodyLength)                    -0.06      0.37    -0.74
## cor(scalebodyLength,diet1:scalebodyLength)           0.03      0.37    -0.67
## cor(Intercept,diet2:scalebodyLength)                -0.15      0.38    -0.79
## cor(diet1,diet2:scalebodyLength)                     0.03      0.38    -0.70
## cor(diet2,diet2:scalebodyLength)                    -0.15      0.37    -0.79
## cor(scalebodyLength,diet2:scalebodyLength)           0.04      0.38    -0.68
## cor(diet1:scalebodyLength,diet2:scalebodyLength)    -0.02      0.38    -0.72
##                                                  u-95% CI Rhat Bulk_ESS
## sd(Intercept)                                       21.36 1.00     1599
## sd(diet1)                                            7.89 1.00     1650
## sd(diet2)                                            8.22 1.00     1506
## sd(scalebodyLength)                                  5.32 1.00     2439
## sd(diet1:scalebodyLength)                            5.78 1.00     2221
## sd(diet2:scalebodyLength)                            5.98 1.00     2163
## cor(Intercept,diet1)                                 0.75 1.00     5263
## cor(Intercept,diet2)                                 0.88 1.00     4551
## cor(diet1,diet2)                                     0.62 1.00     3261
## cor(Intercept,scalebodyLength)                       0.59 1.00     5336
## cor(diet1,scalebodyLength)                           0.72 1.00     4651
## cor(diet2,scalebodyLength)                           0.66 1.00     3869
## cor(Intercept,diet1:scalebodyLength)                 0.70 1.00     6345
## cor(diet1,diet1:scalebodyLength)                     0.74 1.00     4627
## cor(diet2,diet1:scalebodyLength)                     0.65 1.00     4318
## cor(scalebodyLength,diet1:scalebodyLength)           0.71 1.00     3501
## cor(Intercept,diet2:scalebodyLength)                 0.60 1.00     5048
## cor(diet1,diet2:scalebodyLength)                     0.73 1.00     4024
## cor(diet2,diet2:scalebodyLength)                     0.61 1.00     4521
## cor(scalebodyLength,diet2:scalebodyLength)           0.74 1.00     3584
## cor(diet1:scalebodyLength,diet2:scalebodyLength)     0.70 1.00     2956
##                                                  Tail_ESS
## sd(Intercept)                                        1698
## sd(diet1)                                            2254
## sd(diet2)                                            1850
## sd(scalebodyLength)                                  2350
## sd(diet1:scalebodyLength)                            2538
## sd(diet2:scalebodyLength)                            1896
## cor(Intercept,diet1)                                 2858
## cor(Intercept,diet2)                                 3480
## cor(diet1,diet2)                                     3377
## cor(Intercept,scalebodyLength)                       2747
## cor(diet1,scalebodyLength)                           3037
## cor(diet2,scalebodyLength)                           3611
## cor(Intercept,diet1:scalebodyLength)                 2775
## cor(diet1,diet1:scalebodyLength)                     3029
## cor(diet2,diet1:scalebodyLength)                     2766
## cor(scalebodyLength,diet1:scalebodyLength)           3280
## cor(Intercept,diet2:scalebodyLength)                 3217
## cor(diet1,diet2:scalebodyLength)                     2959
## cor(diet2,diet2:scalebodyLength)                     3098
## cor(scalebodyLength,diet2:scalebodyLength)           3590
## cor(diet1:scalebodyLength,diet2:scalebodyLength)     3481
## 
## Population-Level Effects: 
##                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                50.01      4.35    41.32    58.68 1.00     1159
## diet1                   -15.06      1.72   -18.34   -11.23 1.00     2522
## diet2                    15.50      1.82    11.66    18.86 1.00     2537
## scalebodyLength          -0.23      1.30    -2.80     2.42 1.00     4099
## diet1:scalebodyLength    -0.56      1.55    -3.54     2.50 1.00     2753
## diet2:scalebodyLength     1.73      1.53    -1.48     4.68 1.00     2898
##                       Tail_ESS
## Intercept                 1220
## diet1                     2243
## diet2                     2100
## scalebodyLength           2824
## diet1:scalebodyLength     2250
## diet2:scalebodyLength     2588
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    10.78      0.35    10.11    11.48 1.00     6093     2712
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
1
2
describe_posterior(model.diet.length, ci=.95, rope_ci=.95,
                   test=c('pd', 'p_map', 'rope', 'bf'))
1
2
3
4
5
6
7
8
9
10
## # Description of Posterior Distributions
## 
## Parameter             |  Median |             95% CI | p_MAP |      pd |        95% ROPE | % in ROPE |        BF |  Rhat |      ESS
## -----------------------------------------------------------------------------------------------------------------------------------
## Intercept             |  49.960 | [ 40.958,  58.301] | 0.000 | 100.00% | [-2.301, 2.301] |     0.000 | 3.204e+06 | 1.004 | 1064.403
## diet1                 | -15.133 | [-18.368, -11.330] | 0.000 | 100.00% | [-2.301, 2.301] |     0.000 | 15268.750 | 1.000 | 2357.941
## diet2                 |  15.593 | [ 11.707,  18.894] | 0.000 | 100.00% | [-2.301, 2.301] |     0.000 | 28342.766 | 1.002 | 2425.588
## scalebodyLength       |  -0.240 | [ -3.000,   2.199] | 0.973 |  57.83% | [-2.301, 2.301] |    96.212 |     0.054 | 1.000 | 4010.814
## diet1.scalebodyLength |  -0.607 | [ -3.375,   2.650] | 0.845 |  66.72% | [-2.301, 2.301] |    89.766 |     0.069 | 1.000 | 2464.699
## diet2.scalebodyLength |   1.739 | [ -1.383,   4.732] | 0.418 |  87.30% | [-2.301, 2.301] |    65.588 |     0.140 | 1.000 | 2818.100
1
pp_check(model.diet.length)

Since the BFs for diet are > 10, it looks like the effect of diet is still significant after we control for body length. We can also see that there seems to be no main effect of body length or interactions between diet and body length, since their corresponding BFs are < .14.

Let’s plot the output again:

1
2
3
4
5
6
7
8
9
10
dragons %>%
    data_grid(diet, bodyLength=seq_range(bodyLength, 50)) %>%
    add_fitted_draws(model.diet.length, re_formula=NA) %>%
    median_hdi() %>%
    ggplot(aes(x=bodyLength, y=.value, color=diet, fill=diet)) +
    geom_line(size=2) +
    geom_ribbon(aes(ymin=.lower, ymax=.upper), alpha=0.3) +
    ylim(0, NA) +
    xlab("Body Length") + ylab("Test Score") +
    theme_minimal()

Here we can see that indeed, the effect of body length seems to be close to zero for all three diet types.




Hmm…. did adding body length to our model make it better? Since we’re Bayesian statisticians, we can compare distributions of adjusted R2 values for each model to answer this question:

1
2
3
4
5
6
7
8
9
R2 <- data.frame(model.diet=loo_R2(model.diet, summary=FALSE)[,1],
                 model.diet.length=loo_R2(model.diet.length, summary=FALSE)[,1]) %>%
    mutate(diff=model.diet - model.diet.length) %>%
    pivot_longer(model.diet:diff) %>%
    mutate(name=factor(name, levels=c('model.diet', 'model.diet.length', 'diff')))

ggplot(R2, aes(x=name, y=value)) +
    stat_halfeye(point_interval=median_hdi) +
    theme_minimal()

As we can see, both models have an R2 of about 0.75, so adding body length doesn’t seem to help. We can also test their difference using LOO-IC:

1
loo(model.diet, model.diet.length)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
## Output of model 'model.diet':
## 
## Computed from 4000 by 480 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo  -1832.9 21.6
## p_loo        24.4  4.0
## looic      3665.9 43.2
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     474   98.8%   1097      
##  (0.5, 0.7]   (ok)         4    0.8%   230       
##    (0.7, 1]   (bad)        2    0.4%   87        
##    (1, Inf)   (very bad)   0    0.0%   <NA>      
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model.diet.length':
## 
## Computed from 4000 by 480 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo  -1838.1 21.4
## p_loo        30.4  4.1
## looic      3676.2 42.7
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     474   98.8%   1005      
##  (0.5, 0.7]   (ok)         5    1.0%   365       
##    (0.7, 1]   (bad)        1    0.2%   140       
##    (1, Inf)   (very bad)   0    0.0%   <NA>      
## See help('pareto-k-diagnostic') for details.
## 
## Model comparisons:
##                   elpd_diff se_diff
## model.diet         0.0       0.0   
## model.diet.length -5.1       1.8

Again we see that the model without diet does just as well, if not slightly better than, the model with both diet and body length.




Your turn! Try modifying the model above to test whether color is related to testScore, and whether color interacts with diet or bodyLength.

1
2
3
4
5
6
7
8
#Build your model here


#View the output of the model here



#Plot your results below:




Bayesian Multilevel Logistic Regression

Okay, let’s test a new question. Test scores are boring. I actually want to know about which dragons breathe fire. This has way more important practical implications, and is more likely to get me grant funding.

Good news: we have data on fire breathing!

Bad news: It’s a binary variable, so we need to change our model.

With lmer, you need to switch over to the glmer function to gain access to bernoulli models. But in brms, you just need to specify the proper noise distribution family:

1
2
3
4
5
logit_model <- brm(breathesFire ~ color + (1+color|mountainRange),
                   data=dragons, file='fire', family=bernoulli,
                   prior=prior(normal(0, 2)))

summary(logit_model, prior=TRUE)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: breathesFire ~ color + (1 + color | mountainRange) 
##    Data: dragons (Number of observations: 480) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Priors: 
## b ~ normal(0, 2)
## Intercept ~ student_t(3, 0, 2.5)
## L ~ lkj_corr_cholesky(1)
## sd ~ student_t(3, 0, 2.5)
## 
## Group-Level Effects: 
## ~mountainRange (Number of levels: 8) 
##                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)             0.40      0.27     0.03     1.04 1.00     1346
## sd(color1)                0.35      0.28     0.01     1.04 1.00     1597
## sd(color2)                1.03      0.45     0.39     2.13 1.00     1367
## cor(Intercept,color1)    -0.23      0.48    -0.94     0.77 1.00     2969
## cor(Intercept,color2)     0.15      0.44    -0.72     0.88 1.00     1228
## cor(color1,color2)       -0.07      0.49    -0.88     0.84 1.00     1132
##                       Tail_ESS
## sd(Intercept)             1543
## sd(color1)                1163
## sd(color2)                1103
## cor(Intercept,color1)     2705
## cor(Intercept,color2)     1289
## cor(color1,color2)        1570
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.03      0.22    -0.40     0.48 1.00     2196     2212
## color1        2.22      0.26     1.71     2.74 1.00     3108     2349
## color2        0.35      0.42    -0.49     1.17 1.00     1817     1923
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).


Now let’s plot the proportion of dragons that breathe fire by color:

1
2
3
4
5
6
7
8
dragons %>% data_grid(color) %>%
    add_fitted_draws(logit_model, re_formula=NA, scale='response') %>%
    ggplot(aes(x=color, y=.value, fill=color)) +
    stat_halfeye(point_interval=median_hdi, show.legend=FALSE) +
    scale_fill_manual(values=c('blue', 'red', 'yellow')) +
    xlab("Color") +
    ylab("Proportion that Breathes Fire") +
    theme_minimal()

Looks like most blue dragons breathe fire, red dragons are only slightly more likely to breathe fire than not, and few yellow dragons breathe fire.


Your turn! Test whether other variables predict breathesFire.

1
2
3
4
5
#Build your model here:



#Plot your results:

Bayesian Multinomial Regression

One serious advantage of going Bayesian is that you get immediate access to lots of other model types that require significant effort to get from existing Frequentist packages. One example is multinomial regression- this is like logistic regression, but with more than two unordered categories.

Instead of using color to predict body length, let’s see if body length predicts color:

1
2
3
4
5
multi_model <- brm(color ~ bodyLength,
                   data=dragons, file='color', family=categorical,
                   prior=prior(normal(0, 2)))

summary(multi_model, prior=TRUE)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
##  Family: categorical 
##   Links: muRed = logit; muYellow = logit 
## Formula: color ~ bodyLength 
##    Data: dragons (Number of observations: 480) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Priors: 
## b_muRed ~ normal(0, 2)
## b_muYellow ~ normal(0, 2)
## 
## Population-Level Effects: 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## muRed_Intercept         5.62      1.46     2.76     8.52 1.00     3441     3214
## muYellow_Intercept     -7.12      1.67   -10.29    -3.86 1.00     3229     2857
## muRed_bodyLength       -0.03      0.01    -0.04    -0.01 1.00     3399     3110
## muYellow_bodyLength     0.03      0.01     0.02     0.05 1.00     3289     2852
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).


Let’s plot the proportion of dragons of each color by body length:

1
2
3
4
5
6
7
8
9
10
dragons %>% data_grid(bodyLength=seq_range(bodyLength, 100)) %>%
    add_fitted_draws(multi_model, scale='response') %>%
    ggplot(aes(x=bodyLength, y=.value)) +
    stat_lineribbon(aes(color=.category, fill=.category),
                    alpha=0.4, show.legend=FALSE) +
    scale_color_manual(values=c('blue', 'red', 'yellow')) +
    scale_fill_manual(values=c('blue', 'red', 'yellow')) +
    xlab("Body Length") +
    ylab("Probability of Color") +
    theme_minimal()

We can see that while short dragons are likely to be red, long dragons are likely to be yellow, and mid-length dragons are likely to be blue.

Bayesian Multilevel Unequal Variances Regression

Another extremely useful example of a model that integrates seamlessly in brms is unequal variance regression. In our model using diet to predict test scores, we never checked that the variances in test scores between dragons with different diets was equal. Since this is an assumption of that model, it could be bad if that assumption doesn’t hold. An easy way to get around this in Bayesian regression is to simply estimate both a mean and a variance for each diet:

1
2
3
4
5
6
model.diet.uv <- brm(bf(testScore ~ diet + (1+diet|mountainRange),
                        sigma ~ diet + (1+diet|mountainRange)),
                     data=dragons, file='diet_uv',
                     prior=set_prior(paste0('normal(0, ', sd(dragons$testScore), ')')))

summary(model.diet.uv, prior=TRUE)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: testScore ~ diet + (1 + diet | mountainRange) 
##          sigma ~ diet + (1 + diet | mountainRange)
##    Data: dragons (Number of observations: 480) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Priors: 
## b ~ normal(0, 23.0088071104418)
## Intercept ~ student_t(3, 51, 26.1)
## Intercept_sigma ~ student_t(3, 0, 2.5)
## L ~ lkj_corr_cholesky(1)
## sd ~ student_t(3, 0, 26.1)
## sd_sigma ~ student_t(3, 0, 26.1)
## 
## Group-Level Effects: 
## ~mountainRange (Number of levels: 8) 
##                                  Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                       10.78      3.40     6.15    19.37 1.00
## sd(diet1)                            2.33      1.75     0.11     6.64 1.00
## sd(diet2)                            4.47      1.95     1.22     9.15 1.00
## sd(sigma_Intercept)                  0.13      0.09     0.01     0.35 1.00
## sd(sigma_diet1)                      0.22      0.13     0.02     0.54 1.00
## sd(sigma_diet2)                      0.15      0.11     0.01     0.41 1.00
## cor(Intercept,diet1)                 0.03      0.48    -0.85     0.87 1.00
## cor(Intercept,diet2)                 0.74      0.25     0.05     0.99 1.00
## cor(diet1,diet2)                    -0.28      0.45    -0.93     0.74 1.00
## cor(sigma_Intercept,sigma_diet1)     0.08      0.46    -0.80     0.87 1.00
## cor(sigma_Intercept,sigma_diet2)    -0.25      0.46    -0.93     0.73 1.00
## cor(sigma_diet1,sigma_diet2)        -0.42      0.45    -0.96     0.68 1.00
##                                  Bulk_ESS Tail_ESS
## sd(Intercept)                        1699     2357
## sd(diet1)                            1709     2025
## sd(diet2)                            1530     1015
## sd(sigma_Intercept)                  1036     1492
## sd(sigma_diet1)                      1042     1528
## sd(sigma_diet2)                      1094     1931
## cor(Intercept,diet1)                 3565     2868
## cor(Intercept,diet2)                 2202     2767
## cor(diet1,diet2)                     3147     2916
## cor(sigma_Intercept,sigma_diet1)     2046     2306
## cor(sigma_Intercept,sigma_diet2)     2546     2460
## cor(sigma_diet1,sigma_diet2)         2248     2533
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept          49.20      4.04    41.03    57.44 1.00     1289     1798
## sigma_Intercept     2.39      0.08     2.25     2.54 1.00     2288     2156
## diet1             -15.53      1.62   -18.63   -12.27 1.00     1900     2589
## diet2              15.41      2.09    11.08    19.37 1.00     1885     2209
## sigma_diet1         0.16      0.11    -0.05     0.41 1.00     2260     2164
## sigma_diet2         0.02      0.09    -0.14     0.22 1.00     2306     1814
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Let’s plot the predicted mean and standard deviation of test scores by diet:

1
2
3
4
5
6
7
8
9
draws <- dragons %>%
    data_grid(diet) %>%
    add_fitted_draws(model.diet.uv, dpar='sigma', re_formula=NA)

draws %>%
    ggplot(aes(x=diet, y=.value)) +
    stat_halfeye(point_interval=median_hdi) + ylim(0, NA) +
    xlab("Diet") + ylab("Test Score") +
    theme_minimal()

1
2
3
4
5
draws %>%
    ggplot(aes(x=diet, y=sigma)) +
    stat_halfeye(point_interval=median_hdi) + ylim(0, NA) +
    xlab("Diet") + ylab("Std Deviation Test Score") +
    theme_minimal()

It looks like the variation in test scores might be smaller for vegetarian dragons than for carnivorous dragons, though the difference is small. Does this unequal variance model do better than our old model?

1
loo(model.diet, model.diet.uv)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
## Output of model 'model.diet':
## 
## Computed from 4000 by 480 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo  -1832.9 21.6
## p_loo        24.4  4.0
## looic      3665.9 43.2
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     474   98.8%   1097      
##  (0.5, 0.7]   (ok)         4    0.8%   230       
##    (0.7, 1]   (bad)        2    0.4%   87        
##    (1, Inf)   (very bad)   0    0.0%   <NA>      
## See help('pareto-k-diagnostic') for details.
## 
## Output of model 'model.diet.uv':
## 
## Computed from 4000 by 480 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo  -1832.7 21.8
## p_loo        42.3  6.7
## looic      3665.3 43.7
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     469   97.7%   627       
##  (0.5, 0.7]   (ok)         4    0.8%   176       
##    (0.7, 1]   (bad)        7    1.5%   18        
##    (1, Inf)   (very bad)   0    0.0%   <NA>      
## See help('pareto-k-diagnostic') for details.
## 
## Model comparisons:
##               elpd_diff se_diff
## model.diet.uv  0.0       0.0   
## model.diet    -0.3       6.6

In this case, it looks like the model without equal variances does about as well as the model with unequal variances. That means we can probably assume equal variances without any problems. But this isn’t always the case, and many research programs are dedicated to explaining why different groups have different variances!




Conclusions

In this tutorial we demonstrated the basics of using brms to run Bayesian regressions that directly parallel what you’re likely used to running with lm and lmer. We also demonstrated how to run Bayesian significant tests with bayestestR and how to plot results from these models using tidybayes. Though there are many more details to specifying, running, and interpreting Bayesian models, hopefully this tutorial convinced you that making the first step is easier than you imagined!




Convergence

Some parting tips and tricks:

Like with Frequentist multilevel models, one of the biggest concerns with convergence is whether you have enough data for your model structure. Specifically, you need enough observations for every combination of population-level and group-level effect in order for your model to be well-defined. Unlike lmer, brms will try to fit a model even if you don’t have enough data. But in this case, you will either get unstable results, or your model will only be informed by your priors.

There are a few different convergence problems that brms will tell you about, however. If Rhat for any of your parameters is greater than 1.1, then your model has not converged and the samples from your model don’t match the true posterior. Here, brms will likely instruct you to increase the number of samples and/or turn on thinning, which reduces the autocorrelation between samples. If you have divergent transitions, then brms will tell you to increase the adapt_delta parameter, which lowers the step size of the estimator. Finally, you might get warnings about max_treedepth or ESS, both of which brms will give you recommendations for.

In my experience, however, one of the main reasons that these warnings might persist is that your priors are underspecified. If you have flat or extremely wide priors, then brms has to search a massive parameter space to find the most likely values. But if you constrain your prior to a reasonable degree, not only will your model fit faster, but it will also have less convergence issues. No matter what, you want to make sure that your priors truly reflect your expectations and that your results aren’t too heavily influenced by your specific choice of prior.