Warning message about the value of pareto k

[Juntao Ni]

Dear group,

First of all, thank you for creating this excellent package for us investigator in decision-making field. This help me a lot in my study.

Then, my questions is abou the warning message as following shows, the value of pareto k is too high, much more than 0.7. I don't know if I need to pay attention to it or I can ignore it. I have browsed the web and find the similar question, but I didn't fully understand. Do I need to fix it through adjusting some parameters?

Best wishes,

Ni

[Nathaniel Haines]

Hi Ni,

Thanks for your interest! 

So, the reason we get these warnings is related to how we are computing the log likelihood for each participant. Instead of saving the likelihood for each trial, we sum over the trials within participants, and then compute LOOIC on the resulting participant by MCMC sample matrix. The rationale is that this measure is an approximation to a leave one participant out score, rather than a leave one trial out score. 

The problem seems to come up because, for some tasks/models, participants very quite a bit on how well the model captures their choices, causing the importance sampling procedure used to estimate LOO to produce high pareto K diagnostics, which indicate that the LOOIC estimate may not be reliable. There are some details on this exact problem here: https://groups.google.com/g/stan-users/c/kew2-ODKYwQ/m/UBJ1JEiBAQAJ. Overall, the ideal solution is to perform an actual leave-one-participant out cross validation procedure, but this can be rather computationally intensive and we do not have a default functionality for it in the package at the moment. Another solution is to modify the stan code to keep the trial-level likelihood and compute LOOIC on the resulting matrix, but this changes the interpretation a bit. In my own experience, such a change does remedy the pareto k diagnostics, and I have found that both tend to identify the same model as "best", although I am sure this will depend on the model/task/data.

In my opinion, the best approach in this case is to treat the model comparison more comprehensively—e.g., do posterior predictive checks, parameter recovery analyses, or whatever else, and the LOOIC metric is then just another metric to use to help determine which model is best for your given problem. Ideally, you could also perform other types of cross validation checks. We are still thinking about a  more general solution, but hopefully the above can be helpful for now.

Best,

Nate 

[Juntao Ni]

Hi Nate,

Thank you for providing the solutions for me. As you said,  "Another solution is to modify the stan code to keep the trial-level likelihood and compute LOOIC on the resulting matrix...", but I’m not good at programming and not sure what to do about this solution. If there any paper, book, or manual advised for me to refer to and learn. 

Can I ignore the Pareto k warning, doing posterior predictive checks instead to make sure if the model is best?

Best,

Ni

[Haines, Nathaniel B.]

Hi Ni,

I do not have a particular source handy, but I have done this before in some work that has code available online. For example, lines 80-84 and line 110 here: https://github.com/CCS-Lab/Haines_2020_CPS/blob/master/Code/Stan/trait_explanatory_model_testset.stan. The extension is to set up a matrix/array that has rows for participants, and then columns for trials within participants. Then, instead of summing the likelihood over trials within participants, you can save the log likelihood for each trial separately. 

If you take the above approach, one the model is done fitting, you can extract the log likelihood array and then actually compute LOOIC for each participant and model independently, or you can calculate LOOIC on the whole array simultaneously, which is akin to something like leave-one-trial of one-participant out. 

If you do not go this route (or even if you do), I think doing the posterior predictions or other forms of simulation to critique the model would be best. Typically, statistical model fit is not the only thing worth considering, and using a holistic set of approaches may be best to answer your scientific question.

Sorry for the delay—let me know if this answers your question!

Best,

Nate

---

From: hbayesd...@googlegroups.com <hbayesd...@googlegroups.com> on behalf of Juntao Ni <juntao...@gmail.com>
Sent: Saturday, October 10, 2020 6:46 AM
To: hbayesdm-users <hbayesd...@googlegroups.com>
Subject: Re: Warning message about the value of pareto k

 

--
You received this message because you are subscribed to the Google Groups "hbayesdm-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hbayesdm-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/hbayesdm-users/101ac6bc-459d-46dd-91b3-5060aca3e981n%40googlegroups.com.

[Juntao Ni]

Hi Nate,

Thank you very much for providing me the code in detail. I will have a try on my own.

Here, I've got some new questions.

1.  I did the PPC, and the result shows as following, I'm not sure whether this indicated that the model was best fitted, and the task is the Iowa gambling task. I was thinking the reason that simulated data won't go up or down to the value 4 or 1, maybe because this was the average data? This is my inference, I don't know if it means any mistakes or not.
2. Another question is that there were three groups in my study and the results showed that the ORL model for IGT was the best in patient groups, while the VPP model was the best in the healthy group, so can I also use the parameter information of the ORL model in the healthy group to help me compare the characteristics of these three groups? I'm not sure how to make interpretations in this situation. 

Thanks again for answering my questions so kindly and patiently. 

Best wishes,

Ni

[Haines, Nathaniel B.]

Hi Ni,

Apologies for the long delay—for the IGT, I tend to prefer plotting the empirical versus posterior predicted probability of selecting each deck separately over trials. For example, you can access the posterior predictions through the hBayesDM object like so: 

fit$parVals$y_pred

This should be an array of dimension N_samples x N_subjects x N_trials, where each sample will be a value of 1, 2, 3, or 4, indicating a selection of that deck on the given trial for the given participant. 

You can use the R package bayesplot to nicely plot out the probability of selecting each deck across trials. For example, something like: 

library(dplyr)

y_pp <- fit$parVals$y_pred

y_true <- fit$raw_data %>%

     group_by(Trial) %>%

     summarize(prA = mean(choice==1, na.rm=T),

                        prB = mean(choice==2, na.rm=T),

                        prC = mean(choice==3, na.rm=T),

                        prD = mean(choice==4, na.rm=T))

should give you the posterior predictions and the proportion of subjects who selected each deck on a given trial. Note that you should make sure that the y_true$Trial value is actually ordered in the data.frame y_true by Trial number. I find sometimes that dplyr will set the order to something like 1, 10, etc., so just check that before proceeding with the next step. 

Then, you can do:

library(bayesplot)

plot_deckA <- ppc_intervals(

  y = y_true$prA,

  yrep = apply(y_pp, c(1,3), function(x) mean(x==1, na.rm=T)),

  x = y_true$Trial,

  prob = 0.5,

  prob_outer = .95

) +

  labs(

    x = "Trial",

    y = "Pr(choose A)",

    title = "Deck A"

  )

You can do the same to make a plot for deck B, C, and D, and then combine them like so: 

library(patchwork)

(plot_deckA | plot_deckB) /

    (plot_deckC | plot_deckD)  

This should give you a 2x2 grid, where the x-axis for each panel indicates the trial number, and y-axis the proportion of participants who selected the deck on each trial, along with the posterior predictions of the same values.

Sorry again for the delay, and let me know if this works out!

Best,

Nate 

---

From: hbayesd...@googlegroups.com <hbayesd...@googlegroups.com> on behalf of Juntao Ni <juntao...@gmail.com>

Sent: Wednesday, October 21, 2020 10:30 PM

To view this discussion on the web visit https://groups.google.com/d/msgid/hbayesdm-users/ebb41451-f3dc-41dd-af75-e465bfc3751en%40googlegroups.com.

[Juntao Ni]

Hi, Nate,

Thank you for replying to me with your detailed codes. I have tried with my output, but something seems wrong when I ran this code. I don't know how to solve it

Best,

Ni

[Juntao Ni]

It seems the image was not showing, so the error code shown as following

> y_true <- fit$raw_data %>%

+     group_by(trial) %>%

+     summarize(prA = mean(choice==1, na.rm=T),

+               prB = mean(choice==2, na.rm=T),

+               prC = mean(choice==3, na.rm=T),

+               prD = mean(choice==4, na.rm=T))

Error in UseMethod("group_by_") : 

  no applicable method for 'group_by_' applied to an object of class "NULL"

[Haines, Nathaniel B.]

Hi Ni,

It is hard to tell without the actual code and data, but perhaps the raw_data from "fit" is not coming through correctly? FYI—you will need to replace "fit" with the variable name assigned to that hBayesDM object. 

In case that isn't clear, could you post the code you used to fit the model in the first place? 

Best,

Nate 

---

From: hbayesd...@googlegroups.com <hbayesd...@googlegroups.com> on behalf of Juntao Ni <juntao...@gmail.com>

Sent: Thursday, December 10, 2020 3:13 AM

To view this discussion on the web visit https://groups.google.com/d/msgid/hbayesdm-users/9bd7c0cc-cf3a-4c6e-b5d2-ee0f0637fecfn%40googlegroups.com.

[Juntao Ni]

Hi Nate,

Sorry for the delay. This is my data. Thank you for your patient guidance.

My code

output <- igt_orl(

    data = SA100, niter = 4000, nwarmup = 1000, nchain = 4, ncore = 4, max_treedepth = 15, adapt_delta = 0.99, inc_postpred = TRUE)

Best，

Ni