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: http://www-bcf.usc.edu/~gareth/ISL/data.html
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.
In [1]:
import pandas as pd
In [3]:
ad = pd.read_csv("../datasets/advertising.csv", index_col=0)
In [4]:
ad.info()
In [5]:
ad.describe()
Out[5]:
In [6]:
ad.head()
Out[6]:
In [7]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.scatter(ad.TV, ad.Sales, color='blue', label="TV")
plt.scatter(ad.Radio, ad.Sales, color='green', label='Radio')
plt.scatter(ad.Newspaper, ad.Sales, color='red', label='Newspaper')
plt.legend(loc="lower right")
plt.title("Sales vs. Advertising")
plt.xlabel("Advertising [1000 $]")
plt.ylabel("Sales [Thousands of units]")
plt.grid()
plt.show()
In [8]:
ad.corr()
Out[8]:
In [9]:
plt.imshow(ad.corr(), cmap=plt.cm.Blues, interpolation='nearest')
plt.colorbar()
tick_marks = [i for i in range(len(ad.columns))]
plt.xticks(tick_marks, ad.columns, rotation='vertical')
plt.yticks(tick_marks, ad.columns)
Out[9]:
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.
In [10]:
import statsmodels.formula.api as sm
In [11]:
modelAll = sm.ols('Sales ~ TV + Radio + Newspaper', ad).fit()
These are the beta coefficients calculated:
In [12]:
modelAll.params
Out[12]:
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.
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
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.
In [13]:
y_pred = modelAll.predict(ad)
In [14]:
import numpy as np
RSS = np.sum((y_pred - ad.Sales)**2)
RSS
Out[14]:
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.
In [15]:
y_mean = np.mean(ad.Sales) # mean of sales
In [16]:
TSS = np.sum((ad.Sales - y_mean)**2)
TSS
Out[16]:
The F-statistic is the ratio between (TSS-RSS)/p and RSS/(n-p-1)
In [17]:
p=3 # we have three predictors: TV, Radio and Newspaper
n=200 # we have 200 data points (input samples)
In [18]:
F = ((TSS-RSS)/p) / (RSS/(n-p-1))
F
Out[18]:
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.
Once we have rejected the null hypothesis in favor of the alternative hypothesis, it is natural to want 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).
In [19]:
RSE = np.sqrt((1/(n-2))*RSS);
RSE
Out[19]:
In [20]:
np.mean(ad.Sales)
Out[20]:
In [21]:
R2 = 1 - RSS/TSS;
R2
Out[21]:
RSE is 1.68 units while the mean value for the response is 14.02, indicating a percentage error of roughly 12%.
Second, 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.
statsmodels has a handy function that provides the above metrics in one single table:
In [22]:
modelAll.summary()
Out[22]:
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 shows 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.
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 /feature selection.
Ideally, we could perform the variable selection by trying out a lot of different models, each containing a different subset of the features.
We can 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.
In [23]:
def evaluateModel (model):
print("RSS = ", ((ad.Sales - model.predict())**2).sum())
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 exponentially.
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 try here the forward selection.
We start with a null model (no features), we then fit three (p=3) simple linear regressions and add to the null model the variable that results in the lowest RSS.
In [24]:
modelTV = sm.ols('Sales ~ TV', ad).fit()
modelTV.summary().tables[1]
Out[24]:
In [25]:
evaluateModel(modelTV)
The model containing only TV as a predictor had an RSS=2103 and an R2 of 0.61
In [26]:
modelRadio = sm.ols('Sales ~ Radio', ad).fit()
modelRadio.summary().tables[1]
Out[26]:
In [27]:
evaluateModel(modelRadio)
In [28]:
modelPaper = sm.ols('Sales ~ Newspaper', ad).fit()
modelPaper.summary().tables[1]
Out[28]:
In [29]:
evaluateModel(modelPaper)
The lowest RSS and the highest R2 are for the TV medium.
Now we have a best model M1 which contains TV advertising.
We then add to this M1 model the variable that results
in the lowest RSS for the new two-variable model.
This approach is continued until some stopping rule is satisfied.
In [30]:
modelTVRadio = sm.ols('Sales ~ TV + Radio', ad).fit()
modelTVRadio.summary().tables[1]
Out[30]:
In [31]:
evaluateModel(modelTVRadio)
In [32]:
modelTVPaper = sm.ols('Sales ~ TV + Newspaper', ad).fit()
modelTVPaper.summary().tables[1]
Out[32]:
In [33]:
evaluateModel(modelTVPaper)
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 [34]:
evaluateModel(modelAll)
M3 is slightly better than M2 (but remember that R2 always increases when adding new variables) so we call the approach completed and decide that the M2 model with TV and Radio is the good 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 [35]:
modelTVRadio.summary()
Out[35]:
The M2 model has two variables therefore can be plotted as a plane in a 3D chart.
In [36]:
modelTVRadio.params
Out[36]:
The M2 model can be described by this equation:
Sales = 0.19 Radio + 0.05 TV + 2.9 which I can write as:
0.19x + 0.05y - 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)
In [37]:
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[0]*x - normal[1]*y - d)*1./normal[2]
Let's plot the actual values as red points and the model predictions as a cyan plane:
In [38]:
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
fig.suptitle('Regression: Sales ~ Radio + TV Advertising')
ax = Axes3D(fig)
ax.set_xlabel('Radio')
ax.set_ylabel('TV')
ax.set_zlabel('Sales')
ax.scatter(ad.Radio, ad.TV, ad.Sales, c='red')
ax.plot_surface(x,y,z, color='cyan', alpha=0.3)
Out[38]:
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 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 simple model may be incorrect. Suppose that spending money on radio advertising actually increases the effectiveness of TV advertising, so that the slope term for TV 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.
In [39]:
modelSynergy = sm.ols('Sales ~ TV + Radio + TV*Radio', ad).fit()
modelSynergy.summary().tables[1]
Out[39]:
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 [40]:
evaluateModel(modelSynergy)
The R2 for this model is 96.8 %, compared to only 89.7% for the model M2 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 [41]:
modelSynergy.params
Out[41]:
We can interpret β3 as the increase in the effectiveness of TV advertising for a one unit increase in radio advertising (or vice-versa).
In [ ]: