Dear STAN community,
I am using a STAN code for censored hieararchical shifted Wald model provided by Faulkenberry(Faulkenberry, T. (2023, June). A hierarchical Bayesian shifted Wald model with censoring. Paper presented at Virtual MathPsych/ICCM 2023. Via mathpsych.org/presentation/1283)
I am encountering several issues, mostly related to the estimation of non-decision time (theta; H). Here is the original code:
functions{
// shifted Wald (log) PDF
real shiftedWald_lpdf(real X, real gamma, real alpha, real theta){
real tmp1;
real tmp2;
tmp1 = alpha / (sqrt(2 * pi() * (pow((X - theta), 3))));
tmp2 = exp(-1 * (pow((alpha - gamma * (X-theta)),2)/(2*(X-theta))));
return log(tmp1*tmp2) ;
}
// shifted Wald (log) survivor function
real shiftedWald_lS(real X, real gamma, real alpha, real theta){
real tmp1;
real tmp2;
real tmp3;
// compute log of the survival function 1-F for swald
tmp1 = Phi((gamma*(X-theta) - alpha)/sqrt(X-theta));
tmp2 = exp(2*alpha*gamma);
tmp3 = Phi(-1*(gamma*(X-theta) + alpha)/sqrt(X-theta));
return log(1-tmp1-tmp2*tmp3);
}
}
data{
int ns; // number of subjects
int nrt; // number of trials per subject
real rt[ns,nrt]; // matrix of response times
real D[ns,nrt]; // censoring matrix (0 = correct RT, 1 = incorrect RT)
}
parameters{
// noncentered parameters
vector[ns] G_raw;
vector[ns] A_raw;
vector[ns] H_raw;
// define uniform hyperpriors from Matzke & Wagenmakers (2009)
real<lower=0.85,upper=7.43> g;
real<lower=0.67,upper=2.35> a;
real<lower=0,upper=0.82> h;
real<lower=0,upper=1.899> gS;
real<lower=0,upper=0.485> aS;
real<lower=0,upper=0.237> hS;
}
transformed parameters{
// noncentered parameterization (Matt trick) (Betancourt & Girolami, 2003)
vector[ns] G;
vector[ns] A;
vector[ns] H;
for(i in 1:ns){
G[i] = g + gS*G_raw[i];
A[i] = a + aS*A_raw[i];
H[i] = h + hS*H_raw[i];
}
}
model{
for(i in 1:ns){
// group level prior
G_raw[i] ~ normal(0,1);
A_raw[i] ~ normal(0,1);
H_raw[i] ~ normal(0,1);
for(j in 1:nrt){
target += (1-D[i,j])*shiftedWald_lpdf(rt[i,j]| G[i], A[i], H[i]) + D[i,j]*shiftedWald_lS(rt[i,j], G[i], A[i], H[i]);
}
}
}
I encountered many divergent samples due to proposed non-decision time values exceeding observed reaction time. In reaction to that, I introduced to the data block a vector of minimum reaction times for each participant
data {
.
.
.
vector [ns] min_rt;
}
Additionally, to ensure non-decision time is smaller than the minimum reaction time of the given person, I somehow clumsily constrained this parameter in the following way ( which I know is not ideal):
transformed parameters {
.
.
for (i in 1:ns) {
.
.
H[i] = fmin(h + hS * H_raw[i], min_rt[i]-.01);
}
}
However, in several instances, the non-decision time (theta; H) was negative, but it has to be strictly positive.
Could you please provide me with a solution on how to ensure positivity of non-decision time, and at the same time, to make it smaller than the minimum reaction time for each participant?
Thank you in advance.
Tom