# bayes\_R2 estimation
**URL:** https://discourse.mc-stan.org/t/bayes-r2-estimation/7420
**Category:** brms
**Created:** [January 24, 2019, 2:43pm UTC](https://discourse.mc-stan.org/t/bayes-r2-estimation/7420 "2019-01-24T14:43:12Z")
**Posts on this page:** 16
**Page:** 1
### Author:  [@helene](https://discourse.mc-stan.org/u/helene)
#### Post date: [January 24, 2019, 2:43pm UTC](https://discourse.mc-stan.org/t/bayes-r2-estimation/7420/1 "2019-01-24T14:43:12Z")
I am trying to understand the bayes\_R2 function to estimate the proportion of deviance explained for new data and I have been through the information from Gelman et al…
My concern was raised when analyzing my own data. I will illustrate my issue, using the model suggested in the vignette from Paul Bürkner on custom response distributions with the study case using the cbpp data (from the lme4 package).
The model that he suggests in a first step is the following:
fit1 \<- brm(incidence | trials(size) ~ period + (1|herd), data = cbpp, family = binomial())
To play around with the data, I wanted to perform the same model using bernouilli as the probability distribution instead of the Binomial distribution. I therefore performed a similar model but on a transformed version of the initial dataset in a way that each row is one individual observed:
fit2 \<- brm(incidence ~ period + (1|herd), data = cbpp\_individual, family = bernouilli())
the 2 models give the same estimates and errors on the parameters, which is obviously what I would have expected. However, the bayes\_R2 give completely different results for fit1=0.64 and fit2=0.10.
I cannot figure out what may explains this gap. I would be delighted to be enlightened on this point.
---
### Author:  [@avehtari](https://discourse.mc-stan.org/u/avehtari)
#### Post date: [January 24, 2019, 6:30pm UTC](https://discourse.mc-stan.org/t/bayes-r2-estimation/7420/2 "2019-01-24T18:30:50Z")
> [@helene](#):
>
> I am trying to understand the bayes\_R2 function
Can you specify which version? From the paper online appendix or brms package?
I also think that R2 is not valid for comparing different observation models (e.g Bernoulli vs. Binomial). I recommend using loo for that.
---
### Author:  [@helene](https://discourse.mc-stan.org/u/helene)
#### Post date: [January 24, 2019, 6:54pm UTC](https://discourse.mc-stan.org/t/bayes-r2-estimation/7420/3 "2019-01-24T18:54:48Z")
I use the bayes\_R2 from the bmrs 2.7.0 but I assume the function is built as written in the online paper from Gelman et al.
---
### Author:  [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)
#### Post date: [January 24, 2019, 8:26pm UTC](https://discourse.mc-stan.org/t/bayes-r2-estimation/7420/4 "2019-01-24T20:26:12Z")
Yes, Bayes\_R2 in brms is based on the online paper and should yield the same results.
I agree with Aki that using different observation models make the results of Bayes\_R2 incomparable.
You should use the same response vector in both models to make Bayes\_R2 comparisons meaningful.
---
### Author:  [@helene](https://discourse.mc-stan.org/u/helene)
#### Post date: [January 24, 2019, 8:42pm UTC](https://discourse.mc-stan.org/t/bayes-r2-estimation/7420/5 "2019-01-24T20:42:50Z")
my feeling in that the bayes\_R2 is valid for the model with a Bernouilli distribution but not with the Binomial. Indeed, doing the same kind of exploration on my own data, I found for bayesian models, a R2=39% with a binomial distribution, dropping to 7% with a bernouilli distribution. However, if I do the same model in a frequentist framework, I find about the same R2 as with a bernouilli, around 7%.
I am trying to understand the mathematics that may explain that difference. I am thinking loud now, but could not it be linked to the fact that, as you wrote Paul, “A drawback of the binomial model is that – after taking into account the linear predictor – its variance is fixed to Var(yi)=Tipi(1−pi). All variance exceeding this value cannot be not taken into account by the model.”?
---
### Author:  [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)
#### Post date: [January 25, 2019, 8:44am UTC](https://discourse.mc-stan.org/t/bayes-r2-estimation/7420/6 "2019-01-25T08:44:47Z")
One potential problem I see with binomial models and Bayes\_R2 is that responses are not on the same scale if `trials` is not constant across observations. I will have to take a closer look at this again, but for now I would say you should probably indeed prefer the bernoulli version out of the two.
---
### Author:  [@helene](https://discourse.mc-stan.org/u/helene)
#### Post date: [January 25, 2019, 8:50am UTC](https://discourse.mc-stan.org/t/bayes-r2-estimation/7420/7 "2019-01-25T08:50:39Z")
thank you both for your contributions, I will keep an eye open to follow the resolution of this problem.
---
### Author:  [@avehtari](https://discourse.mc-stan.org/u/avehtari)
#### Post date: [January 25, 2019, 9:45am UTC](https://discourse.mc-stan.org/t/bayes-r2-estimation/7420/8 "2019-01-25T09:45:13Z")
> [@paul.buerkner](#):
>
> Bayes\_R2 in brms is based on the online paper
The code in the online appendix does not support Binomial in all cases. It’s better to copy the code from @jonah’s more careful implementation in rstantools (I think Jonah did commit that in github in December).
---
### Author:  [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)
#### Post date: [January 25, 2019, 9:46am UTC](https://discourse.mc-stan.org/t/bayes-r2-estimation/7420/9 "2019-01-25T09:46:13Z")
Thanks! I will take a look!
---
### Author:  [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)
#### Post date: [January 25, 2019, 12:50pm UTC](https://discourse.mc-stan.org/t/bayes-r2-estimation/7420/10 "2019-01-25T12:50:14Z")
The default method in rstantools has no special handling for binomial models, but the method in rstanarm for stanreg objects has. However, this special handling is done in brms anyway (which is not visible in bayes\_R2.brmsfit as it happens on a deeper level), so that brms and rstanarm yield the same results of bayes\_R2 for binomial models.
---
### Author:  [@jonah](https://discourse.mc-stan.org/u/jonah)
#### Post date: [May 26, 2019, 4:15am UTC](https://discourse.mc-stan.org/t/bayes-r2-estimation/7420/11 "2019-05-26T04:15:20Z")
We’ll finally have the updated version in the next rstanarm release coming soon. It’s been implemented on GitHub for a while. The implementation (including for binomial models) is here:
>
and follows the method in the published version of our paper (not my personal GitHub repo):
Journal: [https://www.tandfonline.com/doi/abs/10.1080/00031305.2018.1549100](https://www.tandfonline.com/doi/abs/10.1080/00031305.2018.1549100)
Free on Andrew’s website: [http://www.stat.columbia.edu/~gelman/research/published/bayes\_R2\_v3.pdf](http://www.stat.columbia.edu/~gelman/research/published/bayes_R2_v3.pdf)
---
### Author:  [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)
#### Post date: [May 26, 2019, 7:08am UTC](https://discourse.mc-stan.org/t/bayes-r2-estimation/7420/12 "2019-05-26T07:08:44Z")
Thanks Jonah!
While looking at the stanreg code I noticed that there might be an inefficient pattern in it in the form of `mu_pred <- mu_pred %*% diag(trials)`. If we have lots of observations, `diag(trials)` will be huge and the matrix multiplication will be extremely slow. Here is a quick and dirty benchmark:
```
library(Matrix)
library(pryr)
N <- 10000
S <- 1000
set.seed(1234)
trials <- rpois(N, 1)
object_size(trials)
trials_diag <- diag(trials)
object_size(trials_diag)
trials_mat <- matrix(trials, nrow = S, ncol = N, byrow = TRUE)
object_size(trials_mat)
eta <- matrix(rnorm(S * N), nrow = S, ncol = N)
object_size(eta)
# current approach in bayes_R2.stanreg
system.time(
res1 <- eta %*% trials_diag
)
object_size(res1)
# used by brms
system.time(
res2 <- eta * trials_mat
)
object_size(res2)
# using the Matrix package
system.time(
res3 <- eta %*% Diagonal(x = trials)
)
object_size(res3)
# compare results
all.equal(res1, res2)
all.equal(res1, as.matrix(res3), check.attributes = FALSE)
```
---
### Author:  [@jonah](https://discourse.mc-stan.org/u/jonah)
#### Post date: [May 26, 2019, 2:51pm UTC](https://discourse.mc-stan.org/t/bayes-r2-estimation/7420/13 "2019-05-26T14:51:45Z")
Thanks Paul, that’s a good point! I’ll update the code.
---
### Author:  [@renatoamorais](https://discourse.mc-stan.org/u/renatoamorais)
#### Post date: [June 17, 2020, 1:40am UTC](https://discourse.mc-stan.org/t/bayes-r2-estimation/7420/14 "2020-06-17T01:40:05Z")
Any chance we could get an implementation of the Bayesian R2 for gamma and beta error distributions?
---
### Author:  [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)
#### Post date: [June 17, 2020, 6:16am UTC](https://discourse.mc-stan.org/t/bayes-r2-estimation/7420/15 "2020-06-17T06:16:12Z")
This should work already in brms. Any evidence to the contrary?
---
### Author:  [@renatoamorais](https://discourse.mc-stan.org/u/renatoamorais)
#### Post date: [June 17, 2020, 6:38am UTC](https://discourse.mc-stan.org/t/bayes-r2-estimation/7420/16 "2020-06-17T06:38:23Z")
Thanks for the indication, I got the error when trying the generic function on an rstanarm object.
It worked fine on a brms object, and now I see that the methods are independent… bayes\_R2.brmsfit does work indeed. Thank you very much!