knitr::opts_chunk$set(echo = TRUE)

1. Introduction

Football is a sport increasingly driven by data. While traditional statistics such as goals and assists provide useful information about player performance, they do not always tell the full story. One increasingly popular metric is Expected Goals (xG), which estimates the probability of a shot resulting in a goal based on factors such as its location, angle and type.

But what makes some players generate more expected goals than others? Is it simply the number of shots they take? Does their playing position matter? And do these relationships differ across Europe’s major football leagues?

Regression analysis allows us to investigate these questions by examining relationships between variables. It helps us understand how an outcome changes in relation to one or more predictors.

This analysis uses player statistics from Europe’s top five football leagues between 2017/18 and 2023/24 to investigate the factors associated with expected goals per 90 minutes.

Three regression approaches will be explored:

Simple Linear Regression: Examining the relationship between expected goals and a single predictor.
Multiple Linear Regression: Introducing additional variables to explain differences in expected goals.
Interaction Effects: Investigating whether relationships between variables change depending on player characteristics.

Finally, the models will be compared using statistical measures such as R-squared, AIC, BIC and Root Mean Squared Error (RMSE) to evaluate their performance.

The objective is to demonstrate how increasingly complex regression models can provide a more detailed understanding of football player performance, while recognising that statistical associations do not necessarily imply causation.

2. The Data

The dataset contains individual player statistics from Europe’s five major football leagues: the English Premier League, Spanish La Liga, Italian Serie A, German Bundesliga and French Ligue 1.

Seven seasons of data are included, covering 2017/18 to 2023/24. Each observation represents a player’s statistics for a particular club during a season.

The dataset contains a wide range of performance statistics, including goals, assists, shots, passing accuracy, defensive contributions and expected goals.

Our main outcome variable is Expected Goals per 90 minutes (xG/90). Unlike total expected goals, xG/90 accounts for differences in playing time, making comparisons between players more meaningful.

2.1 Setting Up R

Before importing the data, we need to load several R packages:

tidyverse: Used for data preparation, manipulation and visualisation.
lme4: Used to estimate multilevel regression models.
performance: Used to evaluate model assumptions and compare model performance.
sjPlot: Used to present regression results in tables.

The packages can be installed using install.packages() if they are not already available.

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

2.2 Importing the Data

The seven seasons are stored in separate CSV files. We import each file using read.csv(), before combining them into one dataset using bind_rows().

# Importing player statistics for each season
data_2017 <- read.csv("archive/cleaned_2017-18.csv")
data_2018 <- read.csv("archive/cleaned_2018-19.csv")
data_2019 <- read.csv("archive/cleaned_2019-20.csv")
data_2020 <- read.csv("archive/cleaned_2020-21.csv")
data_2021 <- read.csv("archive/cleaned_2021-22.csv")
data_2022 <- read.csv("archive/cleaned_2022-23.csv")
data_2023 <- read.csv("archive/cleaned_2023-24.csv")

# Combining all seven seasons
Data <- bind_rows(
  data_2017, data_2018, data_2019,
  data_2020, data_2021, data_2022,
  data_2023
)

The bind_rows() function combines datasets vertically, matching variables by their column names. This allows us to analyse player performances across all seven seasons rather than examining each season individually.

2.3 Examining the Dataset

Before conducting any analysis, it is important to understand the structure of the dataset.

The glimpse() function displays the number of observations, variables and their data types. We can also use table() to examine how observations are distributed across seasons and competitions.

# Examine the dataset structure
glimpse(Data)

## Rows: 18,243
## Columns: 65
## $ rk                           <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13…
## $ player                       <chr> "Patrick van Aanholt", "Rolando Aarons", …
## $ nation                       <chr> "Netherlands", "England", "England", "Ita…
## $ pos                          <chr> "DF", "MF", "MF", "DF", "DF", "DF", "MF",…
## $ squad                        <chr> "Crystal Palace", "Newcastle Utd", "Hella…
## $ comp                         <chr> "Premier League", "Premier League", "Seri…
## $ age                          <int> 26, 21, 21, 30, 27, 19, 24, 31, 19, 27, 2…
## $ born                         <int> 1990, 1995, 1995, 1986, 1989, 1997, 1992,…
## $ Matches.Played               <int> 28, 4, 11, 17, 8, 1, 19, 27, 31, 12, 38, …
## $ Avg.Mins.per.Match           <int> 2184, 139, 517, 1057, 499, 77, 1176, 2302…
## $ Goals                        <int> 5, 0, 0, 1, 0, 0, 0, 0, 5, 0, 0, 1, 0, 0,…
## $ Assists                      <int> 1, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 0, 0, 2,…
## $ Goals...Assists              <int> 6, 0, 0, 1, 0, 0, 1, 2, 6, 0, 0, 1, 0, 2,…
## $ Non.Penalty.Goals            <int> 5, 0, 0, 1, 0, 0, 0, 0, 5, 0, 0, 1, 0, 0,…
## $ Penalty.Kicks.Made           <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ Expected.Goals               <dbl> 3.1, 0.1, 0.2, 0.2, 0.1, 0.0, 1.0, 0.4, 6…
## $ Exp.NPG                      <dbl> 3.1, 0.1, 0.2, 0.2, 0.1, 0.0, 1.0, 0.4, 6…
## $ Progressive.Carries          <int> 46, 7, 9, 20, 2, 2, 20, 17, 19, 8, 35, 24…
## $ Progressive.Passes           <int> 92, 3, 17, 81, 12, 1, 49, 103, 20, 38, 14…
## $ Goals.p.90                   <dbl> 0.21, 0.00, 0.00, 0.09, 0.00, 0.00, 0.00,…
## $ Assists.p.90                 <dbl> 0.04, 0.00, 0.00, 0.00, 0.00, 0.00, 0.08,…
## $ Tackles.attempted            <int> 47, 4, 13, 20, 7, 2, 23, 40, 10, 35, 51, …
## $ Tackles.Won                  <int> 32, 4, 8, 17, 3, 1, 15, 32, 8, 32, 32, 6,…
## $ X..Dribbles.tackled          <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ Shots.blocked                <int> 5, 0, 1, 2, 5, 0, 6, 15, 1, 5, 21, 1, 0, …
## $ Passes.blocked               <int> 19, 3, 2, 20, 0, 1, 5, 14, 7, 17, 21, 7, …
## $ Interceptions                <int> 47, 1, 2, 8, 4, 3, 18, 40, 1, 19, 88, 16,…
## $ Clearances                   <int> 64, 0, 0, 29, 20, 3, 25, 141, 11, 20, 255…
## $ Errors.made                  <int> 2, 0, 0, 0, 0, 0, 2, 1, 0, 0, 1, 2, 2, 2,…
## $ Goals.Against                <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ Goals.against.p.90           <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ Saves                        <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ Saves..                      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ Clean.Sheets                 <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ X..Clean.sheets              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ X..Penalty.saves             <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ Passes.Completed             <int> 884, 29, 87, 625, 310, 29, 513, 1088, 206…
## $ Passes.Attempted             <int> 1176, 44, 120, 776, 333, 38, 605, 1321, 3…
## $ Pass.completion..            <dbl> 75.2, 65.9, 72.5, 80.5, 93.1, 76.3, 84.8,…
## $ Progressive.passes.distance  <int> 6422, 77, 325, 4535, 1557, 167, 2345, 847…
## $ X..Short.pass.completed      <dbl> 90.0, 80.0, 83.3, 92.9, 98.4, 89.5, 87.5,…
## $ X..Medium.passes.completed   <dbl> 72.2, 40.0, 71.4, 86.2, 95.5, 83.3, 89.4,…
## $ X..Long.passes.completed     <dbl> 48.7, 66.7, 50.0, 48.3, 70.6, 40.0, 71.0,…
## $ Key.passes                   <int> 18, 0, 3, 10, 0, 1, 7, 7, 12, 3, 10, 6, 6…
## $ X1.3                         <int> 63, 2, 8, 55, 8, 0, 33, 77, 16, 22, 106, …
## $ Passes.into.penalty.area     <int> 28, 1, 7, 20, 0, 1, 1, 3, 6, 6, 2, 2, 7, …
## $ touches_def_pen              <int> 90, 0, 3, 35, 27, 4, 30, 139, 12, 23, 299…
## $ Take.ons.attempted           <int> 31, 6, 15, 1, 1, 1, 20, 10, 47, 5, 23, 24…
## $ X..Successful.take.ons       <dbl> 83.9, 66.7, 53.3, 0.0, 100.0, 100.0, 70.0…
## $ Times.tackled.during.take.on <dbl> 16.1, 33.3, 46.7, 100.0, 0.0, 0.0, 30.0, …
## $ carries_prgc                 <int> 46, 7, 9, 20, 2, 2, 20, 17, 19, 8, 35, 24…
## $ carries.final.3rd            <int> 33, 3, 11, 18, 1, 2, 13, 5, 16, 7, 38, 18…
## $ carries.penalty.area         <int> 7, 3, 2, 0, 0, 0, 0, 0, 8, 1, 0, 2, 2, 0,…
## $ Possessions.lost             <int> 18, 6, 13, 6, 1, 0, 21, 5, 51, 19, 9, 14,…
## $ Goals.Scored                 <int> 5, 0, 0, 1, 0, 0, 0, 0, 5, 0, 0, 1, 0, 0,…
## $ Total.Shots                  <int> 33, 2, 3, 4, 2, 0, 21, 7, 42, 5, 15, 9, 9…
## $ X..Shots.on.target           <dbl> 33.3, 0.0, 0.0, 50.0, 50.0, 0.0, 28.6, 14…
## $ Shots.p.90                   <dbl> 0.45, 0.00, 0.00, 0.17, 0.18, 0.00, 0.46,…
## $ Goals.per.shot               <dbl> 0.15, 0.00, 0.00, 0.25, 0.00, 0.00, 0.00,…
## $ Goals.per.shot.on.target     <dbl> 0.45, 0.00, 0.00, 0.50, 0.00, 0.00, 0.00,…
## $ X..Aerial.Duels.won          <dbl> 54.5, 25.0, 37.5, 55.6, 50.0, 33.3, 44.8,…
## $ Shot.creating.actions.p.90   <dbl> 1.90, 0.65, 0.87, 2.30, 0.36, 2.34, 1.46,…
## $ Goal.creating.actions.p.90   <dbl> 0.16, 0.00, 0.00, 0.34, 0.00, 0.00, 0.08,…
## $ Crosses.Stopped              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ season                       <chr> "2017-2018", "2017-2018", "2017-2018", "2…

# Number of observations and variables
dim(Data)

## [1] 18243    65

# Observations by season
table(Data$season)

## 
## 2017-2018 2018-2019 2019-2020 2020-2021 2021-2022 2022-2023 2023-2024 
##      2549      2503      2558      2634      2689      2687      2623

# Observations by league
table(Data$comp)

## 
##     Bundesliga        La Liga        Ligue 1 Premier League        Serie A 
##           3255           3857           3709           3506           3916

# Check for missing values
sum(is.na(Data))

## [1] 0

These checks help identify potential inconsistencies between seasons and whether any missing values need to be addressed before modelling.

2.4 Preparing the Data

Players accumulate different amounts of playing time throughout a season. To ensure more meaningful comparisons, this analysis focuses on outfield players who have played at least 1,000 minutes, equivalent to approximately 11 complete matches.

This reduces the influence of players with limited appearances, whose per-90 statistics may be less representative of their typical performance. Goalkeepers are also excluded because their responsibilities and attacking statistics differ considerably from those of outfield players.

# Filter players with at least 1,000 minutes
Data <- Data |>
  filter(`Avg.Mins.per.Match` >= 1000, pos != "GK")

Although the variable is named Avg Mins per Match, it represents total minutes played in the supplied data.

Next, we calculate expected goals per 90 minutes:

\[xG/90 = \frac{\text{Total Expected Goals}}{\text{Minutes Played}} \times 90\]

# Calculate expected goals per 90 minutes
Data <- Data |>
  mutate(xg_per90 = Expected.Goals /Avg.Mins.per.Match * 90)

# Summary of expected goals per 90
summary(Data$xg_per90)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.03666 0.07692 0.13807 0.19030 1.15608

# Transforming Carries into the penalty area into per 90 minutes
Data <- Data |>
  mutate(carries_into_penalty_per90 = carries.penalty.area/Avg.Mins.per.Match * 90)

We now have a dataset containing outfield players with at least 1,000 minutes of playing time, alongside a standardised measure of their expected goals.

This provides the foundation for our exploratory analysis and regression models.

3. Exploring the Data

Before building regression models, we need to understand the data. This process is known as Exploratory Data Analysis (EDA) and involves examining distributions, identifying unusual observations and investigating relationships between variables.

Our main outcome is Expected Goals per 90 minutes (xG/90). We will explore how this varies across playing positions and leagues, before examining its relationship with potential predictors.

3.1 Descriptive Statistics

We can begin by using summary() to examine the distribution of xG/90.

summary(Data$xg_per90)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.03666 0.07692 0.13807 0.19030 1.15608

sd(Data$xg_per90, na.rm = TRUE)

## [1] 0.1494009

The summary() function provides the minimum, maximum, mean, median and quartiles, while sd() calculates the standard deviation.

Together, these statistics help us understand the typical xG/90 of a player and how much variation exists across the dataset.

We can also compare expected goals between playing positions.

Data |>
  group_by(pos) |>
  summarise(
    Players = n(),
    Mean_xG = mean(xg_per90, na.rm = TRUE),
    Median_xG = median(xg_per90, na.rm = TRUE),
    SD_xG = sd(xg_per90, na.rm = TRUE))

## # A tibble: 3 × 5
##   pos   Players Mean_xG Median_xG  SD_xG
##   <chr>   <int>   <dbl>     <dbl>  <dbl>
## 1 DF       4216  0.0491    0.0404 0.0390
## 2 FW       2245  0.337     0.314  0.164 
## 3 MF       3394  0.117     0.0919 0.0934

Here, group_by() separates players into positional categories, while summarise() calculates statistics for each group.

3.2 Visualising Expected Goals

A histogram allows us to examine the distribution of expected goals per 90 minutes.

ggplot(Data, aes(x = xg_per90)) +
  geom_histogram(bins = 100, fill = "steelblue", colour = "white") +
  labs(title = "Distribution of Expected Goals per 90 Minutes",
    x = "Expected Goals per 90", y = "Number of Player-Season Observations") +
  theme_classic()

The histogram shows whether xG/90 values are concentrated around a typical value or spread across a wider range. Because attacking opportunities are not distributed equally between players, the distribution is right-skewed, with many players recording relatively low xG/90 and a smaller number recording high values.

3.3 Comparing Playing Positions

A boxplot allows us to compare expected goals across different playing positions.

ggplot(Data, aes(x = pos, y = xg_per90, fill = pos)) +
  geom_boxplot() +
  labs(title = "Expected Goals per 90 by Playing Position",
    x = "Playing Position", y = "Expected Goals per 90") +
  theme_classic() +
  theme(legend.position = "none")

The box represents the middle 50% of observations, while the horizontal line inside each box represents the median. Points beyond the whiskers indicate potentially unusual observations.

This comparison is important because forwards, midfielders and defenders perform different tactical roles, which influences their expected goals.

3.4 Comparing European Leagues

We can also investigate whether expected goals differ between Europe’s five major leagues.

ggplot(Data, aes(x = comp, y = xg_per90, fill = comp)) +
  geom_boxplot() +
  labs(title = "Expected Goals per 90 Across European Leagues",
    x = "League", y = "Expected Goals per 90") +
  theme_classic() +
  theme(legend.position = "none", axis.text.x = element_text(angle = 20, hjust = 1))

3.5 Relationships Between Variables

One potential explanation for differences in expected goals is shot frequency. Players who attempt more shots may generate more expected goals.

We can examine this relationship using a scatterplot.

ggplot(Data, aes(x = `Shots.p.90`, y = xg_per90)) +
  geom_point(alpha = 0.33) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(title = "Relationship Between Shots and Expected Goals",
    x = "Shots per 90 Minutes", y = "Expected Goals per 90") +
  theme_classic()

## `geom_smooth()` using formula = 'y ~ x'

Each point represents a player-season observation, while geom_smooth(method = "lm") adds a linear regression line.

This relationship is expected because xG is calculated from individual shooting opportunities. However, shot frequency alone does not account for the quality of those opportunities.

3.6 Correlation Analysis

Correlation measures the strength and direction of a linear relationship between two numerical variables.

The Pearson correlation coefficient ranges from -1 to +1:

+1: A perfect positive linear relationship.
0: No linear relationship.
-1: A perfect negative linear relationship.

We can calculate correlations between xG/90 and several potential predictors.

# Select variables for correlation analysis
cor_data <- Data |>
  select(xg_per90, `Shots.p.90`, `Goals.p.90`, `Assists.p.90`, `Shot.creating.actions.p.90`)

# Calculate correlation matrix
cor_matrix <- cor(cor_data, use = "pairwise.complete.obs")

round(cor_matrix, 2)

##                            xg_per90 Shots.p.90 Goals.p.90 Assists.p.90
## xg_per90                       1.00       0.89       0.91         0.40
## Shots.p.90                     0.89       1.00       0.86         0.46
## Goals.p.90                     0.91       0.86       1.00         0.39
## Assists.p.90                   0.40       0.46       0.39         1.00
## Shot.creating.actions.p.90     0.45       0.58       0.44         0.68
##                            Shot.creating.actions.p.90
## xg_per90                                         0.45
## Shots.p.90                                       0.58
## Goals.p.90                                       0.44
## Assists.p.90                                     0.68
## Shot.creating.actions.p.90                       1.00

To make these results easier to interpret, we can visualise their correlations with xG/90.

# Create a dataframe of correlations with xG/90
cor_df <- data.frame(
  Variable = names(cor_matrix[, "xg_per90"]),
  Correlation = cor_matrix[, "xg_per90"]) |>
  filter(Variable != "xg_per90")

# correlations plot
ggplot(cor_df, aes(x = reorder(Variable, Correlation),
  y = Correlation)) +
  geom_col(fill = "steelblue") +
  coord_flip() +
  labs(title = "Correlations with Expected Goals per 90", x = "", y = "") +
  theme_classic()

Variables with larger positive correlations tend to increase alongside xG/90, while negative correlations indicate an inverse relationship.

It is important to recognise that correlation does not imply causation. Furthermore, some football statistics are closely related by definition. For example, goals are connected to expected goals, so the strong correlation is not necessarily surprising.

These findings provide a starting point for selecting predictors and developing our regression models.

4. Simple Linear Regression

Our first question is: Do football players who attempt more shots generate higher expected goals?

Simple linear regression examines the relationship between one predictor variable and one outcome variable. In this case:

Outcome (Y): Expected Goals per 90 minutes (xG/90).
Predictor (X): Shots per 90 minutes.

A simple linear regression model can be expressed as:

\[xG/90 = \beta_0 + \beta_1(\text{Shots per 90}) + \varepsilon\]

The intercept (\(\beta_0\)) represents the predicted xG/90 when a player takes zero shots, while the slope (\(\beta_1\)) represents the expected change in xG/90 for each additional shot per 90 minutes.

4.1 Building the Model

In R, we use the lm() function to fit a linear regression model. The formula follows the structure outcome ~ predictor.

# Model 1: Simple Linear Regression
m1 <- lm(xg_per90 ~ Shots.p.90, data = Data)

# Examine regression results
summary(m1)

## 
## Call:
## lm(formula = xg_per90 ~ Shots.p.90, data = Data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34042 -0.03335 -0.00242  0.02450  0.53721 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.002006   0.001012  -1.982   0.0475 *  
## Shots.p.90   0.338668   0.001778 190.496   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06904 on 9853 degrees of freedom
## Multiple R-squared:  0.7865, Adjusted R-squared:  0.7864 
## F-statistic: 3.629e+04 on 1 and 9853 DF,  p-value: < 2.2e-16

4.2 Interpreting the Regression Results

We can also use tab_model() from the sjPlot package to present the regression coefficients in a more readable format.

# Display regression coefficients and confidence intervals
tab_model(m1)

	xg per 90
Predictors	Estimates	CI	p
(Intercept)	-0.00	-0.00 – -0.00	0.048
Shots p 90	0.34	0.34 – 0.34	<0.001
Observations	9855
R² / R² adjusted	0.786 / 0.786

# Calculate 95% confidence intervals
confint(m1)

##                    2.5 %        97.5 %
## (Intercept) -0.003989815 -2.189041e-05
## Shots.p.90   0.335182733  3.421525e-01

4.2 Interpreting the Regression Results

The simple linear regression identified a strong positive relationship between shots per 90 minutes and expected goals per 90 minutes.

The estimated regression equation is:

\[\widehat{xG/90} = -0.002 + 0.339(\text{Shots per 90})\]

The coefficient for shots per 90 was 0.339 (p < 0.001), indicating that each additional shot attempted per 90 minutes is associated with an increase of approximately 0.339 expected goals per 90 minutes.

The model produced an R-squared of 0.787, meaning that shot frequency alone explains approximately 78.7% of the observed variation in xG/90 across the 9,855 player-season observations.

The intercept was -0.002, representing the predicted xG/90 for a player attempting zero shots. Although statistically significant at the 5% level, its magnitude is extremely small and has limited practical importance.

Overall, these results suggest that shot frequency is strongly associated with expected goals. However, this relationship is partly expected because xG is calculated from individual shooting opportunities. The model also does not account for differences in shot quality, playing position or other player characteristics.

5. Multiple Linear Regression

Our simple linear regression found that shots per 90 minutes explained approximately 78.7% of the variation in expected goals per 90 minutes. However, shot frequency alone does not capture every aspect of a player’s attacking performance.

Multiple linear regression allows us to include several predictors in the same model. This helps us investigate whether additional variables are associated with expected goals after accounting for the other predictors.

For this model, we will include:

Shots per 90: How frequently a player attempts shots.
Shot-creating actions per 90: Actions that contribute to creating shooting opportunities.
Carries into the penalty area per 90: How frequently a player carries the ball into the opposition’s penalty area.
Playing position: Whether forwards, midfielders and defenders generate different levels of expected goals.

A multiple linear regression model can be expressed as:

\[Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_3 + \cdots + \varepsilon\]

Each coefficient represents the estimated relationship between a predictor and the outcome, holding the other predictors constant.

5.1 Building the Model

We can add additional predictors to the lm() function using the + symbol.

Before fitting the model, we convert playing position into a factor. This tells R to treat position as a categorical variable rather than a numerical measurement.

# Convert playing position into a categorical variable
Data <- Data |>
  mutate(pos = factor(pos))

# Model 2: Multiple Linear Regression
m2 <- lm(xg_per90 ~ Shots.p.90 + Shot.creating.actions.p.90 + carries_into_penalty_per90 + pos, data = Data)

# Examine regression results
summary(m2)

## 
## Call:
## lm(formula = xg_per90 ~ Shots.p.90 + Shot.creating.actions.p.90 + 
##     carries_into_penalty_per90 + pos, data = Data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34354 -0.02927 -0.00340  0.02303  0.58244 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 0.0123305  0.0013566   9.089  < 2e-16 ***
## Shots.p.90                  0.3055093  0.0028767 106.200  < 2e-16 ***
## Shot.creating.actions.p.90 -0.0066971  0.0007763  -8.627  < 2e-16 ***
## carries_into_penalty_per90 -0.0097358  0.0017718  -5.495 4.01e-08 ***
## posFW                       0.0722780  0.0026633  27.139  < 2e-16 ***
## posMF                       0.0052953  0.0017876   2.962  0.00306 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06545 on 9849 degrees of freedom
## Multiple R-squared:  0.8082, Adjusted R-squared:  0.8081 
## F-statistic:  8299 on 5 and 9849 DF,  p-value: < 2.2e-16

5.2 Interpreting the Regression Results

The multiple regression model explained 80.8% of the variation in expected goals per 90 minutes, compared with 78.7% for the simple linear regression. This represents an improvement of approximately 2.2 percentage points after introducing additional predictors.

We can use tab_model() from the sjPlot package to compare the regression coefficients from both models.

# Compare regression coefficients
tab_model(m1, m2)

	xg per 90			xg per 90
Predictors	Estimates	CI	p	Estimates	CI	p
(Intercept)	-0.00	-0.00 – -0.00	0.048	0.01	0.01 – 0.01	<0.001
Shots p 90	0.34	0.34 – 0.34	<0.001	0.31	0.30 – 0.31	<0.001
Shot creating actions p 90				-0.01	-0.01 – -0.01	<0.001
carries into penalty per90				-0.01	-0.01 – -0.01	<0.001
pos [FW]				0.07	0.07 – 0.08	<0.001
pos [MF]				0.01	0.00 – 0.01	0.003
Observations	9855			9855
R² / R² adjusted	0.786 / 0.786			0.808 / 0.808

The coefficient for shots per 90 decreased from 0.339 in Model 1 to 0.306 (p < 0.001) in Model 2. This indicates that each additional shot per 90 minutes is associated with approximately 0.306 additional expected goals per 90, holding the other predictors constant.

Playing position was also significantly associated with expected goals. Compared with defenders, forwards generated approximately 0.072 more xG per 90 (p < 0.001) after accounting for the other variables. Midfielders generated approximately 0.005 more xG per 90 (p = 0.003).

Interestingly, two attacking statistics produced negative coefficients:

Shot-creating actions per 90: -0.0067 (p < 0.001).
Carries into the penalty area per 90: -0.0097 (p < 0.001).

These negative coefficients do not necessarily suggest that creating shooting opportunities or carrying the ball into the penalty area reduces expected goals. Instead, they describe the relationship between each variable and xG/90 when the other predictors are held constant.

For example, players who frequently create opportunities for teammates may have different attacking responsibilities from players who primarily attempt shots themselves, indicating that the jobs are somewhat independent of each other.

This demonstrates an important feature of multiple regression: a variable’s relationship with the outcome can change after accounting for other predictors.

5.3 Comparing Model Performance

Adding predictors can increase a model’s R-squared, but this does not automatically mean the model is substantially better.

We can compare our two models using compare_performance() from the performance package.

# Compare model performance
compare_performance(m1, m2)

## # Comparison of Model Performance Indices
## 
## Name | Model |    AIC (weights) |   AICc (weights) |    BIC (weights) |    R2
## -----------------------------------------------------------------------------
## m1   |    lm | -24714.4 (<.001) | -24714.4 (<.001) | -24692.8 (<.001) | 0.786
## m2   |    lm | -25762.9 (>.999) | -25762.9 (>.999) | -25712.5 (>.999) | 0.808
## 
## Name | R2 (adj.) |  RMSE | Sigma
## --------------------------------
## m1   |     0.786 | 0.069 | 0.069
## m2   |     0.808 | 0.065 | 0.065

6. Interaction Effects

Our previous models assumed that each additional shot was associated with the same increase in expected goals, regardless of playing position.

We can test whether this relationship differs between defenders, midfielders and forwards by introducing an interaction effect.

6.1 Building an Interaction Model

In R, the * operator includes both variables and their interaction.

# Model 3: Interaction between shots and position
m3 <- lm(xg_per90 ~ Shots.p.90 * pos + Shot.creating.actions.p.90 + carries_into_penalty_per90, data = Data)

# Examine interaction coefficients
summary(m3)

## 
## Call:
## lm(formula = xg_per90 ~ Shots.p.90 * pos + Shot.creating.actions.p.90 + 
##     carries_into_penalty_per90, data = Data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.37162 -0.02683 -0.00629  0.02198  0.58407 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 0.0219182  0.0016981  12.907  < 2e-16 ***
## Shots.p.90                  0.2186484  0.0077983  28.038  < 2e-16 ***
## posFW                       0.0166298  0.0041850   3.974 7.13e-05 ***
## posMF                       0.0004097  0.0025259   0.162    0.871    
## Shot.creating.actions.p.90 -0.0038571  0.0007847  -4.915 9.01e-07 ***
## carries_into_penalty_per90 -0.0093372  0.0017464  -5.347 9.16e-08 ***
## Shots.p.90:posFW            0.1271490  0.0085459  14.878  < 2e-16 ***
## Shots.p.90:posMF            0.0552411  0.0084813   6.513 7.71e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06451 on 9847 degrees of freedom
## Multiple R-squared:  0.8137, Adjusted R-squared:  0.8136 
## F-statistic:  6144 on 7 and 9847 DF,  p-value: < 2.2e-16

The interaction coefficients show how the relationship between shots and expected goals changes relative to defenders, who are the reference category.

6.2 Interpreting Interaction Effects

The interaction model explained 81.4% of the variation in expected goals per 90 minutes, improving on the multiple regression model.

Both interaction terms were statistically significant (p < 0.001), suggesting that the relationship between shots and expected goals differs across playing positions.

The estimated shot coefficients were:

Defenders: 0.219 xG per additional shot per 90.
Midfielders: 0.274 xG per additional shot per 90.
Forwards: 0.346 xG per additional shot per 90.

These results suggest that additional shots are associated with larger increases in expected goals for forwards than for midfielders or defenders, holding the other predictors constant.

6.3 Visualising the Interaction

We can visualise these relationships by plotting separate regression lines for each playing position.

ggplot(Data, aes(x = Shots.p.90, y = xg_per90, colour = pos)) +
  geom_point(alpha = 0.15) +
  geom_smooth(method = "lm", se = FALSE) +
  facet_wrap(~ pos) +
  labs(title = "Shots and Expected Goals by Playing Position",
    x = "Shots per 90 Minutes", y = "Expected Goals per 90 Minutes", colour = "Position") +
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'

6.4 Comparing Model Performance

We can compare all three models and use ANOVA to test whether adding the interaction terms significantly improves model fit.

# Compare all three regression models
compare_performance(m1, m2, m3)

## # Comparison of Model Performance Indices
## 
## Name | Model |    AIC (weights) |   AICc (weights) |    BIC (weights) |    R2
## -----------------------------------------------------------------------------
## m1   |    lm | -24714.4 (<.001) | -24714.4 (<.001) | -24692.8 (<.001) | 0.786
## m2   |    lm | -25762.9 (<.001) | -25762.9 (<.001) | -25712.5 (<.001) | 0.808
## m3   |    lm | -26047.3 (>.999) | -26047.3 (>.999) | -25982.5 (>.999) | 0.814
## 
## Name | R2 (adj.) |  RMSE | Sigma
## --------------------------------
## m1   |     0.786 | 0.069 | 0.069
## m2   |     0.808 | 0.065 | 0.065
## m3   |     0.814 | 0.064 | 0.065

# Test the improvement from interaction effects
anova(m2, m3)

## Analysis of Variance Table
## 
## Model 1: xg_per90 ~ Shots.p.90 + Shot.creating.actions.p.90 + carries_into_penalty_per90 + 
##     pos
## Model 2: xg_per90 ~ Shots.p.90 * pos + Shot.creating.actions.p.90 + carries_into_penalty_per90
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1   9849 42.193                                  
## 2   9847 40.976  2    1.2169 146.22 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Metric	Model 1	Model 2	Model 3
R²	0.786	0.808	0.814
AIC	-24714.4	-25762.9	-26047.3
BIC	-24692.8	-25712.5	-25982.5
RMSE	0.069	0.065	0.064

Model 3 achieved the highest R² and lowest AIC, BIC and RMSE.

The ANOVA test also confirmed that the interaction terms significantly improved model fit (F = 146.22, p < 0.001).

Overall, the results demonstrate that allowing the relationship between shots and expected goals to vary by playing position improves the model, although the increase in explained variation is relatively small.

7. Key Findings

Shot frequency is strongly associated with expected goals. Even the simple linear regression explained approximately 78.7% of the variation in xG/90.
Additional attacking statistics and playing position improve model fit. Multiple regression increased R² to 80.8%, although the improvement was relatively modest.
Shot-creating actions and carries into the penalty area have negative conditional relationships with xG. Both variables produced statistically significant negative coefficients after accounting for shot frequency and playing position. This may reflect differences in attacking responsibilities, where players who frequently create opportunities or carry the ball into dangerous areas are not necessarily those taking the shots themselves. It also suggests that greater attacking involvement does not automatically translate into higher individual expected goals when shot volume is held constant.
The relationship between shots and expected goals differs by position. The interaction model showed that an additional shot per 90 was associated with approximately 0.219 xG/90 for defenders, 0.274 for midfielders and 0.346 for forwards, holding the other predictors constant.

Regression Analysis: Exploring Football Player Performance Across Europe’s Top Five Leagues

Fred Quartermaine

2026-10-12

1. Introduction

2. The Data

2.1 Setting Up R

2.2 Importing the Data

2.3 Examining the Dataset

2.4 Preparing the Data

3. Exploring the Data

3.1 Descriptive Statistics

3.2 Visualising Expected Goals

3.3 Comparing Playing Positions

3.4 Comparing European Leagues

3.5 Relationships Between Variables

3.6 Correlation Analysis

4. Simple Linear Regression

4.1 Building the Model

4.2 Interpreting the Regression Results

4.2 Interpreting the Regression Results

5. Multiple Linear Regression

5.1 Building the Model

5.2 Interpreting the Regression Results

5.3 Comparing Model Performance

6. Interaction Effects

6.1 Building an Interaction Model

6.2 Interpreting Interaction Effects

6.3 Visualising the Interaction

6.4 Comparing Model Performance

7. Key Findings