What is this class?

• Part 4 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-3

• 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

What will we learn in this lesson?

• Beta regression for proportion data
• Zero-inflated Poisson regression for count data that include zeros
• Hurdle models for continuous and proportion data that include zeros

Big-picture concepts

• What are mixture distributions, and how they make stat models more flexible
• How to allow each separate parameter of a statistical distribution to vary as a function of fixed and random effects
• Reminder: not all these things have to be done Bayesian, but it’s often the only way that works in practice!

Regression with proportional data

• In Lesson 2, we modeled binary data with a logistic GLMM instead of fitting a model to the proportions
• But sometimes you just have the proportions and not the counts, or you truly have a continuous proportion not based on counts
• Normally distributed error won’t work because proportions are bounded between 0 and 1
  • It can work OK if most of your data is between ~0.25 and ~0.75, but that often isn’t the case in ag science

Beta distribution to the rescue!

• Can be a bell-like curve, a skewed hump, or even a horseshoe, but it always goes from 0 to 1

Parameters of the beta distribution

• Many sources use parameters \(\alpha\) and \(\beta\) (Wikipedia)
• R function dbeta() uses those parameters but calls them shape1 and shape2
• In brms, it is rewritten in terms of different parameters, \(\mu\) (mu) and \(\phi\) (phi)
  • mu is the mean parameter
  • phi is the precision parameter

Mean parameterization

• Most distributions in brms are parameterized so that the first parameter is the mean or expected value of the distribution
• Easy to separate effects on the mean from effects on other components of the distribution
• You can put any combination of fixed and/or random effects on any parameter of any distribution in brms

Example dataset of proportions

Split-plot design without replication

• New randomization done each year at a different site

Load packages and set options

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

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

options(mc.cores = 4, brms.file_refit = 'on_change') # If you don't have CmdStan
options(mc.cores = 4, brms.file_refit = 'on_change', brms.backend = 'cmdstanr') # If you do have CmdStan

Pre-process data

• Beta distribution requires proportions, not percentages
• Convert gen, strain, and year to factor data type

data(snijders.fusarium)

snijders.fusarium <- snijders.fusarium %>% 
  mutate(
    across(c(gen, strain, year), factor),
    p = y/100
  )

Histograms: distribution of data for each strain

• These shapes are typical for proportion data, and resemble beta distributions

ggplot(snijders.fusarium, aes(x = p, fill = strain)) + 
  geom_histogram(bins = 15) +
  facet_wrap(~ strain) +
  colorblind_friendly +
  theme(legend.position = 'none')

Strip plot: raw data and medians for each strain

• Showing all the data is often preferable to a boxplot, if the dataset is small

ggplot(snijders.fusarium, aes(x = year, y = p, color = strain)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.1, jitter.height = 0, dodge.width = 0.5), alpha = 0.5) +
  stat_summary(fun = 'median', geom = 'point', size = 15, shape = 95, position = position_dodge(width = 0.5)) +
  colorblind_friendly

Model formula and response family

• p ~ strain * year + (1 | gen) + (1 | gen:year)
  • Response p on the left
  • Fixed effects for strain (categorical), year (categorical), and their interaction
  • Random intercept for each genotype and each genotype-year combination
• family = Beta(link = 'logit', link_phi = 'log')
  • Logit link for the mean response
  • Log link for the precision parameter phi

Fit the beta regression model

fusarium_beta_fit <- brm(
  p ~ strain * year + (1 | gen) + (1 | gen:year),
  family = Beta(link = 'logit', link_phi = 'log'),
  data = snijders.fusarium,
  prior = c(
    prior(normal(0, 5), class = b)
  ),
  chains = 4, iter = 2000, warmup = 1000, seed = 555,
  file = 'fits/fusarium_beta_fit'
)

Posterior predictive check: is it a good fit?

pp_check(fusarium_beta_fit)

Note on bf() syntax

• We will revisit shortly when we learn about zero-inflated and hurdle models
• You could even use the distributional regression bf() syntax in a typical general linear model with normally distributed error
  • Standard deviation sigma can vary with fixed or random effects, if you want to model unequal variance among groups
• See the brms vignette on estimating distributional models

pp_check() and model comparison

pp_check(fusarium_beta_fit_varyphi)

fusarium_beta_fit <- add_criterion(fusarium_beta_fit, 'loo')
fusarium_beta_fit_varyphi <- add_criterion(fusarium_beta_fit_varyphi, 'loo')

loo_compare(fusarium_beta_fit, fusarium_beta_fit_varyphi)

What can we conclude?

Making predictions with the beta model

• Use the simpler model with a fixed phi
• Estimate marginal means for each combination of strain and year
  • type = 'response' means use inverse link scale (the same scale as the original data)
• Take contrasts between each pair of strains in each year
  • Use regrid(fusarium_emms, transform = 'log') to get the contrasts in the form of ratios, not odds ratios

(fusarium_emms <- emmeans(fusarium_beta_fit, ~ strain | year, type = 'response'))
(fusarium_contrasts <- contrast(regrid(fusarium_emms, transform = 'log'), 'pairwise'))

Credible intervals and hypothesis tests

• Use describe_posterior() from the bayestestR package to get credible intervals for the means and pairwise differences
  • These are quantiles of the 4000 posterior samples we generated for each of those quantities
• \(p_{MAP}\) and \(p_{d}\) are Bayesian analogues to p-values, if you are keen on testing a null hypothesis

(fusarium_emms_CI <- describe_posterior(fusarium_emms, ci = c(0.66, 0.95), ci_method = 'eti', test = NULL))
(fusarium_contrasts_CI <- describe_posterior(fusarium_contrasts, ci = c(0.66, 0.95), ci_method = 'eti', test = c('pd', 'p_map'), null = 1, as_p = TRUE))

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
• Table

Methods section, continued

We obtained posterior estimates of the marginal means for each combination of strain and year, back-transformed from the logit scale, and took pairwise ratios between strain means within each year. We characterized the posterior distributions of means and ratios using the median and equal-tailed credible intervals. We teested the null hypothesis that each ratio was equal to 1 using two different Bayesian p-value analogues: probability of direction \(p_d\) and maximum a posteriori p-value \(p_{MAP}\).

Results section

Incidence of Fusarium infection was overall low except for strain F39 in 1986 and 1987, with estimated mean infection proportions 0.24 with 95% CI [0.16, 0.34] and 0.28 [0.19, 0.40] respectively. The evidence is consistent with a greater proportion infection for F39 relative to other strains, with approximately a tenfold difference in 1987 and smaller differences in the twofold to fivefold range in 1986 and 1988 (see Table).

Manuscript-quality figure

• Modeled posterior means, credible intervals, and the raw data
• Simultaneously shows that the model is a good fit to the data, and shows the differences between the predicted means
• 66% and 95% credible intervals because that is about ± 1 and 2 SD, if the posterior distribution of the parameter is roughly normal

gather_emmeans_draws(fusarium_emms) %>%
  mutate(.value = plogis(.value)) %>%
  ggplot(aes(x = year, y = .value, group = interaction(strain, year))) +
    stat_interval(aes(color = strain, color_ramp = after_stat(level)), .width = c(.66, .95), position = ggpp::position_dodgenudge(width = .5, x = .03)) +
    stat_summary(fun = 'median', geom = 'point', size = 2, position = ggpp::position_dodgenudge(width = .5, x = .03)) +
    geom_point(aes(y = p, color = strain), data = snijders.fusarium, position = ggpp::position_dodgenudge(width = .5, x = -.03), alpha = .5) +
    colorblind_friendly +
    labs(y = 'proportion infected', color_ramp = 'CI level')

Manuscript-quality table

• Pairwise ratios of proportion infection (posterior estimates and 95% CIs), with \(p_{MAP}\) values

fusarium_contrasts_CI %>%
  filter(CI == 0.95) %>%
  dplyr::select(-c(p, CI)) %>%
  gt(groupname_col = 'year') %>%
  fmt_number(n_sigfig = 3) %>%
  cols_label(Median = 'ratio', CI_low = '95% CI lower bound', CI_high = '95% CI upper bound', p_MAP = md('p~MAP~')) %>%
  sub_zero(columns = p_MAP, zero_text = '< 0.001') %>%
  tab_options(row_group.background.color = 'slateblue3',
              column_labels.font.weight = 'bold')

Zero-inflated and hurdle models

Images (c) Philippa Willitts and Kirby Lee/USA Today Sports

• Gamma and Beta distributions only support positive values (>0)
• So far we’ve used datasets that don’t contain any zeros, so these distributions have worked very well
• But many proportion, count, and non-negative continuous datasets contain zeros, sometimes lots of them!

“Old school” way: get rid of the zeros and transform

• Force the dataset to fit in the “normal” box
• Add a little fudge factor (constant), then log transform
  • Sometimes constant is reality-based: LOD/2
  • Sometimes pulled out of thin air: log(x + 1)? log(x + 0.0000001)?
• Even without the fudge factor, a bias is introduced when transforming and back-transforming
• Or you can simply remove the zeros from the dataset before fitting a model (acceptable if very few zeros)

Better way: mixture distributions

• Mixture distribution: a statistical distribution made up of two or more distributions stacked on top of each other
• Weighted sum of their density functions so that the area under the PDF curve = 1
• You can mix distributions from different families, at different proportions

Application of mixture distributions to datasets with zeros

• We can use mixture distributions to model excess zeros (and sometimes ones) that can’t be fit by a standard statistical distrbution
• Let’s start with an example from the great book Statistical Rethinking by Richard McElreath

Images from Li livres dou santé (13th c.) and Miracles de Notre Dame (15th c.)

• Medieval monks take a long time to copy manuscripts
• Most days they work on copying manuscripts
• On some of those days, they finish one manuscript, some days two or three, some days none at all
• But on other days, they take a break and enjoy a beer, always finishing zero manuscripts

• To model this data, we don’t need to know which zeros are which
• Just assume that there is a greater proportion of zeros in the data than you would expect for a (Poisson) count distribution with a given mean

Zero-inflated models

• For discrete distributions like Poisson that already have some finite probability of being zero.
• ‘Inflate’ the probability of zero with an extra parameter

Hurdle models

• For continuous distributions like Gamma and Beta
• Continuous distributions do not have a probability of being exactly zero, or any other number
• Mixture of two distributions
  • “Point mass” distribution that determines the probability of zero
  • Continuous distribution that gives us the shape of the nonzero values

Count data that contains zeros: Zero Inflated Poisson (ZIP)

• Many datasets that count the abundance of organisms have lots of zero observations
• Poisson and negative binomial distributions are both discrete distributions that work well with count data
• But they may not fit datasets well that have a lot of zeros
• They already have some probability of zero, we just need to add an extra parameter to inflate that number
• Two kinds of zeros in the model: some from the main discrete distribution and others from the extra inflation component

Example count dataset with zeros