bayes_R2 and conditioning on random effects

Modeling
#1

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).

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?

2 Likes
#2

(bumping)

#3

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.

1 Like
#4

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)


# 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")

Created on 2021-01-31 by the reprex package (v0.3.0)

#5

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?

#6

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

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)


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

Created on 2021-02-01 by the reprex package (v0.3.0)

The original arg explanation 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).

#7

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
1 Like
#8

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

#9

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.

#10 
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?

#11

@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?

#12

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.

#13

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)

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:

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])
#14

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

#15

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? 😁

5 Likes
#16

Excellent that you found the cause for the difference

This is the “residual based R^2” in Bayesian R2 and LOO-R2 and R-squared for Bayesian regression models.

This is the “model based R^2” in Bayesian R2 and LOO-R2 and R-squared for Bayesian regression models.

Here https://github.com/paul-buerkner/brms/issues/new 😉

#17

@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).

#18

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?

#19

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)?

1 Like
#20

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.