Final Project Sports Analytics
Exploratory Data Analysis
setwd("E:/S9510/CAP4936")
mlb_stats <- read.csv("MLB Players-hittingstats-ss.csv", header = TRUE)
#Data structure
str(mlb_stats)## 'data.frame': 47 obs. of 17 variables:
## $ Player : chr "Trea Turner" "Bo Bichette" "Amed Rosario" "Xander Bogaerts" ...
## $ Pos : chr "SS" "SS" "SS" "SS" ...
## $ Team : chr "LAD" "TOR" "CLE" "BOS" ...
## $ GS : int 160 158 151 148 161 133 148 151 129 138 ...
## $ AB : int 652 652 637 557 630 522 591 593 481 563 ...
## $ H : int 194 189 180 171 170 152 150 145 135 134 ...
## $ X2B : int 39 43 26 38 25 24 31 24 22 31 ...
## $ X3B : int 4 1 9 0 5 1 6 1 5 0 ...
## $ HR : int 21 24 11 15 26 22 20 33 10 31 ...
## $ RBI : int 100 93 71 73 107 64 80 83 55 98 ...
## $ AVG : num 0.298 0.29 0.283 0.307 0.27 0.291 0.254 0.245 0.281 0.238 ...
## $ OBP : num 0.343 0.333 0.312 0.377 0.339 0.366 0.294 0.317 0.327 0.298 ...
## $ SLG : num 0.466 0.469 0.403 0.456 0.449 0.467 0.428 0.455 0.41 0.458 ...
## $ OPS : num 0.809 0.802 0.715 0.833 0.788 0.834 0.722 0.772 0.736 0.756 ...
## $ WAR : num 4.84 3.44 3.95 5.42 5.4 5.55 1.05 4.04 4.5 4.42 ...
## $ Cash2023: chr "$27,272,727 " "$6,100,000 " "$7,800,000 " "$30,000,000 " ...
## $ Age : int 29 24 26 29 28 27 22 28 25 26 ...names(mlb_stats)## [1] "Player" "Pos" "Team" "GS" "AB" "H"
## [7] "X2B" "X3B" "HR" "RBI" "AVG" "OBP"
## [13] "SLG" "OPS" "WAR" "Cash2023" "Age"#Turn the variable Cash2023 numeric
mlb_stats$Cash2023 <- as.numeric(gsub("[$, ]", "", mlb_stats$Cash2023))
str(mlb_stats)## 'data.frame': 47 obs. of 17 variables:
## $ Player : chr "Trea Turner" "Bo Bichette" "Amed Rosario" "Xander Bogaerts" ...
## $ Pos : chr "SS" "SS" "SS" "SS" ...
## $ Team : chr "LAD" "TOR" "CLE" "BOS" ...
## $ GS : int 160 158 151 148 161 133 148 151 129 138 ...
## $ AB : int 652 652 637 557 630 522 591 593 481 563 ...
## $ H : int 194 189 180 171 170 152 150 145 135 134 ...
## $ X2B : int 39 43 26 38 25 24 31 24 22 31 ...
## $ X3B : int 4 1 9 0 5 1 6 1 5 0 ...
## $ HR : int 21 24 11 15 26 22 20 33 10 31 ...
## $ RBI : int 100 93 71 73 107 64 80 83 55 98 ...
## $ AVG : num 0.298 0.29 0.283 0.307 0.27 0.291 0.254 0.245 0.281 0.238 ...
## $ OBP : num 0.343 0.333 0.312 0.377 0.339 0.366 0.294 0.317 0.327 0.298 ...
## $ SLG : num 0.466 0.469 0.403 0.456 0.449 0.467 0.428 0.455 0.41 0.458 ...
## $ OPS : num 0.809 0.802 0.715 0.833 0.788 0.834 0.722 0.772 0.736 0.756 ...
## $ WAR : num 4.84 3.44 3.95 5.42 5.4 5.55 1.05 4.04 4.5 4.42 ...
## $ Cash2023: num 27272727 6100000 7800000 30000000 27000000 ...
## $ Age : int 29 24 26 29 28 27 22 28 25 26 ...#Eliminating non-numeric columns
mlb_stats_num <- mlb_stats[sapply(mlb_stats, is.numeric)]
str(mlb_stats_num)## 'data.frame': 47 obs. of 14 variables:
## $ GS : int 160 158 151 148 161 133 148 151 129 138 ...
## $ AB : int 652 652 637 557 630 522 591 593 481 563 ...
## $ H : int 194 189 180 171 170 152 150 145 135 134 ...
## $ X2B : int 39 43 26 38 25 24 31 24 22 31 ...
## $ X3B : int 4 1 9 0 5 1 6 1 5 0 ...
## $ HR : int 21 24 11 15 26 22 20 33 10 31 ...
## $ RBI : int 100 93 71 73 107 64 80 83 55 98 ...
## $ AVG : num 0.298 0.29 0.283 0.307 0.27 0.291 0.254 0.245 0.281 0.238 ...
## $ OBP : num 0.343 0.333 0.312 0.377 0.339 0.366 0.294 0.317 0.327 0.298 ...
## $ SLG : num 0.466 0.469 0.403 0.456 0.449 0.467 0.428 0.455 0.41 0.458 ...
## $ OPS : num 0.809 0.802 0.715 0.833 0.788 0.834 0.722 0.772 0.736 0.756 ...
## $ WAR : num 4.84 3.44 3.95 5.42 5.4 5.55 1.05 4.04 4.5 4.42 ...
## $ Cash2023: num 27272727 6100000 7800000 30000000 27000000 ...
## $ Age : int 29 24 26 29 28 27 22 28 25 26 ...#mean
average <- mean(mlb_stats_num$Cash2023, na.rm = TRUE)
paste("average:", format(round(average, 0), scientific = FALSE, big.mark = ","))## [1] "average: 6,855,709"#variance
variance <- var(mlb_stats_num$Cash2023, na.rm = TRUE)
paste("variance:", format(round(variance, 0), scientific = FALSE, big.mark = ","))## [1] "variance: 93,571,304,657,024"#standard deviation
stdv <- sd(mlb_stats_num$Cash2023, na.rm = TRUE)
paste("standard deviation:", format(round(stdv, 0), scientific = FALSE, big.mark = ","))## [1] "standard deviation: 9,673,226"#median
median1 <- median(mlb_stats_num$Cash2023, na.rm = TRUE)
paste("median:", format(round(median1, 0), scientific = FALSE, big.mark = ","))## [1] "median: 2,000,000"#min number
minimum <- min(mlb_stats_num$Cash2023, na.rm = TRUE)
paste("minimum:", format(round(minimum, 0), scientific = FALSE, big.mark = ","))## [1] "minimum: 410,326"#max number
maximum <- max(mlb_stats_num$Cash2023, na.rm = TRUE)
paste("maximum:", format(round(maximum, 0), scientific = FALSE, big.mark = ","))## [1] "maximum: 36,000,000"#range of the numbers
range1 <- range(mlb_stats_num$Cash2023, na.rm = TRUE)
paste("range:", format(round(range1, 0), scientific = FALSE, big.mark = ","))## [1] "range: 410,326" "range: 36,000,000"#difference between max and min
difference <- diff(range(mlb_stats_num$Cash2023, na.rm = TRUE))
paste("difference:", format(round(difference, 0), scientific = FALSE, big.mark = ","))## [1] "difference: 35,589,674"#IQR
IQR1 <- IQR(mlb_stats_num$Cash2023, na.rm = TRUE)
paste("IQR:", format(round(IQR1, 0), scientific = FALSE, big.mark = ","))## [1] "IQR: 7,525,400"#quantile
quant <- quantile(mlb_stats_num$Cash2023, na.rm = TRUE)
paste("quantile:", format(round(quant, 0), scientific = FALSE, big.mark = ","))## [1] "quantile: 410,326" "quantile: 724,600" "quantile: 2,000,000"
## [4] "quantile: 8,250,000" "quantile: 36,000,000"names(mlb_stats_num)## [1] "GS" "AB" "H" "X2B" "X3B" "HR"
## [7] "RBI" "AVG" "OBP" "SLG" "OPS" "WAR"
## [13] "Cash2023" "Age"#correlation of the variables
cor(mlb_stats_num)## GS AB H X2B X3B HR
## GS 1.0000000 0.9878832 0.9387030 0.8804607 0.366752978 0.7339877
## AB 0.9878832 1.0000000 0.9741364 0.9156261 0.389106791 0.7840134
## H 0.9387030 0.9741364 1.0000000 0.9356879 0.375379829 0.7724508
## X2B 0.8804607 0.9156261 0.9356879 1.0000000 0.293089007 0.7399494
## X3B 0.3667530 0.3891068 0.3753798 0.2930890 1.000000000 0.1576920
## HR 0.7339877 0.7840134 0.7724508 0.7399494 0.157692045 1.0000000
## RBI 0.8833421 0.9243508 0.9298298 0.8896238 0.315943566 0.8973054
## AVG 0.2527911 0.3350714 0.4999607 0.4540903 0.115910827 0.2808246
## OBP 0.1430023 0.1935521 0.3283211 0.2978383 0.001986825 0.2055263
## SLG 0.1628489 0.2577958 0.3830350 0.3966754 0.102071214 0.5288822
## OPS 0.1672546 0.2520635 0.3908039 0.3880660 0.069776549 0.4423801
## WAR 0.7585507 0.7838553 0.8122353 0.7383154 0.308012465 0.7154940
## Cash2023 0.4710708 0.5099290 0.5628422 0.4634702 0.030049924 0.6053461
## Age 0.2663626 0.2643532 0.2483782 0.2137952 -0.098017716 0.1988217
## RBI AVG OBP SLG OPS WAR
## GS 0.8833421 0.2527911 0.143002251 0.1628489 0.16725465 0.7585507
## AB 0.9243508 0.3350714 0.193552094 0.2577958 0.25206353 0.7838553
## H 0.9298298 0.4999607 0.328321103 0.3830350 0.39080388 0.8122353
## X2B 0.8896238 0.4540903 0.297838268 0.3966754 0.38806597 0.7383154
## X3B 0.3159436 0.1159108 0.001986825 0.1020712 0.06977655 0.3080125
## HR 0.8973054 0.2808246 0.205526332 0.5288822 0.44238007 0.7154940
## RBI 1.0000000 0.3961871 0.285046646 0.4542913 0.42255309 0.7653890
## AVG 0.3961871 1.0000000 0.807340495 0.7975364 0.86254468 0.4335819
## OBP 0.2850466 0.8073405 1.000000000 0.7032172 0.87380311 0.3843565
## SLG 0.4542913 0.7975364 0.703217214 1.0000000 0.96018521 0.4288262
## OPS 0.4225531 0.8625447 0.873803113 0.9601852 1.00000000 0.4439624
## WAR 0.7653890 0.4335819 0.384356543 0.4288262 0.44396237 1.0000000
## Cash2023 0.5789837 0.3434524 0.373311756 0.3875104 0.41241558 0.6341681
## Age 0.2358260 0.1043517 0.054119374 0.1016015 0.09065536 0.1841738
## Cash2023 Age
## GS 0.47107079 0.26636259
## AB 0.50992895 0.26435320
## H 0.56284220 0.24837822
## X2B 0.46347021 0.21379517
## X3B 0.03004992 -0.09801772
## HR 0.60534606 0.19882166
## RBI 0.57898374 0.23582598
## AVG 0.34345236 0.10435169
## OBP 0.37331176 0.05411937
## SLG 0.38751040 0.10160148
## OPS 0.41241558 0.09065536
## WAR 0.63416813 0.18417379
## Cash2023 1.00000000 0.44225191
## Age 0.44225191 1.00000000If pairs are correlating with each other at 0.7+, that’s multicollinearity, and including both in the same regression could make the coefficients unstable or misleading. In that case, it’s often better to pick just one or two representative predictors rather than throwing all four in together — WAR is a strong single choice since it already tries to summarize overall value in one number.
#box plot chart
options(scipen = 999)
boxplot(mlb_stats_num$Cash2023, main="Boxplot of Salaries", ylab="Price ($)")
#histogram chart
options(scipen = 999)
hist(mlb_stats_num$Cash2023, main = "Histogram of Player Prices", xlab = "Price ($)")
#table
table(mlb_stats_num$Cash2023)##
## 410326 520429 536130 541940 632766 654193 661941 720000
## 1 1 1 1 1 1 1 1
## 720100 722000 723200 724200 725000 727600 730000 734500
## 1 1 1 1 1 1 1 1
## 738600 745750 754900 850000 950000 1800000 2000000 2525000
## 1 1 1 1 1 1 2 1
## 2662000 3000000 5000000 5585000 6000000 6100000 6500000 7000000
## 1 1 1 1 2 1 1 1
## 7800000 8700000 9000000 10000000 10250000 12500000 16000000 22000000
## 1 1 1 1 1 1 1 1
## 27000000 27272727 30000000 35000000 36000000
## 1 1 1 1 1# scatterplot
plot(x = mlb_stats_num$RBI, y = mlb_stats_num$Cash2023,
main = "Scatterplot of RBI vs. Salary",
xlab = "RBI",
ylab = "Price ($)")
#scatter plot against Cash2023
par(mfrow = c(3, 5))
for (col in names(mlb_stats_num)) {
if (col != "Cash2023") {
plot(mlb_stats_num[[col]], mlb_stats_num$Cash2023,
main = paste(col, "vs. Salary"),
xlab = col,
ylab = "Price ($)")
}
}
par(mfrow = c(1, 1)) 
Prediction Model
- Check for multicollinearity among predictors
pairs(mlb_stats_num[, c("WAR", "OPS", "HR", "RBI", "AVG", "OBP", "Age")])
If two predictors are highly correlated with each other (like AVG and
OBP are), keeping both adds little and can destabilize the model.
- Split into training and test sets
set.seed(123)
n <- nrow(mlb_stats_num)
train_idx <- sample(1:n, size = 0.8 * n)
train <- mlb_stats[train_idx, ]
test <- mlb_stats[-train_idx, ]- Choose the best model
Examine model results using summary function
Best Subset Selection
#install.packages("leaps", repos = "https://cran.r-project.org")
library(leaps)
#Best Subset Selection on the train dataset
best_subset <- regsubsets(Cash2023 ~ H + HR + RBI + WAR + AVG + OBP + SLG + OPS + Age + GS + AB + X2B + X3B,
data = train, nvmax = 10) # max number of predictors to consider
summary(best_subset)## Subset selection object
## Call: regsubsets.formula(Cash2023 ~ H + HR + RBI + WAR + AVG + OBP +
## SLG + OPS + Age + GS + AB + X2B + X3B, data = train, nvmax = 10)
## 13 Variables (and intercept)
## Forced in Forced out
## H FALSE FALSE
## HR FALSE FALSE
## RBI FALSE FALSE
## WAR FALSE FALSE
## AVG FALSE FALSE
## OBP FALSE FALSE
## SLG FALSE FALSE
## OPS FALSE FALSE
## Age FALSE FALSE
## GS FALSE FALSE
## AB FALSE FALSE
## X2B FALSE FALSE
## X3B FALSE FALSE
## 1 subsets of each size up to 10
## Selection Algorithm: exhaustive
## H HR RBI WAR AVG OBP SLG OPS Age GS AB X2B X3B
## 1 ( 1 ) " " " " " " "*" " " " " " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " "*" " " " " " " " " "*" " " " " " " " "
## 3 ( 1 ) " " "*" " " "*" " " " " " " " " "*" " " " " " " " "
## 4 ( 1 ) " " "*" " " "*" " " "*" " " " " "*" " " " " " " " "
## 5 ( 1 ) "*" "*" " " "*" " " " " " " " " "*" " " "*" " " " "
## 6 ( 1 ) "*" "*" " " " " " " " " "*" "*" "*" " " "*" " " " "
## 7 ( 1 ) "*" "*" " " "*" " " " " "*" "*" "*" " " "*" " " " "
## 8 ( 1 ) "*" "*" " " "*" " " "*" "*" "*" "*" " " "*" " " " "
## 9 ( 1 ) "*" "*" " " "*" " " "*" "*" "*" "*" " " "*" " " "*"
## 10 ( 1 ) "*" "*" " " "*" " " "*" "*" "*" "*" "*" "*" " " "*"#best values for the linear regression model
results <- summary(best_subset)
results$adjr2 # adjusted R-squared for each model size## [1] 0.3187883 0.4179625 0.4608279 0.4641609 0.4661644 0.4773790 0.4948131
## [8] 0.4884056 0.4858591 0.4688147which.max(results$adjr2) # which size is best## [1] 7coef(best_subset, which.max(results$adjr2)) # variables in that best model## (Intercept) H HR WAR SLG
## -17178396.63 215714.47 862406.66 1382414.00 -132522548.11
## OPS Age AB
## 67932216.23 941931.66 -76017.15#training the liner regression on the train dataset
lm_best <- lm(Cash2023 ~ H + HR + WAR + OBP + SLG + Age + AB, data = train)
summary(lm_best)##
## Call:
## lm(formula = Cash2023 ~ H + HR + WAR + OBP + SLG + Age + AB,
## data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11059212 -3469942 -107570 3723122 13742238
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -17008173 15143618 -1.123 0.27060
## H 216624 131236 1.651 0.10960
## HR 860552 289931 2.968 0.00595 **
## WAR 1381210 970707 1.423 0.16544
## OBP 67032528 46117600 1.454 0.15682
## SLG -64234967 33931047 -1.893 0.06836 .
## Age 940637 318519 2.953 0.00618 **
## AB -76193 39047 -1.951 0.06074 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6204000 on 29 degrees of freedom
## Multiple R-squared: 0.5923, Adjusted R-squared: 0.4939
## F-statistic: 6.018 on 7 and 29 DF, p-value: 0.0002162#Test R squared result on the test dataset
lm_pred <- predict(lm_best, newdata = test)
R_squared <- 1 - sum((test$Cash2023 - lm_pred)^2) / sum((test$Cash2023 - mean(train$Cash2023))^2)
paste("R-squared:", round(R_squared, 4))## [1] "R-squared: 0.7843"prediction_examples <- data.frame(Player = mlb_stats$Player[as.numeric(rownames(test))],
Actual = test$Cash2023,
Predicted = round(lm_pred, 0),
Difference = round(test$Cash2023 - lm_pred, 0)
)
prediction_examples# RMSE formula
rmse <- function(actual, predicted) sqrt(mean((actual - predicted)^2))
# Generate Predictions
lm_pred <- predict(lm_best, newdata = test)
cat("Linear Regression RMSE: $", rmse(test$Cash2023, lm_pred), "\n")## Linear Regression RMSE: $ 5848489Linear Regression: $5,848,489 average prediction error
Conclusion: The linear regression model’s predictions are, on average, about $5.8 million off from the actual salary.
mape <- function(actual, predicted) mean(abs((actual - predicted) / actual)) * 100
lm_mape <- mape(test$Cash2023, lm_pred)
cat("Linear Regression MAPE:", lm_mape,"%\n")## Linear Regression MAPE: 190.7629 %