hBayesDMでお手軽にベイジアン認知モデリングをやってみよう！
この記事は，Stan Advent Calendar 2017の16日目の記事です。
Stan Advent Calendar 2017をみると，Stanを覚えると色々なことができるなあと感想をもたれるのではないでしょうか？しかし，同時に「Stan難しそう」，「長くコードを書く自信がない」，「最初の一歩を踏み出すのが難しい」などの感想もあるかなと思います。そこで，この記事では，「Stanとかとりあえずいいから，気軽にベイジアン認知モデリングをやってみよう！」ということを目指します。
心理学や神経科学では，なんらかの認知課題を参加者にしてもらって，そのパフォーマンス（行動データ）を評価することが多いです。行動データの評価は重要ですが，その行動データが生成された背景にはなんらかのプロセスが存在する可能性があります。認知モデリングは，行動データが生成される過程を数理モデルで表現し，行動データから数理モデルのパラメータを推定する方法です。そして，ベイジアン認知モデリングでは，ベイズ統計学を用いて認知モデリングを行います。今回の記事では，ベイジアン認知モデリングを簡単にやってくれるhBayesDMパッケージを紹介します。
hBayesDMとは？
hBayesDMは，Woo-Young Ahn博士1の研究チームが開発したRパッケージです。hBayesDMは，hierarchical Bayesian modeling of Decision-Making tasksを略したものです。つまり，意思決定課題などの心理学・神経科学で用いられる認知課題を階層ベイズモデリングするパッケージです。認知課題の階層ベイズモデリングは，決して簡単とは言い難いStanコードを書く必要があるのですが，hBayesDMはとてもお手軽に階層ベイズモデリングをしてくれます（1行のコードでベイジアン認知モデリングができちゃいます！）。ここでは，階層ベイズについての説明はしません。日本社会心理学会 第4回春の方法論セミナーで簡単に説明していますので，参照ください。また，hBayesDMの詳細は，Ahn博士の作成した“hBayesDM Getting Started”を参照ください。
hBayesDMで扱うモデル
hBayesDM v0.4.0では，以下のモデルが使えます（今後も増える予定のようです）。たくさんのモデルが準備されていますね！
| 実験課題 | モデル | 関数名 |
|---|---|---|
| 選択反応時間課題 | Drift diffusion model(Ratcliff, 1978) | choiceRT_ddm |
| Linear Ballistic Accumulator model(Brown & Heathcote,2008) | choiceRT_lba | |
| 遅延割引課題 | Constant-Sensitivity (CS) model(Ebert & Prelec,2007) | dd_cs |
| Exponential model(Samuelson, 1937) | dd_exp | |
| Hyperbolic model(Mazur,1987) | dd_hyp | |
| アイオワギャンブル課題 | Prospect Valence Learning-DecayRI (Ahn et al.,2011;2014) | igt_pvl_decay |
| Prospect Valence Learning-Delta (Ahn et al., 2008) | igt_pvl_delta | |
| Value-Plus-Perseverance (Worthy et al., 2013) | igt_vpp | |
| 直交GO/NoGo課題 | RW+noise(Guitart-Masip et al., 2012) | gng_m1 |
| RW+noise+go bias(Guitart-Masip et al., 2012) | gng_m2 | |
| RW+noise+go bias+Pav.bias(Guitart-Masip et al., 2012) | gng_m3 | |
| M5(Cavanagh et al.,2013の表１) | gng_m4 | |
| 確率的逆転学習課題 | Experience-Weighted Attraction(Ouden et al., 2013) ) | prl_ewa |
| Fictitious update(Gläscher et al., 2009) | prl_fictitious | |
| Reward-Punishment(Ouden et al., 2013) | prl_rp | |
| リスク回避課題 | Prospect Theory(Sokol-Hessner et al., 2009) | ra_prospect |
| PT without loss aversion | ra_noLA | |
| PT without risk aversion(Tom et al., 2007) | ra_noRA | |
| ２腕バンディット課題 | Rescorla-Wagner (delta) model(Erev et al., 2010; Hertwig et al., 2004) | bandit2arm_delta |
| ４腕バンディット課題 | Fictive upd.+rew/pun sens.(Seymour et al., 2012) | bandit4arm_4par |
| Fictive upd.+rew/pun sens.+lapse(Seymour et al., 2012) | bandit4arm_lapse | |
| 最後通牒課題 | Ideal Bayesian observer model(Xiang et al., 2013) | ug_bayes |
| Rescorla-Wagner (delta) model(Gu et al., 2015) | ug_delta |
hBayesDMを使う準備
前準備
hBayesDMを使うには，事前に以下を準備しておく必要があります。
また，ggplot2，loo，mail，modeestなどパッケージも必要です（dependencies=TRUEにしておけば，hBayesDMのインストール時に入ります）。
インストール
Windowsユーザーは，以下でインストールしてください。
install.packages("hBayesDM", dependencies=TRUE)MacやLinuxユーザーは，devtoolsパッケージをインストールした上で，以下のようにgithub経由でインストールしてください。
devtools::install_github("CCS-Lab/hBayesDM")hBayesDMで認知モデリングに挑戦！
さくっと使ってみましょう。まずは，hBayesDMをロードします。tidyverseもロードしておきます。tidyverseをインストールされてない方は，install.packages(“tidyverse”, dependencies=TRUE)でインストールください。
library(hBayesDM)
library(tidyverse)遅延割引課題
hBayesDMでは，複数の実験課題に対して複数のモデルが準備されています。今回は，その中から遅延割引についてのモデルを用います。
「今日の1000円」と「1年後の1050円」の２つの選択肢があるとすれば，あなたは，どちらを選ぶでしょうか？
どちらを選ぶかは，人によりますが，「今日の1000円」を選んだ人も多いのではないでしょうか？いつ貰おうと報酬の客観的な金額は変化しないのですが，遅れて貰える報酬は，主観的には価値が割り引かれる傾向があります。このような遅延してもらえる報酬の価値が割り引かれる現象を遅延割引と呼びます。遅延するにつれて価値が割り引かれる様子は，数理モデルによって表現することができます。hBayesDMでは，Constant-Sensitivity model, Exponential model, Hyperbolic modelの３つの数理モデルを準備しています。ここでは，Hyperbolic modelを取り上げて，説明します。
Hyperbolic modelとは？
Mazur(1987)2が提案したHyperbolic modelは，以下のような数理モデルです。Dは遅延になりますので，V(D)は，遅延した報酬の価値を意味します（なお，Vは0から1の範囲の値になります）。V(D)は，1を1+kDで割ったものになりますので，遅延するほど（Dが大きくなるほど），小さな値をとります。
\[ V(D) = \frac{1}{1 + kD} \]
遅延するほど小さな値をとる程度は，割引率（k）によって制御されます。例えば，下図のように，k = 0.001とk = 0.1の場合，kの値が大きいほど遅延による価値割引の程度が大きくなるようになります。
遅延割引課題では，割引率(k)を推定するために，「今日の1000円」と「明日の1100円」のどちらを選ぶかといった選択を何度もしてもらいます。これらの選択肢に対する個人の反応（行動データ）から，kを推定します。なお，kの推定にあたり，選択肢の価値にどのくらい従った反応をするのかにかかわる逆温度（β）パラメータも推定します。このように，認知モデリングでは，時点の異なる選択肢における個人の選択行動の背景にあるプロセスを推定します。
データの準備
さて，遅延割引のHyperbolic modelについての理解が深まったところで，早速推定をしていきましょう！hBayesDMは，複雑なStanコードを書かなくても階層化されたベイジアン認知モデリングをやってくれるのですが，それはある程度，型にはまった形式のデータをセットを使うことで実現されています。
「hBayesDM Getting Started」の「How to use hBayesDM」の「1. Prepare the data」にある「here」というリンクをクリックするとサンプルデータセットを入手できます。今回は，そのサンプルデータセットの中の遅延割引に関するdd_exampleData.txtを用います。では，さっそく，読み込んでみましょう。
(data <- read_tsv("dd_exampleData.txt"))## # A tibble: 2,160 x 7
## subjID trial delay_later amount_later delay_sooner amount_sooner choice
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 6 10.5 0 10 1
## 2 1 2 170 38.3 0 10 1
## 3 1 3 28 13.4 0 10 1
## 4 1 4 28 31.4 0 10 1
## 5 1 5 85 30.9 0 10 1
## 6 1 6 28 21.1 0 10 1
## 7 1 7 28 13 0 10 1
## 8 1 8 1 21.3 0 10 1
## 9 1 9 28 21.1 0 10 1
## 10 1 10 15 30.1 0 10 1
## # … with 2,150 more rowshBayesDMでは，データはテキストファイル(.txt)で保存したものを用います（タブ区切り）。hBayesDMは，簡単にベイジアン認知モデリングができる代わりに，データセットは決められた通りに作成しないといけません。例えば，遅延割引課題の場合，変数名はsubjID, delay_later, amount_later, delay_sooner, amount_sooner, choiceで作成します。このようにデータセットを作らないとパラメータ推定ができません。そういう意味では，hBayesDMを使う前のデータハンドリングが重要になるかもしれない。なお，データは，参加者ごとに試行数が異なっても良いのですが，欠測値(NAも）はデータに含めてはいけません。データハンドリングの段階で，欠測値は除外しておく必要があります。
遅延割引課題のパラメータ推定
遅延割引課題データをhyperbolic modelを用いてパラメータ推定する場合は，dd_hyperbolic()を用います。dd_hyperbolic()の主な引数は以下の通りです。データを適切に作成して，以下の引数を設定するだけで，パラメータ推定できちゃいます。なお，以下では省略していますが，mail引数にメールアドレスを指定しておくと，推定が終了したらメールで通知してくれます。
dd_hyperbolic(data = "データ名.txt",
niter = 反復回数,
nwarmup = ワームアップに指定する反復回数,
nchain = 連鎖の数,
ncore = 並列化の際に使用するコア数,
nthin = 間引き間隔,
inits = 初期値("fixed", "random", 指定),
indPars = MCMCサンプルの要約方法("mean", "median", "mode"),
以下略)では，さっそく，サンプルデータで推定をしてみましょう！以下の１行のコードだけで，hyperbolic modelで階層ベイズモデリングをやってくれます。すごい簡単ですね！
output <- dd_hyperbolic(data = "dd_exampleData.txt", niter = 2000, nwarmup = 1000,
nchain = 4, ncore = 4, nthin = 2)##
## Model name = dd_hyperbolic
## Data file = dd_exampleData.txt
##
## Details:
## # of chains = 4
## # of cores used = 4
## # of MCMC samples (per chain) = 2000
## # of burn-in samples = 1000
## # of subjects = 20
## # of (max) trials per subject = 108
##
##
## ****************************************
## ** Use VB estimates as initial values **
## ****************************************
## Chain 1: ------------------------------------------------------------
## Chain 1: EXPERIMENTAL ALGORITHM:
## Chain 1: This procedure has not been thoroughly tested and may be unstable
## Chain 1: or buggy. The interface is subject to change.
## Chain 1: ------------------------------------------------------------
## Chain 1:
## Chain 1:
## Chain 1:
## Chain 1: Gradient evaluation took 0.001102 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 11.02 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Begin eta adaptation.
## Chain 1: Iteration: 1 / 250 [ 0%] (Adaptation)
## Chain 1: Iteration: 50 / 250 [ 20%] (Adaptation)
## Chain 1: Iteration: 100 / 250 [ 40%] (Adaptation)
## Chain 1: Iteration: 150 / 250 [ 60%] (Adaptation)
## Chain 1: Success! Found best value [eta = 10] earlier than expected.
## Chain 1:
## Chain 1: Begin stochastic gradient ascent.
## Chain 1: iter ELBO delta_ELBO_mean delta_ELBO_med notes
## Chain 1: 100 -3637.296 1.000 1.000
## Chain 1: 200 -1855.146 0.980 1.000
## Chain 1: 300 -1686.434 0.687 0.961
## Chain 1: 400 -5003.714 0.681 0.961
## Chain 1: 500 -1530.761 0.998 0.961
## Chain 1: 600 -1406.043 0.847 0.961
## Chain 1: 700 -1170.984 0.755 0.663
## Chain 1: 800 -997.959 0.682 0.663
## Chain 1: 900 -1012.172 0.608 0.201
## Chain 1: 1000 -1075.240 0.553 0.201
## Chain 1: 1100 -1041.363 0.456 0.173
## Chain 1: 1200 -1027.903 0.361 0.100
## Chain 1: 1300 -1052.742 0.354 0.089
## Chain 1: 1400 -1065.566 0.289 0.059
## Chain 1: 1500 -1047.460 0.063 0.033
## Chain 1: 1600 -1018.781 0.057 0.028
## Chain 1: 1700 -990.140 0.040 0.028
## Chain 1: 1800 -1074.774 0.031 0.028
## Chain 1: 1900 -1006.949 0.036 0.029
## Chain 1: 2000 -1116.032 0.040 0.029
## Chain 1: 2100 -1024.234 0.046 0.029
## Chain 1: 2200 -995.219 0.047 0.029
## Chain 1: 2300 -1113.816 0.056 0.067
## Chain 1: 2400 -1010.520 0.065 0.079
## Chain 1: 2500 -1018.333 0.064 0.079
## Chain 1: 2600 -1165.126 0.073 0.090
## Chain 1: 2700 -981.910 0.089 0.098
## Chain 1: 2800 -1016.077 0.085 0.098
## Chain 1: 2900 -980.237 0.082 0.098
## Chain 1: 3000 -1057.570 0.079 0.090
## Chain 1: 3100 -1025.682 0.073 0.073
## Chain 1: 3200 -987.209 0.074 0.073
## Chain 1: 3300 -990.753 0.064 0.039
## Chain 1: 3400 -986.857 0.054 0.037
## Chain 1: 3500 -1059.822 0.060 0.039
## Chain 1: 3600 -1008.952 0.053 0.039
## Chain 1: 3700 -1012.091 0.034 0.037
## Chain 1: 3800 -1044.237 0.034 0.037
## Chain 1: 3900 -999.452 0.035 0.039
## Chain 1: 4000 -996.331 0.028 0.031
## Chain 1: 4100 -986.455 0.026 0.031
## Chain 1: 4200 -991.453 0.022 0.010
## Chain 1: 4300 -999.335 0.023 0.010
## Chain 1: 4400 -992.752 0.023 0.010
## Chain 1: 4500 -1006.534 0.018 0.010
## Chain 1: 4600 -986.990 0.014 0.010
## Chain 1: 4700 -978.329 0.015 0.010
## Chain 1: 4800 -1021.563 0.016 0.010
## Chain 1: 4900 -988.887 0.015 0.010
## Chain 1: 5000 -986.848 0.015 0.010
## Chain 1: 5100 -990.629 0.014 0.009 MEDIAN ELBO CONVERGED
## Chain 1:
## Chain 1: Drawing a sample of size 1000 from the approximate posterior...
## Chain 1: COMPLETED.## Warning: Pareto k diagnostic value is 2.12. Resampling is disabled. Decreasing
## tol_rel_obj may help if variational algorithm has terminated prematurely.
## Otherwise consider using sampling instead.##
## ************************************
## **** Model fitting is complete! ****
## ************************************パラメータ推定の収束
推定できたので，早速，収束を確認します。hBayesDMのplot関数やrhat()を使って，簡単にトレースプロットや\(\hat{R}\)が確認できます。収束は問題なさそうですね。
plot(output, type = 'trace')
rhat(output)## # A tibble: 65 x 1
## Rhat
## <dbl>
## 1 0.999
## 2 1.00
## 3 1.00
## 4 1.01
## 5 0.999
## 6 0.999
## 7 1.00
## 8 1.00
## 9 0.999
## 10 0.998
## # … with 55 more rows推定されたパラメータ
hBayesDMのplot関数は，デフォルトの場合は，ハイパーパラメータの事後分布をプロットします。
plot(output)
もし，個人ごとのパラメータを見たい場合は，plotInd()が便利です。ここでは，参加者20名のkの事後分布を示します。
plotInd(output,"k")
モデル適合
hBayesDMでは，printFitで情報量基準も出してくれます。デフォルトはLOOICで，指定すればWAICも出力してくれます。
printFit(output)## # A tibble: 1 x 3
## Model LOOIC `LOOIC Weights`
## <chr> <dbl> <dbl>
## 1 dd_hyperbolic 1883. 1printFit(output, ic="waic")## # A tibble: 1 x 3
## Model WAIC `WAIC Weights`
## <chr> <dbl> <dbl>
## 1 dd_hyperbolic 1871. 1最後に
いかがでしたか？思った以上に簡単にベイジアン認知モデリングが体験できたのではないでしょうか？hBayesDMで扱っている実験課題を研究で用いている人は，ぜひとも試してみてください！なお，hBayesDMパッケージを使った場合は，Ahn博士の論文3を引用ください。
Enjoy!
2011年くらいから，なんとなく認知モデリングで階層ベイズを使おうと思っていました。しかし，何のとっかかりがないので，ただボーっと過ごしていました。たまたま2013年に，Ahn博士のモデルベースfMRI解析における階層ベイズ論文を読んで，「これだっ！」と思って直接メールをしました。そのメールのやりとりをしていたら，Ahn博士がStanの存在を教えてくれました。Stanにいち早く気づかせてくださったAhn博士に心から感謝しています。ただ，そこから，またぼーっとしてしまって，とてもゆっくりStanの勉強を進めています。ノロノロとやっているうちに詳しい人が周りに多くなってきました。ありがたいですね。↩︎
Mazur, J. E. (1987). An adjustment procedure for studying delayed reinforcement.↩︎
Ahn, W.-Y., Haines, N., & Zhang, L. (2017). Revealing neuro-computational mechanisms of reinforcement learning and decision-making with the hBayesDM package. Computational Psychiatry. 1:1. https://doi.org/10.1101/064287↩︎