Predicting Wins for Baseball Game Using Multiple Linear Regression

Loading Data:

The data set provided has been uploaded to github repository from where it has been loaded into the markdown using the code chunk below. The reason why it has uploaded to the github is to keep the work reproducible.

training <- read.csv("https://raw.githubusercontent.com/Umerfarooq122/Data_sets/main/moneyball-training-data.csv")

Now our data set has been loaded into the environment. Let’s display the first few rows of our data set.

knitr::kable(head(training))

Everything looks fine but on a first glance we can see that there is an INDEX column which will not be used in the analysis and will be removed from the data set later on in the data preparation stage.

Dimension of Data Set:

Let’s check out the dimension of our raw data set:

dim(training)

## [1] 2276   17

So we can see that we have got 2276 observations and 17 columns which contains INDEX, our target variable TARGET_WINS and all the features. Now these dimension might not remain the same depending on the steps we carry out during data preparation

Descriptive Summary Statistics:

Before we go ahead with plot the relational plot between target and features. let’s check out the descriptive summary for each column. Descriptive summary will give us the mean, median, min, max and quartiles for each column. Descriptive statistics help us to simplify large amounts of data in a sensible way. For instance, finding the average of a column by going through all the 2276 observations could be very hectic thus we use descriptive summary statistics. The code chunk below uses summary() function from base R to give all the summary statistics of each column

knitr::kable(summary(training))

As we can see that the table above gave us the summary for even INDEX column too. So before moving and waiting for data preparation remove that from our data set. The code chunk below remove the redundant INDEX column from our data set.

training <- training[,-1]

Plotting the Target Column:

Before moving ahead with relational plot and correlation between target and features. Let’s check out the distribution of of target column TARGET_WINS.

ggplot()+
  geom_histogram(data = training, mapping = aes(x=TARGET_WINS), bins = 50, color = "black", fill = "grey")+geom_vline(xintercept=mean(training$TARGET_WINS), color='red')+labs(x="Target Wins", y="Count",title ="Distribution of Target Wins")+theme_bw()

As we can see that the distribution of TARGET_WINS is fairly normal with mean around 80.8 (represented by red vertical line), min at 0 and maximum at 146. We can observe min = 0 is a worry-some, since it is very rare for a team to go win less so there might be something that needs to be corrected in the data. We can also revisit the summary statistics for TARGET_WINS column to confirm the figures above.

summary(training$TARGET_WINS)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   71.00   82.00   80.79   92.00  146.00

Relational Plots Between Target and Feature column:

Now we can go ahead and plot the relational plot and we can also check the correlation between target and features. One thing must be kept in mind that the correlation might change also carry out data preparation since we still have the uncleaned, dirty data set. As we know that there are at least 15 features in the data set and it will be impractical to plot the relational graph between the target column and each feature column so we will use pairs.panel from psych library to plot multiple feature on the same graphic. In the graph below our focus will be on the first row and first column. The first row tell us the correlation between TARGET_WINS and feature column and similarly the first column show the relational plot between the same columns. On the diagonal we have the distribution of each column.

Scatter_Matrix <- pairs.panels(training[, c(1, 2:6)], main = "Scatter Plot Matrix for Training Dataset")

In the matrix plot above the value in the first row after TARGET_WINS distribution plot shows the correlation between TARGET_WINS and TEAM_BATTING_H which is .39. Similarly the scatter plot below the TARGET_WINS distribution is the plot between TARGET_WINS and TEAM_BATTING_H with fitted regression line in red color. In the matrix above we have only considered the frist 5 features because to avoid over crowdness. For the remaining features we have the matrix plot below:

pairs.panels(training[, c(1, 7:11)], main = "Scatter Plot Matrix for Training Dataset")

pairs.panels(training[, c(1, 11:16)], main = "Scatter Plot Matrix for Training Dataset")

From the above matrix plots we can see that diagonals have distributions of each features and we can see that some of these distributions contains outliers that needs to be address in our data preparation stage.

Missing values:

Let’s check out the missing values in the columns. We will check the number of missing values in each column using the code check below:

knitr::kable(colSums(is.na(training)))

we can see that some the columns has missing values which also needs to be addressed in our data preparation stage.

Fixing Missing Values and Outliers:

training <- training[, !names(training) %in% c('TEAM_BATTING_HBP','TEAM_BASERUN_CS','TEAM_FIELDING_DP')]

dim(training)

## [1] 2276   13

knitr::kable(colSums(is.na(training)))

training$TEAM_BATTING_SO[is.na(training$TEAM_BATTING_SO)] <- median(training$TEAM_BATTING_SO, na.rm = TRUE)
training$TEAM_BASERUN_SB[is.na(training$TEAM_BASERUN_SB)] <- median(training$TEAM_BASERUN_SB, na.rm = TRUE)
training$TEAM_PITCHING_SO[is.na(training$TEAM_PITCHING_SO)] <- median(training$TEAM_PITCHING_SO, na.rm = TRUE)

dim(training)

## [1] 2276   13

training$TEAM_PITCHING_H[training$TEAM_PITCHING_H > 3*sd(training$TEAM_PITCHING_H)] <- median(training$TEAM_PITCHING_H)
training$TEAM_PITCHING_BB[training$TEAM_PITCHING_BB > 3*sd(training$TEAM_PITCHING_BB)] <- median(training$TEAM_PITCHING_BB)
training$TEAM_PITCHING_SO[training$TEAM_PITCHING_SO > 3*sd(training$TEAM_PITCHING_SO)] <- median(training$TEAM_PITCHING_SO)
training$TEAM_FIELDING_E[training$TEAM_FIELDING_E > 3*sd(training$TEAM_FIELDING_E)] <- median(training$TEAM_FIELDING_E)

ggplot(melt(training), aes(x=value)) + geom_histogram(color = 'black', fill = 'grey') + facet_wrap(~variable, scale='free') + labs(x='', y='Frequency')+theme_bw()

Model 1:

Even though we already know that some the features like TEAM_PITCHING_BB, TEAM_PITCHING_HR has weak correlation with our target which TARGET_WINS but we in our first model we will try to fit all the features in the model and then check there statistical significance. The code chunk below will create our first model

m1 <- lm(TARGET_WINS ~., training)

We can check the intercept and coefficients of the model above using summary function. The summary function also give us the p-value of each coefficient which shows the statistical significance of the that coefficient. So let’s apply summary function to our model:

summary(m1)

## 
## Call:
## lm(formula = TARGET_WINS ~ ., data = training)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -57.862  -8.582   0.353   8.578  56.341 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -1.414986   5.607021  -0.252  0.80079    
## TEAM_BATTING_H    0.038461   0.003769  10.204  < 2e-16 ***
## TEAM_BATTING_2B  -0.019783   0.009431  -2.098  0.03604 *  
## TEAM_BATTING_3B   0.102623   0.017633   5.820 6.72e-09 ***
## TEAM_BATTING_HR   0.075863   0.027554   2.753  0.00595 ** 
## TEAM_BATTING_BB   0.035305   0.004431   7.968 2.52e-15 ***
## TEAM_BATTING_SO   0.011331   0.004418   2.565  0.01039 *  
## TEAM_BASERUN_SB   0.030714   0.004590   6.692 2.76e-11 ***
## TEAM_PITCHING_H   0.006600   0.001173   5.625 2.08e-08 ***
## TEAM_PITCHING_HR -0.036080   0.024315  -1.484  0.13798    
## TEAM_PITCHING_BB -0.018563   0.007410  -2.505  0.01230 *  
## TEAM_PITCHING_SO -0.007138   0.003582  -1.993  0.04643 *  
## TEAM_FIELDING_E  -0.022174   0.003612  -6.138 9.82e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.49 on 2263 degrees of freedom
## Multiple R-squared:   0.27,  Adjusted R-squared:  0.2661 
## F-statistic: 69.74 on 12 and 2263 DF,  p-value: < 2.2e-16

As we cans see that TEAM_PITCHING_HR has a very high p-value which means that it is not statistically significant and it will be in our best interest to remove that from our model. Below are the plots which shows residuals vs fitted values, QQ plot, residuals distributions, standardized residuals and residuals vs leverage.

par(mfrow=c(2,2))
plot(m1)

hist(resid(m1), main="Histogram of Residuals")

Model 2:

As we saw in model 1 that TEAM_PITCHING_HR has a high p-value. Let’s create another model without TEAM_PITCHING_HRand see how that performs.

m2<- lm(TARGET_WINS~TEAM_BATTING_H+TEAM_BATTING_2B+TEAM_BATTING_3B+TEAM_BATTING_HR+TEAM_BATTING_BB+TEAM_BATTING_SO+TEAM_BASERUN_SB+TEAM_PITCHING_H+TEAM_PITCHING_BB+TEAM_FIELDING_E, training)

Now our model is ready let’s check the summary of our model.

summary(m2)

## 
## Call:
## lm(formula = TARGET_WINS ~ TEAM_BATTING_H + TEAM_BATTING_2B + 
##     TEAM_BATTING_3B + TEAM_BATTING_HR + TEAM_BATTING_BB + TEAM_BATTING_SO + 
##     TEAM_BASERUN_SB + TEAM_PITCHING_H + TEAM_PITCHING_BB + TEAM_FIELDING_E, 
##     data = training)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -66.792  -8.559   0.415   8.600  58.424 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       2.296543   5.513356   0.417 0.677053    
## TEAM_BATTING_H    0.037750   0.003753  10.059  < 2e-16 ***
## TEAM_BATTING_2B  -0.017073   0.009419  -1.813 0.070013 .  
## TEAM_BATTING_3B   0.098129   0.017593   5.578 2.73e-08 ***
## TEAM_BATTING_HR   0.039122   0.009739   4.017 6.08e-05 ***
## TEAM_BATTING_BB   0.040121   0.004191   9.573  < 2e-16 ***
## TEAM_BATTING_SO   0.003934   0.002289   1.719 0.085833 .  
## TEAM_BASERUN_SB   0.030545   0.004591   6.653 3.58e-11 ***
## TEAM_PITCHING_H   0.004893   0.001021   4.790 1.78e-06 ***
## TEAM_PITCHING_BB -0.025057   0.007169  -3.495 0.000483 ***
## TEAM_FIELDING_E  -0.021481   0.003610  -5.950 3.10e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.52 on 2265 degrees of freedom
## Multiple R-squared:  0.2661, Adjusted R-squared:  0.2628 
## F-statistic: 82.11 on 10 and 2265 DF,  p-value: < 2.2e-16

As we can that after removing the TEAM_PITCHING_HR from the model our r-squared has dropped down a bit but at the same time features like TEAM_BATTING_SO and TEAM_BATTING_2B has lost their significance too, with p-values way over .05. So let’s remove these from our analysis and create a new model.

par(mfrow=c(2,2))
plot(m2)

hist(resid(m2), main="Histogram of Residuals")

Model 3:

In our model 3 we have removed all the three weakly correlated and low statistically significant features from our consideration. The following code chunk contains our model.

m3<- lm(TARGET_WINS~TEAM_BATTING_H+TEAM_BATTING_3B+TEAM_BATTING_HR+TEAM_BATTING_BB+TEAM_BASERUN_SB+TEAM_PITCHING_H+TEAM_PITCHING_BB+TEAM_FIELDING_E, training)

Let’s check the summary of our model 3.

summary(m3)

## 
## Call:
## lm(formula = TARGET_WINS ~ TEAM_BATTING_H + TEAM_BATTING_3B + 
##     TEAM_BATTING_HR + TEAM_BATTING_BB + TEAM_BASERUN_SB + TEAM_PITCHING_H + 
##     TEAM_PITCHING_BB + TEAM_FIELDING_E, data = training)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -68.259  -8.493   0.383   8.484  56.343 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      10.000771   3.926302   2.547 0.010927 *  
## TEAM_BATTING_H    0.031665   0.002537  12.482  < 2e-16 ***
## TEAM_BATTING_3B   0.098564   0.017475   5.640 1.91e-08 ***
## TEAM_BATTING_HR   0.047017   0.007669   6.131 1.03e-09 ***
## TEAM_BATTING_BB   0.038505   0.004128   9.328  < 2e-16 ***
## TEAM_BASERUN_SB   0.033250   0.004347   7.650 2.96e-14 ***
## TEAM_PITCHING_H   0.004579   0.001012   4.526 6.33e-06 ***
## TEAM_PITCHING_BB -0.024725   0.007151  -3.458 0.000555 ***
## TEAM_FIELDING_E  -0.021393   0.003558  -6.013 2.12e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.53 on 2267 degrees of freedom
## Multiple R-squared:  0.2645, Adjusted R-squared:  0.2619 
## F-statistic: 101.9 on 8 and 2267 DF,  p-value: < 2.2e-16

As we can see that all the features used in the model 3 are statistically significant. Features like TEAM_BATTING_H, TEAM_BATTING_3B, TEAM_BATTING_HR, TEAM_BATTING_BB, TEAM_BASERUN_SB and TEAM_PITCHING_H which are Base hits by batters, Triples by batters, homeruns by batters, walks by batters, stolen bases, and hits allowed, respectively, are all contributing positively towards the wins. Similarly, TEAM_PITCHING_BB and TEAM_FIELDING_E which are walks allowed and errors, respectively, are contributing negatively towards the wins.

Even though r-square went a touch high but all the features are statistically significant. Below are the plots which shows residuals vs fitted values, QQ plot, residuals distributions, standardized residuals and residuals vs leverage.

par(mfrow=c(2,2))
plot(m3)

hist(resid(m3), main = "Histogram of Residuals")

Selecting Models:

We are going with model 3 for prediction since all the feature all statistically significant and there is no weak correlation between the target and features. Another reason why we are not selecting m1 which has all the features and a better r-squared value is the result from analysis of variance (ANOVA). The code chunk below show ANOVA between the biggest model and the samllest model in terms of features used.

anova(m1, m3)

## Analysis of Variance Table
## 
## Model 1: TARGET_WINS ~ TEAM_BATTING_H + TEAM_BATTING_2B + TEAM_BATTING_3B + 
##     TEAM_BATTING_HR + TEAM_BATTING_BB + TEAM_BATTING_SO + TEAM_BASERUN_SB + 
##     TEAM_PITCHING_H + TEAM_PITCHING_HR + TEAM_PITCHING_BB + TEAM_PITCHING_SO + 
##     TEAM_FIELDING_E
## Model 2: TARGET_WINS ~ TEAM_BATTING_H + TEAM_BATTING_3B + TEAM_BATTING_HR + 
##     TEAM_BATTING_BB + TEAM_BASERUN_SB + TEAM_PITCHING_H + TEAM_PITCHING_BB + 
##     TEAM_FIELDING_E
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)   
## 1   2263 412102                               
## 2   2267 415188 -4   -3085.9 4.2365 0.00203 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

As we can see that the p-value is .00203 which below .05, meaning that the null hypothesis can not be rejected. Which in turn means that the results produce by each model i.e. m1 and m3 are not that varied from each other. So we will go with the smaller model in order to save computing power. Below are some of the main features of model 3 (m3)

sum_m3 <- summary(m3)
RSS <- c(crossprod(m3$residuals))
MSE <- RSS/length(m3$residuals)
print(paste0("Mean Squared Error: ", MSE))

## [1] "Mean Squared Error: 182.419913764793"

print(paste0("Root MSE: ", sqrt(MSE)))

## [1] "Root MSE: 13.5062916363002"

print(paste0("Adjusted R-squared: ", sum_m3$adj.r.squared))

## [1] "Adjusted R-squared: 0.261903459682437"

print(paste0("F-statistic: ",sum_m3$fstatistic[1]))

## [1] "F-statistic: 101.90657831176"

Here is the distribution of the residuals which show no pattern accross the horizontal line.

plot(resid(m3))
abline(h=0, col=2)

Predicting the Data:

Now we have selected our model so we can go ahead and load the testing data set to predict the TARGET_WINS using m3 model. The following code chunk load the data in our environment.

testing <- read.csv("https://raw.githubusercontent.com/Umerfarooq122/Data_sets/main/moneyball-evaluation-data.csv")

Here is the first few row of our testing data set:

knitr::kable(head(testing))

testing <- testing[, !names(testing) %in% c('INDEX','TEAM_BATTING_HBP','TEAM_BASERUN_CS','TEAM_FIELDING_DP')]

knitr::kable(head(testing))

dim(testing)

## [1] 259  12

knitr::kable(colSums(is.na(testing)))

testing <- mice(testing, m=5, maxit = 5, method = 'pmm')
testing <- complete(testing)

testing$TEAM_PITCHING_H[testing$TEAM_PITCHING_H > 3*sd(testing$TEAM_PITCHING_H)] <- median(testing$TEAM_PITCHING_H)
testing$TEAM_PITCHING_BB[testing$TEAM_PITCHING_BB > 3*sd(testing$TEAM_PITCHING_BB)] <- median(testing$TEAM_PITCHING_BB)
testing$TEAM_PITCHING_SO[testing$TEAM_PITCHING_SO > 3*sd(testing$TEAM_PITCHING_SO)] <- median(testing$TEAM_PITCHING_SO)
testing$TEAM_FIELDING_E[testing$TEAM_FIELDING_E > 3*sd(testing$TEAM_FIELDING_E)] <- median(testing$TEAM_FIELDING_E)

final <- predict(m3, newdata = testing, interval="prediction")
knitr::kable(head(final,3))