A linear regression with one untransformed predictor fits a straight line. Transformations allow curved relationships, while multiple regression adds predictors. All of these models use lm().

Statcast Fly Ball Data

The files statcast_models_train.csv and statcast_models_test.csv contain a subset of fly balls from the 2019 MLB regular season with complete Statcast measurements. The subset includes exit velocities from 66 to 111 miles per hour and launch angles from 20 to 60 degrees. The training data cover March through July, and the testing data cover August and September.

Each row is one fly ball. The models use three variables:

  • exit_velocity: speed off the bat, in miles per hour
  • launch_angle: vertical angle off the bat, in degrees
  • hit_distance: projected distance traveled, in feet

We will fit each model with the training data and leave the testing data unopened until the candidates are fixed. Testing RMSE will then measure how well the fitted relationships carry forward to later observations. Let’s load the training data first.

statcast_train <- read_csv(
  "data/statcast_models_train.csv",
  show_col_types = FALSE
)

Let’s visualize a random sample of 10,000 fly balls so that the point cloud remains visible. We will still fit the models and calculate RMSE using every training fly ball.

statcast_plot <- statcast_train %>%
  slice_sample(n = 10000)
ggplot(statcast_plot, aes(x = exit_velocity, y = hit_distance)) +
  geom_point(alpha = 0.15) +
  labs(x = "Exit Velocity (mph)", y = "Hit Distance (feet)") +
  theme_minimal()

Fly balls hit harder generally travel farther, but a straight line is only one possible specification. We will now consider logarithmic and quadratic alternatives.

Exit Velocity Transformations

The linear model assumes that one additional mile per hour produces the same fitted difference at 70 mph as at 100 mph. The logarithmic and quadratic models allow that difference to change. Let’s fit all three specifications with the training data.

lin_velo_model <- lm(
  hit_distance ~ exit_velocity,
  data = statcast_train
)

log_velo_model <- lm(
  hit_distance ~ log(exit_velocity),
  data = statcast_train
)

quad_velo_model <- lm(
  hit_distance ~ exit_velocity + I(exit_velocity^2),
  data = statcast_train
)

The logarithm compresses larger exit velocities more than smaller ones, so its fitted curve flattens as exit velocity increases. The quadratic model allows the slope to increase or decrease.

The linear, logarithmic, and quadratic forms are a + b * x, a + b * log(x), and a + b1 * x + b2 * x^2. Inside a model formula, I() tells R to perform ordinary arithmetic.

Comparing the Exit Velocity Curves

To compare the fitted shapes directly, let’s predict each model at the same evenly spaced exit velocity values. Keeping the exit velocity values fixed makes the differences among the curves easier to see.

speed_grid <- tibble(
  exit_velocity = seq(
    min(statcast_train$exit_velocity),
    max(statcast_train$exit_velocity),
    by = 0.5
  )
) %>%
  mutate(
    lin_velo = predict(lin_velo_model, newdata = .),
    log_velo = predict(log_velo_model, newdata = .),
    quad_velo = predict(quad_velo_model, newdata = .)
  ) %>%
  pivot_longer(
    cols = c(lin_velo, log_velo, quad_velo),
    names_to = "model",
    values_to = "fitted_distance"
  ) %>%
  mutate(
    model = factor(
      model,
      levels = c("lin_velo", "log_velo", "quad_velo"),
      labels = c("Lin. velo.", "Log velo.", "Quad. velo.")
    )
  )

pivot_longer() places the three sets of predictions in one column and records their source in model. The resulting long table lets us map model to color and draw all three curves together.

ggplot() +
  geom_point(
    data = statcast_plot,
    mapping = aes(x = exit_velocity, y = hit_distance),
    color = "gray40",
    alpha = 0.1,
    inherit.aes = FALSE
  ) +
  geom_line(
    data = speed_grid,
    mapping = aes(
      x = exit_velocity,
      y = fitted_distance,
      color = model
    ),
    linewidth = 1
  ) +
  labs(
    x = "Exit Velocity (mph)",
    y = "Hit Distance (feet)",
    color = NULL
  ) +
  theme_minimal()

The linear curve has a constant slope, the logarithmic curve flattens, and the quadratic curve becomes steeper.

The curves compare shape, not accuracy. Let’s calculate training RMSE, the typical prediction error in feet; lower is better. We will add fitted values for every training fly ball, define rmse(), and use .by = model inside reframe() to calculate one score per model.

statcast_train <- statcast_train %>%
  mutate(
    lin_velo = predict(lin_velo_model),
    log_velo = predict(log_velo_model),
    quad_velo = predict(quad_velo_model)
  )
rmse <- function(predicted, actual) {
  sqrt(mean((actual - predicted)^2))
}

speed_training_scores <- statcast_train %>%
  pivot_longer(
    cols = c(lin_velo, log_velo, quad_velo),
    names_to = "model",
    values_to = "prediction"
  ) %>%
  reframe(
    rmse = rmse(prediction, hit_distance),
    .by = model
  )

speed_training_scores
## # A tibble: 3 × 2
##   model      rmse
##   <chr>     <dbl>
## 1 lin_velo  31.70
## 2 log_velo  32.68
## 3 quad_velo 30.69

lm() chooses coefficients to reduce error on the training data. The quadratic model can reproduce the linear model by setting the coefficient on the squared term to zero, so its training RMSE cannot be larger. Testing data determine whether that added flexibility improves predictions on later observations.

Adding Launch Angle

Training RMSE ranks the three exit velocity specifications, but it does not show what the models are missing. Let’s plot the residuals from lin_velo_model against launch angle. Each residual is observed hit distance minus fitted hit distance. A positive residual means the fly ball traveled farther than fitted; a negative residual means it traveled less far.

statcast_residual_plot <- statcast_train %>%
  mutate(residual_distance = hit_distance - lin_velo) %>%
  slice_sample(n = 10000)
ggplot(
  statcast_residual_plot,
  aes(x = launch_angle, y = residual_distance)
) +
  geom_hline(yintercept = 0, color = "gray50") +
  geom_point(alpha = 0.15) +
  geom_smooth(
    method = "lm",
    formula = y ~ x + I(x^2),
    se = FALSE,
    color = "lightcoral",
    linewidth = 1
  ) +
  labs(
    x = "Launch Angle (degrees)",
    y = "Residual Hit Distance (feet)"
  ) +
  theme_minimal()

The residuals rise and then fall. The model using only exit velocity underpredicts many fly balls near the middle of the angle range and overpredicts many fly balls with high launch angles.

This pattern motivates a quadratic launch angle specification using launch_angle and I(launch_angle^2). Let’s add that same angle specification to each exit velocity model. Holding the angle terms fixed keeps the comparison focused on the exit velocity transformation.

lin_velo_quad_angle_model <- lm(
  hit_distance ~ exit_velocity + launch_angle + I(launch_angle^2),
  data = statcast_train
)

log_velo_quad_angle_model <- lm(
  hit_distance ~ log(exit_velocity) + launch_angle + I(launch_angle^2),
  data = statcast_train
)

quad_velo_quad_angle_model <- lm(
  hit_distance ~ exit_velocity + I(exit_velocity^2) +
    launch_angle + I(launch_angle^2),
  data = statcast_train
)

These are multiple regressions because each model uses exit velocity and launch angle. The exit velocity terms compare fly balls with the same launch angle. The two launch angle terms allow fitted distance to rise and then fall as launch angle increases. Let’s add fitted values from these models to the training data.

statcast_train <- statcast_train %>%
  mutate(
    lin_velo_quad_angle = predict(
      lin_velo_quad_angle_model
    ),
    log_velo_quad_angle = predict(
      log_velo_quad_angle_model
    ),
    quad_velo_quad_angle = predict(
      quad_velo_quad_angle_model
    )
  )

Each multiple regression contains its counterpart using only exit velocity when both launch angle coefficients are zero. Adding launch angle therefore cannot increase training RMSE. Testing RMSE determines whether the improvement carries forward.

Selecting a Model with Testing Data

The six candidates are now fixed: three use exit velocity alone, and three add quadratic launch angle. Let’s load the later months and generate predictions without refitting. Every testing prediction will come from a relationship learned through July.

statcast_test <- read_csv(
  "data/statcast_models_test.csv",
  show_col_types = FALSE
) %>%
  mutate(
    lin_velo = predict(lin_velo_model, newdata = .),
    log_velo = predict(log_velo_model, newdata = .),
    quad_velo = predict(quad_velo_model, newdata = .),
    lin_velo_quad_angle = predict(
      lin_velo_quad_angle_model,
      newdata = .
    ),
    log_velo_quad_angle = predict(
      log_velo_quad_angle_model,
      newdata = .
    ),
    quad_velo_quad_angle = predict(
      quad_velo_quad_angle_model,
      newdata = .
    )
  )

To compare training and testing error directly, let’s stack the two samples, reshape the six prediction columns, and calculate RMSE for each combination of sample and model. We will then plot the training and testing scores next to each other.

model_comparison <- bind_rows(
  Training = statcast_train,
  Testing = statcast_test,
  .id = "sample"
) %>%
  pivot_longer(
    cols = c(
      lin_velo, log_velo, quad_velo,
      lin_velo_quad_angle, log_velo_quad_angle,
      quad_velo_quad_angle
    ),
    names_to = "model",
    values_to = "prediction"
  ) %>%
  reframe(
    rmse = rmse(prediction, hit_distance),
    .by = c(sample, model)
  ) %>%
  mutate(
    sample = factor(sample, levels = c("Training", "Testing")),
    model = factor(
      model,
      levels = c(
        "lin_velo",
        "log_velo",
        "quad_velo",
        "lin_velo_quad_angle",
        "log_velo_quad_angle",
        "quad_velo_quad_angle"
      ),
      labels = c(
        "Lin. velo.",
        "Log velo.",
        "Quad. velo.",
        "Lin. velo. + quad. angle",
        "Log velo. + quad. angle",
        "Quad. velo. + quad. angle"
      )
    )
  )

model_comparison
## # A tibble: 12 × 3
##    sample   model                      rmse
##    <fct>    <fct>                     <dbl>
##  1 Training Lin. velo.                31.70
##  2 Training Log velo.                 32.68
##  3 Training Quad. velo.               30.69
##  4 Training Lin. velo. + quad. angle  15.14
##  5 Training Log velo. + quad. angle   15.46
##  6 Training Quad. velo. + quad. angle 15.12
##  7 Testing  Lin. velo.                31.85
##  8 Testing  Log velo.                 32.84
##  9 Testing  Quad. velo.               30.82
## 10 Testing  Lin. velo. + quad. angle  14.96
## 11 Testing  Log velo. + quad. angle   15.29
## 12 Testing  Quad. velo. + quad. angle 14.95
ggplot(
  model_comparison,
  aes(x = model, y = rmse, fill = sample)
) +
  geom_col(position = position_dodge(width = 0.9)) +
  geom_text(
    aes(label = scales::number(rmse, accuracy = 0.01)),
    position = position_dodge(width = 0.9),
    vjust = -0.4,
    size = 3
  ) +
  scale_x_discrete(
    labels = scales::label_wrap(18),
    guide = guide_axis(n.dodge = 2)
  ) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  labs(x = NULL, y = "RMSE (feet)", fill = NULL) +
  theme_minimal()

Training RMSE measures fit to the observations used to estimate a model. Testing RMSE measures errors on later observations and provides the basis for selection.

Adding the curved launch angle relationship reduces testing RMSE from more than 30 feet to less than 15.3 feet. Launch angle explains substantial, repeatable information left by the models using only exit velocity.

Once launch angle enters the model, the exit velocity shape matters little. The quadratic velocity model has the smallest testing RMSE, but it beats the linear velocity model by only about 0.015 feet. The large improvement comes from adding launch angle, not from choosing a more complicated exit velocity shape.