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
radon_like = pm.Normal('radon_like', mu=radon_est, sigma=eps, observed=data.log_radon)
with hierarchical_model:
hierarchical_trace = pm.sample(draws=2000, n_init=1000)
What would be a similar set up in pyro? Can I pass a distribution inside of another distribution, as done here? This is done to pull the individual clusters towards a common mean.
And just to be clear, the array looks like this