LAN Tutorial
This tutorial is a rather comprehensive introduction to HDDM with
focus on the new LAN extension. The methods behind the new
HDDMnn()
, HDDMnnRegressor()
and HDDMnnStimcoding()
classes
can be found in our original dedicated
publication. These are
new featues. Please let us know on the HDDM forum and/or via github
reports regarding bugs or other limitations and we will do our best to
help as soon as we can.
Things to look out for:
Networks were trained over a fairly wide range of parameters which hopefully capture the scope of common empirical data. The networks will not accurately report likelihoods outside that range, so we explicitly limit the range of parameters that can be sampled from. If you find that your posterior samples reach and get stuck at the allowed parameter bounds (which you will see in the posterior plots), please notify us and we will do our best to provide improved networks over time.
You may encounter more print output than with standard HDDM. These are sanity checks and the verbosity will vanish progressively.
Section 0: Colab Prep (Optional)
Reminder
In the upper left menu click on Runtime, then Change runtime type and select GPU as hardware accelerator
INSTALLATION COLAB: INSTALL SUPPORT LIBRARIES
# Note: Usually colab has all other packages which we may use already installed
# The basic configuration of colabs does change over time, so you may have to add
# some install commands here if imports below don't work for package xyz
!pip install scikit-learn
!pip install cython
!pip install pymc
Requirement already satisfied: scikit-learn in /Users/afengler/opt/miniconda3/envs/hddmnn_tutorial/lib/python3.7/site-packages (0.24.2)
Requirement already satisfied: numpy>=1.13.3 in /Users/afengler/opt/miniconda3/envs/hddmnn_tutorial/lib/python3.7/site-packages (from scikit-learn) (1.19.1)
Requirement already satisfied: scipy>=0.19.1 in /Users/afengler/opt/miniconda3/envs/hddmnn_tutorial/lib/python3.7/site-packages (from scikit-learn) (1.7.2)
Requirement already satisfied: joblib>=0.11 in /Users/afengler/opt/miniconda3/envs/hddmnn_tutorial/lib/python3.7/site-packages (from scikit-learn) (1.0.1)
Requirement already satisfied: threadpoolctl>=2.0.0 in /Users/afengler/opt/miniconda3/envs/hddmnn_tutorial/lib/python3.7/site-packages (from scikit-learn) (2.1.0)
Requirement already satisfied: cython in /Users/afengler/opt/miniconda3/envs/hddmnn_tutorial/lib/python3.7/site-packages (0.29.24)
Requirement already satisfied: pymc in /Users/afengler/opt/miniconda3/envs/hddmnn_tutorial/lib/python3.7/site-packages (2.3.8)
INSTALLATION COLAB: INSTALL HDDM
Collecting git+https://github.com/hddm-devs/hddm
Cloning https://github.com/hddm-devs/hddm to /private/var/folders/gx/s43vynx550qbypcxm83fv56dzq4hgg/T/pip-req-build-xzqqwrcn
Running command git clone -q https://github.com/hddm-devs/hddm /private/var/folders/gx/s43vynx550qbypcxm83fv56dzq4hgg/T/pip-req-build-xzqqwrcn
Running command git submodule update --init --recursive -q
Building wheels for collected packages: HDDM
Building wheel for HDDM (setup.py) ... [?25l|
Imports
# MODULE IMPORTS ----
# warning settings
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
# Data management
import pandas as pd
import numpy as np
import pickle
# Plotting
import matplotlib.pyplot as plt
import matplotlib
import seaborn as sns
# Stats functionality
from statsmodels.distributions.empirical_distribution import ECDF
# HDDM
import hddm
from hddm.simulators.hddm_dataset_generators import simulator_h_c
Section 1: Model Info / Simulation / Basic Plotting
The main concern of this notebook is to present the extended
capabilities of the HDDM toolbox as a result of the new HDDMnn
classes.
Primarily we are interested in the additional models we can now be fit to data. So let’s take stock of the models that were added to standard HDDM.
2-Choice Models
ANGLE
A model with a linearly collapsing angle. Adds a parameter \(\theta\), which specifies the angle of the bound.
WEIBULL
A model that includes a collapsing bound parameterized as the scaled cdf of a Weibull distribution. This adds two parameters to the standard DDM, \(\alpha\) and \(\beta\).
LEVY
The Levy model is essentially a standard DDM where noise is not driven by a Gaussian distribution, but the noise process is now parameterized by the new parameter \(\alpha\), which interpolates between a Gausian \(\alpha = 2\) and a Cauchy (heavy tailed) \(\alpha = 1\).
ORNSTEIN
This model implements the 2-choice LCA, which includes a an inhibition / excitation parameter \(g\).
Find more details on these models in our companion paper.
3 / 4-Choice Models
NOTE
The addition of 3 choice and 4 choice models, comes with slightly more limited functionality as compared to 2 choice models. Specifically, not all plot-concepts currently standard in HDDM translate immediately to models with more choice options. We are trying to align this functionality going forward.
LCA (Leaky Competing Accumulator)
Please find the original description in this paper.
RACE
Race models simply take out the mutual and self-inhibition of LCAs.
ANGLE versions of LCA / RACE
Implements an linearly collapsing bound as above under the respective 2 choice models
1.1 Access Meta-Data
Let’s first take a look at some of the useful metadata we can use to set
up our models and simulators. If we type
hddm.simulators.model_config
, we get back a dictionary that stores a
bunch of information for each of the models that are currently
implemented in HDDM. It lists,
A
doc
string that gives some information about the status of the model as it pertains to it’s usability as well as some potential usage tips. Please read thedoc
string before using any of the new models.The parameter names under
params
,The parameter bounds that where used for training the network under
param_bounds
The boundary_function (
boundary
)Default parameter values (
params_default
).Slice sampler settings by parameter (
slice_widths
)Under
params_trans
you can choose parameters which will be logit transformed for sampling (order as inparams
)choices
determines valid choice options under the modelUnder
hddm_include
, it lists the parameters which we want to include when initializing our HDDM Model with one of the sequential sampling models available.
You won’t need most of these options if you are getting started, but they do provide you with useful information and a couple extra degrees of freedom when it comes to optimizing your sampler.
# List the models currently available
hddm.model_config.model_config.keys()
dict_keys(['ddm_hddm_base', 'full_ddm_hddm_base', 'ddm', 'angle', 'weibull', 'levy', 'full_ddm', 'ornstein', 'ddm_sdv', 'gamma_drift', 'gamma_drift_angle', 'ds_conflict_drift', 'ds_conflict_drift_angle', 'ddm_par2', 'ddm_par2_no_bias', 'ddm_par2_angle_no_bias', 'ddm_par2_weibull_no_bias', 'ddm_seq2', 'ddm_seq2_no_bias', 'ddm_seq2_angle_no_bias', 'ddm_seq2_weibull_no_bias', 'ddm_mic2_adj', 'ddm_mic2_adj_no_bias', 'ddm_mic2_adj_angle_no_bias', 'ddm_mic2_adj_weibull_no_bias', 'race_no_bias_3', 'race_no_bias_angle_3', 'race_no_bias_4', 'race_no_bias_angle_4', 'lca_no_bias_3', 'lca_no_bias_angle_3', 'lca_no_bias_4', 'lca_no_bias_angle_4', 'weibull_cdf', 'full_ddm2'])
NOTE
You find two kinds of extra models which were not mentioned in the model listing above:
Experimental models, which eventually will be fully documented (or dropped)
hddm_base
models are used predominantly with the basicHDDM()
classes. These models are not to be used with theHDDMnn()
classes.
Now taking a closer look at the angle
model
# Metadata
model = 'ddm'
n_samples = 1000
# Config for our current model
hddm.model_config.model_config[model]
{'doc': 'Basic DDM. Meant for use with the LAN extension. nNote that the boundaries here are coded as -a, and a in line with all other models meant for the LAN extension. nTo compare model fits between standard HDDM and HDDMnn when using the DDM model, multiply the boundary (a) parameter by 2. nWe recommend using standard HDDM if you are interested in the basic DDM, but you might want to use this for testing.', 'params': ['v', 'a', 'z', 't'], 'params_trans': [0, 0, 1, 0], 'params_std_upper': [1.5, 1.0, None, 1.0], 'param_bounds': [[-3.0, 0.3, 0.1, 0.001], [3.0, 2.5, 0.9, 2.0]], 'boundary': <function hddm.simulators.boundary_functions.constant(t=0)>, 'params_default': [0.0, 1.0, 0.5, 0.001], 'hddm_include': ['z'], 'choices': [-1, 1], 'slice_widths': {'v': 1.5, 'v_std': 1, 'a': 1, 'a_std': 1, 'z': 0.1, 'z_trans': 0.2, 't': 0.01, 't_std': 0.15}}
# Looking at the doc string before using the model
print(hddm.model_config.model_config[model]['doc'])
Basic DDM. Meant for use with the LAN extension.
Note that the boundaries here are coded as -a, and a in line with all other models meant for the LAN extension.
To compare model fits between standard HDDM and HDDMnn when using the DDM model, multiply the boundary (a) parameter by 2.
We recommend using standard HDDM if you are interested in the basic DDM, but you might want to use this for testing.
1.2 Generate Data
Let’s start by generating some data from the angle
model. For this
you have available the simulators
module, specifically we will start
with the simulator_h_c
function. If you are curious about all the
capabilities of this function, please check the help()
function for
it.
from hddm.simulators.hddm_dataset_generators import simulator_h_c
data, full_parameter_dict = simulator_h_c(n_subjects = 1,
n_trials_per_subject = n_samples,
model = model,
p_outlier = 0.00,
conditions = None,
depends_on = None,
regression_models = None,
regression_covariates = None,
group_only_regressors = False,
group_only = None,
fixed_at_default = None)
A quick look into what the simulator spits out (you can also read about
it in the docs). We get back a tuple
of two:
First, a DataFrame which holds a
rt
, aresponse
and asubj_idx
column as well as trial-by-trial ground truth parameters.Second a parameter dictionary which has parameter names in accordance with
HDDM()
trace names. This is useful for some of our plots.
data
rt | response | subj_idx | v | a | z | t | |
---|---|---|---|---|---|---|---|
0 | 1.956185 | 1.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
1 | 1.035191 | 0.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
2 | 1.004191 | 0.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
3 | 1.510186 | 0.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
4 | 1.164191 | 0.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
... | ... | ... | ... | ... | ... | ... | ... |
995 | 1.697184 | 0.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
996 | 1.520186 | 1.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
997 | 1.552186 | 0.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
998 | 1.038191 | 0.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
999 | 0.932191 | 1.0 | 0 | -0.481731 | 0.655642 | 0.439841 | 0.887191 |
1000 rows × 7 columns
# Here unspectacularly, parameter names are unchanged
# (single subject fits do not need any parameter name augmentation)
full_parameter_dict
{'v': -0.48173086489284433,
'a': 0.6556418306610691,
't': 0.8871907031605131,
'z': 0.4398408702789776}
1.2 First Plot
Now that we have our simulated data, we look to visualise it. Let’s look at a couple of plots that we can use for this purpose.
The HDDM.plotting
module includes the plot_from_data
function,
which allows you to plot subsets from a dataset, according to a grouping
specified by the groupby
argument.
The plot creates a matplotlib.axes
object for each subset, and you
can provide a function to manipulate this axes object. Some of these
axes manipulators are provided your you. Here we focus on the
_plot_func_model
axes manipulator supplied under the plot_func
argument.
Check out the arguments of plot_from_data
and _plot_func_model
using the help()
function. You have quite some freedom in styling
these plots.
We will refer to this plot as the model cartoon plot
.
The top histogram refers to the probability of choosing option \(1\) across time.
The bottom (upside-down) histogram refers to the probability of choosing option \(-1\) (may be coded as \(0\) as well) across time.
hddm.plotting.plot_from_data(df = data,
generative_model = model,
columns = 1,
groupby = ['subj_idx'],
figsize = (4, 3),
value_range = np.arange(0, 5, 0.1),
plot_func = hddm.plotting._plot_func_model,
**{'alpha': 1.,
'ylim': 3,
'add_data_rts': True,
'add_data_model': False})
plt.show()
subj_idx(0)
If we set add_model = True
, this will add a cartoon of the model on
top of the histograms.
CAUTION
This model cartoon plot
will only work for 2-choice models for
now.
Moreover, often useful for illustration purposes, we can include a bunch
of simulations trajectories into the model plot (note the corresponding
arguments). Common to all models currently included is their conceptual
reliance on there particle trajectories. Reaction times and choices are
simulated as boundary crossings of these particles. If you don’t want
to include these trajectories, just set show_trajectories = False
.
hddm.plotting.plot_from_data(df = data,
generative_model = model,
columns = 1,
groupby = ['subj_idx'],
figsize = (4, 3),
value_range = np.arange(0, 5, 0.1),
plot_func = hddm.plotting._plot_func_model,
**{'alpha': 1.,
'ylim': 3,
'add_data_rts': True,
'add_data_model': True})
plt.show()
subj_idx(0)
If you are interested, you can use this plot to investigate the behavior of models across different parameters setups.
Section 2: Single Subject (or collapsed) Data
Now, we try to fit these models to data! Let’s start with an simple dataset. In other words, we have one single participant who provides \(n\) datatpoints (reaction times and choices) from some two alternative forced choice task paradigm.
Note
In this demo we fit to simulated data. This serves as a template, and you can easily adapt it to your needs.
# Metadata
nmcmc = 1500
model = 'angle'
n_samples = 1000
includes = hddm.model_config.model_config[model]['hddm_include']
Note
When defining includes
, you can also pick only as subset of the
parameters suggested under hddm.model_config.model_config
.
# Generate some simulatred data
from hddm.simulators.hddm_dataset_generators import simulator_h_c
data, full_parameter_dict = simulator_h_c(n_subjects = 1,
n_trials_per_subject = n_samples,
model = model,
p_outlier = 0.00,
conditions = None,
depends_on = None,
regression_models = None,
regression_covariates = None, # need this to make initial covariate matrix from which to use dmatrix (patsy)
group_only_regressors = False,
group_only = None,
fixed_at_default = None)
data
rt | response | subj_idx | v | a | z | t | theta | |
---|---|---|---|---|---|---|---|---|
0 | 2.096904 | 0.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
1 | 2.154903 | 0.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
2 | 1.862907 | 0.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
3 | 1.847907 | 0.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
4 | 1.927906 | 0.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
... | ... | ... | ... | ... | ... | ... | ... | ... |
995 | 2.260902 | 1.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
996 | 1.895906 | 0.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
997 | 1.782908 | 0.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
998 | 1.864907 | 0.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
999 | 1.812907 | 0.0 | 0 | -1.688219 | 1.945201 | 0.543195 | 1.33591 | 0.875893 |
1000 rows × 8 columns
# Define the HDDM model
hddmnn_model = hddm.HDDMnn(data,
informative = False,
include = includes,
p_outlier = 0.01,
w_outlier = 0.1,
model = model,)
Supplied model_config specifies params_std_upper for z as None.
Changed to 10
# Sample
hddmnn_model.sample(nmcmc,
burn = 500)
[-----------------100%-----------------] 1500 of 1500 complete in 103.2 sec
<pymc.MCMC.MCMC at 0x141b8e410>
2.1 Visualization
The plot_caterpillar()
function below displays parameterwise,
as a blue tick-mark the ground truth.
as a thin black line the \(1 - 99\) percentile range of the posterior distribution
as a thick black line the \(5-95\) percentile range of the posterior distribution
Again use the help()
function to learn more.
# Caterpillar Plot: (Parameters recovered ok?)
hddm.plotting.plot_caterpillar(hddm_model = hddmnn_model,
ground_truth_parameter_dict = full_parameter_dict,
figsize = (8, 5),
columns = 3)
plt.show()
2.1.1 Posterior Predictive (via model cartoon plot
)
Another way to examine whether or not our recovery was satisfactory is to perform posterior predictive checks. Essentially, we are looking to simulate datasets from the trace and check whether it aligns with the ground truth participant data. This answers the question of whether or not these parameters that you recovered can actually reproduce the data.
Use the plot_posterior_predictive()
function in the plotting
module for this. It is structured just like the plot_from_data()
function, but instead of providing a dataset, you supply a hddm
model.
Use the help()
function to check out all the functionality.
hddm.plotting.plot_posterior_predictive(model = hddmnn_model,
columns = 1,
groupby = ['subj_idx'],
figsize = (6, 4),
value_range = np.arange(0, 5, 0.1),
plot_func = hddm.plotting._plot_func_model,
parameter_recovery_mode = True,
**{'alpha': 0.01,
'ylim': 3,
'samples': 200})
plt.show()
passing
A small note on convergence:
Note that the MCMC algorithm requires the chain to converge. There are
many heuristics that help you identifying problems with convergence,
such as the trace plot, auto correlation plot, and marginal posterior
histogram. In the trace plots, there might be a problem if you see large
jumps. In the autocorrelation plot, there might be a problem if it does
not drop rapidly. The HDDMnn()
classes support the computation of
the Gelman-Rubin, r-hat statistic, as you would with any hddm
model. Generally, by extracting the traces, you are free to compute any
convergence statistics you want of course.
# TAKING A LOOK AT THE POSTERIOR TRACES
hddmnn_model.plot_posteriors(hddm.simulators.model_config[model]['params'])
plt.show()
Plotting v
Plotting a
Plotting z
Plotting t
Plotting theta
hddm.plotting.plot_posterior_pair(hddmnn_model, save = False,
parameter_recovery_mode = True,
samples = 500,
figsize = (6, 6))
Section 3: Hierarchical Models
The ‘h’ in hddm
stands for hierarchical, so let’s do it! If we have
data from multiple participants and we assume that the parameters of
single participants are drawn from respective group or global
distributions, we can model this explicitly in hddm
by specifying
is_group_model = True
.
Implicitly we are fitting a model of the following kind,
where (let’s say for the angle model),
\(\theta_j = \{v_j, a_j, z_j, t_j, \theta_j \}\), are the model parameters for subject j.
\(\theta_g = \{v_g^{\mu}, a_g^{\mu}, z_g^{\mu}, t_g^{\mu}, \theta_g^{\mu}, v_g^{\sigma}, a_g^{\sigma}, z_g^{\sigma}, t_g^{\sigma}, \theta_g^{\sigma} \}\) (scary, but for completeness), are the mean and variance parameters for our group level normal distributions, and \(\{ \theta_h \}\) are fixed hyperparameters.
\(x_i^j = \{rt_i^j, c_i^j \}\), are the choice and reaction time of subject j during trial i.
In words, the right hand side of the equation tells us that we have a global parameter distribution with certain means and variances for each parameter (we want to figure these means and variances out), from which the subject level parameters are drawn and finally subject level datapoints follow the likelihood distribution of our ddm / angle / weibull / you name it mdoels.
# Metadata
nmcmc = 1000
model = 'angle'
n_trials_per_subject = 200
n_subjects = 10
# test regressors only False
# add p_outliers to the generator !
data, full_parameter_dict = simulator_h_c(data = None,
n_subjects = n_subjects,
n_trials_per_subject = n_trials_per_subject,
model = model,
p_outlier = 0.00,
conditions = None,
depends_on = None,
regression_models = None,
regression_covariates = None,
group_only_regressors = False,
group_only = None,
fixed_at_default = None)
hddmnn_model = hddm.HDDMnn(data,
model = model,
informative = False,
is_group_model = True,
include = hddm.simulators.model_config[model]['hddm_include'],
p_outlier = 0.0)
{'v': 1.5, 'v_std': 1, 'a': 1, 'a_std': 1, 'z': 0.1, 'z_trans': 0.2, 't': 0.01, 't_std': 0.15, 'theta': 0.1, 'theta_std': 0.2}
Supplied model_config specifies params_std_upper for z as None.
Changed to 10
hddmnn_model.sample(nmcmc,
burn = 100) # if you want to save the model specify extra arguments --> dbname='traces.db', db='pickle'. # hddmnn_model.save('test_model')
[-----------------100%-----------------] 1000 of 1000 complete in 339.0 sec
<pymc.MCMC.MCMC at 0x14bb81390>
# Caterpillar Plot: (Parameters recovered ok?)
hddm.plotting.plot_caterpillar(hddm_model = hddmnn_model,
ground_truth_parameter_dict = full_parameter_dict,
figsize = (8, 5),
columns = 3)
plt.show()
hddm.plotting.plot_posterior_predictive(model = hddmnn_model,
columns = 3,
figsize = (10, 7),
groupby = ['subj_idx'],
value_range = np.arange(0, 5, 0.1),
plot_func = hddm.plotting._plot_func_model,
parameter_recovery_mode = True,
**{'alpha': 0.01,
'ylim': 3,
'add_posterior_mean_rts': True,
'add_posterior_mean_model': True,
'add_posterior_uncertainty_rts': False,
'add_posterior_uncertainty_model': False,
'samples': 200,
'legend_fontsize': 7.})
Section 4: Parameter varies by Condition
An important aspect of these posterior analysis, is the consideration of experiment design. We may have an experiment in which subject are exposed to a variety of conditions, such as for example different degrees of difficulty of the same task
It is often reasonable to assume that all but the conceptually relevant parameters are common across conditions.
As a by-product, such experiment designs can help us with the recovery of the constant parameters, by probing those static aspects of the model across varying kinds of datasets (driven by targeted manipulation of variable aspects of the model).
Implicitly we fit the following kind of model,
Where \(\theta_c\) is the condition dependent part of the parameter space, and \(\theta\) forms the portion of parameters which remain constant across condtions.
To give a more concrete example involving the weibull model, consider a dataset for a single participant, who went through four conditions of an experiment. Think of the conditions as manipulating the payoff structure of the experiment to incentivize / disincentivize accuracy in favor of speed. We operationalize this by treating the \(a\) parameter, the initial boundary separation, as affected by the manipulation, while the rest of the parameters are constant across all experiment conditions.
The resulting model would be of the form,
# Metadata
nmcmc = 1000
model = 'angle'
n_trials_per_subject = 500
# We allow the boundary conditions to vary
depends_on = {'a': ['c_one']}
# They will depend on a fictious column 'c_one' that specifies
# levels / conditions
conditions = {'c_one': ['low', 'medium', 'high']}
data, full_parameter_dict = simulator_h_c(n_subjects = 1,
n_trials_per_subject = n_trials_per_subject,
model = model,
p_outlier = 0.00,
conditions = conditions,
depends_on = depends_on,
regression_models = None,
regression_covariates = None,
group_only_regressors = False,
group_only = None,
fixed_at_default = None)
depends_on is: {'a': ['c_one']}
# Let's check the resulting parameter vector
full_parameter_dict
{'theta': 0.7406253194726012,
'v': 1.464554358239174,
'z': 0.6206249211841304,
't': 1.534252965986117,
'a(high)': 1.0519165572885651,
'a(low)': 1.2561997135872933,
'a(medium)': 0.9265856569938499}
# Make HDDM Model
hddmnn_model = hddm.HDDMnn(data,
model = model,
informative = False,
include = hddm.simulators.model_config[model]['hddm_include'],
p_outlier = 0.0,
is_group_model = False,
depends_on = depends_on)
# Sample
hddmnn_model.sample(nmcmc, burn = 100)
[-----------------100%-----------------] 1001 of 1000 complete in 129.4 sec
<pymc.MCMC.MCMC at 0x14f44c690>
# Caterpillar Plot: (Parameters recovered ok?)
hddm.plotting.plot_caterpillar(hddm_model = hddmnn_model,
ground_truth_parameter_dict = full_parameter_dict,
figsize = (8, 5),
columns = 3)
plt.show()
hddm.plotting.plot_posterior_predictive(model = hddmnn_model,
columns = 1,
groupby = ['subj_idx'],
figsize = (4, 4),
value_range = np.arange(0, 5, 0.1),
plot_func = hddm.plotting._plot_func_model,
parameter_recovery_mode = True,
**{'alpha': 0.01,
'ylim': 3,
'add_posterior_uncertainty_rts': True,
'add_posterior_uncertainty_model': True,
'samples': 200})
plt.show()
passing
passing
passing
4.1 Combine Hierarchical and Condition data
# Metadata
nmcmc = 1500
model = 'angle'
n_subjects = 5
n_trials_per_subject = 500
data, full_parameter_dict = simulator_h_c(n_subjects = n_subjects,
n_trials_per_subject = n_trials_per_subject,
model = model,
p_outlier = 0.00,
conditions = {'c_one': ['low', 'medium', 'high']}, #, 'c_three': ['low', 'medium', 'high']},
depends_on = {'v': ['c_one']}, # 'theta': ['c_two']}, # 'theta': ['c_two']}, #regression_models = None, #
regression_models = None, #regression_covariates = None,
regression_covariates = None, # need this to make initial covariate matrix from which to use dmatrix (patsy)
group_only_regressors = False,
group_only = None,
fixed_at_default = None)
depends_on is: {'v': ['c_one']}
# Make HDDM Model
hddmnn_model = hddm.HDDMnn(data,
model = model,
informative = False,
include = hddm.simulators.model_config[model]['hddm_include'],
p_outlier = 0.0,
is_group_model = True,
depends_on = {'v': 'c_one'})
hddmnn_model.sample(nmcmc, burn = 100)
[-----------------100%-----------------] 1500 of 1500 complete in 919.0 sec
<pymc.MCMC.MCMC at 0x14e324a90>
# Caterpillar Plot: (Parameters recovered ok?)
hddm.plotting.plot_caterpillar(hddm_model = hddmnn_model,
ground_truth_parameter_dict = full_parameter_dict,
figsize = (8, 8),
columns = 3)
plt.show()
hddm.plotting.plot_posterior_predictive(model = hddmnn_model,
columns = 2, # groupby = ['subj_idx'],
figsize = (8, 6),
value_range = np.arange(1, 2.5, 0.1),
plot_func = hddm.plotting._plot_func_model,
parameter_recovery_mode = True,
**{'alpha': 0.01,
'ylim': 3,
'add_posterior_uncertainty_rts': True,
'add_posterior_uncertainty_model': True,
'samples': 200,
'legend_fontsize': 7})
plt.show()
passing
passing
passing
passing
passing
passing
passing
passing
passing
passing
passing
passing
passing
passing
passing
Section 5: Regressors
This section provides a simple working example using the Neural Networks with the Regression backend. The regression back-end allows linking parameters to trial-by-trial covariates via a (general) linear model.
# Metadata
nmcmc = 1000
model = 'angle'
n_samples_by_subject = 500
from hddm.simulators.hddm_dataset_generators import simulator_h_c
data, full_parameter_dict = simulator_h_c(n_subjects = 5,
n_samples_by_subject = n_samples_by_subject,
model = model,
p_outlier = 0.00,
conditions = None,
depends_on = None,
regression_models = ['t ~ 1 + covariate_name', 'v ~ 1 + covariate_name'],
regression_covariates = {'covariate_name': {'type': 'continuous', 'range': (0, 1)}},
group_only_regressors = False,
group_only = None,
fixed_at_default = None)
# Set up the regressor a regressor:
reg_model_v = {'model': 'v ~ 1 + covariate_name', 'link_func': lambda x: x}
reg_model_t = {'model': 't ~ 1 + covariate_name', 'link_func': lambda x: x}
reg_descr = [reg_model_t, reg_model_v]
# Make HDDM model
hddmnn_reg = hddm.HDDMnnRegressor(data,
reg_descr,
include = hddm.simulators.model_config[model]['hddm_include'],
model = model,
informative = False,
p_outlier = 0.0)
Supplied model_config specifies params_std_upper for z as None.
Changed to 10
# Sample
hddmnn_reg.sample(nmcmc, burn = 100)
[-----------------100%-----------------] 1001 of 1000 complete in 369.4 sec
<pymc.MCMC.MCMC at 0x149a07890>
# Caterpillar Plot: (Parameters recovered ok?)
hddm.plotting.plot_caterpillar(hddm_model = hddmnn_reg,
ground_truth_parameter_dict = full_parameter_dict,
figsize = (8, 8),
columns = 3)
plt.show()
Section 6: Stim Coding
You can read more about stimulus coding in the documentation.
Here just an example.
# Metadata
nmcmc = 300
model = 'ddm'
n_samples_by_condition = 500
split_param = 'v'
sim_data_stimcoding, parameter_dict = hddm.simulators.simulator_stimcoding(model = model,
split_by = split_param,
drift_criterion = 0.3,
n_trials_per_condition = 500)
sim_data_stimcoding
rt | response | stim | v | a | z | t | subj_idx | |
---|---|---|---|---|---|---|---|---|
0 | 3.190470 | 1.0 | 1 | 0.834704 | 2.426857 | 0.417932 | 1.507448 | none |
1 | 3.942454 | 1.0 | 1 | 0.834704 | 2.426857 | 0.417932 | 1.507448 | none |
2 | 4.186436 | 1.0 | 1 | 0.834704 | 2.426857 | 0.417932 | 1.507448 | none |
3 | 2.205442 | 1.0 | 1 | 0.834704 | 2.426857 | 0.417932 | 1.507448 | none |
4 | 4.669401 | 1.0 | 1 | 0.834704 | 2.426857 | 0.417932 | 1.507448 | none |
... | ... | ... | ... | ... | ... | ... | ... | ... |
495 | 11.207737 | 0.0 | 2 | -0.234704 | 2.426857 | 0.417932 | 1.507448 | none |
496 | 10.334385 | 1.0 | 2 | -0.234704 | 2.426857 | 0.417932 | 1.507448 | none |
497 | 7.077227 | 0.0 | 2 | -0.234704 | 2.426857 | 0.417932 | 1.507448 | none |
498 | 8.740107 | 1.0 | 2 | -0.234704 | 2.426857 | 0.417932 | 1.507448 | none |
499 | 2.621444 | 0.0 | 2 | -0.234704 | 2.426857 | 0.417932 | 1.507448 | none |
1000 rows × 8 columns
parameter_dict
{'v': -0.5347036843723503,
'a': 2.426856838254428,
'z': 0.4179319892615798,
't': 1.5074477307220377,
'dc': 0.3}
hddmnn_model = hddm.HDDMnnStimCoding(sim_data_stimcoding,
include = hddm.simulators.model_config[model]['hddm_include'],
model = model,
stim_col = 'stim',
p_outlier = 0.0,
split_param = split_param,
informative = False,
drift_criterion = True)
hddmnn_model.sample(nmcmc, burn = 100)
[-----------------100%-----------------] 300 of 300 complete in 32.4 sec
<pymc.MCMC.MCMC at 0x14e56a850>
hddmnn_model.gen_stats()
mean | std | 2.5q | 25q | 50q | 75q | 97.5q | mc err | |
---|---|---|---|---|---|---|---|---|
v | -0.539954 | 0.0155719 | -0.572298 | -0.549469 | -0.540782 | -0.528279 | -0.508065 | 0.00112267 |
a | 2.49136 | 0.00895635 | 2.47002 | 2.48866 | 2.49422 | 2.49799 | 2.49988 | 0.000770613 |
z | 0.4031 | 0.0118855 | 0.37898 | 0.397112 | 0.40288 | 0.409666 | 0.431941 | 0.000958058 |
t | 1.48852 | 0.035718 | 1.41497 | 1.46917 | 1.48828 | 1.51355 | 1.56092 | 0.00286112 |
dc | 0.348321 | 0.0202259 | 0.30703 | 0.334826 | 0.349703 | 0.361927 | 0.386422 | 0.00170049 |
# Caterpillar Plot: (Parameters recovered ok?)
hddm.plotting.plot_caterpillar(hddm_model = hddmnn_model,
ground_truth_parameter_dict = parameter_dict,
figsize = (8, 5),
columns = 3)
plt.show()
NOTE:
The hddm.plotting.plot_posterior_predictive()
does not yet accept
stimcoding data. This will be updated as soon as possible.
Section 7: Model Recovery
A crucial exercise in statistical modeling concern model comparison.
We are going to look at model recovery, in this section: Attempt to recover which model generated a given dataset from a set of candidate models.
For the little model recovery study we conduct here, we generate data from the weibull model and fit the data once each to the weibull, angle and ddm models.
We inspect the fits visually and then use the DIC (Deviance information criterion, lower is better :)), to check if we can recover the true model.
# Metadata
model = 'weibull'
n_samples = 300
# test regressors only False
# add p_outliers to the generator !
data, full_parameter_dict = simulator_h_c(n_subjects = 1,
n_samples_by_subject = n_samples,
model = model,
p_outlier = 0.00,
conditions = None,
depends_on = None,
regression_models = None,
regression_covariates = None,
group_only_regressors = False,
group_only = None,
fixed_at_default = None)
data
rt | response | subj_idx | v | a | z | t | alpha | beta | |
---|---|---|---|---|---|---|---|---|---|
0 | 4.204582 | 0.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
1 | 4.269577 | 1.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
2 | 4.404568 | 0.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
3 | 2.960620 | 1.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
4 | 2.223596 | 1.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
95 | 2.304595 | 1.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
96 | 3.067625 | 1.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
97 | 2.379594 | 1.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
98 | 3.991597 | 1.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
99 | 2.904617 | 0.0 | 0 | 0.246969 | 1.470066 | 0.451724 | 1.397603 | 3.268501 | 4.603728 |
100 rows × 9 columns
# Now we fit for each model:
hddmnn_model_weibull = hddm.HDDMnn(data,
informative = False,
model = 'weibull',
p_outlier = 0.0,
include = hddm.simulators.model_config['weibull_cdf']['hddm_include'],
is_group_model = False)
hddmnn_model_angle = hddm.HDDMnn(data,
model = 'angle',
informative = False,
p_outlier = 0.0,
include = hddm.simulators.model_config['angle']['hddm_include'],
is_group_model = False)
hddmnn_model_ddm = hddm.HDDMnn(data,
informative = False,
model = 'ddm',
p_outlier = 0.0,
include = hddm.simulators.model_config['ddm']['hddm_include'],
is_group_model = False)
nmcmc = 1000
hddmnn_model_weibull.sample(nmcmc,
burn = 200)
hddmnn_model_angle.sample(nmcmc,
burn = 200)
hddmnn_model_ddm.sample(nmcmc,
burn = 200)
[-----------------100%-----------------] 1000 of 1000 complete in 23.0 sec
<pymc.MCMC.MCMC at 0x14d606c90>
7.1 Checking Model Fits Visually
Posterior Predictive: Do the ‘Posterior Models’ also make sense?
# WEIBULL
hddm.plotting.plot_posterior_predictive(model = hddmnn_model_weibull,
columns = 1,
groupby = ['subj_idx'],
figsize = (4, 4),
value_range = np.arange(0, 5, 0.1),
plot_func = hddm.plotting._plot_func_model,
parameter_recovery_mode = True,
**{'alpha': 0.01,
'ylim': 3,
'add_posterior_uncertainty_model': True,
'add_posterior_uncertainty_rts': False,
'add_posterior_mean_rts': True,
'samples': 200})
plt.show()
passing
# ANGLE
hddm.plotting.plot_posterior_predictive(model = hddmnn_model_angle,
columns = 1,
groupby = ['subj_idx'],
figsize = (4, 4),
value_range = np.arange(0, 5, 0.1),
plot_func = hddm.plotting._plot_func_model,
parameter_recovery_mode = False,
**{'alpha': 0.01,
'ylim': 3,
'add_posterior_uncertainty_model': True,
'add_posterior_uncertainty_rts': False,
'add_posterior_mean_rts': True,
'samples': 200})
plt.show()
# DDM
hddm.plotting.plot_posterior_predictive(model = hddmnn_model_ddm,
columns = 1,
groupby = ['subj_idx'],
figsize = (4, 4),
value_range = np.arange(0, 5, 0.1),
plot_func = hddm.plotting._plot_func_model,
parameter_recovery_mode = False,
**{'alpha': 0.01,
'ylim': 3,
'add_posterior_uncertainty_model': True,
'add_posterior_uncertainty_rts': False,
'add_posterior_mean_rts': True,
'samples': 200})
plt.show()
7.2 Comparing DIC’s
hddmnn_model_weibull.dic
414.65114936828616
hddmnn_model_angle.dic
415.8001557922363
hddmnn_model_ddm.dic
418.04479835510256
Fingers crossed (this was a random run after all), the DIC usually gives us a result that conforms with the intuition we get from looking at the model plots.
Section 8: Real Data!
# Metadata
nmcmc = 1000
burn = 500
model = 'angle'
8.1 Load and Pre-process dataset
# Load one of the datasets shipping with HDDM
cav_data = hddm.load_csv(hddm.__path__[0] + '/examples/cavanagh_theta_nn.csv')
cav_data
subj_idx | stim | rt | response | theta | dbs | conf | |
---|---|---|---|---|---|---|---|
0 | 0 | LL | 1.210 | 1.0 | 0.656275 | 1 | HC |
1 | 0 | WL | 1.630 | 1.0 | -0.327889 | 1 | LC |
2 | 0 | WW | 1.030 | 1.0 | -0.480285 | 1 | HC |
3 | 0 | WL | 2.770 | 1.0 | 1.927427 | 1 | LC |
4 | 0 | WW | 1.140 | 0.0 | -0.213236 | 1 | HC |
... | ... | ... | ... | ... | ... | ... | ... |
3983 | 13 | LL | 1.450 | 0.0 | -1.237166 | 0 | HC |
3984 | 13 | WL | 0.711 | 1.0 | -0.377450 | 0 | LC |
3985 | 13 | WL | 0.784 | 1.0 | -0.694194 | 0 | LC |
3986 | 13 | LL | 2.350 | 0.0 | -0.546536 | 0 | HC |
3987 | 13 | WW | 1.250 | 1.0 | 0.752388 | 0 | HC |
3988 rows × 7 columns
8.2 Basic Condition Split Model
hddmnn_model_cav = hddm.HDDMnn(cav_data,
model = model,
informative = False,
include = hddm.simulators.model_config[model]['hddm_include'],
p_outlier = 0.05,
is_group_model = False,
depends_on = {'v': 'stim'})
hddmnn_model_cav.sample(nmcmc, burn = burn)
[-----------------100%-----------------] 1000 of 1000 complete in 243.4 sec
<pymc.MCMC.MCMC at 0x154a35350>
hddm.plotting.plot_posterior_predictive(model = hddmnn_model_cav,
columns = 1,
figsize = (4, 4),
value_range = np.arange(0, 5, 0.1),
plot_func = hddm.plotting._plot_func_model,
parameter_recovery_mode = False,
**{'alpha': 0.01,
'ylim': 3,
'add_posterior_uncertainty_model': True,
'add_posterior_uncertainty_rts': False,
'add_posterior_mean_rts': True,
'samples': 200})
plt.show()
8.3 Basic Hierarchical Model
hddmnn_model_cav = hddm.HDDMnn(cav_data,
model = model,
informative = False,
include = hddm.simulators.model_config[model]['hddm_include'],
is_group_model = True,
p_outlier = 0.05)
hddmnn_model_cav.sample(nmcmc, burn = burn)
[-----------------100%-----------------] 1000 of 1000 complete in 471.3 sec
<pymc.MCMC.MCMC at 0x14cca5090>
# Caterpillar Plot: (Parameters recovered ok?)
hddm.plotting.plot_caterpillar(hddm_model = hddmnn_model_cav,
figsize = (8, 8),
columns = 3)
plt.show()
hddm.plotting.plot_posterior_predictive(model = hddmnn_model_cav,
columns = 3,
figsize = (10, 10),
value_range = np.arange(0, 5, 0.1),
plot_func = hddm.plotting._plot_func_model,
parameter_recovery_mode = False,
**{'alpha': 0.01,
'ylim': 3,
'add_posterior_uncertainty_model': True,
'add_posterior_uncertainty_rts': False,
'add_posterior_mean_rts': True,
'samples': 200,
'legend_fontsize': 7,
'subplots_adjust': {'top': 0.9, 'hspace': 0.3, 'wspace': 0.3}})
plt.show()
Note
This is just an example. The angle model might not be the best choice here, and we are moreover ignoring the supplied conditions.
Section 9: Accessing the Neural Network Directly
The network_inspectors
module allows you to inspect the LANs
directly.
9.1 Direct access to batch predictions
You can use the hddm.network_inspectors.get_torch_mlp()
function to
access network predictions.
model = 'angle'
lan_angle = hddm.network_inspectors.get_torch_mlp(model = model)
Let’s predict some likelihoods !
# Make some random parameter set
parameter_df = hddm.simulators.make_parameter_vectors_nn(model = model,
param_dict = None,
n_parameter_vectors = 1)
parameter_matrix = np.tile(np.squeeze(parameter_df.values), (200, 1))
# Initialize network input
network_input = np.zeros((parameter_matrix.shape[0], parameter_matrix.shape[1] + 2)) # Note the + 2 on the right --> we append the parameter vectors with reaction times (+1 columns) and choices (+1 columns)
# Add reaction times
network_input[:, -2] = np.linspace(0, 3, parameter_matrix.shape[0])
# Add choices
network_input[:, -1] = np.repeat(np.random.choice([-1, 1]), parameter_matrix.shape[0])
# Convert to float
network_input = network_input.astype(np.float32)
# Show example output
print(lan_angle(network_input)[:10]) # printing the first 10 outputs
print(lan_angle(network_input).shape) # original shape of output
[[-2.9323568]
[ 2.078088 ]
[ 0.4104141]
[-0.5943402]
[-1.1136726]
[-1.6901499]
[-2.3512228]
[-3.080151 ]
[-3.8215086]
[-4.4257374]]
(200, 1)
9.2 Plotting Utilities
HDDM provides two plotting function to investigate the network outputs
directly. The kde_vs_lan_likelihoods()
plot and the
lan_manifold()
plot.
9.2.1 kde_vs_lan_likelihoods()
The kde_vs_lan_likelihoods()
plot allows you to check the
likelihoods produced by a LAN against Kernel Density Estimates (KDEs)
from model simulations. You can supply a panda DataFrame
that holds
parameter vectors as rows.
# Make some parameters
parameter_df = hddm.simulators.make_parameter_vectors_nn(model = model,
param_dict = None,
n_parameter_vectors = 10)
parameter_df
v | a | z | t | theta | |
---|---|---|---|---|---|
0 | 2.149626 | 1.684902 | 0.232222 | 0.641663 | -0.070030 |
1 | 1.817911 | 0.776330 | 0.535083 | 0.006625 | 1.069452 |
2 | -0.908428 | 0.654107 | 0.301445 | 1.560911 | 0.396448 |
3 | -0.022136 | 1.140235 | 0.479664 | 0.757727 | 1.316409 |
4 | 2.281230 | 0.366558 | 0.409224 | 1.908211 | 1.059872 |
5 | 1.067632 | 1.228020 | 0.337573 | 1.447155 | 0.083665 |
6 | 2.022131 | 1.254037 | 0.262336 | 0.416462 | 0.512724 |
7 | -1.974657 | 0.793536 | 0.791707 | 0.591319 | 1.036441 |
8 | -2.002436 | 1.382722 | 0.442411 | 0.074192 | 0.356522 |
9 | -2.757462 | 0.402900 | 0.738999 | 0.755093 | 1.334423 |
hddm.network_inspectors.kde_vs_lan_likelihoods(parameter_df = parameter_df,
model = model,
cols = 3,
n_samples = 2000,
n_reps = 10,
show = True)
1 of 10
2 of 10
3 of 10
4 of 10
5 of 10
6 of 10
7 of 10
8 of 10
9 of 10
10 of 10
9.2.2 lan_manifold()
Lastly, you can use the lan_manifold()
plot to investigate the LAN
likelihoods over a range of parameters.
The idea is to use a base parameter vector and vary one of the parameters in a prespecificed range.
This plot can be informative if you would like to understand better how a parameter affects model behavior.
# Make some parameters
parameter_df = hddm.simulators.make_parameter_vectors_nn(model = model,
param_dict = None,
n_parameter_vectors = 1)
parameter_df
v | a | z | t | theta | |
---|---|---|---|---|---|
0 | -2.218164 | 0.889863 | 0.254979 | 0.707028 | 0.040745 |
# Now plotting
hddm.network_inspectors.lan_manifold(parameter_df = parameter_df,
vary_dict = {'v': np.linspace(-2, 2, 20)},
model = model,
n_rt_steps = 300,
fig_scale = 1.0,
max_rt = 5,
save = True,
show = True)
Using only the first row of the supplied parameter array !
Hopefully this tutorial proves as a useful starting point for your application.