init

Author: DominiqueMakowski

NAMESPACE

@@ -118,6 +118,8 @@ S3method(check_predictions,brmsfit)
 S3method(check_predictions,default)
 S3method(check_predictions,lme)
 S3method(check_predictions,stanreg)
+S3method(check_reliability,default)
+S3method(check_reliability,estimate_grouplevel)
 S3method(check_residuals,DHARMa)
 S3method(check_residuals,default)
 S3method(check_residuals,performance_simres)
@@ -570,6 +572,7 @@ export(check_normality)
 export(check_outliers)
 export(check_overdispersion)
 export(check_predictions)
+export(check_reliability)
 export(check_residuals)
 export(check_singularity)
 export(check_sphericity)

R/check_reliability.R

@@ -0,0 +1,95 @@
+#' Variability-Over-Uncertainty Ratio (D-vour) for Random Effects Reliability
+#'
+#' @description TODO: Add description.
+#'
+#' @param x A model object.
+#' @param ... Currently not used.
+#'
+#
+#'
+#' @details TODO: Add details.
+#'
+#' @references TODO.
+#'
+#' @family functions to check model assumptions and and assess model quality
+#'
+#' @examplesIf require("lme4")
+#' # Add groups to the data
+#' data <- iris
+#' data$Place <- as.factor(rep(c("P1", "P2", "P3", "P4", "P5", "P6"), each = 25))
+#'
+#' # lme4
+#' m <- lme4::lmer(Sepal.Width ~ Petal.Width + (Petal.Width | Place), data = data)
+#' check_reliability(m)
+#'
+#' @export
+check_reliability <- function(x, ...) {
+  UseMethod("check_reliability")
+}
+
+
+
+
+#' @export
+check_reliability.default <- function(x, ...) {
+  check_reliability(modelbased::estimate_grouplevel(x, ...), ...)
+}
+
+
+#' @export
+check_reliability.estimate_grouplevel <- function(x, ...) {
+
+  coefname <- attributes(x)$coef_name
+  dispname <- names(x)[grep("SE|SD|MAD", names(x))]
+
+  # Sanity checks
+  if(length(unique(x$Level)) <= 3) {
+    warning(paste0("The number of random levels (N = ",
+                   length(unique(x$Level)),
+                   ") might be too low to reliably estimate the variability."))
+  }
+
+  if(length(dispname) == 0) {
+    stop(paste0("This function requires an index of variability of each random ",
+    "effect (e.g., SE) but none was found. Try running check_reliability() on the",
+    " output of modelbased::estimate_grouplevel(model), and make sure the latter ",
+    "returns a table with an index of dispersion."))
+  }
+
+  if(length(dispname) > 1) {
+    warning(paste0("Multiple indices of variability were found (",
+                   paste(dispname, collapse = ", "),
+                   "). Using the first one."))
+    dispname <- dispname[1]
+  }
+
+
+  # Compute reliability
+  if (!"Component" %in% names(x)) x$Component <- "TEMP"
+
+  reliability <- data.frame()
+  for(c in unique(x$Component)) {
+    for(g in unique(x$Group)) {
+      for(p in unique(x$Parameter)) {
+        d <- x[x$Component == c & x$Group == g & x$Parameter == p,]
+        rez <- data.frame(
+          Component = c,
+          Group = g,
+          Parameter = p,
+          Variability = sd(d[[coefname]]),
+          Uncertainty = mean(d[[dispname]])
+        )
+        rez$Reliability <- rez$Variability / rez$Uncertainty
+
+        reliability <- rbind(reliability, rez)
+      }
+    }
+  }
+
+  # Clean-up output
+  if(length(unique(reliability$Component)) == 1 && unique(reliability$Component) == "TEMP") {
+    reliability$Component <- NULL
+  }
+
+  reliability
+}

man/check_reliability.Rd

@@ -0,0 +1,48 @@
+test_that("check_reliability - frequentist", {
+
+  skip_if_not_installed("lme4")
+  skip_if_not_installed("glmmTMB")
+
+  data <- iris
+  data$Place <- as.factor(rep(c("P1", "P2", "P3", "P4", "P5", "P6"), each = 25))
+
+  m <- lme4::lmer(Sepal.Width ~ Petal.Width + (Petal.Width | Place), data = data)
+  out <- check_reliability(m)
+  testthat::expect_true(all(dim(out) == c(2, 5)))
+
+  m <- glmmTMB::glmmTMB(Sepal.Width ~ Petal.Width + (Petal.Width | Place), data = data)
+  out <- check_reliability(m)
+  testthat::expect_true(all(dim(out) == c(2, 5)))
+})
+
+test_that("check_reliability - Bayesian", {
+
+  skip_if_offline()
+  skip_if_not_installed("curl")
+  skip_if_not_installed("httr2")
+  skip_if_not_installed("rstanarm")
+  skip_if_not_installed("rstantools")
+  skip_if_not_installed("brms")
+
+  library(rstanarm)
+
+
+  # m <- rstanarm::stan_lmer(
+  #   Sepal.Width ~ Petal.Width + (Petal.Width | Group),
+  #   data = d,
+  #   refresh = 0
+  # )
+  # m <- brms::brm(
+  #   Sepal.Width ~ Petal.Width + (Petal.Width | Group),
+  #   data = d,
+  #   refresh = 0,
+  #   backend="cmdstanr"
+  # )
+  m <- insight::download_model("brms_mixed_10")
+  out <- check_reliability(m)
+  testthat::expect_true(all(dim(out) == c(2, 5)))
+
+  m <- insight::download_model("brms_sigma_3")
+  out <- check_reliability(m)
+  testthat::expect_true(all(dim(out) == c(4, 6)))
+})