Stanford STATS191 in Python, Lecture 8 : Analysis of variance:(Part 1)
Stanford Stats 191¶
Introduction¶
This is a re-creation of the Stanford Stats 191 course (see https://web.stanford.edu/class/stats191/), using Python eco-system tools, instead of R. This is lecture "Analysis of variance:", Part 1
Initial Notebook Setup¶
watermark documents the current Python and package environment, black is my preferred Python formatter
%load_ext watermark
%load_ext lab_black
%matplotlib inline
import pandas as pd
import numpy as np
import seaborn as sn
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy import stats
from statsmodels.formula.api import ols
from statsmodels.formula.api import rlm
import statsmodels.api as sm
from patsy import dmatrices
from patsy.contrasts import Treatment
from patsy.contrasts import Sum
from patsy.contrasts import Diff
from statsmodels.sandbox.regression.predstd import (
wls_prediction_std,
)
import warnings
Load Datset¶
In this post, we examine a range of datasets, where the explanatory variables are all categorical.
We load the first example which records the time to recover after medical treatment, and where patients fitness levels have been recorded
data = pd.read_csv('../data/rehab.csv')
Get the pandas descriptive quantative summary
data.describe()
Visualize Dataset¶
We use pandas and seaborn to draw boxplots. It does appear that increased Fitness will reduce the recovery Time
fig, ax = plt.subplots(figsize=(8, 8))
data.boxplot(by='Fitness', ax=ax)
sn.boxplot(x='Fitness', y='Time', data=data)
Fit Model by OLS¶
We first declare the Fitness level to be a categorical (not numeric) variable, and then perform an OLS best for of a model where Time depend on Fitness. As we expect, the coefficients for Fitness levels greater than 1 are negative, predicting decreased recovery Time, the more Fitness the patient has
data['Fitness'] = data['Fitness'].astype('category')
data.dtypes
res1 = ols('Time ~ Fitness ', data=data).fit()
res1.summary()
Under the best-fit model, the prediction for Time for a given Fitness will equal the mean for all actual values of Time with theat level of Fitness.
First, we get the means for all levels of Fitness
tr_means = data.groupby('Fitness').mean()
tr_means
tr_counts = data.groupby('Fitness').count()
tr_counts
And now we see the predicted values for all levels of Fitness
res1.predict()
res1.predict({'Fitness': [1, 2, 3]})
Explicit computation of F Statistic¶
#
# get sum of squares of the residual error (predicted - actual)
sse = sum(
(data['Time'] - res1.predict())
* (data['Time'] - res1.predict())
)
mean_sse = sse / (24 - 2 - 1)
mean_sse
#
# get sum of squares of the residual around the mean
ssr = sum(
(data['Time'].mean() - res1.predict())
* (data['Time'].mean() - res1.predict())
)
mean_ssr = ssr / (3 - 1)
mean_ssr
#
# get sum of squares of the actual around the mean
sum(
(data['Time'] - data['Time'].mean())
* (data['Time'] - data['Time'].mean())
)
F = mean_ssr / mean_sse
F
Using Statsmodels ANOVA¶
sm.stats.anova_lm(res1)
Explicit computation of seeing an F value this big, with these degrees of freedom
df_tr = 2
df_res = 21
1 - stats.f.cdf(F, df_tr, df_res)
We can aso use the Survival Function, which is more accurate for small values (here, very close)
stats.f.sf(F, df_tr, df_res)
x = np.linspace(0.001, 30, 100)
y = stats.f.pdf(x, df_tr, df_res)
plt.plot(x, y, label='F pdf')
plt.axhline(0, color='k')
plt.axvline(F, color='r', label='F value')
plt.axvline(0, color='k')
plt.legend(loc='best')
Data Set 2¶
The second dataset again has two categorical variables (Weight, and Duration). The number of Days under care is converted to a log value, maybe because we are interested in the ratio of recovery time, not the actual numbers
| Variable | Description |
|---|---|
Days| Duration of hospital stay Weight| How much weight is gained? Duration| How long under treatment for kidney problems? (two levels)
Read and Explore Dataset¶
data = pd.read_fwf('../data/kidneytable.txt')
len(data)
data.head()
data['Duration'].unique()
data['Weight'].unique()
data.hist(column='Days')
Set datatype for analysis¶
We convert Days to LogDays, compute a Group variable (ie which group an observation is in, given its Weight and Duration values), and declare Weight and Duration to be categorical variables
data['Group'] = (data['Duration'] - 1) * 3 + data['Weight']
data['Duration'] = data['Duration'].astype('category')
data['Weight'] = data['Weight'].astype('category')
data['LogDays'] = np.log(data['Days'] + 1)
data.head()
data.hist(column='LogDays')
Interaction Plot¶
We can assess if there are inter-variable interactions by the interaction plot. If the lines are not parallel, then interactions may be present
_ = sm.graphics.interaction_plot(
data['Weight'],
data['Duration'],
data['LogDays'],
markers=['d', '^'],
markersize=10,
)
Build a Model¶
The model we build is different to previous models (see equation below). Here, alpha_i is the coefficient of the i_th level of the first categorical variable (Duration), beta_j is the coefficient of the j_th level of the second categorical variable, and alphabeta_ij is the coefficient of the i_th level of the first categorical variable (Duration) being present with the j_th level of the second categorical variable
The Null Hypothesis is that there are no main effects (alpha_i = 0 for all i, beta_j = 0 for all j, and alphabeta_ij = 0 for all i,j)
The default way statsmodels OLS builds the matrix of observations for categorical variables is to create new variable (one for each level of the category, except the base level). These new variables are 1 if the observation has that value of the category, and zero otherwise.
In this example, we want the sum of the coefficients to be zero for both categorical variables. We declare this in the ols call, by specifying we want the new variable to be coded so the sum is zero. I must confess I don't know if the C stands for Contrast, Coding or Categorical (I think the latter)
res2 = ols(
'LogDays ~ C(Duration, Sum) * C(Weight, Sum)', data=data
).fit()
res2.summary()
It appears from the above that the interaction effects are not significant. The effect of Duration is significant (at the 5% level). The first coded Weight variable has a significant effect, but not the other coded Weight variable.
Get Design Matrix¶
The Design Matrix is the numerical, matrix representation of our model. Note that we no longer have 0,1 in our new variables, but negative values as well
y, X = dmatrices(
'LogDays ~ C(Duration, Sum) * C(Weight, Sum)',
data=data,
return_type='dataframe',
)
X.head()
X.tail()
Get Coding Rules¶
We can get the new coding rules by calling the Sum object. For example, Weight has three level, so we have two new variables C(Weight, Sum)[S.1] and C(Weight, Sum)[S.2]. The value of each of these new variables has (as Weight ranges over it three level), is given by the columns of the Weight Coding matrix show below. Similarly, the values of the C(Duration, Sum)[S.1] variable (for the two levels of Duration) is in the Duration Coding matrix below
All this matters, because if we want to an explicit calculation of a predicted value given observed Duration and Weight, we have to calculate the value of these newly created variables. The coded_observation function below computes the vector of new variable values, to use with the coefficients returned by the ols call, to get predicted value.
w_levels = [1, 2, 3]
w_contrast = Sum().code_without_intercept(w_levels)
print('Weight Coding\n', w_contrast.matrix)
d_levels = [1, 2]
d_contrast = Sum().code_without_intercept(d_levels)
print('Duration Coding\n', d_contrast.matrix)
d1_codes = {1: 1, 2: -1}
w1_codes = {1: 1, 2: 0, 3: -1}
w2_codes = {1: 0, 2: 1, 3: -1}
def coded_observation(D, W, d1_codes, w1_codes, w2_codes):
vec = [
1, # constant
d1_codes[D], # D1
w1_codes[W], # W1
w2_codes[W], # W2
d1_codes[D] * w1_codes[W], # D1:W1
d1_codes[D] * w2_codes[W], # D1:W2
]
return vec
# end coded_observation
Compare the explicitly calculated prediction with the RegressionResults.predict() method
for D in [1, 2]:
for W in [1, 2, 3]:
dd = {'Duration': D, 'Weight': W}
pred1 = res2.predict(dd)[0]
vec1 = coded_observation(
D, W, d1_codes, w1_codes, w2_codes
)
pred2 = sum(vec1 * res2.params.values)
print(
f'Explicit Prediction vs Method Call: {pred2:10.5f} vs {pred1:10.5f}'
)
# end for
# end for
data['Num'] = list(range(len(data)))
fig, ax = plt.subplots(figsize=(10, 10))
sn.scatterplot(
'Num',
'LogDays',
data=data,
hue='Weight',
alpha=1,
style='Duration',
cmap='rainbow',
palette=['r', 'g', 'b'],
ax=ax,
s=200,
markers=['d', '^'],
)
sn.scatterplot(
'Num',
res2.predict(),
data=data,
hue='Weight',
alpha=1,
style='Duration',
cmap='rainbow',
palette=['r', 'g', 'b'],
ax=ax,
s=200,
markers=['+', 'x'],
legend=False,
)
Plot Residuals against observation number
data['Num'] = list(range(len(data)))
fig, ax = plt.subplots(figsize=(10, 10))
sn.scatterplot(
'Num',
res2.resid,
data=data,
hue='Weight',
alpha=1,
style='Duration',
cmap='rainbow',
palette=['r', 'g', 'b'],
ax=ax,
s=200,
markers=['d', '^'],
)
Plot Residuals Grouped
data['Num'] = list(range(len(data)))
fig, ax = plt.subplots(figsize=(10, 10))
sn.scatterplot(
'Group',
res2.resid,
data=data,
hue='Weight',
alpha=1,
style='Duration',
cmap='rainbow',
palette=['r', 'g', 'b'],
ax=ax,
s=200,
markers=['d', '^'],
)
Predictions¶
Supposed we are interested if the group with Duration level 1, Weight level 3 is significantly different from Duration level 2, Weight level 2.
We compute the predicted value of LogDays for each group
g1 = {'Duration': 1, 'Weight': 3}
g2 = {'Duration': 2, 'Weight': 2}
g1_pred = res2.predict(g1)[0]
g2_pred = res2.predict(g2)[0]
print(g1_pred, g2_pred)
Get the estimated value of sigma, the standard deviation of the noise
sigma_hat = np.sqrt(res2.mse_resid)
sigma_hat
Because we are looking at the difference of means of samples of size n=10, the variance of the estimated difference is as below
We the upper and lower bounds of our estimate, by using the T distribution. 54 is our degrees of freedom of our residuals (from the RegressionResults.summary() output), and we want the 5% confidence interval
n_obs = 10
lower = stats.t.ppf(0.025, 54) * np.sqrt(
2 * sigma_hat ** 2 / n_obs
)
upper = stats.t.ppf(0.975, 54) * np.sqrt(
2 * sigma_hat ** 2 / n_obs
)
lower, upper
Get the difference between the two groups
g_diff = g1_pred - g2_pred
g_diff
Show the upper and lower bounds. Based upon these values, we would say that at the 5% level, these two groups have a different LogDays value
g_diff + lower, g_diff + upper
The code below shows how to get Confidence Intervals for a Prediction, using the PredictionResults object
gp = res2.get_prediction(
{'Duration': 1, 'Weight': 1}
) # PredictionResults Object
gp.summary_frame()
ANOVA¶
The effect of the interactions can be assessed by a anova_lm call. We see that the two categorical variables explain significantly (at the 5% level) more variation that would be predicted by random chance, but not so the interaction between the two variables.
with warnings.catch_warnings():
warnings.simplefilter("ignore")
anv = sm.stats.anova_lm(res2)
# end with
anv
We can also compare a simple model (with no interactions) with out interaction model. The difference between the two models would be seen about 50% of the time by random chance
with warnings.catch_warnings():
warnings.simplefilter("ignore")
anv = sm.stats.anova_lm(
ols('LogDays ~ Duration + Weight', data=data).fit(),
res2,
)
# end with
anv
Environment¶
%watermark -h -iv
%watermark