1. Introduction

Do basic stats such as marks, kicks, handballs and inside 50s have a big impact on score? This tutorial will utilised AFL data from 2023 to 2026 to investigate the relationship between basic stats and score.

2. Set up R

Install these packages if you haven’t already, then load them.

install.packages(c("tidyverse", "fitzRoy", "janitor", "skimr", "GGally",
                   "performance", "sjPlot"))
library(tidyverse)
library(fitzRoy)
library(janitor)
library(GGally)
library(sjPlot)

3. Get AFL Data

To get the data we used the fitzRoy package and collected the player stats for each game between 2023 and 2026. We filtered out the data that we would not use. Always check your data when you import it by using head() or any other function with the same result.

afl_data <-  fetch_player_stats(season = 2023:2026, source = "afltables") |>
  clean_names() |>                       # using janitor package it makes column names consistent, e.g. home_team
  mutate(date = as_date(date))

afl_data |>
  select(season, round, playing_for, kicks, handballs, inside_50s, marks) |>
  head()
## # A data frame: 6 × 7
##   season round playing_for kicks handballs inside_50s marks
## *  <int> <chr> <chr>       <int>     <int>      <int> <int>
## 1   2023 1     Richmond       10         5          1     5
## 2   2023 1     Richmond        6         2          1     3
## 3   2023 1     Richmond       15         3          7     6
## 4   2023 1     Richmond       11         2          2     7
## 5   2023 1     Richmond        9         9          3     3
## 6   2023 1     Richmond       11         9          5     1

4. Make one row for one team for each game

Add up all the player data so each team only has row per game. This will help with our modelling later

team_scores <- afl_data |>
  mutate(round_num = as.numeric(round)) |>
  filter(!is.na(round_num)) |>       # dropping finals, so only regular season games are kept
  mutate(score = if_else(playing_for == home_team, home_score, away_score)) |>
  group_by(season, round, date, team = playing_for, score) |>
  summarise(across(c("kicks", "handballs", "inside_50s", "marks"), ~ sum(.x, na.rm = TRUE)), .groups = "drop")

head(team_scores)
## # A tibble: 6 × 9
##   season round date       team           score kicks handballs inside_50s marks
##    <int> <chr> <date>     <chr>          <int> <int>     <int>      <int> <int>
## 1   2023 1     2023-03-16 Carlton           58   221       120         45    97
## 2   2023 1     2023-03-16 Richmond          58   228       121         66    95
## 3   2023 1     2023-03-17 Collingwood      125   229       143         62   105
## 4   2023 1     2023-03-17 Geelong          103   182       133         46    69
## 5   2023 1     2023-03-18 Brisbane Lions    72   162       102         40    52
## 6   2023 1     2023-03-18 Gold Coast        61   225       141         44    56

It’s important to do sanity checks on your data throughout your modelling. Here we see that each season has 414/2 = 207 games.

team_scores |> count(season)
## # A tibble: 4 × 2
##   season     n
##    <int> <int>
## 1   2023   414
## 2   2024   414
## 3   2025   414
## 4   2026   414

5. Data Exploration

How much are the variables related to each other? Using the ggpairs() function we can see the correlation between the chosen variables. The numbers in the top right are how much they are correlated from -1 to 1 and the bottom left shows the scatter plots.

team_scores |>
  select(score, kicks, handballs, inside_50s, marks) |>
  ggpairs()

You can see that inside 50s has the highest correlation to score by far. We should expect that inside 50s will have the biggest impact.

6. Plotting Variables

It is important to visualise the data that you are going to model. Here we’ll replot the scatterplots from above in a way that is easier to read.

team_scores |>
  pivot_longer(c("kicks", "handballs", "inside_50s", "marks"), names_to = "stat", values_to = "count") |>
  ggplot(aes(count, score)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap(~ stat, scales = "free_x") +
  labs(x = "Count", y = "Team score (points)",
       title = "Team score against each stat") +
  theme_minimal()

Inside 50s had the steepest upwards trend. This backs up the analysis from the correlation plot. It makes the most sense out of these variables as you need an inside 50 to score. By being in an attacking position provides more chances. It is expected that handballs and kicks show shallower trends as the team having possession of the ball doesn’t related to using it effectively.

7. Linear Modelling

Since we expect that inside 50s has the biggest impact our first model will only use inside 50s as a predictor.

lm_1 <- lm(score ~ inside_50s, data = team_scores)

summary(lm_1)
## 
## Call:
## lm(formula = score ~ inside_50s, data = team_scores)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -79.013 -14.345  -0.756  13.827  73.535 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -11.97551    3.28074   -3.65  0.00027 ***
## inside_50s    1.86286    0.06204   30.03  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 20.53 on 1654 degrees of freedom
## Multiple R-squared:  0.3528, Adjusted R-squared:  0.3524 
## F-statistic: 901.6 on 1 and 1654 DF,  p-value: < 2.2e-16

The estimate of 1.86 means that for every inside 50 it is associated with roughly 1.86 more points. With an intercept of -11.98, that means that when inside 50s = 0 then the score is negative. This isn’t meaningful however it’s unrealistic that an AFL team will have 0 inside 50s. Multiple R-squared being 0.3528 means that inside 50s explain roughly 35% of variation in team scores. With a p-value being near zero this mean that it’s significantly unlikely to have these results due to chance. Residual standard error = 20.53 means that predictions are typically off by about 20.5 points, there’s a lot of spread.

Now we’ll model using four predictors and compare it with the previous model.

lm_2 <- lm(score ~ kicks + handballs + inside_50s + marks, data = team_scores)

summary(lm_2)
## 
## Call:
## lm(formula = score ~ kicks + handballs + inside_50s + marks, 
##     data = team_scores)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -78.462 -14.445   0.062  14.049  71.109 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -48.13614    5.96583  -8.069 1.35e-15 ***
## kicks         0.09868    0.03835   2.573   0.0102 *  
## handballs     0.09279    0.02174   4.268 2.08e-05 ***
## inside_50s    1.63542    0.07071  23.129  < 2e-16 ***
## marks         0.14645    0.03695   3.964 7.68e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.92 on 1651 degrees of freedom
## Multiple R-squared:  0.3922, Adjusted R-squared:  0.3907 
## F-statistic: 266.3 on 4 and 1651 DF,  p-value: < 2.2e-16
sjPlot::tab_model(lm_1, lm_2, show.aic = TRUE)
  score score
Predictors Estimates CI p Estimates CI p
(Intercept) -11.98 -18.41 – -5.54 <0.001 -48.14 -59.84 – -36.43 <0.001
inside 50s 1.86 1.74 – 1.98 <0.001 1.64 1.50 – 1.77 <0.001
kicks 0.10 0.02 – 0.17 0.010
handballs 0.09 0.05 – 0.14 <0.001
marks 0.15 0.07 – 0.22 <0.001
Observations 1656 1656
R2 / R2 adjusted 0.353 / 0.352 0.392 / 0.391
AIC 14712.312 14614.414

The R-squared increase to 39.2% so we know that the new variables are genuinely adding information rather than adding more noise. The kicks, handballs and marks effects are quite small, but these stats come in much larger numbers per game, so they add up. For example, 20 extra kicks is worth about 2 points. The t values show that by far inside 50s have the biggest impact and overwhelming evidence of an effect. Handballs and marks show it’s strongly unlikely to be chance.

The lower AIC shows that the four predictor model lm_2 performed better than lm_1