Problem with multivariate hierarchical linear model

Questions
#1

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 :slight_smile:

1 Like
Sampling semantics of multiple observed variables
#2

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.

#3

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?

#4

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