Bayesian Statistical Modelling with Stan and brms (BMSB01)

This session transitions from analytical Bayesian inference to numerical methods via MCMC. We begin by explaining why analytical approaches, while pedagogically valuable, are only possible in restricted cases with conjugate priors and simple models. Markov Chain Monte Carlo methods are introduced as a general numerical solution that can be applied to virtually any Bayesian model. Rather than calculating posterior distributions analytically, MCMC generates samples from the posterior distribution, which can then be used to approximate any posterior quantity. We introduce Stan as a state-of-the-art MCMC implementation and brms as a high-level R interface to Stan that allows fitting models with familiar R formula syntax. To demonstrate the connection with Day 1, we first re-analyze the Bernoulli model using brms, showing that MCMC recovers the same posterior we calculated analytically. We then fit our first Bayesian linear regression using brm and compare results with lm, examining similarities and differences. The session covers understanding brms output, basic visualization of MCMC samples, and interpreting posterior summaries.

Break – 01:00:00

Session 5 – 02:00:00 – Bayesian Linear Regression

This session provides in-depth coverage of Bayesian linear regression. We work through multiple regression models with continuous and categorical predictors, examining the posterior distribution over regression coefficients and comparing Bayesian credible intervals with frequentist confidence intervals. MCMC diagnostics are covered in detail: trace plots showing the sampling path of each chain, effective sample size measuring how many independent samples we have, and Rhat statistics indicating convergence across chains. Participants learn to recognize when MCMC has not converged and what to do about it. We explore varying intercept and varying slope formulations with categorical predictors, understanding how these relate to interaction terms. The brms functions for visualization are introduced: mcmc_plot for posterior distributions, pp_check for posterior predictive checks. The stancode function is used to examine the underlying Stan model, helping participants understand what brms is doing behind the scenes. Through multiple worked examples with real datasets, participants develop fluency in fitting and interpreting Bayesian regression models.

This session addresses two critical aspects of Bayesian analysis: prior specification and model comparison. We begin by examining default priors in brms using prior_summary, understanding that these are weakly informative priors designed to regularize without dominating the likelihood. We then cover how to specify custom priors using set_prior, including different priors for different parameters (intercepts, slopes, variance parameters). Prior predictive checks are introduced: simulating data from the prior before seeing real data to ensure priors are sensible. Prior sensitivity analysis demonstrates how to assess whether posterior inference is robust to prior specification. The session then turns to model comparison, covering multiple approaches: WAIC (Widely Applicable Information Criterion) and LOO (Leave-One-Out cross-validation) for comparing predictive accuracy, and Bayes factors for formal hypothesis comparison. We discuss when each method is appropriate and how to interpret differences between models. Practical workflows for comparing multiple models are demonstrated, including comparing models that differ in fixed effects, interaction terms, or distributional assumptions.

This session demonstrates how Bayesian methods easily extend beyond the restrictive assumptions of normal linear models. Classical normal linear models assume normally distributed residuals with constant variance, but these assumptions often fail with real data. We show how to relax these assumptions in brms. First, robust regression using Student-t distributions is covered: replacing the normal distribution for residuals with a t-distribution allows heavier tails and reduces sensitivity to outliers. The degrees of freedom parameter can be estimated from the data. Second, distributional regression is introduced using the bf (brmsformula) syntax to model sigma (residual standard deviation) as a function of predictors, handling heteroskedasticity directly. We examine real datasets where these extensions improve model fit, using pp_check to visualize how well models capture the data-generating process. Comparing normal, t-distributed, and heteroskedastic models using WAIC demonstrates the practical value of these extensions. This session illustrates a key advantage of Bayesian methods: the ease of modifying model assumptions to match data properties.

This session extends the Bayesian framework to generalized linear models for non-normal response variables. We begin with binary logistic regression, the most widely applicable GLM. The logit link function is explained, connecting linear predictors to probabilities. We fit logistic regression models using brm with family = bernoulli(), interpreting coefficients on the log-odds scale and transforming to probabilities for interpretation. Model comparison for logistic models is demonstrated using LOO. We then briefly cover other GLM families to show the breadth of models available in brms: Poisson regression for count data, negative binomial regression for overdispersed counts, and ordinal regression for ordered categorical outcomes. Throughout, we emphasize that the Bayesian workflow remains the same regardless of model family: specify the model, examine MCMC diagnostics, check priors, evaluate fit via posterior predictive checks, and compare models using information criteria. Real datasets from diverse research domains illustrate these models in practice.

This session introduces Bayesian approaches to mixed effects (multilevel/hierarchical) models for grouped or correlated data. We begin by explaining when mixed models are needed: repeated measures, clustered data, or any situation where observations are grouped and share common features. The varying intercept model is introduced, allowing each group to have its own baseline level while estimating an overall population mean and between-group variance. We extend this to varying slope models, where the effect of a predictor can vary across groups. The brms syntax for mixed models (using the lme4-style formula with random effects specified as (1|group) for varying intercepts and (predictor|group) for varying slopes) is covered in detail. We compare brms output with lme4/lmer output to understand similarities and differences. A key advantage of Bayesian mixed models is demonstrated: they often converge where lme4 struggles, particularly with complex random effects structures, small sample sizes, or binary/count data. Multiple real datasets with hierarchical structure are analyzed, including cases with both linear and generalized linear mixed models. The session concludes with practical guidance on specifying, fitting, and interpreting Bayesian mixed models, and pointers to further learning for more complex hierarchical structures.