| name | statistical-modeling |
| description | Best-practice suggestions for statistical / inferential modelling in Python — OLS and logistic regression with statsmodels, robust standard errors, side-by-side regression tables with stargazer, confidence and prediction intervals, hypothesis tests, and significance reporting. Use when interpreting coefficients, building explanatory regressions, comparing nested models, reporting confidence/prediction intervals, or testing whether an effect is statistically significant. |
Statistical Modelling Best Practices
These are suggestions, not absolute rules. The user is industry-experienced but wants to get sharper at the more academic side of modelling — confidence intervals, prediction intervals, significance tests, and clean regression-table reporting. The patterns here cover that gap.
The patterns are inspired by Békés–Kézdi (Data Analysis for Business, Economics, and Policy) and CEU coursework:
This skill is the inferential / explanatory counterpart to the ml-modeling skill. Use this when the question is "how big is the effect, and are we sure?". Use ml-modeling when the question is "how accurately can we predict?". They overlap on significance testing of model differences — see reference/inference_vs_prediction.md.
When to Use
Auto-apply when the task involves:
- Estimating and interpreting regression coefficients (not just predicting).
- Confidence intervals on parameters or fitted values.
- Prediction intervals for individual observations.
- Side-by-side comparison of nested or competing regression specifications.
- Hypothesis tests (t-test, F-test for joint significance, Wald, likelihood ratio).
- Robust / heteroscedasticity-consistent standard errors.
- Linear splines, log-log specifications, fixed effects on cross-sections.
- Reading a colleague's regression notebook and asking "is this defensible?".
Default Stack
statsmodels — smf.ols, smf.logit, sm.OLS. Use the formula API (smf) by default; falls back to sm.OLS when you need a design matrix.stargazer — side-by-side regression tables (the package used across da_case_studies and CEU Coding 2).scipy.stats — t-tests, F-tests, Wilcoxon, etc. when you need raw tests.pandas, numpy, matplotlib, seaborn for the rest.ydata-profiling (formerly pandas-profiling) for automated EDA reports — especially useful for the exploratory phase before specifying a model.missingno for missing-data visualisation.
For predictive workflows (cross-validation, held-out test sets, leaderboards) → switch to the ml-modeling skill and use scikit-learn. For project setup (folder structure, environments) → see the analytics-project-setup skill.
Core Principles
- Always report standard errors and confidence intervals, never just point estimates.
- Use robust standard errors by default for cross-sectional data —
cov_type="HC3" (or HC0/HC1 to match older textbooks). Don't trust homoscedasticity unless you've checked.
- Distinguish confidence interval from prediction interval. CI is uncertainty around the expected value; PI is uncertainty around an individual observation and is wider.
- Stargazer tables make comparisons honest. When showing >1 specification, put them side by side with the same dependent variable.
- Test significance explicitly — don't infer from "the coefficient is large". A coefficient is significantly different from zero only when the test says so.
- Distinguish statistical from practical significance. A
p < 0.001 effect of $0.02 may be irrelevant. Always report effect size in the units stakeholders understand.
- Pre-specify the model. If you ran 50 specifications and report the 3 with the smallest p-values, you have nothing.
OLS — the workhorse
import statsmodels.formula.api as smf
reg = smf.ols("price ~ distance + rating", data=hotels).fit(cov_type="HC3")
print(reg.summary())
print(reg.params)
print(reg.conf_int(alpha=0.05))
print("p-values:", reg.pvalues)
print(reg.f_test("distance = 0, rating = 0"))
Common specification choices:
| Goal | Specification |
|---|
| Linear | y ~ x |
| Log-level | np.log(y) ~ x (interpret coef as ~% change in y per unit x) |
| Level-log | y ~ np.log(x) (interpret coef as change in y per 1% change in x ≈ coef/100) |
| Log-log (elasticity) | np.log(y) ~ np.log(x) (coef = elasticity directly) |
| Polynomial | y ~ x + I(x**2) + I(x**3) |
| Linear spline | y ~ lspline(x, [knot1, knot2]) (helper from da_case_studies ch00) |
| Categorical | y ~ C(category) |
| Interaction | y ~ x * z (includes x, z, x:z) |
The lspline / knot_ceil helpers live in snippets/spline_helpers.py, adapted from py_helper_functions.py.
Logistic Regression (when the dependent is binary)
logit = smf.logit("default ~ income + age + C(region)", data=df).fit(cov_type="HC3")
print(logit.summary())
import numpy as np
ors = np.exp(logit.params)
ci = np.exp(logit.conf_int())
print(pd.concat([ors.rename("OR"), ci], axis=1))
For pure prediction with logistic regression at scale, switch to sklearn.linear_model.LogisticRegression (ml-modeling skill).
Side-by-side Tables with Stargazer
The pattern from ch13-used-cars-reg:
from stargazer.stargazer import Stargazer
reg1 = smf.ols("price ~ distance", data=hotels).fit(cov_type="HC3")
reg2 = smf.ols("price ~ distance + rating", data=hotels).fit(cov_type="HC3")
reg3 = smf.ols("price ~ distance + rating + C(stars)", data=hotels).fit(cov_type="HC3")
sg = Stargazer([reg1, reg2, reg3])
sg.rename_covariates({"distance": "Distance to centre", "rating": "User rating"})
sg.show_model_numbers(False)
sg.title("Hotel price models")
sg
A worked template lives in snippets/stargazer_table.py.
Confidence vs Prediction Intervals
This trip-up is the single most common one. From ch09-hotels-europe-stability:
- Confidence interval (CI): uncertainty around the fitted value (i.e.
E[Y | X = x]). Narrow.
- Prediction interval (PI): uncertainty around a single new observation — adds the residual variance on top. Wider.
pred = reg.get_prediction(new_X)
summary = pred.summary_frame(alpha=0.05)
When plotting fitted values, draw both intervals as differently shaded bands and label them. See snippets/intervals_plot.py.
Hypothesis Tests You'll Actually Use
| Question | Test |
|---|
| Is this single coefficient ≠ 0? | t-test (built into summary()) |
| Are these k coefficients jointly 0? | F-test: reg.f_test("a = 0, b = 0") |
| Is the bigger model actually better than the nested one? | Likelihood-ratio test or compare adj-R²/BIC; for nested OLS use anova_lm(small, big) |
| Are two means / two groups different? | scipy.stats.ttest_ind (independent) or ttest_rel (paired) |
| Same as above but non-normal | scipy.stats.mannwhitneyu or wilcoxon |
| Are two proportions different? | statsmodels.stats.proportion.proportions_ztest or chi-square |
| Is the residual heteroscedastic? | Breusch–Pagan or White: het_breuschpagan, het_white |
| Is there autocorrelation? | Durbin–Watson (already in summary()) or acorr_ljungbox |
A code catalogue with one-liners for each lives in reference/tests_catalogue.md.
When You Compare Models (not coefficients)
This is the bridge to ml-modeling:
- Nested OLS models (one is a subset of the other): use F-test via
statsmodels.stats.anova.anova_lm(small, big).
- Non-nested: compare on adj-R², AIC, BIC. The model with lower AIC/BIC wins; there's no formal test.
- Predictive ML models on the same CV folds: use the paired t-test on per-fold scores pattern from the
ml-modeling skill.
Significance Reporting Etiquette
When reporting an effect, give the reader four things in one sentence:
The coefficient on experience is 0.034 log-points (95% CI: 0.022–0.046; p < 0.001), implying an additional year of experience is associated with roughly 3.4% higher wages.
That's: estimate, CI, p-value, plain-English interpretation. Always include the unit and the practical magnitude. A p-value alone is not a finding.
Anti-Patterns to Flag
- Reporting only the point estimate, no SE / CI.
- Using default homoscedastic SE on cross-sectional data without checking residuals.
- Reporting
R² = 0.87 without showing test/holdout performance — high R² on training data is uninformative.
- Picking the "best" specification after looking at p-values across many runs (HARKing / p-hacking).
- Calling a coefficient "significant" without saying at what α and with what test.
- Confusing confidence interval with prediction interval in plots.
- Comparing models by R² alone — penalise complexity (adj-R², BIC) or compare on a held-out set.
- Treating logistic-regression coefficients as probabilities (they're log-odds).
Automated EDA / Profiling Before Modelling
Before specifying a model, run an automated profile to catch data quality issues:
from ydata_profiling import ProfileReport
profile = ProfileReport(df, minimal=True)
profile.to_file('eda_report.html')
For quick missing-data diagnostics: import missingno as msno; msno.matrix(df) — shows nullity patterns at a glance.
Data Preparation Utilities
Standardise column names early (df.columns = df.columns.str.replace('-', '_').str.replace(' ', '_').str.lower()). For full column-locking helpers (name_lock, type_lock), see the analytics-project-setup skill.
Code Snippets
In snippets/:
spline_helpers.py — knot_ceil and lspline from py_helper_functions.py for use inside smf.ols("y ~ lspline(x, [...])", ...).stargazer_table.py — full stargazer table template with renamed covariates, custom rows, BIC/AIC footer.intervals_plot.py — fitted line + confidence band + prediction band for OLS, properly labeled.regression_diagnostics.py — residual plots, heteroscedasticity tests, influence diagnostics in one call.
Data Sources Used in the Inspiration Courses
Almost all worked examples across da_case_studies, da_data_repo, and CEU Coding 2
da_data_repo, and da_case_studies come from Békés & Kézdi, Data Analysis for Business,
Economics, and Policy (Cambridge, 2021). The companion data is openly available and is the
single best collection of cleanly-documented small/medium datasets for inferential modeling practice:
Specific datasets you'll keep returning to (all small, all interpretable):
Other free, well-suited sources for inferential exercises:
| Source | What | URL |
|---|
| Statsmodels built-ins | sm.datasets.{anes96, fair, longley, macrodata, statecrime, ...} for quick OLS/Logit demos | import statsmodels.api as sm; sm.datasets.longley.load_pandas().data |
| Seaborn datasets | tips, diamonds, mpg, penguins — small enough to teach, big enough to find effects | import seaborn as sns; sns.load_dataset("diamonds") |
| OpenML | Tabular datasets with metadata; great for replication exercises | https://www.openml.org/ |
| World Bank Open Data | Country-year panels for cross-country regressions | https://data.worldbank.org/ |
| FRED (Federal Reserve Economic Data) | Macro time series for ARIMA / interrupted time-series | https://fred.stlouisfed.org/ via pandas_datareader |
When loading from OSF, the canonical pattern from the da_* repos works directly in Python:
import pandas as pd
hotels = pd.read_csv("https://osf.io/r6uqb/download")
Further Reference
Inspiration repos (check these for full worked examples):
- gabors-data-analysis/da_case_studies — most relevant chapters: ch07-hotels-simple-reg, ch08-hotels-nonlinear, ch09-hotels-europe-stability, ch10-hotels-multiple-reg, ch11-smoking-health-risk, ch13-used-cars-reg, ch17-predicting-firm-exit.
- CEU Coding 2 — gentle intro with
smf.ols, robust SE, log-log, splines, stargazer; private coursework.
- gabors-data-analysis/da_data_repo — companion datasets for the Békés–Kézdi book.
Companion skills:
ml-modeling — for predictive evaluation (cross-validation, leaderboards, threshold tuning).analytics-project-setup — for project scaffolding and AGENTS.md.data-warehousing — for bronze/silver/gold pipeline patterns.
External:
- Békés & Kézdi, Data Analysis for Business, Economics, and Policy (the textbook the
da_* repos accompany).
- Géron, Hands-On Machine Learning with Scikit-Learn, Keras and TensorFlow, 3rd ed. — ch. 2–4 cover complementary data preparation and model training practices.
- statsmodels user guide.
- Stargazer for Python.
Suggestions, not gospel. When in doubt, show the standard errors and the practical magnitude, not just the point estimate and the p-value.
