Skip to content

Use Pearson residuals for glmmTMB/MixMod heteroscedasticity plot - #460

Open
jbogomolovas2 wants to merge 1 commit into
easystats:mainfrom
jbogomolovas2:fix-glmmtmb-homogeneity
Open

Use Pearson residuals for glmmTMB/MixMod heteroscedasticity plot#460
jbogomolovas2 wants to merge 1 commit into
easystats:mainfrom
jbogomolovas2:fix-glmmtmb-homogeneity

Conversation

@jbogomolovas2

Copy link
Copy Markdown

Companion to easystats/performance#927, which @strengejacke asked me to open.

The problem

plot.see_check_heteroscedasticity() carries its own copy of the branch fixed
in performance#927, including a private copy of .sigma_glmmTMB_nonmixed(). So
plot(check_heteroscedasticity(m)) is affected the same way: for glmmTMB and
MixMod models the residuals are divided by a single scalar, which is only
correct when V() does not depend on mu. For non-mixed binomial and Poisson
models that scalar is 1, so no standardization happens and the panel bows for
correctly specified models.

The full write-up, including the per-family table showing this affects nine of
eleven glmmTMB families rather than just binomial, is in performance#927.

The change

Same fix, ported. Try residuals(type = "pearson"), keep the existing
expression as a fallback. see has no .safe() helper, so this uses
tryCatch() directly, matching the surrounding code.

Effect

Correctly specified glmmTMB binomial fit (n = 600, size = 20). Left: patched.
Right: current behaviour, where the trend line humps in the middle purely
because binomial spread is largest near p = 0.5.
see

Also worth flagging,

.sigma_glmmTMB_nonmixed() in R/plot.check_heteroscedasticity.R reads

betad <- model$fit$par["betad"]

glmmTMB renamed that parameter to betadisp (glmmTMB commit 472ccbe0);
performance's copy of this helper was updated, this one was not.

[ does not partial-match, so the lookup returns NA, exp(NA) is NA, and
dividing by it makes the whole panel NA. Nothing errors, so the surrounding
tryCatch does not catch it and the plot comes out silently blank.

Affects gaussian and Gamma glmmTMB fits. binomial/poisson/truncated_poisson
short-circuit to 1 before reaching that line.

Confirmed on glmmTMB 1.1.15:

names(glmmTMB(y ~ x, family = gaussian, data = d)$fit$par)
#> "beta" "betadisp"    -- no "betad"

Fix is the rename, in all three places in the switch(). Happy to PR it, but
noting it separately from # since it is a different failure mode.

Notes

  • Version bumped and NEWS item added, following the same convention as
    performance#927.
  • Added myself to DESCRIPTION.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant