統計数理（2017）
第 65 巻 第 2 号 235–249©2017 統計数理研究所

特集「スポーツ統計科学の新たな挑戦」［原著論文］

項目反応理論を用いた野球選手の能力評価指標の提案
阿部 興1 ・作村 建紀2 ・鎌倉 稔成2（受付 2016 年 12 月 26 日；改訂 2017 年 3 月 1 日；採択 4 月 25 日）

要

旨

本研究は野球選手の，特に打撃の能力に関する新しい評価指標の提案を行う．打者の能力を評価する指標として，最も広く使用されるのは打率である．しかし，打率を一つの統計モデルと捉える場合，これは現実的でない仮定のもとに成立するものだと言える．打率はヒットを打つ確率がピッチャーの能力に関わらず常に一定であると見て計算されている．提案モデルはピッチャーによってヒットを打つ難しさが異なり，対戦するピッチャーが選手ごとに異なるという状況の下で，打撃の能力の選手間での比較を可能にする．我々が提案するモデルは，項目反応理論で用いられる 1 パラメータロジスティックモデル（ラッシュモデル）の拡張である．
ラッシュモデルでは，潜在的な能力パラメータは各人一つとされているが，本論文では打者の能力に対応するマルコフパラメータを導入した．すなわち，我々は打席結果にマルコフ性を仮定して調子の波を表現する．これにより調子の波が前の打席の結果を受けて生じると解釈できる．モデルのパラメータの推定にはハミルトニアン・モンテカルロ法を用いる．パラメータ推定の安定性はシミュレーションによって評価する．提案手法の有益性は日本プロ野球の実際のデータを分析することで示す．
キーワード：階層ベイズモデル，MCMC，セイバーメトリクス，ロジスティックモ
デル．
1.

はじめに

本研究は野球選手の，特に打撃の能力に関する新しい評価指標の提案を行う．選手の能力を公平かつ客観的に評価することは，球団を運営する上で重要な課題である．打者の能力を評価する指標として，最も広く使用されるのは打率である．次いで有名な指標として，長打率（OBP），OPS などがある（Albert and Benett, 2003）．長打率は単打に 1，二塁
（SLG），出塁率
打に 2，三塁打に 3，本塁打に 4 の重みをつけた，ヒット数の重み付き平均である．出塁率はヒットと四球と死球の数を足したものを，打数と四球と死球と犠飛の数を足したもので割った指標である．OPS は出塁率と長打率を足し合わせた指標である．Albert and Benett（2003）では，これらの指標をチームごとの平均的な打撃評価指標と得点との相関という観点から検討し，整理している．
より分析的な打者の能力評価指標として，Albert（2008）によるものがある．一般に打者には1

中央大学大学院 理工学研究科：〒 112–8551 東京都文京区春日 1–13–27

2

中央大学 理工学部：〒 112–8551 東京都文京区春日 1–13–27

236

統計数理

第 65 巻 第 2 号 2017

「調子の波」が存在すると言われている．Albert（2008）はベータ・二項モデルを用いて調子の波を表現した．このベータ・二項モデルは 20 打席ごとのヒットの数が二項分布に従うとし，二項分布の成功確率パラメータ p がベータ分布に従うとしたベイズモデルである．ここでのベータ分布は以下のようにパラメタライズされる．
1(1.1)
pKη−1 (1 − p)K(1−η)−1 (K > 0, 0 < η < 1).
B(Kη, K(1 − η))
ここで B(·) はベータ関数である．η は分布の中心を表す．K は精度パラメータで，大きいほど散らばりが小さい．すなわち，K が小さく推定された選手ほど，調子の波が激しいと解釈できる．このモデルは成功確率 p の変化を表現し，これまで評価の難しかった調子の波を定量的に評価することを可能にした．
上述の指標はいずれもピッチャーの能力を考慮していない．ピッチャーによってヒットを打つ難しさが異なり，対戦するピッチャーが選手ごとに異なるにも関わらず，打撃の能力を選手間で比較したい場合，どのようにすればよいか．項目反応理論（item response theory; IRT）の研究成果は，この問題に一つの解を与える．IRT は，教育における達成度評価のような，対象
となる人物
（被験者）
にある課題
（負荷）
が与えられたときに得られる反応から能力測定を行うために考案されたモデルである（Lord, 1952）．IRT の大きな特徴の一つは，被験者に課された課題への反応から，被験者の能力と課題の難しさを同時に評価することである．つまり，被験者の能力を示す能力パラメータと課された課題のレベルを示す項目パラメータを分離して評価するため, 有益な情報を提供し得る．反応を表すモデルは，課題が持つパラメータの数および得られる反応値の種類によってさまざまなものが考案されている（Hambleton et al., 1991;
De Ayala, 2008; Baker, 1992）．リンク関数にはロジスティック関数やプロビット関数が用いられる．またその推定手法についても，周辺最尤推定法（Bock and Aitkin, 1981）やマルコフ連を用いたベイズ推定（Patz and Junker, 1999a, 1999b）など，さまざ
鎖モンテカルロ法（MCMC）
まな手法が確立されており, それを実装したソフトウェアやパッケージなども多い（Bilog-MG,．もともとは教育工学や心理学の分野で発展した理論
2005; Rizopoulos, 2006; Chalmers, 2012）
であるが，近年ではマーケティングへの応用も報告されている（Jong et al., 2008; Raykov and．
Calantone, 2014）
この手法を野球へ応用し，打席ごとのヒット・アウトという応答の確率を，打者の能力を表すパラメータと，投手の難易度（投球の打たれやすさ）
を表すパラメータに分離して推定することが可能である．IRT における能力パラメータが打者の能力パラメータに，項目パラメータが投手の難易度パラメータにそれぞれ対応する．我々が提案するモデルは，IRT で用いられる
ラッシュモデル（1 パラメータロジスティックモデル）の拡張である．ラッシュモデルでは，潜在的な能力パラメータは各人一つとされているが，本論文では打者の能力に対応するマルコフ
パラメータを導入した．ラッシュモデルはその単純さから，解釈が容易であるとされる．我々のモデルもオッズ比を用いることで，選手間の比較を容易なものにできる．
と同様に，成功確率 p が変化することを仮定する．ただし，上述した
我々は Albert（2008）
ベータ二項モデルにおける p の変化は，これまでの打席結果や，打席の状況に依存しないため，解釈がむずかしい．そこで我々は打席結果にマルコフ性を仮定して調子の波を表現する．
これにより調子の波が前の打席の結果を受けて生じると解釈できる．3 節では，MCMC を用いたパラメータの推定方法について述べる．4 節で，パラメータ推定が可能であることを確かめ
る．5.1 節では，2013 年の日本プロ野球の実際のデータに対して分析を行う．

項目反応理論を用いた野球選手の能力評価指標の提案

2.

237

提案モデル

yi,j を打者 j (j = 1, . . . , n) の打席 i (i = 1, . . . , lj ) での結果
（アウトならば 0，ヒットならば1 の値を取る）とする．ここで lj は打者 j の合計の打席数である．xi,j を打者 j が打席 i で対戦した投手とする．xi,j は (1, 2, . . . , m) のいずれかの値をとる．m はリーグ内の投手の数である．表 1 は，選手 j について記録されるデータを，模式的に示したものである．
yi,j はパラメータ pi,j のベルヌーイ分布に従うと仮定し（pi,j = Pr(yi,j = 1)），pi,j に以下のような回帰型の構造を仮定する．
(2.1)

logit(pi,j ) = β0,j + β1,j yi−1,j + τ (xi,j ).

ここで，τ (x) は，(2.2)

τ (x) =

⎧
⎪
b1
⎪
⎪
⎪
⎪
⎨b2
..
⎪
⎪
⎪.

⎪
⎪
⎩b

m

x=1
x=2
..
.
x=m

なる関数である．投手の難易度パラメータ b（
は平均が 0 になるよう基準化されk k = 1, . . . , m）
2ているとする．事前分布として，bk は平均 0，分散 σ の正規分布に従うとする．β0,j は平均
μ，分散 ξ 2 の正規分布を仮定する．β1,j には，無情報階層事前分布として，区間 (−∞, ∞) の一様分布を仮定する．階層事前分布として, σ, ξ には区間 [0, ∞), μ には区間 (−∞, ∞) の一様分布を仮定する．
β1,j はヒットを打つ確率に対する，直前の打席の影響を表すパラメータである．我々は β1,jの 90% 信用区間を求め，これが 0 を含まない選手を「調子の波がある選手」，0 を含む選手を「調子の波がない選手」とする．
2.1 オッズ比ロジスティックモデルはオッズ比による解釈が容易である．本研究では 3 種類のパラメータに対してオッズ比を評価する．一つ目は打者同士のベースラインとなる能力を比較するためのオッズ比である．直前の打席を失敗とし，ピッチャーを特定の選手に固定した場合にヒットを打つ確率に着目すると，打者 j の打者 h に対するオッズ比は，(2.3)

exp(β0,j )/ exp(β0,h )である．二つ目は投手の能力を比較するためのオッズ比である．投手難易度パラメータは平均が 0 になるよう基準化されていることを仮定したため，ここでは 0 を基準としたオッズ比,
表 1．解析対象となるデータの一例．

統計数理

238

第 65 巻 第 2 号 2017

(2.4)

exp(bj )を考える．このオッズ比は平均的なピッチャーと比べて，どの程度ヒットという事象が起こりやすいかという指標であり，値が小さいほど能力の高いピッチャーと解釈できる．三つ目は打者の調子の波を把握するためのオッズ比である．
(2.5)

exp(β1,j )これは直前の打席結果がアウトだったときを基準として，直前の打席結果がヒットだったときはどの程度ヒットを打ちやすいかという指標である．5.1 節の事例研究ではこれらを使用して選手を評価する．
2.2 定常分布
本モデルでは，投手の難易度パラメータの平均が 0 になるよう基準化されていることを仮定した．打者 j が平均的な投手と対戦した場合のヒットの確率は，
P (yi,j = 1) = logit−1 (β0,j + β1,j yi−1,j )．

である．一打席前でヒットを打てなかった場合にヒットを打つ確率を，
1
p0,j =(2.6)
1 + exp(−β0,j )一打席前でヒットを打った場合にヒットを打つ確率を，
1
p1,j =(2.7)
1 + exp(−(β0,j + β1,j ))と表すことにする．これは遷移行列
Pj =

(2.8)



1 − p0,j
1 − p1,j

p0,j
p1,j



の 2 状態マルコフ連鎖を形作る．このマルコフ連鎖は既約かつ非周期的（Karlin, 1969）であり，以下を満たす行ベクトル π j が定義できる．
πj P j = πj

(2.9)πj = (π0j , π1j ) について解くと,(2.10)π0j =1 − p1j,
1 + p0j − p1j

π1j =

p0j
1 + p0j − p1jである．πj は定常分布と呼ばれる．π1j は十分な時間が経過したときの，選手 j がヒットを打つ割合と解釈でき，初期の状態には影響されない．これも指標の一つとして 5.1 節の事例研究で使用する．
3.

パラメータの推定方法

式（2.1）で示したモデルの事後分布の統計量を求めるために，本研究ではハミルトニアン・モ．
ンテカルロ法（ハイブリッドモンテカルロ法ともいう）による MCMC を用いた（Bishop, 2006）
MCMC は推定の対象となる事後分布に従う乱数を，マルコフ連鎖を用いてシミュレートする方法である．近年，IRT の分野では MCMC が積極的に利用されており，MCMC による推定で議論されているが，ハミルトニアン・モンテカルロ法を用いる方法にの性質は Fox（2010）

項目反応理論を用いた野球選手の能力評価指標の提案

239

ついては記載がない．MCMC 法のバリエーションとしてギブスサンプリング法やメトロポリ
ス・ヘイスティングス法はよく知られており，IRT の分野でも盛んに利用されている．ギブス
サンプリング法ではパラメータごとに全条件付き分布を導出し，変数ごとにサンプリングを行う．メトロポリス・ヘイスティングス法ではパラメータごとに適切な提案分布を選ぶ必要がある．本解析のようにパラメータが多い場合，ハミルトニアン・モンテカルロ法が便利である．
MCMC による IRT モデルのパラメータ推定法の中でも特に有名な手法として Albert（1992）により提案された，データ拡大（data augmetation）アルゴリズムがある．0 または 1の値を取る打席結果データ yi,j に対して，連続値を取る潜在変数 z を仮定する．z は平均
ηi,j = β0,j + β1,j yi−1,j + τ (xi,j )，分散 1 の正規分布に従うとし，zi,j > 0 ならば yi,j = 1，zi,j ≤ 0
ならば yi,j = 0 となる．言い換えると，yi,j = 1 という観測は 0 で右打切りされた zi,j の観測（Turnbull, 1976）された zi,j の観測値である．こ値であり yi,j = 0 という観測は 0 で左打切りの場合尤度関数は
lj
n 


(3.1)

Φ(0|ηi,j , 1)yi,j (1 − Φ(0|ηi,j , 1))1−yi,jj=1 i=2と 書 け る ．本 研 究 で は リ ン ク 関 数 に ロ ジ ッ ト 関 数 を 用 い る ．こ の 場 合 も Albert（1992）と 同 様 の 議 論 に よ り ，潜 在 変 数 z を 仮 定 す る と z は ロ ケ ー シ ョ ン パ ラ メ ー タ が ηij =β0,j + β1,j yi−1,j + τ (xi,j ) のスケールパラメータが 1 のロジスティック分布に従うと考える．F (z|η) をスケールパラメータ 1，ロケーションパラメータ η のロジスティック分布の分布関数とすると，尤度関数は,
lj
N 


(3.2)

F (0|ηi,j )yi,j (1 − F (0|ηi,j ))1−yi,jj=1 i=2と書ける．これはベルヌーイ分布にロジットリンク関数を仮定した場合の尤度と全く同一である．事後分布は,(3.3) π(β0 , β1 , b|y) = C

m

i=1

φ(bi , 0, σ)

m

i=1

φ(β0,i , μ, ξ)

n lj



F (0|ηij )yij {1 − F (0|ηij )}1−yij

j=1 i=2となる．ここで φ は正規分布の密度関数， C は規格化定数である．Albert（1992）は各パラ
メータについての全条件付き分布を導出しているが，ハミルトニアン・モンテカルロ法ではそが分かればサンプリングできる．
の必要がなく，尤度関数と事前分布の積
（3.3）
4.

シミュレーション

IRT の分野では，項目パラメータと能力パラメータが同時に推定される際，最尤推定法や条件付き最尤推定法を用いた場合，項目パラメータは一致推定量にならないことが知られている（豊田，2005）．そのため，パラメータが安定して推定されるか検討が必要である．ここでは，シミュレーションを通じて，パラメータ推定の性質を調べる．シミュレーションでは打者を 75 人，投手を 165 人に設定し，打数は 250 で固定した．これらの設定は事例研究で用いるデータと近い数字にしている．打者能力パラメータ β0,j , β1,j (j = 1, . . . , 75) と，投手難
易度パラメータ bk (k = 1, . . . , 165) を標準正規乱数で生成した．ピッチャー xi,j は各打者の各打席ごとに一様乱数でランダムに選んだ．選択されたピッチャーに対応するパラメータ bk
から確率 logit−1 (β0,j + bk ) で 1，確率 1 − logit−1 (β0,j + bk ) で 0 を生成し，続けて成功確率

240

統計数理

第 65 巻 第 2 号 2017

図 1．推定値と真値の比較．横軸がシミュレーションで仮定した真値，縦軸が推定値である．
（b）打者能力パラメータ β0 ，（c）打者能力パラメータ β1 ．
（a）投手難易度パラメータ b，logit−1 (β0,j + β1,j yi−1,j + bk ) のベルヌーイ試行を繰り返して yi,j を生成した．こうして生成した yi,j と xi,j から元のパラメータを推定する．MCMC の試行回数は 2000 回とし，うち前半の 1000 をバーン・イン期間として捨てた．連鎖の構成数は 3 とした．事後期待値をパラメー（a）–1（c）に示す．こタの点推定値とした．推定値と元のパラメータを比較したプロットを図 1
れらの図は，横軸にシミュレーションで生成した真値，縦軸に推定値を取ってプロットした．
両者でパラメータが一致していれば，散布図の各点が直線状に並ぶが，点はおおむね直線状に分布している．選手間の能力の相対的比較という観点からは，パラメータの大小の順序の関係が保たれていることが必要であり，そのことは相関係数で評価できる．推定値と真値の相関係数は，b については 0.98, β0 が 0.97, β1 が 0.95 と高い値になった．決定係数はそれぞれ 0.95,
0.94, 0.89 であった．
また，区間推定には 90% 信用区間を使用した．事後分布の両端から確率 α/2 の部分を切り取って残った，中央部の確率 (1 − α) に対応する区間を 100 × (1 − α) 信用区間と呼ぶ．ベイズモデルではパラメータが確率変数であると考えるので，90% 信用区間は真値を含む確率が90% となる区間である．本解析では，MCMC サンプルの標本 5% 点と 95% 点をとり，それぞれを信用下限，信用上限とすることで，数値的な 90% 信用区間を得た．表 2 より，信用区間がその中にシミュレーションで仮定した真値を含む割合を計算すると，90% に近い値になっている．このことが，推定の信頼性を高めている．
投手難易度パラメータの事前分布のパラメータ σ については，事後期待値が 1.05，打者能項目反応理論を用いた野球選手の能力評価指標の提案

241

表 2．推定された 90% 信用区間が真値を含んだ回数と割合．

力パラメータの事前分布のパラメータについては μ が −0.11, ξ が 0.98 であった．それぞれシ
ミュレーションで仮定した真値 σ = 1, μ = 0, ξ = 1 に近い値が推定された．
5.

事例研究

ここで扱うデータは，2013 年のセ・リーグ，パ・リーグの試合の打席結果である．ただし，
四球，死球，敬遠，犠打，犠飛は打席から除いた．年間を通じて 100 回以上打席に立った打者と，その選手に対して投球を行った投手を解析対象とした．セ・リーグとパ・リーグは分けて分析を行った．対象となる打者はセ・リーグが 75 名，パ・リーグが 74 名，対象となる投手
はセ・リーグが 165 名，パ・リーグが 162 名である．条件を満たす打者の打席数の平均はセ・リーグで 264.75，パ・リーグで 284.36 であった．
5.1 推定結果の概要
これまでに述べたように，推定にはハミルトニアン・モンテカルロ法を用いた．MCMC の
試行回数は 10000 回とし，うち前半の 5000 をバーン・イン期間として捨てた．連鎖の構成数は 3 とした．
収束の判定はトレースラインの目視による確認と Gelman（1996）の収束判定指標 R̂ によりでは，R̂ が 1.2 または 1.1 より小さければ収束したと判断してよいと行った．Gelman（1996）
されている．本解析においてはすべてのパラメータについて，R̂ が 1.1 より小さいことが確認された．
投手難易度パラメータ bk は値が小さいほど，打つのが難しい投手である．事後期待値の値が小さい順に上位 10 名の選手の推定結果を表 3 に示す．90% 信用区間とオッズ比，被安打率
もあわせて表示した．被安打率はピッチャーの能力を表す指標のひとつで，被安打数を打数で割ったものである．また bk の事後期待値を全投手についてヒストグラムで図示した（図 2（a）–2（b））．図の点線は bk の事後期待値の平均である．これらの図から bi は 0 付近を中心に分布していることがわかる．
表 3 を見ると，被安打率が 0.2 の金子選手が 1 位にランキングされ，被安打率が 0.16 のウィリアムス選手が 10 位にランキングされるという逆転が見られる．ただし，ウィリアムス投手
は打席数が少ないため，信用区間の幅は広く推定されている．この結果がなぜ生まれたかを検討する．金子選手とウィリアムス投手に関して，打率三割以上の打者と対戦したときの被安打率と，打率三割未満の打者と対戦したときの被安打率を集計した（表 4）．ウィリアムス投手は打率の低い選手に対しての成績が好調な反面，打率の高い選手には多くヒットを打たれている．対照的に金子投手はウィリアムス投手と比べ，打率の低い選手に対しての被安打率は低くないが，打率の高い選手のヒットを少なく抑えている．そのため，被安打率と比較したときに特異な値を取る．本モデルは打者の能力と投手の能力を同時に評価しているため，このような場合，能力の高い打者を多くアウトにしている投手が高く評価される．このように提案モデル
は，被安打率だけでは見落としてしまうおそれのある優秀な（能力の高い打者に強い）投手をピックアップできる．

統計数理

242

（a）

第 65 巻 第 2 号 2017（b）

図 2．投手難易度パラメータ bi の事後期待値の分布（（a）セ・リーグ，（b）パ・リーグ）．
表 3．投手難易度パラメータの推定結果．

表 4．打率で層別した被安打率．

打者能力パラメータ β0,j は，直前の打席結果がアウトだったときのヒットを打つ能力を表す．野球ではヒットよりもアウトのほうが起こりやすいため，β0,j は選手 j のベースラインとなる能力と解釈できる．解釈しやすくなるよう，逆ロジット関数で変換した β0,j の事後期待値に加えてその 90% 信用区間，オッズ比，2.2 節で検討した π1j の β0,j , β1,j からのプラグイン
推定値を算出した．オッズ比は β0,j の値を降順に並べて 10 番目の選手を基準としたオッズ比である．表 5 がセ・リーグの選手，表 6 がパ・リーグの選手についての結果である．これらは項目反応理論を用いた野球選手の能力評価指標の提案

243

表 5．打者能力パラメータの推定結果（セ・リーグ）．

表 6．打者能力パラメータの推定結果（パ・リーグ）．

いわば対戦した投手の効果を補正した打率である．
また打者能力パラメータ β1,j は値が正であれば，直前の打席でヒットを打つと続けてヒットを打つ確率が上がり，負であれば直前の打席でヒットを打てなかった場合にヒットを打つ確率が上がると解釈できる．MCMC サンプルから β1,j の 90% 信用区間を求め，これが 0 を含まないとき，調子の波が存在する選手とする．そのような選手を表 7–8 に示した．
表 5 の logit−1 (β0 ) を見ると，阿部選手は打率と比較して β0 が大きく推定されており，梶谷
選手に関しては打率と比較して β0 が小さく推定されている．それに対応して，阿部選手は β1が小さい値を取り，梶谷選手は β1 が大きい値を取っている．β0 と β1 を総合した指標である
π1 は，打率に近い値を取っている．表 7 より，阿部選手，梶谷選手は「調子に波がある」
プレー
ヤーである．
5.2 カイ二乗検定を用いた考察
ここでは，古典的なカイ二乗検定を用いて調子の波について分析を行い，提案モデルの方法
が既存の結果と矛盾しないことを示す．調子の波を定量的に評価するもう一つの方法は，ある打席での成功・失敗と，その直前の打席での成功・失敗との関係を調べることである．このアが用いたものである．いま 2 × 2 分割表の k 行 l 列目のセルを ckl とプローチは Albert（2008）
する（k = 1, 2, l = 1, 2）．選手 j についての分割表を，c11 = i I{yi−1,j = 0, yi,j = 0}，c12 =I{yi−1,j = 0, yi,j = 1}，c21 = i I{yi−1,j = 1, yi,j = 0}，c22 = i I{yi−1,j = 1, yi,j = 1} と
i

244

統計数理

第 65 巻 第 2 号 2017

表 7．打者能力パラメータ β1 の推定結果（セ・リーグ）．

表 8．打者能力パラメータ β1 の推定結果（パ・リーグ）．

表 9．阿部慎之助選手についての分割表．

して作成する．I は指標関数であり，等式を満たすとき 1，他の場合 0 の値を取る．具体例として阿部慎之助選手についての分割表を，表 9 に示した．
この表に対して独立性についてのカイ二乗検定を行った．検定の結果，帰無仮説が棄却された選手は，直前の打席での成功・失敗が，次の打席結果に影響を与えると解釈できる．10% 水準で有意になった選手を「調子の波がある」選手とする．
「調子の波がある」選手を表 10 にまとめた．理解を容易にするために，クラメールの連関係数も併記した．カイ二乗検定では有意と「調子の波がある」選手については，並べて記載しならない選手についても，表 7–8 に上げた
た．有意となった選手には，p 値の右に傍点を振った．
連関係数が正の値をとるときは，直前の打席結果がヒットだった場合，次の打席でのヒットが出やすくなると解釈でき，負の値をとるとき，直前の打席結果がアウトだった場合，次の打「調
席でのヒットが出やすくなると解釈できる．表 7–8 と表 10 を比較すると，カイ二乗検定で子の波がある」と判断された選手は，提案手法で「調子の波がある」と判断された選手にすべて含まれている．このことは提案手法が従来手法と矛盾しないことを示している．逆に提案手法
で「調子の波がある」と判断された選手が，カイ二乗検定で「調子の波がある」と判断されるとは限らない．これは提案手法がヒット・アウトの系列の情報や相手投手の情報を利用して，より感度の高い検出を行っているためである．また仮説検定では，打席結果の情報を，
「有意」
か
「有
意でない」
かという二値の解釈に落としている．提案手法は各パラメータから調子の波の強さ

項目反応理論を用いた野球選手の能力評価指標の提案

245

表 10．カイ二乗検定の結果 10% 水準で有意になった打者一覧（表 7–8 に上げた選手もあわせて記載）．

や，その信用区間といった，より多くの情報を与えることができる．
5.3 ヒットの予測提案モデルでは打者 j とその直前の打席結果，および投手 k が決まれば，ヒットを打つ確率を予測できる．直前の打席結果がアウトだった場合，打者 j と投手 k の対戦でヒットが出る確
率は，
(5.1)

p0,j,k = logit−1 (β0,j + bk ).

直前の打席結果がヒットだった場合，(5.2)

p1,j,k = logit−1 (β0,j + β1,j + bk )と表す．直前の打席結果がアウトだった場合の，打者 j と投手 k の対戦回数を n0,j,k ，そのヒットの回数を h0,j,k と表す．直前の打席結果がヒットだった場合の，打者 j と投手 k の対戦回数
を n1,j,k ，そのヒットの回数を h1,j,k と表す．実際のヒットの割合 h0,j,k /n0,j,k , h1,j,k /n1,j,k と，提案モデルによる確率の推定値 p0,j,k , p1,j,k を比較する．以降，単に実際の割合といった場合
は h0,j,k /n0,j,k または h1,j,k /n1,j,k ，推定確率といった場合は p0,j,k または p1,j,k を指すものと（a）
および図 3（b）に示す．n0,j,k または n1,j,k が大きする．推定確率と実際の割合の差を，図 3
くなるにつれ，実際の割合と予測確率の差は小さくなる傾向が見られた．
n0,j,k または n1,j,k が一定以上大きいにも関わらず，予測確率と実際の割合が大きく食い違
う場合，その打者と投手の組み合わせに固有の特殊な事情がある可能性を考える．n0,j,k または n1,j,k が 10 を超える組み合わせについて，実際の割合と予測確率の差が 0.2 を超えたものを表 11，12 にまとめた．
（表 5）優れた
阿部慎之助選手はベースラインとなる能力パラメータ β0 が大きく推定された選手だが，スタンリッジ投手とは 13 回対戦して 1 回しかヒットを打てていない．いわばスタ
ンリッジ投手は阿部慎之助選手の天敵といえるだろう．本モデルは打者の能力と投手の能力を測る基準となるため，このような特異な関係性を見出すことに利用できる．

統計数理

246

（a）

第 65 巻 第 2 号 2017（b）

図 3．推定確率と実際の割合の差．横軸に n，縦軸に差を取った（（a）セ・リーグ，（b）パ・リーグ）．
表 11．予測確率と実際の割合の比較（セ・リーグ）．

表 12．予測確率と実際の割合の比較（パ・リーグ）．

5.4 パラメータ推定の安定性
4 節ではシミュレーションデータからパラメータ推定の安定性を議論したが，ここでは実データを用いて改めて検討する．一般に項目反応理論では項目の困難度パラメータを正確に推定するためには，多くの被験者が必要であるとされている．被験者の数を少なくとった場合でも，被験者の数を多くした場合と同様の推定結果が得られるならば，モデルとそのパラメータはより信頼できるものとなる．そこでパラメータ推定の対象となる打者が約半数となるように

項目反応理論を用いた野球選手の能力評価指標の提案

（a）

247

（b）

図 4．投手難易度パラメータの推定結果の比較（（a）セ・リーグ，（b）パ・リーグ）．

ランダムサンプリングで打者を選び直した新たなデータセットに対してモデルのフィッティン
グを行い，5.1 節で推定した投手難易度パラメータと，新たなデータセットから推定した投手（a），パ・リーグについての比
難易度パラメータを比較した．セ・リーグについての比較が図 4（b）である．これらの図は縦軸にランダムサンプリングしたデータセットから推定し較が図 4た投手難易度パラメータ，横軸に 5.1 節で推定した投手難易度パラメータをとりプロットした．
点はおおむね直線状に分布しており，相関係数はセ・リーグで 0.81，パ・リーグで 0.77 と高い値になった．決定係数はそれぞれ 0.77, 0.48 であった．このことはサンプルサイズが小さくても結果に矛盾が生じないことを意味しており，安定したパラメータの推定が行えた見込みが高（b）においては，線形回帰した場合に傾きが 1 より小さくなるよういといえる．ただし，図 4な系統的な差が生まれている．ベイズモデルではサンプルサイズが小さい場合，事前分布の影
響で推定値がシュリンクするためである．順序の関係は保たれているため，選手間の相対的比
較をする上では大きな問題とならないが，提案モデルを利用する際は，この性質を認識する必
要がある．
6.

まとめ

本研究では打席結果を 0 または 1 にコード化したデータから，投手の能力と打者の能力を推定するモデルを提案した．ハミルトニアン・モンテカルロ法を用いたパラメータ推定の安定性はシミュレーションと事例研究で確かめた．事例研究では，従来用いられる被安打率では低く評価されるおそれのある優れた選手を見つけることができた．また 0 または 1 のプロセスにマルコフ性を仮定することで，打者の「調子の波」を捉えることができた．
謝

辞

本研究の貸与データは情報・システム研究機構の新領域融合研究プロジェクト『社会コミュ
ニケーション』
データ中心科学リサーチコモンズ事業『人間・社会データ』の支援を受けたもの
No.15K21379）の助成を受けたものです．
です．また，本研究の一部は科学科研費
（若手研究
（B）
心より感謝申し上げます．

248

統計数理

参

第 65 巻 第 2 号 2017

考

文

献

Albert, J. (1992). Bayesian estimation of normal ogive item response curve using Gibbs sampling,
Journal of Educational Statistics, 17(3), 251–269.
Albert, J. (2008). Streaky hitting in baseball, Journal of Quantitative Analysis in Sports, 4(1), DOI:
10.2202/1559–0410.1085.
Albert, J. and Bennett, J. (2003). Curve Ball, Springer-Verlag, New York. （加藤貴昭, 後藤寿彦 訳(2004).『メジャーリーグの数理科学〈上〉』，シュプリンガーフェアラーク，東京．
）
Baker, F. B. (1992). Item Response Theory: Parameter Estimation Techniques, Marcel Dekker, NewYork.
Bilog-MG (2005). Scientific Software International, http://www.ssicentral.com/irt/index.html.
Bishop, C. (2006). Pattern Recognition and Machine Learning (Information Science and Statistics),
Springer, Cambridge.（元田浩, 栗田多喜夫, 樋口知之, 松本裕治, 村田昇 訳 (2008).『パターン認識と機械学習 下（ベイズ理論による統計的予測）』，丸善出版．東京．）
Bock, R. D. and Aitkin, M. (1981). Marginal maximum likelihood estimation of item parameters:
Application of an EM algorithm, Psychometrika, 46(4), 443–459.
Chalmers. R. (2012). Mirt: A multidimensional item response theory package for the r environment,
Journal of Statistical Software, 48(1), 1–29.
De Ayala, R. J. (2008). The Theory and Practice of Item Response Theory, The Guilford Press, NewYork.
Fox, J. P. (2010). Bayesian Item Response Modeling Theory and Applications, Springer-Verlag, NewYork.
Gelman, A. (1996). Markov Chain Monte Carlo in Practice, Chapman & Hall/CRC Interdisciplinary
Statistics, London.
Hambleton R. K., Swaminathan, H. and Rogers, H. J. (1991). Fundamentals of Item Response Theory,
Vol. 2, Sage Publications, New York.
Jong, M. G. D., Steenkamp, J. B. E. M., Fox, J. P. and Baumgartner, Hans (2008). Using item response theory to measure extreme response style in marketing research: A global investigation,
Journal of Marketing Research, 45(1), 104–115.
Karlin, S. (1969). A First Course in Stochastic Processes, Academic Press, Cambridge.（佐藤健一, 佐）藤由身子 訳 (1974).『数理解析とその周辺 3 確率過程講義』，産業図書，東京．
Lord, F. M. (1952). A theory of test scores, Psychometric Monographs, No. 7, Psychometric Corporation, Richmond.
Patz, R. J. and Junker, B. W. (1999a). A straightforward approach to Markov chain Monte Carlo
methods for item response models, Journal of Educational and Behavioral Statistics, 24(2),146–178.
Patz, R. J. and Junker, B. W. (1999b). Applications and extensions of MCMC in IRT: Multiple item
types, missing data, and rated responses, Journal of Educational and Behavioral Statistics,
24(4), 342–366.
Raykov, T. and Calantone, R. J. (2014). The utility of item response modeling in marketing research,
Journal of the Academy of Marketing Science, 42(4), 337–360.
Rizopoulos, D. (2006). ltm: An R package for latent variable modeling and item response analysis,
Journal of Statistical Software, 17(1), 1–25.
豊田秀樹 (2005).『項目反応理論［理論編］』，朝倉書店，東京．
Turnbull, B. W. (1976). The empirical distribution function with arbitrarily grouped, censored andtruncated data, Journal of the Royal Statistical Society, Series B (Methodological), 38, 290–
295.

Proceedings of the Institute of Statistical Mathematics Vol. 65, No. 2, 235–249 (2017)

249

Measurements of Baseball Players’ Batting Abilities
Ko Abe1 , Takenori Sakumura2 and Toshinari Kamakura2
1 Graduate School of Science and Engineering, Chuo University

2 Department of Industrial and Systems Engineering, Chuo University

Statistics of player performance is an important part of baseball. Many stats havebeen proposed to measure a batter’s performance, including batting average, on-base percentage, and slugging percentage. In the ﬁeld of baseball analytics, the “streakiness” ofbatter’s ability is often discussed using a binary sequence of hitting outcomes for a player
during a season.
Unlike previous studies, which use data from the batter, we take a diﬀerent approach.
To analyze a batter’s performance, we simultaneously model the pitcher and batter’s ability. To model a batter’s streakiness, we employ an extension of a one-parameter logistic
item response model. Item response theory (IRT) estimates both the subject’s ability
and item diﬃculty. In this study, the ability parameter and item diﬃculty parametercorrespond to the batter’s ability and pitcher’s ability, respectively. Although simplicity is thought to make the one-parameter logistic model easy to interpret, our model
incorporates numerous parameters. However, using the odds ratio allows athletes to becompared.
We express streakiness by the interactions of previous at bats and imposing theMarkov property on batting data. Speciﬁcally, we use MCMC in the Hamiltonian Monte
Carlo method (also called the hybrid Monte Carlo method). The computation of Gibbs
sampling is complex and time consuming, but the Hamiltonian Monte Carlo method iseasily computed once the prior distribution and the likelihood function are deﬁned. Our
simulation study shows that the true and estimated values agree well. Additionally, thecalculated proportion of times that the credible interval contains the true value is close
to the nominal value.
To demonstrate the usefulness of our proposed method, we applied it to analyze actual
data from Japanese professional baseball. Two-way tables can measure the dependence ofthe previous success and the current success by the Pearson chi-square statistic and thecorresponding p-value of the test of independence. The results provide more information
and are consistent with the results of chi-square test. Because comparing streakiness inthe hypothesis test is diﬃcult, we ranked streaky players from the credible intervals andthe posterior means. IRT requires many subjects to estimate item diﬃculty parameters.
Although we estimated the parameters using fewer batters, the results from our method
are similar to those from IRT.

Key words: Bayesian hierarchical model, MCMC, sabermetrics, logistic model.