Bayesian Statistical Analysis

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

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) or Beta(2, 2)

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:

# Model with different priors
priors = [
    ('weakly_informative', pm.Normal.dist(0, 1)),
    ('diffuse', pm.Normal.dist(0, 10)),
    ('informative', pm.Normal.dist(0.5, 0.3))
]

results = {}
for name, prior in priors:
    with pm.Model():
        effect = pm.Normal('effect', mu=prior.mu, sigma=prior.sigma)
        # ... rest of model
        trace = pm.sample()
        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. Not dependent on sampling intentions (no "peeking" problem)
3. Directly quantifies evidence
4. Can be updated with more data

Python calculation:

# Note: Limited BF support in Python
# Better options: R packages (BayesFactor), JASP software

# Approximate BF from t-statistic
# Using Jeffreys-Zellner-Siow prior
from scipy import stats

def bf_from_t(t, n1, n2, r_scale=0.707):
    """
    Approximate Bayes Factor from t-statistic
    r_scale: Cauchy prior scale (default 0.707 for medium effect)
    """
    # This is simplified; use dedicated packages for accurate calculation
    df = n1 + n2 - 2
    # Implementation requires numerical integration
    pass

---

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:

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

# HDI
hdi = az.hdi(posterior_samples, hdi_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:

# Posterior plot with HDI
az.plot_posterior(trace, hdi_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:

# 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, return_inferencedata=True)

# 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
az.plot_posterior(trace, var_names=['diff'], ref_val=0)

---

Bayesian ANOVA

Purpose: Compare three or more groups

Model:

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, return_inferencedata=True)

# 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:

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

    # Convert to covariance matrix
    cov_matrix = pm.math.stack([[1, rho],
                                [rho, 1]])

    # Likelihood (bivariate normal)
    obs = pm.MvNormal('obs',
                     mu=[0, 0],
                     cov=cov_matrix,
                     observed=np.column_stack([x, y]))

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

# 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:

with pm.Model() as regression_model:
    # 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, beta)

    # Likelihood
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y)

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

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

az.plot_ppc(ppc)

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

---

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
• Lower is better
• Accounts for effective number of parameters

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

• Estimates out-of-sample prediction error
• Lower is better
• More robust than WAIC

Python calculation:

# Calculate WAIC and LOO
waic = az.waic(trace)
loo = az.loo(trace)

print(f"WAIC: {waic.elpd_waic:.2f}")
print(f"LOO: {loo.elpd_loo:.2f}")

# Compare multiple models
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
• ESS > 400 per chain recommended

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
az.plot_ppc(ppc, num_pp_samples=100)

# Quantitative checks
obs_mean = np.mean(observed_data)
pred_means = [np.mean(sample) for sample in ppc.posterior_predictive['y_obs']]
p_value = np.mean(pred_means >= obs_mean)  # Bayesian p-value

---

Reporting Bayesian Results

Example T-Test Report

"A Bayesian inde