実験参加者内要因の条件間でパラメータを比較する(パラメータ間の相関を考慮した比較)
はじめに
本記事は、Stan Advent Calendar 2019 22日目の記事です。
本記事の内容はたのしいベイズモデリング 事例で拓く研究のフロンティア第15章の内容とほぼ同じですが、ちょっとだけ発展的な内容を含んでいます。
次のような手続きで実験を行ったとしましょう。PCのディスプレイ上に、「O」と「D」がランダムに呈示されます。実験参加者は、
- 「O」が出たと思ったら右のボタンを押す
- 「D」が出たと思ったら左のボタンを押す
という風に反応することを求められます。刺激が呈示されてから、実験参加者がボタンを押すまでの反応時間が記録されます。一見簡単に見えますが、これをずーーーっと続けていると、時折間違えた反応をしてしまうこともあるでしょう。また、反応が徐々に遅くなっていくと予想されます。すなわち反応の正確性や反応時間は、「情報処理が円滑に行われている程度」を反映する指標と考えられます。
モデリングの方針
このように、実験参加者に二通りの反応が許されているとき(今回は、左右いずれかのボタンを押す)、その結果として得られる反応時間のデータを、drift diffusion modelによってモデリングしてみましょう。そのようなモデリングを行う利点については、Stan Advent Calendar 2019 16日目の記事を参照してください。
本記事の目的
複数の条件を、異なる実験参加者グループが経験している場合(実験参加者「間」計画)と、同じ実験参加者が経験している場合(実験参加者「内」計画)で、モデリングの仕方が変わります。
実験参加者「内」計画の場合、同じ実験参加者がすべての条件を経験するので、条件同士で反応に相関関係が認められると予想されます。例えば反応時間を測定する実験の場合、慎重な人はすべての条件で反応が遅いかもしれないし、せっかちな人はすべての条件で反応が早いかもしれません。
本記事ではこのような実験参加者内計画における反応時間のデータを、条件間の相関を考慮したうえでモデリングする方法を紹介します。
実践
サンプルデータの作成
drift diffusion modelでは、反応時間がWiener分布に従うと考えます。そこで、Wiener分布に従う乱数を発生させましょう。brms
パッケージには、brms::rwiener()
という関数が用意されています。これは内部ではRWiener
というパッケージを利用しています。もし以下のコードを実行したとき「RWiener
というパッケージは無い」というエラーが出た場合には、追加でRWiener
パッケージをインストールしておいてください。
# load packages ---------------------------- library(tidyverse) library(brms) library(rstan) # sample data ------------------------------ trial = 1000 # 一つの条件における試行数 N_subject = 5 # 実験参加者数 tau = 0.5 # Wiener分布の非決定時間パラメータ beta = 0.5 # Wiener分布の開始点パラメータ set.seed(123) # 以下は順に、実験参加者1の条件1、条件2、実験参加者2の条件1、...、実験参加者5の条件2 dat_all = rbind(brms::rwiener(n = trial, alpha = 0.6, tau = tau, beta = beta, delta = 2), brms::rwiener(n = trial, alpha = 0.75, tau = tau, beta = beta, delta = 2.5), brms::rwiener(n = trial, alpha = 0.5, tau = tau, beta = beta, delta = 1.7), brms::rwiener(n = trial, alpha = 0.55, tau = tau, beta = beta, delta = 2.4), brms::rwiener(n = trial, alpha = 0.65, tau = tau, beta = beta, delta = 2.2), brms::rwiener(n = trial, alpha = 0.8, tau = tau, beta = beta, delta = 2.5), brms::rwiener(n = trial, alpha = 0.5, tau = tau, beta = beta, delta = 1.8), brms::rwiener(n = trial, alpha = 0.6, tau = tau, beta = beta, delta = 2.2), brms::rwiener(n = trial, alpha = 0.7, tau = tau, beta = beta, delta = 1.6), brms::rwiener(n = trial, alpha = 0.75, tau = tau, beta = beta, delta = 2.1)) %>% dplyr::mutate(ID = rep(x = 1:N_subject, each = 2 * trial), condition = rep(1:2, each = trial, time = N_subject))
Wiener分布は4つのパラメータを持ちます。このうち、非決定時間tau
と、情報蓄積の開始点beta
については、今回はすべて共通としています。
反応が起こるために必要な情報蓄積の閾値alpha
と、情報蓄積の速さdelta
が、実験参加者や条件によって異なるサンプルデータです。このサンプルデータには、以下のように4つの変数があります。
- q
- 反応時間(秒)
- resp
- 右のボタンを押したか(1)、左のボタンを押したか(0)
- ID
- 実験参加者のID。5人いるので、1~5。
- condition
- 「O」条件(1)、「D」条件(2)
head(dat_all) q resp ID condition 1 0.5569206 1 1 1 2 0.6966415 0 1 1 3 0.5378940 1 1 1 4 0.5503652 1 1 1 5 0.6447536 1 1 1 6 0.7017261 0 1 1 tail(dat_all) q resp ID condition 9995 0.6605005 0 5 2 9996 0.6017658 1 5 2 9997 0.6058292 1 5 2 9998 0.5311541 1 5 2 9999 0.7406046 1 5 2 10000 0.5725912 0 5 2
Stanコード
Stanにはwiener_lpdf()
という分布関数が用意されているので、基本的にはここに目的変数とパラメータを書き入れるだけです。なお16日目の記事とはやや違う書き方をしています。
パラメータ名は以下のように名付けなおしています。
- 閾値
alpha
→ boundary parametera
- 非決定時間
tau
→t_raw
及びt_new
- これら二つの違いは後述
- 情報蓄積の速さ
delta
→ drift rate parameterd
data { int<lower=0> N; //全データ数 int<lower=0> N_ID; //被験者数 int<lower=0> N_cond; //被験者内要因の条件数 int<lower=0> ID[N]; //被験者を識別するIDが格納されたベクトル int<lower=0> cond[N]; //被験者内条件を識別するラベルが格納されたベクトル vector<lower=0>[N] RT; //反応時間ベクトル real<lower=0> minRT; //全被験者のなかでの最小RT(推定されたtがこれを超えるとエラーになるため) int Resp[N]; //反応のベクトル real<lower=0, upper = 1> beta; //starting point (0.5に固定) } parameters { vector[N_cond] a[N_ID]; //各条件のboundary parameterを計算する材料。被験者内条件数の長さで、そのベクトルが被験者数ある。 vector[N_cond] mean_mu_a; //各条件のboundary parameterが従う、多変量正規分布の平均ベクトル corr_matrix[N_cond] omega_mu_a; //各条件のboundary parameterが従う、多変量正規分布の、分散共分散行列を作るための、相関行列 vector<lower=0>[N_cond] tau_mu_a; //各条件のboundary parameterが従う、多変量正規分布の、分散共分散行列を作るための、標準偏差ベクトル real t_raw; //全被験者に共通するnon-decision time parameterを求めるための材料 vector[N_cond] d[N_ID]; //各条件のdrift rate parameterを計算する材料。被験者内条件数の長さで、そのベクトルが被験者数ある。 vector[N_cond] mean_mu_d; //各条件のdrift rate parameterが従う、多変量正規分布の平均ベクトル corr_matrix[N_cond] omega_mu_d; //各条件のdrift rate parameterが従う、多変量正規分布の、分散共分散行列を作るための、相関行列 vector<lower=0>[N_cond] tau_mu_d; //各条件のdrift rate parameterが従う、多変量正規分布の、分散共分散行列を作るための、標準偏差ベクトル } transformed parameters{ real<lower=0, upper=1> t_new; //全被験者に共通するnon-decision time parameter cov_matrix[N_cond] mu_a_covmatrix; //各条件のboundary parameterが従う、多変量正規分布の、分散共分散行列 cov_matrix[N_cond] mu_d_covmatrix; //各条件のdrift rate parameterが従う、多変量正規分布の、分散共分散行列 t_new = inv_logit(t_raw) * minRT; //ロジスティック変換で0~1に変換し、それに最小RTを掛けているため、 //Wiener分布から推定するt_newは最小RTを超えないように工夫 mu_a_covmatrix = quad_form_diag(omega_mu_a, tau_mu_a); //各条件のboundary parameterが従う、多変量正規分布の、分散共分散行列 mu_d_covmatrix = quad_form_diag(omega_mu_d, tau_mu_d); //各条件のdrift rate parameterが従う、多変量正規分布の、分散共分散行列 } model { // model ---------------------------------------------------- for(n in 1:N){ if(Resp[n] == 1){ target += wiener_lpdf(RT[n] | exp(a[ID[n], cond[n]]), t_new, beta, exp(d[ID[n], cond[n]])); //exp(a_raw), exp(v_raw)により、両parameterが多変量「対数」正規分布から生成されたことを表現 } else{ target += wiener_lpdf(RT[n] | exp(a[ID[n], cond[n]]), t_new, beta, -exp(d[ID[n], cond[n]])); } } // boundary parameter target += multi_normal_lpdf(a | mean_mu_a, mu_a_covmatrix); // drift rate parameter target += multi_normal_lpdf(d | mean_mu_d, mu_d_covmatrix); // priors ----------------------------------------------------- target += normal_lpdf(mean_mu_a | 0, 5); target += lkj_corr_lpdf(omega_mu_a | 2); target += student_t_lpdf(tau_mu_a | 3, 0, 3) - student_t_lccdf(0 | 3, 0, 3); target += normal_lpdf(t_raw | 0, 5); target += normal_lpdf(mean_mu_d | 0, 5); target += lkj_corr_lpdf(omega_mu_d | 2); target += student_t_lpdf(tau_mu_d | 3, 0, 3) - student_t_lccdf(0 | 3, 0, 3); }
ポイントは以下の点です。
wiener_lpdf()
のなかで、a
とd
のパラメータに指数変換を施している- これらのパラメータは正の値を取ります。ということは、両パラメータに
lower=0
と下限を指定するだけでもいいような気がしますが、そうすると収束が不安定になりがちでした。 a
やd
自体には下限を指定せず、exp()
によって正の値に変換しています。
- これらのパラメータは正の値を取ります。ということは、両パラメータに
- 非決定時間パラメータ
- 反応時間の最小値を超えてしまうと推定できなくなるので、それを抑える工夫です。詳しくはコード内のコメントを見てください。
- 各条件の
a
とd
のパラメータが、多変量正規分布に従うと仮定している
1.と2.は推定効率を高めるための工夫ですが、3.が本記事の目的である、相関のあるパラメータを推定するために必要な部分です。
実験参加者「内」条件間でパラメータを比較したいとき、この多変量(対数)正規分布の平均ベクトルを比較することが出来ると思います。今回のモデルではmean_mu_a
やmean_mu_d
です。
もっとも、そういったやり方ではなく、条件間でパラメータが等しいというモデルと、パラメータが異なることを仮定するモデルを、例えばベイズファクター等を用いてモデル比較するという方法もあると思います。
推定
Rコードは以下の通りです。自己相関が高かったので、thin = 10
としています。もうちょっとモデルを工夫できるのかもしれません。
# model ------------------- sm = rstan::stan_model("ddm.stan") # data -------------------- sdata = list(N = nrow(dat_all), N_ID = dat_all %>% dplyr::distinct(ID) %>% nrow(), N_cond = dat_all %>% dplyr::distinct(condition) %>% nrow(), ID = dat_all$ID, cond = dat_all$condition, RT = dat_all$q, Resp = dat_all$resp, minRT = min(dat_all$q), z = 0.5) # sampling -------------------- fit = rstan::sampling(object = sm, data = sdata, seed = 1234, iter = 8000, warmup = 3000, thin = 10)
結果
パラメータが多いので、ここでは結果の要約は載せず、パラメータリカバリできたかどうかに注目します。推定されたa
やd
を指数変換したものが、Wiener分布のパラメータであることに注意してください。
直線上に点が載っていたら、パラメータの真値と推定結果(パラメータの事後期待値)が一致していることを意味します。閾値a
についてはうまくリカバリできたようです。
しかし情報蓄積の速さd
については、真値に近い値は得られているものの、ちょっと誤差があります。このパラメータの推定はもうちょっと工夫が要るのかもしれません。
おわりに
今回は2肢選択事態における反応時間データを、drift diffusion modelによってモデリングしました。さらにその際に、実験参加者内要因、つまり同じ実験参加者が全ての実験条件を経験したことにより、条件間でパラメータに相関が仮定される状況を取り上げました。
実験心理学ではこのような実験計画を組むことが多いので、参考になれば幸いです。
Enjoy!