Somehow I missed the initial question. ArviZ expects the array data to either have labeled dimensions with two of them being chain and draw or to follow the convention chain, draw, *shape if dimensions are not labeled.
To get the result of sample_posterior_predictive to generate arrays following ArviZ convention, you can use the keep_size argument.
There are also some examples of using hdi with multilevel models, even calculating hdi on groupby objects in the “radon” notebook: A Primer on Bayesian Methods for Multilevel Modeling — PyMC3 3.10.0 documentation (hdi is used several times so I am not linking to any specific section, search hdi within the page to see all the occurences)