# bayes\_R2 and conditioning on random effects

**URL:** <https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461>\
**Category:** Modeling\
**Tags:** rstanarm\
**Created:** [January 26, 2021, 11:41pm UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461 "2021-01-26T23:41:46Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![strengejacke](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/strengejacke/32/559_2.png) [@strengejacke](https://discourse.mc-stan.org/u/strengejacke)\
**Post date:** [January 26, 2021, 11:41pm UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/1 "2021-01-26T23:41:46Z")

</div>

I have a question regarding `bayes_R2()`. From the docs, I understand that `re.form` can be used to condition or refrain from conditioning on random effects in mixed models:

> The default, NULL, indicates that all estimated group-level parameters are conditioned on.

However, I don’t quite understand the output. Similar to a frequentist framework, I would expect the conditional R2, i.e. the R2 for the model taking random effects into account, to be higher than the “marginal” R2 that only considers fixed effects. In other word, I would expect the point estimate for `bayes_R2(m, re.form = NULL)` to be higher than for `bayes_R2(m, re.form = NA)`.

```no-highlight
library(rstanarm)
m <- stan_lmer(Sepal.Length ~ Petal.Length + (1 | Species), data = iris, refresh = 0, seed = 333)

# The default, NULL, indicates that all estimated 
# group-level parameters are conditioned on.
mean(bayes_R2(m, re.form = NULL))
#> [1] 0.8278072

# "marginal" R2?
mean(bayes_R2(m, re.form = NA))
#> [1] 0.9515886

```

Is there anything I am missing?

---

<div class="post-metadata">

**Author:** ![strengejacke](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/strengejacke/32/559_2.png) [@strengejacke](https://discourse.mc-stan.org/u/strengejacke)\
**Post date:** [January 29, 2021, 10:21am UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/2 "2021-01-29T10:21:43Z")

</div>

_(bumping)_

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [January 30, 2021, 12:02pm UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/3 "2021-01-30T12:02:03Z")

</div>

I know about R2, but I’m clueless what re.form does. Can you make scatter plots of Sepal.Length vs. the predictions with “condition or refrain from conditioning on random effects”? That would help to see the the different behaviors.

---

<div class="post-metadata">

**Author:** ![DominiqueMakowski](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/dominiquemakowski/32/2315_2.png) [@DominiqueMakowski](https://discourse.mc-stan.org/u/DominiqueMakowski)\
**Post date:** [January 31, 2021, 1:59am UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/4 "2021-01-31T01:59:41Z")

</div>

Hey! Here’s the output :)

```
library(rstanarm)

model <- rstanarm::stan_lmer(Sepal.Length ~ Petal.Length + (1|Species), data=iris, refresh=0)

# The default (includes random factor)
ypred <- rstanarm::posterior_predict(model, re.form = NULL)
ypred <- sapply(as.data.frame(ypred), median)
plot(x=iris$Sepal.Length, y=ypred)

```

![](https://i.imgur.com/eNYlDs0.png)

```
# Do not take into account the Species
ypred <- rstanarm::posterior_predict(model, re.form = NA)
ypred <- sapply(as.data.frame(ypred), median)
plot(x=iris$Sepal.Length, y=ypred, col="red")

```

![](https://i.imgur.com/Nnu7IhM.png)

Created on 2021-01-31 by the [reprex package](https://reprex.tidyverse.org) (v0.3.0)

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [January 31, 2021, 5:08pm UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/5 "2021-01-31T17:08:08Z")

</div>

I still don’t understand what re.form does, but the predictions in the second plot are surprising as group of predictions go below 4, while the target value in the data doesn’t go below 4.4 (and then ypred has unexpected high variance). Can you explain that?

---

<div class="post-metadata">

**Author:** ![DominiqueMakowski](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/dominiquemakowski/32/2315_2.png) [@DominiqueMakowski](https://discourse.mc-stan.org/u/DominiqueMakowski)\
**Post date:** [February 1, 2021, 12:29am UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/6 "2021-02-01T00:29:54Z")

</div>

I’m not exactly sure of what it is _supposed_ to do, but at least it seems to give similar results to lme4:

```r
library(lme4)

model <- lmer(Sepal.Length ~ Petal.Length + (1|Species), data=iris)

ypred <- predict(model, re.form = NULL)
plot(x=iris$Sepal.Length, y=ypred)

```

![](https://i.imgur.com/7ZqO6MC.png)

```r
ypred <- predict(model, re.form = NA)
plot(x=iris$Sepal.Length, y=ypred, col="red")

```

![](https://i.imgur.com/xfA6Phj.png)

Created on 2021-02-01 by the [reprex package](https://reprex.tidyverse.org) (v0.3.0)

The [original arg explanation](https://www.rdocumentation.org/packages/lme4/versions/1.1-26/topics/predict.merMod) says that:

> re.form specify which random effects to condition on when predicting. If `NULL` , include all random effects; if `NA` or `~0` , include no random effects.

If I’m not wrong, the predictions are mostly driven by the population effects than the individual levels (see [this thread](https://stats.stackexchange.com/questions/262277/why-would-you-predict-from-a-mixed-effect-model-without-including-random-effects)).

---

<div class="post-metadata">

**Author:** ![strengejacke](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/strengejacke/32/559_2.png) [@strengejacke](https://discourse.mc-stan.org/u/strengejacke)\
**Post date:** [February 1, 2021, 8:07pm UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/7 "2021-02-01T20:07:53Z")

</div>

In terms of the R2, because a mixed model has different levels of variance, I would expect a higher R2 when conditioning on random effects - that R2 takes _all_ sources of uncertainty into account and thus has a “better explained variance”. `bayes_R2(m, re.form = NA)` should only explain the variance from the fixed effects, and hence should be lower. This is what you also see from functions that compute the R2 for mixed models:

```
library(lme4)
#> Loading required package: Matrix
m <- lmer(Sepal.Length ~ Petal.Length + (1 | Species), data = iris)
performance::r2(m)
#> # R2 for Mixed Models
#> 
#> Conditional R2: 0.969
#> Marginal R2: 0.658

```

---

<div class="post-metadata">

**Author:** ![DominiqueMakowski](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/dominiquemakowski/32/2315_2.png) [@DominiqueMakowski](https://discourse.mc-stan.org/u/DominiqueMakowski)\
**Post date:** [February 22, 2021, 12:19am UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/8 "2021-02-22T00:19:30Z")

</div>

Maybe @paul.buerkner who worked with mixed models might have some idea?

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [February 22, 2021, 2:28pm UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/9 "2021-02-22T14:28:59Z")

</div>

re.form = NA sets all random effects to 0, that is, is still a “conditional” R2 but not condition on the estimates random effects but conditioned on an “average” random effects level. This is _not_, in general a marginal formulation because marginal would mean integrated over the random effects distritbution. Only in linear models, that conditional-on-zero and marginal formulations fall together.

---

<div class="post-metadata">

**Author:** ![DominiqueMakowski](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/dominiquemakowski/32/2315_2.png) [@DominiqueMakowski](https://discourse.mc-stan.org/u/DominiqueMakowski)\
**Post date:** [February 26, 2021, 12:46am UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/10 "2021-02-26T00:46:23Z")

</div>

```
library(rstanarm)
m <- stan_lmer(Sepal.Length ~ Petal.Length + (1 | Species), data = iris, refresh = 0, seed = 333)

mean(bayes_R2(m, re.form = NA))
#> [1] 0.9515886

```

Thanks Paul! If I understand you correctly, the above R2 is conditioned on an average random effect level. But how should we interpret it, in comparison with `mean(bayes_R2(m, re.form = NULL)) = 0.82`?

In the frequentist marginal/conditional R2 framework, the two R2s can be interpreted as the explanatory power of the full model and that of only the fixed effects, from which one can derive the part of the random effects. What is the correct interpretation for the two Bayesian R2s?

---

<div class="post-metadata">

**Author:** ![mattansb](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mattansb/32/20280_2.png) [@mattansb](https://discourse.mc-stan.org/u/mattansb)\
**Post date:** [March 21, 2021, 6:32am UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/11 "2021-03-21T06:32:35Z")

</div>

@paul.buerkner sorry to bump this again, but it’s still not clear why in such a case the Bayesian (stan\_lmer) and frequentist (lmer) R2 are so different. Specifically, why `re.form = NULL` \< `re.form = NA`?

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [March 22, 2021, 9:54am UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/12 "2021-03-22T09:54:33Z")

</div>

The bayes\_R2 presumably uses a very different formula to come up with an R2 estimate. That is, bayes\_R2 and performance::r2 yield incomparable results regardless of what one does with re.form. I have no detailed knowledge of performance::r2 though. If you want to rule out a bug in bayes\_R2 for stan\_lmer, please also try the same thing with brms.

---

<div class="post-metadata">

**Author:** ![strengejacke](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/strengejacke/32/559_2.png) [@strengejacke](https://discourse.mc-stan.org/u/strengejacke)\
**Post date:** [March 22, 2021, 3:50pm UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/13 "2021-03-22T15:50:56Z")

</div>

> [@paul.buerkner](#):
>
> That is, bayes\_R2 and performance::r2 yield incomparable results regardless of what one does with re.form.

No, actually results should be identical, because internally, `performance::r2()` just calls `bayes_R2()`, and for _brms_, results are expected and go into the same direction as for lme4 (i.e. conditional R2 \> marginal R2)

```r
library(brms)
library(performance)

set.seed(123)
model <- brm(Sepal.Length ~ Petal.Length + (1 | Species), data = iris)

bayes_R2(model, robust = FALSE, re.form = NULL)
#> Estimate Est.Error Q2.5 Q97.5
#> R2 0.8327651 0.01103648 0.807663 0.8499841
bayes_R2(model, robust = FALSE, re.form = NA)
#> Estimate Est.Error Q2.5 Q97.5
#> R2 0.7370286 0.01437964 0.7105168 0.7663097

r2(model, robust = FALSE)
#> # Bayesian R2 with Standard Error
#> 
#> Conditional R2: 0.833 (89% CI [0.817, 0.849])
#> Marginal R2: 0.737 (89% CI [0.713, 0.759])

bayes_R2(model, robust = TRUE, re.form = NULL)
#> Estimate Est.Error Q2.5 Q97.5
#> R2 0.834307 0.01035303 0.807663 0.8499841
bayes_R2(model, robust = TRUE, re.form = NA)
#> Estimate Est.Error Q2.5 Q97.5
#> R2 0.7365296 0.01460497 0.7105168 0.7663097

r2(model, robust = TRUE)
#> # Bayesian R2 with Standard Error
#> 
#> Conditional R2: 0.834 (89% CI [0.817, 0.849])
#> Marginal R2: 0.737 (89% CI [0.713, 0.759])

```

however, for _rstanarm_, results are a bit counter-intuitive:

```r
library(rstanarm)
library(performance)

set.seed(123)
model <- stan_lmer(Sepal.Length ~ Petal.Length + (1 | Species), data = iris, refresh = 0)

mean(bayes_R2(model, robust = FALSE, re.form = NULL))
#> [1] 0.828523
mean(bayes_R2(model, robust = FALSE, re.form = NA))
#> [1] 0.951746

r2(model, robust = FALSE)
#> # Bayesian R2 with Standard Error
#> 
#> Conditional R2: 0.829 (89% CI [0.799, 0.863])
#> Marginal R2: 0.952 (89% CI [0.938, 0.966])

```

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [March 22, 2021, 5:03pm UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/14 "2021-03-22T17:03:08Z")

</div>

OK. thanks for clarifying. in this case I would say it is reasonable to assume there is a bug in rstanarm.

---

<div class="post-metadata">

**Author:** ![mattansb](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mattansb/32/20280_2.png) [@mattansb](https://discourse.mc-stan.org/u/mattansb)\
**Post date:** [March 22, 2021, 8:27pm UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/15 "2021-03-22T20:27:37Z")

</div>

It seems like the “bug” (see bellow) is specifically in `rstanarm::bayes_R2()`, as it is the only place where `re.form = NA > re.form = NULL`:

```

y <- iris$Sepal.Length

## lme4::lmer ----

m <- lme4::lmer(Sepal.Length ~ Petal.Length + (1 | Species), data = iris)

p1 <- predict(m, re.form = NULL)
p2 <- predict(m, re.form = NA)

c(
  "with random" = cor(p1, y),
  "no random" = cor(p2, y)
)
#> with random no random 
#> 0.9146892 0.8717538

## nlme::lme ----

m <- nlme::lme(Sepal.Length ~ Petal.Length, random = ~ 1 | Species, data = iris)

p1 <- predict(m, level = 1)
p2 <- predict(m, level = 0)

c(
  "with random" = cor(p1, y),
  "no random" = cor(p2, y)
)
#> with random no random 
#> 0.9146892 0.8717538

## brms::brm -----

m <- brms::brm(Sepal.Length ~ Petal.Length + (1 | Species), data = iris,
               chains = 1, refresh = 0)

p1 <- rstantools::posterior_epred(m, re.form = NULL)
p1 <- apply(p1, 2, mean)
p2 <- rstantools::posterior_epred(m, re.form = NA)
p2 <- apply(p2, 2, mean)

c(
  "with random" = cor(p1, y),
  "no random" = cor(p2, y)
)
#> with random no random 
#> 0.9144018 0.8718073

c(
  "with random" = median(rstantools::bayes_R2(m, re.form = NULL)),
  "no random" = median(rstantools::bayes_R2(m, re.form = NA))
)
#> with random no random 
#> 0.820947 0.724174

## rstanarm::stan_(g)lmer -----

m <- rstanarm::stan_glmer(Sepal.Length ~ Petal.Length + (1 | Species), data = iris,
                          chains = 1, refresh = 0)

p1 <- rstantools::posterior_epred(m, re.form = NULL)
p1 <- apply(p1, 2, mean)
p2 <- rstantools::posterior_epred(m, re.form = NA)
p2 <- apply(p2, 2, mean)

c(
  "with random" = cor(p1, y),
  "no random" = cor(p2, y)
)
#> with random no random 
#> 0.9145250 0.8717538

c(
  "with random" = median(rstantools::bayes_R2(m, re.form = NULL)),
  "no random" = median(rstantools::bayes_R2(m, re.form = NA))
)
#> with random no random 
#> 0.8276229 0.9516522

```

## The cause!

This seems to be because while `brms::bayes_R2.brmsfit()` estimates R2 as:  
`Var(pred) / (Var(pred) + Var(pred - y))`

`rstanarm:::bayes_R2.stanreg()` estimates R2 as:  
`Var(pred) / (Var(pred) + sigma^2)`

If we use the sigma method for `lme4::lmer`, we get the same results as `rstanarm:::bayes_R2.stanreg()`!

```
## lme4::lmer ----

m <- lme4::lmer(Sepal.Length ~ Petal.Length + (1 | Species), data = iris)

p1 <- predict(m, re.form = NULL)
p2 <- predict(m, re.form = NA)

c(
  "with random" = var(p1) / (var(p1) + sigma(m)^2),
  "no random" = var(p2) / (var(p2) + sigma(m)^2)
)
#> with random no random 
#> 0.8317347 0.9533562

```

* * *

Where do I collect my trophy? 😁

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [March 23, 2021, 2:51pm UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/16 "2021-03-23T14:51:43Z")

</div>

Excellent that you found the cause for the difference

> [@mattansb](#):
>
> This seems to be because while `brms::bayes_R2.brmsfit()` estimates R2 as:  
> `Var(pred) / (Var(pred) + Var(pred - y))`

This is the “residual based R^2” in [Bayesian R2 and LOO-R2](https://avehtari.github.io/bayes_R2/bayes_R2.html) and [R-squared for Bayesian regression models](http://www.stat.columbia.edu/~gelman/research/unpublished/bayes_R2_v3.pdf).

> [@mattansb](#):
>
> `rstanarm:::bayes_R2.stanreg()` estimates R2 as:  
> `Var(pred) / (Var(pred) + sigma^2)`

This is the “model based R^2” in [Bayesian R2 and LOO-R2](https://avehtari.github.io/bayes_R2/bayes_R2.html) and [R-squared for Bayesian regression models](http://www.stat.columbia.edu/~gelman/research/unpublished/bayes_R2_v3.pdf).

> [@mattansb](#):
>
> Where do I collect my trophy? 😁

Here [Sign in to GitHub · GitHub](https://github.com/paul-buerkner/brms/issues/new) 😉

---

<div class="post-metadata">

**Author:** ![mattansb](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mattansb/32/20280_2.png) [@mattansb](https://discourse.mc-stan.org/u/mattansb)\
**Post date:** [March 23, 2021, 4:37pm UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/17 "2021-03-23T16:37:30Z")

</div>

@avehtari But is this an issue with `brms` though?

- `rstantools::bayes_r2.default()` uses **residual** based R^2
- `brms:::bayes_R2.brmsfit()` uses **residual** based R^2
- `rstanarm:::bayes_R2.stanreg()` uses _ **model** _ based R^2

(Also specifically for mixed models, wouldn’t the model-based R^2 need to also take into account other – random – sources of variance other than `sigma`? I think this is why the results of this method seem to us a bit wonky / don’t lineup with Nakagawa’s R^2 method for mixed models, which is also model based).

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [March 23, 2021, 5:05pm UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/18 "2021-03-23T17:05:50Z")

</div>

I would favor model based R^2, so based on that it would also be rstantools issue. Unless I missed something and you have justification why model based doesn’t work in some case?

---

<div class="post-metadata">

**Author:** ![mattansb](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mattansb/32/20280_2.png) [@mattansb](https://discourse.mc-stan.org/u/mattansb)\
**Post date:** [March 24, 2021, 6:10am UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/19 "2021-03-24T06:10:38Z")

</div>

Model-based is fine, but it should contain all the relevant sources of variance. In a mixed model these are the variances from the fixed effects from the random effects and the level 1 variance (\sigma^2).

So for mixed effects, I would expect a model-based R^2 to look like this:

\frac{V\_{predicted}}{V\_{fixed} + V\_{random} + \sigma^2}

Then if V\_{predicted} accounts from random effects, this R^2 gives the total variance accounted for the my model, and it V\_{predicted} does not account for random effects, it is the R^2 accounted for only by the fixed effects.

Whereas currently the model-based R^2 is:

\frac{V\_{predicted}}{V\_{predicted} + \sigma^2}

Which is fine (I think) for when the prediction accounts for the random variance (`re.form = NULL`), but I’m not sure how to interpret this method when the prediction does _not_ account for random effects (`re.form = NA` or fixed in some other way) - it is the R^2 in an assumed world where random effects are 0? Or some sort of partial-R^2 (with random effects partialled out)?

---

<div class="post-metadata">

**Author:** ![fweber144](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/fweber144/32/15592_2.png) [@fweber144](https://discourse.mc-stan.org/u/fweber144)\
**Post date:** [March 24, 2021, 4:14pm UTC](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461/20 "2021-03-24T16:14:30Z")

</div>

Also, wouldn’t \sigma^2 have to be re-estimated in case of `re.form = NA`, i.e. when the random effects are set to zero? I’m not sure if this is currently done in the code.

[Next page](https://discourse.mc-stan.org/t/bayes-r2-and-conditioning-on-random-effects/20461.md?page=2)