Scaling GLM Covariance Matrices

Why Scale a Covariance Matrix?

A quasi-binomial or quasi-Poisson model allows the response variance to be larger or smaller than the variance specified by the corresponding ordinary GLM. That change matters for uncertainty: if the variance is multiplied by a common dispersion parameter, the coefficient covariance estimate is multiplied by that parameter as well. The fitted coefficients, however, need not change. The reason is that the common factor appears in every component of the quasi-score and cancels from its estimating equation.

This tutorial derives that relationship and distinguishes the covariance of the score from the covariance of the coefficient estimator. That distinction is important: the score covariance depends on how the score is scaled, while the resulting estimator covariance is the same under either consistent scaling.

The Model-Based Information Matrix

Consider independent responses \(Y_1,\ldots,Y_n\), with means \(\mu_i(\beta)\) determined by a coefficient vector \(\beta\) of dimension \(p\). Suppose their variances have the quasi-GLM form \(\operatorname{Var}(Y_i)=\phi V(\mu_i)\), where \(\phi>0\) is a common dispersion parameter and \(V\) is the variance function. Write \(D\) for the \(n\) by \(p\) matrix whose \(i\)th row is the derivative of \(\mu_i(\beta)\) with respect to \(\beta^\mathsf{T}\), and write \(W=\operatorname{diag}\{1/V(\mu_i)\}\).

Definition: The unscaled model-based information matrix is

H(\beta)=D^\mathsf{T}WD =\sum_{i=1}^n \frac{1}{V(\mu_i)} \frac{\partial\mu_i}{\partial\beta} \frac{\partial\mu_i}{\partial\beta^\mathsf{T}}.

When \(H(\beta)\) is nonsingular, the ordinary model-based covariance estimate at a fitted value \(\hat\beta\) is \(H(\hat\beta)^{-1}\). Under a common dispersion \(\phi\), its quasi-model counterpart is \(\phi H(\hat\beta)^{-1}\).

For a GLM with linear predictor \(\eta_i=x_i^\mathsf{T}\beta\) and link function \(g(\mu_i)=\eta_i\), the derivative row is \((d\mu_i/d\eta_i)x_i^\mathsf{T}\). Thus \(H\) has the familiar weighted-design form, with weight \((d\mu_i/d\eta_i)^2/V(\mu_i)\). The definition above also applies when the mean is expressed directly as a differentiable function of \(\beta\).

Worked Example: Scaling a Given Information Matrix

Suppose the unscaled information matrix at a fitted model is \(H=\begin{pmatrix}4&2\\2&3\end{pmatrix}\), and the estimated dispersion is \(\hat\phi=2.5\). Find the ordinary and scaled covariance estimates, then calculate the scaled standard error for the second coefficient.

The determinant of \(H\) is \(4(3)-2(2)=8\). Therefore,

H^{-1}=\frac{1}{8} \begin{pmatrix}3&-2\\-2&4\end{pmatrix} = \begin{pmatrix}3/8&-1/4\\-1/4&1/2\end{pmatrix}.

The quasi-model covariance estimate is

\hat\phi H^{-1} = \begin{pmatrix}15/16&-5/8\\-5/8&5/4\end{pmatrix}.

The variance estimate for the second coefficient is its diagonal entry, \(5/4=1.25\). Its standard error is \(\sqrt{5/4}\approx1.118\). The corresponding ordinary standard error would be \(\sqrt{1/2}\approx0.707\). The covariance matrix, including its off-diagonal entries, has been multiplied by \(2.5\); standard errors, being square roots of variances, have been multiplied by \(\sqrt{2.5}\), not by \(2.5\).

Deriving the Dispersion Scaling

Use the quasi-score with the common factor \(1/\phi\) included:

U_\phi(\beta)=\frac{1}{\phi}D^\mathsf{T}W(Y-\mu),

where \(Y=(Y_1,\ldots,Y_n)^\mathsf{T}\) and \(\mu=(\mu_1,\ldots,\mu_n)^\mathsf{T}\). For the next derivation, suppose the responses are independent, have the stated means and variances, and the mean functions are differentiable. Assume \(V(\mu_i)>0\) and \(H(\beta)\) is nonsingular at the parameter value being considered.

Theorem: At the true parameter value, the sensitivity matrix and covariance matrix of \(U_\phi\) are both \(H/\phi\). Consequently, the first-order covariance matrix for the quasi-score estimator is \(\phi H^{-1}\). This is an asymptotic, model-based covariance result; it is not in general an exact finite-sample covariance for the estimator.

Proof. The sensitivity matrix is the negative expected derivative of the estimating function with respect to \(\beta^\mathsf{T}\). Differentiating \(U_\phi\) produces a term from differentiating \(D^\mathsf{T}W\), multiplied by \(Y-\mu\), and a term from differentiating \(-\mu\). The expectation of the first term is zero because \(E[Y-\mu]=0\). The second term gives

A_\phi =-E\left[\frac{\partial U_\phi}{\partial\beta^\mathsf{T}}\right] =\frac{1}{\phi}D^\mathsf{T}WD =\frac{H}{\phi}.

For the score covariance, independence makes the covariance matrix of \(Y-\mu\) diagonal, with entries \(\phi V(\mu_i)\). Hence

\begin{aligned} B_\phi &=\operatorname{Cov}\{U_\phi(\beta)\}\\ &=\frac{1}{\phi^2}D^\mathsf{T}W \operatorname{diag}\{\phi V(\mu_i)\}WD\\ &=\frac{1}{\phi}D^\mathsf{T}WD =\frac{H}{\phi}, \end{aligned}

because \(W\operatorname{diag}\{\phi V(\mu_i)\}W=\phi W\). The estimating-equation sandwich covariance is \(A_\phi^{-1}B_\phi A_\phi^{-\mathsf{T}}\). Substituting the two matrices just derived gives

A_\phi^{-1}B_\phi A_\phi^{-\mathsf{T}} =(\phi H^{-1})(H/\phi)(\phi H^{-1}) =\phi H^{-1}.

The sandwich expression is the first-order covariance for the estimator under the usual estimating-equation regularity conditions: the estimator is consistent, the sensitivity is nonsingular, the estimating functions are differentiable, and a central limit theorem applies to their sum. In a sequence of models, these conditions require suitable control of the contributions from individual observations as sample size grows. In practice, the result motivates substituting fitted means and an estimated dispersion into the covariance formula. \(\square\)

There is an equally valid alternative score scaling. Define \(U_1=\phi U_\phi=D^\mathsf{T}W(Y-\mu)\), which removes the common factor. Its sensitivity is \(H\), and its covariance is \(\phi H\). Its sandwich covariance is therefore \(H^{-1}(\phi H)H^{-1}=\phi H^{-1}\). The score covariance is not the same under the two definitions: \(U_\phi\) has covariance \(H/\phi\), whereas \(U_1\) has covariance \(\phi H\). Both lead to the same coefficient covariance when sensitivity and score covariance are paired consistently.

Worked Example: Quasi-Poisson Covariance From Fitted Means

Consider four Poisson log-link observations with covariates \(x=(0,0,1,1)\) and counts \(Y=(0,4,1,5)\). The fitted means are \(2,2,3,3\), giving \(\hat\beta_0=\log 2\) and \(\hat\beta_1=\log(3/2)\). Verify the fit, estimate dispersion by the Pearson statistic divided by residual degrees of freedom, and calculate the coefficient covariance estimate.

For a Poisson log-link model, \(V(\mu)=\mu\) and \(d\mu/d\eta=\mu\). Thus each observation contributes weight \(\mu_i\) to \(H=X^\mathsf{T}\operatorname{diag}(\mu_i)X\). The residual sums within the two covariate groups are \((0+4)-2(2)=0\) and \((1+5)-2(3)=0\); consequently both components of the score are zero at the stated fit. The information matrix is

H= \begin{pmatrix} 2+2+3+3&3+3\\ 3+3&3+3 \end{pmatrix} = \begin{pmatrix}10&6\\6&6\end{pmatrix}.

Its determinant is \(60-36=24\), so

H^{-1}=\frac{1}{24} \begin{pmatrix}6&-6\\-6&10\end{pmatrix} = \begin{pmatrix}1/4&-1/4\\-1/4&5/12\end{pmatrix}.

The Pearson statistic is the sum of squared residuals divided by fitted means. For the first group the two contributions are \(2^2/2=2\) and \((-2)^2/2=2\), summing to \(4\). For the second they are \((-2)^2/3=4/3\) and \(2^2/3=4/3\), summing to \(8/3\). With \(4-2=2\) residual degrees of freedom, \(\hat\phi=(4+8/3)/2=10/3\). Therefore,

\widehat{\operatorname{Cov}}(\hat\beta) =\frac{10}{3}H^{-1} = \begin{pmatrix}5/6&-5/6\\-5/6&25/18\end{pmatrix}.

The slope standard error is \(\sqrt{25/18}=5/(3\sqrt{2})\approx1.179\). The fitted coefficients are still \(\log 2\) and \(\log(3/2)\): the common dispersion changes the covariance estimate, not the score equation’s solution.

Why the Coefficients Do Not Change

The invariance of the coefficient estimate follows directly from the estimating equation. Suppose \(\phi\) is a positive scalar common to all observations. Then, at any \(\beta\),

U_\phi(\beta)=\frac{1}{\phi}U_1(\beta).

Since \(1/\phi\) is nonzero, \(U_\phi(\beta)=0\) if and only if \(U_1(\beta)=0\). Thus the two equations have the same roots, provided they use the same mean model and the same variance function apart from the common scalar. The conclusion applies when a positive dispersion estimate is substituted after solving for \(\beta\) as well.

This argument does not say that every change to a variance specification leaves coefficient estimates unchanged. If the variance changes in a way that alters the observation-specific relative weights, the estimating equation may have different roots. The invariance here depends on a single common multiplicative factor.

Worked Example: Quasi-Binomial Intercept-Only Fit

Two groups each contain 10 trials, with 2 successes in the first group and 6 in the second. Fit an intercept-only quasi-binomial model with a logit link. Use the Pearson estimate of dispersion and find the scaled standard error for the intercept.

The fitted probability is the total number of successes divided by the total number of trials: \(\hat p=(2+6)/(10+10)=0.4\). Therefore \(\hat\beta=\log\{\hat p/(1-\hat p)\}=\log(2/3)\approx-0.4055\). The information is the sum of the binomial contributions \(m_i\hat p(1-\hat p)\):

H=10(0.4)(0.6)+10(0.4)(0.6)=4.8=\frac{24}{5}, \qquad H^{-1}=\frac{5}{24}.

The expected success count in each group is \(10(0.4)=4\). The Pearson contributions are \((2-4)^2/2.4=5/3\) and \((6-4)^2/2.4=5/3\). There is one residual degree of freedom, so \(\hat\phi=10/3\). The scaled variance and standard error are

\widehat{\operatorname{Var}}(\hat\beta) =\frac{10}{3}\frac{5}{24} =\frac{25}{36}, \qquad \operatorname{SE}(\hat\beta)=\sqrt{\frac{25}{36}}=\frac{5}{6}\approx0.833.

For comparison, the unscaled standard error is \(\sqrt{5/24}\approx0.456\). The fitted probability and intercept remain \(0.4\) and \(\log(2/3)\); only their estimated uncertainty is enlarged by the dispersion factor.

Using the Plug-In Covariance Estimate

In an analysis, both the information matrix and the dispersion are usually estimated. Evaluate \(H\) at the fitted coefficients and use the chosen dispersion estimate, such as a Pearson-based or deviance-based estimate:

\widehat{\operatorname{Cov}}(\hat\beta) =\hat\phi H(\hat\beta)^{-1}.

For any coefficient \(\hat\beta_j\), its estimated standard error is \(\sqrt{\hat\phi}\sqrt{\{H(\hat\beta)^{-1}\}_{jj}}\). More generally, for a fixed vector \(a\), the estimated variance of the linear combination \(a^\mathsf{T}\hat\beta\) is \(\hat\phi\,a^\mathsf{T}H(\hat\beta)^{-1}a\). These are covariance calculations; how to form inferential tests using estimated dispersion is a separate question.

Common Mistakes

Multiplying standard errors by the dispersion. The covariance matrix is multiplied by \(\hat\phi\), so standard errors are multiplied by \(\sqrt{\hat\phi}\).
Mixing score scalings. If the score includes \(1/\phi\), its sensitivity and covariance are both \(H/\phi\). If that factor is removed, they are \(H\) and \(\phi H\). Combining the sensitivity from one definition with the covariance from the other gives an inconsistent sandwich calculation.
Claiming an exact finite-sample estimator covariance. The sandwich derivation gives the first-order covariance, and the plug-in matrix estimates it. A nonlinear estimating equation does not generally make this the estimator’s exact finite-sample covariance.
Assuming any variance change preserves coefficients. A common scalar cancels from the score equation; changes to relative observation weights need not.

Key takeaway: For a common dispersion parameter, the quasi-model covariance estimate is the ordinary model-based covariance multiplied by the estimated dispersion. This scales every variance and covariance by \(\hat\phi\), every standard error by \(\sqrt{\hat\phi}\), and leaves the common-score coefficient estimates unchanged.

Check Your Understanding

Use the definitions and derivations above to answer these questions.

If \(\hat\phi=4\), by what factor do the covariance matrix and standard errors change relative to their unscaled values?
For \(U_\phi=U_1/\phi\), what are the score sensitivity and score covariance under the assumptions of the theorem?
Why do \(U_\phi(\beta)=0\) and \(U_1(\beta)=0\) have the same solutions when \(\phi\) is a common positive scalar?
If a coefficient has unscaled variance estimate \(0.36\) and \(\hat\phi=2.25\), calculate its scaled variance and standard error.
Would a change that alters relative observation-specific weights necessarily preserve the coefficient estimates? Explain.