Simple “Dumb” Expected Touchdowns Model

library(tidyverse)
library(nflfastR)

# Define the range of years you want to load (e.g., 2003 to 2023)
years <- 1999:2024

# Initialize an empty list to store the play-by-play data for each year
pbp_list <- list()

# Loop through each year and load the play-by-play data
for (year in years) {
  # Load the data and add the year column
  pbp_list[[as.character(year)]] <- load_pbp(year) %>%
    mutate(year = year)
}

Warning: URL 'https://objects.githubusercontent.com/github-production-release-asset-2e65be/452908115/bed379f4-234b-4ba6-99e8-e2bd2d143a44?X-Amz-Algorithm=AWS4-HMAC-SHA256&X-Amz-Credential=releaseassetproduction%2F20240919%2Fus-east-1%2Fs3%2Faws4_request&X-Amz-Date=20240919T190029Z&X-Amz-Expires=300&X-Amz-Signature=848e0a3834a16aca63449de978e697267d6774dee764a885fee276519c2df0a8&X-Amz-SignedHeaders=host&response-content-disposition=attachment%3B%20filename%3Dplay_by_play_2017.rds&response-content-type=application%2Foctet-stream': status was 'Failure when receiving data from the peer'Warning: Failed to readRDS from <]8;;https://github.com/nflverse/nflverse-data/releases/download/pbp/play_by_play_2017.rdshttps://github.com/nflverse/nflverse-data/releases/download/pbp/play_by_play_2017.rds]8;;>

# Combine all the yearly play-by-play data into one data frame using bind_rows
pbp_data <- bind_rows(pbp_list)

let’s build an easy expected model based only on yards to go.

first, let’s look at all the different play types

play_type_counts <- pbp_data %>%
  group_by(play_type) %>%
  summarise(count = n()) %>%
  arrange(desc(count))

play_type_counts

let’s looks at the playtypes by year:

# Summarize the number of plays by play_type and year
play_type_year_summary <- pbp_data %>%
  group_by(play_type, year) %>%
  summarise(count = n(), .groups = 'drop')


# Reshape the data to have years as columns
play_type_year_pivot <- play_type_year_summary %>%
  pivot_wider(names_from = year, values_from = count, values_fill = 0)

print(play_type_year_pivot)

Let’s quickly look at a visual to see if there are any evident trends

library(ggplot2)
library(ggrepel)

# Pivot the data to long format so that we can plot it easily with ggplot2
play_type_long <- play_type_year_pivot %>%
  pivot_longer(cols = -play_type, names_to = "year", values_to = "count") %>%
  mutate(year = as.numeric(year))  # Convert year to numeric for the plot

# Create the line plot
ggplot(play_type_long, aes(x = year, y = count, color = play_type)) +
  geom_line(size = 1) +  # Line size for better visibility
  labs(title = "Number of Plays by Play Type Over the Years",
       x = "Year",
       y = "Number of Plays",
       color = "Play Type") +  # Labels and title
  theme_minimal() +  # Clean theme
  theme(
    legend.position = "right",  # Place legend on the right
    plot.title = element_text(hjust = 0.5)  # Center the title
  )

Now let’s look at the

# Summarize the number of plays by play_type and year
touchdown_year_summary <- pbp_data %>%
  filter(play_type == "pass" | play_type =="run" | play_type== "punt" | play_type=="kickoff") %>%
  group_by(yardline_100, touchdown) %>%
  summarise(count = n(), .groups = 'drop')

# Reshape the data to have years as columns
play_touchdown_year_pivot <- touchdown_year_summary %>%
  pivot_wider(names_from = yardline_100, values_from = count, values_fill = 0)

play_touchdown_year_pivot <- play_touchdown_year_pivot %>%
  filter(!is.na(touchdown))

print(play_touchdown_year_pivot)

NA

Let’s calculate probabilities

# Calculate the total number of plays and the number of touchdowns for each yardline
probability_data <- play_touchdown_year_pivot %>%
  pivot_longer(cols = -touchdown, names_to = "yardline_100", values_to = "count") %>%
  mutate(yardline_100 = as.numeric(gsub("X", "", yardline_100))) %>%
  group_by(yardline_100) %>%
  summarise(
    total_plays = sum(count),
    touchdowns = sum(count[touchdown == 1]),
    .groups = 'drop'
  ) %>%
  mutate(probability_td = touchdowns / total_plays) %>%
  select(yardline_100, probability_td)

# Create a new dataframe with probabilities for scoring (1) and not scoring (0)
probability_scores <- probability_data %>%
  mutate(probability_no_td = 1 - probability_td) %>%
  select(yardline_100, probability_td, probability_no_td)

# View the resulting probability scores
print(probability_scores)

Let’s join in the probability data

combined_data <- pbp_data %>%
  left_join(probability_scores, by = "yardline_100") 

Next let’s create a “James Expected TD” column

combined_data <- combined_data %>%
  mutate(
    James_xTD = ifelse(play_type %in% c("pass", "run", "punt", "kickoff"), probability_td, NA)
      )

Let’s summarize by player

expected_td_summary <-combined_data %>%
  group_by(fantasy_player_name) %>%
  summarise(count = n(),
            touchdowns = sum(touchdown),
            James_xTD = sum(James_xTD), .groups = 'drop')

expected_td_summary_by_year <-combined_data %>%
  group_by(fantasy_player_name, year) %>%
  summarise(touchdowns = sum(touchdown),
            James_xTD = sum(James_xTD), .groups = 'drop')

Okay, now let’s regress James_xTF on Touchdowns to see how well the model performs:

# Perform linear regression
model <- lm(touchdowns ~ James_xTD, data = expected_td_summary)

# Summary of the model to check the details of regression output
summary(model)

Call:
lm(formula = touchdowns ~ James_xTD, data = expected_td_summary)

Residuals:
    Min      1Q  Median      3Q     Max 
-44.759  -1.454  -1.099  -0.014  72.583 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.086054   0.129503   8.386   <2e-16 ***
James_xTD   0.914747   0.008156 112.149   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.521 on 3025 degrees of freedom
  (611 observations deleted due to missingness)
Multiple R-squared:  0.8061,    Adjusted R-squared:  0.8061 
F-statistic: 1.258e+04 on 1 and 3025 DF,  p-value: < 2.2e-16

Now let’s visualize:

# Plotting the regression line along with the data points
ggplot(expected_td_summary, aes(x = James_xTD, y = touchdowns)) +
  geom_point(aes(color = touchdowns), alpha = 0.5) +  # Data points with semi-transparency
  geom_smooth(method = "lm", se = TRUE, color = "blue") +  # Linear regression line with confidence interval
  labs(title = "Linear Regression of Touchdowns on James_xTD",
       x = "James_xTD",
       y = "Touchdowns") +
  theme_minimal()  # Clean theme

# Load the ggfortify library for diagnostics
library(ggfortify)

# Generate diagnostic plots
autoplot(model)

We can see, the residual variance is not equal….and possibly violating normality assumptions…. let’s try again breaking the data down by year:

# Perform linear regression
model <- lm(touchdowns ~ James_xTD, data = expected_td_summary_by_year)

# Summary of the model to check the details of regression output
summary(model)

Call:
lm(formula = touchdowns ~ James_xTD, data = expected_td_summary_by_year)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.7570  -0.7924  -0.5391   0.4250  18.2975 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.527055   0.021956    24.0   <2e-16 ***
James_xTD   0.845996   0.005701   148.4   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.013 on 11889 degrees of freedom
  (1738 observations deleted due to missingness)
Multiple R-squared:  0.6494,    Adjusted R-squared:  0.6494 
F-statistic: 2.202e+04 on 1 and 11889 DF,  p-value: < 2.2e-16

Now let’s visualize:

# Plotting the regression line along with the data points
ggplot(expected_td_summary_by_year, aes(x = James_xTD, y = touchdowns)) +
  geom_point(aes(color = touchdowns), alpha = 0.5) +  # Data points with semi-transparency
  geom_smooth(method = "lm", se = TRUE, color = "blue") +  # Linear regression line with confidence interval
  labs(title = "Linear Regression of Touchdowns on James_xTD",
       x = "James_xTD",
       y = "Touchdowns") +
  theme_minimal()  # Clean theme

# Load the ggfortify library for diagnostics
library(ggfortify)

# Generate diagnostic plots
autoplot(model)

expected_td_summary_by_year <- expected_td_summary_by_year %>%
  mutate(
   xTD_TD_deficit = touchdowns - James_xTD
  ) %>%
  arrange(desc(xTD_TD_deficit))

expected_td_summary_by_year

expected_td_summary_by_year_2024 <- expected_td_summary_by_year %>%
  filter(year==2024)

expected_td_summary_by_year_2024 <- expected_td_summary_by_year_2024 %>%
  filter(!is.na(xTD_TD_deficit))

expected_td_summary_by_year_2024

Let’s look at just 2024 so far

library(DT)
# Prepare the top 10 and bottom 10 with rounding
top_10 <- expected_td_summary_by_year %>%
  arrange(desc(xTD_TD_deficit)) %>%
  head(10) %>%
  mutate(
    James_xTD = sprintf("%.2f", James_xTD),
    xTD_TD_deficit = sprintf("%.2f", xTD_TD_deficit)
  )

# Prepare and format bottom 10 Unluckiest Players
bottom_10 <- expected_td_summary_by_year %>%
  arrange(xTD_TD_deficit) %>%
  head(10) %>%
  mutate(
    James_xTD = sprintf("%.2f", James_xTD),
    xTD_TD_deficit = sprintf("%.2f", xTD_TD_deficit)
  )

# Prepare and format top 10 Best Seasons xTD
top_10_xTD <- expected_td_summary_by_year %>%
  arrange(desc(James_xTD)) %>%
  head(10) %>%
  mutate(
    James_xTD = sprintf("%.2f", James_xTD),
    xTD_TD_deficit = sprintf("%.2f", xTD_TD_deficit)
  )

# Prepare and format bottom 10 Worst Seasons xTD
bottom_10_xTD <- expected_td_summary_by_year %>%
  arrange(James_xTD) %>%
  head(10) %>%
  mutate(
    James_xTD = sprintf("%.2f", James_xTD),
    xTD_TD_deficit = sprintf("%.2f", xTD_TD_deficit)
  )

Okay, let’s visualize a little better, first lets create some data frames for just want we want to look at:

# Rename columns for all data frames
top_10 <- top_10 %>%
  rename(
    Player = fantasy_player_name,
    Year = year,
    Touchdowns = touchdowns,
    James_xTD = James_xTD,  # Assuming you want to keep this name, otherwise change as needed
    "Expected TD - TD Diff" = xTD_TD_deficit
  )

bottom_10 <- bottom_10 %>%
  rename(
    Player = fantasy_player_name,
    Year = year,
    Touchdowns = touchdowns,
    James_xTD = James_xTD,  # Change as needed
    "Expected TD - TD Diff" = xTD_TD_deficit
  )

top_10_xTD <- top_10_xTD %>%
  rename(
    Player = fantasy_player_name,
    Year = year,
    Touchdowns = touchdowns,
    James_xTD = James_xTD,  # Change as needed
    "Expected TD - TD Diff" = xTD_TD_deficit
  )

let’s rename the columns

# Create a datatable for the Luckiest Players
luckiest_table <- datatable(
  top_10,  # Assuming 'top_10' contains the top 10 Luckiest players
  options = list(
    pageLength = 10,
    dom = 't',  # Keep the table only; no controls
    ordering = FALSE
  ),
  caption = htmltools::tags$caption(
    style = 'caption-side: top; font-weight: bold; text-align: center; font-size: 16px;',
    'Top 10 Luckiest Players since 1999')
) %>%
  formatStyle(
    columns = c('Player', 'Year', 'Touchdowns', 'James_xTD', 'Expected TD - TD Diff'),
    fontWeight = 'normal',
    fontSize = '14px',
    backgroundColor = 'transparent',
    borderBottom = '1px solid #DDD'
  )

# Create a datatable for the Unluckiest Players
unluckiest_table <- datatable(
  bottom_10,  # Assuming 'bottom_10' contains the bottom 10 Unluckiest players
  options = list(
    pageLength = 10,
    dom = 't',  # Keep the table only; no controls
    ordering = FALSE
  ),
  caption = htmltools::tags$caption(
    style = 'caption-side: top; font-weight: bold; text-align: center; font-size: 16px;',
    'Top 10 Unluckiest Players Since 1999')
) %>%
  formatStyle(
    columns = c('Player', 'Year', 'Touchdowns', 'James_xTD', 'Expected TD - TD Diff'),
    fontWeight = 'normal',
    fontSize = '14px',
    backgroundColor = 'transparent',
    borderBottom = '1px solid #DDD'
  )

#Best Seasons xTD
# Create a datatable for the Unluckiest Players
best_table <- datatable(
  top_10_xTD,  # Assuming 'bottom_10' contains the bottom 10 Unluckiest players
  options = list(
    pageLength = 10,
    dom = 't',  # Keep the table only; no controls
    ordering = FALSE
  ),
  caption = htmltools::tags$caption(
    style = 'caption-side: top; font-weight: bold; text-align: center; font-size: 16px;',
    'Top 10 xTD Seasons by Players Since 1999')
) %>%
  formatStyle(
    columns = c('Player', 'Year', 'Touchdowns', 'James_xTD', 'Expected TD - TD Diff'),
    fontWeight = 'normal',
    fontSize = '14px',
    backgroundColor = 'transparent',
    borderBottom = '1px solid #DDD'
  )

# Render the tables in R Markdown or an R HTML widget
luckiest_table

unluckiest_table

best_table

NA

expected_td_summary_team <- combined_data %>%
  group_by(posteam) %>%
  summarise(
    touchdowns = sum(touchdown, na.rm = TRUE),  # Ensure NAs are ignored
    James_xTD = sum(James_xTD, na.rm = TRUE),   # Ignore NAs
    .groups = 'drop'
  ) %>% 
  mutate(
    xTD_TD_deficit = touchdowns - James_xTD,  # Calculate the deficit
    James_xTD = sprintf("%.2f", James_xTD),   # Format James_xTD to 2 decimal places
    xTD_TD_deficit = sprintf("%.2f", xTD_TD_deficit)  # Format the deficit
  )

# Group by posteam and year, and calculate sum, handling NAs
expected_td_summary_by_year_team <- combined_data %>%
  group_by(posteam, year) %>%
  summarise(
    touchdowns = sum(touchdown, na.rm = TRUE),  # Ensure NAs are ignored
    James_xTD = sum(James_xTD, na.rm = TRUE),   # Ignore NAs
    .groups = 'drop'
  ) %>%
  filter(!is.na(posteam)) %>%  # Ensure posteam is not NA
  mutate(
    xTD_TD_deficit = touchdowns - James_xTD,  # Calculate the deficit
    James_xTD = sprintf("%.2f", James_xTD),   # Format James_xTD to 2 decimal places
    xTD_TD_deficit = sprintf("%.2f", xTD_TD_deficit)  # Format the deficit
  )

# View the result
expected_td_summary_by_year_team

NA

expected_td_summary_by_year_team_2024 <- expected_td_summary_by_year_team %>%
  filter(year==2024) %>%
  arrange(desc(xTD_TD_deficit))

expected_td_summary_by_year_team_2024

NA

expected_td_summary_by_year_team_2024 <- expected_td_summary_by_year_team %>%
  filter(year==2024) %>%
  arrange(desc(xTD_TD_deficit))

expected_td_summary_by_year_team_2024