Omitted variable inference after HSSM update

[SaschaFroelich]

I updated HSSM from 0.2.1 to 0.2.3 (and numpyro: 0.14.0 -> 0.15.0, pymc 5.15.1 -> 5.16.1). Now the model omits inference of the different jokerconditions of the z parameter. Here's the model:

model_reg_v_angle_hier = hssm.HSSM(
    data = ddmdata,
    model = "angle",
    hierarchical = False,
    categorical = 'jokercondition',
    include=[
        {
            "name": "v",
            "formula": "v ~ 1 + jokercondition",
            "prior": {
                "Intercept": {"name": "Normal", "mu": 1, "sigma": 2, "initval": 1},
            },
            "link": "identity",
        },
        {
            "name": "z",
            "formula": "z ~ 1 + jokercondition",
            "prior": {
                "Intercept": {
                    "name": "Uniform",
                    "lower": 0.3,
                    "upper": 0.7,
                    "initval": 0.5},
            },
            "link": "identity",
        }
    ],
)

This is how the inference result would look before (only 10 samples for testing, so no convergence):

                        mean     sd  hdi_3%  ...  ess_bulk  ess_tail  r_hat
a                      1.180  0.125   0.967  ...      12.0      20.0   2.27
t                      0.060  0.032   0.028  ...      11.0      20.0   2.43
theta                  1.125  0.074   0.973  ...      11.0      20.0   2.71
v_Intercept            2.406  0.278   1.929  ...      11.0      20.0   2.39
v_jokercondition[2.0]  0.040  0.099  -0.101  ...      13.0      64.0   1.95
v_jokercondition[3.0] -0.109  0.191  -0.335  ...      12.0      15.0   2.29
z_Intercept            0.461  0.021   0.436  ...      14.0      61.0   1.71
z_jokercondition[2.0]  0.049  0.020   0.026  ...      16.0      35.0   1.46
z_jokercondition[3.0] -0.055  0.030  -0.108  ...      21.0      20.0   1.26

This is how the inference comes out now (again only 10 samples):

                        mean     sd  hdi_3%  ...  ess_bulk  ess_tail  r_hat
a                      1.101  0.137   0.950  ...      11.0      20.0   3.19
t                      0.095  0.024   0.055  ...      12.0      20.0   2.14
theta                  1.136  0.063   1.033  ...      12.0      28.0   2.28
v_Intercept            2.555  0.294   1.974  ...      11.0      20.0   2.31
v_jokercondition[2.0]  0.084  0.051  -0.012  ...      16.0      47.0   1.43
v_jokercondition[3.0] -0.022  0.117  -0.251  ...      12.0      28.0   2.12

And then results for the individual z in each trial. So basically, z_Intercept, z_jokercondition[2.0]and z_jokercondition[3.0]are suddenly not inferred anymore?

Best, Sascha

[digicosmos86]

Hi Sascha,

Can you also include a printout of the model? You can do print(model_reg_v_angle_hier) and `print(model_reg_v_angle_hier.model). It would be helpful to know if hssm built the model correctly. One other thing you can try is to specify C(jokercondition) in the formulas to indicate that jokercondition is a categorical variable. I am not 100% clear how the categorical parameter works in bambi, so explicitly using C() notation to indicate categorical variables provides one additional guarantee.

Thanks! Paul

[SaschaFroelich]

Hi Paul,

print(model_reg_v_angle_hier)
Hierarchical Sequential Sampling Model
Model: angle

Response variable: rt,response
Likelihood: approx_differentiable
Observations: 611

Parameters:

v:
    Formula: v ~ 1 + jokercondition
    Priors:
        v_Intercept ~ Normal(mu: 1.0, sigma: 2.0, initval: 1.0)
        v_jokercondition ~ Normal(mu: 0.0, sigma: 0.25)
    Link: identity
    Explicit bounds: (-3.0, 3.0)

a:
    Prior: Uniform(lower: 0.3, upper: 3.0)
    Explicit bounds: (0.3, 3.0)

z:
    Formula: z ~ 1 + jokercondition
    Priors:
        z_Intercept ~ Uniform(lower: 0.3, upper: 0.7, initval: 0.5)
        z_jokercondition ~ Normal(mu: 0.0, sigma: 0.25)
    Link: identity
    Explicit bounds: (0.1, 0.9)

t:
    Prior: Uniform(lower: 0.001, upper: 2.0)
    Explicit bounds: (0.001, 2.0)

theta:
    Prior: Uniform(lower: -0.1, upper: 1.3)
    Explicit bounds: (-0.1, 1.3)

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

Specifying C(jokercondition) does not seem to make any difference.

Best, Sascha

[digicosmos86]

How about print(model_reg_v_angle_hier.model)?

[SaschaFroelich]

Oh sorry, I forgot! Here it is:

print(model_reg_v_angle_hier.model)
       Formula: c(rt, response) ~ 1 + jokercondition
                z ~ 1 + jokercondition
        Family: SSM Family
          Link: v = identity
                z = identity
  Observations: 611
        Priors: 
    target = v
        Common-level effects
            Intercept ~ Normal(mu: 1.0, sigma: 2.0, initval: 1.0)
            jokercondition ~ Normal(mu: 0.0, sigma: 0.25)

        Auxiliary parameters
            theta ~ Uniform(lower: -0.1, upper: 1.3)
            p_outlier ~ 0.05
            a ~ Uniform(lower: 0.3, upper: 3.0)
            t ~ Uniform(lower: 0.001, upper: 2.0)
    target = z
        Common-level effects
            z_Intercept ~ Uniform(lower: 0.3, upper: 0.7, initval: 0.5)
            z_jokercondition ~ Normal(mu: 0.0, sigma: 0.25)
------
* To see a plot of the priors call the .plot_priors() method.
* To see a summary or plot of the posterior pass the object returned by .fit() to az.summary() or az.plot_trace()

[digicosmos86]

Hi Sascha,

Thank you for providing the information! It seems to me that the formulas specified are still the same in these outputs. Do you mean that after specifying the formulas for v and z as v ~ 1 + C(jokercondition) and z ~ 1 + C(jokercondition), nothing has changed?

A few other things that might be helpful to try:

1. Actually code jokercondition in strings, which will ensure that HSSM will not treat it as continuous.
2. Check if the levels for jokercondition exist in the InferenceData object. I'd suspect that they probably don't since the summary doesn't show it, but it might still be a good idea to check.
3. Specify prior_settings to None, for a simple model that you are trying out

Please let me know if any of these helps

[SaschaFroelich]

Hi Paul,

thanks for getting back. Your suggestion to look in the InferenceData object was spot on. Indeed, z_jokercondition was there, and I can print its inference results with print(az.summary(infer_data_trace, var_names = [ 'z_jokercondition'])). Simply printing the summary print(az.summary(infer_data_trace)) puts it at the very end (after the inferences for the individual trials), which was new and confusing. This is something I should have been able to solve by myself, so my apologies for wasting your time. And thank you for your assistance.

Best, Sascha

[digicosmos86]

Hi Paul,

thanks for getting back. Your suggestion to look in the InferenceData object was spot on. Indeed, z_jokercondition was there, and I can print its inference results with print(az.summary(infer_data_trace, var_names = [ 'z_jokercondition'])). Simply printing the summary print(az.summary(infer_data_trace)) puts it at the very end (after the inferences for the individual trials), which was new and confusing. This is something I should have been able to solve by myself, so my apologies for wasting your time. And thank you for your assistance.

Best, Sascha

Hi Sascha,

I am glad everything worked! In fact, we have a wrapped .summary() method for the HSSM class that has some built-in filters to help folks get rid of the deterministics in the InferenceData. You can just call model.summary() after sampling. You don't even need to import arviz :)