Estimated non-decision times larger than reaction times with wiener function in RStan

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
  • How can I fix this?

Any help is appreciated.
Thank you in advance.