Final Project Sports Analytics — 2025 Shortstop Market
Same structure as the 2023 notebook, re-run on the updated dataset
1. Load the Data
# CHANGED from 2023: new filename, and fileEncoding is required because the
# 2025 file contains the names Peña and Báez (Latin-1 encoded).
mlb_stats <- read.csv("SS_hittingstats_2025_cleaned.csv",
header = TRUE, fileEncoding = "latin1")
# Data structure
str(mlb_stats)## 'data.frame': 25 obs. of 17 variables:
## $ Player : chr "Francisco Lindor" "Willy Adames" "Zach Neto" "Trevor Story" ...
## $ Pos : chr "SS" "SS" "SS" "SS" ...
## $ Team : chr "NYM" "SF" "LAA" "BOS" ...
## $ GS : int 160 160 128 157 159 157 162 71 102 150 ...
## $ AB : int 644 591 502 612 590 623 629 255 380 589 ...
## $ H : int 172 133 129 161 144 184 166 61 103 152 ...
## $ X2B : int 35 22 29 29 24 47 31 9 19 23 ...
## $ X3B : int 0 2 1 0 3 6 7 1 0 2 ...
## $ HR : int 31 30 26 25 24 23 22 21 21 20 ...
## $ RBI : int 86 87 62 96 77 88 86 55 50 82 ...
## $ AVG : num 0.267 0.225 0.257 0.263 0.244 0.295 0.264 0.239 0.271 0.258 ...
## $ OBP : num 0.346 0.318 0.319 0.308 0.3 0.351 0.336 0.311 0.373 0.326 ...
## $ SLG : num 0.466 0.421 0.474 0.433 0.417 0.501 0.44 0.529 0.487 0.406 ...
## $ OPS : num 0.812 0.739 0.793 0.741 0.717 0.852 0.776 0.84 0.86 0.732 ...
## $ WAR : num 5.8 3.7 5.3 4.1 4.5 7.1 3.6 3.3 6.2 4.9 ...
## $ Salary2025: int 27050000 21000000 780000 22500000 27000000 7100000 770000 355483 32000000 22000000 ...
## $ Age : int 31 29 24 32 31 25 23 23 31 32 ...names(mlb_stats)## [1] "Player" "Pos" "Team" "GS" "AB"
## [6] "H" "X2B" "X3B" "HR" "RBI"
## [11] "AVG" "OBP" "SLG" "OPS" "WAR"
## [16] "Salary2025" "Age"Note that R automatically renames the doubles and triples columns to
X2B and X3B, because a column name cannot
begin with a number.
2. Data Cleaning
# CHANGED from 2023: the salary column already arrives NUMERIC in this file,
# so the gsub("[$, ]", "", ...) step used on Cash2023 is no longer needed.
# We verify instead of clean:
class(mlb_stats$Salary2025)## [1] "integer"# Standard data quality checks
colSums(is.na(mlb_stats))## Player Pos Team GS AB H X2B
## 0 0 0 0 0 0 0
## X3B HR RBI AVG OBP SLG OPS
## 0 0 0 0 0 0 0
## WAR Salary2025 Age
## 0 0 0sum(duplicated(mlb_stats))## [1] 0# Eliminating non-numeric columns
mlb_stats_num <- mlb_stats[sapply(mlb_stats, is.numeric)]
str(mlb_stats_num)## 'data.frame': 25 obs. of 14 variables:
## $ GS : int 160 160 128 157 159 157 162 71 102 150 ...
## $ AB : int 644 591 502 612 590 623 629 255 380 589 ...
## $ H : int 172 133 129 161 144 184 166 61 103 152 ...
## $ X2B : int 35 22 29 29 24 47 31 9 19 23 ...
## $ X3B : int 0 2 1 0 3 6 7 1 0 2 ...
## $ HR : int 31 30 26 25 24 23 22 21 21 20 ...
## $ RBI : int 86 87 62 96 77 88 86 55 50 82 ...
## $ AVG : num 0.267 0.225 0.257 0.263 0.244 0.295 0.264 0.239 0.271 0.258 ...
## $ OBP : num 0.346 0.318 0.319 0.308 0.3 0.351 0.336 0.311 0.373 0.326 ...
## $ SLG : num 0.466 0.421 0.474 0.433 0.417 0.501 0.44 0.529 0.487 0.406 ...
## $ OPS : num 0.812 0.739 0.793 0.741 0.717 0.852 0.776 0.84 0.86 0.732 ...
## $ WAR : num 5.8 3.7 5.3 4.1 4.5 7.1 3.6 3.3 6.2 4.9 ...
## $ Salary2025: int 27050000 21000000 780000 22500000 27000000 7100000 770000 355483 32000000 22000000 ...
## $ Age : int 31 29 24 32 31 25 23 23 31 32 ...3. Descriptive Statistics — Salary
options(scipen = 999)
average <- mean(mlb_stats_num$Salary2025, na.rm = TRUE)
paste("average:", format(round(average, 0), scientific = FALSE, big.mark = ","))## [1] "average: 11,370,970"variance <- var(mlb_stats_num$Salary2025, na.rm = TRUE)
paste("variance:", format(round(variance, 0), scientific = FALSE, big.mark = ","))## [1] "variance: 136,836,313,145,017"stdv <- sd(mlb_stats_num$Salary2025, na.rm = TRUE)
paste("standard deviation:", format(round(stdv, 0), scientific = FALSE, big.mark = ","))## [1] "standard deviation: 11,697,705"median1 <- median(mlb_stats_num$Salary2025, na.rm = TRUE)
paste("median:", format(round(median1, 0), scientific = FALSE, big.mark = ","))## [1] "median: 7,100,000"minimum <- min(mlb_stats_num$Salary2025, na.rm = TRUE)
paste("minimum:", format(round(minimum, 0), scientific = FALSE, big.mark = ","))## [1] "minimum: 355,483"maximum <- max(mlb_stats_num$Salary2025, na.rm = TRUE)
paste("maximum:", format(round(maximum, 0), scientific = FALSE, big.mark = ","))## [1] "maximum: 32,000,000"range1 <- range(mlb_stats_num$Salary2025, na.rm = TRUE)
paste("range:", format(round(range1, 0), scientific = FALSE, big.mark = ","))## [1] "range: 355,483" "range: 32,000,000"difference <- diff(range(mlb_stats_num$Salary2025, na.rm = TRUE))
paste("difference:", format(round(difference, 0), scientific = FALSE, big.mark = ","))## [1] "difference: 31,644,517"IQR1 <- IQR(mlb_stats_num$Salary2025, na.rm = TRUE)
paste("IQR:", format(round(IQR1, 0), scientific = FALSE, big.mark = ","))## [1] "IQR: 21,720,000"quant <- quantile(mlb_stats_num$Salary2025, na.rm = TRUE)
paste("quantile:", format(round(quant, 0), scientific = FALSE, big.mark = ","))## [1] "quantile: 355,483" "quantile: 780,000" "quantile: 7,100,000"
## [4] "quantile: 22,500,000" "quantile: 32,000,000"Interpretation. The mean (about $11.4M) sits well above the median ($7.1M), so salaries remain right-skewed: a handful of very large contracts pull the average upward. Compared with the 2023 dataset (mean $6.9M, median $2.0M), both measures are higher here — but note this file contains 25 regular starters rather than 47 players including part-time roles, so the 2025 sample is tilted toward established players.
4. Correlations
cor(mlb_stats_num)## GS AB H X2B X3B HR
## GS 1.000000000 0.93291948 0.7104934 0.5361154 0.3493666 0.23129082
## AB 0.932919479 1.00000000 0.8800305 0.6846461 0.3729589 0.36320721
## H 0.710493363 0.88003052 1.0000000 0.7871896 0.3748252 0.23249354
## X2B 0.536115408 0.68464614 0.7871896 1.0000000 0.4180387 0.19528543
## X3B 0.349366576 0.37295886 0.3748252 0.4180387 1.0000000 0.04226220
## HR 0.231290825 0.36320721 0.2324935 0.1952854 0.0422622 1.00000000
## RBI 0.576429730 0.69195837 0.6518262 0.4941619 0.2363763 0.60840723
## AVG 0.007388361 0.24415260 0.6714545 0.5193807 0.1466831 -0.09427099
## OBP 0.057793238 0.24104563 0.5725906 0.4061784 0.1062309 0.08002260
## SLG -0.369284372 -0.10130869 0.1972917 0.2974679 0.1647190 0.51883230
## OPS -0.224702310 0.04057509 0.3939761 0.3865093 0.1601505 0.38888098
## WAR 0.162504465 0.35058230 0.5443251 0.3992738 0.2498754 0.48322847
## Salary2025 0.088957287 0.14075398 0.1668609 -0.0240011 -0.1624172 0.25778798
## Age 0.157165583 0.15971402 0.1451401 -0.1020861 -0.2541877 0.05657190
## RBI AVG OBP SLG OPS WAR
## GS 0.57642973 0.007388361 0.05779324 -0.3692844 -0.22470231 0.1625045
## AB 0.69195837 0.244152595 0.24104563 -0.1013087 0.04057509 0.3505823
## H 0.65182624 0.671454497 0.57259063 0.1972917 0.39397614 0.5443251
## X2B 0.49416193 0.519380677 0.40617835 0.2974679 0.38650933 0.3992738
## X3B 0.23637628 0.146683067 0.10623091 0.1647190 0.16015054 0.2498754
## HR 0.60840723 -0.094270987 0.08002260 0.5188323 0.38888098 0.4832285
## RBI 1.00000000 0.226377258 0.19120746 0.2803547 0.27730016 0.3886038
## AVG 0.22637726 1.000000000 0.80936786 0.5339299 0.73020723 0.5610692
## OBP 0.19120746 0.809367862 1.00000000 0.5385084 0.81985092 0.7483527
## SLG 0.28035467 0.533929867 0.53850843 1.0000000 0.92396157 0.6412525
## OPS 0.27730016 0.730207234 0.81985092 0.9239616 1.00000000 0.7754389
## WAR 0.38860375 0.561069179 0.74835265 0.6412525 0.77543886 1.0000000
## Salary2025 0.18851753 0.144808469 0.14592666 0.1043562 0.13715170 0.2621518
## Age 0.04968468 0.074862051 0.06757985 -0.1492361 -0.07073278 0.1757184
## Salary2025 Age
## GS 0.08895729 0.15716558
## AB 0.14075398 0.15971402
## H 0.16686085 0.14514009
## X2B -0.02400110 -0.10208611
## X3B -0.16241717 -0.25418766
## HR 0.25778798 0.05657190
## RBI 0.18851753 0.04968468
## AVG 0.14480847 0.07486205
## OBP 0.14592666 0.06757985
## SLG 0.10435623 -0.14923610
## OPS 0.13715170 -0.07073278
## WAR 0.26215181 0.17571839
## Salary2025 1.00000000 0.91591552
## Age 0.91591552 1.00000000# Which variables move most closely with salary?
cash_cor <- cor(mlb_stats_num)[, "Salary2025"]
cash_cor[abs(cash_cor) >= 0.5]## Salary2025 Age
## 1.0000000 0.9159155Key finding — the biggest change from 2023. In the 2023 dataset this filter returned five performance statistics (AB, H, HR, RBI and WAR, with WAR strongest at 0.63). In the 2025 dataset, only Age survives the 0.5 cutoff, at roughly 0.92. Every performance statistic is weak: WAR about 0.26, HR 0.26, RBI 0.19, SLG only 0.10.
Interpreted directly: in this sample the market prices age and service time, not current on-field production.
Important caution about this screen. A correlation filter examines one variable at a time, so it can miss a variable that only becomes useful in combination with another. Section 8 shows exactly that happening with SLG.
5. Charts
options(scipen = 999)
boxplot(mlb_stats_num$Salary2025,
main = "Boxplot of Salaries", ylab = "Salary ($)")
options(scipen = 999)
hist(mlb_stats_num$Salary2025,
main = "Histogram of Player Salaries", xlab = "Salary ($)")
table(mlb_stats_num$Salary2025)##
## 355483 760000 764450 770000 770200 778500 780000 780500
## 1 1 1 1 1 1 1 1
## 780600 782300 879000 4100000 7100000 7550000 10000000 16500000
## 1 1 1 1 1 1 1 1
## 21000000 22000000 22500000 25000000 27000000 27000500 27050000 27272727
## 1 1 1 1 1 1 1 1
## 32000000
## 1# CHANGED from 2023: the 2023 notebook plotted RBI vs Salary. In the 2025
# data that relationship has essentially disappeared, so Age vs Salary is
# the informative scatterplot.
plot(x = mlb_stats_num$Age, y = mlb_stats_num$Salary2025,
main = "Scatterplot of Age vs. Salary",
xlab = "Age",
ylab = "Salary ($)")
abline(lm(Salary2025 ~ Age, data = mlb_stats_num), col = "blue", lwd = 2)
# kept for comparison with the 2023 notebook
plot(x = mlb_stats_num$RBI, y = mlb_stats_num$Salary2025,
main = "Scatterplot of RBI vs. Salary",
xlab = "RBI",
ylab = "Salary ($)")
# every numeric variable plotted against salary
par(mfrow = c(3, 5))
for (col in names(mlb_stats_num)) {
if (col != "Salary2025") {
plot(mlb_stats_num[[col]], mlb_stats_num$Salary2025,
main = paste(col, "vs. Salary"),
xlab = col,
ylab = "Salary ($)")
}
}
par(mfrow = c(1, 1))
6. Multicollinearity Check
pairs(mlb_stats_num[, c("WAR", "OPS", "HR", "RBI", "AVG", "OBP", "SLG", "Age")])
round(cor(mlb_stats_num[, c("WAR","OPS","HR","RBI","AVG","OBP","SLG","Age")]), 2)## WAR OPS HR RBI AVG OBP SLG Age
## WAR 1.00 0.78 0.48 0.39 0.56 0.75 0.64 0.18
## OPS 0.78 1.00 0.39 0.28 0.73 0.82 0.92 -0.07
## HR 0.48 0.39 1.00 0.61 -0.09 0.08 0.52 0.06
## RBI 0.39 0.28 0.61 1.00 0.23 0.19 0.28 0.05
## AVG 0.56 0.73 -0.09 0.23 1.00 0.81 0.53 0.07
## OBP 0.75 0.82 0.08 0.19 0.81 1.00 0.54 0.07
## SLG 0.64 0.92 0.52 0.28 0.53 0.54 1.00 -0.15
## Age 0.18 -0.07 0.06 0.05 0.07 0.07 -0.15 1.00OPS is by construction OBP + SLG, so those three should never all appear in the same model. In the 2023 notebook the final model contained both SLG and OPS, and the result was a SLG coefficient of about -132 million alongside an OPS coefficient of +68 million — signs that make no baseball sense. That is the classic symptom of multicollinearity, and we avoid it here.
7. Train / Test Split
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, ]
nrow(train); nrow(test)## [1] 20## [1] 5Limitation to state plainly. With 25 players, the training set holds roughly 20 rows and the test set only 5. Every test metric below is therefore noisy, and a single unusual player can swing it. The 2023 notebook had 37 training and 10 test rows.
8. Best Subset Selection
#install.packages("leaps", repos = "https://cran.r-project.org")
install.packages(c("leaps", "rpart", "randomForest"))
library(leaps)
# CHANGED from 2023: nvmax is capped at 5 rather than 10. With only ~20
# training rows, allowing 10 predictors would mean fitting 11 parameters to
# 20 observations, which guarantees overfitting.
best_subset <- regsubsets(Salary2025 ~ H + HR + RBI + WAR + AVG + OBP + SLG +
OPS + Age + GS + AB + X2B + X3B,
data = train, nvmax = 5)## Reordering variables and trying again:results <- summary(best_subset)
results## Subset selection object
## Call: regsubsets.formula(Salary2025 ~ H + HR + RBI + WAR + AVG + OBP +
## SLG + OPS + Age + GS + AB + X2B + X3B, data = train, nvmax = 5)
## 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
## Age FALSE FALSE
## GS FALSE FALSE
## AB FALSE FALSE
## X2B FALSE FALSE
## X3B FALSE FALSE
## OPS FALSE FALSE
## 1 subsets of each size up to 6
## 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 ) "*" " " "*" "*" " " " " "*" " " "*" " " " " "*" " "# Compare the two selection criteria
results$adjr2 # adjusted R-squared for each model size## [1] 0.8317741 0.9297060 0.9497740 0.9541701 0.9658326 0.9662991results$bic # BIC for each model size## [1] -30.73883 -46.33869 -51.27851 -51.40546 -55.66293 -54.42428which.max(results$adjr2) # size chosen by adjusted R-squared## [1] 6which.min(results$bic) # size chosen by BIC## [1] 5Why we select on BIC rather than adjusted R-squared. Adjusted R-squared applies only a mild penalty for extra predictors, and with a sample this small it keeps rising as variables are added — on the full 2025 data it does not peak until the largest model allowed. BIC penalizes complexity much more strongly and is the appropriate criterion for a small sample.
coef(best_subset, which.min(results$bic))## (Intercept) H RBI WAR SLG GS
## 49435008.2 32355.9 112473.3 2345496.7 -37178315.7 -304651.2Expected selection: Age and SLG. Across the full dataset, BIC selects the two-variable model containing Age and slugging percentage, and this choice is stable — repeating the selection across many different random training splits returns the same pair the large majority of the time.
This is precisely the case the correlation screen in Section 4 missed: SLG has a correlation with salary of only about 0.10 on its own, so a filter of 0.5 discards it. Yet once Age is in the model, SLG becomes clearly significant. A variable can be uninformative alone and informative in combination — which is the reason subset selection is worth running at all.
9. Linear Regression Model
lm_best <- lm(Salary2025 ~ Age + SLG, data = train)
summary(lm_best)##
## Call:
## lm(formula = Salary2025 ~ Age + SLG, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5597505 -2186914 492726 2108004 4589941
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -111763400 9107320 -12.272 0.0000000007138 ***
## Age 3237166 212870 15.207 0.0000000000249 ***
## SLG 81745745 16007918 5.107 0.0000877282834 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3206000 on 17 degrees of freedom
## Multiple R-squared: 0.9371, Adjusted R-squared: 0.9297
## F-statistic: 126.6 on 2 and 17 DF, p-value: 0.00000000006141# Predictions on the held-out test players
lm_pred <- predict(lm_best, newdata = test)
R_squared <- 1 - sum((test$Salary2025 - lm_pred)^2) /
sum((test$Salary2025 - mean(train$Salary2025))^2)
paste("Test R-squared:", round(R_squared, 4))## [1] "Test R-squared: 0.6302"prediction_examples <- data.frame(
Player = test$Player,
Actual = test$Salary2025,
Predicted = round(lm_pred, 0),
Difference = round(test$Salary2025 - lm_pred, 0)
)
prediction_examples## Player Actual Predicted Difference
## 2 Willy Adames 21000000 16529372 4470628
## 4 Trevor Story 22500000 27221819 -4721819
## 7 Elly De La Cruz 770000 -1340454 2110454
## 8 Colson Montgomery 355483 5934917 -5579434
## 16 Jeremy Peña 4100000 14632802 -10532802rmse <- function(actual, predicted) sqrt(mean((actual - predicted)^2))
mape <- function(actual, predicted) mean(abs((actual - predicted) / actual)) * 100
cat("Linear Regression RMSE: $", rmse(test$Salary2025, lm_pred), "\n")## Linear Regression RMSE: $ 6145016cat("Linear Regression MAPE:", mape(test$Salary2025, lm_pred), "%\n")## Linear Regression MAPE: 428.5586 %Watch for two familiar patterns in this output. First, RMSE and MAPE can disagree sharply: a miss of a few million dollars is minor for a $30M veteran and catastrophic for a $400K rookie, so MAPE is usually far worse than RMSE suggests. Second, a linear model can produce negative predicted salaries for the cheapest players, which is impossible in reality — the 2023 notebook produced two such cases. Both are consequences of forcing a straight line through a market with a hard floor at the league minimum.
10. Decision Tree
library(rpart)
dt_model <- rpart(Salary2025 ~ H + HR + RBI + WAR + AVG + OBP + SLG + OPS +
Age + GS + AB + X2B + X3B, data = train)
summary(dt_model)## Call:
## rpart(formula = Salary2025 ~ H + HR + RBI + WAR + AVG + OBP +
## SLG + OPS + Age + GS + AB + X2B + X3B, data = train)
## n= 20
##
## CP nsplit rel error xerror xstd
## 1 0.8703883 0 1.0000000 1.106419 0.1672385
## 2 0.0100000 1 0.1296117 1.106419 0.1672385
##
## Variable importance
## Age AB H HR RBI WAR
## 54 15 8 8 8 8
##
## Node number 1: 20 observations, complexity param=0.8703883
## mean=1.177744e+07, MSE=1.388803e+14
## left son=2 (13 obs) right son=3 (7 obs)
## Primary splits:
## Age < 30.5 to the left, improve=0.8703883, (0 missing)
## AB < 581 to the left, improve=0.2142085, (0 missing)
## AVG < 0.2575 to the left, improve=0.1570155, (0 missing)
## WAR < 4.15 to the left, improve=0.1508251, (0 missing)
## SLG < 0.3945 to the left, improve=0.1495154, (0 missing)
## Surrogate splits:
## AB < 424.5 to the right, agree=0.75, adj=0.286, (0 split)
## H < 110.5 to the right, agree=0.70, adj=0.143, (0 split)
## HR < 19.5 to the left, agree=0.70, adj=0.143, (0 split)
## RBI < 53.5 to the right, agree=0.70, adj=0.143, (0 split)
## WAR < 4.15 to the left, agree=0.70, adj=0.143, (0 split)
##
## Node number 2: 13 observations
## mean=3709658, MSE=2.356423e+13
##
## Node number 3: 7 observations
## mean=2.676046e+07, MSE=7.66789e+12printcp(dt_model)##
## Regression tree:
## rpart(formula = Salary2025 ~ H + HR + RBI + WAR + AVG + OBP +
## SLG + OPS + Age + GS + AB + X2B + X3B, data = train)
##
## Variables actually used in tree construction:
## [1] Age
##
## Root node error: 2777605125949242/20 = 138880256297462
##
## n= 20
##
## CP nsplit rel error xerror xstd
## 1 0.87039 0 1.00000 1.1064 0.16724
## 2 0.01000 1 0.12961 1.1064 0.16724plotcp(dt_model)
Check which variables the tree actually used for splitting — in the 2023 data it chose only AB (playing time) and Age, ignoring WAR, OBP and SLG entirely. Because AB is strongly correlated with H, HR and RBI, the tree simply picks whichever of that group splits cleanest first, after which the rest become redundant.
11. Random Forest
library(randomForest)
rf_model <- randomForest(Salary2025 ~ H + HR + RBI + WAR + AVG + OBP + SLG +
OPS + Age + GS + AB + X2B + X3B, data = train)
rf_model##
## Call:
## randomForest(formula = Salary2025 ~ H + HR + RBI + WAR + AVG + OBP + SLG + OPS + Age + GS + AB + X2B + X3B, data = train)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 4
##
## Mean of squared residuals: 79333758602349
## % Var explained: 42.88importance(rf_model)## IncNodePurity
## H 74765229114916
## HR 98713123196636
## RBI 157448501837517
## WAR 94423835753779
## AVG 104137390383972
## OBP 79611490977152
## SLG 116579442909353
## OPS 157065779223547
## Age 1089796844911076
## GS 138386978112315
## AB 208448324403096
## X2B 93298010350373
## X3B 7416252601903612. Model Comparison
dt_pred <- predict(dt_model, newdata = test)
rf_pred <- predict(rf_model, newdata = test)
cat("Linear Regression RMSE: $", rmse(test$Salary2025, lm_pred), "\n")## Linear Regression RMSE: $ 6145016cat("Decision Tree RMSE: $", rmse(test$Salary2025, dt_pred), "\n")## Decision Tree RMSE: $ 8211600cat("Random Forest RMSE: $", rmse(test$Salary2025, rf_pred), "\n\n")## Random Forest RMSE: $ 7782340cat("Linear Regression MAPE:", mape(test$Salary2025, lm_pred), "%\n")## Linear Regression MAPE: 428.5586 %cat("Decision Tree MAPE: ", mape(test$Salary2025, dt_pred), "%\n")## Decision Tree MAPE: 287.2238 %cat("Random Forest MAPE: ", mape(test$Salary2025, rf_pred), "%\n")## Random Forest MAPE: 887.4576 %The model with the lowest RMSE predicts salary most accurately in dollar terms. Report MAPE alongside it, because the two can point to different winners: linear regression often wins on RMSE while the tree-based models handle low-salary players far better and never predict a negative salary.
13. ANOVA — Is the Simpler Model Enough?
anova(lm_best)## Analysis of Variance Table
##
## Response: Salary2025
## Df Sum Sq Mean Sq F value Pr(>F)
## Age 1 2334932809198462 2334932809198462 227.216 0.00000000002864 ***
## SLG 1 267975882513040 267975882513040 26.077 0.00008772828338 ***
## Residuals 17 174696434237740 10276260837514
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Compare the selected model against Age alone
lm_simple <- lm(Salary2025 ~ Age, data = train)
anova(lm_simple, lm_best)## Analysis of Variance Table
##
## Model 1: Salary2025 ~ Age
## Model 2: Salary2025 ~ Age + SLG
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 18 442672316750780
## 2 17 174696434237740 1 267975882513040 26.077 0.00008773 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1If this test returns p < 0.05, SLG earns its place and the two-variable model is justified. If p >= 0.05, Age alone is sufficient and the simpler model should be preferred and reported.
# Save the trained model
saveRDS(lm_best, "linear_regression_model_2025.rds")Findings
The 2025 shortstop market prices age, not current production. Age correlates with salary at about 0.92 and explains roughly 84% of salary variation on its own. WAR — the single best summary of a player’s total value — correlates at only 0.26 and is not statistically significant once Age is in the model.
Best subset selection identifies Age and slugging percentage. Power is the one performance dimension the market appears to pay for. Total value, which includes defense and baserunning, is not priced.
This creates systematic bargains among young players. Several shortstops aged 23 to 25 produce at or above the level of veterans earning twenty to thirty times more.
Limitations. The sample is small (25 players, roughly 5 test rows) and covers a single season. Salary is normally set by a contract signed years earlier, so part of what Age captures is a career track record this dataset does not contain. Some of the young-player discount is also driven by league contract rules rather than by market error — meaning the bargain is a time-limited window rather than a permanent inefficiency.