Stanford STATS191 in Python, Lecture 4 : Simple Diagnostics for Linear Models
Stanford Stats 191¶
Introduction¶
This is a re-creation of the Stanford Stats 191 course, using Python eco-system tools, instead of R. This is lecture "Diagnostics for simple linear regression" ( see https://web.stanford.edu/class/stats191/notebooks/Simple_diagnostics.html)
We look at how we can assert that a linear model is appropriate (or not, as the case may be)
Initial Notebook Setup¶
watermark documents the Python and package environment, black is my chosen 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.sandbox.regression.predstd import (
wls_prediction_std,
)
from statsmodels.stats.stattools import jarque_bera
from statsmodels.stats.diagnostic import het_white
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.diagnostic import het_goldfeldquandt
import os
Create Data¶
To illustrate the basic concepts of linear regression, we create a dataset, where y depends upon x, and has added Gaussion noise.
x = np.linspace(0, 20, 21)
y = 0.5 * x + 1 + np.random.normal(0, 1, 21)
df_dict = {'x': x, 'y': y}
data = pd.DataFrame(df_dict)
res = ols('y ~ x', data=data).fit()
Plot the data. We choose to use high-level matplotlib methods, at the expense of some plotting elegance
_ = plt.plot(x, y, 'o')
Plot the data and line of best fit
_ = plt.plot(x, y, 'o')
_ = plt.plot(x, res.predict(), 'r-')
_ = plt.axhline(y.mean(), color='k')
The graphs below illustrate the definitions of various sums of squares; SST is Sun Of Squares Total (deviation from the mean of the actual data), SSE is the Sum of Squares of Errors (deviation from the actuals of the predictions), SSR is the Sum of Squares of the Residuals (deviation from the mean of the predictions)
_ = plt.plot(x, y, 'o')
_ = plt.plot(x, res.predict(), 'r-')
_ = plt.axhline(y.mean(), color='k')
for i in range(len(x)):
plt.plot([x[i], x[i]], [y[i], y.mean()], 'y-')
# end for
plt.title('Total Sum of Squares')
_ = plt.plot(x, y, 'o')
_ = plt.plot(x, res.predict(), 'r-')
_ = plt.plot(x, res.predict(), 'r+')
_ = plt.axhline(y.mean(), color='k')
for i in range(len(x)):
plt.plot([x[i], x[i]], [y[i], res.predict()[i]], 'y-')
# end for
plt.title('Error Sum of Squares')
_ = plt.plot(x, y, 'o')
_ = plt.plot(x, res.predict(), 'r-')
_ = plt.plot(x, res.predict(), 'r+')
_ = plt.axhline(y.mean(), color='k')
for i in range(len(x)):
plt.plot(
[x[i], x[i]], [y.mean(), res.predict()[i]], 'y-'
)
# end for
plt.title('Regression Sum of Squares')
data = pd.read_csv('../data/wage.csv')
Explore Data¶
data.head()
_ = plt.plot(data['education'], data['logwage'], 'o')
In order to more clearly see where the data points lie, we add "jitter" so they don't overlap as much. We also adjust alpha
jitter = 0.1
_ = plt.plot(
data['education']
+ np.random.normal(0, 1, len(data['education']))
* jitter,
data['logwage'],
'o',
alpha=0.2,
)
Perform OLS Linear Best Fit¶
res = ols('logwage ~ education', data=data).fit()
Show the line of best fit
jitter = 0.1
_ = plt.plot(
data['education']
+ np.random.normal(0, 1, len(data['education']))
* jitter,
data['logwage'],
'o',
alpha=0.2,
)
data['predict'] = res.predict()
data_std = data.sort_values(by='education')
_ = plt.plot(
data_std['education'], data_std['predict'], 'r-'
)
Compute the Sums of Squares defined above. We use np.ones() to create vectors of constants
n_data = len(data['education'])
SSE = sum(res.resid * res.resid)
sst_array = (
data['logwage']
- np.ones(n_data) * data['logwage'].mean()
)
SST = sum(sst_array * sst_array)
sse_array = (
np.ones(n_data) * data['logwage'].mean() - res.predict()
)
SSR = sum(sse_array * sse_array)
SST, SSE + SSR
Compute the F statistic
F = (SSR / 1) / (SSE / res.df_resid)
F
Check this manually calculated value with the OLS summary
res.summary()
We get the interval that contains 90% of the F distribution function, given the degrees of freedom we have
stats.f.interval(0.90, 1, res.df_resid)
We clearly have a F statistic much larger than this,
so we would reject the Null Hypothesis that there is no linear dependence of logwages on education.
When a Linear Model is Wrong¶
The Anscombe quartet is a famous quartet of datasets, all having the same descriptive numerical values (mean, etc), but wildly different shapes when plotted. We will use the dataset that is a quadratic.
data = pd.read_csv('../data/anscombe.csv')
data.head()
Add some Gaussian noise
data['y2'] = (
data['y2']
+ np.random.normal(0, 1, len(data['y2'])) * 0.45
)
Display the dataset
_ = plt.plot(data['x2'], data['y2'], 'o')
Fit a Linear Model¶
We get the best OLS line, and display it
res = ols('y2 ~ x2', data=data).fit()
_ = plt.plot(data['x2'], data['y2'], 'bo')
_ = plt.plot(data['x2'], res.predict(), 'r-')
When we plot the residuals, we can see that they are apparently not distributed at random (negative at the end, positive in the middle)
_ = plt.plot(data['x2'], res.resid, 'ro')
_ = plt.axhline(0, color='black')
Fit a Quadratic Model¶
We add a term X^2 to our linear model
data['x2sq'] = data['x2'] * data['x2']
data_std = data.sort_values(by='x2')
res2 = ols('y2 ~ x2 + x2sq', data=data_std).fit()
The fit is clearly much better, and the residuals look to be more randomly distributed
_ = plt.plot(data['x2'], data['y2'], 'bo')
_ = plt.plot(data_std['x2'], res2.resid, 'ro')
_ = plt.plot(data_std['x2'], res2.predict(), 'g-')
_ = plt.axhline(0, color='black')
Just plotting the residuals, we see they are smaller than those of the the pure linear model
_ = plt.plot(data['x2'], res2.resid, 'ro')
_ = plt.axhline(0, color='black')
QQ Plots¶
One technique for finding a model that is bad is to look at the QQ Plots for the residuals of eac h OLS fit (linear and quadratic). This is a test for normality; in this case, the difference is very subtle, and could not be used to reject one model or the other
_ = sm.qqplot(res.resid, line='r')
_ = sm.qqplot(res2.resid, line='r')
OLS Summary Comparison¶
A comparison of the RegressionResults summary() output shows that the quadratic model does a much better job of explaining the spread of Y2 values (R Squared of 0.962 versus 0.094). Further, the Jarque-Bera statistic (a test for normality of the residuals) indicates the the residuals are NOT normally distributed in the linear model, and are probably normally distributed in the quadratic model.
The Omnibus statistic also indicates that the residuals are not normal in the pure linear case, and probably are normal in the quadratic case.
res.summary()
res2.summary()
Heteroscedasticity Tests¶
Read and Explore a Medical Dataset¶
data = pd.read_csv('../data/hiv.txt', sep=' ')
data.head()
_ = plt.plot(data['GSS'], data['VL'], 'bo')
Data Cleaning¶
There is clearly an outlier value, so we exclude it.
data2 = data[data['VL'] < 200_000].copy()
len(data), len(data2)
Fit a Linear Model¶
res = ols('VL ~ GSS', data=data2).fit()
Display the line of best fit
_ = plt.plot(data['GSS'], data['VL'], 'bo')
_ = plt.plot(data2['GSS'], res.predict(), 'g-')
Plotting the residuals indicates the the error variance is NOT constant
_ = plt.plot(data2['GSS'], res.resid, 'ro')
_ = plt.axhline(0, color='black')
Looking at the Jarque-Bera statistic in the summary() output, we seec that there is a vanishingly small chance that the residuals have constant variance
res.summary()
Additional Heteroscedasticity Tests¶
statsmodels has a number of tests for constant variance residuals. These tests not addressed in the lecture (but almost certainly R supports them). We show the results for two of them below. They confirm that the residuals are very unlikely to be distributed with constant variance.
x = sm.add_constant(data2['GSS'])
lm, lm_pvalue, fvalue, f_pvalue = het_white(res.resid, x)
print(
'White Test for Heteroscedasticity\n',
f'lagrange multiplier statistic. = {lm}, \n lagrange multiplier statistic = {lm_pvalue},\n',
f'f-statistic of the hypothesis that the error variance does not depend on x = {fvalue}, \n',
f'p-value for the f-statistic = {f_pvalue}',
)
lm, lm_pvalue, fvalue, f_pvalue = het_breuschpagan(
res.resid, x
)
print(
'Breusch-Pagan Test for Heteroscedasticity\n',
f'lagrange multiplier statistic. = {lm}, \n lagrange multiplier statistic = {lm_pvalue},\n',
f'f-statistic of the hypothesis that the error variance does not depend on x = {fvalue}, \n',
f'p-value for the f-statistic = {f_pvalue}',
)
Reproducibility¶
%watermark -h -iv
%watermark