Stats for
Menu

Generalized and Multilevel Regression

Generalized Estimating Equations (GEE)

GEE estimates average effects on outcomes measured repeatedly or in clusters.

Generalized estimating equations take an outcome measured repeatedly or in clusters, such as yes/no answers at several check-ins, and test how the population's average outcome differs with the predictors, allowing for correlation within clusters.

Example

Does a four-week sleep-hygiene program raise the chance of a good night's sleep? We randomly assigned 150 adults with poor sleep to the program or a waitlist, and each answered “Did you sleep well last night?” at four weekly check-ins.

Repeated or Clustered Measurements + Binary Outcome + One or More Predictors → Generalized Estimating Equations (GEE)
Tests
Population-averaged (marginal) regression coefficients, here as odds ratios
Observations
Nested
Parametric
No
Headline Statistic
OR
Effect Size
Population-Averaged Odds Ratio

The Basics

When to Use It

Use it when each person (or couple, classroom or clinic) contributes several observations of the same outcome, often yes/no, and you want to know how the average outcome differs with the predictors across the population rather than for any one person. With the right link, GEE also handles counts and continuous scores. For example:

  • Does a four-week sleep-hygiene program raise the chance of a good night's sleep?
  • Across 12 weekly sessions, are clients who complete their CBT homework more likely to report a week without a panic attack?
  • In classrooms that adopt a mindfulness curriculum, are students less likely to be sent out of class for disruptive behavior?

The question GEE answers is marginal, or population-averaged: how much higher are the odds of a good night across everyone in the program? A multilevel logistic regression answers the conditional one: how much higher are the odds for a given person in the program? With yes/no outcomes those are different numbers, and the right one depends on whether you're describing a population or individuals (Zeger et al., 1988). Because GEE doesn't model the person effects at all, it also makes fewer assumptions than a mixed model (McNeish et al., 2017). It asks more of the data in two ways: missing answers must be missing completely at random, where a mixed model only needs missingness to depend on answers you did observe, and it needs enough people (or groups) for its robust standard errors to be trusted.

GEE describes only the average outcome and how answers vary and correlate, not a full probability model for each person's answers. With no likelihood, it tests coefficients with Wald tests and compares models with the quasi-likelihood criterion QIC rather than likelihood-ratio tests and AIC (Pan, 2001).

Vocabulary You Need

Definition Population-Averaged Effect

The difference in the outcome, averaged over everyone in the population, between people who differ on a predictor: for example, how much higher the share of good outcomes is among all the people who got a treatment than among all who didn't. Also called a marginal effect, it's what generalized estimating equations estimate. Its counterpart is the subject-specific (conditional) effect that mixed models estimate: the difference among people who share the same random effect. For a continuous outcome in a linear model the two coincide. For a yes/no outcome on the odds-ratio scale they don't: averaging many people's S-shaped curves flattens them, so the population-averaged odds ratio sits closer to 1. Neither is wrong; they answer different questions. Glossary →

Definition Working Correlation (\alpha)

The pattern of correlation that generalized estimating equations assume among the observations within one cluster, such as one person's answers at different times, used to decide how much weight each observation gets. Common choices are independence (no correlation), exchangeable (the same correlation, α, between any two observations in a cluster) and first-order autoregressive, AR-1 (α between neighboring times, α² two steps apart, and so on). The size of α is estimated from the model's residuals. It's called “working” because it needn't be right: as long as the model for the average is right, the coefficients stay consistent whichever structure you choose, robust standard errors stay valid, and a closer guess mainly makes the estimates more precise. The exception is a predictor that changes over time and may itself depend on earlier outcomes; then a structure other than independence can bias the coefficients, so use independence. Glossary →

Definition Robust Standard Error (SE_{\text{robust}})

A standard error built from the residuals themselves rather than from what the model assumes about how observations vary and correlate. It's also called the sandwich estimator, after its formula: the model's own idea of how precise the estimates are on both sides (the “bread”), and the spread the residuals actually show, added up cluster by cluster, in the middle (the “filling”). In generalized estimating equations it keeps the standard errors valid even when the working correlation is wrong, as long as the clusters are independent and there are enough of them. With few clusters, a few dozen or fewer, it tends to run too small, and bias-corrected versions exist. In ordinary regression the same idea gives standard errors that hold up when the spread of the outcome isn't the same everywhere (heteroscedasticity-consistent, or HC, standard errors). Glossary →

Assumptions

1. Independent Observations The people (or other units) in your study don't influence each other.
What It Means
Each person's data comes from a separate, unconnected source, so knowing one person's scores tells you nothing about anyone else's. In designs that measure people more than once, their own scores are meant to be related; it's different people who must be independent. Methods built for clustered data, such as multilevel models and GEE, also allow members of the same group to be related, as long as the analysis names the grouping; then it's the groups that must be independent.
Applies To
Observations
How to Check
By design, not by a test. Did each person contribute their own data, with no clusters such as roommates, couples, therapy groups or classrooms whose members influence each other? If there are clusters, does the analysis model every one of them?
How Much It Matters
A lot, but only between clusters. GEE is built for one person's answers being related, so those can correlate however they like. What it needs is for different people to be independent. If people come in groups that share something, such as couples or members of the same online cohort, the robust standard errors miss that sharing and come out too small.

If It's Violated

  • Correction Use the group (the couple, the cohort) as the cluster instead of the person. The robust standard errors then allow for any correlation within the group, including a person's own repeated answers.
  • Design Change Deliver the program individually, or randomize and analyze whole groups.
2. Linearity of the Logit Each continuous predictor relates to the log-odds of the outcome in a straight line.
What It Means
The model assumes that each one-unit rise in a continuous predictor changes the log-odds of the outcome by the same amount, whether it happens low or high on the predictor's range. On the probability scale that's an S-shaped curve, so the probabilities themselves aren't expected to rise in a straight line. A U shape (risk that is high at both ends), a threshold or a leveling off breaks it. Yes/no predictors meet it automatically, since they have only two values.
Applies To
Relationship
How to Check
Add a bend and see whether it helps: a squared term, a spline, or time entered as categories, compared with the straight-line model by a Wald test. GEE has no likelihood, so there's no likelihood-ratio test. A plot helps more than any p: plot the log-odds of the outcome at each value or band of the predictor (here, each week in each group); the points should fall near a line.
How Much It Matters
A lot, because the mean model is the one thing GEE needs right. The robust standard errors forgive a wrong working correlation, not a wrong model for the average. Week is the only continuous predictor, and the waitlist's weekly percentages wobble (54.7% at week 3, 50.7% at week 4), but entering week as four categories instead of a straight line doesn't fit reliably better, Wald χ²(2) = 3.39, p = .18, and the program coefficient barely moves (0.545 against 0.543). The mean model also needs the right terms. If good nights rose in the program group and not in the waitlist, a model without a program × week interaction would give one program effect that fits neither week 1 nor week 4. Here the gap narrows: adding a program × week term puts the program's advantage at about 0.72 log-odds at week 1, shrinking by 0.12 a week, but that change is well within chance (p = .35). The data can't tell whether the effect really fades, so we report the average and show the weekly percentages.

If It's Violated

  • Correction Add a program × week interaction to let the groups change at different rates, or enter week as a set of categories so the trend needn't be a straight line on the log-odds scale.
  • Proceed with Caveat If the program's effect varies a little across weeks, report the single odds ratio as the average over the four check-ins, and show the weekly percentages alongside it.
3. Enough Independent Clusters You have enough separate people or groups, not just enough observations, for the estimates you need.
What It Means
Anything estimated from differences between clusters rests on the number of clusters, not the number of observations inside them. In a GEE or any analysis with robust (sandwich) standard errors, the standard errors are learned from the clusters themselves; with a few dozen clusters or fewer they tend to be too small, so p-values and confidence intervals are too optimistic. In a multilevel (mixed) model, the standard errors come from the model instead, but the random-effect variances and any predictor measured at the cluster level are still estimated from differences between clusters, so with few clusters they are imprecise, the variances can be biased (Moineddin et al., 2007, found this for yes/no outcomes) and their standard errors come out too small.
Applies To
Design
How to Check
Count the independent clusters (people when each person gives repeated answers, groups when people are nested in groups), not the rows. Suggested minimums vary by method, from about 30 for sandwich standard errors (Halekoh et al., 2006) to about 50 for a multilevel model's cluster-level standard errors (Maas & Hox, 2005). With robust standard errors and only a few dozen clusters, use a small-sample correction such as Mancl and DeRouen's bias-corrected estimator; in their simulations with 10 to 40 clusters, the jackknife, like the uncorrected sandwich, could give false-positive rates two to three times the intended level. With a mixed model, add intervals that don't lean on large-sample normality, such as profile-likelihood or parametric-bootstrap intervals, and for a yes/no or count outcome check whether the results hold with more quadrature points.
How Much It Matters
Little with 150 people. The robust standard errors are trustworthy only when there are many independent clusters; with a few dozen or fewer they run too small and p comes out too small (Mancl & DeRouen, 2001). geepack's authors put the danger zone at about 30 clusters or fewer (Halekoh et al., 2006). Adding check-ins per person doesn't help: it's the number of people that counts.

If It's Violated

  • Correction Use Mancl and DeRouen's (2001) bias-corrected sandwich estimator, which kept false-positive rates close to 5% with as few as 10 clusters in their simulations. statsmodels offers it as cov_type="bias_reduced".
  • Correction A jackknife variance estimate leaves out one cluster at a time; in R's geepack, set std.err = "jack" (Halekoh et al., 2006). Treat it as a second opinion, not a fix: in Mancl and DeRouen's (2001) simulations with 10 to 40 clusters and a yes/no outcome, it too gave false-positive rates up to 2–3 times 5%.
  • Design Change Recruit more people (or more groups) rather than more measurements per person.
4. Missing Completely at Random Whether a value (an answer, a code, a score) is missing has nothing to do with what it would have been.
What It Means
Methods that use only the answers people gave, such as GEE, stay consistent (close to right in large samples) when missing answers are missing completely at random: the chance an answer is missing is unrelated to the answer itself and to anything else in the model. If people tend to skip check-ins when things are going badly, the answers that remain paint too rosy a picture.
Applies To
Observations
How to Check
By design and by looking. Count the missing observations by group and time, and compare people with and without gaps on their earlier answers and background. Missingness that tracks earlier answers rules this out; nothing in the data can prove it holds.
How Much It Matters
None here, since everyone answered all four check-ins. In general, a lot: GEE's estimates stay consistent when answers are missing completely at random (Halekoh et al., 2006). If people who slept badly were more likely to skip a check-in, the remaining answers would overstate how often people slept well, and the program effect could be biased.

If It's Violated

  • Alternative Test Multilevel Logistic Regression. A mixed model fitted by maximum likelihood uses every available answer and stays valid when missingness depends on answers you did observe, though it estimates person-specific rather than population-averaged effects.
  • Correction Model each answer's chance of being observed from the person's earlier answers and predictors, and weight each observed answer by the inverse of that chance. The estimates stay population-averaged, and they're consistent if what drives dropout is in that model.
  • Correction Impute the missing answers several times from each person's other answers and predictors, fit the GEE to each completed dataset, and combine the results.
  • Proceed with Caveat Report how many check-ins are missing in each group and whether the people who missed them had been sleeping worse beforehand.

Limitations

This Is Not Your Test If…

What It Doesn't Tell You

  • About any one person. The odds ratio of 1.72 is an average over the population; it doesn't say each person's odds rose that much, or that everyone benefited.
  • Why. Random assignment supports reading the gap as the program's effect, but a waitlist doesn't rule out expectation: people who knew they were in a program may have judged their nights more kindly.
  • Whether it matters. Significance doesn't say whether 13 more good nights in every 100 is worth four weeks of changing bedtime habits.
  • Who else. The result describes adults with poor sleep like these 150, and only as well as one yes/no question captures a night's sleep.

Common Mistakes

  • Don'tRead a GEE odds ratio as the change in any one person's odds.

    InsteadDescribe it as a population average, such as “across all participants, the odds of a good night were 1.72 times higher.”

  • Don'tCompare a GEE odds ratio with a mixed model's odds ratio as if they measured the same thing.

    InsteadExpect the mixed model's person-specific odds ratio to be further from 1, and say which kind each study reports.

  • Don'tStack everyone's repeated yes/no answers into an ordinary logistic regression.

    InsteadUse GEE or a mixed model. Ignoring the clustering here would have shrunk the program's standard error by a quarter and turned p = .015 into p = .001.

  • Don'tReport the model-based (“naive”) standard errors from a GEE.

    InsteadReport the robust (sandwich) standard errors, which don't depend on the working correlation being right.

  • Don'tHunt through working correlations for the smallest p.

    InsteadChoose one from the design before fitting, report it, and show the others only as a sensitivity check.

  • Don'tDescribe a regression coefficient as the effect of raising the predictor.

    InsteadSay the prediction is that much higher per unit of the predictor, on the model's scale (the outcome itself, or its log-odds or log rate in a logistic, ordinal or Poisson model). Causal language needs a design that supports it, such as an experiment.

  • Don'tRead an odds ratio of 3 as “three times as likely.”

    InsteadSay “three times the odds,” and give the two percentages so readers can see the actual chances.

  • Don'tRead a p-value as the chance that the null hypothesis is true.

    InsteadRead it as: if the null were true, how often data at least this extreme would turn up.

Result Decoder

Significance and Effect Size Guidelines

Benchmarks Significance (p Value)
  • p < .001

    Very strong evidence against the null hypothesis.

  • .001 ≤ p < .01

    Strong evidence against the null hypothesis.

  • .01 ≤ p < .05

    Statistically significant at α = .05. Moderate evidence against the null hypothesis.

  • .05 ≤ p < .10

    Not significant. Weak evidence; sometimes called a “trend.” Treat it as inconclusive.

  • p ≥ .10

    Not significant. Little evidence against the null, which is not the same as evidence of no effect.

These evidence labels are one common convention for describing p in graded language rather than as significant or not; the boundaries are rough guides, not rules, and other authors draw them differently. Muff et al. (2022) argue for this kind of language of evidence over a single cutoff. α is chosen before the analysis, and p never says how big or important an effect is.

More about statistical significance ↓

Benchmarks Effect Size (Population-Averaged Odds Ratio)
  • 0.60 < OR < 1.68 Below Small

    Below the odds ratio that matches d = 0.2 when 1% of the comparison group has the outcome.

  • 1.68 ≤ OR < 3.47 (0.29 < OR ≤ 0.60) Small

    About d = 0.2 when 1% of the comparison group has the outcome.

  • 3.47 ≤ OR < 6.71 (0.15 < OR ≤ 0.29) Medium

    About d = 0.5 when 1% of the comparison group has the outcome.

  • OR ≥ 6.71 (OR ≤ 0.15) Large

    About d = 0.8 or more when 1% of the comparison group has the outcome.

These benchmarks were worked out for odds ratios comparing two separate groups when the outcome is rare (Chen et al., 2010), so treat them as a rough guide. A good night is common (46.7% of the waitlist's check-ins), and a common outcome needs a smaller odds ratio for the same Cohen's d, so the label can understate the effect: at a 46.7% rate, an odds ratio of 1.72 matches d ≈ 0.34, still small but well above the d = 0.2 where “small” starts. These benchmarks compare whole groups, as a population-averaged odds ratio does; a mixed model's person-specific odds ratio for the same data sits further from 1 and would look larger against them. Below 1, invert the odds ratio before comparing.

More about effect size ↓

How It's Reported

Equation Form

Example

b = 0.54, SErobust = 0.22, OR = 1.72, 95% CI [1.11, 2.66], z = 2.44, p = .015

Table 1. GEE Logistic Model Predicting a Good Night's Sleep

Example
bSE_{\text{robust}}zpOR95% CI for OR
Intercept −0.340.18−1.86.0630.71[0.50, 1.02]
Program (vs waitlist) 0.540.222.44.0151.72[1.11, 2.66]
Week (per week) 0.140.062.21.0271.15[1.02, 1.29]

Note. N = 150 adults, 600 check-ins. Exchangeable working correlation, α̂ = .27; robust (sandwich) standard errors. Program is coded 1 = program, 0 = waitlist; week is counted from the first check-in (0 to 3), so the intercept is the waitlist's log-odds at week 1. OR = population-averaged odds ratio. The program row is the one worked through in the Napkin Math.

Table 2. Estimated Working Correlation between Check-Ins

Example
Week 1Week 2Week 3Week 4
Week 1 1.00.27.27.27
Week 2 .271.00.27.27
Week 3 .27.271.00.27
Week 4 .27.27.271.00

Note. The exchangeable structure assumes the same correlation between any two of a person's check-ins, estimated from the residuals as α̂ = .27. The observed correlations between weeks' answers range from .20 to .33, with no sign of fading as check-ins get further apart, which fits the exchangeable choice.

Table 3. The Program Effect under Three Working Correlations

Example
\hat{\alpha}b_{\text{program}}SE_{\text{robust}}ORpQIC
Independence 0 (fixed)0.5420.22291.72.015823.59
Exchangeable .270.5430.22291.72.015823.59
AR-1 .270.5550.22391.74.013823.73

Note. Same model, same robust standard errors; only the working correlation changes. α̂ is the estimated correlation (for AR-1, between neighboring weeks, using statsmodels' lag-1 moment estimator, Autoregressive(grid=True); statsmodels' default AR-1 estimator and R's geepack fit α to all pairs of weeks and get α̂ ≈ .36, b ≈ 0.56, SE ≈ 0.225). QIC = quasi-likelihood information criterion (Pan, 2001), computed with R geepack's formula and the binomial scale fixed at 1 (statsmodels' qic() uses a different constant, so its values differ); smaller is better, but differences under 1 point, as here, can't tell the structures apart. A fourth common choice, unstructured, gives every pair of weeks its own correlation (6 here); it suits designs with few times and many people.

Full Text (APA Style)

Example

Across the four weekly check-ins, adults in the sleep-hygiene program reported a good night's sleep 60.0% of the time and those on the waitlist 46.7% of the time. In a GEE logistic model with an exchangeable working correlation (α̂ = .27) and robust standard errors, controlling for week, the population-averaged odds of a good night were 1.72 times as high in the program as on the waitlist, b = 0.54, SErobust = 0.22, OR = 1.72, 95% CI [1.11, 2.66], small by convention (equivalent to about d = 0.34 at this base rate). The difference was statistically significant, z = 2.44, p = .015. Good nights also became more common from week to week, OR = 1.15 per week, 95% CI [1.02, 1.29], p = .027, mainly on the waitlist (38.7% at week 1 and 50.7% at week 4, against 58.7% and 61.3% in the program).

  • Say it's a GEE, and name the link (logit), the working correlation and the kind of standard errors (robust). Readers need all three to judge the result.
  • Give the observed percentages too. An odds ratio of 1.72 means little until readers see 47% against 60%.
  • Call the odds ratio population-averaged, or say “across participants,” so no one reads it as the change for each person.
  • With more than one predictor, put the coefficients in a table like Table 1: b, its robust SE, z (or the Wald χ²), p, and the odds ratio with its CI.

Reading Each Piece

The equation form above, piece by piece, in this order: how big, how certain, then how surprising.

b = 0.54 Across all the adults, the program raises the log-odds of a good night by 0.54 at any given week. Log-odds are hard to picture, which is why the odds ratio below exists.
α̂ = .27 After allowing for program and week, any two of one person's check-ins correlate about .27. That's enough that the four answers count for less than four separate people would.
SErobust = 0.22 The program coefficient would typically wobble by about 0.22 log-odds from one study of 150 people to the next. Ignoring that people's answers go together would have claimed 0.17.
OR = 1.72 Across the population, the odds of a good night are 1.72 times as high in the program, at any given week: at week 1, about 55% against 42%. That's small by the benchmarks above, which were built for rare outcomes; at this outcome's base rate it matches a Cohen's d of about 0.34, still small. Whether 13 more good nights in 100 is worth four weeks of sleep hygiene is a question about sleep, not statistics. Benchmarks ↑
95% CI [1.11, 2.66] The population's odds ratio is plausibly anywhere from 1.11, barely an edge, to 2.66, a sizable one. The whole range is above 1.
z = 2.44 The program's coefficient is 2.44 robust standard errors above zero. R and SPSS report the same test as a Wald χ²(1) = 5.93.
p = .015 If the program made no difference to the population's odds, a coefficient at least this far from zero, in either direction, would turn up about 1.5% of the time. That's below .05, so the result is statistically significant. Benchmarks ↑

Calculations

The Running Example

Example

The sleep-hygiene study: 150 adults with poor sleep, 75 in a four-week sleep-hygiene program and 75 on a waitlist, each answering yes or no to “Did you sleep well last night?” at four weekly check-ins, 600 answers in all.

Null Hypothesis:
\beta_{\text{program}} = 0 \quad (\text{OR} = 1)
Alternative (Two-Sided):
\beta_{\text{program}} \neq 0 \quad (\text{OR} \neq 1)

\beta_{\text{program}} (“beta”) is the population-averaged program coefficient: how much higher the log-odds of a good night are, across all adults like these, in the program than on the waitlist, at the same week. The sample's 0.54 estimates it. A coefficient of 0 is an odds ratio of 1, no difference. The same model also tests the week coefficient (Table 1), but the program is the study's question.

What Your Data Should Look Like

Long format: one row per observation, so one row per person per check-in. Each row carries the cluster ID (whose answer it is), the outcome (1 = yes, 0 = no), the time and the predictors. Keep each cluster's rows together and in time order; a missed check-in simply has no row.

Example
participantconditionprogramweekcellgood_night
1control01control_11
1control02control_20
1control03control_31
1control04control_41
2control01control_10
2control02control_20
… 594 more rows
nYes%
Waitlist, Week 1752938.7%
Waitlist, Week 2753242.7%
Waitlist, Week 3754154.7%
Waitlist, Week 4753850.7%
Program, Week 1754458.7%
Program, Week 2754154.7%
Program, Week 3754965.3%
Program, Week 4754661.3%

Simulated data.

Download all 600 rows (data.csv)

Napkin Math

The whole test, one idea at a time. Under each heading, a line says why the step matters, and the green note reads the equation in plain English with the sleep-hygiene study's numbers. Collapse all to read just the headings: they tell the whole story on their own.

  1. 1. Count How Often Each Group Slept Well

    Start with the plain percentages. Each group gave 300 answers, 4 from each of 75 people, and the share of yes answers is what the model will describe.

    \hat{p} = \dfrac{x}{n} =\dfrac{140}{300} \text{ and } \dfrac{180}{300}

    \hat{p}_{\text{waitlist}} = .47

    \hat{p}_{\text{program}} = .60

    x is the number of yes answers and n the 300 answers in each group. The waitlist said yes 140 times, 140 ÷ 300 = .47; the program group 180 times, 180 ÷ 300 = .60. That's 13 more good nights in every 100.

  2. 2. Write a Model for the Population's Chance of a Good Night

    A yes/no outcome is modeled on the log-odds scale. The odds of a good night are its chance divided by the chance of a bad one: a .60 chance is odds of .60 ÷ .40 = 1.5. Log-odds are the natural log of the odds; they run from minus to plus infinity, with 0 at a 50% chance, so effects can add up in a straight line. GEE fits this model like logistic regression but weights each person's four answers together, instead of treating them as four strangers.

    \ln \dfrac{P(y_{ij} = 1)}{1 - P(y_{ij} = 1)} = b_0 + b_1 x_{1i} + b_2 x_{2ij} =-0.34 + 0.54\,x_{\text{program}} + 0.14\,x_{\text{week}}

    b_0 = −0.34

    b_1 = 0.54 log-odds

    b_2 = 0.14

    x_{\text{program}} is 1 in the program and 0 on the waitlist, and x_{\text{week}} counts weeks from the first check-in, 0 to 3. Software finds the three b's by alternating: fit the coefficients, estimate from the residuals how alike each person's answers are (step 3), refit with weights that reflect it, and repeat until nothing changes. We take its numbers. b_1 = 0.54: at any given week, the program's log-odds of a good night are 0.54 higher.

  3. 3. Estimate How Alike One Person's Nights Are

    Someone who sleeps well one week tends to sleep well the next, so a person's four answers aren't four independent pieces of evidence. GEE estimates that overlap from the residuals and uses it to weight each person's answers. This estimate and the coefficients in step 2 are found together, each update feeding the other.

    \hat{\alpha} = \dfrac{\sum_i \sum_{j<k} r_{ij}\, r_{ik}}{\hat{\phi}\,(M - p)}

    = .27

    r_{ij} is how far answer j from person i sits from the model's prediction, in standard units. Multiplying the residuals for every pair of one person's check-ins, 6 pairs each and M = 900 in all, and averaging (scaled by \hat{\phi} = 1.01, with p = 3) gives \hat{\alpha} = .27.

  4. 4. Measure the Uncertainty from How Much People Actually Varied

    The standard error has to reflect that each person's answers travel together. The robust standard error gets it from the residuals, person by person, so it stays right even if the assumed correlation is off.

    SE_{\text{robust}}(b_j) = \sqrt{\left[A^{-1} B A^{-1}\right]_{jj}}

    = 0.223

    SE_{\text{ignoring clusters}} = 0.166

    A holds the information the model thinks it has and B the scatter the 150 people actually showed. Their sandwich gives 0.223 for the program coefficient. Treating all 600 answers as independent would have given 0.166, too small; multiplying it by √(1 + 3 × .27) = 1.35, the inflation for 4 answers correlating .27, brings it back to about 0.223. That shortcut works for a predictor like the program, which is the same for all of a person's answers, not for one that changes within people, like week.

  5. 5. Compare the Program's Coefficient with Its Uncertainty

    A coefficient only means something relative to its wobble. z counts how many robust standard errors it sits from zero.

    z = \dfrac{b_j}{SE_{\text{robust}}} =\dfrac{0.543}{0.223}

    = 2.44

    \chi^2_{\text{Wald}} = 5.93

    b_j is the program's 0.543 log-odds and SE_{\text{robust}} is 0.223. 0.543 ÷ 0.223 = 2.44 standard errors above zero. Squared, that's the Wald χ² of 5.93 that R prints.

  6. 6. Turn the Coefficient into an Odds Ratio

    Exponentiating a log-odds difference gives a ratio of odds, a size you can compare across studies. Converting back to chances shows what it means in good nights.

    \text{OR} = e^{\,b_j} =e^{0.543}

    = 1.72

    \hat{P}_{\text{waitlist, week 1}} = .42

    \hat{P}_{\text{program, week 1}} = .55

    e^{b_j} with the program's 0.543 is 1.72: the program's odds of a good night are 1.72 times the waitlist's. Turning log-odds back into a chance uses 1 / (1 + e^{-\text{log-odds}}): at week 1 the waitlist's log-odds are −0.34, a chance of .42; the program adds 0.54, giving 0.20 and a chance of .55.

  7. 7. Give a Range of Plausible Values for the Odds Ratio

    Build the range on the log-odds scale, where the estimate's uncertainty is symmetric, then convert both ends into odds ratios.

    e^{\,b_j \pm z^{*} SE_{\text{robust}}} =e^{\,0.543 \pm 1.96 \times 0.223}

    = [1.11, 2.66]

    z^{*} = 1.96 holds the middle 95% of a normal curve. 1.96 × 0.223 = 0.437, so the log-odds range is 0.106 to 0.980, and e^{0.106} = 1.11 and e^{0.980} = 2.66.

  8. 8. Ask How Surprising the Result Would Be If the Program Made No Difference

    Finally, ask how often chance alone would put the coefficient at least this many standard errors from zero, in either direction, if the program changed nothing on average.

    p = 2 \times P(Z \ge |z|) =2 \times P(Z \ge 2.44)

    = .015

    P(Z \ge 2.44) is the chance a standard normal value lands at or above 2.44. Doubling it counts both directions: .015, about 1 in 67.

Run It in Software

The code below reads the running example's data from data.csv; save it in the same folder you run the code from.

SPSS (version 30, as documented by IBM)
  1. Transform → Compute Variable: make week0 = week − 1, so the intercept is week 1.
  2. Analyze → Generalized Linear Models → Generalized Estimating Equations.
  3. Repeated tab: move participant into Subject variables and week into Within-subject variables. Set Structure to Exchangeable, and leave Covariance Matrix on Robust Estimator.
  4. Type of Model tab: choose Binary logistic.
  5. Response tab: move good_night into Dependent Variable, click Reference Category and choose First (lowest value), so the model predicts a 1 (a good night).
  6. Predictors tab: move program and week0 into Covariates. Model tab: move both into Model.
  7. Statistics tab: check Include exponential parameter estimates and Working correlation matrix, then click OK.

Or paste this syntax:

COMPUTE week0 = week - 1.
EXECUTE.
GENLIN good_night (REFERENCE=FIRST) WITH program week0
  /MODEL program week0 INTERCEPT=YES
    DISTRIBUTION=BINOMIAL LINK=LOGIT
  /REPEATED SUBJECT=participant WITHINSUBJECT=week SORT=YES
    CORRTYPE=EXCHANGEABLE COVB=ROBUST
  /PRINT CPS DESCRIPTIVES MODELINFO SUMMARY SOLUTION (EXPONENTIATED) WORKINGCORR.

Reading the Output

Read the Parameter Estimates table: first Exp(B) for program, the population-averaged odds ratio (1.72), then its 95% Wald Confidence Interval for Exp(B) (1.11 to 2.66), then B and its Std. Error (robust, since COVB=ROBUST), and last the Wald Chi-Square and Sig. SPSS tests each coefficient with a Wald χ² on 1 df, the square of the page's z (2.44² = 5.93).

Parameter Estimates

Parameter B Std. Error Wald Chi-Square df Sig. Exp(B) 95% CI Lower Upper
(Intercept) -.339 .1821 3.460 1 .063 .713 .499 1.018
program .543 .2229 5.931 1 .015 1.721 1.112 2.664
week0 .137 .0618 4.900 1 .027 1.147 1.016 1.294

The Working Correlation Matrix shows the same value, about .27, in every off-diagonal cell. SPSS's estimator of that correlation may differ from statsmodels' and R's in the third decimal; the coefficients and standard errors change only in the digits past those shown here.

Output laid out from this page's computed results, not captured from a live run; labels and layout may differ slightly in your version.

R 
library(geepack)

sleep_hygiene <- read.csv("data.csv")
sleep_hygiene$week0 <- sleep_hygiene$week - 1   # 0 = first check-in

fit <- geeglm(good_night ~ program + week0, id = participant,
              data = sleep_hygiene, family = binomial,
              corstr = "exchangeable")
summary(fit)

# Odds ratios with 95% Wald CIs from the robust SEs
est <- coef(summary(fit))
round(exp(cbind(OR = est$Estimate,
                lower = est$Estimate - 1.96 * est$Std.err,
                upper = est$Estimate + 1.96 * est$Std.err)), 3)

Call:
geeglm(formula = good_night ~ program + week0, family = binomial,
    data = sleep_hygiene, id = participant, corstr = "exchangeable")

 Coefficients:
            Estimate Std.err  Wald Pr(>|W|)
(Intercept)  -0.3388  0.1821 3.460   0.0629 .
program       0.5429  0.2229 5.931   0.0149 *
week0         0.1368  0.0618 4.900   0.0269 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation structure = exchangeable
Estimated Scale Parameters:

            Estimate Std.err
(Intercept)        1 0.01871
  Link = identity

Estimated Correlation Parameters:
      Estimate Std.err
alpha   0.2708 0.04363
Number of clusters:   150  Maximum cluster size: 4

        OR lower upper
[1,] 0.713 0.499 1.018
[2,] 1.721 1.112 2.664
[3,] 1.147 1.016 1.294

Read the program row of the odds ratios first (1.72, 95% CI [1.11, 2.66]), then the coefficient and its robust Std.err, then Wald and Pr(>|W|). geepack's Wald is χ² with 1 df, the square of z. Its α̂ is .2708 rather than this page's .2703 because geepack's moment estimates leave out a small degrees-of-freedom correction; nothing that's displayed changes. The data must be sorted by participant, as data.csv is.

Output from a real run: R 4.6.1 (geepack 1.3.13), run 2026-09-29.

Python 
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

sleep_hygiene = pd.read_csv("data.csv")
sleep_hygiene["week0"] = sleep_hygiene["week"] - 1

model = smf.gee("good_night ~ program + week0", groups="participant",
                data=sleep_hygiene, family=sm.families.Binomial(),
                cov_struct=sm.cov_struct.Exchangeable())
result = model.fit()
print(result.summary())
print(model.cov_struct.summary())
print(np.exp(result.params).round(3))
print(np.exp(result.conf_int()).round(3))

The coefficient table from result.summary() (its header block, with the number of clusters, 150, and the working correlation, is left out here), then the other three prints:

==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -0.3388      0.182     -1.860      0.063      -0.696       0.018
program        0.5429      0.223      2.435      0.015       0.106       0.980
week0          0.1368      0.062      2.214      0.027       0.016       0.258
==============================================================================
The correlation between two observations in the same cluster is 0.270
Intercept    0.713
program      1.721
week0        1.147
dtype: float64
               0      1
Intercept  0.499  1.018
program    1.112  2.664
week0      1.016  1.294

The summary's intervals are for the coefficients (log-odds); exponentiate them for the odds ratio's interval, 1.11 to 2.66. std err is the robust (sandwich) standard error by default. statsmodels' α̂ (.270) is the estimator this page uses. For a small number of clusters, model.fit(cov_type="bias_reduced") gives Mancl and DeRouen's corrected standard errors.

Output from a real run: Python 3.9.6 (statsmodels 0.14.6, pandas 2.3.3), run 2026-09-29.

Impact

Effect Size

The population-averaged odds ratio says how many times higher the odds of the outcome are, across the whole population, for a one-unit rise in a predictor. It comes straight from the coefficient, on the same scale as the odds ratios logistic regression and 2 × 2 tables report, so it can be set beside theirs.

\text{OR} = e^{\,b_{\text{program}}} = e^{0.543} = 1.72

Across all 150 adults, the odds of a good night were 1.72 times as high in the program as on the waitlist, at any given week: about 55% against 42% at week 1. That's small by convention. The convention was built for a rare outcome, though: at a good night's base rate, an odds ratio of 1.72 matches a Cohen's d of about 0.34, still small but well above the d = 0.2 where “small” starts.

Report the percentages next to it. An odds ratio compares odds, not chances, and when an outcome is as common as a good night, the odds ratio sits further from 1 than the ratio of the two percentages.

The odds ratio also depends on which question the model asks. A mixed model with a random intercept for each person estimates the conditional odds ratio, the change for a given person, and for yes/no outcomes that's further from 1 than the population-averaged one (Zeger et al., 1988).

Fitted to the same data, a random-intercept logistic model puts the program's odds ratio at 2.09, against GEE's 1.72. Neither is wrong: 2.09 describes a given person's odds, 1.72 the population's.

Benchmarks are conventions, and the odds-ratio benchmarks were worked out for separate groups and a rare outcome (Chen et al., 2010), so treat them as a rough guide. Whether 13 more good nights in 100 is worth the program is a question about sleep, not statistics.

Benchmarks Population-Averaged Odds Ratio at a Glance
Below Small
0.60 < OR < 1.68
Small
1.68 ≤ OR < 3.47 (0.29 < OR ≤ 0.60)
Medium
3.47 ≤ OR < 6.71 (0.15 < OR ≤ 0.29)
Large
OR ≥ 6.71 (OR ≤ 0.15)

Full benchmarks and caveats ↑

Statistical Significance

p answers one narrow question: if the predictor's population-averaged coefficient were zero, how often would its coefficient land at least this many robust standard errors from zero, in either direction? Below the conventional significance level, α = .05, the result is called statistically significant. It leans on the robust standard error, and so on having enough independent clusters (Mancl & DeRouen, 2001).

If the program made no difference, a coefficient at least 2.44 robust standard errors from zero, in either direction, would turn up about 1.5% of the time, roughly 1 in 67. That's below .05, so it's significant.

Ignoring the clustering makes p far too small for a predictor that differs between people, like the program. Treated as 600 independent answers, the same data give the program a standard error of 0.166 instead of 0.223, and p = .001 instead of .015. For a predictor that changes within people, like week, it can go the other way: here week's p would rise from .027 to .065. The choice of working correlation, on the other hand, barely matters here: independence, exchangeable and AR-1 give p = .015, .015 and .013 (Table 3). That's the robust standard error doing its job: as long as the model for the average is right and there are enough clusters, it stays valid even when the working correlation is wrong (Halekoh et al., 2006), so a better working correlation buys precision, not validity.

p is not the probability that the null hypothesis is true, and it says nothing about how large or important the program's effect is. That's what the odds ratio and its interval are for.

Benchmarks p at a Glance
Very Strong
p < .001
Strong
.001 ≤ p < .01
Significant
.01 ≤ p < .05
Weak
.05 ≤ p < .10
Little
p ≥ .10

Strength of evidence against the null hypothesis.

Full benchmarks and caveats ↑

Power and Sample Size

Power is the chance that a study detects an effect that's really there. A real effect doesn't guarantee a significant result: with too few people, it can easily come out non-significant. For GEE, power depends on the size of the population-averaged effect, the significance level (α, usually .05; not the working correlation), how common the outcome is, the number of clusters, the number of observations in each, and how strongly those observations correlate. Extra check-ins help less than they seem to, because correlated answers repeat some of the same information: with 4 answers correlating .27, each person tells us about as much about the program as 4 ÷ (1 + 3 × .27) ≈ 2.2 independent answers would.

Suppose, before collecting any data, we expected effects like the ones we found. We simulated 2,000 studies at each size from the random-intercept logistic model fitted to our data (log-odds −0.46 on the waitlist at week 1, +0.74 for the program, +0.19 per week, and a person SD of 1.41, which reproduces GEE's population-averaged program effect of about 0.54 and the within-person correlation), with half the people in each group, four check-ins and no missing answers, and analyzed each with the same exchangeable GEE. With 150 adults, the program came out significant about 70% of the time; with 200, about 82%. Reaching 80% takes roughly 190 to 200 people.

Researchers usually aim for at least 80% power. Plan it before collecting data, by simulation or with software for clustered designs, using the smallest population-averaged effect that would matter and a realistic within-person correlation. This site doesn't have a calculator for GEE yet.

Other Info

Follow-Up and Post Hoc Tests

None needed. This test answers a single question, so there are no planned contrasts or post hoc tests to run.

References

  • Liang, K.-Y., & Zeger, S. L. (1986). Longitudinal data analysis using generalized linear models. Biometrika, 73(1), 13–22. https://doi.org/10.1093/biomet/73.1.13
    “We introduce a class of estimating equations that give consistent estimates of the regression parameters and of their variance under mild assumptions about the time dependence.” (Summary)
  • Zeger, S. L., Liang, K.-Y., & Albert, P. S. (1988). Models for longitudinal data: A generalized estimating equation approach. Biometrics, 44(4), 1049–1060. https://doi.org/10.2307/2531734
    “population-averaged (PA) models in which the aggregate response for the population is the focus” (Summary)
  • Halekoh, U., Højsgaard, S., & Yan, J. (2006). The R package geepack for generalized estimating equations. Journal of Statistical Software, 15(2), 1–11. https://doi.org/10.18637/jss.v015.i02
    “only correct specification of marginal means is needed for the parameter estimator to be consistent and asymptotically normal” (p. 2)
    “The estimate Σ̂ is a consistent estimate of Σ even if the working correlation matrices R_i(α) are misspecified.” (p. 5)
    “For a small number of clusters (K ≤ 30) the sandwich variance estimator exhibits bias and Paik (1988) recommended using the jackknife variance estimator.” (p. 6)
    “In case of missing values, the GEE estimates are consistent if the values are missing completely at random (Rubin 1976).” (p. 6)
  • McNeish, D., Stapleton, L. M., & Silverman, R. D. (2017). On the unnecessary ubiquity of hierarchical linear modeling. Psychological Methods, 22(1), 114–140. https://doi.org/10.1037/met0000078
    “These alternative methods do not model random effects and thus make a smaller number of assumptions and are interpreted identically to single-level methods with the benefit that estimates are adjusted to reflect clustering of observations.” (Abstract)
  • Mancl, L. A., & DeRouen, T. A. (2001). A covariance estimator for GEE with improved small-sample properties. Biometrics, 57(1), 126–134. https://doi.org/10.1111/j.0006-341X.2001.00126.x
    “Hypothesis tests using the robust covariance estimator can have inflated size when the number of independent clusters is small.” (Summary)
    “The bias-corrected covariance estimator gave tests with sizes close to the nominal level even when the number of subjects was 10 and cluster sizes were unequal, whereas the robust and jackknife covariance estimators gave tests with sizes that could be 2–3 times the nominal level.” (Summary)
  • Pan, W. (2001). Akaike's information criterion in generalized estimating equations. Biometrics, 57(1), 120–125. https://doi.org/10.1111/j.0006-341X.2001.00120.x
    “We propose a modification to AIC, where the likelihood is replaced by the quasi-likelihood and a proper adjustment is made for the penalty term.” (Summary)
  • Chen, H., Cohen, P., & Chen, S. (2010). How big is a big odds ratio? Interpreting the magnitudes of odds ratios in epidemiological studies. Communications in Statistics—Simulation and Computation, 39(4), 860–864. https://doi.org/10.1080/03610911003650383
    “OR = 1.68, 3.47, and 6.71 are equivalent to Cohen's d = 0.2 (small), 0.5 (medium), and 0.8 (large), respectively, when disease rate is 1% in the nonexposed group” (Abstract)
  • Muff, S., Nilsen, E. B., O'Hara, R. B., & Nater, C. R. (2022). Rewriting results sections in the language of evidence. Trends in Ecology & Evolution, 37(3), 203–210. https://doi.org/10.1016/j.tree.2021.10.009
    “we suggest language of evidence that allows for a more nuanced approach to communicate scientific findings as a simple and intuitive alternative to statistical significance testing” (Abstract)
  • Box, G. E. P., & Tidwell, P. W. (1962). Transformation of the independent variables. Technometrics, 4(4), 531–550. https://doi.org/10.1080/00401706.1962.10490038
  • Hosmer, D. W., Lemeshow, S., & Sturdivant, R. X. (2013). Applied logistic regression (3rd ed.). Wiley. https://doi.org/10.1002/9781118548387
  • Maas, C. J. M., & Hox, J. J. (2005). Sufficient sample sizes for multilevel modeling. Methodology, 1(3), 86–92. https://doi.org/10.1027/1614-2241.1.3.86
    “only a small sample size at level two (meaning a sample of 50 or less) leads to biased estimates of the second-level standard errors” (Abstract)
  • Moineddin, R., Matheson, F. I., & Glazier, R. H. (2007). A simulation study of sample size for multilevel logistic regression models. BMC Medical Research Methodology, 7, 34. https://doi.org/10.1186/1471-2288-7-34
    “The estimates of the variance covariance components are slightly biased even with 100 groups and group size of 50.” (Abstract)

Spot an error? Something still confusing? Tell us at support@statsforpsychology.com.