I am dealing with a demography problem over three dimensions, i.e., I am trying to reconstruct a number density across these three dimensions (they do not factorize). I am applying a hierarchical model with selection effects.

I have 2000 objects, each with their position in the three dimensions and additionally 5 nuisance parameters. The data are low count Poisson over ~800 bins, over which a low-level model is defined. One of the three parameters is a shift in related to the normalisation, another in a shift in bin offset, the third affects the model shape. The remaining nuisance parameters also affect the shape slightly.

The bin offset shift parameter has a informative prior from previous measurements, which is for some objects a delta function, for some Gaussian-ish, and for others unknown and thus uniform. The model shape parameter is poorly constrained in many low count datasets. All three are degenerate and frequently have bi-modal solutions (which are meaningful).

For some plots Figure 1 in [1501.02805] Obscuration-dependent evolution of Active Galactic Nuclei has a plot of the distribution of objects, here one posterior sample was drawn for each object assuming a flat population distribution. The three population dimensions are redshift (x-axis), luminosity (y-axis), and NH (inset histogram). But each object’s posterior is really **a cloud** in this 3d space. To get a feel for this, see Figure B.2 in [1402.0004] X-ray spectral modelling of the AGN obscuring region in the CDFS: Bayesian model selection and catalogue for a typical per-object posterior under flat priors, and Figure 5 for a fairly high-count data set.

My problem is that some of the data sets contain bad fits – something went wrong with the low-level model or the nuisance parameters, or the data are badly calibrated. Yet the objects individually are quite uninformative. Therefore posterior predictive checks in the data space of one object is not helpful to diagnose the population inference. Perhaps a projection could do it. Is there a principled way to construct one?

So I am looking for ideas how to visualize such hierarchical models, perhaps with some projections across parameters and/or objects, that are optimized to tell me:

- why the population distribution comes out with parameters A, and not with parameters B
- which objects are driving this result, i.e., where in the data and parameter space some population distribution shape emerges.

The approach I took in 2015 was to infer the population distribution non-parametrically, via a AR(1) inspired model over the 3d parameter space. This imposes some smoothness, while allowing one to visualise where the data wants to go. Nevertheless, it does not give insight which objects are pushing the shape of the population distribution, and why a different shape is disfavored.

Recently, I saw two papers on related topics.

- [2109.02712] Bayesian data selection – this looks quite related, but honestly, I did not quite understand how to apply this in practice.
- [2109.04702] Latent space projection predictive inference by @AlejandroCatalina @paul.buerkner @avehtari – it seems perhaps the projections are related, but the goal (and thus optimization target) seems different?

Pointers are very welcome.

PS: I guess a simple old-school approach is given parameter vector A and B, to color data points according to their likelihood ratio contribution.