--- title: "NFL Load Management - Scratching the Surface" output: html_document: code_folding: hide df_print: paged number_sections: yes toc: yes --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = FALSE) ``` # Introduction Nobody likes seeing players get injured. Players don't like them, as they obviously cause tremendous physical pain and the guilt of feeling like they have let down their teammates. Coaches don't like them, because a player in their care has gotten hurt and they now need to make adjustments. Fans don't like them, as the product worsens when the best players cannot take the field, and even opposing fans don't like them, as the 'what if?' scenarios mean there is built in excuses across the board for whatever the result. Unfortunately, injuries are going to be ever present, but that does not mean they cannot be mitigiated. In the last 5 years, the NFL has instituted reform to protect quarterbacks, defenceless receivers and special teams players, some of the most susceptible players in football. The lower leg injury, however, is far more difficult to protect against. Everyone has seen the clips of a player sprinting down the sideline, only to pull up gingerly and hobble to the sideline. Two days later, that player is out for the season. The seemingly unpredictable nature of non-contact, lower leg injuries makes them infuriating, but also creates the need to try and understand what can be done. Random injuries will never be fully eliminated, but we should do what we can to protect those from what we can. In this notebook, I aim to investigate the factor that influence these lower leg injuries in the NFL. This will be done by firstly: 1. Exploring the 1st and Future dataset 2. Find factors that cause injury 3. If possible, build a model to predict if an injury will occur in a given game # Data Preprocessing ```{r Loading Packages, echo = FALSE} suppressPackageStartupMessages(library(plyr)) suppressPackageStartupMessages(library(tidyverse)) suppressPackageStartupMessages(library(janitor)) suppressPackageStartupMessages(library(reshape2)) suppressPackageStartupMessages(library(gridExtra)) suppressPackageStartupMessages(library(zoo)) suppressPackageStartupMessages(library(RcppRoll)) suppressPackageStartupMessages(library(DMwR)) suppressPackageStartupMessages(library(MLmetrics)) ``` ## Simplification of Categorical Variables A number of steps were taken in the processing of the data. Firstly, a number of the categorical variables in the play data was simplified: - Stadium types were simplified into 2 categories: *indoor* and *outdoor*. Stadiums with retractable roofs were classified according to if the roof was open or shut - Weather was simplified into 5 categories: *dry*, *rain_chance*, *rain*, *snow* and *indoor*. rain_chance could have been used as overcast or simply classified as rain for simplicity. ## Temperature Games played indoors did not have a temperature measurement in the data. In this cases, I made assumed the temperature to be 65 degrees. This was not based on any knowledge, but seemed appropriate as a climate controlled environment. ## Removal of Pre and Post Play Tracking Data Due to the large size of the tracking data file, it was not feasible to use the full dataset while conducting analysis. Therefore, I limited the dataset to include only what I thought was necessary. This meant the removal of all tracking data measured before and after the play. For example, most tracking data starts in the huddle (*huddle_break_offense*), with the players moving to the line (*line_set*), the ball being snapped, the play ending, and then some level of post play action, be it moving back to the huddle, subbing out etc. I identified the events that started and finished a play (listed below) and included only data that fell between these events: - Start of play events: - ball_snap - kickoff - free_kick - End of play events: - qb_spike - qb_kneel - touchback - touchdown - fair_catch - field_goal - field_goal_missed - extra_point - extra_point_missed - tackle - pass_outcome_incomplete - out_of_bounds - qb_sack This reduced the from 76,366,748 rows to 13,577,040 rows, which was small enough to begin conducting analysis using RStudio. ## Speed and Acceleration Speed and acceleration were manually calculated from the xy coordinates in the data. The distance travelled per each NGS observation was calculated from the xy coordintes. Then, the speed was calculated by taking this distance travelled and dividing it by the size of the interval (this was generally but not always 0.1 seconds). Finally, the acceleration was then computed as the change in speed over that time period. ```{r Configs, echo = FALSE} # Pre-processing stadium_types = c("Outdoor","Indoors","Oudoor", "Outdoors","Open","Closed Dome", "Domed, closed","Dome", "Indoor","Domed","Retr. Roof-Closed", "Outdoor Retr Roof-Open","Retractable Roof","Ourdoor", "Indoor, Roof Closed","Retr. Roof - Closed","Bowl", "Outddors","Retr. Roof-Open","Dome, closed", "Indoor, Open Roof","Domed, Open","Domed, open" , "Heinz Field","Cloudy","Retr. Roof - Open", "Retr. Roof Closed","Outdor","Outside") stadium_cats = c("outdoor", "indoor", "outdoor", "outdoor", "outdoor", "indoor", "indoor", "indoor", "indoor", "indoor", "indoor", "outdoor", "indoor", "outdoor", "indoor", "indoor", "indoor", "outdoor", "outdoor", "indoor", "outdoor", "outdoor", "outdoor", "outdoor", "outdoor", "outdoor", "indoor", "outdoor", "outdoor") weather_types = c("Clear and warm", "Mostly Cloudy", "Sunny", "Clear", "Cloudy", "Cloudy, fog started developing in 2nd quarter", "Rain", "Partly Cloudy", "Mostly cloudy", "Cloudy and cold", "Cloudy and Cool", "Rain Chance 40%", "Controlled Climate", "Sunny and warm", "Partly cloudy", "Clear and Cool", "Clear and cold", "Sunny and cold", "Indoor", "Partly Sunny", "N/A (Indoors)", "Mostly Sunny", "Indoors", "Clear Skies", "Partly sunny", "Showers", "N/A Indoor", "Sunny and clear", "Snow", "Scattered Showers", "Party Cloudy", "Clear skies", "Rain likely, temps in low 40s.", "Hazy", "Partly Clouidy", "Sunny Skies", "Overcast", "Cloudy, 50% change of rain", "Fair", "Light Rain", "Partly clear", "Mostly Coudy", "10% Chance of Rain", "Cloudy, chance of rain", "Heat Index 95", "Sunny, highs to upper 80s", "Sun & clouds", "Heavy lake effect snow", "Mostly sunny", "Cloudy, Rain", "Sunny, Windy", "Mostly Sunny Skies", "Rainy", "30% Chance of Rain", "Cloudy, light snow accumulating 1-3\"", "cloudy", "Clear and Sunny", "Coudy", "Clear and sunny", "Clear to Partly Cloudy", "Cloudy with periods of rain, thunder possible. Winds shifting to WNW, 10-20 mph.", "Rain shower", "Cold") #clear, overcast, rain, snow, indoor weather_cats = c("dry", "dry", "dry", "dry", "dry", "dry", "rain", "dry", "dry", "dry", "dry", "rain_chance", "indoor", "dry", "dry", "dry", "dry", "dry", "indoor", "dry", "indoor", "dry", "indoor", "dry", "dry", "rain", "dry", "dry", "snow", "rain", "dry", "dry", "rain_chance", "dry", "dry", "dry", "dry", "rain_chance", "dry", "rain", "dry", "dry", "rain", "rain_chance", "dry", "dry", "dry", "snow", "dry", "rain", "dry", "dry", "rain", "rain_chance", "snow", "dry", "dry", "dry", "dry", "dry", "rain", "rain", "dry") # Start of play and end of play events START_OF_PLAY = c("ball_snap", "kickoff", "free_kick") END_OF_PLAY = c("qb_spike", "qb_kneel", "touchback", "touchdown", "fair_catch", "field_goal", "field_goal_missed", "extra_point", "extra_point_missed", "tackle", "pass_outcome_incomplete", "out_of_bounds", "qb_sack") ``` ```{r Import data, echo = FALSE} # Load the data files load("../input/preprocessing-script/preprocessing.rda") ``` # Exploratory Data Analysis ## The Plays and Games Firstly, I conducted essentially a basic stocktake of the dataset. How many plays do we have? How many games do we have? Below examines these two counts and splits them across a number of categorical variables, to see how well covered our sample is. ```{r Initial eda, fig.width = 10, fig.height = 10} # Play by play level player_play_count = play_list %>% group_by(player_key) %>% summarize(number_of_plays = n()) %>% arrange(desc(number_of_plays)) position_play_count = play_list %>% group_by(position_group) %>% summarize(number_of_plays = n()) %>% arrange(desc(number_of_plays)) stadium_play_count = play_list %>% group_by(stadium_simp) %>% summarize(number_of_plays = n()) %>% arrange(desc(number_of_plays)) weather_play_count = play_list %>% group_by(weather_simp) %>% summarize(number_of_plays = n()) %>% arrange(desc(number_of_plays)) surface_play_count = play_list %>% group_by(field_type) %>% summarize(number_of_plays = n()) %>% arrange(desc(number_of_plays)) number_of_players = play_list %>% group_by(player_key, position_group) %>% summarize(count = n()) %>% group_by(player_key) %>% summarize(primary_position = position_group[count == max(count)]) %>% group_by(primary_position) %>% summarize(players = n()) %>% arrange(desc(players)) player_play_plot = ggplot(player_play_count, aes(x = factor(player_key, levels = player_key), y = number_of_plays)) + geom_col(fill = "darkorchid3") + theme_minimal() + theme(axis.text.x = element_text(angle = 90, size = 4)) + labs(title = "Player per Player", x = "Player Key", y = "Number of Plays") position_play_plot = ggplot(position_play_count, aes(x = factor(position_group, levels = position_group), y = number_of_plays)) + geom_col(fill = "darkorchid3", colour = "black") + theme_minimal() + theme(axis.text.x = element_text(angle = 45, size = 8)) + labs(title = "Plays per Position Group", x = "Position Group", y = "Number of Plays") stadium_play_count_plot = ggplot(stadium_play_count, aes(x = factor(stadium_simp, levels = stadium_simp), y = number_of_plays)) + geom_col(fill = "darkorchid3", colour = "black") + theme_minimal() + theme(axis.text.x = element_text(angle = 45, size = 8)) + labs(title = "Plays per Stadium Type", x = "Stadium Type", y = "Number of Plays") weather_play_count_plot = ggplot(weather_play_count, aes(x = factor(weather_simp, levels = weather_simp), y = number_of_plays)) + geom_col(fill = "darkorchid3", colour = "black") + theme_minimal() + theme(axis.text.x = element_text(angle = 45, size = 8)) + labs(title = "Plays per Weather Type", x = "Weather Type", y = "Number of Plays") surface_play_count_plot = ggplot(surface_play_count, aes(x = factor(field_type, levels = field_type), y = number_of_plays)) + geom_col(fill = "darkorchid3", colour = "black") + theme_minimal() + theme(axis.text.x = element_text(angle = 45, size = 8)) + labs(title = "Plays per Surface Type", x = "Surface Type", y = "Number of Plays") # Plots for the game level player_game_count = game_list %>% group_by(player_key) %>% summarize(number_of_games = n()) %>% arrange(desc(number_of_games)) position_game_count = game_list %>% group_by(position_group) %>% summarize(number_of_games = n()) %>% arrange(desc(number_of_games)) stadium_game_count = game_list %>% group_by(stadium_simp) %>% summarize(number_of_games = n()) %>% arrange(desc(number_of_games)) weather_game_count = game_list %>% group_by(weather_simp) %>% summarize(number_of_games = n()) %>% arrange(desc(number_of_games)) surface_game_count = game_list %>% group_by(field_type) %>% summarize(number_of_games = n()) %>% arrange(desc(number_of_games)) player_game_count_plot = ggplot(player_game_count, aes(x = factor(player_key, levels = player_key), y = number_of_games)) + geom_col(fill = "goldenrod3") + theme_minimal() + theme(axis.text.x = element_text(angle = 90, size = 4)) + labs(title = "Games per Player", x = "Player Key", y = "Number of Games") position_game_count_plot = ggplot(position_game_count, aes(x = factor(position_group, levels = position_group), y = number_of_games)) + geom_col(fill = "goldenrod3", colour = "black") + theme_minimal() + theme(axis.text.x = element_text(angle = 45, size = 8)) + labs(title = "Games per Position Type", x = "Position Group", y = "Number of Games") stadium_game_count_plot = ggplot(stadium_game_count, aes(x = factor(stadium_simp, levels = stadium_simp), y = number_of_games)) + geom_col(fill = "goldenrod3", colour = "black") + theme_minimal() + theme(axis.text.x = element_text(angle = 45, size = 8)) + labs(title = "Games per Stadium Type", x = "Stadium Type", y = "Number of Games") weather_game_count_plot = ggplot(weather_game_count, aes(x = factor(weather_simp, levels = weather_simp), y = number_of_games)) + geom_col(fill = "goldenrod3", colour = "black") + theme_minimal() + theme(axis.text.x = element_text(angle = 45, size = 8)) + labs(title = "Games per Weather Type", x = "Weather Type", y = "Number of Games") surface_game_count_plot = ggplot(surface_game_count, aes(x = factor(field_type, levels = field_type), y = number_of_games)) + geom_col(fill = "goldenrod3", colour = "black") + theme_minimal() + theme(axis.text.x = element_text(angle = 45, size = 8)) + labs(title = "Games per Surface Type", x = "Surface Type", y = "Number of Games") grid.arrange(player_play_plot, player_game_count_plot, position_play_plot, position_game_count_plot, stadium_play_count_plot, stadium_game_count_plot, weather_play_count_plot, weather_game_count_plot, surface_play_count_plot, surface_game_count_plot, nrow = 5) ``` Naturally, the play level data largely reflects the game level data, as of course the number of plays in a game is relatively constant. However, there are some differences in relation to the player positions. For example, despite having more plays for Offensive Linemen, we actually have more recorded games for Wide Receivers, likely due to the fact that Linemen play more snaps in a game than receivers. ## The Injuries Given that a major focus of this project is injuries, we should also look at our Injury Record. The data provided contains 105 recorded injuries on 104 plays (poor player 47307 hurt both his knee and ankle on the same play). ```{r, fig.height = 6, fig.width = 10} injury_games = game_list %>% filter(injury_in_game == 1) injury_position_game_count = injury_games %>% group_by(position_group) %>% summarize(number_of_games = n()) %>% arrange(desc(number_of_games)) injury_stadium_game_count = injury_games %>% group_by(stadium_simp) %>% summarize(number_of_games = n()) %>% arrange(desc(number_of_games)) injury_weather_game_count = injury_games %>% group_by(weather_simp) %>% summarize(number_of_games = n()) %>% arrange(desc(number_of_games)) injury_surface_game_count = injury_games %>% group_by(field_type) %>% summarize(number_of_games = n()) %>% arrange(desc(number_of_games)) injury_position_game_count_plot = ggplot(injury_position_game_count, aes(x = factor(position_group, levels = position_group), y = number_of_games)) + geom_col(fill = "aquamarine4", colour = "black") + theme_minimal() + theme(axis.text.x = element_text(angle = 45, size = 8)) + labs(title = "Games per Position Type", subtitle = "Games where Injuries Occurred", x = "Position Group", y = "Number of Games") injury_stadium_game_count_plot = ggplot(injury_stadium_game_count, aes(x = factor(stadium_simp, levels = stadium_simp), y = number_of_games)) + geom_col(fill = "aquamarine4", colour = "black") + theme_minimal() + theme(axis.text.x = element_text(angle = 45, size = 8)) + labs(title = "Games per Stadium Type", subtitle = "Games where Injuries Occurred", x = "Stadium Type", y = "Number of Games") injury_weather_game_count_plot = ggplot(injury_weather_game_count, aes(x = factor(weather_simp, levels = weather_simp), y = number_of_games)) + geom_col(fill = "aquamarine4", colour = "black") + theme_minimal() + theme(axis.text.x = element_text(angle = 45, size = 8)) + labs(title = "Games per Weather Type", subtitle = "Games where Injuries Occurred", x = "Weather Type", y = "Number of Games") injury_surface_game_count_plot = ggplot(injury_surface_game_count, aes(x = factor(field_type, levels = field_type), y = number_of_games)) + geom_col(fill = "aquamarine4", colour = "black") + theme_minimal() + theme(axis.text.x = element_text(angle = 45, size = 8)) + labs(title = "Games per Surface Type", subtitle = "Games where Injuries Occurred", x = "Surface Type", y = "Number of Games") grid.arrange(injury_position_game_count_plot, injury_stadium_game_count_plot, injury_weather_game_count_plot, injury_surface_game_count_plot, nrow = 2) ``` There are some important things to note here: - We do not have injuries for every position group. For example, no Quarterback injury was recorded in the sample. The NFL rules have changed in recent years, but not so much to fully prevent Quarterback injuries (as much as pass rushers might think they have). Any model that comes out of this should be smart enough to realise this. - The weather and stadium information largely resemble the distributions of the entire dataset. - The majority of injuries took place in games played on a synthetic surface, despite the fact that more games in the data set were played on natural grass. This leads credence to the thought that synthetic turf does cause injury at a higher rate. It should also be noted, injuries were recorded in only 104 of about 267000 plays! Injuries are still a fairly rare event, and this class imbalance will be addressed later on, watch this space. ## Tracking Data Play level data is great, but the NGS tracking data is the real meat of this project. Next we dive into the player tracking to gain insights into the lengths of plays and player movement in a range of situations. ### Length of Plays We firstly look at the length of each play: ```{r Build play summaries df} YPS_5 = 5 YPS_6 = 6 YPS_7 = 7 YPS_8 = 8 ACC_1 = 1 ACC_1.5 = 1.5 ACC_2 = 2 ACC_2.5 = 2.5 ACC_3 = 3 play_summaries = player_track_data %>% group_by(play_key) %>% summarize(play_outcome = first(play_outcome), play_length = max(in_play_time), max_speed_on_play = max(s, na.rm = TRUE), first_time_of_max_speed = first(in_play_time[max_speed_on_play == s]), distance_travelled_on_play = sum(dis), avg_acceleration_first_second = mean(a[in_play_time <= 1], na.rm = TRUE), avg_acceleration_two_seconds = mean(a[in_play_time <= 2], na.rm = TRUE), avg_acceleration_three_seconds = mean(a[in_play_time <= 3], na.rm = TRUE), avg_acceleration_four_seconds = mean(a[in_play_time <= 4], na.rm = TRUE), avg_acceleration_five_seconds = mean(a[in_play_time <= 5], na.rm = TRUE), mag_avg_acceleration_first_second = mean(abs(a[in_play_time <= 1]), na.rm = TRUE), mag_avg_acceleration_two_seconds = mean(abs(a[in_play_time <= 2]), na.rm = TRUE), mag_avg_acceleration_three_seconds = mean(abs(a[in_play_time <= 3]), na.rm = TRUE), mag_avg_acceleration_four_seconds = mean(abs(a[in_play_time <= 4]), na.rm = TRUE), mag_avg_acceleration_five_seconds = mean(abs(a[in_play_time <= 5]), na.rm = TRUE), distance_per_second_on_play = distance_travelled_on_play/play_length, time_spent_faster_5 = sum(s > YPS_5, na.rm = TRUE) * 0.1, time_spent_faster_6 = sum(s > YPS_6, na.rm = TRUE) * 0.1, time_spent_faster_7 = sum(s > YPS_7, na.rm = TRUE) * 0.1, time_spent_faster_8 = sum(s > YPS_8, na.rm = TRUE) * 0.1, time_spent_acc_1 = sum(a > ACC_1, na.rm = TRUE) * 0.1, time_spent_acc_1.5 = sum(a > ACC_1.5, na.rm = TRUE) * 0.1, time_spent_acc_2 = sum(a > ACC_2, na.rm = TRUE) * 0.1, time_spent_acc_2.5 = sum(a > ACC_2.5, na.rm = TRUE) * 0.1, time_spent_acc_3 = sum(a > ACC_3, na.rm = TRUE) * 0.1, time_spent_dec_1 = sum(a < -ACC_1, na.rm = TRUE) * 0.1, time_spent_dec_1.5 = sum(a < -ACC_1.5, na.rm = TRUE) * 0.1, time_spent_dec_2 = sum(a < -ACC_2, na.rm = TRUE) * 0.1, time_spent_dec_2.5 = sum(a < -ACC_2.5, na.rm = TRUE) * 0.1, time_spent_dec_3 = sum(a < -ACC_3, na.rm = TRUE) * 0.1, max_acceleration = max(a, na.rm=TRUE), max_deceleration = -min(a, na.rm=TRUE)) %>% left_join(play_list %>% select(game_id, player_key, play_key, play_type, player_day, player_game, player_game_play, position_group, field_type, temperature, stadium_simp, weather_simp, injury_in_game, injury_on_play), by = "play_key") %>% arrange(player_key, player_game, player_game_play) ``` ```{r Play lengths by type of play, fig.width = 10, fig.height = 4} # Play lengths total_play_length = ggplot(play_summaries, aes(x = play_length)) + geom_density(fill = "darkorchid3") + labs(title = "Length of Each Play in the Dataset", x = "Length of Play from Snap to Play Ending Event", y = "Density") + theme_minimal() play_length_by_type = ggplot(play_summaries, aes(x = play_length)) + geom_density(aes(fill = play_type, y = ..density..)) + labs(title = "Length of Play by Play Type", x = "Length of Play from Snap to Play Ending Event", y = "Density") + facet_wrap(~play_type, nrow = 3) + guides(fill = FALSE, colour = FALSE) + theme_minimal() grid.arrange(total_play_length, play_length_by_type, ncol = 2) ``` From the ball snap to the play ending event, the mean play length was 5.1 seconds, with the majority of plays being less than 10 seconds. There are of course a few outliers - the longest play in the data set was 43.8 seconds, occurring 3 times, with each case being end of game lateral plays (34243-7-52, 42346-7-63, 43490-7-77). I suspect this was the same play. Similarly, there are some extremely short plays < 2 seconds. These are generally spikes and kneeldowns. When broken up by the type of play, we see the difference in density plots across each type. Some observations: - Extra Points are extremely short as expected - most extra points are extremely routine and I suspect it is the two point conversions giving a skew - Kickoffs and Punts are multi-modal, likely skewed by the size of the return - Rushing plays appear on average to take less time than passing plays, with a small second mode for rushing plays that are QB kneels - There were two other play types in the data, "0" and NA. I did not have time to clean these, but left them here anyway. They do not appear to have any defining characteristics. ### Maximum Speed Next, I looked at the maximum speed that was achieved by each player on each play: ```{r Maximum speed reached on play, fig.width = 10, fig.height = 4} max_speed_by_play_plot = ggplot(play_summaries, aes(x = max_speed_on_play)) + geom_density(aes(fill = play_type, y = ..density..)) + labs(title = "Maximum Speed on Play", subtitle = "By Play Type", x = "Max Speed (yards/second)", y = "Density") + guides(fill = FALSE) + theme_minimal() + facet_wrap(~play_type, nrow = 3) max_speed_by_position_plot = ggplot(play_summaries, aes(x = max_speed_on_play)) + geom_density(aes(fill = position_group, y = ..density..)) + labs(title = "Maximum Speed on Play", subtitle = "By Position Group", x = "Max Speed (yards/second)", y = "Density") + guides(fill = FALSE) + theme_minimal() + facet_wrap(~position_group, nrow = 3) grid.arrange(max_speed_by_play_plot, max_speed_by_position_plot, nrow = 1) ``` Players appear to achieve faster top speeds on special teams plays such as Kickoffs and Punts. This makes sense, as these plays generally feature at least one side sprinting as hard as they possibly can to cover further down the field. We additionally see the difference in speed between the positions. Wide Receivers, Running Backs and Defensive Backs appear to reach a range of values, while Offensive Linemen and Quarterbacks consistently reach lower top speeds which makes sense given their roles in the game. That is not to say that all players of the same position are created equally. The following plots show differences between individual players in the same position group. We then look at the number of plays each player in the dataset played, combined with the distance they travelled. ```{r Number of plays, fig.width = 10, fig.height = 4} # Count the number of plays per player play_count = play_summaries %>% group_by(player_key, position_group) %>% summarise(plays = n(), total_distance_travelled = sum(distance_travelled_on_play), average_distance_per_play = mean(distance_travelled_on_play), distance_per_play_variance = var(distance_travelled_on_play), average_max_speed_per_play = mean(max_speed_on_play), average_max_speed_variance = var(max_speed_on_play)) distance_per_play_plot = ggplot(play_count, aes(x = position_group, y = average_distance_per_play)) + geom_violin(aes(fill = position_group)) + geom_point(aes(fill = position_group), size = 0.2) + labs(title = "Average Distance Covered on Play by Player", subtitle = "Facetted by Position", x = "Position Group", y = "Average Yards Travelled per Play") + theme_minimal() speed_per_play_plot = ggplot(play_count, aes(x = position_group, y = average_max_speed_per_play)) + geom_violin(aes(fill = position_group)) + geom_point(aes(fill = position_group), size = 0.2) + labs(title = "Max Speed Per Play by Player", x = "Position Group", y = "Average Yards/Second") + theme_minimal() grid.arrange(distance_per_play_plot, speed_per_play_plot, ncol = 2) ``` There are no surprises here, with skill players covering more distance and reaching higher speeds than players rooted to the line of scrimmage. ### Differences in Movement on Different Surfaces We now compare the movement on a natural surface to that on a synthetic. Does it look like there is a difference? ```{r, fig.width=10, fig.height=4} max_speed_surface_plot = ggplot(play_summaries, aes(x = max_speed_on_play, group = field_type)) + geom_density(aes(fill = field_type), alpha = 0.6) + labs(title = "Maximum Speed on Play", subtitle = "Synthetic vs Natural Grass", x = "Maximum Speed (yards/sec)", y = "Density", fill = "Field Type") + theme_minimal() time_max_speed_surface_plot = ggplot(play_summaries, aes(x = first_time_of_max_speed, group = field_type)) + geom_density(aes(fill = field_type), alpha = 0.6) + labs(title = "Time to Reach Max Speed", subtitle = "Synthetic vs Natural Grass", x = "Time of Max Speed", y = "Density", fill = "Field Type") + theme_minimal() avg_acc_second_1_plot = ggplot(play_summaries, aes(x = mag_avg_acceleration_first_second, group = field_type)) + geom_density(aes(fill = field_type), alpha = 0.6) + labs(title = "Average Acceleration over First Second", subtitle = "Synthetic vs Natural Grass", x = "Avg. Acceleration", y = "Density", fill = "Field Type") + theme_minimal() avg_acc_second_2_plot = ggplot(play_summaries, aes(x = mag_avg_acceleration_two_seconds, group = field_type)) + geom_density(aes(fill = field_type), alpha = 0.6) + labs(title = "Average Acceleration over First 2 Seconds", subtitle = "Synthetic vs Natural Grass", x = "Time of Max Speed", y = "Density", fill = "Field Type") + theme_minimal() avg_acc_second_3_plot = ggplot(play_summaries, aes(x = mag_avg_acceleration_three_seconds, group = field_type)) + geom_density(aes(fill = field_type), alpha = 0.6) + labs(title = "Average Acceleration over First 3 Seconds", subtitle = "Synthetic vs Natural Grass", x = "Time of Max Speed", y = "Density", fill = "Field Type") + theme_minimal() avg_acc_second_4_plot = ggplot(play_summaries, aes(x = mag_avg_acceleration_four_seconds, group = field_type)) + geom_density(aes(fill = field_type), alpha = 0.6) + labs(title = "Average Acceleration over First 4 Seconds", subtitle = "Synthetic vs Natural Grass", x = "Time of Max Speed", y = "Density", fill = "Field Type") + theme_minimal() grid.arrange(max_speed_surface_plot, time_max_speed_surface_plot, ncol = 2) ``` The maximum speed reached on each play appears the same, as well as the time needed to get there. However, I think acceleration is of particular interest. One of the perceived advantages of a synthetic turf is increased friction with the players foot, meaning in theory they should be able to push off faster. Here I look at the average acceleration of each player over the first 1 and 2 seconds of each play to determine if there is this initial pushoff. ```{r, fig.width=10, fig.height = 4} grid.arrange(avg_acc_second_1_plot, avg_acc_second_2_plot, ncol = 2) ``` Again, there appears to be minimal difference. Finally, I tried separating these across position to see if perhaps some positions are affected and others are not. ```{r, fig.width = 10, fig.height = 4} # are these differences present across position? max_speed_pos = ggplot(play_summaries %>% filter(max_speed_on_play <= 12), aes(x = max_speed_on_play, group = field_type)) + geom_density(aes(fill = field_type), alpha = 0.6) + labs(title = "Maximum Speed on Play", subtitle = "Across Position, Synthetic vs Natural Grass", x = "Maximum Speed (yards/sec)", y = "Density", fill = "Field Type") + theme_minimal() + facet_wrap(~position_group, ncol = 2) time_max_speed_pos = ggplot(play_summaries %>% filter(first_time_of_max_speed <= 15), aes(x = first_time_of_max_speed, group = field_type)) + geom_density(aes(fill = field_type), alpha = 0.6) + labs(title = "Time to Reach Max Speed", subtitle = "Across Position, Synthetic vs Natural Grass", x = "Time of Max Speed", y = "Density", fill = "Field Type") + theme_minimal() + facet_wrap(~position_group, ncol = 2) pos_avg_acc_second_1_plot = ggplot(play_summaries, aes(x = mag_avg_acceleration_first_second, group = field_type)) + geom_density(aes(fill = field_type), alpha = 0.6) + labs(title = "Average Acceleration over First Second", subtitle = "Across Position, Synthetic vs Natural Grass", x = "Avg. Acceleration", y = "Density", fill = "Field Type") + theme_minimal() + facet_wrap(~position_group, ncol = 2) pos_avg_acc_second_2_plot = ggplot(play_summaries, aes(x = mag_avg_acceleration_two_seconds, group = field_type)) + geom_density(aes(fill = field_type), alpha = 0.6) + labs(title = "Average Acceleration over First 2 Seconds", subtitle = "Across Position, Synthetic vs Natural Grass", x = "Time of Max Speed", y = "Density", fill = "Field Type") + theme_minimal() + facet_wrap(~position_group, ncol = 2) grid.arrange(max_speed_pos, time_max_speed_pos, ncol = 2) grid.arrange(pos_avg_acc_second_1_plot, pos_avg_acc_second_2_plot, ncol = 2) ``` There could be slight differences across quarterbacks, running backs and wide receivers, although these differences look slight. In all, in terms of speed and acceleraiton there appears to be minimal difference between the synthetic and natural surfaces. That is not to say these effects do not exist (more on that in a second), but they are not overtly obvious. To be fair, if they were obvious we probably would have found them already. # Exploring Factors that Influence Injury ## Research There has been a lot of research into the causes of injury across a number of sports. It is currently an ongoing topic of debate in the NBA, with numerous teams citing 'load management' as a reason to rest players. One particular paper was by Cummins et. al (2019), available [here](https://www.ncbi.nlm.nih.gov/pubmed/30651223). I found helpful focussed on Rugby League, the sport popular in Australia. Although Rugby League is an inherently different game to American Football (continuous plays, less substitutions, less difference between positions), their approach in trying to measure injury risk seemed interesting. They largely focussed on quantifying the level of physical exertion of a player over rolling historical windows, incorporating both training and game information. This concept of measuring a rolling level of physical exertion provided the backbone of the approach I would take with the training data. Similarly, much research has been conducted into the effect of the playing surface with some suggesting that synthetic grass does contribute to injury. As a result, I aimed to combine these two pieces of information relating to causes of injury. ## Approach The approach I decided to take was to build a logistic regression model that would predict the probability that an injury would be sustained in a particular game for a particular player. Although it would be amazing if these predictions were actually accurate in predicting if a player could be injured, the main reasoning for this is to use the parameter estimates to determine which features have a significant influence on the target. If it was statistically significant, then that feature does have an effect, either positively or negatively, on the onset of injury. The sign of the coefficient can be used to determine whether it was positive or negative. The reasoning for this approach was as follows: - **It would have practical use for a player or team.** - Being able to accurately quantify the risk of a player being injured could affect how they are used in the game - It is not particularly helpful to be given the data for a play and predict if an injury happened on that play after the fact. A player in the middle of a play will not be thinking, *"I've been accelerating a lot on this play. I should slow down and let this guy run past me and score a touchdown"* Doing it at the game level allows for a plan to be put in place by a staff *before* the game. - **We can incorporate numeric features** - Much of the approach in the literature has involved using chi-squared tests for categorical data. ie. `Injury Occurred vs Did Not Occur` with `Synthetic vs Grass`. Using a linear model allows us to incorpoate these effects. - - **We did not have complete injury data at the play level** - In the `PlayList` data frame, we had 105 separate injuries, but only 77 injuries had the recorded play where the injury occurred. - Injuries are such rare events that trying to predict if an injury occurred on the play would have been extremely difficult. In building the model, I ignore both the type of injury and the size of the injury sustained. This was to increase the sample of injuries (so they were not separate), and also to keep the problem as simple as possible. ## Methodology The model chosen to do this was a generalised linear model with a logit link function (logistic regression). Although black box methods such as XGBoost or Neural Networks would provide more accurate predictions, the purpose of this project is *interpretability*. Not only do we want to assess the risk of injury, but how much each individual factor is contributing, which is possible within this framework. There were two categories of features that were used to predict the likelihood of injuries: *Environmental Variables*, those referring to the context of where the game was played, and *Load Variables*, referring to the physical toll a player had sustained in previous games. Although we did not have access to the training information of players, I hoped that the game information would be sufficient in measuring a players workload. The location variables were fairly straightforward to devise: **Environmental Variables:** - Field type - Stadium type (indoor or outdoor) - Weather type (dry, rain, rain chance, snow, indoor) - Temperature **Load Variables:** Rolling windows looking back at the previous game, previous two games, and previous four games aim to capture both the long term lingering effects of exertion as well as more recent effects. I wanted to look at how much time a player had spent running at faster than 5 and 7 yards per second, accelerating at faster than 1.5 and 2.5 yards per second per second, and decelerating at faster than 1.5 and 2.5 yards per second per second. 5 yards/sec and 7 yards/sec were arbitrarily chosen to measure high and very high running speeds, as they are somewhat close to the 20km/hr and 25km/hr benchmarks used in the Rugby League paper. Similarly, the acceleration values of 1.5 and 2.5 yards per second per second were arbitrarily chosen, and these values could certainly be tuned in future to find a critical point where injuries occur. Additionally, I included the average distance travelled per play of the player, and the number of games they had played on a synthetic surface over that window. So, the final features were: - Average distance per play in last 1, 2, 4 games - Time spent running at > 5 yards per second in last 1, 2, 4 games - Time spent running at > 7 yards per second in last 1, 2, 4 games - Time spent with acceleration > 1.5 yards per second per second in last 1, 2, 4 games - Time spent with acceleration > 2.5 yards per second per second in last 1, 2, 4 games - Time spent with deceleration > 1.5 yards per second per second in last 1, 2, 4 games - Time spent with deceleration > 2.5 yards per second per second in last 1, 2, 4 games - Number of games played on synthetic turf in last 1, 2, 4 games ## Model Fits ```{r, fig.height = 4, fig.width = 10, warning = FALSE} YPS_5 = 5 YPS_6 = 6 YPS_7 = 7 YPS_8 = 8 ACC_1 = 1 ACC_1.5 = 1.5 ACC_2 = 2 ACC_2.5 = 2.5 ACC_3 = 3 play_summaries = player_track_data %>% group_by(play_key) %>% summarize(play_outcome = first(play_outcome), play_length = max(in_play_time), max_speed_on_play = max(s, na.rm = TRUE), first_time_of_max_speed = first(in_play_time[max_speed_on_play == s]), distance_travelled_on_play = sum(dis), avg_acceleration_first_second = mean(a[in_play_time <= 1], na.rm = TRUE), avg_acceleration_two_seconds = mean(a[in_play_time <= 2], na.rm = TRUE), avg_acceleration_three_seconds = mean(a[in_play_time <= 3], na.rm = TRUE), avg_acceleration_four_seconds = mean(a[in_play_time <= 4], na.rm = TRUE), avg_acceleration_five_seconds = mean(a[in_play_time <= 5], na.rm = TRUE), mag_avg_acceleration_first_second = mean(abs(a[in_play_time <= 1]), na.rm = TRUE), mag_avg_acceleration_two_seconds = mean(abs(a[in_play_time <= 2]), na.rm = TRUE), mag_avg_acceleration_three_seconds = mean(abs(a[in_play_time <= 3]), na.rm = TRUE), mag_avg_acceleration_four_seconds = mean(abs(a[in_play_time <= 4]), na.rm = TRUE), mag_avg_acceleration_five_seconds = mean(abs(a[in_play_time <= 5]), na.rm = TRUE), distance_per_second_on_play = distance_travelled_on_play/play_length, time_spent_faster_5 = sum(s > YPS_5, na.rm = TRUE) * 0.1, time_spent_faster_6 = sum(s > YPS_6, na.rm = TRUE) * 0.1, time_spent_faster_7 = sum(s > YPS_7, na.rm = TRUE) * 0.1, time_spent_faster_8 = sum(s > YPS_8, na.rm = TRUE) * 0.1, time_spent_acc_1 = sum(a > ACC_1, na.rm = TRUE) * 0.1, time_spent_acc_1.5 = sum(a > ACC_1.5, na.rm = TRUE) * 0.1, time_spent_acc_2 = sum(a > ACC_2, na.rm = TRUE) * 0.1, time_spent_acc_2.5 = sum(a > ACC_2.5, na.rm = TRUE) * 0.1, time_spent_acc_3 = sum(a > ACC_3, na.rm = TRUE) * 0.1, time_spent_dec_1 = sum(a < -ACC_1, na.rm = TRUE) * 0.1, time_spent_dec_1.5 = sum(a < -ACC_1.5, na.rm = TRUE) * 0.1, time_spent_dec_2 = sum(a < -ACC_2, na.rm = TRUE) * 0.1, time_spent_dec_2.5 = sum(a < -ACC_2.5, na.rm = TRUE) * 0.1, time_spent_dec_3 = sum(a < -ACC_3, na.rm = TRUE) * 0.1) %>% left_join(play_list %>% select(game_id, player_key, play_key, play_type, player_day, player_game, player_game_play, position_group, field_type, temperature, stadium_simp, weather_simp, injury_in_game, injury_on_play), by = "play_key") %>% arrange(player_key, player_game, player_game_play) # Construct game level summaries # Construct game level summaries game_level_injury_data = play_summaries %>% group_by(game_id, player_key, player_game, player_day, position_group, field_type, stadium_simp, weather_simp, temperature, injury_in_game) %>% summarize(plays_in_game = n(), distance_travelled_in_game = sum(distance_travelled_on_play), distance_per_play = distance_travelled_in_game/plays_in_game, average_max_speed = mean(max_speed_on_play), avg_distance_per_second = mean(distance_per_second_on_play), total_time_faster_5 = sum(time_spent_faster_5), total_time_faster_6 = sum(time_spent_faster_6), total_time_faster_7 = sum(time_spent_faster_7), total_time_faster_8 = sum(time_spent_faster_8), total_time_acc_1.5 = sum(time_spent_acc_1.5), total_time_acc_2.5 = sum(time_spent_acc_2.5), total_time_dec_1.5 = sum(time_spent_dec_1.5), total_time_dec_2.5 = sum(time_spent_dec_2.5), is_synthetic_game = first(ifelse(field_type == "Synthetic", 1, 0))) %>% arrange(player_key, player_game, player_day) training_data = game_level_injury_data %>% group_by(player_key) %>% mutate(distance_per_play_last_1_game = lag(roll_sum(distance_per_play, n = 1, fill = NA, align = "right")), distance_per_play_last_2_game = lag(roll_sum(distance_per_play, n = 2, fill = NA, align = "right")), distance_per_play_last_4_game = lag(roll_sum(distance_per_play, n = 4, fill = NA, align = "right")), time_faster_5_last_1_game = lag(roll_sum(total_time_faster_5, n = 1, fill = NA, align = "right")), time_faster_5_last_2_game = lag(roll_sum(total_time_faster_5, n = 2, fill = NA, align = "right")), time_faster_5_last_4_game = lag(roll_sum(total_time_faster_5, n = 4, fill = NA, align = "right")), time_faster_7_last_1_game = lag(roll_sum(total_time_faster_7, n = 1, fill = NA, align = "right")), time_faster_7_last_2_game = lag(roll_sum(total_time_faster_7, n = 2, fill = NA, align = "right")), time_faster_7_last_4_game = lag(roll_sum(total_time_faster_7, n = 4, fill = NA, align = "right")), time_acc_1.5_last_1_game = lag(roll_sum(total_time_acc_1.5, n = 1, fill = NA, align = "right")), time_acc_1.5_last_2_game = lag(roll_sum(total_time_acc_1.5, n = 2, fill = NA, align = "right")), time_acc_1.5_last_4_game = lag(roll_sum(total_time_acc_1.5, n = 4, fill = NA, align = "right")), time_acc_2.5_last_1_game = lag(roll_sum(total_time_acc_2.5, n = 1, fill = NA, align = "right")), time_acc_2.5_last_2_game = lag(roll_sum(total_time_acc_2.5, n = 2, fill = NA, align = "right")), time_acc_2.5_last_4_game = lag(roll_sum(total_time_acc_2.5, n = 4, fill = NA, align = "right")), time_dec_1.5_last_1_game = lag(roll_sum(total_time_dec_1.5, n = 1, fill = NA, align = "right")), time_dec_1.5_last_2_game = lag(roll_sum(total_time_dec_1.5, n = 2, fill = NA, align = "right")), time_dec_1.5_last_4_game = lag(roll_sum(total_time_dec_1.5, n = 4, fill = NA, align = "right")), time_dec_2.5_last_1_game = lag(roll_sum(total_time_dec_2.5, n = 1, fill = NA, align = "right")), time_dec_2.5_last_2_game = lag(roll_sum(total_time_dec_2.5, n = 2, fill = NA, align = "right")), time_dec_2.5_last_4_game = lag(roll_sum(total_time_dec_2.5, n = 4, fill = NA, align = "right")), last_game_synthetic_1 = lag(roll_sum(is_synthetic_game, n = 1, fill = NA, align = "right")), last_game_synthetic_2 = lag(roll_sum(is_synthetic_game, n = 2, fill = NA, align = "right")), last_game_synthetic_4 = lag(roll_sum(is_synthetic_game, n = 4, fill = NA, align = "right"))) %>% # Doing this basically removes the first 4 games for each player. Approximate based on these group_by(player_key) %>% mutate(distance_per_play_last_1_game = ifelse(is.na(distance_per_play_last_1_game), 0, distance_per_play_last_1_game), distance_per_play_last_2_game = ifelse(is.na(distance_per_play_last_2_game), distance_per_play_last_1_game, distance_per_play_last_2_game), distance_per_play_last_4_game = ifelse(is.na(distance_per_play_last_4_game), distance_per_play_last_1_game, distance_per_play_last_4_game), time_faster_5_last_1_game = ifelse(is.na(time_faster_5_last_1_game), 0, time_faster_5_last_1_game), time_faster_5_last_2_game = ifelse(is.na(time_faster_5_last_2_game), time_faster_5_last_1_game, time_faster_5_last_2_game), time_faster_5_last_4_game = ifelse(is.na(time_faster_5_last_4_game), time_faster_5_last_1_game, time_faster_5_last_4_game), time_faster_7_last_1_game = ifelse(is.na(time_faster_7_last_1_game), 0, time_faster_7_last_1_game), time_faster_7_last_2_game = ifelse(is.na(time_faster_7_last_2_game), time_faster_7_last_1_game, time_faster_7_last_2_game), time_faster_7_last_4_game = ifelse(is.na(time_faster_7_last_4_game), time_faster_7_last_1_game, time_faster_7_last_4_game), time_acc_1.5_last_1_game = ifelse(is.na(time_acc_1.5_last_1_game), 0, time_acc_1.5_last_1_game), time_acc_1.5_last_2_game = ifelse(is.na(time_acc_1.5_last_2_game), time_acc_1.5_last_1_game, time_acc_1.5_last_2_game), time_acc_1.5_last_4_game = ifelse(is.na(time_acc_1.5_last_4_game), time_acc_1.5_last_1_game, time_acc_1.5_last_4_game), time_acc_2.5_last_1_game = ifelse(is.na(time_acc_2.5_last_1_game), 0, time_acc_2.5_last_1_game), time_acc_2.5_last_2_game = ifelse(is.na(time_acc_2.5_last_2_game), time_acc_2.5_last_1_game, time_acc_2.5_last_2_game), time_acc_2.5_last_4_game = ifelse(is.na(time_acc_2.5_last_4_game), time_acc_2.5_last_1_game, time_acc_2.5_last_4_game), time_dec_1.5_last_1_game = ifelse(is.na(time_dec_1.5_last_1_game), 0, time_dec_1.5_last_1_game), time_dec_1.5_last_2_game = ifelse(is.na(time_dec_1.5_last_2_game), time_dec_1.5_last_1_game, time_dec_1.5_last_2_game), time_dec_1.5_last_4_game = ifelse(is.na(time_dec_1.5_last_4_game), time_dec_1.5_last_1_game, time_dec_1.5_last_4_game), time_dec_2.5_last_1_game = ifelse(is.na(time_dec_2.5_last_1_game), 0, time_dec_2.5_last_1_game), time_dec_2.5_last_2_game = ifelse(is.na(time_dec_2.5_last_2_game), time_dec_2.5_last_1_game, time_dec_2.5_last_2_game), time_dec_2.5_last_4_game = ifelse(is.na(time_dec_2.5_last_4_game), time_dec_2.5_last_1_game, time_dec_2.5_last_4_game), last_game_synthetic_1 = ifelse(is.na(last_game_synthetic_1), 0, last_game_synthetic_1), last_game_synthetic_2 = ifelse(is.na(last_game_synthetic_2), last_game_synthetic_1, last_game_synthetic_2), last_game_synthetic_4 = ifelse(is.na(last_game_synthetic_4), last_game_synthetic_1, last_game_synthetic_4)) MODEL_PREDICTORS_CATEGORICAL = c("field_type", "position_group", "stadium_simp", "weather_simp", "position_group") MODEL_PREDICTORS_NUMERIC = c("temperature", "distance_per_play_last_1_game", "distance_per_play_last_2_game", "distance_per_play_last_4_game", "time_faster_5_last_1_game", "time_faster_5_last_2_game", "time_faster_5_last_4_game", "time_faster_7_last_1_game", "time_faster_7_last_2_game", "time_faster_7_last_4_game", "time_acc_1.5_last_1_game", "time_acc_1.5_last_2_game", "time_acc_1.5_last_4_game", "time_acc_2.5_last_1_game", "time_acc_2.5_last_2_game", "time_acc_2.5_last_4_game", "time_dec_1.5_last_1_game", "time_dec_1.5_last_2_game", "time_dec_1.5_last_4_game", "time_dec_2.5_last_1_game", "time_dec_2.5_last_2_game", "time_dec_2.5_last_4_game", "last_game_synthetic_1", "last_game_synthetic_2", "last_game_synthetic_4") MODEL_PREDICTORS = c(MODEL_PREDICTORS_CATEGORICAL, MODEL_PREDICTORS_NUMERIC) ``` Using this rolling window created a number of NAs, as if it is the first game of the season, we obviously will not know their time spent accelerating in their previous 4 games. Mean imputation was used here across individual players to fill these NAs, so players were given an appropriate value for themselves (as we have seen the stats can vary hugely across different players, even within the same position). I also created a 80%-20% train test split at this stage, to evaluate how good the predictions were. ```{r} PREDICTOR_FORMULA = paste(MODEL_PREDICTORS, collapse = " + ") MODEL_FORMULA = as.formula(paste0("injury_in_game ~ ", PREDICTOR_FORMULA)) set.seed(101) test_rows = sample(1:nrow(training_data), 0.2 * nrow(training_data)) orig_training_data = training_data[-test_rows, ] orig_test_data = training_data[test_rows, ] first_glm = glm(MODEL_FORMULA, data = orig_training_data, family = binomial()) summary(first_glm) ``` This is, to use a techincal term, not great. Synthetic turf flags as significant in determining what causes injuries with a positive coefficient, indicating that it increases the likelihood. The only other significant factor is the time spent accelerating at greater than 1.5 yards/sec/sec, which actually has a negative effect on injury likelihood. However interestingly, just missing the cutoff of significant was the time spent decelerating. From an amateur athletes point of view, it makes sense that deceleration would cause injury as that places a significant strain on the legs. It appears only a tenuous link in the model, however. To see if the acceleration features become more significant if other insignificant features are reduced, I removed all numeric features except temperature (which interestingly appeared significant, with warmer climates causing more injury?) and the acceleration features to assess the model fit: ```{r} acceleration_only_glm = glm(injury_in_game ~ field_type + time_acc_1.5_last_1_game + time_acc_1.5_last_2_game + time_acc_1.5_last_4_game + time_acc_2.5_last_1_game + time_acc_2.5_last_2_game + time_acc_2.5_last_4_game + time_dec_1.5_last_1_game + time_dec_1.5_last_2_game + time_dec_1.5_last_4_game + time_dec_2.5_last_1_game + time_dec_2.5_last_2_game + time_dec_2.5_last_4_game, data = orig_training_data, family = binomial()) summary(acceleration_only_glm) ``` We see no changes in the significance here, although for what its worth the AIC has improved. Even though the model does not appear to be performing super well, we will still look at how it performs on the test set: ```{r} orig_test_data = orig_test_data %>% drop_na() predictions = predict(acceleration_only_glm, orig_test_data, "response") pred_vs_actual = data.frame(predicted_prob = round(predictions, 5), actual_value = orig_test_data$injury_in_game) %>% arrange(desc(predicted_prob)) %>% mutate(num = 1:length(predictions)) ggplot(pred_vs_actual, aes(x = num, y = predicted_prob)) + geom_col(aes(fill = as.factor(actual_value))) + scale_fill_manual(values = c("black", "yellow")) + theme_minimal() + labs(title = "Predicted Injury Probabilities in Test Set", x = "Predicted Prob", fill = "Injury in Game") possible = seq(0.0001, 0.05, by = 0.0001) results = matrix(0, ncol = 2, nrow = length(possible), dimnames = list(1:length(possible), c("Sensitivity", "1 - Specificity"))) for (i in 1:length(possible)) { current = possible[i] caret_df = pred_vs_actual %>% mutate(actual_value = factor(actual_value, levels = c("0", "1")), predicted_class = ifelse(predicted_prob >= current, 1, 0), predicted_class = factor(predicted_class, levels = c("0", "1"))) conf_matrix = caret::confusionMatrix(data = caret_df$predicted_class, reference = caret_df$actual_value) results[i, 1] = conf_matrix$byClass["Sensitivity"] results[i, 2] = 1 - conf_matrix$byClass["Specificity"] } results_df = as.data.frame(results) ggplot(results_df, aes(x = `1 - Specificity`, y = Sensitivity)) + geom_point() + geom_abline(intercept = 0, slope = 1, linetype = "dashed") + labs(title = "ROC Curve") + theme_minimal() print("Confusion matrix with 0.015 chosen as threshold for positive:") # Confusion matrix using threshold of 0.015 caret_df = pred_vs_actual %>% mutate(actual_value = factor(actual_value, levels = c("0", "1")), predicted_class = ifelse(predicted_prob >= 0.015, 1, 0), predicted_class = factor(predicted_class, levels = c("0", "1"))) conf_matrix = caret::confusionMatrix(data = caret_df$predicted_class, reference = caret_df$actual_value) conf_matrix ``` The model performance here is extremely poor as seen in the confusion matrix, continually overpredicting injuries where none occurred. It did nail 14 injuries, however 5 went missed. One of the causes of this poor result may be the significant class disparity. We only have 104 games containing injuries, with even less in this training set. One potential remedy to this is to use the SMOTE algorithm (from the `DMrW` package in R) to generate synthetic samples of games with injuries. The SMOTE created an approximately 80-20 class split, which means we oversampled the minority class 10x, and it should be noted this may lead to overfitting. Regardless, the results of the model fit are below: ```{r} print("The class imbalance before using SMOTE") table(orig_training_data$injury_in_game)/nrow(orig_training_data) training_subset = orig_training_data %>% ungroup() %>% select(injury_in_game, MODEL_PREDICTORS) %>% mutate(injury_in_game = as.factor(injury_in_game), field_type = as.factor(field_type), stadium_simp = as.factor(stadium_simp), weather_simp = as.factor(weather_simp), position_group = as.factor(position_group)) %>% as.data.frame() %>% drop_na() smoted_df = SMOTE(injury_in_game ~ . , data = training_subset, perc.over = 1000, perc.under = 500) print ("The class imbalance after using SMOTE") table(smoted_df$injury_in_game)/nrow(smoted_df) smoted_model = glm(injury_in_game ~ field_type + stadium_simp + weather_simp + temperature + position_group + distance_per_play_last_1_game + distance_per_play_last_2_game + distance_per_play_last_4_game + time_faster_5_last_1_game + time_faster_5_last_2_game + time_faster_5_last_4_game + time_faster_7_last_1_game + time_faster_7_last_2_game + time_faster_7_last_4_game + time_acc_1.5_last_1_game + time_acc_1.5_last_2_game + time_acc_1.5_last_4_game + time_acc_2.5_last_1_game + time_acc_2.5_last_2_game + time_acc_2.5_last_4_game + time_dec_1.5_last_1_game + time_dec_1.5_last_2_game + time_dec_1.5_last_4_game + time_dec_2.5_last_1_game + time_dec_2.5_last_2_game + time_dec_2.5_last_4_game + last_game_synthetic_1 + last_game_synthetic_2 + last_game_synthetic_4, data = smoted_df, family = binomial()) summary(smoted_model) ``` We see some significance here! In terms of field conditions, significant factors that positively increased the risk of injury were synthetic turf, indoor fields, rain and the chance of rain. Indoor fields were likely positive due to the fact many indoor fields contain synthetic turf. For the numeric factors, at least one of each set of three factors was significant. So, the load accrued over history does have a significant impact. However, the signs within each set of three point in different directions, likely due to the highly correlated nature of these moving average stats. We can also compare on the same test set: ```{r} predictions = predict(smoted_model, orig_test_data, "response") pred_vs_actual = data.frame(predicted_prob = round(predictions, 5), actual_value = orig_test_data$injury_in_game) %>% arrange(desc(predicted_prob)) %>% mutate(num = 1:length(predictions)) ggplot(pred_vs_actual, aes(x = num, y = predicted_prob)) + geom_col(aes(fill = as.factor(actual_value))) + scale_fill_manual(values = c("black", "yellow")) + theme_minimal() + labs(title = "Predicted Injury Probabilities in Test Set", x = "Predicted Prob", fill = "Injury in Game") results = matrix(0, ncol = 2, nrow = 100, dimnames = list(1:100, c("Sensitivity", "1 - Specificity"))) possible = seq(0.01, 1, by = 0.01) for (i in 1:length(possible)) { current = possible[i] caret_df = pred_vs_actual %>% mutate(actual_value = factor(actual_value, levels = c("0", "1")), predicted_class = ifelse(predicted_prob >= current, 1, 0), predicted_class = factor(predicted_class, levels = c("0", "1"))) conf_matrix = caret::confusionMatrix(data = caret_df$predicted_class, reference = caret_df$actual_value) results[i, 1] = conf_matrix$byClass["Sensitivity"] results[i, 2] = 1 - conf_matrix$byClass["Specificity"] } results_df = as.data.frame(results) ggplot(results_df, aes(x = `1 - Specificity`, y = Sensitivity)) + geom_point() + geom_abline(intercept = 0, slope = 1, linetype = "dashed") + labs(title = "ROC Curve") + theme_minimal() print("Confusion matrix with 0.12 chosen as threshold for positive:") caret_df = pred_vs_actual %>% mutate(actual_value = factor(actual_value, levels = c("0", "1")), predicted_class = ifelse(predicted_prob >= 0.12, 1, 0), predicted_class = factor(predicted_class, levels = c("0", "1"))) conf_matrix = caret::confusionMatrix(data = caret_df$predicted_class, reference = caret_df$actual_value) conf_matrix ``` The prediction metrics see some slight improvement, however are still generally poor. I go into reasons why I think this might be the case in the next section. # Conclusion ## Limitations and Shortcomings Admittedly, there are shortcomings within this analysis: * **The data is not fully clean:** I started on this project fairly late and as a result was not able to clean the data as thoroughly as I would have liked. For instance, there are many records where the stadium type is unknown, and this could be imputed by looking at a players individual indoor/outdoor splits. Given a player plays in the same building for half their games, this value could be imputed in some way. * **Overfitting from SMOTE:** Because we had so few samples to begin with, our synthetic samples were being generated from the same 104 original samples. As a result, I think it is extremely likely we overfit to these * **Practice matters!** This analysis only takes into account games, however NFL players spend the majority of their week practicing. This practice information is vital in knowing how over-exerted a player might be, with the Rugby League paper I used as a basis incorporating both practice and game information. Having this game information would be extremely valuable in being able to predict injury. * **Classifying the start and end of plays:** In order to reduce the size of the data, I removed the tracking data from before the snap and after the play ending event. This was necessary in order to be able to work with the data, but did result in discarding a lot of information. For example, a wide receivers pre-snap motion on a jet sweep is a hard sprint that would definitely have an effect on a players injury load. However, it is not considered here. Similarly, on an incomplete pass, players do not stop running as soon as the ball hits the ground. The deceleration after the play may heighten the chance of injury, particularly if it is on the sideline or near the endzone as the player has to stop suddenly before hitting a wall. With more care, plays like this could be included but again in the interest of time, was not feasible. * **Time:** Truthfully, if you thought everything felt a bit brief, its because it was... I procrastinated and started this project too late. With more time, I would look more at the critical values of speed and acceleration (ie. which are the potential tipping points for injuries), rather than the values I arbitrarily chose. I would also try and do more around the positions, perhaps trying to scale the values for each position to some standard value so they could all be compared to each other. ## Final Remarks This report has explored the NFL 1st and Future data to try and discover links between various factors and injury, as well as differences across the different playing surfaces. I found that injuries appear to be caused by a combination of two categories: environmental factors, such as the use of synthetic grass and playing indoors, as well as load factors: repeated high speeds and high accelerations and decelerations accrued over a season. If you made it this far, thanks for reading. Geaux Tigers! # Bibliography Cummins C., Welch M., Inkster B., Cupples B., Weaving D., Jones B., King D., Murphy A. Modelling the relationships between volume, intensity and injury-risk in professional rugby league players. J. Sci. Med. Sport. 2019;22:653–660. doi: 10.1016/j.jsams.2018.11.028