Bayesian RMSE?

Modeling
1

Hi, is there some way how to get posterior draws for a model RMSE? Similar to how bayes_R2() works? I’ve tried looking around, but all I can find is a single point estimate using the rmse() function from the performance package. I also tried customizing the bayes_R2() function to compute RMSE draws, however, I get a right-skewed distribution with some RMSE values being extremely high (my training data RMSE was ~0.65, some of the RMSE draws were as large as ~3, this is with standardized response + predictor variables). I’m not sure if my code is wrong or if this is an expected/desirable behaviour & I’m just thinking about the problem wrong.

Thank you for your help!

2 Likes
2

This is pretty easy to calculate if you have the predictive posterior distribution (usually function posterior_predict). That will give you a matrix of posterior draws for the outcome/response y_{rep}, and then you can calculate RMSE relative to the original y by iterating over the rows of the posterior predictive matrix for each posterior draw d:

\hat{RMSE} = \sqrt{\frac{\sum_{n=1}^N (y_{repn}-y_n)^2}{N}}

Then you’ve got RMSE with uncertainty based on the number of posterior draws.

3 Likes
3

Can you post your customized code so we can check if it’s correct?

2 Likes
4

The code I used was from the appendix of Gelman, Goodrich, Gabry, & Ali (2017). I commented out the lines that I didn’t use:

bayes_rmse <- function(fit) {
   y <- get_y(fit)
   ypred <- 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)

se <- apply(e, 1, `^`, 2)
mse <- apply(se, 1, mean)
rmse <- sqrt(mse)

rmse

}
5

For the purpose of calculating the predictive errors e=y-ypred for RMSE you’ll need to use posterior_predict instead of posterior_linpred for ypred. The reason is that posterior_predict will give you draws from the posterior predictive distribution of the outcome (instead of just draws from the posterior distribution of the (transformed) linear predictor a + bx). Using “ypred” with posterior_linpred was probably a poor choice of name on our part. Sorry for the confusion!

3 Likes
6

Where did you find that? I’m asking as there is better newer version Gelman, Goodrich, Gabry, & Vehtari (2018) https://doi.org/10.1080/00031305.2018.1549100
with better and more code at
Bayesian R2 and LOO-R2, and it would be good to know if we still have somewhere link to the old version.

1 Like
7

I got it from here, one of the first articles that pops up when I google “bayes R2”: http://www.stat.columbia.edu/~gelman/research/unpublished/bayes_R2.pdf

I’ll check out the newer version, thanks!

8

Finally I got something working which looks approximately like what I expected:

bayes_rmse <- function(object, newdata, y) {
  
  ypred <- posterior_linpred(object, xnew = newdata)
  e <- apply(ypred, 1, `-`, y)
  
  se <- apply(e, 1, `^`, 2)
  mse <- apply(se, 1, mean)
  rmse <- sqrt(mse)
  
  rmse
  
}

The distribution of RMSEs looks a bit more reasonable. Before, it was heavily right-skewed with some values going up to 3, now it’s more similar to the distribution of R2’s I get with bayes_R2(). Weirdly enough, my “observed” RMSE (~ 0.658) using mean posterior point estimates is just outside of the range of the distribution.

r2rmses
r2rmses707×550 5.28 KB

However, one thing I noticed is that the range of the RMSE’s is quite restricted. Does that issue come from using posterior_linpred()? I tried using posterior_predict() as @jonah suggested, however, then I get a symmetrical distribution centered at 0.96 (with range ~ 0.88 - 1.06), which does not match the “observed” RMSE or how well the model seems to predict test data.

Does it make sense to get a posterior draws for RMSE in the first place? Or am I on a wild goose chase here?

9

Thanks. I’ve asked @andrewgelman to replace that with the latest version, but it seems I need to remind him that the old version is still causing confusion.

10

Your code still looks quite different compared to the new code, which make me suspect you didn’t yet look at the latest version of the paper and appendix.

Look at the online appendix Bayesian R2 and LOO-R2. See the difference between residual based and model based R^2. Recognise which term corresponds to squared error.
See the code for “Bayes-R2 function using residuals” and “Bayes-R2 function using modelled (approximate) residual variance”. In the latest version we prefer the modelled version. You may be surprised to see how simple the final solution for RMSE is. I don’t now copy the material from the online appendix to here, but let me know if you don’t see the solution and I can add more explicit comment and add more detailed explanation here.

It is very popular approach.

11

Oh…would it just be sigma from the model?

1 Like
12

Yes! I’m happy you could see that. Inspired by this discussion, I will also add explicit comment about that next time I’m writing about R^2 and RMSE.

1 Like
13

Cheers Aki, thank you for gently nudging me so that I could blunder onto the answer on my own!

2 Likes