この記事はStan Advent Calendar 2019の15日目のエントリー記事です。この記事では，「あなたは何%の試行で方略Aを用いましたか？」といった質問項目を実施したときに得られるパーセントデータの分析事例を紹介します。
1. 分析に用いるデータ
今回分析するデータのヒストグラムを図1に示します。
このデータは，筆者が以前行った研究 (Muto, Matsushita, & Morikawa, 2019)で収集したデータの一部です (OSFで公開しています)。詳細は省きますが，この実験では60名の実験参加者 (12名×実験5つ) に2種類の認知課題を遂行してもらったあとで，それぞれの課題において，研究者があらかじめが想定していた方略Aと方略Bについて実験参加者が使用した割合を0%〜100%の範囲で評定してもらいました。測定には図2のような視覚アナログスケール (visual analogue scale) を使用しました。参加者は，0〜100%の線分上の当てはまる位置に縦線を引くことで割合を回答します。データは0〜100の範囲の整数値として記録されました。
60名の参加者が課題1・課題2のそれぞれについて方略A・方略Bの使用頻度を評定したため方略Aのパーセントデータは60名×2=120個の要素から成る反復測定データですが，今回は簡単のために対応のないデータとして扱います (つまり，120名の参加者が1回ずつ評定した場合と同じとみなす)。図1のヒストグラムはこれらの120個の値のヒストグラムです。 (ちなみに，ここでいう方略Aとは論文中でegocentric transformation strategyと呼んでいる方略ですが，この記事では単に方略Aと呼ぶことにします。方略Bことmental scanning strategyに関してはこの記事では扱いませんが，方略Bで同じ分析をしても類似した結果が得られると予想されます。)
図1のヒストグラムの両端 (濃い棒) は評定値がちょうど0%または100%であった参加者の割合を表し，それ以外の薄い棒は[1%, 99%]の範囲のデータについてbinの幅を9としてそれぞれの割合をプロットしたものです。右端の100%のバーがその隣の[91%, 99%]のバーよりも低いのはbinの幅が異なるためであって，100%と回答した参加者が99%と回答した参加者よりも少なかったことを意味する訳ではないので注意。このヒストグラムを見ると，方略Aが用いられた割合は0%付近と100%付近にピークを持ち，かつ0%と100%の値を含んでいることが分かります。つまり，方略Aを使う参加者はほとんどの試行で方略Aを用いるが，方略Aを使わない参加者は一貫して使わない，といった傾向がありそうです (※実際には，参加者の違いではなく実験状況の違いによって方略の使用頻度が二極化すると考えられます)。このようなデータの生成メカニズムを考えるときに真っ先に思い浮かぶのがベータ分布でしょう (図3)。ベータ分布は(0,1)の範囲の実数を実現値として返す確率分布で，非負のパラメータ α と β によって分布の形が様々に変化します。特に 0 < α < 1 かつ 0 < β < 1 のときにはU字型の分布となり，α = β であれば左右対称の分布，α < β であれば0%側に偏った分布，α > β であれば100%側に偏った分布になります。
この分布を使えば今回のデータをうまく説明できそうな気がしますが，残念ながら両端の0と1に対しては確率密度を定義できない (無限大に発散してしまう) ため，今回のように0%と100%を含むデータに対してはそのまま当てはめることができません。代替案としてゼロワン過剰ベータ分布を使えば恐らくフィッティングはうまくいきますが解釈が困難になります (単純にパラメータが多くなり，しかも評定値が0%になるメカニズムと1%になるメカニズムが質的に異なることを仮定してしまうので，今回のデータには恐らく不適切)。今回のデータをうまく説明するためには一工夫必要なようです。このようなデータの背後にある生成過程を考えてモデルを構築してみましょう。
2. モデル
今回用いる評定値データは，その前に行った320試行から成る認知課題において方略Aを用いた割合の評定値です。そこでまず，ある実験参加者が認知課題において方略Aを用いる確率を θ と置くと，実際に方略Aを用いた試行の数 ν は，二項分布を用いて
と表現できます。ν も θ も直接観測することができない未知のパラメータです。さて，もし実験参加者が，自分が方略Aを用いた試行の数 (= ν) と総試行数 (= 320) を正確に知っていれば，実験参加者の評定値は ν / 320 になるはずです。ところが現実的には記憶の曖昧さや評定値の記入の不正確さなどに由来する誤差が生じるはずなので，実際に観測される評定値 Y は
という正規分布に従うと考えることにします。このモデルには離散パラメータの ν が含まれるため，今回はStanではなくJAGSを用いてパラメータの推定を試みます。が，回答誤差のパラメータ σ がうまく収束せず，試行錯誤もむなしく本記事の執筆までに解決することができなかったので，今回は σ を推定せず，σ = 0.0001と決め打ちして α と β を推定することにします (回答の誤差がほとんどないという強い仮定を置くことを意味します)。うまく収束しなかった原因として考えられるのは，(1) σ だけでなく方略Aが用いられる確率 θ にも測定誤差に関する情報が反映されてしまって識別できなかった，(2) 正規分布が誤差分布として適切ではなかった (例えば両端では誤差は小さく50%付近で誤差が最大となるといった仮定を含めたり，[0,1]の範囲に収まるように工夫したりする必要がある？)，(3) 筆者がJAGSの仕様をきちんと理解していないためモデルの書き方が悪かった，などがあり得ると思います (原因が分かる方がおられましたら教えて頂けると嬉しいです)。とりあえずこの点には目をつむり，以下のJAGSコードを書いて分析してみました。ちなみにJAGSはStanと比べると収束効率が悪いので，バーンイン期間を (Stanで実装するときよりも) 若干長めにしたり，初期値を明示したりといった工夫をしたほうがよい場合があります。今回はn.chains=4, n.iter=5000, n.burnin=1000, n.thin = 2とし，αの初期値を(0.1,0.3)の一様分布，βの初期値を(0.2,0.4)の一様分布から発生させました。
jagsmodel.txt
model{
#prior
alpha ~ dunif(0.0001,2)
beta ~ dunif(0.0001,2)
sigma <- 0.0001
tau <- 1/pow(sigma,2)
#model
for (i in 1:N){
theta[i] ~ dbeta(alpha, beta)T(0.00001,0.99999)
nu[i] ~ dbinom(theta[i], 320)
Y[i] ~ dnorm(nu[i] / 320, tau)
}
#pred
theta_pred ~ dbeta(alpha, beta)
nu_pred ~ dbinom(theta_pred, 320)
Y_pred ~ dnorm(nu_pred / 320, tau)
}
3. 分析結果
図4は，今回のデータのヒストグラムに事後予測分布を重ねたプロットです。予測分布の[73%,90%]のあたりが実測値よりも若干低めですが，全体的なU字のパターンはうまく表現できているようです。パラメータのEAP推定値はα = 0.15 (95% CI = [0.08, 0.22])，β = 0.26 (95% CI = [0.18, 0.36])でした。α < β なので，100%付近の参加者よりも0%付近の参加者のほうが人数が多いと解釈できます。
4. 発展編：実験条件間の比較
これだけで終わってしまうのは物足りないので，もうちょっと踏み込んだ分析もしてみました。今回の実験では，理論的に方略Aが優勢になると予想される実験条件Aと，方略Aではなく方略Bが優勢になると予想される実験条件Bの2通りを設定していました。そこで，実際にこの予想が当たっているかどうかを確かめるために，実験条件Aと実験条件Bのそれぞれについて別々に α と β を推定してみました。コードは以下の通り。n.chains=4, n.iter=5000, n.burnin=1000, n.thin = 2とし，初期値は特に明示しませんでした。
jagsmodel_2cond.txt
model{
#prior
alpha1 ~ dunif(0.0001,2)
beta1 ~ dunif(0.0001,2)
alpha2 ~ dunif(0.0001,2)
beta2 ~ dunif(0.0001,2)
sigma <- 0.0001
tau <- 1/pow(sigma,2)
#model
for (i in 1:N1){
theta1[i] ~ dbeta(alpha1, beta1)T(0.00001,0.99999)
nu1[i] ~ dbinom(theta1[i], 320)
Y1[i] ~ dnorm(nu1[i] / 320, tau)
}
for (i in 1:N2){
theta2[i] ~ dbeta(alpha2, beta2)T(0.00001,0.99999)
nu2[i] ~ dbinom(theta2[i], 320)
Y2[i] ~ dnorm(nu2[i] / 320, tau)
}
#pred
theta_pred1 ~ dbeta(alpha1, beta1)
nu_pred1 ~ dbinom(theta_pred1, 320)
Y_pred1 ~ dnorm(nu_pred1 / 320, tau)
theta_pred2 ~ dbeta(alpha2, beta2)
nu_pred2 ~ dbinom(theta_pred2, 320)
Y_pred2 ~ dnorm(nu_pred2 / 320, tau)
#generated quantities
phi1 <- log(alpha1 / beta1)
phi2 <- log(alpha2 / beta2)
}
実験条件AとBのそれぞれにおける評定値のヒストグラムと予測分布を以下の図5に示しました。ヒストグラムから，条件Aでは方略Aの使用頻度のピークが100%付近にあるのに対し，条件Bでは逆になっていることが分かります。事後予測分布もその特徴をうまく捉えているように見えます。EAP推定値は，条件Aでは α = 1.18 (95% CI = [0.70, 1.77])，β = 0.36 (95% CI = [0.21, 0.52])，条件Bでは α = 0.06 (95% CI = [0.004, 0.14])，β = 0.49 (95% CI = [0.30, 0.75])でした。ここで，α と β の大小関係から分布の偏り度合いが分かるということを利用して，φ = ln(α/β) という生成量を作ってみます。この生成量 φ は，α = β のとき，すなわち方略Aの使用頻度が0%付近と100%付近の両方に等しくピークを持ち分布が左右対称となるときに0と一致します。また，φ が0より小さければ小さいほど0%側に，0より大きければ大きいほど100%側に偏った分布を表します。EAP推定値は条件Aで φ = 1.18 (95% CI = [0.77,1.61])，条件Bで φ = -2.32 (95% CI = [-4.50, -1.33]) でした。これらの結果から，先ほどの予想通り，条件Aでは多くの参加者が方略Aを一貫して用いるのに対し，条件Bでは多くの参加者が方略Aをほとんど使わないということを定量的に示すことができました。
5. まとめ
という訳で，この記事では0%と100%を含むパーセントデータのモデリングを試みました。方略Aが課題中に用いられる確率がベータ分布に従うと仮定し，実際に方略Aが何試行で使われたかのメタ認知によって評定値が決まるという発想です。フィッティングするだけならゼロワン過剰ベータ分布や切断正規分布の混合分布などを使うやり方もありますし，例えば0%と100%のデータをそれぞれ1%と99%に置き換えるといった小細工をすれば単なるベータ分布モデルでも似たような結果を得ることはできるのですが，ドメイン知識を利用することで倹約的かつ解釈しやすいモデルを無理なく作ることができうることをこの記事でお示しできたのではないかと思います。σ の推定がうまくいかないという問題もありますし，まだまだ改善の余地はあるはずなので，もう少し吟味してみようと思います。
以前，某研究会で講演者の方がご自身の論文で報告したデータを再分析してモデリングされていたのを見て，楽しそうだなーと思ったのがこの記事を執筆するきっかけになりました。元の論文 (Muto et al., 2019) では，今回の評定値データはメインの行動データの解釈に役立てるための補助的なデータという位置づけに過ぎなかったので，条件間の平均値差の検定だけしてお茶を濁していました。それはそれで問題ないとは思うのですが，もっと上手な分析方法があるのではないかとずっともやもやしていたというのもあって，この機会に再分析＆モデリングを試みたという訳です。この論文のデータはオープンにしているので実は誰でもモデリングをすることが可能です。オープンデータの営みがもっと普及すればモデリングの可能性も広がる気がしてワクワクしますね。
そしてこの記事を書き終えてから気付いたのですが，この記事は「Stan」 Advent Calendarの記事のはずなのにStan全く使ってねえ……。StanでできないことをJAGSにやってもらったということで許してクレメンス。
引用文献：
Muto, H., Matsushita, S., & Morikawa, K. (2019). Object's symmetry alters spatial perspective-taking processes. Cognition, 191, 103987. https://doi.org/10.1016/j.cognition.2019.05.024
