## Best Subsets Regression Essentials in R

The **best subsets regression** is a model selection approach that consists of testing all possible combination of the predictor variables, and then selecting the best model according to some statistical criteria.

In this chapter, we’ll describe how to compute best subsets regression using R.

Contents:

## Loading required R packages

`tidyverse`

for easy data manipulation and visualization`caret`

for easy machine learning workflow`leaps`

, for computing best subsets regression

```
library(tidyverse)
library(caret)
library(leaps)
```

## Example of data

We’ll use the built-in R `swiss`

data, introduced in the Chapter @ref(regression-analysis), for predicting fertility score on the basis of socio-economic indicators.

```
# Load the data
data("swiss")
# Inspect the data
sample_n(swiss, 3)
```

## Computing best subsets regression

The R function `regsubsets()`

[`leaps`

package] can be used to identify different best models of different sizes.
You need to specify the option `nvmax`

, which represents the maximum number of predictors to incorporate in the model. For example, if `nvmax = 5`

, the function will return up to the best 5-variables model, that is, it returns the best 1-variable model, the best 2-variables model, …, the best 5-variables models.

In our example, we have only 5 predictor variables in the data. So, we’ll use `nvmax = 5`

.

```
models <- regsubsets(Fertility~., data = swiss, nvmax = 5)
summary(models)
```

```
## Subset selection object
## Call: st_build()
## 5 Variables (and intercept)
## Forced in Forced out
## Agriculture FALSE FALSE
## Examination FALSE FALSE
## Education FALSE FALSE
## Catholic FALSE FALSE
## Infant.Mortality FALSE FALSE
## 1 subsets of each size up to 5
## Selection Algorithm: exhaustive
## Agriculture Examination Education Catholic Infant.Mortality
## 1 ( 1 ) " " " " "*" " " " "
## 2 ( 1 ) " " " " "*" "*" " "
## 3 ( 1 ) " " " " "*" "*" "*"
## 4 ( 1 ) "*" " " "*" "*" "*"
## 5 ( 1 ) "*" "*" "*" "*" "*"
```

The function `summary()`

reports the best set of variables for each model size. From the output above, an asterisk specifies that a given variable is included in the corresponding model.

For example, it can be seen that the best 2-variables model contains only Education and Catholic variables (`Fertility ~ Education + Catholic`

). The best three-variable model is (`Fertility ~ Education + Catholic + Infant.mortality`

), and so forth.

A natural question is: which of these best models should we finally choose for our predictive analytics?

## Choosing the optimal model

To answer to this questions, you need some statistical metrics or strategies to compare the overall performance of the models and to choose the best one. You need to estimate the prediction error of each model and to select the one with the lower prediction error.

### Model selection criteria: Adjusted R2, Cp and BIC

The `summary()`

function returns some metrics - Adjusted R2, Cp and BIC (see Chapter @ref(regression-model-accuracy-metrics)) - allowing us to identify the best overall model, where best is defined as the model that maximize the adjusted R2 and minimize the prediction error (RSS, cp and BIC).

The adjusted R2 represents the proportion of variation, in the outcome, that are explained by the variation in predictors values. the higher the adjusted R2, the better the model.

The best model, according to each of these metrics, can be extracted as follow:

```
res.sum <- summary(models)
data.frame(
Adj.R2 = which.max(res.sum$adjr2),
CP = which.min(res.sum$cp),
BIC = which.min(res.sum$bic)
)
```

```
## Adj.R2 CP BIC
## 1 5 4 4
```

There is no single correct solution to model selection, each of these criteria will lead to slightly different models. Remember that, “All models are wrong, some models are useful”.

Here, adjusted R2 tells us that the best model is the one with all the 5 predictor variables. However, using the BIC and Cp criteria, we should go for the model with 4 variables.

So, we have different “best” models depending on which metrics we consider. We need additional strategies.

Note also that the adjusted R2, BIC and Cp are calculated on the training data that have been used to fit the model. This means that, the model selection, using these metrics, is possibly subject to overfitting and may not perform as well when applied to new data.

A more rigorous approach is to select a models based on the prediction error computed on a new test data using k-fold cross-validation techniques (Chapter @ref(cross-validation)).

### K-fold cross-validation

The **k-fold Cross-validation** consists of first dividing the data into k subsets, also known as k-fold, where k is generally set to 5 or 10. Each subset (10%) serves successively as test data set and the remaining subset (90%) as training data. The average cross-validation error is computed as the model prediction error.

The k-fold cross-validation can be easily computed using the function `train()`

[`caret`

package] (Chapter @ref(cross-validation)).

Here, we’ll follow the procedure below:

- Extract the different model formulas from the
`models`

object - Train a linear model on the formula using k-fold cross-validation (with k= 5) and compute the prediction error of each model

We start by defining two helper functions:

`get_model_formula()`

, allowing to access easily the formula of the models returned by the function`regsubsets()`

. Copy and paste the following code in your R console:

```
# id: model id
# object: regsubsets object
# data: data used to fit regsubsets
# outcome: outcome variable
get_model_formula <- function(id, object, outcome){
# get models data
models <- summary(object)$which[id,-1]
# Get outcome variable
#form <- as.formula(object$call[[2]])
#outcome <- all.vars(form)[1]
# Get model predictors
predictors <- names(which(models == TRUE))
predictors <- paste(predictors, collapse = "+")
# Build model formula
as.formula(paste0(outcome, "~", predictors))
}
```

For example to have the best 3-variable model formula, type this:

`get_model_formula(3, models, "Fertility")`

```
## Fertility ~ Education + Catholic + Infant.Mortality
##
```

`get_cv_error()`

, to get the cross-validation (CV) error for a given model:

```
get_cv_error <- function(model.formula, data){
set.seed(1)
train.control <- trainControl(method = "cv", number = 5)
cv <- train(model.formula, data = data, method = "lm",
trControl = train.control)
cv$results$RMSE
}
```

Finally, use the above defined helper functions to compute the prediction error of the different best models returned by the `regsubsets()`

function:

```
# Compute cross-validation error
model.ids <- 1:5
cv.errors <- map(model.ids, get_model_formula, models, "Fertility") %>%
map(get_cv_error, data = swiss) %>%
unlist()
cv.errors
```

`## [1] 9.42 8.45 7.93 7.68 7.92`

```
# Select the model that minimize the CV error
which.min(cv.errors)
```

`## [1] 4`

It can be seen that the model with 4 variables is the best model. It has the lower prediction error. The regression coefficients of this model can be extracted as follow:

`coef(models, 4)`

```
## (Intercept) Agriculture Education Catholic
## 62.101 -0.155 -0.980 0.125
## Infant.Mortality
## 1.078
```

## Discussion

This chapter describes the best subsets regression approach for choosing the best linear regression model that explains our data.

Note that, this method is computationally expensive and becomes unfeasible for a large data set with many variables. A better alternative is provided by the **stepwise regression** method. See Chapter @ref(stepwise-regression).