Approach 1: Weighted Percentile Rank

The first approach I’ll present is a handcrafted “ranking” method.

1. Calculate the proportional difference between the pre-target and target season outperformance ratios– and respectively–for all players .

2. Weight by the player’s accumulated in prior seasons.⁴

3. Calculate the the underperforming unlikeliness as a percentile rank of ascending , i.e. more negative values correspond to a lower percentile.⁵⁶

This is pretty straightforward to calculate once you’ve got the data prepared in the right format.

## `u` for "underperforming unlikelihood"
all_u_approach1 <- all_players_to_evaluate |> 
  dplyr::transmute(
    player,
    prior_o,
    target_o,
    prior_xg,
    weighted_delta_o = prior_shots * (target_o - prior_o) / prior_o,
    u = dplyr::percent_rank(weighted_delta_o)
  ) |> 
  dplyr::arrange(u)

maddison_u_approach1 <- all_u_approach1 |> 
  dplyr::filter(player == 'James Maddison')

maddison_u_approach1 |> dplyr::select(player, prior_o, target_o, u)
#> # A tibble: 1 × 4
#>   player         prior_o target_o      u
#>   <chr>            <dbl>    <dbl>  <dbl>
#> 1 James Maddison    1.38    0.847 0.0321

This approach finds Maddison’s 2023/24 of 0.847 to be about a 3rd percentile outcome. Among the 593 players evaluated, Maddison’s 2023/24 ranks as the 20th most unlikely.

For context, here’s a look at the top 10 most unlikely outcomes for the 2023/24 season.

all_u_approach1 |> head(10) |> dplyr::select(player, prior_o, target_o, u)
#> # A tibble: 10 × 4
#>    player             prior_o target_o       u
#>    <chr>                <dbl>    <dbl>   <dbl>
#>  1 Ciro Immobile        1.23     0.383 0      
#>  2 Giovanni Simeone     1.03     0.306 0.00169
#>  3 Nabil Fekir          1.14     0.5   0.00338
#>  4 Wahbi Khazri         1.11     0.322 0.00507
#>  5 Téji Savanier        1.42     0.282 0.00676
#>  6 Adrien Thomasson     1.18     0.282 0.00845
#>  7 Timo Werner          0.951    0.543 0.0101 
#>  8 Gaëtan Laborde       1.02     0.546 0.0118 
#>  9 Kevin Volland        1.18     0.450 0.0135 
#> 10 Robert Lewandowski   1.01     0.759 0.0152 

Ciro Immobile tops the list, with several other notable attacking players who had less than stellar seasons.

Overall, I’d say that this methodology seems to generate fairly reasonable results, but it’s hard to pinpoint why exactly this approach is defensible other than it may lead to intuitive results. Despite the intuitively appealing outcomes, it’s important to scrutinize key factors affecting this methodology’s robustness and validity.

• Subjectivity in Weighting: The choice to weight the difference in performance ratios by pre-2023/24 xG is inherently subjective. While it’s important to have some form of weighting so as to avoid disproportionately emphasizing players with minimal shooting opportunities, lternative weighting strategies could lead to significantly different rankings.
• Sensitivity to Player Pool: The percentile ranking of a player’s unlikeliness is highly sensitive to the comparison group. For instance, comparing forwards to defenders could skew results due to generally higher variability in defenders’ goal-to-expected goals ratios. Moreover, if a season unusually favors players outperforming their expected goals, even average performers might appear as outliers. This potential for selection bias underlines the importance of carefully choosing the comparison pool.
• Assumption of Uniform Distribution: The methodology assumes a uniform distribution of performance unlikeliness across players. We presume that there must be, for example, 1% of players at the 1st percentile of underperformance. Although this assumption may hold over many seasons, it could misrepresent individual seasons where extremes are less pronounced or absent.

Approach 2: Resampling from Prior History of Shots

There’s only so much you can do with player-season-level data. We need to dive into shot-level data if we want to more robustly understand the uncertainty of outcomes.

Here’s a “resampling” approach to quantify the underperforming unlikeliness of a player in the target season:

1. Sample shots (with replacement⁷) from a player’s past shots . Repeat this for resamples.⁸
2. Count the number of resamples in which the outperformance ratio of the sampled shots is less than or equal to the observed in the target season for the player.⁹ The proportion represents the unlikeness of a given player’s observed (or worse) in the target season.

Here’s how that looks in code.

R <- 1000
resample_player_shots <- function(
    shots, 
    n_shots_to_sample, 
    n_sims = R,
    replace = TRUE,
    seed = 42
) {
  
  withr::local_seed(seed)
  purrr::map_dfr(
    1:n_sims,
    \(.sim) {
      sampled_shots <- shots |> 
        slice_sample(n = n_shots_to_sample, replace = replace)
      
      list(
        sim = .sim,
        xg = sum(sampled_shots$xg),
        g = sum(sampled_shots$g),
        o = sum(sampled_shots$g) / sum(sampled_shots$xg)
      )
    }
  )
}

resample_one_player_o <- function(shots, target_season_end_year) {
  target_shots <- shots |>
    dplyr::filter(season_end_year == target_season_end_year)
  
  prior_shots <- shots |>
    dplyr::filter(season_end_year < target_season_end_year)
  
  prior_shots |> 
    resample_player_shots(
      n_shots_to_sample = nrow(target_shots)
    )
}

resample_player_o <- function(shots, players, target_season_end_year = TARGET_SEASON_END_YEAR) {
  purrr::map_dfr(
    players,
    \(.player) {
      shots |> 
        dplyr::filter(player == .player) |> 
        resample_one_player_o(
          target_season_end_year = target_season_end_year
        ) |> 
        dplyr::mutate(
          player = .player
        )
    }
  )
}

maddison_resampled_o <- np_shots |> 
  resample_player_o(
    players = 'James Maddison'
  ) |> 
  dplyr::inner_join(
    wide_player_np_shots |> 
      dplyr::select(
        player,
        prior_o,
        target_o
      ),
    by = dplyr::join_by(player)
  ) |> 
  dplyr::arrange(player)

maddison_u_approach2 <- maddison_resampled_o |>
  dplyr::summarize(
    .by = c(player, prior_o, target_o),
    u = sum(o <= target_o) / n()
  ) |> 
  dplyr::arrange(player)

maddison_u_approach2 |> dplyr::select(player, prior_o, target_o, u)
#> # A tibble: 1 × 4
#>   player         prior_o target_o     u
#>   <chr>            <dbl>    <dbl> <dbl>
#> 1 James Maddison    1.38    0.847 0.163

The plot below should provide a bit of visual intuition as to what’s going on.

These results imply that Maddison’s 2023/24 ratio of 0.847 (or worse) occurs in 16.3% of simulations, i.e. a 16th percentile outcome. That’s a bit higher than what the first approach showed.

How can we feel more confident about this approach? Well, in the first approach, we assumed that the underperforming unlikelihood percentages should be uniform across all players, hence the percentile ranking. I think that’s a good assumption, so we should see if the same bears out with this second approach.

The plot below shows a histogram of the underperforming unlikelihood across all players, where each player’s estimated unlikelihood is grouped into a decile.

all_resampled_o <- np_shots |> 
  resample_player_o(
    players = all_players_to_evaluate$player
  ) |> 
  dplyr::inner_join(
    wide_player_np_shots |> 
      ## to make sure we just one Rodri, Danilo, and Nicolás González 
      dplyr::filter(player_id %in% all_players_to_evaluate$player_id) |> 
      dplyr::select(
        player,
        prior_o,
        target_o,
        prior_shots,
        target_shots
      ),
    by = dplyr::join_by(player)
  ) |> 
  dplyr::arrange(player, player)

all_u_approach2 <- all_resampled_o |>
  dplyr::summarize(
    .by = c(player, prior_o, target_o, prior_shots, target_shots),
    u = sum(o <= target_o) / n()
  ) |> 
  dplyr::arrange(u)

Indeed, the histogram shows a fairly uniform distribution, with a bit of irregularity at the very edges.

Looking at who is in the lower end of the leftmost decile, we see some of the same names–Immobile and Savanier–among the ten underperformers. (Withholding judgment on the superiority of any methodology, we can find some solace in seeing some of the same names among the most unlikely underperformers here as we did with approach 1.)

all_u_approach2 |> head(10) |> dplyr::select(player, prior_o, target_o, u)
#> # A tibble: 10 × 4
#>    player                    prior_o target_o     u
#>    <chr>                       <dbl>    <dbl> <dbl>
#>  1 Erling Haaland               1.26    0.791 0.004
#>  2 Amine Harit                  1.27    0.262 0.01 
#>  3 Pierre-Emerick Aubameyang    1.07    0.606 0.01 
#>  4 Téji Savanier                1.42    0.282 0.01 
#>  5 Antonio Sanabria             1.05    0.386 0.014
#>  6 Alex Baena                   1.54    0.348 0.016
#>  7 Elye Wahi                    1.38    0.753 0.017
#>  8 Kevin Behrens                1.39    0.673 0.019
#>  9 Ciro Immobile                1.23    0.383 0.02 
#> 10 Ansu Fati                    1.31    0.430 0.024

One familiar face in the printout above is Manchester City’s striker Erling Haaland, whose underperformance this season has been called among fans and the media. His sub-par performance this year ranked as a 7th percentile outcome by approach 1, which is very low, but not quite as low as what this approach finds.

Here are some parting thoughts on this methodology before we look at another:

• Assumption of Shot Profile Consistency: We assume that a player’s past shot behavior accurately predicts their future performance. This generally holds unless a player changes their role or team, or is recovering from an injury. Deviations in these contexts need to be considered when applying this methodology.
• Non-Parametric Nature: This method does not assume any specific distribution for a player’s performance ratios; instead, it relies on the stability of a player’s performance over time. The resampling process itself shapes the outcome distribution, which can vary significantly between players with different shooting behaviors, such as a forward versus a defender.
• Computational Demands: The resampling approach requires relatively more computational resources than the prior approach, especially without parallel processing. Even a relatively small number of resamples, such as , can take a few second per player to compute.

Approach 3: Evaluating a Player-Specific Cumulative Distribution Function (CDF)

If we assume that the set of goals-to-xG ratios come from a Gamma data-generating process, then we can leverage the properties of a player-level Gamma distribution to assess the unlikelihood of a players ratio.

To calculate the underperforming unlikeliness :

1. Estimate a Gamma distribution to model a player’s true outperformance ratio across all prior shots, excluding those in the target season–.
2. Calculate the probability that is less than or equal to the player’s observed in the target season using the Gamma distribution’s cumulative distribution function (CDF).

While that may sound daunting, I promise that it’s not (well, aside from a bit of “magic” in estimating a reasonable Gamma distribution per player).

N_SIMS <- 10000

SHOT_TO_SHAPE_MAPPING <- list(
  'from' = c(50, 750),
  'to' = c(1, 25)
)
estimate_one_gamma_distributed_o <- function(
    shots,
    target_season_end_year
) {
  player_np_shots <- shots |> 
    dplyr::mutate(is_target = season_end_year == target_season_end_year)
  
  prior_player_np_shots <- player_np_shots |> 
    dplyr::filter(!is_target)
  
  target_player_np_shots <- player_np_shots |> 
    dplyr::filter(is_target)
  

  agg_player_np_shots <- player_np_shots |>
    dplyr::summarize(
      .by = c(is_target),
      shots = dplyr::n(),
      dplyr::across(c(g, xg), \(.x) sum(.x))
    ) |> 
    dplyr::mutate(o = g / xg)
  
  agg_prior_player_np_shots <- agg_player_np_shots |> 
    dplyr::filter(!is_target)
  
  agg_target_player_np_shots <- agg_player_np_shots |> 
    dplyr::filter(is_target)

  shape <- dplyr::case_when(
    agg_prior_player_np_shots$shots < SHOT_TO_SHAPE_MAPPING$from[1] ~ SHOT_TO_SHAPE_MAPPING$to[2],
    agg_prior_player_np_shots$shots > SHOT_TO_SHAPE_MAPPING$from[2] ~ SHOT_TO_SHAPE_MAPPING$to[2],
    TRUE ~ scales::rescale(
      agg_prior_player_np_shots$shots, 
      from = SHOT_TO_SHAPE_MAPPING$from, 
      to = SHOT_TO_SHAPE_MAPPING$to
    )
  )
  list(
    'shape' = shape,
    'rate' = shape / agg_prior_player_np_shots$o
  )
}

estimate_gamma_distributed_o <- function(
    shots,
    players,
    target_season_end_year
) {
  
  purrr::map_dfr(
    players,
    \(.player) {
      params <- shots |> 
        dplyr::filter(player == .player) |> 
        estimate_one_gamma_distributed_o(
          target_season_end_year = target_season_end_year
        )
      
      list(
        'player' = .player,
        'params' = list(params)
      )
    }
  )
}

select_gamma_o <- np_shots |> 
  estimate_gamma_distributed_o(
    players = 'James Maddison',
    target_season_end_year = TARGET_SEASON_END_YEAR
  ) |> 
  dplyr::inner_join(
    wide_player_np_shots |> 
      dplyr::select(
        player,
        prior_o,
        target_o
      ),
    by = dplyr::join_by(player)
  ) |> 
  dplyr::arrange(player)

maddison_u_approach3 <- select_gamma_o |> 
  dplyr::mutate(
    u = purrr::map2_dbl(
      target_o,
      params,
      \(.target_o, .params) {
        pgamma(
          .target_o, 
          shape = .params$shape, 
          rate = .params$rate,
          lower.tail = TRUE
        )
      }
    ),
    ou = 1 - u
  ) |> 
  tidyr::unnest_wider(params)

maddison_u_approach3 |> dplyr::select(player, prior_o, target_o, u)
#> # A tibble: 1 × 4
#>   player         prior_o target_o     u
#>   <chr>            <dbl>    <dbl>  <dbl>
#> 1 James Maddison    1.38    0.847 0.0679

We see that Maddison’s 2023/24 ratio of 0.847 (or worse) is about a 7th percentile outcome given his prior shot history. The 7th percentile outcome estimated is about not far from the 3rd percentile outcome estimated with approach 1.

To gain some intuition around this approach, we can plot out the Gamma distributed estimate of Maddison’s . The result is a histogram that looks not all that dissimilar to the one from before with resampled shots, just much smoother (since this is a “parametric” approach).

As with approach 2, we should check to see what the distribution of underperforming unlikeliness looks like–we should expect to see a somewhat uniform distribution.

all_gamma_o <- np_shots |> 
  estimate_gamma_distributed_o(
    players = all_players_to_evaluate$player,
    target_season_end_year = TARGET_SEASON_END_YEAR
  ) |> 
  dplyr::inner_join(
    wide_player_np_shots |> 
      dplyr::filter(
        player_id %in% all_players_to_evaluate$player_id
      ) |> 
      dplyr::select(
        player,
        prior_o,
        target_o
      ),
    by = dplyr::join_by(player)
  ) |> 
  dplyr::arrange(player)

all_u_approach3 <- all_gamma_o |> 
  dplyr::mutate(
    u = purrr::map2_dbl(
      target_o,
      params,
      \(.target_o, .params) {
        pgamma(
          .target_o, 
          shape = .params$shape, 
          rate = .params$rate,
          lower.tail = TRUE
        )
      }
    )
  ) |> 
  tidyr::unnest_wider(params) |> 
  dplyr::arrange(u)

This histogram has a bit more distortion than our resampling approach, so perhaps it’s a little less calibrated.

Looking at the top 10 strongest underperformers, 3 of the names here–Immobile, Savanier, and Sanabria–are shared with approach 2’s top 10, and 7 are shared with approach 1’s top 10.

all_u_approach3 |> head(10) |> dplyr::select(player, prior_o, target_o, u)
#> # A tibble: 10 × 4
#>    player           prior_o target_o         u
#>    <chr>              <dbl>    <dbl>      <dbl>
#>  1 Ciro Immobile       1.23    0.383 0.00000533
#>  2 Téji Savanier       1.42    0.282 0.000127  
#>  3 Giovanni Simeone    1.03    0.306 0.000248  
#>  4 Adrien Thomasson    1.18    0.282 0.000346  
#>  5 Wahbi Khazri        1.11    0.322 0.000604  
#>  6 Fabián Ruiz Peña    1.67    0.510 0.00271   
#>  7 Nabil Fekir         1.14    0.5   0.00308   
#>  8 Kevin Volland       1.18    0.450 0.00408   
#>  9 Antonio Sanabria    1.05    0.386 0.00877   
#> 10 Nicolò Barella      1.11    0.417 0.0138

We can visually check the consistency of the results from this method with the prior two with scatter plots of the estimated underperforming unlikeliness from each.

If two of the approaches were perfectly in agreement, then each point, representing one of the 593 evaluated players, would fall along the 45-degree slope 1 line.

With that in mind, we can see that approach 3 is slightly more precise in relation to approach 1, but approach 3 tends to assign slightly higher percentiles to players on the whole. The results from approaches 2 and 3 also have a fair degree of agreement, and the results are more uniformly calibrated.

Stepping back from the results, what can we say about the principles of the methodology?

• Parametric Nature: The reliance on a Gamma distribution for modeling a player’s performance is both a strength and a limitation. The Gamma distribution is apt for positive, skewed continuous variables, making it suitable for modeling goals-to-xG ratios. However, the dependency on a single distribution type may restrict the scope of analysis.
• Sensitivity to Distribution Parameters: The outcomes of this methodology are highly sensitive to the parameters defining each player’s Gamma distribution. Small adjustments in shape or rate parameters can significantly alter the distribution, causing substantial shifts in the percentile outcomes of player performances. This sensitivity underscores the need for careful parameter selection and calibration.
• Flexibility of the Model: Despite its sensitivity, the Gamma distribution offers considerable flexibility. It allows for fine-tuning of the model to better fit the data, which can be advantageous for capturing the nuances of different players’ shot profiles.