Bradley-Terryモデルとその応用まとめ - RigelのR言語メモであーる(主にpython)

スポーツの勝敗データを想定する。いろんなチーム(あるいは個人)のうち2チームが戦って勝敗を記録したデータ。
このデータから各チームの潜在的な「強さ」を推定してどのくらいの確率でどっちが勝つのかを推定するというお話。

Bradley-Terryモデルとは

Bradley-Terryモデルは提案されたのは1952年となかなか古い統計モデル。早速モデル化していく。
$n$ チームが一対一で何回か対戦する。引き分けは考えない。
チーム $i$ がチーム $j$ に勝つ確率を $p_{ij}$ を
${ \displaystyle p_{ij} = \frac{\pi_i}{\pi_i+\pi_j}\ (i\not=j;i,j=1,2,...,n) }$
と表す。 $\pi_i,\pi_j$ はそれぞれチーム $i$ の「強さ」、チーム $j$ の「強さ」とする。

以上のシンプルなモデルである。

もうちょっと式変形。
引き分けはないものとしているので、チーム $j$ がチーム $i$ に勝つ確率を $p_{ji}$ は
${ \displaystyle p_{ji} = 1-p_{ij}=\frac{\pi_j}{\pi_i+\pi_j} }$
となる。
式変形の見通しをよくするために、 $\pi_{i}=\exp(V_i)$ と表すと
${ \displaystyle \log(\frac{p_{ij}}{1-p_{ij}})=\log(\frac{\frac{\pi_i}{\pi_i+\pi_j}}{\frac{\pi_j}{\pi_i+\pi_j}}) =\log(\frac{\pi_i}{\pi_j})=V_i-V_j }$
と、一般化線形モデル(ロジスティック回帰)の形に帰着する。
こうしてみると、こんなシンプルなモデル(ただの確率の推定)なら、1952年とかに提案されててもおかしくないなぁと思う。
パラメータ $\pi_i\ (i=1,2,...,n)$ の推定はラグランジュの未定乗数法から出発して、連立方程式を数値計算で近似的に求める(詳細はパス)。

BradleyTerry2パッケージで実際にやってみる。

BradleyTerry2はCRANにある。
パッケージにあるサンプルデータを見てみる。

> data(baseball)
> head(baseball)
  home.team away.team home.wins away.wins
1 Milwaukee   Detroit         4         3
2 Milwaukee   Toronto         4         2
3 Milwaukee  New York         4         3
4 Milwaukee    Boston         6         1
5 Milwaukee Cleveland         4         2
6 Milwaukee Baltimore         6         0

左から、ホームチーム名、アウェイチーム名、ホームチーム勝数、アウェイチーム負数である。
パラメータ推定してみる。

> baseballModel1 <- BTm(outcome = cbind(home.wins, away.wins), 
+                       player1 = home.team, player2 = away.team, 
+                       data = baseball, id = "team")
> baseballModel1
Bradley Terry model fit by glm.fit 

Call:  BTm(outcome = cbind(home.wins, away.wins), player1 = home.team, 
    player2 = away.team, id = "team", data = baseball)

Coefficients:
   teamBoston  teamCleveland    teamDetroit  teamMilwaukee   teamNew York    teamToronto  
       1.1077         0.6839         1.4364         1.5814         1.2476         1.2945  

Degrees of Freedom: 42 Total (i.e. Null);  36 Residual
Null Deviance:	    78.02 
Residual Deviance: 44.05 	AIC: 140.5
> BTabilities(baseballModel1)
            ability      s.e.
Baltimore 0.0000000 0.0000000
Boston    1.1076977 0.3338779
Cleveland 0.6838528 0.3318764
Detroit   1.4364084 0.3395682
Milwaukee 1.5813559 0.3432557
New York  1.2476178 0.3358606
Toronto   1.2944851 0.3366691

一番強いのはMilwaukeeと推定された。
推定されたパラメータは $V_i$ なので、
例えば一番弱いと推定されたBaltimore( $i$ )がMilwaukee( $j$ )に勝つ確率は、
${\begin{eqnarray} p_{ij}&=&\frac{\pi_i}{\pi_i+\pi_j}\\ &=&\frac{\exp(V_i)}{\exp(V_i)+\exp(V_j)}\\ &=&\frac{\exp(0.0000000)}{\exp(0.0000000)+\exp(1.5813559)}\\ &=&0.1706035 \end{eqnarray}}$
となる。
データからBaltimore対Milwaukeeは13戦していてBaltimoreの2勝なので、いいぐらいかなと思う。

次にホームチームにバイアスがかかることを考慮したモデルをやってみる。
ホームチームの勝率が上がるというモデル。
${ \displaystyle \log(\frac{p_{ij}}{1-p_{ij}})=V_i-V_j+\alpha z }$
というモデル。ここで、 $z$ はチーム $i$ がホームチームの時に $1$ 、チーム $i$ がアウェイチームの時に $-1$ となる変数で $\alpha$ がその係数。
ホームチーム $i$ がアウェイチーム $j$ に勝つ確率は、
${ \displaystyle p_{ij} = \frac{\exp(V_i+\alpha)}{\exp(V_i+\alpha)+\exp(V_j)} }$
アウェイチーム $i$ がホームチーム $j$ に勝つ確率は、
${ \displaystyle p_{ij} = \frac{\exp(V_i)}{\exp(V_i)+\exp(V_j+\alpha)} }$
ホームチーム $i$ がアウェイチーム $j$ に負ける確率は、
${ \displaystyle p_{ji} = \frac{\exp(V_j)}{\exp(V_i+\alpha)+\exp(V_j)} }$
アウェイチーム $i$ がホームチーム $j$ に負ける確率は、
${ \displaystyle p_{ji} = \frac{\exp(V_j+\alpha)}{\exp(V_i)+\exp(V_j+\alpha)} }$
です。（ちょっと自信ない）
パラメータ推定してみる。

> baseball$home.team <- data.frame(team = baseball$home.team, at.home = 1)
> baseball$away.team <- data.frame(team = baseball$away.team, at.home = 0)
> baseballModel2 <- BTm(outcome = cbind(home.wins, away.wins), 
+                       player1 = home.team, player2 = away.team, 
+                       data = baseball, id = "team", formula = ~ team + at.home)
> baseballModel2
Bradley Terry model fit by glm.fit 

Call:  BTm(outcome = cbind(home.wins, away.wins), player1 = home.team, 
    player2 = away.team, formula = ~team + at.home, id = "team", 
    data = baseball)

Coefficients:
   teamBoston  teamCleveland    teamDetroit  teamMilwaukee   teamNew York    teamToronto        at.home  
       1.1438         0.7047         1.4754         1.6196         1.2813         1.3271         0.3023  

Degrees of Freedom: 42 Total (i.e. Null);  35 Residual
Null Deviance:	    78.02 
Residual Deviance: 38.64 	AIC: 137.1
> BTabilities(baseballModel2)
            ability      s.e.
Baltimore 0.0000000 0.0000000
Boston    1.1438027 0.3378422
Cleveland 0.7046945 0.3350014
Detroit   1.4753572 0.3445518
Milwaukee 1.6195550 0.3473653
New York  1.2813404 0.3404034
Toronto   1.3271104 0.3403222

AICは140.5から137.1と改善した。
$\alpha=0.3023$ と推定された。
つまり、 $\exp(\alpha)=1.352967$ となり、ホームチームの $V_i$ が1.35倍されるという解釈。
逆にアウェイチームの $V_i$ が1/1.35倍されるという解釈でも良い。

Bradley-Terryモデルの発展系

Bradley-Terryモデルの発展系はいっぱいあるみたい。
例えば、複数チームが対戦するもの、引き分けを考慮するもの、などなど。
私は、強さのパラメータが時系列に変化するモデルを探している。
昨年は調子悪かったけど、今年はだんだん調子が上がってきている、とかがわかると嬉しい。
2013年のJRSSのシリーズCに動的Bradley–Terryモデルというものが提案されているが、難しくてよくわからない。。。
http://wrap.warwick.ac.uk/54660/wrap.warwick.ac.uk
詳しい方教えてください。