Generalised estimating equations

Ordinary regression assumes every observation carries independent information. That assumption fails the moment the same person is measured twice, or two patients share a ward. When it fails, the coefficients are usually still reasonable but the uncertainty around them is wrong, and wrong in the direction of making findings look stronger than they are. The scale of the error is set by how much the observations repeat each other, which is covered in detail in the guide to clustered data and the unit of analysis.

Generalised estimating equations, introduced by Liang and Zeger in 1986, extend the generalised linear model to this situation. You still choose a link function and an outcome distribution as you would for logistic regression or Poisson regression. GEE adds a statement about how observations within a cluster are correlated, uses it to weight the data more sensibly, and then computes standard errors in a way that protects you if that statement was wrong.

This is the part most often got wrong in write-ups, and it is not a technicality. GEE answers a question about populations. A mixed effects model answers a question about individuals.

Suppose a service change is associated with better outcomes. The GEE coefficient says what happens to the average outcome across the whole population of wards if every ward adopts the change. The mixed model coefficient says what happens within a given ward when that ward adopts it. Both are legitimate. They are different numbers answering different questions, and a paper that reports one while describing the other has misstated its own finding.

For continuous outcomes with an identity link the two coincide, which is why the distinction is easy to forget. For binary outcomes they do not. The marginal odds ratio from GEE is pulled towards 1 compared with the conditional odds ratio from a mixed model, and the gap widens as variation between clusters increases. Neither is wrong. They are measuring different things.

Fitting a GEE means stating three things. Every one of them belongs in your methods section.

A fourth decision arrives if you are selecting between models. The usual likelihood-based criteria do not apply, because GEE is not a likelihood method. QIC, proposed by Pan in 2001, is the quasi-likelihood analogue of AIC and is what model selection should be based on.

The working correlation describes how observations within a cluster hang together. Choose it to match the design rather than by trying several and keeping the one that gives the best p-value.

The robust or sandwich estimator, also called the Huber-White or empirical variance estimator, computes uncertainty from how much the clusters actually vary rather than from the correlation structure you assumed. That is what lets you be wrong about the structure without being wrong about the answer.

The practical consequence is that your confidence intervals come from the observed between-cluster variation. If you fit a GEE with model-based standard errors instead, you have thrown away the protection and are relying on the working correlation being correct, which it almost never is. Software defaults vary, so check which one you have actually used rather than assuming.

There is one condition. The sandwich estimator is a large-sample result, and "large" here means the number of clusters, not the number of rows. Fifty thousand observations from twelve clinics is a small sample for this purpose.

Forty clinics, fifteen patients each, so 600 observations. The exposure is applied at clinic level and the intraclass correlation is 0.30. Here is what happens to the standard error depending on whether the clustering is respected.

Step 1. The naive analysis

Treat all 600 observations as independent and compare the two groups.

standard error = 0.082

Step 2. The analysis that respects the design

Account for the fact that the 600 observations came from only 40 independently allocated units.

standard error = 0.188

Step 3. Compare, and check against theory

0.188 / 0.082 = 2.29

The honest standard error is 2.29 times larger. Now check that against the design effect:

design effect = 1 + (m − 1) × ICC = 1 + 14 × 0.30 = 5.20
√5.20 = 2.28

The two agree, as they must. Standard errors inflate by the square root of the design effect, and that is precisely the correction the robust estimator applies.

Step 4. What that does to a p-value

A t statistic is the estimate divided by its standard error, so dividing the standard error by 2.29 multiplies the t statistic by 2.29. An effect that would sit at t = 1.3 and be reported as non-significant becomes t = 2.98 under the naive analysis, with a p-value around 0.003. The naive model does not produce a slightly optimistic answer. It produces a confident one that is wrong.

Take a true subject-specific model with a log-odds ratio of 0.80, so a conditional odds ratio of 2.23, and vary how much clusters differ from each other. Averaging the probabilities across the population gives the marginal effect a GEE would estimate.

Reading the table

With no variation between clusters the two coincide exactly. As clusters become more different, the marginal odds ratio is pulled towards 1 while the conditional one stays where it is. At a between-cluster variance of 1.0, a conditional odds ratio of 2.23 corresponds to a marginal one of 1.95, which is an attenuation of about 17% on the log scale.

Why neither is wrong

The conditional estimate answers: if this clinic adopted the intervention, what happens in this clinic. The marginal estimate answers: if every clinic adopted it, what happens to the population average. Those are different quantities and both are useful, which is why a paper must say which one it is reporting.

Note that this divergence only arises for non-linear links. With an identity link on a continuous outcome the two coincide regardless of how much clusters differ, which is why the distinction is easy to forget until a binary outcome appears.

Around forty clusters is the figure usually quoted as the point below which the robust estimator becomes unreliable, and below it the standard errors are biased downwards, which defeats the entire purpose of using them.

If you have fewer, you have three reasonable options. Apply a small-sample correction such as Mancl and DeRouen, Kauermann and Carroll, or Fay and Graubard, all available in standard software. Use a mixed effects model instead, which handles small cluster counts more gracefully. Or aggregate to cluster level and analyse the summaries, which is crude but almost impossible to get wrong.

The honest answer is that your research question decides, not your software or your training.

A reasonable rule of thumb for applied health and social research: if the question is about what a policy would do across a system, GEE fits. If it is about what happens to a patient, a pupil or a ward, a mixed model fits.

This is the most important limitation and the one least often disclosed. Standard GEE is valid when data are missing completely at random, a stronger assumption than the missing at random condition under which likelihood methods including mixed models remain valid.

In practice that assumption is often indefensible. People drop out of studies for reasons connected to their outcomes, and app users stop logging when they feel better or when they feel far worse. If dropout relates to anything you observed, standard GEE can be biased where a mixed model would not be. The repair is weighted GEE using inverse probability weights, which needs a model for why data went missing and should be stated as the extra assumption it is.

R

geepack is the usual choice. The data must be sorted by the cluster identifier before fitting, or the correlation structure is applied to the wrong groupings.

library(geepack)
df <- df[order(df$clinic), ]
fit <- geeglm(outcome ~ treatment + age, id = clinic,
    family = binomial, corstr = "exchangeable", data = df)
summary(fit)

geeglm reports robust standard errors by default, which is what you want. Model selection uses QIC rather than AIC: MuMIn::QIC(fit). For a small number of clusters, saws or geesmv provide small-sample corrections to the sandwich estimator.

SPSS

GEE sits under Analyze, Generalized Linear Models, Generalized Estimating Equations. The subject variable defines the clusters and the within-subject variable defines the ordering for structures such as AR(1).

GENLIN outcome (REFERENCE = FIRST) BY treatment WITH age
  /MODEL treatment age DISTRIBUTION = BINOMIAL LINK = LOGIT
  /REPEATED SUBJECT = clinic CORRTYPE = EXCHANGEABLE
    COVB = ROBUST
  /PRINT CPS DESCRIPTIVES MODELINFO FIT SUMMARY SOLUTION.

COVB = ROBUST is the important line. SPSS will accept MODEL instead, which gives model-based standard errors and removes the protection the method exists to provide.

Stata

xtgee is built in and needs the panel structure declared first.

xtset clinic
xtgee outcome treatment age, family(binomial) link(logit) corr(exchangeable) vce(robust)
xtcorr

vce(robust) requests the sandwich estimator. xtcorr prints the estimated working correlation matrix, which is worth looking at: if an exchangeable structure returns a correlation near zero, the clustering may be negligible, and if it returns something very high your effective sample size is far smaller than the row count suggests.

Reviewers want evidence that each choice was deliberate. Five sentences usually suffice.

In R this is geepack::geeglm with id and corstr; in Stata xtgee with vce(robust); in SPSS the GEE option within Generalized Linear Models; in SAS PROC GENMOD with a REPEATED statement. If an analysis has already been written up and you want it checked before submission, that is what an independent statistical review covers.