Predicting the Winner of Super Bowl LV

TL;DR

• Using Pythagorean expectation we should expect the Baltimore Ravens to be Super Bowl Champions
• Using a Bradley-Terry model we should expect the Kansas City Chiefs to be Super Bowl champions
• Seems like it will be a good year for the AFC

It’s Playoff Time in the NFL!. While my team has unfortunately missed the playoffs, I wanted to take advantage of the season to try to predict who will win the Super Bowl this year through two different mechanisms:

1. Pythagorean Expectation
2. Simulation using Bradley-Terry Models

Getting the Data

While ideally having more historical data would be better, I’m going to keep this exercise quick and dirty by only using the data from the 2020 NFL Regular Season which recently concluded. Data for this season can be easily imported using the nflfastR package. By using the fast_scraper_schedules function, I can quickly get all the games and their results for the 2020 season.

library(tidyverse)
library(nflfastR)
library(scales)

#Get Season 2020 Schedule and results
nfl_games <- fast_scraper_schedules(2020) %>% 
  #Weeks Beyond Week 17 Are the Playoffs
  filter(week <= 17)

knitr::kable(head(nfl_games, 3))

The package returned both the data I’m looking for, but also a lot of additional data that could be used if necessary (day of week, dome vs. outdoor, etc.).

Method 1: Pythagorean expectation

Pythagorean expectation was developed by Bill James for Baseball and estimates the % of games that a team “should win” based on runs scored and runs allowed.

It was adapted for Pro Football by Football Outsiders to use the following formula:

Football Outside Almanac in 2011 stated that “From 1988 through 2004, 11 of 16 Super Bowls were won by the team that led the NFL in Pythagorean wins, while only seven were won by the team with the most actual victories”

There needs to be a little data manipulation to get the NFL schedule data into a format to calculate the pythagorean expectation. Most notably splitting each game into two rows of data to capture information on both the home team and away teams.

p_wins <- nfl_games %>% 
  pivot_longer(
    cols = c(contains('team')),
    names_to = "category",
    values_to = 'team'
  ) %>% 
  mutate(points_for = (category=='home_team')*home_score+
           (category=='away_team')*away_score,
         points_against = (category=='away_team')*home_score+
           (category=='home_team')*away_score
  ) %>% 
  group_by(team) %>%
  summarize(pf = sum(points_for, na.rm = T),
            pa = sum(points_against, na.rm = T),
            actual_wins = sum(points_for > points_against, na.rm = T),
            .groups = 'drop'
  ) %>% 
  mutate(p_expectation = pf^2.37/(pf^2.37+pa^2.37)*16)

By pythagorean expectation the top 3 teams in the NFL are:

According to Pythagorean Expectation, the Baltimore Ravens are the best team in the NFL while the formula would say that the Kansas City Chiefs, the team with the most actual wins, “should” have only had 10.5 wins vs. the 14 actual wins they had.

library(ggrepel)
p_wins %>% 
  mutate(diff_from_exp = actual_wins - p_expectation) %>% 
  ggplot(aes(x = actual_wins, y = p_expectation, fill = diff_from_exp)) + 
    geom_label_repel(aes(label = team)) + 
    geom_abline(lty = 2) + 
    annotate("label", x = 1, y = 10, hjust = 'left', label = "Underachievers") +
    annotate("label", x = 10, y = 5, hjust = 'left', label = "Overachievers") +
    labs(x = "Actual Wins", y = "Expected Wins", 
         title = "What NFL Teams Over/Under Performed?", 
         caption = "Expected Wins Based on Pythagorian Expectation") + 
    scale_fill_gradient2(guide = F) + 
    cowplot::theme_cowplot()

Method #2: Simulation with Bradley-Terry Models

Bradley-Terry Models are probability models to predict the outcomes of paired comparisons (such as sporting events or ranking items in a competition).

In this case, to predict the future winner of Super Bowl LV. I’ll be using regular season data to estimate “ability parameters” for each team and then using those parameters to run simulations to estimate the winners of the NFL Playoff Match-ups.

The Bradley-Terry Model can be fit using the BradleyTerry2 package.

BTm(cbind(win1, win2), team1, team2, ~ team, id = "team", data = sports.data)

#Get List of All Teams
all_teams <- sort(unique(nfl_games$home_team))

nfl_shaped <- nfl_games %>%
  mutate(
    home_team = factor(home_team, levels = all_teams),
    away_team = factor(away_team, levels = all_teams),
    home_wins = if_else(home_score > away_score, 1, 0),
    away_wins = if_else(home_score < away_score, 1, 0) 
  ) %>% 
  group_by(home_team, away_team) %>% 
  summarize(home_wins = sum(home_wins),
            away_wins = sum(away_wins),
            .groups= 'drop') 

knitr::kable(head(nfl_shaped, 3), align = 'c')

library(BradleyTerry2)
base_model <- BTm(cbind(home_wins, away_wins), home_team, away_team,
                  data = nfl_shaped, id = "team")

base_abilities <- qvcalc(BTabilities(base_model)) %>% 
  .[["qvframe"]] %>% 
  as_tibble(rownames = 'team') %>% 
  janitor::clean_names()

knitr::kable(base_abilities %>% 
               mutate(across(where(is.numeric), round, 2)) %>% 
               head(3),
             align = 'c')

Step 4: Simulating Playoff Matchups

To determine each team’s likelihood of winning their match-up I run 1,000 simulations pulling from a distribution of the ability scores using team ability and standard error as parameters. The percent of those 1,000 simulations won by each each represents the likelihood of winning that match-up.

To generate the 1,000 simulations I use the tidyr::crossing() function to replicate each row 1,000 times; then using dplyr to summarize over all simulations.

Since running this for any arbitrary combination of teams isn’t too time consuming, I’ll generate every combination of playoff team across the NFC and AFC even though at least half of these comparisons will be impossible in practice.

playoff_teams = c('BAL', 'BUF', 'CHI', 'CLE', 'GB', 'IND', 'KC', 'LA', 'NO',
                  'PIT', 'SEA', 'TB', 'TEN', 'WAS')

comparisons <- base_abilities %>% 
  filter(team %in% playoff_teams)

#Generate All Potential Combination of Playoff Teams
comparisons <- comparisons %>% 
  rename_with(~paste0("t1_", .x)) %>% 
  crossing(comparisons %>% rename_with(~paste0("t2_", .x)))  %>% 
  filter(t1_team != t2_team)

#Run 1000 Simulations per comparison
set.seed(20210107)

#Draw from Ability Distribution
simulations <- comparisons %>% 
  crossing(simulation = 1:1000) %>% 
  mutate(
    t1_val = rnorm(n(), t1_estimate, t1_quasi_se),
    t2_val = rnorm(n(), t2_estimate, t2_quasi_se),
    t1_win = t1_val > t2_val,
    t2_win = t2_val > t1_val
  )

#Roll up the 1000 Results
sim_summary <- simulations %>% 
  group_by(t1_team, t2_team, t1_estimate, t2_estimate) %>% 
  summarize(t1_wins_pct = mean(t1_win), #Long-Term Average Winning % for Team 1
            t2_wins_pct = mean(t2_win), #Long-Term Average Winning % for Team 2
            .groups = 'drop') %>% 
  mutate(
    #Create a label for the winner
    winner = if_else(t1_wins_pct > t2_wins_pct, t1_team, t2_team)
  )

Step 5: And the winner is….

Now since we have all potential combinations we can step through each of the games on the schedule to determine the likelihood of winning that match-up. For rounds after the initial wild-card round, the teams are re-seeded so the #1 seed will play whatever the lowest winning seed is (can be anywhere from #4 to #7). While initially I wanted to look at each team’s likelihood of winning the Super Bowl, I couldn’t quite figure out how to easily determine the probability of each scenario given the re-seeding process. So I will just step through each round based on the result of the previous round.

For simplicity I define a function to take in the two teams and return the ability scores from the simulations above.

winners <- function(t1, t2){
  dt = sim_summary %>% filter(t1_team == t1 & t2_team == t2) %>% 
    inner_join(
      nflfastR::teams_colors_logos %>% 
        filter(team_abbr == t1) %>% 
        select(t1_team = team_abbr, t1_name = team_name),
      by = "t1_team"
    ) %>% 
    inner_join(
      nflfastR::teams_colors_logos %>% 
        filter(team_abbr == t2) %>% 
        select(t2_team = team_abbr, t2_name = team_name),
      by = "t2_team"
    )
  
  return(
     list(
       team1 = dt$t1_name,
       team1_prob = dt$t1_wins_pct,
       team2 = dt$t2_name,
       team2_prob = dt$t2_wins_pct,
       winner = if_else(dt$winner == dt$t1_team, dt$t1_name, dt$t2_name)
     )
  )
}