Capstone - Complete Analytics System

# Understanding PSxG vs xG
library(dplyr)
library(StatsBombR)

# Load data with freeze frames (for GK position)
matches <- FreeMatches(FreeCompetitions() %>%
                        filter(competition_id == 43))
events <- get.matchFree(matches[1, ])

# StatsBomb freeze frames include shot end location
# We can approximate PSxG concepts

shots <- events %>%
  filter(type.name == "Shot") %>%
  select(player.name, shot.statsbomb_xg, shot.outcome.name,
         shot.end_location.x, shot.end_location.y, shot.end_location.z,
         location.x, location.y)

# Shots on target only have PSxG (off-target = 0)
shots_on_target <- shots %>%
  filter(shot.outcome.name %in% c("Goal", "Saved", "Saved Off Target",
                                   "Saved to Post"))

# Approximate shot quality by end location
# Corners of goal are harder to save
shots_on_target <- shots_on_target %>%
  mutate(
    # Goal is at x=120, y=36-44 (8m wide), z=0-2.44m (8ft high)
    end_y = shot.end_location.y,
    end_z = shot.end_location.z,

    # Distance from goal center
    y_from_center = abs(40 - end_y),  # 40 is goal center
    z_height = ifelse(is.na(end_z), 1, end_z),

    # Corner shots are harder to save
    is_corner = (y_from_center > 2) & (z_height > 1.5),

    # Simple PSxG approximation
    psxg_estimate = case_when(
      is_corner & shot.outcome.name == "Goal" ~ 0.9,
      is_corner ~ 0.6,
      z_height > 2 ~ 0.5,  # High shots
      TRUE ~ shot.statsbomb_xg * 1.2  # Central = closer to pre-shot
    )
  )

print("Shot Placement Analysis:")
print(shots_on_target %>%
        select(player.name, shot.statsbomb_xg, shot.outcome.name,
               y_from_center, z_height, is_corner) %>%
        head(10))

# Calculate Goals Saved Above Expected
library(dplyr)

# For goalkeeper evaluation, we need shots AGAINST their team
# and the PSxG of those shots

# Simulate goalkeeper data (in real analysis, use actual PSxG)
goalkeeper_data <- data.frame(
  goalkeeper = c("Alisson", "Ederson", "De Gea", "Ramsdale", "Pickford"),
  team = c("Liverpool", "Man City", "Man United", "Arsenal", "Everton"),
  shots_faced = c(95, 72, 130, 105, 145),
  saves = c(68, 52, 92, 78, 98),
  goals_conceded = c(27, 20, 38, 27, 47),
  psxg_faced = c(32.5, 24.8, 42.1, 31.2, 52.3)
)

# Calculate GSAx
goalkeeper_data <- goalkeeper_data %>%
  mutate(
    # Goals Saved Above Expected
    GSAx = psxg_faced - goals_conceded,

    # Per 90 (assuming 38 matches)
    GSAx_per_90 = GSAx / 38,

    # Save percentage
    save_pct = saves / shots_faced * 100,

    # PSxG save percentage (saves / shots that should have been saved)
    psxg_save_pct = (shots_faced - goals_conceded) / shots_faced * 100
  ) %>%
  arrange(desc(GSAx))

print("Goalkeeper GSAx Rankings:")
print(goalkeeper_data %>%
        select(goalkeeper, team, goals_conceded, psxg_faced,
               GSAx, GSAx_per_90, save_pct))

# Comprehensive goalkeeper analysis
library(ggplot2)

# Visualize save% vs GSAx
ggplot(goalkeeper_data, aes(x = save_pct, y = GSAx)) +
  geom_point(aes(size = shots_faced), alpha = 0.7) +
  geom_text(aes(label = goalkeeper), vjust = -1, size = 3) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_vline(xintercept = mean(goalkeeper_data$save_pct),
             linetype = "dashed", color = "blue") +
  labs(title = "Save Percentage vs Goals Saved Above Expected",
       subtitle = "Point size = shots faced",
       x = "Save Percentage (%)",
       y = "Goals Saved Above Expected (GSAx)",
       size = "Shots Faced") +
  theme_minimal() +
  annotate("text", x = 75, y = 8, label = "Elite\n(High GSAx, avg save%)",
           size = 3, color = "darkgreen") +
  annotate("text", x = 65, y = -8, label = "Poor\n(Low GSAx, low save%)",
           size = 3, color = "red")

# Shot quality faced comparison
print("\nShot Quality Faced (PSxG per shot):")
goalkeeper_data %>%
  mutate(psxg_per_shot = psxg_faced / shots_faced) %>%
  select(goalkeeper, team, shots_faced, psxg_per_shot, GSAx) %>%
  arrange(desc(psxg_per_shot))

# Analyze finishing skill over multiple seasons
library(dplyr)

# Simulated multi-season data (use actual data in practice)
multi_season <- data.frame(
  player = rep(c("Haaland", "Salah", "Kane", "Son", "Rashford"), each = 3),
  season = rep(c("2021-22", "2022-23", "2023-24"), 5),
  shots = c(82, 120, 110,  # Haaland
            95, 88, 92,    # Salah
            105, 98, 85,   # Kane
            78, 82, 76,    # Son
            88, 95, 72),   # Rashford
  goals = c(22, 36, 32,    # Haaland
            23, 18, 20,    # Salah
            17, 30, 22,    # Kane
            23, 17, 15,    # Son
            12, 17, 8),    # Rashford
  xG = c(18.5, 32.5, 28.2,   # Haaland
         21.8, 17.8, 19.5,   # Salah
         20.2, 25.5, 18.8,   # Kane
         15.2, 14.5, 12.8,   # Son
         14.3, 14.8, 10.2)   # Rashford
)

# Calculate finishing metrics
finishing_analysis <- multi_season %>%
  group_by(player) %>%
  summarise(
    seasons = n(),
    total_shots = sum(shots),
    total_goals = sum(goals),
    total_xG = sum(xG),
    goals_minus_xG = total_goals - total_xG,
    conversion = total_goals / total_shots * 100,
    xG_per_shot = total_xG / total_shots,
    finishing_rate = goals_minus_xG / total_xG * 100  # % over/under xG
  )

# Consistency check: did they beat xG every season?
consistency <- multi_season %>%
  mutate(beat_xG = goals > xG) %>%
  group_by(player) %>%
  summarise(
    seasons_beat_xG = sum(beat_xG),
    avg_over_xG = mean(goals - xG)
  )

finishing_analysis <- finishing_analysis %>%
  left_join(consistency, by = "player") %>%
  arrange(desc(goals_minus_xG))

print("Multi-Season Finishing Analysis:")
print(finishing_analysis)

cat("\nInterpretation:\n")
cat("- Haaland: Elite finisher (consistent overperformance)\n")
cat("- Son: Good finisher (beats xG most seasons)\n")
cat("- Salah: Average finisher (fluctuates around xG)\n")
cat("- Kane: Average (high xG = good positions, not elite finishing)\n")
cat("- Rashford: Below average finisher\n")

# Prepare shot data for xG modeling
library(dplyr)
library(StatsBombR)

# Load multiple matches for training data
comps <- FreeCompetitions() %>%
  filter(competition_id == 43)  # World Cups
matches <- FreeMatches(comps)
events <- free_allevents(MatchesDF = matches[1:50, ])

# Extract shot features
shots <- events %>%
  filter(type.name == "Shot") %>%
  mutate(
    # Target variable
    is_goal = as.numeric(shot.outcome.name == "Goal"),

    # Location features
    x = location.x,
    y = location.y,
    distance = sqrt((120 - x)^2 + (40 - y)^2),
    angle = atan2(8, distance) * 180 / pi,  # Simplified angle

    # Categorical features
    body_part = shot.body_part.name,
    shot_type = shot.type.name,
    technique = shot.technique.name,

    # Binary features
    is_header = as.numeric(body_part == "Head"),
    is_penalty = as.numeric(shot_type == "Penalty"),
    is_first_time = as.numeric(shot.first_time == TRUE),
    is_one_on_one = as.numeric(shot.one_on_one == TRUE),

    # Under pressure
    under_pressure = as.numeric(under_pressure == TRUE)
  ) %>%
  select(is_goal, x, y, distance, angle, is_header, is_penalty,
         is_first_time, is_one_on_one, under_pressure) %>%
  na.omit()

print(sprintf("Training data: %d shots, %.1f%% goals",
              nrow(shots),
              mean(shots$is_goal) * 100))

print("\nFeature summary:")
summary(shots)

# Train xG model using logistic regression
library(caret)

# Split data
set.seed(42)
train_idx <- createDataPartition(shots$is_goal, p = 0.8, list = FALSE)
train_data <- shots[train_idx, ]
test_data <- shots[-train_idx, ]

# Train logistic regression
xg_model <- glm(is_goal ~ distance + angle + is_header + is_penalty +
                  is_first_time + is_one_on_one + under_pressure,
                data = train_data,
                family = binomial(link = "logit"))

# Model summary
print(summary(xg_model))

# Predictions
train_data$xG_pred <- predict(xg_model, train_data, type = "response")
test_data$xG_pred <- predict(xg_model, test_data, type = "response")

# Evaluate: Log loss and calibration
library(MLmetrics)
train_logloss <- LogLoss(train_data$xG_pred, train_data$is_goal)
test_logloss <- LogLoss(test_data$xG_pred, test_data$is_goal)

cat(sprintf("\nTrain Log Loss: %.4f\n", train_logloss))
cat(sprintf("Test Log Loss: %.4f\n", test_logloss))

# Compare to StatsBomb xG (if available)
cat("\nModel captures basic patterns, but commercial models include:")
cat("\n- Defender positions from freeze frames")
cat("\n- Goalkeeper position")
cat("\n- Pass leading to shot characteristics")
cat("\n- Historical player finishing rates")

# Check model calibration
library(ggplot2)

# Bin predictions
test_data$xG_bin <- cut(test_data$xG_pred,
                        breaks = seq(0, 1, 0.1),
                        include.lowest = TRUE)

calibration <- test_data %>%
  group_by(xG_bin) %>%
  summarise(
    n = n(),
    avg_xG = mean(xG_pred),
    actual_rate = mean(is_goal)
  ) %>%
  filter(n >= 5)  # Minimum sample size

# Plot calibration
ggplot(calibration, aes(x = avg_xG, y = actual_rate)) +
  geom_point(aes(size = n)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
  scale_x_continuous(limits = c(0, 1)) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(title = "xG Model Calibration",
       subtitle = "Perfect calibration = points on diagonal line",
       x = "Predicted xG (binned average)",
       y = "Actual Conversion Rate",
       size = "Shots") +
  theme_minimal()

print("Calibration by xG bin:")
print(calibration)