Tommy John Surgery and its Relationship to MLB Pitcher Career Trajectory

Age distribution

In understanding career trajectory and longevity, we should first understand the age distribution of players in the MLB. Using all pitchers that played in the 2021 season, we get the following age distribution and see most pitchers are aged between their low twenties and low thirties.

seasonal_complete %>%
  filter(Year == 2021) %>%
  ggplot(aes(x=Age))+
  geom_density(fill = "lightcyan3")+
  theme_classic()+
  labs(title = "Age distribution of pitchers in the MLB", x="Age", y="")+
  theme(plot.title.position = "plot", 
        plot.title = element_text(family = "mono", size = 11), 
        axis.title.x = element_text(family = "mono", size = 11))

Career length

In trying to understand the impact of having Tommy John surgery on a player’s career, we should seek to understand if there are differences in career length for players that receive the surgery. To do this, we compare the career length in seasons of players that did not have surgery at any point in their career to those that had surgery between their draft and debut date. We remove players that had surgery in the MLB for this graph - the fact that they made it far enough in their MLB career to have a surgery during it is an indicator of career length. Additionally, we do not include players that had surgery before being drafted as scouts deemed them worthy enough to be drafted despite their surgery. When looking at the graph, there do not seem to be any significant differences.

seasonal_complete %>%
  filter(surgery1time %!in% c("after debut", "before draft")) %>%
  mutate(surgery1time = case_when(surgery1time == "between draft and debut" ~ "between draft and debut", 
                                  TRUE ~ "no surgery")) %>%
  group_by(Player, surgery1time) %>%
  summarize(n = n(), lastyear = max(Year)) %>%
  filter(lastyear < 2020) %>%
  ggplot(aes(x=n, fill=surgery1time))+
  geom_density(alpha = 0.5)+
  theme_classic()+
  scale_fill_manual(values = c("deepskyblue4", "sienna2"))+
  labs(x="career length (seasons)", y="", title = "MLB career length (seasons) averages similar length for players that<strong><span style='color:sienna2'> didn't have Tommy John</span></strong></b>", subtitle = "<strong><span style='color:sienna2'> surgery before debuting </span></strong></b>in comparison to <strong><span style='color:deepskyblue4'>those that had surgery between draft and debut</span></strong></b>")+
  theme(plot.title = element_markdown(size = 9, family = "mono"),
        plot.subtitle = element_markdown(size = 9, family = "mono"), 
        plot.title.position = "plot",
        legend.position = "none", 
        axis.title.x = element_text(family = "mono", size = 9))

Performance

For players having their first Tommy John surgery following their MLB debut, it could be interesting to understand how their performance compares before versus after their surgery. In the following graph, we can notice that players average more innings following their surgery if it is done while they are younger. On the other hand, as players age into their thirties they almost always average fewer innings after the surgery. Players that did not play following their surgery are shown in red, but this trend holds even just looking at the blue dots, which represent players that play again after their surgery.

diffIP <- seasonal_complete %>%
  filter(surgery1time == "after debut") %>%
  group_by(Player, statsTime) %>%
  summarize(avg_IP = round(mean(IP),2), 
            surgery_age = time_length(difftime(`Surgery 1`, birth_date), "years")) %>%
  select(Player, surgery_age, statsTime, avg_IP) %>%
  distinct() %>%
  pivot_wider(id_cols = Player:surgery_age, names_from = statsTime, values_from = avg_IP) 

diffIP$after <- replace_na(diffIP$after, 0)

diffIP %>%
  mutate(difference = after - before, 
         leaveMLB = case_when(after == 0 ~ "leave MLB", 
                              TRUE ~ "continue playing")) %>%
  ggplot(aes(x=surgery_age, y= difference, color = leaveMLB))+
  geom_point()+
  theme_classic()+
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") %>%
  labs(title = "Pitchers tend to average fewer seasonal innings post Tommy John surgery than prior to it", subtitle = "Players that leave the MLB completely are shown in<strong><span style='color:red'> red</span></strong></b>", x = "surgery age", y = "Average seasonal IP after - before surgery")+
  theme(plot.title = element_markdown(size = 9, family = "mono"),
        plot.subtitle = element_markdown(size = 9, family = "mono"),
        plot.title.position = "plot", 
        axis.title.x = element_markdown(size = 8, family = "mono"), 
        axis.title.y = element_markdown(size = 8, family = "mono"),
        legend.position = "none")+
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  scale_color_manual(values = c("navy", "red"))

Exploratory Data Analysis

If a team has a player decide to undergo Tommy John surgery, or is debating signing a player needing surgery, it would be helpful to understand what might be correlated with a player returning to the majors following their surgery.

A few possible performance based variables to look into could be the pitcher’s average strikeout rate and average WHIP, as well as total innings pitched prior to the surgery. We see that pitchers that return to the MLB tend to have slightly lower (better) WHIPs and slightly higher (better) strikeout rates. Additionally, the innings pitched variable appears highly informative. 75.29% of the pitchers that did not return to the MLB pitched fewer than 200 innings prior to their surgery. In comparison, only 51.5% of the pitchers that played after pitched less than 200 innings.

#understand relationship between 200 innings pitched prior to surgery and return to MLB rate
stats %>%
  mutate(manyIP = case_when(totalIP >= 200 ~ "high", 
                            TRUE ~ "low")) %>%
  group_by(manyIP, played_after) %>%
  count()

#create graph for relationship between average WHIP and playing after surgery
stats %>%
  filter(avgWHIP < 2.5, avgWHIP > 0.5) %>%
  ggplot(aes(x=avgWHIP, fill = played_after))+
  geom_density(alpha=0.5) +
  theme_classic()+
  labs(x= "average WHIP", y = "", title = "Average WHIP prior to Tommy John surgery comparing pitchers", subtitle =  "that <strong><span style='color:navy'>play again</span></strong></b> versus <strong><span style='color:red'>don't</span></strong></b>")+
  scale_fill_manual(values = c("navy", "red"))+
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(), 
        plot.title.position = "plot", 
        plot.title = element_markdown(family = "mono", size =14),
        plot.subtitle = element_markdown(family = "mono", size =14),
        legend.position = "none", 
        axis.title.x = element_markdown(family="mono", size=11))

#create graph for relationship between average strikeout rate and playing after surgery
stats %>%
  ggplot(aes(x=kRate, fill = played_after))+
  geom_density(alpha=0.5) +
  theme_classic()+
  labs(x= "strikeout rate per inning", y = "", title = "Average K rate prior to Tommy John surgery comparing pitchers",  subtitle = "that <strong><span style='color:navy'>play again</span></strong></b> versus <strong><span style='color:red'>don't</span></strong></b>")+
  scale_fill_manual(values = c("navy", "red"))+
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(), 
        plot.title.position = "plot", 
        plot.title = element_markdown(family = "mono", size =14), 
        plot.subtitle = element_markdown(family = "mono", size =14), 
        legend.position = "none", 
        axis.title.x = element_markdown(family="mono", size=11))

#create graph for relationship between total innings pitched and playing after surgery
stats %>%
  ggplot(aes(x=totalIP, fill = played_after))+
  geom_density(alpha=0.5) +
  theme_classic()+
  labs(x= "total IP", y = "", title = "Total IP prior to Tommy John surgery comparing players that <strong><span style='color:navy'>play again</span></strong></b> versus <strong><span style='color:red'>don't</span></strong></b>")+
  scale_fill_manual(values = c("navy", "red"))+
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(), 
        plot.title.position = "plot", 
        plot.title = element_markdown(family = "mono", size =10), 
        legend.position = "none", 
        axis.title.x = element_markdown(family="mono", size=9))

In addition to performance metrics, we can see if there is any relationship between a pitcher playing after their surgery and their role (starter or reliever) as well as whether they were a high school draftee.

# create mosiac plot showing relationship between starter versus reliever and results of playing again
stats %>%
  mutate(starter = recode(starter, "0" = "reliever", "1" = "starter")) %>%
  ggplot()+
  geom_mosaic(aes(x = product(played_after,starter), fill = played_after), alpha = 0.7)  +
  labs(x= "", y = "", title = "Starters <strong><span style='color:navy'>play after surgery</span></strong></b> more often than relievers")+
  scale_fill_manual(values = c("navy", "red")) + 
  theme_classic()+
  theme(axis.text.y = element_blank(), 
        axis.ticks.y = element_blank(), 
        plot.title = element_markdown(family = "mono", size =13), 
        legend.position = "none", 
        axis.text.x = element_markdown(family = "mono", size = 11))

# create mosiac plot showing relationship between high school draftee indicator and results of playing again
stats %>%
  mutate(hs_draftee = case_when(hs_draftee == TRUE ~ "high school", 
                                hs_draftee == FALSE ~ "college/international")) %>%
  ggplot()+
  geom_mosaic(aes(x = product(played_after,hs_draftee), fill = played_after), alpha = 0.7)  +
  labs(x= "", y = "", title = "High school draftees <strong><span style='color:navy'>play after surgery</span></strong></b> slightly more often", subtitle = "than college or international picks")+
  scale_fill_manual(values = c("navy", "red")) + 
  theme_classic()+
  theme(axis.text.y = element_blank(), 
        axis.ticks.y = element_blank(), 
        plot.title = element_markdown(family = "mono", size =13), 
        plot.subtitle = element_markdown(family = "mono", size =13),
        legend.position = "none", 
        axis.text.x = element_markdown(family = "mono", size = 11))

A final variable to investigate is the age at which the player is having the Tommy John surgery. Below, we see that pitchers that pitch again tend to be slightly younger in age at the time of their surgery in comparison to those that do not play again. However, this difference is not quite as substantial as one might have originally thought.

stats %>%
  ggplot(aes(x=surgery_age, fill = played_after))+
  geom_density(alpha=0.5) +
  theme_classic()+
  labs(x= "surgery age", y = "", title = "Surgery age comparing players that <strong><span style='color:navy'>play again</span></strong></b> versus <strong><span style='color:red'>don't</span></strong></b>")+
  scale_fill_manual(values = c("navy", "red"))+
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(), 
        plot.title.position = "plot", 
        plot.title = element_markdown(family = "mono", size =12), 
        legend.position = "none", 
        axis.title.x = element_markdown(family="mono", size = 11))

Model Introduction and Priors

Having a basic understanding of some of these variables, we can now build a logistic regression model. In building this model, we will let Y be a binary indicator of if the pitcher plays again at the MLB level following his Tommy John surgery, which occurs at probability π. Understanding the results of the exploratory data analysis, this model will utilize the following predictors:

• X1: indicator variable for if a player is a starter (1) or reliever (0)
  
• X2: continuous variable representing the total number of innings pitched by a player prior to their surgery
  
• X3: continuous variable representing the player’s average WHIP prior to their surgery
  
• X4: indicator variable of if a player was a high school draftee
  
• X5: continuous variable representing a pitcher’s strikeout rate prior to their surgery
  
• X6: continuous variable for the age at which a player got surgery
  

Starting with the centered intercept β0c, we recall our prior understanding that for the average player there is roughly a 73% chance of a player returning to the MLB, i.e. π≈0.73. As a result, we can set the prior mean for β0c on the log(odds) scales to 1:

Setting a range for this normal prior, we have a vague understanding of the log(odds of returning) ranges from roughly 0.3 to 1.7 (1±2∗0.35). In other words, the odds of a typical pitcher playing after at the MLB level following surgery could be somewhere between:

We have only a very basic understanding of each predictor’s relationship with a player’s likelihood of returning to pitch at the MLB level and thus will utilize weakly informative priors for these coefficients. As a result, their coefficients are centered at 0 and they have a fairly wide standard deviation. As a whole, the model can be summarized as:

stats_clean <- stats %>%
  mutate(played_after = relevel(played_after, ref='0'), 
         starter = relevel(starter, ref='0')) %>%
  filter(avgWHIP < 3.5)
played_after_mod <- stan_glm(played_after ~ starter + totalIP + avgWHIP + hs_draftee + kRate + surgery_age,
                             data = stats_clean, family = binomial,
                             prior_intercept = normal(1, 0.35),
                             prior = normal(0, 2.5, autoscale = TRUE),
                             chains = 4, iter = 5000*2, seed = 84735, refresh=0)

prior_summary(played_after_mod)$prior_intercept
prior_summary(played_after_mod)$prior

This model uses 315 out of the 318 pitchers that had Tommy John surgery following their major league debut. Three pitchers - Joe Mantiply, Jack Egbert, and Erik Bedard - were removed as their WHIPs were large outliers (all over 4) due to each pitching fewer than three MLB innings before their Tommy John surgery. With these three included in the model, the 80% credible interval on WHIP was almost entirely positive, indicating that as pitchers perform worse their likelihood of returning to the MLB after surgery rises. While Mantiply, Egbert, and Bedard interestingly enough all did return to pitch at the MLB level following their Tommy John surgery, the overall trend that we saw in the density plots indicates that better performing pitchers return more after surgery and thus our model should capture that.

# check out the players that had outlier WHIPs
stats %>%
  filter(avgWHIP > 3.5)

Check the diagnostics

Before getting into the results of the model, we should check the diagnostics of the simulation. The model above was run using Markov Chain Monte Carlo simulation techniques. This process enables us to estimate the posterior of our classification model and get a better understanding of what might contribute to a pitcher returning to the MLB following their surgery. However, not all simulations are perfect. We are looking for Markov chains that roam the sample space of π and in the end imitate a random sample that converges to the posterior. To assess if our simulation does this, we can look at trace plots, compare parallel chains, and assess autocorrelation.

The set of plots below on the left shows the trace plots for our simulation. These plots have no discernible trends or patterns, which is what we want - this implies that our Markov chains are stable. The set of plots below on the right compare parallel chains. In these plots we are looking for consistency across the four chains, which we definitely have. This provides evidence that our simulation is stable and long enough. Additionally, the autocorrelation plots show our Markov chain sufficiently mimics the behavior of an independent sample (which is what we want).

mcmc_trace(played_after_mod)
mcmc_dens_overlay(played_after_mod)+
  ylab("density")

mcmc_acf(played_after_mod)

Model results

Now that our simulation has proven itself reliable, we can check out the coefficients of our model. Using an 80% confidence interval, we see that a pitcher’s total innings pitched, strikeout rate, and surgery age are significant variables in predicting whether they will pitch again at the MLB level (when controlling for the other predictors in the model). On the other hand, the high school draftee indicator variable and average WHIP are not significant variables. Whether or not a player is a starter or reliever is borderline significant as the lowest end of the confidence interval is just slightly negative, indicating it is reasonable to assume that starters are more likely to return to play the MLB level holding all other variables constant.

tidy(played_after_mod, effects = "fixed", conf.int = TRUE, conf.level = 0.80) %>%
  mutate(estimate = round(estimate, 6), 
         std.error = round(std.error, 4), 
         conf.low = round(conf.low, 6), 
         conf.high = round(conf.high, 6)) %>%
  gt() %>%
  tab_header(
    title = md("**Model results**")
  ) 

Looking at the significant coefficients, we understand that:

• For every additional 100 total innings a player has pitched prior to his surgery, there is an 80% posterior chance that the odds of him returning to the pitch at the MLB level following the surgery increase by somewhere between 17.8% and 43% holding all other variables constant: (e0.001635∗100,e0.003577∗100)=(1.178,1.43).
  
• For every additional strikeout per nine inning pitched a player has prior to his surgery, there’s an 80% posterior chance the odds of him returning to pitch at the MLB level following the surgery increase by somewhere between 7.7% and 28.6% holding all other variables constant: (e0.665150∗19,e2.260704∗19)=(1.077,1.286).
  
• For every year older that a player gets when he has surgery, we estimate there is an 80% posterior chance that the odds of him returning to pitch at the MLB level are decreased by a factor of 75.68% to 88.84% holding all other variables constant: (e−0.278644,e−0.118299)=(0.7568,0.8884).

Using the model

With a model, we can now use it in context. For example, suppose we wanted to predict whether or not Tyler Glasnow will return to play again at the MLB level following his Tommy John surgery that occurred on August 4, 2021. Glasnow is currently a starter for the Tampa Bay Rays. Up until his surgery, he has averaged a WHIP of 1.248, 1.255 strikeouts per inning, and has pitched 403 innings. He was drafted out of high school and had surgery at age 27.948. To predict whether Glasnow will return to the MLB, we can approximate the posterior predictive model for the binary outcomes of Y, whether the pitcher plays again, where

We can simulate 20,000 posterior predictions of Y and learn 91.085% of them call for Glasnow to return to the MLB, making it reasonable to predict he will return after recovering from his surgery.

# Posterior predictions of binary outcome
set.seed(84735)
binary_prediction <- posterior_predict(
  played_after_mod, newdata = data.frame(starter="1", totalIP = 403, avgWHIP=1.248, hs_draftee = TRUE, kRate=1.255, surgery_age = 27.948))

table(binary_prediction)
colMeans(binary_prediction)

Model evaluation

To evaluate our classification model, we will ask three questions.

First, how fair is the model?
While some of this information seems pretty private (e.g. when a pitcher has surgery and his age), this is all public data available and players know that their information being public is part of what comes with being a MLB player. Thus, while we cannot ask all players whether they are okay with their information being used in this model, we will assume it okay.

Second, how wrong is the model?
We can complete a posterior predictive check to justify that data simulated from our posterior logistic regression model has characteristics similar to the original data. If this is the case, we can confirm that the assumptions behind our Bayesian logistic regression model (normal priors are appropriate and there is independence between priors) are reasonable. To complete this posterior predictive check, we record the proportion of outcomes Y that are 1 (the proportion of players that play again) from each of 100 posterior simulated datasets. The histogram below of these simulated “play after surgery” rates confirms that they are consistent with the original data.

proportion_again <- function(x){mean(x == 1)}
pp_check(played_after_mod, nreps = 100,
         plotfun = "stat", stat = "proportion_again") + 
  xlab("probability of playing after surgery")

Third and finally, how accurate are our posterior classifications?
In the classification setting where we are predicting whether a pitcher returns to pitch again at the MLB following his surgery, we are interested if our classifications of Y are right or wrong, specifically how often we are right. Thinking about the perspective of a scout of a team, we would probably prefer a model that is highly accurate at identifying players that will not return to play again at the major league level after surgery, as signing these players could lead to large losses as the team gets to return out of them. In that case, we would want a model that favors specificity over sensitivity, meaning we more accurately predict the result correctly for players that do not play again than for players that do. However, it is important to recognize making this sacrifice will hurt our overall model accuracy a bit.

Taking this all into account, our final model will predict a pitcher will return to the MLB after their surgery only if their posterior predicted probability is greater than 72%. Evaluating our model with this classification cutoff on cross-validated data, we see that our model is overall 64% accurate with a specificity of 68.7% and a sensitivity of 63.56%. Unfortunately, this overall accuracy is lower than the no-information rate of our data. The no-information rate is the accuracy we would get if we simply always predicted the majority outcome. In this case, that would mean predicting that a pitcher always returns to play again in the MLB and would give us an accuracy of 73%. The negative of this is our sensitivity (the proportion of players who do not return to the MLB that we predicted correctly) would be 0. So white our model is not as accurate as we might hope, we will go with it and hopefully make further improvements in the future.

set.seed(84735)
cv_accuracy <- classification_summary_cv(
  model = played_after_mod, data = stats_clean, cutoff = 0.72, k = 10)
cv_accuracy$cv

Exploratory Data Analysis

We will first start by looking at how slugging percent against (will be abbreviated SLG from now on) varies with age for starting pitchers. As seen in the plot below, there appears to be a very weak relationship between pitcher age and SLG, with SLG slightly increasing as pitcher age increases. However, this plot has a few issues. First, it pools all seasons together and ignores the data’s grouped structure: an individual pitcher’s SLG one season will be have some correlation to their SLG the following seasons. Second, there is some kind of selection bias happening with age. The few pitchers that are pitching in their late 30s and early 40s are likely doing so because they’ve had a very successful career!

starters19 %>%
  ggplot(aes(x=Age, y=SLG)) +
  geom_point()+
  geom_smooth(se = FALSE)+
  theme_classic()+
  labs(x = "pitcher age", y="SLG", title = "Pitcher age versus slugging percerntage")+
  theme(plot.title.position = "plot", 
        plot.title = element_markdown(family = "mono", size =12), 
        legend.position = "none", 
        axis.title.x = element_markdown(family="mono", size = 11), 
        axis.title.y = element_markdown(family="mono", size = 11))

To minimize these issues, we will look at the data in its grouped format. We can start by taking a sample of 80 starters and plotting the smoothed linear trend in their SLG over the course of their career. Looking at the plot, it appears that for the majority of players SLG tends to increase over the course of their career, indicating their performance is worsening with age.

slopes <-starters19 %>% 
    group_by(Player) %>% 
    do({
      mod = lm(SLG ~ Age, data = .)
      data.frame(Intercept = coef(mod)[1],
                 Slope = coef(mod)[2])
    })

starters19 %>%
  left_join(slopes, by = "Player") %>%
  arrange(Player) %>%
  group_by(Player) %>%
  mutate(group_id = cur_group_id(), 
         positive_slope = as.factor(case_when(Slope >0 ~ 1, 
                                    TRUE ~ 0))) %>%
  filter(group_id <= 80) %>%
  ggplot(aes(x=Age, y=SLG, group = Player))+
  geom_smooth(aes(color = positive_slope), method = "lm", se = FALSE, size = 0.5)+
  theme_classic()+
  labs(title = "SLG seems to <strong><span style='color:red'>worsen</span></strong></b> with age more often than <strong><span style='color:navy'>improve</span></strong></b>")+
  theme(plot.title.position = "plot", 
        plot.title = element_markdown(family = "mono", size =12), 
        legend.position = "none", 
        axis.title.x = element_markdown(family="mono", size = 11), 
        axis.title.y = element_markdown(family="mono", size = 11))+
  scale_color_manual(values = c("red", "navy"))

With some very basic information about SLG trajectory overtime, we will now look into the impact on SLG of undergoing Tommy John surgery at a point in a player’s career. For the 515 players we are analyzing in this study, 85 of them underwent their first Tommy John surgery at some point after their debut. For these players, we will create an categorical variable seasonTime for each season they played. This variable will take one of three values.
- “no surgery”: any season taking place prior to their surgery
- “short after”: any season taking place within 3 years after the year of their surgery
- “long after”: any season taking place at least 3 years after the year of their surgery
For players that did not undergo surgery, all of their seasons will be “no surgery” seasons.

# number of players that underwent Tommy John surgery after their debut
starters19 %>%
  filter(surgery1time == "after debut") %>%
  group_by(Player) %>%
  count() %>%
  nrow()

# create the seasonTime variable
starters19 <- starters19 %>%
  mutate(seasonTime = case_when(surgery1time == "after debut" & (Year - surgery1Year) <= 3 & (Year - surgery1Year) > 0 ~ "short after",
                               surgery1time == "after debut" & (Year - surgery1Year) > 3 ~ "long after",
                               TRUE ~ "no surgery")) 

For the 85 starters that underwent surgery, the average number of seasons prior the surgery was 5.15. This distribution of number of seasons before surgery can be seen in the density plot below. Of these 85 starters, 25 did not return to play again the MLB level following their surgery, and only 35 played “long after” their surgery. 33 players pitched in all three time periods.

breakdown <- starters19 %>%
  filter(surgery1time == "after debut") %>%
  group_by(Player, seasonTime) %>%
  count() %>%
  pivot_wider(id_cols = Player, names_from = seasonTime, values_from = n) %>%
  mutate(`short after` = coalesce(`short after`, 0),
         `long after` = coalesce(`long after`, 0)) %>%
  ungroup() 

# number not playing after surgery
breakdown %>%
  mutate(after = `short after` + `long after`) %>%
  group_by(`after`) %>%
  count() %>%
  head(1)

# number not playing long after surgery
breakdown %>%
  group_by(`long after`) %>%
  count() %>%
  head(1)

# density plot of seasons prior to surgery
breakdown %>%
  ggplot(aes(x=`no surgery`))+
  geom_density(fill = "lightcyan3")+
  labs(x= "seasons prior to surgery", title = "Number of seasons prior to surgery for 85 MLB starters")+
  theme_classic()+
  theme(plot.title.position = "plot", 
        plot.title = element_markdown(family = "mono", size = 12), 
        axis.title.x = element_markdown(family="mono", size = 11), 
        axis.title.y = element_markdown(family="mono", size = 11))

For the 33 starters that played in all three time periods, we plot the trend in SLG across the three periods below, grouped by player. 30 players are included that pitched at least 10 innings in each of the three time periods. For most players, SLG appears to increase slightly as they move into the period of “long after” their surgery, indicating poorer performance. We see two major outliers in SLG prior to having surgery. Matt Riley is on the high end with a SLG of 0.756, though this number is coming from only pitching 11 innings prior to his surgery. On the other hand, Jose Fernandez is on the low end with a phenomenal SLG against of 0.272 over 223.4 innings pitched in two seasons.

breakdown %>%
  filter(`long after` > 0 & `short after` >0) %>%
  left_join(seasonal_complete_starters, by = "Player") %>%
  select(-c(`long after`, `no surgery`, `short after`)) %>%
   mutate(seasonTime = case_when(surgery1time == "after debut" & (Year - surgery1Year) <= 3 & (Year - surgery1Year) > 0 ~ 2,
                               surgery1time == "after debut" & (Year - surgery1Year) > 3 ~ 3,
                               TRUE ~ 1)) %>%
  group_by(Player, seasonTime) %>%
  summarize(avgSLG = weighted.mean(SLG, w=IP), totalIP = sum(IP)) %>%
  filter(Player %!in% c("Erik Bedard", "Joe Wieland", "Shawn Hill")) %>% # did not pitch at least 10 innings per period 
  ggplot(aes(x=seasonTime, y=avgSLG, group = Player))+
  geom_point()+
  geom_smooth(color = "navy")+
  labs(title = "SLG against across for 30 starters over career", y = "SLG against")+
  scale_x_continuous(breaks = 1:3, labels = c("no surgery", "short after", "long after"))+
  theme_classic()+
  theme(plot.title.position = "plot", 
        plot.title = element_markdown(family = "mono", size = 12), 
        axis.title.x = element_markdown(family="mono", size = 11), 
        axis.title.y = element_markdown(family="mono", size = 11))

Model Intro and Priors

The goal with this model will be to predict a player’s seasonal SLG as a function of their age and the seasonTime variable. Again, this variable will remain “no surgery” for all seasons for players that never had Tommy John surgery. This model will use data described earlier that contains information on 481 starters, but filters out seasons with less than 15 IP due to the abnormal SLGs that only few innings pitched often give. Originally, I had thought of including IP as a predictor to account for these abnormal SLGs, but the goal of this analysis is to help scouts predict a pitcher’s success prior to when the season starts - thus, only including variables that are known prior to the season starting is necessary.

The following predictors will be used in this model:
- X1: Pitcher age at season
- X2: Pitcher age squared at season
- X3: Indicator for if the season is happening “short after” surgery (reference level “no surgery”)
- X4: Indicator for if the season in happening “long after” surgery (reference level “no surgery”)

The model built will take the form of a hierarchical model with varying intercepts, allowing the intercept (baseline SLG) to vary from pitcher to pitcher. I really contemplated a model with varying slopes in addition to intercepts, which would allow the rate at which SLG varies with age to change from pitcher to pitcher, but ultimately the runtime and complexity of this model swayed me to stick with just varying intercepts.

#get numbers for prior on intercept
starters19 %>% 
  ungroup() %>%
  filter(IP > 15) %>%
  summarize(mean(SLG), sd(SLG))

#get data ready - this model does not use all players as 
slgData <- starters19 %>%
  ungroup() %>%
  filter(IP > 15)  %>%
  mutate(age2 = Age^2, 
         seasonTime = case_when(surgery1time == "after debut" & (Year - surgery1Year) <= 3 & (Year - surgery1Year) > 0 ~ 1,
                               surgery1time == "after debut" & (Year - surgery1Year) > 3 ~ 2,
                               TRUE ~ 0), 
         seasonTime = as.factor(seasonTime), 
         seasonTime = relevel(seasonTime, ref='0')) 

slgMod<- stan_glmer(
  SLG ~ Age + I(Age^2) + seasonTime + (1 | Player),
  data = slgData, family = gaussian,
  prior_intercept = normal(0.44, 0.069),
  prior = normal(0, 2.5, autoscale = TRUE),
  prior_aux = exponential(1, autoscale = TRUE),
  prior_covariance = decov(reg = 1, conc = 1, shape = 1, scale = 1),
  chains = 4, iter = 5000*2, seed = 84735, adapt_delta = 0.99999
)
#write_rds(slgMod, "data/slgMod.rds")

slgMod <- read_rds("data/slgMod.rds")

# Check the prior specifications
prior_summary(slgMod)$prior_intercept
prior_summary(slgMod)$prior

Checking the diagnostics

Just like with the logistic regression model, we must check the diagnostics of the simulation. Plotting the trace plots for the global parameters below on the left, we see no discernible trends or patterns, which is what we want, implying that our Markov chains are stable. Additionally, the autocorrelation plots shown below on the right shows that our Markov chain sufficiently mimics the behavior of an independent sample. We also learn from a test of R-hat that there is evidence of parallel chains, as variability between chains is roughly equal to that within chains for all parameters.

columns <- as.data.frame(slgMod)[,c(1:5, 445:446)]
mcmc_trace(columns)
mcmc_acf(columns)

rhat(slgMod) %>%
  head(5)

Posterior results

tidyMod <- tidy(slgMod, effects = "fixed",
                       conf.int = TRUE, conf.level = 0.80) 
tidyMod %>%
  mutate(estimate = round(estimate, 6), 
         std.error = round(std.error, 4), 
         conf.low = round(conf.low, 6), 
         conf.high = round(conf.high, 6)) %>%
  gt() %>%
  tab_header(
    title = md("**Model results**")
  ) 

B0 <- tidyMod$estimate[1]
B1 <- tidyMod$estimate[2]
B2 <- tidyMod$estimate[3]

post_med <- function(x){B0 + B1 * x + B2 * x^2}

slgData %>%
  add_fitted_draws(slgMod, n = 200, re_formula = NA) %>%
  ggplot(aes(x = Age, y = SLG)) +
    geom_line(aes(y = .value, group = .draw), alpha = 0.1) +
    geom_function(fun = post_med, color = "blue")+
  labs(title = "Possible posterior relationship between Age and SLG")+
  theme_classic()+
  theme(plot.title.position = "plot", 
        plot.title = element_markdown(family = "mono", size = 12), 
        axis.title.x = element_markdown(family="mono", size = 11), 
        axis.title.y = element_markdown(family="mono", size = 11), 
        legend.title = element_markdown(family="mono", size = 11), 
        legend.text = element_markdown(family="mono", size = 11))+
  xlim(19,36)

mcmc_areas(slgMod, pars = vars(starts_with("season")), prob = 0.8) + 
  geom_vline(xintercept = 0)+
  scale_y_discrete(labels = c("short after surgery", "long after surgery"))

slgData %>%
  group_by(Player) %>%
  mutate(nseasons = n()) %>%
  filter(nseasons == 17) %>%
  ggplot(aes(x=Age, y=SLG, color = Player))+
  geom_point()+
  scale_color_manual(values = c("navy", "red", "darkgreen", "goldenrod"))+
  labs(title = "Variation in SLG from season to season is high")+
  theme_classic()+
  theme(plot.title.position = "plot", 
        plot.title = element_markdown(family = "mono", size = 12), 
        axis.title.x = element_markdown(family="mono", size = 11), 
        axis.title.y = element_markdown(family="mono", size = 11), 
        legend.title = element_markdown(family="mono", size = 11), 
        legend.text = element_markdown(family="mono", size = 11))

Group specific relationships

With this model, we can compare posterior results between pitchers. For example, we will look at the posterior results of 2019 Hall of Fame inductee Ryan Dempster and former starter Jason Hammel. The summary below compares the typical SLGs of these pitchers to the average pitcher, b0j=β0j−β0. When controlling for age and surgery status, Dempster’s seasonal SLG tends to be significantly lower (better) than the average pitchers’s, whereas Hammel’s tends to lean slightly better than average but not significantly. For example, there’s an 80% posterior chance that Dempster’s typical seasonal SLG is between 0.057 and 0.02 points better than average. This is not surprising given that Dempster was admitted to the Hall of Fame!

tidy(slgMod, effects = "ran_vals",
     conf.int = TRUE, conf.level = 0.80) %>% 
  filter(level %in% c("Ryan_Dempster", "Jason_Hammel")) %>% 
  select(level, estimate, conf.low, conf.high) %>%
  rename(Player = level) %>%
  gt() %>%
  tab_header(
    title = md("**Posterior results for Dempster and Hammel**")
  )

If Dempster and Hammel were to come back and play in the 2022 season, we could predict their next season’s SLG using their age and surgery status. More importantly, we can predict SLG for pitchers that are currently active and that we do not have prior MLB seasonal statistics on. For example, we can try to predict the 2022 seasonal SLG for MLB 2021 first round draft pick Jack Leiter. Given that we don’t have any information on prior seasonal SLGs, his posterior predicted interval is much wider. In a future model, a variable that we could add would be draft round, which might help explain more of the baseline variation in SLG.

# Simulate posterior predictive models for the 3 pitchers
set.seed(84735)
predict_next_season <- posterior_predict(
  slgMod,
  newdata = data.frame(
    Player = c("Ryan Dempster", "Jason Hammel", "Jack Leiter"), 
    Age = c(44, 39, 21), age2 = c(1936,1521,21), seasonTime = c("2", "0", "0")))

# Posterior predictive model plots
mcmc_areas(predict_next_season, prob = 0.8) +
 ggplot2::scale_y_discrete(
   labels = c("Ryan Dempster", "Jason Hammel", "Jack Leiter"))

tidy_sigma <- tidy(slgMod, effects = "ran_pars")
tidy_sigma

sigma_0 <- tidy_sigma[1,3]
sigma_y <- tidy_sigma[2,3]
sigma_0^2 / (sigma_0^2 + sigma_y^2)
sigma_y^2 / (sigma_0^2 + sigma_y^2)

Model evaluation

Just as with our logistic regression model, to evaluate this model we will ask three questions.

First, how fair is the model?
The answer to this question is pretty much the same as with the other model - this model uses what some may consider private information, but these players understand that this information being made public is part of what comes with being an MLB player. Thus, while we cannot ask all players whether they are okay with their information being used in this model, we will assume it okay.

Second, how wrong is the model? The posterior predictive check below shows that our model produces posterior simulated datasets of seasonal SLG that are fairly consistent with the main features in the original data the model was run with. This is good news, and allows us to move forward to question 3.

pp_check(slgMod)+
   xlab("SLG")

Finally, what’s the predictive accuracy of this model?

We answer this question by seeing how well our model did predicting the seasonal SLGs for the pitchers that were part of our data. The observed seasonal SLGs tend to be 0.033 points, or 0.539 standard deviations, from their posterior mean predictions. Knowing the randomness and high variability of SLG, this amount of error does not seem unreasonable. The model also has a good proportion of observed SLGs (60.3%) that fall within their respective 50% posterior predictive interval.

set.seed(84735)
prediction_summary(model = slgMod, data = slgData)

Unfortunately, I was unable to obtain cross-validated metrics of posterior predictive accuracy due to the model’s complexity and large amount of data. These metrics would help us understand how well our model predicts the seasonal SLGs of players that weren’t included in our sample data. Despite this unfortunate circumstance, the not cross-validated model results proved to be encouraging.