How the Expected Points Model Is Built (Part II)

Part I established what expected points is and showed the model’s output. This article is about how you get there — why the problem is shaped the way it is, why the obvious approaches fail, and how the shipped model works.

The through-line: predicting expected points is not a regression problem. It looks like one, because the answer is a number. It is not, and understanding why explains every design decision that follows.

The problem is multiclass, not continuous

The naive framing is “predict how many points the offense scores next.” Try to fit that directly and the problems appear immediately.

The target is not continuous. The next score in a half can only be one of seven things, and their values are −7, −3, −2, 0, 2, 3 or 7. There is no such thing as a 4.5-point outcome. A model that predicts 4.5 is not making a claim any football game can satisfy — it is averaging over outcomes it should be enumerating.

So the target is categorical, and expected points is a quantity you compute afterwards, from the class probabilities:

$\mathrm{EP} = \sum_{i=1}^{7} p_i \cdot v_i$

That single reframing is the whole conceptual step. Everything below is about estimating the $p_i$ well.

Why not linear regression

Linear regression assumes the dependent variable is continuous and unbounded, that the relationship between predictors and target is linear, that observations are independent, and that predictors are not collinear.

Point values violate the first assumption outright. But even setting that aside, predictions are unbounded — a linear model will happily tell you a situation is worth 9.2 points, which no single scoring event can produce. And the relationship is emphatically not linear: expected points changes far faster near either goal line than it does at midfield.

Why not binary logistic regression

Logistic regression fixes the boundedness problem. Model “does the offense score a touchdown next, yes or no”, and predictions live in $[0, 1]$ where probabilities belong.

But we need seven probabilities that are mutually exclusive and sum to one. Seven separately-fitted binary models give you seven numbers that do not sum to anything in particular. You cannot form a valid expectation from them.

Multinomial logistic regression

The classical fix is to fit the classes simultaneously against a reference class, and pass the results through a softmax so they are normalised to sum to one. Choose “No Score” as the reference, fit six sets of coefficients against it, and the seventh falls out.

This is exactly what the 2020 model did — and what Yurko, Ventura and Horowitz proposed for football in 2017 (Part III covers that lineage). It is a sound, interpretable and well-understood approach.

It is also where the 96 variables came from.

The cost of linearity: 96 variables

A multinomial logit is linear in its predictors. It can only represent an interaction if you build that interaction by hand, as a column. Football is almost entirely interactions — 3rd and 2 is not “3rd down plus 2 yards to go”, it is its own situation — so the feature list grew to cover them:

Group	Variables
Distance factors (yards to goal, log distance, goal-to-go, and their interaction)	24
Down, and down × distance, down × field position	54
Time factors, under-two-minute indicator, intercepts	18
Total	96

Each bullet is six variables because six classes are fitted against the reference. Every interaction anyone thought of had to be specified. Every interaction nobody thought of was invisible to the model.

There were two further consequences worth naming:

The model could not see the score. Game state was handled through observation weights — plays were weighted by score differential and by how many drives away the next score was — rather than as a feature.
Field goals needed a separate model, a GAM on smoothed kick distance, because the main model had no way to represent “this is makeable.”

What changed: gradient boosting

The current model is an XGBoost multi:softprob classifier. It answers the same seven-class question and produces probabilities the same way. What changes is how the structure is found.

Boosted trees split on feature values, so an interaction is representable without being specified: a tree that splits on down_3, then on distance, then on yards_to_goal is a three-way interaction, discovered from the data. The hand-built columns become unnecessary.

	2020	2026
Estimator	`nnet::multinom`	XGBoost `multi:softprob`
Inputs	96 constructed variables	8 raw features
Sees the score	no (weights only)	yes (`pos_score_diff_start`)
Training rows	2014–2019, non-overtime	2004–2025, 2,219,971 plays
Field goals	separate GAM	separate gradient-boosted `fg_model`

The eight features are simply:

c("TimeSecsRem", "yards_to_goal", "distance",
  "down_1", "down_2", "down_3", "down_4", "pos_score_diff_start")

Trained with eta = 0.025, max_depth = 5 and 525 boosting rounds — which, because there are seven classes, means the shipped booster contains 3,675 trees. That is a useful sanity check when you load it:

# xgb.dump() returns one element per NODE line, so length() is far larger than the
# tree count. Each tree is introduced by a `booster[i]` header -- count those.
sum(grepl("^booster\\[", xgboost::xgb.dump(ep_model)))   # 3675 = 525 rounds x 7 classes

What the model actually leans on

Gain-based importance, computed from the shipped booster:

Feature	Gain
`TimeSecsRem`	30.9%
`yards_to_goal`	21.9%
`down_4`	14.5%
`down_3`	9.4%
`down_1`	7.9%
`down_2`	6.9%
`pos_score_diff_start`	5.0%
`distance`	3.6%

imp <- xgboost::xgb.importance(model = ep_model)
imp[order(-imp$Gain), c("Feature", "Gain")]

Time is the largest single contributor, which is not obvious until you remember what the target is: the next score in this half. With thirty minutes left almost every drive has a successor; with thirty seconds left, “No Score” becomes the overwhelming favourite regardless of field position. The model spends most of its capacity on that.

Note also that down_4 carries more gain than the other three downs combined with each other — fourth down is where the possession itself is at risk, and the distribution shifts hardest.

A worked example, end to end

One situation, all the way through. First and 10 at your own 25, tied, thirty minutes left in the half:

x <- matrix(
  c(1800, 75, 10, 1, 0, 0, 0, 0), nrow = 1,
  dimnames = list(NULL, c("TimeSecsRem", "yards_to_goal", "distance",
                          "down_1", "down_2", "down_3", "down_4",
                          "pos_score_diff_start"))
)
# as.numeric(): predict() may return a 1 x 7 matrix here, and handing a matrix to
# data.frame() below would create seven columns against a seven-ROW outcome vector.
p <- as.numeric(predict(ep_model, x))
data.frame(
  outcome = c("TD", "Opp_TD", "FG", "Opp_FG", "Safety", "Opp_Safety", "No_Score"),
  prob    = round(p, 4),
  value   = c(7, -7, 3, -3, 2, -2, 0)
)

Outcome	Probability	Value	Contribution
Touchdown	0.3965	+7	+2.7753
Opponent touchdown	0.3182	−7	−2.2274
Field goal	0.1564	+3	+0.4691
Opponent field goal	0.1185	−3	−0.3555
Safety	0.0034	+2	+0.0068
Opponent safety	0.0024	−2	−0.0048
No score	0.0047	0	0.0000
			EP = +0.6635

Two thirds of a point, from a near-coin-flip between the two touchdown outcomes. That is what “expected points” means concretely: not a prediction that 0.66 points will be scored, but the average over seven futures, six of which are worth something.

The score-margin feature, and a caveat

Because the current model sees the score, EP at a fixed spot varies with the margin:

Score margin	EP at own 25
−21	−1.38
−7	−0.18
Tied	+0.66
+7	+1.71
+21	+2.71

Read this carefully. It does not mean leading causes you to score more from your own 25. It means teams that are up 21 are, on average, better teams playing worse opponents — the margin is partly a proxy for team quality, and the model has no separate team-strength input to separate the two.

This is a real limitation and worth stating plainly: EP is not opponent- or quality-adjusted. Any team-level metric built by summing EPA inherits that. Adjusted versions, like the ones on Game on Paper, apply opponent and garbage-time corrections on top of raw EPA precisely because the model itself does not.

One artifact, two languages

The largest structural change is not statistical. The models are published once as the cfb_model_artifacts release and read by both cfbfastR and sportsdataverse-py (#138). The same play gets the same EPA in R and in Python, and a retrain updates both by publishing, without a release of either package.

Making that true required pinning down one thing that is easy to get wrong. XGBoost emits classes in its training order, which is not the order cfbfastR has always reported:

manifest$ep_class_contract$class_order
#> "TD" "Opp_TD" "FG" "Opp_FG" "Safety" "Opp_Safety" "No_Score"

manifest$ep_class_contract$cfbfastR_lev_order
#> "No_Score" "FG" "Opp_FG" "Opp_Safety" "Opp_TD" "Safety" "TD"

manifest$ep_class_contract$permutation_to_cfbfastR_lev_1based
#> 7 3 4 6 2 5 1

The permutation is published in the manifest and read at runtime, not hard-coded. cfbfastR keeps a fallback copy, and a test asserts the two agree — so if the bundle is ever retrained with a different class order, the build fails loudly instead of silently reporting a touchdown probability as a safety probability. There is no fixed point between the two orderings, which is to say every position moves, which is to say a mistake here would be catastrophic and completely invisible in the output.

If you score the booster yourself, apply the point values in the bundle’s order, as in the worked example above.

Where the estimates are weakest

Being explicit about limits, since the point of this series is transparency:

No team-quality input. See the score-margin caveat above.
PATs are treated as given. The extra point is folded into the touchdown value of 7 rather than modelled.
Kickoffs assume a touchback baseline; returns past the 25 are points added and shorter returns points lost.
Punts are not modelled separately. Their effect is field position, which the EP model already handles well.
Overtime is not represented in the next-score-in-half framing.
Era effects are handled in the FG, two-point, QBR and fourth-down models through explicit era0–era3 features, but the EP model itself has no era term — it pools 2004 through 2025.

Part III traces where these ideas came from — Virgil Carter’s 1970 paper, the field position and down-distance models that followed, and the line that runs through nflscrapR and cfbscrapR to the shared bundle described here.

Data and artifacts

Models — cfb_model_artifacts
Season play-by-play — cfbfastR-data
Source — cfbfastR · sportsdataverse-py
Applied, live — Game on Paper
Raw API — CollegeFootballData.com, courtesy of @CFB_data

Citation

Gilani, S., Easwaran, A., Lee, J., and Hess, E. (2026). cfbfastR: Access College
Football Play by Play Data. R package version 3.0.0.9000.
https://cfbfastr.sportsdataverse.org

Authors, contributors and related SportsDataverse packages are listed on the package home page.

Saiem Gilani, Akshay Easwaran, Jared Lee

2026-09-03