Bayesian Workflow

Bayesian Workflow

Workflow overview

Every Bayesian analysis follows this sequence. Do not skip steps -- especially model criticism.

1. Formulate — Define the generative story. What underlying process, that we're precisely trying to model, created the data? 2. Specify priors — See references/priors.md 3. Implement in PyMC — Write the model. Prefer PyMC 5+ syntax. Use the latest version possible. 4. Run prior predictive checks — pm.sample_prior_predictive(). Verify priors produce plausible data ranges before fitting 5. Inference — pm.sample(nuts_sampler="nutpie"). Always use nutpie for speed (the nutpie python package provides cutting-edge sampling). Don't hardcode the number of chains — let the sampler pick the best default for the platform. 6. Diagnose convergence — Use arviz_stats.diagnose(idata) as the first check (requires arviz-stats >= 1.0.0). It covers R-hat, ESS, divergences, tree depth, and E-BFMI in one call. See references/diagnostics.md 7. Criticize the model — See references/model-criticism.md 8. Check prior sensitivity — Run psense_summary(idata) to verify conclusions are robust to prior choices. Visualize with plot_psense_dist(idata) from arviz_plots. Requires log_likelihood and log_prior in the InferenceData — compute them after sampling if needed. See references/sensitivity.md 9. Compare models (if applicable) — See references/model-comparison.md 10. Report results — See references/reporting.md. When the user asks for a report or mentions a non-technical audience, generate a standalone markdown report file (not just code comments) using the template in reporting.md. Adapt the language to the audience — if they're new to Bayesian stats, include a glossary and plain-language explanations of key concepts.

Installation

Prefer conda-forge / mamba-forge to install PyMC and its dependencies — pip can cause issues with compiled backends (nutpie, JAX). Example:

mamba install -c conda-forge pymc nutpie arviz arviz-stats preliz

PyMC model template

import pymc as pm
import arviz as az
import numpy as np

RANDOM_SEED = sum(map(ord, "churn-logistic-v1"))
rng = np.random.default_rng(RANDOM_SEED)

# always use dimensions and coordinates in PyMC models
with pm.Model(coords=coords) as model:
    # use Data containers when working on a PyMC model
    data = pm.Data("data", df["y"].to_numpy(), dims="obs")

    # --- Priors ---
    # Always document WHY each prior was chosen
    mu = pm.Normal("mu", mu=0, sigma=10)  # Weakly informative: allows wide range

    # --- Data model ---
    pm.Normal("obs", mu=mu, sigma=1, observed=data, dims="obs")

    # --- Prior predictive check ---
    prior_pred = pm.sample_prior_predictive(random_seed=rng)

    # --- Inference ---
    idata = pm.sample(nuts_sampler="nutpie", random_seed=rng)
    idata.extend(prior_pred)

    # --- Posterior predictive check ---
    idata.extend(pm.sample_posterior_predictive(idata, random_seed=rng))

    # --- Compute log-likelihood and log-prior for sensitivity checks & LOO ---
    pm.compute_log_likelihood(idata, model=model)
    pm.compute_log_prior(idata, model=model)

    # --- Save immediately after sampling ---
    # Late crashes can destroy valid results. Save to disk before any post-processing.
    idata.to_netcdf("model_output.nc")

Critical rules

• Always run prior predictive checks before sampling. If prior predictions span implausible ranges, fix priors first. If you have issues or doubts for some parameters, use the PreliZ package to elicit priors from the user.
• Always check convergence before interpreting results. R-hat > 1.01 or ESS < 100 * nbr_chains means the results are unreliable.
• Always run posterior predictive checks. A model that fits well numerically but cannot reproduce the data is useless.
• Always run calibration checks (PIT / coverage). Use ArviZ's plot_ppc_pit for this — it handles all data types (continuous, binary, count) correctly. See references/model-criticism.md.
• Document every prior choice with a brief justification in a code comment.
• Never report point estimates alone. Always include credible intervals (default: 94% HDI).
• Use `arviz_stats.diagnose(idata)` as the first diagnostic on every model (arviz-stats >= 1.0.0). It checks R-hat, ESS, divergences, tree depth saturation, and E-BFMI in one call. Follow up with az.plot_trace(idata, kind="rank_vlines") for visual inspection.
• Don't hardcode number of chains. Let PyMC / nutpie choose the optimal default for the user's platform. Just call pm.sample() without specifying chains.
• Use reproducible, descriptive seeds. Never use magic numbers like 42. Instead, derive a seed from the analysis name: RANDOM_SEED = sum(map(ord, "my-analysis-name")). Pass it to pm.sample(random_seed=rng), pm.sample_prior_predictive(random_seed=rng), and numpy via rng = np.random.default_rng(RANDOM_SEED).
• Save InferenceData immediately after sampling with idata.to_netcdf("model_output.nc"). Late crashes or kernel restarts can destroy valid MCMC results — save before any post-processing.
• Use ArviZ for all plots and calibration. Don't write custom plotting code when ArviZ already handles it — including for binary data, count data, and calibration. ArviZ developers have thought through edge cases so you don't have to.
• Prefer xarray over numpy for InferenceData operations. InferenceData and DataTree objects are backed by xarray — use xarray's labeled indexing (.sel(), .mean(dim=...), etc.) instead of converting to numpy arrays. This preserves dimension labels, avoids shape bugs, and makes code more readable. Fall back to numpy only when xarray can't do what you need.
• Always generate analysis notes alongside code. When producing a model script, also produce a companion markdown file (analysis_notes.md or similar) that interprets the results — what the diagnostics mean, what the posteriors tell us, what the calibration plots show. Code without interpretation is incomplete.
• Always use the posterior mean (not median) for predictive probabilities. The proper Bayesian predictive distribution averages over the posterior: P(Y=k|x) = (1/S) Σ P(Y=k|x,θₛ). This is the mean, not the median. The median does not correspond to the posterior predictive distribution, can violate probability coherence (probabilities may not sum to 1), and biases calibration due to Jensen's inequality. In code: use np.mean(probs, axis=sample_axis), never np.median(...).
• Use `pm.set_data()` + `pm.sample_posterior_predictive()` for out-of-sample predictions. Don't manually extract posterior samples and recompute predictions — let PyMC propagate uncertainty properly. Define predictors as pm.Data(...) during model building, then swap in new data:

# After fitting the model:
with model:
    pm.set_data({"X": X_new, "group_idx": group_idx_new})
    oos_preds = pm.sample_posterior_predictive(idata, predictions=True, random_seed=rng)

• Check model identifiability before interpreting components. If two model components always appear together in the likelihood (e.g., a league intercept and a home advantage term when every observation is from home perspective), their individual posteriors reflect prior assumptions, not data signal — only their sum is identified. Use az.plot_pair() to check for strong posterior correlations between components. If correlation is near ±1, the components are not separately identifiable — either merge them or restructure the data.

Common model families

Utility scripts

Run diagnose_model.py after sampling to get a structured convergence + diagnostics report:

python scripts/diagnose_model.py --idata path/to/inference_data.nc

Run calibration_check.py to generate calibration plots:

python scripts/calibration_check.py --idata path/to/inference_data.nc

See scripts/ for all available utilities.

Common gotchas

These are battle-tested lessons that save hours of debugging:

• nutpie silently ignores `idata_kwargs` for log_likelihood and log_prior. Always compute them explicitly after sampling: pm.compute_log_likelihood(idata, model=model) (needed for LOO-CV) and pm.compute_log_prior(idata, model=model) (needed for prior sensitivity checks). Don't assume they're stored automatically.
• `az.plot_khat()` requires the LOO object, not InferenceData. Pass the output of az.loo(idata, pointwise=True) to it.
• Flat priors on scale parameters (HalfCauchy, HalfFlat) cause funnels in hierarchical models. Use Gamma(2, ...) or Exponential — these avoid the near-zero region that creates sampling problems. If there's no group-level variation to detect, you don't need the hierarchy.
• Python conditionals in models (if x > 0) don't work inside PyMC. Use pm.math.switch or pytensor.tensor.where instead.
• Forgetting to standardize predictors makes shared priors inappropriate and slows sampling. Always standardize before fitting, then back-transform for interpretation.
• Horseshoe priors create a double-funnel geometry that standard NUTS can struggle with. Always use the regularized (Finnish) horseshoe (Piironen & Vehtari, 2017), which adds a slab component that smooths the geometry. Set target_accept=0.95 or higher. If you see divergences with a horseshoe model, this is almost certainly the cause.
• `np.median` on posterior predictive probabilities is a silent bug. It does not produce the Bayesian predictive distribution and can yield probabilities that don't sum to 1 across categories. Always use np.mean over the posterior samples dimension.

When things go wrong

bayesian-workflow

An opinionated Agent Skill for building, diagnosing, and reporting on Bayesian statistical models using PyMC and ArviZ.

Compatible with Claude Code, Kimi Code, Cursor, Gemini CLI, and any agent that supports the Agent Skills spec.

Full breakdown here.

What it does

Guides your coding agent through the full Bayesian workflow:

1. Formulate the generative story 2. Specify priors with documented justifications 3. Implement in PyMC using modern best practices (coords, dims, nutpie) 4. Prior predictive checks before fitting 5. Inference via MCMC (nutpie sampler by default) 6. Convergence diagnostics (R-hat, ESS, trace plots, energy plots) 7. Model criticism (posterior predictive checks, LOO-PIT calibration) 8. Prior/likelihood sensitivity (power-scaling via PSIS) 9. Model comparison (LOO-CV, ELPD, stacking weights) 10. Reporting with companion analysis notes and audience-adapted reports

The skill enforces guardrails that agents won't apply on their own: 94% HDI, mandatory calibration checks, prior/likelihood sensitivity checks, non-centered parameterizations for hierarchical models, reproducible descriptive seeds, immediate save-to-disk after sampling, and xarray-first data manipulation.

Install

Claude Code

Clone and copy the skill into your personal skills directory:

git clone https://github.com/Learning-Bayesian-Statistics/baygent-skills.git /tmp/baygent-skills
mkdir -p ~/.claude/skills
cp -r /tmp/baygent-skills/bayesian-workflow ~/.claude/skills/

For project-level installation (available only in that project), copy into .claude/skills/ at the project root instead.

Other compatible agents (Kimi Code, Cursor, etc.)

Clone the repo and copy the skill folder into your agent's skills directory:

git clone https://github.com/Learning-Bayesian-Statistics/baygent-skills.git /tmp/baygent-skills
cp -r /tmp/baygent-skills/bayesian-workflow/ ~/.config/agents/skills/bayesian-workflow/

PyMC installation

The skill recommends conda-forge / mamba-forge for PyMC and its dependencies:

mamba install -c conda-forge pymc nutpie arviz arviz-stats preliz

Example prompts

Once installed, just ask your agent naturally:

• "I have customer churn data with binary outcome plus age, tenure, and monthly spend. Build me a Bayesian logistic regression with uncertainty estimates."
• "My model has 47 divergences and R-hat of 1.03. What do I do?"
• "I have test scores for 200 students across 15 schools. Some schools only have 5 students. Help me build a hierarchical model."
• "Compare these two models and tell me which one to use. I have the InferenceData objects."
• "I need to present Bayesian results to my boss who has no stats background."

What's included

bayesian-workflow/
├── SKILL.md                          # Main workflow instructions
├── main.py                           # Entrypoint for programmatic use
├── pyproject.toml                    # Package metadata
├── references/
│   ├── priors.md                     # Prior selection guide
│   ├── diagnostics.md                # Convergence diagnostics
│   ├── model-criticism.md            # PPC, calibration, LOO-PIT
│   ├── model-comparison.md           # LOO-CV, ELPD, stacking weights
│   ├── hierarchical.md               # Partial pooling, non-centered parameterization
│   ├── sensitivity.md                # Prior/likelihood sensitivity analysis (power-scaling)
│   └── reporting.md                  # Report templates, audience adaptation
└── scripts/
    ├── diagnose_model.py             # Post-sampling diagnostics report
    └── calibration_check.py          # Calibration plots from InferenceData

License

MIT - see LICENSE.

Convergence Diagnostics

Contents

• Quick diagnostic checklist
• R-hat
• Effective sample size (ESS)
• Divergences
• Trace plots and rank plots
• Energy diagnostics
• Automated diagnostics workflow

Quick diagnostic checklist

Run this immediately after sampling. If any check fails, do NOT interpret results.

# you just ran this:
# idata = pm.sample(nuts_sampler="nutpie", random_seed=rng)

import arviz_stats as azs

# 1. One-call diagnostics: R-hat, ESS, divergences, tree depth, E-BFMI
# Prints a human-readable summary and returns True if any check failed.
# Defaults already enforce our thresholds (rhat_max=1.01, ess=100*chains, bfmi>=0.3).
has_errors = azs.diagnose(idata)

# 2. If you need the detailed breakdown programmatically:
has_errors, diagnostics = azs.diagnose(idata, return_diagnostics=True, show_diagnostics=False)
# diagnostics contains: "divergent", "treedepth", "bfmi", "ess", "rhat"

# 3. Visual check
az.plot_trace(idata, kind="rank_vlines")

If azs.diagnose is not available (arviz-stats < 1.0.0), fall back to the manual checks below.

R-hat

Measures agreement across chains. Uses the rank-normalized split-R-hat (ArviZ default).

Common causes of high R-hat:

• Insufficient warmup (increase tune)
• Multimodal posterior (reparameterize or use different sampler)
• Model misspecification creating ridges/funnels

Effective sample size (ESS)

The number of independent-equivalent draws. Two flavors matter:

• ESS bulk: Reliability of central tendency estimates (mean, median)
• ESS tail: Reliability of tail estimates (credible intervals, quantiles)

Improving ESS:

• Increase number of draws
• Reparameterize (non-centered for hierarchical models)
• Reduce posterior correlations
• Increase target_accept (trades speed for better exploration)

Divergences

Divergent transitions indicate the sampler encountered regions of high curvature it could not navigate. Even a few divergences (starting from 10+) can bias results.

# Count divergences
n_div = idata.sample_stats["diverging"].sum().item()

# Visualize where divergences occur
# sometimes, when the model is high dimensional and / or has lots of parameters,
# this plot becomes impossible and you want to check just some pairs of potentially
# problematic parameters (e.g population standard deviation in hierarchical models).
az.plot_pair(idata, var_names=["param1", "param2"], divergences=True)

Fix divergences in this order:

1. Increase `target_accept`: pm.sample(target_accept=0.95) — try up to 0.99 2. Reparameterize: Non-centered parameterization for hierarchical models:

# CENTERED (can cause funnel divergences)
mu = pm.Normal("mu", mu=0, sigma=sigma_group)

# NON-CENTERED (usually fixes the funnel)
mu_offset = pm.Normal("mu_offset", mu=0, sigma=1)
mu = pm.Deterministic("mu", mu_offset * sigma_group)

3. Stronger priors on scale parameters: Tight prior on group-level SD can eliminate the funnel, especially avoiding the regions near 0, which don't really mean much in practice anyways (if there is no group-level variation, you don't need to model it!) 4. Marginalize discrete parameters: If possible, integrate out discrete variables analytically

Trace plots and rank plots

# Rank plots (preferred over raw trace plots)
az.plot_trace(idata, kind="rank_vlines")

# What to look for:
# - Rank plots should look uniform (no spikes or gaps)
# - Traces should be "well-mixed" — all chains overlapping
# - No chains stuck in different regions
# - No obvious trends or slow drift

Energy diagnostics

Energy plots detect problems the other diagnostics may miss (e.g., incomplete exploration of the typical set).

az.plot_energy(idata)

# What to look for:
# - Marginal energy and energy transition distributions should overlap
# - Large gap between them indicates poor exploration

Automated diagnostics workflow

The recommended approach is arviz_stats.diagnose() (see Quick diagnostic checklist above). It checks R-hat, ESS, divergences, tree depth saturation, and E-BFMI in one call — with sensible defaults that match our thresholds.

For a script-based workflow, use diagnose_model.py:

python scripts/diagnose_model.py --idata model_output.nc

Or inline (this last approach doesn't require arviz_stats.diagnose()):

def run_diagnostics(idata):
    """Run all convergence diagnostics. Returns dict of results."""
    summary = az.summary(idata)
    num_chains = idata.posterior.sizes["chain"]
    results = {
        "rhat_max": float(summary["r_hat"].max()),
        "rhat_ok": bool((summary["r_hat"] <= 1.01).all()),
        "ess_bulk_min": int(summary["ess_bulk"].min()),
        "ess_tail_min": int(summary["ess_tail"].min()),
        "ess_ok": bool((summary["ess_bulk"] >= 100 * num_chains).all() and (summary["ess_tail"] >= 100 * num_chains).all()),
        "n_divergences": int(idata.sample_stats["diverging"].sum()),
        "divergences_ok": int(idata.sample_stats["diverging"].sum()) == 0,
    }
    results["all_ok"] = results["rhat_ok"] and results["ess_ok"] and results["divergences_ok"]
    return results

Hierarchical (Multilevel) Models

Contents

• When to use hierarchical models
• Partial pooling intuition
• Centered vs. non-centered parameterization
• Common hierarchical structures
• Diagnostics specific to hierarchical models

When to use hierarchical models

Use hierarchical models when data has grouped structure — observations nested within units (students in schools, games in seasons, patients in hospitals, items in categories). Careful: time series data is not hierarchical because timestamps are not interchangeable (they have an order). Tell users that if they have time series data, they should use time series models instead.

This is a great tutorial.

The key question: Do groups share information? If group-level parameters are related (e.g., batting averages across players), hierarchical models borrow strength across groups through partial pooling.

Partial pooling intuition

Three approaches to grouped data:

• Complete pooling: Ignore groups, fit one model. Misses group-level variation. Maximum bias.
• No pooling: Fit separate models per group. Overfits small groups. Maximum variance.
• Partial pooling (hierarchical): Groups share a common distribution. Small groups shrink toward the global mean; large groups retain their own estimate and influence the global population. Trades off worse in-sample coverage for better out-of-sample performance.

Partial pooling is almost always the right choice. It naturally handles imbalanced group sizes.

Centered vs. non-centered parameterization

This is the most common source of divergences in hierarchical models.

# CENTERED — works well when groups have lots of data
# always use dimensions and coordinates in PyMC models
with pm.Model(coords=coords) as centered:
    # use Data containers when working on a PyMC model
    group_idx = pm.Data("group_idx", group_id, dims="obs")
    y = pm.Data("y", df["y"].to_numpy(), dims="obs")
    
    mu_global = pm.Normal("mu_global", mu=0, sigma=10)
    sigma_global = pm.Gamma("sigma_global", 2, 2)
    
    mu_group = pm.Normal("mu_group", mu=mu_global, sigma=sigma_global, dims="group")
    sigma_obs = pm.Gamma("sigma_obs", 2, 2)
    
    # or whatever the data model is
    pm.Normal("obs", mu=mu_group[group_idx], sigma=sigma_obs, observed=y, dims="obs")

# NON-CENTERED — works well when groups have little data
# This is a reparameterization that eliminates the "funnel" geometry
with pm.Model() as noncentered:
    group_idx = pm.Data("group_idx", group_id, dims="obs")
    y = pm.Data("y", df["y"].to_numpy(), dims="obs")
    
    mu_global = pm.Normal("mu_global", mu=0, sigma=10)
    sigma_global = pm.Gamma("sigma_global", 2, 2)
    
    mu_raw = pm.Normal("mu_raw", mu=0, sigma=1, dims="group")
    mu_group = pm.Deterministic("mu_group", mu_global + mu_raw * sigma_global, dims="group")
    
    # or whatever the data model is
    pm.Normal("obs", mu=mu_group[group_idx], sigma=sigma_obs, observed=y, dims="obs")

Rule of thumb: Start with non-centered. Switch to centered only if non-centered shows poor ESS AND groups have substantial data (50+ observations each).

Common hierarchical structures

Varying intercepts

Each group has its own baseline, partially pooled toward a global mean.

# This is the centered parameterization
mu_group = pm.Normal("mu_group", mu=mu_global, sigma=sigma_global, dims="group")
pm.Normal("obs", mu=mu_group[group_idx], sigma=sigma_obs, observed=y, dims="obs")

Varying intercepts and slopes

Each group has its own baseline AND its own effect of a predictor.

# Correlated varying effects (preferred)
# the `n` depends on how many parameters you're relating to each other
# most of the time, it's gonna be intercept + one slope, like here
chol, corr, stds = pm.LKJCholeskyCov(
    "chol", n=2, eta=2.0, sd_dist=pm.Exponential.dist(1.0)
)

# prior for average intercept:
mu_intercept = pm.Normal("mu_intercept", mu=0.0, sigma=5.0)
# prior for average slope:
mu_slope = pm.Normal("mu_slope", mu=0.0, sigma=1.0)

# population of varying effects:
effects = pm.MvNormal(
    "effects", 
    mu=[mu_intercept, mu_slope], 
    chol=chol,
    dims=("group", "param"),
)

# then, continue as usual:
# expected value per group:
# mu_group = effects[group_idx, 0] + effects[group_idx, 1] * slope_data ...

Nested hierarchy

Groups within groups (students in classrooms in schools). Don't go overboard with this, as models become unwieldy and hard to sample and interpret with too many hierarchies.

# School level
mu_school = pm.Normal("mu_school", mu=mu_global, sigma=sigma_global, dims="school")
# Classroom level (nested within school)
mu_class = pm.Normal("mu_class", mu=mu_school[school_idx], sigma=sigma_school, dims="class")
# Student level
pm.Normal("y", mu=mu_class[class_idx], sigma=sigma_student, observed=data, dims="obs")

Diagnostics specific to hierarchical models

In addition to standard diagnostics (reference/diagnostics.md), check:

1. Shrinkage plot: Visualize how much each group is pulled toward the global mean 2. Group-level SD posterior: If sigma_group or sigma_global posterior piles up near zero, the data may not support group-level variation (partial pooling → complete pooling) 3. Funnel plot: Plot group-level means vs. group-level SD. Funnels indicate centered parameterization problems

# Shrinkage plot
group_means_posterior = idata.posterior["mu_group"].mean(dim=["chain", "draw"]).to_numpy()
group_means_obs = [y[group_idx == g].mean() for g in range(n_groups)]

plt.scatter(group_means_raw, group_means_posterior)
plt.plot([min(group_means_obs), max(group_means_obs)],
         [min(group_means_obs), max(group_means_obs)], "r--", label="No pooling")
plt.axhline(np.mean(y), color="gray", linestyle=":", label="Complete pooling")
plt.xlabel("Observed group mean")
plt.ylabel("Posterior group mean")
plt.legend()
plt.title("Shrinkage toward global mean")

---

Identifiability checks

A model component is identifiable only if the data can distinguish its effect from other components. In hierarchical models, identifiability failures are common and subtle — the model samples fine, diagnostics pass, but individual component posteriors reflect prior assumptions, not data signal.

Common identifiability traps

1. Separate intercept + offset when the offset is always active: If every observation has a characteristic (e.g., every match row is from the home team's perspective), you cannot separately estimate a "baseline" intercept and a "home advantage" offset — only their sum is identified. Merge them into a single intercept.

2. Overparameterized intercepts: Group-level intercepts + a global intercept without a sum-to-zero constraint creates a "trading" pattern where the global intercept can shift arbitrarily as group intercepts compensate. Use pm.ZeroSumNormal for group intercepts or drop the global intercept.

3. Collinear covariates at the group level: If a group-level predictor is nearly constant within groups (e.g., all teams in a league share the same league-level feature), it is confounded with the group intercept.

How to detect

# Check posterior correlations between suspect components
az.plot_pair(
    idata,
    var_names=["alpha", "delta"],
    divergences=True,
    marginals=True,
)
# If correlation is near ±1 → components are not separately identifiable

What to do

• Merge confounded components into a single term (honest about what the data can tell you)
• Add identifying variation to the data (e.g., include away-team observations to separate home advantage from league baseline)
• Accept the constraint: if you need both components for interpretability, acknowledge that their individual posteriors are prior-driven and report only their sum

Model Comparison

Contents

• When to compare models
• LOO-CV comparison
• Stacking weights
• WAIC (and when to prefer LOO)
• Reporting comparisons

When to compare models

Compare models when you have genuinely different modeling assumptions — not for variable selection. Bayesian model comparison answers: "Which model predicts unseen data better?" LOO-CV works across different data models — they do not need to share the same observation distribution (see CV-FAQ).

For variable selection, use a BART model, or the kulprit package.

Common comparison scenarios:

• Linear vs. nonlinear trend
• Different hierarchical structures (varying intercepts vs. varying slopes)
• Different covariate sets guided by domain knowledge

LOO-CV comparison

The primary comparison tool. Uses PSIS-LOO via ArviZ.

# Fit multiple models, store InferenceData for each
models = {"m1": idata_1, "m2": idata_2, "m3": idata_3}

# you need to store the pointwise log-likelihood in the InferenceData for each model,
# either upstream with `idata_kwargs={"log_likelihood": True}` when sampling,
# or downstream with:
with m1:
    pm.compute_log_likelihood(idata_1)
with m2:
    pm.compute_log_likelihood(idata_2)
with m3:
    pm.compute_log_likelihood(idata_3)

# Compare
comparison = az.compare(models)
print(comparison)

# Visualize
az.plot_compare(comparison)

Reading the comparison table:

• elpd_loo: Higher is better (less negative = better predictive accuracy)
• se: Standard error of ELPD estimate
• elpd_diff: Difference from best model
• dse: Standard error of the difference
• weight: Stacking weight (see below)
• warning: True if high Pareto k values exist

Interpreting differences:

• If elpd_diff < 2×dse → Models are practically indistinguishable. Prefer the simpler one.
• If elpd_diff > 4×dse → Strong evidence for the better model.
• Between 2–4×dse → Moderate evidence. Consider domain knowledge.

Stacking weights

When no single model is clearly best, use Bayesian stacking to combine predictions:

comparison = az.compare(models, method="stacking")
# The 'weight' column gives optimal combination weights

Stacking weights minimize expected log predictive density loss. They often outperform selecting a single best model. Report stacking weights alongside ELPD differences — they give a more nuanced picture.

WAIC (and when to prefer LOO)

WAIC (Widely Applicable Information Criterion) is asymptotically equivalent to LOO but less robust in practice.

waic = az.waic(idata)

Prefer LOO over WAIC in most cases because:

• LOO provides the Pareto k diagnostic (you know when to trust it)
• WAIC can silently give unreliable results with no warning
• LOO is better calibrated for small samples

Use WAIC only when LOO is computationally infeasible (very rare with PSIS-LOO).

Reporting comparisons

When reporting model comparisons, always include:

1. Table of ELPD values with standard errors 2. ELPD differences with their standard errors 3. Stacking weights 4. Note any high Pareto k warnings, what they mean, and what to do about it 5. The substantive interpretation — what does the better model imply about the phenomenon? Be careful to NOT make causal claims based on model comparison -- it only tells us about predictive accuracy.

Template:

## Model comparison

| Model | ELPD (LOO) | SE | ΔELPD | ΔSE | Weight |
|-------|------------|-----|-------|------|--------|
| Model 1 | -234.5 | 12.3 | 0.0 | — | 0.72 |
| Model 2 | -241.2 | 11.8 | -6.7 | 3.1 | 0.28 |

Model 1 is slightly preferred by LOO (ΔELPD = 6.7, ~2.2× the SE of the difference),
suggesting moderate evidence. Consider domain knowledge to choose one model, 
or use stacking weights to combine predictions.
No observations had Pareto k > 0.7.

Model Criticism

Model criticism answers: "Is this model any good?" Convergence diagnostics (reference/diagnostics.md) only tell you the sampler worked -- they say nothing about whether the model is appropriate for the data.

Contents

• Posterior predictive checks (PPC)
• Leave-one-out cross-validation (LOO-CV)
• Calibration assessment
• Simulation-based calibration (SBC)
• Residual analysis
• Decision workflow

Posterior predictive checks (PPC)

The most important model criticism tool. Simulate data from the fitted model and compare to observed data.

with model:
    idata.extend(pm.sample_posterior_predictive(idata, random_seed=rng))

# Visual check: do simulated datasets resemble the real data?
az.plot_ppc(idata)

What to look for:

• Posterior predictive distribution should envelop the observed data
• Check shape, spread, and key features (skewness, multimodality, tails)
• Systematic departures indicate model misspecification

Targeted PPCs — Check specific data features the model should capture, using az.plot_ppc to isolate parameters of interest (e.g the posterior standard deviation, to check if model captures the observed variance).

Choose test statistics relevant to your problem: mean, variance, skewness, max, proportion above threshold, autocorrelation (for time series), etc.

Leave-one-out cross-validation (LOO-CV)

Estimates out-of-sample predictive accuracy using Pareto-smoothed importance sampling (PSIS-LOO). This is the primary tool for model comparison but also useful for single-model criticism.

loo = az.loo(idata, pointwise=True)
print(loo)

Key outputs:

• elpd_loo: Expected log pointwise predictive density. Higher (less negative) is better.
• p_loo: Effective number of parameters. If p_loo >> actual parameter count, the model may be misspecified or the priors are too weak.
• pareto_k: Per-observation diagnostic. Flags influential or poorly-fit observations.

Pareto k diagnostic — Critical for trusting LOO results:

# Find problematic observations
pareto_k = loo.pareto_k.values
bad_obs = np.where(pareto_k > 0.7)[0]
print(f"Observations with high Pareto k: {bad_obs}")

# Visualize
az.plot_khat(loo)

High Pareto k observations are often outliers or observations the model fits poorly. Investigate them — they may reveal model misspecification.

Calibration assessment

Calibration is mandatory for every model, not optional. A well-calibrated model's X% credible intervals should contain the true value about X% of the time. Run this even for binary and count data — ArviZ handles all data types correctly.

How to run calibration

Always use ArviZ for calibration plots. Don't write custom calibration code — ArviZ's plot_ppc_pit handles continuous, binary, and count data correctly out of the box:

# ArviZ 1.0+ (arviz_plots)
import arviz_plots as azp

# PPC-PIT: compares posterior predictive to observed
azp.plot_ppc_pit(idata)

# LOO-PIT: leave-one-out calibration (more robust, preferred when LOO is available)
azp.plot_ppc_pit(idata, loo_pit=True)

Refer to this guide for detailed coverage interpretation — it's a treasure trove for the whole Bayesian workflow.

Coverage calibration

Interpretation:

• If empirical coverage ≈ nominal → well-calibrated
• If the difference is positive, the model is under-confident: the predictions have a wider spread than the data – they are too uncertain.
• If the difference is negative, the model is over-confident: the predictions have a narrower spread than the data – they are too certain.

PIT histograms (probability integral transform)

A sharper calibration check. If the model is calibrated, PIT values should be uniform. Refer to this section for how to do it, using the new ArviZ.

Patterns:

• U-shaped → underdispersed (intervals too narrow)
• Inverted-U → overdispersed (intervals too wide)
• Skewed → systematic bias in location
• Uniform → well-calibrated

Simulation-based calibration (SBC)

SBC validates that the entire inference pipeline is correct — priors, data model, sampler, and code. It simulates data from the prior, fits the model, and checks that posterior rank statistics are uniform.

This is the gold standard for validating a new model implementation. Run it once per model specification, if you have doubts about the model, since SBC is computationally expensive.

Use the simuk package, either directly, or as inspiration to adapt your own code.

Interpretation:

• Uniform ranks → inference pipeline is correct
• Systematic patterns → implementation bug, wrong prior, sampler failure
• SBC failures mean the model code has a bug — fix before interpreting results

When to run SBC:

• Developing a new model you'll reuse
• Complex hierarchical models where bugs are easy to introduce
• Custom data models or potentials
• Not necessary for routine analyses with standard model families

Residual analysis

For regression-style models, check residuals for patterns:

# Posterior predictive mean
pp_mean = idata.posterior_predictive["obs"].mean(dim=["chain", "draw"])
residuals = y_obs - pp_mean

# Residuals vs. fitted
plt.scatter(pp_mean, residuals, alpha=0.5)
plt.axhline(0, color="red", linestyle="--")
plt.xlabel("Fitted values")
plt.ylabel("Residuals")
plt.title("Residuals vs. Fitted")

# Residuals vs. predictors (check for missed nonlinearity)
for j, name in enumerate(predictor_names):
    plt.figure()
    plt.scatter(X[:, j], residuals, alpha=0.5)
    plt.axhline(0, color="red", linestyle="--")
    plt.xlabel(name)
    plt.ylabel("Residuals")

Look for: trends (missed nonlinearity), fans (heteroscedasticity), clusters (missing grouping variable).

Classification and ordinal model evaluation

Standard PPC and calibration checks apply to classification models — always run `plot_ppc_pit` first (see Calibration assessment above). The metrics below supplement PPC-PIT with classification-specific numeric summaries. Note: sklearn.metrics.brier_score_loss exists but is binary-only; there is no standard package for multiclass ECE or categorical RPS, so we provide lightweight helpers:

Metrics for categorical/ordinal outcomes

def expected_calibration_error(pred_probs, actuals, n_bins=10):
    """Confidence-based ECE: are predicted probabilities well-calibrated?"""
    # Standard ECE: bin on the model's confidence (max predicted probability),
    # compare to accuracy within each bin.
    confidences = np.max(pred_probs, axis=1)
    predictions = np.argmax(pred_probs, axis=1)
    accuracies = (predictions == actuals).astype(float)
    bin_edges = np.linspace(0, 1, n_bins + 1)
    ece = 0
    for i in range(n_bins):
        in_bin = (confidences >= bin_edges[i]) & (confidences < bin_edges[i + 1])
        if in_bin.sum() > 0:
            avg_confidence = confidences[in_bin].mean()
            avg_accuracy = accuracies[in_bin].mean()
            ece += np.abs(avg_confidence - avg_accuracy) * in_bin.sum()
    return ece / len(actuals)

def ranked_probability_score(pred_probs, actuals, n_classes):
    """RPS: gold standard for ordinal outcomes. Penalizes being 'far off' more than Brier."""
    rps = 0
    for i, actual in enumerate(actuals):
        pred_cum = np.cumsum(pred_probs[i])
        actual_cum = np.zeros(n_classes)
        actual_cum[int(actual):] = 1
        rps += np.sum((pred_cum - actual_cum) ** 2)
    return rps / len(actuals)

Per-class calibration plots

For classification models, always check calibration per class, not just overall. A model can be well-calibrated on average but poorly calibrated for specific outcomes (e.g., good at predicting home wins but overconfident on draws).

fig, axes = plt.subplots(1, n_classes, figsize=(5 * n_classes, 4))
for k, ax in enumerate(axes):
    pred_k = pred_probs[:, k]
    actual_k = (actuals == k).astype(float)
    # Bin predictions and compute observed frequency
    bin_edges = np.linspace(0, 1, 11)
    bin_centers, bin_means = [], []
    for i in range(10):
        mask = (pred_k >= bin_edges[i]) & (pred_k < bin_edges[i + 1])
        if mask.sum() > 5:
            bin_centers.append(pred_k[mask].mean())
            bin_means.append(actual_k[mask].mean())
    ax.plot([0, 1], [0, 1], "k--", alpha=0.3)
    ax.scatter(bin_centers, bin_means)
    ax.set_title(f"Class {k}")
    ax.set_xlabel("Predicted P")
    ax.set_ylabel("Observed frequency")

Key metrics summary

Temporal and out-of-sample evaluation

For panel data or time-series-adjacent data (e.g., multiple seasons), always evaluate per time period:

for period in sorted(data_oos["season"].unique()):
    mask = data_oos["season"] == period
    period_metrics = evaluate(pred_probs[mask], actuals[mask])
    print(f"{period}: {period_metrics}")

This reveals temporal degradation — a model that works well on 2020 but poorly on 2023 may be overfitting to historical patterns. If you see degradation, consider whether the data-generating process has changed (concept drift) or whether the training window needs expanding.

Feature importance for shrinkage models

When using sparsity priors (horseshoe, R2-D2), summarize feature relevance via probability of practical significance:

beta_samples = idata.posterior["beta"].stack(samples=("chain", "draw")).values
threshold = 0.05  # on the standardized coefficient scale

importance = pd.DataFrame({
    "feature": features,
    "posterior_mean": beta_samples.mean(axis=-1),
    "posterior_sd": beta_samples.std(axis=-1),
    "P(|beta|>threshold)": (np.abs(beta_samples) > threshold).mean(axis=-1),
}).sort_values("P(|beta|>threshold)", ascending=False)

This is more informative than just looking at posterior means — it tells you the probability that each feature has a practically meaningful effect, which is the natural Bayesian answer to "which features matter?"

Decision workflow

After running diagnostics:

1. Convergence OK?  (reference/diagnostics.md)
   NO  → Fix sampler issues first. Do NOT proceed.
   YES ↓

2. Posterior predictive check pass?
   NO  → Model misspecification. Revise data model or add complexity.
   YES ↓

3. LOO-CV: any high Pareto k?
   YES → Investigate flagged observations. Consider K-fold CV.
   NO  ↓

4. Calibration OK?  (coverage + PIT)
   NO  → Model is mis-calibrated. Check priors, data model, missing predictors.
   YES ↓

5. Residual patterns?
   YES → Missing structure. Add predictors, nonlinearity, or hierarchical effects.
   NO  ↓

→ Model is ready for interpretation and reporting.

Prior Selection Guide

Contents

• Philosophy: why priors matter
• Weakly informative priors (the default)
• Prior families by parameter type
• Sparsity priors (high-dimensional regression)
• Prior predictive checking workflow
• Common mistakes

Philosophy: why priors matter

Priors encode domain knowledge and constrain the model to plausible regions of parameter space. The goal is NOT to be "non-informative" — it is to be honestly informative while avoiding undue influence on the posterior when data is sufficient.

Every prior should have a justification. If you cannot articulate why a prior is reasonable, it is not a good prior.

Weakly informative priors (the default)

When in doubt, use weakly informative priors. These place most mass on plausible values while still allowing the data to dominate.

Principle: A good weakly informative prior rules out nonsense values but does not strongly favor any particular reasonable value.

The PreliZ package is your friend when it comes to choosing priors.

Here are some general rules of thumb:

Regression coefficients

# if you standardize predictors first -- makes priors comparable
beta = pm.Normal("beta", mu=0, sigma=2.5)  # On standardized scale

# if on raw scale with known range: sigma ≈ expected_range / 4
beta_raw = pm.Normal("beta_raw", mu=0, sigma=(plausible_max - plausible_min) / 4)

Scale parameters (standard deviations)

# Gamma avoiding near-zero values: good default
sigma = pm.Gamma("sigma", alpha=2, beta=2)

# Exponential: when you want to favor smaller values
sigma = pm.Exponential("sigma", lam=1)

Intercepts

# Center on observed data mean when possible
intercept = pm.Normal("intercept", mu=y_mean, sigma=2 * y_std)

Correlation matrices (hierarchical models)

# LKJ: the standard choice. eta=1 is uniform, eta=2 pulls toward identity
chol, corr, stds = pm.LKJCholeskyCov(
    "chol", n=k, eta=2.0, sd_dist=pm.Exponential.dist(1.0)
)

Prior families by parameter type

Sparsity priors (high-dimensional regression)

When you have many features and expect only a subset to be relevant, use a sparsity-inducing prior instead of a shared Normal. These adaptively shrink irrelevant coefficients toward zero while preserving signal from important ones.

When to use sparsity priors

Regularized Horseshoe (Finnish Horseshoe)

The go-to sparsity prior. Shrinks irrelevant features toward zero ("spike") while allowing strong signals through ("slab"). Always use the regularized variant (Piironen & Vehtari, 2017) — the original horseshoe has a double-funnel geometry that causes divergences.

import pytensor.tensor as pt

D = X.shape[1]  # number of features
N = X.shape[0]  # number of observations
p0 = 5  # prior guess for number of relevant features

# Global shrinkage — controls overall sparsity
tau = pm.HalfStudentT("tau", nu=2, sigma=p0 / (D - p0) / np.sqrt(N))

# Local shrinkage — per-feature
lam = pm.HalfStudentT("lam", nu=5, shape=D)

# Slab — regularizes large coefficients (prevents double-funnel)
c2 = pm.InverseGamma("c2", alpha=1, beta=1)

# Effective shrinkage
lam_tilde = pt.sqrt(c2 * lam**2 / (c2 + tau**2 * lam**2))

# Coefficients
beta_raw = pm.Normal("beta_raw", mu=0, sigma=1, shape=D)
beta = pm.Deterministic("beta", beta_raw * tau * lam_tilde, dims="feature")

Key hyperparameter: p0 (expected number of relevant features). This controls the global shrinkage tau. Be honest about your prior belief — setting p0 too high defeats the purpose of sparsity.

Practical tip: With horseshoe models, always use target_accept=0.95 or higher.

Feature importance with horseshoe: After fitting, compute P(|β| > threshold) per feature to rank practical significance:

beta_samples = idata.posterior["beta"].values  # (chains, draws, D)
threshold = 0.05  # adjust based on scale
prob_relevant = (np.abs(beta_samples) > threshold).mean(axis=(0, 1))

R2-D2 prior

An alternative where you specify prior beliefs about the total R² (variance explained) rather than per-feature shrinkage. More interpretable when you have a prior sense of overall model fit.

import pytensor.tensor as pt  # if not already imported above

# Prior on R² (proportion of variance explained)
R2 = pm.Beta("R2", alpha=1, beta=1)  # uniform on [0, 1] — adjust based on domain

# Concentration parameter controls how R² is distributed across features
# Higher concentration = more uniform; lower = more sparse
phi = pm.Dirichlet("phi", a=np.ones(D), shape=D)

# Coefficient variances derived from R²
sigma2_y = pm.HalfNormal("sigma2_y", sigma=1)  # residual variance
tau2 = R2 / (1 - R2) * sigma2_y * phi

beta = pm.Normal("beta", mu=0, sigma=pt.sqrt(tau2), dims="feature")

Use R2-D2 when:

• You have domain knowledge about how much variance the model should explain
• You want a more interpretable parameterization than the horseshoe
• The horseshoe's funnel causes persistent divergences even with regularization

Prior predictive checking workflow

This is mandatory. Never skip it.

with model:
    prior_pred = pm.sample_prior_predictive()

# Visualize
az.plot_ppc(prior_pred, group="prior", num_pp_samples=100)

# Check: do simulated datasets look plausible?
# - Are values in a reasonable range?
# - Is the spread of outcomes reasonable?
# - Are there impossible values (negative counts, proportions > 1)?

Decision rule:

• If >10% of prior predictive samples are clearly implausible → tighten priors
• If prior predictions are extremely narrow → priors may be too informative, consider loosening
• If prior predictions are reasonable → proceed to inference

Common mistakes

1. Flat / diffuse priors (e.g., Normal(0, 1000)): These are NOT "non-informative". They place excessive mass on extreme, implausible values and can cause sampling issues. Use weakly informative priors instead.

2. Ignoring scale: A Normal(0, 10) prior means very different things depending on the scale of the data. Always consider the units.

3. Forgetting to standardize predictors: Without standardization, coefficients live on different scales, making shared priors inappropriate and slowing sampling. This is not always true, but it is a common mistake.

4. No prior predictive check: The single most common source of modeling errors. Always visualize what your priors imply before fitting.

5. Informative priors without justification: If you use a tight prior, you need a clear reason (previous study, physical constraint, domain expertise). Document it.

Reporting Bayesian Analyses

Contents

• Reporting principles
• Analysis report template
• Presentation template
• Visualization standards
• Common reporting mistakes

Reporting principles

Bayesian results are inherently richer than frequentist results -- use that richness.

1. Report full posteriors, not just point estimates. Use HDI (Highest Density Interval) as default, 94% unless context demands otherwise (94% avoids the false precision of 95% and conveys its arbitrary nature). 2. Visualize uncertainty. Every parameter estimate should have a visual representation of its posterior. 3. Show the model. Include a model specification section -- readers should know exactly what was assumed. Use pm.model_to_graphviz to visualize the model. 4. Report diagnostics. Convergence and model criticism results build trust. 5. Use probability language, not p-value language. "There is a 94% probability that θ lies in [a, b]" — not "the interval [a, b] is significant."

Analysis report template

Use this structure for written reports. Adapt sections as needed.

# [Analysis Title]

## Executive summary
[2-3 sentence summary of key findings with credible intervals]

## Data description
- Source: [where the data came from]
- Sample size: N = [n]
- Key variables: [list with descriptions]
- Notable features: [missingness, outliers, grouping structure]

## Model specification

### Generative story
[Plain-language description of the assumed data-generating process]

### Mathematical notation
[Model equations using standard notation]

### Prior choices
| Parameter | Prior | Justification |
|-----------|-------|---------------|
| β | Normal(0, 2.5) | Weakly informative on standardized predictors |
| σ | Gamma(2, 2) | Allows wide range of residual variation |

### Model graph
[Figure: model graph, probably done with `pm.model_to_graphviz`]

### Prior predictive check
[Figure: prior predictive distribution vs. plausible data range]
[Brief assessment: "Priors generate data in the range [a, b], consistent with domain knowledge."]

## Results

### Convergence diagnostics
- R-hat: all ≤ 1.01 ✓
- ESS (bulk): minimum [X] ✓
- ESS (tail): minimum [X] ✓
- Divergences: [N] [✓ or ✗]
[Figure: trace/rank plots for key parameters]

### Parameter estimates
| Parameter | Mean | SD | 94% HDI |
|-----------|------|----|---------|
| β₁ | 0.45 | 0.12 | [0.22, 0.68] |
| σ | 1.23 | 0.08 | [1.08, 1.38] |

[Figure: forest plot of key parameters]

### Model criticism
- **Posterior predictive check**: [Figure + assessment]
- **LOO-CV**: ELPD = [X] (SE = [Y]), p_loo = [Z]
- **Calibration**: [Figure + assessment]
- **Pareto k**: all < 0.5 ✓ (or list problematic observations)

### Model comparison (if applicable)
| Model | ELPD | SE | ΔELPD | Weight |
|-------|------|----|-------|--------|
[comparison table]

### Prior sensitivity
| Parameter | Prior CJS | Likelihood CJS | Diagnostic |
|-----------|-----------|----------------|------------|
| β₁ | 0.02 | 0.01 | Low sensitivity ✓ |
| σ | 0.08 | 0.03 | Strong prior / weak likelihood — justified by [domain rationale] |

[Brief interpretation: which parameters are sensitive, whether this affects conclusions,
and justification for any retained informative priors.]

## Interpretation
[What do the results mean substantively? Discuss effect sizes, practical significance,
and how uncertainty affects conclusions. Be explicit about what the model does NOT tell us.]

## Limitations
[Model assumptions that may not hold. Sensitivity to prior choices. Data limitations.]

## Appendix
[Full ArviZ summary table. Additional diagnostic plots. Code repository link.]

Presentation template

For slide-based presentations, use this structure:

Slide 1: Title + one-sentence finding
Slide 2: The question -- what are we trying to learn?
Slide 3: Data overview (1 figure, minimal text)
Slide 4: Model diagram or generative story (visual, not equations)
Slide 5: Prior predictive check ("Our assumptions produce plausible data")
Slide 6: Key results -- posterior distributions (forest plot or ridgeplot)
Slide 7: Posterior predictive check ("The data could have been credibly produced by the model")
Slide 8: Practical implications -- translate posteriors into decisions
Slide 9: Limitations and next steps

Presentation rules:

• One idea per slide
• Visualize posteriors, do not just show tables
• Use ridgeplots or forest plots for multiple parameters
• Show uncertainty in predictions (fan charts, spaghetti plots)
• For non-technical audiences: translate credible intervals into natural language ("There is a 5-in-6 chance that the effect is between X and Y")

Visualization standards

Parameter posteriors

# Forest plot (multiple parameters)
az.plot_forest(idata, var_names=["beta"], combined=True)

# Ridge plot (distributions)
az.plot_forest(idata, kind='ridgeplot', var_names=["beta"], combined=True)

# Pair plot (correlations between parameters)
az.plot_pair(idata, var_names=["beta", "sigma"])

Predictions with uncertainty

The following are very simple examples (most of the time you will use ArviZ's built-in functions, as well as work directly with the InferenceData object through xarray). Nevertheless, the following will give you an idea of the concepts we're interested in:

# Fan chart for time series predictions
percentiles = [5, 25, 50, 75, 95]
for lo, hi in [(5, 95), (25, 75)]:
    plt.fill_between(
        x,
        np.percentile(preds, lo, axis=0),
        np.percentile(preds, hi, axis=0),
        alpha=0.3,
    )
plt.plot(
    x,
    np.percentile(preds, 50, axis=0),
    color="blue",
    label="Median",
)
plt.scatter(
    x_obs,
    y_obs,
    color="black",
    s=10,
    label="Observed",
)

# Spaghetti plot (individual posterior draws)
for i in range(50):
    plt.plot(x, preds[i], alpha=0.05, color="blue")

Calibration plot

Use the new ArviZ's (1.0+) plot_ppc_pit function. Refer to this guide for guidance about coverage interpretation, but also across the whole Bayesian workflow in general -- it's a treasure trove!

Adapting for your audience

The report template above is for a technical audience. When the user mentions a non-technical audience (a boss, a medical board, stakeholders, executives), or says they're new to Bayesian stats, adapt the report:

For non-technical audiences:

• Generate a standalone markdown report file — not just code with inline comments. The report should be readable on its own, without looking at any code.
• Replace jargon with plain language: "There is roughly a 19-in-20 chance the effect is between X and Y" instead of "94% HDI: [X, Y]"
• Include a glossary defining terms like posterior, credible interval, prior, MCMC in everyday language
• Lead with the practical conclusion ("The drug lowers blood pressure by about 10 mmHg"), then support it with the statistical evidence
• Move technical details (convergence diagnostics, model equations) to an appendix
• Use the section "How the Analysis Works (Plain Language)" to explain the generative story as a narrative, not equations

For Bayesian beginners:

• Explain why each step matters, not just what it does: "We check our assumptions first (prior predictive check) to make sure our model doesn't predict impossible values like negative blood pressure"
• Contrast with frequentist approaches where helpful: "Unlike a p-value, a credible interval directly tells you where the parameter probably lies"
• Define acronyms on first use: "HDI (Highest Density Interval — the narrowest range containing 94% of the plausible values)"

Common reporting mistakes

1. Reporting only posterior means: Always include credible intervals. The uncertainty IS the result. 2. Using frequentist language: Avoid "significant", "p < 0.05", "fail to reject". Use "probability", "credible interval", "posterior probability of direction." 3. Hiding diagnostics: If convergence was imperfect, say so. If you had to fix divergences, describe how. 4. Ignoring practical significance: A posterior that excludes zero is not automatically important. Discuss effect sizes in context. 5. Not showing prior sensitivity: Run psense_summary(idata) and report the results — especially for policy-relevant or controversial conclusions. Show the CJS values table, flag any parameters with CJS > 0.05, and briefly explain whether the sensitivity affects your conclusions. If you have an intentionally informative prior that flags as sensitive, justify it explicitly rather than hiding the diagnostic. Readers should know which conclusions depend on prior choices and which are robust. 6. Skipping the generative story: The model specification should make clear what process is assumed to have generated the data.

Prior/Likelihood Sensitivity Analysis

Power-scaling sensitivity analysis checks whether your posterior conclusions are robust to reasonable changes in prior (or likelihood) strength — without refitting the model. It uses Pareto-smoothed importance sampling (PSIS) to simulate what would happen if you made your priors stronger or weaker.

Contents

• Requirements
• Running sensitivity checks
• Interpreting results
• Which variables to check
• Visual diagnostics
• Key principle

Requirements

The InferenceData object must contain both log_likelihood and log_prior groups.

Important: nutpie silently ignores idata_kwargs={"log_likelihood": True} and idata_kwargs={"log_prior": True}. You must compute these after sampling:

# After pm.sample():
pm.compute_log_likelihood(idata, model=model)
pm.compute_log_prior(idata, model=model)

Verify they're present before running sensitivity checks:

assert "log_likelihood" in idata, "Missing log_likelihood — run pm.compute_log_likelihood(idata, model=model)"
assert "log_prior" in idata, "Missing log_prior — run pm.compute_log_prior(idata, model=model)"

Running sensitivity checks

from arviz_stats import psense_summary

summary = psense_summary(idata)
summary

This returns a per-variable table with Cumulative Jensen-Shannon (CJS) divergence values for both prior and likelihood perturbations. CJS > 0.05 flags sensitivity (roughly equivalent to a 0.3 standard deviation shift in the posterior mean).

Interpreting results

Four diagnostic patterns:

Which variables to check

Not every parameter needs sensitivity analysis. Focus on what matters:

• Check: interpretable coefficients, effect sizes, predictions, derived quantities (e.g., Bayesian R², contrasts)
• Skip: group-specific parameters in hierarchical models (power-scale only the top-level hyperpriors), spline/GP basis coefficients, variance components you don't interpret directly

For hierarchical models, sensitivity of the hyperprior (e.g., the group-level standard deviation) is more informative than sensitivity of individual group effects.

# Check only specific variables
summary = psense_summary(idata, var_names=["beta", "sigma"])

Visual diagnostics

plot_psense_dist — How the posterior shifts

Shows posterior marginals at three power-scaling levels (α = 0.8, 1.0, 1.25). Use this to see the direction and magnitude of the shift, not just whether it exceeds the threshold.

from arviz_plots import plot_psense_dist

plot_psense_dist(idata, var_names=["beta"])

Requires arviz-plots (pip install arviz-plots).

plot_psense_quantities — Sensitivity of derived quantities

Shows how predictions or summary statistics shift under perturbation. Use this when you care more about predictive robustness than individual parameter sensitivity.

from arviz_plots import plot_psense_quantities

plot_psense_quantities(idata)

Key principle

Sensitivity warnings are not automatic problems. An intentionally informative prior — grounded in domain knowledge or previous studies — will legitimately flag as sensitive. That's expected: if you have a strong prior and modest data, the prior should matter.

The correct response to a sensitivity flag is: 1. Document the flag and its magnitude 2. Justify why the prior is appropriate (or acknowledge it isn't) 3. Report the sensitivity transparently — readers should know which conclusions depend on prior choices

Do not reflexively loosen priors to silence diagnostics. A well-justified informative prior that flags sensitivity is better science than a vague prior that passes all checks but encodes no knowledge.