I need to compare the performance of my bayesian model against some others built by different methods. One of the metrics I hope to use is the R-squared. Libraries like sklearn allow you to calculate the R-squared using the actual and predicted values of a hold-out set; Sklearn R-Squared Docs

My model was fit using brms. I know brms has the bayes_R2() function to calculate the bayesian R2 using the fitted residuals, but is there an equivalent that allows the actual and predicted values of a second hold-out set as input?

I found Vehtari et al.'s notebook on calculating bayesian R2 https://avehtari.github.io/bayes_R2/bayes_R2.html, which I found helpful. I adapted their method to make use of actual and predicted values from a hold out setā¦

```
bayes_R2_res <- function(mod) {
# method obtained from https://avehtari.github.io/bayes_R2/bayes_R2.html
y <- mod$datahold$avg_sales. # this is the target in the hold-out set
ypred <- predict(mod, mod$datahold, allow_new_levels = T, summary=F) # predict avg_sales in holdout
e <- -1 * sweep(ypred, 2, y)
var_ypred <- apply(ypred, 1, var)
var_e <- apply(e, 1, var)
r2 <- var_ypred / (var_ypred + var_e)
return(as.data.frame(r2))
}
```

For reference, this produces a median R2 of ~0.38

However, when I compare this to the more typical R2, I get a vastly different result

```
y.predicted <- predict(mod, mod$datahold, allow_new_levels = T)[,1] # take the summary column for simplicity's sake
y.actual <- mod$datahold$avg_sales
errors <- (y.actual - y.predicted)
1 - sum(errors^2)/sum((y.actual-mean(y.actual))^2)
```

This returns a value of -79.26 (?!)

taking the square of the correlation coefficient also gives a very different estimate;

```
cor(y.actual, y.predicted)^2 # returns 0.0004
```

I thought I understood R-squared as a basic metric. I must be making a simple error somewhere but I feel like Iām taking crazy pills!

Appreciate it if anyone can either a) tell me where Iām going wrong, or b) point me to a function that will allow me to calculate R-squared from actual and predicted values, or using a brms model object and a new data set.

Thanks