statistical-modeling

# 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: - [gabors-data-analysis/da_case_studies](https://github.com/gabors-data-analysis/da_case_studies) — chapter notebooks for OLS, splines, multiple regression, and prediction ([ch07](https://github.com/gabors-data-analysis/da_case_studies/tree/master/ch07-hotels-simple-reg)–[ch10](https://github.com/gabors-data-analysis/da_case_studies/tree/master/ch10-hotels-multiple-reg), [ch09-hotels-europe-stability](https://github.com/gabors-data-analysis/da_case_studies/tree/master/ch09-hotels-europe-stability), [ch13-used-cars-reg](https://github.com/gabors-data-analysis/da_case_studies/tree/master/ch13-used-cars-reg)). Heavy use of `statsmodels` + `stargazer`. - [da_data_repo on OSF](https://osf.io/3u5em/) — companion datasets for the textbook. - CEU Coding 2 — `intro_to_regression.ipynb` (`smf.ols`, robust SE, `stargazer` side-by-side tables, linear splines); private coursework. - [ch17-predicting-firm-exit](https://github.com/gabors-data-analysis/da_case_studies/tree/master/ch17-predicting-firm-exit) for the link to predictive evaluation under business loss. The inferential-versus-predictive distinction is explained in [reference/inference_vs_prediction.md](reference/inference_vs_prediction.md). ## 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](https://github.com/gabors-data-analysis/da_case_studies) and CEU Coding 2). - `scipy.stats` — t-tests, F-tests, Wilcoxon, etc. when you need raw tests. - `pandas`, `numpy`, and `matplotlib` for data preparation and diagnostics. ## Core Principles 1. **Always report standard errors and confidence intervals**, never just point estimates. 2. **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. 3. **Distinguish confidence interval from prediction interval.** CI is uncertainty around the *expected value*; PI is uncertainty around an *individual observation* and is wider. 4. **Stargazer tables make comparisons honest.** When showing >1 specification, put them side by side with the same dependent variable. 5. **Test significance explicitly** — don't infer from "the coefficient is large". A coefficient is *significantly different from zero* only when the test says so. 6. **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. 7. **Pre-specify the model.** If you ran 50 specifications and report the 3 with the smallest p-values, you have nothing. ## OLS — the workhorse ```python import statsmodels.formula.api as smf reg = smf.ols("price ~ distance + rating", data=hotels).fit(cov_type="HC3") print(reg.summary()) # Coefficients with 95% CI print(reg.params) print(reg.conf_int(alpha=0.05)) # Coefficient-level test print("p-values:", reg.pvalues) # Joint test of multiple coefficients (F-test) 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](https://github.com/gabors-data-analysis/da_case_studies/tree/master/ch00-tech-prep)) | | Categorical | `y ~ C(category)` | | Interaction | `y ~ x * z` (includes `x`, `z`, `x:z`) | The `lspline` / `knot_ceil` helpers live in [snippets/spline_helpers.py](snippets/spline_helpers.py), adapted from [da_helper_functions.py](https://github.com/gabors-data-analysis/da_case_studies/blob/master/ch00-tech-prep/da_helper_functions.py). ## Logistic Regression (when the dependent is binary) ```python logit = smf.logit("default ~ income + age + C(region)", data=df).fit(cov_type="HC3") print(logit.summary()) # Odds ratios with 95% CI 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](https://github.com/gabors-data-analysis/da_case_studies/tree/master/ch13-used-cars-reg): ```python 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 # in Jupyter, renders HTML; .render_latex() for LaTeX ``` A worked template lives in [snippets/stargazer_table.py](snippets/stargazer_table.py). ## Confidence vs Prediction Intervals This trip-up is the single most common one. From [ch09-hotels-europe-stability](https://github.com/gabors-data-analysis/da_case_studies/tree/master/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. ```python pred = reg.get_prediction(new_X) summary = pred.summary_frame(alpha=0.05) # columns: mean, mean_se, mean_ci_lower, mean_ci_upper, ← CI on E[Y|X] # obs_ci_lower, obs_ci_upper ← PREDICTION interval on Y ``` When plotting fitted values, draw both intervals as differently shaded bands and label them. See [snippets/intervals_plot.py](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](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). ## 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/](snippets/): - `spline_helpers.py` — `knot_ceil` and `lspline` from [da_helper_functions.py](https://github.com/gabors-data-analysis/da_case_studies/blob/master/ch00-tech-prep/da_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](https://github.com/gabors-data-analysis/da_case_studies), [da_data_repo on OSF](https://osf.io/3u5em/), and CEU Coding 2 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: | Resource | What | URL | |---|---|---| | Book site (datasets, code in R/Python/Stata, case studies) | Hub for everything | <https://gabors-data-analysis.com/> | | `da_data_repo` (clean datasets; no raw data) | Hosted on OSF — the GitHub mirror is gone | <https://osf.io/3u5em/> · parent project <https://osf.io/7epdj/> | | `da_case_studies` (chapter notebooks) | End-to-end workflows per chapter | <https://github.com/gabors-data-analysis/da_case_studies> | Specific datasets you'll keep returning to (all small, all interpretable): | Dataset | Typical use | Direct URL | |---|---|---| | `hotels-europe` (Vienna subset) | OLS on hotel price; log specs; splines; CIs vs PIs | <https://osf.io/r6uqb/> | | `cps-earnings` (CPS earnings) | Mincer-style log-wage regression, gender gap, robust SE | <https://osf.io/g8p9j/> | | `wms-management-survey` | Cross-country firm management practices regressions | <https://osf.io/uzpce/> | | `used-cars` | Hedonic price regression with feature engineering | <https://osf.io/7gvz9/> | | `share-health` | Logistic regression on health outcomes | <https://osf.io/yhdte/> | | `bisnode-firms` | Probability prediction (probit/logit), classification thresholds with business loss | <https://osf.io/3qyut/> | | `football-managers` | Difference-in-differences, event studies | <https://osf.io/r2psv/> | 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/> |