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

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

データセット

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

#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)

大きい黒丸がコンビごとの平均点で、上下の線が $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()