Discrepancy in residual variance between VarCorr() and get_variance_r…

…esidual() Fixes #929
easystats · Nov 5, 2024 · c736fd0 · c736fd0
1 parent a4d1003
commit c736fd0
Show file tree

Hide file tree

Showing 2 changed files with 15 additions and 2 deletions.
diff --git a/R/get_sigma.R b/R/get_sigma.R
@@ -213,7 +213,7 @@ get_sigma <- function(x, ci = NULL, verbose = TRUE) {
       dat <- as.data.frame(x)
       sigma_column <- grep("sigma", colnames(dat), fixed = TRUE)
       if (length(sigma_column) == 1) {
-        mean(dat[[sigma_column]][1])
+        mean(dat[[sigma_column]])
       } else if (length(sigma_column)) {
         # if more than one sigma column,
         # there isn't a traditional sigma for the model
@@ -226,7 +226,6 @@ get_sigma <- function(x, ci = NULL, verbose = TRUE) {
       NULL
     }
   )
-
   # compute sigma manually ---------------
   if (is_empty_object(s)) {
     # default sigma ---------------

diff --git a/tests/testthat/test-brms.R b/tests/testthat/test-brms.R
@@ -921,3 +921,17 @@ test_that("get_variance works", {
   # make sure it's a matrix
   # expect_true(is.matrix(get_modelmatrix(null_model(mdl))))
 })
+
+
+# get variance
+test_that("get_variance aligns with get_sigma", {
+  skip_if_not_installed("lme4")
+  mdl <- brms::brm(mpg ~ hp + (1 | cyl), data = mtcars, seed = 123)
+  VC <- lme4 <- VarCorr(mdl)
+  out1 <- VC$residual__$sd[1, 1]^2 # Residual variance
+  out2 <- get_variance(mdl)$var.residual
+  out3 <- get_sigma(mdl)^2
+  expect_equal(out1, out2, tolerance = 1e-3, ignore_attr = TRUE)
+  expect_equal(out1, out3, tolerance = 1e-3, ignore_attr = TRUE)
+  expect_equal(out2, out3, tolerance = 1e-3, ignore_attr = TRUE)
+})