Working Correlation Matrices and Standardized Residuals

Why Separate Marginal Variation from Association?

In a generalized estimating equation, responses within a cluster can have different marginal variances and can also be associated with one another. These are distinct features of the covariance model. A response with a larger marginal variance has a larger natural scale of fluctuation; association describes how responses vary together after their scales are taken into account.

The working covariance makes this separation explicit. It combines the marginal variance specified by the mean model with a working correlation matrix. Standardized residuals divide each response residual by its marginal standard deviation. This removes differences in marginal scale, but it does not make the standardized residuals independent.

Working Covariance from Marginal Variances

Let \(Y_i=(Y_{i1},\ldots,Y_{in_i})^\mathsf{T}\) be the response vector for cluster \(i\), where \(i=1,\ldots,m\) and \(n_i\) is the cluster size. Write \(\mu_{ij}=E[Y_{ij}]\) for the marginal mean of response \(j\) in cluster \(i\). Suppose the marginal variance model gives \(\operatorname{Var}(Y_{ij})=\sigma_{ij}^2=\phi v(\mu_{ij})\), where \(\phi>0\) is a dispersion parameter and \(v(\mu_{ij})>0\) is the variance function evaluated at the mean. The formulas below also apply when each marginal variance is specified directly as a positive number \(\sigma_{ij}^2\).

Definition: The diagonal matrix of marginal standard deviations for cluster \(i\) is \(D_i=\operatorname{diag}(\sigma_{i1},\ldots,\sigma_{in_i})\). A working correlation matrix \(R_i(\alpha)\) is a symmetric matrix with diagonal entries equal to 1, indexed by the responses in cluster \(i\), and depending, if needed, on working parameter(s) \(\alpha\). The cluster working covariance is \(V_i=D_iR_i(\alpha)D_i\). For the usual GEE inverse \(V_i^{-1}\) to exist, \(R_i(\alpha)\) is taken to be positive definite.

Because \(D_i\) is diagonal, the entry in row \(j\), column \(k\) of the working covariance is \(\sigma_{ij}\sigma_{ik}R_{i,jk}(\alpha)\). In particular, its diagonal entry is \(\sigma_{ij}^2\), since \(R_{i,jj}=1\). Thus the diagonal of \(V_i\) preserves the marginal variances, while its off-diagonal entries scale the working correlations by the corresponding marginal standard deviations.

V_{i,jk}=\sigma_{ij}\sigma_{ik}R_{i,jk}(\alpha), \qquad V_{i,jj}=\sigma_{ij}^2.

The word working matters. The factorization defines the covariance used in the GEE calculation; it does not by itself assert that the true response covariance has this form. If the working correlation is misspecified, the resulting \(V_i\) is still the chosen working covariance, not a claim that the true pairwise correlations equal its entries after scaling. The marginal variance model and the working association structure are separate inputs.

Worked Example: Scaling a Gaussian Working Correlation

Consider two responses in one cluster with marginal means \(\mu_{i1}=10\) and \(\mu_{i2}=20\), and Gaussian marginal variances \(4\) and \(9\). Their marginal standard deviations are \(2\) and \(3\), so \(D_i=\operatorname{diag}(2,3)\). Suppose the chosen working correlation matrix is \(R_i=\begin{pmatrix}1&0.4\\0.4&1\end{pmatrix}\).

First multiply \(D_i\) by \(R_i\), then by \(D_i\) on the right:

V_i = \begin{pmatrix}2&0\\0&3\end{pmatrix} \begin{pmatrix}1&0.4\\0.4&1\end{pmatrix} \begin{pmatrix}2&0\\0&3\end{pmatrix} = \begin{pmatrix}4&2.4\\2.4&9\end{pmatrix}.

The diagonal entries are the specified marginal variances. The off-diagonal covariance is \(2(3)(0.4)=2.4\), and its implied correlation is \(2.4/(\sqrt{4}\sqrt{9})=2.4/6=0.4\). This checks both parts of the scaling.

If the observed residuals \(Y_{ij}-\mu_{ij}\) are \(4\) and \(-3\), the standardized residuals are \(4/2=2\) and \(-3/3=-1\). This observed pair is a pair of residual values; it is not itself a correlation of \(0.4\).

Standardized Residuals and Their Covariance

Define the marginally standardized residual for response \(j\) in cluster \(i\) by \(r_{ij}=(Y_{ij}-\mu_{ij})/\sigma_{ij}\). Collecting these values gives the residual vector \(r_i=D_i^{-1}(Y_i-\mu_i)\), where \(\mu_i=(\mu_{i1},\ldots,\mu_{in_i})^\mathsf{T}\). These are often called standardized Pearson residuals when the marginal standard deviations come from the GLM variance model.

Theorem: Suppose \(E[Y_i]=\mu_i\), every \(\sigma_{ij}>0\), and the working covariance is \(V_i=D_iR_iD_i\). If the covariance of \(Y_i\) is \(V_i\), then \(E[r_i]=0\) and \(\operatorname{Cov}(r_i)=R_i\). More generally, if the true covariance is \(\Sigma_i\), then \(\operatorname{Cov}(r_i)=D_i^{-1}\Sigma_iD_i^{-1}\); this equals \(R_i\) only when the true covariance has the stated working form.

Proof. Since \(D_i\) is invertible, linearity of expectation gives \(E[r_i]=D_i^{-1}(E[Y_i]-\mu_i)=0\). For any deterministic matrix \(A\) and random vector \(Z\) with finite second moments, \(\operatorname{Cov}(AZ)=A\operatorname{Cov}(Z)A^\mathsf{T}\). Apply this identity with \(A=D_i^{-1}\) and \(Z=Y_i-\mu_i\). The covariance of \(Y_i-\mu_i\) is the covariance of \(Y_i\), so \(\operatorname{Cov}(r_i)=D_i^{-1}\Sigma_i(D_i^{-1})^\mathsf{T}\). The diagonal matrix \(D_i^{-1}\) is symmetric, giving the general formula. Under the stated working covariance model, substitute \(\Sigma_i=V_i=D_iR_iD_i\):

\operatorname{Cov}(r_i) =D_i^{-1}(D_iR_iD_i)D_i^{-1} =R_i.

The calculation uses only means, covariances, and linear transformations; it does not require the response vector to be normally distributed. \(\square\)

This result describes a covariance, not a realized residual vector. Even when the working covariance is correct, one cluster’s standardized residuals will not generally have empirical correlation exactly equal to the corresponding entry of \(R_i\). The theorem concerns the distribution across repeated realizations under the covariance model.

Worked Example: Standardizing Binary Responses

Suppose two binary responses have marginal means \(0.2\) and \(0.5\). A Bernoulli variable with mean \(\mu\) has variance \(\mu(1-\mu)\), so the marginal variances are \(0.2(0.8)=0.16\) and \(0.5(0.5)=0.25\). The standard deviations are \(0.4\) and \(0.5\). Choose the working correlation \(R_i=\begin{pmatrix}1&0.3\\0.3&1\end{pmatrix}\).

The working covariance is therefore

V_i = \begin{pmatrix}0.4&0\\0&0.5\end{pmatrix} \begin{pmatrix}1&0.3\\0.3&1\end{pmatrix} \begin{pmatrix}0.4&0\\0&0.5\end{pmatrix} = \begin{pmatrix}0.16&0.06\\0.06&0.25\end{pmatrix},

because the off-diagonal covariance is \(0.4(0.5)(0.3)=0.06\). For the observed outcome pair \((1,0)\), the raw residuals are \(1-0.2=0.8\) and \(0-0.5=-0.5\). Dividing by their own standard deviations gives \(r_i=(0.8/0.4,-0.5/0.5)^\mathsf{T}=(2,-1)^\mathsf{T}\).

The value \(-2\) obtained by multiplying these two observed residuals is not an estimate of the correlation. It is one realized product. The working correlation \(0.3\) is the off-diagonal entry of the assumed covariance structure after marginal standardization.

Positive Definiteness and Inversion

The GEE uses inverse working covariances, so it is important that the scaling preserve positive definiteness. Here positive definiteness means that a symmetric matrix \(M\) satisfies \(x^\mathsf{T}Mx>0\) for every nonzero vector \(x\) of the appropriate dimension.

Theorem: If every marginal standard deviation in \(D_i\) is positive, then \(V_i=D_iR_iD_i\) is positive definite if and only if \(R_i\) is positive definite. When these matrices are positive definite, \(V_i^{-1}=D_i^{-1}R_i^{-1}D_i^{-1}\).

Proof. For any vector \(x\), \(x^\mathsf{T}V_ix=x^\mathsf{T}D_iR_iD_ix=(D_ix)^\mathsf{T}R_i(D_ix)\). Since all diagonal entries of \(D_i\) are positive, \(D_i\) is invertible, and \(D_ix\ne 0\) whenever \(x\ne 0\). Therefore, if \(R_i\) is positive definite, the last quadratic form is positive for every nonzero \(x\), and \(V_i\) is positive definite. Conversely, if \(V_i\) is positive definite, take any nonzero vector \(z\) and set \(x=D_i^{-1}z\), which is nonzero. Then \(z^\mathsf{T}R_iz=x^\mathsf{T}V_ix>0\), so \(R_i\) is positive definite. Finally, multiplication verifies the inverse formula:

(D_iR_iD_i)(D_i^{-1}R_i^{-1}D_i^{-1})=I.

The factors in the inverse appear in reverse order, as required for a matrix product. \(\square\)

Worked Example: Checking the Inverse by Scaling

Use the Gaussian example with \(D_i=\operatorname{diag}(2,3)\) and \(R_i=\begin{pmatrix}1&0.4\\0.4&1\end{pmatrix}\). The determinant of \(R_i\) is \(1-0.4^2=0.84>0\); its first diagonal entry is positive, so this symmetric \(2\times2\) matrix is positive definite. Its inverse is

R_i^{-1} = \frac{1}{0.84} \begin{pmatrix}1&-0.4\\-0.4&1\end{pmatrix} = \begin{pmatrix}1.190476&-0.476190\\-0.476190&1.190476\end{pmatrix}.

Since \(D_i^{-1}=\operatorname{diag}(0.5,1/3)\), the covariance inverse is

V_i^{-1} =D_i^{-1}R_i^{-1}D_i^{-1} \approx \begin{pmatrix}0.297619&-0.079365\\-0.079365&0.132275\end{pmatrix}.

The unequal marginal standard deviations affect the inverse on both sides. Using \(R_i^{-1}\) alone would omit the marginal scaling and would not be the inverse of \(V_i\).

What Standardization Does—and Does Not Do

Standardization puts residuals measured on different marginal scales into units of their own standard deviations. If the mean and marginal variance models are appropriate, the standardized residuals have mean zero and unit marginal variances. Their off-diagonal covariances, however, remain. Under the working covariance model, those covariances are the entries of \(R_i\).

This separation is useful in checking a model: unusually large standardized residuals may indicate observations far from their fitted marginal means relative to their modeled variability, while patterns across responses within clusters may indicate association not represented well by the working structure. Such checks are descriptive; residuals calculated from fitted parameters are not independent observations of a correlation matrix.

Common Mistakes

Scaling on only one side. Writing \(D_iR_i\) instead of \(D_iR_iD_i\) generally gives a nonsymmetric matrix, which cannot be a covariance matrix. Each covariance entry needs both marginal standard deviations.
Calling the working correlation the true correlation. It is the assumed association structure used in the GEE. It equals the true standardized covariance only if that covariance model is correct.
Assuming standardization removes dependence. Dividing by marginal standard deviations changes the units; it does not force off-diagonal covariances to zero.
Reading a residual product as a correlation. One pair of observed standardized residuals does not determine the population or working correlation.

Key takeaway: The factorization \(V_i=D_iR_iD_i\) assigns marginal scales to \(D_i\) and within-cluster association to \(R_i\). Standardized residuals express departures from fitted means in marginal standard-deviation units; their covariance reflects association rather than marginal scale.

Check Your Understanding

Use the definitions and results above to answer these questions.

If two responses have standard deviations \(2\) and \(5\), and their working correlation is \(0.2\), what is their working covariance?
State the covariance of the standardized residual vector if the true response covariance is \(\Sigma_i\), and give the condition under which it equals \(R_i\).
Why does positive definiteness of \(R_i\) imply positive definiteness of \(D_iR_iD_i\) when all marginal standard deviations are positive?
For a binary response with mean \(0.7\), calculate its marginal variance and standard deviation.
Explain why the product of two observed standardized residuals is not, by itself, the working correlation for that pair.