イントロ

ある打球がホームランになるかどうかは基本的に,

  • 打球の飛距離
  • 打球の方向
  • 球場の形状

によって決まると考えられる.

また打球の飛距離は,

  • 打球速度
  • 打球角度
  • 打球の回転
  • ボールの空気力学的性質
  • 気候 (主に気温)

などによって決まると考えられる.

ここでは, 2016-2017年の打球の速度と角度を調整した状態での本塁打確率を推定し, それを2018年の打球に適用することで打球の速度と角度以外の要因がどの程度変動しているかをぼんやりと調べる.

データは基本的にStatcast 2015-2018レギュラーシーズンを用いる.

# dataは18-10-02取得
require(mgcv)
require(knitr)
require(tidyverse)
df2 <- df2 %>% 
  mutate(HR_FL = ifelse(events == "home_run", 1, 0),
         month =  str_sub(game_date, 6, 7))

typeX <- df2 %>%
  filter(type == "X")%>%
  filter(! grepl("bunt", des)) # bunt除く

HR数/年度 or HR数/月

まずは前提として, MLB15-18の年度毎のHR数, 及びHR/wOBA分母 (概ねPA, IBB入ってないとか色々) を計算する.

df2 %>%
  group_by(game_year)%>%
  dplyr::summarise(HR =  sum(HR_FL, na.rm = TRUE),
                   wOBA_denom =  sum(woba_denom, na.rm = TRUE),
                   HR_rate = HR / wOBA_denom)%>%
  mutate_if(is.numeric, funs(round), 4)%>%
  kable()
game_year HR wOBA_denom HR_rate
2015 4910 181500 0.0271
2016 5611 182637 0.0307
2017 6105 183394 0.0333
2018 5585 183377 0.0305

18年の本塁打率は, 17年に比べれば低下したが, 16年とは同程度を維持している.

月ごとにプロットして見てみる. 3, 10月は試合数が少ないので除いた.

df2 %>%
  group_by(game_year, month)%>%
  dplyr::summarise(HR =  sum(HR_FL, na.rm = TRUE),
                   wOBA_denom =  sum(woba_denom, na.rm = TRUE),
                   HR_rate = HR / wOBA_denom)%>%
  filter(month != "10", month != "03")%>%
ggplot(aes(x = month,  y = HR_rate,
           colour = game_year, group =  game_year,
           linetype =  game_year)) +
  geom_line() + 
  geom_point(size = 1) +
  theme_bw(base_family = "HiraKakuPro-W3") +
  theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +
  labs(title = "月ごとのHR率.",
       subtitle = "MLB15-18.",
       x = "month", y = "HR/wOBA_denom"
  )

2018年は年間を通してだいたい2016年と同程度というところ. 基本的にはシーズン当初と終わりはHRが出にくいが, これはおそらく主に気温が低下するため. 2015年の不自然なシーズン後半での増加が目立つが, これはMLBボールの仕様の変化だろう.

各年度での打球の速度と角度を確認.

df2 %>%
  group_by(game_year)%>%
  dplyr::summarise(speed_mean = mean(launch_speed, na.rm = TRUE),
                   speed_sd = sd(launch_speed, na.rm = TRUE),
                   angle_mean = mean(launch_angle, na.rm = TRUE),
                   angle_sd =  sd(launch_angle, na.rm = TRUE))%>%
  mutate_if(is.numeric, funs(round), 1)%>%
  kable()
game_year speed_mean speed_sd angle_mean angle_sd
2015 86.5 14.4 11.1 27.1
2016 83.5 15.0 15.6 29.0
2017 81.5 16.4 15.9 29.1
2018 83.0 14.8 16.0 29.0

ここ数年で打球速度や角度には顕著な変化が見られる. これ自体も上で示したHR/wOBA_denomに影響があるはずだが, それ以外の要素も影響があったはずである.

速度と角度を調整した状態でのHRの出やすさの年度間比較

HRの出やすさを, 打球の速度+角度について調整し, それ以外 (ボールの空力的な性質, 気温 etc.) の要因の影響をざっくり見積もりたい.

各年度のデータにおいて, それぞれの打球の速度と角度を説明変数として, 本塁打率 (この場合, 打球の速度と角度が記録されている必要があるので分母は打球数BBE) を推定するモデルを一般化加法モデルGAMでつくる. それら各年度のモデルを利用して, 2018年のそれぞれの打球での推定HR確率を計算し, 推定HR数を集計する. これにより, 例えば2017年度と同じ条件であれば何本のHRであったか, などを推定可能となる. 有り体に言えばAlbertのパクリ. 2015年前半に関してはStatcastのデータは信頼性に疑問がある (リンク先2.8)ため, 以下では2016年以後について調べる.

# 下をコメントアウトして計算

# # モデルをたてる
# fit16 <- gam(factor(HR_FL) ~ s(launch_speed, launch_angle),
#            data = typeX %>%
#              filter(game_year ==  2016),
#            family = binomial)
# fit17 <- gam(factor(HR_FL) ~ s(launch_speed, launch_angle),
#              data = typeX %>%
#              filter(game_year ==  2017),
#              family = binomial)
# fit18 <- gam(factor(HR_FL) ~ s(launch_speed, launch_angle),
#              data = typeX %>%
#              filter(game_year ==  2018),
#              family = binomial)
# # 2018年のデータでHR確率を推定させる
# typeX.2018 <- typeX%>%
#   filter(game_year ==  2018)
# 
# typeX.2018 <- typeX.2018%>%
#   mutate(pred_16 =  predict(fit16,
#                              newdata = typeX.2018,
#                              type = "response"),
#          pred_17 =  predict(fit17,
#                              newdata = typeX.2018,
#                              type = "response"),
#          pred_18 =  predict(fit18,
#                              newdata = typeX.2018,
#                              type = "response"))
# # 適当な打球速度, 角度について推定させる
# # 図で示すのに使う
# new.data <- expand.grid(launch_angle = seq(10, 40, 0.1),
#                         launch_speed = seq(95, 110, 5))
# hr.prob <- new.data%>%
#   mutate(pred_16 =  predict(fit16,
#                              newdata = new.data,
#                              type = "response"),
#          pred_17 =  predict(fit17,
#                              newdata = new.data,
#                              type = "response"),
#          pred_18 =  predict(fit18,
#                              newdata = new.data,
#                              type = "response"))
# save(list = c("typeX.2018", "hr.prob"), file = "estimated_hr_2018.Rdata")

load("estimated_hr_2018.Rdata")

まず各年度のモデルから推定されたHR確率を図で示しておく.

hr.prob%>%
  gather(key = year, value = estimated_prob, pred_16, pred_17, pred_18)%>%
ggplot(aes(x =  launch_angle,
           y = estimated_prob,
           group = year,
           colour =  year,
           linetype = year)) +
  geom_line()+
  facet_wrap(~launch_speed)+
  theme_bw(base_family = "HiraKakuPro-W3")+
  labs(title = "打球速度+打球角度による推定HR確率の変化",
       subtitle =  "グラフの上の数字は打球速度 (mph).")

打球速度が100以上では, 2018年は一貫してHR率が16-17年よりも低下しているように見える. 2017年は特に確率が高かったようだ.

打球速度100 mph, 角度25°でのそれぞれの年度での推定HR確率を示す.

hr.prob%>%
  filter(launch_angle == 25, launch_speed == 100)
##   launch_angle launch_speed   pred_16   pred_17   pred_18
## 1           25          100 0.3544047 0.4064717 0.3012768

18年度は5-10%程度低く推定されている. 打球10-20に一つHRが違うというのはそれなりに大きな違いだろう.

実際のデータから24°以上26°以下, 速度98以上102以下の打球でHR確率を計算する.

typeX%>%
  filter(game_year !=  2015)%>%
  filter(launch_angle >= 24, launch_angle <= 26, launch_speed >= 98, launch_speed <= 102)%>%
  group_by(game_year)%>%
  dplyr::summarise(hr_bbe = mean(HR_FL, na.rm = TRUE))
## # A tibble: 3 x 2
##   game_year hr_bbe
##   <fct>      <dbl>
## 1 2016       0.363
## 2 2017       0.403
## 3 2018       0.307

モデルと多少の乖離があるが, 年度間での差の大きさ自体はほぼ一致している.

各年度のモデルを適用することによって推定された, 2018年の総HR数を示す.

typeX.2018 %>%
  dplyr::summarise(HR = sum(HR_FL, na.rm = TRUE),
                   Est_HR_16_model =  sum(pred_16, na.rm = TRUE),
                   Est_HR_17_model =  sum(pred_17, na.rm = TRUE),
                   Est_HR_18_model =  sum(pred_18, na.rm = TRUE))
## # A tibble: 1 x 4
##      HR Est_HR_16_model Est_HR_17_model Est_HR_18_model
##   <dbl>           <dbl>           <dbl>           <dbl>
## 1  5585           6032.           6573.           5585.

17年のモデルからは18年は6573本のHRがあったと推定されたが, 実際にはそれよりも15%以上少ない. この予想された6573本は2017年に記録された6105本よりもかなり多い (17年度のデータに17年度のモデルを適用すると6105本になるはず). また18年においては, 16年と同程度のHRが記録されたが, 16年度のモデルから予想される18年の推定HR数よりは実測値は少ない. これらの結果は, 18年度の打球は過去のデータに照らして非常に本塁打が出やすいような速度や角度を持っていたにも関わらず, その他の要因の変化によって本塁打数が抑制されたことを示唆している. 要因としてはボール自体の空力的な部分の変化, あるいは気温の変化がありうるかもしれない. 打球の方向や回転の違いが関与している可能性もあるかもしれないが, 全体では影響は小さいのではないか.

距離の推定

ついでに打球の速度と角度を使って, 各年度で打球距離を推定するモデルも作ってみる. hit_distance_scにはおかしな値があるので注意. HRの打球距離などを示す.

typeX%>%
  filter(events ==  "home_run")%>%
  arrange(hit_distance_sc)%>%
  select(game_date, events, hit_distance_sc, launch_speed, launch_angle)%>%
  head(20)
## # A tibble: 20 x 5
##    game_date  events   hit_distance_sc launch_speed launch_angle
##    <date>     <chr>              <int>        <dbl>        <dbl>
##  1 2016-08-24 home_run               0        102.          35.9
##  2 2016-09-03 home_run               0        114.          17.8
##  3 2016-09-12 home_run               0        106.          22.5
##  4 2016-09-27 home_run               0        102.          29.6
##  5 2017-05-02 home_run               0         98           39.4
##  6 2017-07-04 home_run               0        104.          27.9
##  7 2017-09-22 home_run               0        106           37.7
##  8 2015-09-02 home_run             238         74.5         30.8
##  9 2018-04-28 home_run             239         92.2         12.3
## 10 2015-07-08 home_run             255         81.4         23.2
## 11 2015-09-25 home_run             262         87.8         14.7
## 12 2018-05-29 home_run             266         89.7         16.9
## 13 2018-07-22 home_run             266         91.7         15.6
## 14 2015-07-30 home_run             286        109.          13.4
## 15 2018-03-29 home_run             293         80.5         22.9
## 16 2017-08-22 home_run             294         97.2         16.5
## 17 2017-07-29 home_run             302         90.4         38.7
## 18 2017-06-22 home_run             311         95.6         41.5
## 19 2017-04-21 home_run             312         99.8         17.7
## 20 2016-04-12 home_run             314         98.2         20.1

HRであるにも関わらず距離が0となっているデータがある. HRに限らず距離0は不自然であり, とりあえず0は除くことにする. トラッキングが完全にできていない場合, 推定値が入っている可能性もある. 推定値は異常なピークとして現れる? 全打球の距離をヒストグラムで確認する.

typeX%>%
ggplot(aes(x = hit_distance_sc)) +
  geom_histogram(binwidth = 1, 
                 alpha = 0.5)+
  facet_wrap(~game_year)

特定の距離に異常に集積しているということは無いのかもしれない. 角度は明らかに推定値を入れている異常なピークがわかりやすい (ここでは示していない) が, 主に注目したいHRになるような打球角度の範囲からはそのようなピークはかなり離れているはずなので, 除かなくてもあまり変わらないだろう. とりあえず, 残りの全打球で計算する.

各年度ごとにGAMで距離を推定させるモデルをつくる.

# 下をコメントアウトして計算
# fit16.dist <- gam(hit_distance_sc ~ s(launch_speed, launch_angle),
#                  data = typeX %>% filter(game_year ==  2016, hit_distance_sc > 0))
# fit17.dist <- gam(hit_distance_sc ~ s(launch_speed, launch_angle),
#                  data = typeX %>% filter(game_year ==  2017, hit_distance_sc > 0))
# fit18.dist <- gam(hit_distance_sc ~ s(launch_speed, launch_angle),
#                  data = typeX %>% filter(game_year ==  2018, hit_distance_sc > 0))
# 
# new.data <- expand.grid(launch_angle = seq(20, 40, 0.1),
#                         launch_speed = seq(95, 110, 5))
# dist.pred <- new.data%>%
#   mutate(pred_16 =  predict(fit16.dist,
#                              newdata = new.data),
#          pred_17 =  predict(fit17.dist,
#                              newdata = new.data),
#          pred_18 =  predict(fit18.dist,
#                              newdata = new.data))
# save(dist.pred, file = "estimated_distance_2018.Rdata")

load("estimated_distance_2018.Rdata")

HRになるような打球について, 打球速度と角度と距離の関係を示す.

dist.pred%>%
  gather(key = year, value = estimated_distance, pred_16, pred_17, pred_18)%>%
ggplot(aes(x =  launch_angle,
           y = estimated_distance,
           group = year,
           colour =  year,
           linetype = year)) +
  theme_bw(base_family = "HiraKakuPro-W3")+
  geom_line()+
  facet_wrap(~launch_speed)+
  labs(title = "打球角度+打球速度による推定飛距離の変化",
       subtitle =  "グラフの上の数字は打球速度 (mph).",
       y = "推定飛距離 (ft)")

少なくともこの速度帯では, 18年は打球速度と打球角度が同じであればボールが飛びにくかったのかもしれない. これは上でHR確率が低かったことと一致している.

打球速度100 mph, 角度25°での推定距離を示す.

dist.pred%>%
  filter(launch_angle == 25, launch_speed == 100)
##   launch_angle launch_speed  pred_16  pred_17  pred_18
## 1           25          100 382.7173 387.1339 379.4498

モデルが正しいのであればだいたい3-7 feet程度距離が違っていたということになる.

実測値と比較しておく.

typeX%>%
  filter(game_year !=  2015)%>%
  filter(launch_angle >= 24, launch_angle <= 26, launch_speed >= 98, launch_speed <= 102)%>%
  group_by(game_year)%>%
  dplyr::summarise(actual_distance = mean(hit_distance_sc, na.rm = TRUE))
## # A tibble: 3 x 2
##   game_year actual_distance
##   <fct>               <dbl>
## 1 2016                 382.
## 2 2017                 384.
## 3 2018                 378.

概ね一致している.

一応書いておくと, これは2018年のボールが飛びにくかったということには必ずしもならない (反発係数は打球速度を介して距離に影響を及ぼすため).