statistics-advanced

Advanced Statistics

Overview

This skill extends the basic statistical-analysis skill with advanced methods for complex data structures: hierarchical/multilevel data, Bayesian inference, resampling-based inference, spatial data, and missing data. Use the basic skill first for standard hypothesis tests.

When to Use This Skill

• Data has hierarchical or nested structure (students in schools, repeated measures per subject)
• You want a full posterior distribution, not just a point estimate and p-value
• Small sample sizes where asymptotic normality doesn't hold
• Data has spatial autocorrelation
• Substantial missing data that cannot be ignored

1. Bayesian Inference with PyMC

When to Use Bayesian Analysis

• Small n (< 30 per group) where priors help regularize
• You have genuine prior knowledge about parameter ranges
• You need full uncertainty quantification (not just confidence intervals)
• Hierarchical / multilevel models are needed

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

# Example: Bayesian t-test
np.random.seed(42)
control = np.random.normal(10, 2, 20)
treatment = np.random.normal(12, 2.5, 20)

with pm.Model() as model:
    # Priors (weakly informative)
    mu_ctrl = pm.Normal("mu_ctrl", mu=10, sigma=5)
    mu_treat = pm.Normal("mu_treat", mu=10, sigma=5)
    sigma_ctrl = pm.HalfNormal("sigma_ctrl", sigma=3)
    sigma_treat = pm.HalfNormal("sigma_treat", sigma=3)

    # Effect size (Cohen's d)
    diff = pm.Deterministic("difference", mu_treat - mu_ctrl)
    pooled_sigma = pm.Deterministic("pooled_sigma",
                                    pm.math.sqrt((sigma_ctrl**2 + sigma_treat**2) / 2))
    effect_size = pm.Deterministic("effect_size", diff / pooled_sigma)

    # Likelihood
    obs_ctrl = pm.Normal("obs_ctrl", mu=mu_ctrl, sigma=sigma_ctrl, observed=control)
    obs_treat = pm.Normal("obs_treat", mu=mu_treat, sigma=sigma_treat, observed=treatment)

    # Sample
    trace = pm.sample(2000, chains=4, target_accept=0.9, random_seed=42)

# Diagnostics
summary = az.summary(trace, var_names=["difference", "effect_size"])
print(summary)

# Key diagnostics to check:
# r_hat < 1.01 → chains converged
# ess_bulk > 400 → enough effective samples
az.plot_trace(trace, var_names=["difference"])

# Posterior predictive check
with model:
    ppc = pm.sample_posterior_predictive(trace)