Bayesian Hierarchical Modeling - report coefficients for multiple levels

This line will not give you what you predictions per se, because you have declared that variable to be observed (and is redundant with the preceding line):

pred = pm.Poisson('Y_pred_city', mu=beta_0_city[7] + (beta_1_city[7] * X['act1']), observed=y)`

It’s likely that you will want to interrogate the trace object generated during sampling to generate your predictions, though the details will depend on exactly what you are trying to predict (new data from city 7? counterfactual conditions in which city 7 had different characteristics?). To just get the predicted mean for city 7, for example, you can just grab sampled values of beta_0_city[7] and beta_1_city[7] and plug them into scipy.stats.poisson. We had some related discussion here. You might also take a look at the documentation for pymc3.model.set_data(). That’s sometimes cleaner/easier depending on your needs (e.g., when trying to predict out of sample).