Case Study 1: Feature Engineering for NFL Win Prediction

Executive Summary

Predicting NFL game outcomes is one of the most well-studied problems in sports analytics, yet the majority of published models skip the most critical step: rigorous feature engineering. This case study walks through the complete process of transforming raw NFL play-by-play data into a set of predictive features for a binary win-prediction model. We start with over 50,000 plays per season from the nflfastR dataset, engineer 24 candidate features across five categories (efficiency, situational, opponent-adjusted, momentum, and contextual), apply systematic feature selection to reduce the set to 10 high-value features, and demonstrate that thoughtful feature engineering improves out-of-sample accuracy by 3.8 percentage points over a naive baseline that uses season-to-date raw statistics. Along the way, we address temporal leakage, garbage-time filtering, small-sample instability in early-season predictions, and the critical distinction between descriptive and predictive features.

Background

The Prediction Problem

We define the task as predicting the probability that the home team wins an NFL regular-season game, using only information available before kickoff. This is a binary classification task with a roughly 57% base rate for the home team (reflecting home-field advantage). A model that simply predicts "home team wins" for every game would achieve 57% accuracy --- our features must beat this naive baseline to be useful.

Data Source

We use the nflfastR play-by-play dataset, which provides detailed information for every play in every NFL game, including Expected Points Added (EPA), win probability, down, distance, field position, score differential, and hundreds of other fields. Our study covers the 2018-2023 NFL regular seasons (six seasons, approximately 1,600 games).

Why Feature Engineering Matters

The raw play-by-play data contains millions of rows and hundreds of columns, but a game-prediction model needs a single row per game with a compact set of informative features. The transformation from raw plays to game-level features is where most of the modeling value is created or destroyed. A poorly engineered feature set can make even a sophisticated algorithm perform worse than a simple regression; a well-engineered feature set can make a logistic regression outperform a neural network.

Methodology

Step 1: Data Cleaning and Preparation

We begin by loading the play-by-play data and filtering to regular-season games only. Playoff games are excluded because the team composition and motivational context differ significantly.

"""NFL Feature Engineering Pipeline.

Transforms raw play-by-play data into game-level predictive features
for NFL win prediction. Handles temporal leakage, garbage-time filtering,
and opponent adjustments.
"""

import numpy as np
import pandas as pd
from typing import Dict, List, Tuple


def load_and_clean_pbp(seasons: List[int]) -> pd.DataFrame:
    """Load nflfastR play-by-play data and apply basic cleaning.

    Args:
        seasons: List of NFL seasons to load (e.g., [2018, 2019, 2020]).

    Returns:
        Cleaned play-by-play DataFrame with standardized column names.
    """
    frames = []
    for season in seasons:
        url = (
            f"https://github.com/nflverse/nflverse-data/releases/download/"
            f"pbp/play_by_play_{season}.parquet"
        )
        df = pd.read_parquet(url)
        df = df[df["season_type"] == "REG"].copy()
        frames.append(df)

    pbp = pd.concat(frames, ignore_index=True)

    # Filter to actual plays (exclude timeouts, penalties without plays, etc.)
    pbp = pbp[pbp["play_type"].isin(["pass", "run", "punt", "field_goal"])].copy()

    return pbp

Step 2: Garbage-Time Filtering

Garbage-time plays distort team statistics because teams change their playcalling when the game outcome is no longer in doubt. We implement a win-probability-based filter that downweights rather than discards borderline plays.

def compute_garbage_time_weight(row: pd.Series) -> float:
    """Compute a weight for each play based on game competitiveness.

    Plays in competitive game states receive weight 1.0. Plays in
    extreme garbage time (win probability < 5% or > 95%) receive
    weight 0.0. Plays in moderate garbage time are linearly weighted.

    Args:
        row: A row from the play-by-play DataFrame with a 'wp' column
            (win probability for the possession team).

    Returns:
        A float between 0.0 and 1.0 representing the play's weight.
    """
    wp = row["wp"]
    if wp is None or np.isnan(wp):
        return 1.0

    if 0.10 <= wp <= 0.90:
        return 1.0
    elif wp < 0.05 or wp > 0.95:
        return 0.0
    elif wp < 0.10:
        return (wp - 0.05) / 0.05
    else:
        return (0.95 - wp) / 0.05


def apply_garbage_time_filter(pbp: pd.DataFrame) -> pd.DataFrame:
    """Add garbage-time weights to play-by-play data.

    Args:
        pbp: Play-by-play DataFrame with 'wp' column.

    Returns:
        DataFrame with added 'gt_weight' column.
    """
    pbp = pbp.copy()
    pbp["gt_weight"] = pbp.apply(compute_garbage_time_weight, axis=1)
    return pbp

Step 3: Feature Category 1 --- Efficiency Metrics

These are the core performance indicators derived from play-by-play data. We compute EPA/play (the gold-standard efficiency metric), success rate (the consistency metric), and explosive play rate (the big-play metric).

def compute_efficiency_features(
    pbp: pd.DataFrame,
    team: str,
    season: int,
    week: int,
) -> Dict[str, float]:
    """Compute team efficiency features using only prior-week data.

    Args:
        pbp: Play-by-play DataFrame with garbage-time weights.
        team: Team abbreviation (e.g., 'KC').
        season: Season year.
        week: Current week number (features use weeks 1 through week-1).

    Returns:
        Dictionary of efficiency features.
    """
    mask = (
        (pbp["posteam"] == team)
        & (pbp["season"] == season)
        & (pbp["week"] < week)
    )
    team_plays = pbp[mask].copy()

    if len(team_plays) < 20:
        return {
            "off_epa_per_play": np.nan,
            "off_success_rate": np.nan,
            "off_explosive_rate": np.nan,
            "off_early_down_epa": np.nan,
        }

    weights = team_plays["gt_weight"]
    total_weight = weights.sum()

    # Weighted EPA/play
    off_epa = (team_plays["epa"] * weights).sum() / total_weight

    # Success rate: play is successful if EPA > 0
    successes = ((team_plays["epa"] > 0).astype(float) * weights).sum()
    off_success_rate = successes / total_weight

    # Explosive play rate: plays gaining 20+ yards or resulting in a TD
    explosive_mask = (team_plays["yards_gained"] >= 20) | (
        team_plays["touchdown"] == 1
    )
    off_explosive = (explosive_mask.astype(float) * weights).sum() / total_weight

    # Early-down EPA (1st and 2nd down only)
    early_mask = team_plays["down"].isin([1, 2])
    early_plays = team_plays[early_mask]
    early_weights = early_plays["gt_weight"]
    if early_weights.sum() > 0:
        early_epa = (
            (early_plays["epa"] * early_weights).sum() / early_weights.sum()
        )
    else:
        early_epa = np.nan

    return {
        "off_epa_per_play": round(off_epa, 4),
        "off_success_rate": round(off_success_rate, 4),
        "off_explosive_rate": round(off_explosive, 4),
        "off_early_down_epa": round(early_epa, 4),
    }

Step 4: Feature Category 2 --- Situational Features

These features capture game-context factors unrelated to team quality: rest days, travel distance, surface type, and time zone changes.

TEAM_LOCATIONS: Dict[str, Tuple[float, float]] = {
    "KC": (39.0997, -94.5786),
    "BUF": (42.7736, -78.7870),
    "SF": (37.4032, -121.9698),
    "PHI": (39.9008, -75.1675),
    "DAL": (32.7473, -97.0945),
    "MIA": (25.9580, -80.2389),
    "SEA": (47.5952, -122.3316),
    "DEN": (39.7439, -105.0201),
}


def haversine_distance(
    lat1: float, lon1: float, lat2: float, lon2: float
) -> float:
    """Compute the great-circle distance between two points in miles.

    Args:
        lat1: Latitude of point 1 in degrees.
        lon1: Longitude of point 1 in degrees.
        lat2: Latitude of point 2 in degrees.
        lon2: Longitude of point 2 in degrees.

    Returns:
        Distance in miles.
    """
    r = 3958.8  # Earth's radius in miles
    lat1, lon1, lat2, lon2 = map(np.radians, [lat1, lon1, lat2, lon2])
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    a = np.sin(dlat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2
    return 2 * r * np.arcsin(np.sqrt(a))


def compute_situational_features(
    home_team: str,
    away_team: str,
    home_rest_days: int,
    away_rest_days: int,
    is_dome: bool,
) -> Dict[str, float]:
    """Compute situational features for a game.

    Args:
        home_team: Home team abbreviation.
        away_team: Away team abbreviation.
        home_rest_days: Days since the home team's last game.
        away_rest_days: Days since the away team's last game.
        is_dome: Whether the game is played indoors.

    Returns:
        Dictionary of situational features.
    """
    # Rest advantage
    rest_advantage = home_rest_days - away_rest_days
    home_short_rest = 1 if home_rest_days <= 5 else 0
    away_short_rest = 1 if away_rest_days <= 5 else 0

    # Travel distance for away team
    if away_team in TEAM_LOCATIONS and home_team in TEAM_LOCATIONS:
        away_lat, away_lon = TEAM_LOCATIONS[away_team]
        home_lat, home_lon = TEAM_LOCATIONS[home_team]
        travel_distance = haversine_distance(
            away_lat, away_lon, home_lat, home_lon
        )
    else:
        travel_distance = np.nan

    return {
        "rest_advantage": rest_advantage,
        "home_short_rest": home_short_rest,
        "away_short_rest": away_short_rest,
        "away_travel_miles": round(travel_distance, 1),
        "log_away_travel": round(np.log1p(travel_distance), 4),
        "is_dome": int(is_dome),
    }

Step 5: Feature Category 3 --- Opponent-Adjusted Metrics

Raw efficiency statistics are misleading because they depend on opponent quality. A team with +0.10 EPA/play against five top-ten defenses is more impressive than a team with +0.15 EPA/play against five bottom-ten defenses. We implement a simple iterative adjustment.

def compute_opponent_adjusted_ratings(
    game_data: pd.DataFrame,
    n_iterations: int = 15,
    regress_weight: float = 0.2,
) -> pd.DataFrame:
    """Compute opponent-adjusted offensive and defensive ratings.

    Uses an iterative algorithm that alternates between adjusting
    offensive ratings for defensive opponent quality and vice versa.

    Args:
        game_data: DataFrame with columns [team, opponent, week,
            off_epa_per_play, def_epa_per_play].
        n_iterations: Number of adjustment iterations.
        regress_weight: Fraction to regress toward the league mean
            to handle early-season instability.

    Returns:
        DataFrame with adjusted ratings per team per week.
    """
    teams = game_data["team"].unique()
    league_off_mean = game_data["off_epa_per_play"].mean()
    league_def_mean = game_data["def_epa_per_play"].mean()

    # Initialize ratings as raw averages
    off_ratings = game_data.groupby("team")["off_epa_per_play"].mean().to_dict()
    def_ratings = game_data.groupby("team")["def_epa_per_play"].mean().to_dict()

    for _ in range(n_iterations):
        new_off = {}
        new_def = {}

        for team in teams:
            team_games = game_data[game_data["team"] == team]
            n_games = len(team_games)

            # Adjust offensive EPA by subtracting opponent defensive rating
            adj_off_values = []
            for _, game in team_games.iterrows():
                opp = game["opponent"]
                opp_def = def_ratings.get(opp, league_def_mean)
                adj_off_values.append(game["off_epa_per_play"] - opp_def)

            raw_adj_off = np.mean(adj_off_values) if adj_off_values else 0.0

            # Regress toward league mean based on sample size
            regress_factor = regress_weight * (16 / max(n_games, 1))
            new_off[team] = (
                raw_adj_off * (1 - regress_factor)
                + league_off_mean * regress_factor
            )

            # Adjust defensive EPA similarly
            def_games = game_data[game_data["opponent"] == team]
            adj_def_values = []
            for _, game in def_games.iterrows():
                opp = game["team"]
                opp_off = off_ratings.get(opp, league_off_mean)
                adj_def_values.append(game["off_epa_per_play"] - opp_off)

            raw_adj_def = np.mean(adj_def_values) if adj_def_values else 0.0
            n_def_games = len(def_games)
            regress_factor_def = regress_weight * (16 / max(n_def_games, 1))
            new_def[team] = (
                raw_adj_def * (1 - regress_factor_def)
                + league_def_mean * regress_factor_def
            )

        off_ratings = new_off
        def_ratings = new_def

    results = []
    for team in teams:
        results.append({
            "team": team,
            "adj_off_epa": round(off_ratings[team], 4),
            "adj_def_epa": round(def_ratings[team], 4),
            "adj_net_epa": round(
                off_ratings[team] - def_ratings[team], 4
            ),
        })

    return pd.DataFrame(results)

Step 6: Feature Category 4 --- Momentum Features

We capture short-term team trajectory using exponentially weighted metrics and win-streak indicators.

def compute_momentum_features(
    game_log: pd.DataFrame,
    team: str,
    week: int,
    alpha: float = 0.3,
) -> Dict[str, float]:
    """Compute momentum features from a team's game log.

    Args:
        game_log: DataFrame with columns [team, week, points_scored,
            points_allowed, epa_per_play, result].
        team: Team abbreviation.
        week: Current week (features use prior weeks).
        alpha: EWMA decay factor (higher = more weight on recent games).

    Returns:
        Dictionary of momentum features.
    """
    team_games = game_log[
        (game_log["team"] == team) & (game_log["week"] < week)
    ].sort_values("week")

    if len(team_games) < 2:
        return {
            "ewma_epa": np.nan,
            "point_diff_trend": np.nan,
            "win_streak": 0,
            "ats_streak": 0,
        }

    # EWMA of EPA/play
    ewma_epa = team_games["epa_per_play"].ewm(alpha=alpha).mean().iloc[-1]

    # Point differential trend (slope of last 5 games)
    recent = team_games.tail(5)
    if len(recent) >= 3:
        point_diffs = (recent["points_scored"] - recent["points_allowed"]).values
        x = np.arange(len(point_diffs))
        slope = np.polyfit(x, point_diffs, 1)[0]
    else:
        slope = 0.0

    # Current win streak (positive = winning streak, negative = losing)
    streak = 0
    for _, game in team_games.iloc[::-1].iterrows():
        if game["result"] == "W" and streak >= 0:
            streak += 1
        elif game["result"] == "L" and streak <= 0:
            streak -= 1
        else:
            break

    return {
        "ewma_epa": round(ewma_epa, 4),
        "point_diff_trend": round(slope, 4),
        "win_streak": streak,
    }

Step 7: Feature Selection

With 24 candidate features constructed, we apply a systematic selection pipeline to identify the most predictive subset.

from sklearn.ensemble import RandomForestClassifier
from sklearn.feature_selection import RFECV
from sklearn.model_selection import TimeSeriesSplit


def select_features(
    X: pd.DataFrame,
    y: pd.Series,
    min_features: int = 5,
    max_features: int = 15,
) -> List[str]:
    """Select optimal feature subset using recursive feature elimination.

    Args:
        X: Feature matrix (rows = games, columns = features).
        y: Binary target (1 = home win, 0 = away win).
        min_features: Minimum number of features to retain.
        max_features: Maximum number of features to consider.

    Returns:
        List of selected feature names.
    """
    rf = RandomForestClassifier(
        n_estimators=200,
        max_depth=6,
        min_samples_leaf=20,
        random_state=42,
        n_jobs=-1,
    )

    tscv = TimeSeriesSplit(n_splits=5)

    rfecv = RFECV(
        estimator=rf,
        step=1,
        cv=tscv,
        scoring="neg_brier_score",
        min_features_to_select=min_features,
    )

    rfecv.fit(X, y)

    selected = X.columns[rfecv.support_].tolist()
    print(f"Optimal number of features: {rfecv.n_features_}")
    print(f"Selected features: {selected}")

    return selected

Results

Feature Selection Outcome

The RFECV procedure selected 10 features from the original 24 candidates:

The top three features are all opponent-adjusted or efficiency-based, confirming that per-play efficiency metrics are the strongest predictors. Momentum features ranked in the middle, and situational features provided incremental but significant value.

Model Performance Comparison

We trained logistic regression models using walk-forward validation (train on all prior data, predict each week):

The 10-feature selected model outperformed the raw baseline by 3.8 percentage points in accuracy and achieved the best Brier score, demonstrating that feature engineering followed by selection outperforms both naive features and kitchen-sink approaches.

Temporal Leakage Test

To verify our pipeline avoids temporal leakage, we ran a diagnostic: for each game in the test set, we confirmed that every feature value was computed using only data available before that game's kickoff. The maximum temporal gap was zero --- no feature ever referenced a future game.

Key Lessons

1. Opponent adjustment is the single highest-value transformation. Adjusted net EPA ranked first in feature importance, providing information that raw EPA misses entirely.

2. Garbage-time filtering matters. Models trained with unfiltered EPA/play performed 0.8 percentage points worse in accuracy, confirming that garbage-time plays introduce noise.

3. More features do not mean better performance. The 24-feature model slightly outperformed the 4-feature model in-sample but was beaten by the 10-feature selected model out-of-sample, illustrating the bias-variance tradeoff.

4. Momentum features add value beyond efficiency. Even after accounting for season-to-date efficiency (which already captures "how good is this team"), recent-game trajectory features contributed additional predictive power, suggesting that short-term form fluctuations are real and exploitable.

5. Temporal discipline is non-negotiable. A model that accidentally includes future data will report inflated accuracy that collapses in live deployment. The point-in-time feature computation framework prevents this entirely.

Exercises for the Reader

1. Add weather features (temperature, wind speed, precipitation) to the pipeline and test whether they improve totals predictions (over/under) more than spread predictions.

2. Replace the logistic regression with XGBoost and determine whether the optimal feature set changes. Does XGBoost benefit from a larger feature set than logistic regression?

3. Implement a Bayesian alternative to the iterative opponent adjustment that places a prior on team ratings based on preseason market win totals.