Fix overdispersion plot for glmmTMB and mixed models: correct residual scale and family variance formulas

[Copilot]

check_model() produced nonsensical overdispersion plots for glmmTMB models because of two compounding bugs: incorrect residual scaling and the wrong family variance formula being applied to nbinom1 models.

Changes

• .model_diagnostic_overdispersion() (R/check_model_diagnostics.R): For glmmTMB and mixed models (when DHARMa is available), uses DHARMa to obtain predictions and then computes response-scale residuals (observedResponse - fittedPredictedResponse):

  • Res2 = squared response-scale residuals — same scale as the family variance V
  • V = family-specific expected variance from .expected_variance() (e.g. μ(1+φ) for nbinom1)
  • StdRes = residuals / sqrt(V) (Pearson-like standardized residuals for the type-2 plot)
  • Falls back to existing Pearson-residual path if DHARMa is unavailable or simulation fails

• .expected_variance() (R/check_model_diagnostics.R): Refactored to correctly dispatch on glmmTMB family names before the general is_negbin branch. Previously, the general is_negbin branch caught nbinom1 models and applied the nbinom2 formula (μ(1 + μ/φ)). Formulas now match the glmmTMB documentation:

  • nbinom1: V = μ(1 + φ)
  • nbinom2: V = μ(1 + μ/φ)
  • Zero-inflated variants use "zprob" (glmmTMB) or "zero" (other packages) correctly

• Removed .new_diag_overdispersion() — unused prior attempt that mixed DHARMa residuals with the old variance scale, producing an inconsistent plot

• R/check_model.R: Pass ... through to .model_diagnostic_overdispersion() so user-provided DHARMa args (e.g. iterations) propagate correctly

model <- glmmTMB(count ~ time * lake, family = nbinom1,
                 control = glmmTMBControl(rank_check = "adjust"),
                 offset = log(min), data = dat2)
check_model(model)  # overdispersion panel now shows correctly scaled Res2 vs V

[strengejacke]

@copilot I think the residuals using DHARMa are ok now, but I'm not sure about the variance. I refactored the code, and the variance is now calculated in .expected_variance(). I think the variance should match the documentation of the families in the glmmTMB package, here is a screenshot:

And you check and - if needed - update the calculation of the expected variance in .expected_variance()?

[Copilot]

@copilot I think the residuals using DHARMa are ok now, but I'm not sure about the variance. I refactored the code, and the variance is now calculated in .expected_variance(). I think the variance s...

Done in commit 178a471. Two issues were fixed:

1. Scale mismatch in DHARMa path: The previous code used qnorm(DHARMa_uniform_residuals) giving Res2 ≈ 1, but V from .expected_variance() was the raw family variance (~500 for nbinom1). The DHARMa path now uses response-scale residuals (observedResponse - fittedPredictedResponse) so both Res2 and V are on the same scale. StdRes is now residuals / sqrt(V) (Pearson-like).

2. Wrong formula for glmmTMB nbinom1: The else if (faminfo$is_negbin && !faminfo$is_zero_inflated) branch in .expected_variance() was catching glmmTMB nbinom1 models before they reached the correct faminfo$family == "nbinom1" check, so the nbinom2 formula μ(1 + μ/φ) was applied instead of the correct μ(1 + φ). The function is now restructured so that the explicit nbinom1/nbinom2 family checks come first, matching the glmmTMB documentation: nbinom1: V=μ(1+φ), nbinom2: V=μ(1+μ/φ).