ベイズマックスへようこそ。
このブログでは,心理学と統計が大好きなmutopsyが,
統計に関するあれやこれやを自由気ままに投稿します。

2017年12月14日

実験心理学者のためのベイジアンモデリング──正答率データをどうやって分析するか──

 こんにちは。mutopsyです。この記事は,Stan Advent Calendar 2017の14日目のエントリー記事です。この記事では,実験心理学でよく扱われる正答率・誤答率(あるいは正答数・誤答数)データの分析を例に,よくある素朴な分析(平均値の点推定・区間推定およびt検定)の背後にあるモデルを可視化した上で,そのモデルを改善した別の3つのモデルを紹介します。また,情報量基準によるモデル比較も行います。そして最後に,実験心理学におけるベイジアンモデリングの有用性について述べたいと思います。

1. はじめに:正答率データを素朴に分析してみる

 改めまして,こんにちは。mutopsyこと武藤拓之です。大阪大学で認知心理学の研究をしています。この記事では,実験心理学でもそれ以外の分野でもお馴染みの,正答率(正答数)データを分析する方法について考えてみたいと思います。昨日はフォン・ミーゼス分布に関する若干マニアックな記事を投稿しましたが,この記事ではもう少しだけ需要のありそうな話をするつもりです。

 さて,いきなりですが,下のグラフ(図1)をご覧ください。このグラフは,筆者が過去に行った記憶に関する実験のデータを要約したものです1。この実験では,実験参加者さんに2枚の画像を順番に見てもらい,2枚目の画像に写っている対象が1枚目の画像と同じであったか異なっていたかを回答してもらいました。この試行を繰り返すと,実験参加者さん一人ひとりの正答数のデータを得ることができます。正答数を合計試行数で割れば正答率になるので,正答率のデータが実験参加者の人数分だけ手に入ると考えることもできます。図1のグラフの青い点は,実験参加者一人ひとりの正答率をプロットしたものです。この実験では難易度の異なる2つの条件が設定されており,条件1と条件2にはそれぞれ異なる実験参加者が割り当てられたと考えてください。なお,この記事の内容は,正答数と合計試行数が分かっているデータであればどんなデータにも当てはまる話なので,読者の皆様にとって分かりやすいデータをイメージして頂ければ良いかと思います。

ac_rawdata.png
図1. この記事で分析に用いる正答率データの平均値とその95%信頼区間。青い点は個人ごとの正答率。

 このデータのcsvファイルはこちらからダウンロードできます。→ accuracy_data.csv

 このデータは,各行が実験参加者1人分のデータに対応しており,以下の5つの変数を含んでいます。

  • ID: 実験参加者のID(1から24までの整数)
  • N.Corr: 正解した試行数(1から6までの整数)
  • N.Total: 合計試行数(6で固定)
  • Corr.Rate: 正答率(「正解した試行数」を「合計試行数」で割った値)
  • Group: 条件群(1または2)

 実験心理学の文脈では,このようなデータが得られたら,まずは平均値の点推定値と誤差範囲を計算するのが一般的だと思います。ここでは各条件における平均正答率に関心があるので,正答率の算術平均(=母平均の点推定値)と母平均の95%信頼区間を計算してみました。これらの値は先ほどの図1に折れ線グラフとエラーバーで示されています。赤色の四角形が各条件の平均正答率,エラーバーがその95%信頼区間です。信頼区間ではなく標準誤差をプロットする場合も多いですが,このようにして算術平均とエラーバーをグラフに描くのが一般的かと思います。データの傾向を確認したら,例えば2つの条件間で平均正答率に差があるか否かを検定したり(Welchのt検定),あるいはある条件の平均正答率がチャンスレベル(ここでは50%)よりも有意に高いと言えるかどうかを検定したり(一標本のt検定),研究仮説に合った分析を行います。実際にこのデータに対してWelchのt検定を行うと,p = .024 という結果が得られ,条件2の平均正答率は条件1の平均正答率よりも平均して23.6%高い,という結論が得られます。

 ただ,改めて図1のグラフを見ると,気になる点がいくつかあります。まず,条件2の平均正答率の95%信頼区間をよく見ると,区間の上側が1.0 (=100%) の値よりも上にあります。正答率は理論上0%から100%の範囲内に収まるはずなのに,平均正答率の母数が100%を超過しているかもしれないという奇妙な推定結果が得られてしまいました。また,実験参加者一人ひとりの正答率の散らばり方も気になります。例えば,条件1でも条件2でも正答率が100%の実験参加者がそれなりにいるように見えます。少なくとも正規分布に従っているようには見えません。たとえデータそのものが正規分布に従っていなくても,サンプルサイズが十分に大きければ平均値の分布は中心極限定理により正規分布に近付きますが,今回のように高々 n = 12 の場合には微妙な気もします。それに,やはり平均値の信頼区間が100%をまたいでしまうのは何となく気持ち悪いです。

 機械的に信頼区間を求めたりt検定を行ったり,あるいは分散分析したり回帰分析したり,ということを行ってしまうと,これらの分析が結局何をやっているのかが見えにくくなってしまいます。そこでまず,上述した分析の背後にあるモデルを可視化してみようと思います。そうすれば,先の分析の違和感の正体が掴めるはずです。

2. モデル1:制約のない正規分布モデル

 この節では,1節で説明したような素朴な分析が何をやっているのかを可視化するために,先の分析の背後にあるモデルをベイズ的に表現してみようと思います。具体的には平均値の区間推定の背後にあるモデルをStanで記述してみます(平均値差の推定については6節で触れます)。「頻度主義 vs. ベイジアン」という枠組みで考えてしまうと違和感を持たれるかもしれませんが,以下のモデルはいわゆる頻度主義の枠組みで想定されているモデルと(ほぼ)等価なモデルであり,後で確認するように,推定結果も(ほぼ)一致します。なお,この記事ではパラメータの解釈に関する頻度主義とベイズの違いについては触れません。前節で説明したようなよく使われる分析方法の背後にあるモデルは,以下に示す「モデル1:制約のない正規分布モデル」として記述することができます2。特に重要な行を赤色でハイライトしています。

Stanコード01: 【モデル1】制約のない正規分布モデル
data {
  int N;  //データの総数
  real CorrRate[N];  //正答率
  int Group[N];  //実験条件(群)
}

parameters {
  real mu[2];  //平均正答率 (μ)
  real<lower=0> sigma[2];  //標準偏差 (σ)
}

model {
  //モデル
  for (n in 1:N){
    CorrRate[n] ~ normal(mu[Group[n]],sigma[Group[n]]);
    //↑正答率データが正規分布に従う
  }
  
  //事前分布
  for (k in 1:2){
    mu[k] ~ uniform(-9999999,9999999); 
    sigma[k] ~ uniform(0,9999999);
  }
}

generated quantities{
  real diff_mu;
  diff_mu = mu[2] - mu[1];  //2群の平均正答率の差
}

 このモデルでは,条件1と条件2の正答率データがそれぞれ独立な正規分布から生成されると仮定しています。正規分布は平均値μと標準偏差σという2つのパラメータを持ちますが,今回の分析ではこの平均値μの分布が求めたい平均正答率の分布に相当します。以下,コードの解説です。

  • 1-5行目: dataブロックでは分析に用いるデータを宣言します。ここではデータの総数(今回はデータが24行なのでN=24),正答率,実験条件(群)を使用します。
  • 7-10行目: parametersプロックでは,このモデルで推定するパラメータを宣言します。ここでは正規分布のパラメータμとσを条件ごとに推定します。σは非負の値なので下限を0に指定しておきます(これを指定しないとうまく推定できません)。
  • 15行目: 正規分布から一人ひとりの正答率データが生成されることを表現しています。この部分がこのモデルの核となります。
  • 21-22行目: 事前分布を指定しています。この部分を省略した場合,Stanは十分に幅の広い一様分布を事前分布として採用します。なので本当はこの2行は不要なのですが,モデルの仮定を示すためにあえて明示的に書いています。伝統的な分析による推定結果は,-∞から+∞の範囲の一様分布3を事前分布としてベイズ推定した時の結果と一致します。ここで指定しているuniform(-9999999,9999999)は,-9999999から+9999999までの範囲の一様分布を表します。つまり,平均値がいくつになるか全く分からないという事前の信念を表現しています。本当はuniform(-Inf,Inf)とでも書きたいところですが,Stanではこのような書き方は(おそらく)できないので,あえて極端に範囲の広い事前分布を指定してみました。
  • 26-29行目: 条件間の平均値差を計算しています。単にmu[0]とmu[1]の引き算をしているだけです。

 このモデルに先ほどのデータを当てはめて,推定したμの事後モード(MAP推定値)と95%最高密度区間を求めると,以下の図2のような結果が得られます。個人のデータのプロットはローデータそのものなので,値は図1と厳密に同じです(描画の都合上,横位置は若干ずれています)。折れ線グラフで表された平均値(μ)の点推定値(MAP推定値)とその95%最高密度区間は上記のコードで推定したものです。推定上の誤差が多少あるとはいえ,図1とほとんど同じ結果になっていることが分かるかと思います。

ac_model1.png
図2. 制約のない正規分布モデル(モデル1)で推定した平均正答率。エラーバーは95%最高密度区間を表す。

 改めてこのモデルを見てみると,まず事前分布の箇所が気になるかと思います。正答率の平均値は理論上0から1の範囲に収まるはずなのに,事前知識としてほぼ-∞から+∞の範囲の一様分布を用いるのは明らかに不自然です。素朴な信頼区間の計算やt検定の背後にはこのようなモデルが想定されているのです。そこで次のモデルでは,「平均正答率(μ)は0から1までの値しか取らない」という事前知識をパラメータの事前分布として反映させてみようと思います。

3. モデル2:μの範囲を制限した正規分布モデル

 モデル2は,正規分布の平均値パラメータμについて「0から1までの値しか取らない」という事前知識を与えたモデルとなります。先ほどのモデル1を僅かに修正するだけでこのモデルを作ることができます。以下のコードの赤色でハイライトされた部分のみがモデル1と異なる部分です。

Stanコード02: 【モデル2】μの範囲を制限した正規分布モデル
data {
  int N;  //データの総数
  real CorrRate[N];  //正答率
  int Group[N];  //実験条件(群)
}

parameters {
  real<lower=0, upper=1> mu[2];  //平均正答率 (μ)
  real<lower=0> sigma[2];  //標準偏差 (σ)
}

model {
  //モデル
  for (n in 1:N){
    CorrRate[n] ~ normal(mu[Group[n]],sigma[Group[n]]);
    //↑正答率データが正規分布に従う
  }
  
  //事前分布
  for (k in 1:2){
    mu[k] ~ uniform(0,1); 
    sigma[k] ~ uniform(0,9999999);
  }
}

generated quantities{
  real diff_mu;
  diff_mu = mu[2] - mu[1];  //2群の平均正答率の差
  }
}

  • 8行目: parametersブロックで平均値パラメータを宣言する時に,あらかじめ下限0と上限1と指定しておきます。
  • 21行目: 平均値パラメータの事前分布として,0から1の範囲の一様分布を指定しています。parametersブロックで平均値パラメータの下限と上限を指定しているので,この行を省略しても自動的に範囲[0,1]の一様分布が事前分布として指定されます。しかし,この上限と下限の制約がパラメータの事前分布であるということを強調するために,あえてここでは明示的に書いています。

 このモデルによる推定結果は以下の図3のようになります。目論見通り,条件2の区間推定幅の上端が100%を越えなくなりました。なお,この記事の下部にある図6にこの記事で用いた全てのモデルの結果を示しているので,比較の際にはこちらをご覧いただくと分かりやすいかと思います。

ac_model2.png
図3. μの範囲を限定した正規分布モデル(モデル2)で推定した平均正答率。エラーバーは95%最高密度区間を表す。

 さて,平均値の範囲の問題は一応解決しましたが,そもそも正規分布を用いたモデルで本当に良いのだろうかという疑問は残ります。次の節では,正規分布ではなく二項分布を用いたモデルについて説明します。

4. モデル3:θを固定した二項分布モデル

 表が出る確率がθであるコインをn回投げた時に表が出る回数kは,pとnをパラメータに持つ二項分布で表現できます。もちろんコイン投げ以外のケースにも二項分布は使えます。例えば,記憶実験のある試行で正答する確率をθ,合計試行数をnと置けば,正答する回数kはやはり二項分布に従います(二項分布を使ったモベイズ推定については過去の記事でも紹介しています)。二項分布はnが大きくなるにつれて正規分布に近づきますが,今回分析に使用しているデータの合計試行数はたったの6なので,正規分布ではなく二項分布を使ったモデルのほうが適切であるように思われます(正規分布を用いるということは,二項分布モデルにおいてnが∞であると仮定することと一致します)。また,kとnは正の整数で,k > nとなることはなく,確率θは0から1の範囲に収まります。つまり,二項分布を用いる以上はモデル2のようにわざわざ範囲を制限する必要がありません(モデル上では制限をかけないと収束しませんが,この制限は事前の信念に基づく制限ではなく,モデルから自然に要求される制限であるという点で異なります)。

 それでは二項分布を使ったモデルについて考えてみましょう。上述したように,ある試行で正答する確率θ,合計試行数nの二項分布から正答回数が生じるモデルを想定しますが,まずは,θが条件のみによって決定するというモデルを考えます。つまり,条件1における正答確率をθ1と置くと,条件1の実験参加者は,どの参加者であっても任意の試行においてθ1の確率で正答し,1-θ1の確率で誤答すると考えます。このパラメータθは平均正答率ではなく各試行で共通の正答確率を表すという点でモデル1とモデル2のμとは厳密には異なります。ただし,このモデルでは同じ条件の参加者でθが共通しているという仮定を置いているため,各条件における正答確率θの参加者間平均は正答確率θと一致します(θ * 12 / 12 = θ)。そのため,ここではθをμと同じように解釈することにします(研究仮説において両者を区別する必要が特になければ問題ありません)。このモデルのStanコードは以下のようになります。推定するパラメータがmuからthetaに変わっていることと,データが正規分布ではなく二項分布に従うと仮定していることが先のモデルとの違いです。

Stanコード03: 【モデル3】θを固定した二項分布モデル
data {
  int N;  //データの総数
  int CorrResp[N];  //正答数
  int TotalResp[N]; //合計試行数
  int Group[N]; //実験条件(群)
}

parameters {
  real<lower=0, upper=1> theta[2]; //正答確率 (θ)
}

model {
  //モデル
  for (n in 1:N){
    CorrResp[n] ~ binomial(TotalResp[n],theta[Group[n]]);
    //↑正答数データが二項分布に従う
  }
  
  //事前分布
  for (k in 1:2){
    theta[k] ~ uniform(0,1);
  }
}

generated quantities{
  real diff_theta;
  diff_theta = theta[2] - theta[1];  //2群の正答確率の差
}

  • 3-4行目: 今度は正答率ではなく正答数と合計試行数をデータとして与えます。今回のデータでは合計試行数は6で固定ですが,この値が実験参加者によって違っていてもこのモデルは走ります。
  • 9行目: 正規分布の時とは推定するパラメータが違うので,名前もmuではなくthetaにしています。thetaの範囲は0から1に指定しておかないとうまく収束しません。
  • 15行目: 二項分布から正答数データが生成されるモデルを表現しています。

 このモデルを当てはめると以下の図4のような結果が得られます(他のモデルとの比較は図6参照)。縦軸が平均正答率(μ)ではなく正答確率(θ)になっていることに注意してください。パラメータのMAP推定値はモデル1・モデル2とあまり変わっていませんが,これまでのモデルよりも条件1における区間推定幅がかなり狭くなっていることが読み取れます。また,当然ながら区間推定幅は100%をまたぎません。

ac_model3.png
図4. θを固定した二項分布モデル(モデル3)で推定した正答確率。エラーバーは95%最高密度区間を表す。

 この節では二項分布を使ったモデルを紹介しましたが,「θが実験参加者にかかわらず一定である」という仮定はなんとなく直感に反するように思えます。記憶成績が良い人も悪い人もいる,というように,θにも個人差があると考える方が自然でしょう。そこで,次の節ではθが人によって異なるモデルの一例を示します。

5. モデル4:θが個人ごとに異なる二項分布モデル

 この節では,二項分布のパラメータθが人によって異なることを仮定したモデルの一例を紹介します。以下にそのStanコードを示します。赤色でハイライトされた行はモデル3と異なる箇所です。

Stanコード04: 【モデル4】θが個人ごとに異なる二項分布モデル
data {
  int N;  //データの総数
  int CorrResp[N];  //正答数
  int TotalResp[N]; //合計試行数
  int Group[N]; //実験条件(群)
}

parameters {
  real<lower=0, upper=1> theta[N]; //正答確率 (θ)
}

model {
  //モデル
  for (n in 1:N){
    CorrResp[n] ~ binomial(TotalResp[n],theta[n]);
    //↑正答数データが二項分布に従う
  }
  
  //事前分布
  for (n in 1:N){
    theta[n] ~ uniform(0,1);
  }
}

generated quantities{
  real diff_mu;
  real mu[2];
  mu[1] = 0;
  mu[2] = 0;
  for (n in 1:N) mu[Group[n]] = mu[Group[n]] + theta[n]/12;
  diff_mu = mu[2] - mu[1];
}

  • 9行目: 前節のモデル3ではパラメータthetaは条件によって異なると仮定していため,2つのthetaを推定すれば良かったのですが,今回のモデルでは実験参加者の人数分だけthetaを推定します。
  • 15行目: 条件ごとのthetaではなく実験参加者ごとのthetaを二項分布のパラメータに指定します。
  • 20-21行目: 事前分布も全てのthetaに対して指定します。
  • 27-30行目: 推定したthetaの平均値muを条件ごとに計算しています。このmuは正規分布のパラメータμではなく,ただの算術平均です(thetaが分布するのでmuも分布します)。

 このモデルでパラメータを推定した結果を以下の図5に示します(モデルの比較は図6を参照)。縦軸は正答確率 (θ) の参加者間平均です。モデル3に比べると,条件2のパラメータの推定値が低くなっています。

ac_model4.png
図5. θが個人ごとに異なる二項分布モデル(モデル4)で推定した正答確率の平均。エラーバーは95%最高密度区間を表す。

6. 4つのモデルを比較する

 さて,これまで4つのモデルを使って正答率データのモデリングを行ってきましたが,結局どのモデルを使うのが良いのでしょうか。図6は4つのモデルによる推定結果を並べたものですが,モデルによって結果が異なっている様子が見て取れます。

modelcomp.png
図6. 4つのモデルによる正答率(正答確率)の推定結果の比較。

 ここではモデル比較のために,情報量基準のWAICを用いることにします。WAICは予測の良さを表す指標で,値が小さいモデルほど真のモデルに近いことを意味します4。それぞれのモデルのWAICは図6に示した通りです5。この結果からは,モデルの良さは「モデル1 < モデル2 < モデル3 < モデル4」であると言えそうです。やはり,素朴な分析の背後にあるモデルと等価な正規分布モデルはあまり適切ではないようです。この中ではθが人によって異なる二項分布モデルが最も良いモデルであるといえそうです。

 これまでは触れてきませんでしたが,各モデルにおける条件間の(平均)正答率の差(diff_muまたはdiff_theta)の推定結果も見てみましょう。以下の図7に結果を示しました。これを見ると,どのモデルにおいても(平均)正答率の差の区間推定幅には0が含まれていませんが,幅の広さにはだいぶ違いが見られます。モデル1の95%最高密度区間はWelchのt検定の計算過程で用いられる95%信頼区間と一致します。どうやら二項分布モデルのほうが推定精度が良さそうです。

ac_diff.png
図7. 4つのモデルによる正答率の条件間差の推定結果。

 ただし,可能なモデルは他にもあるので(例えばベータ二項分布モデルのような,θに分布を仮定するモデルや,感度と反応バイアスを区別する信号検出理論モデル,項目困難度と能力を区別する項目反応理論モデルなど),今回のモデル4があらゆるモデルの中で最良であるとは限りません。研究仮説に応じて適切なモデルを選択することも必要でしょう。

7. まとめ

 この記事では,素朴な分析(平均値の点推定・区間推定およびt検定)の背後にあるモデルを可視化した上で,そのモデルを修正した別の3つのモデルを紹介しました。また,モデルごとの推定結果の違いを概観し,WAICを用いたモデル比較も行いました。

 この記事で最も主張したかったことは,分析の背後にあるモデルを意識することの重要性についてです。機械的に分析を実行することに慣れてしまうと,明らかにデータに合わないモデルを無意識に用いてしまう危険性があります。Stanを使ってベイズ的にモデルを記述すれば,そのモデルが何を仮定しているのかは簡単に分かりますし,この記事で実際にお見せしたように,モデルの修正も容易に行えます。この記事で紹介したモデル1-4のパラメータは,ベイズを用いないで最尤法で解くことも可能だとは思いますが,ベイズ推定で用いられるMCMCアルゴリズムは強力かつ汎用的で,またStanを使えばモデルの記述も修正も拡張も比較的簡単に行えるという点で,非常に一般性が高いです。共通の文法を用いてあらゆるモデルを記述しパラメータを効率的に推定できるベイジアンモデリングはとても魅力的であるように思います。ベイズというと,パラメータの解釈における頻度主義との違いが強調されることが多いですが,個人的にはモデリングの柔軟性こそがベイズの良さであるように感じています。あなたもベイジアンモデラーになってみませんか。

脚注

1本当はもっと多くの条件を設定していたのですが,ここでは2つの条件のデータのみを抽出しています。また,実際にはこの2条件は対応のあるデータなのですが,分かりやすさのために,ここでは対応のないデータとして扱うことにします。とはいえこの記事で書かれている内容は,もっと条件数の多い実験で得られたデータや対応のあるデータに対してももちろん適用可能です(Stanコードは適宜修正しなければなりませんが)。

2分かりやすさを優先し,各実験参加者の正答率が平均μ(=平均正答確率)の正規分布に従うモデル,すなわち平均正答確率を直接推定するモデルを書きました。しかし,正答率ではなく正答数が正規分布に従うモデルを書き,generated quantitiesブロックで平均正答数を合計試行数で割って平均正答率を計算するほうが,モデルの適用範囲は広がります。また,渡すデータが後の二項分布モデル(モデル3・モデル4)と同じになるのでそのままWAICを計算できるようになります。

3厳密に言えば,この関数は全区間で積分しても1にならないので確率密度関数ではない。

4情報量基準は万能ではなく,あくまでも1つの判断基準でしかないのでご注意ください。

5モデル比較のためにWAICを計算する際には,どのモデルにも同じデータを与える必要があります。上述のモデル1とモデル2では正答率データを渡していたのに対し,モデル3とモデル4では正答数と合計試行数のデータを渡していたので,そのままWAICを計算すると,モデル1とモデル2のWAICが極端に小さくなってしまいます。そのため,モデル1とモデル2のWAICを計算する際は,正答確率ではなく正答数が正規分布に従うモデルに修正しなければなりません。この時,修正したモデルには合計試行数のデータを渡す必要はありません(正規分布モデルは合計試行数が∞の二項分布モデルと等しいため)。このことは関西学院大学の清水先生に教えて頂きました!

2018年1月25日:csvデータへのリンク切れを修正しました。

posted by mutopsy at 01:45 | 統計