Why Use a Marginal Logit Model?
When a binary outcome is measured repeatedly on a person, or on several people within a site, outcomes from the same cluster may be dependent. A logistic regression mean model can describe the probability of the outcome, but treating every observation as independent does not account for that dependence in its uncertainty. Generalized estimating equations (GEE) address this by specifying a marginal mean model and using a working correlation to account for within-cluster dependence when estimating its coefficients.
For a binary outcome, the mean is a probability. A logit link models that probability directly as a function of covariates. The resulting regression coefficients describe changes in the population-averaged log odds, rather than effects conditional on a cluster-specific random effect.
The word marginal distinguishes this model from a logistic mixed model, whose mean is specified conditional on random effects. In a GEE model, the probability in the equation is the population-level mean for observations with the stated covariates, without conditioning on a latent cluster effect. The association among observations in a cluster is handled separately through a working correlation.
The Binary GEE and Its Working Correlation
Assume the clusters are independent of one another, outcomes are observed as specified by the analysis, and the marginal mean model is correctly specified. For each cluster, collect its outcomes into \(Y_i\) and its means into \(\mu_i\). Since a Bernoulli variable with mean \(\mu_{ij}\) has variance \(\mu_{ij}(1-\mu_{ij})\), define
The matrix \(R_i(\alpha)\) is a chosen working correlation matrix, indexed by nuisance parameter(s) \(\alpha\). It might use independence, exchangeable, AR(1), or unstructured structure, depending on the design and plausible within-cluster pattern. Its role is to describe a useful working covariance for the estimating procedure. It need not be the true correlation matrix for the regression estimates to be consistent when the mean model is correct and the regularity conditions for GEE hold.
Let \(D_i=\partial\mu_i/\partial\beta^T\), the derivative matrix of the mean vector with respect to the coefficients. The GEE estimator solves the following system:
For the logit link, the derivative of \(\mu_{ij}\) with respect to its linear predictor is \(\mu_{ij}(1-\mu_{ij})\). Thus each row of \(D_i\) is \(\mu_{ij}(1-\mu_{ij})X_{ij}^{T}\). Estimation generally proceeds iteratively because the means, variances, and possibly the working correlation depend on unknown parameters.
The working correlation affects how observations contribute to the coefficient estimates and their efficiency. It is not a license to ignore the binary outcome's mean-variance relationship: the diagonal variance still follows the Bernoulli model. Nor does choosing a correlation structure turn the GEE into a fully specified joint probability model for the binary response vector.
Worked Example: Choosing a Working Correlation for Repeated Outcomes
Suppose a hypothetical study records a binary symptom at baseline, month 1, and month 2 for each participant. A marginal logit model includes treatment and visit indicators. The research team expects measurements close in time to be more alike, and the visits occur at equal intervals.
An AR(1) working correlation is a reasonable starting choice because it represents correlation that declines with the number of time steps separating measurements. If the working lag-one correlation were \(\alpha=0.50\), its three-visit matrix would be
The lag-two entry is \(0.50^2=0.25\). This matrix is a working description, not a claim that the actual within-participant correlations equal these values. If the visits were not equally spaced, or if outcomes at every pair of visits seemed similarly associated, another structure could be more defensible. The choice should also be estimable from the available clusters.
For illustration, if the fitted marginal probabilities at the three visits were \(0.20\), \(0.30\), and \(0.40\), the Bernoulli variances would be \(0.20(0.80)=0.16\), \(0.30(0.70)=0.21\), and \(0.40(0.60)=0.24\). The working covariance is \(V_i=A_i^{1/2}R_iA_i^{1/2}\), so its baseline-to-month-1 covariance entry would be \(\sqrt{0.16}\,(0.50)\sqrt{0.21}\approx0.09165\). This calculation combines the marginal variances with the assumed working correlation; it does not change the fitted marginal means.
Interpreting Population-Averaged Odds Ratios
For a binary covariate \(x\) with coefficient \(\beta_x\), holding the other covariates fixed, a one-unit increase in \(x\) changes the modeled marginal log odds by \(\beta_x\). Exponentiating gives the corresponding population-averaged odds ratio, \(\exp(\beta_x)\). For a treatment indicator, it compares the marginal odds of the outcome for treatment with the marginal odds for the reference group at the same values of the other modeled covariates.
This is an odds ratio, not a risk ratio or a difference in probabilities. Its interpretation also depends on the covariate specification: with a treatment-by-time interaction, for example, the treatment odds ratio varies by time. A marginal odds ratio should not be called a conditional odds ratio given a random effect. Because odds ratios are non-collapsible, marginal and conditional odds ratios can differ even when there is no confounding; the two model types target different quantities.
Worked Example: Interpreting a Treatment Coefficient
Consider a hypothetical GEE for a binary recovery outcome, with treatment coded \(1\) for the new program and \(0\) for usual care. Suppose the fitted treatment coefficient is \(\hat\beta_{\text{treat}}=0.60\), with no treatment interaction in the model.
The population-averaged odds ratio is
At the same values of the other covariates in the model, the estimated marginal odds of recovery are about 1.82 times as high under the program as under usual care. This does not mean that recovery probability is 1.82 times as high, and it does not establish that the program caused a change. The interpretation is about the modeled population-averaged odds contrast.
If the standard error were \(0.20\), a Wald 95% confidence interval on the log-odds scale would be \(0.60\pm1.96(0.20)=(0.208,0.992)\). Exponentiating its endpoints gives \((1.231,2.696)\). Under the model and large-sample approximation, this interval procedure has 95% repeated-sampling coverage for the population-averaged odds ratio; it is not a 95% probability statement about this particular interval.
Why the Estimating Equation Works
The GEE combines residuals \(Y_i-\mu_i\) across independent clusters. The derivative matrix \(D_i\) maps those residuals to information about changes in \(\beta\), while \(V_i^{-1}\) weights the residual pattern using the working covariance. This makes the equation a natural extension of a score equation for a correlated response vector.
A simple derivation is available when the working correlation is independence. Then \(V_i=A_i\), and the contribution of observation \(j\) to the equation is
because \(\partial\mu_{ij}/\partial\beta=X_{ij}\mu_{ij}(1-\mu_{ij})\) for a logit link. Summing gives the familiar logistic-regression score equation. With a non-independence working correlation, residuals within a cluster are instead weighted jointly through \(V_i^{-1}\). The mean remains the same marginal logistic model; it is the weighting that changes.
When the working correlation is not correct, model-based standard errors that treat it as correct may be unreliable. GEE commonly uses the sandwich variance estimator. In broad form, it is the inverse sensitivity of the estimating equation multiplied by the empirical variability of cluster-level estimating contributions and then by the inverse sensitivity again. Because clusters, not individual observations, provide the independent units for this calculation, a large number of independent clusters is important for reliable large-sample inference.
Worked Example: An Intercept-Only Independence GEE
For arithmetic illustration, suppose four hypothetical independent clusters each have two binary observations, with observed pairs \((1,1)\), \((1,0)\), \((0,0)\), and \((0,0)\). Fit an intercept-only marginal logit model using independence as the working correlation. There are \(3\) outcomes equal to \(1\) among \(8\) observations, so the fitted common probability is \(\hat\mu=3/8=0.375\).
The intercept estimate is the logit of that probability:
For the independence working variance, the bread (the sum of \(D_i^TV_i^{-1}D_i\)) is \(8\hat\mu(1-\hat\mu)=8(0.375)(0.625)=1.875\). Each cluster's estimating contribution is its sum of residuals. The four contributions are \(1.25\), \(0.25\), \(-0.75\), and \(-0.75\), since each cluster has size two. The empirical meat is therefore
The uncorrected sandwich variance is \(2.75/(1.875)^2\approx0.7822\), giving a sandwich standard error of \(\sqrt{0.7822}\approx0.8844\) for the intercept. This example demonstrates the cluster-level calculation, not a recommendation to make reliable large-sample inferences from only four clusters. The point estimate is the overall observed proportion's logit here because the model has only an intercept; that simple equality need not hold for a general covariate model or a non-independence working correlation.
Assumptions, Diagnostics, and Alternatives
For interpretation and inference to be credible, the marginal mean model must represent the outcome probabilities adequately, clusters must be independent sampling units, and the number of independent clusters must support the chosen asymptotic inference. The binary variance formula follows from the Bernoulli outcome. A working correlation must be a valid correlation structure for the observations it represents, and the data design must contain enough information to estimate its parameters. Robust variance protects against some forms of working-correlation misspecification; it does not repair an incorrect mean model, dependence between clusters, or severe small-cluster problems.
Useful checks include comparing observed and fitted outcome proportions across important covariate and visit groups, checking that fitted probabilities remain between zero and one, and assessing whether the mean model omits important nonlinearities or interactions. Inspect the estimated working correlation for implausible values or boundary estimates, and check that the assumed structure matches the measurement schedule. If robust and model-based standard errors differ substantially, that signals sensitivity to the working covariance; it does not by itself identify which standard error is correct.
A practical comparison is with ordinary logistic regression and a logistic GLMM. Ordinary logistic regression uses the same link but, if it treats clustered observations as independent, generally does not provide valid uncertainty for dependence within clusters. A GLMM instead introduces random effects and estimates effects conditional on those effects under its distributional assumptions. GEE directly models marginal means and typically uses robust standard errors, without requiring a full random-effects distribution. The choice depends on the scientific estimand and design, not merely on which method returns a smaller standard error.
GEE also has limitations. Robust inference can perform poorly with few independent clusters; a robust standard error is not automatically reliable simply because it is called robust. Results may depend on the marginal mean specification, and the fitted model does not generally provide a full joint distribution for within-cluster binary outcomes. Further, a population-averaged odds ratio does not directly give an absolute probability difference: predicted probabilities or contrasts on the probability scale are needed for that question.
Common Mistake
Wrong reasoning: “I chose an exchangeable working correlation, so the GEE estimates a conditional treatment effect within each cluster.” A working correlation describes how the estimation weights residuals within a cluster; it does not condition the mean model on a cluster-specific random effect. Correction: interpret the logit coefficients from this model as population-averaged coefficients for the specified marginal mean. A conditional interpretation requires a model that explicitly defines a conditional mean, such as a GLMM.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- A treatment coefficient in a marginal logit GEE is \(-0.40\). What is the treatment odds ratio, and how should it be interpreted?
- In a binary GEE, what does the working correlation specify, and what does it not specify?
- Why does a sandwich variance calculation use cluster-level estimating contributions?
- A researcher selects an independence working correlation despite expecting within-person association. Does that choice automatically invalidate the marginal mean estimate? State an important qualification.
Solutions. (1) The odds ratio is \(\exp(-0.40)\approx0.6703\). At the same modeled covariate values, the estimated marginal odds in the treatment group are about 0.67 times the odds in the reference group. (2) It specifies a working pattern for within-cluster correlations used to form the covariance weights; it does not specify the full joint distribution of binary outcomes or turn the model into a conditional random-effects model. (3) Clusters are the independent units under the GEE setup, so the empirical variability of the estimating equation must be calculated across independent clusters rather than treating dependent observations as independent units. (4) Not automatically: with a correctly specified marginal mean model and the required regularity conditions, a working-correlation misspecification need not invalidate the coefficient estimator, though it can affect efficiency. Robust uncertainty and a sufficient number of independent clusters remain important.