Pearson Goodness-of-Fit Statistics

What the Pearson Statistic Measures

A fitted GLM gives a mean for each response, but agreement between observed and fitted values must be judged relative to the response’s variability. A difference of two counts may be unremarkable when the expected count is large and conspicuous when it is small. The Pearson statistic measures each difference in units of its fitted standard deviation, squares it, and adds the contributions across observations.

For independent binomial or Poisson responses, the variance is determined by the fitted mean and the model. Thus the aggregate statistic combines the model’s predictions with its variance assumptions. It is an absolute-fit summary: it measures discrepancies between observed responses and the fitted model, rather than comparing two nested models.

Definition: For independent responses \(Y_i\) with fitted means \(\hat{\mu}_i\) and fitted variances \(\widehat{\operatorname{Var}}(Y_i)\), the Pearson contribution for observation \(i\) is \(r_{Pi}^2=(Y_i-\hat{\mu}_i)^2/\widehat{\operatorname{Var}}(Y_i)\). The aggregate Pearson statistic is \(X_P^2=\sum_{i=1}^{n}r_{Pi}^2\). For the binomial and Poisson GLMs considered here, the dispersion is fixed at one.

The subscript \(P\) distinguishes this statistic from other ways of summarizing fit. For a binomial observation \(Y_i\sim\operatorname{Binomial}(m_i,p_i)\), its mean is \(m_ip_i\) and its variance is \(m_ip_i(1-p_i)\). For a Poisson observation \(Y_i\sim\operatorname{Poisson}(\mu_i)\), its mean and variance are both \(\mu_i\). In each case, the fitted variance goes in the denominator of the Pearson contribution.

Definition: For grouped binomial data, with \(y_i\) successes in \(m_i\) trials and fitted probability \(\hat p_i\), \(X_P^2=\sum_i (y_i-m_i\hat p_i)^2/\{m_i\hat p_i(1-\hat p_i)\}\). For Poisson data with observed counts \(y_i\) and fitted means \(\hat\mu_i\), \(X_P^2=\sum_i (y_i-\hat\mu_i)^2/\hat\mu_i\).

Calculating the Statistic in Fitted Models

The fitted probabilities or means must come from the model being assessed. For example, in a GLM with linear predictor \(\eta_i=x_i^T\beta\) and link \(g\), the mean satisfies \(g(\mu_i)=\eta_i\). After fitting, substitute the fitted parameter values to obtain \(\hat{\mu}_i\), and then calculate each contribution using the corresponding binomial or Poisson variance.

Worked Example: An Intercept-Only Grouped Binomial Model

Three groups have the following numbers of successes and trials: \((y_i,m_i)=(4,10),(11,20),(9,15)\). Fit an intercept-only binomial model, so all groups have a common probability \(p\), and calculate \(X_P^2\).

The maximum-likelihood estimate is the total number of successes divided by the total number of trials: \(\hat p=(4+11+9)/(10+20+15)=24/45=8/15\). The fitted success counts are therefore \(10(8/15)=16/3\), \(20(8/15)=32/3\), and \(15(8/15)=8\). The fitted variances are \(m_i\hat p(1-\hat p)\), with \(1-\hat p=7/15\).

For the first group, the residual is \(4-16/3=-4/3\), and the fitted variance is \(10(8/15)(7/15)=112/45\). Its Pearson contribution is \((16/9)/(112/45)=5/7\).

For the second group, the residual is \(11-32/3=1/3\), and the fitted variance is \(20(8/15)(7/15)=224/45\). Its contribution is \((1/9)/(224/45)=5/224\).

For the third group, the residual is \(9-8=1\), and the fitted variance is \(15(8/15)(7/15)=56/15\). Its contribution is \(1/(56/15)=15/56\). Adding the three gives \(X_P^2=5/7+5/224+15/56=225/224\approx1.0045\).

There are three grouped responses and one fitted parameter, the common probability. The usual large-sample reference would therefore have \(3-1=2\) degrees of freedom. This calculation gives the statistic; whether the reference is reliable depends on the sampling regime, discussed below.

Worked Example: An Intercept-Only Poisson Model

Suppose four independent Poisson counts are \(8,10,12,14\), and fit an intercept-only Poisson model. The fitted mean is the sample mean: \(\hat\mu=(8+10+12+14)/4=11\).

The four residuals are \(-3,-1,1,3\). Since the fitted variance for each count is 11, the contributions are \(9/11,1/11,1/11,9/11\), respectively. Thus \(X_P^2=(9+1+1+9)/11=20/11\approx1.8182\).

One parameter, the common Poisson mean, has been estimated from four responses, so the conventional reference uses \(4-1=3\) degrees of freedom. The statistic is not itself a probability or a goodness-of-fit conclusion; it is compared with a reference distribution only when the approximation is appropriate.

Worked Example: A Grouped Binomial Logistic Model

Consider three groups with \(m_i=20\) trials each and covariate values \(x_i=0,1,2\). Fit a logistic model with \(\operatorname{logit}(p_i)=\beta_0+\beta_1x_i\). Suppose the observed successes are \(5,8,17\), and the fitted probabilities are \(0.2,0.5,0.8\).

These fitted probabilities correspond to \(\hat\beta_0=\log(0.2/0.8)=-\log 4\) and \(\hat\beta_1=\log 4\): the fitted logits at \(x=0,1,2\) are \(-\log4,0,\log4\). The fitted success counts are \(4,10,16\), so the residuals are \(1,-2,1\).

The fitted variances are \(20(0.2)(0.8)=3.2\), \(20(0.5)(0.5)=5\), and \(20(0.8)(0.2)=3.2\). The Pearson statistic is \(X_P^2=1^2/3.2+(-2)^2/5+1^2/3.2=0.3125+0.8+0.3125=1.425\).

To check that these are indeed the fitted probabilities, the logistic-regression score equations require the residuals to sum to zero and their covariate-weighted sum to be zero. Here \(1-2+1=0\) and \(0(1)+1(-2)+2(1)=0\). The design has rank two, and the binomial log-likelihood is strictly concave in these coefficients, so these score equations identify its maximum. There are three groups and two fitted coefficients, giving the conventional reference with one degree of freedom. With only three groups, this example is chiefly an arithmetic illustration, not evidence that a reference approximation is dependable.

Why a Chi-Square Reference Can Arise

At the true parameter, each Pearson contribution has expectation one, provided its variance is positive and the model’s mean and variance are correct. This is an exact expectation statement before fitting; it does not say that the sum has a chi-square distribution.

Result: If \(Y_i\) has mean \(\mu_i\) and positive variance \(v_i\), then \(E[(Y_i-\mu_i)^2/v_i]=1\). Consequently, for \(G\) independent responses with known true means and variances, the sum of their Pearson contributions has expectation \(G\).

The proof is the definition of variance: \(E[(Y_i-\mu_i)^2]=\operatorname{Var}(Y_i)=v_i\). Dividing by \(v_i\) gives expectation one, and linearity of expectation gives \(G\) for the sum. When parameters are estimated, fitted residuals are constrained by the model’s score equations. Under suitable large-sample conditions, this fitting removes approximately one independent direction for each estimated coefficient. This leads to the familiar reference with \(G-q\) degrees of freedom, where \(q\) is the rank of the fitted model’s design matrix.

Theorem: Consider a fixed number \(G\) of independent grouped binomial or Poisson responses and a GLM with \(q\) identifiable regression parameters, where \(q<G\). Suppose the true parameter is interior, the model is correctly specified, the mean function is smooth, the design has rank \(q\), and the maximum-likelihood estimator is consistent and regular. Suppose also that information increases so each binomial group has \(m_i\to\infty\), with probabilities bounded away from zero and one, or each Poisson mean tends to infinity. Under the usual nonsingularity conditions and a nondegenerate limiting fitted design, the fitted Pearson statistic satisfies \(X_P^2\xrightarrow{d}\chi^2_{G-q}\).

Here is the projection argument behind the result. Write \(\mu\) for the vector of true means, \(V\) for the diagonal matrix of true variances, and \(D\) for the \(G\)-by-\(q\) Jacobian matrix whose \(i\)th row is the derivative of \(\mu_i\) with respect to the parameter. Define the standardized response error \(Z=V^{-1/2}(Y-\mu)\), and set \(A=V^{-1/2}D\). The increasing-information assumptions give a joint normal limit for \(Z\), with limiting covariance the identity matrix. A first-order likelihood expansion gives the fitted parameter displacement \(\hat\beta-\beta\approx(D^TV^{-1}D)^{-1}D^TV^{-1}(Y-\mu)\).

Multiplying the fitted mean displacement by \(V^{-1/2}\) shows that the standardized fitted residual vector is, to first order, \(Z-A(A^TA)^{-1}A^TZ=(I-P)Z\), where \(P=A(A^TA)^{-1}A^T\). The matrix \(P\) is the orthogonal projection onto the \(q\)-dimensional column space of \(A\); hence \(I-P\) has rank \(G-q\). Replacing true variances by fitted variances does not change this limiting quadratic form under the stated regularity conditions. Therefore \(X_P^2\) converges in distribution to \(Z^T(I-P)Z\).

To identify that limit, choose an orthonormal basis in which \(I-P\) has \(G-q\) diagonal entries equal to one and all remaining entries equal to zero. A standard normal vector remains standard normal after an orthogonal change of coordinates. Its quadratic form in that basis is the sum of squares of \(G-q\) independent standard normal variables, which has a \(\chi^2_{G-q}\) distribution. This proves the stated reference limit. \(\square\)

Conditions and Interpretation

The theorem describes a particular asymptotic regime: a fixed collection of groups becomes increasingly informative. It does not establish a chi-square approximation merely because a dataset has many rows. In particular, a large number of individual binary observations does not automatically behave like a fixed set of binomial groups whose trial counts all grow. The distributional approximation must fit the way information accumulates in the data.

Even when the reference is reasonable, a Pearson statistic is an aggregate. A small value indicates that the total standardized squared discrepancy is small relative to the reference, while a large value indicates more discrepancy than expected under the model. The sum does not identify which observations drive it, establish why discrepancies occur, or prove that the model is correct. Inspecting individual contributions may help locate discrepancies, but the statistic alone is not a complete diagnostic analysis.

Common Mistake

Common mistake: “The Pearson statistic is approximately chi-square because it is a sum of squared residuals, so any dataset with many observations can be compared with \(\chi^2_{n-q}\).”

Correction: Squaring and adding residuals does not by itself produce a chi-square distribution. The reference follows from an appropriate normal approximation and the effect of estimating parameters. In the fixed-group result above, the number of groups is \(G\), the fitted model has rank \(q\), and information increases within groups. The relevant reference is then approximately \(\chi^2_{G-q}\), not automatically \(\chi^2_{n-q}\) for any representation of the data.

Key takeaway: Calculate each Pearson contribution using the model’s fitted variance, then sum. For grouped binomial or Poisson GLMs, a \(\chi^2_{G-q}\) reference is justified in a suitable large-information regime; it is an approximation, not a universal property of the statistic.

Check Your Understanding

Use the definitions and derivations above to answer the following questions.

For a binomial group with \(m_i\) trials and fitted probability \(\hat p_i\), write its Pearson contribution and explain why its denominator includes both \(\hat p_i\) and \(1-\hat p_i\).
For Poisson observations \(y_i\) with fitted means \(\hat\mu_i\), what is the aggregate Pearson statistic?
In a fixed set of \(G\) groups with a fitted model of rank \(q\), what degrees of freedom appear in the usual large-sample reference, and why does fitting affect them?
Why does the exact expectation of one for each true-parameter Pearson contribution not prove that the aggregate statistic is chi-square?
What feature of the sampling regime is required by the fixed-group asymptotic result?