Mallows's Cp
Mallows's Cp
Main page
M

Mallows's Cp

logo
Community Hub0 subscribers

Wikipedia

Mallows's Cp
View on Wikipedia
from Wikipedia

In statistics, Mallows's ,[1][2] named for Colin Lingwood Mallows, is used to assess the fit of a regression model that has been estimated using ordinary least squares. It is applied in the context of model selection, where a number of predictor variables are available for predicting some outcome, and the goal is to find the best model involving a subset of these predictors. A small value of means that the model is relatively precise.

Mallows's is 'essentially equivalent'[3] to the Akaike information criterion in the case of linear regression. This equivalence is only asymptotic; Akaike[4] notes that requires some subjective judgment in the choice of the variance estimate associated with each response in the linear model (typically denoted as ).

Definition and properties

[edit]

Mallows's addresses the issue of overfitting, in which model selection statistics such as the residual sum of squares always get smaller as more variables are added to a model. Thus, if we aim to select the model giving the smallest residual sum of squares, the model including all variables would always be selected. Instead, the statistic calculated on a sample of data estimates the sum squared prediction error (SSPE) as its population target

where is the fitted value from the regression model for the ith case, E(Yi | Xi) is the expected value for the ith case, and is the error variance (assumed constant across the cases). The mean squared prediction error (MSPE) will not automatically get smaller as more variables are added. The optimum model under this criterion is a compromise influenced by the sample size, the effect sizes of the different predictors, and the degree of collinearity between them.

If p regressors are selected from a set of k regressors, with k > p, the statistic for that particular set of regressors is defined as:

where

  • is the error sum of squares for the model with p regressors,
  • Ypi is the predicted value of the ith observation of Y from the p regressors,
  • S2 is the estimation of residuals variance after regression on the complete set of k regressors and can be estimated by ,[1]
  • and N is the sample size.

Alternative definition

[edit]

Given a linear model such as:

where:

  • are coefficients for predictor variables
  • represents error

An alternate version of can also be defined as:[5]

where

  • RSS is the residual sum of squares on a training set of data
  • p is the number of predictors
  • and refers to an estimate of the variance associated with each response in the linear model (estimated on a model containing all predictors)

Note that this version of the does not give equivalent values to the earlier version, but the model with the smallest from this definition will also be the same model with the smallest from the earlier definition.

Limitations

[edit]

The criterion suffers from two main limitations[6]

  1. the approximation is only valid for large sample size;
  2. the ' cannot handle complex collections of models as in the variable selection (or feature selection) problem.[6]

Practical use

[edit]

The statistic is often used as a stopping rule for various forms of stepwise regression. Mallows proposed the statistic as a criterion for selecting among many alternative subset regressions. Under a model not suffering from appreciable lack of fit (bias), has expectation nearly equal to p; otherwise the expectation is roughly P plus a positive bias term. Nevertheless, even though it has expectation greater than or equal to p, there is nothing to prevent Cp < p or even in extreme cases. It is suggested that one should choose a subset that has approaching p,[7] from above, for a list of subsets ordered by increasing p. In practice, the positive bias can be adjusted for by selecting a model from the ordered list of subsets, such that .

Since the sample-based statistic is an estimate of the MSPE, using for model selection does not completely guard against overfitting. For instance, it is possible that the selected model will be one in which the sample was a particularly severe underestimate of the MSPE.

Model selection statistics such as are generally not used blindly, but rather information about the field of application, the intended use of the model, and any known biases in the data are taken into account in the process of model selection.

See also

[edit]

References

[edit]

Further reading

[edit]

Grokipedia

Mallows's Cp

View on Grokipedia
from Grokipedia
Mallows's C_p is a statistical criterion developed for selecting optimal subset models in multiple linear regression analysis, providing an unbiased estimate of the expected mean squared prediction error relative to the true error variance. Named after the statistician Colin Lingwood Mallows (1930–2023), it balances the trade-off between model bias (from excluding relevant predictors) and variance (from including irrelevant ones), facilitating the identification of parsimonious models with good predictive performance.[1] The statistic was first proposed by Mallows in 1964 during a presentation on variable selection in linear regression and was formally detailed in his 1973 paper, where it was positioned as a tool for assessing model adequacy beyond traditional fit measures like the coefficient of determination.[1] In the paper, Mallows emphasized its role in subset regression, where numerous potential predictor combinations are evaluated to avoid overfitting while minimizing prediction error.[1] The formula for C_p is given by
Cp=RSSpσ^2+2p−n, C_p = \frac{RSS_p}{\hat{\sigma}^2} + 2p - n,
where RSSpRSS_p is the residual sum of squares for the subset model with pp parameters (including the intercept), nn is the number of observations, and σ^2\hat{\sigma}^2 is an unbiased estimate of the error variance, typically derived from the full model's mean squared error.[1] This expression approximates the scaled total mean squared error, adjusting for degrees of freedom to account for model complexity.[2] Interpretation of C_p relies on its value relative to pp: models where Cp≈pC_p \approx p indicate no substantial bias and are deemed adequate for prediction, as the estimated prediction error matches the model's parametric complexity.[1] Values of CpC_p significantly larger than pp suggest bias from model underspecification, while values below pp (though possible due to sampling variability) imply no bias but potential overfitting risks.[2] Mallows introduced C_p-plots, graphing the statistic against pp or the number of predictors, to visually identify stable models with low C_p values near the line Cp=pC_p = p.[1] In practice, C_p is widely applied in best subsets regression algorithms, where it ranks candidate models generated by exhaustive or stepwise searches over predictor combinations, prioritizing those that minimize prediction error without excessive parameters.[2] It offers advantages over metrics like adjusted R2R^2, which focus on explained variance but ignore bias in predictions, by directly estimating the impact of omitted variables on future observations.[2] Extensions of C_p have been developed for biased regression techniques, such as ridge regression, and robust variants to handle outliers, broadening its utility in modern statistical modeling.[1] Despite its strengths, C_p assumes ordinary least squares estimation and known error variance structure, limiting its direct applicability to nonlinear or high-dimensional settings without adaptation.[2]

Background

Linear Regression Context

In ordinary least squares (OLS) regression, the model posits a linear relationship between a response vector Y∈RnY \in \mathbb{R}^n and a design matrix X∈Rn×pX \in \mathbb{R}^{n \times p} of predictors, expressed as Y=Xβ+ϵY = X\beta + \epsilon, where β∈Rp\beta \in \mathbb{R}^p is the vector of unknown coefficients and ϵ∈Rn\epsilon \in \mathbb{R}^n is the error term typically assumed to follow a multivariate normal distribution ϵ∼N(0,σ2In)\epsilon \sim N(0, \sigma^2 I_n) with mean zero, constant variance σ2\sigma^2, and identity covariance matrix InI_n.[3] The OLS estimator minimizes the sum of squared residuals to obtain β^=(XTX)−1XTY\hat{\beta} = (X^T X)^{-1} X^T Y, providing unbiased and efficient estimates under the model's assumptions.[3] The validity of OLS inferences relies on four key assumptions: linearity, which requires the conditional expectation E(Y∣X)=XβE(Y|X) = X\beta to hold; independence, ensuring observations are independent such that residuals are uncorrelated; homoscedasticity, meaning the error variance Var(ϵi)=σ2\text{Var}(\epsilon_i) = \sigma^2 is constant across all levels of the predictors; and normality, where errors are normally distributed for exact finite-sample inference, though this is less critical for large samples due to the central limit theorem.[4] Violations of these assumptions, such as nonlinearity or heteroscedasticity, can lead to biased estimates, inefficient predictions, or invalid hypothesis tests.[5] In multiple linear regression, incorporating numerous predictors allows for capturing complex relationships but introduces the risk of overfitting, where the full model including all pp predictors fits the training data excessively well at the expense of generalizability.[6] Submodels, which select a subset of predictors, aim to balance fit and parsimony, reducing variance and mitigating overfitting by excluding irrelevant variables that inflate estimation error on new data.[6] This tradeoff is particularly acute when pp approaches or exceeds the sample size nn, as the full model may capture noise rather than signal.[6] The development of subset selection methods in linear regression gained prominence in the 1960s and 1970s, coinciding with advances in computational power that made exhaustive searches feasible beyond manual calculations.[7] Prior to the 1970s, computing all 2p2^p subsets was prohibitively expensive, limiting analyses to simple or stepwise approaches, but improved algorithms and hardware enabled more systematic exploration of model spaces.[7] This era marked a shift toward automated variable selection to handle growing datasets in fields like econometrics and environmental modeling.

Model Selection Challenges

In multiple linear regression, selecting an appropriate subset of predictor variables from a larger set poses significant challenges, as the goal is to identify a model that accurately captures the underlying relationships in the data while avoiding excessive complexity or oversimplification. This process is complicated by the need to balance the desire for a good fit to the observed data with the risk of poor performance on new, unseen data, particularly under the standard assumptions of linearity, independence, and homoscedasticity in residuals.[8] One primary issue is overfitting, where including too many predictors results in a model that fits the training data excessively well, achieving low residual error on the sample but failing to generalize to independent data due to capturing noise rather than true signal. Conversely, underfitting arises when relevant predictors are excluded, leading to a model that is too simplistic and exhibits high bias by systematically missing key patterns in the data, thereby producing consistently poor predictions across both training and test sets.[8] These risks highlight the delicate nature of subset selection, where the choice of predictors directly impacts the model's reliability for inference and prediction. At the heart of these challenges lies the bias-variance tradeoff, a fundamental principle in model selection that quantifies how increasing model complexity reduces bias (systematic error from model misspecification) but increases variance (sensitivity to fluctuations in the training data), often leading to an optimal complexity level that minimizes overall prediction error. In linear regression contexts, this tradeoff necessitates criteria that penalize excessive parameters to prevent variance-dominated models, ensuring the selected subset provides a stable estimate of the true regression surface without unnecessary variability.[9] A further complication is the computational explosion inherent in exhaustive subset selection, as evaluating all possible combinations of p predictors requires fitting 2p distinct models, rendering the approach infeasible for even moderately large p (e.g., beyond 20-30 variables) due to the exponential growth in time and resources.[10] This scalability issue has historically driven the development of heuristic alternatives, though it underscores the need for efficient, principled methods to navigate the vast search space without exhaustive enumeration. Early approaches to address these problems, such as stepwise regression (including forward, backward, and bidirectional variants), attempt to approximate optimal subsets through greedy algorithms that iteratively add or remove predictors based on statistical criteria like p-values or information measures.[11] However, these methods suffer from inherent limitations, including their greedy nature, which commits to local decisions without reconsidering global optimality, potentially leading to suboptimal models, biased parameter estimates, and inconsistencies across different algorithmic implementations.[11] Additionally, stepwise procedures can inflate Type I errors by repeatedly testing variables on the same dataset, exacerbating overfitting risks in high-dimensional settings.[12]

Mathematical Formulation

Original Definition

Mallows's CpC_p statistic was introduced by Colin L. Mallows in 1973 as a criterion for selecting subset regression models in linear regression. The original formulation defines CpC_p for a submodel with pp parameters (including the intercept) as
Cp=RSSps2−n+2p, C_p = \frac{\mathrm{RSS}_p}{s^2} - n + 2p,
where RSSp\mathrm{RSS}_p is the residual sum of squares for the submodel, s2s^2 is an unbiased estimate of the error variance σ2\sigma^2 obtained from the full model, and nn is the sample size.[1] This statistic serves as an estimate of the expected prediction error, denoted E(Jp)E(J_p), which measures the scaled sum of squared prediction errors for the submodel. Under ideal conditions where the submodel is unbiased relative to the true model, the expected value E(Cp)E(C_p) approximates pp, providing a baseline for model comparison.[1] Mallows's motivation for developing CpC_p centered on correcting the bias in estimates of the mean squared error of prediction (MSEP), particularly when comparing submodels of varying complexity to avoid overfitting. By incorporating a penalty term 2p2p, the criterion balances goodness-of-fit (via RSSp\mathrm{RSS}_p) with the increased variance from additional parameters, facilitating graphical analysis through CpC_p plots. The concept originated in Mallows's 1973 paper "Some Comments on CpC_p" published in Technometrics.[1]

Key Properties and Derivation

Mallows's CpC_p statistic is derived as an estimator of the mean squared error of prediction (MSEP) for a linear regression model with pp parameters, where MSEP quantifies the expected prediction error for new observations. Specifically, for a new response YnewY_{\text{new}} at covariates XnewX_{\text{new}}, the MSEP is given by
MSEP=E[(Ynew−Y^p)2∣Xnew]=σ2+Bias2(Y^p)+Var(Y^p), \text{MSEP} = E[(Y_{\text{new}} - \hat{Y}_p)^2 \mid X_{\text{new}}] = \sigma^2 + \text{Bias}^2(\hat{Y}_p) + \text{Var}(\hat{Y}_p),
where σ2\sigma^2 is the irreducible error variance, Bias2(Y^p)=(E[Y^p]−μ(Xnew))2\text{Bias}^2(\hat{Y}_p) = (E[\hat{Y}_p] - \mu(X_{\text{new}}))^2 measures systematic deviation from the true mean μ(Xnew)\mu(X_{\text{new}}), and Var(Y^p)\text{Var}(\hat{Y}_p) captures the variability of the predictor. To estimate this in practice, CpC_p adjusts the in-sample residual sum of squares (RSSp_p) to account for both the optimism of training error (which underestimates test error) and model complexity. Under the assumption of a linear model Y=Xβ+ϵY = X\beta + \epsilon with ϵ∼N(0,σ2In)\epsilon \sim N(0, \sigma^2 I_n), the expected sum of squared prediction errors over the training design points, E[∑i=1n(Yi−Y^p,i)2]E[\sum_{i=1}^n (Y_i - \hat{Y}_{p,i})^2], equals nσ2+∑i=1nBias2(Y^p,i)+σ2pn\sigma^2 + \sum_{i=1}^n \text{Bias}^2(\hat{Y}_{p,i}) + \sigma^2 p, while E[RSSp]=σ2(n−p)E[\text{RSS}_p] = \sigma^2 (n - p) for an unbiased model (zero bias). The difference arises from the degrees-of-freedom adjustment, leading to the optimism term 2σ2p2\sigma^2 p. Thus, RSSp+2σ^2p\text{RSS}_p + 2\hat{\sigma}^2 p unbiasedly estimates the expected prediction error sum when the model is correct and σ^2\hat{\sigma}^2 is an unbiased estimator of σ2\sigma^2 (typically from the full model). Scaling by σ^2\hat{\sigma}^2 and recentering yields CpC_p as an estimator of the scaled prediction risk minus nn, specifically E[Cp]=p+1σ2∑Bias2E[C_p] = p + \frac{1}{\sigma^2} \sum \text{Bias}^2 for correct submodels, providing an unbiased assessment of bias and variance contributions. The unbiasedness property of CpC_p holds in the sense that, for submodels nested within the true model (zero bias), E[Cp]=pE[C_p] = p, indicating that CpC_p centers around the number of parameters without systematic over- or underestimation of the scaled risk. This follows directly from the expectation E[RSSp/σ2]=n−pE[\text{RSS}_p / \sigma^2] = n - p, so substituting into Cp=RSSp/σ^2+2p−nC_p = \text{RSS}_p / \hat{\sigma}^2 + 2p - n (with σ^2\hat{\sigma}^2 unbiased) yields the result. For biased submodels (omitting relevant predictors), E[Cp]>pE[C_p] > p, as the positive bias term inflates the estimate, penalizing underfitting. This property makes CpC_p a reliable gauge for model adequacy relative to complexity, outperforming raw RSS by correcting for the tendency of training error to decrease with added parameters. Asymptotically, as the sample size n→∞n \to \infty with fixed pp and correct model specification (zero bias), CpC_p converges in probability to pp, and the normalized estimator Cp/nC_p / n converges to the scaled MSEP limit of 1, reflecting that prediction variance diminishes as O(p/n)→0\mathcal{O}(p/n) \to 0. Under misspecification, CpC_p consistently estimates the inflated risk due to persistent bias, ensuring that suboptimal models are identifiable in large samples. This behavior supports CpC_p's use in selection by favoring models where the minimum Cp≈pC_p \approx p, signaling negligible bias and optimal bias-variance trade-off.

Alternative Expressions

Residual Sum of Squares Form

An equivalent formulation of Mallows's CpC_p expresses the criterion directly in terms of the residual sum of squares (RSS) from the candidate subset model and the full model, eliminating the need for an explicit estimate of the error variance σ2\sigma^2 beyond what is inherent in the RSS values. Specifically, for a subset model with pp predictors (totaling p+1p+1 parameters including the intercept), the RSS form is given by
Cp=(n−kfull)⋅RSSpRSSfull+2(p+1)−n, C_p = (n - k_\mathrm{full}) \cdot \frac{\mathrm{RSS}_p}{\mathrm{RSS}_\mathrm{full}} + 2(p + 1) - n,
where nn is the sample size, RSSp\mathrm{RSS}_p is the residual sum of squares for the subset model, RSSfull\mathrm{RSS}_\mathrm{full} is the RSS for the full model encompassing all candidate predictors, and kfullk_\mathrm{full} is the number of parameters (including the intercept) in the full model.[13][2] This form is advantageous because it is self-contained, relying solely on observable quantities from the fitted models without requiring a separate unbiased estimate of σ2\sigma^2, provided the full model is assumed to be unbiased for the true error variance.[2] In contexts where the full model's RSS serves as a reliable proxy for scaled variance, this avoids additional computational steps or assumptions about external variance estimates.[13] The equivalence to the original definition of Cp=RSSpσ^2+2(p+1)−nC_p = \frac{\mathrm{RSS}_p}{\hat{\sigma}^2} + 2(p+1) - n follows from substituting the plug-in estimate σ^2=RSSfulln−kfull\hat{\sigma}^2 = \frac{\mathrm{RSS}_\mathrm{full}}{n - k_\mathrm{full}}.[13][2] In finite samples, this RSS-based approach can reduce variability in the CpC_p estimate compared to using a potentially noisy standalone σ^2\hat{\sigma}^2, as the full model's RSS leverages the assumption of minimal bias to stabilize the variance proxy, particularly when the full model includes sufficient predictors to capture the underlying structure.[2]

Relationship to Prediction Error

Mallows's CpC_p criterion is closely connected to the predicted residual sum of squares (PRESS) statistic, which serves as a direct measure of a model's predictive accuracy through leave-one-out cross-validation. The PRESS is defined as the sum of squared prediction errors obtained by omitting each observation in turn and predicting it using the model fitted to the remaining data:
PRESSp=∑i=1n(yi−y^(−i))2, \text{PRESS}_p = \sum_{i=1}^n (y_i - \hat{y}_{(-i)})^2,
where y^(−i)\hat{y}_{(-i)} denotes the predicted value for the ii-th observation based on the model excluding that observation. Equivalently, PRESS can be computed using the residuals ei=yi−y^ie_i = y_i - \hat{y}_i and the diagonal elements hiih_{ii} of the hat matrix H=X(XTX)−1XTH = X(X^T X)^{-1} X^T:
PRESSp=∑i=1n(ei1−hii)2. \text{PRESS}_p = \sum_{i=1}^n \left( \frac{e_i}{1 - h_{ii}} \right)^2.
This formulation quantifies the out-of-sample prediction error, making PRESS a rigorous but computationally intensive metric for assessing model performance, especially in scenarios with many candidate subsets.[14] CpC_p functions as an efficient proxy for PRESS by approximating the scaled prediction error without the need for repeated model refitting. Specifically, under the assumption of normally distributed errors,
Cp≈nn−p−1⋅PRESSpσ2, C_p \approx \frac{n}{n - p - 1} \cdot \frac{\text{PRESS}_p}{\sigma^2},
where nn is the sample size, pp is the number of parameters (including the intercept), and σ2\sigma^2 is the error variance. This approximation adjusts PRESS for model complexity and sample size, allowing CpC_p to estimate the total expected prediction error (bias plus variance) in a single computation using the residual sum of squares from the full model. Models with CpC_p values close to p+1p + 1 indicate low bias relative to the true error variance, mirroring the goal of minimizing PRESS. The link relies on the normality assumption for exactness, as deviations can affect the bias-variance decomposition underlying both metrics.[14] In the 1970s context of subset selection in multiple regression, Mallows introduced CpC_p to contrast its computational efficiency with the full PRESS calculation, particularly beneficial when dealing with large numbers of predictors where leave-one-out methods become prohibitive. This efficiency enabled practitioners to evaluate numerous models graphically via CpC_p plots, prioritizing predictive reliability over exhaustive validation in resource-limited settings.[14]

Computation and Estimation

Step-by-Step Calculation

To compute Mallows's CpC_p statistic for model selection in linear regression, begin by fitting the full model that includes all candidate predictors using ordinary least squares (OLS). Let nn denote the number of observations, and let the full model have kk parameters (including the intercept). Compute the residual sum of squares for this full model, denoted RSSfullRSS_{\text{full}}, and estimate the error variance σ^2\hat{\sigma}^2 as σ^2=RSSfull/(n−k)\hat{\sigma}^2 = RSS_{\text{full}} / (n - k). This plug-in estimate σ^2\hat{\sigma}^2 accounts for the uncertainty in the true variance σ2\sigma^2, providing a stable denominator for subsequent calculations across submodels.[13] Next, for each candidate submodel with pp effective parameters (including the intercept), fit the model via OLS on the selected predictors while setting coefficients for excluded variables to zero. Obtain the residual sum of squares RSSpRSS_p for this submodel, where RSSp=∑i=1n(yi−y^p,i)2RSS_p = \sum_{i=1}^n (y_i - \hat{y}_{p,i})^2 and y^p,i\hat{y}_{p,i} are the fitted values. The parameter count pp reflects the dimensionality of the submodel, typically ranging from 1 (intercept only) to kk (full model).[13] Finally, substitute these values into the CpC_p formula:
Cp=RSSpσ^2+2p−n. C_p = \frac{RSS_p}{\hat{\sigma}^2} + 2p - n.
This expression adjusts for the bias introduced by estimating σ^2\hat{\sigma}^2 from the full model and penalizes model complexity via the 2p2p term; when the submodel is adequate, the expected value of CpC_p is approximately pp due to its unbiasedness property under the assumption that the true model is contained in the full model. Compute CpC_p for each submodel and compare values, often plotting them against pp to identify minima near the line Cp=pC_p = p. The adjustment in the formula, incorporating +2+2 relative to a naive 2(p−1)−n2(p-1) - n if pp excludes the intercept, ensures consistency with standard regression conventions where pp includes it.[13][15] When σ^2\hat{\sigma}^2 is estimated from the full model as S2=RSSfull/(n−kfull)S^2 = RSS_{\text{full}} / (n - k_{\text{full}}), where kfullk_{\text{full}} is the number of parameters in the full model, this approach mitigates variability in variance estimates across submodels, promoting reliable comparisons even if the full model overfits slightly. In practice, if the full model's degrees of freedom are low, sensitivity analyses may verify robustness, but the plug-in method remains the standard.[2]

Illustrative Example

Consider a hypothetical dataset with n=10n = 10 observations and 3 candidate predictors (X1,X2,X3X_1, X_2, X_3), so the full model has k=4k = 4 parameters (including intercept). Suppose fitting the full model yields RSSfull=30RSS_{\text{full}} = 30, giving σ^2=30/(10−4)=5\hat{\sigma}^2 = 30 / (10 - 4) = 5. For submodel 1 using only X1X_1 (p=2p = 2), OLS fitting gives RSSp=2=40RSS_{p=2} = 40. Then,
Cp=2=405+2(2)−10=8+4−10=2. C_{p=2} = \frac{40}{5} + 2(2) - 10 = 8 + 4 - 10 = 2.
For submodel 2 using X1X_1 and X2X_2 (p=3p = 3), RSSp=3=35RSS_{p=3} = 35, so
Cp=3=355+2(3)−10=7+6−10=3. C_{p=3} = \frac{35}{5} + 2(3) - 10 = 7 + 6 - 10 = 3.
Both submodels yield Cp=pC_p = p, suggesting reasonable fits relative to their complexity, while a submodel with Cp>pC_p > p would indicate potential bias.

Software Implementation

Mallows's CpC_p is commonly implemented in statistical software for subset selection in linear regression, leveraging exhaustive or heuristic searches to evaluate models based on the criterion. In R, the leaps package provides an efficient branch-and-bound algorithm for best subset regression, computing CpC_p values across possible subsets.[16] For example, the regsubsets function can be used as follows: 
library(leaps)
model <- regsubsets(x = predictors, y = response, nvmax = 10, method = "exhaustive")
cp_values <- summary(model)$cp
plot(cp_values, main = "Mallows's Cp by Number of Predictors")
This approach identifies the subset size where CpC_p is minimized, typically near the number of predictors plus one.[16] In Python, there is no built-in equivalent to leaps in scikit-learn or statsmodels for exhaustive subset selection with CpC_p, due to the combinatorial explosion in high dimensions; instead, custom implementations combine itertools.combinations for subset enumeration with statsmodels.OLS for fitting and CpC_p calculation.[17] A representative custom function might iterate over subsets up to a specified size, compute the residual sum of squares (RSS), and derive Cp=RSSs2+2p−nC_p = \frac{\text{RSS}}{s^2} + 2p - n, where s2s^2 is the mean squared error from the full model: 
import itertools
import statsmodels.api as sm
import numpy as np

def mallows_cp(X, y, max_features=5):
    full_model = sm.OLS(y, sm.add_constant(X)).fit()
    s2 = full_model.mse_resid
    n, p = X.shape
    best_cp = float('inf')
    for k in range(1, min(max_features + 1, p + 1)):
        for combo in itertools.combinations(range(p), k):
            X_subset = X[list(combo)]
            model = sm.OLS(y, sm.add_constant(X_subset)).fit()
            rss = model.ssr
            cp = (rss / s2) + 2 * (k + 1) - n
            if cp < best_cp:
                best_cp = cp
    return best_cp
This method is practical for small pp (e.g., <20) but scales poorly.[17] Commercial software like SAS and Stata also support CpC_p-based selection. In SAS, PROC REG with the SELECTION=CP option in the MODEL statement performs forward, backward, or stepwise selection minimizing CpC_p, outputting selected models and statistics. For instance: 
proc reg data=mydata;
   model response = var1 var2 var3 / selection=cp;
run;
In Stata, CpC_p is not directly built into regress but can be computed post-estimation using the user-contributed rsquare command, which calculates CpC_p alongside other criteria after fitting models. Post-2020 developments have enhanced workflow integration and scalability. In R's tidymodels ecosystem, the broom package tidies outputs from leaps::regsubsets, enabling seamless piping with dplyr and ggplot2 for automated model selection and visualization within machine learning pipelines.[18] For high-dimensional data where p>20p > 20, exact computation becomes infeasible, and approximations via penalized regression in the glmnet package use LASSO paths to mimic subset selection, with CpC_p analogs like generalized cross-validation (GCV) for tuning. The exact exhaustive search for CpC_p has computational complexity O(2pn3)O(2^p n^3) in the naive case, due to 2p2^p subsets and O(n3)O(n^3) least squares fits each, though optimized algorithms like branch-and-bound reduce the constant factors but retain exponential dependence on pp.[19]

Interpretation and Use

Guideline for Model Choice

Mallows's CpC_p statistic serves as a key tool for selecting among candidate regression models by estimating their predictive accuracy. The primary rule for model choice is to identify the model with the smallest CpC_p value near pp, where pp denotes the number of parameters including the intercept, as this indicates the lowest estimated mean squared prediction error with minimal bias relative to the full model. Models where Cp≈pC_p \approx p are particularly preferred, suggesting an unbiased estimate of the prediction error variance.[2] A common diagnostic approach involves plotting CpC_p against pp, the number of parameters in each candidate model. In such plots, desirable models appear near the line Cp=pC_p = p, often forming a diagonal pattern for well-separated subsets; analysts look for a flat minimum or an "elbow" where CpC_p stabilizes close to this line, signaling adequate fit without excessive complexity. If multiple models exhibit similar low CpC_p values forming distinct minima, ambiguity arises, and the parsimonious model (with fewer parameters) is typically selected to favor simplicity while maintaining comparable predictive performance.[20] Thresholds provide further guidance on bias: models with CpC_p significantly greater than pp are indicative of substantial bias, suggesting that the omitted variables contribute meaningfully to prediction error. In contrast, values of CpC_p near or below pp imply little to no bias from model underspecification.[15] CpC_p can be combined with F-tests to assess the significance of differences between candidate and full models.

Practical Examples

In a simulated dataset with 100 observations generated from the model Y = b₁ + b₂X + b₃X³ + ε (with b₁ = b₂ = b₃ = 0.1, X ~ Exp(1) - 1, and ε ~ N(0,1)), Mallows's Cp is applied to compare the constant-only model (p=1) against the model with X (p=2). The Cp value for the p=2 model is 2.28, while for p=1 it is 2.84; an adjusted Cp variant (Cp⊕) selects the p=1 model in 52% of simulations compared to 25% for standard Cp, illustrating differences in model selection performance for submodels of a nonlinear true model.[21] A classic real-world application appears in the Boston Housing dataset from the 1970s, comprising 506 observations on median house prices (MEDV) and 13 predictors including crime rate (CRIM), average number of rooms (RM), and property tax rate (TAX). Using best subsets regression on a subset of five predictors (industry proportion, rooms, age, tax, crime), Cp values are computed across model sizes: 283.9 for p=1, 77.5 for p=2, 23.3 for p=3, 4.5 for p=4, and 6.0 for the full p=5 model. The optimal model with p=4 parameters (three predictors: rooms, age, tax, crime) is selected due to its minimum Cp of 4.5, approximately equal to p=4, indicating low bias relative to the full model with higher Cp, which suggests unnecessary inclusion of the industry variable increases estimated prediction error. A typical Cp plot for this analysis visualizes the curve dipping to a minimum at p=4 before rising, guiding selection toward parsimony while retaining key factors like rooms and crime that explain housing price variance. In post-2000 econometrics applications, such as wage prediction using panel data from the National Longitudinal Survey of Youth (1979–1993) on log hourly wages for white males, Mallows's Cp informs model averaging across seven predictors including education, experience, and family background. Bayesian model averaging weighted by Cp-like criteria selects combinations emphasizing ability and education measures, yielding lower predictive loss via cross-validation compared to single models, as the approach balances bias and variance in high-dimensional economic data.[22]

Limitations and Extensions

Known Shortcomings

Mallows's CpC_p criterion assumes an unbiased estimate of the error variance σ2\sigma^2 derived from the full model encompassing all potential predictors. However, this estimate can become unstable if the full model includes many irrelevant variables, which reduces the degrees of freedom and increases variability, or if multicollinearity is severe, leading to unreliable CpC_p values and potentially suboptimal model selections. Multicollinearity exacerbates this issue by destabilizing coefficient estimates in the full model, which increases the variability of the variance estimate and undermines the reliability of CpC_p for assessing model adequacy.[23] In high-dimensional scenarios where the number of predictors pp approaches or exceeds the sample size nn, CpC_p yields unstable estimates and fails to consistently identify the true model, as its derivation relies on low-dimensional assumptions. This instability was highlighted in 1980s analyses of variable selection procedures, where small sample sizes relative to pp led to erratic model choices.[24] Additionally, CpC_p does not incorporate parameter uncertainty, treating estimates as fixed rather than accounting for variability in coefficients, in contrast to Bayesian approaches that integrate such uncertainty for more robust inference.[25] Empirical simulations from the 1990s demonstrate that under non-normal errors, such as those with outliers or heteroscedasticity, CpC_p frequently selects overly complex models, performing suboptimally compared to criteria like AIC by more frequently incorporating unnecessary predictors in such conditions.[26] This sensitivity arises because CpC_p relies on least squares residuals, which are sensitive to deviations from the assumed error structure, such as non-normality, outliers, or heteroscedasticity, leading to degraded performance when errors deviate from Gaussianity. Furthermore, CpC_p, originally developed for low-dimensional settings where p≪np \ll n, requires modifications to be effective in high-dimensional big data regimes where p≫np \gg n, rendering it ineffective without adaptations for sparse or high-dimensional structures.[27]

Modern Adaptations

In the early 2000s, extensions of Mallows's Cp were developed to accommodate generalized linear models (GLMs), where the response variable follows a distribution from the exponential family. The generalized Cp (GCp) adapts the criterion to use the deviance as a measure of fit, defined as
Cp=deviancepϕ+2p−n, C_p = \frac{\text{deviance}_p}{\phi} + 2p - n,
where deviancep\text{deviance}_p is the model's deviance, ϕ\phi is the dispersion parameter, pp is the number of parameters, and nn is the sample size. This formulation provides an unbiased estimate of prediction error in GLMs, such as logistic or Poisson regression, by scaling the deviance analogously to the residual sum of squares in linear models. Cantoni et al. (2005) introduced GCp for marginal longitudinal GLMs fitted via generalized estimating equations, demonstrating its effectiveness in variable selection under correlated data.[28] Robust variants of Cp emerged to mitigate sensitivity to outliers, replacing the least squares objective with bounded influence functions like the Huber loss. Ronchetti and Staudte (1994) proposed a robust Cp based on M-estimators for linear models, which minimizes a robust loss while penalizing complexity. In the 2010s, these ideas were extended to econometrics, where heteroskedasticity and outliers are common; for instance, Liu and Okui (2013) developed a heteroskedasticity-robust Cp for model averaging in linear regression, improving selection consistency under non-constant variance. Such adaptations maintain the predictive focus of original Cp but enhance reliability in contaminated datasets. High-dimensional settings, where the number of predictors greatly exceeds observations, prompted relaxed Cp criteria incorporating sparsity penalties like those in LASSO. These versions apply coordinate descent to approximate Cp for tuning the penalty parameter, balancing bias reduction with variable selection. Meinshausen (2007) introduced the relaxed LASSO, which performs initial selection followed by debiasing, and subsequent work since 2015 has integrated Cp-like penalties for structured high-dimensional models, such as trees or graphs, ensuring computational feasibility.[29] Bayesian analogs of Cp incorporate prior distributions to regularize model complexity, often merging with the Bayesian Information Criterion (BIC) for posterior model probabilities. BIC, proposed by Schwarz (1978), penalizes with log⁡(n)\log(n) times the number of parameters, serving as a Bayesian counterpart to Cp by favoring parsimonious models under unit information priors. This integration allows Cp-style selection in Bayesian frameworks, enhancing uncertainty quantification in sparse or hierarchical models. Recent software supports these adaptations, with R's Dyn4cast package enabling Cp computation for GLMs and non-linear models, including robust extensions. For high-dimensional sparse models, the glmnet package implements LASSO and relaxed variants, allowing Cp-based tuning via cross-validation proxies, with updates through 2023.

References

  1. Some Comments on C P - jstor
  2. 11.3 - Best Subsets Regression, Adjusted R-Sq, Mallows Cp
  3. [PDF] OLS in Matrix Form
  4. The Four Assumptions of Linear Regression - Statology
  5. Statistical notes for clinical researchers: simple linear regression 3
  6. Overfitting Regression Models: Problems, Detection, and Avoidance
  7. [PDF] The Variable Selection Problem
  8. [PDF] UVA CS 4774: Machine Learning Lecture 9: Bias-Variance Tradeoff
  9. [PDF] Linear Model Selection and Regularization
  10. [PDF] Selection of Subsets of Regression Variables
  11. Why do we still use stepwise modelling in ecology and behaviour?
  12. Step away from stepwise | Journal of Big Data | Full Text
  13. [PDF] Some Comments on C P C. L. Mallows Technometrics, Vol. 15, No ...
  14. Some Comments on C p : Technometrics - Taylor & Francis Online
  15. What is Mallows' Cp? (Defintion & Example) - Statology
  16. Mallows' Cp: Definition, Example - Statistics How To
  17. [PDF] leaps: Regression Subset Selection
  18. Choosing the optimal model: Subset selection - Data Blog
  19. Tidy a(n) regsubsets object - broom
  20. [PDF] Understanding Best Subset Selection: A Tale of Two C(omplex)ities
  21. [PDF] Chapter 9 Model Selection - rafalab
  22. [PDF] Mallows Cp for Out-of-sample Prediction
  23. [PDF] Economic variable selection arXiv:1903.02136v4 [stat.ME] 18 Mar ...
  24. The interpretation of Mallows's Cp-statistic - jstor
  25. [PDF] Comments on Mallows' c, Statistics And Multicollinearity Effects
  26. [PDF] Selecting Important Independent Variables in Linear Regression ...
  27. [PDF] Bayesian Model Averaging: A Tutorial - Colorado State University
  28. A Robust Version of Mallows's Cp - jstor
  29. Choosing a penalty for model selection in heteroscedastic regression
  30. https://academic.oup.com/biometrics/article/61/2/507/7306505
  31. Information criteria for structured parameter selection in high ... - arXiv
User Avatar
No comments yet.