Guides / mixed effects models
Mixed effects models
A mixed effects model lets each group in your data have its own baseline, and sometimes its own slope, while still estimating one overall effect. It is the standard answer for repeated measures and nested data. This guide covers random intercepts and slopes, what the variance components mean, how many clusters you need, why models fail to converge, and when GEE answers your question better.
Written and reviewed by Sidra Khan, Lead BiostatisticianA regression model with two kinds of term in it. Fixed effects are the things you want to estimate for everyone, such as whether a treatment works. Random effects allow each person, ward or school to sit at its own level around that overall picture. It is called mixed because it contains both. The same method is also called a multilevel model, a hierarchical model or a random effects model, and the four names mean the same thing.
The problem
Why you need one
Ordinary regression assumes every row carries independent information. That fails as soon as the same person appears twice, or two pupils share a teacher. When it fails the coefficients usually look sensible but the uncertainty around them is understated, which makes weak findings look strong. How badly depends on how much your observations repeat each other, covered in the guide to clustered data and the unit of analysis.
There are simpler repairs. You can average within each group and analyse the averages, or keep your model and apply cluster-robust standard errors. Both are valid. A mixed effects model earns its extra complexity when you want something those cannot give you: an estimate of how much groups differ, an effect that is allowed to vary between them, or predictions for a particular group.
A mixed model is not automatically the better choice. If you only need honest standard errors and have no interest in the groups themselves, robust errors or aggregation will do the job with far less that can go wrong.
The basic version
Random intercepts
A random intercept gives every group its own starting point. If some wards simply have better outcomes than others before anything is done, a random intercept for ward absorbs that, so the treatment effect is estimated from comparisons within wards rather than being contaminated by which ward a patient happened to be on.
This is the model you should fit first. It handles most clustering, it rarely fails to converge, and it produces the intraclass correlation directly as a by-product. In R it is written as (1 | ward); in Stata as a level in the mixed syntax.
The next step
Random slopes
A random slope lets the effect itself differ between groups. A random intercept says wards start in different places. A random slope says the treatment works better on some wards than others.
Whether you need one is a question about your design, not about model fit. If a predictor varies within groups and you have reason to think its effect differs between them, it belongs in the model. If a predictor is constant within a group, for example the ward's staffing model, it cannot have a random slope at that level.
Random slopes are where mixed models usually break. Every slope you add costs parameters, needs the groups to be numerous enough to estimate variation in that effect, and introduces a correlation between intercept and slope that also has to be estimated. Add them deliberately, one at a time, not as a complete set by default.
The output
Reading the variance components
The part of the output people skip is often the most interesting. Variance components tell you how much of the variation sits between groups and how much within them.
| Quantity | What it tells you | What to do with it |
|---|---|---|
| Between-group variance | How much groups differ from each other once fixed effects are accounted for | Report it. A large value is a finding, not a nuisance. |
| Residual variance | Variation left within groups | Compare with the above to judge where the action is. |
| ICC | Between-group variance as a share of the total | Report it and use it when planning the next study. |
| Slope variance | How much the effect itself differs between groups | If near zero, the random slope is probably not earning its place. |
Two estimation methods appear in every software package. REML gives less biased variance components and is the right default when the variances are what you care about. ML is required if you want to compare models that differ in their fixed effects, because REML likelihoods from different fixed-effect structures are not comparable.
Worked example one
Reading the variance components, and the ICC
A random intercept model of a patient outcome across wards returns this.
| Random effects | Variance | SD |
|---|---|---|
| Ward (intercept) | 4.0 | 2.00 |
| Residual | 16.0 | 4.00 |
Overall mean outcome: 50.0
Step 1. The intraclass correlation
ICC = 4.0 / (4.0 + 16.0) = 0.20
One fifth of the variation lies between wards and four fifths within them. That single number is also what you need to plan the next study, because it drives the design effect.
Step 2. What the between-ward SD means in context
A standard deviation of 2.00 on the ward intercepts means roughly 95% of wards sit within about four points either side of 50, so from 46 to 54. That is a more useful sentence for a reader than "the variance component was 4.0", and it is one line of arithmetic away.
Step 3. Why these numbers are a finding, not a nuisance
If wards vary that much after adjusting for case mix, something about wards matters. The variance components are often the most actionable output of the whole model, and they are routinely omitted from papers that report only the fixed effects.
Worked example two
Shrinkage, calculated
Shrinkage is usually described in words. It is easier to trust once you have seen the arithmetic, because the amount of pull is entirely determined by how much data a group contributed.
The formula
For a group with n observations, the weight given to its own mean is
λ = σ²between / (σ²between + σ²within / n)
The model's estimate for that group is then the overall mean plus λ times the distance from it.
Four wards, same model
| Ward | n | Observed mean | λ | Model estimate |
|---|---|---|---|---|
| A | 4 | 60.0 | 0.500 | 55.0 |
| B | 25 | 60.0 | 0.862 | 58.6 |
| C | 100 | 60.0 | 0.962 | 59.6 |
| D | 4 | 42.0 | 0.500 | 46.0 |
Reading it
Wards A, B and C all observed a mean of 60. The model reports 55.0, 58.6 and 59.6. Ward C contributed enough data to be believed; ward A did not, so it is pulled halfway back towards the overall mean.
Ward D shows the same mechanism in the other direction. An apparently poor ward with only four patients is pulled from 42.0 up to 46.0, because four patients is not enough evidence to conclude it is genuinely worse.
Without shrinkage, the top and bottom of a provider ranking fill up with the smallest units, because small samples produce extreme means by chance. Shrinkage removes most of that artefact, which is why mixed models are the standard approach for comparing schools, wards, surgeons or practices. If you are ranking anything, this is the reason to use one.
A property worth understanding
Shrinkage, and why it is a feature
Group estimates from a mixed model are pulled towards the overall average, and the pull is strongest for groups with little data. A ward with four patients gets dragged close to the average; a ward with four hundred barely moves.
This surprises people who expect the model to report each group's own mean. It is doing something better. A small group's apparent extremeness is mostly noise, and shrinkage prevents you from acting on it. It is why mixed models are preferred for league tables and provider comparisons: without shrinkage the top and bottom of any ranking fill up with the smallest units.
Sample size
How many clusters you need
The number that matters is groups, not rows. Fifty thousand observations from eight clinics is eight pieces of information about how clinics vary.
| Groups | What is realistic |
|---|---|
| Fewer than 10 | Treat group as a fixed effect instead, or aggregate. Variance components will not be trustworthy. |
| 10 to 30 | Random intercepts are usable. Variance estimates are shaky and random slopes usually are not feasible. |
| 30 or more | Random intercepts are comfortable. Random slopes become possible with enough data within each group. |
| 50 or more | Variance components reasonably estimated, and cross-level interactions become realistic. |
Before choosing a model, check how much independent information you have. The free calculator turns your groups, observations each and ICC into a design effect, an effective sample size and the ceiling no amount of extra data can pass.
Open the effective sample size calculatorTroubleshooting
When the model will not converge
Convergence warnings and singular fits are the most common complaint about mixed models, and the cause is nearly always the same: the random structure asks for more than the data can support.
Never respond to a convergence warning by reporting the result anyway. An unconverged model has not estimated anything.
Choosing
Mixed model or GEE?
They answer different questions, and the choice should follow from yours rather than from habit. GEE estimates the population-averaged effect and treats clustering as a nuisance. A mixed model estimates the effect within a group and models the clustering explicitly.
| Choose a mixed model when | Choose GEE when |
|---|---|
| You want the effect within a group | You want the population-averaged effect |
| Variation between groups is itself of interest | Clustering is only a nuisance to correct |
| You have relatively few groups | You have plenty of groups |
| You need predictions for a specific group | You do not |
| Missing data is plausibly at random given what you observed | Missing data is plausibly completely at random |
For continuous outcomes the two usually give similar coefficients. For binary outcomes they do not, and the difference is real: the mixed model odds ratio is larger in magnitude than the GEE one. Reporting one and describing the other misstates the finding.
Doing it
In R, SPSS and Stata
R
lme4 is the standard package. lmer for continuous outcomes, glmer for binary and counts.
library(lme4)
m1 <- lmer(outcome ~ treatment + age + (1 | ward), data = df)
summary(m1)
m2 <- lmer(outcome ~ treatment + age + (1 + treatment | ward), data = df)
anova(m1, m2)
(1 | ward) is a random intercept and (1 + treatment | ward) adds a random slope. performance::icc(m1) returns the ICC directly, and ranef(m1) gives the shrunken group estimates. Note that lme4 deliberately prints no p-values; use lmerTest for Satterthwaite degrees of freedom if a reviewer insists.
SPSS
Linear mixed models are under Analyze, Mixed Models, Linear.
MIXED outcome BY treatment WITH age
/FIXED = treatment age | SSTYPE(3)
/RANDOM = INTERCEPT | SUBJECT(ward) COVTYPE(VC)
/METHOD = REML
/PRINT = SOLUTION TESTCOV.
TESTCOV prints the variance components, which you need for the ICC and which SPSS omits by default. For binary outcomes use GENLINMIXED instead. Add the slope to the RANDOM line as INTERCEPT treatment with COVTYPE(UN) to allow intercept and slope to correlate.
Stata
mixed outcome treatment age || ward:
estat icc
mixed outcome treatment age || ward: treatment, covariance(unstructured)
estat icc computes the intraclass correlation after any mixed fit. For binary outcomes use melogit, and predict, reffects gives the shrunken group estimates. Stata defaults to maximum likelihood, so add , reml when the variance components are what you care about.
Use REML by default, because it gives less biased variance components. Switch to ML when comparing models whose fixed effects differ, since REML likelihoods from different fixed-effect structures are not comparable and any likelihood ratio test between them is meaningless. Comparing random structures is fine under REML.
Avoid these
Six mistakes that cost marks
| Mistake | Why it matters |
|---|---|
| Reporting only the fixed effects | The variance components say how much groups differ, which is often the most useful finding in the model. |
| Fitting random slopes by default | Every slope costs parameters and needs enough groups to estimate variation in that effect. It is the commonest cause of convergence failure. |
| Reporting an unconverged model | It has not estimated anything. Simplify the random structure rather than changing optimiser until a result appears. |
| Comparing fixed effects under REML | REML likelihoods from different fixed structures are not comparable. Refit with ML before any likelihood ratio test on fixed effects. |
| Too few groups | Below ten, variance components are unreliable. Use a fixed effect for group or aggregate instead. |
| Calling a conditional estimate a population effect | For binary outcomes the mixed model odds ratio is larger in magnitude than the marginal one a GEE would report. They answer different questions. |
Writing it up
How to report it
In R this is lme4::lmer and glmer; in Stata mixed and melogit; in SPSS MIXED and GENLINMIXED; in SAS PROC MIXED and GLIMMIX. If the analysis is already written up and needs checking before it is submitted, that is what an independent statistical review covers.
Common questions
Frequently asked questions
What is a mixed effects model?
A regression model containing both fixed effects, which are estimated for everyone, and random effects, which let each group sit at its own level. It is used when observations are grouped, for example repeated measures within people or patients within wards. Multilevel model, hierarchical model and random effects model are other names for the same thing.
What is the difference between a random intercept and a random slope?
A random intercept gives each group its own baseline. A random slope lets the effect of a predictor differ between groups. Start with a random intercept; add slopes only where a predictor varies within groups and you expect its effect to vary between them.
How many groups do I need?
Fewer than ten and you should use fixed effects for group or aggregate instead. Ten to thirty supports random intercepts with shaky variance estimates. Thirty or more is comfortable, and fifty or more is needed before variance components and cross-level interactions are well estimated. The number of groups matters, not the number of rows.
Should I use REML or ML?
REML by default, because it gives less biased variance components. Switch to ML when comparing models that differ in their fixed effects, since REML likelihoods are not comparable across different fixed-effect structures.
My model will not converge. What should I do?
Simplify the random effects first, usually by removing random slopes or the intercept-slope correlation. Centre and scale continuous predictors. Check you have enough groups. Changing optimiser is a last resort, and a result that appears only under one optimiser should not be trusted. Never report an unconverged model.
What does a singular fit mean?
That a variance component has been estimated as essentially zero, usually because the random structure asks for more than the data supports. It is informative rather than an error: the data does not support that random effect, and the model should be simplified.
Why are my group estimates pulled towards the average?
That is shrinkage, and it is deliberate. Groups with little data are pulled hardest because their apparent extremeness is mostly noise. It is why mixed models are preferred for comparing providers or schools, where small units otherwise dominate both ends of any ranking.
Keep reading
Related
Generalised estimating equations
The population-averaged alternative, and how to choose between the two.
Clustered data and the unit of analysis
Why correlated observations shrink your sample, and what the design effect does.
Effective sample size calculator
Design effect, effective sample size and the ceiling you cannot pass. Free.