<jupyter_start><jupyter_text>Multiple Factor ANOVAThe course document uses the [Warpbreaks](http://vincentarelbundock.github.io/Rdatasets/doc/datasets/warpbreaks.html) dataset for this example.Explanatory variables are `wool` with two values `A` and `B`, and `tension`, with 3 values `Low`, `Medium`, and `High`. Response varibale is the number of breaks.<jupyter_code>import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as pm
import scipy.stats as stats
import statsmodels.api as sm

from scipy.special import expit

%matplotlib inline
import warnings
warnings.filterwarnings("ignore")<jupyter_output><empty_output><jupyter_text>Load Data<jupyter_code>try:
    warpbreaks_df = pd.read_csv("warpbreaks.csv", index=False)
except:
    warpbreaks_df = pd.read_csv("http://vincentarelbundock.github.io/Rdatasets/csv/datasets/warpbreaks.csv")
    warpbreaks_df.to_csv("warpbreaks.csv")
    
warpbreaks_df.head()<jupyter_output><empty_output><jupyter_text>Visualize Data* Tension is categorical with 3 levels Low (L), Medium (M), and High (H).* Wool is categorical with 2 levels A and B.* Variance between groups grouped by (wool + tension) is less if we consider log(breaks).<jupyter_code>pd.crosstab(warpbreaks_df["wool"], warpbreaks_df["tension"])
data, labels = [], []
for tension in ["L", "M", "H"]:
    for wool in ["A", "B"]:
        data.append(warpbreaks_df[(warpbreaks_df["tension"] == tension) & 
                                  (warpbreaks_df["wool"] == wool)]["breaks"].values)
        labels.append("".join([wool, tension]))

plt.boxplot(data, labels=labels)
plt.ylabel("breaks")
_ = plt.show()
log_data = []
for d in data:
    logd = np.log(np.array(d))
    log_data.append(logd)
    
plt.boxplot(log_data, labels=labels)
plt.ylabel("log(breaks)")
_ = plt.show()<jupyter_output><empty_output><jupyter_text>One Way ModelUsing `tension` as the explanatory variable only.<jupyter_code>warpbreaks1_df = pd.concat((warpbreaks_df,
                            pd.get_dummies(warpbreaks_df["tension"], prefix="tension", drop_first=True)),
                           axis=1)
warpbreaks1_df.head()
Xdata = warpbreaks1_df.values[:, 4:]
Xone = np.ones(Xdata.shape[0]).reshape(-1, 1)
X = np.hstack((Xone, Xdata)).astype(np.float32)

y = np.log(warpbreaks1_df["breaks"].values)

Xdata.shape, Xone.shape, X.shape, y.shape
model_1 = pm.Model()
with model_1:
    coeffs = pm.Normal("coeffs", mu=0, sigma=1e6, shape=3)
    prec = pm.InverseGamma("prec", alpha=2.5, beta=0.4)

    sd = np.sqrt(prec)
    preds = (coeffs[0] * X[:, 0] +
             coeffs[1] * X[:, 1] + 
             coeffs[2] * X[:, 2])
    
    y_obs = pm.Normal("y_obs", mu=preds, sigma=sd, observed=y)
    
    trace_1 = pm.sample(5000, tune=1000)
_ = pm.traceplot(trace_1, combined=True, compact=False)
pm.summary(trace_1)<jupyter_output><empty_output><jupyter_text>Two-way Additive ModelWe will first fit the additive model which treats the two factors `tension` and `wool` separately with no interaction. There are 6 possible values for the 2 categorical variables taken together.**NOTE:** Model shown here is slightly different from the one in the course, it uses both levels of `wool` (A and B). This is because the model with 1 level of `wool` was not converging properly, possibly because it tries to describe the base state for two variables `tension` and `wool` in the intercept term. In this model, the intercept term contains only the base state for `tension` (H).<jupyter_code>warpbreaks2_df = pd.concat((warpbreaks_df,
                            pd.get_dummies(warpbreaks_df["tension"], prefix="tension", drop_first=True),
                            pd.get_dummies(warpbreaks_df["wool"], prefix="wool")),
                           axis=1)
warpbreaks2_df.head()
Xdata = warpbreaks2_df.values[:, -4:]
Xones = np.ones(Xdata.shape[0]).reshape(-1, 1)
X = np.hstack((Xones, Xdata)).astype(np.float32)

y = np.log(warpbreaks2_df["breaks"].values)

Xdata.shape, Xones.shape, X.shape, y.shape
model_2 = pm.Model()
with model_2:
    coeffs = pm.Normal("coeffs", mu=0, sigma=1e3, shape=5)
    prec = pm.InverseGamma("prec", alpha=2.5, beta=0.4)
    
    sd = np.sqrt(prec)
    pred = (coeffs[0] * X[:, 0] +
            coeffs[1] * X[:, 1] + 
            coeffs[0] * X[:, 2] +
            coeffs[1] * X[:, 3] + 
            coeffs[2] * X[:, 4])
    
    y_obs = pm.Normal("y_obs", mu=pred, sigma=sd, observed=y)

    trace_2 = pm.sample(5000, tune=1000)
_ = pm.traceplot(trace_2, combined=True, compact=False)
pm.summary(trace_2)<jupyter_output><empty_output><jupyter_text>Cell-Means ModelIn the additive model, the effects of `wool_A (coeffs[3])` and `wool_B (coeffs[4])` are both negative. with `wool_A` being more negative than `wool_B`. This does not seem to be borne out by the box plot, since for `tension_M`, `wool_A` seems to have lower `log(breaks)` than `wool_B`, but it is the opposite for the other `tension` classes. This leads us to suspect that there might be some interaction effects. The cell means model tries to account for interaction features, by giving each of the 6 `wool*tension` subgroups its own prior.<jupyter_code>warpbreaks3_df = pd.DataFrame.copy(warpbreaks_df, deep=True)
warpbreaks3_df.head()
for wool in ["A", "B"]:
    for tension in ["L", "M", "H"]:
        colname = "wt_{:s}{:s}".format(wool, tension)
        warpbreaks3_df[colname] = ((warpbreaks3_df["wool"] == wool) & 
                                   (warpbreaks3_df["tension"] == tension)).astype(int)
warpbreaks3_df.head()
X = warpbreaks3_df.values[:, -6:].astype(np.float32)
y = np.log(warpbreaks3_df["breaks"].values)

X.shape, y.shape
model_3 = pm.Model()
with model_3:
    prec = pm.InverseGamma("prec", alpha=1.5, beta=0.67)
    coeffs = pm.Normal("coeffs", mu=0, sigma=1e3, shape=6)

    sd = np.sqrt(prec)
    preds = (coeffs[0] * X[:, 0] +
             coeffs[1] * X[:, 1] + 
             coeffs[2] * X[:, 2] + 
             coeffs[3] * X[:, 3] + 
             coeffs[4] * X[:, 4] + 
             coeffs[5] * X[:, 5])
    
    y_obs = pm.Normal("y_obs", mu=preds, sigma=sd, observed=y)

    trace_3 = pm.sample(5000, tune=1000)
_ = pm.traceplot(trace_3, combined=True, compact=False)
pm.summary(trace_3)<jupyter_output><empty_output><jupyter_text>Find best treatment groupWe look at the trace of the posteriors, and for each generated sample, we compute the one with lowest value (the lowest `log(break)`. Counting these minimum values will give us the "best treatment group".Here out best treatment group looks like `wool_B` and `tension_H`, the same as that found in the course document.<jupyter_code>coeffs = trace_3.get_values("coeffs", combine=True)
unique, counts = np.unique(np.argmin(coeffs, axis=1), return_counts=True)
[(u, c) for u, c in zip(unique, counts)]<jupyter_output><empty_output><jupyter_text>Compare ModelsFindings:* In terms of WAIC, `model_1` appears to be the best model, but only marginally compared to `model_2`. The course document finds that `model_2` is the best in terms of DIC.<jupyter_code>compare_dict = {
    "model_1": trace_1,
    "model_2": trace_2,
    "model_3": trace_3
}
comp = pm.compare(compare_dict)
comp
_ = pm.compareplot(comp)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13821815154038303, 'fraction_numerical': 0.0202608936996947, 'mean_word_length': 3.9772099447513813, 'pattern_counts': {'":': 3, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '577446', 'license': 'Apache-2.0', 'n_tokens_mistral': 2667, 'n_tokens_neox': 2597, 'n_words': 656, 'path': 'techniques-and-models/w03-99-multifactor-anova.ipynb', 'repo_name': 'sujitpal/bayesian-stats-examples', 'size': '7206'}