How to recreate a model in Bambi in pymc

Questions v5
1

I came across the Bambi library and have been playing with it using the sleepstudy data.

Using a relatively simple mode: bmb.Model("Reaction ~ 0 + Days + (1 | Subject)", sleep), how might I recreate it using base pymc? The bambi model produces the following graph

Screenshot 2023-03-25 at 8.02.18 PM
Screenshot 2023-03-25 at 8.02.18 PM1578×950 44.6 KB

which, I can’t seem to figure out how to reproduce in base pymc. In particular, I’m confused how the deterministic block functions in this case.

2

The deterministic block is implementing the “non-centered parameterization” of a hierarchical prior. There’s a nice discussion of this in the radon example notebook. McElreath also has a detailed discussion of why you want to do this in this video.

2 Likes
3

This is the source code of a talk I gave showing how to do the same model with both PyMC and Bambi

Let me know if that helps :slight_smile:

Edit: This may also help unfinished-ideas/sparse.ipynb at main · bambinos/unfinished-ideas · GitHub

Notice the goal there was to replace some indexing operation with a sparse dot product. But there you have how to reproduce the model with PyMC

1 Like
4

Just to be safe, would this be the corresponding pymc model to the bambi model from above?
Note that it has an intercept (Reaction ~ 1 + days + (1 | subject ))

with pm.Model(coords=coords) as model_pm:
    # Priors for the intercept
    intercept_mu = pm.Normal("intercept_mu", mu=0, sigma=250)
    intercept_sigma = pm.HalfNormal("intercept_sigma", sigma=200)
    intercept_offset = pm.Normal("intercept_offset", dims="subjects")
    intercept = pm.Deterministic("intercept", intercept_mu + intercept_sigma * intercept_offset)

    # Priors for the Days coefficient
    days_mu = pm.Normal("days_mu", mu=0, sigma=50)
    days_sigma = pm.HalfNormal("days_sigma", sigma=50)

    sigma = pm.HalfStudentT("sigma", nu=4, sigma=50)

    # Linear model
    mu = intercept[subject_idx] + days_mu * data.Days.values
    mu = pm.Deterministic("mu", mu)

    # likelihood
    reaction = pm.Normal("reaction", mu=mu, sigma=sigma, observed=data.Reaction.values)

    # Sample
    idata_pm = pm.sample(return_inferencedata=True)