##男子サッカーの世界ランキング
サッカーの世界最強国がどこか、日本は何番目なのか、というのはサッカーファンにとって尽きない話題です。
サッカーの国別ランキングとして最もポピュラーなのはFIFAランクですが、これには常に信頼性が微妙という批判がつきまとっています。というのも、算出方法が頻繁に変わっているからです。
過去の試合ごとのポイント制のシステムでは試合をすればするほど有利で、一時(1998年)は日本が9位(https://fifaranking.net/nations/jpn/) になるなどしてしまいました。それはあまりにも理不尽ということで、ポイントの計算方法が変更されたのですが、今度はヨーロッパ南米以外が著しく不利になる傾向が出てきました。
現在は将棋やチェスにも使われるイロレーティングに基づくシステムとなり、比較的信頼性が増していると言われています。それでも、イロレーティングは初期値に依存し収束に時間がかかるため、2018年夏から運用開始された現在の方式が適正値に収束しているとは言えないかもしれません。そこで、今回はイロレーティングやポイントに基づかない方式でサッカーの世界ランクを計算したいと思います。

##ベイズ統計モデリングによる強さの推定
StanとRでベイズ統計モデリングの将棋棋士の強さ推定を参考に、各国の強さを推定するモデルを作成しようと思います。このモデルは以前Jリーグにホームアドバンテージはあるかで書いたものと基本的に同じです。

各チーム$i$が固有の「強さ」$\mu_i$と「勝負ムラ」$\sigma_i>0$を持っていると仮定します。$\mu_i$が大きいほど強いチームで、$\sigma_i$が大きいほどパフォーマンスにムラがあるチームです。
ある試合において、ホームチーム$h$とアウェイチーム$a$が対戦する時、両チームの「パフォーマンス」$p$は、

\begin{align}
p_h =& \mu_h + \varepsilon_h + adv \\
p_a =& \mu_a + \varepsilon_a
\end{align}

で得られるものとします。ただし、$\varepsilon_i \sim N(0, \sigma_i)$とします。中立地での試合の場合、$adv$項はないものとします。
試合結果は、$p_h-p_a>1$のとき、ホームチームの勝ち。$p_a-p_h>1$のときアウェイチームの勝ち。それ以外では引き分けとします。つまり、パフォーマンスに差が1以上つけば勝ち、そうでなければ引き分けになるものとします。
これだけだと強さの基準が一意に定まらないため、$\mu \sim N(0, \sigma_\mu)$つまり強さの平均は0という仮定もおきます。

fifa.stan

data {
  int<lower=0> G_home;
  int<lower=0> G_away;
  int<lower=0> G_draw;
  int<lower=0> G_n_draw;
  int<lower=0> G_n_decided;
  int<lower=0> N;
  int<lower=1, upper=N> home_win[G_home];
  int<lower=1, upper=N> home_lose[G_home];
  int<lower=1, upper=N> away_win[G_away];
  int<lower=1, upper=N> away_lose[G_away];
  int<lower=1, upper=N> draw_home[G_draw];
  int<lower=1, upper=N> draw_away[G_draw];
  int<lower=1, upper=N> n_win[G_n_decided];
  int<lower=1, upper=N> n_lose[G_n_decided];
  int<lower=1, upper=N> n_draw_1[G_n_draw];
  int<lower=1, upper=N> n_draw_2[G_n_draw];
}

parameters {
  vector[N] mu; //各チームの強さ
  real<lower=0> s_mu; //muの標準偏差
  vector<lower=0>[N] s_pf; //勝負ムラ
  real<lower=0> home_adv; //ホームアドバンテージ
  
  ordered[2] performance_home[G_home];
  ordered[2] performance_away[G_away];
  real<lower=-1, upper=1> performance_draw_diff[G_draw];
  vector[G_draw] performance_draw_h;
  ordered[2] performance_n[G_n_decided];
  real<lower=-1, upper=1> performance_n_draw_diff[G_n_draw];
  vector[G_n_draw] performance_n_draw_1;
}

transformed parameters {
  vector[G_home] performance_home_h;
  vector[G_home] performance_home_a;
  vector[G_away] performance_away_h;
  vector[G_away] performance_away_a;
  vector[G_draw] performance_draw_a;
  vector[G_n_draw] performance_n_draw_2;
  
  for (g in 1:G_home){
    performance_home_h[g] = performance_home[g,2]+1-home_adv;
    performance_home_a[g] = performance_home[g,1];
  }
  for (g in 1:G_away){
    performance_away_a[g] = performance_away[g,2]+1+home_adv;
    performance_away_h[g] = performance_away[g,1];
  }
  for (g in 1:G_draw){
    performance_draw_a[g] = performance_draw_diff[g] + performance_draw_h[g] + home_adv;
  }
  for (g in 1:G_n_draw){
    performance_n_draw_2[g] = performance_n_draw_diff[g] + performance_n_draw_1[g];
  }
  
}

model {
  mu ~ normal(0, s_mu);
  for (g in 1:G_home){
    performance_home_h[g] ~ normal(mu[home_win[g]], s_pf[home_win[g]]);
    performance_home_a[g] ~ normal(mu[home_lose[g]], s_pf[home_lose[g]]);
  }
  for (g in 1:G_away){
    performance_away_a[g] ~ normal(mu[away_win[g]], s_pf[away_win[g]]);
    performance_away_h[g] ~ normal(mu[away_lose[g]], s_pf[away_lose[g]]);
  }
  for (g in 1:G_draw){
    performance_draw_h[g] ~ normal(mu[draw_home[g]], s_pf[draw_home[g]]);
    performance_draw_a[g] ~ normal(mu[draw_away[g]], s_pf[draw_away[g]]);
  }
  for (g in 1:G_n_decided){
    performance_n[g,1] ~ normal(mu[n_lose[g]], s_pf[n_lose[g]]);
    performance_n[g,2] ~ normal(mu[n_win[g]], s_pf[n_win[g]]);
  }
  for (g in 1:G_n_draw){
    performance_n_draw_1[g] ~ normal(mu[n_draw_1[g]], s_pf[n_draw_1[g]]);
    performance_n_draw_2[g] ~ normal(mu[n_draw_2[g]], s_pf[n_draw_2[g]]);
  }
}

fifa.r

library(rstan)
library(dplyr)

from_year <- 2016

results <- read.csv("./fifa_results.csv", stringsAsFactors = FALSE)
results$year <- substr(results$date,1,4)
results <- subset(results, year>=from_year)

team <- names(table(c(results$home_team, results$away_team)))
team <- data.frame(team)
team$id <- row_number(team)

results <- dplyr::inner_join(results, team, by=c("home_team"="team"))
results <- dplyr::inner_join(results, team, by=c("away_team"="team"))
results$Result <- (results$home_score > results$away_score)-(results$home_score < results$away_score)
results$Win <- results$id.x
results$Lose <- results$id.y
for(i in 1:nrow(results)){
  if(results$Result[i]==-1){
    results$Win[i] <- results$id.y[i]
    results$Lose[i] <- results$id.x[i]
  }
}

match <- results[c("Win","Lose","Result","neutral")]
match_n <- subset(match, neutral==TRUE)
match <- subset(match, neutral==FALSE)

N <- nrow(team)
match_draw <- subset(match, Result == 0)
match_home <- subset(match, Result == 1)
match_away <- subset(match, Result == -1)
match_n_draw <- subset(match_n, Result == 0)
match_n_decided <- subset(match_n, Result != 0)

G_draw <- nrow(match_draw)
G_home <- nrow(match_home)
G_away <- nrow(match_away)
G_n_draw <- nrow(match_n_draw)
G_n_decided <- nrow(match_n_decided)

data <- list(G_home=G_home, G_away=G_away, G_draw=G_draw, G_n_draw=G_n_draw, G_n_decided=G_n_decided,  N=N, home_win=match_home$Win, home_lose=match_home$Lose, away_win=match_away$Win, away_lose=match_away$Lose, draw_home=match_draw$Win, draw_away=match_draw$Lose, n_win=match_n_decided$Win, n_lose=match_n_decided$Lose, n_draw_1=match_n_draw$Win, n_draw_2=match_n_draw$Lose)

fit <- stan(file="./fifa.stan", data=data, seed=1234, init_r = 5, chains=4, iter=20000, warmup = 5000)
fit
ms <- rstan::extract(fit)
strength <- apply(ms$mu, MARGIN=2, FUN=mean)
strength_sd <- apply(ms$mu, MARGIN=2, FUN=sd)
p_sd <- apply(ms$s_pf, MARGIN=2, FUN=mean)
home_adv <- mean(ms$home_adv)
team$strength <- strength
team$strength_sd <- strength_sd
team$rank <- N - rank(team$strength) + 1
team$p_sd <- p_sd

match_count <- data.frame(table(c(results$home_team, results$away_team)))
colnames(match_count) <- c("team", "count")
team <- inner_join(team, match_count, by="team")

ranking <- team[order(team$rank),]
home_adv <- mean(ms$home_adv)

データはKaggleのInternational football results(2020年8月取得)を使い、2016年以降の試合結果を対象としました。
テクニカルには試合結果のパフォーマンスの不等式条件を実装するため、ordered型とlower,upper条件を駆使しています。もっとよい実装があればご教授ください。。

##結果
以下に、上位30か国の結果を示します。推定強さは、$\mu$の推定分布の平均値です。
FIFAランクには2020年8月16日現在のものを示します。

国	推定強さランキング	FIFAランキング
フランス	1	2
ブラジル	2	3
ベルギー	3	1
スペイン	4	8
オランダ	5	14
ドイツ	6	15
イングランド	7	4
アルゼンチン	8	9
ポルトガル	9	7
コロンビア	10	10
イタリア	11	13
クロアチア	12	6
メキシコ	13	11
スイス	14	12
ウルグアイ	15	5
イラン	16	33
ウクライナ	17	24
デンマーク	18	16
チリ	19	17
セルビア	20	29
ペルー	21	21
ベネズエラ	22	25
ポーランド	23	19
スウェーデン	24	17
韓国	25	40
セネガル	26	20
ボスニア・ヘルツェコビナ	27	49
ナイジェリア	28	31
モロッコ	29	43
日本	30	28

現行FIFAランクのうちtop20は推定強さランキングでもtop30入りしてきました。
FIFAランクはワールドカップの結果を重視するため、オランダやドイツといった2018年ロシアワールドカップで活躍できない/予選敗退した強豪国が低く評価される一方、ウルグアイやクロアチアといった活躍した国のランキングが高く出ています。
今回のランキングはすべての試合を平等に扱っているため、より実感/イメージに近いランキングになっていると思います。
一方、FIFAランクでは評価の低いイランが16位にランクインしました。これはイランがあまり他大陸の強豪相手に試合を行わ(え)ず、近隣の比較的弱い相手に白星を積み重ねた結果、強さ推定値が高く出るものと思われます。こういった方法でロンダリングが可能になってしまうのが、この統計モデリングによる方法の欠点と言えるかもしれません。

##まとめ
以上、サッカーの世界ランキングをベイズ統計モデリングで算出してみました。
今回はすべての試合を平等に扱ったので、もしこのアルゴリズムを運用するとなると、フルメンバーの試合と強化試合の差別化や、対象とする期間などに議論の余地が残るのかもしれません。