Best approach to derive Bayesian p-values in HSSM

[cuddenberg]

Hi all!

I'm new to HSSM and using it to model data from a two forced choice task (data from the experiment is published in papers such as this, building off of the project mentioned in previous post #473. Similar to #421, #287 and #314, I'm hoping to confirm we understand proper HSSM experimental inference, i.e., the right way to use HSSM to test for significant differences in DDM parameters among study factors.

We’re modelling 2 factors in a value-based decision-making experiment, where participants decide whether to keep (response = -1) or discard (response = 1) items:

1. 2 VALUE categories (items displayed can either be ‘TRASH’ or ‘VALUE’; this is a within-subjects factor where everyone has both conditions)

2. 2 GROUP categories (’HD’ or ‘HC’; a between-groups factor)

Previously, we used HDDM on this dataset. Now, we would like to use HSSM to validate those results in preparation for publication, while also capitalizing on HSSM’s ability to test for mixed between/within effects. So our questions here are to make sure we’re correctly testing for significant differences between our two study groups, for experimental conditions, or the between-within interaction effects.

In order to model the data, I used contrast coding to define the levels of each categorical variable: GROUP: HD (-0.5) and HC (0.5); VALUE: VALUE (0.5) and TRASH (-0.5).

Question 1: Is this model instantiated correctly to test for the effects of VALUE and GROUP on each DDM parameter?

lvl = [-0.5, 0.5]

model1 = hssm.HSSM(

model="ddm",

data=hd_data_subset,

extra_namespace={"lvl1": lvl1},

prior_settings = "safe",

include=[

{

"name": "v",

"formula": "v ~ C(group, levels=lvl1) * C(value, levels=lvl1)",

"link": "identity",

},

{

"name": "a",

"formula": "a ~ C(group, levels=lvl1) * C(value, levels=lvl1)",

"link": "identity",

},

{

"name": "z",

"formula": "z ~ C(group, levels=lvl1) * C(value, levels=lvl1)",

"link": "identity",

},

],

)

If anyone sees problems with the model, it would be great to correct our usage here before moving to inference. Next, to extract Bayesian P values for the GROUP, VALUE and interaction effects after sampling, we use the resulting model inference object to test whether the posterior distribution of the interaction term coefficient is > 0. (Previously in HDDM, we subtracted the two traces for each parameter (e.g., TRASH - VALUE) and tested whether the difference > 0.)

Question 2: Is testing the posterior distribution of the coefficient term this way the best approach to derive Bayesian p-values in HSSM? If so, is below the right code to compute this for the terms (e.g., in this example, the interaction)?

model1_idata = az.extract(model1_sample)

v_interaction_term = model1_idata['v_C(group, levels=lvl1):C(value, levels=lvl1)'].values

bayesian_p = (v_interaction_term > 0).mean()

Question 3: Are there any useful code snippets/other ways using HSSM/ArviZ/Bambi to conduct other Bayesian inference tests like testing whether one parameter is significant from 0 (e.g., az.plot_bpv)? Or, just anything we’re missing?

To expand, we foresee a need to maybe do post hoc checks on some individual conditions to make useful points (e.g., whether or not z significantly differs from .5 for each of the groups, etc)… things that just add interpretive context. Because HSSM adds capabilities relative to HDDM, we didn’t want to assume that the old “HDDM way” of testing for effects was the only approach, if indeed there are new and useful alternatives. And, having examples of successful syntax tailored to HSSM for various approaches can be very helpful.

As soon as we get confidence we’re approaching inference correctly, we plan to move on to using the model parameters to simulate data and do posterior predictive checks as our last step before writing this all up. So we might follow up with further questions. Thanks so much for your time!

[cuddenberg]

Hi, I just wanted to follow up on this issue -- any update on this? Thanks so much!

[AlexanderFengler]

Hey @cuddenberg sorry for the delay,

some initial responses,

1. Concerning Bayesian t-tests, you might find this tutorial from the HSSM docs useful.

2. Your model looks ok to me, but to make sure, could you send what it prints when you do print(model)?

3. One approach here is to run the posterior_predictive() with new data, and you set the prediction mode to mean (I will see if I can dig up a code snippet for that and/or add that to the existing tutorial I refer to above). Then you can pass a small dataframe that is 'representative' for your conditions, and you get the condition wise DDM parameter posteriors, from which you can then apply the Bayesian t-test logic.

Thanks for you interested in HSSM. We are trying hard to expand the coverage of our docs with more and more extended examples.

If you figure out a workflow that works for you independently, please know that we are very happy to receive contributed tutorials and will acknowledge you appropriately!

Best,
Alex

[cuddenberg]

Hi Alexander,

Thanks so much for getting back to me and for your help -- I answered each of your points below:

1. Thank you so much for this helpful tutorial! For our dataset, is one approach preferable over the other—fitting separate models for each condition versus using a single regression model for all conditions? Are there specific situations where one approach is better than the other for conducting Bayesian p-values? We want to ensure that we present these results properly when preparing for publication and peer review.

2. This is what print(model) outputs:

Hierarchical Sequential Sampling Model
Model: ddm
Response variable: rt,response
Likelihood: analytical Observations: 6727
Parameters:

v:

      Formula: v ~ C(group, levels=lvl1) * C(value, levels=lvl1)
      Priors:
            v_Intercept ~ Normal(mu: 0.0, sigma: 4.211900234222412)
            v_C(group, levels=lvl1) ~ Normal(mu: 0.0, sigma: 5.150599956512451)
            v_C(value, levels=lvl1) ~ Normal(mu: 0.0, sigma: 5.0)
            v_C(group, levels=lvl1):C(value, levels=lvl1) ~ Normal(mu: 0.0, sigma: 6.3805999755859375)
      Link: identity
      Explicit bounds: (-inf, inf)

a:
      Formula: a ~ C(group, levels=lvl1) * C(value, levels=lvl1)
      Priors:
            a_Intercept ~ Normal(mu: 0.0, sigma: 0.25)            
a_C(group, levels=lvl1) ~ Normal(mu: 0.0, sigma: 0.25)
            a_C(value, levels=lvl1) ~ Normal(mu: 0.0, sigma: 0.25)
            a_C(group, levels=lvl1):C(value, levels=lvl1) ~ Normal(mu: 0.0, sigma: 0.25)
      Link: identity
      Explicit bounds: (0.0, inf)

z:
      Formula: z ~ C(group, levels=lvl1) * C(value, levels=lvl1)
      Priors:
            z_Intercept ~ Normal(mu: 0.0, sigma: 0.25)             
z_C(group, levels=lvl1) ~ Normal(mu: 0.0, sigma: 0.25)
            z_C(value, levels=lvl1) ~ Normal(mu: 0.0, sigma: 0.25)
            z_C(group, levels=lvl1):C(value, levels=lvl1) ~ Normal(mu: 0.0, sigma: 0.25)
      Link: identity
      Explicit bounds: (0.0, 1.0)
t:
      Prior: HalfNormal(sigma: 2.0)
      Explicit bounds: (0.0, inf)

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

3. Just to clarify: so the idea is to use .sample_posterior_predictive() with kind='response_params' and a small dataframe that captures conditions of interest (e.g., different combinations of group and value) to generate new posterior samples for the parameters given that new representative data/predictors? And that returns posterior distributions for the parameters which we then apply Bayesian t-test logic to derive a Bayesian p-value?

I tried writing code to do this, let me know if this makes sense:

I first created a dataframe where each row of group and value are all the combinations of my conditions group and value: test_df = pd.DataFrame({
Group = [0.5, 0.5, -0.5, -0.5],
Value = [0.5, -0.5, 0.5, -0.5]

})
Then I passed this dataframe to .sample_posterior_predictive(kind = 'response_params') to generate what the model thinks the parameters would be for these conditions
model1.sample_posterior_predictive(data=test_df, kind = 'response_params', inplace = True, safe_mode = True)
                  
                  
Am I correct that the code below would give me the Bayesian p value for group = 0.5 in each of the two value conditions? (so for group =0.5, value = 0.5 then group = 0.5, value = -0.5)?

                 a_stacked = model1_groupxvalue_sample.posterior['a']

                 a_combined = a_stacked.stack(sample=("chain", "draw"))
                 a_diff = a_combined[:,1] - a_combined[:,2]
                 bayesian_p = np.mean(a_diff > 0)

Another question I had was with regards to .sample_prior_predictive() -- how is this check/function used in the Bayesian/HSSM workflow?

Thanks so much again for all your help!!

[cuddenberg]

Hi Alex, just following up on this—any updates? Thanks again for your help!

[johnmdusel]

3. One approach here is to run the posterior_predictive() with new data, and you set the prediction mode to mean (I will see if I can dig up a code snippet for that and/or add that to the existing tutorial I refer to above). Then you can pass a small dataframe that is 'representative' for your conditions, and you get the condition wise DDM parameter posteriors, from which you can then apply the Bayesian t-test logic.

It was a different modeling scanario, but FWIW, back in v0.2.3 I had difficulty understanding the output of sample_posterior_predictive with kind='mean'. I ended up working directly with posterior samples (kind='pps'), with fine results.