Population-Averaged Alternatives: GLS, GEE, and Where the Linear Model's Luck Runs Out

Abstract

A marginal model does not explain where within-person correlation came from; it states the covariance directly and moves on. This post takes that refusal seriously. It separates the two clinical questions hiding behind one printout, the subject-specific comparison and the population-averaged one, then works through GLS and GEE: a mean model, a working correlation structure and a sandwich variance estimator, whose point estimate stays consistent when only the working correlation is wrong. With an identity link the two families agree exactly, which is why continuous outcomes forgive the distinction. With a logit link they do not: the same simulated model yields a conditional odds ratio of 6.05 and an exact population-averaged odds ratio of 3.72. Choose by estimand and missing-data assumption, not by habit.

---

Where we are. This series is about studies in which one participant contributes several rows of data — four skin sites here, several visits elsewhere — so the rows cannot be treated as independent. The usual repair is a mixed model: every person is given one personal offset, added unchanged to each of their rows, and that single number is what generates the correlation between a person's own measurements. Post 6 ended on a model that refuses to do this at all: gls states the marginal covariance matrix directly and never asks where the correlation came from. This post takes that refusal seriously and follows it to the end, because it is a real third option and not a degenerate one. For a continuous outcome the refusal costs almost nothing; for a binary outcome it changes the coefficient itself.

Nothing below assumes you have read Posts 1 to 6, and each term is re-introduced where it is first needed. If you would rather start from the gentler overview, the map of the whole series is here: series guide.

Every number in this series comes from simulated teaching data; later posts, including this one, introduce clearly labelled variant simulations. It is not an empirical finding, it describes no real cohort, and it must not be cited as evidence about skin physiology.

The ward round

Two clinicians are standing over the same simulated printout. The first asks about the patient in front of her: a woman with moderate atopic dermatitis whose barrier, from her own earlier readings, is already known to run leaky. She wants to know how much drier this patient's skin reads than that of a person whose barrier runs equally leaky but who does not have the disease. The second clinician asks a wider question: if the whole atopic-dermatitis population in the programme were set against the non-atopic population, how much drier would the skin be on average? Both are legitimate clinical questions, and neither of them is a counterfactual about an individual. The first compares two people held at the same person-level starting point; it does not compare this patient with a version of herself who never developed the disease. Both can be answered from the same simulated dataset.

They do not always have the same answer.

Simulated outcome	Conditional (subject-specific)	Population-averaged (marginal)
TEWL at the hand, AD − control (g·m⁻²·h⁻¹)	+8.0	+8.0
Xerosis at a site, AD vs control (odds ratio)	6.05	3.72

Every entry in that table is an exact property of the simulated model that generated it, including both odds ratios. The familiar attenuation formula, worked through later in this post, returns approximately 4.0 for the population-averaged odds ratio; integrating the same model exactly returns 3.72, and the difference between those two numbers belongs to the formula rather than to the model.

Nothing in that table is an error, and no row is a corrected version of the other. The first row is the coincidence that makes linear mixed models feel forgiving. The second row is where the coincidence expires.

Two philosophies, one table

Before the hard part: The table below sets two ways of modelling the same data side by side. The first builds the within-person correlation out of something you can name and estimate — a person effect, one number per patient. The second simply writes down how strongly a person's measurements travel together and moves on. It is the difference between explaining why two readings from the same patient resemble each other and merely allowing for the fact that they do. Only the first hands you a between-person spread to report; only the second is indifferent to where that spread came from.

Every model in this series so far has built the marginal covariance out of latent parts: one number per person, a variance for that number, and a residual variance. Post 4 wrote the recipe once and for all as $V_i = Z_i G Z_i' + R_i$. A marginal model declines to give a recipe. It writes $V_i$ down.

	Random-effects (conditional)	Marginal (population-averaged)
What you specify	$Z_i$, $G$, $R_i$	$V_i$ itself, or a working version of it
Where the correlation comes from	a shared latent quantity $b_i$	nowhere — it is simply asserted
What $\beta$ means	the change for a person, holding $b_i$ fixed	the change in the population mean
Person-level output	$\widehat{\tau}^{2}$, ICC, BLUPs, subject-specific prediction	none
Estimation	likelihood (REML/ML)	least squares or estimating equations

The two columns are not rivals for the same job. They are two answers to the question what is the correlation for? — a scientific quantity to be estimated, or a nuisance to be accommodated so that the mean structure is estimated honestly.

gls: saying the covariance out loud

Before the hard part: The code below does exactly one thing. It tells the software what the four-by-four table of variances and correlations among one person's four sites is allowed to look like, and it never mentions a person effect at all. Everything the dependence needs is supplied by hand, in two arguments.

The most direct marginal specification is the one Post 6 already introduced, and it needs no new machinery:

# unstructured marginal covariance, no latent person effect (Posts 6, 7)
fitU <- gls(tewl ~ site * group + age + sex + phototype,
            correlation = corSymm(form = ~ siteN | id),
            weights     = varIdent(form = ~ 1 | site), data = skin)

Read the two arguments as one sentence: corSymm supplies the six free correlations among the four sites, and varIdent supplies the four site-specific variances. Together they spend exactly the ten free elements of a symmetric $4\times4$ covariance matrix — the ceiling counted in Post 6. What this fit estimates is the fixed-effect vector $\beta$ and the ten covariance quantities; what it does not estimate is $\tau^2$, because there is no $b_i$ in the model to own it. The within-person dependence is fully accounted for. It is simply accounted for by assertion rather than by decomposition.

GEE: a mean model, a working correlation, and a sandwich

Before the hard part: A GEE is simpler than its notation. You write down what the average outcome should look like — the part you actually care about. You then make a guess about how one patient's repeated measurements hang together. Finally you let the observed spread of the data correct the standard errors in case that guess was poor. The equations below are only the formal way of saying those three things. If the algebra starts to blur, hold on to one sentence: a wrong guess about the correlation costs precision, not correctness, provided the average was modelled correctly.

Generalized estimating equations generalise that idea to outcomes that are not continuous. A GEE has exactly three components, and stating them separately is the whole of the method.

First, a mean model: a link function and a linear predictor, exactly as in an ordinary generalized linear model. Second, a working correlation structure, which describes how observations within a cluster are assumed to hang together. Third, a sandwich (robust) variance estimator, which is what makes the second component forgiving. The parameter estimates solve

\[ \sum_{i=1}^{N} D_i' V_i^{-1}\left(Y_i - \mu_i\right) = 0, \qquad D_i = \frac{\partial\mu_i}{\partial\beta'} \]

which is not a likelihood but a set of equations whose solution is consistent for $\beta$ whenever the mean model is correctly specified. Read the pieces from right to left: $Y_i - \mu_i$ is how far this cluster sits from where the mean model put it, $V_i^{-1}$ decides how heavily each of that cluster's residuals should count, and $D_i$ translates the result back onto the coefficient scale. The index $i$ runs over the $N$ independent clusters, which in the running example means the 80 participants and not the 320 rows. Every large-sample property of the method is stated in terms of that count rather than in terms of the number of observations. The working correlation lives only inside $V_i$, and $V_i$ appears only as a weight.

That is why a wrong working correlation does not break the point estimate. It costs efficiency, and it corrupts the model-based standard errors — which is where the third component earns its name:

\[ \widehat{\operatorname{Var}}(\widehat{\beta}) = \underbrace{\left(\textstyle\sum D_i'V_i^{-1}D_i\right)^{-1}}_{\text{bread}} \underbrace{\left(\textstyle\sum D_i'V_i^{-1}\widehat{r}_{i}\widehat{r}_{i}' V_i^{-1}D_i\right)}_{\text{meat}} \underbrace{\left(\textstyle\sum D_i'V_i^{-1}D_i\right)^{-1}}_{\text{bread}} \]

The two outer factors are what you would have used had the working correlation been right; the middle factor replaces the assumed covariance with the empirical cross-products of the observed residuals $\widehat{r}_{i} = Y_i - \widehat{\mu}_{i}$, cluster by cluster. The data are allowed to overrule the assumption about the variance. Whether that robustness is worth having in a given analysis is a judgement about design, not a property of the estimator — with 80 simulated participants it may reasonably be relied upon, and with a handful of clusters it may not.

Choosing a working correlation

Before the hard part: The working correlation is nothing more than your opening guess about how a patient's repeated readings resemble one another: not at all, all pairs alike, fading as the readings get further apart in time, or every pair on its own terms. The four options below are those four guesses. The caution that follows them is about how many patients you need before the sandwich can be trusted to rescue a bad guess.

Four choices cover almost all practice, and for the running example the arithmetic of Post 2 tells you which one is natural.

1. Independence — work as though the rows inside a cluster were unrelated, which is itself a working assumption rather than the absence of one, and let the sandwich do all the repair work. Valid, and often surprisingly serviceable, but wasteful when the correlation is strong.
2. Exchangeable — one correlation for every pair. This is the working analogue of the compound-symmetry matrix from Post 2, where every off-diagonal cell sat at 0.625.
3. AR-1 — correlation decaying with distance. It requires an ordered dimension; four exchangeable body sites do not provide one, so it belongs to the visit sub-study of Post 5 and not here.
4. Unstructured — six free correlations for $k=4$, estimated without further structure imposed on them, though still constrained to keep the working matrix positive definite, at the cost of six parameters and a larger demand on the number of clusters.

The honest caveat attaches to the sandwich rather than to the choice. The robust variance estimator is a large-cluster-count promise: it is consistent as the number of independent clusters grows, and with few clusters it is known to be anti-conservative, producing intervals that are narrower than they should be. Thus, when the number of clusters is small, a bias-corrected robust variance — available in most implementations — should be considered, and a likelihood-based mixed model may be the more defensible option.

The linear coincidence

Before the hard part: The whole section fits into one line. Some patients read high and some read low, and the person effect is built so that these two departures cancel on average. Average a straight line over departures that cancel, and the line comes back unchanged. That is why, for a continuous outcome, the "for a patient like this" answer and the "in the population" answer are the same number. The algebra below only checks that claim.

Post 3 promised that a random intercept moves the covariance and not the mean. That promise carried a clause — for a linear model — and this is the post where the clause falls due.

With an identity link, the conditional mean of a measurement is $\mu_{j,g} + b_i$ and the marginal mean is obtained by averaging over the person effect. Since $E(b_i)=0$ by construction,

\[ E(Y_{ij}) = E_b\!\left[\mu_{j,g} + b_i\right] = \mu_{j,g} + E(b_i) = \mu_{j,g} \]

so the fixed-effect coefficients of the conditional and the marginal model describe the same quantity on the same scale. That is the whole of the coincidence, and it is arithmetic rather than luck of the data: averaging a straight line over a mean-zero perturbation returns the same straight line. The person-level covariates that the code blocks in this series carry — age, sex and phototype — are constant within a person, sit entirely in the mean structure, and leave that argument untouched. On the simulated continuous outcome, the AD − control difference at the hand is +8.0 g·m⁻²·h⁻¹ whether you fit lme with random = ~ 1 | id, gls with an unstructured covariance, or GEE with an exchangeable working correlation; the three procedures agree on the estimand and differ only in how the correlation is book-kept. They need not return identical standard errors, because the variance machinery is genuinely different.

This is the reason so many analysts have never had to think about the distinction. Continuous outcomes forgive it.

Where the luck runs out

Before the hard part: Now move from a measured value to a yes/no outcome, and the coincidence breaks. Risks and odds do not average the way means do. Picture four patients whose personal risks of visible dry skin are 5%, 20%, 60% and 95%: the average of those four risks is not the risk of the average patient, because the risk scale is compressed at both ends. Averaging a curved relationship across people therefore flattens it, and a flatter relationship is a smaller odds ratio. Everything below is a careful accounting of how much smaller.

Now change the outcome. Extension C of the simulation keeps the same 80 participants and the same four sites, but records whether visible dry skin — xerosis — was present at each site. The generating model is a logistic random-intercept model:

\[ \operatorname{logit}\big[P(\text{xerosis}_{ij}=1 \mid b_i)\big] = \beta_0 + \beta_1\,\text{AD}_i + b_i, \qquad b_i \sim N(0,\ \tau^2) \]

with $\beta_0 = -1.00$, $\beta_1 = +1.80$ and $\tau^2 = 2.00$ in the simulation. The subject-specific structure is exactly the structure of Post 1 — one number per person, added unchanged to every row that person contributes — now added on the logit scale. Exponentiating $\beta_1$ gives a conditional odds ratio of 6.05.

The averaging step that was harmless above is no longer harmless, because the inverse-logit function is not linear:

\[ E_b\!\left[\operatorname{expit}(\eta + b)\right] \;\neq\; \operatorname{expit}\!\left(\eta + E_b[b]\right) \]

Averaging an S-shaped curve over a spread of person effects does not return the same S-shaped curve; it returns a flatter one. The population-averaged curve is therefore shallower than the subject-specific curve, and the population-averaged coefficient is smaller in absolute value. The standard approximation to the shrinkage factor, given by Zeger, Liang and Albert, is

\[ \beta^{M} \approx \beta^{C}\left(1 + c^2\tau^2\right)^{-1/2}, \qquad c = \frac{16\sqrt{3}}{15\pi} \approx 0.588 \]

which converts the conditional coefficient into its population-averaged counterpart using nothing but the between-person variance. Work the arithmetic through for the simulated values, one step at a time:

\[ c^2 = 0.588^2 = 0.346, \qquad 1 + 0.346\times 2.00 = 1.692, \qquad \sqrt{1.692} = 1.30 \]

\[ \beta^{M} \approx \frac{1.80}{1.30} \approx 1.38 \qquad\Rightarrow\qquad \text{marginal OR} = e^{1.38} \approx 4.0 \]

The same attenuation applied to the intercept gives the marginal risks that the formula implies: $-1.00/1.30 = -0.77$ for the control group and $(-1.00+1.80)/1.30 = 0.62$ for the AD group, which on the probability scale are approximately 32% and 65%. Reassembling those two risks into an odds ratio returns 4.0 again. It is worth saying plainly that this is not a check of anything. Both routes divide the same linear predictor by the same $\sqrt{1+c^2\tau^2}$, so the second collapses algebraically to $e^{\beta_1/\sqrt{1+c^2\tau^2}}$ — the first result written out on the probability scale and then read back. The agreement is forced by construction rather than earned. A series that keeps insisting output cannot validate the assumption that produced it should not accept a tautology as confirmation.

The genuine check is to stop approximating and integrate. Averaging $\operatorname{expit}(\beta_0 + \beta_1\text{AD}_i + b)$ over $b \sim N(0,\ 2.00)$ numerically gives exact marginal risks of 32.5% and 64.2%, marginal odds of 0.481 and 1.792, and an exact population-averaged odds ratio of 3.72. The approximation therefore overstates the marginal odds ratio by roughly 7% at this value of $\tau^2$ — small enough to leave the teaching point standing, and large enough that 4.0 and 3.72 must never be tabulated as though they were the same quantity.

One simulated truth. Two coefficients: 6.05 and 3.72.

Read both out loud

The two numbers stop being confusing the moment they are spoken as clinical sentences.

The conditional odds ratio says: for two patients with the same underlying tendency to dry skin, the one with atopic dermatitis has about six times the odds of visible xerosis at a given site. The comparison is held inside a person-level stratum, which is why it is called subject-specific.

The population-averaged odds ratio says: averaged over the whole simulated population, the atopic-dermatitis group has about 3.7 times the odds of visible xerosis. The comparison is between two groups of people, and the person-level heterogeneity has been integrated away rather than held fixed.

Two features of the arithmetic are worth stating flatly, because they remove most of the mystique. First, the marginal coefficient is always closer to the null than the conditional coefficient, and within the approximation the gap is governed entirely by $\tau^2$. When $\tau^2 = 0$ the factor becomes 1 and the two answers coincide exactly, which the exact calculation confirms. Second, computed exactly, the size of the attenuation depends not only on $\tau^2$ but also on where on the logistic curve the comparison sits — that is, on $\beta_0$ as well. That dependence is why 4.0 and 3.72 are not the same number. Attenuation here is not bias, not confounding, and not a small-sample artefact; it is what averaging a non-linear function does.

Which estimand does your question want?

This is the question to settle before opening any software, and it is a question about the research question rather than about the data.

Population-level and policy comparisons lean marginal. If the claim you intend to publish is of the form this group has a higher prevalence than that group, or if the whole clinic population were treated differently, the average outcome would move by this much, then the estimand you want is the population-averaged one, and a marginal model targets it directly. Individual prediction, variance partitioning, the ICC, BLUPs and any question in which the between-person heterogeneity is itself the scientific object require the random-effects formulation, because those quantities do not exist in a marginal model at all.

The subject-specific estimand may also be preferred when the clinical statement is genuinely within-person — a step-down decision for a patient already under follow-up, for instance, rather than a comparison of two populations. Where the two estimands are close to each other, as they will be whenever $\tau^2$ is small, the choice may reasonably be made on convenience. Where they are far apart, the choice must be made on the question, and it should be pre-specified rather than settled after seeing which number is larger.

What a GEE cannot give you

Four things are simply absent, and it is better to know that in advance than to discover it while writing the results section. First, there is no $\tau^2$: a working correlation parameter is not a variance decomposition, so no between-person and within-person split is available. Second, there is consequently no ICC in the variance-partition sense of Post 2. Third, there are no BLUPs and therefore no subject-specific predictions. Finally, there is no likelihood, which means no likelihood-ratio test and no AIC; model comparison uses QIC instead.

If the ICC is the scientific target rather than a nuisance parameter, a marginal model is the wrong tool — and it is worth remembering that the reliability ICC of a measurement study is a different quantity again, sharing a formula with the clustering ICC while answering a different question.

The missing-data caution

Before the hard part: One practical difference often settles the choice before any argument about estimands does, and it is about who drops out. Ask the clinical question directly: are the patients who stopped attending recognisable from what is already in the chart — a worse reading last visit, the sicker group, the older ones? If they are, a likelihood-based mixed model copes with that on its own, whereas a plain GEE needs the dropout to be essentially unrelated to anything. The paragraphs below state that precisely.

One further difference is frequently the deciding one, and it has nothing to do with estimands.

A likelihood-based mixed model is valid under the missing-at-random assumption, meaning that missingness may depend on the observed data. GEE with the usual sandwich variance requires the stronger missing-completely-at-random assumption unless inverse-probability weights are added. In longitudinal clinical studies dropout is very often related to observed history: the participants who stop attending are the ones whose earlier measurements were worst. That pattern is missingness at random rather than missingness completely at random, so the weaker assumption is usually the more plausible one. Thus, when a design is expected to lose participants in a way that depends on what has already been recorded, the likelihood-based mixed model or a properly weighted GEE may be the safer specification, and the reason should be written into the analysis plan rather than defended afterwards.

Dropout that depends on the unobserved value itself — missingness not at random, which is what the word informative is conventionally reserved for — is repaired by neither specification. It requires an explicit assumption about the missingness mechanism and a sensitivity analysis around that assumption, and the choice between a mixed model and a GEE is no longer the operative decision.

What this post does not license

Three over-readings are easy to reach from here, and each is wrong.

First, the identity-link coincidence is not a general rule about non-linear models. It is specific to the identity link, and the closest thing to an exception — the log link, where the person effect factors out multiplicatively and only the intercept absorbs the shift — proves the point rather than softening it. For the logit link the coefficients genuinely differ, and 6.05 and 3.72 are both correct.

Second, the attenuation formula is an approximation and is written with $\approx$ for that reason. On the simulated values it overstates the population-averaged odds ratio by about 7%, which is accurate enough for moderate $\tau^2$ to make the teaching point and to sanity-check a fitted model, and it should not be used to convert a published conditional odds ratio into a marginal one and then reported as if it had been estimated.

Third, "population-averaged" is a statement about averaging over the random effect, and nothing else. It does not mean marginal in the causal sense, it does not imply that confounding has been addressed, and a population-averaged association from an observational design remains an association. Any causal reading still requires the design and the adjustment set to carry it.

What to do in your own analysis

First, write the estimand down as a clinical sentence before choosing the software, and check that the sentence contains either the words for a patient or the words in the population. Second, for a continuous outcome from a balanced design, either specification is defensible, so state which one you fitted and why rather than leaving the reader to infer it. Third, for a binary or count outcome, label every odds ratio in the results table explicitly as conditional or population-averaged, because the two are routinely tabulated side by side across papers as though they were the same quantity. Fourth, if you report a conditional odds ratio, report $\widehat{\tau}^{2}$ alongside it, since the distance between the two scales is a function of that number and a reader cannot reconstruct it otherwise. Fifth, when dropout is expected to depend on what has already been observed, prefer the likelihood-based mixed model or a weighted GEE and say so in the protocol, and where dropout may depend on the unobserved value itself, plan the sensitivity analysis at the same time. Finally, when the number of clusters is small, consider a bias-corrected robust variance and treat the nominal coverage of an uncorrected sandwich interval with caution.

Key takeaways

• A marginal model does not explain where the within-person correlation came from; it states the covariance directly, and the dependence is fully accounted for either way.
• A GEE has exactly three components — a mean model, a working correlation structure, and a sandwich variance estimator — and the point estimate remains consistent when only the working correlation is wrong.
• With an identity link the conditional and marginal fixed effects coincide because $E(b_i)=0$, which is why continuous outcomes forgive the distinction entirely.
• With a logit link they do not coincide: in the simulated example the same generating model yields a conditional OR of 6.05 and an exact population-averaged OR of 3.72, with marginal risks of about 32% and 64%. The Zeger–Liang–Albert approximation puts that second number at about 4.0, some 7% too high, and the two should be labelled accordingly.
• Within that approximation the attenuation is governed by $\tau^2$ alone, it always moves the marginal estimate toward the null, and it vanishes exactly when $\tau^2 = 0$.
• GEE returns no $\tau^2$, no ICC, no BLUPs and no likelihood, so the choice between the two families should be settled by the estimand and by the missing-data assumption, not by habit.

Next. Post 8 takes everything the series has built and points it at one output column: four identical standard errors, and the professor who wants to know whether the software is broken.

References

1. Laird NM, Ware JH. Random-effects models for longitudinal data. Biometrics. 1982;38(4):963–974.
2. Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika. 1986;73(1):13–22.
3. Zeger SL, Liang KY, Albert PS. Models for longitudinal data: a generalized estimating equation approach. Biometrics. 1988;44(4):1049–1060.
4. Verbeke G, Molenberghs G. Linear Mixed Models for Longitudinal Data. Springer; 2000.
5. Diggle PJ, Heagerty P, Liang KY, Zeger SL. Analysis of Longitudinal Data. 2nd ed. Oxford University Press; 2002.
6. Hubbard AE, Ahern J, Fleischer NL, et al. To GEE or not to GEE: comparing population average and mixed models for estimating the associations between neighborhood risk factors and health. Epidemiology. 2010;21(4):467–474.