Gamma Dispersion Estimation By Deviance And Likelihood

Why Estimate Gamma Dispersion?

A gamma response model specifies how the response mean changes with predictors, but it also needs to describe variation around that mean. With mean \(\mu_i\) and dispersion \(\phi>0\), the gamma model has variance \(\operatorname{Var}(Y_i)=\phi\mu_i^2\). Dispersion therefore controls the relative, rather than absolute, size of variation: the standard deviation is \(\sqrt{\phi}\mu_i\). Estimating \(\phi\) tells us how much scatter remains around the fitted means.

Assume throughout that \(Y_1,\ldots,Y_n\) are independent gamma responses with means \(\mu_i>0\) and a common dispersion \(\phi\). Write \(\hat\mu_i\) for the fitted mean. The estimates below differ in how they summarize residual variation: a squared-residual moment, the Pearson statistic, the deviance, or the full gamma likelihood.

Method of Moments and the Pearson Estimate

If the means were known, the variance relationship immediately gives a moment estimate. Define \(R_i=(Y_i-\mu_i)/\mu_i\). Its expected value is zero and its second moment is

Matching the average of the observed squared relative residuals to this expected value gives the method-of-moments estimate. With unknown means, substitute the fitted means for the unknown \(\mu_i\). This yields the raw residual moment below. A degrees-of-freedom adjustment is often used when the means have been estimated.

The approximation behind \(n-p\) can be seen by linearizing the fitted residuals near the true means. In standardized coordinates, suppose the residual vector is approximated by \((I-H)e\), where \(e\) has mean zero and covariance matrix \(\phi I\), and \(H\) is a symmetric idempotent projection matrix of rank \(p\). Then

This calculation motivates the Pearson correction. In a general gamma GLM, the fitted means and residuals are nonlinear functions of the data, so the projection argument is a local approximation rather than an exact finite-sample result.

The Deviance Estimate

The gamma unit deviance measures the difference between the fitted model and a saturated model, which assigns each observation its own mean \(y_i\). For a common dispersion \(\phi\), the gamma log likelihood can be written, up to terms independent of the fitted mean, with mean-dependent contribution \(-\phi^{-1}[\log\mu_i+y_i/\mu_i]\). Comparing the fitted mean \(\hat\mu_i\) with the saturated mean \(y_i\) gives

The deviance estimate is \(\hat\phi_D=D/(n-p)\). Its degrees-of-freedom denominator is motivated by the fact that, near a good fit, deviance and Pearson variation agree to first order. To see the relationship, let \(z=y/\hat\mu\) and write \(z=1+u\). For small \(u\), Taylor expansion gives

The leading term \(u^2\) is the squared Pearson residual. Consequently, when fitted ratios are close to one, deviance and Pearson estimates tend to be similar. For large residuals, the higher-order terms matter and the estimates need not agree.

Maximum Likelihood for Dispersion

The gamma likelihood uses more than a single residual sum of squares. Let \(k=1/\phi\) be the gamma shape, and suppose the fitted means are fixed at \(\hat\mu_i\). Define \(z_i=y_i/\hat\mu_i\). For positive observations, the part of the log likelihood that depends on \(k\) is

Here \(\Gamma\) is the gamma function and \(\psi(k)=\frac{d}{dk}\log\Gamma(k)\) is the digamma function. Differentiating the displayed log likelihood gives the shape score. Setting it to zero and rearranging gives the likelihood equation for dispersion.

To verify the equation, differentiate:

The score is zero precisely when the stated equation holds. If the fitted means are estimated by maximum likelihood through a model \(\mu_i(\beta)\), the mean-dependent part of the score for \(\beta\) is multiplied by \(k\); its maximizing value of \(\beta\) does not depend on \(k\). Thus profiling over the mean parameters leads to the same dispersion equation, provided the fitted means are positive and the mean-model maximum exists.

There is at most one positive solution when \(D>0\). Indeed, the function \(g(k)=\log k-\psi(k)\) decreases strictly from infinity to zero as \(k\) goes from zero to infinity. For strict decrease, its derivative is \(1/k-\psi_1(k)<0\), since \(\psi_1(k)=\sum_{j=0}^\infty(k+j)^{-2}>\int_0^\infty(k+x)^{-2}\,dx=1/k\). The endpoint limits follow from the standard digamma limits \(\psi(k)\sim-1/k\) as \(k\) decreases to zero and \(\psi(k)=\log k-1/(2k)+O(k^{-2})\) as \(k\) increases. Also, \(z-1-\log z\geq0\) for \(z>0\), with equality only at \(z=1\). Hence \(D>0\) gives a unique root, while \(D=0\) places the likelihood maximum at the boundary \(\phi\downarrow0\). \(\square\)

When the Estimates Are Close—and When They Are Not

The estimates share a common purpose but not an identical calculation. The raw moment estimate matches the average squared relative residual to the gamma variance relationship. The Pearson estimate uses that same squared-residual statistic with a degrees-of-freedom adjustment. The deviance estimate uses the likelihood difference from a saturated fit and applies the same conventional adjustment. Maximum likelihood instead solves a nonlinear equation involving the full gamma likelihood.

Common Pitfall and Practical Interpretation

A frequent mistake is to call \(X_P^2/(n-p)\) an exactly unbiased gamma dispersion estimate in every fitted model. The exact moment identity applies when the true means are known: \(E[(Y_i-\mu_i)^2/\mu_i^2]=\phi\). Once means are fitted from the same responses, residuals are dependent on those estimates, and the simple projection argument is generally only approximate. The correction is to describe the Pearson and deviance estimates as conventional residual-degrees-of-freedom estimates, and to state the model and approximation when interpreting them.

A second mistake is to substitute \(D/(n-p)\) into the likelihood equation. The deviance estimate is a useful residual-based estimator; the likelihood estimate instead solves \(\log k-\psi(k)=D/(2n)\). The denominator \(n\) in this score equation is part of the likelihood derivation, not a residual degrees-of-freedom choice.