Dear Stan community,

I am currently fitting a shifted Wald (or shifted inverse gaussian) to reaction time data from an experiment

using rstan. I already asked a question about the same model, but it’s very different from this one, so I figured I should start a new thread, I hope that’s the way I should have proceeded.

I would like to evaluate the effect of experimental conditions A, B and their interaction on the parameters of this distribution.

The model converges nicely if I’m only using fixed-effects, but given the large inter-individual variability among my participants, I would really like to add random effects for all (or some of) my predictors.

Another feature I wanted to implement, is an across-trial variation of the parameter gamma of my distribution (varying drift as in a drift diffusion process).

When I try running the model (below) with everything at once, I can’t get it to start, with the following output:

Chain 1: Rejecting initial value:

Chain 1: Log probability evaluates to log(0), i.e. negative infinity.

Chain 1: Initialization between (-2, 2) failed after 100 attempts.

Chain 1: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the >model.

[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”

[1] “error occurred during calling the sampler; sampling not done”

I can get the model to start after several tries, if I only add random effects by participant (no u_drift) and set the range of initial values to 10e-15, which seems really odd.

My dataset is of about 4400 datapoints split across 20 participants, two subconditions A and two subconditions B (so an average of 55 points per subcondition per participant).

I’ve pieced this model together using tutorials, but I’m probably doing something wrong.

I would be grateful for suggestions or tips pointing me in the direction of what is causing this error, and how to solve it.

Thank you very much!

```
functions {
real shifted_Wald_lpdf(real x, real gamma, real alpha, real theta){
real tmp1;
real tmp2;
tmp1 = log(alpha / (sqrt(2 * pi() * (pow((x - theta), 3)))));
tmp2 = -1 * ((alpha - gamma * (x-theta))^2/(2*(x-theta)));
return tmp1+tmp2 ;
}
}
data{
int<lower=0> N_obs; // observed rts
int<lower=1> J; // n of subjects
vector<lower=0>[N_obs] rt_obs;
int<lower=1, upper=2> soa[N_obs]; // condition A
int<lower=0,upper=1> cg[N_obs]; // condition B
int<lower=1,upper=J> id[N_obs]; // subj identifier
int<lower=0> Tnum[N_obs]; // trial number
}
parameters {
// fixed-effects parameters
vector[12] b;
// random effects parameters
vector<lower=0>[4] sigma_u_g; // random effects for gamma (drift)
vector<lower=0>[4] sigma_u_a; // random effects for alpha (boundary)
vector<lower=0>[4] sigma_u_t; // random effects for theta (non-decision time)
cholesky_factor_corr[4] L_u_g;
cholesky_factor_corr[4] L_u_a;
cholesky_factor_corr[4] L_u_t;
matrix[4,J] z_u_g; // random effect matrix
matrix[4,J] z_u_a;
matrix[4,J] z_u_t;
vector<lower=0>[2] sigma_drift; // for variable drift rates across trials
cholesky_factor_corr[2] L_drift;
matrix[2,N_obs] z_drift;
}
transformed parameters {
matrix[4,J] u_g;
matrix[4,J] u_a;
matrix[4,J] u_t;
matrix[2, N_obs] u_drift;
u_g = diag_pre_multiply(sigma_u_g, L_u_g) * z_u_g;
u_a = diag_pre_multiply(sigma_u_a, L_u_a) * z_u_a;
u_t = diag_pre_multiply(sigma_u_t, L_u_t) * z_u_t;
u_drift = diag_pre_multiply(sigma_drift, L_drift) * z_drift;
}
model{
// model parameters
real gamma; // drift
real alpha; // boundary
real theta; // ndt
//priors
L_u_g ~ lkj_corr_cholesky(2);
L_u_a ~ lkj_corr_cholesky(2);
L_u_t ~ lkj_corr_cholesky(2);
to_vector(z_u_g) ~ normal(0,1);
to_vector(z_u_a) ~ normal(0,1);
to_vector(z_u_t) ~ normal(0,1);
sigma_u_g ~ cauchy(0, 1);
sigma_u_a ~ cauchy(0, 1);
sigma_u_t ~ cauchy(0, 1);
sigma_drift ~ cauchy(0,1);
L_drift ~ lkj_corr_cholesky(2);
to_vector(z_drift) ~ normal(0,1);
b ~ normal(0, 100); // diffuse prior for beta coefficients
for(i in 1 : N_obs) {
gamma = b[1] + u_g[1,id[i]] + u_drift[1,Tnum[i]] + (b[2] + u_g[2,id[i]])*cg[i] + (b[3] + u_g[3,id[i]] + u_drift[2,Tnum[i]] + (b[4] + u_g[4,id[i]])*cg[i])*soa[i];
alpha = b[5] + u_a[1,id[i]] + (b[6] + u_a[2,id[i]])*cg[i] + (b[7] + u_a[3,id[i]] + (b[8] + u_a[4,id[i]])*cg[i])*soa[i];
theta = b[9] + u_t[1,id[i]] + (b[10] + u_t[2,id[i]])*cg[i] + (b[11] + u_t[3,id[i]] + (b[12] + u_t[4,id[i]])*cg[i])*soa[i];
target += shifted_Wald_lpdf(rt_obs[i]|gamma,alpha,theta);
}
}
```