Does this comparison of the quasipoisson model to the poisson model make sense?

Question

I'm messing around with the use of the quasi likelihood method so this is my first attempt. I think I'm misunderstanding a key component, since I do understand why the parameter estimates are the same (as they both maximise the same log likelihood function) and the standard errors are different because of the overdispersion however, why is the deviance the same?

When my professor saw my code, he said you cannot really compare these as they aren't 'nested' but he also said to look at the residual sum of squares. The code is as follows:

> library(ggplot2)
> set.seed(42)
> 
> n <- 100
> x <- runif(n, 0, 10)
> beta_0 <- 1.0
> beta_1 <- 0.5
> mu <- exp(beta_0 + beta_1 * x)
> phi <- 2
> y_overdispersed <- rnbinom(n, size = phi, mu = mu)
> data <- data.frame(x = x, y = y_overdispersed)
> 
> model_quasi <- glm(y ~ x, family = quasipoisson(link = "log"), 
                     data = data)
> data$pred_quasi <- predict(model_quasi, type = "response")
> 
> model_poisson <- glm(y ~ x, family = poisson(link = "log"), 
                   data = data)
> data$pred_poisson <- predict(model_poisson, type = "response")
> 
> cat("Quasi-Poisson Model Summary:\n")
Quasi-Poisson Model Summary:
> summary(model_quasi)

Call:
glm(formula = y ~ x, family = quasipoisson(link = "log"), 
         data = data)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.81145    0.30410   2.668  0.00892 ** 
x            0.50345    0.03593  14.012  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 29.85742)

    Null deviance: 12464.7  on 99  degrees of freedom
Residual deviance:  2696.6  on 98  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

> 
> cat("\nRegular Poisson Model Summary:\n")

Regular Poisson Model Summary:
> summary(model_poisson)

Call:
glm(formula = y ~ x, family = poisson(link = "log"), data = data)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) 0.811453   0.055653   14.58   <2e-16 ***
x           0.503446   0.006575   76.56   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 12464.7  on 99  degrees of freedom
Residual deviance:  2696.6  on 98  degrees of freedom
AIC: 3203.9

Number of Fisher Scoring iterations: 5

Answer 1

The deviance is the same because the code forces it to be the same. The deviance is computed in stats::glm.fit by

dev <- sum(dev.resids(y, mu, weights))

and the dev.resids function is the same for poisson() and quasipoisson() families.

> poisson()$dev.resids
function (y, mu, wt) 
{
    r <- mu * wt
    p <- which(y > 0)
    r[p] <- (wt * (y * log(y/mu) - (y - mu)))[p]
    2 * r
}
> quasipoisson()$dev.resids
function (y, mu, wt) 
{
    r <- mu * wt
    p <- which(y > 0)
    r[p] <- (wt * (y * log(y/mu) - (y - mu)))[p]
    2 * r
}

I think this is how quasilikelihood deviances were defined by McCullagh & Nelder, but I don't have the book at hand to check.

Quasilikelihood-ratio tests involve dividing the deviance difference by the dispersion parameter, and there will be a difference at that point.