# Multilevel modeling, constraint for positive values

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 ;
}

}

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

}

parameters {
// fixed-effects parameters
vector b;                   // fixed effect matrix
vector<lower=0> sigma_u_g;             // random effects standard deviations
vector<lower=0> sigma_u_a;
vector<lower=0> sigma_u_t;
cholesky_factor_corr L_u_g;
cholesky_factor_corr L_u_a;
cholesky_factor_corr 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{
// model parameters
real gamma;                // drift
real alpha;                // boundary
real theta;                // ndt

//priors
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 + u_g[1,id[i]] + (b + u_g[2,id[i]])*cg[i] + (b + u_g[3,id[i]] + (b + u_g[4,id[i]])*cg[i])*soa[i];
alpha = b + u_a[1,id[i]] + (b + u_a[2,id[i]])*cg[i] + (b + u_a[3,id[i]] + (b + u_a[4,id[i]])*cg[i])*soa[i];
theta = b + u_t[1,id[i]] + (b + u_t[2,id[i]])*cg[i] + (b + u_t[3,id[i]] + (b + u_t[4,id[i]])*cg[i])*soa[i];

target += shifted_Wald_lpdf(rt_obs[i] | gamma, alpha, theta + 0.150);
}
}