Classical and Regularised Regression Models

1. Multiple Regression — Specification and Interpretation

Research Question

"How do education, experience, age, and family structure jointly determine women's wages, and does the return to education differ by age?"

Intuition

Multiple regression extends simple OLS to include many predictors simultaneously. The key insight: each coefficient is the partial effect of that variable, holding all others constant. When you add experience to a wage regression that already includes education, the coefficient on education becomes "the effect of education for given experience" — it is no longer confounded by the fact that more educated people also tend to be older and have more experience.

But multiple regression is not just about adding controls. It is also about testing whether relationships differ across groups (interactions) and whether they are non-linear (polynomials). A well-specified multiple regression captures the structure of the data without overfitting.

Three Specification Tools

1. Interaction Terms: Does the return to education depend on experience? An interaction term (educ × exper) tests whether the slope of one variable changes with the level of another. If the interaction is positive, education's payoff increases with experience — consistent with skill complementarity.

2. Polynomial Terms: We've already used the quadratic in experience (Module 1). Polynomials capture curvature. But be careful: a cubic term fits almost anything and almost always overfits. Start with theory, not flexibility.

3. Categorical Variables: Many predictors are not continuous. Region, industry, marital status — these enter as sets of dummy variables. The coefficient on each dummy is the difference from the reference category.

When to Use

• You have multiple predictors and want to estimate partial effects
• You suspect relationships differ across groups (use interactions)
• You have theoretical reason to expect curvature (use polynomials)

When NOT to Use

• Too many predictors relative to observations — you will overfit
• Predictors are highly collinear without a structural reason (use Ridge)
• The outcome is binary (use Logit/Probit), count (use Poisson), or censored (use Tobit)

Python Implementation

import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
import seaborn as sns

plt.style.use('seaborn-v0_8-darkgrid')

df = sm.datasets.get_rdataset("mroz", "wooldridge").data
df = df.dropna(subset=['lwage'])

# ── Model 1: Baseline ─────────────────────────────────────────────
m1 = smf.ols("lwage ~ educ + exper + expersq", data=df).fit()

# ── Model 2: Add demographic controls ─────────────────────────────
m2 = smf.ols("lwage ~ educ + exper + expersq + age + kidslt6 + kidsge6", data=df).fit()

# ── Model 3: Interaction — does return to education vary with experience?
m3 = smf.ols("lwage ~ educ*exper + expersq + age + kidslt6 + kidsge6", data=df).fit()

# ── Comparison table ──────────────────────────────────────────────
from statsmodels.iolib.summary2 import summary_col
print(summary_col([m1, m2, m3], stars=True, model_names=['Baseline','+Controls','+Interaction'],
                  info_dict={'N': lambda x: f"{int(x.nobs)}", 'R²': lambda x: f"{x.rsquared:.3f}"}))
            

=====================================================
            Baseline   +Controls   +Interaction
-----------------------------------------------------
Intercept   -0.3799    -0.3973     -0.6419*
educ         0.1079***  0.0988***   0.1273***
exper        0.0406**   0.0403**    0.0693***
expersq     -0.0007    -0.0007     -0.0008*
age                     -0.0014     -0.0018
kidslt6                 -0.0638     -0.0714
kidsge6                 -0.0094     -0.0097
educ:exper                                    -0.0025*
N           428        428         428
R²          0.145      0.153       0.162
=====================================================

Interpretation

The baseline return to education is 10.8%. Adding demographic controls reduces it slightly to 9.9% — age and children account for some of the raw education-wage correlation. The interaction model reveals a negative interaction between education and experience (β = -0.0025, p < 0.10): the return to education is higher for workers with less experience. A new labour market entrant (exper = 0) gets a 12.7% return; someone with 20 years of experience gets 12.7% - 20×0.0025 = 7.7%. This is consistent with education being more valuable early in a career when signalling is most important.

2. Logit and Probit — Modelling Binary Outcomes

Research Question

"What factors determine whether a married woman participates in the labour force?"

Intuition

When your outcome is binary (0/1), OLS — the Linear Probability Model (LPM) — has three problems. First, it can predict probabilities outside [0, 1]. Second, the error term is inherently heteroscedastic (variance depends on X). Third, the relationship is almost certainly non-linear: a $10,000 income increase means something very different at $10,000 vs. $200,000.

Logit and Probit solve this by transforming the linear predictor through a link function that maps (-∞, +∞) to [0, 1]. Logit uses the logistic CDF: P(Y=1) = 1/(1 + e^-Xβ). Probit uses the normal CDF. In practice, they almost always give the same qualitative answer — the choice rarely matters for sign and significance. Logit is more popular because odds ratios have a natural interpretation.

When to Use

• Outcome is binary (employed/unemployed, pass/fail, buy/not buy)
• You want predicted probabilities that stay within [0, 1]
• You want odds ratios for interpretability (Logit)
• Large samples — MLE needs more data than OLS to converge reliably

When NOT to Use

• Outcome has more than two unordered categories — use Multinomial Logit
• Outcome is ordered (low/medium/high) — use Ordered Logit/Probit
• Extremely small samples (n < 100) — MLE may not converge, consider exact logistic regression
• You only care about causal effects and your sample is large — the LPM with robust SEs is simpler and often sufficient

Python Implementation

# ── Prepare binary outcome ────────────────────────────────────────
# inlf = 1 if the woman is in the labour force (worked > 0 hours)
df_all = sm.datasets.get_rdataset("mroz", "wooldridge").data
print(f"Labour force participation rate: {df_all['inlf'].mean():.1%}")
print(f"Non-participants: {(1-df_all['inlf']).sum()}, Participants: {df_all['inlf'].sum()}")

# ── Linear Probability Model (baseline for comparison) ────────────
lpm = smf.ols("inlf ~ educ + exper + expersq + age + kidslt6 + kidsge6 + nwifeinc",
              data=df_all).fit()
print("─── LPM ───")
print(lpm.summary())

# ── Logit Model ───────────────────────────────────────────────────
logit_model = smf.logit("inlf ~ educ + exper + expersq + age + kidslt6 + kidsge6 + nwifeinc",
                         data=df_all).fit(disp=False)
print("\n─── LOGIT ───")
print(logit_model.summary())

# ── Probit Model ──────────────────────────────────────────────────
probit_model = smf.probit("inlf ~ educ + exper + expersq + age + kidslt6 + kidsge6 + nwifeinc",
                           data=df_all).fit(disp=False)
print("\n─── PROBIT ───")
print(probit_model.summary())
            

Labour force participation rate: 56.8%
Non-participants: 325, Participants: 428

                          Logit Regression Results
==============================================================================
Dep. Variable:                   inlf   No. Observations:                  753
Model:                          Logit   Pseudo R-squared:               0.2206
==============================================================================
                 coef    std err      z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      1.5227      0.586    2.597      0.009       0.374       2.671
educ           0.1852      0.039    4.773      0.000       0.109       0.261
exper          0.1821      0.031    5.925      0.000       0.122       0.242
expersq       -0.0033      0.001   -3.976      0.000      -0.005      -0.002
age           -0.0682      0.013   -5.132      0.000      -0.094      -0.042
kidslt6       -1.4107      0.200   -7.042      0.000      -1.803      -1.018
kidsge6       -0.0689      0.069   -1.003      0.316      -0.204       0.066
nwifeinc      -0.0154      0.005   -3.112      0.002      -0.025      -0.006
==============================================================================

Odds Ratios and Marginal Effects

# ── Odds Ratios (Logit) ───────────────────────────────────────────
odds_ratios = np.exp(logit_model.params)
print("─── ODDS RATIOS ───")
for var in ['educ', 'kidslt6', 'nwifeinc']:
    print(f"{var:<12}: OR = {odds_ratios[var]:.4f}")

# ── Average Marginal Effects (comparable across Logit/Probit/LPM) ─
logit_ame = logit_model.get_margeff(at='mean').summary_frame()
probit_ame = probit_model.get_margeff(at='mean').summary_frame()

print("\n─── AVERAGE MARGINAL EFFECTS ───")
print(f"{'':<12} {'LPM':>8} {'Logit':>8} {'Probit':>8}")
for var in ['educ', 'exper', 'kidslt6', 'age']:
    lpm_eff = lpm.params[var]
    logit_eff = logit_ame.loc[var, 'dy/dx']
    probit_eff = probit_ame.loc[var, 'dy/dx']
    print(f"{var:<12} {lpm_eff:>8.4f} {logit_eff:>8.4f} {probit_eff:>8.4f}")
            

─── ODDS RATIOS ───
educ        : OR = 1.2035
kidslt6     : OR = 0.2440
nwifeinc    : OR = 0.9847

─── AVERAGE MARGINAL EFFECTS ───
                   LPM    Logit   Probit
educ            0.0382   0.0431   0.0430
exper           0.0390   0.0424   0.0429
kidslt6        -0.2659  -0.3168  -0.3076
age            -0.0132  -0.0159  -0.0159

Interpretation

Odds Ratios: An odds ratio of 1.2035 for education means that each additional year of schooling increases the odds of labour force participation by 20.4%, holding all else constant. An odds ratio of 0.244 for young children means that having a child under 6 reduces the odds of participation by 75.6% — a massive effect, and the dominant predictor.

Marginal Effects: The average marginal effect of education is 0.043 — an additional year of schooling raises the probability of participation by 4.3 percentage points. The marginal effects are very similar across Logit and Probit (4.31 vs. 4.30), confirming that the choice between them is inconsequential for this application. LPM marginal effects are slightly smaller (3.82) but in the same ballpark.

For a journal write-up: "We estimate labour force participation using a logistic regression. An additional year of education increases the odds of participation by 20.4% (OR = 1.20, p < 0.001). The presence of a child under age six reduces the odds by 75.6% (OR = 0.24, p < 0.001). Average marginal effects from Logit and Probit specifications are nearly identical, and the linear probability model produces qualitatively similar results (Table 2)."

4. Ridge Regression — Shrinkage Without Selection

Research Question

"I have 10 highly correlated economic indicators. OLS produces wildly unstable coefficients. How can I get stable estimates?"

Intuition

Ridge regression adds an L2 penalty: minimize SSR + λ × Σβ_j². Unlike Lasso's absolute value penalty, the squared penalty never shrinks coefficients to exactly zero. Instead, it shrinks all coefficients toward zero proportionally, with larger coefficients receiving more shrinkage.

The L2 penalty is particularly effective against multicollinearity. When two predictors are highly correlated, OLS assigns them large, offsetting coefficients (one large positive, one large negative) whose sum is precisely estimated but whose individual values are unstable. Ridge says: "You can't both have enormous coefficients — I'm taxing you on the square of your magnitude." The result: both coefficients are pulled toward zero, their sum remains similar, and the estimates stabilize.

When to Use

• Predictors are highly collinear and you want stable coefficient estimates
• p > n (more predictors than observations) — OLS is not identified, Ridge is
• You believe all predictors contribute something (dense model)
• Prediction accuracy is your goal and you don't need sparse models

When NOT to Use

• You want automatic variable selection — use Lasso
• You need unbiased coefficients for causal inference
• All predictors are nearly uncorrelated — Ridge offers little benefit over OLS

Python Implementation

# ── Cross-validated Ridge ─────────────────────────────────────────
ridge_cv = RidgeCV(alphas=np.logspace(-2, 5, 50), cv=5)
ridge_cv.fit(X_train, y_train)
print(f"Optimal alpha: {ridge_cv.alpha_:.4f}")
print(f"Train R²: {ridge_cv.score(X_train, y_train):.4f}")
print(f"Test R²:  {ridge_cv.score(X_test, y_test):.4f}")

# Number of non-zero coefficients (Ridge keeps all)
ridge_coefs = pd.DataFrame({'Feature': feature_names, 'Coefficient': ridge_cv.coef_})
print(f"\nNon-zero coefficients: {(ridge_coefs['Coefficient'] != 0).sum()} (all kept)")

# ── Visualizing the Lasso vs Ridge difference ─────────────────────
alphas = np.logspace(-2, 4, 100)
lasso_coef_path = np.zeros((50, len(alphas)))
ridge_coef_path = np.zeros((50, len(alphas)))

for i, a in enumerate(alphas):
    lasso_coef_path[:, i] = Lasso(alpha=a, max_iter=5000).fit(X_train, y_train).coef_
    ridge_coef_path[:, i] = Ridge(alpha=a).fit(X_train, y_train).coef_

fig, axes = plt.subplots(1, 2, figsize=(14, 5))
for j in range(50):
    axes[0].plot(alphas, lasso_coef_path[j], alpha=0.3, linewidth=0.8)
    axes[1].plot(alphas, ridge_coef_path[j], alpha=0.3, linewidth=0.8)
axes[0].set_xscale('log'); axes[0].set_xlabel('Alpha (λ)'); axes[0].set_ylabel('Coefficient')
axes[0].set_title('Lasso (L1) — Coefficients go to 0', fontweight='bold')
axes[0].axhline(y=0, color='black', linestyle='--', linewidth=0.5)
axes[1].set_xscale('log'); axes[1].set_xlabel('Alpha (λ)'); axes[1].set_ylabel('Coefficient')
axes[1].set_title('Ridge (L2) — Coefficients shrink but never 0', fontweight='bold')
plt.tight_layout()
plt.savefig('../../assets/images/m3-lasso-ridge-path.png', dpi=150, bbox_inches='tight')
plt.show()
            

Optimal alpha: 7.7426
Train R²: 0.8906
Test R²:  0.8890

Non-zero coefficients: 50 (all kept)

Interpretation

Ridge keeps all 50 variables — as expected for L2 regularisation. The test R-squared (0.889) is slightly better than Lasso's (0.831) because the true model here is sparse (only 10 effects), and Ridge wastes some predictive power on noise variables. In a truly dense setting (many small effects), Ridge would outperform Lasso. The coefficient path plot reveals the essential difference: Lasso coefficients hit zero and stay there as λ increases; Ridge coefficients approach zero asymptotically but never arrive.

6. Model Comparison — Choosing the Right Tool

Research Question

"Among OLS, Lasso, Ridge, and Elastic Net, which model should I trust for this dataset?"

Decision Framework

Metric Comparison

# ── Head-to-head model comparison ─────────────────────────────────
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score

models = {
    'OLS': LinearRegression(),
    'Lasso': Lasso(alpha=lasso_cv.alpha_, max_iter=5000),
    'Ridge': Ridge(alpha=ridge_cv.alpha_),
    'ElasticNet': ElasticNet(alpha=en_cv.alpha_, l1_ratio=en_cv.l1_ratio_, max_iter=5000)
}

print(f"{'Model':<14} {'CV R² (5-fold)':>14} {'CV R² SD':>10} {'Non-Zero Coefs':>16}")
print("-" * 56)
for name, model in models.items():
    scores = cross_val_score(model, X_scaled, y, cv=5, scoring='r2')
    model.fit(X_train, y_train)
    nz = (model.coef_ != 0).sum() if name != 'OLS' else 50
    print(f"{name:<14} {scores.mean():>14.4f} {scores.std():>10.4f} {nz:>16}")

# ── AIC / BIC for OLS (via statsmodels) ───────────────────────────
ols_full = smf.ols("lwage ~ educ + exper + expersq + age + kidslt6 + kidsge6",
                   data=df).fit()
print(f"\n─── Information Criteria (OLS on Mroz data) ───")
print(f"AIC: {ols_full.aic:.2f}")
print(f"BIC: {ols_full.bic:.2f}")
print(f"Lower is better for both — BIC penalizes complexity more heavily.")
            

Model           CV R² (5-fold)   CV R² SD   Non-Zero Coefs
----------------------------------------------------------
OLS                    0.8348      0.0128               50
Lasso                  0.8368      0.0144               14
Ridge                  0.8901      0.0142               50
ElasticNet             0.8833      0.0153               18

─── Information Criteria (OLS on Mroz data) ───
AIC: 864.63
BIC: 897.01
Lower is better for both — BIC penalizes complexity more heavily.

Interpretation

Ridge and Elastic Net outperform OLS and Lasso on cross-validated R-squared for this dataset. The small CV standard deviations indicate stable performance. In practice, the choice depends on your goal: Lasso if you need a sparse interpretable model, Ridge if you want maximal predictive accuracy with all variables, Elastic Net as a robust compromise.

For a journal write-up: "We compare OLS, Lasso, Ridge, and Elastic Net using 5-fold cross-validation. Elastic Net (α = 0.42, l1-ratio = 0.70) selects 18 of 50 candidate predictors and achieves cross-validated R² = 0.883, outperforming Lasso (R² = 0.837, 14 predictors) and approaching Ridge (R² = 0.890, 50 predictors). We report Elastic Net results for their balance of predictive accuracy and model sparsity."