I’m trying to create a multivariable linear regression model with multiple groups. The closest I could find online was this resource from one of the authors:

http://twiecki.github.io/blog/2014/03/17/bayesian-glms-3/

Although this one creates multiple groups there is only a single metric predictor. I think the code below creates what I want, but it seems kind of contrived since I end up using list comprehension.

```
model_input_x = [theano.shared(X_train[:,i].reshape(-1,1)) for i in range(0,X_train.shape[1])]
model_input_cat = theano.shared(np.asarray(X_train[:,-1], dtype='int32').reshape(-1,1))
model_output = theano.shared(y_train)
with pm.Model() as hierarchical_model:
mu_beta = pm.Normal('mu_beta', mu=0., sd=1, shape=(num_features))
sigma_beta = pm.HalfCauchy('sigma_beta', beta=1, shape=(num_features))
beta = [pm.Normal('beta_'+str(i), mu=mu_beta[i], sd=sigma_beta[i], shape=(num_hierarchy)) for i in range(0,num_features)]
alpha = pm.Normal('alpha', mu=0, sd=1, shape=(1))
epsilon = pm.HalfCauchy('eps', beta=1) # Model error prior
y_est = alpha + sum([beta_i[model_input_cat]*model_input_x[i] for i,beta_i in enumerate(beta)])
y_like = pm.Normal('y_like', mu=y_est, sd=epsilon, observed=model_output)
```

Is this the best way to build a hierarchical model where each feature/group combination has it’s own slope and slopes across groups for a given feature come from a common distribution?

I think this is the model I’m creating with this code (given 2 features). The pm.model_to_graphviz tool output is also shown