I'm currently trying to fit a shifted inverse gaussian to reaction time data (always postive). My paramterization of the model includes 3 parameters, alpha, gamma and tau, which must always be positive also. The model expression is appropriate and comes from this paper: The shifted Wald distribution for response time data analysis

My problem is, that I expect my effects to be additive, and I'm not sure how I can constrain my parameters (with sensible prior choices?) to always be positive in this case. Of course, a log-transform works quite well, but leads to a different interpretation of the data which I am less interested in because the effects of my predictors are multiplicative in this case.

I have two questions:

1) Is it a major issue for Bayesian hierarchical modeling to allow for nonsensical values? If I run my unconstrained model with appropriate starting values the diagnostics are good, and estimated values make sense. But is it interpretable as is? I assumed the sampling process would just avoid nonsensical regions of parameter space and that it didn't matter, but I'm probably wrong?

2) If it is a major issue, how can I constrain my parameters to be on the positive half-line, in an additive model that makes sense?

I'm using Stan, the code is pasted below. I'm adding 0.150s to theta (the minimum time before response) because it's its lower bound due to the way data was collected and preprocessed.

Thank you!

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)))));
    tmp1 = log(alpha) - 0.5 * (log(2) + log(pi()) + 3 * log(x - theta));
    tmp2 = -1 * ((alpha - gamma * (x-theta))^2/(2*(x-theta)));
    return tmp1+tmp2 ;



 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


parameters {
  // fixed-effects parameters
  vector[12] b;                   // fixed effect matrix
  vector<lower=0>[4] sigma_u_g;             // random effects standard deviations
  vector<lower=0>[4] sigma_u_a;
  vector<lower=0>[4] sigma_u_t;
  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;

transformed parameters {
  matrix[4,J] u_g;
  matrix[4,J] u_a;
  matrix[4,J] u_t;
  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; 

// model parameters
  real gamma;                // drift
  real alpha;                // boundary
  real theta;                // ndt

    b ~ normal(0, 100); // diffuse prior
  L_u_g ~ lkj_corr_cholesky(2);   // LKJ prior for the correlation matrix
  L_u_a ~ lkj_corr_cholesky(2);   
  L_u_t ~ lkj_corr_cholesky(2);   
  to_vector(z_u_g) ~ normal(0,1); // random effects are (initially)     normal variates with SD=1
  to_vector(z_u_a) ~ normal(0,1); 
  to_vector(z_u_t) ~ normal(0,1); 
  sigma_u_g ~ cauchy(0, 2);       // SD of random effects (vectorized)
  sigma_u_a ~ cauchy(0, 2);       
  sigma_u_t ~ cauchy(0, 2);       

  for(i in 1 : N_obs) {

  gamma = b[1] + u_g[1,id[i]] + (b[2] + u_g[2,id[i]])*cg[i] + (b[3] + u_g[3,id[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 + 0.150);

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.