# How can I solve the problem of high Pareto-k values?

**URL:** <https://discourse.mc-stan.org/t/how-can-i-solve-the-problem-of-high-pareto-k-values/31774>\
**Category:** Modeling\
**Created:** [June 18, 2023, 7:11am UTC](https://discourse.mc-stan.org/t/how-can-i-solve-the-problem-of-high-pareto-k-values/31774 "2023-06-18T07:11:51Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![asyy](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/asyy/32/16522_2.png) [@asyy](https://discourse.mc-stan.org/u/asyy)\
**Post date:** [June 18, 2023, 7:11am UTC](https://discourse.mc-stan.org/t/how-can-i-solve-the-problem-of-high-pareto-k-values/31774/1 "2023-06-18T07:11:51Z")

</div>

My data is the data of reversal learning task consisting of 4 blocks. There are two choices(1 and 2) and the values of two choices reversed when a block is switched to the next block. I have 52 participants and in most cases, participants completed 120 trials. I am trying to fit a hierarchical reinforcement learning model to my data using stan. The problem is that when I use loo function to compute the looic value, most of Pareto k diagnostic values are bad or very bad as follows:

## Computed from 12000 by 52 log-likelihood matrix Estimate SE elpd\_loo -2282.4 155.6 p\_loo 65.4 3.4 looic 4564.8 311.1

Monte Carlo SE of elpd\_loo is NA.

Pareto k diagnostic values:  
Count Pct. Min. n\_eff  
(-Inf, 0.5] (good) 0 0.0%   
(0.5, 0.7] (ok) 8 15.4% 112  
(0.7, 1] (bad) 36 69.2% 13  
(1, Inf) (very bad) 8 15.4% 12  
See help(‘pareto-k-diagnostic’) for details.

**Does this mean that my model is misspecified? Is there any problem in my code? I attached my data and code.**

Thank you for helping me.

[RLdata.txt](https://discourse.mc-stan.org/uploads/short-url/uzGjJ3mY081Uw0BhuPxovDF1QfW.txt) (171.7 KB)

> data {  
> int\<lower=1\> N; // Number of subjects
> 
> int\<lower=1\> B; // Maximum number of blocks across subjects  
> int\<lower=1\> Bsubj[N]; // Number of blocks for each subject
> 
> int\<lower=0\> T; // Maximum number of trials across subjects  
> int\<lower=0, upper=T\> Tsubj[N, B]; // Number of trials/blocks for each subject
> 
> int\<lower=-1, upper=2\> choice[N, B, T]; // The choices subjects made  
> real outcome[N, B, T]; // The outcome  
> }
> 
> transformed data {  
> // Default value for (re-)initializing parameter vectors  
> vector[2] initV;  
> initV = rep\_vector(0.0, 2);  
> }
> 
> // Declare all parameters as vectors for vectorizing  
> parameters {  
> // Hyper(group)-parameters  
> vector[2] mu\_pr;  
> vector\<lower=0\>[2] sigma;
> 
> // Subject-level raw parameters (for Matt trick)  
> vector[N] alpha\_pr; // learning rate  
> vector[N] beta\_pr; // inverse temperature  
> }
> 
> transformed parameters {  
> // Transform subject-level raw parameters  
> vector\<lower=0, upper=1\>[N] alpha;  
> vector\<lower=0, upper=10\>[N] beta;
> 
> for (i in 1:N) {  
> alpha[i] = Phi\_approx(mu\_pr[1] + sigma[1] \* alpha\_pr[i]);  
> beta[i] = Phi\_approx(mu\_pr[2] + sigma[2] \* beta\_pr[i]) \* 10;  
> }  
> }
> 
> model {  
> // Hyperparameters  
> mu\_pr ~ normal(0, 1);  
> sigma ~ normal(0, 1);
> 
> // individual parameters  
> alpha\_pr ~ normal(0, 1);  
> beta\_pr ~ normal(0, 1);
> 
> for (i in 1:N) {  
> for (bIdx in 1:Bsubj[i]) { // new  
> // Define Values  
> vector[2] ev; // Expected value  
> real PE; // Prediction error
> 
> ```
> // Initialize values
> ev = initV; // Initial ev values
> 
> for (t in 1:Tsubj[i, bIdx]) {
> // Softmax choice
> choice[i, bIdx, t] ~ categorical_logit(ev * beta[i]);
> 
> // Prediction Error
> PE = outcome[i, bIdx, t] - ev[choice[i, bIdx, t]];
> 
> // Update expected value of chosen stimulus
> ev[choice[i, bIdx, t]] += alpha[i] * PE;
>     
> }
> }
> 
> ```
> 
> }  
> }
> 
> generated quantities {  
> // For group level parameters  
> real\<lower=0, upper=1\> mu\_alpha;  
> real\<lower=0, upper=10\> mu\_beta;
> 
> // For log likelihood calculation  
> real log\_lik[N];
> 
> // For model regressors  
> real ev\_c[N, B, T]; // Expected value of the chosen option  
> real ev\_nc[N, B, T]; // Expected value of the non-chosen option  
> real pe[N, B, T]; // Prediction error
> 
> // For posterior predictive check  
> real y\_pred[N, B, T];
> 
> // Initialize all the variables to avoid NULL values  
> for (i in 1:N) {  
> for (b in 1:B) {  
> for (t in 1:T) {  
> ev\_c[i, b, t] = 0;  
> ev\_nc[i, b, t] = 0;  
> pe[i, b, t] = 0;
> 
> ```
> y_pred[i, b, t] = -1;
> }
> }
> 
> ```
> 
> }
> 
> mu\_alpha = Phi\_approx(mu\_pr[1]);  
> mu\_beta = Phi\_approx(mu\_pr[2]) \* 10;
> 
> { // local section, this saves time and space  
> for (i in 1:N) {
> 
> ```
> log_lik[i] = 0;
> 
> for (bIdx in 1:Bsubj[i]) { // new
> // Define values
> vector[2] ev; // Expected value
> real PE; // prediction error
> 
> // Initialize values
> ev = initV; // initial ev values
> 
> for (t in 1:Tsubj[i, bIdx]) {
> // Softmax choice
> log_lik[i] += categorical_logit_lpmf(choice[i, bIdx, t] | ev * beta[i]);
> 
> // generate posterior prediction for current trial
> y_pred[i, bIdx, t] = categorical_rng(softmax(ev * beta[i]));
> 
> // Prediction Error
> PE = outcome[i, bIdx, t] - ev[choice[i, bIdx, t]];
> 
> // Store values for model regressors
> ev_c[i, bIdx, t] = ev[choice[i, bIdx, t]];
> ev_nc[i, bIdx, t] = ev[1 - choice[i, bIdx, t]];
> pe[i, bIdx, t] = PE;
> 
> // Update expected value of chosen stimulus
> ev[choice[i, bIdx, t]] += alpha[i] * PE;
> }
> }
> }
> 
> ```
> 
> }  
> }

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [June 21, 2023, 2:51pm UTC](https://discourse.mc-stan.org/t/how-can-i-solve-the-problem-of-high-pareto-k-values/31774/2 "2023-06-21T14:51:20Z")

</div>

You are doing leave-one-subject-out cross-validation and at the same time have two parameters for each subject. Removing the all observations for one subject makes the posterior for those two parameters to change a lot. See [LOO package glossary — loo-glossary • loo](https://mc-stan.org/loo/reference/loo-glossary.html#pareto-k-estimates-1) and [Roaches cross-validation demo](https://avehtari.github.io/modelselection/roaches.html#4_Poisson_model_with_%E2%80%9Crandom_effects%E2%80%9D) Section 4. You could solve this by using integrated loo as in [Roaches cross-validation demo](https://avehtari.github.io/modelselection/roaches.html#5_Poisson_model_with_%E2%80%9Crandom_effects%E2%80%9D_and_integrated_LOO) Section 5. As you have two parameters per subject, it is a bit trickier to do the integration but someone did report getting it to work in Stan discourse (sorry, I don’t have time to search that thread)

---

<div class="post-metadata">

**Author:** ![asyy](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/asyy/32/16522_2.png) [@asyy](https://discourse.mc-stan.org/u/asyy)\
**Post date:** [June 21, 2023, 8:52pm UTC](https://discourse.mc-stan.org/t/how-can-i-solve-the-problem-of-high-pareto-k-values/31774/3 "2023-06-21T20:52:07Z")

</div>

Thank you for the reply. I will try your suggestion.  
While reading [Roaches cross-validation demo](https://avehtari.github.io/modelselection/roaches.html#5_Poisson_model_with_%E2%80%9Crandom_effects%E2%80%9D_and_integrated_LOO) Section 5, I was confused whether I can apply the integration to a model with more than two parameters. My other models have 3 to 4 parameters and they also have the high pareto k problem. Can I do the integration when a model has more than 2 parameters?

Before I received your answer, I found that someone who had similar problem solved this problem of high pareto k values by modifying the code to do leave-one trial of one participant-out cross validation. I tried it as below:

```plaintext

data {
  int<lower=1> N; // Number of subjects

  int<lower=1> B; // Maximum number of blocks across subjects
  int<lower=1> Bsubj[N]; // Number of blocks for each subject

  int<lower=0> T; // Maximum number of trials across subjects
  int<lower=0, upper=T> Tsubj[N, B]; // Number of trials/blocks for each subject

  int<lower=-1, upper=2> choice[N, B, T]; // The choices subjects made
  real outcome[N, B, T]; // The outcome
}

transformed data {
  // Default value for (re-)initializing parameter vectors
  vector[2] initV;
  initV = rep_vector(0.0, 2);
}

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

  // Subject-level raw parameters (for Matt trick)
  vector[N] alpha_pr; // learning rate 
  vector[N] beta_pr; // inverse temperature
}

transformed parameters {
  // Transform subject-level raw parameters
  vector<lower=0, upper=1>[N] alpha;
  vector<lower=0, upper=10>[N] beta;

  for (i in 1:N) {
    alpha[i] = Phi_approx(mu_pr[1] + sigma[1] * alpha_pr[i]);
    beta[i] = Phi_approx(mu_pr[2] + sigma[2] * beta_pr[i]) * 10;
  }
}

model {
  // Hyperparameters
  mu_pr ~ normal(0, 1);
  sigma ~ normal(0, 1);

  // individual parameters
  alpha_pr ~ normal(0, 1);
  beta_pr ~ normal(0, 1);

  for (i in 1:N) {
    for (bIdx in 1:Bsubj[i]) { // new
      // Define Values
      vector[2] ev;// Expected value
      real PE; // Prediction error

      // Initialize values
      ev = initV; // Initial ev values

      for (t in 1:Tsubj[i, bIdx]) {
        // Softmax choice
        choice[i, bIdx, t] ~ categorical_logit(ev * beta[i]);

        // Prediction Error
        PE = outcome[i, bIdx, t] - ev[choice[i, bIdx, t]];

        // Update expected value of chosen stimulus
        ev[choice[i, bIdx, t]] += alpha[i] * PE;
        
      }
    }
  }
}

generated quantities {
  // For group level parameters
  real<lower=0, upper=1> mu_alpha;
  real<lower=0, upper=10> mu_beta;

  // For log likelihood calculation
  real log_lik[N, B, T];

  // For model regressors
  real ev_c[N, B, T]; // Expected value of the chosen option
  real ev_nc[N, B, T]; // Expected value of the non-chosen option
  real pe[N, B, T]; // Prediction error

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

  // Initialize all the variables to avoid NULL values
  for (i in 1:N) {
    for (b in 1:B) {
      for (t in 1:T) {
        ev_c[i, b, t] = 0;
        ev_nc[i, b, t] = 0;
        pe[i, b, t] = 0;

        y_pred[i, b, t] = -1;
        log_lik[i, b, t] = 0;
      }
    }
  }

  mu_alpha = Phi_approx(mu_pr[1]);
  mu_beta = Phi_approx(mu_pr[2]) * 10;

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

      for (bIdx in 1:Bsubj[i]) { // new
        // Define values
        vector[2] ev;// Expected value
        real PE;   

        // Initialize values
        ev = initV; // initial ev values

        for (t in 1:Tsubj[i, bIdx]) {
          // Softmax choice
          log_lik[i, bIdx, t] = categorical_logit_lpmf(choice[i, bIdx, t] | ev * beta[i]);

          // generate posterior prediction for current trial
          y_pred[i, bIdx, t] = categorical_rng(softmax(ev * beta[i]));

          // Prediction Error
          PE = outcome[i, bIdx, t] - ev[choice[i, bIdx, t]];

          // Store values for model regressors
          ev_c[i, bIdx, t] = ev[choice[i, bIdx, t]];
          ev_nc[i, bIdx, t] = ev[1 - choice[i, bIdx, t]];
          pe[i, bIdx, t] = PE;

          // Update expected value of chosen stimulus
          ev[choice[i, bIdx, t]] += alpha[i] * PE;
        }
      }
    }
  }
}

```

Now the pareto k result is better than before, but still there are some very bad k values. And the warning says: Can’t fit generalized Pareto distribution because all tail values are the same.  
I cannot understand why this warning appears. Could you please help me if you know why this warning appears and how to solve this problem?

Thank you so much

```plaintext
Computed from 12000 by 4920 log-likelihood matrix

         Estimate SE
elpd_loo -1872.4 41.6
p_loo 51.6 2.4
looic 3744.8 83.2
------
Monte Carlo SE of elpd_loo is NA.

Pareto k diagnostic values:
                         Count Pct. Min. n_eff
(-Inf, 0.5] (good) 4756 96.7% 4149      
 (0.5, 0.7] (ok) 0 0.0% <NA>      
   (0.7, 1] (bad) 0 0.0% <NA>      
   (1, Inf) (very bad) 164 3.3% 6000      
See help('pareto-k-diagnostic') for details.

```

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [June 26, 2023, 8:01am UTC](https://discourse.mc-stan.org/t/how-can-i-solve-the-problem-of-high-pareto-k-values/31774/4 "2023-06-26T08:01:30Z")

</div>

> [@asyy](#):
>
> My other models have 3 to 4 parameters and they also have the high pareto k problem. Can I do the integration when a model has more than 2 parameters?

In theory yes, but in practice it gets very difficult as 1) the number of evaluation points in the quadrature grow exponentially with the number of dimensions and in this case it is likely that some adaptive quadrature is needed so it’s not easy to keep the number of quadrature points small, 2) the implementation in Stan using nested 1D quadrature would be cumbersome.

> [@asyy](#):
>
> solved this problem of high pareto k values by modifying the code to do leave-one trial of one participant-out cross validation

That is often more stable as you are leaving less data out, and if you don’t need leave-one-person-out, then go for it.

> [@asyy](#):
>
> And the warning says: Can’t fit generalized Pareto distribution because all tail values are the same.

Look at the trials corresponding to the high Pareto k’s and warnings, and it is likely that you see that those trials are somehow special or some persons have only one trial.