From Cluster Responses to a Robust Covariance
A GEE analysis brings together two distinct calculations: the mean-model score and the dependence structure used to weight that score. The working correlation affects the weighting, while the sandwich covariance uses the observed cluster scores to quantify their variability. A numerical example makes both roles visible—and helps prevent a common error when estimating correlation from residuals.
Suppose there are \(K\) independent clusters. Cluster \(i\) has response vector \(Y_i\), mean vector \(\mu_i(\beta)\), and derivative matrix \(D_i=\partial\mu_i/\partial\beta^\mathsf{T}\). In a Gaussian identity-link model with common marginal variance \(\phi\) and exchangeable working correlation \(R_i(\rho)\), write \(V_i=\phi R_i(\rho)\). The GEE and its cluster contribution are
After solving the GEE, let \(\hat U_i\) be each fitted cluster score. The empirical sandwich covariance has bread \(H=\sum_i D_i^\mathsf{T}V_i^{-1}D_i\) and meat \(M=\sum_i\hat U_i\hat U_i^\mathsf{T}\):
The cluster, not the individual response, is the unit whose score is squared and summed in \(M\). This calculation is therefore different from treating all measurements as independent when estimating coefficient uncertainty.
A Useful Two-Measurement Simplification
For a cluster with two measurements, the exchangeable working correlation matrix and its inverse are
The condition \(|\rho|<1\) ensures the displayed inverse exists. If the covariates are constant within each cluster, write the cluster’s design matrix as \(D_i=\mathbf{1}_2 x_i^\mathsf{T}\), where \(x_i\) is the covariate vector and \(\mathbf{1}_2=(1,1)^\mathsf{T}\). Multiplication gives \(R(\rho)^{-1}\mathbf{1}_2=\mathbf{1}_2/(1+\rho)\). Thus the cluster score is proportional to the sum of its two residuals.
When all clusters have two observations and the same \(\phi\) and \(\rho\), the factor \(1/\{\phi(1+\rho)\}\) is common to every score. It does not change the root of the GEE. This special simplification need not hold with covariates that vary within a cluster or with different cluster sizes.
Worked Example: Fixed Versus Estimated Marginal Variance
Consider four residual pairs \((1,1),(1,0),(0,1),(0,0)\). This small calculation illustrates why the denominator in a residual-moment estimate depends on whether the marginal variance is fixed or estimated.
The sum of the four cross-products is \(1+0+0+0=1\). The total sum of the eight squared residuals is \(2+1+1+0=4\).
If the marginal variance is fixed at one, the exchangeable moment equation is \(E[e_{i1}e_{i2}]=\rho\). The sample average cross-product therefore gives \(\hat\rho=1/4\).
If the common marginal variance is estimated by its average residual square, the estimated variance moment is \(4/(2\cdot4)=1/2\). The average cross-product is \(1/4\), so the ratio estimate is \(\hat\rho=(1/4)/(1/2)=1/2\). These are different calculations because the variance assumption differs.
Worked Gaussian GEE
There are four independent clusters, each with two responses. Let \(x_i=(1,z_i)^\mathsf{T}\), where \(z_i=0\) denotes control and \(z_i=1\) denotes treatment. The model is \(E[Y_{ij}\mid z_i]=\beta_0+\beta_1z_i\). The observed vectors are:
| Cluster | Group | \(Y_i\) |
|---|---|---|
| 1 | Control | \((5,4)^\mathsf{T}\) |
| 2 | Control | \((3,4)^\mathsf{T}\) |
| 3 | Treatment | \((7,7)^\mathsf{T}\) |
| 4 | Treatment | \((5,5)^\mathsf{T}\) |
The group averages are \(4\) for control and \(6\) for treatment. The fitted coefficients are therefore \(\hat\beta_0=4\) and \(\hat\beta_1=6-4=2\). In this balanced cluster-level design, the common score factor just derived does not change those estimates.
Worked Example: Estimating the Working Correlation
At \(\hat\beta=(4,2)^\mathsf{T}\), the fitted residual vectors are \(e_1=(1,0)^\mathsf{T}\), \(e_2=(-1,0)^\mathsf{T}\), \(e_3=(1,1)^\mathsf{T}\), and \(e_4=(-1,-1)^\mathsf{T}\). Their cross-products sum to \(0+0+1+1=2\), and their squared residuals sum to \(1+1+2+2=6\).
Use a ratio of residual moments: the average pair cross-product is \(2/4=1/2\), while the average marginal residual square is \(6/(2\cdot4)=3/4\). Thus
For the scale in the working covariance, use the residual sum of squares divided by the residual degrees of freedom: \(\hat\phi=6/(8-2)=1\). This degrees-of-freedom adjustment is a scale convention; it is not a reason to replace the moment denominator in the correlation ratio. Here the estimated correlation is strictly between \(-1\) and \(1\), so the two-by-two working correlation inverse exists.
Aggregating Cluster Scores and Calculating the Sandwich
For these fitted residuals, the cluster residual sums are \(1,-1,2,-2\). Let \(q=1/\{\hat\phi(1+\hat\rho)\}=3/5\). The two control clusters have covariate vector \((1,0)^\mathsf{T}\); the two treatment clusters have vector \((1,1)^\mathsf{T}\). Applying the derived score formula gives
The bread follows from \(H=\sum_i D_i^\mathsf{T}V_i^{-1}D_i\). Each cluster contributes \(2q x_ix_i^\mathsf{T}\), since \(D_i=\mathbf{1}_2x_i^\mathsf{T}\). Summing two control and two treatment clusters gives
For the meat, each score contributes its outer product. The two control scores together contribute \(2q^2\begin{pmatrix}1&0\\0&0\end{pmatrix}\); the two treatment scores together contribute \(8q^2\begin{pmatrix}1&1\\1&1\end{pmatrix}\). Hence
Worked Example: Standard Error for the Treatment Effect
The treatment effect is the second coefficient, so its contrast vector is \(c=(0,1)^\mathsf{T}\). Its robust estimated variance is \(c^\mathsf{T}\widehat{\operatorname{Cov}}_{\mathrm{rob}}(\hat\beta)c=5/8\), giving standard error \(\sqrt{5/8}\approx0.791\).
For comparison, the model-based covariance is \(H^{-1}=(1/q)C=(5/3)C\); its treatment-effect variance is \(5/6\), with standard error \(\sqrt{5/6}\approx0.913\). The two estimates need not agree: the robust calculation uses the observed cluster-score outer products, whereas the model-based calculation uses the working covariance directly.
Common Mistake and Interpretation
A frequent error is to say, “the marginal variance is fixed at one,” and then divide the average residual cross-product by the empirical average residual square. That division estimates a ratio of moments when the variance is also estimated; it does not follow from the fixed-variance moment equation. With fixed unit variance, the moment estimate is simply the average cross-product. The distinction is explicit in the four-pair example: the fixed-variance estimate is \(1/4\), while the ratio estimate is \(1/2\).
The worked analysis also illustrates why cluster scores must remain intact when forming the sandwich meat: each cluster contributes one outer product, including dependence among responses in that cluster. Finally, four clusters are far too few to treat the displayed robust standard error as dependable large-sample inference. The calculation demonstrates the mechanics; it does not overcome the limited information in a small number of independent clusters.
Check Your Understanding
Use the definitions and calculations above to check that you can distinguish the mean fit, working correlation, and robust uncertainty calculation.
- For two observations per cluster and cluster-constant covariates, why does the exchangeable GEE score depend on the sum of the two residuals?
- If the marginal variance is fixed at one, what sample moment estimates \(\rho\) from residual pairs?
- In the four-cluster analysis, verify the treatment and control fitted means from the response vectors.
- Why is the sandwich meat formed from \(\hat U_i\hat U_i^\mathsf{T}\) for each cluster rather than from individual residual squares alone?
- Which entry of the robust covariance matrix gives the estimated variance of the treatment coefficient, and what is its standard error?