Chapter 7: Problems with the predictors¶
Load the packages:
In [1]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
import statsmodels.api as sm
import statsmodels.formula.api as smf
%precision 7
Out[1]:
Errors in the predictors¶
Load in the data:
In [2]:
cars = pd.read_csv("data/cars.csv")
cars.head()
Out[2]:
Set up the plot and the random seed
In [3]:
fig, ax = plt.subplots()
ax.scatter(cars.speed, cars.dist)
plt.xlabel("Speed")
plt.ylabel("Distance")
xr = np.array(ax.get_xlim())
np.random.seed(123)
Do the four fits with increasing amounts of error. Since we only need univariate regression np.polyfit is adequate.
In [4]:
est = np.polyfit(cars.speed, cars.dist, 1)
ax.plot(xr, est[1] + est[0] * xr)
est1 = np.polyfit(cars.speed + np.random.normal(size=50), cars.dist, 1)
ax.plot(xr, est1[1] + est1[0] * xr, c='red')
est2 = np.polyfit(cars.speed + np.random.normal(scale=2,size=50), cars.dist, 1)
ax.plot(xr, est2[1] + est2[0] * xr, c='green')
est5 = np.polyfit(cars.speed + np.random.normal(scale=5,size=50), cars.dist, 1)
ax.plot(xr, est5[1] + est5[0] * xr, c='cyan')
plt.show()
In [5]:
est, est1, est2, est5
Out[5]:
In [6]:
vv = np.repeat(np.array([0.1, 0.2, 0.3, 0.4, 0.5]), [1000, 1000, 1000, 1000, 1000])
slopes = np.zeros(5000)
for i in range(5000):
slopes[i] = np.polyfit(cars.speed+np.random.normal(scale=np.sqrt(vv[i]), size=50), cars.dist, 1)[0]
In [7]:
betas = np.reshape(slopes, (5, 1000)).mean(axis=1)
betas = np.append(betas,est[0])
variances = np.array([0.6, 0.7, 0.8, 0.9, 1.0, 0.5])
gv = np.polyfit(variances, betas,1)
gv
Out[7]:
In [8]:
plt.scatter(variances, betas)
xr = np.array([0,1])
plt.plot(xr, np.array(gv[1] + gv[0]*xr))
plt.plot([0], [gv[1]], marker='x', markersize=6)
plt.show()
Do not believe there is a simex package for Python.
Changes of Scale¶
Read in the data and fit the model:
In [9]:
savings = pd.read_csv("data/savings.csv",index_col=0)
lmod = smf.ols(formula = 'sr ~ pop15 + pop75 + dpi + ddpi', data=savings).fit()
lmod.summary()
Out[9]:
Rescale one of the predictors:
In [10]:
lmod = smf.ols(formula = 'sr ~ pop15 + pop75 + I(dpi/1000) + ddpi', data=savings).fit()
lmod.summary()
Out[10]:
Standardize all the variables:
In [11]:
scsav = savings.apply(sp.stats.zscore)
lmod = smf.ols(formula = 'sr ~ pop15 + pop75 + dpi + ddpi', data=scsav).fit()
lmod.summary()
Out[11]:
Set up data frame and make a plot of the estimates and confidence intervals:
In [12]:
edf = pd.concat([lmod.params, lmod.conf_int()],axis=1).iloc[1:,]
edf.columns = ['estimate','lb','ub']
npreds = edf.shape[0]
fig, ax = plt.subplots()
ax.scatter(edf.estimate,np.arange(npreds))
for i in range(npreds):
ax.plot([edf.lb[i], edf.ub[i]], [i, i])
ax.set_yticks(np.arange(npreds))
ax.set_yticklabels(edf.index)
ax.axvline(0)
plt.show()
Do Gelman's rescaling:
In [13]:
savings['age'] = np.where(savings.pop15 > 35, 0, 1)
savings['dpis'] = sp.stats.zscore(savings.dpi)/2
savings['ddpis'] = sp.stats.zscore(savings.ddpi)/2
smf.ols(formula='sr ~ age + dpis + ddpis', data=savings).fit().summary()
Out[13]:
Collinearity¶
Read in the data:
In [14]:
seatpos = pd.read_csv("data/seatpos.csv")
seatpos.head()
Out[14]:
Fit a model with all the predictors:
In [15]:
lmod = smf.ols(formula = 'hipcenter ~ Age+Weight+HtShoes+Ht+Seated+Arm+Thigh+Leg', data=seatpos).fit()
lmod.summary()
Out[15]:
Look at the correlations of the predictors:
In [16]:
seatpos.iloc[:,:-1].corr().round(3)
Out[16]:
Do the eigendecomposition. Gives similar but not identical results to R calculation.
In [17]:
X = seatpos.iloc[:,:-1]
XTX = np.dot(X.T,X)
evals, evecs = np.linalg.eig(XTX)
evals
Out[17]:
Calculate the condition numbers;
In [18]:
evals[0]/evals
Out[18]:
Compute the variance inflation factors:
In [19]:
from patsy import dmatrix
X = dmatrix("Age+Weight+HtShoes+Ht+Seated+Arm+Thigh+Leg", data=seatpos, return_type='dataframe')
lmod = sm.OLS(X['Age'],X.drop('Age',axis=1)).fit()
lmod.rsquared
Out[19]:
In [20]:
1/(1-lmod.rsquared)
Out[20]:
Get them all:
In [21]:
from statsmodels.stats.outliers_influence import variance_inflation_factor
vif = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
pd.Series(vif, X.columns)
Out[21]:
Change in coefficients due to perturbing the response
In [22]:
seatpos['hiperb'] = seatpos.hipcenter+np.random.normal(scale=10,size=38)
lmod = smf.ols(formula = 'hipcenter ~ Age+Weight+HtShoes+Ht+Seated+Arm+Thigh+Leg', data=seatpos).fit()
lmodp = smf.ols(formula = 'hiperb ~ Age+Weight+HtShoes+Ht+Seated+Arm+Thigh+Leg', data=seatpos).fit()
pd.DataFrame([lmod.params, lmodp.params])
Out[22]:
Change in R-squared
In [23]:
lmod.rsquared, lmodp.rsquared
Out[23]:
Correlations of the length variables
In [24]:
pd.DataFrame.corr(X.iloc[3:,3:]).round(3)
Out[24]:
Recommended smaller model:
In [25]:
smf.ols(formula = 'hipcenter ~ Age+Weight+Ht', data=seatpos).fit().summary()
Out[25]:
In [26]:
%load_ext version_information
%version_information pandas, numpy, matplotlib, seaborn, scipy, patsy, statsmodels
Out[26]: