Ranking penalty takers

Since Brazil was eliminated from the World Cup on Sunday, there’s been some discussion on why Ancelotti picked Bruno Guimarães to shoot that first penalty. Supposedly, it was due to Bruno’s historical conversion rate. Let’s investigate this claim.

We begin by plotting penalty data for all current Brazil national team players in recent years prior to the Brazil vs Norway match (hover over each penalty to see match details):

Code

import altair as alt
import polars as pl

df = (
    pl.read_csv("data.csv")
    .filter(pl.col("date").is_between(pl.date(2023, 1, 1), pl.date(2026, 7, 4)))
    .with_columns(
        pl.col("scored").cast(pl.Boolean),
        pl.col("scored").sum().over("player").alias("player_goals"),
        pl.col("scored").count().over("player").alias("player_shots"),
        pl.row_index().over("player", order_by=["date"]).alias("shot"),
    )
    .with_columns(
        pl.when(pl.col("scored"))
        .then(pl.lit("Scored"))
        .otherwise(pl.lit("Missed"))
        .alias("result"),
        pl.format(
            "{} ({}/{})", "player", pl.col("player_goals"), pl.col("player_shots")
        ).alias("label"),
    )
)

order = df.sort(["player_goals", "player_shots"], descending=True)["label"].to_list()
(
    alt.Chart(df)
    .mark_point(filled=True, size=120)
    .encode(
        x=alt.X("shot:O", title="Shot", axis=None),
        y=alt.Y("label:N", sort=order, title=None),
        color=alt.Color(
            "result:N",
            title="Result",
            scale=alt.Scale(range=["#e0e0e0", "steelblue"]),
        ),
        tooltip=[
            alt.Tooltip("player:N", title="Player"),
            alt.Tooltip("date:N", title="Date"),
            alt.Tooltip("competition:N", title="Competition"),
            alt.Tooltip("club:N", title="Team"),
            alt.Tooltip("opponent:N", title="Opponent"),
            alt.Tooltip("result:N", title="Result"),
            alt.Tooltip("goalkeeper:N", title="Goalkeeper"),
            alt.Tooltip("scored:N", title="Scored"),
        ],
    )
    .properties(
        title=alt.TitleParams(
            text="Penalty shots by Brazilian NT players",
            subtitle=[
                "Each dot is one penalty shot between Jan 1st, 2023 and July 4th, 2026."
            ],
        ),
        width="container",
        height=600,
    )
)

How can we quantify the uncertainty around the conversion rate and pick the best shooter for a single penalty? There are many factors, such as the league and goalie of each shot. For simplicity, let’s stick to the player-level shot data. That still leaves us with the uncertainty around the player’s own rate. Intuitively, one might hesitate to pick players such as Éderson, who only shot a single penalty. One might also be wary of picking a player like Fabinho, with a low conversion rate.

A natural way to encode both intuitions is to use a hierarchical partial pooling. The idea is to model a population of penalty takers, and use that estimate to regularize each player’s individual rate. Players with few shots get pulled toward the population mean, while players with many shots are allowed to deviate from it.

Concretely, we assume each player i has a conversion rate θ_i, drawn from a Beta(ϕκ, (1 − ϕ)κ) distribution, where ϕ ∈ [0, 1] is the mean conversion rate for a top-level player and κ > 0 controls how tightly players cluster around it. We place a Beta(8, 2) prior on ϕ, so that the mean conversion rate of the population is centered on 80%. We place a non-informative Pareto(1.5) prior on κ. Finally, each player’s goal count comes from a Binomial(n_i, θ_i).

Code

import arviz as az
import jax
import jax.numpy as jnp
import numpyro
import numpyro.distributions as dist
from numpyro.infer import MCMC, NUTS

numpyro.set_platform("cpu")
numpyro.set_host_device_count(4)

summary = df.unique(["label", "player", "player_goals", "player_shots"]).sort("player")
labels = summary["label"].to_list()
players = summary["player"].to_list()
shots = jnp.array(summary["player_shots"].to_numpy())
goals = jnp.array(summary["player_goals"].to_numpy())
N = len(players)


def model():
    phi = numpyro.sample("phi", dist.Beta(8, 2))
    kappa_log = numpyro.sample("kappa_log", dist.Exponential(1.5))
    kappa = numpyro.deterministic("kappa", jnp.exp(kappa_log))

    with numpyro.plate("players", N):
        theta = numpyro.sample("theta", dist.Beta(phi * kappa, (1.0 - phi) * kappa))

    numpyro.sample("y", dist.Binomial(total_count=shots, probs=theta), obs=goals)


kernel = NUTS(model, target_accept_prob=0.95)
mcmc = MCMC(
    kernel,
    num_warmup=2000,
    num_samples=4000,
    num_chains=4,
    progress_bar=False,
)
seed = jax.random.PRNGKey(0)
mcmc.run(seed)

idata = az.from_numpyro(
    mcmc,
    coords={"players": players},
    dims={"theta": ["players"]},
)

assert float(az.rhat(idata).to_dataset().to_array().max()) < 1.001

Code

theta_samples = idata.posterior["theta"].values.reshape(-1, N)
theta_hdi = az.hdi(idata, var_names=["theta"], prob=0.9)["theta"]
hdi_lo = theta_hdi.sel(ci_bound="lower").values
hdi_hi = theta_hdi.sel(ci_bound="upper").values
means = theta_samples.mean(axis=0)

inference = pl.DataFrame(
    {
        "label": labels,
        "player": players,
        "shots": shots.tolist(),
        "goals": goals.tolist(),
        "mean": means.tolist(),
        "hdi_lo": hdi_lo.tolist(),
        "hdi_hi": hdi_hi.tolist(),
    }
)

base = alt.Chart(inference)

order = inference.sort("mean", descending=True)["label"].to_list()
bars = base.mark_rule(strokeWidth=2).encode(
    x=alt.X(
        "hdi_lo:Q",
        title="Conversion rate",
        axis=alt.Axis(format=".0%"),
        scale=alt.Scale(zero=False),
    ),
    x2="hdi_hi:Q",
    y=alt.Y("label:N", sort=order, title=None),
    color=alt.value("steelblue"),
)

points = base.mark_point(filled=True, size=120).encode(
    x=alt.X("mean:Q"),
    y=alt.Y("label:N", sort=order, title=None),
    color=alt.value("steelblue"),
    tooltip=[
        alt.Tooltip("player:N", title="Player"),
        alt.Tooltip("mean:Q", title="Posterior mean", format=".1%"),
        alt.Tooltip("hdi_lo:Q", title="90% HDI low", format=".1%"),
        alt.Tooltip("hdi_hi:Q", title="90% HDI high", format=".1%"),
    ],
)

(bars + points).properties(
    title=alt.TitleParams(
        text="Penalty conversion posteriors",
        subtitle=["Dot is the posterior mean, bar is the 90% HDI."],
    ),
    width="container",
    height=600,
)

Ranking by posterior mean, Raphinha ends up at the top, with a tight interval. Rayan follows in second place, but his wider interval reflects the weaker evidence. Players with few shots get pulled toward the population mean, shrinking extreme rates like Éderson’s.

Guimarães ends up in third place, tied with Cunha and above Vini and Neymar, despite the large uncertainty around his estimate. We could be more conservative and rank by the lower bound of the posterior, in which case Paquetá would rank third, above Guimarães and Cunha.

Rayan ranks better than Guimarães in both mean and lower bound, so he would have probably been a better pick. Raphinha would have ranked even higher, but he didn’t play in the match against Norway.