Tables, Random Effects and Plots in Multilevel Analysis

Workshop Guide

Authors

Kamarul Imran Musa

Afiq Amsyar

Published

September 10, 2025

Table, random effects and plots

Packages

A few useful packages to further improve the presentation of mixed model results.

Let’s load some packages

lme4 to analyze mixed models
tidyverse to wrangle data and make plots

Code

library(lme4)
library(tidyverse)

Dataset

We we use the built in sleepstudy dataset

Code

data("sleepstudy")

Random intercept model

We start of with random intercept model

Code

model <- lmer(Reaction ~ Days + (1 | Subject), 
              data = sleepstudy)

Random slope model

Next, we estimate random slope model

Code

model2 <- lmer(Reaction ~ Days + (Days | Subject), 
               data = sleepstudy)

Plots

Load additional packages to make plots for multilevel objects

Code

library(sjPlot)
library(lattice)

Plot random intercept model using sjPlot

Plot random effect

Code

plot_model(model, type = "re", 
           show.values = TRUE, 
           value.offset = .3,
           value.size = 3)

Plot fixed effect

Code

plot_model(model, type = "est", 
           show.values = TRUE, 
           value.size = 3)

Plot random slope model using sjPlot

Code

plot_model(model2, type = "re", 
           show.values = TRUE, value.offset = .4,
           value.size = 3)

Combining random effect and predicted value in a plot

Code

plot_model(
  model,
  type = "pred",
  terms = c("Days", "Subject"),
  pred.type = "re",
  ci.lvl = NA
)

Code

plot_model(
  model2,
  type = "pred",
  terms = c("Days", "Subject"),
  pred.type = "re",
  ci.lvl = NA,
  title = "My model 2",
  axis.title = "Reaction (in seconds)"
)

Caterpillar plots

random intercept (model)

Code

dotplot(ranef(model))

$Subject

random slope (model2)

Code

dotplot(ranef(model2))

$Subject

Tables

Useful packages include

Code

library(modelsummary)

using modelsummary to make results in a table

Code

modelsummary(model)

	(1)
(Intercept)	251.405
	(9.747)
Days	10.467
	(0.804)
SD (Intercept Subject)	37.124
SD (Observations)	30.991
Num.Obs.	180
R2 Marg.	0.280
R2 Cond.	0.704
AIC	1794.5
BIC	1807.2
ICC	0.6
RMSE	29.41

Code

modelsummary(model2)

	(1)
(Intercept)	251.405
	(6.825)
Days	10.467
	(1.546)
SD (Intercept Subject)	24.741
SD (Days Subject)	5.922
Cor (Intercept~Days Subject)	0.066
SD (Observations)	25.592
Num.Obs.	180
R2 Marg.	0.279
R2 Cond.	0.799
AIC	1755.6
BIC	1774.8
ICC	0.7
RMSE	23.44

we can combine the results into a table

Code

modelsummary(list(model, model2))

	(1)	(2)
(Intercept)	251.405	251.405
	(9.747)	(6.825)
Days	10.467	10.467
	(0.804)	(1.546)
SD (Intercept Subject)	37.124	24.741
SD (Days Subject)		5.922
Cor (Intercept~Days Subject)		0.066
SD (Observations)	30.991	25.592
Num.Obs.	180	180
R2 Marg.	0.280	0.279
R2 Cond.	0.704	0.799
AIC	1794.5	1755.6
BIC	1807.2	1774.8
ICC	0.6	0.7
RMSE	29.41	23.44

using tab_model

Code

tab_model(model)

	Reaction
Predictors	Estimates	CI	p
(Intercept)	251.41	232.17 – 270.64	<0.001
Days	10.47	8.88 – 12.05	<0.001
Random Effects
σ²	960.46
τ₀₀ _Subject	1378.18
ICC	0.59
N _Subject	18
Observations	180
Marginal R² / Conditional R²	0.280 / 0.704

Code

tab_model(model2)

	Reaction
Predictors	Estimates	CI	p
(Intercept)	251.41	237.94 – 264.87	<0.001
Days	10.47	7.42 – 13.52	<0.001
Random Effects
σ²	654.94
τ₀₀ _Subject	612.10
τ₁₁ _Subject.Days	35.07
ρ₀₁ _Subject	0.07
ICC	0.72
N _Subject	18
Observations	180
Marginal R² / Conditional R²	0.279 / 0.799

combining two tables

Code

tab_model(model, model2)

	Reaction			Reaction
Predictors	Estimates	CI	p	Estimates	CI	p
(Intercept)	251.41	232.17 – 270.64	<0.001	251.41	237.94 – 264.87	<0.001
Days	10.47	8.88 – 12.05	<0.001	10.47	7.42 – 13.52	<0.001
Random Effects
σ²	960.46			654.94
τ₀₀	1378.18 _Subject			612.10 _Subject
τ₁₁				35.07 _Subject.Days
ρ₀₁				0.07 _Subject
ICC	0.59			0.72
N	18 _Subject			18 _Subject
Observations	180			180
Marginal R² / Conditional R²	0.280 / 0.704			0.279 / 0.799

Hands-on

Create a new project and name the project
Install packages
- General package
- Data wrangling
- Read data

Code

library(tidyverse)
library(haven)
library(here)

For multilevel
For presentation

Code

library(lme4)
library(sjPlot)
library(modelsummary)

Read data

Code

immi <- read_dta('imm10_2.dta')

EDA and data wrangling

Code

glimpse(immi)

Rows: 260
Columns: 9
$ schid    <dbl> 7472, 7472, 7472, 7472, 7472, 7472, 7472, 7472, 7472, 7472, 7…
$ math     <dbl> 48, 48, 53, 42, 43, 57, 33, 64, 36, 56, 48, 48, 44, 35, 50, 3…
$ ses      <dbl> -0.13, -0.39, -0.80, -0.72, -0.74, -0.58, -0.83, -0.51, -0.56…
$ white    <dbl+lbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ parented <dbl+lbl> 2, 2, 2, 2, 2, 2, 2, 3, 2, 3, 2, 3, 2, 1, 2, 3, 3, 1, 3, …
$ public   <dbl+lbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
$ sctype   <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ sex      <dbl+lbl> 2, 1, 1, 1, 2, 2, 2, 1, 2, 2, 2, 1, 1, 2, 2, 2, 2, 2, 2, …
$ schnum   <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…

Code

immi <-
  immi |> 
  mutate(across(where(is.labelled), as_factor),
         schnum = factor(schnum))

linear mixed effect model
- random intercept
- and we assume math increase similarly based on ses for all schools

Code

ri_immi <- 
  lmer(math ~ ses +  white + parented + (1 | schnum), 
       data = immi, REML = FALSE)

random slope
and we assume math increase differently based on ses for all schools

Code

rs_immi <- 
  lmer(math ~ ses +  white + parented + 
         (1 + ses | schnum), data = immi, REML = FALSE)

Presentation
- Tables

Code

tab_model(list(ri_immi, rs_immi))

	Math score			Math score
Predictors	Estimates	CI	p	Estimates	CI	p
(Intercept)	43.93	38.82 – 49.03	<0.001	43.45	38.44 – 48.46	<0.001
Socioecnonomic Status	1.43	-1.16 – 4.02	0.276	1.34	-1.45 – 4.12	0.345
white: White	1.07	-2.05 – 4.19	0.500	0.14	-2.63 – 2.91	0.919
parented: Level 2	1.74	-2.02 – 5.51	0.363	2.11	-1.66 – 5.88	0.271
parented: Level 3	6.16	1.99 – 10.32	0.004	7.06	2.88 – 11.23	0.001
parented: Level 4	10.85	4.57 – 17.13	0.001	11.52	5.22 – 17.82	<0.001
parented: Level 5	9.66	2.54 – 16.77	0.008	9.54	2.38 – 16.71	0.009
parented: Level 6	8.21	0.28 – 16.14	0.042	7.31	-0.74 – 15.36	0.075
Random Effects
σ²	62.81			62.90
τ₀₀	8.21 _schnum			7.88 _schnum
τ₁₁				2.49 _schnum.ses
ρ₀₁				1.00 _schnum
ICC	0.12			0.13
N	10 _schnum			10 _schnum
Observations	260			260
Marginal R² / Conditional R²	0.266 / 0.351			0.241 / 0.341

Code

modelsummary(list(
  "random intercept" = ri_immi, 
  "random slope" = rs_immi),
             fmt = fmt_decimal(digits = 2, pdigits = 2))

	random intercept	random slope
(Intercept)	43.93	43.45
	(2.59)	(2.54)
ses	1.43	1.34
	(1.31)	(1.41)
whiteWhite	1.07	0.14
	(1.58)	(1.41)
parentedLevel 2	1.74	2.11
	(1.91)	(1.91)
parentedLevel 3	6.16	7.06
	(2.12)	(2.12)
parentedLevel 4	10.85	11.52
	(3.19)	(3.20)
parentedLevel 5	9.66	9.54
	(3.61)	(3.64)
parentedLevel 6	8.21	7.31
	(4.03)	(4.09)
SD (Intercept schnum)	2.86	2.81
SD (ses schnum)		1.58
Cor (Intercept~ses schnum)		1.00
SD (Observations)	7.93	7.93
Num.Obs.	260	260
R2 Marg.	0.266	0.268
R2 Cond.	0.351
AIC	1827.2	1830.1
BIC	1862.8	1872.8
ICC	0.1
RMSE	7.81	7.83

Code

modelsummary(list(
  "random intercept" = ri_immi, 
  "random slope" = rs_immi),
             fmt = fmt_decimal(digits = 2, pdigits = 2), 
             output = "table.docx")

Plots

Code

plot_model(ri_immi, 
           type = 'est')

Code

plot_model(ri_immi, 
           type = 'pred',
           terms = 'ses')

Code

plot_model(ri_immi, 
           type = 'pred',
           terms = c('ses' , 'white', 'parented'))

Code

plot_model(ri_immi, 
           type = 'pred',
           terms = c('ses' , 'white', 'schnum'),
           pred.type = 're', ci.lvl=NA)

Code

plot_model(ri_immi, 
           type = 're')

Code

plot_model(ri_immi, 
           type = 'diag')

[[1]]


[[2]]
[[2]]$schnum



[[3]]


[[4]]

Code

plot_model(rs_immi, 
           type = 'est')

Code

plot_model(rs_immi, 
           type = 're')

Code

plot_model(rs_immi, 
           type = 'pred',
           terms = 'ses')

Code

plot_model(rs_immi, 
           type = 'pred',
           terms = c('ses', 'white', 'parented'))

Code

plot_model(rs_immi, 
           type = 'pred',
           terms = c('ses' , 'schnum', 'white'),
           pred.type = 're', ci.lvl=NA)

Code

plot_model(rs_immi, 
           type = 'diag')

[[1]]


[[2]]
[[2]]$schnum



[[3]]


[[4]]