Your model looks mostly correct to me, though you may want to check that your predictors are on the same scale and if not standardise them (your model does not seem to be built to deal with this since priors for slopes seem to have the same sigma for every predictor). If I am not mistaken, you also defined a separate intercept for every predictor. In my experience the more common way is to define a single intercept and multiple slopes (which would be multivariate linear regression basically). See for instance Chapter 18 in Kruscheke talks about these. There is a notebook which has translated the code there to pymc:
https://nbviewer.org/github/cluhmann/DBDA-python/tree/master/Notebooks/
For doing predictions see the following link (I think it has linear regression examples, though not hierarchical, but should easily generalize to your model)
For checking the model, usually I would do some checks roughly in the following order (some of them would be only useful if you do 4 or more chains):
1- You can construct a diagram of your model:
https://www.pymc.io/projects/docs/en/v5.6.1/api/generated/pymc.model_to_graphviz.html
This would make it easier to eyeball your model and see if it is defined the way you imagined it.
2- Prior and posterior checks: Real test of performance for Bayesian models is the posterior predictive check:
https://www.pymc.io/projects/docs/en/stable/learn/core_notebooks/posterior_predictive.html
Prior checks can also be used, to see if the priors you have chosen give results within reasonable values. This is also shown in the link above.
Note that when you do
sample_posterior_predictive(trace, extend_inference=True, predictions=True)
then your trace object will contain sampled predictions in trace[“predictions”][“Y_obs”]. This will have the shape number of draws x number of chains x Y_obs.shape. If you want to you can take the mean across the first two dimensions and do a residual plot maybe with error bands. However see 4 below.
3- Use arviz to produce posterior plots and trace plots for your parameters. Your posteriors from different chains should produce similar results. And also I think in a simple linear regression setting one would expect posteriors to look like unimodal normal distributions.
https://python.arviz.org/en/latest/api/generated/arviz.plot_posterior.html
https://python.arviz.org/en/latest/api/generated/arviz.plot_trace.html
4- loo_pit, This is again some manner of comparing your observed distribution to predictive samples. This will in the end produce a transformed distribution which tells you on (weighted) average how well your observation looks like the sampled predictions. weights come from the log probability of each sample. So I would say this is more akin to perhaps a distribution-like version of RMSE over multiple samples. Though, in this case. you want your transformed distribution to be as close to a uniform distribution as possible but using ecdf=True, you can get something which tells your model is good if the resulting statistics is close to 0. A link, again doing this for linear regression, is here:
Here it also actually presents a way to transform your posterior predictive samples for linear regression so posterior predictive plots make more sense.
5- If you want to test different priors and do model selection amongst them:
https://www.pymc.io/projects/docs/en/stable/learn/core_notebooks/model_comparison.html