{"cells":[{"metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"## Importing packages\nlibrary(tidyverse) # metapackage with lots of helpful functions\n## Import the data\ninj_data <- read.csv('../input/nfl-playing-surface-analytics/InjuryRecord.csv')\nplay_data <- read.csv('../input/nfl-playing-surface-analytics/PlayerTrackData.csv')\nplay_info <- read.csv('../input/nfl-playing-surface-analytics/PlayList.csv')\n##  Load data from the data cleaning kernel processing  to save time\nXY_offline <- read.csv('../input/nfl-injury-phases/XY.csv')\nphase_1_offline <- read.csv('../input/nfl-injury-phases/phase_1.csv')\nphase_2_offline <- read.csv('../input/nfl-injury-phases/phase_2.csv')\n## Load the data for standardizing Phase 1 and Phase 2\nstandardize_mean_p1 <- read.csv('../input/nfl-injury-data-standardization/standardize_mean_p1.csv')\nstandardize_mean_p2 <- read.csv('../input/nfl-injury-data-standardization/standardize_mean_p2.csv')\nstandardize_sd_p1 <- read.csv('../input/nfl-injury-data-standardization/standardize_sd_p1.csv')\nstandardize_sd_p2 <- read.csv('../input/nfl-injury-data-standardization/standardize_sd_p2.csv')\n## Break up the data in the healthy games and games with injuries\nall_games <- unique(play_info$GameID)\ninj_games <- unique(inj_data$GameID)\nhealthy_games <- all_games[-which(all_games %in% inj_games)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"play_info$StadiumType[play_info$Temperature == -999] <- \"Indoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Closed Dome\"] <- \"Indoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Dome, closed\"] <- \"Indoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Domed, closed\"] <- \"Indoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Indoor, Roof Closed\"] <- \"Indoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Indoors\"] <- \"Indoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Retr. Roof-Closed\"] <- \"Indoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Retr. Roof - Closed\"] <- \"Indoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Retr. Roof Closed\"] <- \"Indoor\"\n\nplay_info$StadiumType[play_info$StadiumType == \"Bowl\"] <- \"Outdoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Domed, open\"] <- \"Outdoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Domed, Open\"] <- \"Outdoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Heinz Field\"] <- \"Outdoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Indoor, Open Roof\"] <- \"Outdoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Open\"] <- \"Outdoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Oudoor\"] <- \"Outdoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Ourdoor\"] <- \"Outdoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Outddors\"] <- \"Outdoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Outdoor Retr Roof-Open\"] <- \"Outdoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Outdoors\"] <- \"Outdoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Outdor\"] <- \"Outdoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Outside\"] <- \"Outdoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Retr. Roof-Open\"] <- \"Outdoor\"\nplay_info$StadiumType[play_info$StadiumType == \"Retr. Roof - Open\"] <- \"Outdoor\"\n\nlevels(play_info$Weather)[which(grepl(\"Rain\",levels(as.factor(play_info$Weather)),fixed=TRUE))] <- \"Rain\"\nlevels(play_info$Weather)[which(grepl(\"rain\",levels(as.factor(play_info$Weather)),fixed=TRUE))] <- \"Rain\"\nlevels(play_info$Weather)[which(grepl(\"Shower\",levels(as.factor(play_info$Weather)),fixed=TRUE))] <- \"Rain\"\n\nlevels(play_info$Weather)[which(grepl(\"Snow\",levels(as.factor(play_info$Weather)),fixed=TRUE))] <- \"Snow\"\nlevels(play_info$Weather)[which(grepl(\"snow\",levels(as.factor(play_info$Weather)),fixed=TRUE))] <- \"Snow\"\n\nplay_info$Temperature[play_info$Temperature==-999] <- mean(play_info$Temperature[play_info$Temperature!=-999])\n\nlevels(play_info$PlayType)[which(grepl(\"Punt\",levels(as.factor(play_info$PlayType)),fixed=TRUE))] <- \"Return\"\nlevels(play_info$PlayType)[which(grepl(\"Kickoff\",levels(as.factor(play_info$PlayType)),fixed=TRUE))] <- \"Return\"\n\nlevels(play_info$PlayType)[which(grepl(\"Extra Point\",levels(as.factor(play_info$PlayType)),fixed=TRUE))] <- \"Kick\"\nlevels(play_info$PlayType)[which(grepl(\"Field Goal\",levels(as.factor(play_info$PlayType)),fixed=TRUE))] <- \"Kick\"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Movement and Injury in the NFL\n\n# Introduction\n\nThis notebook contains my analysis for the 2019 NFL 1st and Future Analytics comptetion. \n\n### Objective\n\n> Your challenge is to characterize any differences in player movement between the playing surfaces and identify specific scenarios (e.g., field surface, weather, position, play type, etc.) that interact with player movement to present an elevated risk of injury.\n\n## Summary of Results\n\nIn the following results a few conclusions will be argued:\n\n1. There are specific movement patterns which are related to increased injury probability, and there are indications that are also interactions between game-specific variables and some movement patterns which together lead to an even higher increase in injury probability\n2. There are also game-specific variables (field surface, weather, and play type) which also are related to increased injury probability\n3. There are slight differences in movement between synthetic and natural playing surfaces, however these differences have less of an affect than the field surface itself\n\nIn order to arrive at these conclusions, a generalized linear model to predict the injury probability of any given play taking into account game-specific and movement conditions (as well as interactions between all possible combinations of these variables) is created. The development of the model will be shown below. The final model predicted below average injury probability for 74% of plays with no injury, and above average injury probability for 75% of plays with an injury.\n\n\n"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plot(c(0,2),c(0,0),type=\"l\",lty=2,col=1,ylim=c(-100,100),ylab=\"Percent Predicted Injury Probability Above/Below Average\",axes=FALSE,xlab=\"\")\naxis(1,labels=c(\"No Injury\",\"Injury\"), at=c(.5,1.5),las=1,cex=.25)\naxis(2)\npolygon(c(.1,.1,.9,.9),c(0,-74,-74,0),col=rgb(0,1,0,.75),border=FALSE)\npolygon(c(.1,.1,.9,.9),c(0,26,26,0),col=rgb(0,1,0,.4),border=FALSE)\npolygon(c(1.1,1.1,1.9,1.9),c(0,75,75,0),col=rgb(1,0,0,.75),border=FALSE)\npolygon(c(1.1,1.1,1.9,1.9),c(0,-25,-25,0),col=rgb(1,0,0,.4),border=FALSE)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From this model we can summarize some of our results.\n\n### 1. Movement Patterns which Increase Injury Probability\n\nTwo main movement patterns emerged in the analysis which increase the probability of injury: quick acceleration and immediate deceleration, and just quick acceleration. Further, an interaction between the quick acceleration movement pattern and the synthetic field combined to increase the probability of injury even more. In order to demonstrate the effect of these movement patterns, one play is chosen which shows how that movement type affects the predicted injury probability by the final model. More infromation on these plays will be discussed later.\n\nThe first movement type is a quick acceleration and immediate quick deceleration. Predicted injury probability on one of the plays on which a player was injured who exhibited this pattern is shown below."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plot(c(-1,3),c(0.000288951,0.0002889591),type=\"l\",lty=2,ylim=c(0,2e-3),axes=FALSE,ylab=\"Predicted Injury Probability\",xlab=\"\")\nbarplot(c(5.854532e-04,1.598115e-03),names=c(\"Without Pattern\",\"With Pattern\"),add=TRUE,col=c(3,2))\ntext(-.45,.00032,\"Average Injury Probability\",cex=.7)\ntext(1.9,1.675e-3,\"2.7x Increase\",col=2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This play had a predicted injury probability well above the average injury probability across all plays in the dataset, and the influence of factors pertaining to the first movement pattern lead to 2.7 times increase in injury probability.\n\nThe second movement type is a quick acceleration. This movement type is also notable because the model suggests that when this movement type happens on a synthetic playing surface, the effect of the movment type to increase injury probability is greatly increased."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plot(c(-1,3),c(0.000288951,0.0002889591),type=\"l\",lty=2,ylim=c(0,1.5e-2),axes=FALSE,ylab=\"Predicted Injury Probability\",xlab=\"\")\nbarplot(c(3.764652e-03,1.239269e-02),names=c(\"Without Pattern\",\"With Pattern\"),add=TRUE,col=c(3,2))\ntext(-.45,.0005,\"Average Injury Probability\",cex=.7)\ntext(1.9,1.3e-2,\"3.3x Increase\",col=2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This play showed an even higher predicted injury probability than the previous play. The combination of the play happening on a synthetic playing surface and the second movement type lead to a 3.3 times increase in the predicted injury probability.\n\n### 2. Game-specific Variables which Increase Injury Probability\n\nA variety of game-specific varaibles and interactions between game-specific variables and movement patterns increase the probability of injury according to the model. A sample of these how these effects compare to the effect of a synthetic surface are shown below."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"ptemp <- 0.3191/0.6076*100\nprain <- 0.282/0.6076*100\npreturn <- 1.0632/.6076*100\npfields22 <- 0.3346/0.6076*100\nprainmaxdv <- 1.2311/.6076*100\n\nbarplot(c(100,pfields22,ptemp,prain,prainmaxdv,preturn),names=c(\"Synthetic\",\"Synth:s22\",\"Temp\",\"Rain\",\"Rain:maxdv\",\"Return\"),col=c(2,2,2,2,2,2),ylab=\"Percentage of Effect of Synthetic Turf\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that not only playing on a synthetic surface increases injury probability, but also playing at higher temperatures, playing in the rain, or playing on a kickoff or punt return. In fact the maximum effect of a single game-specific variable is playing on a return play. \n\nThere are also intractions, such as a large increase in injury probability if it is raining and the player is rapidly changing movement speed and/or direction. We can also see that according to the model, the second movement type discussed above increases the effect of the synthetic field by an extra 50%.\n\n### Differences in Movement Across Playing Surfaces\n\nIt would be theoretically possible that an increase in injury occurence on syntheitic playing surfaces is a secondary effect from players moving differently on the different surface. However, for the main movement types which affect the injury probability, the differences between player movement on each playing surface would account for less than 5% of the effect of the synthetic field itself. This shows that there is indeed an increase in injury probability simply by playing a play on a synthetic field."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"p32 <- -1.78153\np22 <- -0.7925787\npmaxdv <- 0.6512364\npmaxs <- 4.72725\npv3 <- 2.295924\n\nbarplot(c(100,p22,p32,pmaxdv,pmaxs,pv3),names=c(\"Synthetic\",\"s22\",\"s32\",\"maxdv\",\"maxs\",\"v3\"),col=c(2,3,3,2,2,2),ylab=\"Percentage of Effect of Synthetic Turf\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n\n# Analysis\n\nThe main analysis in this notebook revolves around using generalized linear models to predict the probability of a player getting injured on any specific play. This analysis is done in two phases:\n\nPhase 1 includes all game-specific variables (e.g., field surface, weather, position, play type, etc.) along with autocorrelations of movement variables including the speed, movement direction, and the change in movement vector (taking into account both speed and direction simultaneously). The autocorrelation function (ACF) is used in Phase 1 because a play could last any length of time, and the exact time of the injury within that timeframe is unknown. \n\nIn Phase 2, lessons learned from Phase 1 are leveraged. Subsets of 3 seconds of each injury play are used to find the time period with the highest injury probability as predicted by Phase 1. The time histories of the movement variables during this time period then become factors for the general linear model in the second phase.\n\n### Singular Value Decomposition for Pattern Recognition\n\nIn both phases, the Singular Value Decomposition (SVD) is used to turn movement patterns during known injury plays into a smaller number of templates which can be used in our model fitting. \n\nThe first step, is to create a matrix, A, which has a movement pattern vector for each of the 76 known injury plays as a column. In Phase 1, for example, we create a matrix each for the speed, direction change, and movement vector change where each column is 3 seconds of the ACF. A graphical example of this matrix is shown below for the speed variable in Phase 1 where 10 of the ACF vectors for the first 10 injury plays are shown as lines."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"## Find all the plays where we know there is an injury\ninj_plays <- unique(inj_data$PlayKey[inj_data$PlayKey!=\"\"])\nlag_max <- 30;\nP <- length(inj_plays);\nAv <- matrix(0,lag_max+1,P); As <- matrix(0,lag_max+1,P); Adir <- matrix(0,lag_max+1,P)\nfor(k in 1:P){\n  sing_play_data <- play_data[which(play_data$PlayKey == as.character(inj_plays[k])),]\n  N <- nrow(sing_play_data)\n  xv <- sing_play_data$s*sin(sing_play_data$dir*pi/180)\n  yv <- sing_play_data$s*cos(sing_play_data$dir*pi/180)\n  dv <- sqrt((xv[2:N]-xv[1:(N-1)])^2+(yv[2:N]-yv[1:(N-1)])^2)\n  ds <- sing_play_data$s[2:(N)]-sing_play_data$s[1:(N-1)]\n  dang <- yv/xv; dang[is.nan(dang)] <- 0\n  ddir <- atan(dang)\n  rv<-acf(dv,lag.max=lag_max,type=\"covariance\",plot=FALSE)\n  rs<-acf(ds,lag.max=lag_max,type=\"covariance\",plot=FALSE)\n  rdir<-acf(ddir,lag.max=lag_max,type=\"covariance\",plot=FALSE)\n  Av[,k] <- rv$acf\n  As[,k] <- rs$acf\n  Adir[,k] <- rdir$acf\n}\n\nplot(.5+As[,1]/max(As[,1])*2/3,seq(from=0,to=3,by=.1),type=\"l\",xlim=c(0,11),xlab=\"Injury Plays\",ylab=\"Autocorrelation Time (s)\",xaxt=\"n\")\nfor(k in 1:10){\n  lines(c(k,k)-.5,c(-1,4),lty=2,col=\"gray\")\n}\nfor(k in 2:10){\n  lines(k-.5+As[,k]/max(As[,k])*2/3,seq(from=0,to=3,by=.1),col=k)\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"If we call the matrix created with the movement indicators from the injury plays $A$, then with the SVD we break this matrix into three matrices\n\n$$ A = U \\Sigma V^* $$\n\nwhere $ U $ is a matrix of the singular vectors. The singular vectors represent the vectors which explain the most variance from the original matrix $A$. Therefore, we take first five singular vectors for this matrix then represent the dominant patterns from these matrices: meaning that the SVD allows us to pull out patterns which happen are common to all injury plays. The singular vectors for each movement type in Phase 1 are shown below."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":false},"cell_type":"code","source":"Sv <- svd(Av)\nSs <- svd(As)\nSdir <- svd(Adir)\n#par(mfrow=c(3,1))\nplot(Ss$u[,1],type=\"l\",ylim=c(-.5,.5),ylab=\"Speed ACF Patterns\",xlab=\"Lag Index\")\nfor(v in 2:5){lines(Ss$u[,v],col=v)}; legend(\"topleft\",legend=c(\"1\",\"2\",\"3\",\"4\",\"5\"),horiz=TRUE,lty=1,col=1:5,cex=1)\nplot(Sv$u[,1],type=\"l\",ylim=c(-.5,.5),ylab=\"Vector ACF Patterns\",xlab=\"Lag Index\")\nfor(v in 2:5){lines(Sv$u[,v],col=v)}; legend(\"topleft\",legend=c(\"1\",\"2\",\"3\",\"4\",\"5\"),horiz=TRUE,lty=1,col=1:5,cex=1)\nplot(Sdir$u[,1],type=\"l\",ylim=c(-.5,.5),ylab=\"Direction ACF Patterns\",xlab=\"Lag Index\")\nfor(v in 2:5){lines(Sdir$u[,v],col=v)}; legend(\"topleft\",legend=c(\"1\",\"2\",\"3\",\"4\",\"5\"),horiz=TRUE,lty=1,col=1:5,cex=1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Data and Model Formulation\n\nA total of 46 variables which may affect the injury probability of an NFL player are used in the final model of this analysis. These effects are in three main categories: game-specific variables, Phase 1, and Phase 2 variables. All the variables for each type are:\n\n### Game Specific Variables\n\n| Variable | Description |\n| -- | -- |\n| Field | Field type, either \"Synthetic\" or \"Natural\" |\n| Return | Boolean for if the play type is a kickoff or punt return |\n| Kick | Boolean for if the play type is a field goal or extra point |\n| Rush | Boolean for if the play type is a rushing play |\n| Pass | Boolean for if the play type is a passing play |\n| Snow | Boolean for if the weather report indicated snow |\n| Rain | Boolean for if the weather report indicated rain |\n| Precipitation | Boolean for if the weather report indicated snow and/or rain |\n| Indoor | Boolean for if the game is played indoors (or with a retractable roof closed) |\n| Outdoor | Boolean for if the game is played outdoors (or with a retractable roof open) |\n| Defense | Boolean for if the player's position is on the defense |\n| Offense | Boolean for if the player's position is on the offense |\n| Line | Boolean for if the player's position is on the offensive or defensive line |\n| Temperature | The temperature recorded for the game (temperature for indoors games is set to the mean of all other games) |\n\n### Phase 1 Movement Variables\n\n| Variable | Description |\n| -- | -- |\n| maxs | The maximum speed of the player during the play |\n| maxdv | The maximum movement vector change during the play |\n| v1-v5 | Correlation with the first five singular vectors for movement vector changes during injury plays |\n| s1-s5 | Correlation with the first five singular vectors for speed changes during injury plays |\n| d1-d5 | Correlation with the first five singular vectors for direction changes during injury plays |\n\n### Phase 2 Movement Variables\n\n| Variable | Description |\n| -- | -- |\n| v12-v52 | Inner product with the five movement vector change time history templates |\n| s12-s52 | Inner product with the five speed time history templates |\n| d12-d52 | Inner product with the five direction change time history templates |\n\n### Differences between Phase 1 and Phase 2\n\nThere are two main differences between the Phase 1 and Phase 2 of the analysis. Phase 1 only takes into account the autocorrelation of the movement variables since the time of injury is unknown. This means that only patterns in the overall history can be accounted for (are there quick changes, are there periodic changes), while Phase 2 can take into account specific movement patterns (did the player accelerate, accelerate then quickly decelerate etc.). The second difference is that since there are no time history templates available in Phase 1, the magnitude of the autcorrelation patterns is not accounted for, while in Phase 2 the magnitude can be accounted for (allowing a difference betweeen acceleration to a slow speed and acceleration to a high top speed).\n\n### Model Formulation\n\nThe model which will be used is a binomial generalized linear model to predict the injury probability of each play. The model can be expressed as:\n\n$$ p_i = e^{x_i \\beta } $$\n\nwhere $p_i$ is the injury probaility of play $i$, $x_i$ are the relevant play factors (such as the playing surface type, or the magnitude of movement types) for the play, and $\\beta$ are the model coefficients for each factor. When fitting the model, interactions between all combinations of factors are considered, and the important factors are selected by step-wise fitting with BIC as the decision criteria. This ensures that the most important factors are in our final model while also ensuring model parsimony."},{"metadata":{},"cell_type":"markdown","source":"# Phase 1\n\nIn phase 1, the autocorrelation singular vectors discussed above are used in the model fitting. Following the step-wise fitting, we arrive at the phase 1 model below."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"XY1 <- cbind(XY_offline,phase_1_offline)\n## In order to save time, we initialize our step-wise fitting at the correct model.\n# null1 <- suppressWarnings(glm(inj~1,data=XY1,family=\"binomial\"))\n# full1 <- suppressWarnings(glm(inj~.^2,data=XY1[1:1000,],family=\"binomial\"))\n# fit1 <- suppressWarnings(step(null1,scope=formula(full1),direction=\"both\",trace=1,k=2))\nfit1 <- suppressWarnings(glm(inj~1+maxs+return+d2+v4+temperature+field+kick+v3+s4+v1+rain+maxs*s3+maxs*v3+maxs*v1+return*v4+maxs*rain+v3*rain+v4*s3+temperature*v1,data=XY1,family=\"binomial\"))\nsummary(fit1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see which variables were selected as important by the step-wise fitting. Both movement variables and game-specific variables are chosen as important, along with some interactions between these. The cumulative distributions of the predicted injury probability for both injury and non-injury plays is shown below. Even in phase 1, the model is capable of predicting above average injury probability for 68% of injury plays and below average injury probability for 74% of non-injury plays."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"pred1<-predict(fit1,newdata=XY1,type=\"response\")\n\nplot(c(76/nrow(XY1),76/nrow(XY1)),c(0,1),type=\"l\",lty=2,xlim=c(-.000,0.003),ylim=c(0,1.05),main=\"\",xlab=\"Predicted Injury Probability\",ylab=\"Cumulative Distribution\")\n\np0 <- sum(pred1[XY1$inj==0]<76/nrow(XY1))/nrow(XY1)\np1 <- sum(pred1[XY1$inj==1]<76/nrow(XY1))/76\n\npolygon(c(-1,-1,76/nrow(XY1),76/nrow(XY1)),c(0,p0,p0,0),col=rgb(0,1,0,.1),border=FALSE)\npolygon(c(76/nrow(XY1),76/nrow(XY1),.004,004),c(p1,1,1,p1),col=rgb(1,0,0,.1),border=FALSE)\nlines(ecdf(pred1[XY1$inj==0]),col=3)\nlines(ecdf(pred1[XY1$inj==1]),col=2,cex=0)\n\nlines(c(-1,76/nrow(XY1)),c(p0,p0),col=3,lty=2)\n\nlines(c(76/nrow(XY1),1),c(p1,p1),col=2,lty=2)\ntext(76/nrow(XY1)+1/10000,.625,\"68% ABOVE Average\",col=2,adj=0,cex=.7,srt=40)\ntext(76/nrow(XY1)-1.3/10000,.65,\"74% BELOW Average\",col=3,adj=1,cex=.7,srt=80)\ntext(.00025,.9,\"Average Injury Probability\",cex=.7,adj=.5,srt=90)\nlegend(.00038,.15,legend=c(\"No Injury Play\",\"Injury Play\"),col=c(3,2),lty=1,cex=.7,box.lty=0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Phase 2\n\nIn order to extend from the lessons learned in phase 1 and find specific time-history templates for movement variables which can be more easily interpreted we turn to phase 2 of our analysis. First, we take each injury play and consider 3 second subsections of the play and use the phase 1 model to predict what 3 second window of the play was the most likley to cause injury (this is shown below as the shaded area). The 3 second time history of the speed, direction change, and movement vector change during the most likely injury time were taken for each play and the same SVD process as before was used to get movement templates for phase 2."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"s_patterns <- matrix(0,(lag_max+1),length(inj_plays))\ndv_patterns <- matrix(0,(lag_max+1),length(inj_plays))\ndir_patterns <- matrix(0,(lag_max+1),length(inj_plays))\nfor(j in 1:length(inj_plays)){\nsing_play_data <- play_data[which(play_data$PlayKey == as.character(inj_plays[j])),]\nN <- nrow(sing_play_data)\n      #xv <- sing_play_data$x[2:(N)]-sing_play_data$x[1:(N-1)]\n      #yv <- sing_play_data$y[2:(N)]-sing_play_data$y[1:(N-1)]\n      xv <- sing_play_data$s*sin(sing_play_data$dir*pi/180)\n      yv <- sing_play_data$s*cos(sing_play_data$dir*pi/180)\n      #n <- N-1\n      n <- N\n      dv <- sqrt((xv[2:n]-xv[1:(n-1)])^2+(yv[2:n]-yv[1:(n-1)])^2)\n      ds <- sing_play_data$s[2:(N)]-sing_play_data$s[1:(N-1)]\n      ddir <- abs(sin((sing_play_data$dir[2:N]-sing_play_data$dir[1:(N-1)])*pi/180))\n      \nv4 <- rep(0,N-32); v3 <- rep(0,N-32); d2 <- rep(0,N-32); maxs <- rep(0,N-32); maxdv <- rep(0,N-32)\nfor(k in 1:(N-32)){\n  maxs[k] <- (max(sing_play_data$s[(1+k):(1+lag_max+k)]) - 4.848755)/2.020227\n  maxdv[k] <- (max(dv[k:(lag_max+k)]) - 0.5017628)/0.4935381\n  rv <- acf(dv[k:(lag_max+k)],lag.max=lag_max,type=\"covariance\",plot=FALSE)\n  rd <- acf(ddir[k:(lag_max+k)],lag.max=lag_max,type=\"covariance\",plot=FALSE)\n  if(sd(dv[k:(lag_max+k)])>0){\n  v3[k] <- (cor(rv$acf,Sv$u[,3]) - 0.001636858)/0.1722944\n  v4[k] <- (cor(rv$acf,Sv$u[,4]) - 0.03662055)/0.1869571}\n  if(sd(ddir[k:(lag_max+k)])>0){d2[k] <- (cor(rd$acf,Sdir$u[,2]) + 0.8219374)/0.09065627}\n  \n}\nprob <- exp(-8.7375 + 0.8808*maxs + 0.3094*d2 - 0.9219*maxdv -0.3066*v4 + 0.3673*v3 +0.4158*v3*maxdv -0.3917*v4*maxdv)\nif(j == 75){\n    plot(prob/max(prob),type=\"l\",xlab=\"Sample Number\",ylab=\"Normalized Value\")\n    lines(sing_play_data$s/max(sing_play_data$s),col=2,lty=2)\n    lines(ddir/max(ddir),col=3,lty=2)\n    lines(dv/max(dv),col=4,lty=2)\n    ind1 <- which.max(prob)\n    polygon(c(ind1,ind1,ind1+lag_max,ind1+lag_max),c(0,1,1,0),col=rgb(0,0,0,.1),border=FALSE)\n    legend(\"topright\",legend=c(\"Injury Probability\",\"Speed\",\"Direction Change\",\"Movement Vector                               .\"),col=c(1,2,3,4),lty=c(1,2,2,2),cex=0.8)\n}\n    \nplot_inds <- which.max(prob):(which.max(prob)+lag_max)\n\ns_patterns[,j] <- sing_play_data$s[plot_inds+1]\ndv_patterns[,j] <- dv[plot_inds]\ndir_patterns[,j] <- ddir[plot_inds]\n\n#if(j==1){plot(sing_play_data$s[plot_inds]/max(sing_play_data$s),col=2,type=\"l\",ylim=c(0,1))}\n#if(j>1){lines(sing_play_data$s[plot_inds]/max(sing_play_data$s),col=2)}\n#lines(ddir[plot_inds]/max(ddir),col=3)\n#lines(dv[plot_inds]/max(dv),col=4)\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The first five singular vectors for use in phase 2 are shown below:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"Ss2 <- svd(s_patterns)\nSv2 <- svd(dv_patterns)\nSdir2 <- svd(dir_patterns)\n#par(mfrow=c(3,1))\nplot(Ss2$u[,1],type=\"l\",ylim=c(-.5,.5),ylab=\"Speed Time History Patterns\",xlab=\"Lag Index\")\nfor(v in 2:5){lines(Ss2$u[,v],col=v)}; legend(\"topleft\",legend=c(\"1\",\"2\",\"3\",\"4\",\"5\"),horiz=TRUE,lty=1,col=1:5,cex=1)\n\nplot(Sv2$u[,1],type=\"l\",ylim=c(-.5,.5),ylab=\"Vector Time History Patterns\",xlab=\"Lag Index\")\nfor(v in 2:5){lines(Sv2$u[,v],col=v)}; legend(\"topleft\",legend=c(\"1\",\"2\",\"3\",\"4\",\"5\"),horiz=TRUE,lty=1,col=1:5,cex=1)\n\nplot(Sdir2$u[,1],type=\"l\",ylim=c(-.5,.5),ylab=\"Direction Time History Patterns\",xlab=\"Lag Index\")\nfor(v in 2:5){lines(Sdir2$u[,v],col=v)}; legend(\"topleft\",legend=c(\"1\",\"2\",\"3\",\"4\",\"5\"),horiz=TRUE,lty=1,col=1:5,cex=1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we can fit our Phase 2 model with all variables from Phase 1, and also the templates from Phase 2. For each play, the 3 second time period with the highest inner product with each singular vector tells us the magnitude with which the play contains that template. The final phase 2 model is then fitted in the same step-wise fashion:"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"XY2 <- cbind(XY_offline,phase_1_offline,phase_2_offline)\n## In order to save time, we initialize our step-wise fitting at the correct model\n# null2 <- suppressWarnings(glm(inj~1,data=XY2[1:1000,],family=\"binomial\"))\n# full2 <- suppressWarnings(glm(inj~.^2,data=XY2[1:1000,],family=\"binomial\"))\n# fit2 <- suppressWarnings(step(null2,scope=formula(full2),direction=\"both\",trace=1,k=2))\nfit2 <- suppressWarnings(glm(inj~1+s32+d2+maxdv+temperature+field+kick+v3+v4+rain+return+s2+d32+s22+v1+maxs+s32*rain+s32*v3+s32*v4+d2*d32+maxdv*s22+field*s22+v4*d32+maxdv*v1+maxdv*v3+maxdv*rain+maxdv*s2+s2*v1,data=XY2,family=\"binomial\"))\nsummary(fit2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that Phase 2 has added useful information since the Phase 2 model is capable of predicting above average injury probability for 75% of injury plays and below average injury probability for 74% of non-injury plays."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"pred2<-predict(fit2,newdata=XY2,type=\"response\")\n\nplot(c(76/nrow(XY2),76/nrow(XY2)),c(0,1),type=\"l\",lty=2,xlim=c(-.000,0.003),ylim=c(0,1.05),main=\"\",xlab=\"Predicted Injury Probability\",ylab=\"Cumulative Distribution\")\n\np0 <- sum(pred2[XY2$inj==0]<76/nrow(XY2))/nrow(XY2)\np1 <- sum(pred2[XY2$inj==1]<76/nrow(XY2))/76\n\npolygon(c(-1,-1,76/nrow(XY2),76/nrow(XY2)),c(0,p0,p0,0),col=rgb(0,1,0,.1),border=FALSE)\npolygon(c(76/nrow(XY2),76/nrow(XY2),.004,004),c(p1,1,1,p1),col=rgb(1,0,0,.1),border=FALSE)\nlines(ecdf(pred2[XY2$inj==0]),col=3)\nlines(ecdf(pred2[XY2$inj==1]),col=2,cex=0)\n\nlines(c(-1,76/nrow(XY2)),c(p0,p0),col=3,lty=2)\n\nlines(c(76/nrow(XY2),1),c(p1,p1),col=2,lty=2)\ntext(76/nrow(XY2)+1/10000,.5,\"75% ABOVE Average\",col=2,adj=0,cex=.7,srt=40)\ntext(76/nrow(XY2)-1.3/10000,.65,\"74% BELOW Average\",col=3,adj=1,cex=.7,srt=80)\ntext(.00025,.9,\"Average Injury Probability\",cex=.7,adj=.5,srt=90)\nlegend(.00038,.15,legend=c(\"No Injury Play\",\"Injury Play\"),col=c(3,2),lty=1,cex=.7,box.lty=0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Movement Patterns which Increase Injury Probability\n\nWith the results from our Phase 2 model, we can begin to make inferences about what movement patterns increase injury probability.\n\n### Quick Acceleration and Deceleration\n\nThe first movement pattern which the model suggests increases injury probability is \"s32\", which is a speed template meaning that the player had a quick acceleration and immediately following quick deceleration."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plot(Ss2$u[,3],type=\"l\",ylim=c(-.5,.5),ylab=\"Speed Time History Pattern\",xlab=\"Lag Index\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"An example of a play which shows a significant influence of this movement pattern on the pedicted injury probability is shown below."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#play_key_all <- unique(play_data$PlayKey)\n#play_key_all <- c(as.character(play_key_all[sample(1:length(play_key_all),1000,replace=FALSE)]),as.character(play_keys))\nfor(k in 27){\nplay_num <- k\ncur_info <- play_info[which(play_info$PlayKey==as.character(inj_plays[k])),]\nsing_play_data <- play_data[which(play_data$PlayKey == as.character(inj_plays[k])),]\nN <- nrow(sing_play_data)  \nN_plays <- N - lag_max - 1  \ninj <- rep(1,N_plays); \nv1 <- rep(0,N_plays); v2 <- rep(0,N_plays); v3 <- rep(0,N_plays); v4 <- rep(0,N_plays); v5 <- rep(0,N_plays);\ns1 <- rep(0,N_plays); s2 <- rep(0,N_plays); s3 <- rep(0,N_plays); s4 <- rep(0,N_plays); s5 <- rep(0,N_plays);\nd1 <- rep(0,N_plays); d2 <- rep(0,N_plays); d3 <- rep(0,N_plays); d4 <- rep(0,N_plays); d5 <- rep(0,N_plays);\n#\nv12 <- rep(0,N_plays); v22 <- rep(0,N_plays); v32 <- rep(0,N_plays); v42 <- rep(0,N_plays); v52 <- rep(0,N_plays);\ns12 <- rep(0,N_plays); s22 <- rep(0,N_plays); s32 <- rep(0,N_plays); s42 <- rep(0,N_plays); s52 <- rep(0,N_plays);\nd12 <- rep(0,N_plays); d22 <- rep(0,N_plays); d32 <- rep(0,N_plays); d42 <- rep(0,N_plays); d52 <- rep(0,N_plays);\n#\nmaxs <- rep(0,N_plays); maxdv <- rep(0,N_plays); field <- rep(\"\",N_plays);\n# Play Types\nreturn <- rep(0,N_plays); kick <- rep(0,N_plays); rush <- rep(0,N_plays); pass <- rep(0,N_plays);\n# Weather\nsnow <- rep(0,N_plays); rain <- rep(0,N_plays); precipitation <- rep(0,N_plays);\n# Stadium\nindoor <- rep(0,N_plays); outdoor <- rep(0,N_plays);\n# Temp\ntemperature <- rep(0,N_plays)\n# Defense or Offense\ndefense <- rep(0,N_plays); offense <- rep(0,N_plays); line <- rep(0,N_plays)\n\n  #xv <- sing_play_data$x[2:(N)]-sing_play_data$x[1:(N-1)]\n  #yv <- sing_play_data$y[2:(N)]-sing_play_data$y[1:(N-1)]\n  xv <- sing_play_data$s*sin(sing_play_data$dir*pi/180)\n  yv <- sing_play_data$s*cos(sing_play_data$dir*pi/180)\n  #n <- N-1\n  n <- N\n  dv <- sqrt((xv[2:n]-xv[1:(n-1)])^2+(yv[2:n]-yv[1:(n-1)])^2)\n  ds <- sing_play_data$s[2:(N)]-sing_play_data$s[1:(N-1)]\n  ddir <- abs(sin((sing_play_data$dir[2:N]-sing_play_data$dir[1:(N-1)])*pi/180))\n  \n  s<- sing_play_data$s\n  #\n  \n  for(j in 1:N_plays){\n    rv<-acf(dv,lag.max=lag_max,type=\"covariance\",plot=FALSE)\n    rs<-acf(ds,lag.max=lag_max,type=\"covariance\",plot=FALSE)\n    rdir<-acf(ddir,lag.max=lag_max,type=\"covariance\",plot=FALSE)\n    \n    v12[j] <- dv[j:(j+lag_max)] %*% Sv2$u[,1]\n    v22[j] <- dv[j:(j+lag_max)] %*% Sv2$u[,2]\n    v32[j] <- dv[j:(j+lag_max)] %*% Sv2$u[,3]\n    v42[j] <- dv[j:(j+lag_max)] %*% Sv2$u[,4]\n    v52[j] <- dv[j:(j+lag_max)] %*% Sv2$u[,5]\n    #\n    s12[j] <- s[j:(j+lag_max)] %*% Ss2$u[,1]\n    s22[j] <- s[j:(j+lag_max)] %*% Ss2$u[,2]\n    s32[j] <- s[j:(j+lag_max)] %*% Ss2$u[,3]\n    s42[j] <- s[j:(j+lag_max)] %*% Ss2$u[,4]\n    s52[j] <- s[j:(j+lag_max)] %*% Ss2$u[,5]\n    #\n    d12[j] <- ddir[j:(j+lag_max)] %*% Sdir2$u[,1]\n    d22[j] <- ddir[j:(j+lag_max)] %*% Sdir2$u[,2]\n    d32[j] <- ddir[j:(j+lag_max)] %*% Sdir2$u[,3]\n    d42[j] <- ddir[j:(j+lag_max)] %*% Sdir2$u[,4]\n    d52[j] <- ddir[j:(j+lag_max)] %*% Sdir2$u[,5]\n  \n  #\n  v1[j] <- cor(rv$acf,Sv$u[,1])\n  v2[j] <- cor(rv$acf,Sv$u[,2])\n  v3[j] <- cor(rv$acf,Sv$u[,3])\n  v4[j] <- cor(rv$acf,Sv$u[,4])\n  v5[j] <- cor(rv$acf,Sv$u[,5])\n  #\n  s1[j] <- cor(rs$acf,Ss$u[,1])\n  s2[j] <- cor(rs$acf,Ss$u[,2])\n  s3[j] <- cor(rs$acf,Ss$u[,3])\n  s4[j] <- cor(rs$acf,Ss$u[,4])\n  s5[j] <- cor(rs$acf,Ss$u[,5])\n  #\n  d1[j] <- cor(rdir$acf,Sdir$u[,1])\n  d2[j] <- cor(rdir$acf,Sdir$u[,2])\n  d3[j] <- cor(rdir$acf,Sdir$u[,3])\n  d4[j] <- cor(rdir$acf,Sdir$u[,4])\n  d5[j] <- cor(rdir$acf,Sdir$u[,5])\n  maxs[j] <- max(sing_play_data$s[j:(j+lag_max)])\n  maxdv[j] <- max(dv[j:(j+lag_max)])\n  \n  field[j] <- as.character(cur_info$FieldType)\n  # Play Types\n  if(cur_info$PlayType == \"Return\"){return[j] <- 1}\n  if(cur_info$PlayType == \"Kick\"){kick[j] <- 1}\n  if(cur_info$PlayType == \"Rush\"){rush[j] <- 1}\n  if(cur_info$PlayType == \"Pass\"){pass[j] <- 1}\n  # Weather Types\n  if(cur_info$Weather == \"Snow\"){snow[j] <- 1; precipitation[j] <- 1}\n  if(cur_info$Weather == \"Rain\"){rain[j] <- 1; precipitation[j] <- 1}\n  # Stadium Types\n  if(cur_info$StadiumType == \"Indoor\"){indoor[j] <- 1}\n  if(cur_info$StadiumType == \"Outdoor\"){outdoor[j] <- 1}\n  # Temperature\n  temperature[j] <- cur_info$Temperature\n  # Position Types\n  cur_pos <- cur_info$PositionGroup\n  if(cur_pos==\"DB\"|cur_pos==\"DL\"|cur_pos==\"LB\"){defense[j] <- 1}\n  if(cur_pos==\"RB\"|cur_pos==\"OL\"|cur_pos==\"QB\"|cur_pos==\"TE\"|cur_pos==\"WR\"){offense[j] <- 1}\n  if(cur_pos==\"DL\"|cur_pos==\"OL\"){line[j] <- 1}\n  }\n\nplay_XY <- data.frame(inj=inj,field=as.factor(field),return=return,kick=kick,rush=rush,pass=pass,snow=snow,rain=rain,precipitation=precipitation,indoor=indoor,outdoor=outdoor,defense=defense,offense=offense,line=line)\nplay_phase_1 <- data.frame(temperature=temperature,maxs=maxs,maxdv=maxdv,v1=v1,v2=v2,v3=v3,v4=v4,v5=v5,s1=s1,s2=s2,s3=s3,s4=s4,s5=s5,d1=d1,d2=d2,d3=d3,d4=d4,d5=d5)\nplay_phase_2 <- data.frame(v12=v12,v22=v22,v32=v32,v42=v42,v52=v52,s12=s12,s22=s22,s32=s32,s42=s42,s52=s52,d12=d12,d22=d22,d32=d32,d42=d42,d52=d52)\n\nfor(k in 1:ncol(play_phase_1)){\n  play_phase_1[,k] <- (play_phase_1[,k] - standardize_mean_p1[k,1])/standardize_sd_p1[k,1]\n}\n\nfor(k in 1:ncol(play_phase_2)){\n  play_phase_2[,k] <- (play_phase_2[,k] - standardize_mean_p2[k,1])/standardize_sd_p2[k,1]\n}\n\nplay_XY2 <- cbind(play_XY,play_phase_1,play_phase_2)\n\nprob <- predict(fit2,play_XY2,type=\"response\")\n\nplot(prob/max(prob),type=\"l\",xlab=\"Sample Number\",ylab=\"Normalized Value\",ylim=c(0,1))\n    lines(sing_play_data$s/max(sing_play_data$s),col=2,lty=2)\n    lines(ddir/max(ddir),col=3,lty=2)\n    lines(dv/max(dv),col=4,lty=2)\n    ind1 <- which.max(prob)\n    polygon(c(ind1,ind1,ind1+lag_max,ind1+lag_max),c(0,1,1,0),col=rgb(0,0,0,.1),border=FALSE)\n    legend(\"topleft\",legend=c(\"Injury Probability\",\"Speed\",\"Direction Change\",\"Movement Vector Change\"),col=c(1,2,3,4),lty=c(1,2,2,2),cex=0.8,box.lty=0)\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here we can see once again how the injury probability changes over the course of the play, and we can see that the time period with the highest injury probability exhibits the expected acceleration and deceleration pairing. We can go even further, showing how all coefficients in the model combine to give the final injury probability, and we can see that the coefficients which are related to \"s32\" lead to a 2.7 times increase in the injury probability."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"coeffs <- fit2$coefficients\n#coeffs <- coeffs[c(1,7,11,5,16,3,12,20,13,15,8,9,2,18,19,23,28,24,17,10,25,26,4,27,14,21,22,6)]\n#coeffs <- coeffs[c(1,14,21,22,6)]\ncoeffs <- coeffs[c(1,7,11,5,16,3,12,20,13,15,8,9,23,28,24,10,25,26,4,27,14,21,22,6,2,18,19,17)]\n\nXY2m <- model.matrix(inj~.^2,data=XY2)\ncoeff_strings <- strsplit(names(coeffs),\":\");\ncols <- rep(0,length(coeff_strings))\nfor(k in 1:length(coeff_strings)){\n  search_str <- coeff_strings[[k]]\n  if(length(lengths(coeff_strings[[k]])) == 1){\n    cols[k] <- which(is.finite(match(colnames(XY2m),search_str)))\n  }\n  if(length(lengths(coeff_strings[[k]])) == 2){\n    search_str1 <- paste(search_str[1],search_str[2],sep=\":\")\n    search_str2 <- paste(search_str[2],search_str[1],sep=\":\")\n    colk_1 <- which(is.finite(match(colnames(XY2m),search_str1)))\n    colk_2 <- which(is.finite(match(colnames(XY2m),search_str2)))\n    if(length(colk_1)>0){cols[k] <- colk_1}\n    if(length(colk_2)>0){cols[k] <- colk_2}\n  }\n}\nfor(j in 27){\nxi <- XY2m[j,cols]; probs <- rep(0,length(xi))\nfor(k in 1:length(xi)){\n  probs[k] <- exp(xi[1:k] %*% coeffs[1:k])\n}\n\nplot(c(1,length(xi)),c(76/nrow(XY2),76/nrow(XY2)),type=\"l\",lty=2,ylim=c(0,max(probs)),axes=FALSE,xlab=\"\",ylab=\"Predicted Injury Probability\")\naxis(1,labels=rep(\"\",28), at=1:28,las=2,cex=.25)\naxis(2)\n\nend_point = .25 + 28 #this is the line which does the trick (together with barplot \"space = 1\" parameter)\ntext(seq(1.25,end_point,by=1), par(\"usr\")[3]-max(probs)/100, labels = names(coeffs), srt = 45, adj = c(1.1,1.1), xpd = TRUE, cex=0.65)\n\n#lines(probs,col=rgb(1,0,0,.2))\nfor(k in 1:length(xi)){\n  col_use <- rgb(0,0,0,.2)\n  #if(grepl(\"field\",names(coeffs)[k])|grepl(\"s22\",names(coeffs)[k])){col_use <- rgb(1,0,0,.2)}\n  if(grepl(\"s32\",names(coeffs)[k])){col_use <- rgb(1,0,0,.2)}\n  if(k < (length(xi)+1)){\n  lines(.5+c(k,k-1),c(probs[k],probs[k]),col=1)\n  }\n  #points(k,probs[k],col=col_use)\n  if(k>1){\n    #points(c(k,k),c(probs[k-1],probs[k]),col=rgb(1,0,0,.2))\n    #lines(c(k,k),c(probs[k-1],probs[k]),col=col_use)\n    polygon(.5+c((k-1),k,k,(k-1)),c(probs[k-1],probs[k-1],probs[k],probs[k]),col=col_use,border=FALSE)\n  }\n}}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Quick Acceleration\n\nA second movement type which the model suggests increases injury probability is quick acceleration. The singular vector accounting for this is \"s22\"."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plot(Ss2$u[,2],type=\"l\",ylim=c(-.5,.5),ylab=\"Speed Time History Pattern\",xlab=\"Lag Index\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And a play which exhibits the effect of this pattern on injury probability is shown below."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#play_key_all <- unique(play_data$PlayKey)\n#play_key_all <- c(as.character(play_key_all[sample(1:length(play_key_all),1000,replace=FALSE)]),as.character(play_keys))\nfor(k in 62){\nplay_num <- k\ncur_info <- play_info[which(play_info$PlayKey==as.character(inj_plays[k])),]\nsing_play_data <- play_data[which(play_data$PlayKey == as.character(inj_plays[k])),]\nN <- nrow(sing_play_data)  \nN_plays <- N - lag_max - 1  \ninj <- rep(1,N_plays); \nv1 <- rep(0,N_plays); v2 <- rep(0,N_plays); v3 <- rep(0,N_plays); v4 <- rep(0,N_plays); v5 <- rep(0,N_plays);\ns1 <- rep(0,N_plays); s2 <- rep(0,N_plays); s3 <- rep(0,N_plays); s4 <- rep(0,N_plays); s5 <- rep(0,N_plays);\nd1 <- rep(0,N_plays); d2 <- rep(0,N_plays); d3 <- rep(0,N_plays); d4 <- rep(0,N_plays); d5 <- rep(0,N_plays);\n#\nv12 <- rep(0,N_plays); v22 <- rep(0,N_plays); v32 <- rep(0,N_plays); v42 <- rep(0,N_plays); v52 <- rep(0,N_plays);\ns12 <- rep(0,N_plays); s22 <- rep(0,N_plays); s32 <- rep(0,N_plays); s42 <- rep(0,N_plays); s52 <- rep(0,N_plays);\nd12 <- rep(0,N_plays); d22 <- rep(0,N_plays); d32 <- rep(0,N_plays); d42 <- rep(0,N_plays); d52 <- rep(0,N_plays);\n#\nmaxs <- rep(0,N_plays); maxdv <- rep(0,N_plays); field <- rep(\"\",N_plays);\n# Play Types\nreturn <- rep(0,N_plays); kick <- rep(0,N_plays); rush <- rep(0,N_plays); pass <- rep(0,N_plays);\n# Weather\nsnow <- rep(0,N_plays); rain <- rep(0,N_plays); precipitation <- rep(0,N_plays);\n# Stadium\nindoor <- rep(0,N_plays); outdoor <- rep(0,N_plays);\n# Temp\ntemperature <- rep(0,N_plays)\n# Defense or Offense\ndefense <- rep(0,N_plays); offense <- rep(0,N_plays); line <- rep(0,N_plays)\n\n  #xv <- sing_play_data$x[2:(N)]-sing_play_data$x[1:(N-1)]\n  #yv <- sing_play_data$y[2:(N)]-sing_play_data$y[1:(N-1)]\n  xv <- sing_play_data$s*sin(sing_play_data$dir*pi/180)\n  yv <- sing_play_data$s*cos(sing_play_data$dir*pi/180)\n  #n <- N-1\n  n <- N\n  dv <- sqrt((xv[2:n]-xv[1:(n-1)])^2+(yv[2:n]-yv[1:(n-1)])^2)\n  ds <- sing_play_data$s[2:(N)]-sing_play_data$s[1:(N-1)]\n  ddir <- abs(sin((sing_play_data$dir[2:N]-sing_play_data$dir[1:(N-1)])*pi/180))\n  \n  s<- sing_play_data$s\n  #\n  \n  for(j in 1:N_plays){\n    rv<-acf(dv,lag.max=lag_max,type=\"covariance\",plot=FALSE)\n    rs<-acf(ds,lag.max=lag_max,type=\"covariance\",plot=FALSE)\n    rdir<-acf(ddir,lag.max=lag_max,type=\"covariance\",plot=FALSE)\n    \n    v12[j] <- dv[j:(j+lag_max)] %*% Sv2$u[,1]\n    v22[j] <- dv[j:(j+lag_max)] %*% Sv2$u[,2]\n    v32[j] <- dv[j:(j+lag_max)] %*% Sv2$u[,3]\n    v42[j] <- dv[j:(j+lag_max)] %*% Sv2$u[,4]\n    v52[j] <- dv[j:(j+lag_max)] %*% Sv2$u[,5]\n    #\n    s12[j] <- s[j:(j+lag_max)] %*% Ss2$u[,1]\n    s22[j] <- s[j:(j+lag_max)] %*% Ss2$u[,2]\n    s32[j] <- s[j:(j+lag_max)] %*% Ss2$u[,3]\n    s42[j] <- s[j:(j+lag_max)] %*% Ss2$u[,4]\n    s52[j] <- s[j:(j+lag_max)] %*% Ss2$u[,5]\n    #\n    d12[j] <- ddir[j:(j+lag_max)] %*% Sdir2$u[,1]\n    d22[j] <- ddir[j:(j+lag_max)] %*% Sdir2$u[,2]\n    d32[j] <- ddir[j:(j+lag_max)] %*% Sdir2$u[,3]\n    d42[j] <- ddir[j:(j+lag_max)] %*% Sdir2$u[,4]\n    d52[j] <- ddir[j:(j+lag_max)] %*% Sdir2$u[,5]\n  \n  #\n  v1[j] <- cor(rv$acf,Sv$u[,1])\n  v2[j] <- cor(rv$acf,Sv$u[,2])\n  v3[j] <- cor(rv$acf,Sv$u[,3])\n  v4[j] <- cor(rv$acf,Sv$u[,4])\n  v5[j] <- cor(rv$acf,Sv$u[,5])\n  #\n  s1[j] <- cor(rs$acf,Ss$u[,1])\n  s2[j] <- cor(rs$acf,Ss$u[,2])\n  s3[j] <- cor(rs$acf,Ss$u[,3])\n  s4[j] <- cor(rs$acf,Ss$u[,4])\n  s5[j] <- cor(rs$acf,Ss$u[,5])\n  #\n  d1[j] <- cor(rdir$acf,Sdir$u[,1])\n  d2[j] <- cor(rdir$acf,Sdir$u[,2])\n  d3[j] <- cor(rdir$acf,Sdir$u[,3])\n  d4[j] <- cor(rdir$acf,Sdir$u[,4])\n  d5[j] <- cor(rdir$acf,Sdir$u[,5])\n  maxs[j] <- max(sing_play_data$s[j:(j+lag_max)])\n  maxdv[j] <- max(dv[j:(j+lag_max)])\n  \n  field[j] <- as.character(cur_info$FieldType)\n  # Play Types\n  if(cur_info$PlayType == \"Return\"){return[j] <- 1}\n  if(cur_info$PlayType == \"Kick\"){kick[j] <- 1}\n  if(cur_info$PlayType == \"Rush\"){rush[j] <- 1}\n  if(cur_info$PlayType == \"Pass\"){pass[j] <- 1}\n  # Weather Types\n  if(cur_info$Weather == \"Snow\"){snow[j] <- 1; precipitation[j] <- 1}\n  if(cur_info$Weather == \"Rain\"){rain[j] <- 1; precipitation[j] <- 1}\n  # Stadium Types\n  if(cur_info$StadiumType == \"Indoor\"){indoor[j] <- 1}\n  if(cur_info$StadiumType == \"Outdoor\"){outdoor[j] <- 1}\n  # Temperature\n  temperature[j] <- cur_info$Temperature\n  # Position Types\n  cur_pos <- cur_info$PositionGroup\n  if(cur_pos==\"DB\"|cur_pos==\"DL\"|cur_pos==\"LB\"){defense[j] <- 1}\n  if(cur_pos==\"RB\"|cur_pos==\"OL\"|cur_pos==\"QB\"|cur_pos==\"TE\"|cur_pos==\"WR\"){offense[j] <- 1}\n  if(cur_pos==\"DL\"|cur_pos==\"OL\"){line[j] <- 1}\n  }\n\nplay_XY <- data.frame(inj=inj,field=as.factor(field),return=return,kick=kick,rush=rush,pass=pass,snow=snow,rain=rain,precipitation=precipitation,indoor=indoor,outdoor=outdoor,defense=defense,offense=offense,line=line)\nplay_phase_1 <- data.frame(temperature=temperature,maxs=maxs,maxdv=maxdv,v1=v1,v2=v2,v3=v3,v4=v4,v5=v5,s1=s1,s2=s2,s3=s3,s4=s4,s5=s5,d1=d1,d2=d2,d3=d3,d4=d4,d5=d5)\nplay_phase_2 <- data.frame(v12=v12,v22=v22,v32=v32,v42=v42,v52=v52,s12=s12,s22=s22,s32=s32,s42=s42,s52=s52,d12=d12,d22=d22,d32=d32,d42=d42,d52=d52)\n\nfor(k in 1:ncol(play_phase_1)){\n  play_phase_1[,k] <- (play_phase_1[,k] - standardize_mean_p1[k,1])/standardize_sd_p1[k,1]\n}\n\nfor(k in 1:ncol(play_phase_2)){\n  play_phase_2[,k] <- (play_phase_2[,k] - standardize_mean_p2[k,1])/standardize_sd_p2[k,1]\n}\n\nplay_XY2 <- cbind(play_XY,play_phase_1,play_phase_2)\n\nprob <- predict(fit2,play_XY2,type=\"response\")\n\nplot(prob/max(prob),type=\"l\",xlab=\"Sample Number\",ylab=\"Normalized Value\",ylim=c(0,1))\n    lines(sing_play_data$s/max(sing_play_data$s),col=2,lty=2)\n    lines(ddir/max(ddir),col=3,lty=2)\n    lines(dv/max(dv),col=4,lty=2)\n    ind1 <- which.max(prob)\n    polygon(c(ind1,ind1,ind1+lag_max,ind1+lag_max),c(0,1,1,0),col=rgb(0,0,0,.1),border=FALSE)\n    legend(\"topleft\",legend=c(\"Injury Probability\",\"Speed\",\"Direction Change\",\"Movement Vector Change\"),col=c(1,2,3,4),lty=c(1,2,2,2),cex=0.8,box.lty=0)\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In this case, we know that the interaction this movement type and the synthetic field is important. In fact, together, the effects of they synthetic field and this movment pattern lead to a 3.3 times increase in injury probability."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"coeffs <- fit2$coefficients\ncoeffs <- coeffs[c(1,7,11,5,16,3,12,20,13,15,8,9,2,18,19,23,28,24,17,10,25,26,4,27,14,21,22,6)]\n#coeffs <- coeffs[c(1,14,21,22,6)]\n#coeffs <- coeffs[c(1,7,11,5,16,3,12,20,13,15,8,9,23,28,24,10,25,26,4,27,14,21,22,6,2,18,19,17)]\n\nXY2m <- model.matrix(inj~.^2,data=XY2)\ncoeff_strings <- strsplit(names(coeffs),\":\");\ncols <- rep(0,length(coeff_strings))\nfor(k in 1:length(coeff_strings)){\n  search_str <- coeff_strings[[k]]\n  if(length(lengths(coeff_strings[[k]])) == 1){\n    cols[k] <- which(is.finite(match(colnames(XY2m),search_str)))\n  }\n  if(length(lengths(coeff_strings[[k]])) == 2){\n    search_str1 <- paste(search_str[1],search_str[2],sep=\":\")\n    search_str2 <- paste(search_str[2],search_str[1],sep=\":\")\n    colk_1 <- which(is.finite(match(colnames(XY2m),search_str1)))\n    colk_2 <- which(is.finite(match(colnames(XY2m),search_str2)))\n    if(length(colk_1)>0){cols[k] <- colk_1}\n    if(length(colk_2)>0){cols[k] <- colk_2}\n  }\n}\nfor(j in 62){\nxi <- XY2m[j,cols]; probs <- rep(0,length(xi))\nfor(k in 1:length(xi)){\n  probs[k] <- exp(xi[1:k] %*% coeffs[1:k])\n}\n\nplot(c(1,length(xi)),c(76/nrow(XY2),76/nrow(XY2)),type=\"l\",lty=2,ylim=c(0,max(probs)),axes=FALSE,xlab=\"\",ylab=\"Predicted Injury Probability\")\naxis(1,labels=rep(\"\",28), at=1:28,las=2,cex=.25)\naxis(2)\n\nend_point = .25 + 28 #this is the line which does the trick (together with barplot \"space = 1\" parameter)\ntext(seq(1.25,end_point,by=1), par(\"usr\")[3]-max(probs)/100, labels = names(coeffs), srt = 45, adj = c(1.1,1.1), xpd = TRUE, cex=0.65)\n\n#lines(probs,col=rgb(1,0,0,.2))\nfor(k in 1:length(xi)){\n  col_use <- rgb(0,0,0,.2)\n  if(grepl(\"field\",names(coeffs)[k])|grepl(\"s22\",names(coeffs)[k])){col_use <- rgb(1,0,0,.2)}\n  #if(grepl(\"s32\",names(coeffs)[k])){col_use <- rgb(1,0,0,.2)}\n  if(k < (length(xi)+1)){\n  lines(.5+c(k,k-1),c(probs[k],probs[k]),col=1)\n  }\n  #points(k,probs[k],col=col_use)\n  if(k>1){\n    #points(c(k,k),c(probs[k-1],probs[k]),col=rgb(1,0,0,.2))\n    #lines(c(k,k),c(probs[k-1],probs[k]),col=col_use)\n    polygon(.5+c((k-1),k,k,(k-1)),c(probs[k-1],probs[k-1],probs[k],probs[k]),col=col_use,border=FALSE)\n  }\n}}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Movement Differences between Playing Surfaces\n\nAs shown by the cumulative distributions of important movement factors and the t-tests, there are some slight but statistically significant differences between how players move on natural and synthetic fields."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plot(c(0,0),c(0,1),type=\"l\",lty=2,xlim=c(-3,3),ylim=c(0,1),main=\"\",ylab=\"Cumulative Distribution\",xlab=\"Standardized Factor\")\nlines(c(-4,4),c(0,0),lty=2,col=\"gray\")\nlines(c(-4,4),c(1,1),lty=2,col=\"gray\")\n\nNatural <- XY2$field==\"Natural\"\nSynthetic <- XY2$field==\"Synthetic\"\n\ncdf1 <- ecdf(XY2$s22[Natural])\ncdf2 <- ecdf(XY2$s22[Synthetic])\nx_seq1 <- seq(from=-4,to=4,by=.1)\nx_seq2 <- seq(from=4,to=-4,by=-.1)\npolygon(c(x_seq1,x_seq2),c(cdf1(x_seq1),cdf2(x_seq2)),col=rgb(0,1,0,.1),border=3)\n\ncdf1 <- ecdf(XY2$s32[Natural])\ncdf2 <- ecdf(XY2$s32[Synthetic])\nx_seq1 <- seq(from=-4,to=4,by=.1)\nx_seq2 <- seq(from=4,to=-4,by=-.1)\npolygon(c(x_seq1,x_seq2),c(cdf1(x_seq1),cdf2(x_seq2)),col=rgb(1,0,0,.1),border=2)\n\ncdf1 <- ecdf(XY2$maxs[Natural])\ncdf2 <- ecdf(XY2$maxs[Synthetic])\nx_seq1 <- seq(from=-4,to=4,by=.1)\nx_seq2 <- seq(from=4,to=-4,by=-.1)\npolygon(c(x_seq1,x_seq2),c(cdf1(x_seq1),cdf2(x_seq2)),col=rgb(0,0,1,.1),border=4)\n\ncdf1 <- ecdf(XY2$maxdv[Natural])\ncdf2 <- ecdf(XY2$maxdv[Synthetic])\nx_seq1 <- seq(from=-4,to=4,by=.1)\nx_seq2 <- seq(from=4,to=-4,by=-.1)\npolygon(c(x_seq1,x_seq2),c(cdf1(x_seq1),cdf2(x_seq2)),col=rgb(0,1,1,.1),border=5)\n\ncdf1 <- ecdf(XY2$v3[Natural])\ncdf2 <- ecdf(XY2$v3[Synthetic])\nx_seq1 <- seq(from=-4,to=4,by=.1)\nx_seq2 <- seq(from=4,to=-4,by=-.1)\npolygon(c(x_seq1,x_seq2),c(cdf1(x_seq1),cdf2(x_seq2)),col=rgb(1,1,0,.1),border=6)\n\nlegend(\"topleft\",legend=c(\"s22\",\"s32\",\"maxdv\",\"maxs\",\"v3\"),col=c(3,2,4,5,6),lty=1)\n\nt.test(XY2$s22[Natural],XY2$s22[Synthetic],alternative = \"two.sided\", var.equal = FALSE)\nt.test(XY2$s32[Natural],XY2$s32[Synthetic],alternative = \"two.sided\", var.equal = FALSE)\nt.test(XY2$maxdv[Natural],XY2$maxdv[Synthetic],alternative = \"two.sided\", var.equal = FALSE)\nt.test(XY2$maxs[Natural],XY2$maxs[Synthetic],alternative = \"two.sided\", var.equal = FALSE)\nt.test(XY2$v3[Natural],XY2$v3[Synthetic],alternative = \"two.sided\", var.equal = FALSE)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"However, the average differences between the movement variables in our final model do not lead to a large effect. The most significant difference is that players appear to run a little faster on synthetic fields versus natural fields. However, the average difference of speed on the two fields only leads to less than 5% of the effect of the field itself in the model. This tells us that although there are some statistically significant differences between our important movement variables across playing surfaces, they are not the reason more injuries happen on synthetic fields. It appears that there is something more fundamental to the synthetic field itself that leads to increased injury probability."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"p32 <- (mean(XY2$s32[XY2$field==\"Synthetic\"])-mean(XY2$s32[XY2$field==\"Natural\"]))*0.492/0.6076*100\np22 <- (mean(XY2$s22[XY2$field==\"Synthetic\"])-mean(XY2$s22[XY2$field==\"Natural\"]))*-.518/0.6076*100\npmaxdv <- (mean(XY2$maxdv[XY2$field==\"Synthetic\"])-mean(XY2$maxdv[XY2$field==\"Natural\"]))*-.9855/0.6076*100\npmaxs <- (mean(XY2$maxs[XY2$field==\"Synthetic\"])-mean(XY2$maxs[XY2$field==\"Natural\"]))*.5725/0.6076*100\npv3 <- (mean(XY2$v3[XY2$field==\"Synthetic\"])-mean(XY2$v3[XY2$field==\"Natural\"]))*-.4494/0.6076*100\n\nbarplot(c(100,p22,p32,pmaxdv,pmaxs,pv3),names=c(\"Synthetic\",\"s22\",\"s32\",\"maxdv\",\"maxs\",\"v3\"),col=c(2,3,3,2,2,2),ylab=\"Percentage of Effect of Synthetic Turf\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Final Conclusions\n\nIn this analysis, a framework for investigating the effect and interaction of game-specific and movement variables on injury probability in NFL plays has been presented. A generalized linear model which is able to predict 74% of non-injury plays to have below average injury probability and 75% of injury plays to have above average injury probability was created. Investigation of this model gave our three main conclusions:\n\n1. There are specific movement patterns which are related to increased injury probability, and there are indications that are also interactions between game-specific variables and some movement patterns which together lead to an even higher increase in injury probability\n2. There are also game-specific variables (field surface, weather, and play type) which also are related to increased injury probability\n3. There are slight differences in movement between synthetic and natural playing surfaces, however these differences have less of an affect than the field surface itself\n\n### How can we reduce injuries?\n\nFootball is a fast moving and hard hitting sport. It is not possible to completely remove the possibility of injury, especially when maintaining the integrity of the game. For example, according to this analysis some of the most effective methods for decreasing the non-contact injury probability are to never play in the rain and to eliminate all kickoff and punt returns. However, these changes would have large fundamental changes on how the game is played. One of the next most effective ways of decreasing non-contact injury probability would be to decrease the number of games played on synthetic playing surfaces in favor of the less dangerous natural playing surfaces. It may not be possible to completely remove the synthetic playing surfaces, however, as some field which are strictly indoors may not be able to grow natural surface playing fields as easily.\n\nA realistic path forward would be to incentivice teams to play on natural playing surfaces. For example, if a field is outdoors or can have an the roof opened in order to grow natural turf then those fields should not have synthetic surfaces. Additionally, it could be suggested that new fields being built should have retractable roofs so that natural playing surfaces can be grown, while the roof can be closed to reduce effects from rain and temperature which also increase injury probability. \n\nFurther research should also be put into why synthetic fields increase non-contact injury probability. It has been suggested that cleat-surface interactions could be in part responsible. From this analysis it appears that the interaction between synthetic fields and acceleration are significant, which could give further direction to research into cleat-surface intreactions that could be responsible, hopefully leading to newly designed synthetic surfaces which do not increase injury probability.\n\n# Future Work\n\nThe dataset provided is incredibly rich. In this work, consideration was not given to injury location or injury duration. Instead, this work focused on setting up an analysis framework which allows us to glean information about movement patterns which influence injury probability and allowing a method for estimating the specific time during the play which was most likely to cause the injury. With data from more injury plays, and re-running a similar analysis for specific injury locations or injury durations it is possible that even more important conclusions can be reached to further limit the number of non-contact injuries happening in the NFL."}],"metadata":{"kernelspec":{"display_name":"R","language":"R","name":"ir"},"language_info":{"mimetype":"text/x-r-source","name":"R","pygments_lexer":"r","version":"3.4.2","file_extension":".r","codemirror_mode":"r"}},"nbformat":4,"nbformat_minor":1}