easystats / bayestestR

:ghost: Utilities for analyzing Bayesian models and posterior distributions
https://easystats.github.io/bayestestR/
GNU General Public License v3.0
571 stars 55 forks source link

Convert RE variances from cov to cor/sd #536

Open strengejacke opened 2 years ago

strengejacke commented 2 years ago

you see that correlations are not bounded between [-1, 1]. This is something to be fixed in bayestestR, I think, but do you have any idea how to transform the posterior in order to have a proper "point estimate" here?

I think the values from m1 are covariances and variances, not correlations and standard deviation. So you would have to:

library(rstanarm)
library(lme4)
data(cake)
set.seed(123)

m1 <- stan_lmer(angle ~ temp + (temp | replicate),
                data = cake,
                chains = 1)

bayestestR::describe_posterior(m1, effects = "all")
#> ...
#> 
#> # Random effects SD/Cor: replicate
#> 
#> Parameter          |   Median |         95% CI |     pd |          ROPE | % in ROPE |  Rhat |    ESS
#> ----------------------------------------------------------------------------------------------------
#> (Intercept)        | 2.90e-03 | [ 0.00,  0.70] |   100% | [-0.82, 0.82] |      100% | 0.999 | 612.00
#> temp ~ (Intercept) | 1.86e-04 | [-0.01,  0.01] | 59.20% | [-0.82, 0.82] |      100% | 0.999 | 395.00
#> temp               | 1.12e-03 | [ 0.00,  0.00] |   100% | [-0.82, 0.82] |      100% | 1.000 | 325.00
#> 
#> ...

posteriors <- as.data.frame(m1)

rho <- with(posteriors, {
  `Sigma[replicate:temp,(Intercept)]` / 
    (sqrt(`Sigma[replicate:temp,temp]`) * sqrt(`Sigma[replicate:(Intercept),(Intercept)]`))
})

# Cor
bayestestR::describe_posterior(rho)
#> Summary of Posterior Distribution
#> 
#> Parameter | Median |        95% CI |     pd |          ROPE | % in ROPE
#> -----------------------------------------------------------------------
#> Posterior |   0.19 | [-0.92, 0.95] | 59.20% | [-0.10, 0.10] |     9.05%

# Slopes
a <- with(posteriors, sqrt(`Sigma[replicate:temp,temp]`))

bayestestR::describe_posterior(a)
#> Summary of Posterior Distribution
#> 
#> Parameter | Median |       95% CI |   pd |          ROPE | % in ROPE
#> --------------------------------------------------------------------
#> Posterior |   0.03 | [0.02, 0.05] | 100% | [-0.10, 0.10] |      100%

# Compare to:
lme4::lmer(angle ~ temp + (temp | replicate),
           data = cake) |> 
  lme4::VarCorr()
#>  Groups    Name        Std.Dev. Corr 
#>  replicate (Intercept) 2.604527      
#>            temp        0.027073 0.134
#>  Residual              4.825927

Originally posted by @mattansb in https://github.com/easystats/parameters/issues/720#issuecomment-1143318379

strengejacke commented 2 weeks ago

This only affects stanreg, not brms. And we may have to revise the wording in insight::clean_parameters(), and change "SD/COR" to "var/cov" for stanreg-models.