<jupyter_start><jupyter_text>---title: "Connecting a paired t-test with a Bayesian multilevel model (DRAFT)"toc: truetoc_sticky: truemathjax: true---I wondered how to get the same numerical answer for the t-statistic when calculating by hand, using a `scipy.stats` function, and using a Bayesian multilevel model.Dr. McElreath states:> I think about it like this: ANOVA etc are procedures applied tospecific linear models. Those models are neither Bayesian nornon-Bayesian. One can justify those models from differentperspectives: measurement or causal identification. And then how toconstruct an estimator for the piece of the model that you want is yetanother step.> Many non-Bayesian estimators can be understood as approximations ofthe posterior distribution. But lots of non-Bayesian estimators arejust quantities derived from specific models. There is no generaltheory for them. A paired t-test for example is a particular estimatorfor a specific kind of mixed model.<jupyter_code>import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as pm
import scipy.stats as stats
import seaborn as sns
import daft
from causalgraphicalmodels import CausalGraphicalModel

from scipy.optimize import curve_fit
%load_ext nb_black
%config InlineBackend.figure_format = 'retina'
%load_ext watermark
RANDOM_SEED = 8927
np.random.seed(RANDOM_SEED)
az.style.use("arviz-darkgrid")
az.rcParams["stats.hdi_prob"] = 0.89  # sets default credible interval used by arviz
sns.set_context("talk")
def standardize(x):
    x = (x - np.mean(x)) / np.std(x)
    return x<jupyter_output><empty_output><jupyter_text>A simple dataset taken from WikipediaThe data I'll use is taken from [Wikipedia's entry on the paired t-test](https://en.wikipedia.org/wiki/Student%27s_t-testDependent_t-test_for_paired_samples). It's a simple set with four subjects taking a test twice, an example of repeated measures.  Interestingly, in the Chapter 1 introduction and in Chapter 13, Dr. McElreath highlights multi-level models as a way to get "improved estimates for repeat sampling".<jupyter_code>df_tests = pd.DataFrame(
    {
        "Name": ["Mike", "Melanie", "Melissa", "Mitchell"],
        "Test_1": [35, 50, 90, 78],
        "Test_2": [67, 46, 86, 91],
    }
)
df_tests<jupyter_output><empty_output><jupyter_text>Calculate the t-statisticThe paired t-test is calculating whether the average of the difference is significantly different. We can get this quickly with a Python function.  T-statistic using scipy.stats<jupyter_code>stats.ttest_rel(df_tests["Test_1"], df_tests["Test_2"])<jupyter_output><empty_output><jupyter_text>T-statistic calculated manuallyThe formula for the t-statistic is given on the same Wikipedia page. $$t = \frac{\bar{X}_D - \mu_0}{s_D / \sqrt{n}}$$where $\bar{X}_D$ is the average of the differences between all pairs and $s_D$ is the standard deviation of the differences. We can set $mu_0$ to 0 if we want to test whether the average is significantly different. This is straightforward to calculate.<jupyter_code>diff = df_tests["Test_1"] - df_tests["Test_2"]
t_stat = diff.mean() / (diff.std() / np.sqrt(len(df_tests)))
print("Manual calculated t-statistic: ", t_stat)<jupyter_output>Manual calculated t-statistic:  -1.0784834690588145<jupyter_text>We get exactly the same answer as `scipy.stats` so that is comforting.  Calculate the t-statistic with a Bayesian approachNow the fun part. We can set this up as a Bayesian multi-level linear model. A non-Bayesian approach for connecting the paired t-test to a linear mixed model was described [here](https://vasishth.github.io/Freq_CogSci/from-the-paired-t-test-to-the-linear-mixed-model.html).First, let's do some table reformatting, so that we can use it in our model. We'll represent name and times into numerical codes and standardize the test scores.<jupyter_code>df_tests2 = pd.melt(df_tests, id_vars="Name")
df_tests2.columns = ["Name", "Time", "Score"]
df_tests2["Name_code"] = pd.Categorical(df_tests2["Name"]).codes
df_tests2["Time_code"] = pd.Categorical(df_tests2["Time"]).codes
df_tests2["Score_std"] = standardize(df_tests2["Score"])
df_tests2<jupyter_output><empty_output><jupyter_text>There are different "clusters" here. On one level, is time, whether it was the first or second test. On another level are the individual subjects.  Non multi-level model with time as an indexed, categorical variableLet's start without a multi-level model. Since the main question is whether there are differences between the first and second tests, we'll use only time as a cluster for now and ignore the subject cluster. This is an intercept only model, with time represented as an indexed categorical variable. Not sure if I set this up right. I think we'd still use a Normal here but should I use a Student t distribution instead. Also I used a flat prior since non-Bayesians wouldn't regularize . $$\text{score}_i \text{ ~ Normal}(\mu_i, \sigma) $$$$\mu_i = \alpha_{\text{TIME}}$$$$\alpha_j \text{ ∼ Normal}(0, 10)$$$$\sigma \text{ ~ Exp}(1)$$<jupyter_code>with pm.Model() as m1:

    # prior for SD of testers
    sigma = pm.Exponential("sigma", 1.0)

    # regularizing prior
    a = pm.Normal("a", 0, 10, shape=2)  # two time points

    # mu is deterministic, equivalent to alpha indexed by time
    mu = a[df_tests2["Time_code"]]

    # likelihood
    score = pm.Normal("score", mu=mu, sd=sigma, observed=df_tests2["Score_std"])

    trace_m1 = pm.sample(
        draws=1000, random_seed=19, return_inferencedata=True, progressbar=False
    )
az.summary(trace_m1)<jupyter_output><empty_output><jupyter_text>Not sure if I set this up right and even if I did, I'm not sure what to do next. Now let's calculate the contrast between the alpha parameters representing the two timepoints (??)<jupyter_code># Put the trace object in a dataframe so we can pull out the posteriors
trace_m1_df = trace_m1.to_dataframe()

# calculate contrast
post_diff = (
    trace_m1_df[("posterior", "a[0]", 0)] - trace_m1_df[("posterior", "a[1]", 1)]
)<jupyter_output><empty_output><jupyter_text>I'm just using the t-statistic equation here.<jupyter_code>post_diff.mean() / (post_diff.std() / np.sqrt(4))<jupyter_output><empty_output><jupyter_text>With a pretty flat prior (what non-Bayesians would do), I got pretty close to the answer, but not sure if I just lucky.  Multi-level modelWe can now recognize both clusters and let the model do adaptive regularization. Not sure if I set this up right. I'm unsure what the top line sigma would be for example. $$\text{score}_i \text{ ~ Normal}(\mu_i, \sigma) $$$$\mu_i = \alpha_{\text{TIME}} + \gamma_{\text{SUBJECT}} $$$$\alpha_j \text{ ∼ Normal}(\bar{\alpha}, \sigma_\alpha)$$$$\gamma_j \text{ ∼ Normal}(0, \sigma_\gamma) $$$$\bar{\alpha} \text{ ∼ Normal}(0, 10)$$$$\sigma_\alpha \text{ ~ Exp}(1)$$$$\sigma_\gamma \text{ ~ Exp}(1)$$Alpha is using an adaptive prior for time (j=1,2) and gamma is using adaptive prior for subject (j = 1..4). We'll use one global mean parameter ("bar_alpha") since each varying intercept type will be added to the same linear prediction. This global mean parameter is essentially the average across all tests and subjects. Like in the previous model, it is pretty wide to represent the non-Bayesian perspective. Each cluster also gets its own sigma.<jupyter_code>with pm.Model() as m2:
    # Top line sigma (not sure if this is right)
    sigma = pm.Exponential("sigma", 1.0)

    # prior for average person and timepoint
    a_bar = pm.Normal("a_bar", 0.0, 10)

    # prior for SD of testers and time
    sigma_a = pm.Exponential("sigma_a", 1.0)
    sigma_g = pm.Exponential("sigma_g", 1.0)

    # adaptive priors?
    a = pm.Normal("a", a_bar, sigma_a, shape=2)  # 2 time points
    g = pm.Normal("g", 0, sigma_g, shape=4)  # 4 subjects

    mu = a[df_tests2["Time_code"]] + g[df_tests2["Name_code"]]

    # likelihood
    score = pm.Normal("score", mu=mu, sd=sigma, observed=df_tests2["Score_std"])

    trace_m2 = pm.sample(
        draws=1000, random_seed=19, return_inferencedata=True, progressbar=False
    )<jupyter_output>Auto-assigning NUTS sampler...
INFO:pymc3:Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
INFO:pymc3:Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
INFO:pymc3:Multiprocess sampling (4 chains in 4 jobs)
NUTS: [g, a, sigma_g, sigma_a, a_bar, sigma]
INFO:pymc3:NUTS: [g, a, sigma_g, sigma_a, a_bar, sigma]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 14 seconds.
INFO:pymc3:Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 14 seconds.
There were 88 divergences after tuning. Increase `target_accept` or reparameterize.
ERROR:pymc3:There were 88 divergences after tuning. Increase `target_accept` or reparameterize.
There were 118 divergences after tuning. Increase `target_accept` or reparameterize.
ERROR:pymc3:There were 118 divergences after tuning. Increase `target_accept` or reparameterize.
There were 96 divergences after tuning. Increas[...]<jupyter_text>The divergences indicate I need to re-paramaterize, but let me know if my overall approach is right<jupyter_code>az.summary(trace_m2)
# Put the trace object in a dataframe so we can pull out the posteriors
trace_m2_df = trace_m2.to_dataframe()

# calculate contrast
post_diff2 = (
    trace_m2_df[("posterior", "a[0]", 0)] - trace_m2_df[("posterior", "a[1]", 1)]
)<jupyter_output><empty_output><jupyter_text>Again, I'm just using the t-statistic equation here.<jupyter_code>post_diff2.mean() / (post_diff2.std() / np.sqrt(4))<jupyter_output><empty_output><jupyter_text>I get an answer further than the previous model but in the right direction. Again, not sure if I had done something incorrect Appendix: Environment and system parameters<jupyter_code>%watermark -n -u -v -iv -w<jupyter_output>Last updated: Wed Aug 04 2021

Python implementation: CPython
Python version       : 3.8.6
IPython version      : 7.20.0

seaborn   : 0.11.1
scipy     : 1.6.0
matplotlib: 3.3.4
arviz     : 0.11.1
daft      : 0.1.0
pandas    : 1.2.1
pymc3     : 3.11.0
json      : 2.0.9
numpy     : 1.20.1

Watermark: 2.1.0
{'fraction_non_alphanumeric': 0.10390755068063853, 'fraction_numerical': 0.02467926745666438, 'mean_word_length': 4.7727529677784055, 'pattern_counts': {'":': 3, '<': 46, '<?xml version=': 0, '>': 48, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '398555', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3494, 'n_tokens_neox': 3235, 'n_words': 1164, 'path': '_posts_drafts/2021-08-04-stats_rethinking_paired_ttest.ipynb', 'repo_name': 'benslack19/benslack19.github.io', 'size': '10211'}