0

私は JAGS で階層ドリフト拡散モデルを実行しています。このモデルは、 Pleskac et al., 2018から採用しました。この実験は、2 択の反応時間タスクでした。

私の目標は、3 つの実験条件について個別のパラメーター推定値 (muAlpha、muBeta、muDelta、および muNDT) を取得することです。さらに、muAlpha と muBeta の推定値は con1 (factor1) の被験者操作に依存し、muDelta と muNDT は con2 (factor1*factor2) の被験者操作に依存します。

私の問題は、最初の実験条件の事後推定 (muAlpha、muBeta、muDelta、および muNDT) のみがもっともらしい形状に見え、2 番目と 3 番目の実験条件の事後推定は一様分布のように見えることです。

同様の問題に関する別の質問を見つけましたが、私が理解している限り、モデル内のネストされたインデックス作成に誤りがありました。残念ながら、私は自分のモデルで何を誤って指定したかを理解できませんでした (以下を参照)。誰​​かが私を助けてくれれば幸いです!

#For demonstration I generated a data set of 30 participants (10 in each condition) with 100 trials each:

library(rjags)
library(runjags)
library(tidyr)    
library(retimes) #for generating reaction time data
load.module("wiener") #for using the wiener distribution in JAGS (see info below model)

#Generate an ex-Gaussian population for RT-values:
rt <- rexgauss(3000, 300, 10, 20, positive = T)

#Generate responses
x = rep(0:1, length(rt))
resp <- sample(x, length(rt))

#Generate factor1 and factor2 within-subject conditions
y = rep(1:2, length(rt))
f1 <- sample(y, length(rt))
f2 <- sample(y, length(rt))

#Generate condition & subject ID
cond <- rep(1:3, each = 100, times = 10)
sub <- rep(1:30, each = 100)

#bind in data frame
df <- cbind(sub, cond, f1, f2, resp, rt)
df <- as.data.frame(df)

#Change "0" response times to negative for JAGS wiener module
df$rt <- ifelse(df$rt < 100, NA, df$rt)
y <- ifelse(df$resp == 0, (df$rt*-1)/1000, df$rt/1000)
y <- ifelse(abs(y) < .1, NA, y) #Remove data below 100ms (0)
y <- ifelse(abs(y) > 3, NA, y) #Remove data above 3000ms (0)

#alpha and beta vary by f2
con1 <- as.numeric(df$f2)

#tau and delta vary by f1 and f2
con2 <- as.numeric(with(df, interaction(f1, f2)))

#indicator vector of between subject condition 
btwnCon <- as.numeric(df$cond)

#Create numbers for looping purposes
nData <- nrow(df)
nSub <- length(unique(df$sub))
ncon1 <- max(con1)
ncon2 <- max(con2)
nBtwn <- max(btwnCon)

#Create a list of the data for JAGS
datalist <- list(y = y, con1 = con1, con2 = con2, sub = df$sub, nData = nData, 
                 nSub = nSub, ncon1 = ncon1, ncon2 = ncon2, btwnCon = btwnCon, nBtwn = nBtwn)

#These are the model specifications:

#initial values
initfunction <- function(chain){
  return(list(
    muAlpha = matrix(runif(ncon1,.2,4.9),nrow=ncon1,ncol=nBtwn,byrow=TRUE), #threshold mean
    muBeta = matrix(runif(ncon1, .15, .85),nrow=ncon1,ncol=nBtwn,byrow=TRUE), #start point mean
    muNDT = matrix(runif(ncon2, .01, .05),nrow=ncon2,ncol=nBtwn,byrow=TRUE), #non-decision time mean
    muDelta = matrix(runif(ncon2, -4.9, 4.9),nrow=ncon2,ncol=nBtwn,byrow=TRUE), #drift rate mean
    tauAlpha = runif(nBtwn, .01, 100), #threshold precision 
    tauBeta = runif(nBtwn, .01, 100), #start point precision
    tauNDT = runif(nBtwn, .01, 100), #non-decision precision
    tauDelta = runif(nBtwn, .01, 100), #drift rate precision
    .RNG.name = "lecuyer::RngStream",
    .RNG.seed = sample.int(1e10, 1, replace = F)))
}

#And this is the model:

model {
#hyperpriors
for (bIdx1 in 1:nBtwn) { #between condition index
    for (bCon2 in 1:ncon2) { #within conditions index con2
        muDelta[bCon2,bIdx1] ~ dunif(-5,5)
        muNDT[bCon2,bIdx1] ~ dunif(.001,1)
        #sample hyper priors for inference tests
        muDeltaPrior[bCon2,bIdx1] ~ dunif(-5,5)
        muNDTPrior[bCon2,bIdx1] ~ dunif(.001,1)
    }
}
for (bIdx2 in 1:nBtwn) { #between condition index
    for (bCon1 in 1:ncon1) { #within condition index con1
        muBeta[bCon1,bIdx2] ~ dunif(.1,.9)
        muBetaPrior[bCon1,bIdx2] ~ dunif(.1,.9)
        muAlpha[bCon1,bIdx2] ~ dunif(.1,5)
        muAlphaPrior[bCon1,bIdx2] ~ dunif(.1,5)
    }
}
for (bIdx3 in 1:nBtwn) { #between condition index
    tauDelta[bIdx3] ~ dunif(0,1000)
    tauBeta[bIdx3] ~ dunif(0,1000)
    tauAlpha[bIdx3] ~ dunif(0,1000)
    tauNDT[bIdx3] ~ dunif(0,1000)
}

#PRIORS
for (sIdx in 1:nSub) { #subject index
    for (sCon1 in 1:ncon1) { #within subject index con1
        beta[sCon1,sIdx] ~ dnorm(muBeta[sCon1, btwnCon[sIdx]], tauBeta[btwnCon[sIdx]]) 
        alpha[sCon1,sIdx] ~ dnorm(muAlpha[sCon1, btwnCon[sIdx]], tauAlpha[btwnCon[sIdx]]) 
    }
    for (sCon2 in 1:ncon2) { #within subject index con2
        delta[sCon2,sIdx] ~ dnorm(muDelta[sCon2, btwnCon[sIdx]], tauDelta[btwnCon[sIdx]])
    ndt[sCon2,sIdx] ~ dnorm(muNDT[sCon2, btwnCon[sIdx]], tauNDT[btwnCon[sIdx]])
}}

for (dIdx in 1:nData) {
    y[dIdx] ~ dwiener(alpha[con1[dIdx], sub[dIdx]], ndt[con2[dIdx], sub[dIdx]], beta[con1[dIdx], sub[dIdx]], delta[con2[dIdx], sub[dIdx]])
}
}

(モデルを実行する場合は、JAGS のJAGS Wiener モジュール拡張機能をインストールする必要があります。)

4

0 に答える 0