Working Covariance as Residual Weighting

Why the Working Covariance Changes Residual Weighting

In a generalized estimating equation, a cluster contributes through its residual vector and the inverse of its working covariance. The inverse does more than adjust for different marginal variances: when responses are associated, it also combines residual information across observations within the cluster. Understanding that combination makes the GEE contribution more interpretable and helps catch matrix-ordering errors.

Consider $K$ independent clusters. Cluster $i$ contains $m_i$ responses, collected in the vector $\mathbf Y_i$, with marginal mean vector $\boldsymbol\mu_i=E[\mathbf Y_i]$. Write the residual vector as $\mathbf r_i=\mathbf Y_i-\boldsymbol\mu_i$. Let $\mathbf D_i=\partial\boldsymbol\mu_i/\partial\boldsymbol\beta^T$ be the $m_i$ by $p$ matrix of mean derivatives with respect to the $p$-dimensional parameter $\boldsymbol\beta$. The cluster GEE contribution has the form $\mathbf D_i^T\mathbf V_i^{-1}\mathbf r_i$, where $\mathbf V_i$ is the working covariance matrix for that cluster.

The covariance factorization separates two sources of scale: the marginal variance of each response and the working correlation among responses. That separation lets us see what each part does to the residuals before they are combined with the mean derivatives.

Factorization and Standardized Residuals

Definition: Suppose the working marginal variance for response $j$ in cluster $i$ is $\nu_{ij}>0$. Define $\mathbf A_i=\operatorname{diag}(\nu_{i1},\ldots,\nu_{im_i})$ and $\mathbf S_i=\mathbf A_i^{1/2}=\operatorname{diag}(\sqrt{\nu_{i1}},\ldots,\sqrt{\nu_{im_i}})$. A working correlation matrix $\mathbf R_i$ is symmetric, positive definite, and has diagonal entries equal to one. The corresponding working covariance is $$ \mathbf V_i=\mathbf S_i\mathbf R_i\mathbf S_i. $$ The standardized residual vector is $\mathbf e_i=\mathbf S_i^{-1}\mathbf r_i$.

In a generalized linear model, $\nu_{ij}$ is often specified from the marginal mean and variance function, with any dispersion factor included when appropriate. The factorization says that $\mathbf S_i$ places each response on its marginal standard-deviation scale, while $\mathbf R_i$ describes the working association on the standardized scale. The entries of $\mathbf R_i$ are correlations, not covariances.

Positive marginal variances make $\mathbf S_i$ invertible. Positive definiteness of $\mathbf R_i$ makes $\mathbf V_i$ invertible, so the inverse-covariance GEE contribution is defined. The same assumptions also ensure that the quadratic form induced by the inverse covariance is positive for every nonzero residual vector.

Result: Under these conditions, $$ \mathbf V_i^{-1}=\mathbf S_i^{-1}\mathbf R_i^{-1}\mathbf S_i^{-1}, \qquad \mathbf V_i^{-1}\mathbf r_i=\mathbf S_i^{-1}\mathbf R_i^{-1}\mathbf e_i. $$ Thus the standardized residuals are first combined by the inverse working correlation and then rescaled by the inverse marginal standard deviations.

The order of the factors matters. Matrix multiplication is not generally commutative: one cannot move $\mathbf S_i^{-1}$ through $\mathbf R_i^{-1}$ or replace the inverse by a componentwise reciprocal of the covariance entries. The factorization gives an exact way to compute the inverse-vector product without losing track of which scale is being applied.

Worked Example: Two Responses With Unequal Marginal Scales

Suppose a cluster has marginal standard deviations 2 and 1, and working correlation 0.5. Thus $$ \mathbf S=\begin{pmatrix}2&0\\0&1\end{pmatrix}, \qquad \mathbf R=\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}. $$ The working covariance and its inverse are $$ \mathbf V=\mathbf S\mathbf R\mathbf S =\begin{pmatrix}4&1\\1&1\end{pmatrix}, \qquad \mathbf V^{-1} =\frac{1}{3}\begin{pmatrix}1&-1\\-1&4\end{pmatrix}. $$ The inverse follows directly from the two-by-two inverse formula: the determinant of $\mathbf V$ is $4(1)-1(1)=3$.

Take observed residuals $\mathbf r=(2,1)^T$. Standardizing gives $\mathbf e=\mathbf S^{-1}\mathbf r=(1,1)^T$. Also, $$ \mathbf R^{-1}=\frac{1}{0.75}\begin{pmatrix}1&-0.5\\-0.5&1\end{pmatrix} =\begin{pmatrix}4/3&-2/3\\-2/3&4/3\end{pmatrix}. $$ Therefore $\mathbf R^{-1}\mathbf e=(2/3,2/3)^T$, and $$ \mathbf S^{-1}\mathbf R^{-1}\mathbf e =\begin{pmatrix}1/3\\2/3\end{pmatrix} =\mathbf V^{-1}\mathbf r. $$ The two calculations agree. Each residual is first expressed in standard-deviation units, and the inverse correlation then combines those standardized values.

Proof of the Inverse and Its Residual Action

The inverse formula follows by multiplying in the correct order. Since $\mathbf S_i$ and $\mathbf R_i$ are invertible, the proposed inverse exists. Multiplying it by the covariance gives $$ (\mathbf S_i\mathbf R_i\mathbf S_i) (\mathbf S_i^{-1}\mathbf R_i^{-1}\mathbf S_i^{-1}) =\mathbf S_i\mathbf R_i\mathbf R_i^{-1}\mathbf S_i^{-1} =\mathbf I. $$ Multiplication in the reverse order also gives the identity: $$ (\mathbf S_i^{-1}\mathbf R_i^{-1}\mathbf S_i^{-1}) (\mathbf S_i\mathbf R_i\mathbf S_i) =\mathbf S_i^{-1}\mathbf R_i^{-1}\mathbf R_i\mathbf S_i =\mathbf I. $$ So the proposed matrix is the inverse of $\mathbf V_i$. Applying it to the residual vector and using $\mathbf e_i=\mathbf S_i^{-1}\mathbf r_i$ yields $$ \mathbf V_i^{-1}\mathbf r_i =\mathbf S_i^{-1}\mathbf R_i^{-1}\mathbf S_i^{-1}\mathbf r_i =\mathbf S_i^{-1}\mathbf R_i^{-1}\mathbf e_i. $$ This proves both identities. $\square$

To see how this acts in the estimating equation, define $\mathbf q_i=\mathbf R_i^{-1}\mathbf e_i$. Then $$ \mathbf D_i^T\mathbf V_i^{-1}\mathbf r_i =\mathbf D_i^T\mathbf S_i^{-1}\mathbf q_i. $$ If $\mathbf d_{ij}^T$ is row $j$ of $\mathbf D_i$, the contribution can also be written as $$ \sum_{j=1}^{m_i}\frac{\mathbf d_{ij}}{\sqrt{\nu_{ij}}}\,q_{ij}. $$ Here $q_{ij}$ is generally a combination of standardized residuals from several observations, not just the residual from observation $j$. The mean-derivative row and the marginal scale then determine how that combined residual enters the parameter estimating function.

Worked Example: A Cluster GEE Contribution

Suppose there are two observations, with working standard deviations 1 and 2 and working correlation 0.25. Let the residual vector be $\mathbf r=(1,-2)^T$, and let the mean-derivative matrix be $$ \mathbf D=\begin{pmatrix}1&0\\1&1\end{pmatrix}. $$ The standardized residuals are $\mathbf e=(1,-1)^T$. Since $1-0.25^2=15/16$, $$ \mathbf R^{-1} =\frac{16}{15}\begin{pmatrix}1&-0.25\\-0.25&1\end{pmatrix}. $$ Thus $$ \mathbf q=\mathbf R^{-1}\mathbf e =\frac{16}{15}\begin{pmatrix}1+0.25\\-0.25-1\end{pmatrix} =\begin{pmatrix}4/3\\-4/3\end{pmatrix}. $$ Rescaling by $\mathbf S^{-1}=\operatorname{diag}(1,1/2)$ gives $\mathbf V^{-1}\mathbf r=(4/3,-2/3)^T$. The cluster contribution is therefore $$ \mathbf D^T\mathbf V^{-1}\mathbf r =\begin{pmatrix}1&1\\0&1\end{pmatrix} \begin{pmatrix}4/3\\-2/3\end{pmatrix} =\begin{pmatrix}2/3\\-2/3\end{pmatrix}. $$ In this calculation, the inverse correlation first mixes the two standardized residuals; the standard-deviation rescaling and derivative matrix then determine their effect on the two components of the estimating function.

The Residual Quadratic Form

The same factorization describes the size assigned to a residual vector by inverse-covariance weighting. This quadratic form is not itself the GEE contribution, which is a vector, but it helps explain the geometry of the weights.

Result: For any residual vector $\mathbf r_i$, $$ \mathbf r_i^T\mathbf V_i^{-1}\mathbf r_i =\mathbf e_i^T\mathbf R_i^{-1}\mathbf e_i. $$ If $\mathbf R_i=\mathbf L_i\mathbf L_i^T$ is a Cholesky factorization, then $$ \mathbf e_i^T\mathbf R_i^{-1}\mathbf e_i =\left\|\mathbf L_i^{-1}\mathbf e_i\right\|^2. $$

For the first equality, substitute $\mathbf V_i^{-1}=\mathbf S_i^{-1}\mathbf R_i^{-1}\mathbf S_i^{-1}$. Because $\mathbf S_i$ is symmetric and $\mathbf e_i=\mathbf S_i^{-1}\mathbf r_i$, the result is $$ \mathbf r_i^T\mathbf S_i^{-1}\mathbf R_i^{-1}\mathbf S_i^{-1}\mathbf r_i =(\mathbf S_i^{-1}\mathbf r_i)^T\mathbf R_i^{-1}(\mathbf S_i^{-1}\mathbf r_i) =\mathbf e_i^T\mathbf R_i^{-1}\mathbf e_i. $$ For the second equality, positive definiteness of $\mathbf R_i$ guarantees a Cholesky factorization with nonsingular $\mathbf L_i$. Since $\mathbf R_i^{-1}=\mathbf L_i^{-T}\mathbf L_i^{-1}$, $$ \mathbf e_i^T\mathbf R_i^{-1}\mathbf e_i =\mathbf e_i^T\mathbf L_i^{-T}\mathbf L_i^{-1}\mathbf e_i =(\mathbf L_i^{-1}\mathbf e_i)^T(\mathbf L_i^{-1}\mathbf e_i) =\|\mathbf L_i^{-1}\mathbf e_i\|^2. $$ This is nonnegative and is zero only when $\mathbf e_i=\mathbf 0$, equivalently $\mathbf r_i=\mathbf 0$. $\square$

The Cholesky expression shows that correlation adjustment can be understood as a transformation of the standardized residual vector followed by an ordinary squared length. It does not mean that GEE requires a multivariate likelihood for the cluster: the working covariance specifies weights in the estimating equation, and the algebra above describes those weights.

Worked Example: A Three-Observation Correlation Pattern

Consider the working correlation matrix $$ \mathbf R=\begin{pmatrix}1&0.5&0\\0.5&1&0.5\\0&0.5&1\end{pmatrix}. $$ Its determinant is $1-0.25-0.25=0.5$, and its leading principal minors are $1$, $0.75$, and $0.5$, all positive. It is therefore positive definite. Its inverse is $$ \mathbf R^{-1} =\begin{pmatrix}1.5&-1&0.5\\-1&2&-1\\0.5&-1&1.5\end{pmatrix}. $$ Suppose $\mathbf S=\operatorname{diag}(1,2,0.5)$ and $\mathbf r=(1,2,0)^T$. Then $\mathbf e=\mathbf S^{-1}\mathbf r=(1,1,0)^T$, so $$ \mathbf q=\mathbf R^{-1}\mathbf e =\begin{pmatrix}1.5-1\\-1+2\\0.5-1\end{pmatrix} =\begin{pmatrix}0.5\\1\\-0.5\end{pmatrix}. $$ The inverse-covariance-weighted residual is $$ \mathbf V^{-1}\mathbf r=\mathbf S^{-1}\mathbf q =\begin{pmatrix}0.5\\0.5\\-1\end{pmatrix}. $$ For a scalar mean parameter whose derivative across the three observations is $(1,2,1)^T$, the cluster contribution is $$ (1,2,1)\begin{pmatrix}0.5\\0.5\\-1\end{pmatrix} =0.5+1-1=0.5. $$ Although the third raw residual is zero, its weighted residual is $-1$: the inverse correlation combines information across the cluster. A response’s weighted residual need not have the same sign as its raw residual.

Common Pitfall and Practical Interpretation

A frequent mistake is to assume that inverse-covariance weighting assigns a separate positive weight to each residual. That is true only in special cases where the working correlation is diagonal. In general, $\mathbf R_i^{-1}$ has off-diagonal entries, and $\mathbf q_i=\mathbf R_i^{-1}\mathbf e_i$ mixes residuals. Its entries may be negative, so it is not generally meaningful to describe the operation as multiplying each residual by an independent positive weight.

A related error is to treat the entries of $\mathbf V_i^{-1}$ as reciprocal covariance entries. Matrix inversion involves the whole matrix, including its off-diagonal structure. The safe calculation is to use $\mathbf V_i^{-1}=\mathbf S_i^{-1}\mathbf R_i^{-1}\mathbf S_i^{-1}$, or to solve the corresponding linear system, while preserving the order of the factors.

The working covariance therefore affects how a cluster’s residual evidence is combined in its estimating function. Marginal standard deviations set the units of the residuals; the inverse correlation adjusts their joint contribution. These are weights in an estimating equation, not a claim that the specified covariance is the true joint distribution. Keeping that distinction clear helps make sense of both the algebra and the interpretation of a GEE fit.

Check Your Understanding

Use the factorization and the residual-weighting interpretation to answer these questions.

State the conditions under which $\mathbf V_i=\mathbf S_i\mathbf R_i\mathbf S_i$ is invertible.
Why must the standardized residual vector be $\mathbf S_i^{-1}\mathbf r_i$, rather than $\mathbf S_i\mathbf r_i$?
For a non-diagonal working correlation, why can a weighted residual have a different sign from its raw residual?
Explain how $\mathbf D_i^T\mathbf S_i^{-1}\mathbf R_i^{-1}\mathbf e_i$ separates residual combination from mean-derivative weighting.
What does the Cholesky form of the residual quadratic show, and what does it not imply about the GEE model?