crgagne / volatility_paper_elife

10 stars 6 forks source link

Errors of importing fitting data #1

Open lilihub opened 2 years ago

lilihub commented 2 years ago

Hi Crgagne,

This is Lili, a PhD student from Dublin City University. I'm really interested in the paper you published on eLief and am so impressed by the modeling analysis you made.

It was said in the left plot of figure 3 that the group mean learning rates are in logit space. I'm confused because I thought if it's in logit, the values should all be in the range 0 to 1. So I am trying to replicate your plots in figure 3. However, it gave me errors when I was trying to import the pickle files of the fitting results (please see the following error codes in the _Figures_3-4_Behavioral_Model_Exp1Results.ipynb file and the screenshot of the error information). I have followed the tutorial created the conda environment and installed the related packages. Could not figure out what was wrong on my side. Could you please help with this? Or at least could you please explain a bit the meaning of the logit space? I'm doing the modeling using Rstan apart from replicating your results with your pymc3 codes, so wanna make sure the Rstan models are correct by comparing the posterior distribution reported in your paper.

model_name = 'model=11_covariate=Bi3itemCDM_date=2020_4_27_samples=2000_seed=3_exp=1.pkl' with open('../fitting_behavioral_model/model_fits/'+model_name, "rb" ) as buff: model_output = pickle.load(buff)

image image

crgagne commented 2 years ago

Hi Lili,

Thanks for reaching out and for your interest in our work. To answer the question about the logits in Figure 3.

Sincerely, Chris

lilihub commented 2 years ago

Hi Chris,

Sorry for taking your first name wrong :) and thank you so much for your reply! I really appreciate your time explaining this.

Does it mean you used the posterior samples of the baseline and other learning rates directly without doing any transformation to plot them in Figure 3? In the modeling files, e.g. the following sentence in models_1.py under model_code folder, lr_t is the learning rate which was mapped to (0,1) using inverse-logit transform, but lr_baseline, lr_goodbas, etc were not converted to logit range, so they are actually in the R space. You took the posterior sample of these parameters directly and plotted them in Figure 3. Am I understanding it correctly? This question is important to me because I am trying to replicate your models in Rstan and testing them on a smaller sample. The range of the posterior distribution I obtained for the baseline learning rate is much larger. But in your results, all learning rates components are perfectly constrained in (0,1) or (-1,0), which makes me doubt the correctness of my implementation. _lr_t = lr_baseline + \ outcome_valence_tlr_goodbad + \ stabvol_tlr_stabvol + \ rewpain_tlr_rewpain + \ outcome_valence_tstabvol_tlr_goodbad_stabvol + \ outcome_valence_trewpain_tlr_rewpain_goodbad + \ stabvol_trewpain_tlr_rewpain_stabvol + \ outcome_valence_tstabvol_trewpain_tlr_rewpain_goodbad_stabvol lr_t = pm.invlogit(lrt)

Another question about model 1 is that it was mentioned in the paper that there were only three parameters in model 1 but in model 1.py, there are four extra parameters _lr_ct, _epst, _Amixt, _Bct and I could not find the risk parameter.

Thanks for the shared link. I will check it and see if I can figure the error out.

Very sorry that I have so many questions and bothering you again.

Best, Lili

lilihub commented 2 years ago

Hi Chris,

These are the group-level posterior distributions (without doing any transformation on the posterior samples) of the learning rate components I obtained for a much smaller sample. As you see, their ranges are much larger, so that's why I'm wondering if you did any transformation (didn't see any transformation in the plotting code though). image

Thanks again, Lili

crgagne commented 2 years ago

Hi Lili,

No problem at all for asking questions. I'm happy that you are interested in the models/code. In looking into your question about model_1.py, it looks like the file was misnamed in the version of the repo that you have. My sincerest apologies for that -- when I uploaded the repo to github for sharing, I changed the file names to be more intuitive and accidentally swapped the names for models_1.py and models2thr9_11.py. I adjusted this and re-pushed the repo, so the now current version is correct.

I'm running the model on my computer now and it seems to be running without issue. What's your email? I'll send you the environment.yml in case you want to try installing the conda environment that way. It won't upload here for some reason.

For fIgure 3a, I'm not doing any transformations. For figure 3b, I add the samples for the components up, using +1 -1 depending on the condition, with the formula:
lr_t = lr_baseline + outcome_valence_tlr_goodbad + stabvol_tlr_stabvol + rewpain_tlr_rewpain + ....

And then inverse-logit transform.

There's a few reasons that you might be seeing larger ranges of values than for you model. One might be because we also have separate components (e.g. lr_baseline, etc.) that get multiplied by the factor scores, so if you don't have these, then the basic components might have a larger range. The range is also probably highly depend on the number of participants -- with more participants, the range tends to be smaller for these group level parameters, because there is more evidence for the mean across participants (e.g. for lr_baseline).

lilihub commented 2 years ago

Hi Chris,

Thank you so much for your quick reply. It will be great if you can share the environment file. Thanks for updating the repo! my email: lili.zhang27@mail.dcu.ie

"lr_baseline multiplied by the factor score" refers to the following equation, right?

image

I do have sth similar to this in my model. Here A_base_cov corresponds to the weight parameters (beta_g, beta_d, etc) and cov corresponds to the clinical factors (not the same factors as in your paper, they are the clinical scores for the Parkinson subjects in our study instead. mu_pr corresponds to the mu in your equation.

image

These clinical factors are standardized before fitting to the model). I do agree with the second possible reason. I only have 20 participants tested on this model, which is much smaller than that of yours. I'm so happy to confirm with you about this. At least I got a bit more confidence about my implementation:).

Have a nice day!

Best, Lili