diff --git a/Python/hbayesdm/models/__init__.py b/Python/hbayesdm/models/__init__.py index 19d91f0b..b2cd71ec 100644 --- a/Python/hbayesdm/models/__init__.py +++ b/Python/hbayesdm/models/__init__.py @@ -40,10 +40,13 @@ from ._peer_ocu import peer_ocu from ._prl_ewa import prl_ewa from ._prl_fictitious import prl_fictitious +from ._prl_fictitious_cnc import prl_fictitious_cnc +from ._prl_fictitious_cnc_multipleB import prl_fictitious_cnc_multipleB from ._prl_fictitious_multipleB import prl_fictitious_multipleB from ._prl_fictitious_rp import prl_fictitious_rp from ._prl_fictitious_rp_woa import prl_fictitious_rp_woa from ._prl_fictitious_woa import prl_fictitious_woa +from ._prl_fictitious_woa_multipleB import prl_fictitious_woa_multipleB from ._prl_rp import prl_rp from ._prl_rp_multipleB import prl_rp_multipleB from ._pstRT_ddm import pstRT_ddm @@ -106,10 +109,13 @@ 'peer_ocu', 'prl_ewa', 'prl_fictitious', + 'prl_fictitious_cnc', + 'prl_fictitious_cnc_multipleB', 'prl_fictitious_multipleB', 'prl_fictitious_rp', 'prl_fictitious_rp_woa', 'prl_fictitious_woa', + 'prl_fictitious_woa_multipleB', 'prl_rp', 'prl_rp_multipleB', 'pstRT_ddm', diff --git a/Python/hbayesdm/models/_prl_fictitious_cnc.py b/Python/hbayesdm/models/_prl_fictitious_cnc.py new file mode 100644 index 00000000..1163b73d --- /dev/null +++ b/Python/hbayesdm/models/_prl_fictitious_cnc.py @@ -0,0 +1,244 @@ +from typing import Sequence, Union, Any +from collections import OrderedDict + +from numpy import Inf, exp +import pandas as pd + +from hbayesdm.base import TaskModel +from hbayesdm.preprocess_funcs import prl_preprocess_func + +__all__ = ['prl_fictitious_cnc'] + + +class PrlFictitiousCnc(TaskModel): + def __init__(self, **kwargs): + super().__init__( + task_name='prl', + model_name='fictitious_cnc', + model_type='', + data_columns=( + 'subjID', + 'choice', + 'outcome', + ), + parameters=OrderedDict([ + ('eta_c', (0, 0.5, 1)), + ('eta_nc', (0, 0.5, 1)), + ('alpha', (-Inf, 0, Inf)), + ('beta', (0, 1, 10)), + ]), + regressors=OrderedDict([ + ('ev_c', 2), + ('ev_nc', 2), + ('pe_c', 2), + ('pe_nc', 2), + ('dv', 2), + ]), + postpreds=['y_pred'], + parameters_desc=OrderedDict([ + ('eta_c', 'learning rate for chosen option'), + ('eta_nc', 'learning rate for non-chosen option'), + ('alpha', 'indecision point'), + ('beta', 'inverse temperature'), + ]), + additional_args=OrderedDict([ + + ]), + additional_args_desc=OrderedDict([ + + ]), + **kwargs, + ) + + _preprocess_func = prl_preprocess_func + + +def prl_fictitious_cnc( + data: Union[pd.DataFrame, str, None] = None, + niter: int = 4000, + nwarmup: int = 1000, + nchain: int = 4, + ncore: int = 1, + nthin: int = 1, + inits: Union[str, Sequence[float]] = 'vb', + ind_pars: str = 'mean', + model_regressor: bool = False, + vb: bool = False, + inc_postpred: bool = False, + adapt_delta: float = 0.95, + stepsize: float = 1, + max_treedepth: int = 10, + **additional_args: Any) -> TaskModel: + """Probabilistic Reversal Learning Task - Fictitious Update Model, with chosen/not-chosen learning rates + + Hierarchical Bayesian Modeling of the Probabilistic Reversal Learning Task + using Fictitious Update Model, with chosen/not-chosen learning rates [Glascher2009]_ with the following parameters: + "eta_c" (learning rate for chosen option), "eta_nc" (learning rate for non-chosen option), "alpha" (indecision point), "beta" (inverse temperature). + + + + + .. [Glascher2009] Glascher, J., Hampton, A. N., & O'Doherty, J. P. (2009). Determining a Role for Ventromedial Prefrontal Cortex in Encoding Action-Based Value Signals During Reward-Related Decision Making. Cerebral Cortex, 19(2), 483-495. https://doi.org/10.1093/cercor/bhn098 + + .. codeauthor:: Jinwoo Jeong + + User data should contain the behavioral data-set of all subjects of interest for + the current analysis. When loading from a file, the datafile should be a + **tab-delimited** text file, whose rows represent trial-by-trial observations + and columns represent variables. + + For the Probabilistic Reversal Learning Task, there should be 3 columns of data + with the labels "subjID", "choice", "outcome". It is not necessary for the columns to be + in this particular order; however, it is necessary that they be labeled + correctly and contain the information below: + + - "subjID": A unique identifier for each subject in the data-set. + - "choice": Integer value representing the option chosen on that trial: 1 or 2. + - "outcome": Integer value representing the outcome of that trial (where reward == 1, and loss == -1). + + .. note:: + User data may contain other columns of data (e.g. ``ReactionTime``, + ``trial_number``, etc.), but only the data within the column names listed + above will be used during the modeling. As long as the necessary columns + mentioned above are present and labeled correctly, there is no need to + remove other miscellaneous data columns. + + .. note:: + + ``adapt_delta``, ``stepsize``, and ``max_treedepth`` are advanced options that + give the user more control over Stan's MCMC sampler. It is recommended that + only advanced users change the default values, as alterations can profoundly + change the sampler's behavior. See [Hoffman2014]_ for more information on the + sampler control parameters. One can also refer to 'Section 34.2. HMC Algorithm + Parameters' of the `Stan User's Guide and Reference Manual`__. + + .. [Hoffman2014] + Hoffman, M. D., & Gelman, A. (2014). + The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. + Journal of Machine Learning Research, 15(1), 1593-1623. + + __ https://mc-stan.org/users/documentation/ + + Parameters + ---------- + data + Data to be modeled. It should be given as a Pandas DataFrame object, + a filepath for a data file, or ``"example"`` for example data. + Data columns should be labeled as: "subjID", "choice", "outcome". + niter + Number of iterations, including warm-up. Defaults to 4000. + nwarmup + Number of iterations used for warm-up only. Defaults to 1000. + + ``nwarmup`` is a numerical value that specifies how many MCMC samples + should not be stored upon the beginning of each chain. For those + familiar with Bayesian methods, this is equivalent to burn-in samples. + Due to the nature of the MCMC algorithm, initial values (i.e., where the + sampling chains begin) can have a heavy influence on the generated + posterior distributions. The ``nwarmup`` argument can be set to a + higher number in order to curb the effects that initial values have on + the resulting posteriors. + nchain + Number of Markov chains to run. Defaults to 4. + + ``nchain`` is a numerical value that specifies how many chains (i.e., + independent sampling sequences) should be used to draw samples from + the posterior distribution. Since the posteriors are generated from a + sampling process, it is good practice to run multiple chains to ensure + that a reasonably representative posterior is attained. When the + sampling is complete, it is possible to check the multiple chains for + convergence by running the following line of code: + + .. code:: python + + output.plot(type='trace') + ncore + Number of CPUs to be used for running. Defaults to 1. + nthin + Every ``nthin``-th sample will be used to generate the posterior + distribution. Defaults to 1. A higher number can be used when + auto-correlation within the MCMC sampling is high. + + ``nthin`` is a numerical value that specifies the "skipping" behavior + of the MCMC sampler. That is, only every ``nthin``-th sample is used to + generate posterior distributions. By default, ``nthin`` is equal to 1, + meaning that every sample is used to generate the posterior. + inits + String or list specifying how the initial values should be generated. + Options are ``'fixed'`` or ``'random'``, or your own initial values. + ind_pars + String specifying how to summarize the individual parameters. + Current options are: ``'mean'``, ``'median'``, or ``'mode'``. + model_regressor + Whether to export model-based regressors. For this model they are: "ev_c", "ev_nc", "pe_c", "pe_nc", "dv". + vb + Whether to use variational inference to approximately draw from a + posterior distribution. Defaults to ``False``. + inc_postpred + Include trial-level posterior predictive simulations in + model output (may greatly increase file size). Defaults to ``False``. + adapt_delta + Floating point value representing the target acceptance probability of a new + sample in the MCMC chain. Must be between 0 and 1. See note below. + stepsize + Integer value specifying the size of each leapfrog step that the MCMC sampler + can take on each new iteration. See note below. + max_treedepth + Integer value specifying how many leapfrog steps the MCMC sampler can take + on each new iteration. See note below. + **additional_args + Not used for this model. + + Returns + ------- + model_data + An ``hbayesdm.TaskModel`` instance with the following components: + + - ``model``: String value that is the name of the model ('prl_fictitious_cnc'). + - ``all_ind_pars``: Pandas DataFrame containing the summarized parameter values + (as specified by ``ind_pars``) for each subject. + - ``par_vals``: OrderedDict holding the posterior samples over different parameters. + - ``fit``: A PyStan StanFit object that contains the fitted Stan model. + - ``raw_data``: Pandas DataFrame containing the raw data used to fit the model, + as specified by the user. + - ``model_regressor``: Dict holding the extracted model-based regressors. + + Examples + -------- + + .. code:: python + + from hbayesdm import rhat, print_fit + from hbayesdm.models import prl_fictitious_cnc + + # Run the model and store results in "output" + output = prl_fictitious_cnc(data='example', niter=2000, nwarmup=1000, nchain=4, ncore=4) + + # Visually check convergence of the sampling chains (should look like "hairy caterpillars") + output.plot(type='trace') + + # Plot posterior distributions of the hyper-parameters (distributions should be unimodal) + output.plot() + + # Check Rhat values (all Rhat values should be less than or equal to 1.1) + rhat(output, less=1.1) + + # Show the LOOIC and WAIC model fit estimates + print_fit(output) + """ + return PrlFictitiousCnc( + data=data, + niter=niter, + nwarmup=nwarmup, + nchain=nchain, + ncore=ncore, + nthin=nthin, + inits=inits, + ind_pars=ind_pars, + model_regressor=model_regressor, + vb=vb, + inc_postpred=inc_postpred, + adapt_delta=adapt_delta, + stepsize=stepsize, + max_treedepth=max_treedepth, + **additional_args) diff --git a/Python/hbayesdm/models/_prl_fictitious_cnc_multipleB.py b/Python/hbayesdm/models/_prl_fictitious_cnc_multipleB.py new file mode 100644 index 00000000..2ec35997 --- /dev/null +++ b/Python/hbayesdm/models/_prl_fictitious_cnc_multipleB.py @@ -0,0 +1,246 @@ +from typing import Sequence, Union, Any +from collections import OrderedDict + +from numpy import Inf, exp +import pandas as pd + +from hbayesdm.base import TaskModel +from hbayesdm.preprocess_funcs import prl_multipleB_preprocess_func + +__all__ = ['prl_fictitious_cnc_multipleB'] + + +class PrlFictitiousCncMultipleb(TaskModel): + def __init__(self, **kwargs): + super().__init__( + task_name='prl', + model_name='fictitious_cnc', + model_type='multipleB', + data_columns=( + 'subjID', + 'block', + 'choice', + 'outcome', + ), + parameters=OrderedDict([ + ('eta_c', (0, 0.5, 1)), + ('eta_nc', (0, 0.5, 1)), + ('alpha', (-Inf, 0, Inf)), + ('beta', (0, 1, 10)), + ]), + regressors=OrderedDict([ + ('ev_c', 3), + ('ev_nc', 3), + ('pe_c', 3), + ('pe_nc', 3), + ('dv', 3), + ]), + postpreds=['y_pred'], + parameters_desc=OrderedDict([ + ('eta_c', 'learning rate for chosen option'), + ('eta_nc', 'learning rate for non-chosen option'), + ('alpha', 'indecision point'), + ('beta', 'inverse temperature'), + ]), + additional_args=OrderedDict([ + + ]), + additional_args_desc=OrderedDict([ + + ]), + **kwargs, + ) + + _preprocess_func = prl_multipleB_preprocess_func + + +def prl_fictitious_cnc_multipleB( + data: Union[pd.DataFrame, str, None] = None, + niter: int = 4000, + nwarmup: int = 1000, + nchain: int = 4, + ncore: int = 1, + nthin: int = 1, + inits: Union[str, Sequence[float]] = 'vb', + ind_pars: str = 'mean', + model_regressor: bool = False, + vb: bool = False, + inc_postpred: bool = False, + adapt_delta: float = 0.95, + stepsize: float = 1, + max_treedepth: int = 10, + **additional_args: Any) -> TaskModel: + """Probabilistic Reversal Learning Task - Fictitious Update Model, with chosen/not-chosen learning rates + + Multiple-Block Hierarchical Bayesian Modeling of the Probabilistic Reversal Learning Task + using Fictitious Update Model, with chosen/not-chosen learning rates [Glascher2009]_ with the following parameters: + "eta_c" (learning rate for chosen option), "eta_nc" (learning rate for non-chosen option), "alpha" (indecision point), "beta" (inverse temperature). + + + + + .. [Glascher2009] Glascher, J., Hampton, A. N., & O'Doherty, J. P. (2009). Determining a Role for Ventromedial Prefrontal Cortex in Encoding Action-Based Value Signals During Reward-Related Decision Making. Cerebral Cortex, 19(2), 483-495. https://doi.org/10.1093/cercor/bhn098 + + .. codeauthor:: Jinwoo Jeong + + User data should contain the behavioral data-set of all subjects of interest for + the current analysis. When loading from a file, the datafile should be a + **tab-delimited** text file, whose rows represent trial-by-trial observations + and columns represent variables. + + For the Probabilistic Reversal Learning Task, there should be 4 columns of data + with the labels "subjID", "block", "choice", "outcome". It is not necessary for the columns to be + in this particular order; however, it is necessary that they be labeled + correctly and contain the information below: + + - "subjID": A unique identifier for each subject in the data-set. + - "block": A unique identifier for each of the multiple blocks within each subject. + - "choice": Integer value representing the option chosen on that trial: 1 or 2. + - "outcome": Integer value representing the outcome of that trial (where reward == 1, and loss == -1). + + .. note:: + User data may contain other columns of data (e.g. ``ReactionTime``, + ``trial_number``, etc.), but only the data within the column names listed + above will be used during the modeling. As long as the necessary columns + mentioned above are present and labeled correctly, there is no need to + remove other miscellaneous data columns. + + .. note:: + + ``adapt_delta``, ``stepsize``, and ``max_treedepth`` are advanced options that + give the user more control over Stan's MCMC sampler. It is recommended that + only advanced users change the default values, as alterations can profoundly + change the sampler's behavior. See [Hoffman2014]_ for more information on the + sampler control parameters. One can also refer to 'Section 34.2. HMC Algorithm + Parameters' of the `Stan User's Guide and Reference Manual`__. + + .. [Hoffman2014] + Hoffman, M. D., & Gelman, A. (2014). + The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. + Journal of Machine Learning Research, 15(1), 1593-1623. + + __ https://mc-stan.org/users/documentation/ + + Parameters + ---------- + data + Data to be modeled. It should be given as a Pandas DataFrame object, + a filepath for a data file, or ``"example"`` for example data. + Data columns should be labeled as: "subjID", "block", "choice", "outcome". + niter + Number of iterations, including warm-up. Defaults to 4000. + nwarmup + Number of iterations used for warm-up only. Defaults to 1000. + + ``nwarmup`` is a numerical value that specifies how many MCMC samples + should not be stored upon the beginning of each chain. For those + familiar with Bayesian methods, this is equivalent to burn-in samples. + Due to the nature of the MCMC algorithm, initial values (i.e., where the + sampling chains begin) can have a heavy influence on the generated + posterior distributions. The ``nwarmup`` argument can be set to a + higher number in order to curb the effects that initial values have on + the resulting posteriors. + nchain + Number of Markov chains to run. Defaults to 4. + + ``nchain`` is a numerical value that specifies how many chains (i.e., + independent sampling sequences) should be used to draw samples from + the posterior distribution. Since the posteriors are generated from a + sampling process, it is good practice to run multiple chains to ensure + that a reasonably representative posterior is attained. When the + sampling is complete, it is possible to check the multiple chains for + convergence by running the following line of code: + + .. code:: python + + output.plot(type='trace') + ncore + Number of CPUs to be used for running. Defaults to 1. + nthin + Every ``nthin``-th sample will be used to generate the posterior + distribution. Defaults to 1. A higher number can be used when + auto-correlation within the MCMC sampling is high. + + ``nthin`` is a numerical value that specifies the "skipping" behavior + of the MCMC sampler. That is, only every ``nthin``-th sample is used to + generate posterior distributions. By default, ``nthin`` is equal to 1, + meaning that every sample is used to generate the posterior. + inits + String or list specifying how the initial values should be generated. + Options are ``'fixed'`` or ``'random'``, or your own initial values. + ind_pars + String specifying how to summarize the individual parameters. + Current options are: ``'mean'``, ``'median'``, or ``'mode'``. + model_regressor + Whether to export model-based regressors. For this model they are: "ev_c", "ev_nc", "pe_c", "pe_nc", "dv". + vb + Whether to use variational inference to approximately draw from a + posterior distribution. Defaults to ``False``. + inc_postpred + Include trial-level posterior predictive simulations in + model output (may greatly increase file size). Defaults to ``False``. + adapt_delta + Floating point value representing the target acceptance probability of a new + sample in the MCMC chain. Must be between 0 and 1. See note below. + stepsize + Integer value specifying the size of each leapfrog step that the MCMC sampler + can take on each new iteration. See note below. + max_treedepth + Integer value specifying how many leapfrog steps the MCMC sampler can take + on each new iteration. See note below. + **additional_args + Not used for this model. + + Returns + ------- + model_data + An ``hbayesdm.TaskModel`` instance with the following components: + + - ``model``: String value that is the name of the model ('prl_fictitious_cnc_multipleB'). + - ``all_ind_pars``: Pandas DataFrame containing the summarized parameter values + (as specified by ``ind_pars``) for each subject. + - ``par_vals``: OrderedDict holding the posterior samples over different parameters. + - ``fit``: A PyStan StanFit object that contains the fitted Stan model. + - ``raw_data``: Pandas DataFrame containing the raw data used to fit the model, + as specified by the user. + - ``model_regressor``: Dict holding the extracted model-based regressors. + + Examples + -------- + + .. code:: python + + from hbayesdm import rhat, print_fit + from hbayesdm.models import prl_fictitious_cnc_multipleB + + # Run the model and store results in "output" + output = prl_fictitious_cnc_multipleB(data='example', niter=2000, nwarmup=1000, nchain=4, ncore=4) + + # Visually check convergence of the sampling chains (should look like "hairy caterpillars") + output.plot(type='trace') + + # Plot posterior distributions of the hyper-parameters (distributions should be unimodal) + output.plot() + + # Check Rhat values (all Rhat values should be less than or equal to 1.1) + rhat(output, less=1.1) + + # Show the LOOIC and WAIC model fit estimates + print_fit(output) + """ + return PrlFictitiousCncMultipleb( + data=data, + niter=niter, + nwarmup=nwarmup, + nchain=nchain, + ncore=ncore, + nthin=nthin, + inits=inits, + ind_pars=ind_pars, + model_regressor=model_regressor, + vb=vb, + inc_postpred=inc_postpred, + adapt_delta=adapt_delta, + stepsize=stepsize, + max_treedepth=max_treedepth, + **additional_args) diff --git a/Python/hbayesdm/models/_prl_fictitious_woa_multipleB.py b/Python/hbayesdm/models/_prl_fictitious_woa_multipleB.py new file mode 100644 index 00000000..f4ebd40b --- /dev/null +++ b/Python/hbayesdm/models/_prl_fictitious_woa_multipleB.py @@ -0,0 +1,242 @@ +from typing import Sequence, Union, Any +from collections import OrderedDict + +from numpy import Inf, exp +import pandas as pd + +from hbayesdm.base import TaskModel +from hbayesdm.preprocess_funcs import prl_multipleB_preprocess_func + +__all__ = ['prl_fictitious_woa_multipleB'] + + +class PrlFictitiousWoaMultipleb(TaskModel): + def __init__(self, **kwargs): + super().__init__( + task_name='prl', + model_name='fictitious_woa', + model_type='multipleB', + data_columns=( + 'subjID', + 'block', + 'choice', + 'outcome', + ), + parameters=OrderedDict([ + ('eta', (0, 0.5, 1)), + ('beta', (0, 1, 10)), + ]), + regressors=OrderedDict([ + ('ev_c', 3), + ('ev_nc', 3), + ('pe_c', 3), + ('pe_nc', 3), + ('dv', 3), + ]), + postpreds=['y_pred'], + parameters_desc=OrderedDict([ + ('eta', 'learning rate'), + ('beta', 'inverse temperature'), + ]), + additional_args=OrderedDict([ + + ]), + additional_args_desc=OrderedDict([ + + ]), + **kwargs, + ) + + _preprocess_func = prl_multipleB_preprocess_func + + +def prl_fictitious_woa_multipleB( + data: Union[pd.DataFrame, str, None] = None, + niter: int = 4000, + nwarmup: int = 1000, + nchain: int = 4, + ncore: int = 1, + nthin: int = 1, + inits: Union[str, Sequence[float]] = 'vb', + ind_pars: str = 'mean', + model_regressor: bool = False, + vb: bool = False, + inc_postpred: bool = False, + adapt_delta: float = 0.95, + stepsize: float = 1, + max_treedepth: int = 10, + **additional_args: Any) -> TaskModel: + """Probabilistic Reversal Learning Task - Fictitious Update Model, without alpha (indecision point) + + Multiple-Block Hierarchical Bayesian Modeling of the Probabilistic Reversal Learning Task + using Fictitious Update Model, without alpha (indecision point) [Glascher2009]_ with the following parameters: + "eta" (learning rate), "beta" (inverse temperature). + + + + + .. [Glascher2009] Glascher, J., Hampton, A. N., & O'Doherty, J. P. (2009). Determining a Role for Ventromedial Prefrontal Cortex in Encoding Action-Based Value Signals During Reward-Related Decision Making. Cerebral Cortex, 19(2), 483-495. https://doi.org/10.1093/cercor/bhn098 + + .. codeauthor:: Jinwoo Jeong + + User data should contain the behavioral data-set of all subjects of interest for + the current analysis. When loading from a file, the datafile should be a + **tab-delimited** text file, whose rows represent trial-by-trial observations + and columns represent variables. + + For the Probabilistic Reversal Learning Task, there should be 4 columns of data + with the labels "subjID", "block", "choice", "outcome". It is not necessary for the columns to be + in this particular order; however, it is necessary that they be labeled + correctly and contain the information below: + + - "subjID": A unique identifier for each subject in the data-set. + - "block": A unique identifier for each of the multiple blocks within each subject. + - "choice": Integer value representing the option chosen on that trial: 1 or 2. + - "outcome": Integer value representing the outcome of that trial (where reward == 1, and loss == -1). + + .. note:: + User data may contain other columns of data (e.g. ``ReactionTime``, + ``trial_number``, etc.), but only the data within the column names listed + above will be used during the modeling. As long as the necessary columns + mentioned above are present and labeled correctly, there is no need to + remove other miscellaneous data columns. + + .. note:: + + ``adapt_delta``, ``stepsize``, and ``max_treedepth`` are advanced options that + give the user more control over Stan's MCMC sampler. It is recommended that + only advanced users change the default values, as alterations can profoundly + change the sampler's behavior. See [Hoffman2014]_ for more information on the + sampler control parameters. One can also refer to 'Section 34.2. HMC Algorithm + Parameters' of the `Stan User's Guide and Reference Manual`__. + + .. [Hoffman2014] + Hoffman, M. D., & Gelman, A. (2014). + The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. + Journal of Machine Learning Research, 15(1), 1593-1623. + + __ https://mc-stan.org/users/documentation/ + + Parameters + ---------- + data + Data to be modeled. It should be given as a Pandas DataFrame object, + a filepath for a data file, or ``"example"`` for example data. + Data columns should be labeled as: "subjID", "block", "choice", "outcome". + niter + Number of iterations, including warm-up. Defaults to 4000. + nwarmup + Number of iterations used for warm-up only. Defaults to 1000. + + ``nwarmup`` is a numerical value that specifies how many MCMC samples + should not be stored upon the beginning of each chain. For those + familiar with Bayesian methods, this is equivalent to burn-in samples. + Due to the nature of the MCMC algorithm, initial values (i.e., where the + sampling chains begin) can have a heavy influence on the generated + posterior distributions. The ``nwarmup`` argument can be set to a + higher number in order to curb the effects that initial values have on + the resulting posteriors. + nchain + Number of Markov chains to run. Defaults to 4. + + ``nchain`` is a numerical value that specifies how many chains (i.e., + independent sampling sequences) should be used to draw samples from + the posterior distribution. Since the posteriors are generated from a + sampling process, it is good practice to run multiple chains to ensure + that a reasonably representative posterior is attained. When the + sampling is complete, it is possible to check the multiple chains for + convergence by running the following line of code: + + .. code:: python + + output.plot(type='trace') + ncore + Number of CPUs to be used for running. Defaults to 1. + nthin + Every ``nthin``-th sample will be used to generate the posterior + distribution. Defaults to 1. A higher number can be used when + auto-correlation within the MCMC sampling is high. + + ``nthin`` is a numerical value that specifies the "skipping" behavior + of the MCMC sampler. That is, only every ``nthin``-th sample is used to + generate posterior distributions. By default, ``nthin`` is equal to 1, + meaning that every sample is used to generate the posterior. + inits + String or list specifying how the initial values should be generated. + Options are ``'fixed'`` or ``'random'``, or your own initial values. + ind_pars + String specifying how to summarize the individual parameters. + Current options are: ``'mean'``, ``'median'``, or ``'mode'``. + model_regressor + Whether to export model-based regressors. For this model they are: "ev_c", "ev_nc", "pe_c", "pe_nc", "dv". + vb + Whether to use variational inference to approximately draw from a + posterior distribution. Defaults to ``False``. + inc_postpred + Include trial-level posterior predictive simulations in + model output (may greatly increase file size). Defaults to ``False``. + adapt_delta + Floating point value representing the target acceptance probability of a new + sample in the MCMC chain. Must be between 0 and 1. See note below. + stepsize + Integer value specifying the size of each leapfrog step that the MCMC sampler + can take on each new iteration. See note below. + max_treedepth + Integer value specifying how many leapfrog steps the MCMC sampler can take + on each new iteration. See note below. + **additional_args + Not used for this model. + + Returns + ------- + model_data + An ``hbayesdm.TaskModel`` instance with the following components: + + - ``model``: String value that is the name of the model ('prl_fictitious_woa_multipleB'). + - ``all_ind_pars``: Pandas DataFrame containing the summarized parameter values + (as specified by ``ind_pars``) for each subject. + - ``par_vals``: OrderedDict holding the posterior samples over different parameters. + - ``fit``: A PyStan StanFit object that contains the fitted Stan model. + - ``raw_data``: Pandas DataFrame containing the raw data used to fit the model, + as specified by the user. + - ``model_regressor``: Dict holding the extracted model-based regressors. + + Examples + -------- + + .. code:: python + + from hbayesdm import rhat, print_fit + from hbayesdm.models import prl_fictitious_woa_multipleB + + # Run the model and store results in "output" + output = prl_fictitious_woa_multipleB(data='example', niter=2000, nwarmup=1000, nchain=4, ncore=4) + + # Visually check convergence of the sampling chains (should look like "hairy caterpillars") + output.plot(type='trace') + + # Plot posterior distributions of the hyper-parameters (distributions should be unimodal) + output.plot() + + # Check Rhat values (all Rhat values should be less than or equal to 1.1) + rhat(output, less=1.1) + + # Show the LOOIC and WAIC model fit estimates + print_fit(output) + """ + return PrlFictitiousWoaMultipleb( + data=data, + niter=niter, + nwarmup=nwarmup, + nchain=nchain, + ncore=ncore, + nthin=nthin, + inits=inits, + ind_pars=ind_pars, + model_regressor=model_regressor, + vb=vb, + inc_postpred=inc_postpred, + adapt_delta=adapt_delta, + stepsize=stepsize, + max_treedepth=max_treedepth, + **additional_args) diff --git a/Python/tests/test_prl_fictitious_cnc.py b/Python/tests/test_prl_fictitious_cnc.py new file mode 100644 index 00000000..6118555a --- /dev/null +++ b/Python/tests/test_prl_fictitious_cnc.py @@ -0,0 +1,12 @@ +import pytest + +from hbayesdm.models import prl_fictitious_cnc + + +def test_prl_fictitious_cnc(): + _ = prl_fictitious_cnc( + data="example", niter=10, nwarmup=5, nchain=1, ncore=1) + + +if __name__ == '__main__': + pytest.main() diff --git a/Python/tests/test_prl_fictitious_cnc_multipleB.py b/Python/tests/test_prl_fictitious_cnc_multipleB.py new file mode 100644 index 00000000..1a3f0a9c --- /dev/null +++ b/Python/tests/test_prl_fictitious_cnc_multipleB.py @@ -0,0 +1,12 @@ +import pytest + +from hbayesdm.models import prl_fictitious_cnc_multipleB + + +def test_prl_fictitious_cnc_multipleB(): + _ = prl_fictitious_cnc_multipleB( + data="example", niter=10, nwarmup=5, nchain=1, ncore=1) + + +if __name__ == '__main__': + pytest.main() diff --git a/Python/tests/test_prl_fictitious_woa_multipleB.py b/Python/tests/test_prl_fictitious_woa_multipleB.py new file mode 100644 index 00000000..86bec63e --- /dev/null +++ b/Python/tests/test_prl_fictitious_woa_multipleB.py @@ -0,0 +1,12 @@ +import pytest + +from hbayesdm.models import prl_fictitious_woa_multipleB + + +def test_prl_fictitious_woa_multipleB(): + _ = prl_fictitious_woa_multipleB( + data="example", niter=10, nwarmup=5, nchain=1, ncore=1) + + +if __name__ == '__main__': + pytest.main() diff --git a/R/DESCRIPTION b/R/DESCRIPTION index c2d00ccd..b7d12a12 100644 --- a/R/DESCRIPTION +++ b/R/DESCRIPTION @@ -109,6 +109,9 @@ Collate: 'prl_fictitious_rp.R' 'prl_fictitious_rp_woa.R' 'prl_fictitious_woa.R' + 'prl_fictitious_woa_multipleB.R' + 'prl_fictitious_cnc.R' + 'prl_fictitious_cnc_multipleB.R' 'prl_rp.R' 'prl_rp_multipleB.R' 'pstRT_ddm.R' diff --git a/R/NAMESPACE b/R/NAMESPACE index 0f29aefe..36d6db4b 100644 --- a/R/NAMESPACE +++ b/R/NAMESPACE @@ -55,6 +55,9 @@ export(prl_fictitious_multipleB) export(prl_fictitious_rp) export(prl_fictitious_rp_woa) export(prl_fictitious_woa) +export(prl_fictitious_woa_multipleB) +export(prl_fictitious_cnc) +export(prl_fictitious_cnc_multipleB) export(prl_rp) export(prl_rp_multipleB) export(pstRT_ddm) diff --git a/R/R/prl_fictitious_cnc.R b/R/R/prl_fictitious_cnc.R new file mode 100644 index 00000000..f744df07 --- /dev/null +++ b/R/R/prl_fictitious_cnc.R @@ -0,0 +1,51 @@ +#' @templateVar MODEL_FUNCTION prl_fictitious_cnc +#' @templateVar CONTRIBUTOR \href{https://github.com/bugoverdose}{Jinwoo Jeong} <\email{jwjeong96@@gmail.com}> +#' @templateVar TASK_NAME Probabilistic Reversal Learning Task +#' @templateVar TASK_CODE prl +#' @templateVar TASK_CITE +#' @templateVar MODEL_NAME Fictitious Update Model, with chosen/not-chosen learning rates +#' @templateVar MODEL_CODE fictitious_cnc +#' @templateVar MODEL_CITE (Glascher et al., 2009) +#' @templateVar MODEL_TYPE Hierarchical +#' @templateVar DATA_COLUMNS "subjID", "choice", "outcome" +#' @templateVar PARAMETERS \code{eta_c} (learning rate for chosen option), \code{eta_nc} (learning rate for non-chosen option), \code{alpha} (indecision point), \code{beta} (inverse temperature) +#' @templateVar REGRESSORS "ev_c", "ev_nc", "pe_c", "pe_nc", "dv" +#' @templateVar POSTPREDS "y_pred" +#' @templateVar LENGTH_DATA_COLUMNS 3 +#' @templateVar DETAILS_DATA_1 \item{subjID}{A unique identifier for each subject in the data-set.} +#' @templateVar DETAILS_DATA_2 \item{choice}{Integer value representing the option chosen on that trial: 1 or 2.} +#' @templateVar DETAILS_DATA_3 \item{outcome}{Integer value representing the outcome of that trial (where reward == 1, and loss == -1).} +#' @templateVar LENGTH_ADDITIONAL_ARGS 0 +#' +#' @template model-documentation +#' +#' @export +#' @include hBayesDM_model.R +#' @include preprocess_funcs.R + +#' @references +#' Glascher, J., Hampton, A. N., & O'Doherty, J. P. (2009). Determining a Role for Ventromedial Prefrontal Cortex in Encoding Action-Based Value Signals During Reward-Related Decision Making. Cerebral Cortex, 19(2), 483-495. https://doi.org/10.1093/cercor/bhn098 +#' + + +prl_fictitious_cnc <- hBayesDM_model( + task_name = "prl", + model_name = "fictitious_cnc", + model_type = "", + data_columns = c("subjID", "choice", "outcome"), + parameters = list( + "eta_c" = c(0, 0.5, 1), + "eta_nc" = c(0, 0.5, 1), + "alpha" = c(-Inf, 0, Inf), + "beta" = c(0, 1, 10) + ), + additional_args = NULL, + regressors = list( + "ev_c" = 2, + "ev_nc" = 2, + "pe_c" = 2, + "pe_nc" = 2, + "dv" = 2 + ), + postpreds = c("y_pred"), + preprocess_func = prl_preprocess_func) diff --git a/R/R/prl_fictitious_cnc_multipleB.R b/R/R/prl_fictitious_cnc_multipleB.R new file mode 100644 index 00000000..557a3b03 --- /dev/null +++ b/R/R/prl_fictitious_cnc_multipleB.R @@ -0,0 +1,52 @@ +#' @templateVar MODEL_FUNCTION prl_fictitious_cnc_multipleB +#' @templateVar CONTRIBUTOR \href{https://github.com/bugoverdose}{Jinwoo Jeong} <\email{jwjeong96@@gmail.com}> +#' @templateVar TASK_NAME Probabilistic Reversal Learning Task +#' @templateVar TASK_CODE prl +#' @templateVar TASK_CITE +#' @templateVar MODEL_NAME Fictitious Update Model, with chosen/not-chosen learning rates +#' @templateVar MODEL_CODE fictitious_cnc +#' @templateVar MODEL_CITE (Glascher et al., 2009) +#' @templateVar MODEL_TYPE Multiple-Block Hierarchical +#' @templateVar DATA_COLUMNS "subjID", "block", "choice", "outcome" +#' @templateVar PARAMETERS \code{eta_c} (learning rate for chosen option), \code{eta_nc} (learning rate for non-chosen option), \code{alpha} (indecision point), \code{beta} (inverse temperature) +#' @templateVar REGRESSORS "ev_c", "ev_nc", "pe_c", "pe_nc", "dv" +#' @templateVar POSTPREDS "y_pred" +#' @templateVar LENGTH_DATA_COLUMNS 4 +#' @templateVar DETAILS_DATA_1 \item{subjID}{A unique identifier for each subject in the data-set.} +#' @templateVar DETAILS_DATA_2 \item{block}{A unique identifier for each of the multiple blocks within each subject.} +#' @templateVar DETAILS_DATA_3 \item{choice}{Integer value representing the option chosen on that trial: 1 or 2.} +#' @templateVar DETAILS_DATA_4 \item{outcome}{Integer value representing the outcome of that trial (where reward == 1, and loss == -1).} +#' @templateVar LENGTH_ADDITIONAL_ARGS 0 +#' +#' @template model-documentation +#' +#' @export +#' @include hBayesDM_model.R +#' @include preprocess_funcs.R + +#' @references +#' Glascher, J., Hampton, A. N., & O'Doherty, J. P. (2009). Determining a Role for Ventromedial Prefrontal Cortex in Encoding Action-Based Value Signals During Reward-Related Decision Making. Cerebral Cortex, 19(2), 483-495. https://doi.org/10.1093/cercor/bhn098 +#' + + +prl_fictitious_cnc_multipleB <- hBayesDM_model( + task_name = "prl", + model_name = "fictitious_cnc", + model_type = "multipleB", + data_columns = c("subjID", "block", "choice", "outcome"), + parameters = list( + "eta_c" = c(0, 0.5, 1), + "eta_nc" = c(0, 0.5, 1), + "alpha" = c(-Inf, 0, Inf), + "beta" = c(0, 1, 10) + ), + additional_args = NULL, + regressors = list( + "ev_c" = 3, + "ev_nc" = 3, + "pe_c" = 3, + "pe_nc" = 3, + "dv" = 3 + ), + postpreds = c("y_pred"), + preprocess_func = prl_multipleB_preprocess_func) diff --git a/R/R/prl_fictitious_woa_multipleB.R b/R/R/prl_fictitious_woa_multipleB.R new file mode 100644 index 00000000..82a773a8 --- /dev/null +++ b/R/R/prl_fictitious_woa_multipleB.R @@ -0,0 +1,50 @@ +#' @templateVar MODEL_FUNCTION prl_fictitious_woa_multipleB +#' @templateVar CONTRIBUTOR \href{https://github.com/bugoverdose}{Jinwoo Jeong} <\email{jwjeong96@@gmail.com}> +#' @templateVar TASK_NAME Probabilistic Reversal Learning Task +#' @templateVar TASK_CODE prl +#' @templateVar TASK_CITE +#' @templateVar MODEL_NAME Fictitious Update Model, without alpha (indecision point) +#' @templateVar MODEL_CODE fictitious_woa +#' @templateVar MODEL_CITE (Glascher et al., 2009) +#' @templateVar MODEL_TYPE Multiple-Block Hierarchical +#' @templateVar DATA_COLUMNS "subjID", "block", "choice", "outcome" +#' @templateVar PARAMETERS \code{eta} (learning rate), \code{beta} (inverse temperature) +#' @templateVar REGRESSORS "ev_c", "ev_nc", "pe_c", "pe_nc", "dv" +#' @templateVar POSTPREDS "y_pred" +#' @templateVar LENGTH_DATA_COLUMNS 4 +#' @templateVar DETAILS_DATA_1 \item{subjID}{A unique identifier for each subject in the data-set.} +#' @templateVar DETAILS_DATA_2 \item{block}{A unique identifier for each of the multiple blocks within each subject.} +#' @templateVar DETAILS_DATA_3 \item{choice}{Integer value representing the option chosen on that trial: 1 or 2.} +#' @templateVar DETAILS_DATA_4 \item{outcome}{Integer value representing the outcome of that trial (where reward == 1, and loss == -1).} +#' @templateVar LENGTH_ADDITIONAL_ARGS 0 +#' +#' @template model-documentation +#' +#' @export +#' @include hBayesDM_model.R +#' @include preprocess_funcs.R + +#' @references +#' Glascher, J., Hampton, A. N., & O'Doherty, J. P. (2009). Determining a Role for Ventromedial Prefrontal Cortex in Encoding Action-Based Value Signals During Reward-Related Decision Making. Cerebral Cortex, 19(2), 483-495. https://doi.org/10.1093/cercor/bhn098 +#' + + +prl_fictitious_woa_multipleB <- hBayesDM_model( + task_name = "prl", + model_name = "fictitious_woa", + model_type = "multipleB", + data_columns = c("subjID", "block", "choice", "outcome"), + parameters = list( + "eta" = c(0, 0.5, 1), + "beta" = c(0, 1, 10) + ), + additional_args = NULL, + regressors = list( + "ev_c" = 3, + "ev_nc" = 3, + "pe_c" = 3, + "pe_nc" = 3, + "dv" = 3 + ), + postpreds = c("y_pred"), + preprocess_func = prl_multipleB_preprocess_func) diff --git a/R/inst/plotting/plot_functions.R b/R/inst/plotting/plot_functions.R index dbcf2b3e..f800c6e2 100644 --- a/R/inst/plotting/plot_functions.R +++ b/R/inst/plotting/plot_functions.R @@ -178,7 +178,7 @@ plot_prl_fictitious_rp_woa <- function(obj, fontSize = 10, ncols = 3, binSize = return(h_all) } -plot_prl_fictitious_woa <- function(obj, fontSize = 10, ncols = 2, binSize = 30) { +plot_prl_fictitious_woa <- plot_prl_fictitious_woa_multipleB <- function(obj, fontSize = 10, ncols = 2, binSize = 30) { pars = obj$parVals h1 = plotDist(sample = pars$mu_eta, fontSize = fontSize, binSize = binSize, xLim = c(0,1), xLab = expression(paste(eta, " (Learning Rate)"))) h2 = plotDist(sample = pars$mu_beta, fontSize = fontSize, binSize = binSize, xLab = expression(paste(beta, " (Inverse Temp.)"))) @@ -186,6 +186,15 @@ plot_prl_fictitious_woa <- function(obj, fontSize = 10, ncols = 2, binSize = 30) return(h_all) } +plot_prl_fictitious_cnc <- plot_prl_fictitious_cnc_multipleB <- function(obj, fontSize = 10, ncols = 3, binSize = 30) { + pars = obj$parVals + h1 = plotDist(sample = pars$mu_eta_c, fontSize = fontSize, binSize = binSize, xLim = c(0,1), xLab = expression(paste(eta_c, " (Chosen. Learning Rate)"))) + h2 = plotDist(sample = pars$mu_eta_nc, fontSize = fontSize, binSize = binSize, xLim = c(0,1), xLab = expression(paste(eta_nc, " (Not-Chosen. Learning Rate)"))) + h3 = plotDist(sample = pars$mu_beta, fontSize = fontSize, binSize = binSize, xLab = expression(paste(beta, " (Inverse Temp.)"))) + h_all = multiplot(h1, h2, h3, cols = ncols) + return(h_all) +} + plot_prl_rp <- plot_prl_rp_multipleB <- function(obj, fontSize = 10, ncols = 3, binSize = 30) { pars = obj$parVals h1 = plotDist(sample = pars$mu_Apun, fontSize = fontSize, binSize = binSize, xLim = c(0,1), xLab = expression(paste(A[pun], " (Pun. Learning Rate)"))) diff --git a/R/tests/testthat/test_prl_fictitious_cnc.R b/R/tests/testthat/test_prl_fictitious_cnc.R new file mode 100644 index 00000000..05359466 --- /dev/null +++ b/R/tests/testthat/test_prl_fictitious_cnc.R @@ -0,0 +1,10 @@ +context("Test prl_fictitious_cnc") +library(hBayesDM) + +test_that("Test prl_fictitious_cnc", { + # Do not run this test on CRAN + skip_on_cran() + + expect_output(prl_fictitious_cnc( + data = "example", niter = 10, nwarmup = 5, nchain = 1, ncore = 1)) +}) diff --git a/R/tests/testthat/test_prl_fictitious_cnc_multipleB.R b/R/tests/testthat/test_prl_fictitious_cnc_multipleB.R new file mode 100644 index 00000000..c0b887a1 --- /dev/null +++ b/R/tests/testthat/test_prl_fictitious_cnc_multipleB.R @@ -0,0 +1,10 @@ +context("Test prl_fictitious_cnc_multipleB") +library(hBayesDM) + +test_that("Test prl_fictitious_cnc_multipleB", { + # Do not run this test on CRAN + skip_on_cran() + + expect_output(prl_fictitious_cnc_multipleB( + data = "example", niter = 10, nwarmup = 5, nchain = 1, ncore = 1)) +}) diff --git a/R/tests/testthat/test_prl_fictitious_woa_multipleB.R b/R/tests/testthat/test_prl_fictitious_woa_multipleB.R new file mode 100644 index 00000000..df8f6db9 --- /dev/null +++ b/R/tests/testthat/test_prl_fictitious_woa_multipleB.R @@ -0,0 +1,10 @@ +context("Test prl_fictitious_woa_multipleB") +library(hBayesDM) + +test_that("Test prl_fictitious_woa_multipleB", { + # Do not run this test on CRAN + skip_on_cran() + + expect_output(prl_fictitious_woa_multipleB( + data = "example", niter = 10, nwarmup = 5, nchain = 1, ncore = 1)) +}) diff --git a/commons/models/prl_fictitious_cnc.yml b/commons/models/prl_fictitious_cnc.yml new file mode 100644 index 00000000..0c3d9c12 --- /dev/null +++ b/commons/models/prl_fictitious_cnc.yml @@ -0,0 +1,47 @@ +task_name: + code: prl + desc: Probabilistic Reversal Learning Task + cite: +model_name: + code: fictitious_cnc + desc: Fictitious Update Model, with chosen/not-chosen learning rates + cite: + - Glascher, J., Hampton, A. N., & O'Doherty, J. P. (2009). Determining a Role for + Ventromedial Prefrontal Cortex in Encoding Action-Based Value Signals During Reward-Related + Decision Making. Cerebral Cortex, 19(2), 483-495. https://doi.org/10.1093/cercor/bhn098 +model_type: + code: + desc: Hierarchical +notes: +contributors: + - name: Jinwoo Jeong + email: jwjeong96@gmail.com + link: https://github.com/bugoverdose +data_columns: + subjID: A unique identifier for each subject in the data-set. + choice: "Integer value representing the option chosen on that trial: 1 or 2." + outcome: + Integer value representing the outcome of that trial (where reward == 1, + and loss == -1). +parameters: + eta_c: + desc: learning rate for chosen option + info: [0, 0.5, 1] + eta_nc: + desc: learning rate for non-chosen option + info: [0, 0.5, 1] + alpha: + desc: indecision point + info: [-Inf, 0, Inf] + beta: + desc: inverse temperature + info: [0, 1, 10] +regressors: + ev_c: 2 + ev_nc: 2 + pe_c: 2 + pe_nc: 2 + dv: 2 +postpreds: + - y_pred +additional_args: diff --git a/commons/models/prl_fictitious_cnc_multipleB.yml b/commons/models/prl_fictitious_cnc_multipleB.yml new file mode 100644 index 00000000..d57c2b68 --- /dev/null +++ b/commons/models/prl_fictitious_cnc_multipleB.yml @@ -0,0 +1,48 @@ +task_name: + code: prl + desc: Probabilistic Reversal Learning Task + cite: +model_name: + code: fictitious_cnc + desc: Fictitious Update Model, with chosen/not-chosen learning rates + cite: + - Glascher, J., Hampton, A. N., & O'Doherty, J. P. (2009). Determining a Role for + Ventromedial Prefrontal Cortex in Encoding Action-Based Value Signals During Reward-Related + Decision Making. Cerebral Cortex, 19(2), 483-495. https://doi.org/10.1093/cercor/bhn098 +model_type: + code: multipleB + desc: Multiple-Block Hierarchical +notes: +contributors: + - name: Jinwoo Jeong + email: jwjeong96@gmail.com + link: https://github.com/bugoverdose +data_columns: + subjID: A unique identifier for each subject in the data-set. + block: A unique identifier for each of the multiple blocks within each subject. + choice: "Integer value representing the option chosen on that trial: 1 or 2." + outcome: + Integer value representing the outcome of that trial (where reward == 1, + and loss == -1). +parameters: + eta_c: + desc: learning rate for chosen option + info: [0, 0.5, 1] + eta_nc: + desc: learning rate for non-chosen option + info: [0, 0.5, 1] + alpha: + desc: indecision point + info: [-Inf, 0, Inf] + beta: + desc: inverse temperature + info: [0, 1, 10] +regressors: + ev_c: 3 + ev_nc: 3 + pe_c: 3 + pe_nc: 3 + dv: 3 +postpreds: + - y_pred +additional_args: diff --git a/commons/models/prl_fictitious_woa_multipleB.yml b/commons/models/prl_fictitious_woa_multipleB.yml new file mode 100644 index 00000000..9a76a88c --- /dev/null +++ b/commons/models/prl_fictitious_woa_multipleB.yml @@ -0,0 +1,42 @@ +task_name: + code: prl + desc: Probabilistic Reversal Learning Task + cite: +model_name: + code: fictitious_woa + desc: Fictitious Update Model, without alpha (indecision point) + cite: + - Glascher, J., Hampton, A. N., & O'Doherty, J. P. (2009). Determining a Role for + Ventromedial Prefrontal Cortex in Encoding Action-Based Value Signals During Reward-Related + Decision Making. Cerebral Cortex, 19(2), 483-495. https://doi.org/10.1093/cercor/bhn098 +model_type: + code: multipleB + desc: Multiple-Block Hierarchical +notes: +contributors: + - name: Jinwoo Jeong + email: jwjeong96@gmail.com + link: https://github.com/bugoverdose +data_columns: + subjID: A unique identifier for each subject in the data-set. + block: A unique identifier for each of the multiple blocks within each subject. + choice: "Integer value representing the option chosen on that trial: 1 or 2." + outcome: + Integer value representing the outcome of that trial (where reward == 1, + and loss == -1). +parameters: + eta: + desc: learning rate + info: [0, 0.5, 1] + beta: + desc: inverse temperature + info: [0, 1, 10] +regressors: + ev_c: 3 + ev_nc: 3 + pe_c: 3 + pe_nc: 3 + dv: 3 +postpreds: + - y_pred +additional_args: diff --git a/commons/stan_files/prl_fictitious_cnc.stan b/commons/stan_files/prl_fictitious_cnc.stan new file mode 100644 index 00000000..3489146b --- /dev/null +++ b/commons/stan_files/prl_fictitious_cnc.stan @@ -0,0 +1,179 @@ +#include /pre/license.stan + +/** + * Probabilistic Reversal Learning (PRL) Task + * + * Fictitious Update Model (Glascher et al., 2008, Cerebral Cortex) + */ + +data { + int N; // Number of subjects + int T; // Maximum number of trials across subjects + int Tsubj[N]; // Number of trials/blocks for each subject + int choice[N, T]; // The choices subjects made + real outcome[N, 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[4] mu_pr; + vector[4] sigma; + + // Subject-level raw parameters (for Matt trick) + vector[N] eta_c_pr; // learning rate + vector[N] eta_nc_pr; // learning rate + vector[N] alpha_pr; // indecision point + vector[N] beta_pr; // inverse temperature +} + +transformed parameters { + // Transform subject-level raw parameters + vector[N] eta_c; + vector[N] eta_nc; + vector[N] alpha; + vector[N] beta; + + for (i in 1:N) { + eta_c[i] = Phi_approx(mu_pr[1] + sigma[1] * eta_c_pr[i]); + eta_nc[i] = Phi_approx(mu_pr[2] + sigma[2] * eta_nc_pr[i]); + beta[i] = Phi_approx(mu_pr[4] + sigma[4] * beta_pr[i]) * 10; + } + alpha = mu_pr[3] + sigma[3] * alpha_pr; +} + +model { + // Hyperparameters + mu_pr ~ normal(0, 1); + sigma[1] ~ normal(0, 0.2); + sigma[2] ~ normal(0, 0.2); + sigma[3] ~ cauchy(0, 1.0); + sigma[4] ~ normal(0, 0.2); + + // Individual parameters + eta_c_pr ~ normal(0, 1); + eta_nc_pr ~ normal(0, 1); + alpha_pr ~ normal(0, 1); + beta_pr ~ normal(0, 1); + + for (i in 1:N) { + // Define values + vector[2] ev; // expected value + vector[2] prob; // probability + real prob_1_; + + real PE; // prediction error + real PEnc; // fictitious prediction error (PE-non-chosen) + + // Initialize values + ev = initV; // initial ev values + + for (t in 1:(Tsubj[i])) { + // Compute action probabilities + prob[1] = 1 / (1 + exp(beta[i] * (alpha[i] - (ev[1] - ev[2])))); + prob_1_ = prob[1]; + prob[2] = 1 - prob_1_; + choice[i, t] ~ categorical(prob); + + // Prediction error + PE = outcome[i, t] - ev[choice[i, t]]; + PEnc = -outcome[i, t] - ev[3-choice[i, t]]; + + // Value updating (learning) + ev[choice[i, t]] += eta_c[i] * PE; + ev[3-choice[i, t]] += eta_nc[i] * PEnc; + } + } +} + +generated quantities { + // For group level parameters + real mu_eta_c; + real mu_eta_nc; + real mu_alpha; + real mu_beta; + + // For log likelihood calculation + real log_lik[N]; + + // For model regressors + real ev_c[N, T]; // Expected value of the chosen option + real ev_nc[N, T]; // Expected value of the non-chosen option + + real pe_c[N, T]; //Prediction error of the chosen option + real pe_nc[N, T]; //Prediction error of the non-chosen option + real dv[N, T]; //Decision value = PE_chosen - PE_non-chosen + + // For posterior predictive check + real y_pred[N, T]; + + // Set all posterior predictions, model regressors to 0 (avoids NULL values) + for (i in 1:N) { + for (t in 1:T) { + ev_c[i, t] = 0; + ev_nc[i, t] = 0; + + pe_c[i, t] = 0; + pe_nc[i, t] = 0; + dv[i, t] = 0; + + y_pred[i, t] = -1; + } + } + + mu_eta_c = Phi_approx(mu_pr[1]); + mu_eta_nc = Phi_approx(mu_pr[2]); + mu_alpha = mu_pr[3]; + mu_beta = Phi_approx(mu_pr[4]) * 10; + + { // local section, this saves time and space + for (i in 1:N) { + // Define values + vector[2] ev; // expected value + vector[2] prob; // probability + real prob_1_; + + real PE; // prediction error + real PEnc; // fictitious prediction error (PE-non-chosen) + + // Initialize values + ev = initV; // initial ev values + + log_lik[i] = 0; + + for (t in 1:(Tsubj[i])) { + // compute action probabilities + prob[1] = 1 / (1 + exp(beta[i] * (alpha[i] - (ev[1] - ev[2])))); + prob_1_ = prob[1]; + prob[2] = 1 - prob_1_; + + log_lik[i] += categorical_lpmf(choice[i, t] | prob); + + // generate posterior prediction for current trial + y_pred[i, t] = categorical_rng(prob); + + // prediction error + PE = outcome[i, t] - ev[choice[i, t]]; + PEnc = -outcome[i, t] - ev[3-choice[i, t]]; + + // Store values for model regressors + ev_c[i, t] = ev[choice[i, t]]; + ev_nc[i, t] = ev[3 - choice[i, t]]; + + pe_c[i, t] = PE; + pe_nc[i, t] = PEnc; + dv[i, t] = PE - PEnc; + + // value updating (learning) + ev[choice[i, t]] += eta_c[i] * PE; + ev[3-choice[i, t]] += eta_nc[i] * PEnc; + } + } + } +} diff --git a/commons/stan_files/prl_fictitious_cnc_multipleB.stan b/commons/stan_files/prl_fictitious_cnc_multipleB.stan new file mode 100644 index 00000000..bf0519f3 --- /dev/null +++ b/commons/stan_files/prl_fictitious_cnc_multipleB.stan @@ -0,0 +1,193 @@ +#include /pre/license.stan + +/** + * Probabilistic Reversal Learning (PRL) Task + * + * Fictitious Update Model (Glascher et al., 2008, Cerebral Cortex) + */ + +data { + int N; // Number of subjects + + int B; // Max number of blocks across subjects + int Bsubj[N]; // Number of blocks for each subject + + int T; // Max number of trials across subjects + int Tsubj[N, B]; // Number of trials/block for each subject + + int choice[N, B, T]; // Choice for each subject-block-trial + real outcome[N, B, T]; // Outcome (reward/loss) for each subject-block-trial +} + +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[4] mu_pr; + vector[4] sigma; + + // Subject-level raw parameters (for Matt trick) + vector[N] eta_c_pr; // learning rate + vector[N] eta_nc_pr; // learning rate + vector[N] alpha_pr; // indecision point + vector[N] beta_pr; // inverse temperature +} + +transformed parameters { + // Transform subject-level raw parameters + vector[N] eta_c; + vector[N] eta_nc; + vector[N] alpha; + vector[N] beta; + + for (i in 1:N) { + eta_c[i] = Phi_approx(mu_pr[1] + sigma[1] * eta_c_pr[i]); + eta_nc[i] = Phi_approx(mu_pr[2] + sigma[2] * eta_nc_pr[i]); + beta[i] = Phi_approx(mu_pr[4] + sigma[4] * beta_pr[i]) * 10; + } + alpha = mu_pr[3] + sigma[3] * alpha_pr; +} + +model { + // Hyperparameters + mu_pr ~ normal(0, 1); + sigma[1] ~ normal(0, 0.2); + sigma[2] ~ normal(0, 0.2); + sigma[3] ~ cauchy(0, 1.0); + sigma[4] ~ normal(0, 0.2); + + // Individual parameters + eta_c_pr ~ normal(0, 1); + eta_nc_pr ~ normal(0, 1); + 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 + vector[2] prob; // probability + real prob_1_; + + real PE; // prediction error + real PEnc; // fictitious prediction error (PE-non-chosen) + + // Initialize values + ev = initV; // initial ev values + + for (t in 1:(Tsubj[i, bIdx])) { // new + // compute action probabilities + prob[1] = 1 / (1 + exp(beta[i] * (alpha[i] - (ev[1] - ev[2])))); + prob_1_ = prob[1]; + prob[2] = 1 - prob_1_; + choice[i, bIdx, t] ~ categorical(prob); + //choice[i, t] ~ bernoulli(prob); + + // prediction error + PE = outcome[i, bIdx, t] - ev[choice[i, bIdx, t]]; //new + PEnc = -outcome[i, bIdx, t] - ev[3-choice[i, bIdx, t]]; //new + + // value updating (learning) + ev[choice[i, bIdx, t]] += eta_c[i] * PE; //new + ev[3-choice[i, bIdx, t]] += eta_nc[i] * PEnc; //new + } // end of t loop + } // end of bIdx loop + } // end of i loop +} + +generated quantities { + // For group level parameters + real mu_eta_c; + real mu_eta_nc; + real mu_alpha; + real 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_c[N, B, T]; //Prediction error of the chosen option + real pe_nc[N, B, T]; //Prediction error of the non-chosen option + real dv[N, B, T]; //Decision value = PE_chosen - PE_non-chosen + + // For posterior predictive check + real y_pred[N, B, T]; + + // Set all posterior predictions, model regressors to 0 (avoids 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_c[i, b, t] = 0; + pe_nc[i, b, t] = 0; + dv[i, b, t] = 0; + + y_pred[i, b, t] = -1; + } + } + } + + mu_eta_c = Phi_approx(mu_pr[1]); + mu_eta_nc = Phi_approx(mu_pr[2]); + mu_alpha = mu_pr[3]; + mu_beta = Phi_approx(mu_pr[4]) * 10; + + { // local section, this saves time and space + for (i in 1:N) { + + log_lik[i] = 0; + + for (bIdx in 1:Bsubj[i]) { + // Define values + vector[2] ev; // expected value + vector[2] prob; // probability + real prob_1_; + + real PE; // prediction error + real PEnc; // fictitious prediction error (PE-non-chosen) + + // Initialize values + ev = initV; // initial ev values + + for (t in 1:(Tsubj[i, bIdx])) { + // compute action probabilities + prob[1] = 1 / (1 + exp(beta[i] * (alpha[i] - (ev[1] - ev[2])))); + prob_1_ = prob[1]; + prob[2] = 1 - prob_1_; + + log_lik[i] += categorical_lpmf(choice[i, bIdx, t] | prob); //new + + // generate posterior prediction for current trial + y_pred[i, bIdx, t] = categorical_rng(prob); + + // prediction error + PE = outcome[i, bIdx, t] - ev[choice[i, bIdx, t]]; //new + PEnc = -outcome[i, bIdx, t] - ev[3-choice[i, bIdx, t]]; //new + + // Store values for model regressors + ev_c[i, bIdx, t] = ev[choice[i, bIdx, t]]; + ev_nc[i, bIdx, t] = ev[3 - choice[i, bIdx, t]]; + + pe_c[i, bIdx, t] = PE; + pe_nc[i, bIdx, t] = PEnc; + dv[i, bIdx, t] = PE - PEnc; + + // value updating (learning) + ev[choice[i, bIdx, t]] += eta_c[i] * PE; //new + ev[3-choice[i, bIdx, t]] += eta_nc[i] * PEnc; //new + } // end of t loop + } // end of bIdx loop + } + } +} + diff --git a/commons/stan_files/prl_fictitious_woa_multipleB.stan b/commons/stan_files/prl_fictitious_woa_multipleB.stan new file mode 100644 index 00000000..df62626c --- /dev/null +++ b/commons/stan_files/prl_fictitious_woa_multipleB.stan @@ -0,0 +1,178 @@ +#include /pre/license.stan + +/** + * Probabilistic Reversal Learning (PRL) Task + * + * Fictitious Update Model (Glascher et al., 2008, Cerebral Cortex) without alpha (indecision point) + */ + +data { + int N; // Number of subjects + + int B; // Max number of blocks across subjects + int Bsubj[N]; // Number of blocks for each subject + + int T; // Max number of trials across subjects + int Tsubj[N, B]; // Number of trials/block for each subject + + int choice[N, B, T]; // Choice for each subject-block-trial + real outcome[N, B, T]; // Outcome (reward/loss) for each subject-block-trial +} + +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[2] sigma; + + // Subject-level raw parameters (for Matt trick) + vector[N] eta_pr; // learning rate + vector[N] beta_pr; // inverse temperature +} + +transformed parameters { + // Transform subject-level raw parameters + vector[N] eta; + vector[N] beta; + + for (i in 1:N) { + eta[i] = Phi_approx(mu_pr[1] + sigma[1] * eta_pr[i]); + beta[i] = Phi_approx(mu_pr[2] + sigma[2] * beta_pr[i]) * 10; + } +} +model { + // Hyperparameters + mu_pr ~ normal(0, 1); + sigma[1] ~ normal(0, 0.2); + sigma[2] ~ normal(0, 0.2); + + // individual parameters + eta_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 + vector[2] prob; // probability + real prob_1_; + + real PE; // prediction error + real PEnc; // fictitious prediction error (PE-non-chosen) + + // Initialize values + ev = initV; // initial ev values + + for (t in 1:(Tsubj[i, bIdx])) { // new + // compute action probabilities + prob[1] = 1 / (1 + exp(beta[i] * (ev[2] - ev[1]))); + prob_1_ = prob[1]; + prob[2] = 1 - prob_1_; + choice[i, bIdx, t] ~ categorical(prob); + //choice[i, t] ~ bernoulli(prob); + + // prediction error + PE = outcome[i, bIdx, t] - ev[choice[i, bIdx, t]]; //new + PEnc = -outcome[i, bIdx, t] - ev[3-choice[i, bIdx, t]]; //new + + // value updating (learning) + ev[choice[i, bIdx, t]] += eta[i] * PE; //new + ev[3-choice[i, bIdx, t]] += eta[i] * PEnc; //new + } // end of t loop + } // end of bIdx loop + } // end of i loop +} + +generated quantities { + // For group level parameters + real mu_eta; + real 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_c[N, B, T]; //Prediction error of the chosen option + real pe_nc[N, B, T]; //Prediction error of the non-chosen option + real dv[N, B, T]; //Decision value = PE_chosen - PE_non-chosen + + // For posterior predictive check + real y_pred[N, B, T]; + + // Set all posterior predictions, model regressors to 0 (avoids 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_c[i, b, t] = 0; + pe_nc[i, b, t] = 0; + dv[i, b, t] = 0; + + y_pred[i, b, t] = -1; + } + } + } + + mu_eta = 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]) { + // Define values + vector[2] ev; // expected value + vector[2] prob; // probability + real prob_1_; + + real PE; // prediction error + real PEnc; // fictitious prediction error (PE-non-chosen) + + // Initialize values + ev = initV; // initial ev values + + for (t in 1:(Tsubj[i, bIdx])) { + // compute action probabilities + prob[1] = 1 / (1 + exp(beta[i] * (ev[2] - ev[1]))); + prob_1_ = prob[1]; + prob[2] = 1 - prob_1_; + + log_lik[i] += categorical_lpmf(choice[i, bIdx, t] | prob); //new + + // generate posterior prediction for current trial + y_pred[i, bIdx, t] = categorical_rng(prob); + + // prediction error + PE = outcome[i, bIdx, t] - ev[choice[i, bIdx, t]]; //new + PEnc = -outcome[i, bIdx, t] - ev[3-choice[i, bIdx, t]]; //new + + // Store values for model regressors + ev_c[i, bIdx, t] = ev[choice[i, bIdx, t]]; + ev_nc[i, bIdx, t] = ev[3 - choice[i, bIdx, t]]; + + pe_c[i, bIdx, t] = PE; + pe_nc[i, bIdx, t] = PEnc; + dv[i, bIdx, t] = PE - PEnc; + + // value updating (learning) + ev[choice[i, bIdx, t]] += eta[i] * PE; //new + ev[3-choice[i, bIdx, t]] += eta[i] * PEnc; //new + } // end of t loop + } // end of bIdx loop + } + } +} +