Mixed Effects Models in R

June 14, 2026 · 10 min read ·Updated June 14, 2026 ·intermediate

rlme4nlmemixed-effectsstatisticsregressionglmer

Mixed effects models: when to use them

Mixed effects models in R solve a specific problem: observations are not independent. Subjects are measured repeatedly, students sit in classrooms, batches come from the same machine, and rows of the data frame are correlated in ways ordinary least squares cannot see. The standard errors that lm() produces under those conditions are too small, and the p-values are too good.

The conceptual payoff is partial pooling. A complete-pooling model fits one grand mean and ignores the groups; a no-pooling model fits a separate mean per group and ignores the population. A mixed model sits between them: each group gets its own mean, but those group means are shrunk toward the population mean in proportion to how much information each group carries. That shrinkage is what makes a model with 10 observations per subject behave sensibly when you have 18 subjects instead of 18,000.

The two packages that do this in R are lme4 and nlme. They overlap a lot, but they make different bets about what models are easy to express.

What this guide covers

lme4 vs nlme
The lmer formula syntax
Random intercepts, random slopes, nested and crossed designs
Generalized linear mixed models with glmer
Structured residuals with nlme::lme
REML vs ML, singular fits, and convergence
Diagnostics and reporting

lme4 vs nlme: which one to pick

Need	Pick
Speed on large or sparse designs	`lme4`
Generalized linear mixed models (binomial, Poisson, Gamma)	`lme4` (`glmer`)
Nonlinear mixed models	`lme4` (`nlmer`)
Correlated within-group errors (AR(1), compound symmetry)	`nlme` (`lme` with `corAR1()`)
Heteroscedastic variance functions	`nlme`
P-values with Satterthwaite df out of the box	`lmerTest` (wraps `lme4`)
Bayesian mixed model	`brms`, `rstanarm`, `MCMCglmm`

If you don’t need correlated residuals or variance functions, default to lme4. It is faster, scales to larger data, and has the cleaner glmer interface.

The lmer formula syntax

Random effects go in parentheses with the form (expr | grouping). The expr is the random-effects design matrix and the grouping is the cluster variable.

Formula	Meaning
`(1 \| g)`	Random intercept per level of `g`
`(x \| g)`	Random intercept and random slope of `x`, correlation estimated
`(x \|\| g)`	Random intercept and random slope, correlation fixed to 0
`(0 + x \| g)`	Random slope of `x` only, no random intercept
`(1 \| g1) + (1 \| g2)`	Crossed random effects (every level of `g1` sees every level of `g2`)
`(1 \| g1/g2)`	`g2` nested within `g1`

A worked example on the built-in sleepstudy data, which has 180 reaction-time measurements from 18 subjects across 10 days of sleep deprivation:

library(lme4)
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
summary(fm1)
# Linear mixed model fit by REML. t-tests use Satterthwaite's method [lmerTest]
# ...
# Fixed effects:
#             Estimate Std. Error t value
# (Intercept)  251.405      6.632   37.91
# Days          10.467      1.502    6.97
#
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr
#  Subject  (Intercept) 612.09   24.74
#           Days         35.07    5.92    0.07
#  Residual             654.94   25.59
# Number of obs: 180, groups:  Subject, 18

Three numbers carry the bulk of the interpretation. The fixed-effect intercept (251 ms) is the population-average reaction time on day 0. The random-intercept SD (24.7 ms) says a typical subject’s day-0 reaction time deviates from the population by about 25 ms. The random-slope SD (5.9 ms/day) says subjects differ in how quickly they slow down under sleep deprivation, with a SD of about 6 ms per additional day. The correlation between intercept and slope is small (0.07), so a subject’s baseline tells you very little about their slope.

The grouping factor must be a factor. If Subject came in as a character vector it would still work in newer lme4 versions, but converting it explicitly is safer. See the factor reference for the conversion.

Reading the summary

summary() on a merMod object prints four blocks. Fixed effects give you the population coefficients. Random effects give you the variance components and, when the model has correlated random terms, the correlation matrix. The correlation of fixed effects block is the multicollinearity check for the fixed-effect design. The model fit block prints the AIC, BIC, log-likelihood, REML criterion, and the number of groups per random term.

The extractor generics are how you pull these numbers back out as data:

fixef(fm1)                # named numeric: population fixed effects
ranef(fm1)$Subject        # subject-level random intercepts and slopes
coef(fm1)$Subject         # fixed + per-subject random combined
VarCorr(fm1)              # variance components
confint(fm1, method = "Wald")
predict(fm1, newdata = new_df)

ranef() returns the conditional modes of the random effects, the best linear unbiased predictors in older terminology. coef() is the same plus the fixed effects, so the per-subject intercept in coef(fm1)$Subject[, "(Intercept)"] is the actual fitted day-0 reaction time for that subject.

Random slopes and why they matter

fm1 estimated both a random intercept and a random slope for Days. If you fit a random-intercept-only model you get a different answer for the Days coefficient because you forced the slope to be the same across subjects:

fm_int <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy)
anova(fm1, fm_int)

The likelihood-ratio test between these two models asks whether the additional random-slope variance is justified. With only 18 subjects the test is conservative; with hundreds, the difference can be dramatic.

A common refinement is to force the intercept-slope correlation to zero with ||:

fm_uncorr <- lmer(Reaction ~ Days + (Days || Subject), sleepstudy)

This estimates two independent variance components instead of the full 2x2 covariance matrix. Use it when the model is over-fit, the correlation is poorly estimated, or summary() reports a singular fit.

Nested and crossed random effects

The distinction between nested and crossed factors is a frequent source of confusion. If your class IDs are unique across schools (class_1, class_2 in school A and the same names in school B) you should still write (1 | school/class). lme4 will treat this as nested even if the inner IDs are not unique. If you instead write (1 | school) + (1 | class) you are telling lme4 that the levels of class are shared across schools, which inflates the number of groups and changes the variance decomposition.

Crossed random effects appear in designs where every subject sees every item, or every patient gets every drug:

# 30 subjects, 50 items, each subject sees about half the items
fm_crossed <- lmer(
  RT ~ condition + (1 | subject) + (1 | item),
  data  = my_df,
  REML  = FALSE
)

The subject and item terms are crossed. The grouping factor on the left of | is independent of the one on the right.

Generalized linear mixed models with glmer

For non-Gaussian responses, lme4::glmer extends lmer with a family argument that takes any family() object. The classic example is the cbpp dataset, a 56-row frame of bovine tuberculosis counts from 15 herds across 4 time periods:

data(cbpp, package = "lme4")
gm1 <- glmer(
  cbind(incidence, size - incidence) ~ period + (1 | herd),
  data   = cbpp,
  family = binomial(link = "logit")
)
summary(gm1)
# Fixed effects:
#             Estimate Std. Error z value Pr(>|z|)
# (Intercept)   -1.386      0.239   -5.80   7.2e-09 ***
# period2       -0.535      0.253   -2.12   0.0342 *
# period3       -0.778      0.266   -2.93   0.0034 **
# period4       -1.118      0.290   -3.86   0.0001 ***

The (Intercept) row is the log-odds of incidence in period 1. Each periodN row is the change in log-odds from period 1. The random-intercept SD across herds absorbs the herd-level variation in baseline incidence.

Group proportions are useful to keep in view while you read the model:

library(dplyr)
cbpp |> summarise(p = sum(incidence) / sum(size), .by = period)
#   period          p
# 1      1 0.07142857
# 2      2 0.04901961
# 3      3 0.03571429
# 4      4 0.02105263

The summarise() call here is the dplyr summarise pattern. For count data, pass family = poisson or, with overdispersion, family = MASS::negative.binomial(theta = ...). For proportions bounded on (0, 1) with no 0/1 values, family = betareg::Beta() inside a glmmTMB call is the usual route.

Structured residuals with nlme

lme4 is silent on the within-group error structure: residuals are assumed i.i.d. normal. When that assumption fails (variance grows with the fitted value, or adjacent residuals are correlated), nlme::lme lets you model the structure directly. The Orthodont dataset has growth measurements on 27 children at ages 8, 10, 12, and 14. Heteroscedasticity by sex is a common question in this kind of data:

library(nlme)
fm_lme <- lme(
  fixed   = distance ~ age + Sex,
  data    = Orthodont,
  random  = ~ age | Subject,
  weights = varPower()
)
summary(fm_lme)

varPower() says the residual variance is proportional to the absolute value of the fitted value raised to a power 2*delta. The estimated delta is reported in the summary. Other useful variance functions are varIdent() (different SDs per group), varExp(), and varConstPower(). For temporal autocorrelation, pass correlation = corAR1() to the lme() call.

The output of lme is an lme object, not a merMod. Most of the same extractor generics work (fixef, ranef, VarCorr, predict, residuals), but the anova() method follows nlme’s own conventions and the broom.mixed::tidy() interface may need different arguments.

REML vs ML and the model comparison trap

lmer defaults to REML = TRUE. Restricted maximum likelihood is unbiased for variance components, which is what you usually want for reporting a single model’s parameters. It is also the wrong likelihood to use when you compare models with different fixed effects, because the REML criterion depends on the fixed-effect design.

If you want to LRT a fixed effect, refit both candidates with REML = FALSE:

m1 <- lmer(Reaction ~ 1 + (Days | Subject), sleepstudy, REML = FALSE)
m2 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy, REML = FALSE)
anova(m1, m2)

For random-effects comparisons, REML is fine, and the test is more accurate. For fixed-effects comparisons, ML is required.

Diagnostics and common errors

The error message you will see most often is boundary (singular) fit: see help(isSingular). It means one of the variance components (or a correlation) was estimated at the boundary, usually 0 or ±1. The fix is to either drop the offending term, force the correlation to 0 with ||, or collect more data. The performance package has a check_model() helper that runs the standard residual, normality, and influential-observation plots in one call:

performance::check_model(fm1)

Other messages worth knowing:

Model failed to converge: try lmerControl(optimizer = "bobyqa") (the default in current lme4) or fall back to "Nelder_Mead". Rescale continuous predictors; simplify the random structure.
number of levels of grouping factor must be > 1: the grouping factor is constant. Check for typos or subsetting that wiped out the variation.
contrasts can be applied only to factors with 2 or more levels: a factor collapsed to one level after subsetting. Coerce with droplevels() before fitting.
fitted probabilities occur numerically 0 or 1 in glmer: complete separation. Try blme::bglmer for a weakly-informative prior on the fixed effects, or fit a Firth-corrected model.

A practical habit is to centre continuous predictors before fitting models with random slopes. Centring moves the random-intercept correlation off the boundary of the parameter space and makes the intercept interpretable as the response at the predictor’s mean.

Reporting results

For a paper, the minimum set of numbers is: the fixed-effect estimates with standard errors, the random-effect variance components, the number of groups and observations, and a fit statistic (AIC, BIC, or marginal R² from performance::r2_nakagawa(fm1)). The papaja package will format these inline in an R Markdown or Quarto document, including the random-effects table.

For a model with both fixed and random effects, broom.mixed::tidy(fm1, effects = "fixed") returns a tibble ready for gt or kable:

broom.mixed::tidy(fm1, effects = "fixed", conf.int = TRUE)
# # A tibble: 2 × 7
#   term        estimate std.error statistic   p.value conf.low conf.high
#   <chr>          <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
# 1 (Intercept)    251.       6.63     37.9  2.67e- 24    238.      265.
# 2 Days            10.5      1.50      6.97 5.05e-  9      7.51     13.4

For random effects, broom.mixed::tidy(fm1, effects = "ran_vals") or broom.mixed::tidy(fm1, effects = "ran_modes") give the variance components and the per-group modes respectively.

If you want the raw population variance as a sanity check, var is the base function to look at; the VarCorr() output for a random intercept can be cross-checked against var(tapply(response, group, mean)).