M1の審査員の採点が適正だったのか異常検知により分析する

先日開かれたM1において一部の審査員の採点に問題があったのかどうかについて話題になっているそうです。
たまたま私が最近異常検知という分析手法について勉強しているところだったので、これを使って審査員の採点が異常であったかどうかについて分析してみたいと思います。
なお得られる標本データが少ないので流行りの機械学習による分析ではなく、純粋統計に近い分析になっております。

データセット

M1の決勝における各審査員の各漫才コンビに対する採点を分析対象のデータセットとします。
審査員は7人、漫才コンビは10組の計70のレコードです。

#R
> library(tidyverse)
> m1 <- read_csv(
	       "m1_2018data.csv",
	       col_names=c("combi", "order", "judge", "cat", "score"),
	       skip=1) %>%
    select(combi, judge, score)
> m1
# A tibble: 70 x 3
   combi      judge  score
   <chr>      <chr>  <int>
 1 かまいたち 巨人      89
 2 かまいたち 礼二      92
 3 かまいたち 塙        92
 4 かまいたち 志らく    88
 5 かまいたち 富澤      91
 6 かまいたち 松本      90
 7 かまいたち 上沼      94
 8 ギャロップ 巨人      87
 9 ギャロップ 礼二      90
10 ギャロップ 塙        89
# ... with 60 more rows

まずはデータの概観をざっくり確認しておきます。

#R
> m1 %>%
    ggplot(aes(combi, score)) +
      stat_summary(fun.data = function(x)mean_se(x, mult=2), size = 1) +
      geom_point() +
      geom_text(aes(label = judge), hjust = -1.4)
M1決勝採点データの要約

大きい黒丸がコンビごとの平均点で、上下の線が2\sigmaを表します。
2\sigmaから外れれば平均からは大きくズレたと言えそうなので、確かに一部の審査員の採点が平均から外れているようにも見えます。
しかし個別の採点を見るだけで全体を異常であると判断するのは危険です。
審査員のすべての採点を同時に評価して判断をくださなければなりません。

問題設定

知りたいことはこの平均からのズレがある審査員について大きく偏っているかどうかです。
よって我々の問題を次のように設定することにしましょう。

それぞれの漫才コンビを変数とし、それぞれの審査員の採点を測定とします。
つまり審査員が各漫才コンビの真の点数を測定していると考えます。
測定には誤差が入ります。
ある測定において多くの変数について誤差があまりに大きく、通常の状態では発現しない水準を上回っていれば、その測定は異常であったと判断します。

また分析を単純にするためにそれぞれの得点は正規分布に従うこととします。

採点の異常検知

ここでホテリング法による多変量異常検知を実施します。
ホテリング法ではデータが従う分布から負の対数尤度を計算してこの数字を異常値とします。
そしてこの統計量が従う分布から発生確率を求め、それが非常に小さい確率であれば何か異常事態が発生したのであろうと判定します。

分布は正規分布に従うとしたので確率密度関数は下のようになります。


\displaystyle p(x)=\frac{|\Sigma|^{-1/2}}{(2\pi)^{M/2}}\exp\left\{-\frac{1}{2}(x-\mu)^\mathrm T\Sigma^{-1}(x-\mu)\right\}

ただし\muは母平均、\Sigmaは母共分散行列、Mは変数の数(今回はコンビ数)です。

教科書通りにやるのであれば標本平均を\muの、標本共分散行列を\Sigmaの最尤推定値として採用し尤度を計算します。


\displaystyle\hat\mu=\frac1N\sum_i x_i
\displaystyle\hat\Sigma=\frac1N\sum_i(x_i-\hat\mu)(x_i-\hat\mu)^\mathrm T

しかしここで問題が発生します。
今回のデータは変数よりも測定の方が少ないため、この場合には標本共分散行列は特異になり逆行列を持ちません。
よって教科書通りに多変量正規分布の共分散行列を最尤推定することができません。

変数のほうが多い分散共分散行列の推定についてWikipediaを見ると小難しいことがたくさん書いてあり、なんとなく色々な課題があって単純な答えが存在しない雰囲気です。
こういう場合は付け焼刃で仕掛けるよりは推奨されているパッケージを無心に使うのが良い気がします。
よってWikipediaで推奨されているcorpcorというRのライブラリを使って分散共分散行列(の逆行列)の推定値を得ることにしました。
不安が無いこともないので詳しい方はアドバイスいただけると助かります。

corpcor: Efficient Estimation of Covariance and (Partial) Correlation

#R
# scoreをコンビごとに平均値からの差分にしてから行列へと変換
> m1_val <- m1 %>%
    group_by(combi) %>%
    mutate(score = score - mean(score)) %>%
    ungroup %>%
    spread(combi, score)
> m1_mat <- m1_val %>%
    select(-judge) %>%
    as.matrix
# corpcorを使って共分散行列の推定
> install.packages("corpcor")
> library(corpcor)
> sigma_inv <- invcov.shrink(m1_mat)

ともかくも分散共分散行列の逆行列が推定できれば測定xの異常値a(x)は負の対数尤度で計算することができます。


a(x) = (x-\hat\mu)^\mathrm T\hat\Sigma^{-1}(x-\hat\mu)

ここで\hat\muは変数ごとの標本平均、\hat\Sigma^{-1}は推定した分散共分散行列の逆行列です。

#R
# 異常値を負の対数尤度によって計算
> a <- rowSums((m1_mat %*% sigma_inv) * m1_mat
> m1_val %>% mutate(anomaly = a) %>% select(judge, anomaly)
# A tibble: 7 x 2
  judge    anomaly
  <chr>     <dbl>
1 上沼       9.39
2 塙         4.88
3 富澤       4.41
4 巨人       4.22
5 志らく     9.67
6 松本       4.95
7 礼二       7.09
> m1_val %>%
     mutate(annomaly = a) %>%
     ggplot(aes(judge, anomaly)) +
       geom_col()
審査員ごとの採点の異常度

採点の異常度の一位は志らく、二位は上沼、少し開けて三位が礼二になりました。
残りの四人は同程度の異常度で、最も異常度が低かったのは巨人でした。
審査員のクラスターは一位二位のトンガリグループ、四位以下の普通グループ、礼二だけその中間に浮いている、と分類できないこともないです。

ただこれらが異常な水準と言えるかというとそうではなく、F分布を用いて確率計算すると一位の志らくですら異常確率12%程度なので異常な測定であるという仮説は棄却されます。

なお、本当は母集団の推定を行うときには異常度の評価対象の測定値を省くべきなのですが面倒くさくて省略しました。すいません。
サンプル数が少ないので実際には大きく影響するのかもしれません。

まとめ

異常検知により審査員の採点が異常であったかを統計的に推定しました。
その結果、異常であるという仮説は棄却されました。
なお、私は番組を見ていないので採点か妥当であったかは全くわかりません。すいません。

コメントを残す