Multilevel Logistic Regression

The Basics

When to Use It

Use it when your outcome is yes/no, your people come in groups whose members resemble each other (clustered data: patients in clinics, students in classrooms, clients in therapy groups), and your predictors are measured on the people, the groups or both. For example:

Do clinics that send text reminders have more patients come back for their follow-up session?
Does a school's anti-bullying program lower the chance that a student reports being bullied this term, once grade level is taken into account?
Are clients in group CBT less likely to drop out when their group's therapist has more years of experience, given their own baseline anxiety score?

If the question is whether a group-level predictor changes how a person-level one relates to the outcome (does therapist experience weaken the link between baseline anxiety and dropout?), add a cross-level interaction, with a random slope for the person-level predictor.

The model builds on logistic regression; start there if odds, log-odds and odds ratios are new.

Vocabulary You Need

Definition Log-Odds (

\text{logit}(p)

)

The natural log of the odds, ln(p ÷ (1 − p)), also called the logit. A probability is stuck between 0 and 1 and odds between 0 and infinity, but log-odds run from minus to plus infinity, with 0 at a probability of .50: negative below it, positive above it, and symmetric, so .20 and .80 are −1.39 and +1.39. That unlimited scale is why logistic models are built on it: each predictor adds a fixed amount to the log-odds. A coefficient is a change in log-odds; e raised to it is an odds ratio, and 1 ÷ (1 + e^(−log-odds)) turns log-odds back into a probability. In a multinomial model the log-odds compare a category with the baseline category, ln(P(j) ÷ P(base)), rather than with everything else. Glossary →

Definition Random Intercept (

u_{0j}

)

A random effect that lets each group have its own intercept: its own level of the outcome (on the model's scale, such as the log-odds in a logistic model), after the predictors are taken into account. Group j's intercept is the overall intercept plus u₀ⱼ (often written u_j), its deviation, and the model estimates how much those deviations vary (their variance, written τ₀₀ or τ²). A random-intercept model assumes the predictors' slopes are the same in every group, so the groups' regression lines are parallel, shifted up or down. It's the simplest multilevel model; if a predictor's slope varies between groups, that predictor usually also needs a random slope for its standard error to be right. Glossary →

Definition Cluster-Specific Effect

What a multilevel model's coefficients mean: the difference on the outcome's model scale (such as the log-odds) for people who share a cluster, or whose clusters have the same random shift. For a person-level predictor it compares two people in the same cluster (if the predictor is unrelated to the cluster shifts); for a cluster-level one, two clusters with the same random shift, which is not a match on anything else about them. In a linear model it equals the population-averaged effect; with a yes/no outcome it is further from no effect, and the gap grows with the variance between clusters. Also called the conditional or subject-specific effect. Glossary →

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 Latent-Variable ICC (

\rho_{\text{latent}}

)

An intraclass correlation for a yes/no outcome. Picture the outcome as a hidden continuous propensity that tips into “yes” past a threshold. With a logit link that propensity's within-cluster variance is fixed at π²/3 ≈ 3.29, the variance of the standard logistic distribution, so the between-cluster share is τ² ÷ (τ² + π²/3). It describes the propensity, not the observed 0s and 1s, so it isn't comparable with an ICC from a linear model of the same outcome. The latent-variable view suits outcomes that can be thought of as a cutoff on a continuum; for others its meaning is looser. Also called the variance partition coefficient. 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 at a different level. The model allows patients in the same clinic to be alike; it assumes the clinics are independent of each other, and that once a clinic's shift is known, one patient's attendance tells you nothing more about another's. If patients are also grouped by therapist within clinics, or clinics share a manager, the standard errors come out too small for predictors at the level that's been left out.

If It's Violated

Correction Add a random intercept for the level you left out, such as therapists within clinics: (1 | clinic/therapist) in R.
Design Change Recruit clinics that don't share staff, owners or referral sources.

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, or a restricted cubic spline, for each continuous predictor, compared with the straight-line model by a likelihood-ratio test (Hosmer et al., 2013). The Box-Tidwell check adds x × ln(x) instead, for predictors with only positive values. A plot helps more than any p: cut the predictor into bands (such as deciles), and plot the log-odds of the outcome in each band against the band's average; the points should fall near a line.
How Much It Matters: A lot, for what the severity coefficient means. The model gives PHQ-9 one straight-line effect on the log-odds; if attendance holds steady up to moderate symptoms and then drops, one slope blurs the two. Check by grouping severity into bands and plotting each band's log-odds of attending, or by adding a curve and seeing whether it helps.

If It's Violated

Correction Add a squared term or a spline for PHQ-9, centering it first.
Transformation Log-transform a predictor whose effect levels off, such as number of prior sessions. The odds ratio is then per unit of the log.

3. Normally Distributed Random Effects The groups' deviations from the average are roughly bell-shaped.

What It Means: A multilevel model treats each group's intercept (and any random slope) as a draw from a normal distribution with a mean of zero and a variance the model estimates. That's what lets it estimate a variance component from a sample of groups and pull unreliable group estimates toward the average. Groups whose deviations are skewed, bimodal or include a few extreme groups break it.
Applies To: Covariance structure
How to Check: Plot the model's predicted group effects (the empirical Bayes estimates, ranef() in R) in a histogram or Q-Q plot and look for strong skew, two humps or a group far from the rest; then look at what makes any outlying group different. With few groups, the plot says little, so reason about how the groups were sampled too.
How Much It Matters: Usually modest for the fixed effects, more for the variance. The clinic shifts are assumed to follow a normal curve; a few clinics far from the rest pull on τ² and on anything built from it (the ICC, the median odds ratio, the population-averaged effect). With only 30 clinics the shape is hard to judge; plot the predicted clinic shifts and look for one or two far from the others.

If It's Violated

Proceed with Caveat Refit without the one or two most extreme clinics and report whether the conclusions change.
Alternative Test Generalized Estimating Equations (GEE). A GEE estimates the population-averaged effect without assuming any distribution for the clinic shifts.

4. Predictors Unrelated to the Clinic Shifts The clinics' random shifts aren't related to the predictors in the model.

What It Means: The model treats each clinic's shift as unrelated to its patients' predictor values. If clinics with sicker patients also differ in attendance for other reasons, the PHQ-9 slope mixes two relationships: how severity relates to attendance among patients in the same clinic, and how a clinic's average severity relates to its attendance. It then no longer compares two patients in the same clinic.
Applies To: Predictors
How to Check: Add each clinic's mean PHQ-9 as a predictor alongside the patient's own score and see how far its slope is from 0; equivalently, enter PHQ-9 centered within each clinic plus the clinic mean, and compare the within-clinic and between-clinic slopes.
How Much It Matters: Little here, for the severity slope. With the clinic mean added, the within-clinic slope is −0.064, close to the model's −0.061, and the clinic-mean term (the difference between the between-clinic and within-clinic slopes, 0.28 log-odds) can't be told from 0, p = .09. With 30 clinics whose mean PHQ-9 runs only from 10.3 to 14.8, that test has little power, so it's a check, not proof. The reminder odds ratio barely moves (2.14 against 2.11).

If It's Violated

Correction Enter the predictor centered within each cluster plus the cluster mean, and report both slopes. The within-cluster slope then compares people in the same cluster whatever the clusters' shifts.

5. Same Severity Slope in Every Clinic Severity relates to the log-odds of attending about equally in every clinic.

What It Means: A random-intercept model lets clinics differ in their level of attendance but gives every clinic the same PHQ-9 slope, so the clinics' lines are parallel on the log-odds scale. If the slope really varies between clinics and the model leaves that out, the standard error of the average PHQ-9 slope comes out too small.
Applies To: Relationship
How to Check: Add a random slope, (1 + phq9_c | clinic) in R, and compare the two models with a likelihood-ratio test. Keep the slope if it improves the fit, or if theory says slopes should vary and the data can support it; a singular-fit warning (a variance at 0 or a correlation at ±1) says they can't.
How Much It Matters: Little here. A random PHQ-9 slope doesn't improve the fit, χ²(2) = 0.44, p = .80 (Laplace fits), and the fit is singular: the slope's standard deviation is 0.017 log-odds with a correlation of 1.0, so the random-intercept model stands. Reminders vary only between clinics, so they can't have a random slope across clinics.

If It's Violated

Correction Add a random slope for the person-level predictor whose slope varies, and a cross-level interaction if a clinic-level predictor might predict how much it varies.

6. 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 clinics, not the patients. Suggested minimums come mostly from linear multilevel models: with 50 or fewer clusters, the standard errors of the cluster-level variances come out too small, though the fixed effects' are usually fine (Maas & Hox, 2005). For a multilevel logistic model, Moineddin et al. (2007) recommend at least 50 groups of at least 50 people. With fewer, add intervals that don't lean on large-sample normality, such as profile-likelihood or parametric-bootstrap intervals, and check whether the results hold with more quadrature points.
How Much It Matters: A lot, for clinic-level predictors and the variance. Reminders vary only between clinics, so their standard error depends mostly on the 30 clinics, not the 662 patients. With few or small clinics, τ² is estimated imprecisely and with some bias, and its standard error tends to come out too small (Moineddin et al., 2007), and a clinic-level effect has little power. In our 30 clinics, checks that don't lean on the Wald z agree: a likelihood-ratio test of reminders gives χ²(1) = 5.02, p = .025, and the profile-likelihood 95% interval for its odds ratio is [1.11, 4.08] (see Run It in Software).

If It's Violated

Design Change Recruit more clinics rather than more patients per clinic; each extra clinic adds more to the precision of a clinic-level effect.
Correction With small clinics, fit with adaptive quadrature (nAGQ = 10 or more) rather than the Laplace approximation, and compare.

Limitations

This Is Not Your Test If…

Your outcome is a score, not yes/no Multilevel Regression (HLM) →
Your people aren't grouped in clinics, schools or families Logistic Regression →
You want the average effect across the whole population Generalized Estimating Equations (GEE) →
Your outcome has three or more unordered categories Multinomial Regression →
Your outcome is a set of ordered categories Ordinal Regression →
Your outcome is a count, like sessions attended (fit a Poisson GEE, or this model with family = poisson) Generalized Estimating Equations (GEE) →

What It Doesn't Tell You

Why. Clinics chose whether to send reminders, so reminder clinics may differ in other ways (better-organized front desks, less deprived neighborhoods) that raise attendance on their own. The model adjusts for PHQ-9 severity, not for those.
Whether it matters. Significance doesn't make the difference worth the cost. The useful question is whether a 15-point gap in attendance, if reminders caused it, would justify a reminder system, and the interval runs from barely any difference to a large one.
Who else. The results describe clinics like these 30 and their patients. A clinic-level effect estimated from 30 clinics is fragile, and one unusual clinic can move it.
Which effect, for whom. The odds ratio of 2.11 compares clinics with the same random shift, a comparison no one can observe; across the whole population the implied odds ratio is 1.97, and pairs of real clinics vary widely around both (the interval odds ratio runs from 0.58 to 7.65). None of these tells you what reminders would do at any one clinic.

Common Mistakes

Don'tDescribe the multilevel odds ratio as the effect averaged over the whole population.

InsteadSay it compares patients in the same clinic, or clinics with the same random shift. Report a population-averaged odds ratio if that's the question.
Don'tTest the clinic variance with the ordinary χ² p-value.

InsteadHalve it. The null value τ² = 0 sits on the boundary, so the usual p is about twice too large.
Don'tCount patients as your sample size for a clinic-level predictor.

InsteadPlan and interpret its precision by the number of clinics, which limits it far more.
Don'tCompare coefficients or τ² across models with different predictors as if the scale were fixed.

InsteadRemember the latent residual variance is fixed at π²/3, so adding predictors rescales everything else.
Don'tAccept the Laplace default without checking when clusters are small.

InsteadRefit with adaptive quadrature (nAGQ = 10 or more) and report the method you used.
Don'tCompare a latent-variable ICC with an ICC from a linear model or from the observed 0/1 scores.

InsteadSay which scale it's on. The latent ICC fixes the within-cluster variance at π²/3, so it's a different quantity.
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 bigger χ² as a stronger effect.

InsteadReport an effect size suited to the design, such as Cramér's V, an odds ratio or a model's coefficients. χ² grows with sample size even when the effect doesn't.
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 (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.

For an odds ratio below 1, use its reciprocal (1 ÷ OR), so 0.25 counts the same as 4; the bands in brackets are the same cutoffs below 1. Chen, Cohen and Chen (2010) matched these odds ratios to Cohen's d when 1% of the comparison group has the outcome. A more common outcome needs a smaller odds ratio for the same d, so for common outcomes these labels understate the effect. They are rough conventions, not measures of practical or clinical importance.

More about effect size ↓

How It's Reported

Equation Form

Example

OR = 2.11, 95% CI [1.12, 3.94], z = 2.33, p = .020; clinics: τ² = 0.51, latent ICC = .13, MOR = 1.97, LR χ² = 27.40, p < .001

Table 1. Multilevel Logistic Regression Predicting Attendance at the 4-Week Session

Example

	$b$	$SE$	$z$	$p$	OR	95% CI (OR)
Intercept	0.524	0.222	2.36	.018	1.69	[1.09, 2.61]
PHQ-9 (per point)	−0.061	0.019	−3.22	.001	0.94	[0.91, 0.98]
Text reminders	0.744	0.320	2.33	.020	2.11	[1.12, 3.94]

Note. N = 662 patients in 30 clinics. Maximum likelihood, adaptive Gauss-Hermite quadrature (25 points). PHQ-9 is centered at 12, so the intercept's odds are for a patient scoring 12 at a typical clinic without reminders. b is in log-odds; OR = odds ratio, with its 95% Wald interval. Odds ratios are cluster-specific (see Table 3).

Table 2. Variation between Clinics

Example

	$\tau^2$	$\tau$	Latent ICC	MOR	LR $\chi^2$	$p$
Clinic (intercept)	0.51	0.71	.13	1.97	27.40	< .001

Note. τ² is the variance of the clinics' random intercepts (log-odds scale) and τ its standard deviation. The latent-variable ICC is τ² ÷ (τ² + π²/3). MOR = median odds ratio. The likelihood-ratio χ² compares the model with ordinary logistic regression; its p is halved because τ² = 0 lies on the boundary.

Table 3. The Reminder Effect, Three Ways

Example

	OR	$SE$	95% CI
Multilevel model (cluster-specific)	2.11	0.320	[1.12, 3.94]
Population average implied by the multilevel model (PHQ-9 = 12)	1.97
GEE, independence (robust SE)	1.93	0.286	[1.10, 3.38]
Ordinary logistic regression (naive SE)	1.93	0.171	[1.38, 2.70]

Note. The multilevel model's odds ratio compares clinics with the same random shift. The population-averaged odds ratio it implies is for a patient with a PHQ-9 score of 12, averaging over the normal distribution of clinic shifts (it changes a little with the PHQ-9 score chosen; no standard error is computed for it here). GEE with an independence working correlation gives the same estimate as ordinary logistic regression with a cluster-robust standard error; the naive standard error treats the 662 patients as independent. SEs are for the log-odds.

Full Text (APA Style)

Example

Of 662 patients in 30 clinics, 75.2% attended their 4-week follow-up session at the 15 clinics that sent text reminders and 61.2% at the 15 that didn't. In a multilevel logistic regression with a random intercept for clinic, fitted by maximum likelihood with adaptive quadrature, patients at reminder clinics had about twice the odds of attending as patients with the same baseline severity at clinics without reminders that had the same random shift, OR = 2.11, 95% CI [1.12, 3.94]. For a patient with a PHQ-9 score of 12 at a typical clinic, the predicted chance of attending was 78% with reminders and 63% without. The difference was statistically significant, z = 2.33, p = .020. Patients with higher baseline PHQ-9 scores were less likely to attend: each additional point went with about 6% lower odds, OR = 0.94, 95% CI [0.91, 0.98], z = −3.22, p = .001. Clinics differed considerably beyond these predictors, τ² = 0.51 (latent ICC = .13, median odds ratio = 1.97), and a likelihood-ratio test showed the variation was significant, χ² = 27.40, p < .001 (halved for the boundary). Averaged over clinics, the model implies a population-averaged odds ratio of 1.97 for a patient with a PHQ-9 score of 12.

Name the estimation method (Laplace, or adaptive quadrature with how many points) and the software. Different methods give slightly different numbers.
Say that the odds ratios are cluster-specific. If readers need the population-wide effect, report a population-averaged odds ratio too, and say how you got it.
Report the clinic variance (τ²), not just the fixed effects: give τ² or τ, and the latent ICC or median odds ratio to put it in context.
Give a predicted probability or two alongside each important odds ratio, and say which predictor values they're for.

Reading Each Piece

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

p̂ = .63 without reminders, .78 with them	At a typical clinic, about 63 in 100 patients with a PHQ-9 of 12 would be expected to come back without reminders and 78 in 100 with them: 15 percentage points more.
OR = 2.11	Patients at reminder clinics have about twice the odds of attending, comparing clinics with the same random shift. That's odds, not chances: at a typical clinic it's 78% against 63%. By the benchmarks above, which were written for rare outcomes, 2.11 is small and its interval runs from below small to medium; a 15-point gap in attendance is not small for a clinic. Benchmarks ↑
95% CI [1.12, 3.94]	The true clinic-specific odds ratio is plausibly anywhere from barely above 1 (12% higher odds) to almost 4. The data point to higher attendance at reminder clinics but can't pin down how much higher, because only 30 clinics carry the comparison.
τ² = 0.51, latent ICC = .13	About 13% of the variation in the propensity to attend lies between clinics, beyond what reminders and severity predict; in a model with no predictors it's 15%. It's on the latent scale, so it can't be compared with the ICC from a linear model or from the 0s and 1s themselves.
MOR = 1.97	Take two patients with the same PHQ-9 at two randomly chosen clinics (both with or both without reminders): in half of such pairs, the one at the higher-attendance clinic has less than 1.97 times the other's odds, in half more. Clinics differ from each other by about as much as reminder clinics differ from the rest.
OR_PA = 1.97	Across the whole population of clinics, patients at reminder clinics have 1.97 times the odds of attending of patients at the rest. It's a different question from the clinic-specific 2.11, not a correction of it: a GEE would estimate this one.
z = 2.33	The reminder coefficient is 2.33 standard errors above zero. A Wald z is judged against the normal curve, so there are no degrees of freedom to report.
LR χ² = 27.40	Adding a separate shift for each clinic improves the fit far more than chance would. Report it with the halved p, and say so.
p = .020	If patients at reminder and no-reminder clinics with the same random shift had the same odds of attending, a z at least this far from 0, in either direction, would turn up about 2% of the time. That's below .05, so it's statistically significant, if not by a wide margin. Benchmarks ↑
p < .001	If clinics didn't differ at all beyond reminders and severity, a likelihood gain at least this large would turn up far less than 0.1% of the time. The data give strong evidence that clinics differ. Benchmarks ↑

Calculations

The Running Example

Example

The follow-up study: 662 patients in 30 clinics (15 send text reminders), a yes/no outcome (attended the 4-week session) and each patient's baseline PHQ-9 score.

	Null Hypothesis	Alternative
Reminders	$\beta_2 = 0$	$\beta_2 \neq 0$
Clinic Variation	$\tau^2 = 0$	$\tau^2 > 0$

$\beta_2$ is the difference in the log-odds of attending between a clinic that sends reminders and one that doesn't but has the same random shift, for patients with the same PHQ-9 score. $\tau^2$ is the variance of the clinics' random intercepts: how much clinics differ in attendance beyond what reminders and severity predict. The second alternative is one-directional because a variance can't be negative.

What Your Data Should Look Like

One row per person, in long format: a 0/1 outcome column, a column naming the person's cluster (clinic 1 to 30), and a column for each predictor. A cluster-level predictor, like whether the clinic sends reminders, repeats the same value on every row from that clinic. Clusters can have different numbers of people. The box below shows attendance by reminders; the model uses every patient's row.

Example

patient	clinic	reminders	phq9	attended
1	1	0	11	1
2	1	0	14	1
3	1	0	21	1
4	1	0	11	0
5	1	0	16	1
6	1	0	12	1
… 656 more rows

	n	Yes	%
No Reminders	327	200	61.2%
Text Reminders	335	252	75.2%

Simulated data.

Download all 662 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 follow-up study's numbers. Collapse all to read just the headings: they tell the whole story on their own.

1. Write Down How Attendance Depends on Severity, Reminders and the Clinic

Odds are the chance of attending divided by the chance of not: a 75% chance is odds of 3 to 1. Logistic regression predicts the log of the odds, which can take any value, so predictors can add to it in straight lines. The multilevel version adds one term, each clinic's own shift up or down, drawn from a normal distribution whose spread the model estimates too.

$\ln\dfrac{p_{ij}}{1 - p_{ij}} = \beta_0 + \beta_1 x_{ij} + \beta_2 w_j + u_j, \quad u_j \sim N(0,\, \tau^2)$ = $0.524 - 0.061\,(\text{PHQ-9} - 12) + 0.744\,\text{reminders} + u_j$

$b_0$ = 0.524
$b_1$ = −0.061
$b_2$ = 0.744

$p_{ij}$ is patient i's chance of attending at clinic j, $x_{ij}$ their PHQ-9 minus 12, $w_j$ is 1 if clinic j sends reminders and $u_j$ is clinic j's own shift; $\tau^2$ is the variance of the $u_j$ (some books write $u_{0j}$ and $\tau_{00}$ ). The $\beta$ 's are the true values and the $b$ 's our estimates of them. PHQ-9 is centered at 12, near the sample mean, so $b_0$ describes a patient with typical severity; centering at one number for everyone doesn't separate differences within clinics from differences between them (see Predictors Unrelated to the Clinic Shifts). No napkin can find the $b$ 's: software averages each clinic's likelihood over every possible $u_j$ (adaptive quadrature, 25 points per clinic) and picks the values that make the observed attendance most likely.
2. Turn the Reminder Coefficient into an Odds Ratio

A coefficient is a difference in log-odds, which nobody thinks in. Exponentiating it gives how many times higher the odds of attending are at reminder clinics.

$\text{OR} = e^{b}$ = $e^{0.7444}$

= 2.11

$b$ is reminders' 0.7444 log-odds. $e^{0.7444}$ = 2.11: at a clinic that sends reminders, a patient's odds of attending are about twice those of a patient with the same PHQ-9 at a clinic that doesn't but has the same random shift. The same random shift means the two clinics sit at the same point in the spread of clinics; it isn't a match on staffing, neighborhood or anything else about them.
3. Give a Range of Plausible Odds Ratios

Build the range on the log-odds scale, where the estimate is roughly normal, then exponentiate both ends. That's why the interval isn't symmetric around 2.11.

$\left[\,e^{\,b - z^{*} SE},\;\; e^{\,b + z^{*} SE}\,\right]$ = $e^{\,0.744 \pm 1.960 \times 0.320}$

= [1.12, 3.94]

$SE_{b_2}$ = 0.320

$SE$ is 0.320, from the curvature of the likelihood at its peak. 1.960 × 0.320 = 0.628, so the log-odds range is 0.117 to 1.372, and $e^{0.117}$ = 1.12 to $e^{1.372}$ = 3.94.
4. Test Whether the Reminder Coefficient Differs from 0

Divide the coefficient by its standard error and compare the result with the normal curve. With a random intercept, that standard error reflects how much clinics differ, not just how many patients there are.

$z = \dfrac{b}{SE_b}$ = $\dfrac{0.744}{0.320}$

= 2.33

$p$ = .020

$b$ = 0.744 and $SE_b$ = 0.320, so the coefficient is 2.33 standard errors above 0. The chance of a z at least that far from 0, either way, is .020.
5. Measure How Much Clinics Differ

The model estimates the variance of the clinic shifts on the log-odds scale. To read it as a share, compare it with the variance of the underlying propensity within a clinic, which a logit model fixes at π²/3.

$\rho_{\text{latent}} = \dfrac{\tau^2}{\tau^2 + \pi^2/3}$ = $\dfrac{0.506}{0.506 + 3.290}$

= .13

$\tau^2$ = 0.506
$\text{MOR}$ = 1.97
$\tau^2_{\text{Laplace}}$ = 0.50

$\tau^2$ is the clinics' variance, 0.506, and $\pi^2/3$ = 3.290. 0.506 ÷ 3.796 = .13: about 13% of the variation in patients' underlying tendency to attend lies between clinics, beyond what reminders and severity predict. On the odds scale, the median odds ratio $e^{\sqrt{2 \times 0.506} \times 0.6745}$ = 1.97. The Laplace approximation, lme4's default, puts $\tau^2$ at 0.50 instead.
6. Test Whether Clinics Really Differ

Compare how probable the data are with and without the clinic intercepts. A variance can't go below 0, so the usual chi-square p-value is too large; halve it.

$\chi^2 = 2(\ln L_1 - \ln L_0), \quad p = \tfrac{1}{2}\,P(\chi^2_1 \ge \chi^2)$ = $2\,\big(-388.39 - (-402.09)\big)$

= 27.40

$\ln L_1$ = −388.39
$\ln L_0$ = −402.09
$p$ < .001

$\ln L_1$ = −388.39 with the clinic intercepts and $\ln L_0$ = −402.09 without them (ordinary logistic regression). 2 × 13.70 = 27.40. Half the chance of a χ² with 1 df reaching 27.40 is far below .001.
7. Predict Attendance at a Typical Clinic

Odds ratios are hard to feel. Plug values into the equation, set the clinic's shift to 0 (a typical clinic) and turn the log-odds back into a probability.

$\hat{p} = \dfrac{1}{1 + e^{-(b_0 + b_1 x + b_2 w)}}$ = $\dfrac{1}{1 + e^{-(0.524 + 0.744)}}, \quad \dfrac{1}{1 + e^{-0.524}}$

= .63

$\hat{p}_{\text{reminders}}$ = .78

For a patient with PHQ-9 = 12, $b_1 x$ is 0, so the log-odds are $b_0$ = 0.524 without reminders and 0.524 + 0.744 = 1.268 with them. 1 ÷ (1 + $e^{-0.524}$ ) = .63, and 1 ÷ (1 + $e^{-1.268}$ ) = .78.
8. Average over All Clinics to Get the Population-Averaged Effect

The odds ratio so far compares clinics with the same shift. To compare all reminder clinics with all others, average each group's probability over the whole spread of clinics first. Averaging flattens the S-shaped curve, so the odds ratio shrinks toward 1.

$\bar{p}_w = \int \dfrac{1}{1 + e^{-(b_0 + b_2 w + u)}}\, \phi(u;\, 0, \tau^2)\, du, \quad \text{OR}_{PA} = \dfrac{\bar{p}_1 / (1 - \bar{p}_1)}{\bar{p}_0 / (1 - \bar{p}_0)}$ = $\dfrac{.7590 / .2410}{.6155 / .3845}$

= 1.97

$\bar{p}_0$ = .616
$\bar{p}_1$ = .759

Averaged over clinic shifts with $\tau^2$ = 0.506, a PHQ-9 12 patient's chance of attending is $\bar{p}_1$ = .759 across reminder clinics and $\bar{p}_0$ = .616 across the rest. Their odds, 3.15 and 1.60, give an odds ratio of 1.97, below the clinic-specific 2.11.
9. See What Ignoring the Clinics Would Do

Ordinary logistic regression treats the 662 patients as independent. Its odds ratio lands near the population-averaged one, but its standard error is far too small. A cluster-robust standard error, as GEE uses, corrects most of that, though with 30 clinics it still tends to come out a little too small.

$\text{OR}_{\text{ordinary}} = e^{\,b_{\text{ordinary}}}$ = $e^{0.656}$

= 1.93

$b_{\text{ordinary}}$ = 0.656
$SE_{\text{naive}}$ = 0.171
$SE_{\text{robust}}$ = 0.286

Ignoring clinics, reminders' coefficient is 0.656, so $e^{0.656}$ = 1.93, with a standard error of 0.171 log-odds, as if 662 separate patients had each been assigned reminders. Summing the scatter clinic by clinic instead gives 0.286, and the multilevel model's is 0.320: the reminder comparison really rests on 30 clinics.

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)

Set attended to Nominal in Variable View, and compute the centered severity: Transform → Compute Variable, phq9_c = phq9 - 12.
Analyze → Mixed Models → Generalized Linear. On Data Structure, drag clinic into Subjects.
Fields and Effects → Target: choose attended, then Custom: Binomial distribution, Logit link. Under More, set the reference category to 0, so the model predicts attending.
Fixed Effects: add phq9_c and reminders, with the intercept. Random Effects: Add Block, check Include intercept, subject combination clinic.
Build Options: leave the defaults, then Run.

Or paste this syntax:

COMPUTE phq9_c = phq9 - 12.
EXECUTE.
GENLINMIXED
  /DATA_STRUCTURE SUBJECTS=clinic
  /FIELDS TARGET=attended
  /TARGET_OPTIONS REFERENCE=0 DISTRIBUTION=BINOMIAL LINK=LOGIT
  /FIXED EFFECTS=phq9_c reminders USE_INTERCEPT=TRUE
  /RANDOM USE_INTERCEPT=TRUE SUBJECTS=clinic COVARIANCE_TYPE=VARIANCE_COMPONENTS
  /BUILD_OPTIONS TARGET_CATEGORY_ORDER=ASCENDING INPUTS_CATEGORY_ORDER=ASCENDING
    MAX_ITERATIONS=100 CONFIDENCE_LEVEL=95 DF_METHOD=RESIDUAL COVB=MODEL.

Reading the Output

SPSS shows the results in the Model Viewer. Read the Fixed Coefficients table first: each row's Exp(Coefficient) (the odds ratio) and its 95% interval, then Coefficient (log-odds), t and Sig. Then open Random Effect Covariance for the clinic intercept's variance (τ²).

Expect SPSS's numbers to differ a little from R's. GENLINMIXED doesn't maximize the likelihood by quadrature: IBM's algorithms documentation says it estimates GLMMs "using linearization-based methods, also called the pseudo likelihood approach." It also tests coefficients with t on residual degrees of freedom rather than z. Neither method is wrong, but report which you used; for an exact maximum likelihood fit with adaptive quadrature, use R's glmer as below.

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(lme4)
followup <- read.csv("data.csv")
followup$phq9_c <- followup$phq9 - 12   # center severity at 12

fit <- glmer(attended ~ phq9_c + reminders + (1 | clinic), data = followup,
             family = binomial, nAGQ = 25)
summary(fit)

# Odds ratios with 95% Wald intervals
b <- fixef(fit); se <- sqrt(diag(vcov(fit)))
round(exp(cbind(OR = b, lower = b - qnorm(0.975) * se, upper = b + qnorm(0.975) * se)), 3)

# Clinic variation: latent-variable ICC and median odds ratio
tau2 <- VarCorr(fit)$clinic[1]
c(tau2 = tau2, ICC = tau2 / (tau2 + pi^2 / 3), MOR = exp(sqrt(2 * tau2) * qnorm(0.75)))

# Likelihood-ratio test of the clinic variance, p halved for the boundary
fit0 <- glm(attended ~ phq9_c + reminders, data = followup, family = binomial)
lr <- as.numeric(2 * (logLik(fit) - logLik(fit0)))
c(LR = lr, p = pchisq(lr, df = 1, lower.tail = FALSE) / 2)

# The Laplace approximation (glmer's default) for comparison
VarCorr(update(fit, nAGQ = 1))$clinic[1]

# With only 30 clinics, check the reminder effect without the Wald z:
# a likelihood-ratio test and a profile-likelihood interval for its odds ratio
no_reminders <- update(fit, . ~ . - reminders)
anova(no_reminders, fit)
round(exp(confint(fit, parm = "reminders", method = "profile")), 3)

Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 25) [glmerMod]
 Family: binomial  ( logit )
Formula: attended ~ phq9_c + reminders + (1 | clinic)
   Data: followup

      AIC       BIC    logLik -2*log(L)  df.resid
    784.8     802.8    -388.4     776.8       658

Scaled residuals:
    Min      1Q  Median      3Q     Max
-2.9236 -0.9355  0.4673  0.6528  1.5545

Random effects:
 Groups Name        Variance Std.Dev.
 clinic (Intercept) 0.5065   0.7117
Number of obs: 662, groups:  clinic, 30

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.52357    0.22226   2.356  0.01849 *
phq9_c      -0.06082    0.01890  -3.218  0.00129 **
reminders    0.74446    0.32019   2.325  0.02007 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
          (Intr) phq9_c
phq9_c    -0.023
reminders -0.686 -0.022

               OR lower upper
(Intercept) 1.688 1.092 2.610
phq9_c      0.941 0.907 0.976
reminders   2.105 1.124 3.943

     tau2       ICC       MOR
0.5064702 0.1334102 1.9715890

          LR            p
2.740184e+01 8.263677e-08

[1] 0.4961049

Data: followup
Models:
no_reminders: attended ~ phq9_c + (1 | clinic)
fit: attended ~ phq9_c + reminders + (1 | clinic)
             npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)
no_reminders    3 787.80 801.29 -390.90    781.80
fit             4 784.78 802.76 -388.39    776.78 5.0213  1    0.02504 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Computing profile confidence intervals ...
          2.5 % 97.5 %
reminders 1.105  4.082

Read the odds ratios and their intervals first, then each row of Fixed effects (Estimate in log-odds, z value, Pr(>|z|)), then the clinic's Variance under Random effects with the ICC and MOR computed from it. The first line names the method: without nAGQ, glmer uses the Laplace approximation, which here gives a clinic variance of 0.496 instead of 0.506. The last digits of τ² differ slightly from the page's (0.5064702 against 0.5064658) because the optimizer stops at a tolerance; every value agrees to the precision shown on the page. Use logLik(fit) for the likelihood-ratio test: deviance() on a glmer fit returns a different quantity. The last two blocks check the reminder effect without the Wald z, which matters with only 30 clinics: the likelihood-ratio test gives χ²(1) = 5.02, p = .025, and the profile-likelihood interval for the odds ratio, [1.11, 4.08], is a little wider at the top than the Wald [1.12, 3.94]. Both agree with the Wald result.

If glmer warns of a singular (boundary) fit, a variance is estimated at or near 0, or a correlation at ±1: the data can't support that part of the random structure, so simplify it. If it warns that the model failed to converge, check with allFit(fit), which refits with every available optimizer; if they agree, the warning can usually be set aside.

Output from a real run: R 4.6.1 (lme4 2.0-6), run 2026-10-03.

Python

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

followup = pd.read_csv("data.csv")
followup["phq9_c"] = followup["phq9"] - 12

# statsmodels has no maximum likelihood fit for a multilevel logistic model
# (BinomialBayesMixedGLM is a variational Bayes approximation). A GEE gives the
# population-averaged comparison with cluster-robust standard errors.
gee = smf.gee("attended ~ phq9_c + reminders", groups="clinic", data=followup,
              family=sm.families.Binomial(), cov_struct=sm.cov_struct.Independence()).fit()
print(gee.summary())
print(np.exp(gee.params).round(3))

                               GEE Regression Results
===================================================================================
Dep. Variable:                    attended   No. Observations:                  662
Model:                                 GEE   No. clusters:                       30
Method:                        Generalized   Min. cluster size:                  14
                      Estimating Equations   Max. cluster size:                  28
Family:                           Binomial   Mean cluster size:                22.1
Dependence structure:         Independence   Num. iterations:                     2
Date:                     Tue, 29 Sep 2026   Scale:                           1.000
Covariance type:                    robust   Time:                         13:36:32
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      0.4583      0.228      2.008      0.045       0.011       0.906
phq9_c        -0.0489      0.017     -2.936      0.003      -0.081      -0.016
reminders      0.6563      0.286      2.296      0.022       0.096       1.217
==============================================================================
Skew:                         -0.7417   Kurtosis:                      -1.2857
Centered skew:                -0.6099   Centered kurtosis:             -0.9036
==============================================================================
Intercept    1.581
phq9_c       0.952
reminders    1.928
dtype: float64

These are population-averaged estimates, the third row of Table 3: reminders' odds ratio is 1.93 with a robust standard error of 0.286. For the multilevel model itself, fit glmer in R rather than reading BinomialBayesMixedGLM's approximate estimates as maximum likelihood ones.

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

Impact

Effect Size

A multilevel logistic regression has effect sizes at two levels. For a predictor, the odds ratio says how many times higher the odds of the outcome are per unit, comparing people in the same cluster or clusters with the same shift. For the clusters themselves, the latent-variable ICC gives the share of the variation between clusters, and the median odds ratio puts that variation on the odds ratio scale.

\text{OR} = e^{b_2} = e^{0.744} = 2.11, \qquad \text{MOR} = e^{\sqrt{2 \times 0.506} \times 0.6745} = 1.97

Patients at reminder clinics have about twice the odds of attending, comparing clinics with the same random shift. Two clinics picked at random typically differ by about as much: reminder clinics differ from the rest by about as much as clinics differ from each other.

Odds ratios overstate how much the chances differ when an outcome is common, as attending is here. At a typical clinic, the predicted chance for a patient with a PHQ-9 of 12 is 63% without reminders and 78% with them, 15 percentage points apart; give numbers like these alongside the odds ratio.

For a clinic-level predictor like reminders, the odds ratio is harder to read than it looks. It compares two clinics with the same random shift, a comparison no one can observe, because each clinic has only one shift and either sends reminders or doesn't. Austin and Merlo (2017) suggest adding the interval odds ratio: the middle 80% of the odds ratios you'd get comparing a patient at a random reminder clinic with an identical patient at a random no-reminder clinic. Here it runs from 0.58 to 7.65. Because that range includes 1, the variation between clinics is large next to the reminder difference: in about 23% of such pairs, the patient at the no-reminder clinic has the higher odds.

The odds ratio also depends on the question. The multilevel model's 2.11 is cluster-specific; averaged over the whole spread of clinics, the same model implies a population-averaged odds ratio of 1.97 for a patient with a PHQ-9 score of 12, and a GEE estimates 1.93. For a yes/no outcome the population-averaged effect is closer to 1 whenever clusters vary, and the gap grows with the variance between clusters (Austin & Merlo, 2017).

The odds ratio benchmarks above come from Chen, Cohen and Chen (2010) and assume a rare outcome; by them 2.11 is small, and its interval, 1.12 to 3.94, runs from below small to medium. Benchmarks are conventions. Whether a 15-point gap in attendance would matter, if reminders caused it, is a question about clinical care, not statistics.

Benchmarks 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

A multilevel logistic regression answers two kinds of significance question. Each coefficient's Wald z asks: if this predictor were unrelated to the odds (with the other predictors and the cluster's shift held constant), how often would its estimate land at least this far from zero, in either direction? The likelihood-ratio test of the random intercept asks: if clusters didn't differ at all, how often would adding a shift for each cluster improve the fit this much? Because a variance can't be negative, that test's ordinary chi-square p is about twice too large and is halved (likelihood-ratio test; Self & Liang, 1987).

If patients at reminder and no-reminder clinics with the same random shift had the same odds of attending, a z at least this far from 0, in either direction, would turn up about 2% of the time (p = .020). If clinics didn't differ at all, a likelihood gain of 27.40 or more would turn up far less than 0.1% of the time (p < .001). Ignoring the clinics would have made the reminder difference look far more certain than it is: p = .0001, because the naive standard error of 0.171 treats 662 patients as 662 independent tests of reminders.

The p-values depend on how the likelihood was approximated. Here the Laplace approximation (lme4's default) and 25-point adaptive quadrature agree closely, with 14 to 28 patients per clinic; with small clusters or rare outcomes they can differ, so use quadrature and say which you used.

p is not the probability that reminder and no-reminder clinics have the same odds, and it says nothing about how large the difference is. With only 30 clinics, the reminder test has modest power, so a non-significant result would have said little. That's what the odds ratio, its interval and the predicted probabilities 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 clusters, it can easily come out non-significant. In a multilevel logistic regression, power depends on the size of the effect, how common the outcome is, how much clusters differ (τ²), how many people each cluster has, and, above all for a cluster-level predictor, how many clusters there are.

There is no simple formula, so power is usually worked out by simulation: generate many data sets from the model you expect, fit the model to each, and count how often the effect comes out significant.

Suppose the truth were exactly what we estimated: an odds ratio of 2.11 for reminders, a PHQ-9 slope of −0.061 and τ² = 0.51. We simulated 2,000 studies for each design, with PHQ-9 scores drawn from a normal curve (mean 12, SD 5, rounded and kept within 0–27), reminders at half the clinics and every clinic the same size. Each study was fitted with glmer (Laplace) and the reminder effect tested with its Wald z. With 30 clinics of 22 patients, the reminder effect came out significant about 68% of the time. Doubling the patients per clinic to 44 helped only a little (76%); doubling the clinics to 60 raised it to 92%. With 2,000 studies, each figure is good to about ±2 percentage points.

The extra patients add little because reminders vary only between clinics: once each clinic's attendance rate is known reasonably well, more patients from the same clinics don't add new comparisons. Estimates need enough clusters too. In Moineddin, Matheson and Glazier's (2007) simulations, variance components were somewhat biased and their standard errors too small even with many groups, and they recommend at least 50 groups of at least 50 people for valid estimates. In their results the fixed effects' standard errors were unbiased while the variance components' were too small, so with our 30 clinics of 14 to 28 patients, read τ², the ICC and the median odds ratio as approximate in particular.

Researchers usually aim for at least 80% power. Plan it before collecting data, by simulation, for the smallest effect that would matter, and add clusters rather than people when the key predictor is measured at the cluster level. There's no calculator for this design on the site; a short R loop that simulates data sets and refits glmer to each runs the simulation.

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.

Related Methods

Related DesignLogistic Regression: logistic regression is this model with the clinic variance fixed at 0. Add a random intercept and the coefficients become cluster-specific, usually a little further from 0.
Related DesignMultilevel Regression (HLM): the same random-intercept model for a continuous outcome. There the within-cluster and population-averaged effects are the same number; with a logit link they aren't.
Related DesignGeneralized Estimating Equations (GEE): a GEE answers the population-averaged question for the same clustered yes/no data, with robust standard errors and no model for the clinic shifts.
Related DesignCochran's Q: Cochran's Q compares a yes/no outcome across three or more occasions in the same people. It uses only people measured every time and takes no covariates; a multilevel logistic model keeps incomplete cases and describes a typical person.
Related DesignMarkov Models: For repeated categorical states, a multilevel model of the next state, with random effects for people, is the regression route to a mixed Markov model in which each person has their own transition chances.

References

Wong, G. Y., & Mason, W. M. (1985). The hierarchical logistic regression model for multilevel analysis. Journal of the American Statistical Association, 80(391), 513–524. https://doi.org/10.1080/01621459.1985.10478148
Austin, P. C., & Merlo, J. (2017). Intermediate and advanced topics in multilevel logistic regression analysis. Statistics in Medicine, 36(20), 3257–3277. https://doi.org/10.1002/sim.7336
“the estimated covariates from a population‐average model are closer to zero than are those from the cluster‐specific or conditional model” (Section 4.3, Marginal or population-average regression coefficients)

“When outcomes are binary, the subject‐level variance is fixed when using the latent variable formulation (being π²/3 when using the logit link and 1 when using the probit link).” (Section 4.6, Proportional change in cluster variation)

“The variance of a logistic distribution with scale parameter equal to one is π²/3.” (Section 4.5.1, The variance partition coefficient)

“If one were to repeatedly sample at random two subjects with the same covariates from different clusters, then the MOR is the median odds ratio between the subject at higher risk of the outcome and the subject at the lower risk of the outcome” (Section 4.5.2, Median odds ratio)

“An advantage to use of the MOR for quantifying the contextual effect is that it is on the same scale as that used for estimating measures of association” (Section 4.5.2, Median odds ratio)
Goldstein, H., Browne, W., & Rasbash, J. (2002). Partitioning variation in multilevel models. Understanding Statistics, 1(4), 223–231. https://doi.org/10.1207/S15328031US0104_02
Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01
Bates, D., Maechler, M., Bolker, B., & Walker, S. (2026). lme4: Linear mixed-effects models using Eigen and S4 (R package version 2.0-6) [Computer software]. https://CRAN.R-project.org/package=lme4
“Defaults to 1, corresponding to the Laplace approximation. Values greater than 1 produce greater accuracy in the evaluation of the log-likelihood at the expense of speed.” (glmer help page, argument nAGQ)
Bolker, B., et al. (2025, July 19). GLMM FAQ. https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html
“in the simplest case (single random effect variance), the p-value is approximately twice as large as it should be (Pinheiro and Bates 2000)” (Testing hypotheses)
Self, S. G., & Liang, K.-Y. (1987). Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. Journal of the American Statistical Association, 82(398), 605–610. https://doi.org/10.1080/01621459.1987.10478472
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
“We recommend using a minimum group size of 50 with at least 50 groups to produce valid estimates for multi-level logistic regression models.” (Abstract)

“The estimates of the variance covariance components are slightly biased even with 100 groups and group size of 50.” (Abstract)
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)
IBM Corp. (2022). IBM SPSS Statistics 29 algorithms. IBM Corp. https://www.ibm.com/docs/en/SSLVMB_29.0.0/pdf/IBM_SPSS_Statistics_Algorithms.pdf
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)
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
Schielzeth, H., Dingemanse, N. J., Nakagawa, S., Westneat, D. F., Allegue, H., Teplitsky, C., Réale, D., Dochtermann, N. A., Garamszegi, L. Z., & Araya-Ajoy, Y. G. (2020). Robustness of linear mixed-effects models to violations of distributional assumptions. Methods in Ecology and Evolution, 11(9), 1141–1152. https://doi.org/10.1111/2041-210X.13434
“Model estimates were usually robust to violations of assumptions, with the exception of slight upward biases in estimates of random effect variance if the generating distribution was bimodal but was modelled by Gaussian error distributions.” (Abstract)
MacCallum, R. C., Roznowski, M., & Necowitz, L. B. (1992). Model modifications in covariance structure analysis: The problem of capitalization on chance. Psychological Bulletin, 111(3), 490–504. https://doi.org/10.1037/0033-2909.111.3.490
“The use of alternative a priori models is recommended as a preferred strategy.” (Abstract)
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
“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)
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
“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)