Bayesian Statistical Analysis

This document provides guidance on conducting and interpreting Bayesian statistical analyses, which offer an alternative framework to frequentist (classical) statistics.

Contents

• Bayesian vs. Frequentist Philosophy
• Bayes' Theorem
• Prior Distributions
• Bayesian Hypothesis Testing
• Bayesian Estimation
• Common Bayesian Analyses
• Hierarchical (Multilevel) Models
• Model Comparison
• Checking Bayesian Models
• Reporting Bayesian Results
• Advantages and Limitations
• Key Python Packages
• When to Use Bayesian Methods

Bayesian vs. Frequentist Philosophy

Fundamental Differences

Key Question Difference

Frequentist: "If the null hypothesis is true, what is the probability of observing data this extreme or more extreme?"

Bayesian: "Given the observed data, what is the probability that the hypothesis is true?"

The Bayesian question is more intuitive and directly addresses what researchers want to know.

Bayes' Theorem

Formula:

P(θ|D) = P(D|θ) × P(θ) / P(D)

In words:

Posterior = Likelihood × Prior / Evidence

Where:

• θ (theta): Parameter of interest (e.g., mean difference, correlation)
• D: Observed data
• P(θ|D): Posterior distribution (belief about θ after seeing data)
• P(D|θ): Likelihood (probability of data given θ)
• P(θ): Prior distribution (belief about θ before seeing data)
• P(D): Marginal likelihood/evidence (normalizing constant)

Prior Distributions

Types of Priors

1. Informative Priors

When to use: When you have substantial prior knowledge from:

• Previous studies
• Expert knowledge
• Theory
• Pilot data

Example: Meta-analysis shows effect size d ≈ 0.40, SD = 0.15

• Prior: Normal(0.40, 0.15)

Advantages:

• Incorporates existing knowledge
• More efficient (smaller samples needed)
• Can stabilize estimates with small data

Disadvantages:

• Subjective (but subjectivity can be strength)
• Must be justified and transparent
• May be controversial if strong prior conflicts with data

2. Weakly Informative Priors

When to use: Default choice for most applications

Characteristics:

• Regularizes estimates (prevents extreme values)
• Has minimal influence on posterior with moderate data
• Prevents computational issues

Example priors:

• Effect size: Normal(0, 1) or Cauchy(0, 0.707)
• Variance: Half-Cauchy(0, 1)
• Correlation: Uniform(-1, 1), a rescaled Beta on (-1, 1) (e.g. 2×Beta(2, 2)−1; plain Beta(2, 2) has support [0, 1]), or an LKJ prior for correlation matrices

Advantages:

• Balances objectivity and regularization
• Computationally stable
• Broadly acceptable

3. Non-Informative (Flat/Uniform) Priors

When to use: When attempting to be "objective"

Example: Uniform(-∞, ∞) for any value

⚠️ Caution:

• Can lead to improper posteriors
• May produce non-sensible results
• Not truly "non-informative" (still makes assumptions)
• Often not recommended in modern Bayesian practice

Better alternative: Use weakly informative priors

Prior Sensitivity Analysis

Always conduct: Test how results change with different priors

Process:

1. Fit model with default/planned prior
2. Fit model with more diffuse prior
3. Fit model with more concentrated prior
4. Compare posterior distributions

Reporting:

• If results are similar: Evidence is robust
• If results differ substantially: Data are not strong enough to overwhelm prior

Python example:

import pymc as pm

prior_specs = [
    ('weakly_informative', 0, 1),
    ('diffuse', 0, 10),
    ('informative', 0.5, 0.3),
]

results = {}
for name, mu_prior, sigma_prior in prior_specs:
    with pm.Model() as model:
        effect = pm.Normal('effect', mu=mu_prior, sigma=sigma_prior)
        # ... likelihood and observed data
        trace = pm.sample(2000, tune=1000)
        results[name] = trace

Bayesian Hypothesis Testing

Bayes Factor (BF)

What it is: Ratio of evidence for two competing hypotheses

Formula:

BF₁₀ = P(D|H₁) / P(D|H₀)

Interpretation:

Advantages over p-values:

1. Can provide evidence for null hypothesis
2. Less dependent on sampling intentions than p-values — but only Bayes factors with fixed priors are relatively insensitive to optional stopping; posterior-based decision rules are still affected, and transparency requires reporting the stopping rule
3. Directly quantifies evidence
4. Can be updated with more data

Python calculation:

# Pingouin 0.5+: BF10 for independent two-sided t-tests; one-sided BF removed.
import pingouin as pg

result = pg.ttest(group1, group2, correction=False)
bf10 = result['BF10'].values[0]

# Rigorous Bayes Factors: BayesFactor (R), JASP, or PyMC model comparison (see pymc skill)

Region of Practical Equivalence (ROPE)

Purpose: Define range of negligible effect sizes

Process:

1. Define ROPE (e.g., d ∈ [-0.1, 0.1] for negligible effects)
2. Calculate % of posterior inside ROPE
3. Make decision:

   • 95% in ROPE: Accept practical equivalence

   • 95% outside ROPE: Reject equivalence

   • Otherwise: Inconclusive

Advantage: Directly tests for practical significance

Python example:

# Define ROPE
rope_lower, rope_upper = -0.1, 0.1

# Calculate % of posterior in ROPE
in_rope = np.mean((posterior_samples > rope_lower) &
                  (posterior_samples < rope_upper))

print(f"{in_rope*100:.1f}% of posterior in ROPE")

Bayesian Estimation

Credible Intervals

What it is: Interval containing parameter with X% probability

95% Credible Interval interpretation:

"There is a 95% probability that the true parameter lies in this interval."

This is what people THINK confidence intervals mean (but don't in frequentist framework)

Types:

Equal-Tailed Interval (ETI)

• 2.5th to 97.5th percentile
• Simple to calculate
• May not include mode for skewed distributions

Highest Density Interval (HDI)

• Narrowest interval containing 95% of distribution
• Always includes mode
• Better for skewed distributions

Python calculation:

import arviz as az

# Equal-tailed interval
eti = np.percentile(posterior_samples, [2.5, 97.5])

# HDI (ArviZ 1.x renamed the keyword hdi_prob= to prob=)
hdi = az.hdi(posterior_samples, prob=0.95)

Posterior Distributions

Interpreting posterior distributions:

1. Central tendency:

   • Mean: Average posterior value
   • Median: 50th percentile
   • Mode: Most probable value (MAP - Maximum A Posteriori)

2. Uncertainty:

   • SD: Spread of posterior
   • Credible intervals: Quantify uncertainty

3. Shape:

   • Symmetric: Similar to normal
   • Skewed: Asymmetric uncertainty
   • Multimodal: Multiple plausible values

Visualization:

import matplotlib.pyplot as plt
import arviz as az

# Posterior plot with 95% credible interval
# (ArviZ 1.x replaced plot_posterior with plot_dist and hdi_prob= with ci_prob=)
az.plot_dist(trace, ci_prob=0.95)

# Trace plot (check convergence)
az.plot_trace(trace)

# Forest plot (multiple parameters)
az.plot_forest(trace)

Common Bayesian Analyses

Bayesian T-Test

Purpose: Compare two groups (Bayesian alternative to t-test)

Outputs:

1. Posterior distribution of mean difference
2. 95% credible interval
3. Bayes Factor (BF₁₀)
4. Probability of directional hypothesis (e.g., P(μ₁ > μ₂))

Python implementation:

import pymc as pm
import arviz as az

# Bayesian independent samples t-test
with pm.Model() as model:
    # Priors for group means
    mu1 = pm.Normal('mu1', mu=0, sigma=10)
    mu2 = pm.Normal('mu2', mu=0, sigma=10)

    # Prior for pooled standard deviation
    sigma = pm.HalfNormal('sigma', sigma=10)

    # Likelihood
    y1 = pm.Normal('y1', mu=mu1, sigma=sigma, observed=group1)
    y2 = pm.Normal('y2', mu=mu2, sigma=sigma, observed=group2)

    # Derived quantity: mean difference
    diff = pm.Deterministic('diff', mu1 - mu2)

    # Sample posterior
    trace = pm.sample(2000, tune=1000)

# Analyze results
print(az.summary(trace, var_names=['mu1', 'mu2', 'diff']))

# Probability that group1 > group2
prob_greater = np.mean(trace.posterior['diff'].values > 0)
print(f"P(μ₁ > μ₂) = {prob_greater:.3f}")

# Plot posterior (ArviZ 1.x: plot_posterior was replaced by plot_dist;
# add a reference line at 0 with matplotlib if needed)
az.plot_dist(trace, var_names=['diff'])

Bayesian ANOVA

Purpose: Compare three or more groups

Model:

import pymc as pm

with pm.Model() as anova_model:
    # Hyperpriors
    mu_global = pm.Normal('mu_global', mu=0, sigma=10)
    sigma_between = pm.HalfNormal('sigma_between', sigma=5)
    sigma_within = pm.HalfNormal('sigma_within', sigma=5)

    # Group means (hierarchical)
    group_means = pm.Normal('group_means',
                            mu=mu_global,
                            sigma=sigma_between,
                            shape=n_groups)

    # Likelihood
    y = pm.Normal('y',
                  mu=group_means[group_idx],
                  sigma=sigma_within,
                  observed=data)

    trace = pm.sample(2000, tune=1000)

# Posterior contrasts
contrast_1_2 = trace.posterior['group_means'][:,:,0] - trace.posterior['group_means'][:,:,1]

Bayesian Correlation

Purpose: Estimate correlation between two variables

Advantage: Provides distribution of correlation values

Python implementation:

import numpy as np
import pymc as pm

# Standardize both variables first: with z-scored data the bivariate normal
# can fix mu = [0, 0] and unit variances, leaving rho as the only free
# parameter (correlation is unchanged by linear rescaling). Alternatively,
# model the means and SDs as parameters (or use pm.LKJCholeskyCov).
xz = (x - x.mean()) / x.std()
yz = (y - y.mean()) / y.std()

with pm.Model() as corr_model:
    # Prior on correlation
    rho = pm.Uniform('rho', lower=-1, upper=1)

    # Correlation (= covariance) matrix for standardized data
    cov_matrix = pm.math.stack([[1, rho],
                                [rho, 1]])

    # Likelihood (bivariate normal on standardized data)
    obs = pm.MvNormal('obs',
                     mu=[0, 0],
                     cov=cov_matrix,
                     observed=np.column_stack([xz, yz]))

    trace = pm.sample(2000, tune=1000)

# Summarize correlation
print(az.summary(trace, var_names=['rho']))

# Probability that correlation is positive
prob_positive = np.mean(trace.posterior['rho'].values > 0)

Bayesian Linear Regression

Purpose: Model relationship between predictors and outcome

Advantages:

• Uncertainty in all parameters
• Natural regularization (via priors)
• Can incorporate prior knowledge
• Credible intervals for predictions

Python implementation:

import pymc as pm

with pm.Model() as regression_model:
    # Mutable data container: required for pm.set_data() to swap in new
    # predictors later (a raw array would make set_data fail)
    X_data = pm.Data('X', X)

    # Priors for coefficients
    alpha = pm.Normal('alpha', mu=0, sigma=10)  # Intercept
    beta = pm.Normal('beta', mu=0, sigma=10, shape=n_predictors)
    sigma = pm.HalfNormal('sigma', sigma=10)

    # Expected value
    mu = alpha + pm.math.dot(X_data, beta)

    # Likelihood (shape=mu.shape lets predictions resize with new data)
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y, shape=mu.shape)

    trace = pm.sample(2000, tune=1000)

# Posterior predictive checks
with regression_model:
    ppc = pm.sample_posterior_predictive(trace)

az.plot_ppc_dist(ppc)  # ArviZ 1.x: plot_ppc was replaced by plot_ppc_dist

# Predictions with uncertainty
with regression_model:
    pm.set_data({'X': X_new})
    posterior_pred = pm.sample_posterior_predictive(trace, predictions=True)

Hierarchical (Multilevel) Models

When to use:

• Nested/clustered data (students within schools)
• Repeated measures
• Meta-analysis
• Varying effects across groups

Key concept: Partial pooling

• Complete pooling: Ignore groups (biased)
• No pooling: Analyze groups separately (high variance)
• Partial pooling: Borrow strength across groups (Bayesian)

Example: Varying intercepts:

with pm.Model() as hierarchical_model:
    # Hyperpriors
    mu_global = pm.Normal('mu_global', mu=0, sigma=10)
    sigma_between = pm.HalfNormal('sigma_between', sigma=5)
    sigma_within = pm.HalfNormal('sigma_within', sigma=5)

    # Group-level intercepts
    alpha = pm.Normal('alpha',
                     mu=mu_global,
                     sigma=sigma_between,
                     shape=n_groups)

    # Likelihood
    y_obs = pm.Normal('y_obs',
                     mu=alpha[group_idx],
                     sigma=sigma_within,
                     observed=y)

    trace = pm.sample()

Model Comparison

Methods

1. Bayes Factor

• Directly compares model evidence
• Sensitive to prior specification
• Can be computationally intensive

2. Information Criteria

WAIC (Widely Applicable Information Criterion):

• Bayesian analog of AIC
• Reported on the elpd (expected log pointwise predictive density) scale: HIGHER elpd is better (only on the deviance scale, −2 × elpd, is lower better)
• Accounts for effective number of parameters
• az.waic was removed in ArviZ 1.x — use LOO

LOO (Leave-One-Out Cross-Validation):

• Estimates out-of-sample prediction error
• Also on the elpd scale: higher elpd is better
• More robust than WAIC

Python calculation:

import arviz as az
import pymc as pm

# LOO needs pointwise log-likelihoods
with model:
    pm.compute_log_likelihood(trace)

loo = az.loo(trace)
print(f"LOO elpd: {loo.elpd:.2f}")  # higher is better

# Compare multiple models: az.compare ranks them correctly (rank 0 = best)
comparison = az.compare({
    'model1': trace1,
    'model2': trace2,
    'model3': trace3
})
print(comparison)

Checking Bayesian Models

1. Convergence Diagnostics

R-hat (Gelman-Rubin statistic):

• Compares within-chain and between-chain variance
• Values close to 1.0 indicate convergence
• R-hat < 1.01: Good
• R-hat > 1.05: Poor convergence

Effective Sample Size (ESS):

• Number of independent samples
• Higher is better
• Bulk-ESS > 400 in total across all chains recommended (Vehtari et al., 2021); also check tail-ESS

Trace plots:

• Should look like "fuzzy caterpillar"
• No trends, no stuck chains

Python checking:

# Automatic summary with diagnostics
print(az.summary(trace, var_names=['parameter']))

# Visual diagnostics
az.plot_trace(trace)
az.plot_rank(trace)  # Rank plots

2. Posterior Predictive Checks

Purpose: Does model generate data similar to observed data?

Process:

1. Generate predictions from posterior
2. Compare to actual data
3. Look for systematic discrepancies

Python implementation:

with model:
    ppc = pm.sample_posterior_predictive(trace)

# Visual check (ArviZ 1.x: plot_ppc was replaced by plot_ppc_dist)
az.plot_ppc_dist(ppc, num_samples=100)

# Quantitative check: compute the statistic per posterior draw over the
# observation dimension (do NOT iterate the array directly - that loops
# over chains, not draws)
obs_mean = np.mean(observed_data)
pp = ppc.posterior_predictive['y_obs']                # dims: (chain, draw, obs)
pred_means = pp.mean(dim=pp.dims[-1]).values.ravel()  # one mean per draw
p_value = np.mean(pred_means >= obs_mean)  # Bayesian (posterior predictive) p-value

Reporting Bayesian Results

Example T-Test Report

"A Bayesian independent samples t-test was conducted to compare groups A and B. Weakly informative priors were used: Normal(0, 1) for the mean difference and Half-Cauchy(0, 1) for the pooled standard deviation. The posterior distribution of the mean difference had a mean of 5.2 (95% CI [2.3, 8.1]), indicating that Group A scored higher than Group B. The Bayes Factor BF₁₀ = 23.5 provided strong evidence for a difference between groups, and there was a 99.7% probability that Group A's mean exceeded Group B's mean."

Example Regression Report

"A Bayesian linear regression was fitted with weakly informative priors (Normal(0, 10) for coefficients, Half-Cauchy(0, 5) for residual SD). The model explained substantial variance (R² = 0.47, 95% CI [0.38, 0.55]). Study hours (β = 0.52, 95% CI [0.38, 0.66]) and prior GPA (β = 0.31, 95% CI [0.17, 0.45]) were credible predictors (95% CIs excluded zero). Posterior predictive checks showed good model fit. Convergence diagnostics were satisfactory (all R-hat < 1.01, ESS > 1000)."

Advantages and Limitations

Advantages

1. Intuitive interpretation: Direct probability statements about parameters
2. Incorporates prior knowledge: Uses all available information
3. Flexible: Handles complex models easily
4. Less sensitive to optional stopping: Bayes factors with fixed priors are relatively robust to analyzing data as it arrives, but posterior-based decision rules are still affected — always report the stopping rule
5. Quantifies uncertainty: Full posterior distribution
6. Small samples: Works at any sample size (but small-n posteriors are prior-dominated — report a prior-sensitivity check)

Limitations

1. Computational: Requires MCMC sampling (can be slow)
2. Prior specification: Requires thought and justification
3. Complexity: Steeper learning curve
4. Software: Fewer tools than frequentist methods
5. Communication: May need to educate reviewers/readers

Key Python Packages

Install with uv (see SKILL.md). ArviZ requires Python >= 3.10. ArviZ 1.x is the current line and is a breaking rewrite: az.summary defaults to 89% intervals and takes ci_prob= (the old hdi_prob= keyword is gone), az.hdi takes prob=, az.plot_posterior/az.plot_ppc were replaced by az.plot_dist/az.plot_ppc_dist, and az.waic was removed (use az.loo).

• PyMC (pymc>=5): Full Bayesian modeling framework
• ArviZ (arviz>=1.0): Visualization and diagnostics (docs)
• Bambi: High-level interface for regression models (uv pip install bambi)
• cmdstanpy: Python interface to Stan (use instead of the discontinued PyStan)
• TensorFlow Probability: Bayesian inference with TensorFlow

When to Use Bayesian Methods

Use Bayesian when:

• You have prior information to incorporate
• You want direct probability statements
• Sample size is small
• Model is complex (hierarchical, missing data, etc.)
• You want to update analysis as data arrives

Frequentist may be sufficient when:

• Standard analysis with large sample
• No prior information
• Computational resources limited
• Reviewers unfamiliar with Bayesian methods