Plug-In Dispersion and Correlation Estimates

Why Estimate the Nuisance Parameters?

A generalized estimating equation (GEE) models a marginal mean while using a working covariance to weight the cluster-level residuals. Its regression parameter is the main target, but the weights typically involve nuisance quantities: a dispersion parameter and one or more working-correlation parameters. These are usually not known in advance. A practical GEE therefore estimates them from residuals and substitutes the estimates into the next regression update and the covariance calculation.

The key question is whether this substitution changes the large-sample inference for the mean parameter. Under regularity conditions, consistent nuisance estimates preserve the first-order limit of the GEE estimator and its sandwich covariance. The reason is not that estimated weights have no effect in finite samples. Rather, the estimating function remains centered at the true mean parameter across working-correlation choices, so estimating the weights contributes no additional first-order term.

Residual Moments for Dispersion and Correlation

Let \(Y_i=(Y_{i1},\ldots,Y_{im_i})^\mathsf{T}\) denote the response vector for independent cluster \(i\), for \(i=1,\ldots,n\). The marginal mean model is \(\mu_i(\beta)=E[Y_i]\), with regression parameter \(\beta\), and \(D_i(\beta)=\partial\mu_i(\beta)/\partial\beta^\mathsf{T}\). Suppose the marginal variance of component \(j\) is \(\phi v_{ij}(\beta)\), where \(v_{ij}(\beta)\) is specified by the mean-variance relationship and \(\phi>0\) is a dispersion parameter. Write \(A_i(\beta)=\operatorname{diag}(v_{i1}(\beta),\ldots,v_{im_i}(\beta))\), and let \(R_i(\alpha)\) be a working correlation matrix indexed by nuisance parameter \(\alpha\). The working covariance is

V_i(\beta,\phi,\alpha) =\phi A_i(\beta)^{1/2}R_i(\alpha)A_i(\beta)^{1/2}.

For a provisional estimate \(\tilde\beta\), define the marginally standardized residual \(r_{ij}(\tilde\beta)=\{Y_{ij}-\mu_{ij}(\tilde\beta)\}/\sqrt{v_{ij}(\tilde\beta)}\). If the mean model is correct and the parameters are evaluated at their true values, then \(E[r_{ij}^2]=\phi\), and for distinct components \(j,k\) in the same cluster, \(E[r_{ij}r_{ik}]=\phi R_{i,jk}(\alpha)\). These moment identities motivate the estimators; they do not require the full joint response distribution to be specified.

Definition: A moment estimator of dispersion matches an empirical average of squared standardized residuals to its marginal expectation \(\phi\). A moment estimator of working correlation matches empirical products of distinct standardized residuals to their expectations \(\phi R_{i,jk}(\alpha)\). In practice the residuals use a fitted mean, and \(\phi\) and \(\alpha\) may be estimated jointly or successively.

A simple Pearson-type dispersion estimate is

\hat\phi =\frac{1}{N_n}\sum_{i=1}^n\sum_{j=1}^{m_i} r_{ij}(\hat\beta)^2,\qquad N_n=\sum_{i=1}^n m_i.

A degrees-of-freedom adjustment, such as replacing \(N_n\) by \(N_n-p\) for a \(p\)-dimensional \(\beta\), is sometimes used. It is not generally an exact unbiasedness correction for correlated, fitted residuals: fitting the mean changes their second moments, and the amount of change depends on the design and covariance. For large-sample consistency, the difference between these denominators is negligible when \(p\) is fixed and \(N_n\) grows.

For an exchangeable working correlation, \(R_{i,jk}(\alpha)=\rho\) for every distinct within-cluster pair. When there are at least some clusters with two or more observations, a direct moment estimate is

\hat\rho =\frac{\displaystyle\sum_{i=1}^n\sum_{j<k}r_{ij}(\hat\beta)r_{ik}(\hat\beta)} {\displaystyle\hat\phi\sum_{i=1}^n {m_i\choose 2}}.

The denominator scales the average residual cross-product by the estimated marginal dispersion. More generally, estimate \(\alpha\) by solving moment equations or minimizing a discrepancy between empirical residual products and \(\hat\phi R_{i,jk}(\alpha)\). The estimate must be kept in a parameter region for which each working correlation matrix is positive definite. If the true parameter is in the interior and the unconstrained moment estimate is consistent, a suitable constraint has no asymptotic effect; near a boundary, constraint and small-sample behavior require more care.

Worked Example: Pearson Dispersion and Exchangeable Correlation

Suppose three clusters each have two responses. At a provisional mean fit, the standardized residual pairs are \((1,-1)\), \((2,0)\), and \((-1,1)\). For this calculation, treat the fitted mean and variance functions as fixed. Calculate the Pearson dispersion estimate and the exchangeable moment estimate.

There are \(N_n=6\) residuals. Their squared values sum to \(1+1+4+0+1+1=8\), so \(\hat\phi=8/6=4/3\).

There is one distinct pair per cluster. The cross-products are \(-1\), \(0\), and \(-1\), summing to \(-2\); the number of pairs is \(3\). Therefore,

\hat\rho=\frac{-2}{(4/3)(3)}=-\frac{1}{2}.

The estimate is in the valid exchangeable range for clusters of size two, \(-1<\rho<1\). This is a residual-moment calculation, not a claim that the estimate is unbiased after fitting a mean model. With estimated means, residual dependence on \(\hat\beta\) must be accounted for in the large-sample argument.

When Are the Moment Estimates Consistent?

A useful basic setting has independent clusters, a fixed-dimensional mean parameter, uniformly bounded cluster sizes, and a correctly specified marginal mean and variance function. Assume the cluster observations have finite moments sufficient for laws of large numbers for the residual squares and cross-products. Also assume \(\hat\beta\xrightarrow{p}\beta_0\), and that the mean and variance functions are continuous and have bounded derivatives in a neighborhood of \(\beta_0\). These conditions make replacing \(\beta_0\) by \(\hat\beta\) negligible in averages of the residual moments.

For example, with bounded cluster sizes, the difference between a fitted and true standardized residual is \(O_p(\|\hat\beta-\beta_0\|)\) within each cluster under local smoothness. The difference between their squares is bounded by the residual difference times the sum of the absolute residuals. Finite second moments and averaging over independent clusters then make the difference in the empirical average converge to zero in probability. The same argument applies to cross-products. Consequently, the Pearson estimate converges to \(\phi_0\), and a well-identified solution of the correlation moment equations converges to the true working-correlation parameter. Identification matters: if the design supplies no within-cluster pairs, for example, an exchangeable correlation parameter cannot be estimated from these moments.

Worked Example: Why a Correlation Moment Needs Pairs

Suppose every cluster has size one. The dispersion moment can still use the standardized residual squares. But the exchangeable correlation formula contains \(\sum_i {m_i\choose 2}\) in its denominator. Here each term is \({1\choose 2}=0\), so there are no within-cluster products and the denominator is zero.

This is not a numerical inconvenience that can be fixed by choosing a different divisor. With one observation per cluster, the data contain no within-cluster pair information about correlation. The correlation parameter is not identified by these residual moments. An independence working correlation is available, but a nonzero within-cluster correlation cannot be estimated from this design.

Why Plugging In Correlation Preserves First-Order GEE Inference

Write the cluster estimating function at working parameter \(\alpha\) as

\psi_i(\beta,\alpha) =D_i(\beta)^\mathsf{T}V_i(\beta,\phi,\alpha)^{-1} \{Y_i-\mu_i(\beta)\},\qquad U_n(\beta,\alpha)=\sum_{i=1}^n\psi_i(\beta,\alpha).

For this argument, treat dispersion as fixed or consistently estimated, and suppose that for every \(\alpha\) in a neighborhood of \(\alpha_0\), \(E[\psi_i(\beta_0,\alpha)]=0\). This centering follows when the marginal mean is correct, since \(E[Y_i-\mu_i(\beta_0)]=0\); it does not require the working covariance to be the true covariance. Assume the relevant derivatives and second moments are finite, the matrices defining the estimating function are nonsingular, and the usual laws of large numbers, central limit theorem, and stochastic equicontinuity hold. In particular, assume \(n^{-1/2}\partial U_n(\beta_0,\alpha)/\partial\alpha\) is bounded in probability uniformly near \(\alpha_0\). These conditions are stated to make the nuisance-substitution step precise.

Theorem: Suppose the GEE based on the fixed working parameter \(\alpha_0\) has a consistent root \(\tilde\beta\), with \(\sqrt n(\tilde\beta-\beta_0)\xrightarrow{d}N(0,H^{-1}JH^{-\mathsf T})\), where \(H\) is the probability limit of the normalized negative derivative of \(U_n\) with respect to \(\beta\), and \(J\) is the probability limit of the normalized sum of cluster score outer products. Under the centering and regularity conditions above, if \(\hat\alpha\xrightarrow{p}\alpha_0\), the consistent root \(\hat\beta\) of \(U_n(\beta,\hat\alpha)=0\) has the same first-order limit.

Compare the estimating functions at the true mean parameter.
Apply a Taylor expansion in \(\alpha\) between \(\alpha_0\) and \(\hat\alpha\). The difference, divided by \(\sqrt n\), is bounded in probability by a term of order \(\|\hat\alpha-\alpha_0\|\), and hence converges to zero.

Expand in the mean parameter.
For either root, expand \(U_n\) around \(\beta_0\). The normalized negative derivative converges to \(H\), which is nonsingular, so the root's scaled error is determined by \(H^{-1}n^{-1/2}U_n(\beta_0,\alpha_0)\), up to a term converging to zero.

Use the same limiting score.
The score difference from Step 1 is negligible on the \(\sqrt n\) scale. Thus both roots have the same limiting distribution, with covariance \(H^{-1}JH^{-\mathsf T}\).

For completeness, the first step follows from the mean-value expansion

\frac{U_n(\beta_0,\hat\alpha)-U_n(\beta_0,\alpha_0)}{\sqrt n} = \left\{\frac{1}{\sqrt n} \frac{\partial U_n(\beta_0,\bar\alpha)}{\partial\alpha^\mathsf{T}}\right\} (\hat\alpha-\alpha_0),

for an intermediate \(\bar\alpha\). The bracketed term is \(O_p(1)\) uniformly near \(\alpha_0\), and the last factor is \(o_p(1)\), so the product is \(o_p(1)\). The root expansion and the assumed central limit theorem then give the stated common limit. This proves the result. \(\square\)

The same reasoning applies to a consistent dispersion estimate when dispersion affects the estimating weights smoothly and the estimating function remains centered at \(\beta_0\) for nearby dispersion values. If dispersion is a common scalar factor in \(V_i\), there is an even simpler fact: it multiplies the entire GEE by \(1/\phi\), so it does not change the root at all. Correlation generally changes relative weights and can change the finite-sample coefficient estimate, but consistent plug-in estimation does not change its first-order limit under the theorem’s conditions.

Worked Example: A Common Dispersion Factor Cancels From the GEE

Suppose \(V_i(\beta,\phi,\alpha)=\phi W_i(\beta,\alpha)\), where \(W_i\) does not contain \(\phi\). Then \(V_i^{-1}=\phi^{-1}W_i^{-1}\), and the estimating equation becomes

U_n(\beta,\phi,\alpha) =\frac{1}{\phi}\sum_{i=1}^n D_i^\mathsf{T}W_i^{-1}(Y_i-\mu_i).

For any positive \(\phi\), \(U_n(\beta,\phi,\alpha)=0\) exactly when the sum on the right without \(1/\phi\) is zero. Thus replacing \(\phi\) by \(\hat\phi>0\) leaves the coefficient root unchanged, not merely unchanged to first order.

The covariance calculation still uses scale: the bread based on this score is multiplied by \(1/\phi\), and the score outer-product meat by \(1/\phi^2\). In the sandwich, the inverse bread appears on both sides, so these factors cancel. A model-based covariance does not generally have this cancellation: it uses the assumed covariance itself and therefore depends on \(\phi\). The distinction is why dispersion estimation remains relevant even when it does not affect the coefficient root.

Plug-In Sandwich Covariance and Practical Checks

Let \(\hat\psi_i=\psi_i(\hat\beta,\hat\alpha)\), with any dispersion estimate included if needed. The empirical sandwich uses a normalized bread and meat, for example

\hat H=-\frac{1}{n}\sum_{i=1}^n \frac{\partial\psi_i(\hat\beta,\hat\alpha)}{\partial\beta^\mathsf{T}}, \qquad \hat J=\frac{1}{n}\sum_{i=1}^n\hat\psi_i\hat\psi_i^\mathsf{T}, \qquad \widehat{\operatorname{Var}}(\hat\beta) =\frac{1}{n}\hat H^{-1}\hat J\hat H^{-\mathsf{T}}.

If \(\hat\beta\) and the nuisance estimates converge to their probability limits, and the bread, score outer products, and their derivatives satisfy the corresponding continuity and law-of-large-numbers conditions, then \(\hat H\) and \(\hat J\) converge to their population limits. The plug-in sandwich therefore estimates the first-order covariance. This does not make it reliable in every finite sample: the number of independent clusters, influential clusters, and the quality of the mean model still matter.

A useful distinction is that a correctly specified mean is what centers the GEE. The working correlation controls weighting and efficiency, not that centering property. If the correlation moments are noisy or converge to a different admissible working parameter, a sandwich can still support inference when the mean estimating equation remains centered and the required limits exist. That broader issue should not be confused with the particular claim here: consistent nuisance estimates preserve the first-order limit associated with their target working parameter.

Common Pitfall

Wrong reasoning: “The correlation is estimated from the same residuals as the regression, so its estimation must add a first-order variance term to \(\hat\beta\).”
Correction: Dependence alone does not establish a first-order contribution. Under the stated regularity conditions, the estimating function is centered at \(\beta_0\) for each nearby working correlation. The derivative with respect to the nuisance parameter is therefore centered as well, and a consistent nuisance substitution changes the score by \(o_p(\sqrt n)\). The conclusion can fail if the mean-centering property, consistency, smoothness, or nonsingularity conditions fail.
Another error: treating a Pearson dispersion estimate with a degrees-of-freedom divisor as exactly unbiased in every correlated GEE. Fitted, correlated residuals do not generally have the simple independent-regression residual degrees-of-freedom calculation.

Key takeaway: Residual squares estimate marginal dispersion, and within-cluster residual products estimate working correlation. With identifiable moments and regularity conditions, consistent nuisance estimates can be plugged into the GEE and sandwich without changing first-order inference. A common scalar dispersion factor cancels from the GEE root, although it remains relevant to model-based covariance calculations.

Check Your Understanding

Answer each question by stating the relevant moment identity or asymptotic condition.

Why does the expectation of a squared marginally standardized residual identify \(\phi\)?
What information is missing when every cluster has size one and an exchangeable correlation is to be estimated from residual products?
Under what centering property does consistent plug-in estimation of \(\alpha\) leave the first-order GEE limit unchanged?
Why can a common dispersion factor leave the coefficient root unchanged but still matter for a model-based covariance?
Why should a Pearson estimate with a degrees-of-freedom adjustment not automatically be described as exactly unbiased for correlated fitted residuals?