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
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.
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.
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 
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
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)