How to fit a model

To fit a model, there are 3 main things to specify:

1. Data: The data, which should be in the form of a pandas DataFrame.

Different model classes might require different columns in the data. You should check in the API Reference of each model class (or using model.fit?) what the required data columns are.

2. The priors (optional): You can decide whether to use the default priors (which you can see after initializing the model) or whether you want to change the mean or SD of the prior or hyper-prior distributions. Whether you changed the priors or not, they are always printed out when the model starts fitting.

3. Sampling parameters: The sampling parameters (number of chains, iterations, warmups, thinning, etc.) are the arguments to the pystan.StanModel.sampling() function, and we simply refer to the pystan documentation for a better overview.

Additional learning parameters: While all sequential sampling models (DDM and race models) without a learning component only require a data argument, all models with a learning components (RL models, RLDDMs, and RL+race models) also require a K argument, which is the total number of different options in a learning block (note that this can be different from the number of options presented in each trial), and initial_value_learning, which is the initial Q value (before learning).

import rlssm

Non-learning example (non-hierarchical, simulated data)

model_ddm = rlssm.DDModel(hierarchical_levels=1)

Using cached StanModel

# simulate some DDM data:
from rlssm.random import simulate_ddm
data_ddm = simulate_ddm(
    n_trials=400,
    gen_drift=.8,
    gen_threshold=1.3,
    gen_ndt=.23)

For the simple, non-hierarchical DDM, it is only necessary to have rt and accuracy columns:

• rt, response times in seconds.

• accuracy, 0 if the incorrect option was chosen, 1 if the correct option was chosen.

data_ddm.head()

		drift	rel_sp	threshold	ndt	rt	accuracy
participant	trial
1	1	0.8	0.5	1.3	0.23	1.134	0.0
	2	0.8	0.5	1.3	0.23	0.433	1.0
	3	0.8	0.5	1.3	0.23	0.702	1.0
	4	0.8	0.5	1.3	0.23	0.695	1.0
	5	0.8	0.5	1.3	0.23	0.351	1.0

# Run 2 chains, with 2000 samples each, 1000 of which warmup, with custom priors:
model_fit_ddm = model_ddm.fit(
    data_ddm,
    drift_priors={'mu':.5, 'sd':1},
    threshold_priors={'mu':0, 'sd':.5},
    ndt_priors={'mu':0, 'sd':.1},
    chains=2,
    iter=2000,
    warmup=1000,
    thin=1)

Fitting the model using the priors:
drift_priors {'mu': 0.5, 'sd': 1}
threshold_priors {'mu': 0, 'sd': 0.5}
ndt_priors {'mu': 0, 'sd': 0.1}

WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.
To run all diagnostics call pystan.check_hmc_diagnostics(fit)

Checks MCMC diagnostics:
n_eff / iter looks reasonable for all parameters
0.0 of 2000 iterations ended with a divergence (0.0%)
0 of 2000 iterations saturated the maximum tree depth of 10 (0.0%)
E-BFMI indicated no pathological behavior

Learning example (hierarchical, real data)

model_rl = rlssm.RLModel_2A(hierarchical_levels = 2)

Using cached StanModel

# import some example data:
data_rl = rlssm.load_example_dataset(hierarchical_levels = 2)

data_rl.head()

	participant	block_label	trial_block	f_cor	f_inc	cor_option	inc_option	times_seen	rt	accuracy
0	1	1	1	43	39	2	1	1	1.244082	0
1	1	1	2	60	50	4	3	1	1.101821	1
2	1	1	3	44	36	4	2	2	1.029923	0
3	1	1	4	55	55	4	3	2	1.368007	0
4	1	1	5	52	49	4	3	3	1.039329	1

Since this learning model is only fit on choices, rt are not required.

Other columns/indexes that should be included are:

• accuracy, 0 if the incorrect option was chosen, 1 if the correct option was chosen.

• trial_block, the number of trial in a learning session. Should be integers starting from 1.

• f_cor, the output from the correct option in the presented pair (the option with higher outcome on average).

• f_inc, the output from the incorrect option in the presented pair (the option with lower outcome on average).

• cor_option, the number identifying the correct option in the presented pair (the option with higher outcome on average).

• inc_option, the number identifying the incorrect option in the presented pair(the option with lower outcome on average).

• block_label, the number identifying the learning session. Should be integers starting from 1. Set to 1 in case there is only one learning session.

If the model is hierarchical, also include:

• participant, the participant number. Should be integers starting from 1.

If increasing_sensitivity is True, also include:

• times_seen, average number of times the presented options have been seen in a learning session.

data_rl.head()

	participant	block_label	trial_block	f_cor	f_inc	cor_option	inc_option	times_seen	rt	accuracy
0	1	1	1	43	39	2	1	1	1.244082	0
1	1	1	2	60	50	4	3	1	1.101821	1
2	1	1	3	44	36	4	2	2	1.029923	0
3	1	1	4	55	55	4	3	2	1.368007	0
4	1	1	5	52	49	4	3	3	1.039329	1

# Run 2 chains, with 3000 samples each, 1000 of which warmup, with thinning and custom priors:
model_fit_rl = model_rl.fit(
    data_rl,
    K=4,
    initial_value_learning=27.5,
    alpha_priors={'mu_mu':-.3, 'sd_mu':.1, 'mu_sd':0, 'sd_sd':.1},
    sensitivity_priors={'mu_mu':-.1, 'sd_mu':.1, 'mu_sd':0, 'sd_sd':.1},
    chains=2,
    iter=3000,
    warmup=1000,
    print_diagnostics=False, # (not suggested, see below)
    thin=2)

Fitting the model using the priors:
alpha_priors {'mu_mu': -0.3, 'sd_mu': 0.1, 'mu_sd': 0, 'sd_sd': 0.1}
sensitivity_priors {'mu_mu': -0.1, 'sd_mu': 0.1, 'mu_sd': 0, 'sd_sd': 0.1}

WARNING:pystan:Maximum (flat) parameter count (1000) exceeded: skipping diagnostic tests for n_eff and Rhat.
To run all diagnostics call pystan.check_hmc_diagnostics(fit)

Diagnostics

As you can see, the MCMC diagnostics are already printed by default (if you do not want this, you can set print_diagnostics to False). I refer to https://mc-stan.org/users/documentation/case-studies/divergences_and_bias.html for an excellent explanation of what these diagnostics actually mean and how to assess them.

On top of these, you can also check the convergence of the chains and the WAIC:

model_fit_ddm.rhat

	rhat	variable
0	1.000630	drift
1	0.999847	threshold
2	0.999237	ndt

model_fit_rl.rhat.describe()

	rhat
count	58.000000
mean	1.000025
std	0.000724
min	0.999071
25%	0.999492
50%	0.999797
75%	1.000427
max	1.001968

model_fit_ddm.waic

{'lppd': -201.4814621446239,
 'p_waic': 3.304311118929816,
 'waic': 409.5715465271074,
 'waic_se': 45.969396831225254}

model_fit_rl.waic

{'lppd': -2632.8620973496854,
 'p_waic': 53.04664551553182,
 'waic': 5371.817485730435,
 'waic_se': 94.06444827215812}

If you want to also see the point-wise WAIC, you can set pointwise_waic to True.

Save the results

By default, the model fit results are saved in the same folder, using the model_label as filename. you can specify a different location using the filename argument.

model_fit_ddm.to_pickle()

Saving file as: /Users/laurafontanesi/git/rlssm/docs/notebooks/DDM.pkl

model_fit_rl.to_pickle()

Saving file as: /Users/laurafontanesi/git/rlssm/docs/notebooks/hierRL_2A.pkl

Re-load previously saved results

model_fit_ddm = rlssm.load_model_results('/Users/laurafontanesi/git/rlssm/docs/notebooks/DDM.pkl')

model_fit_rl = rlssm.load_model_results('/Users/laurafontanesi/git/rlssm/docs/notebooks/hierRL_2A.pkl')

The data the model was fit on are stored in data_info:

model_fit_rl.data_info['data']

	index	participant	block_label	trial_block	f_cor	f_inc	cor_option	inc_option	times_seen	rt	accuracy
0	0	1	1	1	43	39	2	1	1	1.244082	0
1	1	1	1	2	60	50	4	3	1	1.101821	1
2	2	1	1	3	44	36	4	2	2	1.029923	0
3	3	1	1	4	55	55	4	3	2	1.368007	0
4	4	1	1	5	52	49	4	3	3	1.039329	1
...	...	...	...	...	...	...	...	...	...	...	...
6459	6459	27	3	76	37	36	2	1	39	1.875327	1
6460	6460	27	3	77	58	41	4	2	39	1.696957	1
6461	6461	27	3	78	64	49	4	3	38	2.059956	1
6462	6462	27	3	79	44	37	3	1	39	1.623731	1
6463	6463	27	3	80	51	51	4	3	40	1.115363	1

6464 rows × 11 columns

The priors are stored in priors:

model_fit_ddm.priors

{'drift_priors': {'mu': 0.5, 'sd': 1},
 'threshold_priors': {'mu': 0, 'sd': 0.5},
 'ndt_priors': {'mu': 0, 'sd': 0.1}}

model_fit_rl.priors

{'alpha_priors': {'mu_mu': -0.3, 'sd_mu': 0.1, 'mu_sd': 0, 'sd_sd': 0.1},
 'sensitivity_priors': {'mu_mu': -0.1, 'sd_mu': 0.1, 'mu_sd': 0, 'sd_sd': 0.1}}

And different parameter information are stored in parameter_info:

model_fit_ddm.parameters_info

{'hierarchical_levels': 1,
 'n_parameters_individual': 3,
 'n_parameters_trial': 0,
 'n_posterior_samples': 2000,
 'parameters_names': ['drift', 'threshold', 'ndt'],
 'parameters_names_transf': ['transf_drift', 'transf_threshold', 'transf_ndt'],
 'parameters_names_all': ['drift', 'threshold', 'ndt']}

model_fit_rl.waic

{'lppd': -2632.8620973496854,
 'p_waic': 53.04664551553182,
 'waic': 5371.817485730435,
 'waic_se': 94.06444827215812}

And, of course, you can inspect the model’s posteriors, see How to inspect a model for more details on this.