Why don't the default residuals from sdmTMB match glm() or mgcv::gam()?

[seananderson]

I've seen a couple versions of the following question:

I wanted to compare how a GLM created with stats::glm() compares to the one created with sdmTMB. When I evaluate the predictive performances of these two models (confusion matrix, auc) the results are identical. However, when I plot residual vs. fitted values, the plots are different.

[seananderson]

Answer:

residuals.glm() is returning deviance residuals by default. residuals.sdmTMB() is returning randomized quantile residuals by default. Randomized-quantile residuals are transformed to be a standard normal distribution if the model is consistent with the data. We don't have deviance residuals in sdmTMB currently, although we probably should. The type = "response" should match between the two. An sdmTMB model with spatial = "off" (and spatiotemporal = "off") should match a call to glm() exactly. We have many internal tests to make sure that's the case.

I made the default randomized quantile because I thought they'd be most useful for my work, but that was before I planned on others using the package. So, it's unfortunate that the defaults are different. I may change that eventually or add a message when using the default residuals.

There are many descriptions of randomized quantile residuals out there (also known as probability integral transform [PIT] residuals). This description by Anders Nielsen is pretty good (although the ones in sdmTMB by default are not 'one-step-ahead').

Classic references:

Dunn, P.K. & Smyth, G.K. (1996). Randomized Quantile Residuals. Journal of Computational and Graphical Statistics, 5, 236–244.

Smith, J.Q. (1985). Diagnostic checks of non-standard time series models. Journal of Forecasting, 4, 283–291.

Illustration:

library(sdmTMB)
m1 <- sdmTMB(present ~ 1, family = binomial(), data = pcod, spatial = "off")
m2 <- glm(present ~ 1, family = binomial(), data = pcod)

r1 <- residuals(m1)
r2 <- residuals(m2)

par(mfrow = c(1, 2))
qqnorm(r1, main = "sdmTMB randomized quantile residuals");qqline(r1)
qqnorm(r2, main = "glm deviance residuals");qqline(r2)

r1_resp <- residuals(m1, type = "response")
r2_resp <- residuals(m2, type = "response")
qqnorm(r1_resp, main = "sdmTMB response residuals");qqline(r1_resp)
qqnorm(r2_resp, main = "glm response residuals");qqline(r2_resp)

^{Created on 2023-02-01 with reprex v2.0.2}