野球っぽいデータで多重比較に伴う問題を眺める
Author: @sleep_in_nmbrs
2019-10-12
1 背景と目的
このポストでp値の多重性の補正について言及した. この点について, ランダムに生成したデータで具体的に適用してみる. 多重性の問題は対処法はp値で特に広く使われているが, 実際の問題としてはp値に限ったことではなく, 頻度論的な方法でよく見られる問題である.
p値は
- 仮に特定の統計モデルのもとで帰無仮説が正しい (ざっくりいうと比較したいものと差がない) と仮定した場合に,
検定方法の前提 (標本分布, データの取得方法/データ収集停止規則 etc.) が完全に満たされている
非実在状況において,得られた結果と同じかそれ以上に起こりにくい結果になる確率
を示す. 典型的には何らかの有意性基準, 例えば0.05を定めて, これよりも小さい場合, つまり起こりにくいことが起こった場合に有意であると考えることが多い. 野球のデータ分析で用いられることは少ないが, 例えばFiveThirtyEightのこの記事で示されていたり, あるいはdeltagraphsのこの記事でも示されている.
p値はそれ自体の意味も正確に理解されにくいが, さらに多重に比較することで解釈はより誤解を生みやすくなる. 例えば, 打者の得点圏とそれ以外について出塁率を比較して, その値に違いがあるかをp値で判断する場合を考える.
仮に世界に打者が１人であり, この打者の得点圏とそうでない場合の出塁能力が同一である場合を考える. この場合, 少なくとも一人の打者が偶然にpが0.05以下となるような違いを持つ確率は,
1 - 0.95## [1] 0.05で5%である. 十分に低いと言えるかもしれない (思えないならもっと基準を厳しくすればいいだろう).
次に世界に300人の打者が存在し, それらの打者の得点圏とそうでない場合の出塁能力が同一である場合を考える. この場合, 期待値としては15人が偶然によって生じた誤差によって, この基準をクリアする.
300 * 0.05## [1] 15このようなp値をクリアする人数の期待値と, 現実のデータを比較することも古くから行われている. Pete Palmerは得点圏打率 (pdf; pp.6) でTom Tangoは投手BABIPなどで, ほぼ同じものである (p値へと直接変換可能な) 各選手のZ scoreを用いてこのような観点から結果を評価している.
あるいは, 同じく300人の世界で少なくとも一人の打者が偶然にpが0.05以下となる確率は,
1 - 0.95^300## [1] 0.9999998で99.99998%となり, ほぼ間違いなくp値が0.05以下をクリアする選手が1人は存在することになる.
pの値が同一であっても, 比較する対象 (この場合では選手数) が多くなると, 偶然にp値の基準をクリアする確率が上がってしまうため, 選手が1人でp値の基準をクリアする場合に比べてその意味は大きく異なる. これは多重比較の問題として古くから知られており, 様々な対処法が考えられている. 近年ではその問題は広く知られるようになっており, 少なくとも学術誌では, 多重性について何らかの補正がなされたp値が報告されることが多くなっている.
1.1 目的
ここでは乱数から結果を取得し, 得点圏とそうでない場合の出塁率っぽい数値を生成する. このデータを用いて, p値と, 補正されたp値を計算し, それらを判断基準として使った場合に, どの程度正しい, あるいは間違った結論が出されるのかを眺める. p値を使った有意性判定を使って調べていくが, 比較する仮説の数が増えればどこかで偶然極端な結果が得られる確率が上がることは, p値以外の仮説を設定する方法でも同じことである.
一応書いておくと, p値を使うことが望ましいとかそういうような話ではない. あくまで, 多重性の問題がどういうものであり, それを補正するという行為がどういうことかを眺めるのが目的である. 上で引用した野球の例はそれぞれなにがしかの指標が持つ効果を推定しようとする試みだが, 近年ではこのような場合, 指標の信頼性という形で整理されている. 例えば, 勝負強さについてであれば6000 PA程度で半分が本人のタレントと推定できるといった具合である. このような方法は全体の傾向を取り込んでいるため, 多重性の問題が生じにくい. ちなみに学術系の業界でも, 多くの統計学者がp値のみに基づく判断は大きな問題を引き起こしていると批判を強めており (ASA声明翻訳版pdf), 某有力週刊誌にコメントとして掲載されていたりする (Amrhein et al., 2019).
ちなみに中の人は統計がよくわかっていないので, 以下の内容が正しいかどうかは各自の責任で判断をお願いしたい.
2 結果
2.1 状況に関して差がない場合
各選手について, 異なる2つの状況の成績を比較したい. ここではその2つの状況において, 成績には実際には差がないというシナリオをランダムサンプリングでモデル化する.
二項分布から300人分のデータを2つの異なるサンプルサイズ (打席数) をもつ状況について取得する. 便宜的に2つの状況は, 得点圏 (RISP) とそれ以外 (notRISP) とし, 適当にそれぞれ100, 500打席分とする. 出塁率0.33として出塁率として考えることにする. 人数をサンプル数, 打席数をサンプルサイズと呼ぶ. サンプル数の数だけ検定が繰り返されることになる.
require(tidyverse)
require(gt)
# 疑似乱数の種を設定
set.seed(2019)
# 打席数の設定
ab.RISP <- 100
ab.not.RISP <- 500
# 300人分データを二項分布で取り出す
# もっと多くしたいところだが, ここでは2つの理由で現実に近づけておく
# 1. 300人程度でも偶然に有意性の基準を超えるものが十分に現れることを示す
# 2. ここでサンプル数を増やすとあとで示す多重性の補正が現実に比べて強くなりすぎる
player.num <- 300
data <- data.frame(player = 1:player.num,
notRISP = rbinom(n = player.num, size = ab.not.RISP, prob = 0.33),
RISP = rbinom(n = player.num, size = ab.RISP , prob = 0.33))p値の計算には, 2群の比率の差の検定を用いる. 有意水準αは慣習に従って0.05とする (気に食わないなら好きに変えれば良い).
各仮想選手ごとに成績をまとめ, p値が低い順に示す.
# 頻度の検定を行いp値を返すための関数
return.p <- function(x){
p <- prop.test(c(x[1], x[2]), c(ab.not.RISP, ab.RISP ), correct=F)[3]
p <- as.numeric(p)
return(p)
}
p.values <- data%>%select(-1)%>%
apply(1, return.p) # return.pを各行に対して適用
data.prop.test <- data %>%
mutate(p = p.values,
notRISP = notRISP /ab.not.RISP, # rateに変換
RISP = RISP /ab.RISP )
data.prop.test %>%
arrange(p)%>%
head(10)%>%
mutate_if(is.numeric, funs(round), 3)%>%
gt()| player | notRISP | RISP | p |
|---|---|---|---|
| 24 | 0.318 | 0.47 | 0.003 |
| 110 | 0.324 | 0.18 | 0.004 |
| 79 | 0.360 | 0.22 | 0.007 |
| 283 | 0.302 | 0.44 | 0.007 |
| 212 | 0.324 | 0.46 | 0.009 |
| 127 | 0.298 | 0.43 | 0.010 |
| 19 | 0.348 | 0.22 | 0.013 |
| 264 | 0.294 | 0.42 | 0.013 |
| 3 | 0.346 | 0.22 | 0.014 |
| 275 | 0.308 | 0.43 | 0.018 |
この集団はRISPとnotRISPで十分に大きいサンプルサイズでは同じ成績になるように取り出しているにもかかわらず, pが0.05以下となったデータが予想通りそれなりに存在することがわかる.
pが0.05以下となった頻度を計算する. ここでは差がない状況を考えているので, 理想的には0.05 (5%) になるはずである.
mean(data.prop.test$p <= 0.05)## [1] 0.056666670.05よりも少し大きくなったが, サンプルを増やせば0.05に近づいた (本稿末尾のAppendix 1参照). ここで計算されたp値に基づいて判断した場合, 状況によって差がない仮想選手の一部に対して差があるという判断を下すことになる.
p値の分布を示す.
data.prop.test%>%
ggplot(aes(x = p)) +
geom_histogram(alpha = 0.7,
binwidth = 0.02,
colour="black", fill="gray")+
theme_bw(base_family = "HiraKakuPro-W3")+
theme(axis.text.x = element_text(size=10),
axis.text.y = element_text(size=10),
axis.title = element_text(size = 12)) +
labs(title = "p値のヒストグラム")
実際の数値を散布図で示す. 色と形でpが0.05以下であるかどうかを示す.
data.prop.test%>%
mutate(p = ifelse(p <= 0.05, "<= 0.05", "> 0.05"))%>%
ggplot(aes(x = notRISP,
y = RISP,
color = p,
shape = p)) +
geom_point(size = 1)+
coord_equal() +
theme_bw(base_family = "HiraKakuPro-W3")+
theme(axis.text.x = element_text(size=10),
axis.text.y = element_text(size=10),
axis.title = element_text(size = 12)) +
labs(x = "非得点圏 (notRISP)", y = "得点圏 (RISP)")
図の縦軸と横軸のスケールは一致させている. RISPではサンプルサイズが小さいため, 結果のばらつきが大きい.
p値を補正する. p.adjust()を用いて, Bonferroniの方法とHolmの方法による補正値を計算する. これらの方法は広いデータに対して適用可能な, 代表的な多重比較法である.
data.prop.test <- data.prop.test %>%
mutate(p_bonferroni = p.adjust(p, "bonferroni"),
p_holm = p.adjust(p, "holm"))
data.prop.test %>%
arrange(p)%>%
head(10)%>%
mutate_if(is.numeric, funs(round), 3)%>%
gt()| player | notRISP | RISP | p | p_bonferroni | p_holm |
|---|---|---|---|---|---|
| 24 | 0.318 | 0.47 | 0.003 | 1 | 1 |
| 110 | 0.324 | 0.18 | 0.004 | 1 | 1 |
| 79 | 0.360 | 0.22 | 0.007 | 1 | 1 |
| 283 | 0.302 | 0.44 | 0.007 | 1 | 1 |
| 212 | 0.324 | 0.46 | 0.009 | 1 | 1 |
| 127 | 0.298 | 0.43 | 0.010 | 1 | 1 |
| 19 | 0.348 | 0.22 | 0.013 | 1 | 1 |
| 264 | 0.294 | 0.42 | 0.013 | 1 | 1 |
| 3 | 0.346 | 0.22 | 0.014 | 1 | 1 |
| 275 | 0.308 | 0.43 | 0.018 | 1 | 1 |
Bonferroni法による補正値はp_bonferroni, Holm法による補正値はp_holmである. いずれの方法でも, 差があると判定された選手はいないという正しい結論となった.
Bonferroni法ではすべてのp値に仮説の数 k を掛け合わせて補正する. つまりこの場合, 選手の数 (k = 300) を掛ける. Holm法も似たような補正をするのだが, p値が小さい順からk, k-1, k-2, …を掛けるためBonferroniよりもやや補正が弱い. 理屈に興味がある方は教科書, 例えば永田と吉田, 統計的多重比較法の基礎あたりを参照していただきたい.
2.2 状況に関して差がある場合
2.2.1 現実的なサンプル数 (選手数) の場合
次に, 上で使用したデータの1割 (30人分) をRISPとnotRISPで成績に差がある選手 (clutchタイプと呼ぶ) に入れ替えて, 同様に計算していく. これは, 状況別成績に差がある選手を見つけられるかどうかを調べるためである.
- 差があるものを差があると判定すること
- 差がないものを差がないと判定すること
これらいずれも正しく判定することであるが, 問題として異なることに注意.
効果 (差) の大きさは0.04とした. これは決して現実においてそのようなものがありうると言っているのではなく, 単にそうだった場合にどうなるかを調べるというだけである.
打者のタイプ (clutch vs non_clutch) ごとに, 成績の平均やSDを示す.
RISP.effect <- 0.04
player.num2 <- 30
data2 <- data.frame(player = 1:player.num2,
notRISP = rbinom(n = player.num2, size = ab.not.RISP, prob = 0.33),
RISP = rbinom(n = player.num2, size = ab.RISP , prob = 0.33 + RISP.effect))
data2 <- data%>%
mutate(is_clutch_hitter = "non_clutch")%>%
slice(1:270)%>%
bind_rows(data2 %>%
mutate(is_clutch_hitter = "clutch"))
p.values2 <- data2%>%select(-c(1,4))%>%
apply(1, return.p)
data.prop.test2 <- data2 %>%
mutate(p = p.values2,
notRISP = notRISP /ab.not.RISP,
RISP = RISP /ab.RISP )
data.prop.test2%>%
group_by(is_clutch_hitter)%>%
dplyr::summarise(N = n(),
Mean_notRISP = mean(notRISP),
Mean_RISP = mean(RISP),
SD_notRISP = sd(notRISP),
SD_RISP = sd(RISP))%>%
mutate_if(is.numeric, funs(round), 4)%>%
gt()| is_clutch_hitter | N | Mean_notRISP | Mean_RISP | SD_notRISP | SD_RISP |
|---|---|---|---|---|---|
| clutch | 30 | 0.3327 | 0.3680 | 0.0214 | 0.0472 |
| non_clutch | 270 | 0.3302 | 0.3264 | 0.0192 | 0.0489 |
差がある選手 (clutch) ではRISPがで0.035ほど高くなっていたようだ. 0.04より少し低くなっているのは, 30人だけなので設定した効果からの誤差がある程度残っている結果だろう.
補正されていないp値が0.05以下となったデータの頻度を示す.
data.prop.test2%>%
group_by(is_clutch_hitter)%>%
summarise(N = sum(p <= 0.05))%>%
mutate(freq = N / sum(N))%>%
mutate_if(is.numeric, funs(round), 3)%>%
gt()| is_clutch_hitter | N | freq |
|---|---|---|
| clutch | 4 | 0.235 |
| non_clutch | 13 | 0.765 |
freqはpが0.05以下となったデータに占めるそれぞれの比率示している. 差がある選手は全体では10%だが, この例ではRISPにおいて能力差があるためにclutchタイプの濃縮がかかり23%程度まで増加した.
実際の数値を散布図で示す. 色と形でpが0.05以下であるかどうかを示す.
data.prop.test2%>%
mutate(p = ifelse(p <= 0.05, "<= 0.05", "> 0.05"))%>%
ggplot(aes(x = notRISP,
y = RISP,
color = p,
shape = p)) +
geom_point(size = 1)+
facet_wrap(~ is_clutch_hitter)+
coord_equal() +
theme_bw(base_family = "HiraKakuPro-W3")+
theme(axis.text.x = element_text(size=10),
axis.text.y = element_text(size=10),
axis.title = element_text(size = 12)) +
labs(x = "非得点圏 (notRISP)", y = "得点圏 (RISP)")
clutchタイプはRISPで全体的に数値が高くなっていることが確認できる.
同様に補正されたp値を計算する.
data.prop.test2 <- data.prop.test2 %>%
mutate(p_bonferroni = p.adjust(p, "bonferroni"),
p_holm = p.adjust(p, "holm"))
data.prop.test2 %>%
arrange(p)%>%
head(10)%>%
mutate_if(is.numeric, funs(round), 3)%>%
gt()| player | notRISP | RISP | is_clutch_hitter | p | p_bonferroni | p_holm |
|---|---|---|---|---|---|---|
| 8 | 0.288 | 0.44 | clutch | 0.003 | 0.833 | 0.833 |
| 24 | 0.318 | 0.47 | non_clutch | 0.003 | 1.000 | 1.000 |
| 110 | 0.324 | 0.18 | non_clutch | 0.004 | 1.000 | 1.000 |
| 79 | 0.360 | 0.22 | non_clutch | 0.007 | 1.000 | 1.000 |
| 212 | 0.324 | 0.46 | non_clutch | 0.009 | 1.000 | 1.000 |
| 127 | 0.298 | 0.43 | non_clutch | 0.010 | 1.000 | 1.000 |
| 7 | 0.318 | 0.45 | clutch | 0.011 | 1.000 | 1.000 |
| 19 | 0.348 | 0.22 | non_clutch | 0.013 | 1.000 | 1.000 |
| 264 | 0.294 | 0.42 | non_clutch | 0.013 | 1.000 | 1.000 |
| 3 | 0.346 | 0.22 | non_clutch | 0.014 | 1.000 | 1.000 |
補正されたp値で0.05以下をクリアしたデータは存在しなかった. 前項では, これらの補正を利用することで, 差がないものを差がないと判断する正しい結論を出すことはできた. 一方, 差があるものを見つけ出すことはできなかった.
しかし, これは必ずしも方法の問題であるとは限らない. そもそも補正されていないp値を見ても上位はほとんど差がないタイプのデータ (non_clutch) となっている. これは, 設定された差の大きさ (集団の10%の選手に+0.04の効果) を見いだせるほどのサンプルサイズ (打席数) がそもそも与えられていない可能性を示唆している.
そこで効果の大きさをそのままにして, サンプルサイズを大きくする. NPBの国内FA取得年に合わせて8倍にする.
# 打席数を増やして検出力を上げる
ab.RISP <- 800
ab.not.RISP <- 4000
n1 <- 270
n2 <- 30
data <- data.frame(player = 1:n1,
notRISP = rbinom(n = n1, size = ab.not.RISP, prob = 0.33),
RISP = rbinom(n = n1, size = ab.RISP , prob = 0.33))
data2 <- data.frame(player = 1:n2,
notRISP = rbinom(n = n2, size = ab.not.RISP, prob = 0.33),
RISP = rbinom(n = n2, size = ab.RISP , prob = 0.33 + 0.04))
data2 <- data%>%
mutate(is_clutch_hitter = "non_clutch")%>%
bind_rows(data2 %>%
mutate(is_clutch_hitter = "clutch"))
data.prop.test2 <- data2 %>%
mutate(p = apply(data2%>%select(- c(1,4)),
1, # 各行に対して適用
return.p),
notRISP = notRISP /ab.not.RISP,
RISP = RISP /ab.RISP)
data.prop.test2%>%
group_by(is_clutch_hitter)%>%
dplyr::summarise(N = n(),
Mean_notRISP = mean(notRISP),
Mean_RISP = mean(RISP),
SD_notRISP = sd(notRISP),
SD_RISP = sd(RISP))%>%
mutate_if(is.numeric, funs(round), 3)%>%
gt()| is_clutch_hitter | N | Mean_notRISP | Mean_RISP | SD_notRISP | SD_RISP |
|---|---|---|---|---|---|
| clutch | 30 | 0.333 | 0.372 | 0.009 | 0.016 |
| non_clutch | 270 | 0.330 | 0.331 | 0.007 | 0.017 |
pが0.05以下となったデータの頻度を示す.
data.prop.test2%>%
group_by(is_clutch_hitter)%>%
summarise(N = sum(p <= 0.05))%>%
mutate(freq = N / sum(N))%>%
mutate_if(is.numeric, funs(round), 3)%>%
gt()| is_clutch_hitter | N | freq |
|---|---|---|
| clutch | 18 | 0.45 |
| non_clutch | 22 | 0.55 |
実際の数値を散布図で示す. 色と形でpが0.05以下であるかどうかを示す.
data.prop.test2%>%
mutate(p = ifelse(p <= 0.05, "<= 0.05", "> 0.05"))%>%
ggplot(aes(x = notRISP,
y = RISP,
color = p,
shape = p)) +
geom_point(size = 1)+
facet_wrap(~ is_clutch_hitter)+
coord_equal() +
theme_bw(base_family = "HiraKakuPro-W3")+
theme(axis.text.x = element_text(size=10),
axis.text.y = element_text(size=10),
axis.title = element_text(size = 12)) +
labs(x = "非得点圏", y = "得点圏")
補正なしの状態において, 差のある選手については, 半数以上において差があると判定されるようになった.
ここではnon_clutchは状況に関連なく0.33の確率で出塁するように設定されている. しかし, これだけ打席を増やした状態で, 非得点圏 (4000 PA) ですら偶然生じる誤差がそれなりに残っている. もちろんより打席数が少ない得点圏 (800 PA) では誤差はさらに大きい.
同様に補正されたp値を計算する.
data.prop.test2 <- data.prop.test2 %>%
mutate(p_bonferroni = p.adjust(p, "bonferroni"),
p_holm = p.adjust(p, "holm"),
p_BH = p.adjust(p, "BH"))
data.prop.test2 %>%
arrange(p)%>%
head(10)%>%
mutate_if(is.numeric, funs(round), 3)%>%
gt()| player | notRISP | RISP | is_clutch_hitter | p | p_bonferroni | p_holm | p_BH |
|---|---|---|---|---|---|---|---|
| 20 | 0.323 | 0.401 | clutch | 0.000 | 0.005 | 0.005 | 0.004 |
| 17 | 0.332 | 0.410 | clutch | 0.000 | 0.008 | 0.008 | 0.004 |
| 28 | 0.325 | 0.399 | clutch | 0.000 | 0.017 | 0.017 | 0.006 |
| 155 | 0.350 | 0.282 | non_clutch | 0.000 | 0.063 | 0.062 | 0.016 |
| 16 | 0.319 | 0.385 | clutch | 0.000 | 0.092 | 0.091 | 0.018 |
| 1 | 0.312 | 0.375 | clutch | 0.001 | 0.159 | 0.156 | 0.026 |
| 22 | 0.327 | 0.389 | clutch | 0.001 | 0.224 | 0.219 | 0.032 |
| 51 | 0.322 | 0.378 | non_clutch | 0.002 | 0.733 | 0.715 | 0.092 |
| 19 | 0.328 | 0.380 | clutch | 0.005 | 1.000 | 1.000 | 0.164 |
| 173 | 0.320 | 0.370 | non_clutch | 0.006 | 1.000 | 1.000 | 0.187 |
サンプルサイズを大きくしたことによって, 補正なしのp値は差があるタイプ (clutch) が占めるようになった. しかし , Bonferroni (p_bonferroni) でもHolm (p_holm) でも, 0.05基準をクリアした選手はこの例では3人だけだった. 実のところ, ここで示した方法は多重比較法の中でも広く使えるが補正が強すぎることでお馴染みである. 教科書を参照すれば, 特定のデータタイプに対して利用可能な感度がより高い他の方法も紹介されているだろう.
その他の補正方法として, ここではBH法の結果も示した. 補正した値を比較すると, p_BHはBonferroniやHolmよりも低めに出ている (補正の強さが弱い). この方法は偽発見率 (False Discovery Rate: FDR) の調整を目的とする方法であり, 例えばBH法でpが0.05以下のデータを集めれば, FDRが0.05となる, ということである. この例で言えばFDRは
- 基準をクリアした選手の中で,
- 差がない選手 / (差がある選手 + 差がない選手)
から計算される.
p_BHが0.1以下のデータで頻度を集計する.
compute.FDR <- function(df,
Q = 0.1){
df%>%
filter(p_BH <= Q)%>%
group_by(is_clutch_hitter)%>%
dplyr::summarise(N = n())%>%
mutate(freq = N / sum(N))
}
data.prop.test2%>%
compute.FDR(Q = 0.1)%>%
mutate_if(is.numeric, funs(round), 3)%>%
gt()| is_clutch_hitter | N | freq |
|---|---|---|
| clutch | 6 | 0.75 |
| non_clutch | 2 | 0.25 |
FDRは0.25となり0.1とは多少乖離がある. そもそもサンプル数が少ない (差がない選手データは2つ) ため信頼できる結果であるかは疑問である.
2.2.2 非現実的に多いサンプル数 (選手数) でFDRを調べる
FDRがこの例でうまく働いているのかどうか調べるために, サンプル数を増やす. 状況による能力の有無の比率は変化させず, 差がないデータを90000, 差があるデータを10000とする.
実際の数値を散布図で示す. 色と形でpが0.05以下であるかどうかを示す.
ab.RISP <- 800
ab.not.RISP <- 4000
n1 <- 90000
n2 <- 10000
data <- data.frame(player = 1:n1,
notRISP = rbinom(n = n1, size = ab.not.RISP, prob = 0.33),
RISP = rbinom(n = n1, size = ab.RISP , prob = 0.33))
data2 <- data.frame(player = 1:n2,
notRISP = rbinom(n = n2, size = ab.not.RISP, prob = 0.33),
RISP = rbinom(n = n2, size = ab.RISP , prob = 0.33 + 0.04))
data2 <- data%>%
mutate(is_clutch_hitter = "non_clutch")%>%
bind_rows(data2 %>%
mutate(is_clutch_hitter = "clutch"))
data.prop.test2 <- data2 %>%
mutate(p = apply(data2%>%select(- c(1,4)),
1, # 各行に対して適用
return.p),
notRISP = notRISP /ab.not.RISP,
RISP = RISP /ab.RISP)data.prop.test2%>%
sample_n(10000)%>%
mutate(p = ifelse(p <= 0.05, "<= 0.05", "> 0.05"))%>%
ggplot(aes(x = notRISP,
y = RISP,
color = p,
shape = p)) +
geom_point(size = 1)+
facet_wrap(~ is_clutch_hitter)+
coord_equal() +
theme_bw(base_family = "HiraKakuPro-W3")+
theme(axis.text.x = element_text(size=10),
axis.text.y = element_text(size=10),
axis.title = element_text(size = 12)) +
labs(x = "非得点圏", y = "得点圏")
同様に補正されたp値を計算する.
data.prop.test2 <- data.prop.test2 %>%
mutate(p_bonferroni = p.adjust(p, "bonferroni"),
p_holm = p.adjust(p, "holm"),
p_BH = p.adjust(p, "BH"))
data.prop.test2 %>%
arrange(p)%>%
head(10)%>%
mutate_if(is.numeric, funs(round), 3)%>%
gt()| player | notRISP | RISP | is_clutch_hitter | p | p_bonferroni | p_holm | p_BH |
|---|---|---|---|---|---|---|---|
| 6493 | 0.320 | 0.429 | clutch | 0 | 0.000 | 0.000 | 0.000 |
| 6962 | 0.320 | 0.425 | clutch | 0 | 0.001 | 0.001 | 0.000 |
| 798 | 0.324 | 0.430 | clutch | 0 | 0.001 | 0.001 | 0.000 |
| 6536 | 0.313 | 0.418 | clutch | 0 | 0.001 | 0.001 | 0.000 |
| 4658 | 0.318 | 0.420 | clutch | 0 | 0.002 | 0.002 | 0.000 |
| 1684 | 0.319 | 0.420 | clutch | 0 | 0.003 | 0.003 | 0.000 |
| 830 | 0.310 | 0.410 | clutch | 0 | 0.003 | 0.003 | 0.000 |
| 6047 | 0.326 | 0.425 | clutch | 0 | 0.006 | 0.006 | 0.001 |
| 3246 | 0.315 | 0.414 | clutch | 0 | 0.007 | 0.007 | 0.001 |
| 9134 | 0.324 | 0.422 | clutch | 0 | 0.007 | 0.007 | 0.001 |
この例でのBonferoni, Holmで有意水準をクリアした頻度を示す.
data.prop.test2 %>%
group_by(is_clutch_hitter)%>%
dplyr::summarise(p_bonferroni = mean(p_bonferroni <= 0.05),
p_holm = mean(p_holm <= 0.05))%>%
mutate_if(is.numeric, funs(round), 3)%>%
gt()| is_clutch_hitter | p_bonferroni | p_holm |
|---|---|---|
| clutch | 0.003 | 0.003 |
| non_clutch | 0.000 | 0.000 |
いずれもclutchタイプの3%程度だけを検出している. これらの方法では, このデータで存在する程度の小さい効果を検出するのは困難のようだ. 一応書いておくと, 野球で仮に個人間で比較して年間成績で出塁率に0.04の差があれば, かなりの得点能力の違いにつながる大きな差であり, 得点圏だけで記録されるとしても多少影響はあるだろう.
予想されるFDRが0.5, 0.4, 0.2, 0.1, 0.05の基準を満たすデータでの, 差がある選手と差がない選手の比率を示す.
data.prop.test2%>%
compute.FDR(Q = 0.5)%>%
mutate_if(is.numeric, funs(round), 3)%>%
gt()| is_clutch_hitter | N | freq |
|---|---|---|
| clutch | 6055 | 0.548 |
| non_clutch | 4992 | 0.452 |
data.prop.test2%>%
compute.FDR(Q = 0.4)%>%
mutate_if(is.numeric, funs(round), 3)%>%
gt()| is_clutch_hitter | N | freq |
|---|---|---|
| clutch | 5188 | 0.64 |
| non_clutch | 2914 | 0.36 |
data.prop.test2%>%
compute.FDR(Q = 0.2)%>%
mutate_if(is.numeric, funs(round), 3)%>%
gt()| is_clutch_hitter | N | freq |
|---|---|---|
| clutch | 3266 | 0.819 |
| non_clutch | 724 | 0.181 |
data.prop.test2%>%
compute.FDR(Q = 0.1)%>%
mutate_if(is.numeric, funs(round), 3)%>%
gt()| is_clutch_hitter | N | freq |
|---|---|---|
| clutch | 1919 | 0.913 |
| non_clutch | 182 | 0.087 |
data.prop.test2%>%
compute.FDR(Q = 0.05)%>%
mutate_if(is.numeric, funs(round), 3)%>%
gt()| is_clutch_hitter | N | freq |
|---|---|---|
| clutch | 1074 | 0.96 |
| non_clutch | 45 | 0.04 |
完全には一致しないが, 概ね期待通りのFDRとなっていることがわかる. FDRの調整を目的とした方法としてはStoreyらによるQ-valueもあり, これはBioconductorパッケージのqvalueで計算可能らしい (使ったことがない).
このようにして得られた候補は一定の間違いを許容した上で, 効果がある選手の候補と言える. しかし, そもそも“効果がある”, あるいは“勝負強い打者である”などと仮に主張したとして, それにどういう意味があるのかはここではそれほど明らかではない. 例えばFAの選手と契約するために評価を行う際に, このような手法が役立つだろうかと問われれば, 無いよりはマシというところだろう. これはこの主張が効果がある/ないという二分法 (dichotomy) をベースとしているためである. 本来ここで我々が知りたいことは, “効果があるかどうか”ではなく, “効果の大きさはどの程度か”という効果の大きさ (size) である.
2.3 回帰で調べる
効果の大きさを表現する方法の一つは平均への回帰の大きさを考慮する方法である. 同じようなデータを2セット生成し, それぞれのセットでの効果 (出塁率の差) の回帰量を調べる. サンプルサイズも大きく (4年分相当 x2), サンプル数も多い (差がない: 90000, 差がある: 10000).
散布図で示す.
RISP.effect <- 0.04
ab.RISP <- 400
ab.not.RISP <- 2000
n1 <- 90000
n2 <- 10000
data <- data.frame(player = 1:n1,
notRISP.x = rbinom(n = n1, size = ab.not.RISP, prob = 0.33) /ab.not.RISP ,
RISP.x = rbinom(n = n1, size = ab.RISP , prob = 0.33) / ab.RISP,
notRISP.y = rbinom(n = n1, size = ab.not.RISP, prob = 0.33) / ab.not.RISP ,
RISP.y = rbinom(n = n1, size = ab.RISP , prob = 0.33) / ab.RISP
)%>%
mutate(diff.x = RISP.x - notRISP.x,
diff.y = RISP.y - notRISP.y)
data2 <- data.frame(player = 1:n2,
notRISP.x = rbinom(n = n2, size = ab.not.RISP, prob = 0.33) /ab.not.RISP,
RISP.x = rbinom(n = n2, size = ab.RISP , prob = 0.33 + RISP.effect) / ab.RISP,
notRISP.y = rbinom(n = n2, size = ab.not.RISP, prob = 0.33) /ab.not.RISP,
RISP.y = rbinom(n = n2, size = ab.RISP , prob = 0.33 + RISP.effect) / ab.RISP)%>%
mutate(diff.x = RISP.x - notRISP.x,
diff.y = RISP.y - notRISP.y)
data2 <- data%>%
mutate(is_clutch_hitter = "non_clutch")%>%
bind_rows(data2 %>%
mutate(is_clutch_hitter = "clutch"))
data2%>%
sample_n(1000)%>%
ggplot(aes(x = diff.x, y = diff.y))+
geom_point()+
theme_bw(base_family = "HiraKakuPro-W3")+
labs(x = "効果x", y = "効果y",
subtitle = "1000人をランダムにサンプリングして示した.")
これは期間xと期間yにおいて見られた効果を比較している. この例における相関係数を求める. p値ではなくr族の効果量で考える, と言ってもいいかもしれない.
cor.test(data2$diff.x, data2$diff.y)##
## Pearson's product-moment correlation
##
## data: data2$diff.x and data2$diff.y
## t = 56.006, df = 99998, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1683783 0.1803972
## sample estimates:
## cor
## 0.1743942この条件では相関はかなり小さい. この相関係数と計算に利用したサンプルサイズから, 相関係数が0.5となるサンプルサイズ (A) を推定する.
R <- as.numeric(cor.test(data2$diff.x, data2$diff.y)$estimate)
N <- 2000 + 400
A <- (1- R) / R * N
A## [1] 11361.92得点圏の比率を固定し, 打席数を11000 * 2で相関係数を再計算してみる.
ab.RISP <- 11000 * 1/5
ab.not.RISP <- 11000 * 4/5
n1 <- 90000
n2 <- 10000
data <- data.frame(player = 1:n1,
notRISP.x = rbinom(n = n1, size = ab.not.RISP, prob = 0.33) /ab.not.RISP ,
RISP.x = rbinom(n = n1, size = ab.RISP , prob = 0.33) / ab.RISP,
notRISP.y = rbinom(n = n1, size = ab.not.RISP, prob = 0.33) / ab.not.RISP ,
RISP.y = rbinom(n = n1, size = ab.RISP , prob = 0.33) / ab.RISP
)%>%
mutate(diff.x = RISP.x - notRISP.x,
diff.y = RISP.y - notRISP.y)
data2 <- data.frame(player = 1:n2,
notRISP.x = rbinom(n = n2, size = ab.not.RISP, prob = 0.33) /ab.not.RISP,
RISP.x = rbinom(n = n2, size = ab.RISP , prob = 0.33 + 0.04) / ab.RISP,
notRISP.y = rbinom(n = n2, size = ab.not.RISP, prob = 0.33) /ab.not.RISP,
RISP.y = rbinom(n = n2, size = ab.RISP , prob = 0.33 + 0.04) / ab.RISP)%>%
mutate(diff.x = RISP.x - notRISP.x,
diff.y = RISP.y - notRISP.y)
data2 <- data%>%
mutate(is_clutch_hitter = "non_clutch")%>%
bind_rows(data2 %>%
mutate(is_clutch_hitter = "clutch"))
data2%>%
sample_n(1000)%>%
ggplot(aes(x = diff.x, y = diff.y))+
geom_point()+
theme_bw(base_family = "HiraKakuPro-W3")+
labs(x = "効果x", y = "効果y",
subtitle = "1000人をランダムにサンプリングして示した.")
cor.test(data2$diff.x, data2$diff.y)##
## Pearson's product-moment correlation
##
## data: data2$diff.x and data2$diff.y
## t = 198.29, df = 99998, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5267931 0.5356905
## sample estimates:
## cor
## 0.5312564概ね0.5となった. これはこの世界では「11000程度のPAがあればその時点でのRISPとnotRISPの差の半分程度を, 実際のスキルであると推定できる.」と記述できそうだ,というイメージである. 仮にこの世界がMLBやNPBと同様に600 PA程度で選手寿命が十数年, また加齢による能力低下があるのであれば, 選手起用や補強において情報として利用することは難しいだろう. あるいは仮に, 1年に10000 PAほどある世界であれば, 能力として十分存在すると言えるかもしれない. しかし, 他の指標の信頼性も改善するため, 結局の所, 利用価値自体は大きくは変わらないだろう.
実際にはこの例は非常に不自然な, 状況別成績における差があるか無いかがはっきり区別された例であるため, これはあまり妥当な結論ではない可能性が高い (散布図を見れば明らかだろう). しかし, 現実では多くの指標ではその効果に個人差があり, 差がない~(検出できないほど) 極めて小さい差~(検出できる程度の) 小さい差~大きな差といった違いがあると考えられ, 回帰の推定は妥当になることが予想される.
この回帰分析の例では, 全体の傾向を考えてRISPによる効果のサイズを推定している. このため多重比較の問題は小さい. 例えば, ここから特定の特徴を共有した選手を見つけて, それらを取り出して個別に回帰分析をすれば, 多重比較的な問題が起こるかもしれない.
3 まとめと議論
出塁率っぽい数値について, 得点圏っぽい状況において差がない場合と差がある場合についてデータを生成し, Bonferroni法, Holm法, BH法を適用し, p値への多重比較補正を試みた.
補正されていないp値は, 実際には差がない場合でも, それなりの数のデータが有意性基準をクリアした. これはp値の定義上, 当然の結果である. 比較する仮説が多い場合, 多重比較の問題のため通常採用される有意水準よりもp値が小さくても, 差がないものを差があると判断するリスクが大きくなる. Bonferroni法, Holm法を適用して, 全体のサンプル数を考慮に入れることで, このような偽陽性 (type I error) を抑制することができた. p値を補正しない問題点と, 補正する利点が現れていると言えるだろう (ここで使った例でFamilywise error, つまり差がないものから一つでも差があるとされるか, がどれくらい起こるかについてはAppendix 2を参照).
しかし, Bonferroni法, Holm法は, ここでの例でデータに含めた差のある選手については, 少数しか検出することはできなかった. 効果の大きさ, サンプルサイズ, α (有意水準), β (1 - 検出力) の間には密接な関係がある. ここでβは「差があるものを差がないという誤った判断 (偽陰性; type II error) をする確率」である. 効果の大きさ, サンプルサイズが一定であれば, αを小さくするとβが大きくなる (つまり検出力が低下する) 性質がある. 個別の比較に関して, Bonferroni法, Holm法はαを小さくすることで偽陽性を抑制するため, 偽陰性の問題はある程度避けがたいコストとなっているだろうと思われる. 特におそらくここで示したような仮説の数 (つまり選手の数) に対して効果が小さい場合では, 問題が非常に大きくなる. 一般的にはα = 0.05, β = 0.2が慣習として用いられることが多い. これはαが扱う誤り, つまり差がないものを差があると判断してしまう誤り (type I error) は, βが対象とする差があるものを差がないとする誤り (type II error) に比べて, 結果として生じる問題が大きいためにこのようなバランスとなっていると考えられる. データ自体の量や質を改善せずに, 差があるものを差がないとする誤りを減らそうとすれば, 帰結として差がないものを差があると判断してしまう誤りが増えることになる. これは一般的に問題が小さいとみなされるような誤りを減らして, 問題が大きいとみなされる誤りを増やすことになるかもしれない.
多重性の補正を考慮したp値は検出力が低い. このため, これらの補正をかけたp値で一般的に採用される有意差が見いだせなかったからと言って, 差がないという結論を出すこともできない (判定は保留される, というぼんやりした結論となることになっている). 実際のところ, 現実にNPBでデータを調べるときに数百人の選手のp値のなかで偽陽性を1つも出さない確率を0.05にする, というのもやや荒唐無稽な目標に思えるかもしれない.
FDRの調整を目的とするBH法も試した. この方法は偽発見を許容することで, p値の補正の程度は小さくなる. このため, FDRを基準として用いることで, βの上昇が抑制されより多くの差がある候補を見つけることができる. しかし, あくまで偽発見を許容した結果であり, 例えば0.05を基準として使えば, 期待値としては得られた候補の5%程度で誤った判断となる. このため, 科学的研究の中ではFDRを基準として得られた候補は, 何らかの独立した実験や調査の対象とされ結論が出されることが多い (ような気がする). このような追加的な調査は野球のデータでは難しいかもしれない. しかし, ある注目している選手が, どの程度の不確実性を許容すれば候補となるか, はある意味で不確実性の記述になっているとも言えるだろう (ただし, これも効果の強さについては何も言えないというべきだと思うが).
野球のデータで多重比較が関連する状況では, 多くの場合, 選手やチーム数が多いことが関与すると思われる. このような状況は回帰の大きさによって表現することが可能であることが多い. この方法は, “効果があるかないか”という閾値を使った二分法ではなく, “効果がどれくらいあるか”という大きさを扱えるため, 仮に効果が小さい状況でもどれくらい小さいのかの表現が可能となる. 野球でこのような効果がどの程度あるかを推測する手段としては, 多くの場合において回帰の大きさを利用する方法が可能である. 典型的には年度間相関が利用されるが, より正確なサンプルサイズを調整する方法によって偶然の影響と能力差の影響を考慮する方法も利用されている. この場合, 注目する選手の結果を全体の傾向で補正することで, 多重性の問題に対処する. この方法では妥当性のある点推定値を提示することが可能となり, 平均的なケースに比べてどの程度の差があるかというサイズの議論が容易になる. 特に比較の数が多くなるために偽陰性率が上昇する状況においては, “効果があるかないか”に注目するp値やそれに類する考え方よりも, 回帰やベイズ統計的な方法で“どれぐらいの効果がありうるのか”の表現を検討することが望ましいだろう.
野球の統計においてこのようなサイズの議論を有効に利用した例としてはTangoらによるThe Bookが挙げられる. The Bookでは平均への回帰を考慮して, 状況別成績 (左右別, あるいは勝負への影響が大きい場面) における個人成績を推定している. また, この推定された効果における個人差の大きさを, 全体的な成績の個人差の大きさと比較することで, 状況別成績の利用価値について議論している.
4 参考文献
@sleep_in_nmbrs (2019). wOBAの誤差範囲: 年度レベルのSD
Rob Arthur and Greg Matthews (2017). Madison Bumgarner Rode A Hot Streak To Greatness, And We Know Who Could Be Next, FiveThirtyEight.
佐藤 文彦（Student) (2018). 捕手の構えどおりに投げることは被安打リスクを下げるか
Christie Ashwanden (2015). Not Even Scientists Can Easily Explain P-values, FiveThirtyEight.
山本陵平. Clinical Journal Club 1. 多重比較.
Cristea IA & Ioannidis JP (2018). P values in display items are ubiquitous and almost invariably significant: A survey of top science journals, PLOS one
Pete Palmer (1990). Clutch Hitting One More Time, BTN.
Tango. DIPS Bands.
蛭川 皓平 (baseballconcrete) (2011). 指標の信頼性と平均への回帰.
Tango (2009). The color of clutch, THT.
Wasserstein RL, Lazar NA., 邦訳 佐藤俊哉 (2016). Editorial: The ASA’s statement on p-values: Context, process, and purpose., The American Statistician, 70, 129-133
Valentin Amrhein, Sander Greenland, and Blake McShane (2019). Retire statistical significance, Nature, 567, 305.
舟尾 暢男. 66. カテゴリデータの検定.
永田 靖, 吉田 道弘 (1997). 「統計的多重比較法の基礎」, サイエンティスト社.
水本篤, 竹内理. 効果量と検定力分析入門 ―統計的検定を正しく使うために―, 『より良い外国語教育研究のための方法』(pp. 47–73).
@sleep_in_nmbrs (2019). wOBAの誤差範囲: スキルの推定 (95% CI, 回帰)
Tango (2010). Regression equations for pitcher events.
Tango, MGL, and Dolphin (2007). The Book, Potmac Books.
5 Appendix
5.1 サンプルを増やしてpが0.05以下となった頻度を計算
冒頭のp値を数を増やして示す.
# 打席数の設定
ab.RISP <- 100
ab.not.RISP <- 500
player.num <- 40000
data <- data.frame(player = 1:player.num,
notRISP = rbinom(n = player.num, size = ab.not.RISP, prob = 0.33),
RISP = rbinom(n = player.num, size = ab.RISP , prob = 0.33))
data.prop.test <- data %>%
mutate(p = apply(data%>%select(-1),
1, # 各行に対して適用
return.p),
notRISP = notRISP /ab.not.RISP,
RISP = RISP /ab.RISP )
data.prop.test %>%
arrange(p)%>%
head(10)%>%
mutate_if(is.numeric, funs(round), 3)## player notRISP RISP p
## 1 17623 0.294 0.52 0
## 2 14327 0.284 0.49 0
## 3 6145 0.304 0.51 0
## 4 19881 0.292 0.49 0
## 5 9038 0.344 0.15 0
## 6 18426 0.338 0.54 0
## 7 11591 0.286 0.48 0
## 8 15137 0.288 0.48 0
## 9 5833 0.362 0.17 0
## 10 31464 0.290 0.48 0pが0.05以下となった頻度を計算する.
mean(data.prop.test$p <= 0.05)## [1] 0.049325だいたい5%となった.
5.2 Familywise error rate
Familywise Error (FWE) rate, つまりこの場合でいうと, 差がない選手の数だけ比較していって, 一つでも差があるという判断をする確率を計算して確認する. 補正方法はBonferroniとHolmを用いる.
return.res <- function(ab.RISP = 100,
ab.not.RISP = 500,
player.num = 300,
success_rate = 0.33){
data <- data.frame(player = 1:player.num,
notRISP = rbinom(n = player.num, size = ab.not.RISP, prob = success_rate),
RISP = rbinom(n = player.num, size = ab.RISP , prob = success_rate))
# 頻度の検定を行いp値を返すための関数
return.p <- function(x){
p <- prop.test(c(x[1], x[2]), c(ab.not.RISP, ab.RISP ), correct=F)[3]
p <- as.numeric(p)
return(p)
}
p.values <- data%>%select(-1)%>%
apply(1, # 各行に対して適用
return.p)
data.prop.test <- data %>%
mutate(p = p.values,
notRISP = notRISP /ab.not.RISP,
RISP = RISP /ab.RISP )
data.prop.test <- data.prop.test %>%
mutate(p_bonferroni = p.adjust(p, "bonferroni"),
p_holm = p.adjust(p, "holm"))
out <- data.prop.test %>%
dplyr::summarise(N= n(),
p_bonf_significant = sum(p_bonferroni <= 0.05),
p_holm_significant = sum(p_holm <= 0.05))%>%
mutate(FWE_bonf = as.numeric(p_bonf_significant != 0),
FWE_holm = as.numeric(p_holm_significant != 0))
return(out)}10000回繰り返してFWEの割合を計算する.
# eval = FALSE
res <- 1:10000%>%
map_df(~return.res())%>%
dplyr::summarise_all(funs(sum))
save(res, file = "Familywise_error_result.rdata")load("Familywise_error_result.rdata")
res%>%
mutate(Trial = 10000,
FWE_rate_bonf = FWE_bonf / Trial,
FWE_rate_holm = FWE_holm / Trial)%>%
select(Trial, everything())## Trial N p_bonf_significant p_holm_significant FWE_bonf FWE_holm
## 1 10000 3000000 492 492 471 471
## FWE_rate_bonf FWE_rate_holm
## 1 0.0471 0.0471概ね期待される0.05だが, 少し低いかもしれない.