How to choose shape parameter for gamma distribution when using generalized linear model?

Question

I guess we can figure out coefficients without setting specific shape parameter $\alpha$ for gamma distribution because its effect is cancelled out when we use IRLS:

$$\beta^{(m)}=(X^TWX)^{-1}X^TWz^{(m)}$$

Here, the $i$th diagonal element of weight matrix is $W_{ii}=\frac{\alpha}{\mu_i^2 g'(\mu_i)}$ and the $g$ is the link function. So if we write $W=\alpha Q$, the equation above will be $\beta^{(m)}=(X^TQX)^{-1}X^TQz^{(m)}$ and we can do IRLS not knowing $\alpha$.

My question is that, after we find and calculate the coefficients, which value of $\alpha$ should we choose?

In linear regression case, we could estimate error by using sum of squre of errors. I think there should be a similar thing (like,maybe function of deviance?) to estimate the shape parameter.

It seems like the function glm in R shows dispersion parameter which is reciprocal of the shape parameter automatically. (Dispersion parameter for Gamma family). I wonder how it works.

Answer 1

What R does by default is to compute the dispersion as the residual deviance divided by the residual degrees of freedom. The deviance is analogous to the residual sum of squares (so the dispersion is analogous to the residual standard deviation); as in the linked question, the shape parameter for the Gamma (1/CV^2, or the reciprocal of the variance scaled by the square of the mean) is the reciprocal of the dispersion. MASS::gamma.shape() improves on this estimate:

A glm fit for a Gamma family correctly calculates the maximum likelihood estimate of the mean parameters but provides only a crude estimate of the dispersion parameter. This function takes the results of the glm fit and solves the maximum likelihood equation for the reciprocal of the dispersion parameter, which is usually called the shape (or exponent) parameter.

This should be described in more detail in Venables and Ripley's book (Modern Applied Statistics with S), but a quick Google Books search can't find it (my copy is not in the same place I am right now). If you're brave you can look at the code of MASS:::gamma.shape.glm to see what it's actually doing.

Answer 2

Thanks to Ben Bolker, I have searched based on his answer. And I would like to add it here for someone else like me studying generalized linear model. The method

to compute the dispersion as the residual deviance divided by the residual degrees of freedom

is called Deviance method and it turned out there are other methods to estimate the dispersion parameter. One of them is Pearson method using the sum of the squares of the residuals divided by the residual degrees of freedom.

Here I add some source code to run in R.

set.seed(1)
#making artificial data
x<-c(1:20)
eta<-1-.5*x
y<-c()
for(i in 1:20){
y[i]<-rgamma(1,shape=4,scale=(exp(eta[i])/4))
}
df<-data.frame(x=x,y=y)

#glm fitting
model<-glm(y~x,family=Gamma(link='log'),data=df)
summary(model) # result

#Dispersion parameter estimation
#1.Deviance method
(model$deviance/model$df.residual)  
#2.Pearson method
sum(resid(model,type='pear')^2)/model$df.residual # estimation used in summary(model)

The result of running this code is below:

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

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.21011  -0.36473   0.06017   0.33754   0.42223  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.1796     0.1845   6.392  5.1e-06 ***
x            -0.5168     0.0154 -33.548  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 0.1577925)

    Null deviance: 115.4624  on 19  degrees of freedom
Residual deviance:   3.8942  on 18  degrees of freedom
AIC: -143.07

Number of Fisher Scoring iterations: 5

> 
> #Dispersion parameter estimation
> #1.Deviance method
> (model$deviance/model$df.residual)  
[1] 0.2163464
> #2.Pearson method
> #resid(model,type='pear')
> sum(resid(model,type='pear')^2)/model$df.residual # estimation used in summary(model)
[1] 0.1577925
> 

And I guess the discrepancy shows they are 'crude' estimates.