Vignette.qmd

---
title: "Vignette - NFL Game Outcome Prediction Modeling"
author: "Anshi Arora, Joshua Charfauros, Christina Cui, Sean Reagan"
format: html
editor: visual
---

```{r setup, echo=FALSE, warning=FALSE, message=FALSE}
library(nflfastR)
library(knitr)
library(tidyverse)
library(dplyr)
library(ggplot2)
library(ggrepel)
library(nflreadr)
library(ggpath)
library(nflplotR)
library(caret)
library(keras)
library(fastDummies)
library(reticulate)
library(randomForest)
library(visdat)
```

# Introduction

In the NFL world, predicting game outcomes is a highly sought after accomplishment for a wide range of stakeholders like teams, analysts, fans, etc. A wide variety of approaches and methods have been utilized in sports forecasting but the application of random forests to this field is not a widespread practice. By leveraging historical play-by-play data, which includes detailed information on team performance, player actions, game context, and situational factors, we aim to build a model that can forecast the winner of an NFL game with a high degree of accuracy.

# Overview of the NFLFastR Library and Dataset

The NFLFastR Library is a package in R that provides access to detailed play-by-play data for NFL games. Statistics about every play, including yardages, player stats, game situations, win probability, down, etc., are included in the datasets. Overall information about each game (final score, precipitation, home team, weekday, etc.) is also included.

The package includes helper functions to pull specific data. We will specifically be using the following function:

load_pbp(): downloads play-by-play data for a given season (or multiple)

The usage of this function to pull play-by-play data for years 2022-2024 is displayed below. A small subset of the dataframe is also shown.

```{r, cache = T, message = F, warning = F}
pbp_example <- load_pbp(2022:2024)
kable(head(pbp_example[,c(1, 2, 4, 11, 15, 26, 13)], 5))
```

As you can see the first few observations are different plays from the same game.

The dataset for each year's play-by-play data has 372 variables in total. The definitions for each variable can be found in the library's directory. Here is the link to look through the variables:

[www.nflfastr.com/articles/stats_variables.html](https://www.nflfastr.com/articles/stats_variables.html)

# Objectives and Methodology

The objective of this vignette is to use a multitude of variables to predict binary win/loss outcomes of a game.

To determine which variables have strong correlations with game win, and thereby likely will serve as strong predictors, we will be conducting some exploratory data analysis. Then, we will train a random forest model on the data. After making the model, we can evaluate its accuracy on the test set and account for any issues that arise. We will also calculate variable importance scores to determine which predictors serve the largest roles in determining the prediction. This model is further developed by adding training controls.

More information on why/how random forests were used can be found in their respective sections below.

Here is a visualization of this vignette's methodological steps:

### Steps

```{mermaid}
flowchart TD
    A["Review NFLFastR Library: 
    Get Familiar with the Functions, Structure of Data, Variables"]
    B["Exploratory Data Analysis:
    Looking at Trends over Time, Correlations between Variables etc."]
    C["Data Preprocessing:
    Cleaning, Formatting, Choosing Variables, Creating Train/Test Sets"]
    E["Building & Training the Random Forest Model"]
    F["Testing model & Calculating Accuracy"]
    G["Evaluate Variable Importance Scores and Adjust RF Model"]
    H["Retrain Model and Calculate Test Accuracy"]
    A --> 
    B --> 
    C -->
    E --> 
    F -->
    G -->
    H
```

# Exploratory Data Analysis

Below is a graph displaying the teams that have the highest Offense EPA (Expected Points Added) and Defense EPA. Here we are looking for the teams in the top right quadrant of the graph which means that they have a negative Defense EPA (take points from their opponent) and a positive Offense EPA (add points to their own team). The top teams by this metric are the Bills, 49ers, Chiefs, and Eagles which were all really good teams during this time period.

```{r, cache=TRUE}
EPA <- readRDS('RDS files/EPA_2024.rds')
EPA
```

Another interesting component to look at is the type of stadium that the teams are playing in. Home teams are favored by both Vegas and their actual results when playing in a dome or outdoors. This advantage flips the other way for closed or open stadiums. For some reason, stadiums with retractable roofs make it less likely for their home team to take home a win.

```{r, cache = T}
LV_vs_ACT = readRDS('RDS files/LV_vs_Act.rds')
LV_vs_ACT
```

This following chart is attempting to see the strength of point differential on wins. With the chart, we can see that most teams that have a positive point differential go on to win 8 or more games. This would make sense as those teams are at least even in game or have a positive record. Some outliers from the positive side of win percentage drop down below 0 in the point differential category. This means that the relationship between number of wins and point differential is not 1-to-1.

```{r, warning=FALSE, cache=TRUE}
ptdiff = readRDS('RDS files/pt_diff.rds')
ptdiff
```

Here is a graph about quarterback efficiency. Notice anything? Those same teams as before populate our top spots. Teams like the Chiefs, 49ers, Bills, and Eagles have the most efficient or near the most efficient quarterbacks in the entire league. This proves to be another good indicator of a teams overall success. Which is good quarterback play.

```{r warning=FALSE, cache=TRUE}
QB_eff = readRDS('RDS files/QB_efficiency.rds')
QB_eff
```

This final graph addresses the point totals at the end of games and whether or not the home team is favored based off of that. From the stacked bar chart we can see that the home team is always favored. However, for relatively low or high scoring point totals, that advantage drops slightly. Games that fall within the 38-49 point total typically favor the home team a little bit more than usual.

```{r, cache=TRUE}
Total_wp = readRDS('RDS files/Stacked_Total.rds')
Total_wp
```

# Preprocessing

NFLFastR already has a pre built, schedule data frame which has each observation as a an NFL game. It was built by Lee Sharpe and is explained in this article [https://www.nflfastr.com/articles/beginners_guide.html#real-life-example-lets-make-a-win-total-model](https://www.nflfastr.com/articles/beginners_guide.html#real-life-example-lets-make-a-win-total-model). To access this, we use the load_schedules function from NFLFastR.

```{r}
games <- nflreadr::load_schedules(2021:2024)
```

Sharpe's set, however, captures a lot of data we don't need, things like quarterback and coach name etc, and doesn't capture plenty of things we do care about, such as quantitative performance metrics. To solve this, we need to manipulate our play-by-play data to have single games as observations, scrape our new game-by-game data for our useful metrics, and then merge this data set with Sharpe's schedule data set. Once this is complete, we will have time series data where each time increment is one game. Breaking this down step-by-step, first load the play-by-play data, group it by game_id, and slice all but the last observation for that game.

```{r}
pbp <- load_pbp(2021:2024) 
performances <- pbp %>% group_by(game_id) %>% slice_tail(n=1)
```

At this point, our 'performances' data set is each game from the 2021-2024 seasons. Now, we will pull from 'performances,' the useful quantitative metrics we would like. Anything that says EPA means expected points added, and anything that says WPA means win percentage added.

```{r}
performances <- performances %>% select(game_id, total_home_epa, total_home_rush_epa, total_home_pass_epa, total_home_comp_air_epa, total_home_raw_air_epa, total_home_comp_yac_epa, total_home_comp_air_wpa, total_home_comp_yac_wpa, total_home_pass_wpa, total_home_raw_air_wpa, total_home_rush_wpa, total_home_raw_yac_epa, total_home_raw_yac_wpa, total_away_comp_air_epa, total_away_comp_air_wpa, total_away_comp_yac_epa, total_away_comp_yac_wpa, total_away_epa, total_away_pass_epa, total_away_pass_wpa, total_away_raw_air_epa, total_away_raw_air_wpa, total_away_raw_yac_epa, total_away_raw_yac_wpa, total_away_rush_epa, total_away_rush_wpa)
```

Now that are 'performances' data has everything we want, we will left merge it with the 'games' data (which is Sharpe's original schedule set), and finally remove any unimportant predictors, such as names and unneeded ID's.

```{r}
games <- games %>% left_join(performances, by = 'game_id')

games <- games %>% select(-old_game_id, -nfl_detail_id, -pfr, -pff, -espn, -ftn, -away_qb_id, -home_qb_id, -away_qb_name, -home_qb_name, -away_coach, -home_coach, -referee, -stadium, -stadium_id, -location)

games$home_win <- ifelse(games$result > 0, 1,0)

games <- games %>% arrange(home_team)
```

Finally, we have our desired data set. 'games' now has each observation as a unique NFL game from the 2021-2024 seasons, with many quantitative measures and void of the non-important predictors.

We then want to get rid of any missing data throughout the dataset. The visualization below can help us identify where and how many NA's are in the dataframe. 
```{r}
vis_miss(games)
```

We can then go ahead and remove the rows with missing values

```{r}
# Specifically check the `home_win` column
sum(is.na(games$home_win))

# Remove Rows with Missing Values
games <- games %>% drop_na()

games <- games %>%
  mutate(across(where(is.numeric), ~ ifelse(is.na(.), mean(., na.rm = TRUE), .)))

games <- games %>%
  mutate(across(where(is.factor), ~ ifelse(is.na(.), "Unknown", .)))
```

# Predictive Modeling

#### Overview of Random Forests

A Random forest is a collection of binary classification trees. A large number of trees are built using random bootstrap samples and subsets of predictor (x) variables. Each of these trees classifies the the data and a majority vote is taken across the forest to reach the final conclusion. 

#### Building the RF Model

First, we want to begin by splitting the data into train and test sets. Then, we split each train and test set into its predictor and target variables. 


```{r, cache = T}
# Split the dataset into training and testing sets
set.seed(666)
train_index <- createDataPartition(games$home_win, p = 0.8, list = FALSE)
train_data <- games[train_index, ]
test_data <- games[-train_index, ]

# Separate predictors and target variable
x_train <- train_data %>% select(-home_win, -result)
y_train <- train_data$home_win
x_test <- test_data %>% select(-home_win, -result)
y_test <- test_data$home_win
```


We can then train the Random Forest model on this data. We will be using 100 trees for our forest. This can be adjusted later if needed. 


```{r, cache = T}

# Train the Random Forest model
rf_model <- randomForest(
  x = x_train,
  y = factor(y_train), # Target variable as factor
  ntree = 100,         # Number of trees
  importance = TRUE,   # Calculate feature importance
  proximity = TRUE     # Proximity measure (optional, useful for visualization)
)

# View the model summary
print(rf_model)
```


As we can see by the model summary, the random forest is a good prediction model for this dataset. The confusion matrix presents extremely low error rates. This can be shown as most game outcomes were predicted correctly (all except 7).

#### Testing and Accuracy

We can then apply this model to the test set. 

```{r, cache = T}
# Predict on the test data
rf_predictions <- predict(rf_model, newdata = x_test)

# Confusion Matrix
conf_matrix <- confusionMatrix(
  factor(rf_predictions), 
  factor(y_test), 
  positive = "1"
)
print(conf_matrix)
```

Again, the confusion matrix shows that most games were corectly predicted to be either win/loss (all except 4). This is reflected in the accuracy score which is .9661. The other statistics also support the conclusion that this is a solid model. 


# Variable Importance Scores

Variable Importance scores measure the average change in node-homogenity, entropy, or DecreaseAccuracy across all nodes associated with a given predictor. This provides a measure of how important the predictor is in the end decision. 

We can calculate these scores for the model to identify the variables that play the largest role in the prediction. 

```{r, cache = T}
# Extract feature importance
importance <- importance(rf_model)
importance_df <- as.data.frame(importance)

# Remove columns with NA names
importance_df <- importance_df[, !is.na(colnames(importance_df))]

# Sort by MeanDecreaseAccuracy
importance_df <- importance_df %>%
  arrange(desc(MeanDecreaseAccuracy))

# Create a data frame from the importance object with feature names
importance_df <- as.data.frame(importance)

# Add feature names as a proper column
importance_df$Feature <- rownames(importance_df)

# Reset row names to avoid confusion
rownames(importance_df) <- NULL

# Limit the number of features to top 20 for better visibility
top_features <- importance_df %>%
  top_n(20, wt = MeanDecreaseAccuracy)

ggplot(top_features, aes(x = reorder(Feature, MeanDecreaseAccuracy), y = MeanDecreaseAccuracy)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  coord_flip() +  # Flip axes for better readability
  labs(
    title = "Top 20 Feature Importance (Random Forest)",
    x = "Features",
    y = "Mean Decrease in Accuracy"
  ) +
  theme_minimal() +
  theme(
    axis.text.y = element_text(size = 10),  # Increase y-axis font size
    axis.title = element_text(size = 12),
    plot.title = element_text(size = 14, hjust = 0.5, face = "bold")
  )
```

The graph above shows the top 20 predictors with the highest variable importance scores. One potential method of improving the model would be to limit the predictors to decrease the complexity of the model without significantly impacting the accuracy. 

# Adjusting the Initial Model

We then build on the previous model and use a more developed method of building and training the random forest. This model uses training control. This involves cross-validation, 5 folds, and random search for hyperparameters. We then use the train function from the caret library. 

```{r, cache = T}
# Ensure x_train and x_test are data frames
x_train_df <- as.data.frame(x_train)
x_test_df <- as.data.frame(x_test)

# Ensure y_train is a factor
y_train <- as.factor(y_train)

# Ensure y_test is a factor
y_test <- as.factor(y_test)

# Define training control
train_control <- trainControl(
  method = "cv",        # Cross-validation
  number = 5,           # Number of folds
  search = "random"     # Random search for hyperparameters
)

# Explicitly use caret's train function
tuned_rf_model <- caret::train(
  x = x_train_df,
  y = y_train,        # Ensure this is a factor
  method = "rf",
  trControl = train_control,
  tuneLength = 5
)

# Print tuned model results
print(tuned_rf_model)
```

As we can see this model increases the accuracy of the previous model to .983. The Graph below shows the variable/feature importance scores. As you can notice, many of the variables with high importance scores are common between this model and the previous one.

```{r}
# Feature importance
varImp <- varImp(tuned_rf_model)

# Plot feature importance
plot(varImp, top = 20, main = "Top 20 Feature Importance (Random Forest)")
```

# References/Further Resources
[https://www.nflfastr.com/index.html](https://www.nflfastr.com/index.html)
[https://nflreadr.nflverse.com/](https://nflreadr.nflverse.com/)
[https://www.nflfastr.com/articles/beginners_guide.html#real-life-example-lets-make-a-win-total-model](https://www.nflfastr.com/articles/beginners_guide.html#real-life-example-lets-make-a-win-total-model)
[www.nflfastr.com/articles/stats_variables.html](https://www.nflfastr.com/articles/stats_variables.html)