Chapter 15: One Factor Models¶
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
Coagulation example¶
In [2]:
coagulation = pd.read_csv("data/coagulation.csv", index_col=0)
coagulation.head()
Out[2]:
pandas version of boxplot
In [3]:
coagulation.boxplot('coag',by='diet')
plt.suptitle("")
plt.show()
seaborn version of boxplot
In [4]:
sns.boxplot(x="diet", y="coag", data=coagulation)
plt.show()
Maybe seaborn plot is preferable (although without the colors). Pandas version is annoying over the title and failure to label the y-axis.
swarmplot is good when there are replicated points:
In [5]:
sns.swarmplot(x="diet", y="coag", data=coagulation)
plt.show()
Fit the model:
In [6]:
lmod = smf.ols("coag ~ diet", coagulation).fit()
lmod.summary()
Out[6]:
Examine the coefficients:
In [7]:
lmod.params.round(1)
Out[7]:
Look at the design matrix:
In [8]:
import patsy
patsy.dmatrix('~ diet', coagulation)
Out[8]:
Test the significance of diet.
In [9]:
sm.stats.anova_lm(lmod)
Out[9]:
Fit without an intercept term:
In [10]:
lmodi = smf.ols("coag ~ diet-1", coagulation).fit()
lmodi.summary()
Out[10]:
Do the ANOVA explicitly:
In [11]:
%%capture --no-display
lmodnull = smf.ols("coag ~ 1", coagulation).fit()
sm.stats.anova_lm(lmodnull, lmod)
Out[11]:
Do it with sum contrasts:
In [12]:
from patsy.contrasts import Sum
lmods = smf.ols("coag ~ C(diet,Sum)", coagulation).fit()
lmods.summary()
Out[12]:
Diagnostics¶
In [13]:
sns.residplot(lmod.fittedvalues, lmod.resid)
plt.xlabel("Fitted values")
plt.ylabel("Residuals")
plt.show()
In [14]:
fig = sm.qqplot(lmod.resid, line="q")
plt.show()
Do the Levene test:
In [15]:
coagulation['meds'] = coagulation.groupby('diet').transform(np.median)
coagulation['mads'] = abs(coagulation.coag - coagulation.meds)
coagulation.head()
Out[15]:
In [16]:
lmodb = smf.ols('mads ~ diet', coagulation).fit()
sm.stats.anova_lm(lmodb)
Out[16]:
In [17]:
coagulation.coag[coagulation.diet == "A"]
Out[17]:
Do the Bartlett test:
In [18]:
stats.bartlett(coagulation.coag[coagulation.diet == "A"],
coagulation.coag[coagulation.diet == "B"],
coagulation.coag[coagulation.diet == "C"],
coagulation.coag[coagulation.diet == "D"])
Out[18]:
Pairwise comparisons¶
Assemble the pieces necessary for pairwise comparisons:
In [19]:
lmod.bse
Out[19]:
In [20]:
lmod.params
Out[20]:
In [21]:
stats.t.ppf(0.975,20)
Out[21]:
In [22]:
lmod.params[1] + np.array([-1, 1]) * stats.t.ppf(0.975,20) * lmod.bse[1]
Out[22]:
Can do the Tukey pairwise:
In [23]:
thsd = sm.stats.multicomp.pairwise_tukeyhsd(coagulation.coag, coagulation.diet)
thsd.summary()
Out[23]:
In [24]:
fig = thsd.plot_simultaneous()
plt.show()
False discovery rate¶
Read in the data:
In [25]:
jsp = pd.read_csv("data/jsp.csv", index_col=0)
jsp.head()
Out[25]:
In [26]:
jsp['mathcent'] = jsp.math - np.mean(jsp.math)
Try two different plotting methods
In [27]:
jsp.boxplot('mathcent',by='school')
plt.suptitle("")
plt.xticks(fontsize=8, rotation=90)
plt.show()
In [28]:
sns.boxplot(x="school", y="mathcent", data=jsp)
plt.xticks(fontsize=8, rotation=90)
plt.show()
In [29]:
lmod = smf.ols("mathcent ~ C(school) - 1", jsp).fit()
lmod.summary()
Out[29]:
In [30]:
sm.stats.anova_lm(smf.ols("mathcent ~ C(school)", jsp).fit())
Out[30]:
In [31]:
lmod.pvalues
Out[31]:
In [32]:
from statsmodels.sandbox.stats.multicomp import multipletests
reject, padj, _, _ = multipletests(lmod.pvalues, method="bonferroni")
These schools are rejected by Bonferroni.
In [33]:
lmod.params[reject]
Out[33]:
In [34]:
selsch = np.argsort(lmod.pvalues)[np.sort(lmod.pvalues) < np.arange(1,50)*0.05/49]
lmod.params.index[selsch]
Out[34]:
Schools identified using false discovery rate.
In [35]:
reject, padj, _, _ = multipletests(lmod.pvalues, method="fdr_bh")
lmod.params[reject]
Out[35]:
In [36]:
%load_ext version_information
%version_information pandas, numpy, matplotlib, seaborn, scipy, patsy, statsmodels
Out[36]: