Chapter 39: Regression & Generalized Linear Models

The Population Model and OLS Estimator

Linear regression posits that the expected value of a response $Y$ is a linear function of predictors $\mathbf{x}$:

In matrix notation this is $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$ where $\mathbf{X}$ is the $n \times (p+1)$ design matrix (the first column is all ones for the intercept).

The Ordinary Least Squares (OLS) estimator minimises the residual sum of squares:

The Gauss-Markov theorem guarantees that among all linear unbiased estimators, OLS has the smallest variance — it is BLUE (Best Linear Unbiased Estimator) — provided the four LINE assumptions hold:

• Linearity: <!--MATHBLOCK10--> is a linear function of <!--MATHBLOCK11-->.
• Independence: observations are mutually independent.
• Normal errors: <!--MATHBLOCK12--> (needed for exact inference; asymptotically relaxable).
• Equal variance (homoscedasticity): <!--MATHBLOCK13--> for all <!--MATHBLOCK14-->.

Reading the Coefficient Summary

The fitted model returns much more than point estimates. The estimated standard error of each coefficient is:

Under $H_0: \beta_j = 0$, the test statistic $t_j = \hat{\beta}_j / \text{SE}(\hat{\beta}_j)$ follows a $t_{n-p-1}$ distribution, yielding the p-value shown in every regression output. A 95% confidence interval is $\hat{\beta}_j \pm t_{0.025, n-p-1} \cdot \text{SE}(\hat{\beta}_j)$.

The overall F-test checks whether the model explains any variance beyond the intercept-only baseline:

And the coefficient of determination:

Adjusted $R^2$ penalises model complexity and should always be preferred when comparing models of different sizes.

A Worked Example with statsmodels

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

rng = np.random.default_rng(42)
n = 200
df = pd.DataFrame({
    'sqft': rng.normal(1500, 400, n),
    'bedrooms': rng.integers(1, 6, n).astype(float),
    'age': rng.uniform(0, 50, n),
})
# True DGP: price = 50 + 0.12*sqft + 8*bedrooms - 0.5*age + noise
df['price'] = (
    50 + 0.12 * df['sqft'] + 8 * df['bedrooms']
    - 0.5 * df['age'] + rng.normal(0, 15, n)
)

model = smf.ols('price ~ sqft + bedrooms + age', data=df).fit()
print(model.summary())

The summary() output contains three tables. The top table reports model-level statistics (R², F-statistic, AIC, BIC, degrees of freedom). The middle table — the one practitioners spend most time on — shows each coefficient's estimate, standard error, t-statistic, p-value, and 95% confidence interval. The bottom table warns about potential issues such as high condition number (multicollinearity) or skewness in residuals.

Pitfall: interpreting coefficients without units. The coefficient on sqft of 0.12 means a one-square-foot increase in floor area is associated with a \$120 increase in price, holding bedrooms and age constant. The phrase "holding other predictors constant" is the ceteris-paribus clause that distinguishes multiple regression from simple bivariate correlation.

The Intercept

The intercept $\hat{\beta}_0$ is the model's predicted value when all predictors equal zero. In many practical models (e.g., age = 0, square footage = 0) this is extrapolation and the intercept has no substantive interpretation. Centering predictors ($x_j' = x_j - \bar{x}_j$) makes the intercept interpretable as the predicted response at the mean of all predictors and also improves numerical conditioning.

Why Diagnostics Matter

Fitting a linear model is easy; trusting its inferences requires verifying the LINE assumptions. Violations do not always invalidate predictions, but they do invalidate the standard errors, p-values, and confidence intervals — the very quantities that justify regression-based decisions. This section presents a systematic diagnostic workflow.

The Residual Toolkit

All diagnostics begin with the residuals $e_i = y_i - \hat{y}_i$ and their standardised variants.

The internally studentised residual is $r_i = e_i / (\hat{\sigma}\sqrt{1 - h_{ii}})$, where $h_{ii}$ is the $i$-th diagonal element of the hat matrix $\mathbf{H} = \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top$. The value $h_{ii}$ — called leverage — measures how far observation $i$ is from the centre of the predictor space. A point with high leverage has the potential to pull the regression line toward itself.

Cook's Distance combines leverage and residual size into a single influence metric:

A rule of thumb flags $D_i > 4/n$ as potentially influential.

The Four Diagnostic Plots

1. Residuals vs. Fitted. Plot $e_i$ against $\hat{y}_i$. A random scatter around zero indicates linearity and homoscedasticity. A U-shape suggests a missing quadratic term. A funnel shape (residuals spreading as $\hat{y}$ grows) signals heteroscedasticity — the most common real-world violation.

2. Normal Q-Q Plot. Plot the quantiles of standardised residuals against theoretical normal quantiles. Heavy tails or S-curves indicate non-normality. Mild non-normality matters little for large samples (CLT); severe skew or outliers warrant concern.

3. Scale-Location (Spread-Location). Plot $\sqrt{|r_i|}$ against $\hat{y}_i$. A flat LOESS smoother confirms homoscedasticity; a rising smoother confirms heteroscedasticity. This plot amplifies the funnel pattern.

4. Residuals vs. Leverage. Overlaying Cook's distance contours reveals which observations are both influential (high leverage) and poorly fitted (large residual). Points in the upper-right or lower-right corners of this plot warrant investigation.

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.graphics.gofplots as gof
import matplotlib.pyplot as plt
from statsmodels.stats.outliers_influence import OLSInfluence

rng = np.random.default_rng(7)
n = 150
df = pd.DataFrame({'x': rng.uniform(1, 10, n)})
# Introduce heteroscedasticity: variance grows with x
df['y'] = 3 + 2 * df['x'] + rng.normal(0, 0.5 * df['x'], n)

fit = smf.ols('y ~ x', data=df).fit()
infl = OLSInfluence(fit)

fig, axes = plt.subplots(2, 2, figsize=(10, 8))

# 1. Residuals vs Fitted
axes[0, 0].scatter(fit.fittedvalues, fit.resid, alpha=0.5, s=15)
axes[0, 0].axhline(0, color='red', lw=1)
axes[0, 0].set(xlabel='Fitted', ylabel='Residuals', title='Residuals vs Fitted')

# 2. Q-Q plot
gof.qqplot(fit.resid, line='s', ax=axes[0, 1], alpha=0.5)
axes[0, 1].set_title('Normal Q-Q')

# 3. Scale-Location
stud_resid = infl.resid_studentized_internal
axes[1, 0].scatter(fit.fittedvalues, np.sqrt(np.abs(stud_resid)), alpha=0.5, s=15)
axes[1, 0].set(xlabel='Fitted', ylabel='√|Studentized Resid|',
               title='Scale-Location')

# 4. Residuals vs Leverage
lev = infl.hat_matrix_diag
axes[1, 1].scatter(lev, stud_resid, alpha=0.5, s=15)
axes[1, 1].axhline(0, color='red', lw=1)
axes[1, 1].set(xlabel='Leverage', ylabel='Studentized Residuals',
               title='Residuals vs Leverage')

plt.tight_layout()
plt.savefig('diagnostics.png', dpi=120)

Formal Tests

While plots are primary, several formal tests supplement visual inspection:

• Breusch-Pagan / White test for heteroscedasticity: regresses <!--MATHBLOCK14--> on the predictors; a significant F-statistic rejects homoscedasticity.
• Shapiro-Wilk / Anderson-Darling for normality of residuals.
• Durbin-Watson for serial autocorrelation (values near 2 are ideal; < 1.5 or > 2.5 suggest correlation).

from statsmodels.stats.diagnostic import het_breuschpagan, linear_harvey_collier
from statsmodels.stats.stattools import durbin_watson

lm_stat, p_bp, _, _ = het_breuschpagan(fit.resid, fit.model.exog)
print(f'Breusch-Pagan p-value: {p_bp:.4f}')   # p < 0.05 → heteroscedasticity
dw = durbin_watson(fit.resid)
print(f'Durbin-Watson: {dw:.3f}')             # near 2 → no autocorrelation

Remedies

Once a violation is identified, standard remedies are:

• Heteroscedasticity: use HC3 robust standard errors (cov_type='HC3' in statsmodels), or transform <!--MATHBLOCK15--> (log, square root), or switch to Weighted Least Squares.
• Non-linearity: add polynomial terms, splines, or interaction terms.
• Influential outliers: investigate whether the observation is a data error or a genuinely extreme case; never delete automatically.
• Non-normality in small samples: use bootstrap confidence intervals.

The Multicollinearity Problem

Multicollinearity occurs when two or more predictors are highly linearly correlated. It does not bias $\hat{\boldsymbol{\beta}}$, but it inflates standard errors dramatically — sometimes to the point where no individual coefficient is statistically significant even though the model as a whole is highly significant. Geometrically, when $\mathbf{X}^\top\mathbf{X}$ is nearly singular, small changes in the data cause large swings in the coefficient estimates.

Variance Inflation Factor

The Variance Inflation Factor for predictor $j$ is:

where $R^2_j$ is the $R^2$ from regressing $x_j$ on all other predictors. A $\text{VIF}_j = 10$ means the variance of $\hat{\beta}_j$ is 10 times larger than it would be if $x_j$ were uncorrelated with the rest. Common thresholds: VIF > 5 is a warning; VIF > 10 is severe.

import pandas as pd
import numpy as np
from statsmodels.stats.outliers_influence import variance_inflation_factor
import statsmodels.formula.api as smf

rng = np.random.default_rng(3)
n = 200
z = rng.normal(0, 1, n)
# x2 is nearly a linear function of x1
df = pd.DataFrame({
    'x1': z + rng.normal(0, 0.05, n),
    'x2': z + rng.normal(0, 0.05, n),   # very collinear with x1
    'x3': rng.normal(0, 1, n),
})
df['y'] = 1 + 2*df['x1'] + 3*df['x3'] + rng.normal(0, 1, n)

exog = smf.ols('y ~ x1 + x2 + x3', data=df).fit().model.exog
vif = pd.DataFrame({
    'feature': ['intercept', 'x1', 'x2', 'x3'],
    'VIF': [variance_inflation_factor(exog, i) for i in range(exog.shape[1])]
})
print(vif.round(1))

Expected output reveals enormous VIFs for x1 and x2 while x3 remains low.

Remedies for Multicollinearity

• Drop redundant predictors when theory supports it.
• Combine collinear predictors into a composite (e.g., principal components).
• Collect more data — collinearity is partly a small-sample problem.
• Regularization — the preferred remedy in high-dimensional settings.

Ridge Regression

Ridge adds an $L_2$ penalty to the OLS objective:

The closed-form solution is:

Adding $\lambda\mathbf{I}$ to the (potentially near-singular) $\mathbf{X}^\top\mathbf{X}$ matrix conditions it away from singularity. Ridge shrinks all coefficients toward zero but never sets them exactly to zero. Crucially, predictors must be standardised before applying ridge, or the penalty is applied unequally.

Lasso Regression

Lasso uses an $L_1$ penalty:

Because the $L_1$ ball has corners at the coordinate axes, lasso solutions often have exactly zero components — lasso performs automatic feature selection. It is the preferred method when sparsity (many irrelevant predictors) is expected.

Elastic Net interpolates: $\lambda[\alpha\|\boldsymbol{\beta}\|_1 + (1-\alpha)\|\boldsymbol{\beta}\|^2]$, inheriting lasso's sparsity and ridge's grouping of correlated features.

import numpy as np
import pandas as pd
from sklearn.linear_model import RidgeCV, LassoCV, ElasticNetCV
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)
n, p = 200, 20
X = rng.normal(0, 1, (n, p))
# Only 5 of 20 features matter
true_beta = np.zeros(p)
true_beta[:5] = [2, -1.5, 1, -0.8, 0.6]
y = X @ true_beta + rng.normal(0, 1, n)

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)
sc = StandardScaler()
Xt = sc.fit_transform(X_train); Xv = sc.transform(X_test)

ridge = RidgeCV(alphas=np.logspace(-3, 3, 50), cv=5).fit(Xt, y_train)
lasso = LassoCV(alphas=np.logspace(-3, 1, 50), cv=5, random_state=0).fit(Xt, y_train)

print(f'Ridge best alpha: {ridge.alpha_:.4f},  Test R²: {ridge.score(Xv, y_test):.3f}')
print(f'Lasso best alpha: {lasso.alpha_:.4f},  Test R²: {lasso.score(Xv, y_test):.3f}')
print(f'Lasso non-zero coefs: {(lasso.coef_ != 0).sum()}')

Interpreting Regularised Coefficients

A critical warning: regularised coefficient estimates are biased and their standard errors from OLS theory are invalid. Do not report p-values from a lasso fit. For inference after variable selection, use post-selection inference methods (e.g., the selectiveInference R package) or report coefficients from a refit OLS on the selected subset — with appropriate multiple-testing caution.

From Linear to Generalized Linear Models

Linear regression assumes $Y \sim \mathcal{N}(\mu, \sigma^2)$. Real data rarely cooperates: binary outcomes violate this assumption fundamentally (predictions can exceed [0,1]), count data have non-negative integer support, and proportions are bounded. The Generalized Linear Model (GLM) unifies these cases with a single elegant framework.

A GLM has three components:

1. Random component: <!--MATHBLOCK5--> follows a distribution from the exponential family with mean <!--MATHBLOCK6-->.
2. Systematic component: a linear predictor <!--MATHBLOCK7-->.
3. Link function: a monotone differentiable function <!--MATHBLOCK8--> such that <!--MATHBLOCK9-->.

The GLM is estimated by Maximum Likelihood via Iteratively Reweighted Least Squares (IRLS), not OLS. Consequently, goodness-of-fit is measured by deviance rather than RSS:

The null deviance (intercept only) and residual deviance (fitted model) appear in every GLM summary, and their difference, divided by the degrees-of-freedom difference, follows a $\chi^2$ distribution under $H_0$.

Logistic Regression

For a binary outcome $Y_i \in \{0, 1\}$, the natural distributional family is Bernoulli with $\mu_i = P(Y_i = 1 | \mathbf{x}_i)$. The canonical link is the logit:

Inverting, the predicted probability is:

Interpreting Logistic Regression Coefficients

The coefficient $\hat{\beta}_j$ is the change in the log-odds of the event per one-unit increase in $x_j$, holding others constant. Exponentiated, $e^{\hat{\beta}_j}$ is the multiplicative change in the odds ratio:

An odds ratio of 1.5 means the odds of the event are 50% higher for each one-unit increase in $x_j$. A crucial caution: odds ratios are not probability ratios (relative risks), and the distinction matters clinically when the baseline probability is not near zero.

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from sklearn.datasets import load_breast_cancer

bc = load_breast_cancer()
df = pd.DataFrame(bc.data[:, :5], columns=bc.feature_names[:5])
df.columns = ['radius', 'texture', 'perimeter', 'area', 'smoothness']
df['malignant'] = (bc.target == 0).astype(int)  # 1 = malignant

# Standardise for numerical stability
for col in df.columns[:-1]:
    df[col] = (df[col] - df[col].mean()) / df[col].std()

fit = smf.logit(
    'malignant ~ radius + texture + perimeter + area + smoothness',
    data=df
).fit()
print(fit.summary())

# Odds ratios with 95% CI
odds = pd.DataFrame({
    'OR': np.exp(fit.params),
    'lower': np.exp(fit.conf_int()[0]),
    'upper': np.exp(fit.conf_int()[1]),
}).round(3)
print('\nOdds Ratios:')
print(odds)

Model Evaluation for Logistic Regression

Unlike linear regression, logistic $R^2$ analogues (McFadden, Nagelkerke) are not intuitive. Prefer:

• Likelihood-ratio test: compare nested models with <!--MATHBLOCK19-->.
• AUC-ROC: the area under the Receiver Operating Characteristic curve; 0.5 = random, 1.0 = perfect.
• Brier score: mean squared error of predicted probabilities; lower is better.
• Calibration: do predicted probabilities match observed proportions (use a calibration plot or Hosmer-Lemeshow test)?

from sklearn.metrics import roc_auc_score, brier_score_loss
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import LogisticRegression

X = df[['radius', 'texture', 'perimeter', 'area', 'smoothness']].values
y = df['malignant'].values

lr = LogisticRegression(max_iter=500)
print('CV AUC: ', cross_val_score(lr, X, y, cv=5, scoring='roc_auc').mean().round(3))
print('CV Brier:', cross_val_score(lr, X, y, cv=5,
      scoring='neg_brier_score').mean().__neg__().round(3))

The Wald Test vs. the Likelihood-Ratio Test

The summary table reports Wald z-statistics ($\hat{\beta}_j / \text{SE}(\hat{\beta}_j)$). These can be unreliable when the sample is small or when separation (a predictor perfectly predicts the outcome) causes coefficients to diverge. The likelihood-ratio test (LRT) is more reliable: fit the full and reduced models, compute $\Delta D$, and compare to $\chi^2$.

fit_reduced = smf.logit(
    'malignant ~ radius + texture + area + smoothness',  # drop perimeter
    data=df
).fit(disp=0)

lr_stat = fit_reduced.llr_pvalue  # full model vs null
delta_d = 2 * (fit.llf - fit_reduced.llf)  # compare full vs reduced
import scipy.stats as st
p_lrt = st.chi2.sf(delta_d, df=1)
print(f'LRT for perimeter: Δdeviance={delta_d:.2f}, p={p_lrt:.4f}')

When the Response is a Count

Count data — number of insurance claims, hospital readmissions, website visits, mutations per genome — are non-negative integers. Applying OLS to counts violates homoscedasticity (variance grows with mean for counts) and can produce negative predictions. The natural starting point is Poisson regression.

Poisson Regression

The Poisson GLM assumes $Y_i \sim \text{Poisson}(\mu_i)$ with the canonical log link:

Since $\text{Var}(Y_i) = \mu_i = E[Y_i]$ for the Poisson distribution, the model constrains variance to equal the mean — a strong assumption.

Interpretation of Poisson Coefficients

Exponentiating: $\mu_i = e^{\mathbf{x}_i^\top \boldsymbol{\beta}}$. A one-unit increase in $x_j$ multiplies the expected count by $e^{\hat{\beta}_j}$. This is the incidence rate ratio (IRR). An IRR of 1.3 means 30% more events per unit increase in $x_j$.

Exposure and Offsets

When counts are collected over different periods or population sizes, the expected count should be proportional to the exposure $E_i$:

The term $\log(E_i)$ is added as an offset — a predictor whose coefficient is fixed at 1.0 rather than estimated.

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

rng = np.random.default_rng(5)
n = 500
df = pd.DataFrame({
    'days_observed': rng.integers(30, 365, n),   # exposure varies
    'age_group': rng.choice(['young', 'middle', 'old'], n),
    'treatment': rng.integers(0, 2, n),
})
log_mu = (
    np.log(df['days_observed'])
    + np.where(df['age_group'] == 'young', 0.0,
        np.where(df['age_group'] == 'middle', 0.4, 0.9))
    - 0.5 * df['treatment']
)
df['events'] = rng.poisson(np.exp(log_mu))

fit = smf.poisson(
    'events ~ age_group + treatment + offset(np.log(days_observed))',
    data=df
).fit()
print(fit.summary())
print('\nIRR (exp(coef)):')
print(np.exp(fit.params).round(3))

Overdispersion: The Achilles Heel of Poisson

Real count data almost always exhibit overdispersion: $\text{Var}(Y_i) > E[Y_i]$. Overdispersion inflates the Poisson deviance and deflates standard errors, producing overconfident (spuriously significant) inferences. Detect it by comparing the residual deviance to its degrees of freedom — a ratio $\gg 1$ is a warning sign — or with a formal test.

# Quasi-Poisson: scales SE by sqrt(dispersion) — corrects SEs but not fit
fit_quasi = smf.glm(
    'events ~ age_group + treatment + offset(np.log(days_observed))',
    data=df,
    family=smf.families.family.Poisson()
).fit()

dispersion = fit_quasi.pearson_chi2 / fit_quasi.df_resid
print(f'Dispersion parameter: {dispersion:.2f}')  # ≈1 is fine; >>1 = overdispersed

Negative Binomial Regression

The Negative Binomial model addresses overdispersion by introducing a per-observation random effect $\nu_i \sim \text{Gamma}(\theta, \theta)$, leading to:

As $\theta \to \infty$, the negative binomial converges to the Poisson. Smaller $\theta$ means more overdispersion. This model is preferable to quasi-Poisson because it is a proper likelihood, allowing AIC comparison and model selection.

import statsmodels.api as sm

# Simulate overdispersed count data
rng = np.random.default_rng(11)
n = 800
df2 = pd.DataFrame({
    'x1': rng.normal(0, 1, n),
    'x2': rng.choice([0, 1], n),
})
true_mu = np.exp(0.5 + 0.8 * df2['x1'] - 0.4 * df2['x2'])
# Negative binomial with theta=2 (moderately overdispersed)
df2['y'] = rng.negative_binomial(2, 2 / (2 + true_mu))

exog = sm.add_constant(df2[['x1', 'x2']])
nb_fit = sm.NegativeBinomial(df2['y'], exog).fit(disp=0)
print(nb_fit.summary())
print(f'\nAlpha (1/theta): {nb_fit.params["alpha"]:.3f}')  # alpha = 1/theta
print(f'IRR for x1: {np.exp(nb_fit.params["x1"]):.3f}')

Zero-Inflated Models

Some count datasets have excess zeros beyond what any Poisson or NB model can explain — e.g., the number of cigarettes smoked per day when many participants are non-smokers. Zero-Inflated Poisson (ZIP) and Zero-Inflated Negative Binomial (ZINB) models mix a point mass at zero with a count component. The statsmodels ZeroInflatedPoisson and ZeroInflatedNegativeBinomialP classes handle these.

Model Selection Checklist

• Start with Poisson and check residual deviance / degrees of freedom.
• If dispersion <!--MATHBLOCK16-->, switch to Negative Binomial or use quasi-Poisson standard errors.
• If the zero count is notably higher than the Poisson/NB prediction, consider zero-inflated or hurdle models.
• Use AIC to compare non-nested models (Poisson vs. NB); use LRT (<!--MATHBLOCK17--> on <!--MATHBLOCK18--> deviance) for nested models.

A Unified Workflow for GLM Analysis

Fitting a regression model is only the beginning. A principled workflow moves through six stages: exploratory data analysis, model specification, estimation, diagnostics, comparison/selection, and communication of results. This section ties those stages together and addresses the most common practical issues.

Information Criteria: AIC and BIC

When comparing non-nested models or models of different complexity, likelihood-ratio tests are not valid. Information criteria penalise log-likelihood by model complexity:

where $k$ is the number of estimated parameters. Lower is better. BIC penalises complexity more severely and tends toward sparser models; AIC minimises prediction error on new data (in expectation). For model selection in a predictive context, prefer AIC or cross-validation; for identifying the true data-generating model when it is in the candidate set, prefer BIC.

Warning: never interpret the absolute value of AIC — only differences between models on the same dataset are meaningful.

Cross-Validation for Regression Models

For predictive models, k-fold cross-validation on a held-out metric (RMSE, log-loss, AUC) is more reliable than information criteria because it estimates generalisation error directly and makes no parametric assumptions about model fit.

import numpy as np
import pandas as pd
from sklearn.linear_model import Ridge, Lasso, LogisticRegression
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import KFold, cross_val_score
from sklearn.datasets import load_diabetes, load_breast_cancer

# Regression comparison: Ridge vs. Lasso vs. OLS (alpha → 0)
X, y = load_diabetes(return_X_y=True)
kf = KFold(n_splits=10, shuffle=True, random_state=42)

for name, est in [
    ('Ridge',  Pipeline([('sc', StandardScaler()), ('m', Ridge(alpha=1.0))])),
    ('Lasso',  Pipeline([('sc', StandardScaler()), ('m', Lasso(alpha=0.5))])),
    ('OLS',    Pipeline([('sc', StandardScaler()), ('m', Ridge(alpha=1e-6))])),
]:
    scores = cross_val_score(est, X, y, cv=kf, scoring='neg_root_mean_squared_error')
    print(f'{name:6s}: RMSE = {(-scores.mean()):.2f} ± {scores.std():.2f}')

Feature Engineering Inside the Model Formula

Both statsmodels (via Patsy) and R's formula interface support inline transformations, making the specification concise and reproducible.

import statsmodels.formula.api as smf
import pandas as pd, numpy as np

rng = np.random.default_rng(0)
n = 300
df = pd.DataFrame({
    'age':    rng.uniform(20, 70, n),
    'income': rng.lognormal(10, 0.5, n),
    'gender': rng.choice(['M', 'F'], n),
})
df['y'] = (
    0.02 * df['age']**2 - 2 * df['age']
    + 0.3 * np.log(df['income'])
    + np.where(df['gender'] == 'F', 5, 0)
    + rng.normal(0, 8, n)
)

# Patsy handles np.log, I() for arithmetic, and C() for categorical
fit = smf.ols(
    'y ~ I(age**2) + age + np.log(income) + C(gender, Treatment("M"))',
    data=df
).fit()
print(fit.summary().tables[1])

Communicating Results

A regression table published in a paper or technical report should include:

• Coefficient estimates with standard errors (not just p-values).
• Confidence intervals.
• Sample size and degrees of freedom.
• Model fit statistics (R², adjusted R², AIC, or deviance for GLMs).
• For logistic/count models: effect sizes on the natural scale (ORs, IRRs) with 95% CIs.

The table is not the end. Supplement with effect plots showing the predicted response as one predictor varies while others are held at their median values — these communicate non-linearity and interaction effects that coefficient tables obscure.

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

rng = np.random.default_rng(8)
n = 400
df = pd.DataFrame({
    'dose':    rng.uniform(0, 10, n),
    'weight':  rng.normal(70, 15, n),
    'treated': rng.integers(0, 2, n),
})
logit_p = -3 + 0.6 * df['dose'] - 0.02 * df['weight'] + 1.5 * df['treated']
df['response'] = rng.binomial(1, 1 / (1 + np.exp(-logit_p)))

fit = smf.logit('response ~ dose + weight + treated', data=df).fit(disp=0)

# Effect plot: response probability vs. dose, at median weight, both treatment arms
dose_grid = np.linspace(0, 10, 200)
med_weight = df['weight'].median()

fig, ax = plt.subplots(figsize=(6, 4))
for t, color, label in [(0, 'steelblue', 'Control'), (1, 'tomato', 'Treated')]:
    pred_df = pd.DataFrame({'dose': dose_grid,
                            'weight': med_weight,
                            'treated': t})
    preds = fit.get_prediction(pred_df).summary_frame()
    ax.plot(dose_grid, preds['mean'], color=color, label=label)
    ax.fill_between(dose_grid, preds['mean_ci_lower'], preds['mean_ci_upper'],
                    alpha=0.2, color=color)
ax.set_xlabel('Dose'); ax.set_ylabel('P(Response)')
ax.set_title('Effect Plot: Dose vs. Response Probability')
ax.legend()
plt.tight_layout()
plt.savefig('effect_plot.png', dpi=120)

Common Pitfalls Summarised

• Confusing statistical significance with practical significance: a tiny effect can be highly significant with large <!--MATHBLOCK2-->. Always report effect sizes.
• Data leakage in cross-validation: standardisation scalers and imputers must be fit only on the training fold.
• Ignoring the scale: never compare AIC across datasets or models fit on different samples.
• Extrapolation: the model's predictions outside the observed covariate range are unreliable.
• Post-selection inference: selecting variables based on p-values and then reporting those p-values as valid is anti-conservative. Use Bonferroni corrections, FWER, or FDR procedures, or pre-register your model.