---
title: "NFL Analytics"
author: "Alan Hartung"
date: "January 9, 2019"
output: html_document
---


## Intro

This competition peaked my interest, because I have a unique combination of data science background and football experience. In my five years playing college football competitively, I received two concussions: one after a helmet-to-helmet hit while catching a pass over the middle; one while blocking as I lead with my head to initiate a block. If the rules as regulation today were implied then, such as hitting a defenseless player with excessive force and lowering the head to initiate contact, my injuries may have been prevented.  

## Analysis Approach

Special teams are overlooked by most, however coaches and players know special team plays hold critical opportunities within the game to change a team’s momentum or. My analysis supports changes both from a player-safety standpoint and from a competitive aspect.  The primary based on the high rate of penalties and injuries on punt plays. 

When ask about the prevalence of penalties on punt plays, Rich McKay, the head of the competition committee, said.  "We know our fans are not big fans of every play ending in a penalty, and that's one thing we definitely want to look at,"

Illegal blocks above the waist and offensive holding are considered along with concussions.  

I chose to include Illegal blocks above the waist and offensive holding as other events to consider in player safety.  Most of these penalty events will not end up in a concussion.  However, from my experience, these behaviors and techniques are leading indicators for other serious player injuries.

The prevalence of these events will be examined through metrics including: speed of players, event time durations, and distance from returner. 

I was looking for the most subtle rule enhancement that could be implemented that would reduce concussions, reduce penalties, and preserve gameplay competition.  

A link to my most recent analysis presentation
[Presentation On Google Docs](https://drive.google.com/open?id=1nG-h5B3oExqjaeQKptYJ4meYzDXsCarJ)


# Rule Enhancement Suggestions
##Suggested Rule Enhancement # 1

RULE 9 SCRIMMAGE KICK
SECTION 1 KICK FROM SCRIMMAGE
ARTICLE 3.  DEFENSIVE TEAM FORMATION.

1. Enhancement - A maximum of one player for Team A and Team B may line up between the inbounds lines and the yard-line number on each side of the field. Remaining players must line up between each yard-line numbers.


## Suggested Rule Enhancement # 2

RULE 9 SCRIMMAGE KICK 
SECTION 1 KICK FROM SCRIMMAGE
ARTICLE 3.  DEFENSIVE TEAM FORMATION. 

2. Enhancement - Team A player lined up inside the numbers can not initiate contact with Team B player lined up between the inbounds lines and the yard-line number until after 15 yards past the original line of scrimmage.


## Potential Rule Enhancement Next Year

RULE 10 
SECTION 2 FAIR CATCH 
ARTICLE 4 PUTTING BALL IN PLAY AFTER FAIR CATCH

Enhancement 3 - After a fair catch is made, the receiving team will be awarded an additional 5 yards from the spot of the fair catch.




# Start of Analysis
## Initial Setup - Importing the Data 

First I look to import the data and standardize it a bit.  After some initial inspection of the data, I chose to standardize a few things like: lower case, time formats, etc

When I inspected the NGS Data, there are many observations for each player.  Observations before the play and after the play.  I chose to focus on observations that occurred between the snap of the play (event = 'ball_snap' and the last event recorded before the play ended. (Event with the maximum time that is not 'play_submit' )  
This cut down the overall observations to a more manageable dataset.  Again, the primary focus of my analysis occurs between the snap and 4 or 5 seconds after the returner has caught the ball.

I wrote a little function to handle this for each file.


```{r echo=FALSE}
options(contrasts = rep("contr.treatment", 2))
options(stringsAsFactors = FALSE)

suppressMessages(suppressWarnings({
  library("data.table")
  library("dplyr")
  library("stringr")
  library("tidytext")
  library("ggplot2")
  library("fasttime")
  library("scales")
  library("MASS")
}))

# I prefer to standardize a few things first, like lower case characters, names and using data.table
# This will help when merging tables later

data_import <- function(filename){
  newdata  <- fread(filename)
  newdata <- mutate_if(newdata,.predicate = is.character, .funs=tolower)
  names(newdata) <- tolower(names(newdata))
  setDT(newdata)
  newdata
}

player_punt_data <- data_import(filename = "../input/player_punt_data.csv")
player_role_data <- data_import(filename ="../input/play_player_role_data.csv")
play_information <- data_import(filename ="../input/play_information.csv")
game_data <- data_import(filename ="../input/game_data.csv")
video_review <- data_import(filename ="../input/video_review.csv")
video_footage <- data_import(filename ="../input/video_footage-injury.csv")

# Next Gen Data has millions of observations 

ngs_import <- function(filename){
  
  newdata  <- fread(filename,colClasses = c("Event"="character"))
  newdata <- mutate_if(newdata,.predicate = is.character, .funs=tolower)
  names(newdata) <- tolower(names(newdata))
  setDT(newdata)
  
  newdata[,time_num := fastPOSIXct(time)]
  newdata[event=="",event := NA_character_]
  
  events <- newdata[!is.na(event),]
  snap_time <- events[event=='ball_snap',list(time_num_snap = min(time_num)), by=list(season_year,gamekey,playid)]
  event_time_max <- events[ event!='play_submit', list(time_num_max = max(time_num)),by=list(season_year,gamekey,playid)]
  
  newdata <- merge(newdata,snap_time,by=c("season_year","gamekey","playid"))
  newdata <- merge(newdata,event_time_max,by=c("season_year","gamekey","playid"))
  newdata <- newdata[time_num >= time_num_snap & time_num <= time_num_max , ]
  
  newdata
}

ngs_16_pre <- ngs_import(filename = "../input/NGS-2016-pre.csv")
ngs_16_reg1 <- ngs_import(filename ="../input/NGS-2016-reg-wk1-6.csv")
ngs_16_reg2 <- ngs_import(filename ="../input/NGS-2016-reg-wk7-12.csv")
ngs_16_reg3 <- ngs_import(filename ="../input/NGS-2016-reg-wk13-17.csv")
ngs_16_post <- ngs_import(filename ="../input/NGS-2016-post.csv")
ngs_17_pre <- ngs_import(filename ="../input/NGS-2017-pre.csv")
ngs_17_reg1 <- ngs_import(filename ="../input/NGS-2017-reg-wk1-6.csv")
ngs_17_reg2 <- ngs_import(filename ="../input/NGS-2017-reg-wk7-12.csv")
ngs_17_reg3 <- ngs_import(filename ="../input/NGS-2017-reg-wk13-17.csv")
ngs_17_post <- ngs_import(filename ="../input/NGS-2017-post.csv")

list_nextgen <- list(ngs_16_pre,ngs_16_reg1,ngs_16_reg2,ngs_16_reg3,ngs_16_post,ngs_17_pre,ngs_17_reg1,ngs_17_reg2,ngs_17_reg3,ngs_17_post)
rm(list = ls(pattern = "ngs"))
setattr(list_nextgen, 'names', c("pre","reg","reg","reg","post","pre","reg","reg","reg","post"))
ngs_all <- rbindlist(list_nextgen, use.names=TRUE, fill=TRUE, idcol="season_type")
rm(list_nextgen)


```


## Data Wrangling

Here I will look at the player_punt_data, video_review, and do some data wrangling

```{r}


player_punt_data[, number := as.numeric(str_replace_all(number,"[a-z]",""))]
video_review[, primary_partner_gsisid := as.integer(primary_partner_gsisid) ]
player_punt_data <- unique(player_punt_data ,by="gsisid" )
video_review[,primary_partner_gsisid := as.integer(primary_partner_gsisid)]

setnames(video_footage,"season","season_year")
video <- merge(video_footage,video_review,by=c("season_year","gamekey","playid"),all.x = TRUE)

video <- merge(video,player_role_data,by=c("season_year","gamekey","playid","gsisid"),all.x = TRUE)
setnames(video,"role","primary_role")

# Doing the same for concussion partner
video <- merge(video,player_role_data,by.x =c("season_year","gamekey","playid","primary_partner_gsisid"),
               by.y =c("season_year","gamekey","playid","gsisid"),all.x = TRUE)
setnames(video,"role","partner_role")

video <- merge(video,unique(player_punt_data ,by='gsisid'),by=c("gsisid"))
video[,concussion := 1]


game_data <- game_data[,list(season_year,gamekey,game_date,game_day,turf,stadiumtype,temperature)]
game_data[, date := as.Date(game_date)]
game_data[,month := as.character(month(date))]

game_data[ str_detect(turf,'grass|natur'), turf_type := "grass"]
game_data[ turf %in% c('astroturf gameday grass 3d'), turf_type := "turf"]
game_data[is.na(turf_type), turf_type := "turf"]
game_data[stadiumtype=="", stadiumtype := NA_character_]

game_data[ , stadium_type := if_else(str_detect(stadiumtype,'outdoor')==TRUE,'uncontrolled','controlled','uncontrolled')]
game_data[,turf := NULL]
game_data[,game_date := NULL]
game_data[,date := NULL]
game_data[is.na(temperature) & stadium_type=='controlled', temperature := 68]

```


## All Punt Plays Data

This dataset has all punt plays.  I will look to build features and extract game situations and conditions for each punt play. This will help give insight into whether particular game circumstances play a role with concussions, penalties, and outcomes

I want to standardize the yard line so we know location of line of scrimmage.
NYG 40 yrd line translates into a different position on the field depending on who has possession of the ball when the ball is snapped for the punt.  Essentially, looking to see if the punter has 60 yards to punt until he reaches the end zone or 40 yards.  Converting to 0 to 100 yards

Using the possession team, make a defense / receiving team
For validation, we'll look at a histogram of punt line of scrimmage


```{r}
all_punts <- merge(play_information,video[,list(season_year,gamekey,playid,concussion)],
                   by=c("season_year","gamekey","playid"),all.x = TRUE)
all_punts[is.na(concussion), concussion := 0]

all_punts[,def_team := str_replace(home_team_visit_team,poss_team,'')]
all_punts[,def_team := trimws( str_replace(def_team,'-',''))]

# If possession team has the ball on their side of the 50, then use the yard line, otherwise add yardage over 50
all_punts[str_detect(yardline,poss_team), yardline_num := as.numeric(str_replace(yardline,poss_team,'')) ]
all_punts[str_detect(yardline,def_team) & is.na(yardline_num), yardline_num := (50-as.numeric(str_replace(yardline,def_team,''))) + 50 ]

all_punts[,yardline_grp := cut(yardline_num,breaks = c(-Inf,20,50,Inf)) ]

ggplot(data=all_punts, aes(yardline_num)) +
  ggtitle(label = "Punt Location on Field") +
  geom_histogram(binwidth = 1,aes(y=..density..), color="black", fill="white") +
  geom_vline(aes(xintercept=mean(yardline_num)),color="blue", linetype="dashed", size=1) +
  geom_vline(aes(xintercept=mean(yardline_num) - sd(yardline_num)),color="blue", linetype="dashed", size=1) +
  geom_vline(aes(xintercept=mean(yardline_num) + sd(yardline_num)),color="blue", linetype="dashed", size=1) +
  geom_density(alpha=.2, fill="#FF6666") +
  scale_x_continuous(name = "Yard Line" , breaks = seq(0,100,5),limits = c(0,100)) +
  scale_y_continuous(name = "Punt Play Density") 

```


Most punts happen between 20 and 50 yard line.  The average line scrimmage is at the 34 yard line.

## Information from play description 
The play description has tons of useful information. I will look to build more play features like:
Penalties, Punted Yards, Fair Catch, Touchback, Yard Returned, etc

I like using the this function first to get a sense of common words or combinations of words
Using tidytext, Let’s look at some common 3 word combinations to make some features

```{r}
library("tidytext")
trigram_data1 <- all_punts %>% unnest_tokens(trigram, playdescription, token = "ngrams", n = 3) %>% 
  count(trigram, sort = TRUE)
head(trigram_data1)
```


Using some of the events from description, lets build some features
First I'll set up some feature defaults and then scan the desription for key word
Also, using some regex to get punt yards

```{r }
all_punts[, `:=`(fair_catch = 0, touchback = 0, punt_yards = 0)]
# Removing commas from description for ease of use with extraction
all_punts[, playdescription := str_replace_all(playdescription,",",'')]
# Plays with a fair catch
all_punts[str_detect(playdescription,"fair catch"), fair_catch := 1]#
# Plays with touchback
all_punts[str_detect(playdescription,"touchback"), touchback := 1]
# Extract punted yards 
all_punts[, punt_yards :=
            as.numeric(trimws(str_replace_all(str_extract(playdescription,"punts\\s?\\d+\\s?yards"),'punts|yards','')))]
# Get overall penalty
all_punts[, penalty := if_else(str_detect(playdescription,'penalty')==TRUE,1,0,0) ]
#Look for offensive holding
all_punts[, penalty_off_holding := if_else(str_detect(playdescription,'offensive holding')==TRUE,1,0,0) ]
#Look for illegal block above the waist
all_punts[, penalty_ill_block := if_else(str_detect(playdescription,
                                                    'illegal block above the waist|illegal blindside')==TRUE,1,0,0) ]
# Major penalties are illegal block above the waist and holding
all_punts[, penalty_ill_block_holding := if_else(str_detect(playdescription,
                                        'illegal block above the waist|offensive holding|illegal blindside')==TRUE,1,0,0) ]

all_punts[, minute := as.numeric(str_sub(game_clock,1,2)) ]
all_punts[, total_minute := (quarter*15) - minute]

all_punts[,home_team :=str_replace(str_extract(home_team_visit_team,"\\w+\\s?(\\-)"),'-','')]
all_punts[,visit_team :=str_replace(str_extract(home_team_visit_team,"(\\-)\\s?\\w+"),'-','')]

# make a home team score and visiting team score
all_punts[,home_team_score := as.numeric(str_replace( str_extract(score_home_visiting,"\\d+\\s?(\\-)"),'-',''))]
all_punts[,visit_team_score := as.numeric(str_replace( str_extract(score_home_visiting,"(\\-)\\s?\\d+"),'-',''))]

all_punts[,poss_team_score := ifelse(poss_team==home_team,home_team_score,visit_team_score) ]
all_punts[,def_team_score := ifelse(def_team==home_team,home_team_score,visit_team_score) ]
all_punts[,poss_team_score_differential := poss_team_score - def_team_score ]
all_punts[,poss_team_score_differential_grp := cut(poss_team_score_differential,breaks = c(-Inf,-17,-9,-4,-1,0,3,8,17,Inf),
   labels = c("down by 17+","down by 9 tp 16","down by 4 to 8","down by 3 or less","tied","up by 3 or less","up by 4 to 8","up by 9 to 16","up by 17+")) ]

all_punts[,return_yards := as.numeric(str_replace_all(str_extract(playdescription,'for\\s?\\d+\\s?yards|for\\s?\\d+\\s?yds'),'for|yards|yds',''))]
all_punts[is.na(return_yards), return_yards := 0]



ggplot(data=all_punts[punt_yards %between% c(10,80),], aes(punt_yards)) +
  ggtitle(label = "Distance of Punt") +
  geom_histogram(binwidth = 1,aes(y=..density..), color="black", fill="white") +
  geom_vline(aes(xintercept=mean(punt_yards)),color="blue", linetype="dashed", size=1) +
  geom_vline(aes(xintercept=mean(punt_yards) - sd(punt_yards)),color="blue", linetype="dashed", size=1) +
  geom_vline(aes(xintercept=mean(punt_yards) + sd(punt_yards)),color="blue", linetype="dashed", size=1) +
  geom_density(alpha=.2, fill="#FF6666") +
  scale_x_continuous(name = "Yards" , breaks = seq(10,80,5),limits = c(10,80) ) +
  scale_y_continuous(name = "Punt Play Density") 

```

Most punts are 37 to 54 yards long


## High Level Insights

Here are some insights looking at Rates for Concussion, Touchback, Fair Catch, Illegal Block in the Back, and
Offensive Holding Rates by different groupings.


```{r echo=FALSE}


#############################################################################
##### Graph
temp <- all_punts[,list(
                #concussion_rate = mean(concussion)*100,
                fair_catch_rate = mean(fair_catch)*100,
                touchback_rate = mean(touchback)*100,
                illegal_block_rate = mean(penalty_ill_block)*100,
                offensive_holding_rate = mean(penalty_off_holding)*100),by=list(season_type)]

temp <- melt(temp,id.vars = "season_type",value.name = "rate",variable.name = "rate_metric")
ggplot(data=temp, aes(x=season_type,y=rate,fill=rate_metric)) +
  ggtitle(label = "Metrics by Season Type") +
  geom_bar(stat = "identity", color="blue",position=position_dodge()) +
  geom_text(aes(label=sprintf("%1.2f%%",rate) ), vjust=-1, color="black",
            position = position_dodge(.95), size=3) +
   scale_fill_brewer(palette="Paired") +
  scale_y_continuous(name = "Percent" ,limits = c(0,max(temp$rate)*1.1)) +
  xlab("Season Type")+
  theme_minimal()+
  labs(fill="Rate Metric")


#############################################################################
##### Graph
temp <- all_punts[,list(
                concussion_rate = mean(concussion)*100
                #fair_catch_rate = mean(fair_catch)*100,
                #touchback_rate = mean(touchback)*100,
                #illegal_block_rate = mean(penalty_ill_block)*100,
                #offensive_holding_rate = mean(penalty_off_holding)*100
  ),by=list(season_type)]

temp <- melt(temp,id.vars = "season_type",value.name = "rate",variable.name = "rate_metric")
ggplot(data=temp, aes(x=season_type,y=rate,fill=rate_metric)) +
  ggtitle(label = "Metrics by Season Type") +
  geom_bar(stat = "identity", color="blue",position=position_dodge()) +
  geom_text(aes(label=sprintf("%1.2f%%",rate) ), vjust=-1, color="black",
            position = position_dodge(.95), size=3) +
   scale_fill_brewer(palette="Paired") +
  scale_y_continuous(name = "Percent" ,limits = c(0,max(temp$rate)*1.1)) +
    xlab("Season Type")+
  theme_minimal()+
  labs(fill="Rate Metric")



#############################################################################
##### Graph
temp <- all_punts[,list(
                #concussion_rate = mean(concussion)*100,
                fair_catch_rate = mean(fair_catch)*100,
                touchback_rate = mean(touchback)*100,
                illegal_block_rate = mean(penalty_ill_block)*100,
                offensive_holding_rate = mean(penalty_off_holding)*100
                ),by=list(yardline_grp)]

temp <- melt(temp,id.vars = "yardline_grp",value.name = "rate",variable.name = "rate_metric")
ggplot(data=temp, aes(x=yardline_grp,y=rate,fill=rate_metric)) +
  ggtitle(label = "Metrics by Punting Team Line of Scrimmage") +
  geom_bar(stat = "identity", color="blue",position=position_dodge()) +
  geom_text(aes(label=sprintf("%1.2f%%",rate) ), vjust=-1, color="black",
            position = position_dodge(1), size=3) +
   scale_fill_brewer(palette="Paired") +
  scale_y_continuous(name = "Percent" ,limits = c(0,max(temp$rate)*1.1)) +
      xlab("Yard Line Group") +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))+
  labs(fill="Rate Metric")


#############################################################################
##### Graph
temp <- all_punts[,list(
                concussion_rate = mean(concussion)*100
                #fair_catch_rate = mean(fair_catch)*100,
                #touchback_rate = mean(touchback)*100,
                #illegal_block_rate = mean(penalty_ill_block)*100,
               # offensive_holding_rate = mean(penalty_off_holding)*100
                ),by=list(yardline_grp)]

temp <- melt(temp,id.vars = "yardline_grp",value.name = "rate",variable.name = "rate_metric")
ggplot(data=temp, aes(x=yardline_grp,y=rate,fill=rate_metric)) +
  ggtitle(label = "Metrics by Punting Team Line of Scrimmage") +
  geom_bar(stat = "identity", color="blue",position=position_dodge()) +
  geom_text(aes(label=sprintf("%1.2f%%",rate) ), vjust=-1, color="black",
           position = position_dodge(1), size=3) +
   scale_fill_brewer(palette="Paired") +
  scale_y_continuous(name = "Percent" ,limits = c(0,max(temp$rate)*1.1)) +
  xlab("Yard Line Group") +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))+
  labs(fill="Rate Metric")


```


## Some insights here :

1. Concussion Rates are higher in the pre-season.  This make sense since there are more in-experienced players trying to make the team. NFL teams carry 90 or so players during the pre-season but only 53 during the regular and post-season. These extra players are trying to make the "Big Play" on special teams since that is where most of the open spots tend to be initially for new players.  Sometimes "Big Play" = High Risk. 

2. As Fair Catch and Touch Back Rates increase, concussion and penalty rates decrease.  This is more a function of the line of scrimmage.  Once the offense reaches the 40 to 50 yard line, space is limited to attempt a punt return. I plan to show additional information with the NGS data that take player distance from punt returner into account

3. Major punt penalties such as offensive holding and illegal block in the back decreases as the line of scrimmage moves closer to the opposing teams goal line.  This is due to the amount of space available to return and the position of the players relative to punt returner. As well as focusing on concussions, I will include analysis to decrease the amount of penalties on punt plays. 
When ask about the prevalence of penalties on punt plays, Rich McKay, the head of the competition committee, said.  "We know our fans are not big fans of every play ending in a penalty, and that's one thing we definitely want to look at,"
I will show additional insights with NGS data that may guide rule enhancements.

4. Interestingly enough, punts within the 5 yard line (punting out of their own end zone) show a 0 concussion rate.  This is most likely a function of a sample size of 100 or so, however the spacing for the punt team is condensed.  Punters typically position themselves about 14 yards behind the line of scrimmage.  This may trigger different formations as well as a quicker snap and punt to get a successful punt off.


## Player and Team Roles

After inspecting the player_role_data , the "role" of each id is very useful. I chose to categorize these roles into
which team they are on " Punting Team" or "Receiving Team" = role_team.
Also, I chose to categorize player roles by their primary assignment and formation alignment. 
"Gunner","Backfield","Line" = role_pos. I will use these player grouping to generalize metrics and behaviors.

```{r echo=FALSE}

player_role_data[str_detect(role,"gr|gl"), `:=`( role_team = "pnt" , role_pos = "gunner") ]
player_role_data[str_detect(role,"vr|vl"), `:=`( role_team = "rec" , role_pos = "gunner") ]
player_role_data[str_detect(role,"pd"), `:=`( role_team = "rec" , role_pos = "line") ]
player_role_data[role %in% c("pr"), `:=`( role_team = "rec" , role_pos = "returner") ]
player_role_data[role %in% c("p"), `:=`( role_team = "pnt" , role_pos = "punter") ]
player_role_data[str_detect(role,"pll|plm|plr|pfb"), `:=`( role_team = "rec" , role_pos = "backfield") ]
player_role_data[str_detect(role,"pp|plw|prw|pc"), `:=`( role_team = "pnt" , role_pos = "backfield") ]
player_role_data[role %in% c("plt","plg","pls","prg","prt"), `:=`( role_team = "pnt" , role_pos = "line") ]
player_role_data[str_detect(role,"l"), role_side := "left" ]
player_role_data[str_detect(role,"r"), role_side := "right" ]
player_role_data[str_detect(role,"pls"), role_side := "middle" ]
player_role_data[is.na(role_side), role_side := "middle" ]


```


# NGS Data

Using Next Gen Stats data by player will provide the best insight for proposed rule enhancements

Main approach with NGS Data:
I will look to build features that I can summarize to a play level.  Within the NGS data, there are events flagged to indicate what punt event just happened.  I will use these events, along with other markers to build metrics.

Here are the events I will consider:


1. Snap of the Punt Play

2. Punt - when the punter punts the ball

3. I will make other "events" to indicate the time since the ball was punted.  These will be in 1 second intervals

4. Fair Catch, Punt Received, Punt Muffed, etc.  Looking to get the event when the punt is near the ground.

5. Similar to #3, extend the 1 second intervals to after the catch

Essentially I would like to know these attributes at each "event"

1. Speed of the player

2. Distance from the punt returner

3. Time Duration between major events

I will use this information to exhibit a few punt play dynamics that increase the likelihood of concussions and major punt penalties 

```{r echo=FALSE}

master <- copy(all_punts)

ngs_events_playid <- ngs_all[!is.na(event),list(time_num = min(time_num)), by=list(season_year,gamekey,playid,event)]
t1 <-  ngs_events_playid[,.N,by=list(event)]

# Get the time of the punt to append for each play
dt_time_punt <- ngs_events_playid[event=='punt',list(time_num_punt = min(time_num)), by=list(season_year,gamekey,playid)]

# Get the time when the put was near the ground or caught
dt_time_punt_after <- ngs_events_playid[event %in% c("fair_catch","punt_received","punt_land","punt_muffed","punt_downed","touchback","out_of_bounds"),
                                list(time_num_punt_after = min(time_num)), by=list(season_year,gamekey,playid)]


# Make a dataset that has generated time events.  NGS data has snap, punt, fair catch etc.
# Looking make events based on seconds from punt
# Get the time events every second after punt to calculate metrics
dt_time_punt1 <- ngs_events_playid[event=='punt',][, `:=`( time_num = time_num+1 , event = 'punt1')][]
dt_time_punt2 <- ngs_events_playid[event=='punt',][, `:=`( time_num = time_num+2 , event = 'punt2')][]
dt_time_punt3 <- ngs_events_playid[event=='punt',][, `:=`( time_num = time_num+3 , event = 'punt3')][]
dt_time_punt4 <- ngs_events_playid[event=='punt',][, `:=`( time_num = time_num+4 , event = 'punt4')][]
dt_time_punt5 <- ngs_events_playid[event=='punt',][, `:=`( time_num = time_num+5 , event = 'punt5')][]
dt_time_punt6 <- ngs_events_playid[event=='punt',][, `:=`( time_num = time_num+6 , event = 'punt6')][]
dt_time_punt7 <- ngs_events_playid[event=='punt',][, `:=`( time_num = time_num+7 , event = 'punt7')][]
dt_time_punt8 <- ngs_events_playid[event=='punt',][, `:=`( time_num = time_num+8 , event = 'punt8')][]
dt_time_punt9 <- ngs_events_playid[event=='punt',][, `:=`( time_num = time_num+9 , event = 'punt9')][]



ngs_events_playid <- rbind(ngs_events_playid,dt_time_punt1,dt_time_punt2,dt_time_punt3,dt_time_punt4,
                                 dt_time_punt5,dt_time_punt6,dt_time_punt7,dt_time_punt8,dt_time_punt9)
ngs_events_playid[,new_event := event]
ngs_events_playid <- unique(ngs_events_playid,by=c("season_year","gamekey","playid","time_num"))

ngs_events_playid[event %in% c("fair_catch","punt_received","punt_land","punt_muffed"), new_event := "punt_after"]
rm(dt_time_punt1,dt_time_punt2,dt_time_punt3,dt_time_punt4,dt_time_punt5,dt_time_punt6,dt_time_punt7,dt_time_punt8,dt_time_punt9)

# Merge in the new events to the NGS data
ngs_all_new <- merge(ngs_all,ngs_events_playid[,list(season_year,gamekey,playid,time_num,new_event)],by=c("season_year","gamekey","playid","time_num"),all.x = TRUE)

#Keep only the NGS obsservations from events
ngs_all_new <- ngs_all_new[new_event %in% c('ball_snap','punt','punt1','punt2','punt3','punt4',
                                            'punt5','punt6','punt7','punt8','punt9','punt_after')]

#Bring in the Punt time and Punt After times for each play
ngs_all_new <- merge(ngs_all_new,dt_time_punt,by=c("season_year","gamekey","playid"))
ngs_all_new <- merge(ngs_all_new,dt_time_punt_after,by=c("season_year","gamekey","playid"))


#Bring in Role data
ngs_all_new <- merge(ngs_all_new, player_role_data,by=c("season_year","gamekey","playid","gsisid"))
#Want to know where the returner is and how other player are located relative to the returner by each time event
# Looking to calculate player distance from returner
ret_loc <- ngs_all_new[role_pos %in% c("returner"), list(season_year,gamekey,playid,new_event,time_num,x,y)]
setnames(ret_loc,c("x","y"),c("ret_x","ret_y"))
ret_loc <- unique(ret_loc)

# Merge the retuner location at each event back to the data
ngs_all_new <- merge(ngs_all_new, ret_loc,by=c("season_year","gamekey","playid","new_event","time_num"),all.x = TRUE)

# Making some metrics like speed, event duration, and distance
ngs_all_new[,speed_meter := (dis/.9144)/.1]
ngs_all_new[,dist_from_ret := sqrt((x-ret_x)^2 + (y-ret_y)^2)]
ngs_all_new[,snap_duration := as.numeric(difftime(time_num_punt,time_num_snap,units = "sec"))]
ngs_all_new[,hang_time_duration := as.numeric(difftime(time_num_punt_after,time_num_punt,units = "sec"))]
ngs_all_new[,change_pos_duration := snap_duration+hang_time_duration]
ngs_all_new[,time_since_punt := as.numeric(difftime(time_num,time_num_punt,units = "sec"))]

ngs_all_new <- ngs_all_new[order(season_year,gamekey,playid,time_num)]

# Get the count of players by role position at the time of the snap
snap_position <- ngs_all_new[event=='ball_snap' & !is.na(role_team),]
snap_position <- snap_position[,list(count=as.numeric(.N)),by=list(season_year,gamekey,playid,role_team,role_pos)]
snap_position <- melt(snap_position,id.vars = c("season_year","gamekey","playid","role_team","role_pos"))
snap_position <- dcast(snap_position, season_year+gamekey+playid ~ paste0("snapform_",role_team,"_",role_pos,"_",variable),value.var = "value",fun.aggregate = mean)
snap_position[is.na(snap_position)] <- 0

ngs_all_new <- merge(ngs_all_new,snap_position, by=c("season_year","gamekey","playid"))

duration_metrics <- ngs_all_new[,list(snap_duration=min(snap_duration,na.rm = TRUE),
                                      hang_time_duration=min(hang_time_duration,na.rm = TRUE),
                                      change_pos_duration=min(change_pos_duration,na.rm = TRUE)) ,by=list(season_year,gamekey,playid)]



speed_data <- dcast(ngs_all_new, season_year+gamekey+playid ~ str_c("speed",new_event,role_team,role_pos,sep = "_"), value.var = "speed_meter",fun.aggregate = mean,fill = NA)
distance_data <- dcast(ngs_all_new, season_year+gamekey+playid ~ str_c("distance",new_event,role_team,role_pos,sep = "_"), value.var = "dist_from_ret",fun.aggregate = mean,fill = NA)
speed_data[] <- lapply(speed_data, function(x) ifelse(is.na(x), median(x, na.rm = TRUE), x))
distance_data[] <- lapply(distance_data, function(x) ifelse(is.na(x), median(x, na.rm = TRUE), x))

final <- merge(master,speed_data,by=c("season_year","gamekey","playid"),all.x = TRUE)
final <- merge(final,distance_data,by=c("season_year","gamekey","playid"),all.x = TRUE)
final <- merge(final,snap_position,by=c("season_year","gamekey","playid"),all.x = TRUE)
final <- merge(final,duration_metrics,by=c("season_year","gamekey","playid"),all.x = TRUE)

```


## Coverage Waves

Conventional football wisdom would likely conclude that the distance from the punt returner after 3 seconds or so would give good indication on whether a punt returner will attempt to return the ball or call for a fair catch.

The punt returners decision is based on how close the opposing team is, how fast they are approaching, and if they are unblocked.  


```{r echo=FALSE}

pop <- final[!(season_type %in% c("pre")) & yardline_num %between% c(0,49) ,]
temp <- merge(pop[,list(season_year,gamekey,playid)],ngs_all_new,by=c("season_year","gamekey","playid"))
temp[role_team=='pnt', role_team :='Punting Team']
temp[role_team=='rec', role_team :='Receiving Team']
temp[role_pos=='gunner', role_pos :='Gunner']
temp[role_pos=='line', role_pos :='Line / Backfield']
temp[role_pos=='backfield', role_pos :='Line / Backfield']
t1<- temp[ ,list(metric= mean(dist_from_ret,na.rm = TRUE),
                 time = mean(time_since_punt,na.rm = TRUE), obs=.N),by=list(new_event,role_team,role_pos)]

total_plays <- pop[,.N]

ggplot(data=t1[!(role_pos %in% c("punter","returner")) & time > -2.5,], aes(x = time, 
                                              color = paste0(role_team," : ",role_pos) ) ) +
  ggtitle(label = "Avg Distance by Punt Team Unit (mps)" , 
          subtitle = paste0("Sample of ",total_plays," Punt Plays") ) +
  stat_smooth(aes(y=metric, x=time), method = lm, formula = y ~ poly(x,4), se = FALSE) +
  scale_linetype_manual(name="Team" , values=c("solid","dotdash")) +
  scale_color_manual(name="Position", values=c("red","green","blue","yellow")) +
  geom_vline(aes(xintercept=-2),color="black", linetype="dotted", size=.75) +
  geom_vline(aes(xintercept=0),color="black", linetype="dotted", size=.75) +
  geom_vline(aes(xintercept=mean(pop$hang_time_duration,na.rm = TRUE)),color="black", linetype="dotted", size=.75) +
  scale_x_continuous(name = "Time (sec.) Since Punt" , breaks = seq(-3,8,1)) +
  scale_y_continuous(name = "Player Distance from Punt Returner (yds)") +
  geom_text(aes(y=max(t1$metric)*1.1, x=-1.5 ,label="Snap"), color="black", size=3.5) +
  geom_text(aes(y=max(t1$metric)*1.1, x=.5 ,label="Punt"), color="black", size=3.5) +
  geom_text(aes(y=max(t1$metric)*1.1, x=(mean(pop$hang_time_duration,na.rm = TRUE)-1) ,label="Punt Catch"), color="black", size=3.5) 


```


The Gunner Group are fist down the field, about 15 yards ahead of the line and back field

As you can see, the punt coverage units have 2 waves of units going down the field.
This is much different than kickoffs since they typically have 1 wave due to everyone moving down the field at the same time.


```{r echo=FALSE}


pop <- final[!(season_type %in% c("pre")) & yardline_num %between% c(0,49) ,]
temp <- merge(pop[,list(season_year,gamekey,playid)],ngs_all_new,by=c("season_year","gamekey","playid"))
temp[role_team=='pnt', role_team :='Punting Team']
temp[role_team=='rec', role_team :='Receiving Team']
temp[role_pos=='gunner', role_pos :='Gunner']
temp[role_pos=='line', role_pos :='Line / Backfield']
temp[role_pos=='backfield', role_pos :='Line / Backfield']
t1<- temp[ ,list(metric= mean(dist_from_ret,na.rm = TRUE),
                 time = mean(time_since_punt,na.rm = TRUE), obs=.N),by=list(new_event,role_team,role_pos)]

total_plays <- pop[,.N]

t1<- temp[ ,list(metric= median(speed_meter,na.rm = TRUE),
                        time = median(time_since_punt,na.rm = TRUE), obs=.N),by=list(new_event,role_team,role_pos)]

ggplot(data=t1[!(role_pos %in% c("punter","returner")),], aes(x = time,  color = paste0(role_team," : ",role_pos) ) ) +
  ggtitle(label = "Avg Speed by Punt Team Unit (mps) " , 
          subtitle = paste0("Sample of ",total_plays," Punt Plays") ) +
  stat_smooth(aes(y=metric, x=time), method = lm, formula = y ~ poly(x,4), se = FALSE) +
  scale_linetype_manual(name="Team" , values=c("solid","dotdash")) +
  scale_color_manual(name="Position", values=c("red","green","blue","yellow")) +
  geom_vline(aes(xintercept=-2),color="black", linetype="dotted", size=.75) +
  geom_vline(aes(xintercept=0),color="black", linetype="dotted", size=.75) +
  geom_vline(aes(xintercept=mean(pop$hang_time_duration,na.rm = TRUE)),color="black", linetype="dotted", size=.75) +
  scale_x_continuous(name = "Time (sec.) Since Punt" , breaks = seq(-3,8,1)) +
  scale_y_continuous(name = "Player Speed (mps)") +
  geom_text(aes(y=max(t1$metric)*1.1, x=-1.5 ,label="Snap"), color="black", size=3.5) +
  geom_text(aes(y=max(t1$metric)*1.1, x=.5 ,label="Punt"), color="black", size=3.5) +
  geom_text(aes(y=max(t1$metric)*1.1, x=(mean(pop$hang_time_duration,na.rm = TRUE)-1) ,label="Punt Catch"), color="black", size=3.5) 

```


Gunner Units max out their speed at around 2 seconds after the ball has been kicked.
Line and Back field unit max out their speed at around 4 seconds after the ball has been kicked.

With 2 different waves of players moving down the field moving at different rates of speed at different times is a chaotic environment to be in.

## Punt Team Gunners and Receiving Team Gunners

In the NFL, the two outside punting team players can advance downfield before the ball is punted.  The remaining punt team players must wait until the punter has punted the ball to advance down field.
(see image below)

Next I will examine the relationship of the punt team gunners on our key metrics.  For this analysis, I chose to focus on Receiving Team Gunner Counts.
There are varying levels on how to defend against the punt team gunners.
For this analysis, I chose to focus on Receiving Team Gunner Counts of 2, 3, and 4.

```{r echo=FALSE}

#############################################################################
##### Graph

temp <- final[!(season_type %in% c("pre")) & yardline_num >= 50,]
total_plays <- temp[,.N]
temp <- temp[,list(
                concussion_rate = mean(concussion)*100,
                fair_catch_rate = mean(fair_catch)*100,
                touchback_rate = mean(touchback)*100,
                illegal_block_rate = mean(penalty_ill_block)*100,
                offensive_holding_rate = mean(penalty_off_holding)*100
                #obs=.N
                ),by=list(snapform_rec_gunner_count)]

temp <- melt(temp[snapform_rec_gunner_count %between% c(2,4),],id.vars = "snapform_rec_gunner_count",value.name = "rate",variable.name = "rate_metric")
ggplot(data=temp, aes(x=snapform_rec_gunner_count,y=rate,fill=rate_metric)) +
  ggtitle(label = "Metrics by # of Receiving Team Gunner Coverage \n Punts past the 50 yard line",
          subtitle = paste0("Sample of ",total_plays," Reg/Post Season Punt Plays")) +
  geom_bar(stat = "identity", color="blue",position=position_dodge()) +
  geom_text(aes(label=sprintf("%1.2f%%",rate) ), vjust=-1, color="black",
           position = position_dodge(1), size=3) +
   scale_fill_brewer(palette="Paired") +
  scale_y_continuous(name = "Percent of Punt Plays" ,limits = c(0,max(temp$rate)*1.1)) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1)) + xlab("# of Receiving Team Gunner Coverage") +
  labs(fill="Rate Metric")

#############################################################################
##### Graph

temp <- final[!(season_type %in% c("pre")) & yardline_num %between% c(0,49),]
total_plays <- temp[,.N]
temp <- temp[,list(
                concussion_rate = mean(concussion)*100,
                fair_catch_rate = mean(fair_catch)*100,
                touchback_rate = mean(touchback)*100,
                illegal_block_rate = mean(penalty_ill_block)*100,
                offensive_holding_rate = mean(penalty_off_holding)*100
                #obs=.N
                ),by=list(snapform_rec_gunner_count)]

temp <- melt(temp[snapform_rec_gunner_count %between% c(2,4),],id.vars = "snapform_rec_gunner_count",value.name = "rate",variable.name = "rate_metric")
ggplot(data=temp, aes(x=snapform_rec_gunner_count,y=rate,fill=rate_metric)) +
  ggtitle(label = "Metrics by # of Receiving Team Gunner Coverage \n Punts between the 1 and 49 yard line",
          subtitle = paste0("Sample of ",total_plays," Reg/Post Season Punt Play")) +
  geom_bar(stat = "identity", color="blue",position=position_dodge()) +
  geom_text(aes(label=sprintf("%1.2f%%",rate) ), vjust=-1, color="black",
           position = position_dodge(1), size=3) +
   scale_fill_brewer(palette="Paired") +
  scale_y_continuous(name = "Percent of Punt Plays" ,limits = c(0,max(temp$rate)*1.1)) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1)) + xlab("# of Receiving Team Gunner Coverage")+
  labs(fill="Rate Metric")


```



#Insight #1
Appears to be a telling story here.  Punt plays that only have 2 Receiving Gunner Coverage players typically have decreased concussions and penalties while increasing the liklihood of a fair catch.

Lets look at the percentage increase for each additional Receiving Gunner Coverage Player


```{r echo=FALSE}

#############################################################################
##### Graph


temp <- final[!(season_type %in% c("pre")) & yardline_num %between% c(0,49),]
total_plays <- temp[,.N]
temp <- temp[,list(
  concussion_rate = mean(concussion)*100,
  fair_catch_rate = mean(fair_catch)*100,
  touchback_rate = mean(touchback)*100,
  illegal_block_rate = mean(penalty_ill_block)*100,
  offensive_holding_rate = mean(penalty_off_holding)*100
  #obs=.N
),by=list(snapform_rec_gunner_count)]
temp <- melt(temp[snapform_rec_gunner_count %between% c(2,4),],id.vars = "snapform_rec_gunner_count",value.name = "rate",variable.name = "rate_metric")
base_temp <- temp[snapform_rec_gunner_count==2,list(rate_metric,rate)]
setnames(base_temp,"rate","base_rate")
  
temp <- merge(temp,base_temp,by = "rate_metric")
temp1 <- temp[snapform_rec_gunner_count!=2,rate := ((rate-base_rate)/base_rate)*100][snapform_rec_gunner_count!=2,][]
temp1[is.nan(rate),rate := 0]
temp1[,snapform_rec_gunner_count := as.factor(snapform_rec_gunner_count)]
ggplot(data=temp1, aes(x=snapform_rec_gunner_count,y=rate,fill=rate_metric)) +
  ggtitle(label = "% Increase of Metrics with Additional Gunner Coverage \n Punts between the 1 and 49 yard line",
          subtitle = paste0("Sample of ",total_plays," Reg/Post Season Punt Play")) +
  geom_bar(stat = "identity", color="blue",position=position_dodge()) +
  geom_text(aes(label=sprintf("%1.0f%%",rate) ), vjust=-1, color="black",
            position = position_dodge(1), size=3) +
  scale_fill_brewer(palette="Paired") +
  scale_y_continuous(name = "Percent Increase : Baseline 2 Gunners" ,limits = c(min(temp1$rate),max(temp1$rate)*1.1)) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1)) + xlab("# of Receiving Team Gunner Coverage")+
  labs(fill="Rate Metric")


```


```{r echo=FALSE}

#############################################################################
##### Graph
pop <- final[!(season_type %in% c("pre")) & yardline_num %between% c(0,49) ,]
temp <- merge(pop[,list(season_year,gamekey,playid)],ngs_all_new,by=c("season_year","gamekey","playid"))
temp[role_team=='pnt', role_team :='Punting Team']
temp[role_team=='rec', role_team :='Receiving Team']
temp[role_pos=='gunner', role_pos :='Gunner']
temp[role_pos=='line', role_pos :='Line / Backfield']
temp[role_pos=='backfield', role_pos :='Line / Backfield']
t1<- temp[role_pos =='Gunner' & role_team=='Punting Team' & snapform_rec_gunner_count %between% c(2,4) ,list(metric= mean(dist_from_ret,na.rm = TRUE),
                 time = mean(time_since_punt,na.rm = TRUE), obs=.N),by=list(new_event,role_team,snapform_rec_gunner_count)]

total_plays <- pop[,.N]

ggplot(data=t1[ time > -2.5,], aes(x = time, color = as.factor(snapform_rec_gunner_count)) ) +
  ggtitle(label = "Avg Punt Team Gunner Distance from Punt Returner" , 
          subtitle = paste0("Sample of ",total_plays," Punt Plays") ) +
  stat_smooth(aes(y=metric, x=time), method = lm, formula = y ~ poly(x,4), se = FALSE) +
  scale_linetype_manual(name="Team" , values=c("solid","dotdash")) +
  scale_color_manual(name="Gunner Coverage Count", values=c("red","green3","blue")) +
  geom_vline(aes(xintercept=-2),color="black", linetype="dotted", size=.75) +
  geom_vline(aes(xintercept=0),color="black", linetype="dotted", size=.75) +
  geom_vline(aes(xintercept=mean(pop$hang_time_duration,na.rm = TRUE)),color="black", linetype="dotted", size=.75) +
  scale_x_continuous(name = "Time (sec.) Since Punt" , breaks = seq(-3,8,1)) +
  scale_y_continuous(name = "Punt Team Gunner Distance from Punt Returner (yds)") +
  geom_text(aes(y=max(t1$metric)*1.1, x=-1.5 ,label="Snap"), color="black", size=3.5) +
  geom_text(aes(y=max(t1$metric)*1.1, x=.5 ,label="Punt"), color="black", size=3.5) +
  geom_text(aes(y=max(t1$metric)*1.1, x=(mean(pop$hang_time_duration,na.rm = TRUE)-1) ,label="Punt Catch"), color="black", size=3.5) 

```


```{r echo=FALSE}

t1 <- t1[new_event=="ball_snap", time := -2.1]
t1 <- t1[new_event=="punt_after", time := 4.5]
t1 <- t1[, role_team := 'gunner']
t2 <- dcast(t1,  new_event+snapform_rec_gunner_count+time ~ role_team, value.var = "metric")
base_temp <- t2[snapform_rec_gunner_count==2,list(new_event,time, gunner)]
setnames(base_temp,"gunner","base_metric")

temp <- merge(t2,base_temp,by = c("new_event","time"))
temp1 <- temp[snapform_rec_gunner_count!=2,metric := (base_metric-gunner)][snapform_rec_gunner_count!=2,][]
temp1[is.nan(metric),metric := 0]
temp1[,snapform_rec_gunner_count := as.factor(snapform_rec_gunner_count)]

total_plays <- pop[,.N]
ggplot(data=temp1[,], aes(x = time, color = as.factor(snapform_rec_gunner_count)) ) +
  ggtitle(label = "Avg Difference Compared to Only 2 Gunner Coverage Count" , subtitle = paste0("Sample of ",total_plays," Punt Plays") ) +
  stat_smooth(aes(y=metric, x=time), method = lm, formula = y ~ poly(x,2), se = FALSE) +
  scale_linetype_manual(name="Team" , values=c("solid","dotdash")) +
  scale_color_manual(name="Gunner Coverage Count", values=c("red","green3","blue")) +
  geom_vline(aes(xintercept=-2),color="black", linetype="dotted", size=.75) +
  geom_vline(aes(xintercept=0),color="black", linetype="dotted", size=.75) +
  geom_vline(aes(xintercept=mean(pop$hang_time_duration,na.rm = TRUE)),color="black", linetype="dotted", size=.75) +
  scale_x_continuous(name = "Time (sec.) Since Punt" , breaks = seq(-3,8,1)) +
  scale_y_continuous(name = "Punt Team Gunner Distance from Punt Returner (yds)") +
  geom_text(aes(y=max(temp1$metric)*1, x=-1.5 ,label="Snap"), color="black", size=3.5) +
  geom_text(aes(y=max(temp1$metric)*1, x=.5 ,label="Punt"), color="black", size=3.5) +
  geom_text(aes(y=max(temp1$metric)*1, x=(mean(pop$hang_time_duration,na.rm = TRUE)-1) ,label="Punt Catch"), color="black", size=3.5) 

```


```{r echo=FALSE}

pop <- final[!(season_type %in% c("pre")) & yardline_num %between% c(0,49) ,]
temp <- merge(pop[,list(season_year,gamekey,playid)],ngs_all_new,by=c("season_year","gamekey","playid"))
temp[role_team=='pnt', role_team :='Punting Team']
temp[role_team=='rec', role_team :='Receiving Team']
temp[role_pos=='gunner', role_pos :='Gunner']
temp[role_pos=='line', role_pos :='Line / Backfield']
temp[role_pos=='backfield', role_pos :='Line / Backfield']
t1<- temp[role_pos =='Gunner' & role_team=='Punting Team' & snapform_rec_gunner_count %between% c(2,4) ,list(metric= mean(speed_meter,na.rm = TRUE),
                   time = mean(time_since_punt,na.rm = TRUE), obs=.N),by=list(new_event,role_team,snapform_rec_gunner_count)]

total_plays <- pop[,.N]

ggplot(data=t1[ time > -2.5,], aes(x = time, color = as.factor(snapform_rec_gunner_count)) ) +
  ggtitle(label = "Avg Punt Gunner Speed (mps) by Gunner Coverage Count" , subtitle = paste0("Sample of ",total_plays," Punt Plays") ) +
  stat_smooth(aes(y=metric, x=time), method = lm, formula = y ~ poly(x,4), se = FALSE) +
  scale_linetype_manual(name="Team" , values=c("solid","dotdash")) +
  scale_color_manual(name="Gunner Coverage Count", values=c("red","green3","blue")) +
  geom_vline(aes(xintercept=-2),color="black", linetype="dotted", size=.75) +
  geom_vline(aes(xintercept=0),color="black", linetype="dotted", size=.75) +
  geom_vline(aes(xintercept=mean(pop$hang_time_duration,na.rm = TRUE)),color="black", linetype="dotted", size=.75) +
  scale_x_continuous(name = "Time (sec.) Since Punt" , breaks = seq(-3,8,1)) +
  scale_y_continuous(name = "Punt Gunner Speed (mps)") +
  geom_text(aes(y=max(t1$metric)*1.1, x=-1.5 ,label="Snap"), color="black", size=3.5) +
  geom_text(aes(y=max(t1$metric)*1.1, x=.5 ,label="Punt"), color="black", size=3.5) +
  geom_text(aes(y=max(t1$metric)*1.1, x=(mean(pop$hang_time_duration,na.rm = TRUE)-1) ,label="Punt Catch"), color="black", size=3.5) 


```


```{r echo=FALSE}

t1 <- t1[new_event=="ball_snap", time := -2.1]
t1 <- t1[new_event=="punt_after", time := 4.5]
t1 <- t1[, role_team := 'gunner']
t2 <- dcast(t1,  new_event+snapform_rec_gunner_count+time ~ role_team, value.var = "metric")
base_temp <- t2[snapform_rec_gunner_count==2,list(new_event,time, gunner)]
setnames(base_temp,"gunner","base_metric")

temp <- merge(t2,base_temp,by = c("new_event","time"))
temp1 <- temp[snapform_rec_gunner_count!=2,metric := (gunner-base_metric)][snapform_rec_gunner_count!=2,][]
temp1[is.nan(metric),metric := 0]
temp1[,snapform_rec_gunner_count := as.factor(snapform_rec_gunner_count)]

total_plays <- pop[,.N]
ggplot(data=temp1[,], aes(x = time, color = as.factor(snapform_rec_gunner_count)) ) +
  ggtitle(label = "Avg Difference Compared to Only 2 Gunner Coverage Count" , subtitle = paste0("Sample of ",total_plays," Punt Plays") ) +
  stat_smooth(aes(y=metric, x=time), method = lm, formula = y ~ poly(x,2), se = FALSE) +
  scale_linetype_manual(name="Team" , values=c("solid","dotdash")) +
  scale_color_manual(name="Gunner Coverage Count", values=c("red","green3","blue")) +
  geom_vline(aes(xintercept=-2),color="black", linetype="dotted", size=.75) +
  geom_vline(aes(xintercept=0),color="black", linetype="dotted", size=.75) +
  geom_vline(aes(xintercept=mean(pop$hang_time_duration,na.rm = TRUE)),color="black", linetype="dotted", size=.75) +
  scale_x_continuous(name = "Time (sec.) Since Punt" , breaks = seq(-3,8,1)) +
  scale_y_continuous(name = "Punt Gunner Speed (mps)") +
  geom_text(aes(y=max(temp1$metric)*1, x=-1.5 ,label="Snap"), color="black", size=3.5) +
  geom_text(aes(y=max(temp1$metric)*1, x=.5 ,label="Punt"), color="black", size=3.5) +
  geom_text(aes(y=max(temp1$metric)*1, x=(mean(pop$hang_time_duration,na.rm = TRUE)-1) ,label="Punt Catch"), color="black", size=3.5) 

```


Overall Graph Conclusions on the Increase of Receiving Gunner Count
* Concussions increased
* Fair Catches decreased 
* Touchbacks decreased 
* Illegal Blocks increased 
* Offensive Holding increased


# Predicted Probabilities by Punt Line of Scrimmage, Punt Duration, and Receiving Team Gunner Count
 
 The above information is good discovery on where to focus the attention moving forward.
 Next I intend to predict the probability of concussions, major penalties, and fair catches/touchbacks given certain play attributes.
 
 The main drivers from the analysis above were:
 1. Punt Line of Scrimmage - Location on the field where the punt play started.
 2. Change of Possession Duration - The amount of time that occurs between the snap of the ball and when the returner attempts to catch the ball
 3. Receiving Team Gunner Count - The number Receiving Team players designated to defend the punt team gunners
 4. Speed - This was a function of Punt Duration and Receiving Team Gunner Count

 Below I will build predictive models and build a dataset that leverages the model information to produces some visuals
 
```{r echo=FALSE}


library("MASS")
library("scales")


model_data <- copy(final)
model_data[,penalty_major := if_else(penalty_off_holding+penalty_ill_block >0,1,0,0)]
model_data[,fair_or_touchback := if_else(fair_catch+touchback >0,1,0,0)]

lower<-quantile(model_data$change_pos_duration,probs = .05,na.rm = TRUE)
upper<-quantile(model_data$change_pos_duration,probs = .95,na.rm = TRUE)
model_data <- model_data[change_pos_duration %between% c(lower,upper),]

form <- c("yardline_num", "snapform_rec_gunner_count","snapform_rec_backfield_count","snapform_rec_line_count","change_pos_duration")
form <- paste(form,collapse="+") 

fit1 <- glm( paste("concussion~",form) , data = model_data, family = "binomial")
summary(fit1)

fit3 <- glm( paste("penalty_major~",form) , data = model_data, family = "binomial")
summary(fit3)

fit2 <- glm( paste("fair_or_touchback~",form)  ,data = model_data, family = "binomial")
summary(fit2)


#Fitting Final Models
fit1 <- glm( concussion~yardline_num+snapform_rec_gunner_count+change_pos_duration , data = model_data, family = "binomial")
fit2 <- glm(penalty_major~yardline_num+snapform_rec_gunner_count+change_pos_duration  ,data = model_data, family = "binomial")
fit3 <- glm(fair_or_touchback~yardline_num+snapform_rec_gunner_count+change_pos_duration  ,data = model_data, family = "binomial")
fit4 <- glm(log(return_yards+1)~yardline_num+snapform_rec_gunner_count+change_pos_duration+fair_catch, 
            data=model_data[])


##### Build data to get predicted probabilites
yardline_num_vec  <- as.numeric(seq(5,60,1))
snapform_rec_gunner_count_vec <- as.numeric(seq(2,5,1))
change_pos_duration_vec <- as.numeric(seq(4,9,1))
dt_pop <- CJ(yardline_num = yardline_num_vec, 
             snapform_rec_gunner_count = snapform_rec_gunner_count_vec,
             change_pos_duration = change_pos_duration_vec,
             unique = TRUE)

dt_pop[,prob_concussion := predict(fit1,newdata = .SD,type = "response")]
dt_pop[,prob_penalty_major := predict(fit2,newdata = .SD,type = "response")]
dt_pop[,prob_fair_or_touchback := predict(fit3,newdata = .SD,type = "response")]
dt_pop[,fair_catch :=1]
dt_pop[,exp_return_yds_fair_catch := exp(predict(fit4,newdata = .SD))]
dt_pop[,fair_catch :=0]
dt_pop[,exp_return_yds_nofair_catch := exp(predict(fit4,newdata = .SD))]
dt_pop[,fair_catch_diff := exp_return_yds_fair_catch - exp_return_yds_nofair_catch ]

```


```{r echo=FALSE}

# yardline_num,snapform_rec_gunner_count, prob_concussion
gdata <- dt_pop[,list(metric=mean(prob_concussion)),by=list(yardline_num,snapform_rec_gunner_count)] 
setnames(gdata,c("yardline_num","snapform_rec_gunner_count"),c("xvar","yvar"))
ggplot(data=gdata, aes(x = xvar, color= as.factor(yvar)) ) +
  ggtitle(label = "Probability of a Concussion Event" , subtitle = paste0("by Line of Scrimmage \nand # of Receiving Gunner Coverage Players") ) +
  geom_line(aes(y=metric)) +
  scale_x_continuous(name = "Line of Scrimmage : (60 = opposing to 40 yd line)" , breaks = seq(0,60,10)) +
  scale_y_continuous(name = "Probability of a Concussion Event",labels = percent) +
  labs(color="# of Receiving\n Gunner Coverage \nPlayers")


# change_pos_duration,snapform_rec_gunner_count, prob_concussion
gdata <- dt_pop[,list(metric=mean(prob_concussion)),by=list(change_pos_duration,snapform_rec_gunner_count)] 
setnames(gdata,c("change_pos_duration","snapform_rec_gunner_count"),c("xvar","yvar"))
ggplot(data=gdata, aes(x = xvar, color= as.factor(yvar)) ) +
  ggtitle(label = "Probability of a Concussion Event" , subtitle = paste0("by Change of Possession Duration \nand # of Receiving Gunner Coverage Players") ) +
  geom_line(aes(y=metric)) +
  scale_x_continuous(name = "Time : Seconds Since the Snap" ) +
  scale_y_continuous(name = "Probability of a Concussion Event",labels = percent) +
  labs(color="# of Receiving \nGunner Coverage \nPlayers")


```

Graph Insight 1:
1. More Receiving Gunner Coverage Players increases the likelihood of a concussion.
2. As the line of scrimmage moves closer to the opponent’s goal line, the likelihood of a concussion decreases
3. Increasing the Receiving Gunner Coverage Players from 2 to 4 or 5 doubles the likelihood of a concussion

Graph Insight 2:
1. More Receiving Gunner Coverage Players increases the likelihood of a concussion.
2. As play duration increases, the likelihood of a slightly increases as well



```{r echo=FALSE}

# yardline_num,snapform_rec_gunner_count, prob_penalty_major
gdata <- dt_pop[,list(metric=mean(prob_penalty_major)),by=list(yardline_num,snapform_rec_gunner_count)] 
setnames(gdata,c("yardline_num","snapform_rec_gunner_count"),c("xvar","yvar"))
ggplot(data=gdata, aes(x = xvar, color= as.factor(yvar)) ) +
  ggtitle(label = "Probability of a Penalty Event" , subtitle = paste0("by Line of Scrimmage \nand # of Receiving Gunner Coverage Players") ) +
  geom_line(aes(y=metric)) +
  scale_x_continuous(name = "Line of Scrimmage : (60 = opposing to 40 yd line)" , breaks = seq(0,60,10)) +
  scale_y_continuous(name = "Probability of a Punt Penalty Event",labels = percent) +
  labs(color="# of Receiving\n Gunner Coverage \nPlayers")

# change_pos_duration,snapform_rec_gunner_count, prob_penalty_major
gdata <- dt_pop[,list(metric=mean(prob_penalty_major)),by=list(change_pos_duration,snapform_rec_gunner_count)] 
setnames(gdata,c("change_pos_duration","snapform_rec_gunner_count"),c("xvar","yvar"))
ggplot(data=gdata, aes(x = xvar, color= as.factor(yvar)) ) +
  ggtitle(label = "Probability of a Penalty Event" , subtitle = paste0("by Change of Possession Duration \nand # of Receiving Gunner Coverage Players") ) +
  geom_line(aes(y=metric)) +
  scale_x_continuous(name = "Time : Seconds since the Snap" ) +
  scale_y_continuous(name = "Probability of Punt Penalty Events",labels = percent) +
  labs(color="# of Receiving\n Gunner Coverage \nPlayers")

```

Graph Insight 3:
1. More Receiving Gunner Coverage Players increase the likelihood of a punt penalty (Offensive Holding or any type of Illegal Block).
2. As the line of scrimmage moves closer to the opponent’s goal line, the likelihood of a punt penalty decreases
3. Increasing the Receiving Gunner Coverage Players from 2 to 4 or 5 increases the likelihood of a punt penalty by 50% to 100%

Graph Insight 4:
1. More Receiving Gunner Coverage Players increases the likelihood of a punt penalty.
2. As play duration increases, the likelihood of a punt penalty decreases



```{r echo=FALSE}

# yardline_num,snapform_rec_gunner_count, prob_penalty_major
gdata <- dt_pop[,list(metric=mean(prob_fair_or_touchback)),by=list(yardline_num,snapform_rec_gunner_count)] 
setnames(gdata,c("yardline_num","snapform_rec_gunner_count"),c("xvar","yvar"))
ggplot(data=gdata, aes(x = xvar, color= as.factor(yvar)) ) +
  ggtitle(label = "Probability of a Fair Catch/Touchback Event" , subtitle = paste0("by Line of Scrimmage \n and # of Receiving Gunner Coverage Players") ) +
  geom_line(aes(y=metric)) +
  scale_x_continuous(name = "Line of Scrimmage : (60 = opposing to 40 yd line)" , breaks = seq(0,60,10)) +
  scale_y_continuous(name = "Probability of a Fair Catch/Touchback Event",labels = percent) +
  labs(color="# of Receiving\n Gunner Coverage \nPlayers")

# change_pos_duration,snapform_rec_gunner_count, prob_penalty_major
gdata <- dt_pop[,list(metric=mean(prob_fair_or_touchback)),by=list(change_pos_duration,snapform_rec_gunner_count)] 
setnames(gdata,c("change_pos_duration","snapform_rec_gunner_count"),c("xvar","yvar"))
ggplot(data=gdata, aes(x = xvar, color= as.factor(yvar)) ) +
  ggtitle(label = "Probability of a Fair Catch/Touchback Event" , subtitle = paste0("by Change of Possession Duration \nand # of Receiving Gunner Coverage Players") ) +
  geom_line(aes(y=metric)) +
  scale_x_continuous(name = "Time : Seconds since the Snap" ) +
  scale_y_continuous(name = "Probability of a Fair Catch/Touchback Event",labels = percent) +
  labs(color="# of Receiving\n Gunner Coverage \nPlayers")

```
 

Graph Insight 5:
1. More Receiving Gunner Coverage Players decrease the likelihood of a fair catch or touchback
2. As the line of scrimmage moves closer to the opponent’s goal line, the likelihood of a fair catch or touchback increases
3. Increasing the Receiving Gunner Coverage Players from 2 to 4 or 5 decreases the likelihood of a fair catch or touchback by -50% to -100%

Graph Insight 6:
1. More Receiving Gunner Coverage Players decrease the likelihood of a fair catch or touchback
2. As play duration increases, the likelihood of a fair catch or touchback increases


## Expected Punt Return Yard Loss

As we have seen above, certain formation alignments impact concussions, penalties, and fair catch rates.
Next I will briefly look at how Loss of Field Position with Fair Catch increases.


```{r echo=FALSE}

# yardline_num,snapform_rec_gunner_count, prob_penalty_major
gdata <- dt_pop[,list(metric=mean(fair_catch_diff)),by=list(yardline_num,snapform_rec_gunner_count)] 
setnames(gdata,c("yardline_num","snapform_rec_gunner_count"),c("xvar","yvar"))
ggplot(data=gdata, aes(x = xvar, color= as.factor(yvar)) ) +
  ggtitle(label = "Loss of Field Position Due to Fair Catch" , 
          subtitle = paste0("by Line of Scrimmage \nand # of Receiving Gunner Coverage Players") ) +
  geom_line(aes(y=metric)) +
  scale_x_continuous(name = "Line of Scrimmage : (60 = opposing to 40 yd line)" , breaks = seq(0,60,10)) +
  scale_y_continuous(name = "Punt Return Yards Lost",breaks = seq(-7,0,1),limits = c(-7,0)) +
  labs(color="# of Receiving\n Gunner Coverage \nPlayers")


gdata <- dt_pop[,list(metric=mean(fair_catch_diff)),by=list(change_pos_duration,snapform_rec_gunner_count)] 
setnames(gdata,c("change_pos_duration","snapform_rec_gunner_count"),c("xvar","yvar"))
ggplot(data=gdata, aes(x = xvar, color= as.factor(yvar)) ) +
  ggtitle(label = "Loss of Field Position Due to Fair Catch" , 
          subtitle = paste0("by Change of Possession Duration \nand # of Receiving Gunner Coverage Players") ) +
  geom_line(aes(y=metric)) +
  scale_x_continuous(name = "Time : Seconds since the Snap" ) +
  scale_y_continuous(name = "Punt Return Yards Lost",breaks = seq(-7,0,1),limits = c(-7,0)) +
  labs(color="# of Receiving \nGunner Coverage \nPlayers")

```



Graph Insight 7:
1. A team will expect to lose field position due to a fair catch as the number of Receiving Gunner Coverage Players increase.  (I.e Since having more Receiving Gunner Coverage Players reduce fair catch events, they stand to lose more potential yards)
2. As the line of scrimmage moves closer to the opponent’s goal line, receiving team loss of field position is decreased.
3. Increasing the Receiving Gunner Coverage Players from 2 to 4 or 5 results in about 1 to 2 yards of field position loss.

Graph Insight 8:
1. Similar trends as in Graph 7, as play duration increases, receiving team loss of field position is decreased.



# Suggested Rule Enhancement # 1

RULE 9 SCRIMMAGE KICK
SECTION 1 KICK FROM SCRIMMAGE
ARTICLE 3.  DEFENSIVE TEAM FORMATION.

Enhancement 1 - A maximum of one player for Team A and Team B may line up between the inbounds lines and the yard-line number on each side of the field. Remaining players must line up between each yard-line numbers.


# Suggested Rule Enhancement # 2

RULE 9 SCRIMMAGE KICK 
SECTION 1 KICK FROM SCRIMMAGE
ARTICLE 3.  DEFENSIVE TEAM FORMATION. 

Enhancement 2 - Team A player lined up inside the numbers can not initiate contact with Team B player lined up between the inbounds lines and the yard-line number until after 15 yards past the original line of scrimmage.


# Potential Rule Enhancement Next Year

RULE 10 SECTION 2 FAIR CATCH ARTICLE 4 PUTTING BALL IN PLAY AFTER FAIR CATCH

Enhancement 3 - After a fair catch is made, the receiving team will be awarded an additional 5 yards from the spot of the fair catch.

# Comments on Rule Enhancements
The rule enhancement promotes 1 on 1 blocking alignment in lieu of double team blocking alignment or triple team blocking down field. This will also promote a more balanced blocking ratio with interior line players as well.  

Data analysis and visuals show sufficient evidence that reducing the gunner blocking ratio decrease the likelihood of a concussions occurring and reducing the likelihood of offensive holding and illegal blocking penalties. 

The rule enhancement will most likely increase fair catch rates and increase the speed of the gunners.  However, I believe the subtle rule enhancements will not impact overall gameplay integrity.  The above enhancements are also consistent with Kickoff Rule changes in 2018 that promoted more predictable blocking and tackling scenarios with formation alignment and designated blocking area

The rule enhancements are expected to impact receiving team by about 1 to 2 yards of field position loss.
At this point I feel awarding the receiving team an additional 5 yards to be a little pre-mature due to the expected loss of field position