Somehow my google search didn’t yield many results for how to compute R2 for Turing models.

StatsBase.jl defines a few generic methods, but I was in particular interested in the so-called Bayesian R2 (Gelman et al., 2019), defined as:

Would that be easy to implement in Julia/Turing? Its R implementation seems pretty straightforward:

```
bayes_R2 <- function(fit) {
y <- rstanarm::get_y(fit)
ypred <- rstanarm::posterior_linpred(fit, transform = TRUE)
if (family(fit)$family == "binomial" && NCOL(y) == 2) {
trials <- rowSums(y)
y <- y[, 1]
ypred <- ypred %*% diag(trials)
}
e <- -1 * sweep(ypred, 2, y)
var_ypred <- apply(ypred, 1, var)
var_e <- apply(e, 1, var)
var_ypred / (var_ypred + var_e)
}
```