1. Introduction

Imagine you are an AFL analyst. You have the average number of disposals per game for 25 players from each of the 18 teams, along with the percentage of game time each player spends on the ground. We want to know how time on the ground is related to disposals.

There is a problem using ordinary linear regression in this situation. Players on the same team share a coach, a game plan, and a playing style. This way, players from the same team tend to be more similar to each other than to players from the other teams. Regression assumes every player is independent, so using it here would make our results look much better than they should be.

A linear mixed model fixes this. It has two parts:

  • Fixed effects are the average effects for everyone, like the coefficients in regression (e.g. the average change in disposals for each extra 1% of time on the ground).
  • Random effects allow each group (each team in this case) to differ from that average. A random intercept gives each team its own average number of disposals.

By the end of this guide, you will be able to fit a mixed model, work out how much of the variation is between teams, and see how a mixed model differs from ordinary regression.

Set up. Install the packages below (only once in the console). Then load them each time you start R.

install.packages("tidyverse")
install.packages("dplyr")
install.packages("sjPlot")
install.packages("performance")
install.packages("lme4")
install.packages("patchwork") # needed for check model
install.packages("see") # needed for check model
library(tidyverse)
library(dplyr)
library(sjPlot)
library(performance)
library(lme4)

2. The data

We simulate data for the 450 players (25 from each of the 18 teams), so no downloading or importing of data is needed. (With real data you will need to import your own file wth read.csv()).

  • team: a number identifying each team (1 to 18).
  • time_on_ground: the percentage of game time the player spends on the ground.
  • disposals: Outcome -> the average number of disposals (kicks and handballs) per game.

The code below builds the data. You don’t need to memorise it. Each comment explains what each part does, the actual analysis will start after this chunk.

set.seed(123) # gives us the same random numbers/data every time
n_teams <- 18 # number of teams
players_per_team <- 25 # number of players per team

# Each team gets its own boost (or drop) in disposals due to its game plan.

# This is the random intercept
team_effect <- rnorm(n_teams, mean = 0, sd = 2.5)

# Build the data: 25 rows (players) for each team
# Disposals = 2 + 0.2 x time on ground + the team's effect + random variation.
# rnorm(n(), 0, 2.5) adds random variation, as no two players are identical.

# rep() repeats each team number 25 times, so every player gets a team
players <- data.frame(team = rep(1:n_teams, each = players_per_team)) |> 
  mutate(time_on_ground = round(runif(n(), 55, 95)),
         disposals = round(2 + 0.2 * time_on_ground + team_effect[team] + rnorm(n(), 0, 2.5))) |> 
  as_tibble()

Inspect the data with str:

str(players)
## tibble [450 × 3] (S3: tbl_df/tbl/data.frame)
##  $ team          : int [1:450] 1 1 1 1 1 1 1 1 1 1 ...
##  $ time_on_ground: num [1:450] 85 64 68 64 61 72 72 70 61 61 ...
##  $ disposals     : num [1:450] 15 13 19 13 9 13 16 14 11 12 ...

We can see that team is stored as a number, but we know it’s a label for a group. For mixed modelling it needs to be a factor, so we convert it with mutate().

players <- players |> 
  mutate(team = as.factor(team))

# View it again to check the change to factors
str(players)
## tibble [450 × 3] (S3: tbl_df/tbl/data.frame)
##  $ team          : Factor w/ 18 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ time_on_ground: num [1:450] 85 64 68 64 61 72 72 70 61 61 ...
##  $ disposals     : num [1:450] 15 13 19 13 9 13 16 14 11 12 ...

3. Exploring the data

Before building any model, we explore the data. First, a histogram of the outcome is made to see how disposals are spread:

ggplot(players, aes(x = disposals)) +
  geom_histogram()

Next, a scatter plot of disposals against time on the ground, for all players together, as well as with a regression line:

ggplot(players, aes(x = time_on_ground, y = disposals)) + 
  geom_point() +
  geom_smooth(method = "lm")

Finally, the same plot with facet_wrap(), which gives each team its own panel:

ggplot(players, aes(x = time_on_ground, y = disposals)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap("team")

From these visualisation plots, we can highlight that players who spend more time on the ground have more disposals. But if we look at the team panels, some teams sit higher on the plot than others. Each team also has its own average level of disposals, which is exactly what a random intercept captures.

Factors such as team game style and coaching can vary across teams and can therefore be a reason to why each team can be different to one another.

4. Ordinary regression (ignoring teams)

First, we fit an ordinary linear regression, treating all 450 players as independent. This is what we learnt in the previous regression guide. (Check out the ‘Learning to Performing Linear regression in R’ first if you haven’t!).

m1 <-  lm(disposals ~ time_on_ground, data = players)

tab_model(m1)
  disposals
Predictors Estimates CI p
(Intercept) 1.57 -0.58 – 3.71 0.152
time on ground 0.21 0.18 – 0.24 <0.001
Observations 450
R2 / R2 adjusted 0.322 / 0.320

tab_model() shows the estimate for time_on_ground, its confidence interval (CI), and R-squared. The intercept is the predicted disposals at 0% time on the ground, which isn’t a realistic player, so we focus on the slope.

Results: Each extra 1% of time on ground is linked to about 0.21 more disposals per game. Time on the ground alone explains about 32% of the differences between players.

5. The unconditional means model or the ‘NULL’ model

The simplest mixed model has no predictors (NULL). It only has a random intercept, written (1 | team), which gives each team its own average. This model tells us how much of the variation in disposals is due to the differences between teams. We fit mixed models with lmer().

m2 <- lmer(disposals ~ 1 + (1 | team), data = players)

tab_model(m2)
  disposals
Predictors Estimates CI p
(Intercept) 17.28 16.06 – 18.49 <0.001
Random Effects
σ2 11.29
τ00 team 6.47
ICC 0.36
N team 18
Observations 450
Marginal R2 / Conditional R2 0.000 / 0.364

tab_model() shows the average disposals for all players (the intercept) and the variance split between teams and within teams. It also gives the ICC (intraclass correlation coefficient), which is the amount of total variation that is due to differences between teams. A higher ICC implies players from the same team are more similar to each other.

Results: The average number of disposals across players is about 17.28 per game. The ICC is 0.36, which means about 36% of the differences in disposals are between teams. This is not small, so players from the same team are related and a mixed model is needed.

6. Adding time on the ground

Now we add time_on_ground to the mixed model. Each team still gets its own average (the random intercept), but all teams share the same fixed slope for the time on the ground.

m3 <- lmer(disposals ~ time_on_ground + (1 | team), data = players)

tab_model(m1, m3)
  disposals disposals
Predictors Estimates CI p Estimates CI p
(Intercept) 1.57 -0.58 – 3.71 0.152 2.32 0.39 – 4.26 0.019
time on ground 0.21 0.18 – 0.24 <0.001 0.20 0.18 – 0.22 <0.001
Random Effects
σ2   6.18
τ00   5.98 team
ICC   0.49
N   18 team
Observations 450 450
R2 / R2 adjusted 0.322 / 0.320 0.294 / 0.641

tab_model(m1, m3) puts the ordinary regression and the mixed model side by side. Now we are able to compare the time_on_ground estimate and its CI in the two models.

Results: In the mixed model, each extra 1% of time on the ground is linked to about 0.20 more disposals per game (95% CI: 0.18 to 0.22), compared with 0.21 in the ordinary regression. The estimate is similar, but the CI is narrower in the mixed model as it accounts for the differences between teams. The team variance is 5.98, which shows teams still differ in their average disposals after time on the ground is accounted for.

plot_model() with type = "re" plots each team’s random intercept. Each dot shows how far that teams’ average is above or below the overall average, so teams to the right have more disposals than average.

plot_model(m3, type = "re")

As with regression, we then check the model’s assumptions. check_model() draws the diagnostics plots, and each panel explains what a good plot looks like.

check_model(m3)

7. Comparing models

Which mixed model is better, the one with no predictors (m2) or the one with time on the ground (m3)? We compare them with the AICc, where lower is better.

tab_model(m2, m3,
          show.aicc = TRUE)
  disposals disposals
Predictors Estimates CI p Estimates CI p
(Intercept) 17.28 16.06 – 18.49 <0.001 2.32 0.39 – 4.26 0.019
time on ground 0.20 0.18 – 0.22 <0.001
Random Effects
σ2 11.29 6.18
τ00 6.47 team 5.98 team
ICC 0.36 0.49
N 18 team 18 team
Observations 450 450
Marginal R2 / Conditional R2 0.000 / 0.364 0.294 / 0.641
AICc 2421.261 2167.624

anova() also compares the models and tests whether the second fits better than the first. A small p-value means the extra term was worth adding.

anova(m2, m3)
## refitting model(s) with ML (instead of REML)
## Data: players
## Models:
## m2: disposals ~ 1 + (1 | team)
## m3: disposals ~ time_on_ground + (1 | team)
##    npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)    
## m2    3 2422.1 2434.4 -1208.0    2416.1                         
## m3    4 2161.0 2177.5 -1076.5    2153.0 263.04  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Results: The model with time on the ground (m3) has a lower AICc than the model without it (2167.624 compared to 2421.261), and the p-value from anova() is small, so adding time on the ground made the model better.

8. Key points

  • Players on the same team are related, so ordinary regression isn’t appropriate.
  • Fixed effects are the average effects for everyone. Random effects let each group (here, each team) vary from the average.
  • A random intercept (1 | team) gives each team its own average.
  • The ICC tells us how much of the variation is between groups.
  • Compare models using AICc, and anova().