With \(p\) predictors, there are \(2^p\) possible linear models. The goal of variable selection is to find the subset of predictors that best predicts \(Y\) while avoiding overfitting. Including too many variables increases variance and reduces interpretability; including too few leads to bias.
Three main approaches exist:
The most direct approach: fit all \(2^p\) possible models and choose the best one according to cross-validation error or an information criterion.
For each number of predictors \(k = 0, 1, \ldots, p\):
Then select the single best model across all sizes using CV, AIC, BIC, or adjusted \(R^2\).
Computational cost: fitting \(2^p\) models is infeasible for large \(p\). For \(p = 40\), this means over \(10^{12}\) models. Practical algorithms like branch-and-bound (implemented in the leaps R package) prune the search space but are still limited to roughly \(p \leq 40\).
Start with the null model (intercept only). At each step, add the variable whose inclusion most improves the fit (e.g., largest reduction in RSS). Stop when no variable provides a statistically significant improvement or when the chosen criterion no longer improves.
This is a greedy algorithm. It requires fitting:
\[ 1 + \sum_{k=0}^{p-1}(p - k) = 1 + \frac{p(p+1)}{2} \in O(p^2) \text{ models}\]
Much cheaper than best subset, but not guaranteed to find the globally optimal model.
Start with the full model (all \(p\) predictors). At each step, remove the variable whose removal causes the smallest increase in RSS. Continue until a stopping rule is met.
Requirement: \(n > p\) (the full model must be estimable). Same \(O(p^2)\) computational cost as forward selection.
Combines forward and backward steps: at each iteration, variables can be added or removed. This allows correction of earlier greedy decisions. Most implementations of step() in R use this approach by default.
These criteria penalize model complexity to avoid overfitting. They allow model comparison without a separate validation set.
\[ \text{AIC} = -2 \ln \hat{L} + 2k \]
where \(\hat{L}\) is the maximized likelihood and \(k\) is the number of estimated parameters. The AIC is asymptotically equivalent to leave-one-out cross-validation. Select the model with the lowest AIC.
\[ \text{BIC} = -2 \ln \hat{L} + k \ln(n) \]
Heavier penalty than AIC (since \(\ln(n) > 2\) for \(n \geq 8\)), leading to sparser models. BIC is consistent: as \(n \to \infty\), it selects the true model with probability 1 (assuming the true model is in the candidate set).
\[ R^2_{\text{adj}} = 1 - \frac{RSS / (n - p - 1)}{TSS / (n - 1)} \]
Unlike raw \(R^2\), the adjusted version penalizes additional predictors. A new variable only increases \(R^2_{\text{adj}}\) if its contribution to fit outweighs the loss of a degree of freedom. Select the model with the highest adjusted \(R^2\).
\[ C_p = \frac{RSS}{\hat{\sigma}^2} - n + 2p \]
where \(\hat{\sigma}^2\) is estimated from the full model. \(C_p\) estimates the mean squared prediction error. For an unbiased model, \(C_p \approx p\). Select the model with the lowest \(C_p\) (or where \(C_p \approx p\)).
Filter methods rank features independently of any predictive model. Each feature is scored individually against the target, and a threshold or top-\(k\) rule is applied.
Common scoring functions:
Advantages: very fast, scales to high-dimensional data, model-agnostic.
Disadvantages: ignores feature interactions and redundancies. A feature uninformative alone may be powerful in combination with others.
Wrapper methods use a specific model's performance (e.g., cross-validated accuracy) to evaluate feature subsets. The model acts as a "black box" that scores each candidate subset.
Examples:
Advantages: accounts for feature interactions; tailored to the specific model.
Disadvantages: computationally expensive (multiple model fits); risk of overfitting to the specific model choice.
Feature selection is built directly into the model training procedure.
The LASSO solves:
\[ \min_{\beta} \left\{ \frac{1}{2n} \sum_{i=1}^{n}(y_i - x_i^T \beta)^2 + \lambda \sum_{j=1}^{p} |\beta_j| \right\} \]
The \(L_1\) penalty drives coefficients exactly to zero, performing automatic variable selection. The regularization parameter \(\lambda\) controls sparsity — larger \(\lambda\) means fewer selected variables. Typically chosen by cross-validation.
Combines \(L_1\) and \(L_2\) penalties:
\[ \min_{\beta} \left\{ \frac{1}{2n} \|y - X\beta\|_2^2 + \lambda \left( \alpha \|\beta\|_1 + \frac{1-\alpha}{2}\|\beta\|_2^2 \right) \right\} \]
When predictors are correlated, LASSO tends to select one and ignore others. Elastic Net groups correlated features together. Setting \(\alpha = 1\) gives LASSO; \(\alpha = 0\) gives Ridge (no selection).
Decision trees and ensembles (Random Forest, Gradient Boosting) naturally rank features by how much they reduce impurity (Gini, entropy) or prediction error. A common approach:
| Scenario | Recommended approach |
|---|---|
| Low \(p\) (\(< 20\)) | Best subset selection or stepwise + information criteria |
| Moderate \(p\) (20–100) | Stepwise with BIC/AIC, or LASSO |
| High \(p\) (\(> 100\), possibly \(p > n\)) | LASSO, Elastic Net, or filter + wrapper pipeline |
| Correlated predictors | Elastic Net or grouped LASSO |
| Non-linear relationships | Tree-based importance or mutual information filters |
Key principles: