Today we’re going to bootstrap NFL games and learn about Pythagorean Expectation.
Let’s start by loading the data, our favorite packages, and taking a look:
library(dplyr)
library(ggplot2)
nfl <- read.csv('/home/rstudioshared/shared_files/data/nflgames2016.csv')
head(nfl)We can calculate 2016 standings using dplyr:
standings <- nfl %>% group_by(team) %>%
summarize(pts = sum(points), opp_pts = sum(opp_points),
wins = sum(win), losses = sum(loss), ties=sum(tie),
wp=(wins+0.5*ties)/16)
head(standings)Bill James devised a formula to calculate how often a team would be expected to win (Win%) based on the the numbers of rusn they scored (RS) and the number or runs they allowed (RA).
\[Win\ \% = \frac{RS^2}{RS^2 + RA^2}\]
This formula can be extended to work in other sports by simply changing the exponent. In the NFL, the following formula will work:
\[Win\ \% = \frac{points\ scored^{2.37}}{points\ scored^{2.37} + points\ allowed^{2.37}}\]
Let’s add “Pythagorean” wins, losses, and winning percentage to our standings:
xp <- 2.37
standings <- standings %>% mutate(pythag_wins = 16*pts^xp/(pts^xp + opp_pts^xp),
pythag_losses = 16*opp_pts^xp/(pts^xp + opp_pts^xp),
pythag_wp = pts^xp/(pts^xp + opp_pts^xp)
) Now, let’s take a good look at how the Pythagorean standings and the real standings compare:
View(standings)We can compare teams pythagorean winning percentages to their actual winning percentages graphically:
ggplot(standings, aes(pythag_wp, wp))+geom_point() + geom_smooth(method="lm")We might also be wondering, how confident we can be that the teams that performed the best, are the best. Let’s ltry to answer this question by boostrapping. Here’s one alternate reality created by sampling 16 games, with replacement for each team.
nfl.sample <- nfl %>% group_by(team) %>% sample_n(size = 16, replace=TRUE)
sample.standings <- nfl.sample %>% group_by(team) %>%
summarize(pts = sum(points), opp_pts = sum(opp_points),
wins = sum(win), losses = sum(loss), ties=sum(tie), wp=(wins+0.5*ties)/16) %>%
mutate(pythag_wins = 16*pts^xp/(pts^xp + opp_pts^xp),
pythag_losses = 16*opp_pts^xp/(pts^xp + opp_pts^xp),
pythag_wp = pts^xp/(pts^xp + opp_pts^xp)
)
sample.standings %>% arrange(desc(pythag_wp))Just as in our previously lab, we don’t want to create only one alternate reality but many. Let’s create 100 bootstrapped samples of the season. First we’ll create a place holder, sim.results, and then we’ll use a for loop to repeat our sampling code 100 times, each time adding our results to the sim.results table.
sim.results <- standings[FALSE, ]
for (i in 1:100){
nfl.sample <- nfl %>% group_by(team) %>% sample_n(size = 16, replace=TRUE)
sample.standings <- nfl.sample %>% group_by(team) %>%
summarize(pts = sum(points), opp_pts = sum(opp_points),
wins = sum(win), losses = sum(loss), ties=sum(tie), wp=(wins+0.5*ties)/16) %>%
mutate(pythag_wins = 16*pts^xp/(pts^xp + opp_pts^xp),
pythag_losses = 16*opp_pts^xp/(pts^xp + opp_pts^xp),
pythag_wp = pts^xp/(pts^xp + opp_pts^xp)
)
sim.results <- rbind(sim.results, sample.standings)
}How often did teams have perfect seasons and which teams accomplished this?
sim.results %>% filter(wins==16)How many teams went winless in at least one simulation and how often did they go winless?
sim.results %>% filter(wins==0) %>% group_by(team) %>% summarize(times.winless = length(team))Finally, let’s calculate the best and worst seasons for each team and the standard deviations in both winning percentage and Pythagorean winning percentages:
sim.summary <- sim.results %>% group_by(team) %>% summarize(
min.wins = min(wins), max.wins=max(wins), sd(wins), sd(pythag_wins)
)View(sim.summary)How wide are the distributions of outcomes for each team? How do the standard deviations in winning percentage compare to the standard deviations in pythagorean winning percentage?