MCMC chain did not coverge in fitting ITCH model

Modeling
,
6

To constrain the beta to be either positive or negative when using the non-centered parameterisation, you need to place the appropriate lower= or upper= on the respective _raw variable.

For your beta variables, if you want to constrain them to be positive, then you need to declare the respective beta_raw variable with the lower bound: <lower = -mean_raw/sd_raw> (more background on that derivation in this post). If you want them to be negative, then that constraint simply becomes <upper= -mean_raw/sd_raw>.

Based on your original code, I’ll assume that beta[1] (the intercept) can be either positive or negative, beta[2] and beta[3] should be positive, beta[4] should be negative, and beta[5] should be positive. But if that’s not the case, then you can just modify the code accordingly.

Here’s how that construction would look:

data{
  int<lower=1> n_s; // number of subjects
  int<lower=1> nTrials; // number of trials
  int<lower=0, upper=2> choice[n_s, nTrials]; // 0 for ss choice, 1 for ll choice
  vector<lower=0>[nTrials] ss_money[n_s];
  vector<lower=0>[nTrials] ss_time[n_s];
  vector<lower=1>[nTrials] ll_money[n_s];
  vector<lower=1>[nTrials] ll_time[n_s];
 
} // end data 

parameters{
// group raw parameters
  // mean
  vector[5] mean_raw;  
  // SD
  vector <lower=0> [5] sd_raw;
// subject raw parameters
  vector[n_s] beta_raw1;
  vector<lower = -mean_raw[2]/sd_raw[2]>[n_s] beta_raw2; // implies strictly positive beta[2]
  vector<lower = -mean_raw[3]/sd_raw[3]>[n_s] beta_raw3; // implies strictly positive beta[3]
  vector<upper = -mean_raw[4]/sd_raw[4]>[n_s] beta_raw4; // implies strictly negative beta[4]
  vector<lower = -mean_raw[5]/sd_raw[5]>[n_s] beta_raw5; // implies strictly positive beta[5]
}

transformed parameters{
  // real subject parameters
  vector[n_s] beta[5];

  beta[1] = mean_raw[1] + beta_raw1 * sd_raw[1];
  beta[2] = mean_raw[2] + beta_raw2 * sd_raw[2];
  beta[3] = mean_raw[3] + beta_raw3 * sd_raw[3];
  beta[4] = mean_raw[4] + beta_raw4 * sd_raw[4];
  beta[5] = mean_raw[5] + beta_raw5 * sd_raw[5];

}

model{
  mean_raw ~ std_normal();
  sd_raw ~ cauchy(0,3);

  beta_raw1 ~ std_normal();
  beta_raw2 ~ std_normal();
  beta_raw3 ~ std_normal();
  beta_raw4 ~ std_normal();
  beta_raw5 ~ std_normal();
  
  for (i in 1:n_s){
      vector[nTrials] utility;
      vector[nTrials] x_star = (ll_money[i] + ss_money[i]) / 2;
      vector[nTrials] t_star = (ll_time[i] + ss_time[i]) / 2;

      utility = beta[1][i] + 
                beta[2][i] * (ll_money[i] - ss_money[i]) +
                beta[3][i] * ((ll_money[i] - ss_money[i]) ./ x_star) +
                beta[4][i] * (ll_time[i] - ss_time[i]) +
                beta[5][i] * ((ll_time[i] - ss_time[i]) ./ t_star);
      choice[i] ~ bernoulli_logit(utility);
  }
  
}
1 Like