Hierarchical modeling¶

This tutorial demonstrates how to take advantage of HSSM's hierarchical modeling capabilities. We will cover the following:

• How to define a mixed-effect regression
• How to define a hierarchical HSSM model
• How to fit the model and read the group-level summary

Run this tutorial¶

On Colab, uncomment and run the installation cell below once, then restart the runtime. For local setup, GPU extras, and troubleshooting see the Installation guide.

In [1]:

# %pip install hssm

# %pip install hssm

Import Modules¶

In [2]:

import arviz as az

import hssm

%matplotlib inline

import arviz as az import hssm %matplotlib inline

1. Defining Regressions¶

Under the hood, HSSM uses bambi for model creation. bambi takes inspiration from the lme4 package in R and supports the definition of generalized linear mixed-effect models through R-like formulas and concepts such as link functions. This makes it possible to create arbitrary mixed-effect regressions in HSSM, which is one advantage of HSSM over HDDM. Now let's walk through the ways to define a parameter with a regression in HSSM.

Specifying fixed- and random-effect terms¶

Suppose that we want to define a parameter v that has a regression defined. There are two ways to define such a parameter - either through a dictionary or through a hssm.Param object:

# The following code are equivalent,
# including the definition of the formula.

# The dictionary way:
param_v = {
    "name": "v",
    "formula": "v ~ x + y + x:y + (1|participant_id)",
    "link": "identity",
    "prior": {
        "Intercept": {"name": "Normal", "mu": 0.0, "sigma": 0.25},
        "1|participant_id": {
            "name": "Normal",
            "mu": 0.0,
            "sigma": {"name": "HalfNormal", "sigma": 0.2},  # this is a hyperprior
        },
        "x": {"name": "Normal", "mu": 0.0, "sigma": 0.25},
    },
}

# The object-oriented way
param_v = hssm.Param(
    "v",
    formula="v ~ 1 + x*y + (1|participant_id)",
    link="identity",
    prior={
        "Intercept": hssm.Prior("Normal", mu=0.0, sigma=0.25),
        "1|participant_id": hssm.Prior(
            "Normal",
            mu=0.0,
            sigma=hssm.Prior("HalfNormal", sigma=0.2), # this is a hyperprior
        ),
        "x": hssm.Prior("Normal", mu=0.0, sigma=0.25),
    },
)

The formula "v ~ x + y + x:y + (1|participant_id)" defines a random-intercept model. Like R, unless otherwise specified, a fixed-effect intercept term is added to the formula by default. You can make this explicit by adding a 1 to the formula. Or, if your regression does not have an intercept. you can explicitly remove the intercept term by using a 0 in the place of 1: "v ~ 0 + x * y + (1|participant_id)". We recommend that the random effect terms should be specified after the fixed effect terms.

Other fixed effect covariates are x, y, and the interaction term x:y. When all three terms are present, you can use the shortcut x * y in place of the three terms.

The only random effect term in this model is 1|participant_id. It is a random-intercept term with participant_id indicating the grouping variable. You can add another random-effect term in a similar way: "v ~ x + y + x:y + (1|participant_id) + (x|participant_id)", or more briefly, "v ~ x + y + x:y + (1 + x|participant_id)".

Specifying priors for fixed- and random-effect terms:¶

As demonstrated in the above code, you can specify priors of each term through a dictionary, with the key being the name of each term, and the corresponding value being the prior specification, etiher through a dictionary, or a hssm.Prior object. There are a few things to note:

• The prior of fixed-effect intercept is specified with "Intercept", capitalized.
• For random effects, you can specify hyperpriors for the parameters of of their priors.

Specifying the link functions:¶

Link functions is another concept in frequentist generalized linear models, which defines a transformation between the linear combination of the covariates and the response variable. This is helpful especially when the response variable is not normally distributed, e.g. in a logistic regression. In HSSM, the link function is identity by default. However, since some parameters of SSMs are defined on (0, inf) or (0, 1), link function can be helpful in ensuring the result of the regression is defined for these parameters. We will come back to this later.

2. Defining a hierarchical HSSM model¶

HSSM does not differentiate between hierarchical and non-hierarchical models — a hierarchical model is simply a model in which one or more parameters are defined as regressions with group-specific terms.

Applying a formula to every parameter with global_formula¶

When specified, global_formula adds the given formula to every parameter that is not explicitly defined by the user. Setting it to "y ~ 1 + (1|participant_id)" therefore gives every parameter a random intercept per participant — the classic fully hierarchical setup. y stands for the name of each parameter; HSSM is agnostic to the name left of the ~ sign, y is simply customary. You can supply any formula you want applied across parameters.

Note

In HSSM, the default grouping variable is `participant_id`.

In [3]:

# Load a package-supplied dataset
cav_data = hssm.load_data("cavanagh_theta")

# Load a package-supplied dataset cav_data = hssm.load_data("cavanagh_theta")

In [4]:

# Define a basic non-hierarchical model
model_non_hierarchical = hssm.HSSM(data=cav_data)
model_non_hierarchical

# Define a basic non-hierarchical model model_non_hierarchical = hssm.HSSM(data=cav_data) model_non_hierarchical

Model initialized successfully.

Out[4]:

Hierarchical Sequential Sampling Model
Model: ddm

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

Parameters:

v:
    Prior: Normal(mu: 0.0, sigma: 2.0)
    Explicit bounds: (-inf, inf)

a:
    Prior: HalfNormal(sigma: 2.0)
    Explicit bounds: (0.0, inf)

z:
    Prior: Uniform(lower: 0.0, upper: 1.0)
    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)

In [5]:

# Give every parameter a random intercept per participant
model_hierarchical = hssm.HSSM(
    data=cav_data, global_formula="y ~ 1 + (1|participant_id)"
)
model_hierarchical

# Give every parameter a random intercept per participant model_hierarchical = hssm.HSSM( data=cav_data, global_formula="y ~ 1 + (1|participant_id)" ) model_hierarchical

Model initialized successfully.

Out[5]:

Hierarchical Sequential Sampling Model
Model: ddm

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

Parameters:

v:
    Formula: v ~ 1 + (1|participant_id)
    Priors:
        v_Intercept ~ Normal(mu: 2.0, sigma: 3.0)
        v_1|participant_id ~ Normal(mu: 0.0, 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)

3. Fitting and checking the hierarchical model¶

Specification is only half the story — let's sample the hierarchical model we just built and look at the group-level summary. (omit_offsets=True keeps the non-centered offset nuisance variables out of the results.)

In [6]:

infer_data = model_hierarchical.sample(
    draws=500, tune=500, chains=4, target_accept=0.9, omit_offsets=True
)

infer_data = model_hierarchical.sample( draws=500, tune=500, chains=4, target_accept=0.9, omit_offsets=True )

Using default initvals. 

Initializing NUTS using adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [v_Intercept, v_1|participant_id_sigma, v_1|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]

Sampling 4 chains for 500 tune and 500 draw iterations (2_000 + 2_000 draws total) took 243 seconds.

In [7]:

az.summary(
    model_hierarchical.traces, var_names=["_Intercept", "_sigma"], filter_vars="like"
)

az.summary( model_hierarchical.traces, var_names=["_Intercept", "_sigma"], filter_vars="like" )

Out[7]:

	mean	sd	eti89_lb	eti89_ub	ess_bulk	ess_tail	r_hat	mcse_mean	mcse_sd
v_Intercept	0.448	0.08	0.32	0.57	701	1011	1.00	0.003	0.0023
z_1|participant_id_sigma	0.0198	0.0098	0.0059	0.037	809	1027	1.01	0.00033	0.00027
t_Intercept	0.478	0.031	0.43	0.53	627	767	1.01	0.0012	0.001
a_1|participant_id_sigma	0.167	0.037	0.12	0.23	903	1148	1.01	0.0012	0.0012
v_1|participant_id_sigma	0.281	0.06	0.2	0.39	826	922	1.00	0.002	0.0017
z_Intercept	0.4909	0.008	0.48	0.5	1811	1476	1.00	0.00019	0.00015
t_1|participant_id_sigma	0.116	0.026	0.083	0.16	749	981	1.00	0.00095	0.00099
a_Intercept	0.988	0.043	0.92	1.1	564	764	1.01	0.0018	0.0015

The rows to read first are the group-level standard deviations (the _sigma terms): they quantify how much each parameter varies across participants — the whole point of going hierarchical.

Next steps¶

• Hierarchical DDM — the advanced deep-dive: within-subject, between-subject, and interaction designs mapped to formulas, with full parameter recovery.
• Centered vs. non-centered parameterizations and per-parameter control — when hierarchical models sample badly, this is the first knob to understand.
• How to specify priors and fix parameters — including the prior_settings/link_settings safeguards for hierarchical models.
• Variational inference — when MCMC on a large hierarchical model is too slow.