ベイズ統計モデリングでプレミアリーグ各チームの強さを検証する

2021 - 2022シーズンのサッカープレミアリーグの結果を基に、各チームの強さをスコア化し、信頼区間を求めて順位の妥当性を検証してみます。
こちらのサイトのデータを使用しました。
https://www.football-data.co.uk/englandm.php

### モデル
単純な回帰モデルで検証します。
強さを測る指標は勝ち点$Y$
各チームのシュート数$X_1$、被シュート数$X_2$の時、各効果量を$b_1$、$b_2$とした時、切片$a$を加えた以下のモデルで強さを算出します。

$Y = a + b_1X_1 + b_1X_2$

$a$、$b_1$、$b_2$は全チーム共通の値です。

### 検証

#### 使用データ
引用したデータを整形し、各チーム毎にシュート数 (Shots)、被シュート数 (GivenShots)、勝ち点 (Points)をまとめます。
分析結果の振り返りのため、各チームのカップ戦の出場権 (CL, EL, ECL)、プレミアリーグ残留 (Keep)、下部リーグ降格 (Relegate)の条件を付したstatusも加えました。

```r
> d_team
# A tibble: 20 x 4
   Team           Shots GivenShots Points status  
   <chr>          <int>      <int>  <dbl> <chr>   
 1 Man City         715        237     93 CL      
 2 Liverpool        730        296     92 CL      
 3 Chelsea          592        341     74 CL      
 4 Tottenham        491        481     71 CL      
 5 Arsenal          588        425     69 EL      
 6 Man United       508        511     58 EL      
 7 West Ham         433        548     56 ECL     
 8 Leicester        435        558     52 Keep    
 9 Brighton         490        476     51 Keep    
10 Wolves           402        511     51 Keep    
11 Newcastle        448        520     49 Keep    
12 Crystal Palace   412        443     48 Keep    
13 Brentford        442        502     46 Keep    
14 Aston Villa      454        477     45 Keep    
15 Southampton      482        497     40 Keep    
16 Everton          438        522     39 Keep    
17 Leeds            485        582     38 Keep    
18 Burnley          407        603     35 Relegate
19 Watford          396        561     23 Relegate
20 Norwich          374        631     22 Relegate
```
### stanモデリング

stanコードは以下の通りです。

```stan:hierarchical_bayesian_model.stan
data {
  int<lower=0> N;
  int<lower=0> S[N];
  int<lower=0> G_S[N];
  int<lower=0> P[N];
}

parameters {
  real<lower=0> sigma;
  real a;
  real b1;
  real b2;
}

transformed parameters{
  real mu[N];
  for (n in 1:N) {
    mu[n] = a + b1*S[n] + b2*G_S[n];
  }
}

model {
  for (n in 1:N) {
    P[n] ~ normal(mu[n], sigma);
  }
}

```

インプットするデータはチーム数$N$、シュート数$S$、被シュート数$G\\_S$、勝ち点$P$
推定するパラメータは切片$a$、シュート数効果$b_1$、被シュート数効果$b_2$、勝ち点のばらつき$sigma$です。

#### 実行結果

```r
dat <- list(N = nrow(d_team), S=d_team$Shots, G_S=d_team$GivenShots, P=d_team$Points)
fit <- stan(file = 'hierarchical_bayesian_model.stan', data = dat, seed = 1234)

> fit
Inference for Stan model: hierarchical_bayesian_model.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

         mean se_mean    sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
sigma    9.26    0.06  1.82   6.52   7.97   8.98  10.24  13.65   919    1
a       48.18    1.43 44.72 -41.40  18.74  46.95  77.00 134.17   978    1
b1       0.10    0.00  0.05   0.00   0.07   0.10   0.13   0.19  1030    1
b2      -0.09    0.00  0.05  -0.18  -0.12  -0.09  -0.06   0.01  1021    1
mu[1]   96.23    0.10  5.81  84.87  92.41  96.23 100.04 107.88  3579    1
mu[2]   92.54    0.12  5.96  80.60  88.82  92.42  96.54 104.54  2544    1
mu[3]   75.38    0.08  3.90  67.91  72.78  75.35  77.95  83.41  2244    1
mu[4]   53.52    0.03  2.10  49.29  52.16  53.49  54.85  57.77  3901    1
mu[5]   67.70    0.07  3.41  60.77  65.55  67.68  69.95  74.54  2146    1
mu[6]   52.54    0.07  2.99  46.47  50.62  52.57  54.45  58.46  1728    1
mu[7]   42.13    0.04  2.47  37.23  40.52  42.16  43.75  47.03  4012    1
mu[8]   41.45    0.05  2.64  36.12  39.74  41.50  43.20  46.58  3357    1
mu[9]   53.86    0.04  2.12  49.62  52.47  53.84  55.21  58.15  3420    1
mu[10]  42.37    0.10  3.69  35.19  39.97  42.32  44.71  49.90  1345    1
mu[11]  46.00    0.04  2.25  41.53  44.54  46.01  47.45  50.67  3815    1
mu[12]  49.24    0.18  5.81  37.73  45.42  49.04  53.14  60.58  1038    1
mu[13]  46.99    0.06  2.57  41.91  45.35  46.93  48.62  52.36  1920    1
mu[14]  50.31    0.07  2.85  44.71  48.45  50.22  52.14  56.09  1480    1
mu[15]  51.26    0.03  2.12  46.93  49.87  51.30  52.62  55.52  3944    1
mu[16]  44.87    0.04  2.37  40.19  43.34  44.84  46.34  49.77  3030    1
mu[17]  44.17    0.14  4.89  34.25  41.02  44.27  47.42  53.87  1266    1
mu[18]  34.85    0.07  3.42  27.91  32.61  34.94  37.06  41.54  2632    1
mu[19]  37.45    0.05  2.89  31.78  35.58  37.45  39.23  43.36  3062    1
mu[20]  29.25    0.07  3.79  21.50  26.71  29.28  31.71  36.62  3134    1
lp__   -51.42    0.05  1.58 -55.43 -52.27 -51.05 -50.23 -49.45  1043    1

```

#### 可視化

分析結果を可視化します。

**パラメータb1, b2**
シュート数に関するパラメータb1は95%s信頼区間で0をまたぐことなく正の値を取っています。
一方被シュート数に関するパラメータb2は95%s信頼区間で0をまたいでいます。
シュート数が勝ち点にポジティブに与える影響と被シュート数が勝ち点にネガティブに与える影響では、前者のほうが信頼性が高そうです。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1031677/75d9b76e-2e15-d71c-4646-fec3275dd718.png)

**勝ち点**

各チームの勝ち点推定値の中央値、95%信頼区間を可視化します。
マンチェスターシティ、リヴァプールの上位2チームは文句なしに見えます。
惜しくもCL出場権を逃したアーセナルは実力的には4位に食い込む力がありそうですが、他の要因もあり順位を落としてしまったのでしょうか。
9位のブライトンも推定値ではカップ戦出場権を勝ち得ていてもおかしくないですね。シーズン通した試合の戦い方や運用の仕方が肝となりそうです。
降格してしまった下位3チームは仕方ないといったところでしょうか。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1031677/63d9333f-1fe7-c0ee-db45-27f23f75d7a7.png)


### まとめ

シュート数と被シュート数だけでもおおむね納得感のある結果が出ました。しかしカップ戦がかかるくらい上位の戦いになると、それ以外の要因も大きくかかわってきそうです。今年のプレミアリーグも楽しみですね。

### 参考資料
サッカーデータまとめサイト
https://www.football-data.co.uk/englandm.php

アヒル本
[StanとRでベイズ統計モデリング](https://www.kyoritsu-pub.co.jp/book/b10003786.html)