# Applied Bayesian Modelling: Regression, Hierarchical Models, and Clinical Applications {#sec-applied-bayesian}
```{r}
#| include: false
library(tidyverse)
library(knitr)
```
## Introduction
The previous chapter introduced the conceptual foundations of Bayesian inference: priors, likelihoods, posteriors, and MCMC. This chapter puts those ideas to work. We will fit Bayesian regression models, learn the critical model-checking steps that distinguish careful Bayesian analysis from careless button-pressing, and then tackle the most powerful application of the Bayesian framework in clinical research: **hierarchical (multilevel) models**.
Hierarchical models deserve special attention because they solve a problem that arises constantly in medicine. Clinical data are almost always grouped: patients within hospitals, measurements within patients, sites within multi-centre trials. Frequentist mixed-effects models handle this, but the Bayesian hierarchical framework does it more flexibly and with clearer interpretation. By the end of this chapter, you will understand partial pooling, shrinkage, and why these ideas matter for clinical decision-making.
## Bayesian Linear Regression
### How It Differs from Frequentist Regression
In frequentist ordinary least squares (OLS), we obtain point estimates of regression coefficients and standard errors. A 95% confidence interval is a statement about long-run coverage.
In Bayesian linear regression, each coefficient has a **full posterior distribution**. We obtain not just an estimate and an interval, but the entire shape of our uncertainty. This is the same linear model:
$$
y_i = \beta_0 + \beta_1 x_{1i} + \cdots + \beta_p x_{pi} + \epsilon_i, \quad \epsilon_i \sim \text{Normal}(0, \sigma^2)
$$
Reading this equation in words: the outcome for patient $i$ (say, their blood pressure, $y_i$) is built up from a baseline value plus a contribution from each predictor, plus some leftover noise.
- $y_i$ is the outcome we want to model for patient $i$.
- $\beta_0$ (the **intercept**) is the predicted outcome when every predictor equals zero --- a baseline level.
- $\beta_1, \ldots, \beta_p$ (the **slopes** or coefficients) say how much the outcome changes for a one-unit increase in each predictor $x_{1}, \ldots, x_p$. For example, $\beta_1 = 0.5$ for age means the outcome rises by 0.5 units for each additional year.
- $x_{1i}, \ldots, x_{pi}$ are the predictor values measured on patient $i$ (age, BMI, and so on).
- $\epsilon_i$ is the **residual** --- the part of patient $i$'s outcome that the predictors do not explain. We assume it is random noise centred on zero with standard deviation $\sigma$, which captures how much patients scatter around the line.
The symbol $\sim \text{Normal}(0, \sigma^2)$ simply means "is drawn from a bell-shaped (normal) distribution centred at 0 with spread $\sigma$." This is the same linear model you would fit with `lm()`; nothing about the structure changes.
What *is* new in the Bayesian version is that we place priors on all parameters: $\beta_0, \beta_1, \ldots, \beta_p, \sigma$. A prior is just a starting belief about plausible values before seeing the data. The MCMC sampler then combines those priors with the data and generates draws from the joint posterior $P(\beta_0, \beta_1, \ldots, \beta_p, \sigma \mid \mathbf{y})$ --- our updated belief about every parameter once the data have been taken into account.
### Clinical Example: Predicting Systolic Blood Pressure
We will model systolic blood pressure (SBP) as a function of age, BMI, and whether the patient is on antihypertensive medication.
::: {.callout-important}
## Load your packages first
Every R example in this chapter (and the exercises) assumes you have loaded the **tidyverse** --- this is where functions like `tibble()`, `mutate()`, the pipe `%>%`, and `ggplot()` come from. If you skip it, R will not recognise these functions. We show `library(tidyverse)` at the top of each code block as a reminder; you only need to run it once per session. The Bayesian models additionally need `library(brms)`, and some plots need `library(bayesplot)`.
:::
::: {.panel-tabset}
### R
```{r}
#| label: bayesian-lm-r
#| message: false
#| warning: false
library(tidyverse) # provides tibble(), mutate(), %>%, ggplot()
library(brms) # provides brm() for Bayesian models
# Simulate clinical data
set.seed(123)
n <- 200
bp_data <- tibble(
age = round(rnorm(n, 55, 12)),
bmi = round(rnorm(n, 28, 5), 1),
on_meds = rbinom(n, 1, 0.4),
sbp = round(90 + 0.5 * age + 0.8 * bmi - 8 * on_meds + rnorm(n, 0, 12))
)
# Fit Bayesian linear regression with weakly informative priors
fit_bp <- brm(
sbp ~ age + bmi + on_meds,
data = bp_data,
family = gaussian(),
prior = c(
prior(normal(0, 20), class = "b"), # weakly informative for slopes
prior(normal(120, 30), class = "Intercept"), # centred on typical SBP
prior(exponential(0.1), class = "sigma") # positive scale parameter
),
chains = 4,
iter = 2000,
warmup = 1000,
seed = 42,
silent = 2,
refresh = 0
)
summary(fit_bp)
```
```{r}
#| label: fig-posterior-coefs-r
#| fig-cap: "Posterior distributions of regression coefficients for the blood pressure model."
#| fig-width: 9
#| fig-height: 5
library(tidyverse)
library(bayesplot) # provides mcmc_areas() for posterior plots
mcmc_areas(
fit_bp,
pars = c("b_age", "b_bmi", "b_on_meds"),
prob = 0.95,
prob_outer = 0.99
) +
labs(
title = "Posterior Distributions of Regression Coefficients",
subtitle = "Shaded region = 95% credible interval"
) +
theme_minimal(base_size = 13)
```
### Python
```{python}
#| label: bayesian-lm-py
#| message: false
import numpy as np
import pandas as pd
import pymc as pm
import arviz as az
np.random.seed(123)
n = 200
bp_data = pd.DataFrame({
'age': np.round(np.random.normal(55, 12, n)),
'bmi': np.round(np.random.normal(28, 5, n), 1),
'on_meds': np.random.binomial(1, 0.4, n)
})
bp_data['sbp'] = np.round(
90 + 0.5 * bp_data['age'] + 0.8 * bp_data['bmi']
- 8 * bp_data['on_meds'] + np.random.normal(0, 12, n)
)
with pm.Model() as bp_model:
# Priors
intercept = pm.Normal("Intercept", mu=120, sigma=30)
b_age = pm.Normal("b_age", mu=0, sigma=20)
b_bmi = pm.Normal("b_bmi", mu=0, sigma=20)
b_meds = pm.Normal("b_on_meds", mu=0, sigma=20)
sigma = pm.Exponential("sigma", lam=0.1)
# Linear model
mu = intercept + b_age * bp_data['age'].values + \
b_bmi * bp_data['bmi'].values + \
b_meds * bp_data['on_meds'].values
# Likelihood
y_obs = pm.Normal("sbp", mu=mu, sigma=sigma, observed=bp_data['sbp'].values)
# Sample
trace_bp = pm.sample(1000, tune=1000, chains=4, random_seed=42,
progressbar=False)
print(az.summary(trace_bp, var_names=["Intercept", "b_age", "b_bmi",
"b_on_meds", "sigma"]))
```
```{python}
#| label: fig-posterior-coefs-py
#| fig-cap: "Posterior distributions of regression coefficients for the blood pressure model."
az.plot_forest(trace_bp, var_names=["b_age", "b_bmi", "b_on_meds"],
combined=True, hdi_prob=0.95, figsize=(8, 4))
import matplotlib.pyplot as plt
plt.title("Posterior Distributions of Regression Coefficients\n95% HDI")
plt.tight_layout()
plt.show()
```
:::
The output shows the **posterior mean**, **standard deviation** (the Bayesian analogue of the standard error), and the **95% credible interval** for each coefficient. Notice that the on-medication effect has a credible interval that is entirely below zero, providing direct evidence that medication lowers SBP.
::: {.callout-note}
## What are `chains`, `iter`, and `warmup`?
These three arguments control the MCMC sampler and appear in every `brm()` call in this chapter, so it is worth pausing on them.
- **`chains = 4`** runs the sampler four separate times from different starting points. If all four chains end up exploring the same region, that is reassuring evidence the sampler has found the right answer (this is what the $\hat{R}$ diagnostic checks).
- **`iter = 2000`** is the total number of steps each chain takes.
- **`warmup = 1000`** (sometimes called *burn-in*) is the number of *initial* steps that are **thrown away**. When a chain starts, it begins at an arbitrary point that may be far from the posterior, and it also needs time to tune its internal step size. Those early, unreliable samples would distort our estimates, so we discard them and keep only what comes after.
So `iter = 2000` with `warmup = 1000` means each chain takes 2000 steps but we *keep* only the last 1000 --- giving $4 \times 1000 = 4000$ usable posterior draws. A useful analogy: warm-up is like letting a car engine idle until it reaches operating temperature before you start the timed lap. The warm-up laps still happen; you simply do not record them.
:::
## Bayesian Logistic Regression
For binary outcomes (e.g., 30-day readmission, treatment response), we use Bayesian logistic regression. The model is identical to the frequentist version except that all parameters receive priors:
$$
\text{logit}(P(y_i = 1)) = \beta_0 + \beta_1 x_{1i} + \cdots + \beta_p x_{pi}
$$
::: {.panel-tabset}
### R
```{r}
#| label: bayesian-logistic-r
#| message: false
#| warning: false
library(tidyverse)
library(brms)
# Simulate readmission data
set.seed(456)
n <- 300
readmit_data <- tibble(
age = round(rnorm(n, 65, 10)),
charlson = rpois(n, 2),
los = rpois(n, 5) + 1, # length of stay in days
readmit_30 = rbinom(
n,
1,
plogis(
-3 +
0.02 * round(rnorm(n, 65, 10)) +
0.3 * rpois(n, 2) +
0.1 * (rpois(n, 5) + 1)
)
)
)
fit_readmit <- brm(
readmit_30 ~ age + charlson + los,
data = readmit_data,
family = bernoulli(link = "logit"),
prior = c(
prior(normal(0, 2.5), class = "b"), # weakly informative on log-odds
prior(normal(0, 5), class = "Intercept")
),
chains = 4,
iter = 2000,
warmup = 1000,
seed = 42,
silent = 2,
refresh = 0
)
# Posterior summary on the odds ratio scale
posterior_samples <- as_draws_df(fit_readmit)
or_summary <- posterior_samples %>%
transmute(
OR_age = exp(b_age),
OR_charlson = exp(b_charlson),
OR_los = exp(b_los)
) %>%
summarise(across(
everything(),
list(mean = mean, q2.5 = ~ quantile(., 0.025), q97.5 = ~ quantile(., 0.975))
))
cat("Posterior odds ratios (mean [95% CrI]):\n")
cat(sprintf(
" Age: %.3f [%.3f, %.3f]\n",
or_summary$OR_age_mean,
or_summary$OR_age_q2.5,
or_summary$OR_age_q97.5
))
cat(sprintf(
" Charlson: %.3f [%.3f, %.3f]\n",
or_summary$OR_charlson_mean,
or_summary$OR_charlson_q2.5,
or_summary$OR_charlson_q97.5
))
cat(sprintf(
" LOS: %.3f [%.3f, %.3f]\n",
or_summary$OR_los_mean,
or_summary$OR_los_q2.5,
or_summary$OR_los_q97.5
))
```
### Python
```{python}
#| label: bayesian-logistic-py
#| message: false
import numpy as np
import pandas as pd
import pymc as pm
import arviz as az
np.random.seed(456)
n = 300
readmit_data = pd.DataFrame({
'age': np.round(np.random.normal(65, 10, n)),
'charlson': np.random.poisson(2, n),
'los': np.random.poisson(5, n) + 1
})
from scipy.special import expit
lp = -3 + 0.02 * readmit_data['age'] + 0.3 * readmit_data['charlson'] + \
0.1 * readmit_data['los']
readmit_data['readmit_30'] = np.random.binomial(1, expit(lp))
with pm.Model() as readmit_model:
intercept = pm.Normal("Intercept", mu=0, sigma=5)
b_age = pm.Normal("b_age", mu=0, sigma=2.5)
b_charl = pm.Normal("b_charlson", mu=0, sigma=2.5)
b_los = pm.Normal("b_los", mu=0, sigma=2.5)
logit_p = intercept + b_age * readmit_data['age'].values + \
b_charl * readmit_data['charlson'].values + \
b_los * readmit_data['los'].values
y_obs = pm.Bernoulli("readmit", logit_p=logit_p,
observed=readmit_data['readmit_30'].values)
trace_readmit = pm.sample(1000, tune=1000, chains=4, random_seed=42,
progressbar=False)
# Compute odds ratios from posterior
posterior = az.extract(trace_readmit)
for var in ["b_age", "b_charlson", "b_los"]:
or_vals = np.exp(posterior[var].values)
print(f"OR {var}: {or_vals.mean():.3f} "
f"[{np.percentile(or_vals, 2.5):.3f}, "
f"{np.percentile(or_vals, 97.5):.3f}]")
```
:::
## Prior Predictive Checks
A prior predictive check asks: **before seeing any data, does the model generate plausible outcomes?** This is a critical step that is often skipped. If your priors allow the model to predict systolic blood pressures of 500 mmHg or negative values, your priors are too vague.
The procedure:
1. Sample parameter values from the priors (not the posterior).
2. Simulate data from the model using those parameter values.
3. Plot the simulated data and check whether it covers a plausible range.
::: {.panel-tabset}
### R
```{r}
#| label: fig-prior-predictive-r
#| fig-cap: "Prior predictive check: simulated SBP distributions under the chosen priors."
#| fig-width: 6.5
#| fig-height: 4
#| message: false
#| warning: false
library(tidyverse)
library(brms)
# Generate prior predictive samples
pp_check_prior <- brm(
sbp ~ age + bmi + on_meds,
data = bp_data,
family = gaussian(),
prior = c(
prior(normal(0, 20), class = "b"),
prior(normal(120, 30), class = "Intercept"),
prior(exponential(0.1), class = "sigma")
),
sample_prior = "only",
chains = 4,
iter = 2000,
warmup = 1000,
seed = 42,
silent = 2,
refresh = 0
)
pp_samples <- posterior_predict(pp_check_prior)
# Plot density of prior predictive samples vs observed data
tibble(
simulated = as.vector(pp_samples[sample(nrow(pp_samples), 50), ]),
type = "Prior predictive"
) %>%
ggplot(aes(x = simulated)) +
geom_density(fill = "steelblue", alpha = 0.3) +
geom_density(data = bp_data, aes(x = sbp), fill = "firebrick", alpha = 0.3) +
labs(
x = "Systolic Blood Pressure (mmHg)",
y = "Density",
title = "Prior Predictive Check",
subtitle = "Blue = simulated from priors; Red = observed data"
) +
xlim(-100, 350) +
theme_minimal(base_size = 13)
```
### Python
```{python}
#| label: fig-prior-predictive-py
#| fig-cap: "Prior predictive check: simulated SBP distributions under the chosen priors."
#| fig-width: 6.5
#| fig-height: 4
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
n_sim = 2000
# Simulate from priors
intercepts = np.random.normal(120, 30, n_sim)
b_ages = np.random.normal(0, 20, n_sim)
b_bmis = np.random.normal(0, 20, n_sim)
b_meds_arr = np.random.normal(0, 20, n_sim)
sigmas = np.random.exponential(10, n_sim)
# Pick a "typical" patient: age=55, bmi=28, on_meds=0
sim_sbp = intercepts + b_ages * 55 + b_bmis * 28 + b_meds_arr * 0
sim_sbp += np.random.normal(0, sigmas)
fig, ax = plt.subplots(figsize=(8, 5))
ax.hist(sim_sbp, bins=80, density=True, alpha=0.4, color="steelblue",
label="Prior predictive", range=(-500, 800))
ax.hist(bp_data['sbp'].values, bins=30, density=True, alpha=0.4,
color="firebrick", label="Observed data")
ax.set_xlabel("Systolic Blood Pressure (mmHg)")
ax.set_ylabel("Density")
ax.set_title("Prior Predictive Check\nBlue = simulated from priors; Red = observed")
ax.set_xlim(-200, 400)
ax.legend()
plt.tight_layout()
plt.show()
```
:::
**How to read this plot.** The blue curve is the range of blood pressures the model considers possible *based on the priors alone*, before it has seen a single patient. The red curve is the real data, shown only for comparison. What we are looking for is a sanity check: the blue curve should comfortably cover the clinically realistic range (roughly 80--200 mmHg) without putting serious weight on impossible values such as negative pressures or readings above 400 mmHg. Here the blue curve is wide --- it does not yet know the typical SBP --- but it is centred in a sensible place and is not absurd, so these priors are acceptable. It is meant to be broader than the data; the data will sharpen it later.
If the prior predictive distribution is wildly broader than any plausible clinical range, consider tightening the priors. If it is too narrow and excludes plausible values, widen them. This iterative process happens **before** fitting the model to data.
## Posterior Predictive Checks
After fitting, we ask: **does the fitted model generate data that look like the observed data?** This is the Bayesian equivalent of residual diagnostics.
The procedure:
1. Draw parameter values from the posterior.
2. Simulate new data from the model using those posterior draws.
3. Compare the distribution of simulated data to the observed data.
If the model is well specified, the simulated data should be nearly indistinguishable from the real data in terms of summary statistics (mean, variance, range, shape).
**How to read this plot.** Each light blue line is one fake dataset the *fitted* model generates; the dark line is the observed data. Unlike the prior predictive check, here the model has already learned from the data, so we expect a close match. If the dark line sits comfortably within the spread of blue lines --- same centre, same width, same overall shape --- the model is reproducing the data well. If the dark line falls outside the blue cloud (for example, the data are skewed but the model only produces symmetric curves), the model is misspecified and the likelihood or priors need revisiting.
::: {.panel-tabset}
### R
```{r}
#| label: fig-posterior-predictive-r
#| fig-cap: "Posterior predictive check: the model-generated data (blue) should resemble the observed data (dark line)."
#| fig-width: 6.5
#| fig-height: 4
#| message: false
library(tidyverse)
library(brms)
pp_check(fit_bp, type = "dens_overlay", ndraws = 100) +
labs(
title = "Posterior Predictive Check",
subtitle = "Light blue = simulated from posterior; Dark = observed"
) +
theme_minimal(base_size = 13)
```
### Python
```{python}
#| label: fig-posterior-predictive-py
#| fig-cap: "Posterior predictive check for the blood pressure model."
#| fig-width: 6.5
#| fig-height: 4
import pymc as pm
import arviz as az
import matplotlib.pyplot as plt
# Using the previously fitted model
with bp_model:
ppc = pm.sample_posterior_predictive(trace_bp, random_seed=42,
progressbar=False)
az.plot_ppc(az.from_pymc3(posterior_predictive=ppc,
observed_data={"sbp": bp_data['sbp'].values}),
num_pp_samples=100, figsize=(8, 5))
plt.title("Posterior Predictive Check")
plt.tight_layout()
plt.show()
```
:::
## Bayesian Hierarchical (Multilevel) Models
### Why They Matter
Consider a multi-site clinical trial testing a new antihypertensive drug across 12 hospitals. The treatment effect might vary from site to site due to differences in patient populations, clinical protocols, or local practice patterns. How should we estimate the treatment effect at each site?
There are three naive approaches:
1. **Complete pooling**: ignore the site structure and estimate a single treatment effect. This assumes all sites are identical --- clearly wrong.
2. **No pooling**: estimate a separate treatment effect at each site. This ignores that the sites are studying the same drug. Small sites will have wildly imprecise estimates.
3. **Partial pooling (hierarchical model)**: the Bayesian approach. Each site has its own treatment effect, but those effects are drawn from a shared distribution. Sites with small samples are **pulled toward the overall mean**, while sites with large samples retain their individual estimates.
This pulling toward the grand mean is called **shrinkage**, and it is one of the most powerful ideas in applied statistics.
### The Model Structure
A Bayesian hierarchical model for treatment effects across $J$ hospitals:
$$
y_{ij} = \beta_0 + \beta_1 \text{treatment}_{ij} + u_j + \epsilon_{ij}
$$
Reading the equation: $y_{ij}$ is the outcome for patient $i$ in hospital $j$. The first two terms, $\beta_0 + \beta_1 \text{treatment}_{ij}$, are the *average* picture across all hospitals --- the same intercept and treatment effect we would get from an ordinary regression. The new piece is $u_j$, a small "nudge" that is specific to hospital $j$. It captures how that hospital departs from the average. The terms are:
- $\beta_0$ is the overall baseline outcome (the **fixed intercept**, shared by all hospitals).
- $\beta_1$ is the **average treatment effect** across all hospitals.
- $u_j$ is the **random effect** for hospital $j$: how much hospital $j$'s own treatment effect sits above or below the average $\beta_1$. A hospital with $u_j = +2$ responds 2 units more than the typical site; one with $u_j = -2$ responds 2 units less.
- $\epsilon_{ij}$ is the residual error for an individual patient --- the usual unexplained patient-to-patient noise.
**What does "random effect" mean?** Instead of estimating each hospital's deviation $u_j$ in complete isolation, we assume the deviations themselves come from a shared distribution:
$$
u_j \sim \text{Normal}(0, \tau^2).
$$
This single line is what makes the model "hierarchical." It says the hospital effects scatter around zero (zero meaning "exactly average") with a spread of $\tau$. Because all the $u_j$ are tied together by this common distribution, information is shared across hospitals: a hospital with very few patients borrows strength from the others rather than relying only on its own noisy data. That borrowing is exactly the partial pooling and shrinkage we explore below.
**The role of $\tau$.** The parameter $\tau$ is the **between-hospital standard deviation** --- it answers the clinically important question "how much do hospitals actually differ?"
- If $\tau$ is close to **zero**, the hospitals are nearly identical, every $u_j$ is tiny, and the model collapses toward a single shared treatment effect (complete pooling).
- If $\tau$ is **large**, hospitals differ substantially, and each site keeps much of its own estimate (closer to no pooling).
Crucially, $\tau$ is not fixed in advance; it is **estimated from the data**, so the model learns how much pooling is warranted rather than assuming it.
In the Bayesian framework, $\tau$ also gets a prior. Because $\tau$ is a standard deviation it cannot be negative, so we use a prior defined only on positive values. Two common choices are the **half-normal** and the **half-Cauchy** --- the right-hand halves of the ordinary normal and Cauchy distributions. Both place most of their weight near zero (gently regularising toward "hospitals are similar") while still permitting larger values. The difference is in the tail: the half-Cauchy has a much heavier tail, so it is more willing to allow a large between-hospital spread if the data demand it. This makes the half-Cauchy a popular, mildly conservative default for group-level standard deviations, especially when the number of groups is small.
::: {.panel-tabset}
### R
```{r}
#| label: fig-tau-priors-r
#| fig-cap: "Two common priors for the between-hospital standard deviation tau. Both concentrate near zero, but the half-Cauchy has a heavier tail, allowing larger between-group variation."
#| fig-width: 6.5
#| fig-height: 4
#| message: false
#| warning: false
library(tidyverse)
tau_grid <- seq(0, 20, length.out = 400)
prior_df <- bind_rows(
tibble(
tau = tau_grid,
density = 2 * dnorm(tau_grid, 0, 5), # half-normal(0, 5)
prior = "Half-Normal(0, 5)"
),
tibble(
tau = tau_grid,
density = 2 * dcauchy(tau_grid, 0, 5), # half-Cauchy(0, 5)
prior = "Half-Cauchy(0, 5)"
)
)
ggplot(prior_df, aes(x = tau, y = density, colour = prior)) +
geom_line(linewidth = 1) +
scale_colour_manual(
values = c(
"Half-Normal(0, 5)" = "steelblue",
"Half-Cauchy(0, 5)" = "firebrick"
)
) +
labs(
x = expression(tau ~ "(between-hospital SD)"),
y = "Prior density",
colour = "Prior",
title = "Priors for a Group-Level Standard Deviation",
subtitle = "Both favour small tau; the half-Cauchy keeps a heavier tail"
) +
theme_minimal(base_size = 13) +
theme(legend.position = "bottom")
```
### Python
```{python}
#| label: fig-tau-priors-py
#| fig-cap: "Two common priors for the between-hospital standard deviation tau."
#| fig-width: 6.5
#| fig-height: 4
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
tau_grid = np.linspace(0, 20, 400)
half_normal = 2 * stats.norm.pdf(tau_grid, 0, 5) # half-normal(0, 5)
half_cauchy = 2 * stats.cauchy.pdf(tau_grid, 0, 5) # half-Cauchy(0, 5)
fig, ax = plt.subplots(figsize=(6.5, 4))
ax.plot(tau_grid, half_normal, color="steelblue", linewidth=2,
label="Half-Normal(0, 5)")
ax.plot(tau_grid, half_cauchy, color="firebrick", linewidth=2,
label="Half-Cauchy(0, 5)")
ax.set_xlabel(r"$\tau$ (between-hospital SD)")
ax.set_ylabel("Prior density")
ax.set_title("Priors for a Group-Level Standard Deviation\n"
"Both favour small tau; the half-Cauchy keeps a heavier tail")
ax.legend()
plt.tight_layout()
plt.show()
```
:::
### Clinical Example: Treatment Effects Across Hospitals
::: {.panel-tabset}
### R
```{r}
#| label: fig-hierarchical-r
#| fig-cap: "Shrinkage in a Bayesian hierarchical model: hospital-specific estimates are pulled toward the grand mean."
#| fig-width: 10
#| fig-height: 6
#| message: false
#| warning: false
library(tidyverse) # tibble(), map2_dfr(), %>%, ggplot()
library(brms)
# Simulate multi-site trial data
set.seed(789)
n_hospitals <- 12
patients_per_hospital <- c(15, 20, 25, 30, 35, 40, 50, 60, 80, 100, 120, 150)
true_grand_effect <- -8 # mmHg reduction in SBP
tau_true <- 3 # between-hospital SD
hospital_effects <- rnorm(n_hospitals, 0, tau_true)
trial_data <- map2_dfr(1:n_hospitals, patients_per_hospital, function(j, nj) {
tibble(
hospital = factor(j),
treatment = rep(0:1, length.out = nj),
sbp = 140 +
(true_grand_effect + hospital_effects[j]) * treatment +
rnorm(nj, 0, 15)
)
})
# Fit hierarchical model
fit_hier <- brm(
sbp ~ treatment + (treatment | hospital),
data = trial_data,
prior = c(
prior(normal(0, 20), class = "b"),
prior(normal(140, 30), class = "Intercept"),
prior(cauchy(0, 5), class = "sd"),
prior(exponential(0.1), class = "sigma")
),
chains = 4,
iter = 2000,
warmup = 1000,
seed = 42,
silent = 2,
refresh = 0,
control = list(adapt_delta = 0.95)
)
# Extract hospital-specific treatment effects (partial pooling)
ranefs <- ranef(fit_hier)$hospital[, , "treatment"]
grand_mean <- fixef(fit_hier)["treatment", "Estimate"]
# No-pooling estimates (separate OLS per hospital)
no_pool <- trial_data %>%
group_by(hospital) %>%
summarise(
no_pool_est = coef(lm(sbp ~ treatment))[2],
n = n(),
.groups = "drop"
) %>%
mutate(
hospital_num = as.numeric(hospital),
partial_pool_est = grand_mean + ranefs[, "Estimate"],
# Label each hospital with its sample size so the reader can see
# which sites are small (and therefore shrink most)
hospital_label = paste0("Hospital ", hospital, " (n = ", n, ")")
)
# Plot shrinkage
ggplot(no_pool, aes(y = reorder(hospital_label, n))) +
geom_point(aes(x = no_pool_est, colour = "No pooling"), size = 3) +
geom_point(aes(x = partial_pool_est, colour = "Partial pooling"), size = 3) +
geom_vline(xintercept = grand_mean, linetype = "dashed", colour = "grey40") +
annotate(
"text",
x = grand_mean + 0.5,
y = 12.5,
label = paste0("Grand mean = ", round(grand_mean, 1)),
hjust = 0,
size = 3.5
) +
geom_segment(
aes(
x = no_pool_est,
xend = partial_pool_est,
yend = reorder(hospital_label, n)
),
arrow = arrow(length = unit(0.15, "cm")),
colour = "grey60"
) +
scale_colour_manual(
values = c("No pooling" = "steelblue", "Partial pooling" = "firebrick")
) +
labs(
x = "Treatment Effect (mmHg change in SBP)",
y = "Hospital (ordered by sample size)",
colour = "Estimate Type",
title = "Shrinkage in a Hierarchical Model",
subtitle = "Arrows show how partial pooling pulls estimates toward the grand mean"
) +
theme_minimal(base_size = 13) +
theme(legend.position = "bottom")
```
### Python
```{python}
#| label: fig-hierarchical-py
#| fig-cap: "Shrinkage in a Bayesian hierarchical model."
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
np.random.seed(789)
n_hospitals = 12
patients = [15, 20, 25, 30, 35, 40, 50, 60, 80, 100, 120, 150]
grand_effect = -8
tau = 3
hospital_effects = np.random.normal(0, tau, n_hospitals)
rows = []
for j in range(n_hospitals):
nj = patients[j]
trt = np.tile([0, 1], nj // 2 + 1)[:nj]
sbp = 140 + (grand_effect + hospital_effects[j]) * trt + \
np.random.normal(0, 15, nj)
for i in range(nj):
rows.append({'hospital': j, 'treatment': trt[i], 'sbp': sbp[i]})
trial_data = pd.DataFrame(rows)
# No-pooling estimates
no_pool = []
for j in range(n_hospitals):
df_j = trial_data[trial_data['hospital'] == j]
trt_mean = df_j[df_j['treatment'] == 1]['sbp'].mean()
ctl_mean = df_j[df_j['treatment'] == 0]['sbp'].mean()
no_pool.append(trt_mean - ctl_mean)
no_pool = np.array(no_pool)
grand_mean_est = np.mean(no_pool)
# Simulate partial pooling (shrinkage toward grand mean)
# Shrinkage factor depends on sample size: larger sites shrink less
shrinkage_factor = np.array([15 / (15 + nj) for nj in patients])
partial_pool = grand_mean_est + (1 - shrinkage_factor) * (no_pool - grand_mean_est)
# Plot
fig, ax = plt.subplots(figsize=(10, 6))
order = np.argsort(patients)
y_pos = np.arange(n_hospitals)
for i, idx in enumerate(order):
ax.plot([no_pool[idx], partial_pool[idx]], [i, i],
color="grey", linewidth=1, zorder=1)
ax.annotate("", xy=(partial_pool[idx], i), xytext=(no_pool[idx], i),
arrowprops=dict(arrowstyle="->", color="grey"))
ax.scatter(no_pool[order], y_pos, color="steelblue", s=60, zorder=2,
label="No pooling")
ax.scatter(partial_pool[order], y_pos, color="firebrick", s=60, zorder=2,
label="Partial pooling")
ax.axvline(grand_mean_est, linestyle="--", color="grey", linewidth=1)
ax.set_yticks(y_pos)
ax.set_yticklabels([f"Hospital {order[i]+1}\n(n={patients[order[i]]})"
for i in range(n_hospitals)], fontsize=9)
ax.set_xlabel("Treatment Effect (mmHg)")
ax.set_title("Shrinkage in a Hierarchical Model\n"
"Arrows show estimates pulled toward the grand mean")
ax.legend()
plt.tight_layout()
plt.show()
```
:::
### Understanding Shrinkage
The figure above illustrates the key behaviour of hierarchical models:
- **Small hospitals** (few patients) are pulled strongly toward the grand mean. Their individual data are noisy, so the model relies more on the shared information from other sites.
- **Large hospitals** (many patients) retain estimates close to their no-pooling values. Their data are informative enough to override the group-level information.
- **The grand mean** itself is estimated from all sites, so every site contributes to it.
This is **partial pooling** --- a principled compromise between treating all sites as identical (complete pooling) and treating them as completely independent (no pooling). It is particularly valuable when some sites have very few patients, where no-pooling estimates would be unreliable.
::: {.callout-note}
## Shrinkage is not bias
Shrinkage toward the grand mean might seem like it introduces bias. In a narrow sense it does: individual site estimates are biased toward the mean. But this bias is traded for reduced variance, and the overall mean squared error is lower. This is the bias-variance trade-off that we discussed in earlier chapters, now appearing naturally in the Bayesian framework.
:::
### Varying Slopes
The model above includes a varying (random) slope for treatment by hospital. This means we are estimating not just "does the treatment effect vary across hospitals?" but also "by how much?" The between-hospital standard deviation $\tau$ directly quantifies this heterogeneity. If the 95% credible interval for $\tau$ excludes zero, there is evidence of meaningful variation across sites.
## The Practical Bayesian Workflow
A complete Bayesian analysis follows a structured workflow:
### Step 1: Specify the Model
- Write down the likelihood (data-generating process).
- Choose priors for all parameters.
- Justify your choices: are the priors weakly informative? Are they based on prior studies?
::: {.callout-tip}
## You don't have to write a likelihood from scratch
"Specify the likelihood" sounds like it requires deriving probability formulas by hand. In practice, with `brms` (R) or `bambi`/PyMC (Python) you almost never do. The likelihood is determined by the **type of outcome** you are modelling, and you select it with a single `family =` argument. The mapping is the same one you already use for choosing a regression in everyday practice:
| Your outcome | Example | `family` to use | Underlying likelihood |
|---|---|---|---|
| Continuous, roughly symmetric | Blood pressure, BMI | `gaussian()` | Normal |
| Binary (yes/no) | 30-day readmission, mortality | `bernoulli()` | Bernoulli (logistic) |
| Count | Number of admissions | `poisson()` | Poisson |
| Time-to-event | Survival time | `cox()` / `weibull()` | Survival |
So the practical question is not "what is the mathematical form of the likelihood?" but the familiar clinical question "**what kind of measurement is my outcome?**" Pick the matching family, and the software writes the likelihood for you. The equations in this chapter are there so you understand *what the software is doing* --- not because you need to type them out.
:::
### Step 2: Prior Predictive Check
- Simulate data from the priors alone.
- Verify that the implied range of outcomes is clinically plausible.
- Iterate if necessary: tighten or widen priors.
### Step 3: Fit the Model
- Run MCMC (typically 4 chains, 1000--2000 post-warmup iterations each).
- Check for warnings (divergent transitions, max treedepth warnings).
### Step 4: Diagnose Convergence
- Inspect trace plots: should look like "hairy caterpillars".
- Check $\hat{R} < 1.01$ for all parameters.
- Check ESS > 400 for both bulk and tail.
- If any diagnostic fails, do **not** interpret the results. Fix the model first.
### Step 5: Posterior Predictive Check
- Simulate data from the fitted posterior.
- Compare to observed data: distributions, summary statistics, patterns.
- If the model is badly misspecified, revise the likelihood or priors.
### Step 6: Summarise and Report
- Report posterior means (or medians), credible intervals, and relevant posterior probabilities.
- Include sensitivity analyses: how do results change under alternative priors?
- Report MCMC diagnostics alongside the results.
::: {.callout-tip}
## Reporting Template
"We estimated the treatment effect using a Bayesian [model type] with [prior specification]. The posterior mean treatment effect was X (95% credible interval: [L, U]). The posterior probability that the treatment effect exceeds the clinically meaningful threshold of Y was Z%. All chains converged ($\hat{R} < 1.01$, bulk ESS > [value], tail ESS > [value])."
:::
## Implementation: Software Choices
### R Ecosystem
**brms** (Bayesian Regression Models using Stan) is the recommended package for applied Bayesian analysis in R. It uses the familiar R formula syntax and compiles models in Stan behind the scenes. Essentially, anything you can fit with `lme4` can be fit with `brms`, but with priors and full posterior inference.
**rstanarm** provides pre-compiled Stan models for common regression types. It is faster to start (no compilation time) but less flexible than brms.
Both produce Stan code that you can inspect, which is excellent for learning.
### Python Ecosystem
**PyMC** (version 5+) is the most mature probabilistic programming library in Python. It uses the JAX or NumPy backends and provides NUTS sampling out of the box.
**Bambi** (BAyesian Model-Building Interface) is the Python analogue of brms: it provides a formula-based interface on top of PyMC, making it easy to specify complex models without writing raw PyMC code.
::: {.panel-tabset}
### R (brms)
```{r}
#| label: brms-workflow-example
#| eval: false
library(brms)
# Specify and fit
fit <- brm(
outcome ~ predictor1 + predictor2 + (1 | group),
data = my_data,
family = gaussian(),
prior = c(
prior(normal(0, 10), class = "b"),
prior(normal(0, 20), class = "Intercept"),
prior(cauchy(0, 5), class = "sd"),
prior(exponential(0.1), class = "sigma")
),
chains = 4,
iter = 2000,
warmup = 1000,
seed = 42
)
# Diagnose
plot(fit) # trace plots
summary(fit) # R-hat, ESS
pp_check(fit) # posterior predictive check
conditional_effects(fit) # visualise effects
# Posterior probability of meaningful effect
hypothesis(fit, "predictor1 > 0.5")
```
### Python (bambi)
```{python}
#| label: bambi-workflow-example
#| eval: false
import bambi as bmb
import arviz as az
# Specify and fit
model = bmb.Model("outcome ~ predictor1 + predictor2 + (1 | group)",
data=my_data, family="gaussian")
results = model.fit(draws=1000, tune=1000, chains=4, random_seed=42)
# Diagnose
az.plot_trace(results) # trace plots
print(az.summary(results)) # R-hat, ESS
az.plot_ppc(results) # posterior predictive check
# Posterior probability
posterior = az.extract(results)
prob = (posterior['predictor1'] > 0.5).mean()
print(f"P(predictor1 > 0.5 | data) = {prob:.3f}")
```
:::
## Exercises
::: {.callout-tip title="Exercise 1: Bayesian Logistic Regression for ICU Mortality"}
A dataset contains 500 ICU admissions with the following variables: age, APACHE II score, mechanical ventilation status (yes/no), and 28-day mortality (outcome).
a. Simulate a dataset with plausible parameter values.
b. Fit a Bayesian logistic regression model using weakly informative priors.
c. Perform a prior predictive check: do the priors imply plausible mortality rates?
d. Report posterior odds ratios with 95% credible intervals.
e. Calculate $P(\text{OR}_{\text{APACHE}} > 1.10 \mid \text{data})$ --- the probability that each unit increase in APACHE score increases the odds of death by more than 10%.
::: {.callout-note collapse="true" title="Solution"}
::: {.panel-tabset}
#### R
```{r}
#| eval: false
#| file: ../solutions/R/applied_bayesian_modeling_ex1.R
```
#### Python
```{python}
#| eval: false
#| file: ../solutions/python/applied_bayesian_modeling_ex1.py
```
:::
:::
:::
::: {.callout-tip title="Exercise 2: Hierarchical Model for Multi-Site Drug Trial"}
Twelve hospitals participate in a trial comparing a new statin to standard care. The primary outcome is change in LDL cholesterol (mg/dL) at 12 weeks. Hospitals vary in patient volume from 20 to 200 patients.
a. Simulate data where the true average treatment effect is --25 mg/dL with between-hospital SD of 5 mg/dL.
b. Fit a Bayesian hierarchical model with random intercepts and random slopes for treatment by hospital.
c. Create a shrinkage plot showing how hospital-specific estimates are pulled toward the grand mean.
d. Compare the hierarchical model estimates to separate hospital-by-hospital OLS regressions. Which approach is more appropriate and why?
::: {.callout-note collapse="true" title="Solution"}
::: {.panel-tabset}
#### R
```{r}
#| eval: false
#| file: ../solutions/R/applied_bayesian_modeling_ex2.R
```
#### Python
```{python}
#| eval: false
#| file: ../solutions/python/applied_bayesian_modeling_ex2.py
```
:::
:::
:::
::: {.callout-tip title="Exercise 3: Prior Sensitivity for Rare Events"}
A new surgical technique is tested in 40 patients. Zero patients experience a major adverse event.
a. Compute the posterior distribution for the adverse event rate using three priors: Beta(1,1), Beta(0.5, 0.5) (Jeffreys prior), and Beta(1, 9) (informative: prior belief ~10% rate).
b. Report the posterior mean and 95% upper credible bound for each prior.
c. A frequentist would report the point estimate as 0/40 = 0. Explain why the Bayesian estimates are more useful for regulatory decision-making.
::: {.callout-note collapse="true" title="Solution"}
::: {.panel-tabset}
#### R
```{r}
#| eval: false
#| file: ../solutions/R/applied_bayesian_modeling_ex3.R
```
#### Python
```{python}
#| eval: false
#| file: ../solutions/python/applied_bayesian_modeling_ex3.py
```
:::
:::
:::
## Summary
Bayesian regression models extend familiar linear and logistic regression by placing priors on parameters and yielding full posterior distributions rather than point estimates. The practical workflow --- specify, check priors, fit, diagnose, check posteriors, report --- ensures rigorous and transparent analysis. Hierarchical models, with their partial pooling and shrinkage properties, are especially powerful for multi-site clinical studies where borrowing information across groups improves estimation. Modern software (brms in R, PyMC/bambi in Python) makes these methods accessible to applied researchers.
::: {.callout-tip}
## Key Takeaways
1. Bayesian regression yields full posterior distributions for all coefficients, not just point estimates.
2. Prior and posterior predictive checks are essential quality-control steps.
3. Hierarchical models perform partial pooling: small groups shrink toward the grand mean, large groups retain their individual estimates.
4. Shrinkage reduces mean squared error by trading a small amount of bias for a large reduction in variance.
5. Always report MCMC diagnostics ($\hat{R}$, ESS, trace plots) alongside posterior summaries.
:::
## References and Further Reading
- For applied Bayesian modelling, see @gelman2007data, @gelman2013bda, and @mcelreath2020statistical.
- For Bayesian clinical trial design, see @fda2010bayesian and @berry2010bayesian.
- For brms software, see @buerkner2017brms.
- For hierarchical shrinkage and modern Bayesian workflows, see @efron1977stein and @martin2024bayesian (a practical guide to PyMC for applied researchers).
::: {.sectionrefs}
:::