Chapter 16 - Metric-Predicted Variable on One or Two Groups¶
In [1]:
import arviz as az
import numpy as np
import pandas as pd
import pymc as pm
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.stats import norm, t
from IPython.display import Image
plt.style.use('seaborn-white')
color = '#87ceeb'
In [2]:
%load_ext watermark
%watermark -p matplotlib,numpy,pandas,pymc,seaborn,scipy
matplotlib: 3.5.1 numpy : 1.23.1 pandas : 1.4.3 pymc : 5.0.0 seaborn : 0.12.2 scipy : 1.8.1
Data¶
In [3]:
df = pd.read_csv('data/TwoGroupIQ.csv', dtype={'Group':'category'})
df.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 120 entries, 0 to 119 Data columns (total 2 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 Score 120 non-null int64 1 Group 120 non-null category dtypes: category(1), int64(1) memory usage: 1.3 KB
In [4]:
df.head()
Out[4]:
| Score | Group | |
|---|---|---|
| 0 | 102 | Smart Drug |
| 1 | 107 | Smart Drug |
| 2 | 92 | Smart Drug |
| 3 | 101 | Smart Drug |
| 4 | 110 | Smart Drug |
In [5]:
# Mean and standard deviation
df.groupby('Group').agg(['mean', 'std'])
Out[5]:
| Score | ||
|---|---|---|
| mean | std | |
| Group | ||
| Placebo | 100.035088 | 17.894497 |
| Smart Drug | 107.841270 | 25.445201 |
In [6]:
fig, axes = plt.subplots(nrows=1, ncols=2)
for i, group in enumerate(["Placebo", "Smart Drug"]):
sns.histplot(x="Score", color="#87ceeb", data=df[df["Group"] == group], ax=axes[i])
axes[i].title.set_text(group)
plt.tight_layout()
In [7]:
# we are only interested in the scores of group "Smart Drug"
y = df[df["Group"] == "Smart Drug"]["Score"]
16.1 - Estimating the mean and standard deviation of a normal distribution¶
Model (Kruschke, 2015)¶
In [8]:
Image('images/fig16_2.png', width=300)
Out[8]:
In [9]:
with pm.Model() as model:
mu = pm.Normal("mu", y.mean(), sigma=y.std())
sigma = pm.Uniform("sigma", y.std() / 1000, y.std() * 1000)
# PyMC's Normal likelihood can take either precision or standard deviation as an argument.
likelihood = pm.Normal("likelihood", mu, sigma=sigma, observed=y)
In [10]:
pm.model_to_graphviz(model)
Out[10]:
In [11]:
with model:
idata = pm.sample(5000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (4 chains in 4 jobs) NUTS: [mu, sigma]
Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 2 seconds.
In [12]:
az.plot_trace(idata);
Figure 16.3¶
In [13]:
fig, [(ax1, ax2), (ax3, ax4)] = plt.subplots(2, 2, figsize=(10, 6))
font_d = {"size": 16}
# Upper left
pm.plot_posterior(
idata.posterior["mu"], point_estimate="mode", ref_val=100, ax=ax1, color=color
)
ax1.set_xlabel("$\mu$", fontdict=font_d)
ax1.set_title("Mean", fontdict=font_d)
# Upper right
tr_len = idata.posterior.sizes["draw"]
# Plot only 20 posterior prediction curves.
n_curves = 20
# Create an index of length 20 with which we step through the trace.
stepIdxVec = np.arange(0, tr_len, tr_len // n_curves)
x_range = np.arange(y.min(), y.max())
x = np.tile(x_range.reshape(-1, 1), (1, n_curves))
# x = np.linspace(y.min(), y.max(), n_curves)
ax2.hist(y, bins=25, density=True, color="steelblue")
ax2.plot(
x,
norm.pdf(
x,
idata.posterior["mu"].sel(chain=0, draw=stepIdxVec),
idata.posterior["sigma"].sel(chain=0, draw=stepIdxVec),
),
c=color,
)
ax2.set_xlabel("y", fontdict=font_d)
ax2.set_title("Data w. Post. Pred.\nN=63")
[ax2.spines[spine].set_visible(False) for spine in ["left", "right", "top"]]
ax2.yaxis.set_visible(False)
# Lower left
pm.plot_posterior(
idata.posterior["sigma"], point_estimate="mode", ref_val=15, ax=ax3, color=color
)
ax3.set_xlabel("$\sigma$", fontdict=font_d)
ax3.set_title("Std. Dev.", fontdict=font_d)
# Lower right
pm.plot_posterior(
(idata.posterior["mu"] - 100) / idata.posterior["sigma"],
point_estimate="mode",
ref_val=0,
ax=ax4,
color=color,
)
ax4.set_title("Effect Size", fontdict=font_d)
ax4.set_xlabel("$(\mu - 100)/\sigma$", fontdict=font_d)
plt.tight_layout()
In [14]:
az.summary(idata)
Out[14]:
| mean | sd | hdi_3% | hdi_97% | mcse_mean | mcse_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|
| mu | 107.851 | 3.288 | 101.713 | 114.086 | 0.025 | 0.017 | 18022.0 | 13482.0 | 1.0 |
| sigma | 26.015 | 2.401 | 21.637 | 30.531 | 0.019 | 0.013 | 16746.0 | 13980.0 | 1.0 |
Mean of mu is 108ish, mean of sigma is 26ish. What is the mean effect size (e.g., (mean-100)/sigma)? About 0.3 or so? Did you forget that mu and sigma are simply 2 dimensions of a multivariate posterior? To determine our posterior for the effect size, we need to calculate the effect size for each mu-sigma pair that we have in our trace (our logbook). Then we can take the mean of these effect sizes.
In [15]:
((idata.posterior["mu"] - 100) / idata.posterior["sigma"]).mean()
Out[15]:
<xarray.DataArray ()> array(0.30422986)
Fortunately (for you), your initial guess wasn't that far off because mu and sigma are distributed more or less independently. But next time you may not be so lucky!
In [16]:
# test mu parameter against 100 (definitionally 'average' within the population)
(idata.posterior["mu"] > 100).mean()
Out[16]:
<xarray.DataArray 'mu' ()> array(0.99045)
In [17]:
# test mu parameter against 15 (definitionally the SD of the population)
(idata.posterior["sigma"] > 15).mean()
Out[17]:
<xarray.DataArray 'sigma' ()> array(1.)
In [18]:
# test the standardized effect size against zero
# i.e., what is the probability that the mean is > 0 SDs above 100?
(((idata.posterior['mu']-100)/idata.posterior['sigma']) > 0).mean()
Out[18]:
<xarray.DataArray ()> array(0.99045)
In [19]:
# test effect size against 0.1 (e.g., top of Kruschke's ROPE)
# i.e., what is the probability that the mean is > 0.1 SDs above 100?
(((idata.posterior['mu']-100)/idata.posterior['sigma'])>0.1).mean()
Out[19]:
<xarray.DataArray ()> array(0.94305)
16.2 - Outliers and robust estimation: the t distribution¶
Model¶
In [20]:
with pm.Model() as model2:
mu = pm.Normal("mu", y.mean(), sigma=y.std())
sigma = pm.Uniform("sigma", y.std() / 1000, y.std() * 1000)
nu_minus1 = pm.Exponential("nu_minus1", 1 / 29)
nu = pm.Deterministic("nu", nu_minus1 + 1)
likelihood = pm.StudentT("likelihood", nu, mu, sigma=sigma, observed=y)
In [21]:
pm.model_to_graphviz(model2)
Out[21]:
In [22]:
with model2:
idata2 = pm.sample(5000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (4 chains in 4 jobs) NUTS: [mu, sigma, nu_minus1]
Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 4 seconds.
In [23]:
az.plot_trace(idata2)
plt.tight_layout();
Figure 16.9¶
In [24]:
fig, [(ax1, ax2), (ax3, ax4), (ax5, ax6)] = plt.subplots(3, 2, figsize=(10, 8))
# Upper left
pm.plot_posterior(
idata2.posterior["mu"], point_estimate="mode", ref_val=100, ax=ax1, color=color
)
ax1.set_xlabel("$\mu$", fontdict=font_d)
ax1.set_title("Mean", fontdict=font_d)
# Upper right
tr_len = idata2.posterior.sizes["draw"]
n_curves = 20
stepIdxVec = np.arange(0, tr_len, tr_len // n_curves)
x_range = np.arange(y.min(), y.max())
x = np.tile(x_range.reshape(-1, 1), (1, 20))
ax2.hist(y, bins=25, density=True, color="steelblue")
ax2.plot(
x,
norm.pdf(
x,
idata2.posterior["mu"].sel(chain=0, draw=stepIdxVec),
idata2.posterior["sigma"].sel(chain=0, draw=stepIdxVec),
),
c="#87ceeb",
)
ax2.set_xlabel("y", fontdict=font_d)
ax2.set_title("Data w. Post. Pred.")
[ax2.spines[spine].set_visible(False) for spine in ["left", "right", "top"]]
ax2.yaxis.set_visible(False)
# Middle left
pm.plot_posterior(
idata2.posterior["sigma"], point_estimate="mode", ref_val=15, ax=ax3, color=color
)
ax3.set_xlabel("$\sigma$", fontdict=font_d)
ax3.set_title("Scale", fontdict=font_d)
# Middle right
pm.plot_posterior(
(idata2.posterior["mu"] - 100) / idata2.posterior["sigma"],
point_estimate="mode",
ref_val=0,
ax=ax4,
color=color,
)
ax4.set_title("Effect Size", fontdict=font_d)
ax4.set_xlabel("$(\mu - 100)/\sigma$", fontdict=font_d)
# Lower left
pm.plot_posterior(
np.log10(idata2.posterior["nu"]), point_estimate="mode", ax=ax5, color=color
)
ax5.set_title("Normality", fontdict=font_d)
ax5.set_xlabel(r"log10($\nu$)", fontdict=font_d)
plt.tight_layout()
ax6.set_visible(False)
Note that sigma does not control the standard deviation of the t-distribution. This means that comparing sigma against a value of 15 is not meaningful. More importantly, this means that (mu - 100)/sigma is not a conventional standardized effect size.
In [25]:
az.summary(idata2)
Out[25]:
| mean | sd | hdi_3% | hdi_97% | mcse_mean | mcse_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|
| mu | 107.175 | 2.841 | 101.677 | 112.391 | 0.027 | 0.019 | 11349.0 | 11429.0 | 1.0 |
| sigma | 19.890 | 3.467 | 13.600 | 26.493 | 0.035 | 0.025 | 9411.0 | 9254.0 | 1.0 |
| nu_minus1 | 8.857 | 13.173 | 0.075 | 27.319 | 0.133 | 0.094 | 9485.0 | 10887.0 | 1.0 |
| nu | 9.857 | 13.173 | 1.075 | 28.319 | 0.133 | 0.094 | 9485.0 | 10887.0 | 1.0 |
In [26]:
# test mu parameter against 100 (definitionally 'average' within the population)
(idata2.posterior["mu"] > 100).mean()
Out[26]:
<xarray.DataArray 'mu' ()> array(0.99355)
In [27]:
# lower alpha to see where the vast majority of samples are
alpha = 0.5
#alpha = 0.005
az.plot_pair(
idata2,
var_names=["mu", "sigma", "nu"],
scatter_kwargs={"alpha": alpha},
figsize=(7, 7),
)
Out[27]:
array([[<AxesSubplot:ylabel='sigma'>, <AxesSubplot:>],
[<AxesSubplot:xlabel='mu', ylabel='nu'>,
<AxesSubplot:xlabel='sigma'>]], dtype=object)
16.2 - Two Groups¶
Model (Kruschke, 2015)¶
In [28]:
Image('images/fig16_11.png', width=400)
Out[28]:
In [29]:
# compute indicies for study groups
grp_idx, grp_codes = pd.factorize(df["Group"])
n_grps = df["Group"].nunique()
with pm.Model(coords = {"group":["Smart Drug","Placebo"]}) as model3:
mu = pm.Normal("mu", df["Score"].mean(), sigma=df["Score"].std(), dims="group")
sigma = pm.Uniform(
"sigma", df["Score"].std() / 1000, df["Score"].std() * 1000, dims="group"
)
nu_minus1 = pm.Exponential("nu_minus1", 1 / 29)
nu = pm.Deterministic("nu", nu_minus1 + 1)
likelihood = pm.StudentT(
"likelihood", nu, mu[grp_idx], sigma=sigma[grp_idx], observed=df["Score"]
)
In [30]:
pm.model_to_graphviz(model3)
Out[30]:
In [31]:
with model3:
idata3 = pm.sample(5000)
Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Initializing NUTS using jitter+adapt_diag... Multiprocess sampling (4 chains in 4 jobs) NUTS: [mu, sigma, nu_minus1]
Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 4 seconds.
In [32]:
az.plot_trace(idata3)
plt.tight_layout();
Figure 16.12¶
In [33]:
tr3_mu0 = idata3.posterior["mu"].sel(group="Placebo")
tr3_mu1 = idata3.posterior["mu"].sel(group="Smart Drug")
tr3_sigma0 = idata3.posterior["sigma"].sel(group="Placebo")
tr3_sigma1 = idata3.posterior["sigma"].sel(group="Smart Drug")
tr3_nu = np.log10(idata3.posterior["nu"])
fig, axes = plt.subplots(5, 2, figsize=(12, 12))
# Left column figs
l_trace_vars = (tr3_mu0, tr3_mu1, tr3_sigma0, tr3_sigma1, tr3_nu)
l_axes_idx = np.arange(5)
l_xlabels = ("$\mu_0$", "$\mu_1$", "$\sigma_0$", "$\sigma_1$", r"log10($\nu$)")
l_titles = (
"Placebo Mean",
"Smart Drug Mean",
"Placebo Scale",
"Smart Drug Scale",
"Normality",
)
for var, ax_i, xlabel, title in zip(l_trace_vars, l_axes_idx, l_xlabels, l_titles):
pm.plot_posterior(var, point_estimate="mode", ax=axes[ax_i, 0], color=color)
axes[ax_i, 0].set_xlabel(xlabel, font_d)
axes[ax_i, 0].set_title(title, font_d)
# Right column figs
tr_len = idata3.posterior.sizes["draw"]
n_curves = 20
stepIdxVec = np.arange(0, tr_len, tr_len // n_curves)
x_range = np.arange(df["Score"].min(), df["Score"].max())
x = np.tile(x_range.reshape(-1, 1), (1, n_curves))
# 1
axes[0, 1].hist(
df[df["Group"] == "Placebo"]["Score"], bins=25, density=True, color="steelblue"
)
axes[0, 1].plot(
x,
t.pdf(
x,
loc=tr3_mu0.sel(chain=0, draw=stepIdxVec),
scale=tr3_sigma0.sel(chain=0, draw=stepIdxVec),
df=tr3_nu.sel(chain=0, draw=stepIdxVec),
),
c="#87ceeb",
)
axes[0, 1].set_xlabel("y", font_d)
[axes[0, 1].spines[spine].set_visible(False) for spine in ["left", "right", "top"]]
axes[0, 1].yaxis.set_visible(False)
axes[0, 1].set_title("Data for Placebo w. Post. Pred.", font_d)
# 2
axes[1, 1].hist(
df[df["Group"] == "Smart Drug"]["Score"], bins=25, density=True, color="steelblue"
)
axes[1, 1].plot(
x,
t.pdf(
x,
loc=tr3_mu1.sel(chain=0, draw=stepIdxVec),
scale=tr3_sigma1.sel(chain=0, draw=stepIdxVec),
df=tr3_nu.sel(chain=0, draw=stepIdxVec),
),
c="#87ceeb",
)
axes[1, 1].set_xlabel("y", font_d)
[axes[1, 1].spines[spine].set_visible(False) for spine in ["left", "right", "top"]]
axes[1, 1].yaxis.set_visible(False)
axes[1, 1].set_title("Data for Smart Drug w. Post. Pred.", font_d)
# 3-5
r_vars = (
tr3_mu1 - tr3_mu0,
tr3_sigma1 - tr3_sigma0,
(tr3_mu1 - tr3_mu0) / np.sqrt((tr3_sigma0**2 + tr3_sigma1**2) / 2),
)
r_axes_idx = np.arange(start=2, stop=5)
r_xlabels = (
"$\mu_1 - \mu_0$",
"$\sigma_1 - \sigma_0$",
r"$\frac{(\mu_1-\mu_0)}{\sqrt{(\sigma_0^2+\sigma_1^2)/2}}$",
)
r_titles = ("Difference of Means", "Difference of Scales", "Effect Size")
for var, ax_i, xlabel, title in zip(r_vars, r_axes_idx, r_xlabels, r_titles):
pm.plot_posterior(
var, point_estimate="mode", ref_val=0, ax=axes[ax_i, 1], color=color
)
axes[ax_i, 1].set_xlabel(xlabel, font_d)
axes[ax_i, 1].set_title(title, font_d)
plt.tight_layout()
Again, the lower right panel is not our belief about conventional measures of standardized effect size. We are using a t-distribution, which means that sigma is not the SD.
In [34]:
# mean of the "Smart Drug" group's location parameter
idata3.posterior["mu"].mean(dim=["draw", "chain"]).sel(group="Smart Drug")
Out[34]:
<xarray.DataArray 'mu' ()>
array(107.09970567)
Coordinates:
group <U10 'Smart Drug'In [35]:
# mean of the "Placebo" group's location parameter
idata3.posterior["mu"].mean(dim=["draw", "chain"]).sel(group="Placebo")
Out[35]:
<xarray.DataArray 'mu' ()>
array(99.27939086)
Coordinates:
group <U10 'Placebo'In [36]:
# belief about the difference in location parameters (i.e., the unstandardized effect size)
(
idata3.posterior["mu"].sel(group="Smart Drug")
- idata3.posterior["mu"].sel(group="Placebo")
).mean()
Out[36]:
<xarray.DataArray 'mu' ()> array(7.82031481)
In [37]:
# probability smart drug group location > placebo group location
(
idata3.posterior["mu"].sel(group="Smart Drug")
> idata3.posterior["mu"].sel(group="Placebo")
).mean()
Out[37]:
<xarray.DataArray 'mu' ()> array(0.99255)
