Minimal Linear Regression Graphics
Minimal Visualization of Regression Results¶
Introduction¶
I was working my way through An Introduction to Statistical Learning (https://www.statlearning.com/), and (as usual) the examples in the book were in R (to be expected from professional statisticians). My response should be "I'll learn R to replicate these!", but it is usually "How would I do this in Python?".
This post is how I replicated some graphics associated with linear regression in Python. The situation is that we have data on sales of some product in various areas, with details of money spent in advertising in three channels (TV, newspapers, and radio). We try to model the effectiveness of the three channels, based upon the sales data, and then visualize the models.
Implementation¶
%matplotlib inline
%load_ext watermark
%load_ext lab_black
# all imports should go here
import pandas as pd
import numpy as np
import statsmodels.api as sm
from statsmodels.formula.api import ols
import statsmodels.graphics.api as smg
# housekeeping imports
import sys
import os
import subprocess
import datetime
import platform
import datetime
# graphic imports
import seaborn as sns
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# path management
from pathlib import Path
data_dir_path = Path('d:/IntroToStatLearning/')
data_path = data_dir_path / 'Advertising.csv'
ads = pd.read_csv(data_path)
ads.head()
Analysis¶
In the analysis phase. we first try to fit a linear relationship between the individual spend in advertsing channels, and the sales. We start with the TV spend.
TV Alone¶
We peform a Ordinary Least Squares fit of a linear relationship between TV spend, and sales, and print s asummary of the results.
res1 = ols('sales ~ TV ', data=ads).fit()
res1.summary()
We note that there is (almost certainly) a relationship between TV spend and sales, but with R^2 = 0.612, there is a lot of variation in the sales data left unexplained by this simple model.
Visualization¶
Now we visualize the regression model we have just created.
We start by ploting the sales data against TV spend (in practise, we would probably do this as part of initial data exploration). We use matplotlib quick and minimal plot facilities, but will add more detail later.
plt.plot(
ads['TV'], ads['sales'], 'r+',
)
We can see that the variation of the sales data is not uniform, but there is a definite upwards trend.
Now we plot the OLS line overv the range of the TV spend data. We use the Regression Results params member to get the coefficients of the model.
x_range = np.asarray([min(ads['TV']), max(ads['TV'])])
plt.plot(
x_range,
res1.params['Intercept'] + res1.params['TV'] * x_range,
'g-',
)
Now we put these two plots together. For each data point, we draw a small red dot ('ro', ms=2), and then draw a thin (lw=1) faint (alpha=0.2) black line to the regression line. Then we draw the regression line in green, and label the X and Y axis.
fig, ax = plt.subplots(figsize=(12, 8),)
for x, y_seen, y_fit in zip(
ads['TV'], ads['sales'], res1.fittedvalues
):
ax.plot(x, y_seen, 'ro', ms=2)
ax.plot([x, x], [y_seen, y_fit], 'k-', lw=1, alpha=0.2)
# end for
x_range = np.asarray([min(ads['TV']), max(ads['TV'])])
ax.plot(
x_range,
res1.params['Intercept'] + res1.params['TV'] * x_range,
'g-',
)
ax.set_xlabel('TV')
ax.set_ylabel('Sales')
There are clearly issues with the model as depicted above. The errors (actual-predicted) are not uniform, and for small values of TV spend, the error is consistently negative.
For interest, we show the mean and observation 95% Confidence Intervals (CIs) below. We get a linear spread of values of TV spend, and use the get_prediction method to get a DataFrame with CI values for these input values
We then add these Confidence Intervals to the graphic we had before.
x_ci = np.linspace(0, 300, 20)
gp = res1.get_prediction({'TV': x_ci},)
pred_df = gp.summary_frame()
fig, ax = plt.subplots(figsize=(12, 8),)
for x, y_seen, y_fit in zip(
ads['TV'], ads['sales'], res1.fittedvalues
):
ax.plot(x, y_seen, 'ro', ms=2)
ax.plot([x, x], [y_seen, y_fit], 'k-', lw=1, alpha=0.2)
# end for
ax.plot(
x, y_seen, 'ro', ms=2, label='Actual',
)
x_range = np.asarray([min(ads['TV']), max(ads['TV'])])
ax.plot(
x_range,
res1.params['Intercept'] + res1.params['TV'] * x_range,
'g-',
label='Fitted Line',
)
ax.set_xlabel('TV')
ax.set_ylabel('Sales')
gp = res1.get_prediction({'TV': x_ci})
pred_df = gp.summary_frame()
ax.plot(x_ci, pred_df['mean_ci_upper'], 'b-')
ax.plot(
x_ci, pred_df['mean_ci_lower'], 'b-', label='Mean CI',
)
ax.plot(x_ci, pred_df['obs_ci_upper'], 'b:')
ax.plot(
x_ci, pred_df['obs_ci_lower'], 'b:', label='Obs. CI',
)
ax.legend()
The summary_frame of the predictions looks like:
pred_df.head()
Further simple models¶
Next we construct individual linear models for the relationship between sales and spend on radio and newspapers. In summary, these explain little of the variance in the sales dataset (R^2 values).
res1 = ols('sales ~ radio ', data=ads).fit()
res1.summary()
res1 = ols('sales ~ newspaper ', data=ads).fit()
res1.summary()
Linear model with all variables¶
If we fit a linear model with all variables included, we get a much better explanation of the variation in the salses dataset (R^2 = 0.897). Unexpectedly, we find a weakly negative relationship between newspaper spend and sales! In fact, the 95% CI for the newspaper coefficient includes zero.
res4 = ols('sales ~ TV + newspaper + radio', data=ads).fit()
res4.summary()
We can use numpy to get the correlation between the various variables, and pandas to turn this into a DataFrame for a more polished presentation. We find that the newspaper spend is most strongly correlated with the radio spend.
cm = np.corrcoef(
ads[['TV', 'radio', 'newspaper', 'sales']],
rowvar=False,
)
cm
We show the matrix in a more civilized format
cm_df = pd.DataFrame(
data=cm,
columns=['TV', 'radio', 'newspaper', 'sales'],
index=['TV', 'radio', 'newspaper', 'sales'],
)
cm_df
Final linear model¶
In our final purely linear model, we drop the newspaper spend (and find no reduction in the explained variation in the sales dataset).
res5 = ols('sales ~ TV + radio', data=ads).fit()
res5.summary()
The parameters of this model can be accessed via the params attribute of the regression results.
res5.params
n_points will be the number of points we use in plotting results
n_points = 20
We now plot the actuals vs predicted, using the 3D features of matplotlib
fig = plt.figure(figsize=(12, 8),)
ax = fig.add_subplot(111, projection='3d')
# draw the raw data points
ax.scatter(
ads['TV'],
ads['radio'],
ads['sales'],
c='red',
alpha=0.5,
)
# label the X, Y, Z axis
ax.set_xlabel('TV')
ax.set_ylabel('radio')
ax.set_zlabel('sales')
# plot faint vertical lines from the data points to X,Y plane
for x, y, z in zip(ads['TV'], ads['radio'], ads['sales']):
ax.plot(
[x, x], [y, y], [0, z], 'b-', alpha=0.1,
)
# end for
# plot thicker red lines from each data point to predicted value for that data point
for x, y, z1, z2 in zip(
ads['TV'], ads['radio'], ads['sales'], res5.fittedvalues
):
ax.plot(
[x, x], [y, y], [z1, z2], 'r-', alpha=0.8,
)
# end for
# get the linear span of the X and Y axis values
tv_range = np.linspace(
min(ads['TV']), max(ads['TV']), n_points
)
radio_range = np.linspace(
min(ads['radio']), max(ads['radio']), n_points
)
# get datapoints for lines in the vertical planes for X = 0, X= Max(X) and Y=0, Y=Max(Y)
line_tv = (
tv_range * res5.params['TV'] + res5.params['Intercept']
)
line_radio = (
radio_range * res5.params['radio']
+ res5.params['Intercept']
)
line2_tv = (
tv_range * res5.params['TV']
+ res5.params['Intercept']
+ np.ones(n_points)
* max(radio_range)
* res5.params['radio']
)
line2_radio = (
radio_range * res5.params['radio']
+ res5.params['Intercept']
+ np.ones(n_points) * max(tv_range) * res5.params['TV']
)
# plot lines in the vertical planes for X = 0, X= Max(X) and Y=0, Y=Max(Y)
ax.plot(tv_range, np.zeros(n_points), line_tv, 'g-')
ax.plot(np.zeros(n_points), radio_range, line_radio, 'g-')
ax.plot(
tv_range,
np.ones(n_points) * max(radio_range),
line2_tv,
'g-',
)
ax.plot(
np.ones(20) * max(tv_range),
radio_range,
line2_radio,
'g-',
)
# plot the predicted values as a mesh grid, and as a colored surface
X, Y = np.meshgrid(tv_range, radio_range)
Z = (
X * res5.params['TV']
+ Y * res5.params['radio']
+ res5.params['Intercept']
)
surf = ax.plot_wireframe(X, Y, Z, color='green', alpha=0.4)
surfs = ax.plot_surface(X, Y, Z, color='green', alpha=0.1)
# set the viewing angle
ax.view_init(elev=20, azim=60)
We can see that the prediction errors are not distributed uniformly across the range of the datasets (e.g. the prediction errors are uniformly positive at the left and right corners of the prediction plane, as shown above (predicted > actual).
We can also plot the raw data values in the X=0, and Y=0 planes, as below.
fig = plt.figure(figsize=(12, 8),)
ax = fig.add_subplot(111, projection='3d')
# do 3d scatter plot
ax.scatter(
ads['TV'],
ads['radio'],
ads['sales'],
c='red',
alpha=0.5,
)
# set X,Y,Z axis labels
ax.set_xlabel('TV')
ax.set_ylabel('radio')
ax.set_zlabel('sales')
# draw faint line from scatter point to TV/radio plane (sales==0)
for x, y, z in zip(ads['TV'], ads['radio'], ads['sales']):
ax.plot(
[x, x], [y, y], [0, z], 'b-', alpha=0.1,
)
# end for
# show raw datapoints in radio=0 plane
ax.scatter(
ads['TV'],
np.zeros(len(ads)),
ads['sales'],
c='blue',
alpha=0.3,
)
# show raw datapoints in TV=0 plane
ax.scatter(
np.zeros(len(ads)),
ads['radio'],
ads['sales'],
c='blue',
alpha=0.3,
)
# draw red line from scatter point to fitted point
for x, y, z1, z2 in zip(
ads['TV'], ads['radio'], ads['sales'], res5.fittedvalues
):
ax.plot(
[x, x], [y, y], [z1, z2], 'r-', alpha=0.8,
)
# end for
# get range of X, Y (or TV, radio) values for plotting
tv_range = np.linspace(
min(ads['TV']), max(ads['TV']), n_points
)
radio_range = np.linspace(
min(ads['radio']), max(ads['radio']), n_points
)
# get fitted line values on radio==0 plane, and tv==0 plane
line_tv = (
tv_range * res5.params['TV'] + res5.params['Intercept']
)
line_radio = (
radio_range * res5.params['radio']
+ res5.params['Intercept']
)
# get fitted line values on radio=max(radio) plane
line2_tv = (
tv_range * res5.params['TV']
+ res5.params['Intercept']
+ np.ones(n_points)
* max(radio_range)
* res5.params['radio']
)
# get fitted line values on tv=max(tv) plane
line2_radio = (
radio_range * res5.params['radio']
+ res5.params['Intercept']
+ np.ones(n_points) * max(tv_range) * res5.params['TV']
)
# draw fitted lines in the tv==0 , and the radio==0 planes
ax.plot(tv_range, np.zeros(n_points), line_tv, 'g-')
ax.plot(np.zeros(n_points), radio_range, line_radio, 'g-')
# draw fitted lines in the tv== max(tv), and the radio==max(radio) planes
ax.plot(
tv_range,
np.ones(n_points) * max(radio_range),
line2_tv,
'g-',
)
ax.plot(
np.ones(20) * max(tv_range),
radio_range,
line2_radio,
'g-',
)
X, Y = np.meshgrid(tv_range, radio_range)
Z = (
X * res5.params['TV']
+ Y * res5.params['radio']
+ res5.params['Intercept']
)
# draw wireframe and surface
surf = ax.plot_wireframe(X, Y, Z, color='green', alpha=0.4)
surfs = ax.plot_surface(X, Y, Z, color='green', alpha=0.1)
ax.view_init(elev=50, azim=50)
Interaction Effects¶
In our final model, we consider an interaction effect. The idea is that maybe spending on TV increases the effectiveness of radio, and vice versa. We find an interaction effect that is statistically greater than zero. R^2 has gone from 0.897 to 0.968, so almost all of the variation in the ads dataset is explained by this model.
res5 = ols('sales ~ TV * radio', data=ads).fit()
res5.summary()
Final model visualized¶
When we visualize the model, we can see that the distribution of the residuals (actual-predicted) is now much more uniform.
fig = plt.figure(figsize=(12, 8),)
ax = fig.add_subplot(111, projection='3d')
# do 3d scatter plot
ax.scatter(
ads['TV'],
ads['radio'],
ads['sales'],
c='red',
alpha=0.5,
)
# set X,Y,Z axis labels
ax.set_xlabel('TV')
ax.set_ylabel('radio')
ax.set_zlabel('sales')
# draw faint line from scatter point to TV/radio plane (sales==0)
for x, y, z in zip(ads['TV'], ads['radio'], ads['sales']):
ax.plot(
[x, x], [y, y], [0, z], 'b-', alpha=0.1,
)
# end for
# draw line from scatter point to fitted point
for x, y, z1, z2 in zip(
ads['TV'], ads['radio'], ads['sales'], res5.fittedvalues
):
ax.plot(
[x, x], [y, y], [z1, z2], 'r-', alpha=0.8,
)
# end for
X, Y = np.meshgrid(tv_range, radio_range)
Z = (
X * res5.params['TV']
+ Y * res5.params['radio']
+ X * Y * res5.params['TV:radio']
+ res5.params['Intercept']
)
# draw wireframe and surface
surf = ax.plot_wireframe(X, Y, Z, color='green', alpha=0.4)
surfs = ax.plot_surface(X, Y, Z, color='green', alpha=0.1)
Summary¶
The comparion with the graphics in ITSL is interesting (I assume the graphics in ITSL are done with R). By default, matplotlib 3D places ticks and tick labels, and shows the 3D aspect by a grid in each of the figure 'back planes'. The ITSL graphics are very much paired back, with no ticks or tick labels, and no grids. I can achieve the same result, as below.
fig = plt.figure(figsize=(12, 8),)
ax = fig.add_subplot(111, projection='3d')
# clear ticks and labels
ax.set_zticks([])
ax.set_yticks([])
ax.set_xticks([])
# clear background panes
ax.xaxis.pane.fill = False
ax.yaxis.pane.fill = False
ax.zaxis.pane.fill = False
# turn off wireframe at back of 3D plot
ax.zaxis.pane.set_edgecolor('white')
ax.yaxis.pane.set_edgecolor('white')
ax.xaxis.pane.set_edgecolor('white')
# do 3d scatter plot
ax.scatter(
ads['TV'],
ads['radio'],
ads['sales'],
c='red',
alpha=0.5,
)
# set X,Y,Z axis labels
ax.set_xlabel('TV')
ax.set_ylabel('radio')
ax.set_zlabel('sales')
# draw faint line from scatter point to TV/radio plane (sales==0)
for x, y, z in zip(ads['TV'], ads['radio'], ads['sales']):
ax.plot(
[x, x], [y, y], [0, z], 'b-', alpha=0.1,
)
# end for
# draw line from scatter point to fitted point
for x, y, z1, z2 in zip(
ads['TV'], ads['radio'], ads['sales'], res5.fittedvalues
):
ax.plot(
[x, x], [y, y], [z1, z2], 'r-', alpha=0.8,
)
# end for
X, Y = np.meshgrid(tv_range, radio_range)
Z = (
X * res5.params['TV']
+ Y * res5.params['radio']
+ X * Y * res5.params['TV:radio']
+ res5.params['Intercept']
)
# draw wireframe and surface
surf = ax.plot_wireframe(X, Y, Z, color='green', alpha=0.4)
surfs = ax.plot_surface(X, Y, Z, color='green', alpha=0.1)
Conclusion¶
Very similar minimalist graphics can be achieved in Python, to match those shown in the ITSL book.
Reproducibility¶
Notebook version status¶
theNotebook = 'ISLR-LinReg'
# show info to support reproducibility
def python_env_name():
envs = subprocess.check_output(
'conda env list'
).splitlines()
# get unicode version of binary subprocess output
envu = [x.decode('ascii') for x in envs]
active_env = list(
filter(lambda s: '*' in str(s), envu)
)[0]
env_name = str(active_env).split()[0]
return env_name
# end python_env_name
print('python version : ' + sys.version)
print('python environment :', python_env_name())
print('pandas version : ' + pd.__version__)
print('current wkg dir: ' + os.getcwd())
print('Notebook name: ' + theNotebook)
print(
'Notebook run at: '
+ str(datetime.datetime.now())
+ ' local time'
)
print(
'Notebook run at: '
+ str(datetime.datetime.utcnow())
+ ' UTC'
)
print('Notebook run on: ' + platform.platform())
%watermark
%watermark -h -iv
import matplotlib
matplotlib.__version__