Chapter 19: Machine Learning for Soccer

"The goal is to turn data into information, and information into insight." --- Carly Fiorina

Machine learning has fundamentally reshaped how soccer clubs, broadcasters, and governing bodies extract value from the ever-growing volume of match data. From expected goals models that quantify finishing quality to clustering algorithms that discover latent player archetypes, ML techniques now underpin decisions worth hundreds of millions of dollars each transfer window. This chapter provides a rigorous, practitioner-oriented treatment of the ML methods most relevant to soccer analytics. We assume familiarity with the statistical foundations covered in Chapters 3 and 7, and build toward production-ready pipelines that can be integrated into club workflows.

19.1 ML Fundamentals for Soccer Applications

19.1.1 The Machine Learning Landscape

Machine learning algorithms learn patterns from data rather than following explicitly programmed rules. In the soccer domain we encounter three canonical paradigms:

This chapter focuses on supervised and unsupervised methods, which account for the vast majority of applied work in soccer analytics today.

19.1.2 The Supervised Learning Pipeline

Every supervised learning project in soccer follows a common workflow:

1. Problem formulation --- Define the target variable and the decision the model will inform.
2. Data collection --- Aggregate event data (e.g., StatsBomb, Opta), tracking data (e.g., Second Spectrum, SkillCorner), or both.
3. Feature engineering --- Transform raw events into informative predictor variables.
4. Train/validation/test split --- Partition data while respecting temporal ordering (see below).
5. Model selection and tuning --- Compare candidate algorithms; optimize hyper-parameters.
6. Evaluation --- Assess on the held-out test set using task-appropriate metrics.
7. Deployment and monitoring --- Serve predictions; detect model drift.

Callout: Temporal Splitting in Soccer

Standard random train/test splits can leak future information in soccer data. A match played on matchday 30 should never appear in the training set when the target is a matchday 25 event. Always split by date or season:

• Training set: Seasons 2017/18 -- 2020/21
• Validation set: Season 2021/22
• Test set: Season 2022/23

This mirrors the real-world deployment scenario where the model must generalize to unseen future matches.

19.1.3 Data Types in Soccer ML

Soccer data comes in several modalities, each requiring different preprocessing:

Event data contains discrete on-the-ball actions (passes, shots, tackles) with $(x, y)$ coordinates, timestamps, and categorical qualifiers. A single Premier League season produces roughly 500,000--700,000 events.

Tracking data records the $(x, y)$ position of every player and the ball at 25 Hz, yielding approximately 3.3 million frames per match. Tracking data enables features like defensive line height, pressing intensity, and off-ball movement metrics.

Aggregated statistics summarize per-match or per-season totals (goals, assists, progressive passes). These are the coarsest but most widely available data type.

19.1.4 The Bias-Variance Trade-Off

A model's generalization error decomposes as:

$$ \text{Error} = \text{Bias}^2 + \text{Variance} + \text{Irreducible Noise} $$

In soccer contexts:

• High bias (underfitting): A logistic regression predicting match outcomes from possession percentage alone misses the complex interactions that determine results.
• High variance (overfitting): A deep decision tree memorizes specific scorelines from training matches but fails on new fixtures.

The practitioner's task is to navigate this trade-off through model complexity control, regularization, and ensemble methods.

19.1.5 Evaluation Metrics for Soccer ML

Different soccer tasks demand different evaluation criteria:

Callout: Calibration Matters More Than Discrimination

In expected goals modeling, a well-calibrated model --- one whose predicted probabilities match observed frequencies --- is often more valuable than one with a higher AUC. A model that assigns $\hat{p} = 0.15$ to a class of shots should see roughly 15% of those shots result in goals. Always inspect calibration curves alongside discrimination metrics.

19.1.6 Cross-Validation Strategies for Temporal Soccer Data

For soccer data, we recommend time-series cross-validation (also called walk-forward validation):

Fold 1:  Train [S1, S2, S3]  |  Val [S4]
Fold 2:  Train [S1, S2, S3, S4]  |  Val [S5]
Fold 3:  Train [S1, S2, S3, S4, S5]  |  Val [S6]

where $S_i$ denotes season $i$. This ensures:

• No future data leaks into training.
• The training set grows over time, mimicking production conditions.
• Validation performance across folds reveals how the model degrades or improves as more data becomes available.

from sklearn.model_selection import TimeSeriesSplit

tscv = TimeSeriesSplit(n_splits=5)
for train_idx, val_idx in tscv.split(X):
    X_train, X_val = X.iloc[train_idx], X.iloc[val_idx]
    y_train, y_val = y.iloc[train_idx], y.iloc[val_idx]
    # Fit and evaluate model

An alternative approach for soccer data is grouped time-series cross-validation, where the grouping variable is the match date or matchday number. This ensures that all events from a single match are either entirely in the training set or entirely in the validation set, preventing leakage from within-match correlations (e.g., if a team scores one goal, the probability of scoring another in the same match is not independent).

For match-level prediction tasks, another useful strategy is leave-one-season-out cross-validation, where each fold uses a single season as the validation set and all other seasons as training data. This provides a natural measure of how well the model generalizes across seasons and captures any season-to-season variation in playing styles, rule changes, or data quality:

import numpy as np

def leave_one_season_out_cv(X, y, season_labels, model_class, **model_kwargs):
    """Leave-one-season-out cross-validation for soccer data.

    Args:
        X: Feature matrix.
        y: Target vector.
        season_labels: Array of season identifiers for each observation.
        model_class: Sklearn-compatible model class.
        **model_kwargs: Keyword arguments for the model constructor.

    Returns:
        Dictionary mapping season to validation score.
    """
    results = {}
    for season in np.unique(season_labels):
        val_mask = season_labels == season
        train_mask = ~val_mask

        model = model_class(**model_kwargs)
        model.fit(X[train_mask], y[train_mask])
        score = model.score(X[val_mask], y[val_mask])
        results[season] = score

    return results

Callout: Beware of Within-Match Leakage

Even with proper temporal splitting, within-match leakage can inflate performance estimates. If your training set includes some events from a match, and your validation set includes other events from the same match, shared match-level features (score state, team formations, weather) create a subtle form of information leakage. Always ensure that the split boundary falls between matches, not within them. For shot-level xG models, this means assigning all shots from a single match to the same fold.

19.2 Classification Problems in Soccer

Classification is the most common supervised learning task in soccer analytics. We predict discrete outcomes: goal or no goal, home win or away win, successful pass or turnover.

19.2.1 Binary Classification: Goal Prediction

The canonical example is the expected goals (xG) model, which estimates the probability that a shot results in a goal. We treat this as binary classification where $y \in \{0, 1\}$.

Feature set for an xG model:

Logistic Regression Baseline:

$$ P(y=1 \mid \mathbf{x}) = \sigma(\mathbf{w}^T \mathbf{x} + b) = \frac{1}{1 + e^{-(\mathbf{w}^T \mathbf{x} + b)}} $$

where $\sigma$ is the sigmoid function. Logistic regression provides a strong, interpretable baseline. Features like distance_to_goal and angle_to_goal alone yield an AUC of approximately 0.76--0.78 on typical event data.

from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score, brier_score_loss

model = LogisticRegression(C=1.0, max_iter=1000)
model.fit(X_train, y_train)

y_prob = model.predict_proba(X_test)[:, 1]
print(f"AUC: {roc_auc_score(y_test, y_prob):.4f}")
print(f"Brier Score: {brier_score_loss(y_test, y_prob):.4f}")

19.2.2 Multi-Class Classification: Match Outcome Prediction

Predicting match outcomes (home win, draw, away win) is a three-class problem. Common approaches include:

• Multinomial logistic regression (softmax regression)
• Gradient-boosted trees with multi-class log-loss
• Ordinal regression, treating draw as an intermediate outcome

The key challenge is that draws are inherently difficult to predict. In most top-five leagues, draws occur 23--27% of the time, and models struggle to distinguish draws from narrow wins.

Feature engineering for match prediction:

• Rolling averages of xG, xGA over the last $n$ matches (typically $n = 5$ or $n = 10$).
• Elo or Pi-rating differentials.
• Home advantage adjustment.
• Days since last match (fatigue proxy).
• Head-to-head historical record.

19.2.3 Pass Success Prediction

Predicting whether a pass will be completed is valuable for assessing a player's decision-making quality. The model estimates $P(\text{success} \mid \text{pass context})$, and a player who consistently completes passes with low predicted success rates demonstrates above-average passing ability.

Key features include:

• Pass distance and direction.
• Packing rate (number of defenders bypassed).
• Receiver's movement speed and direction.
• Pressure on the passer (nearest opponent distance).
• Pitch zone (defensive third, middle third, attacking third).

Callout: Class Imbalance in Soccer Classification

Many soccer classification tasks exhibit significant class imbalance:

• Shots to goals: Only ~10% of shots result in goals.
• Tackles resulting in fouls: Approximately 25--30%.
• Red card events: Extremely rare (<0.5% of matches for a given player).

Strategies to handle imbalance:

1. Use probability-based metrics (log-loss, Brier score) rather than accuracy.
2. Apply class weights inversely proportional to class frequency.
3. Use SMOTE or other oversampling techniques cautiously --- they can distort calibration.
4. Ensure stratified splits preserve class ratios across folds.

19.2.4 Handling Imbalanced Classes: Goals Are Rare Events

The class imbalance problem deserves deeper treatment because it is so pervasive in soccer analytics. Goals are scored on roughly 10% of shots, meaning that a naive model predicting "no goal" for every shot achieves 90% accuracy while being entirely useless.

Cost-sensitive learning assigns different misclassification costs to different classes. In scikit-learn, this is implemented via the class_weight parameter:

from sklearn.linear_model import LogisticRegression

# Automatically weight classes inversely proportional to frequency
model = LogisticRegression(class_weight="balanced", max_iter=1000)
model.fit(X_train, y_train)

The "balanced" option sets weights to $w_c = \frac{n}{k \cdot n_c}$ where $n$ is the total number of samples, $k$ is the number of classes, and $n_c$ is the number of samples in class $c$. For a dataset with 90% non-goals and 10% goals, the goal class receives approximately 9x the weight of the non-goal class.

SMOTE (Synthetic Minority Over-sampling Technique) generates synthetic examples of the minority class by interpolating between existing minority samples:

from imblearn.over_sampling import SMOTE

smote = SMOTE(random_state=42)
X_resampled, y_resampled = smote.fit_resample(X_train, y_train)

Callout: SMOTE and Calibration

While SMOTE can improve recall for the minority class, it often degrades model calibration. A model trained on SMOTE-resampled data will tend to overestimate the probability of the minority class (goals). If calibrated probabilities are important for your application --- and in soccer analytics, they almost always are --- prefer cost-sensitive learning over resampling, or apply Platt scaling or isotonic regression after training to recalibrate the model's outputs.

Threshold tuning is another approach: rather than using the default 0.5 decision threshold, choose a threshold that optimizes a task-relevant metric (e.g., F1 score, precision at a given recall level). For xG models, the threshold is less relevant because we use the predicted probabilities directly rather than converting them to binary predictions.

Focal loss, introduced by Lin et al. (2017), down-weights easy examples and focuses learning on hard, misclassified examples:

$$ FL(p_t) = -\alpha_t (1 - p_t)^\gamma \log(p_t) $$

where $\gamma > 0$ is a focusing parameter and $\alpha_t$ is a class-balancing weight. With $\gamma = 2$, a well-classified example with $p_t = 0.9$ has its loss reduced by a factor of 100 compared to standard cross-entropy.

19.2.5 Decision Boundaries and Non-Linearity

Linear classifiers assume that the decision boundary is a hyperplane in feature space. In soccer, many relationships are non-linear:

• A shot from 6 yards is far more likely to be scored than one from 7 yards, but the difference between 30 and 31 yards is negligible.
• Header accuracy drops off sharply beyond certain distances but is less sensitive to angle.

Tree-based models and kernel methods capture these non-linearities naturally. We explore these in detail in Section 19.5.

19.2.6 Model Selection: Logistic Regression vs Tree-Based vs Neural Networks

The choice of model architecture depends on the specific task, available data volume, and the relative importance of interpretability versus predictive performance.

Logistic regression remains the gold standard for interpretable soccer models. Its coefficients have clear interpretations: a coefficient of $-0.08$ on distance_to_goal means that each additional meter reduces the log-odds of scoring by 0.08. For tasks where stakeholders need to understand why the model makes specific predictions --- such as presenting xG methodology to coaching staff --- logistic regression is often the best choice.

Tree-based models (decision trees, random forests, gradient boosting) are the workhorses of modern soccer ML. They automatically capture non-linear relationships, feature interactions, and heterogeneous effects across subgroups. A gradient-boosted model can learn that headers from crosses behave differently from headers from corners without requiring the analyst to manually engineer interaction features. The trade-off is reduced interpretability, though tools like SHAP values partially mitigate this.

Neural networks offer the greatest flexibility but require substantially more data and computational resources. In soccer analytics, neural networks are most valuable when processing sequential data (sequences of match events) or spatial data (pitch heatmaps, tracking data frames). For tabular data with fewer than 100,000 training examples --- the typical regime for most soccer tasks --- gradient boosting generally outperforms neural networks.

Callout: The Tabular Data Paradox

Despite the hype around deep learning, for tabular soccer data (the most common format in the industry), gradient-boosted trees consistently match or outperform neural networks in benchmarks. A well-tuned XGBoost or LightGBM model is almost always the right starting point for a new soccer ML project. Reserve neural networks for tasks that genuinely require processing raw sequential or spatial data (e.g., learning from full tracking data frames or event sequences), where the structure of the data naturally suits neural architectures.

19.3 Regression Applications

Regression models predict continuous outcomes. In soccer analytics, regression is used for player valuation, performance rating systems, and continuous expected-value metrics.

19.3.1 Expected Threat (xT) as a Regression Problem

Expected Threat assigns a value to every zone on the pitch representing the probability that possession in that zone leads to a goal within the next $n$ actions. While the original formulation uses a Markov chain, a regression approach can incorporate richer context:

$$ \text{xT}(x, y, \text{context}) = f(x, y, \text{game state}, \text{time}, \text{player positions}) $$

A gradient-boosted regression model can predict the continuous xT value for each action, conditioned on spatial and contextual features.

19.3.2 Player Rating Models

VAEP (Valuing Actions by Estimating Probabilities) and similar frameworks decompose player contributions into offensive and defensive value. The regression targets are:

$$ \Delta P_{\text{score}} = P(\text{score} \mid a_t) - P(\text{score} \mid a_{t-1}) $$

$$ \Delta P_{\text{concede}} = P(\text{concede} \mid a_t) - P(\text{concede} \mid a_{t-1}) $$

where $a_t$ is the action at time $t$. These probability changes are estimated by regression models trained on sequences of match events.

19.3.3 Salary and Transfer Fee Prediction

Predicting player market value is a regression task with high commercial relevance. Features include:

• Age, contract length remaining, current league.
• Performance metrics: goals, assists, xG, xA, progressive carries.
• Market factors: selling club's financial situation, buying club's budget.
• Categorical: position, nationality, agent representation.

Regularization is critical because of multicollinearity among performance features. Ridge regression ($L_2$ penalty) and Lasso ($L_1$ penalty) help:

$$ \mathcal{L}_{\text{Ridge}} = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^{p} w_j^2 $$

$$ \mathcal{L}_{\text{Lasso}} = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^{p} |w_j| $$

Lasso has the added benefit of performing automatic feature selection by shrinking irrelevant coefficients to zero.

from sklearn.linear_model import Ridge, Lasso
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

pipeline = Pipeline([
    ("scaler", StandardScaler()),
    ("ridge", Ridge(alpha=1.0))
])
pipeline.fit(X_train, y_train)

19.3.4 Handling Non-Linear Relationships

Many soccer regression targets exhibit non-linear dependencies. For example, the relationship between a player's age and market value follows an inverted-U shape peaking around age 27. Options for capturing non-linearity include:

1. Polynomial features: Add $x^2$, $x^3$, or interaction terms.
2. Splines: Fit piecewise polynomials with smooth joins.
3. Tree-based models: Gradient-boosted regressors automatically capture non-linearity and interactions.
4. Generalized Additive Models (GAMs): Fit smooth functions of each predictor individually.

19.3.5 Evaluation of Regression Models

Always report multiple metrics, as each tells a different story about model performance.

19.4 Clustering for Player Roles

19.4.1 Why Cluster Players?

Traditional position labels (striker, centre-back, left winger) are increasingly inadequate for describing modern football roles. A "centre-forward" might be a target man, a false nine, or a pressing forward --- three roles with vastly different statistical profiles.

Clustering algorithms discover data-driven role definitions from player performance metrics, enabling:

• Scouting: Find replacement players who match a departing player's statistical profile.
• Tactical analysis: Identify how a team's roles differ from the league norm.
• Player development: Track a young player's role evolution over time.

19.4.2 Feature Selection for Clustering

Feature selection is critical for meaningful clusters. We recommend:

1. Per-90-minute normalization to control for playing time differences.
2. Standardization (z-scores) so that features with larger scales do not dominate the distance metric.
3. Dimensionality reduction (PCA) when using many correlated features.

A typical feature set for outfield player clustering:

19.4.3 K-Means Clustering

K-means partitions $n$ observations into $k$ clusters by minimizing the within-cluster sum of squares:

$$ \min_{\{C_1, \ldots, C_k\}} \sum_{j=1}^{k} \sum_{\mathbf{x}_i \in C_j} \|\mathbf{x}_i - \boldsymbol{\mu}_j\|^2 $$

where $\boldsymbol{\mu}_j$ is the centroid of cluster $C_j$.

Choosing $k$: Use the elbow method (plot inertia vs. $k$) and the silhouette score. For outfield player roles in top-five leagues, $k = 7$ to $k = 12$ typically yields interpretable clusters.

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

inertias = []
silhouettes = []
K_range = range(3, 16)

for k in K_range:
    km = KMeans(n_clusters=k, n_init=10, random_state=42)
    labels = km.fit_predict(X_scaled)
    inertias.append(km.inertia_)
    silhouettes.append(silhouette_score(X_scaled, labels))

19.4.4 Hierarchical Clustering

Agglomerative hierarchical clustering builds a tree (dendrogram) of nested clusters by iteratively merging the closest pair of clusters. Linkage criteria include:

• Ward's method: Minimizes the total within-cluster variance. Tends to produce compact, equally sized clusters.
• Complete linkage: Uses the maximum pairwise distance between clusters.
• Average linkage: Uses the mean pairwise distance.

Hierarchical clustering is particularly useful in soccer because the dendrogram reveals the granularity of role distinctions. For example, at a coarse level "attackers" and "defenders" separate; at a finer level, "ball-playing centre-backs" split from "traditional centre-backs."

19.4.5 Gaussian Mixture Models

Gaussian Mixture Models (GMMs) provide a probabilistic alternative to K-means. Each cluster is modeled as a multivariate Gaussian:

$$ p(\mathbf{x}) = \sum_{j=1}^{k} \pi_j \, \mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}_j, \boldsymbol{\Sigma}_j) $$

GMMs offer soft assignments: a versatile midfielder might have 60% probability of belonging to the "box-to-box" cluster and 40% to the "advanced playmaker" cluster. This reflects the reality that player roles exist on a continuum.

19.4.6 Discovering Tactical Patterns with Unsupervised Learning

Beyond player role clustering, unsupervised learning can discover tactical patterns at the team level:

Formation detection. K-means or GMM clustering applied to the average positions of all 10 outfield players across a match (or within specific game states) can automatically detect the formation a team is playing. By clustering the 20-dimensional vector of average $(x, y)$ positions, the algorithm discovers that some matches cluster around a 4-3-3 shape, others around a 4-4-2 shape, and so on. This is particularly valuable for detecting in-game formation changes, where a team shifts from a 4-3-3 to a 3-5-2 after a substitution.

Pressing pattern classification. Clustering the spatial configuration of the defending team at the moment of ball recovery reveals distinct pressing archetypes: high-press recoveries (ball recovered in the attacking third), mid-block recoveries (middle third), and low-block recoveries (defensive third). Further sub-clustering within each zone reveals whether the team uses man-oriented or space-oriented pressing.

Set-piece grouping. Clustering the spatial configuration of players at corner kicks or free kicks reveals different set-piece routines. Teams typically have 3--5 distinct corner kick setups, and clustering can automatically identify which routine was used for each corner, enabling analysis of the success rate of each approach.

19.4.7 Validating Clusters

Cluster validation combines quantitative metrics with domain expertise:

Callout: The Importance of Domain Validation

A clustering solution with a high silhouette score but uninterpretable clusters is useless in practice. Always visualize cluster centroids using radar charts and present results to domain experts (coaches, scouts) for validation. The best clustering is the one that tells scouts something they did not already know while remaining believable.

19.4.8 Dimensionality Reduction for Visualization

High-dimensional player profiles need to be projected into 2D for visualization. Common techniques:

• PCA (Principal Component Analysis): Linear projection preserving maximum variance. Fast and deterministic. PCA is particularly useful in soccer analytics because the principal components often have interpretable meanings: the first component frequently separates attacking from defensive players, while the second separates central from wide players.
• t-SNE: Non-linear embedding that preserves local neighborhood structure. Good for visualization but not stable across runs. The perplexity parameter (typically 20--50) controls the balance between local and global structure. For player clustering visualizations, perplexity values around 30 often produce the most interpretable plots.
• UMAP (Uniform Manifold Approximation and Projection): Similar to t-SNE but faster and better at preserving global structure. UMAP tends to produce tighter, more separated clusters than t-SNE, making it particularly effective for visualizing player role clusters. It is also deterministic when a random seed is set, unlike t-SNE.

from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

pca = PCA(n_components=2)
X_2d = pca.fit_transform(X_scaled)

plt.scatter(X_2d[:, 0], X_2d[:, 1], c=cluster_labels, cmap="tab10", alpha=0.7)
plt.xlabel(f"PC1 ({pca.explained_variance_ratio_[0]:.1%} variance)")
plt.ylabel(f"PC2 ({pca.explained_variance_ratio_[1]:.1%} variance)")
plt.title("Player Roles via PCA")
plt.colorbar(label="Cluster")
plt.show()

Callout: PCA Component Interpretation

When using PCA for player analysis, examine the loadings (component coefficients) to understand what each principal component represents. In a typical outfield player analysis, the first component often loads heavily on attacking metrics (xG, shots, touches in the box) with negative loadings on defensive metrics (tackles, interceptions). The second component frequently separates central players (high pass completion, progressive passes) from wide players (high dribbles, crosses). Naming your components based on their loadings --- e.g., "Attacking vs. Defensive Orientation" and "Central vs. Wide Profile" --- makes the visualization immediately interpretable to non-technical stakeholders.

19.5 Ensemble Methods and Model Stacking

19.5.1 Why Ensembles?

Ensemble methods combine multiple models to achieve better predictive performance than any single model. The theoretical justification comes from Condorcet's jury theorem and the bias-variance decomposition:

• Bagging (Bootstrap Aggregating) reduces variance by averaging predictions from models trained on bootstrapped samples.
• Boosting reduces bias by sequentially training models that focus on the errors of their predecessors.
• Stacking combines heterogeneous models using a meta-learner.

19.5.2 Random Forests

A random forest is an ensemble of decision trees trained via bagging with random feature subsampling at each split.

Key hyperparameters:

from sklearn.ensemble import RandomForestClassifier

rf = RandomForestClassifier(
    n_estimators=500,
    max_depth=12,
    min_samples_leaf=20,
    max_features="sqrt",
    random_state=42,
    n_jobs=-1
)
rf.fit(X_train, y_train)

Feature importance from random forests provides a useful (if imperfect) ranking of predictor relevance. For xG models, distance_to_goal and angle_to_goal consistently rank highest.

19.5.3 Gradient Boosting

Gradient boosting builds an additive model by sequentially fitting trees to the negative gradient of the loss function. At iteration $m$:

$$ F_m(\mathbf{x}) = F_{m-1}(\mathbf{x}) + \eta \cdot h_m(\mathbf{x}) $$

where $\eta$ is the learning rate and $h_m$ is the tree fitted to the pseudo-residuals of $F_{m-1}$.

Modern implementations:

from sklearn.ensemble import GradientBoostingClassifier

gbc = GradientBoostingClassifier(
    n_estimators=300,
    learning_rate=0.05,
    max_depth=5,
    subsample=0.8,
    random_state=42
)
gbc.fit(X_train, y_train)

19.5.4 Hyperparameter Tuning

For gradient boosting, the most impactful hyperparameters are:

1. Learning rate ($\eta$): Lower values require more trees but often generalize better. Start with 0.05--0.1.
2. Number of trees (n_estimators): Use early stopping on validation loss to determine the optimal number.
3. Tree depth (max_depth): 3--8 for boosting (shallower than in random forests).
4. Subsampling rate (subsample): 0.7--0.9 introduces stochasticity that reduces overfitting.
5. Regularization (reg_alpha, reg_lambda): L1/L2 penalties on leaf weights.

Use RandomizedSearchCV or Bayesian optimization (e.g., Optuna) for efficient hyperparameter search:

from sklearn.model_selection import RandomizedSearchCV
from scipy.stats import uniform, randint

param_dist = {
    "n_estimators": randint(100, 800),
    "learning_rate": uniform(0.01, 0.19),
    "max_depth": randint(3, 9),
    "subsample": uniform(0.6, 0.4),
    "min_samples_leaf": randint(10, 50)
}

search = RandomizedSearchCV(
    GradientBoostingClassifier(random_state=42),
    param_distributions=param_dist,
    n_iter=50,
    scoring="neg_log_loss",
    cv=TimeSeriesSplit(n_splits=4),
    random_state=42,
    n_jobs=-1
)
search.fit(X_train, y_train)
print(f"Best log-loss: {-search.best_score_:.4f}")
print(f"Best params: {search.best_params_}")

19.5.5 Model Stacking

Stacking (stacked generalization) trains a meta-learner on the out-of-fold predictions of multiple base models:

Level 0 (base models): - Logistic regression - Random forest - Gradient boosting - K-nearest neighbors

Level 1 (meta-learner): - Logistic regression (or another simple model) trained on the stacked predictions

The key is to use out-of-fold predictions for the Level 0 features to avoid information leakage:

from sklearn.ensemble import StackingClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from sklearn.neighbors import KNeighborsClassifier

estimators = [
    ("lr", LogisticRegression(max_iter=1000)),
    ("rf", RandomForestClassifier(n_estimators=200, random_state=42)),
    ("gb", GradientBoostingClassifier(n_estimators=200, random_state=42)),
    ("knn", KNeighborsClassifier(n_neighbors=15))
]

stacking_clf = StackingClassifier(
    estimators=estimators,
    final_estimator=LogisticRegression(),
    cv=5,
    stack_method="predict_proba"
)
stacking_clf.fit(X_train, y_train)

Callout: When to Use Stacking

Stacking provides marginal gains (often 0.5--1.5% improvement in AUC) at the cost of significant complexity. In soccer analytics:

• Use stacking for high-stakes models (e.g., xG models used in broadcast graphics or recruitment decisions).
• Avoid stacking for exploratory analysis or when interpretability is paramount.
• Always benchmark against a well-tuned single gradient boosting model first.

19.5.6 Comparing Model Performance

When comparing models, use paired statistical tests or confidence intervals rather than single-point estimates:

from sklearn.model_selection import cross_val_score

models = {
    "Logistic Regression": LogisticRegression(max_iter=1000),
    "Random Forest": RandomForestClassifier(n_estimators=300, random_state=42),
    "Gradient Boosting": GradientBoostingClassifier(n_estimators=300, random_state=42),
    "Stacking": stacking_clf
}

for name, model in models.items():
    scores = cross_val_score(
        model, X_train, y_train,
        cv=TimeSeriesSplit(n_splits=5),
        scoring="neg_log_loss"
    )
    print(f"{name}: Log-loss = {-scores.mean():.4f} (+/- {scores.std():.4f})")

19.6 Feature Engineering and Selection

19.6.1 The Art of Feature Engineering

Feature engineering is often the difference between a mediocre model and a state-of-the-art one. In soccer, domain knowledge is the primary driver of good features.

Categories of engineered features:

1. Spatial features: Distance to goal, angle to goal, pitch zone indicators, distance to nearest defender.
2. Temporal features: Time remaining in match, time since last event, possession duration.
3. Sequential features: Previous $n$ actions (action type, location, outcome), possession chain length.
4. Contextual features: Score differential, home/away, match importance (league stage, knockout round).
5. Aggregated features: Rolling averages (xG over last 5 matches), season totals, career statistics.

19.6.2 Feature Engineering from Soccer Event Data

Event data provides a rich source of engineerable features, but extracting maximum value requires understanding both the data structure and the game itself:

Action sequence features capture the buildup to an event. For an xG model, the sequence of actions leading to a shot --- was it preceded by a cross, a through ball, a dribble, or a set piece? --- significantly affects the goal probability. Encoding the last 3--5 actions as a feature vector (action type, location, success/failure) adds substantial predictive power. A common approach is to create binary indicators for the most informative preceding actions: preceded_by_cross, preceded_by_through_ball, preceded_by_cutback, preceded_by_individual_play.

Spatial derivative features go beyond raw $(x, y)$ coordinates to capture spatial relationships:

• Distance to nearest defender: Computed from tracking data or estimated from event data qualifiers.
• Angle bisector features: For shots, the angle between the direction of approach and the line to goal captures whether the shooter is running toward goal (favorable) or across the face of goal (less favorable).
• Zone transition indicators: Binary flags indicating whether the action moved the ball from one pitch zone to another (e.g., middle third to attacking third).

Possession-level aggregations summarize the characteristics of the entire possession chain up to the current action:

• Total distance the ball has traveled during the possession.
• Number of passes in the possession.
• Maximum and minimum $x$-coordinates reached (how far forward and backward the possession has gone).
• Duration of the possession.
• Number of progressive actions (passes or carries that move the ball significantly toward goal).

def engineer_shot_features(shot_event, preceding_events, tracking_frame=None):
    """Engineer features for a single shot event.

    Args:
        shot_event: Dictionary with shot attributes.
        preceding_events: List of events in the possession chain before the shot.
        tracking_frame: Optional tracking data at the moment of the shot.

    Returns:
        Dictionary of engineered features.
    """
    features = {}

    # Spatial features
    goal_center = (105.0, 34.0)  # Assuming standard coordinates
    dx = goal_center[0] - shot_event["x"]
    dy = goal_center[1] - shot_event["y"]
    features["distance_to_goal"] = np.sqrt(dx**2 + dy**2)
    features["angle_to_goal"] = np.abs(np.arctan2(dy, dx))

    # Sequence features
    if preceding_events:
        last_action = preceding_events[-1]
        features["preceded_by_cross"] = int(last_action["type"] == "cross")
        features["preceded_by_through_ball"] = int(last_action["type"] == "through_ball")
        features["possession_length"] = len(preceding_events)
        features["possession_duration"] = (
            shot_event["timestamp"] - preceding_events[0]["timestamp"]
        )
    else:
        features["preceded_by_cross"] = 0
        features["preceded_by_through_ball"] = 0
        features["possession_length"] = 0
        features["possession_duration"] = 0

    # Contextual features
    features["score_diff"] = shot_event.get("score_diff", 0)
    features["minute"] = shot_event.get("minute", 45)
    features["is_home"] = int(shot_event.get("is_home", True))

    return features

19.6.3 Encoding Categorical Variables

Soccer data contains many categorical features (body part, action type, team name, formation). Encoding strategies:

from sklearn.preprocessing import OneHotEncoder
from sklearn.compose import ColumnTransformer

categorical_features = ["body_part", "shot_type", "prev_action"]
numerical_features = ["distance_to_goal", "angle_to_goal", "gk_distance"]

preprocessor = ColumnTransformer([
    ("num", StandardScaler(), numerical_features),
    ("cat", OneHotEncoder(drop="first", sparse_output=False), categorical_features)
])

19.6.4 Feature Selection Methods

With dozens of potential features, selection prevents overfitting and improves interpretability:

Filter methods rank features by statistical criteria independent of the model:

• Mutual information: Captures non-linear dependencies between feature and target.
• Chi-squared test: For categorical features vs. categorical target.
• Correlation analysis: Identify and remove highly correlated feature pairs ($|r| > 0.85$).

Wrapper methods evaluate feature subsets by training models:

• Recursive Feature Elimination (RFE): Iteratively removes the least important feature.
• Forward/backward selection: Greedily adds or removes features.

Embedded methods perform selection during model training:

• Lasso ($L_1$ regularization): Drives irrelevant coefficients to zero.
• Tree-based importance: Features that appear frequently in splits and reduce impurity the most are ranked highest.

from sklearn.feature_selection import mutual_info_classif, SelectKBest

selector = SelectKBest(mutual_info_classif, k=15)
X_selected = selector.fit_transform(X_train, y_train)

# Get selected feature names
selected_mask = selector.get_support()
selected_features = X_train.columns[selected_mask].tolist()
print(f"Selected features: {selected_features}")

19.6.5 Feature Engineering Pipeline

A production-ready feature engineering pipeline should be:

1. Reproducible: All transformations are defined in code, not manual spreadsheets.
2. Fitted on training data only: Scalers, encoders, and imputers learn parameters from training data and transform test data accordingly.
3. Versioned: Feature definitions are tracked alongside model versions.

from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer

full_pipeline = Pipeline([
    ("imputer", SimpleImputer(strategy="median")),
    ("preprocessor", preprocessor),
    ("selector", SelectKBest(mutual_info_classif, k=20)),
    ("classifier", GradientBoostingClassifier(random_state=42))
])

full_pipeline.fit(X_train, y_train)
y_pred = full_pipeline.predict_proba(X_test)[:, 1]

19.6.6 Interaction Features

Some of the most powerful features in soccer ML are interactions:

• Distance x Body Part: Headers from close range are scored at a much higher rate than headers from distance.
• Angle x Defender Count: A wide angle with no defenders is very different from a wide angle with a crowded box.
• Score Differential x Time Remaining: Trailing by one goal with 5 minutes left changes team behavior dramatically.

Tree-based models discover interactions automatically, but explicit interaction features can help linear models and improve interpretability.

19.6.7 Dealing with Missing Data

Soccer datasets frequently have missing values:

• Tracking data coverage may be partial.
• Historical data lacks modern event qualifiers.
• Injury and fitness data may be proprietary and incomplete.

Strategies:

19.7 Model Interpretation: SHAP Values and Feature Importance

19.7.1 Why Interpretability Matters in Soccer

In soccer analytics, model interpretability is not a luxury --- it is a requirement. Coaches and sporting directors will not trust a model they cannot understand. A black-box model that predicts "sign this player" without explaining why is unlikely to influence real decisions. Interpretability also helps analysts debug models, detect data quality issues, and build domain knowledge.

19.7.2 Permutation Feature Importance

Permutation importance measures how much a model's performance degrades when a feature is randomly shuffled. For each feature $j$:

1. Compute the baseline model score $s$ on the validation set.
2. Randomly shuffle feature $j$ across all validation samples.
3. Recompute the model score $s_j^{\text{shuffled}}$.
4. The importance is $I_j = s - s_j^{\text{shuffled}}$.

Repeat this process multiple times and take the mean to reduce variance. Permutation importance is model-agnostic and avoids the biases of tree-based impurity importance (which favors high-cardinality features).

from sklearn.inspection import permutation_importance

result = permutation_importance(
    model, X_val, y_val,
    n_repeats=10,
    random_state=42,
    scoring="neg_log_loss"
)

for i in result.importances_mean.argsort()[::-1]:
    if result.importances_mean[i] > 0.001:
        print(f"{feature_names[i]}: "
              f"{result.importances_mean[i]:.4f} +/- {result.importances_std[i]:.4f}")

19.7.3 SHAP Values

SHAP (SHapley Additive exPlanations) values provide a unified framework for interpreting any machine learning model. Based on cooperative game theory, SHAP assigns each feature a contribution to the prediction for a specific input:

$$ \hat{f}(\mathbf{x}) = \phi_0 + \sum_{j=1}^{p} \phi_j(\mathbf{x}) $$

where $\phi_0$ is the expected model output (the mean prediction) and $\phi_j$ is the SHAP value for feature $j$ --- the marginal contribution of feature $j$ to the prediction for this specific input.

For a shot-level xG prediction, SHAP values answer questions like: "This shot has an xG of 0.35. The average xG is 0.10. The short distance to goal contributes +0.18, the favorable angle contributes +0.05, but the header body part contributes -0.03, and the two defenders in the cone contribute -0.05."

import shap

explainer = shap.TreeExplainer(gbc)
shap_values = explainer.shap_values(X_test)

# Summary plot: global feature importance
shap.summary_plot(shap_values, X_test, feature_names=feature_names)

# Waterfall plot: single prediction explanation
shap.waterfall_plot(shap.Explanation(
    values=shap_values[0],
    base_values=explainer.expected_value,
    data=X_test.iloc[0],
    feature_names=feature_names
))

Callout: SHAP for Scouting Reports

SHAP values can be used to generate automated scouting reports that explain why a model rates a particular player highly. For a transfer target, the report might say: "This player's predicted goals contribution next season is 14.2. The primary drivers are: xG per 90 (+3.5 goals), age 24 (+1.8 goals, pre-peak player), league adjustment from Eredivisie to Premier League (-2.1 goals), and high progressive carrying rate (+1.0 goals)." This format bridges the gap between model output and decision-making by providing actionable explanations that scouts and sporting directors can evaluate using their domain expertise.

19.7.4 Partial Dependence Plots

Partial dependence plots (PDPs) show the marginal effect of a single feature on the model's predictions, averaging over the values of all other features:

$$ \hat{f}_j(x_j) = \frac{1}{n} \sum_{i=1}^{n} \hat{f}(x_j, \mathbf{x}_{i,-j}) $$

For an xG model, a PDP of distance_to_goal reveals the shape of the relationship between distance and goal probability: a steep decline from 0 to 15 meters, followed by a gradual flattening beyond 20 meters. This visualization is intuitive even for non-technical audiences and can be used in presentations to coaching staff.

19.8 Avoiding Overfitting in Small-Sample Soccer Datasets

19.8.1 The Small-Sample Problem

Soccer datasets are often smaller than practitioners realize. While a season contains approximately 380 matches in a major league, the effective sample size for many tasks is much smaller:

• Team-level match prediction: Only 380 matches per season, with each team appearing 38 times. A model with 20+ features risks overfitting severely.
• Player-level season statistics: Only 500--600 players per league meet a reasonable minutes threshold (e.g., 900 minutes). Many features per player create a high-dimensional problem.
• Rare events: Only 20--30 red cards per league per season, fewer than 10 penalty misses, and perhaps 5 own goals. Building models for these events requires extreme care.

19.8.2 Regularization Techniques

Regularization is the primary defense against overfitting:

• L1 regularization (Lasso): Drives irrelevant feature weights to exactly zero, performing automatic feature selection. Useful when you suspect that many features are irrelevant.
• L2 regularization (Ridge): Shrinks all weights toward zero without eliminating any. Useful when features are correlated and all potentially relevant.
• Elastic Net: Combines L1 and L2 penalties, offering a balance between feature selection and coefficient shrinkage.
• Early stopping: For iterative models (gradient boosting, neural networks), monitor validation loss and stop training when it begins to increase.
• Dropout (neural networks): Randomly sets a fraction of neuron activations to zero during training, preventing co-adaptation.

19.8.3 Cross-Validation as an Overfitting Detector

If cross-validation performance is substantially worse than training performance, the model is overfitting. A useful diagnostic is the generalization gap:

$$ \text{Gap} = \text{Training Score} - \text{CV Score} $$

A gap exceeding 5--10% of the training score warrants investigation. Common remedies include reducing model complexity, adding regularization, removing features, or acquiring more training data.

Callout: The "n > p" Rule of Thumb

A rough guideline for soccer datasets: the number of training observations ($n$) should exceed the number of features ($p$) by at least a factor of 10--20. For a match prediction model with 380 matches, this suggests using no more than 19--38 features. Exceeding this ratio without strong regularization is a recipe for overfitting. When you have more features than this guideline suggests, use dimensionality reduction (PCA) or feature selection before training.

19.9 Common Pitfalls in Soccer ML Applications

19.9.1 Survivorship Bias

Many soccer datasets only include players who achieved a certain level of success --- those who played in top leagues, received professional contracts, or were involved in transfers. Models trained on these biased samples may not generalize to the full population of potential players. For example, a transfer success model trained only on completed transfers ignores all the players who were scouted but not signed --- the "near misses" that contain valuable negative examples.

19.9.2 Target Leakage

Target leakage occurs when information that would not be available at prediction time is included in the training features. Common examples in soccer:

• Including a player's season-end statistics when predicting mid-season performance.
• Using the final match score as a feature in a model predicting in-game events.
• Including xG values computed by a different model as features for your own xG model (circular dependency).

19.9.3 Simpson's Paradox in Soccer Statistics

Aggregate statistics can reverse direction when disaggregated by a confounding variable. For example, a player may have a higher pass completion rate than another player overall, but a lower rate in every individual match. This paradox arises when the two players attempt passes with systematically different difficulty distributions. Always condition on relevant confounders (pitch zone, opposition quality, game state) before drawing conclusions from aggregate statistics.

19.9.4 The Curse of Metrics Proliferation

Modern soccer data provides hundreds of metrics per player per match. The temptation to include all of them in a model leads to overfitting, multicollinearity, and interpretability loss. A disciplined feature selection process --- guided by domain knowledge, statistical criteria, and cross-validation performance --- is essential.

Callout: The Five-Feature Rule

Before building a complex model with dozens of features, ask: "Can I build a model with five features that captures 80% of the predictive power?" In many soccer tasks, the answer is yes. For xG: distance, angle, body part, shot type, and whether it followed a cross. For match prediction: home/away, Elo rating difference, recent form, days rest, and key player availability. Start with these parsimonious models and only add complexity when the cross-validated improvement justifies it.

19.10 MLOps: Deploying and Monitoring Soccer ML Models

19.10.1 From Notebook to Production

The majority of soccer ML models never leave a Jupyter notebook. Moving to production requires:

1. Code refactoring: Modularize feature engineering, model training, and prediction into separate functions or classes.
2. Testing: Unit tests for feature calculations, integration tests for the full pipeline.
3. Serialization: Save trained models using joblib or pickle.
4. API serving: Wrap the model in a REST API (Flask, FastAPI) for real-time or batch predictions.

import joblib

# Save the trained pipeline
joblib.dump(full_pipeline, "xg_model_v2.1.joblib")

# Load and predict
loaded_pipeline = joblib.load("xg_model_v2.1.joblib")
predictions = loaded_pipeline.predict_proba(new_data)[:, 1]

19.10.2 Model Versioning

Maintain a model registry that tracks:

19.10.3 Monitoring for Model Drift

Soccer evolves. Rule changes (e.g., VAR introduction, five-substitute rule), tactical trends, and new data sources can cause model drift --- a degradation in model performance over time.

Types of drift:

• Data drift (covariate shift): The distribution of input features changes. Example: a new league season features more shots from outside the box due to tactical trends.
• Concept drift: The relationship between features and the target changes. Example: VAR overturning goals changes the effective conversion rate for certain shot types.
• Label drift: The distribution of the target variable changes.

Monitoring strategies:

1. Track live performance metrics weekly/monthly against a baseline.
2. Population Stability Index (PSI) for detecting distributional shifts in features:

$$ \text{PSI} = \sum_{i=1}^{B} (p_i^{\text{new}} - p_i^{\text{ref}}) \ln\left(\frac{p_i^{\text{new}}}{p_i^{\text{ref}}}\right) $$

where $B$ is the number of bins and $p_i$ is the proportion in bin $i$. A PSI above 0.2 typically indicates significant drift.

3. Calibration monitoring: Plot predicted vs. observed probabilities for rolling windows of matches.
4. Automated alerts: Trigger retraining when performance drops below a threshold.

19.10.4 Retraining Strategies

19.10.5 Ethical Considerations and Fairness

Machine learning models in soccer raise ethical questions:

• Player privacy: Tracking data reveals granular physical performance metrics. Models built on this data must comply with data protection regulations (e.g., GDPR).
• Algorithmic bias: If training data overrepresents certain leagues, models may systematically undervalue players from underrepresented leagues.
• Transparency: When ML models influence multi-million-dollar transfer decisions, stakeholders deserve explanations. Use SHAP values or LIME for local interpretability.

import shap

explainer = shap.TreeExplainer(gbc)
shap_values = explainer.shap_values(X_test)

shap.summary_plot(shap_values, X_test, feature_names=feature_names)

Callout: Algorithmic Bias in Scouting Models

Scouting models trained predominantly on data from Europe's top five leagues may systematically undervalue players from leagues with less comprehensive data coverage (South American leagues, African leagues, Asian leagues). This creates a feedback loop: undervalued players are less likely to be scouted, which means less data is collected on them, which further reduces their model-predicted value. Responsible ML practitioners should audit their models for league-based bias, use league adjustment factors, and supplement quantitative models with qualitative scouting assessments for underrepresented markets.

19.10.6 A Production Architecture

A complete soccer ML system typically includes:

+------------------+     +-------------------+     +------------------+
|  Data Ingestion  | --> | Feature Pipeline  | --> |  Model Training  |
|  (APIs, files)   |     | (Spark/Pandas)    |     |  (sklearn, XGB)  |
+------------------+     +-------------------+     +------------------+
                                                          |
                                                          v
+------------------+     +-------------------+     +------------------+
|  Monitoring &    | <-- |  Prediction API   | <-- |  Model Registry  |
|  Alerting        |     |  (FastAPI)        |     |  (MLflow)        |
+------------------+     +-------------------+     +------------------+

Callout: Start Simple, Iterate

The most common mistake in soccer ML projects is over-engineering the initial solution. Start with:

1. A logistic regression baseline with 5--10 well-chosen features.
2. A simple train/test split by season.
3. A joblib serialized model loaded by a Python script.

Only add complexity (gradient boosting, stacking, real-time APIs, drift monitoring) when you have evidence that the simpler approach is insufficient.

Summary

This chapter has covered the full lifecycle of machine learning in soccer analytics, from problem formulation through deployment and monitoring. The key themes are:

1. Respect the data structure: Soccer data is temporal, spatial, and hierarchical. Standard ML recipes must be adapted accordingly.
2. Start with strong baselines: Logistic regression and simple decision trees are surprisingly competitive for many soccer tasks.
3. Feature engineering is king: Domain-informed features consistently outperform algorithmic complexity.
4. Ensemble methods push the frontier: Gradient boosting and stacking provide the best predictive performance for structured soccer data.
5. Clustering reveals hidden structure: Data-driven player roles complement traditional positional labels.
6. Model interpretation is essential: SHAP values and permutation importance bridge the gap between model outputs and actionable insights for coaches, scouts, and sporting directors.
7. Handle imbalanced classes carefully: Goals, injuries, and other rare events require specialized techniques to model effectively.
8. Avoid common pitfalls: Survivorship bias, target leakage, small-sample overfitting, and metrics proliferation are pervasive in soccer ML.
9. Production readiness matters: A model that cannot be deployed, monitored, and maintained delivers no value.

In the next chapter, we extend these foundations to deep learning approaches that can process raw tracking data and learn representations end-to-end.