Remove/estimate correlation of parameters (or get to know if correlation is real)

Hi all,

I have a basic question. I was playing with [RA-prospect theory model](https://github.com/CCS-Lab/hBayesDM/blob/master/commons/stan_files/ra_prospect.stan, code and data follow) and I’m getting a correlation between my fitted parameters (see Fig.

). I have no or very little experience with this, so I would have several questions and would be very grateful for advice and explanations:

  1. I have been reading other posts (eg., Infinite location parameter warning and correlated parameters, Reducing correlation between parameters of the posterior distribution, …) but I guess I need a bit more explanation and help.
    What I am looking at and what I’m plotting here are the correlations of means of the estimated parameters. A naive question is, is this the same correlation as which would be between the parameters (ie. if I take the distribution of parameters for each subject and look at correlations within a subject)? Why/why not? Where would be the difference?
    Second question would be, what does it mean if the correlation gets more significant if I include more samples? See Fig.

Btw. I’m seeing that eg. \rho and \tau have nice and clear exponential dependency. See points bellow. The plots are just illustrative.

  1. What was proposed by several people were three things:
    A) reparametrise the equation so that the correlation is removed (do some math transformation, as in the first post) – that I am unable to do since the equation I’m fitting is quite complex:

sv = \text{gain}^\rho - \lambda \cdot \text{loss}^\rho,

followed by:

y = \frac{e^{\tau(sv-c^\rho)}}{1+e^{\tau(sv-c^\rho)}},

where parameters are \tau inverse temperature, \lambda loss aversion coefficient, \rho risk aversion coefficient; and c being some given constant, and gains, losses the inputs.

B) To try to estimate the correlation within STAN and then account for it. Here I’d have two questions – how to do that (in STAN, pySTAN) and if I have it, how to use it – what will it tell me and how should I interpret it?

C) To try to use different priors, eg. multinominal/multigaussian. Here I’d have also two questions – how do I select my priors? How should I create them? And then how and why does it work? If I have any correlation of the parameters left after this, does that mean that actually the parameters are correlated or does it still mean that maybe my model is wrong?

Note: I’m checking convergence of the parameters (n_eff, Rhat, divergence, tree_depth), all are ok.


The stan code I’m using is the following (taken from the page mentioned before):

data {
  int<lower=1> N;
  int<lower=1> T;
  int<lower=1, upper=T> Tsubj[N];
  real<lower=0> gain[N, T];
  real<lower=0> loss[N, T];  // absolute loss amount
  real cert[N, T];
  int<lower=-1, upper=1> gamble[N, T];
}
transformed data {
}
parameters {
  vector[3] mu_pr;
  vector<lower=0>[3] sigma; 
  vector[N] rho_pr;
  vector[N] lambda_pr;
  vector[N] tau_pr;
}
transformed parameters {
  vector<lower=0, upper=2>[N]  rho;
  vector<lower=0, upper=5>[N]  lambda;
  vector<lower=0, upper=30>[N] tau;

  for (i in 1:N) {
    rho[i]    = Phi_approx(mu_pr[1] + sigma[1] * rho_pr[i]) * 2;
    lambda[i] = Phi_approx(mu_pr[2] + sigma[2] * lambda_pr[i]) * 5;
    tau[i]    = Phi_approx(mu_pr[3] + sigma[3] * tau_pr[i]) * 30;
  }
}
model {
  // ra_prospect: Original model in Soko-Hessner et al 2009 PNAS
  // hyper parameters
  mu_pr  ~ normal(0, 1.0);
  sigma ~ normal(0, 0.2);

  // individual parameters w/ Matt trick
  rho_pr    ~ normal(0, 1.0);
  lambda_pr ~ normal(0, 1.0);
  tau_pr    ~ normal(0, 1.0);

  for (i in 1:N) {
    for (t in 1:Tsubj[i]) {
      real evSafe;    // evSafe, evGamble, pGamble can be a scalar to save memory and increase speed.
      real evGamble;  // they are left as arrays as an example for RL models.
      real pGamble;

      // loss[i, t]=absolute amount of loss (pre-converted in R)
      evSafe   = pow(cert[i, t], rho[i]);
      evGamble = 0.5 * (pow(gain[i, t], rho[i]) - lambda[i] * pow(loss[i, t], rho[i]));
      pGamble  = inv_logit(tau[i] * (evGamble - evSafe));
      gamble[i, t] ~ bernoulli(pGamble);
    }
  }
}
generated quantities {
  real<lower=0, upper=2>  mu_rho;
  real<lower=0, upper=5>  mu_lambda;
  real<lower=0, upper=30> mu_tau;

  real log_lik[N];

  // For posterior predictive check
  real y_pred[N, T];

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

  mu_rho    = Phi_approx(mu_pr[1]) * 2;
  mu_lambda = Phi_approx(mu_pr[2]) * 5;
  mu_tau    = Phi_approx(mu_pr[3]) * 30;

  { // local section, this saves time and space
    for (i in 1:N) {
      log_lik[i] = 0;
      for (t in 1:Tsubj[i]) {
        real evSafe;    // evSafe, evGamble, pGamble can be a scalar to save memory and increase speed.
        real evGamble;  // they are left as arrays as an example for RL models.
        real pGamble;

        evSafe     = pow(cert[i, t], rho[i]);
        evGamble   = 0.5 * (pow(gain[i, t], rho[i]) - lambda[i] * pow(fabs(loss[i, t]), rho[i]));
        pGamble    = inv_logit(tau[i] * (evGamble - evSafe));
        log_lik[i] += bernoulli_lpmf(gamble[i, t] | pGamble);

        // generate posterior prediction for current trial
        y_pred[i, t] = bernoulli_rng(pGamble);
      }
    }
  }
}

The code I’m using for fitting in python is the following:

# define conditions and path
import os, sys
import pystan

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

import time
from datetime import timedelta

import scipy.stats as stats

data_path = '/home/user/LossAversion/'
model_folder = "/home/user/LossAversion/stan_files/"

#prepare model
model = 'ra_prospect.stan'
model_path = os.path.join(model_folder, model)

summary_dict = None
start = time.time()

# Compile the model
start_compile = time.time()
stan_model_RA = pystan.StanModel(file=model_path)
end_compile = time.time()
print('Duration of compilation: ', timedelta(seconds=end_compile - start_compile))

# run the model
start_model = time.time()
fit_ra_cleared = stan_model_RA.sampling(data=RA_prep_test, 
                    chains=4, iter=2000, warmup=1000, thin=1, init='random', verbose=True, 
                    control = {"adapt_delta":0.95, "stepsize":1, "max_treedepth":10}, n_jobs=-1)
end_model = time.time()
print('Duration of the model run is: ', timedelta(seconds=end_model - start_model))

# Create an ordered dictionary
summary_dict = fit_ra_cleared.summary()

# Create a data frame out of the output
df_cleared = pd.DataFrame(summary_dict['summary'], 
                  columns=summary_dict['summary_colnames'], 
                  index=summary_dict['summary_rownames'])
    
end = time.time()

print('Duration of the entire process is: ', timedelta(seconds=end - start))

The subsequent analysis is just reshaping the results, creating a new Data Frame with columns being the parameters and rows subjects and plotting using seaborn and PairGrid. Statistics are done using scipy.stats and the Pearson and Spearman correlation from them. I can post this code as well if needed (might be), I do not now because it is already quite heavy.

The data are accessible here, data_LA_exp1_STAN_correlation.txt (160.5 KB), and contain a subset of my data.

Thank you all in advance.

I believe these are the same thing. The idea is that you currently have three vectors, rho_pr, lambda_pr, and tau_pr that are modeled as sampled independently from a normal distribution, and instead you can model them as sampled from a multivariate distribution to provide explicit inference on their correlations:

...
parameters {
  vector[3] mu_pr;
  vector<lower=0>[3] sigma; 
  corr_matrix[3] cors;
  vector[3] rho_lambda_tau_pr[N]; 
}
transformed parameters {
  vector<lower=0, upper=2>[N]  rho;
  vector<lower=0, upper=5>[N]  lambda;
  vector<lower=0, upper=30>[N] tau;

  for (i in 1:N) {
    rho[i]    = Phi_approx(rho_lambda_tau_pr[i][1]) * 2;
    lambda[i] = Phi_approx(rho_lambda_tau_pr[i][2]) * 5;
    tau[i]    = Phi_approx(rho_lambda_tau_pr[i][3]) * 30;
  }
}
model {
  mu_pr  ~ normal(0, 1.0);
  sigma ~ normal(0, 0.2);
  rho_lambda_tau_pr ~ multi_normal(mu_pr,quad_form_diag(cors, sigma));
  ...
}

To simplify things, the above is the “centered” version of the hierarchical model, so if you’ve found that a non-centered version samples better for your data you’ll want to tweak the above to implement the non-centered version. You should also take a look at the “Multivariate priors for hierarchical models” section of the Stan Users Guide and implement the more efficient cholesky factorization version (which can be applied to either the centered or non-centered variants)

Thanks for your reply and the reference. I have read it briefly and I have tried, I got the following error after compilation:

“”"
Traceback (most recent call last):
File “/home/user/anaconda3/lib/python3.7/multiprocessing/pool.py”, line 121, in worker
result = (True, func(*args, **kwds))
File “/home/user/anaconda3/lib/python3.7/multiprocessing/pool.py”, line 44, in mapstar
return list(map(*args))
File “stanfit4anon_model_1be6737ea3855e86a262b997b3445e8b_3706375220267016912.pyx”, line 373, in stanfit4anon_model_1be6737ea3855e86a262b997b3445e8b_3706375220267016912._call_sampler_star
File “stanfit4anon_model_1be6737ea3855e86a262b997b3445e8b_3706375220267016912.pyx”, line 406, in stanfit4anon_model_1be6737ea3855e86a262b997b3445e8b_3706375220267016912._call_sampler
RuntimeError: Initialization failed.

The above exception was the direct cause of the following exception:

RuntimeError                              Traceback (most recent call last)
<ipython-input-11-8484c5f09526> in <module>
     20 fit_ra_cleared_multi = stan_model_RA.sampling(data=LA_data_prep_cleared, 
     21                     chains=4, iter=2000, warmup=1000, thin=1, init='random', verbose=True,
---> 22                     control = {"adapt_delta":0.95, "stepsize":1, "max_treedepth":10}, n_jobs=-1)
     23 end_model = time.time()
     24 print('Duration of the model run is: ', timedelta(seconds=end_model - start_model))

~/anaconda3/lib/python3.7/site-packages/pystan/model.py in sampling(self, data, pars, chains, iter, warmup, thin, seed, init, sample_file, diagnostic_file, verbose, algorithm, control, n_jobs, **kwargs)
    811         call_sampler_args = izip(itertools.repeat(data), args_list, itertools.repeat(pars))
    812         call_sampler_star = self.module._call_sampler_star
--> 813         ret_and_samples = _map_parallel(call_sampler_star, call_sampler_args, n_jobs)
    814         samples = [smpl for _, smpl in ret_and_samples]
    815 

~/anaconda3/lib/python3.7/site-packages/pystan/model.py in _map_parallel(function, args, n_jobs)
     83         try:
     84             pool = multiprocessing.Pool(processes=n_jobs)
---> 85             map_result = pool.map(function, args)
     86         finally:
     87             pool.close()

~/anaconda3/lib/python3.7/multiprocessing/pool.py in map(self, func, iterable, chunksize)
    266         in a list that is returned.
    267         '''
--> 268         return self._map_async(func, iterable, mapstar, chunksize).get()
    269 
    270     def starmap(self, func, iterable, chunksize=None):

~/anaconda3/lib/python3.7/multiprocessing/pool.py in get(self, timeout)
    655             return self._value
    656         else:
--> 657             raise self._value
    658 
    659     def _set(self, i, obj):

RuntimeError: Initialization failed.

I would also like to ask about the generated quantities. They have not shown it in the example, how are they generated? What need to be changed?

So the model compiled fine but you get an error when you try to sample?

yes, sorry, was above this message:

INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_1be6737ea3855e86a262b997b3445e8b NOW.

Duration of compilation:  0:01:19.850109

---------------------------------------------------------------------------
RemoteTraceback                           Traceback (most recent call last)
RemoteTraceback: 

Hm, that looks like an error in how you’re packaging the data in python. I don’t have any experience with PyStan. @ahartikainen have any insight on what’s going wrong here?

Hi, looks like your initialization failed.

There are couple of things that might cause this.

Some parameter is not initialized, e.g. subitem in vector / matrix. This is quite common reason. Solution: make sure you initialize all the parameters, check for common error of using wrong index (e.g. N instead of n, where N is the total and n is loop variable) (edit: Here I mean, e.g. you create an empty matrix and then don’t add anything on the last column)

Some priors are too tight or are too far from the default range and sampler has hard time to find typical set from random init (unconstrained space, -2 to 2). Solution: Use better priors or if the priors are ok, use user given inits.

1 Like

Thank you both for your comments. But I’m not sure I understand. The initialisation, ie. all what @ahartikainen writes should be done in STAN. All what has changed is that now the parametes come from multivariate and not three separate normals. So the priors did not change (still the same bound, tight, center,…). And also the initialisation should be the same, souldn’t it? To verify that, I run the exact same code just with the previous model (STAN prospect, as before). The sama data, same commands. It runs. So where can be a problem?

Can you please post your latest model code?

Yes, of course, sorry, should have posted already before…

data {
  int<lower=1> N;
  int<lower=1> T;
  int<lower=1, upper=T> Tsubj[N];
  real<lower=0> gain[N, T];
  real<lower=0> loss[N, T];  // absolute loss amount
  real cert[N, T];
  int<lower=-1, upper=1> gamble[N, T];
}
transformed data {
}
parameters {
  vector[3] mu_pr;
  vector<lower=0>[3] sigma; 
  corr_matrix[3] cors;
  vector[3] rho_lambda_tau_pr[N]; 
}
transformed parameters {
  vector<lower=0, upper=2>[N]  rho;
  vector<lower=0, upper=5>[N]  lambda;
  vector<lower=0, upper=30>[N] tau;

  for (i in 1:N) {
    rho[i]    = Phi_approx(rho_lambda_tau_pr[i][1]) * 2;
    lambda[i] = Phi_approx(rho_lambda_tau_pr[i][2]) * 5;
    tau[i]    = Phi_approx(rho_lambda_tau_pr[i][3]) * 30;
  }
}
model {
  mu_pr  ~ normal(0, 1.0);
  sigma ~ normal(0, 0.2);

// individual parameters w/ Matt trick
  rho_lambda_tau_pr ~ multi_normal(mu_pr,quad_form_diag(cors, sigma));

  for (i in 1:N) {
    for (t in 1:Tsubj[i]) {
      real evSafe;    // evSafe, evGamble, pGamble can be a scalar to save memory and increase speed.
      real evGamble;  // they are left as arrays as an example for RL models.
      real pGamble;

      // loss[i, t]=absolute amount of loss (pre-converted in R)
      evSafe   = pow(cert[i, t], rho[i]);
      evGamble = 0.5 * (pow(gain[i, t], rho[i]) - lambda[i] * pow(loss[i, t], rho[i]));
      pGamble  = inv_logit(tau[i] * (evGamble - evSafe));
      gamble[i, t] ~ bernoulli(pGamble);
    }
  }
}
generated quantities {
  real<lower=0, upper=2>  mu_rho;
  real<lower=0, upper=5>  mu_lambda;
  real<lower=0, upper=30> mu_tau;

  real log_lik[N];

  // For posterior predictive check
  real y_pred[N, T];

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

  mu_rho    = Phi_approx(mu_pr[1]) * 2;
  mu_lambda = Phi_approx(mu_pr[2]) * 5;
  mu_tau    = Phi_approx(mu_pr[3]) * 30;

  { // local section, this saves time and space
    for (i in 1:N) {
      log_lik[i] = 0;
      for (t in 1:Tsubj[i]) {
        real evSafe;    // evSafe, evGamble, pGamble can be a scalar to save memory and increase speed.
        real evGamble;  // they are left as arrays as an example for RL models.
        real pGamble;

        evSafe     = pow(cert[i, t], rho[i]);
        evGamble   = 0.5 * (pow(gain[i, t], rho[i]) - lambda[i] * pow(fabs(loss[i, t]), rho[i]));
        pGamble    = inv_logit(tau[i] * (evGamble - evSafe));
        log_lik[i] += bernoulli_lpmf(gamble[i, t] | pGamble);

        // generate posterior prediction for current trial
        y_pred[i, t] = bernoulli_rng(pGamble);
      }
    }
  }
}

Try putting this in the model block:

    cors ~ lkj_corr(2) ;

(Edited with correct spelling of lkj_corr)

If you mean:

model {
  mu_pr  ~ normal(0, 1.0);
  sigma ~ normal(0, 0.2);
  cors ~ lkj_cor(2);

// individual parameters w/ Matt trick
  rho_lambda_tau_pr ~ multi_normal(mu_pr,quad_form_diag(cors, sigma));

  for (i in 1:N) {
...

then I’m getting the following error:

ValueError: Failed to parse Stan model 'anon_model_d125e4bdad2924d8f45350df09870368'. Error message:
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Probability function must end in _lpdf or _lpmf. Found distribution family = lkj_cor with no corresponding probability function lkj_cor_lpdf, lkj_cor_lpmf, or lkj_cor_log
 error in 'ra_prospect_multivariate.stan' at line 32, column 20
  -------------------------------------------------
    30:   mu_pr  ~ normal(0, 1.0);
    31:   sigma ~ normal(0, 0.2);
    32:   cors ~ lkj_cor(2);
                           ^
    33: 
  -------------------------------------------------

Add r in cor --> corr

Ah, thanks, that worked, model compiled, however, initialisation failed again:

INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_ce77c6561d58b5b1ce594dc10c5aa2c6 NOW.

Duration of compilation:  0:01:21.795654

---------------------------------------------------------------------------
RemoteTraceback                           Traceback (most recent call last)
RemoteTraceback: 
"""
Traceback (most recent call last):
  File "/home/jan/anaconda3/lib/python3.7/multiprocessing/pool.py", line 121, in worker
    result = (True, func(*args, **kwds))
  File "/home/jan/anaconda3/lib/python3.7/multiprocessing/pool.py", line 44, in mapstar
    return list(map(*args))
  File "stanfit4anon_model_ce77c6561d58b5b1ce594dc10c5aa2c6_1688361003633273438.pyx", line 373, in stanfit4anon_model_ce77c6561d58b5b1ce594dc10c5aa2c6_1688361003633273438._call_sampler_star
  File "stanfit4anon_model_ce77c6561d58b5b1ce594dc10c5aa2c6_1688361003633273438.pyx", line 406, in stanfit4anon_model_ce77c6561d58b5b1ce594dc10c5aa2c6_1688361003633273438._call_sampler
RuntimeError: Initialization failed.
"""

The above exception was the direct cause of the following exception:

RuntimeError                              Traceback (most recent call last)
<ipython-input-13-8484c5f09526> in <module>
     20 fit_ra_cleared_multi = stan_model_RA.sampling(data=LA_data_prep_cleared, 
     21                     chains=4, iter=2000, warmup=1000, thin=1, init='random', verbose=True,
---> 22                     control = {"adapt_delta":0.95, "stepsize":1, "max_treedepth":10}, n_jobs=-1)
     23 end_model = time.time()
     24 print('Duration of the model run is: ', timedelta(seconds=end_model - start_model))

~/anaconda3/lib/python3.7/site-packages/pystan/model.py in sampling(self, data, pars, chains, iter, warmup, thin, seed, init, sample_file, diagnostic_file, verbose, algorithm, control, n_jobs, **kwargs)
    811         call_sampler_args = izip(itertools.repeat(data), args_list, itertools.repeat(pars))
    812         call_sampler_star = self.module._call_sampler_star
--> 813         ret_and_samples = _map_parallel(call_sampler_star, call_sampler_args, n_jobs)
    814         samples = [smpl for _, smpl in ret_and_samples]
    815 

~/anaconda3/lib/python3.7/site-packages/pystan/model.py in _map_parallel(function, args, n_jobs)
     83         try:
     84             pool = multiprocessing.Pool(processes=n_jobs)
---> 85             map_result = pool.map(function, args)
     86         finally:
     87             pool.close()

~/anaconda3/lib/python3.7/multiprocessing/pool.py in map(self, func, iterable, chunksize)
    266         in a list that is returned.
    267         '''
--> 268         return self._map_async(func, iterable, mapstar, chunksize).get()
    269 
    270     def starmap(self, func, iterable, chunksize=None):

~/anaconda3/lib/python3.7/multiprocessing/pool.py in get(self, timeout)
    655             return self._value
    656         else:
--> 657             raise self._value
    658 
    659     def _set(self, i, obj):

RuntimeError: Initialization failed.

try model.sampling(..., init_r=0.1) or something similar (small value)

Ok, I think that model should not have correlation matrix (bounded from -1 to 1), but covariance matrix. I mean it is used in multi_normal with non-normalized diagonal.

I added print(quad_form_diag(cors, sigma)) to model block to see the values.

Changing init_r did not help.

I have tried to add it in the following way:

data {
  int<lower=1> N;
  int<lower=1> T;
  int<lower=1, upper=T> Tsubj[N];
  real<lower=0> gain[N, T];
  real<lower=0> loss[N, T];  // absolute loss amount
  real cert[N, T];
  int<lower=-1, upper=1> gamble[N, T];
}
transformed data {
}
parameters {
  vector[3] mu_pr;
  vector<lower=0>[3] sigma; 
  corr_matrix[3] cors;
  vector[3] rho_lambda_tau_pr[N]; 
}
transformed parameters {
  vector<lower=0, upper=2>[N]  rho;
  vector<lower=0, upper=5>[N]  lambda;
  vector<lower=0, upper=30>[N] tau;
  cov_matrix[3] Sigma; 
  Sigma = quad_form_diag(cors, sigma); 

  for (i in 1:N) {
    rho[i]    = Phi_approx(rho_lambda_tau_pr[i][1]) * 2;
    lambda[i] = Phi_approx(rho_lambda_tau_pr[i][2]) * 5;
    tau[i]    = Phi_approx(rho_lambda_tau_pr[i][3]) * 30;
  }
}
model {
  mu_pr  ~ normal(0, 1.0);
  sigma ~ normal(0, 0.2);
  cors ~ lkj_corr(3);
  
// individual parameters w/ Matt trick
  rho_lambda_tau_pr ~ multi_normal(mu_pr,Sigma);

  for (i in 1:N) {
    for (t in 1:Tsubj[i]) {
      real evSafe;    // evSafe, evGamble, pGamble can be a scalar to save memory and increase speed.
      real evGamble;  // they are left as arrays as an example for RL models.
      real pGamble;

      // loss[i, t]=absolute amount of loss (pre-converted in R)
      evSafe   = pow(cert[i, t], rho[i]);
      evGamble = 0.5 * (pow(gain[i, t], rho[i]) - lambda[i] * pow(loss[i, t], rho[i]));
      pGamble  = inv_logit(tau[i] * (evGamble - evSafe));
      gamble[i, t] ~ bernoulli(pGamble);
    }
  }
}
generated quantities {
  real<lower=0, upper=2>  mu_rho;
  real<lower=0, upper=5>  mu_lambda;
  real<lower=0, upper=30> mu_tau;

  real log_lik[N];

  // For posterior predictive check
  real y_pred[N, T];

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

  mu_rho    = Phi_approx(mu_pr[1]) * 2;
  mu_lambda = Phi_approx(mu_pr[2]) * 5;
  mu_tau    = Phi_approx(mu_pr[3]) * 30;

  { // local section, this saves time and space
    for (i in 1:N) {
      log_lik[i] = 0;
      for (t in 1:Tsubj[i]) {
        real evSafe;    // evSafe, evGamble, pGamble can be a scalar to save memory and increase speed.
        real evGamble;  // they are left as arrays as an example for RL models.
        real pGamble;

        evSafe     = pow(cert[i, t], rho[i]);
        evGamble   = 0.5 * (pow(gain[i, t], rho[i]) - lambda[i] * pow(fabs(loss[i, t]), rho[i]));
        pGamble    = inv_logit(tau[i] * (evGamble - evSafe));
        log_lik[i] += bernoulli_lpmf(gamble[i, t] | pGamble);

        // generate posterior prediction for current trial
        y_pred[i, t] = bernoulli_rng(pGamble);
      }
    }
  }
}

as I have found in this post:
http://stla.github.io/stlapblog/posts/StanLKJprior.html
It compiles but still fails to initialise.

EDIT:
I have also tried to implement the Cholesky decomposition in the following way (with the same result, compiles but fails to initialise):

data {
  int<lower=1> N;            // num of subjects
  int<lower=1> T;            // num of trials per subject
  int<lower=1, upper=T> Tsubj[N];  // maximum possible num of trials
  real<lower=0> gain[N, T];
  real<lower=0> loss[N, T];  // absolute loss amount
  real cert[N, T];
  int<lower=-1, upper=1> gamble[N, T];
}
transformed data {
}
parameters {
  vector[3] mu_pr;
  vector<lower=0>[3] sigma; 
  cholesky_factor_corr[3] L_Omega;
  //corr_matrix[3] cors;
  vector[3] rho_lambda_tau_pr[N]; 
  
}
transformed parameters {
  vector<lower=0, upper=2>[N]  rho;
  vector<lower=0, upper=5>[N]  lambda;
  vector<lower=0, upper=30>[N] tau;

  for (i in 1:N) {
    rho[i]    = Phi_approx(rho_lambda_tau_pr[i][1]) * 2;
    lambda[i] = Phi_approx(rho_lambda_tau_pr[i][2]) * 5;
    tau[i]    = Phi_approx(rho_lambda_tau_pr[i][3]) * 30;
  }
}
model {
  mu_pr  ~ normal(0, 1.0);
  sigma ~ normal(0, 0.2);
  //cors ~ lkj_cor(2);
  L_Omega ~ lkj_corr_cholesky(3);

// individual parameters w/ Matt trick
  rho_lambda_tau_pr ~  multi_normal_cholesky(mu_pr,quad_form_diag(L_Omega, sigma));

  for (i in 1:N) {
    for (t in 1:Tsubj[i]) {
      real evSafe;    // evSafe, evGamble, pGamble can be a scalar to save memory and increase speed.
      real evGamble;  // they are left as arrays as an example for RL models.
      real pGamble;

      // loss[i, t]=absolute amount of loss (pre-converted in R)
      evSafe   = pow(cert[i, t], rho[i]);
      evGamble = 0.5 * (pow(gain[i, t], rho[i]) - lambda[i] * pow(loss[i, t], rho[i]));
      pGamble  = inv_logit(tau[i] * (evGamble - evSafe));
      gamble[i, t] ~ bernoulli(pGamble);
    }
  }
}
generated quantities {
  real<lower=0, upper=2>  mu_rho;
  real<lower=0, upper=5>  mu_lambda;
  real<lower=0, upper=30> mu_tau;
  
  matrix[3,3] Omega;
  matrix[3,3] Sigma;

  real log_lik[N];

  // For posterior predictive check
  real y_pred[N, T];

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

  mu_rho    = Phi_approx(mu_pr[1]) * 2;
  mu_lambda = Phi_approx(mu_pr[2]) * 5;
  mu_tau    = Phi_approx(mu_pr[3]) * 30;

  Omega = multiply_lower_tri_self_transpose(L_Omega);
  Sigma = quad_form_diag(Omega, sigma); 

  { // local section, this saves time and space
    for (i in 1:N) {
      log_lik[i] = 0;
      for (t in 1:Tsubj[i]) {
        real evSafe;    // evSafe, evGamble, pGamble can be a scalar to save memory and increase speed.
        real evGamble;  // they are left as arrays as an example for RL models.
        real pGamble;

        evSafe     = pow(cert[i, t], rho[i]);
        evGamble   = 0.5 * (pow(gain[i, t], rho[i]) - lambda[i] * pow(fabs(loss[i, t]), rho[i]));
        pGamble    = inv_logit(tau[i] * (evGamble - evSafe));
        log_lik[i] += bernoulli_lpmf(gamble[i, t] | pGamble);

        // generate posterior prediction for current trial
        y_pred[i, t] = bernoulli_rng(pGamble);
      }
    }
  }
}

I have no experience with R so this might have even more problems but I have tried the following code also in R. And I got a compilation error as well so maybe the problem is in the STAN code, any ideas @ahartikainen, @mike-lawrence?

> fit <- stan(file = stan_model1, data = stan_data, chains=4, iter=2000, warmup=1000, thin=1, init='random')
Error in compileCode(f, code, language = language, verbose = verbose) : 
  Compilation ERROR, function(s)/method(s) not created! file25313aa703c2.cpp:6:9: warning: ISO C++11 requires whitespace after the macro name
    6 | #define STAN__SERVICES__COMMAND_HPP#include <boost/integer/integer_log2.hpp>
      |         ^~~~~~~~~~~~~~~~~~~~~~~~~~~
In file included from /apps/RH7U2/gnu/R/3.5.0/lib64/R/library/RcppEigen/include/Eigen/Core:383,
                 from /apps/RH7U2/gnu/R/3.5.0/lib64/R/library/RcppEigen/include/Eigen/Dense:1,
                 from /apps/RH7U2/gnu/R/3.5.0/lib64/R/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4,
                 from <command-line>:
/apps/RH7U2/gnu/R/3.5.0/lib64/R/library/RcppEigen/include/Eigen/src/Core/arch/SSE/PacketMath.h:60:39: warning: ignoring attributes on template argument ‘__m128’ {aka ‘__vector(4) float’} [-Wignored-attributes]
   60 | template<> struct is_arithmetic<__m128>  { enum { value = true }; };
      |                                       ^
/apps/RH7U2/gnu/R/3.5.0/lib64/R
In addition: Warning message:
In system(cmd, intern = !verbose) :
  running command '/apps/RH7U2/gnu/R/3.5.0/lib64/R/bin/R CMD SHLIB file25313aa703c2.cpp 2> file25313aa703c2.cpp.err.txt' had status 1
Error in sink(type = "output") : invalid connection

The code what I used was:

library(rstan)
stanc("ra_prospect_multivariate.stan")
stan_model1 <-"ra_prospect_multivariate.stan"
stan_data <- list(N = 1,T = 256,Tsubj = c(256),gain = c(10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 10, 12, 14, 16, 
                                                        18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 10,12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 10, 12, 14, 16, 18, 20, 
                                                        22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 10, 12, 14, 16, 18, 20, 22, 24,
                                                        26, 28, 30, 32, 34, 36, 38, 40, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 
                                                        30, 32, 34, 36, 38, 40, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 
                                                        34, 36, 38, 40, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40),loss = c(40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40,40, 40, 40, 38, 38, 38, 38, 38, 38, 38, 38, 38, 38,38, 38, 38, 
                                                                                                                                                 38, 38, 38, 36, 36, 36, 36, 36, 36, 36,36, 36, 36, 36, 36, 36, 36, 36, 36, 34, 34, 34, 34,34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 32,32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32,32, 32, 30,
                                                                                                                                                 30, 30, 30, 30, 30, 30, 30, 30, 30, 30,30, 30, 30, 30, 30, 28, 28, 28, 28, 28, 28, 28, 28,28, 28, 28, 28, 28, 28, 28, 28, 26, 26, 26, 26, 26,26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 24, 24,24, 24, 24, 
                                                                                                                                                 24, 24, 24, 24, 24, 24, 24, 24, 24, 24,24, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22,22, 22, 22, 22, 20, 20, 20, 20, 20, 20, 20, 20, 20,20, 20, 20, 20, 20, 20, 20, 18, 18, 18, 18, 18, 18,18, 18, 18, 
                                                                                                                                                 18, 18, 18, 18, 18, 18, 18, 16, 16, 16,16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16,14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14,14, 14, 14, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12,12, 12, 12, 
                                                                                                                                                 12, 12, 12, 10, 10, 10, 10, 10, 10, 10,10, 10, 10, 10, 10, 10, 10, 10, 10),cert = rep(5, 256),gamble = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                                                                                                                                                                                                                          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                                                                                                                                                                                                                          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                                                                                                                                                                                                                          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                                                                                                                                                                                                                          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                                                                                                                                                                                                                          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                                                                                                                                                                                                                          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                                                                                                                                                                                                                          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                                                                                                                                                                                                                          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                                                                                                                                                                                                                          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                                                                                                                                                                                                                          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                                                                                                                                                                                                                                                          0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

fit <- stan(file = stan_model1, data = stan_data, chains=4, iter=2000, warmup=1000, thin=1, init='random')

Notice that I used just one subject unlike in the python code. Just to be clear that the data is a bit different.

Yeah, that looks like something is wrong with how you installed Stan. What OS are you using? How did you install Stan?