Today, we will return to the batting dataset from Problem Set 3, reshape it in a way that is useful for plotting and modeling, and then introduce the ideas that lead to linear regression.

If you have not yet completed the challenge at the end of Problem Set 3, you can download the final batting_2014_2015 dataset here and save it to your data folder.

load("data/batting_2014_2015.RData")
batting_2014_2015 <- batting_2014_2015 %>% ungroup()
batting_2014_2015
## # A tibble: 140 × 3
##    playerID  yearID    BA
##    <chr>      <int> <dbl>
##  1 abreujo02   2014 0.317
##  2 abreujo02   2015 0.290
##  3 altuvjo01   2014 0.341
##  4 altuvjo01   2015 0.313
##  5 andruel01   2014 0.263
##  6 andruel01   2015 0.258
##  7 aybarer01   2014 0.278
##  8 aybarer01   2015 0.270
##  9 bautijo02   2014 0.286
## 10 bautijo02   2015 0.250
## # ℹ 130 more rows

Tidy Data

Many analyses become easier once the data are in a tidy format:

  1. Each variable must have its own column.
  2. Each observation must have its own row.
  3. Each value must have its own cell.

Right now, our tbl has one row per player-season. To compare a player’s 2014 and 2015 batting averages side by side, it will be more useful to have one row per player and a separate column for each season. We can do that with pivot_wider(). The names_from argument tells R which column should supply the new column names, and values_from tells R which values should fill those new columns.

batting_2014_2015 <- batting_2014_2015 %>%
  pivot_wider(
    names_from = yearID,
    values_from = BA,
    names_prefix = "BA_"
  )
head(batting_2014_2015)
## # A tibble: 6 × 3
##   playerID  BA_2014 BA_2015
##   <chr>       <dbl>   <dbl>
## 1 abreujo02   0.317   0.290
## 2 altuvjo01   0.341   0.313
## 3 andruel01   0.263   0.258
## 4 aybarer01   0.278   0.270
## 5 bautijo02   0.286   0.250
## 6 beltrad01   0.324   0.287

Sometimes the long format is more useful for plotting. We can return to that structure with pivot_longer(), which takes several columns and stacks them into a pair of columns.

batting_2014_2015_long <- batting_2014_2015 %>%
  pivot_longer(
    cols = starts_with("BA_"),
    names_to = "yearID",
    values_to = "BA"
  ) %>%
  mutate(yearID = sub("BA_", "", yearID))

head(batting_2014_2015_long)
## # A tibble: 6 × 3
##   playerID  yearID    BA
##   <chr>     <chr>  <dbl>
## 1 abreujo02 2014   0.317
## 2 abreujo02 2015   0.290
## 3 altuvjo01 2014   0.341
## 4 altuvjo01 2015   0.313
## 5 andruel01 2014   0.263
## 6 andruel01 2015   0.258

We can use this long table to make a side-by-side bar plot of batting averages in 2014 and 2015.

batting_2014_2015_slice <- batting_2014_2015_long %>%
  slice_head(n = 20)

ggplot(batting_2014_2015_slice, aes(x = playerID, y = BA, fill = yearID)) +
  geom_col(position = "dodge", color = "black", alpha = 0.8) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

To model the relationship between the two seasons, we will return to the wide version of the data, where each row is one player.

Fitting Batting Averages

We know each player’s batting average in both 2014 and 2015. We will use BA_2014 to group players with similar earlier-season performance, then assign each player a fitted 2015 value based on the observed BA_2015 values in that group. Comparing the fitted values with the observed 2015 values measures how closely a method fits these players. It does not measure how accurately the method would forecast a future season or a new player.

The simplest model ignores BA_2014 and assigns every player the overall mean of BA_2015. This mean is 0.273, whereas the mean 2014 batting average is 0.272. We will save these means as avg_2014 and avg_2015 for use in our plots. We can do this by using the $ operator to feed an entire column into mean().

batting_2014_2015 <-
  batting_2014_2015 %>% 
  mutate(yhat_0 = mean(BA_2015))
batting_2014_2015
## # A tibble: 70 × 4
##    playerID  BA_2014 BA_2015 yhat_0
##    <chr>       <dbl>   <dbl>  <dbl>
##  1 abreujo02   0.317   0.290  0.273
##  2 altuvjo01   0.341   0.313  0.273
##  3 andruel01   0.263   0.258  0.273
##  4 aybarer01   0.278   0.270  0.273
##  5 bautijo02   0.286   0.250  0.273
##  6 beltrad01   0.324   0.287  0.273
##  7 blackch02   0.288   0.287  0.273
##  8 bogaexa01   0.240   0.320  0.273
##  9 brantmi02   0.327   0.310  0.273
## 10 braunry02   0.266   0.285  0.273
## # ℹ 60 more rows
avg_2014 = mean(batting_2014_2015$BA_2014)
avg_2015 = mean(batting_2014_2015$BA_2015)

The overall mean provides a baseline fit in which every player receives the same fitted value. Plotting the 2015 batting averages against the 2014 batting averages allows us to assess visually whether there is a relationship that a more flexible model could capture. We have added dashed red horizontal and vertical lines at the overall means of the 2014 and 2015 data, respectively.

ggplot(batting_2014_2015) +
  geom_point(aes(x = BA_2014, y = BA_2015), col = 'black', size = 1) +
  geom_hline(yintercept = avg_2015, col = 'red', lty = 2) +
  geom_vline(xintercept = avg_2014, col = 'red', lty = 2)

It certainly looks like there is a relationship! Using the 2014 batting averages may allow us to fit the observed 2015 values more closely than assigning every player the same 2015 average.

Looking carefully at the plot, we see that most players with below average batting averages in 2014 tended to also have below average batting averages in 2015. Similarly, most players with above average batting averages in 2014 tended to have above average batting averages in 2015.

One way to improve on the baseline fit would be as follows:

  1. Divide the players into two groups, one for those with above average BA in 2014 and one for those with below average BA in 2014.
  2. Average the 2015 BA within each group and use that average as the fitted value for every member of the group.

In order to do this, we will need to mutate our tbl batting_2014_2015 with a column indicating to which group each player belongs. Then we can pass this column to group_by() and compute the average BA_2015 within each group. To create the column indicating group membership, we will use the powerful cut() function, which divides the range of a numerical vector into intervals and recodes the numerical values according to which interval they fall.

The following code chunk does precisely that with two more steps: once we have calculated the fitted values, we ungroup the tbl and drop the column indicating the interval in which each observation falls.

batting_2014_2015 <- 
  batting_2014_2015 %>%
  mutate(bins = cut(BA_2014, breaks = c(0.15, 0.272, 0.40))) %>%
  group_by(bins) %>%
  mutate(yhat_1 = mean(BA_2015)) %>%
  ungroup() %>%
  select(-bins)
batting_2014_2015
## # A tibble: 70 × 5
##    playerID  BA_2014 BA_2015 yhat_0 yhat_1
##    <chr>       <dbl>   <dbl>  <dbl>  <dbl>
##  1 abreujo02   0.317   0.290  0.273  0.281
##  2 altuvjo01   0.341   0.313  0.273  0.281
##  3 andruel01   0.263   0.258  0.273  0.265
##  4 aybarer01   0.278   0.270  0.273  0.281
##  5 bautijo02   0.286   0.250  0.273  0.281
##  6 beltrad01   0.324   0.287  0.273  0.281
##  7 blackch02   0.288   0.287  0.273  0.281
##  8 bogaexa01   0.240   0.320  0.273  0.265
##  9 brantmi02   0.327   0.310  0.273  0.281
## 10 braunry02   0.266   0.285  0.273  0.265
## # ℹ 60 more rows

When we run this code and print out our tbl, we see that there is a new column called yhat_1 that contains our new fitted values. Before proceeding, we should talk a bit about the syntax used in cut(). The first argument is the variable we want to discretize. The next argument, breaks =, is a vector that tells R where the endpoints of these intervals are. These are often called “cut points.” In this particular case, we wanted to divide the players into those with below-average batting averages in 2014 and above-average batting averages in 2014. The first element of the cut point vector, 0.15, is much less than the smallest BA_2014 value, whereas the second element, 0.272, is the overall mean of the BA_2014 values. The last element, 0.40, is much greater than the largest BA_2014 value.

Now that we have two different ways of fitting the 2015 batting averages, let us compare them visually.

ggplot(batting_2014_2015) +
  geom_point(aes(x = BA_2014, y = BA_2015), col = 'black', size = 2) +
  geom_point(aes(x = BA_2014, y = yhat_0), col = 'lightcoral', size = 2) + 
  geom_point(aes(x = BA_2014, y = yhat_1), col = 'lightblue', size = 2)

Visually, the blue points (corresponding to yhat_1) appear a bit closer to the observed values than the red points (corresponding to yhat_0). Dividing the players into two bins according to their 2014 batting average therefore fits these data better than assigning every player the overall average BA_2015 value.

Of course, we can continue with this line of reasoning and divide the players into even more bins. When we do that, instead of hand-coding the vector of cut points, we can use the function seq(), which generates a vector of equally spaced numbers. To demonstrate, suppose we wanted to divide the interval [0,1] into 10 equally sized intervals: (0,0.1], (0.1, 0.2], …, (0.9, 1]. To get the vector of cut points, we need to tell seq() either how many points we wanted or the spacing between the points:

seq(from = 0, to = 1, length = 11)
##  [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
seq(from = 0, to = 1, by = 0.1)
##  [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Suppose we divide the 2014 batting averages into intervals of length 0.05 and use the average BA_2015 within each bin as its fitted value. We could run the following:

batting_2014_2015 <- 
  batting_2014_2015 %>%
  mutate(bins = cut(BA_2014, breaks = seq(from = 0.15, to = 0.4, by = 0.05))) %>%
  group_by(bins) %>%
  mutate(yhat_2 = mean(BA_2015)) %>%
  ungroup() %>%
  select(-bins)
batting_2014_2015
## # A tibble: 70 × 6
##    playerID  BA_2014 BA_2015 yhat_0 yhat_1 yhat_2
##    <chr>       <dbl>   <dbl>  <dbl>  <dbl>  <dbl>
##  1 abreujo02   0.317   0.290  0.273  0.281  0.300
##  2 altuvjo01   0.341   0.313  0.273  0.281  0.300
##  3 andruel01   0.263   0.258  0.273  0.265  0.271
##  4 aybarer01   0.278   0.270  0.273  0.281  0.271
##  5 bautijo02   0.286   0.250  0.273  0.281  0.271
##  6 beltrad01   0.324   0.287  0.273  0.281  0.300
##  7 blackch02   0.288   0.287  0.273  0.281  0.271
##  8 bogaexa01   0.240   0.320  0.273  0.265  0.257
##  9 brantmi02   0.327   0.310  0.273  0.281  0.300
## 10 braunry02   0.266   0.285  0.273  0.265  0.271
## # ℹ 60 more rows

We can also visualize our new fitted values, this time with green points.

ggplot(batting_2014_2015) +
  geom_point(aes(x = BA_2014, y = BA_2015), size = 2) + 
  geom_point(aes(x = BA_2014, y = yhat_0), col = 'lightcoral', size = 2, alpha = 0.8) + 
  geom_point(aes(x = BA_2014, y = yhat_1), col = 'lightblue',  size = 2, alpha = 0.8) +
  geom_point(aes(x = BA_2014, y = yhat_2), col = 'darkseagreen3', size = 2, alpha = 0.8)

The fitted value exactly matches the observed 2015 batting average of the player with the lowest 2014 batting average. Why do you think this happened?

Let’s add two more fitted-value columns by dividing the 2014 batting averages into bins of length 0.005 and 0.001. The code below saves them as yhat_3 and yhat_4, then adds them to the plot in pink and orange.

batting_2014_2015 <- 
  batting_2014_2015 %>%
  mutate(bins = cut(BA_2014, breaks = seq(from = 0.15, to = 0.4, by = 0.005))) %>%
  group_by(bins) %>%
  mutate(yhat_3 = mean(BA_2015)) %>%
  ungroup() %>%
  select(-bins)

batting_2014_2015 <- 
  batting_2014_2015 %>%
  mutate(bins = cut(BA_2014, breaks = seq(from = 0.15, to = 0.4, by = 0.001))) %>%
  group_by(bins) %>%
  mutate(yhat_4 = mean(BA_2015)) %>%
  ungroup() %>%
  select(-bins)

ggplot(batting_2014_2015) +
  geom_point(aes(x = BA_2014, y = BA_2015), size = 2) + 
  geom_point(aes(x = BA_2014, y = yhat_0), col = 'lightcoral', size = 2) + 
  geom_point(aes(x = BA_2014, y = yhat_1), col = 'lightblue', size = 2) +
  geom_point(aes(x = BA_2014, y = yhat_2), col = 'darkseagreen3', size = 2) +
  geom_point(aes(x = BA_2014, y = yhat_3), col = 'hotpink', size = 2) + 
  geom_point(aes(x = BA_2014, y = yhat_4), col = 'darkorange', size = 1.5, alpha = .7)

Here, we can see that some of the orange points are directly on top of the actual values.

Assessing Fit with RMSE

We now have several ways of fitting BA_2015. Qualitatively, the orange fitted values (formed by binning BA_2014 into very small intervals) appear to fit the observed data better than the green, blue, and red fitted values. To assess fit quantitatively, we often rely on root mean square error, or RMSE. RMSE is the square root of the mean square error (MSE), which averages the squared differences between the observed and fitted values.

reframe(batting_2014_2015, 
          rmse_0 = sqrt(mean((BA_2015 - yhat_0)^2)),
          rmse_1 = sqrt(mean((BA_2015 - yhat_1)^2)),
          rmse_2 = sqrt(mean((BA_2015 - yhat_2)^2)),
          rmse_3 = sqrt(mean((BA_2015 - yhat_3)^2)),
          rmse_4 = sqrt(mean((BA_2015 - yhat_4)^2)))
## # A tibble: 1 × 5
##   rmse_0 rmse_1 rmse_2 rmse_3 rmse_4
##    <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
## 1 0.0257 0.0244 0.0226 0.0185 0.0113

The orange fitted values have the lowest RMSE, followed by the pink, green, blue, and red fitted values. Binning into smaller intervals gives the model more flexibility to follow the observed data. Because these RMSE values are calculated on the same players used to form the bin averages, they measure training fit rather than performance on new players or a future season.

Why not keep shrinking the bins? A bin containing one player reproduces that player’s 2015 batting average exactly, reducing training error without establishing a pattern that must hold for anyone else. Empty bins also provide no fitted value. A straight line avoids these abrupt jumps and lets all of the observed players inform the fit across the range of 2014 batting averages. In Lecture 6, we will use linear regression to fit that line.