DDM with hierarchical regressions¶

This tutorial walks you through fitting hierarchical Drift Diffusion Models (DDMs) with HSSM. If you're new to HSSM, please click this link. For background on hierarchical modeling in HSSM, please click this link. We'll move through four progressively more complex scenarios, each demonstrating a different real-world experimental design and the kind of regression formula you'd use to capture it.

Here are all 4 scenarios:

Variant	Scenario	v formula
A	A within-subject experiment where each trial varies in difficulty.	v ~ 0 + difficulty + (1 + difficulty \| participant_id)
B	Same as A, but participants also differ in age (a between-subject factor). Age affects drift rate on top of the within-subject difficulty effect.	v ~ 0 + age + difficulty + (1 + difficulty \| participant_id)
C	Same as B, but with an interaction: the effect of difficulty on drift rate also affects age.	v ~ 0 + age + difficulty + age:difficulty + (1 + difficulty \| participant_id)
D	Same v regression as A, but now every DDM parameter (a, z, t) gets its own per-participant random intercept — this is "fully hierarchical" case.	(v same as A; full hierarchy on a, z, t)

In Variants A–C the other parameters (a, z, t) are estimated globally (param ~ 1) — this keeps the focus on the v regression and makes sampling fast. Variant D shows what happens when we add hierarchy to every parameter, demonstrating the fully hierarchical case.

Note¶

Each variant follows the same pattern: simulate synthetic data with known parameter values, fit the model, then check whether we recovered the true parameters. This is parameter recovery — the standard sanity check before applying a model to real data.

If you already have your own data to estimate, you can skip the simulation steps and go straight to the model fitting. But running the simulation to recovery loop on synthetic data first is always a good idea — it confirms your modelling pipeline can actually recover the parameters it's supposed to estimate before you trust its inferences on real data.

Layout of each variant¶

1. Pick ground-truth coefficients and simulate synthetic data with the do-operator
2. Fit a HSSM model to that data
3. Model graph and posterior summary
4. Check parameter recovery and posterior predictive fit

In [77]:

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

import hssm

rng = np.random.default_rng(0)

sample_kwargs = dict(
    draws=1000,
    tune=500,
    chains=2,
    target_accept=0.9,
    omit_offsets=True,
    idata_kwargs={"log_likelihood": False},
)

import arviz as az import matplotlib.pyplot as plt import numpy as np import pandas as pd import pymc as pm import hssm rng = np.random.default_rng(0) sample_kwargs = dict( draws=1000, tune=500, chains=2, target_accept=0.9, omit_offsets=True, idata_kwargs={"log_likelihood": False}, )

Covariates¶

We'll have 10 participants doing 250 trials each, giving us 2500 trials total.

Two covariates are tracked throughout:

• difficulty — varies trial-by-trial (within-subject)
• age — one value per participant (between-subject)

difficulty is simulated on a 0-10 scale, and age is represented as participant age in years. The rt and response columns start as placeholders and get filled in after simulation.

In this synthetic setup, larger difficulty values correspond to trials with stronger or clearer evidence, so drift rate increases as difficulty increases. age is represented as participant age in years; older participants are assumed to have higher drift rates.

In [ ]:

n_participants = 10
n_trials = 250

pid = np.repeat(np.arange(n_participants), n_trials)

# Trial-level difficulty on a readable 0-10 scale.
difficulty = rng.uniform(0, 10, size=n_participants * n_trials)

# Participant-level age in years.
age_by_p = rng.uniform(20, 75, size=n_participants)
age = age_by_p[pid]

covariates = pd.DataFrame(
    {
        "participant_id": pid,
        "difficulty": difficulty,
        "age": age,
        "rt": np.nan,
        "response": np.nan,
    }
)
covariates.head()

n_participants = 10 n_trials = 250 pid = np.repeat(np.arange(n_participants), n_trials) # Trial-level difficulty on a readable 0-10 scale. difficulty = rng.uniform(0, 10, size=n_participants * n_trials) # Participant-level age in years. age_by_p = rng.uniform(20, 75, size=n_participants) age = age_by_p[pid] covariates = pd.DataFrame( { "participant_id": pid, "difficulty": difficulty, "age": age, "rt": np.nan, "response": np.nan, } ) covariates.head()

Out[ ]:

	participant_id	difficulty	age	rt	response
0	0	6.369617	60.357299	NaN	NaN
1	0	2.697867	60.357299	NaN	NaN
2	0	0.409735	60.357299	NaN	NaN
3	0	0.165276	60.357299	NaN	NaN
4	0	8.132702	60.357299	NaN	NaN

Simulating synthetic data with the do-operator¶

What is the do-operator doing?¶

It intervenes on the model's parameters by fixing them to specific values, then sample the likelihood to generate fresh rt / response data. The result is synthetic data that we know was generated from that exact parameter setting.

In code, this looks like model.sample_do(params={...}). We pass it a dictionary of parameter values and it returns a fresh dataset.

For a deeper walkthrough of sample_do, see the HSSM do-operator tutorial.

Variant A: within-subject¶

Each trial varies in difficulty. Each participant has their own baseline drift rate and their own sensitivity to difficulty.

Drift rate v depends on difficulty, and each participant gets their own intercept and slope.

The other parameters (a, z, t) are kept as global scalars - no per-participant variation, no regression.

Step 1: Simulate the data¶

We pick ground-truth values for the regression coefficients (e.g. beta_v_difficulty = 0.08), determine the per-participant random effects, and pass everything to the model via the do-operator. Outcomes are synthetic rt/response data that we know was generated from this exact parameter setting.

Because a, z, t are global scalars in this variant, we simulate them as scalars too.

In [79]:

# true coefficients for variant A
beta_v_difficulty = 0.08
sigma_v_int, sigma_v_slope = 0.3, 0.2

v_int_p = rng.normal(0, sigma_v_int, n_participants)
v_slope_p = rng.normal(0, sigma_v_slope, n_participants)
a_true, z_true, t_true = 1.2, 0.5, 0.25

# Noncentered offsets for the random effects
v_int_offset = v_int_p / sigma_v_int
v_slope_offset = v_slope_p / sigma_v_slope

# Build HSSM model with the regression formula
dummy_data = covariates.copy()
dummy_data["rt"] = 1.0
dummy_data["response"] = 1.0

model_A_gen = hssm.HSSM(
    data=dummy_data,
    model="ddm",
    global_formula="y ~ 1",
    include=[
        {
            "name": "v",
            "formula": "v ~ 0 + difficulty + (1 + difficulty | participant_id)",
        }
    ],
    noncentered=True,
)

# true coefficients for variant A beta_v_difficulty = 0.08 sigma_v_int, sigma_v_slope = 0.3, 0.2 v_int_p = rng.normal(0, sigma_v_int, n_participants) v_slope_p = rng.normal(0, sigma_v_slope, n_participants) a_true, z_true, t_true = 1.2, 0.5, 0.25 # Noncentered offsets for the random effects v_int_offset = v_int_p / sigma_v_int v_slope_offset = v_slope_p / sigma_v_slope # Build HSSM model with the regression formula dummy_data = covariates.copy() dummy_data["rt"] = 1.0 dummy_data["response"] = 1.0 model_A_gen = hssm.HSSM( data=dummy_data, model="ddm", global_formula="y ~ 1", include=[ { "name": "v", "formula": "v ~ 0 + difficulty + (1 + difficulty | participant_id)", } ], noncentered=True, )

No common intercept. Bounds for parameter v is not applied due to a current limitation of Bambi. This will change in the future.
Model initialized successfully.

/Users/mayan/HSSM/src/hssm/hssm.py:504: UserWarning: You set choices to be (-1, 1), but [-1] are missing from your dataset.
  self._post_check_data_sanity()

To see how HSSM relates to PyMC, please click this link.

In [80]:

# Use do-operator
synth_idata_A, synth_model_A = model_A_gen.sample_do(
    params={
        "v_difficulty": beta_v_difficulty,
        "v_1|participant_id_mu": 0.0,
        "v_1|participant_id_sigma": sigma_v_int,
        "v_1|participant_id_offset": v_int_offset,
        "v_difficulty|participant_id_mu": 0.0,
        "v_difficulty|participant_id_sigma": sigma_v_slope,
        "v_difficulty|participant_id_offset": v_slope_offset,
        "a_Intercept": a_true,
        "z_Intercept": z_true,
        "t_Intercept": t_true,
    },
    draws=1,
    var_names=["rt,response"],
    return_model=True,
)

# Use do-operator synth_idata_A, synth_model_A = model_A_gen.sample_do( params={ "v_difficulty": beta_v_difficulty, "v_1|participant_id_mu": 0.0, "v_1|participant_id_sigma": sigma_v_int, "v_1|participant_id_offset": v_int_offset, "v_difficulty|participant_id_mu": 0.0, "v_difficulty|participant_id_sigma": sigma_v_slope, "v_difficulty|participant_id_offset": v_slope_offset, "a_Intercept": a_true, "z_Intercept": z_true, "t_Intercept": t_true, }, draws=1, var_names=["rt,response"], return_model=True, )

Sampling: [rt,response]

The graph below shows the response-generating part of the synthetic model after the do-operator has fixed the parameters used for simulation. For details on how PyMC renders model graphs, see the pm.model_to_graphviz documentation.

In [81]:

pm.model_to_graphviz(synth_model_A, var_names=["rt,response"])

pm.model_to_graphviz(synth_model_A, var_names=["rt,response"])

Out[81]:

The below cell converts simulated prior predictive outputs from the HSSM model into a DataFrame and extracts the generated reaction times (rt) and responses. It then merges these simulated outcomes with the original covariates to create a complete synthetic dataset and saves it as a CSV file.

In [82]:

synth_df_A = hssm.utils.predictive_idata_to_dataframe(
    synth_idata_A, predictive_group="prior_predictive"
)

data_A = covariates.copy()
data_A["rt"] = synth_df_A["rt"].values
data_A["response"] = synth_df_A["response"].values

data_A.to_csv("sim_variant_A.csv", index=False)
data_A.head()

synth_df_A = hssm.utils.predictive_idata_to_dataframe( synth_idata_A, predictive_group="prior_predictive" ) data_A = covariates.copy() data_A["rt"] = synth_df_A["rt"].values data_A["response"] = synth_df_A["response"].values data_A.to_csv("sim_variant_A.csv", index=False) data_A.head()

Out[82]:

	participant_id	difficulty	age	rt	response
0	0	6.369617	60.357299	0.684816	1.0
1	0	2.697867	60.357299	2.162073	1.0
2	0	0.409735	60.357299	1.049083	-1.0
3	0	0.165276	60.357299	2.346251	-1.0
4	0	8.132702	60.357299	0.543396	1.0

Step 2: Fit a fresh HSSM model¶

Now we pretend we don't know the true parameters and try to recover them from the data. We reload the saved CSV, build a new HSSM model with the same regression structure, and explore the posterior.

In [84]:

data_A = pd.read_csv("sim_variant_A.csv")

# Bin difficulty into 3 levels (low/med/high) for PPC plots
data_A["difficulty_level"] = pd.qcut(
    data_A["difficulty"], q=3, labels=["low", "med", "high"]
)

model_A = hssm.HSSM(
    data=data_A,
    model="ddm",
    global_formula="y ~ 1",
    include=[
        {
            "name": "v",
            "formula": "v ~ 0 + difficulty + (1 + difficulty | participant_id)",
        }
    ],
    noncentered=True,
)

data_A = pd.read_csv("sim_variant_A.csv") # Bin difficulty into 3 levels (low/med/high) for PPC plots data_A["difficulty_level"] = pd.qcut( data_A["difficulty"], q=3, labels=["low", "med", "high"] ) model_A = hssm.HSSM( data=data_A, model="ddm", global_formula="y ~ 1", include=[ { "name": "v", "formula": "v ~ 0 + difficulty + (1 + difficulty | participant_id)", } ], noncentered=True, )

No common intercept. Bounds for parameter v is not applied due to a current limitation of Bambi. This will change in the future.
Model initialized successfully.

Sanity check before sampling¶

print(model) shows every parameter's formula in text. The PyMC/Graphviz call below draws the response-generating subgraph, which is useful for checking that the regression structure and participant-level hierarchy feed into rt,response as expected. For details on PyMC model graphs, see the pm.model_to_graphviz documentation.

For more example graphs across HSSM models, see the Scientific Workflow tutorial.

In [85]:

print(model_A)
pm.model_to_graphviz(model_A.pymc_model, var_names=["rt,response"])

print(model_A) pm.model_to_graphviz(model_A.pymc_model, var_names=["rt,response"])

Hierarchical Sequential Sampling Model
Model: ddm

Response variable: rt,response
Likelihood: analytical
Observations: 2500

Parameters:

v:
    Formula: v ~ 0 + difficulty + (1 + difficulty | participant_id)
    Priors:
        v_difficulty ~ Normal(mu: 0.0, sigma: 0.25)
        v_1|participant_id ~ Normal(mu: Normal(mu: 2.0, sigma: 3.0), sigma: HalfNormal(sigma: 2.0))
        v_difficulty|participant_id ~ Normal(mu: Normal(mu: 0.0, sigma: 0.25), sigma: Weibull(alpha: 1.5, beta: 0.3))
    Link: identity
    Explicit bounds: (-inf, inf)

a:
    Formula: a ~ 1
    Priors:
        a_Intercept ~ Gamma(mu: 1.5, sigma: 0.75)
    Link: identity
    Explicit bounds: (0.0, inf)

z:
    Formula: z ~ 1
    Priors:
        z_Intercept ~ Beta(alpha: 10.0, beta: 10.0)
    Link: identity
    Explicit bounds: (0.0, 1.0)

t:
    Formula: t ~ 1
    Priors:
        t_Intercept ~ Gamma(mu: 0.2, sigma: 0.2)
    Link: identity
    Explicit bounds: (0.0, inf)


Lapse probability: 0.05
Lapse distribution: Uniform(lower: 0.0, upper: 20.0)

Out[85]:

In [86]:

idata_A = model_A.sample(**sample_kwargs)

idata_A = model_A.sample(**sample_kwargs)

Using default initvals. 

Initializing NUTS using adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [v_difficulty, v_1|participant_id_mu, v_1|participant_id_sigma, v_1|participant_id_offset, v_difficulty|participant_id_mu, v_difficulty|participant_id_sigma, v_difficulty|participant_id_offset, a_Intercept, z_Intercept, t_Intercept]
/Users/mayan/.local/share/uv/python/cpython-3.12.12-macos-aarch64-none/lib/python3.12/multiprocessing/popen_fork.py:66: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.
  self.pid = os.fork()

Output()

/Users/mayan/.local/share/uv/python/cpython-3.12.12-macos-aarch64-none/lib/python3.12/multiprocessing/popen_fork.py:66: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.
  self.pid = os.fork()

Sampling 2 chains for 500 tune and 1_000 draw iterations (1_000 + 2_000 draws total) took 252 seconds.
There was 1 divergence after tuning. Increase `target_accept` or reparameterize.
We recommend running at least 4 chains for robust computation of convergence diagnostics

Step 3: Inspect what got estimated¶

First look at all estimated variables (no var_names filter) — this is how you discover the participant-wise parameters needed for recovery checks.

In [87]:

az.summary(idata_A)

az.summary(idata_A)

Out[87]:

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
v_difficulty|participant_id_mu	0.001	0.245	-0.442	0.478	0.004	0.005	3073.0	1667.0	1.00
v_difficulty|participant_id[0]	-0.159	0.072	-0.292	-0.029	0.003	0.002	664.0	669.0	1.00
v_difficulty|participant_id[1]	0.266	0.078	0.122	0.416	0.003	0.002	747.0	782.0	1.00
v_difficulty|participant_id[2]	0.149	0.077	0.004	0.293	0.003	0.002	663.0	681.0	1.01
v_difficulty|participant_id[3]	-0.078	0.074	-0.220	0.055	0.003	0.002	654.0	741.0	1.00
v_difficulty|participant_id[4]	0.216	0.078	0.056	0.349	0.003	0.002	656.0	694.0	1.00
v_difficulty|participant_id[5]	0.160	0.078	0.013	0.306	0.003	0.002	668.0	807.0	1.00
v_difficulty|participant_id[6]	0.222	0.077	0.078	0.370	0.003	0.002	662.0	717.0	1.00
v_difficulty|participant_id[7]	-0.110	0.075	-0.232	0.040	0.003	0.002	743.0	771.0	1.00
v_difficulty|participant_id[8]	-0.056	0.076	-0.201	0.080	0.003	0.002	722.0	667.0	1.00
v_difficulty|participant_id[9]	-0.100	0.074	-0.233	0.041	0.003	0.002	669.0	674.0	1.00
v_1|participant_id_mu	1.996	2.924	-3.510	7.447	0.057	0.062	2672.0	1589.0	1.00
v_1|participant_id_sigma	0.329	0.100	0.167	0.515	0.004	0.003	732.0	1165.0	1.00
z_Intercept	0.507	0.014	0.481	0.532	0.000	0.000	1801.0	1319.0	1.00
v_difficulty|participant_id_sigma	0.198	0.058	0.114	0.307	0.003	0.002	585.0	709.0	1.00
t_Intercept	0.253	0.004	0.246	0.259	0.000	0.000	1941.0	1546.0	1.00
v_difficulty	0.757	0.070	0.618	0.878	0.003	0.002	617.0	624.0	1.01
v_1|participant_id[0]	-0.490	0.140	-0.749	-0.234	0.003	0.003	2008.0	1612.0	1.00
v_1|participant_id[1]	-0.201	0.154	-0.485	0.077	0.003	0.003	2241.0	1273.0	1.00
v_1|participant_id[2]	0.067	0.141	-0.202	0.326	0.003	0.003	2496.0	1439.0	1.00
v_1|participant_id[3]	0.004	0.129	-0.237	0.248	0.003	0.004	2184.0	1336.0	1.00
v_1|participant_id[4]	0.004	0.142	-0.254	0.272	0.003	0.004	2444.0	1178.0	1.00
v_1|participant_id[5]	0.190	0.144	-0.094	0.455	0.003	0.004	2610.0	1435.0	1.00
v_1|participant_id[6]	-0.108	0.137	-0.362	0.145	0.003	0.003	2777.0	1754.0	1.00
v_1|participant_id[7]	-0.538	0.128	-0.776	-0.302	0.003	0.003	2314.0	1356.0	1.00
v_1|participant_id[8]	0.195	0.124	-0.044	0.418	0.002	0.003	2791.0	1597.0	1.00
v_1|participant_id[9]	0.022	0.133	-0.247	0.256	0.003	0.003	2631.0	1536.0	1.00
a_Intercept	1.213	0.025	1.168	1.262	0.000	0.001	2691.0	1202.0	1.00

Now filter to the key global coefficients.

In [88]:

az.summary(idata_A, var_names=["~offset", "~participant_id"], filter_vars="regex")

az.summary(idata_A, var_names=["~offset", "~participant_id"], filter_vars="regex")

/Users/mayan/HSSM/.venv/lib/python3.12/site-packages/arviz/utils.py:146: UserWarning: Items starting with ~: ['offset'] have not been found and will be ignored
  warnings.warn(

Out[88]:

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
z_Intercept	0.507	0.014	0.481	0.532	0.000	0.000	1801.0	1319.0	1.00
t_Intercept	0.253	0.004	0.246	0.259	0.000	0.000	1941.0	1546.0	1.00
v_difficulty	0.757	0.070	0.618	0.878	0.003	0.002	617.0	624.0	1.01
a_Intercept	1.213	0.025	1.168	1.262	0.000	0.001	2691.0	1202.0	1.00

az.plot_trace shows two views of each parameter side by side:

• Left column — the posterior density, one line per chain.
• Right column — the trace: parameter value across posterior draws, shown separately for each chain.

For readability, the trace plots below focus on global/group-level parameters and omit participant-specific random effects; those participant-level effects are checked separately in the recovery plots.

For more on plotting in HSSM, please click this link.

In [89]:

az.plot_trace(idata_A, var_names=["~offset", "~participant_id"], filter_vars="regex")
plt.tight_layout()

az.plot_trace(idata_A, var_names=["~offset", "~participant_id"], filter_vars="regex") plt.tight_layout()

/Users/mayan/HSSM/.venv/lib/python3.12/site-packages/arviz/utils.py:146: UserWarning: Items starting with ~: ['offset'] have not been found and will be ignored
  warnings.warn(

Step 4: Posterior predictive check¶

We sample many parameter sets from the posterior, forward-simulate rt/response data from each, and compare to what we observed. If the predicted distribution overlaps the observed one, the model is doing its job. If they diverge, the model is missing something.

In [90]:

ax = hssm.plotting.plot_predictive(model_A, col_wrap=5)
ax.set_xlim(-7, 7)

ax = hssm.plotting.plot_predictive(model_A, col_wrap=5) ax.set_xlim(-7, 7)

No posterior_predictive samples found. Generating posterior_predictive samples using the provided InferenceData object and the original data. This will modify the provided InferenceData object, or if not provided, the traces object stored inside the model.

Out[90]:

(-7.0, 7.0)

In [91]:

g = hssm.plotting.plot_predictive(model_A, col="participant_id", col_wrap=5)
g.set(xlim=(-7, 7))

g = hssm.plotting.plot_predictive(model_A, col="participant_id", col_wrap=5) g.set(xlim=(-7, 7))

Out[91]:

<seaborn.axisgrid.FacetGrid at 0x140ea92b0>

In [92]:

# Participant-wise parameter recovery for v (random intercepts and slopes)
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 5), sharey=True)

for ax, (param_name, truth, label) in zip(
    axes,
    [
        ("v_1|participant_id", v_int_p, "v random intercept"),
        ("v_difficulty|participant_id", v_slope_p, "v random slope on difficulty"),
    ],
):
    post = idata_A.posterior[param_name]
    mean = post.mean(dim=["chain", "draw"]).values
    hdi = az.hdi(post, hdi_prob=0.95)[param_name].values

    yvals = np.arange(n_participants)
    ax.hlines(yvals, hdi[:, 0], hdi[:, 1], color="blue", alpha=0.5)
    ax.plot(mean, yvals, "|", color="blue", label="posterior mean")
    ax.plot(truth, yvals, "x", color="red", label="ground truth")
    ax.set_title(label)
    ax.set_yticks(yvals)
    ax.set_ylabel("participant_id")
    ax.legend()

plt.tight_layout()

# Participant-wise parameter recovery for v (random intercepts and slopes) fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 5), sharey=True) for ax, (param_name, truth, label) in zip( axes, [ ("v_1|participant_id", v_int_p, "v random intercept"), ("v_difficulty|participant_id", v_slope_p, "v random slope on difficulty"), ], ): post = idata_A.posterior[param_name] mean = post.mean(dim=["chain", "draw"]).values hdi = az.hdi(post, hdi_prob=0.95)[param_name].values yvals = np.arange(n_participants) ax.hlines(yvals, hdi[:, 0], hdi[:, 1], color="blue", alpha=0.5) ax.plot(mean, yvals, "|", color="blue", label="posterior mean") ax.plot(truth, yvals, "x", color="red", label="ground truth") ax.set_title(label) ax.set_yticks(yvals) ax.set_ylabel("participant_id") ax.legend() plt.tight_layout()

In [93]:

# Posterior predictive: split by difficulty level
g = hssm.plotting.plot_predictive(model_A, col="difficulty_level", col_wrap=3)
g.set(xlim=(-7, 7))
g.fig.suptitle("")

# Posterior predictive: split by difficulty level g = hssm.plotting.plot_predictive(model_A, col="difficulty_level", col_wrap=3) g.set(xlim=(-7, 7)) g.fig.suptitle("")

Out[93]:

Text(0.5, 0.98, '')

Variant B: between and within¶

Same within-subject difficulty effect as in Variant A, but now we also let age (a between-subject covariate) affect drift rate.

The a, z, t parameters remain global scalars.

Step 1: Simulate the data¶

In [102]:

# true coefficients for variant B
beta_v_difficulty = 0.08
beta_v_age = 0.005

v_int_p = rng.normal(0, sigma_v_int, n_participants)
v_slope_p = rng.normal(0, sigma_v_slope, n_participants)

v_int_offset = v_int_p / sigma_v_int
v_slope_offset = v_slope_p / sigma_v_slope

# Build HSSM model with the regression formula
dummy_data = covariates.copy()
dummy_data["rt"] = 1.0
dummy_data["response"] = 1.0

model_B_gen = hssm.HSSM(
    data=dummy_data,
    model="ddm",
    global_formula="y ~ 1",
    include=[
        {
            "name": "v",
            "formula": "v ~ 0 + age + difficulty + (1 + difficulty | participant_id)",
        }
    ],
    noncentered=True,
)

# Use do-operator: intervene on regression coefficients and random effects
synth_idata_B, synth_model_B = model_B_gen.sample_do(
    params={
        "v_age": beta_v_age,
        "v_difficulty": beta_v_difficulty,
        "v_1|participant_id_mu": 0.0,
        "v_1|participant_id_sigma": sigma_v_int,
        "v_1|participant_id_offset": v_int_offset,
        "v_difficulty|participant_id_mu": 0.0,
        "v_difficulty|participant_id_sigma": sigma_v_slope,
        "v_difficulty|participant_id_offset": v_slope_offset,
        "a_Intercept": a_true,
        "z_Intercept": z_true,
        "t_Intercept": t_true,
    },
    draws=1,
    var_names=["rt,response"],
    return_model=True,
)

pm.model_to_graphviz(synth_model_B, var_names=["rt,response"])

# true coefficients for variant B beta_v_difficulty = 0.08 beta_v_age = 0.005 v_int_p = rng.normal(0, sigma_v_int, n_participants) v_slope_p = rng.normal(0, sigma_v_slope, n_participants) v_int_offset = v_int_p / sigma_v_int v_slope_offset = v_slope_p / sigma_v_slope # Build HSSM model with the regression formula dummy_data = covariates.copy() dummy_data["rt"] = 1.0 dummy_data["response"] = 1.0 model_B_gen = hssm.HSSM( data=dummy_data, model="ddm", global_formula="y ~ 1", include=[ { "name": "v", "formula": "v ~ 0 + age + difficulty + (1 + difficulty | participant_id)", } ], noncentered=True, ) # Use do-operator: intervene on regression coefficients and random effects synth_idata_B, synth_model_B = model_B_gen.sample_do( params={ "v_age": beta_v_age, "v_difficulty": beta_v_difficulty, "v_1|participant_id_mu": 0.0, "v_1|participant_id_sigma": sigma_v_int, "v_1|participant_id_offset": v_int_offset, "v_difficulty|participant_id_mu": 0.0, "v_difficulty|participant_id_sigma": sigma_v_slope, "v_difficulty|participant_id_offset": v_slope_offset, "a_Intercept": a_true, "z_Intercept": z_true, "t_Intercept": t_true, }, draws=1, var_names=["rt,response"], return_model=True, ) pm.model_to_graphviz(synth_model_B, var_names=["rt,response"])

No common intercept. Bounds for parameter v is not applied due to a current limitation of Bambi. This will change in the future.

/Users/mayan/HSSM/src/hssm/hssm.py:504: UserWarning: You set choices to be (-1, 1), but [-1] are missing from your dataset.
  self._post_check_data_sanity()

Model initialized successfully.

Sampling: [rt,response]

Out[102]:

In [103]:

synth_df_B = hssm.utils.predictive_idata_to_dataframe(
    synth_idata_B, predictive_group="prior_predictive"
)

data_B = covariates.copy()
data_B["rt"] = synth_df_B["rt"].values
data_B["response"] = synth_df_B["response"].values
data_B.to_csv("sim_variant_B.csv", index=False)

synth_df_B = hssm.utils.predictive_idata_to_dataframe( synth_idata_B, predictive_group="prior_predictive" ) data_B = covariates.copy() data_B["rt"] = synth_df_B["rt"].values data_B["response"] = synth_df_B["response"].values data_B.to_csv("sim_variant_B.csv", index=False)

Step 2: Fit a fresh HSSM model¶

In [104]:

data_B = pd.read_csv("sim_variant_B.csv")

# Bin difficulty (3 levels) and age (3 levels) for faceted PPC plots
data_B["difficulty_level"] = pd.qcut(
    data_B["difficulty"], q=3, labels=["low", "med", "high"]
)
data_B["age_level"] = pd.qcut(data_B["age"], q=3, labels=["young", "middle", "older"])

model_B = hssm.HSSM(
    data=data_B,
    model="ddm",
    global_formula="y ~ 1",
    include=[
        {
            "name": "v",
            "formula": "v ~ 0 + age + difficulty + (1 + difficulty | participant_id)",
        }
    ],
    noncentered=True,
)

data_B = pd.read_csv("sim_variant_B.csv") # Bin difficulty (3 levels) and age (3 levels) for faceted PPC plots data_B["difficulty_level"] = pd.qcut( data_B["difficulty"], q=3, labels=["low", "med", "high"] ) data_B["age_level"] = pd.qcut(data_B["age"], q=3, labels=["young", "middle", "older"]) model_B = hssm.HSSM( data=data_B, model="ddm", global_formula="y ~ 1", include=[ { "name": "v", "formula": "v ~ 0 + age + difficulty + (1 + difficulty | participant_id)", } ], noncentered=True, )

No common intercept. Bounds for parameter v is not applied due to a current limitation of Bambi. This will change in the future.

Model initialized successfully.

In [105]:

print(model_B)
pm.model_to_graphviz(model_B.pymc_model, var_names=["rt,response"])

print(model_B) pm.model_to_graphviz(model_B.pymc_model, var_names=["rt,response"])

Hierarchical Sequential Sampling Model
Model: ddm

Response variable: rt,response
Likelihood: analytical
Observations: 2500

Parameters:

v:
    Formula: v ~ 0 + age + difficulty + (1 + difficulty | participant_id)
    Priors:
        v_age ~ Normal(mu: 0.0, sigma: 0.25)
        v_difficulty ~ Normal(mu: 0.0, sigma: 0.25)
        v_1|participant_id ~ Normal(mu: Normal(mu: 2.0, sigma: 3.0), sigma: HalfNormal(sigma: 2.0))
        v_difficulty|participant_id ~ Normal(mu: Normal(mu: 0.0, sigma: 0.25), sigma: Weibull(alpha: 1.5, beta: 0.3))
    Link: identity
    Explicit bounds: (-inf, inf)

a:
    Formula: a ~ 1
    Priors:
        a_Intercept ~ Gamma(mu: 1.5, sigma: 0.75)
    Link: identity
    Explicit bounds: (0.0, inf)

z:
    Formula: z ~ 1
    Priors:
        z_Intercept ~ Beta(alpha: 10.0, beta: 10.0)
    Link: identity
    Explicit bounds: (0.0, 1.0)

t:
    Formula: t ~ 1
    Priors:
        t_Intercept ~ Gamma(mu: 0.2, sigma: 0.2)
    Link: identity
    Explicit bounds: (0.0, inf)


Lapse probability: 0.05
Lapse distribution: Uniform(lower: 0.0, upper: 20.0)

Out[105]:

In [106]:

idata_B = model_B.sample(**sample_kwargs)

idata_B = model_B.sample(**sample_kwargs)

Using default initvals. 

Initializing NUTS using adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [v_age, v_difficulty, v_1|participant_id_mu, v_1|participant_id_sigma, v_1|participant_id_offset, v_difficulty|participant_id_mu, v_difficulty|participant_id_sigma, v_difficulty|participant_id_offset, a_Intercept, z_Intercept, t_Intercept]
/Users/mayan/.local/share/uv/python/cpython-3.12.12-macos-aarch64-none/lib/python3.12/multiprocessing/popen_fork.py:66: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.
  self.pid = os.fork()

Output()

/Users/mayan/.local/share/uv/python/cpython-3.12.12-macos-aarch64-none/lib/python3.12/multiprocessing/popen_fork.py:66: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.
  self.pid = os.fork()

Sampling 2 chains for 500 tune and 1_000 draw iterations (1_000 + 2_000 draws total) took 190 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
100%|██████████| 2000/2000 [00:01<00:00, 1275.31it/s]

Step 3: Inspect what got estimated¶

In [154]:

az.summary(idata_B)

az.summary(idata_B)

Out[154]:

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
v_difficulty|participant_id_mu	-0.001	0.254	-0.470	0.487	0.006	0.006	1969.0	1464.0	1.00
v_age	0.006	0.002	0.003	0.010	0.000	0.000	911.0	1119.0	1.01
v_difficulty|participant_id[0]	-0.073	0.097	-0.247	0.113	0.005	0.003	330.0	687.0	1.01
v_difficulty|participant_id[1]	-0.415	0.099	-0.595	-0.227	0.005	0.003	351.0	794.0	1.01
v_difficulty|participant_id[2]	-0.153	0.098	-0.337	0.029	0.005	0.003	315.0	804.0	1.01
v_difficulty|participant_id[3]	-0.245	0.097	-0.417	-0.052	0.005	0.003	327.0	692.0	1.01
v_difficulty|participant_id[4]	0.380	0.099	0.185	0.557	0.005	0.003	351.0	775.0	1.01
v_difficulty|participant_id[5]	0.207	0.098	0.018	0.388	0.005	0.003	340.0	750.0	1.01
v_difficulty|participant_id[6]	0.061	0.098	-0.113	0.252	0.005	0.003	336.0	773.0	1.01
v_difficulty|participant_id[7]	0.477	0.100	0.290	0.657	0.005	0.003	362.0	866.0	1.01
v_difficulty|participant_id[8]	0.235	0.097	0.065	0.431	0.005	0.003	334.0	695.0	1.01
v_difficulty|participant_id[9]	-0.267	0.097	-0.444	-0.082	0.005	0.003	314.0	683.0	1.01
v_1|participant_id_mu	2.053	2.984	-3.460	7.537	0.071	0.069	1762.0	1328.0	1.00
v_1|participant_id_sigma	0.244	0.091	0.095	0.402	0.003	0.003	743.0	965.0	1.00
z_Intercept	0.499	0.008	0.484	0.515	0.000	0.000	1881.0	1477.0	1.00
v_difficulty|participant_id_sigma	0.319	0.075	0.191	0.460	0.003	0.002	577.0	1009.0	1.00
t_Intercept	0.259	0.005	0.248	0.268	0.000	0.000	1552.0	1361.0	1.00
v_difficulty	0.094	0.095	-0.083	0.270	0.005	0.003	294.0	647.0	1.01
v_1|participant_id[0]	0.075	0.144	-0.190	0.371	0.004	0.003	1136.0	1366.0	1.01
v_1|participant_id[1]	-0.219	0.128	-0.458	0.016	0.003	0.002	1594.0	1453.0	1.00
v_1|participant_id[2]	0.279	0.167	-0.015	0.603	0.006	0.004	922.0	1149.0	1.01
v_1|participant_id[3]	0.004	0.135	-0.229	0.265	0.004	0.003	1131.0	1316.0	1.01
v_1|participant_id[4]	-0.104	0.144	-0.382	0.165	0.004	0.003	1371.0	1465.0	1.01
v_1|participant_id[5]	-0.230	0.124	-0.439	0.022	0.003	0.003	1972.0	1562.0	1.00
v_1|participant_id[6]	0.188	0.129	-0.034	0.451	0.003	0.003	1730.0	1482.0	1.00
v_1|participant_id[7]	-0.143	0.155	-0.457	0.121	0.004	0.003	1307.0	1332.0	1.00
v_1|participant_id[8]	0.016	0.114	-0.209	0.221	0.003	0.002	1736.0	1691.0	1.00
v_1|participant_id[9]	-0.241	0.122	-0.475	-0.020	0.003	0.002	1864.0	1534.0	1.00
a_Intercept	1.165	0.016	1.136	1.196	0.000	0.000	1591.0	1595.0	1.00

In [108]:

az.summary(idata_B, var_names=["~offset", "~participant_id"], filter_vars="regex")

az.summary(idata_B, var_names=["~offset", "~participant_id"], filter_vars="regex")

Out[108]:

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
v_age	0.430	0.129	0.187	0.669	0.005	0.004	707.0	857.0	1.0
v_difficulty	0.758	0.132	0.496	0.987	0.005	0.004	606.0	864.0	1.0
a_Intercept	1.224	0.016	1.193	1.252	0.000	0.000	1778.0	1331.0	1.0
z_Intercept	0.509	0.007	0.494	0.522	0.000	0.000	1986.0	1414.0	1.0
t_Intercept	0.258	0.008	0.244	0.271	0.000	0.000	1671.0	1349.0	1.0

In [168]:

az.plot_trace(idata_B, var_names=["~offset", "~participant_id"], filter_vars="regex")
plt.tight_layout()

az.plot_trace(idata_B, var_names=["~offset", "~participant_id"], filter_vars="regex") plt.tight_layout()

/Users/mayan/HSSM/.venv/lib/python3.12/site-packages/arviz/utils.py:146: UserWarning: Items starting with ~: ['offset'] have not been found and will be ignored
  warnings.warn(

Step 4: Posterior predictive check¶

In [169]:

g = hssm.plotting.plot_predictive(model_B, col="participant_id", col_wrap=5)
g.set(xlim=(-7, 7))

g = hssm.plotting.plot_predictive(model_B, col="participant_id", col_wrap=5) g.set(xlim=(-7, 7))

Out[169]:

<seaborn.axisgrid.FacetGrid at 0x13fd224b0>

In [171]:

# Participant-wise parameter recovery for v (random intercepts and slopes)
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 5), sharey=True)

for ax, (param_name, truth, label) in zip(
    axes,
    [
        ("v_1|participant_id", v_int_p, "v random intercept"),
        ("v_difficulty|participant_id", v_slope_p, "v random slope on difficulty"),
    ],
):
    post = idata_B.posterior[param_name]
    mean = post.mean(dim=["chain", "draw"]).values
    hdi = az.hdi(post, hdi_prob=0.95)[param_name].values

    yvals = np.arange(n_participants)
    ax.hlines(yvals, hdi[:, 0], hdi[:, 1], color="blue", alpha=0.5)
    ax.plot(mean, yvals, "|", color="blue", label="posterior mean")
    ax.plot(truth, yvals, "x", color="red", label="ground truth")
    ax.set_title(label)
    ax.set_yticks(yvals)
    ax.set_ylabel("participant_id")
    ax.legend()

plt.tight_layout()

# Participant-wise parameter recovery for v (random intercepts and slopes) fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 5), sharey=True) for ax, (param_name, truth, label) in zip( axes, [ ("v_1|participant_id", v_int_p, "v random intercept"), ("v_difficulty|participant_id", v_slope_p, "v random slope on difficulty"), ], ): post = idata_B.posterior[param_name] mean = post.mean(dim=["chain", "draw"]).values hdi = az.hdi(post, hdi_prob=0.95)[param_name].values yvals = np.arange(n_participants) ax.hlines(yvals, hdi[:, 0], hdi[:, 1], color="blue", alpha=0.5) ax.plot(mean, yvals, "|", color="blue", label="posterior mean") ax.plot(truth, yvals, "x", color="red", label="ground truth") ax.set_title(label) ax.set_yticks(yvals) ax.set_ylabel("participant_id") ax.legend() plt.tight_layout()

In [170]:

# Posterior predictive: split by age level (between-subject factor)
g = hssm.plotting.plot_predictive(model_B, col="age_level", col_wrap=3)
g.set(xlim=(-7, 7))
g.fig.suptitle("")

# Posterior predictive: split by age level (between-subject factor) g = hssm.plotting.plot_predictive(model_B, col="age_level", col_wrap=3) g.set(xlim=(-7, 7)) g.fig.suptitle("")

Out[170]:

Text(0.5, 0.98, '')

In [172]:

# Posterior predictive: split by difficulty level (within-subject factor)
g = hssm.plotting.plot_predictive(model_B, col="difficulty_level", col_wrap=3)
g.set(xlim=(-7, 7))
g.fig.suptitle("")

# Posterior predictive: split by difficulty level (within-subject factor) g = hssm.plotting.plot_predictive(model_B, col="difficulty_level", col_wrap=3) g.set(xlim=(-7, 7)) g.fig.suptitle("")

Out[172]:

Text(0.5, 0.98, '')

Variant C: interaction¶

What if older and younger participants don't just differ in their baseline drift rate, but also in how much difficulty affects them? That's an interaction between age and difficulty.

To capture this, we add an age:difficulty term to the v regression. Same setup as Variant B otherwise.

Step 1: Simulate the data¶

In [155]:

# true coefficients for variant C
beta_v_difficulty = 0.08
beta_v_age = 0.005
beta_v_age_difficulty = 0.0005

v_int_p = rng.normal(0, sigma_v_int, n_participants)
v_slope_p = rng.normal(0, sigma_v_slope, n_participants)

v_int_offset = v_int_p / sigma_v_int
v_slope_offset = v_slope_p / sigma_v_slope

# Build HSSM model with the regression formula
dummy_data = covariates.copy()
dummy_data["rt"] = 1.0
dummy_data["response"] = 1.0

model_C_gen = hssm.HSSM(
    data=dummy_data,
    model="ddm",
    global_formula="y ~ 1",
    include=[
        {
            "name": "v",
            "formula": (
                "v ~ 0 + age + difficulty + age:difficulty + "
                "(1 + difficulty | participant_id)"
            ),
        }
    ],
    noncentered=True,
)

# Use do-operator: intervene on regression coefficients and random effects
synth_idata_C, synth_model_C = model_C_gen.sample_do(
    params={
        "v_age": beta_v_age,
        "v_difficulty": beta_v_difficulty,
        "v_age:difficulty": beta_v_age_difficulty,
        "v_1|participant_id_mu": 0.0,
        "v_1|participant_id_sigma": sigma_v_int,
        "v_1|participant_id_offset": v_int_offset,
        "v_difficulty|participant_id_mu": 0.0,
        "v_difficulty|participant_id_sigma": sigma_v_slope,
        "v_difficulty|participant_id_offset": v_slope_offset,
        "a_Intercept": a_true,
        "z_Intercept": z_true,
        "t_Intercept": t_true,
    },
    draws=1,
    var_names=["rt,response"],
    return_model=True,
)

pm.model_to_graphviz(synth_model_C, var_names=["rt,response"])

# true coefficients for variant C beta_v_difficulty = 0.08 beta_v_age = 0.005 beta_v_age_difficulty = 0.0005 v_int_p = rng.normal(0, sigma_v_int, n_participants) v_slope_p = rng.normal(0, sigma_v_slope, n_participants) v_int_offset = v_int_p / sigma_v_int v_slope_offset = v_slope_p / sigma_v_slope # Build HSSM model with the regression formula dummy_data = covariates.copy() dummy_data["rt"] = 1.0 dummy_data["response"] = 1.0 model_C_gen = hssm.HSSM( data=dummy_data, model="ddm", global_formula="y ~ 1", include=[ { "name": "v", "formula": ( "v ~ 0 + age + difficulty + age:difficulty + " "(1 + difficulty | participant_id)" ), } ], noncentered=True, ) # Use do-operator: intervene on regression coefficients and random effects synth_idata_C, synth_model_C = model_C_gen.sample_do( params={ "v_age": beta_v_age, "v_difficulty": beta_v_difficulty, "v_age:difficulty": beta_v_age_difficulty, "v_1|participant_id_mu": 0.0, "v_1|participant_id_sigma": sigma_v_int, "v_1|participant_id_offset": v_int_offset, "v_difficulty|participant_id_mu": 0.0, "v_difficulty|participant_id_sigma": sigma_v_slope, "v_difficulty|participant_id_offset": v_slope_offset, "a_Intercept": a_true, "z_Intercept": z_true, "t_Intercept": t_true, }, draws=1, var_names=["rt,response"], return_model=True, ) pm.model_to_graphviz(synth_model_C, var_names=["rt,response"])

No common intercept. Bounds for parameter v is not applied due to a current limitation of Bambi. This will change in the future.

/Users/mayan/HSSM/src/hssm/hssm.py:504: UserWarning: You set choices to be (-1, 1), but [-1] are missing from your dataset.
  self._post_check_data_sanity()

Model initialized successfully.

Sampling: [rt,response]

Out[155]:

In [156]:

synth_df_C = hssm.utils.predictive_idata_to_dataframe(
    synth_idata_C, predictive_group="prior_predictive"
)

data_C = covariates.copy()
data_C["rt"] = synth_df_C["rt"].values
data_C["response"] = synth_df_C["response"].values
data_C.to_csv("sim_variant_C.csv", index=False)

synth_df_C = hssm.utils.predictive_idata_to_dataframe( synth_idata_C, predictive_group="prior_predictive" ) data_C = covariates.copy() data_C["rt"] = synth_df_C["rt"].values data_C["response"] = synth_df_C["response"].values data_C.to_csv("sim_variant_C.csv", index=False)

Step 2: Fit a fresh HSSM model¶

In [157]:

data_C = pd.read_csv("sim_variant_C.csv")

# Bin difficulty (3 levels) and age (3 levels) for faceted PPC plots
data_C["difficulty_level"] = pd.qcut(
    data_C["difficulty"], q=3, labels=["low", "med", "high"]
)
data_C["age_level"] = pd.qcut(data_C["age"], q=3, labels=["young", "middle", "older"])

model_C = hssm.HSSM(
    data=data_C,
    model="ddm",
    global_formula="y ~ 1",
    include=[
        {
            "name": "v",
            "formula": (
                "v ~ 0 + age + difficulty + age:difficulty + "
                "(1 + difficulty | participant_id)"
            ),
        }
    ],
    noncentered=True,
)

data_C = pd.read_csv("sim_variant_C.csv") # Bin difficulty (3 levels) and age (3 levels) for faceted PPC plots data_C["difficulty_level"] = pd.qcut( data_C["difficulty"], q=3, labels=["low", "med", "high"] ) data_C["age_level"] = pd.qcut(data_C["age"], q=3, labels=["young", "middle", "older"]) model_C = hssm.HSSM( data=data_C, model="ddm", global_formula="y ~ 1", include=[ { "name": "v", "formula": ( "v ~ 0 + age + difficulty + age:difficulty + " "(1 + difficulty | participant_id)" ), } ], noncentered=True, )

No common intercept. Bounds for parameter v is not applied due to a current limitation of Bambi. This will change in the future.
Model initialized successfully.

In [158]:

print(model_C)
pm.model_to_graphviz(model_C.pymc_model, var_names=["rt,response"])

print(model_C) pm.model_to_graphviz(model_C.pymc_model, var_names=["rt,response"])

Hierarchical Sequential Sampling Model
Model: ddm

Response variable: rt,response
Likelihood: analytical
Observations: 2500

Parameters:

v:
    Formula: v ~ 0 + age + difficulty + age:difficulty + (1 + difficulty | participant_id)
    Priors:
        v_age ~ Normal(mu: 0.0, sigma: 0.25)
        v_difficulty ~ Normal(mu: 0.0, sigma: 0.25)
        v_age:difficulty ~ Normal(mu: 0.0, sigma: 0.25)
        v_1|participant_id ~ Normal(mu: Normal(mu: 2.0, sigma: 3.0), sigma: HalfNormal(sigma: 2.0))
        v_difficulty|participant_id ~ Normal(mu: Normal(mu: 0.0, sigma: 0.25), sigma: Weibull(alpha: 1.5, beta: 0.3))
    Link: identity
    Explicit bounds: (-inf, inf)

a:
    Formula: a ~ 1
    Priors:
        a_Intercept ~ Gamma(mu: 1.5, sigma: 0.75)
    Link: identity
    Explicit bounds: (0.0, inf)

z:
    Formula: z ~ 1
    Priors:
        z_Intercept ~ Beta(alpha: 10.0, beta: 10.0)
    Link: identity
    Explicit bounds: (0.0, 1.0)

t:
    Formula: t ~ 1
    Priors:
        t_Intercept ~ Gamma(mu: 0.2, sigma: 0.2)
    Link: identity
    Explicit bounds: (0.0, inf)


Lapse probability: 0.05
Lapse distribution: Uniform(lower: 0.0, upper: 20.0)

Out[158]:

In [159]:

idata_C = model_C.sample(**sample_kwargs)

idata_C = model_C.sample(**sample_kwargs)

Using default initvals. 

Initializing NUTS using adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [v_age, v_difficulty, v_age:difficulty, v_1|participant_id_mu, v_1|participant_id_sigma, v_1|participant_id_offset, v_difficulty|participant_id_mu, v_difficulty|participant_id_sigma, v_difficulty|participant_id_offset, a_Intercept, z_Intercept, t_Intercept]
/Users/mayan/.local/share/uv/python/cpython-3.12.12-macos-aarch64-none/lib/python3.12/multiprocessing/popen_fork.py:66: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.
  self.pid = os.fork()

Output()

/Users/mayan/.local/share/uv/python/cpython-3.12.12-macos-aarch64-none/lib/python3.12/multiprocessing/popen_fork.py:66: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.
  self.pid = os.fork()

Sampling 2 chains for 500 tune and 1_000 draw iterations (1_000 + 2_000 draws total) took 856 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
100%|██████████| 2000/2000 [00:01<00:00, 1386.90it/s]

Step 3: Inspect what got estimated¶

In [160]:

az.summary(idata_C)

az.summary(idata_C)

Out[160]:

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
v_difficulty|participant_id_mu	0.007	0.249	-0.464	0.440	0.005	0.006	2583.0	1448.0	1.00
v_age	0.004	0.003	-0.001	0.010	0.000	0.000	872.0	813.0	1.00
v_difficulty|participant_id[0]	-0.101	0.090	-0.271	0.068	0.004	0.002	621.0	1132.0	1.00
v_difficulty|participant_id[1]	-0.166	0.080	-0.318	-0.012	0.003	0.002	745.0	855.0	1.01
v_difficulty|participant_id[2]	-0.129	0.116	-0.359	0.077	0.005	0.003	604.0	980.0	1.00
v_difficulty|participant_id[3]	-0.120	0.071	-0.258	0.016	0.003	0.002	692.0	1052.0	1.00
v_difficulty|participant_id[4]	0.117	0.082	-0.035	0.280	0.003	0.002	693.0	1042.0	1.00
v_difficulty|participant_id[5]	0.034	0.091	-0.141	0.204	0.003	0.003	729.0	1004.0	1.00
v_difficulty|participant_id[6]	0.366	0.081	0.208	0.515	0.003	0.002	780.0	889.0	1.00
v_difficulty|participant_id[7]	0.176	0.100	-0.018	0.357	0.004	0.003	607.0	1008.0	1.00
v_difficulty|participant_id[8]	-0.062	0.086	-0.220	0.105	0.003	0.003	739.0	938.0	1.00
v_difficulty|participant_id[9]	-0.077	0.088	-0.234	0.099	0.003	0.003	730.0	889.0	1.00
v_1|participant_id_mu	1.995	2.857	-3.795	7.025	0.058	0.063	2466.0	1450.0	1.00
v_1|participant_id_sigma	0.379	0.125	0.177	0.610	0.004	0.004	761.0	1005.0	1.00
z_Intercept	0.505	0.007	0.491	0.518	0.000	0.000	1917.0	1196.0	1.00
v_difficulty|participant_id_sigma	0.204	0.057	0.116	0.316	0.002	0.002	626.0	931.0	1.00
v_age:difficulty	0.000	0.003	-0.005	0.006	0.000	0.000	627.0	807.0	1.00
t_Intercept	0.260	0.007	0.247	0.274	0.000	0.000	1800.0	1261.0	1.00
v_difficulty	0.059	0.141	-0.221	0.319	0.005	0.004	689.0	788.0	1.00
v_1|participant_id[0]	0.015	0.190	-0.355	0.383	0.006	0.005	1036.0	1202.0	1.00
v_1|participant_id[1]	0.741	0.146	0.477	1.010	0.004	0.003	1679.0	1889.0	1.00
v_1|participant_id[2]	-0.223	0.216	-0.628	0.186	0.007	0.006	1045.0	1076.0	1.00
v_1|participant_id[3]	0.104	0.174	-0.209	0.443	0.005	0.005	1265.0	1139.0	1.00
v_1|participant_id[4]	0.170	0.189	-0.155	0.553	0.005	0.005	1238.0	1270.0	1.00
v_1|participant_id[5]	0.069	0.118	-0.151	0.290	0.003	0.002	1757.0	1504.0	1.00
v_1|participant_id[6]	0.113	0.137	-0.144	0.383	0.003	0.003	1581.0	1437.0	1.00
v_1|participant_id[7]	-0.376	0.197	-0.752	-0.012	0.006	0.006	969.0	905.0	1.00
v_1|participant_id[8]	0.066	0.122	-0.179	0.287	0.003	0.003	1479.0	1371.0	1.00
v_1|participant_id[9]	-0.080	0.121	-0.310	0.133	0.003	0.002	1588.0	1506.0	1.00
a_Intercept	1.195	0.014	1.170	1.222	0.000	0.000	1781.0	1341.0	1.00

In [161]:

az.summary(idata_C, var_names=["~offset", "~participant_id"], filter_vars="regex")

az.summary(idata_C, var_names=["~offset", "~participant_id"], filter_vars="regex")

/Users/mayan/HSSM/.venv/lib/python3.12/site-packages/arviz/utils.py:146: UserWarning: Items starting with ~: ['offset'] have not been found and will be ignored
  warnings.warn(

Out[161]:

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
v_age	0.004	0.003	-0.001	0.010	0.000	0.000	872.0	813.0	1.0
z_Intercept	0.505	0.007	0.491	0.518	0.000	0.000	1917.0	1196.0	1.0
v_age:difficulty	0.000	0.003	-0.005	0.006	0.000	0.000	627.0	807.0	1.0
t_Intercept	0.260	0.007	0.247	0.274	0.000	0.000	1800.0	1261.0	1.0
v_difficulty	0.059	0.141	-0.221	0.319	0.005	0.004	689.0	788.0	1.0
a_Intercept	1.195	0.014	1.170	1.222	0.000	0.000	1781.0	1341.0	1.0

In [162]:

az.plot_trace(idata_C, var_names=["~offset", "~participant_id"], filter_vars="regex")
plt.tight_layout()

az.plot_trace(idata_C, var_names=["~offset", "~participant_id"], filter_vars="regex") plt.tight_layout()

/Users/mayan/HSSM/.venv/lib/python3.12/site-packages/arviz/utils.py:146: UserWarning: Items starting with ~: ['offset'] have not been found and will be ignored
  warnings.warn(

Step 4: Posterior predictive check¶

In [163]:

g = hssm.plotting.plot_predictive(model_C, col="participant_id", col_wrap=5)
g.set(xlim=(-7, 7))

g = hssm.plotting.plot_predictive(model_C, col="participant_id", col_wrap=5) g.set(xlim=(-7, 7))

No posterior_predictive samples found. Generating posterior_predictive samples using the provided InferenceData object and the original data. This will modify the provided InferenceData object, or if not provided, the traces object stored inside the model.

In [164]:

# Participant-wise parameter recovery for v (random intercepts and slopes)
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 5), sharey=True)

for ax, (param_name, truth, label) in zip(
    axes,
    [
        ("v_1|participant_id", v_int_p, "v random intercept"),
        ("v_difficulty|participant_id", v_slope_p, "v random slope on difficulty"),
    ],
):
    post = idata_C.posterior[param_name]
    mean = post.mean(dim=["chain", "draw"]).values
    hdi = az.hdi(post, hdi_prob=0.95)[param_name].values

    yvals = np.arange(n_participants)
    ax.hlines(yvals, hdi[:, 0], hdi[:, 1], color="blue", alpha=0.5)
    ax.plot(mean, yvals, "|", color="blue", label="posterior mean")
    ax.plot(truth, yvals, "x", color="red", label="ground truth")
    ax.set_title(label)
    ax.set_yticks(yvals)
    ax.set_ylabel("participant_id")
    ax.legend()

plt.tight_layout()

# Participant-wise parameter recovery for v (random intercepts and slopes) fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 5), sharey=True) for ax, (param_name, truth, label) in zip( axes, [ ("v_1|participant_id", v_int_p, "v random intercept"), ("v_difficulty|participant_id", v_slope_p, "v random slope on difficulty"), ], ): post = idata_C.posterior[param_name] mean = post.mean(dim=["chain", "draw"]).values hdi = az.hdi(post, hdi_prob=0.95)[param_name].values yvals = np.arange(n_participants) ax.hlines(yvals, hdi[:, 0], hdi[:, 1], color="blue", alpha=0.5) ax.plot(mean, yvals, "|", color="blue", label="posterior mean") ax.plot(truth, yvals, "x", color="red", label="ground truth") ax.set_title(label) ax.set_yticks(yvals) ax.set_ylabel("participant_id") ax.legend() plt.tight_layout()

In [165]:

# Posterior predictive: faceted by age level (between-subject factor)
g = hssm.plotting.plot_predictive(model_C, col="age_level", col_wrap=3)
g.set(xlim=(-7, 7))
g.fig.suptitle("")

# Posterior predictive: faceted by age level (between-subject factor) g = hssm.plotting.plot_predictive(model_C, col="age_level", col_wrap=3) g.set(xlim=(-7, 7)) g.fig.suptitle("")

Out[165]:

Text(0.5, 0.98, '')

In [166]:

# Posterior predictive: faceted by difficulty level (within-subject factor)
g = hssm.plotting.plot_predictive(model_C, col="difficulty_level", col_wrap=3)
g.set(xlim=(-7, 7))
g.fig.suptitle("")

# Posterior predictive: faceted by difficulty level (within-subject factor) g = hssm.plotting.plot_predictive(model_C, col="difficulty_level", col_wrap=3) g.set(xlim=(-7, 7)) g.fig.suptitle("")

Out[166]:

Text(0.5, 0.98, '')

In [167]:

# Posterior predictive: 3x3 interaction grid (age x difficulty)
g = hssm.plotting.plot_predictive(model_C, row="age_level", col="difficulty_level")
g.set(xlim=(-7, 7))
g.set_titles(template="age={row_name} | difficulty={col_name}")
g.fig.suptitle("")

# Posterior predictive: 3x3 interaction grid (age x difficulty) g = hssm.plotting.plot_predictive(model_C, row="age_level", col="difficulty_level") g.set(xlim=(-7, 7)) g.set_titles(template="age={row_name} | difficulty={col_name}") g.fig.suptitle("")

Out[167]:

Text(0.5, 0.98, '')

Variant D: full hierarchy on every parameter¶

In Variants A–C we kept a, z, t as global scalars to keep the spotlight on the v regression. In real experiments, however, every DDM parameter usually shows individual differences — some people have higher boundary separations, some have longer non-decision times, etc.

Variant D shows what the fully hierarchical case looks like: same v regression as Variant A, but now global_formula="y ~ 1 + (1 | participant_id)" adds participant-level random intercepts to a, z, t as well. The simulator generates per-participant values for these to match.

This is the most realistic specification for real data, but it's also slower to fit (more parameters to sample).

Step 1: Simulate the data¶

In [173]:

# Variant D: simulate with per-participant a, z, t
beta_v_difficulty = 0.08

v_int_p = rng.normal(0, sigma_v_int, n_participants)
v_slope_p = rng.normal(0, sigma_v_slope, n_participants)

# Per-participant a, z, t (centered on the scalar means used in A/B/C)
sigma_a, sigma_z, sigma_t = 0.15, 0.05, 0.05
a_p = a_true + rng.normal(0, sigma_a, n_participants)
z_p = np.clip(z_true + rng.normal(0, sigma_z, n_participants), 0.05, 0.95)
# t clipped to >= 0.001 (relaxed from the old 0.05 bound)
t_p = np.clip(t_true + rng.normal(0, sigma_t, n_participants), 0.001, None)

# Noncentered offsets, derived from the (possibly clipped) per-participant values
v_int_offset = v_int_p / sigma_v_int
v_slope_offset = v_slope_p / sigma_v_slope
a_int_offset = (a_p - a_true) / sigma_a
z_int_offset = (z_p - z_true) / sigma_z
t_int_offset = (t_p - t_true) / sigma_t

# Build HSSM model with full hierarchy
dummy_data = covariates.copy()
dummy_data["rt"] = 1.0
dummy_data["response"] = 1.0

model_D_gen = hssm.HSSM(
    data=dummy_data,
    model="ddm",
    global_formula="y ~ 1 + (1 | participant_id)",
    include=[
        {
            "name": "v",
            "formula": "v ~ 0 + difficulty + (1 + difficulty | participant_id)",
        }
    ],
    noncentered=True,
)

# Use do-operator: intervene on every regression coefficient and random effect
synth_idata_D, synth_model_D = model_D_gen.sample_do(
    params={
        "v_difficulty": beta_v_difficulty,
        "v_1|participant_id_mu": 0.0,
        "v_1|participant_id_sigma": sigma_v_int,
        "v_1|participant_id_offset": v_int_offset,
        "v_difficulty|participant_id_mu": 0.0,
        "v_difficulty|participant_id_sigma": sigma_v_slope,
        "v_difficulty|participant_id_offset": v_slope_offset,
        "a_Intercept": a_true,
        "a_1|participant_id_sigma": sigma_a,
        "a_1|participant_id_offset": a_int_offset,
        "z_Intercept": z_true,
        "z_1|participant_id_sigma": sigma_z,
        "z_1|participant_id_offset": z_int_offset,
        "t_Intercept": t_true,
        "t_1|participant_id_sigma": sigma_t,
        "t_1|participant_id_offset": t_int_offset,
    },
    draws=1,
    var_names=["rt,response"],
    return_model=True,
)

pm.model_to_graphviz(synth_model_D, var_names=["rt,response"])

# Variant D: simulate with per-participant a, z, t beta_v_difficulty = 0.08 v_int_p = rng.normal(0, sigma_v_int, n_participants) v_slope_p = rng.normal(0, sigma_v_slope, n_participants) # Per-participant a, z, t (centered on the scalar means used in A/B/C) sigma_a, sigma_z, sigma_t = 0.15, 0.05, 0.05 a_p = a_true + rng.normal(0, sigma_a, n_participants) z_p = np.clip(z_true + rng.normal(0, sigma_z, n_participants), 0.05, 0.95) # t clipped to >= 0.001 (relaxed from the old 0.05 bound) t_p = np.clip(t_true + rng.normal(0, sigma_t, n_participants), 0.001, None) # Noncentered offsets, derived from the (possibly clipped) per-participant values v_int_offset = v_int_p / sigma_v_int v_slope_offset = v_slope_p / sigma_v_slope a_int_offset = (a_p - a_true) / sigma_a z_int_offset = (z_p - z_true) / sigma_z t_int_offset = (t_p - t_true) / sigma_t # Build HSSM model with full hierarchy dummy_data = covariates.copy() dummy_data["rt"] = 1.0 dummy_data["response"] = 1.0 model_D_gen = hssm.HSSM( data=dummy_data, model="ddm", global_formula="y ~ 1 + (1 | participant_id)", include=[ { "name": "v", "formula": "v ~ 0 + difficulty + (1 + difficulty | participant_id)", } ], noncentered=True, ) # Use do-operator: intervene on every regression coefficient and random effect synth_idata_D, synth_model_D = model_D_gen.sample_do( params={ "v_difficulty": beta_v_difficulty, "v_1|participant_id_mu": 0.0, "v_1|participant_id_sigma": sigma_v_int, "v_1|participant_id_offset": v_int_offset, "v_difficulty|participant_id_mu": 0.0, "v_difficulty|participant_id_sigma": sigma_v_slope, "v_difficulty|participant_id_offset": v_slope_offset, "a_Intercept": a_true, "a_1|participant_id_sigma": sigma_a, "a_1|participant_id_offset": a_int_offset, "z_Intercept": z_true, "z_1|participant_id_sigma": sigma_z, "z_1|participant_id_offset": z_int_offset, "t_Intercept": t_true, "t_1|participant_id_sigma": sigma_t, "t_1|participant_id_offset": t_int_offset, }, draws=1, var_names=["rt,response"], return_model=True, ) pm.model_to_graphviz(synth_model_D, var_names=["rt,response"])

No common intercept. Bounds for parameter v is not applied due to a current limitation of Bambi. This will change in the future.
Model initialized successfully.

/Users/mayan/HSSM/src/hssm/hssm.py:504: UserWarning: You set choices to be (-1, 1), but [-1] are missing from your dataset.
  self._post_check_data_sanity()
Sampling: [rt,response]

Out[173]:

In [174]:

synth_df_D = hssm.utils.predictive_idata_to_dataframe(
    synth_idata_D, predictive_group="prior_predictive"
)

data_D = covariates.copy()
data_D["rt"] = synth_df_D["rt"].values
data_D["response"] = synth_df_D["response"].values
data_D.to_csv("sim_variant_D.csv", index=False)
data_D.head()

synth_df_D = hssm.utils.predictive_idata_to_dataframe( synth_idata_D, predictive_group="prior_predictive" ) data_D = covariates.copy() data_D["rt"] = synth_df_D["rt"].values data_D["response"] = synth_df_D["response"].values data_D.to_csv("sim_variant_D.csv", index=False) data_D.head()

Out[174]:

	participant_id	difficulty	age	rt	response
0	0	6.369617	60.357299	1.303465	1.0
1	0	2.697867	60.357299	0.989182	1.0
2	0	0.409735	60.357299	1.040517	-1.0
3	0	0.165276	60.357299	2.973009	1.0
4	0	8.132702	60.357299	1.392237	1.0

Step 2: Fit a fresh HSSM model¶

In [175]:

data_D = pd.read_csv("sim_variant_D.csv")

model_D = hssm.HSSM(
    data=data_D,
    model="ddm",
    global_formula="y ~ 1 + (1 | participant_id)",
    include=[
        {
            "name": "v",
            "formula": "v ~ 0 + difficulty + (1 + difficulty | participant_id)",
        }
    ],
    noncentered=True,
)

data_D = pd.read_csv("sim_variant_D.csv") model_D = hssm.HSSM( data=data_D, model="ddm", global_formula="y ~ 1 + (1 | participant_id)", include=[ { "name": "v", "formula": "v ~ 0 + difficulty + (1 + difficulty | participant_id)", } ], noncentered=True, )

No common intercept. Bounds for parameter v is not applied due to a current limitation of Bambi. This will change in the future.
Model initialized successfully.

In [ ]:

print(model_D)
pm.model_to_graphviz(model_D.pymc_model, var_names=["rt,response"])

print(model_D) pm.model_to_graphviz(model_D.pymc_model, var_names=["rt,response"])

Hierarchical Sequential Sampling Model
Model: ddm

Response variable: rt,response
Likelihood: analytical
Observations: 2500

Parameters:

v:
    Formula: v ~ 0 + x + (1 + x | participant_id)
    Priors:
        v_x ~ Normal(mu: 0.0, sigma: 0.25)
        v_1|participant_id ~ Normal(mu: Normal(mu: 2.0, sigma: 3.0), sigma: HalfNormal(sigma: 2.0))
        v_x|participant_id ~ Normal(mu: Normal(mu: 0.0, sigma: 0.25), sigma: Weibull(alpha: 1.5, beta: 0.3))
    Link: identity
    Explicit bounds: (-inf, inf)

a:
    Formula: a ~ 1 + (1 | participant_id)
    Priors:
        a_Intercept ~ Gamma(mu: 1.5, sigma: 0.75)
        a_1|participant_id ~ Normal(mu: 0.0, sigma: Weibull(alpha: 1.5, beta: 0.3))
    Link: identity
    Explicit bounds: (0.0, inf)

z:
    Formula: z ~ 1 + (1 | participant_id)
    Priors:
        z_Intercept ~ Beta(alpha: 10.0, beta: 10.0)
        z_1|participant_id ~ Normal(mu: 0.0, sigma: Weibull(alpha: 1.5, beta: 0.3))
    Link: identity
    Explicit bounds: (0.0, 1.0)

t:
    Formula: t ~ 1 + (1 | participant_id)
    Priors:
        t_Intercept ~ Gamma(mu: 0.2, sigma: 0.2)
        t_1|participant_id ~ Normal(mu: 0.0, sigma: Weibull(alpha: 1.5, beta: 0.3))
    Link: identity
    Explicit bounds: (0.0, inf)


Lapse probability: 0.05
Lapse distribution: Uniform(lower: 0.0, upper: 20.0)

Out[ ]:

In [176]:

idata_D = model_D.sample(**sample_kwargs)

idata_D = model_D.sample(**sample_kwargs)

Using default initvals. 

Initializing NUTS using adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [v_difficulty, v_1|participant_id_mu, v_1|participant_id_sigma, v_1|participant_id_offset, v_difficulty|participant_id_mu, v_difficulty|participant_id_sigma, v_difficulty|participant_id_offset, a_Intercept, a_1|participant_id_sigma, a_1|participant_id_offset, z_Intercept, z_1|participant_id_sigma, z_1|participant_id_offset, t_Intercept, t_1|participant_id_sigma, t_1|participant_id_offset]
/Users/mayan/.local/share/uv/python/cpython-3.12.12-macos-aarch64-none/lib/python3.12/multiprocessing/popen_fork.py:66: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.
  self.pid = os.fork()

Output()

/Users/mayan/.local/share/uv/python/cpython-3.12.12-macos-aarch64-none/lib/python3.12/multiprocessing/popen_fork.py:66: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.
  self.pid = os.fork()

Sampling 2 chains for 500 tune and 1_000 draw iterations (1_000 + 2_000 draws total) took 415 seconds.
We recommend running at least 4 chains for robust computation of convergence diagnostics
100%|██████████| 2000/2000 [00:01<00:00, 1349.72it/s]

Step 3: Inspect what got estimated¶

In [177]:

az.summary(idata_D)

az.summary(idata_D)

Out[177]:

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
v_difficulty|participant_id_mu	-0.002	0.255	-0.444	0.498	0.005	0.006	2760.0	1335.0	1.0
a_1|participant_id_sigma	0.162	0.051	0.080	0.252	0.002	0.002	658.0	929.0	1.0
v_difficulty|participant_id[0]	0.043	0.078	-0.108	0.188	0.004	0.003	347.0	427.0	1.0
v_difficulty|participant_id[1]	0.238	0.080	0.083	0.384	0.004	0.003	372.0	455.0	1.0
v_difficulty|participant_id[2]	-0.091	0.078	-0.238	0.062	0.004	0.003	332.0	464.0	1.0
...	...	...	...	...	...	...	...	...	...
t_1|participant_id[6]	0.057	0.024	0.012	0.103	0.001	0.001	1314.0	1464.0	1.0
t_1|participant_id[7]	0.037	0.026	-0.011	0.087	0.001	0.001	1450.0	1253.0	1.0
t_1|participant_id[8]	-0.047	0.022	-0.091	-0.005	0.001	0.001	1153.0	1317.0	1.0
t_1|participant_id[9]	-0.014	0.025	-0.059	0.034	0.001	0.001	1206.0	967.0	1.0
a_Intercept	1.211	0.056	1.114	1.323	0.002	0.002	700.0	1033.0	1.0

61 rows × 9 columns

In [178]:

az.summary(idata_D, var_names=["~offset", "~participant_id"], filter_vars="regex")

az.summary(idata_D, var_names=["~offset", "~participant_id"], filter_vars="regex")

/Users/mayan/HSSM/.venv/lib/python3.12/site-packages/arviz/utils.py:146: UserWarning: Items starting with ~: ['offset'] have not been found and will be ignored
  warnings.warn(

Out[178]:

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
z_Intercept	0.491	0.023	0.451	0.539	0.001	0.001	950.0	1293.0	1.0
t_Intercept	0.265	0.017	0.234	0.299	0.001	0.001	974.0	1035.0	1.0
v_difficulty	0.091	0.076	-0.052	0.242	0.004	0.003	323.0	436.0	1.0
a_Intercept	1.211	0.056	1.114	1.323	0.002	0.002	700.0	1033.0	1.0

In [179]:

az.plot_trace(idata_D, var_names=["~offset", "~participant_id"], filter_vars="regex")
plt.tight_layout()

az.plot_trace(idata_D, var_names=["~offset", "~participant_id"], filter_vars="regex") plt.tight_layout()

/Users/mayan/HSSM/.venv/lib/python3.12/site-packages/arviz/utils.py:146: UserWarning: Items starting with ~: ['offset'] have not been found and will be ignored
  warnings.warn(

Step 4: Posterior predictive check¶

In [181]:

g = hssm.plotting.plot_predictive(model_D, col="participant_id", col_wrap=5)
g.set(xlim=(-7, 7))

g = hssm.plotting.plot_predictive(model_D, col="participant_id", col_wrap=5) g.set(xlim=(-7, 7))

No posterior_predictive samples found. Generating posterior_predictive samples using the provided InferenceData object and the original data. This will modify the provided InferenceData object, or if not provided, the traces object stored inside the model.

Out[181]:

<seaborn.axisgrid.FacetGrid at 0x138eb5730>

In [180]:

# Participant-wise parameter recovery for all DDM parameters
fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(16, 5), sharey=True)

for ax, (param_name, truth, label) in zip(
    axes,
    [
        ("v_1|participant_id", v_int_p, "v (random intercept)"),
        ("a_1|participant_id", a_p - a_true, "a (random intercept)"),
        ("z_1|participant_id", z_p - z_true, "z (random intercept)"),
        ("t_1|participant_id", t_p - t_true, "t (random intercept)"),
    ],
):
    post = idata_D.posterior[param_name]
    mean = post.mean(dim=["chain", "draw"]).values
    hdi = az.hdi(post, hdi_prob=0.95)[param_name].values

    yvals = np.arange(n_participants)
    ax.hlines(yvals, hdi[:, 0], hdi[:, 1], color="blue", alpha=0.5)
    ax.plot(mean, yvals, "|", color="blue", label="posterior mean")
    ax.plot(truth, yvals, "x", color="red", label="ground truth")
    ax.set_title(label)
    ax.set_yticks(yvals)
    if ax is axes[0]:
        ax.set_ylabel("participant_id")
        ax.legend()

plt.tight_layout()

# Participant-wise parameter recovery for all DDM parameters fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(16, 5), sharey=True) for ax, (param_name, truth, label) in zip( axes, [ ("v_1|participant_id", v_int_p, "v (random intercept)"), ("a_1|participant_id", a_p - a_true, "a (random intercept)"), ("z_1|participant_id", z_p - z_true, "z (random intercept)"), ("t_1|participant_id", t_p - t_true, "t (random intercept)"), ], ): post = idata_D.posterior[param_name] mean = post.mean(dim=["chain", "draw"]).values hdi = az.hdi(post, hdi_prob=0.95)[param_name].values yvals = np.arange(n_participants) ax.hlines(yvals, hdi[:, 0], hdi[:, 1], color="blue", alpha=0.5) ax.plot(mean, yvals, "|", color="blue", label="posterior mean") ax.plot(truth, yvals, "x", color="red", label="ground truth") ax.set_title(label) ax.set_yticks(yvals) if ax is axes[0]: ax.set_ylabel("participant_id") ax.legend() plt.tight_layout()