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

library(tidyverse)
library(nflfastR)
library(nflplotR)
library(gt)

# Load team-level data
teams <- load_teams() %>%
  filter(season >= 2015, season <= 2023)

# Load schedules and results
schedules <- load_schedules(2015:2023) %>%
  filter(game_type == "REG")

# Calculate team records and point differential
team_records <- schedules %>%
  filter(!is.na(result)) %>%
  # Home team perspective
  bind_rows(
    schedules %>%
      select(
        season, week, team = home_team, opponent = away_team,
        points = home_score, points_allowed = away_score, result
      ),
    schedules %>%
      select(
        season, week, team = away_team, opponent = home_team,
        points = away_score, points_allowed = home_score
      ) %>%
      mutate(result = -result)
  ) %>%
  filter(!is.na(points)) %>%
  group_by(season, team) %>%
  summarise(
    games = n(),
    wins = sum(result > 0),
    losses = sum(result < 0),
    ties = sum(result == 0),
    points = sum(points),
    points_allowed = sum(points_allowed),
    point_diff = points - points_allowed,
    .groups = "drop"
  )

# Calculate Pythagorean expected wins
# Formula: W_exp = G * (PF^2.37 / (PF^2.37 + PA^2.37))
team_records <- team_records %>%
  mutate(
    pythag_wins = games * (points^2.37 / (points^2.37 + points_allowed^2.37)),
    pythag_diff = wins - pythag_wins,
    win_pct = wins / games
  )

# Display recent examples
team_records %>%
  filter(season == 2023) %>%
  arrange(desc(pythag_diff)) %>%
  select(team, wins, losses, pythag_wins, pythag_diff) %>%
  head(10) %>%
  gt() %>%
  cols_label(
    team = "Team",
    wins = "Wins",
    losses = "Losses",
    pythag_wins = "Expected Wins",
    pythag_diff = "Over/Under"
  ) %>%
  fmt_number(columns = c(pythag_wins, pythag_diff), decimals = 2) %>%
  tab_header(
    title = "Teams Exceeding Pythagorean Expectation",
    subtitle = "2023 Regular Season"
  )

#| label: schedule-strength-r
#| message: false
#| warning: false

# Calculate opponent win percentage for each team
schedule_strength <- schedules %>%
  filter(!is.na(result)) %>%
  # For each game, get opponent's season record
  left_join(
    team_records %>% select(season, team, win_pct),
    by = c("season", "away_team" = "team")
  ) %>%
  rename(away_win_pct = win_pct) %>%
  left_join(
    team_records %>% select(season, team, win_pct),
    by = c("season", "home_team" = "team")
  ) %>%
  rename(home_win_pct = win_pct) %>%
  # Calculate schedule strength for each team
  group_by(season, home_team) %>%
  summarise(
    sos_home = mean(away_win_pct, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  rename(team = home_team) %>%
  bind_rows(
    schedules %>%
      filter(!is.na(result)) %>%
      left_join(
        team_records %>% select(season, team, win_pct),
        by = c("season", "home_team" = "team")
      ) %>%
      rename(home_win_pct = win_pct) %>%
      group_by(season, away_team) %>%
      summarise(
        sos_away = mean(home_win_pct, na.rm = TRUE),
        .groups = "drop"
      ) %>%
      rename(team = away_team)
  ) %>%
  group_by(season, team) %>%
  summarise(
    sos = mean(c(sos_home, sos_away), na.rm = TRUE),
    .groups = "drop"
  )

# Join with team records
team_records_adj <- team_records %>%
  left_join(schedule_strength, by = c("season", "team")) %>%
  mutate(
    # Adjust wins for schedule strength
    # Teams with harder schedules get bonus
    sos_adj_wins = wins + (sos - 0.500) * games,
    sos_diff = sos_adj_wins - wins
  )

# Show 2023 schedule strength
team_records_adj %>%
  filter(season == 2023) %>%
  arrange(desc(sos)) %>%
  select(team, wins, losses, sos, sos_adj_wins) %>%
  head(10) %>%
  gt() %>%
  cols_label(
    team = "Team",
    wins = "Actual Wins",
    losses = "Losses",
    sos = "SOS",
    sos_adj_wins = "Adjusted Wins"
  ) %>%
  fmt_number(columns = c(sos, sos_adj_wins), decimals = 2) %>%
  tab_header(
    title = "Schedule Strength Adjustments",
    subtitle = "2023 Regular Season - Hardest Schedules"
  )

#| label: fig-expected-wins-r
#| fig-cap: "Actual wins vs Pythagorean expected wins by team (2023)"
#| fig-width: 10
#| fig-height: 8
#| message: false
#| warning: false

team_records_adj %>%
  filter(season == 2023) %>%
  ggplot(aes(x = pythag_wins, y = wins)) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "gray40") +
  geom_nfl_logos(aes(team_abbr = team), width = 0.05, alpha = 0.8) +
  scale_x_continuous(breaks = seq(0, 17, 2), limits = c(0, 17)) +
  scale_y_continuous(breaks = seq(0, 17, 2), limits = c(0, 17)) +
  labs(
    title = "Actual Wins vs Pythagorean Expected Wins",
    subtitle = "2023 NFL Regular Season | Teams above line exceeded expectations",
    x = "Pythagorean Expected Wins",
    y = "Actual Wins",
    caption = "Data: nflfastR | Teams above diagonal exceeded expectations based on points scored/allowed"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    plot.subtitle = element_text(size = 11),
    panel.grid.minor = element_blank()
  )

#| label: fourth-down-analysis-r
#| message: false
#| warning: false
#| cache: true

# Load play-by-play data
pbp <- load_pbp(2023)

# Analyze 4th down decisions
fourth_downs <- pbp %>%
  filter(down == 4, !is.na(ydstogo)) %>%
  mutate(
    decision = case_when(
      punt_attempt == 1 ~ "Punt",
      field_goal_attempt == 1 ~ "FG Attempt",
      play_type %in% c("pass", "run") ~ "Go For It",
      TRUE ~ "Other"
    ),
    # Optimal decision based on EPA
    should_go = if_else(go_boost > 0.5, TRUE, FALSE),
    correct_decision = case_when(
      decision == "Go For It" & should_go ~ TRUE,
      decision != "Go For It" & !should_go ~ TRUE,
      TRUE ~ FALSE
    )
  ) %>%
  filter(decision != "Other")

# Calculate decision quality by team
team_4th_down <- fourth_downs %>%
  group_by(posteam) %>%
  summarise(
    fourth_downs = n(),
    go_attempts = sum(decision == "Go For It"),
    go_rate = mean(decision == "Go For It"),
    should_go_situations = sum(should_go),
    correct_decisions = sum(correct_decision),
    decision_rate = mean(correct_decision),
    .groups = "drop"
  ) %>%
  arrange(desc(decision_rate))

# Display top decision-makers
team_4th_down %>%
  head(10) %>%
  gt() %>%
  cols_label(
    posteam = "Team",
    fourth_downs = "4th Downs",
    go_attempts = "Go Attempts",
    go_rate = "Go Rate",
    correct_decisions = "Correct",
    decision_rate = "Correct %"
  ) %>%
  fmt_percent(columns = c(go_rate, decision_rate), decimals = 1) %>%
  tab_header(
    title = "Fourth Down Decision Quality",
    subtitle = "2023 Regular Season - Best Decision Makers"
  )

#| label: challenge-analysis-r
#| message: false
#| warning: false

# Analyze replay challenges
# Note: Challenge data is limited in pbp data
challenges <- pbp %>%
  filter(!is.na(replay_or_challenge), replay_or_challenge == 1) %>%
  mutate(
    challenge_successful = if_else(
      !is.na(replay_or_challenge_result) &
      str_detect(tolower(replay_or_challenge_result), "upheld|confirmed"),
      FALSE,
      TRUE
    )
  )

# Summarize by team (possessing team when challenge occurred)
team_challenges <- challenges %>%
  filter(!is.na(posteam)) %>%
  group_by(posteam) %>%
  summarise(
    challenges = n(),
    successful = sum(challenge_successful, na.rm = TRUE),
    success_rate = mean(challenge_successful, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  filter(challenges >= 3) %>%  # Minimum sample
  arrange(desc(success_rate))

if(nrow(team_challenges) > 0) {
  team_challenges %>%
    head(10) %>%
    gt() %>%
    cols_label(
      posteam = "Team",
      challenges = "Challenges",
      successful = "Successful",
      success_rate = "Success Rate"
    ) %>%
    fmt_percent(columns = success_rate, decimals = 1) %>%
    tab_header(
      title = "Challenge Success Rates",
      subtitle = "2023 Regular Season (min. 3 challenges)"
    )
} else {
  cat("Challenge data not available in sufficient detail\n")
}

#| label: timeout-analysis-r
#| message: false
#| warning: false

# Analyze timeout usage
timeout_analysis <- pbp %>%
  filter(!is.na(timeout), timeout == 1) %>%
  mutate(
    # Categorize timeouts by situation
    timeout_category = case_when(
      qtr <= 2 & game_seconds_remaining < 120 ~ "End of 1st Half",
      qtr == 4 & game_seconds_remaining < 120 ~ "End of Game",
      abs(score_differential) <= 3 ~ "Close Game",
      TRUE ~ "Other"
    ),
    # Was it effective? (led to positive EPA on next play)
    timeout_team = timeout_team
  )

# Count timeouts by team and situation
team_timeouts <- timeout_analysis %>%
  filter(!is.na(timeout_team)) %>%
  group_by(timeout_team, timeout_category) %>%
  summarise(
    timeouts = n(),
    .groups = "drop"
  ) %>%
  pivot_wider(
    names_from = timeout_category,
    values_from = timeouts,
    values_fill = 0
  ) %>%
  mutate(
    total_timeouts = `Close Game` + `End of 1st Half` + `End of Game` + Other
  ) %>%
  arrange(desc(total_timeouts))

# Display
team_timeouts %>%
  head(10) %>%
  gt() %>%
  cols_label(
    timeout_team = "Team",
    `Close Game` = "Close Game",
    `End of 1st Half` = "End Half",
    `End of Game` = "End Game",
    Other = "Other",
    total_timeouts = "Total"
  ) %>%
  tab_header(
    title = "Timeout Usage by Situation",
    subtitle = "2023 Regular Season"
  )

#| label: oc-analysis-r
#| message: false
#| warning: false

# Calculate offensive efficiency metrics
offense_metrics <- pbp %>%
  filter(!is.na(posteam), !is.na(epa), play_type %in% c("pass", "run")) %>%
  group_by(season, posteam) %>%
  summarise(
    plays = n(),
    epa_per_play = mean(epa),
    success_rate = mean(epa > 0),
    explosive_rate = mean(epa > 1.5),
    # Situational efficiency
    third_down_conv = mean(third_down_converted[down == 3], na.rm = TRUE),
    red_zone_td_rate = mean(
      touchdown[yardline_100 <= 20],
      na.rm = TRUE
    ),
    .groups = "drop"
  )

# Add team colors for visualization
offense_metrics <- offense_metrics %>%
  left_join(
    load_teams() %>%
      select(team_abbr, team_color, team_color2) %>%
      distinct(),
    by = c("posteam" = "team_abbr")
  )

# Display 2023 top offenses
offense_metrics %>%
  filter(season == 2023) %>%
  arrange(desc(epa_per_play)) %>%
  select(posteam, plays, epa_per_play, success_rate, explosive_rate,
         third_down_conv, red_zone_td_rate) %>%
  head(10) %>%
  gt() %>%
  cols_label(
    posteam = "Team",
    plays = "Plays",
    epa_per_play = "EPA/Play",
    success_rate = "Success Rate",
    explosive_rate = "Explosive Rate",
    third_down_conv = "3rd Down Conv",
    red_zone_td_rate = "RZ TD Rate"
  ) %>%
  fmt_number(columns = epa_per_play, decimals = 3) %>%
  fmt_percent(columns = c(success_rate, explosive_rate, third_down_conv, red_zone_td_rate),
              decimals = 1) %>%
  tab_header(
    title = "Offensive Coordinator Efficiency",
    subtitle = "2023 Regular Season - Top 10 Offenses by EPA"
  )

#| label: dc-analysis-r
#| message: false
#| warning: false

# Calculate defensive efficiency metrics
defense_metrics <- pbp %>%
  filter(!is.na(defteam), !is.na(epa), play_type %in% c("pass", "run")) %>%
  group_by(season, defteam) %>%
  summarise(
    plays = n(),
    epa_allowed = mean(epa),  # Lower is better
    success_rate_allowed = mean(epa > 0),  # Lower is better
    explosive_rate_allowed = mean(epa > 1.5),  # Lower is better
    # Pressure and turnovers
    pressure_rate = mean(qb_hit == 1 | qb_hurry == 1 | sack == 1, na.rm = TRUE),
    turnover_rate = mean(interception == 1 | fumble_lost == 1, na.rm = TRUE),
    # Situational
    third_down_stop_rate = mean(!third_down_converted[down == 3], na.rm = TRUE),
    red_zone_td_rate = mean(touchdown[yardline_100 <= 20], na.rm = TRUE),
    .groups = "drop"
  )

# Display 2023 top defenses
defense_metrics %>%
  filter(season == 2023) %>%
  arrange(epa_allowed) %>%  # Lower is better
  select(defteam, plays, epa_allowed, success_rate_allowed,
         pressure_rate, turnover_rate, third_down_stop_rate) %>%
  head(10) %>%
  gt() %>%
  cols_label(
    defteam = "Team",
    plays = "Plays",
    epa_allowed = "EPA Allowed",
    success_rate_allowed = "Success Rate",
    pressure_rate = "Pressure Rate",
    turnover_rate = "Turnover Rate",
    third_down_stop_rate = "3rd Down Stops"
  ) %>%
  fmt_number(columns = epa_allowed, decimals = 3) %>%
  fmt_percent(columns = c(success_rate_allowed, pressure_rate, turnover_rate,
                          third_down_stop_rate), decimals = 1) %>%
  tab_header(
    title = "Defensive Coordinator Efficiency",
    subtitle = "2023 Regular Season - Top 10 Defenses by EPA Allowed"
  ) %>%
  tab_style(
    style = cell_fill(color = "lightgreen"),
    locations = cells_body(columns = epa_allowed, rows = epa_allowed < -0.05)
  )

#| label: fig-coordinator-scatter-r
#| fig-cap: "Offensive and defensive EPA by team (2023)"
#| fig-width: 10
#| fig-height: 8
#| message: false
#| warning: false

# Combine offense and defense metrics
combined_metrics <- offense_metrics %>%
  filter(season == 2023) %>%
  select(posteam, off_epa = epa_per_play) %>%
  inner_join(
    defense_metrics %>%
      filter(season == 2023) %>%
      select(defteam, def_epa = epa_allowed),
    by = c("posteam" = "defteam")
  )

# Create scatter plot
combined_metrics %>%
  ggplot(aes(x = off_epa, y = -def_epa)) +  # Negative def_epa so up is better
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_nfl_logos(aes(team_abbr = posteam), width = 0.05, alpha = 0.8) +
  labs(
    title = "Offensive and Defensive Coordinator Efficiency",
    subtitle = "2023 NFL Regular Season | Top-right quadrant = elite on both sides",
    x = "Offensive EPA per Play (Better →)",
    y = "Defensive EPA per Play (Better ↑)",
    caption = "Data: nflfastR | Defensive EPA inverted so up = better defense"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    plot.subtitle = element_text(size = 11),
    panel.grid.minor = element_blank()
  )

#| label: player-development-r
#| message: false
#| warning: false
#| cache: true

# Load multiple seasons
pbp_multi <- load_pbp(2021:2023)

# Calculate QB development (EPA improvement by experience)
qb_development <- pbp_multi %>%
  filter(
    play_type %in% c("pass", "run"),
    !is.na(passer_player_id),
    !is.na(epa)
  ) %>%
  group_by(season, passer_player_id, passer_player_name) %>%
  summarise(
    dropbacks = n(),
    epa_per_play = mean(epa),
    success_rate = mean(epa > 0),
    .groups = "drop"
  ) %>%
  filter(dropbacks >= 100) %>%  # Minimum sample
  # Calculate year-over-year improvement
  group_by(passer_player_id) %>%
  arrange(season) %>%
  mutate(
    seasons_played = row_number(),
    epa_change = epa_per_play - lag(epa_per_play),
    success_change = success_rate - lag(success_rate)
  ) %>%
  ungroup()

# Identify top improvers (2022 to 2023)
top_improvers <- qb_development %>%
  filter(season == 2023, !is.na(epa_change)) %>%
  arrange(desc(epa_change)) %>%
  select(passer_player_name, season, dropbacks, epa_per_play, epa_change)

# Display
top_improvers %>%
  head(10) %>%
  gt() %>%
  cols_label(
    passer_player_name = "Quarterback",
    season = "Season",
    dropbacks = "Dropbacks",
    epa_per_play = "EPA/Play",
    epa_change = "EPA Change"
  ) %>%
  fmt_number(columns = c(epa_per_play, epa_change), decimals = 3) %>%
  tab_header(
    title = "Quarterback Development",
    subtitle = "2023 Season - Biggest YoY Improvements in EPA"
  ) %>%
  tab_style(
    style = cell_fill(color = "lightgreen"),
    locations = cells_body(columns = epa_change, rows = epa_change > 0.05)
  )

#| label: team-development-r
#| message: false
#| warning: false

# Calculate team-level player development
# Focus on young players (first 3 seasons)
young_player_development <- pbp_multi %>%
  filter(
    play_type %in% c("pass", "run"),
    !is.na(epa)
  ) %>%
  # Get passer or rusher
  mutate(
    player_id = coalesce(passer_player_id, rusher_player_id),
    player_name = coalesce(passer_player_name, rusher_player_name)
  ) %>%
  filter(!is.na(player_id)) %>%
  group_by(season, player_id, player_name, posteam) %>%
  summarise(
    plays = n(),
    epa_per_play = mean(epa),
    .groups = "drop"
  ) %>%
  filter(plays >= 50) %>%
  # Calculate improvement
  group_by(player_id) %>%
  arrange(season) %>%
  mutate(
    seasons = n(),
    epa_change = epa_per_play - lag(epa_per_play)
  ) %>%
  filter(seasons >= 2, seasons <= 3, !is.na(epa_change)) %>%
  ungroup()

# Aggregate by team
team_development <- young_player_development %>%
  group_by(posteam) %>%
  summarise(
    players = n(),
    avg_improvement = mean(epa_change),
    positive_improvement_rate = mean(epa_change > 0),
    .groups = "drop"
  ) %>%
  filter(players >= 5) %>%
  arrange(desc(avg_improvement))

# Display
team_development %>%
  head(10) %>%
  gt() %>%
  cols_label(
    posteam = "Team",
    players = "Players",
    avg_improvement = "Avg EPA Improvement",
    positive_improvement_rate = "% Improving"
  ) %>%
  fmt_number(columns = avg_improvement, decimals = 3) %>%
  fmt_percent(columns = positive_improvement_rate, decimals = 1) %>%
  tab_header(
    title = "Team Player Development",
    subtitle = "2021-2023 | Young Players (Years 1-3)"
  )

#| label: halftime-adjustments-r
#| message: false
#| warning: false

# Analyze first half vs second half performance
halftime_analysis <- pbp %>%
  filter(!is.na(epa), play_type %in% c("pass", "run")) %>%
  mutate(
    half = if_else(qtr <= 2, "First Half", "Second Half")
  ) %>%
  group_by(season, posteam, half) %>%
  summarise(
    plays = n(),
    epa_per_play = mean(epa),
    success_rate = mean(epa > 0),
    .groups = "drop"
  ) %>%
  pivot_wider(
    names_from = half,
    values_from = c(plays, epa_per_play, success_rate),
    names_glue = "{half}_"
  ) %>%
  mutate(
    epa_improvement = `Second Half_epa_per_play` - `First Half_epa_per_play`,
    success_improvement = `Second Half_success_rate` - `First Half_success_rate`
  )

# Display 2023 best adjusters
halftime_analysis %>%
  filter(season == 2023) %>%
  arrange(desc(epa_improvement)) %>%
  select(posteam, `First Half_epa_per_play`, `Second Half_epa_per_play`,
         epa_improvement, success_improvement) %>%
  head(10) %>%
  gt() %>%
  cols_label(
    posteam = "Team",
    `First Half_epa_per_play` = "1st Half EPA",
    `Second Half_epa_per_play` = "2nd Half EPA",
    epa_improvement = "EPA Improvement",
    success_improvement = "Success Rate Δ"
  ) %>%
  fmt_number(columns = c(`First Half_epa_per_play`, `Second Half_epa_per_play`,
                         epa_improvement), decimals = 3) %>%
  fmt_number(columns = success_improvement, decimals = 3) %>%
  tab_header(
    title = "Halftime Adjustments",
    subtitle = "2023 Regular Season - Best Second Half Improvements"
  ) %>%
  tab_style(
    style = cell_fill(color = "lightgreen"),
    locations = cells_body(columns = epa_improvement, rows = epa_improvement > 0.05)
  )