Marginal Contrasts with Robust Covariance

Why Contrasts Need the Robust Covariance

A fitted GEE model usually reports estimates for individual regression coefficients, but scientific questions often concern combinations of them. For example, a treatment effect may depend on visit, so the treatment effect at a particular visit is the sum of a treatment coefficient and a treatment-by-visit interaction coefficient. To quantify uncertainty in that sum, it is not enough to look at the two coefficient standard errors separately: their covariance also matters.

GEE coefficients describe the specified marginal mean model. In particular, for a response $Y_{ij}$ in cluster $i$, the model specifies its mean given the covariates, $\mu_{ij}=E[Y_{ij}\mid X_{ij}]=g^{-1}(X_{ij}^{T}\beta)$. The link $g$ connects that population-average mean to the linear predictor. This tutorial focuses on contrasts of $\beta$, their cluster-robust uncertainty, and the careful translation of those contrasts to the response scale.

Definition: A linear contrast of a coefficient vector $\beta\in\mathbb{R}^{p}$ is $c^{T}\beta$, where $c\in\mathbb{R}^{p}$ is a fixed, specified vector. Its GEE estimate is $c^{T}\hat{\beta}$. A contrast is a scalar quantity; a contrast matrix $C$ can specify several such quantities at once, giving $C\beta$.

Let $\widehat{\operatorname{Cov}}_{\mathrm{rob}}(\hat{\beta})$ denote the cluster-robust sandwich covariance estimate. It accounts for the estimated variability across independent clusters and does not require the working within-cluster correlation to be correct for robust large-sample inference, provided the marginal mean model and the relevant regularity conditions are appropriate. A short way to remember the calculation is that the entire covariance matrix, including its off-diagonal entries, enters a contrast's variance.

Result: The estimated robust variance and standard error of a fixed linear contrast are $$ \widehat{\operatorname{Var}}_{\mathrm{rob}}(c^{T}\hat{\beta}) =c^{T}\widehat{\operatorname{Cov}}_{\mathrm{rob}}(\hat{\beta})c, \qquad \widehat{\operatorname{SE}}_{\mathrm{rob}}(c^{T}\hat{\beta}) =\sqrt{c^{T}\widehat{\operatorname{Cov}}_{\mathrm{rob}}(\hat{\beta})c}. $$

Deriving the Contrast Variance

This formula follows directly from the covariance rule for linear transformations. It is useful to see why that rule requires the full covariance matrix rather than only the diagonal entries.

Derivation: For any random vector $Z$ with finite covariance matrix $\Sigma$ and any fixed vector $c$, $$ \operatorname{Var}(c^{T}Z) =E\left[\left\{\sum_{r=1}^{p}c_r(Z_r-E[Z_r])\right\}^{2}\right] =\sum_{r=1}^{p}\sum_{s=1}^{p}c_rc_s\operatorname{Cov}(Z_r,Z_s) =c^{T}\Sigma c. $$ Apply this identity to the estimated coefficient vector and estimate its covariance with the robust sandwich matrix. The diagonal terms give each coefficient's variance contribution; the off-diagonal terms give the covariance contributions. In practice the sandwich covariance is estimated, so this calculation gives an estimated variance, not an exact finite-sample variance. $\square$

If there are $m$ independent clusters and the fitted mean model and regularity conditions support the usual GEE asymptotics, then $\hat{\beta}$ is approximately multivariate normal for sufficiently large $m$. The robust Wald interval for a scalar contrast is consequently $$ c^{T}\hat{\beta}\ \mathbin{\pm}\ z_{1-\alpha/2} \sqrt{c^{T}\widehat{\operatorname{Cov}}_{\mathrm{rob}}(\hat{\beta})c}, $$ where $z_{1-\alpha/2}$ is the standard normal quantile. This is a large-sample approximation, not a guarantee of accurate inference with few clusters. The number of independent clusters, not merely the total number of repeated observations, is central to its reliability.

Worked Example: A Treatment Effect at One Visit

Suppose a GEE linear predictor includes treatment $T$, a visit indicator $V$, and their interaction: $\eta=\beta_0+\beta_TT+\beta_VV+\beta_{TV}TV$. Let $T=1$ denote treatment and $V=1$ the visit of interest. Assume the fitted estimates are $\hat{\beta}_T=0.30$ and $\hat{\beta}_{TV}=-0.10$, with robust variances $0.0400$ and $0.0225$, respectively, and robust covariance $-0.0150$.

At $V=1$, the treatment-versus-control difference in the linear predictor is $\beta_T+\beta_{TV}$. Thus $c=(1,1)^T$ for the coefficient subvector $(\beta_T,\beta_{TV})^T$, and the estimated contrast is $0.30-0.10=0.20$.

Its robust variance includes both covariance terms: $0.0400+0.0225+2(-0.0150)=0.0325$. Therefore its robust standard error is $\sqrt{0.0325}=0.1803$. A 95% large-sample Wald interval on the linear-predictor scale is $0.20\pm1.96(0.1803)$, or approximately $(-0.153,0.553)$.

If this is a log-link model, exponentiating the contrast gives a population-average mean ratio of $\exp(0.20)=1.221$: at this visit, the fitted treatment mean is about 1.22 times the fitted control mean, under the model. Exponentiating the interval endpoints gives approximately $(0.858,1.739)$. The negative covariance reduced the variance here; ignoring it would give $0.0625$, not $0.0325$.

From a Coefficient Contrast to a Response-Scale Effect

A contrast $c^{T}\beta$ is a difference on the linear-predictor scale. Its response-scale interpretation depends on the link. Under an identity link, a difference in linear predictors is a difference in means. Under a log link, it is a log ratio of means, so exponentiating gives a mean ratio. Under a logit link, it is a log odds ratio, so exponentiating gives an odds ratio. In the latter two cases, the coefficient contrast is not itself a difference in expected responses.

For a nonlinear inverse link, a response-scale difference at two covariate patterns $x_A$ and $x_B$ is $h(\beta)=g^{-1}(x_A^{T}\beta)-g^{-1}(x_B^{T}\beta)$. This is generally a nonlinear function of the coefficients, even though the difference of linear predictors, $(x_A-x_B)^T\beta$, is a linear contrast. Its uncertainty can be approximated by the delta method.

Result: Suppose $\hat{\beta}\xrightarrow{p}\beta$, $\sqrt{m}(\hat{\beta}-\beta)\xrightarrow{d}N(0,\Omega)$, and $h$ is differentiable in a neighborhood of $\beta$. If a consistent covariance estimate for $\hat{\beta}$ is available, then $$ \widehat{\operatorname{Var}}\{h(\hat{\beta})\} \ \approx\ \nabla h(\hat{\beta})^{T} \widehat{\operatorname{Cov}}_{\mathrm{rob}}(\hat{\beta}) \nabla h(\hat{\beta}). $$

For the two-pattern response difference, write $\mu_A=g^{-1}(x_A^{T}\beta)$ and $\mu_B=g^{-1}(x_B^{T}\beta)$. By the chain rule its gradient is $$ \nabla h(\beta) =\left.\frac{d\,g^{-1}(\eta)}{d\eta}\right|_{\eta=x_A^{T}\beta}x_A -\left.\frac{d\,g^{-1}(\eta)}{d\eta}\right|_{\eta=x_B^{T}\beta}x_B. $$ This gradient, evaluated at $\hat{\beta}$, is used in the sandwich quadratic form above.

Justification of the delta-method variance: Differentiability at $\beta$ gives the first-order expansion $$ h(\hat{\beta})-h(\beta) =\nabla h(\beta)^{T}(\hat{\beta}-\beta)+r_m, \qquad \frac{|r_m|}{\|\hat{\beta}-\beta\|}\xrightarrow{p}0. $$ The assumed asymptotic normality implies $\hat{\beta}-\beta=O_p(m^{-1/2})$, so $\sqrt{m}r_m\xrightarrow{p}0$. Therefore the first-order term determines the limiting distribution of $\sqrt{m}\{h(\hat{\beta})-h(\beta)\}$; its asymptotic variance is $\nabla h(\beta)^T\Omega\nabla h(\beta)$. Replacing the gradient and covariance by consistent estimates yields the displayed plug-in variance estimate. This reasoning requires differentiability and the stated large-sample conditions; it does not make the approximation exact in a small sample. $\square$

Worked Example: A Log-Link Contrast as a Mean Ratio

Suppose a Poisson GEE uses a log link and includes treatment $T$ and a treatment-by-period indicator $V$. At $V=1$, the treatment contrast is $L=\beta_T+\beta_{TV}$. Let $\hat{\beta}_T=0.25$, $\hat{\beta}_{TV}=0.15$, and suppose the robust covariance submatrix is $\begin{pmatrix}0.0100&0.0020\\0.0020&0.0064\end{pmatrix}$.

The estimated log mean ratio is $\hat L=0.25+0.15=0.40$. With $c=(1,1)^T$, its robust variance is $0.0100+0.0064+2(0.0020)=0.0204$, and its standard error is $\sqrt{0.0204}=0.1428$. The 95% interval for $L$ is $0.40\pm1.96(0.1428)=(0.120,0.680)$, to three decimals.

The response-scale mean ratio is $\exp(0.40)=1.492$, with transformed interval $(\exp(0.120),\exp(0.680))=(1.128,1.974)$. Thus the fitted population-average mean at the specified period is about 1.49 times as large under treatment as under control. This multiplicative interpretation follows from the log link; it is not a claim that the mean increases by 0.40 response units.

Worked Example: A Logistic GEE Probability Difference

Consider a logistic GEE with $\eta=\beta_0+\beta_TT$, where $T=0$ is control and $T=1$ is treatment. Suppose $\hat{\beta}_0=-1.00$, $\hat{\beta}_T=0.80$, and the robust covariance matrix is $\begin{pmatrix}0.09&-0.015\\-0.015&0.04\end{pmatrix}$. The coefficient contrast for treatment versus control is $\beta_T$, with standard error $\sqrt{0.04}=0.20$. It is a log odds ratio; the fitted odds ratio is $\exp(0.80)=2.23$.

For the response-scale probability difference, first calculate the two fitted probabilities. The control linear predictor is $-1.00$, giving $\hat p_0=1/(1+\exp(1))=0.2689$. The treatment linear predictor is $-1.00+0.80=-0.20$, giving $\hat p_1=1/(1+\exp(0.20))=0.4502$. Thus $\hat p_1-\hat p_0=0.1812$, or about 18.1 percentage points.

For $h(\beta)=p_1-p_0$, the gradient with respect to $(\beta_0,\beta_T)^T$ is $(p_1(1-p_1)-p_0(1-p_0),\,p_1(1-p_1))^T$. At the fitted values this is approximately $(0.0509,0.2475)^T$. The delta-method variance is $0.09(0.0509)^2+0.04(0.2475)^2+2(-0.015)(0.0509)(0.2475)=0.00231$, so the standard error is approximately $0.0480$. A rough 95% Wald interval for the probability difference is $0.1812\pm1.96(0.0480)=(0.087,0.275)$. The negative covariance contributes to this uncertainty calculation; it is not optional even though the final estimand is a probability difference.

What “Population-Average” Means Here

In GEE, the fitted mean $g^{-1}(x^T\hat{\beta})$ is a marginal mean for responses at covariate pattern $x$, rather than a mean conditional on an unobserved cluster-specific random effect. That is why the log-link ratio and logistic-model probabilities in the examples are described as population-average model quantities. They are still conditional on the stated observed covariates.

If the scientific target is an effect averaged over a specified covariate distribution, a prediction at one pattern is not enough. For example, with prespecified weights $w_k$ summing to one over patterns $x_k$, a standardized mean difference is $\sum_k w_k\{g^{-1}(x_{k,1}^{T}\beta)-g^{-1}(x_{k,0}^{T}\beta)\}$. Its gradient is the same weighted sum of the pattern-specific gradients, and its robust delta-method variance uses that gradient and the full robust covariance matrix. State the target patterns or weights so the reported response-scale effect has an unambiguous population interpretation.

Common Pitfalls

Dropping covariance terms. The incorrect calculation $\operatorname{Var}(\hat{\beta}_1+\hat{\beta}_2)=\operatorname{Var}(\hat{\beta}_1)+\operatorname{Var}(\hat{\beta}_2)$ omits $2\operatorname{Cov}(\hat{\beta}_1,\hat{\beta}_2)$. Use the full robust covariance matrix.
Using the wrong contrast for an interaction model. A treatment coefficient alone is the treatment effect at the reference value of an interacting variable. At another value, form the appropriate sum of coefficients.
Calling a logit contrast a probability difference. A treatment coefficient in a logistic GEE is a log odds ratio. Calculate fitted probabilities and, if needed, their difference separately.
Treating a robust standard error as a small-sample guarantee. Sandwich inference is asymptotic. With few independent clusters, its normal approximation can be poor.
Leaving the target population implicit. A response-scale difference at one covariate pattern and a difference averaged over a covariate distribution are different estimands.

Key takeaway: Define the contrast from the scientific question, calculate its variance as $c^T\widehat{\operatorname{Cov}}_{\mathrm{rob}}(\hat{\beta})c$, and interpret it through the link. For nonlinear response-scale effects, calculate the predictions and use the gradient with the same robust covariance matrix.

Check Your Understanding

For each question, identify the relevant contrast or response-scale quantity and the uncertainty calculation required.

A model has treatment and treatment-by-visit coefficients. What contrast represents the treatment effect at a nonreference visit?
Why does the robust variance of a sum of two coefficients include twice their covariance?
Under a log link, what response-scale quantity corresponds to a difference of two linear predictors?
In a logistic GEE, why is exponentiating a treatment coefficient not the same as calculating a treatment probability difference?
What additional specification is needed to interpret a response-scale effect standardized across covariate patterns?