A crash course in Bayesian mixed models with brms (Lesson 5)

Introduction and review

What is this class?

Part 5 of a practical introduction to fitting Bayesian multilevel models in R and Stan
Uses the brms package (Bayesian Regression Models using Stan)

How to follow the course

Slides and text version of lessons are online
Fill in code in the worksheet (replace ... with code)
You can always copy and paste code from text version of lesson if you fall behind

What we’ve learned in Lessons 1-4

Relationship between prior, likelihood, and posterior
How to sample from the posterior with Hamiltonian Monte Carlo
How to use posterior distributions from our models to make predictions and test hypotheses
GLMMs with non-normal response distributions and random intercepts and slopes
Models with residuals correlated in space and time
Spatial generalized additive mixed models
Beta regression for proportion data
Zero-inflated and hurdle models for data with zeros

What will we learn in this lesson?

Cumulative logistic mixed models for ordered categorical data
Nonlinear mixed models

Big-picture concepts

Explore in more detail how link functions work in mixed models
Discover how every parameter in a mixed model can have fixed and random effects

Part 1. Cumulative logistic mixed models for ordinal data

Load packages and set options

library(brms)
library(agridat)
library(tidyverse)
library(emmeans)
library(tidybayes)
library(ggdist)
library(bayestestR)
library(gt)

theme_set(theme_bw())
colorblind_friendly <- scale_color_manual(values = unname(palette.colors(5)[-1]), aesthetics = c('colour', 'fill'))

# IF YOU ARE RUNNING THIS CODE ON YOUR OWN MACHINE
options(mc.cores = 4, brms.file_refit = 'on_change')

# IF YOU ARE RUNNING THIS CODE ON THE POSIT CLOUD SERVER (or you have configured CmdStan on your own machine)
options(brms.backend = 'cmdstanr', mc.cores = 4, brms.file_refit = 'on_change')

Categorical data

Image (c) Allison Horst

We’ve seen binary response data in lesson 2
Nominal data (for example: blood group) has statistical models but we won’t cover them today

Ordered categorical data

Lots of data in ag science are on an ordinal rating scale
Common in phenotyping: capture a complex phenotype, or a quantitative one that is too time-consuming to measure directly, with a simple score
- Rate disease severity on a scale from 1 to 5
- Instead of counting hairs on a stem, rate low/medium/high

Shortcut at your own risk

It is faster to score phenotypes on a rating scale than to physically measure the underlying variable directly
But this speed has a price: more data points are required for similarly precise estimates, than if you had a continuous variable
Some people try to have it both ways and treat the numerical ratings as continuous variables with normally distributed error
This “quick and dirty” shortcut often works OK
But why do that when you have statistical models for ordinal response data?!?

Review: modeling binary data

For binary data we use a logistic model: binomial response distribution with logit link function
We’re trying to predict a probability, which has to be between 0 and 1
But our linear predictor is a combination of fixed and random terms, multiplied by coefficients and added together, which can be any value, positive or negative
Logit, or log-odds, link function takes a probability between 0 and 1 and maps it to a scale that can have any value, positive or negative
- $\text{logit } p = \log \frac{p}{1-p}$

Binomial distribution for binary data

Logit model is modeling a binary outcome, we only need to predict single probability
Probability of success + probability of failure = 1
If we know one, we know the other
Model this process as following a binomial distribution which only has one parameter
Example: model coin flips as following a binomial distribution with one parameter: $\theta$, the probability of flipping heads

Multinomial distribution, from coins to dice

If we have $n$ ordered categories, we have to estimate $n - 1$ probabilities
They all have to be between 0 and 1
Instead of a binomial distribution, use a multinomial distribution.
$n - 1$ parameters, the probability of each outcome except for the last
You can model a dice roll as following a multinomial distribution with equal probability, 1/6, of the first five outcomes
- If we know probability of rolling 1, 2, 3, 4, and 5, we know the probability of rolling a 6

Cumulative logit link

We still use a logit link function but it’s a bit more complicated
Example: we have four categories, 1, 2, 3, 4
- Estimate probability of category 1, probability of 2 or less, and probability of 3 or less
- (Probability of 4 or less = 1)
These are cumulative probabilities so we’ll fit a cumulative logit model

Link functions in a cumulative logit model

$n - 1$ link functions, one for each category except the last
Map our cumulative probabilities for each category to the linear predictor scale
If probabilities are $p_{y \le 1}$, $p_{y \le 2}$, and $p_{y \le 3}$, the link functions are:
- $\log \frac {p_{y \le 1}}{1-p_{y \le 1}}$
- $\log \frac {p_{y \le 1} + p_{y \le 2}}{1-(p_{y \le 1} + p_{y \le 2})}$
- $\log \frac {p_{y \le 1} + p_{y \le 2} + p_{y \le 3}}{1-(p_{y \le 1} + p_{y \le 2} + p_{y \le 3})}$
Estimate fixed and random effects not just for one outcome, but for multiple outcomes: the cumulative probability of each class

Example dataset: turfgrass color ratings

karcher.turfgrass from agridat package
Effect of management style and nitrogen level on turfgrass color
Completely randomized design
- Four management levels are different ways of applying N (at the surface or injected at different depths)
- Two nitrogen levels: low and high N addition

Experimental plots rated for four weeks for their color quality
Each week each combination within each replicate was rated a variable number of times
Rating scale: Poor < Average < Good < Excellent

Import and examine example dataset

data(karcher.turfgrass)
head(karcher.turfgrass)
table(karcher.turfgrass$rating)

Process example data

Pool Good and Excellent to “Good+” and use ordered() to set ascending order
Data for a brms multinomial model need to be in long format (single outcome per row) instead of counts in each row
Convert week to a discrete factor variable

karcher_long <- karcher.turfgrass %>%
  mutate(rating = if_else(rating %in% c('Good', 'Excellent'), 'Good+', rating) %>% ordered(levels = c('Poor', 'Average', 'Good+'))) %>%
  group_by(week, rep, manage, nitro) %>%
  reframe(rating = rep(rating, count)) %>%
  mutate(week = factor(week))

Plot example data

ggplot(karcher_long, aes(x = interaction(nitro, manage))) +
  geom_bar(aes(fill = rating), stat = 'count', position = 'stack') +
  facet_wrap(~ week) +
  labs(x = 'Treatment combination')

Difference between management treatments?
Interaction between management and nitrogen?
Change over time?

Fit the model

rating ~ nitro * manage + week
- Model interactive effect of nitro and manage on cumulative probability of each rating class
- Include week as categorical fixed effect
(1 | rep:nitro:manage)
- Account for repeated measures over time within plots
family = cumulative(link = 'logit', threshold = 'flexible')
- Fit a cumulative logit model
- flexible means allow the threshold between each consecutive category to have a different intercept
- So each link function has different intercept but the same fixed and random effects
Weakly informative normal(0, 2) prior on fixed effects (reasonable for log-odds scale)

karcher_clmm <- brm(rating ~ nitro * manage + week + (1 | rep:nitro:manage),
                    family = cumulative(link = 'logit', threshold = 'flexible'),
                    prior = c(
                      prior(normal(0, 2), class = b)
                    ),
                    data = karcher_long,
                    seed = 928, file = 'fits/karcher_clmm')

Diagnostics

Posterior predictive check

pp_check(karcher_clmm, type = 'bars_grouped', group = 'manage')
pp_check(karcher_clmm, type = 'bars_grouped', group = 'nitro')

Quick way to look at model predictions (effect of management within each level of nitrogen)

conditional_effects(karcher_clmm, categorical = TRUE, effects = 'manage', conditions = data.frame(nitro = c('N1', 'N2')))

Model predictions

In previous lessons we have used emmeans() to estimate and compare treatment marginal means
Unfortunately emmeans does not fully support family = cumulative in brms
We will need to write our own code to get the estimates we want

Step 1: get expected values of posterior predictive

Create prediction grid with all treatment and time combinations
Generate expected value of posterior predictive distribution for each posterior draw for each combination
- Random effects are set to 0 with re_formula = ~ 0, just as emmeans() does
Name category column 'rating, posterior estimate column called '.epred' by default

pred_grid <- with(karcher_long, expand_grid(nitro = unique(nitro), manage = unique(manage), week = unique(week)))

karcher_epred <- pred_grid %>%
  add_epred_draws(karcher_clmm, re_formula = ~ 0, category = 'rating')

Get predictions averaged over time
Call this 'overall' in the week column and combine with the other predictions

karcher_epred_avgtime <- karcher_epred %>%
  group_by(nitro, manage, .draw, rating) %>%
  summarize(.epred = mean(.epred)) 

karcher_epred <- bind_rows(karcher_epred, karcher_epred_avgtime) %>%
  mutate(week = forcats::fct_na_value_to_level(week, level = 'overall'))

Step 2. Make a plot of the model predictions

ggplot(karcher_epred, aes(x = interaction(nitro, manage), y = .epred, group = rating)) +
  stat_interval(aes(color = rating, color_ramp = after_stat(level)), .width = c(.5, .8, .95), position = position_dodge(width = .5)) + 
  stat_summary(fun = 'median', geom = 'point', position = position_dodge(width = .5)) +
  facet_wrap(~ week, ncol = 2) +
  colorblind_friendly +
  labs(x = 'Treatment combination', color_ramp = 'CI width')

Mean class prediction: a sketchy assumption?

Use modeled probabilities of each category to estimate the mean predicted category for each treatment combination
This does assume that decimal values between the scores are meaningful
But at least the underlying model was fit assuming they were discrete ordered categories
Weighted average of class values, weighted by probabilities:
- Example: assume Poor = 1, Average = 2, and Good+ = 3
- For one posterior draw, we have P(Poor) = 0.05, P(Average) = 0.20, and P(Good+) = 0.75
- Estimated mean class is 0.05 * 1 + 0.20 * 2 + 0.75 * 3 or 2.70.

karcher_epred_meanclass <- karcher_epred %>%
  group_by(nitro, manage, week, .draw) %>%
  summarize(mean_class = weighted.mean(x = as.numeric(rating), w = .epred))

Plot mean predicted classes

ggplot(karcher_epred_meanclass, aes(x = interaction(nitro, manage), y = mean_class)) +
  stat_interval(.width = c(.5, .8, .95)) + 
  stat_summary(fun = 'median', geom = 'point') +
  facet_wrap(~ week, ncol = 2) +
  scale_color_brewer(palette = 'Blues') +
  labs(x = 'Treatment combination', y = 'Predicted mean class', color = 'CI width')

Step 3. Testing whether mean class differs between treatment combinations

In the past we’ve used contrast() from emmeans to do this but now we have to do it longhand
Custom function to take all pairwise differences by group
Apply the function to contrast mean predicted classes, averaged over weeks
- Difference between the two nitrogen levels within each management level
- Difference between each pair of management levels within each nitrogen level
Other comparisons are possible

pairwise_diffs <- function(values, group_names, reverse = TRUE) {
  if (reverse) {
    values <- rev(values)
    group_names <- rev(group_names)
  }
  diffs <- apply(combn(values, 2), 2, function(x) x[1] - x[2])
  contrast_names <- apply(combn(as.character(group_names), 2), 2, paste, collapse = '-')
  data.frame(contrast = contrast_names, value = diffs)
}

karcher_mgmt_contrasts <- karcher_epred_meanclass %>% 
  filter(week == 'overall') %>%
  group_by(nitro, .draw) %>%
  reframe(pairwise_diffs(mean_class, manage, reverse = TRUE))

karcher_nitro_contrasts <- karcher_epred_meanclass %>% 
  filter(week == 'overall') %>%
  group_by(manage, .draw) %>%
  reframe(pairwise_diffs(mean_class, nitro, reverse = TRUE))

Making inference based on the contrasts

describe_posterior() generates table of medians and credible intervals for contrasts
$p_d$ and $p_{MAP}$: “Bayesian p-values”
Note we have 4000 posterior samples so pd = 1 is really pd > 3999/4000
$p_{MAP}$ is not as strictly limited but you do need a lot of samples to characterize the tail of a distribution

(karcher_mgmt_contrast_post <- karcher_mgmt_contrasts %>%
  group_by(nitro, contrast) %>%
  summarize(describe_posterior(value, ci_method = 'eti', ci_width = c(.66, .95), test = c('pd', 'p_map'))))

(karcher_nitro_contrast_post <- karcher_nitro_contrasts %>%
  group_by(manage, contrast) %>%
  summarize(describe_posterior(value, ci_method = 'eti', ci_width = c(.66, .95), test = c('pd', 'p_map'))))

Presentation ideas: How to report this analysis in a manuscript?

Description of statistical model for methods section (including model comparison results)
Verbal description of results
Figure

Methods section

We fit a generalized linear mixed model to the turfgrass color rating data, pooling Good and Excellent to a single category for a total of three categories. The model assumed a multinomial response distribution with cumulative logit link function. We allowed the intercept of each linear predictor to vary by rating category but estimated a single set of fixed effects (nitrogen treatment, management treatment, week as a discrete variable, and all their interactions) and random effects (random intercept for replicate nested within treatment combination) for each category threshold. We put weakly informative Normal(0, 2) priors on all fixed effects. Using Hamiltonian Monte Carlo, we sampled 1000 warmup iterations (discarded) and 1000 sampling iterations for each of four chains.

Methods section, continued

We generated posterior distributions of the estimated marginal probability of each rating category, and the mean category prediction, for each combination of treatment and week and averaged over weeks. We compared the mean category prediction for each pair of treatments averaged over weeks, testing against a null hypothesis of no difference using the Bayesian probability of direction $p_d$ and maximum a posteriori p-value $p_{MAP}$. We characterized the posterior distributions using the medians and equal-tailed credible intervals.

Results section

The probability of higher color ratings tended to be higher in the M3 and M4 management treatments relative to M1 and M2, and higher in the high N addition treatment relative to the low N addition treatment. There were only minimal interactions between management treatment and nitrogen addition level, and the relative differences between the treatments varied very little over time.

(At this point I would refer to specific tables and figures because it would become cumbersome to repeat many probabilities in the text).

Manuscript-quality figure

Only present means averaged over weeks (as we ignored interactions with time anyway)
Present the probabilities instead of the (sketchy) mean class predictions
Separate into a grid of panels by nitrogen treatment and rating class

ggplot(karcher_epred %>% filter(week == 'overall'), aes(x = manage, y = .epred)) +
  stat_interval(aes(color = rating, color_ramp = after_stat(level)), .width = c(.5, .8, .95)) + 
  stat_summary(fun = 'median', geom = 'point') +
  facet_grid(rating ~ nitro,
             labeller = labeller(
               nitro = c('N1' = 'low N addition', 'N2' = 'high N addition')
             )) +
  colorblind_friendly +
  labs(x = 'Management treatment', y = 'Modeled class probability', color_ramp = 'CI width')

Part 2. Nonlinear mixed models

Nonlinear model versus GAM

In Lesson 3 we used GAMs (semiparametric models) to model spatial variation
GAMs use a bunch of polynomials to approximate a nonlinear relationship with linear models
But what if we have a strong theoretical expectation about how two variables are related?
We can fit a fully parametric nonlinear mixed model with brms
We can use any error distribution we want and put fixed and/or random effects on any of the parameters

Microbial growth curves

A population of microbial cells grows exponentially at an increasing rate until they start to crowd each other and growth levels off
This looks like an S-curve
Different mathematical equations exist to approximate this relationship
Based on existing knowledge and theory we can hypothesize which functional form will fit the data best

Gompertz equation

\[f(t) = a e^{-b e^{-ct}}\]

One of many equations that looks kind of like an S
Has been shown to model microbial biomass growth well
Population size as a function of time $t$, with three parameters:
- $a$: asymptotic population size, or carrying capacity of environment
- $b$: parameter to offset the curve to left or right along time axis (length of lag phase)
- $c$: determines growth rate (how fast population grows toward carrying capacity)

Example Gompertz curves

Example data

Small subset of a real ARS dataset, with noise added
How much do pyrrocidines inhibit the growth of Fusarium mold?
Bioscreen: “microbial growth curve analysis system”
- Machine with a growth plate inside, which has many wells
Three replicate plates (I am only giving you a subset of the treatments on each plate here)
Three compounds: Pyrrocidine A, Pyrrocidine B, and an equal mixture of A+B
Fusarium growth measured in the presence of each compound, at regular time intervals

Example data, continued

Not possible to directly measure biomass
Reader shines a light through each well and measures how much light the stuff in the well blocks
Optical density (OD) correlates very closely with the amount of biomass in the well

Import data

Read from CSV in the data folder on cloud server

bioscreen <- read_csv('https://usda-ree-ars.github.io/SEAStats/brms_crash_course/bioscreen.csv')

If not on cloud server, read from GitHub URL

bioscreen <- read_csv('data/bioscreen.csv')

Visualize data

bioscreen <- bioscreen %>%
  mutate(across(c(rep, well, compound), factor))

head(bioscreen)

ggplot(bioscreen, aes(x = HOUR, y = OD, color = compound, group = interaction(compound, rep, well))) +
  geom_line(alpha = 0.25, linewidth = 0.5) + 
  stat_summary(aes(group = compound), fun = 'median', geom = 'line', linewidth = 1.5) +
  scale_y_log10() +
  colorblind_friendly

Note logarithmic scale on y-axis

Adjust OD values

bioscreen <- bioscreen %>%
  group_by(well, rep) %>%
  mutate(adjusted_log_OD = log(OD) - min(log(OD)))

Gompertz equation requires positive response values
log(OD) is negative for many values in our raw data
Transform data by taking the logarithm and adjusting so that the minimum value in each well is 0

ggplot(bioscreen, aes(x = HOUR, y = adjusted_log_OD, color = compound, group = interaction(compound, rep, well))) +
  geom_line(alpha = 0.25, linewidth = 0.5) + 
  stat_summary(aes(group = compound), fun = 'median', geom = 'line', linewidth = 1.5) +
  colorblind_friendly

Determine priors for nonlinear parameters

nls(adjusted_log_OD ~ A * exp(-B * exp(-C * HOUR)), data = bioscreen, start = list(A = 3, B = 1, C = 0.1))

Fit frequentist nonlinear least squares model to the data, ignoring all fixed and random effects
Use the parameters of this model for means of the prior distributions of the intercepts for the parameters of the curve in our Bayesian model

Fixed and random effects on nonlinear parameters

Just like slopes and intercepts in a mixed model can have fixed and random effects, so can any of the nonlinear parameters in our model
Use bf() function with nl = TRUE to create multi-line formula
Each nonlinear parameter gets its own formula defining its fixed and random effects

Nonlinear formula with `bf()`

First line is the nonlinear formula
- Here it is Gompertz equation adjusted_log_OD ~ a * exp(-b * exp(-c * HOUR))
Next we allow asymptote a to vary as a function of compound (fixed) and well nested within rep (random)
- a ~ compound + (1 | well:rep)
Fit a single b and c parameter across all compounds and wells
- b ~ 1 and c ~ 1

Model syntax

bioscreen_gompertz_fit <- brm(
  bf(
    adjusted_log_OD ~ a * exp(-b * exp(-c * HOUR)),
    a ~ compound + (1 | well:rep),
    b ~ 1,
    c ~ 1,
    nl = TRUE
  ),
  prior = c(
    prior(normal(0, 1), class = b, nlpar = a),
    prior(normal(3, 1), class = b, nlpar = a, coef = Intercept),
    prior(normal(20, 5), class = b, nlpar = b, coef = Intercept),
    prior(normal(0, 1), class = b, nlpar = c, coef = Intercept)
  ),
  data = bioscreen,
  iter = 3500, warmup = 1000, chains = 4,
  seed = 950, file = 'fits/bioscreen_gompertz_fit'
)

Assume normally distributed error (default family)
Normal priors on each parameter’s intercept with means determined previously
normal(0, 1) prior on fixed effects on a (conservative given range of log OD)
All priors must include nlpar argument
Total of 10000 sampling iterations

Posterior predictive check plot

pp_check(bioscreen_gompertz_fit)

Model captures bimodal distribution (expected from S curve)
Mismatch may be due to our assumption of normally distributed error

Increase model complexity

bioscreen_gompertz_fit_morepars <- brm(
  bf(
    adjusted_log_OD ~ a * exp(-b * exp(-c * HOUR)),
    a ~ compound + (1 | well:rep),
    b ~ compound,
    c ~ compound,
    nl = TRUE
  ),
  prior = c(
    prior(normal(0, 1), class = b, nlpar = a),
    prior(normal(3, 1), class = b, nlpar = a, coef = Intercept),
    prior(normal(0, 1), class = b, nlpar = b),
    prior(normal(20, 5), class = b, nlpar = b, coef = Intercept),
    prior(normal(0, 1), class = b, nlpar = c),
    prior(normal(0, 1), class = b, nlpar = c, coef = Intercept)
  ),
  data = bioscreen,
  iter = 3500, warmup = 1000, chains = 4,
  seed = 951, file = 'fits/bioscreen_gompertz_fit_morepars'
)

Allow b and c to vary with compound also
Only fixed effects on b and c
normal(0, 1) priors on fixed effects

Compare the two models

Posterior predictive check plot for more complex models

pp_check(bioscreen_gompertz_fit_morepars)

Compare the two models using leave-one-out cross-validation

bioscreen_gompertz_fit <- add_criterion(bioscreen_gompertz_fit, 'loo')
bioscreen_gompertz_fit_morepars <- add_criterion(bioscreen_gompertz_fit_morepars, 'loo')

loo_compare(bioscreen_gompertz_fit, bioscreen_gompertz_fit_morepars)

Model predictions

bioscreen_pred <- expand_grid(
  HOUR = 0:120,
  compound = unique(bioscreen$compound)
) %>%
  add_epred_draws(bioscreen_gompertz_fit_morepars, re_formula = ~ 0, ndraws = 1000, seed = 1)

Create prediction grid with combinations of time and compound
add_epred_draws() to get expected value of posterior predictive
- re_formula = ~ 0 sets random effects to zero
- ndraws = 1000 only gets 1000 of the 10000 posterior draws (plenty for a quick plot)

Plot model predictions

ggplot(bioscreen_pred, aes(x = HOUR, y = .epred)) +
  stat_lineribbon(aes(fill = compound, color = compound, fill_ramp = after_stat(level)), alpha = 0.5) +
  colorblind_friendly
  
ggplot(bioscreen_pred, aes(x = HOUR, y = .epred, fill = compound, color = compound)) +
  geom_line(aes(y = adjusted_log_OD, group = interaction(compound, rep, well)), data = bioscreen, alpha = 0.25, linewidth = 0.5) +
  stat_lineribbon(aes(fill_ramp = after_stat(level)), alpha = 0.5) +
  colorblind_friendly

Second plot is messier but shows the reality of the data better

Estimate marginal means of the nonlinear parameters

Test hypothesis that any of the parameters’ means differs between compounds
emmeans() supports nonlinear brms models using nlpar argument
pairwise ~ compound estimates means and takes pairwise differences between compounds

emm_a <- emmeans(bioscreen_gompertz_fit_morepars, pairwise ~ compound, nlpar = 'a')
emm_b <- emmeans(bioscreen_gompertz_fit_morepars, pairwise ~ compound, nlpar = 'b')
emm_c <- emmeans(bioscreen_gompertz_fit_morepars, pairwise ~ compound, nlpar = 'c')

Inference from posterior distributions of parameter means

emm_a_post <- describe_posterior(emm_a$emmeans, ci_method = 'eti', ci = c(0.5, 0.8, 0.95), test = NULL)
emm_b_post <- describe_posterior(emm_b$emmeans, ci_method = 'eti', ci = c(0.5, 0.8, 0.95), test = NULL)
emm_c_post <- describe_posterior(emm_c$emmeans, ci_method = 'eti', ci = c(0.5, 0.8, 0.95), test = NULL)

contr_a_post <- describe_posterior(emm_a$contrasts, ci_method = 'eti', ci = c(0.5, 0.8, 0.95), test = c('pd', 'p_map'))
contr_b_post <- describe_posterior(emm_b$contrasts, ci_method = 'eti', ci = c(0.5, 0.8, 0.95), test = c('pd', 'p_map'))
contr_c_post <- describe_posterior(emm_c$contrasts, ci_method = 'eti', ci = c(0.5, 0.8, 0.95), test = c('pd', 'p_map'))

Equal-tailed quantile credible intervals of means in $emmeans objects
Credible intervals and Bayesian hypothesis tests in $contrasts objects

How do parameters differ between compounds? (means)

emm_a_post
emm_b_post
print(emm_c_post, digits = 4)

contr_a_post
contr_b_post
print(contr_c_post, digits = 4)

Presentation ideas: How to report this analysis in a manuscript?

Description of statistical model for methods section (including model comparison results)
Verbal description of results
Figure

Methods section

We fit a nonlinear mixed model to the OD data. The response variable was adjusted by taking the logarithm and subtracting the minimum value of each well so that the lowest value in all wells was zero. We assumed the growth rate of OD with time followed a Gompertz relationship $f(t) = a e^{-b e^{-ct}}$. A fixed effect for compound was put on each of the three nonlinear parameters: $a$, the asymptote; $b$, the offset, and $c$, the growth rate. In addition, we put a random intercept for each well nested within replicate on the asymptote parameter. We put weakly informative Normal(0, 1) priors on the fixed effects and priors for the intercepts of each nonlinear parameter which were normally distributed with means derived from the point estimates from a least squares nonlinear fit.

Methods section, continued

We sampled the posterior distribution of model parameters using the Hamiltonian Monte Carlo algorithm with four chains, each run for 1000 warmup iterations (discarded) and 2500 sampling iterations (retained). We compared this model with a similar model that did not allow the $b$ and $c$ parameters to vary by compound, and found the more complex model to be a better fit; we present predictions from only that model in what follows.

We obtained posterior estimates of the marginal mean for each nonlinear parameter for each compound and compared the estimates between each pair of compounds. We characterized the posterior distributions using the medians and equal-tailed credible intervals, and tested the contrasts against a null value of zero using the Bayesian probability of direction ($p_d$) and maximum a posteriori p-value ($p_{MAP}$).

Results section

The model predicted a similar asymptote value for adjusted log OD for each of the three compounds (A: 3.06, 95% CI [2.97, 3.15]; B: 3.19 [3.10, 3.28]; AB: 3.15 [3.06, 3.24]). The model predicted a slightly higher growth rate for B (0.063 [0.062, 0.065]) compared to A (0.061 [0.059, 0.063]) and AB (0.059 [0.058, 0.061]). This indicates that the A compound may have been slightly more effective at initially slowing the growth rate but that the long-term inhibition differs very little between the compounds.

Manuscript-quality figure

Previous figure recreated with more publication-quality options

ggplot(bioscreen_pred, aes(x = HOUR, y = .epred, fill = compound, color = compound)) +
  geom_line(aes(y = adjusted_log_OD, group = interaction(compound, rep, well)), data = bioscreen, alpha = 0.25, linewidth = 0.5) +
  stat_lineribbon(aes(fill_ramp = after_stat(level)), alpha = 0.5) +
  colorblind_friendly +
  labs(fill_ramp = 'CI width', x = 'time (h)', y = 'adjusted log(OD)') +
  theme(
    panel.grid = element_blank(),
    legend.position = 'inside',
    legend.position.inside = c(0, 1),
    legend.justification = c(0, 1),
    legend.background = element_blank()
  )

Conclusion

You now have a lot of techniques in your Bayesian toolkit! Today you learned how to:

Fit a statistical model when you have ordered categorical response data
Make different kinds of predictions and comparisons using a Bayesian cumulative logit-link mixed model
Fit a nonlinear model with fixed and random effects
Compare parameters of a nonlinear model between treatment groups in a Bayesian framework

Well done!

See text version of lesson for further reading and useful resources
Please provide feedback to help me improve this course!
- See course homepage for link to feedback form

A crash course in Bayesian mixed models with brms (Lesson 5)

Introduction and review

What is this class?

How to follow the course

What we’ve learned in Lessons 1-4

What will we learn in this lesson?

Big-picture concepts

Part 1. Cumulative logistic mixed models for ordinal data

Load packages and set options

Categorical data

Ordered categorical data

Shortcut at your own risk

Review: modeling binary data

Binomial distribution for binary data

Multinomial distribution, from coins to dice

Cumulative logit link

Link functions in a cumulative logit model

Example dataset: turfgrass color ratings

Import and examine example dataset

Process example data

Plot example data

Fit the model

Diagnostics

Model predictions

Step 1: get expected values of posterior predictive

Step 2. Make a plot of the model predictions

Mean class prediction: a sketchy assumption?

Plot mean predicted classes

Step 3. Testing whether mean class differs between treatment combinations

Making inference based on the contrasts

Presentation ideas: How to report this analysis in a manuscript?

Methods section

Methods section, continued

Results section

Manuscript-quality figure

Part 2. Nonlinear mixed models

Nonlinear model versus GAM

Microbial growth curves

Gompertz equation

Example Gompertz curves

Example data

Example data, continued

Import data

Visualize data

Adjust OD values

Determine priors for nonlinear parameters

Fixed and random effects on nonlinear parameters

Nonlinear formula with bf()

Model syntax

Posterior predictive check plot

Increase model complexity

Compare the two models

Model predictions

Plot model predictions

Estimate marginal means of the nonlinear parameters

Inference from posterior distributions of parameter means

How do parameters differ between compounds? (means)

Presentation ideas: How to report this analysis in a manuscript?

Methods section

Methods section, continued

Results section

Manuscript-quality figure

Conclusion

Nonlinear formula with `bf()`