MLB Contract Valuation analysis and discovering underrated traits
STAT 420, Kyungmo Eric Ryu
Abstract
Baseball has provided a wealth of statistical data that many have
analyzed over the years to answer various questions. There’s a lot of
specific questions one can answer with this data, such as “when
was the last time a pitcher who allowed 5 home runs in a game and still
won?”, to more general questions like how good is a player overall.
Such a general question has been addressed by various metrics such as
the Wins Above Replacement, or WAR, of which countless have
tried using to predict how players will perform. Indeed, the baseball
industry is the forerunner in modern sports data analysis which is well
demonstrated in the renowned book and film Moneyball. While
there have been many studies done and many debates waged over more
general statistics as to how they should be calculated and interpreted,
there is so much depth in the amount of data in baseball that there are
still ways to conduct a creative study. Among many possible explorable
areas, this study would like to accomplish at least two goals. First, to
find a model that predicts player success based on raw batted ball data.
A ‘hit’ in baseball depends on many factors-most notably what direction
and launch angle the ball was hit, and even a bit of luck as it depends
on whether the opponent fielder is able to field the ball. In this
exercise, we would like to compare a player’s BABIP(Batting Average on
Balls In Play) and see how it compares to his raw batted ball data such
as Hard hit %, Pull %, Line drive %, etc. BABIP does not tell the full
story of a hitter’s performance as a single isn’t an equal performance
to a home run, but it is the best response variable that shows the ratio
of successful plate appearances to all plate appearances. Secondly, this
study attempts to examine how player success and predictors of success
compare to player salaries.
Data information
There are a number of major websites that provide baseball statistics which include: Baseball Prospectus, Baseball Reference, and Fangraphs. This study will focus mostly on Fangraphs when performing our analyis but will also make use of Baseball Reference as they have contract information for players which is not available on Fangraphs.
Fangraphs provides a large amount of variables, including most typical baseball statistical values as well as more unique and specific values and probabilities. Of particular interest are the three true outcomes (strikeouts, home runs, and walks) and batted ball data (line drive%, hard hit%, etc). Please see the Appendix for definitions of terms in datatables.
Batter data was collected from Fangraphs’s Split Leaderboard tool for
viewing batter statistics. Only those with PA greater than
150 were collected for the 2018 and 2019 seasons. This is about 150
players for each season.
Salary data was collected from Baseball-Reference via a python web
scraper script for the 2018 and 2019 seasons. Further python scripts
were made to refine and combine the salary data with the data collected
from Fangraphs. These aggregate datasets are contained in
Combined-2018.csv and Combined-2019.csv. See
appendix for samples scripts used for the 2019 season.
Methods
BABIP Investigations
We will use the 2018 fangraphs data as the training set, and the 2019 data as the testing set. Some columns such as player name and counting statistics with low variability are deleted for a more efficient analysis.
#import and refine data
fangraphs_2019 = read.csv("Combined-2019.csv")
fangraphs_2018 = read.csv("Combined-2018.csv")
fangraphs_2018 = subset(fangraphs_2018, select = -c(BUH, BUH., IFH, Name, Team, playerid))
fangraphs_2019 = subset(fangraphs_2019, select = -c(s.age, s.lifespan, Name, Team, playerid))Now we want to fit a mulitiple linear regression model that predicts a player’s BABIP(Batting Average on Balls In Play) and ultimately their performance from it. BABIP is defined as : \[BABIP = (H-HR)/(AB-K-HR+SF)\] We want to find a model that predicts BABIP from raw batted ball data.
full_batted_mod = lm(BABIP ~ GB.FB+LD.+GB.+FB.+IFFB.+HR.FB+IFH.+Pull.+Cent.+Oppo.+Soft.+Med.+Hard.+Swing.+Contact., data = fangraphs_2018)
#AIC search
aic_trn = step(full_batted_mod, direction = "both", trace = FALSE)
#BIC search
n = nrow(fangraphs_2018)
bic_trn = step(full_batted_mod, direction = "both", k = log(n), trace = FALSE)
summary(bic_trn)$coef## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.011522 0.0291637 0.3951 6.934e-01
## LD. 0.005703 0.0006832 8.3483 7.277e-14
## GB. 0.001851 0.0003363 5.5047 1.799e-07
## IFH. 0.005005 0.0007828 6.3939 2.443e-09
## Hard. 0.001499 0.0003786 3.9587 1.213e-04babip_mod = lm(BABIP ~ LD. + GB. + IFH. + Hard., data = fangraphs_2019)
war_mod = lm(WAR ~ BABIP + BB. + K. + BsR + Spd, data = fangraphs_2019)
wpa_mod = lm(WPA ~ BABIP + BB. + K. + BsR + Spd, data = fangraphs_2019)Salary Investigations
All tests conducted assume a \(\alpha = 0.05\).
Predicting Pay Level Based on Performance Metrics
# extract only salaried players from 2018 and 2019 data.
bstats2018 = read.csv("Combined-2018.csv", stringsAsFactors = FALSE, strip.white = TRUE)
bstats2019 = read.csv("Combined-2019.csv", stringsAsFactors = FALSE, strip.white = TRUE)
bstats2018 = subset(bstats2018, agreement == "CS")
bstats2019 = subset(bstats2019, agreement == "CS")
bstats2018$salary = as.numeric(bstats2018$salary)/1000000
bstats2019$salary = as.numeric(bstats2019$salary)/1000000
bstats2019$s.age = as.numeric(bstats2019$s.age)
bstats2019$s.lifespan = as.numeric(bstats2019$s.lifespan)
# Add column for Pay Level (Low = <5 mil, High = >5 mil, the median of the salaries)
bstats2018 = transform(bstats2018, pay.level=ifelse(salary<median(bstats2018$salary), 0, 1))
bstats2019 = transform(bstats2019, pay.level=ifelse(salary<median(bstats2019$salary), 0, 1))
# select only parameters we care about
bstats2018 = subset(bstats2018, select = c("salary", "WAR", "WPA", "BsR", "Hard.", "LD.", "BABIP", "pay.level"))
bstats2019 = subset(bstats2019, select = c("salary", "WAR", "WPA", "BsR", "Hard.", "LD.", "BABIP", "pay.level"))
# Run pairs to see if any correlations between salary and other values
pairs(bstats2018)
# Train Classifiers on 2018 data
nilmodel = glm(pay.level ~ 1, data = bstats2018, family = binomial)
model_add = glm(pay.level ~ WAR + WPA + BsR + Hard. + LD. + BABIP, data = bstats2018, family = binomial)
model_add_log_hard = glm(pay.level ~ WAR + WPA + BsR + log(Hard.) + LD. + BABIP, data = bstats2018, family = binomial)
model_int = glm(pay.level ~ WAR * WPA * BsR * Hard. * LD. * BABIP, data = bstats2018, family = binomial)
model_int_log_hard = glm(pay.level ~ WAR * WPA * BsR * log(Hard.) * LD. * BABIP, data = bstats2018, family = binomial)
forward_aic_add = step(nilmodel, scope = pay.level ~ WAR + WPA + BsR + Hard. + LD. + BABIP, direction="forward", trace = FALSE)
forward_aic_add_log_hard = step(nilmodel, scope = pay.level ~ WAR + WPA + BsR + log(Hard.) + LD. + BABIP, direction="forward", trace = FALSE)
summary(forward_aic_add)##
## Call:
## glm(formula = pay.level ~ Hard. + BsR, family = binomial, data = bstats2018)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.6611 1.2280 -2.17 0.030 *
## Hard. 0.0723 0.0323 2.24 0.025 *
## BsR -0.0854 0.0506 -1.69 0.091 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 181.60 on 130 degrees of freedom
## Residual deviance: 172.15 on 128 degrees of freedom
## AIC: 178.2
##
## Number of Fisher Scoring iterations: 4summary(forward_aic_add_log_hard)##
## Call:
## glm(formula = pay.level ~ log(Hard.) + BsR, family = binomial,
## data = bstats2018)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.1426 4.2209 -2.17 0.030 *
## log(Hard.) 2.5450 1.1653 2.18 0.029 *
## BsR -0.0846 0.0505 -1.67 0.094 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 181.60 on 130 degrees of freedom
## Residual deviance: 172.32 on 128 degrees of freedom
## AIC: 178.3
##
## Number of Fisher Scoring iterations: 4model_hard = glm(pay.level ~ Hard., data = bstats2018, family = binomial)
model_log_hard = glm(pay.level ~ log(Hard.), data = bstats2018, family = binomial)
model_bsr = glm(pay.level ~ BsR, data = bstats2018, family = binomial)
model_int_hard_bsr = glm(pay.level ~ Hard. * BsR, data = bstats2018, family = binomial)
model_war = glm(pay.level ~ WAR, data = bstats2018, family = binomial)
model_wpa = glm(pay.level ~ WPA, data = bstats2018, family = binomial)
model_babip = glm(pay.level ~ BABIP, data = bstats2018, family = binomial)
# check if BsR is actually helpful
anova( model_hard, forward_aic_add, test="LRT")## Analysis of Deviance Table
##
## Model 1: pay.level ~ Hard.
## Model 2: pay.level ~ Hard. + BsR
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 129 175
## 2 128 172 1 2.95 0.086 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# check if Hard. is actually helpful
anova( model_bsr, forward_aic_add, test="LRT")## Analysis of Deviance Table
##
## Model 1: pay.level ~ BsR
## Model 2: pay.level ~ Hard. + BsR
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 129 178
## 2 128 172 1 5.35 0.021 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# check if log(Hard.) is actually helpful
anova( forward_aic_add, forward_aic_add_log_hard, test="LRT")## Analysis of Deviance Table
##
## Model 1: pay.level ~ Hard. + BsR
## Model 2: pay.level ~ log(Hard.) + BsR
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 128 172
## 2 128 172 0 -0.167# check if aic model is actually better than overfitted additive model
anova( forward_aic_add, model_add, test="LRT")## Analysis of Deviance Table
##
## Model 1: pay.level ~ Hard. + BsR
## Model 2: pay.level ~ WAR + WPA + BsR + Hard. + LD. + BABIP
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 128 172
## 2 124 171 4 1.4 0.84# check if interaction between Hard. and BsR is helpful
anova( forward_aic_add, model_int_hard_bsr, test="LRT")## Analysis of Deviance Table
##
## Model 1: pay.level ~ Hard. + BsR
## Model 2: pay.level ~ Hard. * BsR
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 128 172
## 2 127 172 1 0.203 0.65# compare best model to a WAR model
anova( model_war, forward_aic_add, test="LRT")## Analysis of Deviance Table
##
## Model 1: pay.level ~ WAR
## Model 2: pay.level ~ Hard. + BsR
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 129 182
## 2 128 172 1 9.33 0.0022 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# compare best model to a WPA model
anova( model_wpa, forward_aic_add, test="LRT")## Analysis of Deviance Table
##
## Model 1: pay.level ~ WPA
## Model 2: pay.level ~ Hard. + BsR
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 129 179
## 2 128 172 1 7.22 0.0072 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# compare best model to a WAR model
anova( model_babip, forward_aic_add, test="LRT")## Analysis of Deviance Table
##
## Model 1: pay.level ~ BABIP
## Model 2: pay.level ~ Hard. + BsR
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 129 181
## 2 128 172 1 9.06 0.0026 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test classifiers on 2019 data
hard_pred = ifelse(predict(model_hard, bstats2019) > 0, "Low", "High")
bsr_pred = ifelse(predict(model_bsr, bstats2019) > 0, "Low", "High")
war_pred = ifelse(predict(model_war, bstats2019) > 0, "Low", "High")
wpa_pred = ifelse(predict(model_wpa, bstats2019) > 0, "Low", "High")
babip_pred = ifelse(predict(model_babip, bstats2019) > 0, "Low", "High")
aic_pred = ifelse(predict(forward_aic_add, bstats2019) > 0, "Low", "High")
aiclog_pred = ifelse(predict(forward_aic_add_log_hard, bstats2019) > 0, "Low", "High")
conf_mat_hard = make_conf_mat(predicted=hard_pred, actual=bstats2019$pay.level)
conf_mat_bsr = make_conf_mat(predicted=bsr_pred, actual=bstats2019$pay.level)
conf_mat_war = make_conf_mat(predicted=war_pred, actual=bstats2019$pay.level)
conf_mat_wpa = make_conf_mat(predicted=wpa_pred, actual=bstats2019$pay.level)
conf_mat_babip = make_conf_mat(predicted=babip_pred, actual=bstats2019$pay.level)
conf_mat_aic = make_conf_mat(predicted=aic_pred, actual=bstats2019$pay.level)
conf_mat_aiclog = make_conf_mat(predicted=aiclog_pred, actual=bstats2019$pay.level)
hard_mis = mean(hard_pred != bstats2019$pay.level)
bsr_mis = mean(bsr_pred != bstats2019$pay.level)
war_mis = mean(war_pred != bstats2019$pay.level)
wpa_mis = mean(wpa_pred != bstats2019$pay.level)
babip_mis = mean(babip_pred != bstats2019$pay.level)
aic_mis = mean(aic_pred != bstats2019$pay.level)
aiclog_mis = mean(aiclog_pred != bstats2019$pay.level)
# Analyze Classifiers (sensitivity and specificity)
sens_hard = get_sens(conf_mat_hard)
sens_bsr = get_sens(conf_mat_bsr)
sens_war = get_sens(conf_mat_war)
sens_wpa = get_sens(conf_mat_wpa)
sens_babip = get_sens(conf_mat_babip)
sens_aic = get_sens(conf_mat_aic)
sens_aiclog = get_sens(conf_mat_aiclog)
spec_hard = get_spec(conf_mat_hard)
spec_bsr = get_spec(conf_mat_bsr)
spec_war = get_spec(conf_mat_war)
spec_wpa = get_spec(conf_mat_wpa)
spec_babip = get_spec(conf_mat_babip)
spec_aic = get_spec(conf_mat_aic)
spec_aiclog = get_spec(conf_mat_aiclog)Predicting Salary Lifespan Based on Performance Metrics
# extract only salaried players from 2019 data.
bstats2019 = read.csv("Combined-2019.csv", stringsAsFactors = FALSE, strip.white = TRUE)
bstats2019 = subset(bstats2019, agreement == "CS")
bstats2019$salary = as.numeric(bstats2019$salary)/1000000
bstats2019$s.age = as.numeric(bstats2019$s.age)
bstats2019$s.lifespan = as.numeric(bstats2019$s.lifespan)
# Add column for Salary Lifespan Type (Short = <=1 year, Long = >1)
bstats2019 = transform(bstats2019, s.lifespan.type=ifelse(s.lifespan<=1, 0, 1))
# Divide dataset into two equal halves randomly
set.seed(42)
trn_idx = sample(nrow(bstats2019), nrow(bstats2019)/2)
bstats_trn = bstats2019[trn_idx, ]
bstats_tst = bstats2019[-trn_idx, ]
# select only parameters we care about
bstats_trn = subset(bstats2019, select = c("salary", "WAR", "WPA", "BsR", "Hard.", "LD.", "BABIP", "s.lifespan.type"))
bstats_tst = subset(bstats2019, select = c("salary", "WAR", "WPA", "BsR", "Hard.", "LD.", "BABIP", "s.lifespan.type"))
# Train Classifiers on training data
nilmodel = glm(s.lifespan.type ~ 1, data = bstats_trn, family = binomial)
model_add = glm(s.lifespan.type ~ WAR + WPA + BsR + Hard. + LD. + BABIP, data = bstats_trn, family = binomial)
model_add_log_hard = glm(s.lifespan.type ~ WAR + WPA + BsR + log(Hard.) + LD. + BABIP, data = bstats_trn, family = binomial)
model_int = glm(s.lifespan.type ~ WAR * WPA * BsR * Hard. * LD. * BABIP, data = bstats_trn, family = binomial)
model_int_log_hard = glm(s.lifespan.type ~ WAR * WPA * BsR * log(Hard.) * LD. * BABIP, data = bstats_trn, family = binomial)
n = length(bstats_trn$salary)
backward_bic_add = step(model_int, direction="backward", trace = FALSE, k = log(n))
forward_aic_add_squares = step(nilmodel, scope = s.lifespan.type ~ poly(WAR, 2, raw = TRUE)
+ poly(WPA, 2, raw = TRUE)
+ poly(BsR, 2, raw = TRUE)
+ poly(Hard., 2, raw = TRUE)
+ poly(LD. , 2, raw = TRUE)
+ poly(BABIP, 2, raw = TRUE), direction="forward", trace = FALSE)
forward_aic_add = step(nilmodel, scope = s.lifespan.type ~ WAR + WPA + BsR + Hard. + LD. + BABIP, direction="forward", trace = FALSE)
summary(forward_aic_add)##
## Call:
## glm(formula = s.lifespan.type ~ Hard., family = binomial, data = bstats_trn)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.4065 1.2529 -1.92 0.055 .
## Hard. 0.0494 0.0308 1.60 0.109
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 197.20 on 146 degrees of freedom
## Residual deviance: 194.55 on 145 degrees of freedom
## AIC: 198.5
##
## Number of Fisher Scoring iterations: 4summary(backward_bic_add)##
## Call:
## glm(formula = s.lifespan.type ~ WAR * WPA * BsR * Hard. * LD. *
## BABIP, family = binomial, data = bstats_trn)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.08e+17 1.41e+10 -35982926 <2e-16 ***
## WAR 2.43e+17 8.48e+09 28643375 <2e-16 ***
## WPA -3.96e+16 1.04e+10 -3819615 <2e-16 ***
## BsR -2.17e+17 1.02e+10 -21356952 <2e-16 ***
## Hard. 1.18e+16 3.52e+08 33583869 <2e-16 ***
## LD. 2.38e+16 6.87e+08 34637182 <2e-16 ***
## BABIP 1.97e+18 4.85e+10 40688250 <2e-16 ***
## WAR:WPA 6.20e+16 4.55e+09 13627226 <2e-16 ***
## WAR:BsR 6.73e+16 4.76e+09 14122927 <2e-16 ***
## WPA:BsR 2.36e+17 7.62e+09 31037694 <2e-16 ***
## WAR:Hard. -5.09e+15 2.14e+08 -23798658 <2e-16 ***
## WPA:Hard. -1.37e+15 2.53e+08 -5428367 <2e-16 ***
## BsR:Hard. 4.15e+15 2.48e+08 16728621 <2e-16 ***
## WAR:LD. -9.51e+15 4.10e+08 -23199025 <2e-16 ***
## WPA:LD. 4.90e+14 4.83e+08 1014677 <2e-16 ***
## BsR:LD. 1.06e+16 4.86e+08 21772536 <2e-16 ***
## Hard.:LD. -5.56e+14 1.71e+07 -32607207 <2e-16 ***
## WAR:BABIP -9.54e+17 2.81e+10 -33903102 <2e-16 ***
## WPA:BABIP 2.62e+17 3.39e+10 7728112 <2e-16 ***
## BsR:BABIP 7.24e+17 3.40e+10 21260243 <2e-16 ***
## Hard.:BABIP -4.68e+16 1.21e+09 -38661083 <2e-16 ***
## LD.:BABIP -9.29e+16 2.34e+09 -39746255 <2e-16 ***
## WAR:WPA:BsR -7.01e+16 2.75e+09 -25519877 <2e-16 ***
## WAR:WPA:Hard. -8.37e+14 1.06e+08 -7878924 <2e-16 ***
## WAR:BsR:Hard. -1.08e+15 1.15e+08 -9399000 <2e-16 ***
## WPA:BsR:Hard. -5.31e+15 1.90e+08 -27922280 <2e-16 ***
## WAR:WPA:LD. -3.02e+15 2.10e+08 -14375383 <2e-16 ***
## WAR:BsR:LD. -3.39e+15 2.31e+08 -14655092 <2e-16 ***
## WPA:BsR:LD. -1.05e+16 3.47e+08 -30322039 <2e-16 ***
## WAR:Hard.:LD. 2.01e+14 1.03e+07 19588983 <2e-16 ***
## WPA:Hard.:LD. 9.60e+13 1.19e+07 8098416 <2e-16 ***
## BsR:Hard.:LD. -2.03e+14 1.19e+07 -17144881 <2e-16 ***
## WAR:WPA:BABIP -1.86e+17 1.46e+10 -12773856 <2e-16 ***
## WAR:BsR:BABIP -2.45e+17 1.61e+10 -15276073 <2e-16 ***
## WPA:BsR:BABIP -8.36e+17 2.60e+10 -32117043 <2e-16 ***
## WAR:Hard.:BABIP 2.07e+16 7.09e+08 29234547 <2e-16 ***
## WPA:Hard.:BABIP 1.70e+15 8.25e+08 2054410 <2e-16 ***
## BsR:Hard.:BABIP -1.38e+16 8.31e+08 -16627540 <2e-16 ***
## WAR:LD.:BABIP 3.90e+16 1.34e+09 29038226 <2e-16 ***
## WPA:LD.:BABIP -7.69e+15 1.58e+09 -4877793 <2e-16 ***
## BsR:LD.:BABIP -3.58e+16 1.62e+09 -22062173 <2e-16 ***
## Hard.:LD.:BABIP 2.21e+15 5.81e+07 38034569 <2e-16 ***
## WAR:WPA:BsR:Hard. 1.43e+15 6.42e+07 22246582 <2e-16 ***
## WAR:WPA:BsR:LD. 3.28e+15 1.28e+08 25671475 <2e-16 ***
## WAR:WPA:Hard.:LD. 3.99e+13 4.95e+06 8069608 <2e-16 ***
## WAR:BsR:Hard.:LD. 5.68e+13 5.60e+06 10140822 <2e-16 ***
## WPA:BsR:Hard.:LD. 2.32e+14 8.66e+06 26742237 <2e-16 ***
## WAR:WPA:BsR:BABIP 2.62e+17 9.31e+09 28119264 <2e-16 ***
## WAR:WPA:Hard.:BABIP 2.22e+15 3.40e+08 6534527 <2e-16 ***
## WAR:BsR:Hard.:BABIP 4.12e+15 3.87e+08 10633237 <2e-16 ***
## WPA:BsR:Hard.:BABIP 1.90e+16 6.48e+08 29264377 <2e-16 ***
## WAR:WPA:LD.:BABIP 9.02e+15 6.64e+08 13580108 <2e-16 ***
## WAR:BsR:LD.:BABIP 1.27e+16 7.78e+08 16385937 <2e-16 ***
## WPA:BsR:LD.:BABIP 3.74e+16 1.18e+09 31716530 <2e-16 ***
## WAR:Hard.:LD.:BABIP -8.60e+14 3.37e+07 -25526637 <2e-16 ***
## WPA:Hard.:LD.:BABIP -1.82e+14 3.85e+07 -4723835 <2e-16 ***
## BsR:Hard.:LD.:BABIP 6.89e+14 3.96e+07 17399051 <2e-16 ***
## WAR:WPA:BsR:Hard.:LD. -6.64e+13 3.00e+06 -22147768 <2e-16 ***
## WAR:WPA:BsR:Hard.:BABIP -5.45e+15 2.17e+08 -25107092 <2e-16 ***
## WAR:WPA:BsR:LD.:BABIP -1.24e+16 4.33e+08 -28669614 <2e-16 ***
## WAR:WPA:Hard.:LD.:BABIP -1.06e+14 1.56e+07 -6791720 <2e-16 ***
## WAR:BsR:Hard.:LD.:BABIP -2.23e+14 1.88e+07 -11873022 <2e-16 ***
## WPA:BsR:Hard.:LD.:BABIP -8.34e+14 2.94e+07 -28364684 <2e-16 ***
## WAR:WPA:BsR:Hard.:LD.:BABIP 2.57e+14 1.01e+07 25351353 <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: 197.2 on 146 degrees of freedom
## Residual deviance: 2811.4 on 83 degrees of freedom
## AIC: 2939
##
## Number of Fisher Scoring iterations: 25summary(forward_aic_add_squares)##
## Call:
## glm(formula = s.lifespan.type ~ poly(Hard., 2, raw = TRUE), family = binomial,
## data = bstats_trn)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -17.09636 8.18684 -2.09 0.037 *
## poly(Hard., 2, raw = TRUE)1 0.80120 0.40940 1.96 0.050 .
## poly(Hard., 2, raw = TRUE)2 -0.00945 0.00508 -1.86 0.063 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 197.20 on 146 degrees of freedom
## Residual deviance: 190.27 on 144 degrees of freedom
## AIC: 196.3
##
## Number of Fisher Scoring iterations: 4# forward_aic_add is the same as type vs Hard
model_bsr = glm(s.lifespan.type ~ BsR, data = bstats_trn, family = binomial)
model_war = glm(s.lifespan.type ~ WAR, data = bstats_trn, family = binomial)
model_wpa = glm(s.lifespan.type ~ WPA, data = bstats_trn, family = binomial)
model_babip = glm(s.lifespan.type ~ BABIP, data = bstats_trn, family = binomial)
model_hard_squared = glm(s.lifespan.type ~ I(Hard.^2), data = bstats_trn, family = binomial)
# check if linear Hard term is needed
anova( model_hard_squared, forward_aic_add_squares, test="LRT")## Analysis of Deviance Table
##
## Model 1: s.lifespan.type ~ I(Hard.^2)
## Model 2: s.lifespan.type ~ poly(Hard., 2, raw = TRUE)
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 145 195
## 2 144 190 1 4.89 0.027 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1so we do want to keep the linear Hard term
# check if squared Hard term is necessary
anova( forward_aic_add, forward_aic_add_squares, test="LRT")## Analysis of Deviance Table
##
## Model 1: s.lifespan.type ~ Hard.
## Model 2: s.lifespan.type ~ poly(Hard., 2, raw = TRUE)
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 145 194
## 2 144 190 1 4.28 0.039 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1So we do want to keep the squared Hard term
# check if the big model is any better than our squared model
anova( forward_aic_add_squares, backward_bic_add, test="LRT")## Analysis of Deviance Table
##
## Model 1: s.lifespan.type ~ poly(Hard., 2, raw = TRUE)
## Model 2: s.lifespan.type ~ WAR * WPA * BsR * Hard. * LD. * BABIP
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 144 190
## 2 83 2811 61 -2621Clearly the big model is super over-fitted. The squared Hard model appears to be best.
# Test classifiers on other half of 2019 data
hard_pred_2 = ifelse(predict(forward_aic_add, bstats_tst) > 0, "Single", "Multi")
bsr_pred_2 = ifelse(predict(model_bsr, bstats_tst) > 0, "Single", "Multi")
war_pred_2 = ifelse(predict(model_war, bstats_tst) > 0, "Single", "Multi")
wpa_pred_2 = ifelse(predict(model_wpa, bstats_tst) > 0, "Single", "Multi")
babip_pred_2 = ifelse(predict(model_babip, bstats_tst) > 0, "Single", "Multi")
back_bic_pred = ifelse(predict(backward_bic_add, bstats_tst) > 0, "Single", "Multi")
for_add_squares_pred = ifelse(predict(forward_aic_add_squares, bstats_tst) > 0, "Single", "Multi")
conf_mat_hard_2 = make_conf_mat(predicted=hard_pred_2, actual=bstats_tst$s.lifespan.type)
conf_mat_bsr_2 = make_conf_mat(predicted=bsr_pred_2, actual=bstats_tst$s.lifespan.type)
conf_mat_war_2 = make_conf_mat(predicted=war_pred_2, actual=bstats_tst$s.lifespan.type)
conf_mat_wpa_2 = make_conf_mat(predicted=wpa_pred_2, actual=bstats_tst$s.lifespan.type)
conf_mat_babip_2 = make_conf_mat(predicted=babip_pred_2, actual=bstats_tst$s.lifespan.type)
conf_mat_back_bic = make_conf_mat(predicted=back_bic_pred, actual=bstats_tst$s.lifespan.type)
conf_mat_for_add_squares = make_conf_mat(predicted=for_add_squares_pred, actual=bstats_tst$s.lifespan.type)
hard_mis_2 = mean(hard_pred_2 != bstats_tst$s.lifespan.type)
bsr_mis_2 = mean(bsr_pred_2 != bstats_tst$s.lifespan.type)
war_mis_2 = mean(war_pred_2 != bstats_tst$s.lifespan.type)
wpa_mis_2 = mean(wpa_pred_2 != bstats_tst$s.lifespan.type)
babip_mis_2 = mean(babip_pred_2 != bstats_tst$s.lifespan.type)
back_bic_mis = mean(back_bic_pred != bstats_tst$s.lifespan.type)
for_add_squares_mis = mean(for_add_squares_pred != bstats_tst$s.lifespan.type)
# Analyze Classifier (sensitivity and specificity)
sens_hard_2 = get_sens(conf_mat_hard_2)
sens_bsr_2 = get_sens(conf_mat_bsr_2)
sens_war_2 = get_sens(conf_mat_war_2)
sens_wpa_2 = get_sens(conf_mat_wpa_2)
sens_babip_2 = get_sens(conf_mat_babip_2)
sens_back_bic = get_sens(conf_mat_back_bic)
sens_for_add_squares = get_sens(conf_mat_for_add_squares)
spec_hard_2 = get_spec(conf_mat_hard_2)
spec_bsr_2 = get_spec(conf_mat_bsr_2)
spec_war_2 = get_spec(conf_mat_war_2)
spec_wpa_2 = get_spec(conf_mat_wpa_2)
spec_babip_2 = get_spec(conf_mat_babip_2)
spec_back_bic = get_spec(conf_mat_back_bic)
spec_for_add_squares = get_spec(conf_mat_for_add_squares)Results
BABIP Investigations
summary(full_batted_mod)##
## Call:
## lm(formula = BABIP ~ GB.FB + LD. + GB. + FB. + IFFB. + HR.FB +
## IFH. + Pull. + Cent. + Oppo. + Soft. + Med. + Hard. + Swing. +
## Contact., data = fangraphs_2018)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.05800 -0.01686 -0.00223 0.01545 0.06028
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.260538 7.384561 0.17 0.86
## GB.FB -0.017278 0.014633 -1.18 0.24
## LD. 0.024211 0.045143 0.54 0.59
## GB. 0.021244 0.045041 0.47 0.64
## FB. 0.018578 0.045098 0.41 0.68
## IFFB. 0.000211 0.000839 0.25 0.80
## HR.FB 0.000618 0.000543 1.14 0.26
## IFH. 0.004915 0.000817 6.02 1.8e-08 ***
## Pull. -0.064287 0.045399 -1.42 0.16
## Cent. -0.063640 0.045293 -1.41 0.16
## Oppo. -0.062810 0.045306 -1.39 0.17
## Soft. 0.031913 0.044206 0.72 0.47
## Med. 0.033493 0.044235 0.76 0.45
## Hard. 0.034019 0.044194 0.77 0.44
## Swing. -0.000499 0.000453 -1.10 0.27
## Contact. -0.000260 0.000490 -0.53 0.60
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.024 on 124 degrees of freedom
## Multiple R-squared: 0.543, Adjusted R-squared: 0.488
## F-statistic: 9.83 on 15 and 124 DF, p-value: 5.8e-15summary(babip_mod)##
## Call:
## lm(formula = BABIP ~ LD. + GB. + IFH. + Hard., data = fangraphs_2019)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.06628 -0.02235 -0.00048 0.01823 0.08969
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.075616 0.035890 2.11 0.037 *
## LD. 0.005834 0.000820 7.11 5.0e-11 ***
## GB. 0.001638 0.000401 4.08 7.4e-05 ***
## IFH. 0.003572 0.000878 4.07 7.8e-05 ***
## Hard. 0.000217 0.000442 0.49 0.625
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0289 on 143 degrees of freedom
## Multiple R-squared: 0.349, Adjusted R-squared: 0.331
## F-statistic: 19.2 on 4 and 143 DF, p-value: 1.24e-12#Breusch-Pagan test
bptest(babip_mod)##
## studentized Breusch-Pagan test
##
## data: babip_mod
## BP = 4.7, df = 4, p-value = 0.3#Shapiro test
shapiro.test(resid(babip_mod))##
## Shapiro-Wilk normality test
##
## data: resid(babip_mod)
## W = 0.98, p-value = 0.08#large leverage
sum(hatvalues(babip_mod) > 2 * mean(hatvalues(babip_mod)))## [1] 9#influential points
cd_babip_mod = cooks.distance(babip_mod)
large_cd_babip = cd_babip_mod > 4 / length(cd_babip_mod)
cd_babip_mod[large_cd_babip]## 5 34 56 97 145
## 0.03070 0.03343 0.06568 0.02867 0.03650par(mfrow = c(2, 2))
plot(babip_mod)
summary(wpa_mod)$r.squared## [1] 0.4154summary(war_mod)$r.squared## [1] 0.4707par(mfrow = c(2, 2))
plot(war_mod)
par(mfrow = c(2, 2))
plot(wpa_mod)
Salary Investigations
Predicting Pay Level Based on Performance Metrics
| . | Hard. | BsR | WAR | WPA | BABIP | Hard.BsR | log.Hard…BsR |
|---|---|---|---|---|---|---|---|
| Missclassification Rate | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 |
| Sensitivity | 0.7838 | 0.6081 | 0.3649 | 0.5541 | 0.5000 | 0.7432 | 0.7838 |
| Specifity | 0.3699 | 0.3836 | 0.6986 | 0.6164 | 0.3836 | 0.3836 | 0.3699 |
Predicting Salary Lifespan Based on Performance Metrics
| . | Hard. | BsR | WAR | WPA | BABIP | Overfitted.Model | Hard..Hard..2 |
|---|---|---|---|---|---|---|---|
| Missclassification Rate | 1.0000 | 1 | 1 | 1.0000 | 1 | 1.0000 | 1.0000 |
| Sensitivity | 0.0345 | 0 | 0 | 0.0517 | 0 | 0.8448 | 0.7838 |
| Specifity | 0.9663 | 1 | 1 | 0.9888 | 1 | 0.6629 | 1.0000 |
Discussion
BABIP Investigations
As discussed briefly in the abstract, BABIP is not a perfect
indicator of success at a plate appearance. The main difference between
batting average and BABIP is that batting average also
takes into account of strikeouts and base on balls. Strikeouts and base
on balls are important batter skills, but they are what is called ‘dead
balls’, which are plate appearance outcomes that end without putting the
ball in play. This means that approximately 10-40% of every player’s
plate appearances will not be affected by any raw batted ball skill.
Therefore, this study attempts to predict BABIP, which
attempts to only measure what percentage of balls that were put in play
were successful, from batted ball data predictors. Then we combine
BABIP, percentage of base on balls(BB),
percentage of strikeouts(K), baserunning skills(BsR and
Spd), to predict batter performance according to metrics known as
WAR(wins above replacement) and WPA(win probability added). Note that
this exercise still isn’t perfect, because there may be better
predictors that aren’t publically available for analysis and the minimum
number of plate appearances (PA)for each player is set to
150 in the data, which may not be enough.
Looking at the BABIP model that we created using the
2018 training data, the four predictors are line drive rate, groundball
rate, and infield hit rate, and hard hit rate. Immediately we can see
that these predictors are all important variables that professional
scouts and coaches often use to determine a good hitter. Line drives are
batted balls that fly at a medium angle which gives the opponent
fielders the least time to field. Also note that infield hit rates are a
measure of infield ground balls that end up being a hit, and both
IFH% and GB% are chosen predictors in the
model. This means that the percentage of fly balls(batted balls that fly
at a high angle) are relatively unimportant. I was also surprised to see
that contact% and swing% weren’t significant
predictors for BABIP. It is generally believed that players
who have high contact%(therefore are good at getting the
ball to the field instead of getting struck out) have good batting
skill, but in this excercise we see that putting more balls in play does
not necessarily mean that it is correlated to higher or lower
BABIP. In fact, in the full model we see that the
coefficient is negative and shows that hitters with high contact% are
often reaching out for worse balls and therefore running a lower
BABIP.
Running the Breusch-Pagan test, Shapiro test, and looking for high
leverage and influential points, we see that the model looks okay. The
QQ plot and other residual plots look great with the exception of a few
points. I’ve found that point 56 is Andrew Benintendi of the Boston Red
Sox, who has a .360 BABIP(compared to his .328 career
average) and has suddenly become a more flyball heavy hitter compared to
his past years.
Now we model BABIP along with BB,
K, BsR, and Spd to predict a
player’s WAR and WPA. WAR(wins
above replacement) is a simple numeric value calculated by giving a
value for each plate appearance, baserunning, and defensive outcome and
adding it all up over a season. WPA(Win probability added)
is similar, except that this time it assigns the added win probability
within a game after each plate appearance, baserunning, and defensive
outcome. Unfortunately our model is unable to add defensive statistics
and therefore incomplete, but it is still sufficient to provide a good
model for predicting WAR and WPA from each
predictor variable.
Salary Investigations
For the Salary Investigations we had attempted creating binary classifiers to help address these two questions:
- does overall performance necessarily translate to a higher pay rate, and if so what metrics might determine this?
- does overall performance necessarily translate to a longer lasting contract, and if so what metrics might determine this?
For this portion of the study we looked at only general metrics of
performance, such as BsR, WAR,
WPA, and BABIP. But we also included
Hard% as it appeared like it might be worth investigating.
And interestingly, out of the other performance metrics
Hard% appeared to show the most promise in playing a role
in a player’s salary.
For the first question we used data collected from the 2018 season as
the training data set and the 2019 season as the testing data set. We
added a new column to the dataset, pay.level, which was a
factor variable describing the pay level of the players. It’s just a
simple division using the median, those below the median yearly pay were
labeled as “low” and those receiving a yearly pay above the median were
labeled as “high”. We then used R’s pairs function to see
if there was any apparent correlation between the salary
parameter and any of the metrics mentioned above. Overall, there didn’t
really look like there was any correlation between anything with
salary. Though Hard% did appear to have some
correlation, maybe even an exponential relationship with
salary. Using this information we constructed a number of
training models using the 2018 season data. In particular we tried using
forward and backward stepwise AIC model building. However, attempting to
use the backward method resulted in errors so we were unable to use it.
The forward model building method yielded a model relying on
Hard% and BsR. Out of interest, a number of
models were constructed using the log(Hard%) as well to see
if there may indeed be an exponential relationship between
Hard% and salary. Plus simple models where
each metric of performance was used in making a classifier. We used
several anova tests to compare each of these classifier
models together. After all of these anova comparisons we came to
conclude that the pay.level vs Hard% +
BsR model was the best looking so far. We also tried a
pay.level vs log(Hard%) + BsR
model to account for the possibility of the exponential relationship in
Hard%. We then calculated the missclassification rate,
sensitivity, and specificity of each of the classifiers made and overall
the results were not that great. In the results section one will find a
summary of the 7 classifiers we created. The first five are where the
corresponding column title was the only predictor used to make the
classifier, then the last two columns are using the best model
(Hard% + BsR as predictors) and with the
log(Hard%) variant. All columns have similar values for
missclassification rate, sensitivity, and specificity. This may not be
so surprising since we divided the players between the median salary.
Although anything involving Hard% tended to have a high
Specificity, so how hard a player can hit a ball may play a factor in
how a player’s salary is determined. Though is definitely in question as
anything using Hard% as a predictor does not have a great
Sensitivity. Indeed, looking at the sensitivity of the other metrics
like WAR and WPA there are much better metrics
for determining how heavily paid a player is besides Hard%.
As for the log(Hard%), we again see not much of a
difference with just the regular Hard%, so there likely is
no exponential relationship. In hindsight it might have been good to
look at polynomial relationships, which is what was examined in the next
question studied.
For the second question we could only use the 2019 season data as
that was the only season for which Baseball-Reference provided contract
age and lifespan data. So we randomly divided the 2019 season data into
two equally sized groups for the training and testing data sets. Then we
labeled each player as having either a “Single” year contract or a
“Multi” year contract. These two labels were decided on when looking
over the collected salary data. It was observed that most contracted
players have just a single year on their contract. So we were interested
in seeing if player performance has any role in it. For this procedure
we followed a similar procedure to the first question, fitting various
models, using the forward and backward AIC model building techniques,
and using anova to compare models. This time we were able
to use a backward BIC model, but the resulting model was very large and
likely overfitted as demonstrated by the anova tests
conducted. Either way though we included this model in classifier
comparisons later on (the column labled “Overfitted model” is this large
model) Additionally, this time we looked at square terms. Interestingly,
the forward AIC model building method yielded a model using
Hard% and Hard%^2 as the predictors. Using a
\(\alpha = 0.05\), the
anova tests seemed to confirm that this was indeed the best
classifier model. However when we started calculating missclassification
rate, sensitivity and specificity, we saw some problems. These poor
calculations as seen in the results section are likely due to too small
of a sample size. There’s only about 150 players in our data set so this
likely may be why there’s these poor numbers. The best numbers are for
the Overfitted model generatd from the backwards BIC model generation
method. But that’s to be expected for an overfitted model. Although we
do see that anything using Hard% seems to have a great
sensitivity. Despite the questionable numbers, this does seem to suggest
that Hard% may play a role in how long their contract gets
written out for. Although these data sets are based on players having at
least 150 PA, so there’s more chances for them to get hard
hits. This would make sense as teams have an incentive to put longer
contracted players into play.
Overall, the results are inconclusive. We would need to use larger
data sets to make any conclusions involving a player’s salary. Although
Hard% did appear like it might play some role, so further
study on that metric with salary might be worth looking into.
Additionally, binary classifiers were probably not the best choice for
studying salary. It may be worth trying some typical linear regression
models in future studies.
Appendix
Quick Baseball Stat Definitions
| Term | Definition | Equation |
|---|---|---|
| H | Number of Hits Where Batter Reaches a Base Safely | |
| HR | Home Runs | |
| AB | Sum of Hits, Outs, and Times Reached By Error | |
| K | Strikeouts | |
| PA | Plate Appearances | |
| BB | Walks | |
| SF | Flyball out that results in a score | |
| BABIP | Batting Average on Balls In Play | \[ \frac{H - HR}{AB - K - HR + SF}\] |
| LD% | Line Drive Percentage | \[ \frac{LD}{H - HR}\] |
| GB% | Ground Ball Percentage | \[ \frac{GB}{H - HR}\] |
| FB% | Fly Ball Percentage | \[ \frac{FB}{H - HR}\] |
| IFFB% | Infield Fly Ball Percentage | \[ \frac{IFFB}{FB}\] |
| GB/FB | Ratio of GB to FB | |
| HR/FB | Ratio of Home Runs to Fly Balls | |
| IFH | Infield Hits | |
| IFH% | Percentage of ground balls that are infield hits | \[ \frac{IFH}{GB}\] |
| BUH | Hits that are bunts | |
| BUH% | Percentage of Bunts that result in a Hit | \[ \frac{BUH}{BU}\] |
| Pull% | Percentage of Hits to Batter’s Side of Field | |
| Cent% | Percentage of Hits to Center Field | |
| Oppo% | Percentage of Hits to Side of Field Opposite to Batter | |
| Soft% | Percentage of Hits that are Soft in Intensity | |
| Med% | Percentage of Hits that are Medium in Intensity | |
| Hard% | Percentage of Hits that are Hard in Intensity | |
| BB% | Percentage of PA that result in a Walk | \[\frac{BB}{PA}\] |
| K% | Percentage of PA that result in a Strikeout | \[\frac{K}{PA}\] |
| Spd | How fast the Batter can run | |
| Swing% | Percentage of Pithces the Batter Swung at | \[\frac{Swing}{PIT}\] |
| Contact% | Percentage of Pitches the Batter made contact with | \[\frac{Contact}{PIT}\] |
| BsR | Fangraph’s all encompassing base running statistic (see link for more details on calculation) | \[wSB + UBR + wGDP\] |
| WAR | Wins Above Replacement (Fangraph’s definition) | See Link |
| WPA | Win Probability Added | See Link |
Terms Specific to This study
| Term | Definition | Equation |
|---|---|---|
| Agreement | What kind of deal player has with organization | |
| CS | Player with a Contracted Salary | |
| PAE | Player who is Pre-Arbitration Eligible | |
| MLD | Player playing with a Minor League Deal | |
| NA | No payment information available for this player | |
| S1 | Player who is 1st-year Arbitration Eligible | |
| S2+ | Player who is 1st-year Arbitration Eligible (Super 2) | |
| S2 | Player who is 2nd-year Arbitration Eligible | |
| S3 | Player who is 3rd-year Arbitration Eligible | |
| Salary | How much the player makes in 1 year in USD; salaried players only (does not take into account other benefits, like vesting options, team options, etc) | \[\frac{Total Contract Amount}{S.lifespan}\] |
| S.age | How old the contract is in years | |
| S.lifespan | How long the contract lasts | |
| pay.level | Low or High, whether Salary is below or above the median | |
| s.lifespan.type | Single = 1 year contract or Multi = >1 year contract |
R Function Definitions
make_conf_mat = function(predicted, actual) {
table(predicted = predicted, actual = actual)
}
get_sens = function(conf_mat) {
conf_mat[2, 2] / sum(conf_mat[, 2])
}
get_spec = function(conf_mat) {
conf_mat[1, 1] / sum(conf_mat[, 1])
}
calc_loocv_rmse = function(model) {
sqrt(mean((resid(model) / (1 - hatvalues(model))) ^ 2))
}Contract Data Mining & Manipulation Python Scripts
Baseball-Reference Contract Data Scraper
A script that collects salary information for each player on each team. Note that the resulting Team csv files are not comma delimited, but are instead vertical bar delimited.
import urllib.request
from bs4 import BeautifulSoup
teams = ['NYY', 'TBR', 'BOS', 'TOR', 'BAL', 'MIN', 'CLE', 'CHW', 'KCR', 'DET', 'HOU', 'OAK', 'LAA', 'TEX', 'SEA',
'ATL', 'WSN', 'PHI', 'NYM', 'MIA', 'CHC', 'STL', 'MIL', 'CIN', 'PIT', 'LAD', 'SFG', 'ARI', 'SDP', 'COL']
for team in teams:
csv = open(team + ".csv", 'w')
url = urllib.request.urlopen('https://www.baseball-reference.com/teams/' + team + '/2019-roster.shtml')
html = url.read()
soup = BeautifulSoup(html, 'html.parser')
table = soup.find(id="all_the40man")
br_stem = 'https://www.baseball-reference.com'
for l in table.find_all('a'):
player_address = br_stem + l.get('href')
player_req = urllib.request.urlopen(player_address)
player_html = player_req.read()
player_soup = BeautifulSoup(player_html, 'html.parser')
info = player_soup.find(id='info')
player_name = info.find(itemprop='name').get_text()
for item in info.find_all('p'):
if item.get_text().find('2019 Contract Status:') != -1:
contract = item.get_text().split(":")[1].strip()
break
csv.write(player_name + "|" + contract + '\n')
csv.close()
print("completed: " + team)Qualifying Contract Data Script
Baseball-reference provides detailed descriptions of how much a player is making, such as what kind of agreement they currently have, how long it lasts, and how money much in total. For this study we were interested in yearly earnings, so this script was used to take a simple average of each contracted value over the contract’s lifespan.
from decimal import Decimal
teams = ['NYY', 'TBR', 'BOS', 'TOR', 'BAL', 'MIN', 'CLE', 'CHW', 'KCR', 'DET', 'HOU', 'OAK', 'LAA', 'TEX', 'SEA',
'ATL', 'WSN', 'PHI', 'NYM', 'MIA', 'CHC', 'STL', 'MIL', 'CIN', 'PIT', 'LAD', 'SFG', 'ARI', 'SDP', 'COL']
heading = "first,last,team,agreement,salary,s.age,s.lifespan\n"
all_csv = open("Salaries_Detailed.csv", "w")
all_csv.write(heading)
for team in teams:
csv = open(team + ".csv", 'r')
lines = csv.readlines()
csv.close()
for line in lines:
items = line.split("|")
name = items[0].split()
first = name[0]
last = name[1]
contract = items[1].strip()
all_csv.write(first + ",")
all_csv.write(last + ",")
all_csv.write(team + ",")
agreement = "ERR"
salary = "ERR"
salary_age = "ERR"
salary_lifespan = "ERR"
if (contract.find("Signed thru") != -1):
#salary calculations
agreement = "CS"
contract_items = contract.split(',')
ind = 0
while ind < len(contract_items) and contract_items[ind].find("$") == -1:
ind += 1
contract_dets = contract_items[ind].split()
years = contract_dets[0]
price_string = contract_dets[1].split('$')[1]
factor = 1
if (price_string[-1] == "k" or price_string[-1] == "K"):
factor = 1000
if (price_string[-1] == "m" or price_string[-1] == "M"):
factor = 1000000
if (price_string[-1] == "b" or price_string[-1] == "B"):
factor = 1000000000
price_dec = Decimal(price_string[:-1]) * factor
salary = str(price_dec/int(years))
year_span = contract_dets[2][1:-1].split("-")
curr_year = 19
salary_age = str(curr_year - int(year_span[0]))
salary_lifespan = years
if (contract == "Pre-Arb Eligible"):
#PAE
agreement = "PAE"
salary = "PAE"
salary_age = "PAE"
salary_lifespan = "PAE"
if (contract == "Minor League Deal"):
#MLD
agreement = "MLD"
salary ="MLD"
salary_age ="MLD"
salary_lifespan ="MLD"
if (contract == "Not Updated"):
#NA
agreement = "NA"
salary ="NA"
salary_age ="NA"
salary_lifespan ="NA"
if (contract == "2nd-Year Arb Eligible"):
#S2
agreement = "S2"
salary ="S2"
salary_age ="S2"
salary_lifespan ="S2"
if (contract == "1st-Year Arb Eligible"):
#S1
agreement = "S1"
salary ="S1"
salary_age ="S1"
salary_lifespan ="S1"
if (contract == "3rd-Year Arb Eligible"):
#S3
agreement = "S3"
salary ="S3"
salary_age ="S3"
salary_lifespan ="S3"
if (contract == "1st Year Arb Eligible (Super 2)"):
#S2+
agreement = "S2+"
salary ="S2+"
salary_age ="S2+"
salary_lifespan ="S2+"
all_csv.write(agreement + ",")
all_csv.write(salary + ",")
all_csv.write(salary_age + ",")
all_csv.write(salary_lifespan + "\n")
all_csv.close()Combining Contract Data with Fangraphs Data Script
# need unidecode to strip accented characters
import unidecode
fan = open("FanGraphs-Leaderboard.csv", "r")
sal = open("Salaries_Detailed.csv", "r")
fan_lines = fan.readlines()
sal_lines = sal.readlines()
fan.close()
sal.close()
comb = open("Combined-2019.csv", "w")
# column headings for combined dataset
comb_head = fan_lines[0].strip() + sal_lines[0][15:]
comb.write(comb_head)
line_count = 1
player_to_line_number = {}
for line in fan_lines[1:]:
# split line and grab first value (FirstName LastName), strip away the " characters
player = line.split(",")[0].replace('"', '')
# strip name suffixes
player = player.split()[0] + " " + player.split()[1]
# map player name to FG data line number
player_to_line_number[player] = line_count
line_count += 1
comb_lines = {}
for line in sal_lines[1:]:
# strip away accented characters and combine first last name separated by space.
player = unidecode.unidecode(line.split(",")[0] + " " + line.split(",")[1])
try:
# combine FG and BR data into single row, store in map by line number
comb_lines[player_to_line_number[player]] = fan_lines[player_to_line_number[player]].strip() \
+ "," + ",".join(line.split(",")[3:])
except:
pass
# now write the combined data rows into a file.
for c in range(1, len(fan_lines)):
try:
comb.write(comb_lines[c])
except:
# used to spot players with different names between FG and BR (EG, Peter vs Pete)
print()
comb.close()Citations and Readings
- fangraphs.com
- baseball-reference.com
- MLB Standard Stat Definitions
- Wikipedia’s page on WAR
- Sabermetrics
- Fangraph’s The beginner’s Guide to Using Statistics Properly
- What is WAR?
- WAR for Position Players
- WAR for Pitchers
- Calculating Position Player WAR, A Complete Example
- Win Values Explained: Part Six
- How to do baseball research: statistical databases and websites
- Wikipedia’s page on baseball statistics