# How to implement Multilevel/Hierarchical Regressions?

Hello,

I am trying to run a hierarchical regression model (sometimes multilevel, or mixed effects) and I am having trouble getting off the ground.

I am hoping to be able to set it up in the spirit of Pymc3, as in this post:

https://docs.pymc.io/notebooks/GLM-hierarchical.html

Here is the code example I am thinking of, in Pymc3. Here county_idx is an array indicating which group/cluster each intercept should belong to. This would be considered a random intercepts regression with random slopes as well.

``````county_idx = data.county_code.values
n_counties = len(data.county.unique())

with pm.Model() as hierarchical_model:
# Hyperpriors for group nodes
mu_a = pm.Normal('mu_a', mu=0., sigma=100**2)
sigma_a = pm.HalfCauchy('sigma_a', 5)
mu_b = pm.Normal('mu_b', mu=0., sigma=100**2)
sigma_b = pm.HalfCauchy('sigma_b', 5)

# Intercept for each county, distributed around group mean mu_a
# Above we just set mu and sd to a fixed value while here we
# plug in a common group distribution for all a and b (which are
# vectors of length n_counties).
a = pm.Normal('a', mu=mu_a, sigma=sigma_a, shape=n_counties)
# Intercept for each county, distributed around group mean mu_a
b = pm.Normal('b', mu=mu_b, sigma=sigma_b, shape=n_counties)

# Model error
eps = pm.HalfCauchy('eps', 5)

radon_est = a[county_idx] + b[county_idx] * data.floor.values

# Data likelihood