The puck is cylindrical

There is a German phrase about soccer: Der Ball ist rund translating into "the ball is round". Anything can happen, the game resists certainty, and it is always a probability function without absolutes. Hockey takes that same premise with an amplifier. It slides across a slick surface, flips, stands on edges, disappears under skates, and deflects off of pads before going in the net.

That makes hockey analytics a credit assignment problem with a long list of extenuating circumstances. The sport is shift-oriented and deployment-oriented and succombs to a randomness that makes outcomes for less talented teams more likely. Goalies contort and flail their impossibly nimble limbs. Power plays and empty nets create games within the game. The points system attributes goal scorers and at most two assistors, but this is often enabled by a dozen invisible decisions that never touch a scoresheet... The defenseman who killed a rush thirty seconds earlier, the smooth skater who did the clean zone entry upon the prior shift leading to an offensive zone draw, the winger whose forecheck forced the turnover, the center who lost the faceoff but won the board battle after.

Plus-minus, a naive approach for on-ice credit, is so influenced by teammates, opponents, deployment, and goaltending that it mostly measures whose line you're on. The modern answer is RAPM, Regularized Adjusted Plus-Minus, which stops asking "who scored" and starts asking "what changes when you're on the ice, after accounting for everyone else who is also on the ice." It helps isolate impact for specific player contribution, but can also be misleading for players that always are on this ice together or players with sheltered offensive zone deployments. Answering that question properly requires reconstructing every second of a season, which is a data engineering problem before it's ever a statistics problem.

So let's try to make sense of the ice rink. nhl_2526 is a namespace within our oleander oleander (yes, we have accountId 1 of course) account that pulls the entire 2025-26 NHL season from the league's public API into oleander, builds an expected-goals model, chops 1,312 games into 518,282 stints, and solves one very large, very sparse ridge regression. This post walks through the whole thing, including the actual math on one Celebrini shift.

Pulling down a season

Everything starts at the NHL's public API (api-web.nhle.com). Most of the ingestion layer was prompt-driven from the NHL API documentation: oleander Python was written for each endpoint to conform to the API's actual response shape, then checked against representative responses and iterated until the data landed cleanly. The result is a series of inspectable raw tables - raw_players, raw_schedule, raw_skater_game_logs, raw_goalie_game_logs, the NHL Edge tables, raw_pbp_shots, and raw_shifts - rather than one prematurely modeled dataset. A single entrypoint.py orchestrates the PySpark jobs that populate those tables in the oleander.nhl_2526 namespace; the joins, feature engineering, and models begin downstream of that raw boundary.

The ingestion pattern is deliberately boring. We fetched JSON with a 1/2 second delay between requests and exponential backoff on HTTP 429, then Spark takes over for everything downstream:

# nhl/ingest_pbp.py (abridged)
for game_id in game_ids:
    payload = fetch_with_backoff(
        f"https://api-web.nhle.com/v1/gamecenter/{game_id}/play-by-play"
    )
    rows.extend(parse_shot_events(payload))

df = spark.createDataFrame(rows, schema=PBP_SHOT_SCHEMA)
df.writeTo("oleander.nhl_2526.raw_pbp_shots").createOrReplace()

Deployment to oleander's managed Spark is a Makefile and two commands. make zips the nhl module and packs a virtualenv, then:

# the agent ran this, not me
oleander spark jobs upload entrypoint.py \
  --py-files out/pyfiles.zip., \
  --virtualenv out/environment.tar.gz

oleander spark jobs submit entrypoint.py \
  --namespace nhl_2526 --name nhl-full-pipeline

Agents and the oleander MCP built and ran this entire workflow looping through a couple of errors until it happily validated full game and roster collections. Because every job runs through oleander, each table write emits lineage events automatically, so when a downstream WAR number looks suspicious you can walk the lineage graph from skater_percentiles back through rapm_skater_results to the exact ingest_shifts run that fed it, with logs and traces attached.

The full pipeline is five phases, and order matters, because everything compounds:

Phase 1: Ingest          rosters, schedule, game logs, Edge, PBP, shifts
Phase 2: xG + stints     build_xg (expected goals), build_stint_events, game context
Phase 3: RAPM            ridge regression on 5v5 stints
Phase 4: Player models   GAx, (goals above expected) GSAx, (goals saved above expected, for goalies) WAR components (wins above replacement)
Phase 5: Rankings        percentiles within F / D / G

Stints: the atomic unit of on-ice credit

A shift is one player's trip over the boards. A stint is the intersection of everybody's shifts: a maximal slice of game time during which the ten skaters and two goalies (for 5v5 play) on the ice do not change. Any substitution, even one winger swapping on the fly, ends the current stint and starts a new one. The ingest_shifts job takes 986,183 raw shift records, turns every shift start and end into a boundary, and slices each game along all of them.

The 2025-26 season decomposes into 518,282 stints, of which 445,916 are 5v5, totaling 1,064 hours of even-strength hockey. Their durations are wildly skewed:

The median 5v5 stint is 3 seconds long. That sounds broken until you watch a line change: players trickle on and off one at a time, so a single change produces several short-lived intermediate units. The regression handles this gracefully because every stint is weighted by its duration; a 2-second stint contributes almost nothing.

One shift under the microscope

Here is Macklin Celebrini's fifth shift from San Jose at Detroit on January 16, 2026. He hops over the boards at 9:04 of the first period with Will Smith and Pavol Regenda, Dmitry Orlov and John Klingberg behind them. The Sharks keep the same five skaters out for 50 seconds while Detroit changes all three forwards halfway through:

San Jose Sharksaway · attacking unit

Pavol Regenda

LW

Macklin Celebrini

C

Will Smith

C

Dmitry Orlov

D

John Klingberg

D

Yaroslav Askarov

G

VS

Detroit Red Wingshome · final defending unit

Patrick Kane

RW

Andrew Copp

C

Alex DeBrincat

RW

Moritz Seider

D

Simon Edvinsson

D

John Gibson

G

The puck stays at one end. Celebrini tips the first look wide at 9:06. Smith misses a sharp-angle backhand six seconds later. After Detroit changes forwards, resulting in a new stint, Celebrini puts a snap shot on John Gibson at 9:50. Four seconds later Smith taps in a one-foot backhand, Celebrini gets the primary assist, and the shift and stint ends. The rink view reconstructs all four attempts:

The official NHL clip shows how the sequence ends:

Four attempts, 1.138 expected goals, and a wonderous clean zone entry. The first three chances total only 0.188 xG. Smith's final backhand is worth 0.95 by itself: one foot from the line, no angle, and no plausible world in which the goalie should recover. The shooting calibration mechanics explained later will dive into some of the nuances here. San Jose's teenage savior is on the ice for all four attempts, takes two, and creates the last one.

This is also a clean demonstration of why a shift is not a stint. Celebrini's 50 seconds become two 25-second regression observations because Detroit changes from Larkin–Raymond–van Riemsdyk to Kane–Copp–DeBrincat around 9:29. The final stint's corrected ledger is below; columns are from Detroit's home perspective:

duration_seconds = 25          strength_state = 5v5
xgf_home = 0.0000              xga_home = 1.0069
gf_home  = 0                   ga_home  = 1
cf_home  = 0                   ca_home  = 2

The 1.0069 is Celebrini's 0.0569 snap shot plus Smith's 0.95 goal.

The xG math

Every unblocked shot attempt gets a goal probability from a Spark ML logistic regression (xg_model = logistic_v2 in the current build), trained on all attempts with a goalie in net:

P(goal | shot) = σ(β₀ + β·f)     where σ(z) = 1 / (1 + e^(-z))

f = [ log(1 + distance_ft),      # distance decays roughly exponentially
      angle_deg,                 # sharper angle, smaller net
      score_diff (clamped ±3),   # score effects on shot selection
      period,                    # period of game - more likely to score in 2nd period due to longer shift change
      is_rebound,                # secondary shot
      is_rush,                   # in transition, odd numbered sitution
      crossed_royal_road,        # lateral pass before shot
      is_behind_net,
      is_tip,                    # shot that is tipped in front of the goal, changing direction
      is_slap,
      is_wraparound,             # shooting with body from behind the net
      is_pp,                     # is power play
      is_sh ]

Celebrini's snap shot and Smith's one-foot backhand show the tradeoff between angle and distance. Celebrini's attempt is nearly straight on at 1.3 degrees, which helps, but it is still 44 feet from the goal. The model sees log(1 + 44.0) = 3.81, so the shot stays in low-probability territory at 0.057 xG. Smith's finish is the opposite case: no angle penalty, no meaningful distance, and the puck already at the edge of the crease. After calibration, that kind of chance is allowed to hit the model ceiling at 0.950 xG.

Raw logistic output isn't the final xG. It gets calibrated twice:

1. Per strength bucket (5v5, 5v4, 4v5, other): probabilities are scaled so each bucket's total xG equals its total goals. Power plays finish better than the same shot at evens, and the base model underrates that.
2. Globally on shots-on-goal with a goalie in net, so that league-wide Goals Saved Above Expected sums to approximately zero. If you skip this, every goalie in the league looks collectively brilliant or collectively cursed, and goalie WAR inherits the bias.

The receipts, by decile:

Empty-net situations skip the model entirely and use a multiplicative heuristic, because a logistic regression trained on goalie-in-net shots has no idea what to do when there's no goalie:

xg_raw = exp(-distance * 0.04) * angle_factor * shot_type_factor
       * rebound_factor * rush_factor * royal_road_factor

Only 1,822 of 112,888 attempts take this path. And where do the other 111,066 come from? Everywhere, but not uniformly:

Toggle that chart to "finishing" and the xG model stops being abstract: shot volume piles up at the blue line where little goes in, and conversion concentrates in a tight cone in front of the crease. Smith's goal sits at the brightest edge of that cone; Celebrini's preceding snap shot came from the high-volume, low-yield interior where goalies are expected to makes saves from. Distance and angle carry most of the model because distance and angle carry most of hockey.

The RAPM iteration

Now the main event. RAPM asks: holding everyone else constant, how does xG flow change when a given player is on the ice? Mechanically, every stint becomes two rows in a design matrix, one per attacking perspective. Each row is really just a roster: five attacking skaters, five defending skaters, the defending goalie, team context, and home ice. The final 25 seconds of this shift say that San Jose generated 120.8 xGF/60 with the unit shown above, while those same twelve players pointed the other way produced zero. The targets are per-60 rates, and this stint is short enough that the safety valve in the conversion fires:

rate_factor = 3600 / max(duration, 30)         # 30s floor for tiny stints

duration = 25s → rate_factor = 3600 / 30 = 120   (floor applied)

Row 1 (SJS attacking):  y_xgf = 1.0069 × 120 = 120.8 xGF/60
                        y_gf  = 1      × 120 = 120.0 GF/60
Row 2 (DET attacking):  y_xgf = 0      × 120 = 0.0 xGF/60
                        y_gf  = 0.0

weight w = 25 / 60 ≈ 0.42 (stint duration in minutes)

Without the floor, the goal would extrapolate to 144 GF/60, as if this unit scores every 25 seconds. Even with it, 120.8 xGF/60 is absurd; a good team sustains about 2.7 over a full game. That is fine. This row carries about 0.42 minutes of weight in a regression built from 1,064 hours of hockey. Individually it is an anecdote; collectively, 891,832 of these rows are a season.

Assemble every 5v5 stint the same way and you get the full system: X is 891,832 rows by 2,043 columns (940 skaters twice, once as offense and once as defense, plus 98 goalies, 32 teams twice, and one home-ice column), W is a diagonal matrix of stint durations in minutes, and y is the per-60 rate vector. Ordinary least squares dies here, because teammates who always play together are nearly collinear and low-minute players produce wild coefficients. So we solve the ridge-regularized weighted normal equations instead:

β = (XᵀWX + λI)⁻¹ XᵀWy

λ = 800  for the xGF/60 target
λ = 2500 for the GF/60 target
λ × 0.25 for team and home-ice columns

Each choice is doing a specific job:

• W (minute weights): a 25-second observation deserves more trust than a 2-second one, and a 60-second one deserves more than either. Weighting in minutes keeps the ridge penalty on the same scale as the sample size a regular sees over a season. Longer stints are also less contaminated by line-change noise.
• λ = 800 vs λ = 2500: shrinkage strength tracks noise. Shot generation (xGF) is a high-frequency, fairly repeatable process; goals (GF) are rare events ruled by shooting and save percentage luck, so the GF model gets more than three times the shrinkage. Smith's actual goal is worth 120 GF/60 in row 1, but λ = 2500 ensures no player's reputation rides on which bounces happen to go in.
• λ × 0.25 on team/home controls: the nuisance columns get a lighter penalty on purpose. If team effects were shrunk as hard as player effects, systematic team context (coaching, systems, schedule) would leak into every player's coefficient instead of being absorbed by the team column.

Those λ values were not picked by vibes alone. A chronological 80/20 game holdout plus a λ sweep lands nearest λ = 800 for xGF (about 0.22% skill versus a stint-mean baseline). Nearby λ values shuffle the board only a little; the absolute lift is small because most stint variance is noise.

After solving, offensive and defensive coefficient blocks are centered on their means, defense is sign-flipped so positive always means good, and each skater gets:

rapm_off_xgf_60 = β_off[i] − mean(β_off)       # xG created above average
rapm_def_xgf_60 = −(β_def[i] − mean(β_def))    # xG suppressed above average
rapm_xgf_60     = rapm_off_xgf_60 + rapm_def_xgf_60

Drive the multiplication yourself

Reading about a dot product is one thing; turning the crank is better. The widget below computes one row of Xβ for our stint using each player's actual published 2025-26 coefficient. Step through the terms one at a time, click any player to bench them (their x flips from 1 to 0 and their term drops out), and flip the perspective to compute the other row. The published coefficients are centered on the league mean, so the team, home-ice, and centering constants are folded into a single baseline term of 2.61, the league-average 5v5 xGF/60:

RAPM contribution waterfall

SJS attackingDET attacking

How the lineup moves the prediction

Start at league context, then add each player's isolated RAPM effect.

ŷ = 2.98 xGF/60

raises predictionlowers predictionClick a player bar to toggle their X entry

Show the underlying matrix math

X (1 × 12) × β (12 × 1) = ŷ (1 × 1)

One attacking row of the design matrix multiplied by its coefficient vector.

2.98 xGF/60

X: design row

Baseline

Regenda

Celebrini

Smith

Orlov

Klingberg

Kane

Copp

DeBrincat

Seider

Edvinsson

Gibson

1

1

1

1

1

1

1

1

1

1

1

1

×

β: coefficients

Baseline

Regenda

Celebrini

Smith

Orlov

Klingberg

Kane

Copp

DeBrincat

Seider

Edvinsson

Gibson

+2.61

+0.06

+0.13

+0.06

+0.09

+0.12

+0.01

−0.06

+0.01

−0.15

−0.00

+0.10

=

ŷ

2.98

xGF/60

The production row has 2,043 columns and 14 non-zero entries. This compact view shows the ten player columns, goalie column, and one baseline term that folds together the two team columns, home ice, and centering constants.

Lineup coefficients · click a player to toggle x between 1 and 0

β₀Baselineleague-average context: team, home ice, and centering constantsx = 12.61+2.61RegendaLW · offense columnx = 1β_off +0.06+0.06CelebriniC · offense columnx = 1β_off +0.13+0.13SmithC · offense columnx = 1β_off +0.06+0.06OrlovD · offense columnx = 1β_off +0.09+0.09KlingbergD · offense columnx = 1β_off +0.12+0.12KaneRW · defense columnx = 1β_def −0.01+0.01CoppC · defense columnx = 1β_def +0.06−0.06DeBrincatRW · defense columnx = 1β_def −0.01+0.01SeiderD · defense columnx = 1β_def +0.15−0.15EdvinssonD · defense columnx = 1β_def +0.00−0.00GibsonG · goalie columnx = 1β_g −0.10+0.10

Multiply next termShow allReset

ŷ = 2.98 xGF/60

model expectation ŷ

3.0

observed y in this stint (w = 25 of 3.8M seconds)

120.8

With everyone on the ice, the model expects this San Jose unit to generate about 2.98 xGF/60. Start at the 2.61 league baseline: the five Sharks' offensive terms add 0.46, Detroit's defensive terms subtract 0.19, and Gibson gives 0.10 back. The observed row is 120.8. That residual looks ridiculous because it is, and then it is multiplied by only 0.42 minutes of weight out of roughly 128,000. One spectacular shift argues with the model; it does not win.

So what did the final 25 seconds in Detroit actually do? They nudged five Sharks offensive coefficients and five Red Wings defensive coefficients toward “ice tilts when these ones are out there,” by an amount proportional to 0.42 minutes of weight and inversely proportional to everything else those ten players did across the rest of the season. One stint moves a coefficient by a rounding error. A thousand shifts of consistently winning your minutes moves it to the top of this chart:

What the model actually thinks

A few things that are interesting in the results, in the spirit of arguing with a spreadsheet and ML model:

There are two species of RAPM leader. Brandon Hagel (+0.33) and Jesper Bratt (+0.32) get there mostly by creating: their offensive component dwarfs their defensive one. Jordan Kyrou (+0.38) and Moritz Seider (+0.30) get to the same neighborhood with a more even offense/defense split. Classic plus-minus could never distinguish these players from their linemates; the whole point of the design matrix is that a defensive column gets credit only for suppression that persists after goalie, partners, and team context are accounted for.

The leaderboard is still a little weird. Alongside Kyrou, Hagel, and Bratt sit Jordan Spence and Linus Karlsson. Some of that is real (players who quietly win their minutes), and some of it is single-season RAPM being noisy for unusual deployment. This is why the WAR model downstream blends RAPM at 80% with individual shot creation rather than swallowing it whole, and why many sports-oriented data scientists pool multiple seasons. Season-long RAPM can bias toward a very specific set of linemates, coaching strategies, and deployments.

xGF and GF RAPM disagree, and that's the design working. Bratt's rapm_gf_60 (+0.01) is a fraction of his xG-based number (+0.32). Somewhere between shot quality and goals, finishing luck in his minutes ran cold. The paired targets with different λ values let you see the luck rather than average it away.

From coefficients to player cards

RAPM is deliberately narrow: it estimates how 5v5 chance flow changes with a player on the ice. WAR asks a broader question: how many wins did the player add relative to a replacement-level skater? In this build, replacement level for skaters is the 25th-percentile component rate within position group, forwards and defensemen separately, among skaters with at least 10 games.

The version here is a small, transparent build rather than a claim to have solved hockey. It turns several player signals into the same currency, wins, then adds them:

Each raw signal is compared within a position group, translated into goal value, then divided by 5.15 goals per win. So a component worth roughly +2.6 goals becomes about +0.5 WAR, and a +5.15-goal edge becomes +1.0 WAR. After that conversion, the 25th-percentile skater baseline is shifted to zero, so positive numbers are wins above replacement and negative numbers are wins below that replacement line. That is why a player can sit near average in 5v5 RAPM and still grade out as excellent if the other columns are doing enough work.

For the three Sharks in the featured shift, the accounting looks like this:

Additive WAR ramp

Every player starts at zero WAR. Each component moves the running total; teal steps add value, red steps subtract value.

positive WARnegative WAR

Macklin Celebrini

C · 82 games

+4.41 WAR

final total

Will Smith

C · 69 games

+2.10 WAR

final total

William Eklund

LW · 78 games

+1.48 WAR

final total

Each step is already measured in wins above replacement, so the components can be added left to right. Celebrini's 4.41 WAR is not a claim that his 5v5 RAPM alone made the season. His even-strength RAPM is solidly above average; the model still credits a 99th-percentile season mostly because elite playmaking, power-play value, and finishing do the heavy lifting. Eklund's card tells the inverse finishing story: transition and playmaking survive a negative shooting component.

Choose any skater from the 2025–26 season to render their card:

Choose a playerMCMacklin Celebrini· SJS · C

Loading the full 2025–26 player list…

Macklin Celebrini

SJS · C · 82 games

off_def_xgf_team_ctx_v2 · 1313 min 5v5

4.41

WAR

99

WAR percentile

99

Composite

EV offense+0.50 WAR

EV defense+0.04 WAR

Transition+0.33 WAR

Power play+0.94 WAR

Playmaking+1.15 WAR

Shooting+1.12 WAR

Penalties+0.33 WAR

−2.0 WAR0+2.0 WAR

Shared scale for every component and every player

A 99th-percentile WAR season: above-average even-strength RAPM offense, with the bulk still coming from playmaking, power-play creation, and finishing.

Goalie cards

Goalies need their own accounting. A skater's value can be split across offense, defense, transition, power play, playmaking, shooting, and penalties. A goalie's job is quite a bit narrower: face a sequence of shots, inherit the chaos of screens and rebounds and rushes, and turn expected goals against into actual goals against. It borrows from the xG models that we created before and used for skaters, so let's use it again here.

The goalie card starts with GSAx, Goals Saved Above Expected. For every shot with a goalie in net, the xG model estimates the chance of a goal. Add those probabilities up and compare them with the goals actually allowed:

GSAx = expected goals against − actual goals against

Then the model shrinks GSAx toward zero based on shots faced, because a goalie with 200 shots faced should not get the same certainty as one with 1,800. That shrunk GSAx is divided by 5.15 goals per win, blended with goalie RAPM, and only then shifted above the 35th-percentile replacement goalie. This is intentionally different from the skater baseline: skaters use the 25th percentile within position group, while goalies use the 35th percentile among qualified goalies. The actual sequence is:

GSAx WAR        = shrunk GSAx / 5.15
goalie RAPM WAR = clamp((rapm_xga_60 − league_mean_rate) × TOI_hours / 5.15)

goalie WAR raw  = 0.85 × GSAx WAR + 0.15 × goalie RAPM WAR
goalie WAR      = goalie WAR raw − p35(goalie WAR raw)

The RAPM piece asks whether expected goals against changes when that goalie is in net after accounting for the skaters, teams, and home ice. The GSAx piece still does most of the work, but RAPM keeps the card from treating every environment as equally difficult since some teams and shooters are better than others. Replacement is applied once at the end so low-TOI backups do not inherit a full-season baseline bonus.

Choose a goalieISIlya Sorokin· NYI · G

Loading the full 2025–26 goalie list…

Ilya Sorokin

NYI · G · 55 games

4.33

Goalie WAR

100

WAR percentile

90

Composite

GSAx WAR+4.98 WAR

RAPM WAR+2.00 WAR

How the WAR is built

0.85 × +4.98 + 0.15 × +2.00 = 4.54

4.54 − 0.20 p35 shift = 4.33 WAR

Adj SV%

.907

Raw SV%

.906

GSAx

+28.0

xGA

172.0

Ilya Sorokin's goalie card starts with shot quality faced, shrinks GSAx for sample size, blends in goalie RAPM, and then asks how far above a replacement-level goalie that result sits.

This work is inspired by the excellent content at Evolving Hockey and HockeyStats, simplified to a single season. RAPM's regularized-regression lineage traces back to basketball's adjusted plus-minus work, and the shrinkage story is the same in any sport: when your signal is rare events distributed across shared credit, the honest move is to pull every estimate toward zero and let sustained evidence drag it back out. Celebrini's 50-second shift produced four attempts, 1.138 expected goals, and the winning geometry of Smith's one-foot finish. The model enjoyed the goal as much as anyone, but it doesn't treat this exceptional shift as the norm.