Why Estimate Dispersion From Residuals?
In a quasi-likelihood model, the mean structure and variance function specify how the response varies around its mean, but the dispersion parameter \(\phi\) sets the overall scale of that variation. It is usually unknown. A residual-based estimate is useful for describing the observed variation relative to the model’s variance function and, where appropriate, for scaling model-based standard errors.
Two common estimates use the Pearson statistic and the deviance. Both compare observed responses with fitted means. Neither is an exact finite-sample formula for \(\phi\) in a general regression model: fitting the mean parameters reduces residual variation, and the familiar degrees-of-freedom correction compensates for that reduction approximately.
Why the Denominator Is Usually \(n-p\)
If the means were known rather than estimated, then each standardized squared residual would have expectation \(\phi\): $$ E\left[\frac{(Y_i-\mu_i)^2}{V(\mu_i)}\right] =\frac{\operatorname{Var}(Y_i)}{V(\mu_i)} =\phi. $$ The sum would therefore have expectation \(n\phi\). But fitted means are chosen to make the residuals smaller. In a regular regression, the reduction is approximately equivalent to using up \(p\) residual degrees of freedom.
The next result gives the first-order justification. It is an approximation from local linearization, not an exact identity for arbitrary finite-sample generalized linear models.
To see the calculation, let \(z_i=(Y_i-\mu_i)/\sqrt{V(\mu_i)}\), so the vector \(z\) has mean zero and covariance matrix \(\phi I_n\), where \(I_n\) is the \(n\)-dimensional identity matrix. Write \(\eta_i\) for the linear predictor and let \(x_i\) be the row of covariates. A first-order change in the fitted mean is \((d\mu_i/d\eta_i)x_i^\mathsf{T}(\hat\beta-\beta)\). After standardizing by \(\sqrt{V(\mu_i)}\), the local design matrix \(A\) has rows \((d\mu_i/d\eta_i)x_i^\mathsf{T}/\sqrt{V(\mu_i)}\). The linearized likelihood or quasi-likelihood score equation is \(A^\mathsf{T}\{z-A(\hat\beta-\beta)\}=0\), which gives \(\hat\beta-\beta\approx(A^\mathsf{T}A)^{-1}A^\mathsf{T}z\).
The first-order fitted standardized residual vector is \((I_n-H)z\), where \(H=A(A^\mathsf{T}A)^{-1}A^\mathsf{T}\) is the projection matrix onto the columns of \(A\). Since \(H\) is symmetric, idempotent, and has rank \(p\), $$ E\left[z^\mathsf{T}(I_n-H)z\right] =\operatorname{tr}\left((I_n-H)E[zz^\mathsf{T}]\right) =\phi\operatorname{tr}(I_n-H) =\phi(n-p). $$ This is the first-order expected sum of squared standardized residuals. Replacing the unknown means and parameter values by their fitted estimates gives the usual Pearson statistic, so division by \(n-p\) yields the stated moment estimate. The approximation relies on the local linearization; curvature, high leverage, or a poor mean model can make it inaccurate. \(\square\)
Worked Example: Pearson Dispersion for Poisson Counts
Suppose four independent Poisson counts are modeled with a common mean \(\mu\). The observed counts are \(1,2,3,4\). The model has one fitted parameter, the common mean, so \(n=4\), \(p=1\), and the residual degrees of freedom are \(4-1=3\).
The fitted mean is the sample average: $$ \hat\mu=\frac{1+2+3+4}{4}=2.5. $$ For a Poisson response, \(V(\mu)=\mu\). Thus $$ X^2=\frac{(1-2.5)^2}{2.5}+\frac{(2-2.5)^2}{2.5} +\frac{(3-2.5)^2}{2.5}+\frac{(4-2.5)^2}{2.5} =\frac{2.25+0.25+0.25+2.25}{2.5}=2. $$ The Pearson estimate is \(\hat\phi_P=2/3\approx0.667\). It is below one, but with only three residual degrees of freedom it is highly variable and should not be read as precise evidence for underdispersion.
For comparison, the Poisson unit deviance is \(d(y,\mu)=2\{y\log(y/\mu)-(y-\mu)\}\), with \(0\log(0/\mu)=0\). Since the fitted mean equals the observed average, the terms involving \(y-\hat\mu\) sum to zero: $$ D=2\sum_{i=1}^4 y_i\log(y_i/2.5) =2\{1\log(0.4)+2\log(0.8)+3\log(1.2)+4\log(1.6)\} \approx 2.1288. $$ The deviance estimate is \(D/3\approx0.7096\). The two estimates differ, as they need not coincide in a small sample.
Deviance-Based Dispersion Estimates
The deviance is another measure of residual discrepancy. In a unit-dispersion exponential-family model, the unit deviance for a response \(y\) and mean \(\mu\) is twice the log-likelihood difference between the saturated model for that observation and the model with mean \(\mu\). A saturated model allows the fitted mean for the observation to equal its observed value.
The qualification about unit dispersion matters. In an exponential-dispersion model with scale \(\phi\), the log-likelihood difference is scaled by \(1/\phi\): twice the log-likelihood difference is the unit deviance divided by \(\phi\). For example, in the Gaussian model with variance \(\phi\), the unit deviance is \((y-\mu)^2\), whereas twice the log-likelihood difference is \((y-\mu)^2/\phi\).
For a specified variance function \(V\), the quasi-deviance contribution can be defined by $$ d(y,\mu)=2\int_\mu^y \frac{y-t}{V(t)}\,dt, $$ when the integral is well-defined over the interval between \(\mu\) and \(y\). Summing these contributions at the fitted means gives \(D=\sum_i d(Y_i,\hat\mu_i)\). For a compatible exponential-family model this agrees with the usual unit deviance. The common deviance moment estimate is $$ \hat\phi_D=\frac{D}{n-p}, \qquad n>p. $$
For the derivation, put \(\delta=y-\mu\) and change variables to \(t=\mu+s\): $$ d(y,\mu)=2\int_0^\delta \frac{\delta-s}{V(\mu+s)}\,ds. $$ Continuity and positivity of \(V\) imply \(1/V(\mu+s)=1/V(\mu)+o(1)\) as \(s\) tends to zero. The leading term is therefore $$ \frac{2}{V(\mu)}\int_0^\delta(\delta-s)\,ds =\frac{2}{V(\mu)}\cdot\frac{\delta^2}{2} =\frac{\delta^2}{V(\mu)}. $$ The remaining contribution is \(o(\delta^2)\), since the difference \(1/V(\mu+s)-1/V(\mu)\) tends uniformly to zero on the shrinking interval between \(0\) and \(\delta\). This proves the expansion. Applied observation by observation, it explains why the two residual statistics often track one another when deviations from fitted means are modest. \(\square\)
Worked Example: Gaussian Regression and the Dispersion Scale
Consider four independent Gaussian observations with a common fitted mean and variance \(\phi\). Let the observed values be \(1,2,2,5\). The fitted mean is $$ \hat\mu=\frac{1+2+2+5}{4}=2.5, $$ and there is one fitted parameter, so \(n-p=4-1=3\).
For the Gaussian variance function \(V(\mu)=1\), the Pearson statistic is the residual sum of squares: $$ X^2=(1-2.5)^2+(2-2.5)^2+(2-2.5)^2+(5-2.5)^2 =2.25+0.25+0.25+6.25=9. $$ Thus \(\hat\phi_P=9/3=3\). The Gaussian unit deviance is also the squared residual, so \(D=9\) and \(\hat\phi_D=3\).
This example also checks the scale convention. At a single observation, the unit deviance is \((y-\mu)^2\), but with Gaussian variance \(\phi\), twice the log-likelihood difference from the saturated fit is \((y-\mu)^2/\phi\). They agree only when \(\phi=1\); the deviance itself is not the scaled log-likelihood difference.
Worked Example: Grouped Binomial Dispersion Estimates
Four independent groups each contain 10 trials. The observed success counts are \(2,4,6,8\), and the model has a common success probability. For a count \(Y_i\) from 10 trials, the variance function at success probability \(\mu_i/10\) gives variance \(10p(1-p)\) when \(\phi=1\). There is one fitted probability parameter, and the fitted probability is $$ \hat p=\frac{2+4+6+8}{4\cdot10}=\frac{20}{40}=0.5. $$
Each fitted expected count is 5, and the unit-dispersion variance of each count is \(10(0.5)(0.5)=2.5\). The Pearson statistic and estimate are $$ X^2=\frac{(2-5)^2+(4-5)^2+(6-5)^2+(8-5)^2}{2.5} =\frac{9+1+1+9}{2.5}=8, \qquad \hat\phi_P=\frac{8}{4-1}=\frac{8}{3}\approx2.667. $$
For a grouped binomial count \(y\) out of 10, the unit deviance contribution is $$ d(y,10\hat p)=2\left\{y\log\left(\frac{y}{5}\right)+(10-y)\log\left(\frac{10-y}{5}\right)\right\}. $$ The contributions for \(y=2\) and \(y=8\) are equal: \(2\{2\log(0.4)+8\log(1.6)\}\approx3.8549\) each. Those for \(y=4\) and \(y=6\) are also equal: \(2\{4\log(0.8)+6\log(1.2)\}\approx0.4027\) each. Hence \(D\approx2(3.8549)+2(0.4027)=8.5152\), and $$ \hat\phi_D=\frac{8.5152}{3}\approx2.8384. $$ Both estimates suggest extra variation relative to the binomial variance, but the sample has only three residual degrees of freedom. The numerical difference between the estimates is not surprising.
When the Estimates Can Be Unreliable
The \(n-p\) adjustment is a useful first-order convention, not a guarantee of unbiasedness. A model with many parameters relative to the sample size leaves little residual information; when \(n-p\) is small, either estimate can vary widely. If the local design is rank deficient, counting nominal coefficients overstates the number of independently estimated parameters. The relevant rank is the rank of the local design.
Sparse responses, fitted means near a boundary, and highly influential observations can also make the two statistics unstable or lead them to tell different stories. A single large Pearson residual contributes its squared standardized residual directly to \(X^2\); one unusual count can therefore dominate \(\hat\phi_P\). Deviance contributions use a different scale and can respond differently, particularly when observations are sparse. These estimates also cannot distinguish genuine extra variation from a misspecified mean structure or dependence among responses.
In practice, inspect the residual contributions as well as the aggregate estimates, and consider the Pearson and deviance results together. A large estimate should prompt investigation of the model and data, not an automatic conclusion that the variance scale alone is responsible. Conversely, a small estimate from a low-degree-of-freedom fit is weak evidence that extra variation is absent.
Common Pitfall
Check Your Understanding
Use the definitions and first-order reasoning above to answer the following questions.
- For \(n=40\) independent responses and a full-rank model with \(p=5\) fitted regression parameters, what denominator is used in the usual Pearson dispersion estimate?
- Why does fitting the mean parameters reduce the expected Pearson residual sum from approximately \(n\phi\) to approximately \((n-p)\phi\)?
- In a Gaussian model with variance \(\phi\), distinguish the unit deviance from twice the log-likelihood difference for one observation.
- What local condition makes Pearson and deviance contributions approximately equal?
- Name two features of a fitted data set or model that can make residual-based dispersion estimates unreliable, and explain why either matters.