Everything you’ve (n)ever wanted to know about penalty kicks: statistical models

After the descriptives of part 1, we turn to the statistical models. I came up with six questions that were of interest to me:

1. Where in the goal does a penalty have the best chance of going in?¹
2. Which takers are best at sending the goalkeeper the wrong way, and which keepers read them best?
3. Who are the best and worst penalty takers and penalty-stopping goalkeepers overall?
4. How good could the very best penalty takers be?
5. How good are teams at selecting their penalty takers?
6. How many penalties does a player need to take before their own placement history tells us more than the global average?

NoteAt a glance

Sending keeper wrong way

• The keepers dive the wrong way in 56% of penalties.
• Max Kruse is the best at sending the goalkeeper the wrong way, sending the keeper the wrong way in 25 of his 28 penalties (89%). The model estimates that he manages this 12.7pp more often versus the average goalkeeper than the average player versus the average goalkeeper.
• Ante Budimir and Neymar are best at this from players selected for the World Cup 2026.
• Darren Bennet is the worst (4/19 scored, for -10.6pp)
• Youri Tielemans is the worst among players active at WC 2026, sending the goalkeeper the wrong way only on 5 (29%) of his 17 attempts (+5.2pp).

Penalty taking skill

• Harry Kane is the best penalty taker (97 scored out of 107 (91%)) by considerable margin (1.3pp more added than the second best). The model estimates that he adds 8.2pp above the average penalty taker against the average goalkeeper.
• Chicarito is the worst penalty taker (8 scored out of 17 (47%)). Because of regularisation the model believes his ‘true’ conversion rate is -6.4pp than the average penalty taker.
• Uncertainty remains large: for only 3 players the 95% credible interval excluding the average.
• Uncertainty is especially large for the worst penalty takers: teams do not wait for complete certainty as switching your taker is so cheap.

Penalty stopping skill

• Tim Melia is the best goalkeeper (26/53 not scored, 49%). Also by a considerable margin (1.4pp more added). The model estimates that the average penalty taker scores 9.4pp less of their penalties against him.
• Among the active goalkeepers, Dayne St Clair is the best (16/35 saved) for 5.9pp less scoring chance versus average taker. He is in goal for Canada this World Cup.
• Roman Bürki is the worst penalty stopper, conceding 44 from 47 penalty kicks (94%). Part of this was bad luck as the model estimates that the average penalty taker scores 4.6pp more from their penalties versus him.

How good could the very best be?

• The estimated spread in taker skill implies that takers with a true conversion rate of ≈89% (+12pp above average) exist: roughly 1 in 1,000, so about 6 of the 5,934 takers in the data. The model just cannot say who they are as a taker’s own record only outweighs the population average after ~66 penalties, a bar only 9 takers clear.
• A true 95% taker would be about one in 538 million takers, although this far out in the tail the uncertainty is enormous (one in 1.4 million to one in 10¹⁴) and such a large number may point at model misspecification (e.g. fat-tailed model may be more appropiate).
• The best conceivable goalkeepers matter more than the best conceivable takers: a 99.9th-percentile keeper prevents 14.9pp, more than the +12.2pp the equivalent taker adds.

Placement preference

• A taker’s own placement history only outpredicts the global average (informed by foot) after about 30 penalties.
• Players seem to differ from one another mainly in how often they go down the middle, not in how often they use their weak side.

Team selection quality

• Teams are good at picking penalty takers: the correlation between a player’s likelihood of being selected and their scoring skill is estimated at 0.71 using data from the English top level and cup games.

A simple expected penalty goals model based on shot placement and kicker’s foot

We model whether an on-target penalty is scored from the kicker’s foot and where the ball crosses the goal line:

\[ \begin{aligned} \text{Goal}_i &\sim \operatorname{Bernoulli}(p_i) \\ \operatorname{logit}(p_i) &= \beta_0 + \beta_1\,\mathbb{1}[\text{foot}_i = \text{Right}] + f_{\text{foot}_i}(x_i, y_i), \end{aligned} \]

where \((x_i, y_i)\) is where the ball crosses the goal line and \(f_{\text{Left}}\), \(f_{\text{Right}}\) are thin-plate spline surfaces fit separately per foot (the by = kick_foot term). The index \(i\) runs over on-target shots only — scored or saved — since shots that miss or strike the frame are deterministically no-goal and carry no placement information here.

library(tidyverse)
source("functions_features.R")
source("functions_plotting.R")

df_male <- nanoparquet::read_parquet(
  "data/penalties_ws.parquet"
) |>
  convert_opta_to_meters() |>
  add_features() |>
  filter(!is_female_league)

# fit model
psxg_interaction <- brms::brm(
  formula = is_goal ~ kick_foot +
    s(shot_x_meters, shot_y_meters, by = kick_foot),
  # only shots on goal: post and wide are deterministically no-goal here
  data = df_male |> filter(outcome %in% c("Goal", "Saved")),
  family = bernoulli(link = "logit"),
  backend = 'cmdstanr',
  cores = 4,
  chains = 4,
  file = "models/brms_psxg_interaction"
)

# prepare marginal effects calculation
x_range <- c(post_left, post_right)
y_range <- c(0, crossbar_height)

# n_y is a resolution parameter. the heatmap is a raster, so we sample the fitted
# surface on a grid of cells and colour each one. we could sample every cm of the
# goal (~730x240 cells), but there are two reasons not to. first, it buys nothing:
# s(x, y) is a *smooth* spline with no cm-scale detail to resolve, so a fine raster
# looks identical to a coarse one (and heatmaps read as smooth long before per-pixel
# resolution). second, it is expensive in a way the cell count hides: predictions()
# evaluates each cell on *every posterior draw* and then summarises -- that is what
# produces the credible interval -- so a cm grid with ~4000 draws is hundreds of
# millions of evaluations to hold in memory. 80 rows is smooth on screen while
# keeping that posterior summary cheap; bump n_y if you ever want it crisper.
#
# keep cells roughly square: the goal is ~3x wider than tall
n_y <- 80
n_x <- round(n_y * diff(x_range) / diff(y_range))

# evaluate at cell *centres* so the raster fills the goal frame exactly and
# doesn't bleed half a cell past the posts/line into the grass
cell_centres <- function(rng, n) {
  edges <- seq(rng[1], rng[2], length.out = n + 1)
  utils::head(edges, -1) + diff(edges) / 2
}

# calculate marginal effects
psxg_grid <- marginaleffects::predictions(
  psxg_interaction,
  newdata = marginaleffects::datagrid(
    shot_x_meters = cell_centres(x_range, n_x),
    shot_y_meters = cell_centres(y_range, n_y),
    kick_foot = c("Left", "Right")
  )
) |>
  tibble::as_tibble() |>
  dplyr::select(
    kick_foot,
    shot_x_meters,
    shot_y_meters,
    prob = estimate,
    prob_lower = conf.low,
    prob_upper = conf.high
  )

To turn the fitted model into a picture, we lay a grid of points across the goal and read off the predicted scoring probability at each one, separately for a left- and right-footed taker. Every cell is one marginaleffects prediction (the posterior mean and its 95% credibility interval); the heatmap below simply colours them.

library(plotly)

# reversed so red marks low goal probability (where you don't want to shoot)
pal_psxg <- rev(wesanderson::wes_palette(
  "Zissou1Continuous",
  100,
  type = "continuous"
))

x_vals <- sort(unique(psxg_grid$shot_x_meters))
y_vals <- sort(unique(psxg_grid$shot_y_meters))

# reshape one foot's predictions into a z-matrix (rows = y, cols = x) plus a
# matching matrix of hover labels carrying the credible interval
foot_layer <- function(foot) {
  d <- psxg_grid |>
    dplyr::filter(kick_foot == foot) |>
    dplyr::arrange(shot_y_meters, shot_x_meters)
  to_mat <- function(v) {
    matrix(v, nrow = length(y_vals), ncol = length(x_vals), byrow = TRUE)
  }
  list(
    z = to_mat(d$prob),
    text = to_mat(sprintf(
      paste0(
        "Goal probability: <b>%.1f%%</b><br>",
        "95%% credible interval: %.1f%%–%.1f%%<br>",
        "x: %.2f m   y: %.2f m"
      ),
      d$prob * 100,
      d$prob_lower * 100,
      d$prob_upper * 100,
      d$shot_x_meters,
      d$shot_y_meters
    ))
  )
}

left <- foot_layer("Left")
right <- foot_layer("Right")
z_range <- range(psxg_grid$prob)

# colorbar ticks: pretty interior values, but always pin a tick at the data min
# and max so *both* ends of the scale are labelled (plotly's auto ticks label
# the top but not the minimum)
cb_tickvals <- {
  interior <- pretty(z_range)
  interior <- interior[interior > z_range[1] & interior < z_range[2]]
  c(z_range[1], interior, z_range[2])
}

# reversed Zissou1 as a plotly colorscale (red = low probability)
colorscale <- Map(
  function(i, col) list((i - 1) / (length(pal_psxg) - 1), col),
  seq_along(pal_psxg),
  pal_psxg
)

# goal frame, goal line and dive zones as static shapes
r <- diameter_post / 2
post_rect <- function(x_centre) {
  list(
    type = "rect",
    layer = "above",
    x0 = x_centre - r,
    x1 = x_centre + r,
    y0 = 0,
    y1 = crossbar_height + r,
    fillcolor = "black",
    line = list(width = 0)
  )
}
dive_line <- function(x) {
  list(
    type = "line",
    layer = "above",
    x0 = x,
    x1 = x,
    y0 = 0,
    y1 = crossbar_height,
    line = list(color = "gray50", width = 1, dash = "dash")
  )
}
goal_shapes <- list(
  post_rect(post_left),
  post_rect(post_right),
  list(
    type = "rect",
    layer = "above",
    x0 = post_left - r,
    x1 = post_right + r,
    y0 = crossbar_height - r,
    y1 = crossbar_height + r,
    fillcolor = "black",
    line = list(width = 0)
  ),
  list(
    type = "line",
    layer = "above",
    x0 = -(post_offset + 0.3),
    x1 = post_offset + 0.3,
    y0 = 0,
    y1 = 0,
    line = list(color = "darkgreen", width = 2)
  ),
  dive_line(dive_zone_left_offset),
  dive_line(dive_zone_right_offset)
)

heatmap_trace <- function(p, layer, visible) {
  add_trace(
    p,
    x = x_vals,
    y = y_vals,
    z = layer$z,
    text = layer$text,
    type = "heatmap",
    hoverinfo = "text",
    visible = visible,
    colorscale = colorscale,
    zmin = z_range[1],
    zmax = z_range[2],
    colorbar = list(
      title = list(text = "Goal probability", side = "top"),
      tickformat = ".0%",
      tickmode = "array",
      tickvals = cb_tickvals,
      orientation = "h",
      x = 0.5,
      xanchor = "center",
      y = 1.05,
      yanchor = "bottom",
      thickness = 10,
      len = 0.3,
      lenmode = "fraction"
    )
  )
}

# a ball that follows the cursor
# so the hover doubles as a to-scale sense of how big the ball is in the goal. it
# is appended as the next shape after the goal frame and moved on hover via JS.
ball_shape_index <- length(goal_shapes)
ball_js <- paste0(
  "function(el, x) {\n",
  "  var r = ",
  ball_radius,
  ";\n",
  "  var idx = ",
  ball_shape_index,
  ";\n",
  "  var ball = {type: 'circle', xref: 'x', yref: 'y',\n",
  "    x0: 0, x1: 0, y0: 0, y1: 0, visible: false, layer: 'above',\n",
  "    fillcolor: 'rgba(255,255,255,0.92)', line: {color: '#222222', width: 1.5}};\n",
  "  var init = {}; init['shapes[' + idx + ']'] = ball;\n",
  "  Plotly.relayout(el, init);\n",
  "  el.on('plotly_hover', function(d) {\n",
  "    var pt = d.points[0], u = {};\n",
  "    u['shapes[' + idx + '].x0'] = pt.x - r;\n",
  "    u['shapes[' + idx + '].x1'] = pt.x + r;\n",
  "    u['shapes[' + idx + '].y0'] = pt.y - r;\n",
  "    u['shapes[' + idx + '].y1'] = pt.y + r;\n",
  "    u['shapes[' + idx + '].visible'] = true;\n",
  "    Plotly.relayout(el, u);\n",
  "  });\n",
  "  el.on('plotly_unhover', function(d) {\n",
  "    var u = {}; u['shapes[' + idx + '].visible'] = false;\n",
  "    Plotly.relayout(el, u);\n",
  "  });\n",
  "}"
)

# height is pinned in pixels: with the 1:1 scaleanchor lock, plotly otherwise
# stretches the y-range to fill a tall container, leaving huge empty bands above
# and below the (wide, short) goal. width stays responsive (out-width 100%).
plot_ly(height = 400) |>
  heatmap_trace(left, visible = TRUE) |>
  heatmap_trace(right, visible = FALSE) |>
  layout(
    title = list(
      text = "Post-shot goal probability across the goal mouth",
      y = 0.97,
      yref = "container",
      yanchor = "top",
      x = 0.5
    ),
    hovermode = "closest",
    # symmetric l/r so the (centred) data area lines up under the centred title
    margin = list(t = 95, b = 35, l = 45, r = 45),
    xaxis = list(
      title = "",
      range = c(-(post_offset + 0.05), post_offset + 0.05),
      zeroline = FALSE,
      showgrid = FALSE,
      constrain = "domain"
    ),
    yaxis = list(
      title = "",
      range = c(-0.25, crossbar_height + 0.05),
      zeroline = FALSE,
      showgrid = FALSE,
      scaleanchor = "x",
      scaleratio = 1
    ),
    shapes = goal_shapes,
    updatemenus = list(list(
      type = "buttons",
      direction = "right",
      x = 0.02,
      xanchor = "left",
      y = 1.05,
      yanchor = "bottom",
      buttons = list(
        list(
          method = "restyle",
          label = "Left-footed",
          args = list("visible", list(TRUE, FALSE))
        ),
        list(
          method = "restyle",
          label = "Right-footed",
          args = list("visible", list(FALSE, TRUE))
        )
      )
    ))
  ) |>
  htmlwidgets::onRender(ball_js)