# Features selection for linear regression

We have seen how to fit a linear regression model, with multiple features, how we can interpret the results and how accurately we can predict.

One issue in linear regression with multiple features was which one to use. Not all available variables should be used: beside the overfitting problem, each new variable requires more data. We will see now how to appropriately choose variables.

Following is an example taken from the masterpiece book Introduction to Statistical Learning  by Hastie, Witten, Tibhirani, James. It is based on an Advertising Dataset, available on the accompanying web site or on Kaggle.

The dataset contains statistics about the sales of a product in 200 different markets, together with advertising budgets in each of these markets for different media channels: TV, radio and newspaper.

Imaging being the Marketing responsible and you need to prepare a new advertising plan for next year. You may be interested in answering questions such as:
Which media contribute to sales?
Do all three media—TV, radio, and newspaper—contribute to sales, or do just one or two of the media contribute?

Suppose that in your role are asked to suggest, on the basis of these data, a marketing plan for next year that will result in high product sales. What information would be useful in order to provide such a recommendation?

As usual, this example is also available as Jupyter notebook in Github.

Note: I will use interchangeably the words feature, variable or predictor to indicate the x_i inputs of the linear regression.

```import pandas as pd

import matplotlib.pyplot as plt

plt.show()```

# ## Is there a relationship betwen sales and advertising?

First of all, we fit a regression line using the Ordinary Least Square algorithm, i.e. the line that minimises the squared differences between the actual Sales and the line itself.

The multiple linear regression model takes the form:
Sales = β0 + β1 * TV + β2 * Radio + β3 * Newspaper + ε, where Beta are the regression coefficients we want to find and epsilon is the error that we want to minimise.
For this we use the statsmodels package and its ols() function.

### Fit the LR model

```import statsmodels.formula.api as sm

These are the beta coefficients calculated:

```In: modelAll.params
Out:
Intercept 2.938889
TV 0.045765
Newspaper -0.001037
dtype: float64```

We interpret these results as follows: for a given amount of TV and newspaper advertising, spending an additional 1000 dollars on radio advertising leads to an increase in sales by approximately 189 units.
In contrast, the coefficient for newspaper represents the average effect (negligible) of increasing newspaper spending by 1000 dollars while holding TV and radio fixed.

## Is at least one of the features useful in predicting Sales?

We use a hypothesis test to answer this question.

The most common hypothesis test involves testing the null hypothesis of:
H0: There is no relationship between the media and sales versus the alternative hypothesis
Ha: There is some relationship between the media and sales.

Mathematically, this corresponds to testing
H0: β1 = β2 = β3 = β4 = 0
versus
Ha: at least one βi is non-zero.

This hypothesis test is performed by computing the F-statistic.

### The F-statistic

We need first of all the Residual Sum of Squares (RSS), i.e. the sum of all squared errors (differences between actual sales and predictions from the regression line). Remember this is the number that the regression is trying to minimise.

```y_pred = modelAll.predict(ad)

import numpy as np

Out: 556.825262902187```

Now we need the Total Sum of Squares (TSS): the total variance in the response Y, and can be thought of as the amount of variability inherent in the response before the regression is performed.
The distance from any point in a collection of data, to the mean of the data, is the deviation.

```y_mean = np.mean(ad.Sales) # mean of sales

In: TSS
Out: 5417.148749999998```

```p=3 # we have three predictors: TV, Radio and Newspaper
n=200 # we have 200 data points (input samples)```
```F = ((TSS-RSS)/p) / (RSS/(n-p-1))

In: F
Out: 570.2707036590941```

When there is no relationship between the response and predictors, one would expect the F-statistic to take on a value close to 1.
On the other hand, if Ha is true, then we expect F to be greater than 1.

In this case, F is far larger than 1: at least one of the three advertising media must be related to sales.

## How strong is the relationship?

Once we have rejected the null hypothesis in favour of the alternative hypothesis, the next question is to quantify the extent to which the model fits the data.

The quality of a linear regression fit is typically assessed using two related quantities: the residual standard error (RSE) and the R2 statistic (the square of the correlation of the response and the variable, when close to 1 means high correlation).

```RSE = np.sqrt((1/(n-2))*RSS)
In: RSE
Out: 1.6769760888385672

Out: 14.022500000000003

In: R2
Out: 0.8972106381789521```

RSE is 1.68 units while the mean value for the response is 14,02 indicating a percentage error of roughly 12%.
The R2 statistic records the percentage of variability in the response that is explained by the predictors. The predictors explain almost 90% of the variance in sales.

## Summary

statsmodels has a handy function that provides the above metrics in one single table:

```In: modelAll.summary()
Out:```
Dep. Variable: R-squared: Sales 0.897 OLS 0.896 Least Squares 570.3 Mon, 08 May 2017 1.58e-96 22:45:01 -386.18 200 780.4 196 793.6 3 nonrobust
coef std err t P>|t| [95.0% Conf. Int.] 2.9389 0.312 9.422 0.000 2.324 3.554 0.0458 0.001 32.809 0.000 0.043 0.049 0.1885 0.009 21.893 0.000 0.172 0.206 -0.0010 0.006 -0.177 0.860 -0.013 0.011
 Omnibus: Durbin-Watson: 60.414 2.084 0 151.241 -1.327 1.44e-33 6.332 454

One thing to note is that R2 (R-squared above) will always increase when more variables are added to the model, even if those variables are only weakly associated with the response.
Therefore an adjusted R2 is provided, which is R2 adjusted by the number of predictors.

Another thing to note is that the summary table provides also a t-statistic and a p-value for each single feature.
These provide information about whether each individual predictor is related to the response (high t-statistic or low p-value).

But be careful looking only at these individual p-values instead of looking at the overall F-statistic. It seems likely that if any one of the p-values for the individual features is very small, then at least one of the predictors is related to the response. However, this logic is flawed, especially when you have many predictors; statistically about 5 % of the p-values will be below 0.05 by chance (this is the effect infamously leveraged by the so-called p-hacking).
The F-statistic does not suffer from this problem because it adjusts for the number of predictors.

## Which media contribute to sales?

To answer this question, we could examine the p-values associated with each predictor’s t-statistic. In the multiple linear regression above, the p-values for TV and radio are low, but the p-value for newspaper is not. This suggests that only TV and radio are related to sales.

But as just seen, if p is large then we are likely to make some false discoveries.

The task of determining which predictors are associated with the response, in order to fit a single model involving only those predictors, is referred to as variable or feature selection.

Ideally, we could perform the feature selection by trying out a lot of different models – each containing a different subset of the features – and then select the best model out of all of the models that we have considered.
For example, the model with the smallest RSS and the biggest R2. Other used metrics are the Mallow’s Cp, Akaike information criterion (AIC), Bayesian information criterion (BIC), and adjusted R2. All of them are visible in the summary model.

```def evaluateModel (model):
print("R2 = ", model.rsquared)```

Unfortunately, there are a total of 2^p models that contain subsets of p variables.
For three predictors, it would still be manageable, only 8 models to fit and evaluate but as p increases, the number of models grows esponentially.

Instead, we can use other approaches.
The three classical ways are the forward selection (start with no features and add one after the other until a threshold is reached); the backward selection (start with all features and remove one by one) and the mixed selection (a combination of the two).

We show here the forward selection.

### Forward selection

We start with a null model (no features), we then fit three (p=3) different simple linear regressions (one for each feature) and add to the null model the variable that results in the lowest RSS.

```modelTV = sm.ols('Sales ~ TV', ad).fit()

In: modelTV.summary().tables
Out:```
coef std err t P>|t| [95.0% Conf. Int.] 7.0326 0.458 15.360 0.000 6.130 7.935 0.0475 0.003 17.668 0.000 0.042 0.053
```In: evaluateModel(modelTV)
Out:
R2 = 0.61187505085```

The model containing only TV as a predictor had an RSS=2103 and an R2 of 0.61

```modelRadio = sm.ols('Sales ~ Radio', ad).fit()
coef std err t P>|t| [95.0% Conf. Int.] 9.3116 0.563 16.542 0.000 8.202 10.422 0.2025 0.020 9.921 0.000 0.162 0.243
```In: evaluateModel(modelRadio)
Out:
R2 = 0.332032455445

modelPaper = sm.ols('Sales ~ Newspaper', ad).fit()
modelPaper.summary().tables
```
coef std err t P>|t| [95.0% Conf. Int.] 12.3514 0.621 19.876 0.000 11.126 13.577 0.0547 0.017 3.300 0.001 0.022 0.087
```In: evaluateModel(modelPaper)
Out:

R2 = 0.0521204454443```

The lowest RSS and the highest R2 are for the TV medium.
Now we have a best model M1 which contains one single feature: TV advertising.
Next step is to add to this M1 model the variable that results in the lowest RSS among the fitted two-variable model (TV + another variable).
This approach is continued until some stopping rule is satisfied.

Now fitting and evaluating the models with combinations of the TV feature plus Radio or Newspaper:

```modelTVRadio = sm.ols('Sales ~ TV + Radio', ad).fit()
Out:
```
coef std err t P>|t| [95.0% Conf. Int.] 2.9211 0.294 9.919 0.000 2.340 3.502 0.0458 0.001 32.909 0.000 0.043 0.048 0.1880 0.008 23.382 0.000 0.172 0.204
```In: evaluateModel(modelTVRadio)
Out:
R2 = 0.897194261083

modelTVPaper = sm.ols('Sales ~ TV + Newspaper', ad).fit()
In: modelTVPaper.summary().tables```
coef std err t P>|t| [95.0% Conf. Int.] 5.7749 0.525 10.993 0.000 4.739 6.811 0.0469 0.003 18.173 0.000 0.042 0.052 0.0442 0.010 4.346 0.000 0.024 0.064
```In: evaluateModel(modelTVPaper)
Out:
R2 = 0.645835493829```

Well, the model with TV AND Radio greatly decreased RSS and increased R2, so that will be our M2 model.

Now, we have only three variables here. We can decide to stop at M2 or use an M3 model with all three variables.
Recall that we already fitted and evaluated a model with all features, just at the beginning:

```In: evaluateModel(modelAll)
Out:
R2 = 0.897210638179```

M3 is *slightly* better than M2 (but remember that the R2 metric always increases when adding new variables) so we call the approach completed and decide that the M2 model with TV and Radio is the best compromise. Adding the newspaper could possibly overfits on new test data.
Next year no budget for newspaper advertising and that amount will be used for TV and Radio instead!

`In: modelTVRadio.summary()`
Dep. Variable: R-squared: Sales 0.897 OLS 0.896 Least Squares 859.6 Tue, 09 May 2017 4.83e-98 17:26:12 -386.20 200 778.4 197 788.3 2 nonrobust
coef std err t P>|t| [95.0% Conf. Int.] 2.9211 0.294 9.919 0.000 2.340 3.502 0.0458 0.001 32.909 0.000 0.043 0.048 0.1880 0.008 23.382 0.000 0.172 0.204
 Omnibus: Durbin-Watson: 60.022 2.081 0 148.679 -1.323 5.19e-33 6.292 425

## Plotting the model

The M2 model has two variables therefore can be plotted as a plane in a 3D chart.

```In: modelTVRadio.params
Out:
Intercept    2.921100
TV           0.045755
dtype: float64```

The M2 model can be described therefore by this equation:
Sales = 0.19 * Radio + 0.05 * TV + 2.9 which I can write as:
0.19*x + 0.05*y – z + 2.9 = 0
Its normal is (0.19, 0.05, -1)
and a point on the plane is (-2.9/0.19,0,0) = (-15.26,0,0)

```normal = np.array([0.19,0.05,-1])
point = np.array([-15.26,0,0])
# a plane is a*x + b*y +c*z + d = 0
# [a,b,c] is the normal.
# Thus, we have to calculate d and we're set
d = -np.sum(point*normal) # dot product
# create x,y
x, y = np.meshgrid(range(50), range(300))
# calculate corresponding z
z = (-normal*x - normal*y - d)*1./normal```

Let’s plot the actual values as red points and the model predictions as a cyan plane: ## Is there synergy among the advertising media?

Adding radio to the model leads to a substantial improvement in R2. This implies that a model that uses TV and radio expenditures to predict sales is substantially better than one that uses only TV advertising.

In our previous analysis of the Advertising data, we concluded that both TV and radio seem to be associated with sales.  The linear models that formed the basis for this conclusion assumed that the effect on sales of increasing one advertising medium is independent of the amount spent on the other media. For example, the linear model M2 states that the average effect on sales of a one-unit increase in TV is always β1, regardless of the amount spent on radio.
However, this assumption may be incorrect. Suppose that spending money on radio advertising actually increases the effectiveness of TV advertising, so that the slope term for TV (β1) should increase as radio increases. In this situation, given a fixed budget of \$100K spending half on radio and half on TV may increase sales more than allocating the entire amount to either TV or to radio. In marketing, this is known as a synergy effect.

The figure above suggests that such an effect may be present in the advertising data. Notice that when levels of either TV or radio are low, then the true sales are lower than predicted by the linear model. But when advertising is split between the two media, then the model tends to underestimate sales.

It is possible to fit a model that integrates a synergy effect via ols():

```modelSynergy = sm.ols('Sales ~ TV + Radio + TV*Radio', ad).fit()

In: modelSynergy.summary().tables
Out:```
coef std err t P>|t| [95.0% Conf. Int.] 6.7502 0.248 27.233 0.000 6.261 7.239 0.0191 0.002 12.699 0.000 0.016 0.022 0.0289 0.009 3.241 0.001 0.011 0.046 0.0011 5.24e-05 20.727 0.000 0.001 0.001

The results strongly suggest that the model that includes the interaction term is superior to the model that contains only main effects. The p-value for the interaction term, TV×radio, is extremely low, indicating that there is strong evidence for Ha : β3 not zero. In other words, it is clear that the true relationship is not additive.

```In: evaluateModel(modelSynergy)
Out:
R2 = 0.967790549848```

The R2 for the model is now 96.8 %, compared to only 89.7% for the M2 model that predicts sales using TV and radio without an interaction term. This means that (96.8 − 89.7)/(100 − 89.7) = 69 % of the variability in sales that remains after fitting the additive model has been explained by the interaction term.
A linear model that uses radio, TV, and an interaction between the two to predict sales takes the form:
sales = β0 + β1 × TV + β2 × radio + β3 × (radio×TV) + ε

```In:  modelSynergy.params
Out:
Intercept 6.750220
TV 0.019101