#| label: load-libraries-r
#| message: false
#| warning: false

# Load required libraries for data manipulation and visualization
library(tidyverse)      # Core data manipulation (dplyr, ggplot2, etc.)
library(nflfastR)       # NFL play-by-play data
library(nflplotR)       # NFL-specific plotting utilities
library(gt)             # Grammar of tables for nice table formatting
library(gtExtras)       # Additional table styling options
library(patchwork)      # Combining multiple plots

# Set consistent plot theme for all visualizations
theme_set(theme_minimal(base_size = 12))

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

# Load play-by-play data from 2010-2023
# This provides 14 seasons: 5 pre-rule change (2010-2014) and 9 post-rule change (2015-2023)
pbp <- load_pbp(2010:2023)

# Filter for extra point and two-point conversion attempts
# We need both types to compare rates and success probabilities
conversion_attempts <- pbp %>%
  filter(
    # Include plays that are either two-point or extra point attempts
    two_point_attempt == 1 | extra_point_attempt == 1,
    # Exclude plays with missing team information (data quality check)
    !is.na(posteam)
  ) %>%
  mutate(
    # Create readable attempt type label
    attempt_type = case_when(
      two_point_attempt == 1 ~ "Two-Point",
      extra_point_attempt == 1 ~ "Extra Point"
    ),
    # Define success for each attempt type
    # Two-point: Check if conversion result was "success"
    # Extra point: Check if result was "good"
    success = case_when(
      two_point_attempt == 1 ~ two_point_conv_result == "success",
      extra_point_attempt == 1 ~ extra_point_result == "good"
    ),
    # Flag for post-rule change era (2015 onward)
    post_rule_change = season >= 2015
  )

cat("Loaded", nrow(conversion_attempts), "conversion attempts from 2010-2023\n")

#| label: fig-historical-trends-r
#| fig-cap: "Two-point conversion attempt rates and success rates over time. The top panel shows that coaches increased two-point attempts after the 2015 rule change, recognizing the changed expected value calculus. The bottom panel reveals that extra point success dropped dramatically after 2015, while two-point conversion rates remained relatively stable, fundamentally altering the strategic landscape."
#| fig-width: 12
#| fig-height: 8
#| message: false
#| warning: false

# Calculate annual statistics for attempt rates and success rates
# This aggregation will show trends over time
annual_stats <- conversion_attempts %>%
  group_by(season, attempt_type) %>%
  summarise(
    # Count total attempts of each type
    attempts = n(),
    # Calculate success rate (proportion of successful attempts)
    success_rate = mean(success, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  group_by(season) %>%
  mutate(
    # Calculate total attempts across both types
    total_attempts = sum(attempts),
    # Calculate percentage of attempts that were two-point conversions
    # This shows how aggressive coaches were
    pct_of_attempts = attempts / total_attempts
  ) %>%
  ungroup()

# Plot 1: Two-point attempt rate over time
# This reveals strategic evolution in coaching decisions
p1 <- annual_stats %>%
  filter(attempt_type == "Two-Point") %>%
  ggplot(aes(x = season, y = pct_of_attempts)) +
  geom_line(color = "#0077BE", size = 1.2) +
  geom_point(size = 3, color = "#0077BE") +
  # Add vertical line at rule change
  geom_vline(xintercept = 2015, linetype = "dashed", color = "red", alpha = 0.5) +
  annotate("text", x = 2015, y = max(annual_stats$pct_of_attempts[annual_stats$attempt_type == "Two-Point"]),
           label = "Rule Change", vjust = -0.5, color = "red") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 0.1)) +
  scale_x_continuous(breaks = seq(2010, 2023, 2)) +
  labs(
    title = "Two-Point Conversion Attempt Rate",
    subtitle = "Percentage of post-TD attempts that are two-point conversions",
    x = NULL,
    y = "% of Attempts"
  ) +
  theme(plot.title = element_text(face = "bold"))

# Plot 2: Success rates over time for both attempt types
# This shows the impact of the rule change on outcomes
p2 <- annual_stats %>%
  ggplot(aes(x = season, y = success_rate, color = attempt_type)) +
  geom_line(size = 1.2) +
  geom_point(size = 3) +
  # Add vertical line at rule change
  geom_vline(xintercept = 2015, linetype = "dashed", color = "red", alpha = 0.5) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_x_continuous(breaks = seq(2010, 2023, 2)) +
  scale_color_manual(
    values = c("Two-Point" = "#0077BE", "Extra Point" = "#FF6B35"),
    name = "Attempt Type"
  ) +
  labs(
    title = "Success Rates Over Time",
    subtitle = "Extra point success dropped after 2015 rule change",
    x = "Season",
    y = "Success Rate"
  ) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "top"
  )

# Combine plots vertically for easy comparison
p1 / p2 +
  plot_annotation(
    caption = "Data: nflfastR | Rule change: 2015 (XP moved from 2-yard line to 15-yard line)"
  )

#| label: calculate-ev-r
#| message: false
#| warning: false

# Calculate success rates pre and post rule change
# We'll compute expected values and confidence intervals for each era
ev_comparison <- conversion_attempts %>%
  # Create era labels for grouping
  mutate(period = if_else(post_rule_change, "Post-2015", "Pre-2015")) %>%
  # Group by era and attempt type
  group_by(period, attempt_type) %>%
  summarise(
    # Count attempts and successes
    attempts = n(),
    successes = sum(success, na.rm = TRUE),
    # Calculate success rate (proportion of successes)
    success_rate = mean(success, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(
    # Calculate expected value based on attempt type
    # Extra point: success_rate * 1 point
    # Two-point: success_rate * 2 points
    expected_value = case_when(
      attempt_type == "Extra Point" ~ success_rate * 1,
      attempt_type == "Two-Point" ~ success_rate * 2
    ),
    # Calculate standard error using binomial formula: sqrt(p(1-p)/n)
    std_error = sqrt(success_rate * (1 - success_rate) / attempts),
    # Calculate 95% confidence interval bounds
    # Using normal approximation: estimate ± 1.96 * SE
    ev_lower = case_when(
      attempt_type == "Extra Point" ~ (success_rate - 1.96 * std_error) * 1,
      attempt_type == "Two-Point" ~ (success_rate - 1.96 * std_error) * 2
    ),
    ev_upper = case_when(
      attempt_type == "Extra Point" ~ (success_rate + 1.96 * std_error) * 1,
      attempt_type == "Two-Point" ~ (success_rate + 1.96 * std_error) * 2
    )
  )

# Display table with formatting
ev_comparison %>%
  select(period, attempt_type, attempts, success_rate, expected_value) %>%
  arrange(period, desc(attempt_type)) %>%
  gt() %>%
  cols_label(
    period = "Period",
    attempt_type = "Attempt Type",
    attempts = "Attempts",
    success_rate = "Success Rate",
    expected_value = "Expected Value"
  ) %>%
  # Format numbers appropriately
  fmt_number(
    columns = c(success_rate, expected_value),
    decimals = 3
  ) %>%
  fmt_number(
    columns = attempts,
    decimals = 0,
    use_seps = TRUE
  ) %>%
  # Highlight maximum expected value in each period
  tab_style(
    style = cell_fill(color = "#E8F4F8"),
    locations = cells_body(
      columns = expected_value,
      rows = expected_value == max(expected_value, na.rm = TRUE)
    )
  ) %>%
  tab_header(
    title = "Expected Value Comparison",
    subtitle = "Pre-2015 vs Post-2015 Rule Change"
  ) %>%
  tab_source_note(
    source_note = "95% confidence intervals calculated using normal approximation"
  )

#| label: situational-success-r
#| message: false
#| warning: false

# Analyze two-point success by various factors
# Focus on post-2015 era for current relevance
two_pt_situations <- conversion_attempts %>%
  filter(two_point_attempt == 1, season >= 2015) %>%
  mutate(
    # Capture score differential before the conversion attempt
    # This is the score after TD but before the conversion
    score_diff_before = score_differential,
    # Categorize quarter and time situations
    quarter_situation = case_when(
      qtr <= 2 ~ "First Half",
      qtr == 3 ~ "Third Quarter",
      qtr == 4 & game_seconds_remaining > 300 ~ "Early 4th Q",
      qtr == 4 & game_seconds_remaining <= 300 ~ "Late 4th Q",
      TRUE ~ "Overtime"
    ),
    # Classify play type from play description
    # This is imperfect but gives us directional insight
    play_type_cat = case_when(
      grepl("pass", desc, ignore.case = TRUE) ~ "Pass",
      grepl("rush|run", desc, ignore.case = TRUE) ~ "Run",
      TRUE ~ "Other"
    )
  )

# Calculate success by quarter situation
quarter_success <- two_pt_situations %>%
  group_by(quarter_situation) %>%
  summarise(
    attempts = n(),
    success_rate = mean(success, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(desc(attempts))  # Order by sample size

# Success by play type
play_type_success <- two_pt_situations %>%
  filter(play_type_cat != "Other") %>%
  group_by(play_type_cat) %>%
  summarise(
    attempts = n(),
    success_rate = mean(success, na.rm = TRUE),
    .groups = "drop"
  )

# Display quarter situation table
quarter_success %>%
  gt() %>%
  cols_label(
    quarter_situation = "Game Situation",
    attempts = "Attempts",
    success_rate = "Success Rate"
  ) %>%
  fmt_number(
    columns = success_rate,
    decimals = 1,
    scale_by = 100,
    pattern = "{x}%"
  ) %>%
  fmt_number(
    columns = attempts,
    decimals = 0,
    use_seps = TRUE
  ) %>%
  tab_header(
    title = "Two-Point Conversion Success by Game Situation",
    subtitle = "2015-2023 Seasons"
  ) %>%
  tab_source_note(
    source_note = "Data: nflfastR"
  )

#| label: two-pt-model-r
#| message: false
#| warning: false

# Prepare modeling data with relevant features
model_data <- conversion_attempts %>%
  filter(
    two_point_attempt == 1,     # Only two-point conversion attempts
    season >= 2015,               # Post-rule change era
    !is.na(success)              # Must have known outcome
  ) %>%
  mutate(
    # Create binary indicator for home attempts
    home_attempt = if_else(posteam == home_team, 1, 0),
    # Score differential (positive = leading, negative = trailing)
    score_diff = score_differential,
    # Late game indicator (4th quarter, <5 minutes)
    is_late_game = if_else(qtr == 4 & game_seconds_remaining < 300, 1, 0),
    # Simplify roof types to three categories
    roof_type = case_when(
      roof == "outdoors" ~ "Outdoor",
      roof == "dome" ~ "Dome",
      TRUE ~ "Retractable"
    )
  ) %>%
  # Select modeling variables and remove any rows with missing data
  select(success, home_attempt, score_diff, is_late_game, qtr, roof_type) %>%
  filter(complete.cases(.))

# Fit logistic regression model
# Logistic regression models binary outcomes (success/failure)
# It estimates log-odds of success as a linear function of predictors
model <- glm(
  success ~ home_attempt + score_diff + is_late_game + qtr + roof_type,
  data = model_data,
  family = binomial(link = "logit")  # Logistic regression specification
)

# Display model summary
summary(model)

# Generate predicted probabilities for each attempt
# type = "response" converts log-odds to probabilities
model_data$predicted_prob <- predict(model, type = "response")

# Evaluate model performance
library(pROC)

# ROC curve and AUC measure discrimination ability
# AUC of 1.0 = perfect predictions, 0.5 = no better than random
roc_obj <- roc(model_data$success, model_data$predicted_prob)
auc_value <- auc(roc_obj)

# Brier score measures calibration (lower is better)
# It's the mean squared error of predicted probabilities
brier_score <- mean((model_data$success - model_data$predicted_prob)^2)

cat("\nModel Performance:\n")
cat("AUC:", round(auc_value, 3), "\n")
cat("Brier Score:", round(brier_score, 4), "\n")

# Create interpretable coefficient table
coef_df <- data.frame(
  variable = names(coef(model)),
  coefficient = coef(model),
  odds_ratio = exp(coef(model)),  # Exponentiate to get odds ratios
  p_value = summary(model)$coefficients[, 4]
) %>%
  filter(variable != "(Intercept)")

# Display coefficient table
coef_df %>%
  gt() %>%
  cols_label(
    variable = "Variable",
    coefficient = "Coefficient",
    odds_ratio = "Odds Ratio",
    p_value = "P-value"
  ) %>%
  fmt_number(
    columns = c(coefficient, odds_ratio),
    decimals = 3
  ) %>%
  fmt_number(
    columns = p_value,
    decimals = 4
  ) %>%
  tab_header(
    title = "Two-Point Conversion Success Model",
    subtitle = "Logistic Regression Coefficients"
  )

#| label: two-point-chart-r
#| message: false
#| warning: false
#| fig-width: 10
#| fig-height: 8

# Create two-point decision chart for trailing scenarios
# This function encodes analytical guidance for when to go for two

create_two_point_chart <- function(p_2pt = 0.48, p_xp = 0.94) {
  # Score differentials to consider (when trailing before scoring TD)
  score_diffs <- 1:14

  # Calculate optimal decision for each scenario
  decisions <- tibble(
    score_diff_after_td = score_diffs - 6,  # After scoring TD, before conversion
    go_for_2 = case_when(
      # Down 8 (now down 2): Go for 2 to tie immediately
      score_diffs == 8 ~ TRUE,
      # Down 14 (now down 8): Go for 2 to set up second 2PT for tie
      score_diffs == 14 ~ TRUE,
      # Down 7 (now down 1): Kick to tie (most risk-averse option)
      score_diffs == 7 ~ FALSE,
      # Down 1-6: Generally kick to take lead or close gap
      score_diffs <= 6 ~ FALSE,
      # Down 9-13: Complex situations, generally kick
      TRUE ~ FALSE
    ),
    # Provide reasoning for each decision
    reasoning = case_when(
      score_diffs == 8 ~ "Go for tie",
      score_diffs == 14 ~ "Set up 2nd 2PT to tie",
      score_diffs == 7 ~ "Take the tie",
      score_diffs <= 6 ~ "Take lead/close gap",
      TRUE ~ "Kick XP"
    )
  )

  return(decisions)
}

# Generate chart data
chart_data <- create_two_point_chart()

# Create visual representation
chart_data %>%
  mutate(
    decision_label = if_else(go_for_2, "GO FOR 2", "KICK XP"),
    score_label = paste0("Down ", -score_diff_after_td)
  ) %>%
  ggplot(aes(x = reorder(score_label, -score_diff_after_td),
             y = 1,
             fill = go_for_2)) +
  geom_tile(color = "white", size = 2) +
  # Add decision text
  geom_text(aes(label = decision_label),
            size = 5, fontface = "bold", color = "white") +
  # Add reasoning text below decision
  geom_text(aes(label = reasoning),
            y = 0.6, size = 3.5, color = "white") +
  scale_fill_manual(
    values = c("TRUE" = "#D50032", "FALSE" = "#0076CE"),
    guide = "none"
  ) +
  labs(
    title = "Two-Point Conversion Decision Chart",
    subtitle = "Optimal decisions when TRAILING after scoring a touchdown",
    x = "Score Differential (After TD, Before Conversion Attempt)",
    y = NULL,
    caption = "Assumes: 48% 2PT success, 94% XP success\nChart represents general strategy; late-game situations require dynamic analysis"
  ) +
  theme_minimal() +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    panel.grid = element_blank(),
    plot.title = element_text(face = "bold", size = 16),
    plot.subtitle = element_text(size = 11),
    axis.text.x = element_text(size = 11, face = "bold")
  ) +
  coord_fixed(ratio = 3)

#| label: dynamic-programming-r
#| message: false
#| warning: false

# Dynamic programming for optimal two-point decisions
# This is a simplified model for demonstration purposes
# A full model would incorporate possession probabilities, opponent scoring rates, etc.

calculate_win_probability_2pt <- function(
  score_diff_after_td,  # Score after TD, before conversion
  time_remaining_seconds,
  p_2pt = 0.48,
  p_xp = 0.94
) {
  # This simplified model uses score differential as proxy for win probability
  # In practice, you'd use a full win probability model incorporating:
  # - Time remaining
  # - Timeouts
  # - Possession dynamics
  # - Team quality
  # - Field position

  # Calculate score after each possible outcome
  score_after_2pt_success <- score_diff_after_td + 2
  score_after_2pt_fail <- score_diff_after_td
  score_after_xp_success <- score_diff_after_td + 1
  score_after_xp_fail <- score_diff_after_td

  # Simple win probability based on score differential
  # Using normal CDF centered at 0 as approximation
  # Positive score differential → higher win probability
  wp_2pt_success <- pnorm(score_after_2pt_success / 7)
  wp_2pt_fail <- pnorm(score_after_2pt_fail / 7)
  wp_xp_success <- pnorm(score_after_xp_success / 7)
  wp_xp_fail <- pnorm(score_after_xp_fail / 7)

  # Calculate expected win probability for each decision
  ev_2pt <- p_2pt * wp_2pt_success + (1 - p_2pt) * wp_2pt_fail
  ev_xp <- p_xp * wp_xp_success + (1 - p_xp) * wp_xp_fail

  return(list(
    ev_2pt = ev_2pt,
    ev_xp = ev_xp,
    optimal = if_else(ev_2pt > ev_xp, "Go for 2", "Kick XP"),
    advantage = ev_2pt - ev_xp
  ))
}

# Create decision matrix for different score/time scenarios
scenarios <- expand_grid(
  score_diff = seq(-8, -1, 1),  # Trailing by 1-8 after TD
  time_remaining = c(600, 300, 120, 60)  # 10:00, 5:00, 2:00, 1:00
) %>%
  rowwise() %>%
  mutate(
    # Calculate optimal decision for each scenario
    result = list(calculate_win_probability_2pt(score_diff, time_remaining)),
    ev_2pt = result$ev_2pt,
    ev_xp = result$ev_xp,
    optimal = result$optimal,
    advantage = result$advantage
  ) %>%
  ungroup() %>%
  mutate(
    # Create readable time labels
    time_label = case_when(
      time_remaining == 600 ~ "10:00",
      time_remaining == 300 ~ "5:00",
      time_remaining == 120 ~ "2:00",
      time_remaining == 60 ~ "1:00"
    )
  )

# Visualize decision matrix as heatmap
scenarios %>%
  ggplot(aes(x = factor(score_diff), y = time_label, fill = advantage)) +
  geom_tile(color = "white", size = 1) +
  # Add optimal decision text to each cell
  geom_text(aes(label = optimal), size = 4, fontweight = "bold") +
  scale_fill_gradient2(
    low = "#0076CE",      # Blue for "kick XP"
    mid = "white",         # White for neutral
    high = "#D50032",      # Red for "go for 2"
    midpoint = 0,
    name = "WP Advantage\n(2PT - XP)"
  ) +
  labs(
    title = "Dynamic Two-Point Conversion Decisions",
    subtitle = "Optimal strategy by score and time remaining (simplified model)",
    x = "Score Differential (After TD, Before Conversion)",
    y = "Time Remaining",
    caption = "Assumes: 48% 2PT success, 94% XP success | Model simplified for illustration"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    legend.position = "right"
  )

#| label: play-calling-analysis-r
#| message: false
#| warning: false
#| fig-width: 10
#| fig-height: 6

# Analyze play-calling on two-point attempts
# We'll classify plays as pass or run and track success rates

two_pt_plays <- conversion_attempts %>%
  filter(
    two_point_attempt == 1,
    season >= 2015
  ) %>%
  mutate(
    # Classify play type from play description text
    # This is imperfect but provides directional insight
    play_type_clean = case_when(
      grepl("pass", desc, ignore.case = TRUE) ~ "Pass",
      grepl("rush|run", desc, ignore.case = TRUE) ~ "Run",
      grepl("scramble", desc, ignore.case = TRUE) ~ "Scramble",
      TRUE ~ "Other"
    )
  ) %>%
  # Focus on standard pass and run plays
  filter(play_type_clean %in% c("Pass", "Run"))

# Calculate annual trends in play-calling and success
play_calling_trends <- two_pt_plays %>%
  group_by(season, play_type_clean) %>%
  summarise(
    attempts = n(),
    success_rate = mean(success, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  group_by(season) %>%
  mutate(
    # Calculate percentage of two-point attempts using each play type
    pct_of_attempts = attempts / sum(attempts)
  )

# Plot 1: Distribution of play types over time
p1 <- play_calling_trends %>%
  ggplot(aes(x = season, y = pct_of_attempts, color = play_type_clean)) +
  geom_line(size = 1.2) +
  geom_point(size = 3) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_color_manual(values = c("Pass" = "#0077BE", "Run" = "#FF6B35")) +
  labs(
    title = "Two-Point Attempt Play Type Distribution",
    x = NULL,
    y = "% of Attempts",
    color = "Play Type"
  ) +
  theme(legend.position = "top")

# Plot 2: Success rates by play type over time
p2 <- play_calling_trends %>%
  ggplot(aes(x = season, y = success_rate, color = play_type_clean)) +
  geom_line(size = 1.2) +
  geom_point(size = 3) +
  # Add reference line at overall 48% success rate
  geom_hline(yintercept = 0.48, linetype = "dashed", alpha = 0.5) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_color_manual(values = c("Pass" = "#0077BE", "Run" = "#FF6B35")) +
  labs(
    title = "Success Rate by Play Type",
    x = "Season",
    y = "Success Rate",
    color = "Play Type"
  ) +
  theme(legend.position = "top")

# Combine plots vertically
p1 / p2 +
  plot_annotation(
    caption = "Data: nflfastR | Dashed line shows overall 48% success rate"
  )