Hello Stan community,
I am trying to fit a drift-diffusion model (DDM) on some reaction time (RT) data. Participants had to choose between a 20 cents (lower boundary) and 60 cents (upper boundary) bet.
During data preprocessing, I removed all trials faster than 150 ms (a frequently used lower-bound for RTs in the literature).
I used a hierarchical model with the basic Wiener function:
data {
int <lower=1> nSubs;
int <lower=1> maxTrials;
int <lower=1> subTrials[nSubs];
real <lower=-9> betsize[nSubs, maxTrials];
real <lower=-9> betRT[nSubs, maxTrials];
}
parameters {
// Group level
//// Non-decision time
real <lower=0.01, upper=2> mu_tau;
real <lower=0.01, upper=2> sd_tau;
//// Boundary separation
real <lower=0.1> mu_alpha;
real <lower=0.01, upper=5> sd_alpha;
//// Response bias
real <lower=-4, upper=4> mu_beta;
real <lower=0> sd_beta;
//// Drift rate
real mu_delta;
real <lower=0> sd_delta;
// Subject level
//// Non-decision time
vector <lower=0.01, upper=2> [nSubs] tau;
//// Boundary separation
vector <lower=0.1> [nSubs] alpha;
//// Response bias
vector <lower=-4, upper=4> [nSubs] beta;
//// Drift rate
vector [nSubs] delta;
}
model {
// Define values
real betaval[nSubs];
// Sample parameter values
//// Group level
mu_tau ~ uniform(0.01, 2);
sd_tau ~ uniform(0.01, 2);
mu_alpha ~ uniform(0.1, 5);
sd_alpha ~ uniform(0.01, 5);
mu_beta ~ uniform(-4, 4);
sd_beta ~ uniform(0, 1);
mu_delta ~ normal(0, 10);
sd_delta ~ uniform(0, 10);
//// Subject level
for (subi in 1:nSubs) {
tau[subi] ~ normal(mu_tau, sd_tau)T[0.01, ];
alpha[subi] ~ normal(mu_alpha, sd_alpha)T[0.1, ];
beta[subi] ~ normal(mu_beta, sd_beta);
delta[subi] ~ normal(mu_delta, sd_delta);
}
for (subi in 1:nSubs) {
// Parameter values
betaval[subi] = Phi(beta[subi]);
for (triali in 1:subTrials[subi]) {
if (betsize[subi, triali] > 0) {
// Sample RT
if (betsize[subi, triali] == 60) {
betRT[subi, triali] ~ wiener(
alpha[subi], tau[subi], betaval[subi], delta[subi]
);
} else {
betRT[subi, triali] ~ wiener(
alpha[subi], tau[subi], 1 - betaval[subi], -delta[subi]
);
}
}
}
}
}
generated quantities {
// Define values
real betaval[nSubs];
real log_lik[nSubs];
for (subi in 1:nSubs) {
betaval[subi] = Phi(beta[subi]);
log_lik[subi] = 0;
for (triali in 1:subTrials[subi]) {
if (betsize[subi, triali] > 0) {
// Log likelihood
if (betsize[subi, triali] == 60) {
log_lik[subi] = log_lik[subi] + wiener_lpdf(
betRT[subi, triali] | alpha[subi], tau[subi], betaval[subi],
delta[subi]
);
} else {
log_lik[subi] = log_lik[subi] + wiener_lpdf(
betRT[subi, triali] | alpha[subi], tau[subi], 1 - betaval[subi],
-delta[subi]
);
}
}
}
}
}
When I fit this model and do some posterior checks, I find that several participants have RTs faster than the mean estimated non-decision time (tau
) at the subject level. This does not make sense. It should not be possible to make a decision faster than the non-decision time.
I have tried many different boundaries for the priors of tau
(e.g., lower=0.001, upper=1
). The only time the model estimated tau
< min(RT) is when I kept removed RTs < 250 ms. However, this is not great because then I am removing a lot of data and I have no theoretical justification for a cutoff of 250 ms.
I have several questions about this and I hope you can help me.
- Why does this happen?
- Why is Stan not “Rejecting initial value”?
- It should be unable to estimate an RT faster than
tau
- It should be unable to estimate an RT faster than
- How can I fix this?
Any help is appreciated.
Thank you in advance.