--- title: 'Isolating and Analyzing Punt Returner Performance' author: 'Henry Gise (University of Pittsburgh)' date: '`r Sys.Date()`' output: html_document: number_sections: true toc: true --- # Introduction Success on special teams, especially in the punt game, is often overlooked and undervalued. Field position heavily influences drive success and is often correlated with winning and losing. Much of that influence is shouldered by the punter and the punt returner, each working to maximize its team's chances of success via optimizing field position. The punt returner only controls so much of a play, which introduces the importance of context. A punt returner has to work within his circumstances to maximize his return yardage and, thus, optimize field position. When this context is controlled for, we isolate the performance of the punt returner over the course of a return. That's what we attempt to do here. In this Notebook, a model is created that predicts return yardage at each frame of a given punt return. This will facilitate two major determinations: * The influence of individual factors on punt return success * Which returners have been the most successful ```{r setup, include = FALSE} knitr::opts_chunk$set(echo = TRUE) require(tidyverse) require(ggthemes) require(future) plan("multisession") #turning off warnings options(warn=-1) require(deldir) require(pracma) require(glmnet) require(brnn) require(randomForest) require(xgboost) require(gt) remotes::install_github("jthomasmock/gtExtras") require(gtExtras) require(nflfastR) require(ggimage) require(scales) require(gganimate) require(ggmap) require(GGally) require(magick) ``` ```{r load in the data, include = FALSE} # includes schedule info for games df_games <- read_csv("../input/nfl-big-data-bowl-2022/games.csv", col_types = cols()) # includes play-by-play info on specific plays df_plays <- read_csv("../input/nfl-big-data-bowl-2022/plays.csv", col_types = cols()) # includes background info for players df_players <- read_csv("../input/nfl-big-data-bowl-2022/players.csv", col_types = cols()) ##Reading tracking data (iterated over three years) #weeks of NFL season years <- seq(2018,2020) #blank dataframe to store tracking data df_tracking <- data.frame() #iterating through all years for (y in years){ #temporary dataframe used for reading week for given iteration df_tracking_temp <- read_csv(paste0("../input/nfl-big-data-bowl-2022/tracking",y,".csv"), col_types = cols()) #storing temporary dataframe in full dataframe df_tracking <- bind_rows(df_tracking_temp,df_tracking) } #Standardizing tracking data so its always in direction of offense vs raw on-field coordinates. df_tracking <- df_tracking %>% mutate(x = ifelse(playDirection == "left",120 - x, x), y = ifelse(playDirection == "left",160/3 - y, y), dir = ifelse(playDirection == "left", ifelse(dir - 180 < 0, dir + 180, dir - 180), dir), o = ifelse(playDirection == "left", ifelse(o - 180 < 0, o + 180, o - 180), o)) ##Reading scouting data df_scouting <- read_csv("../input/nfl-big-data-bowl-2022/PFFScoutingData.csv", col_types = cols()) ``` # Data Overview ```{r games, include = FALSE} # data from every game 2018-2020 df_games %>% dim(); df_games %>% names() ``` ```{r plays, include = FALSE} # play-level data from every game 2018-2020 df_plays %>% dim(); df_plays %>% names() ``` ```{r players, include = FALSE} # player data 2018-2020 df_players %>% dim(); df_players %>% names() ``` ```{r tracking, include = FALSE} # tracking data from every play 2018-2020 df_tracking %>% dim(); df_tracking %>% names() ``` ```{r scouting, include = FALSE} # PFF scouting data from every play 2018-2020 df_scouting %>% dim(); df_scouting %>% names() ``` All necessary infromation for each each punt return is available, including the locations and speeds of every player on the field at 0.1 second intervals. ```{r clean the data, echo = FALSE, message = FALSE} # create a dataframe from df_plays that includes only punts that included returns df_punt <- df_plays %>% filter(specialTeamsPlayType == "Punt", specialTeamsResult == "Return") %>% left_join(df_scouting, by = c("gameId","playId")) %>% select(gameId,playId,playDescription,possessionTeam,returnerId,penaltyCodes,kickReturnYardage) %>% left_join(df_games) df_punt %>% left_join(df_tracking, by = c("gameId","playId")) %>% arrange(gameId, playId, frameId) %>% select(frameId, displayName, jerseyNumber, position, x, y, s, dir, o) %>% head(6) %>% knitr::kable(format = "html", caption = "(9:20) C.Johnston punts 56 yards to ATL 36, Center-R.Lovato. J.Hardy to ATL 41...") ``` Detailed explanations of data can be found [here](https://www.kaggle.com/c/nfl-big-data-bowl-2022/data). # Variable Selection and Feature Engineering ```{r load df_return, include=FALSE} # loading data to save time df_return <- read_csv("../input/df-return/df_return.csv") ``` ```{r create df_return, eval = FALSE, include = FALSE} # creating df_return --- takes ~ 3.5 hours # establish a vector of penalties that are considered okay to occur on a play fine_penalties <- c("DOF","DSQd","FMM","HC","IDK","ILF","ILM","ISH","POK","POK;UNR","RNK","TAU","TRP","TRPd","UNR","UNRd","UNS","UNSd","UOHd",NA) # create the dataframe with all desired features df_return <- df_punt %>% # do not need to track the football filter(displayName != "football", # get rid of plays where there is a return but no returnerId !is.na(returnerId)) %>% # creating the groups for each player (returner, returnTeam, kickingTeam) mutate(group = case_when( nflId == returnerId ~ "returner", (homeTeamAbbr == possessionTeam & team == "home") | (homeTeamAbbr != possessionTeam & team == "away") ~ "kickingTeam", TRUE ~ "returnTeam" )) %>% # filter out fair catches, muffed punts, and penalties that aren't considered OK for analysis mutate(dummy = ifelse(event %in% c("fair_catch","punt_muffed") | (!penaltyCodes %in% fine_penalties),1,0)) %>% group_by(gameId,playId) %>% mutate(dummy = sum(dummy)) %>% ungroup() %>% filter(dummy == 0) %>% # filter out frames that are before punt_received mutate(dummy = ifelse(event == "punt_received" & group == "returner",frameId,0)) %>% group_by(playId,gameId) %>% mutate(dummy = sum(dummy)) %>% ungroup() %>% filter(frameId >= dummy) %>% # filter for plays in which there is a punt_received event filter(dummy != 0) %>% #filter for plays that ended in a tackle, fumble, touchdown, out of bounds, or punt downed mutate(dummy = ifelse(event %in% c("tackle","fumble","touchdown","out_of_bounds","punt_downed"),1,0), # mark the frame of the end of the play last_frame = ifelse(event %in% c("tackle","fumble","touchdown","out_of_bounds","punt_downed") & group == "returner",frameId,0), # mark plays that ended in touchdowns td_marker = ifelse(event == "touchdown" & group == "returner",1,0)) %>% group_by(gameId,playId) %>% arrange(frameId) %>% mutate(cumulative = cumsum(dummy)) %>% ungroup() %>% # filter out plays that have multiple "end" events (e.g., scoop and recovery should end on the fumble) filter(cumulative < 22 | dummy*cumulative == 22) %>% group_by(gameId,playId) %>% mutate(dummy = sum(dummy), last_frame = sum(last_frame), td_marker = sum(td_marker)) %>% ungroup() %>% # this takes out all frames that occur after the specified ending frame, ensuring the last frame is the one corresponding the ending "event" filter(dummy > 0, frameId <= last_frame) %>% select(-cumulative) %>% group_by(playId,gameId,frameId) %>% # feature engineering mutate(adj_x = x - 10, # denoting how many yards away from the target endzone adj_y = y - 80/3, # denoting how many yards from the middle of the field (positive values denote left side) adj_dir = ifelse(dir + 90 > 360, dir - 270, dir + 90), #direction player is facing w.r.t. target endzone, 0 indicates facing endzone, returner_adj_x = ifelse(group == "returner", adj_x, 0), #how many yards away returner is from target endzone returner_adj_x = sum(returner_adj_x, na.rm = TRUE), returner_adj_y = ifelse(group == "returner", adj_y, 0), #how many yards returner is from the middle of the field (positive values denote left side) returner_adj_y = sum(returner_adj_y, na.rm = TRUE), dist_from_returner = sqrt((adj_x - returner_adj_x)^2 + (adj_y - returner_adj_y)^2), #euclidean distance from returner x_s = s * cos(adj_dir * 0.0174532925199)) %>% # speed in the direction of the endzone # using kickReturnYardage and td_marker to denote the location of the end of the play group_by(gameId,playId) %>% mutate( returner_adj_x_end = ifelse(event == "punt_received" & group == "returner", returner_adj_x - kickReturnYardage, 0), returner_adj_x_end = ifelse(td_marker == 1,0,sum(returner_adj_x_end))) %>% # filtering out frames where returner is already past the goal line filter(returner_adj_x >= 0) %>% # arranging players by their distance from returner group_by(playId,gameId,frameId,group) %>% arrange(dist_from_returner) %>% mutate(proximity = ifelse(group == "returner",0,row_number())) %>% # voronoi tessellations group_by(playId,gameId,frameId) %>% mutate( # assert that as long as the runner is in bounds, we can calculate his area, otherwise it's 0 voronoi_returner_area_noBlockers = ifelse(abs(returner_adj_y) <= 80/3, v_function_area(adj_x,adj_y,group,12), 0), # area in front without blockers blockers = ifelse(abs(returner_adj_y) <= 80/3, v_function_blockers(adj_x,adj_y,group,12), 0), # number of blockers within area in front voronoi_returner_close_adj = ifelse(abs(returner_adj_y) <= 80/3, v_function_close(adj_x,adj_y,group,12), 0), # closest point on region to endzone voronoi_returner_area_in_front = ifelse(abs(returner_adj_y) <= 80/3, v_function_area(adj_x,adj_y,group,22), 0)) %>% # area in front with blockers #### filtering out those further than 2 people away from returner (on each team) filter(group %in% c("returner","kickingTeam")) %>% group_by(playId,gameId,frameId,group) %>% arrange(dist_from_returner) %>% mutate(proximity = ifelse(group == "returner",0,row_number())) %>% filter(proximity < 3) %>% ### pivot out the data so that each frame has one single observation pivot_wider(id_cols = c(week,gameId,playId,frameId,returner_adj_x_end,voronoi_returner_close_adj,voronoi_returner_area_in_front,voronoi_returner_area_noBlockers,blockers), names_from = c(group,proximity), values_from = c(adj_x, adj_y, s, x_s, dist_from_returner)) %>% group_by(gameId,playId) %>% filter(abs(adj_y_returner_0) < 80/3) %>% # getting rid of returner out of bounds plays # calculating return yards mutate(return_yards = ifelse(frameId == min(frameId), adj_x_returner_0 - returner_adj_x_end, 0)) %>% mutate(return_yards = sum(return_yards)) %>% ungroup() %>% # filter out punt that included deflection filter(gameId != 2018100711 | playId != 253, !is.na(returner_adj_x_end)) # re-join season df_return <- left_join(df_return, df_punt %>% select(gameId, playId, season)) ``` ```{r creating training and testing sets, include=FALSE} # randomly shuffle rows in df_return #### set.seed(123) rows <- sample(nrow(df_return)) # randomly remove 3,862 frames of returns > 30 yards in order to remove skewness #### df_return_jumbled <- (df_return[rows,] %>% mutate(big_return = case_when( return_yards > 20 ~ "big", return_yards <= 20 ~ "small")) %>% arrange(big_return))[c(5718:77344),] # create training set from adjusted data (3,862 frames removed) #### rows <- sample(nrow(df_return_jumbled)) df_return_initial <- df_return_jumbled[rows,] %>% # create y speed mutate(y_s_returner_0 = sqrt(s_returner_0^2 - x_s_returner_0^2)) %>% # add features to capture the rates of change in variables (e.g., rate of decrease of distance of closest defender) group_by(gameId, playId) %>% arrange(frameId) %>% mutate( # rate of change area of the area in front of the returner (w/o blockers) D_vor_area_12 = voronoi_returner_area_noBlockers - lag(voronoi_returner_area_noBlockers), D_vor_area_12 = ifelse(is.na(D_vor_area_12), lead(D_vor_area_12), D_vor_area_12), # rate of change of the distance of the closest defender D_dist_1 = dist_from_returner_kickingTeam_1 - lag(dist_from_returner_kickingTeam_1), D_dist_1 = ifelse(is.na(D_dist_1), lead(D_dist_1), D_dist_1), # rate of change of the distance of the second-closest defender D_dist_2 = dist_from_returner_kickingTeam_2 - lag(dist_from_returner_kickingTeam_2), D_dist_2 = ifelse(is.na(D_dist_2), lead(D_dist_2), D_dist_2)) %>% filter(!is.na(D_dist_2)) %>% # this removes a play that only included one frame and thus cannot be used select(season, returner_adj_x_end, # dependent variable voronoi_returner_close_adj, voronoi_returner_area_in_front, voronoi_returner_area_noBlockers, blockers, adj_x_returner_0, s_returner_0, x_s_returner_0, y_s_returner_0, dist_from_returner_kickingTeam_1, dist_from_returner_kickingTeam_2, D_vor_area_12, D_dist_1, D_dist_2, gameId, playId, week, frameId) %>% # rename some columns for ease of interpretation rename(end_of_play = returner_adj_x_end, vor_close_point = voronoi_returner_close_adj, vor_area_22 = voronoi_returner_area_in_front, vor_area_12 = voronoi_returner_area_noBlockers, adj_x_returner = adj_x_returner_0, s_returner = s_returner_0, x_s_returner = x_s_returner_0, y_s_returner = y_s_returner_0, dist_1 = dist_from_returner_kickingTeam_1, dist_2 = dist_from_returner_kickingTeam_2) # create testing set using unadjusted data (includes all frames) #### df_return_analysis <- df_return %>% # create y speed mutate(y_s_returner_0 = sqrt(s_returner_0^2 - x_s_returner_0^2)) %>% # add features to capture the rates of change in variables (e.g., rate of decrease of distance of closest defender) group_by(gameId, playId) %>% arrange(frameId) %>% mutate( # rate of change area of the area in front of the returner (w/o blockers) D_vor_area_12 = voronoi_returner_area_noBlockers - lag(voronoi_returner_area_noBlockers), D_vor_area_12 = ifelse(is.na(D_vor_area_12), lead(D_vor_area_12), D_vor_area_12), # rate of change of the distance of the closest defender D_dist_1 = dist_from_returner_kickingTeam_1 - lag(dist_from_returner_kickingTeam_1), D_dist_1 = ifelse(is.na(D_dist_1), lead(D_dist_1), D_dist_1), # rate of change of the distance of the second-closest defender D_dist_2 = dist_from_returner_kickingTeam_2 - lag(dist_from_returner_kickingTeam_2), D_dist_2 = ifelse(is.na(D_dist_2), lead(D_dist_2), D_dist_2)) %>% filter(!is.na(D_dist_2)) %>% # this removes a play that only included one frame and thus cannot be used select(season, returner_adj_x_end, # dependent variable voronoi_returner_close_adj, voronoi_returner_area_in_front, voronoi_returner_area_noBlockers, blockers, adj_x_returner_0, adj_y_returner_0, s_returner_0, x_s_returner_0, y_s_returner_0, dist_from_returner_kickingTeam_1, dist_from_returner_kickingTeam_2, D_vor_area_12, D_dist_1, D_dist_2, gameId, playId, week, frameId) %>% # rename some columns for ease of interpretation rename(end_of_play = returner_adj_x_end, vor_close_point = voronoi_returner_close_adj, vor_area_22 = voronoi_returner_area_in_front, vor_area_12 = voronoi_returner_area_noBlockers, adj_x_returner = adj_x_returner_0, adj_y_returner = adj_y_returner_0, s_returner = s_returner_0, x_s_returner = x_s_returner_0, y_s_returner = y_s_returner_0, dist_1 = dist_from_returner_kickingTeam_1, dist_2 = dist_from_returner_kickingTeam_2) ``` The model will predict the **distance from the end zone at which a play ends**, whether the end result be a tackle, fumble, out of bounds, or touchdown. ## Data Wrangling Each observation is one frame (occurring every 0.1 seconds) of one punt return. For this reason, skewness presents a minor problem. Returns that go for more yards have more frames because they generally occur over longer periods of time, which tells the model that there are more instances of long returns than there actually are. ```{r plot skewness, echo=FALSE} fills <- c("Frames" = "#FFB81C", "Plays" = "#003594") ggplot() + geom_density(data = df_return_analysis, aes(x = end_of_play, fill = "Frames"), alpha = 0.5) + geom_density(data = df_return_analysis %>% group_by(gameId, playId) %>% filter(frameId == min(frameId)), aes(x = end_of_play, fill = "Plays"), alpha = 0.5) + theme_fivethirtyeight() + labs(title = "Play-level Density and Frame-level Density", caption = "By Henry Gise | Big Data Bowl 2022", x = "End of Play Distance from Endzone (Yds)", y = "Density") + scale_fill_manual(values = fills) + theme(plot.title = element_text(size = 14, hjust = 0.5, face = "bold"), plot.caption = element_text(size = 8), axis.title = element_text(hjust = 0.5, size = 12, face = "bold"), axis.text = element_text(size = 8), legend.title = element_blank(), legend.text = element_text(size = 8), legend.key.size = unit(.5,"cm"), legend.position = "bottom") ``` We address this by doing the following. Shown below is the `Frames per Play` for returns shorter and greater than 20 yards. ```{r skewed frames per play table, echo = FALSE} df_return %>% select(gameId, playId, frameId, return_yards, returner_adj_x_end) %>% group_by(gameId, playId) %>% # group by play mutate(frames = n()) %>% ungroup() %>% # calculate the total # of frames in each play mutate(`Return Length` = case_when( return_yards > 20 ~ "> 20 yards", # longer returns return_yards <= 20 ~ "<= 20 yards")) %>% #shorter returns select(-frameId) %>% distinct() %>% group_by(`Return Length`) %>% summarise("Total Plays" = n(), "Total Frames" = sum(frames)) %>% mutate("Frames per Play" = `Total Frames`/`Total Plays`) %>% # total frames per play for long vs. short returns knitr::kable(format = "html") ``` By randomly removing enough long-return frames from `df_return` to ensure that `Frames per Play` is equal among shorter and longer returns, some bias is removed from the model. In this case **5,717** rows are removed from the training data, cutting the `Total Frames` for these long returns down to **5,606**. ```{r short returns} # Frames per play for returns <= 20 yards 66021 / 1790 ``` ```{r long returns} # Frames per play for returns > 20 yards (11323 - 5717) / 152 ``` Notably, the data used for testing is not modified, just the training data. ## Feature Engineering The features for this model include the following, with examples from the first frame of Dwayne Harris's 99-yard TD: ```{r introduction to features, echo = FALSE} tibble( "Variable" = names(df_return_initial)[c(2,7:12,14:15)], "Definition" = c("OUTCOME - The point at which the play ends", "Returner's current distance from endzone", "Overall speed of the returner", "Speed of the returner in the direction of the endzone", "Horizontal speed of the returner (in either direction)", "Distance between returner and closest defender", "Distance between returner and second-closest defender", "Rate of change of dist_1", "Rate of change of dist_2"), "Example" = as.vector(t((df_return_analysis %>% filter(gameId == 2018122400, playId == 241, frameId == min(frameId)))[c(2,7,9:13,15:16)])), "Units" = c("yds","yds","yds/second","yds/second","yds/second","yds","yds","yds/decisecond","yds/decisecond") ) %>% knitr::kable(format = 'html') ``` Beyond that, **Voronoi tessellations** ([Yurko, et al.](https://www.degruyter.com/document/doi/10.1515/jqas-2019-0056/html)) were used to extract specific features relating to returner spacing on the field. A Voronoi tessellation partitions a plane (a football field) into regions representing the set of points on a plane which are closer to one specified object (player) than any other. Thus, a single player's Voronoi region is the set of all points on the football field that are closer to him than to anyone else. Two separate Voronoi tessellations were created from each frame of each play --- one including all 22 players, and one excluding the returner's teammates. From these tessellations, the following features are extracted for each frame in each return. ```{r explantion of Voronoi, echo = FALSE} tibble( "Variable" = names(df_return_initial)[c(4:5,3,6,13)], "Definition" = c("Area in front of returner within his Voronoi region (all 22 players)", "Area in front of returner within his Voronoi region (without teammates)", "Distance of closest point within vor_area_12 to the endzone", "The number of teammates located within vor_area_12", "Rate of change of vor_area_12"), "Example" = as.vector(t((df_return_analysis %>% filter(gameId == 2018122400, playId == 241, frameId == min(frameId)))[c(4:5,3,6,14)])), "Units" = c("sq. yds","sq. yds","yds","players","sq. yds/decisecond") ) %>% knitr::kable(format = "html") ``` ```{r create voronoi functions, eval = FALSE, include = FALSE} field <- c(0,110,-80/3,80/3) v_function_area <- function(x,y,group,blocking){ if (blocking != 22){ players <- data.frame(cbind(x,y,group)) %>% filter(group %in% c("returner","kickingTeam")) x <- as.numeric(players$x) y <- as.numeric(players$y) } delTest <- deldir(x,y,rw = field) # create tessellations df_area <- data.frame( "x_area" = c((delTest$dirsgs %>% filter(ind1 == 1 | ind2 == 1))$x1, (delTest$dirsgs %>% filter(ind1 == 1 | ind2 == 1))$x2), "y_area" = c((delTest$dirsgs %>% filter(ind1 == 1 | ind2 == 1))$y1, (delTest$dirsgs %>% filter(ind1 == 1 | ind2 == 1))$y2)) %>% distinct() %>% mutate(sideline_checker = ifelse(x_area < x[1], 1, 0)) # create initial dataframe of returner's region's vertices if (prod(df_area$sideline_checker) == 1){ df_area <- df_area %>% rbind(data.frame( "x_area" = c(x[1], x[1]), "y_area" = c(-80/3, 80/3), "sideline_checker" = 0)) } # in case returner has two open sidelines next to him, add points on the sidelines x_center <- mean(df_area$x_area) y_center <- mean(df_area$y_area) df_area <- df_area %>% select(-sideline_checker) %>% mutate(angle = atan2(y_area - mean(y_area), x_area - mean(x_area))*180/3.14159) # add angle for sorting points if (NROW(df_area %>% filter(y_area <= -26.66666)) == 1){ # the case where only one point touches the lower y axis if (NROW(df_area %>% filter(x_area == 0)) == 1){ # the case where a point touches opponent's endzone df_area <- df_area %>% rbind(data.frame( "x_area" = 0, "y_area" = -80/3) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159))} if (NROW(df_area %>% filter(x_area == 110)) == 1){ # the case where a point touches own endzone df_area <- df_area %>% rbind(data.frame( "x_area" = 110, "y_area" = -80/3) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159))} } # adding points on the front or back pylons if necessary (- sideline) if (NROW(df_area %>% filter(y_area >= 26.66666)) == 1){ # the case where only one point touches the upper y axis if (NROW(df_area %>% filter(x_area == 0)) == 1){ # the case where a point touches opponent's endzone df_area <- df_area %>% rbind(data.frame( "x_area" = 0, "y_area" = 80/3) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159)) } if (NROW(df_area %>% filter(x_area == 110)) == 1){ # the case where a point touches own endzone df_area <- df_area %>% rbind(data.frame( "x_area" = 110, "y_area" = 80/3) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159)) } } # adding points on the front or back pylons if necessary (+ sideline) if (NROW(df_area %>% filter(y_area <= -26.66666)) == 1 & NROW(df_area %>% filter(y_area >= 26.66666)) == 1){ if (x[1] == min(x)){ df_area <- df_area %>% rbind(data.frame( "x_area" = c(0,0), "y_area" = c(-80/3,80/3)) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159))} # if the returner is in front of the entire defense if (NROW(df_area %>% filter(abs(y_area) >= 26.66666, x_area < x[1])) == 2){ if (x[1] == max(x)){ df_area <- df_area %>% rbind(data.frame( "x_area" = c(x[1],x[1]), "y_area" = c(-80/3,80/3)) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159))} # if the returner is behind the entire defense } # if both sideline points are in front of the returner if (NROW(df_area %>% filter(abs(y_area) >= 26.66666, x_area < x[1])) == 1){ temp_x <- (df_area %>% filter(abs(y_area) >= 26.66666))$x_area temp_y <- (df_area %>% filter(abs(y_area) >= 26.66666))$y_area if (temp_x[1] < x[1]){ df_area <- df_area %>% rbind(data.frame( "x_area" = x[1], "y_area" = temp_y[1]) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159))} # if the first point is in front of returner if (temp_x[2] < x[1]){ df_area <- df_area %>% rbind(data.frame( "x_area" = x[1], "y_area" = temp_y[2]) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159))} # if the second point is in front of returner } # if only one sideline point is in front of the returner } # if there are two sideline points df_area <- df_area %>% arrange(angle) # sort points df_area[NROW(df_area) + 1,] <- df_area[1,] # add extra point to allow for slope calculations df_area <- df_area %>% mutate(m = (lead(y_area) - y_area) / (lead(x_area) - x_area)) # add slope for calculating intersections crossers <- df_area %>% filter((x_area < x[1] & lead(x_area) > x[1]) | (x_area > x[1] & lead(x_area) < x[1])) %>% mutate(intersection = m*(x[1] - x_area) + y_area, intersection = ifelse(intersection > 80/3, 80/3, ifelse(intersection < -80/3, -80/3, intersection))) # calculating intersections min_int <- -Inf max_int <- Inf if (NROW(crossers) > 0){ min_int <- ifelse(min(crossers$intersection) > y[1], -Inf, min(crossers$intersection)) max_int <- ifelse(max(crossers$intersection) < y[1], Inf, max(crossers$intersection))} # find the min and max y points on region final_area <- df_area %>% select(x_area, y_area) %>% rbind(data.frame( "x_area" = rep(x[1], length(crossers$intersection)), "y_area" = crossers$intersection)) %>% mutate(x_area = ifelse(x_area > x[1], x[1], x_area)) %>% filter(x_area < x[1] | (y_area <= max_int & y_area >= min_int)) %>% distinct() %>% mutate(angle = atan2(y_area - mean(y_area), x_area - mean(x_area))*180/3.14159) %>% arrange(angle) area_in_front <- polyarea(final_area$x_area, final_area$y_area) # area in front of returner } # this function returns the area in front of the returner (all 22) or (no blockers) v_function_close <- function(x,y,group,blocking){ if (blocking != 22){ players <- data.frame(cbind(x,y,group)) %>% filter(group %in% c("returner","kickingTeam")) x <- as.numeric(players$x) y <- as.numeric(players$y) } delTest <- deldir(x,y,rw = field) # create tessellations df_area <- data.frame( "x_area" = c((delTest$dirsgs %>% filter(ind1 == 1 | ind2 == 1))$x1, (delTest$dirsgs %>% filter(ind1 == 1 | ind2 == 1))$x2), "y_area" = c((delTest$dirsgs %>% filter(ind1 == 1 | ind2 == 1))$y1, (delTest$dirsgs %>% filter(ind1 == 1 | ind2 == 1))$y2)) %>% distinct() %>% mutate(sideline_checker = ifelse(x_area < x[1], 1, 0)) # create initial dataframe of returner's region's vertices if (prod(df_area$sideline_checker) == 1){ df_area <- df_area %>% rbind(data.frame( "x_area" = c(x[1], x[1]), "y_area" = c(-80/3, 80/3), "sideline_checker" = 0)) } # in case returner has two open sidelines next to him, add points on the sidelines x_center <- mean(df_area$x_area) y_center <- mean(df_area$y_area) df_area <- df_area %>% select(-sideline_checker) %>% mutate(angle = atan2(y_area - mean(y_area), x_area - mean(x_area))*180/3.14159) # add angle for sorting points if (NROW(df_area %>% filter(y_area <= -26.66666)) == 1){ # the case where only one point touches the lower y axis if (NROW(df_area %>% filter(x_area == 0)) == 1){ # the case where a point touches opponent's endzone df_area <- df_area %>% rbind(data.frame( "x_area" = 0, "y_area" = -80/3) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159))} if (NROW(df_area %>% filter(x_area == 110)) == 1){ # the case where a point touches own endzone df_area <- df_area %>% rbind(data.frame( "x_area" = 110, "y_area" = -80/3) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159))} } # adding points on the front or back pylons if necessary (- sideline) if (NROW(df_area %>% filter(y_area >= 26.66666)) == 1){ # the case where only one point touches the upper y axis if (NROW(df_area %>% filter(x_area == 0)) == 1){ # the case where a point touches opponent's endzone df_area <- df_area %>% rbind(data.frame( "x_area" = 0, "y_area" = 80/3) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159)) } if (NROW(df_area %>% filter(x_area == 110)) == 1){ # the case where a point touches own endzone df_area <- df_area %>% rbind(data.frame( "x_area" = 110, "y_area" = 80/3) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159)) } } # adding points on the front or back pylons if necessary (+ sideline) if (NROW(df_area %>% filter(y_area <= -26.66666)) == 1 & NROW(df_area %>% filter(y_area >= 26.66666)) == 1){ if (x[1] == min(x)){ df_area <- df_area %>% rbind(data.frame( "x_area" = c(0,0), "y_area" = c(-80/3,80/3)) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159))} # if the returner is in front of the entire defense if (NROW(df_area %>% filter(abs(y_area) >= 26.66666, x_area < x[1])) == 2){ if (x[1] == max(x)){ df_area <- df_area %>% rbind(data.frame( "x_area" = c(x[1],x[1]), "y_area" = c(-80/3,80/3)) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159))} # if the returner is behind the entire defense } # if both sideline points are in front of the returner if (NROW(df_area %>% filter(abs(y_area) >= 26.66666, x_area < x[1])) == 1){ temp_x <- (df_area %>% filter(abs(y_area) >= 26.66666))$x_area temp_y <- (df_area %>% filter(abs(y_area) >= 26.66666))$y_area if (temp_x[1] < x[1]){ df_area <- df_area %>% rbind(data.frame( "x_area" = x[1], "y_area" = temp_y[1]) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159))} # if the first point is in front of returner if (temp_x[2] < x[1]){ df_area <- df_area %>% rbind(data.frame( "x_area" = x[1], "y_area" = temp_y[2]) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159))} # if the second point is in front of returner } # if only one sideline point is in front of the returner } # if there are two sideline points df_area <- df_area %>% arrange(angle) # sort points df_area[NROW(df_area) + 1,] <- df_area[1,] # add extra point to allow for slope calculations df_area <- df_area %>% mutate(m = (lead(y_area) - y_area) / (lead(x_area) - x_area)) # add slope for calculating intersections crossers <- df_area %>% filter((x_area < x[1] & lead(x_area) > x[1]) | (x_area > x[1] & lead(x_area) < x[1])) %>% mutate(intersection = m*(x[1] - x_area) + y_area, intersection = ifelse(intersection > 80/3, 80/3, ifelse(intersection < -80/3, -80/3, intersection))) # calculating intersections min_int <- -Inf max_int <- Inf if (NROW(crossers) > 0){ min_int <- ifelse(min(crossers$intersection) > y[1], -Inf, min(crossers$intersection)) max_int <- ifelse(max(crossers$intersection) < y[1], Inf, max(crossers$intersection))} # find the min and max y points on region final_area <- df_area %>% select(x_area, y_area) %>% rbind(data.frame( "x_area" = rep(x[1], length(crossers$intersection)), "y_area" = crossers$intersection)) %>% mutate(x_area = ifelse(x_area > x[1], x[1], x_area)) %>% filter(x_area < x[1] | (y_area <= max_int & y_area >= min_int)) %>% distinct() %>% mutate(angle = atan2(y_area - mean(y_area), x_area - mean(x_area))*180/3.14159) %>% arrange(angle) close <- min(final_area$x_area) } # this function returns the closest point from target endzone on returner's tessellation (no blockers) v_function_blockers <- function(x_all,y_all,group,blocking){ if (blocking != 22){ players <- data.frame(cbind(x_all,y_all,group)) %>% filter(group %in% c("returner","kickingTeam")) x <- as.numeric(players$x) y <- as.numeric(players$y) } delTest <- deldir(x,y,rw = field) # create tessellations df_area <- data.frame( "x_area" = c((delTest$dirsgs %>% filter(ind1 == 1 | ind2 == 1))$x1, (delTest$dirsgs %>% filter(ind1 == 1 | ind2 == 1))$x2), "y_area" = c((delTest$dirsgs %>% filter(ind1 == 1 | ind2 == 1))$y1, (delTest$dirsgs %>% filter(ind1 == 1 | ind2 == 1))$y2)) %>% distinct() %>% mutate(sideline_checker = ifelse(x_area < x[1], 1, 0)) # create initial dataframe of returner's region's vertices if (prod(df_area$sideline_checker) == 1){ df_area <- df_area %>% rbind(data.frame( "x_area" = c(x[1], x[1]), "y_area" = c(-80/3, 80/3), "sideline_checker" = 0)) } # in case returner has two open sidelines next to him, add points on the sidelines x_center <- mean(df_area$x_area) y_center <- mean(df_area$y_area) df_area <- df_area %>% select(-sideline_checker) %>% mutate(angle = atan2(y_area - mean(y_area), x_area - mean(x_area))*180/3.14159) # add angle for sorting points if (NROW(df_area %>% filter(y_area <= -26.66666)) == 1){ # the case where only one point touches the lower y axis if (NROW(df_area %>% filter(x_area == 0)) == 1){ # the case where a point touches opponent's endzone df_area <- df_area %>% rbind(data.frame( "x_area" = 0, "y_area" = -80/3) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159))} if (NROW(df_area %>% filter(x_area == 110)) == 1){ # the case where a point touches own endzone df_area <- df_area %>% rbind(data.frame( "x_area" = 110, "y_area" = -80/3) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159))} } # adding points on the front or back pylons if necessary (- sideline) if (NROW(df_area %>% filter(y_area >= 26.66666)) == 1){ # the case where only one point touches the upper y axis if (NROW(df_area %>% filter(x_area == 0)) == 1){ # the case where a point touches opponent's endzone df_area <- df_area %>% rbind(data.frame( "x_area" = 0, "y_area" = 80/3) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159)) } if (NROW(df_area %>% filter(x_area == 110)) == 1){ # the case where a point touches own endzone df_area <- df_area %>% rbind(data.frame( "x_area" = 110, "y_area" = 80/3) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159)) } } # adding points on the front or back pylons if necessary (+ sideline) if (NROW(df_area %>% filter(y_area <= -26.66666)) == 1 & NROW(df_area %>% filter(y_area >= 26.66666)) == 1){ if (x[1] == min(x)){ df_area <- df_area %>% rbind(data.frame( "x_area" = c(0,0), "y_area" = c(-80/3,80/3)) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159))} # if the returner is in front of the entire defense if (NROW(df_area %>% filter(abs(y_area) >= 26.66666, x_area < x[1])) == 2){ if (x[1] == max(x)){ df_area <- df_area %>% rbind(data.frame( "x_area" = c(x[1],x[1]), "y_area" = c(-80/3,80/3)) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159))} # if the returner is behind the entire defense } # if both sideline points are in front of the returner if (NROW(df_area %>% filter(abs(y_area) >= 26.66666, x_area < x[1])) == 1){ temp_x <- (df_area %>% filter(abs(y_area) >= 26.66666))$x_area temp_y <- (df_area %>% filter(abs(y_area) >= 26.66666))$y_area if (temp_x[1] < x[1]){ df_area <- df_area %>% rbind(data.frame( "x_area" = x[1], "y_area" = temp_y[1]) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159))} # if the first point is in front of returner if (temp_x[2] < x[1]){ df_area <- df_area %>% rbind(data.frame( "x_area" = x[1], "y_area" = temp_y[2]) %>% mutate(angle = atan2(y_area - y_center, x_area - x_center)*180/3.14159))} # if the second point is in front of returner } # if only one sideline point is in front of the returner } # if there are two sideline points df_area <- df_area %>% arrange(angle) # sort points df_area[NROW(df_area) + 1,] <- df_area[1,] # add extra point to allow for slope calculations df_area <- df_area %>% mutate(m = (lead(y_area) - y_area) / (lead(x_area) - x_area)) # add slope for calculating intersections crossers <- df_area %>% filter((x_area < x[1] & lead(x_area) > x[1]) | (x_area > x[1] & lead(x_area) < x[1])) %>% mutate(intersection = m*(x[1] - x_area) + y_area, intersection = ifelse(intersection > 80/3, 80/3, ifelse(intersection < -80/3, -80/3, intersection))) # calculating intersections min_int <- -Inf max_int <- Inf if (NROW(crossers) > 0){ min_int <- ifelse(min(crossers$intersection) > y[1], -Inf, min(crossers$intersection)) max_int <- ifelse(max(crossers$intersection) < y[1], Inf, max(crossers$intersection))} # find the min and max y points on region final_area <- df_area %>% select(x_area, y_area) %>% rbind(data.frame( "x_area" = rep(x[1], length(crossers$intersection)), "y_area" = crossers$intersection)) %>% mutate(x_area = ifelse(x_area > x[1], x[1], x_area)) %>% filter(x_area < x[1] | (y_area <= max_int & y_area >= min_int)) %>% distinct() %>% mutate(angle = atan2(y_area - mean(y_area), x_area - mean(x_area))*180/3.14159) %>% arrange(angle) blockers <- sum((data.frame(cbind(x_all,y_all,group)) %>% filter(group == "returnTeam") %>% mutate(within = ifelse( as.numeric(x_all) >= min(final_area$x_area) & as.numeric(x_all) <= max(final_area$x_area) & as.numeric(y_all) >= min(final_area$y_area) & as.numeric(y_all) <= max(final_area$y_area), 1, 0)))$within) } # this function returns the total no. of blockers inside returner's in-front region ``` These extracted features are illustrated in context below. ``` {r print images, echo = FALSE} magick::image_read('../input/d-harris-22-annot/d_harris_22_annot.png') magick::image_read('../input/d-harris-12-annot/d_harris_12_annot.png') ``` After feature engineering, we have the following data: ```{r glimpse at analysis, echo = FALSE} glimpse(df_return_analysis) ``` # Modeling Six different models were trained using leave-one-season-out (LOSO) cross-validation and root-mean-square-error (RMSE) as the error metric. ```{r OLS, include = FALSE} MSE_set_lm <- data.frame() results_lm <- data.frame() # LOSO cross-validation for (year in 2018:2020){ # create training set #### train <- df_return_initial %>% filter(season != year) # create testing set #### test <- df_return_analysis %>% filter(season == year) # model #### lm_fit <- lm(end_of_play ~ . - playId, data = train %>% ungroup() %>% select(-week,-gameId,-frameId)) # predict #### y_test_results_lm <- as.data.frame(predict.glm(lm_fit, test)) y_test_results_lm$`predict.glm(lm_fit, test)`[y_test_results_lm$`predict.glm(lm_fit, test)` < 0] = 0 # post process #### # this creates a dataset with the predicted values, actual values, and square errors subtract_lm <- cbind(y_test_results_lm,test$end_of_play) %>% rename(predicted = `predict.glm(lm_fit, test)`, actual = `test$end_of_play`) %>% mutate(SE = (predicted - actual)^2, over_expected = predicted - actual) # this combined all of the datasets for all seasons results_lm <- rbind(results_lm, subtract_lm %>% cbind(test)) %>% group_by(gameId, playId) %>% arrange(frameId) %>% mutate(pred_change = predicted - lag(predicted), pred_change = ifelse(is.na(pred_change), 0, pred_change), avg_change = sum(abs(pred_change))/(n() - 1)) # calculate MSE MSE = mean(subtract_lm$SE) # create a list of all MSE's for each season MSE_set_lm <- rbind(MSE_set_lm,MSE) # display RMSE for each season print(paste0(year, " RMSE: ", sqrt(MSE))) } print(paste0("RMSE: ", sqrt(mean(MSE_set_lm[,1])))) ``` ```{r Ridge, include=FALSE} MSE_set_ridge <- data.frame() results_ridge <- data.frame() lambda_seq <- 10^seq(2, -2, by = -.1) # set range of lambdas for (year in 2018:2020){ # create training set #### train <- df_return_initial %>% filter(season != year) # create testing set #### test <- df_return_analysis %>% filter(season == year) # create matrices #### x_train <- data.matrix(train[,3:15]) # make x values x_test <- data.matrix(test[,c(3:7,9:16)]) # make x values y_train <- data.matrix(train[,c("end_of_play")]) # make y values y_test <- data.matrix(test[,c("end_of_play")]) # make y values # model #### ridge_cv <- cv.glmnet(x_train, y_train, alpha = 0, lambda = lambda_seq) # cross validation best_lambda <- ridge_cv$lambda.min # get lambda value that offers minimum best_ridge <- glmnet(x_train, y_train, alpha = 0, lambda = best_lambda) # build final model # predict #### y_test_results <- as.data.frame(predict.glmnet(best_ridge, x_test, s = best_lambda)) y_test_results$`1`[y_test_results$`1` < 0] = 0 # this creates a dataset with the predicted values, actual values, and square errors subtract_ridge <- cbind(y_test_results,test$end_of_play) %>% rename(predicted = `1`, actual = `test$end_of_play`) %>% mutate(SE = (predicted - actual)^2, over_expected = predicted - actual) # this combined all of the datasets for all seasons results_ridge <- rbind(results_ridge, subtract_ridge %>% cbind(test)) %>% group_by(gameId, playId) %>% arrange(frameId) %>% mutate(pred_change = predicted - lag(predicted), pred_change = ifelse(is.na(pred_change), 0, pred_change), avg_change = sum(abs(pred_change))/(n() - 1)) # calculate MSE MSE = mean(subtract_ridge$SE) # create a list of all MSE's for each season MSE_set_ridge <- rbind(MSE_set_ridge,MSE) # display RMSE and value of lambda for each season print(paste0(year, " RMSE: ", sqrt(MSE), " lambda: ", best_lambda)) } print(paste0("RMSE: ", sqrt(mean(MSE_set_ridge[,1])))) ``` ```{r LASSO, include=FALSE} MSE_set_lasso <- data.frame() results_lasso <- data.frame() lambda_seq <- 10^seq(2, -2, by = -.1) # set range of lambdas for (year in 2018:2020){ # create training set #### train <- df_return_initial %>% filter(season != year) # create testing set #### test <- df_return_analysis %>% filter(season == year) # create matrices #### x_train <- data.matrix(train[,3:15]) # make x values x_test <- data.matrix(test[,c(3:7,9:16)]) # make x values y_train <- data.matrix(train[,c("end_of_play")]) # make y values y_test <- data.matrix(test[,c("end_of_play")]) # make y values # model #### lasso_cv <- cv.glmnet(x_train, y_train, alpha = 1, lambda = lambda_seq) # cross validation best_lambda <- lasso_cv$lambda.min # get lambda value that offers minimum best_lasso <- glmnet(x_train, y_train, alpha = 1, lambda = best_lambda) # build final model # predict #### y_test_results <- as.data.frame(predict.glmnet(best_lasso, x_test, s = best_lambda)) y_test_results$`1`[y_test_results$`1` < 0] = 0 # this creates a dataset with the predicted values, actual values, and square errors subtract_lasso <- cbind(y_test_results,test$end_of_play) %>% rename(predicted = `1`, actual = `test$end_of_play`) %>% mutate(SE = (predicted - actual)^2, over_expected = predicted - actual) # this combined all of the datasets for all seasons results_lasso <- rbind(results_lasso, subtract_lasso %>% cbind(test)) %>% group_by(gameId, playId) %>% arrange(frameId) %>% mutate(pred_change = predicted - lag(predicted), pred_change = ifelse(is.na(pred_change), 0, pred_change), avg_change = sum(abs(pred_change))/(n() - 1)) # calculate MSE MSE = mean(subtract_lasso$SE) # create a list of all MSE's for each season MSE_set_lasso <- rbind(MSE_set_lasso,MSE) # display RMSE and value of lambda for each season print(paste0(year, " RMSE: ", sqrt(MSE), " lambda: ", best_lambda)) } print(paste0("RMSE: ", sqrt(mean(MSE_set_lasso[,1])))) ``` ```{r Bayesian Regularized Neural Network (brnn), include=FALSE} install.packages("brnn") MSE_set_nn <- data.frame() results_nn <- data.frame() # LOSO cross-validation for (year in 2018:2020){ # create training set #### train <- df_return_initial %>% filter(season != year) # create testing set #### test <- df_return_analysis %>% filter(season == year) # model #### best_nn <- brnn::brnn(end_of_play ~ . - playId, data = train %>% ungroup() %>% select(-week,-gameId,-frameId), normalize = TRUE) # predict #### y_test_results_nn <- as.data.frame(brnn::predict.brnn(best_nn, test)) y_test_results_nn$`brnn::predict.brnn(best_nn, test)`[y_test_results_nn$`brnn::predict.brnn(best_nn, test)` < 0] = 0 # post process #### # this creates a dataset with the predicted values, actual values, and square errors subtract_nn <- cbind(y_test_results_nn,test$end_of_play) %>% rename(predicted = `brnn::predict.brnn(best_nn, test)`, actual = `test$end_of_play`) %>% mutate(SE = (predicted - actual)^2, over_expected = predicted - actual) # this combined all of the datasets for all seasons results_nn <- rbind(results_nn, subtract_nn %>% cbind(test)) %>% group_by(gameId, playId) %>% arrange(frameId) %>% mutate(pred_change = predicted - lag(predicted), pred_change = ifelse(is.na(pred_change), 0, pred_change), avg_change = sum(abs(pred_change))/(n() - 1)) # calculate MSE MSE = mean(subtract_nn$SE) # create a list of all MSE's for each season MSE_set_nn <- rbind(MSE_set_nn,MSE) # display RMSE for each season print(paste0(year, " RMSE: ", sqrt(MSE))) } print(paste0("RMSE: ", sqrt(mean(MSE_set_nn[,1])))) ``` ```{r xgBoost, include=FALSE} MSE_set_boost <- data.frame() results_boost <- data.frame() for (year in 2018:2020){ # create training set #### train <- df_return_initial %>% filter(season != year) # create testing set #### test <- df_return_analysis %>% filter(season == year) # create matrices #### x_train <- data.matrix(train[,3:15]) # make x values x_test <- data.matrix(test[,c(3:7,9:16)]) # make x values y_train <- data.matrix(train[,c("end_of_play")]) # make y values y_test <- data.matrix(test[,c("end_of_play")]) # make y values # model #### best_boost <- xgboost(data = x_train, label = y_train, nrounds = 19, objective = "reg:squarederror", verbose = 0) # predict #### y_test_results_boost <- as.data.frame(predict(best_boost, x_test)) y_test_results_boost$`predict(best_boost, x_test)`[y_test_results_boost$`predict(best_boost, x_test)` < 0] = 0 # this creates a dataset with the predicted values, actual values, and square errors subtract_boost <- cbind(y_test_results_boost,test$end_of_play) %>% rename(predicted = `predict(best_boost, x_test)`, actual = `test$end_of_play`) %>% mutate(SE = (predicted - actual)^2, over_expected = predicted - actual) # this combined all of the datasets for all seasons results_boost <- rbind(results_boost, subtract_boost %>% cbind(test)) %>% group_by(gameId, playId) %>% arrange(frameId) %>% mutate(pred_change = predicted - lag(predicted), pred_change = ifelse(is.na(pred_change), 0, pred_change), avg_change = sum(abs(pred_change))/(n() - 1)) # calculate MSE MSE = mean(subtract_boost$SE) # create a list of all MSE's for each season MSE_set_boost <- rbind(MSE_set_boost,MSE) # display RMSE and value of lambda for each season print(paste0(year, " RMSE: ", sqrt(MSE))) } print(paste0("RMSE: ", sqrt(mean(MSE_set_boost[,1])))) ``` ```{r Random Forest, include = FALSE} MSE_set_rf <- data.frame() results_rf <- data.frame() for (year in 2018:2020){ # create training set #### train <- df_return_initial %>% filter(season != year) # create testing set #### test <- df_return_analysis %>% filter(season == year) # create matrices #### x_train <- data.matrix(train[,3:15]) # make x values y_train <- data.matrix(train[,c("end_of_play")]) # make y values # model #### best_rf <- randomForest(x_train, y_train, ntree = 64) # build model # predict #### y_test_results_rf <- as.data.frame(predict(best_rf, test)) y_test_results_rf$`predict(best_rf, test)`[y_test_results_rf$`predict(best_rf, test)` < 0] = 0 # this creates a dataset with the predicted values, actual values, and square errors subtract_rf <- cbind(y_test_results_rf,test$end_of_play) %>% rename(predicted = `predict(best_rf, test)`, actual = `test$end_of_play`) %>% mutate(SE = (predicted - actual)^2, over_expected = predicted - actual) # this combined all of the datasets for all seasons results_rf <- rbind(results_rf, subtract_rf %>% cbind(test)) %>% group_by(gameId, playId) %>% arrange(frameId) %>% mutate(pred_change = predicted - lag(predicted), pred_change = ifelse(is.na(pred_change), 0, pred_change), avg_change = sum(abs(pred_change))/(n() - 1)) # calculate MSE MSE = mean(subtract_rf$SE) # create a list of all MSE's for each season MSE_set_rf <- rbind(MSE_set_rf,MSE) # display RMSE and value of lambda for each season print(paste0(year, " RMSE: ", sqrt(MSE))) } print(paste0("RMSE: ", sqrt(mean(MSE_set_rf[,1])))) ``` Performance is summarized below. ```{r RMSE table, echo = FALSE} tibble( "Model" = c("OLS", "Ridge", "Lasso", "BRNN", "xgBoost", "RF"), "RMSE" = c(sqrt(mean(MSE_set_lm[,1])), sqrt(mean(MSE_set_ridge[,1])), sqrt(mean(MSE_set_lasso[,1])), sqrt(mean(MSE_set_nn[,1])), sqrt(mean(MSE_set_boost[,1])), sqrt(mean(MSE_set_rf[,1])))) %>% knitr::kable(format = "html") ``` # Analysis From here, the **Random Forest** model is selected, given it produces the lowest RMSE on the test data. Feature importance can then be visualized, illustrating the degree to which different within-play factors contribute to return yardage. ```{r feature importance, echo = FALSE} ggplot(data.frame(best_rf$importance) %>% rownames_to_column(var = "Feature") %>% rename(Importance = IncNodePurity)) + geom_col(aes(x = reorder(Feature, Importance), y = Importance), fill = "lightblue", color = "navyblue") + theme_fivethirtyeight() + labs(title = "Feature Importance From Random Forest Model", caption = "By Henry Gise | Big Data Bowl 2022", x = "Feature", y = "Importance (RSS)") + theme_fivethirtyeight() + coord_flip() + theme(plot.title = element_text(size = 14, hjust = 0.5, face = "bold"), plot.caption = element_text(size = 8), axis.title = element_text(hjust = 0.5, size = 12, face = "bold"), axis.text = element_text(size = 8), legend.title = element_blank(), legend.text = element_text(size = 8), legend.key.size = unit(.5,"cm"), legend.position = "bottom") ``` It's clear that the returner's position (`adj_x_returner`) and his `vor_close_point`, which attempts to measure the space in front of him, make up the bulk of the importance. The significance of the `vor_close_point` feature is apparent, and a big reason is its ability to mark when a returner "breaks free" during a return. When the closest point on his Voronoi region is the endzone, there is effectively nobody on the coverage team closer to the endzone than he (at least for some stretch of the goal line). This effect can be illustrated using Dwayne Harris's 99-yard return TD. ```{r Dwayne Harris map, echo = FALSE} colors <- c("vor_close_point" = "darkgreen", "adj_x_returner" = "blue", "Expected Play End" = "red", "Actual Play End" = "black") ggplot(results_rf %>% arrange(frameId) %>% filter(gameId == 2018122400, playId == 241), aes(x = frameId - min(frameId) + 1)) + geom_path(aes(y = actual, color = "Actual Play End"), size = 3, alpha = .5) + geom_path(aes(y = vor_close_point, color = "vor_close_point"), size = 3, alpha = .5) + geom_path(aes(y = adj_x_returner, color = "adj_x_returner"), size = 3, alpha = .5) + geom_path(aes(y = predicted, color = "Expected Play End"), size = 3, alpha = .5) + scale_color_manual(values = colors) + labs(title = "Return Yardage Prediction Through Dwayne Harris's 99-yd TD", y = "Distance from Endzone (yds)", x = "Frames into Return", color = "Legend") + theme(plot.title = element_text(size = 14, hjust = 0.5, face = "bold"), plot.caption = element_text(size = 8), axis.title = element_text(hjust = 0.5, size = 12, face = "bold"), axis.text = element_text(size = 8), legend.title = element_blank(), legend.text = element_text(size = 8), legend.key.size = unit(.5,"cm"), legend.position = "bottom") ``` We can clearly see the point at which Harris broke away from the coverage team and was off to the races. Further, we can utilize this model for its most widely accepted purpose, which is to compute individual **return yards over expected**. This is done by determining an expectation during the first frame of the play and taking the difference from the result. ```{r return yards over expected, include = FALSE} df_ryoe_player <- results_rf %>% group_by(gameId, playId) %>% filter(season == 2020, frameId == min(frameId)) %>% left_join(df_plays) %>% left_join(df_players %>% mutate(nflId = as.character(nflId)), by = c("returnerId" = "nflId")) %>% group_by(displayName, returnerId) %>% mutate(success = ifelse(over_expected > 0, 1, 0), boom = ifelse(over_expected >= 2.3428, 1, 0), # 25th %tile bust = ifelse(over_expected <= -6.2811, 1, 0)) %>% #75th %tile summarise(ret = n(), x_yds = sum(adj_x_returner - predicted), yds = sum(adj_x_returner - actual), ryoe = yds - x_yds, yds_ret = yds/ret, ryoe_ret = ryoe/ret, succ_rate = sum(success)/ret, boom_rate = sum(boom)/ret, bust_rate = sum(bust)/ret) %>% filter(ret > 6) %>% left_join( results_rf %>% filter(season == 2020) %>% left_join(df_plays) %>% left_join(df_players %>% mutate(nflId = as.character(nflId)), by = c("returnerId" = "nflId")) %>% mutate(D_close_point = -vor_close_point + lag(vor_close_point), D_adj_x = -adj_x_returner + lag(adj_x_returner), D_x_s = x_s_returner - lag(x_s_returner)) %>% group_by(displayName) %>% summarise(x_s_returner = mean(x_s_returner), D_x_s = mean(D_x_s, na.rm = TRUE), s_returner = mean(s_returner), D_dist_1 = mean(D_dist_1), D_dist_2 = mean(D_dist_2), D_vor_area_12 = mean(D_vor_area_12), vor_area_12 = mean(vor_area_12), D_close_point = mean(D_close_point, na.rm = TRUE), D_adj_x = mean(D_adj_x, na.rm = TRUE))) %>% left_join(nflfastR::fast_scraper_roster(2020) %>% select(full_name, headshot_url), by = c("displayName" = "full_name")) ``` ```{r stability, eval = FALSE, include = FALSE} df_ryoe_stable <- results_rf %>% filter(frameId == min(frameId)) %>% left_join(df_plays) %>% left_join(df_players %>% mutate(nflId = as.character(nflId)), by = c("returnerId" = "nflId")) %>% group_by(displayName, season) %>% mutate(success = ifelse(over_expected > 0, 1, 0), boom = ifelse(over_expected >= 2.3428, 1, 0), bust = ifelse(over_expected <= -6.2811, 1, 0)) %>% #75th %tile summarise(ret = n(), x_yds = sum(adj_x_returner - predicted), yds = sum(adj_x_returner - actual), ryoe = yds - x_yds, yds_ret = yds/ret, ryoe_ret = ryoe/ret, succ_rate = sum(success)/ret, boom_rate = sum(boom)/ret, bust_rate = sum(bust)/ret) %>% filter(ret > 6) %>% group_by(displayName) %>% arrange(season) %>% mutate(ryoe_n_1 = lead(ryoe_ret)) %>% filter(!is.na(ryoe_n_1)) ``` Let's summarize the following key advanced stats for all NFL players who returned at least seven punts in 2020: * RYOE (Return Yards Over Expected) * RYOE/Return: RYOE averaged over all returns * Success Rate: % of returns with **positive** RYOE * Boom Rate: % of returns that went for more than **2.3** RYOE (75th percentile) * Non-Bust Rate: % of returns that DID NOT go for less than **-6.2** RYOE (25th percentile) ``` {r paste ryoe table, echo = FALSE} magick::image_read('../input/tables/ryoe_players.png') ``` Now that we know who were the best returners in 2020, we can look a little deeper. What led to Gunner Olszewski's and Diontae Spencer's dominance? How did Trent Taylor have such a high rate of success? What made CeeDee Lamb so bad? We can start by looking at what leads to a high success rate (using this as opposed to RYOE will be used to rein in skewness, under the assumption that TD scorers are not as heavily rewarded for achieving such a high RYOE on one play). Thus, if we take a returner with a higher success rate to be a more efficient returner, what breeds efficiency? Let's look at the average of our model features for each returner, and how they change over the course of a return. We'll create a correlation matrix with `ggpairs` and include `D_close_point` to capture the rate at which a returner moves his `vor_close_point` closer to the end zone. ```{r ggpairs success rate, echo = FALSE} ggpairs(df_ryoe_player %>% ungroup() %>% select(succ_rate, x_s_returner, s_returner, D_vor_area_12, vor_area_12, D_close_point)) ``` So, interestingly, we see slightly more correlation between a the rate of change of the returner's `vor_close_point` than that of his actual distance from the end zone (`x_s_returner`). Intuitively, we presume that those who "get north and south" see more success in the return game. This is generally true, and the model supports it, but there's an argument to be made for the moves that sacrifice yardage for *space* or *potential* yardage. To dig deeper into this notion, let's look at every 1.5 second interval from each return between 2018 and 2020. There were **49,010**. Further, let's group those returns into four categories: Those that... 1. *improved* **both** `adj_x_returner` and `vor_close_point` 2. *improved* `adj_x_returner` but *worsened* `vor_close_point` 3. *worsened* `adj_x_returner` but *improved* `vor_close_point` 4. *worsened* **both** `adj_x_returner` and `vor_close_point` For the sake of simplicity, I'll refer to `adj_x_returner` as "Position" and `vor_close_point` as "Space" (both measured in distance from the end zone), while `adj_s_returner` (speed towards end zone) and `D_vor_close_point` are just the rates of change of position and space. ```{r diving deeper into adj_x and close_point, include = FALSE} results_rf %>% filter(max(frameId) - min(frameId) >= 15) %>% mutate(adj_x_change = adj_x_returner - lead(adj_x_returner, 15), close_point_change = vor_close_point - lead(vor_close_point, 15), pred_change = predicted - lead(predicted, 15), adj_x_returner = ifelse(adj_x_change > 0, "forwards", ifelse(adj_x_change < 0, "backwards", NA)), vor_close_point = ifelse(close_point_change >= 0, "better", ifelse(close_point_change < 0, "worse", NA))) %>% filter(!is.na(adj_x_returner) & !is.na(vor_close_point)) %>% group_by(adj_x_returner, vor_close_point) %>% summarise(count = n(), avg_pred_change = mean(pred_change)) %>% ungroup() %>% mutate(freq = count/sum(count)) %>% arrange(desc(adj_x_returner)) ``` ``` {r paste space position ryoe table, echo = FALSE} magick::image_read('../input/tables/space_position.png') ``` It's evident that the only pairs, on average, that improve expected return yards are those that improve **space** (if you're curious, this holds true for other reasonable time intervals as well as for individual seasons). Just 17% of instances of backwards movement increased the amount of room in front of the returner, though, so while a sound way to ensure a good return, it's a tough feat. Likewise, 74% of instances that involved forward movement resulted in a better `vor_close_point`, which is expected as they're correlated. Let's animate an example. In this case, Jo-Jo Natson runs backwards and secures enough space in front of him to make for a solid 35-yard return. The red line represents his expectation, while the green line shows the location of his `vor_close_point`. ![](https://media.giphy.com/media/zXYtgvicS3o0OSPTqo/giphy.gif) So, the art of punt returning is not so much a secret, at least from the returner's perspective. As any special teams coach will say, it comes down to three steps: 1. Secure the ball 2. Make a guy miss 3. Run like hell We have, however, uncovered unexpected value in the rate of increase of the "space" that a returner has, relative to the returner's speed. The 17% effectiveness of backwards movement as it relates to space is not substantial, but it facilitates a narrative that supports movement away from the desired end zone in certain cases. Beyond that, speed kills. ```{r plot good returners, echo = FALSE} ggplot(df_ryoe_player) + geom_image(aes(x = x_s_returner, y = succ_rate, image = headshot_url), size = .08, asp = 16/9) + #geom_label(aes(x = x_s_returner, y = succ_rate, label = displayName), vjust = 1.5) + theme_fivethirtyeight() + labs(title = "The Effect of Speed on Punt Returner Success", subtitle = "2020 NFL Season, 7+ Returns", caption = "By Henry Gise | Big Data Bowl 2022", x = "Average Speed Towards Endzone (yds/s)", y = "Success Rate") + scale_y_continuous(labels = percent_format(2), limits = c(0,1), breaks = seq(0,1,.2)) + theme(plot.title = element_text(size = 14, hjust = 0.5, face = "bold"), plot.subtitle = element_text(size = 12), plot.caption = element_text(size = 8), axis.title = element_text(hjust = 0.5, size = 12, face = "bold"), axis.text = element_text(size = 8), aspect.ratio = 9/16) ggplot(df_ryoe_player) + geom_image(aes(x = D_close_point, y = succ_rate, image = headshot_url), size = .08, asp = 16/9) + #geom_label(aes(x = D_close_point, y = succ_rate, label = displayName), vjust = 1.5) + theme_fivethirtyeight() + labs(title = "The Effect of Space on Punt Returner Success", subtitle = "2020 NFL Season, 7+ Returns", caption = "By Henry Gise | Big Data Bowl 2022", x = "Rate of Improvement of vor_close_point", y = "Success Rate") + scale_y_continuous(labels = percent_format(2), limits = c(0,1), breaks = seq(0,1,.2)) + theme(plot.title = element_text(size = 14, hjust = 0.5, face = "bold"), plot.subtitle = element_text(size = 12), plot.caption = element_text(size = 8), axis.title = element_text(hjust = 0.5, size = 12, face = "bold"), axis.text = element_text(size = 8), aspect.ratio = 9/16) ``` # Limitations 1. None of the advanced return metrics are stable on a season-to-season basis. Thus, we can't establish that outperforming expectation indicates returner value, but rather just performance on a play-by-play basis. 2. The only variable that attempts to captture the direction of the returner is `x_s_returner`, which uses his direction of movement and his speed to calculate speed towards the endzone. What we could do is incorporate the direction when we calculate `vor_area_12`, which will penalize a returner's for facing away from the endzone and rewards those are are facing north. 3. We do not capture the horizontal coordinate of `vor_close_point`, meaning the model sees no difference between `vor_close_point` being in front of the returner and being on the other side of the field. Thank you very much for reading.