MLB

Since the Pythagorean winning percentage is used most commonly in baseball, I wanted to look at baseball first. As previously mentioned, the ideal exponent has been found to be 1.83, so I will compare the result I get from here to that exponent.

The data to calculate the exponent comes from the Lahman package, which contains baseball statistics all the way from 1871 through the present. Since some statistics were recorded early on, any year after 1900 was used. The run differential and winning percentage was calculated for each team every season. Linear regression was then used to calculate what the value of the exponent would be. Note that from this formula, there is no intercept in the model, which is why -1 is included in the call to the lm function.

The resulting value given is 1.854, which is close to the ideal value of 1.83 and can be derived simply using only linear regression.

library(tidyverse)
library(Lahman)
mlb <- Teams %>% 
    filter(yearID > 1900) %>%
    select(teamID, yearID, lgID, G, W, L, R, RA) %>%
    mutate(RD = R - RA, 
           Wpct = W / (W + L))

mlb_pyt <- mlb %>%
    mutate(log_w_ratio = log(W/L),
           log_r_ratio = log(R / RA))
pyt_fit <- lm(log_w_ratio ~ -1 + log_r_ratio, data = mlb_pyt)
round(coef(pyt_fit), digits = 3)

## log_r_ratio 
##       1.854

NHL

For NHL data, I used the nhlapi package, which uses the open NHL API to get data. The standings for all seasons from 1990 through 2022 were pulled, where the point differential could be calculated using the goals scored and the goals against.

Linear regression is used again in the same way that is was for the MLB data and the resulting value for the exponent is 2.35. A study of Pythagorean winning percentage by Kevin D. Dayaratna and Steven J. Miller for the NHL found an ideal exponent of just above 2.²

library(nhlapi)

nhl <- nhl_standings(seasons = 1990:2022)
nhl_data <- nhl$teamRecords %>%
    bind_rows() %>%
    mutate(point_diff = goalsScored - goalsAgainst,
           Win_pct = leagueRecord.wins/leagueRecord.losses,
           log_win_ratio = log(leagueRecord.wins/leagueRecord.losses),
           log_point_ratio = log(goalsScored/goalsAgainst))
pyt_fit_nhl <- lm(log_win_ratio ~ -1 + log_point_ratio, data = nhl_data)
coef(pyt_fit_nhl)

## log_point_ratio 
##        2.354845

NFL

The NFL data was scraped using the espnscrapeR package. The package can be used to scrape the NFL standings for each year. The earliest year that provided the points for and against a team was 2002 so the data for the model goes from 2002-2021.

One adjustment I had to make for the NFL data was I had to use log1p() rather than log() for the wins and losses of each team. Since there have been a few instances of a team winning or losing zero games, the natural logarithm would have been undefined. log1p() calculates \(log(1 + x)\), which adjusts those values that would have been zero. There may be other ways to handle those specific seasons, but log1p() was a quick way that I thought of to cover that. I also used map_df() to get the standings for each season and return a data frame without using a for loop. The map functions in the purrr package are very useful in a variety of scenarios.

The output from the linear regression model gives an exponent of 2.396. Football Outsiders has previously used an exponent of 2.37, which they credited to Daryl Morey. They had also suggested a new exponent that adjusts for offensive environment.³.

library(espnscrapeR)
library(purrr)

nfl_data <- map_df(2002:2021, get_nfl_standings)

nfl_pyt <- nfl_data %>%
    mutate(log_w_ratio = log1p(wins) - log1p(as.numeric(losses)),
           log_point_ratio = log(pts_for/pts_against))

pyt_fit_nfl <- lm(log_w_ratio ~ -1 + log_point_ratio, data = nfl_pyt)
coef(pyt_fit_nfl)

## log_point_ratio 
##        2.396146

NBA

The NBA data was scraped using the hoopr package, which also uses ESPN to scrape the yearly standings for teams. Data was taken from the 1994-1995 season through the 2021-2022 season. The average points for and against each team were givent in the standings, so the total number of points scored and allowed for the season were calculated by multiplying by the total number of games in the season (this is generally 82, but there are a few lockout seasons that were shortened and the COVID-19 shortened season as well).

The exponent given from linear regression here is 14.085. Daryl Morey is credited with deriving the exponent for basketball and he got an exponent of 13.91.

library(hoopR)
nba_data <- map_df(1994:2021, espn_nba_standings)
nba_data <- nba_data %>%
    mutate(wins = as.numeric(wins),
           losses = as.numeric(losses),
           avgpointsfor = as.numeric(avgpointsfor),
           avgpointsagainst = as.numeric(avgpointsagainst),
           games = wins + losses,
           points_for = games*avgpointsfor,
           points_against = games*avgpointsagainst,
           point_diff = points_for - points_against,
           log_w_ratio = log(wins/losses),
           log_point_ratio = log(points_for/points_against))
pyt_fit_nba <- lm(log_w_ratio ~ -1 + log_point_ratio, data = nba_data)
coef(pyt_fit_nba)

## log_point_ratio 
##        14.08544