Bambi's hierarchical modeling

[MatteoBenaud]

Hello everyone, I have a question regarding the theory behing bambi and specificly hierarchical modeling.

In traditional hierarchical Bayesian modeling, it is customary to incorporate hyper-priors to ensure that the resulting prior adequately captures the average effect across all categories.

However, upon reviewing the implementation in Bambi, for example here, I noticed that there is only one prior specified for sigma. This raised a question regarding the absence of a hyper-prior on mu.

Could you help me elucidate the reason of this choice ?
Is it possible to implement an hyperprior on mu ?

Thanks a lot !

[tomicapretto]

Hi @MatteoBenaud this is a great question and it's actually something I wanted to incorporate into our docs for a long time. I'll make an attempt at explaining how it works.

Bambi, as most (or all?) formula based interfaces, have some frequentist flavor when doing hierarchical (mixed, multi-level, etc.) modeling. In frequentist stats a mixed effects models for the mean (fixed and random effects) is written as follows

$$ \boldsymbol{\mu} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{Z}\boldsymbol{u} $$

The first dot product represents the fixed effects and the second dot product represents the random effects. Bambi is always following this same structure.

Let's present a simpler example. Let's imagine we have J groups and we want partially pooled intercepts. Using distributions we can write this model as

$$ \begin{aligned} Y_i &\sim \text{Normal}(m_{j[i]}, s) \\ m_j &\sim \text{Normal}(\mu, \sigma) & \text{for all } j = 1 \dots J\\ \mu &\sim \text{Normal}(...) \\ \sigma &\sim \text{HalfNormal}(...) \\ s &\sim \text{Something}(...) \end{aligned} $$

So each group has its own mean $m_j$, and these means all come from the same normal hyperprior with mean $\mu$ and standard deviation $\sigma$.

However, this is not the only way we can write this model. Have a look at the following:

$$ \begin{aligned} Y_i &\sim \text{Normal}(\mu + m'_{j[i]}, s) \\ m'_j &\sim \text{Normal}(0, \sigma) & \text{for all } j = 1 \dots J \\ \mu &\sim \text{Normal}(...) \\ \sigma &\sim \text{HalfNormal}(...) \\ s &\sim \text{Something}(...) \end{aligned} $$

Notice the $m'_j$ are now centered at 0. These models are equivalent. However, the second one allows for a different "view". We can think of $\mu$ as a common mean and $m'_j$ as a group specific deflection around. This is exactly what Bambi represents with 1 and (1 | Pig) in the model. 1 is the common intercept, and (1 | Pig) is the pig specific deviation around that common intercept. As in the math notation, the sum of these two components give the intercept for each pig.

This is the reason why you see (1 | Pig) ~ Normal(0, HalfNormal(...)). The group-specific effect is a deviation around the common effect. The same idea applies to slopes. Ah, and 1 goes into the X * beta part, and 1 | Pig goes into the Z * u part in the frequentist notation.

I don't want to extend too much here. I think the main idea is summarized above, but please feel free to ask for further clarification :)

[peeeffchang]

One example would be like this:
https://www.pymc-labs.com/blog-posts/causal-analysis-with-pymc-answering-what-if-with-the-new-do-operator/

[tomicapretto]

This one?

with pm.Model(coords_mutable={"i": [0]}) as model_generative:
    # priors on Y <- C -> Z
    beta_y0 = pm.Normal("beta_y0")
    beta_cy = pm.Normal("beta_cy")
    beta_cz = pm.Normal("beta_cz")
    # priors on Z -> Y causal path
    beta_z0 = pm.Normal("beta_z0")
    beta_zy = pm.Normal("beta_zy")
    # observation noise on Y
    sigma_y = pm.HalfNormal("sigma_y")
    # core nodes and causal relationships
    c = pm.Normal("c", mu=0, sigma=1, dims="i")
    z = pm.Bernoulli("z", p=pm.invlogit(beta_z0 + beta_cz * c), dims="i")
    y_mu = pm.Deterministic("y_mu", beta_y0 + (beta_zy * z) + (beta_cy * c), dims="i")
    y = pm.Normal("y", mu=y_mu, sigma=sigma_y, dims="i")

pm.model_to_graphviz(model_generative)

[peeeffchang]

Yep that's the one describing the system

[tomicapretto]

I'm afraid it's not possible to implement that model in Bambi because of:

• The multiple likelihoods (c, z, and y)
• The dependency between them (e.g. z and c go into y_mu).

It's possible one could write a custom family that hacks through Bambi internals and implements this PyMC model. But if one is willing to do such amount of work, it's just simpler to use plain PyMC.

[peeeffchang]

@tomicapretto Noted thanks