実験参加者内要因の条件間でパラメータを比較する(パラメータ間の相関を考慮した比較)

はじめに

本記事は、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 parameter a
  • 非決定時間tau → t_raw及びt_new
    • これら二つの違いは後述
  • 情報蓄積の速さdelta → drift rate parameter d
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;       //ロジスティック変換で01に変換し、それに最小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);
}

ポイントは以下の点です。

  1. wiener_lpdf()のなかで、adのパラメータに指数変換を施している
    • これらのパラメータは正の値を取ります。ということは、両パラメータにlower=0と下限を指定するだけでもいいような気がしますが、そうすると収束が不安定になりがちでした。
    • ad自体には下限を指定せず、exp()によって正の値に変換しています。
  2. 非決定時間パラメータ
    • 反応時間の最小値を超えてしまうと推定できなくなるので、それを抑える工夫です。詳しくはコード内のコメントを見てください。
  3. 各条件のadのパラメータが、多変量正規分布に従うと仮定している
    • 今回は2条件しかないので、実質的に2変量正規分布
    • 多変量正規分布は、分散共分散行列をパラメータに持つので、ここで条件間の相関を表現できる
    • なお今回は、adを最終的に指数変換しているため、結果的にこれらのパラメータは多変量対数正規分布に従うというモデルになっている

1.と2.は推定効率を高めるための工夫ですが、3.が本記事の目的である、相関のあるパラメータを推定するために必要な部分です。

実験参加者「内」条件間でパラメータを比較したいとき、この多変量(対数)正規分布の平均ベクトルを比較することが出来ると思います。今回のモデルではmean_mu_amean_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)

結果

パラメータが多いので、ここでは結果の要約は載せず、パラメータリカバリできたかどうかに注目します。推定されたadを指数変換したものが、Wiener分布のパラメータであることに注意してください。

f:id:das_Kino:20191222215445p:plain

直線上に点が載っていたら、パラメータの真値と推定結果(パラメータの事後期待値)が一致していることを意味します。閾値aについてはうまくリカバリできたようです。

しかし情報蓄積の速さdについては、真値に近い値は得られているものの、ちょっと誤差があります。このパラメータの推定はもうちょっと工夫が要るのかもしれません。

おわりに

今回は2肢選択事態における反応時間データを、drift diffusion modelによってモデリングしました。さらにその際に、実験参加者内要因、つまり同じ実験参加者が全ての実験条件を経験したことにより、条件間でパラメータに相関が仮定される状況を取り上げました。

実験心理学ではこのような実験計画を組むことが多いので、参考になれば幸いです。

Enjoy!