library(ggplot2)
library(forcats)

###############
# Import Data #
###############

nfl_final <- read.csv('../input/cleaning-game-player-data-py/nfl_final.csv', header = T)

#####################
# Null Hypothesises #
#####################

# Null Hypothesis 1: Players with ball possession and players with no ball position have the same association for increased concussions during punt plays.
# Null Hypothesis 2: NFL games on turf fields and grass fields have the same association for increased concussions.
# Null Hypothesis 3: 'Gunners' have the same association for increased concussion as other punt posiitons.
# Null Hypothesis 4: Players on the receiving side and kicking side of the play have the same association for increased concussions.
# Null Hypothesis 4.1: Punt returners have the same association for increased concussion as other punt posiitons.

###########################
# Sample Size Calculation #
###########################

# logic regression requires a large sample size; below is the calculation to determine sample size based on number of dependent variables
# sample size = ((10 * number of independent variables)/probability of least frequent outcome)
sample_size.1 <- as.integer(((10*1)/(37/146573)))
sample_size.2 <- as.integer(((10*2)/(37/146573)))
sample_size.3 <- as.integer(((10*3)/(37/146573)))
sample_size.4 <- as.integer(((10*4)/(37/146573)))

# Due to our sample size, we can have 3 dependent variables in a model
# Sample size needed for 3 dependent variables: n=118842

####################################
# Data Analysis - Whole Population #
####################################

## Null Hypothesis 1: Players with ball possession and players with no ball position have the same association for increased concussions during punt plays.
# Defining players with ball possession as Punter and Kick Returner
# Y value: 'concussion_status'; whether or not an individual has a concussion
# X value: 'ball_posession'; was the player a Punter (P) or a Kick Returner (PR)
# Stat Test: logisitc regression

nfl_ball <- nfl_final

# Creating variable 'ball_posession'
nfl_ball$ball_posession <- ifelse(nfl_ball$Role == "P" | nfl_ball$Role == "PR", 1, 0)

# logistic regression
ball_glm <- glm(concussion_status ~ ball_posession, data = nfl_ball, family = binomial)
summary(ball_glm)

# p > 0.05, cannot reject null hypothesis

# does accounting for confounders change things?

# Multiple logistic regression to adddress confounders
# Controlling for GameWeather and Quarter
ball_glm <- glm(concussion_status ~ 
                  ball_posession +
                  GameWeather +
                  Quarter,
                data = nfl_ball,
                family = binomial
                )
summary(ball_glm)

# p still > 0.05, cannot reject null hypothesis

## Null Hypothesis 2: NFL games on turf fields and grass fields have the same association for increased concussions.
# Y value: 'concussion_status'; whether or not an individual has a concussion
# X value: 'Turf'; was the playing field grass or turf
# Stat Test: logisitc regression

nfl_turf <- nfl_final

#Turning all variable labels into a coded number
nfl_turf$Turf <- ifelse(nfl_final$Turf == "Turf", 1, 0)

# logistic regression
turf_glm <- glm(concussion_status ~ Turf, data = nfl_turf, family = binomial)
summary(turf_glm)

# p > 0.05, cannot reject null hypothesis

# does accounting for confounders change things?

# Multiple logistic regression to adddress confounders
# Controlling for GameWeather
turf_glm <- glm(concussion_status ~ 
                  Turf +
                  GameWeather +
                  Turf*GameWeather,
                data = nfl_turf,
                family = binomial
                )
summary(turf_glm)

# p still > 0.05 and worse than previous test

## Null Hypothesis 3: 'Gunners' have the same association for increased concussion as other punt posiitons.
# Y value: 'concussion_status'; whether or not an individual has a concussion
# X value: 'gunner'; is the NFL player a 'punt returner'gunner' punt position
# Stat Test: logisitc regression

nfl_gunner <- nfl_final

# coding whether or not player is a gunner (GL or GR)
nfl_gunner$gunner <- ifelse(nfl_gunner$Role == "GL" | nfl_gunner$Role == "GR", 1, 0)

# logistic regression
gunner_glm <- glm(concussion_status ~ gunner, data = nfl_gunner, family = binomial)
summary(gunner_glm)

# p > 0.05, cannot reject null hypothesis

# Does controlling for observed confounders change anything?

# Multiple logistic regression to adddress confounders
# Controlling for Quarter and Stadium Type (i.e. indoor or outdoor)
gunner_glm <- glm(concussion_status ~
                gunner +
                Quarter +
                StadiumType,
              data = nfl_gunner,
              family = binomial
              )
summary(gunner_glm)

# p > 0.05, cannot reject null hypothesis

## Null Hypothesis 4: Players on the receiving side and kicking side of the play have the same association for increased concussions.
# Y value: 'concussion_status'; whether or not an individual has a concussion
# X value: 'side'; is the NFL player on recieving side or kicking side
# Stat Test: logisitc regression

nfl_side <- nfl_final

# coding whether or not player is on kicking or receiving side
# 0 = kicking, 1 = receiving
nfl_side$side <- ifelse(nfl_side$Role == "GL" |
                          nfl_side$Role == "PLW" |
                          nfl_side$Role == "PLT" |
                          nfl_side$Role == "PLG" |
                          nfl_side$Role == "PLS" |
                          nfl_side$Role == "PRG" |
                          nfl_side$Role == "PRT" |
                          nfl_side$Role == "PRW" |
                          nfl_side$Role == "PC" |
                          nfl_side$Role == "PPR" |
                          nfl_side$Role == "P" |
                          nfl_side$Role == "GR",
                        0, 1
                        )

# logistic regression
side_glm <- glm(concussion_status ~ side, data = nfl_side, family = binomial)
summary(side_glm)
confint(side_glm)

# p value = 0.0039, results are statistically significant

# converting log odds (i.e. logistic regression coefficient) to probability
(exp(1.0694)/(1 + exp(1.0694)))

# probability = 0.74
# In other words, being a being a NFL punt recieving side increases your liklihood of a concussion by 50%

# Does controlling for observed confounders change anything?

# Multiple logistic regression to adddress confounders
# Controlling for Quarter and Stadium Type (i.e. indoor or outdoor)
side_glm <- glm(concussion_status ~
                side +
                Quarter +
                StadiumType,
              data = nfl_side,
              family = binomial
              )

summary(side_glm)
options(warn=-1) # turns off warning message
confint(side_glm)
options(warn=-0) # turns on warning message

# p value still 0.0039, results are statistically significant

# With this, let's go further into understanding this increase for the punt kick recieivng side

## Null Hypothesis 4.1: Punt returners have the same association for increased concussion as other punt posiitons.
# Y value: 'concussion_status'; whether or not an individual has a concussion
# X value: 'PR'; is the NFL player a punt returner
# Stat Test: logisitc regression

nfl_pr <- nfl_final

# coding whether or not player is a punt returner (PR)
nfl_pr$pr <- ifelse(nfl_pr$Role == "PR", 1, 0)

# logistic regression
pr_glm <- glm(concussion_status ~ pr, data = nfl_pr, family = binomial)
summary(pr_glm)
confint(pr_glm)

# p value = 0.0139, results are statistically significant

# converting log odds (i.e. logistic regression coefficient) to probability
(exp(-1.1832)/(1 + exp(-1.1832)))

# probability = 0.23
# In other words, being a being a NFL punt returner increases your liklihood of a concussion by 50%

# Does controlling for observed confounders change anything?

# Multiple logistic regression to adddress confounders
# Controlling for Quarter and Stadium Type (i.e. indoor or outdoor)
pr_glm <- glm(concussion_status ~
                pr +
                Quarter +
                StadiumType,
              data = nfl_pr,
              family = binomial
              )
summary(pr_glm)
options(warn=-1) # turns off warning message
confint(pr_glm)
options(warn=-0) # turns on warning message

# p value still 0.0139, results are statistically significant

##################################################
# Bar Charts - Injury Descriptive Statistics  #
##################################################

# Will look at descreptive statistics of the NFL concussed punt player population
# Though we can't show causality from just looking at the descreptive statistics, it adds pieces of info to the larger picture

########################
# Subsetting Dataframe #
########################

# subsetting to dataframe to only players with concussion
# still keeping column 'concussion_status'
nfl_concussion <- subset(nfl_final, subset = concussion_status == 'concussion', select = c(Season_Year:concussion_status))

####################
# Visualizing Data #
####################

## Bar Chart: Y = Number of Concussions, X = Covariate

## Season_Type

season_t.plt <- 
  ggplot(nfl_concussion, aes(Season_Type)) +
  geom_bar() + 
  labs(x = "Season Type",
       y = 'Number of Concussions',
       title = 'Total Number of Punt Play Concussions by Season Type',
       subtitle = '2016-2017 Seasons Combined'
       ) +
  scale_y_continuous(breaks = seq(1, 30, by = 2))

season_t.plt

## StadiumType

stadium_t.plt <- 
  ggplot(nfl_concussion, aes(StadiumType)) +
  geom_bar() + 
  labs(x = "Stadium Type",
       y = 'Number of Concussions',
       title = 'Total Number of Punt Play Concussions by Stadium Type',
       subtitle = '2016-2017 Seasons Combined'
       ) +
  scale_y_continuous(breaks = seq(1, 30, by = 2))

stadium_t.plt

## Turf

turf.plt <- 
  ggplot(nfl_concussion, aes(Turf)) +
  geom_bar() + 
  labs(x = "Turf Type",
       y = 'Number of Concussions',
       title = 'Total Number of Punt Play Concussions by Turf Type',
       subtitle = '2016-2017 Seasons Combined'
       ) +
  scale_y_continuous(breaks = seq(1, 30, by = 2))

turf.plt

## Player_Activity_Derived

nfl_concussion$Player_Activity_Derived <- fct_infreq(nfl_concussion$Player_Activity_Derived)
player_activity.plt <- 
  ggplot(nfl_concussion, aes(Player_Activity_Derived)) +
  geom_bar() + 
  labs(x = "Player Activity",
       y = 'Number of Concussions',
       title = 'Total Number of Punt Play Concussions by Player Activity',
       subtitle = '2016-2017 Seasons Combined'
  ) +
  scale_y_continuous(breaks = seq(1, 30, by = 2))

player_activity.plt

## Primary_Partner_Activity_Derived

nfl_concussion$Primary_Partner_Activity_Derived <- fct_infreq(nfl_concussion$Primary_Partner_Activity_Derived)
partner_activity.plt <- 
  ggplot(nfl_concussion, aes(Primary_Partner_Activity_Derived)) +
  geom_bar() + 
  labs(x = "Partner Activity",
       y = 'Number of Concussions',
       title = 'Total Number of Punt Play Concussions by Partner Activity',
       subtitle = '2016-2017 Seasons Combined'
  ) +
  scale_y_continuous(breaks = seq(1, 30, by = 2))

partner_activity.plt

## Punt_Position

nfl_concussion$Role <- fct_infreq(nfl_concussion$Role)
punt_position.plt <- 
  ggplot(nfl_concussion, aes(Role)) +
  geom_bar() + 
  labs(x = "Punt Position",
       y = 'Number of Concussions',
       title = 'Total Number of Punt Play Concussions by Punt Position',
       subtitle = '2016-2017 Seasons Combined'
  ) +
  scale_y_continuous(breaks = seq(1, 10, by = 1))

punt_position.plt

## Primary_Impact_Type

nfl_concussion$Primary_Impact_Type <- fct_infreq(nfl_concussion$Primary_Impact_Type)
impact_t.plt <- 
  ggplot(nfl_concussion, aes(Primary_Impact_Type)) +
  geom_bar() + 
  labs(x = "Impact Type",
       y = 'Number of Concussions',
       title = 'Total Number of Punt Play Concussions by Impact Type',
       subtitle = '2016-2017 Seasons Combined'
  ) +
  scale_y_continuous(breaks = seq(1, 30, by = 2))

impact_t.plt

## Friendly_Fire

nfl_concussion$Friendly_Fire <- fct_infreq(nfl_concussion$Friendly_Fire)
friend_fire.plt <- 
  ggplot(nfl_concussion, aes(Friendly_Fire)) +
  geom_bar() + 
  labs(x = "Friendly Fire",
       y = 'Number of Concussions',
       title = 'Total Number of Punt Play Concussions by Friendly Fire',
       subtitle = '2016-2017 Seasons Combined'
  ) +
  scale_y_continuous(breaks = seq(1, 30, by = 2))

friend_fire.plt

## Week

week.plt <- 
  ggplot(nfl_concussion, aes(Week)) +
  geom_bar() + 
  labs(x = "Game Week",
       y = 'Number of Concussions',
       title = 'Total Number of Punt Play Concussions by Week',
       subtitle = '2016-2017 Seasons Combined'
       ) +
  scale_x_continuous(breaks = seq(1, 17, by = 1)) +
  scale_y_continuous(breaks = seq(1, 8, by = 1))

week.plt

###################
# End of Analysis #
###################