Predições de Resultados de Futebol

Vamos documentar uma abordagem de predição de partidas de futebol utilizando modelos de regressões de Poisson.

Data de Publicação

20 de dezembro de 2022

Introdução

Durante a Copa do Mundo de futebol de 2022, o interesse por predições de resultados de futebol aumentou consideravelmente. Neste sentido, este post irá demonstrar uma abordagem simples que pode gerar resultados razoáveis. O modelo utilizado neste post é baseado no trabalho The Methodology of Red, a Football Forecasting Model from theUniversity of Reading de J. James Reade.

Basicamente, o autor utiliza um modelo de Poisson para prever as probabilidades de gols marcados por cada time. A partir dessas probabilidades, é possível estimar também a chance de vitória, empate ou derrota de cada time. Como variáveis, iremos utilizar o Elo Rating de cada time e uma variável que indica se o jogo ocorreu em campo neutro ou não. O modelo a ser estimado para cada time é: $Y_{i j t}^{n} = \exp (α^{n} + β_{1}^{n} R_{i t} + β_{2}^{n} R_{j t} + β_{3}^{n} N_{i j t}) ε_{i j t}^{n}, n \in {t i m e 1, t i m e 2},$ em que $Y_{i j t}^{n}$ representa o número de gols marcados pelo time $n$ (time 1 ou time2) em uma partida entre o time $i$ (time 1) e o time $j$ (time 2) na data $t$ . As variáveis $R_{i t}$ e $R_{j t}$ representam, respectivamente, os ratings dos times $i$ e $j$ antes da partida de data $t$ . Por último, a variável $N_{i j t}$ é um variável dummy que indica se a partida foi realizada em um campo neutro. Reforçando, iremos estimar dois modelos: um para o time 1 e outro para o time 2.

Pacotes

Os pacotes utilizados estão listados abaixo.

library(EloRating)
library(tidymodels)
library(tidyverse)

# Tema para os gráficos
theme_set(theme_bw())

Dados

Serão utilizados dados de partidas entre seleções que estão disponíveis neste repositório: martj42/international_results. A primeira partida disponível na base data de 30/11/1872.

data <- read_csv(
  file = "https://github.com/martj42/international_results/raw/master/results.csv",
  show_col_types = FALSE
)

# Remove partidas sem resultados e ajusta os nomes das variáveis de home e away
# para team1 e team2
data <- data |>
  filter(!is.na(home_score)) |>
  rename(
    team1 = home_team,
    team2 = away_team,
    team1_score = home_score,
    team2_score = away_score
  )

head(data)

# A tibble: 6 × 9
  date       team1    team2    team1_score team2_score tournament city   country
  <date>     <chr>    <chr>          <dbl>       <dbl> <chr>      <chr>  <chr>  
1 1872-11-30 Scotland England            0           0 Friendly   Glasg… Scotla…
2 1873-03-08 England  Scotland           4           2 Friendly   London England
3 1874-03-07 Scotland England            2           1 Friendly   Glasg… Scotla…
4 1875-03-06 England  Scotland           2           2 Friendly   London England
5 1876-03-04 Scotland England            3           0 Friendly   Glasg… Scotla…
6 1876-03-25 Scotland Wales              4           0 Friendly   Glasg… Scotla…
# ℹ 1 more variable: neutral <lgl>

Elo Rating

O código abaixo irá calcular os scores (rating) do times utilizando o pacote EloRating. A fórmula de atualização dos scores depende um fator atualização $k$ . Quanto maior esse valor, maior será o tamanho do ajuste dos scores dos times envolvidos em cada partida. Aqui, será utilizado o valor 20 igual ao adotado por J. James Reade.

# Determina o vencendor e o perdedor
# Em caso de empate, não importa a ordem, mas é preciso indicar que é um empate
data <- data |>
  mutate(
    winner = case_when(
      team1_score >= team2_score ~ team1,
      TRUE ~ team2
    ),
    loser = case_when(
      team1_score < team2_score ~ team1,
      TRUE ~ team2
    ),
    draw = team1_score == team2_score,
    match_id = 1:n()
  )

# Calcula os ratings
elo_fit <- elo.seq(
  winner = data$winner,
  loser = data$loser,
  Date = data$date,
  draw = data$draw,
  k = 20
)

# Adiciona os ratings ao data.frame.
# Note que utilizamos o rate do dia anterior à partida para que seja refletido
# o score de cada time antes da partida.
data <- data |>
  mutate(
    team1_elo = extract_elo(
      eloobject = elo_fit,
      extractdate = pmax(min(elo_fit$truedates), date - 1),
      IDs = team1
    ),
    team2_elo = extract_elo(
      eloobject = elo_fit,
      extractdate = pmax(min(elo_fit$truedates), date - 1),
      IDs = team2
    )
  )

Treino, validação e teste

Como não vamos realizar nenhuma seleção de hiperparâmetros, já que o modelo é bastante simples, iremos dividir os dados apenas em dois conjuntos: treino e teste.

data <- data |>
  mutate(friendly = tournament == "Friendly") |>
  filter(date >= as.Date("1930-01-01"))

# Split - Treino e Teste
train_test <- initial_time_split(data, prop = 0.8)
train <- training(train_test)
test <- testing(train_test)

Resultados

Neste primeiro bloco de código, os modelos para os times 1 e 2 são estimados usando a função glm() com a opção family = poisson.

fit1 <- glm(
  formula = team1_score ~ team1_elo + team2_elo + neutral,
  data = train,
  family = "poisson"
)

summary(fit1)

pred1 <- predict(fit1, test, type = "response")

fit2 <- glm(
  formula = team2_score ~ team1_elo + team2_elo + neutral,
  data = train,
  family = "poisson"
)
summary(fit2)

pred2 <- predict(fit2, test, type = "response")


Call:
glm(formula = team1_score ~ team1_elo + team2_elo + neutral, 
    family = "poisson", data = train)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.412e+00  4.211e-02  33.528  < 2e-16 ***
team1_elo    2.005e-03  3.357e-05  59.716  < 2e-16 ***
team2_elo   -2.876e-03  3.315e-05 -86.761  < 2e-16 ***
neutralTRUE -7.586e-02  9.508e-03  -7.979 1.48e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 58873  on 35068  degrees of freedom
Residual deviance: 50166  on 35065  degrees of freedom
AIC: 119035

Number of Fisher Scoring iterations: 5


Call:
glm(formula = team2_score ~ team1_elo + team2_elo + neutral, 
    family = "poisson", data = train)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  5.695e-01  5.158e-02   11.04   <2e-16 ***
team1_elo   -2.847e-03  4.023e-05  -70.75   <2e-16 ***
team2_elo    2.324e-03  4.139e-05   56.14   <2e-16 ***
neutralTRUE  2.990e-01  1.089e-02   27.45   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 53631  on 35068  degrees of freedom
Residual deviance: 46925  on 35065  degrees of freedom
AIC: 100215

Number of Fisher Scoring iterations: 5

Neste segundo bloco, são computadas as probabilidades de vitória do time 1, empate e vitória do time 2. Para isso, são utilizadas as predições anteriores como o parâmetro ( $λ$ ) da distribuição de Poisson e são calculadas as probabilidades de cada time marcar 0 a 10 gols. Multiplicando as probabilidades são obtidas as probabilidades de cada placar possível. Por último, a probabilidade do time 1 ganhar é calculada como a soma das probabilidades dos placares nos quais o time 1 marca mais gols que o time2. As probabilidades de empate e vitória do time 2 são calculadas de forma similar.

preds <- map_df(seq_len(nrow(test)), ~ {
  lambda1 <- pred1[.x]
  lambda2 <- pred2[.x]

  prob1 <- dpois(0:10, lambda1)
  prob2 <- dpois(0:10, lambda2)
  probs <- outer(prob1, prob2)
  probs <- probs / sum(probs)

  prob_team1 <- sum(lower.tri(probs) * probs)
  prob_team2 <- sum(upper.tri(probs) * probs)
  prob_draw <- sum(diag(probs))

  data.frame(
    .pred_team1 = prob_team1,
    .pred_draw = prob_draw,
    .pred_team2 = prob_team2
  )
})

O código abaixo computa algumas métricas (log-loss e acurácia) para a base de testes.

# Junta as predições aos dados de teste, cria a variável de target e a predição
# por classe.
test <- test |>
  bind_cols(preds) |>
  mutate(
    target = case_when(
      team1_score > team2_score ~ "team1",
      team1_score == team2_score ~ "draw",
      team1_score < team2_score ~ "team2"
    ),
    target = factor(target, c("team1", "draw", "team2")),
    # predição por classe
    pred_class = case_when(
      .pred_team1 == pmax(.pred_team1, .pred_team2, .pred_draw) ~ "team1",
      .pred_draw == pmax(.pred_team1, .pred_team2, .pred_draw) ~ "draw",
      .pred_team2 == pmax(.pred_team1, .pred_team2, .pred_draw) ~ "team2"
    ),
    pred_class = factor(pred_class, c("team1", "draw", "team2"))
  )

score_acc <- accuracy(
  data = test,
  truth = target,
  estimate = pred_class
)

score_ll <- mn_log_loss(
  data = test,
  truth = target,
  contains(".pred")
)

bind_rows(score_acc, score_ll)

# A tibble: 2 × 3
  .metric     .estimator .estimate
  <chr>       <chr>          <dbl>
1 accuracy    multiclass     0.588
2 mn_log_loss multiclass     0.899

E qual é o desempenho do modelo para as partidas da Copa do Mundo de 2022?

score_acc_wc2022 <- accuracy(
  data = test |>
    filter(date >= as.Date("2022-01-01"), tournament == "FIFA World Cup"),
  truth = target,
  estimate = pred_class
)

score_ll_wc2022 <- mn_log_loss(
  data = test |>
    filter(date >= as.Date("2022-01-01"), tournament == "FIFA World Cup"),
  truth = target,
  contains(".pred")
)

bind_rows(score_acc_wc2022, score_ll_wc2022)

# A tibble: 2 × 3
  .metric     .estimator .estimate
  <chr>       <chr>          <dbl>
1 accuracy    multiclass     0.531
2 mn_log_loss multiclass     1.03

Verifica-se uma acurácia de 53,1% e log-loss de 1.03. Estes resultados indicam uma acurácia menor do que aquela obtida para o conjunto total de teste. Por fim, para dar uma sensibilidade da qualidade dos resultados obtidos por esse modelo, deixo aqui uma comparação realizada pelo Octosport:

The FIFA world cup has just ended. It is time to share our model performance over the 64 games. We reached an accuracy of 56.2% and a log-loss of -0.982 on the final-time winner prediction (90 minutes).For comparison, famous 538 did 53.1% with a log-loss of -1.031 while kickoff.ai did 54.6% and -1.028.