Qualitative Predictors, Interaction Effects, Structural Breaks, Polynomial Splines & Model Selection

In classical regression modeling, the columns of the design matrix $\mathbf{X} \in \mathbb{R}^{n \times p}$ must be numerical vectors, yet many important explanatory factors in economics, biostatistics, sociology, and quality engineering are **qualitative (categorical or nominal) variables** possessing $G \ge 2$ discrete, unordered categories or levels $\{C_1, C_2, \dots, C_G\}$. Assigning arbitrary consecutive integers $1, 2, \dots, G$ to a nominal variable inside a single regression column would impose a false metric assumption that the difference between category $1$ and category $2$ is identical in magnitude and direction to the difference between category $2$ and category $3$. To incorporate a $G$-level qualitative factor into a linear regression model without imposing any artificial ordering or metric spacing among the categories, we construct **indicator (dummy) variables**. For a qualitative factor with $G$ mutually exclusive and exhaustive levels, we introduce $G - 1$ binary indicator variables $D_1, D_2, \dots, D_{G-1}$ alongside the model intercept $\beta_0$. Under the standard **Reference Cell (or Corner-Point) Coding** convention, one category—say $C_G$ (or $C_1$)—is designated as the **baseline reference group** and assigned zeros across all $G - 1$ indicators, while each remaining category $C_g$ ($g = 1, \dots, G - 1$) activates its corresponding indicator: $$D_{ig} = \begin{cases} 1, & \text{if observation } i \text{ belongs to category } C_g \\ 0, & \text{otherwise} \end{cases} \qquad (g = 1, 2, \dots, G - 1)$$ When these $G - 1$ dummy variables are included alongside a continuous covariate $X_i$ in an additive multiple regression model: $$Y_i = \beta_0 + \beta_1 x_i + \gamma_1 D_{i1} + \gamma_2 D_{i2} + \dots + \gamma_{G-1} D_{i, G-1} + \varepsilon_i$$ evaluating the conditional expectation function $\mathbb{E}[Y_i \mid x_i, \text{Group}]$ for each category reveals the exact structural interpretation of every parameter: - **Reference Category $C_G$ ($D_{i1} = \dots = D_{i, G-1} = 0$)**: $\mathbb{E}[Y_i \mid x_i, C_G] = \beta_0 + \beta_1 x_i$. Thus, $\beta_0$ is the baseline intercept for the reference category $C_G$. - **Category $C_g$ ($D_{ig} = 1$, all other $D_{im} = 0$)**: $\mathbb{E}[Y_i \mid x_i, C_g] = (\beta_0 + \gamma_g) + \beta_1 x_i$. Thus, $\gamma_g = \mathbb{E}[Y \mid x, C_g] - \mathbb{E}[Y \mid x, C_G]$ measures the **vertical parallel shift (differential intercept)** between category $C_g$ and the reference category $C_G$ holding the continuous predictor $x$ constant! Moreover, the contrast between any two non-reference categories $C_g$ and $C_m$ is simply $\gamma_g - \gamma_m$, with standard error $s\sqrt{c_{gg} + c_{mm} - 2c_{gm}}$. ### The Dummy Variable Trap (Exact Linear Collinearity) Why must a model with an intercept $\beta_0$ include **at most $G - 1$** dummy variables for a $G$-level factor rather than all $G$ indicators $D_1, \dots, D_G$? Because the $G$ categories are mutually exclusive and exhaustive, every observation $i = 1, \dots, n$ belongs to exactly one category, which forces the pointwise algebraic identity: $$\sum_{g=1}^G D_{ig} = D_{i1} + D_{i2} + \dots + D_{iG} = 1 = x_{i0} \quad \forall i = 1, 2, \dots, n$$ In vector notation, the sum of all $G$ indicator columns equals the intercept column of ones: $\sum_{g=1}^G \mathbf{d}_{(g)} = \mathbf{1}_n$. If an analyst includes both the intercept column $\mathbf{1}_n$ and all $G$ dummy columns in $\mathbf{X}$, the columns of $\mathbf{X}$ are linearly dependent ($\text{rank}(\mathbf{X}) = G < G + 1$), the Gram matrix $\mathbf{X}^T\mathbf{X}$ is strictly singular ($\det(\mathbf{X}^T\mathbf{X}) = 0$), and unique Ordinary Least Squares estimation collapses! This exact collinearity pathology is universally known in econometrics as the **Dummy Variable Trap**. There are three mathematically equivalent ways to avoid the Dummy Variable Trap: 1. **Reference Cell Coding ($1$ Intercept $+$ $G - 1$ Dummies)**: Retain $\beta_0$ and omit one reference dummy $D_G$ (equivalent to imposing the corner-point constraint $\gamma_G = 0$). 2. **Cell-Means Coding (No Intercept $+$ All $G$ Dummies)**: Omit the intercept $\beta_0$ and include all $G$ indicator columns $D_1, \dots, D_G$: $$Y_i = \mu_1 D_{i1} + \mu_2 D_{i2} + \dots + \mu_G D_{iG} + \beta_1 x_i + \varepsilon_i$$ Because the $G$ dummy columns are mutually orthogonal ($\mathbf{d}_{(g)}^T \mathbf{d}_{(m)} = 0$ for $g \neq m$ and $\mathbf{d}_{(g)}^T \mathbf{d}_{(g)} = n_g$), this parameterization is full rank and estimates each group's intercept $\mu_g = \beta_0 + \gamma_g$ directly. 3. **Effects (Sum-to-Zero / Deviation) Coding**: Retain the grand intercept $\mu_0$ and define $G - 1$ contrast variables $E_{ig} \in \{-1, 0, +1\}$ corresponding to the ANOVA zero-sum restriction $\sum_{g=1}^G \alpha_g = 0$ (so $\alpha_G = -\sum_{g=1}^{G-1} \alpha_g$): $$E_{ig} = D_{ig} - D_{iG} = \begin{cases} +1, & \text{if observation } i \in C_g \\ -1, & \text{if observation } i \in C_G \\ 0, & \text{otherwise} \end{cases} \qquad (g = 1, \dots, G - 1)$$ Under effects coding, the intercept $\mu_0 = \frac{1}{G}\sum_{g=1}^G \mu_g$ represents the unweighted grand mean across all $G$ categories, and each coefficient $\alpha_g = \mu_g - \mu_0$ measures the deviation of category $C_g$ from the unweighted grand mean rather than from a single reference cell. Crucially, because all three coding schemes span the exact same $G$-dimensional column space $\mathcal{C}(\mathbf{d}_{(1)}, \dots, \mathbf{d}_{(G)})$, the projection matrix $\mathbf{H}$, the fitted values $\hat{\mathbf{Y}}$, the residuals $\mathbf{e}$, and the $\text{SSE}$ are strictly identical across all three parameterizations!

The additive dummy-variable model of Section 7.1 assumes that the marginal effect $\beta_1 = \frac{\partial \mathbb{E}[Y \mid x, D]}{\partial x}$ of the continuous predictor $X$ is identical across all $G$ qualitative groups, producing $G$ parallel regression lines that differ only in their vertical intercepts. In statistical biometry and experimental design, this parallel-lines specification is the classical Analysis of Covariance (ANCOVA) model: by adjusting for a continuous baseline concomitant variable $X$ (such as pre-treatment blood pressure or initial plot fertility) that is unaffected by the treatment, ANCOVA removes within-group error variance, reduces $\text{MSE}$, and adjusts treatment group comparisons $\hat{\gamma}_g = (\bar{Y}_{g\cdot} - \bar{Y}_{G\cdot}) - \hat{\beta}_1 (\bar{x}_{g\cdot} - \bar{x}_{G\cdot})$ for baseline covariate imbalances.

However, in many scientific settings, the slope with respect to $X$ varies across qualitative categories—an effect known as a qualitative-by-quantitative interaction (or moderation / non-parallelism). Geometrically, an interaction means that the partial derivative of $\mathbb{E}[Y \mid X_1, X_2]$ with respect to one predictor depends on the level of another predictor. Forming the cross-product interaction terms $(D_{ig} x_i)$ for $g = 1, \dots, G - 1$ yields the Unrestricted Separate-Regressions (Interaction) Model:

Evaluating the conditional expectation function within each category reveals how both the intercept and the slope shift across regimes:

• Reference Group $C_G$: $\mathbb{E}[Y_i \mid x_i, C_G] = \beta_0 + \beta_1 x_i$ (intercept $\beta_0$, slope $\beta_1$).
• Group $C_g$ ($g = 1, \dots, G - 1$): $\mathbb{E}[Y_i \mid x_i, C_g] = (\beta_0 + \gamma_g) + (\beta_1 + \delta_g) x_i$ (intercept $\beta_0 + \gamma_g$, slope $\beta_1 + \delta_g$).

Here, $\gamma_g$ is the differential intercept and $\delta_g$ is the differential slope of group $g$ relative to the reference group $G$. This nested hierarchy enables three distinct General Linear Hypothesis $F$-tests:

1. Test for Homogeneity of Slopes (Parallelism / No Interaction): $H_0: \delta_1 = \delta_2 = \dots = \delta_{G-1} = 0$ ($q = G - 1$ degrees of freedom), comparing the Full interaction model against the Reduced parallel-lines ANCOVA model.

2. Test for Equal Group Intercepts Given Parallelism (ANCOVA Treatment Effect): Assuming parallelism holds, $H_0: \gamma_1 = \dots = \gamma_{G-1} = 0$ ($q = G - 1$ degrees of freedom).

3. Test for Complete Coincidence of All $G$ Regression Lines: $H_0: \gamma_1 = \dots = \gamma_{G-1} = \delta_1 = \dots = \delta_{G-1} = 0$ ($q = 2(G - 1)$ degrees of freedom), testing whether a single pooled regression line $Y_i = \beta_0 + \beta_1 x_i + \varepsilon_i$ holds across all $G$ groups.

The Chow Test for Structural Stability Across Regimes ($p$ Parameters)

In macroeconomics, finance, and policy evaluation, an investigator often divides a sample of size $n = n_1 + n_2$ into two regimes—such as pre-policy vs. post-policy periods ($t = 1, \dots, n_1$ and $t = n_1 + 1, \dots, n$), or urban vs. rural households—and wishes to test whether the entire $p \times 1$ regression parameter vector remains structurally stable across the two regimes:

Assuming homoscedastic, independent normal disturbances $\varepsilon_i \overset{\text{iid}}{\sim} \mathcal{N}(0, \sigma^2)$ across both regimes with $n_1 > p$ and $n_2 > p$, Gregory Chow (1960) showed that testing $H_0: \boldsymbol{\beta}_{(1)} = \boldsymbol{\beta}_{(2)}$ is a General Linear $F$-test with $q = p$ restrictions:

• Reduced (Pooled) Model: Fit a single regression to all $n = n_1 + n_2$ observations pooled together to obtain the pooled residual sum of squares $\text{SSE}_R = \text{SSE}_{\text{pooled}}$ with $\text{df}_R = n_1 + n_2 - p$ degrees of freedom.
• Full (Unrestricted Separate Regimes) Model: Fit the regression separately to Regime 1 (yielding $\text{SSE}_1$ with $n_1 - p$ df) and separately to Regime 2 (yielding $\text{SSE}_2$ with $n_2 - p$ df). Because the two regimes have zero overlapping parameters in the unrestricted block-diagonal model $\begin{pmatrix} \mathbf{Y}_1 \\ \mathbf{Y}_2 \end{pmatrix} = \begin{pmatrix} \mathbf{X}_1 & \mathbf{0} \\ \mathbf{0} & \mathbf{X}_2 \end{pmatrix} \begin{pmatrix} \boldsymbol{\beta}_{(1)} \\ \boldsymbol{\beta}_{(2)} \end{pmatrix} + \boldsymbol{\varepsilon}$, the unrestricted Full-model residual sum of squares is simply the sum $\text{SSE}_F = \text{SSE}_1 + \text{SSE}_2$ with $\text{df}_F = (n_1 - p) + (n_2 - p) = n - 2p$ degrees of freedom!

Notice that this block-diagonal model is algebraically identical to defining a regime dummy $D_i \in \{0, 1\}$ and interacting $D_i$ with the intercept and all $k = p - 1$ predictors in a single pooled equation! Therefore, the exact Chow Structural Break $F$-Statistic is:

(Note: If the second regime has fewer observations than parameters, $n_2 \le p$, $\text{SSE}_2$ cannot be estimated separately; in that case, Chow's Predictive Stability Test compares $\text{SSE}_{\text{pooled}} - \text{SSE}_1$ against $\text{SSE}_1 / (n_1 - p)$ using an $F_{n_2, \, n_1 - p}$ distribution).

When theoretical considerations or residual plots indicate that the conditional expectation function $\mu(x) = \mathbb{E}[Y \mid X = x]$ is a smooth curve rather than a straight line, the classical Weierstrass Approximation Theorem motivates approximating $\mu(x)$ on a bounded interval by a $d$-th degree Polynomial Regression Model:

Even though $\mu(x)$ is a nonlinear curve in $x$, the polynomial regression model is strictly linear in the unknown parameter vector $\boldsymbol{\beta} = (\beta_0, \beta_1, \dots, \beta_d)^T \in \mathbb{R}^{d+1}$! Defining the $n \times (d + 1)$ Vandermonde Design Matrix $\mathbf{X}$ with entries $X_{ij} = x_i^j$ ($j = 0, 1, \dots, d$) casts polynomial regression directly into the standard multiple linear regression apparatus $\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$.

Two Critical Statistical and Numerical Principles in Polynomial Regression

1. The Hierarchy (Marginalitya) Principle: Because transforming the origin of $x$ via $x^* = x - c$ expands $(x^* + c)^2 = (x^*)^2 + 2c x^* + c^2$, a quadratic polynomial without a linear term ($\beta_1 = 0$) in one coordinate system acquires a non-zero linear term in every shifted coordinate system! To ensure that the fitted curve is invariant to linear shifts of scale and origin, polynomial models must respect the hierarchical inclusion principle: whenever a higher-order power $x^d$ is retained in the model, all lower-order powers $1, x, x^2, \dots, x^{d-1}$ should be retained even if one of their individual $t$-statistics is small, and model degree $d$ should be selected via sequential (Type I) sums of squares $\text{SSR}(x^d \mid x, \dots, x^{d-1})$.

2. Severe Multicollinearity and Ill-Conditioning of the Raw Vandermonde Basis: When $x_i > 0$ (for instance, years $x_i \in [2000, 2025]$ or temperatures $x_i \in [300, 400]$), the raw power columns $x_i, x_i^2, x_i^3$ are almost perfectly linearly correlated ($r(x, x^2) > 0.99$)! Even on $[0, 1]$, the continuous limit of $\frac{1}{n}\mathbf{X}^T\mathbf{X}$ as $n \to \infty$ is the notorious Hilbert matrix $H_{jk} = \int_0^1 x^{j+k} dx = \frac{1}{j + k + 1}$, whose condition number grows exponentially ($\kappa(\mathbf{H}_5) \approx 4.8 \times 10^5$), causing severe variance inflation and catastrophic numerical instability for $d \ge 4$.

Centering and Forsythe–Gram–Schmidt Orthogonal Polynomials

A simple first remedy is mean-centering the predictor via $\tilde{x}_i = x_i - \bar{x}$ (or scaling to $[-1, 1]$) before raising to powers: $Y_i = \alpha_0 + \alpha_1 \tilde{x}_i + \alpha_2 \tilde{x}_i^2 + \dots + \alpha_d \tilde{x}_i^d + \varepsilon_i$. For any symmetric design around $\bar{x}$, odd and even powers of $\tilde{x}_i$ are strictly orthogonal ($\sum \tilde{x}_i^{2m-1} \tilde{x}_i^{2r} = 0$), immediately eliminating the correlation between $\tilde{x}_i$ and $\tilde{x}_i^2$!

To achieve complete mutual orthogonality across all degrees $j = 0, 1, \dots, d$ for any arbitrary spacing of sample points $x_1, \dots, x_n$, we replace the raw monomials $\{1, x, x^2, \dots, x^d\}$ by Discrete Orthogonal Polynomials $\{P_0(x), P_1(x), \dots, P_d(x)\}$, where $P_j(x)$ is a monic polynomial of exact degree $j$ satisfying the sample orthogonality condition:

Using Forsythe's (1957) Three-Term Recurrence Relation (discrete Gram–Schmidt orthogonalization), the orthogonal polynomials are generated recursively without ever forming ill-conditioned moments:

• Degree $0$: $P_0(x_i) = 1$ for all $i = 1, \dots, n$.
• Degree $1$: $P_1(x_i) = x_i - \bar{x}$ (since $\sum P_0(x_i) P_1(x_i) = \sum (x_i - \bar{x}) = 0$).
• Degree $j \ge 2$:

where the recurrence constants $a_j$ and $b_j$ are chosen so that $\sum_{i=1}^n P_j(x_i) P_{j-1}(x_i) = 0$ and $\sum_{i=1}^n P_j(x_i) P_{j-2}(x_i) = 0$ (and by induction $\sum P_j(x_i) P_m(x_i) = 0$ for all $m < j$):

In the orthogonal polynomial regression model $Y_i = \gamma_0 P_0(x_i) + \gamma_1 P_1(x_i) + \dots + \gamma_d P_d(x_i) + \varepsilon_i$, the Gram matrix $\mathbf{P}^T\mathbf{P} = \text{diag}(\sum P_0^2, \sum P_1^2, \dots, \sum P_d^2)$ is strictly diagonal! Consequently:

Because the columns are mutually orthogonal, adding a higher-degree polynomial $P_{d+1}(x)$ never alters the previously computed estimates $\hat{\gamma}_0, \hat{\gamma}_1, \dots, \hat{\gamma}_d$ or their sums of squares $\text{SSR}(P_j)$!

High-degree global polynomials ($d \ge 4$) suffer from a well-known analytical pathology known as Runge's Phenomenon: because a single global polynomial is infinitely rigid (analytic everywhere), a local perturbation or outlier near one end of the predictor domain induces wild, non-local oscillations across the entire interval, and polynomial tails diverge rapidly to $\pm\infty$ at the boundaries. When the underlying physical or biological mechanism changes slope smoothly or abruptly at interior threshold values, piecewise polynomials and regression splines provide flexible, locally adaptable curves within the exact linear regression framework.

1. Piecewise Linear (Broken-Stick / Segmented) Regression at a Known Knot $\xi$

Suppose the linear slope of $Y$ on $x$ changes at a known threshold (knot or join point) $x = \xi$ (for example, material yield above a critical stress $\xi$, or electricity demand above a comfort temperature $\xi = 22^\circ\text{C}$). If we fitted two completely separate regressions on $x \le \xi$ and $x > \xi$ using a step dummy $D_i = \mathbb{I}(x_i > \xi)$, the fitted function would suffer an unphysical vertical jump discontinuity at $x = \xi$ unless constrained. To guarantee exact $C^0$ continuity at the knot $x = \xi$ while allowing the slope to change from $\beta_1$ (for $x < \xi$) to $\beta_1 + \beta_2$ (for $x > \xi$), we define the hinge function (or degree-$1$ truncated power basis function, identical to the ReLU activation function in modern deep learning):

The Continuous Piecewise Linear Regression Model is then a standard two-predictor linear regression:

Let us verify continuity and slope interpretation on both sides of the knot $\xi$:

• Left Regime ($x \le \xi$): $(x - \xi)_+ = 0 \implies \mathbb{E}[Y \mid x] = \beta_0 + \beta_1 x$ (slope is $\beta_1$). Approaching the knot from the left gives $\lim_{x \to \xi^-} \mathbb{E}[Y \mid x] = \beta_0 + \beta_1 \xi$.
• Right Regime ($x > \xi$): $(x - \xi)_+ = x - \xi \implies \mathbb{E}[Y \mid x] = (\beta_0 - \beta_2 \xi) + (\beta_1 + \beta_2) x$ (slope is $\beta_1 + \beta_2$). Approaching the knot from the right gives $\lim_{x \to \xi^+} \mathbb{E}[Y \mid x] = \beta_0 + \beta_1 \xi + \beta_2(\xi - \xi) = \beta_0 + \beta_1 \xi$.

Both limits match at $x = \xi$, and testing $H_0: \beta_2 = 0$ via an ordinary $t$-test directly tests whether the slope changes at $\xi$!

2. Cubic Regression Splines and the Truncated Power Basis ($C^2$ Smoothness)

In general, a degree-$d$ regression spline with $K$ interior knots $\xi_1 < \xi_2 < \dots < \xi_K$ is a piecewise degree-$d$ polynomial on each subinterval $[\xi_k, \xi_{k+1}]$ that is constrained to be continuous together with its first $d - 1$ derivatives ($C^{d-1}$ smooth) at every knot $\xi_k$. The most widely used choice is the Cubic Regression Spline ($d = 3$), which is continuous in value ($f(\xi_k^-) = f(\xi_k^+)$), first derivative/slope ($f'(\xi_k^-) = f'(\xi_k^+)$), and second derivative/curvature ($f''(\xi_k^-) = f''(\xi_k^+)$) at every knot, making the joins imperceptible to the human eye.

How many free linear parameters does a cubic spline with $K$ knots possess? In the Truncated Power Basis, we start with the global cubic polynomial $1, x, x^2, x^3$ ($4$ parameters) and add one cubic truncated power function $(x - \xi_k)_+{}^3 = [\max(0, x - \xi_k)]^3$ at each knot $k = 1, \dots, K$:

Why does adding $(x - \xi_k)_+{}^3$ automatically preserve continuity of $f, f',$ and $f''$ while allowing the third derivative $f'''$ to jump by $6\theta_k$ at $x = \xi_k$? Because the function $h_k(x) = (x - \xi_k)_+{}^3$ has $h_k(\xi_k) = 0$, first derivative $h_k'(x) = 3(x - \xi_k)_+{}^2 \implies h_k'(\xi_k) = 0$, and second derivative $h_k''(x) = 6(x - \xi_k)_+ \implies h_k''(\xi_k) = 0$, while $h_k'''(x) = 6 \mathbb{I}(x > \xi_k)$ jumps by $6$! Thus, a cubic spline with $K$ interior knots is a linear regression model with exact dimension:

In computational practice (Course Stat H-312), the truncated power basis columns $(x - \xi_k)_+{}^3$ become collinear for large $x$, so software replaces them with Curry–Schoenberg B-Splines (Basis Splines)—an algebraically equivalent basis spanning the exact same $(K + 4)$-dimensional spline subspace, where each B-spline basis function $B_m(x)$ is non-zero only on a compact local support of at most $d + 2 = 5$ adjacent knots, yielding a well-conditioned banded Gram matrix! Furthermore, Natural Cubic Splines impose the boundary constraint that the curve is linear ($f''(x) = f'''(x) = 0$) beyond the boundary knots $\xi_1$ and $\xi_K$, saving $4$ degrees of freedom (reducing dimension from $K + 4$ to $K$) and eliminating erratic boundary tail variance.

When diagnostic residual plots reveal that the response variable $Y > 0$ violates both linearity and homoscedasticity/normality simultaneously—such as positive skewness coupled with an outward megaphone fan where $\text{Var}(Y_i)$ increases with the conditional mean $\mu_i = \mathbb{E}[Y_i \mid \mathbf{x}_i]$—transforming the response variable $Y \mapsto g(Y)$ often rectifies non-normality, heteroscedasticity, and nonlinearity in a single stroke.

1. First-Order Delta-Method Derivation of Variance-Stabilizing Transformations

Suppose $\mathbb{E}[Y_i] = \mu_i$ and the error standard deviation is a known function of the mean: $\text{Var}(Y_i) = \sigma^2 [h(\mu_i)]^2$. By a first-order Taylor expansion of a smooth transformation $g(Y_i)$ around $\mu_i$, $g(Y_i) \approx g(\mu_i) + g'(\mu_i)(Y_i - \mu_i)$, so the approximate variance of the transformed response is:

To make $\text{Var}(g(Y_i))$ approximately constant (independent of $\mu_i$), we choose $g(\mu)$ so that $g'(\mu) h(\mu) = 1$, i.e., $g(\mu) = \int \frac{d\mu}{h(\mu)}$:

• Poisson-Like Count Variance ($\text{Var}(Y_i) \propto \mu_i \implies h(\mu) = \mu^{1/2}$): $g(\mu) = \int \mu^{-1/2} d\mu = 2\sqrt{\mu} \implies$ Square-Root Transformation $Y^* = \sqrt{Y}$.
• Constant Coefficient of Variation / Multiplicative Error ($\text{Var}(Y_i) \propto \mu_i^2 \implies h(\mu) = \mu$): $g(\mu) = \int \mu^{-1} d\mu = \ln(\mu) \implies$ Logarithmic Transformation $Y^* = \ln(Y)$.
• Quadratic Variance ($\text{Var}(Y_i) \propto \mu_i^4 \implies h(\mu) = \mu^2$): $g(\mu) = \int \mu^{-2} d\mu = -\mu^{-1} \implies$ Reciprocal Transformation $Y^* = 1/Y$.

2. The Box–Cox Parametric Power Transformation Family (1964)

Rather than guessing a power transformation by trial and error, George Box and Sir David Cox (1964) embedded all power and logarithmic transformations of a strictly positive response $Y_i > 0$ into a single continuous parametric family indexed by $\lambda \in \mathbb{R}$:

Notice why subtracting $1$ and dividing by $\lambda$ is essential: by L'Hôpital's rule, $\lim_{\lambda \to 0} \frac{y^\lambda - 1}{\lambda} = \lim_{\lambda \to 0} \frac{y^\lambda \ln(y)}{1} = \ln(y)$, making $y^{(\lambda)}$ a continuous function of $\lambda$ at $\lambda = 0$! Special cases include $\lambda = 1$ (linear $y - 1$), $\lambda = 0.5$ (square root $2(\sqrt{y} - 1)$), $\lambda = 0$ (natural log $\ln y$), $\lambda = -0.5$ (inverse square root), and $\lambda = -1$ (reciprocal $1 - y^{-1}$).

How is the optimal transformation parameter $\hat{\lambda}$ estimated by Maximum Likelihood? Box and Cox postulate that for the true parameter $\lambda$, the transformed vector $\mathbf{y}^{(\lambda)}$ satisfies the classical normal linear model $\mathbf{y}^{(\lambda)} \sim \mathcal{N}_n(\mathbf{X}\boldsymbol{\beta}, \sigma^2 \mathbf{I}_n)$. However, we cannot simply pick $\lambda$ to minimize $\text{SSE}(\lambda) = \|(\mathbf{I}_n - \mathbf{H})\mathbf{y}^{(\lambda)}\|_2^2$, because changing $\lambda$ changes the scale of $y_i^{(\lambda)}$ (for instance, shrinking $\lambda$ can change $\text{SSE}(\lambda)$ purely by compressing the numerical units of $y^{(\lambda)}$)! Instead, we must write the exact joint density of the original untranformed data $(Y_1, \dots, Y_n)$ by multiplying the normal density of $\mathbf{y}^{(\lambda)}$ by the Jacobian determinant of the transformation $y_i \mapsto y_i^{(\lambda)}$:

Profiling out $\hat{\boldsymbol{\beta}}(\lambda) = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T \mathbf{y}^{(\lambda)}$ and $\hat{\sigma}_{\text{MLE}}^2(\lambda) = \frac{1}{n}\text{SSE}(\lambda)$ yields the Box–Cox Profile Log-Likelihood:

Even more conveniently, let $\dot{y} = \left(\prod_{i=1}^n y_i\right)^{1/n} = \exp\left(\frac{1}{n}\sum_{i=1}^n \ln y_i\right)$ denote the geometric mean of the original responses, and define the Jacobian-normalized Box–Cox response $z_i^{(\lambda)}$:

Because the Jacobian of $y_i \mapsto z_i^{(\lambda)}$ is identically $1$ ($\prod_{i=1}^n (y_i / \dot{y})^{\lambda - 1} = 1$), maximizing $\ell_p(\lambda)$ is algebraically equivalent to minimizing the Residual Sum of Squares $\text{SSE}_z(\lambda)$ of the normalized response $z_i^{(\lambda)}$:

An approximate $100(1 - \alpha)\%$ likelihood-ratio confidence interval for $\lambda$ consists of all values of $\lambda$ satisfying $2[\ell_p(\hat{\lambda}) - \ell_p(\lambda)] \le \chi_{\alpha, 1}^2$ (i.e., $\ell_p(\lambda) \ge \ell_p(\hat{\lambda}) - \frac{1}{2}\chi_{\alpha, 1}^2 = \ell_p(\hat{\lambda}) - 1.921$ for $95\%$ confidence), from which the analyst selects a nearby interpretable value such as $\lambda \in \{-1, -0.5, 0, 0.5, 1\}$.

Suppose an analyst has collected a pool of $K$ candidate explanatory variables ($P = K + 1$ total parameters including the intercept), giving rise to $2^K$ possible subset regression models. Why not simply include all $K$ predictors in a single "kitchen-sink" model? As proved in Problem 2.2, while omitting relevant variables ($\beta_j \neq 0$) induces omitted-variable bias, including extraneous or redundant variables ($\beta_j = 0$) inflates the sampling variances $\text{Var}(\hat{\beta}_m)$ of the remaining coefficients, inflates the out-of-sample prediction variance $\text{Var}(\hat{Y}_0) = \sigma^2 \mathbf{x}_0^T(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{x}_0$, and wastes error degrees of freedom. To balance goodness of fit (low $\text{SSE}_p$) against model parsimony (small parameter count $p = k + 1 \le P$) across candidate subset models of dimension $p$, statistical theory employs five primary **Model Selection Criteria**: ### 1. Adjusted Coefficient of Determination ($R_{a,p}^2$) and Minimum $\text{MSE}_p$ As established in Section 2.6, maximizing $R_{a,p}^2 = 1 - \frac{\text{SSE}_p / (n - p)}{\text{SST} / (n - 1)}$ is strictly equivalent to **minimizing the residual mean square $\text{MSE}_p = \frac{\text{SSE}_p}{n - p}$**. ### 2. Colin Mallows' $C_p$ Criterion (1973) Mallows designed $C_p$ to estimate the standardized **Total Scaled Mean Squared Prediction Error** $\Gamma_p$ of a $p$-parameter subset model across all $n$ sample design points: $$\Gamma_p = \frac{1}{\sigma^2} \sum_{i=1}^n \mathbb{E}\left[ (\hat{Y}_{i(p)} - \mu_i)^2 \right] = \frac{1}{\sigma^2} \left( \sum_{i=1}^n \text{Var}(\hat{Y}_{i(p)}) + \sum_{i=1}^n \left[\text{Bias}(\hat{Y}_{i(p)})\right]^2 \right)$$ Because $\sum_{i=1}^n \text{Var}(\hat{Y}_{i(p)}) = \sigma^2 \text{tr}(\mathbf{H}_p) = p \sigma^2$, and $\mathbb{E}[\text{SSE}_p] = (n - p)\sigma^2 + \sum_{i=1}^n [\text{Bias}(\hat{Y}_{i(p)})]^2$, solving for the squared bias sum $\sum [\text{Bias}(\hat{Y}_{i(p)})]^2 = \mathbb{E}[\text{SSE}_p] - (n - p)\sigma^2$ and substituting into $\Gamma_p$ yields: $$\Gamma_p = p + \frac{\mathbb{E}[\text{SSE}_p] - (n - p)\sigma^2}{\sigma^2} = \frac{\mathbb{E}[\text{SSE}_p]}{\sigma^2} - (n - 2p)$$ Replacing $\mathbb{E}[\text{SSE}_p]$ by the observed subset $\text{SSE}_p$ and replacing the unknown true variance $\sigma^2$ by the unbiased variance estimator $\hat{\sigma}_{\text{full}}^2 = \text{MSE}_P = \frac{\text{SSE}_P}{n - P}$ from the **full $P$-parameter model** defines **Mallows' $C_p$ Statistic**: $$C_p = \frac{\text{SSE}_p}{\text{MSE}_P} - (n - 2p) = p + (P - p) \left( F_{P - p} - 1 \right)$$ where $F_{P-p} = \frac{(\text{SSE}_p - \text{SSE}_P)/(P - p)}{\text{MSE}_P}$ is the partial $F$-statistic for testing whether the $P - p$ omitted variables have zero coefficients! Consequently: - If the $p$-parameter subset model has **zero specification bias** ($\mathbb{E}[\text{SSE}_p] = (n - p)\sigma^2$), then $\mathbb{E}[C_p] \approx p$. - For the full model ($p = P$), $C_P = P$ identically by algebra. - Good subset models exhibit **both small $p$ and $C_p \le p$** (sitting near or below the $45^\circ$ reference line $C_p = p$). ### 3. Akaike Information Criterion ($\text{AIC}$), Hurvich–Tsai $\text{AIC}_c$, and Schwarz Bayesian Information Criterion ($\text{BIC}$) Information criteria penalize twice the negative maximized log-likelihood $-2 \ell(\hat{\boldsymbol{\beta}}_p, \hat{\sigma}_{p,\text{MLE}}^2) = n \ln(\text{SSE}_p / n) + n(1 + \ln 2\pi)$ by a complexity penalty proportional to $p$. Dropping additive constants that depend only on $n$, for a linear regression model of dimension $p$: $$\text{AIC}(p) = n \ln\left(\frac{\text{SSE}_p}{n}\right) + 2p, \qquad \text{AIC}_c(p) = \text{AIC}(p) + \frac{2(p+1)(p+2)}{n - p - 2}, \qquad \text{BIC}(p) = n \ln\left(\frac{\text{SSE}_p}{n}\right) + p \ln(n)$$ When $K$ is too large for exhaustive **All-Possible-Regressions** ($2^K$ subsets, accelerated via the Furnival–Wilson leaps-and-bounds branch-and-bound algorithm), greedy heuristic searches—**Forward Selection**, **Backward Elimination**, and **Stepwise Regression** (which alternates forward additions with backward partial-$F$ deletions to remove variables rendered redundant by newly entered collinear predictors)—are employed, though analysts must remember that post-selection $p$-values from stepwise searches suffer from severe selection bias (Freedman's paradox) and should be validated on held-out test data or replaced by LASSO regularization.