library(nflreadr)
library(dplyr)
library(ggplot2)
library(pROC)
library(tidyverse)
library(patchwork)
#Player directory: draft info + birth date, one row per player
player_info <- load_players() |>
#just pulling biographical info about a player
#id, birth date, draft info, and year of their rookie season (undrafted players don't have a draft year)
select(gsis_id, birth_date, draft_year, draft_round, draft_pick, rookie_season) |>
#keeping only one instance of each player
distinct(gsis_id, .keep_all = TRUE) |>
#change the birth_date variable to a datetime variable
mutate(birth_date = as.Date(birth_date))
#want to get juset wr stats over the last 20 years, and just pulling offense players for the sake of memory and data processing time
wr_stats <- load_player_stats(seasons = 2006:2025) |>
#just want regular season stats from WR
filter(season_type == "REG", position == "WR") |>
#grouping by season, id, name, display name, in order to add other stats to individual player season totals (probably only need to group by season and player_id but this is fine too)
group_by(season, player_id, player_name, player_display_name) |>
#for each unqique player and unique season, find number games, targets, reception, receiving yards, receiving tds, air yards, rec_epa
summarise(
games = n_distinct(week),
targets = sum(targets, na.rm = TRUE),
receptions = sum(receptions, na.rm = TRUE),
rec_yards = sum(receiving_yards, na.rm = TRUE),
rec_tds = sum(receiving_tds, na.rm = TRUE),
air_yards = sum(receiving_air_yards, na.rm = TRUE), #reveiving_air_yards is a weird stat, its air yards of the ball on each TARGET including incompletions, this number is frequently bigger that rec_yrds because of that
rec_epa = sum(receiving_epa, na.rm = TRUE),
.groups = "drop" #use .groups = "drop" to ungroup after getting the necessary stats
) |>
#we want receivers that see some field action, 20 total targets in any given season feels like a good amount
filter(targets >= 20) |>
#adding yards, targets, yards per target, for each season
mutate(
yards_per_game = rec_yards / games,
targets_per_game = targets / games,
yards_per_target = rec_yards / targets
) |>
#joining the player_info df (contains draft info, birth day) onto the main stats df
#player_id from main stats df equals gsis_id from player_info df
left_join(player_info, by = c("player_id" = "gsis_id")) |>
#grouping by player_id -- add info for each player in each season
group_by(player_id) |>
#add the first season for a player with min(season) (the df is grouped by player id)
mutate(first_season_in_data = min(season)) |>
#ungroup after pulling a player's first season
ungroup() |>
#addding draft info abount each player
mutate(
#drafted binary logical variable needs to be made BEFORE adding dummy variables to undrafted players
#if a player has a draft_pick, drafted is TRUE, else FALSE
drafted = !is.na(draft_pick),
#coalese() returns the FIRST non-NA value it finds, so draft year will essentially just end up being the first year a player has any recorded stats
draft_year = coalesce(draft_year, rookie_season, first_season_in_data),
draft_pick = if_else(drafted, as.numeric(draft_pick), 263), #assign number for undrafted players
draft_round = if_else(drafted, as.numeric(draft_round), 8), #assign round for undrafted players
age = as.numeric(difftime(as.Date(paste0(season, "-09-01")),
birth_date, units = "days")) / 365.25,
#age floor rounds age down, so that you can group players by age number later on (age is in decimal for some reason)
age_floor = floor(age),
#find number of years in the league for each player each season
yrs_in_league = season - draft_year
) |>
#the rookie_season and first_season_in_data columns were just use to create the draft_year column, they're not needed anymore
select(-rookie_season, -first_season_in_data) |>
#put all these columns after player_id, makes it cleanere
relocate(player_display_name, player_name, birth_date, age, age_floor,
draft_year, draft_round, draft_pick, drafted,
.after = player_id)
#general production stats grouped by age
age_curve <- wr_stats |>
filter(!is.na(age_floor)) |>
group_by(age_floor) |>
summarise(
player_seasons = n(),
avg_targets_pg = mean(targets_per_game),
avg_yards_pg = mean(yards_per_game),
avg_ypt = mean(yards_per_target),
avg_tds = mean(rec_tds),
avg_epa = mean(rec_epa, na.rm = TRUE),
.groups = "drop"
)
age_curve
wr_theme <- theme_minimal(base_family = 'mono', base_size = 14) +
theme(
plot.margin = margin(20, 20, 20, 20),
text = element_text(color = 'black', face = 'bold'),
plot.title = element_text(color = 'black', face = 'bold', size = 16),
axis.title = element_text(color = 'black', face = 'bold'),
axis.text = element_text(color = 'black', face = 'bold')
)
#the goal is to see at what point in a WR's career are they at their best, making extremely quick plots about age and number of years in the league
#how many instances of each age are there in the dataset
age <- wr_stats %>% ggplot(aes(x = age)) +
geom_histogram(binwidth = 1, fill = 'royalblue3', color = 'white') +
wr_theme +
labs(title = 'Counts of WR by Age', x = 'Age', y = 'Counts')
age
#how many instances of career years in the league
years <- wr_stats %>% ggplot(aes(x = yrs_in_league)) +
geom_histogram(binwidth = 1, fill = 'royalblue3', color = 'white') +
wr_theme +
labs(title = 'Counts of WR by Career Years', x = 'Years in League', y = 'Counts')
years
#age for a season vs yards per game
age_v_yards <- wr_stats %>% ggplot(aes(x = age, y = yards_per_game)) +
geom_point(color = 'royalblue3', alpha = .5) +
wr_theme +
labs(title = 'Rec Yards Per Game vs Age', x = 'Age', y = 'Receiving Yards Per Game')
age_v_yards
#number of years spent in the league for a given seaseon vs yards per game
years_v_yards <- wr_stats %>% ggplot(aes(x = yrs_in_league, y = yards_per_game)) +
geom_point(color = 'royalblue3', alpha = .5) +
wr_theme +
labs(title = 'Rec Yards Per Game vs Career Years', x = 'Years in the League', y = 'Receiving Yards per Game')
years_v_yards
combined <- (age | years) / (age_v_yards | years_v_yards)
ggsave(filename = "~/Documents/elite_wr/wr_plot.png", combined,
width = 12, height = 9, dpi = 300, bg = "white")
combined
#elite wrs
#defining an "elite" seaon as having either 1000 receiving yards or 100 catches (rarely is it one or the other however)
elite_wr_stats <- wr_stats %>% filter((rec_yards >= 1000) | (receptions >= 100), drafted == TRUE)
#all recievers who never had either 1000 yards or 100 catches are non-elite
non_elite_wr_stats <- wr_stats %>% filter((rec_yards < 1000) & (receptions < 100), drafted == TRUE)
#making the same quick plots but for elite receivers only
age <- elite_wr_stats %>% ggplot(aes(x=age)) +
geom_histogram(binwidth=1,fill = 'royalblue3', color = 'white') +
wr_theme +
labs(title='Counts of Elite WR by Age', x='Age', y='Counts')
age
years <- elite_wr_stats %>% ggplot(aes(x=yrs_in_league)) +
geom_histogram(binwidth=1, fill = 'royalblue3', color = 'white') +
wr_theme +
labs(title = 'Counts of Elite WR by Career Years', x ='Years in League', y='Counts')
years
age_v_yards <- elite_wr_stats %>% ggplot(aes(x=age, y=yards_per_game)) +
geom_point(color = 'royalblue3', alpha = .5) +
wr_theme +
labs(title = 'Elite WR Rec Yards Per Game vs Age', x ='Age', y='Receiving Yards Per Game')
age_v_yards
years_v_yards <- elite_wr_stats %>% ggplot(aes(x=yrs_in_league, y=yards_per_game)) +
geom_point(color = 'royalblue3', alpha = .5) +
wr_theme +
labs(title = str_wrap('Elite WR Rec Yards Per Game vs Career Years', width=40), x ='Years in the League', y='Receiving Yards per Game')
years_v_yards
combined <- (age | years) / (age_v_yards | years_v_yards)
ggsave(filename = "~/Documents/elite_wr/elite_wr_plot.png", combined, width = 12, height = 9, dpi = 300, bg = "white")
combined
#going to make a histogram that combines elite wrs and non-elite wrs
overlay_df <- bind_rows(
non_elite_wr_stats |> transmute(age = age, group = "Non-Elite"),
elite_wr_stats |> transmute(age = age, group = "Elite")
)
#getting the means ages of elite and non-elite for visualization
#making a df that will contain the mean, the group, and a label that will go on the final graph
means <- overlay_df |>
#the group column is "Elite" or "Non-Elite"
group_by(group) |>
summarise(mean_val = mean(age, na.rm = TRUE), .groups = "drop") |>
#adding a label column that will read "[group name] mean: [value]"
mutate(label = paste0(group, " mean age: ", round(mean_val, 1)))
#(new) pick a color for each group here, any R color name or hex code works
group_colors <- c("Elite" = "firebrick3", "Non-Elite" = "royalblue3")
#using the overlay_df data for making the overlayed histogram
#in aes() x is age, and fill color will be by group
p <- ggplot(overlay_df, aes(x = age, fill = group)) +
#after_stat(density) tells ggplot to maek the y-axis density (within group) after it's done the binning
#position = "identity" tells ggplot to keep bars overlapped, "stack" would make them go on top of each other
geom_histogram(aes(y = after_stat(density)), # density so unequal group sizes are comparable
#color "white" makes the outlines white so that bars ear easy to see
position = "identity", alpha = 0.6, bins = 30, color = "white") +
#add the vertical line with the means of each group, use the means df from before
geom_vline(data = means, aes(xintercept = mean_val, color = group),
#show.legend = FALSE because we're labeling the lines directly
linewidth = 1, linetype = "dashed", show.legend = FALSE) +
#add the labels to the mean lines
geom_label(data = means,
#y = Inf so that the line goes up to the top of the graph
#use data = means and pull the label from the label column
aes(x = mean_val, y = Inf, label = label, color = group),
vjust = c(1.5, 3.2), hjust = -0.05, # staggered so labels don't collide
fill = "white", show.legend = FALSE) +
#(new) fill colors the bars, color colors the mean lines and label text
scale_fill_manual(values = group_colors) +
scale_color_manual(values = group_colors) +
labs(title = 'Proportion of Elite and Non-Elite WR Seasons by Age', x = "Age", y = "Density", fill = NULL) +
#using the mono font, looks better
wr_theme
p
ggsave(filename = "~/Documents/elite_wr/elite_v_non-e_age.png", p, width = 12, height = 9, dpi = 300, bg = "white")
#welch's T-test compare 2 means from different populations without assuming they have equal variances
t.test(age ~ group, data = overlay_df)
##
## Welch Two Sample t-test
##
## data: age by group
## t = 1.7881, df = 497.11, p-value = 0.07437
## alternative hypothesis: true difference in means between group Elite and group Non-Elite is not equal to 0
## 95 percent confidence interval:
## -0.03274627 0.69556750
## sample estimates:
## mean in group Elite mean in group Non-Elite
## 26.41623 26.08482
#primary output: t = 1.6944, df = 489.97, p-value = 0.09082
#the p-value is low but not low enough to be significant, fails to reject the null hypothesis that the mean ages between elite and non-elite wr seasons are equal
#doing everything again but instead of age, using yrs_in_league
overlay_df <- bind_rows(
non_elite_wr_stats |> transmute(yrs_in_league = yrs_in_league, group = "Non-Elite"),
elite_wr_stats |> transmute(yrs_in_league = yrs_in_league, group = "Elite")
)
means <- overlay_df |>
group_by(group) |>
summarise(mean_val = mean(yrs_in_league, na.rm = TRUE), .groups = "drop") |>
mutate(label = paste0(group, " mean years in league: ", round(mean_val, 1)))
group_colors <- c("Elite" = "firebrick3", "Non-Elite" = "royalblue3")
p <- ggplot(overlay_df, aes(x = yrs_in_league, fill = group)) +
geom_histogram(aes(y = after_stat(density)), # density so unequal group sizes are comparable
position = "identity", alpha = 0.5, binwidth = 1, color = "white") +
geom_vline(data = means, aes(xintercept = mean_val, color = group),
linewidth = 1, linetype = "dashed", show.legend = FALSE) +
geom_label(data = means,
aes(x = mean_val, y = Inf, label = label, color = group),
vjust = c(1.5, 3.2), hjust = -0.05, # staggered so labels don't collide
fill = "white", show.legend = FALSE) +
labs(title = str_wrap('Proportion of Elite and Non-Elite WR Seasons by Career Years', width=40), x = "Years in League", y = "Density", fill = NULL) +
scale_fill_manual(values = group_colors) +
scale_color_manual(values = group_colors) +
wr_theme
p
ggsave(filename = "~/Documents/elite_wr/elite_v_non-e_years.png", p, width = 12, height = 9, dpi = 300, bg = "white")
t.test(yrs_in_league ~ group, data = overlay_df)
##
## Welch Two Sample t-test
##
## data: yrs_in_league by group
## t = 3.779, df = 501.6, p-value = 0.0001764
## alternative hypothesis: true difference in means between group Elite and group Non-Elite is not equal to 0
## 95 percent confidence interval:
## 0.3252914 1.0297958
## sample estimates:
## mean in group Elite mean in group Non-Elite
## 4.203125 3.525581
#primary results: t = 3.6503, df = 494.75, p-value = 0.0002898
#the p-value is significant and rejects the null hypothesis that the mean years in the league between elite and non-elite WR seasons are equal
#now compare elite vs non-elite draft position
overlay_df <- bind_rows(
non_elite_wr_stats |> transmute(draft_pick = draft_pick, group = "Non-Elite"),
elite_wr_stats |> transmute(draft_pick = draft_pick, group = "Elite")
)
means <- overlay_df |>
group_by(group) |>
summarise(mean_val = mean(draft_pick, na.rm = TRUE), .groups = "drop") |>
mutate(label = paste0(group, " mean draft position: ", round(mean_val, 1)))
group_colors <- c("Elite" = "firebrick3", "Non-Elite" = "royalblue3")
p <- ggplot(overlay_df, aes(x = draft_pick, fill = group)) +
geom_histogram(aes(y = after_stat(density)), # density so unequal group sizes are comparable
position = "identity", alpha = 0.5, binwidth = 10, color = "white") +
geom_vline(data = means, aes(xintercept = mean_val, color = group),
linewidth = 1, linetype = "dashed", show.legend = FALSE) +
geom_label(data = means,
aes(x = mean_val, y = Inf, label = label, color = group),
vjust = c(1.5, 3.2), hjust = 0.03, # staggered so labels don't collide
fill = "white", show.legend = FALSE) +
labs(title = str_wrap('Proportion of Elite and Non-Elite WR Seasons by Draft Pick', width=40), x = "Draft Pick", y = "Density", fill = NULL) +
scale_fill_manual(values = group_colors) +
scale_color_manual(values = group_colors) +
wr_theme
p
ggsave(filename = "~/Documents/elite_wr/elite_v_non-e_pick.png", p, width = 12, height = 9, dpi = 300, bg = "white")
t.test(draft_pick ~ group, data = overlay_df)
##
## Welch Two Sample t-test
##
## data: draft_pick by group
## t = -8.5347, df = 505.41, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Elite and group Non-Elite is not equal to 0
## 95 percent confidence interval:
## -38.35059 -23.99812
## sample estimates:
## mean in group Elite mean in group Non-Elite
## 57.59375 88.76811
#pimrary results: t = -8.341, df = 498.54, p-value = 7.257e-16
#this is the most significant of the three, rejects the null hypothesis
#now compare elite vs non-elite draft position
overlay_df <- bind_rows(
non_elite_wr_stats |> transmute(draft_round = draft_round, group = "Non-Elite"),
elite_wr_stats |> transmute(draft_round = draft_round, group = "Elite")
)
means <- overlay_df |>
group_by(group) |>
summarise(mean_val = mean(draft_round, na.rm = TRUE), .groups = "drop") |>
mutate(label = paste0(group, " mean draft round: ", round(mean_val, 1)))
group_colors <- c("Elite" = "firebrick3", "Non-Elite" = "royalblue3")
p <- ggplot(overlay_df, aes(x = draft_round, fill = group)) +
geom_histogram(aes(y = after_stat(density)), # density so unequal group sizes are comparable
position = "identity", alpha = 0.5, binwidth = 1, color = "white") +
geom_vline(data = means, aes(xintercept = mean_val, color = group),
linewidth = 1, linetype = "dashed", show.legend = FALSE) +
geom_label(data = means,
aes(x = mean_val, y = Inf, label = label, color = group),
vjust = c(1.5, 3.2), hjust = -0.05, # staggered so labels don't collide
fill = "white", show.legend = FALSE) +
labs(title = str_wrap('Proportion of Elite and Non-Elite WR Seasons by Draft Round', width=40), x = "Draft Round", y = "Density", fill = NULL) +
scale_fill_manual(values = group_colors) +
scale_color_manual(values = group_colors) +
wr_theme
p
ggsave(filename = "~/Documents/elite_wr/elite_v_non-e_round.png", p, width = 12, height = 9, dpi = 300, bg = "white")
t.test(draft_round ~ group, data = overlay_df)
##
## Welch Two Sample t-test
##
## data: draft_round by group
## t = -8.7989, df = 507.15, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Elite and group Non-Elite is not equal to 0
## 95 percent confidence interval:
## -1.0783374 -0.6846817
## sample estimates:
## mean in group Elite mean in group Non-Elite
## 2.228125 3.109635
#Modeling and inference
#what is the probability of an elite season given draft position and years in the league
#making a df that will be used to build regression models of elite wrs
model_df <- wr_stats |>
#only using wrs that were draft, all the undrafted players skew the data
filter(drafted, !is.na(yrs_in_league), !is.na(draft_pick)) |>
mutate(
#adds binary logical variable for whether a wr is elite or not base on the definition we've been using
elite = as.integer(rec_yards >= 1000 | receptions >= 100),
#effect of draft pick isn't linear (as shown by the graphs of elite WR by draft pick/round)
#massive drop off afer round 1, then after round 2/3, flattens in later rounds; difference between pick 1 and 10 is important, difference between pick 210 and 220 is negligible
#will be doing LOGISTIC regression as a result of this and compare to a standard linear regression
log_pick = log(draft_pick)
)
#Linear version, regressing elite (TRUE/FALSE) on yrs_in_league, and draft_pick (family=binomial b/c elite is a binomial variable)
m1 <- glm(elite ~ yrs_in_league + draft_pick, data = model_df, family = binomial)
#Better version: career arcs are curved, and pick value is steep early then flat
#I(yrs_in_league^2) creates a squared term for years in the league, fits a curved relationship (WR careers usually follow a parabolic curve, not a straight line)
m2 <- glm(elite ~ yrs_in_league + I(yrs_in_league^2) + log_pick,
data = model_df, family = binomial)
summary(m2)
##
## Call:
## glm(formula = elite ~ yrs_in_league + I(yrs_in_league^2) + log_pick,
## family = binomial, data = model_df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.09209 0.25639 -0.359 0.719
## yrs_in_league 0.37513 0.06556 5.722 1.05e-08 ***
## I(yrs_in_league^2) -0.03227 0.00622 -5.187 2.14e-07 ***
## log_pick -0.55743 0.05760 -9.677 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1694.5 on 1824 degrees of freedom
## Residual deviance: 1560.4 on 1821 degrees of freedom
## AIC: 1568.4
##
## Number of Fisher Scoring iterations: 5
#yrs_in_league has positive correlation (.3642), and I(yrs_in_league^2) has a slight negative correlation (-.03), this means the curve of this is an inverted U (what we expect)
#log_pick correlation is -.539, doubling of draft pick number multiplies odds of an elite season: exp(-0.538 * log(2)) = exp(-0.373) ≈ 0.69; double the draft pick means likelihood of elite season is cut by 30% (from the original draft position, so if pick 50 is 15%, then pick 100 is 10%)
#each predictor is significant
out <- capture.output(summary(m2))
txt <- paste(out, collapse = "\n")
p <- ggplot() +
annotate("text", x = 0, y = 1, label = txt, hjust = 0, vjust = 1,
family = "mono", size = 3.5, color = "#d4d4d4") +
xlim(0, 1) + ylim(0, 1) +
theme_void() +
theme(plot.background = element_rect(fill = "#1e1e1e", color = NA),
plot.margin = margin(20, 20, 20, 20))
ggsave(filename = "~/Documents/elite_wr/m2_summary.png", p, width = 8, height = 6, dpi = 300)
AIC(m1, m2) # lower is better
#using the second model to predict an elite season from a 3rd year wr that was drafted 15 overall
predict(m2,
newdata = data.frame(yrs_in_league = 3, log_pick = log(15)),
type = "response") #type = "response" makes the output a normal probability instead of a log-odds ('link' returns the log-odds from glm)
## 1
## 0.3172072
#expand.grid() makes a df of every combination of values you give
#in this case, every combination of yrs_in_league (0 to 12) and draft_pick(5, 15, 32, 64, 100, 150, 220), then for each combination of year in the league and draft pick we calculate probability of an elite season with the model
grid <- expand.grid(yrs_in_league = 0:12, draft_pick = c(5, 15, 32, 64, 100, 150, 220)) |>
mutate(log_pick = log(draft_pick))
#add probability of an elite season to the grid df
grid$p_elite <- predict(m2, grid, type = "response")
p <- ggplot(grid, aes(yrs_in_league, p_elite, color = factor(draft_pick))) +
geom_line(linewidth = 1) +
labs(title = str_wrap("Probability of Elite WR Season \nGiven Year in League and Draft Pick", width=40), x = "Years in league", y = "P(elite season)", color = "Draft pick") +
wr_theme
p
ggsave(filename = "~/Documents/elite_wr/elite_prob_m1.png", p, width = 12, height = 9, dpi = 300, bg = "white")
roc(model_df$elite, fitted(m2)) # AUC; 0.5 = coin flip, 1 = perfect
##
## Call:
## roc.default(response = model_df$elite, predictor = fitted(m2))
##
## Data: fitted(m2) in 1505 controls (model_df$elite 0) < 320 cases (model_df$elite 1).
## Area under the curve: 0.7005
#a random elite season gets a higher predicted probability than a random non-elite season .6955 of the time
#new model adds INTERACTION between draft position and career curve; shape of the years-in-league curve can differ by draft pick
#in m2 every player shares one curve, one years where career peaks the most across all draft picks (about year 6)
#earlier draft picks may get to peak earlier in career, later picks may take longer
m4 <- glm(elite ~ (yrs_in_league + I(yrs_in_league^2)) * log_pick,
data = model_df, family = binomial)
summary(m4)
##
## Call:
## glm(formula = elite ~ (yrs_in_league + I(yrs_in_league^2)) *
## log_pick, family = binomial, data = model_df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.661e+00 4.969e-01 3.343 0.000828 ***
## yrs_in_league -2.307e-01 1.999e-01 -1.154 0.248438
## I(yrs_in_league^2) 8.864e-05 1.715e-02 0.005 0.995875
## log_pick -1.048e+00 1.385e-01 -7.566 3.86e-14 ***
## yrs_in_league:log_pick 1.677e-01 5.572e-02 3.010 0.002616 **
## I(yrs_in_league^2):log_pick -8.737e-03 4.854e-03 -1.800 0.071836 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1694.5 on 1824 degrees of freedom
## Residual deviance: 1541.3 on 1819 degrees of freedom
## AIC: 1553.3
##
## Number of Fisher Scoring iterations: 5
out <- capture.output(summary(m4))
txt <- paste(out, collapse = "\n")
p <- ggplot() +
annotate("text", x = 0, y = 1, label = txt, hjust = 0, vjust = 1,
family = "mono", size = 3.5, color = "#d4d4d4") +
xlim(0, 1) + ylim(0, 1) +
theme_void() +
theme(plot.background = element_rect(fill = "#1e1e1e", color = NA),
plot.margin = margin(20, 20, 20, 20))
ggsave(filename = "~/Documents/elite_wr/m4_summary.png", p, width = 8, height = 6, dpi = 300)
grid$p_int <- predict(m4, grid, type = "response")
grid |>
mutate(tier = cut(draft_pick, c(0, 10, 32, 64, 100, 256))) |>
group_by(tier) |>
slice_max(p_int, n = 1, with_ties = FALSE)
p <- ggplot(grid, aes(yrs_in_league, p_int, color = factor(draft_pick))) +
geom_line(linewidth = 1) +
labs(title = str_wrap("Probability of Elite WR Season \nGiven Year in League and Draft Pick (Interaction Model)", width=40), x = "Years in league", y = "P(elite season)", color = "Draft pick") +
wr_theme
ggsave(filename = "~/Documents/elite_wr/elite_prob_m2.png", p, width = 12, height = 9, dpi = 300, bg = "white")
p
#see which model performs better
AIC(m2, m4)
anova(m2, m4, test = "Chisq")
#m2 residual deviance: 1589
#m4 residual deviance: 1571
#large residual deviance drop-off
#rows of same players correlate with each other, can create bias
#m4 fits better than m2
#heat map of probability for an elite season from m2
heat_grid <- expand.grid(
yrs_in_league = 0:12,
draft_pick = 1:256
) |>
mutate(log_pick = log(draft_pick))
heat_grid$p_elite <- predict(m2, heat_grid, type = "response")
p <- ggplot(heat_grid, aes(yrs_in_league, draft_pick, fill = p_elite)) +
geom_tile() +
scale_y_reverse(breaks = seq(0, 250, 25)) +
scale_x_continuous(breaks = 0:12) +
scale_fill_gradientn(
colors = c("#08306b", "#4393c3", "#f7f7f7", "#f4a582", "#b2182b"),
labels = scales::percent, name = "P(elite)"
) +
labs(title = 'Probability of Elite Season by Draft Pick', x = "Years in league", y = "Draft pick") +
wr_theme
p
ggsave(filename = "~/Documents/elite_wr/m1_heatmap.png", p, width = 12, height = 9, dpi = 300, bg = "white")
#heat map of probability for an elite season from m4
heat_grid <- expand.grid(
yrs_in_league = 0:12,
draft_pick = 1:256
) |>
mutate(log_pick = log(draft_pick))
heat_grid$p_elite <- predict(m4, heat_grid, type = "response")
p <- ggplot(heat_grid, aes(yrs_in_league, draft_pick, fill = p_elite)) +
geom_tile() +
scale_y_reverse(breaks = seq(0, 250, 25)) +
scale_x_continuous(breaks = 0:12) +
scale_fill_gradientn(
colors = c("#08306b", "#4393c3", "#f7f7f7", "#f4a582", "#b2182b"),
labels = scales::percent, name = "P(elite)"
) +
labs(title=str_wrap('Probability of Elite Season by Draft Pick (Interaction Model)', width=40), x = "Years in league", y = "Draft pick") +
wr_theme
p
ggsave(filename = "~/Documents/elite_wr/m2_heatmap.png", p, width = 12, height = 9, dpi = 300, bg = "white")
model_df |> group_by(yrs_in_league) |> summarise(rate = mean(elite), n = n())
#raw counts of elite seasons for each draft pick and year
#based on eye test, looks like it agrees more with m4 than m2
elite_counts <- model_df |>
filter(elite == 1) |>
count(yrs_in_league, draft_pick)
p <- ggplot(elite_counts, aes(yrs_in_league, draft_pick, fill = n)) +
geom_tile() +
scale_y_reverse(breaks = seq(0, 250, 25)) +
scale_x_continuous(breaks = 0:15) +
scale_fill_gradientn(
colors = c("#f4a582", "#b2182b"),
name = "Elite seasons",
breaks = scales::breaks_pretty()
) +
labs(title='Counts of Elite Seasons by Draft Pick', x = "Years in league", y = "Draft pick") +
wr_theme
p
ggsave(filename = "~/Documents/elite_wr/counts_heatmap.png", p, width = 12, height = 9, dpi = 300, bg = "white")