wOBAの誤差範囲: スキルの推定 (95% CI, 回帰)
Author: @sleep_in_nmbrs
2019-09-22
1 背景と目的
以前のポストで多項分布の真のパラメータが定まった状態で, ランダムサンプリングによって得られたwOBAがどの程度変動するかを調べた. このような変動の存在は, 特にサンプルサイズが限定された状況において, 実際に測定されたデータから真の確率を推定する場合に不確実性が伴い, 単純な平均値を考えるだけでは不十分であることを示唆する.
1.1 データから確率を推定する
上で触れたポストでは各イベントの確率が定まった状態を使ったが, このような状況を考えることは奇異に映るかもしれない. 現実の状況では真の確率分布を知るすべがあることは非常に珍しく, むしろ得られたデータを利用して真の確率やその分布を推定したいことが多い. サンプルサイズが小さい状態で不確実性が問題になるのはデータから推定する場合でも変わらない. このような推定を可視化が容易な二項分布における尤度を使った議論でも示す.
ある特定の打者が出塁率0.4を記録した場合を考える. この打者の能力に関する事前の情報は手元に無いものとする (これも現実にはずいぶん馬鹿げた仮定だがお付き合い願いたい. Appendixでは, 事前の情報をベイズ統計で考慮した場合の例を示した). ここで, n回試行でk回 (= n * 0.4) 成功が得られたという条件における, 真の確率を表すパラメータ\(\theta\)の確からしさを示す尤度関数\(f(\theta\,|\,n,k)\)は以下のようになる.
\[ f(\theta\,|\,n,k) = \binom{n}{k}\theta^k(1-\theta)^{n-k}\]
ｎ=100の時の, \(\theta\)に対する尤度をプロットする.
require(tidyverse)
require(gt)
L <- function(p, n = 100, k = 40) {
choose(n, k) * p^k * (1 - p) ^ (n-k)
}
par(family="HiraKakuProN-W3")
probs <- seq(0.3, 0.5, by = 0.001)
plot(probs, L(probs, n = 100, k = 40),
xlab = expression(theta),
ylab = "尤度",
type = "l",
main = "N = 100")
これはパラメータ\(\theta\)の値への妥当な確信の程度を示す. 結果の平均値である0.4が最尤推定値であるが, 0.35や0.45などが真の確率\(\theta\)である可能性もそれなりあることがわかる. 複雑な事を言っていると思うかもしれないが, 実のところこれは, 単にそれぞれの\(\theta\)の値の時に二項分布によって期待される確率をプロットして相対的に比較しているだけに過ぎない. \(\binom{n}{k}\theta^k(1-\theta)^{n-k}\)を, 特定の定まった結果を持った条件で\(\theta\)を推定する式と見れば尤度関数であり, 特定の定まった\(\theta\)のもとでn回試行してデータが得られる確率と見れば確率密度関数になる.
n=500について同様にプロットする.
par(family="HiraKakuProN-W3")
plot(probs, L(probs, n = 500, k = 200),
xlab = expression(theta),
ylab = "尤度",
type = "l",
main = "N = 500")
n=500であれば, 尤度関数の曲線はn=100に比べて凸になり, \(\theta\)が0.35や0.45である可能性は, 0.40である可能性に比べてかなり小さく, \(\theta\)は0.40に比較的近いと予測される. 真の確率に関してどの程度の確信を持つことができるかはサンプルサイズに依存することがわかる. これは, 真の確率が既知であった場合において, サンプルサイズが大きくなればその誤差が小さくなることに対応している.
1.2 95%信頼区間と母平均の推定
不確実性が推定に影響を与えるのであれば, 持っている不確実性を表現する方法が必要となる. 多くの場合において, 我々が興味を持っているのは真の平均値, つまり母平均値だと思われる. この母平均値の不確実性を提示する区間推定的な手法の中で, 最も簡単に計算可能なものは95%信頼区間だろう.
下の式は, 平均を中心として, よく使われる (1標本の) 母平均の95%信頼区間 (CI) の範囲を示している.
\[ Mean\,± SE\,*\,t_{critical} \]
ここでSEは標準誤差, \(t_{critical}\)はt値の臨界値でサンプルサイズによって変化する. \(t_{critical}\)については, 野球のデータではよほど細かくデータを分割しない限りサンプルは多いため, 両側5%水準であれば概ね1.96程度であり, 正確に用いることが困難な場合1.96あるいは2を使えば十分に良い近似となる.
このようにして求められた信頼区間は, 一定の条件を満たしていれば, 得られた信頼区間の95%は母平均を含んでいるという性質を持つ. つまり100個の信頼区間を作れば, 期待としてはそれらのうち95個は母平均を含んでいる. 「母平均が95%の確率で含まれる」という説明をたまに見るが, このような性質を持つを区間を得たい場合は, ベイズの信用区間が必要となる (母平均の値は未知であっても定数である前提なので, 信頼区間に入っているか入っていないかのどちらかでしかない).
1.3 平均への回帰とタレントの推定
この母平均の推定は, よりSabermetrics的に表現するなら, そのデータの元となった選手や状況における, スキルやタレント, つまり真の能力的なものを推定していると言えるだろう. Sabermetricsにおいて, タレントの推定に使われる方法としては, 平均への回帰の影響を取り込む方法がある (The BOOKのappendixに詳しい説明がある; 日本語ではbbconcreteさんによるタンゴの回帰式). “平均への回帰”という言葉自体は様々な文脈で広く使われているが, このポストでは以下で述べる式を利用して回帰量や能力を推定する方法を指すことにする.
回帰の大きさが指標によって異なることはよく知られている. それぞれの指標の回帰の大きさを表現するためによく利用される数値として, 平均的な成績からの乖離の半分を能力とみなすことができる打席数がある. この打席数が小さいほど, 少ない打席数で本人の能力が推定できるということになる. この数値は各指標によって大きさが異なり, 眺めることで指標の回帰の大きさの違いをなんとなく見て取ることができる. しかし, この数値の本当の利用価値は, 適切な回帰の量を考慮した母平均値の推定値を簡単な計算で求めることができることにある.
実際に観察されたある選手の数値を\(Observed\ value\), 打席数を\(PA_{batter}\), 回帰すべき平均を\(Average\)とすると, スキルの点推定値\(Estimated\, skill\)は以下の式で求めることができる.
\[Estimated\ skill = Average + (Observed\ value - Average) * \frac{PA_{batter}}{PA_{batter} + PA_{even}}\]
ここで\(PA_{even}\)は前述した平均値からの乖離の半分が打者の能力とみなす事が可能な打席数で, 例えばMLBの数年分のデータなどから計算することができる (必ずしもPAを使うわけではない. 打球数を分母とする指標であれば打球数を単位として使ったほうが自然). 右辺は:
平均値
平均からの乖離の内どれだけが打者の能力とみなすことができるか
の2つの要素で構成されていると解釈できる. \(PA_{batter}\)が\(PA_{even}\)に比べて十分に大きくなれば\(\frac{PA_{batter}}{PA_{batter} + PA_{even}}\)が1に近づくため, \(Estimated\ skill\)はその時点で得られているデータ\(Observed\ value\)とほぼ一致する. PAが0であれば, 推定値は\(Average\)と一致する. 利用可能な情報がない状態では, 平均を予測値とするのは理にかなっているだろう.
計算例を示す. ある指標の平均値が0.33, \(PA_{evet}\)が350とする. ある選手が500 PAで0.4を記録した. この時回帰の大きさを考慮した推定値は以下のように求めることができる.
PA.batter <- 500
PA.even <- 350
0.33 + (0.4 - 0.33) * (PA.batter / (PA.batter + PA.even))## [1] 0.3711765この式は:
事前に持っている情報, つまり回帰する先となる値\(Average\)と, 回帰の大きさを示す\(PA_{even}\), そして
その選手のデータ (記録された値\(Observed\ value\)とサンプルサイズ\(PA_{batter}\))
という要素から構成されており, 事前の情報と, 特定の選手の測定を統合するイメージである. (実際, この式は、メタ分析などで用いられる複数の測定結果を分散の逆数で重み付けする方法の式から導出できる (The BOOK, pp. 380で記述されているがAppendix 2にも書いた; この説明がベイズっぽいと思った方はAppendix 1もどうぞ).).
一つの問題は\(PA_{even}\)をどうやって推定するかだが, 具体的な方法としては, 計測の信頼性 (テスト理論的な意味で) を利用して計算する場合と, 成績表などからランダムサンプリングを利用する場合などがある. 前者としてはThe BOOKのAppendixやSean Dolinarのこの記事が, 後者としてはTangoやpizzacutterこと Russel Carletonの説明が参考になるだろう.
1.4 本稿の目的
95% CI, あるいは平均への回帰を, ランダムサンプリングされたデータや実際のデータに適用してみる. ここでは, 前回のポスト同様wOBAに注目する. これは言うまでもなく, 現状 (おそらくこれからもずっと) wOBAが最も重要な打撃指標であるからである. wOBAは非常に汚い標本分布となるが, この指標で95% CIが使えるのかどうか検討する. まず正規分布から得られたwOBAっぽい数値を使って95% CIの性質を確認し (section 2.1), その後 多項分布からのサンプリング結果 (1打席ごとの結果のシミュレーション) を使って, 母平均値が実際に区間に含まれているかを調べる (section 2.2). その後 (section 2.3), 恣意的に選んだ例に関して実データにも適用してみて, 回帰の大きさから推測する方法と比較して示す. 結論から言えば, 可能なら回帰を利用したほうが色々と利点が多い.
2 結果と考察
2.1 95% CI: 正規分布からのサンプリング
wOBAの話に入る前に, まず95%信頼区間の性質を確認する. 正規分布から取り出した例で, 推定された信頼区間の95%が平均を含むことを示す. せっかくなので野球のデータっぽくする. 以下のコードでは規定到達打者っぽい, 平均0.33, sd = 0.03の正規分布から300サンプリングし, それを100,000回繰り返している. 最初と最後の3行を示す.
seasons = 1:100000
random.data <- vector("list", length(seasons))
for(i in seq_along(seasons)){
temp <- data.frame(wOBA =rnorm(300,
mean = 0.33, # 観測された平均値のイメージ
sd = 0.03)) # 観測されたSDのイメージ
temp$year <- i
random.data[[i]] <- temp
}
random.data <- bind_rows(random.data)
random.data%>%
head(3)%>%
mutate_if(is.numeric, funs(round), 3)%>%
gt()| wOBA | year |
|---|---|
| 0.359 | 1 |
| 0.302 | 1 |
| 0.387 | 1 |
random.data%>%
tail(3)%>%
mutate_if(is.numeric, funs(round), 3)%>%
gt()
| wOBA | year |
|---|---|
| 0.281 | 1e+05 |
| 0.355 | 1e+05 |
| 0.307 | 1e+05 |
このデータの各年度について95%信頼区間を計算する. この場合, 集団の母平均 (あるいは平均的なスキル) がどの範囲にあるかを考える. ランダムに選んだ10年分の計算結果を示す.
# rではt値はqt()で取り出すことができる
# 片側0.025, 自由度はサンプルサイズ N - 1
random.data %>%
group_by(year)%>%
dplyr::summarise(N = n(),
Avg = mean(wOBA),
SD = sd(wOBA),
SE = SD /sqrt(N),
t_crtcl = qt(0.975, df = N-1),
CI_95 = SE * t_crtcl,
CI_upper = Avg + CI_95,
CI_lower = Avg - CI_95)%>%
sample_n(10)%>%
mutate_if(is.numeric, funs(round), 4)%>%
gt()| year | N | Avg | SD | SE | t_crtcl | CI_95 | CI_upper | CI_lower |
|---|---|---|---|---|---|---|---|---|
| 39338 | 300 | 0.3288 | 0.0289 | 0.0017 | 1.9679 | 0.0033 | 0.3321 | 0.3255 |
| 57294 | 300 | 0.3319 | 0.0299 | 0.0017 | 1.9679 | 0.0034 | 0.3353 | 0.3285 |
| 26201 | 300 | 0.3332 | 0.0310 | 0.0018 | 1.9679 | 0.0035 | 0.3367 | 0.3297 |
| 23463 | 300 | 0.3271 | 0.0304 | 0.0018 | 1.9679 | 0.0034 | 0.3305 | 0.3236 |
| 82919 | 300 | 0.3292 | 0.0300 | 0.0017 | 1.9679 | 0.0034 | 0.3326 | 0.3258 |
| 53117 | 300 | 0.3289 | 0.0309 | 0.0018 | 1.9679 | 0.0035 | 0.3324 | 0.3254 |
| 58561 | 300 | 0.3306 | 0.0324 | 0.0019 | 1.9679 | 0.0037 | 0.3343 | 0.3269 |
| 31022 | 300 | 0.3299 | 0.0303 | 0.0017 | 1.9679 | 0.0034 | 0.3333 | 0.3264 |
| 64833 | 300 | 0.3282 | 0.0304 | 0.0018 | 1.9679 | 0.0035 | 0.3316 | 0.3247 |
| 30662 | 300 | 0.3292 | 0.0314 | 0.0018 | 1.9679 | 0.0036 | 0.3328 | 0.3257 |
N=300あるが, それでもAvgが0.33から少し離れているケースもあるかもしれない.
ランダムに選んだ各打者 (100,000x) の平均WOBAっぽい数値の分布をヒストグラムで示す.
random.data%>%
sample_n(100000)%>%
ggplot(aes(x = wOBA))+
geom_histogram(alpha = 0.7,
binwidth = 0.01,
colour="black", fill="gray")+
theme_bw(base_family = "HiraKakuPro-W3")+
labs(y = "生起回数")
これらのデータから計算された, 各年度での信頼区間のうち, 母平均0.33が含まれている割合を計算する.
random.data %>%
group_by(year)%>%
dplyr::summarise(N = n(),
Avg = mean(wOBA),
SD = sd(wOBA),
SE = SD /sqrt(N),
t_crtcl = qt(0.975, df = N-1),
CI_95 = SE * t_crtcl,
CI_upper = Avg + CI_95,
CI_lower = Avg - CI_95)%>%
mutate(within_CI = ifelse(0.33 <= CI_upper & 0.33 >= CI_lower, 1, 0))%>%
dplyr::summarise(within_CI_pct = mean(within_CI) * 100)%>%
gt()| within_CI_pct |
|---|
| 95.077 |
確かに概ね95%になる.
2.2 95% CI: 多項分布からのサンプリング
1打席ごとの結果を多項分布からサンプリングしたwOBAに関しても成立するかどうかを確認する.
関数の設定.
sim.wOBA.PA <- function(df = average.hitter,
PA = 100,
season = 50000){
seasons = 1:season
seasons <- vector("list", length(seasons)) # 空のリスト. これに入れていって最後にbind_rows
for(i in seq_along(seasons)){
the.season <- data.frame(t(rmultinom(PA, 1, df$N))) # 今回は1打席ずつ記録する
the.season$year = i
seasons[[i]] <- the.season
}
seasons <- bind_rows(seasons)
names(seasons) <- c(as.character(df$events), "year")
season.stats <- seasons %>%
group_by(year)%>%
mutate(woba_value = single * 0.9 + double * 1.2 + triple * 1.6 + home_run * 2.0 +
field_error * 0.9 + walk * 0.7 + hit_by_pitch * 0.7 + woba_0 * 0)%>%
dplyr::summarise(N = n(),
Avg = mean(woba_value),
SD = sd(woba_value),
SE = SD /sqrt(N),
t_crtcl = qt(0.975, df = N-1),
CI_95 = SE * t_crtcl,
CI_upper = Avg + CI_95,
CI_lower = Avg - CI_95)
}中身の一部に関して例を示す. 下の例では平均的な打者の1000打席をシミュレーションして, そのうちの10打席だけを示している.
test.season <- data.frame(t(rmultinom(1000, 1, average.hitter$N))) # 設定したPA分を計算
names(test.season) <- as.character(average.hitter$events)
test.season%>%
mutate(woba_value = single * 0.9 + double * 1.2 + triple * 1.6 + home_run * 2.0 +
field_error * 0.9 + walk * 0.7 + hit_by_pitch * 0.7 + woba_0 * 0)%>%
sample_n(10)%>%
gt()| single | double | triple | home_run | field_error | walk | hit_by_pitch | woba_0 | woba_value |
|---|---|---|---|---|---|---|---|---|
| 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0.9 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0.0 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0.0 |
| 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 2.0 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0.0 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0.0 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0.0 |
| 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 2.0 |
| 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1.2 |
| 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0.7 |
wOBA valueを横軸にとり, その分布をヒストグラムで示す.
test.season%>%
mutate(woba_value = single * 0.9 + double * 1.2 + triple * 1.6 + home_run * 2.0 +
field_error * 0.9 + walk * 0.7 + hit_by_pitch * 0.7 + woba_0 * 0)%>%
ggplot(aes(x = woba_value))+
geom_histogram(alpha = 0.7,
binwidth = 0.1,
colour="black", fill="gray")+
theme_bw(base_family = "HiraKakuPro-W3")+
labs(y = "発生回数")
それぞれの打席結果の平均的な価値を与えているため非常に離散的である. また, アウト (value = 0) が多く, ひどく不自然な分布である. この不自然なものの平均値が, この平均的な打者の1000打席のwOBAとなる.
このデータから信頼区間を求める.
test.season%>%
mutate(woba_value = single * 0.9 + double * 1.2 + triple * 1.6 + home_run * 2.0 +
field_error * 0.9 + walk * 0.7 + hit_by_pitch * 0.7 + woba_0 * 0)%>%
dplyr::summarise(N = n(),
Avg = mean(woba_value),
SD = sd(woba_value),
SE = SD /sqrt(N),
t_crtcl = qt(0.975, df = N-1),
CI_95 = SE * t_crtcl,
CI_upper = Avg + CI_95,
CI_lower = Avg - CI_95)%>%
mutate_if(is_numeric, funs(round), 3)%>%
gt()| N | Avg | SD | SE | t_crtcl | CI_95 | CI_upper | CI_lower |
|---|---|---|---|---|---|---|---|
| 1000 | 0.318 | 0.518 | 0.016 | 1.962 | 0.032 | 0.35 | 0.285 |
このような不自然な分布のデータでも95% CIは使い物になるのだろうか? 結論から言うと, ある程度のサンプルサイズがあれば95% CIとしての性質を十分に示す. 少なくともこのシミュレーションでは明らかに中心極限定理 (central limiting theorem) の前提を満たしているため, サンプル数が大きければ, その“平均値の分布”, つまりある選手のwOBAとして測定される値の分布は正規分布に近づくことが容易に予想できる. 正規分布に近似できるのであれば, 基本的には上で示した正規分布からwOBAっぽい値を取り出した例とあまり変わらないはずである.
打席数30, 150, 1000で, 50000人分の信頼区間を計算する. 母平均値はすべての打者のwOBAの平均の平均をとる. 厳密には母平均では無いが, サンプルが大きいので, 近似としては十分だろう.
# avg.hitter.result <- c(30, 150, 1000) %>%
# map(~sim.wOBA.PA(df = average.hitter, PA = ., season = 50000))%>%
# bind_rows()
# save(avg.hitter.result, file="avg_hitter_95CI.Rdata")
load("avg_hitter_95CI.Rdata")
avg.hitter.result <- avg.hitter.result%>%
mutate(Diff = abs(mean(avg.hitter.result$Avg) - Avg),
within_SE = ifelse(Diff <= SE, 1, 0),
within_CI = ifelse(Diff <= SE * t_crtcl, 1, 0))
avg.hitter.result %>%
group_by(N)%>%
dplyr::summarise(Season = n(),
within_SE_pct = mean(within_SE) * 100,
within_CI_pct = mean(within_CI) * 100)%>%
rename(PA = N)%>%
mutate_if(is_numeric, funs(round), 3)%>%
gt()| PA | Season | within_SE_pct | within_CI_pct |
|---|---|---|---|
| 30 | 50000 | 67.072 | 93.512 |
| 150 | 50000 | 67.902 | 94.734 |
| 1000 | 50000 | 68.464 | 95.034 |
95%信頼区間的な性質はそれなりに持っていると言えそうであり, 目安としては悪くないようだ. 打席数が少ない状態では信頼区間に含まれる割合が少し低い.
ここで求めた信頼区間は, サンプルの平均値, つまりwOBA (ここではAvg列) が正規分布していることを仮定している. 中心極限定理は非常に強力で, 様々な標本分布で機能するが, サンプルサイズがある程度大きいことが必要である. 打席数が少ない状態では, ここで求めたwOBAの分布は正規分布から少し離れたものであったかもしれない.
一応density plotを示しておくが, 正規分布からの乖離はこの程度の方法では明らかではない.
avg.hitter.result%>%
mutate(N = factor(N))%>%
ggplot(aes(x = Avg, color = N, fill = N))+
geom_density(stat = "density",
position = "identity",
alpha = 0.5,
adjust = 1)+
theme_bw(base_family = "HiraKakuPro-W3")+
labs(x = "wOBA", y = "確率密度",
title = "wOBAの分布.")
平均wOBAとその信頼区間の幅をプロットする.
avg.hitter.result%>%
sample_n(5000)%>%
ggplot(aes(x = Avg,
y = CI_95))+
geom_point(size = 0.1,
alpha = 0.6)+
facet_wrap(~N)+
theme_bw(base_family = "HiraKakuPro-W3")+
labs(x = "wOBA", y = "信頼区間の幅 (片側)")
30 PAでは極端に高い or 低い結果, 例えばwOBA0.5以上, も偶然によって生じるが, その場合は信頼区間の幅も広いため, 母平均値を含んでいたと考えられる.
ランダムに選んだ信頼区間300とデータ全体の平均の関係.
avg.hitter.result%>%
sample_n(300)%>%
ggplot(aes(x = year,
y = Avg))+
geom_point(size = 0.2)+
geom_errorbar(aes(ymin = Avg - CI_95,
ymax = Avg + CI_95),
alpha=0.6,
size=0.5)+
geom_hline(yintercept = mean(avg.hitter.result$Avg), linetype = 2, colour ="red")+
facet_wrap(~N)+
theme_bw(base_family = "HiraKakuPro-W3")+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
labs(x = "", y = "wOBA",
subtitle ="平均wOBA ± 95%信頼区間. 赤の点線はデータ全体の平均.")
ほとんどの信頼区間はデータ全体の平均を含んでいる. サンプルサイズを増やすほど (SEが小さくなるため) 信頼区間の幅は狭くなり, 不確実性が低下していることが確認できる.
平均を含まない信頼区間もある. 含まないものを300選んで同様に示す.
avg.hitter.result%>%
filter(within_CI == 0)%>%
sample_n(300)%>%
ggplot(aes(x = year,
y = Avg))+
geom_point(size = 0.2)+
geom_errorbar(aes(ymin = Avg - CI_95,
ymax = Avg + CI_95),
alpha=0.5,
size=0.5)+
geom_hline(yintercept = 0.33, linetype = 2, colour ="red")+
facet_wrap(~N)+
theme_bw(base_family = "HiraKakuPro-W3")+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
labs(x = "", y = "wOBA",
subtitle ="平均wOBA ± 95%信頼区間. 赤の点線はデータ全体の平均.
信頼区間にデータ全体の平均が含まれなかった場合の区間を示す.")
全体的に平均の下側にぶれてしまうことが多い. これも分布の歪みによるものだろう. bootstrap CIの方が歪みには強いと思われるので, そちらを試してみても面白いかもしれない (面倒).
2.3 実データ
ここからはStatcastデータを利用して, 実際のデータに95% CIと回帰を適用してみる. 推定における不確実性の問題が特に大きいのは, 打席数が少ない状態である (もう少し厳密に言えば, サンプルサイズの大きさとデータが発生する際に従う分布に依存するデータのばらつきの大きさ, それに対する選手間の能力差, これらの間のバランスで決まる. 能力差が大きいならサンプルサイズが小さくても問題は起きにくい). 典型的には個人の状況別成績において問題が生じる. 以下では, よく知られた打者の状況別成績として:
- 左右別対戦成績
- 得点圏打撃成績
を見ていく. 使う指標は当然wOBAである. ここでは, 大した意味もなく選んだ特定の選手を取り上げて比較する. これは特定の選手の状況別成績を取り上げて, その能力について何らかの結論を下すタイプの記事を真似ようと努力しただけであり, そもそもこのような方法では状況別成績の価値についての評価は困難である. これらの状況別成績がどの程度予測力があるか, あるいは状況別能力における個人差がどの程度と推定されるか, という点についてはThe BOOKが詳しく検討しているのでそちらを参照していただきたい.
2.3.1 左右の相性
MLB 15-18における打席の左右 (stand) と, 投手が投球に使った腕の左右 (p_throws; 利き腕とは限らないし, 両投げもありうるのでこんな怪しげな日本語になっている) ごとに, 打席でのwOBAとその信頼区間を示す.以下の計算ではIBBは含まれていない. またbuntに関連した打席の大部分は除かれている.
load("sc_data_for_CI_calc.Rdata")
data %>%
filter(woba_denom == 1)%>%
filter(! grepl("bunt", des))%>%
group_by(stand, p_throws)%>%
dplyr::summarise(N = n(),
Avg = mean(woba_value),
SD = sd(woba_value),
SE = SD /sqrt(N),
t_crtcl = qt(0.975, df = N-1),
CI_95 = SE * t_crtcl,
CI_upper = Avg + CI_95,
CI_lower = Avg - CI_95)%>%
rename(wOBA = Avg, PA = N)%>%
mutate_if(is.numeric, funs(round), 3)%>%
select(-c(SD, SE, t_crtcl))%>%
gt()| p_throws | PA | wOBA | CI_95 | CI_upper | CI_lower |
|---|---|---|---|---|---|
| L | |||||
| L | 54952 | 0.305 | 0.004 | 0.310 | 0.301 |
| R | 246052 | 0.336 | 0.002 | 0.338 | 0.334 |
| R | |||||
| L | 139276 | 0.336 | 0.003 | 0.339 | 0.333 |
| R | 285581 | 0.322 | 0.002 | 0.324 | 0.320 |
# woba_denom == 1に制限しているためIBBは除かれている
# buntもだいたいは除かれているこれは左右に関してのMLBでの平均的な効果を示す. 左上に示された投手の使用した腕 (p_throws) に関して, テーブルが上下に分割されている. 右投手を基準とすると, 左打者は左投手に対して-0.031 (0.305 - 0.336) のdisadvantageを, 右打者は左投手に対して0.014のadvantageを持つ. 95% CIを見ると, MLB全体4年分であるにも関わらず, 左投手vs左打者は他に比べてやや大きな誤差が残っていることがわかる. これは左投手, 左打者ともに右に比べて比率が少ないことと, 左打者が選択的に右投手に対して起用されていることによって、 サンプルサイズがやや小さくなっているため.
なんとなく目についた個人で適用してみる. 2017年のLorenzo Cain (RHH) の対左右別成績を示す. なぜ彼を選んだかといえば, 単に投手左右別対戦成績で近い成績の選手を適当に選んだだけである.
data %>%
filter(woba_denom == 1)%>%
filter(! grepl("bunt", des))%>%
filter(game_year == "2017")%>%
filter(player_name == "Lorenzo Cain")%>%
group_by(player_name, p_throws)%>%
dplyr::summarise(N = n(),
Avg = mean(woba_value),
SD = sd(woba_value),
SE = SD /sqrt(N),
t_crtcl = qt(0.975, df = N-1),
CI_95 = SE * t_crtcl,
CI_upper = Avg + CI_95,
CI_lower = Avg - CI_95)%>%
rename(wOBA = Avg, PA = N)%>%
mutate_if(is.numeric, funs(round), 3)%>%
select(-c(SD, SE, t_crtcl))%>%
gt()| p_throws | PA | wOBA | CI_95 | CI_upper | CI_lower |
|---|---|---|---|---|---|
| Lorenzo Cain | |||||
| L | 143 | 0.362 | 0.090 | 0.452 | 0.272 |
| R | 501 | 0.371 | 0.044 | 0.415 | 0.327 |
平均wOBAは同程度だが, 左投手のほうが打席数がかなり少ないため信頼区間の上限と下限が広い. 母平均の推定という意味では, 仮に同じwOBAであってもサンプルサイズが違えばその数値の意味は同じではない.
次に, 2017年のMcCutchen (RHH)の結果を示す. なぜ彼を選んだかといえば, 単に成績表を眺めて差が大きい選手を探してきただけである.
| 2017 | |||||
|---|---|---|---|---|---|
| p_throws | PA | wOBA | CI_95 | CI_upper | CI_lower |
| Andrew McCutchen | |||||
| L | 145 | 0.470 | 0.103 | 0.572 | 0.367 |
| R | 500 | 0.342 | 0.046 | 0.387 | 0.296 |
wOBAとして非常に大きな差 (0.128) があるが, これでも信頼区間は重なっている. これは一つにはそもそも左投手が少ないため, 対左の打席が少なく, どうしても推定が甘くなるためである. このように区間があると, 2つの状況について, それぞれの信頼区間の重なりの有無で何らかの意味を見出したくなるのがヒトの性というものかもしれない. が, この信頼区間の重なりがあるかどうかという二分法の意味はよく誤解される. 95%信頼区間の重なりを元にした有意性の判断は頻度論的な帰無仮説検定に比べてかなり保守的である (つまり有意差が出にくい) とされており, 重なりがあるから差がないとはしないほうがいいかもしれない (とはいえ最近は頻度論的な帰無仮説検定のp<.05は緩すぎるので, 例えばp<.005にしろとかいう意見が色々な方面から出ていることを考えると, より保守的で好ましいとも言えるかもしれない; そもそも二分法は可能なら避けたほうがいいというのが正しいツッコミだろう).
ついでなので無駄にさらに球種別に分けてみる.
| wOBA of McCutchen 2017 by p_throws and pitch types | |||||
|---|---|---|---|---|---|
| p_throws | pitch_type | PA | wOBA | CI_upper | CI_lower |
| L | CH | 21 | 0.321 | 0.519 | 0.123 |
| L | CU | 10 | 0.430 | 0.836 | 0.024 |
| L | FC | 6 | 0.483 | 1.350 | -0.383 |
| L | FF | 64 | 0.526 | 0.685 | 0.367 |
| L | FS | 3 | 0.667 | 3.535 | -2.202 |
| L | FT | 22 | 0.441 | 0.722 | 0.160 |
| L | KC | 3 | 0.000 | 0.000 | 0.000 |
| L | SI | 4 | 0.800 | 1.744 | -0.144 |
| L | SL | 11 | 0.509 | 1.049 | -0.030 |
| R | CH | 32 | 0.306 | 0.453 | 0.160 |
| R | CU | 27 | 0.567 | 0.861 | 0.273 |
| R | FC | 32 | 0.323 | 0.480 | 0.167 |
| R | FF | 165 | 0.365 | 0.452 | 0.278 |
| R | FS | 7 | 0.407 | 0.900 | -0.085 |
| R | FT | 70 | 0.381 | 0.504 | 0.259 |
| R | KC | 9 | 0.317 | 0.697 | -0.064 |
| R | KN | 8 | 0.088 | 0.294 | -0.119 |
| R | SI | 43 | 0.287 | 0.418 | 0.156 |
| R | SL | 103 | 0.267 | 0.360 | 0.174 |
左投手に対するwOBAが高いのは主にFFの結果が非常に優れていたことによる (平均wOBAとしてはSLも同じくらい高いが, サンプルサイズが小さいためwOBA全体への影響はFFより少ない; サンプルが足りていない球種では下限はwOBAが取りえない0より小さくなっているが, この辺はご愛嬌というところで). とはいえ, FFもNが小さく, 信頼区間が広いため下限は本人の全打席のwOBAと同程度というところ.
この例では都合よく (?) 信頼区間の幅の広さからこのようなデータの分割が, 推定という観点からは意味があまりなさそうなことがわかるが, このように状況を分けていくことは比較の多重性の問題 (multiple comparison probrem) を引き起こすので, そこも注意が必要である.
McCutchenが左に強いのかどうかというのであれば, データを増やせばより妥当性のある議論が可能になるだろう. とりあえず2015-2018にしてみる.
| 2015-2018 | |||||
|---|---|---|---|---|---|
| p_throws | PA | wOBA | CI_95 | CI_upper | CI_lower |
| Andrew McCutchen | |||||
| L | 627 | 0.390 | 0.044 | 0.433 | 0.346 |
| R | 2041 | 0.357 | 0.023 | 0.380 | 0.335 |
1年分に比べて効果はだいぶ小さくなった. 上で示したようにMLB全体での平均的な効果は0.014であった. これらの事前の情報よりはまだ大きいが, 2017単独に比べれば差は小さくなった. 信頼区間も狭くなったが, やはり区間同士には重なりがある. 実際にはこの場合は簡単に過去データを増やせるので, もっと増やすことでより意味のある判断が可能になるだろう.
この例では, 2017で大きな差があったがサンプルを増やしたことで効果は小さくなった. これは適当に選んで結果的にそうなったということで, いわゆるanecdotalな例である. もちろん常にそうなるわけではないが, 左右差の効果が平均的な打者とJudgeとの間の差より大きいというのはあまりありそうな話ではなく, 平均的な効果に向かって回帰したのは偶然ではないだろう.
The BOOKは左右の効果自体の回帰の大きさにも言及しており, 右打者の場合, 対左投手の打席が2200 (!) ほどで半分が回帰で失われるとしている (pp. 159).
つまりMcCutchenの2017年の結果に適用すると,
PA.batter <- 145
PA.even <- 2200
avg <- 0.014
lr.apparent.effect <- 0.47 - 0.342
avg + (lr.apparent.effect - avg) * (PA.batter / (PA.batter + PA.even))## [1] 0.02104904が回帰を考慮して推定された効果であり, 追加的な効果の大きさは下のようになる.
(lr.apparent.effect - avg) * (PA.batter / (PA.batter + PA.even))## [1] 0.007049041つまり, 平均的な効果 (0.014) にかなり近いところにまで推定される効果 (0.014 + 0.007) は低下した. 半分が回帰で失われる2200打席に比べて, McCuthenの2017の対左打席数は非常に小さく, 見かけ上の効果 (0.47 - 0.342 = 0.128) のほとんどが回帰で失われた.
95%信頼区間と回帰を使った議論を比較したい. 信頼区間では平均値を中心にして母平均がどの程度の範囲にありそうか, という形で不確実性を提示する. 対して, 回帰では観測値を平均方向へとシフトさせる. 回帰を利用した手法では, ある指標が持つ傾向を事前の情報として考慮することで不確実性に対処する. このため, このような平均値へと近づける効果が起こる. 点推定値による議論が可能なことは効果を量的に表現する上で大きな恩恵を持っている. この例で言えば, 対左投手が年間150 PA程度でwOBAで0.007なら, 野球においてそれほど重要な違いとは言えない, といった議論がありうるだろう. このような効果の大きさについての議論は, 95%信頼区間を利用した方法では不可能である.
2.3.2 得点圏
ついでにみんな大好き得点圏 (これはあくまで例のために使っているのであり, 実際にチャンスに強いかどうかを議論したい場合, Leverage Index (LI) などのより多くの要素を考慮して状況を量的に表現した指標を利用したほうが望ましい).
下はなんとなく目についた, 得点圏かどうかでwOBAに差がありそうだった2017年のCarlos Gonzalezの成績を示す.
data <- data %>%
mutate(RISP_FL = ifelse(is.na(on_2b) == FALSE | is.na(on_3b) == FALSE, 1, 0))
data %>%
filter(woba_denom == 1, game_year == "2017")%>%
filter(! grepl("bunt", des))%>%
filter(player_name == "Carlos Gonzalez")%>%
group_by(player_name, RISP_FL)%>%
dplyr::summarise(N = n(),
Avg = mean(woba_value),
SD = sd(woba_value),
SE = SD /sqrt(N),
t_crtcl = qt(0.975, df = N-1),
CI_95 = SE * t_crtcl,
CI_upper = Avg + CI_95,
CI_lower = Avg - CI_95)%>%
rename(wOBA = Avg, PA = N)%>%
mutate_if(is_numeric, funs(round), 3)%>%
gt()%>%
tab_header(title = "2017")| 2017 | ||||||||
|---|---|---|---|---|---|---|---|---|
| RISP_FL | PA | wOBA | SD | SE | t_crtcl | CI_95 | CI_upper | CI_lower |
| Carlos Gonzalez | ||||||||
| 0 | 389 | 0.371 | 0.535 | 0.027 | 1.966 | 0.053 | 0.424 | 0.318 |
| 1 | 141 | 0.241 | 0.429 | 0.036 | 1.977 | 0.071 | 0.313 | 0.170 |
得点圏 (RISP_FL = 1) ではかなりひどい成績になっている. 得点圏での打席数が少ないため信頼区間は広いが, 各状況の信頼区間はかろうじて重なっていない. NPBであれば, スポーツ新聞やなんJには怨嗟の声が満ちていたことだろう.
比較のためMLB全体の結果を計算する (MLB 15-18).
data %>%
filter(woba_denom == 1)%>%
filter(! grepl("bunt", des))%>%
group_by(RISP_FL)%>%
dplyr::summarise(N = n(),
Avg = mean(woba_value),
SD = sd(woba_value),
SE = SD /sqrt(N),
t_crtcl = qt(0.975, df = N-1),
CI_95 = SE * t_crtcl,
CI_upper = Avg + CI_95,
CI_lower = Avg - CI_95)%>%
rename(wOBA = Avg, PA = N)%>%
mutate_if(is_numeric, funs(round), 3)%>%
gt()%>%
tab_header(title = "MLB 2015-2018")| MLB 2015-2018 | ||||||||
|---|---|---|---|---|---|---|---|---|
| RISP_FL | PA | wOBA | SD | SE | t_crtcl | CI_95 | CI_upper | CI_lower |
| 0 | 548711 | 0.327 | 0.520 | 0.001 | 1.96 | 0.001 | 0.328 | 0.326 |
| 1 | 177150 | 0.331 | 0.511 | 0.001 | 1.96 | 0.002 | 0.334 | 0.329 |
全体で見ればほとんど違いが無いことがわかる.
[このパラグラフは本文と関係がない] 一応補足しておくと, 得点圏でわずかにwOBAが高いようにも見えるが, この結果だけから全体的に打者が勝負強い/投手が勝負弱い, などというような議論は困難である. まず, そもそも状況によって投手も打者も能力が同じではないと思われる. 例えば, 得点圏を多くつくる投手は, wOBAが高いだろう. また, 得点圏では1塁が空いている事が多くBBが出やすい (IBBは除かれているがIBB的な要素を持ったBBが発生しやすい; BBとIBBは, 0% IBBと100% IBBの間のどこかに位置する連続的なものと見たほうが良いだろう), あるいは守備位置が制限されるため安打が出やすくなりうる, などの小さい影響がありうる. さらに補足すると, 仮に得点圏などの状況で成績の違いらしきものが見られたとして, 多くの人が持つ概念としての“勝負強さ”と無関係に, そのスキルセットが上で挙げたような状況に対して相性の良い打者/投手が好成績を残している可能性もありうる. いずれにせよ状況に関して注目して個人の能力への影響を調べた場合, MLBでもNPBでもあっても, 影響は非常に小さいことは繰り返し示されている.
当然Carlos Gonzalezについて, サンプルを増やした時の結果を確認する. 2015-2018のCarlos Gonzalezの成績.
| 2015-2018 | |||||
|---|---|---|---|---|---|
| RISP_FL | PA | wOBA | CI_95 | CI_upper | CI_lower |
| Carlos Gonzalez | |||||
| 0 | 1668 | 0.360 | 0.026 | 0.386 | 0.334 |
| 1 | 584 | 0.344 | 0.046 | 0.391 | 0.298 |
差はかなり小さくなった. これもやはりanecdotalな例だが, 得点圏での成績については通算レベルでも全体的には偶然生じる程度の差とあまり変わらない, というよく知られた事前の情報を考えると, 回帰によって効果の大部分が失われたのは偶然と考えるべきではないだろう.
保守的な判断基準であるはずの信頼区間が重なっていない例なのにデータを増やしたら効果がほぼ失われた. これは, 2017のCarlos Gonzalezは得点圏別の状況で差が大きいものを恣意的に選んできたことに関係があるだろう. つまり, たくさんの信頼区間を持つデータの中から, 他の信頼区間 (つまり他の選手) を排除している. いわゆるdata dredging (データの峻別) 的な意味で, (1標本の) 信頼区間を用いた議論の前提が満たされていない (頻度論的な有意差検定でも, 比較の多重性を考慮しないと同じような問題が生じる). この例の場合, デモンストレーションとして意図的に問題のある状況を作りだしたわけだが, 問題を理解していないために意図せずにこのような議論を展開してもこのような平均への回帰的な現象は同様に生じる.
こちらでも2017のデータから回帰を使った点推定をざっくりしてみる. 得点圏ではないがclutchは6250 PA程度で半分が回帰するとされているので, これをそのまま近似として利用してみる.
PA.batter <- 141+389
PA.even <- 6250
avg <- 0.004
risp.apparent.effect <- 0.24 - 0.37
avg + (risp.apparent.effect - avg) * (PA.batter / (PA.batter + PA.even))## [1] -0.006474926一年分の打席数は得点圏の効果を推定する上ではあまり意味がない程度の打席数であるため, 測定された効果の大部分 (90%以上) は回帰によって失われて, 推定された効果は-0.006程度となった. 得点圏は得点価値や勝利価値の振れ幅が大きいため, この程度でも意味は多少あるかもしれない. また, これは1年で記録されたかなり極端な結果を恣意的に取り出した例であり, 打席数を増やして回帰の量を適用してみればより妥当な推定が可能になるだろう. また, 得点圏におけるwOBA差の回帰の大きさを実際に調べて, それを適用するべきだろう (そんな暇があればLIを使って調べた方がいいが).
3 まとめ
サンプルサイズが限定的な状況においては, 測定結果には偶然生じた誤差の影響が含まれている可能性が高い. このため, 実際のスキルがどの程度なのかを推測する場合には, 単純にデータから得られた平均ではなく, 母平均を推定する必要が生じる. 簡単な方法としては95%信頼区間が, Sabermetricsでよく使われる方法としては平均への回帰を利用する方法がある. もちろん, 仮にサンプルサイズがそれほど大きくない状況であっても, 得られた平均値は過去に起こったことを表現する数値としては妥当である. しかし, 得られた平均値は過去に起こったことは考慮する一方で, 過去に起こる可能性があったが起こらなかったことは考慮しないため, 実際に起こったことに過剰に適合している. そのため, 実際の能力を反映していない可能性や, 予測力を欠く可能性があることは注意が必要となる.
平均への回帰の計算では, ある打者の観測値の予測力に影響を与える要因としては下の要素がある:
その選手の結果の平均からの乖離の大きさ
指標自体に内在するばらつきの大きさ
打者の間の能力差
3つの要素のうち後者2つは冒頭の式で言えば\(PA_{even}\)であり, これは平均からの乖離の半分が打者の能力とみなせる打席数ということになる. この値は冒頭の式を介して, サンプル数が少ない状態においても比較的妥当な推定値を与えうる. 95% CIは1や3は考慮しない.
予測力を欠く指標としてよく知られる代表的なものは, 得点圏における成績だろう. 個人の1年単位で見れば極端な成績を残す選手は数多くいるが, 全体的に見ればこれは極めて弱い予測力しかないことが繰り返し示されている (蛭川, 2018). これは得点圏の打席自体が小さいために偶然の影響が大きく, また個人間の能力差がその影響に対して非常に小さいため, 個人の1年単位で非常に極端な成績であってもほとんどの場合偶然によって生じたものを見てしまっているためだと理解できる. このような状況に対して, 95% CIは測定値を中心として広がりの大きさで不確実性を提示する. これは測定から得られる情報のみから作られており個人間の能力差などは考慮しない. 対して, 平均への回帰は, その指標における回帰の大きさを使って偶然の影響を取り込み, その選手の実際の成績と妥当な折り合いをつけた数値を提示する.
多項分布からwOBAを取り出して, その母平均の95% CIが機能するかどうか検討した. wOBAの各打席ごとの結果についての分布はやや特殊だが, 30 PA程度のサンプルサイズで十分に信頼区間的な性質を持っているようだった. 得られた平均値からスキルを推定する際に, 生じている不確実性を理解する, あるいは読者に示すための方法として妥当性はありそうだ. 95%信頼区間は計算自体は非常に容易であり, 平均などを計算するついでに計算できる点が利点だろう.
実際には信頼区間を用いることで問題が起こりうると思われる. 多くの場合, 状況ごとの比較に興味があることが多いが, 信頼区間同士の比較はデータの示し方としては複雑となる. またその解釈にはある程度知識が必要で, さらに多重比較的な問題を含んでいるため, 野球記事の想定読者を考えると混乱を招きやすいかもしれない. 回帰を利用した方法は, 全体的な傾向を前提として考慮に入れるため多重比較的な問題を回避でき, 点推定的にデータを示すだけでもある程度正当性がある. このようなデータの提示は, 誤解が少ないという極めて良い利点があると思われる. また, この点推定で示せることによる利点として, 興味を持っている効果, つまり回帰すべき値との差, が野球においてどの程度の影響があるかについて, 価値への変換が容易となる点も挙げられるだろう (2.3.1 終わり当たり参照).
プロ野球の世界には一見すると偶然ではあり得ないような数値を示すような選手を探すことはそれほど難しくはない. この一つの原因は単純に多くの選手がいるためであり, 多重比較の問題として理解できる. 例えば, 300人に1回しか起こらないイベントも, 規定打席到達者が300人いれば期待値としては1人において起こる. このような稀な出来事に対して, 左右に関しての球種ごとの平均値を出したように, 状況を分割することでどこでその極端な数値が起こったかを調べ, それに基づいてそれらしい話をでっち上げることは可能だろう. しかし, これは打席数の減少による不確実性の増加を伴う. スキルの推定や将来の予測などの目的においては, 単年すべての結果の平均値をそのまま使うことですら妥当性はあまり高くはない. また比較する対象が増えて多重性の問題は深刻になる. 結果として, その分割された状況において測定された平均値から議論することは推定や予測に関してほとんど意味を持たない.
野球の統計において, 平均への回帰を表現する道具として最も一般的なものは, 年度間相関だろう. 年度間相関は計算が非常に容易であり, また回帰の大きさの悪くない桁感を与えてくれる. しかし, 年度間相関はサンプルサイズとの関係が不明瞭であるという欠陥を抱えている. 例えばある指標の年度間相関が0.5であり, 300 PAをサンプルサイズとして持っているある選手がいる状況を考える. この時, 測定値を平均へどれだけ回帰させれば適切な推定値になるかということは残念ながら計算できない. これは年度間相関の計算ではサンプルサイズの違いに起因する分散と, 能力差に起因する分散が分解されていないことによって起こる (サンプルサイズの違いの影響を小さくするために一定の打席数以上に絞るが, サンプルサイズの個人差は残っている). 一方, 半分が回帰するために必要なサンプルサイズを計算する場合では, サンプルサイズの影響が調整されるため, 冒頭で示した式と組み合わせることで妥当な回帰量を得ることができる.
回帰の大きさを利用した方法の問題の一つは, それぞれの効果に対して回帰の量を計算する必要性があることだろう. 回帰を考慮した効果の大きさを推定するのは, 概念的に単純とは言えないかもしれない. ランダムサンプリングなどを利用する場合, 概念的には単純だがやや面倒であるし, 技術的には少しはコーディングができる必要があるだろう. 平均的状況におけるwOBAの回帰の大きさを各状況に適用して, それを使って効果を計算することもできるかもしれないが, それぞれの効果自体の回帰の大きさを考慮した方法の結果と一致するとは限らない. また, MLBの数値をそのままNPBの選手に適用することも指標によっては不適切になりうる. ここではwOBAを扱ってきたが, NPBでのwOBAの回帰の大きさは, MLBよりもかなり小さいようだ (データは示していない; 主に長打力の個人差が大きいことが原因だと思われる). 興味を持っている個別の効果について1つづつ, 効果の回帰の大きさを評価していくほかないだろう.
4 参考文献
大久保街亜 岡田謙介, 2012, 伝えるための心理統計, 勁草書房.
John K. Kruschke, 2015, (邦訳: 前田 和寛, 小杉 考司 ほか, 2017), ベイズ統計モデリング: R, JAGS, Stanによるチュートリアル, 共立出版.
西原史暁, 2017, DeNA に対する第三者委員会の調査報告書での信頼区間の説明.
Steve Slowinski, 2010, Sample Size, Fangraphs.
Tango, MGL, and Dolphin, 2007, The BOOK, Potomac Books.
BellCurve, 2017, 平均への回帰、相関係数―統計学史(2)
蛭川 皓平 (baseballconcrete), 2013, タンゴの回帰式.
Sean Dolinar, 2015, A New Way to Look at Sample Size: Math Supplement, Fangraphs.
Tango, 2010, Regression equations for pitcher events
Russele Carleton, 2012, Baseball Therapy: It’s a Small Sample Size After All, BP.
Inverse-variance weighting (Wikipedia)
小池裕太, 2017, 統計データ解析 (I) 講義資料7
MINAKA Nobuhiro, 2013, Central Limit Theorem using R.
Statcast data are property of MLB advanced media.
Alex Reinhart (西原史暁 訳) ダメな統計学 (6) 有意であるかないかの違いが有意差でない場合, Original ver.
Johnson VE, 2013, Revised standards for statistical evidence, PNAS, 110, 48, pp. 19313-19317.
multiple comparison probrem (Wikipedia)
Tango, 2009, The color of clutch, The Hard Ball Times.
Cyril Morong, Clutch Hitting Links.
蛭川 皓平 (baseballconcrete), 勝負強さの研究 in プロ野球を統計学と客観分析で考える デルタ・ベースボール・リポート2, 2018.
5 Appendix
5.1 出塁率のベイズ推定と, その点推定値と回帰による推定値の比較
冒頭で示した100 PA, 500 PAにおける出塁率0.4のデータが得られた場合について, 事前の情報がある場合の推定を行う. この場合は非常に単純なケースなので, 解析的に取り扱うことができる. 計算方法の詳細については説明しないで結果だけ示す. 参考文献のJohn K. Kruschke, 2015 (邦訳2017) のchapter 6にデータがベルヌーイ分布に従う場合の詳しい説明があるので, そちらなどを参照していただきたい.
事前分布にはベータ分布を使い, 平均は0.33とし, SDはThe BOOKで例として示された出塁率のSD(talent), 0.025に近づける. このときの\(PA_{even}\)は350となる (pp. 380 in The BOOK). ベータ分布における, スケールパラメータa, bと平均値やSDとの関係から, 事前分布のパラメータaとbの値を以下のように設定することにする.
# mean = a /(a+b) なので
# aを33の倍数, bを67の倍数で適当に入れていく
a <- 33 * 3.5
b <- 67 * 3.5
mean <- a /(a+b)
kappa <- a+b
Var <- (a * b)/ (kappa^2 * (kappa+1))
SD <- sqrt(Var)このときの事前分布の平均とSDを示す.
mean## [1] 0.33SD## [1] 0.0250981これらのパラメータを持つbeta分布をプロットする. 横軸がOBPの推定値\(\theta\)を, 縦軸は確率密度を示す.
plot.beta<- function(a, b){
y.txt <- paste0("dbeta(", a," ,", b, ")")
plot(function(x) dbeta(x, a, b),
0.2, 0.5,
xlab = "theta",
ylab = y.txt)
}
plot.beta(a, b)
この事前分布は0.33あたりにピークを持ち, 0.4を超えることはほとんどない. 上位1%をクリアするための数値を確認する.
qbeta(0.01, a, b, lower.tail = FALSE)## [1] 0.3897089これはイメージとしては, 通算成績レベルではOBP 0.39や0.4をクリアする選手はほとんどいない, という実際に観測される全体の傾向, あるいは野球を見ている人間の常識や感覚, のようなものを反映している (実際には偶然の影響で起こるようなばらつきが取り除かれた状態のSDを考えている).
この事前分布のもとで, 100, 500, 7000 PAにおいてOBP = 0.4が得られたときの, 事後のOBPの推定値\(\theta\)をプロットする.
a <- 33 * 3.5
b <- 67 * 3.5
par(mfrow=c(4,1))
plot.beta(a, b)
plot.beta(a+40, b+100-40)
plot.beta(a+200, b+500-200)
plot.beta(a+2800, b+7000-2800)
上から事前分布, 100 PA, 500 PA, 7000 PAにおける, 事後の推定結果を示す. 事前の情報がない場合, つまり尤度関数のみから推定した場合ピークは0.4になるが, 100 PA, 500 PAともに0.4よりはかなり低い. 事前の情報が事後分布を全体的に事前分布へと近づけている. データ (尤度関数) への重み付けは\(\frac{N}{N+a+b}\)であり, Nは増えるほど尤度関数へと近づいていく. これによって, 一番下に示した例のようにPAが非常に大きくなれば, その測定値付近の狭い範囲にあると推定されることになる.
ベータ分布\(beta(a, b)\)の平均値は\(\frac{a}{a+b}\)であるから, 500 PAを観測したあと平均値の推定値は以下で求められる.
(a + 200) /(a + 200 + b + 300)## [1] 0.3711765あるいは最頻値\(\omega\)は\(\frac{a-1}{a+b-2}\)であるから, 500 PAを観測したあとの最大事後確率 (MAP) の推定値は以下で求められる.
(a + 200-1) /(a + 200 + b + 300-2)## [1] 0.3708726平均への回帰からの点推定値と比較する. 500 PAの場合では以下のようになる.
PA.batter <- 500
PA.even <- 350
0.33 + (0.4 - 0.33) * (PA.batter / (PA.batter + PA.even))## [1] 0.3711765ベイズで得られた平均値やMAP推定値とほぼ同じ点推定値を提示していることがわかる.
いずれの方法も特定の選手について成績を推定する際に, その選手のデータに加えて事前の事実や常識のようなものを考慮に入れて推定する. 比較的少ない打席数でOBP 0.4という極端な数値が得られた場合で言えば, この極端な測定値と事前の情報の間の乖離が大きく, それらの間の折り合いをつけるようなイメージになる. 重要なことは, 事前の情報 (事前分布) と, その選手の成績とサンプルサイズを考慮して, “量的に適切な”折り合いをつけようとしていること. 全体の傾向を無視してその選手の成績をそのまま受け入れる, あるいは反対にサンプルサイズが小さいため信用できないとデータを全て捨てる, などといった二分法に比べて, より多くの情報を考慮して結論を得ることが可能になる.
5.2 回帰の大きさを示す式の導出
\[Estimated\ talent = Average + (Observed\ value - Average) * \frac{PA_{batter}}{PA_{batter} + PA_{even}}\]
この式はすでに説明したように, 特定の打者について回帰の大きさを考慮した推定値を得るために用いることができる.
異なる測定結果を統合する方法として, 分散の逆数で重み付けする方法があり, これが上の式を導出する方法の1つとなる. 分散が大きいほど, つまり不確実性が大きいほど, 分散の逆数は小さくなり, 弱い重み付けを持つことになる. 例えば, ある選手に関して2つの出塁率 (OBP) の測定結果がある状況を考える. 測定結果iの成績を\(OBP_i\), その分散を\(Var(OBP_i)\)としてこの成績を統合して, 2つの測定結果が反映された値\(\hat{OBP}\)を得るためには以下のような計算を行えば良い.
\[\hat{OBP} = \frac{\frac{OBP_1}{Var(OBP_1)} + \frac{OBP_2}{Var(OBP_2)}}{\frac{1}{Var(OBP_1)} + \frac{1}{Var(OBP_2)}}\]
The BOOKの具体例を引用する. 仮に全体の平均を0.33とし, SD(talent)が0.025とする. 100 PAにおいて0.45の出塁率を残した選手がいたとする. 全体の能力差の情報にこの追加的な情報を加えて, この選手の出塁率における能力を推定したい. この打者の成績の分散として, 二項分布から期待される大きさを利用して, 分散の逆数を重み付けに使って平均を取ると以下のようになる.
\[\begin{aligned} Estimated\ talent &= \frac{Avg/Var(talent) + Observed\ value/Var(luck)}{1/Var(talent) + 1/Var(luck)}\\ &= \frac{0.33/0.025^2 + 0.45/0.05^2}{1/0.025^2 + 1/0.05^2} \end{aligned}\]
ここで二項分布の公式から\(SD(luck) = \sqrt{0.45*(1-0.45)/100}\), つまり:
sqrt(0.45*(1-0.45)/100)## [1] 0.04974937と\(Var = SD^2\)を利用していることに注意. これを計算すると推定値は以下のようになる.
num <- 0.33/0.025^2 + 0.45/0.05^2
denom <- 1/0.025^2 + 1/0.05^2
num/denom## [1] 0.354基本的な操作としてはこの逆分散を利用した重み付けだけなのだが, 上の式は見通しが悪いのでもう少し変形させる. 平均値からの差を測定する状況を考える. この場合, 無情報における推定値は平均が0であり, SD(talent)やSD(luck)は平均値を考えた場合と変わらない (単に分布をシフトさせただけであるため). ある特定の選手を考えると:
\[\begin{aligned} Estimated\ talent - Avg &= \frac{0/Var(talent) + (Observed\ value - Avg)/Var(luck)}{1/Var(talent) + 1/Var(luck)}\\ &= (Observed\ value - Avg) *\frac{1/Var(luck)}{1/Var(talent) + 1/Var(luck)} \end{aligned}\]
ここでOBPの場合を考えると, \(Var(luck) = 0.33 * (1-0.33) /N ≈ 0.22/PA_{batter}\)で概ね近似できる.
0.33 * (1-0.33)## [1] 0.2211これを使い, さらにAvgを移項して上の式を変形すると以下のようになる.
\[\begin{aligned} Estimated\ talent &= Avg + (Observed\ value - Avg) *\frac{1/Var(luck)}{1/Var(talent) + 1/Var(luck)}\\ &≈ Avg + (Observed\ value - Avg) * \frac{PA_{batter}/0.22}{PA_{batter}/0.22 + 1/Var(talent)}\\ &= Avg + (Observed\ value - Avg) * \frac{PA_{batter}}{PA_{batter} + 0.22/Var(talent)}\\ &≈ Avg + (Observed\ value - Avg) * \frac{PA_{batter}}{PA_{batter} + 350} \end{aligned}\]
最後の350はこの場合のOBPの\(PA_{even}\)であり, ここでは以下によって求めている.
0.22 / (0.025^2) # 0.22 / sd(talnet)^2## [1] 352これにより:
- 全体の平均値
- 特定の選手の成績
- その打席数
という馴染みのある測定値を利用して (実際には\(PA_{even}\)を得るためにその指標の分布の性質が利用されている), さらに右辺に:
- 平均値
- 回帰を考慮しても残った効果
の2つの要素を持つ, 推定値を得るための見通しの良い式が得られた.