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