Chapter 4: Data Wrangling for NFL Analytics

#| label: setup-r
#| message: false
#| warning: false

# Load core tidyverse packages
# The tidyverse is actually a meta-package that loads multiple packages:
# - dplyr (data manipulation)
# - ggplot2 (visualization)
# - tidyr (data reshaping)
# - readr (data import)
# - purrr (functional programming)
# - tibble (modern data frames)
# - stringr (string manipulation)
# - forcats (factor handling)
library(tidyverse)      # Data manipulation and visualization

# Load NFL data packages
# nflfastR provides play-by-play data with advanced metrics pre-calculated
# This includes EPA, WPA, success rate, and 370+ other variables
library(nflfastR)       # NFL play-by-play data

# Load additional packages for specific tasks
# lubridate simplifies date/time operations (parsing, arithmetic, extraction)
library(lubridate)      # Date/time handling

# gt creates publication-quality tables with extensive formatting options
library(gt)             # Professional table formatting

# nflplotR adds NFL team logos, colors, and styling to ggplot2 visualizations
library(nflplotR)       # NFL-specific plotting tools

# Set default ggplot2 theme for consistent appearance
# theme_minimal() provides clean, uncluttered plots suitable for presentations
theme_set(theme_minimal())

# Display confirmation
cat("✓ All R packages loaded successfully\n")
cat("✓ Ready to wrangle football data!\n")

#| label: load-sample-data-r
#| message: false
#| warning: false
#| cache: true

# Load 2023 play-by-play data
# load_pbp() is the primary nflfastR function for accessing play-by-play data
# It accepts a season year or vector of years, e.g., load_pbp(2020:2023)
# The function automatically caches data locally, so subsequent loads are fast
pbp <- load_pbp(2023)

# Load roster data for the 2023 season
# Rosters include player names, positions, height, weight, college, draft info
# This data is essential for player-level analyses and demographics
rosters <- load_rosters(2023)

# Load schedule data
# Schedules include game dates, locations, outcomes, and other game metadata
# This is crucial for aggregating play-level data to game level
schedules <- load_schedules(2023)

# Display summary information about loaded datasets
cat("=" %R% 50, "\n")
cat("LOADED DATA SUMMARY\n")
cat("=" %R% 50, "\n\n")

cat("Play-by-Play Data:\n")
cat(sprintf("  - %s plays (rows)\n", format(nrow(pbp), big.mark = ",")))
cat(sprintf("  - %s variables (columns)\n", ncol(pbp)))
cat(sprintf("  - Memory size: %.1f MB\n",
            object.size(pbp) / 1024^2))
cat(sprintf("  - Weeks: %d through %d\n",
            min(pbp$week, na.rm = TRUE), max(pbp$week, na.rm = TRUE)))

cat("\nRoster Data:\n")
cat(sprintf("  - %s player-season records (rows)\n",
            format(nrow(rosters), big.mark = ",")))
cat(sprintf("  - %s variables (columns)\n", ncol(rosters)))

cat("\nSchedule Data:\n")
cat(sprintf("  - %s games (rows)\n", nrow(schedules)))
cat(sprintf("  - %s variables (columns)\n", ncol(schedules)))

# Display sample of key variables from play-by-play data
cat("\n" %R% "=" %R% 50, "\n")
cat("SAMPLE PLAY-BY-PLAY DATA (First 3 Plays)\n")
cat("=" %R% 50, "\n")
pbp %>%
  select(game_id, posteam, defteam, down, ydstogo, yardline_100,
         desc, yards_gained, epa) %>%
  head(3) %>%
  print(width = 120)

#| label: select-columns-r

# =====================================
# Basic Column Selection
# =====================================

# Select specific columns by name
# This is the most explicit and self-documenting approach
# Use this when you know exactly which columns you need
pbp_basic <- pbp %>%
  select(
    # Identifiers
    game_id, play_id, game_date,

    # Teams
    posteam, defteam,

    # Situation
    down, ydstogo, yardline_100,

    # Outcome
    play_type, yards_gained,

    # Advanced metrics
    epa, wpa, success
  )

cat("Basic selection:", ncol(pbp_basic), "columns\n")

# =====================================
# Pattern-Based Selection
# =====================================

# Select columns matching patterns
# This is useful when you want all columns related to a specific concept
# The tidyverse provides helper functions for pattern matching:
# - starts_with(): columns that start with a prefix
# - ends_with(): columns that end with a suffix
# - contains(): columns that contain a string
# - matches(): columns matching a regular expression
pbp_pass <- pbp %>%
  select(
    game_id, play_id,                # Always include identifiers
    starts_with("pass"),              # All passing-related columns
    ends_with("_epa"),                # All EPA variants
    contains("air_yards")             # All air yards columns
  )

cat("Pass-focused selection:", ncol(pbp_pass), "columns\n")

# =====================================
# Selection by Range
# =====================================

# Select consecutive columns
# Use the colon operator to select all columns between two endpoints
# This is less common because column order might not be meaningful
pbp_range <- pbp %>%
  select(game_id:play_id, down:ydstogo, epa:wpa)

cat("Range selection:", ncol(pbp_range), "columns\n")

# =====================================
# Negative Selection (Dropping Columns)
# =====================================

# Sometimes it's easier to specify what to remove rather than what to keep
# Use the minus sign or - operator to drop columns
pbp_no_ids <- pbp %>%
  select(
    -old_game_id,           # Drop individual columns
    -ends_with("_player_id"), # Drop pattern-matched columns
    -starts_with("td_")     # Drop all touchdown probability columns
  )

cat("After dropping columns:", ncol(pbp_no_ids), "columns\n")

# =====================================
# Selection with Renaming
# =====================================

# Select and rename simultaneously for clearer variable names
pbp_renamed <- pbp %>%
  select(
    game = game_id,
    offense = posteam,
    defense = defteam,
    down = down,
    distance = ydstogo,
    field_position = yardline_100,
    play_description = desc,
    yards = yards_gained,
    expected_points_added = epa
  )

cat("Selection with renaming:", ncol(pbp_renamed), "columns\n")

# Display column names from each selection approach
cat("\n========================================\n")
cat("COLUMN NAMES FROM BASIC SELECTION:\n")
cat("========================================\n")
cat(paste(names(pbp_basic), collapse = ", "), "\n")

cat("\n========================================\n")
cat("SAMPLE OF PASS-FOCUSED COLUMNS:\n")
cat("========================================\n")
cat(paste(head(names(pbp_pass), 15), collapse = ", "), "\n")
cat("... and", ncol(pbp_pass) - 15, "more columns\n")

#| label: filter-basic-r

# =====================================
# Single-Condition Filtering
# =====================================

# Filter to only pass plays
# The == operator tests for equality
# This creates a logical vector (TRUE/FALSE for each row)
# filter() keeps rows where the condition is TRUE
pass_plays <- pbp %>%
  filter(play_type == "pass")

cat("Total plays:", format(nrow(pbp), big.mark = ","), "\n")
cat("Pass plays:", format(nrow(pass_plays), big.mark = ","),
    sprintf("(%.1f%%)\n", 100 * nrow(pass_plays) / nrow(pbp)))

# =====================================
# Multiple-Condition Filtering (AND logic)
# =====================================

# Filter to passes on 3rd down
# Multiple conditions in filter() are combined with AND
# A row must satisfy ALL conditions to be kept
third_down_passes <- pbp %>%
  filter(
    play_type == "pass",    # Must be a pass AND
    down == 3               # Must be 3rd down
  )

cat("\n3rd down passes:", format(nrow(third_down_passes), big.mark = ","), "\n")

# =====================================
# Complex Filtering with Multiple Operators
# =====================================

# Filter to critical situations: 3rd/4th down, medium/long distance,
# 4th quarter, close game
# This demonstrates combining multiple logical operators:
# - %in% tests if value is in a set
# - >= tests greater than or equal
# - <= tests less than or equal
critical_plays <- pbp %>%
  filter(
    down %in% c(3, 4),              # 3rd or 4th down
    ydstogo >= 5,                   # 5+ yards to go
    qtr == 4,                       # 4th quarter
    score_differential >= -8,       # Within 8 points (trailing)
    score_differential <= 8         # Within 8 points (leading)
  )

cat("\nCritical situation plays:", format(nrow(critical_plays), big.mark = ","), "\n")

# =====================================
# Filtering with OR Logic
# =====================================

# Filter to plays that are either touchdowns OR explosive gains (20+ yards)
# Use the | operator for OR logic
# A row is kept if it satisfies ANY of the conditions
explosive_plays <- pbp %>%
  filter(
    touchdown == 1 |                # Touchdown, OR
    yards_gained >= 20,             # 20+ yard gain
    !is.na(epa)                     # Must have valid EPA
  )

cat("\nExplosive plays:", format(nrow(explosive_plays), big.mark = ","), "\n")

# =====================================
# Negative Filtering (Exclusion)
# =====================================

# Filter to non-penalty, regular offensive plays with EPA
# Use ! for negation (NOT)
# !is.na(x) means "is not missing"
# ! (condition) means "not condition"
clean_plays <- pbp %>%
  filter(
    !is.na(epa),                    # EPA must not be missing
    play_type %in% c("pass", "run"), # Must be pass or run
    penalty == 0 | is.na(penalty)   # No penalty (or penalty is NA)
  )

cat("\nClean offensive plays:", format(nrow(clean_plays), big.mark = ","), "\n")

# =====================================
# Filtering by Field Position Zones
# =====================================

# Filter to red zone plays (inside opponent's 20-yard line)
# Field position in nflfastR is yards to opponent's goal
# yardline_100 = 20 means 20 yards from opponent's goal (own 80-yard line)
# yardline_100 = 10 means 10 yards from opponent's goal (opponent's 10)
red_zone <- pbp %>%
  filter(
    yardline_100 <= 20,             # 20 or fewer yards to goal
    yardline_100 > 0,               # Not in end zone
    play_type %in% c("pass", "run"),
    !is.na(epa)
  )

cat("\nRed zone plays:", format(nrow(red_zone), big.mark = ","), "\n")

# =====================================
# Filtering by Team
# =====================================

# Filter to Kansas City Chiefs offensive plays
# posteam is the team with possession (offensive team)
chiefs_offense <- pbp %>%
  filter(
    posteam == "KC",
    play_type %in% c("pass", "run"),
    !is.na(epa)
  )

cat("\nKC offensive plays:", format(nrow(chiefs_offense), big.mark = ","), "\n")

# =====================================
# Demonstrate Filter Logic with Counts
# =====================================

cat("\n========================================\n")
cat("FILTERING IMPACT SUMMARY\n")
cat("========================================\n")

filter_summary <- tibble(
  Filter = c(
    "All plays",
    "Pass plays only",
    "3rd down passes",
    "Critical situations",
    "Explosive plays",
    "Clean offensive plays",
    "Red zone plays",
    "KC offensive plays"
  ),
  Rows = c(
    nrow(pbp),
    nrow(pass_plays),
    nrow(third_down_passes),
    nrow(critical_plays),
    nrow(explosive_plays),
    nrow(clean_plays),
    nrow(red_zone),
    nrow(chiefs_offense)
  )
) %>%
  mutate(
    `Percent of Total` = Rows / nrow(pbp) * 100,
    `Rows Kept` = format(Rows, big.mark = ","),
    `% Kept` = sprintf("%.1f%%", `Percent of Total`)
  ) %>%
  select(Filter, `Rows Kept`, `% Kept`)

print(filter_summary)

#| label: red-zone-analysis-r

# =====================================
# Define and Filter Red Zone Plays
# =====================================

# Red zone is defined as inside opponent's 20-yard line
# In nflfastR, yardline_100 represents yards to opponent's goal:
# - yardline_100 = 20 means you're AT your opponent's 20-yard line
# - yardline_100 = 10 means you're at your opponent's 10-yard line
# - yardline_100 = 1 means you're at the 1-yard line
# We exclude yardline_100 = 0 (plays in end zone, usually PATs)

red_zone <- pbp %>%
  filter(
    # Only pass and run plays (exclude special teams)
    play_type %in% c("pass", "run"),

    # Inside the 20-yard line but not in end zone
    yardline_100 <= 20,
    yardline_100 > 0,

    # Only plays with valid EPA
    !is.na(epa)
  )

cat("Red Zone Analysis\n")
cat("=" %R% 50, "\n")
cat("Total plays in red zone:", format(nrow(red_zone), big.mark = ","), "\n")
cat("Percentage of all offensive plays:",
    sprintf("%.1f%%\n", 100 * nrow(red_zone) / sum(pbp$play_type %in% c("pass", "run") & !is.na(pbp$epa))))

# =====================================
# Calculate Efficiency Metrics by Play Type
# =====================================

# For each play type (pass/run), calculate:
# - Total plays (volume)
# - Touchdowns (scoring outcomes)
# - Touchdown rate (efficiency)
# - Average EPA (comprehensive efficiency measure)
# - Success rate (consistency)

red_zone_summary <- red_zone %>%
  group_by(play_type) %>%
  summarise(
    # Count total plays
    plays = n(),

    # Count touchdowns
    # touchdown variable is 1 for TD, 0 otherwise
    # sum() counts the 1s
    touchdowns = sum(touchdown == 1, na.rm = TRUE),

    # Calculate TD rate
    # This is the primary outcome we care about in red zone
    td_rate = touchdowns / plays,

    # Calculate average EPA
    # EPA captures full value including field position changes
    avg_epa = mean(epa, na.rm = TRUE),

    # Calculate success rate
    # Success = positive EPA
    success_rate = mean(epa > 0, na.rm = TRUE),

    # Calculate median EPA for robustness to outliers
    median_epa = median(epa, na.rm = TRUE),

    .groups = "drop"
  ) %>%
  # Calculate percentage of total red zone plays
  mutate(
    pct_of_rz_plays = plays / sum(plays) * 100
  )

# =====================================
# Display Results in Professional Table
# =====================================

red_zone_summary %>%
  # Create publication-quality table
  gt() %>%
  # Add title and context
  tab_header(
    title = "Red Zone Efficiency by Play Type",
    subtitle = "2023 NFL Regular Season | Plays inside opponent's 20-yard line"
  ) %>%
  # Readable column labels
  cols_label(
    play_type = "Play Type",
    plays = "Plays",
    pct_of_rz_plays = "% of RZ Plays",
    touchdowns = "Touchdowns",
    td_rate = "TD Rate",
    avg_epa = "Avg EPA",
    median_epa = "Median EPA",
    success_rate = "Success Rate"
  ) %>%
  # Format numbers appropriately
  fmt_number(
    columns = c(plays, touchdowns),
    decimals = 0,
    use_seps = TRUE
  ) %>%
  fmt_percent(
    columns = c(td_rate, success_rate, pct_of_rz_plays),
    decimals = 1
  ) %>%
  fmt_number(
    columns = c(avg_epa, median_epa),
    decimals = 3
  ) %>%
  # Add source note
  tab_source_note(
    source_note = "Data: nflfastR | Red zone defined as inside opponent's 20-yard line"
  ) %>%
  # Highlight best values
  # Highest values are best for avg_epa, td_rate, success_rate
  data_color(
    columns = c(avg_epa, td_rate, success_rate),
    palette = "Greens"
  )

# =====================================
# Narrative Summary of Findings
# =====================================

cat("\n\n" %R% "=" %R% 50, "\n")
cat("KEY FINDINGS\n")
cat("=" %R% 50, "\n\n")

pass_data <- red_zone_summary %>% filter(play_type == "pass")
run_data <- red_zone_summary %>% filter(play_type == "run")

cat(sprintf("VOLUME:\n"))
cat(sprintf("- Passes: %s plays (%.1f%% of red zone plays)\n",
            format(pass_data$plays, big.mark = ","),
            pass_data$pct_of_rz_plays))
cat(sprintf("- Runs: %s plays (%.1f%% of red zone plays)\n",
            format(run_data$plays, big.mark = ","),
            run_data$pct_of_rz_plays))

cat(sprintf("\nTOUCHDOWN EFFICIENCY:\n"))
cat(sprintf("- Pass TD rate: %.1f%%\n", pass_data$td_rate * 100))
cat(sprintf("- Run TD rate: %.1f%%\n", run_data$td_rate * 100))
cat(sprintf("- Advantage: %s by %.1f percentage points\n",
            ifelse(pass_data$td_rate > run_data$td_rate, "Pass", "Run"),
            abs(pass_data$td_rate - run_data$td_rate) * 100))

cat(sprintf("\nEPA EFFICIENCY:\n"))
cat(sprintf("- Pass avg EPA: %.3f\n", pass_data$avg_epa))
cat(sprintf("- Run avg EPA: %.3f\n", run_data$avg_epa))
cat(sprintf("- Advantage: %s by %.3f EPA\n",
            ifelse(pass_data$avg_epa > run_data$avg_epa, "Pass", "Run"),
            abs(pass_data$avg_epa - run_data$avg_epa)))

cat(sprintf("\nCONSISTENCY (Success Rate):\n"))
cat(sprintf("- Pass success rate: %.1f%%\n", pass_data$success_rate * 100))
cat(sprintf("- Run success rate: %.1f%%\n", run_data$success_rate * 100))

#| label: mutate-basic-r

# =====================================
# Create Simple Derived Variables
# =====================================

# mutate() creates new variables while keeping all existing ones
# Each line inside mutate() creates one new variable
# Variables can reference previously created variables in the same mutate()

pbp_enhanced <- pbp %>%
  filter(play_type %in% c("pass", "run"), !is.na(epa)) %>%
  mutate(
    # ===== Simple Transformations =====

    # Binary indicator for successful plays (EPA > 0)
    # Logical comparisons (epa > 0) create TRUE/FALSE values
    # You can use these directly or convert to 0/1
    is_success = epa > 0,

    # Alternative: convert to numeric 0/1
    success_numeric = as.numeric(epa > 0),

    # Absolute value of EPA (magnitude regardless of direction)
    epa_magnitude = abs(epa),

    # EPA relative to league average (EPA over expectation)
    # mean(epa, na.rm = TRUE) calculates league average
    epa_over_avg = epa - mean(epa, na.rm = TRUE),

    # ===== Categorizing Continuous Variables =====

    # Down as descriptive text
    # case_when() is like a series of if-else statements
    # Each line is: condition ~ result
    # The first TRUE condition determines the result
    down_label = case_when(
      down == 1 ~ "First Down",
      down == 2 ~ "Second Down",
      down == 3 ~ "Third Down",
      down == 4 ~ "Fourth Down",
      TRUE ~ "Unknown"  # Default/catch-all case
    ),

    # Distance category (short/medium/long yardage)
    # Short: 1-3 yards (power situations)
    # Medium: 4-7 yards (standard situations)
    # Long: 8+ yards (obvious passing situations)
    distance_cat = case_when(
      ydstogo <= 3 ~ "Short (1-3)",
      ydstogo <= 7 ~ "Medium (4-7)",
      ydstogo > 7 ~ "Long (8+)",
      TRUE ~ "Unknown"
    ),

    # Field position zones
    # Field position dramatically affects strategy and EPA
    field_zone = case_when(
      yardline_100 > 80 ~ "Own Territory (1-20)",
      yardline_100 > 50 ~ "Own Side (21-50)",
      yardline_100 > 20 ~ "Opponent Side (31-80)",
      yardline_100 > 0 ~ "Red Zone (1-20)",
      TRUE ~ "Unknown"
    ),

    # ===== Creating Binary Flags =====

    # Flag explosive plays (20+ yards gained)
    # These are high-value plays that dramatically shift field position
    is_explosive = yards_gained >= 20,

    # Flag chunk plays (10+ yards)
    # These advance the offense significantly
    is_chunk = yards_gained >= 10,

    # Flag negative plays (lost yardage)
    # These hurt offense and often create difficult down-distance situations
    is_negative = yards_gained < 0,

    # Flag scoring plays (touchdowns)
    is_touchdown = touchdown == 1,

    # ===== Time-Based Variables =====

    # Convert game_seconds_remaining to minutes for readability
    game_minutes_remaining = game_seconds_remaining / 60,

    # Flag two-minute drill situations (less than 2 minutes in half)
    two_minute_drill = game_seconds_remaining <= 120,

    # ===== Multi-Variable Combinations =====

    # Yards per attempt (could be passing or rushing)
    # This normalizes yards_gained to per-play basis
    ypa = yards_gained,  # In this context, it's the same as yards_gained

    # Yards needed percentage (what % of needed yards were gained)
    # A 5-yard gain on 3rd-and-10 is 50% of needed
    # A 5-yard gain on 3rd-and-3 is 167% of needed (conversion)
    yards_needed_pct = if_else(
      ydstogo > 0,
      yards_gained / ydstogo,
      NA_real_  # Avoid division by zero
    )
  )

# =====================================
# Display Sample of New Variables
# =====================================

cat("=" %R% 50, "\n")
cat("SAMPLE OF NEWLY CREATED VARIABLES\n")
cat("=" %R% 50, "\n\n")

# Show a few plays with original and derived variables
pbp_enhanced %>%
  select(
    # Original variables
    down, ydstogo, yardline_100, yards_gained, epa,
    # Derived variables
    down_label, distance_cat, field_zone,
    is_success, is_explosive, is_chunk
  ) %>%
  head(5) %>%
  print(width = 120)

# =====================================
# Summary Statistics for Derived Variables
# =====================================

cat("\n" %R% "=" %R% 50, "\n")
cat("SUMMARY OF DERIVED VARIABLE DISTRIBUTIONS\n")
cat("=" %R% 50, "\n\n")

# For categorical variables, show counts
cat("Distance Categories:\n")
pbp_enhanced %>%
  count(distance_cat) %>%
  mutate(pct = n / sum(n) * 100) %>%
  arrange(desc(n)) %>%
  print()

cat("\nField Position Zones:\n")
pbp_enhanced %>%
  count(field_zone) %>%
  mutate(pct = n / sum(n) * 100) %>%
  arrange(desc(n)) %>%
  print()

# For binary variables, show TRUE percentages
cat("\nBinary Indicators:\n")
binary_summary <- tibble(
  Variable = c("Success Rate", "Explosive Play Rate",
               "Chunk Play Rate", "Negative Play Rate",
               "Touchdown Rate", "Two-Minute Drill"),
  Percentage = c(
    mean(pbp_enhanced$is_success) * 100,
    mean(pbp_enhanced$is_explosive) * 100,
    mean(pbp_enhanced$is_chunk) * 100,
    mean(pbp_enhanced$is_negative) * 100,
    mean(pbp_enhanced$is_touchdown, na.rm = TRUE) * 100,
    mean(pbp_enhanced$two_minute_drill) * 100
  )
) %>%
  mutate(Percentage = sprintf("%.1f%%", Percentage))

print(binary_summary)