Chapter 17: Experiments with Blocks¶
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
Randomised block design¶
In [2]:
oatvar = pd.read_csv("data/oatvar.csv", index_col=0, dtype={'yield':'int','block':'category','variety':'category'})
oatvar.head()
Out[2]:
Display the data as a table:
In [3]:
oatvar.pivot(index = 'variety', columns='block', values='yield')
Out[3]:
Make plots of the data:
In [4]:
sns.boxplot(x="variety", y="yield", data=oatvar)
plt.show()
In [5]:
sns.boxplot(x="block", y="yield", data=oatvar)
plt.show()
In [6]:
sns.catplot(x='variety',y='yield',hue='block', data=oatvar, kind='point')
plt.show()
In [7]:
sns.catplot(x='block',y='yield',hue='variety', data=oatvar, kind='point')
plt.show()
The obvious way to fit the model fails, apparently because patsy doesn't like the response to be called yield.
In [8]:
try:
lmod = smf.ols('yield ~ variety + block', oatvar).fit()
lmod.summary()
except:
pass
Look at the design matrix:
In [9]:
import patsy
dm = patsy.dmatrix('~ variety + block', oatvar)
np.asarray(dm)
Out[9]:
In [10]:
lmod = sm.OLS(oatvar['yield'], dm).fit()
lmod.summary()
Out[10]:
Another solution is to create a response with a new name:
In [11]:
oatvar['y'] = oatvar['yield']
lmod = smf.ols('y ~ variety + block', oatvar).fit()
sm.stats.anova_lm(lmod).round(3)
Out[11]:
Fit with only block
In [12]:
lmod = smf.ols('y ~ block', oatvar).fit()
sm.stats.anova_lm(lmod).round(3)
Out[12]:
Leave out the first observation:
In [13]:
oat1 = oatvar.iloc[1:,]
lmod = smf.ols('y ~ variety + block', oat1).fit()
sm.stats.anova_lm(lmod).round(3)
Out[13]:
Change the factor order:
In [14]:
lmod = smf.ols('y ~ block + variety', oat1).fit()
sm.stats.anova_lm(lmod).round(3)
Out[14]:
Or use type III ANOVA which is like drop() in R
In [15]:
sm.stats.anova_lm(lmod,typ=3)
Out[15]:
In [16]:
lmod = smf.ols('y ~ variety + block', oatvar).fit()
sns.residplot(lmod.fittedvalues, lmod.resid)
plt.xlabel("Fitted values")
plt.ylabel("Residuals")
plt.show()
In [17]:
res = np.sort(lmod.resid)
n = len(res)
qq = stats.norm.ppf(np.arange(1,n+1)/(n+1))
fig, ax = plt.subplots()
ax.scatter(qq, res)
plt.xlabel("Normal Quantiles")
plt.ylabel("Sorted Residuals")
plt.show()
Do the interaction test:
In [18]:
lmod.params
Out[18]:
In [19]:
varcoefs = np.append([0.],lmod.params[1:8])
varcoefs = np.repeat(varcoefs,5)
blockcoefs = np.append([0.],lmod.params[8:12])
blockcoefs = np.tile(blockcoefs,8)
In [20]:
oatvar['crossp'] = varcoefs * blockcoefs
tmod = smf.ols("y ~ variety + block + crossp", oatvar).fit()
tmod.summary()
Out[20]:
Relative efficiency
In [21]:
lmcrd = smf.ols("y ~ variety", oatvar).fit()
np.sqrt(lmcrd.scale)
Out[21]:
In [22]:
np.sqrt(lmod.scale)
Out[22]:
In [23]:
lmcrd.scale/lmod.scale
Out[23]:
Latin square¶
In [24]:
abrasion = pd.read_csv("data/abrasion.csv", index_col=0,
dtype={'run':'category','position':'category','material':'category','wear':'int'})
abrasion.head()
Out[24]:
See the Latin square layout:
In [25]:
abrasion.pivot(index = 'run', columns='position', values='wear')
Out[25]:
In [26]:
abrasion.pivot(index = 'run', columns='position', values='material')
Out[26]:
In [27]:
sns.catplot(x='run',y='wear',hue='position', data=abrasion, kind='point',scale=0.5)
Out[27]:
In [28]:
sns.catplot(x='run',y='wear',hue='material', data=abrasion, kind='point',scale=0.5)
Out[28]:
In [29]:
lmod = smf.ols('wear ~ run + position + material', abrasion).fit()
sm.stats.anova_lm(lmod,typ=3).round(3)
Out[29]:
In [30]:
%%capture --no-display
lmod.summary()
Out[30]:
Do the Tukey computation:
In [31]:
from statsmodels.sandbox.stats.multicomp import get_tukeyQcrit
get_tukeyQcrit(4,6)*lmod.bse[1]/np.sqrt(2)
Out[31]:
In [32]:
mcoefs = np.append([0],lmod.params[7:10])
nmats = ['A','B','C','D']
In [33]:
import itertools
dp = set(itertools.combinations(range(0,4),2))
dcoef = []
namdiff = []
for cp in dp:
dcoef.append(mcoefs[cp[0]] - mcoefs[cp[1]])
namdiff.append(nmats[cp[0]] + '-' + nmats[cp[1]])
thsd = pd.DataFrame({'Difference':dcoef},index=namdiff)
thsd["lb"] = thsd.Difference - get_tukeyQcrit(4,6) * lmod.bse[1]/np.sqrt(2)
thsd["ub"] = thsd.Difference + get_tukeyQcrit(4,6) * lmod.bse[1]/np.sqrt(2)
thsd.round(2)
Out[33]:
Relative efficiency
In [34]:
lmodr = smf.ols('wear ~ material', abrasion).fit()
lmodr.scale/lmod.scale
Out[34]:
Balanced incomplete block design¶
In [35]:
rabbit = pd.read_csv("data/rabbit.csv", index_col=0,
dtype={'treat':'category','gain':'float','block':'category'})
rabbit.head()
Out[35]:
Show the arrange of the BIB. I replace the NaNs with spaces for clarity.
In [36]:
pt = rabbit.pivot(index = 'treat', columns='block', values='gain')
pt.replace(np.nan," ", regex=True)
Out[36]:
In [37]:
sns.swarmplot(x='block',y='gain',hue='treat',data=rabbit)
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()
In [38]:
sns.swarmplot(x='treat',y='gain',hue='block',data=rabbit)
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()
In [39]:
lmod = smf.ols('gain ~ treat + block', rabbit).fit()
sm.stats.anova_lm(lmod,typ=3).round(3)
Out[39]:
In [40]:
sns.residplot(lmod.fittedvalues, lmod.resid)
plt.xlabel("Fitted values")
plt.ylabel("Residuals")
plt.show()
In [41]:
res = np.sort(lmod.resid)
n = len(res)
qq = stats.norm.ppf(np.arange(1,n+1)/(n+1))
fig, ax = plt.subplots()
ax.scatter(qq, res)
plt.xlabel("Normal Quantiles")
plt.ylabel("Sorted Residuals")
plt.show()
In [42]:
lmod.summary()
Out[42]:
Tukey pairwise differences:
In [43]:
mcoefs = np.append([0],lmod.params[1:6])
nmats = [chr(i) for i in range(ord('a'),ord('f')+1)]
In [44]:
p = len(mcoefs)
dp = set(itertools.combinations(range(0,p),2))
dcoef = []
namdiff = []
for cp in dp:
dcoef.append(mcoefs[cp[0]] - mcoefs[cp[1]])
namdiff.append(nmats[cp[0]] + '-' + nmats[cp[1]])
thsd = pd.DataFrame({'Difference':dcoef},index=namdiff)
thsd["lb"] = thsd.Difference - get_tukeyQcrit(p,lmod.df_resid) * lmod.bse[1]/np.sqrt(2)
thsd["ub"] = thsd.Difference + get_tukeyQcrit(p,lmod.df_resid) * lmod.bse[1]/np.sqrt(2)
thsd.round(2)
Out[44]:
Relative efficiency
In [45]:
lmodr = smf.ols('gain ~ treat', rabbit).fit()
lmodr.scale/lmod.scale
Out[45]:
In [46]:
%load_ext version_information
%version_information pandas, numpy, matplotlib, seaborn, scipy, patsy, statsmodels
Out[46]: