Modelling Player Performance with Linear Mixed Models

This tutorial examines player performance with generated demo data. It progresses through six ordinary linear models, six mixed models with a team random intercept, compares alternative random-effects structures, and checks diagnostics and team-specific fitness slopes. The simulated records make the analysis reproducible without sharing real player data.

Question: How do previous-season performance, fitness, goals, team facilities, and games played relate to current performance, after accounting for players grouped within teams?

The generated values are for learning only. They are not evidence about real athletes, and associations do not prove cause and effect.

1. Get ready

Use RStudio. Create a new R Markdown file, replace its contents with this tutorial, and save it as sport-performance-mixed-models.Rmd. Run chunks in order using each green triangle. Install these packages once in the Console if prompted:

install.packages(c("tidyverse", "lme4", "sjPlot", "performance", "see", "corrplot"))

When all chunks run, click Knit to make an HTML page. lme4 fits mixed models; sjPlot presents model tables; performance and see provide diagnostic plots.

library(tidyverse)
library(lme4)
library(sjPlot)
library(performance)
library(see)

2. Generate a demo player dataset

Generate example players with a facility category and a shared team effect. Players have previous-season scores, fitness, goals, games played, and other descriptors. This gives us data to demonstrate the analysis steps.

set.seed(2026)
n_teams <- 18
players_per_team <- 20
n_players <- n_teams * players_per_team

team_facility <- rep(c("Low", "Medium", "High"), length.out = n_teams)
team_intercept <- rnorm(n_teams, mean = 0, sd = 3)
team_fitness_slope <- rnorm(n_teams, mean = 0, sd = 0.45)
team_number <- rep(seq_len(n_teams), each = players_per_team)

players <- tibble(
  team_id = team_number,
  player_id = seq_len(n_players),
  age = sample(15:25, n_players, replace = TRUE),
  experience = pmin(age - 14, sample(0:8, n_players, replace = TRUE)),
  games_played = pmin(22, rpois(n_players, lambda = 16)),
  goals_scored = rpois(n_players, lambda = 8),
  injury_history = rbinom(n_players, size = 1, prob = 0.22),
  training_hours_weekly = pmax(0, rnorm(n_players, mean = 8, sd = 2)),
  previous_season_score = rnorm(n_players, mean = 55, sd = 8),
  fitness_level = rnorm(n_players, mean = 0, sd = 1),
  team_training_facilities = team_facility[team_number]
) |>
  mutate(
    performance_score = 12 +
      0.62 * previous_season_score +
      1.7 * fitness_level +
      0.20 * goals_scored +
      0.12 * games_played +
      case_when(
        team_training_facilities == "Medium" ~ 1.0,
        team_training_facilities == "High" ~ 2.0,
        TRUE ~ 0.0
      ) +
      team_intercept[team_number] +
      team_fitness_slope[team_number] * fitness_level +
      rnorm(n(), mean = 0, sd = 4)
  )

head(players)
## # A tibble: 6 × 12
##   team_id player_id   age experience games_played goals_scored injury_history
##     <int>     <int> <int>      <dbl>        <dbl>        <int>          <int>
## 1       1         1    23          0           14            9              0
## 2       1         2    21          4           15            6              1
## 3       1         3    18          1           21            8              0
## 4       1         4    23          3           12            8              0
## 5       1         5    20          1           14            1              1
## 6       1         6    18          4            9           10              0
## # ℹ 5 more variables: training_hours_weekly <dbl>, previous_season_score <dbl>,
## #   fitness_level <dbl>, team_training_facilities <chr>,
## #   performance_score <dbl>

set.seed() makes the generated data repeatable. team_id and player_id are identifiers, so convert them to factors before modelling. Set Low as the reference category for facilities so other categories are compared with it.

player_data <- players |>
  mutate(
    team_training_facilities = factor(
      team_training_facilities,
      levels = c("Low", "Medium", "High")
    ),
    across(c(team_id, player_id), factor)
  )

stopifnot(levels(player_data$team_training_facilities)[1] == "Low")
glimpse(player_data)
## Rows: 360
## Columns: 12
## $ team_id                  <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ player_id                <fct> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14…
## $ age                      <int> 23, 21, 18, 23, 20, 18, 23, 15, 22, 20, 18, 2…
## $ experience               <dbl> 0, 4, 1, 3, 1, 4, 6, 1, 1, 6, 4, 1, 1, 2, 8, …
## $ games_played             <dbl> 14, 15, 21, 12, 14, 9, 9, 16, 10, 18, 15, 18,…
## $ goals_scored             <int> 9, 6, 8, 8, 1, 10, 8, 11, 5, 6, 8, 12, 4, 14,…
## $ injury_history           <int> 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, …
## $ training_hours_weekly    <dbl> 9.180716, 8.197151, 7.665934, 9.631859, 5.427…
## $ previous_season_score    <dbl> 45.13630, 58.88130, 60.71854, 63.14258, 59.08…
## $ fitness_level            <dbl> -2.3534674, -0.2702758, -0.6892603, -0.211318…
## $ team_training_facilities <fct> Low, Low, Low, Low, Low, Low, Low, Low, Low, …
## $ performance_score        <dbl> 37.28223, 50.76638, 49.62395, 51.54077, 56.36…
colSums(is.na(player_data))
##                  team_id                player_id                      age 
##                        0                        0                        0 
##               experience             games_played             goals_scored 
##                        0                        0                        0 
##           injury_history    training_hours_weekly    previous_season_score 
##                        0                        0                        0 
##            fitness_level team_training_facilities        performance_score 
##                        0                        0                        0
n_distinct(player_data$team_id)
## [1] 18

3. Explore the data

Check younger players, experience, correlations, and weekly training. These quick checks help understand the variables before fitting models.

sum(player_data$age < 16)
## [1] 33
player_data |>
  filter(age < 16) |>
  select(team_id, player_id, age, experience, games_played,
         goals_scored, injury_history, performance_score) |>
  arrange(age) |>
  slice_head(n = 10)
## # A tibble: 10 × 8
##    team_id player_id   age experience games_played goals_scored injury_history
##    <fct>   <fct>     <int>      <dbl>        <dbl>        <int>          <int>
##  1 1       8            15          1           16           11              0
##  2 1       17           15          1           16            4              1
##  3 1       20           15          0           15            6              0
##  4 2       21           15          0           22           13              0
##  5 2       25           15          1           16            6              0
##  6 2       28           15          1           16            8              0
##  7 2       34           15          1           18           11              1
##  8 2       39           15          0           12           12              1
##  9 3       45           15          1           17            5              0
## 10 3       54           15          1           12            6              0
## # ℹ 1 more variable: performance_score <dbl>

Experience can be interpreted alongside the age when a player began. Correlations describe pairwise linear relationships; they do not control for other variables.

player_data |>
  mutate(debut_age = age - experience) |>
  ggplot(aes(debut_age)) +
  geom_histogram(bins = 30) +
  labs(x = "Approximate debut age", y = "Players") +
  theme_minimal()

cor(player_data$age, player_data$experience)
## [1] 0.4294302
cor_data <- player_data |> select(where(is.numeric))
sort(cor(cor_data)[, "performance_score"])
##        injury_history                   age training_hours_weekly 
##         -0.0264917656         -0.0003949346          0.0117745529 
##            experience          goals_scored          games_played 
##          0.0344136504          0.0615678140          0.1083684696 
##         fitness_level previous_season_score     performance_score 
##          0.2731264830          0.6856019986          1.0000000000
corrplot::corrplot(cor(cor_data), method = "color", type = "upper",
                   addCoef.col = "black", number.cex = 0.6, tl.cex = 0.7)

ggplot(player_data, aes(training_hours_weekly, performance_score)) +
  geom_point(alpha = 0.5) +
  geom_smooth() +
  labs(x = "Weekly training hours", y = "Current performance score") +
  theme_minimal()
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

4. Build six ordinary linear models

Add predictors in stages. Each model uses the same outcome and rows, so information criteria can be compared. AIC and AICc balance model fit and complexity; smaller values are preferred within this candidate set.

m0 <- lm(performance_score ~ 1, data = player_data)
m1 <- lm(performance_score ~ previous_season_score, data = player_data)
m2 <- lm(performance_score ~ previous_season_score + fitness_level,
         data = player_data)
m3 <- lm(performance_score ~ previous_season_score + fitness_level + goals_scored,
         data = player_data)
m4 <- lm(performance_score ~ previous_season_score + fitness_level + goals_scored +
           team_training_facilities, data = player_data)
m5 <- lm(performance_score ~ previous_season_score + fitness_level + goals_scored +
           team_training_facilities + games_played, data = player_data)

tab_model(m0, m1, m2, m3, m4, m5, show.aic = TRUE, show.aicc = TRUE)
  performance score performance score performance score performance score performance score performance score
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) 49.27 48.53 – 50.01 <0.001 15.33 11.54 – 19.11 <0.001 14.36 10.94 – 17.77 <0.001 11.73 8.00 – 15.46 <0.001 13.00 9.19 – 16.82 <0.001 9.69 5.47 – 13.91 <0.001
previous season score 0.62 0.55 – 0.69 <0.001 0.64 0.58 – 0.70 <0.001 0.65 0.59 – 0.71 <0.001 0.64 0.58 – 0.70 <0.001 0.64 0.58 – 0.70 <0.001
fitness level 2.23 1.75 – 2.70 <0.001 2.22 1.75 – 2.69 <0.001 2.14 1.67 – 2.60 <0.001 2.20 1.73 – 2.66 <0.001
goals scored 0.27 0.11 – 0.44 0.001 0.28 0.11 – 0.45 0.001 0.28 0.11 – 0.44 0.001
team training facilities
[Medium] 		-0.87 -2.04 – 0.29 0.142 -0.80 -1.95 – 0.35 0.174
team training facilities
[High] 		-1.95 -3.13 – -0.78 0.001 -1.91 -3.07 – -0.76 0.001
games played 0.22 0.09 – 0.34 0.001
Observations 360 360 360 360 360 360
R2 / R2 adjusted 0.000 / 0.000 0.470 / 0.469 0.572 / 0.570 0.584 / 0.581 0.597 / 0.591 0.609 / 0.603
AIC 2441.198 2214.608 2139.599 2131.263 2124.439 2114.940
AICc 2441.232 2214.675 2139.711 2131.432 2124.757 2115.350

Compare fitted and observed values for the m4 model.

ggplot(m4$model, aes(x = fitted(m4), y = performance_score)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "grey40") +
  geom_point(alpha = 0.4, colour = "steelblue") +
  coord_fixed() +
  labs(x = "Predicted performance score", y = "Observed performance score") +
  theme_minimal()

5. Add team random intercepts

Players within the same team may share context. Fit six matching mixed models with a team-specific random intercept. (1 | team_id) allows each team its own baseline while estimating average predictor associations across teams. Each model uses maximum likelihood (REML = FALSE) to support model comparisons.

mm0 <- lmer(performance_score ~ 1 + (1 | team_id), data = player_data, REML = FALSE)
mm1 <- lmer(performance_score ~ previous_season_score + (1 | team_id),
            data = player_data, REML = FALSE)
mm2 <- lmer(performance_score ~ previous_season_score + fitness_level + (1 | team_id),
            data = player_data, REML = FALSE)
mm3 <- lmer(performance_score ~ previous_season_score + fitness_level + goals_scored +
              (1 | team_id), data = player_data, REML = FALSE)
mm4 <- lmer(performance_score ~ previous_season_score + fitness_level + goals_scored +
              team_training_facilities + (1 | team_id),
            data = player_data, REML = FALSE)
mm5 <- lmer(performance_score ~ previous_season_score + fitness_level + goals_scored +
              team_training_facilities + games_played + (1 | team_id),
            data = player_data, REML = FALSE)

tab_model(mm0, mm1, mm2, mm3, mm4, mm5)
  performance score performance score performance score performance score performance score performance score
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p
(Intercept) 49.27 47.81 – 50.73 <0.001 15.39 11.82 – 18.97 <0.001 14.03 10.86 – 17.20 <0.001 11.86 8.46 – 15.26 <0.001 12.84 9.02 – 16.65 <0.001 10.69 6.57 – 14.82 <0.001
previous season score 0.62 0.56 – 0.68 <0.001 0.65 0.59 – 0.70 <0.001 0.65 0.60 – 0.70 <0.001 0.65 0.60 – 0.70 <0.001 0.65 0.60 – 0.70 <0.001
fitness level 2.18 1.77 – 2.58 <0.001 2.17 1.77 – 2.57 <0.001 2.16 1.76 – 2.56 <0.001 2.19 1.79 – 2.59 <0.001
goals scored 0.23 0.09 – 0.37 0.002 0.23 0.09 – 0.37 0.001 0.23 0.09 – 0.37 0.001
team training facilities
[Medium] 		-0.87 -3.79 – 2.06 0.559 -0.82 -3.67 – 2.03 0.571
team training facilities
[High] 		-1.93 -4.86 – 1.00 0.196 -1.90 -4.75 – 0.94 0.189
games played 0.14 0.03 – 0.25 0.012
Random Effects
σ2 43.25 20.01 15.13 14.72 14.72 14.49
τ00 7.77 team_id 7.03 team_id 6.70 team_id 6.52 team_id 5.89 team_id 5.55 team_id
ICC 0.15 0.26 0.31 0.31 0.29 0.28
N 18 team_id 18 team_id 18 team_id 18 team_id 18 team_id 18 team_id
Observations 360 360 360 360 360 360
Marginal R2 / Conditional R2 0.000 / 0.152 0.470 / 0.608 0.575 / 0.705 0.585 / 0.712 0.603 / 0.716 0.611 / 0.719

The likelihood-ratio comparisons below compare nested models; they do not prove that added predictors cause performance changes.

mixed_comparison <- anova(mm0, mm1, mm2, mm3, mm4, mm5)
mixed_comparison
## Data: player_data
## Models:
## mm0: performance_score ~ 1 + (1 | team_id)
## mm1: performance_score ~ previous_season_score + (1 | team_id)
## mm2: performance_score ~ previous_season_score + fitness_level + (1 | team_id)
## mm3: performance_score ~ previous_season_score + fitness_level + goals_scored + (1 | team_id)
## mm4: performance_score ~ previous_season_score + fitness_level + goals_scored + team_training_facilities + (1 | team_id)
## mm5: performance_score ~ previous_season_score + fitness_level + goals_scored + team_training_facilities + games_played + (1 | team_id)
##     npar    AIC    BIC  logLik -2*log(L)    Chisq Df Pr(>Chisq)    
## mm0    3 2411.2 2422.8 -1202.6    2405.2                           
## mm1    4 2145.7 2161.3 -1068.9    2137.7 267.4587  1  < 2.2e-16 ***
## mm2    5 2050.9 2070.3 -1020.5    2040.9  96.8135  1  < 2.2e-16 ***
## mm3    6 2042.9 2066.2 -1015.4    2030.9  10.0216  1   0.001547 ** 
## mm4    8 2045.3 2076.4 -1014.6    2029.3   1.6117  2   0.446714    
## mm5    9 2041.0 2075.9 -1011.5    2023.0   6.3150  1   0.011972 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6. Compare random-effects structures

Using the mm4 fixed effects, compare a random intercept, a fitness random slope without a random intercept, and a model with both. A random slope lets the fitness association vary across teams.

mm4s <- lmer(
  performance_score ~ previous_season_score + fitness_level + goals_scored +
    team_training_facilities + (0 + fitness_level | team_id),
  data = player_data, REML = FALSE
)

mm4b <- lmer(
  performance_score ~ previous_season_score + fitness_level + goals_scored +
    team_training_facilities + (1 + fitness_level | team_id),
  data = player_data, REML = FALSE
)

anova(mm4, mm4s, mm4b)
## Data: player_data
## Models:
## mm4: performance_score ~ previous_season_score + fitness_level + goals_scored + team_training_facilities + (1 | team_id)
## mm4s: performance_score ~ previous_season_score + fitness_level + goals_scored + team_training_facilities + (0 + fitness_level | team_id)
## mm4b: performance_score ~ previous_season_score + fitness_level + goals_scored + team_training_facilities + (1 + fitness_level | team_id)
##      npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)    
## mm4     8 2045.3 2076.4 -1014.6    2029.3                         
## mm4s    8 2125.0 2156.1 -1054.5    2109.0  0.000  0               
## mm4b   10 2047.8 2086.7 -1013.9    2027.8 81.239  2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(mm4, mm4s, mm4b)
  performance score performance score performance score
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 12.84 9.02 – 16.65 <0.001 12.90 9.15 – 16.65 <0.001 12.84 9.04 – 16.64 <0.001
previous season score 0.65 0.60 – 0.70 <0.001 0.64 0.58 – 0.70 <0.001 0.65 0.60 – 0.70 <0.001
fitness level 2.16 1.76 – 2.56 <0.001 2.17 1.61 – 2.73 <0.001 2.17 1.69 – 2.65 <0.001
goals scored 0.23 0.09 – 0.37 0.001 0.29 0.13 – 0.46 0.001 0.24 0.10 – 0.38 0.001
team training facilities
[Medium] 		-0.87 -3.79 – 2.06 0.559 -0.85 -2.00 – 0.31 0.150 -0.87 -3.82 – 2.07 0.559
team training facilities
[High] 		-1.93 -4.86 – 1.00 0.196 -1.89 -3.05 – -0.73 0.001 -1.91 -4.85 – 1.03 0.203
Random Effects
σ2 14.72 20.11 14.36
τ00 5.89 team_id   5.97 team_id
τ11   0.45 team_id.fitness_level 0.33 team_id.fitness_level
ρ01     0.03 team_id
ICC 0.29 0.02 0.31
N 18 team_id 18 team_id 18 team_id
Observations 360 360 360
Marginal R2 / Conditional R2 0.603 / 0.716 0.597 / 0.607 0.601 / 0.723

The best structure depends on fit, complexity, convergence, and whether estimates make sense. Boundary or singular-fit warnings can mean that the demo data do not support all random-effect terms; do not ignore such warnings in a real analysis.

7. Check model diagnostics

Inspect diagnostics for m4 and mm4b. These plots check residual patterns and influential observations. They are aids for judgement, not automatic pass/fail tests.

performance::check_model(m4)

performance::check_model(mm4b)

8. Inspect team-specific fitness slopes

With random intercepts and slopes, coef(mm4b) combines overall fixed effects with each team’s estimated deviation. The dashed line shows the overall fitness coefficient; points show estimated team-specific slopes.

team_slopes <- data.frame(
  team_id = rownames(coef(mm4b)$team_id),
  slope = coef(mm4b)$team_id$fitness_level
)

ggplot(team_slopes, aes(x = slope, y = reorder(team_id, slope))) +
  geom_vline(xintercept = fixef(mm4b)[["fitness_level"]],
             linetype = "dashed", colour = "grey40") +
  geom_point(size = 2.5, colour = "#2C6E9B") +
  labs(x = "Estimated fitness association", y = "Team") +
  theme_minimal()

What to take away

  • The staged model progression compares predictor sets and team random-effects structures using generated demo data.
  • Demo results depend on the simulated player records and values.
  • A random intercept accounts for shared team baselines; a random slope allows a predictor association to vary by team.
  • Treat model coefficients as associations unless the study design supports causal claims.

Reproduce this resource

  1. Install R and RStudio.
  2. Create an R Markdown document and paste this page into it.
  3. Install the packages listed above once using the Console.
  4. Run chunks from top to bottom, then click Knit.
  5. In the HTML preview, click Publish and choose RPubs.