Assignment1
Yujia Cao
2024-03-25
1.DATA EXPLORATION
The Moneyball training dataset encompasses 2,276 observations across 17 key variables, capturing an array of team performance metrics ranging from batting and pitching to fielding statistics. The dataset’s variables exhibit varying degrees of dispersion. The mean number of wins (TARGET_WINS) sits at 80.79, with a standard deviation of 15.75, reflecting a moderate variability in team performance from season to season. The median wins pegged at 82, suggesting that half of the teams won more than 82 games, while the other half won fewer.
Visual analysis through a bar chart of missing data highlights gaps in several variables, including TEAM_BATTING_SO, TEAM_BASERUN_SB, and particularly TEAM_BATTING_HBP, which has data for only a fraction of the teams. This suggests a necessity for data imputation to ensure a comprehensive analysis. Additionally, a correlation heatmap reveals significant relationships between various variables and the target variable of wins, with some metrics showing stronger correlations than others, indicating their potential predictive value for team success.
While the dataset provides a rich foundation for analyzing factors contributing to baseball team successes, addressing missing values and leveraging correlations will be pivotal in harnessing this data for predictive insights, aligning with the strategic Moneyball philosophy.
# Load necessary libraries
library(ggplot2)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(corrplot)## corrplot 0.92 loadedlibrary(mice)##
## Attaching package: 'mice'## The following object is masked from 'package:stats':
##
## filter## The following objects are masked from 'package:base':
##
## cbind, rbindlibrary(Amelia)## Loading required package: Rcpp## ##
## ## Amelia II: Multiple Imputation
## ## (Version 1.8.1, built: 2022-11-18)
## ## Copyright (C) 2005-2024 James Honaker, Gary King and Matthew Blackwell
## ## Refer to http://gking.harvard.edu/amelia/ for more information
## ### Load the CSV files
data1 <- read.csv("moneyball-training-data (1).csv")
data2 <- read.csv("moneyball-evaluation-data (1).csv")
training_data <- read.csv("moneyball-training-data (1).csv")
head(training_data)summary(training_data)## INDEX TARGET_WINS TEAM_BATTING_H TEAM_BATTING_2B
## Min. : 1.0 Min. : 0.00 Min. : 891 Min. : 69.0
## 1st Qu.: 630.8 1st Qu.: 71.00 1st Qu.:1383 1st Qu.:208.0
## Median :1270.5 Median : 82.00 Median :1454 Median :238.0
## Mean :1268.5 Mean : 80.79 Mean :1469 Mean :241.2
## 3rd Qu.:1915.5 3rd Qu.: 92.00 3rd Qu.:1537 3rd Qu.:273.0
## Max. :2535.0 Max. :146.00 Max. :2554 Max. :458.0
##
## TEAM_BATTING_3B TEAM_BATTING_HR TEAM_BATTING_BB TEAM_BATTING_SO
## Min. : 0.00 Min. : 0.00 Min. : 0.0 Min. : 0.0
## 1st Qu.: 34.00 1st Qu.: 42.00 1st Qu.:451.0 1st Qu.: 548.0
## Median : 47.00 Median :102.00 Median :512.0 Median : 750.0
## Mean : 55.25 Mean : 99.61 Mean :501.6 Mean : 735.6
## 3rd Qu.: 72.00 3rd Qu.:147.00 3rd Qu.:580.0 3rd Qu.: 930.0
## Max. :223.00 Max. :264.00 Max. :878.0 Max. :1399.0
## NA's :102
## TEAM_BASERUN_SB TEAM_BASERUN_CS TEAM_BATTING_HBP TEAM_PITCHING_H
## Min. : 0.0 Min. : 0.0 Min. :29.00 Min. : 1137
## 1st Qu.: 66.0 1st Qu.: 38.0 1st Qu.:50.50 1st Qu.: 1419
## Median :101.0 Median : 49.0 Median :58.00 Median : 1518
## Mean :124.8 Mean : 52.8 Mean :59.36 Mean : 1779
## 3rd Qu.:156.0 3rd Qu.: 62.0 3rd Qu.:67.00 3rd Qu.: 1682
## Max. :697.0 Max. :201.0 Max. :95.00 Max. :30132
## NA's :131 NA's :772 NA's :2085
## TEAM_PITCHING_HR TEAM_PITCHING_BB TEAM_PITCHING_SO TEAM_FIELDING_E
## Min. : 0.0 Min. : 0.0 Min. : 0.0 Min. : 65.0
## 1st Qu.: 50.0 1st Qu.: 476.0 1st Qu.: 615.0 1st Qu.: 127.0
## Median :107.0 Median : 536.5 Median : 813.5 Median : 159.0
## Mean :105.7 Mean : 553.0 Mean : 817.7 Mean : 246.5
## 3rd Qu.:150.0 3rd Qu.: 611.0 3rd Qu.: 968.0 3rd Qu.: 249.2
## Max. :343.0 Max. :3645.0 Max. :19278.0 Max. :1898.0
## NA's :102
## TEAM_FIELDING_DP
## Min. : 52.0
## 1st Qu.:131.0
## Median :149.0
## Mean :146.4
## 3rd Qu.:164.0
## Max. :228.0
## NA's :286# Check for missing values
missing_values <- sapply(training_data, function(x) sum(is.na(x)))
print(missing_values)## INDEX TARGET_WINS TEAM_BATTING_H TEAM_BATTING_2B
## 0 0 0 0
## TEAM_BATTING_3B TEAM_BATTING_HR TEAM_BATTING_BB TEAM_BATTING_SO
## 0 0 0 102
## TEAM_BASERUN_SB TEAM_BASERUN_CS TEAM_BATTING_HBP TEAM_PITCHING_H
## 131 772 2085 0
## TEAM_PITCHING_HR TEAM_PITCHING_BB TEAM_PITCHING_SO TEAM_FIELDING_E
## 0 0 102 0
## TEAM_FIELDING_DP
## 286training_data <- read.csv("moneyball-training-data (1).csv")
# Correlation matrix
correlation_matrix <- cor(training_data[,-1], use="complete.obs")
# Plotting correlation heatmap
corrplot(correlation_matrix, method="circle", type="upper", tl.cex=0.6)
# Boxplot for a key variable
ggplot(training_data, aes(y=TARGET_WINS)) +
geom_boxplot() +
theme_minimal() +
labs(title="Boxplot of Target Wins", y="Wins")
2.DATA PREPARATION
Step 1&2: I address missing values by first creating flags for any missing numeric values and then imputing these missing values with the median of their respective columns. This ensures no data point is lost while also keeping track of which values were originally missing.
Step3: To categorizes teams based on their home run performance into ‘Low’, ‘Medium’, or ‘High’ tiers, enabling easier analysis of performance segments.
Step 4: I applies a log transformation to the strikeouts by pitchers (TEAM_PITCHING_SO), which can help normalize data that may have a skewed distribution, making it more suitable for linear models.
Step 5: I introduces a new variable, ‘batting_efficiency’, which combines hits, walks, and strikeouts to provide a single metric indicative of a team’s offensive performance efficiency.
library(dplyr)
library(tidyr)
df <- read.csv("moneyball-training-data (1).csv")
# 1.Fix missing values by replacing them with the median of their respective columns
df <- df %>%
mutate(across(where(is.numeric), ~ifelse(is.na(.), median(., na.rm = TRUE), .)))
# 2.Create flags for missing values
df <- df %>%
mutate(across(where(is.numeric), ~ifelse(is.na(.), 1, 0), .names = "missing_{.col}"))
# 3.Transform data by bucketing the TEAM_BATTING_HR into Low, Medium, and High groups
df$HR_bucket <- cut(df$TEAM_BATTING_HR,
breaks = quantile(df$TEAM_BATTING_HR, probs = c(0, .33, .67, 1), na.rm = TRUE),
labels = c("Low", "Medium", "High"),
include.lowest = TRUE)
# 4.Mathematical transformation - log transformation of TEAM_PITCHING_SO to normalize the distribution
df$TEAM_PITCHING_SO_log <- log(df$TEAM_PITCHING_SO + 1)
# 5.Combine variables to create new ones - for example, a batting efficiency metric
df$batting_efficiency <- with(df, TEAM_BATTING_H / (TEAM_BATTING_BB + TEAM_BATTING_H + TEAM_BATTING_SO))
# View the transformed dataset
head(df)3.BUILD MODELS
Model 1: Using basic batting metrics
It’s might show that basic offensive metrics like hits and home runs have a positive relationship with wins, which aligns with expectations. The batting efficiency’s coefficient should also be positive, suggesting that efficient teams likely win more games.
# Load necessary library for linear modeling
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericmodel1 <- lm(TARGET_WINS ~ TEAM_BATTING_H + TEAM_BATTING_HR + batting_efficiency, data = df)
summary(model1)##
## Call:
## lm(formula = TARGET_WINS ~ TEAM_BATTING_H + TEAM_BATTING_HR +
## batting_efficiency, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -59.881 -9.067 0.386 9.452 46.025
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.714767 3.121605 6.956 4.55e-12 ***
## TEAM_BATTING_H 0.064581 0.002954 21.864 < 2e-16 ***
## TEAM_BATTING_HR -0.005754 0.007016 -0.820 0.412
## batting_efficiency -64.160085 6.257479 -10.253 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.93 on 2272 degrees of freedom
## Multiple R-squared: 0.2192, Adjusted R-squared: 0.2182
## F-statistic: 212.6 on 3 and 2272 DF, p-value: < 2.2e-16Model 2: Including pitching and fielding errors
it’s introduces a transformed pitching metric and fielding errors. A negative coefficient for TEAM_FIELDING_E would make sense, indicating that more errors are associated with fewer wins. The transformed pitching strikeouts variable should positively influence wins, reflecting the importance of strong pitching.
# Applying a square root transformation to reduce skewness and enhance model fit
df$TEAM_PITCHING_SO_sqrt = sqrt(df$TEAM_PITCHING_SO)
model2 <- lm(TARGET_WINS ~ TEAM_BATTING_HR + TEAM_PITCHING_SO_sqrt + TEAM_FIELDING_E, data = df)
summary(model2)##
## Call:
## lm(formula = TARGET_WINS ~ TEAM_BATTING_HR + TEAM_PITCHING_SO_sqrt +
## TEAM_FIELDING_E, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -75.473 -9.869 0.431 9.683 72.061
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 90.414477 1.760169 51.367 < 2e-16 ***
## TEAM_BATTING_HR 0.044331 0.006864 6.458 1.29e-10 ***
## TEAM_PITCHING_SO_sqrt -0.432181 0.058592 -7.376 2.27e-13 ***
## TEAM_FIELDING_E -0.007905 0.001737 -4.552 5.60e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.27 on 2272 degrees of freedom
## Multiple R-squared: 0.06164, Adjusted R-squared: 0.0604
## F-statistic: 49.75 on 3 and 2272 DF, p-value: < 2.2e-16Model 3: Interaction between home runs and strikeouts
It’s explores the interaction between batting and pitching. A positive interaction term would suggest that teams with both strong batting and pitching outperform others, a reasonable expectation.
model3 <- lm(TARGET_WINS ~ TEAM_BATTING_HR * TEAM_PITCHING_SO_log + TEAM_FIELDING_E, data = df)
summary(model3)##
## Call:
## lm(formula = TARGET_WINS ~ TEAM_BATTING_HR * TEAM_PITCHING_SO_log +
## TEAM_FIELDING_E, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -61.585 -9.496 0.464 9.422 62.389
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 73.104707 4.199675 17.407 < 2e-16 ***
## TEAM_BATTING_HR 0.963055 0.096645 9.965 < 2e-16 ***
## TEAM_PITCHING_SO_log 0.458925 0.577059 0.795 0.42653
## TEAM_FIELDING_E -0.006095 0.001851 -3.293 0.00101 **
## TEAM_BATTING_HR:TEAM_PITCHING_SO_log -0.133766 0.013855 -9.655 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.08 on 2271 degrees of freedom
## Multiple R-squared: 0.08506, Adjusted R-squared: 0.08345
## F-statistic: 52.78 on 4 and 2271 DF, p-value: < 2.2e-164.SELECT MODELS
library(ggplot2)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode# Function to calculate RMSE
calculate_rmse <- function(actual, predicted) {
sqrt(mean((predicted - actual) ^ 2))
}
# Evaluate each model
evaluate_model <- function(model, data) {
# Predictions
predictions <- predict(model, newdata = data)
# Calculating Metrics
rmse <- calculate_rmse(data$TARGET_WINS, predictions)
r_squared <- summary(model)$r.squared
adjusted_r_squared <- summary(model)$adj.r.squared
f_statistic <- summary(model)$fstatistic[1]
list(RMSE = rmse, R_Squared = r_squared, Adjusted_R_Squared = adjusted_r_squared, F_Statistic = f_statistic)
}
# Model evaluations
eval_model1 <- evaluate_model(model1, df)
eval_model2 <- evaluate_model(model2, df)
eval_model3 <- evaluate_model(model3, df)
# Print model evaluations
print("Model 1 Evaluation:")## [1] "Model 1 Evaluation:"print(eval_model1)## $RMSE
## [1] 13.91597
##
## $R_Squared
## [1] 0.2192038
##
## $Adjusted_R_Squared
## [1] 0.2181728
##
## $F_Statistic
## value
## 212.6167print("Model 2 Evaluation:")## [1] "Model 2 Evaluation:"print(eval_model2)## $RMSE
## [1] 15.25559
##
## $R_Squared
## [1] 0.06164112
##
## $Adjusted_R_Squared
## [1] 0.06040209
##
## $F_Statistic
## value
## 49.74949print("Model 3 Evaluation:")## [1] "Model 3 Evaluation:"print(eval_model3)## $RMSE
## [1] 15.06404
##
## $R_Squared
## [1] 0.08505702
##
## $Adjusted_R_Squared
## [1] 0.08344549
##
## $F_Statistic
## value
## 52.78047# Residual plots for the selected model (replace modelx with the chosen model after evaluation)
selected_model <- model1 # Assuming model1 is selected based on evaluation
par(mfrow=c(2,2))
# Plotting diagnostic plots
plot(selected_model)
Residuals vs. Fitted: This plot should ideally show no clear pattern, indicating that the model’s errors are randomly distributed and the linear relationship is properly captured. In this plot, there is a slight pattern, which may suggest that the model isn’t capturing all the predictive structure in the data.
Q-Q Plot (Quantile-Quantile): A perfect alignment along the line would suggest that the residuals are normally distributed. But Q-Q plot seems to indicate a normal distribution of residuals, although there may be slight deviations in the tails.
Scale-Location (Spread-Location): This plot checks for homoscedasticity, where the residuals are spread equally across all levels of the fitted values. The plot shows some signs of heteroscedasticity, as indicated by the non-random pattern where the spread increases for larger fitted values.
Residuals vs. Leverage: This helps identify influential points that could have a strong impact on the model’s predictions. The points outside the Cook’s distance lines are particularly influential. The plot indicates a few points with higher leverage, but they do not appear to be influential enough to be concerning, as they are within the Cook’s distance boundary.
Q-Q plots as best model & apply for each model’s residuals
The Q-Q plot helps identify if the residuals of my model are normally distributed by comparing them to a perfect normal distribution. A model whose residuals closely follow the line in the Q-Q plot indicates that the normality assumption is reasonably met.
# Model 1 Q-Q Plot
qqnorm(residuals(model1))
qqline(residuals(model1), col = "steelblue", lwd = 2)
title("Model 1: Q-Q Plot of Residuals")
# Model 2 Q-Q Plot
qqnorm(residuals(model2))
qqline(residuals(model2), col = "steelblue", lwd = 2)
title("Model 2: Q-Q Plot of Residuals")
# Model 3 Q-Q Plot
qqnorm(residuals(model3))
qqline(residuals(model3), col = "steelblue", lwd = 2)
title("Model 3: Q-Q Plot of Residuals")
Residuals vs Leverage plots as best model & apply for each model’s residuals
Selecting the best model based on the “Residuals vs Leverage” plot is insightful because it helps in identifying if any particular observations have a disproportionate influence on the model fit. Influential points can significantly sway the regression line, leading to a model that doesn’t generalize well to other data. Here’s the “Residuals vs Leverage” plot for each model and proceed with the model selection based on this criterion:
# Function to plot Residuals vs Leverage
plot_residuals_vs_leverage <- function(model) {
plot(influencePlot(model, id.method="identify", main="Residuals vs Leverage"),
sub="Circle size is proportional to Cook's Distance")
}
# Plotting Residuals vs Leverage for each model
par(mfrow=c(3, 1), mar=c(4, 4, 2, 1))
# Model 1
plot_residuals_vs_leverage(model1)## Warning in plot.window(...): "id.method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "id.method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "id.method" is not
## a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "id.method" is not
## a graphical parameter## Warning in box(...): "id.method" is not a graphical parameter## Warning in title(...): "id.method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "id.method" is not a
## graphical parameter
# Model 2
plot_residuals_vs_leverage(model2)## Warning in plot.window(...): "id.method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "id.method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "id.method" is not
## a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "id.method" is not
## a graphical parameter## Warning in box(...): "id.method" is not a graphical parameter## Warning in title(...): "id.method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "id.method" is not a
## graphical parameter
# Model 3
plot_residuals_vs_leverage(model3)## Warning in plot.window(...): "id.method" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "id.method" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "id.method" is not
## a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "id.method" is not
## a graphical parameter## Warning in box(...): "id.method" is not a graphical parameter## Warning in title(...): "id.method" is not a graphical parameter## Warning in plot.xy(xy.coords(x, y), type = type, ...): "id.method" is not a
## graphical parameter
Best Model select & Predict
# Model3 was determined to be the best model
best_model <- model3eval_df <- read.csv("moneyball-training-data (1).csv")
predictions <- predict(best_model)
# Add predictions to eval_df
eval_df$predicted_wins <- predictions
# View the results
head(eval_df)