# Goodness-of-fit glm: Pearson's residuals or deviance residuals?

I want to evaluate the goodness-of-fit (or badness-of-fit) of a negative binomial glm. However, even here within CV, I've seen multiple different approaches for doing so.

Some use the the residual deviance (here, and second answer here), some don't specify which deviance to use (otherwise nice answer here), still others emphasize that you really ought to use the Pearson's residuals (see pg. 13 of these great lecture notes here, see Zuur et al. 2009, see this post, see this post, see update to this post)

Can someone confirm that I ought to be using the Pearson's residuals -- not deviance residuals -- to get a sense of model fit?

EP <- resid(model.nb.step, type = "pearson")
ED <- resid(model.nb.step, type = "deviance")
sum(EP^2) # 59.05
sum(ED^2) # 56.84
model.nb.step\$deviance # 56.84041 (same as above)

## I think this is the *improper* way to do it:
pchisq(model.nb.step$$deviance, df=model.nb.step$$df.residual, lower.tail=FALSE)
# 0.1540072
pchisq(sum(ED^2), df = (51 - 4), lower.tail=FALSE)
# 0.1540072 (same as above)

## I think this is the *proper* way to do it:
pchisq(sum(EP^2), df = (51 - 4), lower.tail=FALSE)
# 0.1115992


with the null being that the model is correctly specified such that p > 0.05 suggest I can retain my model. The results from the two approaches aren't that different, but I want to make sure I'm doing this properly. Thanks!

Model below.

Call:
glm.nb(formula = tally ~ slp100 + shrub_perc + min_live_seed_dist +
offset(log(area)), data = data, maxit = 100, init.theta = 1.318046566,

Deviance Residuals:
Min       1Q   Median       3Q      Max
-2.5082  -0.8702  -0.3185   0.4379   2.6724

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)          7.9815     0.1228  64.984  < 2e-16 ***
slp100              -0.2337     0.1321  -1.770 0.076735 .
shrub_perc          -0.4206     0.1277  -3.294 0.000987 ***
min_live_seed_dist  -0.3991     0.1336  -2.988 0.002806 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(1.318) family taken to be 1)

Null deviance: 74.979  on 50  degrees of freedom
Residual deviance: 56.840  on 47  degrees of freedom
AIC: 589.34

Number of Fisher Scoring iterations: 1

Theta:  1.318
Std. Err.:  0.241

2 x log-likelihood:  -579.345

• I don't think there is a clear right or wrong type of residual. Both types of residuals are used in the literature (see for example Chapter 6.4 "GLM Diagnostics" in Faraway's Extending the Linear Model with R, 2006). Based on your output it also shows that they are indeed not much different. Whichever you decide to use just report which ones you used. And why not reporting both? The more you know about your model the better. Commented Mar 11, 2019 at 16:24
• @Stefan, yes I'm quite aware that both are used in the literature, but I'm asking specifically about applying them in a goodness-of-fit test. There seems to be some strong feelings about deviance residuals not being as good as Pearson's residuals for evaluating fit -- the former perhaps does not approximate X^2 as well. See p. 12 of these lecture notes for an explanation (as well as other links I included above). I'm just seeking more consensus. Commented Mar 11, 2019 at 17:17