# Metric as straightforward as R^2 for Bayesian models

So, the beauty of the $R^2$ in linear models or the deviance-based pseudo-$R^2$ from GLMs is their intuitive interpretation for non-specialists. There's also some nice developments on this front for GLMMs (see Nakagawa & Schielzeth's (2013) methods and extensions from folk like @jslefche's https://github.com/jslefche/piecewiseSEM package).

But - when it comes to funkier models, particularly those fit using Bayesian methods, I'm at a bit more of a loss. While the LOO metrics make sense if you do a nice deep dive into what they mean, for those of us who are collaborating with non-specialists in particular, is there something out there - a test statistic or somesuch - that is as straightforward to communicate as the $R^2$? I've been working a bit with quantile residuals for assessing GLMMs and more funky models a lot lately - https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html - and wonder if there might be something there?

But, how do we make fit metrics for complex Bayesian models that are interpretable to mere mortals?

See also Andrew Gelman's paper literally titled "R-squared for Bayesian regression models":

http://www.stat.columbia.edu/~gelman/research/unpublished/bayes_R2.pdf

In that paper, he proposes $R^2 = Var(\hat{y}) / (Var(\hat{y}) + Var(e))$ which has a similar interpretation to traditional $R^2$. However, since this is bayesian, that value is now a r.v. and needs to be summarised, eg by posterior median.

He also provides a small R function for calculating it.

To check model adequacy, hierarchical Bayesian models are usually evaluated exactly like in DHARMa - you simulated from the fitted model and calculate the quantile residuals. The only difference to DHARMa is that you vary parameters as well while doing the simulations.

The approach is explained in many textbooks. Keywords are "posterior predictive simulations", "posterior model checks", or "Bayesian p-values". A few notes

1. As for mixed models, the question arises how many levels of the hierarchy your re-simulate. Especially state-space population models look great when re-simulating only the final observation model, but if you re-simulate the process model as well, it often looks very very different -> good reason to be skeptical when people only show the final fit of the SSMs, this can look great even if the model is horribly wrong.

2. The distributional theory for these residuals is not quite clear for complex models - for GLMMs, all simulations I did so far showed that simulated quantile residuals are with good approximation flat, but this may not be so for any complex model structure. In doubt, run a test on a properly specified model.

3. DHARMa has a function to input posterior predictive simulations (createDHARMa), so that you can use all tools of DHARMa together with hierarchical Bayesian models as well

However, all this is more a substitute for standard residual diagnostics, not R2, which is a measure of how well you explain the target variable.

For R2, you can simply calculated R2 or RMSE over the posterior predictive distributions, whatever you like best, however, with all the caveats that arise for mixed models, i.e. you have to consider what you mean by "explain", does a random effect "explain" the data or not - this is again particularly apparent in state-space models, where you can get your fitted states very close to the data, even if the model has basically no clue why the data moves this way.