# T-test#

[1]:

import arviz as az
import bambi as bmb
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

[2]:

az.style.use('arviz-darkgrid')
np.random.seed(1234)


In this notebook we demo how to perform a two-sample Bayesian t-test to compare the mean value of two Gaussian populations using Bambi.

## Generate data#

We generate 160 values from a Gaussian with $$\mu=6$$ and $$\sigma=2.5$$ and another 120 values from a Gaussian’ with $$\mu=8$$ and $$\sigma=2$$

[3]:

np.random.seed(2935)
a = np.random.normal(6, 2.5, 160)
b = np.random.normal(8, 2, 120)
df = pd.DataFrame(
{
'Group': ['a'] * 160 + ['b'] * 120,
'Val':np.hstack([a, b])
}
)
az.plot_violin({'a':a, 'b':b});


### Model#

When we carry out a two sample t-test we’re implicitly using the following model:

$\mu_i = \beta_0 + \beta_1 (i) + \epsilon_i$

where $$i = 0$$ represents the population 1, $$i = 1$$ the popoulation 2 and $$\epsilon_i$$ is a random error with mean 0. If we replace the indicator variables for the two groups we have

$\mu_0 = \beta_0 + \epsilon_i$

and

$\mu_1 = \beta_0 + \beta_1 + \epsilon_i$

Thus, we can see that testing whether the mean of the two populations are equal is equivalent to testing whether $$\beta_1$$ is 0.

### Analysis#

We start by instantiating our model and specifying the model previously described.

[4]:

model = bmb.Model(df)
results = model.fit('Val ~ Group')

Auto-assigning NUTS sampler...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [Val_sigma, Group, Intercept]

100.00% [4000/4000 00:01<00:00 Sampling 2 chains, 0 divergences]
Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 3 seconds.


We’ve only specified the formula for the model and Bambi automatically selected priors distributions and values for their parameters. We can inspect both the setup and the priors as following:

[5]:

model.backend.model

[5]:

$$\begin{array}{rcl} \text{Intercept} &\sim & \text{Normal}(\mathit{mu}=6.888841520355615,~\mathit{sigma}=2.77556648438955)\\\text{Group} &\sim & \text{Normal}(\mathit{mu}=0.0,~\mathit{sigma}=3.71960189063118)\\\text{Val_sigma} &\sim & \text{HalfStudentT}(\mathit{nu}=4.0,~\mathit{sigma}=2.272128158801779)\\\text{Val} &\sim & \text{Normal}(\mathit{mu}=f(f(f(),~f(f(array,~f(\text{Intercept})))),~f(f(array,~f(\text{Group})))),~\mathit{sigma}=f(\text{Val_sigma})) \end{array}$$
[6]:

model.plot_priors();


To inspect our posterior and the sampling process we can call az.plot_trace(). The option kind='rank_vlines' gives us a variant of the rank plot that uses lines and dots and helps us to inspect the stationarity of the chains. Since there is no clear pattern or serious deviations from the horizontal lines, we can conclude the chains are stationary.

[7]:

az.plot_trace(results, kind='rank_vlines');

[8]:

az.summary(results)

[8]:

mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_mean ess_sd ess_bulk ess_tail r_hat
Intercept 6.171 0.171 5.856 6.486 0.005 0.004 1115.0 1109.0 1120.0 1222.0 1.0
Group 1.678 0.259 1.201 2.173 0.008 0.005 1176.0 1158.0 1176.0 1262.0 1.0
Val_sigma 2.121 0.091 1.954 2.292 0.002 0.002 1538.0 1538.0 1529.0 1052.0 1.0

In the summary table we can see the 94% highest density interval for $$\beta_1$$ ranges from 1.201 to 2.173. Thus, according to the data and the model used, we conclude the difference between the two population means is somewhere between 1.2 and 2.2 and hence we support the hypotehsis that $$\beta_1 \ne 0$$.

Similar conclusions can be made with the density estimate for the posterior distribution of $$\beta_1$$. As seen in the table, most of the probability for the difference in the mean roughly ranges from 1.2 to 2.2.

[9]:

# Grab just the posterior of the term of interest (group)
group_posterior = results.posterior['Group'].values
ax = az.plot_kde(group_posterior)
ax.axvline(0, 0, 3, linestyle='--', color='k');


Another way to arrive to a similar conclusion is by calculating the probability that the parameter $$\beta_1 > 0$$. This probability, practically equal to 1, tells us that the mean of the two populations are different.

[10]:

# Probabiliy that posterior is > 0
(group_posterior > 0).mean()

[10]:

1.0

[11]:

%load_ext watermark
%watermark -n -u -v -iv -w

numpy  1.19.1
arviz  0.10.0
bambi  0.2.0
pandas 1.1.1
last updated: Fri Dec 11 2020

CPython 3.8.5
IPython 7.18.1
watermark 2.0.2