Models with within-subject parameters

For example, a group of subjects underwent two treatments. The parameters I am interested in are the learning rate (LR) and inverse temperature (IT) of the basic reinforcement learning model. I want to model the within-subject change of parameters LR and IT(Treatment1 and Treatment2).

Since LR should be between 0 and 1, I used a logistic transformation for LR in treatment 1: LR1=1/[1+exp(LR_norm)], in which LR_norm is a vector with elements following a normal distribution.

For LR in treatment 2, I used LR_norm and LR_diff to represent the term within the exp(). Therefore LR2=1/[1+exp(LR_norm+LR_diff)].

I use a similar way to model IT in treatment 1 and treatment 2:

IT1=exp(IT_norm);

IT2=exp(IT_norm+IT_diff);

Do you see anything wrong with this?

The weird thing happens when I do this:

I use 30 LR values and 30 IT values to simulate data for 30 subjects in treatment 1, and the same 30 LR and IT values to simulate data for these subjects in treatment 2. If everything is correct, the estimated LR and IT values in the two treatment conditions shouldn’t be significantly different. However, this model ALWAYS gives highly significant treatment effects( e.g. a paired t-test of LR gives t = -8.0381, df = 29, p-value = 7.278e-09).

The full model and the simulation scripts are attached if you are interested in running the model. You only need to modify directories in the R file:
simulation.R (9.2 KB)
stan_1LR1IT_withinsub.stan (3.9 KB)

The whole model is as follows:

//#####################################################
data {
int T; //number of trials
int S; //number of subjects

int choice[T]; //choice: 1 or 2
vector[T] out; //outcome values
vector[T] block_num; //vector (value per trial) of block number per subject
vector[T] subj; //vector (value per trial) of subject number
}

//#####################################################
parameters {
real LR_m;
real<lower=0> LR_s;
vector[S] LR_raw;

real LR_diff_m;
real<lower=0> LR_diff_s;
vector[S] LR_diff_raw;
//-----------------------------
real IT_m;
real <lower=0> IT_s;
vector[S] IT_raw;

real IT_diff_m;
real<lower=0> IT_diff_s;
vector[S] IT_diff_raw;
}
//#####################################################
transformed parameters {
vector[S] LR_norm;
vector[S] IT_norm;

vector[S] LR1;
vector[S] LR2;
vector[S] IT1;
vector[S] IT2;

vector[S] LR_diff;
vector[S] IT_diff;

LR_norm=LR_m+LR_sLR_raw;
IT_norm=IT_m+IT_s
IT_raw;

LR_diff=LR_diff_m+LR_diff_sLR_diff_raw;
IT_diff=IT_diff_m+IT_diff_s
IT_diff_raw;

LR1=rep_vector(1,S) ./ (1+exp(LR_norm));
LR2=rep_vector(1,S) ./ (1+exp(LR_norm+LR_diff));
));

IT1=exp(IT_norm);
IT2=exp(IT_norm+IT_diff);
}

//#####################################################
model {
vector[T+1] Q1; // Q1 is EV of the better option,size T+1 to contain the last trial
vector[T+1] Q2;
//------------------------------
LR_m ~ normal(0,5);
LR_s ~ cauchy(0,5);
LR_raw ~ normal(0,1);
LR_diff_m ~ normal(0,5);
LR_diff_s ~ cauchy(0,5);
LR_diff_raw ~ normal(0,1);
//------------------------------
IT_m ~ normal(0,5);
IT_s ~ cauchy(0,5);
IT_raw ~ normal(0,1);
IT_diff_m ~ normal(0,5);
IT_diff_s ~ cauchy(0,5);
IT_diff_raw ~ normal(0,1);

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

{int s; //1 to start, subject #
s = 1;
for (t in 1:T) {
// ********************** 1st trial ****************************************
if (t==1) {
Q1[t] = 0;
Q2[t] = 0;}
// ************************** new subject **********************************
else if (subj[t] != subj[t-1]) {
s = s+1; //update subject
Q1[t] = 0;
Q2[t] = 0;}
// ************************** new block ************************************
else if (block_num[t] != block_num[t-1]) { //same subj, different block
Q1[t] = 0;
Q2[t] = 0;}
// ************************** same subject, same block **********************
else {}
//-----------------------------------------------------------------------------------------------------
if (s<S+1) {

      choice[t] ~ bernoulli_logit(IT1[s]*(Q1[t]-Q2[t]));
      
      if (choice[t]==1) {			//EV update, only option selected updated
        Q1[t+1] = Q1[t]+LR1[s]*(out[t]-Q1[t]);
        Q2[t+1] = Q2[t];}
      else if (choice[t]==0) {
        Q2[t+1] = Q2[t]+LR1[s]*(out[t]-Q2[t]);
        Q1[t+1] = Q1[t];}
    } else if (s>S){
      
      choice[t] ~ bernoulli_logit(IT2[s-S]*(Q1[t]-Q2[t]));
      
      if (choice[t]==1) {			//EV update, only option selected updated
        Q1[t+1] = Q1[t]+LR2[s-S]*(out[t]-Q1[t]);
        Q2[t+1] = Q2[t];}
      else if (choice[t]==0) {
        Q2[t+1] = Q2[t]+LR2[s-S]*(out[t]-Q2[t]);
        Q1[t+1] = Q1[t];} 
    }

}	// end of trial 

} //end of group
} // end of model section

There must be something wrong or some bias in the model.
I appreciate it if you can give me any suggestions.

Thanks!
Sheng