Differences in DHARMa::plotQQunif and check_residuals() for glms

[jebyrnes]

So, I've noticed when plotting simulated quantile residuals, the results can be quite different between using DHARMa and performance.

For example, in this code:

library(DHARMa)
library(lme4)
library(performance)

library(DHARMa)
library(lme4)
library(performance)

set.seed(31415)
testData = createData(sampleSize = 10, overdispersion = 0.5, family = poisson())
fittedModel <- glm(observedResponse ~ Environment1,
                     family = "poisson", data = testData)

simulationOutput <- simulateResiduals(fittedModel = fittedModel)

plotQQunif(simulationOutput,
           testDispersion = FALSE,
           testUniformity = FALSE,
           testOutliers = FALSE)

check_residuals(simulationOutput) |> plot()

This produces

While qualitatively these match, there are some big differences at the upper tail with the performance version looking really problematic while the DHARMa version does not.

What's going on here? This seems like there shouldn't be disagreement, and it makes me nervous that there is!

This is a reprex, but in working with real data, the discrepencies are even larger - although the KS tests agree. So, the values are the same? I'm confused.

[strengejacke]

Well, the patterns look very identical?

[strengejacke]

I see that in the lower plot, the line does not start at (0|0)?
@mattansb we use qqplotr in see::plot.see_performance_simres.R, is that package inaccurate, or does it somehow slightly "rotate" the line?

library(DHARMa)
#> This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(performance)

data(Salamanders, package = "glmmTMB")

m <- glm(count ~ mined + spp, family = "poisson", data = Salamanders)

simulationOutput <- simulateResiduals(fittedModel = m)

plotQQunif(simulationOutput,
           testDispersion = FALSE,
           testUniformity = FALSE,
           testOutliers = FALSE)

check_residuals(simulationOutput) |> plot()

^{Created on 2025-03-26 with reprex v2.1.1}

[strengejacke]

Here is the detrended plot, the lower left points now no longer start at (0|0)...

check_residuals(simulationOutput) |> 
  plot(detrend = TRUE) +
  ggplot2::ylim(c(-1, 1))

[mattansb]

This seems to be controlled by the identity= argument in {qqplotr} to identity = TRUE (I have no idea what that does or if it should always be used?)

library(DHARMa)
library(performance)

library(ggplot2)
library(qqplotr)

data(Salamanders, package = "glmmTMB")

m <- glm(count ~ mined + spp, family = "poisson", data = Salamanders)

simulationOutput <- simulateResiduals(fittedModel = m)

plotQQunif(simulationOutput,
           testDispersion = FALSE,
           testUniformity = FALSE,
           testOutliers = FALSE)

ggplot(mapping = aes(sample = simulationOutput$scaledResiduals)) + 
  stat_qq_band(distribution = "unif", identity = TRUE) + 
  stat_qq_line(distribution = "unif", identity = TRUE) + 
  stat_qq_point(distribution = "unif", identity = TRUE)

[strengejacke]

Ah, nice!
The docs say

identity: Logical. Should an identity line be used as the reference line? If TRUE, the identity line is used. If FALSE (default), the commonly-used Q-Q line that intercepts two data quantiles specified in qprobs is used. Please notice that the chosen reference line will also be used for the detrending procedure, if detrend = TRUE.

If I understand, basically it's about whether the line crossed min/max or q25%/75%.

So our approach is the one "commonly used"?

[mattansb]

I guess? 🤷‍♂️

[jebyrnes]

Looks like somehow gap::qqunif - which is what DHARMa uses - is producing a different plot.

library(DHARMa)
library(performance)

data(Salamanders, package = "glmmTMB")

m <- glm(count ~ mined + spp, family = "poisson", data = Salamanders)

a <- simulateResiduals(fittedModel = m)


ggplot(mapping =aes(sample = a$scaledResiduals)) +
  stat_qq_line(distribution  = "unif", identity = FALSE) + 
  stat_qq_point(distribution  = "unif", identity = FALSE)

gap::qqunif(a$scaledResiduals, logscale = FALSE)

versus

[strengejacke]

@florianhartig any thoughts? Which would you consider more appropriate, or is it a matter of taste?

[strengejacke]

(may I re-ping @florianhartig again? 😬 )

[florianhartig]

Sorry for the late reply, hadn't seen this before. @melina-leite, maybe this discussion is also interesting for (just to be informed)

I had never thought about how to define the line, I'm was just using the defaults of gap::qqunif . However, reading the discussion on https://stats.stackexchange.com/questions/601731/how-to-define-the-line-to-fit-in-q-q-plot , it seems to me that the DHARMa default with intercept 0, slope 1 is the best choice for testing against uniform as there shouldn't be any outlier issue etc.

Cheers
Florian