wonky plot from `check_model()` on a `glmmTMB` example

[bbolker]

This is from an nbinom1 model - the "overdispersion" and "normality of residuals" plots both look odd ...

library(glmmTMB)
library(dplyr) ## for mutate_at, %>%
#Build example data
x <- c("A", "B", "C", "D")
(time <- rep(x, each=20, times=3)) #time factor
y <- c("exposed", "ref1", "ref2")
(lake <- rep (y, each=80))  #lake factor
set.seed(123)
min <- runif(n=240, min=4.5, max=5.5) #mins used in model offset
set.seed(123)
(count <- rnbinom(n=240,mu=10,size=100)) #randomly generated negative binomial data

#make data frame
dat <- as.data.frame(cbind(time, lake, min, count)) 
dat <- dat %>% 
   mutate_at(c('min', 'count'), as.numeric)

#remove one combination of factors to make example rank deficient (all observations from time A and lake ref1)
dat2 <- filter(dat, time!="A" | lake !="ref1")

model <-glmmTMB(count~time*lake, family=nbinom1,
                      control = glmmTMBControl(rank_check = "adjust"),
                      offset=log(min), data=dat2)

library(performance)
check_model(model)

[strengejacke]

A quick guess is that inappropriate residuals are calculated. This is the code to detect overdispersion for this specific model:

    d <- data.frame(Predicted = stats::predict(model, type = "response"))
    d$Residuals <- insight::get_response(model) - as.vector(d$Predicted)
    d$Res2 <- d$Residuals^2
    d$V <- d$Predicted * (1 + d$Predicted / insight::get_sigma(model))
    d$StdRes <- insight::get_residuals(model, type = "pearson")

And the qq-plot for glmmTMB simply uses stats::residuals(model).

If you don't have a suggestion for calculating the most appropriate residuals, the best solution is probably to go with simulated residuals and diagnostics based on DHARMa (#643)

[bbolker]

OK, the Q-Q plot should definitely be using stats::residuals(model, type = "pearson"). Or residuals(model, type = "deviance") if that's what the y-axis label says ... (we can't yet get standardized deviance residuals, because we still haven't got machinery to extract the hat matrix ...)

I don't know enough about how the columns in the overdispersion machinery are being used downstream or why you need to compute the variance yourself - it seems like it should be possibly to get it more reliably with built-in glmmTMB machinery (but I haven't dug around there ...)

[strengejacke]

Changing to deviance residuals changes the scale of the qq-plot, but not the pattern itself, hm...

[strengejacke]

Ok, found the issue. For binomial or count models, see used a halfnorm distribution, see https://github.com/easystats/see/blob/3e13fc2a3066201191802f6080f5fa42ff5663e4/R/plot.check_normality.R#L170

This was also used for glmmTMB. I now changed the behaviour that glmmTMB models always (no matter which distribution) use the "regular" qq-plot.

Not sure what the reason was for this decision (i.e. using halfnorm) - I can remember we had good reasons, but can't remember details. @easystats/core-team any ideas? @bbolker what would you suggest, the new solution presented below? (as long as we don't rely on DHARMa)

library(glmmTMB)
library(dplyr) ## for mutate_at, %>%
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
# Build example data
x <- c("A", "B", "C", "D")
time <- rep(x, each = 20, times = 3) # time factor
y <- c("exposed", "ref1", "ref2")
lake <- rep(y, each = 80) # lake factor
set.seed(123)
min <- runif(n = 240, min = 4.5, max = 5.5) # mins used in model offset
set.seed(123)
count <- rnbinom(n = 240, mu = 10, size = 100) # randomly generated negative binomial data

# make data frame
dat <- as.data.frame(cbind(time, lake, min, count))
dat <- dat %>%
  mutate_at(c("min", "count"), as.numeric)

# remove one combination of factors to make example rank deficient (all observations from time A and lake ref1)
dat2 <- filter(dat, time != "A" | lake != "ref1")

model <- glmmTMB(count ~ time * lake,
  family = nbinom1,
  control = glmmTMBControl(rank_check = "adjust"),
  offset = log(min), data = dat2
)
#> dropping columns from rank-deficient conditional model: timeD:lakeref1

library(performance)
check_model(model)
#> `check_outliers()` does not yet support models of class `glmmTMB`.

[strengejacke]

I don't know enough about how the columns in the overdispersion machinery are being used downstream or why you need to compute the variance yourself - it seems like it should be possibly to get it more reliably with built-in glmmTMB machinery (but I haven't dug around there ...)

@bwiernik I think you wrote most/all of the code for the overdisperson plot/check.

[bwiernik]

Not sure what the reason was for this decision (i.e. using halfnorm) - I can remember we had good reasons, but can't remember details. @easystats/core-team any ideas? @bbolker what would you suggest, the new solution presented below? (as long as we don't rely on DHARMa)

Base R switched to using half-normnal for binomial and count models. I would suggest we keep using it for glmmTMB for the relevant families?

[bwiernik]

I don't know enough about how the columns in the overdispersion machinery are being used downstream or why you need to compute the variance yourself - it seems like it should be possibly to get it more reliably with built-in glmmTMB machinery (but I haven't dug around there ...)

@bwiernik I think you wrote most/all of the code for the overdisperson plot/check.

I think I just wrote one set of code there and didn't check if anything was already available for glmmTMB models. It would be good to use existing machinery there if possible.

[strengejacke]

I would suggest we keep using it for glmmTMB for the relevant families?

That would be poisson and binomial, but not nbinom?

[strengejacke]

It would be good to use existing machinery there if possible.

Ok - but I'm not sure how to revise the relevant code section. Happy if someone could make a proposal.

[bwiernik]

I would suggest we keep using it for glmmTMB for the relevant families?

That would be poisson and binomial, but not nbinom?

Base R uses a half-normal with absolute standardized deviance residuals for any family of model fit with glm(). We should probably be a little smarter than that and limit it to non-Gaussian models.

If we can't get standardized residuals, then we probably need to adjust the reference distribution to not be a standard normal/half-normal.

[bwiernik]

Let me dig into these plots a bit and see what the most justifiable default should be.

[strengejacke]

or why you need to compute the variance yourself - it seems like it should be possibly to get it more reliably with built-in glmmTMB machinery (but I haven't dug around there ...)

@bbolker we compute the per-observation variance - not sure how to do this with family$variance()?

[bbolker]

Hmm. In principle it would be nice if it were sigma(model)^2*family(model)$variance(predict(model, type = "response")) (at least for models without ZI or dispersion components ...), but sigma.glmmTMB is wonky - returns different values, not necessarily equal to the sqrt of the scaling factor of the $variance() function, for different families (see ?sigma.glmmTMB): maybe there needs to be another function that will get us what we want (or a flag to sigma.glmmTMB ?)

[bwiernik]

Theory on half-normal plot for deviance residuals https://www.jstatsoft.org/article/view/v081i10

[strengejacke]

This looks better for glmmTMB now, but overdispersion plot for glm.nb looks strange to me. Maybe it's correct, though? See #680

library(glmmTMB)
library(performance)
library(MASS)
library(dplyr) ## for mutate_at, %>%
#> 
#> Attaching package: 'dplyr'
#> The following object is masked from 'package:MASS':
#> 
#>     select
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

x <- c("A", "B", "C", "D")
time <- rep(x, each = 20, times = 3) # time factor
y <- c("exposed", "ref1", "ref2")
lake <- rep(y, each = 80) # lake factor
set.seed(123)
min <- runif(n = 240, min = 4.5, max = 5.5) # mins used in model offset
set.seed(123)
count <- rnbinom(n = 240, mu = 10, size = 100) # randomly generated negative binomial data

# make data frame
dat <- as.data.frame(cbind(time, lake, min, count))
dat <- dat %>%
  mutate_at(c("min", "count"), as.numeric)

# remove one combination of factors to make example rank deficient (all observations from time A and lake ref1)
dat2 <- filter(dat, time != "A" | lake != "ref1")

model <- glmmTMB(count ~ time * lake,
  family = nbinom1,
  control = glmmTMBControl(rank_check = "adjust"),
  offset = log(min), data = dat2
)
#> dropping columns from rank-deficient conditional model: timeD:lakeref1
check_model(model)
#> `check_outliers()` does not yet support models of class `glmmTMB`.

m2 <- glm.nb(count ~ time * lake + offset(log(min)), data = dat2)
check_model(m2)

^{Created on 2024-02-05 with reprex v2.1.0}

[bbolker]

We still need to be careful. m3 <- update(model, family = nbinom2); plot(check_overdispersion(m3)) is also busted (precisely because sigma() has the inverse meaning for nbinom2 ... (haven't dug in enough to see what's up with glm.nb)

[strengejacke]

precisely because sigma() has the inverse meaning for nbinom2

We can add a different handling for nbinom2, but what would be the solution in this case? 1/sigma?

[strengejacke]

If the new code is correct, this would be the result.

First, the new code:

  if (faminfo$is_negbin && !faminfo$is_zero_inflated) {
    if (inherits(model, "glmmTMB")) {
      d <- data.frame(Predicted = stats::predict(model, type = "response"))
      d$Residuals <- insight::get_residuals(model, type = "pearson")
      d$Res2 <- d$Residuals^2
      d$StdRes <- insight::get_residuals(model, type = "pearson")
      if (faminfo$family == "nbinom1") {
        # for nbinom1, we can use "sigma()"
        d$V <- insight::get_sigma(model)^2 * stats::family(model)$variance(d$Predicted)
      } else {
        # for nbinom2, "sigma()" has "inverse meaning" (see #654)
        d$V <- (1 / insight::get_sigma(model)^2) * stats::family(model)$variance(d$Predicted)
      }
    } else {
      ## FIXME: this is not correct for glm.nb models?
      d <- data.frame(Predicted = stats::predict(model, type = "response"))
      d$Residuals <- insight::get_response(model) - as.vector(d$Predicted)
      d$Res2 <- d$Residuals^2
      d$V <- d$Predicted * (1 + d$Predicted / insight::get_sigma(model))
      d$StdRes <- insight::get_residuals(model, type = "pearson")
    }
  }

Result:

model <- glmmTMB(count ~ time * lake,
  family = nbinom1,
  control = glmmTMBControl(rank_check = "adjust"),
  offset = log(min), data = dat2
)
m3 <- update(model, family = nbinom2)

p1 <- plot(check_overdispersion(model))
p2 <- plot(check_overdispersion(m3))

p1 + p2

[bbolker]

Looks good (although obviously glm.nb still needs to be fixed ... moving forward I'd like to try to find a way to fix the inconsistency of sigma() definitions (see also https://github.com/glmmTMB/glmmTMB/issues/951). Maybe we could provide some kind of variance function that gave the predicted variance as a function of mu directly (rather than the predicted scaled variance ...)

[strengejacke]

While we're talking about variance, documentation etc.:

This is the variance-function from beta-families:

glmmTMB::beta_family()$variance
#> function (mu) 
#> {
#>     mu * (1 - mu)
#> }
#> <bytecode: 0x000001dc83b3ae58>
#> <environment: 0x000001dc83b3aa30>

The docs, however, say:

Could you clarify which one is correct? I used the one from the docs in insight::get_variance() instead of family$variance(), and if the latter is correct (not the docs), then this could resolve some issues (https://github.com/easystats/insight/issues?q=is%3Aissue+is%3Aopen+label%3Aget_variance)

[bwiernik]

Looks good (although obviously glm.nb still needs to be fixed ... moving forward I'd like to try to find a way to fix the inconsistency of sigma() definitions (see also https://github.com/glmmTMB/glmmTMB/issues/951). Maybe we could provide some kind of variance function that gave the predicted variance as a function of mu directly (rather than the predicted scaled variance ...)

A "variance" option in predict() would be super useful

[bbolker]

That's the thing: the $variance() function as defined/implemented in base R does not return the variance for a given mean, it returns the scaled variance (even though ?family says it is "the variance as a function of the mean"). For example, gaussian$variance is function (mu) rep.int(1, length(mu)). We probably screwed up when we allowed sigma.glmmTMB to return things that are not scale parameters; it should be the case that

sigma(model)^2*family(model)$variance(predict(model, type = "response"))

works generally ...

It could be worth making breaking changes to sigma.glmmTMB to fix this ... although I guess adding a type = "variance" objection to predict would be the most straightforward/least-breaking solution ...

[strengejacke]

I think this is where my lack of statistical expertise prevents me from understanding how to proceed ;-)

The initial code base in insight::get_variance() was drafted by you, and I tried to enhance for further model families. Based on Nakagawa et al. 2017, you commented:

in general want log(1+var(x)/mu^2)

https://github.com/easystats/insight/blob/5f886703cf2303d8f0d0c44b89b29dd092046360/R/compute_variances.R#L603

This is, e.g., what is done for the beta-family, based on the docs (see my screenshot above):

# 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)
  }
}

This sometimes leads to weird results when computing R2 or ICC for mixed models. For most families, performance::r2_nakagawa() is in line with results from other packages. However, "new" families, which are not yet supported by other packages, sometime produce weird results for performance::r2_nakagawa(). Not sure if the distributional variance (calculated in insight:::..compute_variance_distribution() and insight:::.variance_distributional()) is accurate for all model families.

Is there some way to "validate" the results against something? E.g. against non-mixed or Bayesian models, or some simulated data where the R2 is known? (not sure how to simulate such data, though)

[strengejacke]

Initial issue should be fixed, but re-open for the later discussed points here.

[strengejacke]

Q-Q plot now based on DHARMa (see #643), but still need to think about overdispersion plots.

library(glmmTMB)
library(performance)
library(datawizard)
# Build example data
x <- c("A", "B", "C", "D")
time <- rep(x, each = 20, times = 3) # time factor
y <- c("exposed", "ref1", "ref2")
lake <- rep(y, each = 80) # lake factor
set.seed(123)
min <- runif(n = 240, min = 4.5, max = 5.5) # mins used in model offset
set.seed(123)
count <- rnbinom(n = 240, mu = 10, size = 100) # randomly generated negative binomial data

# make data frame
dat <- as.data.frame(cbind(time, lake, min, count))
dat <- dat |>
  data_modify(.at = c("min", "count"), .modify = as.numeric)

# remove one combination of factors to make example rank deficient (all observations from time A and lake ref1)
dat2 <- data_filter(dat, time != "A" | lake != "ref1")

model <- glmmTMB(count ~ time * lake,
  family = nbinom1,
  control = glmmTMBControl(rank_check = "adjust"),
  offset = log(min), data = dat2
)
#> dropping columns from rank-deficient conditional model: timeD:lakeref1

check_model(model)
#> `check_outliers()` does not yet support models of class `glmmTMB`.

^{Created on 2024-03-16 with reprex v2.1.0}

[strengejacke]

@bwiernik Based on my comments here: https://github.com/easystats/performance/pull/643#issuecomment-2001945347

the question is whether we need to do anything regarding the code of the overdispersion plot? The current code relies on residuals / Pearson residuals. To check for over-/underdispersion in more complex models, we now use simulated residuals based on the DHARMa package. Can these residuals possibly be used for the code to create overdispersion plots?

A draft to play with is .new_diag_overdispersion:

https://github.com/easystats/performance/blob/35b5e19988386b584d91116be542baca1e98f33f/R/check_model_diagnostics.R#L296

I'm not sure if this code really works well? I'm not fully understanding the implementation in .diag_overdispersion:

https://github.com/easystats/performance/blob/35b5e19988386b584d91116be542baca1e98f33f/R/check_model_diagnostics.R#L370

and how this "translated" into a function using simulated residuals?

[strengejacke]

Or is the major concern still the variance function and/or sigma?

[bwiernik]

I think we should be able to use the dharma residuals. Let me take a look