Hi,

I’m trying to build a multivariate hierarchical linear model, but I’m having trouble finding resources on how to do so. I can find examples of hierarchical univariate models and non-hierarchical multivariate models, but not the two together.

Fundamentally, I’m trying to build the following linear model (one per category) with shared hyperpriors between categories:

y_est = a[category] + b[category]*X

where category is an element of [0, 1, 2, 3] and X is an input matrix of num_training_samples by num_pred (10) covariates.

I’ve adapted the scikit learn wrapped linear model implementation at https://github.com/parsing-science/pymc3_models to perform this regression.

In the code below, which isn’t the full code, just the part that constructs the model, the following variables mean:

num_pred - number of predictors or covariates (10)

num_cats - number of categories for which to compute different coefficients (4)

model_cats - vector representing the category for each training sample [0, 1, 2, 3, 0, 1, 2, 3 … etc]

```
model_input = theano.shared(np.zeros([num_training_samples, num_pred]))
model_output = theano.shared(np.zeros(num_training_samples, dtype='float'))
model_cats = theano.shared(np.zeros(num_training_samples, dtype='int'))
model = pm.Model()
with model:
# Hyperpriors for group nodes
mu_alpha = pm.Normal('mu_alpha', mu=0, sd=100)
mu_beta = pm.Normal('mu_beta', mu=0, sd=100)
sigma_alpha = pm.HalfCauchy('sigma_alpha', 5)
sigma_beta = pm.HalfCauchy('sigma_beta', 5)
# model error
eps = pm.HalfCauchy('eps', 5)
# alpha and betas
alpha = pm.Normal('alpha', mu=mu_alpha, sd=sigma_alpha, shape=(num_cats,))
betas = pm.Normal('beta', mu=mu_beta, sd=sigma_beta, shape=(num_cats, num_pred))
c = model_cats
test_est = alpha[c] + T.sum(betas[c] * model_input, 1)
y_est = pm.Normal('y_est', test_est, sd=eps, observed=model_output)
```

So I am getting one alpha per category and all 4 alphas are sharing the hyperprior.

I don’t believe my beta hyperpriors are quite right. Since each beta represents a different covariate, it seems like it’s a mistake to set the same hyperpriors to be used for all of them. I would like to for my beta hyperpriors to be shared between categories, but not between different covariates.

I’ve made the following change to the shape to accomplish this:

```
mu_beta = pm.Normal('mu_beta', mu=0, sd=100, shape=(1, self.num_pred))
sigma_beta = pm.HalfCauchy('sigma_beta', 5, shape=(1, self.num_pred))
```

The model seems to be building and inference seems to work, but I am having some trouble getting it to converge.

Is my concept reasonable?

Is my modification to the code above the correct way to accomplish my goal?