Pairwise Residual Moments for Correlation Estimation

Why Pairwise Residual Products Estimate Correlation

In a generalized estimating equation analysis, the mean model describes the expected response at each measurement, while a working correlation describes dependence between responses from the same cluster. Once the marginal means and variances have been fitted, each response can be put on a standardized residual scale. The product of two such residuals is informative about their correlation: positive products indicate that the two responses tend to lie on the same side of their marginal means, while negative products indicate opposite deviations.

The method-of-moments idea is to match an average of observed residual products to the corresponding model correlation. The important practical detail is deciding which pairs belong in that average. A pair type may be specified by visit labels, time separation, or another feature of the design. Clusters can contribute different numbers of eligible pairs, and fitted marginal means enter every product. These choices determine the numerical estimate.

Standardized Residuals and Pair Types

Let $Y_{ij}$ be the response for measurement $j$ in cluster $i$, where $i=1,\ldots,n$ and cluster $i$ has $m_i$ observed measurements. Write $\mu_{ij}=E[Y_{ij}]$ and $\operatorname{Var}(Y_{ij})=\phi v(\mu_{ij})$, where $v$ is the marginal variance function and $\phi$ is a dispersion parameter when one is needed. The standardized residual based on the true marginal quantities is

Z_{ij}=\frac{Y_{ij}-\mu_{ij}}{\sqrt{\phi v(\mu_{ij})}}.

In practice, the marginal mean and variance are estimated. If $\hat\mu_{ij}$ and $\hat\phi$ are the fitted quantities, the plug-in standardized residual is

\hat Z_{ij}=\frac{Y_{ij}-\hat\mu_{ij}}{\sqrt{\hat\phi v(\hat\mu_{ij})}}.

For a Gaussian model with known common variance $\sigma^2$, this reduces to $(Y_{ij}-\hat\mu_{ij})/\sigma$. The denominator matters: raw residual products are on the response scale, whereas standardized products target a correlation scale.

Definition: A pair type is a specified class of within-cluster measurement pairs that share one correlation parameter under the working correlation model. For a pair type $g$, let $\mathcal P_g$ be the set of eligible triples $(i,j,k)$, with $j<k$, whose two measurements form a pair of type $g$. The number of eligible pairs is $N_g=|\mathcal P_g|$.

The type definition must be fixed by the design or the working correlation structure before the products are pooled. For example, a type might consist of pairs at a specified visit separation. This tutorial uses pair types generally; it does not assume one particular correlation structure.

The Moment Identity

The connection between a correlation and a standardized residual product follows directly from the definition of covariance. It requires finite, positive marginal variances and correctly specified marginal means and variances for the responses in question.

Theorem: Suppose $E[Y_{ij}]=\mu_{ij}$, $0<\operatorname{Var}(Y_{ij})=\phi v(\mu_{ij})<\infty$, and the analogous conditions hold for $Y_{ik}$. Then $$ E[Z_{ij}Z_{ik}]=\operatorname{Corr}(Y_{ij},Y_{ik}). $$ Consequently, if every pair in $\mathcal P_g$ has the same correlation $\rho_g$, then the average true standardized residual product over those pairs has expectation $\rho_g$, provided $\mathcal P_g$ is fixed and nonempty.

Proof. By definition, $E[Z_{ij}]=E[Z_{ik}]=0$. The denominators in $Z_{ij}$ and $Z_{ik}$ are nonrandom marginal standard deviations, so

\begin{aligned} E[Z_{ij}Z_{ik}] &=\frac{E[(Y_{ij}-\mu_{ij})(Y_{ik}-\mu_{ik})]} {\sqrt{\operatorname{Var}(Y_{ij})\operatorname{Var}(Y_{ik})}}\\ &=\frac{\operatorname{Cov}(Y_{ij},Y_{ik})} {\sqrt{\operatorname{Var}(Y_{ij})\operatorname{Var}(Y_{ik})}}\\ &=\operatorname{Corr}(Y_{ij},Y_{ik}). \end{aligned}

If all $N_g$ pairs have correlation $\rho_g$, linearity of expectation gives

E\left[\frac{1}{N_g}\sum_{(i,j,k)\in\mathcal P_g}Z_{ij}Z_{ik}\right] =\frac{1}{N_g}\sum_{(i,j,k)\in\mathcal P_g}\rho_g =\rho_g.

No independence between pairs is needed for this expectation calculation. $\square$

This identity motivates the plug-in method-of-moments estimate

\hat\rho_g^{\,\mathrm{mom}} =\frac{1}{N_g}\sum_{(i,j,k)\in\mathcal P_g}\hat Z_{ij}\hat Z_{ik}, \qquad N_g>0.

It is a pair-weighted average: every eligible pair contributes one product. This formula is a moment estimate, not a claim that the fitted residual products have exactly the same expectation as products formed with the true marginal quantities.

Worked Example: Pooling Products for One Pair Type

Suppose three independent clusters each contribute one eligible pair of the same type. Their true standardized residuals, for the two measurements in each pair, are observed as $(1,-0.5)$, $(0.5,1)$, and $(-1,-0.5)$. Calculate the pairwise moment.

The three products are

1(-0.5)=-0.5,\qquad 0.5(1)=0.5,\qquad (-1)(-0.5)=0.5.

There are $N_g=3$ eligible pairs. Their pooled average is therefore

\hat\rho_g^{\,\mathrm{mom}} =\frac{-0.5+0.5+0.5}{3} =\frac{0.5}{3} =\frac{1}{6}\approx 0.167.

This is the moment estimate from these observed products. It is not necessary for each individual product to resemble a correlation or to lie between $-1$ and $1$; the moment estimate itself can also fall outside that interval in a finite sample.

Unequal Numbers of Eligible Pairs

Unequal cluster sizes affect the calculation through the eligible-pair sets. A cluster with more measurements may contribute more pairs of a given type. Under the pair-weighted formula, its contribution is the sum of its eligible products divided by the total number of eligible pairs across all clusters. It is not automatically given the same total weight as a smaller cluster.

To make this explicit, define $q_{ig}$ as the number of type-$g$ pairs in cluster $i$, and define $U_{ig}$ as their sum of standardized residual products. Then $N_g=\sum_i q_{ig}$, and the pooled estimate can be written

\hat\rho_g^{\,\mathrm{mom}} =\frac{\sum_{i=1}^{n}U_{ig}}{\sum_{i=1}^{n}q_{ig}} =\sum_{i:q_{ig}>0}\frac{q_{ig}}{N_g}\left(\frac{U_{ig}}{q_{ig}}\right).

Thus the cluster-specific average product is weighted by $q_{ig}/N_g$. An alternative calculation that averages cluster-specific averages equally answers a different weighting question, and is not generally the same estimator.

Worked Example: Pair Weighting Versus Equal Cluster Weighting

Two clusters contribute products for the same pair type. Cluster 1 has one eligible pair, with product $0.8$. Cluster 2 has two eligible pairs, with products $0.2$ and $-0.2$. Compare the pooled pair moment with the average of the two cluster averages.

For pooled pair weighting, there are three products in all:

\hat\rho_g^{\,\mathrm{mom}} =\frac{0.8+0.2-0.2}{3} =\frac{0.8}{3} \approx 0.267.

The cluster-specific averages are $0.8/1=0.8$ and $(0.2-0.2)/2=0$. Giving the two clusters equal weight instead produces

\frac{0.8+0}{2}=0.4.

The estimates differ because the second calculation assigns half the total weight to each cluster, even though the clusters supply different numbers of pairs. Neither weighting should be substituted for the other without stating the intended target and calculation.

Theorem: Suppose clusters are independent, the eligible pair sets are fixed, and every type-$g$ pair has correlation $\rho_g$ and correctly standardized marginal residuals. Assume $q_{ig}\le C$ for a fixed finite constant $C$, $\sum_i q_{ig}\ge cn$ for some $c>0$, and $\operatorname{Var}(U_{ig})\le K$ for a fixed finite $K$. Then the average using true residuals converges in probability to $\rho_g$ as $n$ grows.

Proof. The moment identity gives $E[U_{ig}]=q_{ig}\rho_g$. Hence

E\left[\frac{\sum_i U_{ig}}{\sum_i q_{ig}}\right]=\rho_g.

Since clusters are independent, $\operatorname{Var}(\sum_i U_{ig})=\sum_i\operatorname{Var}(U_{ig})\le nK$. The denominator is fixed and at least $cn$, so Chebyshev’s inequality gives, for every $\epsilon>0$,

P\left(\left|\frac{\sum_i U_{ig}}{\sum_i q_{ig}}-\rho_g\right|>\epsilon\right) \le \frac{nK}{\epsilon^2 c^2 n^2} =\frac{K}{\epsilon^2 c^2 n} \longrightarrow 0.

This proves convergence in probability for the average formed with true standardized residuals. $\square$

The result does not itself establish an exact finite-sample property for plug-in residuals. In applications, fitted means and variances replace the true ones; their estimation must be accounted for when assessing finite-sample behavior or establishing a result for a particular fitting procedure.

How Fitted Marginal Means Enter

A fitted mean is measurement-specific whenever the marginal model says the mean varies with covariates or visit. Therefore, each residual in a pair must use the fitted mean and variance for that particular measurement. Using one overall mean for all observations can change both residual signs and the resulting product. The following calculation shows the effect directly.

Worked Example: Residual Products After Fitting Visit Means

Consider two clusters measured at two visits. A Gaussian marginal model has variance $4$ at each visit and a separate mean for each visit. The observed responses are $(5,9)$ in cluster 1 and $(3,7)$ in cluster 2. Calculate the moment for the pair type consisting of the two visits.

The fitted mean at visit 1 is the average $(5+3)/2=4$, and at visit 2 it is $(9+7)/2=8$. The standard deviation is $\sqrt{4}=2$. Thus the standardized residuals are

\begin{array}{c|cc|c} &\hat Z_{i1}&\hat Z_{i2}&\hat Z_{i1}\hat Z_{i2}\\ \hline \text{Cluster 1}&(5-4)/2=0.5&(9-8)/2=0.5&0.25\\ \text{Cluster 2}&(3-4)/2=-0.5&(7-8)/2=-0.5&0.25 \end{array}

There are two eligible pairs, so the pairwise moment is $(0.25+0.25)/2=0.25$. If one instead incorrectly used the overall mean $6$ for both visits, the products would be $(-0.5)(1.5)=-0.75$ and $(-1.5)(0.5)=-0.75$, giving $-0.75$. That is a different residual calculation because it ignores the fitted visit-specific marginal means.

Interpretation and Common Pitfall

A pairwise residual moment is best understood as an empirical match to a model identity. When the marginal means and variances are known, the identity is exact in expectation. In a fitted GEE, the estimated means and possibly estimated dispersion are substituted, so the resulting products are dependent on the data used to fit the mean model. This can affect small-sample behavior. The calculation still provides a transparent way to summarize the dependence information in residual pairs, but an estimate should not be described as exactly unbiased merely because the population identity holds.

Missing measurements also change the eligible-pair set. A pair contributes only when both responses are observed and their pair type is defined. Consequently, both the numerator and $N_g$ must be calculated from the same set of observed pairs. If the set of observed pairs is random, the simple fixed-set expectation and proof above do not automatically cover the resulting estimator.

Common pitfall: Averaging cluster-level pair averages and calling the result the pooled pair moment. This gives each contributing cluster equal weight, even when one cluster supplies several times as many eligible pairs as another. For the pair-weighted moment, sum every eligible product and divide by the total number of eligible pairs. If equal cluster weighting is intended, state that it is a different summary.

Before interpreting a reported moment estimate, check three items: the pair-type definition, the number of contributing pairs from each cluster, and the fitted marginal mean and variance used for every residual. Those details determine what the number summarizes.

Check Your Understanding

Use the definitions and calculations above to answer these questions.

Why does the expectation of a product of two true standardized residuals equal the correlation of the corresponding responses?
A pair type has three eligible products, $0.4$, $-0.1$, and $0.3$. What is the pair-weighted moment estimate?
Two clusters contribute one and three eligible pairs, respectively. How does the pooled pair moment weight their cluster-specific average products?
Why can replacing visit-specific fitted means with one overall mean substantially change the estimated moment?
Does the moment identity imply that a plug-in estimate using fitted means is exactly unbiased in a finite sample? Explain.