Chapter 14: Categorical Predictors¶
Import 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-level factor¶
Read in the data
In [2]:
sexab = pd.read_csv("data/sexab.csv", index_col=0)
sexab.head()
Out[2]:
Construct something similar to the by function summary from LMR.
In [3]:
lfuncs = ['min','median','max']
sexab.groupby('csa').agg({'cpa': lfuncs,'ptsd': lfuncs}).round(1)
Out[3]:
In [4]:
sexab.boxplot('ptsd',by='csa')
plt.suptitle("")
plt.show()
In [5]:
sns.pairplot(x_vars="cpa", y_vars="ptsd", data=sexab, hue="csa",size=5)
plt.show()
In [6]:
stats.ttest_ind(sexab.ptsd[sexab.csa == 'Abused'], sexab.ptsd[sexab.csa == 'NotAbused'])
Out[6]:
Define dummy variables and including an intercept does fit a model but the meaning of the coefficients and standard errors are problematic. Look out for the design matrix is singular warning at the end.
In [7]:
df1 = (sexab.csa == 'Abused').astype(int)
df2 = (sexab.csa == 'NotAbused').astype(int)
X = np.column_stack((df1,df2))
lmod = sm.OLS(sexab.ptsd,sm.add_constant(X)).fit()
lmod.summary()
Out[7]:
Use only one of the dummy variables:
In [8]:
lmod = sm.OLS(sexab.ptsd,sm.add_constant(df2)).fit()
lmod.summary()
Out[8]:
Or we can just not have an intercept term:
In [9]:
lmod = sm.OLS(sexab.ptsd,X).fit()
lmod.summary()
Out[9]:
Using the formula method produces R-like results
In [10]:
lmod = smf.ols('ptsd ~ csa', sexab).fit()
lmod.summary()
Out[10]:
Can view the design matrix: (picking out the first two rows of each level of csa).
In [11]:
import patsy
selcols = [0,1,45,46]
patsy.dmatrix('~ csa', sexab)[selcols,:]
Out[11]:
Can check the types of variables:
In [12]:
sexab.csa.dtype, sexab.ptsd.dtype
Out[12]:
pandas has a way to create the dummy variables:
In [13]:
sac = pd.concat([sexab,pd.get_dummies(sexab.csa)],axis=1)
lmod = smf.ols('ptsd ~ Abused', sac).fit()
lmod.summary()
Out[13]:
Alternatively, we can set the reference level to NotAbused:
In [14]:
lmod = smf.ols('ptsd ~ C(csa,Treatment(reference="NotAbused"))', sexab).fit()
lmod.summary()
Out[14]:
Factors and Quantitative predictors¶
In [15]:
lmod4 = smf.ols('ptsd ~ csa*cpa', sexab).fit()
lmod4.summary()
Out[15]:
Take a look at a chunk from the design matrix:
In [16]:
patsy.dmatrix('~ csa*cpa', sexab)[40:50,]
Out[16]:
Plot the fit by group - seem like a lot of work!
In [17]:
abused = (sexab.csa == "Abused")
plt.scatter(sexab.cpa[abused], sexab.ptsd[abused], marker='x',label="abused")
xl,xu = [-3, 9]
a, b = (lmod4.params[0], lmod4.params[2])
plt.plot([xl,xu], [a+xl*b,a+xu*b])
plt.scatter(sexab.cpa[~abused], sexab.ptsd[~abused], marker='o',label="not abused")
a, b = (lmod4.params[0]+lmod4.params[1], lmod4.params[2]+lmod4.params[3])
plt.plot([xl,xu], [a+xl*b,a+xu*b])
plt.legend()
plt.show()
Can produce essentially the same plot but the regression lines are fit independently to the groups. In this case, the fits will be identical.
In [18]:
sns.lmplot(x="cpa", y="ptsd", hue="csa", data=sexab, ci=None)
plt.show()
In [19]:
lmod3 = smf.ols('ptsd ~ csa+cpa', sexab).fit()
lmod3.summary()
Out[19]:
No shortcut to producing the plot this time.
In [20]:
abused = (sexab.csa == "Abused")
plt.scatter(sexab.cpa[abused], sexab.ptsd[abused], marker='x',label="abused")
xl,xu = [-3, 9]
a, b = (lmod3.params[0], lmod3.params[2])
plt.plot([xl,xu], [a+xl*b,a+xu*b])
plt.scatter(sexab.cpa[~abused], sexab.ptsd[~abused], marker='o',label="not abused")
a, b = (lmod3.params[0]+lmod4.params[1], lmod3.params[2])
plt.plot([xl,xu], [a+xl*b,a+xu*b])
plt.legend()
plt.show()
Get the confidence intervals:
In [21]:
lmod3.conf_int()
Out[21]:
In [22]:
sns.residplot(lmod3.fittedvalues, lmod3.resid, lowess=True)
plt.xlabel("Fitted values")
plt.ylabel("Residuals")
plt.show()
In [23]:
lmod1 = smf.ols('ptsd ~ cpa', sexab).fit()
lmod1.summary()
Out[23]:
Interpretation with interaction terms¶
In [24]:
ws = pd.read_csv("data/whiteside.csv", index_col=0)
ws.head()
Out[24]:
In [25]:
sns.lmplot(x="Temp", y="Gas", col="Insul", data=ws)
plt.show()
In [26]:
lmod = smf.ols('Gas ~ Temp*Insul', ws).fit()
lmod.summary()
Out[26]:
In [27]:
ws.Temp.mean()
Out[27]:
In [28]:
ws['cTemp'] = ws.Temp - ws.Temp.mean()
lmod = smf.ols('Gas ~ cTemp*Insul', ws).fit()
lmod.summary()
Out[28]:
Factors with more than two levels¶
In [29]:
ff = pd.read_csv("data/fruitfly.csv", index_col=0)
ff.head()
Out[29]:
In [30]:
sns.pairplot(x_vars="thorax", y_vars="longevity", data=ff, hue="activity",size=5)
plt.show()
In [31]:
sns.lmplot(x="thorax", y="longevity", data=ff, col="activity",size=5)
plt.show()
Fit the model
In [32]:
lmod = smf.ols('longevity ~ thorax*activity', ff).fit()
lmod.summary()
Out[32]:
Construct and show selected rows of the design matrix (one for each level of activity). DataFrame is just for pretty printing.
In [33]:
mm = patsy.dmatrix('~ thorax*activity', ff)
ii = (1, 25, 49, 75, 99)
pd.DataFrame(mm[ii,:],index=ii,columns=lmod.params.index)
Out[33]:
In [34]:
sns.residplot(lmod.fittedvalues, lmod.resid, lowess=True)
plt.xlabel("Fitted values")
plt.ylabel("Residuals")
plt.show()
In [35]:
sm.qqplot(lmod.resid, line="q")
Out[35]:
In [36]:
sm.stats.anova_lm(lmod)
Out[36]:
In [37]:
lmods = smf.ols('longevity ~ thorax+activity', ff).fit()
sm.stats.anova_lm(lmods,lmod)
Out[37]:
In [38]:
lmod = smf.ols('np.log(longevity) ~ thorax+activity', ff).fit()
lmod.summary()
Out[38]:
In [39]:
sns.residplot(lmod.fittedvalues, lmod.resid, lowess=True)
plt.xlabel("Fitted values")
plt.ylabel("Residuals")
plt.show()
In [40]:
np.exp(lmod.params[1:5])
Out[40]:
In [41]:
lmod = smf.ols('longevity ~ activity', ff).fit()
lmod.summary()
Out[41]:
In [42]:
sm.stats.anova_lm(lmod)
Out[42]:
In [43]:
lmod = smf.ols('np.log(longevity) ~ activity', ff).fit()
lmod.summary()
Out[43]:
Factor coding¶
In [44]:
from patsy.contrasts import Treatment
levels = [1,2,3,4]
contrast = Treatment(reference=0).code_without_intercept(levels)
print(contrast.matrix)
In [45]:
from patsy.contrasts import Sum
contrast = Sum().code_without_intercept(levels)
print(contrast.matrix)
In [46]:
from patsy.contrasts import Helmert
contrast = Helmert().code_without_intercept(levels)
print(contrast.matrix)
In [47]:
lmod = smf.ols('ptsd ~ C(csa,Sum)', sexab).fit()
lmod.summary()
Out[47]:
In [48]:
%load_ext version_information
%version_information pandas, numpy, matplotlib, seaborn, scipy, patsy, statsmodels
Out[48]: