Table of Contents

Subset selection in python¶

This notebook explores common methods for performing subset selection on a regression model, namely

• Best subset selection
• Forward stepwise selection
• Criteria for choosing the optimal model
  • $C_p$, AIC, BIC, $R^2_{adj}$

The figures, formula and explanation are taken from the book "Introduction to Statistical Learning (ISLR)" Chapter 6 and have been adapted in python

Sources:¶

• Introduction to Statistical Learning (ISLR) - Chapter 6
• https://github.com/JWarmenhoven/ISLR-python
• http://www.science.smith.edu/~jcrouser/SDS293/labs/lab8-py.html

Libraries¶

The dataset¶

The credit dataset is a use case for linear regression where some predictors are qualitative.

Note - all datasets from the book are available here

Encoding qualitative data¶

Best subset selection¶

To perform best selection, we fit separate models for each possible combination of the $n$ predictors and then select the best subset. That is we fit:

• All models that contains exactly one predictor
• All models that contain 2 predictors at the second step: $\binom{n}{2}$
• Until reaching the end point where all $n$ predictors are included in the model

This results in $2^n$ possibilities as this is a power set problem. In our case there are $2^{11} = 2048$ possible combinations

Algorithm¶

• Let $\mathcal{M}_0$ denote the null model which contains no predictors, this model simply predicts the sample mean of each observation

  • For $k = 1,2,...,n$

    • Fit all $\binom{n}{k}$ models that contain exactly $k$ predictors
    • Pick the best among these $\binom{n}{k}$ models, and call it $\mathcal{M}_k$. Here the best is defined as having the smallest RSS, or an equivalent measure

  • Select the single best model among $\mathcal{M_0, M_1,...,M_n}$ using cross validated predicton error, $C_p$, BIC, adjusted $R^2$ or any other method.

Helper function for fitting linear regression (Sklearn)¶

Implementing Best subset selection (using itertools.combinations)¶

Finding the best subsets for each number of features¶

Using the smallest RSS value, or the largest R_squared value

Adding columns to the dataframe with RSS and R squared values of the best subset¶

Plotting the best subset selection process¶

Forward stepwise selection¶

For computational reasons, the best subset cannot be applied for any large $n$ due to the $2^n$ complexity. Forward Stepwise begins with a model containing no predictors, and then adds predictors to the model, one at the time. At each step, the variable that gives the greatest additional improvement to the fit is added to the model.

Algorithm¶

Let $\mathcal{M}_0$ denote the null model which contains no predictors

• For $k = 1,2,...,n-1$
  • Consider all $n - k$ models that augment the predictors in $\mathcal{M_k}$ with one additional predictor
  • Choose the \textit{best} among these $n - k$ models, and call it $\mathcal{M_{k+1}}$
• Select the single best model among $\mathcal{M_0, M_1,...,M_n}$ using cross validated predicton error, $C_p$, BIC, adjusted $R^2$ or any other method.

Displaying results of the first 4 steps¶

Comparing models: AIC, BIC, Mallows'CP¶

The training set Mean Squared Error (MSE) is generally an underestimate of the test MSE. This is because when we fit a model to the training data using least squares, we specifically estimate the regression coefficients such that the training RSS is minimized. In particular, the training RSS decreases as we add more features to the model, but the test error may not. Therefore the training RSS and $R^2$ may not be used for selecting the best model unless we adjust for this underestimation.

Mallow's $C_p$¶

Mallow's $C_p$ is named after Colin Lingwood Mallows and is defined as:

$$ C_p = \frac{1}{m} (RSS + 2d\hat\sigma^2)$$

where $\hat\sigma^2$ is an estimate of the variance of the error $\epsilon$ associated with each response measurement. Typically $\hat\sigma^2$ is estimated using the full model containing all predictors.

Essentially,the $C_p$ statistic adds a penalty of $2d\hat\sigma^2$ to the training RSS in order to adjust for the fact that the training error tends to underestimate the test error. Clearly, the penalty increases as the number of predictors in the model increases, and this is intended to adjust for the corresponding decrease in training RSS. It can be shown that if $\hat\sigma^2$ is an unbiased estimate of $\sigma^2$ then $C_p$ is an unbiased estimate of the test MSE, so we choose the model with the smallest $C_p$.

Akaike's Information Criteria (AIC)¶

The AIC criterion is defiend for a large class of models fit by maximum likelihood. In the case of a linear model with Gaussian errors, MLE and least squares are the same thing and the AIC is given by

$$ AIC = \frac{1}{m\hat \sigma^2} (RSS + 2d\hat\sigma^2)$$

Bayesian Information Criteria (BIC)¶

BIC is derived from a Bayesian point of view, and looks similar to the $C_p$ and $AIC$ - it is defined (up to irrelevant constants) as:

$$ BIC = \frac{1}{m\hat \sigma^2} (RSS + \log(m) d \hat\sigma^2)$$

Like $C_p$ and AIC, the BIC will tend to take small values for a model with low test error.

Adjusted $R^2$¶

Since the $R^2$ always increases as more variables are added, the adjusted $R^2$ accounts for that fact and introduces a penalty. The intuition is that once all the correct variables have been included in the model,additional noise variables will lead to a very small decrase in RSS, but an increase in $k$ and hence will decrease the adjusted $R^2$. In effect, we pay a price for the inclusion of unnecessary variables in the model.

$$ R_{a}^2 = 1 - \frac{RSS / (m - k -1)}{TSS / (m - 1)} = 1 - (1 - R^2) \frac{m - 1}{m - k - 1} $$

 Theoretical justification¶

$C_p$, AIC, BIC all have rigorous theoretical justification that rely on asymptotic arguments, i.e. when the sample size $m$ grows very large, whereas the adjusted $R^2$, although quite intuitive, is not as well motivated in statistical theory.

Estimation of $\hat\sigma^2$¶

Using the RSS of the full model with $p$ features, (i.e.the smallest RSS) we estimate $\hat\sigma^2$ as:

$$ \hat\sigma^2 = \frac{RSS}{m - p -1} $$

Combining forward stepwise results into a new DataFrame¶

Computing the C_p, AIC, BIC and R-square adjusted¶

Plotting the computed values as a function of number of features¶

Miscellaneous¶

Pairplot¶

Statsmodel linear regression¶

Least squares coefficient estimates associated with the regression of balance onto ethnicity in the Credit data set. The linear model is given in (3.30). That is, ethnicity is encoded via two dummy variables

Table 3.8 - page 86¶