砂になった人

砂になった人の趣味や勉強についての備忘録

【R/Stan】項目反応理論でM-1を分析してみる

RとStanを使って、過去三回のM-1についての分析を項目反応理論を応用して分析してみました。目次は以下の通りです。

はじめに

こんにちは。最近Stan力上げないとヤバいなと思いネットサーフィンをしていたら、Stanを用いてM-1について分析する記事をいくつか発見しました(といっても2個目は去年のゼミで話題になりました)。

norimune.net

mtokuoka.net

めちゃくちゃ面白くないですか?Stanで柔軟なモデリングをすれば何か漫才についての色々なことが見えてきそうな気がします。

そこで、私自身もStanの練習を兼ねてM-1について分析してみようと思いました。

何を分析するか

この記事では、過去三回のM-1の成績データを使って、

  • スコアではわからない出演者の潜在的な漫才力
  • 審査員の採点に関する潜在的な厳しさ

をStanでモデルを組んで推定しようと思います。

では、なぜ最近3年のデータしか使わないといえば、最近の動向が気になるのはもちろんのことですが、この3年では審査員が固定されているので分析が楽ちんという点があります。この点については後ほど説明します。

どうやって分析するのか

潜在的なものを推定したいのでベイズを使うのですが、ここでは項目反応理論(Item Response Theory, IRT)という手法を応用しようと思います。私自身は政治学の畑の人ですが、卒論で心理学や教育学系の手法であるIRTを使うので、これは格好の機会だなと思い、やってみることにします(実際のところ、政治学でのIRT応用例は結構あります)。

項目反応理論を応用する

IRTについては、いろんなところで解説されつくされていると思うのでWikipediaでもいいので適当に見てください。

ja.wikipedia.org

私は主にこの本で勉強しています。

www.amazon.co.jp

通常のIRTのモデルは、以下のようなロジスティック回帰の時に用いられる形で書かれます(またはプロビット)。

$$ \text{Pr}(Y_i=1) = \frac{\exp(\beta_j\theta_i - \alpha_j)}{1 + \exp(\beta_j \theta_i - \alpha_j)} $$

一般に、$i$をあるテストの被験者とし、$j$をあるテストの項目(問題)とします。$\theta_i$は個人$i$の潜在的な能力とされ、$\alpha_j$*1は項目$j$の困難度、$\beta_j$は識別パラメータと呼ばれ、能力$\theta_i$が高い人とそうでない人とを分けている度合いとして捉えられます(要は曲線の傾き)。

しかし、今回扱うM-1のデータは従属変数$Y_i$がバイナリーでもなければ多値選択変数でもなく、得点というカウント変数です。そこで、カウントといえばポアソン分布の登場です。「上記のロジットをポアソンへと拡張すれば、カウント変数でもIRTできるんじゃね?」と考えました。実際に、Slapin & Proksch (2008) が提案したWordfishというモデルも同じような発想だと思います。つまり、

\begin{align} Y_i &\sim Poisson(\lambda_i),\\ \lambda_i &= \exp[\beta_j(\theta_i - \alpha_j)]. \end{align}

というモデルを考えます。$\beta_j(\theta_i - \alpha_j)$をlinear aggregator内に投入した形です。この発想が正しいかどうかは微妙ですが、一応テストしてみました。

f:id:ukuk1014:20210520194451p:plain f:id:ukuk1014:20210520194510p:plain

上が$\alpha$の値を動かしたときのもので、下が$\beta$の値を動かしたものです。通常の場合と同じような曲線の動きが見て取れると思います。なので、不安ですが概ね大丈夫かと。。。

それでは、モデルの説明をしたいと思います。まず、ここではM-1の分析をしたいので項目$j$を審査員$j$であるとして、困難度である$\alpha_j$をその審査員$j$の採点に関する厳しめ度合いということにします。また、M-1の文脈で$\beta$をどう表そうかと考えてましたが、良い日本語が浮かばなかったので、良い出演者$i$とそうでないものを分けるパラメータとして考えます。つまり、良い出演者にはしっかりといい点数が付き、逆にそうでない出演者には厳しい値が付くというような識別度合いですね。また、$\theta_i$を出演者$i$の潜在的な漫才力を表すとします。

さらに今回推定するデータは、時間が3年間に渡っているため、時間の影響考慮をするためにモデルに時間のパラメータ$\gamma_{t}$を加えることにします。

\begin{align} Y_i &\sim Poisson(\lambda_i),\\ \lambda_i &= \exp[\beta_j (\theta_i - \alpha_j) + \gamma_{t}]. \end{align}

とは言っても、時間に関するパラメータは実際のIRTにはほとんど入ることはないと思います。

それでは、早速分析していきましょう。

パッケージの読み込み

以下のパッケージを読み込みます。{tidybayes}はstan回した後の事後処理に非常に便利です。

library(rstan)
library(tidyverse)
library(tidylog)
library(tidybayes)

また、以下のStanに関する設定をします。これやると速くなる(?)。

options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

データの読み込み

今回用いるデータは、過去三回のM-1(2018 ~ 2020)の成績データです。M-1のサイトやら何やらを参照して作成しました(ミスありましたらスミマセン)。ダウンロードはこちら

# そのまま読み込むと文字化けするのでencodingをセットしてます。
df <- read.csv("data/m1_data.csv", encoding = "shift_jis") %>% 
  as_tibble() # tibbleが好きなので

中身を確認しましょう。

dim(df)
head(df)
tail(df)
[1] 161   4
# A tibble: 6 x 4
   time group      referee    score
  <int> <chr>      <chr>      <int>
1    14 ギャロップ オール巨人    87
2    14 ギャロップ 松本人志      86
3    14 ギャロップ 上沼恵美子    89
4    14 ギャロップ 塙宜之        89
5    14 ギャロップ 富澤たけし    87
6    14 ギャロップ 立川志らく    86
# A tibble: 6 x 4
   time group          referee    score
  <int> <chr>          <chr>      <int>
1    16 東京ホテイソン 松本人志      86
2    16 東京ホテイソン 上沼恵美子    92
3    16 東京ホテイソン 塙宜之        85
4    16 東京ホテイソン 富澤たけし    91
5    16 東京ホテイソン 立川志らく    89
6    16 東京ホテイソン 礼二          88

ここでは、2回連続で出場している組も何組かあり(e.g. 和牛)、分析の際に重複すると面倒くさいので、その組については一番新しい出場データのみを適用しています。今回やりたいことは、審査員$j$をIRTでいうところのItemとして捉えることで各出演者$i$の潜在的な能力$\theta_i$(潜在的な漫才力)をIRTで推定することです。この時、時系列の分析にしてしまうと欠損値の扱いが大変になってしまったりするし、単純にデータを集めてしまうと重複する出演者が出てきてしまうので、それを防ぐために上記の収集ルールにしています。

それでは、このデータを用いてまずは簡単な可視化をします。まず、出演者は全部で何組いるのか確認します。

length(unique(df$group))
[1] 23

通常、M-1決勝は10組ずつなので23組はちょっとおかしいですが、上記のルールを適用したからこうなっているだけなので気にしないようにします。また、審査員も確認します。

length(unique(df$referee))
[1] 7

先ほども説明した通り、過去三回は審査員が全員同じなので、データとしては非常に分析しやすいです。

出演者ごとの総得点を見てみます。ここで、いきなりボンッと得点に関するグラフを貼ってしまうと、どの組も600点以上を取っているのであまり得点の差が見えないようなグラフができてしまいます。そこで、全出演者の総得点から最低総得点を引いて+1した値を採用します。

df %>% 
  group_by(group) %>% 
  summarise(sum = sum(score)) %>%
  mutate(sum = sum - min(sum) + 1) %>% 
  ggplot(aes(x = fct_reorder(group, sum), y = sum)) + # fct_reorder(factor, value)で大きい順に並び変える
  geom_bar(stat = "identity", fill = "#00b3b3") +
  coord_flip() +
  ylab("Total Score (Adjusted)") + xlab("") +
  theme_light()

f:id:ukuk1014:20210520184413p:plain

14 ~ 15回の最終決戦チームが結構高得点を取ってますが、前回は何か微妙な感じですね(笑)。次に、審査員ごとの総得点を見てみます。

df %>% 
  group_by(referee) %>% 
  summarise(sum = sum(score)) %>% 
  mutate(sum = sum - min(sum) + 1) %>% 
  ggplot(aes(x = fct_reorder(referee, sum), y = sum)) + 
  geom_bar(stat = "identity", fill = "#00b3b3") +
  coord_flip() +
  ylab("Total Score (Adjusted)") + xlab("") +
  theme_light()

f:id:ukuk1014:20210520184346p:plain

ここから単純に見れば、巨人さんは厳しめに得点を付けていて上沼さんは甘めに付けているということが何となくですがわかります。しかし、本当にそうなのかはわかりません。

項目反応理論(IRT)をやってみる

モデル

それではIRTの実装に入りたいと思います。先述の通り、ここで行うIRTは従属変数がバイナリーというわけではなく、また多値の選択変数でもないただのカウント変数ため、いわゆるロジットやプロビットで行われる通常の2Pモデルではなく、ポアソンへと応用したモデルを扱います。generative modelは、

\begin{align} Y_i &\sim Poisson(\lambda_i),\\ \lambda_i &= \exp[\beta_j (\theta_i - \alpha_j) + \gamma_{t}]. \end{align}

です。

データの整形

IRTを行う際には、IRTのための行列(以下、IRT行列)を作るのが慣例というか恒例というか義務(?)なので、早速作りたいと思います。IRT行列とは、行に回答者$i$・列に項目$j$、要素は正解or不正解のような形です。ここでの事例に合わせるならば、行は出演者$i$・列は審査員$j$、要素は得点です。以下で加工してみます。

df %>% 
  select(group, referee, score) %>% 
  mutate(score = score - min(df$score)) %>% 
  pivot_wider(names_from = referee, values_from = score) %>% # wider型へ変換
  as.matrix() -> m_irt

得点に関してですが、90点とかだと大きすぎるので、各得点から最低点を引いた値を採用しています。諸々残りをやります。

gname <- m_irt[, 1] # 1列目(出演者名)を抽出
m_irt <- m_irt[, -1] # 1列目を削除
rownames(m_irt) <- gname # 行ラベルとして抽出した出演者名を適用
rname <- colnames(m_irt) # 後ほど使うので審査員名も抽出

できた行列を見てみます。

m_irt
オール巨人 松本人志 上沼恵美子 塙宜之 富澤たけし 立川志らく 礼二
ギャロップ 7 6 9 9 7 6 10
ジャルジャル 13 12 8 13 10 19 13
スーパーマラドーナ 7 5 9 9 9 8 10
トムブラウン 7 11 6 13 9 17 10
ミキ 10 8 18 10 10 9 13
ゆにばーす 4 0 4 2 6 7 11
霜降り明星 13 14 17 18 11 13 16
かまいたち 13 15 15 15 13 15 14
からし蓮根 13 10 14 10 10 9 13
すゑひろがりず 12 9 12 11 10 12 11
ぺこぱ 13 14 16 14 14 11 12
ミルクボーイ 17 17 18 19 17 17 16
和牛 12 12 12 16 11 16 13
アキナ 9 5 12 7 8 10 11
インディアンス 9 10 13 5 9 9 10
ウェストランド 8 10 12 5 11 6 10
おいでやすこが 12 15 14 13 13 16 15
オズワルド 8 8 12 15 11 13 15
ニューヨーク 8 12 14 13 13 11 11
マヂカルラブリー 8 13 14 14 14 10 16
錦鯉 7 9 13 15 12 15 12
見取り図 11 11 15 13 12 13 13
東京ホテイソン 6 6 12 5 11 9 8

出演者の順番はあまり気にしないでください。

次に、Stanへ渡すためのデータを作らなくてはなりません。Stanくんはlist型のデータではないと読み込んでくれないのがちょっとダルい点です。

Y <- as.numeric(m_irt) # 得点
N <- length(Y) # observation
i <- rep(1:dim(m_irt)[1], times = dim(m_irt)[2]) # 出演者の観察数
j <- rep(1:dim(m_irt)[2], each = dim(m_irt)[1]) #審査員の観察数
I <- max(i) # 出演者数
J <- max(j) # 審査員数

# 時間の影響考慮もしたいので抜き出します

df %>% 
  arrange(referee) -> df2

t <- df2$time - min(df2$time) + 1 # 時間の観察数
T <- max(t) # 時間数

data <- list(I = I, J = J, N = N, T = T,  
             i = i, j = j, t = t, Y = Y)

Stanの実装

次に、Stanを書きます。今回のStanコードは以下の通りです(pois_irt.stanという名前で保存しました)。

また、priorとして

$$ \begin{align} \alpha &\sim \mathcal{N}(0, 1), \\\ \beta &\sim \mathcal{N}(0, 10), \\\ \gamma &\sim \mathcal{N}(0, 10), \\\ \theta &\sim \mathcal{N}(0, 1). \end{align} $$

を与えます。

data{
    int<lower=1> N; // # of observations
    int<lower=1> I; // # of groups
    int<lower=1> J; // # of referees
    int<lower=1> T; // # of time
    int<lower=1, upper=I> i[N]; // groups observation
    int<lower=1, upper=J> j[N]; // referees observation
    int<lower=1, upper=T> t[N]; // time observation
    int<lower=0> Y[N]; // score
}

parameters {
    vector[T] gamma; //time parameter
    vector[J] alpha; //difficulty
    vector[J] beta; //discrimination
    vector[I] theta; //latent ability
}

model {
    for (n in 1:N){
        Y[n] ~ poisson(exp(beta[j[n]] * (theta[i[n]] - alpha[j[n]]) + gamma[t[n]])); 
    }

    theta ~ std_normal(); // for standardizing
    alpha ~ std_normal();
    beta ~ normal(0, 10);
    gamma ~ normal(0, 10);
}

できたので早速Stanを回しましょう。最初にモデルを読み込みます。

model <- stan_model("stan/pois_irt.stan")

警告が出る時がありますが無視してかまいません。キック用のコードは以下の通りです。

fit <- sampling(model, 
                data = data, 
                iter = 20000, 
                warmup = 5000, 
                chains = 1,
                seed = 1, 
                thin = 10, 
                control = list(adapt_delta = 0.99,  # controlなしでやったら警告が出たので入れてます
                               max_treedepth = 15))

収束しているかチェックしましょう。

all(summary(fit)$summary[, "Rhat"] <= 1.10, na.rm = TRUE)
[1] TRUE

TRUEが返されているので、すべてのパラメータに関するRhatが1.1以下であったことがわかり、収束していると判断します。上記の便利な判定コードは、以下のslideshareを参照しました。めちゃくちゃ便利ですね。

www.slideshare.net

結果

次に、各パラメータの事後分布を可視化します。まずは、各出演者の潜在的な漫才力$\theta_i$を見てみましょう。この値が高ければ高いほど潜在的な漫才力が高いことを表します。また、推定値の線は細いものから順に95%、90%、80%の信用区間を示し、黒丸は事後分布の中央値です。

group <- tibble(group = gname, n = 1:length(gname)) # 出演者の名前を投入するためにデータフレームを作ります

fit %>% 
  spread_draws(theta[n]) %>% # stanの結果をtidyに変換
  median_qi(.width = c(.95, .9, .8)) %>% # 信用区間やら諸々
  left_join(group, by = "n") %>% 
  ggplot(aes(x = theta, y = fct_reorder(group, theta),
             xmin = .lower, xmax = .upper, 
             size = -.width)) +
  geom_pointinterval(interval_size_range = c(0.4, 1.5)) +
  theme_light() +
  xlab("Estimated Latent Score") + ylab("")

f:id:ukuk1014:20210521000528p:plain

霜降り明星とミルクボーイ強いですねぇ~。得点だけで見たらミルクボーイが圧倒的でしたが、いざ分析してみると霜降りの方が強いとは。個人的にはトムブラウンが割と上位の方にいて嬉しいですね~。でんじゃらすじーさんを見ている感じがしてとても良い。

個人$i$だけでなく項目$j$の推定も見てみましょう。まずは、審査員の厳しめ度合いを表す$\alpha_j$から見てみます。

referee <- tibble(referee = rname, n = 1:length(rname)) # 審査員の名前のデータフレーム

fit %>% 
  spread_draws(alpha[n]) %>%
  median_qi(.width = c(.95, .9, .8)) %>% 
  left_join(referee, by = "n") %>% 
  ggplot(aes(x = alpha, y = fct_reorder(referee, alpha),
             xmin = .lower, xmax = .upper, 
             size = -.width)) +
  geom_pointinterval(interval_size_range = c(0.4, 1.5)) +
  theme_light() +
  xlab("Estimated Difficulties") + ylab("")

f:id:ukuk1014:20210521000547p:plain

結果としては言われている通り巨人さんでが厳しめで、上沼さんは甘めなのでしょうか。しかし、信用区間がかなり広く不確実性が大きいので何とも言えないです。

最後に$\beta$の値を見てみます。解釈が難しいパラメータですが、、、、

fit %>% 
  spread_draws(beta[n]) %>%
  median_qi(.width = c(.95, .9, .8)) %>% 
  left_join(referee, by = "n") %>% 
  ggplot(aes(x = beta, y = fct_reorder(referee, beta),
             xmin = .lower, xmax = .upper, 
             size = -.width)) +
  geom_pointinterval(interval_size_range = c(0.4, 1.5)) +
  theme_light() +
  xlab("Estimated Saliency") + ylab("")

f:id:ukuk1014:20210521000615p:plain

これ、何と呼んでいいかわからず横軸をSaliency(政治学では$\beta$はsailiency parameterと呼ばれることが多い)にしてしまいました。塙さんが一番大きな値を獲得してますね。能力のある芸人にはしっかりと点数を付けて、そうでない人には厳しいといった評価を与えているという風に解釈できますね。そう考えると、礼二さんが一番可もなく不可もなくのような採点をしているのでしょうか。

初めてお笑いの分析をしてみましたが結構面白いですね。今回やったことがいろんな方面に活かせたら結構面白いですね~。

*1:政治学の文脈では$\beta$がdiscriminationで$\alpha$がdifficultyとして扱われることが多いのですが、心理学では逆のことが多いようです。