easystats / insight

:crystal_ball: Easy access to model information for various model objects
https://easystats.github.io/insight/
GNU General Public License v3.0
404 stars 39 forks source link

Discrepancy in residual variance between VarCorr() and get_variance_residual() #929

Closed matschmitz closed 2 weeks ago

matschmitz commented 2 months ago

There seems to be a discrepancy between the residual variance calculated using brms::VarCorr() and insight::get_variance_residual() for a brms model. Apologies if I'm missing something obvious.

library(brms)    # v2.21.0
library(insight) # v0.20.4.3

mdl <- brm(mpg ~ hp + (1 | cyl), data = mtcars, seed = 123)
# Family: gaussian 
#  Links: mu = identity; sigma = identity 
# Formula: mpg ~ hp + (1 | cyl) 
#  Data: mtcars (Number of observations: 32) 
#  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#         total post-warmup draws = 4000
# 
# Multilevel Hyperparameters:
#   ~cyl (Number of levels: 3) 
#               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sd(Intercept)     4.75      2.55     1.33    11.21 1.00     1192     1719
# 
# Regression Coefficients:
#           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept    24.79      3.50    17.26    31.25 1.00     1454     1545
# hp           -0.03      0.02    -0.07    -0.00 1.00     2036     2145
# 
# Further Distributional Parameters:
#       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sigma     3.29      0.46     2.53     4.32 1.00     2045     2157

VC <- VarCorr(mdl)
# $cyl
# $cyl$sd
#           Estimate Est.Error     Q2.5    Q97.5
# Intercept 4.751544  2.549082 1.334246 11.21103
# 
# 
# $residual__
# $residual__$sd
# Estimate Est.Error     Q2.5    Q97.5
# 3.290936 0.4609515 2.531565 4.321534

VC$residual__$sd[1, 1]^2 # Residual variance
# [1] 10.83026

get_variance_residual(mdl)
# var.residual 
#      7.85012 
mariobecerra commented 2 months ago

Seems to be related to the new version. I'm running an earlier version (0.19.8) and I do get the same result in both cases (10.83026).

strengejacke commented 2 weeks ago
library(brms)    # v2.21.0
library(insight) # v0.20.4.3

mdl <- brm(mpg ~ hp + (1 | cyl), data = mtcars, seed = 123)

VC <- VarCorr(mdl)
VC$residual__$sd[1, 1]^2 # Residual variance
#> [1] 10.83026

get_variance(mdl)$var.residual
#> [1] 10.83026
#> attr(,"class")
#> [1] "insight_aux" "numeric"

get_sigma(mdl)^2
#> [1] 10.83026
#> attr(,"class")
#> [1] "insight_aux" "numeric"

get_variance_residual(mdl)
#> var.residual 
#>     10.83026

Created on 2024-11-05 with reprex v2.1.1