Multilevel models

lme4 for MLM regression; lavaan for two-level SEM

Author
Affiliation

George Michaelides

Norwich Business School, University of East Anglia

Published

July 7, 2026

Learning objectives

By the end of this tutorial you will be able to:

  • describe what a multilevel model is and when it is appropriate
  • fit a random-intercept and random-slope model with lme4::lmer()
  • interpret VarCorr() output and compute the intraclass correlation coefficient (ICC)
  • describe how lavaan extends to two-level SEM and recognise the syntax
  • decide when to use multilevel regression versus multilevel SEM

Setup

Setup code 
library(lme4)
library(lavaan)
library(ggplot2)

What is a multilevel model?

Multilevel models (also called hierarchical linear models or mixed-effects models) are for data with a nested structure — observations clustered within higher-level units. Common examples in organisational psychology:

  • Employees nested within teams
  • Teams nested within departments
  • Repeated measures nested within individuals (longitudinal data — see Tutorial 15)

The defining feature is that observations within the same cluster are not independent — two employees in the same team tend to be more similar than two employees from different teams. Standard regression assumes independence, so applying it to nested data gives biased standard errors and incorrect p-values.

Multilevel models handle this by modelling the clustering explicitly. They estimate:

  • Fixed effects — the average relationship across all clusters (analogous to regression coefficients)
  • Random effects — how individual clusters deviate from the average (variance components)

When to use a multilevel model

Use a multilevel model when:

  1. Your data has a nested structure (e.g. employees within teams, repeated measures within individuals)
  2. There is meaningful variation between clusters (test this with the ICC — see below)
  3. You have enough clusters (a rough rule of thumb: 20+ clusters with 5+ observations each)

Do not use a multilevel model if:

  • Your data is single-level (no clustering) — use lm() or lavaan instead
  • You have very few clusters (fewer than 10) — the random effects estimates will be unstable
  • The clustering is incidental (e.g. data collected on two different days, with no substantive meaning to the day) — you can often ignore it

The teaching dataset

The PS dataset we have used throughout this series is single-level — one row per participant, no clustering. To demonstrate multilevel modelling, we use a simulated dataset of employees nested within teams, generated by data-raw/simulate-multilevel.R.

ml <- read.csv("data/multilevel_sim.csv")
head(ml)
##   employee_id team_id     team_x      ind_x        y
## 1       E0001       1 -0.1102855 -0.5817275 1.697048
## 2       E0002       1 -0.1102855  0.7000802 2.313822
## 3       E0003       1 -0.1102855  1.4921766 3.106420
## 4       E0004       1 -0.1102855  0.5265534 1.570512
## 5       E0005       1 -0.1102855  1.0377721 1.971969
## 6       E0006       1 -0.1102855 -1.8607164 2.839831
nrow(ml)
## [1] 348
length(unique(ml$team_id))
## [1] 30

The dataset has 348 employees across 30 teams, with an average team size of about 12. Each row has:

  • employee_id — a unique employee identifier
  • team_id — the team the employee belongs to (1–30)
  • team_x — a team-level predictor (a characteristic of the team, shared by all members)
  • ind_x — an individual-level predictor (a characteristic of the employee)
  • y — the outcome variable

Multilevel regression with lme4

The lme4 package fits multilevel models in R. The main function is lmer() — “linear mixed-effects regression”.

Random intercept model

The simplest multilevel model adds a random intercept for each cluster. This allows each team to have its own baseline level of y, while keeping the slopes (the effects of ind_x and team_x) the same across teams.

The syntax extends the lm() formula with a (1 | group) term:

fit_ri <- lmer(y ~ ind_x + team_x + (1 | team_id), data = ml)
summary(fit_ri)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ ind_x + team_x + (1 | team_id)
##    Data: ml
## 
## REML criterion at convergence: 842.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.88841 -0.65618 -0.03109  0.60856  2.68643 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  team_id  (Intercept) 0.371    0.6091  
##  Residual             0.537    0.7328  
## Number of obs: 348, groups:  team_id, 30
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  2.51366    0.13162  19.098
## ind_x        0.44701    0.04169  10.722
## team_x       0.37428    0.12879   2.906
## 
## Correlation of Fixed Effects:
##        (Intr) ind_x 
## ind_x  -0.013       
## team_x  0.440 -0.029

Reading the output:

  • Fixed effects — the Estimate column gives the coefficients for ind_x and team_x, interpreted exactly as in lm(). The intercept is the predicted y when all predictors are zero, for a team with an average random effect.
  • Random effects — the (Intercept) variance under team_id:(Intercept) tells you how much team means vary around the grand mean. The Residual variance is the within-team variance.

The intraclass correlation coefficient (ICC)

The ICC tells you what proportion of the total variance in y is between teams rather than within teams:

\[ \text{ICC} = \frac{\sigma^2_{\text{between}}}{\sigma^2_{\text{between}} + \sigma^2_{\text{within}}} \tag{1}\]

Extract the variance components:

vc <- VarCorr(fit_ri)
vc
##  Groups   Name        Std.Dev.
##  team_id  (Intercept) 0.60906 
##  Residual             0.73279

Compute the ICC manually:

var_between <- attr(vc$team_id, "stddev")^2
var_within  <- attr(vc, "sc")^2
icc <- var_between / (var_between + var_within)
icc
## (Intercept) 
##   0.4085704

An ICC above .05 suggests meaningful clustering — a multilevel model is warranted. An ICC near 0 means the clustering explains very little variance and a single-level model may suffice.

NoteInterpreting the ICC

The ICC is the proportion of variance in the outcome that is due to differences between clusters rather than differences within clusters. An ICC of .20 means 20% of the variance in y is between teams, and 80% is within teams. Higher ICC means stronger clustering — and a stronger need for a multilevel model.

Random slope model

To allow the effect of ind_x to vary across teams — some teams may show a stronger relationship between ind_x and y than others — add a random slope:

fit_rs <- lmer(y ~ ind_x + team_x + (1 + ind_x | team_id), data = ml)
summary(fit_rs)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ ind_x + team_x + (1 + ind_x | team_id)
##    Data: ml
## 
## REML criterion at convergence: 834.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.68860 -0.67861  0.00018  0.64023  2.43519 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  team_id  (Intercept) 0.36606  0.6050        
##           ind_x       0.03118  0.1766   0.71 
##  Residual             0.50904  0.7135        
## Number of obs: 348, groups:  team_id, 30
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  2.50883    0.12962  19.355
## ind_x        0.43392    0.05245   8.273
## team_x       0.34687    0.11956   2.901
## 
## Correlation of Fixed Effects:
##        (Intr) ind_x 
## ind_x   0.353       
## team_x  0.425 -0.042

The (1 + ind_x | team_id) term says “allow each team to have its own intercept AND its own slope for ind_x”. The random effects section now shows:

  • The variance of the team intercepts ((Intercept))
  • The variance of the team slopes (ind_x)
  • The correlation between intercepts and slopes (Corr)
TipTry it

Fit a random-intercept model with ind_x as the only predictor (no team_x). Compare the ICC to the model with team_x included. Does adding team_x change the ICC? Why?

fit_ri2 <- lmer(y ~ ind_x + (1 | team_id), data = ml)
vc2 <- VarCorr(fit_ri2)
var_b2 <- attr(vc2$team_id, "stddev")^2
var_w2 <- attr(vc2, "sc")^2
var_b2 / (var_b2 + var_w2)
## (Intercept) 
##   0.4718956

Adding team_x (a team-level predictor) should reduce the between-team variance, because some of what looked like “team differences” is actually explained by team_x. The ICC should be smaller in the model with team_x included.

Comparing random-intercept and random-slope models

Compare the two models with a likelihood-ratio test:

anova(fit_ri, fit_rs)
## Data: ml
## Models:
## fit_ri: y ~ ind_x + team_x + (1 | team_id)
## fit_rs: y ~ ind_x + team_x + (1 + ind_x | team_id)
##        npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)  
## fit_ri    5 843.61 862.87 -416.80    833.61                       
## fit_rs    7 838.97 865.93 -412.48    824.97 8.6428  2    0.01328 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A significant result means the random slope improves the model — the effect of ind_x does vary across teams. A non-significant result means the simpler random-intercept model is sufficient.

lavaan for two-level SEM

lavaan can also fit multilevel models, but it goes beyond what lme4 does: it allows latent variables at each level. In a two-level SEM, you can model:

  • A measurement model at Level 1 (within individuals)
  • A structural model at Level 2 (between teams)
  • Different factor structures at each level

This is powerful but complex. We provide a conceptual overview and a small worked example.

The level: syntax

lavaan’s two-level syntax uses the level: keyword to separate within-cluster and between-cluster model specifications:

model_2level <- '
  level: 1
    fw =~ y1 + y2 + y3

  level: 2
    fb =~ y1 + y2 + y3
    fb ~ team_x
'

The level: 1 block specifies the within-cluster model (how individual-level indicators co-vary within teams). The level: 2 block specifies the between-cluster model (how team-level aggregates co-vary across teams). The same observed variables can have different factor structures at each level.

Fitting a two-level model

To fit a two-level model, pass cluster = "team_id" to sem():

fit_2level <- sem(model_2level, data = ml, cluster = "team_id",
                  estimator = "MLR")
WarningTwo-level SEM requires a different dataset

The simulated multilevel_sim.csv has only one outcome (y) per individual, which is not enough to demonstrate a two-level SEM with latent variables at each level. A proper two-level SEM needs at least three indicators per individual to define a within-level factor.

For a real two-level SEM in your dissertation, you would need multivariate data — multiple items per employee, nested within teams. The syntax above shows how such a model would be specified; fitting it requires a dataset with the right structure.

When to use lme4 versus lavaan for multilevel modelling

Use lme4 (lmer()) when… Use lavaan (sem() with cluster =) when…
You have a single outcome variable You have multiple indicators per individual
You are interested in regression coefficients You want to model latent variables at either level
You do not need a measurement model You want to separate within-cluster and between-cluster factor structures
Your data is continuous and roughly normal You need fit indices to compare models

For most MSc dissertation applications, lme4 is sufficient. lavaan’s two-level SEM is a PhD-level technique, but it is worth knowing that the capability exists.

Common errors

Warningboundary (singular) fit: see ?isSingular

The random effects variance is estimated as zero or near-zero — the clusters do not differ enough to justify a random effect. Try removing the random effect and using lm() instead, or simplify the random-effects structure (e.g. remove a random slope if you have one).

number of levels of each grouping factor must be < number of observations

You have more clusters than observations — usually means your grouping variable is unique per row (e.g. using employee_id as the cluster when each employee has only one observation). Check that your grouping variable actually identifies clusters with multiple members.

Convergence warnings

The model is too complex for the data. Reduce the number of random effects, collect more data, or use a simpler model. Convergence problems in lmer() are common with small numbers of clusters.

Exercises

  1. Fit a random-intercept model predicting y from ind_x with team_id as the cluster. Compute the ICC. Is a multilevel model justified?

  2. Add team_x to the model from Exercise 1. Does the between-team variance decrease? Why? Compare the two models with anova().

  3. Fit a random-slope model allowing the effect of ind_x to vary across teams. Compare it to the random-intercept model from Exercise 2 with anova(). Is the random slope justified?

  4. Write a one-paragraph APA-style results summary for the model from Exercise 2, reporting the fixed effects, the variance components, and the ICC.

Exercise 1

fit1 <- lmer(y ~ ind_x + (1 | team_id), data = ml)
summary(fit1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ ind_x + (1 | team_id)
##    Data: ml
## 
## REML criterion at convergence: 848.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9124 -0.6449 -0.0343  0.6069  2.6983 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  team_id  (Intercept) 0.4797   0.6926  
##  Residual             0.5368   0.7327  
## Number of obs: 348, groups:  team_id, 30
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  2.34552    0.13263   17.68
## ind_x        0.45169    0.04173   10.82
## 
## Correlation of Fixed Effects:
##       (Intr)
## ind_x -0.001
vc1 <- VarCorr(fit1)
var_b <- attr(vc1$team_id, "stddev")^2
var_w <- attr(vc1, "sc")^2
icc1 <- var_b / (var_b + var_w)
icc1
## (Intercept) 
##   0.4718956

If the ICC is above .05, a multilevel model is justified — there is meaningful variation between teams.

Exercise 2

fit2 <- lmer(y ~ ind_x + team_x + (1 | team_id), data = ml)
summary(fit2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ ind_x + team_x + (1 | team_id)
##    Data: ml
## 
## REML criterion at convergence: 842.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.88841 -0.65618 -0.03109  0.60856  2.68643 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  team_id  (Intercept) 0.371    0.6091  
##  Residual             0.537    0.7328  
## Number of obs: 348, groups:  team_id, 30
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  2.51366    0.13162  19.098
## ind_x        0.44701    0.04169  10.722
## team_x       0.37428    0.12879   2.906
## 
## Correlation of Fixed Effects:
##        (Intr) ind_x 
## ind_x  -0.013       
## team_x  0.440 -0.029
vc2 <- VarCorr(fit2)
var_b2 <- attr(vc2$team_id, "stddev")^2
var_w2 <- attr(vc2, "sc")^2
var_b2 / (var_b2 + var_w2)
## (Intercept) 
##   0.4085704
anova(fit1, fit2)
## Data: ml
## Models:
## fit1: y ~ ind_x + (1 | team_id)
## fit2: y ~ ind_x + team_x + (1 | team_id)
##      npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)   
## fit1    4 849.51 864.91 -420.75    841.51                        
## fit2    5 843.61 862.87 -416.80    833.61 7.8965  1   0.004953 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Adding team_x should reduce the between-team variance because some of the team-level differences are explained by team_x. The anova() test tells you whether team_x significantly improves the model.

Exercise 3

fit3 <- lmer(y ~ ind_x + team_x + (1 + ind_x | team_id), data = ml)
anova(fit2, fit3)
## Data: ml
## Models:
## fit2: y ~ ind_x + team_x + (1 | team_id)
## fit3: y ~ ind_x + team_x + (1 + ind_x | team_id)
##      npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)  
## fit2    5 843.61 862.87 -416.80    833.61                       
## fit3    7 838.97 865.93 -412.48    824.97 8.6428  2    0.01328 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A significant result means the random slope improves the model — the effect of ind_x varies across teams. A non-significant result means the simpler random-intercept model is sufficient.

Exercise 4

Using the output from Exercise 2, write: “A multilevel model was fitted with employees nested within teams (N = 348, 30 teams). The intraclass correlation coefficient was ICC = …, indicating that …% of the variance in y was between teams. The effect of ind_x was \(b = ...\), \(p < .001\). The effect of team_x was \(b = ...\), \(p = ...\). The between-team variance was … and the within-team variance was ….”

Summary & next tutorial

You can now fit multilevel regression models with lme4::lmer(), interpret the variance components and ICC, and recognise the syntax for two-level SEM in lavaan. For most MSc applications, lme4 is the right tool; lavaan’s two-level SEM extends this to latent variables at each level. In the next tutorial we cover growth curve models — a specific type of multilevel model for longitudinal data, where repeated measures are nested within individuals.

Further reading

  • R Core Team (2024) — see the documentation for nlme and lme4 packages.
  • Wickham (2019) — for the deeper R programming underlying mixed-effects models.
  • Hoyle (2012) — the Handbook of Structural Equation Modeling includes chapters on multilevel SEM.

References

Hoyle, R. H. (2012). Handbook of structural equation modeling.
R Core Team. (2024). R: A language and environment for statistical computing. https://www.R-project.org/
Wickham, H. (2019). Advanced R (2nd ed.). CRC Press.