I have a dataset with 2 covariates (i.e. X.shape = (600, 2)) and 1 dependent variable (y). As my data comes from two populations I want to use a mixed effects / hierarchical model. The non-hierarchical model looks like this (just treating all the data as equivalent):
alpha = pm.Normal('alpha', mu=0, sd=1)
beta = pm.Normal('beta', mu=0, sd=1, shape=2)
# Expected value of outcome
mu = alpha + T.dot(beta, minibatch_x.T)
Y_obs = pm.Bernoulli('Y_obs', p=logistic(mu),
observed=minibatch_y, total_size=ys.shape[0])
Iām struggling to make this into a hierarchical model. Basically the problem is that because I have two covariates the beta of my beta above is already 2. So when I come to defining my random effects, i.e.
b = pm.Normal('b', mu=beta, sd=sigma_b, shape=n_domains)
The shape is messed up. Does anybody know a way of doing what I want to achieve - I feel like it must be a fairly common issue?
Thanks