What is this class?

• Part 2 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

Concepts we covered in Lesson 1

• What is Bayesian inference
• Prior, likelihood, and posterior
• Credible intervals

Techniques we covered in Lesson 1

• Use Hamiltonian Monte Carlo method, implemented in Stan/brms, to fit Bayesian models
• Fit a multilevel (mixed) model with random intercepts and random slopes
• Specify different priors
• Diagnose and deal with convergence problems
• Plot model parameters and predictions
• Compare models with information criteria
• Do hypothesis testing with Bayesian analogues of p-values

What will we learn in this lesson?

• Specify strongly and weakly informative priors and see how they may influence your results
• Fit a Bayesian generalized linear model (GLM), for situations where you can’t assume normally distributed error in your model
• Fit a Bayesian generalized linear mixed model (GLMM), to accommodate non-normal data and a mixture of fixed and random effects

Back to Bayes-ics: Bayes’ Theorem

\[P(A|B) = \frac{P(B|A)P(A)}{P(B)}\]

The probability of A being true given that B is true (\(P(A|B)\)) is equal to the probability that B is true given that A is true (\(P(B|A)\)) times the ratio of probabilities that A and B are true (\(\frac{P(A)}{P(B)}\))

Bayes’ theorem in a statistical context

\[P(model|data) = \frac {P(data|model)P(model)}{P(data)}\]

or

\[posterior \propto likelihood \times prior\]

Bayes’ theorem in statistics, explained

• We have some prior knowledge about the world, formalized with a statistical distribution of what we think the parameters are without having collected any data yet
  • prior, represented by \(P(model)\), or the probability of the model not dependent on the data
• We go out and collect some data and calculate \(P(data|model)\)
  • likelihood, a function that just depends on the data.
  • How likely it was to observe the data we did, given all possible combinations of parameters we could have in our model

• Divide the likelihood, \(P(data|model)\), by the marginal probability, \(P(data)\), to make the likelihood function have an area under the curve of 1 and become a probability distribution
• Multiply the prior and marginalized likelihood together to get \(P(model|data)\)
  • posterior, our new knowledge of the world, after collecting the data, represented by a new statistical distribution for our parameter estimates

The posterior distribution

Priors have a bad reputation

Argument 1: Assumptions are not a bad thing

• Prior distributions are not the only assumption; likelihood always involves assumptions too!
  • Example: a general linear model, Bayesian or not, assumes x-y relationship is linear, and the noise around that is normally distributed
  • That’s a strong assumption that only roughly approximates the real world
  • If you won’t make assumptions, you can’t fit any statistical model at all!

Argument 2: Priors are (often) more practical than philosophical

• If you have a large dataset, the prior distribution usually has relatively little influence on the results
• The prior makes the computations needed to fit the statistical model a little bit easier
• In those cases, only the Bayesian model can give us an answer in the first place

Posterior distributions get closer to likelihood with increasing sample size

• As the amount of data goes up, evidence that the true immunity rate is 70% becomes stronger
• The influence of our prior distribution centered at 50% keeps getting weaker and weaker
• The influence of the likelihood (informed by the data) keeps getting stronger
• The more data you have, the more evidence you have, and the more the likelihood pulls the posterior towards itself and away from the prior

Hands-on example of prior distribution influence

• If you have lots of data, you don’t have to worry that much about the influence of the prior
• But what if you have a very small dataset?
• Let’s do a hands-on example using data on a controversial topic where people may have strong prior beliefs

Setup

Load packages.

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

• Set plotting options (including colorblind-friendly palette)
• Create directory for fitted model objects if it does not exist

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

if (!dir.exists('fits')) dir.create('fits')

• Set brms options
• Cloud server has pre-fit model objects loaded to save time

options(mc.cores = 4, brms.file_refit = 'on_change', brms.backend = 'cmdstanr')

• If you are on a local machine without CmdStanR, do not set brms.backend

options(mc.cores = 4, brms.file_refit = 'on_change')

Meet the example data

• Small sample size: \(n=8\) Thysanoptera abundance measurements on each line

data(gathmann.bt)

ggplot(gathmann.bt, aes(x = gen, y = thysan)) + geom_boxplot()

• Set factor order so that the control line (ISO) is treated as the reference level in the model
• A positive coefficient in the model will mean that Bt has higher Thysanoptera abundance than control

gathmann.bt <- mutate(gathmann.bt, gen = relevel(gen, ref = 'ISO'))

Frequentist analysis: two-sided t-test

t.test(formula = thysan ~ gen, data = gathmann.bt)

• We ask: is there enough evidence to reject the null hypothesis that the two means are equal?
• Only considers the likelihood, which is calculated from the data
• Does not formally incorporate any prior knowledge

Bayesian analyses with different priors

• The response variable ranges from ~4 to 18, so an effect of +5 or -5 is “big”
• We’re going to fit the model with 7 different priors
• Default completely uninformative “flat” prior assigns equal prior probability to any effect size
• Not realistic for most situations
• We usually at least know the order of magnitude of expected effect size, even if we don’t know if it will be positive or negative

Six other possible prior distributions

Comparison of posterior estimates

• Pull out the posterior estimates for the difference between the means
• Plot the median of the posterior with 95% credible interval

get_post <- function(fit) {
  fixedeff <- summary(fit)$fixed
  setNames(fixedeff['genBt', c('Estimate','l-95% CI', 'u-95% CI')], c('median', 'lower95', 'upper95'))
}

gathmann_post_estimates <- cbind(
  prior = c('flat', 'normal(0,5)', 'normal(0,1)', 'normal(5,5)', 'normal(5,1)', 'normal(-5,5)', 'normal(-5,1)'),
  map_dfr(list(m1gathmann, m2gathmann, m3gathmann, m4gathmann, m5gathmann, m6gathmann, m7gathmann), get_post)
) %>%
  mutate(prior = factor(prior, levels = unique(prior)))

ggplot(gathmann_post_estimates, aes(x = prior, y = median, ymin = lower95, ymax = upper95)) +
  geom_pointrange() +
  geom_hline(yintercept = 0, linetype = 'dashed', linewidth = 1) +
  labs(y = 'effect of Bt on Thysanoptera abundance')

Comparison of posterior estimates

• \(\text{Normal}(0,5)\), \(\text{Normal}(5,5)\), and \(\text{Normal}(-5,5)\): all similar inference to the flat prior
• \(\text{Normal}(5, 1)\): similar point estimate to flat prior, but narrower uncertainty interval (informative prior)
• \(\text{Normal}(0, 1)\): informative prior that assigns a very low prior probability to strong positive and strong negative effects
• \(\text{Normal}(-5, 1)\): informative prior that assigns low prior probability to everything except strong negative effects
• The informative priors with low standard deviations influence our conclusions

Comparison of posterior estimates: Bayesian-style p-values

• “maximum a posteriori” p-value, or p_MAP (see Lesson 1)
• Ratio of posterior probability of null value to probability of posterior point estimate
• The lower p_MAP is, the stronger the evidence in favor of the posterior estimate relative to null

gathmann_pvalues <- data.frame(
  prior = c('flat', 'normal(0,5)', 'normal(0,1)', 'normal(5,5)', 'normal(5,1)', 'normal(-5,5)', 'normal(-5,1)'),
  pMAP = c(p_map(m1gathmann)$p_MAP[2], p_map(m2gathmann)$p_MAP[2], p_map(m3gathmann)$p_MAP[2], p_map(m4gathmann)$p_MAP[2], p_map(m5gathmann)$p_MAP[2], p_map(m6gathmann)$p_MAP[2], p_map(m7gathmann)$p_MAP[2])
)

• So far we have fit simple models where Bayesian and frequentist methods both work well
• Let’s explore some territory where Bayesian methods really shine!

Review: general linear mixed models

\[y \sim \text{Normal}(Xb + Zu, \sigma), u \sim \text{Normal}(0, \sigma_u)\]

• This means response \(y\) is a function of fixed effects and random effects
  • \(X\): the predictor variables
  • \(b\): the fixed effect parameters
  • \(Z\): model matrix that says which group (block) each observation belongs to
  • \(u\): the random effect parameters
  • \(\sigma\): normally distributed noise

Review: general linear mixed models

\[y \sim \text{Normal}(Xb + Zu, \sigma), u \sim \text{Normal}(0, \sigma_u)\]

• Random effects \(u\) are normally distributed too
• Some random effects are positive, others are negative, a few are far from zero, most are close to zero
• Same goes for the model residuals

Non-normal data: making predictions

• Predictions made from a general linear mixed model fit to non-normal data often don’t make sense
  • Model fit to binary data may give you predicted probabilities greater than 1 or less than 0
  • Model fit to the skewed data may give you negative predictions
• Even if you don’t care about prediction but only about inference, this is still a problem
• If your model doesn’t make good predictions, it doesn’t fit the data well, and you need a good fit to make inference!

Dealing with non-normal data: transforming the data

• Old school solution: transform the data!
• Not ideal because it introduces bias (possible exception: log transformation)
• We should find a model that fits the data, instead of torturing the data to make it fit the model
• Better alternative: generalized linear models!

Generalized linear mixed models: an introduction

\[g(y) \sim \text{Distribution}(Xb + Zu | \theta), u \sim \text{Normal}(0, \sigma_u)\]

• The distribution can be any of a wide variety of distributions, with parameters \(\theta\)
• The left side now says \(g(y)\) instead of just \(y\): link function
  • Gives us the relationship between the distribution function and the predicted variable
  • Converts variables from the scale of the statistical distribution, to the scale of the data

What is a generalized linear mixed model?

• generalized: it allows for response distributions that aren’t normal
• linear: Y is modeled as a linear combination of effects (multiply each variable by something, then add them together)
• mixed: it includes both fixed and random effects
• This example will be generalized and linear, but not mixed (fixed effects only)

Experimental design

• Many plates, each with a variable number of seeds from one treatment combination
• Completely randomized design: no blocks, so no random effects
  • n: total number of seeds on the plate
  • germ: the number of seeds that germinated

• Compute and plot p (proportion of germinated seeds)

data(crowder.seeds)

crowder.seeds <- mutate(crowder.seeds, p = germ/n)

ggplot(crowder.seeds, aes(x = gen, y = p, color = extract)) +
  geom_boxplot() + colorblind_friendly

Fit a linear model to the proportions

• This is OK in many cases if proportions are between ~0.05 and ~0.95
• Assumption of normally distributed error can still work out
• Because proportions are between 0 and 1, \(\text{Normal}(0, 1)\) prior on categorical fixed effect is completely uninformative (it encompasses any possible effect size)

crowder_lm <- brm(p ~ gen * extract, 
                  prior = c(
                    prior(normal(0, 1), class = b)
                  ),
                  data = crowder.seeds, 
                  seed = 111, file = 'fits/crowder_lm')

Fit a generalized linear model to the binary outcomes

• Outcome is binary (0 = seed did not germinate, 1 = seed germinated)
• Response variable is specified with two parts
  • Number of successes (germinated seeds in each dish)
  • Total number of trials (all seeds in each dish including the ones that didn’t germinate)
• germ | trials(n)

The logit link function

• Linear combination of covariates in a linear model \(Xb + Zu\) can take any value, positive or negative
• But we are trying to model the probability of germination, which can only be between 0 and 1
• Logit link function is used to transform the unbounded value to a scale that ranges from 0 to 1

GLM code

• normal(0, 5) prior is sensible for fixed effects on the log odds scale
• Again we will use all default model fitting options

crowder_glm <- brm(germ | trials(n) ~ gen * extract, 
                   family = binomial(link = 'logit'),
                   prior = c(
                     prior(normal(0, 5), class = b)
                   ),
                   data = crowder.seeds, 
                   seed = 112, file = 'fits/crowder_glm')

Comparing LM and GLM

What do you notice about the posterior predictive check plots?

• thick black line: distribution of the observed data
• each thin blue line: distribution of the predicted values from one random draw from the posterior

pp_check(crowder_lm)
pp_check(crowder_glm)

LM and GLM: Posterior predictive check plots

• Both are not too bad: shape of the fitted values is similar to shape of the observed data
• Scale of x-axis is different
• But notice some predictions from LM are <0 and >1
• Why? Why is this a problem?

LM: Estimating marginal means

• pairwise ~ extract | gen: within each level of genotype (gen), estimate the predicted mean for each level of extract and take the pairwise contrast between them
• LM contrasts are taken by subtracting one modeled proportion from the other

emm_crowder_lm <- emmeans(crowder_lm, pairwise ~ extract | gen)

emm_crowder_lm$emmeans
emm_crowder_lm$contrasts

GLM: Estimating marginal means

• Estimated marginal means on the logit scale
• Contrasts on the log odds ratio scale

emm_crowder_glm <- emmeans(crowder_glm, pairwise ~ extract | gen)

emm_crowder_glm$emmeans
emm_crowder_glm$contrasts

Estimated marginal means on the response scale

• The argument type = 'response' transforms the results back to the data scale
• The logit link function in R is qlogis() and its inverse is plogis()

emmresponse_crowder_glm <- emmeans(crowder_glm, pairwise ~ extract | gen, type = 'response')

emmresponse_crowder_glm$emmeans
emmresponse_crowder_glm$contrasts

Odds ratios

• Now estimated marginal means and credible intervals are probabilities
• Contrasts are odds ratios
  • Ratio of the odds of germination between the two levels of extract
  • An odds ratio of 1 indicates no difference between the odds
  • The ratio is bean / cucumber so we have evidence that bean extract is less effective at stimulating germination than cucumber extract

Hypothesis testing

• What if we want to do a test relative to a null hypothesis?
• Use describe_posterior() from the bayestestR package
• Describe the posterior distributions of the differences between the means of the extract treatments within each genotype for the general linear model
• Describe the posterior distributions of the odds ratios of the means for the generalized linear model

describe_posterior() explained

• We have 1000 sampling iterations from each of 4 chains: total of 4000 posterior draws for the distribution of each contrast
• By default describe_posterior() gives us the median of the 4000 values as the point estimate
• Specify equal-tailed credible interval using ci_method = 'eti', default width 95%
  • We get 2.5% and 97.5% quantiles of the 4000 values

Comparing inference from LM and GLM

• The generalized linear model provides stronger evidence in favor of differences
• It is taking into account the greater precision we can get from having many partially independent binary outcomes (seeds) in each dish
• 200/1000 is more evidence than 2/10 (both are treated the same in the linear model fit to proportions)

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

• Description of the statistical model for the methods section
• Verbal description of the results
• Figure
• Table

Results section

The modeled probability of germination was greater when seeds were treated with cucumber extract versus bean extract for both the O73 genotype (bean extract: posterior median 0.396, 95% CI [0.316, 0.478]; cucumber extract: 0.531 [0.458, 0.614]) and the O75 genotype (bean extract: 0.364 [0.305, 0.421]; cucumber extract: 0.681 [0.629, 0.734]). The odds ratio between cucumber and bean extract was smaller for the O73 genotype than the O75 genotype (ratio 0.46 [0.25, 0.83], p_MAP < 0.001), indicating a greater relative effect of cucumber extract on germination probability in the O75 genotype.

Generalized linear mixed models

• GLMM: non-normal data and a mixture of fixed and random effects
• Bayesian GLMMs really outperform frequentist
• Model syntax requires random terms on right-hand side of formula, and to specify response family and link functione.

Meet the example data

data(mcconway.turnip)
ggplot(mcconway.turnip, aes(x = yield)) + geom_histogram()

Turnip experimental design

• 2×2×4 factorial design experiment, with three factors
  • genotype (gen, either "Barkant" or "Marco")
  • planting date (date, either "21Aug1990" or "28Aug1990")
  • planting density (density, 1, 2, 4, or 8 kg/ha)
• Randomized complete block design with four blocks
• Each block has 16 plots, one for each treatment combination
• A single measurement of yield was taken from each plot

• Treat density as a discrete variable with four levels
• This is easier because we don’t need to worry about whether the effect of planting density on yield is linear

mcconway.turnip <- mutate(mcconway.turnip, density = factor(density))

• Boxplot of the yield values in each treatment combination
• Split the genotypes into separate panels using facet_wrap()
• Arrange the boxplots in increasing order of density on the x-axis
• Show the planting dates for each density with a different color box

ggplot(mcconway.turnip, aes(x = density, y = yield, group = interaction(density, date), color = date)) +
  facet_wrap(~ gen) +
  geom_boxplot() +
  colorblind_friendly +
  labs(x = 'Planting density (kg/ha)', y = 'Yield', color = 'Planting date')

Yield distribution

• Interaction between planting date and density
• Possible interaction with genotype
• Relationship between mean and variance
• Linear model assumes one single standard deviation for all the residuals, so unequal variance can cause bad predictions

Side note: divergent transitions

• You may see a warning about divergent transitions
• In this case they are few enough that we can ignore them
• They often occur when fitting complex models to small datasets
• May indicate numerical problems with the sampler; it might be “missing” parts of the posterior distribution
• Fix this by setting narrower priors or tweaking the sampler arguments
• Make the sampler take smaller steps with each iteration, which reduces chance of divergent transition but is slower

Posterior predictive check plot on linear model

pp_check(mcconway_lmm)

• What is our model getting wrong?
• Also look at estimated marginal means by density

emmeans(mcconway_lmm, ~ density)

Posterior predictive check plot: LMM with log-transformation

pp_check(mcconway_loglmm)

Fit a GLMM with gamma response distribution

Fit a GLMM with gamma response distribution and log link function using family = Gamma(link = 'log').

mcconway_glmm <- brm(yield ~ gen * date * density + (1 | block), 
                     family = Gamma(link = 'log'), 
                     prior = c(
                       prior(normal(0, 2), class = b)
                     ),
                     data = mcconway.turnip, 
                     seed = 115, file = 'fits/mcconway_glmm')

Posterior predictive check: Gamma GLMM

pp_check(mcconway_glmm)

Estimating marginal means: Gamma GLMM

• Marginal means for each combination of fixed effects
• The means and CIs do not include random effect of block
• Expected value for a block that’s “smack-dab in the middle”
• type = 'response' gives predictions on data scale, not log link scale

emm_mcconway <- emmeans(mcconway_glmm, ~ density + date + gen, type = 'response')

emm_mcconway

Contrasts: Gamma GLMM

• For example, pairwise contrasts between density levels within each combination of planting date and genotype
• Contrasts are ratios
• 'revpairwise' so that the higher density is in the numerator

contrast(emm_mcconway, 'revpairwise', by = c('date', 'gen'))

Assessing evidence for interactions between treatments

• Effect sizes look bigger for Marco than Barkant
• Interaction contrasts (contrasts of contrasts) to test this

contrast(emm_mcconway, 'revpairwise', by = 'date', interaction = 'revpairwise')

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

• Description of the statistical model for the methods section
• Text write-up for the results section
• A figure

Manuscript-quality figure

• exp() is used to back-transform the posterior estimates

gather_emmeans_draws(emm_mcconway) %>%
  mutate(.value = exp(.value)) %>%
  ggplot(aes(x = density, y = .value, group = interaction(density, date))) +
    stat_interval(aes(color = date, color_ramp = after_stat(level)), position = position_dodge(width = .5)) +
    stat_summary(fun = 'median', geom = 'point', size = 2, position = position_dodge(width = .5)) +
    facet_wrap(~ gen) +
    colorblind_friendly +
    labs(x = 'Planting density (kg/ha)', y = 'Modeled yield', color_ramp = 'CI width', color = 'Planting date')

We’ve learned . . .

• How to explore different priors and how they may influence (or not influence) your conclusions
• How to fit a Bayesian generalized linear model (GLM) for non-normal response data
• How to fit a Bayesian generalized linear mixed model (GLMM) for non-normal response data with random effects
• Techniques from packages brms, emmeans, tidybayes, and bayestestR

Give yourself a well-deserved round of applause!

• See text version of lesson for further reading and useful resources
• Please provide feedback to help me improve this course!
  • Send feedback to quentin.read@usda.gov