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),
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?
Usually you can make it work by specifying the dimension of the random effect like:
b = pm.Normal('b', mu=beta, sd=sigma_b, shape=(n_domains, 1)), and do matrix multiplication along the correct dimension.
For more information, have a look at my repository https://github.com/junpenglao/GLMM-in-Python and @Jack_Caster’s repository https://github.com/JackCaster/GLM_with_PyMC3
There are quite a few examples in there.
Ok thanks, so it would be defined like this,
beta = pm.Normal('beta', mu=0, sd=1, shape=(1,2))
beta_print = T.printing.Print('beta')(beta.shape)
b = pm.Normal('b', mu=beta, sd=1, shape=(3,2))
b_print = T.printing.Print('b')(b.shape)
The shape of beta is like 1 sample, with 2 covariates. And then for b we just broadcast the 1 sample up to 3?
Not sure what you mean by 1 sample - but just to be clear the shape of the coefficients (b and beta) are determined by the shape of the covariates (i.e., design matrix).