MLB Batting Statistics

Let’s first load the Lahman library used in Lecture 4.

library(Lahman)
batting <- as_tibble(Batting)

In the code below, we’ll create a dataset containing player batting averages from the 2017 and 2018 seasons. We keep players with at least 300 at-bats in both seasons and reshape the data so that each player has one row and each season has its own column.

batting_2017_2018 <- batting %>%
  filter(yearID %in% c(2017, 2018)) %>%
  group_by(playerID, yearID) %>%
  reframe(H = sum(H), AB = sum(AB)) %>%
  filter(AB >= 300) %>%
  mutate(BA = H / AB) %>%
  ungroup() %>%
  group_by(playerID) %>%
  filter(n() == 2) %>%
  ungroup() %>%
  select(playerID, yearID, BA) %>%
  pivot_wider(
    names_from = yearID,
    values_from = BA,
    names_prefix = "BA_"
  ) %>%
  arrange(playerID)

batting_2017_2018
## # A tibble: 196 × 3
##    playerID  BA_2017 BA_2018
##    <chr>       <dbl>   <dbl>
##  1 abreujo02   0.304   0.265
##  2 adamsma01   0.274   0.239
##  3 alonsyo01   0.266   0.25 
##  4 altuvjo01   0.346   0.316
##  5 anderti01   0.257   0.240
##  6 andruel01   0.297   0.256
##  7 arciaor01   0.277   0.236
##  8 arenano01   0.309   0.297
##  9 baezja01    0.273   0.290
## 10 barnhtu01   0.270   0.248
## # ℹ 186 more rows

Predicting Batting Averages

We start by plotting each player’s 2017 batting average against their 2018 batting average.

ba_plot <- batting_2017_2018 %>%
  ggplot(aes(x = BA_2017, y = BA_2018)) +
  geom_point() +
  labs(x = "2017 Batting Average", y = "2018 Batting Average")

ba_plot

Correlation

The correlation coefficient, rr, summarizes the linear relationship between two variables. It ranges from -1 to 1: values near 1 indicate a strong positive linear relationship, values near -1 indicate a strong negative linear relationship, and values near 0 indicate a weak linear relationship.

We can calculate the correlation between the 2017 and 2018 batting averages with cor():

cor(batting_2017_2018$BA_2017, batting_2017_2018$BA_2018)
## [1] 0.4795003

The two seasons have a moderate positive correlation, but there is still a lot of variation from player to player. Remember that correlation measures only a linear relationship; variables can be related in other ways even when their correlation is close to zero.

Fitting a Linear Model

We can describe the linear relationship and use it for prediction by fitting a linear model with lm():

fit <- lm(BA_2018 ~ BA_2017, data = batting_2017_2018)

The first argument is a formula in the form response ~ predictor. Here, BA_2018 is the response we want to predict and BA_2017 is the predictor. The data argument tells R where to find those variables.

The fitted model contains the intercept and slope of the regression line:

fit$coefficients
## (Intercept)     BA_2017 
##   0.1305237   0.4735976

We can use those coefficients to add the fitted line to our scatterplot.

intercept <- fit$coefficients[[1]]
slope <- fit$coefficients[[2]]

ba_plot <- ba_plot +
  geom_abline(intercept = intercept, slope = slope, color = "red") +
  labs(
    title = paste0(
      "Predicted BA = ", round(intercept, 3),
      " + ", round(slope, 3), " x 2017 BA"
    )
  )

ba_plot

Question: How would you interpret the slope? What does the model predict will happen to 2018 batting average when 2017 batting average increases by 0.010?

Predictions and Residuals

The function predict() uses the fitted model to calculate a predicted 2018 batting average for every player in the data. A residual is the difference between a player’s actual and predicted batting average.

batting_2017_2018 <- batting_2017_2018 %>%
  mutate(
    pred = predict(fit),
    resid = BA_2018 - pred
  )

batting_2017_2018
## # A tibble: 196 × 5
##    playerID  BA_2017 BA_2018  pred    resid
##    <chr>       <dbl>   <dbl> <dbl>    <dbl>
##  1 abreujo02   0.304   0.265 0.275 -0.0101 
##  2 adamsma01   0.274   0.239 0.260 -0.0219 
##  3 alonsyo01   0.266   0.25  0.257 -0.00654
##  4 altuvjo01   0.346   0.316 0.294  0.0222 
##  5 anderti01   0.257   0.240 0.252 -0.0125 
##  6 andruel01   0.297   0.256 0.271 -0.0155 
##  7 arciaor01   0.277   0.236 0.262 -0.0259 
##  8 arenano01   0.309   0.297 0.277  0.0199 
##  9 baezja01    0.273   0.290 0.260  0.0307 
## 10 barnhtu01   0.270   0.248 0.259 -0.0107 
## # ℹ 186 more rows

We can draw each residual as the vertical distance between an observed value and the fitted line.

ggplot(batting_2017_2018) +
  geom_segment(
    aes(x = BA_2017, xend = BA_2017, y = BA_2018, yend = pred),
    color = "dodgerblue"
  ) +
  geom_point(aes(x = BA_2017, y = BA_2018)) +
  geom_abline(intercept = intercept, slope = slope, color = "red") +
  labs(x = "2017 Batting Average", y = "2018 Batting Average")

Another way to summarize the fit is R2R^2, the proportion of variation in the response explained by the model. We can obtain it from the model summary:

summary(fit)$r.squared
## [1] 0.2299205

Predicting New Observations

So far, we have predicted values for the same players used to fit the model. Suppose we know three other players’ 2017 batting averages but not their 2018 batting averages. We can pass a new tibble to the newdata argument of predict().

new_data <- tibble(BA_2017 = c(0.241, 0.310, 0.265))

new_data <- new_data %>%
  mutate(pred = predict(fit, newdata = new_data))

new_data
## # A tibble: 3 × 2
##   BA_2017  pred
##     <dbl> <dbl>
## 1   0.241 0.245
## 2   0.31  0.277
## 3   0.265 0.256

Finally, we add the new predictions to the original plot.

ggplot(batting_2017_2018, aes(x = BA_2017, y = BA_2018)) +
  geom_point() +
  geom_abline(intercept = intercept, slope = slope, color = "red") +
  geom_point(
    data = new_data,
    aes(x = BA_2017, y = pred),
    color = "dodgerblue",
    size = 3
  ) +
  labs(x = "2017 Batting Average", y = "2018 Batting Average")