はじめに

Linear Ballistic Accumulator modelを用いた反応時間分析では，１名の参加者の反応時間データをRstanのユーザー定義関数を活用しながら解析する方法について解説を行いました。その際に，試行数は500施行という現実的ではない値を用いていました。今回は，10名が100試行反応するという，心理学実験でありえそうな状況を想定してみます。なお，以下でStanコードを載せていますが，HTML上で表記がおかしくなることがあるので，自分の環境で動かしてみたい方は，ここから適宜RmdファイルをDLください。

使用パッケージ

今回は，以下のパッケージを用います。

rm(list=ls())
library(DiagrammeR)
library(ggplot2)
library(plotly)
library(bayesplot)
library(rstan)
library(loo)
library(rtdists)
library(msm)

10名分のデータの作成

Annis et al. (2016)による Bayesian inference with Stan: A tutorial on adding custom distributionsのコードを参考にしつつ，10名の参加者が３条件の２選択肢の課題をそれぞれ100試行行ったというデータを作ってみる。

パタメータの準備

まずは，人数分のパラメータを準備します。その際に，vは条件によって異なる（変化するor条件の効果がある）という設定にします。

set.seed(123)
# 参加者数，試行数，条件数
numberSubjects <- 10
numberTrials <- 100
numberConditions <- 3

#　ハイパーパラメータの設定
A_mu = .5
b_mu = 1
t0_mu = .5
v1_mu_c1 = 2
v2_mu_c1 = 2.5
v1_mu_c2 = 3
v2_mu_c2 = 3
v1_mu_c3 = 3.5
v2_mu_c3 = 2.5
s = 1

#　ハイパーパラメータをもとにした個人のパラメータ（切断正規分布を利用）
A = rtnorm(numberSubjects,A_mu,.5,0,Inf)
b = rtnorm(numberSubjects,b_mu,.5,A,Inf)
t0 = rtnorm(numberSubjects,t0_mu,.5,0,1)
v1_c1 = rtnorm(numberSubjects,v1_mu_c1,1,0,Inf)
v2_c1 = rtnorm(numberSubjects,v2_mu_c1,1,0,Inf)
v1_c2 = rtnorm(numberSubjects,v1_mu_c2,1,0,Inf)
v2_c2 = rtnorm(numberSubjects,v2_mu_c2,1,0,Inf)
v1_c3 = rtnorm(numberSubjects,v1_mu_c3,1,0,Inf)
v2_c3 = rtnorm(numberSubjects,v2_mu_c3,1,0,Inf)

# パラメータとハイパーパラメータのデータフレーム化
hyperParameters <- data.frame(A_mu,b_mu,t0_mu,v1_mu_c1,v2_mu_c1,v1_mu_c2,v2_mu_c2,v1_mu_c3,v2_mu_c3,s)
print(hyperParameters)

##   A_mu b_mu t0_mu v1_mu_c1 v2_mu_c1 v1_mu_c2 v2_mu_c2 v1_mu_c3 v2_mu_c3 s
## 1  0.5    1   0.5        2      2.5        3        3      3.5      2.5 1

parameters <- data.frame(A,b,t0,v1_c1,v2_c1,v1_c2,v2_c2,v1_c3,v2_c3)
print(parameters)

##            A         b        t0     v1_c1    v2_c1    v1_c2     v2_c2
## 1  0.2197622 1.1799069 0.2330341 2.8215811 3.707962 2.774229 1.9282088
## 2  0.3849113 1.2003857 0.4659625 2.6886403 1.376891 4.516471 3.3035286
## 3  1.2793542 1.6458696 0.2659726 2.5539177 2.097115 1.451247 3.4482098
## 4  0.5352542 1.0553414 0.8578277 1.9380883 2.033345 3.584614 3.0530042
## 5  0.5646439 0.7220794 0.8100644 1.6940373 3.279965 3.123854 3.9222675
## 6  1.3575325 1.8806917 0.4422001 1.6195290 2.416631 3.215942 5.0500847
## 7  0.7304581 1.8934566 0.7989248 1.3052930 2.753319 3.379639 2.5089688
## 8  1.1120409 1.7147624 0.1218993 1.7920827 2.471453 2.497677 0.6908311
## 9  0.1565736 1.2489252 0.5609480 0.7346036 2.457130 2.666793 4.0057385
## 10 0.2771690 1.3506780 0.2065314 4.1689560 3.868602 1.981425 2.2907992
##       v1_c3    v2_c3
## 1  2.811991 2.279513
## 2  4.525571 2.831782
## 3  3.215227 3.596839
## 4  2.279282 2.935181
## 5  3.681303 2.174068
## 6  3.361109 3.648808
## 7  3.505764 3.493504
## 8  3.885280 3.048397
## 9  3.129340 2.738732
## 10 4.144377 1.872094

反応と反応時間データの作成

上記で設定したパラメータをもとに，rtdistsのrLBA関数を用いて反応時間と反応のデータを作ります。作ったデータは，ロング形式にしました。そして，stanLongDataにStanで使えるようにリスト形式でデータを保存しています。

longData <- NULL
for(i in 1:numberSubjects){
    id <- rep(i,numberTrials)
    # 条件１
    tempRt <- rLBA(numberTrials, A=parameters$A[i], b=parameters$b[i], t0=parameters$t0[i], mean_v=c(parameters$v1_c1[i],parameters$v2_c1[i]), sd_v=c(s,s))
    condition <- rep(1,numberTrials)
    indDataC1 <- data.frame(id,condition,tempRt)
    
    # 条件2
    tempRt <- rLBA(numberTrials, A=parameters$A[i], b=parameters$b[i], t0=parameters$t0[i], mean_v=c(parameters$v1_c2[i],parameters$v2_c2[i]), sd_v=c(s,s))
    condition <- rep(2,numberTrials)
    indDataC2 <- data.frame(id,condition,tempRt)
    
    # 条件3
    tempRt <- rLBA(numberTrials, A=parameters$A[i], b=parameters$b[i], t0=parameters$t0[i], mean_v=c(parameters$v1_c3[i],parameters$v2_c3[i]), sd_v=c(s,s))
    condition <- rep(3,numberTrials)
    indDataC3 <- data.frame(id,condition,tempRt)
    indData <- rbind(indDataC1,indDataC2,indDataC3)
    
    # ロングデータ
    longData <- rbind(longData,indData)
}

numberAllTrials <- length(longData$rt)
# Stan用データの作成
stanLongData <- list(rt=longData$rt,response=longData$response,id=longData$id,condition = longData$condition, numberConditions=numberConditions,numberSubjects=numberSubjects,numberChoices=2,numberAllTrials=numberAllTrials)

データの確認

longDataは，id，condition,rt,responseからなります。idは10名，conditionは3条件です。

head(longData)

##   id condition        rt response
## 1  1         1 0.4637995        1
## 2  1         1 0.4430283        1
## 3  1         1 0.6445698        1
## 4  1         1 0.4954407        2
## 5  1         1 0.5488672        2
## 6  1         1 0.5467491        2

全員の結果を確認するわけにもいかないので，適当に，ID=1の条件１，２，３の結果を見てみることにします。

# 参加者IDが1の人の各条件のデータを作る。
dataSub1 <- subset(longData,longData$id==1)
dataSub1Con1 <- subset(dataSub1,dataSub1$condition==1)
dataSub1Con2 <- subset(dataSub1,dataSub1$condition==2)
dataSub1Con3 <- subset(dataSub1,dataSub1$condition==3)

各条件の選択肢別の反応数

# 条件１
table(dataSub1Con1$response)

## 
##  1  2 
## 31 69

# 条件２
table(dataSub1Con2$response)

## 
##  1  2 
## 73 27

# 条件３
table(dataSub1Con3$response)

## 
##  1  2 
## 63 37

条件１の選択肢ごとの反応時間のヒストグラム

# 選択肢1の反応時間
dataSub1Con1Choice1 <- subset(dataSub1Con1,dataSub1Con1$response==1)
p1 <- ggplot(dataSub1Con1Choice1,aes(x = rt)) + geom_histogram()+coord_cartesian(xlim = c(0, 3), ylim = c(0, 10))

# 選択肢2の反応時間
dataSub1Con1Choice2 <- subset(dataSub1Con1,dataSub1Con1$response==2)
p2 <- ggplot(dataSub1Con1Choice2,aes(x = rt)) + geom_histogram()+coord_cartesian(xlim = c(0, 3), ylim = c(0, 10))
subplot(p1, p2)

条件２の選択肢ごとの反応時間のヒストグラム

# 選択肢1の反応時間
dataSub1Con2Choice1 <- subset(dataSub1Con2,dataSub1Con2$response==1)
p1 <- ggplot(dataSub1Con2Choice1,aes(x = rt)) + geom_histogram()+coord_cartesian(xlim = c(0, 3), ylim = c(0, 10))

# 選択肢2の反応時間
dataSub1Con2Choice2 <- subset(dataSub1Con2,dataSub1Con2$response==2)
p2 <- ggplot(dataSub1Con2Choice2,aes(x = rt)) + geom_histogram()+coord_cartesian(xlim = c(0, 3), ylim = c(0, 10))
subplot(p1, p2)

条件３の選択肢ごとの反応時間のヒストグラム

# 選択肢1の反応時間
dataSub1Con3Choice1 <- subset(dataSub1Con3,dataSub1Con3$response==1)
p1 <- ggplot(dataSub1Con3Choice1,aes(x = rt)) + geom_histogram()+coord_cartesian(xlim = c(0, 3), ylim = c(0, 10))

# 選択肢2の反応時間
dataSub1Con3Choice2 <- subset(dataSub1Con3,dataSub1Con3$response==2)
p2 <- ggplot(dataSub1Con3Choice2,aes(x = rt)) + geom_histogram()+coord_cartesian(xlim = c(0, 3), ylim = c(0, 10))
subplot(p1, p2)

階層ベイズモデル

　今回は，10名の参加者が３つの条件にたいしてそれぞれ100試行反応するというデータを作成しました。これを階層ベイズを使って推定します。参加者ごとにLinear Ballistic Accumulator Modelのパラメータがありますが，それらが切断正規分布に従うという事前分布をおいて，階層ベイズによる推定を行います。また，同じ参加者が３つの条件を経験しているので，こちらも階層化しますが，条件の効果がパラメータにかかる場合とかからない場合でモデルを組みます。　条件の効果は，vのみにかかることにします。おそらく，無判断時間や閾値や開始点の上限は，条件によって変わらないのではないかという前提のものとで，vに条件の効果がかかるかどうかを検討します（実際は，もう少し複数のモデルを検討しても良いかと思います）。

条件による効果なしモデル

まず，k,A,v,tauは参加者ごとに異なるが，条件によってパラメータは異ならないという前提で解析を行います。グラフィカルモデルを描くと以下のようになる（Conditionは不要ですが，比較のために描いています）。

グラフィカルモデル

grViz("
  digraph dot {
    graph [splines = line,compound = true, nodesep = .5, ranksep = .25,
           color = black, label='Hierarchical Bayes Model of Linear Ballistic Accumulator']

      node [shape = circle,style = filled,fillcolor = white,color = black,label = 'k@_{n}'] k
      node [label = '&mu;@_{k}'] mk
      node [label = '&sigma;@_{k}'] sk
      node [label = 'A@_{n}'] A
      node [label = '&mu;@_{A}'] mA
      node [label = '&sigma;@_{A}'] sA
      node [label = '&tau;@_{n}'] t
      node [label = '&mu;@_{&tau;}'] mt
      node [label = '&sigma;@_{&tau;}'] st
      node [label = 'v[1]@_{n}'] v1
      node [label = '&mu;@_{v[1]}'] mv1
      node [label = '&sigma;@_{v[1]}'] sv1
      node [label = 'v[2]@_{n}'] v2
      node [label = '&mu;@_{v[2]}'] mv2
      node [label = '&sigma;@_{v[2]}'] sv2
      node [label = 'S'] s
      node [label = 'LBA@_{n,c,t}'] lba
      node [fillcolor = grey,label = 'RT@_{n,c,t}'] rt

        subgraph cluster1 {
          labelloc=b
          label = 'Participants n = 1...N'
            subgraph cluster2 {
              labelloc=b
              label = 'Condition c = 1...C'
                subgraph cluster3 {
                  labelloc=b
                  label = 'Trials t = 1...T'
                    edge [color = black]
                      lba -> rt
                }
              edge [color = black]
            }
          edge [color = black]
            k -> lba
            t -> lba
            A -> lba
            v1 -> lba
            v2 -> lba
        }
      edge [color = black]
        mk -> k
        sk -> k
        mA -> A
        sA -> A
        mt -> t
        st -> t
        mv1 -> v1
        sv1 -> v1
        mv2 -> v2
        sv2 -> v2 [tailport=s,headport=e]
        s -> lba
  }",engine = "dot")

事前分布はAnnisらを参考に以下にしました(ハイパーパラメータの分散の事前分布のみGamma(1,1)からstudent T(4,0,1変更）。 \[ k_{n} \sim Normal(\mu_{k},\sigma_{k})T(0,\infty)\\ A_{n} \sim Normal(\mu_{A},\sigma_{A})T(0,\infty)\\ \tau_{n} \sim Normal(\mu_{\tau},\sigma_{\tau})T(0,\infty)\\ v[1]_{n} \sim Normal(\mu_{v[1]},\sigma_{v[1]})T(0,\infty)\\ v[2]_{n} \sim Normal(\mu_{v[2]},\sigma_{v[2]})T(0,\infty)\\ \mu_{A},\mu_{k} \sim Normal(0.5,1)T(0,\infty)\\ \mu_{\tau} \sim Normal(0.5,0.5)T(0,\infty)\\ \mu_{v[1]},\mu_{v[2]} \sim Normal(2,1)T(0,\infty)\\ \sigma_{A},\sigma_{k}, \sigma_{v[1]},\sigma_{v[2]},\sigma_{\tau} \sim student\,t(4,0,1)T(0,\infty)\\ S = 1 \hspace{30pt}※Sは定数\\ RT[n,c, t] \sim LBA(k_{n},A_{n},\tau_{n},v[1]_{n},s) [選択肢１の時]\\ RT[n,c, t] \sim LBA(k_{n},A_{n},\tau_{n},v[2]_{n},s) [選択肢2の時]\\ \]

Stanコード

基本的には，「RStanによる反応時間の解析：Linear Ballistic Accumulator modelを用いて」で解説したユーザー定義関数を用います。

functions{
     
     real lba_pdf(real t, real b, real A, real v, real s){
          //PDF of the LBA model
          
          real b_A_tv_ts;
          real b_tv_ts;
          real term_1;
          real term_2;
          real term_3;
          real term_4;
          real pdf;
          
          b_A_tv_ts = (b - A - t*v)/(t*s);
          b_tv_ts = (b - t*v)/(t*s);
          term_1 = v*Phi(b_A_tv_ts);
          term_2 = s*exp(normal_lpdf(b_A_tv_ts|0,1));
          term_3 = v*Phi(b_tv_ts);
          term_4 = s*exp(normal_lpdf(b_tv_ts|0,1)); 
          pdf = (1 / A)*(-term_1 + term_2 + term_3 - term_4);
          
          return pdf;
     }
     
     real lba_cdf(real t, real b, real A, real v, real s){
          //CDF of the LBA model
          
          real b_A_tv;
          real b_tv;
          real ts;
          real term_1;
          real term_2;
          real term_3;
          real term_4;
          real cdf; 
          
          b_A_tv = b - A - t*v;
          b_tv = b - t*v;
          ts = t*s;
          term_1 = b_A_tv/A * Phi(b_A_tv/ts);   
          term_2 = b_tv/A   * Phi(b_tv/ts);
          term_3 = ts/A     * exp(normal_lpdf(b_A_tv/ts|0,1)); 
          term_4 = ts/A     * exp(normal_lpdf(b_tv/ts|0,1)); 
          cdf = 1 + term_1 - term_2 + term_3 - term_4;
          
          return cdf;
          
     }
     
     real lba_lpdf(real rt, real res, real k, real A, vector v, real s, real tau){
          
          real t;
          real b;
          real cdf;
          real pdf;     
          real prob;
          real out;
          real prob_neg;

          b = A + k;
          t = rt - tau;
          
          if(t > 0){            
                cdf = 1;
                for(j in 1 : num_elements(v)){
                      if(res == j){
                            pdf = lba_pdf(t, b, A, v[j], s);
                      }else{    
                            cdf = (1 - lba_cdf(t, b, A, v[j], s)) * cdf;
                      }
                }
                
                prob_neg = 1;
                for(j in 1 : num_elements(v)){
                      prob_neg = Phi(-v[j]/s) * prob_neg;    
                }
                prob = pdf*cdf;     
                prob = prob/(1 - prob_neg); 
                if(prob < 1e-10){
                      prob = 1e-10;             
                }
          }else{
              prob = 1e-10;         
          }
          out = log(prob);
          return out;       
     }
     
    vector lba_rng(real k, real A, vector v, real s, real tau){
          
          int get_pos_drift;    
          int no_pos_drift;
          int get_first_pos;
          vector[num_elements(v)] drift;
          int max_iter;
          int iter;
          real start[num_elements(v)];
          real ttf[num_elements(v)];
          int resp[num_elements(v)];
          real rt;
          vector[2] pred;
          real b;
          
          //try to get a positive drift rate
          get_pos_drift = 1;
          no_pos_drift = 0;
          max_iter = 1000;
          iter = 0;
          while(get_pos_drift){
               for(j in 1 : num_elements(v)){
                    drift[j] = normal_rng(v[j],s);
                    if(drift[j] > 0){
                         get_pos_drift = 0;
                    }
               }
               iter = iter + 1;
               if(iter > max_iter){
                    get_pos_drift = 0;
                    no_pos_drift = 1;
               }    
          }
          //if both drift rates are <= 0
          //return an infinite response time
          if(no_pos_drift){
               pred[1] = -1;
               pred[2] = -1;
          }else{
               b = A + k;
               for(i in 1 : num_elements(v)){
                    //start time of each accumulator    
                    start[i] = uniform_rng(0, A);
                    //finish times
                    ttf[i] = (b-start[i])/drift[i];
               }
               //rt is the fastest accumulator finish time  
               //if one is negative get the positive drift
               resp = sort_indices_asc(ttf);
               ttf = sort_asc(ttf);
               get_first_pos = 1;
               iter = 1;
               while(get_first_pos){
                    if(ttf[iter] > 0){
                         pred[1] = ttf[iter] + tau;
                         pred[2] = resp[iter]; 
                         get_first_pos = 0;
                    }
                    iter = iter + 1;
               }
          }
          return pred;  
     }
}

data{
     int numberConditions;
     int numberSubjects;
     int numberChoices;
     int numberAllTrials;
     vector[numberAllTrials] rt;
     vector[numberAllTrials] response;
     int condition[numberAllTrials];
     int id[numberAllTrials];
}

parameters {
     real < lower = 0 > k[numberSubjects];
     real < lower = 0 > A[numberSubjects];
     real < lower = 0 > psi[numberSubjects];
     vector < lower = 0 > [numberChoices] v [numberSubjects];
     real < lower = 0 > k_sigma;
     real < lower = 0 > A_sigma;
     real < lower = 0 > psi_sigma;
     vector < lower = 0 > [numberChoices] v_sigma;
     real < lower = 0 > k_mu;
     real < lower = 0 > A_mu;
     real < lower = 0 > psi_mu;
     vector < lower=0> [numberChoices] v_mu;
}

transformed parameters {
     real s;
     s = 1;
}

model {
     int subIndex;
     // ハイパーパラメータの事前分布の設定
     k_mu ~ normal(.5, 1)T[0, ];
     A_mu ~ normal(.5, 1)T[0, ];
     psi_mu ~ normal(.5, .5)T[0, ];
     k_sigma ~ student_t(4, 0, 1)T[0, ];
     A_sigma ~ student_t(4, 0, 1)T[0, ];
     psi_sigma ~ student_t(4, 0, 1)T[0, ];
     
     for(n in 1 : numberChoices){
        v_mu[n] ~ normal(2, 1)T[0, ];
        v_sigma[n] ~ student_t(4, 0, 1)T[0, ];
      }
     
     // パラメータの事前分布の設定
     for(i in 1 : numberSubjects){      
          k[i] ~ normal(k_mu, k_sigma)T[0, ];
          A[i] ~ normal(A_mu, A_sigma)T[0, ];
          psi[i] ~ normal(psi_mu, psi_sigma)T[0, ];

          for(n in 1 : numberChoices){
              v[i,n] ~ normal(v_mu[n],v_sigma[n])T[0, ];
          }
     }
    
    // 対数事後確率の計算部分
    for(m in 1 : numberAllTrials){
           subIndex = id[m];
          rt[m] ~ lba(response[m], k[subIndex], A[subIndex], v[subIndex], s, psi[subIndex]);
    }
}

generated quantities {
     vector[2] pred[numberSubjects];
     vector[numberAllTrials] log_lik;
     int subIndex;
     for(i in 1 : numberSubjects){
          pred[i] = lba_rng(A[i], k[i], v[i], s, psi[i]);
     }
    for(m in 1 : numberAllTrials){
          subIndex = id[m];
          log_lik[m] = lba_lpdf(rt[m] | response[m], k[subIndex], A[subIndex], v[subIndex], s, psi[subIndex]);
    }
}

条件によるAとvへの効果ありモデル

次は，条件によるv（ドリフト率の平均）への効果を仮定したモデルになります。おそらく条件によって異なるという前提がある場合は，こちらのほうが自然に思います。

グラフィカルモデル

grViz("
  digraph dot {
    graph [splines = line,compound = true, nodesep = .5, ranksep = .25,
           color = black, label='Hierarchical Bayes Model of Linear Ballistic Accumulator']

      node [shape = circle,style = filled,fillcolor = white,color = black,label = 'k@_{n}'] k
      node [label = '&mu;@_{k}'] mk
      node [label = '&sigma;@_{k}'] sk
      node [label = 'A@_{n,c}'] A
      node [label = '&mu;@_{A}'] mA
      node [label = '&sigma;@_{A}'] sA
      node [label = '&tau;@_{n}'] t
      node [label = '&mu;@_{&tau;}'] mt
      node [label = '&sigma;@_{&tau;}'] st
      node [label = 'v[1]@_{n,c}'] v1
      node [label = '&mu;@_{v[1]}'] mv1
      node [label = '&sigma;@_{v[1]}'] sv1
      node [label = 'v[2]@_{n,c}'] v2
      node [label = '&mu;@_{v[2]}'] mv2
      node [label = '&sigma;@_{v[2]}'] sv2
      node [label = 'S'] s
      node [label = 'LBA@_{n,c,t}'] lba
      node [fillcolor = grey,label = 'RT@_{n,c,t}'] rt

        subgraph cluster1 {
          labelloc=b
          label = 'Participants n = 1...N'
            subgraph cluster2 {
              labelloc=b
              label = 'Condition c = 1...C'
                subgraph cluster3 {
                  labelloc=b
                  label = 'Trials t = 1...T'
                    edge [color = black]
                      lba -> rt
                }
              edge [color = black]
                A -> lba
                v1 -> lba
                v2 -> lba
            }
          edge [color = black]
            k -> lba
            t -> lba
        }
      edge [color = black]
        mk -> k
        sk -> k
        mA -> A
        sA -> A
        mt -> t
        st -> t
        mv1 -> v1
        sv1 -> v1
        mv2 -> v2
        sv2 -> v2 [tailport=s,headport=e]
        s -> lba
  }",engine = "dot")

事前分布はAnnisらを参考に以下にしました。 \[ k_{i} \sim Normal(\mu_{k},\sigma_{k})T(0,\infty)\\ A_{i} \sim Normal(\mu_{A},\sigma_{A})T(0,\infty)\\ \tau_{i} \sim Normal(\mu_{\tau},\sigma_{\tau})T(0,\infty)\\ v1_{i} \sim Normal(\mu_{v1},\sigma_{v1})T(0,\infty)\\ v2_{i} \sim| Normal(\mu_{v2},\sigma_{v2})T(0,\infty)\\ \mu^{A},\mu^{k} \sim Normal(0.5,1)T(0,\infty)\\ \mu^{\tau} \sim Normal(0.5,0.5)T(0,\infty)\\ \mu^{v1},\mu^{v2} \sim Normal(2,1)T(0,\infty)\\ \sigma^{A},\sigma^{k}, \sigma^{v1},\sigma^{v2},\sigma^{\tau} \sim student_t(4,0,1)T(0,\infty) [Gamma(1,1)から変更]\\ S = 1 [Sは定数]\\ RT[t, i] \sim LBA(k_{i},A_{i},\tau_{i},v1_{i,j},s) [選択肢１の時]\\ RT[t, i] \sim LBA(k_{i},A_{i},\tau_{i},v2_{i,j},s) [選択肢2の時]\\ k \sim Normal(.5,1)T[0,]\\ A \sim Normal(.5,1)T[0,]\\ \tau \sim Normal(.5,.5)T[0,]\\ v[1],v[2] \sim Normal(2,1)T[0,]\\ S = 1 \hspace{30pt}※Sは定数 \]

Stanコード

functions{
     
     real lba_pdf(real t, real b, real A, real v, real s){
          //PDF of the LBA model
          
          real b_A_tv_ts;
          real b_tv_ts;
          real term_1;
          real term_2;
          real term_3;
          real term_4;
          real pdf;
          
          b_A_tv_ts = (b - A - t*v)/(t*s);
          b_tv_ts = (b - t*v)/(t*s);
          term_1 = v*Phi(b_A_tv_ts);
          term_2 = s*exp(normal_lpdf(b_A_tv_ts|0,1));
          term_3 = v*Phi(b_tv_ts);
          term_4 = s*exp(normal_lpdf(b_tv_ts|0,1)); 
          pdf = (1 / A)*(-term_1 + term_2 + term_3 - term_4);
          
          return pdf;
     }
     
     real lba_cdf(real t, real b, real A, real v, real s){
          //CDF of the LBA model
          
          real b_A_tv;
          real b_tv;
          real ts;
          real term_1;
          real term_2;
          real term_3;
          real term_4;
          real cdf; 
          
          b_A_tv = b - A - t*v;
          b_tv = b - t*v;
          ts = t*s;
          term_1 = b_A_tv/A * Phi(b_A_tv/ts);   
          term_2 = b_tv/A   * Phi(b_tv/ts);
          term_3 = ts/A     * exp(normal_lpdf(b_A_tv/ts|0,1)); 
          term_4 = ts/A     * exp(normal_lpdf(b_tv/ts|0,1)); 
          cdf = 1 + term_1 - term_2 + term_3 - term_4;
          
          return cdf;
          
     }
     
     real lba_lpdf(real rt, real res, real k, real A, vector v, real s, real tau){
          
          real t;
          real b;
          real cdf;
          real pdf;     
          real prob;
          real out;
          real prob_neg;

          b = A + k;
          t = rt - tau;
          
          if(t > 0){            
                cdf = 1;
                for(j in 1 : num_elements(v)){
                      if(res == j){
                            pdf = lba_pdf(t, b, A, v[j], s);
                      }else{    
                            cdf = (1 - lba_cdf(t, b, A, v[j], s)) * cdf;
                      }
                }
                
                prob_neg = 1;
                for(j in 1 : num_elements(v)){
                      prob_neg = Phi(-v[j]/s) * prob_neg;    
                }
                prob = pdf*cdf;     
                prob = prob/(1 - prob_neg); 
                if(prob < 1e-10){
                      prob = 1e-10;             
                }
          }else{
              prob = 1e-10;         
          }
          out = log(prob);
          return out;       
     }
     
    vector lba_rng(real k, real A, vector v, real s, real tau){
          
          int get_pos_drift;    
          int no_pos_drift;
          int get_first_pos;
          vector[num_elements(v)] drift;
          int max_iter;
          int iter;
          real start[num_elements(v)];
          real ttf[num_elements(v)];
          int resp[num_elements(v)];
          real rt;
          vector[2] pred;
          real b;
          
          //try to get a positive drift rate
          get_pos_drift = 1;
          no_pos_drift = 0;
          max_iter = 1000;
          iter = 0;
          while(get_pos_drift){
               for(j in 1 : num_elements(v)){
                    drift[j] = normal_rng(v[j],s);
                    if(drift[j] > 0){
                         get_pos_drift = 0;
                    }
               }
               iter = iter + 1;
               if(iter > max_iter){
                    get_pos_drift = 0;
                    no_pos_drift = 1;
               }    
          }
          //if both drift rates are <= 0
          //return an infinite response time
          if(no_pos_drift){
               pred[1] = -1;
               pred[2] = -1;
          }else{
               b = A + k;
               for(i in 1 : num_elements(v)){
                    //start time of each accumulator    
                    start[i] = uniform_rng(0, A);
                    //finish times
                    ttf[i] = (b-start[i])/drift[i];
               }
               //rt is the fastest accumulator finish time  
               //if one is negative get the positive drift
               resp = sort_indices_asc(ttf);
               ttf = sort_asc(ttf);
               get_first_pos = 1;
               iter = 1;
               while(get_first_pos){
                    if(ttf[iter] > 0){
                         pred[1] = ttf[iter] + tau;
                         pred[2] = resp[iter]; 
                         get_first_pos = 0;
                    }
                    iter = iter + 1;
               }
          }
          return pred;  
     }
}

data{
     int numberConditions;
     int numberSubjects;
     int numberChoices;
     int numberAllTrials;
     vector[numberAllTrials] rt;
     vector[numberAllTrials] response;
     int condition[numberAllTrials];
     int id[numberAllTrials];
}

parameters {
     real < lower = 0 > k[numberSubjects];
     real < lower = 0 > A [numberSubjects];
     real < lower = 0 > psi[numberSubjects];
     vector < lower = 0 > [numberChoices] v [numberSubjects, numberConditions];
     real < lower = 0 > k_sigma;
     real < lower = 0 > A_sigma;
     real < lower = 0 > psi_sigma;
     vector < lower = 0 > [numberChoices] v_sigma [numberConditions];
     real < lower = 0 > k_mu;
     real < lower = 0 > A_mu;
     real < lower = 0 > psi_mu;
     vector < lower=0> [numberChoices] v_mu [numberConditions];
}

transformed parameters {
     real s;
     s = 1;
}

model {
     int subIndex;
     int condIndex;
     // ハイパーパラメータの事前分布の設定
     k_mu ~ normal(.5, 1)T[0, ];
     psi_mu ~ normal(.5, .5)T[0, ];
     k_sigma ~ student_t(4, 0, 1)T[0, ];
     psi_sigma ~ student_t(4, 0, 1)T[0, ];
     A_sigma ~ student_t(4, 0, 1)T[0, ];
     A_mu ~ normal(.5, 1)T[0, ];
     
     for( j in 1 : numberConditions){
          for(n in 1 : numberChoices){
            v_mu[j, n] ~ normal(2, 1)T[0, ];
            v_sigma[j, n] ~ student_t(4, 0, 1)T[0, ];
          }
     }
      
     // パラメータの事前分布の設定
     for(i in 1 : numberSubjects){      
        k[i] ~ normal(k_mu, k_sigma)T[0, ];
        psi[i] ~ normal(psi_mu, psi_sigma)T[0, ];
        A[i] ~ normal(A_mu, A_sigma)T[0, ];
          for( j in 1 : numberConditions){
              for(n in 1 : numberChoices){
                v[i, j, n] ~ normal(v_mu[j, n],v_sigma[j, n])T[0, ];
              }
          }
     }
    
    for(m in 1 : numberAllTrials){
      subIndex = id[m];
      condIndex = condition[m];
      rt[m] ~ lba(response[m], k[subIndex], A[subIndex], v[subIndex, condIndex], s, psi[subIndex]);
    }
}

generated quantities {
     vector[2] pred[numberSubjects, numberConditions];
     vector[numberAllTrials] log_lik;
     int subIndex;
     int condIndex;
     
     for(i in 1 : numberSubjects){
          for(j in 1 : numberConditions){
               pred[i, j] = lba_rng(k[i], A[i], v[i, j], s, psi[i]);
          }
     }

    for(m in 1 : numberAllTrials){
          subIndex = id[m];
          condIndex = condition[m];
          log_lik[m] = lba_lpdf(rt[m] | response[m], k[subIndex], A[subIndex], v[subIndex, condIndex], s, psi[subIndex]);
    }
     
}

Rstanによるサンプリング

さてさて，サンプリングしてみます。まずは，条件による効果なしモデルからサンプリングを行います。

rstan_options(auto_write=TRUE)
options(mc.cores = parallel::detectCores())
hmcIter = 1000
hmcChains = 4
hmcWarmup = 300
hmcThin = 2

fit1 <- sampling(lbaModel1, 
            seed = 1234, 
            data = stanLongData,
            warmup = hmcWarmup, 
            iter = hmcIter,
            chains = hmcChains,
            thin =hmcThin)

次に，条件による効果ありモデルのサンプリングを行います。

fit2 <- sampling(lbaModel2, 
            seed = 1234, 
            data = stanLongData,
            warmup = hmcWarmup, 
            iter = hmcIter,
            chains = hmcChains,
            thin =hmcThin)

収束判定(条件による効果なしモデル)

収束してそうですね。収束判定には，bayesplotを活用しています。このページが参考になりました。なお，定数のSは除外しています。

トレースプロット

posterior11 <- rstan::extract(fit1, inc_warmup = TRUE, permuted = FALSE)
# k[1:5]
color_scheme_set("mix-blue-pink")
pk1 <- mcmc_trace(posterior11,  pars = c("k[1]", "k[2]", "k[3]", "k[4]", "k[5]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pk1 <- pk1+ facet_text(size = 15)
plot(pk1)

# k[6:10]
pk2 <- mcmc_trace(posterior11,  pars = c("k[6]", "k[7]", "k[8]", "k[9]", "k[10]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pk2 <- pk2+ facet_text(size = 15)
plot(pk2)

# A[1:5]
pA1 <- mcmc_trace(posterior11,  pars = c("A[1]", "A[2]", "A[3]", "A[4]", "A[5]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pA1 <- pA1+ facet_text(size = 15)
plot(pA1)

# A[6:10]
pA2 <- mcmc_trace(posterior11,  pars = c( "A[6]", "A[7]", "A[8]", "A[9]", "A[10]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pA2 <- pA2+ facet_text(size = 15)
plot(pA2)

# tau[1:5]
pt <- mcmc_trace(posterior11,  pars = c("psi[1]", "psi[2]", "psi[3]", "psi[4]", "psi[5]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pt <- pt+ facet_text(size = 15)
plot(pt)

# tau[6:10]
pt <- mcmc_trace(posterior11,  pars = c("psi[6]", "psi[7]", "psi[8]", "psi[9]", "psi[10]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pt <- pt+ facet_text(size = 15)
plot(pt)

# v[1:3]
pv1 <- mcmc_trace(posterior11,  pars = c("v[1,1]", "v[1,2]", "v[2,1]", "v[2,2]", "v[3,1]", "v[3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv1 <- pv1+ facet_text(size = 15)
plot(pv1)

# v[4:6]
pv2 <- mcmc_trace(posterior11,  pars = c("v[4,1]", "v[4,2]", "v[5,1]", "v[5,2]", "v[6,1]", "v[6,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv2 <- pv2+ facet_text(size = 15)
plot(pv2)

# v[7:8]
pv3 <- mcmc_trace(posterior11,  pars = c("v[7,1]", "v[7,2]", "v[8,1]", "v[8,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv3 <- pv3+ facet_text(size = 15)
plot(pv3)

# v[9:10]
pv4 <- mcmc_trace(posterior11,  pars = c("v[9,1]", "v[9,2]", "v[10,1]", "v[10,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv4 <- pv4+ facet_text(size = 15)
plot(pv4)

# hyperparameter 1
ph1 <- mcmc_trace(posterior11,  pars = c("k_mu", "A_mu", "psi_mu", "v_mu[1]", "v_mu[2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
ph1 <- ph1+ facet_text(size = 15)
plot(ph1)

# hyperparameter 2
ph2 <- mcmc_trace(posterior11,  pars = c("k_sigma", "A_sigma", "psi_sigma", "v_sigma[1]", "v_sigma[2]","lp__"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
ph2 <- ph2+ facet_text(size = 15)
plot(ph2)

\(\hat{R}\)の確認

bayesplotで確認をしても良いですが，ちょっとパラメータ多いすぎるので，一括で検討します。\(\hat{R}\)が1.10以下ならTRUE，1.10以上のものがあればFALSEになります。

all(stan_rhat(fit1)$data <=1.10,na.rm=T)

## [1] TRUE

有効サンプルサイズ

bayesplotで確認をしても良いですが，ちょっとパラメータ多いすぎるので，一括で検討します。有効サンプルサイズが0.1以上ならTRUE，0.1以下のものがあればFALSEになります。

all(stan_ess(fit1)$data >= 0.1,na.rm=T)

## [1] TRUE

自己相関

posterior12 <- as.matrix(fit1)
mcmc_acf(posterior12, par = c("k[1]","k[2]","k[3]","k[4]","k[5]"), lags = 10)

mcmc_acf(posterior12, par = c("k[6]","k[7]","k[8]","k[9]","k[10]"), lags = 10)

mcmc_acf(posterior12, par = c( "A[1]", "A[2]", "A[3]", "A[4]", "A[5]"), lags = 10)

mcmc_acf(posterior12, par = c("A[6]", "A[7]", "A[8]", "A[9]", "A[10]"), lags = 10)

mcmc_acf(posterior12, par = c("psi[1]", "psi[2]", "psi[3]", "psi[4]", "psi[5]"), lags = 10)

mcmc_acf(posterior12, par = c("psi[6]", "psi[7]", "psi[8]", "psi[9]", "psi[10]"), lags = 10)

mcmc_acf(posterior12, par = c("v[1,1]", "v[1,2]", "v[2,1]", "v[2,2]", "v[3,1]", "v[3,2]"), lags = 10)

mcmc_acf(posterior12, par = c("v[4,1]", "v[4,2]", "v[5,1]", "v[5,2]", "v[6,1]", "v[6,2]"), lags = 10)

mcmc_acf(posterior12, par = c("v[7,1]", "v[7,2]", "v[8,1]", "v[8,2]", "v[9,1]", "v[9,2]"), lags = 10)

mcmc_acf(posterior12, par = c("v[10,1]", "v[10,2]", "k_mu", "A_mu", "psi_mu", "v_mu[1]"), lags = 10)

mcmc_acf(posterior12, par = c("v_mu[2]", "k_sigma", "A_sigma", "psi_sigma", "v_sigma[1]", "v_sigma[2]","lp__"), lags = 10)

収束判定(条件による効果ありモデル)

収束してそうですね。収束判定には，bayesplotを活用しています。このページが参考になりました。なお，定数のSは除外しています。

トレースプロット

posterior21 <- rstan::extract(fit2, inc_warmup = TRUE, permuted = FALSE)
# k[1:5]
color_scheme_set("mix-blue-pink")
pk1 <- mcmc_trace(posterior21,  pars = c("k[1]", "k[2]", "k[3]", "k[4]", "k[5]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pk1 <- pk1+ facet_text(size = 15)
plot(pk1)

# k[6:10]
pk2 <- mcmc_trace(posterior21,  pars = c("k[6]", "k[7]", "k[8]", "k[9]", "k[10]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pk2 <- pk2+ facet_text(size = 15)
plot(pk2)

# A[1:5]
pA1 <- mcmc_trace(posterior21,  pars = c("A[1]", "A[2]", "A[3]", "A[4]", "A[5]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pA1 <- pA1+ facet_text(size = 15)
plot(pA1)

# A[6:10]
pA2 <- mcmc_trace(posterior21,  pars = c("A[1]", "A[2]", "A[3]", "A[4]", "A[5]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pA2 <- pA2+ facet_text(size = 15)
plot(pA2)

# tau[1:5]
pt <- mcmc_trace(posterior21,  pars = c("psi[1]", "psi[2]", "psi[3]", "psi[4]", "psi[5]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pt <- pt+ facet_text(size = 15)
plot(pt)

# tau[6:10]
pt2 <- mcmc_trace(posterior21,  pars = c("psi[6]", "psi[7]", "psi[8]", "psi[9]", "psi[10]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pt2 <- pt2+ facet_text(size = 15)
plot(pt2)

# v[1]
pv1 <- mcmc_trace(posterior21,  pars = c("v[1,1,1]", "v[1,1,2]", "v[1,2,1]", "v[1,2,2]", "v[1,3,1]", "v[1,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv1 <- pv1+ facet_text(size = 15)
plot(pv1)

# v[2]
pv2 <- mcmc_trace(posterior21,  pars = c("v[2,1,1]", "v[2,1,2]", "v[2,2,1]", "v[2,2,2]", "v[2,3,1]", "v[2,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv2 <- pv2+ facet_text(size = 15)
plot(pv2)

# v[3]
pv3 <- mcmc_trace(posterior21,  pars = c("v[3,1,1]", "v[3,1,2]", "v[3,2,1]", "v[3,2,2]", "v[3,3,1]", "v[3,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv3 <- pv3+ facet_text(size = 15)
plot(pv3)

# v[4]
pv4 <- mcmc_trace(posterior21,  pars = c("v[4,1,1]", "v[4,1,2]", "v[4,2,1]", "v[4,2,2]", "v[4,3,1]", "v[4,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv4 <- pv4+ facet_text(size = 15)
plot(pv4)

# v[5]
pv5 <- mcmc_trace(posterior21,  pars = c("v[5,1,1]", "v[5,1,2]", "v[5,2,1]", "v[5,2,2]", "v[5,3,1]", "v[5,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv5 <- pv5+ facet_text(size = 15)
plot(pv5)

# v[6]
pv6 <- mcmc_trace(posterior21,  pars = c("v[6,1,1]", "v[6,1,2]", "v[6,2,1]", "v[6,2,2]", "v[6,3,1]", "v[6,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv6 <- pv6+ facet_text(size = 15)
plot(pv6)

# v[7]
pv7 <- mcmc_trace(posterior21,  pars = c("v[7,1,1]", "v[7,1,2]", "v[7,2,1]", "v[7,2,2]", "v[7,3,1]", "v[7,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv7 <- pv7+ facet_text(size = 15)
plot(pv7)

# v[8]
pv8 <- mcmc_trace(posterior21,  pars = c("v[8,1,1]", "v[8,1,2]", "v[8,2,1]", "v[8,2,2]", "v[8,3,1]", "v[8,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv8 <- pv8+ facet_text(size = 15)
plot(pv8)

# v[9]
pv9 <- mcmc_trace(posterior21,  pars = c("v[9,1,1]", "v[9,1,2]", "v[9,2,1]", "v[9,2,2]", "v[9,3,1]", "v[9,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv9 <- pv9+ facet_text(size = 15)
plot(pv9)

# v[10]
pv10 <- mcmc_trace(posterior21,  pars = c("v[10,1,1]", "v[10,1,2]", "v[10,2,1]", "v[10,2,2]", "v[10,3,1]", "v[10,3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
pv10 <- pv10+ facet_text(size = 15)
plot(pv10)

# hyperparameter 1
ph1 <- mcmc_trace(posterior21,  pars = c("k_mu","k_sigma", "A_mu","A_sigma","psi_mu","psi_sigma","lp__"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
ph1 <- ph1+ facet_text(size = 15)
plot(ph1)

# hyperparameter 2
ph2 <- mcmc_trace(posterior21,  pars = c("v_mu[1,1]", "v_mu[1,2]", "v_mu[2,1]", "v_mu[2,2]", "v_mu[3,1]", "v_mu[3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
ph2 <- ph2+ facet_text(size = 15)
plot(ph2)

# hyperparameter 3
ph3 <- mcmc_trace(posterior21,  pars = c("v_sigma[1,1]", "v_sigma[1,2]", "v_sigma[2,1]", "v_sigma[2,2]", "v_sigma[3,1]", "v_mu[3,2]"), n_warmup = 125,facet_args = list(nrow = 2, labeller = label_parsed))
ph3 <- ph3+ facet_text(size = 15)
plot(ph3)

\(\hat{R}\)の確認

all(stan_rhat(fit2)$data <=1.10,na.rm=T)

## [1] TRUE

有効サンプルサイズ

all(stan_ess(fit2)$data >= 0.1,na.rm=T)

## [1] TRUE

自己相関

posterior22 <- as.matrix(fit2)
mcmc_acf(posterior22, par = c("k[1]", "k[2]", "k[3]", "k[4]", "k[5]"), lags = 10)

mcmc_acf(posterior22, par = c("k[6]", "k[7]", "k[8]", "k[9]", "k[10]"), lags = 10)

mcmc_acf(posterior22, par = c("A[1]","A[2]","A[3]","A[4]","A[5]"), lags = 10)

mcmc_acf(posterior22, par = c("A[6]","A[7]","A[8]","A[9]","A[10]"), lags = 10)

mcmc_acf(posterior22, par = c("psi[1]", "psi[2]", "psi[3]", "psi[4]", "psi[5]"), lags = 10)

mcmc_acf(posterior22, par = c("psi[6]", "psi[7]", "psi[8]", "psi[9]", "psi[10]"), lags = 10)

mcmc_acf(posterior22, par = c("v[1,1,1]", "v[1,1,2]", "v[1,2,1]", "v[1,2,2]", "v[1,3,1]", "v[1,3,2]"), lags = 10)

mcmc_acf(posterior22, par = c("v[2,1,1]", "v[2,1,2]", "v[2,2,1]", "v[2,2,2]", "v[2,3,1]", "v[2,3,2]"), lags = 10)

mcmc_acf(posterior22, par = c("v[3,1,1]", "v[3,1,2]", "v[3,2,1]", "v[3,2,2]", "v[3,3,1]", "v[3,3,2]"), lags = 10)

mcmc_acf(posterior22, par = c("v[4,1,1]", "v[4,1,2]", "v[4,2,1]", "v[4,2,2]", "v[4,3,1]", "v[4,3,2]"), lags = 10)

mcmc_acf(posterior22, par = c("v[5,1,1]", "v[5,1,2]", "v[5,2,1]", "v[5,2,2]", "v[5,3,1]", "v[5,3,2]"), lags = 10)

mcmc_acf(posterior22, par = c("v[6,1,1]", "v[6,1,2]", "v[6,2,1]", "v[6,2,2]", "v[6,3,1]", "v[6,3,2]"), lags = 10)

mcmc_acf(posterior22, par = c("v[7,1,1]", "v[7,1,2]", "v[7,2,1]", "v[7,2,2]", "v[7,3,1]", "v[7,3,2]"), lags = 10)

mcmc_acf(posterior22, par = c("v[8,1,1]", "v[8,1,2]", "v[8,2,1]", "v[8,2,2]", "v[8,3,1]", "v[8,3,2]"), lags = 10)

mcmc_acf(posterior22, par = c("v[9,1,1]", "v[9,1,2]", "v[9,2,1]", "v[9,2,2]", "v[9,3,1]", "v[9,3,2]"), lags = 10)

mcmc_acf(posterior22, par = c("v[10,1,1]", "v[10,1,2]", "v[10,2,1]", "v[10,2,2]", "v[10,3,1]", "v[10,3,2]"), lags = 10)

mcmc_acf(posterior22, par = c("k_mu", "k_sigma","A_mu","A_sigma","psi_mu","psi_sigma","lp__"), lags = 10)

mcmc_acf(posterior22, par = c("v_mu[1,1]", "v_mu[1,2]", "v_mu[2,1]", "v_mu[2,2]", "v_mu[3,1]", "v_mu[3,2]"), lags = 10)

mcmc_acf(posterior22, par = c("v_sigma[1,1]", "v_sigma[1,2]", "v_sigma[2,1]", "v_sigma[2,2]", "v_sigma[3,1]", "v_mu[3,2]"), lags = 10)

情報量基準(WAICとLOOCV)の算出

@berobero11氏のStan Advent 2016を参考にしつつWAICとLOOCVの算出をしてみました。どちらの指標においても，条件の効果を考慮したモデルが当てはまりがう良さそうです。

# WAICとLOOCVの計算
extractFit1 <- rstan::extract(fit1)
extractFit2 <- rstan::extract(fit2)
fitValue <- NULL
fitValue[1] <- loo::loo(extractFit1$log_lik)$looic/(2*numberAllTrials)
fitValue[2] <- loo::waic(extractFit1$log_lik)$waic/(2*numberAllTrials)
fitValue[3] <- loo::loo(extractFit2$log_lik)$looic/(2*numberAllTrials)
fitValue[4] <- loo::waic(extractFit2$log_lik)$waic/(2*numberAllTrials)

# プロットの準備
modelNames <- c("Without Condition Effect","Without Condition Effect","With Condition Effect","With Condition Effect")
valueNames <- c("LOOCV","WAIC","LOOCV","WAIC")
WAICData <- data.frame(modelNames,valueNames,fitValue)

# Bar plotの作成
bp <- ggplot(WAICData,aes(x=modelNames,y=fitValue,fill=valueNames)) + geom_bar(position="dodge",stat="identity")
bp <- bp + xlab("Model")+ylab("Values")+geom_text(aes(label = sprintf("%5.4f",fitValue)), vjust = 1.5, colour = "black", position = position_dodge(0.9),size = 5)+ylim(-1, NA)
plot(bp)

推定結果とパラメータリカバリのチェック

以降では，推定によって，元々設定したパラメータの値をリカバリできているのかチェックをしていきます（推定値の点推定値としてはMAPとEAPを算出しましたが，プロットはEAPのみにしています）。

パラメータA

まあまあ，パラメターリカバリできているかな。

# MAP推定値計算関数を準備
mapEstimate <- function(z){
  density(z)$x[which.max(density(z)$y)]
}

fit1AEap <- apply(extractFit1$A,2,mean)
fit1AMap <- apply(extractFit1$A,2,mapEstimate)
fit2AEap <- apply(extractFit2$A,2,mean)
fit2AMap <- apply(extractFit2$A,2,mapEstimate)
parameterA <-data.frame(parameters$A,fit1AEap,fit1AMap,fit2AEap,fit2AMap)
print(parameterA)

##    parameters.A  fit1AEap   fit1AMap  fit2AEap  fit2AMap
## 1     0.2197622 0.3596942 0.37050078 0.3885457 0.4046331
## 2     0.3849113 0.1382490 0.05222353 0.2153580 0.1267809
## 3     1.2793542 1.2193958 1.24630463 1.3172280 1.3182410
## 4     0.5352542 0.3197556 0.33902718 0.5531708 0.6318820
## 5     0.5646439 0.5233939 0.53986090 0.6136198 0.5923683
## 6     1.3575325 0.8050358 0.89361177 0.9848496 1.0369273
## 7     0.7304581 0.6254361 0.68121676 0.8972380 0.9707640
## 8     1.1120409 1.2784274 1.28359037 1.2380688 1.2450074
## 9     0.1565736 0.4323484 0.46767085 0.4532021 0.4953147
## 10    0.2771690 0.2904164 0.23710658 0.4131040 0.4468086

Apfit1Eap <- ggplot(parameterA, aes(x=parameters.A, y=fit1AEap))+geom_point()+ xlab("Estimated parameter in Without Condition Effect Model")+ylab("Original value of A")
Apfit2Eap <- ggplot(parameterA, aes(x=parameters.A, y=fit2AEap))+geom_point()+ xlab("Estimated parameter in With Condition Effect Model")+ylab("Original value of A")
plot(Apfit1Eap)

plot(Apfit2Eap)

パラメータk

まあまあ，パラメターリカバリできているかな。

# 元々はbからデータを作っていたので，kに変更(k=b-A)
originK <- parameters$b-parameters$A

fit1kEap <- apply(extractFit1$k,2,mean)
fit1kMap <- apply(extractFit1$k,2,mapEstimate)
fit2kEap <- apply(extractFit2$k,2,mean)
fit2kMap <- apply(extractFit2$k,2,mapEstimate)
parameterK <-data.frame(originK,fit1kEap,fit1kMap,fit2kEap,fit2kMap)
print(parameterK)

##      originK  fit1kEap  fit1kMap  fit2kEap  fit2kMap
## 1  0.9601447 0.7052006 0.6612541 0.9209895 0.9375435
## 2  0.8154745 0.4124235 0.4316756 0.7645841 0.7987263
## 3  0.3665154 0.4254030 0.3588894 0.5682246 0.5111807
## 4  0.5200872 0.6047961 0.6145263 0.6440033 0.5416150
## 5  0.1574356 0.1202811 0.1074498 0.1783537 0.1570896
## 6  0.5231592 0.4858275 0.3678424 0.7171649 0.6095679
## 7  1.1629985 0.8687153 0.7451857 1.0639499 0.9486935
## 8  0.6027215 0.5182855 0.4772118 0.8695459 0.8037778
## 9  1.0923517 0.7831610 0.6994084 1.1512616 1.1029487
## 10 1.0735089 0.4516198 0.4804440 0.6854127 0.6061864

kpfit1Eap <- ggplot(parameterK, aes(x=originK, y=fit1kEap))+geom_point()+ xlab("Estimated parameter in Without Condition Effect Model")+ylab("Original value of k")
kpfit2Eap <- ggplot(parameterK, aes(x=originK, y=fit2kEap))+geom_point()+ xlab("Estimated parameter in With Condition Effect Model")+ylab("Original value of k")
plot(kpfit1Eap)

plot(kpfit2Eap)

パラメータTau

まあまあ，パラメターリカバリできているかな。

fit1psiEap <- apply(extractFit1$psi,2,mean)
fit1psiMap <- apply(extractFit1$psi,2,mapEstimate)
fit2psiEap <- apply(extractFit2$psi,2,mean)
fit2psiMap <- apply(extractFit2$psi,2,mapEstimate)
parameterPsi <-data.frame(parameters$t0,fit1psiEap,fit1psiMap,fit2psiEap,fit2psiMap)
print(parameterPsi)

##    parameters.t0 fit1psiEap fit1psiMap fit2psiEap fit2psiMap
## 1      0.2330341  0.2542777  0.2549910 0.23602476 0.23320105
## 2      0.4659625  0.5082014  0.5059192 0.47775354 0.47348916
## 3      0.2659726  0.2417448  0.2575985 0.22742514 0.23692016
## 4      0.8578277  0.8248316  0.8249007 0.83384639 0.84220059
## 5      0.8100644  0.8168498  0.8217065 0.80906909 0.81227905
## 6      0.4422001  0.4138236  0.4476275 0.39606031 0.41790517
## 7      0.7989248  0.8206257  0.8289087 0.81669371 0.83232547
## 8      0.1218993  0.1135934  0.1276360 0.07814652 0.08905505
## 9      0.5609480  0.5635367  0.5681844 0.54558210 0.55010207
## 10     0.2065314  0.2752165  0.2620464 0.26701436 0.27176891

Psipfit1Eap <- ggplot(parameterPsi, aes(x=parameters.t0, y=fit1psiEap))+geom_point()+ xlab("Estimated parameter in Without Condition Effect Model")+ylab("Original value of Tau")
Psipfit2Eap <- ggplot(parameterPsi, aes(x=parameters.t0, y=fit2psiEap))+geom_point()+ xlab("Estimated parameter in With Condition Effect Model")+ylab("Original value of Tau")
plot(Psipfit1Eap)

plot(Psipfit2Eap)

パラメータV(条件１の選択肢１)

当たり前だが，With Condition Effect Modelの方は，パラメータリカバリできている。

fit1v11Eap <- apply(extractFit1$v[,,1],2,mean)
fit1v11Map <- apply(extractFit1$v[,,1],2,mapEstimate)
fit2v11Eap <- apply(extractFit2$v[,,1,1],2,mean)
fit2v11Map <- apply(extractFit2$v[,,1,1],2,mapEstimate)
parameterV11 <-data.frame(parameters$v1_c1,fit1v11Eap,fit1v11Map,fit2v11Eap,fit2v11Map)
print(parameterV11)

##    parameters.v1_c1 fit1v11Eap fit1v11Map fit2v11Eap fit2v11Map
## 1         2.8215811   2.462893   2.472908   2.976195   3.012800
## 2         2.6886403   2.302890   2.272929   2.432274   2.437962
## 3         2.5539177   2.187280   2.240267   2.633322   2.561235
## 4         1.9380883   2.336235   2.293017   2.193843   2.180408
## 5         1.6940373   2.621454   2.441395   1.633005   1.582418
## 6         1.6195290   1.930384   1.940349   1.618144   1.607161
## 7         1.3052930   2.329237   2.306559   1.475796   1.406477
## 8         1.7920827   2.591646   2.554154   2.127964   2.104065
## 9         0.7346036   1.924409   1.934199   1.349092   1.358975
## 10        4.1689560   2.022118   2.060492   3.706664   3.624642

V11pfit1Eap <- ggplot(parameterV11, aes(x=parameters.v1_c1, y=fit1v11Eap))+geom_point()+ xlab("Estimated parameter in Without Condition Effect Model")+ylab("Original value of v in condition 1 & choice1")
V11pfit2Eap <- ggplot(parameterV11, aes(x=parameters.v1_c1, y=fit2v11Eap))+geom_point()+ xlab("Estimated parameter in With Condition Effect Model")+ylab("Original value of v in condition 1 & choice1")
plot(V11pfit1Eap)

plot(V11pfit2Eap)