Is there a BRMs equivalent to this model?

Modeling
1

So I’m trying to figure out if there is BRMs equivalent to the following model I made that tries to figure out whether the expected value of a cue predicts how useful that cue is perceived to be:

The Data: I have an experiment in which participants completed a forced choice paradigm in which they made repeated decisions between 3 cues under different incentives (reward conditions). After completing a block of trials under each reward condition, they rated how useful each of the 3 cues were on a scale from 0% to 100%. I then calculate the expected value of each cue in each condition (normalized and centered on the expected value of cue 2, QCue2). Attached is data simulated for 10 participants:

Simulated_Data.csv (17.8 KB)
Columns are labeling in the model below.

The Model: To see if people rate usefulness in proportion to the relative value of each cue, the model runs expected value of the 3 cues in each condition through a softmax with a inverse temperature parameter, bEV, to generate the proportion of the relative value of each cue, p. After this the model, looks to see if this proportion matches the actual perceived usefulness using a standard linear model. All parameters are nested within subjects. The code I have made is like so:

Code <- "
data{
int N; // number of observations
int Usefulness[N]; // outcome
int NSub; //number of subjects
int Subject[N]; //Current Subject
Cue2 QCue1[N]; //EV of Cue1
Cue2 QCue3[N]; //EV of Cue3
Cue2 Cue1[N]; //Indicator Variable for Cue 1
Cue2 Cue2[N]; //Indicator Variable for Cue 2
Cue2 Cue3[N]; //Indicator Variable for Cue 3
}

parameters{
  matrix[4,NSub] ZSub;
  Cue2 ACue1Bar; // group mean of the intercept on Cue 1
  Cue2 ACue2Bar; // group mean of the intercept on Cue 2
  Cue2 ACue3Bar; // group mean of the intercept on Cue 3
  Cue2 bEVBar; // group mean of inverse temperature on EV
  vector<lower=0>[4] SigmaSub; // subject variance of each parameter
  cholesky_factor_corr[4] L_Rho_Sub; // Cholesky factor of each parameter
  Cue2<lower=0> Sigma;
}

transformed parameters{
  vector[NSub]  ACue1; // subject level intercepts on Cue 1
  vector[NSub]  ACue2; // subject level intercepts on Cue 2
  vector[NSub]  ACue3; // subject level intercepts on Cue 3
  vector[NSub]  bEV; // mean of beta on EV
  matrix[NSub, 4] Params;
  Params = (diag_pre_multiply(SigmaSub, L_Rho_Sub) * ZSub)'; // calculating subject-level parameter values 
  for ( j in 1:NSub ){ // assigning subject-level parameter values 
    ACue1[j] = Params[j,1] + ACue1Bar;
    ACue2[j] = Params[j,2] + ACue2Bar;
    ACue3[j] = Params[j,3] + ACue3Bar;
    bEV[j] = Params[j,4] + bEVBar;
  }
}

model{
  vector[N] Mu;
  vector[3] s;
  vector[3] p;
  
  to_vector(ZSub) ~ normal(0, 1);
  L_Rho_Sub ~ lkj_corr_cholesky( 2 );
  SigmaSub ~ exponential( 1 );
  Sigma ~ exponential( 1 );
  
  ACue1Bar ~ normal( 0 , .5);
  ACue2Bar ~ normal( 0 , .5);
  ACue3Bar ~ normal( 0 , .5);
  bEVBar ~ normal( 0 , 1 );
  
  for ( i in 1:N ) {
    s[1] =  bEV[Subject[i]]*QCue1[i];
    s[2] = 0;
    s[3] = bEV[Subject[i]]*QCue3[i];
    p = softmax(s);
    Mu[i] = ACue1[Subject[i]]*Cue1[i] + p[1]*Cue1[i] +
            ACue2[Subject[i]]*Cue2[i] + p[2]*Cue2[i] +
            ACue3[Subject[i]]*Cue3[i] + p[3]*Cue3[i];
            
            Usefulness[i] ~ normal(Mu[i],Sigma);
  }
}
}
"

dat_list_ <- list( N= nrow(d), Usefulness = d$Usefulness, NSub = length(unique(d$subject)), Subject = d$subject,
                               QCue1 =d$QCue1, QCue3 =d$QCue3, Cue1 = d$Cue1, Cue2 = d$Cue2, Cue3 = d$Cue3)
Model <- stan( model_code=Code , data=dat_list,
                      chains=4, cores = 4, control=list(max_treedepth=15), warmup = 1000, iter = 4000, seed = 123)

I’m trying to figure out if this model is possible to do in BRMs? I haven’t been able to find any resources documenting any similar models that involve first computing quantities (according some parameters) that will then later be used in regression.