Why Work Through Every Subset?
When a regression has only a few candidate predictors, an exhaustive comparison makes the trade-off between fit and model size visible. A model with more predictors cannot have a larger residual sum of squares than one of its nested submodels, because least squares can always set an added coefficient to zero. But improving fit alone is not a sufficient reason to keep every predictor. Best-subset analysis compares the fitted subsets using criteria that balance fit and complexity in different ways.
This tutorial works through a hypothetical example with three candidate predictors. The aim is not to claim that one criterion is universally best. It is to calculate each criterion correctly, identify its preferred subset, and explain why the answers can differ.
Hypothetical Data and Model Summaries
Suppose a hypothetical data set contains \(n=30\) individuals. The response \(Y\) is an exam score in points. The candidate predictors are \(A\), weekly practice hours; \(B\), a baseline quiz score in points; and \(C\), attendance in days. For this example, assume the same 30 complete observations are used in every model, each fit is ordinary least squares with an intercept, and the total sum of squares about the response mean is \(SST=50\) score-points squared.
Assume the following residual sums of squares, \(SSE\), are obtained by fitting all subsets. These summaries constitute the hypothetical analysis data; they are not results from a real study. Because there are three candidates, the full list has \(2^3=8\) models.
| Predictors included | Number of predictors \(k\) | \(SSE\) |
|---|---|---|
| None (intercept only) | 0 | 50 |
| A | 1 | 42 |
| B | 1 | 40 |
| C | 1 | 48 |
| A, B | 2 | 22 |
| A, C | 2 | 30 |
| B, C | 2 | 28 |
| A, B, C | 3 | 20 |
The reported sums of squares obey the basic nesting check: each two-predictor model has no larger \(SSE\) than either of its one-predictor submodels, and the three-predictor model has no larger \(SSE\) than any model with fewer predictors. For example, adding \(C\) to the \(A,B\) model reduces \(SSE\) from 22 to 20. This reduction improves raw fit, but it does not by itself settle whether \(C\) should be retained.
Worked Example: Calculate Fit and Complexity Criteria
Worked Example: Calculate Fit and Complexity Criteria
Calculate adjusted \(R^2\) and Mallows’ \(C_p\) for the \(A,B\) model, and calculate AIC and BIC for the \(C\)-only model.
Adjusted \(R^2\) for \(A,B\). Here \(k=2\), so the residual degrees of freedom are \(n-k-1=30-2-1=27\). The formula is
Substituting \(SSE=22\), \(SST=50\), and \(n=30\),
Mallows’ \(C_p\) for \(A,B\). Let \(p=k+1\) count all estimated regression coefficients, including the intercept. Estimate the error variance from the full model: \(MSE_{\mathrm{full}}=20/(30-3-1)=20/26\). Then
The unrounded full-model estimate \(20/26\) is used in this calculation. Rounding it first could change the final displayed result.
AIC and BIC for \(C\). For this example, use the versions with constants common to all candidate models omitted:
For \(C\), \(k=1\), so \(p=2\), including the intercept. Substitution gives
Both information criteria are minimized, so smaller values are preferred when comparing these models. Omitting terms common to all models does not change their ranking.
What the Criteria Measure
The coefficient of determination is \(R^2=1-SSE/SST\), which measures the fraction of the observed response variation accounted for by the fitted model. Since \(SST\) is fixed across these models, minimizing \(SSE\) maximizes \(R^2\). That makes raw \(R^2\) unsuitable on its own for choosing a subset: it can only stay the same or rise as predictors are added.
Adjusted \(R^2\) corrects for model size by comparing residual variation per residual degree of freedom with total variation per total degree of freedom. Its usual form for an intercept model is
A larger adjusted \(R^2\) is preferred. Adding a predictor can increase or decrease it: the predictor must reduce \(SSE\) enough to offset the loss of a residual degree of freedom.
Mallows’ \(C_p\) compares a subset model’s \(SSE\) with an estimate of the error variance from a specified larger reference model. Here that reference is the full \(A,B,C\) model, for which \(MSE_{\mathrm{full}}=20/26\). In the convention used here,
A relatively small \(C_p\) is desirable; a useful rule of thumb is to look for a value near \(p\), while also checking fit and the purpose of the analysis. The rule is not an automatic selection test. Its meaning depends on the full model providing a suitable estimate of the error variance.
AIC and BIC are likelihood-based criteria. Under the stated Gaussian regression model, their fit term reflects how well the model fits, while their penalty increases with the number of estimated parameters. AIC uses a penalty of \(2p\); BIC uses \(p\ln(n)\). Smaller values are preferred for both. Since \(\ln(30)>2\), BIC’s per-parameter penalty exceeds AIC’s in this example, so BIC can favor a smaller subset. This is a difference in their trade-offs, not a calculation error.
Worked Example: Rank All Eight Models
Worked Example: Rank All Eight Models
Calculate the criteria for each subset and identify the selected model under each rule. The table uses \(p=k+1\), \(df=30-k-1\), and \(MSE_{\mathrm{full}}=20/26\). Values are rounded for display; the rankings use the unrounded calculations.
| Subset | \(k\) | \(SSE\) | Adjusted \(R^2\) | \(C_p\) | AIC | BIC |
|---|---|---|---|---|---|---|
| None | 0 | 50 | 0.000 | 37.000 | 17.325 | 18.726 |
| A | 1 | 42 | 0.130 | 28.600 | 14.094 | 16.897 |
| B | 1 | 40 | 0.171 | 26.000 | 12.630 | 15.433 |
| C | 1 | 48 | 0.006 | 36.400 | 18.100 | 20.903 |
| A, B | 2 | 22 | 0.5274 | 4.600 | -3.305 | 0.899 |
| A, C | 2 | 30 | 0.3556 | 15.000 | 6.000 | 10.204 |
| B, C | 2 | 28 | 0.3985 | 12.400 | 3.930 | 8.134 |
| A, B, C | 3 | 20 | 0.5538 | 4.000 | -4.164 | 1.441 |
For instance, for the full model, \(p=4\) and \(df=26\). Its adjusted coefficient of determination is
Its information criteria are
The smallest \(C_p\) in the table is 4.000 for \(A,B,C\), and the largest adjusted \(R^2\) is also for \(A,B,C\). AIC is smallest for \(A,B,C\), while BIC is smallest for \(A,B\), at approximately 0.899 rather than 1.441. Thus three criteria select the full subset in this hypothetical comparison; BIC selects the two-predictor subset.
Worked Example: Explain the Change in Selected Subset
Worked Example: Explain the Change in Selected Subset
Compare the \(A,B\) and \(A,B,C\) models and explain why BIC chooses the smaller one even though the full model has lower residual error.
Adding \(C\) reduces \(SSE\) from 22 to 20. This improves the raw fit, and adjusted \(R^2\) rises from approximately 0.5274 to 0.5538. The full model’s \(C_p=4.000\) is also a little smaller than the \(A,B\) model’s \(C_p=4.600\). AIC falls from approximately \(-3.305\) to \(-4.164\), so it too favors the full model.
BIC also rewards the improved fit, but charges a larger penalty for one additional parameter. Its value rises from approximately 0.899 for \(A,B\) to 1.441 for \(A,B,C\), so BIC prefers \(A,B\) among these candidates. The criterion is not asserting that \(C\) has no relationship with the outcome. It is applying its particular balance of fit and complexity to this set of models and this sample size.
The choices are therefore criterion-dependent. A careful report would state the criteria considered and their selected subsets, rather than presenting “the best model” as though the data determined one uniquely correct answer without a decision rule.
Common Pitfalls and Practical Interpretation
- Choosing the model with the smallest \(SSE\) automatically. The full model cannot fit worse than its nested submodels, so its smaller \(SSE\) alone does not account for its added complexity. Compare the specified selection criteria.
- Using the wrong parameter count. In the formulas here, \(p\) includes the intercept. A model with two predictors has \(p=3\), not \(p=2\). The residual degrees of freedom are \(n-k-1\).
- Minimizing adjusted \(R^2\), or maximizing AIC/BIC. The preferred direction is the reverse: maximize adjusted \(R^2\), and minimize \(C_p\), AIC, and BIC under the conventions defined here.
- Comparing criteria calculated on different observations. The table’s comparisons depend on every model using the same response observations. If different models use different cases, their fit criteria are not directly comparable in this way.
- Treating a selected subset as proof of causation or certainty. A criterion ranks the candidates supplied to the search under a stated rule. It does not establish that a selected predictor causes score changes, or that an omitted predictor is irrelevant in every setting.
Best-subset analysis is particularly transparent when the number of candidates is small: every model can be fitted and inspected. Its value lies in showing how conclusions depend on the fit-complexity trade-off. For the hypothetical scores here, the \(A,B\) subset is the BIC choice, while \(A,B,C\) is preferred by adjusted \(R^2\), \(C_p\), and AIC. That difference is part of the result and should be reported plainly.
Check Your Understanding
Answer these questions using the definitions and calculations in the tutorial.
- How many subsets must be fitted when there are four candidate predictors, if the intercept-only model is included?
- For a model with \(n=30\) observations and \(k=2\) predictors plus an intercept, what are the residual degrees of freedom and the parameter count \(p\) used in these formulas?
- Why does ordinary \(R^2\) not provide a penalty for adding a predictor?
- Which subset does BIC select in the hypothetical analysis, and what feature of its trade-off helps explain that choice?
- For the \(A,B\) model, verify the adjusted \(R^2\) using \(SSE=22\) and \(SST=50\). What does the value summarize?
- Why should the final conclusion name the selection criterion rather than simply call one subset “the best model”?