廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

[RStan]多項ロジスティックモデルで μ's とAqours の人気の差を調べる

背景の整理

μ's とAqours の人気の差(驚異のアニヲタ社会復帰への道) のデータを使わせていただきます。

あるアニメショップでラブライブのキャラの人気投票をしたとき、μ's のメンバーが、Aqours のメンバーより全体的に上位だったそうです。

そこで 2 グループの人気はどれくらいの差かを考えます。

  • ラブライブとは架空のアイドルグループの奮闘と成長を描いた作品で、アニメ版やゲーム版がある
  • μ's(ミューズ)はラブライブの主人公たちの所属するアイドルグループ
  • Aqours(アクア)はラブライブの派生作品の主人公たちが所属するアイドルグループ
  • μ's を主人公とした作品群のほうが先に制作され、Aqours を主人公とした作品群が後から制作された

参考:ラブライブ! - Wikipedia

モデル

キャラクター i の得票数 y_i はパラメータ p_i の多項分布に従うとします。

ここで p_i は、

\displaystyle p_i = \frac{\exp(gx_i+b_i)}{\sum\exp(gx_i+b_i)}

で、x_i はキャラクター i が μ's に所属しているとき 1、さもなくば 0 の値を取る変数とします。

投票で  gx_i+b_i がキャラクター i が選ばれる「強さ」を表していると解釈します。

係数 g はグループが人気に与える効果の強さを表すと解釈できます。

係数 b_i は個人ごとに異なる人気度を表すと解釈し、平均 0、標準偏差 σ の正規分布に従うと仮定します。

\displaystyle {\rm softmax}(a_i) = \frac{\exp(a_i)}{\sum\exp(a_i)}

はソフトマックス関数とよばれ、p_i の総和が 1 になるようにするために使用しています。

さて、ソフトマックス関数は変数を平行移動しても値は同じになります。

たとえば、

\displaystyle {\rm softmax}(a_i + 1) =  {\rm softmax}(a_i)

ということです。

キャラクターの人気度は相対的なものなので基準をどこかで決めてやらないといけません。

今回は b_i に平均 0 の正規分布を仮定したので、基準は 0 に決まります。

モデルのパラメータから確率の比を求めることができます。

μ's に所属していたときの得票する確率は、

\displaystyle \frac{\exp(g+b_i)}{\sum\exp(gx_i+b_i)}

μ's に所属していないときの得票する確率は、

\displaystyle \frac{\exp(b_i)}{\sum\exp(gx_i+b_i)}

なので比を取ると、

\exp(g)

となります。

\exp(g) はμ's に所属することで、得票する確率がどれだけ上がるかの指標になります。

また仮に所属するグループが同じだった場合の、キャラクター i とキャラクター j の得票する確率の比も同様に求めることができます。

\exp(b_i)/\exp(b_j)

b_i は平均が 0 になるよう基準化されていることを仮定したので、b_j = 0 として、\exp(b_i) を、得票する確率が平均にくらべどれだけ高いかの指標とすることにします。

質問:この \exp(g)\exp(b_i) をリスク比と呼んでもいいのでしょうか?

解析対象

こんなデータです。

f:id:abrahamcow:20170624124125p:plain

横軸が得票数。

結果

Stan で MCMC を2000回して収束しました。

サンプリングは一瞬で終わります。

トレースプロットはこんな感じです。

f:id:abrahamcow:20170624152823p:plain

当てはまりを見るために予測値と実測値をプロットしてみます。

f:id:abrahamcow:20170624203448p:plain

赤が実測値、青が予測値、線分は95%予測区間です。

当てはまりは良さそうです。

次にグループ効果 g の事後分布を見ます。

f:id:abrahamcow:20170624153131p:plain

青い実線が事後期待値、点線が95%信用区間です。

μ's は Aqours にくらべ 2.4倍くらい人気であるといえそうです。

各キャラクターの、グループ効果を補正した人気度 \exp(b_i) はこんな感じです。

f:id:abrahamcow:20170624153511p:plain

丸が事後期待値、線分が95%信用区間です。

矢澤にこ(にこにー)はぼくでも聞いたことあるキャラクターなので、にこにーが一位に来るかな―と思っていたのですが、そうはなりませんでした。

それでもにこにーのベイズ信頼区間の下限は 1 から十分に離れており、平均的なキャラクターよりは人気であると、確信を持って言えそうです。

もし μ's と Aqours が同時期に活動を開始していたら、一番人気は津島善子になっていたのかもしれません。

この結果が背景知識と矛盾しないかは、ラブライブに詳しい方の意見を伺いたいです。

付録

R のコードと Stan のコードを貼り付けます。

#R
library(rstan)
library(cowplot)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
pop <-"name vote group
矢澤にこ 384 muse
西木野真姫 384 muse
南ことり 370 muse
東條希 336 muse
綾瀬絵里 283 muse
園田海未 263 muse
小泉花陽 251 muse
星空凛 226 muse
津島善子 190 aqours
黒澤ルビィ 171 aqours
国木田花丸 155 aqours
渡辺曜 141 aqours
高坂穂乃果 128 muse
松浦果南 110 aqours
黒澤ダイヤ 103 aqours
桜内梨子 93 aqours
小原鞠莉 80 aqours
高海千歌 74 aqours"
pop_dat <-read.table(text = pop,header = TRUE)
ggplot(pop_dat,aes(x=reorder(name,vote),y=vote,fill=group))+
  geom_col()+
  coord_flip()+
  xlab("")+
  theme_cowplot(font_family = "HiraKakuProN-W3")

pop_dat4stan <- list(Member=nrow(pop_dat),Vote=pop_dat$vote,group=as.integer(pop_dat$group=="muse"))

lovelive <-stan_model("lovelive.stan")
smp_pop <-sampling(lovelive,pop_dat4stan)
traceplot(smp_pop,pars=c("b","g","sigma"))
all(summary(smp_pop)$summary[,"Rhat"]<1.1)
ex_pop <-extract(smp_pop)

pop_dat$pred <- get_posterior_mean(smp_pop,pars="pred_boot")[,"mean-all chains"]
pop_dat$lower <- apply(ex_pop$pred_boot,2,function(x)quantile(x,probs = 0.025))
pop_dat$upper <- apply(ex_pop$pred_boot,2,function(x)quantile(x,probs = 0.975))

ggplot(pop_dat,aes(x=reorder(name,vote)))+
  geom_linerange(aes(ymin=lower,ymax=upper,colour="pred"))+
  geom_point(aes(y=pred,colour="pred"),size=3,alpha=0.5)+
  geom_point(aes(y=vote,colour="obs"),size=3,alpha=0.5)+
  ylim(c(0,NA))+
  coord_flip()+
  xlab("")+
  theme_cowplot(font_family = "HiraKakuProN-W3")


hist(exp(ex_pop$g),breaks = "scott",xlab = "exp(g)",main = "")
abline(v=quantile(exp(ex_pop$g),probs = 0.025),col="royalblue",lwd=2,lty=2)
abline(v=quantile(exp(ex_pop$g),probs = 0.975),col="royalblue",lwd=2,lty=2)
abline(v=mean(exp(ex_pop$g)),col="royalblue",lwd=2)

mean(exp(ex_pop$g))

est <-data.frame(name=pop_dat$name,
                 group=pop_dat$group,
                 RR=apply(ex_pop$b,2,function(x)mean(exp(x))),
                 lower=apply(ex_pop$b,2,function(x)quantile(exp(x),probs = 0.025)),
                 upper=apply(ex_pop$b,2,function(x)quantile(exp(x),probs = 0.975)))

ggplot(est,aes(x=reorder(name,RR),y=RR,colour=group))+
  geom_pointrange(aes(ymin=lower,ymax=upper))+
  geom_hline(yintercept=1,linetype=2)+
  coord_flip()+
  xlab("")+ylab("exp(b)")+
  theme_cowplot(font_family = "HiraKakuProN-W3")
#Stan
data{
  int<lower=1> Member;
  int<lower=1> Vote[Member];
  int<lower=0, upper=1> group[Member];
}
parameters{
  real b[Member];
  real g;
  real<lower=0> sigma;
}
transformed parameters{
  vector[Member] alpha;
  for(i in 1:Member){
    alpha[i] = g*group[i] + b[i];
  }
}
model{
  b ~ normal(0,sigma);
  Vote ~ multinomial(softmax(alpha));
}
generated quantities{
  int pred_boot[Member];
  pred_boot = multinomial_rng(softmax(alpha),sum(Vote));
}

関連エントリ

abrahamcow.hatenablog.com