# Intro

This report examines the impact on a team's win probability from our proposed
rule changes to incentivize less returns on punts by NFL teams.

The NFL provides an API for accessing play-by-play data for all games dating back to the 2009 season.
We designed the `R` package [`nflscrapR`](https://github.com/maksimhorowitz/nflscrapR) for accessing and compiling data from this API.

We analyze the subset of all punt plays in 2016-17 seasons pre-, post-, regular from this API. All play-by-play going back to 
2009 is already accessible in one of our [`nflscrapR`-data](https://github.com/ryurko/nflscrapR-data) GitHub repository and is also
available via [Kaggle](https://www.kaggle.com/maxhorowitz/nflplaybyplay2009to2016#NFL%20Play%20by%20Play%202009-2018%20(v5).csv) as well
due to our collaborator Maksim Horowitz.

# Data

We first gather NFL.com play-by-play data accessed via the `nflscrapR` in `R`, 
with all pre-, post-, and regular season data for the 2016 and 2017 seasons available here: https://github.com/ryurko/nflscrapR-data/tree/master/play_by_play_data

Then we only grab punts and plays following punts, making necessary indicators denoting what type of punt outcome there was:
fair catch, downed, touchback, out of bounds, or just a return.

```{r, warning = FALSE, message = FALSE}
# Access tidyverse:
# install.packages("tidyverse")
library(tidyverse)

# Load the play-by-play data from the 2016 and 2017 season available in the repository above:

# Create the combination of year and type of games to access data for:
season_type_data <- expand.grid(c(2016, 2017), 
                              c("pre", "regular", "post")) %>%
  rename(season = Var1, type = Var2)

# Join all of these seasons together:
pbp_data <- map2_dfr(season_type_data$season, as.character(season_type_data$type),
                    .f = function(season, type) {
                      # Shorten abbreviation to reg for regular season:
                      type_abbr <- if_else(type == "regular", "reg", type)
                      
                      # Now load the data:
                      read_csv(paste0(
                        "https://raw.githubusercontent.com/ryurko/nflscrapR-data/master/play_by_play_data/", 
                        type, "_season/", type_abbr, "_pbp_", season, ".csv"),
                        progress = FALSE) %>%
                        select(-fumble_recovery_2_yards)
                    })

# Access all punts as well as any play following the punt in order to have
# proper win probability calculation:
punt_pbp_data <- pbp_data %>%
  group_by(game_id) %>%
  # Create booleans for the certain situations we're interested in making
  # rule changes for:
  mutate(follow_punt = ifelse(lag(play_type == "punt"), TRUE, FALSE),
         follow_fair_catch = ifelse(lag(punt_fair_catch) == 1, TRUE, FALSE),
         follow_punt_oob = ifelse(lag(punt_out_of_bounds) == 1, TRUE, FALSE),
         follow_punt_touchback = ifelse(follow_punt & lag(touchback == 1), 
                                        TRUE, FALSE),
         follow_downed_punt = ifelse(lag(punt_downed) == 1, TRUE, FALSE)) %>%
  filter(play_type == "punt" |
         follow_punt) %>%
  ungroup()

```

# Proposed Rule Change

Since our proposed rule change to lower the concussion rate on punts is to simply
incentivize teams to return less often, it is necessary to view the impact of our
proposed incentives on the game itself. We propose the following two rule changes:

#. Return team gains 5 additional yards on a fair catch (__incentive for return team__).
#. Return team loses 5 yards from the point the ball goes out of bounds in the air (__incentive for the punting team__). 

There are no proposed incentive structures for punts that land in play, regular 
rules will apply.

# Distribution of Return Yards

The histogram below displays the distribution of return yards for all punt 
returns (excluding fair catch, touchbacks, out of bounds, and downed punts) with 
a reference line at our proposed shift of five yards. Roughly `r round(digits = 2, nrow(filter(punt_pbp_data,play_type == "punt", punt_fair_catch == 0 & punt_out_of_bounds == 0 & touchback == 0 & punt_downed == 0, return_yards <= 5)) / nrow(filter(punt_pbp_data,play_type == "punt", punt_fair_catch == 0 & punt_out_of_bounds == 0 & touchback == 0 & punt_downed == 0)) * 100)`% of punt
returns are less than or equal to 5 yards.

```{r, fig.width = 9, fig.height = 7, fig.align = "center"}
punt_pbp_data %>%
  filter(play_type == "punt",
         punt_fair_catch == 0 & punt_out_of_bounds == 0 & 
           touchback == 0 &
           punt_downed == 0) %>%
  ggplot(aes(x = return_yards)) +
  geom_histogram(fill = "darkblue", color = "black",
                 breaks = seq(-15, 100, by = 5)) +
  scale_x_continuous(limits = c(-15, 100), breaks = seq(-10, 100, by = 10)) +
  geom_rug(color = "darkblue", alpha = 0.75) +
  geom_vline(xintercept = 5, color = "darkorange", linetype = "dashed", size = 1.5) +
  labs(x = "Return yards",
       y = "Frequency", 
       title = "Distribution of punt return yards from pre-, post-, and regular season games (2016-17)",
       subtitle = "Only includes punts with returns; dashed vertical line denoting five yard returns",
       caption = "Pelechrinis, Yurko, Ventura (2019)\nData courtesy of NFL.com, accessed with `nflscrapR`") +
  theme_bw() 
```


The empirical CDF below is displaying the same information, but now is it easier
to see that the proposed yardage shift of 5 yards is comparable to the median
of return yards which is `r punt_pbp_data %>% filter(play_type == "punt", punt_fair_catch == 0, punt_out_of_bounds == 0, touchback == 0, punt_downed == 0) %>% pull(return_yards) %>% median()` yards.

```{r, fig.width = 9, fig.height = 7, fig.align = "center"}
punt_pbp_data %>%
  filter(play_type == "punt",
         punt_fair_catch == 0 & punt_out_of_bounds == 0 & 
           touchback == 0 & punt_downed == 0) %>%
  ggplot(aes(x = return_yards)) +
  stat_ecdf(color = "darkblue") +
  scale_x_continuous(limits = c(-15, 100), breaks = seq(-10, 100, by = 10)) +
  geom_vline(xintercept = 5, color = "darkorange", linetype = "dashed") +
  labs(x = "Return yards",
       y = "Cumulative proportion of returns", 
       title = "Empirical CDF of punt return yards from pre-, post-, and regular season games (2016-17)",
       subtitle = "Only includes punts with returns; dashed vertical line denoting five yard returns",
       caption = "Pelechrinis, Yurko, Ventura (2019)\nData courtesy of NFL.com, accessed with `nflscrapR`") +
  theme_bw() 
```

# Impact on Win Probability Added

__Win probability added (WPA)__ is the change in a team's probability of winning that can be attributed to a single play.

#. Ex:  A 2-yard run on 3rd-and-1 is worth more according to WPA than a 2-yard run on 2nd-and-10.
#. Ex:  Cody Parkey's ``double-doink'' missed field goal vs. the Eagles in the 2019 Wild Card Round was worth -89\% in win probability.

We use historical data from the NFL's API to build our WP model. Our model is a generalized additive model that accounts for time 
remaining, expected score differential, down, yard line, yards to go, etc. Full details regarding both our expected points and 
win probability model are [available here on ArXiv](https://arxiv.org/abs/1802.00998).

Using this win probability model in `nflscrapR`, we calculate the 
win probability added from the proposed change for all punt returns with a
fair catch or out of bounds, and compare it to the win probability
added from returns. This means we are comparing the change in win probability
from where the returning team would have hypothetically received the ball given
the current rules to where the return ends or the 5 yard shift from the proposed
rule change. Downed punts and touchbacks are not included in the figure below.

Due to kernel restraints, we made the adjustments to the punts on a dataset available publicly in our GitHub repository here: https://github.com/ryurko/nflscrapR-data/tree/master/R/punt_competition
The code to make the below file `punt_adj_data` is located in this `Rmarkdown` file that is a repeat of this kernel but with the `nflscrapR` functionality
to calculate win probability: https://github.com/ryurko/nflscrapR-data/blob/master/R/punt_competition/rule_change_wpa.Rmd

```{r, warning = FALSE, message = FALSE, fig.width=10, fig.height=10, fig.align="center"}

# Read in the created dataset available in the nflscrapR-data repository with our WPA adjustments
punt_adj_data <- read_csv("https://raw.githubusercontent.com/ryurko/nflscrapR-data/master/R/punt_competition/punt_adj_data.csv")


# Visualize the distributions of WPA due from returns compared to the proposed
# rule changes using histograms:
punt_adj_data %>%
  # Only look at the plays following punts, removing QB kneels and overtime,
  # as well as touchbacks:
  filter(follow_punt, play_type != "qb_kneel", qtr != 5) %>%
  # Reorder and relabel the punt_type:
  mutate(punt_type = fct_relevel(punt_type, "return", "fair_catch", 
                                 "out_of_bounds"),
         punt_type = fct_recode(punt_type, `Historical punt returns` = "return", 
                                              `Rule change #1: 5 yard fair catch adjustment` = "fair_catch",
                                              #`Rule change: touchback` = "touchback",
                                              `Rule change #2: -5 yards out of bounds adjustment` = "out_of_bounds")) %>%
  # Create beeswarm plots with boxplots overtop:
  ggplot(aes(x = rule_change_wpa)) +
  geom_histogram(aes(fill = punt_type), color = "black",
                 breaks = seq(-.05, .35, by = .005)) +
  ggthemes::scale_fill_colorblind(guide = FALSE) +
  facet_wrap(~punt_type, ncol = 1) +
  # geom_density_ridges(aes(y = punt_type), rel_min_height = 0.01, jittered_points = TRUE, color = "white",
  #                     position = position_points_jitter(width = 0, height = 0),
  #                     point_shape = '|', point_size = 1, point_alpha = 0.7, alpha = 0.8,
  #                     fill = "darkblue", point_color = "darkblue") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "darkred", size = 1.5) +
  geom_text(data = data.frame("punt_type" = c("Rule change #1: 5 yard fair catch adjustment",
                                              "Rule change #2: -5 yards out of bounds adjustment"),
                              "text_label" = c("Rule change #1\nfavors return team",
                                               "Rule change #2\nfavors punting team"),
                              "x" = c(.06, .06),
                              "y" = c(450, 450)), aes(x, y, label = text_label),
            color = "darkred", size = 6) +
  scale_x_continuous(limits = c(-0.05, .2), breaks = seq(-.05, .2, by = .05)) +
  labs(title = "Distribution of win probability added from returns compared to proposed rule changes",
       #subtitle = "Calculated with respect to possession team",
       subtitle = "All punt returns from pre-, post-, and regular season games in 2016 and 2017",
       caption = "Pelechrinis, Yurko, Ventura (2019)\nData courtesy of NFL.com, accessed with `nflscrapR`",
       x = "Win probability added (WPA) with respect to return team",
       y = "Frequency") +
  #scale_y_discrete(expand = c(0.01, 0.01)) +
  #scale_x_continuous(limits = c(-.5, .5)) +
  theme_bw() +
  theme(axis.title = element_text(size = 14),
        axis.text = element_text(size = 14),
        plot.title = element_text(size = 16),
        plot.subtitle = element_text(size = 14),
        plot.caption = element_text(size = 12),
        strip.background = element_blank(),
        strip.text = element_text(size = 14))  


```

The figure below just displays the empirical CDF form of the above chart:


```{r, warning = FALSE, message = FALSE,  fig.align="center", fig.width=10, fig.height=10}
ecdf_war_plot <- punt_adj_data %>%
  # Only look at the plays following punts, removing QB kneels and overtime:
  filter(follow_punt, play_type != "qb_kneel", qtr != 5) %>%
  # Reorder and relabel the punt_type:
  mutate(punt_type = fct_relevel(punt_type, "return", "fair_catch",
                                 "out_of_bounds"),
         punt_type = fct_recode(punt_type, `Historical punt returns` = "return", 
                                              `Rule change #1: 5 yard fair catch adjustment` = "fair_catch",
                                              #`Rule change: touchback` = "touchback",
                                              `Rule change #2: -5 yards out of bounds adjustment` = "out_of_bounds")) %>%
  # Create beeswarm plots with boxplots overtop:
  ggplot(aes(x = rule_change_wpa, color = punt_type)) +
  stat_ecdf(alpha = 0.8) +
  ggthemes::scale_color_colorblind() +
  geom_vline(xintercept = 0, linetype = "dashed", color = "darkred", size = 1.5) +
  labs(title = "Empirical CDF of win probability added from returns compared to proposed rule changes",
       #subtitle = "Calculated with respect to possession team",
       subtitle = "All punt returns from pre-, post-, and regular season games in 2016 and 2017",
       caption = "Data courtesy of NFL.com, WPA generated w/ `nflscrapR`",
       x = "Win probability added (WPA) with respect to return team",
       y = "Cumulative proportion of punt return type", 
       color = "Punt return type") +
  #scale_y_discrete(expand = c(0.01, 0.01)) +
  #scale_x_continuous(limits = c(-.5, .5)) +
  theme_bw() +
  theme(axis.title = element_text(size = 10),
        axis.text = element_text(size = 10),
        plot.title = element_text(size = 12),
        plot.subtitle = element_text(size = 10),
        plot.caption = element_text(size = 8),
        strip.background = element_blank(),
        strip.text = element_text(size = 14),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10))  
library(cowplot)
ecdf_legend <- get_legend(ecdf_war_plot)

plot_grid(ecdf_war_plot + theme(legend.position = "none"),
          ecdf_legend, ncol = 1, rel_heights = c(3,1))

```

Our proposed rule changes incentive both teams to behave in ways that will __reduce the number of returned punts__.

__Fewer punt returns leads to fewer concussions__.

