Bio Applied Bayesian Statistics Python

Bayesian Statistics in Python

When to Use

• Need a full posterior distribution (not just a point estimate + p-value) for an effect size, dose-response slope, or group difference.
• Data has a grouped/nested structure (measurements within individuals, sites, species) and you want partial pooling instead of complete pooling or no pooling.
• Small sample sizes where informative priors from prior literature should be incorporated.
• Count data with overdispersion or excess zeros (Poisson doesn't fit; need Negative Binomial or zero-inflated models).
• Comparing competing models by predictive accuracy (LOO-CV, WAIC) rather than nested-model F-tests.

Version Compatibility

PyMC ≥ 5.10, ArviZ ≥ 0.17, Bambi ≥ 0.13, Python ≥ 3.10. PyMC 5.x uses PyTensor (not Theano/Aesara) as its backend — code written for PyMC3 (import pymc3, Theano ops) is not compatible.

Prerequisites

pip install pymc arviz bambi statsmodels palmerpenguins

Familiarity with linear regression and basic probability (prior/likelihood/posterior). For R users: install.packages(c("brms","bayesplot")).

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf
import statsmodels.api as sm
import pymc as pm
import arviz as az
import bambi as bmb

try:
    from palmerpenguins import load_penguins
    penguins = load_penguins().dropna()
except ImportError:
    rng = np.random.default_rng(42)
    n = 333
    flipper = rng.normal(200, 14, n)
    penguins = pd.DataFrame({
        "species": rng.choice(["Adelie", "Chinstrap", "Gentoo"], n),
        "flipper_length_mm": flipper,
        "body_mass_g": 100 * flipper + rng.normal(0, 500, n),
    })

Frequentist vs Bayesian

Frequentist: parameters are fixed unknowns; CIs describe the long-run procedure; p-values test H0. Bayesian: parameters are uncertain random variables. Posterior ∝ likelihood × prior: $$p(\theta \mid y) = \frac{p(y \mid \theta) , p(\theta)}{p(y)}$$ Key outputs: the posterior p(θ|y); the HDI credible interval (a direct probability statement, unlike a CI); and the posterior predictive distribution for model checking.

Goal: compare an OLS estimate of a slope to its Bayesian posterior on the same standardized data. Approach: fit smf.ols for the point estimate/CI, then fit an equivalent PyMC linear model and sample with NUTS.

def fit_bayesian_linear(df, x_col, y_col, draws=1000, tune=500, seed=42):
    """Fit y ~ x with weakly-informative priors on standardized data.

    Returns (ols_result, idata) so frequentist and Bayesian estimates
    can be compared side by side.
    """
    std = (df[[x_col, y_col]] - df[[x_col, y_col]].mean()) / df[[x_col, y_col]].std()
    X, y = std[x_col].values, std[y_col].values

    ols = smf.ols(f"{y_col} ~ {x_col}", data=std).fit()

    with pm.Model() as model:
        alpha = pm.Normal("alpha", mu=0, sigma=2)
        beta = pm.Normal("beta", mu=0, sigma=1)
        sigma = pm.HalfNormal("sigma", sigma=1)
        mu = alpha + beta * X
        pm.Normal("obs", mu=mu, sigma=sigma, observed=y)
        idata = pm.sample(draws, tune=tune, target_accept=0.9,
                           progressbar=False, random_seed=seed)
    return ols, idata

ols, idata_lin = fit_bayesian_linear(penguins, "flipper_length_mm", "body_mass_g")
print(f"OLS beta = {ols.params.iloc[1]:.3f}, 95% CI: {ols.conf_int().iloc[1].values.round(3)}")
print(az.summary(idata_lin, var_names=["alpha", "beta", "sigma"], hdi_prob=0.95)
      [["mean", "hdi_2.5%", "hdi_97.5%", "r_hat"]].round(3))

Prior Specification

Prior predictive check: sample from the prior alone (no data) and confirm simulated values are physically plausible (e.g. no negative body masses) before fitting to real data.

def prior_predictive_check(X, n_samples=200, seed=42):
    """Simulate data from the prior only, to sanity-check prior scale."""
    with pm.Model():
        alpha = pm.Normal("alpha", mu=0, sigma=1)
        beta = pm.Normal("beta", mu=0, sigma=0.5)
        sigma = pm.HalfNormal("sigma", sigma=0.5)
        mu = alpha + beta * X
        pm.Normal("obs", mu=mu, sigma=sigma)
        prior = pm.sample_prior_predictive(n_samples, random_seed=seed)
    return prior.prior_predictive["obs"].values.reshape(-1, len(X))

pp_obs = prior_predictive_check(np.zeros(50))
print(f"Prior predictive range: [{pp_obs.min():.2f}, {pp_obs.max():.2f}] (should cover plausible data range)")

Multiple Regression and Collinearity

VIF > 5-10 flags a collinearity concern: VIF_j = 1 / (1 - R²_j), where R²_j comes from regressing predictor j on all other predictors.

from statsmodels.stats.outliers_influence import variance_inflation_factor

multi_std = (penguins[["body_mass_g", "flipper_length_mm"]]
             .assign(bill_length_mm=penguins.get("bill_length_mm", penguins["flipper_length_mm"] * 0.2))
             .pipe(lambda d: (d - d.mean()) / d.std()))
X_mat = sm.add_constant(multi_std[["flipper_length_mm", "bill_length_mm"]])
vif = pd.DataFrame({
    "feature": X_mat.columns,
    "VIF": [variance_inflation_factor(X_mat.values, i) for i in range(X_mat.shape[1])],
})
print(vif.to_string(index=False))

Model Comparison: LOO-CV and WAIC

Both are estimates of expected log predictive density (ELPD), computed from the pointwise log-likelihood. LOO-CV uses Pareto-smoothed importance sampling (more robust); WAIC is faster but more sensitive to influential points. Only compare models fit to the same observed data with pm.sample_posterior_predictive / log-likelihood stored.

comp = az.compare({"flipper_only": idata_m1, "flipper_plus_bill": idata_m2}, ic="loo")
print(comp[["elpd_loo", "p_loo", "d_loo", "weight"]].round(2))
az.plot_compare(comp, insample_dev=False)
plt.tight_layout(); plt.show()

Hierarchical / Mixed-Effects Models (Bambi)

y ~ x + (1|group) adds a random intercept per group; y ~ x + (x|group) adds a random slope. Random effects use partial pooling — sparse groups borrow strength from the population estimate, shrinking noisy per-group estimates toward the mean.

m_me = bmb.Model("body_mass_g ~ flipper_length_mm + (1|species)",
                  data=penguins, family="gaussian")
idata_me = m_me.fit(draws=800, tune=400, target_accept=0.9,
                     progressbar=False, random_seed=42)
print(az.summary(idata_me, var_names=["flipper_length_mm", "Intercept"], hdi_prob=0.95)
      [["mean", "hdi_2.5%", "hdi_97.5%", "r_hat"]].round(3))
az.plot_forest(idata_me, var_names=["1|species"], combined=True)
plt.tight_layout(); plt.show()

GLM Families for Count Data

def fit_count_glm(counts, x_std, overdispersed=False, seed=1):
    """Fit a Poisson or Negative Binomial GLM: log(mu) = a + b*x."""
    with pm.Model() as model:
        a = pm.Normal("a", 0, 2)
        b = pm.Normal("b", 0, 1)
        mu = pm.math.exp(a + b * x_std)
        if overdispersed:
            alpha = pm.Exponential("alpha", lam=1)
            pm.NegativeBinomial("y", mu=mu, alpha=alpha, observed=counts)
        else:
            pm.Poisson("y", mu=mu, observed=counts)
        idata = pm.sample(800, tune=400, progressbar=False, random_seed=seed)
    return idata

idata_nb = fit_count_glm(counts_negbin, X_temp, overdispersed=True)

Convergence Diagnostics

Always check before trusting a posterior: R-hat should be < 1.01 (chains agree); bulk/tail ESS should be > 400 (enough effective samples); divergences should be ~0 (sampler didn't hit pathological geometry).

summ = az.summary(idata_lin, var_names=["alpha", "beta", "sigma"])
print(summ[["mean", "sd", "ess_bulk", "ess_tail", "r_hat"]].round(3))
n_divergences = int(idata_lin.sample_stats["diverging"].sum())
print(f"Divergences: {n_divergences} (should be 0 or near-0)")
az.plot_trace(idata_lin, var_names=["alpha", "beta", "sigma"], compact=True)
plt.tight_layout(); plt.show()

R Equivalent (brms)

R practitioners typically fit the same hierarchical Bayesian models with brms (a Stan/NUTS front end with the same lme4-style formula syntax as Bambi):

library(brms)
library(bayesplot)

fit <- brm(
  body_mass_g ~ flipper_length_mm + (1 | species),
  data = penguins, family = gaussian(),
  chains = 4, iter = 2000, warmup = 1000, cores = 4, seed = 42
)
summary(fit)          # includes Rhat, Bulk_ESS, Tail_ESS per parameter
mcmc_trace(fit)        # trace plots for convergence
loo(fit)               # LOO-CV, comparable to az.compare in Python

Pitfalls

• Divergences during sampling: indicate the sampler hit pathological posterior geometry (often a "funnel" in hierarchical models). Increase target_accept (e.g. 0.95-0.99) or use a non-centered parameterization before trusting the result.
• Ignoring R-hat/ESS: an R-hat > 1.01 or ESS < 400 means the posterior estimate is unreliable — do not report means/HDIs from a non-converged chain.
• Flat/improper priors: often cause slow mixing or non-identifiability, especially in hierarchical models; use weakly informative priors instead.
• Comparing LOO/WAIC across models fit to different data or likelihoods: az.compare is only valid when all models share the same observed data and response distribution.
• Skipping the prior predictive check: a prior that generates implausible data (negative counts, masses) silently biases the posterior, especially with small samples.
• VIF only flags linear collinearity: it won't catch interaction or nonlinear collinearity — inspect pairwise correlations too.

See Also

• pymc — lower-level PyMC modeling API reference
• statistical-analysis — general frequentist statistical testing
• bio-experimental-design-power-analysis — sample size/power before running a study
• bio-machine-learning-model-validation — cross-validation for predictive (non-Bayesian) models