Bayesian t-tests¶

In this quick tutorial we illustrate how to use posterior samples to perform Bayesian hypothesis testing.

In [ ]:

Copied!





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

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

import hssm

Example 1: Separate models¶

Simulate Data¶

We will simulate a simple dataset that contains two conditions.

In [ ]:

Copied!





# Condition 1
condition_1 = hssm.simulate_data(
    model="ddm", theta=dict(v=0.5, a=1.5, z=0.5, t=0.1), size=500
)

# Condition 2
condition_2 = hssm.simulate_data(
    model="ddm", theta=dict(v=1.0, a=1.5, z=0.5, t=0.1), size=500
)
# Condition 1
condition_1 = hssm.simulate_data(
    model="ddm", theta=dict(v=0.5, a=1.5, z=0.5, t=0.1), size=500
)

# Condition 2
condition_2 = hssm.simulate_data(
    model="ddm", theta=dict(v=1.0, a=1.5, z=0.5, t=0.1), size=500
)

Specify Models¶

We will fit two separate models to the data.

In [ ]:

Copied!

# Model 1
m1 = hssm.HSSM(model="ddm", data=condition_1)

m1.sample(sampler="mcmc", tune=500, draws=500)

# Model 2
m2 = hssm.HSSM(model="ddm", data=condition_2)

m2.sample(sampler="mcmc", tune=500, draws=500)
# Model 1
m1 = hssm.HSSM(model="ddm", data=condition_1)

m1.sample(sampler="mcmc", tune=500, draws=500)

# Model 2
m2 = hssm.HSSM(model="ddm", data=condition_2)

m2.sample(sampler="mcmc", tune=500, draws=500)

Model initialized successfully.
Using default initvals.

Initializing NUTS using adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [z, t, a, v]

Output()

Sampling 4 chains for 500 tune and 500 draw iterations (2_000 + 2_000 draws total) took 6 seconds.
100%|██████████| 2000/2000 [00:00<00:00, 5329.32it/s]

Model initialized successfully.
Using default initvals.

Initializing NUTS using adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [z, t, a, v]

Output()

Sampling 4 chains for 500 tune and 500 draw iterations (2_000 + 2_000 draws total) took 8 seconds.
100%|██████████| 2000/2000 [00:00<00:00, 5427.84it/s]

Out[ ]:

arviz.InferenceData

posterior
log_likelihood
sample_stats
observed_data

Bayesian t-test¶

Now, let's ask the (informed) question:

Is the v parameter higher in condition 1 than condition 2?

The beauty of the Bayesian approach is that this questions boils down to simple counting.

Reformulating the question: Are the v samples of my posterior for condition 1 higher than the samples of my posterior for condition 2?

Let's check for that, we have everything we need in our posterior samples!

In [ ]:

Copied!





# Percent of posterior samples
print(
    f"p(v_1 > v_2) = {np.mean(m1.traces.posterior['v'].values > m2.traces.posterior['v'].values)}"
)
# Percent of posterior samples
print(
    f"p(v_1 > v_2) = {np.mean(m1.traces.posterior['v'].values > m2.traces.posterior['v'].values)}"
)

p(v_1 > v_2) = 0.0

Looks like our inference indicates that there is a $0\%$ chance for v parameter of condition 1 to be higher than the v parameter of condition 2.

Plotting¶

We can also plot the posterior samples to get a sense of the uncertainty.

In [ ]:

Copied!





# Specify a plot with one row and two columns
fig, ax = plt.subplots(1, 1, figsize=(10, 5))

# Plot the posterior samples for condition 1
az.plot_posterior(
    m1.traces.posterior, var_names="v", ax=ax, color="blue", hdi_prob=0.95
)

# Plot the posterior samples for condition 2
az.plot_posterior(m2.traces.posterior, var_names="v", ax=ax, color="red", hdi_prob=0.95)

# Set x-axis limit
ax.set_xlim(0.25, 1.25)

# Create proxy artists for the legend
from matplotlib.lines import Line2D

legend_elements = [
    Line2D([0], [0], color="blue", label="Condition 1"),
    Line2D([0], [0], color="red", label="Condition 2"),
]

# Add a legend
ax.legend(handles=legend_elements)

# Add a title
ax.set_title("Posterior samples for v")

# Add a x-axis label
ax.set_xlabel("v")
# Specify a plot with one row and two columns
fig, ax = plt.subplots(1, 1, figsize=(10, 5))

# Plot the posterior samples for condition 1
az.plot_posterior(
    m1.traces.posterior, var_names="v", ax=ax, color="blue", hdi_prob=0.95
)

# Plot the posterior samples for condition 2
az.plot_posterior(m2.traces.posterior, var_names="v", ax=ax, color="red", hdi_prob=0.95)

# Set x-axis limit
ax.set_xlim(0.25, 1.25)

# Create proxy artists for the legend
from matplotlib.lines import Line2D

legend_elements = [
    Line2D([0], [0], color="blue", label="Condition 1"),
    Line2D([0], [0], color="red", label="Condition 2"),
]

# Add a legend
ax.legend(handles=legend_elements)

# Add a title
ax.set_title("Posterior samples for v")

# Add a x-axis label
ax.set_xlabel("v")

Out[ ]:

Text(0.5, 0, 'v')

No description has been provided for this image

A glance a these posteriors visually corroborates our simple counting analysis. The two posteriors for v are clearly separated.

To appreciate the difference, let us also plot the posteriors of the respective a parameters.

In [ ]:

Copied!





# Specify a plot with one row and two columns
fig, ax = plt.subplots(1, 1, figsize=(10, 5))

# Plot the posterior samples for condition 1
az.plot_posterior(
    m1.traces.posterior, var_names="a", ax=ax, color="blue", hdi_prob=0.95
)

# Plot the posterior samples for condition 2
az.plot_posterior(m2.traces.posterior, var_names="a", ax=ax, color="red", hdi_prob=0.95)

# Set x-axis limit
ax.set_xlim(1.3, 1.7)

# Create proxy artists for the legend
from matplotlib.lines import Line2D

legend_elements = [
    Line2D([0], [0], color="blue", label="Condition 1"),
    Line2D([0], [0], color="red", label="Condition 2"),
]

# Add a legend
ax.legend(handles=legend_elements)

# Add a title
ax.set_title("Posterior samples for a")

# Add a x-axis label
ax.set_xlabel("v")
# Specify a plot with one row and two columns
fig, ax = plt.subplots(1, 1, figsize=(10, 5))

# Plot the posterior samples for condition 1
az.plot_posterior(
    m1.traces.posterior, var_names="a", ax=ax, color="blue", hdi_prob=0.95
)

# Plot the posterior samples for condition 2
az.plot_posterior(m2.traces.posterior, var_names="a", ax=ax, color="red", hdi_prob=0.95)

# Set x-axis limit
ax.set_xlim(1.3, 1.7)

# Create proxy artists for the legend
from matplotlib.lines import Line2D

legend_elements = [
    Line2D([0], [0], color="blue", label="Condition 1"),
    Line2D([0], [0], color="red", label="Condition 2"),
]

# Add a legend
ax.legend(handles=legend_elements)

# Add a title
ax.set_title("Posterior samples for a")

# Add a x-axis label
ax.set_xlabel("v")

Out[ ]:

Text(0.5, 0, 'v')

and correspondingly,

In [ ]:

Copied!

print(
    f"p(a_1 > a_2) = {np.mean(m1.traces.posterior['a'].values > m2.traces.posterior['a'].values)}"
)
print(
    f"p(a_1 > a_2) = {np.mean(m1.traces.posterior['a'].values > m2.traces.posterior['a'].values)}"
)

p(a_1 > a_2) = 0.801

Example 2: Combined Model¶

In [ ]:

Copied!

condition_1["condition"] = "C1"
condition_2["condition"] = "C2"

data = pd.concat([condition_1, condition_2]).reset_index(drop=True)
condition_1["condition"] = "C1"
condition_2["condition"] = "C2"

data = pd.concat([condition_1, condition_2]).reset_index(drop=True)

In [ ]:

Copied!





m_combined = hssm.HSSM(
    model="ddm",
    data=data,
    include=[
        {
            "name": "v",
            "formula": "v ~ 1 + condition",
        }
    ],
)

idata_combined = m_combined.sample(sampler="mcmc", tune=500, draws=500)
m_combined = hssm.HSSM(
    model="ddm",
    data=data,
    include=[
        {
            "name": "v",
            "formula": "v ~ 1 + condition",
        }
    ],
)

idata_combined = m_combined.sample(sampler="mcmc", tune=500, draws=500)

Model initialized successfully.
Using default initvals.

Initializing NUTS using adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [z, t, a, v_Intercept, v_condition]

Output()

Sampling 4 chains for 500 tune and 500 draw iterations (2_000 + 2_000 draws total) took 15 seconds.
100%|██████████| 2000/2000 [00:00<00:00, 2699.50it/s]

Note, now we don't have two distinct idata objects, but instead we have a single idata object with a posterior for v_condition. Let's take a closer look.

In [ ]:

Copied!

m_combined.traces.posterior
m_combined.traces.posterior

Out[ ]:

<xarray.Dataset> Size: 84kB
Dimensions:          (chain: 4, draw: 500, v_condition_dim: 1)
Coordinates:
  * chain            (chain) int64 32B 0 1 2 3
  * draw             (draw) int64 4kB 0 1 2 3 4 5 6 ... 494 495 496 497 498 499
  * v_condition_dim  (v_condition_dim) <U2 8B 'C2'
Data variables:
    z                (chain, draw) float64 16kB 0.485 0.4746 ... 0.4711 0.5008
    a                (chain, draw) float64 16kB 1.46 1.486 1.462 ... 1.457 1.453
    t                (chain, draw) float64 16kB 0.1121 0.1184 ... 0.1261 0.1177
    v_Intercept      (chain, draw) float64 16kB 0.6678 0.5982 ... 0.624 0.6209
    v_condition      (chain, draw, v_condition_dim) float64 16kB 0.3642 ... 0...
Attributes:
    created_at:                  2025-02-16T20:21:04.416816+00:00
    arviz_version:               0.19.0
    inference_library:           pymc
    inference_library_version:   5.19.1
    sampling_time:               12.396317958831787
    tuning_steps:                500
    modeling_interface:          bambi
    modeling_interface_version:  0.15.0

xarray.Dataset

Dimensions:
- chain: 4
- draw: 500
- v_condition_dim: 1
Coordinates: (3)
- chain
  (chain)
  int64
  0 1 2 3
```
array([0, 1, 2, 3])
```
- draw
  (draw)
  int64
  0 1 2 3 4 5 ... 495 496 497 498 499
```
array([  0,   1,   2, ..., 497, 498, 499])
```
- v_condition_dim
  (v_condition_dim)
  <U2
  'C2'
```
array(['C2'], dtype='<U2')
```

Data variables: (5)

(chain, draw)

float64

0.485 0.4746 ... 0.4711 0.5008

array([[0.48504465, 0.47459921, 0.48263437, ..., 0.51856957, 0.47234456,
        0.48974099],
       [0.48741236, 0.49183998, 0.48604471, ..., 0.50843754, 0.47933609,
        0.49490765],
       [0.4969859 , 0.45415678, 0.49671196, ..., 0.47621925, 0.4740677 ,
        0.45323482],
       [0.4791636 , 0.48893406, 0.48301926, ..., 0.50522429, 0.47111393,
        0.50082489]])

(chain, draw)

float64

1.46 1.486 1.462 ... 1.457 1.453

array([[1.46011486, 1.48561764, 1.46230045, ..., 1.45907265, 1.40297597,
        1.50857533],
       [1.43276716, 1.4927821 , 1.49015971, ..., 1.50472936, 1.44479993,
        1.51576519],
       [1.51898333, 1.45648345, 1.42525304, ..., 1.47832583, 1.51056357,
        1.50351921],
       [1.41249163, 1.46813043, 1.49958572, ..., 1.45175015, 1.45700247,
        1.453499  ]])

(chain, draw)

float64

0.1121 0.1184 ... 0.1261 0.1177

array([[0.11209944, 0.11839754, 0.09273701, ..., 0.12320129, 0.13698208,
        0.12163452],
       [0.1493607 , 0.10238487, 0.0876791 , ..., 0.12378699, 0.12613497,
        0.10494174],
       [0.09821629, 0.08315578, 0.13580166, ..., 0.09594787, 0.07434846,
        0.06184704],
       [0.12817953, 0.14456496, 0.11074408, ..., 0.13165505, 0.12613824,
        0.11774354]])

v_Intercept

(chain, draw)

float64

0.6678 0.5982 ... 0.624 0.6209

array([[0.66781977, 0.59817023, 0.57463505, ..., 0.53497314, 0.60374539,
        0.63218083],
       [0.64469056, 0.56995575, 0.6023748 , ..., 0.55816654, 0.66105045,
        0.61612952],
       [0.65376653, 0.66932408, 0.57343211, ..., 0.56230566, 0.6489641 ,
        0.68306346],
       [0.61707905, 0.59860124, 0.65748562, ..., 0.59916945, 0.623997  ,
        0.62091879]])

v_condition

(chain, draw, v_condition_dim)

float64

0.3642 0.3753 ... 0.4057 0.2887

array([[[0.36418892],
        [0.37534147],
        [0.43011859],
        ...,
        [0.48740229],
        [0.32755279],
        [0.34391929]],

       [[0.36340026],
        [0.34181001],
        [0.33541986],
        ...,
        [0.31666528],
        [0.28445689],
        [0.40995799]],

       [[0.32477034],
        [0.3420066 ],
        [0.34262551],
        ...,
        [0.39151984],
        [0.41306592],
        [0.28285667]],

       [[0.4190528 ],
        [0.41912041],
        [0.3021753 ],
        ...,
        [0.37824254],
        [0.40573261],
        [0.28871663]]])

Indexes: (3)

chain

PandasIndex

PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='chain'))

draw

PandasIndex

PandasIndex(Index([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,
       ...
       490, 491, 492, 493, 494, 495, 496, 497, 498, 499],
      dtype='int64', name='draw', length=500))

v_condition_dim

PandasIndex

PandasIndex(Index(['C2'], dtype='object', name='v_condition_dim'))

Attributes: (8)
created_at :
2025-02-16T20:21:04.416816+00:00
arviz_version :
0.19.0
inference_library :
pymc
inference_library_version :
5.19.1
sampling_time :
12.396317958831787
tuning_steps :
500
modeling_interface :
bambi
modeling_interface_version :
0.15.0

Under the hood, Bambi created a model with a dummy variable for the C2 condition (the C1 condition represents the Intercept).

So what do we need to test here...

We still ask the same question: Is the v_condition[C2] posterior higher than the v_condition[C1] posterior?

But now we only need to check if our v_condition[C2] variable is above 0, since it represents the offset from the Intercept directly.

Let's check:

In [ ]:

Copied!

print(
    f"p(v_condition[C2] > v_condition[C1]) = {np.mean(idata_combined.posterior['v_condition'].values > 0)}"
)
print(
    f"p(v_condition[C2] > v_condition[C1]) = {np.mean(idata_combined.posterior['v_condition'].values > 0)}"
)

p(v_condition[C2] > v_condition[C1]) = 1.0

or visually,

In [ ]:

Copied!

az.plot_posterior(idata_combined.posterior, var_names="v_condition", hdi_prob=0.95)
az.plot_posterior(idata_combined.posterior, var_names="v_condition", hdi_prob=0.95)

Out[ ]:

<Axes: title={'center': 'v_condition\nC2'}>

You can use this approach to test any number of complex statements about your parameters. There will essentially always be a way to turn your question into a simple comparison of the posterior samples, a simple counting problem.