translationalneuromodeling / tapas

TAPAS - Translational Algorithms for Psychiatry-Advancing Science
https://translationalneuromodeling.github.io/tapas/
GNU General Public License v3.0
219 stars 90 forks source link

complex values in HGF output #248

Open Heechberri opened 1 year ago

Heechberri commented 1 year ago

Hi HGF experts,

I am a newbie in coding and computational neuroscience, so please be patient with me ;D

I am using tapas_hgf_whatworld, tapas_logrt_linear_whatworld (edited, edits explained below), and tapas_quasinewton_optim to model the Bayesian inference parameters the reaction time from a pSRTT task (for experiment details see below, adapted from Marshall et al, 2016).

About half of the model outputs contain complex numbers. Here is a sample of the output:

Ignored trials: none Irregular trials: 570 Warning: optimization terminated because the maximum number of resets was reached. Parameter estimates for the perceptual model: mu2_0: [1×16 double] sa2_0: [1×16 double] mu3_0: 1 sa3_0: 0.1000 ka: 0.7311 om: -3.2273 + 0.0262i th: 0 Parameter estimates for the observation model: be0: 2.7675 - 0.0042i be1: 0.2245 - 0.0025i be2: 0.3921 - 0.0052i be3: 0.4723 + 0.0962i be4: 0.2618 - 0.0041i be5: 0.0103 - 0.0000i ze: 0.2518 + 0.0098i Model quality: LME (more is better): -178323.7548+332651.3928i AIC (less is better): 356671.2263+665302.0195i BIC (less is better): 356706.9056+665302.0195i AIC and BIC are approximations to -2*LME = 356647.5096-665302.7856i.

Changes made to the response model logrt_linear_whatworld (attached).

1. I am using original RT (milliseconds/100) values instead of log RT. RT in milliseconds is divided by 100 to scale the values of the reaction time down closer to the scale of the Bayesian parameters which are around the scan of 10 to the power of 1 to -2

2. I changed the response variables to RT~ be0 + be1.prediction error at level 1 + be2.precision weights at level 2 + be3.precision weighted prediction error at level 3, and be4.mean of level 3 + be5.post error slowing + ze.

There are two things I couldn't understand about the original script,

1. Why were the trajectories not used directly in the response model instead of the infStates?
I assumed that Marshall et al, 2016 used the linear_logrt_whatworld scripts, however, the variables are different from what is reported in Figure 3 of the paper. From Figure, I assumed that the variable of interest is first-level PE, third-level precision weighted prediction errors, and third-level mean and post error slowing. These perceptual parameters can already be found in the traj variable. I tried to use actual traj variables (attached as whatworld3 family of scripts) and equations from mathys et al. 2011;14 to derive the variables I needed (attached as whatworld2 family of scripts) and both produced complex numbers.

From my understanding, the first equation in linear_logrt_whatworld is calculating Shannon's surprise, i.e. -log(probability seeing that stimulus), but I don't quite understand the other two equations. which brings me to my next question

2. Why were the calculated response variables "transformed to 1st level" and how are they transformed to the first level?

I am well aware that my edits could have caused the weird problems, and I have spent quite some time trying to troubleshoot myself to no avail... Thank you for all the patience and understanding and help!

Vae

tapas_rt_linear_whatworld2.txt tapas_rt_linear_whatworld3.txt tapas_hgf_whatworld2.txt tapas_hgf_whatworld3.txt tapas_hgf_whatworld3_config.txt tapas_rt_linear_whatworld3_config.txt tapas_rt_linear_whatworld2_config.txt tapas_hgf_whatworld2_config.txt

I also turned on verbose (in tapas_quasinewton_optim_config) for one of the subjects just to see what is going on and here is the output if it helps. After about 9 iterations the improvements were really very small (less than 1), I am not sure why the program is still optimizing.

Initial argument: -3.5 2.5 0 0 0 0 0 0.5 Initial value: 676.1826

Argument: -3.5 2.573 0.1162 0.17898 0.0072271 0.071485 0.036284 0.57185 Value: 664.6223 Improvement: 11.5602 Regularizations: 2

Argument: -3.5354 2.1268 0.54092 0.53861 0.060322 -0.36061 0.3556 1.0173 Value: 656.6636 Improvement: 7.9587 Regularizations: 0

Argument: -3.6502 1.8161 1.1133 0.8078 0.12187 -0.65904 0.71337 0.49836 Value: 612.4991 Improvement: 44.1645 Regularizations: 0

Argument: -3.9365 1.7706 1.453 0.71613 0.15424 -0.69966 0.88412 0.39641 Value: 599.1132 Improvement: 13.3859 Regularizations: 1

Argument: -3.7525 1.7894 1.8366 0.46892 0.19065 -0.66942 0.9566 0.393 Value: 595.5051 Improvement: 3.6081 Regularizations: 1

Argument: -3.9411 1.7534 2.1248 0.38109 0.24591 -0.69908 0.025114 0.42057 Value: 593.787 Improvement: 1.7181 Regularizations: 0

Argument: -4.0697 2.0909 1.5708 0.033141 0.23601 -0.35894 0.57925 0.55758 Value: 564.7664 Improvement: 29.0205 Regularizations: 0

Argument: -3.8629 2.2682 1.2781 0.087939 0.013656 -0.41633 0.4871 0.3621 Value: 544.3717 Improvement: 20.3947 Regularizations: 0

Argument: -3.7904 1.8708 1.2785 0.15775 0.87632 -0.14187 0.59426 0.39276 Value: 540.9362 Improvement: 3.4355 Regularizations: 0

Argument: -3.7311 2.0964 1.4282 0.16457 -0.018872 -0.48497 0.65426 0.4135 Value: 539.8702 Improvement: 1.066 Regularizations: 0

Argument: -3.5995 2.1153 1.4102 0.09726 0.17489 -0.43349 0.65159 0.3657 Value: 536.6007 Improvement: 3.2695 Regularizations: 0

Argument: -3.6175 2.0196 1.3826 0.10298 0.1452 -0.31318 0.62911 0.37312 Value: 536.3899 Improvement: 0.21085 Regularizations: 1

Argument: -3.6197 2.099 1.3683 0.10514 0.20577 -0.38499 0.62214 0.37675 Value: 536.3502 Improvement: 0.039706 Regularizations: 3

Argument: -3.5744 2.0976 1.3607 0.11437 0.17808 -0.38471 0.61591 0.38122 Value: 536.1995 Improvement: 0.15066 Regularizations: 0

Argument: -3.3981 2.1026 1.3382 0.10022 0.14635 -0.35944 0.61155 0.38873 Value: 535.9581 Improvement: 0.24145 Regularizations: 0

Argument: -3.2844 2.1042 1.303 0.11465 0.14947 -0.35198 0.61126 0.38404 Value: 535.7392 Improvement: 0.21885 Regularizations: 0

Argument: -3.2012 2.1102 1.3029 0.091889 0.15758 -0.33967 0.61899 0.38211 Value: 535.5285 Improvement: 0.2107 Regularizations: 0

Argument: -3.1274 2.1151 1.2841 0.089102 0.15962 -0.33197 0.62214 0.38135 Value: 535.4554 Improvement: 0.073162 Regularizations: 0

Argument: -2.9519 2.1229 1.2521 0.078532 0.15698 -0.31251 0.62519 0.38072 Value: 535.281 Improvement: 0.17436 Regularizations: 0

Argument: -2.906 2.1238 1.2367 0.079883 0.15553 -0.30579 0.62484 0.381 Value: 535.2619 Improvement: 0.019114 Regularizations: 2

Argument: -2.9256 2.1221 1.2318 0.085255 0.15439 -0.30566 0.62263 0.38179 Value: 535.2587 Improvement: 0.0031995 Regularizations: 1

Regularizations exhausted - resetting algorithm. Initial argument: -2.9831 2.1599 1.1086 0.07673 0.13895 -0.2751 0.56037 0.39361 Initial value: 537.6173

Argument: -2.9816 2.1635 1.1147 0.089637 0.13928 -0.27155 0.56207 0.39046 Value: 537.0565 Improvement: 0.56079 Regularizations: 6

Argument: -2.7648 1.69 1.5178 0.43112 0.18197 -0.7306 0.82792 -0.017462 Value: -10533.1412+895.353906i Improvement: 11070.1977-895.353906i Regularizations: 0

Argument: -2.5774-0.00011243i 1.5976+0.00051045i 1.5093-7.9624e-05i 0.53169+0.00023444i 0.18416-2.053e-05i -0.82141+0.00049676i 0.84283-0.00014167i -0.026891+0.0096193i Value: -155274.9814+54657.77116i Improvement: 144741.8402-53762.41726i Regularizations: 2

Regularizations exhausted - resetting algorithm. Initial argument: -2.6696-0.00010119i 1.6878+0.0004594i 1.3583-7.1661e-05i 0.47852+0.000211i 0.16575-1.8477e-05i -0.73927+0.00044709i 0.75854-0.0001275i 0.025798+0.0086573i Initial value: 68838.5032+23184.7368i

Regularizations exhausted - resetting algorithm. Initial argument: -2.7527-9.1071e-05i 1.7691+0.00041346i 1.2225-6.4495e-05i 0.43067+0.0001899i 0.14917-1.6629e-05i -0.66534+0.00040238i 0.68269-0.00011475i 0.073218+0.0077916i Initial value: 2043.8548+209.61148i

Argument: -3.1527-0.042494i 1.8677+0.011611i 1.2416+0.0025151i 0.18476-0.024652i 0.14718-0.00020585i -0.56746+0.011516i 0.67167-0.0012071i 0.028989-0.063329i Value: 639.66078-1598.8741i Improvement: 1404.1941+1808.4856i Regularizations: 1

Argument: -3.3042-0.052972i 1.9482-0.065636i 1.2713-0.041557i 0.13543-0.11833i 0.14676-0.0011506i -0.48832-0.063752i 0.66753-0.0024727i -0.0021834-0.07045i Value: -901.85163-949.85631i Improvement: 1541.5124-649.01778i Regularizations: 2

Argument: -3.3146-0.088994i 1.9684-0.073807i 1.2616-0.042556i 0.11085-0.10184i 0.14528-0.001165i -0.46848-0.071709i 0.66188-0.00073914i -0.039239+0.033997i Value: -1207.5107+1236.645i Improvement: 305.65908-2186.5013i Regularizations: 3

Argument: -3.3208-0.10383i 1.9653-0.084521i 1.2566-0.03293i 0.11791-0.088846i 0.14479-7.5683e-05i -0.47152-0.08237i 0.66113+0.0030756i -0.050065+0.033993i Value: -1227.9977+832.09881i Improvement: 20.48693+404.5462i Regularizations: 5

Regularizations exhausted - resetting algorithm. Initial argument: -3.3388-0.093447i 2.0187-0.076069i 1.131-0.029637i 0.10612-0.079961i 0.13031-6.8115e-05i -0.42437-0.074133i 0.59502+0.002768i 0.0049413+0.030594i Initial value: 1339.797+2757.5232i

Argument: -3.3411-0.060095i 2.4287+0.22093i 1.3793+0.16562i 0.60643+0.3023i 0.13631+0.0057761i -0.024902+0.21531i 0.60585+0.017233i -0.01441-0.016616i Value: -33607.9946+7012.945i Improvement: 34947.7916-4255.42175i Regularizations: 0

Regularizations exhausted - resetting algorithm. Initial argument: -3.357-0.054085i 2.4359+0.19883i 1.2414+0.14905i 0.54579+0.27207i 0.12268+0.0051984i -0.022412+0.19378i 0.54526+0.01551i 0.037031-0.014955i Initial value: 1941.24079-14534.4361i

Regularizations exhausted - resetting algorithm. Initial argument: -3.3713-0.048677i 2.4423+0.17895i 1.1172+0.13415i 0.49121+0.24486i 0.11041+0.0046786i -0.020171+0.1744i 0.49073+0.013959i 0.083328-0.013459i Initial value: 2061.3925-5128.8501i

Regularizations exhausted - resetting algorithm. Initial argument: -3.3841-0.043809i 2.4481+0.16106i 1.0055+0.12073i 0.44209+0.22038i 0.099369+0.0042107i -0.018154+0.15696i 0.44166+0.012563i 0.12499-0.012113i Initial value: 1431.8639-2647.775i

Regularizations exhausted - resetting algorithm. Initial argument: -3.3957-0.039428i 2.4532+0.14495i 0.90495+0.10866i 0.39788+0.19834i 0.089432+0.0037897i -0.016338+0.14127i 0.3975+0.011307i 0.1625-0.010902i Initial value: 1084.7694-1583.8371i

Argument: -3.4415+0.079807i 2.0657+0.46417i 0.69772+0.29954i -0.06574+0.58144i 0.0862+0.0081696i -0.39375+0.45247i 0.39864+0.019002i 0.07065+0.19578i Value: -3925.8444-2022.2796i Improvement: 5010.6138+438.44257i Regularizations: 0

Regularizations exhausted - resetting algorithm. Initial argument: -3.4474+0.071826i 2.1092+0.41775i 0.62795+0.26959i -0.059166+0.52329i 0.07758+0.0073527i -0.35437+0.40722i 0.35877+0.017102i 0.11359+0.1762i Initial value: -3334.7892-710.2984i

Argument: -3.4386+0.075176i 2.0926+0.44022i 0.61843+0.28387i -0.077293+0.55179i 0.07738+0.0077368i -0.37048+0.42912i 0.35856+0.01804i 0.064059+0.075543i Value: -8459.1434-56.554157i Improvement: 5124.3542-653.74424i Regularizations: 3

Regularizations exhausted - resetting algorithm. Initial argument: -3.4448+0.067658i 2.1334+0.3962i 0.55659+0.25548i -0.069564+0.49661i 0.069642+0.0069631i -0.33343+0.38621i 0.3227+0.016236i 0.10765+0.067989i Initial value: -4873.2916+1750.8074i

Warning: optimization terminated because the maximum number of resets was reached.

chmathys commented 1 year ago

Hi Vae - just vey quickly, scanning your problem, I suspect that at some point you're taking the log (or square root) of a negative number, which gives you a complex result. This probably happens because you're using RT's instead of log-RT's. When you use a Gaussian model with RT's, negative values are possible according to that model, which doesn't make sense. The solution is to use either a log-Gaussian model with RT's or a Gaussian model with log-RT's

Heechberri commented 1 year ago

Hi @chmathys

Thanks for the reply!

I have generated another set of parameters using logrt, and all the parameters look fine :) Thank you for your help.

Another issue is that I would like to include the subjects' actual response (categorical) alongside the RT and stimulus train when using HGF. I came across a poster from your lab at CCN this year discussing this exact application. Is the script for this already available in tapas?

Lastly, I still can't figure out how to get the equations for expected uncertainty and unexpected uncertainty in the logrt_linear_whatworld response model, could you give me some hints? below are the equations that I am referring to:

% mu1 contains the actually occurring transition -> multiply with
% mu1hat to get probability of that transition (other elements are
% zero)
otp    = mu1.*mu1hat; % observed transition probabilities (3-dim)
otps3  = sum(otp, 3, 'omitnan');      % sum over 3rd dim
otps23 = sum(otps3, 2, 'omitnan');    % sum over 2nd dim

surp = -log(otps23);%according to shannon's formula  the information contained in a data set D is given by  -log P (D)
surp(r.irr) = [];

% Expected uncertainty
% ~~~~~~~~~~~~~~~~~~~~
euo    = mu1.*sa2;    % expected uncertainty of observed transition (3-dim)
euos3  = sum(euo, 3, 'omitnan');      % sum over 3rd dim
euos23 = sum(euos3, 2, 'omitnan');    % sum over 2nd dim

to     = mu1.*mu2;    % tendency of observed transition (3-dim)
tos3   = sum(to, 3, 'omitnan');       % sum over 3rd dim
tos23  = sum(tos3, 2, 'omitnan');     % sum over 2nd dim

eu = tapas_sgm(tos23,1).*(1-tapas_sgm(tos23,1)).*euos23; % transform down to 1st level
eu(r.irr) = [];

% Unexpected uncertainty
% ~~~~~~~~~~~~~~~~~~~~~~
ueu = tapas_sgm(tos23,1).*(1-tapas_sgm(tos23,1)).*exp(mu3); % transform down to 1st level
ueu(r.irr) = [];

Thank you so much!!!

Regards, Vae

chmathys commented 1 year ago

Hi Vae - @alexjhess might be able to help you with multimodal responses.

Regarding the uncertainty calculations, much of that is about finding the right values in multidimensional arrays. I suggest using the debugger to look at how this unfolds step by step. The only substantive thing that happens is that once we have the desired values (euos23 and exp(mu3)), we need to transform them down to the 1st level so they are on the same scale as the other predictors in the linear model. An explanation of this transformation is in the supplementary material to Iglesias et al. (2013) https://doi.org/10.1016/j.neuron.2013.09.009.

alexjhess commented 1 year ago

Hi @Heechberri,

All the functionality you need to fit your models to multiple response data modalities is already implemented in tapas. The most straightforward way to build a customized response model for your application is to create a new obs model where you essentially sum the log likelihood of separate response data modalities. There is no concrete example model implemented in the toolbox yet, maybe we can include one in the next release (@chmathys).

Anyways, we hope to be able to upload a preprint of the work you were referring to including example code and models by the end of this month. Feel free to contact me via hess@biomed.ee.ethz.ch in case of delays or if you're in desperate need of example code in the mean time. Hope this is helpful.

All the best, Alex

Heechberri commented 1 year ago

@chmathys

Thank you for referring me to Igelesias et al supplementary, it really helped me understand more of the transforms. Since $dx_1=\frac{dx_1}{dx_2}\cdot dx_2$ only applies to the relationship between $x_1$ and $x_2$ then $s(x)(1-s(x))dx_2$ would only apply to transforms for level 2 right? For level 3 transforms it will be $$dx_1=\frac{dx_1}{dx_2}\cdot\frac{dx_2}{dx_3}\cdot dx_3$$ and that will be $$s(x)(1-s(x))e^{(x_3)}dx_3$$ since the relationship between $x_2$ and $x_3$ is $e^{\kappa +\omega}$?

How about those variables that are calculated from two different levels? If I wanted to transform precision weight at level 3 (which is $\frac{\pi_2}{\pi_3}$) down to level 1, would I apply tapas_sgm(tos23,1).(1-tapas_sgm(tos23,1)).exp(mu3).*psi3 or do I have to do it separately for $\pi_2$ (transform from level 2 to 1) and $\pi_3$ (transform from level 3 to 2) and calculate precision weight at level 3 from these transformed values?

@alexjhess Thanks for your response! I see... do you mean editing fitmodel.m to 1) accept two response model, 2) sum logLl from both response models, and 3) sum negLogJoint from both logLI and their priors as well?

Here are the code sniplets farom fitmodel.m that I am talking about.

trialLogLls = obs_fun(r, infStates, ptrans_obs);
logLl = sum(trialLogLls, 'omitnan');
negLogLl = -logLl;
logObsPriors = -1/2.*log(8*atan(1).*r.c_obs.priorsas(obs_idx)) - 1/2.*(ptrans_obs(obs_idx) - r.c_obs.priormus(obs_idx)).^2./r.c_obs.priorsas(obs_idx);
logObsPrior  = sum(logObsPriors);

negLogJoint = -(logLl + logPrcPrior + logObsPrior);

Hahaha... as I am a novice coder, I am not confident in coding these by myself so I think I will still email you for the sample code soon!

Thank you sooo much for your help! looking forward to the publication!!

Regards, Vae

alexjhess commented 1 year ago

@Heechberri

No, I was thinking of something even simpler than modifying the fitModel.m function, namely creating a new response model that consists of a combination of several response models for different data stream. Your function could look something like this:

function [logp, yhat, res, logp_split] = comb_obs(r, infStates, ptrans)

  % part of obs model for binary choices
  [logp_bin, yhat_bin, res_bin] = tapas_unitsq_sgm(r, infStates, ptrans);

  % part of obs model for continuous response data modality (e.g. RTs)
  [logp_rt, yhat_rt, res_rt] = tapas_logrt_linear_whatworld(r, infStates, ptrans);

  % calculate log data likelihood (assuming independence)
  logp = logp_bin + logp_rt;

  yhat = [yhat_bin yhat_rt];
  res = [res_bin res_rt];
  logp_split = [logp_bin logp_rt];

end

Of course you would need to add the corresponding files '_config.m', '_transp.m', '_namep.m' and '_sim.m' that are used in the fitModel and simModel routines of the HGF Toolbox.

Hope this helps, otherwise shoot me an e-mail. :)

alexjhess commented 9 months ago

Hi @Heechberri,

I just wanted to let you know that we have uploaded a preprint introducing some example HGF response models that simultaneously model binary choices and continuous RTs. You can check it out at https://www.biorxiv.org/content/10.1101/2024.02.19.581001v1 and there are also links to code and data included.

Cheers!

Heechberri commented 9 months ago

Hi! Thanks! This is awesome!

Regards, Vae

On Fri, 23 Feb 2024 at 6:26 PM, alexjhess @.***> wrote:

Hi @Heechberri https://github.com/Heechberri,

I just wanted to let you know that we have uploaded a preprint introducing some example HGF response models that simultaneously model binary choices and continuous RTs. You can check it out at https://www.biorxiv.org/content/10.1101/2024.02.19.581001v1 and there are also links to code and data included.

Cheers!

— Reply to this email directly, view it on GitHub https://github.com/translationalneuromodeling/tapas/issues/248#issuecomment-1961074649, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALSMZKOKJPLPEZ6AAUWAMHLYVBVDZAVCNFSM6AAAAAA5MP3JHOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSNRRGA3TINRUHE . You are receiving this because you were mentioned.Message ID: @.***>