Statistical Modeling

Statistical Modeling for Sports

Overview

This skill's job is classical statistical inference on a user-owned sports table: choose a GLM or test, check assumptions, report effects and uncertainty.

It is not a universal modeling workflow. Do not use it to load data, design walk-forward ML horse races, or certify leakage. Mixed, Bayesian, causal, and power-analysis material below is routed methodology — the bundled scripts only cover univariate assumption checks and an additive numeric binomial GLM.

Conduct the analysis the way a careful analyst would: right model family for the outcome, verified assumptions, honest effect sizes and uncertainty, and a write-up that survives review.

Typical stack: pandas, numpy, scipy, statsmodels, and scikit-learn (+ optional pingouin, seaborn, pymc, and arviz).

When to Use This Skill

Use this skill when:

• Modeling win/loss, margins, points, goals, runs, or other sports outcomes
• Comparing groups (home vs away, pre/post rule change, position groups)
• Fitting logistic, linear, Poisson, or NegBin models on sports panels
• Checking residual diagnostics and reporting effects with uncertainty
• After the GLM path, clustered or partially pooled effects (see hierarchical section)
• Reporting coefficients, odds ratios, effect sizes, and model fit
• Choosing between classical stats models and jumping to ML

Do not use this skill as a substitute for:

• raw data loading → use the appropriate loader skill, then return with a user-owned analysis table
• pure EDA → eda-sports
• walk-forward ML horse races → predictive-modeling
• leakage review of feature pipelines → leakage-audit
• rating systems as the primary object → ratings-strength-models

Installation

Install only the libraries required by the chosen analysis:

bashCopy code
python -m pip install pandas numpy scipy statsmodels scikit-learn
# Optional richer classical statistics and plots:
python -m pip install "pingouin>=0.6" "seaborn>=0.13"
# Optional Bayesian hierarchical work:
python -m pip install pymc arviz

Parquet input to the bundled helpers also needs pyarrow or fastparquet.

Compatibility notes:

• Prefer statsmodels>=0.14 with a recent SciPy.
• Pingouin 0.6+ uses p_val, cohen_d, and CI95.
• Read bayesian_statistics.md before selecting Bayesian priors or diagnostics.

Analysis Workflow

Every sound sports statistical analysis follows this arc. Do not skip.

1. Frame the question before fitting.

   • Outcome (win, margin, goals, player stat)
   • Predictors known at decision time T (if predictive)
   • Unit of analysis (game, team-game, player-game, possession)
   • Estimand: description, prediction, or causal claim (most sports work is predictive/descriptive)

2. Inspect the data (eda-sports).

   • n by season/group
   • missingness
   • score distributions, zero inflation, outliers
   • home/away balance

3. Choose the model family using the selection guide below (and references/test_selection_guide.md + references/sports_glm_guide.md).

4. Enforce time safety for predictive work.

   • Features must be knowable at T
   • Use walk-forward validation (validation-design)

5. Fit a minimal model first.

   • home indicator + one strength feature beats a kitchen sink

6. Check assumptions / diagnostics.

   • Use scripts/assumption_checks.py and scripts/glm_diagnostics.py
   • residual plots, overdispersion, leverage, calibration for probabilities

7. Compare to baselines.

   • constant rate, home-only, simple rating differential (baseline-models)

8. Report completely.

   • family, formula, n, metrics, coefficients with uncertainty, limits
   • templates in references/reporting_standards.md + sports templates below

If the user only needs one step (e.g. “is home advantage real in this sample?”), jump to that section — but still state design assumptions.

Reference routing:

• Read test_selection_guide.md when choosing a comparison or test.
• Read sports_glm_guide.md when mapping an outcome to a family and link.
• Read assumptions_and_diagnostics.md and diagnostics_checklist.md when validating a fitted model.
• Read effect_sizes_and_power.md for effect-size, interval, and power choices.
• Read bayesian_statistics.md for prior, partial-pooling, and posterior-diagnostic decisions.
• Read reporting_standards.md when producing the final analysis.

Model Selection Guide (Sports)

Binary outcomes (win/loss, threshold events)

Continuous outcomes (margin, rating residual)

Counts (goals, runs, sacks, threes made)

Ordered outcomes (win/draw/loss)

• Proportional odds / ordered logit
• Or two-stage: draw model + conditional direction

Hierarchical structure

• Team strength as random intercepts
• Pitcher/batter, goalie, QB effects with partial pooling
• Prefer partial pooling over dummy-every-player with tiny n

Group comparisons (not full predictive models)

Worked Path on a User-Provided Team-Game Table

Start with one row per evaluated decision and document whether that means one row per game or one row per team-game. A useful binary-outcome frame contains:

For symmetric team-game data, either model both perspectives with dependence-aware uncertainty or select one perspective per game. Never treat doubled rows as independent evidence.

Logistic model with robust inference

pythonCopy code
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf

train = model_df[model_df["season"] < 2024].copy()
fit = smf.glm(
    "won ~ is_home + rating_diff + rest_diff",
    data=train,
    family=sm.families.Binomial(),
).fit(cov_type="HC3")

odds_ratios = np.exp(fit.params)
odds_ratio_ci = np.exp(fit.conf_int())
print(fit.summary())
print(odds_ratios)
print(odds_ratio_ci)

HC3 is useful for heteroskedasticity but does not solve team, player, or temporal dependence. Use cluster-robust uncertainty or a hierarchical model when the sampling structure requires it.

Margin model

pythonCopy code
home = train[train["is_home"] == 1].copy()  # one row per game
margin_fit = smf.ols(
    "point_margin ~ rating_diff + rest_diff",
    data=home,
).fit(cov_type="HC3")
print(margin_fit.summary())

Inspect residual shape, changing variance, influential games, and nonlinear patterns. If the scientific question concerns an average contrast rather than a conditional linear effect, report that estimand directly too.

Count model

pythonCopy code
poisson_fit = smf.glm(
    "points_for ~ is_home + rating_diff",
    data=train,
    family=sm.families.Poisson(),
).fit()
dispersion = poisson_fit.pearson_chi2 / poisson_fit.df_resid
print(poisson_fit.summary())
print("Pearson dispersion:", dispersion)

Material overdispersion, excess zeros, or correlated paired scores can require a negative-binomial, hurdle/zero-inflated, or joint score model. Select that extension from diagnostics and the data-generating process—not a better-looking p-value.

Predictive evaluation

Freeze expanding time folds before fitting. Within each fold, estimate every transform and model on earlier rows only; score the later block against a training-rate constant and a simple legal-feature baseline. Report every fold, not only the mean.

For a standalone additive binomial diagnostic report:

bashCopy code
python /path/to/statistical-modeling/scripts/glm_diagnostics.py \
  --input model_frame.csv \
  --formula "won ~ is_home + rating_diff + rest_diff" \
  --out glm_report.json

The helper reports HC3 inference, odds ratios, intervals, AIC, and in-sample calibration bins. Those bins diagnose fit; they are not held-out evidence.

Assumption Checks & Diagnostics

Always check assumptions before interpreting, and report the checks.

Standalone assumption checks

bashCopy code
python /path/to/statistical-modeling/scripts/assumption_checks.py \
  --input observations.csv --value-col margin \
  --group-col venue --out assumptions.json

The helper reports Shapiro-Wilk normality, IQR outliers, and optional median-centered Levene variance testing. Formal tests supplement plots and domain judgment; at large sample sizes they can reject practically harmless departures.

Standalone binomial GLM diagnostics

bashCopy code
python /path/to/statistical-modeling/scripts/glm_diagnostics.py \
  --input model_frame.csv \
  --formula "won ~ is_home + rating_diff + rest_diff" \
  --out glm_report.json

Inputs may be CSV, Parquet, JSON, JSONL, or NDJSON. Formulas support additive numeric predictors joined by +; use --outcome-col when the binary target is not won.

What to check by family

Logistic / classification

• separation / extreme coefficients
• probability calibration (calibration-check)
• influential games (high leverage matchups)
• stability across seasons

Gaussian margin models

• residual vs fitted curvature
• QQ plot of residuals
• residual variance by season or predicted mean
• inspect influence and leverage with the fitted library's diagnostics (for example, statsmodels OLSInfluence), and use Breusch-Pagan or another predeclared variance check when appropriate

The bundled glm_diagnostics.py fits only additive numeric binomial GLMs. It does not diagnose OLS, Poisson, negative-binomial, mixed, or Bayesian models; use the model library's diagnostics and the linked checklist for those families.

Poisson scoring models

• overdispersion: Pearson χ² / df >> 1 → NegBin
• excess zeros

All predictive sports models

• leakage: any feature using same-game outcomes?
• temporal drift: does 2019 model die in 2024?

See references/assumptions_and_diagnostics.md and references/diagnostics_checklist.md.

Effect Sizes & Metrics (Sports Defaults)

p-values are usually the wrong headline for predictive sports work. Prefer:

Odds ratio from logistic coefficient:

pythonCopy code
import numpy as np
or_home = np.exp(fit.params["is_home"])
ci = np.exp(fit.conf_int().loc["is_home"])
print(or_home, ci.tolist())

For classical effect size tables, power analysis, and CI methods, use references/effect_sizes_and_power.md and pingouin:

pythonCopy code
import pingouin as pg
# home vs away margins on home-row frame already encodes home perspective;
# example two-sample on a continuous player metric:
# pg.ttest(group_a, group_b, correction='auto')  # returns cohen_d, CI95, etc.

Hierarchical / Partial Pooling Patterns

Read this section when the estimand is a clustered or partially pooled effect. This is not the default path of the skill, and there is no bundled mixed-model helper. Prefer the GLM/test workflow above unless grouping structure is the question.

When player or team n is uneven, complete-pooling (ignore groups) and no-pooling (dummy each group) both fail. Partial pooling is the default serious approach.

pythonCopy code
import statsmodels.formula.api as smf
# Mixed LM example for continuous margin; requires sufficient group size
md = smf.mixedlm("point_diff ~ is_home", data=train, groups=train["team"])
mdf = md.fit()
print(mdf.summary())

For full Bayesian hierarchical rating models, keep priors explicit and report posterior predictive checks — see references/bayesian_statistics.md.

Group Comparison Examples (Sports)

Home advantage on margins (home rows)

pythonCopy code
home = panel.loc[panel["is_home"].eq(1), ["season", "point_diff"]].dropna()
print(home["point_diff"].describe())  # pooled descriptive game-weighted mean
print(home.groupby("season")["point_diff"].agg(["count", "mean", "std"]))

Do not attach an ordinary one-sample t-test to clustered team/season data by default. Choose the uncertainty method from the sampling structure: a season-block bootstrap for a game-weighted mean, an equal-season estimand on season summaries when seasons are the independent units, or a model with defensible cluster/hierarchical structure. Games connect two teams, so a one-way team cluster can also be inadequate. State the estimand and number of independent clusters; with few clusters, use small-sample-aware methods and weaken the claim.

Correlation: pre-game form vs margin

pythonCopy code
required = {"pre_form_diff", "point_diff"}
missing = required.difference(panel.columns)
if missing:
    raise ValueError(f"build and join legal pre-game form columns first: {sorted(missing)}")
feat = panel.dropna(subset=["pre_form_diff", "point_diff"])
print(pg.corr(feat["pre_form_diff"], feat["point_diff"], method="pearson"))
print(pg.corr(feat["pre_form_diff"], feat["point_diff"], method="spearman"))

The example assumes pre_form_diff was created outside this skill with an as-of/shifted feature workflow such as time-series-sports; no bundled add_pregame_form_features function exists.

Always state that observational sports associations are not automatically causal.

Reporting Templates (Sports)

Predictive logistic win model

textCopy code
Question: Pre-game P(team win) on NFL team-game panel, seasons S0–S1.
Unit: team-game. Decision time: scheduled kickoff.
Model: logistic GLM with HC3 SE.
Formula: won ~ is_home + rating_diff + pre_form_diff
Features: all declared pre-event predictors, with construction and as-of timing documented.
Validation: walk-forward by season; train on past seasons only.
Baselines: constant train win rate; home-only logistic optional.
Results: mean walk-forward log-loss = X.XXX vs constant Y.YYY; accuracy Z.ZZ.
Coefficients (train window): is_home OR=… (95% CI …); …
Diagnostics: …; calibration: …
Limits: no opponent-adjusted EPA; roster/injury not modeled; team renames handled by raw abbreviations.

Explanatory home-advantage estimate

textCopy code
Question: Mean home margin in NFL seasons …
Design: home rows only (one row per game).
Analysis: game-weighted mean point_diff; season-block bootstrap (or another
dependence-aware method justified from the design); 95% interval.
Assumptions: one row per game; cluster definition and number of clusters stated.
Result: M=…, SD=…, dependence-aware 95% CI ….
Limits: era/rule confounds; not a causal estimate of venue alone.

APA-style general templates: references/reporting_standards.md.

Integrity Rules

1. Commit the model family before peeking at test folds.
2. Do not drop awkward seasons after seeing metrics.
3. Always beat a dumb baseline before claiming value.
4. Never use same-game score components as predictors for pre-game models.
5. Separate exploratory coefficient fishing from confirmatory walk-forward evaluation.
6. Report failures and non-improvements.
7. Distinguish confirmatory from exploratory analyses (see statistical integrity notes in references).
8. Make it reproducible — seasons, formula, package versions, seeds.

Bundled Resources

references/

scripts/

Related Skills

Quick Command Card

bashCopy code
python /path/to/statistical-modeling/scripts/assumption_checks.py \
  --input observations.csv --value-col margin --group-col venue \
  --out assumptions.json

python /path/to/statistical-modeling/scripts/glm_diagnostics.py \
  --input model_frame.csv \
  --formula "won ~ is_home + rating_diff + rest_diff" \
  --out glm_report.json