What is the best way to define a multivariate regression in the context of hierarchical bayesian models?

@greghor you could try Bambi to build hierarchical models on top of PyMC3 using a formula interface like the one in R. You have several examples in the Examples section of the documentation.
Also, feel free to ping me if you have any questions about how it works.