How does R find the dispersion parameter in a GLM?

Question

I'm working on a problem involving fitting a GLM to data and I'm curious about how R calculates the dispersion parameter. For example, I have this output for the summary of my GLM.

glm(formula = Lifespan ~ glucose + Temperature, family = 
              Gamma(link = "inverse"), 
data = dat)

Deviance Residuals: 
Min       1Q   Median       3Q      Max  
-2.2337  -0.7800  -0.1906   0.3331   2.1397  

Coefficients:
         Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.455352   0.048656   9.359  < 2e-16 ***
glucose     -0.066062   0.015573  -4.242 3.41e-05 ***
Temperature -0.007778   0.001248  -6.233 2.73e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 0.6810965)

Null deviance: 172.16  on 199  degrees of freedom
Residual deviance: 135.52  on 197  degrees of freedom
AIC: 1207.6

Number of Fisher Scoring iterations: 6

My question is, how does R calculate the dispersion parameter to be equal to 0.6810965? I've tried looking at the documentation for the GLM function but I can't seem to find it.

As a follow-up question too, I tried changing the summary so that the dispersion is one, and this doesn't affect the AIC value or the null/residual deviances, so how can I tell if my model fits better with dispersion 1 or the given dispersion in the output, 0.68?

Answer 1

You can look in the code for summary.glm where you'll see:

sum((object$weights * object$residuals^2)[object$weights > 0])/df.r

• df.r is the "residual degrees of freedom" (number of observations - number of parameters)
• object$residuals (from ?glm):

residuals: the working residuals, that is the residuals in the final iteration of the IWLS fit. Since cases with zero weights are omitted, their working residuals are ‘NA’.

I don't have the definition of the working residuals at the top of my brain right now, but it turns out that this is equivalent to calculating the sum of the squared Pearson residuals ($\sum (Y_i-\mu_i)^2/v_i$, where $v_i$ is the scaled variance predicted by the model, equal to $\mu_i^2$ for the Gamma family) divided by the residual df.

Running example("glm") to get the fitted object glm.D93:

> sum(residuals(glm.D93, "pearson")^2)
[1] 5.173202
> with(glm.D93, sum(weights*residuals^2))
[1] 5.173202

This is basically the sum of squares of the residuals (scaled to make the residuals homoscedastic), divided by the residual df (so it matches up with the way you estimate the standard deviation of a (Gaussian) linear model (mean squared error / residual df)).

One thing to keep in mind is that most GLM estimates of dispersion are approximations (and that even if you went to the trouble of getting a maximum likelihood estimate it would be biased for small sample size). For example, you can also compute sigma(glm.D93), which divides the deviance (or equivalently the sum of squares of the deviance residuals) by the residual df. Venables and Ripley MASS have a comment about the unreliability of these approximations for small sample sizes.

I'm not clear that your question ("does my model fit better with dispersion 1 or the estimated dispersion") makes sense; changing the dispersion does not change the predictions of the model at all.