BABIP

1. Randomly split at bats for each pitcher

I’ll split each ball in play for each pitcher into two random groups. Then we’ll compare the BABIP in each group.

d2 %>% mutate(group12=sample(1:2,n(), T)) %>% group_by(pitcher_name, group12) %>% summarize(ipo=sum(event %in% inplayouts), iph=sum(event %in% inplayhits), N=iph+ipo, BABIP=iph/(iph+ipo))  %>% filter(!is.na(pitcher_name)) %>% group_by(pitcher_name) %>% summarize(BABIP1=BABIP[1], BABIP2=BABIP[2], N1=N[1], N2=N[2], N=N1+N2) %>% filter(N1>80 & N2>80) %>% {ggplot(., aes(BABIP1, BABIP2)) + geom_point(aes(color=N1+N2)) + geom_smooth()}

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

There appears to be no correlation. The correlation coefficient is about 0.04.

d2 %>% mutate(group12=sample(1:2,n(), T)) %>% group_by(pitcher_name, group12) %>% summarize(ipo=sum(event %in% inplayouts), iph=sum(event %in% inplayhits), N=iph+ipo, BABIP=iph/(iph+ipo))  %>% filter(!is.na(pitcher_name)) %>% group_by(pitcher_name) %>% summarize(BABIP1=BABIP[1], BABIP2=BABIP[2], N1=N[1], N2=N[2], N=N1+N2) %>% filter(N1>80 & N2>80) %>% {cor(.$BABIP1, .$BABIP2)}

## [1] 0.03165568

It would be nice to see what this would look like comparing another quantity, such as percentage of balls put in play.

d2 %>% mutate(group12=sample(1:2,n(), T)) %>% group_by(pitcher_name, group12) %>% summarize(inplay=sum(event %in% c(inplayhits, inplayouts)), notinplay=sum(!(event %in% c(inplayouts, inplayhits))), N=inplay+notinplay, InPlayPer=inplay/(inplay+notinplay))  %>% filter(!is.na(pitcher_name)) %>% group_by(pitcher_name) %>% summarize(InPlayPer1=InPlayPer[1], InPlayPer2=InPlayPer[2], N1=N[1], N2=N[2], N=N1+N2) %>% filter(N1>80 & N2>80) %>% {ggplot(., aes(InPlayPer1, InPlayPer2)) + geom_point(aes(color=N1+N2)) + geom_smooth() + geom_text(aes(label=ifelse(N1+N2>850, pitcher_name, '')))}

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

We see from the play that there is decent correlation between the two groups. The correlation coefficient is .61, which is significant.

d2 %>% mutate(group12=sample(1:2,n(), T)) %>% group_by(pitcher_name, group12) %>% summarize(inplay=sum(event %in% c(inplayhits, inplayouts)), notinplay=sum(!(event %in% c(inplayouts, inplayhits))), N=inplay+notinplay, InPlayPer=inplay/(inplay+notinplay))  %>% filter(!is.na(pitcher_name)) %>% group_by(pitcher_name) %>% summarize(InPlayPer1=InPlayPer[1], InPlayPer2=InPlayPer[2], N1=N[1], N2=N[2], N=N1+N2) %>% filter(N1>80 & N2>80) %>% {cor(.$InPlayPer1, .$InPlayPer2)}

## [1] 0.6627211

2. Split up each pitchers results by the two halves of the season

I’ll do the same as before, but now assign groups so the first half of the pitchers in play balls are in one group, the rest in the other.

d2 %>% group_by(pitcher_name) %>% mutate(group12=rep(c(1,2), c(floor(n()/2),n()-floor(n()/2)))) %>% ungroup %>% group_by(pitcher_name, group12) %>% summarize(ipo=sum(event %in% inplayouts), iph=sum(event %in% inplayhits), N=iph+ipo, BABIP=iph/(iph+ipo))  %>% filter(!is.na(pitcher_name)) %>% group_by(pitcher_name) %>% summarize(BABIP1=BABIP[1], BABIP2=BABIP[2], N1=N[1], N2=N[2], N=N1+N2) %>% filter(N1>80 & N2>80) %>% {ggplot(., aes(BABIP1, BABIP2)) + geom_point(aes(color=N1+N2)) + geom_smooth()}

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Again, no clear pattern, and a correlation coefficient near zero.

d2 %>% group_by(pitcher_name) %>% mutate(group12=rep(c(1,2), c(floor(n()/2),n()-floor(n()/2)))) %>% ungroup %>% group_by(pitcher_name, group12) %>% summarize(ipo=sum(event %in% inplayouts), iph=sum(event %in% inplayhits), N=iph+ipo, BABIP=iph/(iph+ipo))  %>% filter(!is.na(pitcher_name)) %>% group_by(pitcher_name) %>% summarize(BABIP1=BABIP[1], BABIP2=BABIP[2], N1=N[1], N2=N[2], N=N1+N2) %>% filter(N1>80 & N2>80) %>% with(cor(BABIP1, BABIP2))

## [1] -0.02318874

Let’s check proportion of balls put in play as before, but with this new grouping method.

d2 %>% group_by(pitcher_name) %>% mutate(group12=rep(c(1,2), c(floor(n()/2),n()-floor(n()/2)))) %>% ungroup %>% group_by(pitcher_name, group12) %>% summarize(inplay=sum(event %in% c(inplayhits, inplayouts)), notinplay=sum(!(event %in% c(inplayouts, inplayhits))), N=inplay+notinplay, InPlayPer=inplay/(inplay+notinplay))  %>% filter(!is.na(pitcher_name)) %>% group_by(pitcher_name) %>% summarize(InPlayPer1=InPlayPer[1], InPlayPer2=InPlayPer[2], N1=N[1], N2=N[2], N=N1+N2) %>% filter(N1>80 & N2>80) %>% {ggplot(., aes(InPlayPer1, InPlayPer2)) + geom_point(aes(color=N1+N2)) + geom_smooth() + geom_text(aes(label=ifelse(N1+N2>850, pitcher_name, '')))}

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Again, there is a pattern with .59 correlation coefficient.

d2 %>%  group_by(pitcher_name) %>% mutate(group12=rep(c(1,2), c(floor(n()/2),n()-floor(n()/2)))) %>% ungroup %>% group_by(pitcher_name, group12) %>% summarize(inplay=sum(event %in% c(inplayhits, inplayouts)), notinplay=sum(!(event %in% c(inplayouts, inplayhits))), N=inplay+notinplay, InPlayPer=inplay/(inplay+notinplay))  %>% filter(!is.na(pitcher_name)) %>% group_by(pitcher_name) %>% summarize(InPlayPer1=InPlayPer[1], InPlayPer2=InPlayPer[2], N1=N[1], N2=N[2], N=N1+N2) %>% filter(N1>80 & N2>80) %>% {cor(.$InPlayPer1, .$InPlayPer2)}

## [1] 0.5909042

I’m going to skip 3 and 4 and go to 5. Those seem like too much work.

5. Logistic regression to predict if ball in play is hit or out

We’re going to need to decide which features to include in the model. I worked backwards and removed all the ones I don’t want, so it’s a mess of code below. Basically I want to keep the location, speed, and spin data.

d2 %>% select(-c(gameday_link.x, gameday_link.y, end_tfs_zulu, event_num.x, event_num.y, event_es, event_es, event2_es, date, event4, pitcher_name, batter_name, event3, event2, next_.y, inning.y, inning_side.y, inning.x, inning_side.x, away_team_runs, home_team_runs, url, des, des_es, atbat_des, atbat_des_es, p_throws, b_height, stand, start_tfs, start_tfs_zulu, pitcher, batter, play_guid, code, count, on_3b, on_2b, on_1b, num, next_.x, type_confidence, id, type, x, y, o, cc, mt, tfs, tfs_zulu, sz_top, sz_bot, x0, y0, z0, vx0, vy0, vz0, sv_id, b, s)) %>% filter(event %in% c(inplayhits, inplayouts)) %>%  mutate(iph=ifelse(event %in% inplayhits, 1, 0))  %>% with(iph) %>% table

## .
##     0     1 
## 83667 35595

# data in play
dip <- d2 %>% select(-c(gameday_link.x, gameday_link.y, end_tfs_zulu, event_num.x, event_num.y, event_es, event_es, event2_es, date, event4, pitcher_name, batter_name, event3, event2, next_.y, inning.y, inning_side.y, inning.x, inning_side.x, away_team_runs, home_team_runs, url, des, des_es, atbat_des, atbat_des_es, p_throws, b_height, stand, start_tfs, start_tfs_zulu, pitcher, batter, play_guid, code, count, on_3b, on_2b, on_1b, num, next_.x, type_confidence, id, type, x, y, o, cc, mt, tfs, tfs_zulu, sz_top, sz_bot, x0, y0, z0, vx0, vy0, vz0, sv_id, b, s)) %>% filter(event %in% c(inplayhits, inplayouts)) %>%  mutate(iph=ifelse(event %in% inplayhits, 1, 0)) %>% select(-event)

I’m going to randomly assign the training and testing.

train_inds <- sample(1:nrow(dip), floor(.7*nrow(dip)), T)
dip_train <- dip[train_inds,]
dip_test  <- dip[-train_inds,]

First I’ll fit a simpler logistc regression model that only uses three of the factors.

mod1 <- glm(iph ~ start_speed + spin_rate + break_length, "binomial", dip_train)
pred1train <- plogis(predict(mod1, dip_train))
qplot(pred1train, dip_train$iph+rnorm(nrow(dip_train), 0, .04), alpha=.0001)

## Warning: Removed 782 rows containing missing values (geom_point).

pred1test <- plogis(predict(mod1, dip_test))
qplot(pred1test, dip_test$iph+rnorm(nrow(dip_test), 0, .02))

## Warning: Removed 586 rows containing missing values (geom_point).

dip_train %>% mutate(p1 = pred1train, iphf=as.factor(iph)) %>% {ggplot(data=., mapping=aes(x=p1, group=iphf, fill=iphf)) + geom_histogram(alpha=.6)}

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 782 rows containing non-finite values (stat_bin).

dip_test %>% mutate(p1 = pred1test, iphf=as.factor(iph)) %>% {ggplot(data=., mapping=aes(x=p1, group=iphf, fill=iphf)) + geom_histogram(alpha=.6)}

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 586 rows containing non-finite values (stat_bin).

The plots all show that the model is not able to find any difference between the two groups. The first and third plots shows that it can’t even find a difference on the training data, which should be easier than testing data. Let’s try another model with more features.

mod2 <- glm(iph ~ start_speed + spin_rate + break_length + pfx_x + pfx_z + px + pz + break_y + nasty, "binomial", dip_train)
pred1train <- plogis(predict(mod2, dip_train))
qplot(pred1train, dip_train$iph+rnorm(nrow(dip_train), 0, .04), alpha=.0001)

## Warning: Removed 782 rows containing missing values (geom_point).

pred1test <- plogis(predict(mod2, dip_test))
qplot(pred1test, dip_test$iph+rnorm(nrow(dip_test), 0, .02))

## Warning: Removed 586 rows containing missing values (geom_point).

dip_test %>% mutate(p1 = pred1test, iphf=as.factor(iph)) %>% {ggplot(data=., mapping=aes(x=p1, group=iphf, fill=iphf)) + geom_histogram(alpha=.6)}

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 586 rows containing non-finite values (stat_bin).

With these additional features, it still can’t really tell the difference between the two groups.

BABIP by zone

I’m curious to see how the zone the pitch is in affects the BABIP.

d2 %>% group_by(zone) %>% summarize(ipo=sum(event %in% inplayouts), iph=sum(event %in% inplayhits), N=iph+ipo, BABIP=iph/(iph+ipo)) %>% filter(N>100) %>% arrange(BABIP)

## # A tibble: 14 x 5
##     zone   ipo   iph     N BABIP
##    <dbl> <int> <int> <int> <dbl>
##  1    11  6024  2211  8235 0.268
##  2    12  5183  1984  7167 0.277
##  3    14  8367  3232 11599 0.279
##  4     1  4204  1627  5831 0.279
##  5     2  4878  1933  6811 0.284
##  6    13  7283  2935 10218 0.287
##  7     3  3716  1604  5320 0.302
##  8     9  6191  2711  8902 0.305
##  9     7  5537  2442  7979 0.306
## 10    NA   790   351  1141 0.308
## 11     4  7557  3386 10943 0.309
## 12     6  7558  3394 10952 0.310
## 13     5  9149  4231 13380 0.316
## 14     8  7230  3554 10784 0.330

zones <- data.frame(
  zone = (c(rep(1:9, each=4), rep(11:14, each=6))),
  x = c(1,2,2,1,2,3,3,2,3,4,4,3,1,2,2,1,2,3,3,2,3,4,4,3,1,2,2,1,2,3,3,2,3,4,4,3,  0,1,1,2.5,2.5,0, 2.5,4,4,5,5,2.5, 2.5,5,5,4,4,2.5, 0,2.5,2.5,1,1,0),
  y = c(3,3,4,4,3,3,4,4,3,3,4,4,2,2,3,3,2,2,3,3,2,2,3,3,1,1,2,2,1,1,2,2,1,1,2,2,  2.5,2.5,4,4,5,5, 4,4,2.5,2.5,5,5, 0,0,2.5,2.5,1,1, 0,0,1,1,2.5,2.5)
)
zones$bs <- ifelse(as.numeric(zones$zone) <=9, "S","B")

babipzonesdf <- d2 %>% group_by(zone) %>% summarize(ipo=sum(event %in% inplayouts), iph=sum(event %in% inplayhits), N=iph+ipo, BABIP=iph/(iph+ipo)) %>% filter(N>100) %>% arrange(BABIP) %>% inner_join(zones, by="zone")
ggplot(babipzonesdf) + geom_polygon(aes(x=x,y=y,fill=BABIP, group=zone), alpha=.8)  + scale_fill_gradientn(colours=c('blue','white','red'))

Here we can get something. The BABIP is highest in the lower two thirds of the strike zone, and lowest outside the strike zone. This isn’t too surprising since it’s hard to get contact on pitches outside the zone. But at the same time, this would mean that pitchers can control their BABIP by pitching to specific zones.

BABIP by nasty

Looks like nastiness is not a good predictor of BABIP by itself.

dip %>% {ggplot(data=., mapping=aes(x=nasty, group=iph, fill=iph)) + geom_histogram(alpha=.6)}

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 1141 rows containing non-finite values (stat_bin).

Conclusion

From some basic data exploration here, I’ve become more convinced in the assumption behind FIP/DIPS: that pitchers can’t control the proportion of balls in play that become hits. This means that saying that some pitchers are able to induce weak contact is not really a thing. I could do some exploration on this. The only thing useful I found was that BABIP is lower outside of the strike zone, which makes sense.