# 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](#bayesian-vs-frequentist-philosophy)
- [Bayes' Theorem](#bayes-theorem)
- [Prior Distributions](#prior-distributions)
- [Bayesian Hypothesis Testing](#bayesian-hypothesis-testing)
- [Bayesian Estimation](#bayesian-estimation)
- [Common Bayesian Analyses](#common-bayesian-analyses)
- [Hierarchical (Multilevel) Models](#hierarchical-multilevel-models)
- [Model Comparison](#model-comparison)
- [Checking Bayesian Models](#checking-bayesian-models)
- [Reporting Bayesian Results](#reporting-bayesian-results)
- [Advantages and Limitations](#advantages-and-limitations)
- [Key Python Packages](#key-python-packages)
- [When to Use Bayesian Methods](#when-to-use-bayesian-methods)

## Bayesian vs. Frequentist Philosophy

### Fundamental Differences

| Aspect | Frequentist | Bayesian |
|--------|-------------|----------|
| **Probability interpretation** | Long-run frequency of events | Degree of belief/uncertainty |
| **Parameters** | Fixed but unknown | Random variables with distributions |
| **Inference** | Based on sampling distributions | Based on posterior distributions |
| **Primary output** | p-values, confidence intervals | Posterior probabilities, credible intervals |
| **Prior information** | Not formally incorporated | Explicitly incorporated via priors |
| **Hypothesis testing** | Reject/fail to reject null | Probability of hypotheses given data |
| **Sample size** | Often requires minimum | Works at any n, but small-n posteriors are prior-dominated — report a prior-sensitivity check |
| **Interpretation** | Indirect (probability of data given H₀) | Direct (probability of hypothesis given data) |

### 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**:
```python
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**:

| BF₁₀ | Evidence |
|------|----------|
| >100 | Decisive for H₁ |
| 30-100 | Very strong for H₁ |
| 10-30 | Strong for H₁ |
| 3-10 | Moderate for H₁ |
| 1-3 | Anecdotal for H₁ |
| 1 | No evidence |
| 1/3-1 | Anecdotal for H₀ |
| 1/10-1/3 | Moderate for H₀ |
| 1/30-1/10 | Strong for H₀ |
| 1/100-1/30 | Very strong for H₀ |
| <1/100 | Decisive for H₀ |

**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**:
```python
# 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**:
```python
# 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**:
```python
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**:
```python
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**:
```python
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**:
```python
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**:
```python
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**:
```python
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**:
```python
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**:
```python
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**:
```python
# 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**:
```python
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](https://python.arviz.org))
- **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