diff --git a/DESCRIPTION b/DESCRIPTION
index c02a0cbdd..892a937ad 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Type: Package
 Package: insight
 Title: Easy Access to Model Information for Various Model Objects
-Version: 0.20.1.1
+Version: 0.20.1.2
 Authors@R: 
     c(person(given = "Daniel",
              family = "Lüdecke",
@@ -210,3 +210,4 @@ Roxygen: list(markdown = TRUE)
 Config/testthat/edition: 3
 Config/testthat/parallel: true
 Config/Needs/website: easystats/easystatstemplate
+Remotes: easystats/performance
diff --git a/NEWS.md b/NEWS.md
index 0a967309b..36f944bb7 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,6 +2,10 @@
 
 ## General
 
+* Massive overhaul of `get_variance()`. The function should be now more
+  accurate for different distributional families, in particular for
+  mixed regression models with Beta family.
+
 * Improved accuracy of singularity-checks in `get_variance()`.
 
 # insight 0.20.1
diff --git a/R/compute_variances.R b/R/compute_variances.R
index ac4a448a7..560e430c6 100644
--- a/R/compute_variances.R
+++ b/R/compute_variances.R
@@ -3,8 +3,10 @@
                                name_fun = NULL,
                                name_full = NULL,
                                verbose = TRUE,
-                               tolerance = 1e-5,
-                               model_component = "conditional") {
+                               tolerance = 1e-8,
+                               model_component = "conditional",
+                               model_null = NULL,
+                               approximation = "lognormal") {
   ## Original code taken from GitGub-Repo of package glmmTMB
   ## Author: Ben Bolker, who used an cleaned-up/adapted
   ## version of Jon Lefcheck's code from SEMfit
@@ -17,6 +19,9 @@
 
   faminfo <- model_info(x, verbose = FALSE)
 
+  # check argument
+  approx_method <- match.arg(approximation, c("lognormal", "delta", "trigamma", "observation_level"))
+
   if (any(faminfo$family == "truncated_nbinom1")) {
     if (verbose) {
       format_warning(sprintf(
@@ -27,6 +32,11 @@
     return(NA)
   }
 
+  # rename lme4 neg-binom family
+  if (startsWith(faminfo$family, "Negative Binomial")) {
+    faminfo$family <- "negative binomial"
+  }
+
   # get necessary model information, like fixed and random effects,
   # variance-covariance matrix etc.
   vals <- .get_variance_information(
@@ -37,13 +47,28 @@
     model_component = model_component
   )
 
+  # we also need necessary model information, like fixed and random effects,
+  # variance-covariance matrix etc. for the null model
+  if (is.null(model_null)) {
+    model_null <- .safe(null_model(x, verbose = FALSE))
+  }
+  vals_null <- .get_variance_information(
+    model_null,
+    faminfo = faminfo,
+    name_fun = name_fun,
+    verbose = verbose,
+    model_component = model_component
+  )
+
   # Test for non-zero random effects ((near) singularity)
   no_random_variance <- FALSE
-  if (performance::check_singularity(x, tolerance = tolerance) && !(component %in% c("slope", "intercept"))) {
+  singular_fit <- isTRUE(.safe(performance::check_singularity(x, tolerance = tolerance)))
+
+  if (singular_fit && !(component %in% c("slope", "intercept"))) {
     if (verbose) {
       format_warning(
-        sprintf("Can't compute %s. Some variance components equal zero. Your model may suffer from singularity (see `?lme4::isSingular` and `?performance::check_singularity`).", name_full),
-        "Solution: Respecify random structure! You may also decrease the `tolerance` level to enforce the calculation of random effect variances."
+        sprintf("Can't compute %s. Some variance components equal zero. Your model may suffer from singularity (see `?lme4::isSingular` and `?performance::check_singularity`).", name_full), # nolint
+        "Solution: Respecify random structure! You may also decrease the `tolerance` level to enforce the calculation of random effect variances." # nolint
       )
     }
     no_random_variance <- TRUE
@@ -52,6 +77,7 @@
   # initialize return values, if not all components are requested
   var.fixed <- NULL
   var.random <- NULL
+  var.random_null <- NULL
   var.residual <- NULL
   var.distribution <- NULL
   var.dispersion <- NULL
@@ -83,6 +109,18 @@
     var.random <- .compute_variance_random(not.obs.terms, x = x, vals = vals)
   }
 
+  # Variance of random effects for NULL model
+  if (!singular_fit && !is.null(vals_null)) {
+    # Separate observation variance from variance of random effects
+    nr <- vapply(vals_null$re, nrow, numeric(1))
+    not.obs.terms_null <- names(nr[nr != n_obs(model_null)])
+    var.random_null <- .compute_variance_random(
+      not.obs.terms_null,
+      x = model_null,
+      vals = vals_null
+    )
+  }
+
   # Residual variance, which is defined as the variance due to
   # additive dispersion and the distribution-specific variance (Johnson et al. 2014)
 
@@ -91,6 +129,9 @@
       x = x,
       var.cor = vals$vc,
       faminfo,
+      model_null = model_null,
+      revar_null = var.random_null,
+      approx_method = approximation,
       name = name_full,
       verbose = verbose
     )
@@ -175,6 +216,11 @@
                                       name_fun = "get_variances",
                                       verbose = TRUE,
                                       model_component = "conditional") {
+  # sanity check
+  if (is.null(x)) {
+    return(NULL)
+  }
+
   # installed?
   check_if_installed("lme4", reason = "to compute variances for mixed models")
 
@@ -258,7 +304,7 @@
     )
     names(vals$re) <- re_names[seq_along(vals$re)]
 
-    # nlme
+    # nlme / glmmPQL
     # ---------------------------
   } else if (inherits(x, "lme")) {
     re_names <- find_random(x, split_nested = TRUE, flatten = TRUE)
@@ -419,6 +465,7 @@
 }
 
 
+# same as above, but for the zero-inflation component
 .collapse_zi <- function(x) {
   if (is.list(x) && "zi" %in% names(x)) {
     x[["zi"]]
@@ -431,8 +478,14 @@
 
 
 
-# fixed effects variance ----
-# ---------------------------
+# fixed effects variance ------------------------------------------------------
+#
+# This is in line with Nakagawa et al. 2017, Suppl. 2
+# see https://royalsocietypublishing.org/action/downloadSupplement?doi=10.1098%2Frsif.2017.0213&file=rsif20170213supp2.pdf
+#
+# However, package MuMIn differs and uses "fitted()" instead, leading to minor
+# deviations
+# -----------------------------------------------------------------------------
 .compute_variance_fixed <- function(vals) {
   with(vals, stats::var(as.vector(beta %*% t(X))))
 }
@@ -441,8 +494,11 @@
 
 
 
-# variance associated with a random-effects term (Johnson 2014) ----
-# ------------------------------------------------------------------
+# variance associated with a random-effects term (Johnson 2014) --------------
+#
+# This is in line with Nakagawa et al. 2017, Suppl. 2, and package MuMIm
+# see https://royalsocietypublishing.org/action/downloadSupplement?doi=10.1098%2Frsif.2017.0213&file=rsif20170213supp2.pdf
+# ----------------------------------------------------------------------------
 .compute_variance_random <- function(terms, x, vals) {
   if (is.null(terms)) {
     return(NULL)
@@ -480,17 +536,29 @@
 
 
 
-# distribution-specific variance (Nakagawa et al. 2017) ----
-# ----------------------------------------------------------
-.compute_variance_distribution <- function(x, var.cor, faminfo, name, verbose = TRUE) {
-  if (inherits(x, "lme")) {
-    sig <- x$sigma
-  } else {
-    sig <- attr(var.cor, "sc")
+# distribution-specific (residual) variance (Nakagawa et al. 2017) ------------
+#
+# This is in line with Nakagawa et al. 2017, Suppl. 2, and package MuMIm
+# see https://royalsocietypublishing.org/action/downloadSupplement?doi=10.1098%2Frsif.2017.0213&file=rsif20170213supp2.pdf
+#
+# There may be small deviations to Nakagawa et al. for the null-model, which
+# despite being correctly re-formulated in "null_model()", returns slightly
+# different values for the log/delta/trigamma approximation.
+# -----------------------------------------------------------------------------
+.compute_variance_distribution <- function(x,
+                                           var.cor,
+                                           faminfo,
+                                           model_null = NULL,
+                                           revar_null = NULL,
+                                           name,
+                                           approx_method = "lognormal",
+                                           verbose = TRUE) {
+  sig <- .safe(get_sigma(x))
+
+  if (is.null(sig)) {
+    sig <- 1
   }
 
-  if (is.null(sig)) sig <- 1
-
   # Distribution-specific variance depends on the model-family
   # and the related link-function
 
@@ -498,82 +566,175 @@
     # linear / Gaussian ----
     # ----------------------
 
-    dist.variance <- sig^2
+    resid.variance <- sig^2
   } else if (faminfo$is_betabinomial) {
     # beta-binomial ----
     # ------------------
 
-    dist.variance <- switch(faminfo$link_function,
+    resid.variance <- switch(faminfo$link_function,
       logit = ,
       probit = ,
       cloglog = ,
-      clogloglink = .variance_distributional(x, faminfo, sig, name = name, verbose = verbose),
+      clogloglink = .variance_distributional(
+        x,
+        faminfo = faminfo,
+        sig = sig,
+        approx_method = approx_method,
+        name = name,
+        verbose = verbose
+      ),
       .badlink(faminfo$link_function, faminfo$family, verbose = verbose)
     )
   } else if (faminfo$is_binomial) {
     # binomial / bernoulli  ----
     # --------------------------
 
-    dist.variance <- switch(faminfo$link_function,
-      logit = pi^2 / 3,
-      probit = 1,
+    # we need this to adjust for "cbind()" outcomes
+    y_factor <- .binomial_response_weight(x)
+
+    # for observation level approximation, when we don't want the "fixed"
+    # residual variance, pi^2/3, but the variance based on the distribution
+    # of the response
+    pmean <- .obs_level_variance(model_null, revar_null)
+
+    # sanity check - clmm-models are "binomial" but have no pmean
+    if (is.null(pmean) && identical(approx_method, "observation_level")) {
+      approx_method <- "lognormal"
+    }
+
+    resid.variance <- switch(faminfo$link_function,
+      logit = switch(approx_method,
+        observation_level = 1 / (y_factor * pmean * (1 - pmean)),
+        pi^2 / (3 * y_factor)
+      ),
+      probit = switch(approx_method,
+        observation_level = 2 * pi / y_factor * pmean * (1 - pmean) * exp((stats::qnorm(pmean) / sqrt(2))^2)^2, # nolint
+        1 / y_factor
+      ),
       cloglog = ,
-      clogloglink = pi^2 / 6,
+      clogloglink = switch(approx_method,
+        observation_level = pmean / y_factor / log(1 - pmean)^2 / (1 - pmean),
+        pi^2 / (6 * y_factor)
+      ),
       .badlink(faminfo$link_function, faminfo$family, verbose = verbose)
     )
+  } else if (faminfo$is_tweedie) {
+    # Tweedie  ----
+    # -------------
+
+    resid.variance <- .variance_distributional(
+      x,
+      faminfo = faminfo,
+      sig = sig,
+      model_null = model_null,
+      revar_null = revar_null,
+      approx_method = approx_method,
+      name = name,
+      verbose = verbose
+    )
   } else if (faminfo$is_count) {
     # count  ----
     # -----------
 
-    dist.variance <- switch(faminfo$link_function,
-      log = .variance_distributional(x, faminfo, sig, name = name, verbose = verbose),
-      sqrt = 0.25,
+    resid.variance <- switch(faminfo$link_function,
+      log = .variance_distributional(
+        x,
+        faminfo = faminfo,
+        sig = sig,
+        model_null = model_null,
+        revar_null = revar_null,
+        approx_method = approx_method,
+        name = name,
+        verbose = verbose
+      ),
+      sqrt = 0.25 * sig,
       .badlink(faminfo$link_function, faminfo$family, verbose = verbose)
     )
   } else if (faminfo$family %in% c("Gamma", "gamma")) {
     # Gamma  ----
     # -----------
 
-    ## TODO needs some more checking - should now be in line with other packages
-    dist.variance <- switch(faminfo$link_function,
+    resid.variance <- switch(faminfo$link_function,
       inverse = ,
-      identity = ,
-      log = stats::family(x)$variance(sig),
-      # log = .variance_distributional(x, faminfo, sig, name = name, verbose = verbose),
+      identity = stats::family(x)$variance(sig),
+      log = switch(approx_method,
+        delta = 1 / sig^-2,
+        trigamma = trigamma(sig^-2),
+        log1p(1 / sig^-2)
+      ),
       .badlink(faminfo$link_function, faminfo$family, verbose = verbose)
     )
   } else if (faminfo$family == "beta") {
     # Beta  ----
     # ----------
 
-    dist.variance <- switch(faminfo$link_function,
-      logit = .variance_distributional(x, faminfo, sig, name = name, verbose = verbose),
-      .badlink(faminfo$link_function, faminfo$family, verbose = verbose)
-    )
-  } else if (faminfo$is_tweedie) {
-    # Tweedie  ----
-    # -------------
-
-    dist.variance <- switch(faminfo$link_function,
-      log = .variance_distributional(x, faminfo, sig, name = name, verbose = verbose),
+    resid.variance <- switch(faminfo$link_function,
+      logit = .variance_distributional(
+        x,
+        faminfo = faminfo,
+        sig = sig,
+        model_null = model_null,
+        revar_null = revar_null,
+        name = name,
+        approx_method = approx_method,
+        verbose = verbose
+      ),
       .badlink(faminfo$link_function, faminfo$family, verbose = verbose)
     )
   } else if (faminfo$is_orderedbeta) {
     # Ordered Beta  ----
     # ------------------
 
-    dist.variance <- switch(faminfo$link_function,
-      logit = .variance_distributional(x, faminfo, sig, name = name, verbose = verbose),
+    resid.variance <- switch(faminfo$link_function,
+      logit = .variance_distributional(
+        x,
+        faminfo = faminfo,
+        sig = sig,
+        model_null = model_null,
+        revar_null = revar_null,
+        approx_method = approx_method,
+        name = name,
+        verbose = verbose
+      ),
       .badlink(faminfo$link_function, faminfo$family, verbose = verbose)
     )
   } else {
-    dist.variance <- sig
+    resid.variance <- sig
   }
 
-  dist.variance
+  resid.variance
 }
 
 
+# helper for .compute_variance_distribution()
+#
+# calculate weight if we have a "cbind()" response for binomial models
+# needed to weight the residual variance
+.binomial_response_weight <- function(x) {
+  # we need this to adjust for "cbind()" outcomes
+  resp_value <- .safe(stats::model.frame(x)[, find_response(x)])
+  # sanity check
+  if (is.null(resp_value)) {
+    resp_value <- .safe(get_data(x, source = "frame")[, find_response(x)])
+  }
+  if (!is.null(resp_value) && !is.null(ncol(resp_value)) && ncol(resp_value) > 1) {
+    mean(rowSums(resp_value, na.rm = TRUE))
+  } else {
+    1
+  }
+}
+
+
+# helper for .compute_variance_distribution()
+#
+# for observation level approximation, when we don't want the "fixed"
+# residual variance, pi^2/3, but the variance based on the distribution
+# of the response
+.obs_level_variance <- function(model_null, revar_null) {
+  fe_null <- .safe(as.numeric(.collapse_cond(lme4::fixef(model_null))))
+  .safe(as.numeric(stats::plogis(fe_null - 0.5 * sum(revar_null) * tanh(fe_null * (1 + 2 * exp(-0.5 * sum(revar_null))) / 6)))) # nolint
+}
+
 
 
 
@@ -594,133 +755,194 @@
 
 
 # This is the core-function to calculate the distribution-specific variance
-# Nakagawa et al. 2017 propose three different methods, here we only rely
-# on the lognormal-approximation.
+# Nakagawa et al. 2017 propose three different methods, which are now also
+# implemented here.
+#
+# This is in line with Nakagawa et al. 2017, Suppl. 2, and package MuMIm
+# see https://royalsocietypublishing.org/action/downloadSupplement?doi=10.1098%2Frsif.2017.0213&file=rsif20170213supp2.pdf
+#
+# There may be small deviations to Nakagawa et al. for the null-model, which
+# despite being correctly re-formulated in "null_model()", returns slightly
+# different values for the log/delta/trigamma approximation.
 #
-.variance_distributional <- function(x, faminfo, sig, name, verbose = TRUE) {
+# what we get here is the following rom the Nakagawa et al. Supplement:
+# VarOdF <- 1 / lambda + 1 / thetaF # the delta method
+# VarOlF <- log(1 + (1 / lambda) + (1 / thetaF)) # log-normal approximation
+# VarOtF <- trigamma((1 / lambda + 1 / thetaF)^(-1)) # trigamma function
+# -----------------------------------------------------------------------------
+.variance_distributional <- function(x,
+                                     faminfo,
+                                     sig,
+                                     model_null = NULL,
+                                     revar_null = NULL,
+                                     name,
+                                     approx_method = "lognormal",
+                                     verbose = TRUE) {
   check_if_installed("lme4", "to compute variances for mixed models")
 
-  # lognormal-approximation of distributional variance,
-  # see Nakagawa et al. 2017
-
-  # in general want log(1+var(x)/mu^2)
-  .null_model <- null_model(x, verbose = verbose)
+  # ------------------------------------------------------------------
+  # approximation of distributional variance, see Nakagawa et al. 2017
+  # in general want somethinh lije log(1+var(x)/mu^2) for log-approximation,
+  # and for the other approximations accordingly
+  # ------------------------------------------------------------------
 
   # check if null-model could be computed
-  if (is.null(.null_model)) {
+  if (is.null(model_null)) {
     mu <- NA
   } else {
-    if (inherits(.null_model, "cpglmm")) {
+    if (inherits(model_null, "cpglmm")) {
       # installed?
       check_if_installed("cplm")
-      null_fixef <- unname(cplm::fixef(.null_model))
+      null_fixef <- unname(cplm::fixef(model_null))
     } else {
-      null_fixef <- unname(.collapse_cond(lme4::fixef(.null_model)))
+      null_fixef <- unname(.collapse_cond(lme4::fixef(model_null)))
     }
     # brmsfit also returns SE and CI, so we just need the first value
-    if (inherits(.null_model, "brmsfit")) {
+    if (inherits(model_null, "brmsfit")) {
       null_fixef <- as.vector(null_fixef)[1]
     }
-
     mu <- null_fixef
   }
 
   if (is.na(mu)) {
     if (verbose) {
       format_warning(
-        "Can't calculate model's distribution-specific variance. Results are not reliable."
+        "Can't calculate model's distribution-specific variance. Results are not reliable.",
+        "A reason can be that the null model could not be computed manually. Try to fit the null model manually and pass it to `null_model`." # nolint
       )
     }
     return(0)
-  } else if (is.null(faminfo$family)) {
-    mu <- exp(mu)
+  }
+
+  # transform expected mean of the null model
+  if (is.null(faminfo$family)) {
+    mu <- link_inverse(x)(mu)
   } else {
     # transform mu
     mu <- switch(faminfo$family,
-      beta = mu,
+      beta = ,
+      betabinomial = ,
       ordbeta = stats::plogis(mu),
-      exp(mu)
-    )
-  }
-
-  # check if mu is too close to zero, but not for beta-distribution
-  if (mu < 6 && verbose && isFALSE(faminfo$family %in% c("beta", "ordbeta"))) {
-    format_warning(
-      sprintf("mu of %0.1f is too close to zero, estimate of %s may be unreliable.", mu, name)
+      poisson = ,
+      quasipoisson = ,
+      nbinom = ,
+      nbinom1 = ,
+      nbinom2 = ,
+      tweedie = ,
+      `negative binomial` = exp(mu + 0.5 * as.vector(revar_null)),
+      link_inverse(x)(mu) ## TODO: check if this is better than "exp(mu)"
+      # exp(mu)
     )
   }
 
   cvsquared <- tryCatch(
     {
-      vv <- switch(faminfo$family,
-
-        # (zero-inflated) poisson ----
-        # ----------------------------
-        `zero-inflated poisson` = ,
-        poisson = .variance_family_poisson(x, mu, faminfo),
-
-        # hurdle-poisson ----
-        # -------------------
-        `hurdle poisson` = ,
-        truncated_poisson = stats::family(x)$variance(sig),
-
-        # Gamma, exponential ----
-        # -----------------------
-        Gamma = stats::family(x)$variance(sig),
-
-        # (zero-inflated) negative binomial ----
-        # --------------------------------------
-        `zero-inflated negative binomial` = ,
-        `negative binomial` = ,
-        genpois = ,
-        nbinom1 = ,
-        nbinom2 = .variance_family_nbinom(x, mu, sig, faminfo),
-        truncated_nbinom2 = stats::family(x)$variance(mu, sig),
-
-        # other distributions ----
-        # ------------------------
-        tweedie = .variance_family_tweedie(x, mu, sig),
-        beta = .variance_family_beta(x, mu, sig),
-        ordbeta = .variance_family_orderedbeta(x, mu),
-        # betabinomial = stats::family(x)$variance(mu, sig),
-        # betabinomial = .variance_family_betabinom(x, mu, sig),
-
-        # default variance for non-captured distributions ----
-        # ----------------------------------------------------
-        .variance_family_default(x, mu, verbose)
-      )
+      if (faminfo$link_function == "tweedie") {
+        vv <- .variance_family_tweedie(x, mu, sig)
+      } else {
+        vv <- switch(faminfo$family,
+
+          # (zero-inflated) poisson ----
+          # ----------------------------
+          `zero-inflated poisson` = .variance_family_poisson(x, mu, faminfo),
+          poisson = 1,
+
+          # hurdle-poisson ----
+          # -------------------
+          `hurdle poisson` = ,
+          truncated_poisson = stats::family(x)$variance(sig),
+
+          # Gamma, exponential ----
+          # -----------------------
+          Gamma = stats::family(x)$variance(sig),
+
+          # (zero-inflated) negative binomial ----
+          # --------------------------------------
+          nbinom = ,
+          nbinom1 = ,
+          nbinom2 = ,
+          quasipoisson = ,
+          `negative binomial` = sig,
+          `zero-inflated negative binomial` = ,
+          genpois = .variance_family_nbinom(x, mu, sig, faminfo),
+          truncated_nbinom2 = stats::family(x)$variance(mu, sig),
+
+          # other distributions ----
+          # ------------------------
+          tweedie = .variance_family_tweedie(x, mu, sig),
+          beta = ,
+          betabinomial = .variance_family_beta(x, mu, sig),
+          ordbeta = .variance_family_orderedbeta(x, mu),
+          # betabinomial = stats::family(x)$variance(mu, sig),
+          # betabinomial = .variance_family_betabinom(x, mu, sig),
+
+          # default variance for non-captured distributions ----
+          # ----------------------------------------------------
+          .variance_family_default(x, mu, verbose)
+        )
+      }
 
       if (vv < 0 && isTRUE(verbose)) {
         format_warning(
           "Model's distribution-specific variance is negative. Results are not reliable."
         )
       }
-      vv / mu^2
+
+      # now compute cvsquared -------------------------------------------
+
+      # for cpglmm with tweedie link, the model is not of tweedie family,
+      # only the link function is tweedie
+      if (faminfo$link_function == "tweedie") {
+        vv
+      } else if (identical(approx_method, "trigamma")) {
+        switch(faminfo$family,
+          nbinom = ,
+          nbinom1 = ,
+          nbinom2 = ,
+          `negative binomial` = ((1 / mu) + (1 / sig))^-1,
+          poisson = ,
+          quasipoisson = mu / vv,
+          vv / mu^2
+        )
+      } else {
+        switch(faminfo$family,
+          nbinom = ,
+          nbinom1 = ,
+          nbinom2 = ,
+          `negative binomial` = (1 / mu) + (1 / sig),
+          poisson = ,
+          quasipoisson = vv / mu,
+          vv / mu^2
+        )
+      }
     },
     error = function(x) {
       if (verbose) {
         format_warning(
-          "Can't calculate model's distribution-specific variance. Results are not reliable."
+          "Can't calculate model's distribution-specific variance. Results are not reliable.",
+          paste("The following error occured: ", x$message)
         )
       }
       0
     }
   )
 
-  log1p(cvsquared)
+  switch(approx_method,
+    delta = cvsquared,
+    trigamma = trigamma(cvsquared),
+    log1p(cvsquared)
+  )
 }
 
 
 
-
-
 # Get distributional variance for poisson-family
 # ----------------------------------------------
 .variance_family_poisson <- function(x, mu, faminfo) {
   if (faminfo$is_zero_inflated) {
     .variance_zip(x, faminfo, family_var = mu)
   } else if (inherits(x, "MixMod")) {
-    return(mu)
+    mu
   } else if (inherits(x, "cpglmm")) {
     .get_cplm_family(x)$variance(mu)
   } else {
@@ -730,22 +952,22 @@
 
 
 
-
-
 # Get distributional variance for beta-family
 # ----------------------------------------------
 .variance_family_beta <- function(x, mu, phi) {
-  if (inherits(x, "MixMod")) {
-    stats::family(x)$variance(mu)
-  } else {
-    mu * (1 - mu) / (1 + phi)
-  }
+  stats::family(x)$variance(mu)
+  # if (inherits(x, "MixMod")) {
+  #   stats::family(x)$variance(mu)
+  # } else {
+  #   # was:
+  #   # mu * (1 - mu) / (1 + phi)
+  #   # but that code is not what "glmmTMB" uses for the beta family
+  #   mu * (1 - mu)
+  # }
 }
 
 
 
-
-
 # Get distributional variance for ordered beta-family
 # ----------------------------------------------
 .variance_family_orderedbeta <- function(x, mu) {
@@ -758,8 +980,6 @@
 
 
 
-
-
 # Get distributional variance for beta-family
 # ----------------------------------------------
 .variance_family_betabinom <- function(x, mu, phi) {
@@ -773,24 +993,25 @@
 
 
 
-
-
 # Get distributional variance for tweedie-family
 # ----------------------------------------------
 .variance_family_tweedie <- function(x, mu, phi) {
-  if ("psi" %in% names(x$fit$par)) {
-    psi <- x$fit$par["psi"] # glmmmTMB >= 1.1.5
+  if (inherits(x, "cpglmm")) {
+    phi <- x@phi
+    p <- x@p - 2
   } else {
-    psi <- x$fit$par["thetaf"]
+    if ("psi" %in% names(x$fit$par)) {
+      psi <- x$fit$par["psi"] # glmmmTMB >= 1.1.5
+    } else {
+      psi <- x$fit$par["thetaf"]
+    }
+    p <- unname(stats::plogis(psi) + 1)
   }
-  p <- unname(stats::plogis(psi) + 1)
   phi * mu^p
 }
 
 
 
-
-
 # Get distributional variance for nbinom-family
 # ----------------------------------------------
 .variance_family_nbinom <- function(x, mu, sig, faminfo) {
@@ -809,8 +1030,6 @@
 
 
 
-
-
 # For zero-inflated negative-binomial models,
 # the distributional variance is based on Zuur et al. 2012
 # ----------------------------------------------
@@ -849,8 +1068,6 @@
 
 
 
-
-
 # For zero-inflated poisson models, the
 # distributional variance is based on Zuur et al. 2012
 # ----------------------------------------------
@@ -872,8 +1089,6 @@
 
 
 
-
-
 # Get distribution-specific variance for general and
 # undefined families / link-functions
 # ----------------------------------------------
@@ -1066,7 +1281,7 @@
   # check if any polynomial / I term in random slopes.
   # we then have correlation among levels
   rs_names <- unique(unlist(lapply(corrs, colnames)))
-  pattern <- paste0("(I|poly)(.*)(", paste0(rnd_slopes, collapse = "|"), ")")
+  pattern <- paste0("(I|poly)(.*)(", paste(rnd_slopes, collapse = "|"), ")")
   poly_random_slopes <- any(grepl(pattern, rs_names))
 
   if (length(rnd_slopes) < 2 && !isTRUE(cat_random_slopes) && !isTRUE(poly_random_slopes)) {
diff --git a/R/get_nested_lme_varcorr.R b/R/get_nested_lme_varcorr.R
index 78834cf6a..13215a995 100644
--- a/R/get_nested_lme_varcorr.R
+++ b/R/get_nested_lme_varcorr.R
@@ -8,6 +8,7 @@
   vcor <- lme4::VarCorr(x)
   class(vcor) <- "matrix"
 
+  ## FIXME: doesn't work for nested RE from MASS::glmmPQL, see Nakagawa example
   re_index <- (which(rownames(vcor) == "(Intercept)") - 1)[-1]
   vc_list <- split(data.frame(vcor, stringsAsFactors = FALSE), findInterval(seq_len(nrow(vcor)), re_index))
   vc_rownames <- split(rownames(vcor), findInterval(seq_len(nrow(vcor)), re_index))
@@ -49,5 +50,9 @@
 
 
 .is_nested_lme <- function(x) {
-  sapply(find_random(x), function(i) any(grepl(":", i, fixed = TRUE)))
+  if (inherits(x, "glmmPQL")) {
+    length(find_random(x, flatten = TRUE)) > 1
+  } else {
+    sapply(find_random(x), function(i) any(grepl(":", i, fixed = TRUE)))
+  }
 }
diff --git a/R/get_sigma.R b/R/get_sigma.R
index 2c6318236..b4fd407ed 100644
--- a/R/get_sigma.R
+++ b/R/get_sigma.R
@@ -74,6 +74,50 @@ get_sigma <- function(x, ci = NULL, verbose = TRUE) {
 # special handling ---------------
 
 
+.get_sigma.glmerMod <- function(x, ...) {
+  check_if_installed("lme4")
+  if (startsWith(stats::family(x)$family, "Negative Binomial(")) {
+    lme4::getME(x, "glmer.nb.theta")
+  } else {
+    stats::sigma(x)
+  }
+}
+
+
+.get_sigma.glmmadmb <- function(x, ...) {
+  check_if_installed("lme4")
+  vc <- lme4::VarCorr(x)
+  s <- attr(vc, "sc")
+  # sanity check
+  if (is.null(s)) {
+    s <- .safe(x$alpha)
+  }
+  s
+}
+
+
+.get_sigma.glmmTMB <- function(x, ...) {
+  # The commented code is what MuMIn returns for sigma for nbinom1 models.
+  # However, I think this is wrong. Nakagawa et al. (2017) used this in their
+  # code because glmmadmb models with nbinom1 family are actually Quasi-Poisson
+  # models (see also Supplement 2). Thus, we revert and just use "sigma()" again.
+  # This will return results for `get_variance()` that are in line with the code
+  # in the Supplement 2 from Nakaawa et al. (2017).
+  # if (stats::family(x)$family == "nbinom1") {
+  #   add_value <- 1
+  # } else {
+  #   add_value <- 0
+  # }
+  # stats::sigma(x) + add_value
+  stats::sigma(x)
+}
+
+
+.get_sigma.merMod <- function(x, ...) {
+  stats::sigma(x)
+}
+
+
 .get_sigma.model_fit <- function(x, verbose = TRUE, ...) {
   .get_sigma(x$fit, verbose = verbose)
 }
@@ -144,6 +188,25 @@ get_sigma <- function(x, ci = NULL, verbose = TRUE) {
 }
 
 
+.get_sigma.lme <- function(x, ...) {
+  .safe(x$sigma)
+}
+
+
+.get_sigma.mjoint <- function(x, ...) {
+  .safe(x$coef$sigma2[[1]])
+}
+
+
+.get_sigma.glmmPQL <- function(x, ...) {
+  switch(x$family$family,
+    gaussian = ,
+    Gamma = x$sigma,
+    x$sigma^2
+  )
+}
+
+
 .get_sigma.brmsfit <- function(x, verbose = TRUE, ...) {
   s <- tryCatch(
     {
@@ -206,15 +269,6 @@ get_sigma <- function(x, ci = NULL, verbose = TRUE) {
     }
   }
 
-  if (is_empty_object(s)) {
-    if (is.null(info)) {
-      info <- model_info(x, verbose = FALSE)
-    }
-    if (!is.null(info) && info$is_mixed) {
-      s <- .safe(sqrt(get_variance_residual(x, verbose = FALSE)))
-    }
-  }
-
   if (is_empty_object(s)) {
     s <- .safe(
       sqrt(get_deviance(x, verbose = verbose) / get_df(x, type = "residual", verbose = verbose))
@@ -254,7 +308,6 @@ get_sigma <- function(x, ci = NULL, verbose = TRUE) {
 as.numeric.insight_aux <- function(x, ...) {
   if (is.null(x) || is.na(x) || is.infinite(x)) {
     return(NULL)
-  } else {
-    mean(x, na.rm = TRUE)
   }
+  mean(x, na.rm = TRUE)
 }
diff --git a/R/get_variances.R b/R/get_variances.R
index c14da343e..3696c7dbd 100644
--- a/R/get_variances.R
+++ b/R/get_variances.R
@@ -10,15 +10,24 @@
 #'
 #' @param x A mixed effects model.
 #' @param component Character value, indicating the variance component that should
-#'   be returned. By default, all variance components are returned. The
-#'   distribution-specific (`"distribution"`) and residual (`"residual"`)
-#'   variance are the most computational intensive components, and hence may
-#'   take a few seconds to calculate.
+#' be returned. By default, all variance components are returned. The
+#' distribution-specific (`"distribution"`) and residual (`"residual"`)
+#' variance are the most computational intensive components, and hence may
+#' take a few seconds to calculate.
 #' @param verbose Toggle off warnings.
 #' @param tolerance Tolerance for singularity check of random effects, to decide
-#'   whether to compute random effect variances or not. Indicates up to which
-#'   value the convergence result is accepted. The larger tolerance is, the
-#'   stricter the test will be. See [performance::check_singularity()].
+#' whether to compute random effect variances or not. Indicates up to which
+#' value the convergence result is accepted. The larger tolerance is, the
+#' stricter the test will be. See [performance::check_singularity()].
+#' @param null_model Optional, a null-model to be used for the calculation of
+#' random effect variances. If `NULL`, the null-model is computed internally.
+#' @param approximation Character string, indicating the approximation method
+#' for the distribution-specific (residual) variance. Only applies to non-Gaussian
+#' models. Can be `"lognormal"` (default), `"delta"` or `"trigamma"`. For binomial
+#' models, can also be `"observation_level"`. See _Nakagawa et al. 2017_ for
+#' details.
+#' @param model_component For models that can have a zero-inflation component,
+#' specify for which component variances should be returned.
 #' @param ... Currently not used.
 #'
 #' @return A list with following elements:
@@ -124,18 +133,10 @@
 #' get_variance_residual(m)
 #' }
 #' @export
-get_variance <- function(x,
-                         component = c(
-                           "all", "fixed", "random", "residual",
-                           "distribution", "dispersion", "intercept",
-                           "slope", "rho01", "rho00"
-                         ),
-                         verbose = TRUE,
-                         ...) {
+get_variance <- function(x, ...) {
   UseMethod("get_variance")
 }
 
-
 #' @export
 get_variance.default <- function(x,
                                  component = c(
@@ -152,6 +153,7 @@ get_variance.default <- function(x,
 }
 
 
+#' @rdname get_variance
 #' @export
 get_variance.merMod <- function(x,
                                 component = c(
@@ -159,8 +161,10 @@ get_variance.merMod <- function(x,
                                   "distribution", "dispersion", "intercept",
                                   "slope", "rho01", "rho00"
                                 ),
+                                tolerance = 1e-8,
+                                null_model = NULL,
+                                approximation = "lognormal",
                                 verbose = TRUE,
-                                tolerance = 1e-5,
                                 ...) {
   component <- match.arg(component)
   .safe(.compute_variances(
@@ -169,7 +173,9 @@ get_variance.merMod <- function(x,
     name_fun = "get_variance",
     name_full = "random effect variances",
     verbose = verbose,
-    tolerance = tolerance
+    tolerance = tolerance,
+    model_null = null_model,
+    approximation = approximation
   ))
 }
 
@@ -205,6 +211,7 @@ get_variance.lme <- get_variance.merMod
 get_variance.brmsfit <- get_variance.merMod
 
 
+#' @rdname get_variance
 #' @export
 get_variance.glmmTMB <- function(x,
                                  component = c(
@@ -212,9 +219,11 @@ get_variance.glmmTMB <- function(x,
                                    "distribution", "dispersion", "intercept",
                                    "slope", "rho01", "rho00"
                                  ),
-                                 verbose = TRUE,
-                                 tolerance = 1e-5,
                                  model_component = NULL,
+                                 tolerance = 1e-8,
+                                 null_model = NULL,
+                                 approximation = "lognormal",
+                                 verbose = TRUE,
                                  ...) {
   component <- match.arg(component)
   .safe(.compute_variances(
@@ -224,7 +233,9 @@ get_variance.glmmTMB <- function(x,
     name_full = "random effect variances",
     verbose = verbose,
     tolerance = tolerance,
-    model_component = model_component
+    model_component = model_component,
+    model_null = null_model,
+    approximation = approximation
   ))
 }
 
@@ -239,8 +250,10 @@ get_variance.mixed <- function(x,
                                  "distribution", "dispersion", "intercept",
                                  "slope", "rho01", "rho00"
                                ),
+                               tolerance = 1e-8,
+                               null_model = NULL,
+                               approximation = "lognormal",
                                verbose = TRUE,
-                               tolerance = 1e-5,
                                ...) {
   component <- match.arg(component)
   .safe(.compute_variances(
@@ -249,7 +262,9 @@ get_variance.mixed <- function(x,
     name_fun = "get_variance",
     name_full = "random effect variances",
     verbose = verbose,
-    tolerance = tolerance
+    tolerance = tolerance,
+    model_null = null_model,
+    approximation = approximation
   ))
 }
 
@@ -268,7 +283,7 @@ get_variance_fixed <- function(x, verbose = TRUE, ...) {
 
 #' @rdname get_variance
 #' @export
-get_variance_random <- function(x, verbose = TRUE, tolerance = 1e-5, ...) {
+get_variance_random <- function(x, verbose = TRUE, tolerance = 1e-8, ...) {
   unlist(get_variance(x, component = "random", verbose = verbose, tolerance = tolerance, ...))
 }
 
diff --git a/R/null_model.R b/R/null_model.R
index f8e81e2e2..09c49111b 100644
--- a/R/null_model.R
+++ b/R/null_model.R
@@ -39,25 +39,23 @@ null_model <- function(model, verbose = TRUE, ...) {
     stats::update(model, location = ~1, scale = ~1)
   } else if (inherits(model, "multinom")) {
     stats::update(model, ~1, trace = FALSE)
+  } else if (is.null(offset_term)) {
+    # stats::update(model, ~1)
+    out <- stats::update(model, ~1, evaluate = FALSE)
+    eval(out, envir = parent.frame())
   } else {
-    if (!is.null(offset_term)) {
-      tryCatch(
-        do.call(stats::update, list(model, ~1, offset = str2lang(offset_term))),
-        error = function(e) {
-          if (verbose) {
-            format_warning(
-              "Model contains offset-terms, which could not be considered in the returned null-model.",
-              "Coefficients might be inaccurate."
-            )
-          }
-          stats::update(model, ~1)
+    tryCatch(
+      do.call(stats::update, list(model, ~1, offset = str2lang(offset_term))),
+      error = function(e) {
+        if (verbose) {
+          format_warning(
+            "Model contains offset-terms, which could not be considered in the returned null-model.",
+            "Coefficients might be inaccurate."
+          )
         }
-      )
-    } else {
-      # stats::update(model, ~1)
-      out <- stats::update(model, ~1, evaluate = FALSE)
-      eval(out, envir = parent.frame())
-    }
+        stats::update(model, ~1)
+      }
+    )
   }
 }
 
@@ -65,10 +63,10 @@ null_model <- function(model, verbose = TRUE, ...) {
 .null_model_mixed <- function(model, offset_term = NULL, verbose = TRUE) {
   if (inherits(model, "MixMod")) {
     nullform <- stats::as.formula(paste(find_response(model), "~ 1"))
-    null.model <- stats::update(model, fixed = nullform)
+    null.model <- suppressWarnings(stats::update(model, fixed = nullform))
   } else if (inherits(model, "cpglmm")) {
     nullform <- find_formula(model, verbose = FALSE)[["random"]]
-    null.model <- stats::update(model, nullform)
+    null.model <- suppressWarnings(stats::update(model, nullform))
   } else {
     f <- stats::formula(model)
     resp <- find_response(model)
@@ -79,21 +77,21 @@ null_model <- function(model, verbose = TRUE, ...) {
     re.terms <- paste0("(", sapply(.findbars(f), safe_deparse), ")")
     nullform <- stats::reformulate(re.terms, response = resp)
     null.model <- tryCatch(
-      if (!is.null(offset_term)) {
-        do.call(stats::update, list(model, formula = nullform, offset = str2lang(offset_term)))
+      if (is.null(offset_term)) {
+        suppressWarnings(stats::update(model, nullform))
       } else {
-        stats::update(model, nullform)
+        suppressWarnings(do.call(stats::update, list(model, formula = nullform, offset = str2lang(offset_term))))
       },
       error = function(e) {
         msg <- e$message
         if (verbose) {
           if (grepl("(^object)(.*)(not found$)", msg)) {
-            print_color("Can't calculate null-model. Probably the data that was used to fit the model cannot be found.\n", "red")
+            print_color("Can't calculate null-model. Probably the data that was used to fit the model cannot be found.\n", "red") # nolint
           } else if (startsWith(msg, "could not find function")) {
-            print_color("Can't calculate null-model. Probably you need to load the package that was used to fit the model.\n", "red")
+            print_color("Can't calculate null-model. Probably you need to load the package that was used to fit the model.\n", "red") # nolint
           }
         }
-        return(NULL)
+        NULL
       }
     )
   }
diff --git a/WIP/nakagawa_suppl.R b/WIP/nakagawa_suppl.R
new file mode 100644
index 000000000..d61385da4
--- /dev/null
+++ b/WIP/nakagawa_suppl.R
@@ -0,0 +1,509 @@
+vsCodeSnippets::load_debug_pkg()
+
+# installing glmmADMB
+library(R2admb)
+library(glmmADMB)
+library(lme4)
+library(cplm)
+library(performance)
+library(insight)
+library(glmmTMB)
+
+
+set.seed(1234)
+# 12 different populations n = 960
+Population <- gl(12, 80, 960)
+# 120 containers (8 individuals in each container)
+Container <- gl(120, 8, 960)
+# Sex of the individuals. Uni-sex within each container (individuals are
+# sorted at the pupa stage)
+Sex <- factor(rep(rep(c("Female", "Male"), each = 8), 60))
+# Habitat at the collection site: dry or wet soil (four indiviudal from each
+# Habitat in each container)
+Habitat <- factor(rep(rep(c("Dry", "Wet"), each = 4), 120))
+# Food treatment at the larval stage: special food ('Exp') or standard food
+# ('Cont')
+Treatment <- factor(rep(c("Cont", "Exp"), 480))
+# Data combined in a data frame
+Data <- data.frame(
+  Population = Population, Container = Container, Sex = Sex,
+  Habitat = Habitat, Treatment = Treatment
+)
+
+# Subset the design matrix (only females lay eggs)
+DataFemale <- Data[Data$Sex == "Female", ]
+# set seed for reproduciblity (this will enable one to get the same data
+# every time)
+set.seed(777)
+# simulation of the underlying random effects (Population and Container with
+# variance of 0.4 and 0.05, respectively)
+PopulationE <- rnorm(12, 0, sqrt(0.4))
+ContainerE <- rnorm(120, 0, sqrt(0.05))
+# generation of response values on latent scale (!) based on fixed effects,
+# random effects and residual errors
+EggL <- with(DataFemale, 1.1 + 0.5 * (as.numeric(Treatment) - 1) + 0.1 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container] + rnorm(480, 0, sqrt(0.1)))
+# data generation (on data scale!) based on Poisson distribution
+DataFemale$Egg <- rpois(length(EggL), exp(EggL))
+
+# Data frame for both sex
+DataAll <- Data
+# simulation of the underlying random effects (Population and Container with
+# variance of 0.5 and 0.8, respectively)
+PopulationE <- rnorm(12, 0, sqrt(0.5))
+ContainerE <- rnorm(120, 0, sqrt(0.8))
+# generation of response values on latent scale (!) based on fixed effects
+# and random effects
+ParasiteL <- with(DataAll, 1.8 + 2 * (-1) * (as.numeric(Sex) - 1) + 0.8 * (-1) * (as.numeric(Treatment) - 1) + 0.7 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container])
+# data generation (on data scale!) based on negative binomial distributions;
+# size = theta
+DataAll$Parasite <- rnbinom(length(ParasiteL), size = 5, mu = exp(ParasiteL))
+
+# simulation of the underlying random effects (Population and Container with
+# variance of 1.3 and 0.3, respectively)
+PopulationE <- rnorm(12, 0, sqrt(1.3))
+ContainerE <- rnorm(120, 0, sqrt(0.3))
+# data generation based on fixed effects, random effects and random
+# residuals errors
+DataAll$BodyL <- 15 + 3 * (-1) * (as.numeric(Sex) - 1) + 0.4 * (as.numeric(Treatment) - 1) + 0.15 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container] + rnorm(960, 0, sqrt(1.2))
+# simulation of the underlying random effects (Population and Container with
+# variance of 0.2 and 0.2, respectively)
+PopulationE <- rnorm(12, 0, sqrt(0.2))
+ContainerE <- rnorm(120, 0, sqrt(0.2))
+# generation of response values on latent scale (!) based on fixed effects
+# and random effects
+ExplorationL <- with(DataAll, 4 + 1 * (-1) * (as.numeric(Sex) - 1) + 2 * (as.numeric(Treatment) - 1) + 0.5 * (-1) * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container])
+# data generation (on data scale!) based on gamma distribution; size = theta
+DataAll$Exploration <- rgamma(length(ExplorationL), shape = exp(ExplorationL) * 0.3, rate = 0.3)
+
+# Subset the design matrix (only males express colour morphs)
+DataMale <- subset(Data, Sex == "Male")
+# simulation of the underlying random effects (Population and Container with
+# variance of 1.2 and 0.2, respectively)
+PopulationE <- rnorm(12, 0, sqrt(1.2))
+ContainerE <- rnorm(120, 0, sqrt(0.2))
+# generation of response values on latent scale (!) based on fixed effects
+# and random effects
+ColourL <- with(DataMale, 0.8 * (-1) + 0.8 * (as.numeric(Treatment) - 1) + 0.5 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container])
+# data generation (on data scale!) based on binomial distribution
+DataMale$Colour <- rbinom(length(ColourL), 1, plogis(ColourL))
+
+
+
+# ==============================================================
+# Quasi-Poisson with log link (page 5) -------------------------
+# glmmadmb -----------------------------------------------------
+# ==============================================================
+
+# Fit null model without fixed effects (but including all random effects)
+fecmodADMBr <- glmmadmb(
+  Egg ~ 1 + (1 | Population) + (1 | Container),
+  family = "nbinom1", data = DataFemale
+)
+# Fit alternative model including fixed and all random effects
+fecmodADMBf <- glmmadmb(
+  Egg ~ Treatment + Habitat + (1 | Population) + (1 | Container),
+  family = "nbinom1", data = DataFemale
+)
+
+# Calculation of the variance in fitted values
+VarF <- var(as.vector(model.matrix(fecmodADMBf) %*% fixef(fecmodADMBf)))
+# getting the observation-level variance Null model
+omegaN <- fecmodADMBr$alpha # overdispersion omega is alpha in glmmadmb
+lambda <- as.numeric(exp(fixef(fecmodADMBr) + 0.5 * (as.numeric(VarCorr(fecmodADMBr)[1]) + as.numeric(VarCorr(fecmodADMBr)[2]))))
+# lambda2 <- mean(DataFemale$Egg) # for lambda we use the mean of all
+# observations
+VarOdN <- omegaN / lambda # the delta method
+VarOlN <- log(1 + omegaN / lambda) # log-normal approximation
+VarOtN <- trigamma(lambda / omegaN) # trigamma function
+# comparing the three
+c(VarOdN = VarOdN, VarOlN = VarOlN, VarOtN = VarOtN)
+
+# Full model
+omegaF <- fecmodADMBf$alpha # overdispersion omega is alpha in glmmadmb
+VarOdF <- omegaF / lambda # the delta method
+VarOlF <- log(1 + omegaF / lambda) # log-normal approximation
+VarOtF <- trigamma(lambda / omegaF) # trigamma function
+# comparing the three
+c(VarOdF = VarOdF, VarOlF = VarOlF, VarOtF = VarOtF)
+
+# R2[GLMM(m)] - marginal R2[GLMM]
+R2glmmM <- VarF / (VarF + sum(as.numeric(VarCorr(fecmodADMBf))) + VarOlF)
+# R2[GLMM(c)] - conditional R2[GLMM] for full model
+R2glmmC <- (VarF + sum(as.numeric(VarCorr(fecmodADMBf)))) / (VarF + sum(as.numeric(VarCorr(fecmodADMBf))) + VarOlF)
+# Raw unadjusted ICC[Population]
+ICCrawPop <- VarCorr(fecmodADMBr)$Population[1] / (sum(as.numeric(VarCorr(fecmodADMBr))) + VarOlN)
+# adjusted ICC[Population]
+ICCadjPop <- VarCorr(fecmodADMBf)$Population[1] / (sum(as.numeric(VarCorr(fecmodADMBf))) + VarOlF)
+# Raw unadjusted ICC[Container]
+ICCrawCont <- VarCorr(fecmodADMBr)$Container[1] / (sum(as.numeric(VarCorr(fecmodADMBr))) + VarOlN)
+# adjusted ICC[Container]
+ICCadjCont <- VarCorr(fecmodADMBf)$Container[1] / (sum(as.numeric(VarCorr(fecmodADMBf))) + VarOlF)
+# comparing the results
+c(R2glmmM = R2glmmM, R2glmmC = R2glmmC, ICCrawPop = ICCrawPop, ICCadjPop = ICCadjPop, ICCrawCont = ICCrawCont, ICCadjCont = ICCadjCont)
+
+performance::r2_nakagawa(fecmodADMBf)
+
+
+# ==============================================================
+# Quasi-Poisson with log link (page 7) -------------------------
+# glmmPQL- -----------------------------------------------------
+# ==============================================================
+
+# Fit null model without fixed effects (but including all random effects)
+fecmodPQLr <- glmmPQL(Egg ~ 1,
+  random = list(~ 1 | Population, ~ 1 | Container),
+  family = "quasipoisson", data = DataFemale
+)
+# Fit alternative model including fixed and all random effects
+fecmodPQLf <- glmmPQL(Egg ~ Treatment + Habitat, random = list(
+  ~ 1 | Population,
+  ~ 1 | Container
+), family = "quasipoisson", data = DataFemale)
+
+# Calculation of the variance in fitted values
+VarF <- var(as.vector(model.matrix(~ Treatment + Habitat, data = DataFemale) %*% fixef(fecmodPQLf)))
+# getting the observation-level variance Null model
+omegaN <- as.numeric(VarCorr(fecmodPQLr)[5, 1]) # overdispersion omega is residual variance in glmmPQL
+lambda <- as.numeric(exp(fixef(fecmodPQLr) + 0.5 * (as.numeric(VarCorr(fecmodPQLr)[2, 1]) + as.numeric(VarCorr(fecmodPQLr)[4, 1]))))
+# lambda2 <- mean(DataFemale$Egg)
+VarOdN <- omegaN / lambda
+
+VarOlN <- log(1 + omegaN / lambda)
+VarOtN <- trigamma(lambda / omegaN)
+# comparing the three
+c(VarOdN = VarOdN, VarOlN = VarOlN, VarOtN = VarOtN)
+
+omegaF <- as.numeric(VarCorr(fecmodPQLf)[5, 1]) # overdispersion omega is residual variance in glmmPQL
+VarOdF <- omegaF / lambda
+VarOlF <- log(1 + omegaF / lambda)
+VarOtF <- trigamma(lambda / omegaF)
+# comparing the three
+c(VarOdF = VarOdF, VarOlF = VarOlF, VarOtF = VarOtF)
+
+# R2[GLMM(m)] - marginal R2[GLMM]
+R2glmmM <- VarF / (VarF + sum(as.numeric(VarCorr(fecmodPQLf)[c(2, 4), 1])) + VarOlF)
+# R2[GLMM(c)] - conditional R2[GLMM] for full model
+R2glmmC <- (VarF + sum(as.numeric(VarCorr(fecmodPQLf)[c(2, 4), 1]))) / (VarF + sum(as.numeric(VarCorr(fecmodPQLf)[c(2, 4), 1])) + VarOlF)
+# Raw unadjusted ICC[Population]
+ICCrawPop <- as.numeric(VarCorr(fecmodPQLr)[2, 1]) / (sum(as.numeric(VarCorr(fecmodPQLr)[c(2, 4), 1])) + VarOlN)
+# adjusted ICC[Population]
+ICCadjPop <- as.numeric(VarCorr(fecmodPQLf)[2, 1]) / (sum(as.numeric(VarCorr(fecmodPQLf)[c(2, 4), 1])) + VarOlF)
+# Raw unadjusted ICC[Container]
+ICCrawCont <- as.numeric(VarCorr(fecmodPQLr)[4, 1]) / (sum(as.numeric(VarCorr(fecmodPQLr)[c(2, 4), 1])) + VarOlN)
+# adjusted ICC[Container]
+ICCadjCont <- as.numeric(VarCorr(fecmodPQLf)[4, 1]) / (sum(as.numeric(VarCorr(fecmodPQLf)[c(2, 4), 1])) + VarOlF)
+# comparing the results
+c(
+  R2glmmM = R2glmmM, R2glmmC = R2glmmC, ICCrawPop = ICCrawPop, ICCadjPop = ICCadjPop,
+  ICCrawCont = ICCrawCont, ICCadjCont = ICCadjCont
+)
+
+performance::r2_nakagawa(fecmodPQLf, null_model = fecmodPQLr)
+
+
+# ==============================================================
+# Neg bin with log link (page 9) -------------------------
+# glmmadmb- -----------------------------------------------------
+# ==============================================================
+
+# Fit null model without fixed effects (but including all random effects)
+parmodADMBr <- glmmadmb(
+  Parasite ~ 1 + (1 | Population) + (1 | Container),
+  family = "nbinom2", data = DataAll
+)
+# Fit alternative model including fixed and all random effects
+parmodADMBf <- glmmadmb(
+  Parasite ~ Sex + Treatment + Habitat + (1 | Population) + (1 | Container),
+  family = "nbinom2", data = DataAll
+)
+
+# Calculation of the variance in fitted values
+VarF <- var(as.vector(model.matrix(parmodADMBf) %*% fixef(parmodADMBf)))
+# getting the observation-level variance Null model
+thetaN <- parmodADMBr$alpha # note that theta is called alpha in glmmadmb
+lambda <- as.numeric(exp(fixef(parmodADMBr) + 0.5 * (as.numeric(VarCorr(parmodADMBr)[1]) + as.numeric(VarCorr(parmodADMBr)[2]))))
+# lambda2 <- mean(DataAll$Parasite)
+VarOdN <- 1 / lambda + 1 / thetaN # the delta method
+VarOlN <- log(1 + (1 / lambda) + (1 / thetaN)) # log-normal approximation
+VarOtN <- trigamma((1 / lambda + 1 / thetaN)^(-1)) # trigamma function
+# comparing the three
+c(VarOdN = VarOdN, VarOlN = VarOlN, VarOtN = VarOtN)
+
+thetaF <- parmodADMBf$alpha # note that theta is called alpha in glmmadmb
+VarOdF <- 1 / lambda + 1 / thetaF # the delta method
+VarOlF <- log(1 + (1 / lambda) + (1 / thetaF)) # log-normal approximation
+VarOtF <- trigamma((1 / lambda + 1 / thetaF)^(-1)) # trigamma function
+# comparing the three
+c(VarOdF = VarOdF, VarOlF = VarOlF, VarOtF = VarOtF)
+
+# R2[GLMM(m)] - marginal R2[GLMM]
+R2glmmM <- VarF / (VarF + sum(as.numeric(VarCorr(parmodADMBf))) + VarOlF)
+# R2[GLMM(c)] - conditional R2[GLMM] for full model
+R2glmmC <- (VarF + sum(as.numeric(VarCorr(parmodADMBf)))) / (VarF + sum(as.numeric(VarCorr(parmodADMBf))) + VarOlF)
+# Raw unadjusted ICC[Population]
+ICCrawPop <- VarCorr(parmodADMBr)$Population[1] / (sum(as.numeric(VarCorr(parmodADMBr))) + VarOlN)
+# adjusted ICC[Population]
+ICCadjPop <- VarCorr(parmodADMBf)$Population[1] / (sum(as.numeric(VarCorr(parmodADMBf))) + VarOlF)
+# Raw unadjusted ICC[Container]
+ICCrawCont <- VarCorr(parmodADMBr)$Container[1] / (sum(as.numeric(VarCorr(parmodADMBr))) + VarOlN)
+# adjusted ICC[Container]
+ICCadjCont <- VarCorr(parmodADMBf)$Container[1] / (sum(as.numeric(VarCorr(parmodADMBf))) + VarOlF)
+# comparing the results
+c(
+  R2glmmM = R2glmmM, R2glmmC = R2glmmC, ICCrawPop = ICCrawPop,
+  ICCadjPop = ICCadjPop, ICCrawCont = ICCrawCont, ICCadjCont = ICCadjCont
+)
+
+performance::r2_nakagawa(parmodADMBf, null_model = parmodADMBr)
+
+
+# ==============================================================
+# Neg bin with log link (page 9) -------------------------
+# glmmTMB ------------------------------------------------------
+# ==============================================================
+
+# Fit null model without fixed effects (but including all random effects)
+glmmTMBr <- glmmTMB(
+  Parasite ~ 1 + (1 | Population) + (1 | Container),
+  family = "nbinom2", data = DataAll, REML = TRUE
+)
+# Fit alternative model including fixed and all random effects
+glmmTMBf <- glmmTMB(
+  Parasite ~ Sex + Treatment + Habitat + (1 | Population) + (1 | Container),
+  family = "nbinom2", data = DataAll, REML = TRUE
+)
+
+# Calculation of the variance in fitted values
+VarF <- var(as.vector(get_modelmatrix(glmmTMBf) %*% fixef(glmmTMBf)$cond))
+# getting the observation-level variance Null model
+thetaN <- sigma(glmmTMBr)
+lambda <- as.numeric(exp(fixef(glmmTMBr)$cond + 0.5 * (as.numeric(VarCorr(glmmTMBr)$cond[1]) + as.numeric(VarCorr(glmmTMBr)$cond[2]))))
+# lambda2 <- mean(DataAll$Parasite)
+VarOdN <- 1 / lambda + 1 / thetaN # the delta method
+VarOlN <- log(1 + (1 / lambda) + (1 / thetaN)) # log-normal approximation
+VarOtN <- trigamma((1 / lambda + 1 / thetaN)^(-1)) # trigamma function
+# comparing the three
+c(VarOdN = VarOdN, VarOlN = VarOlN, VarOtN = VarOtN)
+
+thetaF <- sigma(glmmTMBf) # note that theta is called alpha in glmmadmb
+VarOdF <- 1 / lambda + 1 / thetaF # the delta method
+VarOlF <- log(1 + (1 / lambda) + (1 / thetaF)) # log-normal approximation
+VarOtF <- trigamma((1 / lambda + 1 / thetaF)^(-1)) # trigamma function
+# comparing the three
+c(VarOdF = VarOdF, VarOlF = VarOlF, VarOtF = VarOtF)
+
+# R2[GLMM(m)] - marginal R2[GLMM]
+R2glmmM <- VarF / (VarF + sum(as.numeric(VarCorr(glmmTMBf)$cond)) + VarOlF)
+# R2[GLMM(c)] - conditional R2[GLMM] for full model
+R2glmmC <- (VarF + sum(as.numeric(VarCorr(glmmTMBf)$cond))) / (VarF + sum(as.numeric(VarCorr(glmmTMBf)$cond)) + VarOlF)
+c(R2glmmM = R2glmmM, R2glmmC = R2glmmC)
+
+performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr)
+
+
+# ==============================================================
+# Neg bin1 with log link (page 9) -------------------------
+# glmmTMB ------------------------------------------------------
+# ==============================================================
+
+# Fit null model without fixed effects (but including all random effects)
+glmmTMBr <- glmmTMB(
+  Parasite ~ 1 + (1 | Population) + (1 | Container),
+  family = "nbinom1", data = DataAll, REML = TRUE
+)
+# Fit alternative model including fixed and all random effects
+glmmTMBf <- glmmTMB(
+  Parasite ~ Sex + Treatment + Habitat + (1 | Population) + (1 | Container),
+  family = "nbinom1", data = DataAll, REML = TRUE
+)
+
+# Calculation of the variance in fitted values
+VarF <- var(as.vector(get_modelmatrix(glmmTMBf) %*% fixef(glmmTMBf)$cond))
+# getting the observation-level variance Null model
+thetaN <- sigma(glmmTMBr)
+lambda <- as.numeric(exp(fixef(glmmTMBr)$cond + 0.5 * (as.numeric(VarCorr(glmmTMBr)$cond[1]) + as.numeric(VarCorr(glmmTMBr)$cond[2]))))
+# lambda2 <- mean(DataAll$Parasite)
+VarOdN <- 1 / lambda + 1 / thetaN # the delta method
+VarOlN <- log(1 + (1 / lambda) + (1 / thetaN)) # log-normal approximation
+VarOtN <- trigamma((1 / lambda + 1 / thetaN)^(-1)) # trigamma function
+# comparing the three
+c(VarOdN = VarOdN, VarOlN = VarOlN, VarOtN = VarOtN)
+
+thetaF <- sigma(glmmTMBf) # note that theta is called alpha in glmmadmb
+VarOdF <- 1 / lambda + 1 / thetaF # the delta method
+VarOlF <- log(1 + (1 / lambda) + (1 / thetaF)) # log-normal approximation
+VarOtF <- trigamma((1 / lambda + 1 / thetaF)^(-1)) # trigamma function
+# comparing the three
+c(VarOdF = VarOdF, VarOlF = VarOlF, VarOtF = VarOtF)
+
+# R2[GLMM(m)] - marginal R2[GLMM]
+R2glmmM <- VarF / (VarF + sum(as.numeric(VarCorr(glmmTMBf)$cond)) + VarOlF)
+# R2[GLMM(c)] - conditional R2[GLMM] for full model
+R2glmmC <- (VarF + sum(as.numeric(VarCorr(glmmTMBf)$cond))) / (VarF + sum(as.numeric(VarCorr(glmmTMBf)$cond)) + VarOlF)
+c(R2glmmM = R2glmmM, R2glmmC = R2glmmC)
+
+performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr)
+
+
+# ==============================================================
+# Neg bin with log link (page 12) -------------------------
+# glmer.nb ------------------------------------------------------
+# ==============================================================
+
+# Fit null model without fixed effects (but including all random effects)
+parmodGLMERr <- glmer.nb(Parasite ~ (1 | Population) + (1 | Container),
+  data = DataAll
+)
+# Fit alternative model including fixed and all random effects
+parmodGLMERf <- glmer.nb(
+  Parasite ~ Sex + Treatment + Habitat + (1 | Population) + (1 | Container),
+  data = DataAll
+)
+
+# Calculation of the variance in fitted values
+VarF <- var(as.vector(model.matrix(parmodGLMERf) %*% fixef(parmodGLMERf)))
+# getting the observation-level variance Null model
+thetaN <- getME(parmodGLMERr, "glmer.nb.theta")
+lambda <- as.numeric(exp(fixef(parmodGLMERr) + 0.5 * (as.numeric(VarCorr(parmodGLMERr)$Population) + as.numeric(VarCorr(parmodGLMERr)$Container))))
+# lambda2 <- mean(DataAll$Parasite)
+VarOdN <- 1 / lambda + 1 / thetaN # the delta method
+VarOlN <- log(1 + (1 / lambda) + (1 / thetaN)) # log-normal approximation
+VarOtN <- trigamma((1 / lambda + 1 / thetaN)^(-1)) # trigamma function
+# comparing the three
+c(VarOdN = VarOdN, VarOlN = VarOlN, VarOtN = VarOtN)
+# Full model
+thetaF <- getME(parmodGLMERf, "glmer.nb.theta")
+VarOdF <- 1 / lambda + 1 / thetaF # the delta method
+VarOlF <- log(1 + (1 / lambda) + (1 / thetaF)) # log-normal approximation
+VarOtF <- trigamma((1 / lambda + 1 / thetaF)^-1) # trigamma function
+# comparing the three
+c(VarOdF = VarOdF, VarOlF = VarOlF, VarOtF = VarOtF)
+
+R2glmmM <- VarF / (VarF + sum(as.numeric(VarCorr(parmodGLMERf))) + VarOlF)
+R2glmmC <- (VarF + sum(as.numeric(VarCorr(parmodGLMERf)))) / (VarF + sum(as.numeric(VarCorr(parmodGLMERf)) + VarOlF))
+c(R2glmmM = R2glmmM, R2glmmC = R2glmmC)
+
+performance::r2_nakagawa(parmodGLMERf)
+MuMIn::r.squaredGLMM(parmodGLMERf)
+
+
+
+
+
+
+# Fit null model without fixed effects (but including all random effects)
+glmmTMBr <- glmmTMB::glmmTMB(
+  count ~ (1 | site),
+  family = glmmTMB::nbinom1(),
+  data = Salamanders, REML = TRUE
+)
+glmmTMBf <- glmmTMB::glmmTMB(
+  count ~ mined + spp + (1 | site),
+  family = glmmTMB::nbinom1(),
+  data = Salamanders, REML = TRUE
+)
+# Calculation of the variance in fitted values
+VarF <- var(as.vector(get_modelmatrix(glmmTMBf) %*% fixef(glmmTMBf)$cond))
+# getting the observation-level variance Null model
+thetaN <- sigma(glmmTMBr)
+lambda <- as.numeric(exp(fixef(glmmTMBr)$cond + 0.5 * (as.numeric(VarCorr(glmmTMBr)$cond[1]))))
+# lambda2 <- mean(DataAll$Parasite)
+VarOdN <- 1 / lambda + 1 / thetaN # the delta method
+VarOlN <- log(1 + (1 / lambda) + (1 / thetaN)) # log-normal approximation
+VarOtN <- trigamma((1 / lambda + 1 / thetaN)^(-1)) # trigamma function
+# comparing the three
+c(VarOdN = VarOdN, VarOlN = VarOlN, VarOtN = VarOtN)
+
+thetaF <- sigma(glmmTMBf) # note that theta is called alpha in glmmadmb
+VarOdF <- 1 / lambda + 1 / thetaF # the delta method
+VarOlF <- log(1 + (1 / lambda) + (1 / thetaF)) # log-normal approximation
+VarOtF <- trigamma((1 / lambda + 1 / thetaF)^(-1)) # trigamma function
+# comparing the three
+c(VarOdF = VarOdF, VarOlF = VarOlF, VarOtF = VarOtF)
+
+# R2[GLMM(m)] - marginal R2[GLMM]
+R2glmmM <- VarF / (VarF + sum(as.numeric(VarCorr(glmmTMBf)$cond)) + VarOlF)
+# R2[GLMM(c)] - conditional R2[GLMM] for full model
+R2glmmC <- (VarF + sum(as.numeric(VarCorr(glmmTMBf)$cond))) / (VarF + sum(as.numeric(VarCorr(glmmTMBf)$cond)) + VarOlF)
+c(R2glmmM = R2glmmM, R2glmmC = R2glmmC)
+
+performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr)
+MuMIn::r.squaredGLMM(glmmTMBf)
+get_variance(glmmTMBf, null_model = glmmTMBr)
+
+
+
+
+
+sizemodeGLMERr <- glmer(
+  BodyL ~ 1 + (1 | Population) + (1 | Container),
+  family = Gamma(link = log),
+  data = DataAll
+)
+# Fit alternative model including fixed and all random effects
+sizemodeGLMERf <- glmer(
+  BodyL ~ Sex + Treatment + Habitat + (1 | Population) + (1 | Container),
+  family = Gamma(link = log), data = DataAll
+)
+
+
+VarF <- var(as.vector(model.matrix(sizemodeGLMERf) %*% fixef(sizemodeGLMERf)))
+# getting the observation-level variance Null model
+nuN <- 1 / attr(VarCorr(sizemodeGLMERr), "sc")^2 # note that glmer report 1/nu not nu as resiudal varian
+VarOdN <- 1 / nuN # the delta method
+VarOlN <- log(1 + 1 / nuN) # log-normal approximation
+VarOtN <- trigamma(nuN) # trigamma function
+# comparing the three
+c(VarOdN = VarOdN, VarOlN = VarOlN, VarOtN = VarOtN)
+## VarOdN VarOlN VarOtN
+## 0.008370998 0.008336156 0.008406133
+# Full model
+nuF <- 1 / attr(VarCorr(sizemodeGLMERf), "sc")^2 # note that glmer report 1/nu not nu as resiudal varian
+VarOdF <- 1 / nuF # the delta method
+VarOlF <- log(1 + 1 / nuF) # log-normal approximation
+VarOtF <- trigamma(nuF) # trigamma function
+# comparing the three
+c(VarOdF = VarOdF, VarOlF = VarOlF, VarOtF = VarOtF)
+## VarOdF VarOlF VarOtF
+## 0.006750704 0.006728020 0.006773541
+# R2[GLMM(m)] - marginal R2[GLMM]
+R2glmmM <- VarF / (VarF + sum(as.numeric(VarCorr(sizemodeGLMERf))) + VarOlF)
+# R2[GLMM(c)] - conditional R2[GLMM] for full model
+R2glmmC <- (VarF + sum(as.numeric(VarCorr(sizemodeGLMERf)))) / (VarF + sum(as.numeric(VarCorr(sizemodeGLMERf))) + VarOlF)
+c(R2glmmM = R2glmmM, R2glmmC = R2glmmC)
+performance::r2_nakagawa(sizemodeGLMERf, null_model = sizemodeGLMERr)
+
+
+
+
+# Fit null model without fixed effects (but including all random effects)
+parmodCPr <- cpglmm(
+  Parasite ~ 1 + (1 | Population) + (1 | Container),
+  link = 0,
+  data = DataAll
+)
+# Fit alternative model including fixed and all random effects
+parmodCPf <- cpglmm(
+  Parasite ~ Sex + Treatment + Habitat + (1 | Population) + (1 | Container),
+  link = 0,
+  data = DataAll
+)
+
+# Calculation of the variance in fitted values
+VarF <- var(as.vector(model.matrix(parmodCPf) %*% fixef(parmodCPf)))
+phiF <- parmodCPf@phi # the dispersion parameter
+pF <- parmodCPf@p # the index parameter
+VarOdF <- phiF * mu^(pF - 2) # the delta method
+VarOlF <- log(1 + (phiF * mu^(pF - 2)))
+
+R2glmmM <- VarF / (VarF + sum(as.numeric(VarCorr(parmodCPf))) + VarOlF)
+# R2[GLMM(c)] - conditional R2[GLMM] for full model
+R2glmmC <- (VarF + sum(as.numeric(VarCorr(parmodCPf)))) / (VarF + sum(as.numeric(VarCorr(parmodCPf))) + VarOlF)
+
+performance::r2_nakagawa(parmodCPf, null_model = parmodCPr)
+
+model_info(parmodCPf)
+get_sigma(parmodCPf)
+parse(text = safe_deparse(parmodCPf@call))[[1]]$link
+x <- parmodCPf
+x
diff --git a/man/get_variance.Rd b/man/get_variance.Rd
index a9c6434cc..060964885 100644
--- a/man/get_variance.Rd
+++ b/man/get_variance.Rd
@@ -2,6 +2,8 @@
 % Please edit documentation in R/get_variances.R
 \name{get_variance}
 \alias{get_variance}
+\alias{get_variance.merMod}
+\alias{get_variance.glmmTMB}
 \alias{get_variance_residual}
 \alias{get_variance_fixed}
 \alias{get_variance_random}
@@ -13,10 +15,27 @@
 \alias{get_correlation_slopes}
 \title{Get variance components from random effects models}
 \usage{
-get_variance(
+get_variance(x, ...)
+
+\method{get_variance}{merMod}(
+  x,
+  component = c("all", "fixed", "random", "residual", "distribution", "dispersion",
+    "intercept", "slope", "rho01", "rho00"),
+  tolerance = 1e-08,
+  null_model = NULL,
+  approximation = "lognormal",
+  verbose = TRUE,
+  ...
+)
+
+\method{get_variance}{glmmTMB}(
   x,
   component = c("all", "fixed", "random", "residual", "distribution", "dispersion",
     "intercept", "slope", "rho01", "rho00"),
+  model_component = NULL,
+  tolerance = 1e-08,
+  null_model = NULL,
+  approximation = "lognormal",
   verbose = TRUE,
   ...
 )
@@ -25,7 +44,7 @@ get_variance_residual(x, verbose = TRUE, ...)
 
 get_variance_fixed(x, verbose = TRUE, ...)
 
-get_variance_random(x, verbose = TRUE, tolerance = 1e-05, ...)
+get_variance_random(x, verbose = TRUE, tolerance = 1e-08, ...)
 
 get_variance_distribution(x, verbose = TRUE, ...)
 
@@ -42,20 +61,32 @@ get_correlation_slopes(x, verbose = TRUE, ...)
 \arguments{
 \item{x}{A mixed effects model.}
 
+\item{...}{Currently not used.}
+
 \item{component}{Character value, indicating the variance component that should
 be returned. By default, all variance components are returned. The
 distribution-specific (\code{"distribution"}) and residual (\code{"residual"})
 variance are the most computational intensive components, and hence may
 take a few seconds to calculate.}
 
-\item{verbose}{Toggle off warnings.}
-
-\item{...}{Currently not used.}
-
 \item{tolerance}{Tolerance for singularity check of random effects, to decide
 whether to compute random effect variances or not. Indicates up to which
 value the convergence result is accepted. The larger tolerance is, the
 stricter the test will be. See \code{\link[performance:check_singularity]{performance::check_singularity()}}.}
+
+\item{null_model}{Optional, a null-model to be used for the calculation of
+random effect variances. If \code{NULL}, the null-model is computed internally.}
+
+\item{approximation}{Character string, indicating the approximation method
+for the distribution-specific (residual) variance. Only applies to non-Gaussian
+models. Can be \code{"lognormal"} (default), \code{"delta"} or \code{"trigamma"}. For binomial
+models, can also be \code{"observation_level"}. See \emph{Nakagawa et al. 2017} for
+details.}
+
+\item{verbose}{Toggle off warnings.}
+
+\item{model_component}{For models that can have a zero-inflation component,
+specify for which component variances should be returned.}
 }
 \value{
 A list with following elements:
diff --git a/tests/testthat/test-clmm.R b/tests/testthat/test-clmm.R
index 42b1d827e..93d6710c5 100644
--- a/tests/testthat/test-clmm.R
+++ b/tests/testthat/test-clmm.R
@@ -79,12 +79,9 @@ test_that("link_inverse", {
 })
 
 test_that("get_data", {
-  expect_equal(nrow(get_data(m1)), 72)
-  expect_identical(
-    colnames(get_data(m1)),
-    c("rating", "temp", "contact", "judge")
-  )
-  expect_equal(nrow(get_data(m2)), 1847)
+  expect_identical(nrow(get_data(m1)), 72L)
+  expect_named(get_data(m1), c("rating", "temp", "contact", "judge"))
+  expect_identical(nrow(get_data(m2)), 1847L)
   expect_identical(colnames(get_data(m2)), c("SURENESS", "PROD", "RESP"))
 })
 
@@ -137,8 +134,8 @@ test_that("find_terms", {
 })
 
 test_that("n_obs", {
-  expect_equal(n_obs(m1), 72)
-  expect_equal(n_obs(m2), 1847)
+  expect_identical(n_obs(m1), 72)
+  expect_identical(n_obs(m2), 1847)
 })
 
 test_that("linkfun", {
diff --git a/tests/testthat/test-lme.R b/tests/testthat/test-lme.R
index 2dd2129f9..b32d1002e 100644
--- a/tests/testthat/test-lme.R
+++ b/tests/testthat/test-lme.R
@@ -13,12 +13,12 @@ m1 <- nlme::lme(Reaction ~ Days,
 m2 <- nlme::lme(distance ~ age + Sex, data = Orthodont, random = ~1)
 
 set.seed(123)
-sleepstudy$mygrp <- sample(1:5, size = 180, replace = TRUE)
+sleepstudy$mygrp <- sample.int(5, size = 180, replace = TRUE)
 sleepstudy$mysubgrp <- NA
 for (i in 1:5) {
   filter_group <- sleepstudy$mygrp == i
   sleepstudy$mysubgrp[filter_group] <-
-    sample(1:30, size = sum(filter_group), replace = TRUE)
+    sample.int(30, size = sum(filter_group), replace = TRUE)
 }
 
 m3 <- nlme::lme(Reaction ~ Days,
diff --git a/tests/testthat/test-r2_nakagawa_bernoulli.R b/tests/testthat/test-r2_nakagawa_bernoulli.R
new file mode 100644
index 000000000..acbb0620a
--- /dev/null
+++ b/tests/testthat/test-r2_nakagawa_bernoulli.R
@@ -0,0 +1,266 @@
+skip_on_cran()
+
+skip_if_not_installed("glmmTMB")
+skip_if_not_installed("MuMIn")
+skip_if_not_installed("lme4")
+skip_if_not_installed("performance")
+skip_if_not_installed("datawizard")
+
+
+# ==============================================================================
+# Bernoulli mixed models, glmmTMB ----
+# ==============================================================================
+
+test_that("glmmTMB, bernoulli", {
+  # dataset ---------------------------------
+  set.seed(123)
+  dat <- data.frame(
+    outcome = rbinom(n = 500, size = 1, prob = 0.3),
+    var_binom = as.factor(rbinom(n = 500, size = 1, prob = 0.3)),
+    var_cont = rnorm(n = 500, mean = 10, sd = 7)
+  )
+  dat$var_cont <- datawizard::standardize(dat$var_cont)
+  dat$group <- NA
+  dat$group[dat$outcome == 1] <- sample(
+    letters[1:5],
+    size = sum(dat$outcome == 1),
+    replace = TRUE,
+    prob = c(0.1, 0.2, 0.3, 0.1, 0.3)
+  )
+  dat$group[dat$outcome == 0] <- sample(
+    letters[1:5],
+    size = sum(dat$outcome == 0),
+    replace = TRUE,
+    prob = c(0.3, 0.1, 0.1, 0.4, 0.1)
+  )
+
+  # glmmTMB, no random slope -------------------------------------------------
+  m <- glmmTMB::glmmTMB(
+    outcome ~ var_binom + var_cont + (1 | group),
+    data = dat,
+    family = binomial(link = "logit")
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # glmmTMB, probit, no random slope -----------------------------------------
+  m <- glmmTMB::glmmTMB(
+    outcome ~ var_binom + var_cont + (1 | group),
+    data = dat,
+    family = binomial(link = "probit")
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # glmmTMB, cloglog, no random slope -----------------------------------------
+  m <- glmmTMB::glmmTMB(
+    outcome ~ var_binom + var_cont + (1 | group),
+    data = dat,
+    family = binomial(link = "cloglog")
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # glmmTMB, probit, random slope -------------------------------------------------
+  m <- suppressWarnings(glmmTMB::glmmTMB(
+    outcome ~ var_binom + var_cont + (1 + var_cont | group),
+    data = dat,
+    family = binomial(link = "probit")
+  ))
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m, tolerance = 1e-8)
+  # matches theoretical values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # glmmTMB, random slope -------------------------------------------------
+  m <- glmmTMB::glmmTMB(
+    outcome ~ var_binom + var_cont + (1 + var_cont | group),
+    data = dat,
+    family = binomial(link = "logit")
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # glmmTMB, cloglog, random slope -------------------------------------------------
+  m <- glmmTMB::glmmTMB(
+    outcome ~ var_binom + var_cont + (1 + var_cont | group),
+    data = dat,
+    family = binomial(link = "cloglog")
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+})
+
+
+# ==============================================================================
+# Bernoulli mixed models, lme4 ----
+# ==============================================================================
+
+test_that("lme4, bernoulli", {
+  # dataset ---------------------------------
+  set.seed(123)
+  dat <- data.frame(
+    outcome = rbinom(n = 500, size = 1, prob = 0.3),
+    var_binom = as.factor(rbinom(n = 500, size = 1, prob = 0.3)),
+    var_cont = rnorm(n = 500, mean = 10, sd = 7)
+  )
+  dat$var_cont <- datawizard::standardize(dat$var_cont)
+  dat$group <- NA
+  dat$group[dat$outcome == 1] <- sample(
+    letters[1:5],
+    size = sum(dat$outcome == 1),
+    replace = TRUE,
+    prob = c(0.1, 0.2, 0.3, 0.1, 0.3)
+  )
+  dat$group[dat$outcome == 0] <- sample(
+    letters[1:5],
+    size = sum(dat$outcome == 0),
+    replace = TRUE,
+    prob = c(0.3, 0.1, 0.1, 0.4, 0.1)
+  )
+
+  # lme4, no random slope ----------------------------------------------------
+  m <- lme4::glmer(
+    outcome ~ var_binom + var_cont + (1 | group),
+    data = dat,
+    family = binomial(link = "logit")
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # lme4, probit, no random slope ---------------------------------------------
+  m <- lme4::glmer(
+    outcome ~ var_binom + var_cont + (1 | group),
+    data = dat,
+    family = binomial(link = "probit")
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # lme4, cloglog, no random slope ---------------------------------------------
+  m <- lme4::glmer(
+    outcome ~ var_binom + var_cont + (1 | group),
+    data = dat,
+    family = binomial(link = "cloglog")
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # lme4, random slope -------------------------------------------------
+  m <- lme4::glmer(
+    outcome ~ var_binom + var_cont + (1 + var_cont | group),
+    data = dat,
+    family = binomial(link = "logit")
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # lme4, cloglog, random slope -------------------------------------------------
+  m <- lme4::glmer(
+    outcome ~ var_binom + var_cont + (1 + var_cont | group),
+    data = dat,
+    family = binomial(link = "cloglog")
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+})
+
+
+# ==============================================================================
+# Validate against Nakagawa et al. 2017 paper!
+test_that("glmer, Bernoulli", {
+  # example data from Nakagawa et al. 2017
+  Population <- gl(12, 80, 960)
+  Container <- gl(120, 8, 960)
+  Sex <- factor(rep(rep(c("Female", "Male"), each = 8), 60))
+  Habitat <- factor(rep(rep(c("Dry", "Wet"), each = 4), 120))
+  Treatment <- factor(rep(c("Cont", "Exp"), 480))
+  Data <- data.frame(
+    Population = Population, Container = Container, Sex = Sex,
+    Habitat = Habitat, Treatment = Treatment
+  )
+  DataFemale <- Data[Data$Sex == "Female", ]
+  set.seed(777)
+  PopulationE <- rnorm(12, 0, sqrt(0.4))
+  ContainerE <- rnorm(120, 0, sqrt(0.05))
+  EggL <- with(DataFemale, 1.1 + 0.5 * (as.numeric(Treatment) - 1) + 0.1 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container] + rnorm(480, 0, sqrt(0.1)))
+  DataFemale$Egg <- rpois(length(EggL), exp(EggL))
+  DataAll <- Data
+  PopulationE <- rnorm(12, 0, sqrt(0.5))
+  ContainerE <- rnorm(120, 0, sqrt(0.8))
+  ParasiteL <- with(DataAll, 1.8 + 2 * (-1) * (as.numeric(Sex) - 1) + 0.8 * (-1) * (as.numeric(Treatment) - 1) + 0.7 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container])
+  DataAll$Parasite <- rnbinom(length(ParasiteL), size = 5, mu = exp(ParasiteL))
+  PopulationE <- rnorm(12, 0, sqrt(1.3))
+  ContainerE <- rnorm(120, 0, sqrt(0.3))
+  DataAll$BodyL <- 15 + 3 * (-1) * (as.numeric(Sex) - 1) + 0.4 * (as.numeric(Treatment) - 1) + 0.15 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container] + rnorm(960, 0, sqrt(1.2))
+  PopulationE <- rnorm(12, 0, sqrt(0.2))
+  ContainerE <- rnorm(120, 0, sqrt(0.2))
+  ExplorationL <- with(DataAll, 4 + 1 * (-1) * (as.numeric(Sex) - 1) + 2 * (as.numeric(Treatment) - 1) + 0.5 * (-1) * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container])
+  DataAll$Exploration <- rgamma(length(ExplorationL), shape = exp(ExplorationL) * 0.3, rate = 0.3)
+  DataMale <- subset(Data, Sex == "Male")
+  PopulationE <- rnorm(12, 0, sqrt(1.2))
+  ContainerE <- rnorm(120, 0, sqrt(0.2))
+  ColourL <- with(DataMale, 0.8 * (-1) + 0.8 * (as.numeric(Treatment) - 1) + 0.5 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container])
+  DataMale$Colour <- rbinom(length(ColourL), 1, plogis(ColourL))
+
+  morphmodGLMERr <- lme4::glmer(
+    Colour ~ 1 + (1 | Population) + (1 | Container),
+    family = binomial(),
+    data = DataMale
+  )
+  # Fit alternative model including fixed and all random effects
+  morphmodGLMERf <- lme4::glmer(
+    Colour ~ Treatment + Habitat + (1 | Population) + (1 | Container),
+    family = binomial(),
+    data = DataMale
+  )
+
+  VarF <- var(as.vector(model.matrix(morphmodGLMERf) %*% lme4::fixef(morphmodGLMERf)))
+  VarDS <- pi^2 / 3
+  Vt <- lme4::VarCorr(morphmodGLMERr)$Population + lme4::VarCorr(morphmodGLMERr)$Container
+  pmean <- as.numeric(plogis(as.vector(lme4::fixef(morphmodGLMERr)) - 0.5 * Vt * tanh(as.vector(lme4::fixef(morphmodGLMERr)) * (1 + 2 * exp(-0.5 * Vt)) / 6)))
+  VarOL <- 1 / (pmean * (1 - pmean))
+
+  R2glmmM <- VarF / (VarF + sum(as.numeric(lme4::VarCorr(morphmodGLMERf))) + VarDS)
+  R2glmmC <- (VarF + sum(as.numeric(lme4::VarCorr(morphmodGLMERf)))) / (VarF + sum(as.numeric(lme4::VarCorr(morphmodGLMERf))) + VarDS)
+  out <- performance::r2_nakagawa(morphmodGLMERf)
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+
+  R2glmmM <- VarF / (VarF + sum(as.numeric(lme4::VarCorr(morphmodGLMERf))) + VarOL)
+  R2glmmC <- (VarF + sum(as.numeric(lme4::VarCorr(morphmodGLMERf)))) / (VarF + sum(as.numeric(lme4::VarCorr(morphmodGLMERf))) + VarOL)
+  out <- performance::r2_nakagawa(morphmodGLMERf, approximation = "observation_level")
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+})
diff --git a/tests/testthat/test-r2_nakagawa_beta.R b/tests/testthat/test-r2_nakagawa_beta.R
new file mode 100644
index 000000000..ecbd56c6b
--- /dev/null
+++ b/tests/testthat/test-r2_nakagawa_beta.R
@@ -0,0 +1,27 @@
+skip_on_cran()
+
+skip_if_not_installed("glmmTMB")
+skip_if_not_installed("MuMIn")
+skip_if_not_installed("performance")
+
+
+# ==============================================================================
+# beta mixed models, glmmTMB
+# ==============================================================================
+
+skip_if_not_installed("betareg")
+
+test_that("glmmTMB, beta_family", {
+  # dataset ---------------------------------
+  data(FoodExpenditure, package = "betareg")
+  m <- glmmTMB::glmmTMB(
+    I(food / income) ~ income + (1 | persons),
+    data = FoodExpenditure,
+    family = glmmTMB::beta_family()
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- suppressWarnings(performance::r2_nakagawa(m, verbose = FALSE))
+  # matches theoretical values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+})
diff --git a/tests/testthat/test-r2_nakagawa_binomial.R b/tests/testthat/test-r2_nakagawa_binomial.R
new file mode 100644
index 000000000..660dc83e6
--- /dev/null
+++ b/tests/testthat/test-r2_nakagawa_binomial.R
@@ -0,0 +1,50 @@
+skip_on_cran()
+
+skip_if_not_installed("glmmTMB")
+skip_if_not_installed("MuMIn")
+skip_if_not_installed("lme4")
+skip_if_not_installed("performance")
+
+
+# ==============================================================================
+# Binomial mixed models, lme4 ----
+# ==============================================================================
+
+test_that("lme4, binomial", {
+  # dataset
+  data(cbpp, package = "lme4")
+
+  # lme4, no random slope ----------------------------------------------------
+  m <- lme4::glmer(
+    cbind(incidence, size - incidence) ~ period + (1 | herd),
+    data = cbpp,
+    family = binomial()
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-3)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-3)
+})
+
+
+# ==============================================================================
+# Binomial mixed models, glmmTMB ----
+# ==============================================================================
+
+test_that("glmmTMB, binomial", {
+  # dataset
+  data(cbpp, package = "lme4")
+
+  # lme4, no random slope ----------------------------------------------------
+  m <- glmmTMB::glmmTMB(
+    cbind(incidence, size - incidence) ~ period + (1 | herd),
+    data = cbpp,
+    family = binomial()
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-3)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-3)
+})
diff --git a/tests/testthat/test-r2_nakagawa_gamma.R b/tests/testthat/test-r2_nakagawa_gamma.R
new file mode 100644
index 000000000..48185123d
--- /dev/null
+++ b/tests/testthat/test-r2_nakagawa_gamma.R
@@ -0,0 +1,79 @@
+skip_on_cran()
+
+skip_if_not_installed("MuMIn")
+skip_if_not_installed("performance")
+
+# ==============================================================================
+# Gamma mixed models, glmmTMB ----
+# ==============================================================================
+
+# ==============================================================================
+# Validate against Nakagawa et al. 2017 paper!
+test_that("glmmTMB, Gamma", {
+  # example data from Nakagawa et al. 2017
+  Population <- gl(12, 80, 960)
+  Container <- gl(120, 8, 960)
+  Sex <- factor(rep(rep(c("Female", "Male"), each = 8), 60))
+  Habitat <- factor(rep(rep(c("Dry", "Wet"), each = 4), 120))
+  Treatment <- factor(rep(c("Cont", "Exp"), 480))
+  Data <- data.frame(
+    Population = Population, Container = Container, Sex = Sex,
+    Habitat = Habitat, Treatment = Treatment
+  )
+  DataFemale <- Data[Data$Sex == "Female", ]
+  set.seed(777)
+  PopulationE <- rnorm(12, 0, sqrt(0.4))
+  ContainerE <- rnorm(120, 0, sqrt(0.05))
+  EggL <- with(DataFemale, 1.1 + 0.5 * (as.numeric(Treatment) - 1) + 0.1 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container] + rnorm(480, 0, sqrt(0.1)))
+  DataFemale$Egg <- rpois(length(EggL), exp(EggL))
+  DataAll <- Data
+  PopulationE <- rnorm(12, 0, sqrt(0.5))
+  ContainerE <- rnorm(120, 0, sqrt(0.8))
+  ParasiteL <- with(DataAll, 1.8 + 2 * (-1) * (as.numeric(Sex) - 1) + 0.8 * (-1) * (as.numeric(Treatment) - 1) + 0.7 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container])
+  DataAll$Parasite <- rnbinom(length(ParasiteL), size = 5, mu = exp(ParasiteL))
+  PopulationE <- rnorm(12, 0, sqrt(1.3))
+  ContainerE <- rnorm(120, 0, sqrt(0.3))
+  DataAll$BodyL <- 15 + 3 * (-1) * (as.numeric(Sex) - 1) + 0.4 * (as.numeric(Treatment) - 1) + 0.15 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container] + rnorm(960, 0, sqrt(1.2))
+  PopulationE <- rnorm(12, 0, sqrt(0.2))
+  ContainerE <- rnorm(120, 0, sqrt(0.2))
+  ExplorationL <- with(DataAll, 4 + 1 * (-1) * (as.numeric(Sex) - 1) + 2 * (as.numeric(Treatment) - 1) + 0.5 * (-1) * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container])
+  DataAll$Exploration <- rgamma(length(ExplorationL), shape = exp(ExplorationL) * 0.3, rate = 0.3)
+
+  sizemodeGLMERr <- lme4::glmer(
+    BodyL ~ 1 + (1 | Population) + (1 | Container),
+    family = Gamma(link = log),
+    data = DataAll
+  )
+  # Fit alternative model including fixed and all random effects
+  sizemodeGLMERf <- lme4::glmer(
+    BodyL ~ Sex + Treatment + Habitat + (1 | Population) + (1 | Container),
+    family = Gamma(link = log), data = DataAll
+  )
+
+  VarF <- var(as.vector(model.matrix(sizemodeGLMERf) %*% lme4::fixef(sizemodeGLMERf)))
+  nuF <- 1 / attr(lme4::VarCorr(sizemodeGLMERf), "sc")^2
+  VarOdF <- 1 / nuF # the delta method
+  VarOlF <- log(1 + 1 / nuF) # log-normal approximation
+  VarOtF <- trigamma(nuF) # trigamma function
+
+  # lognormal
+  R2glmmM <- VarF / (VarF + sum(as.numeric(lme4::VarCorr(sizemodeGLMERf))) + VarOlF)
+  R2glmmC <- (VarF + sum(as.numeric(lme4::VarCorr(sizemodeGLMERf)))) / (VarF + sum(as.numeric(lme4::VarCorr(sizemodeGLMERf))) + VarOlF)
+  out <- performance::r2_nakagawa(sizemodeGLMERf, null_model = sizemodeGLMERr)
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+
+  # delta
+  R2glmmM <- VarF / (VarF + sum(as.numeric(lme4::VarCorr(sizemodeGLMERf))) + VarOdF)
+  R2glmmC <- (VarF + sum(as.numeric(lme4::VarCorr(sizemodeGLMERf)))) / (VarF + sum(as.numeric(lme4::VarCorr(sizemodeGLMERf))) + VarOdF)
+  out <- performance::r2_nakagawa(sizemodeGLMERf, null_model = sizemodeGLMERr, approximation = "delta")
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+
+  # trigamma
+  R2glmmM <- VarF / (VarF + sum(as.numeric(lme4::VarCorr(sizemodeGLMERf))) + VarOtF)
+  R2glmmC <- (VarF + sum(as.numeric(lme4::VarCorr(sizemodeGLMERf)))) / (VarF + sum(as.numeric(lme4::VarCorr(sizemodeGLMERf))) + VarOtF)
+  out <- performance::r2_nakagawa(sizemodeGLMERf, null_model = sizemodeGLMERr, approximation = "trigamma")
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+})
diff --git a/tests/testthat/test-r2_nakagawa_linear.R b/tests/testthat/test-r2_nakagawa_linear.R
new file mode 100644
index 000000000..d51ba4ea1
--- /dev/null
+++ b/tests/testthat/test-r2_nakagawa_linear.R
@@ -0,0 +1,93 @@
+skip_on_cran()
+
+skip_if_not_installed("glmmTMB")
+skip_if_not_installed("MuMIn")
+skip_if_not_installed("lme4")
+skip_if_not_installed("performance")
+
+
+# ==============================================================================
+# linear mixed models, glmmTMB ----
+# ==============================================================================
+
+test_that("glmmTMB, linear", {
+  data(sleepstudy, package = "lme4")
+
+  # no random effects
+  m1 <- glmmTMB::glmmTMB(Reaction ~ Days, data = sleepstudy)
+  m2 <- lm(Reaction ~ Days, data = sleepstudy)
+  out1 <- performance::r2(m1)
+  out2 <- performance::r2(m2)
+  expect_equal(out1$R2, out2$R2, tolerance = 1e-4)
+  expect_equal(out1$R2_adjusted, out2$R2_adjusted, tolerance = 1e-1)
+
+  # linear, no random slope
+  m <- glmmTMB::glmmTMB(Reaction ~ Days + (1 | Subject), data = sleepstudy)
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  expect_equal(out1[, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # linear, no random slope, inverse
+  m <- suppressWarnings(glmmTMB::glmmTMB(
+    Reaction ~ Days + (1 | Subject),
+    data = sleepstudy,
+    family = gaussian("inverse")
+  ))
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- suppressWarnings(performance::r2_nakagawa(m))
+  expect_equal(out1[, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # linear, with random slope
+  m <- glmmTMB::glmmTMB(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  expect_equal(out1[, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # linear, random slope, inverse
+  m <- suppressWarnings(glmmTMB::glmmTMB(
+    Reaction ~ Days + (1 + Days | Subject),
+    data = sleepstudy,
+    family = gaussian("inverse")
+  ))
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- suppressWarnings(performance::r2_nakagawa(m))
+  expect_equal(out1[, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # linear, random slope, log
+  m <- glmmTMB::glmmTMB(
+    Reaction ~ Days + (1 + Days | Subject),
+    data = sleepstudy,
+    family = gaussian("log")
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  expect_equal(out1[, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+})
+
+
+# ==============================================================================
+# linear mixed models, lme4 ----
+# ==============================================================================
+
+test_that("lme4, linear", {
+  data(sleepstudy, package = "lme4")
+
+  # linear, no random slope
+  m <- lme4::lmer(Reaction ~ Days + (1 | Subject), data = sleepstudy)
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  expect_equal(out1[, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # linear, with random slope
+  m <- lme4::lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  expect_equal(out1[, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+})
diff --git a/tests/testthat/test-r2_nakagawa_negbin.R b/tests/testthat/test-r2_nakagawa_negbin.R
new file mode 100644
index 000000000..587eb7870
--- /dev/null
+++ b/tests/testthat/test-r2_nakagawa_negbin.R
@@ -0,0 +1,208 @@
+skip_on_cran()
+
+skip_if_not_installed("glmmTMB")
+skip_if_not_installed("MuMIn")
+skip_if_not_installed("lme4")
+skip_if_not_installed("performance")
+
+
+# ==============================================================================
+# neg-binomial mixed models, lme4 ----
+# ==============================================================================
+
+test_that("glmer, negbin", {
+  # dataset ---------------------------------
+  data(Salamanders, package = "glmmTMB")
+
+  # no random slopes
+  m <- lme4::glmer.nb(
+    count ~ mined + spp + (1 | site),
+    data = Salamanders
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[2, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[2, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+})
+
+
+# ==============================================================================
+# neg-binomial1 mixed models, glmmTMB ----
+# ==============================================================================
+
+test_that("glmmTMB, Nbinom1", {
+  # we skip this test for now, because MuMIn might use a wrong computation
+  # of the approximation here. See discussion in #877 for details
+  skip_if(TRUE)
+
+  # dataset ---------------------------------
+  data(Salamanders, package = "glmmTMB")
+
+  # no random slopes
+  m <- glmmTMB::glmmTMB(
+    count ~ mined + spp + (1 | site),
+    family = glmmTMB::nbinom1(),
+    data = Salamanders
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[2, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[2, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # with random slopes
+  m <- suppressWarnings(glmmTMB::glmmTMB(
+    count ~ mined + spp + cover + (1 + cover | site),
+    family = glmmTMB::nbinom1(),
+    data = Salamanders
+  ))
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m, tolerance = 1e-8)
+  # matches theoretical values
+  expect_equal(out1[2, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-1)
+  expect_equal(out1[2, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-1)
+
+  # no random slopes, sqrt
+  m <- glmmTMB::glmmTMB(
+    count ~ mined + spp + (1 | site),
+    family = glmmTMB::nbinom1("sqrt"),
+    data = Salamanders
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches delta values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+})
+
+
+# ==============================================================================
+# Validate against Nakagawa et al. 2017 paper!
+test_that("glmmTMB, Nbinom1", {
+  data(Salamanders, package = "glmmTMB")
+  glmmTMBr <- glmmTMB::glmmTMB(
+    count ~ (1 | site),
+    family = glmmTMB::nbinom1(),
+    data = Salamanders, REML = TRUE
+  )
+  glmmTMBf <- glmmTMB::glmmTMB(
+    count ~ mined + spp + (1 | site),
+    family = glmmTMB::nbinom1(),
+    data = Salamanders, REML = TRUE
+  )
+  # Calculation based on Supplement 2 of Nakagawa et al. 2017
+  VarF <- var(as.vector(get_modelmatrix(glmmTMBf) %*% glmmTMB::fixef(glmmTMBf)$cond))
+  # this is "mu" in insight
+  lambda <- as.numeric(exp(glmmTMB::fixef(glmmTMBr)$cond + 0.5 * (as.numeric(glmmTMB::VarCorr(glmmTMBr)$cond[1]))))
+  # this is "sig" in insight
+  thetaF <- sigma(glmmTMBf) # note that theta is called alpha in glmmadmb
+  # this is what ".variance_distributional()" returns
+  VarOdF <- 1 / lambda + 1 / thetaF # the delta method
+  VarOlF <- log(1 + (1 / lambda) + (1 / thetaF)) # log-normal approximation
+  VarOtF <- trigamma((1 / lambda + 1 / thetaF)^-1) # trigamma function
+
+  # lognormal
+  R2glmmM <- VarF / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOlF)
+  R2glmmC <- (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond))) / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOlF)
+  out <- performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr)
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+
+  # delta
+  R2glmmM <- VarF / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOdF)
+  R2glmmC <- (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond))) / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOdF)
+  out <- performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr, approximation = "delta")
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+
+  # trigamma
+  R2glmmM <- VarF / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOtF)
+  R2glmmC <- (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond))) / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOtF)
+  out <- performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr, approximation = "trigamma")
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+})
+
+
+
+# ==============================================================================
+# neg-binomial2 mixed models, glmmTMB
+# ==============================================================================
+
+test_that("glmmTMB, Nbinom2", {
+  # dataset ---------------------------------
+  data(Salamanders, package = "glmmTMB")
+
+  # no random slopes
+  m <- glmmTMB::glmmTMB(
+    count ~ mined + spp + (1 | site),
+    family = glmmTMB::nbinom2(),
+    data = Salamanders
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[2, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[2, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # with random slopes
+  m <- suppressWarnings(glmmTMB::glmmTMB(
+    count ~ mined + spp + cover + (1 + cover | site),
+    family = glmmTMB::nbinom2(),
+    data = Salamanders
+  ))
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m, tolerance = 1e-8)
+  # matches theoretical values
+  expect_equal(out1[2, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-1)
+  expect_equal(out1[2, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-1)
+})
+
+
+# ==============================================================================
+# Validate against Nakagawa et al. 2017 paper!
+test_that("glmmTMB, Nbinom2", {
+  data(Salamanders, package = "glmmTMB")
+  glmmTMBr <- glmmTMB::glmmTMB(
+    count ~ (1 | site),
+    family = glmmTMB::nbinom2(),
+    data = Salamanders, REML = TRUE
+  )
+  glmmTMBf <- glmmTMB::glmmTMB(
+    count ~ mined + spp + (1 | site),
+    family = glmmTMB::nbinom2(),
+    data = Salamanders, REML = TRUE
+  )
+  # Calculation based on Supplement 2 of Nakagawa et al. 2017
+  VarF <- var(as.vector(get_modelmatrix(glmmTMBf) %*% glmmTMB::fixef(glmmTMBf)$cond))
+  # this is "mu" in insight
+  lambda <- as.numeric(exp(glmmTMB::fixef(glmmTMBr)$cond + 0.5 * (as.numeric(glmmTMB::VarCorr(glmmTMBr)$cond[1]))))
+  # this is "sig" in insight
+  thetaF <- sigma(glmmTMBf) # note that theta is called alpha in glmmadmb
+  # this is what ".variance_distributional()" returns
+  VarOdF <- 1 / lambda + 1 / thetaF # the delta method
+  VarOlF <- log(1 + (1 / lambda) + (1 / thetaF)) # log-normal approximation
+  VarOtF <- trigamma((1 / lambda + 1 / thetaF)^-1) # trigamma function
+
+  # lognormal
+  R2glmmM <- VarF / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOlF)
+  R2glmmC <- (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond))) / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOlF)
+  out <- performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr)
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+
+  # delta
+  R2glmmM <- VarF / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOdF)
+  R2glmmC <- (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond))) / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOdF)
+  out <- performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr, approximation = "delta")
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+
+  # trigamma
+  R2glmmM <- VarF / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOtF)
+  R2glmmC <- (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond))) / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOtF)
+  out <- performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr, approximation = "trigamma")
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+})
diff --git a/tests/testthat/test-r2_nakagawa_negbin_zi.R b/tests/testthat/test-r2_nakagawa_negbin_zi.R
new file mode 100644
index 000000000..5bc50ada6
--- /dev/null
+++ b/tests/testthat/test-r2_nakagawa_negbin_zi.R
@@ -0,0 +1,108 @@
+skip_on_cran()
+
+skip_if_not_installed("glmmTMB")
+skip_if_not_installed("MuMIn")
+skip_if_not_installed("performance")
+
+
+# ==============================================================================
+# neg-binomial1 zero-inflated mixed models, glmmTMB
+# ==============================================================================
+
+test_that("glmmTMB, Nbinom1 zero-inflated", {
+  # we skip this test for now, because MuMIn might use a wrong computation
+  # of the approximation here. See discussion in #877 for details
+  skip_if(TRUE)
+
+  # dataset ---------------------------------
+  data(Salamanders, package = "glmmTMB")
+
+  # no random slopes
+  m <- glmmTMB::glmmTMB(
+    count ~ mined + spp + (1 | site),
+    ziformula = ~mined,
+    family = glmmTMB::nbinom1(),
+    data = Salamanders
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[2, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[2, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # with random slopes
+  m <- suppressWarnings(glmmTMB::glmmTMB(
+    count ~ mined + spp + cover + (1 + cover | site),
+    ziformula = ~mined,
+    family = glmmTMB::nbinom1(),
+    data = Salamanders
+  ))
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m, tolerance = 1e-8)
+  # matches theoretical values
+  expect_equal(out1[2, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-1)
+  expect_equal(out1[2, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-1)
+
+  # no random slopes, sqrt
+  m <- glmmTMB::glmmTMB(
+    count ~ mined + spp + (1 | site),
+    ziformula = ~mined,
+    family = glmmTMB::nbinom1("sqrt"),
+    data = Salamanders
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches delta values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+})
+
+
+# ==============================================================================
+# Validate against Nakagawa et al. 2017 paper!
+test_that("glmmTMB, Nbinom1 zero-inflated", {
+  data(Salamanders, package = "glmmTMB")
+  glmmTMBr <- glmmTMB::glmmTMB(
+    count ~ (1 | site),
+    ziformula = ~1,
+    family = glmmTMB::nbinom1(),
+    data = Salamanders, REML = TRUE
+  )
+  glmmTMBf <- glmmTMB::glmmTMB(
+    count ~ mined + spp + (1 | site),
+    ziformula = ~mined,
+    family = glmmTMB::nbinom1(),
+    data = Salamanders, REML = TRUE
+  )
+  # Calculation based on Supplement 2 of Nakagawa et al. 2017
+  VarF <- var(as.vector(get_modelmatrix(glmmTMBf) %*% glmmTMB::fixef(glmmTMBf)$cond))
+  # this is "mu" in insight
+  lambda <- as.numeric(exp(glmmTMB::fixef(glmmTMBr)$cond + 0.5 * (as.numeric(glmmTMB::VarCorr(glmmTMBr)$cond[1]))))
+  # this is "sig" in insight
+  thetaF <- sigma(glmmTMBf) # note that theta is called alpha in glmmadmb
+  # this is what ".variance_distributional()" returns
+  VarOdF <- 1 / lambda + 1 / thetaF # the delta method
+  VarOlF <- log(1 + (1 / lambda) + (1 / thetaF)) # log-normal approximation
+  VarOtF <- trigamma((1 / lambda + 1 / thetaF)^-1) # trigamma function
+
+  # lognormal
+  R2glmmM <- VarF / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOlF)
+  R2glmmC <- (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond))) / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOlF)
+  out <- performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr)
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+
+  # delta
+  R2glmmM <- VarF / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOdF)
+  R2glmmC <- (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond))) / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOdF)
+  out <- performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr, approximation = "delta")
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+
+  # trigamma
+  R2glmmM <- VarF / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOtF)
+  R2glmmC <- (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond))) / (VarF + sum(as.numeric(glmmTMB::VarCorr(glmmTMBf)$cond)) + VarOtF)
+  out <- performance::r2_nakagawa(glmmTMBf, null_model = glmmTMBr, approximation = "trigamma")
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+})
diff --git a/tests/testthat/test-r2_nakagawa_poisson.R b/tests/testthat/test-r2_nakagawa_poisson.R
new file mode 100644
index 000000000..53043de98
--- /dev/null
+++ b/tests/testthat/test-r2_nakagawa_poisson.R
@@ -0,0 +1,131 @@
+skip_on_cran()
+
+skip_if_not_installed("glmmTMB")
+skip_if_not_installed("MuMIn")
+skip_if_not_installed("lme4")
+skip_if_not_installed("performance")
+
+
+# ==============================================================================
+# Poisson mixed models, glmmTMB ----
+# ==============================================================================
+
+test_that("glmmTMB, Poisson", {
+  # dataset ---------------------------------
+  data(Salamanders, package = "glmmTMB")
+
+  # glmmTMB, no random slope -------------------------------------------------
+  m <- glmmTMB::glmmTMB(count ~ mined + (1 | site),
+    family = poisson(), data = Salamanders
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[2, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[2, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # glmmTMB, sqrt, no random slope -------------------------------------------------
+  m <- glmmTMB::glmmTMB(count ~ mined + (1 | site),
+    family = poisson("sqrt"), data = Salamanders
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches delta values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # glmmTMB, random slope -------------------------------------------------
+  m <- suppressWarnings(glmmTMB::glmmTMB(count ~ mined + cover + (1 + cover | site),
+    family = poisson(), data = Salamanders
+  ))
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m, tolerance = 1e-8)
+  # we have slight differences here: MuMIn uses "var(fitted())" to exctract fixed
+  # effects variances, while insight uses "var(beta %*% t(mm))". The latter gives
+  # different values when random slopes are involved
+  expect_equal(out1[2, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-1)
+  expect_equal(out1[2, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-1)
+
+  # glmmTMB, sqrt, random slope -------------------------------------------------
+  m <- suppressWarnings(glmmTMB::glmmTMB(count ~ mined + cover + (1 + cover | site),
+    family = poisson("sqrt"), data = Salamanders
+  ))
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m, tolerance = 1e-8)
+  # matches delta values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+})
+
+
+# ==============================================================================
+# Validate against Nakagawa et al. 2017 paper!
+test_that("glmmTMB, Poisson, Nakagawa", {
+  # dataset ---------------------------------
+  data(Salamanders, package = "glmmTMB")
+
+  # glmmTMB, no random slope -------------------------------------------------
+  fecmodADMBr <- glmmTMB::glmmTMB(count ~ (1 | site),
+    family = poisson(), data = Salamanders
+  )
+  fecmodADMBf <- glmmTMB::glmmTMB(count ~ mined + (1 | site),
+    family = poisson(), data = Salamanders
+  )
+  VarF <- var(as.vector(model.matrix(fecmodADMBf) %*% glmmTMB::fixef(fecmodADMBf)$cond))
+  lambda <- as.numeric(exp(glmmTMB::fixef(fecmodADMBr)$cond + 0.5 * (as.numeric(glmmTMB::VarCorr(fecmodADMBr)$cond))))
+  omegaF <- sigma(fecmodADMBf) # overdispersion omega is alpha in glmmadmb
+  VarOdF <- omegaF / lambda # the delta method
+  VarOlF <- log(1 + omegaF / lambda) # log-normal approximation
+  VarOtF <- trigamma(lambda / omegaF) # trigamma function
+
+  # lognormal
+  R2glmmM <- VarF / (VarF + sum(as.numeric(glmmTMB::VarCorr(fecmodADMBf)$cond)) + VarOlF)
+  R2glmmC <- (VarF + sum(as.numeric(glmmTMB::VarCorr(fecmodADMBf)$cond))) / (VarF + sum(as.numeric(glmmTMB::VarCorr(fecmodADMBf)$cond)) + VarOlF)
+  out <- performance::r2_nakagawa(fecmodADMBf)
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+
+  # delta
+  R2glmmM <- VarF / (VarF + sum(as.numeric(glmmTMB::VarCorr(fecmodADMBf)$cond)) + VarOdF)
+  R2glmmC <- (VarF + sum(as.numeric(glmmTMB::VarCorr(fecmodADMBf)$cond))) / (VarF + sum(as.numeric(glmmTMB::VarCorr(fecmodADMBf)$cond)) + VarOdF)
+  out <- performance::r2_nakagawa(fecmodADMBf, approximation = "delta")
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+
+  # trigamma
+  R2glmmM <- VarF / (VarF + sum(as.numeric(glmmTMB::VarCorr(fecmodADMBf)$cond)) + VarOtF)
+  R2glmmC <- (VarF + sum(as.numeric(glmmTMB::VarCorr(fecmodADMBf)$cond))) / (VarF + sum(as.numeric(glmmTMB::VarCorr(fecmodADMBf)$cond)) + VarOtF)
+  out <- performance::r2_nakagawa(fecmodADMBf, approximation = "trigamma")
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+})
+
+
+# ==============================================================================
+# Poisson mixed models, lme4 ----
+# ==============================================================================
+
+test_that("lme4, Poisson", {
+  # dataset ---------------------------------
+  data(Salamanders, package = "glmmTMB")
+
+  # lme4, no random slope -------------------------------------------------
+  m <- lme4::glmer(count ~ mined + (1 | site),
+    family = poisson(), data = Salamanders
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[2, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[2, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # lme4, sqrt, no random slope -------------------------------------------------
+  m <- lme4::glmer(count ~ mined + (1 | site),
+    family = poisson("sqrt"), data = Salamanders
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+})
diff --git a/tests/testthat/test-r2_nakagawa_poisson_zi.R b/tests/testthat/test-r2_nakagawa_poisson_zi.R
new file mode 100644
index 000000000..529736531
--- /dev/null
+++ b/tests/testthat/test-r2_nakagawa_poisson_zi.R
@@ -0,0 +1,61 @@
+skip_on_cran()
+
+skip_if_not_installed("glmmTMB")
+skip_if_not_installed("MuMIn")
+skip_if_not_installed("performance")
+
+
+# ==============================================================================
+# Poisson zero-inflated mixed models, glmmTMB
+# ==============================================================================
+
+test_that("glmmTMB, Poisson zero-inflated", {
+  # dataset ---------------------------------
+  data(Salamanders, package = "glmmTMB")
+
+  # glmmTMB, no random slope -------------------------------------------------
+  m <- glmmTMB::glmmTMB(count ~ mined + (1 | site),
+    ziformula = ~mined,
+    family = poisson(), data = Salamanders
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches theoretical values
+  expect_equal(out1[2, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[2, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # glmmTMB, sqrt, no random slope -------------------------------------------------
+  m <- glmmTMB::glmmTMB(count ~ mined + (1 | site),
+    ziformula = ~mined,
+    family = poisson("sqrt"), data = Salamanders
+  )
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m)
+  # matches delta values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+
+  # glmmTMB, random slope -------------------------------------------------
+  m <- suppressWarnings(glmmTMB::glmmTMB(count ~ mined + cover + (1 + cover | site),
+    ziformula = ~mined,
+    family = poisson(), data = Salamanders
+  ))
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m, tolerance = 1e-8)
+  # we have slight differences here: MuMIn uses "var(fitted())" to exctract fixed
+  # effects variances, while insight uses "var(beta %*% t(mm))". The latter gives
+  # different values when random slopes are involved
+  expect_equal(out1[2, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-1)
+  expect_equal(out1[2, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-1)
+
+  # glmmTMB, sqrt, random slope -------------------------------------------------
+  m <- suppressWarnings(glmmTMB::glmmTMB(count ~ mined + cover + (1 + cover | site),
+    ziformula = ~mined,
+    family = poisson("sqrt"), data = Salamanders
+  ))
+  out1 <- suppressWarnings(MuMIn::r.squaredGLMM(m))
+  out2 <- performance::r2_nakagawa(m, tolerance = 1e-8)
+  # matches delta values
+  expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
+  expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
+})
diff --git a/tests/testthat/test-r2_nakagawa_tweedie.R b/tests/testthat/test-r2_nakagawa_tweedie.R
new file mode 100644
index 000000000..f004d5948
--- /dev/null
+++ b/tests/testthat/test-r2_nakagawa_tweedie.R
@@ -0,0 +1,98 @@
+skip_on_cran()
+
+skip_if_not_installed("cplm")
+skip_if_not_installed("glmmTMB")
+skip_if_not_installed("performance")
+
+# cplm::cpglmm doesn't work
+suppressPackageStartupMessages({
+  suppressWarnings(suppressMessages(library(cplm, quietly = TRUE, warn.conflicts = FALSE)))
+})
+
+# ==============================================================================
+# Tweedie mixed models, cplm
+# ==============================================================================
+
+test_that("cplm, tweedie", {
+  set.seed(1234)
+  Population <- gl(12, 80, 960)
+  Container <- gl(120, 8, 960)
+  Sex <- factor(rep(rep(c("Female", "Male"), each = 8), 60))
+  Habitat <- factor(rep(rep(c("Dry", "Wet"), each = 4), 120))
+  Treatment <- factor(rep(c("Cont", "Exp"), 480))
+  Data <- data.frame(
+    Population = Population, Container = Container, Sex = Sex,
+    Habitat = Habitat, Treatment = Treatment
+  )
+
+  DataFemale <- Data[Data$Sex == "Female", ]
+  set.seed(777)
+  PopulationE <- rnorm(12, 0, sqrt(0.4))
+  ContainerE <- rnorm(120, 0, sqrt(0.05))
+  EggL <- with(DataFemale, 1.1 + 0.5 * (as.numeric(Treatment) - 1) + 0.1 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container] + rnorm(480, 0, sqrt(0.1)))
+  DataFemale$Egg <- rpois(length(EggL), exp(EggL))
+  DataAll <- Data
+  PopulationE <- rnorm(12, 0, sqrt(0.5))
+  ContainerE <- rnorm(120, 0, sqrt(0.8))
+  ParasiteL <- with(DataAll, 1.8 + 2 * (-1) * (as.numeric(Sex) - 1) + 0.8 * (-1) * (as.numeric(Treatment) - 1) + 0.7 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container])
+  DataAll$Parasite <- rnbinom(length(ParasiteL), size = 5, mu = exp(ParasiteL))
+  PopulationE <- rnorm(12, 0, sqrt(1.3))
+  ContainerE <- rnorm(120, 0, sqrt(0.3))
+  DataAll$BodyL <- 15 + 3 * (-1) * (as.numeric(Sex) - 1) + 0.4 * (as.numeric(Treatment) - 1) + 0.15 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container] + rnorm(960, 0, sqrt(1.2))
+  PopulationE <- rnorm(12, 0, sqrt(0.2))
+  ContainerE <- rnorm(120, 0, sqrt(0.2))
+  ExplorationL <- with(DataAll, 4 + 1 * (-1) * (as.numeric(Sex) - 1) + 2 * (as.numeric(Treatment) - 1) + 0.5 * (-1) * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container])
+  DataAll$Exploration <- rgamma(length(ExplorationL), shape = exp(ExplorationL) * 0.3, rate = 0.3)
+  DataMale <- subset(Data, Sex == "Male")
+  PopulationE <- rnorm(12, 0, sqrt(1.2))
+  ContainerE <- rnorm(120, 0, sqrt(0.2))
+  ColourL <- with(DataMale, 0.8 * (-1) + 0.8 * (as.numeric(Treatment) - 1) + 0.5 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container])
+  DataMale$Colour <- rbinom(length(ColourL), 1, plogis(ColourL))
+
+  parmodCPr <- cpglmm(
+    Parasite ~ 1 + (1 | Population) + (1 | Container),
+    link = 0,
+    data = DataAll
+  )
+  parmodCPf <- cpglmm(
+    Parasite ~ Sex + Treatment + Habitat + (1 | Population) + (1 | Container),
+    link = 0,
+    data = DataAll
+  )
+
+  VarF <- var(as.vector(model.matrix(parmodCPf) %*% cplm::fixef(parmodCPf)))
+  phiF <- parmodCPf@phi # the dispersion parameter
+  pF <- parmodCPf@p # the index parameter
+  mu <- exp(cplm::fixef(parmodCPr) + 0.5 * (cplm::VarCorr(parmodCPr)$Population[1] + cplm::VarCorr(parmodCPr)$Container[1]))
+  VarOdF <- phiF * mu^(pF - 2) # the delta method
+  VarOlF <- log(1 + (phiF * mu^(pF - 2)))
+
+  # lognormal
+  R2glmmM <- VarF / (VarF + sum(as.numeric(cplm::VarCorr(parmodCPf))) + VarOlF)
+  R2glmmC <- (VarF + sum(as.numeric(cplm::VarCorr(parmodCPf)))) / (VarF + sum(as.numeric(cplm::VarCorr(parmodCPf))) + VarOlF)
+  out <- performance::r2_nakagawa(parmodCPf, null_model = parmodCPr)
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+
+  # delta
+  R2glmmM <- VarF / (VarF + sum(as.numeric(cplm::VarCorr(parmodCPf))) + VarOdF)
+  R2glmmC <- (VarF + sum(as.numeric(cplm::VarCorr(parmodCPf)))) / (VarF + sum(as.numeric(cplm::VarCorr(parmodCPf))) + VarOdF)
+  out2 <- performance::r2_nakagawa(parmodCPf, null_model = parmodCPr, approximation = "delta")
+  expect_equal(out2$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out2$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+
+  # results for glmmTMB are very close to cplm
+  m0 <- glmmTMB::glmmTMB(
+    Parasite ~ (1 | Population) + (1 | Container),
+    family = glmmTMB::tweedie(1),
+    data = DataAll
+  )
+  m <- glmmTMB::glmmTMB(
+    Parasite ~ Sex + Treatment + Habitat + (1 | Population) + (1 | Container),
+    family = glmmTMB::tweedie(1),
+    data = DataAll
+  )
+  out3 <- performance::r2_nakagawa(m, null_model = m0)
+  expect_equal(out$R2_conditional, out3$R2_conditional, tolerance = 1e-2, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, out3$R2_marginal, tolerance = 1e-2, ignore_attr = TRUE)
+})
diff --git a/tests/testthat/test-rlmer.R b/tests/testthat/test-rlmer.R
index 3a9b2cc36..764b2af63 100644
--- a/tests/testthat/test-rlmer.R
+++ b/tests/testthat/test-rlmer.R
@@ -267,7 +267,8 @@ test_that("get_variance", {
       var.slope = c(Subject.Days = 41.8070895953001),
       cor.slope_intercept = c(Subject = -0.0387835013909591)
     ),
-    tolerance = 1e-3
+    tolerance = 1e-3,
+    ignore_attr = TRUE
   )
 
   expect_warning(
@@ -287,7 +288,8 @@ test_that("get_variance", {
         mygrp = 16.1137111496375
       )
     ),
-    tolerance = 1e-3
+    tolerance = 1e-3,
+    ignore_attr = TRUE
   )
 })
 
diff --git a/tests/testthat/test-rstanarm.R b/tests/testthat/test-rstanarm.R
index 8d0f9f350..3f34f7108 100644
--- a/tests/testthat/test-rstanarm.R
+++ b/tests/testthat/test-rstanarm.R
@@ -3,6 +3,7 @@ skip_if_offline()
 skip_if_not_installed("lme4")
 skip_if_not_installed("BayesFactor")
 skip_if_not_installed("rstanarm")
+skip_if_not_installed("httr2")
 
 suppressPackageStartupMessages({
   suppressWarnings(suppressMessages(library(rstanarm, quietly = TRUE, warn.conflicts = FALSE)))
@@ -431,8 +432,8 @@ test_that("get_variance", {
     list(
       var.fixed = 0.36274,
       var.random = 0.5988885,
-      var.residual = 3.28987,
-      var.distribution = 3.28987,
+      var.residual = 0.2188036,
+      var.distribution = 0.2188036,
       var.dispersion = 0,
       var.intercept = c(herd = 0.59889)
     ),
@@ -448,11 +449,11 @@ test_that("get_variance", {
     tolerance = 1e-4
   )
   expect_equal(get_variance_residual(m1),
-    c(var.residual = 3.289868),
+    c(var.residual = 0.2188036),
     tolerance = 1e-4
   )
   expect_equal(get_variance_distribution(m1),
-    c(var.distribution = 3.289868),
+    c(var.distribution = 0.2188036),
     tolerance = 1e-4
   )
   expect_equal(get_variance_dispersion(m1),