Method-of-Moments Estimation for Exchangeable Correlation

From Residual Pairs to a Pooled Moment

An exchangeable working correlation uses one parameter, \(\rho\), for every distinct pair of responses in the same cluster. Estimating that parameter from observed data therefore amounts to pooling information across within-cluster pairs. The key choice is how to pool: give every pair equal weight, or first summarize each cluster and then give every cluster equal weight. These approaches agree when cluster sizes are equal, but need not agree otherwise.

Let cluster \(i\) contain \(m_i\) responses, with marginal mean \(\mu_{ij}\) and variance \(s_{ij}^2\), for \(i=1,\ldots,K\) and \(j=1,\ldots,m_i\). Define the true standardized residual \(U_{ij}=(Y_{ij}-\mu_{ij})/s_{ij}\). In practice, the means and variances are estimated, giving the standardized residual \(\hat r_{ij}=(Y_{ij}-\hat\mu_{ij})/\hat s_{ij}\). For a GLM, the estimated marginal variance may, for example, be written \(\hat s_{ij}^2=\hat\phi V(\hat\mu_{ij})\), where \(\hat\phi\) is a scale estimate and \(V\) is the variance function.

There are \(q_i=\binom{m_i}{2}=m_i(m_i-1)/2\) distinct pairs in cluster \(i\). Clusters of size one have \(q_i=0\), so they provide no within-cluster pair for estimating \(\rho\). Write \(\hat S_i=\sum_{j<k}\hat r_{ij}\hat r_{ik}\). The pooled estimators below use these cross-products; they are moment estimators based on fitted residuals, not exact finite-sample unbiased estimators in general.

Two Ways to Pool the Cross-Products

The pair-weighted estimator combines all observed pair products before dividing by the total number of pairs. Thus a cluster with many pairs contributes more than a cluster with few pairs. The cluster-weighted estimator first divides each cluster’s cross-product sum by its own pair count, then averages those cluster averages equally. Let \(\mathcal I=\{i:m_i\ge2\}\) and \(C=|\mathcal I|\).

The first formula is undefined if there are no pairs in the data; so is the second if \(C=0\). This happens, for example, if every cluster has size one. The pair-weighted formula is equivalently a weighted average of the cluster-specific averages \(\hat S_i/q_i\), with weights \(q_i\). This identity makes the weighting difference explicit.

To verify the result, let \(Q=\sum_{i\in\mathcal I}q_i>0\). For each within-cluster pair, \(E[U_{ij}U_{ik}]=\rho\). Linearity of expectation therefore gives \(E[S_i]=q_i\rho\), without requiring different pair products to be independent. Hence \[ E\left[\frac{\sum_i S_i}{Q}\right]=\frac{\sum_i q_i\rho}{Q}=\rho. \] For the cluster-weighted version, \(E[S_i/q_i]=\rho\) for every \(i\in\mathcal I\), so \[ E\left[\frac{1}{C}\sum_{i\in\mathcal I}\frac{S_i}{q_i}\right] =\frac{1}{C}\sum_{i\in\mathcal I}\rho=\rho. \] These calculations establish the population-moment claim. They do not establish that either estimator has smaller variance, or that either plug-in estimate is exactly unbiased after fitting the mean model. \(\square\)

Admissible Values From the Correlation Matrix

A moment estimate is not automatically a valid exchangeable correlation. For a cluster of size \(m\ge2\), its exchangeable correlation matrix has ones on the diagonal and \(\rho\) off the diagonal. It can be written \(R_m(\rho)=(1-\rho)I_m+\rho J_m\), where \(I_m\) is the identity matrix and \(J_m\) is the matrix of ones.

For a proof, the vector of all ones is an eigenvector of \(J_m\) with eigenvalue \(m\). Every vector orthogonal to it is an eigenvector of \(J_m\) with eigenvalue zero. Consequently, \(R_m(\rho)\) has eigenvalue \(1+(m-1)\rho\) in the all-ones direction and eigenvalue \(1-\rho\) in each of the \(m-1\) orthogonal directions. Positive definiteness requires both eigenvalues to be strictly positive: \[ 1+(m-1)\rho>0 \quad\text{and}\quad 1-\rho>0. \] Solving gives \(\rho>-1/(m-1)\) and \(\rho<1\). Replacing strict inequalities by weak ones gives the positive-semidefinite condition. \(\square\)

For several cluster sizes, a common working correlation must satisfy the bound for every cluster. If \(m_{\max}\) is the largest observed cluster size of at least two, the combined positive-definite range is \[ -\frac{1}{m_{\max}-1}<\rho<1. \] Indeed, \(-1/(m-1)\) increases as \(m\) increases, becoming less negative. Therefore the largest cluster size gives the greatest lower bound, which is the strictest one. Singleton clusters impose no correlation-matrix restriction because they contain no off-diagonal entries.

Interpretation, Practical Use, and a Common Pitfall

The two pooling rules answer slightly different practical questions. Pair weighting estimates the common moment by treating each within-cluster pair as one contribution. Cluster weighting treats each cluster with at least two observations as one contribution, regardless of how many pairs it contains. Under the common-correlation moment model, both target \(\rho\); with finite data, their observed values can differ, as the examples show. The chosen rule should be reported clearly, especially when cluster sizes vary.

In an actual GEE analysis, the residuals are computed using estimated means and possibly an estimated scale. The resulting cross-products are useful plug-in moments, but fitting can induce dependence among residuals and alter their finite-sample moments. Also, the pair products within a cluster are generally dependent because they share responses. The expectation calculations above rely only on linearity and the stated common pair moment; they do not imply independence of products, nor do they by themselves provide a standard error for \(\hat\rho\).

• Common mistake: averaging all cluster averages and calling it pair-weighted. That is the cluster-weighted estimator. Pair weighting divides the total cross-product sum by the total pair count.
• Common mistake: checking only that the estimate lies between \(-1\) and \(1\). That is not enough for an exchangeable matrix when cluster size exceeds two. For size \(m\), the lower bound is \(-1/(m-1)\).
• Correction: Find the largest cluster size, derive its lower bound, and ensure the working value is strictly inside the common admissible interval when positive definiteness is required.