Learning ObjectivesBy the end of this chapter, you will be able to:

  1. Build comprehensive offensive evaluation frameworks combining multiple metrics
  2. Calculate and interpret advanced efficiency metrics (DVOA, EPA/play, success rate)
  3. Analyze drive efficiency and scoring probability
  4. Adjust offensive metrics for pace, schedule strength, and era
  5. Decompose offensive performance by unit (line, skill positions, coordinator)
  6. Compare offensive systems and their effectiveness
  7. Project future offensive performance using historical data
  8. Create professional-grade offensive scouting reports

Introduction

Throughout Part II, we've examined individual components of offensive performance: passing efficiency, rushing effectiveness, EPA generation, success rates, and play-calling tendencies. While each of these analyses provides valuable insights, evaluating an offense comprehensively requires synthesizing multiple metrics into a holistic framework.

This chapter represents the culmination of everything we've learned in Part II. Think of the previous chapters as examining individual instruments in an orchestra—each important in its own right, but incomplete without the others. Here, we learn to hear the full symphony: how passing and rushing efficiency interact, how consistency and explosiveness balance, how context shapes performance, and how to build evaluation frameworks that capture the multi-dimensional nature of offensive football.

Modern NFL front offices don't evaluate offenses based solely on yards per game or points scored. They build sophisticated models combining efficiency metrics, situational performance, opponent adjustments, and trend analysis. A team might rank 5th in yards but 18th in efficiency. Another might score many points due to favorable field position from turnovers, masking underlying offensive struggles. A third might appear mediocre overall but excel in the critical situations that determine wins.

This chapter equips you to conduct professional-grade offensive analysis. We'll learn how to build composite metrics that balance multiple dimensions, adjust for pace and schedule strength, decompose total performance into unit-level contributions, identify coordinator effects, compare offenses across eras, and create comprehensive scouting reports. By the end, you'll be able to evaluate any NFL offense with the same rigor and nuance used by analytics departments across the league.

What is Comprehensive Offensive Efficiency?

Comprehensive offensive efficiency is a multi-dimensional evaluation framework that: - Combines multiple performance metrics (EPA, success rate, explosiveness) - Adjusts for external factors (opponent strength, pace, game situation) - Attributes performance to specific units (QB, line, receivers, coordinator) - Enables fair comparisons across teams, seasons, and eras - Provides actionable insights for team evaluation and decision-making

The Need for Comprehensive Analysis

Limitations of Single Metrics

No single metric tells the complete story of offensive performance. Each traditional and advanced metric captures an important dimension but misses critical context that separates truly elite offenses from those that merely appear successful.

Total Yards: The most visible statistic—total offensive yards—fundamentally misses efficiency. A team might gain 400 yards by running 80 plays (5.0 yards per play), while another gains 350 yards on 55 plays (6.4 yards per play). The second offense is far more efficient, but yards alone won't reveal this. Yardage also ignores field position: 50 yards gained from your own 20 to your 30 is worth less than 50 yards from the opponent's 40 to the end zone.

Points Per Game: Scoring statistics seem definitive—the goal is to score points, after all. But points per game conflates offensive performance with factors beyond the offense's control. Defensive touchdowns inflate point totals without offensive contribution. Turnovers by the defense create short fields, making scoring easier. Teams that play faster get more possessions, inflating per-game statistics. A team averaging 27 points per game on 10 drives looks worse than one averaging 24 points on 8 drives, despite identical 2.7 points per drive efficiency.

EPA/Play: Expected Points Added per play is our best single efficiency metric, capturing field position value and down-and-distance context. But even EPA has blindspots. It doesn't distinguish between an offense that consistently gains 4-5 yards (successful but not explosive) from one that mixes 1-yard gains with frequent 20-yard explosions (same average EPA, very different profile). EPA treats a +0.4 EPA play the same whether it's a critical third-down conversion or a meaningless second down gain, despite different strategic value.

Success Rate: The percentage of plays with positive EPA measures consistency—how reliably an offense moves the chains and avoids negative plays. But consistency alone doesn't win games. An offense with 50% success rate but no explosive plays will struggle to score. Another with 45% success rate but frequent touchdowns may be more dangerous. Success rate also treats all positive plays equally: barely converting third-and-10 counts the same as a 40-yard touchdown.

Common Pitfall: Single-Metric Rankings

Many analysts rank offenses solely by EPA/play or yards per game, then express surprise when rankings don't predict playoff success. Elite offenses typically excel across multiple dimensions simultaneously: efficiency, consistency, explosiveness, and situational performance. A team ranking 3rd in EPA/play but 22nd in explosive play rate has a fundamentally different profile—and different strategic implications—than one ranking 3rd in both metrics.

The Multi-Metric Approach

A comprehensive evaluation framework combines complementary metrics that together capture the full picture of offensive performance. Modern analytics departments typically track 10-15 core metrics spanning several key dimensions:

1. Efficiency: How well the offense generates expected points relative to opportunities. EPA per play remains our foundation, but we can enhance it by examining efficiency in different contexts: early downs versus late downs, neutral game script versus leading/trailing, inside versus outside the red zone. An offense might be efficient overall but inefficient in critical situations, which affects win probability more than raw averages suggest.

2. Consistency: How reliably the offense executes successfully and avoids setbacks. Success rate (percentage of plays with positive EPA) provides the baseline, but we also want to examine negative play rate (plays losing substantial EPA), three-and-out percentage, and drive sustainability. Consistent offenses control the clock, wear down defenses, and protect their own defense by limiting opponent possessions.

3. Explosiveness: Big-play capability that changes games and creates easy scoring opportunities. We measure this through explosive play rate (percentage of plays gaining 20+ passing yards or 10+ rushing yards), EPA from explosive plays specifically, and the ability to generate chunk gains that put the defense on their heels. Some offenses efficiently march down the field; others strike quickly. Both can succeed, but they require different supporting casts and defensive matchups.

4. Sustainability: Drive success and scoring efficiency that convert possessions into points. Points per drive, scoring drive percentage, and touchdown rate in the red zone measure how often the offense finishes drives. An offense might move the ball well (high success rate) but stall in scoring position (low red zone efficiency), scoring many field goals instead of touchdowns. Over a season, these differences compound into significant point differentials.

5. Context: Adjustments for external factors that shape performance but don't reflect offensive quality. Opponent defensive strength, home/away splits, weather conditions, pace of play, and game script all affect raw statistics. A team playing the league's easiest schedule might produce impressive stats while being merely average. Strength-of-schedule and pace adjustments enable fair comparisons across teams and seasons.

Key Insight: Dimensional Trade-offs

The best offenses don't necessarily lead every category. Instead, they find optimal combinations suited to their personnel. A team with an elite quarterback but weak offensive line might emphasize quick passes and explosiveness over consistency. Another with a dominant line but average quarterback might prioritize consistency and sustainable drives. Understanding these trade-offs—why teams excel in certain dimensions and struggle in others—reveals strategic choices and personnel constraints that shape offensive identity.

This multi-dimensional approach mirrors how NFL coaches actually evaluate offenses. When preparing for an opponent, defensive coordinators don't just look at yards per game. They examine third-down efficiency (do we need to stop them on early downs or can we rely on getting them in third-and-long?), explosive play rate (do we need to prevent the big play or focus on consistent pressure?), red zone tendencies (are they more dangerous passing or running near the goal line?), and situational performance (how do they respond when trailing in the fourth quarter?). Our analytics should provide the same comprehensive picture.

Comprehensive Efficiency Metrics

DVOA (Defense-adjusted Value Over Average)

Before building our own composite metrics, let's understand the industry standard. Defense-adjusted Value Over Average (DVOA), pioneered by Football Outsiders, has been the gold standard for comprehensive offensive evaluation since the early 2000s. Unlike raw statistics, DVOA compares each play to the league-average performance in the same situation, then adjusts for opponent strength.

The fundamental insight behind DVOA is that gaining 4 yards on 3rd-and-3 is far more valuable than gaining 4 yards on 3rd-and-15, even though both produce identical yardage. DVOA uses historical success rates for every down-distance-field position combination to establish baseline expectations, then measures how much each play exceeded or fell short of that baseline.

The mathematical representation is:

$$ \text{DVOA} = \frac{\text{Success}_{\text{actual}} - \text{Success}_{\text{expected}}}{\text{Success}_{\text{expected}}} $$

Where success is measured in terms of yards gained relative to expected yards for each down-and-distance situation. A DVOA of +15% means the offense performed 15% better than league average in comparable situations. The opponent adjustment further refines this by accounting for defensive quality—beating an elite defense is worth more than beating a poor one.

DVOA vs EPA: Complementary Approaches

DVOA and EPA measure similar concepts through different methodologies. DVOA uses success rates based on yards gained and first downs, while EPA uses expected point values. DVOA adjusts for opponent strength more directly; EPA captures field position value more precisely. Both correlate strongly with winning (~0.7-0.8 correlation with point differential), and elite offenses typically rank highly in both metrics. We'll use EPA as our foundation because the data is readily available in nflfastR, but the conceptual approach mirrors DVOA's multi-dimensional philosophy.

EPA-Based Composite Metrics

Building on DVOA's principles, we can construct our own comprehensive metric using EPA components. The key is combining efficiency, consistency, and explosiveness into a single rating that balances multiple dimensions of offensive performance.

Our composite offensive rating follows this formula:

$$ \text{Offensive Rating} = w_1 \cdot \text{EPA/play} + w_2 \cdot \text{Success Rate} + w_3 \cdot \text{Explosiveness} $$

Where weights ($w_1, w_2, w_3$) are optimized to predict future scoring or wins. But how do we choose these weights? Several approaches exist:

Equal Weighting: The simplest approach assigns equal importance to each component (33% each). This makes no assumptions about relative importance but may underweight EPA efficiency, which typically correlates most strongly with winning.

Empirical Optimization: We can use regression analysis to find weights that best predict point differential or win rate. Historical analysis suggests EPA/play deserves the highest weight (~40-50%), with success rate and explosiveness splitting the remainder. This approach is data-driven but requires careful validation to avoid overfitting.

Strategic Weighting: Different contexts might require different weights. When evaluating playoff performance, explosive play capability might deserve higher weight since playoff defenses typically limit sustained drives. When evaluating offensive line quality specifically, we might emphasize negative play avoidance. The weights should match your analytical goal.

Standardization: Before combining metrics, we standardize each to z-scores (standard deviations from the mean). This ensures that metrics with different scales contribute proportionally to the final rating. Without standardization, a metric ranging 0-100 would dominate one ranging 0-1, regardless of actual importance.

Best Practice: Validate Your Weights

After selecting weights, validate them using out-of-sample testing: use one season's weights to predict the next season's performance. If your composite rating from Season N correlates strongly with Season N+1 success, your weights capture sustainable performance rather than random variation. Correlation of 0.5+ with next season's EPA/play indicates good predictive validity.

Throughout this chapter, we'll use these weights: 40% EPA/play, 25% success rate, 20% explosiveness, and 15% negative play avoidance. These weights emerged from empirical testing showing the best predictive accuracy for future offensive performance while maintaining face validity (EPA/play as the primary driver makes intuitive sense).

Setting Up Our Environment

#| label: setup-r
#| message: false
#| warning: false

# Load required packages
library(tidyverse)
library(nflfastR)
library(nflplotR)
library(gt)
library(gtExtras)
library(ggridges)
library(patchwork)

# Set theme
theme_set(theme_minimal())

# Load multiple seasons for comprehensive analysis
pbp <- load_pbp(2020:2023)

cat("✓ Loaded", nrow(pbp), "plays from 2020-2023 seasons\n")

#| label: setup-py
#| message: false
#| warning: false

import pandas as pd
import numpy as np
import nfl_data_py as nfl
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from sklearn.preprocessing import StandardScaler

# Set display options
pd.set_option('display.max_columns', None)
pd.set_option('display.width', None)

# Load multiple seasons
pbp = nfl.import_pbp_data(list(range(2020, 2024)))

print(f"✓ Loaded {len(pbp):,} plays from 2020-2023 seasons")

Calculating Comprehensive Offensive Ratings

Basic Efficiency Metrics

Let's start by calculating the fundamental efficiency metrics for each team:

#| label: basic-efficiency-r
#| message: false
#| warning: false

# Calculate comprehensive offensive metrics
offensive_efficiency <- pbp %>%
  filter(season == 2023, !is.na(epa), !is.na(posteam)) %>%
  filter(play_type %in% c("pass", "run")) %>%
  group_by(posteam) %>%
  summarise(
    plays = n(),
    epa_per_play = mean(epa),
    success_rate = mean(epa > 0),
    explosiveness = mean(epa[abs(epa) > 1], na.rm = TRUE),

    # Separate pass and run efficiency
    pass_epa = mean(epa[play_type == "pass"], na.rm = TRUE),
    run_epa = mean(epa[play_type == "run"], na.rm = TRUE),
    pass_success = mean(epa[play_type == "pass"] > 0, na.rm = TRUE),
    run_success = mean(epa[play_type == "run"] > 0, na.rm = TRUE),

    # Big play rates
    explosive_pass_rate = mean(yards_gained[play_type == "pass"] >= 20, na.rm = TRUE),
    explosive_run_rate = mean(yards_gained[play_type == "run"] >= 10, na.rm = TRUE),

    # Negative plays
    negative_play_rate = mean(epa < -0.5),

    .groups = "drop"
  ) %>%
  arrange(desc(epa_per_play))

# Display top 10 offenses
offensive_efficiency %>%
  slice_head(n = 10) %>%
  gt() %>%
  cols_label(
    posteam = "Team",
    plays = "Plays",
    epa_per_play = "EPA/Play",
    success_rate = "Success%",
    explosiveness = "Explosive EPA",
    pass_epa = "Pass EPA",
    run_epa = "Run EPA",
    pass_success = "Pass SR",
    run_success = "Run SR",
    explosive_pass_rate = "Exp Pass%",
    explosive_run_rate = "Exp Run%",
    negative_play_rate = "Neg Play%"
  ) %>%
  fmt_number(
    columns = c(epa_per_play, pass_epa, run_epa, explosiveness),
    decimals = 3
  ) %>%
  fmt_percent(
    columns = c(success_rate, pass_success, run_success,
                explosive_pass_rate, explosive_run_rate, negative_play_rate),
    decimals = 1
  ) %>%
  fmt_number(columns = plays, decimals = 0, use_seps = TRUE) %>%
  data_color(
    columns = c(epa_per_play, success_rate),
    palette = "RdYlGn"
  ) %>%
  tab_header(
    title = "Comprehensive Offensive Efficiency Metrics",
    subtitle = "2023 NFL Season - Top 10 Offenses by EPA/Play"
  ) %>%
  tab_source_note("Data: nflfastR")

#| label: basic-efficiency-py
#| message: false
#| warning: false

# Filter data
pbp_2023 = pbp[(pbp['season'] == 2023) &
               (pbp['epa'].notna()) &
               (pbp['posteam'].notna()) &
               (pbp['play_type'].isin(['pass', 'run']))]

# Calculate metrics for each team
def calculate_efficiency_metrics(df):
    pass_plays = df[df['play_type'] == 'pass']
    run_plays = df[df['play_type'] == 'run']

    return pd.Series({
        'plays': len(df),
        'epa_per_play': df['epa'].mean(),
        'success_rate': (df['epa'] > 0).mean(),
        'explosiveness': df[df['epa'].abs() > 1]['epa'].mean() if len(df[df['epa'].abs() > 1]) > 0 else 0,
        'pass_epa': pass_plays['epa'].mean() if len(pass_plays) > 0 else 0,
        'run_epa': run_plays['epa'].mean() if len(run_plays) > 0 else 0,
        'pass_success': (pass_plays['epa'] > 0).mean() if len(pass_plays) > 0 else 0,
        'run_success': (run_plays['epa'] > 0).mean() if len(run_plays) > 0 else 0,
        'explosive_pass_rate': (pass_plays['yards_gained'] >= 20).mean() if len(pass_plays) > 0 else 0,
        'explosive_run_rate': (run_plays['yards_gained'] >= 10).mean() if len(run_plays) > 0 else 0,
        'negative_play_rate': (df['epa'] < -0.5).mean()
    })

offensive_efficiency = (pbp_2023
    .groupby('posteam')
    .apply(calculate_efficiency_metrics)
    .reset_index()
    .sort_values('epa_per_play', ascending=False)
)

# Display top 10
print("\nTop 10 Offenses by EPA/Play (2023):")
print(offensive_efficiency.head(10).to_string(index=False))

Composite Offensive Rating

Now let's create a composite rating that combines multiple dimensions:

#| label: composite-rating-r
#| message: false
#| warning: false

# Create composite offensive rating
# Standardize metrics to z-scores for fair weighting
offensive_rating <- offensive_efficiency %>%
  mutate(
    # Standardize key metrics
    epa_z = scale(epa_per_play)[,1],
    success_z = scale(success_rate)[,1],
    explosive_z = scale(explosiveness)[,1],
    negative_z = scale(-negative_play_rate)[,1], # Negative sign because lower is better

    # Weighted composite (EPA/play weighted most heavily)
    composite_rating = (
      0.40 * epa_z +           # 40% weight to EPA efficiency
      0.25 * success_z +       # 25% weight to consistency
      0.20 * explosive_z +     # 20% weight to explosiveness
      0.15 * negative_z        # 15% weight to avoiding negative plays
    ),

    # Convert to 0-100 scale for interpretability
    rating_100 = ((composite_rating - min(composite_rating)) /
                  (max(composite_rating) - min(composite_rating))) * 100
  ) %>%
  arrange(desc(rating_100)) %>%
  select(posteam, plays, epa_per_play, success_rate, explosiveness,
         negative_play_rate, composite_rating, rating_100)

# Display rankings with team logos
offensive_rating %>%
  gt() %>%
  cols_label(
    posteam = "Team",
    plays = "Plays",
    epa_per_play = "EPA/Play",
    success_rate = "Success%",
    explosiveness = "Explosive EPA",
    negative_play_rate = "Neg Play%",
    composite_rating = "Z-Score",
    rating_100 = "Rating (0-100)"
  ) %>%
  fmt_number(
    columns = c(epa_per_play, explosiveness, composite_rating),
    decimals = 3
  ) %>%
  fmt_percent(
    columns = c(success_rate, negative_play_rate),
    decimals = 1
  ) %>%
  fmt_number(
    columns = c(plays, rating_100),
    decimals = 0
  ) %>%
  data_color(
    columns = rating_100,
    palette = "RdYlGn",
    domain = c(0, 100)
  ) %>%
  tab_header(
    title = "Composite Offensive Rating",
    subtitle = "Combines EPA efficiency, consistency, explosiveness, and discipline"
  ) %>%
  tab_source_note("Rating: 40% EPA/play, 25% Success Rate, 20% Explosiveness, 15% Negative Play Avoidance")

#| label: composite-rating-py
#| message: false
#| warning: false

from sklearn.preprocessing import StandardScaler

# Standardize metrics
scaler = StandardScaler()
metrics_to_scale = ['epa_per_play', 'success_rate', 'explosiveness', 'negative_play_rate']

# Create copy and calculate z-scores
offensive_rating = offensive_efficiency.copy()
z_scores = scaler.fit_transform(offensive_rating[metrics_to_scale])

offensive_rating['epa_z'] = z_scores[:, 0]
offensive_rating['success_z'] = z_scores[:, 1]
offensive_rating['explosive_z'] = z_scores[:, 2]
offensive_rating['negative_z'] = -z_scores[:, 3]  # Negative because lower is better

# Calculate composite rating
offensive_rating['composite_rating'] = (
    0.40 * offensive_rating['epa_z'] +
    0.25 * offensive_rating['success_z'] +
    0.20 * offensive_rating['explosive_z'] +
    0.15 * offensive_rating['negative_z']
)

# Convert to 0-100 scale
min_rating = offensive_rating['composite_rating'].min()
max_rating = offensive_rating['composite_rating'].max()
offensive_rating['rating_100'] = (
    ((offensive_rating['composite_rating'] - min_rating) /
     (max_rating - min_rating)) * 100
)

# Sort and display
offensive_rating = offensive_rating.sort_values('rating_100', ascending=False)

print("\nComposite Offensive Ratings (2023):")
print(offensive_rating[['posteam', 'epa_per_play', 'success_rate',
                        'explosiveness', 'rating_100']].head(15).to_string(index=False))

Visualizing Multi-Dimensional Performance

#| label: fig-efficiency-scatter-r
#| fig-cap: "Offensive efficiency vs consistency (2023)"
#| fig-width: 12
#| fig-height: 8
#| message: false
#| warning: false

# Create scatter plot with team logos
offensive_efficiency %>%
  ggplot(aes(x = success_rate, y = epa_per_play)) +
  geom_hline(yintercept = mean(offensive_efficiency$epa_per_play),
             linetype = "dashed", color = "gray50", alpha = 0.7) +
  geom_vline(xintercept = mean(offensive_efficiency$success_rate),
             linetype = "dashed", color = "gray50", alpha = 0.7) +
  geom_nfl_logos(aes(team_abbr = posteam), width = 0.05, alpha = 0.8) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_y_continuous(labels = scales::number_format(accuracy = 0.01)) +
  labs(
    title = "Offensive Efficiency: EPA vs Success Rate",
    subtitle = "2023 NFL Season | Lines indicate league average",
    x = "Success Rate (% of plays with positive EPA)",
    y = "EPA per Play",
    caption = "Data: nflfastR | Visualization: nflplotR"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
    plot.subtitle = element_text(size = 11, hjust = 0.5),
    axis.title = element_text(size = 12),
    panel.grid.minor = element_blank()
  ) +
  # Add quadrant labels
  annotate("text", x = 0.35, y = 0.15, label = "Low Efficiency\nLow Consistency",
           color = "gray30", size = 3.5, alpha = 0.6) +
  annotate("text", x = 0.51, y = 0.15, label = "High Consistency\nLow Efficiency",
           color = "gray30", size = 3.5, alpha = 0.6) +
  annotate("text", x = 0.35, y = 0.30, label = "High Efficiency\nLow Consistency",
           color = "gray30", size = 3.5, alpha = 0.6) +
  annotate("text", x = 0.51, y = 0.30, label = "Elite Offense",
           color = "darkgreen", size = 4, fontface = "bold", alpha = 0.7)

#| label: fig-efficiency-scatter-py
#| fig-cap: "Offensive efficiency vs consistency - Python (2023)"
#| fig-width: 12
#| fig-height: 8
#| message: false
#| warning: false

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

# Plot points
scatter = ax.scatter(offensive_efficiency['success_rate'],
                     offensive_efficiency['epa_per_play'],
                     s=200, alpha=0.6, c=offensive_efficiency['rating_100'],
                     cmap='RdYlGn', edgecolors='black', linewidth=1)

# Add team labels
for idx, row in offensive_efficiency.iterrows():
    ax.annotate(row['posteam'],
                (row['success_rate'], row['epa_per_play']),
                fontsize=8, ha='center', va='center', fontweight='bold')

# Add reference lines
ax.axhline(y=offensive_efficiency['epa_per_play'].mean(),
           color='gray', linestyle='--', alpha=0.5, label='League Avg EPA')
ax.axvline(x=offensive_efficiency['success_rate'].mean(),
           color='gray', linestyle='--', alpha=0.5, label='League Avg Success%')

# Labels and formatting
ax.set_xlabel('Success Rate (% of plays with positive EPA)', fontsize=12)
ax.set_ylabel('EPA per Play', fontsize=12)
ax.set_title('Offensive Efficiency: EPA vs Success Rate\n2023 NFL Season',
             fontsize=14, fontweight='bold', pad=20)

# Add colorbar
cbar = plt.colorbar(scatter, ax=ax)
cbar.set_label('Composite Rating', rotation=270, labelpad=20)

# Format axes
ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, p: f'{x:.0%}'))
ax.grid(True, alpha=0.3)
ax.legend(loc='lower right')

plt.tight_layout()
plt.show()

Drive Efficiency and Scoring Probability

Analyzing Drive Success

Drive efficiency measures how effectively an offense converts possessions into points:

#| label: drive-efficiency-r
#| message: false
#| warning: false

# Calculate drive-level metrics
drive_efficiency <- pbp %>%
  filter(season == 2023, !is.na(fixed_drive)) %>%
  group_by(posteam, fixed_drive, game_id) %>%
  summarise(
    drive_plays = n(),
    drive_yards = sum(yards_gained, na.rm = TRUE),
    drive_epa = sum(epa, na.rm = TRUE),
    drive_result = last(series_result),
    scored_td = any(touchdown == 1),
    scored_fg = any(field_goal_result == "made"),
    scored_any = scored_td | scored_fg,
    points_scored = sum(ifelse(touchdown == 1, 7, 0) +
                       ifelse(field_goal_result == "made", 3, 0)),
    .groups = "drop"
  ) %>%
  group_by(posteam) %>%
  summarise(
    drives = n(),
    plays_per_drive = mean(drive_plays),
    yards_per_drive = mean(drive_yards),
    epa_per_drive = mean(drive_epa),
    scoring_drive_rate = mean(scored_any),
    td_drive_rate = mean(scored_td),
    fg_drive_rate = mean(scored_fg),
    points_per_drive = mean(points_scored),
    three_and_out_rate = mean(drive_plays <= 3 & !scored_any),
    .groups = "drop"
  ) %>%
  arrange(desc(points_per_drive))

# Display drive efficiency table
drive_efficiency %>%
  gt() %>%
  cols_label(
    posteam = "Team",
    drives = "Drives",
    plays_per_drive = "Plays/Drive",
    yards_per_drive = "Yards/Drive",
    epa_per_drive = "EPA/Drive",
    scoring_drive_rate = "Score%",
    td_drive_rate = "TD%",
    fg_drive_rate = "FG%",
    points_per_drive = "Pts/Drive",
    three_and_out_rate = "3&Out%"
  ) %>%
  fmt_number(
    columns = c(plays_per_drive, yards_per_drive, epa_per_drive, points_per_drive),
    decimals = 2
  ) %>%
  fmt_percent(
    columns = c(scoring_drive_rate, td_drive_rate, fg_drive_rate, three_and_out_rate),
    decimals = 1
  ) %>%
  fmt_number(columns = drives, decimals = 0) %>%
  data_color(
    columns = c(points_per_drive, scoring_drive_rate),
    palette = "RdYlGn"
  ) %>%
  tab_header(
    title = "Drive Efficiency Metrics",
    subtitle = "2023 NFL Season"
  ) %>%
  tab_source_note("Data: nflfastR")

#| label: drive-efficiency-py
#| message: false
#| warning: false

# Calculate drive-level metrics
pbp_2023_drives = pbp[(pbp['season'] == 2023) & (pbp['fixed_drive'].notna())]

# Group by drive
drive_summary = []
for (team, drive, game), drive_data in pbp_2023_drives.groupby(['posteam', 'fixed_drive', 'game_id']):
    drive_summary.append({
        'posteam': team,
        'drive_plays': len(drive_data),
        'drive_yards': drive_data['yards_gained'].sum(),
        'drive_epa': drive_data['epa'].sum(),
        'scored_td': (drive_data['touchdown'] == 1).any(),
        'scored_fg': (drive_data['field_goal_result'] == 'made').any(),
        'scored_any': (drive_data['touchdown'] == 1).any() or (drive_data['field_goal_result'] == 'made').any(),
        'points_scored': (drive_data[drive_data['touchdown'] == 1].shape[0] * 7 +
                         drive_data[drive_data['field_goal_result'] == 'made'].shape[0] * 3)
    })

drive_df = pd.DataFrame(drive_summary)

# Calculate team-level metrics
drive_efficiency = (drive_df
    .groupby('posteam')
    .agg(
        drives=('drive_plays', 'count'),
        plays_per_drive=('drive_plays', 'mean'),
        yards_per_drive=('drive_yards', 'mean'),
        epa_per_drive=('drive_epa', 'mean'),
        scoring_drive_rate=('scored_any', 'mean'),
        td_drive_rate=('scored_td', 'mean'),
        fg_drive_rate=('scored_fg', 'mean'),
        points_per_drive=('points_scored', 'mean')
    )
    .reset_index()
    .assign(
        three_and_out_rate=lambda x: drive_df.groupby('posteam').apply(
            lambda g: ((g['drive_plays'] <= 3) & (~g['scored_any'])).mean()
        ).values
    )
    .sort_values('points_per_drive', ascending=False)
)

print("\nDrive Efficiency Metrics (2023):")
print(drive_efficiency.head(10).to_string(index=False))

Scoring Probability by Field Position

#| label: fig-scoring-prob-r
#| fig-cap: "Scoring probability by field position (2023)"
#| fig-width: 12
#| fig-height: 7
#| message: false
#| warning: false

# Calculate scoring probability by field position
scoring_by_position <- pbp %>%
  filter(season == 2023, !is.na(yardline_100), !is.na(fixed_drive)) %>%
  filter(down == 1, ydstogo == 10) %>%  # First and 10 for consistency
  mutate(
    field_zone = case_when(
      yardline_100 >= 90 ~ "Own 10",
      yardline_100 >= 80 ~ "Own 20",
      yardline_100 >= 70 ~ "Own 30",
      yardline_100 >= 60 ~ "Own 40",
      yardline_100 >= 50 ~ "Midfield",
      yardline_100 >= 40 ~ "Opp 40",
      yardline_100 >= 30 ~ "Opp 30",
      yardline_100 >= 20 ~ "Opp 20",
      yardline_100 >= 10 ~ "Opp 10",
      TRUE ~ "Red Zone"
    ),
    field_zone = factor(field_zone, levels = c("Own 10", "Own 20", "Own 30",
                                               "Own 40", "Midfield", "Opp 40",
                                               "Opp 30", "Opp 20", "Opp 10",
                                               "Red Zone"))
  ) %>%
  group_by(fixed_drive, game_id, field_zone) %>%
  slice(1) %>%  # First occurrence of each drive in each zone
  ungroup() %>%
  left_join(
    pbp %>%
      filter(season == 2023) %>%
      group_by(fixed_drive, game_id) %>%
      summarise(
        scored_td = any(touchdown == 1, na.rm = TRUE),
        scored_fg = any(field_goal_result == "made", na.rm = TRUE),
        scored_any = scored_td | scored_fg,
        .groups = "drop"
      ),
    by = c("fixed_drive", "game_id")
  ) %>%
  group_by(field_zone) %>%
  summarise(
    drives = n(),
    td_prob = mean(scored_td, na.rm = TRUE),
    fg_prob = mean(scored_fg, na.rm = TRUE),
    any_score_prob = mean(scored_any, na.rm = TRUE),
    avg_epa = mean(epa, na.rm = TRUE),
    .groups = "drop"
  )

# Create visualization
scoring_by_position %>%
  pivot_longer(cols = c(td_prob, fg_prob, any_score_prob),
               names_to = "score_type", values_to = "probability") %>%
  mutate(score_type = factor(score_type,
                             levels = c("any_score_prob", "fg_prob", "td_prob"),
                             labels = c("Any Score", "Field Goal", "Touchdown"))) %>%
  ggplot(aes(x = field_zone, y = probability, fill = score_type)) +
  geom_col(position = "dodge", alpha = 0.8) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),
                    expand = c(0, 0), limits = c(0, 0.8)) +
  scale_fill_manual(values = c("Any Score" = "#4CAF50",
                               "Field Goal" = "#2196F3",
                               "Touchdown" = "#FF5722")) +
  labs(
    title = "Scoring Probability by Starting Field Position",
    subtitle = "1st & 10 situations, 2023 NFL Season",
    x = "Starting Field Position",
    y = "Probability of Scoring",
    fill = "Score Type",
    caption = "Data: nflfastR"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "top",
    panel.grid.major.x = element_blank()
  )

📊 Visualization Output

The code above generates a visualization. To see the output, run this code in your R or Python environment. The resulting plot will help illustrate the concepts discussed in this section.

#| label: fig-scoring-prob-py
#| fig-cap: "Scoring probability by field position - Python (2023)"
#| fig-width: 12
#| fig-height: 7
#| message: false
#| warning: false

# Filter for first and 10
first_and_10 = pbp_2023_drives[
    (pbp_2023_drives['down'] == 1) &
    (pbp_2023_drives['ydstogo'] == 10) &
    (pbp_2023_drives['yardline_100'].notna())
].copy()

# Create field zones
def categorize_field_position(yards):
    if yards >= 90: return "Own 10"
    elif yards >= 80: return "Own 20"
    elif yards >= 70: return "Own 30"
    elif yards >= 60: return "Own 40"
    elif yards >= 50: return "Midfield"
    elif yards >= 40: return "Opp 40"
    elif yards >= 30: return "Opp 30"
    elif yards >= 20: return "Opp 20"
    elif yards >= 10: return "Opp 10"
    else: return "Red Zone"

first_and_10['field_zone'] = first_and_10['yardline_100'].apply(categorize_field_position)

# Get scoring outcomes by drive
drive_scores = (pbp_2023_drives
    .groupby(['game_id', 'fixed_drive'])
    .agg({
        'touchdown': lambda x: (x == 1).any(),
        'field_goal_result': lambda x: (x == 'made').any()
    })
    .reset_index()
)
drive_scores['scored_any'] = drive_scores['touchdown'] | drive_scores['field_goal_result']

# Merge and calculate probabilities
first_and_10_with_scores = first_and_10.merge(
    drive_scores, on=['game_id', 'fixed_drive'], how='left'
)

zone_order = ["Own 10", "Own 20", "Own 30", "Own 40", "Midfield",
              "Opp 40", "Opp 30", "Opp 20", "Opp 10", "Red Zone"]

scoring_by_position = (first_and_10_with_scores
    .groupby('field_zone')
    .agg(
        drives=('field_zone', 'count'),
        td_prob=('touchdown', 'mean'),
        fg_prob=('field_goal_result', 'mean'),
        any_score_prob=('scored_any', 'mean')
    )
    .reset_index()
)
scoring_by_position['field_zone'] = pd.Categorical(
    scoring_by_position['field_zone'],
    categories=zone_order,
    ordered=True
)
scoring_by_position = scoring_by_position.sort_values('field_zone')

# Visualize
fig, ax = plt.subplots(figsize=(12, 7))

x = np.arange(len(scoring_by_position))
width = 0.25

ax.bar(x - width, scoring_by_position['any_score_prob'], width,
       label='Any Score', color='#4CAF50', alpha=0.8)
ax.bar(x, scoring_by_position['fg_prob'], width,
       label='Field Goal', color='#2196F3', alpha=0.8)
ax.bar(x + width, scoring_by_position['td_prob'], width,
       label='Touchdown', color='#FF5722', alpha=0.8)

ax.set_xlabel('Starting Field Position', fontsize=12)
ax.set_ylabel('Probability of Scoring', fontsize=12)
ax.set_title('Scoring Probability by Starting Field Position\n1st & 10 situations, 2023 NFL Season',
             fontsize=14, fontweight='bold', pad=20)
ax.set_xticks(x)
ax.set_xticklabels(scoring_by_position['field_zone'], rotation=45, ha='right')
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda y, p: f'{y:.0%}'))
ax.legend(loc='upper left')
ax.grid(axis='y', alpha=0.3)
ax.set_ylim(0, 0.8)

plt.tight_layout()
plt.show()

Pace-Adjusted Metrics

Understanding Offensive Pace

Offensive pace (plays per game or seconds per play) significantly affects raw counting stats. Adjusting for pace allows fair comparisons:

#| label: pace-adjustment-r
#| message: false
#| warning: false

# Calculate pace-adjusted metrics
pace_adjusted <- pbp %>%
  filter(season == 2023, !is.na(epa), play_type %in% c("pass", "run")) %>%
  group_by(posteam, game_id) %>%
  summarise(
    offensive_plays = n(),
    total_epa = sum(epa, na.rm = TRUE),
    avg_epa = mean(epa, na.rm = TRUE),
    points_scored = sum((touchdown == 1) * 7 +
                       (field_goal_result == "made") * 3, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  group_by(posteam) %>%
  summarise(
    games = n(),
    plays_per_game = mean(offensive_plays),
    epa_per_play = mean(avg_epa),
    epa_per_game = mean(total_epa),
    points_per_game = mean(points_scored),

    # Pace-adjusted to league average (assuming ~64 plays/game)
    pace_adj_epa = epa_per_play * 64,
    pace_adj_points = (points_scored / offensive_plays) * 64,
    .groups = "drop"
  ) %>%
  arrange(desc(epa_per_play))

# Compare raw vs pace-adjusted rankings
pace_comparison <- pace_adjusted %>%
  mutate(
    rank_raw_ppg = rank(desc(points_per_game)),
    rank_pace_adj = rank(desc(pace_adj_points)),
    rank_change = rank_raw_ppg - rank_pace_adj
  ) %>%
  arrange(rank_pace_adj) %>%
  select(posteam, plays_per_game, points_per_game, pace_adj_points,
         rank_raw_ppg, rank_pace_adj, rank_change)

pace_comparison %>%
  gt() %>%
  cols_label(
    posteam = "Team",
    plays_per_game = "Plays/Game",
    points_per_game = "PPG (Raw)",
    pace_adj_points = "PPG (Pace Adj)",
    rank_raw_ppg = "Rank (Raw)",
    rank_pace_adj = "Rank (Adj)",
    rank_change = "Rank Δ"
  ) %>%
  fmt_number(
    columns = c(plays_per_game, points_per_game, pace_adj_points),
    decimals = 1
  ) %>%
  fmt_number(
    columns = c(rank_raw_ppg, rank_pace_adj, rank_change),
    decimals = 0
  ) %>%
  data_color(
    columns = rank_change,
    palette = c("#FF4136", "white", "#2ECC40"),
    domain = c(-5, 5)
  ) %>%
  tab_header(
    title = "Pace-Adjusted Offensive Rankings",
    subtitle = "Comparing raw points per game to pace-adjusted values"
  ) %>%
  tab_source_note("Pace adjusted to league average of 64 plays per game")

#| label: pace-adjustment-py
#| message: false
#| warning: false

# Calculate per-game metrics
game_stats = (pbp_2023
    .groupby(['posteam', 'game_id'])
    .agg(
        offensive_plays=('epa', 'count'),
        total_epa=('epa', 'sum'),
        avg_epa=('epa', 'mean'),
        points_scored=('touchdown', lambda x: (x == 1).sum() * 7 +
                      ((pbp_2023['field_goal_result'] == 'made').sum() * 3))
    )
    .reset_index()
)

# Calculate pace-adjusted metrics
league_avg_plays = 64  # Typical plays per game

pace_adjusted = (game_stats
    .groupby('posteam')
    .agg(
        games=('game_id', 'count'),
        plays_per_game=('offensive_plays', 'mean'),
        epa_per_play=('avg_epa', 'mean'),
        points_per_game=('points_scored', 'mean')
    )
    .reset_index()
)

pace_adjusted['pace_adj_epa'] = pace_adjusted['epa_per_play'] * league_avg_plays
pace_adjusted['pace_adj_points'] = (
    pace_adjusted['points_per_game'] / pace_adjusted['plays_per_game'] * league_avg_plays
)

# Rankings
pace_adjusted['rank_raw_ppg'] = pace_adjusted['points_per_game'].rank(ascending=False)
pace_adjusted['rank_pace_adj'] = pace_adjusted['pace_adj_points'].rank(ascending=False)
pace_adjusted['rank_change'] = pace_adjusted['rank_raw_ppg'] - pace_adjusted['rank_pace_adj']

pace_adjusted = pace_adjusted.sort_values('rank_pace_adj')

print("\nPace-Adjusted Rankings (2023):")
print(pace_adjusted[['posteam', 'plays_per_game', 'points_per_game',
                     'pace_adj_points', 'rank_change']].head(10).to_string(index=False))

Strength of Schedule Adjustments

Opponent-Adjusted Performance

Offensive performance must be evaluated in context of defensive quality faced:

#| label: sos-adjustment-r
#| message: false
#| warning: false

# First, calculate each team's defensive quality
defensive_quality <- pbp %>%
  filter(season == 2023, !is.na(epa), play_type %in% c("pass", "run")) %>%
  group_by(defteam) %>%
  summarise(
    def_epa_allowed = mean(epa),
    .groups = "drop"
  ) %>%
  mutate(
    def_epa_vs_avg = def_epa_allowed - mean(def_epa_allowed),
    def_rating = scale(def_epa_vs_avg)[,1]
  )

# Calculate opponent-adjusted offensive performance
sos_adjusted <- pbp %>%
  filter(season == 2023, !is.na(epa), play_type %in% c("pass", "run")) %>%
  left_join(defensive_quality, by = "defteam") %>%
  group_by(posteam) %>%
  summarise(
    plays = n(),
    raw_epa = mean(epa, na.rm = TRUE),
    avg_opponent_def = mean(def_epa_vs_avg, na.rm = TRUE),

    # Adjust EPA by subtracting opponent defensive quality
    # Positive opponent def EPA = weak defense, so subtract it
    adjusted_epa = raw_epa - avg_opponent_def,

    .groups = "drop"
  ) %>%
  mutate(
    raw_rank = rank(desc(raw_epa)),
    adjusted_rank = rank(desc(adjusted_epa)),
    rank_change = raw_rank - adjusted_rank
  ) %>%
  arrange(desc(adjusted_epa))

sos_adjusted %>%
  gt() %>%
  cols_label(
    posteam = "Team",
    plays = "Plays",
    raw_epa = "Raw EPA/Play",
    avg_opponent_def = "Avg Opp Def",
    adjusted_epa = "Adj EPA/Play",
    raw_rank = "Rank (Raw)",
    adjusted_rank = "Rank (Adj)",
    rank_change = "Rank Δ"
  ) %>%
  fmt_number(
    columns = c(raw_epa, avg_opponent_def, adjusted_epa),
    decimals = 3
  ) %>%
  fmt_number(
    columns = c(plays, raw_rank, adjusted_rank, rank_change),
    decimals = 0
  ) %>%
  data_color(
    columns = rank_change,
    palette = c("#FF4136", "white", "#2ECC40"),
    domain = c(-8, 8)
  ) %>%
  tab_header(
    title = "Strength of Schedule Adjusted Offensive EPA",
    subtitle = "2023 NFL Season"
  ) %>%
  tab_source_note(
    "Positive Avg Opp Def = faced weaker defenses | Adj EPA = Raw EPA - Opp Quality"
  )

#| label: sos-adjustment-py
#| message: false
#| warning: false

# Calculate defensive quality
defensive_quality = (pbp_2023
    .groupby('defteam')
    .agg(def_epa_allowed=('epa', 'mean'))
    .reset_index()
)

defensive_quality['def_epa_vs_avg'] = (
    defensive_quality['def_epa_allowed'] -
    defensive_quality['def_epa_allowed'].mean()
)

# Merge with offensive plays
pbp_with_def = pbp_2023.merge(defensive_quality, on='defteam', how='left')

# Calculate opponent-adjusted performance
sos_adjusted = (pbp_with_def
    .groupby('posteam')
    .agg(
        plays=('epa', 'count'),
        raw_epa=('epa', 'mean'),
        avg_opponent_def=('def_epa_vs_avg', 'mean')
    )
    .reset_index()
)

sos_adjusted['adjusted_epa'] = sos_adjusted['raw_epa'] - sos_adjusted['avg_opponent_def']
sos_adjusted['raw_rank'] = sos_adjusted['raw_epa'].rank(ascending=False)
sos_adjusted['adjusted_rank'] = sos_adjusted['adjusted_epa'].rank(ascending=False)
sos_adjusted['rank_change'] = sos_adjusted['raw_rank'] - sos_adjusted['adjusted_rank']

sos_adjusted = sos_adjusted.sort_values('adjusted_epa', ascending=False)

print("\nStrength of Schedule Adjusted EPA (2023):")
print(sos_adjusted[['posteam', 'raw_epa', 'avg_opponent_def',
                    'adjusted_epa', 'rank_change']].head(10).to_string(index=False))

Decomposing Offensive Contributions

Offensive Line Impact

Isolating offensive line performance from skill position performance:

#| label: oline-impact-r
#| message: false
#| warning: false

# Analyze metrics that isolate O-line performance
oline_metrics <- pbp %>%
  filter(season == 2023, !is.na(epa)) %>%
  group_by(posteam) %>%
  summarise(
    # Run blocking (before contact metrics would be ideal)
    rush_plays = sum(play_type == "run", na.rm = TRUE),
    rush_epa = mean(epa[play_type == "run"], na.rm = TRUE),
    rush_success = mean(epa[play_type == "run"] > 0, na.rm = TRUE),
    stuff_rate = mean(yards_gained[play_type == "run"] <= 0, na.rm = TRUE),

    # Pass blocking
    pass_plays = sum(play_type == "pass", na.rm = TRUE),
    sack_rate = sum(sack == 1, na.rm = TRUE) / pass_plays,
    pressures = sum(qb_hit == 1 | sack == 1, na.rm = TRUE),
    pressure_rate = pressures / pass_plays,

    # Time-related (requires additional data but approximated)
    scramble_rate = sum(qb_scramble == 1, na.rm = TRUE) / pass_plays,

    .groups = "drop"
  ) %>%
  mutate(
    # Composite O-line rating (inverse of bad outcomes)
    oline_rating = scale(-stuff_rate)[,1] * 0.3 +
                   scale(-sack_rate)[,1] * 0.4 +
                   scale(-pressure_rate)[,1] * 0.3,
    oline_rank = rank(desc(oline_rating))
  ) %>%
  arrange(desc(oline_rating))

oline_metrics %>%
  select(posteam, rush_epa, rush_success, stuff_rate,
         sack_rate, pressure_rate, oline_rating, oline_rank) %>%
  gt() %>%
  cols_label(
    posteam = "Team",
    rush_epa = "Rush EPA",
    rush_success = "Rush SR",
    stuff_rate = "Stuff%",
    sack_rate = "Sack%",
    pressure_rate = "Pressure%",
    oline_rating = "O-Line Rating",
    oline_rank = "Rank"
  ) %>%
  fmt_number(columns = c(rush_epa, oline_rating), decimals = 3) %>%
  fmt_percent(columns = c(rush_success, stuff_rate, sack_rate, pressure_rate),
              decimals = 1) %>%
  fmt_number(columns = oline_rank, decimals = 0) %>%
  data_color(
    columns = oline_rating,
    palette = "RdYlGn"
  ) %>%
  tab_header(
    title = "Offensive Line Performance Metrics",
    subtitle = "2023 NFL Season"
  ) %>%
  tab_source_note(
    "O-Line Rating: Composite of run stuffing, sack rate, and pressure rate"
  )

#| label: oline-impact-py
#| message: false
#| warning: false

# Calculate O-line metrics
def calc_oline_metrics(df):
    rush_plays = df[df['play_type'] == 'run']
    pass_plays = df[df['play_type'] == 'pass']

    return pd.Series({
        'rush_plays': len(rush_plays),
        'rush_epa': rush_plays['epa'].mean() if len(rush_plays) > 0 else 0,
        'rush_success': (rush_plays['epa'] > 0).mean() if len(rush_plays) > 0 else 0,
        'stuff_rate': (rush_plays['yards_gained'] <= 0).mean() if len(rush_plays) > 0 else 0,
        'pass_plays': len(pass_plays),
        'sack_rate': (pass_plays['sack'] == 1).sum() / len(pass_plays) if len(pass_plays) > 0 else 0,
        'pressure_rate': ((pass_plays['qb_hit'] == 1) | (pass_plays['sack'] == 1)).sum() / len(pass_plays) if len(pass_plays) > 0 else 0,
        'scramble_rate': (pass_plays['qb_scramble'] == 1).sum() / len(pass_plays) if len(pass_plays) > 0 else 0
    })

oline_metrics = (pbp_2023
    .groupby('posteam')
    .apply(calc_oline_metrics)
    .reset_index()
)

# Standardize for composite rating
scaler = StandardScaler()
oline_metrics[['stuff_z', 'sack_z', 'pressure_z']] = scaler.fit_transform(
    oline_metrics[['stuff_rate', 'sack_rate', 'pressure_rate']]
)

# Composite rating (negative because lower is better)
oline_metrics['oline_rating'] = (
    -oline_metrics['stuff_z'] * 0.3 +
    -oline_metrics['sack_z'] * 0.4 +
    -oline_metrics['pressure_z'] * 0.3
)
oline_metrics['oline_rank'] = oline_metrics['oline_rating'].rank(ascending=False)
oline_metrics = oline_metrics.sort_values('oline_rating', ascending=False)

print("\nOffensive Line Performance Metrics (2023):")
print(oline_metrics[['posteam', 'rush_epa', 'stuff_rate', 'sack_rate',
                     'pressure_rate', 'oline_rating']].head(10).to_string(index=False))

Skill Position Contributions

#| label: skill-positions-r
#| message: false
#| warning: false

# Analyze skill position performance
skill_position_metrics <- pbp %>%
  filter(season == 2023, !is.na(epa)) %>%
  group_by(posteam) %>%
  summarise(
    # Passing game (QB + receivers)
    pass_plays = sum(play_type == "pass" & sack == 0, na.rm = TRUE),
    pass_epa = mean(epa[play_type == "pass" & sack == 0], na.rm = TRUE),
    completion_pct = mean(complete_pass[play_type == "pass"], na.rm = TRUE),
    yards_per_completion = mean(yards_gained[complete_pass == 1], na.rm = TRUE),
    air_yards_per_att = mean(air_yards[play_type == "pass"], na.rm = TRUE),
    yac_per_completion = mean(yards_after_catch[complete_pass == 1], na.rm = TRUE),

    # Explosive plays (skill position excellence)
    explosive_pass_rate = mean(yards_gained[play_type == "pass"] >= 20, na.rm = TRUE),
    explosive_run_rate = mean(yards_gained[play_type == "run"] >= 10, na.rm = TRUE),

    # Red zone (where skill matters most)
    rz_plays = sum(yardline_100 <= 20 & yardline_100 > 0, na.rm = TRUE),
    rz_td_rate = sum(touchdown == 1 & yardline_100 <= 20, na.rm = TRUE) / rz_plays,

    .groups = "drop"
  ) %>%
  mutate(
    # Skill position rating
    skill_rating = scale(pass_epa)[,1] * 0.35 +
                   scale(explosive_pass_rate)[,1] * 0.25 +
                   scale(explosive_run_rate)[,1] * 0.20 +
                   scale(rz_td_rate)[,1] * 0.20,
    skill_rank = rank(desc(skill_rating))
  ) %>%
  arrange(desc(skill_rating))

skill_position_metrics %>%
  select(posteam, pass_epa, completion_pct, yac_per_completion,
         explosive_pass_rate, rz_td_rate, skill_rating, skill_rank) %>%
  gt() %>%
  cols_label(
    posteam = "Team",
    pass_epa = "Pass EPA",
    completion_pct = "Comp%",
    yac_per_completion = "YAC/Comp",
    explosive_pass_rate = "Exp Pass%",
    rz_td_rate = "RZ TD%",
    skill_rating = "Skill Rating",
    skill_rank = "Rank"
  ) %>%
  fmt_number(columns = c(pass_epa, yac_per_completion, skill_rating), decimals = 2) %>%
  fmt_percent(columns = c(completion_pct, explosive_pass_rate, rz_td_rate),
              decimals = 1) %>%
  fmt_number(columns = skill_rank, decimals = 0) %>%
  data_color(
    columns = skill_rating,
    palette = "RdYlGn"
  ) %>%
  tab_header(
    title = "Skill Position Performance Metrics",
    subtitle = "QB + Receivers + RBs, 2023 NFL Season"
  ) %>%
  tab_source_note(
    "Skill Rating: Composite of passing efficiency, explosiveness, and red zone scoring"
  )

#| label: skill-positions-py
#| message: false
#| warning: false

# Calculate skill position metrics
def calc_skill_metrics(df):
    pass_plays = df[(df['play_type'] == 'pass') & (df['sack'] == 0)]
    completions = df[df['complete_pass'] == 1]
    rz_plays = df[(df['yardline_100'] <= 20) & (df['yardline_100'] > 0)]

    return pd.Series({
        'pass_plays': len(pass_plays),
        'pass_epa': pass_plays['epa'].mean() if len(pass_plays) > 0 else 0,
        'completion_pct': (df['complete_pass'] == 1).mean() if len(df[df['play_type'] == 'pass']) > 0 else 0,
        'yac_per_completion': completions['yards_after_catch'].mean() if len(completions) > 0 else 0,
        'explosive_pass_rate': (pass_plays['yards_gained'] >= 20).mean() if len(pass_plays) > 0 else 0,
        'explosive_run_rate': (df[df['play_type'] == 'run']['yards_gained'] >= 10).mean() if len(df[df['play_type'] == 'run']) > 0 else 0,
        'rz_plays': len(rz_plays),
        'rz_td_rate': (rz_plays['touchdown'] == 1).sum() / len(rz_plays) if len(rz_plays) > 0 else 0
    })

skill_metrics = (pbp_2023
    .groupby('posteam')
    .apply(calc_skill_metrics)
    .reset_index()
)

# Standardize and create composite
scaler = StandardScaler()
skill_metrics[['pass_epa_z', 'exp_pass_z', 'exp_run_z', 'rz_z']] = scaler.fit_transform(
    skill_metrics[['pass_epa', 'explosive_pass_rate', 'explosive_run_rate', 'rz_td_rate']]
)

skill_metrics['skill_rating'] = (
    skill_metrics['pass_epa_z'] * 0.35 +
    skill_metrics['exp_pass_z'] * 0.25 +
    skill_metrics['exp_run_z'] * 0.20 +
    skill_metrics['rz_z'] * 0.20
)
skill_metrics['skill_rank'] = skill_metrics['skill_rating'].rank(ascending=False)
skill_metrics = skill_metrics.sort_values('skill_rating', ascending=False)

print("\nSkill Position Performance Metrics (2023):")
print(skill_metrics[['posteam', 'pass_epa', 'explosive_pass_rate',
                     'rz_td_rate', 'skill_rating']].head(10).to_string(index=False))

Coordinator and System Analysis

Identifying Coordinator Effects

#| label: coordinator-effects-r
#| message: false
#| warning: false

# Analyze offensive tendencies by team (proxy for coordinator philosophy)
coordinator_analysis <- pbp %>%
  filter(season == 2023, !is.na(epa), down <= 3) %>%
  filter(play_type %in% c("pass", "run")) %>%
  group_by(posteam) %>%
  summarise(
    plays = n(),

    # Play calling tendencies
    pass_rate = mean(play_type == "pass"),
    pass_rate_neutral = mean(play_type[wp >= 0.2 & wp <= 0.8] == "pass", na.rm = TRUE),

    # Early down aggression
    early_down_pass_rate = mean(play_type[down <= 2] == "pass", na.rm = TRUE),

    # Situational efficiency
    third_down_conv_rate = mean(
      series_result[down == 3] == "First down", na.rm = TRUE
    ),

    # Play action usage
    play_action_rate = sum(play_action == 1, na.rm = TRUE) /
                      sum(play_type == "pass", na.rm = TRUE),

    # Personnel groupings (if available)
    # This would require roster data

    # Creativity metrics
    motion_rate = sum(!is.na(shift_since_snap), na.rm = TRUE) / plays,

    # Tempo
    avg_seconds_to_snap = mean(time_to_snap, na.rm = TRUE),

    # Efficiency by tendency
    pass_epa = mean(epa[play_type == "pass"], na.rm = TRUE),
    run_epa = mean(epa[play_type == "run"], na.rm = TRUE),

    .groups = "drop"
  ) %>%
  mutate(
    tendency = case_when(
      pass_rate_neutral >= 0.65 ~ "Pass-Heavy",
      pass_rate_neutral >= 0.55 ~ "Pass-First",
      pass_rate_neutral >= 0.45 ~ "Balanced",
      pass_rate_neutral >= 0.35 ~ "Run-First",
      TRUE ~ "Run-Heavy"
    )
  ) %>%
  arrange(desc(pass_epa))

coordinator_analysis %>%
  select(posteam, pass_rate_neutral, play_action_rate, third_down_conv_rate,
         pass_epa, run_epa, tendency) %>%
  gt() %>%
  cols_label(
    posteam = "Team",
    pass_rate_neutral = "Neutral Pass%",
    play_action_rate = "Play Action%",
    third_down_conv_rate = "3rd Down%",
    pass_epa = "Pass EPA",
    run_epa = "Run EPA",
    tendency = "Offensive Tendency"
  ) %>%
  fmt_percent(
    columns = c(pass_rate_neutral, play_action_rate, third_down_conv_rate),
    decimals = 1
  ) %>%
  fmt_number(columns = c(pass_epa, run_epa), decimals = 3) %>%
  data_color(
    columns = c(pass_epa, run_epa),
    palette = "RdYlGn"
  ) %>%
  tab_header(
    title = "Coordinator Philosophy and Performance",
    subtitle = "2023 NFL Season"
  ) %>%
  tab_source_note("Neutral situations: Win Probability between 20% and 80%")

#| label: coordinator-effects-py
#| message: false
#| warning: false

# Filter to relevant plays
coord_data = pbp_2023[
    (pbp_2023['down'] <= 3) &
    (pbp_2023['epa'].notna())
]

# Calculate coordinator metrics
def calc_coordinator_metrics(df):
    pass_plays = df[df['play_type'] == 'pass']
    neutral = df[(df['wp'] >= 0.2) & (df['wp'] <= 0.8)]
    third_down = df[df['down'] == 3]

    return pd.Series({
        'plays': len(df),
        'pass_rate': (df['play_type'] == 'pass').mean(),
        'pass_rate_neutral': (neutral['play_type'] == 'pass').mean() if len(neutral) > 0 else 0,
        'early_down_pass_rate': (df[df['down'] <= 2]['play_type'] == 'pass').mean() if len(df[df['down'] <= 2]) > 0 else 0,
        'third_down_conv_rate': (third_down['series_result'] == 'First down').mean() if len(third_down) > 0 else 0,
        'play_action_rate': (pass_plays['play_action'] == 1).sum() / len(pass_plays) if len(pass_plays) > 0 else 0,
        'pass_epa': pass_plays['epa'].mean() if len(pass_plays) > 0 else 0,
        'run_epa': df[df['play_type'] == 'run']['epa'].mean() if len(df[df['play_type'] == 'run']) > 0 else 0
    })

coordinator_analysis = (coord_data
    .groupby('posteam')
    .apply(calc_coordinator_metrics)
    .reset_index()
)

# Categorize tendencies
def categorize_tendency(pass_rate):
    if pass_rate >= 0.65: return "Pass-Heavy"
    elif pass_rate >= 0.55: return "Pass-First"
    elif pass_rate >= 0.45: return "Balanced"
    elif pass_rate >= 0.35: return "Run-First"
    else: return "Run-Heavy"

coordinator_analysis['tendency'] = coordinator_analysis['pass_rate_neutral'].apply(categorize_tendency)
coordinator_analysis = coordinator_analysis.sort_values('pass_epa', ascending=False)

print("\nCoordinator Philosophy and Performance (2023):")
print(coordinator_analysis[['posteam', 'pass_rate_neutral', 'play_action_rate',
                            'pass_epa', 'run_epa', 'tendency']].head(10).to_string(index=False))

Era Adjustments and Historical Comparisons

Adjusting for Rule Changes and Era

Football has evolved significantly over time. Fair historical comparisons require era adjustments:

#| label: era-adjustments-r
#| message: false
#| warning: false
#| cache: true

# Calculate league-average offensive metrics by season
league_averages <- pbp %>%
  filter(!is.na(epa), play_type %in% c("pass", "run")) %>%
  group_by(season) %>%
  summarise(
    league_epa = mean(epa, na.rm = TRUE),
    league_pass_epa = mean(epa[play_type == "pass"], na.rm = TRUE),
    league_run_epa = mean(epa[play_type == "run"], na.rm = TRUE),
    league_pass_rate = mean(play_type == "pass"),
    league_completion_pct = mean(complete_pass[play_type == "pass"], na.rm = TRUE),
    league_yards_per_pass = mean(yards_gained[play_type == "pass"], na.rm = TRUE),
    .groups = "drop"
  )

# Calculate team performance relative to era
era_adjusted_performance <- pbp %>%
  filter(!is.na(epa), play_type %in% c("pass", "run")) %>%
  left_join(league_averages, by = "season") %>%
  group_by(season, posteam) %>%
  summarise(
    plays = n(),
    team_epa = mean(epa, na.rm = TRUE),
    league_epa = first(league_epa),
    epa_vs_avg = team_epa - league_epa,

    team_pass_epa = mean(epa[play_type == "pass"], na.rm = TRUE),
    league_pass_epa = first(league_pass_epa),
    pass_epa_vs_avg = team_pass_epa - league_pass_epa,

    .groups = "drop"
  ) %>%
  group_by(season) %>%
  mutate(
    epa_z_score = scale(epa_vs_avg)[,1],
    pass_epa_z_score = scale(pass_epa_vs_avg)[,1]
  ) %>%
  ungroup()

# Find best offensive seasons (era-adjusted)
best_offenses_all_time <- era_adjusted_performance %>%
  arrange(desc(epa_z_score)) %>%
  slice_head(n = 20) %>%
  select(season, posteam, team_epa, league_epa, epa_vs_avg, epa_z_score)

best_offenses_all_time %>%
  gt() %>%
  cols_label(
    season = "Season",
    posteam = "Team",
    team_epa = "Team EPA",
    league_epa = "League EPA",
    epa_vs_avg = "vs Average",
    epa_z_score = "Z-Score"
  ) %>%
  fmt_number(
    columns = c(team_epa, league_epa, epa_vs_avg, epa_z_score),
    decimals = 3
  ) %>%
  data_color(
    columns = epa_z_score,
    palette = "RdYlGn"
  ) %>%
  tab_header(
    title = "Best Offensive Seasons (Era-Adjusted)",
    subtitle = "2020-2023 NFL Seasons | Adjusted for league-average offensive performance"
  ) %>%
  tab_source_note("Z-score represents standard deviations above league average for that season")

#| label: era-adjustments-py
#| message: false
#| warning: false

# Calculate league averages by season
def calc_league_averages(df):
    pass_plays = df[df['play_type'] == 'pass']
    return pd.Series({
        'league_epa': df['epa'].mean(),
        'league_pass_epa': pass_plays['epa'].mean() if len(pass_plays) > 0 else 0,
        'league_run_epa': df[df['play_type'] == 'run']['epa'].mean() if len(df[df['play_type'] == 'run']) > 0 else 0,
        'league_pass_rate': (df['play_type'] == 'pass').mean(),
        'league_completion_pct': (pass_plays['complete_pass'] == 1).mean() if len(pass_plays) > 0 else 0
    })

league_averages = (pbp[pbp['epa'].notna()]
    .groupby('season')
    .apply(calc_league_averages)
    .reset_index()
)

# Calculate team performance vs era
team_by_season = (pbp[pbp['epa'].notna()]
    .groupby(['season', 'posteam'])
    .agg(
        plays=('epa', 'count'),
        team_epa=('epa', 'mean'),
        team_pass_epa=('epa', lambda x: pbp.loc[x.index][pbp['play_type'] == 'pass']['epa'].mean())
    )
    .reset_index()
)

# Merge with league averages
era_adjusted = team_by_season.merge(league_averages, on='season')
era_adjusted['epa_vs_avg'] = era_adjusted['team_epa'] - era_adjusted['league_epa']
era_adjusted['pass_epa_vs_avg'] = era_adjusted['team_pass_epa'] - era_adjusted['league_pass_epa']

# Calculate z-scores within each season
era_adjusted['epa_z_score'] = era_adjusted.groupby('season')['epa_vs_avg'].transform(
    lambda x: (x - x.mean()) / x.std()
)

# Best offenses (era-adjusted)
best_offenses = era_adjusted.nlargest(20, 'epa_z_score')

print("\nBest Offensive Seasons (Era-Adjusted, 2020-2023):")
print(best_offenses[['season', 'posteam', 'team_epa', 'league_epa',
                     'epa_vs_avg', 'epa_z_score']].to_string(index=False))

Visualizing League Evolution

#| label: fig-league-evolution-r
#| fig-cap: "Evolution of NFL offensive efficiency (2020-2023)"
#| fig-width: 12
#| fig-height: 7
#| message: false
#| warning: false

# Plot league trends over time
league_averages %>%
  select(season, league_epa, league_pass_epa, league_run_epa) %>%
  pivot_longer(cols = -season, names_to = "metric", values_to = "value") %>%
  mutate(
    metric = factor(metric,
                   levels = c("league_epa", "league_pass_epa", "league_run_epa"),
                   labels = c("Overall EPA", "Passing EPA", "Rushing EPA"))
  ) %>%
  ggplot(aes(x = season, y = value, color = metric, group = metric)) +
  geom_line(linewidth = 1.5, alpha = 0.8) +
  geom_point(size = 3) +
  scale_color_manual(values = c("Overall EPA" = "#2C3E50",
                                "Passing EPA" = "#3498DB",
                                "Rushing EPA" = "#E74C3C")) +
  scale_x_continuous(breaks = 2020:2023) +
  labs(
    title = "Evolution of NFL Offensive Efficiency",
    subtitle = "League-average EPA per play by season",
    x = "Season",
    y = "EPA per Play",
    color = "Metric",
    caption = "Data: nflfastR"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 14),
    legend.position = "top",
    panel.grid.minor = element_blank()
  )

#| label: fig-league-evolution-py
#| fig-cap: "Evolution of NFL offensive efficiency - Python (2020-2023)"
#| fig-width: 12
#| fig-height: 7
#| message: false
#| warning: false

# Create evolution plot
fig, ax = plt.subplots(figsize=(12, 7))

seasons = league_averages['season']
ax.plot(seasons, league_averages['league_epa'],
        marker='o', linewidth=2, markersize=8, label='Overall EPA', color='#2C3E50')
ax.plot(seasons, league_averages['league_pass_epa'],
        marker='s', linewidth=2, markersize=8, label='Passing EPA', color='#3498DB')
ax.plot(seasons, league_averages['league_run_epa'],
        marker='^', linewidth=2, markersize=8, label='Rushing EPA', color='#E74C3C')

ax.set_xlabel('Season', fontsize=12)
ax.set_ylabel('EPA per Play', fontsize=12)
ax.set_title('Evolution of NFL Offensive Efficiency\nLeague-average EPA per play by season',
             fontsize=14, fontweight='bold', pad=20)
ax.set_xticks(seasons)
ax.legend(loc='upper left', fontsize=11)
ax.grid(True, alpha=0.3)
ax.axhline(y=0, color='black', linestyle='--', alpha=0.3)

plt.tight_layout()
plt.show()

Building Predictive Models

Predicting Future Offensive Performance

#| label: predictive-model-r
#| message: false
#| warning: false

# Create dataset for prediction
# Use first half of season to predict second half
prediction_data <- pbp %>%
  filter(season == 2023, !is.na(epa), !is.na(week)) %>%
  filter(play_type %in% c("pass", "run")) %>%
  mutate(half = ifelse(week <= 9, "first_half", "second_half")) %>%
  group_by(posteam, half) %>%
  summarise(
    plays = n(),
    epa_per_play = mean(epa, na.rm = TRUE),
    success_rate = mean(epa > 0, na.rm = TRUE),
    explosive_rate = mean(abs(epa) > 1, na.rm = TRUE),
    pass_epa = mean(epa[play_type == "pass"], na.rm = TRUE),
    run_epa = mean(epa[play_type == "run"], na.rm = TRUE),
    .groups = "drop"
  ) %>%
  pivot_wider(
    names_from = half,
    values_from = c(plays, epa_per_play, success_rate, explosive_rate, pass_epa, run_epa)
  )

# Simple linear model
model <- lm(epa_per_play_second_half ~ epa_per_play_first_half +
              success_rate_first_half + explosive_rate_first_half,
            data = prediction_data)

# Model summary
summary(model)

# Predictions vs actuals
prediction_data <- prediction_data %>%
  mutate(
    predicted_epa = predict(model, prediction_data),
    prediction_error = abs(epa_per_play_second_half - predicted_epa)
  ) %>%
  arrange(prediction_error)

prediction_data %>%
  select(posteam, epa_per_play_first_half, epa_per_play_second_half,
         predicted_epa, prediction_error) %>%
  gt() %>%
  cols_label(
    posteam = "Team",
    epa_per_play_first_half = "1st Half EPA",
    epa_per_play_second_half = "2nd Half EPA",
    predicted_epa = "Predicted EPA",
    prediction_error = "Abs Error"
  ) %>%
  fmt_number(decimals = 3) %>%
  data_color(
    columns = prediction_error,
    palette = "Greens",
    reverse = TRUE
  ) %>%
  tab_header(
    title = "Predicting Second Half Performance from First Half",
    subtitle = "2023 NFL Season"
  ) %>%
  tab_source_note(
    sprintf("Model R-squared: %.3f | RMSE: %.3f",
            summary(model)$r.squared,
            sqrt(mean((prediction_data$epa_per_play_second_half - prediction_data$predicted_epa)^2)))
  )

#| label: predictive-model-py
#| message: false
#| warning: false

from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, mean_squared_error

# Create prediction dataset
pbp_2023_pred = pbp_2023[pbp_2023['week'].notna()].copy()
pbp_2023_pred['half'] = pbp_2023_pred['week'].apply(lambda x: 'first_half' if x <= 9 else 'second_half')

# Calculate metrics by half
def calc_half_metrics(df):
    pass_plays = df[df['play_type'] == 'pass']
    run_plays = df[df['play_type'] == 'run']

    return pd.Series({
        'plays': len(df),
        'epa_per_play': df['epa'].mean(),
        'success_rate': (df['epa'] > 0).mean(),
        'explosive_rate': (df['epa'].abs() > 1).mean(),
        'pass_epa': pass_plays['epa'].mean() if len(pass_plays) > 0 else 0,
        'run_epa': run_plays['epa'].mean() if len(run_plays) > 0 else 0
    })

half_stats = (pbp_2023_pred
    .groupby(['posteam', 'half'])
    .apply(calc_half_metrics)
    .reset_index()
)

# Pivot to wide format
prediction_data = half_stats.pivot(index='posteam', columns='half')
prediction_data.columns = ['_'.join(col).strip() for col in prediction_data.columns.values]
prediction_data = prediction_data.reset_index()

# Build model
X = prediction_data[['epa_per_play_first_half', 'success_rate_first_half',
                     'explosive_rate_first_half']].fillna(0)
y = prediction_data['epa_per_play_second_half'].fillna(0)

model = LinearRegression()
model.fit(X, y)

# Predictions
prediction_data['predicted_epa'] = model.predict(X)
prediction_data['prediction_error'] = abs(
    prediction_data['epa_per_play_second_half'] - prediction_data['predicted_epa']
)

r2 = r2_score(y, prediction_data['predicted_epa'])
rmse = np.sqrt(mean_squared_error(y, prediction_data['predicted_epa']))

print("\nPredicting Second Half Performance from First Half (2023):")
print(f"Model R-squared: {r2:.3f}")
print(f"RMSE: {rmse:.3f}\n")

prediction_data_sorted = prediction_data.sort_values('prediction_error')
print(prediction_data_sorted[['posteam', 'epa_per_play_first_half',
                              'epa_per_play_second_half', 'predicted_epa',
                              'prediction_error']].head(10).to_string(index=False))

Complete Offensive Report Card

Building a Comprehensive Evaluation

Let's create a complete offensive scouting report combining all our analyses:

#| label: report-card-r
#| message: false
#| warning: false

# Combine all metrics into comprehensive report card
offensive_report_card <- offensive_efficiency %>%
  left_join(
    drive_efficiency %>% select(posteam, points_per_drive, scoring_drive_rate),
    by = "posteam"
  ) %>%
  left_join(
    oline_metrics %>% select(posteam, sack_rate, stuff_rate),
    by = "posteam"
  ) %>%
  left_join(
    skill_position_metrics %>% select(posteam, explosive_pass_rate, rz_td_rate),
    by = "posteam"
  ) %>%
  left_join(
    offensive_rating %>% select(posteam, rating_100),
    by = "posteam"
  ) %>%
  mutate(
    # Letter grades for each component
    epa_grade = case_when(
      epa_per_play >= quantile(epa_per_play, 0.9) ~ "A",
      epa_per_play >= quantile(epa_per_play, 0.75) ~ "B",
      epa_per_play >= quantile(epa_per_play, 0.5) ~ "C",
      epa_per_play >= quantile(epa_per_play, 0.25) ~ "D",
      TRUE ~ "F"
    ),
    drive_grade = case_when(
      points_per_drive >= quantile(points_per_drive, 0.9, na.rm = TRUE) ~ "A",
      points_per_drive >= quantile(points_per_drive, 0.75, na.rm = TRUE) ~ "B",
      points_per_drive >= quantile(points_per_drive, 0.5, na.rm = TRUE) ~ "C",
      points_per_drive >= quantile(points_per_drive, 0.25, na.rm = TRUE) ~ "D",
      TRUE ~ "F"
    ),
    oline_grade = case_when(
      sack_rate <= quantile(sack_rate, 0.1, na.rm = TRUE) ~ "A",
      sack_rate <= quantile(sack_rate, 0.25, na.rm = TRUE) ~ "B",
      sack_rate <= quantile(sack_rate, 0.5, na.rm = TRUE) ~ "C",
      sack_rate <= quantile(sack_rate, 0.75, na.rm = TRUE) ~ "D",
      TRUE ~ "F"
    ),
    explosive_grade = case_when(
      explosive_pass_rate >= quantile(explosive_pass_rate, 0.9, na.rm = TRUE) ~ "A",
      explosive_pass_rate >= quantile(explosive_pass_rate, 0.75, na.rm = TRUE) ~ "B",
      explosive_pass_rate >= quantile(explosive_pass_rate, 0.5, na.rm = TRUE) ~ "C",
      explosive_pass_rate >= quantile(explosive_pass_rate, 0.25, na.rm = TRUE) ~ "D",
      TRUE ~ "F"
    )
  ) %>%
  arrange(desc(rating_100)) %>%
  select(posteam, rating_100, epa_per_play, epa_grade, points_per_drive,
         drive_grade, sack_rate, oline_grade, explosive_pass_rate, explosive_grade)

# Display report card
offensive_report_card %>%
  gt() %>%
  cols_label(
    posteam = "Team",
    rating_100 = "Overall",
    epa_per_play = "EPA/Play",
    epa_grade = "Grade",
    points_per_drive = "Pts/Drive",
    drive_grade = "Grade",
    sack_rate = "Sack%",
    oline_grade = "Grade",
    explosive_pass_rate = "Explosive%",
    explosive_grade = "Grade"
  ) %>%
  tab_spanner(label = "Efficiency", columns = c(epa_per_play, epa_grade)) %>%
  tab_spanner(label = "Drive Success", columns = c(points_per_drive, drive_grade)) %>%
  tab_spanner(label = "O-Line", columns = c(sack_rate, oline_grade)) %>%
  tab_spanner(label = "Explosiveness", columns = c(explosive_pass_rate, explosive_grade)) %>%
  fmt_number(columns = rating_100, decimals = 0) %>%
  fmt_number(columns = c(epa_per_play, points_per_drive), decimals = 2) %>%
  fmt_percent(columns = c(sack_rate, explosive_pass_rate), decimals = 1) %>%
  data_color(
    columns = rating_100,
    palette = "RdYlGn",
    domain = c(0, 100)
  ) %>%
  data_color(
    columns = c(epa_grade, drive_grade, oline_grade, explosive_grade),
    palette = c("F" = "#FF4136", "D" = "#FF851B", "C" = "#FFDC00",
                "B" = "#7FDBFF", "A" = "#2ECC40")
  ) %>%
  tab_header(
    title = "Complete Offensive Report Card",
    subtitle = "2023 NFL Season | Comprehensive Multi-Metric Evaluation"
  ) %>%
  tab_source_note("Overall rating combines EPA efficiency, consistency, explosiveness, and discipline")

#| label: report-card-py
#| message: false
#| warning: false

# Combine all metrics
report_card = offensive_efficiency.copy()
report_card = report_card.merge(
    drive_efficiency[['posteam', 'points_per_drive', 'scoring_drive_rate']],
    on='posteam', how='left'
)
report_card = report_card.merge(
    oline_metrics[['posteam', 'sack_rate', 'stuff_rate']],
    on='posteam', how='left'
)
report_card = report_card.merge(
    skill_metrics[['posteam', 'explosive_pass_rate', 'rz_td_rate']],
    on='posteam', how='left'
)
report_card = report_card.merge(
    offensive_rating[['posteam', 'rating_100']],
    on='posteam', how='left'
)

# Assign letter grades
def assign_grade(value, percentiles, reverse=False):
    if pd.isna(value):
        return 'N/A'
    if reverse:  # For metrics where lower is better
        if value <= percentiles[0.1]: return 'A'
        elif value <= percentiles[0.25]: return 'B'
        elif value <= percentiles[0.5]: return 'C'
        elif value <= percentiles[0.75]: return 'D'
        else: return 'F'
    else:  # For metrics where higher is better
        if value >= percentiles[0.9]: return 'A'
        elif value >= percentiles[0.75]: return 'B'
        elif value >= percentiles[0.5]: return 'C'
        elif value >= percentiles[0.25]: return 'D'
        else: return 'F'

report_card['epa_grade'] = report_card['epa_per_play'].apply(
    lambda x: assign_grade(x, report_card['epa_per_play'].quantile([0.1, 0.25, 0.5, 0.75, 0.9]))
)
report_card['drive_grade'] = report_card['points_per_drive'].apply(
    lambda x: assign_grade(x, report_card['points_per_drive'].quantile([0.1, 0.25, 0.5, 0.75, 0.9]))
)
report_card['oline_grade'] = report_card['sack_rate'].apply(
    lambda x: assign_grade(x, report_card['sack_rate'].quantile([0.1, 0.25, 0.5, 0.75, 0.9]), reverse=True)
)
report_card['explosive_grade'] = report_card['explosive_pass_rate'].apply(
    lambda x: assign_grade(x, report_card['explosive_pass_rate'].quantile([0.1, 0.25, 0.5, 0.75, 0.9]))
)

# Sort and display
report_card = report_card.sort_values('rating_100', ascending=False)

print("\nComplete Offensive Report Card (2023):")
print("\nTop 15 Teams:")
print(report_card[['posteam', 'rating_100', 'epa_per_play', 'epa_grade',
                   'points_per_drive', 'drive_grade', 'sack_rate', 'oline_grade',
                   'explosive_pass_rate', 'explosive_grade']].head(15).to_string(index=False))

Summary

In this chapter, we built a comprehensive offensive evaluation framework that synthesizes all the concepts from Part II. We learned how to:

  1. Combine Multiple Metrics: Created composite ratings that balance efficiency, consistency, explosiveness, and discipline

  2. Calculate Advanced Metrics: Implemented DVOA-style adjustments and comprehensive efficiency measures

  3. Analyze Drive Efficiency: Evaluated how effectively offenses convert possessions into points

  4. Adjust for Context: Applied pace and strength-of-schedule adjustments for fair comparisons

  5. Decompose Performance: Isolated offensive line and skill position contributions

  6. Identify System Effects: Analyzed coordinator philosophies and play-calling tendencies

  7. Compare Across Eras: Made historical comparisons by adjusting for league-average performance

  8. Build Predictive Models: Used first-half performance to project second-half results

  9. Create Report Cards: Synthesized all analyses into comprehensive offensive evaluations

Key Takeaways

- No single metric tells the complete story—use multiple dimensions - Context matters: Always adjust for opponent, pace, and era - Decomposition helps identify specific strengths and weaknesses - Predictive models require balancing complexity with sample size - Professional-grade analysis combines quantitative metrics with qualitative insight

Exercises

Conceptual Questions

  1. Metric Selection: Explain why a composite rating might better evaluate offensive performance than EPA/play alone. What dimensions does EPA miss?

  2. Context Dependency: A team averages 0.15 EPA/play but faced the league's easiest schedule. Another averages 0.12 EPA/play against the hardest schedule. Which offense is better? Explain your reasoning.

  3. Era Comparisons: The 2007 Patriots averaged 0.29 EPA/play in an era where league average was 0.02. The 2020 Chiefs averaged 0.22 EPA/play when league average was 0.06. Which offense was more dominant relative to their era?

  4. System vs Talent: How would you distinguish between an offense succeeding due to superior talent versus superior scheme/play-calling?

Coding Exercises

Exercise 1: Custom Composite Rating

Create your own composite offensive rating using different weights than those shown in this chapter. a) Experiment with different weight combinations b) Validate your weights by predicting future game outcomes c) Compare your rating to the one presented in the chapter **Bonus**: Use optimization to find weights that maximize predictive accuracy

Exercise 2: Detailed Strength of Schedule

Build a more sophisticated strength-of-schedule adjustment that accounts for: a) Opponent defensive EPA allowed b) Home/away adjustments c) Weather conditions d) Opponent rest (days since last game) Compare raw rankings to fully-adjusted rankings and identify the biggest movers.

Exercise 3: Coordinator Stability Analysis

Analyze teams that changed offensive coordinators during the 2020-2023 period: a) Calculate before/after performance metrics b) Identify which coordinator changes had the biggest impact c) Determine if personnel or scheme was the primary driver of changes **Hint**: You'll need to research which teams changed coordinators and when

Exercise 4: Playoff Performance Prediction

Build a model that predicts playoff offensive performance using regular season metrics: a) Calculate regular season comprehensive ratings for playoff teams b) Compare to actual playoff EPA/play c) Identify which regular season metrics best predict playoff success d) Discuss why certain teams over/under-perform in playoffs **Hint**: Consider that playoff defenses are typically stronger

Exercise 5: Complete Team Scouting Report

Choose an NFL team and create a comprehensive offensive scouting report including: a) Overall efficiency metrics and rankings b) Strengths and weaknesses by unit (line, QB, receivers, RBs) c) Play-calling tendencies and situational performance d) Historical trends (improving or declining?) e) Opponent-adjusted performance f) Recommendations for defensive game planning Present your analysis in a professional format with visualizations.

Further Reading

Academic Papers

  • Yurko, R., Ventura, S., & Horowitz, M. (2019). "nflWAR: A Reproducible Method for Offensive Player Evaluation in Football." Journal of Quantitative Analysis in Sports, 15(3), 163-183.

  • Lopez, M. J., Matthews, G. J., & Baumer, B. S. (2018). "How often does the best team win? A unified approach to understanding randomness in North American sport." The Annals of Applied Statistics, 12(4), 2483-2516.

  • Burke, B. (2009). "Calibration of Win Probability Models." Advanced NFL Stats.

Industry Resources

  • Football Outsiders. (2023). "DVOA Explained." https://www.footballoutsiders.com/info/methods

  • Baldwin, B. (2020). "Fundamentals of EPA: Expected Points and EPA Explained." https://www.nfeloapp.com/analysis/

  • Sharp Football. (2023). "Comprehensive Team Analysis Methodology."

Books

  • Alamar, B. (2013). Sports Analytics: A Guide for Coaches, Managers, and Other Decision Makers. Columbia University Press.

  • Winston, W. (2012). Mathletics: How Gamblers, Managers, and Sports Enthusiasts Use Mathematics. Princeton University Press.

References

:::