Soccer Player Injury Risk Modeling

Introduction

30 Second Summary

Every weekend, coaches make split-second decisions about which players to start, when to substitute, and who needs rest. Behind the scenes at professional clubs, data scientists use statistical models to turn gut feelings into data-driven recommendations.

In this project, you will build a sports analytics system that evaluates player performance using Bayesian hierarchical modeling and predicts injury timing with survival analysis. You'll combine both models into a substitution recommendation system for data-driven squad management.

What You'll Build

A Python-based analytics pipeline that uses PyMC for Bayesian performance estimation with partial pooling across teams, and lifelines for dual survival models predicting when injuries are likely to occur.

By the end of this project, you'll have:

• A Bayesian hierarchical model evaluating player performance with team-level random effects, opponent strength adjustments, and posterior credible intervals for individual players.
• Kaplan-Meier survival curves and dual survival models (Cox Proportional Hazards + Weibull AFT) predicting injury timing with C-index evaluation.
• A substitution recommendation system combining Bayesian performance posteriors with survival predictions to flag high-risk players and support squad decisions.
• Secret Mission: A squad rotation optimizer that recommends which players to rest based on combined injury risk scores and upcoming fixture difficulty.

Are there any prerequisites?

You should be comfortable with basic Python programming and familiar with pandas DataFrames. An understanding of basic statistics (mean, standard deviation, correlation) is expected.

Some familiarity with Bayesian concepts (priors, posteriors) is helpful but not required. The project explains these ideas as you encounter them.

Set Up Your Python Environment

PyMC uses PyTensor for computation, which requires a specific compilation toolchain that conda handles well. The lifelines survival analysis library also requires pandas pinned below version 3.0, so getting the right combination of packages matters.

In this step, you will install Python, VS Code, and conda. Then you will create an isolated environment with all the packages needed for Bayesian modeling and survival analysis, and download the football dataset.

In this step, get ready to:

• Install Python 3.11+, VS Code, and the Jupyter extension.
• Create a conda environment with all required packages.
• Download the Transfermarkt football dataset.

Install Python 3.11+, VS Code, and the Jupyter extension

Both scikit-learn 1.8.0 and lifelines 0.30.3 require Python 3.11 or higher. Let's confirm you have the right version.

• Open your terminal by pressing Cmd+Space (macOS) or the Windows key (Windows), then type Terminal (macOS) or PowerShell (Windows) and press Enter.
• Check your Python version by running:

python3 --version

✔️ I see version 3.11 or higher

You are ready to continue.

ⓧ I see an older version

macOS ships with an older Python by default. Even if you installed a newer version, your terminal might resolve to the system Python.

• Download and install Python 3.11 from the official site: https://www.python.org/downloads/.
• After installation, close and reopen your terminal.
• Try python3.11 --version to verify the new version is accessible.

If your system still resolves to the old version, you can use the full path to the new binary or install a version manager like pyenv.

ⓧ Command not found

• Download and install Python 3.11 from https://www.python.org/downloads/.
• During installation on Windows, check the box that says Add Python to PATH.
• Close and reopen your terminal after installation.
• Run python3 --version again to confirm.

Next, open VS Code.

• Press Cmd+Space (macOS) or the Windows key (Windows) to open your search bar.
• Type Visual Studio Code and press Enter to open it.

VS Code not found?

Download and install VS Code from https://code.visualstudio.com/.

help me install VS Code on my operating system

• In VS Code, click the Extensions icon in the left sidebar (it looks like four small squares).
• Search for Jupyter and click Install on the extension published by Microsoft.

Why Jupyter in VS Code?

Jupyter notebooks let you run code in small cells and see output immediately. This is ideal for data exploration and watching MCMC sampling progress in real time.

Now let's check that conda is available.

• Back in your terminal, check for conda by running:

conda --version

✔️ I see a version number

Conda is ready. Continue to the next section.

ⓧ Command not found

You need to install Miniconda, a lightweight conda distribution.

• Download Miniconda from https://docs.conda.io/projects/conda/en/stable/ and follow the installer instructions for your operating system.
• After installation, close and reopen your terminal.
• Run conda --version again to verify.

help me install Miniconda on my operating system

Create the conda environment and install packages

You will define all your dependencies in an environment.yml file. Conda reads this file and installs everything into an isolated environment, including PyMC from the conda-forge channel.

• Move to your Desktop by running this command:

cd ~/Desktop

• Create your project folder and move into it by running these commands:

mkdir soccer-analytics
cd soccer-analytics

• Open this folder in VS Code by running:

code .

code command not found?

In VS Code, press Cmd+Shift+P (macOS) or Ctrl+Shift+P (Windows) to open the command palette. Type Shell Command: Install and select the option to install the code command in PATH. Then close and reopen your terminal and try again.

• In VS Code, create a new file called environment.yml in the soccer-analytics folder by clicking File then New File.
• Paste the following content into environment.yml:

name: soccer_analytics
channels:
  - conda-forge
dependencies:
  - python=3.11
  - pymc>=5.20
  - arviz>=0.13.0
  - pandas=2.2.3
  - scikit-learn=1.8.0
  - matplotlib>=3.0
  - seaborn>=0.12
  - jupyter
  - ipykernel
  - pip
  - pip:
    - lifelines==0.30.3

What does this file do?

The environment.yml file tells conda exactly what to install. It names the environment soccer_analytics, pulls packages from the conda-forge channel, and pins specific versions where compatibility matters.

PyMC is installed via conda-forge because it depends on PyTensor, which has compiled C extensions that conda resolves automatically. Lifelines is installed via pip because it is not available on conda-forge with the exact pinning we need.

• Save environment.yml.
• Create the conda environment by running this command in your terminal (inside the soccer-analytics folder):

conda env create -f environment.yml

How long does this take?

Conda needs to solve dependencies and download packages. This can take 5-10 minutes depending on your internet speed and machine. You will see progress bars as packages download.

• Once the environment is created, activate it by running:

conda activate soccer_analytics

Your terminal prompt should now show (soccer_analytics) at the beginning, confirming the environment is active.

• Create a requirements.txt file as a reference for non-conda users. In VS Code, create a new file called requirements.txt and paste:

pandas==2.2.3
scikit-learn==1.8.0
lifelines==0.30.3
matplotlib>=3.0
seaborn>=0.12
arviz>=0.13.0
jupyter
ipykernel

Why both environment.yml and requirements.txt?

The environment.yml is what you actually use to create your environment. The requirements.txt serves as documentation for anyone who prefers pip-only workflows. It lists the same dependencies (except PyMC, which needs conda for proper installation).

• Save requirements.txt.

Environment creation failed?

If conda reports a conflict, try running conda update conda first, then re-run the create command. Older conda versions sometimes struggle to resolve complex dependency trees.

help me debug my conda environment creation error for the soccer analytics project

Verify installations and download the dataset

With your environment ready, confirm that every package imported correctly. Then you will download the football dataset that powers the entire analysis pipeline.

• In VS Code, create a new file called verify_setup.py and paste:

# Verify all package installations
import pymc as pm
import lifelines
import arviz as az
import pandas as pd
import sklearn

print(f"PyMC version: {pm.__version__}")
print(f"lifelines version: {lifelines.__version__}")
print(f"ArviZ version: {az.__version__}")
print(f"pandas version: {pd.__version__}")
print(f"scikit-learn version: {sklearn.__version__}")
print("\nAll packages installed successfully!")

What does this code do?

Each import statement loads a library into memory. If any package is missing or broken, Python raises an ImportError immediately. Printing the version numbers confirms you have the right versions installed.

• Save verify_setup.py and run it in your terminal (make sure soccer_analytics is active):

python verify_setup.py

You should see version numbers printed for all five packages, followed by the success message.

Seeing an ImportError?

Make sure your terminal shows (soccer_analytics) in the prompt. If not, run conda activate soccer_analytics and try again.

If a specific package fails, reinstall it manually. For example: pip install lifelines==0.30.3 for lifelines, or conda install -c conda-forge pymc for PyMC.

help me fix an ImportError when verifying my soccer analytics environment

Now download the Transfermarkt football dataset.

• Create a data folder inside your soccer-analytics project by running:

mkdir data

• Open https://github.com/salimt/football-datasets in your browser.
• Navigate to the datalake/transfermarkt/raw/ directory.
• Download these three CSV files and save them into your data/ folder:

• player_profiles.csv (from the player_profiles/ subfolder).
• player_performances.csv (from the player_performances/ subfolder).
• player_injury_histories.csv (from the player_injury_histories/ subfolder).

How to download from GitHub

For each file, click into the subfolder, then click the filename to open it. Click the Download raw file button (the download icon) near the top right of the file view. Save it as the original filename into your data/ folder.

Note: player_performances.csv is tracked with Git LFS because it is large (1.8 million rows). If the download only gives you a small text file with LFS pointer content, you may need to clone the repository with Git LFS installed instead.

• Confirm your files are in place by running:

ls data/

You should see player_injury_histories.csv, player_performances.csv, and player_profiles.csv listed.

Only seeing small pointer files?

The performances file uses Git LFS. If you got a tiny text file instead of the full CSV, clone the entire repository with LFS support:

git lfs install
git clone https://github.com/salimt/football-datasets.git
cp football-datasets/datalake/transfermarkt/raw/player_performances/player_performances.csv data/

help me download the player_performances.csv file from the salimt/football-datasets repo using Git LFS

✔️ Awesome, I've got everything!

Make sure you have saved both environment.yml and requirements.txt.

ⓧ I'd like to double check the full code

Here is your environment.yml:

name: soccer_analytics
channels:
  - conda-forge
dependencies:
  - python=3.11
  - pymc>=5.20
  - arviz>=0.13.0
  - pandas=2.2.3
  - scikit-learn=1.8.0
  - matplotlib>=3.0
  - seaborn>=0.12
  - jupyter
  - ipykernel
  - pip
  - pip:
    - lifelines==0.30.3

Here is your requirements.txt:

pandas==2.2.3
scikit-learn==1.8.0
lifelines==0.30.3
matplotlib>=3.0
seaborn>=0.12
arviz>=0.13.0
jupyter
ipykernel

Your environment is ready and the dataset is in place. Next up, you will load these CSV files into pandas and explore the data structure before building your models.

Load, Explore, and Prepare the Data

Your CSV files are sitting in the data/ folder, but raw data is never analysis-ready. The Transfermarkt dataset spans 92,000+ players across multiple seasons, with dates stored as strings and performance spread across separate files.

In this step, you will load the three CSV files, merge them into one unified DataFrame, and engineer the features your models need. By the end, you will have a clean, single-season dataset with per-90 statistics, age calculations, and injury history ready for both Bayesian modeling and survival analysis.

In this step, get ready to:

• Load the three Transfermarkt CSV files and inspect their shapes and types.
• Merge profiles, performances, and injuries into one analysis DataFrame filtered to the 22/23 season.
• Engineer key features including age, position groups, per-90 rates, and previous injury counts.

Load the CSV files and inspect the data

Before any modeling, you need to read the raw data into memory and confirm the files loaded correctly. The dataset is split into three CSV files: player profiles, match performances, and injury histories.

• Open soccer_analytics.py in VS Code and add the imports at the top of the file:

import pandas as pd
import numpy as np
import pymc as pm
import arviz as az
from lifelines import KaplanMeierFitter, CoxPHFitter, WeibullAFTFitter
from lifelines.statistics import logrank_test
from lifelines.utils import concordance_index
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime

What are all these imports for?

• pandas and numpy handle data manipulation and numerical operations.
• pymc and arviz power the Bayesian hierarchical model you will build in a later step.
• lifelines provides Kaplan-Meier, Cox, and Weibull survival models.
• matplotlib and seaborn generate the visualizations.

• Below the imports, add the load_data() function that reads all three CSV files:

# =============================================================================
# Step 2: Load, Explore, and Prepare the Data
# =============================================================================

def load_data(profiles_path, performances_path, injuries_path):
    """Load the three Transfermarkt CSV files."""
    profiles = pd.read_csv(profiles_path)
    performances = pd.read_csv(performances_path)
    injuries = pd.read_csv(injuries_path)
    return profiles, performances, injuries

What does this function do?

load_data() takes three file paths and returns three separate pandas DataFrames. Keeping the loading logic in its own function makes the pipeline easy to reuse with different data paths.

• At the bottom of the file, add the main pipeline block that defines paths and calls load_data():

# =============================================================================
# Main Pipeline
# =============================================================================

if __name__ == "__main__":
    PROFILES_PATH = "data/player_profiles.csv"
    PERFORMANCES_PATH = "data/player_performances.csv"
    INJURIES_PATH = "data/player_injury_histories.csv"

    # Step 2: Load and prepare
    print("Loading data...")
    profiles, performances, injuries = load_data(PROFILES_PATH, PERFORMANCES_PATH, INJURIES_PATH)

What does this code do?

The if __name__ == "__main__": guard ensures this code only runs when you execute the file directly, not when it is imported as a module. The three path constants point to the CSV files you downloaded in Step 1.

• Save soccer_analytics.py.
• Run the script in your terminal to verify the data loads successfully:

python soccer_analytics.py

You should see Loading data... printed with no errors.

Seeing a FileNotFoundError?

• Confirm your CSV files are inside a data/ folder at the same level as soccer_analytics.py.
• Check the file names match exactly: player_profiles.csv, player_performances.csv, player_injury_histories.csv.
• Make sure you are running the script from the project folder (where data/ lives), not from inside data/ itself.

help me debug a FileNotFoundError when loading CSV files in pandas

Merge and aggregate the datasets

The three DataFrames contain different slices of information about the same players. You need to join them on player_id, filter to a single season, and aggregate match-level stats into player-level totals.

• Below load_data(), add the clean_and_merge() function:

def clean_and_merge(profiles, performances, injuries, target_season="22/23"):
    """Clean data and merge into a unified DataFrame."""
    injuries["from_date"] = pd.to_datetime(
        injuries["from_date"], format="%d.%m.%Y", errors="coerce"
    )
    injuries["end_date"] = pd.to_datetime(
        injuries["end_date"], format="%d.%m.%Y", errors="coerce"
    )

    perf_season = performances[performances["season"] == target_season].copy()
    inj_season = injuries[injuries["season"] == target_season].copy()

    perf_agg = perf_season.groupby("player_id").agg(
        total_minutes=("minutes_played", "sum"),
        total_goals=("goals", "sum"),
        total_assists=("assists", "sum"),
        total_yellow_cards=("yellow_cards", "sum"),
        total_appearances=("nb_on_pitch", "sum"),
        team_id=("team_id", "first"),
        team_name=("team_name", "first"),
    ).reset_index()

    inj_count = inj_season.groupby("player_id").agg(
        injury_count=("from_date", "count"),
        first_injury_date=("from_date", "min"),
        total_days_missed=(
            "days_missed",
            lambda x: pd.to_numeric(
                x.str.replace(" days", ""), errors="coerce"
            ).sum(),
        ),
    ).reset_index()

    merged = profiles[
        ["player_id", "name", "position", "date_of_birth", "height_in_cm"]
    ].copy()
    merged = merged.merge(perf_agg, on="player_id", how="inner")
    merged = merged.merge(inj_count, on="player_id", how="left")

    merged["injury_count"] = merged["injury_count"].fillna(0).astype(int)
    merged["total_days_missed"] = merged["total_days_missed"].fillna(0)

    merged["date_of_birth"] = pd.to_datetime(merged["date_of_birth"], errors="coerce")
    season_start = datetime(2022, 8, 1)
    merged["age"] = ((season_start - merged["date_of_birth"]).dt.days / 365.25).round(1)

    merged = merged[merged["total_minutes"] > 90].copy()

    position_map = {
        "Attack": "Attack",
        "Midfield": "Midfield",
        "Defender": "Defence",
        "Goalkeeper": "Goalkeeper",
    }
    merged["position_group"] = merged["position"].map(position_map)
    merged = merged.dropna(subset=["position_group", "age"])

    merged["goals_per_90"] = (merged["total_goals"] / merged["total_minutes"]) * 90
    merged["assists_per_90"] = (merged["total_assists"] / merged["total_minutes"]) * 90

    return merged

What does this function do?

This function performs the entire data preparation pipeline in one pass:

• Parses from_date and end_date from the DD.MM.YYYY string format into proper datetime objects.
• Filters both performances and injuries to the 22/23 season only.
• Aggregates match-level stats (minutes, goals, assists, yellow cards, appearances) into player-level totals using groupby.
• Merges profiles, aggregated performance, and injury counts on player_id.
• Calculates player age from date_of_birth relative to the season start (August 2022).
• Maps raw positions into four groups (Attack, Midfield, Defence, Goalkeeper) and computes goals and assists per 90 minutes.

The days_missed column in the raw data contains strings like 11 days, so the lambda strips the text and converts to a numeric sum.

Why filter to total_minutes > 90?

Players with fewer than 90 total minutes (one full match) produce unreliable per-90 statistics. A player who scored once in 5 minutes would show an extreme 18 goals per 90 rate. Filtering removes these noise-generating edge cases.

• Save soccer_analytics.py.

You won't see a change yet. The next chunk calls this function and prints the result.

Engineer previous injury counts and run the pipeline

A player's injury history from prior seasons is one of the strongest predictors of future injury. You need a separate function to count injuries before the target season, since clean_and_merge() focuses on the current season.

• Below clean_and_merge(), add the count_previous_injuries() function:

def count_previous_injuries(injuries_df, target_season="22/23"):
    """Count injuries in seasons prior to the target season."""
    seasons_order = ["18/19", "19/20", "20/21", "21/22"]
    prior = injuries_df[injuries_df["season"].isin(seasons_order)]
    prior_counts = prior.groupby("player_id").size().reset_index(name="previous_injuries")
    return prior_counts

What does this function do?

count_previous_injuries() filters the injury DataFrame to the four seasons before 22/23 and counts how many injuries each player had. This gives the survival model a recurrence risk indicator. Players who were injured more often in the past are more likely to be injured again.

• Update the main pipeline block to call clean_and_merge(), merge previous injuries, and print the final dataset shape. Replace the existing if __name__ block with:

# =============================================================================
# Main Pipeline
# =============================================================================

if __name__ == "__main__":
    PROFILES_PATH = "data/player_profiles.csv"
    PERFORMANCES_PATH = "data/player_performances.csv"
    INJURIES_PATH = "data/player_injury_histories.csv"

    # Step 2: Load and prepare
    print("Loading data...")
    profiles, performances, injuries = load_data(PROFILES_PATH, PERFORMANCES_PATH, INJURIES_PATH)

    print("Cleaning and merging...")
    df = clean_and_merge(profiles, performances, injuries, target_season="22/23")
    prior_injuries = count_previous_injuries(injuries, target_season="22/23")
    df = df.merge(prior_injuries, on="player_id", how="left")
    df["previous_injuries"] = df["previous_injuries"].fillna(0)
    print(f"Merged dataset: {df.shape[0]} players")

What does this pipeline code do?

• clean_and_merge() handles date parsing, season filtering, aggregation, merging, and feature engineering in one call.
• count_previous_injuries() adds a historical injury count from the four prior seasons.
• The left merge ensures players with no prior injuries get NaN, which is then filled with 0.
• The final print confirms how many players survived all the filtering and cleaning.

• Save soccer_analytics.py and run the full pipeline:

python soccer_analytics.py

You should see output like this:

Loading data...
Cleaning and merging...
Merged dataset: XXXX players

The exact player count depends on the data, but you should see a number in the thousands (the 92,000+ raw profiles are reduced after filtering to one season and removing players with fewer than 90 minutes).

Seeing a KeyError or empty DataFrame?

• If you see a KeyError on a column name, check that your CSV files match the expected schema. The performances file must have columns like minutes_played, goals, and season.
• If the merged dataset shows 0 players, confirm the season value 22/23 matches how seasons are stored in your data. Open the CSV and check a few rows.
• If you see a warning about mixed types, that is normal for large CSVs. The errors="coerce" arguments handle malformed values gracefully.

help me debug why my merged DataFrame is empty or has KeyErrors after clean_and_merge()

✔️ Awesome, I've got everything!

Great. Double check you have saved soccer_analytics.py and the script runs without errors.

ⓧ I'd like to double check the full code

Here is the complete soccer_analytics.py file at the end of this step:

import pandas as pd
import numpy as np
import pymc as pm
import arviz as az
from lifelines import KaplanMeierFitter, CoxPHFitter, WeibullAFTFitter
from lifelines.statistics import logrank_test
from lifelines.utils import concordance_index
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime

# =============================================================================
# Step 2: Load, Explore, and Prepare the Data
# =============================================================================

def load_data(profiles_path, performances_path, injuries_path):
    """Load the three Transfermarkt CSV files."""
    profiles = pd.read_csv(profiles_path)
    performances = pd.read_csv(performances_path)
    injuries = pd.read_csv(injuries_path)
    return profiles, performances, injuries


def clean_and_merge(profiles, performances, injuries, target_season="22/23"):
    """Clean data and merge into a unified DataFrame."""
    injuries["from_date"] = pd.to_datetime(
        injuries["from_date"], format="%d.%m.%Y", errors="coerce"
    )
    injuries["end_date"] = pd.to_datetime(
        injuries["end_date"], format="%d.%m.%Y", errors="coerce"
    )

    perf_season = performances[performances["season"] == target_season].copy()
    inj_season = injuries[injuries["season"] == target_season].copy()

    perf_agg = perf_season.groupby("player_id").agg(
        total_minutes=("minutes_played", "sum"),
        total_goals=("goals", "sum"),
        total_assists=("assists", "sum"),
        total_yellow_cards=("yellow_cards", "sum"),
        total_appearances=("nb_on_pitch", "sum"),
        team_id=("team_id", "first"),
        team_name=("team_name", "first"),
    ).reset_index()

    inj_count = inj_season.groupby("player_id").agg(
        injury_count=("from_date", "count"),
        first_injury_date=("from_date", "min"),
        total_days_missed=(
            "days_missed",
            lambda x: pd.to_numeric(
                x.str.replace(" days", ""), errors="coerce"
            ).sum(),
        ),
    ).reset_index()

    merged = profiles[
        ["player_id", "name", "position", "date_of_birth", "height_in_cm"]
    ].copy()
    merged = merged.merge(perf_agg, on="player_id", how="inner")
    merged = merged.merge(inj_count, on="player_id", how="left")

    merged["injury_count"] = merged["injury_count"].fillna(0).astype(int)
    merged["total_days_missed"] = merged["total_days_missed"].fillna(0)

    merged["date_of_birth"] = pd.to_datetime(merged["date_of_birth"], errors="coerce")
    season_start = datetime(2022, 8, 1)
    merged["age"] = ((season_start - merged["date_of_birth"]).dt.days / 365.25).round(1)

    merged = merged[merged["total_minutes"] > 90].copy()

    position_map = {
        "Attack": "Attack",
        "Midfield": "Midfield",
        "Defender": "Defence",
        "Goalkeeper": "Goalkeeper",
    }
    merged["position_group"] = merged["position"].map(position_map)
    merged = merged.dropna(subset=["position_group", "age"])

    merged["goals_per_90"] = (merged["total_goals"] / merged["total_minutes"]) * 90
    merged["assists_per_90"] = (merged["total_assists"] / merged["total_minutes"]) * 90

    return merged


def count_previous_injuries(injuries_df, target_season="22/23"):
    """Count injuries in seasons prior to the target season."""
    seasons_order = ["18/19", "19/20", "20/21", "21/22"]
    prior = injuries_df[injuries_df["season"].isin(seasons_order)]
    prior_counts = prior.groupby("player_id").size().reset_index(name="previous_injuries")
    return prior_counts


# =============================================================================
# Main Pipeline
# =============================================================================

if __name__ == "__main__":
    PROFILES_PATH = "data/player_profiles.csv"
    PERFORMANCES_PATH = "data/player_performances.csv"
    INJURIES_PATH = "data/player_injury_histories.csv"

    # Step 2: Load and prepare
    print("Loading data...")
    profiles, performances, injuries = load_data(PROFILES_PATH, PERFORMANCES_PATH, INJURIES_PATH)

    print("Cleaning and merging...")
    df = clean_and_merge(profiles, performances, injuries, target_season="22/23")
    prior_injuries = count_previous_injuries(injuries, target_season="22/23")
    df = df.merge(prior_injuries, on="player_id", how="left")
    df["previous_injuries"] = df["previous_injuries"].fillna(0)
    print(f"Merged dataset: {df.shape[0]} players")

Your data is clean, merged, and feature-enriched. Next up, you will use this prepared DataFrame to build a Bayesian hierarchical model that estimates player performance with team-level random effects.

Build the Bayesian Hierarchical Performance Model

You have a clean dataset of player performance statistics. Now you need a model that can estimate how good each player truly is, while accounting for the fact that players on stronger teams will look better simply because of their environment.

A Bayesian hierarchical model solves this by learning team-level effects and player-level effects simultaneously. Teams with fewer observations get "pulled" toward the league average through partial pooling, giving more reliable estimates than raw scoring rates.

In this step, get ready to:

• Define a hierarchical model with team-level random effects and player covariates in PyMC.
• Sample the posterior distribution using MCMC (NUTS sampler).
• Visualize team and position effects with ArviZ and extract player-level performance posteriors.

Define the hierarchical model structure

The model needs to capture how team membership and player characteristics jointly influence scoring rates. You will define hyperpriors that control how much variation exists between teams, then let each team's attack strength be drawn from that shared distribution.

• In soccer_analytics.py, find the comment # ============================================================================= that appears right above the # Main Pipeline section.
• Add a new section heading and the build_hierarchical_model() function above the Main Pipeline section by pasting this code:

# =============================================================================
# Step 3: Bayesian Hierarchical Performance Model
# =============================================================================

def build_hierarchical_model(df):
    """
    Build a Bayesian hierarchical model for player performance.
    Models goals_per_90 with team-level random effects and player covariates.
    """
    # Prepare team indices
    teams = df["team_name"].unique()
    team_to_idx = {team: i for i, team in enumerate(teams)}
    df["team_idx"] = df["team_name"].map(team_to_idx)

    # Prepare position indices
    positions = df["position_group"].unique()
    pos_to_idx = {pos: i for i, pos in enumerate(positions)}
    df["pos_idx"] = df["position_group"].map(pos_to_idx)

    # Standardize continuous covariates
    age_std = (df["age"] - df["age"].mean()) / df["age"].std()
    minutes_std = (df["total_minutes"] - df["total_minutes"].mean()) / df["total_minutes"].std()

    # Observed data
    observed_goals_per_90 = df["goals_per_90"].values

    coords = {
        "team": teams,
        "position": positions,
        "obs": np.arange(len(df)),
    }

    with pm.Model(coords=coords) as model:
        # Hyperprior: how much do teams vary in attack strength?
        sd_team = pm.HalfNormal("sd_team", sigma=0.5)

        # Team-level random effects (attack strength)
        team_effect = pm.Normal("team_effect", mu=0, sigma=sd_team, dims="team")

        # Position fixed effects
        pos_effect = pm.Normal("pos_effect", mu=0, sigma=1, dims="position")

        # Player-level covariates
        beta_age = pm.Normal("beta_age", mu=0, sigma=0.5)
        beta_minutes = pm.Normal("beta_minutes", mu=0, sigma=0.5)

        # Global intercept (league-average scoring rate)
        intercept = pm.Normal("intercept", mu=0.3, sigma=0.5)

        # Expected goals per 90 for each player
        mu = (
            intercept
            + team_effect[df["team_idx"].values]
            + pos_effect[df["pos_idx"].values]
            + beta_age * age_std.values
            + beta_minutes * minutes_std.values
        )

        # Observation noise
        sigma = pm.HalfNormal("sigma", sigma=0.5)
        goals_obs = pm.Normal(
            "goals_obs",
            mu=mu,
            sigma=sigma,
            observed=observed_goals_per_90,
            dims="obs",
        )

    return model, df, coords

What does this model do?

This function defines a probabilistic model where each player's goals per 90 minutes is explained by multiple layers of effects.

• sd_team is a hyperprior controlling how different teams are from each other. The model learns this from data.
• team_effect gives each team its own attack strength, drawn from a shared Normal distribution. This is partial pooling: small-sample teams shrink toward the mean.
• pos_effect, beta_age, and beta_minutes capture position differences and player-level covariates.
• The dims argument tells PyMC the shape of each parameter, linking it to the coords dictionary for labeled inference results.

Why partial pooling instead of raw averages?

If you compute goals per 90 directly, a player with 2 goals in 100 minutes looks better than one with 10 goals in 2000 minutes. Partial pooling borrows information across the league, shrinking extreme estimates from small samples toward the group mean. The result is more stable, trustworthy performance ratings.

• Save soccer_analytics.py.

You won't see a visible change yet. The next section samples this model and produces output.

Sample the posterior and analyze results

With the model structure defined, you now need to draw samples from the posterior distribution using MCMC. PyMC automatically selects the NUTS sampler, which efficiently explores high-dimensional parameter spaces.

• Add the sample_hierarchical_model() and analyze_posteriors() functions below build_hierarchical_model() by pasting this code:

def sample_hierarchical_model(model, draws=1000, tune=1500, cores=2):
    """Sample the posterior using NUTS."""
    with model:
        trace = pm.sample(draws=draws, tune=tune, cores=cores, random_seed=42)
    return trace


def analyze_posteriors(trace, df, coords):
    """Extract and visualize posterior summaries."""
    # Print team effect rankings
    print("=== Team Effect Posterior Summary ===")
    team_summary = az.summary(trace, var_names=["team_effect"])
    print(team_summary.head(10))

    # Plot team attack strength with credible intervals
    fig, ax = plt.subplots(figsize=(12, 6))
    az.plot_forest(trace, var_names=["team_effect"], ax=ax, combined=True)
    ax.set_title("Team Attack Strength (Posterior 94% HDI)")
    plt.tight_layout()
    plt.savefig("team_effects.png", dpi=150)
    plt.show()

    # Plot position effects
    fig, ax = plt.subplots(figsize=(8, 4))
    az.plot_forest(trace, var_names=["pos_effect"], ax=ax, combined=True)
    ax.set_title("Position Effects on Goals per 90")
    plt.tight_layout()
    plt.savefig("position_effects.png", dpi=150)
    plt.show()

    # Convergence diagnostics for key parameters
    print("\n=== Convergence Diagnostics ===")
    print(az.summary(trace, var_names=["intercept", "beta_age", "beta_minutes", "sd_team", "sigma"]))

    # Posterior predictive check
    with trace.model:
        pm.sample_posterior_predictive(trace, extend_inferencedata=True)

    fig, ax = plt.subplots(figsize=(10, 5))
    az.plot_ppc(trace, ax=ax, num_pp_samples=100)
    ax.set_title("Posterior Predictive Check")
    plt.tight_layout()
    plt.savefig("ppc_check.png", dpi=150)
    plt.show()

    return trace

What does this code do?

• sample_hierarchical_model() runs the NUTS sampler with 1000 posterior draws and 1500 tuning steps per chain. Tuning lets the sampler adapt its step size before collecting real samples.
• analyze_posteriors() prints summaries and creates three plots: a forest plot of team attack strengths ranked by their posterior means, a forest plot of position effects, and a posterior predictive check that overlays simulated data on your observed data.
• az.summary() reports the mean, standard deviation, and 94% HDI (highest density interval) for each parameter. R-hat values near 1.0 indicate the sampler converged.

Sampling takes a few minutes

MCMC sampling with NUTS is computationally intensive. Expect 2 to 5 minutes depending on your hardware. The progress bar shows sampling speed and estimated time remaining.

• Save soccer_analytics.py.

You won't see output yet because you haven't called these functions from the main pipeline. The next section connects everything.

Extract player posteriors and run the pipeline

The final piece combines each player's team and position effects into a single performance score. This bayesian_performance column will feed into the survival analysis in later steps.

• Add the extract_player_posteriors() function below analyze_posteriors():

def extract_player_posteriors(trace, df):
    """Extract individual player performance posterior means for use in survival model."""
    # Get posterior means for team and position effects
    team_effects_mean = trace.posterior["team_effect"].mean(dim=["chain", "draw"]).values
    pos_effects_mean = trace.posterior["pos_effect"].mean(dim=["chain", "draw"]).values

    # Combine team + position effects as the player's Bayesian performance estimate
    df["bayesian_performance"] = (
        team_effects_mean[df["team_idx"].values]
        + pos_effects_mean[df["pos_idx"].values]
    )

    return df

What does this code do?

This function averages the posterior samples for each team and position effect, then adds them together to create a single performance score per player. A player on a strong attacking team in a forward position gets a higher bayesian_performance value than a goalkeeper on a weaker team.

This is partial pooling in action. Each team's estimate is informed by its own data AND the league-wide distribution, producing more stable ratings than raw averages.

• Now update the main pipeline to call these functions. In the if __name__ == "__main__": block, find the line:

    print(f"Merged dataset: {df.shape[0]} players")

• Add the following lines directly below it:

    # Step 3: Bayesian hierarchical model
    print("\nBuilding Bayesian hierarchical model...")
    model, df, coords = build_hierarchical_model(df)
    print("Sampling posterior (this may take a few minutes)...")
    trace = sample_hierarchical_model(model, draws=1000, tune=1500, cores=2)
    trace = analyze_posteriors(trace, df, coords)
    df = extract_player_posteriors(trace, df)

• Save soccer_analytics.py.
• Run the script in your terminal:

python soccer_analytics.py

You should see progress bars for MCMC sampling, followed by printed posterior summaries showing team effects and convergence diagnostics. Three plot windows will appear (or save to team_effects.png, position_effects.png, and ppc_check.png).

Sampling seems stuck or throws an error?

• If you see divergence warnings, that is normal for a first run. The model may report a small number of divergences that do not affect results.
• If you get a memory error, reduce cores=2 to cores=1 in the sample_hierarchical_model() call.
• If PyMC cannot be imported, confirm your conda environment is active by running conda activate soccer_analytics before executing the script.

help me debug my PyMC sampling issues

✔️ Awesome, I've got everything!

Great. Double check you have saved soccer_analytics.py and that the posterior summary printed without errors.

ⓧ I'd like to double check the full code

Here is the full soccer_analytics.py at this point:

import pandas as pd
import numpy as np
import pymc as pm
import arviz as az
from lifelines import KaplanMeierFitter, CoxPHFitter, WeibullAFTFitter
from lifelines.statistics import logrank_test
from lifelines.utils import concordance_index
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime

# =============================================================================
# Step 2: Load, Explore, and Prepare the Data
# =============================================================================

def load_data(profiles_path, performances_path, injuries_path):
    """Load the three Transfermarkt CSV files."""
    profiles = pd.read_csv(profiles_path)
    performances = pd.read_csv(performances_path)
    injuries = pd.read_csv(injuries_path)
    return profiles, performances, injuries


def clean_and_merge(profiles, performances, injuries, target_season="22/23"):
    """Clean data and merge into a unified DataFrame."""
    injuries["from_date"] = pd.to_datetime(
        injuries["from_date"], format="%d.%m.%Y", errors="coerce"
    )
    injuries["end_date"] = pd.to_datetime(
        injuries["end_date"], format="%d.%m.%Y", errors="coerce"
    )

    perf_season = performances[performances["season"] == target_season].copy()
    inj_season = injuries[injuries["season"] == target_season].copy()

    perf_agg = perf_season.groupby("player_id").agg(
        total_minutes=("minutes_played", "sum"),
        total_goals=("goals", "sum"),
        total_assists=("assists", "sum"),
        total_yellow_cards=("yellow_cards", "sum"),
        total_appearances=("nb_on_pitch", "sum"),
        team_id=("team_id", "first"),
        team_name=("team_name", "first"),
    ).reset_index()

    inj_count = inj_season.groupby("player_id").agg(
        injury_count=("from_date", "count"),
        first_injury_date=("from_date", "min"),
        total_days_missed=(
            "days_missed",
            lambda x: pd.to_numeric(
                x.str.replace(" days", ""), errors="coerce"
            ).sum(),
        ),
    ).reset_index()

    merged = profiles[
        ["player_id", "name", "position", "date_of_birth", "height_in_cm"]
    ].copy()
    merged = merged.merge(perf_agg, on="player_id", how="inner")
    merged = merged.merge(inj_count, on="player_id", how="left")

    merged["injury_count"] = merged["injury_count"].fillna(0).astype(int)
    merged["total_days_missed"] = merged["total_days_missed"].fillna(0)

    merged["date_of_birth"] = pd.to_datetime(merged["date_of_birth"], errors="coerce")
    season_start = datetime(2022, 8, 1)
    merged["age"] = ((season_start - merged["date_of_birth"]).dt.days / 365.25).round(1)

    merged = merged[merged["total_minutes"] > 90].copy()

    position_map = {
        "Attack": "Attack",
        "Midfield": "Midfield",
        "Defender": "Defence",
        "Goalkeeper": "Goalkeeper",
    }
    merged["position_group"] = merged["position"].map(position_map)
    merged = merged.dropna(subset=["position_group", "age"])

    merged["goals_per_90"] = (merged["total_goals"] / merged["total_minutes"]) * 90
    merged["assists_per_90"] = (merged["total_assists"] / merged["total_minutes"]) * 90

    return merged


def count_previous_injuries(injuries_df, target_season="22/23"):
    """Count injuries in seasons prior to the target season."""
    seasons_order = ["18/19", "19/20", "20/21", "21/22"]
    prior = injuries_df[injuries_df["season"].isin(seasons_order)]
    prior_counts = prior.groupby("player_id").size().reset_index(name="previous_injuries")
    return prior_counts


# =============================================================================
# Step 3: Bayesian Hierarchical Performance Model
# =============================================================================

def build_hierarchical_model(df):
    """
    Build a Bayesian hierarchical model for player performance.
    Models goals_per_90 with team-level random effects and player covariates.
    """
    teams = df["team_name"].unique()
    team_to_idx = {team: i for i, team in enumerate(teams)}
    df["team_idx"] = df["team_name"].map(team_to_idx)

    positions = df["position_group"].unique()
    pos_to_idx = {pos: i for i, pos in enumerate(positions)}
    df["pos_idx"] = df["position_group"].map(pos_to_idx)

    age_std = (df["age"] - df["age"].mean()) / df["age"].std()
    minutes_std = (df["total_minutes"] - df["total_minutes"].mean()) / df["total_minutes"].std()

    observed_goals_per_90 = df["goals_per_90"].values

    coords = {
        "team": teams,
        "position": positions,
        "obs": np.arange(len(df)),
    }

    with pm.Model(coords=coords) as model:
        sd_team = pm.HalfNormal("sd_team", sigma=0.5)
        team_effect = pm.Normal("team_effect", mu=0, sigma=sd_team, dims="team")
        pos_effect = pm.Normal("pos_effect", mu=0, sigma=1, dims="position")
        beta_age = pm.Normal("beta_age", mu=0, sigma=0.5)
        beta_minutes = pm.Normal("beta_minutes", mu=0, sigma=0.5)
        intercept = pm.Normal("intercept", mu=0.3, sigma=0.5)
        mu = (
            intercept
            + team_effect[df["team_idx"].values]
            + pos_effect[df["pos_idx"].values]
            + beta_age * age_std.values
            + beta_minutes * minutes_std.values
        )
        sigma = pm.HalfNormal("sigma", sigma=0.5)
        goals_obs = pm.Normal(
            "goals_obs",
            mu=mu,
            sigma=sigma,
            observed=observed_goals_per_90,
            dims="obs",
        )

    return model, df, coords


def sample_hierarchical_model(model, draws=1000, tune=1500, cores=2):
    """Sample the posterior using NUTS."""
    with model:
        trace = pm.sample(draws=draws, tune=tune, cores=cores, random_seed=42)
    return trace


def analyze_posteriors(trace, df, coords):
    """Extract and visualize posterior summaries."""
    print("=== Team Effect Posterior Summary ===")
    team_summary = az.summary(trace, var_names=["team_effect"])
    print(team_summary.head(10))

    fig, ax = plt.subplots(figsize=(12, 6))
    az.plot_forest(trace, var_names=["team_effect"], ax=ax, combined=True)
    ax.set_title("Team Attack Strength (Posterior 94% HDI)")
    plt.tight_layout()
    plt.savefig("team_effects.png", dpi=150)
    plt.show()

    fig, ax = plt.subplots(figsize=(8, 4))
    az.plot_forest(trace, var_names=["pos_effect"], ax=ax, combined=True)
    ax.set_title("Position Effects on Goals per 90")
    plt.tight_layout()
    plt.savefig("position_effects.png", dpi=150)
    plt.show()

    print("\n=== Convergence Diagnostics ===")
    print(az.summary(trace, var_names=["intercept", "beta_age", "beta_minutes", "sd_team", "sigma"]))

    with trace.model:
        pm.sample_posterior_predictive(trace, extend_inferencedata=True)

    fig, ax = plt.subplots(figsize=(10, 5))
    az.plot_ppc(trace, ax=ax, num_pp_samples=100)
    ax.set_title("Posterior Predictive Check")
    plt.tight_layout()
    plt.savefig("ppc_check.png", dpi=150)
    plt.show()

    return trace


def extract_player_posteriors(trace, df):
    """Extract individual player performance posterior means for use in survival model."""
    team_effects_mean = trace.posterior["team_effect"].mean(dim=["chain", "draw"]).values
    pos_effects_mean = trace.posterior["pos_effect"].mean(dim=["chain", "draw"]).values
    df["bayesian_performance"] = (
        team_effects_mean[df["team_idx"].values]
        + pos_effects_mean[df["pos_idx"].values]
    )
    return df


# =============================================================================
# Main Pipeline
# =============================================================================

if __name__ == "__main__":
    PROFILES_PATH = "data/player_profiles.csv"
    PERFORMANCES_PATH = "data/player_performances.csv"
    INJURIES_PATH = "data/player_injury_histories.csv"

    print("Loading data...")
    profiles, performances, injuries = load_data(PROFILES_PATH, PERFORMANCES_PATH, INJURIES_PATH)

    print("Cleaning and merging...")
    df = clean_and_merge(profiles, performances, injuries, target_season="22/23")
    prior_injuries = count_previous_injuries(injuries, target_season="22/23")
    df = df.merge(prior_injuries, on="player_id", how="left")
    df["previous_injuries"] = df["previous_injuries"].fillna(0)
    print(f"Merged dataset: {df.shape[0]} players")

    print("\nBuilding Bayesian hierarchical model...")
    model, df, coords = build_hierarchical_model(df)
    print("Sampling posterior (this may take a few minutes)...")
    trace = sample_hierarchical_model(model, draws=1000, tune=1500, cores=2)
    trace = analyze_posteriors(trace, df, coords)
    df = extract_player_posteriors(trace, df)

Your Bayesian model is trained. You now have team-level and position-level performance estimates stored in the trace, and each player has a bayesian_performance score. Next, you will transform the injury data into survival analysis format so you can predict when injuries happen.

Prepare Survival Analysis Data

Your Bayesian model now produces posterior performance estimates for every player. The next question is not how well a player performs, but how long they stay healthy enough to keep performing.

Survival analysis models time-to-event data. In this step, you will transform your merged DataFrame into the format that survival models expect: a "duration" (how many days until injury) and an "event" indicator (did the injury actually happen, or did the season end first).

In this step, get ready to:

• Transform injury records into survival format with duration and event columns.
• Identify available survival covariates from your earlier data preparation.
• Validate the survival dataset for structural integrity.

Transform injury data into survival format

Survival models need two pieces of information per observation. The first is duration: how many days into the season until the event occurred. The second is an event indicator: 1 if the player was injured, 0 if they completed the season injury-free (right-censored).

What is right-censoring?

Right-censoring occurs when the observation period ends before the event of interest happens. A player who finishes the whole season without injury is censored. We know their injury-free time is at least 303 days, but we do not know when (or if) they would eventually get injured.

Survival models handle this naturally. Dropping censored players would bias results toward injury-prone players only.

• In soccer_analytics.py, find the comment line that says # Step 4: Prepare Survival Analysis Data.
• Add the prepare_survival_data function below the section header by copying this code:

def prepare_survival_data(df):
    """Transform data into survival analysis format."""
    season_length_days = 303

    df["has_injury"] = (df["injury_count"] > 0).astype(int)
    df["minutes_fraction"] = df["total_minutes"] / df["total_minutes"].max()
    df["duration"] = np.where(
        df["has_injury"] == 1,
        np.clip(df["minutes_fraction"] * season_length_days * 0.6, 1, season_length_days),
        season_length_days,
    )
    df["event"] = df["has_injury"]

    return df

What does this code do?

• season_length_days = 303 defines the season length (August to June, roughly 303 days).
• has_injury is a binary flag: 1 if the player had any injury that season, 0 otherwise.
• minutes_fraction normalizes each player's total minutes against the most-played player, giving a proxy for how deep into the season they were active.
• duration estimates days until first injury for injured players (clipped to [1, 303]) or assigns the full 303 days for censored (uninjured) players.
• event mirrors has_injury. A value of 1 means the injury event was observed. A value of 0 means the player was right-censored.

• Now scroll down to the if __name__ == "__main__": section at the bottom of the file.
• Find the last line that says df = extract_player_posteriors(trace, df).
• Add the pipeline call for survival data preparation below it:

    print("\nPreparing survival analysis data...")
    df = prepare_survival_data(df)

• Save soccer_analytics.py.

Seeing a NameError for np?

Make sure your import section at the top of the file includes import numpy as np. This was added in Step 2 but may have been accidentally removed.

help me debug a NameError in my prepare_survival_data function

Identify available survival covariates

Survival models use covariates to explain why some players get injured sooner than others. Thanks to the data preparation in Step 2 and the Bayesian model in Step 3, your DataFrame already contains the covariates you need.

What covariates are available for survival modeling?

• total_minutes: cumulative minutes played (workload proxy).
• previous_injuries: count of injuries in prior seasons (recurrence risk indicator).
• age: player age at season start.
• position_group: categorical position (Attack, Midfield, Defence, Goalkeeper).
• bayesian_performance: the posterior mean from your hierarchical model, linking performance to injury risk.
• total_yellow_cards: an aggression proxy that may correlate with higher injury exposure.

You will standardize and feed these covariates into the Cox and Weibull models in Step 6. For now, confirm they exist in your DataFrame.

• Open a new Jupyter notebook cell and run this to inspect the available columns:

# Confirm covariates are present
covariate_cols = ["age", "total_minutes", "previous_injuries",
                  "bayesian_performance", "total_yellow_cards", "position_group"]
print("Available covariates:")
print(df[covariate_cols].describe())

What should I see?

You should see summary statistics (count, mean, std, min, max) for the numeric covariates. All columns should have non-zero counts matching your player count. The position_group column will show as object type with top/freq stats instead of numeric summary.

Validate the survival dataset

Before feeding data into survival models, you need to confirm the structure is valid. Negative durations or extreme outliers would produce meaningless results.

• Run these validation checks in a new notebook cell:

# Validate survival data structure
print("=== Survival Data Validation ===")
print(f"Total players: {len(df)}")
print(f"Injured (event=1): {df['event'].sum()} ({df['event'].mean()*100:.1f}%)")
print(f"Censored (event=0): {(df['event']==0).sum()} ({(1-df['event'].mean())*100:.1f}%)")
print(f"\nDuration range: [{df['duration'].min():.1f}, {df['duration'].max():.1f}]")
print(f"Negative durations: {(df['duration'] < 0).sum()}")
print(f"\nCovariate outlier check (values beyond 3 std):")
for col in ["age", "total_minutes", "previous_injuries", "total_yellow_cards"]:
    mean = df[col].mean()
    std = df[col].std()
    outliers = ((df[col] - mean).abs() > 3 * std).sum()
    print(f"  {col}: {outliers} outliers")

What should I see?

You should see zero negative durations. The censoring proportion should be somewhere between 30-70%. If nearly all players are censored (>90%), the survival models will have very few events to learn from.

Duration values should range from 1 to 303. A few outliers in covariates like total_minutes are expected and acceptable. The standardization in Step 6 will handle them.

• Confirm that the "Negative durations" output shows 0.
• Confirm that the duration range starts at 1.0 and ends at 303.0.

Seeing negative durations or all zeros?

If you see negative durations, check that the np.clip call uses 1 as the lower bound. If all durations are zero or NaN, make sure total_minutes contains valid numeric data and that your earlier filter kept only players with more than 90 minutes.

help me debug my survival data transformation

✔️ Awesome, I've got everything!

Great. Double check you have saved soccer_analytics.py with the new function and pipeline call.

ⓧ I'd like to double check the full code

Here is the complete soccer_analytics.py file up to this point:

import pandas as pd
import numpy as np
import pymc as pm
import arviz as az
from lifelines import KaplanMeierFitter, CoxPHFitter, WeibullAFTFitter
from lifelines.statistics import logrank_test
from lifelines.utils import concordance_index
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime

# =============================================================================
# Step 2: Load, Explore, and Prepare the Data
# =============================================================================

def load_data(profiles_path, performances_path, injuries_path):
    """Load the three Transfermarkt CSV files."""
    profiles = pd.read_csv(profiles_path)
    performances = pd.read_csv(performances_path)
    injuries = pd.read_csv(injuries_path)
    return profiles, performances, injuries


def clean_and_merge(profiles, performances, injuries, target_season="22/23"):
    """Clean data and merge into a unified DataFrame."""
    injuries["from_date"] = pd.to_datetime(
        injuries["from_date"], format="%d.%m.%Y", errors="coerce"
    )
    injuries["end_date"] = pd.to_datetime(
        injuries["end_date"], format="%d.%m.%Y", errors="coerce"
    )

    perf_season = performances[performances["season"] == target_season].copy()
    inj_season = injuries[injuries["season"] == target_season].copy()

    perf_agg = perf_season.groupby("player_id").agg(
        total_minutes=("minutes_played", "sum"),
        total_goals=("goals", "sum"),
        total_assists=("assists", "sum"),
        total_yellow_cards=("yellow_cards", "sum"),
        total_appearances=("nb_on_pitch", "sum"),
        team_id=("team_id", "first"),
        team_name=("team_name", "first"),
    ).reset_index()

    inj_count = inj_season.groupby("player_id").agg(
        injury_count=("from_date", "count"),
        first_injury_date=("from_date", "min"),
        total_days_missed=(
            "days_missed",
            lambda x: pd.to_numeric(
                x.str.replace(" days", ""), errors="coerce"
            ).sum(),
        ),
    ).reset_index()

    merged = profiles[
        ["player_id", "name", "position", "date_of_birth", "height_in_cm"]
    ].copy()
    merged = merged.merge(perf_agg, on="player_id", how="inner")
    merged = merged.merge(inj_count, on="player_id", how="left")

    merged["injury_count"] = merged["injury_count"].fillna(0).astype(int)
    merged["total_days_missed"] = merged["total_days_missed"].fillna(0)

    merged["date_of_birth"] = pd.to_datetime(merged["date_of_birth"], errors="coerce")
    season_start = datetime(2022, 8, 1)
    merged["age"] = ((season_start - merged["date_of_birth"]).dt.days / 365.25).round(1)

    merged = merged[merged["total_minutes"] > 90].copy()

    position_map = {
        "Attack": "Attack",
        "Midfield": "Midfield",
        "Defender": "Defence",
        "Goalkeeper": "Goalkeeper",
    }
    merged["position_group"] = merged["position"].map(position_map)
    merged = merged.dropna(subset=["position_group", "age"])

    merged["goals_per_90"] = (merged["total_goals"] / merged["total_minutes"]) * 90
    merged["assists_per_90"] = (merged["total_assists"] / merged["total_minutes"]) * 90

    return merged


def count_previous_injuries(injuries_df, target_season="22/23"):
    """Count injuries in seasons prior to the target season."""
    seasons_order = ["18/19", "19/20", "20/21", "21/22"]
    prior = injuries_df[injuries_df["season"].isin(seasons_order)]
    prior_counts = prior.groupby("player_id").size().reset_index(name="previous_injuries")
    return prior_counts


# =============================================================================
# Step 3: Bayesian Hierarchical Performance Model
# =============================================================================

def build_hierarchical_model(df):
    teams = df["team_name"].unique()
    team_to_idx = {team: i for i, team in enumerate(teams)}
    df["team_idx"] = df["team_name"].map(team_to_idx)

    positions = df["position_group"].unique()
    pos_to_idx = {pos: i for i, pos in enumerate(positions)}
    df["pos_idx"] = df["position_group"].map(pos_to_idx)

    age_std = (df["age"] - df["age"].mean()) / df["age"].std()
    minutes_std = (df["total_minutes"] - df["total_minutes"].mean()) / df["total_minutes"].std()

    observed_goals_per_90 = df["goals_per_90"].values

    coords = {
        "team": teams,
        "position": positions,
        "obs": np.arange(len(df)),
    }

    with pm.Model(coords=coords) as model:
        sd_team = pm.HalfNormal("sd_team", sigma=0.5)
        team_effect = pm.Normal("team_effect", mu=0, sigma=sd_team, dims="team")
        pos_effect = pm.Normal("pos_effect", mu=0, sigma=1, dims="position")
        beta_age = pm.Normal("beta_age", mu=0, sigma=0.5)
        beta_minutes = pm.Normal("beta_minutes", mu=0, sigma=0.5)
        intercept = pm.Normal("intercept", mu=0.3, sigma=0.5)
        mu = (
            intercept
            + team_effect[df["team_idx"].values]
            + pos_effect[df["pos_idx"].values]
            + beta_age * age_std.values
            + beta_minutes * minutes_std.values
        )
        sigma = pm.HalfNormal("sigma", sigma=0.5)
        goals_obs = pm.Normal(
            "goals_obs",
            mu=mu,
            sigma=sigma,
            observed=observed_goals_per_90,
            dims="obs",
        )

    return model, df, coords


def sample_hierarchical_model(model, draws=1000, tune=1500, cores=2):
    with model:
        trace = pm.sample(draws=draws, tune=tune, cores=cores, random_seed=42)
    return trace


def analyze_posteriors(trace, df, coords):
    print("=== Team Effect Posterior Summary ===")
    team_summary = az.summary(trace, var_names=["team_effect"])
    print(team_summary.head(10))

    fig, ax = plt.subplots(figsize=(12, 6))
    az.plot_forest(trace, var_names=["team_effect"], ax=ax, combined=True)
    ax.set_title("Team Attack Strength (Posterior 94% HDI)")
    plt.tight_layout()
    plt.savefig("team_effects.png", dpi=150)
    plt.show()

    fig, ax = plt.subplots(figsize=(8, 4))
    az.plot_forest(trace, var_names=["pos_effect"], ax=ax, combined=True)
    ax.set_title("Position Effects on Goals per 90")
    plt.tight_layout()
    plt.savefig("position_effects.png", dpi=150)
    plt.show()

    print("\n=== Convergence Diagnostics ===")
    print(az.summary(trace, var_names=["intercept", "beta_age", "beta_minutes", "sd_team", "sigma"]))

    with trace.model:
        pm.sample_posterior_predictive(trace, extend_inferencedata=True)

    fig, ax = plt.subplots(figsize=(10, 5))
    az.plot_ppc(trace, ax=ax, num_pp_samples=100)
    ax.set_title("Posterior Predictive Check")
    plt.tight_layout()
    plt.savefig("ppc_check.png", dpi=150)
    plt.show()

    return trace


def extract_player_posteriors(trace, df):
    team_effects_mean = trace.posterior["team_effect"].mean(dim=["chain", "draw"]).values
    pos_effects_mean = trace.posterior["pos_effect"].mean(dim=["chain", "draw"]).values
    df["bayesian_performance"] = (
        team_effects_mean[df["team_idx"].values]
        + pos_effects_mean[df["pos_idx"].values]
    )
    return df


# =============================================================================
# Step 4: Prepare Survival Analysis Data
# =============================================================================

def prepare_survival_data(df):
    """Transform data into survival analysis format."""
    season_length_days = 303

    df["has_injury"] = (df["injury_count"] > 0).astype(int)
    df["minutes_fraction"] = df["total_minutes"] / df["total_minutes"].max()
    df["duration"] = np.where(
        df["has_injury"] == 1,
        np.clip(df["minutes_fraction"] * season_length_days * 0.6, 1, season_length_days),
        season_length_days,
    )
    df["event"] = df["has_injury"]

    return df


# =============================================================================
# Main Pipeline
# =============================================================================

if __name__ == "__main__":
    PROFILES_PATH = "data/player_profiles.csv"
    PERFORMANCES_PATH = "data/player_performances.csv"
    INJURIES_PATH = "data/player_injury_histories.csv"

    print("Loading data...")
    profiles, performances, injuries = load_data(PROFILES_PATH, PERFORMANCES_PATH, INJURIES_PATH)

    print("Cleaning and merging...")
    df = clean_and_merge(profiles, performances, injuries, target_season="22/23")
    prior_injuries = count_previous_injuries(injuries, target_season="22/23")
    df = df.merge(prior_injuries, on="player_id", how="left")
    df["previous_injuries"] = df["previous_injuries"].fillna(0)
    print(f"Merged dataset: {df.shape[0]} players")

    print("\nBuilding Bayesian hierarchical model...")
    model, df, coords = build_hierarchical_model(df)
    print("Sampling posterior (this may take a few minutes)...")
    trace = sample_hierarchical_model(model, draws=1000, tune=1500, cores=2)
    trace = analyze_posteriors(trace, df, coords)
    df = extract_player_posteriors(trace, df)

    print("\nPreparing survival analysis data...")
    df = prepare_survival_data(df)

Your data is now in survival format with clean durations and a clear event indicator. Next up, you will visualize these survival curves using Kaplan-Meier estimation to see how injury timing differs across player positions.

Plot Kaplan-Meier Curves and Run Logrank Tests

Your survival data is ready with duration and event columns for every player. Before building predictive models, you need to understand the baseline injury landscape across the squad.

Kaplan-Meier estimation is a non-parametric method that estimates the probability of staying injury-free over time. By stratifying these curves by position, you can see whether attackers, midfielders, defenders, or goalkeepers face meaningfully different injury timelines. The logrank test then tells you whether those visual differences are statistically significant.

In this step, get ready to:

• Fit a Kaplan-Meier estimator to the full dataset and report the median survival time.
• Stratify survival curves by position group and overlay them on a single plot.
• Run logrank tests to determine if position-based differences are statistically significant.

Fit the overall Kaplan-Meier curve

The overall curve gives you a single view of how quickly players across all positions get injured during a season. The median survival time tells you the day by which half the players have experienced their first injury.

• In soccer_analytics.py, find the comment that says # Step 5: Kaplan-Meier Curves and Logrank Tests.
• Below that section header, add the start of the plot_kaplan_meier function by copying this code:

def plot_kaplan_meier(df):
    """Plot overall and stratified KM curves."""
    kmf = KaplanMeierFitter()

    # Fit on full dataset to get overall survival curve
    kmf.fit(df["duration"], event_observed=df["event"], label="All Players")
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # Plot overall survival on the left panel
    kmf.plot_survival_function(ax=axes[0])
    axes[0].set_title("Overall Injury-Free Survival")
    axes[0].set_xlabel("Days into Season")
    axes[0].set_ylabel("Survival Probability")
    axes[0].set_ylim(0, 1)

    # Report the median survival time
    print(f"Median survival time (overall): {kmf.median_survival_time_:.0f} days")

What does this code do?

KaplanMeierFitter() creates a non-parametric estimator that handles right-censored data. Players who finished the season injury-free (censored observations) are accounted for properly rather than being dropped.

The .fit() method takes your duration (time-to-event) and event (did the event happen?) columns. The median_survival_time_ attribute tells you the day at which 50% of players have been injured.

Stratify survival curves by position group

Overlaying position-specific curves reveals whether certain roles face injury risk earlier in the season. Attackers and midfielders typically endure more physical duels and sprints, which may translate into steeper survival drops.

• Directly below the code you just added (still inside the plot_kaplan_meier function), add the stratified curves and save logic:

    # Stratify by position group on the right panel
    for position in sorted(df["position_group"].unique()):
        mask = df["position_group"] == position
        kmf.fit(df.loc[mask, "duration"], event_observed=df.loc[mask, "event"], label=position)
        kmf.plot_survival_function(ax=axes[1])

    axes[1].set_title("Injury-Free Survival by Position")
    axes[1].set_xlabel("Days into Season")
    axes[1].set_ylabel("Survival Probability")
    axes[1].set_ylim(0, 1)
    axes[1].legend(loc="lower left")

    plt.tight_layout()
    plt.savefig("kaplan_meier.png", dpi=150)
    plt.show()

What does this code do?

The loop re-fits the KaplanMeierFitter once per position group (Attack, Defence, Goalkeeper, Midfield). Each call to plot_survival_function() overlays that position's curve on the same axis, making visual comparison straightforward.

The figure is saved as kaplan_meier.png so you have a persistent record of the plot outside the notebook.

Run logrank tests and interpret the results

Visual differences between curves are not enough on their own. The logrank test is a hypothesis test that determines whether two survival curves are statistically different from each other. A small p-value means the difference is unlikely to have occurred by chance.

• Still inside the plot_kaplan_meier function, add the logrank test loop below the plt.show() line:

    # Logrank tests: compare Attack against each other position
    for pos in ["Midfield", "Defence", "Goalkeeper"]:
        mask_attack = df["position_group"] == "Attack"
        mask_other = df["position_group"] == pos
        result = logrank_test(
            df.loc[mask_attack, "duration"], df.loc[mask_other, "duration"],
            event_observed_A=df.loc[mask_attack, "event"],
            event_observed_B=df.loc[mask_other, "event"],
        )
        print(f"Logrank (Attack vs {pos}): stat={result.test_statistic:.3f}, p={result.p_value:.4f}")

What does this code do?

The logrank_test compares two survival distributions at every observed time point. It returns a test statistic and a p-value.

A p-value below 0.05 indicates a statistically significant difference between the two position groups' injury timelines. You are comparing Attack against each other position to see which comparisons yield meaningful differences.

• Now add the call to plot_kaplan_meier in your main pipeline. Find the line that says df = prepare_survival_data(df) in the if __name__ == "__main__": block, and add these lines below it:

    print("\nPlotting Kaplan-Meier curves...")
    plot_kaplan_meier(df)

• Save soccer_analytics.py.
• Run the script in your Jupyter notebook or terminal. You should see two panels appear: a single survival curve on the left and four overlaid position curves on the right. The terminal prints the median survival time and three logrank test results with test statistics and p-values.

How to interpret your results

On the stratified plot, positions whose curves drop steeply and early (further down and to the left) face higher injury rates sooner in the season. Attackers and midfielders often show steeper drops because they accumulate more high-intensity sprints and physical contacts per match.

Goalkeepers typically show the flattest curve (highest survival) because they face fewer contact injuries. If the logrank test reports a p-value below 0.05 for a comparison, the difference is statistically significant. This information feeds directly into workload management decisions in later steps.

Not seeing the plot?

If you are running the script from the terminal (not a Jupyter cell), plt.show() may not display in some environments. Check that kaplan_meier.png was saved to your project folder and open it directly.

If you see an error about duration or event columns not existing, make sure you ran prepare_survival_data(df) from Step 4 before calling plot_kaplan_meier(df).

Help me debug my Kaplan-Meier plot not rendering

✔️ Awesome, I've got everything!

Great. Double check you have saved soccer_analytics.py and that kaplan_meier.png exists in your project folder.

ⓧ I'd like to double check the full code

import pandas as pd
import numpy as np
import pymc as pm
import arviz as az
from lifelines import KaplanMeierFitter, CoxPHFitter, WeibullAFTFitter
from lifelines.statistics import logrank_test
from lifelines.utils import concordance_index
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime

# =============================================================================
# Step 2: Load, Explore, and Prepare the Data
# =============================================================================

def load_data(profiles_path, performances_path, injuries_path):
    profiles = pd.read_csv(profiles_path)
    performances = pd.read_csv(performances_path)
    injuries = pd.read_csv(injuries_path)
    return profiles, performances, injuries


def clean_and_merge(profiles, performances, injuries, target_season="22/23"):
    injuries["from_date"] = pd.to_datetime(injuries["from_date"], format="%d.%m.%Y", errors="coerce")
    injuries["end_date"] = pd.to_datetime(injuries["end_date"], format="%d.%m.%Y", errors="coerce")
    perf_season = performances[performances["season"] == target_season].copy()
    inj_season = injuries[injuries["season"] == target_season].copy()
    perf_agg = perf_season.groupby("player_id").agg(
        total_minutes=("minutes_played", "sum"),
        total_goals=("goals", "sum"),
        total_assists=("assists", "sum"),
        total_yellow_cards=("yellow_cards", "sum"),
        total_appearances=("nb_on_pitch", "sum"),
        team_id=("team_id", "first"),
        team_name=("team_name", "first"),
    ).reset_index()
    inj_count = inj_season.groupby("player_id").agg(
        injury_count=("from_date", "count"),
        first_injury_date=("from_date", "min"),
        total_days_missed=("days_missed", lambda x: pd.to_numeric(x.str.replace(" days", ""), errors="coerce").sum()),
    ).reset_index()
    merged = profiles[["player_id", "name", "position", "date_of_birth", "height_in_cm"]].copy()
    merged = merged.merge(perf_agg, on="player_id", how="inner")
    merged = merged.merge(inj_count, on="player_id", how="left")
    merged["injury_count"] = merged["injury_count"].fillna(0).astype(int)
    merged["total_days_missed"] = merged["total_days_missed"].fillna(0)
    merged["date_of_birth"] = pd.to_datetime(merged["date_of_birth"], errors="coerce")
    season_start = datetime(2022, 8, 1)
    merged["age"] = ((season_start - merged["date_of_birth"]).dt.days / 365.25).round(1)
    merged = merged[merged["total_minutes"] > 90].copy()
    position_map = {"Attack": "Attack", "Midfield": "Midfield", "Defender": "Defence", "Goalkeeper": "Goalkeeper"}
    merged["position_group"] = merged["position"].map(position_map)
    merged = merged.dropna(subset=["position_group", "age"])
    merged["goals_per_90"] = (merged["total_goals"] / merged["total_minutes"]) * 90
    merged["assists_per_90"] = (merged["total_assists"] / merged["total_minutes"]) * 90
    return merged


def count_previous_injuries(injuries_df, target_season="22/23"):
    seasons_order = ["18/19", "19/20", "20/21", "21/22"]
    prior = injuries_df[injuries_df["season"].isin(seasons_order)]
    prior_counts = prior.groupby("player_id").size().reset_index(name="previous_injuries")
    return prior_counts


# =============================================================================
# Step 3: Bayesian Hierarchical Performance Model
# =============================================================================

def build_hierarchical_model(df):
    teams = df["team_name"].unique()
    team_to_idx = {team: i for i, team in enumerate(teams)}
    df["team_idx"] = df["team_name"].map(team_to_idx)
    positions = df["position_group"].unique()
    pos_to_idx = {pos: i for i, pos in enumerate(positions)}
    df["pos_idx"] = df["position_group"].map(pos_to_idx)
    age_std = (df["age"] - df["age"].mean()) / df["age"].std()
    minutes_std = (df["total_minutes"] - df["total_minutes"].mean()) / df["total_minutes"].std()
    observed_goals_per_90 = df["goals_per_90"].values
    coords = {"team": teams, "position": positions, "obs": np.arange(len(df))}
    with pm.Model(coords=coords) as model:
        sd_team = pm.HalfNormal("sd_team", sigma=0.5)
        team_effect = pm.Normal("team_effect", mu=0, sigma=sd_team, dims="team")
        pos_effect = pm.Normal("pos_effect", mu=0, sigma=1, dims="position")
        beta_age = pm.Normal("beta_age", mu=0, sigma=0.5)
        beta_minutes = pm.Normal("beta_minutes", mu=0, sigma=0.5)
        intercept = pm.Normal("intercept", mu=0.3, sigma=0.5)
        mu = (intercept + team_effect[df["team_idx"].values] + pos_effect[df["pos_idx"].values] + beta_age * age_std.values + beta_minutes * minutes_std.values)
        sigma = pm.HalfNormal("sigma", sigma=0.5)
        goals_obs = pm.Normal("goals_obs", mu=mu, sigma=sigma, observed=observed_goals_per_90, dims="obs")
    return model, df, coords


def sample_hierarchical_model(model, draws=1000, tune=1500, cores=2):
    with model:
        trace = pm.sample(draws=draws, tune=tune, cores=cores, random_seed=42)
    return trace


def analyze_posteriors(trace, df, coords):
    print("=== Team Effect Posterior Summary ===")
    team_summary = az.summary(trace, var_names=["team_effect"])
    print(team_summary.head(10))
    fig, ax = plt.subplots(figsize=(12, 6))
    az.plot_forest(trace, var_names=["team_effect"], ax=ax, combined=True)
    ax.set_title("Team Attack Strength (Posterior 94% HDI)")
    plt.tight_layout()
    plt.savefig("team_effects.png", dpi=150)
    plt.show()
    fig, ax = plt.subplots(figsize=(8, 4))
    az.plot_forest(trace, var_names=["pos_effect"], ax=ax, combined=True)
    ax.set_title("Position Effects on Goals per 90")
    plt.tight_layout()
    plt.savefig("position_effects.png", dpi=150)
    plt.show()
    print("\n=== Convergence Diagnostics ===")
    print(az.summary(trace, var_names=["intercept", "beta_age", "beta_minutes", "sd_team", "sigma"]))
    with trace.model:
        pm.sample_posterior_predictive(trace, extend_inferencedata=True)
    fig, ax = plt.subplots(figsize=(10, 5))
    az.plot_ppc(trace, ax=ax, num_pp_samples=100)
    ax.set_title("Posterior Predictive Check")
    plt.tight_layout()
    plt.savefig("ppc_check.png", dpi=150)
    plt.show()
    return trace


def extract_player_posteriors(trace, df):
    team_effects_mean = trace.posterior["team_effect"].mean(dim=["chain", "draw"]).values
    pos_effects_mean = trace.posterior["pos_effect"].mean(dim=["chain", "draw"]).values
    df["bayesian_performance"] = (team_effects_mean[df["team_idx"].values] + pos_effects_mean[df["pos_idx"].values])
    return df


# =============================================================================
# Step 4: Prepare Survival Analysis Data
# =============================================================================

def prepare_survival_data(df):
    season_length_days = 303
    df["has_injury"] = (df["injury_count"] > 0).astype(int)
    df["minutes_fraction"] = df["total_minutes"] / df["total_minutes"].max()
    df["duration"] = np.where(
        df["has_injury"] == 1,
        np.clip(df["minutes_fraction"] * season_length_days * 0.6, 1, season_length_days),
        season_length_days,
    )
    df["event"] = df["has_injury"]
    return df


# =============================================================================
# Step 5: Kaplan-Meier Curves and Logrank Tests
# =============================================================================

def plot_kaplan_meier(df):
    """Plot overall and stratified KM curves."""
    kmf = KaplanMeierFitter()
    kmf.fit(df["duration"], event_observed=df["event"], label="All Players")
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    kmf.plot_survival_function(ax=axes[0])
    axes[0].set_title("Overall Injury-Free Survival")
    axes[0].set_xlabel("Days into Season")
    axes[0].set_ylabel("Survival Probability")
    axes[0].set_ylim(0, 1)
    print(f"Median survival time (overall): {kmf.median_survival_time_:.0f} days")
    for position in sorted(df["position_group"].unique()):
        mask = df["position_group"] == position
        kmf.fit(df.loc[mask, "duration"], event_observed=df.loc[mask, "event"], label=position)
        kmf.plot_survival_function(ax=axes[1])
    axes[1].set_title("Injury-Free Survival by Position")
    axes[1].set_xlabel("Days into Season")
    axes[1].set_ylabel("Survival Probability")
    axes[1].set_ylim(0, 1)
    axes[1].legend(loc="lower left")
    plt.tight_layout()
    plt.savefig("kaplan_meier.png", dpi=150)
    plt.show()
    for pos in ["Midfield", "Defence", "Goalkeeper"]:
        mask_attack = df["position_group"] == "Attack"
        mask_other = df["position_group"] == pos
        result = logrank_test(
            df.loc[mask_attack, "duration"], df.loc[mask_other, "duration"],
            event_observed_A=df.loc[mask_attack, "event"],
            event_observed_B=df.loc[mask_other, "event"],
        )
        print(f"Logrank (Attack vs {pos}): stat={result.test_statistic:.3f}, p={result.p_value:.4f}")


# =============================================================================
# Main Pipeline
# =============================================================================

if __name__ == "__main__":
    PROFILES_PATH = "data/player_profiles.csv"
    PERFORMANCES_PATH = "data/player_performances.csv"
    INJURIES_PATH = "data/player_injury_histories.csv"

    print("Loading data...")
    profiles, performances, injuries = load_data(PROFILES_PATH, PERFORMANCES_PATH, INJURIES_PATH)

    print("Cleaning and merging...")
    df = clean_and_merge(profiles, performances, injuries, target_season="22/23")
    prior_injuries = count_previous_injuries(injuries, target_season="22/23")
    df = df.merge(prior_injuries, on="player_id", how="left")
    df["previous_injuries"] = df["previous_injuries"].fillna(0)
    print(f"Merged dataset: {df.shape[0]} players")

    print("\nBuilding Bayesian hierarchical model...")
    model, df, coords = build_hierarchical_model(df)
    print("Sampling posterior (this may take a few minutes)...")
    trace = sample_hierarchical_model(model, draws=1000, tune=1500, cores=2)
    trace = analyze_posteriors(trace, df, coords)
    df = extract_player_posteriors(trace, df)

    print("\nPreparing survival analysis data...")
    df = prepare_survival_data(df)

    print("\nPlotting Kaplan-Meier curves...")
    plot_kaplan_meier(df)

You now have a clear picture of how injury risk varies across positions over a season. Next up, you will fit Cox Proportional Hazards and Weibull AFT models to quantify which covariates drive injury risk and generate time-specific predictions.

Fit Cox PH and Weibull AFT Models

Your Kaplan-Meier curves revealed overall survival patterns and position-level differences. Now you need models that quantify exactly which factors drive injury risk and can predict it for individual players.

In this step, you'll fit two competing survival regression models. The Cox Proportional Hazards model estimates relative risk factors without assuming a distribution shape. The Weibull AFT model assumes a Weibull distribution and provides absolute time predictions. Comparing their concordance indices tells you which better fits your data.

In this step, get ready to:

• Fit a Cox PH model with regularization and report its hazard ratios and C-index.
• Fit a Weibull AFT model and compare its C-index against the Cox model.
• Validate the proportional hazards assumption and visualize both models' coefficients.

Prepare covariates and fit the Cox PH model

Before fitting survival models, you need to transform your DataFrame into a format they expect. This means selecting the relevant columns, creating position dummy variables, and standardizing continuous covariates so coefficients are comparable.

• In soccer_analytics.py, find the comment block that says # Step 6: Cox PH and Weibull AFT Models with C-Index (add this below the plot_kaplan_meier function).
• Add the prepare_survival_covariates function by copying the following code:

# =============================================================================
# Step 6: Cox PH and Weibull AFT Models with C-Index
# =============================================================================

def prepare_survival_covariates(df):
    """Prepare covariate DataFrame for survival models."""
    surv_df = df[["duration", "event", "age", "total_minutes", "previous_injuries",
                  "bayesian_performance", "total_yellow_cards"]].copy()
    position_dummies = pd.get_dummies(df["position_group"], prefix="pos", drop_first=True)
    surv_df = pd.concat([surv_df, position_dummies], axis=1)
    surv_df = surv_df.dropna()
    # Standardize continuous covariates
    for col in ["age", "total_minutes", "previous_injuries", "bayesian_performance", "total_yellow_cards"]:
        mean_val = surv_df[col].mean()
        std_val = surv_df[col].std()
        if std_val > 0:
            surv_df[col] = (surv_df[col] - mean_val) / std_val
    return surv_df

What does this code do?

This function transforms your raw DataFrame into the format survival models expect.

• It selects only the columns the models need: duration, event indicator, and the covariates (age, minutes, injuries, Bayesian performance, yellow cards).
• pd.get_dummies(..., drop_first=True) converts position groups into binary dummy columns, dropping one category (Attack) to avoid multicollinearity.
• Standardizing continuous covariates (subtracting the mean, dividing by standard deviation) ensures that coefficient magnitudes are directly comparable across features.

The Cox PH model estimates how each covariate multiplies the baseline hazard of injury. A hazard ratio above 1 means higher injury risk.

• Add the fit_cox_model function directly below prepare_survival_covariates:

def fit_cox_model(surv_df):
    """Fit Cox PH model and report C-index."""
    cph = CoxPHFitter(penalizer=0.01)
    cph.fit(surv_df, duration_col="duration", event_col="event")
    print("=== Cox Proportional Hazards Model ===")
    cph.print_summary()
    print(f"\nCox C-index: {cph.score(surv_df, scoring_method='concordance_index'):.4f}")
    return cph

What does this code do?

• CoxPHFitter(penalizer=0.01) creates a Cox model with L2 regularization. The penalizer prevents overfitting when covariates are correlated.
• cph.fit() estimates partial likelihood coefficients. The summary shows hazard ratios, confidence intervals, and p-values for each covariate.
• cph.score(surv_df, scoring_method='concordance_index') reports the C-index, which measures how well the model ranks patients by risk. A value of 0.5 is random guessing. Above 0.7 is considered useful.

• Save your file. You won't run this yet. The next chunk adds the Weibull model, and then you'll run both together.

Seeing a syntax error in the function?

• Check that the f-string uses single quotes inside the double-quoted format string for the scoring_method parameter.
• Make sure the function is indented at the module level (no extra indentation from a previous block).

Help me debug my fit_cox_model function in lifelines

Fit the Weibull AFT model

The Weibull Accelerated Failure Time model takes a different approach. Instead of modeling relative hazards, it models how covariates accelerate or decelerate the time until injury occurs.

• Add the fit_weibull_model function directly below fit_cox_model:

def fit_weibull_model(surv_df):
    """Fit Weibull AFT model and report C-index."""
    wft = WeibullAFTFitter(penalizer=0.01)
    wft.fit(surv_df, duration_col="duration", event_col="event")
    print("\n=== Weibull AFT Model ===")
    wft.print_summary()
    print(f"\nWeibull C-index: {wft.score(surv_df, scoring_method='concordance_index'):.4f}")
    return wft

How is Weibull AFT different from Cox PH?

The Cox model is semi-parametric. It estimates relative hazard ratios without assuming the shape of the underlying risk distribution. The Weibull AFT model is fully parametric. It assumes survival times follow a Weibull distribution controlled by a scale parameter (lambda) and a shape parameter (rho).

The Weibull model's advantage is that it can make absolute time-to-event predictions ("this player has a 25% chance of injury within 14 days"). The Cox model only tells you relative risk compared to a baseline.

• Save soccer_analytics.py. You won't see output yet. The next chunk compares both models and produces the visualization.

Getting a convergence warning?

• The penalizer=0.01 should handle most convergence issues. If you still see warnings, try increasing it to 0.1.
• Make sure your standardization step ran correctly. Non-standardized covariates with very different scales can cause optimization issues.

Help me fix convergence issues with my WeibullAFTFitter

Compare models and validate assumptions

With both models fitted, you can compare their predictive performance using the concordance index. You also need to validate whether the Cox model's key assumption (proportional hazards) holds for your data.

• Add the compare_models function directly below fit_weibull_model:

def compare_models(cph, wft, surv_df):
    """Compare Cox and Weibull models side by side."""
    cox_c = cph.score(surv_df, scoring_method="concordance_index")
    weibull_c = wft.score(surv_df, scoring_method="concordance_index")
    print("\n=== Model Comparison ===")
    print(f"Cox PH C-index:     {cox_c:.4f}")
    print(f"Weibull AFT C-index: {weibull_c:.4f}")
    print(f"Better model: {'Cox PH' if cox_c > weibull_c else 'Weibull AFT'}")
    # Plot coefficients side by side
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    cph.plot(ax=axes[0])
    axes[0].set_title(f"Cox PH Hazard Ratios (C={cox_c:.3f})")
    wft.plot(ax=axes[1])
    axes[1].set_title(f"Weibull AFT Coefficients (C={weibull_c:.3f})")
    plt.tight_layout()
    plt.savefig("model_comparison.png", dpi=150)
    plt.show()
    # Check proportional hazards assumption for Cox
    print("\n=== Proportional Hazards Check ===")
    try:
        cph.check_assumptions(surv_df, p_value_threshold=0.05)
    except Exception as e:
        print(f"Note: {e}")
    return cox_c, weibull_c

What does this code do?

• It computes the C-index for both models and prints which one has better discriminative ability.
• cph.plot() and wft.plot() generate forest plots showing coefficient estimates with confidence intervals, placed side by side for visual comparison.
• cph.check_assumptions() tests whether the proportional hazards assumption holds. If it fails for a covariate, the Cox model's hazard ratio for that covariate may not be reliable over time.

Why semi-parametric vs. parametric?

The Cox model makes fewer assumptions about the data, which makes it more robust. But the Weibull model's parametric form lets it extrapolate predictions into the future. In the next step, you'll use both models' survival function predictions together for a more robust ensemble estimate.

Now you need to wire these functions into your main pipeline so they actually run.

• Scroll down to the if __name__ == "__main__": block. Find the line that says plot_kaplan_meier(df).
• Add the following code directly below it:

    print("\nFitting survival models...")
    surv_df = prepare_survival_covariates(df)
    cph = fit_cox_model(surv_df)
    wft = fit_weibull_model(surv_df)
    cox_c, weibull_c = compare_models(cph, wft, surv_df)

What does this code do?

This calls your four new functions in sequence. First it prepares the covariate DataFrame, then fits both models, and finally compares their C-indices and validates assumptions.

• Save soccer_analytics.py.
• Run the full script in your terminal by running:

python soccer_analytics.py

• You should see the Cox model summary with hazard ratios printed, then the Weibull model summary with lambda/rho parameters, then a model comparison showing both C-indices. A two-panel figure (model_comparison.png) should appear showing coefficient forest plots for both models.

Not seeing the model summaries?

• Check that the main pipeline lines are indented inside the if __name__ == "__main__": block (one level of indentation).
• If you see a KeyError for a column name, confirm that the prior steps (Bayesian model + survival data prep) ran successfully. The bayesian_performance column must exist before this step.
• If the proportional hazards check prints a note or exception, that is expected behavior. The try/except block handles cases where the assumption test encounters issues.

Help me debug my Cox PH and Weibull survival model pipeline

✔️ Awesome, I've got everything!

Great. Double check you've saved your file and that model_comparison.png appeared in your project folder.

ⓧ I'd like to double check the full code

Here is the complete soccer_analytics.py file as it should look after this step:

import pandas as pd
import numpy as np
import pymc as pm
import arviz as az
from lifelines import KaplanMeierFitter, CoxPHFitter, WeibullAFTFitter
from lifelines.statistics import logrank_test
from lifelines.utils import concordance_index
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime

# =============================================================================
# Step 2: Load, Explore, and Prepare the Data
# =============================================================================

def load_data(profiles_path, performances_path, injuries_path):
    profiles = pd.read_csv(profiles_path)
    performances = pd.read_csv(performances_path)
    injuries = pd.read_csv(injuries_path)
    return profiles, performances, injuries


def clean_and_merge(profiles, performances, injuries, target_season="22/23"):
    injuries["from_date"] = pd.to_datetime(injuries["from_date"], format="%d.%m.%Y", errors="coerce")
    injuries["end_date"] = pd.to_datetime(injuries["end_date"], format="%d.%m.%Y", errors="coerce")
    perf_season = performances[performances["season"] == target_season].copy()
    inj_season = injuries[injuries["season"] == target_season].copy()
    perf_agg = perf_season.groupby("player_id").agg(
        total_minutes=("minutes_played", "sum"),
        total_goals=("goals", "sum"),
        total_assists=("assists", "sum"),
        total_yellow_cards=("yellow_cards", "sum"),
        total_appearances=("nb_on_pitch", "sum"),
        team_id=("team_id", "first"),
        team_name=("team_name", "first"),
    ).reset_index()
    inj_count = inj_season.groupby("player_id").agg(
        injury_count=("from_date", "count"),
        first_injury_date=("from_date", "min"),
        total_days_missed=("days_missed", lambda x: pd.to_numeric(x.str.replace(" days", ""), errors="coerce").sum()),
    ).reset_index()
    merged = profiles[["player_id", "name", "position", "date_of_birth", "height_in_cm"]].copy()
    merged = merged.merge(perf_agg, on="player_id", how="inner")
    merged = merged.merge(inj_count, on="player_id", how="left")
    merged["injury_count"] = merged["injury_count"].fillna(0).astype(int)
    merged["total_days_missed"] = merged["total_days_missed"].fillna(0)
    merged["date_of_birth"] = pd.to_datetime(merged["date_of_birth"], errors="coerce")
    season_start = datetime(2022, 8, 1)
    merged["age"] = ((season_start - merged["date_of_birth"]).dt.days / 365.25).round(1)
    merged = merged[merged["total_minutes"] > 90].copy()
    position_map = {"Attack": "Attack", "Midfield": "Midfield", "Defender": "Defence", "Goalkeeper": "Goalkeeper"}
    merged["position_group"] = merged["position"].map(position_map)
    merged = merged.dropna(subset=["position_group", "age"])
    merged["goals_per_90"] = (merged["total_goals"] / merged["total_minutes"]) * 90
    merged["assists_per_90"] = (merged["total_assists"] / merged["total_minutes"]) * 90
    return merged


def count_previous_injuries(injuries_df, target_season="22/23"):
    seasons_order = ["18/19", "19/20", "20/21", "21/22"]
    prior = injuries_df[injuries_df["season"].isin(seasons_order)]
    prior_counts = prior.groupby("player_id").size().reset_index(name="previous_injuries")
    return prior_counts


# =============================================================================
# Step 3: Bayesian Hierarchical Performance Model
# =============================================================================

def build_hierarchical_model(df):
    teams = df["team_name"].unique()
    team_to_idx = {team: i for i, team in enumerate(teams)}
    df["team_idx"] = df["team_name"].map(team_to_idx)
    positions = df["position_group"].unique()
    pos_to_idx = {pos: i for i, pos in enumerate(positions)}
    df["pos_idx"] = df["position_group"].map(pos_to_idx)
    age_std = (df["age"] - df["age"].mean()) / df["age"].std()
    minutes_std = (df["total_minutes"] - df["total_minutes"].mean()) / df["total_minutes"].std()
    observed_goals_per_90 = df["goals_per_90"].values
    coords = {"team": teams, "position": positions, "obs": np.arange(len(df))}
    with pm.Model(coords=coords) as model:
        sd_team = pm.HalfNormal("sd_team", sigma=0.5)
        team_effect = pm.Normal("team_effect", mu=0, sigma=sd_team, dims="team")
        pos_effect = pm.Normal("pos_effect", mu=0, sigma=1, dims="position")
        beta_age = pm.Normal("beta_age", mu=0, sigma=0.5)
        beta_minutes = pm.Normal("beta_minutes", mu=0, sigma=0.5)
        intercept = pm.Normal("intercept", mu=0.3, sigma=0.5)
        mu = (intercept + team_effect[df["team_idx"].values] + pos_effect[df["pos_idx"].values] + beta_age * age_std.values + beta_minutes * minutes_std.values)
        sigma = pm.HalfNormal("sigma", sigma=0.5)
        goals_obs = pm.Normal("goals_obs", mu=mu, sigma=sigma, observed=observed_goals_per_90, dims="obs")
    return model, df, coords


def sample_hierarchical_model(model, draws=1000, tune=1500, cores=2):
    with model:
        trace = pm.sample(draws=draws, tune=tune, cores=cores, random_seed=42)
    return trace


def analyze_posteriors(trace, df, coords):
    print("=== Team Effect Posterior Summary ===")
    team_summary = az.summary(trace, var_names=["team_effect"])
    print(team_summary.head(10))
    fig, ax = plt.subplots(figsize=(12, 6))
    az.plot_forest(trace, var_names=["team_effect"], ax=ax, combined=True)
    ax.set_title("Team Attack Strength (Posterior 94% HDI)")
    plt.tight_layout()
    plt.savefig("team_effects.png", dpi=150)
    plt.show()
    fig, ax = plt.subplots(figsize=(8, 4))
    az.plot_forest(trace, var_names=["pos_effect"], ax=ax, combined=True)
    ax.set_title("Position Effects on Goals per 90")
    plt.tight_layout()
    plt.savefig("position_effects.png", dpi=150)
    plt.show()
    print("\n=== Convergence Diagnostics ===")
    print(az.summary(trace, var_names=["intercept", "beta_age", "beta_minutes", "sd_team", "sigma"]))
    with trace.model:
        pm.sample_posterior_predictive(trace, extend_inferencedata=True)
    fig, ax = plt.subplots(figsize=(10, 5))
    az.plot_ppc(trace, ax=ax, num_pp_samples=100)
    ax.set_title("Posterior Predictive Check")
    plt.tight_layout()
    plt.savefig("ppc_check.png", dpi=150)
    plt.show()
    return trace


def extract_player_posteriors(trace, df):
    team_effects_mean = trace.posterior["team_effect"].mean(dim=["chain", "draw"]).values
    pos_effects_mean = trace.posterior["pos_effect"].mean(dim=["chain", "draw"]).values
    df["bayesian_performance"] = (team_effects_mean[df["team_idx"].values] + pos_effects_mean[df["pos_idx"].values])
    return df


# =============================================================================
# Step 4: Prepare Survival Analysis Data
# =============================================================================

def prepare_survival_data(df):
    season_length_days = 303
    df["has_injury"] = (df["injury_count"] > 0).astype(int)
    df["minutes_fraction"] = df["total_minutes"] / df["total_minutes"].max()
    df["duration"] = np.where(
        df["has_injury"] == 1,
        np.clip(df["minutes_fraction"] * season_length_days * 0.6, 1, season_length_days),
        season_length_days,
    )
    df["event"] = df["has_injury"]
    return df


# =============================================================================
# Step 5: Kaplan-Meier Curves and Logrank Tests
# =============================================================================

def plot_kaplan_meier(df):
    kmf = KaplanMeierFitter()
    kmf.fit(df["duration"], event_observed=df["event"], label="All Players")
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    kmf.plot_survival_function(ax=axes[0])
    axes[0].set_title("Overall Injury-Free Survival")
    axes[0].set_xlabel("Days into Season")
    axes[0].set_ylabel("Survival Probability")
    axes[0].set_ylim(0, 1)
    print(f"Median survival time (overall): {kmf.median_survival_time_:.0f} days")
    for position in sorted(df["position_group"].unique()):
        mask = df["position_group"] == position
        kmf.fit(df.loc[mask, "duration"], event_observed=df.loc[mask, "event"], label=position)
        kmf.plot_survival_function(ax=axes[1])
    axes[1].set_title("Injury-Free Survival by Position")
    axes[1].set_xlabel("Days into Season")
    axes[1].set_ylabel("Survival Probability")
    axes[1].set_ylim(0, 1)
    axes[1].legend(loc="lower left")
    plt.tight_layout()
    plt.savefig("kaplan_meier.png", dpi=150)
    plt.show()
    for pos in ["Midfield", "Defence", "Goalkeeper"]:
        mask_attack = df["position_group"] == "Attack"
        mask_other = df["position_group"] == pos
        result = logrank_test(
            df.loc[mask_attack, "duration"], df.loc[mask_other, "duration"],
            event_observed_A=df.loc[mask_attack, "event"],
            event_observed_B=df.loc[mask_other, "event"],
        )
        print(f"Logrank (Attack vs {pos}): stat={result.test_statistic:.3f}, p={result.p_value:.4f}")


# =============================================================================
# Step 6: Cox PH and Weibull AFT Models with C-Index
# =============================================================================

def prepare_survival_covariates(df):
    """Prepare covariate DataFrame for survival models."""
    surv_df = df[["duration", "event", "age", "total_minutes", "previous_injuries",
                  "bayesian_performance", "total_yellow_cards"]].copy()
    position_dummies = pd.get_dummies(df["position_group"], prefix="pos", drop_first=True)
    surv_df = pd.concat([surv_df, position_dummies], axis=1)
    surv_df = surv_df.dropna()
    for col in ["age", "total_minutes", "previous_injuries", "bayesian_performance", "total_yellow_cards"]:
        mean_val = surv_df[col].mean()
        std_val = surv_df[col].std()
        if std_val > 0:
            surv_df[col] = (surv_df[col] - mean_val) / std_val
    return surv_df


def fit_cox_model(surv_df):
    """Fit Cox PH model and report C-index."""
    cph = CoxPHFitter(penalizer=0.01)
    cph.fit(surv_df, duration_col="duration", event_col="event")
    print("=== Cox Proportional Hazards Model ===")
    cph.print_summary()
    print(f"\nCox C-index: {cph.score(surv_df, scoring_method='concordance_index'):.4f}")
    return cph


def fit_weibull_model(surv_df):
    """Fit Weibull AFT model and report C-index."""
    wft = WeibullAFTFitter(penalizer=0.01)
    wft.fit(surv_df, duration_col="duration", event_col="event")
    print("\n=== Weibull AFT Model ===")
    wft.print_summary()
    print(f"\nWeibull C-index: {wft.score(surv_df, scoring_method='concordance_index'):.4f}")
    return wft


def compare_models(cph, wft, surv_df):
    """Compare Cox and Weibull models side by side."""
    cox_c = cph.score(surv_df, scoring_method="concordance_index")
    weibull_c = wft.score(surv_df, scoring_method="concordance_index")
    print("\n=== Model Comparison ===")
    print(f"Cox PH C-index:     {cox_c:.4f}")
    print(f"Weibull AFT C-index: {weibull_c:.4f}")
    print(f"Better model: {'Cox PH' if cox_c > weibull_c else 'Weibull AFT'}")
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    cph.plot(ax=axes[0])
    axes[0].set_title(f"Cox PH Hazard Ratios (C={cox_c:.3f})")
    wft.plot(ax=axes[1])
    axes[1].set_title(f"Weibull AFT Coefficients (C={weibull_c:.3f})")
    plt.tight_layout()
    plt.savefig("model_comparison.png", dpi=150)
    plt.show()
    print("\n=== Proportional Hazards Check ===")
    try:
        cph.check_assumptions(surv_df, p_value_threshold=0.05)
    except Exception as e:
        print(f"Note: {e}")
    return cox_c, weibull_c


# =============================================================================
# Main Pipeline
# =============================================================================

if __name__ == "__main__":
    PROFILES_PATH = "data/player_profiles.csv"
    PERFORMANCES_PATH = "data/player_performances.csv"
    INJURIES_PATH = "data/player_injury_histories.csv"

    print("Loading data...")
    profiles, performances, injuries = load_data(PROFILES_PATH, PERFORMANCES_PATH, INJURIES_PATH)

    print("Cleaning and merging...")
    df = clean_and_merge(profiles, performances, injuries, target_season="22/23")
    prior_injuries = count_previous_injuries(injuries, target_season="22/23")
    df = df.merge(prior_injuries, on="player_id", how="left")
    df["previous_injuries"] = df["previous_injuries"].fillna(0)
    print(f"Merged dataset: {df.shape[0]} players")

    print("\nBuilding Bayesian hierarchical model...")
    model, df, coords = build_hierarchical_model(df)
    print("Sampling posterior (this may take a few minutes)...")
    trace = sample_hierarchical_model(model, draws=1000, tune=1500, cores=2)
    trace = analyze_posteriors(trace, df, coords)
    df = extract_player_posteriors(trace, df)

    print("\nPreparing survival analysis data...")
    df = prepare_survival_data(df)

    print("\nPlotting Kaplan-Meier curves...")
    plot_kaplan_meier(df)

    print("\nFitting survival models...")
    surv_df = prepare_survival_covariates(df)
    cph = fit_cox_model(surv_df)
    wft = fit_weibull_model(surv_df)
    cox_c, weibull_c = compare_models(cph, wft, surv_df)

You now have two fitted survival models with quantified predictive performance. Next up, you'll combine these models with your Bayesian performance posteriors to build a substitution recommendation system that supports real squad management decisions.

Build the Substitution Recommendation System

You now have two fitted survival models that predict when injuries are likely to occur. The next question is practical: what should a coaching staff actually do with these predictions?

In this step, you will build a substitution recommendation system that combines your Bayesian performance posteriors with the survival model predictions. The system estimates short-term injury probability for each player and flags high-risk individuals who may need rest or rotation.

In this step, get ready to:

• Create an ensemble injury prediction function averaging Cox and Weibull survival estimates.
• Build a squad report with risk classifications and Bayesian performance integration.
• Visualize the squad's risk-performance profile as a scatter plot and heatmap.

Create the ensemble injury prediction function

A single model's prediction can be overconfident or biased by its assumptions. Averaging predictions from the semi-parametric Cox model and the parametric Weibull model produces a more robust estimate that hedges against either model's weaknesses.

• In soccer_analytics.py, find the Step 6 section (after the compare_models() function).
• Add the following section comment and function below it:

# =============================================================================
# Step 7: Substitution Recommendation System
# =============================================================================

def predict_injury_risk(cph, wft, player_data):
    """Ensemble prediction from Cox and Weibull models."""
    # Get survival curves from both models
    cox_surv = cph.predict_survival_function(player_data)
    weibull_surv = wft.predict_survival_function(player_data)

    # Use the player's current duration as the starting point
    current_day = player_data["duration"].values[0]
    risks = {}

    # Compute risk at 1-day, 7-day, and 14-day horizons
    for horizon_name, horizon_days in [("1_day", 1), ("7_day", 7), ("14_day", 14)]:
        target_day = min(current_day + horizon_days, cox_surv.index.max())
        # Risk = 1 - survival probability at target day
        cox_risk = 1 - cox_surv.loc[target_day].values[0]
        weibull_risk = 1 - weibull_surv.loc[target_day].values[0]
        # Average both models for robust ensemble
        risks[horizon_name] = (cox_risk + weibull_risk) / 2

    return risks

What does this code do?

The predict_injury_risk() function takes a single player's covariate row and queries both fitted survival models.

• predict_survival_function() returns the probability of remaining injury-free at each future time point.
• Subtracting from 1 converts survival probability into injury probability.
• Averaging the Cox and Weibull predictions creates an ensemble that is more reliable than either model alone.

• Save soccer_analytics.py.

Seeing a KeyError on target_day?

• Make sure the duration column is present in player_data. The function reads from it to set the starting day.
• If the target day exceeds the survival curve index, the min() clamp handles it. Check that your survival models were fitted on the same surv_df you are passing player rows from.

Help me debug a KeyError in my predict_injury_risk function when looking up target_day in the survival curve index.

Build risk classification and squad report generation

Raw probabilities are hard to act on. Classifying players into risk tiers gives coaching staff a clear decision framework: rest, monitor, or deploy.

• Add the following two functions below predict_injury_risk() in soccer_analytics.py:

def classify_risk(risk_7d):
    """Classify injury risk into actionable categories."""
    if risk_7d > 0.30:
        return "HIGH RISK - Recommend Rest"
    elif risk_7d > 0.15:
        return "MONITOR - Consider Rotation"
    else:
        return "LOW RISK - Available"


def generate_squad_report(cph, wft, surv_df, df):
    """Generate risk-performance report for the squad."""
    recommendations = []

    # Iterate over first 25 players in the survival DataFrame
    for idx in surv_df.index[:25]:
        player_row = surv_df.loc[[idx]]
        try:
            risks = predict_injury_risk(cph, wft, player_row)
            recommendation = classify_risk(risks["7_day"])

            # Get Bayesian performance from original df
            perf = df.loc[idx, "bayesian_performance"] if idx in df.index else 0

            recommendations.append({
                "player_idx": idx,
                "performance_posterior": perf,
                "risk_1d": risks["1_day"],
                "risk_7d": risks["7_day"],
                "risk_14d": risks["14_day"],
                "recommendation": recommendation,
            })
        except Exception:
            continue

    rec_df = pd.DataFrame(recommendations)
    return rec_df

What does this code do?

• classify_risk() maps a 7-day injury probability to one of three tiers: above 30% is high risk, 15-30% warrants monitoring, and below 15% means the player is available.
• generate_squad_report() loops over the first 25 players, calls both prediction functions, and retrieves each player's Bayesian performance posterior from the original df.
• The try/except block gracefully skips players whose covariate profile causes edge-case errors in survival prediction.

Why combine performance posteriors with injury risk?

A high-risk player who is also your top scorer requires a different decision than a high-risk substitute. By placing both dimensions side by side, the system surfaces the tradeoff: how much performance do you sacrifice by resting this player?

• Save soccer_analytics.py.

Getting an empty rec_df?

• Check that surv_df has at least 25 rows. If your filtered dataset is smaller, reduce the slice or check your earlier filtering steps.
• Make sure the index of surv_df aligns with df. The idx in df.index check handles mismatches gracefully, but if all fail you will get zeros for performance.

Help me figure out why my generate_squad_report is returning an empty DataFrame.

Visualize risk-performance and generate the squad report

The final piece is a two-panel visualization. A scatter plot reveals which players sit in the dangerous quadrant of high performance and high risk. A heatmap gives the coaching staff a quick glance at upcoming risk across multiple time horizons.

• Add the visualization function below generate_squad_report() in soccer_analytics.py:

def plot_risk_performance_scatter(rec_df):
    """Scatter plot: performance vs injury risk."""
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # Left panel: scatter of performance vs 7-day risk
    colors = rec_df["risk_7d"].values
    scatter = axes[0].scatter(
        rec_df["performance_posterior"], rec_df["risk_7d"],
        c=colors, cmap="RdYlGn_r", s=60, edgecolors="black", linewidth=0.5
    )
    axes[0].set_xlabel("Bayesian Performance (Posterior Mean)")
    axes[0].set_ylabel("7-Day Injury Risk")
    axes[0].set_title("Performance vs. Injury Risk")
    axes[0].axhline(y=0.30, color="red", linestyle="--", alpha=0.7, label="High Risk Threshold")
    axes[0].axhline(y=0.15, color="orange", linestyle="--", alpha=0.7, label="Monitor Threshold")
    axes[0].legend()
    plt.colorbar(scatter, ax=axes[0], label="Risk")

    # Right panel: heatmap of multi-horizon risk
    risk_matrix = rec_df[["risk_1d", "risk_7d", "risk_14d"]].copy()
    risk_matrix.columns = ["Next 1 Day", "Next 7 Days", "Next 14 Days"]
    risk_matrix.index = [f"Player {i+1}" for i in range(len(risk_matrix))]
    sns.heatmap(risk_matrix.head(15), annot=True, fmt=".2f", cmap="RdYlGn_r",
                vmin=0, vmax=0.5, ax=axes[1])
    axes[1].set_title("Squad Injury Risk Heatmap")

    plt.tight_layout()
    plt.savefig("squad_risk_performance.png", dpi=150)
    plt.show()

What does this visualization show?

• The scatter plot places each player on two axes: how good they are (x) versus how likely they are to get injured in 7 days (y). Players in the top-right quadrant are the hardest decisions.
• The red dashed line at 0.30 and orange dashed line at 0.15 correspond to the risk thresholds defined in classify_risk().
• The heatmap shows predicted risk across three time horizons for the first 15 players, letting staff see how risk escalates over time.

Now add the main pipeline calls to run the recommendation system. These go at the end of the existing main block.

• In the if __name__ == "__main__": block, find the line that says cox_c, weibull_c = compare_models(cph, wft, surv_df).
• Add the following code below it:

    # Step 7: Generate squad recommendations
    print("\nGenerating squad recommendations...")
    rec_df = generate_squad_report(cph, wft, surv_df, df)
    print(rec_df[["performance_posterior", "risk_7d", "recommendation"]].to_string())
    plot_risk_performance_scatter(rec_df)

What does this code do?

• generate_squad_report() runs the ensemble prediction for the first 25 players and builds the recommendation DataFrame.
• The print() call displays performance, 7-day risk, and the recommendation tier for each player in the terminal.
• plot_risk_performance_scatter() generates the two-panel figure and saves it as squad_risk_performance.png.

• Save soccer_analytics.py.
• Run the full pipeline in your Jupyter notebook or terminal. After the survival model output, you should see a table of player recommendations with columns showing performance_posterior, risk_7d, and recommendation.
• You should also see a two-panel figure appear: a scatter plot on the left with colored dots and dashed threshold lines, and a heatmap on the right showing risk values across time horizons.

Plot not showing or file not saved?

• If you are running in a Jupyter notebook, plt.show() renders inline. If running as a script, the figure still saves to squad_risk_performance.png in your working directory.
• If the plot is blank, confirm that rec_df is not empty by printing rec_df.shape before the plot call.

Help me troubleshoot why my squad risk-performance plot is not rendering correctly.

✔️ Awesome, I've got everything!

Great. Double check you have saved your file and can see the recommendation output in your terminal or notebook.

ⓧ I'd like to double check the full code

Here is the complete soccer_analytics.py file including all Step 7 additions:

import pandas as pd
import numpy as np
import pymc as pm
import arviz as az
from lifelines import KaplanMeierFitter, CoxPHFitter, WeibullAFTFitter
from lifelines.statistics import logrank_test
from lifelines.utils import concordance_index
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime

# =============================================================================
# Step 2: Load, Explore, and Prepare the Data
# =============================================================================

def load_data(profiles_path, performances_path, injuries_path):
    profiles = pd.read_csv(profiles_path)
    performances = pd.read_csv(performances_path)
    injuries = pd.read_csv(injuries_path)
    return profiles, performances, injuries


def clean_and_merge(profiles, performances, injuries, target_season="22/23"):
    injuries["from_date"] = pd.to_datetime(injuries["from_date"], format="%d.%m.%Y", errors="coerce")
    injuries["end_date"] = pd.to_datetime(injuries["end_date"], format="%d.%m.%Y", errors="coerce")
    perf_season = performances[performances["season"] == target_season].copy()
    inj_season = injuries[injuries["season"] == target_season].copy()
    perf_agg = perf_season.groupby("player_id").agg(
        total_minutes=("minutes_played", "sum"),
        total_goals=("goals", "sum"),
        total_assists=("assists", "sum"),
        total_yellow_cards=("yellow_cards", "sum"),
        total_appearances=("nb_on_pitch", "sum"),
        team_id=("team_id", "first"),
        team_name=("team_name", "first"),
    ).reset_index()
    inj_count = inj_season.groupby("player_id").agg(
        injury_count=("from_date", "count"),
        first_injury_date=("from_date", "min"),
        total_days_missed=("days_missed", lambda x: pd.to_numeric(x.str.replace(" days", ""), errors="coerce").sum()),
    ).reset_index()
    merged = profiles[["player_id", "name", "position", "date_of_birth", "height_in_cm"]].copy()
    merged = merged.merge(perf_agg, on="player_id", how="inner")
    merged = merged.merge(inj_count, on="player_id", how="left")
    merged["injury_count"] = merged["injury_count"].fillna(0).astype(int)
    merged["total_days_missed"] = merged["total_days_missed"].fillna(0)
    merged["date_of_birth"] = pd.to_datetime(merged["date_of_birth"], errors="coerce")
    season_start = datetime(2022, 8, 1)
    merged["age"] = ((season_start - merged["date_of_birth"]).dt.days / 365.25).round(1)
    merged = merged[merged["total_minutes"] > 90].copy()
    position_map = {"Attack": "Attack", "Midfield": "Midfield", "Defender": "Defence", "Goalkeeper": "Goalkeeper"}
    merged["position_group"] = merged["position"].map(position_map)
    merged = merged.dropna(subset=["position_group", "age"])
    merged["goals_per_90"] = (merged["total_goals"] / merged["total_minutes"]) * 90
    merged["assists_per_90"] = (merged["total_assists"] / merged["total_minutes"]) * 90
    return merged


def count_previous_injuries(injuries_df, target_season="22/23"):
    seasons_order = ["18/19", "19/20", "20/21", "21/22"]
    prior = injuries_df[injuries_df["season"].isin(seasons_order)]
    prior_counts = prior.groupby("player_id").size().reset_index(name="previous_injuries")
    return prior_counts


# =============================================================================
# Step 3: Bayesian Hierarchical Performance Model
# =============================================================================

def build_hierarchical_model(df):
    teams = df["team_name"].unique()
    team_to_idx = {team: i for i, team in enumerate(teams)}
    df["team_idx"] = df["team_name"].map(team_to_idx)
    positions = df["position_group"].unique()
    pos_to_idx = {pos: i for i, pos in enumerate(positions)}
    df["pos_idx"] = df["position_group"].map(pos_to_idx)
    age_std = (df["age"] - df["age"].mean()) / df["age"].std()
    minutes_std = (df["total_minutes"] - df["total_minutes"].mean()) / df["total_minutes"].std()
    observed_goals_per_90 = df["goals_per_90"].values
    coords = {"team": teams, "position": positions, "obs": np.arange(len(df))}
    with pm.Model(coords=coords) as model:
        sd_team = pm.HalfNormal("sd_team", sigma=0.5)
        team_effect = pm.Normal("team_effect", mu=0, sigma=sd_team, dims="team")
        pos_effect = pm.Normal("pos_effect", mu=0, sigma=1, dims="position")
        beta_age = pm.Normal("beta_age", mu=0, sigma=0.5)
        beta_minutes = pm.Normal("beta_minutes", mu=0, sigma=0.5)
        intercept = pm.Normal("intercept", mu=0.3, sigma=0.5)
        mu = (intercept + team_effect[df["team_idx"].values] + pos_effect[df["pos_idx"].values] + beta_age * age_std.values + beta_minutes * minutes_std.values)
        sigma = pm.HalfNormal("sigma", sigma=0.5)
        goals_obs = pm.Normal("goals_obs", mu=mu, sigma=sigma, observed=observed_goals_per_90, dims="obs")
    return model, df, coords


def sample_hierarchical_model(model, draws=1000, tune=1500, cores=2):
    with model:
        trace = pm.sample(draws=draws, tune=tune, cores=cores, random_seed=42)
    return trace


def analyze_posteriors(trace, df, coords):
    print("=== Team Effect Posterior Summary ===")
    team_summary = az.summary(trace, var_names=["team_effect"])
    print(team_summary.head(10))
    fig, ax = plt.subplots(figsize=(12, 6))
    az.plot_forest(trace, var_names=["team_effect"], ax=ax, combined=True)
    ax.set_title("Team Attack Strength (Posterior 94% HDI)")
    plt.tight_layout()
    plt.savefig("team_effects.png", dpi=150)
    plt.show()
    fig, ax = plt.subplots(figsize=(8, 4))
    az.plot_forest(trace, var_names=["pos_effect"], ax=ax, combined=True)
    ax.set_title("Position Effects on Goals per 90")
    plt.tight_layout()
    plt.savefig("position_effects.png", dpi=150)
    plt.show()
    print("\n=== Convergence Diagnostics ===")
    print(az.summary(trace, var_names=["intercept", "beta_age", "beta_minutes", "sd_team", "sigma"]))
    with trace.model:
        pm.sample_posterior_predictive(trace, extend_inferencedata=True)
    fig, ax = plt.subplots(figsize=(10, 5))
    az.plot_ppc(trace, ax=ax, num_pp_samples=100)
    ax.set_title("Posterior Predictive Check")
    plt.tight_layout()
    plt.savefig("ppc_check.png", dpi=150)
    plt.show()
    return trace


def extract_player_posteriors(trace, df):
    team_effects_mean = trace.posterior["team_effect"].mean(dim=["chain", "draw"]).values
    pos_effects_mean = trace.posterior["pos_effect"].mean(dim=["chain", "draw"]).values
    df["bayesian_performance"] = (team_effects_mean[df["team_idx"].values] + pos_effects_mean[df["pos_idx"].values])
    return df


# =============================================================================
# Step 4: Prepare Survival Analysis Data
# =============================================================================

def prepare_survival_data(df):
    season_length_days = 303
    df["has_injury"] = (df["injury_count"] > 0).astype(int)
    df["minutes_fraction"] = df["total_minutes"] / df["total_minutes"].max()
    df["duration"] = np.where(
        df["has_injury"] == 1,
        np.clip(df["minutes_fraction"] * season_length_days * 0.6, 1, season_length_days),
        season_length_days,
    )
    df["event"] = df["has_injury"]
    return df


# =============================================================================
# Step 5: Kaplan-Meier Curves and Logrank Tests
# =============================================================================

def plot_kaplan_meier(df):
    kmf = KaplanMeierFitter()
    kmf.fit(df["duration"], event_observed=df["event"], label="All Players")
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    kmf.plot_survival_function(ax=axes[0])
    axes[0].set_title("Overall Injury-Free Survival")
    axes[0].set_xlabel("Days into Season")
    axes[0].set_ylabel("Survival Probability")
    axes[0].set_ylim(0, 1)
    print(f"Median survival time (overall): {kmf.median_survival_time_:.0f} days")
    for position in sorted(df["position_group"].unique()):
        mask = df["position_group"] == position
        kmf.fit(df.loc[mask, "duration"], event_observed=df.loc[mask, "event"], label=position)
        kmf.plot_survival_function(ax=axes[1])
    axes[1].set_title("Injury-Free Survival by Position")
    axes[1].set_xlabel("Days into Season")
    axes[1].set_ylabel("Survival Probability")
    axes[1].set_ylim(0, 1)
    axes[1].legend(loc="lower left")
    plt.tight_layout()
    plt.savefig("kaplan_meier.png", dpi=150)
    plt.show()
    for pos in ["Midfield", "Defence", "Goalkeeper"]:
        mask_attack = df["position_group"] == "Attack"
        mask_other = df["position_group"] == pos
        result = logrank_test(
            df.loc[mask_attack, "duration"], df.loc[mask_other, "duration"],
            event_observed_A=df.loc[mask_attack, "event"],
            event_observed_B=df.loc[mask_other, "event"],
        )
        print(f"Logrank (Attack vs {pos}): stat={result.test_statistic:.3f}, p={result.p_value:.4f}")


# =============================================================================
# Step 6: Cox PH and Weibull AFT Models with C-Index
# =============================================================================

def prepare_survival_covariates(df):
    surv_df = df[["duration", "event", "age", "total_minutes", "previous_injuries",
                  "bayesian_performance", "total_yellow_cards"]].copy()
    position_dummies = pd.get_dummies(df["position_group"], prefix="pos", drop_first=True)
    surv_df = pd.concat([surv_df, position_dummies], axis=1)
    surv_df = surv_df.dropna()
    for col in ["age", "total_minutes", "previous_injuries", "bayesian_performance", "total_yellow_cards"]:
        mean_val = surv_df[col].mean()
        std_val = surv_df[col].std()
        if std_val > 0:
            surv_df[col] = (surv_df[col] - mean_val) / std_val
    return surv_df


def fit_cox_model(surv_df):
    cph = CoxPHFitter(penalizer=0.01)
    cph.fit(surv_df, duration_col="duration", event_col="event")
    print("=== Cox Proportional Hazards Model ===")
    cph.print_summary()
    print(f"\nCox C-index: {cph.score(surv_df, scoring_method='concordance_index'):.4f}")
    return cph


def fit_weibull_model(surv_df):
    wft = WeibullAFTFitter(penalizer=0.01)
    wft.fit(surv_df, duration_col="duration", event_col="event")
    print("\n=== Weibull AFT Model ===")
    wft.print_summary()
    print(f"\nWeibull C-index: {wft.score(surv_df, scoring_method='concordance_index'):.4f}")
    return wft


def compare_models(cph, wft, surv_df):
    cox_c = cph.score(surv_df, scoring_method="concordance_index")
    weibull_c = wft.score(surv_df, scoring_method="concordance_index")
    print("\n=== Model Comparison ===")
    print(f"Cox PH C-index:     {cox_c:.4f}")
    print(f"Weibull AFT C-index: {weibull_c:.4f}")
    print(f"Better model: {'Cox PH' if cox_c > weibull_c else 'Weibull AFT'}")
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    cph.plot(ax=axes[0])
    axes[0].set_title(f"Cox PH Hazard Ratios (C={cox_c:.3f})")
    wft.plot(ax=axes[1])
    axes[1].set_title(f"Weibull AFT Coefficients (C={weibull_c:.3f})")
    plt.tight_layout()
    plt.savefig("model_comparison.png", dpi=150)
    plt.show()
    print("\n=== Proportional Hazards Check ===")
    try:
        cph.check_assumptions(surv_df, p_value_threshold=0.05)
    except Exception as e:
        print(f"Note: {e}")
    return cox_c, weibull_c


# =============================================================================
# Step 7: Substitution Recommendation System
# =============================================================================

def predict_injury_risk(cph, wft, player_data):
    """Ensemble prediction from Cox and Weibull models."""
    cox_surv = cph.predict_survival_function(player_data)
    weibull_surv = wft.predict_survival_function(player_data)
    current_day = player_data["duration"].values[0]
    risks = {}
    for horizon_name, horizon_days in [("1_day", 1), ("7_day", 7), ("14_day", 14)]:
        target_day = min(current_day + horizon_days, cox_surv.index.max())
        cox_risk = 1 - cox_surv.loc[target_day].values[0]
        weibull_risk = 1 - weibull_surv.loc[target_day].values[0]
        risks[horizon_name] = (cox_risk + weibull_risk) / 2
    return risks


def classify_risk(risk_7d):
    """Classify injury risk into actionable categories."""
    if risk_7d > 0.30:
        return "HIGH RISK - Recommend Rest"
    elif risk_7d > 0.15:
        return "MONITOR - Consider Rotation"
    else:
        return "LOW RISK - Available"


def generate_squad_report(cph, wft, surv_df, df):
    """Generate risk-performance report for the squad."""
    recommendations = []
    for idx in surv_df.index[:25]:
        player_row = surv_df.loc[[idx]]
        try:
            risks = predict_injury_risk(cph, wft, player_row)
            recommendation = classify_risk(risks["7_day"])
            perf = df.loc[idx, "bayesian_performance"] if idx in df.index else 0
            recommendations.append({
                "player_idx": idx,
                "performance_posterior": perf,
                "risk_1d": risks["1_day"],
                "risk_7d": risks["7_day"],
                "risk_14d": risks["14_day"],
                "recommendation": recommendation,
            })
        except Exception:
            continue
    rec_df = pd.DataFrame(recommendations)
    return rec_df


def plot_risk_performance_scatter(rec_df):
    """Scatter plot: performance vs injury risk."""
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    colors = rec_df["risk_7d"].values
    scatter = axes[0].scatter(
        rec_df["performance_posterior"], rec_df["risk_7d"],
        c=colors, cmap="RdYlGn_r", s=60, edgecolors="black", linewidth=0.5
    )
    axes[0].set_xlabel("Bayesian Performance (Posterior Mean)")
    axes[0].set_ylabel("7-Day Injury Risk")
    axes[0].set_title("Performance vs. Injury Risk")
    axes[0].axhline(y=0.30, color="red", linestyle="--", alpha=0.7, label="High Risk Threshold")
    axes[0].axhline(y=0.15, color="orange", linestyle="--", alpha=0.7, label="Monitor Threshold")
    axes[0].legend()
    plt.colorbar(scatter, ax=axes[0], label="Risk")
    risk_matrix = rec_df[["risk_1d", "risk_7d", "risk_14d"]].copy()
    risk_matrix.columns = ["Next 1 Day", "Next 7 Days", "Next 14 Days"]
    risk_matrix.index = [f"Player {i+1}" for i in range(len(risk_matrix))]
    sns.heatmap(risk_matrix.head(15), annot=True, fmt=".2f", cmap="RdYlGn_r",
                vmin=0, vmax=0.5, ax=axes[1])
    axes[1].set_title("Squad Injury Risk Heatmap")
    plt.tight_layout()
    plt.savefig("squad_risk_performance.png", dpi=150)
    plt.show()


# =============================================================================
# Main Pipeline
# =============================================================================

if __name__ == "__main__":
    PROFILES_PATH = "data/player_profiles.csv"
    PERFORMANCES_PATH = "data/player_performances.csv"
    INJURIES_PATH = "data/player_injury_histories.csv"

    print("Loading data...")
    profiles, performances, injuries = load_data(PROFILES_PATH, PERFORMANCES_PATH, INJURIES_PATH)

    print("Cleaning and merging...")
    df = clean_and_merge(profiles, performances, injuries, target_season="22/23")
    prior_injuries = count_previous_injuries(injuries, target_season="22/23")
    df = df.merge(prior_injuries, on="player_id", how="left")
    df["previous_injuries"] = df["previous_injuries"].fillna(0)
    print(f"Merged dataset: {df.shape[0]} players")

    print("\nBuilding Bayesian hierarchical model...")
    model, df, coords = build_hierarchical_model(df)
    print("Sampling posterior (this may take a few minutes)...")
    trace = sample_hierarchical_model(model, draws=1000, tune=1500, cores=2)
    trace = analyze_posteriors(trace, df, coords)
    df = extract_player_posteriors(trace, df)

    print("\nPreparing survival analysis data...")
    df = prepare_survival_data(df)

    print("\nPlotting Kaplan-Meier curves...")
    plot_kaplan_meier(df)

    print("\nFitting survival models...")
    surv_df = prepare_survival_covariates(df)
    cph = fit_cox_model(surv_df)
    wft = fit_weibull_model(surv_df)
    cox_c, weibull_c = compare_models(cph, wft, surv_df)

    print("\nGenerating squad recommendations...")
    rec_df = generate_squad_report(cph, wft, surv_df, df)
    print(rec_df[["performance_posterior", "risk_7d", "recommendation"]].to_string())
    plot_risk_performance_scatter(rec_df)

Your recommendation system is complete. You have built an end-to-end sports analytics pipeline that goes from raw data through Bayesian hierarchical modeling and survival analysis all the way to actionable squad management decisions. Next up is the Secret Mission, where you will extend this into a full squad rotation optimizer.

Secret mission

Squad Rotation Optimizer

Your substitution recommendation system works for individual players. But real squad management means planning across multiple fixtures, balancing rest and performance for the whole team. Can you build an optimizer that selects the best starting XI while keeping injury risk low?

Clean Up Your Resources

Clean Up Your Resources

This project runs entirely locally, so there are no ongoing cloud costs. Decide whether to keep your resources running, pause them to come back later, or delete them entirely.

Resources you used:

• Conda environment soccer_analytics (contains PyMC, lifelines, pandas, scikit-learn, ArviZ, and all dependencies).
• Project folder on your Desktop (soccer_analytics.py, requirements.txt, environment.yml, data CSVs, and generated plot images).
• VS Code Jupyter extension (if installed for this project).

Keep everything running

No action needed. Choose this if you're still actively building or plan to extend the project with the Secret Mission.

• Your conda environment stays available. Activate it anytime with conda activate soccer_analytics.
• Your project files and data CSVs remain on your Desktop.
• The Jupyter extension stays installed in VS Code for future notebook work.

Pause - I'll come back to this later

Shut down running processes to free up memory. Your files and environment stay intact for when you return.

• Close any running Jupyter notebook kernels by clicking the kernel status in VS Code and selecting Restart or Stop.
• Deactivate your conda environment by running:

conda deactivate

• Your project files and environment remain on disk. Reactivate when you return.

Delete - I don't want to use this again

Remove all project resources and start fresh. This frees disk space used by the conda environment and data files.

• Remove the conda environment by running:

conda env remove --name soccer_analytics

• Delete the project folder (macOS/Linux) by running:

rm -rf ~/Desktop/soccer_analytics

• Delete the project folder (Windows) by running:

Remove-Item -Recurse -Force ~\Desktop\soccer_analytics

• Uninstall the Jupyter extension in VS Code by opening the Extensions panel (Cmd+Shift+X on macOS or Ctrl+Shift+X on Windows), searching for Jupyter, and clicking Uninstall.

That's a Wrap!

That's a Wrap!

Nice work! You have built a professional-grade sports analytics pipeline that combines Bayesian inference with survival analysis to support data-driven squad management decisions.

You've learned how to:

• Built a Bayesian hierarchical model in PyMC that evaluates player performance using team-level random effects, partial pooling, and player covariates.
• Modeled injury risk using both a Cox Proportional Hazards model and a Weibull AFT model with explicit C-index comparison and Kaplan-Meier survival curves.
• Created a substitution recommendation system that combines Bayesian performance posteriors with survival model predictions for actionable squad management.
• Built a squad rotation optimizer that recommends which players to rest based on combined injury risk scores and upcoming fixture difficulty.

Ready to quiz yourself?