0

我在 JAGS 中运行分层漂移扩散模型，我改编自：Pleskac et al., 2018该实验是一个二选一的反应时间任务。

我的目标是为 3 个实验条件获得单独的参数估计值（muAlpha、muBeta、muDelta 和 muNDT）。此外，估计 muAlpha 和 muBeta 依赖于对 con1 (factor1) 的主体操作，而 muDelta 和 muNDT 取决于对 con2 (factor1*factor2) 的主体操纵。

我的问题是，只有第一个实验条件的后验估计（muAlpha、muBeta、muDelta 和 muNDT）在形状上看起来是合理的，而第二个和第三个实验条件的后验估计看起来像是均匀分布。

我发现另一个类似问题的问题，据我了解，模型中的嵌套索引存在错误。不幸的是，我无法弄清楚我在模型中错误指定的内容（见下文），如果有人能帮助我，我将不胜感激！

#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