Leverage-Adjusted GEE Covariance Estimation

Why Adjust Cluster Contributions?

The ordinary robust GEE covariance estimator treats each cluster’s fitted score as an empirical contribution to the covariance estimate. But fitted residuals are not raw errors: estimation of the regression coefficients has already pulled them toward the fitted mean. This effect can be substantial when a cluster has high leverage, meaning its covariates contribute strongly to determining the fitted coefficients. With few clusters, a handful of such contributions can make the unadjusted sandwich covariance too small or unstable.

Fay–Graubard and CR2-style adjustments address this effect by inflating cluster residuals according to their leverage. The adjustment is matrix-valued when a cluster has multiple observations, and its construction depends on the covariance scale used to define leverage. It is a finite-sample correction to the covariance estimator, not a replacement for the GEE estimating equation or a guarantee of accurate inference.

Whitened Cluster Leverage

For cluster \(i\), let \(Y_i\) be the response vector, \(\mu_i\) the fitted mean vector, \(D_i=\partial\mu_i/\partial\beta^\mathsf{T}\) the mean-derivative matrix, and \(V_i\) the working covariance matrix. Assume \(V_i\) is positive definite and the total GEE bread matrix is nonsingular. Evaluate these quantities at the fitted coefficient vector when estimating covariance.

The whitened cluster score is \(s_i=W_i^\mathsf{T}r_i=D_i^\mathsf{T}V_i^{-1}(Y_i-\mu_i)\). The ordinary empirical sandwich covariance is

Because \(B\) is symmetric, the two bread factors here are also written \(B^{-\mathsf{T}}\) and \(B^{-1}\) in the general sandwich form. When \(W\) has full column rank, \(P\) is a symmetric idempotent projection: \(P^\mathsf{T}=P\), and \(P^2=WB^{-1}(W^\mathsf{T}W)B^{-1}W^\mathsf{T}=P\). Its diagonal blocks need not themselves be projections.

In particular, \(P_{ii}\) is positive semidefinite and has eigenvalues between zero and one, but it is generally not idempotent. The full identity \(P^2=P\) involves cross-cluster blocks as well as diagonal blocks. Also, \(\sum_i\operatorname{tr}(P_{ii})=\operatorname{tr}(P)=p\), where \(p\) is the number of regression coefficients. A large eigenvalue of \(P_{ii}\) indicates a direction in that cluster’s whitened response space that is strongly constrained by fitting.

Fay–Graubard and CR2-Style Adjustments

The central construction is to inflate the whitened residual block in directions where fitting has reduced its variance. For the usual working-covariance target, define a symmetric adjustment matrix for cluster \(i\) by

where \(m_i\) is the cluster size and the inverse square root is the symmetric matrix function: if \(P_{ii}=U_i\Lambda_iU_i^\mathsf{T}\), then \(A_i=U_i(I-\Lambda_i)^{-1/2}U_i^\mathsf{T}\). This is the standard whitened form of a Fay–Graubard-type leverage adjustment and of the basic CR2-style adjustment under a working-covariance target. The adjusted score and covariance estimator are

The factor \(A_i\) acts on the whitened residual, not directly on the unwhitened response residual unless the corresponding covariance transformations are included. Other CR2 implementations can use a specified target covariance and construct an adjustment to match that target. The formula above is the working-covariance version; the target and parameterization must be stated when comparing implementations.

Here is the justification. Consider the whitened linear model \(y=W\beta+\varepsilon\), with \(E[\varepsilon]=0\) and \(\operatorname{Var}(\varepsilon)=I\). The fitted residual vector is \(r=(I-P)y=(I-P)\varepsilon\), because \((I-P)W=0\). Its covariance is

Taking the \(i\)th diagonal block gives \(\operatorname{Var}(r_i)=I_{m_i}-P_{ii}\). Since \(A_i\) is symmetric and is a function of \(P_{ii}\),

This proves the stated residual-covariance property under the assumptions of the linearized model. In GEE, the mean model is generally nonlinear in its coefficients, and the covariance estimate is formed from observed cluster scores. Thus the same construction is a leverage correction motivated by local linearization; it does not make the finite-sample sandwich exact under arbitrary dependence or mean-model misspecification. \(\square\)

Worked Calculations

Conditions, Interpretation, and Common Pitfalls

The inverse square root exists only if \(I_{m_i}-P_{ii}\) is positive definite, or equivalently every eigenvalue of \(P_{ii}\) is strictly below one. If an eigenvalue equals one, the corresponding residual direction has zero fitted-residual variance in the linear projection calculation, and the stated adjustment is undefined. If an eigenvalue is close to one, the adjustment factor in that direction is large, so a small change in the residual can produce a large change in the covariance estimate. Numerical truncation or alternative regularization may be used in some software, but it changes the stated estimator and should be reported rather than silently treated as the same correction.

High leverage is especially concerning when there are few independent clusters, when cluster sizes are very unequal, or when one cluster contains covariate patterns that are poorly represented elsewhere. Under these conditions, the unadjusted robust standard error may be unreliable because the empirical meat is based on fitted residual scores that have been shrunk by estimation. Leverage adjustment can reduce this downward bias in suitable settings, but it cannot create independent information, repair a misspecified mean model, or ensure reliable tests with very few clusters. The correction’s performance also depends on the working covariance and target used to define the whitened geometry.

• Mistake: treating \(P_{ii}\) as a projection. The full \(P\) is idempotent, but a diagonal block \(P_{ii}\) is generally not. Its eigenvalues can be between zero and one without being only zero or one.
• Mistake: setting cross-cluster blocks to zero while keeping arbitrary diagonal blocks. The resulting block-diagonal matrix need not be idempotent and need not represent a fitted-value projection. Any full projection must satisfy all block equations implied by \(P^2=P\).
• Mistake: adjusting residuals but forgetting the design in the meat. The contribution is \(\tilde s_i\tilde s_i^\mathsf{T}\), with \(\tilde s_i=W_i^\mathsf{T}A_ir_i\); in a scalar example it includes \(w_i^2\).
• Mistake: treating a proposed residual block as freely selectable. A stacked fitted residual must satisfy \(Pr=0\), equivalently \(W^\mathsf{T}r=0\). One block equation alone does not establish compatibility.