Stanford STATS191 in Python, Lecture 5 : Multiple Linear Regression
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 "Multiple linear regression" ( see https://web.stanford.edu/class/stats191/notebooks/Multiple_linear_regression.html )
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
from statsmodels.stats.anova import anova_lm
import os
Read and Explore Data¶
Read a medical dataset relating to prostate cancer
Data details are:
| Variable | Description |
|---|---|
| lcavol | (log) Cancer Volume |
| lweight | (log) Weight |
| age | Patient age |
| lbph | (log) Vening Prostatic Hyperplasia |
| svi | Seminal Vesicle Invasion |
| lcp | (log) Capsular Penetration |
| gleason | Gleason score |
| pgg45 | Percent of Gleason score 4 or 5 |
| lpsa | (log) Prostate Specific Antigen |
| train | Label for test / training split |
data = pd.read_csv('../data/prostate.txt', sep='\t')
data.head()
Display Data¶
We use pandas to give use a quick snapshot of all the variables. We could have removed the "Unnamed" column, but this is just initial exploration. If we look at the lpsa row or column, we see that there is apparently linear correlation with some of the other variables.
_ = pd.plotting.scatter_matrix(data, figsize=(10, 10))
Fit Model¶
The model is of the form of the equation below (linear in all parameters of interest, with Gaussian noise added). Under the Null Hypothesis (H0), all the beta-s are zero.
res = ols(
' lpsa ~ lcavol + lweight + age + lbph + svi + lcp + pgg45',
data=data,
).fit()
We list the model parameter values (the betas), and the probability of seeing each value under H0
print(
f'{"Parameter":30} {"Value":>12} {"Probability (under H0)"}'
)
_ = [
print(f'{na:30} {pa:12.5f} {pr:.5g}')
for na, pa, pr in zip(
res.params.index, res.params.values, res.pvalues
)
]
We sort them by probability
p_srtd = pd.DataFrame()
p_srtd['parameter'] = res.params.index
p_srtd['value'] = res.params.values
p_srtd['probsh0'] = res.pvalues.values
temp = p_srtd.sort_values(by='probsh0')
temp
The estimate of the error signal must take into account the number of parameters we have used in the model. More parameters will make the estimate larger for a given Sum Of Squared Errors (SSE). We show how the RegressionResults attribute df_resid is calculated
print(f'{res.df_resid} = {len(data)} - {res.df_model} - 1')
Compute our estimate of the standard deviation of the error noise.
sigma_hat = np.sqrt(
sum(res.resid * res.resid) / res.df_resid
)
sigma_hat
We can check this against the value in the RegressionResults object
np.sqrt(res.mse_resid)
Display Summary of Regression Results¶
We display the results of the regression by OLS, for two different values of alpha (that sets the probability bounds for which we want confidence intervals). The default value for statsmodels is 0.05. The lecture series we are following has alpha at 0.1
Statsmodels default¶
res.summary()
Lecture series default for alpha¶
alpha is the probability used to determine confidence intervals of parameter estimates
res.summary(alpha=0.1)
We display the confidence intervals at the 90% level for the parameters we have estimated
res.conf_int(alpha=0.1)
Explicit Check of Parameter Values¶
We can check (and better understand the summary) by explicitly computing the parameter estimates
We convert pandas data structures into numpy data structures, to use the linear algebra functions of numpy. We have to add a column of ones, to get the intercept
data['ones'] = np.ones(len(data))
X = data[
[
'ones',
'lcavol',
'lweight',
'age',
'lbph',
'svi',
'lcp',
'pgg45',
]
]
y = data['lpsa'].to_numpy()
X.to_numpy()
beta_hat = np.linalg.inv(X.T @ X) @ X.T @ y.T
If we compare the OLS values with the explicit values, we see they match
print(f'{"Linear Alg Value":>20} {"OLS Value":>20}')
_ = [
print(f'{p1:20.5f}, {p2:20.5f}')
for p1, p2 in zip(beta_hat, res.params.values)
]
Adjusted R^2¶
We compute the R^2 values explicitly, and show they match values in the RegressionResults object. Note there seems to be a difference between the abbreviations used in the lectures, and used in the RegressionResults object. In my code, I will use the STATS191 lecture conventions (even though I think the statsmodels conventions make more sense)
SST = sum(
(data['lpsa'] - data['lpsa'].mean())
* (data['lpsa'] - data['lpsa'].mean())
)
print(
f' Explicit Calculation = {SST:15.5f},\n Statsmodels Value = {res.centered_tss:15.5f}'
)
SSR = sum(
(res.predict() - data['lpsa'].mean())
* (res.predict() - data['lpsa'].mean())
)
SSR, res.ess
print(
f' Explicit Calculation = {SSR:15.5f},\n Statsmodels Value = {res.ess:15.5f}'
)
SSE = sum(
(res.predict() - data['lpsa'])
* (res.predict() - data['lpsa'])
)
SSE, res.ssr
print(
f' Explicit Calculation = {SSE:15.5f},\n Statsmodels Value = {res.ssr:15.5f}'
)
r_2 = SSR / SST
r_2, res.rsquared
print(
f' Explicit Calculation = {r_2:15.5f},\n Statsmodels Value = {res.rsquared:15.5f}'
)
If we add more and more explanatory variables, even random ones, R^2 will get smaller. We adjust for this
We compute R2_Adjusted in three different ways; as expected, all three values match
r_adjust_1 = 1 - res.mse_resid / res.mse_total
r_adjust_2 = 1 - (SSE / (len(data) - len(res.params))) / (
SST / (len(data) - 1)
)
print(
f' Explicit Calculation = {r_adjust_2:15.5f},\n Explicit Calculation2= {r_adjust_1:15.5f},\n Statsmodels Value = {res.rsquared_adj:15.5f}'
)
Goodness of Fit¶
We form the F statistic by computing the ratio of:
mean square difference of the predictions from the mean of the dependent variable
to
mean square difference of the predictions from the actual values of the dependent variable
If the Null Hypothesis (H0) holds (no relationship between explanatory and dependent variable), then this should be close to 1.0. The larger the F statistic is, the more the variation of the dependent variable is explained by the explanatory variables.
When we ran our OLS fit, the F statistic was computed, and returned in the RegressionResults object. We check its value, and the RegressionResults object is consistent
F = res.mse_model / res.mse_resid
F2 = (SSR / (len(res.params) - 1)) / (
SSE / (len(data) - len(res.params))
)
print(
f' Explicit Calculation = {F:15.5f},\n Explicit Calculation2= {F2:15.5f},\n Statsmodels Value = {res.fvalue:15.5f}'
)
Compute the degrees of freedom of the numerator and denominator of the F ratio
dfn = len(res.params) - 1
dfd = len(data['lpsa']) - dfn - 1
dfn, dfd
Display a plot of the F probability density function, and the value of F we actually see
x = np.linspace(0.01, 26, 201)
y = stats.f.pdf(x, dfn, dfd)
_ = plt.plot(x, y, 'r-', label='F Prob Dist')
_ = plt.axvline(F2, color='green', label='F Value Seen')
_ = plt.axhline(0, color='black')
_ = plt.legend(loc='best')
_ = plt.grid()
We use the survival function sf rather than 1-cdf(), as it is supposed to be more accurate. Note that the value from the equivalent R function gives a different answer, but this to be expected, given the very small value returned.
Clearly, the F value we see is very small, and so we reject H0.
stats.f.sf(res.fvalue, res.df_model, res.df_resid)
Prediction¶
To get a prediction of the value of lpsa for a given set of values of the explanatory variables in our model, we can either form a vector of values (including 1 for the constant or intercept term), and form the dot product with our vector of parameter values (beta). Alternatively, we can form a dict given the variable values by name, and call the predict() method of the RegressionResults method
a = [1, 1.3, 3.6, 64, 0.1, 0.2, -0.2, 25]
a_dict = {
'lcavol': 1.3,
'lweight': 3.6,
'age': 64,
'lbph': 0.1,
'svi': 0.2,
'lcp': -0.2,
'pgg45': 25,
}
Form the dot product
y_hat = sum(res.params.values * a)
y_hat
Call the predict method, as expected, the answers are the same
res.predict(a_dict)
We can get the confidence intervals for our prediction by the get_prediction method. Note that we can be more confident about the mean of our predictions, than about any one prediction (which must include the error term of our model)
gp = res.get_prediction(a_dict)
gp.summary_frame(alpha=0.05)
Comparing models¶
Having build and fit our model, we may believe that the true parameter values for lcavol and svi are the same.
One way to test this is to use the statsmodels formula constraint test (obvously not covered in the STATS191 lecture)
to_test = 'lcavol = svi'
t_test = res.t_test(to_test)
print(t_test)
We cannot reject the hypothesis the the parameters are the same
A different way to analyse the equivalence of parameters, is to form a different model that has a built in constraint that parameter values for lcavol and svi are the same. We do this by creating a new variable, being the sum of the two variables lcavol and svi, and including that is our new model
data['z'] = data['lcavol'] + data['svi']
res2 = ols(
'lpsa ~ z + lweight + age + lbph + lcp + pgg45',
data=data,
).fit()
We now perform a Analysis of Variance on the two models, to see if one has a markedly better explanation of the variance of the dependent variable lpsa. It appears that the difference we see could arise about 50% of the time, so we are unable to say which model is better. Note we get some spurious diagnostics; however the values match those in the STATS191 lecture
| F | df.num | df.den | p.value |
|---|---|---|---|
| 0.4818657 | 1 | 89 | 0.4893864 |
anova_lm(res2, res)
The final example involved testing if lcavol and svi parameters sum to 1.0
The first way to do this is by the statsmodels formula constraint test. Again, we find that we cannot get markedly better variance explanation by this constrained model
to_test = 'lcavol + svi = 1'
t_test = res.t_test(to_test)
print(t_test)
Again at the 5% level, we cannot say which model is better
We can build a new model that enforces our constraint, fit it by OLS, and compare the variance via the F test. The original model is below
The new model is as below. We form new variable z2 and y2
data['z2'] = data['lcavol'] - data['svi']
data['y2'] = data['lpsa'] - data['svi']
res3 = ols(
'y2 ~ z2 + lweight + age + lbph + lcp + pgg45',
data=data,
).fit()
Perform ANOVA
anova_lm(res3, res)
Again, we cannot reject the difference between models are being more than random chance.
Reproducibility¶
%watermark -h -iv
%watermark