bayes_R2 computation does not match Bayesian R-squared paper eq (3) #1815
Description
"R-squared for Bayesian Regression Models" paper at https://doi.org/10.1080/00031305.2018.1549100 titles Equation (3) as "Bayesian R^2".
brms uses the following implementation
.bayes_R2 <- function(y, ypred, ...) {
e <- -1 * sweep(ypred, 2, y)
var_ypred <- matrixStats::rowVars(ypred)
var_e <- matrixStats::rowVars(e)
as.matrix(var_ypred / (var_ypred + var_e))
}
which matches Equation (2) titled as "alternative R^2". This could be also called as "residual R^2".
The code for the paper is at https://avehtari.github.io/bayes_R2/bayes_R2.html. Modifying the bayes_R2 code for brms would look like
bayes_R2 <- function(fit) {
mupred <- brms::posterior_epred(fit)
var_mupred <- apply(mupred, 1, var)
sigma2 <- as.matrix(fit, variable = c("sigma"))^2
var_mupred / (var_mupred + sigma2)
}
The important part of Bayes-R2 is that it doesn't depend on y.
The biggest difference with the current implementation and actual Bayes-R^2 can be seen when using sampling from the prior with sample_prior = "only". The above fixed bayes_R2 code correctly returns draws from the implied prior on R^2.
When sampling from the posterior there are differences for small n, and with big enough n, the differences get smaller (see examples at https://avehtari.github.io/bayes_R2/bayes_R2.html)
I suggest
- change bayes_r2.brmsfit to match the paper on normal and Bernoulli, for Bernoulli use Tjur's pseudo-variance
- for other families which do not yet have pseudo-variance, give a warning and fall back to residual-R^2
- for other families, add pseudo-variance when adding pseudo-variance support for R2D2 type priors
The benefit of using Bayes-R^2 as in the paper is that we can get the implied prior on R^2 for any priors on coefficients, and in case of R2D2 type priors the implied priors should match the beta prior used.