Stanford STATS191 in Python, Lecture 6 : Diagnostics for multiple linear regression
Stanford Stats 191¶
Introduction¶
This is a re-creation of the Stanford Stats 191 course (see https://web.stanford.edu/class/stats191/), using Python eco-system tools, instead of R. This is lecture "Diagnostics for multiple regression"
Initial Notebook Setup¶
watermark documents the current Python and package environment, black is my preferred Python formatter
%load_ext watermark
%load_ext lab_black
%matplotlib inline
All imports go here
import pandas as pd
import numpy as np
import seaborn as sn
import math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy import stats
from statsmodels.formula.api import ols
from statsmodels.formula.api import rlm
import statsmodels.api as sm
from statsmodels.stats.stattools import jarque_bera
import os
data = pd.read_csv('../data/hills.csv', sep='\t')
data.head()
data.sort_values(by='Distance').tail()
Visualize Data¶
We use Seaborn to visually explore the relationships between the variables
sn.pairplot(data)
Same graph as before, but we also add a regression line to each scatter plot
sn.pairplot(data, kind='reg')
We choose to highlight a particular race that looks unusual. We do this by linking the color of the dot (hue) to the race name
data['RaceName'] = [
'LairigGhru' if r == 'LairigGhru' else '-'
for r in data['Race']
]
sn.pairplot(
data, kind='scatter', diag_kind='hist', hue='RaceName'
)
We can also use pandas to do our scatter plots
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
__ = pd.plotting.scatter_matrix(
data, ax=ax, alpha=0.95, grid=False
)
Perform OLS¶
We build a model where time to complete the race depends linearly upon the distance, and height to climb (and of course, a Gaussian error term)
res = ols('Time ~ Distance + Climb', data=data).fit()
res.summary()
_ = plt.plot(res.resid, 'o', label='Data')
_ = plt.legend(loc='best')
_ = plt.grid()
_ = plt.xlabel('Observation Number')
_ = plt.ylabel('Residual')
_ = plt.legend()
We can also plot our two explanatory variables against the predicted and actual values of time
_ = plt.plot(
data['Distance'],
res.predict(),
'o',
label='Predicted Time',
)
_ = plt.plot(
data['Distance'],
data['Time'],
'r+',
label='Actual Time',
)
_ = plt.legend(loc='best')
_ = plt.grid()
_ = plt.xlabel('Distance')
_ = plt.ylabel('Time')
_ = plt.legend()
plt.plot(
data['Climb'],
res.predict(),
'o',
label='Predicted Time',
)
plt.plot(
data['Climb'], data['Time'], 'r+', label='Actual Time'
)
_ = plt.legend(loc='best')
_ = plt.grid()
_ = plt.xlabel('Climb')
_ = plt.ylabel('Time')
_ = plt.legend()
We can graph a more detailed view of the residuals, and of residuals plotted against prediction
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
ax.plot(
range(len(res.resid)),
res.resid,
'o',
label='OLS Residuals',
)
ax.grid()
ax.set_xlabel('Residual No.')
ax.set_ylabel('Residual')
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.legend(loc='best')
plt.plot(res.predict(), res.resid, 'o', label='Residual')
_ = plt.legend(loc='best')
_ = plt.grid()
_ = plt.xlabel('Predicted Time')
_ = plt.ylabel('Residual')
_ = plt.axhline(0, c='k')
_ = plt.legend()
Display Normality Test Plots¶
One way to check for the normality of the residuals is to plot the actual quantiles against theoretical quantiles (assuming a normal distribution). We see that because of two outliers, we do not get the expected straight line
There are two ways of getting QQ Plots, via statsmodels (qqplot() ), and via scipy.stats (probplot() )
fig = sm.qqplot(res.resid, line='s')
__ = stats.probplot(res.resid, dist='norm', plot=plt)
A Normal distribution should be a straight line, as below
_ = stats.probplot(
np.random.normal(0, 1, 50), dist='norm', plot=plt
)
(osm, osr), (slope, intercept, r) = stats.probplot(
res.resid, dist='norm', plot=None
)
plt.plot(osm, osr, 'r+')
plt.plot(osm, osm * slope + intercept, 'g-')
Plot Studentized Residuals¶
Plot Studentized Residuals against Residuals.
Externally Studentized means taking all other residuals (ie excluding the residual under consideration) to get a 'standard deviation', to scale the residual to get a variable that is distributed as per t distribution.
Internally Studentized means taking all residuals (ie including the residual under consideration) to get a 'standard deviation', to scale the residual to get a variable that is distributed approximately as per t distribution.
The RegressionResults object has a method outlier_test, that returns the two types of "Studentized" residuals
res_df = res.outlier_test()
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
ax.plot(
res.resid,
res_df['student_resid'],
'o',
label='Studentized Resid.',
)
ax.set_xlabel('Residual')
ax.set_ylabel(
'Studentized Residual (= External Studentized Residual = rstudent in R)'
)
ax.grid()
ax.legend(loc='best')
As we see above, we have a smallnumber of residuals that require closer attention
The outlier_test() method returns a pandas DataFrame as below. The column student_resid is the Externally Studentized Residuals as discussed above. The columns unadj_p and bonf(p) columns are discussed below
res_df.head()
We can also get the standardized residuals (both types) from the OLSInfluence object, using a method get_influence() in the RegressionResults object.
influence = res.get_influence()
type(influence)
standardized_residuals_int = (
influence.resid_studentized_internal
)
standardized_residuals_ext = (
influence.resid_studentized_external
)
Outlier_test results (Bonferroni Adjustment)¶
The approach we could take is to look at the residuals for each race (actual time - predicted time), and reject those observations where the residual is unexpectedly large. That is, larger than a pre-chosen threshold (e.g. a residual of this size would only happen at a frequency of 5%, or 1%, or at whatever our chosen cutoff frequency was). The problem with this approach is that the error term in our model is random, and if we have a very large number of valid data points, then we expect 5% of them to be exceed our threshold, if we have set 5% as our threshold.
Carlo Emilio Bonferroni (an Italian mathematician) developed a correction that reduces the number of valid data points being rejected.
The DataFrame returned by outlier_test , holds (for each observation) the probability that we would see a residual of this size by random chance (assuming all residuals are T-distributed). The bonfp() column I understand to be probability that the observation is valid (after the Bonferroni adjustment).
res_df.head()
We can use develop a boolean mask that extracts those observations that are very surprising (i.e. almost certainly not valid)
outlier_m = [r < 0.05 for r in res_df['bonf(p)']]
print(res_df[outlier_m])
print(data[outlier_m])
We can plot the unadjusted and adjusted probabilities, and a 5% cutoff line. We see the one point below the cutoff.
idf = range(1, len(data) + 1)
plt.plot(idf, res_df['bonf(p)'], 'o', label='Adjusted Prob')
plt.plot(idf, res_df['unadj_p'], 'r+', label='Raw Prob')
plt.axhline(0.05, c='g', label='Threshhold')
plt.legend(loc='best')
plt.xlabel('Observation No')
plt.ylabel('Probability')
It turns out that the datset has a known error (http://www.statsci.org/data/general/hills.html)
The data contains a known error - Atkinson (1986) reports that the record for Knock Hill (observation 17) should actually be 18 minutes rather than 78 minutes.
Extended Analysis of Residuals¶
We see that outlier_test() returns the same array as the influence.resid_studentized_external
we will compare the externally and internally standardized residuals for each point; again we see a few suspicious points
res_df.head()
print(
'outlier_test student_resid ',
standardized_residuals_ext[0:5],
)
standardized_residuals_int = (
influence.resid_studentized_internal
)
standardized_residuals_ext = (
influence.resid_studentized_external
)
print(
'External Stud Resid ',
standardized_residuals_ext[0:5],
)
print(
'Internal Stud Resid ',
standardized_residuals_int[0:5],
)
Display the plots of Externally Studentized Residuals versus the Residuals, and Internally Studentized Residuals. Valid points should like close to a straight line. Again we see a few points that are 'suspicious'
fig, ax = plt.subplots(1, 2, figsize=(10, 8))
ax[0].plot(
res.resid,
standardized_residuals_ext,
'ro',
label='Ext Stud Resid',
)
ax[0].set_xlabel('Residual')
ax[0].set_ylabel(
'External Studentized Residual (= rstudent in R)'
)
ax[0].grid()
ax[0].legend(loc='best')
ax[1].plot(
standardized_residuals_int,
standardized_residuals_ext,
'ro',
label='Int Stud Resid',
)
ax[1].set_ylabel(
'External Studentized Residual (= rstudent in R)'
)
ax[1].set_xlabel(
'Internal Studentized Residual (= rstandard in R)'
)
ax[1].grid()
ax[1].legend(loc='best')
fig.suptitle(
'Standard Residuals\n(as per\n'
+ 'https://web.stanford.edu/class/stats191/notebooks/Diagnostics_for_multiple_regression.html'
)
Constant Variance Test¶
We look at what is some sense the 'variance' of each externally standardize residual, and again can identify a few points that seem different from the others. Otherwise the variance appears to be constant
zz = np.sqrt(np.abs(standardized_residuals_ext))
fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(res.predict(), zz, 'o', label='Data')
ax.set_xlabel('Fitted Value')
ax.set_ylabel('sqrt(abs(Externally Studentized Residual))')
ax.grid()
ax.axhline(1, c='g', label='Outlier Cutoff')
ax.legend(loc='best')
We can plot the residuals against the predicted values. Again variance appear constant across the range of predicted values
plt.plot(res.predict(), res.resid, 'o')
_ = plt.xlabel('Predicted Time')
_ = plt.ylabel('Residual')
Measure of Influence¶
DFFITS¶
DFFITS (standing for Difference in Fits) is a measure of how influential an observation is, gauged by seeing how the prediction for an observation changes if the observation is excluded from the regression.
The OLSInfluence object can be used to get this value for each observation
# df[0] is an array of the DFFITS value for each observation
# df[1] is the suggested threshhold for suspicious DFFITS values (= 2 * sqrt(k / n),
# where k is the size of our model [=3, constant, distance, climb])
# where n is the dataset size
df = influence.dffits
df_val = df[0]
df_thresh = df[1]
ldf = len(df_val)
print(f'DFFITS threshhold value = {df_thresh: 15.5f}')
Plot our DFFITS, and the threshhold
idf = range(ldf)
plt.plot(idf, df_val, 'o')
plt.axhline(df_thresh, c='g')
Show the races for which the DFFITS indicates to be investigated
outlier = [r > df_thresh for r in df_val]
data[outlier]
Checking the DFFITS threshhold value returned by the OLSInfluence object
n_feat = 2
thresh = 2 * np.sqrt((n_feat + 1) / len(data['Time']))
print('Outlier Threshhold ', thresh)
Cooks Distance¶
Cooks distance is a measure of how much the whole regression changes, when we delete an observation. We can get this from the OLSInfluence object. Two arrays are returned, the first matches the R convention
cooks = influence.cooks_distance
Plot the Cooks Distance for each observation
plt.plot(idf, cooks[1], 'o', label='cooks[1]')
plt.plot(
idf, cooks[0], 'r+', label='cooks[0] = R'
) # R definition
plt.xlabel('Observation No')
plt.ylabel('Cooks Distance')
plt.legend(loc='best')
Again we get mask in the races with Cooks Distance over a threshhold, and we get our usual suspects
outlier = [r > 0.1 for r in cooks[0]]
data[outlier]
DFBETA¶
DFBETA measures the difference in the estimated coefficients when a single observation is removed from the OLS fit. We can get this from the OLSInfluence object.
Because we are estimating three parameters (constant, climb, distance), we get three arrays back
dfbetas = influence.dfbetas
type(dfbetas)
Plot the influence observations have over the coefficient of Climb
dfb0 = [r[0] for r in dfbetas]
dfb1 = [r[1] for r in dfbetas]
dfb2 = [r[2] for r in dfbetas]
plt.plot(idf, dfb2, 'o')
# Climb
Plot the influence each observation has on the Distance coefficient
plt.plot(idf, dfb1, 'o')
# Distance
Combine the two plots, and draw threshhold at +-2* sqrt(n)
plt.plot(idf, dfb2, 'o') # Climb
plt.plot(idf, dfb1, 'o') # Distance
plt.axhline(2 / np.sqrt(len(data['Time'])), c='g')
plt.axhline(-2 / np.sqrt(len(data['Time'])), c='g')
Investigation of outliers¶
Given that we have a set of unusual races that a number of metrics have fingered, we can use seaborn the show these in the scatter plots. We explicitly name each suspect race in a new DataFrame column, as assign a placeholder name to all the rest.
data['RaceName'] = [
r
if r in ['BensofJura', 'LairigGhru', 'KnockHill']
else '-'
for r in data['Race']
]
sn.pairplot(
data, kind='scatter', diag_kind='hist', hue='RaceName'
)
One observation we can make is that two of our outliers are indeed outliers in input space, being the longest and highest race. Having outliers like this might make us question our linear model. We might ask if runners operate in two biological regimes, sprints and marathons? Should a model that fits short distances be expected to work at long distances?
Outliers like KnockHill we are more likely attribute to human error (as indeed was the case here, reportedly)
Leverage¶
We can get a view as to how far an observation is from the norm in input space by using the leverage concept, defined as below. We can get this from the OLSInfluence object.
hat = influence.hat_matrix_diag
Plot the leverage values
plt.plot(idf, hat, 'o')
plt.xlabel('Observation No.')
plt.ylabel('Leverage')
Print out the outlier races. We get our races at extremes of distance and height
outlier = [r > 0.3 for r in hat]
data[outlier]
We can display those observation with large residuals, and large leverage
plt.plot(hat, standardized_residuals_int, 'o')
plt.xlabel('Leverage')
plt.ylabel('Internally Standardised Residuals')
A slightly better version of the graph above, done with the OO interface to Matplotlib
fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(hat, standardized_residuals_int, 'o', label='Data')
ax.set_xlabel('Influence')
ax.set_ylabel('Standard Residual (Internal)')
ax.grid()
ax.legend(loc='best')
statsmodels has a purpose-build graphic for visualizing observation influence. Once again, we get our usual suspects.
fig, ax = plt.subplots(figsize=(10, 8))
_ = sm.graphics.influence_plot(res, ax=ax)
data.iloc[[17, 6, 10]]
statsmodels has another closely related graphic to plot leverage statistics vs. normalized residuals squared
fig, ax = plt.subplots(figsize=(10, 8))
_ = sm.graphics.plot_leverage_resid2(res, ax=ax)
Excluding Large Residuals¶
The example below shows the impact of Bonferroni correction. We plot standardised residuals against fitted value, and set thresholds based upon the T test (at the 95% level)
fig, ax = plt.subplots(figsize=(10, 8))
n_features = 2
thresh1, thresh2 = stats.t.interval(
0.95, len(data['Time']) - n_features - 2
)
ax.plot(idf, standardized_residuals_ext, 'o', label='Data')
ax.set_xlabel('Fitted Value')
ax.set_ylabel('Externally Studentized Residual)')
ax.grid()
ax.axhline(thresh1, c='g', label='Outlier Cutoff')
ax.axhline(thresh2, c='g', label='Outlier Cutoff')
ax.legend(loc='best')
print(thresh1, thresh2)
Because we are looking at many residuals, even under the Null Hypothesis (all observations valid), some will be rejected. An example below, where all points are genuinely from Normal(0,1), still has some rejected.
sample = np.random.normal(0, 1, 100)
n_features = 1
thresh1, thresh2 = stats.t.interval(
0.95, len(sample) - n_features - 2
)
plt.plot(range(100), sample, 'o')
plt.axhline(thresh1, c='g', label='Outlier Cutoff')
plt.axhline(thresh2, c='g', label='Outlier Cutoff')
print(thresh1, thresh2)
Apply Bonferroni Correction (widen acceptance region)¶
We widen the acceptance region by the Bonferroni correction
fig, ax = plt.subplots(figsize=(10, 8))
n_points = len(data['Time'])
n_features = 2
alpha = 0.05
thresh1, thresh2 = stats.t.interval(
1 - alpha / n_points, len(data['Time']) - n_features - 2
)
ax.plot(idf, standardized_residuals_ext, 'o', label='Data')
ax.set_xlabel('Fitted Value')
ax.set_ylabel('Externally Studentized Residual)')
ax.grid()
ax.axhline(thresh1, c='g', label='Outlier Cutoff')
ax.axhline(thresh2, c='g', label='Outlier Cutoff')
ax.legend(loc='best')
print(thresh1, thresh2)
We see that under our new threshholds, onlyone observation is excluded, and it is the one point that we know (from other sources) is not valid
outlier = [
np.abs(r) > thresh2 for r in standardized_residuals_ext
]
data[outlier]
The data contains a known error - Atkinson (1986) reports that the record for Knock Hill (observation 17) should actually be 18 minutes rather than 78 minutes.
Additional Plots¶
There are two additional plots that can be used to assess the linear model. STATS191 calls them "added-variable and component plus residual plots"; statsmodel calls them "partial regression plot for a set of regressors" and "component and component-plus-residual (CCPR) plots" respectively
Both graphics seem to support our linear model
CCPR¶
fig = plt.figure(figsize=(8, 8))
_ = sm.graphics.plot_ccpr_grid(res, fig=fig)
plt.show()
Partial Regression Plot¶
fig = plt.figure(figsize=(8, 8))
_ = sm.graphics.plot_partregress_grid(res, fig=fig)
Reproducibility¶
%watermark -h -iv
%watermark