Dear Stan community,
I’m new and currently trying to fit a right-censored shifted-Wald distribution to reaction time data using rstan. Using the Wiener distribution makes less sense due to the way the experiment was designed (it’s a go/no go task, so there is only one decision boundary and not two). Mainly, I would like to evaluate the effect of experimental conditions A, B and their interaction on the parameters of this distribution.
I managed to produce the following code with the help of the available literature. It may not be optimal, but the chains converge so far.
It seems like I should use the log pdf of the shifted-Wald, can someone explain whether this is the appropriate choice and why?
I would like to check the quality of the model fit by means of posterior predictive checking, but I’m at a loss as to what I should do within the generated quantities block given my model. Would you mind pointing me to the appropriate documentation if I missed it, or providing me with some help?
Thank you in advance.
DRR
functions {
// shifted Wald pdf (see Miller 2018)
real shifted_Wald_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); // unsure about why I should use the log.
}
real swald_lccdf(real x, real gamma, real alpha, real theta){
real tmp1;
real tmp2;
real tmp3;
// compute the log of the complementary CDF of the shifted wald (see Miller 2018)
tmp1 = Phi((gamma*(x - theta) - alpha)/sqrt(x-theta));
tmp2 = exp(2*alpha*gamma);
tmp3 = Phi((gamma*(x - theta) + alpha)/sqrt(x-theta));
return log(1 - tmp1-(tmp2*tmp3));
}
}
data{
int<lower=0> N_obs; // number of observed RTs
int<lower=1> J; // n of participants
vector<lower=0>[J] N_cens; // number of right-censored trials
vector<lower=0>[N_obs] rt_obs; // RT data
vector<lower=0>[J] U; // censoring point (eventually this will be fixed for all participants when the actual censoring point is retrieved)
int<lower=1, upper=2> soa[N_obs]; // experimental manipulation A (2 levels)
int<lower=0,upper=1> cg[N_obs]; // experimental manipulation B (2 levels)
int<lower=1,upper=J> id[N_obs]; // identifier
}
parameters {
// fixed-effects parameters
vector[12] b;
}
transformed parameters {
}
model{
// parameters
real gamma; // drift
real alpha; // boundary
real theta; // ndt
//priors
b ~ normal(0, 100); // diffuse prior
for(i in 1 : N_obs) {
gamma = b[1] + (b[2])*cg[i] + (b[3] + (b[4])*cg[i])*soa[i];
alpha = b[5] + (b[6])*cg[i] + (b[7] + (b[8])*cg[i])*soa[i];
theta = b[9] + (b[10])*cg[i] + (b[11] + (b[12])*cg[i])*soa[i];
target += shifted_Wald_lpdf(rt_obs[i]|gamma,alpha,theta) + N_cens[id[i]] * swald_lccdf(U[id[i]]|gamma,alpha,theta);
}
}
generated quantities {
vector[N_obs] log_lik;
real gamma;
real alpha;
real theta;
for (i in 1:N_obs){
gamma = b[1] + (b[2])*cg[i] + (b[3] + (b[4])*cg[i])*soa[i];
alpha = b[5] + (b[6])*cg[i] + (b[7] + (b[8])*cg[i])*soa[i];
theta = b[9] + (b[10])*cg[i] + (b[11] + (b[12])*cg[i])*soa[i];
log_lik[i] = shifted_Wald_lpdf(rt_obs[i]|gamma,alpha,theta) + N_cens[id[i]] * swald_lccdf(U[id[i]]|gamma,alpha,theta);
}
}