Broadcasting error when using default intercept (1) on DDM model

[morgbead]

I am currently using HSSM on Google Colab, and have been working on tweaking some analyses over the last month or two using this package. The code below was working to define the model and sample the data up until recently - it is now producing a broadcasting error, specific verbiage being:

[/usr/local/lib/python3.12/dist-packages/pytensor/tensor/elemwise.py](https://localhost:8080/#) in _check_runtime_broadcast(node, inputs)
    740         ):
    741             if any(d != 1 for d, _ in dims_and_bcast) and (1, False) in dims_and_bcast:
--> 742                 raise ValueError(
    743                     "Runtime broadcasting not allowed. "
    744                     "At least one input has a distinct dimension length of 1, but was not marked as broadcastable.\n"

ValueError: Runtime broadcasting not allowed. At least one input has a distinct dimension length of 1, but was not marked as broadcastable.
If broadcasting was intended, use `specify_broadcastable` on the relevant input.

Here is the model definition and sampling code that was previously working:

# Model Definition
model1 = hssm.HSSM(
    model="ddm", )
    noncentered=False, 
    data=dat, # dataset
    prior_settings="safe", 
    include=[
        {
            "name": "v",
            "formula": "v ~ group",
            "link": "identity",
        },
        {
            "name": "a",
            "formula": "a ~ group",
            "link": "identity",
        },
        {
            "name": "z",
            "formula": "z ~ group",
            "link": "identity", 
        },
    ],
)

# Model Sampling
idata_m1 = model1.sample(
    sampler="nuts_numpyro",  
    cores=2,  
    chains=4,  
    draws=8000,  
    tune=4000,  
    target_accept=0.95,  
    idata_kwargs=dict(log_likelihood=True), 
)

The model will sample if I change the default intercept to 0 (in the model definition, alter the formula for each of the 3 parameters to be ~ 0 + group). It will also sample if I do not define the group parameter in the model, and instead use the default model settings, such as - model1 = hssm.HSSM(data = dat, model = "ddm"). It will not sample if I specify the formula to have an intercept of 1 (in the model definition, alter the formula for each of the 3 parameters to be ~ 1 + group).

Given the recent onset of this error on code that was previously working, I suspect that this is due to a recent update. I would be grateful if someone could offer advice on a way to define my model in a way that circumvents the error.

Happy to elaborate on any of the above! Thank you in advance for your time.

[AlexanderFengler]

@morgbead does group have only two values?
In that case you would not want to specify an intercept (so setting to 0 + group is correct).

[morgbead]

@AlexanderFengler group does have two values. My understanding was that the intercept of 1 should be used unless there is a theory-driven reason not to use it, so this is valuable advice. Is there any general guidance about when to use a 0 vs another value for the intercept?

For example, in some iterations of this analysis, we use an interaction of group and another two-value variable - should the intercept of that model also be 0?

[AlexanderFengler]

In general, if you treat the variable as a dummy (so not as categorical, but [0,1] coded as a continuous covariate) you should be ok, because it assigns one beta.

If it is actually a categorical variable, (to make that explicit you can use the C(my_variable) coding in the regression model) and treated as such, you should run into the problem above.

This should translate to the interaction case too.

Definitely suggest to check the Bambi docs as well.

[morgbead]

Hi Alexander,

Thank you for your insight! I tested a few different approaches to the syntax of this model, and I'd appreciate any thoughts or guidance you might have. Apologies for the length of this message, but we wanted to be thorough. While the details in this post are related in process to this thread, the issues I'm looking for guidance on are likely of a more general problem. If this query should be moved to a different part of the list serv, I'd be happy to readjust the location of the post to ensure that others can benefit from any advice.

The example I gave above was a dummy / one-hot variable that is within subject. When it is coded as categorical and the intercept is 1, the error I initially mentioned (the 'broadcasting array' message) appeared and the model failed to sample. Using the same categorical coding for the variable, but defining the model to have an intercept of 0, allows the model to sample - however, it produces the maximum number of divergences when doing so, and the model summary table statistics confirms a complete lack of model convergence (Gelman-Rubins values of NaN). When this variable is coded numerically (0/1) and defined as a numeric/float variable type, an intercept of 0 produces the same pattern seen before (maximum divergences in model sampling, Gelman-Rubin values of NaN). But now, an intercept of 1 appears to be yielding interpretable results. This is accompanied by a large number of divergences (15,138 divergences), but all Gelman-Rubin values in the summary table are <1.1.

When I adjust the model definition to incorporate a between subjects interaction variable, an intercept of 0 produces the same pattern seen above (maximum divergences, Gelman-Rubins values of NaN). The model will run when the intercept for this is a value of 1, as long as one of the two variables has been coded as 0/1 and defined as a numeric/float type (the other can remain a categorical variable). Of interest is that this doesn't change the model output - if both variables are defined as numeric, the model summary details are nearly identical (within 0.001) of the results from when one is categorical and the other is a factor.

The above description is a shortened recap of a few days' worth of careful model testing and retesting - hopefully the details are helpful in contextualizing where we're at, but please let me know if additional information would be productive in that regard. Ultimately, there are 3 items that I'd love further clarification on, which would help us to accurately move forward in using this package for a slew of interesting data that we have lined up. We recently designed a pilot study for a clinical population around DDM-compatible tasks, and are excited to examine some early results. To do so, we want to be confident in our model set-up and evaluation, which these three clarifications will help to provide:

1. Divergences. As noted above, we have produced models that have a concerningly high number of divergences (e.g. 15138/16000 (where 16000 was the total samples for the model); 28802/32000; 71315/80000). I have listed below the series of attempts we've made to decrease divergences. Of note is that all of the models appear to still be converging: the trace plots and marginal posteriors look okay, and the Gelman-Rubin values are all below 1.1. Employing the different approaches to minimize divergences (listed below) does move the needle in the right direction by decreasing total divergences, but they remain concerningly high.
   a. increasing the draws and tuning
   b. increasing the 'target_accept' value (from 0.9-0.99)
   c. changing the model definition parameter of 'noncentered' to 'False' (a recommendation from this thread: hierarchical DDM does not seem to converge at all on 0.2.2 #471 to address funnel-shapes in hierarchical models)
   d. using the ‘thin’ argument, in case autocorrelation issues were the root cause
   e. excluding participants from the data set that had highly skewed choices, decreasing the ‘maximum skew’ from 90/10 to 75/25 and 70/30
   f. selectively excluding trials where participants took a long time (>3 sec) to respond
   g. These attempts have helped to decrease the total number of divergences, but not nullify them. All other model checks appear to be reflecting a sound model
   We’re hoping for some guidance on resolving these issues with divergences for our models – or, alternatively, information on whether this is incompatible with the current toolbox, and any recommendations for where to turn to if that’s the case?

2. Syntax. For the model I've described throughout this post, we are hoping to assess: the main effect of a single within subjects factor; interaction effect for a within and between subjects factor; and an interaction effect of 2 within subjects factors. Our testing has pointed to the use of the value of 1 for an intercept (given that a 0 intercept does not produce model results), is that the correct choice? For the syntax of the model definition, is our approach appropriate? Example syntax below ('group' is our current within-subjects factor, 'between' is an example of a between subjects factor)

model1 = hssm.HSSM(
model="ddm",
noncentered=False,
data=dat,
prior_settings="safe",
include=[
{
"name": "v",
"formula": "v ~ 1 + group * between",
"link": "identity",
},
{
"name": "a",
"formula": "a ~ 1 + group* between",
"link": "identity",
},
{
"name": "z",
"formula": "z ~ 1 + group* between ",
"link": "identity",
},
],
)

3. Coding categorical variables. Is there any detriment to generally encoding binary categorical variables as 0/1 numeric-coded variables for use in this model? And, are there any barrier to interpretation when we do so? I presume the answer is that it is okay to dummy/one-hot encode them generally and there aren't interpretation limitations when doing so (given the similarity of the interaction results that I mentioned earlier when including both variables as numeric VS one numeric and one group) - but I wanted to confirm this before proceeding.

We're torn between whether these issues- particularly the divergences issue -is related to faulty install of HSSM, or (more probably) some simple user error on our part. We're hoping that someone with more insight to the general mechanisms of the package might be able to provide some insight to these issues.

Many thanks,

Morgan