Homogeneity-of-variance panel is not variance-standardized for `glmmTMB`/`MixMod` fits

[jbogomolovas2]

Hi, thanks for performance. I'm implementing a Conway-Maxwell-binomial family for
glmmTMB and noticed the Homogeneity of Variance panel shows a pronounced
inverted-U peaking near p = 0.5 for a correctly specified model. I think the
residuals in that panel aren't standardized for glmmTMB objects.

Where

.model_diagnostic_homogeneity() in R/check_model_diagnostics.R branches by
class. The glm branch uses rstandard(model, type = "pearson"), which is
correctly variance-standardized. The glmmTMB/MixMod branch uses
residuals(model) / residual_sigma, where residuals.glmmTMB() defaults to
type = "response" and residual_sigma is a scalar (either
sqrt(get_variance_residual(model)) or .sigma_glmmTMB_nonmixed()). Neither
varies with mu_i.

For non-mixed fits it's stronger than that: .sigma_glmmTMB_nonmixed() returns
1 outright for binomial, poisson and truncated_poisson, so the panel
plots sqrt(|raw response residuals|) with no standardization at all.

This differs from what ?check_model documents under Residuals for
(Generalized) Linear Models, which says the homogeneity plots use standardized
Pearson's residuals for generalized linear models.

Why it bows

For binomial-type families the response-scale SD is proportional to
sqrt(p(1-p)), maximal at p = 0.5. Since the panel plots sqrt(|r_i|),

E sqrt(|r_i|)  ∝  V(mu_i)^(1/4)

so it has to bow whenever fitted probabilities span a reasonable range, however
well specified the model is.

Reproducible example

Same data, same model, two backends. Only the glmmTMB one bows

library(glmmTMB)
library(performance)

set.seed(1)
n <- 600; size <- 20
x <- runif(n, -3, 3)
d <- data.frame(x = x, y = rbinom(n, size, plogis(-0.5 + 1.2 * x)))
d$f <- size - d$y

m_tmb <- glmmTMB(cbind(y, f) ~ x, family = binomial, data = d)
m_glm <- glm(cbind(y, f) ~ x, family = binomial, data = d)

h_tmb <- performance:::.model_diagnostic_homogeneity(m_tmb)
h_glm <- performance:::.model_diagnostic_homogeneity(m_glm)

par(mfrow = c(1, 2))
plot(h_tmb, main = "glmmTMB"); lines(lowess(h_tmb$x, h_tmb$y), lwd = 2, col = 2)
plot(h_glm, main = "glm");     lines(lowess(h_glm$x, h_glm$y), lwd = 2, col = 2)

performance:::.sigma_glmmTMB_nonmixed(m_tmb, insight::model_info(m_tmb))
#> [1] 1

The glm smooth is flat at about 0.85 (roughly E sqrt|N(0,1)|). The glmmTMB
smooth runs 0.17 to 0.27. The y-axes also differ by an order of magnitude, which
is the same point from a different angle.

Suggested fix

Use Pearson residuals in the glmmTMB/MixMod branch, keeping the current
expression as a fallback for families where glmmTMB can't compute them:

} else if (inherits(model, c("glmmTMB", "MixMod"))) {
  r_pearson <- .safe(stats::residuals(model, type = "pearson"))
  if (!is.null(r_pearson)) {
    r_pearson
  } else {
    residual_sigma <- if (faminfo$is_mixed) {
      sqrt(insight::get_variance_residual(model))
    } else {
      .sigma_glmmTMB_nonmixed(model, faminfo)
    }
    stats::residuals(model) / residual_sigma
  }
}

Note this gives Pearson rather than standardized Pearson residuals, so no
sqrt(1 - h_ii) leverage correction.

One incidental note: for families outside the
c("binomial", "poisson", "truncated_poisson") shortcut,
.sigma_glmmTMB_nonmixed() falls through to exp(model$fit$par[["betadisp"]]).
For a family whose second parameter is a dispersion exponent rather than a scale
(mine, and probably others), that isn't a residual SD.

Happy to open a PR if useful.

[strengejacke]

Thanks for this detailed report! It would be great if you could start with a PR, if you have the time.

[Ix3WUhpgTPE]

It seems that the current implementation in .model_diagnostic_homogeneity() for glmmTMB does not correctly account for variance standardization in the residuals, particularly for binomial families. A concrete solution could be to modify the residual calculation for the glmmTMB branch to use a variance standardization method that accounts for the fitted probabilities, perhaps by implementing a new function that calculates Pearson residuals based on the fitted mu_i. Is there a specific strategy or formula you have in mind for this standardization that aligns with how rstandard() handles it for glm fits?