This is a tutorial for using the HDDMrl module to simultaneously estimate reinforcement learning parameters and decision parameters within a fully hierarchical bayesian estimation framework, including steps for sampling, assessing convergence, model fit, parameter recovery, and posterior predictive checks (model validation). The module uses the reinforcement learning drift diffusion model (RLDDM), a reinforcement learning model that replaces the standard “softmax” choice function with a drift diffusion process. The softmax and drift diffusion process is equivalent for capturing choice proportions, but the DDM also takes RT distributions into account; options are provided to also only fit RL parameters without RT. The RLDDM estimates trial-by-trial drift rate as a scaled difference in expected rewards (expected reward for upper bound alternative minus expected reward for lower bound alternative). Expected rewards are updated with a delta learning rule using either a single learning rate or with separate learning rates for positive and negative prediction errors. The model also includes the standard DDM-parameters. The RLDDM is described in detail in Pedersen, Frank & Biele (2017). (Note this approach differs from Frank et al (2015) who used HDDM to fit instrumental learning but did not allow for simultaneous estimation of learning parameters).
1. Background
2. Installing the module
3. How the RLDDM works
4. Structuring data
5. Running basic model
6. Checking results
7. Posterior predictive checks
8. Parameter recovery
9. Separate learning rates for positive and negative prediction errors
10. depends_on vs. split_by
11. Probabilistic binary outcomes vs. normally distributed outcomes
12. HDDMrlRegressor
13. Regular RL without RT
Traditional RL models typically assume static decision processes (e.g. softmax), and the DDM typically assumes static decision variables (stimuli are modeled with the same drift rate across trials). The RLDDM combines dynamic decision variables from RL and dynamic choice process from DDM by assuming trial-by-trial drift rate that depends on the difference in expected rewards, which are updated on each trial by a rate of the prediction error dependent on the learning rate. The potential benefit of the RLDDM is thus to gain a better insight into decision processes in instrumental learning by also accounting for speed of decision making.
The new version of HDDM (version 0.7.1) that includes the RLDDM is currently not uploaded to conda. So to install you would either have to use pip:
'pip install hddm'
OR
docker: 'pull madslupe/hddm', which runs hddm in jupyter notebook
The main idea of the RLDDM is that reward-based choices can be captured by an accumulation-to-bound process where drift rate is proportional to the difference in expected reward between options, and where expected rewards subsequently are updated in a trial-by-trial basis via reinforcement learning.
drift rate on each trial depends on difference in expected rewards for the two alternatives (q_up and q_low):
drift rate = (q_up - q_low) * scaling
where the scaling parameter describes the weight to put on the difference in q-values.
expeceted reward (q) for chosen option is updated according to delta learning rule :
q_chosen = q_chosen + alpha * (feedback-q_chosen)
where alpha weights the rate of learning on each trial.
So in principle all you need is the Wiener first passage time likelihood-function. The reason why HDDM is useful (and hence also HDDMrl) is that it automates the process of setting up and running your model, which tends to be very time consuming. So after structuring the data it is simple to run a model with HDDMrl. In particular it separates subjects and conditions (using the split_by-column, see next section) so that the updating process works correctly, which can be especially difficult to do for RL models.
The HDDMrl module was created to make it easier to model instrumental learning data with the RLDDM. If you are familiar with using HDDM it shouldn't be a big step to start using HDDMrl. Please refresh your memory by starting with the tutorial for HDDM first (especially critical if you have not used HDDM at all). Running HDDMrl does require a few extra steps compared to HDDM, and because the model includes increased potential for parameter colinearity (typically learning rate and the scaling parameter on drift rate) it is even more important to assess model fit, which will be covered below. Here are the most important steps for structuring your dataframe:
To illustrate how to run the model we will use example data from the learning phase of the probabilistic selection task (PST). During the learning phase of the PST subjects choose between two stimuli presented as Hiragana-letters (here represented as letters from the latin alphabet). There are three conditions with different probabilities of receiving reward (feedback=1) and non-reward (feedback=0). In the AB condition A is rewarded with 80% probability, B with 20%. In the CD condition C is rewarded with 70% probability and D with 30%, while in the EF condition E is rewarded with a 60% probability and F with 40%. The dataset is included in the data-folder in your installation of HDDM.
# import
import pandas as pd
import numpy as np
import hddm
from scipy import stats
import seaborn as sns
import matplotlib.pyplot as plt
import pymc
import kabuki
sns.set(style="white")
%matplotlib inline
from tqdm import tqdm
import warnings
warnings.filterwarnings("ignore", category=np.VisibleDeprecationWarning)
/home/krishn/anaconda3/envs/hddm/lib/python3.6/site-packages/tensorflow/python/framework/dtypes.py:526: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'. _np_qint8 = np.dtype([("qint8", np.int8, 1)]) /home/krishn/anaconda3/envs/hddm/lib/python3.6/site-packages/tensorflow/python/framework/dtypes.py:527: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'. _np_quint8 = np.dtype([("quint8", np.uint8, 1)]) /home/krishn/anaconda3/envs/hddm/lib/python3.6/site-packages/tensorflow/python/framework/dtypes.py:528: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'. _np_qint16 = np.dtype([("qint16", np.int16, 1)]) /home/krishn/anaconda3/envs/hddm/lib/python3.6/site-packages/tensorflow/python/framework/dtypes.py:529: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'. _np_quint16 = np.dtype([("quint16", np.uint16, 1)]) /home/krishn/anaconda3/envs/hddm/lib/python3.6/site-packages/tensorflow/python/framework/dtypes.py:530: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'. _np_qint32 = np.dtype([("qint32", np.int32, 1)]) /home/krishn/anaconda3/envs/hddm/lib/python3.6/site-packages/tensorflow/python/framework/dtypes.py:535: FutureWarning: Passing (type, 1) or '1type' as a synonym of type is deprecated; in a future version of numpy, it will be understood as (type, (1,)) / '(1,)type'. np_resource = np.dtype([("resource", np.ubyte, 1)])
tanh tanh tanh linear
/home/krishn/anaconda3/envs/hddm/lib/python3.6/site-packages/IPython/parallel.py:13: ShimWarning: The `IPython.parallel` package has been deprecated since IPython 4.0. You should import from ipyparallel instead. "You should import from ipyparallel instead.", ShimWarning)
# load data. you will find this dataset in your hddm-folder under hddm/examples/rlddm_data.csv
data = hddm.load_csv("rlddm_data.csv")
# check structure
data.head()
subj_idx | response | cond | rt | trial | split_by | feedback | q_init | |
---|---|---|---|---|---|---|---|---|
0 | 42 | 0.0 | CD | 1.255 | 1.0 | 1 | 0.0 | 0.5 |
1 | 42 | 1.0 | EF | 0.778 | 1.0 | 2 | 0.0 | 0.5 |
2 | 42 | 1.0 | AB | 0.647 | 1.0 | 0 | 1.0 | 0.5 |
3 | 42 | 1.0 | AB | 0.750 | 2.0 | 0 | 1.0 | 0.5 |
4 | 42 | 0.0 | EF | 0.772 | 2.0 | 2 | 1.0 | 0.5 |
The columns in the datafile represent: subj_idx (subject id), response (1=best option), cond (identifies condition, but not used in model), rt (in seconds), 0=worst option), trial (the trial-iteration for a subject within each condition), split_by (identifying condition, used for running the model), feedback (whether the response given was rewarded or not), q_init (the initial q-value used for the model, explained above).
# run the model by calling hddm.HDDMrl (instead of hddm.HDDM for normal HDDM)
m = hddm.HDDMrl(data)
# set sample and burn-in
m.sample(1500, burn=500, dbname="traces.db", db="pickle")
# print stats to get an overview of posterior distribution of estimated parameters
m.print_stats()
{} () No model attribute --> setting up standard HDDM Includes supplied: () printing self.nn False Set model to ddm [-----------------100%-----------------] 1500 of 1500 complete in 143.7 sec mean std 2.5q 25q 50q 75q 97.5q mc err a 1.7433 0.0813138 1.58494 1.68967 1.74494 1.7954 1.9187 0.00292143 a_std 0.406522 0.0610479 0.312902 0.3623 0.401829 0.441929 0.558704 0.00266828 a_subj.3 2.00633 0.10129 1.83542 1.93472 1.99867 2.06634 2.23372 0.00443095 a_subj.4 1.89855 0.0539594 1.79832 1.86133 1.89825 1.93209 2.01204 0.002394 a_subj.5 1.48664 0.0741329 1.34377 1.4389 1.48745 1.53229 1.6436 0.00351514 a_subj.6 2.39376 0.100291 2.19236 2.32632 2.39406 2.4608 2.5925 0.00447612 a_subj.8 2.49353 0.132293 2.23753 2.40396 2.49135 2.57203 2.77695 0.00695635 a_subj.12 1.90207 0.0684117 1.77964 1.85209 1.89715 1.94514 2.05585 0.00313201 a_subj.17 1.40287 0.0822611 1.22774 1.35337 1.40235 1.45487 1.56221 0.00487811 a_subj.18 1.90812 0.125013 1.67796 1.81644 1.90952 1.99562 2.15789 0.00679012 a_subj.19 2.08064 0.085651 1.92067 2.02268 2.08202 2.13764 2.24183 0.00446387 a_subj.20 1.71228 0.0794146 1.55761 1.65782 1.70888 1.76279 1.8681 0.00369196 a_subj.22 1.19307 0.0335414 1.12915 1.17129 1.19207 1.21472 1.26232 0.00158518 a_subj.23 1.93739 0.0895423 1.76888 1.87749 1.93578 1.99651 2.12028 0.00376337 a_subj.24 1.68377 0.0952974 1.514 1.62176 1.68111 1.74522 1.87941 0.00482898 a_subj.26 2.23807 0.0803151 2.08342 2.18239 2.23925 2.29066 2.40204 0.00362931 a_subj.33 1.56624 0.066382 1.44514 1.51945 1.56395 1.61366 1.70245 0.00274963 a_subj.34 1.81803 0.0883775 1.64683 1.75351 1.82055 1.87853 1.99149 0.00376558 a_subj.35 1.83533 0.0975695 1.64512 1.77156 1.83433 1.9008 2.03218 0.00488217 a_subj.36 1.31686 0.0685052 1.19867 1.2683 1.31396 1.35797 1.46606 0.00292263 a_subj.39 1.57041 0.0530064 1.47084 1.53222 1.57151 1.60465 1.68414 0.00183899 a_subj.42 1.72821 0.0757661 1.58521 1.67929 1.7261 1.77948 1.8817 0.00388646 a_subj.50 1.46984 0.0556356 1.3657 1.42899 1.46951 1.51015 1.57257 0.00366184 a_subj.52 2.1227 0.0885946 1.95207 2.06009 2.12394 2.18527 2.29723 0.00371883 a_subj.56 1.47025 0.0482088 1.37467 1.44058 1.46858 1.50214 1.56976 0.00245215 a_subj.59 1.31005 0.0817522 1.16181 1.25224 1.30721 1.36636 1.46954 0.00524295 a_subj.63 1.8639 0.0883703 1.71047 1.801 1.8624 1.92234 2.04382 0.00375231 a_subj.71 1.24626 0.0340768 1.17937 1.22192 1.24647 1.27036 1.31251 0.00140437 a_subj.75 0.97873 0.0365919 0.913123 0.951568 0.977399 1.00419 1.05616 0.00179393 a_subj.80 2.187 0.0965787 2.00696 2.12386 2.18203 2.24945 2.37954 0.00388193 v 4.08209 0.737245 2.64887 3.57663 4.09291 4.5691 5.5414 0.0458143 v_std 3.05484 0.575751 2.07762 2.64448 2.98722 3.40931 4.34097 0.0386874 v_subj.3 5.3036 2.4037 1.33138 3.46797 5.08761 6.83865 10.67 0.180208 v_subj.4 0.936614 0.141934 0.679051 0.839064 0.931719 1.02271 1.23745 0.00588382 v_subj.5 5.06788 2.26694 1.35156 3.39827 4.84072 6.55001 10.3183 0.158873 v_subj.6 4.58494 2.4764 1.19679 2.43696 4.25878 6.21347 9.97977 0.218163 v_subj.8 3.1593 1.1972 1.73287 2.40894 2.89231 3.50444 6.66669 0.0985572 v_subj.12 1.83104 0.224408 1.44052 1.67555 1.81061 1.9754 2.30106 0.0107878 v_subj.17 6.93416 2.22622 3.44608 5.39757 6.56186 8.19763 12.2283 0.168861 v_subj.18 6.12272 1.59534 3.55255 4.97474 5.93287 7.09307 9.95941 0.122813 v_subj.19 11.7997 1.9933 8.60538 10.4711 11.5523 12.8174 16.4139 0.156147 v_subj.20 4.72675 1.23807 2.99018 3.87997 4.53657 5.3181 7.96896 0.104607 v_subj.22 2.19815 0.283286 1.68683 1.9937 2.18769 2.39485 2.73678 0.0128246 v_subj.23 2.4331 1.50181 0.959062 1.4395 1.90909 2.89597 6.60491 0.128198 v_subj.24 5.12042 0.753118 3.67884 4.60547 5.10993 5.62647 6.63844 0.0334464 v_subj.26 0.602543 0.863859 0.24547 0.404374 0.466825 0.555059 1.98844 0.0705156 v_subj.33 2.73658 2.10267 0.282828 1.14428 2.1263 3.89181 8.07948 0.120839 v_subj.34 4.11408 2.32502 1.11361 2.29809 3.65292 5.45819 9.77123 0.15346 v_subj.35 3.14171 1.87384 1.09744 1.80252 2.57169 3.96941 7.90195 0.155472 v_subj.36 3.87701 2.00542 1.52318 2.39526 3.27531 4.79694 9.1923 0.161708 v_subj.39 2.0298 1.62953 0.34513 0.935087 1.45793 2.58893 6.68562 0.107378 v_subj.42 5.16806 1.89934 2.28067 3.90109 4.91005 6.18106 9.81468 0.164759 v_subj.50 5.57649 1.79009 2.85308 4.24585 5.34125 6.59774 9.84759 0.141939 v_subj.52 3.18577 1.84878 1.16424 1.85057 2.71706 3.88391 8.41286 0.159 v_subj.56 2.05017 1.91277 0.258263 0.703241 1.36834 2.86095 7.32663 0.14203 v_subj.59 10.9111 1.82836 7.62548 9.62643 10.8841 12.1672 14.6141 0.133133 v_subj.63 2.22187 1.32504 0.988662 1.4741 1.91208 2.45281 6.4643 0.105276 v_subj.71 2.17006 1.78241 0.603637 1.08317 1.48305 2.50471 7.27104 0.149917 v_subj.75 3.71969 0.552382 2.70805 3.34183 3.69385 4.06496 4.9227 0.0277269 v_subj.80 4.03709 2.17552 0.986179 2.27813 3.66944 5.33457 9.111 0.170157 t 0.432885 0.0415403 0.356521 0.404121 0.431653 0.462548 0.511254 0.00178714 t_std 0.232643 0.0394899 0.169301 0.203101 0.227423 0.257429 0.315877 0.00179352 t_subj.3 1.06317 0.0271769 0.999802 1.04743 1.06651 1.08276 1.10684 0.00129274 t_subj.4 0.527059 0.0179398 0.488522 0.515254 0.528076 0.539991 0.557209 0.000826279 t_subj.5 0.528425 0.0158093 0.493437 0.519202 0.529042 0.538151 0.559482 0.000756408 t_subj.6 0.398728 0.0316425 0.330672 0.379382 0.399754 0.41924 0.461281 0.00145695 t_subj.8 0.663496 0.0406703 0.573784 0.638121 0.667493 0.692995 0.734004 0.00220312 t_subj.12 0.408965 0.0147454 0.378986 0.398904 0.410359 0.41963 0.436432 0.000711611 t_subj.17 0.500941 0.0182598 0.472563 0.491633 0.499139 0.505433 0.563134 0.00142979 t_subj.18 0.433627 0.0282406 0.371316 0.416532 0.438131 0.455078 0.47849 0.00153929 t_subj.19 0.419058 0.0150433 0.387912 0.409382 0.419564 0.429151 0.4471 0.000741052 t_subj.20 0.51628 0.0127393 0.491093 0.507661 0.516726 0.525528 0.538053 0.000610317 t_subj.22 0.33496 0.00518669 0.324114 0.331928 0.335304 0.338463 0.344235 0.000253281 t_subj.23 0.481226 0.0226842 0.431716 0.466665 0.483849 0.497623 0.519431 0.000941342 t_subj.24 0.453762 0.013039 0.425115 0.445645 0.454502 0.462701 0.47683 0.000579678 t_subj.26 0.440285 0.0338985 0.372521 0.414917 0.443772 0.468245 0.49323 0.00164635 t_subj.33 0.194727 0.0156474 0.164383 0.184158 0.194446 0.204649 0.225083 0.000671091 t_subj.34 0.344477 0.0311183 0.276129 0.327923 0.345049 0.361748 0.413056 0.00126463 t_subj.35 0.326443 0.0258584 0.273692 0.309526 0.32619 0.343838 0.377765 0.00119399 t_subj.36 0.44915 0.0134952 0.415288 0.443814 0.452192 0.458502 0.466911 0.000650595 t_subj.39 0.619389 0.0110796 0.595795 0.612534 0.620065 0.627238 0.638565 0.000401916 t_subj.42 0.393542 0.0155789 0.362458 0.383731 0.394086 0.404417 0.422266 0.000830944 t_subj.50 0.518358 0.0214237 0.479266 0.498422 0.524998 0.536913 0.548587 0.00150307 t_subj.52 0.514745 0.0238354 0.46713 0.499788 0.515337 0.530868 0.56368 0.00104662 t_subj.56 0.118365 0.0108243 0.0959725 0.11129 0.119001 0.12587 0.139238 0.00054126 t_subj.59 0.37923 0.0228073 0.332184 0.363319 0.377838 0.39906 0.415741 0.00167972 t_subj.63 0.478244 0.0193531 0.437445 0.466752 0.479564 0.491179 0.513794 0.000800474 t_subj.71 0.162673 0.0054732 0.150777 0.159311 0.162957 0.166633 0.172557 0.000235596 t_subj.75 0.260777 0.00626259 0.244606 0.258207 0.262153 0.264765 0.269822 0.000295097 t_subj.80 0.0417991 0.0176072 0.0114696 0.028693 0.041848 0.0536039 0.0781576 0.000700406 alpha -2.46715 0.356073 -3.157 -2.69308 -2.48242 -2.24392 -1.74121 0.0195181 alpha_std 1.46992 0.275563 1.0209 1.28032 1.45203 1.61981 2.09465 0.0155966 alpha_subj.3 -3.43548 0.740258 -4.48635 -3.90149 -3.56324 -3.08342 -1.45227 0.0608327 alpha_subj.4 0.191786 0.595174 -0.801725 -0.213263 0.126076 0.52319 1.65872 0.0268693 alpha_subj.5 -3.28297 0.771161 -4.52565 -3.7898 -3.39336 -2.93093 -1.27709 0.0562957 alpha_subj.6 -3.69185 0.841373 -4.92836 -4.33939 -3.87057 -3.1469 -1.67019 0.0743001 alpha_subj.8 -2.29198 0.611898 -3.76992 -2.61173 -2.22875 -1.86592 -1.29704 0.047517 alpha_subj.12 -0.780933 0.425905 -1.62084 -1.05694 -0.786673 -0.502726 0.0835042 0.0206346 alpha_subj.17 -2.81596 0.484016 -3.76465 -3.14894 -2.81048 -2.5293 -1.77409 0.0371539 alpha_subj.18 -2.75671 0.392851 -3.44184 -3.03159 -2.77331 -2.50711 -1.88007 0.030011 alpha_subj.19 -3.7218 0.226315 -4.21729 -3.85334 -3.70381 -3.58146 -3.29365 0.0171907 alpha_subj.20 -2.37897 0.438575 -3.25441 -2.65439 -2.3874 -2.09008 -1.51432 0.0362439 alpha_subj.22 -0.975383 0.212324 -1.36486 -1.12002 -0.979027 -0.845113 -0.547454 0.00859101 alpha_subj.23 -2.29157 0.971566 -4.14016 -2.95807 -2.24846 -1.50684 -0.582361 0.0809592 alpha_subj.24 -1.49115 0.231074 -1.95491 -1.64755 -1.47849 -1.33537 -1.05145 0.00874484 alpha_subj.26 0.477848 1.54335 -5.25248 0.166103 0.700563 1.2395 2.52254 0.12061 alpha_subj.33 -3.25956 1.33119 -5.54438 -4.22104 -3.29747 -2.31559 -0.650837 0.0832202 alpha_subj.34 -2.9963 0.924299 -4.55411 -3.65087 -3.1139 -2.40745 -1.05377 0.0637804 alpha_subj.35 -2.68266 0.985176 -4.42859 -3.41642 -2.76104 -1.94719 -0.782148 0.0807771 alpha_subj.36 -2.35113 0.964568 -3.91305 -3.07013 -2.41253 -1.6805 -0.415103 0.0775148 alpha_subj.39 -3.07825 1.23011 -5.28098 -4.03919 -3.06089 -2.16057 -0.719406 0.0801558 alpha_subj.42 -3.14098 0.59356 -4.12935 -3.52327 -3.22374 -2.80894 -1.82638 0.0516964 alpha_subj.50 -3.78284 0.656364 -4.82933 -4.26283 -3.87635 -3.39796 -2.22331 0.0510745 alpha_subj.52 -3.24081 0.75318 -4.72941 -3.78718 -3.26367 -2.68309 -1.89825 0.0642394 alpha_subj.56 -3.88438 1.4086 -6.07973 -4.95683 -4.04305 -2.87882 -0.949276 0.107409 alpha_subj.59 -2.54493 0.225305 -2.98248 -2.69149 -2.55072 -2.39734 -2.06605 0.0149846 alpha_subj.63 -2.04163 0.949775 -4.16 -2.56008 -1.95309 -1.44029 -0.246927 0.0683846 alpha_subj.71 -2.84608 1.54153 -5.48198 -4.07678 -2.63727 -1.8612 0.43124 0.123145 alpha_subj.75 -1.37871 0.44339 -2.24472 -1.67589 -1.36566 -1.06067 -0.558518 0.0223029 alpha_subj.80 -3.22538 0.795376 -4.50105 -3.82306 -3.32231 -2.75346 -1.47653 0.0631666 DIC: 10127.947410 deviance: 10059.166818 pD: 68.780592
Interpreting output from print_stats:
The model estimates group mean and standard deviation parameters and subject parameters for the following latent variables:
a = decision threshold
v = scaling parameter
t = non-decision time
alpha = learning rate, note that it's not bound between 0 and 1. to transform take inverse logit: np.exp(alpha)/(1+np.exp(alpha))
The columns represent the mean, standard deviation and quantiles of the approximated posterior distribution of each parameter
There are a few things to note that is different from the normal HDDM model.
First of all, the estimated learning rate does not necessarily fall between 0 and 1. This is because it is estimated as a normal distribution for purposes of sampling hierarchically and then transformed by an inverse logit function to 0<alpha<1. So to interpret alpha as learning rate you have to transform the samples in the trace back with np.exp(alpha)/(1+np.exp(alpha)). And if you estimate separate learning rates for positive and negative prediction errors (see here) then you get learning rate for negative prediction errors with np.exp(alpha)/(1+np.exp(alpha)) and positive prediction errors with np.exp(pos_alpha)/(1+np.exp(pos_alpha)).
Second, the v-parameter in the output is the scaling factor that is multiplied by the difference in q-values, so it is not the actual drift rate (or rather, it is the equivalent drift rate when the difference in Q values is exactly 1).
# plot the posteriors of parameters
m.plot_posteriors()
Plotting a Plotting a_std Plotting v Plotting v_std Plotting t Plotting t_std Plotting alpha Plotting alpha_std
Fig. The mixing of the posterior distribution and autocorrelation looks ok.
The Gelman-Rubin statistic is a test of whether the chains in the model converges. The Gelman-Ruben statistic measures the degree of variation between and within chains. Values close to 1 indicate convergence and that there is small variation between chains, i.e. that they end up as the same distribution across chains. A common heuristic is to assume convergence if all values are below 1.1. To run this you need to run multiple models, combine them and perform the Gelman-Rubin statistic:
# estimate convergence
from kabuki.analyze import gelman_rubin
models = []
for i in range(3):
m = hddm.HDDMrl(data=data)
m.sample(1500, burn=500, dbname="traces.db", db="pickle")
models.append(m)
gelman_rubin(models)
np.max(list(gelman_rubin(models).values()))
The model seems to have converged, i.e. the Gelman-Rubin statistic is below 1.1 for all parameters. It is important to always run this test, especially for more complex models (as with separate learning rates for positive and negative prediction errors). So now we can combine these three models to get a better approximation of the posterior distribution.
# Combine the models we ran to test for convergence.
m = kabuki.utils.concat_models(models)
Another test of the model is to look at collinearity. If the estimation of parameters is very codependent (correlation is strong) it can indicate that their variance trades off, in particular if there is a negative correlation. The following plot shows there is generally low correlation across all combinations of parameters. It does not seem to be the case for this dataset, but common for RLDDM is a negative correlation between learning rate and the scaling factor, similar to what's usually observed between learning rate and inverse temperature for RL models that uses softmax as the choice rule (e.g. Daw, 2011).
alpha, t, a, v = m.nodes_db.node[["alpha", "t", "a", "v"]]
samples = {"alpha": alpha.trace(), "t": t.trace(), "a": a.trace(), "v": v.trace()}
samp = pd.DataFrame(data=samples)
def corrfunc(x, y, **kws):
r, _ = stats.pearsonr(x, y)
ax = plt.gca()
ax.annotate("r = {:.2f}".format(r), xy=(0.1, 0.9), xycoords=ax.transAxes)
g = sns.PairGrid(samp, palette=["red"])
g.map_upper(plt.scatter, s=10)
g.map_diag(sns.distplot, kde=False)
g.map_lower(sns.kdeplot, cmap="Blues_d")
g.map_lower(corrfunc)
g.savefig("matrix_plot.png")
An important test of the model is its ability to recreate the observed data. This can be tested with posterior predictive checks, which involves simulating data using estimated parameters and comparing observed and simulated results.
The first step then is to extract the traces from the estimated model. The function get_traces() gives you the samples (row) from the approaximated posterior distribution for all of the estimated group and subject parameters (column).
traces = m.get_traces()
traces.head()
Now that we have the traces the next step is to simulate data using the estimated parameters. The RLDDM includes a function to simulate data. Here's an example of how to use the simulation-function for RLDDM. This example explains how to generate data with binary outcomes. See here for an example on simulating data with normally distributed outcomes. Inputs to function:
a = decision threshold
t = non-decision time
alpha = learning rate
pos_alpha = defaults to 0. if given it defines the learning rate for positive prediction errors. alpha then becomes the learning rate_ for negative prediction errors.
scaler = the scaling factor that is multiplied with the difference in q-values to calculate trial-by-trial drift rate
p_upper = the probability of reward for the option represented by the upper boundary. The current version thus only works for outcomes that are either 1 or 0
p_lower = the probability of reward for the option represented by the lower boundary.
subjs = number of subjects to simulate data for.
split_by = define the condition which makes it easier to append data from different conditions.
size = number of trials per subject.
hddm.generate.gen_rand_rlddm_data(
a=1,
t=0.3,
alpha=0.2,
scaler=2,
p_upper=0.8,
p_lower=0.2,
subjs=1,
split_by=0,
size=10,
)
How to interpret columns in the resulting dataframe
q_up = expected reward for option represented by upper boundary
q_low = expected reward for option represented by lower boundary
sim_drift = the drift rate for each trial calculated as: (q_up-q_low)*scaler
response = simulated choice
rt = simulated response time
feedback = observed feedback for chosen option
subj_idx = subject id (starts at 0)
split_by = condition as integer
trial = current trial (starts at 1)
Now that we know how to extract traces and simulate data we can combine this to create a dataset similar to our observed data. This process is currently not automated but the following is an example code using the dataset we analyzed above.
from tqdm import tqdm # progress tracker
# create empty dataframe to store simulated data
sim_data = pd.DataFrame()
# create a column samp to be used to identify the simulated data sets
data["samp"] = 0
# load traces
traces = m.get_traces()
# decide how many times to repeat simulation process. repeating this multiple times is generally recommended,
# as it better captures the uncertainty in the posterior distribution, but will also take some time
for i in tqdm(range(1, 51)):
# randomly select a row in the traces to use for extracting parameter values
sample = np.random.randint(0, traces.shape[0] - 1)
# loop through all subjects in observed data
for s in data.subj_idx.unique():
# get number of trials for each condition.
size0 = len(
data[(data["subj_idx"] == s) & (data["split_by"] == 0)].trial.unique()
)
size1 = len(
data[(data["subj_idx"] == s) & (data["split_by"] == 1)].trial.unique()
)
size2 = len(
data[(data["subj_idx"] == s) & (data["split_by"] == 2)].trial.unique()
)
# set parameter values for simulation
a = traces.loc[sample, "a_subj." + str(s)]
t = traces.loc[sample, "t_subj." + str(s)]
scaler = traces.loc[sample, "v_subj." + str(s)]
alphaInv = traces.loc[sample, "alpha_subj." + str(s)]
# take inverse logit of estimated alpha
alpha = np.exp(alphaInv) / (1 + np.exp(alphaInv))
# simulate data for each condition changing only values of size, p_upper, p_lower and split_by between conditions.
sim_data0 = hddm.generate.gen_rand_rlddm_data(
a=a,
t=t,
scaler=scaler,
alpha=alpha,
size=size0,
p_upper=0.8,
p_lower=0.2,
split_by=0,
)
sim_data1 = hddm.generate.gen_rand_rlddm_data(
a=a,
t=t,
scaler=scaler,
alpha=alpha,
size=size1,
p_upper=0.7,
p_lower=0.3,
split_by=1,
)
sim_data2 = hddm.generate.gen_rand_rlddm_data(
a=a,
t=t,
scaler=scaler,
alpha=alpha,
size=size2,
p_upper=0.6,
p_lower=0.4,
split_by=2,
)
# append the conditions
sim_data0 = sim_data0.append([sim_data1, sim_data2], ignore_index=True)
# assign subj_idx
sim_data0["subj_idx"] = s
# identify that these are simulated data
sim_data0["type"] = "simulated"
# identify the simulated data
sim_data0["samp"] = i
# append data from each subject
sim_data = sim_data.append(sim_data0, ignore_index=True)
# combine observed and simulated data
ppc_data = data[
["subj_idx", "response", "split_by", "rt", "trial", "feedback", "samp"]
].copy()
ppc_data["type"] = "observed"
ppc_sdata = sim_data[
["subj_idx", "response", "split_by", "rt", "trial", "feedback", "type", "samp"]
].copy()
ppc_data = ppc_data.append(ppc_sdata)
ppc_data.to_csv("ppc_data_tutorial.csv")
Now that we have a dataframe with both observed and simulated data we can plot to see whether the simulated data are able to capture observed choice and reaction times. To capture the uncertainty in the simulated data we want to identify how much choice and reaction differs across the simulated data sets. A good measure of this is to calculate the highest posterior density/highest density interval for summary scores of the generated data. Below we calculate highest posterior density with an alpha set to 0.1, which means that we are describing the range of the 90% most likely values.
# for practical reasons we only look at the first 40 trials for each subject in a given condition
plot_ppc_data = ppc_data[ppc_data.trial < 41].copy()
# bin trials to for smoother estimate of response proportion across learning
plot_ppc_data["bin_trial"] = pd.cut(
plot_ppc_data.trial, 11, labels=np.linspace(0, 10, 11)
).astype("int64")
# calculate means for each sample
sums = (
plot_ppc_data.groupby(["bin_trial", "split_by", "samp", "type"])
.mean()
.reset_index()
)
# calculate the overall mean response across samples
ppc_sim = sums.groupby(["bin_trial", "split_by", "type"]).mean().reset_index()
# initiate columns that will have the upper and lower bound of the hpd
ppc_sim["upper_hpd"] = 0
ppc_sim["lower_hpd"] = 0
for i in range(0, ppc_sim.shape[0]):
# calculate the hpd/hdi of the predicted mean responses across bin_trials
hdi = pymc.utils.hpd(
sums.response[
(sums["bin_trial"] == ppc_sim.bin_trial[i])
& (sums["split_by"] == ppc_sim.split_by[i])
& (sums["type"] == ppc_sim.type[i])
],
alpha=0.1,
)
ppc_sim.loc[i, "upper_hpd"] = hdi[1]
ppc_sim.loc[i, "lower_hpd"] = hdi[0]
# calculate error term as the distance from upper bound to mean
ppc_sim["up_err"] = ppc_sim["upper_hpd"] - ppc_sim["response"]
ppc_sim["low_err"] = ppc_sim["response"] - ppc_sim["lower_hpd"]
ppc_sim["model"] = "RLDDM_single_learning"
ppc_sim.to_csv("ppc_choicedata_tutorial.csv")
# plotting evolution of choice proportion for best option across learning for observed and simulated data.
fig, axs = plt.subplots(figsize=(15, 5), nrows=1, ncols=3, sharex=True, sharey=True)
for i in range(0, 3):
ax = axs[i]
d = ppc_sim[(ppc_sim.split_by == i) & (ppc_sim.type == "simulated")]
ax.errorbar(
d.bin_trial,
d.response,
yerr=[d.low_err, d.up_err],
label="simulated",
color="orange",
)
d = ppc_sim[(ppc_sim.split_by == i) & (ppc_sim.type == "observed")]
ax.plot(d.bin_trial, d.response, linewidth=3, label="observed")
ax.set_title("split_by = %i" % i, fontsize=20)
ax.set_ylabel("mean response")
ax.set_xlabel("trial")
plt.legend()
fig.savefig("PPCchoice.pdf")
Fig. The plots display the rate of choosing the best option (response = 1) across learning and condition. The model generates data (orange) that closely follows the observed behavior (blue), with the exception of overpredicting performance early in the most difficult condition (split_by=2). Uncertainty in the generated data is captured by the 90% highest density interval of the means across simulated datasets.
# set reaction time to be negative for lower bound responses (response=0)
plot_ppc_data["reaction time"] = np.where(
plot_ppc_data["response"] == 1, plot_ppc_data.rt, 0 - plot_ppc_data.rt
)
# plotting evolution of choice proportion for best option across learning for observed and simulated data. We use bins of trials because plotting individual trials would be very noisy.
g = sns.FacetGrid(plot_ppc_data, col="split_by", hue="type")
g.map(sns.kdeplot, "reaction time", bw=0.05).set_ylabels("Density")
g.add_legend()
g.savefig("PPCrt_dist.pdf")
Fig. Density plots of observed and predicted reaction time across conditions. RTs for lower boundary choices (i.e. worst option choices) are set to be negative (0-RT) to be able to separate upper and lower bound responses.
To validate the RLDDM we ran a parameter recovery study to test to which degree the model can recover the parameter values used to simulate data. To do this we generated 81 synthetic datasets with 50 subjects performing 70 trials each. The 81 datasets were simulated using all combinations of three plausible parameter values for decision threshold, non-decision time, learning rate and the scaling parameter onto drift rate.
We can plot simulated together with the estimated values to test the models ability to recover parameters, and to see if there are any values that are more difficult to recover than others.
param_recovery = hddm.load_csv("recovery_sim_est_rlddm.csv")
g = sns.catplot(x="a", y="e_a", data=param_recovery, palette="Set1")
g.set_axis_labels("Simulated threshold", "Estimated threshold")
plt.title("Decision threshold")
g.savefig("Threshold_recovery.pdf")
g = sns.catplot(x="alpha", y="e_alphaT", data=param_recovery, palette="Set1")
g.set_axis_labels("Simulated alpha", "Estimated alpha")
plt.title("Learning rate")
g.savefig("Alpha_recovery.pdf")
g = sns.catplot(x="scaler", y="e_v", data=param_recovery, palette="Set1")
g.set_axis_labels("Simulated scaling", "Estimated scaling")
plt.title("Scaling drift rate")
g.savefig("Scaler_recovery.pdf")
g = sns.catplot(x="t", y="e_t", data=param_recovery, palette="Set1")
g.set_axis_labels("Simulated NDT", "Estimated NDT")
plt.title("Non-decision time")
g.savefig("NDT_recovery.pdf")
Fig. The correlation between simulated and estimated parameter values are high, which means recovery is good. There is somewhat worse recovery for the learning rate and scaling parameter, which makes sense given that they to a degree can explain the same variance (see below).
Several studies have reported differences in updating of expected rewards following positive and negative prediction errors (e.g. to capture differences between D1 and D2 receptor function). To model asymmetric updating rates for positive and negative prediction errors you can set dual=True in the model. This will produce two estimated learning rates; alpha and pos_alpha, of which alpha then becomes the estimated learning rate for negative prediction errors.
# set dual=True to model separate learning rates for positive and negative prediction errors.
m_dual = hddm.HDDMrl(data, dual=True)
# set sample and burn-in
m_dual.sample(1500, burn=500, dbname="traces.db", db="pickle")
# print stats to get an overview of posterior distribution of estimated parameters
m_dual.print_stats()
m_dual.plot_posteriors()
Fig. There's more autocorrelation in this model compared to the one with a single learning rate. First, let's test whether it converges.
# estimate convergence
models = []
for i in range(3):
m = hddm.HDDMrl(data=data, dual=True)
m.sample(1500, burn=500, dbname="traces.db", db="pickle")
models.append(m)
# get max gelman-statistic value. shouldn't be higher than 1.1
np.max(list(gelman_rubin(models).values()))
gelman_rubin(models)
Convergence looks good, i.e. no parameters with gelman-rubin statistic > 1.1.
# Create a new model that has all traces concatenated
# of individual models.
m_dual = kabuki.utils.concat_models(models)
And then we can have a look at the joint posterior distribution:
alpha, t, a, v, pos_alpha = m_dual.nodes_db.node[["alpha", "t", "a", "v", "pos_alpha"]]
samples = {
"alpha": alpha.trace(),
"pos_alpha": pos_alpha.trace(),
"t": t.trace(),
"a": a.trace(),
"v": v.trace(),
}
samp = pd.DataFrame(data=samples)
def corrfunc(x, y, **kws):
r, _ = stats.pearsonr(x, y)
ax = plt.gca()
ax.annotate("r = {:.2f}".format(r), xy=(0.1, 0.9), xycoords=ax.transAxes)
g = sns.PairGrid(samp, palette=["red"])
g.map_upper(plt.scatter, s=10)
g.map_diag(sns.distplot, kde=False)
g.map_lower(sns.kdeplot, cmap="Blues_d")
g.map_lower(corrfunc)
g.savefig("matrix_plot.png")
Fig. The correlation between parameters is generally low.
The DIC for this dual learning rate model is better than for the single learning rate model. We can therefore check whether we can detect this improvement in the ability to recreate choice and RT patterns:
# create empty dataframe to store simulated data
sim_data = pd.DataFrame()
# create a column samp to be used to identify the simulated data sets
data["samp"] = 0
# get traces, note here we extract traces from m_dual
traces = m_dual.get_traces()
# decide how many times to repeat simulation process. repeating this multiple times is generally recommended as it better captures the uncertainty in the posterior distribution, but will also take some time
for i in tqdm(range(1, 51)):
# randomly select a row in the traces to use for extracting parameter values
sample = np.random.randint(0, traces.shape[0] - 1)
# loop through all subjects in observed data
for s in data.subj_idx.unique():
# get number of trials for each condition.
size0 = len(
data[(data["subj_idx"] == s) & (data["split_by"] == 0)].trial.unique()
)
size1 = len(
data[(data["subj_idx"] == s) & (data["split_by"] == 1)].trial.unique()
)
size2 = len(
data[(data["subj_idx"] == s) & (data["split_by"] == 2)].trial.unique()
)
# set parameter values for simulation
a = traces.loc[sample, "a_subj." + str(s)]
t = traces.loc[sample, "t_subj." + str(s)]
scaler = traces.loc[sample, "v_subj." + str(s)]
# when generating data with two learning rates pos_alpha represents learning rate for positive prediction errors and alpha for negative prediction errors
alphaInv = traces.loc[sample, "alpha_subj." + str(s)]
pos_alphaInv = traces.loc[sample, "pos_alpha_subj." + str(s)]
# take inverse logit of estimated alpha and pos_alpha
alpha = np.exp(alphaInv) / (1 + np.exp(alphaInv))
pos_alpha = np.exp(pos_alphaInv) / (1 + np.exp(pos_alphaInv))
# simulate data for each condition changing only values of size, p_upper, p_lower and split_by between conditions.
sim_data0 = hddm.generate.gen_rand_rlddm_data(
a=a,
t=t,
scaler=scaler,
alpha=alpha,
pos_alpha=pos_alpha,
size=size0,
p_upper=0.8,
p_lower=0.2,
split_by=0,
)
sim_data1 = hddm.generate.gen_rand_rlddm_data(
a=a,
t=t,
scaler=scaler,
alpha=alpha,
pos_alpha=pos_alpha,
size=size1,
p_upper=0.7,
p_lower=0.3,
split_by=1,
)
sim_data2 = hddm.generate.gen_rand_rlddm_data(
a=a,
t=t,
scaler=scaler,
alpha=alpha,
pos_alpha=pos_alpha,
size=size2,
p_upper=0.6,
p_lower=0.4,
split_by=2,
)
# append the conditions
sim_data0 = sim_data0.append([sim_data1, sim_data2], ignore_index=True)
# assign subj_idx
sim_data0["subj_idx"] = s
# identify that these are simulated data
sim_data0["type"] = "simulated"
# identify the simulated data
sim_data0["samp"] = i
# append data from each subject
sim_data = sim_data.append(sim_data0, ignore_index=True)
# combine observed and simulated data
ppc_dual_data = data[
["subj_idx", "response", "split_by", "rt", "trial", "feedback", "samp"]
].copy()
ppc_dual_data["type"] = "observed"
ppc_dual_sdata = sim_data[
["subj_idx", "response", "split_by", "rt", "trial", "feedback", "type", "samp"]
].copy()
ppc_dual_data = ppc_dual_data.append(ppc_dual_sdata)
# for practical reasons we only look at the first 40 trials for each subject in a given condition
plot_ppc_dual_data = ppc_dual_data[ppc_dual_data.trial < 41].copy()
# bin trials to for smoother estimate of response proportion across learning
plot_ppc_dual_data["bin_trial"] = pd.cut(
plot_ppc_dual_data.trial, 11, labels=np.linspace(0, 10, 11)
).astype("int64")
# calculate means for each sample
sums = (
plot_ppc_dual_data.groupby(["bin_trial", "split_by", "samp", "type"])
.mean()
.reset_index()
)
# calculate the overall mean response across samples
ppc_dual_sim = sums.groupby(["bin_trial", "split_by", "type"]).mean().reset_index()
# initiate columns that will have the upper and lower bound of the hpd
ppc_dual_sim["upper_hpd"] = 0
ppc_dual_sim["lower_hpd"] = 0
for i in range(0, ppc_dual_sim.shape[0]):
# calculate the hpd/hdi of the predicted mean responses across bin_trials
hdi = pymc.utils.hpd(
sums.response[
(sums["bin_trial"] == ppc_dual_sim.bin_trial[i])
& (sums["split_by"] == ppc_dual_sim.split_by[i])
& (sums["type"] == ppc_dual_sim.type[i])
],
alpha=0.1,
)
ppc_dual_sim.loc[i, "upper_hpd"] = hdi[1]
ppc_dual_sim.loc[i, "lower_hpd"] = hdi[0]
# calculate error term as the distance from upper bound to mean
ppc_dual_sim["up_err"] = ppc_dual_sim["upper_hpd"] - ppc_dual_sim["response"]
ppc_dual_sim["low_err"] = ppc_dual_sim["response"] - ppc_dual_sim["lower_hpd"]
ppc_dual_sim["model"] = "RLDDM_dual_learning"
# plotting evolution of choice proportion for best option across learning for observed and simulated data.
fig, axs = plt.subplots(figsize=(15, 5), nrows=1, ncols=3, sharex=True, sharey=True)
for i in range(0, 3):
ax = axs[i]
d = ppc_dual_sim[(ppc_dual_sim.split_by == i) & (ppc_dual_sim.type == "simulated")]
ax.errorbar(
d.bin_trial,
d.response,
yerr=[d.low_err, d.up_err],
label="simulated",
color="orange",
)
d = ppc_sim[(ppc_dual_sim.split_by == i) & (ppc_dual_sim.type == "observed")]
ax.plot(d.bin_trial, d.response, linewidth=3, label="observed")
ax.set_title("split_by = %i" % i, fontsize=20)
ax.set_ylabel("mean response")
ax.set_xlabel("trial")
plt.legend()
Fig. The plots display the rate of choosing the best option (response = 1) across learning and condition. The model generates data (orange) that closely follows the observed behavior (blue), with the exception of performance early in the most difficult condition (split_by=2).
To get a better sense of differences in ability to predict data between the single and dual learning rate model we can plot them together:
# plotting evolution of choice proportion for best option across learning for observed and simulated data. Compared for model with single and dual learning rate.
fig, axs = plt.subplots(figsize=(15, 5), nrows=1, ncols=3, sharex=True, sharey=True)
for i in range(0, 3):
ax = axs[i]
d_single = ppc_sim[(ppc_sim.split_by == i) & (ppc_sim.type == "simulated")]
# slightly move bin_trial to avoid overlap in errorbars
d_single["bin_trial"] += 0.2
ax.errorbar(
d_single.bin_trial,
d_single.response,
yerr=[d_single.low_err, d_single.up_err],
label="simulated_single",
color="orange",
)
d_dual = ppc_dual_sim[
(ppc_dual_sim.split_by == i) & (ppc_dual_sim.type == "simulated")
]
ax.errorbar(
d_dual.bin_trial,
d_dual.response,
yerr=[d_dual.low_err, d_dual.up_err],
label="simulated_dual",
color="green",
)
d = ppc_sim[(ppc_dual_sim.split_by == i) & (ppc_dual_sim.type == "observed")]
ax.plot(d.bin_trial, d.response, linewidth=3, label="observed")
ax.set_title("split_by = %i" % i, fontsize=20)
ax.set_ylabel("mean response")
ax.set_xlabel("trial")
plt.xlim(-0.5, 10.5)
plt.legend()
Fig. The predictions from the model with two learning rates are not very different from the model with single learning rate, and a similar overprediction of performance early on for the most difficult condition (split_by =2).
plot_ppc_data["type_compare"] = np.where(
plot_ppc_data["type"] == "observed",
plot_ppc_data["type"],
"simulated_single_learning",
)
plot_ppc_dual_data["type_compare"] = np.where(
plot_ppc_dual_data["type"] == "observed",
plot_ppc_dual_data["type"],
"simulated_dual_learning",
)
dual_vs_single_pcc = plot_ppc_data.append(plot_ppc_dual_data)
dual_vs_single_pcc["reaction time"] = np.where(
dual_vs_single_pcc["response"] == 1,
dual_vs_single_pcc.rt,
0 - dual_vs_single_pcc.rt,
)
# plotting evolution of choice proportion for best option across learning for observed and simulated data. We use bins of trials because plotting individual trials would be very noisy.
g = sns.FacetGrid(dual_vs_single_pcc, col="split_by", hue="type_compare", height=5)
g.map(sns.kdeplot, "reaction time", bw=0.01).set_ylabels("Density")
g.add_legend()
Fig. Again there's not a big difference between the two models. Both models slightly overpredict performance for the medium (split_by =1) and hard (split_by = 2) conditions, as identified by lower densities for the negative (worst option choices) in the simulated compared to observed data.
To interpret the parameter estimates for alpha and pos_alpha you have to transform them with the inverse logit where learning rate for negative prediction error is alpha and learning rate for positive prediction errors is pos_alpha. For this dataset the learning rate is estimated to be higher for positive than negative prediction errors.
# plot alpha for positive and negative learning rate
traces = m_dual.get_traces()
neg_alpha = np.exp(traces["alpha"]) / (1 + np.exp(traces["alpha"]))
pos_alpha = np.exp(traces["pos_alpha"]) / (1 + np.exp(traces["pos_alpha"]))
sns.kdeplot(
neg_alpha, color="r", label="neg_alpha: " + str(np.round(np.mean(neg_alpha), 3))
)
sns.kdeplot(
pos_alpha, color="b", label="pos_alpha: " + str(np.round(np.mean(pos_alpha), 3))
)
Fig. The positive learning rate is estimated to be stronger than the negative learning rate. Sticky choice, tendencies to repeat choices, could be driving some of this difference. The current model does not allow to test for this, however, but it could be tested in the future if we implement a regression version of RLDDM (similar to HDDMRegressor).
Here's how you would simulate data with a learning rate for positive and negative predictions of 0.2 and 0.4, respectively:
hddm.generate.gen_rand_rlddm_data(
a=1, t=0.3, alpha=0.2, pos_alpha=0.4, scaler=2, p_upper=0.8, p_lower=0.2, size=10
)
HDDMrl can be used to estimate separate parameters just as in the standard HDDM. But in RL you typically estimate the same learning rates and inverse temperature across conditions. That's one reason why you have to specify condition in the split_by-column instead of depends_on. (The other is that if you use depends_on the expected rewards will not get updated properly). But depends_on is still useful, for example if you want to estimate the effect of group on parameters. As an example we can simulate a dataset with two groups that have different decision thresholds:
data1 = hddm.generate.gen_rand_rlddm_data(
a=1, t=0.3, alpha=0.4, scaler=2, p_upper=0.8, p_lower=0.2, subjs=50, size=50
)
data1["group"] = "group1"
data2 = hddm.generate.gen_rand_rlddm_data(
a=2, t=0.3, alpha=0.4, scaler=2, p_upper=0.8, p_lower=0.2, subjs=50, size=50
)
data2["group"] = "group2"
group_data = data1.append(data2)
group_data["q_init"] = 0.5
m = hddm.HDDMrl(
group_data, depends_on={"v": "group", "a": "group", "t": "group", "alpha": "group"}
)
m.sample(1500, burn=500, dbname="traces.db", db="pickle")
# the plot shows that the model was able to recover the different decision threshold across groups.
a_group1, a_group2 = m.nodes_db.node[["a(group1)", "a(group2)"]]
hddm.analyze.plot_posterior_nodes([a_group1, a_group2])
plt.xlabel("decision threshold")
plt.ylabel("Posterior probability")
plt.xlim(0.7, 2.3)
plt.title("Posterior of decision threshold group means")
The examples so far have all been using a task structure where the outcomes are binary and probabilistic. But the model can also be applied to other types of outcomes. Here we show how you can generate and model data with normally distributed outcomes. As you will see you don't have to do any modifications to the model estimation process, but you have to change the input for generating data. Also note that the scaling parameter (v) will scale negatively with the values of the observed outcomes because the combined drift rate needs to be plausible.
# This is how we generated data so far, defining the probability of reward (1) for actions/stimuli associated with upper and lower boundary.
# binary probabilistic outcomes
hddm.generate.gen_rand_rlddm_data(
a=2, t=0.3, scaler=2, alpha=0.2, size=10, p_upper=0.2, p_lower=0.8
)
# If instead the outcomes are drawn from a normal distribution you will have to set binary_outcome to False and instead of p_upper and p_upper define the mean (mu) and sd
# of the normal distribution for both alternatives. Note that we change the initial q-value to 0, and that we reduce the scaling factor.
# normally distributed outcomes
hddm.generate.gen_rand_rlddm_data(
a=2,
t=0.3,
scaler=0.2,
alpha=0.2,
size=10,
mu_upper=8,
mu_lower=2,
sd_upper=1,
sd_lower=1,
binary_outcome=False,
q_init=0,
)
# We can generate a dataset where 30 subjects perform 50 trials each. Note that we set the scaler to be lower than for the binary outcomes as otherwise
# the resulting drift will be unrealistically high.
norm_data = hddm.generate.gen_rand_rlddm_data(
a=2,
t=0.3,
scaler=0.2,
alpha=0.2,
size=50,
subjs=30,
mu_upper=8,
mu_lower=2,
sd_upper=2,
sd_lower=2,
binary_outcome=False,
q_init=0,
)
# and then we can do estimation as usual
# but first we need to define inital q-value
norm_data["q_init"] = 0
m_norm = hddm.HDDMrl(norm_data)
m_norm.sample(1500, burn=500, dbname="traces.db", db="pickle")
m_norm.print_stats()
As of version 0.7.6. HDDM includes a module for estimating the impact of continuous regressor onto RLDDM parameters. The module, called HDDMrlRegressor, works the same way as the HDDMRegressor for the normal DDM. The method allows estimation of the association of e.g. neural measures onto parameters. To illustrate the method we extend the function to generate rlddm_data by adding a normally distributed regressor and including a coefficient called 'neural'.Note that to run the HDDMrlRegressor you need to include alpha when specifying the model. For more information on how to set up regressor models look at the tutorial for HDDM.
# function to generate rlddm-data that adds a neural regressor to decision threshold
def gen_rand_reg_rlddm_data(
a,
t,
scaler,
alpha,
neural,
size=1,
p_upper=1,
p_lower=0,
z=0.5,
q_init=0.5,
split_by=0,
subjs=1,
):
all_data = []
n = size
# set sd for variables to generate subject-parameters from group distribution
sd_t = 0.02
sd_a = 0.1
sd_alpha = 0.1
sd_v = 0.25
# save parameter values as group-values
tg = t
ag = a
alphag = alpha
scalerg = scaler
for s in range(0, subjs):
t = (
np.maximum(0.05, np.random.normal(loc=tg, scale=sd_t, size=1))
if subjs > 1
else tg
)
a = (
np.maximum(0.05, np.random.normal(loc=ag, scale=sd_a, size=1))
if subjs > 1
else ag
)
alpha = (
np.minimum(
np.minimum(
np.maximum(0.001, np.random.normal(loc=alphag, scale=sd_a, size=1)),
alphag + alphag,
),
1,
)
if subjs > 1
else alphag
)
scaler = (
np.random.normal(loc=scalerg, scale=sd_v, size=1) if subjs > 1 else scalerg
)
# create a normalized regressor that is combined with the neural coefficient to create trial-by-trial values for decision threshold
neural_reg = np.random.normal(0, 1, size=n)
q_up = np.tile([q_init], n)
q_low = np.tile([q_init], n)
response = np.tile([0.5], n)
feedback = np.tile([0.5], n)
rt = np.tile([0], n)
rew_up = np.random.binomial(1, p_upper, n).astype(float)
rew_low = np.random.binomial(1, p_lower, n).astype(float)
sim_drift = np.tile([0], n)
subj_idx = np.tile([s], n)
d = {
"q_up": q_up,
"q_low": q_low,
"sim_drift": sim_drift,
"rew_up": rew_up,
"rew_low": rew_low,
"response": response,
"rt": rt,
"feedback": feedback,
"subj_idx": subj_idx,
"split_by": split_by,
"trial": 1,
"neural_reg": neural_reg,
}
df = pd.DataFrame(data=d)
df = df[
[
"q_up",
"q_low",
"sim_drift",
"rew_up",
"rew_low",
"response",
"rt",
"feedback",
"subj_idx",
"split_by",
"trial",
"neural_reg",
]
]
# generate data trial-by-trial using the Intercept (a), regressor (neural_reg) and coefficient (neural) for decision threshold.
data, params = hddm.generate.gen_rand_data(
{
"a": a + neural * df.loc[0, "neural_reg"],
"t": t,
"v": df.loc[0, "sim_drift"],
"z": z,
},
subjs=1,
size=1,
)
df.loc[0, "response"] = data.response[0]
df.loc[0, "rt"] = data.rt[0]
if data.response[0] == 1.0:
df.loc[0, "feedback"] = df.loc[0, "rew_up"]
else:
df.loc[0, "feedback"] = df.loc[0, "rew_low"]
for i in range(1, n):
df.loc[i, "trial"] = i + 1
df.loc[i, "q_up"] = (
df.loc[i - 1, "q_up"] * (1 - df.loc[i - 1, "response"])
) + (
(df.loc[i - 1, "response"])
* (
df.loc[i - 1, "q_up"]
+ (alpha * (df.loc[i - 1, "rew_up"] - df.loc[i - 1, "q_up"]))
)
)
df.loc[i, "q_low"] = (
df.loc[i - 1, "q_low"] * (df.loc[i - 1, "response"])
) + (
(1 - df.loc[i - 1, "response"])
* (
df.loc[i - 1, "q_low"]
+ (alpha * (df.loc[i - 1, "rew_low"] - df.loc[i - 1, "q_low"]))
)
)
df.loc[i, "sim_drift"] = (df.loc[i, "q_up"] - df.loc[i, "q_low"]) * (scaler)
data, params = hddm.generate.gen_rand_data(
{
"a": a + neural * df.loc[i, "neural_reg"],
"t": t,
"v": df.loc[i, "sim_drift"],
"z": z,
},
subjs=1,
size=1,
)
df.loc[i, "response"] = data.response[0]
df.loc[i, "rt"] = data.rt[0]
if data.response[0] == 1.0:
df.loc[i, "feedback"] = df.loc[i, "rew_up"]
else:
df.loc[i, "feedback"] = df.loc[i, "rew_low"]
all_data.append(df)
all_data = pd.concat(all_data, axis=0)
all_data = all_data[
[
"q_up",
"q_low",
"sim_drift",
"response",
"rt",
"feedback",
"subj_idx",
"split_by",
"trial",
"neural_reg",
]
]
return all_data
# Create data with function defined above.
# This will create trial-by-trial values for decision threshold (a) by adding the coefficient neural (here set to 0.2)
# multiplied by a normalized regressor (neural_reg) to the 'Intercept' value of a (here set to 1)
data_neural = gen_rand_reg_rlddm_data(
a=1,
t=0.3,
scaler=2,
alpha=0.2,
neural=0.2,
size=100,
p_upper=0.7,
p_lower=0.3,
subjs=25,
)
data_neural["q_init"] = 0.5
data_neural.head()
# run a regressor model estimating the impact of 'neural' on decision threshold a. This should estimate the coefficient a_neural_reg to be 0.2
# to run the HDDMrlRegressor you need to include alpha
m_reg = hddm.HDDMrlRegressor(
data_neural, "a ~ neural_reg", include=["v", "a", "t", "alpha"]
)
m_reg.sample(1000, burn=250)
m_reg.print_stats()
HDDMrl also includes a module to run an RL-model that uses softmax to transform q-values to probability of choosing options associated with upper (response=1) or lower (response=0) boundary. To run this model you type hddm.Hrl instead of hddm.HDDMrl. The setup is the same as for HDDMrl, and for now, the model won't run if you don't include an rt-column. This will be fixed for a future version, but for now, if you don't have RTs you can just create an rt-column where you set all rts to e.g. 0.5. You can choose to estimate separate learning rates for positive and negative learning rate by setting dual=True (see here for more information). The model will by default estimate posterior distributions for the alpha and v parameters. The probability of choosing upper boundary is captured as:
$p_{up} =(e^{-2*z*d_t}-1)/ (e^{-2*d_t}-1)$,
where ${d_t}=q_{up_t}-q_{low}*v$ and z represents starting point (which for now is fixed to be 0.5).
This calculation is equivalent to soft-max transformation when z=0.5.
# run the model by calling hddm.Hrl (instead of hddm.HDDM for normal model and hddm.HDDMrl for rlddm-model)
m_rl = hddm.Hrl(data)
# set sample and burn-in
m_rl.sample(1500, burn=500, dbname="traces.db", db="pickle")
# print stats to get an overview of posterior distribution of estimated parameters
m_rl.print_stats()
Parameter estimates from the pure RL-model are a bit different compared to the RLDDM. This is to be expected as probability of choice in DDM is dependent both on the decsision threshold and the scaled difference in q-values, whereas the RL model only uses the scaled difference in q-values.
m_rl.plot_posteriors()
Fig. Mixing and autocorrelation looks good.
# estimate convergence
models = []
for i in range(3):
m = hddm.Hrl(data=data)
m.sample(1500, burn=500, dbname="traces.db", db="pickle")
models.append(m)
# get max gelman-statistic value. shouldn't be higher than 1.1
np.max(list(gelman_rubin(models).values()))
Convergence looks good, i.e. no parameters with gelman-rubin statistic > 1.1.
# Create a new model that has all traces concatenated
# of individual models.
m_rl = kabuki.utils.concat_models(models)
alpha, v = m_rl.nodes_db.node[["alpha", "v"]]
samples = {"alpha": alpha.trace(), "v": v.trace()}
samp = pd.DataFrame(data=samples)
def corrfunc(x, y, **kws):
r, _ = stats.pearsonr(x, y)
ax = plt.gca()
ax.annotate("r = {:.2f}".format(r), xy=(0.1, 0.9), xycoords=ax.transAxes)
g = sns.PairGrid(samp, palette=["red"])
g.map_upper(plt.scatter, s=10)
g.map_diag(sns.distplot, kde=False)
g.map_lower(sns.kdeplot, cmap="Blues_d")
g.map_lower(corrfunc)
Fig. The correlation in the posterior distribution for alpha and v/scaling is somewhat negative.
We can also do posterior predictive check on the RL-model by generating new data with hddm.generate.gen_rand_rl_data.
# create empty dataframe to store simulated data
sim_data = pd.DataFrame()
# create a column samp to be used to identify the simulated data sets
data["samp"] = 0
# load traces
traces = m_rl.get_traces()
# decide how many times to repeat simulation process. repeating this multiple times is generally recommended as it better captures the uncertainty in the posterior distribution, but will also take some time
for i in tqdm(range(1, 51)):
# randomly select a row in the traces to use for extracting parameter values
sample = np.random.randint(0, traces.shape[0] - 1)
# loop through all subjects in observed data
for s in data.subj_idx.unique():
# get number of trials for each condition.
size0 = len(
data[(data["subj_idx"] == s) & (data["split_by"] == 0)].trial.unique()
)
size1 = len(
data[(data["subj_idx"] == s) & (data["split_by"] == 1)].trial.unique()
)
size2 = len(
data[(data["subj_idx"] == s) & (data["split_by"] == 2)].trial.unique()
)
# set parameter values for simulation
scaler = traces.loc[sample, "v_subj." + str(s)]
alphaInv = traces.loc[sample, "alpha_subj." + str(s)]
# take inverse logit of estimated alpha
alpha = np.exp(alphaInv) / (1 + np.exp(alphaInv))
# simulate data for each condition changing only values of size, p_upper, p_lower and split_by between conditions.
sim_data0 = hddm.generate.gen_rand_rl_data(
scaler=scaler, alpha=alpha, size=size0, p_upper=0.8, p_lower=0.2, split_by=0
)
sim_data1 = hddm.generate.gen_rand_rl_data(
scaler=scaler, alpha=alpha, size=size1, p_upper=0.7, p_lower=0.3, split_by=1
)
sim_data2 = hddm.generate.gen_rand_rl_data(
scaler=scaler, alpha=alpha, size=size2, p_upper=0.6, p_lower=0.4, split_by=2
)
# append the conditions
sim_data0 = sim_data0.append([sim_data1, sim_data2], ignore_index=True)
# assign subj_idx
sim_data0["subj_idx"] = s
# identify that these are simulated data
sim_data0["type"] = "simulated"
# identify the simulated data
sim_data0["samp"] = i
# append data from each subject
sim_data = sim_data.append(sim_data0, ignore_index=True)
# combine observed and simulated data
ppc_rl_data = data[
["subj_idx", "response", "split_by", "trial", "feedback", "samp"]
].copy()
ppc_rl_data["type"] = "observed"
ppc_rl_sdata = sim_data[
["subj_idx", "response", "split_by", "trial", "feedback", "type", "samp"]
].copy()
ppc_rl_data = ppc_rl_data.append(ppc_rl_sdata)
# for practical reasons we only look at the first 40 trials for each subject in a given condition
plot_ppc_rl_data = ppc_rl_data[ppc_rl_data.trial < 41].copy()
# bin trials to for smoother estimate of response proportion across learning
plot_ppc_rl_data["bin_trial"] = pd.cut(
plot_ppc_rl_data.trial, 11, labels=np.linspace(0, 10, 11)
).astype("int64")
# calculate means for each sample
sums = (
plot_ppc_rl_data.groupby(["bin_trial", "split_by", "samp", "type"])
.mean()
.reset_index()
)
# calculate the overall mean response across samples
ppc_rl_sim = sums.groupby(["bin_trial", "split_by", "type"]).mean().reset_index()
# initiate columns that will have the upper and lower bound of the hpd
ppc_rl_sim["upper_hpd"] = 0
ppc_rl_sim["lower_hpd"] = 0
for i in range(0, ppc_rl_sim.shape[0]):
# calculate the hpd/hdi of the predicted mean responses across bin_trials
hdi = pymc.utils.hpd(
sums.response[
(sums["bin_trial"] == ppc_rl_sim.bin_trial[i])
& (sums["split_by"] == ppc_rl_sim.split_by[i])
& (sums["type"] == ppc_rl_sim.type[i])
],
alpha=0.1,
)
ppc_rl_sim.loc[i, "upper_hpd"] = hdi[1]
ppc_rl_sim.loc[i, "lower_hpd"] = hdi[0]
# calculate error term as the distance from upper bound to mean
ppc_rl_sim["up_err"] = ppc_rl_sim["upper_hpd"] - ppc_rl_sim["response"]
ppc_rl_sim["low_err"] = ppc_rl_sim["response"] - ppc_rl_sim["lower_hpd"]
ppc_rl_sim["model"] = "RL"
# plotting evolution of choice proportion for best option across learning for observed and simulated data. Compared for RL and RLDDM models, both with single learnign rate.
fig, axs = plt.subplots(figsize=(15, 5), nrows=1, ncols=3, sharex=True, sharey=True)
for i in range(0, 3):
ax = axs[i]
d_single = ppc_sim[(ppc_sim.split_by == i) & (ppc_sim.type == "simulated")]
# slightly move bin_trial to avoid overlap in errorbars
d_single["bin_trial"] += 0.2
ax.errorbar(
d_single.bin_trial,
d_single.response,
yerr=[d_single.low_err, d_single.up_err],
label="simulated_RLDDM",
color="orange",
)
ax = axs[i]
d_rl = ppc_rl_sim[(ppc_rl_sim.split_by == i) & (ppc_rl_sim.type == "simulated")]
ax.errorbar(
d_rl.bin_trial,
d_rl.response,
yerr=[d_rl.low_err, d_rl.up_err],
label="simulated_RL",
color="green",
)
ax = axs[i]
d = ppc_sim[(ppc_dual_sim.split_by == i) & (ppc_dual_sim.type == "observed")]
ax.plot(d.bin_trial, d.response, linewidth=3, label="observed")
ax.set_title("split_by = %i" % i, fontsize=20)
ax.set_ylabel("mean response")
ax.set_xlabel("trial")
plt.xlim(-0.5, 10.5)
plt.legend()
Fig. The predicted choice for the RL-model is very similar to what was predicted in the RLDDM. That is not surprising given that they use the same calculation to get the choice likelihood. The difference between them is instead that the RLDDM could potentially detect the unique contribution of the scaling/drift parameter and the decision threshold onto choice.
Another way to visualize this is to look at how the predicted choice misses on the observed across learning, i.e. predicted-observed. As for the other plots we see that the two methods are very similar.
# rl
error_prediction = (
plot_ppc_rl_data.groupby(["split_by", "type", "bin_trial"])["response"]
.mean()
.reset_index()
)
ep = error_prediction.pivot_table(
index=["split_by", "bin_trial"], columns="type", values="response"
).reset_index()
ep["diff"] = ep["simulated"] - ep["observed"]
ep["model"] = "RL"
# rlddm
error_prediction = (
plot_ppc_data.groupby(["split_by", "type", "bin_trial"])["response"]
.mean()
.reset_index()
)
ep_rlddm = error_prediction.pivot_table(
index=["split_by", "bin_trial"], columns="type", values="response"
).reset_index()
ep_rlddm["diff"] = ep_rlddm["simulated"] - ep_rlddm["observed"]
ep_rlddm["model"] = "RLDDM"
# combine
ep = ep.append(ep_rlddm)
# plot
g = sns.relplot(
x="bin_trial",
y="diff",
col="split_by",
hue="model",
kind="line",
ci=False,
data=ep,
palette="Set2_r",
)
g.map(plt.axhline, y=0, ls=":", c=".5")