statcastデータにおけるlaunch_speed, launch_angleの欠測値について（後編）

↓の記事の続きです。

tsuyupon.hatenablog.com

中編では欠測値補完の手法についてまとめました。

後編では、紹介した手法のうち以下の3手法を試し、補完前後におけるlaunch_speed, launch_angleの変化を見ていきたいと思います。

ホットデック法
確率的回帰代入法
多重代入法（EMBアルゴリズム）

使用データ

前編で作成したmissingdataを使います。

※コードはこちらを参照

このlaunch_speed, launch_angleの欠測を様々な手法で補完していきます。

ちょっとだけ前処理

missingdata <- readRDS("missingdata.rds") %>% 
  mutate(hc_x = hc_x - 125.42,
         hc_y = 198.27 - hc_y, # 原点をホームベースにする
         r = sqrt(hc_x^2 + hc_y^2),
         theta = atan(hc_y/hc_x), # 極座標変換
         flg = ifelse(is.na(flg), 0, flg),　# 欠測フラグ
         flg2 = ifelse(is.na(flg2), 0, flg2)) %>%　# 欠測種類フラグ 
  filter(!is.na(bb_type)) # 打球結果が入っていない打球が1球あったので除外…

# データ列を絞っておく
dat <- missingdata %>% 
  select(park, bb_type, events, r, theta, zone, launch_speed, launch_angle)

結果の解釈がしやすいよう打球位置は極座標変換後の値を使いました。

ホットデック法

※手法詳細はこちらを参照ください

library(simputation)
m1 <- impute_shd(dat, launch_speed+launch_angle ~ park+bb_type+events+r+theta+zone, pool = "complete", backend = "VIM")

これだけ。

simputationパッケージにはimpute_shd（同じ区分内における直前の観測値で補完する方法）とimpute_rhd（同じ区分内の観測値からランダムなレコードを選んで補完する方法）がありますが、後者は時間がかかってしまうため、前者を採用しました。

補完した結果として得られた分布がこちらです。

※flg2=0が欠測していなかったレコード

f:id:tsuyu_pon:20200420154935p:plain

欠測の種類別に補完されていて、滑らかにうまく補完できていそうです。

確率的回帰代入法

つづいて、確率的回帰代入法。

※手法詳細はこちらを参照ください

m2 <- impute_lm(dat, launch_speed+launch_angle ~ park+bb_type+events+r+theta+zone, add_residual = "normal")

これだけ。

formulaでは目的変数を2つにしていますがそれぞれに対し、代入モデルを構築し、補完する値を計算しています。

線形回帰で計算しているため、ホットデックの例よりもlaunch_angleの70付近の値が多く補完されているようです。

補完した結果として得られた分布がこちらです。

※flg2=0が欠測していなかったレコード

f:id:tsuyu_pon:20200420162809p:plain

ホットデックではレコード単位で見て、補完値まで決めましたが、回帰代入の方法ではlaunch_speed, launch_angle別々で考えているので、分布ではきれいに表現できてもlaunch_speed, launch_angleの組み合わせは上手くいっていない可能性もあります。

多重代入法（EMBアルゴリズム）

※手法詳細はこちらを参照ください

library(amelia)
set.seed(42)
m3 <- amelia(dat, m=100, noms = c("park", "bb_type", "events", "zone"))
# overimpute(m3, var = 7) # データが重すぎるので断念
# overimpute(m3, var = 8) # データが重すぎるので断念
# disperse(m3, dims = 1, m = 100) # データが重すぎるので断念
# m3_2 <- datlist2mids(m3$imputations)
# densityplot(m3_2) # データが重すぎるので一部のみで実行

# 100通りの代入値から得た平均と標準偏差から正規乱数を発生させる
# 本来はここで100通りの回帰モデルを構築して回帰係数の評価などを行う
ls_m <- rowMeans(bind_cols(lapply(m3$imputations, function(x) x$launch_speed)))
ls_sd <- apply(bind_cols(lapply(m3$imputations, function(x) x$launch_speed)), 1, sd)
set.seed(42)
ls <- mapply(function(mean, sd) rnorm(1, mean, sd), ls_m, ls_sd)

la_m <- rowMeans(bind_cols(lapply(m3$imputations, function(x) x$launch_angle)))
la_sd <- apply(bind_cols(lapply(m3$imputations, function(x) x$launch_angle)), 1, sd)
set.seed(42)
la <- mapply(function(mean, sd) rnorm(1, mean, sd), la_m, la_sd)

dat3 <- missingdata %>% 
  mutate(launch_speed = ls, launch_angle = la)

f:id:tsuyu_pon:20200420172129p:plain

シミュレーションごとのdensity plotを見ると欠測が多かったところに山が出来ているのがわかるので、共変量の選定は筋が悪くなかったようです。

f:id:tsuyu_pon:20200420172035p:plain

分布は、確率的回帰代入のような代入値を作り方をしたので先ほどの結果と似たものになりました。

ちなみに100通りの代入値の平均をとった場合は、次のような分布になります。

※flg2=0が欠測していなかったレコード

f:id:tsuyu_pon:20200420180939p:plain

↑もともとのデータセットのように大きな山がいくつか出てきてしまっています。

f:id:tsuyu_pon:20200402133601j:plain

↑これは前編で載せた図で元のデータセットから出した分布ですが、この図で異常に頻度が高かったあたりに補完後の山が出来ているので、元のデータセットをつくるために使った共変量と今回使った共変量が似ているのかもしれません。

結果の比較

ここまで3手法を試してきましたが、それぞれ違った特徴がありました。

この結果を使って、実際に欠測ありのデータにおける平均値からどれだけ動いたかを様々な粒度で確かめてみました。

Yankees

f:id:tsuyu_pon:20200420221421p:plain

Astros

f:id:tsuyu_pon:20200420221503p:plain

Christian Yelich

f:id:tsuyu_pon:20200420221601p:plain

José Altuve

f:id:tsuyu_pon:20200420221622p:plain

チームごとに集計しても選手ごとに集計しても元のデータセットから算出した平均値との差異は大きくありませんでした。

ただ、全体的に言えることですが、以下のような傾向があるようです。

回数は少ないが、欠測のあった犠牲バントやフィルダースチョイスの動きが大きい
確率的回帰代入と多重代入の結果は比較的似ているが、ホットデック法の結果は極端に動くことが多い
欠測が多かったfield_out（多くはポップフライ）は、動きが少ないが補完結果の動く方向はどの補完方法でも同じで、打球速度が下がり打球角度が上がる

つまり、ポップフライが多かった選手は補完によって、打球速度が下がり打球角度が上がることになります。

どれが正しいかは欠測なのでわかりませんが、手法のメリット・デメリットについて以下のようにまとめることが出来ます。

手法	計算速度	単純集計	分析モデルへの応用
ホットデック法	〇	〇	△
確率的回帰代入法	〇	〇	△
多重代入法	×	〇	〇

ヒットかどうかのモデリング

せっかく多重代入したものが出来たので、打球がヒットになったかどうかを、launch_speed, launch_angleを説明変数（の一部）としてモデリングしてみました。

> # 多重代入法での回帰モデル
> m3_2 <- datlist2mids(m3$imputations)
> modelA <- glm.mids(events %in% c("home_run", "single", "double", "triple") ~ 
+                    launch_speed + launch_angle + bb_type, data = m3_2, family = binomial)
> summary(pool(modelA))
                term     estimate    std.error   statistic        df      p.value
1        (Intercept) -4.458569506 0.0321505225 -138.677980  69446.39 0.000000e+00
2       launch_speed  0.039167041 0.0003261071  120.104846  55910.64 0.000000e+00
3       launch_angle -0.003422507 0.0003618348   -9.458756  27450.93 0.000000e+00
4 bb_typeground_ball -0.105359051 0.0190754693   -5.523274  49009.96 3.344248e-08
5  bb_typeline_drive  1.396924883 0.0124204220  112.470001 172571.68 0.000000e+00
6       bb_typepopup -2.114878632 0.0436746879  -48.423440 329411.25 0.000000e+00


> # 欠測ありデータでの回帰モデル
> modelB <- glm(events %in% c("home_run", "single", "double", "triple") ~ 
+                 launch_speed + launch_angle + bb_type, data = dat, family = binomial)
> summary(modelB)

Call:
glm(formula = events %in% c("home_run", "single", "double", "triple") ~ 
    launch_speed + launch_angle + bb_type, family = binomial, 
    data = dat)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9134  -0.8420  -0.6082   0.9502   3.1933  

Coefficients:
                     Estimate Std. Error  z value Pr(>|z|)    
(Intercept)        -4.4654444  0.0330124 -135.265  < 2e-16 ***
launch_speed        0.0397393  0.0003335  119.151  < 2e-16 ***
launch_angle       -0.0041502  0.0003615  -11.480  < 2e-16 ***
bb_typeground_ball -0.1382680  0.0191057   -7.237 4.59e-13 ***
bb_typeline_drive   1.3637115  0.0124485  109.549  < 2e-16 ***
bb_typepopup       -1.5938762  0.0444717  -35.840  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 444369  on 342670  degrees of freedom
Residual deviance: 381670  on 342665  degrees of freedom
  (34156 observations deleted due to missingness)
AIC: 381682

Number of Fisher Scoring iterations: 5