Who Was the Best NFL Kicker, 2018–2020?
IS 470 Sports Analytics · Assignment 1
Luc Alliaume
October 08, 2026
Introduction
The question
This assignment asks a simple question: who was the best NFL kicker from 2018 to 2020? The obvious answer would be the kicker with the highest make percentage, but that measure is unfair. Kickers don’t get the same kicks. One kicker might spend a season attempting short field goals in a dome, while another attempts 50-yarders into the wind in December. The first will post a better percentage even if the second is the better kicker.
The approach
To compare kickers fairly, this project separates how hard a kick was from how well the kicker kicked it. It first builds a model of what an average NFL kicker would do on every kick, given its distance, angle, surface and weather. Each kicker is then judged only on how much better or worse he did than that average kicker would have on the same kicks.
The method has four steps:
- Kick difficulty. A logistic regression estimates the chance an average kicker makes each kick.
- Accuracy. Each kick’s aim error is compared with the aim error expected in those conditions.
- Points over expected. Each kicker is credited with the points he scored beyond what an average kicker would have scored on the same kicks.
- Standardize and rank. Points over expected is put on a per-kick basis, adjusted for small samples, and turned into z-scores so kickers can be ranked.
Key terms
- \(\hat{p}_i\) (expected make chance): the probability an average kicker makes kick \(i\). Its complement, \(1 - \hat{p}_i\), is the kick’s difficulty.
- Points over expected (POE): points a kicker scored above what an average kicker would have scored from the same kicks.
- Aim Error Over Expected (AEOE): how far off center a kicker’s kicks were, compared with what’s expected in those conditions. Lower is better.
- z-score: how many standard deviations a kicker sits above or below the average qualifying kicker.
The data
Sources
The dataset combines three sources:
- NFL Big Data Bowl 2022: play results, games, players, PFF scouting data and ball tracking for every special-teams play from 2018 to 2020. Tracking gives each kick’s hold spot and where the ball crossed the goal line.
- NFLWeather.com: temperature, wind speed and direction, precipitation and stadium details (surface, roof and compass orientation) for each quarter of each game.
- Three hand-filled lookup files: each kicker’s kicking foot; each stadium’s surface, roof and elevation; and which compass end of each stadium a reference kick was aimed at, found on game video. The last file links the tracking coordinates to real compass directions, so wind can be split into crosswind and head/tailwind.
Building the dataset
The build script keeps every real field goal and extra point attempt, removing blocked kicks, fakes and bad snaps. From the tracking data it measures distance, the hold spot’s sideways position and where the ball crossed the goal line. When tracking stopped just short of the goal line, the ball’s last straight path is extended (these kicks are flagged). Weather is matched to the quarter of each kick, and indoor games have their wind, rain and snow set to zero and their temperature set to 70°F.
Several variables are derived from these measurements:
- Kicking direction \(\varphi_i = \arctan(\tilde{L}_i / D_i)\), where \(\tilde{L}_i\) is the hold spot’s sideways offset, signed toward the kicking-foot side so left- and right-footed kickers compare fairly.
- Aim error \(e_i\), the yards between where the ball crossed and the center of the uprights.
- Crosswind (wind across the kick’s path) and tailwind (wind along it, positive when it helps).
- Cold (degrees below 50°F) and thin air (the drop in air density from elevation and temperature).
This code is shown for reference but isn’t re-run when the page is built, since it reads several gigabytes of tracking data and downloads about 768 web pages.
# =============================================================================
# build_kick_dataset.R
# Builds one row per kick (field goals + extra points, 2018-2020) by combining:
# 1. NFL Big Data Bowl 2022 data (plays, games, players, PFF scouting, tracking)
# 2. Game-by-quarter weather scraped from NFLWeather.com
# 3. Three small hand-filled lookup files (kicker footedness, stadium details,
# and which compass end of each stadium a kick was aimed at)
#
# Forecasts provided by NFLWeather.com (https://www.nflweather.com/).
# Their terms of use require this attribution next to any republished data.
#
# Run from the same working directory as your Rmd (the folder that contains
# NFLBDB2022/). The first run downloads ~820 pages from NFLWeather.com at a
# polite pace (about 20-25 minutes); every page is cached to disk, so later
# runs are fast and never hit the site again.
# =============================================================================
library(data.table)
library(dplyr)
library(stringr)
library(rvest)
# ---- Settings ---------------------------------------------------------------
DATA_DIR <- "nfl-big-data-bowl-2022"
CACHE_DIR <- "nflweather_cache" # downloaded pages are saved here
SEASONS <- 2018:2020
WEEKS <- 1:17 # BDB 2022 covers regular season weeks 1-17
PAUSE_SEC <- 1.5 # wait between downloads (be polite)
ASSUME_RETRACTABLE_CLOSED <- TRUE # roof status per game isn't on the site
WIND_DIR_IS_FROM <- TRUE # TRUE if NFLWeather's wind direction is where the
# wind comes FROM (standard). See the instructions
# doc for the one-time check.
POST_LEFT <- 23.58 # goal-post y positions used in Demo.Rmd
POST_RIGHT <- 29.75
POST_CENTER <- (POST_LEFT + POST_RIGHT) / 2
UPRIGHT_W <- 6.17 # yards between uprights (18 ft 6 in)
HOLD_DEPTH <- 7 # ball is kicked ~7 yards behind the line
MAX_EXTRAP_YDS <- 10 # extend a kick's path to the end line only if
# tracking stopped within this many yards of it
dir.create(CACHE_DIR, showWarnings = FALSE)
# =============================================================================
# PART 1. Kicks from the Big Data Bowl files
# =============================================================================
build_kicks <- function() {
plays <- fread(file.path(DATA_DIR, "plays.csv"))
games <- fread(file.path(DATA_DIR, "games.csv"))
players <- fread(file.path(DATA_DIR, "players.csv"))
pff <- fread(file.path(DATA_DIR, "PFFScoutingData.csv"))
# Check the bad-snap codes once with: table(pff$snapDetail)
kicks <- plays |>
filter(specialTeamsPlayType %in% c("Field Goal", "Extra Point"),
# keeps only real attempts: drops blocked kicks and fakes
specialTeamsResult %in% c("Kick Attempt Good", "Kick Attempt No Good")) |>
left_join(pff |> select(gameId, playId, snapDetail), by = c("gameId", "playId")) |>
filter(is.na(snapDetail) | snapDetail == "OK") |> # drop bad snaps
left_join(games |> select(gameId, season, week, gameDate,
homeTeamAbbr, visitorTeamAbbr), by = "gameId") |>
filter(season %in% SEASONS, week %in% WEEKS) |>
left_join(players |> select(nflId, kicker = displayName),
by = c("kickerId" = "nflId")) |>
mutate(is_xp = as.integer(specialTeamsPlayType == "Extra Point"),
made = as.integer(specialTeamsResult == "Kick Attempt Good"),
points = ifelse(is_xp == 1, 1, 3)) |>
select(gameId, playId, season, week, gameDate, quarter, gameClock,
homeTeamAbbr, visitorTeamAbbr, possessionTeam, playDescription,
kickerId, kicker, is_xp, made, points, kickLength)
as.data.table(kicks)
}
# Ball position at the snap and where it crossed the end line, per kick.
# Reads one tracking file at a time and keeps only the ball on kick plays,
# so memory stays manageable.
build_ball_geometry <- function(kicks) {
ids <- unique(kicks[, .(gameId = as.numeric(gameId), playId = as.numeric(playId))])
read_ball <- function(f) {
tr <- fread(f, select = c("gameId", "playId", "frameId", "team",
"x", "y", "event", "playDirection"))
tr[, `:=`(gameId = as.numeric(gameId), playId = as.numeric(playId))]
tr <- tr[team == "football"][ids, on = .(gameId, playId), nomatch = 0]
# Keep the ORIGINAL direction before flipping. "left" = kick heads toward
# the low-x end (x = 0) of the real stadium; "right" = toward x = 120.
tr[, play_dir_raw := playDirection]
# Same flip as Demo.Rmd: every kick now travels toward x = 0, and
# larger y is the kicker's right. Used only for the left/right geometry.
tr[playDirection == "right", `:=`(x = 120 - x, y = 160 / 3 - y)]
tr
}
files <- file.path(DATA_DIR, paste0("tracking", SEASONS, ".csv"))
ball <- rbindlist(lapply(files, read_ball))
setorder(ball, gameId, playId, frameId)
ball[, {
fs <- which(event == "ball_snap")[1]
if (is.na(fs)) fs <- 1L
ic <- which(seq_along(x) > fs & x < 0)[1] # first frame past end line
y_cross <- NA_real_
if (!is.na(ic) && x[ic] != x[ic - 1]) { # interpolate y at x = 0
t <- (0 - x[ic - 1]) / (x[ic] - x[ic - 1])
y_cross <- y[ic - 1] + t * (y[ic] - y[ic - 1])
}
# Tracking sometimes stops just before the ball reaches the end line
# (common on short made kicks). Extend the ball's last straight-line
# path to x = 0, using up to 5 frames of it moving toward the posts.
e_extrap <- 0L
if (is.na(y_cross)) {
fl <- which(seq_along(x) > fs & x < x[fs] - 1 & c(FALSE, diff(x) < 0))
fl <- tail(fl, 5)
if (length(fl) >= 2) {
a <- fl[1]; b <- fl[length(fl)]
if (x[b] <= MAX_EXTRAP_YDS && x[a] - x[b] > 0.5) {
slope <- (y[b] - y[a]) / (x[b] - x[a])
y_cross <- y[b] + slope * (0 - x[b])
e_extrap <- 1L
}
}
}
list(x_snap = x[fs], y_snap = y[fs], y_cross = y_cross,
e_extrapolated = e_extrap, play_dir_raw = play_dir_raw[1])
}, by = .(gameId, playId)]
}
# =============================================================================
# PART 2. Weather from NFLWeather.com
# =============================================================================
cached_page <- function(url) {
file <- file.path(CACHE_DIR, paste0(gsub("[^A-Za-z0-9]+", "_",
sub("^https?://", "", url)), ".html"))
if (!file.exists(file)) {
Sys.sleep(PAUSE_SEC)
ok <- tryCatch({ download.file(url, file, quiet = TRUE); TRUE },
error = function(e) { message("FAILED: ", url); FALSE })
if (!ok) return(NULL)
}
read_html(file)
}
week_game_links <- function(season, week) {
url <- sprintf("https://www.nflweather.com/week/%d/week-%d", season, week)
page <- cached_page(url)
if (is.null(page)) return(character(0))
hrefs <- html_attr(html_elements(page, "a"), "href")
hrefs <- unique(hrefs[str_detect(hrefs, sprintf("/games/%d/week-%d/", season, week))])
url_absolute(hrefs, "https://www.nflweather.com")
}
num <- function(txt, pattern) as.numeric(str_match(txt, pattern)[, 2])
# One quarter's block of text -> one row of weather
parse_block <- function(b) {
data.table(
temp = num(b, "(-?\\d+)\\s*°F"), # first °F = air temp
feels_like = num(b, "Feels Like:\\s*(-?\\d+)"),
wind_mph = num(b, "(\\d+)\\s*mph"),
wind_dir = str_match(b, "mph\\s+(?:[a-z_]+\\s+)?([NSEW]{1,3})\\b")[, 2],
precip_prob = num(b, "Prec\\. Prob\\.:\\s*(\\d+)"),
cloud = num(b, "Cloud Cover:\\s*(\\d+)"),
humidity = num(b, "Humidity:\\s*(\\d+)"),
dew_point = num(b, "Dew Point:\\s*(-?\\d+)"),
visibility = num(b, "Visibility:\\s*([\\d.]+)"),
rain = as.integer(str_detect(b, regex("rain|drizzle|shower", ignore_case = TRUE))),
snow = as.integer(str_detect(b, regex("snow|flurr|sleet", ignore_case = TRUE)))
)
}
# Full game page text -> 4 rows (Kickoff, Q2, Q3, Q4) plus stadium details
parse_game_text <- function(txt, url = NA_character_, stadium_slug = NA_character_) {
heads <- c("Kickoff", "Q2", "Q3", "Q4")
starts <- sapply(heads, function(h)
str_locate(txt, regex(paste0("^\\W*", h, "\\s*$"), multiline = TRUE))[1, "start"])
if (anyNA(starts)) return(NULL)
stad_at <- str_locate(txt, "Surface:")[1, "start"]
if (is.na(stad_at)) stad_at <- nchar(txt) + 1
# Quarters are not always listed in order on the page (e.g. Kickoff, Q3,
# Q4, Q2), so each block ends where the next block on the page begins.
nxt <- sapply(starts, function(st) min(c(starts[starts > st], stad_at)))
ends <- nxt - 1
wx <- rbindlist(lapply(seq_along(heads), function(k)
parse_block(substr(txt, starts[k], ends[k]))))
wx[, wx_period := heads]
orient <- str_match(txt, "Orientation:\\s*(?:[a-z_]+\\s*)?([NSEW]{1,2})\\s*[–-]\\s*([NSEW]{1,2})")
wx[, `:=`(
url = url,
game_time = str_match(txt, "(\\d{2}/\\d{2}/\\d{2} \\d{2}:\\d{2} [AP]M)")[, 2],
stadium = stadium_slug,
surface = str_trim(str_match(txt, "Surface:\\s*(.*?)\\s*(?:house\\s*)?Location:")[, 2]),
roof = str_trim(str_match(txt, "Type:\\s*(.*?)\\s*(?:explore\\s*)?Orientation:")[, 2]),
orient_end1 = orient[, 2],
orient_end2 = orient[, 3]
)]
wx[]
}
parse_game_page <- function(url) {
page <- cached_page(url)
if (is.null(page)) return(NULL)
hrefs <- html_attr(html_elements(page, "a"), "href")
stad <- str_match(hrefs[str_detect(hrefs, "/stadium/")][1], "/stadium/([^/?#]+)")[, 2]
parse_game_text(html_text2(page), url, stad)
}
build_weather <- function() {
grid <- expand.grid(season = SEASONS, week = WEEKS)
links <- rbindlist(lapply(seq_len(nrow(grid)), function(i) {
l <- week_game_links(grid$season[i], grid$week[i])
data.table(season = grid$season[i], week = grid$week[i], url = l)
}))
message(nrow(links), " game pages found")
wx <- rbindlist(lapply(links$url, parse_game_page), fill = TRUE)
# Some quarters show wind as "TBD". Borrow the nearest quarter of the same
# game (next quarter first, then the previous one) and flag it.
wx[, wind_filled := as.integer(is.na(wind_mph))]
fill_nearest <- function(v) {
ok <- which(!is.na(v))
if (length(ok) == 0) return(v)
v[sapply(seq_along(v), function(k) {
if (!is.na(v[k])) return(k)
nxt <- ok[ok > k]
if (length(nxt)) nxt[1] else tail(ok, 1)
})]
}
wx[, `:=`(wind_mph = fill_nearest(wind_mph), wind_dir = fill_nearest(wind_dir)), by = url]
links[wx, on = "url"]
}
# =============================================================================
# PART 3. Matching keys and lookups
# =============================================================================
# BDB team abbreviation -> nickname used in NFLWeather URLs
abbr_to_nick <- c(
ARI = "cardinals", ATL = "falcons", BAL = "ravens", BUF = "bills",
CAR = "panthers", CHI = "bears", CIN = "bengals", CLE = "browns",
DAL = "cowboys", DEN = "broncos", DET = "lions", GB = "packers",
HOU = "texans", IND = "colts", JAX = "jaguars", JAC = "jaguars",
KC = "chiefs", LA = "rams", LAR = "rams", LAC = "chargers",
LV = "raiders", OAK = "raiders", MIA = "dolphins", MIN = "vikings",
NE = "patriots", NO = "saints", NYG = "giants", NYJ = "jets",
PHI = "eagles", PIT = "steelers", SEA = "seahawks", SF = "49ers",
TB = "buccaneers", TEN = "titans", WAS = "washington"
)
home_nick_from_url <- function(url) {
# Decode first: Washington's 2020 pages use "football%20team"
nick <- tolower(str_match(URLdecode(url), "-at-(.+?)/?$")[, 2])
nick <- str_replace_all(str_trim(nick), "\\s+", "-")
# Washington's URLs use different names by season
ifelse(str_detect(nick, "redskins|washington|football|commanders"), "washington", nick)
}
compass_deg <- c(N = 0, NNE = 22.5, NE = 45, ENE = 67.5, E = 90, ESE = 112.5,
SE = 135, SSE = 157.5, S = 180, SSW = 202.5, SW = 225,
WSW = 247.5, W = 270, WNW = 292.5, NW = 315, NNW = 337.5)
opposite_end <- c(N = "S", NE = "SW", E = "W", SE = "NW",
S = "N", SW = "NE", W = "E", NW = "SE")
# Keeps a hand-filled CSV in sync: writes it the first time; on later runs
# adds rows for any new keys (blank to fill) and never touches your edits.
sync_template <- function(f, fresh, key) {
fresh <- fresh[, lapply(.SD, as.character)]
if (!file.exists(f)) {
fwrite(fresh, f)
message("Wrote ", f, " (", nrow(fresh), " rows) - fill it in, then rerun.")
return(fresh)
}
old <- fread(f, colClasses = "character", na.strings = c("", "NA"))
add <- fresh[!fresh[[key]] %in% old[[key]]]
if (nrow(add) > 0) {
out <- rbind(old, add, fill = TRUE)
fwrite(out, f)
message("Added ", nrow(add), " new rows to ", f, " - fill them in, then rerun.")
return(out)
}
old
}
# Hand-filled file 1: kicker footedness ("R" is pre-filled; change lefties to "L")
kicker_foot_lookup <- function(kicks) {
fresh <- unique(kicks[, .(kickerId, kicker)])[order(kicker)][, foot := "R"]
ft <- sync_template("kicker_foot.csv", fresh, "kickerId")
ft[, .(kickerId = as.numeric(kickerId), foot = toupper(str_trim(foot)))]
}
# Hand-filled file 2: stadium details (fill elevation_ft; check grass, indoor)
stadium_lookup <- function(weather) {
st <- unique(weather[!is.na(stadium), .(stadium, surface, roof)], by = "stadium")
st[, `:=`(
grass = as.integer(str_detect(surface, regex("grass", ignore_case = TRUE)) &
!str_detect(surface, regex("turf|artificial|synthetic|realgrass",
ignore_case = TRUE))),
indoor = as.integer(str_detect(roof, regex("dome|indoor|closed|fixed", ignore_case = TRUE)) |
(ASSUME_RETRACTABLE_CLOSED &
str_detect(roof, regex("retract", ignore_case = TRUE)))),
elevation_ft = NA)]
st <- sync_template("stadiums.csv", st[order(stadium)], "stadium")
st[, .(stadium, grass = as.integer(grass), indoor = as.integer(indoor),
elevation_ft = as.numeric(gsub(",", "", elevation_ft)))]
}
# Hand-filled file 3: which compass end each stadium's low-x end faces.
# For each stadium the script picks two field goals from two different games
# and says which coordinate end each kick went toward. You look the kick up on
# video, note which compass end of the stadium it was kicked toward, and type
# that end (one of the two letters in `orientation`) into a1_kicked_toward and
# a2_kicked_toward. The script works out the rest.
stadium_ends_lookup <- function(d, stad) {
fg <- d[is_xp == 0 & !is.na(stadium) & !is.na(play_dir_raw)]
fg <- fg[order(stadium, gameId, -kickLength)]
fg <- fg[, .SD[1], by = .(stadium, gameId)] # best kick per game
fg <- fg[, head(.SD, 2), by = stadium] # two games per stadium
fg[, n := seq_len(.N), by = stadium]
fg[, desc := sprintf("%s: %s @ %s, Q%s %s - %s", gameDate, visitorTeamAbbr,
homeTeamAbbr, quarter, gameClock, playDescription)]
fg[, toward := ifelse(play_dir_raw == "left", "LOW x end (x = 0)", "HIGH x end (x = 120)")]
wide <- unique(d[!is.na(stadium), .(stadium, orientation = paste0(orient_end1, "-", orient_end2))],
by = "stadium")
for (k in 1:2) {
a <- fg[n == k, .(stadium, gameId, playId, desc, toward)]
setnames(a, c("gameId", "playId", "desc", "toward"),
paste0("a", k, c("_gameId", "_playId", "_kick", "_coord_end")))
wide <- a[wide, on = "stadium"]
}
wide <- stad[, .(stadium, indoor)][wide, on = "stadium"]
wide[, needed := ifelse(indoor == 1, "skip (indoor)", "FILL IN")]
wide[, `:=`(a1_kicked_toward = NA, a2_kicked_toward = NA)]
setcolorder(wide, c("stadium", "needed", "orientation",
"a1_kick", "a1_coord_end", "a1_kicked_toward",
"a2_kick", "a2_coord_end", "a2_kicked_toward",
"a1_gameId", "a1_playId", "a2_gameId", "a2_playId"))
wide[, indoor := NULL]
en <- sync_template("stadium_ends.csv", wide[order(needed, stadium)], "stadium")
# Turn each observation into "which compass end is the low-x end"
low_from <- function(obs, coord_end) {
obs <- toupper(str_trim(obs))
ifelse(is.na(obs) | is.na(coord_end), NA_character_,
ifelse(str_detect(coord_end, "LOW"), obs, opposite_end[obs]))
}
en[, `:=`(low1 = low_from(a1_kicked_toward, a1_coord_end),
low2 = low_from(a2_kicked_toward, a2_coord_end))]
bad <- en[!is.na(low1) & !is.na(low2) & low1 != low2]
if (nrow(bad) > 0)
warning("Anchor kicks disagree at: ", paste(bad$stadium, collapse = ", "),
". Recheck those videos (or the stadium uses different coordinates by game).")
ends <- strsplit(en$orientation, "-")
ok <- mapply(function(l, e) is.na(l) || l %in% e, fcoalesce(en$low1, en$low2), ends)
if (any(!ok))
warning("Entries that aren't one of the two orientation ends at: ",
paste(en$stadium[!ok], collapse = ", "))
en[, low_x_end := fifelse(ok, fcoalesce(low1, low2), NA_character_)]
en[, .(stadium, low_x_end)]
}
# Air density relative to sea level at 59 F (0 = normal air, larger = thinner)
thin_air <- function(elev_ft, temp_f) {
h <- elev_ft * 0.3048
P <- 101325 * (1 - 2.25577e-5 * h)^5.25588
rho <- P / (287.05 * ((temp_f - 32) * 5 / 9 + 273.15))
1 - rho / (101325 / (287.05 * 288.15))
}
# =============================================================================
# PART 4. Put it all together
# =============================================================================
build_kick_dataset <- function() {
kicks <- build_kicks()
geom <- build_ball_geometry(kicks)
weather <- build_weather()
foot <- kicker_foot_lookup(kicks)
stad <- stadium_lookup(weather)
kicks[, `:=`(gameId = as.numeric(gameId), playId = as.numeric(playId),
kickerId = as.numeric(kickerId))]
weather[, home_nick := home_nick_from_url(url)]
kicks[, `:=`(home_nick = abbr_to_nick[homeTeamAbbr],
wx_period = c("Kickoff", "Q2", "Q3", "Q4", "Q4")[pmin(quarter, 5)])]
d <- geom[kicks, on = .(gameId, playId)]
d <- weather[d, on = .(season, week, home_nick, wx_period)]
d <- stad[d, on = "stadium"]
d <- foot[d, on = "kickerId"]
ends <- stadium_ends_lookup(d, stad)
d <- ends[d, on = "stadium"]
d[, `:=`(
# Step 1: kick geometry
D = fifelse(!is.na(kickLength), as.numeric(kickLength), round(x_snap + HOLD_DEPTH)),
L_raw = y_snap - POST_CENTER # + = kicker's right
)]
d[, `:=`(
L_tilde = fifelse(foot == "L", -L_raw, L_raw), # + = kicking-foot side
e = y_cross - POST_CENTER # Step 2: aim error (yards)
)]
d[, `:=`(
theta = (atan((L_tilde + UPRIGHT_W / 2) / D) - atan((L_tilde - UPRIGHT_W / 2) / D)) * 180 / pi,
phi = atan(L_tilde / D) * 180 / pi,
alpha = atan(abs(e) / D) * 180 / pi
)]
# Step 1: weather, with indoor games zeroed out
axis <- compass_deg[d$orient_end1]
wdir <- compass_deg[d$wind_dir]
wind_to <- if (WIND_DIR_IS_FROM) (wdir + 180) %% 360 else wdir # where the air is going
low_deg <- compass_deg[d$low_x_end]
kick_deg <- ifelse(d$play_dir_raw == "left", low_deg, (low_deg + 180) %% 360)
d[, `:=`(
crosswind = fifelse(indoor == 1, 0, wind_mph * abs(sin((wdir - axis) * pi / 180))),
along_wind = fifelse(indoor == 1, 0, wind_mph * abs(cos((wdir - axis) * pi / 180))),
# + = tailwind (helps distance), - = headwind
tailwind = fifelse(indoor == 1, 0, wind_mph * cos((wind_to - kick_deg) * pi / 180)),
wind_mph = fifelse(indoor == 1, 0, wind_mph),
temp_used = fifelse(indoor == 1, 70, temp),
rain = fifelse(indoor == 1, 0L, rain),
snow = fifelse(indoor == 1, 0L, snow)
)]
d[, `:=`(cold = pmax(0, 50 - temp_used),
thin_air = thin_air(fcoalesce(elevation_ft, 0.0), temp_used))]
# Quick checks
message(sum(is.na(d$temp)), " kicks with no weather match (check home_nick / week)")
message(sum(is.na(d$y_snap)), " kicks with no tracking")
message(sum(is.na(d$foot)), " kicks with no footedness")
message(sum(is.na(d$elevation_ft)), " kicks at stadiums with no elevation (treated as sea level)")
message(sum(d$e_extrapolated == 1, na.rm = TRUE), " kicks with e extended to the end line")
message(sum(d$wind_filled == 1, na.rm = TRUE), " kicks with wind borrowed from a nearby quarter")
message(sum(is.na(d$tailwind)), " kicks with no head/tailwind (stadium_ends.csv not filled)")
out <- d[, .(gameId, playId, season, week, quarter, kickerId, kicker, foot,
is_xp, made, points, D, L_raw, L_tilde, theta, phi, e, e_extrapolated, alpha,
grass, indoor, stadium, wind_mph, wind_dir, wind_filled, crosswind, along_wind, tailwind,
play_dir_raw, low_x_end,
temp = temp_used, feels_like, cold, rain, snow, precip_prob,
humidity, dew_point, visibility, cloud, elevation_ft, thin_air)]
setorder(out, gameId, playId)
fwrite(out, "kick_dataset.csv")
message("Saved kick_dataset.csv with ", nrow(out), " kicks")
out
}
if (!exists("SKIP_MAIN")) kick_df <- build_kick_dataset()Creating one model-ready file
A second script keeps the columns the model needs, grouped by role, and adds three helpers: each kick’s point value, each kicker’s number of attempts and a flag for kickers with 30+ attempts.
# =============================================================================
# make_model_data.R
# IS 470 Assignment 1 - Best NFL kicker, 2018-2020
#
# Turns kick_dataset.csv (made by build_kick_dataset.R) into ONE tidy,
# model-ready file, with columns grouped in the order the model uses them:
#
# 1. Identifiers - which kick, which kicker
# 2. Outcomes - what happened (never used as predictors)
# 3. Difficulty (G) - kick geometry and setup
# 4. Weather (E) - conditions and environment
# 5. Data flags - where values were filled in
# 6. Kicker totals - attempts and the 30+ attempt qualifier
#
# Also writes variable_dictionary.csv describing every column.
#
# Run build_kick_dataset.R first, then this, from the same folder.
# Outputs: kick_model_data.csv, variable_dictionary.csv
# =============================================================================
library(data.table)
INPUT <- "kick_dataset.csv"
OUTPUT <- "kick_model_data.csv"
DICTIONARY <- "variable_dictionary.csv"
MIN_ATTEMPTS <- 30 # kickers with at least this many attempts get ranked
if (!file.exists(INPUT))
stop("Can't find ", INPUT, ". Run build_kick_dataset.R first, from this folder.")
d <- fread(INPUT)
# ---- Small additions --------------------------------------------------------
d[, v := fifelse(is_xp == 1, 1L, 3L)] # Step 3 point weight
d[, n_attempts := .N, by = kickerId] # all kicks, FG + XP
d[, qualifies := as.integer(n_attempts >= MIN_ATTEMPTS)] # ranked in Step 4
# ---- Column groups (order = order in the output file) -----------------------
dict <- rbindlist(list(
# 1. Identifiers
list("gameId", "Identifier", "", "Game ID (Big Data Bowl)"),
list("playId", "Identifier", "", "Play ID within the game"),
list("season", "Identifier", "", "Season (2018-2020)"),
list("week", "Identifier", "", "Week (1-17)"),
list("quarter", "Identifier", "", "Quarter (5 = overtime)"),
list("kickerId", "Identifier", "k", "Kicker ID"),
list("kicker", "Identifier", "", "Kicker name"),
list("stadium", "Identifier", "", "Stadium"),
# 2. Outcomes
list("made", "Outcome", "y_i", "1 = kick good, 0 = no good (Step 1 response)"),
list("v", "Outcome", "v_i", "Point weight: 3 for a field goal, 1 for an extra point (Step 3)"),
list("e", "Outcome", "e_i", "Aim error at the goal line, yards (+ = right of center)"),
list("alpha", "Outcome", "alpha_i", "Aim-error angle, degrees = arctan(|e|/D) (Step 2)"),
# 3. Difficulty
list("D", "Difficulty", "D_i", "Kick distance, yards; enters the model as a natural spline f(D)"),
list("theta", "Difficulty", "theta_i", "Kicking window: angle the uprights fill, degrees (nearly determined by D)"),
list("phi", "Difficulty", "phi_i", "Kicking direction: angle off center, degrees, + toward the kicking-foot side"),
list("grass", "Difficulty", "S_i", "1 = natural grass, 0 = artificial turf"),
list("is_xp", "Difficulty", "X_i", "1 = extra point, 0 = field goal"),
list("L_tilde", "Difficulty", "L~_i", "Hold spot's sideways offset from center, yards, + toward the kicking-foot side"),
list("foot", "Difficulty", "", "Kicking foot (L/R); sets the sign of L~ and phi"),
# 4. Weather / environment
list("crosswind", "Weather", "W_i", "Wind across the kick's path, mph (0 indoors)"),
list("tailwind", "Weather", "", "Wind along the kick's path, mph: + helps (tailwind), - hurts (headwind)"),
list("cold", "Weather", "C_i", "Degrees below 50 F, 0 if warmer"),
list("rain", "Weather", "R_i", "1 = rain (0 indoors)"),
list("snow", "Weather", "N_i", "1 = snow (0 indoors)"),
list("thin_air", "Weather", "A_i", "Drop in air density vs. sea level at 59 F (elevation + temperature)"),
list("precip_prob", "Weather", "", "Chance of precipitation, % (candidate)"),
list("temp", "Weather", "", "Temperature, F (70 indoors); source of cold"),
list("feels_like", "Weather", "", "Feels-like temperature, F (near-duplicate of temp)"),
list("humidity", "Weather", "", "Relative humidity, % (candidate)"),
list("dew_point", "Weather", "", "Dew point, F (candidate)"),
list("visibility", "Weather", "", "Visibility (candidate)"),
list("cloud", "Weather", "", "Cloud cover, % (candidate)"),
list("indoor", "Weather", "", "1 = dome or closed roof (wind, rain, snow already zeroed)"),
# 5. Data flags
list("e_extrapolated","Flag", "", "1 = ball tracking ended early; e extended to the end line"),
list("wind_filled", "Flag", "", "1 = wind missing for that quarter; borrowed from the nearest quarter"),
# 6. Kicker totals
list("n_attempts", "Kicker", "n_k", "Kicker's total attempts in the data (FG + XP)"),
list("qualifies", "Kicker", "", paste0("1 = kicker has ", MIN_ATTEMPTS, "+ attempts and is ranked in Step 4"))
))
setnames(dict, c("column", "group", "symbol", "description"))
missing_cols <- setdiff(dict$column, names(d))
if (length(missing_cols) > 0)
stop("kick_dataset.csv is missing: ", paste(missing_cols, collapse = ", "),
". Rerun the latest build_kick_dataset.R.")
out <- d[, dict$column, with = FALSE]
setorder(out, season, week, gameId, playId)
# ---- Quick checks ------------------------------------------------------------
core <- c("made", "D", "phi", "grass", "is_xp", "crosswind", "tailwind",
"cold", "rain", "snow", "thin_air")
na_counts <- sapply(out[, ..core], function(x) sum(is.na(x)))
message(nrow(out), " kicks (", sum(out$is_xp == 0), " field goals, ",
sum(out$is_xp == 1), " extra points)")
message(uniqueN(out$kickerId), " kickers, ",
uniqueN(out[qualifies == 1, kickerId]), " with ", MIN_ATTEMPTS, "+ attempts")
if (any(na_counts > 0)) {
message("Missing values in model columns:")
print(na_counts[na_counts > 0])
} else {
message("No missing values in the core model columns")
}
fwrite(out, OUTPUT)
fwrite(dict, DICTIONARY)
message("Saved ", OUTPUT, " (", ncol(out), " columns) and ", DICTIONARY)The finished dataset
d <- fread("kick_model_data.csv")
data.table(
Item = c("Kicks", "Field goals", "Extra points", "Misses",
"Kickers", "Kickers with 30+ attempts", "Distance range (yards)"),
Value = c(nrow(d), sum(d$is_xp == 0), sum(d$is_xp == 1), sum(d$made == 0),
uniqueN(d$kickerId), uniqueN(d[qualifies == 1, kickerId]),
paste(range(d$D), collapse = "–"))
) |> kable(align = "lr")| Item | Value |
|---|---|
| Kicks | 6055 |
| Field goals | 2604 |
| Extra points | 3451 |
| Misses | 585 |
| Kickers | 60 |
| Kickers with 30+ attempts | 45 |
| Distance range (yards) | 19–67 |
if (file.exists("variable_dictionary.csv")) {
dict <- fread("variable_dictionary.csv")
kable(dict[group %in% c("Outcome", "Difficulty", "Weather")],
col.names = c("Column", "Group", "Symbol", "Description"))
}| Column | Group | Symbol | Description |
|---|---|---|---|
| made | Outcome | y_i | 1 = kick good, 0 = no good (Step 1 response) |
| v | Outcome | v_i | Point weight: 3 for a field goal, 1 for an extra point (Step 3) |
| e | Outcome | e_i | Aim error at the goal line, yards (+ = right of center) |
| alpha | Outcome | alpha_i | Aim-error angle, degrees = arctan(|e|/D) (Step 2) |
| D | Difficulty | D_i | Kick distance, yards; enters the model as a natural spline f(D) |
| theta | Difficulty | theta_i | Kicking window: angle the uprights fill, degrees (nearly determined by D) |
| phi | Difficulty | phi_i | Kicking direction: angle off center, degrees, + toward the kicking-foot side |
| grass | Difficulty | S_i | 1 = natural grass, 0 = artificial turf |
| is_xp | Difficulty | X_i | 1 = extra point, 0 = field goal |
| L_tilde | Difficulty | L~_i | Hold spot’s sideways offset from center, yards, + toward the kicking-foot side |
| foot | Difficulty | Kicking foot (L/R); sets the sign of L~ and phi | |
| crosswind | Weather | W_i | Wind across the kick’s path, mph (0 indoors) |
| tailwind | Weather | Wind along the kick’s path, mph: + helps (tailwind), - hurts (headwind) | |
| cold | Weather | C_i | Degrees below 50 F, 0 if warmer |
| rain | Weather | R_i | 1 = rain (0 indoors) |
| snow | Weather | N_i | 1 = snow (0 indoors) |
| thin_air | Weather | A_i | Drop in air density vs. sea level at 59 F (elevation + temperature) |
| precip_prob | Weather | Chance of precipitation, % (candidate) | |
| temp | Weather | Temperature, F (70 indoors); source of cold | |
| feels_like | Weather | Feels-like temperature, F (near-duplicate of temp) | |
| humidity | Weather | Relative humidity, % (candidate) | |
| dew_point | Weather | Dew point, F (candidate) | |
| visibility | Weather | Visibility (candidate) | |
| cloud | Weather | Cloud cover, % (candidate) | |
| indoor | Weather | 1 = dome or closed roof (wind, rain, snow already zeroed) |
Interpreting the data
Extra points make up 57% of all kicks, and only 9.7% of kicks were missed. Most kicks are close to automatic, so the differences between kickers come from a fairly small number of hard kicks and rare misses. That is why the method weighs each kick by its difficulty instead of simply counting makes: a missed extra point and a missed 55-yarder are very different events, and make percentage treats them the same.
# Model settings
KNOTS <- c(25, 33, 40, 50, 55) # inner knots for the distance spline
ALT_KNOTS <- list("30, 40, 50" = c(30, 40, 50),
"30, 36, 42, 48, 55" = c(30, 36, 42, 48, 55))
SHRINK <- 25 # POE rate = POE / (n + SHRINK)
MIN_ATTEMPTS <- 30 # kickers ranked in Step 4
G_VARS <- c("phi", "grass", "is_xp") # kick difficulty
E_VARS <- c("crosswind", "tailwind", "cold", "rain", "snow", "thin_air") # weather
d <- d[complete.cases(d[, c("made", "D", G_VARS, E_VARS), with = FALSE])]
BOUNDARY <- range(d$D) # outer spline knots at the shortest and longest kick
make_formula <- function(response, knots) {
as.formula(paste0(
response, " ~ ns(D, knots = c(", paste(knots, collapse = ", "),
"), Boundary.knots = c(", BOUNDARY[1], ", ", BOUNDARY[2], ")) + ",
paste(c(G_VARS, E_VARS), collapse = " + ")))
}Step 1: How hard was each kick?
Theory
The first step estimates the chance an average kicker makes each kick. Because the outcome is made or missed, a logistic regression is used: it models the log-odds of a make as a sum of effects, then converts that into a probability between 0 and 1. The model is fit on every kick by every kicker, so it describes the league as a whole rather than any one kicker. Only the conditions of a kick enter the model, never its outcome.
\[\hat{p}_i = \frac{1}{1 + e^{-\left(\beta_0 + f(D_i) + G_i + E_i\right)}}\]
The \(\beta\) coefficients are estimated from the data, so the model decides how much each factor matters rather than us guessing.
Distance: a natural spline
Distance is by far the most important factor, but its effect isn’t a straight line: kicks are almost automatic up to about 30 yards, then get harder faster and faster. Rather than forcing a fixed shape like a quadratic, distance enters as a natural cubic spline \(f(D_i)\). A spline joins smooth cubic curves at points called knots, here at 25, 33, 40, 50 and 55 yards, so the curve can bend where difficulty changes. “Natural” means it becomes a straight line beyond the shortest and longest kicks, which keeps it from swinging wildly where data is thin.
Kick difficulty and weather
\[G_i = \beta_1\varphi_i + \beta_2 S_i + \beta_3 X_i\]
- \(\varphi_i\): kicking direction, the angle of the hold spot off center, signed toward the kicking-foot side.
- \(S_i\): 1 for natural grass, 0 for turf.
- \(X_i\): 1 for an extra point, 0 for a field goal.
\[E_i = \beta_4 W_i + \beta_5 T_i + \beta_6 C_i + \beta_7 R_i + \beta_8 N_i + \beta_9 A_i\]
- \(W_i\): crosswind and \(T_i\): tailwind, both in mph (0 indoors).
- \(C_i\): degrees below 50°F.
- \(R_i\), \(N_i\): 1 for rain or snow.
- \(A_i\): thin-air measure from elevation and temperature, which matters in Denver and Mexico City.
The kicking-window angle (how wide the uprights look from the hold spot) is left out because distance almost fully determines it, and including both would make the model unstable.
Fitting the model
m1 <- glm(make_formula("made", KNOTS), data = d, family = binomial)
d[, p_hat := predict(m1, type = "response")]
co1 <- as.data.table(summary(m1)$coefficients, keep.rownames = "term")
setnames(co1, c("term", "estimate", "std_error", "z_value", "p_value"))
co1[, odds_ratio := exp(estimate)]
co1[, term := sub("^ns\\(D.*\\)\\)([0-9]+)$", "spline basis \\1", term)]
kable(co1[, .(term, estimate, std_error, odds_ratio, p_value)], digits = 4,
col.names = c("Term", "Estimate", "Std. error", "Odds ratio", "p-value"))| Term | Estimate | Std. error | Odds ratio | p-value |
|---|---|---|---|---|
| (Intercept) | 7.2332 | 1.9198 | 1384.6997 | 0.0002 |
| spline basis 1 | -3.8275 | 1.8193 | 0.0218 | 0.0354 |
| spline basis 2 | -5.8223 | 2.0080 | 0.0030 | 0.0037 |
| spline basis 3 | -6.0906 | 1.8991 | 0.0023 | 0.0013 |
| spline basis 4 | -5.5806 | 1.1830 | 0.0038 | 0.0000 |
| spline basis 5 | -11.1013 | 3.9652 | 0.0000 | 0.0051 |
| spline basis 6 | -6.6024 | 1.2481 | 0.0014 | 0.0000 |
| phi | -0.0045 | 0.0132 | 0.9955 | 0.7303 |
| grass | -0.1543 | 0.0969 | 0.8570 | 0.1113 |
| is_xp | -0.2406 | 0.1917 | 0.7862 | 0.2095 |
| crosswind | -0.0174 | 0.0131 | 0.9828 | 0.1838 |
| tailwind | 0.0066 | 0.0110 | 1.0066 | 0.5487 |
| cold | 0.0024 | 0.0100 | 1.0024 | 0.8065 |
| rain | 0.0065 | 0.2021 | 1.0065 | 0.9744 |
| snow | 12.0955 | 229.1683 | 179061.4706 | 0.9579 |
| thin_air | 1.4885 | 1.3190 | 4.4304 | 0.2591 |
The six spline terms have no meaning on their own; together they form the distance curve. An odds ratio above 1 means a factor makes a kick more likely to be made, and below 1 means less likely.
Testing distance and weather as groups
The p-values in the table test each column on its own, but distance enters as six spline columns together and weather as six separate factors. The proper test for a group is a likelihood ratio test: refit the model without the group and check whether the fit gets worse by more than chance would explain.
f_noD <- as.formula(paste("made ~", paste(c(G_VARS, E_VARS), collapse = " + ")))
f_noE <- update(formula(m1), paste(". ~ . -", paste(E_VARS, collapse = " - ")))
lr <- function(small, label) {
a <- anova(glm(small, data = d, family = binomial), m1, test = "LRT")
data.table(Group = label, df = a$Df[2], `Chi-square` = a$Deviance[2], `p-value` = a$`Pr(>Chi)`[2])
}
lr_tab <- rbind(lr(f_noD, "Distance spline f(D)"), lr(f_noE, "Weather factors (all six)"))
lr_p <- lr_tab[["p-value"]]
fmt_p <- function(p) ifelse(p < 0.001, "< 0.001", sprintf("= %.3f", p))
lr_tab[, `p-value` := sub("= ", "", fmt_p(`p-value`))]
kable(lr_tab, digits = 1)| Group | df | Chi-square | p-value |
|---|---|---|---|
| Distance spline f(D) | 6 | 362.4 | < 0.001 |
| Weather factors (all six) | 6 | 6.0 | 0.424 |
Interpreting the model
Distance dominates. Removing the spline worsens the fit by a chi-square of 362.4 on 6 degrees of freedom (p < 0.001). Nothing else in the model comes close, which matches intuition: how far a kick is matters far more than anything else about it. The individual spline coefficients look extreme, but they have no meaning one at a time; only the curve they form together does, shown below.
The other factors are small. Of the nine factors besides distance, 0 are statistically significant at the 5% level, and the six weather factors together do not significantly improve the fit (p = 0.424). Their directions still mostly make physical sense. Each mph of crosswind lowers the odds of a make by about 1.7%, each mph of tailwind raises them by about 0.7%, and natural grass lowers them by about 14.3% compared with turf.
Why the weather effects are weak. Two things work against them. First, distance explains most misses, and with only 585 misses there is little left over to detect small effects. Second, the weather comes from per-quarter forecasts rather than conditions at the moment of the kick, and that measurement error pulls estimated effects toward zero. The factors are kept anyway: the goal of Step 1 is a fair baseline, and each one has a clear physical reason to affect a kick.
Snow can’t be estimated. Every snow kick in the data was made, so the model has no misses to learn from and pushes the snow coefficient toward infinity, which is why its odds ratio and standard error are enormous. This is known as separation; the coefficient should be ignored rather than interpreted.
Visualization: the distance spline
The fitted spline is drawn for typical field goal conditions, against the actual make rate in 5-yard bins (larger dots have more kicks). Dotted lines mark the knots.
typical <- d[is_xp == 0, lapply(.SD, median), .SDcols = c(G_VARS, E_VARS)]
typical[, is_xp := 0]
grid <- data.table(D = seq(BOUNDARY[1], BOUNDARY[2], by = 0.5))[, (names(typical)) := typical]
grid[, fit := predict(m1, newdata = grid, type = "response")]
binned <- d[is_xp == 0, .(rate = mean(made), n = .N),
by = .(D5 = 5 * floor(D / 5) + 2.5)][order(D5)]
plot(grid$D, grid$fit, type = "l", lwd = 2, ylim = c(0, 1),
xlab = "Kick distance (yards)", ylab = "Chance of making it",
main = "Field goals: fitted make chance vs. actual rate")
points(binned$D5, binned$rate, pch = 19, cex = pmax(0.6, sqrt(binned$n) / 8), col = "grey30")
abline(v = KNOTS, lty = 3, col = "grey60")
legend("bottomleft", c("Spline model (typical conditions)", "Actual rate, 5-yard bins", "Knots"),
lty = c(1, NA, 3), pch = c(NA, 19, NA), lwd = c(2, NA, 1),
col = c("black", "grey30", "grey60"), bty = "n")
Interpreting the spline
For a field goal in typical conditions, an average kicker makes about 98% of kicks from 25 yards, 85% from 40, 70% from 50 and 43% from 60. The curve is nearly flat through the short range, then bends downward from the high 30s and falls steeply past 50 yards. That changing shape is exactly why a spline is used: a straight line would badly overstate short-kick difficulty, and a quadratic would force the same bend everywhere.
The dots follow the curve closely from 20 to 50 yards, where most kicks are. At the longest distances the actual rates tend to sit above the curve. That’s expected: coaches only attempt very long kicks with strong kickers and in favorable conditions, so those kicks are easier than “a 60-yarder in typical conditions.” The model accounts for each kick’s actual conditions, which is why the calibration check below still holds.
Visualization: the other factors
Each factor’s odds ratio is shown with its 95% confidence interval. Factors whose interval crosses 1 (the dashed line) can’t be separated from no effect. Continuous factors are shown per unit (one degree of angle, one mph, one degree of cold); thin air is shown per 0.01 change in air density.
ci <- confint.default(m1)
or <- data.table(term = rownames(ci), est = coef(m1), lo = ci[, 1], hi = ci[, 2])
or <- or[term %in% c(G_VARS, E_VARS)]
or[term == "thin_air", `:=`(est = est / 100, lo = lo / 100, hi = hi / 100)]
or[, label := c(phi = "Kicking angle (per degree)", grass = "Natural grass",
is_xp = "Extra point", crosswind = "Crosswind (per mph)",
tailwind = "Tailwind (per mph)", cold = "Cold (per degree)",
rain = "Rain", snow = "Snow", thin_air = "Thin air (per 0.01)")[term]]
or <- or[is.finite(hi) & hi - lo < 20] # drops terms with no usable estimate
or <- or[order(est)]
par(mar = c(4.5, 13, 2, 1))
plot(exp(or$est), seq_len(nrow(or)), log = "x", pch = 19, yaxt = "n", ylab = "",
xlim = range(exp(c(or$lo, or$hi))), xlab = "Odds ratio (log scale)",
main = "Effect of each factor on the chance of a make")
segments(exp(or$lo), seq_len(nrow(or)), exp(or$hi), seq_len(nrow(or)), lwd = 2)
abline(v = 1, lty = 2, col = "grey50")
axis(2, at = seq_len(nrow(or)), labels = or$label, las = 1)
Any factor missing from this chart (such as snow) had no misses to learn from, so its effect can’t be estimated.
Interpreting the factor chart
A dot to the right of the dashed line means the factor makes kicks easier, and to the left means harder. The bar around each dot is the range of plausible values; when it crosses the dashed line, the data can’t rule out no effect at all. Narrow bars that sit near 1 (such as wind per mph) mean the effect is measured precisely and is small, while wide bars (such as rain) mean there are too few kicks in those conditions to say much either way.
Step 2: How accurate was each kick?
Theory
Making a kick is a yes-or-no outcome, but accuracy shows more detail: a kick that splits the uprights shows more skill than one that sneaks inside the post. Aim error is measured as an angle rather than in yards, so a 2-yard miss on a short kick counts as a bigger mistake than the same miss from 55 yards:
\[\alpha_i = \arctan\left(\frac{|e_i|}{D_i}\right)\]
A linear regression on the same variables as Step 1 gives the expected angle \(\hat{\alpha}_i\) for each kick. A kicker’s Aim Error Over Expected is his average gap between actual and expected:
\[\text{AEOE}_k = \frac{1}{n_k}\sum_{i \in k}\left(\alpha_i - \hat{\alpha}_i\right)\]
Lower AEOE means a more accurate kicker. The standard deviation of those gaps measures consistency: how much a kicker’s accuracy varies from kick to kick.
m2 <- lm(make_formula("alpha", KNOTS), data = d[!is.na(alpha)])
d[, alpha_hat := predict(m2, newdata = d)]
d[, aim_resid := alpha - alpha_hat] # + = less accurate than expectedInterpreting the accuracy model
The conditions explain only 1.3% of the variation in aim error (R² = 0.013). That might look like a weak model, but for this purpose it’s good news: it means how close to center a kick flies depends very little on distance, wind or surface, and mostly on the kicker and on luck. So AEOE mainly measures the kicker himself.
Accuracy also captures something make-or-miss can’t. Two kickers can make the same kicks while one consistently splits the uprights and the other keeps sneaking inside the post. Over time, the second kicker is more likely to start missing, so AEOE works as a second, independent view of skill.
Step 3: Points over expected
Theory
This step turns difficulty into a score. Each kick earns its point value times the gap between what happened and what an average kicker would be expected to do:
\[\text{POE}_k = \sum_{i \in k} v_i\left(y_i - \hat{p}_i\right)\]
- \(y_i\): 1 if the kick was made, 0 if missed.
- \(v_i\): the kick’s point value, 3 for a field goal and 1 for an extra point.
- The sum runs over all of kicker \(k\)’s kicks.
How the credit works
Take a 50-yard field goal that an average kicker makes 70% of the time. Making it earns \(3 \times (1 - 0.70) = +0.9\) points; missing it costs \(3 \times (0 - 0.70) = -2.1\) points. An extra point made 94% of the time earns only \(+0.06\) when made, because almost everyone makes it, but costs \(-0.94\) when missed. Hard makes earn a lot and easy misses cost a lot, which is exactly what make percentage ignores. An average kicker finishes near zero by design.
Interpreting points over expected
Added up over every kick in the league, POE comes to +0.0 points, close to zero. That is a useful check: the model’s “average kicker” really is the league average, so a kicker’s POE measures him against his peers. Among kickers with 30+ attempts, totals range from +34.4 (Justin Tucker) to -10.7 (Adam Vinatieri). Because POE is measured in points, it translates directly into what a kicker is worth on the scoreboard.
Step 4: Standardize and rank
Theory
Total POE grows with the number of kicks, so a busy kicker would rank above an equally good kicker with fewer chances. POE is therefore put on a per-kick basis, with a small adjustment for sample size:
\[\text{POE rate}_k = \frac{\text{POE}_k}{n_k + 25}\]
The +25 pulls kickers with few attempts slightly toward zero, so a kicker who made a handful of lucky kicks can’t top the list. Its effect fades as attempts grow.
The rates are then standardized against the kickers with 30+ attempts, the only ones ranked:
\[z_k = \frac{\text{POE rate}_k - \bar{x}}{s}\]
where \(\bar{x}\) and \(s\) are the mean and standard deviation of their POE rates. A \(z\) of +2 means far better than a typical kicker; 0 means exactly average. AEOE and consistency are standardized the same way, signed so that positive always means better.
For readability, each POE rate is also converted into points added per season: the rate times the number of kicks an average team attempts in a season.
k <- d[, .(
attempts = .N,
fg_att = sum(is_xp == 0),
fg_made = sum(made[is_xp == 0]),
xp_att = sum(is_xp == 1),
xp_made = sum(made[is_xp == 1]),
POE = sum(poe),
AEOE = mean(aim_resid, na.rm = TRUE),
consistency_sd = sd(aim_resid, na.rm = TRUE)
), by = .(kickerId, kicker)]
k[, POE_rate := POE / (attempts + SHRINK)]
k[, qualifies := as.integer(attempts >= MIN_ATTEMPTS)]
zscore <- function(x, q) (x - mean(x[q == 1])) / sd(x[q == 1])
k[, z_POE := zscore(POE_rate, qualifies)]
k[, z_accuracy := -zscore(AEOE, qualifies)]
k[, z_consistency := -zscore(consistency_sd, qualifies)]
k[qualifies == 0, c("z_POE", "z_accuracy", "z_consistency") := NA_real_]
setorder(k, -qualifies, -z_POE)
k[qualifies == 1, rank := seq_len(.N)]
# Average kicks per team-season, using the kicking team from plays.csv
DATA_DIR <- "nfl-big-data-bowl-2022"
plays_csv <- file.path(DATA_DIR, "plays.csv")
has_plays <- file.exists(plays_csv)
if (has_plays) {
plays <- fread(plays_csv, select = c("gameId", "playId", "possessionTeam"))
kicks <- plays[d[, .(gameId, playId, kickerId, season)], on = .(gameId, playId)]
SEASON_KICKS <- nrow(kicks) / uniqueN(kicks[!is.na(possessionTeam), .(possessionTeam, season)])
} else {
SEASON_KICKS <- nrow(d) / (32 * uniqueN(d$season)) # 32 teams per season
}
k[, POE_per_season := POE_rate * SEASON_KICKS]
# Team colors: each kicker gets the color of the team he kicked for most
TEAM_COLORS <- c(
ARI = "#97233F", ATL = "#A71930", BAL = "#241773", BUF = "#00338D",
CAR = "#0085CA", CHI = "#C83803", CIN = "#FB4F14", CLE = "#311D00",
DAL = "#003594", DEN = "#FB4F14", DET = "#0076B6", GB = "#203731",
HOU = "#03202F", IND = "#002C5F", JAX = "#006778", KC = "#E31837",
LA = "#003594", LAC = "#0080C6", LV = "#000000", OAK = "#000000",
MIA = "#008E97", MIN = "#4F2683", NE = "#002244", NO = "#D3BC8D",
NYG = "#0B2265", NYJ = "#125740", PHI = "#004C54", PIT = "#FFB612",
SEA = "#69BE28", SF = "#AA0000", TB = "#D50A0A", TEN = "#4B92DB",
WAS = "#5A1414")
q <- k[qualifies == 1]
if (has_plays) {
teams <- kicks[!is.na(possessionTeam), .N, by = .(kickerId, possessionTeam)][order(kickerId, -N)]
q <- teams[, .(main_team = possessionTeam[1],
all_teams = paste0(possessionTeam, " (", N, ")", collapse = ", ")),
by = kickerId][q, on = "kickerId"]
} else {
q[, `:=`(main_team = NA_character_, all_teams = "not available")]
}
q[, color := fifelse(main_team %in% names(TEAM_COLORS), TEAM_COLORS[main_team], "#888888")]
q[, hover := sprintf(paste0(
"<b>%s</b> (#%d)<br>Team: %s<br>Skill: %+.2f SD from average<br>",
"Points added per kick: %+.3f<br>Points added per %.0f-kick season: %+.1f<br>",
"Total points over expected: %+.1f<br>Accuracy: %+.2f SD<br>",
"FG %d/%d | XP %d/%d | %d attempts"),
kicker, rank, all_teams, z_POE, POE_rate, SEASON_KICKS, POE_per_season,
POE, z_accuracy, fg_made, fg_att, xp_made, xp_att, attempts)]
best <- q[order(rank)][1]An average team attempts 63.1 kicks per season, which is the season length used below.
Results
The top 15 kickers
q[order(rank)][1:min(15, .N),
.(Rank = rank, Kicker = kicker, Attempts = attempts,
FG = paste0(fg_made, "/", fg_att), XP = paste0(xp_made, "/", xp_att),
`Skill (z)` = sprintf("%+.2f", z_POE),
`Points per kick` = sprintf("%+.3f", POE_rate),
`Points per season` = sprintf("%+.1f", POE_per_season),
`Total POE` = sprintf("%+.1f", POE),
`Accuracy (z)` = sprintf("%+.2f", z_accuracy),
`Consistency (z)` = sprintf("%+.2f", z_consistency))] |>
kable(align = "rlrrrrrrrrr")| Rank | Kicker | Attempts | FG | XP | Skill (z) | Points per kick | Points per season | Total POE | Accuracy (z) | Consistency (z) |
|---|---|---|---|---|---|---|---|---|---|---|
| 1 | Justin Tucker | 228 | 87/92 | 133/136 | +2.60 | +0.136 | +8.6 | +34.4 | +2.60 | +1.88 |
| 2 | Josh Lambo | 105 | 54/57 | 45/48 | +2.42 | +0.127 | +8.0 | +16.4 | +1.68 | +1.08 |
| 3 | Graham Gano | 96 | 42/45 | 46/51 | +1.86 | +0.097 | +6.1 | +11.7 | +0.47 | +0.25 |
| 4 | Brandon McManus | 164 | 72/82 | 78/82 | +1.36 | +0.070 | +4.4 | +13.2 | +0.37 | +0.40 |
| 5 | Wil Lutz | 231 | 77/85 | 144/146 | +1.26 | +0.065 | +4.1 | +16.5 | +1.14 | +1.78 |
| 6 | Nick Folk | 82 | 38/41 | 38/41 | +1.25 | +0.064 | +4.1 | +6.9 | +0.96 | +1.11 |
| 7 | Jason Myers | 197 | 73/80 | 106/117 | +1.22 | +0.063 | +3.9 | +13.9 | -0.06 | -0.17 |
| 8 | Mason Crosby | 201 | 63/71 | 124/130 | +0.99 | +0.050 | +3.2 | +11.4 | +1.20 | +0.74 |
| 9 | Jason Sanders | 173 | 71/83 | 88/90 | +0.88 | +0.044 | +2.8 | +8.7 | +0.84 | +0.84 |
| 10 | Dustin Hopkins | 152 | 73/85 | 63/67 | +0.78 | +0.039 | +2.5 | +6.9 | +0.47 | -0.02 |
| 11 | Harrison Butker | 238 | 77/84 | 143/154 | +0.70 | +0.035 | +2.2 | +9.2 | +0.79 | +0.78 |
| 12 | Younghoe Koo | 101 | 52/57 | 40/44 | +0.51 | +0.025 | +1.5 | +3.1 | +0.73 | +0.71 |
| 13 | Randy Bullock | 160 | 66/77 | 79/83 | +0.46 | +0.022 | +1.4 | +4.1 | +0.78 | +0.18 |
| 14 | Kai Forbath | 37 | 16/18 | 17/19 | +0.31 | +0.014 | +0.9 | +0.9 | -2.07 | -2.10 |
| 15 | Greg Zuerlein | 198 | 76/92 | 104/106 | +0.29 | +0.013 | +0.8 | +2.9 | +0.17 | +0.28 |
Interpreting the rankings
Justin Tucker ranks first, 2.60 standard deviations above the average qualifying kicker, ahead of Josh Lambo at 2.42. Only 7 of 45 kickers sit more than one standard deviation above average, while 33 fall within one standard deviation either side. Once difficulty is accounted for, most NFL kickers are close to interchangeable, and a small group stands out.
Two things are worth weighing when reading the table. First, kickers with fewer attempts have less certain positions: a couple of kicks either way can move them several places, even with the +25 adjustment. Second, the accuracy and consistency columns are independent of the ranking. A kicker who ranks high on both points and accuracy has stronger evidence behind his position than one who ranks high on points alone.
Visualization: the kicker skill curve
This is the main visualization. Each dot is a kicker with 30+ attempts, placed on a bell curve by how many standard deviations his points over expected per kick sits from the average kicker. Dots are colored by the team he kicked for most, and the top three are labeled. Hover over any dot to see the kicker’s name, team and stats.
q[, y := dnorm(z_POE)]
xmax <- max(3.2, ceiling(max(abs(q$z_POE)) * 10) / 10 + 0.3)
curve <- data.table(x = seq(-xmax, xmax, length.out = 400))[, y := dnorm(x)]
top <- q[order(rank)][1:min(3, .N)]
plot_ly() |>
add_lines(data = curve, x = ~x, y = ~y, hoverinfo = "none",
line = list(color = "rgba(90,90,90,0.8)", width = 2),
fill = "tozeroy", fillcolor = "rgba(150,150,150,0.12)") |>
add_markers(data = q, x = ~z_POE, y = ~y, text = ~hover, hoverinfo = "text",
marker = list(color = ~color, size = 15, opacity = 0.95,
line = list(color = "white", width = 1.5))) |>
layout(
title = list(text = "<b>NFL kicker skill, 2018–2020</b>"),
xaxis = list(title = "Standard deviations from the average kicker (POE z-score)",
range = c(-xmax, xmax), dtick = 1),
yaxis = list(title = "", showticklabels = FALSE, showgrid = FALSE, range = c(0, 0.46)),
shapes = lapply(c(-2, -1, 1, 2), function(s) list(
type = "line", x0 = s, x1 = s, y0 = 0, y1 = dnorm(s),
line = list(color = "rgba(120,120,120,0.5)", dash = "dot", width = 1))),
annotations = lapply(seq_len(nrow(top)), function(i) list(
x = top$z_POE[i], y = top$y[i], text = top$kicker[i], showarrow = TRUE,
arrowhead = 0, ax = 0, ay = -28 - 14 * (i - 1))),
showlegend = FALSE, hoverlabel = list(align = "left"), margin = list(t = 60))Visualization: points added per season
The same ranking in points: how many more points each of the top 15 kickers adds in an average season than an average kicker taking the same kicks.
t15 <- q[order(rank)][1:min(15, .N)]
plot_ly(t15, x = ~POE_per_season, y = ~factor(kicker, levels = rev(kicker)),
type = "bar", orientation = "h", text = ~hover, hoverinfo = "text",
textposition = "none", marker = list(color = ~color)) |>
layout(xaxis = list(title = "Points added per season over an average kicker"),
yaxis = list(title = ""), margin = list(l = 140))Visualization: points vs. accuracy
Points over expected and accuracy are measured independently: one from whether kicks went in, the other from how close to center they flew. Kickers in the upper right are strong on both, which makes their ranking more convincing than either measure alone.
plot_ly(q, x = ~z_accuracy, y = ~z_POE, type = "scatter", mode = "markers",
text = ~hover, hoverinfo = "text",
marker = list(color = ~color, size = ~pmax(8, sqrt(attempts) * 1.6), opacity = 0.9,
line = list(color = "white", width = 1))) |>
layout(xaxis = list(title = "Accuracy (z-score, + = more accurate)", zeroline = TRUE),
yaxis = list(title = "Points over expected (z-score)", zeroline = TRUE),
annotations = lapply(seq_len(nrow(top)), function(i) list(
x = top$z_accuracy[i], y = top$z_POE[i], text = top$kicker[i],
showarrow = FALSE, yshift = 16)))Dot size reflects the number of attempts, so larger dots have more reliable placements.
Model checks
Calibration
Every score is built on \(\hat{p}_i\), so it matters most that the predicted chances match what actually happened. Kicks are split into 10 equal groups by predicted chance, and each group’s average prediction is compared with its actual make rate. Points on the dashed line mean the model is predicting accurately.
d[, cal_bin := cut(p_hat, breaks = unique(quantile(p_hat, 0:10 / 10)),
include.lowest = TRUE, labels = FALSE)]
cal <- d[, .(Kicks = .N, Predicted = mean(p_hat), Actual = mean(made)), by = cal_bin][order(cal_bin)]
lims <- range(c(cal$Predicted, cal$Actual))
plot(cal$Predicted, cal$Actual, pch = 19, xlim = lims, ylim = lims,
xlab = "Predicted make chance", ylab = "Actual make rate",
main = "Calibration: predicted vs. actual (10 groups)")
abline(0, 1, lty = 2, col = "grey40")
| Group | Kicks | Predicted | Actual |
|---|---|---|---|
| 1 | 606 | 0.664 | 0.665 |
| 2 | 605 | 0.801 | 0.795 |
| 3 | 606 | 0.907 | 0.914 |
| 4 | 605 | 0.935 | 0.931 |
| 5 | 606 | 0.940 | 0.936 |
| 6 | 605 | 0.944 | 0.944 |
| 7 | 605 | 0.948 | 0.954 |
| 8 | 606 | 0.952 | 0.962 |
| 9 | 605 | 0.957 | 0.945 |
| 10 | 606 | 0.985 | 0.988 |
The Brier score (the average squared gap between prediction and result, lower is better) is 0.0779.
Interpreting calibration
Across the 10 groups, the largest gap between predicted and actual make rate is 1.1 percentage points, from the hardest kicks (about 66% predicted) to the easiest (about 99%). The model’s difficulty estimates match reality, which is the most important property for this method, because every kicker’s credit is built on them.
The Brier score is 0.0779, compared with 0.0873 for simply guessing the league make rate on every kick, an improvement of 11%. That sounds modest, but most kicks are near-certain makes, so there’s little room to improve on them; the gain comes from the hard kicks, which is where it matters.
Multicollinearity
When predictors overlap heavily, their coefficients become unstable. This is checked with the generalized VIF: for a term with \(m\) columns, \(\text{GVIF}^{1/(2m)}\) is compared with \(\sqrt{10} \approx 3.16\), the equivalent of the usual VIF > 10 rule.
gvif <- function(model) {
X <- model.matrix(model)[, -1, drop = FALSE]
assign <- attr(model.matrix(model), "assign")[-1]
labels <- attr(terms(model), "term.labels")
keep <- apply(X, 2, sd) > 0
X <- X[, keep, drop = FALSE]; assign <- assign[keep]
R <- cor(X); detR <- det(R)
rbindlist(lapply(unique(assign), function(j) {
a <- which(assign == j)
g <- det(R[a, a, drop = FALSE]) * det(R[-a, -a, drop = FALSE]) / detR
data.table(Term = labels[j], df = length(a), `Adjusted GVIF` = g^(1 / (2 * length(a))))
}))
}
gv <- gvif(m1)
gv[Term %like% "^ns\\(D", Term := "f(D) spline"]
kable(gv, digits = 2)| Term | df | Adjusted GVIF |
|---|---|---|
| f(D) spline | 6 | 1.10 |
| phi | 1 | 1.00 |
| grass | 1 | 1.04 |
| is_xp | 1 | 1.76 |
| crosswind | 1 | 1.06 |
| tailwind | 1 | 1.00 |
| cold | 1 | 1.15 |
| rain | 1 | 1.02 |
| snow | 1 | 1.01 |
| thin_air | 1 | 1.14 |
Interpreting multicollinearity
The largest adjusted GVIF is 1.76, well below the 3.16 cutoff. The predictors don’t overlap enough to destabilize each other, so each coefficient reflects its own factor. This is partly by design: the kicking-window angle was dropped because it nearly duplicated distance, and kicking direction was used in its place.
Knot placement
To check that the results don’t depend on exactly where the knots are, the model is refit with two alternative knot sets and compared by AIC, which rewards fit but penalizes extra knots. Differences under 2 are effectively a tie.
knot_sets <- c(list(KNOTS), ALT_KNOTS)
names(knot_sets)[1] <- paste0(paste(KNOTS, collapse = ", "), " (chosen)")
aics <- sapply(knot_sets, function(kn)
AIC(glm(make_formula("made", kn), data = d, family = binomial)))
data.table(Knots = names(aics), AIC = aics, `vs. chosen` = aics - aics[1]) |>
kable(digits = 1)| Knots | AIC | vs. chosen |
|---|---|---|
| 25, 33, 40, 50, 55 (chosen) | 3367.8 | 0.0 |
| 30, 40, 50 | 3366.0 | -1.8 |
| 30, 36, 42, 48, 55 | 3368.6 | 0.8 |
Interpreting the knot comparison
The three knot sets differ by at most 1.8 AIC points. That is a tie: all three fit equally well, so the results don’t depend on exactly where the knots are placed. The chosen knots are kept because they sit at meaningful distances, including 33 yards, the length of an extra point.
Sensitivity: extended aim errors
Some kicks had their goal-line crossing extended because tracking stopped early. Step 2 is rerun without them to see whether that changes the accuracy rankings.
d_s <- d[e_extrapolated == 0 & !is.na(alpha)]
m2s <- lm(make_formula("alpha", KNOTS), data = d_s)
d_s[, aim_resid := alpha - predict(m2s, newdata = d_s)]
cmp <- merge(k[qualifies == 1, .(kickerId, AEOE)],
d_s[, .(AEOE_s = mean(aim_resid)), by = kickerId], by = "kickerId")
rho <- cor(cmp$AEOE, cmp$AEOE_s, method = "spearman")Removing the 93 extended kicks gives a rank correlation of 0.971 with the full accuracy ranking, where 1 would mean no change at all.
Interpreting the sensitivity check
A rank correlation this close to 1 means the extended kicks barely move anyone’s accuracy ranking. Extending the ball’s path for the few kicks where tracking stopped early doesn’t drive the results, so keeping them is safe.
Conclusion
By this method, the best NFL kicker from 2018 to 2020 was Justin Tucker. Across 228 attempts, he scored +34.4 points over what an average kicker would have from the same kicks, about 8.6 points per season, placing him 2.60 standard deviations above the average qualifying kicker. His accuracy z-score of +2.60 shows whether his kicks also flew closer to center than expected, a second, independent check on the result.
The strength of this approach is that it compares every kicker on the same footing: each is judged only against what an average kicker would have done with his exact kicks, so a heavy workload of long or windy kicks no longer counts against him.
Assumptions and limitations
- Weather is forecast, not measured. It comes from per-quarter forecasts, with wind direction to 8 compass points, which likely weakens the weather effects.
- Retractable roofs are assumed closed, since roof status for each game isn’t available.
- Kicking foot sets the sign of \(\varphi_i\), so it must be correct for every kicker.
- Rare conditions can’t be estimated: if a condition (such as snow) has no misses, the model can’t measure its effect.
- Small samples: kickers with fewer attempts have less certain rankings, even with the +25 adjustment.
Data: NFL Big Data Bowl 2022. Forecasts provided by NFLWeather.com.