{"cells":[{"metadata":{"_uuid":"ce71ae32f92e4fe32bd20cdb204d898d8e5690ed","_execution_state":"idle","trusted":true,"_kg_hide-output":true,"_kg_hide-input":true},"cell_type":"code","source":"###############################################################################################\n### Reducing Injury Rates on Punts by Incentivizing Fewer Punt Returns                      ###\n# Author: Konstantinos Pelechrinis, Ron Yurko, Sam Venture                                    #\n# Date: 01/07/2019                                                                            #\n###############################################################################################\n\nlibrary(tidyverse) # metapackage with lots of helpful functions\nlibrary(dplyr)\nlibrary(qdapRegex)\nlibrary(stringr)\nlibrary(LearnGeom)\nlibrary(latex2exp)\nlibrary(ggplot2)\noptions(warn = -1)\n## Running code\n\n# In a notebook, you can run a single code cell by clicking in the cell and then hitting \n# the blue arrow to the left, or by clicking in the cell and pressing Shift+Enter. In a script, \n# you can run code by highlighting the code you want to run and then clicking the blue arrow\n# at the bottom of this window.\n\nplayer_play_punts <- read.csv(\"../input/play_player_role_data.csv\")\ninjuries <- read.csv(\"../input/video_review.csv\")\npunts <- read.csv(\"../input/play_information.csv\")\n# We will create a key to identify uniquely every punt; this is the combination of the game ID and the Play ID \npunts$key = paste0(punts$GameKey,\"-\",punts$PlayID)\ninjuries$key = paste0(injuries$GameKey,\"-\",injuries$PlayID)\nplayer_play_punts$key = paste0(player_play_punts$GameKey,\"-\",player_play_punts$PlayID)\n\n### Find how many injury plays were 1. fair catches 2. out-of-bounds 3. touchbacks 4. downed punts and the same for all the punts in the dataset\n## calculate the rate of concussions for each type of punts \n\npunts_data = read.csv(\"../input/play_information.csv\")\npunts_data$key = paste0(punts_data$GameKey,\"-\",punts_data$PlayID)\n\npunts_data$PlayDescription = as.character(punts_data$PlayDescription)\nfc_punts = length(which(grepl(\"fair catch\",punts_data$PlayDescription)))\noob_punts = length(which(grepl(\"out of bounds\",punts_data$PlayDescription))) # these do not include returns where the returner was pushed out-of-bounds since those include \"pushed ob\" \ntback_punts = length(which(grepl(\"Touchback\",punts_data$PlayDescription)))\ndowned_punts = length(which(grepl(\"downed\",punts_data$PlayDescription)))\n\ninjuries_data = read.csv(\"../input/video_footage-injury.csv\")\ninjuries_data$key = paste0(injuries_data$gamekey,\"-\",injuries_data$playid)\ninjuries_data$PlayDescription= as.character(injuries_data$PlayDescription)\n\nfc_injury = length(which(grepl(\"fair catch\",injuries_data$PlayDescription)))\noob_injury = length(which(grepl(\"out of bounds\",injuries_data$PlayDescription)))\ntback_injury = length(which(grepl(\"Touchback\",injuries_data$PlayDescription)))\ndowned_injury = length(which(grepl(\"downed\",injuries_data$PlayDescription)))\n\nsprintf(\"================Concussion rates================\")\nsprintf(\"Fair Catch: %f%%\",100*(fc_injury/fc_punts))\nsprintf(\"Out-of-bounds: %f%%\",100*(oob_injury/oob_punts))\nsprintf(\"Touchbacks: %f%%\",100*(tback_injury/tback_punts))\nsprintf(\"Downed: %f%%\",100*(downed_injury/downed_punts))\nsprintf(\"Returned: %f%%\",100*((dim(injuries_data)[1]-(fc_injury+oob_injury+tback_injury+downed_injury))/(dim(punts_data)[1]-(fc_punts+oob_punts+tback_punts+downed_punts))))\n\ntype = c(\"Fair Catch\",\"Out-of-bounds\",\"Touchback\",\"Downed\",\"Returned\")\nrate = c(100*(fc_injury/fc_punts),100*(oob_injury/oob_punts),100*(tback_injury/tback_punts),100*(downed_injury/downed_punts),100*((dim(injuries_data)[1]-(fc_injury+oob_injury+tback_injury+downed_injury))/(dim(punts_data)[1]-(fc_punts+oob_punts+tback_punts+downed_punts))))\nrate.lower = 100*c(prop.test(x=fc_injury,n=fc_punts)$conf.int[1],prop.test(x=oob_injury,n=oob_punts)$conf.int[1],prop.test(x=tback_injury,n=tback_punts)$conf.int[1],prop.test(x=downed_injury,n=downed_punts)$conf.int[1],prop.test(x=(dim(injuries_data)[1]-(fc_injury+oob_injury+tback_injury+downed_injury)),n=(dim(punts_data)[1]-(fc_punts+oob_punts+tback_punts+downed_punts)))$conf.int[1])\nrate.upper = 100*c(prop.test(x=fc_injury,n=fc_punts)$conf.int[2],prop.test(x=oob_injury,n=oob_punts)$conf.int[2],prop.test(x=tback_injury,n=tback_punts)$conf.int[2],prop.test(x=downed_injury,n=downed_punts)$conf.int[2],prop.test(x=(dim(injuries_data)[1]-(fc_injury+oob_injury+tback_injury+downed_injury)),n=(dim(punts_data)[1]-(fc_punts+oob_punts+tback_punts+downed_punts)))$conf.int[2])\nsamples = c(fc_punts,oob_punts,tback_punts,downed_punts,(dim(punts_data)[1]-(fc_punts+oob_punts+tback_punts+downed_punts)))\n\nbar.data <- data.frame(Type=type,Rate=rate,lower = rate.lower, upper = rate.upper, SampleSize = samples)\n\nggplot(bar.data,aes(Type,Rate,fill=Rate))+geom_col(color=\"black\")+geom_errorbar(aes(ymin=lower, ymax=upper), colour=\"black\", width=.1)+scale_fill_gradient(low=\"grey\",high=\"red\")+labs(x=\"Punt Type\",y=\"Concussion Rate (%)\",title=\"Concussion Rates (%) for Different Types of Punts\",caption=\"Pelechrinis, Yurko, Ventura (2019)\")+theme_bw(base_size=17)+geom_text(aes(label=SampleSize),vjust=-0.35,hjust = 1.5)\n\n## Aggregate returned/non-returned punts \n\naggr.data <- data.frame(Type = c(\"Returned\",\"Non Returned\"), Rate = c(100*((dim(injuries_data)[1]-(fc_injury+oob_injury+tback_injury+downed_injury))/(dim(punts_data)[1]-(fc_punts+oob_punts+tback_punts+downed_punts))),100*(fc_injury+oob_injury+tback_injury+downed_injury)/(fc_punts+oob_punts+tback_punts+downed_punts)), lower = c(100*prop.test(x=(dim(injuries_data)[1]-(fc_injury+oob_injury+tback_injury+downed_injury)),n=(dim(punts_data)[1]-(fc_punts+oob_punts+tback_punts+downed_punts)))$conf.int[1],100*prop.test(x=fc_injury+oob_injury+tback_injury+downed_injury,n=fc_punts+oob_punts+tback_punts+downed_punts)$conf.int[1]),upper = c(100*prop.test(x=(dim(injuries_data)[1]-(fc_injury+oob_injury+tback_injury+downed_injury)),n=(dim(punts_data)[1]-(fc_punts+oob_punts+tback_punts+downed_punts)))$conf.int[2],100*prop.test(x=fc_injury+oob_injury+tback_injury+downed_injury,n=fc_punts+oob_punts+tback_punts+downed_punts)$conf.int[2]))\nggplot(aggr.data,aes(Type,Rate,fill=Rate))+geom_col(color=\"black\")+geom_errorbar(aes(ymin=lower, ymax=upper), colour=\"black\", width=.1)+scale_fill_gradient(low=\"grey\",high=\"red\")+labs(x=\"Punt Type\",y=\"Concussion Rate (%)\",title=\"Concussion Rates (%) for Returned and Not-Returned Punts\",caption=\"Pelechrinis, Yurko, Ventura (2019)\")+theme_bw(base_size=17)\n\n#### Extract the NGS for the returned punts and find y-coordinate at the time the returner receives the punt\n'%!in%' <- function(x,y)!('%in%'(x,y))\n\n## store the NGS data from all the different files to a single data frame\n# This is not the most efficient but makes the code easier to follow and support\n# the following might take some time to run\n\nNGS_files = list.files(path=\"../input\",pattern = \"NGS*\")\n\nngs.data = read.csv(paste0(\"../input/\",NGS_files[1]))\n\nfor (i in 2:length(NGS_files)){\n        ngs.tmp.f = read.csv(paste0(\"../input/\",NGS_files[i]))\n        ngs.data <- rbind(ngs.data,ngs.tmp.f)\n\n}\n\n\nngs.data$key = paste0(ngs.data$GameKey,\"-\",ngs.data$PlayID)\nind <- c(which(grepl(\"fair catch\",punts_data$PlayDescription)), which(grepl(\"out of bounds\",punts_data$PlayDescription)), which(grepl(\"Touchback\",punts_data$PlayDescription)),which(grepl(\"downed\",punts_data$PlayDescription)))\n\nloc.punt.rec <- data.frame(key=c(),y=c())\n\nfor (i in 1:dim(punts_data)[1]){\n\n        if (i %!in% ind){\n\n                ngs.tmp <- ngs.data[which(ngs.data$key == punts_data[i,]$key),]\n                if ( dim(ngs.tmp)[1] > 0){\n                        # find the returer\n                        punts.roles.tmp <- player_play_punts[which(player_play_punts$key == punts_data[i,]$key),]\n                        pr = punts.roles.tmp[which(punts.roles.tmp$Role == \"PR\"),]$GSISID\n                        # find his y-coord during the reception of the punt \"punt_received\"\n                        y = ngs.tmp[which(as.integer(as.numeric(as.character(ngs.tmp$GSISID))) == pr & as.character(ngs.tmp$Event) == \"punt_received\"),]$y\n                        if (length(as.numeric(as.character(y))) > 0){ # making sure the coordinates are not \"NAs\" -- just double checking to avoid cases where the punts have not been correctly annotated\n                                loc.punt.rec <- rbind(loc.punt.rec,data.frame(key=punts_data[i,]$key,y=as.numeric(as.character(y))))\n                        }\n                }\n\n        }\n\n}\n\n# find the closest sideline \n\nd.side <- data.frame(key=c(),d = c())\n\nfor (i in 1:dim(loc.punt.rec)[1]){\n\n        d = min(abs(53.3-loc.punt.rec[i,]$y),abs(loc.punt.rec[i,]$y-0)) # use abs because in 2 cases the y-coord is greater than 53.3 most probably due to either measurement error\n        d.side <- rbind(d.side,data.frame(key=loc.punt.rec[i,]$key,d = d))\n\n}\n\nd.side <- d.side %>% mutate(inj = ifelse(key %in% injuries$key,1,0))\n\ndsidinj.mod <- glm(inj~d,data=d.side,family=\"binomial\")\n\ndistance <- c(0:25)\n\nglm.plot <- data.frame(d=distance, prob = 100*predict(dsidinj.mod,newdata=data.frame(d=distance),type=\"response\"),se = 100*predict(dsidinj.mod,newdata=data.frame(d=distance),type=\"response\",se=T)$se.fit)\n\nannotations <- data.frame(loc = c(13,22),type=c(\"Numbers\",\"Hashmarks\"))\nggplot(glm.plot,aes(d,prob))+geom_line(color=\"blue\")+geom_errorbar(aes(ymin=prob-se, ymax=prob+se), colour=\"black\", width=.1)+labs(x=\"Distance from the closest sideline (yards)\",y=\"Concussion Incident Probability of a Returned Punt(%)\",title=\"Concussion Probability Based on the Distance from the Nearest Sideline\",caption=\"Pelechrinis, Yurko, Ventura (2019)\")+theme_bw(base_size=17)+geom_vline(data=annotations, mapping=aes(xintercept=loc), color=\"red\",linetype=\"dashed\")+geom_text(data=annotations, mapping=aes(x=loc, y=0, label=type), size=6, angle=90, vjust=-0.4, hjust=0)\n\n########################### Punt Angles ######################### \n\n### Find the angle from the 0 degree line for the punts (returned, fair catches -- for downed punts we cannot really tell the angle since we do not know the place where the ball landed first, while for out of bounds punts we do not know which sideline they were out from) behind a team's own 30 yard line\n## Steps: find the punts of interest, find the location of the punter at the time of punt and the location of the punt returner at the time of fair-catch/punt reception\n\nind.punt.RetFC <- c(which(grepl(\"out of bounds\",punts_data$PlayDescription)), which(grepl(\"Touchback\",punts_data$PlayDescription)),which(grepl(\"downed\",punts_data$PlayDescription)))\n\nloc.punt.RetFC <- data.frame(key=c(),yrdline=c(),x1=c(),y1=c(),x2=c(),y2=c(),theta=c())\n\n## add an indicator in the punts_data on whether the punt is on the own territory behind the 30 yard line \n\npunts_data = cbind(punts_data, read.table(text = as.character(punts_data$YardLine), sep = \" \"))\npunts_data <- punts_data %>% mutate(own30 = ifelse(as.numeric(punts_data$V2)<=30,1,0))\n\nfor (i in 1:dim(punts_data)[1]){\n\n        if (i %!in% ind.punt.RetFC & punts_data[i,]$own30 == 1){\n\n                ngs.tmp <- ngs.data[which(ngs.data$key == punts_data[i,]$key),]\n                if ( dim(ngs.tmp)[1] > 0){\n                        # find the returner\n                        punts.roles.tmp <- player_play_punts[which(player_play_punts$key == punts_data[i,]$key),]\n                        pr = punts.roles.tmp[which(punts.roles.tmp$Role == \"PR\"),]$GSISID\n                        # find his xy-coord during the reception of the punt \"punt_received\" or the fair catch \"fair_catch\"\n                        x2 = ngs.tmp[which(as.integer(as.numeric(as.character(ngs.tmp$GSISID))) == pr & (as.character(ngs.tmp$Event) == \"punt_received\" |as.character(ngs.tmp$Event) == \"fair_catch\" )),]$x\n                        y2 = ngs.tmp[which(as.integer(as.numeric(as.character(ngs.tmp$GSISID))) == pr & (as.character(ngs.tmp$Event) == \"punt_received\"|as.character(ngs.tmp$Event) == \"fair_catch\" )),]$y\n                        # find the punter \n                        punter = punts.roles.tmp[which(punts.roles.tmp$Role == \"P\"),]$GSISID\n                        # find his xy-coord during the punt (Event: \"punt\")\n                        x1 = ngs.tmp[which(as.integer(as.numeric(as.character(ngs.tmp$GSISID))) == punter & as.character(ngs.tmp$Event) == \"punt\"),]$x\n                        y1 = ngs.tmp[which(as.integer(as.numeric(as.character(ngs.tmp$GSISID))) == punter & as.character(ngs.tmp$Event) == \"punt\"),]$y\n                        if (length(as.numeric(as.character(y1))) > 0 & length(as.numeric(as.character(y2))) ){ # making sure the coordinates are not \"NAs\" -- just double checking to avoid cases where the punts have not been correctly annotated \n                                # find the angle between the straight line defined by the 0 degree line from the punter [(x1,y1),(x2,y1)] and the actual line of the punt [(x1,y1),(x2,y2)]\n                                A = c(as.numeric(as.character(x2)),as.numeric(as.character(y2)))\n                                B = c(as.numeric(as.character(x1)),as.numeric(as.character(y1)))\n                                C = c(as.numeric(as.character(x2)),as.numeric(as.character(y1)))\n                                theta = Angle(A, B,C)[[1]]\n                                loc.punt.RetFC <- rbind(loc.punt.RetFC,data.frame(key=punts_data[i,]$key,yrdline=as.numeric(punts_data[i,]$V2),x1=as.numeric(as.character(x1)),y1=as.numeric(as.character(y1)),x2=as.numeric(as.character(x2)),y2=as.numeric(as.character(y2)),theta=theta))\n                        }\n\n                }\n\n        }\n\n}\n\n\n####### Possible SIDEFFECTS - Muffed Punts\n\n### One of the possible \"side-effects\" is that a punt returner might try to catch on the fly more punts for a fair catch that other wise he would have let land and downed by the covering team. This might lead to more muffed punts so we need to examine what is the rate of concussions in muffed punts.\n\nmuffed_punts = length(which(grepl(\"MUFFS\",punts_data$PlayDescription)))\nmuffed_injury = length(which(grepl(\"MUFFS\",injuries_data$PlayDescription)))\nmuffed_rate.lower = c(prop.test(x=muffed_injury,n=muffed_punts)$conf.int[1])\nmuffed_rate.upper = c(prop.test(x=muffed_injury,n=muffed_punts)$conf.int[2])\n\nsprintf(\"================Concussion rates================\")\nsprintf(\"Muffed Punts: %f%%\",100*(muffed_injury/muffed_punts))\nsprintf(\"Muffed Punts 95%% CI: [%f%%, %f%%]\",100*muffed_rate.lower,100*muffed_rate.upper)\n\n# but some muffed punts were signaled for fair catch and some were attempted to be returned. The NGS have annotations for the fair catch signal so we will find which muffed punts were signaled for fair catch and which were not (i.e., they would be attempted to be returned).\n\nmuffed_keys = punts_data[which(grepl(\"MUFFS\",punts_data$PlayDescription)),]$key\n\nmuffed_punts.data = data.frame(key=c(),fc=c())\n\nfor (i in 1:length(muffed_keys)){\n        ngs.tmp <- ngs.data[which(ngs.data$key == muffed_keys[i]),]\n        e = which(as.character(ngs.tmp$Event) == \"fair_catch\")\n        if (length(e) > 0){\n                muffed_punts.data = rbind(muffed_punts.data,data.frame(key=muffed_keys[i],fc = 1))\n        }else{\n                muffed_punts.data = rbind(muffed_punts.data,data.frame(key=muffed_keys[i],fc = 0))\n        }\n}\n\nmuffed_inj_keys = injuries_data[which(grepl(\"MUFFS\",injuries_data$PlayDescription)),]$key\n\nmuffed_punts.data <- muffed_punts.data %>% mutate(inj = ifelse(as.character(key) %in% muffed_inj_keys,1,0))\n\nmuffed_punts_fc = length(which(muffed_punts.data$fc == 1))\nmuffed_punts_fc_inj = length(which(muffed_punts.data$fc == 1 & muffed_punts.data$inj == 1))\nmuffed_punts_nfc = length(which(muffed_punts.data$fc == 0))\nmuffed_punts_nfc_inj = length(which(muffed_punts.data$fc == 0 & muffed_punts.data$inj == 1))\n\nmuffed_punts_fc.lower = c(prop.test(x=muffed_punts_fc_inj,n=muffed_punts_fc)$conf.int[1])\nmuffed_punts_fc.upper = c(prop.test(x=muffed_punts_fc_inj,n=muffed_punts_fc)$conf.int[2])\nmuffed_punts_nfc.lower = c(prop.test(x=muffed_punts_nfc_inj,n=muffed_punts_nfc)$conf.int[1])\nmuffed_punts_nfc.upper = c(prop.test(x=muffed_punts_nfc_inj,n=muffed_punts_nfc)$conf.int[2])\n\nmuffed.rates.dat <- data.frame(type=c(\"All Muffed\",\"Muffed FC\",\"Muffed Non-FC\"),rate=c(muffed_injury/muffed_punts,muffed_punts_fc_inj/muffed_punts_fc,muffed_punts_nfc_inj/muffed_punts_nfc),lower = c(muffed_rate.lower,muffed_punts_fc.lower,muffed_punts_nfc.lower),upper=c(muffed_rate.upper,muffed_punts_fc.upper,muffed_punts_nfc.upper))\n\nmuffed.rates.dat$Rate = 100*muffed.rates.dat$rate\n\nggplot(muffed.rates.dat,aes(type,Rate,fill=Rate))+geom_col(color=\"black\")+geom_errorbar(aes(ymin=100*lower, ymax=100*upper), colour=\"black\", width=.1)+scale_fill_gradient(low=\"grey\",high=\"red\")+labs(x=\"Muffed Punt Type\",y=\"Concussion Rate (%)\",title=\"Concussion Rates (%) for Different Types of Muffed Punts\",caption=\"Pelechrinis, Yurko, Ventura (2019)\")+theme_bw(base_size=17)+ylim(c(0,5))\n\n\n####### SIMULATE the impact of the proposed rule change #######\n### we want to get an estimate of the number of injuries that will be prevented by this rule\n## we will perform simulations based on assumptions for the percentage of punts to not be returned and for the closeness to the sideline of the kicks that will be returned \n## this will provide us with a range of possibilities and then we can obtain a fairly realistic view of what to expect\n## Currently ~47% of the punts are returned. Based on the simulated reduction in punts returned this will be adjusted accordingly\n\nreduction_rate = c(0.05,0.1,0.15,0.2) # this is an assumed reduction rate for the returned punts. Realistically we should not expect more than 20% reduction\nyards_closer = c(2.5,5,7.5,10) # this is the assumed expected shift of the punt towards the sideline. Assume a normal distribution with variance 7 yards (this is the variance of the sideline distance in the punts observed) \n\ninjuries_simulations = data.frame(red_rate = c(),yrds_closer = c(), xCR = c(), lower = c())\n \nfor (r in reduction_rate){\n \n         for (y in yards_closer){\n                 # we estimate the rate per 1,000 exposures, so we assume 1,000 punts\n                 # based on the base rates 470 of them would be expected to be returned and 530 to not be returned \n                 xPR = 470 - round(r*470) \n                 xNPR = 530 + round(r*470)\n                 xInjNPR = aggr.data[2,]$Rate*0.01*xNPR # aggr.data data frame has the expected concussion rates - as well as upper and lower bounds - for returned and non-returned punts as calculated above\n                 xInjNPR_lower = aggr.data[2,]$lower*0.01*xNPR\n                 boot_samples = rep(0,500)\n                 boot_samples_lower = rep(0,500)\n                 for (b in 1:500){\n                         returned_dis_sim = sample(d.side$d,xPR,replace=T)-rnorm(xPR,y,7)\n                         returned_dis_sim[which(returned_dis_sim<0)] = 1\n                         returned_dis_sim[which(returned_dis_sim>26)] = 26\n                         xInjRet = sum(predict(dsidinj.mod,newdata=data.frame(d=returned_dis_sim),type=\"response\"))\n                         xInjRet_lower = xInjRet-sum(predict(dsidinj.mod,newdata=data.frame(d=returned_dis_sim),type=\"response\",se.fit=T)$se.fit)\n                         boot_samples[b] = xInjRet\n                         boot_samples_lower[b] = xInjRet_lower\n                 }\n                 injuries_simulations <- rbind(injuries_simulations,data.frame(red_rate = r, yrds_closer= y,xCR = xInjNPR+mean(boot_samples),lower=xInjNPR_lower+mean(boot_samples_lower)))\n \n         }\n \n}\n\ncurrent_concussion_rate_per1Kexposures = (dim(injuries)[1]/dim(punts_data)[1])*1000\n\nggplot(injuries_simulations, aes(x = red_rate, y = xCR,color=as.factor(yrds_closer)))+geom_line(size=2) + geom_hline(yintercept=current_concussion_rate_per1Kexposures,linetype=\"dashed\",size=2)+scale_color_manual(TeX(\"$s_d$\"),values=c(\"blue\",\"brown\",\"red\",\"orange\",\"purple\"))+labs(x=TeX(\"r\"),y=\"# of Concussions per 1,000 Exposures (Punts)\",title=\"Our Suggested Rule Changes Expected Impact on Concussion Rates\",caption=\"Pelechrinis, Yurko, Ventura (2019)\")+theme_bw(base_size=17)+geom_text(mapping=aes(x=0.075, y=5.6, label=\"Current # of Concussions per 1,000 Exposures\"), size=6,hjust=0,vjust=-0.05,color=\"black\")\n\nggplot(injuries_simulations[which(injuries_simulations$red_rate == 0.05),], aes(x = yrds_closer))+geom_line(aes(y=xCR),size=2,color=\"red\")+geom_line(aes(y=lower),size=2,color=\"brown\",linetype=\"dashed\") + geom_hline(yintercept=current_concussion_rate_per1Kexposures,linetype=\"dashed\",size=2)+labs(x=TeX(\"$s_d$\"),y=\"# of Concussions per 1,000 Exposures (Punts)\",title=\"Our Suggested Rule Changes Expected Impact on Concussion Rates\",caption=\"Pelechrinis, Yurko, Ventura (2019)\")+theme_bw(base_size=17)+geom_text(mapping=aes(x=2.075, y=5.6, label=\"Current # of Concussions per 1,000 Exposures\"), size=6,hjust=0,vjust=-0.05,color=\"black\")+xlim(2,11)+geom_text(mapping=aes(x=3.5, y=4, label=\"This is the region of possible conussion rates to be observed with our suggested rule change\"), size=6,hjust=0,vjust=-0.05,color=\"black\")+geom_text(mapping=aes(x=3.5,y = 5.1,label=\"This is the expected number of concussions per 1,000 exposures with our suggested rule change\"), size=6,hjust=0,vjust=-0.05,color=\"red\")+geom_text(mapping=aes(x=3.5,y = 2.4,label=\"This is a conservative lower bound number of concussions per 1,000 exposures with our suggested rule change\"), size=6,hjust=0,vjust=-0.05,color=\"brown\")\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d90cda2ad532a7c328c40bfb887c847d2c85f270"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}