Capstone - Complete Analytics System

# Understanding xG features in StatsBomb data
library(StatsBombR)
library(dplyr)

# Load sample match
matches <- FreeMatches(FreeCompetitions() %>%
                        filter(competition_id == 43, season_id == 106))
events <- get.matchFree(matches[1, ])

# Examine shot features
shots <- events %>%
  filter(type.name == "Shot") %>%
  select(player.name, location.x, location.y,
         shot.statsbomb_xg, shot.body_part.name,
         shot.type.name, shot.technique.name,
         shot.outcome.name) %>%
  mutate(
    # Calculate distance to goal center (120, 40)
    distance = sqrt((120 - location.x)^2 + (40 - location.y)^2),
    # Calculate angle (simplified)
    angle = atan2(8, distance) * 180 / pi  # Goal is 8 yards wide
  )

# View relationship between features and xG
print(shots %>%
        select(player.name, distance, angle, shot.statsbomb_xg,
               shot.body_part.name, shot.outcome.name) %>%
        arrange(desc(shot.statsbomb_xg)))

# Average xG by body part
print(shots %>%
        group_by(shot.body_part.name) %>%
        summarise(
          shots = n(),
          avg_xG = mean(shot.statsbomb_xg, na.rm = TRUE),
          goals = sum(shot.outcome.name == "Goal")
        ))

# Visualize xG distribution
library(ggplot2)

# Load more matches for better distribution
events <- free_allevents(MatchesDF = matches[1:10, ])
shots <- events %>% filter(type.name == "Shot")

# xG histogram
ggplot(shots, aes(x = shot.statsbomb_xg)) +
  geom_histogram(bins = 50, fill = "#1B5E20", color = "white", alpha = 0.8) +
  geom_vline(xintercept = mean(shots$shot.statsbomb_xg, na.rm = TRUE),
             color = "red", linetype = "dashed", size = 1) +
  annotate("text", x = 0.2, y = 150,
           label = paste("Mean xG:", round(mean(shots$shot.statsbomb_xg, na.rm = TRUE), 3)),
           color = "red") +
  labs(title = "Distribution of Expected Goals",
       subtitle = "Most shots have low xG; high-quality chances are rare",
       x = "xG Value", y = "Number of Shots") +
  theme_minimal() +
  scale_x_continuous(breaks = seq(0, 1, 0.1))

# xG categories
shots %>%
  mutate(xg_category = case_when(
    shot.statsbomb_xg >= 0.5 ~ "Big Chance (0.5+)",
    shot.statsbomb_xg >= 0.2 ~ "Good Chance (0.2-0.5)",
    shot.statsbomb_xg >= 0.1 ~ "Reasonable (0.1-0.2)",
    shot.statsbomb_xg >= 0.05 ~ "Low Quality (0.05-0.1)",
    TRUE ~ "Very Low (<0.05)"
  )) %>%
  group_by(xg_category) %>%
  summarise(
    shots = n(),
    pct = n() / nrow(shots) * 100,
    avg_conversion = mean(shot.outcome.name == "Goal") * 100
  ) %>%
  arrange(desc(pct))

# Working with StatsBomb xG
library(StatsBombR)
library(dplyr)

# Load competition data
comps <- FreeCompetitions() %>%
  filter(competition_id == 43, season_id == 106)  # World Cup 2022
matches <- FreeMatches(comps)
events <- free_allevents(MatchesDF = matches)

# Extract all shots with xG
shots <- events %>%
  filter(type.name == "Shot") %>%
  select(match_id, team.name, player.name, minute,
         location.x, location.y,
         shot.statsbomb_xg, shot.outcome.name,
         shot.body_part.name, shot.type.name,
         shot.first_time, shot.one_on_one)

# Team xG totals
team_xg <- shots %>%
  group_by(team.name) %>%
  summarise(
    matches = n_distinct(match_id),
    shots = n(),
    goals = sum(shot.outcome.name == "Goal"),
    xG = sum(shot.statsbomb_xg, na.rm = TRUE),
    npxG = sum(shot.statsbomb_xg[shot.type.name != "Penalty"], na.rm = TRUE)
  ) %>%
  mutate(
    xG_per_match = round(xG / matches, 2),
    goals_minus_xG = goals - xG,
    conversion = round(goals / shots * 100, 1)
  ) %>%
  arrange(desc(xG))

print("World Cup 2022 Team xG:")
print(head(team_xg, 10))

# Player xG leaders
player_xg <- shots %>%
  group_by(player.name, team.name) %>%
  summarise(
    shots = n(),
    goals = sum(shot.outcome.name == "Goal"),
    xG = sum(shot.statsbomb_xg, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(goals_minus_xG = goals - xG) %>%
  arrange(desc(xG))

print("\nTop 10 Players by xG:")
print(head(player_xg, 10))

# Working with Understat xG via understatr
library(understatr)
library(dplyr)

# Get team season data
epl_teams <- get_league_teams_stats(league_name = "EPL", year = 2023)

# View team xG data
team_xg <- epl_teams %>%
  select(team_name, matches, scored, missed, xG, xGA, xpts, pts) %>%
  mutate(
    goals_minus_xG = scored - xG,
    conceded_minus_xGA = missed - xGA,
    pts_minus_xpts = pts - xpts
  ) %>%
  arrange(desc(xG))

print("EPL 2023-24 Team xG:")
print(team_xg)

# Get player data
player_data <- get_league_players_stats(league_name = "EPL", year = 2023)

# Top scorers by xG
top_xg <- player_data %>%
  select(player_name, team_name, games, goals, xG, shots) %>%
  mutate(
    goals_minus_xG = goals - xG,
    xG_per_shot = xG / shots
  ) %>%
  arrange(desc(xG)) %>%
  head(15)

print("\nTop 15 Players by xG:")
print(top_xg)

# Scraping FBref xG data with worldfootballR
library(worldfootballR)
library(dplyr)

# Get league-wide player stats
epl_stats <- fb_big5_advanced_season_stats(
  season_end_year = 2024,
  stat_type = "shooting",
  team_or_player = "player"
)

# Filter for EPL and analyze xG
epl_players <- epl_stats %>%
  filter(Comp == "Premier League") %>%
  select(Player, Squad, Min_Playing, Gls_Standard, xG_Expected,
         npxG_Expected, Sh_Standard, SoT_Standard) %>%
  mutate(
    nineties = Min_Playing / 90,
    goals_minus_xG = Gls_Standard - xG_Expected,
    xG_per_90 = xG_Expected / nineties,
    shots_per_90 = Sh_Standard / nineties
  ) %>%
  filter(Min_Playing >= 900) %>%  # Minimum 10 matches
  arrange(desc(xG_Expected))

print("EPL Top Scorers by xG (min 900 mins):")
print(head(epl_players, 15))

# Biggest overperformers
print("\nBiggest Overperformers (Goals - xG):")
print(epl_players %>%
        arrange(desc(goals_minus_xG)) %>%
        select(Player, Squad, Gls_Standard, xG_Expected, goals_minus_xG) %>%
        head(10))

# Simulate goal outcomes from xG
library(dplyr)

simulate_goals <- function(xg_values, n_simulations = 10000) {
  # Each shot is a Bernoulli trial with p = xG
  goals_per_sim <- sapply(1:n_simulations, function(i) {
    sum(runif(length(xg_values)) < xg_values)
  })

  return(goals_per_sim)
}

# Example: Team had shots with these xG values
shot_xgs <- c(0.75, 0.35, 0.12, 0.08, 0.05, 0.03, 0.02, 0.02)
total_xg <- sum(shot_xgs)  # 1.42 xG

# Simulate 10,000 times
simulated_goals <- simulate_goals(shot_xgs)

# What percentage of simulations result in each goal count?
goal_distribution <- table(simulated_goals) / 10000 * 100

cat(sprintf("Total xG: %.2f\n", total_xg))
cat("\nGoal distribution from 10,000 simulations:\n")
print(round(goal_distribution, 1))

cat(sprintf("\nMost likely outcome: %d goals (%.1f%%)\n",
            as.numeric(names(which.max(goal_distribution))),
            max(goal_distribution)))
cat(sprintf("Probability of 0 goals: %.1f%%\n",
            goal_distribution["0"]))

# Analyze over/underperformance
library(dplyr)

# Calculate player finishing skill
player_finishing <- shots %>%
  group_by(player.name, team.name) %>%
  summarise(
    shots = n(),
    goals = sum(shot.outcome.name == "Goal"),
    xG = sum(shot.statsbomb_xg, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  filter(shots >= 10) %>%  # Minimum sample size
  mutate(
    goals_minus_xG = goals - xG,
    finishing_pct = (goals - xG) / xG * 100,  # % over/under xG
    conversion = goals / shots * 100,
    xG_per_shot = xG / shots
  )

# Overperformers (potential regression candidates)
cat("\nBiggest Overperformers (may regress):\n")
print(player_finishing %>%
        filter(goals_minus_xG > 0) %>%
        arrange(desc(goals_minus_xG)) %>%
        select(player.name, shots, goals, xG, goals_minus_xG) %>%
        head(5))

# Underperformers (potential bounce-back candidates)
cat("\nBiggest Underperformers (may improve):\n")
print(player_finishing %>%
        filter(goals_minus_xG < 0) %>%
        arrange(goals_minus_xG) %>%
        select(player.name, shots, goals, xG, goals_minus_xG) %>%
        head(5))

# Note: True finishing skill requires multi-season analysis
cat("\nNote: Single tournament data is noisy.")
cat("\nMulti-season analysis needed for true finishing skill.")

# Calculate npxG
player_npxg <- shots %>%
  group_by(player.name, team.name) %>%
  summarise(
    total_shots = n(),
    penalties = sum(shot.type.name == "Penalty"),
    non_pen_shots = sum(shot.type.name != "Penalty"),
    goals = sum(shot.outcome.name == "Goal"),
    pen_goals = sum(shot.type.name == "Penalty" & shot.outcome.name == "Goal"),
    np_goals = goals - pen_goals,
    xG = sum(shot.statsbomb_xg, na.rm = TRUE),
    npxG = sum(shot.statsbomb_xg[shot.type.name != "Penalty"], na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(
    np_goals_minus_npxG = np_goals - npxG
  ) %>%
  filter(non_pen_shots >= 5)

# Compare with and without penalties
print("Player xG vs npxG Comparison:")
print(player_npxg %>%
        arrange(desc(xG)) %>%
        select(player.name, goals, xG, np_goals, npxG, penalties) %>%
        head(10))

# Calculate team xG difference
library(dplyr)

# Get xG for and against per team
team_xg_analysis <- events %>%
  filter(type.name == "Shot") %>%
  group_by(match_id, team.name) %>%
  summarise(xG_for = sum(shot.statsbomb_xg, na.rm = TRUE),
            goals_for = sum(shot.outcome.name == "Goal"),
            .groups = "drop")

# Get opponent xG for each team in each match
match_xg <- team_xg_analysis %>%
  group_by(match_id) %>%
  mutate(
    xG_against = sum(xG_for) - xG_for,
    goals_against = sum(goals_for) - goals_for
  ) %>%
  ungroup()

# Aggregate to team level
team_xgd <- match_xg %>%
  group_by(team.name) %>%
  summarise(
    matches = n(),
    xG_for = sum(xG_for),
    xG_against = sum(xG_against),
    goals_for = sum(goals_for),
    goals_against = sum(goals_against)
  ) %>%
  mutate(
    xGD = xG_for - xG_against,
    actual_GD = goals_for - goals_against,
    xGD_per_match = round(xGD / matches, 2),

    # Performance vs expectation
    goals_vs_xG = goals_for - xG_for,
    conceded_vs_xGA = goals_against - xG_against
  ) %>%
  arrange(desc(xGD))

print("Team xG Difference Rankings:")
print(team_xgd %>%
        select(team.name, matches, xG_for, xG_against, xGD,
               actual_GD, xGD_per_match) %>%
        head(10))

# Calculate expected points from match xG
calculate_xpts <- function(xg_for, xg_against, n_sims = 10000) {
  # Simulate goals using Poisson distribution
  goals_for <- rpois(n_sims, xg_for)
  goals_against <- rpois(n_sims, xg_against)

  # Calculate points: 3 for win, 1 for draw, 0 for loss
  points <- ifelse(goals_for > goals_against, 3,
                   ifelse(goals_for == goals_against, 1, 0))

  # Return expected points and win/draw/loss probabilities
  return(list(
    xPts = mean(points),
    win_prob = mean(goals_for > goals_against),
    draw_prob = mean(goals_for == goals_against),
    loss_prob = mean(goals_for < goals_against)
  ))
}

# Example: A match where Team A had 2.1 xG and Team B had 0.8 xG
result <- calculate_xpts(2.1, 0.8)
cat(sprintf("Team A (2.1 xG vs 0.8 xG):\n"))
cat(sprintf("  Expected Points: %.2f\n", result$xPts))
cat(sprintf("  Win Probability: %.1f%%\n", result$win_prob * 100))
cat(sprintf("  Draw Probability: %.1f%%\n", result$draw_prob * 100))
cat(sprintf("  Loss Probability: %.1f%%\n", result$loss_prob * 100))

# Calculate xPts for all matches
match_xpts <- match_xg %>%
  rowwise() %>%
  mutate(
    xPts = calculate_xpts(xG_for, xG_against)$xPts,
    actual_pts = case_when(
      goals_for > goals_against ~ 3,
      goals_for == goals_against ~ 1,
      TRUE ~ 0
    )
  ) %>%
  ungroup()

# Team xPts totals
team_xpts <- match_xpts %>%
  group_by(team.name) %>%
  summarise(
    matches = n(),
    xPts = sum(xPts),
    actual_pts = sum(actual_pts),
    pts_difference = actual_pts - xPts
  ) %>%
  arrange(desc(xPts))

print("\nTeam Expected Points:")
print(team_xpts)

# Create xG shot map with size encoding
library(ggplot2)
library(ggsoccer)

match_shots <- events %>%
  filter(type.name == "Shot", match_id == matches$match_id[1])

# xG shot map
ggplot(match_shots) +
  annotate_pitch(colour = "white", fill = "#1a472a") +
  geom_point(aes(x = location.x, y = location.y,
                 size = shot.statsbomb_xg,
                 color = shot.outcome.name == "Goal"),
             alpha = 0.8) +
  scale_size_continuous(range = c(2, 12), name = "xG") +
  scale_color_manual(values = c("FALSE" = "#CCCCCC", "TRUE" = "#FFD700"),
                     labels = c("No Goal", "Goal"), name = "Result") +
  coord_flip(xlim = c(60, 120)) +
  theme_pitch() +
  facet_wrap(~team.name, ncol = 2) +
  labs(title = "Match xG Shot Map",
       subtitle = "Point size represents expected goals value") +
  theme(legend.position = "bottom",
        strip.text = element_text(size = 12, face = "bold"))

from statsbombpy import sb
import matplotlib.pyplot as plt
from mplsoccer import VerticalPitch
import numpy as np

# Load World Cup Final
events = sb.events(match_id=3869685)
shots = events[events["type"] == "Shot"].copy()

# Extract coordinates
shots["x"] = shots["location"].apply(lambda loc: loc[0])
shots["y"] = shots["location"].apply(lambda loc: loc[1])
shots["is_goal"] = shots["shot_outcome"] == "Goal"

# Create figure with two half-pitches
fig, axes = plt.subplots(1, 2, figsize=(16, 10))

teams = ["Argentina", "France"]
for idx, team in enumerate(teams):
    pitch = VerticalPitch(
        pitch_type="statsbomb", half=True,
        pitch_color="#1a472a", line_color="white", linewidth=1
    )
    pitch.draw(ax=axes[idx])

    team_shots = shots[shots["team"] == team]

    # Create scatter with xG color gradient
    scatter = pitch.scatter(
        team_shots["x"], team_shots["y"],
        s=team_shots["shot_statsbomb_xg"] * 800 + 100,
        c=team_shots["shot_statsbomb_xg"],
        cmap="RdYlBu_r",
        edgecolors="white",
        linewidth=1.5,
        alpha=0.85,
        ax=axes[idx],
        vmin=0, vmax=0.8
    )

    # Mark goals with stars
    goals = team_shots[team_shots["is_goal"]]
    pitch.scatter(
        goals["x"], goals["y"],
        s=300, marker="*", c="gold",
        edgecolors="black", linewidth=1,
        ax=axes[idx], zorder=5
    )

    # Add xG total
    total_xg = team_shots["shot_statsbomb_xg"].sum()
    goals_scored = team_shots["is_goal"].sum()
    axes[idx].set_title(f"{team}\n{goals_scored} Goals | {total_xg:.2f} xG",
                        color="white", fontsize=14, fontweight="bold", pad=10)

# Add colorbar
cbar = fig.colorbar(scatter, ax=axes, orientation="horizontal",
                    fraction=0.05, pad=0.08, aspect=40)
cbar.set_label("xG Value", color="white", fontsize=12)
cbar.ax.xaxis.set_tick_params(color="white")
plt.setp(plt.getp(cbar.ax.axes, "xticklabels"), color="white")

fig.suptitle("xG Shot Map: World Cup 2022 Final", fontsize=18,
             fontweight="bold", color="white", y=0.98)
fig.patch.set_facecolor("#1a472a")
plt.tight_layout(rect=[0, 0.05, 1, 0.95])
plt.savefig("xg_shot_map.png", dpi=150, bbox_inches="tight", facecolor="#1a472a")
plt.show()

import matplotlib.pyplot as plt
import numpy as np

# Aggregate player data
player_xg = shots.groupby(["player", "team"]).agg(
    shots_count=("shot_statsbomb_xg", "count"),
    goals=("shot_outcome", lambda x: (x == "Goal").sum()),
    xG=("shot_statsbomb_xg", "sum")
).reset_index()

player_xg = player_xg[player_xg["shots_count"] >= 3]

# Create scatter plot
fig, ax = plt.subplots(figsize=(10, 8))

# Identity line
ax.plot([0, 4], [0, 4], "k--", alpha=0.5, linewidth=2, label="Expected")

# Scatter by team
colors = {"Argentina": "#75AADB", "France": "#002654"}
for team in colors:
    team_data = player_xg[player_xg["team"] == team]
    ax.scatter(
        team_data["xG"], team_data["goals"],
        s=team_data["shots_count"] * 30 + 50,
        c=colors[team], alpha=0.7, edgecolors="white",
        linewidth=1.5, label=team
    )

# Add labels for key players
for _, row in player_xg[player_xg["goals"] >= 2].iterrows():
    last_name = row["player"].split()[-1]
    ax.annotate(last_name, (row["xG"], row["goals"]),
                xytext=(5, 5), textcoords="offset points",
                fontsize=10, fontweight="bold")

# Add region labels
ax.text(0.5, 2.8, "Overperforming", color="#4CAF50",
        fontsize=11, fontstyle="italic")
ax.text(2.5, 0.5, "Underperforming", color="#F44336",
        fontsize=11, fontstyle="italic")

ax.set_xlabel("Expected Goals (xG)", fontsize=12)
ax.set_ylabel("Actual Goals", fontsize=12)
ax.set_title("Goals vs xG: World Cup 2022 Final\n" +
             "Points above dashed line = finishing above expectation",
             fontsize=14, fontweight="bold")
ax.legend(loc="upper left")
ax.set_xlim(0, 3.5)
ax.set_ylim(0, 3.5)
ax.set_aspect("equal")
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("xg_vs_goals.png", dpi=150, bbox_inches="tight")
plt.show()

# Simplified version for single match timeline with game phases
import matplotlib.pyplot as plt
import numpy as np

# Calculate cumulative xG for both teams
fig, ax = plt.subplots(figsize=(14, 7))

for team, color in [("Argentina", "#75AADB"), ("France", "#002654")]:
    team_shots = shots[shots["team"] == team].sort_values(["minute", "second"])
    team_shots["cumulative_xG"] = team_shots["shot_statsbomb_xg"].cumsum()

    # Add starting point
    minutes = [0] + team_shots["minute"].tolist()
    cum_xg = [0] + team_shots["cumulative_xG"].tolist()

    ax.step(minutes, cum_xg, where="post", linewidth=2.5,
            color=color, label=f"{team} ({cum_xg[-1]:.2f} xG)", alpha=0.9)

    # Mark goals
    goals = team_shots[team_shots["is_goal"]]
    for _, goal in goals.iterrows():
        ax.scatter(goal["minute"], goal["cumulative_xG"],
                   marker="*", s=400, c=color, edgecolors="gold",
                   linewidth=2, zorder=5)

# Add match phase indicators
ax.axvline(x=45, color="gray", linestyle="--", alpha=0.5, linewidth=1.5)
ax.axvline(x=90, color="gray", linestyle="--", alpha=0.5, linewidth=1.5)
ax.axvline(x=105, color="gray", linestyle=":", alpha=0.5, linewidth=1.5)

ax.text(22.5, ax.get_ylim()[1]*0.95, "1st Half", ha="center",
        fontsize=10, color="gray")
ax.text(67.5, ax.get_ylim()[1]*0.95, "2nd Half", ha="center",
        fontsize=10, color="gray")
ax.text(112.5, ax.get_ylim()[1]*0.95, "ET", ha="center",
        fontsize=10, color="gray")

ax.set_xlabel("Minute", fontsize=12)
ax.set_ylabel("Cumulative xG", fontsize=12)
ax.set_title("Cumulative xG Timeline: World Cup 2022 Final\n" +
             "Stars indicate goals scored",
             fontsize=14, fontweight="bold")
ax.legend(loc="upper left", fontsize=11)
ax.set_xlim(0, 125)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("cumulative_xg_timeline.png", dpi=150, bbox_inches="tight")
plt.show()