Use Pearson residuals for glmmTMB/MixMod heteroscedasticity plot

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: see
 Title: Model Visualisation Toolbox for 'easystats' and 'ggplot2'
-Version: 0.14.1.2
+Version: 0.14.1.3
 Authors@R:
     c(person(given = "Daniel",
              family = "Lüdecke",
@@ -43,6 +43,10 @@ Authors@R:
              role = "ctb",
              email = "jeffrey.r.stevens@gmail.com",
              comment = c(ORCID = "0000-0003-2375-1360")),
+      person(given = "Julius",
+             family = "Bogomolovas",
+             role = "ctb",
+             comment = c(ORCID = "0000-0001-8344-1909")),
       person(given = "Matthew",
              family = "Smith",
              role = "rev",

NEWS.md

@@ -5,6 +5,15 @@
 * Argument `linewidth` was renamed to `size_line`, to be consistent across the
   easystats-ecosystem.
 
+## Bug fixes
+
+* `plot()` for `check_heteroscedasticity()` now uses Pearson residuals for
+  `glmmTMB` and `MixMod` models. Previously these residuals were divided by a
+  single scalar, which is only correct when the variance function does not
+  depend on the mean; for non-mixed binomial and Poisson models that scalar was
+  1, so the plot could suggest heteroscedasticity for correctly specified
+  models. See easystats/performance#926.
+
 # see 0.14.1
 
 ## Changes

R/plot.check_heteroscedasticity.R

@@ -39,12 +39,25 @@ plot.see_check_heteroscedasticity <- function(
     if (inherits(model, "merMod")) {
       stats::residuals(model, scaled = TRUE)
     } else if (inherits(model, c("glmmTMB", "MixMod"))) {
-      sig <- if (faminfo$is_mixed) {
-        sqrt(insight::get_variance_residual(model))
+      ## Pearson residuals are scaled by the family's variance function V(mu_i),
+      ## which varies across observations. The fallback below divides by a single
+      ## scalar, which is only correct when V() does not depend on mu (e.g.
+      ## gaussian). For non-mixed binomial/poisson models `.sigma_glmmTMB_nonmixed()`
+      ## returns 1, i.e. no standardization at all.
+      r_pearson <- tryCatch(
+        stats::residuals(model, type = "pearson"),
+        error = function(e) NULL
+      )
+      if (is.null(r_pearson) || all(is.na(r_pearson))) {
+        sig <- if (faminfo$is_mixed) {
+          sqrt(insight::get_variance_residual(model))
+        } else {
+          .sigma_glmmTMB_nonmixed(model, faminfo)
+        }
+        stats::residuals(model, type = "response") / sig
       } else {
-        .sigma_glmmTMB_nonmixed(model, faminfo)
+        r_pearson
       }
-      stats::residuals(model) / sig
     } else if (inherits(model, "glm")) {
       stats::rstandard(model, type = "pearson")
     } else {