LA Dodgers Analysis
2024-03-20
Firstly we import all of the necessary packages for the analysis.
Then the data sets are imported.
lad_batting <- read.csv('LAD_batting.csv')
lad_pitching <- read.csv('LAD_pitching.csv')For batting analysis, to remove bias from small sample sizes, we filter by over 100 plate appearances for batting, and over 100 batters faced for pitchers.
lad_batting <- lad_batting %>%
filter(Plate_Appearances > 100)
lad_pitching <- lad_pitching %>%
filter(Batters_Faced > 100)Next we formulate a new column called Position Group in order to break down the numerous positions into some smaller categories.
position_group <- function(position){
if (position %in% c('C')){
return('Catcher')}
else if (position %in% c('1B', '2B', '3B','UT', 'SS','MI', 'CF','IF','CI')){
return('Infielder')}
else if(position %in% c('RF', 'LF','CF','OF')){
return('Outfielder')}
else if(position %in% c('DH')){
return('Designated Hitter')}
else if(position %in% c('P')){
return('Pitcher')}
else{
return('Unknown')}
}
lad_batting <- lad_batting %>%
mutate(PositionGroup = sapply(Position, position_group))The mean of the batting average and rbi per team is calculated and saved as a new data frame.
team_batting_average_per_year <- lad_batting %>%
group_by(Year) %>%
summarise(LeagueAverage = mean(Batting_Average))
team_rbi_per_year <- lad_batting%>%
group_by(Year) %>%
summarise(LeagueRBI = mean(Runs_Batted_In))Q1. Carefully constructed numerical and graphical summaries (using ggplot) of 5 relevant variables.
The first variable is a histogram of batting average, commonly used to evaluate batting performance. It is calculated as total hits/total at bats. A higher batting average indicates a better batting performance.
summary(lad_batting$Batting_Average)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0640 0.2340 0.2600 0.2557 0.2820 0.3960The team median batting average is 0.2600 which aligns with data from the league.
ggplot(lad_batting, aes(x=Batting_Average)) + geom_histogram(binwidth = 0.01, fill=("blue")) + labs(x="Player Batting Average",y="Frequency", title="Histogram of Player Batting Average")
Batting average within the data set is roughly normally distributed. The highest batting average with over 100 plate appearances was left fielder Manny Ramirez in the 2008 season, who had 74 hits from his 187 at bats.
Slugging percentage is a measure of how many ‘bases’ a player gains per at bat. It is calculated as (1(singles) + 2(doubles) + 3(triples) + 4(home runs))/at bats. It is a measure of how hard a player hits the ball, and how often.
summary(lad_batting$Slugging_Percentage)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0910 0.3330 0.3830 0.3886 0.4500 0.7430Team median slugging percentage is 0.3830.
ggplot(lad_batting, aes(x=Slugging_Percentage))+geom_boxplot(fill="red") + labs(x="Player Slugging Percentage")
Whilst the inter-quartile range is quite small, the range of all data points is quite large. Some players have very high or very low slugging percentages.
The next visualisation is a breakdown of players with over 100 plate appearances by position.
table(lad_batting$Position)##
## 1B 2B 3B C CF CI DH IF LF MI OF P RF SS UT
## 86 91 100 129 93 11 3 26 95 11 67 27 86 92 52This instead breaks it down by position group to reduce the number of levels.
table(lad_batting$PositionGroup)##
## Catcher Designated Hitter Infielder Outfielder
## 129 3 562 248
## Pitcher
## 27There are more infielders than any other position group who reached 100 plate appearances. There are very few players who are listed as a designated hitter, as normally the designated hitter in each game is listed as another position in the data set. Also, whilst there are more pitchers in the data set than any other position by far, when specifying at least 100 plate appearances most of these pitchers are removed.
Earned run average is a statistic used to grade pitching performance, and is a measure of how many earned runs (runs scored without the aid of a fielding error) they allow per 9 innings. A smaller ERA value indicates a better pitching performance.
summary(lad_pitching$Earned_Run_Average)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.650 2.760 3.485 3.651 4.240 13.140ggplot(lad_pitching, aes(x=Earned_Run_Average)) + geom_histogram(binwidth=0.2, fill="red") + labs(x="Earned Run Average", y= "Count")
The data for earned run average is approximately normal. There is one extreme outlier, Orel Hershiser in the 2000 season who recorded a woeful 13.14 ERA. He retired after the season. This was a bad end to an otherwise stellar 17 year MLB career, with an average career ERA of 3.48.
The final visualisation is fielding independent pitching. It only measures outcomes that come directly from the interactions between the hitter and the pitcher, eliminating any effects from the quality of the fielding team.
summary(lad_pitching$Fielding_Independent_Pitching)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.860 3.110 3.675 3.749 4.250 9.330ggplot(lad_pitching, aes(x=Fielding_Independent_Pitching)) + geom_boxplot(fill="blue") + labs(x="Fielding Independent Pitching")
The best recorded Fielding Independent Pitching figure was recorded by Éric Gagné in the 2003 season, a season in which he was awarded the Cy Young award for the best pitcher in baseball, and set a record by converting 84 consecutive save opportunities.
Q2. a. A scatterplot with two continuous variables, where points are coloured based on a categorical variable. Whilst both statistics quantify positive batting actions, on base percentage and slugging percentage represent slightly different things. Slugging percentage is a metric that is based on a players ability to power hit, whereas on base percentage includes walks and hit by pitches.
ggplot(lad_batting, aes(x=On_Base_Percentage,y=Slugging_Percentage, color=PositionGroup)) + geom_point() + labs(x="On Base Percentage",y="Slugging Percentage", color="Position Group", title="On Base Percentage against Slugging Percentage, Coloured based on position") 
There is a clear correlation between these two values, indicating that players who are better at power hitting are also better at getting on base. You can also see grouping of the pitchers towards the bottom left corner, indicating they are generally worse hitters. There is no real distinction between catchers, infielders, designated hitters and outfielders within the data set.
Q2b. A bar chart with error bars. A pitcher earns a win for a few reasons: - They pitch at least 5 innings and leave when their team is in the lead, and the team goes on to win. - They enter the game with the team behind or drawing, and the team goes on to win the game.
The bars represent the mean winning percentage of the Dodgers pitchers for each year. The error bars represent the maximum and minimum recorded winning percentage by a single pitcher for each year.
winning_percentage <- lad_pitching %>%
group_by(Year) %>%
filter(Year > 2010) %>%
summarise(League_Winning_Percentage = mean(Winning_Percentage),
Min_Winning_Percentage = min(Winning_Percentage),
Max_Winning_Percentage = max(Winning_Percentage))
ggplot(winning_percentage, aes(x=Year,y=League_Winning_Percentage)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin=Min_Winning_Percentage,ymax=Max_Winning_Percentage), width=0.2) + labs(title="Mean Winning Percentage with Error Bars (2010-2023)", x ="Year", y="Winning Percentage")
A lot of the years have very large ranges, indicating both pitchers with very high and very low winning percentage are common within each pitching group.
Q2c. A plot showing overlapping histograms for two different groups. Ensure the transparency setting allows for visibility of both distributions.
We plot the number of innings pitched between dominant hands
ggplot(lad_pitching, aes(x = Innings_Pitched, fill = Dominant_Hand)) +
geom_histogram(alpha = 0.5, position = "identity", bins = 30) +
labs(title = "Overlapping Histograms of Innings Pitched by Dominant Hand",x = "Innings Pitched",y = "Frequency",fill = "Pitching Hand") + scale_fill_manual(values = c("blue", "red")) +theme_minimal()
This graphic shows that the distributions for left and right handed pitchers are remarkably similar. This suggests that the number of innings pitched is entirely independent from dominant hand. It can also be seen that less innings are pitched across the board by left handed pitchers, showing that left handed pitchers are less common. However, the data suggests that the percentage of MLB pitchers who are left handed is much higher than the global average of left-handedness which is around 10%.
Q3. Fit, evaluate, interpret, and compare two logistic models. Both models should focus on the same response variable and include at least two predictors.
For this section, I have chosen to use the same response variable and predictors, and vary the type of linear model fitted. A home run is the best outcome from an individual at bat, where the ball is hit into the playing area and all four bases are rounded in one go, scoring a run. Often this comes as the result of the ball being hit very hard and a long way into the stands, and home runs are some of the most impressive events in baseball. We would like to use some batting statistics as a predictor of number of home runs hit.
ggplot(lad_batting, aes(x=Batting_Average,y=Home_Runs)) + geom_point() + labs(x="Batting Average", y="Home Runs")
ggplot(lad_batting, aes(x=Slugging_Percentage,y=Home_Runs)) + geom_point() + labs(x="Slugging Percentage",y="Home Runs")
Here we can see that both batting average and slugging percentage appear to be positively correlated with the number of home runs hit. Therefore, we use these two statistics as the predictors in our model. The first model uses a Gaussian family, modelling the response variable (home runs) as continuous.
linmod1 <- glm(Home_Runs ~ Batting_Average + Slugging_Percentage, family=gaussian, data=lad_batting)
summary(linmod1)##
## Call:
## glm(formula = Home_Runs ~ Batting_Average + Slugging_Percentage,
## family = gaussian, data = lad_batting)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -20.9226 -3.0557 -0.2545 2.7022 21.0366
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.258 1.138 -9.016 <2e-16 ***
## Batting_Average -86.050 6.477 -13.285 <2e-16 ***
## Slugging_Percentage 106.652 2.945 36.216 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 30.13261)
##
## Null deviance: 80377 on 968 degrees of freedom
## Residual deviance: 29108 on 966 degrees of freedom
## AIC: 6054.9
##
## Number of Fisher Scoring iterations: 2The summary of this gaussian model indicates that batting average and slugging percentage are very good indicators of number of home runs, with p scores of <2e-16 suggesting the chance of observing this result by random chance is almost zero. To better analyse the fit of the data we plot the residuals.
plot(lad_batting$Batting_Average,linmod1$residuals, main="Residuals vs Batting Average", xlab="Batting Average",ylab="Residuals") +
abline(h = 0, col = "blue")
## integer(0)This residual plot shows that batting average fits the model well. The residual points are pretty uniformly distributed, with clustering centrally near 0 residual.
plot(lad_batting$Slugging_Percentage,linmod1$residuals, main="Residuals vs Slugging Percentage",xlab="Slugging Percentage",ylab="Residuals") +
abline(h = 0, col = "blue")
## integer(0)This plot shows a worse fit of slugging percentage to the data. The plot exhibits heteroscedastic behaviour, suggesting that as the slugging percentage value grows, the fit of the model is worse.
qqnorm(linmod1$residuals)
The qq plot for this model is pretty strong. It is approximately the desired straight diagonal line.
The second linear model takes the same response and predictor variables, but instead fits a Poisson model. This is appropriate for home runs as they are discrete, only taking integer values and unable to be measured to an increased resolution. Also, the probability of hitting a home run in a single at bat is very low, however over a sufficient time period such as an entire season, that several home runs will be observed. This behavior fits a poisson distribution.
linmod2 <- glm(Home_Runs ~ Batting_Average + Slugging_Percentage, family=poisson ,data=lad_batting)
summary(linmod2)##
## Call:
## glm(formula = Home_Runs ~ Batting_Average + Slugging_Percentage,
## family = poisson, data = lad_batting)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -9.8800 -1.4893 -0.4646 0.9634 4.1913
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.17544 0.07472 -2.348 0.0189 *
## Batting_Average -7.40226 0.39703 -18.644 <2e-16 ***
## Slugging_Percentage 10.31139 0.16570 62.227 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 8198.0 on 968 degrees of freedom
## Residual deviance: 2991.7 on 966 degrees of freedom
## AIC: 6292.2
##
## Number of Fisher Scoring iterations: 5Again, we observe that batting average and slugging percentage are very good indicators for the number of home runs hit.
plot(lad_batting$Batting_Average,linmod2$residuals, main="Residuals vs Batting Average", xlab="Batting Average", ylab="Residuals") + abline(h=0,col="blue")
## integer(0)This residual plot is relatively normal, with a pretty even distribution that is clustered at the center. However, there is a line of -1.0 residuals at the bottom of the plot. I would estimate that this means the model is poor at dealing with individuals with very low numbers of home runs.
plot(lad_batting$Slugging_Percentage,linmod2$residuals, main="Residuals vs Slugging Percentage",xlab="Slugging Percentage",ylab="Residuals") + abline(h=0,col="blue")
## integer(0)The general fit for the residuals is very good, the data is randomly spread with some clustering towards the centre of the graph. Again, we can observe a residual line caused by players with low numbers of home runs.
qqnorm(linmod2$residuals)
This inability to model players with very low numbers of home runs affects the qqplot, which exhibits the same straight line behaviour. The model would be a much better fit with a minimum number of home runs required.
Comparisons: Both models found batting average and slugging percentage to be highly significant indicators of the number of home runs a player hits. The poisson model is generally a much better fit to the data, with lower residual deviance, at 2991.7 against 29108. The poisson model also fits the slugging percentage better, as the distribution of the residual plot is much more uniform. However, the Gaussian model is better at handling the very low home run values, which the poisson model struggles with. I would advise that a poisson model is used to model home runs, but with a minimum home run requirement to reduce some of the residuals.
Q4. Fit, evaluate, interpret, and compare two linear mixed effect models. Both models should focus on the same response variable and include at least one predictor and one random effect.
We want model the number of stolen bases. A stolen base is when a player runs between bases whilst the pitching team is in possession of the ball, not following a hit by a teammate. The idea is to reach the base before the fielding player can throw the ball there, requiring the hitter to be fast and time their steal attempt correctly. The player’s name has a big impact on the number of stolen bases due to each players individual tenancies, but a lot of levels, so it seems appropriate to include it as a random effect.
lmem1 <- lmer(Stolen_Bases ~ (1|Name)+On_Base_Percentage, data=lad_batting)
lmem1## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: Stolen_Bases ~ (1 | Name) + On_Base_Percentage
## Data: lad_batting
## REML criterion at convergence: 6921.659
## Random effects:
## Groups Name Std.Dev.
## Name (Intercept) 7.078
## Residual 7.020
## Number of obs: 969, groups: Name, 335
## Fixed Effects:
## (Intercept) On_Base_Percentage
## -10.24 50.16mean(lad_batting$Stolen_Bases)## [1] 7.154799There are 335 individuals within the data set. The mean number of stolen bases is 7.154799. The standard deviation of the random effect Name is 7.078, so it seems very important.
names1 <- ranef(lmem1)$Name
hist(names1$`(Intercept)`, breaks=30,xlab="Intercept",main="Intercepts based on the random effect 'Name'")
Higher values for intercept suggests that the player is better at stealing bases.
Now we will have a look at the residuals.
plot(lad_batting$On_Base_Percentage,lmem1@resp$wtres,xlab="On Base Percentage",ylab="Residuals")
The residuals look pretty good, bar some extreme positive outliers. Before moving onto predictions, we will outline the other model, which uses a poisson family instead. This is appropriate for stolen bases as they are discrete, and cannot be measured to a higher resolution.
lmem2 <- glmer(Stolen_Bases ~ (1|Name)+On_Base_Percentage, data=lad_batting, family="poisson")
lmem2## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Stolen_Bases ~ (1 | Name) + On_Base_Percentage
## Data: lad_batting
## AIC BIC logLik deviance df.resid
## 5760.560 5775.189 -2877.280 5754.560 966
## Random effects:
## Groups Name Std.Dev.
## Name (Intercept) 1.283
## Number of obs: 969, groups: Name, 335
## Fixed Effects:
## (Intercept) On_Base_Percentage
## -1.783 8.231The standard deviation for Name in the poisson model is 1.283, so it suggests that name has less of an impact on the number of stolen bases than the previous gaussian model.
names2 <- ranef(lmem2)$Name
hist(names2$`(Intercept)`, breaks=30)
Here, the intercepts are pretty normally distributed. This provides further evidence that the names variable is less important in the poisson model.
plot(lad_batting$On_Base_Percentage,lmem2@resp$wtres,xlab="On Base Percentage",ylab="Residuals")
The poisson model seems like a better fit for the On Base Percentage, with smaller values for the residuals.
Now we move onto comparing predictions for the two models. Current Dodgers right fielder Mookie Betts is proficient at stealing bases, so we will compare the model predictions for his on base percentages for each season he has played, and the number of bases he stole. In 2023, Betts stole 14 bases with OBP of 0.408. In 2022 he stole 12 bases with OBP of 0.340. In 2021, he stole 10 with OBP 0.367. In 2020 he stole 10 with OBP 0.366.
mookie_betts <- subset(lad_batting, Name == "Mookie Betts")
predict(lmem1,newdata=data.frame(Name="Mookie Betts",On_Base_Percentage=0.408))## 1
## 12.76841predict(lmem2,newdata=data.frame(Name="Mookie Betts",On_Base_Percentage=0.408))## 1
## 2.717263In 2023, model 1 predicts 12.76841 stolen bases which is pretty accurate, and model 2 is a long way off.
predict(lmem1,newdata=data.frame(Name="Mookie Betts",On_Base_Percentage=0.340))## 1
## 9.357347predict(lmem2,newdata=data.frame(Name="Mookie Betts",On_Base_Percentage=0.340))## 1
## 2.157562Again, model 1 is pretty good, and model 2 is way off.
predict(lmem1,newdata=data.frame(Name="Mookie Betts",On_Base_Percentage=0.367))## 1
## 10.71174predict(lmem2,newdata=data.frame(Name="Mookie Betts",On_Base_Percentage=0.367))## 1
## 2.379796predict(lmem1,newdata=data.frame(Name="Mookie Betts",On_Base_Percentage=0.366))## 1
## 10.66158predict(lmem2,newdata=data.frame(Name="Mookie Betts",On_Base_Percentage=0.366))## 1
## 2.371565It is clear to see through all of these that model 1 is a pretty robust model for predicting stolen bases, although on a very small sample size. Model 2, whilst a good fit for the on base percentage residuals, is not a strong model. I would certainly recommend using model 1 for analysis of stolen bases from the name of the player.
Q5. Fit, evaluate, interpret a Lasso penalised regression model.
We use a Lasso penalised regression model to highlight the main pitching predictors for ERA.
ERA <- as.vector(lad_pitching$Earned_Run_Average)
lad_pitching_q5 <- lad_pitching %>%
select(-Earned_Run_Average, -Earned_Runs,-Name,-Position,-Year,-Dominant_Hand,-Switch_Hitter,-Adjusted_Earned_Run_Average)
era_predictors <- as.matrix(lad_pitching_q5)This eliminates all of the string values from the data set and saves it as a new matrix of predictors.
lprmod <- glmnet(era_predictors,ERA,alpha=1)
plot(lprmod, xvar="dev")
This visualisation of the fractional deviance shows that 80% of the deviance can be explained by 12 parameter values. To find a lambda value for this, we analyse the data frame.
lprmod##
## Call: glmnet(x = era_predictors, y = ERA, alpha = 1)
##
## Df %Dev Lambda
## 1 0 0.00 1.02100
## 2 1 10.22 0.93010
## 3 2 19.44 0.84740
## 4 2 27.93 0.77220
## 5 2 34.98 0.70360
## 6 3 40.98 0.64110
## 7 3 46.32 0.58410
## 8 3 50.75 0.53220
## 9 3 54.43 0.48490
## 10 3 57.48 0.44190
## 11 3 60.02 0.40260
## 12 4 62.46 0.36680
## 13 4 64.59 0.33420
## 14 4 66.37 0.30460
## 15 4 67.84 0.27750
## 16 4 69.07 0.25280
## 17 4 70.08 0.23040
## 18 4 70.92 0.20990
## 19 4 71.62 0.19130
## 20 5 72.29 0.17430
## 21 5 72.87 0.15880
## 22 6 73.36 0.14470
## 23 6 73.84 0.13180
## 24 6 74.24 0.12010
## 25 7 74.59 0.10950
## 26 7 74.92 0.09973
## 27 7 75.20 0.09087
## 28 7 75.43 0.08280
## 29 7 75.62 0.07544
## 30 7 75.78 0.06874
## 31 8 75.94 0.06263
## 32 8 76.13 0.05707
## 33 10 76.34 0.05200
## 34 10 76.71 0.04738
## 35 10 77.01 0.04317
## 36 11 77.48 0.03933
## 37 12 78.10 0.03584
## 38 12 79.09 0.03266
## 39 12 79.90 0.02976
## 40 12 80.57 0.02711
## 41 13 81.24 0.02470
## 42 14 81.82 0.02251
## 43 15 82.34 0.02051
## 44 15 82.86 0.01869
## 45 16 83.29 0.01703
## 46 16 83.66 0.01551
## 47 16 83.97 0.01414
## 48 17 84.23 0.01288
## 49 18 84.57 0.01174
## 50 18 84.86 0.01069
## 51 19 85.09 0.00974
## 52 19 85.29 0.00888
## 53 19 85.46 0.00809
## 54 20 85.60 0.00737
## 55 21 85.73 0.00672
## 56 21 85.83 0.00612
## 57 21 85.92 0.00558
## 58 21 85.99 0.00508
## 59 22 86.07 0.00463
## 60 23 86.23 0.00422
## 61 23 86.36 0.00384
## 62 23 86.47 0.00350
## 63 23 86.56 0.00319
## 64 23 86.63 0.00291
## 65 23 86.70 0.00265
## 66 22 86.74 0.00241
## 67 22 86.78 0.00220
## 68 22 86.81 0.00200
## 69 22 86.83 0.00183
## 70 23 86.85 0.00166
## 71 23 86.87 0.00152
## 72 23 86.88 0.00138
## 73 24 86.90 0.00126
## 74 24 86.92 0.00115
## 75 24 86.93 0.00104
## 76 24 86.95 0.00095
## 77 24 86.96 0.00087
## 78 25 86.97 0.00079
## 79 25 86.98 0.00072
## 80 27 86.98 0.00066
## 81 27 86.99 0.00060
## 82 27 87.05 0.00054
## 83 27 87.10 0.00050
## 84 27 87.10 0.00045So we select a lambda value of 0.02976.
era_12 <- coef(lprmod, s=0.02976)
era_12@Dimnames[[1]][1+era_12@i]## [1] "(Intercept)" "Rank"
## [3] "Winning_Percentage" "Games_Played"
## [5] "Hits_Allowed" "Runs_Allowed"
## [7] "Walks_Allowed" "Intentional_Walks_Allowed"
## [9] "Fielding_Independent_Pitching" "Walks_Hits_Per_Inning_Pitched"
## [11] "Hits_Per_Nine_Innings" "Home_Runs_Per_Nine_Innings"
## [13] "Strikeouts_Per_Nine_Innings"This highlights 12 main predictors for the ERA. Now, we want to select an optimal range of values for lambda.
era_cv <- cv.glmnet(era_predictors, ERA)
plot(era_cv)
This suggets a log(lambda) value between -7.5 and -5.2 is optimal. Firstly, we analyse the higher value of lambda to find the most significant predictors in the model.
era_min <- coef(lprmod,s=era_cv$lambda.1se)
setdiff(era_min@Dimnames[[1]][1+era_min@i],era_12@Dimnames[[1]][1+era_12@i])## [1] "Wins" "Losses" "Complete_Games"
## [4] "Shutouts" "Saves" "Innings_Pitched"
## [7] "Home_Runs_Allowed" "Strikeouts" "Hit_By_Pitch"
## [10] "Wild_Pitches" "Strikeout_Walk_Ratio"Then we analyse the maximum end of the ideal range.
era_max <- coef(lprmod,s=era_cv$lambda.min)
era_max@Dimnames[[1]][1+era_max@i]## [1] "(Intercept)" "Rank"
## [3] "Age" "Wins"
## [5] "Losses" "Winning_Percentage"
## [7] "Games_Played" "Games_Started"
## [9] "Complete_Games" "Shutouts"
## [11] "Saves" "Innings_Pitched"
## [13] "Hits_Allowed" "Runs_Allowed"
## [15] "Home_Runs_Allowed" "Walks_Allowed"
## [17] "Intentional_Walks_Allowed" "Strikeouts"
## [19] "Hit_By_Pitch" "Balks"
## [21] "Wild_Pitches" "Batters_Faced"
## [23] "Fielding_Independent_Pitching" "Walks_Hits_Per_Inning_Pitched"
## [25] "Hits_Per_Nine_Innings" "Home_Runs_Per_Nine_Innings"
## [27] "Strikeouts_Per_Nine_Innings" "Strikeout_Walk_Ratio"Now that we have highlighted the significant parameters for modelling ERA, we fit and evaluate the final model.
set.seed(321)
training.samples <- createDataPartition(lad_pitching$Earned_Run_Average, p=0.8, list=FALSE)
train.data <- lad_pitching[training.samples, ]
test.data <- lad_pitching[-training.samples, ]For the training model we use the 12 main predictors from earlier. Rank is modelled as a random effect as it is a discrete value which takes a lot of levels, whilst still significant in predicting the ERA.
train.model <- lmer(Earned_Run_Average ~ (1|Rank) + Winning_Percentage + Games_Played + Hits_Allowed + Runs_Allowed + Home_Runs_Allowed + Walks_Allowed + Intentional_Walks_Allowed + Fielding_Independent_Pitching + Walks_Hits_Per_Inning_Pitched + Home_Runs_Per_Nine_Innings + Strikeouts_Per_Nine_Innings, data=train.data)
predictions <- predict(train.model, test.data)
data.frame( R2 = R2(predictions, test.data$Earned_Run_Average),
RMSE = RMSE(predictions, test.data$Earned_Run_Average),
MAE = MAE(predictions, test.data$Earned_Run_Average))## R2 RMSE MAE
## 1 0.8831565 0.4650469 0.3556758This produces an R2 value of 0.8831565 which is excellent and shows a very strong correlation of 0.93976.
Q6. Using a clustering method seen in class, explore whether the data tend to form clusters and describe your findings.
We want to investigate to see if any clusters form within the batting data set.
lad_batting$Switch_Hitter <- as.factor(lad_batting$Switch_Hitter)
lad_batting$PositionGroup <- as.factor(lad_batting$PositionGroup)
lad_batting$Dominant_Hand <- as.factor(lad_batting$Dominant_Hand)
lad_batting_q6 <- lad_batting %>%
select(-Name,-Position) %>%
mutate(nPosition_Group = as.numeric(recode(PositionGroup,
"Catcher"=1,
"Infielder"=2,
"Outfielder"=3,
"Designated Hitter"=4,
"Pitcher"=5))) %>%
mutate(nDominant_Hand = as.numeric(recode(Dominant_Hand,
"Left"=1,
"Right"=2))) %>%
mutate(nSwitch_Hitter = as.numeric(recode(Switch_Hitter,
"No"=1,
"Yes"=2)))The discrete string variables are re-coded as numerical so that they can be analysed in the model. The ordering of these numbers is not important.
gower.dist <- vegdist(model.matrix(~.,lad_batting_q6[,]),method="gower")## Warning in decostand(x, "range", 2, na.rm = TRUE, ...): result contains NaN, perhaps due to impossible mathematical
## operationaggl.clust.c <- hclust(gower.dist, method="complete")
plot(aggl.clust.c, main="Agglomerative, complete linkages")
abline(h=0.59,col=2)
abline(h=0.43,col=3)
The red cut-line forms 2 clusters, whilst the green cut-line forms 4
clusters.
lad_batting_q6$red_clusters <- as.factor(cutree(aggl.clust.c,h=0.59))
lad_batting_q6$green_clusters <- as.factor(cutree(aggl.clust.c,h=0.43))
summary(lad_batting_q6$red_clusters)## 1 2
## 204 765The first cluster has 204 entries, and the second has 765.
summary(lad_batting_q6$green_clusters)## 1 2 3 4
## 204 337 289 139Now we want to have a look at the breakdown of these clusters based on the discrete variables that we re-coded earlier. This is a breakdown of the clusters by position group.
prop.table(table(lad_batting_q6$nPosition_Group,lad_batting_q6$red_clusters))##
## 1 2
## 1 0.017543860 0.115583075
## 2 0.136222910 0.443756450
## 3 0.054695562 0.201238390
## 4 0.002063983 0.001031992
## 5 0.000000000 0.027863777prop.table(table(lad_batting_q6$nPosition_Group,lad_batting_q6$green_clusters))##
## 1 2 3 4
## 1 0.017543860 0.039215686 0.065015480 0.011351909
## 2 0.136222910 0.167182663 0.162022704 0.114551084
## 3 0.054695562 0.141382869 0.042311662 0.017543860
## 4 0.002063983 0.000000000 0.001031992 0.000000000
## 5 0.000000000 0.000000000 0.027863777 0.000000000As a reminder, “Catcher”=1, “Infielder”=2,“Outfielder”=3,“Designated Hitter”=“4”,“Pitcher”=“5”. Some interesting points are that all of the pitchers belong to cluster number 3, and all the designated hitters belong to clusters 1 and 3. The rest of hte positions are spread relatively evenly across all four.
Now we look at whether a player is a switch hitter or not. A switch hitter can hit with both left and right hand.
prop.table(table(lad_batting_q6$nSwitch_Hitter,lad_batting_q6$red_clusters))##
## 1 2
## 1 0.2105263 0.6460268
## 2 0.0000000 0.1434469prop.table(table(lad_batting_q6$nSwitch_Hitter,lad_batting_q6$green_clusters))##
## 1 2 3 4
## 1 0.2105263 0.3477812 0.2982456 0.0000000
## 2 0.0000000 0.0000000 0.0000000 0.1434469This clustering shows that all of the switch hitters belong to cluster number 4, whilst none of the non switch hitters belong to cluster 4.
Now we have a look at dominant hand.
prop.table(table(lad_batting_q6$nDominant_Hand,lad_batting_q6$red_clusters))##
## 1 2
## 1 0.1001032 0.2156863
## 2 0.1104231 0.5737874prop.table(table(lad_batting_q6$nDominant_Hand,lad_batting_q6$green_clusters))##
## 1 2 3 4
## 1 0.10010320 0.13415893 0.08152735 0.00000000
## 2 0.11042312 0.21362229 0.21671827 0.14344685Here we see that cluster number 4 contains exclusively right handed players. This shows an interesting point in the data, that all switch hitters are recorded as right handers. It also means that no pitchers or designated hitters are recorded as switch hitters in the data set. It also means that clusters 1, 2 and 3 contain a larger percentage of left handed players then observed in the data set.
Now we see the silhouette plots of the red and green clusters.
sil=silhouette(as.numeric(lad_batting_q6$red_clusters),gower.dist)
plot(sil,col=1:2,border=NA)
abline(v=mean(sil[,3]),lty=2)
sil=silhouette(as.numeric(lad_batting_q6$green_clusters),gower.dist)
plot(sil,col=1:4,border=NA)
abline(v=mean(sil[,3]),lty=2)
To better visualise the effect of adding clusters, we show a mean
silohette plot.
silav<-2:18
for(i in 2:18){
sili=silhouette(cutree(aggl.clust.c,k=i),gower.dist)
silav[i]<-mean(sili[,3])
}
plot(2:18,silav[2:18],type="l",ylab="Correlation Score Between Data Points",xlab="Number of Clusters")
The mean silhouette plot shows that correlation between data points is
at a maximum with 2 clusters, and falls to minimum with 3 clusters.
However, the correlation increases up until 15 clusters. This suggests
that with 15 clusters, there is quite a good level of correlation
between data points in each cluster.