Pearson-Based Tests for Overdispersion

Why a Pearson Test Needs More Than a Variance Estimate

A Pearson statistic compares observed responses with their fitted means, after scaling each discrepancy by its model-based variance. A large value can signal that the observed variation exceeds the equidispersion level. But its reference distribution is not determined just by the fact that the variance function has been specified: it also depends on how the means were fitted and on the distributional behavior of the standardized residuals.

The key degrees-of-freedom question is therefore precise: when does fitting $q$ mean parameters remove exactly $q$ directions of variation from the Pearson statistic? The answer is not “whenever $q$ parameters were estimated.” The familiar $m-q$ calibration follows when the fitted mean changes, to first order, by the orthogonal projection of the standardized data onto a $q$-dimensional model space.

The Statistic and Its Null Hypothesis

Suppose $Y_1,\ldots,Y_m$ are independent responses with means $\mu_i$, and suppose their variances are modeled as $\operatorname{Var}(Y_i)=\phi V(\mu_i)$, where $V(\mu_i)>0$ is known up to the mean parameters and $\phi$ is a dispersion multiplier. Equidispersion means $\phi=1$; overdispersion corresponds to $\phi>1$. Let $\hat{\mu}_i$ be fitted means under the specified mean model.

Definition: The Pearson goodness-of-fit statistic for the fitted mean model is $$ Q=\sum_{i=1}^m \frac{(Y_i-\hat{\mu}_i)^2}{V(\hat{\mu}_i)}. $$ A Pearson-based upper-tail test of equidispersion compares $Q$ with a chi-square reference distribution. Its degrees of freedom depend on the residual projection induced by the fitting procedure, not just on the number of estimated parameters.

The statistic is dimensionless: each squared discrepancy is divided by its model-based variance at dispersion $\phi=1$. If the means were known and standardized errors were approximately independent standard normal variables, each summand would be approximately a squared standard normal. Estimating the means changes that argument because the fitted residuals are no longer the original standardized errors.

How Fitting the Mean Changes the Residuals

Write $\boldsymbol{\mu}=(\mu_1,\ldots,\mu_m)^\mathsf{T}$, and define the standardized error vector under the equidispersion null by $$ z_i=\frac{Y_i-\mu_i}{\sqrt{V(\mu_i)}},\qquad i=1,\ldots,m. $$ Let $\beta$ denote the $q$-dimensional mean parameter, with true value $\beta_0$. Let $D$ be the $m\times q$ derivative matrix whose $i$th row is $\partial\mu_i/\partial\beta^\mathsf{T}$, evaluated at $\beta_0$, and let $W=\operatorname{diag}\{1/V(\mu_i)\}$. The standardized local changes in the mean lie in the columns of $$ B=W^{1/2}D. $$ Assume $B$ has full column rank $q$. The matrix projecting orthogonally onto this column space is $$ H=B(B^\mathsf{T}B)^{-1}B^\mathsf{T}. $$ In particular, the assertion that fitting removes the component $Hz$ requires an appropriate fitting rule. Weighted least squares has this property directly in its linearized problem; a regular quasi-score fit has the same first-order projection under the usual local linearization conditions. A generic consistent mean estimator need not have it.

Theorem: Suppose $m$ and $q$ are fixed, the mean model is correctly specified, $B$ has rank $q$, and under $\phi=1$ the standardized error vector satisfies $z\xrightarrow{d}N_m(0,I_m)$. Assume the fitted standardized residual vector $r$, with components $r_i=(Y_i-\hat{\mu}_i)/\sqrt{V(\hat{\mu}_i)}$, satisfies $$ r=(I-H)z+o_p(1). $$ Then $Q=r^\mathsf{T}r\xrightarrow{d}\chi^2_{m-q}$. Under a corresponding dispersion sequence with $\phi>0$, if the standardized errors converge in distribution to $N_m(0,\phi I_m)$ and the same projection applies, then $Q\xrightarrow{d}\phi\chi^2_{m-q}$.

The assumptions describe the calibration, not merely the construction of $Q$. In particular, the result requires both an approximately normal standardized error vector and a first-order orthogonal-projection representation of the fitted residuals. A quasi-likelihood variance specification by itself need not imply the first condition.

Proof of the Chi-Square Calibration

First consider the linearized fitting step. If the standardized mean change is $B\delta$, weighted least squares chooses $\delta$ to minimize $\|z-B\delta\|^2$. Differentiating this quadratic objective gives the normal equations $$ B^\mathsf{T}B\hat{\delta}=B^\mathsf{T}z. $$ Because $B$ has full column rank, $B^\mathsf{T}B$ is invertible, and hence $\hat{\delta}=(B^\mathsf{T}B)^{-1}B^\mathsf{T}z$. The fitted standardized change is $B\hat{\delta}=Hz$, so the linearized residual is $(I-H)z$. This derives the projection for that fitting rule; it does not establish it for every possible estimator.

By the theorem's assumed residual expansion, $r=(I-H)z+o_p(1)$. Since $m$ is fixed and $z\xrightarrow{d}N_m(0,I_m)$, $z=O_p(1)$. Consequently, the difference between the squared norms of $r$ and $(I-H)z$ is $o_p(1)$: expanding the square gives a cross-product of an $O_p(1)$ vector with an $o_p(1)$ vector, plus an $o_p(1)$ squared norm. Thus $$ Q=z^\mathsf{T}(I-H)z+o_p(1), $$ where $I-H$ is symmetric and idempotent.

To identify the limit, choose an orthogonal matrix $U$ that diagonalizes $I-H$. Since $H$ is an orthogonal projection of rank $q$, $I-H$ has $m-q$ eigenvalues equal to 1 and $q$ equal to 0. Also $U^\mathsf{T}z\xrightarrow{d}N_m(0,I_m)$, because an orthogonal transformation preserves the standard multivariate normal distribution. Therefore the quadratic form $z^\mathsf{T}(I-H)z$ converges in distribution to the sum of squares of $m-q$ independent standard normal variables, which is $\chi^2_{m-q}$. Slutsky's theorem gives the asserted limit for $Q$.

If, under the corresponding dispersion sequence, the standardized errors converge in distribution to $N_m(0,\phi I_m)$, the same argument writes the limiting standardized error vector as $\sqrt{\phi}\,Z$, where $Z\sim N_m(0,I_m)$. The quadratic form is then $\phi Z^\mathsf{T}(I-H)Z$, giving the scaled chi-square limit. $\square$

The same projection argument explains the degrees of freedom: $H$ has rank $q$, so $I-H$ retains $m-q$ orthogonal directions. If the design or derivative matrix has rank $r_B<q$, the projection has rank $r_B$, and the corresponding residual degrees of freedom are $m-r_B$, provided the other assumptions still hold.

Worked Calculations

Worked Example: Known Poisson Means

Suppose four independent Poisson counts have specified means $100,100,100,100$, and the observed counts are $110,90,105,95$. The means are known rather than fitted, so no mean parameters are estimated. For Poisson responses, $V(\mu_i)=\mu_i$, and the Pearson statistic is $$ Q=\frac{(110-100)^2}{100}+\frac{(90-100)^2}{100}+\frac{(105-100)^2}{100}+\frac{(95-100)^2}{100} =1+1+0.25+0.25=2.5. $$

Here $q=0$, so the asymptotic reference degrees of freedom are $m-q=4$. The upper-tail probability is $P(\chi^2_4\geq2.5)=e^{-1.25}(1+1.25)\approx0.6446$. This is not evidence of unusually large Pearson variation. The chi-square calibration is approximate for these Poisson counts; it becomes more credible when standardized counts are close to normal.

Worked Example: Estimating a Common Poisson Mean

Now suppose five independent Poisson observations share one unknown mean and have values $80,120,100,140,60$. The maximum-likelihood estimate is their average: $$ \hat{\mu}=\frac{80+120+100+140+60}{5}=100. $$ There is one fitted parameter, so the projection argument gives $m-q=5-1=4$ residual degrees of freedom.

The Pearson statistic is $$ Q=\frac{(80-100)^2}{100}+\frac{(120-100)^2}{100}+\frac{(100-100)^2}{100}+\frac{(140-100)^2}{100}+\frac{(60-100)^2}{100} =4+4+0+16+16=40. $$

Using $\chi^2_4$, the upper-tail probability is $e^{-40/2}(1+40/2)=21e^{-20}\approx4.33\times10^{-8}$. Under the stated approximation, this is strong evidence against equidispersion. The subtraction of one degree of freedom is justified here by the common-mean fit's projection onto the one-dimensional intercept space; it is not a universal consequence of having estimated one parameter.

Worked Example: Checking the Degrees of Freedom and Tail Probability

Consider ten independent normal responses with known unit variance and a two-parameter linear mean, fitted by least squares. For a concrete data set, take $x_i=i$ for $i=1,\ldots,10$, and define $$ Y_i=10+2x_i+r_i,\qquad (r_1,\ldots,r_{10})=\sqrt{6}(1,-1,0,0,0,0,0,0,-1,1). $$ The residual vector is orthogonal to both the intercept and the $x$ column: its entries sum to zero, and its inner product with $x$ is $\sqrt{6}(1-2-9+10)=0$. Thus the least-squares fitted line is $10+2x_i$, with exactly these residuals.

The Pearson statistic, which here is the residual sum of squares because every variance is one, is $$ Q=\sum_{i=1}^{10}r_i^2=4(6)=24. $$ There are $m=10$ responses and $q=2$ fitted mean parameters, so the residual degrees of freedom are $8$. For this normal model the $\chi^2_8$ calibration is exact. Its upper-tail probability is $$ P(\chi^2_8\geq24)=e^{-12}\left(1+12+\frac{12^2}{2}+\frac{12^3}{6}\right)=373e^{-12}\approx0.00229179. $$ Thus the probability is about $0.00229$, or $0.0023$ to two significant figures, under the equidispersion model.

When the Degrees-of-Freedom Shortcut Fails

It is tempting to say that $m$ squared residuals contribute $m$ degrees of freedom and each fitted mean parameter removes one. That reasoning skips the central issue: which directions does the fitting procedure remove? The $m-q$ result follows when the residual operation is the orthogonal projection $I-H$ in standardized coordinates and the standardized errors are asymptotically spherical normal. A nonlinear mean model can still have this local projection if its estimating equations produce it, but mere consistency of $\hat{\beta}$ is not enough.

For example, a different estimator can have a different first-order influence on the fitted means. Its residual vector may have a covariance other than $I-H$; a quadratic form in such residuals generally has a weighted sum of chi-square variables rather than a single chi-square distribution with $m-q$ degrees of freedom. Dependence among observations can also invalidate the spherical covariance assumption. In these cases, report that the usual Pearson chi-square calibration is not justified without deriving an estimator-specific reference distribution or using another suitable calibration method.

A second pitfall is interpreting a large $Q$ as proof of overdispersion. A wrong mean structure, influential observations, dependence, or an incorrect variance function can also produce a large statistic. The test assesses compatibility with the specified mean and variance assumptions together; it does not identify the source of a discrepancy. Nor does a quasi-likelihood variance specification alone supply a full finite-sample distribution for $Q$.

Key takeaway: Under a valid standardized-residual projection and an appropriate normal approximation, $Q$ has an asymptotic $\chi^2_{m-q}$ distribution under equidispersion. Subtract the number of fitted mean parameters only when the fitting rule actually removes those orthogonal standardized directions.

Check Your Understanding

Use the assumptions of the projection result when deciding whether a chi-square calibration applies.

If there are $m=18$ responses and a full-rank projection onto a mean space of dimension $q=4$, what are the residual degrees of freedom?
Why does estimating $q$ parameters by itself not prove that the Pearson statistic has $m-q$ chi-square degrees of freedom?
What additional distributional condition on standardized errors is used in the asymptotic calibration?
For the ten-response example, verify that the displayed residual vector is orthogonal to both the intercept and slope columns.
Give two reasons, other than overdispersion, that a Pearson statistic might be unusually large.