Tutorials › Exam SRM - Statistics for Risk Modeling › Why Mallows’ Cp Estimates Prediction Error

Exam SRM · Information Criteria · Tutorial 182 of 926

Why Mallows’ Cp Estimates Prediction Error

See how residual error and a complexity adjustment combine to estimate a candidate model’s prediction error under classical linear-model assumptions.

Intermediate 10 min read

What You'll Learn

  • State the classical assumptions behind the Mallows’ Cp prediction-error relationship
  • Distinguish prediction error for the mean response from error for a future observation
  • Derive the relationship between candidate-model RSS and expected prediction error
  • Explain why the full-model residual variance estimate is used in Cp
  • Count model dimension correctly, including the intercept
  • Identify when a plug-in Cp value is an estimate rather than an exact unbiased quantity

Why a Complexity Adjustment Estimates Prediction Error

A candidate regression model can have a small residual sum of squares simply because it uses many coefficients to fit the observed responses. Prediction error asks a different question: how far, on average, will predictions from the fitted model be from the true mean response, or from a future response? Mallows’ \(C_p\) adjusts the candidate model’s residual sum of squares for the optimism created by fitting coefficients to the same data.

The key comparison is between training residual error and prediction error at the same design points. A fitted model has already adapted to its training responses, so its training residuals do not fully reveal the error expected on new data. Under the classical linear-model assumptions, the adjustment is \(2p\) error variances, where \(p\) is the candidate model’s dimension. The reference error variance is estimated from an adequately specified full model.

Definition: Let \(n\) be the number of observations and let a candidate least-squares model have design matrix \(X_M\) of rank \(p\). The rank \(p\) is the number of fitted coefficients, including the intercept when present. Its residual sum of squares is \(\operatorname{RSS}_M\). Given an estimate \(\widehat{\sigma}^2\) of the error variance from a full model, Mallows’ \(C_p\) is $$ C_p=\frac{\operatorname{RSS}_M}{\widehat{\sigma}^2}-n+2p. $$

In the classical setup, responses at fixed design points satisfy \(Y=\mu+\varepsilon\), where the true mean vector \(\mu\) is represented by an adequate full linear model and the errors are independent with mean zero and common variance \(\sigma^2\). The usual normal-error assumption, \(\varepsilon\sim N(0,\sigma^2 I)\), supports the standard inferential setting; the expectation calculations below require only the stated means, variances, and independence. Candidate models are fitted by least squares. If the candidate omits part of the true mean structure, that omission contributes prediction bias; it is not assumed away.

The full-model estimate is commonly \(\widehat{\sigma}^2=\operatorname{RSS}_F/(n-p_F)\), where \(p_F\) is the rank of the full-model design matrix. This estimate is unbiased for \(\sigma^2\) when the full model represents the true mean and the error assumptions hold. It supplies the scale for comparing candidate models: the candidate’s residual error is measured in units of the estimated error variance.

Deriving the Prediction-Error Relationship

Write \(H_M=X_M(X_M^\mathsf{T}X_M)^{-1}X_M^\mathsf{T}\) for the candidate model’s projection matrix when its columns are linearly independent. More generally, \(H_M\) is the projection onto the column space of \(X_M\). It is symmetric and idempotent, with \(\operatorname{tr}(H_M)=p\). The fitted response vector is \(H_MY\), and the residual sum of squares is \(Y^\mathsf{T}(I-H_M)Y\).

Let \(b_M=(I-H_M)\mu\). This is the part of the true mean vector that the candidate model cannot represent. Its squared length, \(b_M^\mathsf{T}b_M\), is the candidate’s squared approximation bias over the observed design points. The expected training residual sum of squares is therefore:

$$ \begin{aligned} E[\operatorname{RSS}_M] &=E\!\left[(\mu+\varepsilon)^\mathsf{T}(I-H_M)(\mu+\varepsilon)\right]\\ &=b_M^\mathsf{T}b_M+\sigma^2\operatorname{tr}(I-H_M)\\ &=b_M^\mathsf{T}b_M+(n-p)\sigma^2. \end{aligned} $$

The cross-product has expectation zero because \(E[\varepsilon]=0\). The variance term follows from \(E[\varepsilon^\mathsf{T}A\varepsilon]=\sigma^2\operatorname{tr}(A)\) for a fixed symmetric matrix \(A\). Here, \(\operatorname{tr}(I-H_M)=n-p\).

Now consider prediction of the mean response \(\mu\) at these same design points, using a fitted model based on the training responses. The prediction error vector is \(\mu-H_MY=b_M-H_M\varepsilon\). The bias and random-error cross-product again has expectation zero. Since \(H_M^2=H_M\), the expected squared prediction error is:

$$ \begin{aligned} E\!\left[\|\mu-H_MY\|^2\right] &=b_M^\mathsf{T}b_M+ E[\varepsilon^\mathsf{T}H_M\varepsilon]\\ &=b_M^\mathsf{T}b_M+\sigma^2\operatorname{tr}(H_M)\\ &=b_M^\mathsf{T}b_M+p\sigma^2. \end{aligned} $$

Subtract the expected training RSS from this prediction-error expression. The bias term cancels, leaving a difference of \(2p\sigma^2-n\sigma^2\). Equivalently, the prediction error is the training RSS plus \(2p\sigma^2\), in expectation. This is the reason for the complexity adjustment:

$$ E\!\left[\|\mu-H_MY\|^2\right] =E[\operatorname{RSS}_M]+2p\sigma^2-n\sigma^2. $$

Thus, if \(\sigma^2\) were known, the quantity \(\operatorname{RSS}_M+2p\sigma^2-n\sigma^2\) would be an unbiased estimator of the candidate model’s expected squared prediction error for the mean response at these design points. Dividing by \(\sigma^2\) gives the familiar \(C_p\) form with known variance. In practice, \(\sigma^2\) is replaced by the full-model estimate \(\widehat{\sigma}^2\), so \(C_p\) is a plug-in estimate. Because that estimate is itself random, the plug-in ratio is not generally exactly unbiased in finite samples.

If the target is a future response vector \(Y^0=\mu+\varepsilon^0\), independent of the training errors and with the same variance \(\sigma^2\), the expected squared prediction error includes an additional \(n\sigma^2\) of irreducible future-response noise. Therefore, the same Cp adjustment estimates the part of total prediction error beyond that noise floor. It is important to say which prediction target is being discussed.

Worked Examples

Worked Example: Translating Cp Into Prediction-Error Units

A hypothetical insurer models annual claim counts for \(n=20\) policy groups using a candidate linear model with \(p=3\) fitted coefficients, including the intercept. Its residual sum of squares is 70. An adequate full model gives an error-variance estimate of \(\widehat{\sigma}^2=4\). Calculate \(C_p\), then interpret it in squared-response units. Assume the classical linear-model setup described above.

Method. Use the plug-in definition of \(C_p\), then multiply the score by the estimated variance to recover the estimated prediction-error adjustment in the original squared units.

Step 1: Scale the candidate RSS. \(70/4=17.5\).

Step 2: Apply the sample-size and dimension terms. \(C_p=17.5-20+2(3)=17.5-20+6=3.5\).

Step 3: Return to squared-response units. The estimated prediction error for the mean response beyond the \(n\widehat{\sigma}^2\) noise baseline is \(C_p\widehat{\sigma}^2=3.5(4)=14\). The corresponding plug-in estimate for prediction of a future response vector is \(14+20(4)=94\).

Check. Directly, the adjusted estimate in squared units is \(70+2(3)(4)-20(4)=70+24-80=14\), agreeing with the Cp calculation. The common mistake avoided here is calling \(C_p=3.5\) an error of 3.5 squared-response units; it is scaled by the variance estimate.

Worked Example: How Omitted Signal Enters the Relationship

Suppose a fixed-design study of \(n=12\) hypothetical service regions has independent normal errors with variance \(\sigma^2=3\). For candidate model A, the squared omitted-mean component is \(b_A^\mathsf{T}b_A=8\) and its dimension is \(p_A=2\). For candidate model B, these quantities are 1 and \(p_B=4\). Find each model’s expected RSS, expected mean-response prediction error, and expected \(C_p\) if the true variance is known.

Method. Apply the two expectation formulas derived above. In this example the variance is known, so the theoretical Cp expression is used without a plug-in variance estimate.

Model A: expected RSS. \(8+(12-2)(3)=8+30=38\).

Model A: mean-response prediction error. \(8+2(3)=14\).

Model A: expected Cp. The theoretical expression is \(E[\operatorname{RSS}_A]/\sigma^2-n+2p_A\), giving \(38/3-12+4=38/3-8=14/3\), or approximately 4.6667.

Model B: expected RSS. \(1+(12-4)(3)=1+24=25\).

Model B: mean-response prediction error. \(1+4(3)=13\).

Model B: expected Cp. \(25/3-12+8=25/3-4=13/3\), or approximately 4.3333.

Check. Multiplying each expected Cp by \(\sigma^2=3\) gives 14 for A and 13 for B, exactly the respective expected mean-response prediction errors. Model B’s larger dimension raises its variance contribution, but its smaller omitted-signal bias more than offsets that increase. The common mistake avoided is assuming the candidate must contain the true mean: the bias term remains in the derivation when a candidate omits signal.

Worked Example: Estimating Error Variance From a Full Model

For \(n=10\) hypothetical portfolio segments, a full linear model has rank \(p_F=4\) and residual sum of squares 12. A candidate model has rank \(p=2\) and residual sum of squares 22. Calculate the full-model variance estimate and the candidate’s plug-in \(C_p\). Interpret the result, assuming the full model adequately represents the mean and the classical error assumptions hold.

Method. First estimate the error variance by dividing the full-model RSS by its residual degrees of freedom. Then use this estimate in the candidate Cp formula.

Step 1: Full-model residual degrees of freedom. \(n-p_F=10-4=6\).

Step 2: Estimate the variance. \(\widehat{\sigma}^2=12/6=2\).

Step 3: Calculate Cp. The candidate’s scaled RSS is \(22/2=11\). Therefore \(C_p=11-10+2(2)=11-10+4=5\).

Step 4: Interpret in squared units. The estimated mean-response prediction error beyond the noise baseline is \(5(2)=10\). Including the future-response noise baseline gives \(10+10(2)=30\) for a future response vector.

Check. The direct adjusted expression is \(22+2(2)(2)-10(2)=22+8-20=10\). This agrees with \(C_p\widehat{\sigma}^2\). The common mistake avoided is dividing the full-model RSS by \(n\), rather than by the full model’s residual degrees of freedom \(n-p_F\).

Assumptions, Interpretation, and a Common Trap

The argument depends on the candidate least-squares fit being a projection with rank \(p\), and on a common error variance across observations. The full model must adequately represent the true mean if its residual mean square is to estimate the error variance without systematic bias. If the full model omits important mean structure, its residual variation can include signal as well as random error, making \(\widehat{\sigma}^2\) an unsuitable reference scale. Unequal variances, dependence, or a different prediction design also change the calculation and require an appropriate criterion.

A frequent error is to use \(C_p=\operatorname{RSS}/\widehat{\sigma}^2+2p\) and omit the \(-n\) term. For the data in the first worked example, that incorrect calculation gives \(17.5+6=23.5\), not 3.5. The missing subtraction matters because the scaled residual sum of squares includes a baseline contribution associated with \(n\) observations. The full expression is \(\operatorname{RSS}/\widehat{\sigma}^2-n+2p\).

Another counting error is to let \(p\) mean only the number of slope predictors when an intercept is fitted. The projection rank counts every estimated coefficient. If a model has an intercept and two slopes, then \(p=3\); using \(p=2\) understates the adjustment by 2 on the Cp scale. Finally, the estimated variance is the full-model error variance, not the candidate model’s residual mean square. Using candidate residual variance would make the scale change with each candidate and would no longer implement the comparison derived here.

Key takeaway: Under the classical linear-model assumptions, \(C_p\) applies a \(2p\) complexity adjustment to the candidate RSS and uses the full-model error variance as its scale. With known variance, its expectation equals mean-response prediction error divided by \(\sigma^2\); with estimated variance, it is a plug-in estimate whose accuracy depends on the full-model variance estimate.

Check Your Understanding

For each problem, use the stated model dimension, count an intercept when present, and distinguish a Cp-scale value from squared-response units.

  1. A candidate model has \(n=15\), \(p=3\), \(\operatorname{RSS}=48\), and known error variance \(\sigma^2=4\). Calculate \(C_p\) and the corresponding estimated mean-response prediction error in squared units. (Answer: \(C_p=3\); prediction error \(=12\).)
  2. A full model has \(n=18\), \(p_F=5\), and \(\operatorname{RSS}_F=26\). A candidate has \(p=2\) and \(\operatorname{RSS}=32\). Calculate the full-model variance estimate, candidate \(C_p\), and adjusted prediction error in squared units. (Answer: \(\widehat{\sigma}^2=2\), \(C_p=2\), adjusted prediction error \(=4\).)
  3. With \(n=20\) and \(\widehat{\sigma}^2=3\), a candidate has \(p=4\) and \(\operatorname{RSS}=41\). Calculate the plug-in Cp to two decimal places and its adjusted prediction error in squared units. (Answer: \(C_p=1.67\); adjusted prediction error \(=5\).)
  4. Two models use the same \(n=16\) observations and a full-model variance estimate of 2. Model A has \(p=2\) and \(\operatorname{RSS}=30\); model B has \(p=4\) and \(\operatorname{RSS}=25\). Calculate both Cp values and state which has the lower adjusted prediction-error estimate. (Answer: \(C_{p,A}=3\), \(C_{p,B}=4.5\); model A is lower.)
  5. At fixed design points, the error variance is \(\sigma^2=2\). Candidate A has \(p=2\) and squared omitted-mean component 12; candidate B has \(p=3\) and squared omitted-mean component 2. Find each expected mean-response prediction error and identify the lower one. (Answer: A has \(12+2(2)=16\), B has \(2+3(2)=8\); model B is lower.)
Clinical Biostats