Chapter 16: Models with several factors¶
Load the packages
In [1]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels.formula.api as smf
import seaborn as sns
from scipy import stats
Two factors with no replication¶
Load in the data
In [2]:
composite = pd.read_csv("data/composite.csv", index_col=0)
composite
Out[2]:
Create the interaction plots:
In [3]:
sns.catplot(x="laser", y="strength", hue="tape", data=composite, kind="point")
Out[3]:
In [4]:
sns.catplot(x="tape", y="strength", hue="laser", data=composite, kind="point")
Out[4]:
Fit a saturated model. Don't show the error messages (get complaints about small datasets)
In [5]:
%%capture --no-display
lmod = smf.ols("strength ~ tape * laser", composite).fit()
lmod.summary()
Out[5]:
Fit main effects only model:
In [6]:
lmod = smf.ols("strength ~ tape + laser", composite).fit()
lmod.params
Out[6]:
In [7]:
lasercoefs = np.tile([0, 6.57, 12.21],3)
lasercoefs
Out[7]:
In [8]:
tapecoefs = np.repeat([1.48, 5.52, 0], 3)
tapecoefs
Out[8]:
In [9]:
%%capture --no-display
composite['crossp'] = tapecoefs * lasercoefs
tmod = smf.ols("strength ~ tape + laser + crossp", composite).fit()
tmod.summary()
Out[9]:
In [10]:
sm.stats.anova_lm(tmod)
Out[10]:
Fix the order of the categories
In [11]:
cat_type = pd.api.types.CategoricalDtype(categories=['slow','medium','fast'],ordered=True)
composite['tape'] = composite.tape.astype(cat_type)
In [12]:
%%capture --no-display
from patsy.contrasts import Poly
lmod = smf.ols("strength ~ C(tape,Poly) + C(laser,Poly)", composite).fit()
lmod.summary()
Out[12]:
Check out the design matrix
In [13]:
import patsy
dm = patsy.dmatrix('~ C(tape,Poly) + C(laser,Poly)', composite)
np.asarray(dm).round(2)
Out[13]:
In [14]:
composite['Nlaser'] = np.tile([40,50,60],3)
composite['Ntape'] = np.repeat([6.42,13,27], 3)
In [15]:
%%capture --no-display
lmodn = smf.ols("strength ~ np.log(Ntape) + I(np.log(Ntape)**2) + Nlaser", composite).fit()
lmodn.summary()
Out[15]:
Two factors with replication¶
Read in the data while making sure the predictors are treated as categorical variables:
In [16]:
pvc = pd.read_csv("data/pvc.csv", index_col=0, dtype={'psize':'float','operator':'category','resin':'category'})
pvc.head()
Out[16]:
Compute the group means
In [17]:
pvc.groupby('operator')[['psize']].mean()
Out[17]:
In [18]:
sns.catplot(x='resin',y='psize',hue='operator', data=pvc, kind="point",ci=None)
Out[18]:
In [19]:
ax = sns.catplot(x='operator',y='psize',hue='resin', data=pvc, ci=None,scale=0.5,kind="point")
ax = sns.swarmplot(x='operator',y='psize',hue='resin', data=pvc)
ax.legend_.remove()
Create the ANOVA table
In [20]:
lmod = smf.ols('psize ~ operator*resin', pvc).fit()
sm.stats.anova_lm(lmod).round(3)
Out[20]:
Now remove the interaction term:
In [21]:
lmod = smf.ols('psize ~ operator+resin', pvc).fit()
sm.stats.anova_lm(lmod).round(3)
Out[21]:
Make the diagnostic plots:
In [22]:
fig = sm.qqplot(lmod.resid, line="q")
In [23]:
sns.residplot(lmod.fittedvalues, lmod.resid)
plt.xlabel("Fitted values")
plt.ylabel("Residuals")
plt.show()
In [24]:
sns.swarmplot(pvc['operator'], lmod.resid)
plt.axhline(0)
plt.xlabel("Operator")
plt.ylabel("Residuals")
plt.show()
Test for equality of variance:
In [25]:
lmod = smf.ols('psize ~ operator*resin', pvc).fit()
ii = np.arange(0,48,2)
pvce = pvc.iloc[ii,:].copy()
pvce['res'] = np.sqrt(abs(lmod.resid.loc[ii+1]))
vmod = smf.ols('res ~ operator+resin', pvce).fit()
sm.stats.anova_lm(vmod,typ=3).round(3)
Out[25]:
In [26]:
lmod = smf.ols('psize ~ operator+resin', pvc).fit()
lmod.summary()
Out[26]:
Get the Tukey critical value:
In [27]:
from statsmodels.sandbox.stats.multicomp import get_tukeyQcrit
get_tukeyQcrit(3,38)
Out[27]:
Get the first pairwise difference:
In [28]:
lmod.params[1] + np.array([-1,1]) * get_tukeyQcrit(3,38) * lmod.bse[1] / np.sqrt(2)
Out[28]:
Two factors with an interaction¶
Read in the data:
In [29]:
warpbreaks = pd.read_csv("data/warpbreaks.csv", index_col=0, dtype={'breaks':'int','wool':'category','tension':'category'})
warpbreaks.head()
Out[29]:
In [30]:
sns.boxplot(x="wool", y="breaks", data=warpbreaks)
plt.show()
In [31]:
sns.catplot(x='wool',y='breaks',hue='tension', data=warpbreaks, ci=None,kind="point")
Out[31]:
In [32]:
warpbreaks['nwool'] = warpbreaks['wool'].cat.codes
In [33]:
ax = sns.catplot(x='wool',y='breaks',hue='tension', data=warpbreaks, ci=None,kind="point")
ax = sns.swarmplot(x='wool',y='breaks',hue='tension', data=warpbreaks)
ax.legend_.remove()
In [34]:
ax = sns.catplot(x='tension',y='breaks',hue='wool', data=warpbreaks, ci=None, kind="point")
ax = sns.swarmplot(x='tension',y='breaks',hue='wool', data=warpbreaks)
ax.legend_.remove()
In [35]:
lmod = smf.ols('breaks ~ wool * tension', warpbreaks).fit()
sns.regplot(lmod.fittedvalues, lmod.resid,scatter_kws={'alpha':0.3}, ci=None)
plt.xlabel("Fitted values")
plt.ylabel("Residuals")
Out[35]:
In [36]:
lmod = smf.ols('np.sqrt(breaks) ~ wool * tension', warpbreaks).fit()
sns.regplot(lmod.fittedvalues, lmod.resid,scatter_kws={'alpha':0.3}, ci=None)
plt.xlabel("Fitted values")
plt.ylabel("Residuals")
Out[36]:
In [37]:
sm.stats.anova_lm(lmod)
Out[37]:
In [38]:
lmod.summary()
Out[38]:
In [39]:
lmod = smf.ols('np.sqrt(breaks) ~ wool:tension-1', warpbreaks).fit()
lmod.summary()
Out[39]:
Width of Tukey HSD bands is ( qt/sqrt(2) se sqrt(2))
In [40]:
get_tukeyQcrit(6,48) * lmod.bse[0]
Out[40]:
Compute all the pairwise differences:
In [41]:
import itertools
dp = set(itertools.combinations(range(0,6),2))
dcoef = []
namdiff = []
for cp in dp:
dcoef.append(lmod.params[cp[0]] - lmod.params[cp[1]])
namdiff.append(lmod.params.index[cp[0]] + '-' + lmod.params.index[cp[1]])
thsd = pd.DataFrame({'Difference':dcoef},index=namdiff)
thsd["lb"] = thsd.Difference - get_tukeyQcrit(6,48) * lmod.bse[0]
thsd["ub"] = thsd.Difference + get_tukeyQcrit(6,48) * lmod.bse[0]
thsd
Out[41]:
In [42]:
lmod.params.index
Out[42]:
Larger factorial experiments¶
Read in the data
In [43]:
speedo = pd.read_csv("data/speedo.csv", index_col=0)
speedo
Out[43]:
In [44]:
%%capture --no-display
lmod = smf.ols('y ~ h+d+l+b+j+f+n+a+i+e+m+c+k+g+o', speedo).fit()
lmod.summary()
Out[44]:
Look at the design matrix:
In [45]:
dm = patsy.dmatrix('~ h+d+l+b+j+f+n+a+i+e+m+c+k+g+o', speedo)
np.asarray(dm)
Out[45]:
Make a QQ plot of the coefficients:
In [46]:
ii = np.argsort(lmod.params[1:])
scoef = np.array(lmod.params[1:])[ii]
lcoef = (speedo.columns[:-1])[ii]
n = len(scoef)
qq = stats.norm.ppf(np.arange(1,n+1)/(n+1))
fig, ax = plt.subplots()
ax.scatter(qq, scoef,s=1)
for i in range(len(qq)):
ax.annotate(lcoef[i], (qq[i],scoef[i]))
plt.xlabel("Normal Quantiles")
plt.ylabel("Sorted coefficients")
plt.show()
In [47]:
n = len(scoef)
qq = stats.norm.ppf((n + np.arange(1,n+1))/(2*n+1))
acoef = np.abs(lmod.params[1:])
ii = np.argsort(acoef)
acoef = acoef[ii]
lcoef = (speedo.columns[:-1])[ii]
fig, ax = plt.subplots()
ax.scatter(qq, acoef,s=1)
for i in range(len(qq)):
ax.annotate(lcoef[i], (qq[i],acoef[i]))
plt.xlabel("Normal Half Quantiles")
plt.ylabel("Sorted absolute coefficients")
plt.show()
In [48]:
%load_ext version_information
%version_information pandas, numpy, matplotlib, seaborn, scipy, patsy, statsmodels
Out[48]: