---
title: "2022 NFL Big Data Bowl Submission"
subtitle: "Regularized Adjusted Field Position Contribution"
author: "Matt Johnson, UNC Chapel Hill"
date: "January 6, 2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(reactable)
library(htmltools)
library(knitr)
library(Matrix)
library(kableExtra)
```

```{r, include = F}

Sparse.Df <- readMM("../input/bdb-final-data2/Sparse.Matrix.txt")

punts <- read_csv("../input/bdb-final-data2/clean.plays.csv")

player.index <- read_csv("../input/bdb-final-data2/player.index.csv")

#first lets attach the unique newIds to the Sparse.Df as it's own column

#we need to convert the sparse matrix into a large data frame

Sparse.tib <- as.data.frame(as.matrix(Sparse.Df))

#now we can figure out which column corresponds to which player with the player.index

col.names <- player.index %>%
    arrange(ColIdx)  %>%
    distinct(nflId, ColIdx) %>%
    select(nflId)

names(Sparse.tib) <- as.character(unlist(col.names))

#now we just put in all the variables we need

Sparse.tib.named <- Sparse.tib %>% 
  mutate(newId = unique(player.index$newId)) %>% 
  mutate(NYG = punts$NetYardsGained) %>% 
  mutate(Field.Pos = punts$field.pos) %>% 
  mutate(Pen.Yrds = punts$penalty.yards.clean) %>% 
  select(newId,Field.Pos, everything())
```

# Introduction

This paper introduces a new metric for assessing player contribution for both offensive and defensive players during a punt called Regularized Adjusted Field Position Contribution (RAFPC). I adapt a metric historically used in basketball called Adjusted Plus-Minus (APM). This metric adjusts the conventional Plus-Minus metric to control for every player on the court during a stint where the same 10 players are on the court. In the NFL punt return case, this is simply the 22 players on the field during a punt. In classic APM for basketball, the variable of interest is the points scored during the stint, this will not work as well for football since points are rarely scored during punts. Instead, I build a field position metric to represent success during a punt return. Field position is the difference between where the team punted from and where the ball ended up in the resulting return. I think this metric encapsulates the success of a punt in terms of the final field position. This new metric can be used to rank players regardless of their position on how much they positively contribute towards the success of their team in punting situations. The definition of contribution is discussed later in the methodology section, but for now it is a positive number indicating how much that individual helped their team in terms of the resulting field position.

# Data
The final dataframe used for analysis had 5697 rows and 2007 columns. Each row is a play from 2018-2020 where one team successfully punted the ball away to the other team. Each column represents a player that could have been on the field during any punt for which there is tracking data. Each cell corresponding to a player is a 1 if the player's team is receiving the punt, a -1 if the player's team is punting, and a 0 if the player is not on the field. The reason I code the punting team as -1 is so they will still have positive coefficients in the regression analysis. This dataframe is inherently rank deficient and therefore I cannot use standard linear regression, rather I must use a penalized approach so the sample covariance matrix is invertible. The dataframe is rank deficient because if I removed a players column, I would still know whether he was on the field on offense or defense simply using the rest of the columns. According to this, the matrix does not have full rank.

Below is a glimpse of the dataframe. 

```{r, echo = F}
kable(Sparse.tib.named[1:10, 1:10]) %>% 
  kable_styling(full_width = T)
```

#### Exploratory Data Analysis

I complete an exploratory data analysis to help visualize the dataframe used for the modeling section. Below is a histogram and a boxplot for the response variable: Field Position. Field Position seems to be approximately normally distributed with mean -40.

```{r, echo = F, message = F}
Sparse.tib.named %>% 
  rename("Field Position" = Field.Pos) %>% 
  pivot_longer(cols = c(`Field Position`), 
               names_to = "Variable", values_to = "Value") %>% 
  select(Variable, Value) %>% 
  ggplot(aes(x = Value)) +
  geom_histogram(fill='#A4A4A4', color="purple") +
  theme_bw() +
  labs(y = "Count", 
       title = "Field Position Histogram")
```

```{r, echo = F}
Sparse.tib.named %>% 
  rename("Field Position" = Field.Pos) %>% 
  pivot_longer(cols = c(`Field Position`), 
               names_to = "Variable", values_to = "Value") %>% 
  select(Variable, Value) %>% 
  ggplot(aes(x = Value, y = Variable)) +
  geom_boxplot(fill='#A4A4A4', color="purple") +
  geom_jitter(color = "purple", alpha = .1) +
  theme_bw() +
  coord_flip() +
  labs(y = element_blank(), 
       title = "Field Position Boxplot")
```

# Methodology

#### Variable Construction

To build the Field Position response variable I first normalize the field so in every punting situation the endzone behind the punter was the 0 yardline. From left to right the yardline increases from 0 to 100. I then take the difference between where the ball was kicked and where the ball was at the end of the play. Below I use the gg_field function developed by Marschall Furman to visualize this variable. The function documentation can be found [HERE](https://github.com/mlfurman3/gg_field). 

```{r, include = F}
## gg_field function - set up as a list of annotations
gg_field <- function(yardmin=0, yardmax=120, buffer=5, direction="horiz",
                     field_color="forestgreen",line_color="white",
                     sideline_color=field_color, endzone_color="darkgreen"){
  
  ## field dimensions (units=yards)
  xmin <- 0
  xmax <- 120
  
  ymin <- 0
  ymax <- 53.33
  
  
  ## distance from sideline to hash marks in middle (70 feet, 9 inches)
  hash_dist <- (70*12+9)/36
  
  ## yard lines locations (every 5 yards) 
  yd_lines <- seq(15,105,by=5)
  
  ## hash mark locations (left 1 yard line to right 1 yard line)
  yd_hash <- 11:109
  
  ## field number size
  num_size <- 5
  
  ## rotate field numbers with field direction
  ## first element is for right-side up numbers, second for upside-down
  angle_vec <- switch(direction, "horiz" = c(0, 180), "vert" = c(270, 90))
  num_adj <- switch(direction, "horiz" = c(-1, 1), "vert" = c(1, -1))
  
  ## list of annotated geoms
  p <- list(
    
    ## add field background 
    annotate("rect", xmin=xmin, xmax=xmax, ymin=ymin-buffer, ymax=ymax+buffer, 
             fill=field_color),
    
    ## add end zones
    annotate("rect", xmin=xmin, xmax=xmin+10, ymin=ymin, ymax=ymax, fill=endzone_color),
    annotate("rect", xmin=xmax-10, xmax=xmax, ymin=ymin, ymax=ymax, fill=endzone_color),
    
    ## add yardlines every 5 yards
    annotate("segment", x=yd_lines, y=ymin, xend=yd_lines, yend=ymax,
             col=line_color),
    
    ## add thicker lines for endzones, midfield, and sidelines
    annotate("segment",x=c(0,10,60,110,120), y=ymin, xend=c(0,10,60,110,120), yend=ymax,
             lwd=1.3, col=line_color),
    annotate("segment",x=0, y=c(ymin, ymax), xend=120, yend=c(ymin, ymax),
             lwd=1.3, col=line_color) ,
    
    ## add field numbers (every 10 yards)
    ## field numbers are split up into digits and zeros to avoid being covered by yard lines
    ## numbers are added separately to allow for flexible ggplot stuff like facetting
    
    ## 0
    annotate("text",x=seq(20,100,by=10) + num_adj[2], y=ymin+12, label=0, angle=angle_vec[1],
             col=line_color, size=num_size),
    
    ## 1
    annotate("text",label=1,x=c(20,100) + num_adj[1], y=ymin+12, angle=angle_vec[1],
             colour=line_color, size=num_size),
    ## 2
    annotate("text",label=2,x=c(30,90) + num_adj[1], y=ymin+12, angle=angle_vec[1],
             colour=line_color, size=num_size),
    ## 3
    annotate("text",label=3,x=c(40,80) + num_adj[1], y=ymin+12, angle=angle_vec[1],
             colour=line_color, size=num_size),
    ## 4
    annotate("text",label=4,x=c(50,70) + num_adj[1], y=ymin+12, angle=angle_vec[1],
             colour=line_color, size=num_size),
    ## 5
    annotate("text",label=5,x=60 + num_adj[1], y=ymin+12, angle=angle_vec[1],
             colour=line_color, size=num_size),
    
    
    ## upside-down numbers for top of field
    
    ## 0
    annotate("text",x=seq(20,100,by=10) + num_adj[1], y=ymax-12, angle=angle_vec[2],
             label=0, col=line_color, size=num_size),
    ## 1
    annotate("text",label=1,x=c(20,100) + num_adj[2], y=ymax-12, angle=angle_vec[2],
             colour=line_color, size=num_size),
    ## 2
    annotate("text",label=2,x=c(30,90) + num_adj[2], y=ymax-12, angle=angle_vec[2],
             colour=line_color, size=num_size),
    ## 3
    annotate("text",label=3,x=c(40,80) + num_adj[2], y=ymax-12, angle=angle_vec[2],
             colour=line_color, size=num_size),
    ## 4
    annotate("text",label=4,x=c(50,70) + num_adj[2], y=ymax-12, angle=angle_vec[2],
             colour=line_color, size=num_size),
    ## 5
    annotate("text",label=5,x=60 + num_adj[2], y=ymax-12, angle=angle_vec[2],
             colour=line_color, size=num_size),
    
    
    ## add hash marks - middle of field
    annotate("segment", x=yd_hash, y=hash_dist - 0.5, xend=yd_hash, yend=hash_dist + 0.5,
             color=line_color),
    annotate("segment", x=yd_hash, y=ymax - hash_dist - 0.5, 
             xend=yd_hash, yend=ymax - hash_dist + 0.5,color=line_color),
    
    ## add hash marks - sidelines
    annotate("segment", x=yd_hash, y=ymax, xend=yd_hash, yend=ymax-1, color=line_color),
    annotate("segment", x=yd_hash, y=ymin, xend=yd_hash, yend=ymin+1, color=line_color),
    
    ## add conversion lines at 2-yard line
    annotate("segment",x=12, y=(ymax-1)/2, xend=12, yend=(ymax+1)/2, color=line_color),
    annotate("segment",x=108, y=(ymax-1)/2, xend=108, yend=(ymax+1)/2, color=line_color),
    
    ## cover up lines outside of field with sideline_color
    annotate("rect", xmin=0, xmax=xmax, ymin=ymax, ymax=ymax+buffer, fill=sideline_color),
    annotate("rect",xmin=0, xmax=xmax, ymin=ymin-buffer, ymax=ymin, fill=sideline_color),
    
    ## remove axis labels and tick marks
    labs(x="", y=""),
    theme(axis.text.x = element_blank(),axis.text.y = element_blank(),
          axis.ticks = element_blank()),
    
    ## clip axes to view of field
    if(direction=="horiz"){
      coord_cartesian(xlim=c(yardmin, yardmax), ylim = c(ymin-buffer,ymax+buffer), 
                      expand = FALSE)
      
    } else if (direction=="vert"){
      ## flip entire plot to vertical orientation
      coord_flip(xlim=c(yardmin, yardmax), ylim = c(ymin-buffer,ymax+buffer), expand = FALSE)
      
    }
  )
  
  return(p)

}
```

```{r, echo = F}
ggplot() +
  gg_field() +
  geom_segment(aes(x = 20, y = 27, xend = 80, yend = 27),
                  arrow = arrow(length = unit(0.5, "cm")), size = 2, color = "black") +
  geom_label(aes(label = "10", x = 20, y = 30)) +
  geom_label(aes(label = "70", x = 80, y = 30)) +
  geom_label(aes(label = "Field Position = -60", x = 35, y = 18))
```

#### Modeling

I employ a similar methodology as in Dan Rosenbaum's paper about adjusted plus minus in basketball found [HERE](http://www.82games.com/comm30.htm). In order to quantify a player's contribution while controlling for everyone on his team as well as everyone playing against him I use a penalized Ridge regression and use the beta coefficients as the value for "Contribution". I use a penalized regression to combat the rank deficiency of the dataframe described in the **Data** section. Specifically, I want to use a Ridge regression because I do not want to shrink anyone player's coefficient to 0 as I would not be able to rank their contribution. By running this type of regression I am able to evaluate players regardless of their position and assess how much they contribute during the play. Ridge regression will also remedy most multicollinearity issues. I expect there to be some since many special teams players will be on the field in pairs, e.g., Long Snapper and Punter.

*A Note on Ridge Regression:* Ridge regression is a form of multiple regression that adds a penalty factor for the $\beta$ estimates. Ridge regression allows for the sample covariance matrix to be non-invertible because of the addition of the penalty parameter. Ridge also shrinks coefficients and controls the trade-off between fit and magnitude of $\beta{}s$. In Ordinary Least Squares Regression, 
$$\hat\beta = (X^tX)^{-1}X^ty$$ 

while in Ridge, 
$$\hat\beta_{ridge} = (X^tX + \lambda I_p)^{-1}X^ty$$ 

where $I_p$ is the pxp Identity matrix and $\lambda$ is the penalty factor. I use cross validation to find the optimal $\lambda$ that minimizes the error. Ridge regression can also be written in the form of an optimization program: 
$$\min_{\beta_0, \beta} [\sum_{i=1}^N(y_i-\beta_0-x_i^t\beta)^2] \ \ \ \ s.t.\sum_{j=1}^p|\beta_j|_2 \le t$$

Where $|\beta_j|_2$ is the L2 norm of $\beta_j$.

This $\hat\beta_{ridge}$ defines the Contribution variable. 

# Results

Below is the final table showing the top 20 players with regard to positive contribution towards field position (min 25 snaps). I adjust the table to account for snap-count as there are a few players only on the field during a single return that happened to result in a big play. The percentile metric shows where that player stands with respect to every player regardless of their snap--counts. It is interesting to see that the model finds players on both sides of the ball that usually would not get any credit during the return. It also finds some specialty return players that have a direct impact on the ball moving down-field. 

```{r, include = F}
#bring in data
DF <- read_csv("../input/bdb-final-data2/top20.csv") %>%
mutate(Contribution = round(Contribution, digits = 4))

#they make us build this reactable in the RMD
#lets try to use a reactable

#first build a function to make a nice barchart, we want there to be a nice barchart in this table

bar_chart <- function(label, width = "100%", height = "16px", fill = "#00bfc4", background = NULL) {
  bar <- div(style = list(background = fill, width = width, height = height))
  chart <- div(style = list(flexGrow = 1, marginLeft = "8px", background = background), bar)
  div(style = list(display = "flex", alignItems = "center"), label, chart)
}

#ok this is a great graph, lets save it

reactable1 <- DF %>% 
  select(nflId, Name, Position, Contribution, Snaps, Percentile) %>% 
  reactable(
    columns = list(
      Contribution = colDef(name = "Contribution", align = "left", cell = function(value) {
        width <- paste0(value / max(DF$Contribution) * 100, "%")
        bar_chart(value, width = width, fill = "#7c5295", background = "#e1e1e1")
      }), 
      Snaps = colDef(name = "Snaps", align = "left", cell = function(value) {
        width <- paste0(value / max(DF$Snaps) * 100, "%")
        bar_chart(value, width = width, fill = "#56a0d3", background = "#e1e1e1")
      }), 
      Percentile = colDef(align = "center"), 
      Position = colDef(align = "center")
    )
  )
```

```{r, echo = F}
reactable1
```

# Conclusions

#### Results

I have successfully developed a new metric to assess any player's contribution during punting plays: RAFPC. This type of metric is useful for on-field decision makers as it can be used to decide which players should be on the field during a punt and it makes it easy to compare across different positions. It can be used to evaluate other teams special teams players and potentially find weaknesses. Also, it can be easily adjusted for any type of play, not only punt returns.

I found players in the top 10% of positive contribution seem to make up the top 20 players after adjusting for snap count. Ideally, we only want to include players that consistently play. It is interesting to see that we find players of all types on both offense and defense, not only punters who have the highest snap counts. With over 2000 players in the data set and 5000+ plays, it is encouraging to see the model found players that seemingly have impact while having minimal snaps. These players could be difference makers when it comes to special teams.  

#### Future Work

1. Future work should apply this model to other play types. In theory, one should be able to apply this model to any play since the same 22 players are always on the field during a play. A new success metric would need to be quantified for other play types.
2. Future work should include analysis on fumbles. This analysis removed fumbles because the Field Position variable could not be calculated when the punting team recovered the ball. This is a major drawback of this analysis since in reality fumbles on special teams plays are some of the most impactful and can instantly change a game. 
3. Future work should build cumulative team totals to assess which teams have best special teams in terms of field position. This could be useful in ranking teams. 
4. Future work should consider individual players salaries in an attempt to find undervalued players on special teams. These players could be key players in the draft, during free agency, and in trades. 

In conclusion, RAFPC can be a useful metric for ranking players across different positions for their contribution during punt plays with regard to relative field position.

# Appendix

Data cleaning and preliminary EDA can be found on [github](https://github.com/mattymo18/2022-NFL-Big-Data-Bowl). 


