Hierarchical Bayes Model with Multiple Between-Subjects Conditions

Hi,

I am trying to fit a Hierarchical Bayes model emulating the hBayesDM approach. It is working and I have three studies, and two of the studies have multiple different between subjects conditions. Right now, for any between-subjects inference, I just fit the models separately for each condition and then compare their group distributions. This is what hBayesDM vignette advises as well.

However, I am thinking it may be better to assume they are all conditions drawn from the same population distribution, right? And treat each distribution as it’s own separate group-level distribution within the population-level distribution? How might I accomplish this?

Below is my code:

data {
  int<lower=1> nSubjects; // number of subjects total
  int<lower = 1> maxTrials; //number of testing/generalization trials
  int<lower = 1> maxTrain; // max number of training trials
  int<lower = 1> nTrain[nSubjects]; // per participant number of training trials
  int<lower = 1> nTrials[nSubjects]; //per participant number of testing/generalization trials
  int<lower = 0, upper = 2> groupChoice[nSubjects,maxTrials]; // which group chosen
  real sg[nSubjects, maxTrials]; // self-weighted similarity
  
  vector[maxTrain] prevSim[nSubjects, maxTrials]; // 3-dimensional matrix of subject's training traits similarities to testing traits
  //int<lower = 0, upper = 7> prevSelf[nSubjects,maxTrain]; // matrix of training self-evaluations
  vector[maxTrain] prevSelf[nSubjects]; // matrix of training self-evaluations
  
}


parameters {
  // Declare all parameters as vectors for vectorizing
  // Hyper(group)-parameters
  vector[3] mu_pr;
  vector<lower=0>[3] sigma;

  // Subject-level raw parameters (for Matt trick)
  vector[nSubjects] tau_pr;  // inverse temperature
  vector[nSubjects] m_pr;  // slope for ingroup
  vector[nSubjects] bias_pr;  // slope for bias
  
}

transformed parameters {
  // subject-level parameters
  vector<lower=0, upper=10>[nSubjects] tau; # temperature
  vector<lower=0, upper=10>[nSubjects] m; # growth rate
  vector<lower=0, upper=1>[nSubjects] bias; # bias

  for (i in 1:nSubjects) {
    
    tau[i] = Phi_approx(mu_pr[1] + sigma[1] * tau_pr[i]) * 10; 
    m[i] = Phi_approx(mu_pr[2] + sigma[2] * m_pr[i]) * 10; 
    bias[i] = Phi_approx(mu_pr[3] + sigma[3] * bias_pr[i]); 
    
  }
  
}

model {
  // Hyperparameters
  mu_pr  ~ normal(0, 1);
  sigma ~ normal(0, 0.2);
  //sigma[3:4] ~ cauchy(0, 0.35);

  // individual parameters
  tau_pr ~ normal(0, 1);
  m_pr ~ normal(0, 1);
  bias_pr ~ normal(0, 1);

  for (s in 1:nSubjects) {

    vector[2] simW;
    vector[2] prob;
    vector[nTrain[s]] curSelf;
    vector[nTrain[s]] GPin;
    vector[nTrain[s]] GPout;
    vector[nTrain[s]] PS;

    GPin[1:nTrain[s]] = rep_vector(1,nTrain[s])./(1 + exp((-m[s])*(prevSelf[s,1:nTrain[s]]-4))); # logistic function for transforming eval ratings to ingroup
    GPout[1:nTrain[s]] = rep_vector(1,nTrain[s])./(1 + exp((m[s])*(prevSelf[s,1:nTrain[s]]-4))); # logistic function for transforming eval ratings to outgroup
    
    for (t in 1:nTrials[s]) {
      
      PS[1:nTrain[s]] = prevSim[s,t,1:nTrain[s]]; // similarities from training to current generalization trait
      
      simW[1] = dot_product(GPout[1:nTrain[s]],PS); # dot-product of ingroup and sim
      simW[2] = dot_product(GPin[1:nTrain[s]],PS); #dot-product of ingroup and sim
      
      prob[1] = ( (1 - bias[s]) * pow( ( simW[1] ) ,tau[s] ) ) / ( ( (1 - bias[s]) *  pow( simW[1]  , tau[s] ) ) + ( (bias[s]) * pow( ( simW[2] ) ,tau[s] ) ) ); // convert to probabilities
      prob[2] = ( (bias[s]) * pow( ( simW[2] ) ,tau[s] ) ) / ( ( (1 - bias[s]) *  pow( simW[1]  , tau[s] ) ) + ( (bias[s]) * pow( ( simW[2] ) ,tau[s] ) ) ); // convert to probabilities
      
      groupChoice[s,t] ~ categorical( prob );
      
    }
  }    
}

generated quantities {
  // For group level parameters
  real<lower=0, upper=10> mu_tau;
  real<lower=0, upper=10> mu_m;
  real<lower=0, upper=1> mu_bias;

  // For log likelihood calculation
  real log_lik[nSubjects];

  // For posterior predictive check
  real y_pred[nSubjects, maxTrials];

  // Set all posterior predictions to 0 (avoids NULL values)
  for (i in 1:nSubjects) {
    for (t in 1:maxTrials) {
      y_pred[i, t] = -1;
    }
  }

  //mu_A   = Phi_approx(mu_pr[1]);
  mu_tau = Phi_approx(mu_pr[1]) * 10;
  mu_m   = Phi_approx(mu_pr[2]) * 10;
  mu_bias   = Phi_approx(mu_pr[3]);

  { // local section, this saves time and space
    
    for (s in 1:nSubjects) {

    vector[2] simW;
    vector[2] prob;
    vector[nTrain[s]] curSelf;
    vector[nTrain[s]] GPin;
    vector[nTrain[s]] GPout;
    vector[nTrain[s]] PS;
    
    log_lik[s] = 0;
    
    GPin[1:nTrain[s]] = rep_vector(1,nTrain[s])./(1 + exp((-m[s])*(prevSelf[s,1:nTrain[s]]-4)));
    GPout[1:nTrain[s]] = rep_vector(1,nTrain[s])./(1 + exp((m[s])*(prevSelf[s,1:nTrain[s]]-4)));
    
    for (t in 1:nTrials[s]) {

      PS[1:nTrain[s]] = prevSim[s,t,1:nTrain[s]]; // similarities from training to current generalization trait
      
      simW[1] = dot_product(GPout[1:nTrain[s]],PS);
      simW[2] = dot_product(GPin[1:nTrain[s]],PS);
      
      prob[1] = ( (1 - bias[s]) * pow( ( simW[1] ) ,tau[s] ) ) / ( ( (1 - bias[s]) *  pow( simW[1]  , tau[s] ) ) + ( (bias[s]) * pow( ( simW[2] ) ,tau[s] ) ) ); // convert to probabilities
      prob[2] = ( (bias[s]) * pow( ( simW[2] ) ,tau[s] ) ) / ( ( (1 - bias[s]) *  pow( simW[1]  , tau[s] ) ) + ( (bias[s]) * pow( ( simW[2] ) ,tau[s] ) ) ); // convert to probabilities
        
        // compute log likelihood of current trial
        log_lik[s] += categorical_lpmf( groupChoice[s, t] | prob );

        // generate posterior prediction for current trial
        y_pred[s, t] = categorical_rng(prob);
        
      }

  }   
    
  }
}

I am thinking maybe I could modify it so that vector[3] mu_pr denoting the group-level estimates for the parameters is the overall population, and then there are vector[nGroups] sub_mu_pr denoting sub groups for as many conditions/studies as there were.

Does anyone have advice on how to implement a hierarchical bayes model where instead of estimating the model separately each time for each condition/study, it’s estimated simultaneously for all groups and treats conditons/between-subjects groups as being drawn from an overall group distribution?

Thank you for any advice!