Best subset selection

Best subset selection searches through different models by simply considering every possible combination of the available variables. We then choose the model which optimizes some performance criterion. We’ll start with \(R^2_{adj}\) as our model performance criterion, so we want to choose the combination of explanatory variables that maximizes \(R^2_{adj}\) when predicting site EUI.

Best subset selection can be run with the regsubsets function in the leaps package. The following code runs best subset selection over models with up to 3 explanatory variables (nvmax = 3).

model_select <- regsubsets(site_eui ~ ., data = large_residential, 
                           nvmax = 3, method="exhaustive")

Once we’ve run best subset selection, we can look at the results.

• To get the number of variables in the model which maximizes \(R^2_{adj}\), run

which.max(summary(model_select)$adjr2)

• To plot \(R^2_{adj}\) against the number of parameters in the model, run

plot(summary(model_select)$adjr2)

• To see which variables are included in the model which maximizes \(R^2_{adj}\), run

summary(model_select)$which[which.max(summary(model_select)$adjr2),]

• Best subset selection considers all models with up to nvmax variables. To visualize which variables are included in the best model of each size, run

plot(model_select, scale="adjr2")

In this plot, \(R^2_{adj}\) is shown on the y-axis, and each variable included at that \(R^2_{adj}\) is shaded black.

If our client is specifically interested in building age, we need to make sure it is included in the model. We can ensure a variable is included with the force.in argument in the regsubsets function. For example, if age is the 2nd explanatory variable listed in my dataset, then I would do something like

model_select <- regsubsets(site_eui ~ ., data = large_residential, 
                           nvmax = 10, force.in = 2, method="exhaustive")

Different optimality criteria

Instead of maximizing \(R^2_{adj}\), we can optimize other performance criteria. Two common metrics are Mallows’ \(C_p\) and the Bayesian Information Criterion (BIC). Like \(R^2_{adj}\), both metrics penalize the sum of squared residuals to account for the number of parameters in the model. The difference is in how the penalty term is defined.

Mallows’ \(C_p\): If we have \(p\) parameters in the model (including the intercept) and \(n\) observations, then

\[C_p = \dfrac{SSE_p}{MSE} - n + 2p\]

where \(SSE_p\) is the sum of squared residuals for the model with \(p\) parameters, and \(MSE\) is the mean squared error for the model with all possible variables.

BIC: If we have \(p\) parameters in the model (including the intercept) and \(n\) observations, then

\[BIC = n \log (SSE_p/n) + p \log(n)\]

Instead of maximizing \(R^2_{adj}\), we could choose to minimize \(C_p\) or BIC. For example, the number of parameters in the model which minimizes \(C_p\) is which.min(summary(model_select)$cp), while for minimizing BIC it is which.min(summary(model_select)$bic).

Stepwise selection

Considering all possible variable combinations is computationally expensive when we have many observations and variables. An alternative to best subset selection is stepwise selection, in which we consider adding or removing one variable at a time. Backward stepwise selection starts with a full model and removes variables, whereas forward stepwise selection starts with an empty model and adds variables.

In R, we can run forward stepwise selection with

model_select <- regsubsets(site_eui ~ ., data = large_residential, 
                           nvmax = 10, method="forward")

and backward stepwise selection with

model_select <- regsubsets(site_eui ~ ., data = large_residential, 
                           nvmax = 10, method="backward")

Citation