Formatting data for radd¶

Required columns¶

• idx: Subject ID number
• ttype: Trial-Type ('go' if no stop-signal, 'stop' if stop-signal trial)
• response: Trial Response (1 if response recorded, 0 if no response)
• acc: Choice Accuracy (1 if correct, 0 if error)
• rt: Choice Response-Time (in seconds, can be any value no-response trials)
• ssd: Stop-Signal Delay (in milliseconds, 1000 on go trials)

Optional columns¶

• input dataframe can contain columns for experimental conditions of interest (choose any name)
• in the dataframe above, the Cond column contains 'bsl' and 'pnl'
  • in the 'bsl' or "Baseline" condition, errors on go and stop trials are equally penalized
  • in the 'pnl' or "Caution" condition, penalties are doubled for stop trial errors (e.g., response=1)
• See below for fitting models with conditional parameter dependencies
  • e.g., drift-rate depends on levels of 'Cond'

Building a model¶

Animation of Dependent Process Model¶

Header of observed dataframe (model.observedDF)¶

• idx: subject ID
• Cond: Baseline(bsl)/Caution(pnl) (could be any experimental condition of interest)
• Acc: Accuracy on "go" trials
• sacc: Mean accuracy on "stop" trials (mean condition SSD used during simulations)
• c10 - c90: 10th - 90th RT quantiles for correct responses
• e10 - e90: 10th - 90th RT quantiles for error responses

Bounded Global & Local Optimization¶

Stochastic Global Optimization (Basinhopping w/ bounds)¶

set_basinparams() method gives control over low-level parameters used for global opt

• xtol = ftol = tol: error tolerance of global optimization (default=1e-20)

• stepsize (default=.05): set basinhopping initial step-size

  • see HopStep class in radd.fit for more details
  • see get_stepsize_scalars() in radd.theta for parameter-specific step-sizes

• nsamples (default=5000): number of parameter sets to sample

  • optimization is performed on best ninits (see below)

• ninits (default=5): number of initial parameter sets to perform global optimization on

  • evaluate all nsamples parameter sets, select ninits with the lowest error
  • perform global optimization using each parameter set (out of the ninits chosen)
  • pass the optimized parameters with the lowest global minimum to local optimizer

• nsuccess (default=40): criterion number of successful steps without finding new global minimum to exit basinhopping

• interval (default=10): number of steps before adaptively updating the stepsize

• T (default=1.0): set the basinhopping "temperature"

  • higher T will result in accepted steps with larger changes in function value (larger changes in model error)

Local Optimization (Nelder-Mead Simplex w/ bounds)¶

set_fitparams() method gives control over low-level parameters used for local opt. Local optimization polishes parameter estimates passed from global optimization step.

• xtol = ftol = tol (default=1e-30): error tolerance of optimization

• maxfev (default=2000): max number of func evaluations to perform

• ntrials (default=20000): num. simulated trials per condition

• quantiles (default=np.array([.10, .20, ... .90]): quantiles of RT distribution to fit

Setting parameters example¶

Set basinparams¶

model.set_basinparams(tol=1e-25, stepsize=.05, nsamples=2000, ninits=3)

Set fitparams¶

quantiles=np.arange(.1, 1.,.05)
model.set_fitparams(tol=1e-35, ntrials=30000, maxfev=2500, quantiles=quantiles)

Steps in Fitting Routine¶

Step 1. Flat Fits¶

• All models are initially fit by optimizing the full set of parameters to the "flattened" data (flat meaning the average data collapsing across all conditions of interest). For instance, at this stage fitting the dependent process model involves finding the parameter values for each included parameter that minimizes the cost-function cost function:

  $$Cost = \sum [\omega * (\hat{Y} - Y)]^2$$

• $Y$ is an array of observed data (e.g., accuracy, RT quantiles, etc.)

• $\hat{Y}$ is an equal length array of corresponding model-predicted values given parameters $\theta$

• The error between the predicted and the observed data ($\hat{Y} - Y$) is weighted by an equal length array of scalars $\omega$ proportional to the inverse of the variance in each value of $Y$.

• The array of weighted differences is summed, yielding a single cost value equal to the summed squared-error ($SSE$).

Step 1a:¶

• Global optimzation on flat data (average values collapsing across experimental conditions)

Step 1b:¶

• Local optimzation using parameters passed from global optimizer as starting values

Step 2. Conditional Fits¶

• Conditional models can be fit in which all parameters from Step 1 are held constant except for one or more designated conditional parameters which is free to vary across levels of an experimental condition of interest. Global and local optimization are performed by minimizing the cost-function:

• where $\sum[\omega_i*(\hat{Y_i} - Y_i)]^2$ gives the $Cost$ for level $i$ of condition $C$

• the total $Cost$ is equal to the $SSE$ across all $N_c$ levels of that condition

• Specifying parameter dependencies is done using depends_on --> {parameter_id : condition_name}.

• For instance, in Dunovan et al., (2015) subjects performed two versions of a stop-signal task

  • Baseline ("bsl") condition: errors on go and stop trials are equally penalized
  • Caution ("pnl") condition: penalties 2x higher for stop trial errors (e.g., response=1)

• Different hypotheses about the mechanism underlying behavioral adaptation to assymetric penalties

• To test the hypothesis that a change in decision threshold (a) provides a better account of the data than a change in the drift-rate (v), where 'Cond' is the name of a column header in your data file

  | How-To |Code | |:---:|:---:| | initiate boundary model |model_a = build.Model(data, depends_on={'a':'Cond'}) | | initiate drift model | model_v = build.Model(data, depends_on={'v':'Cond'}) | | fit boundary model | model_a.optimize() | | fit drift model | model_v.optimize() | | Is boundary model better? | model_a.finfo['AIC'] < model_v.finfo['AIC'] |

Step 2a¶

• Global optimzation of conditional parameters

Step 2b:¶

• Local optimzation of conditional parameters passed from global optimizer

How to access...¶

Typically...¶

• ...you'll have specific hypotheses (H) about which parameters will depend on various task conditions so there's no reason to run step 1 and 2 separately.
• In this case, simply initialize the model with the depends_on dictionary corresponding to your H...

  model = build.Model(kind='xdpm', depends_on={'v': 'Cond'})
  

...and hit the optimize button to fit the full model (Step 1 and 2) py model.optimize(progress=True)

Nested optimization of alternative models¶

• Nested optimization allows alternative models to be optimized using a single initial parameter set.

  m = build.Model(kind='xdpm', data=data, depends_on={'v': 'Cond'})
  nested_models = ['xb', 'a', 'tr']
  m.optimize_nested_models(models=nested_models)
  

• After fitting the model with depends_on={'v': 'Cond'}, 'v' is replaced by the next parameter in the nested_models list ('xb' in this case) and optimized using the same init params as the original model

• As a result, model selection is less likely to be biased by the initial parameter set

• Also, because Step 1 takes significantly longer than Step 2, nested optimization of multiple models is significantly faster than optimizing each model individually

• Run the cell below and go take a leisurely coffee break

Model Comparison¶

Examine Best-Fit Model¶

Troubleshooting Ugly Fits¶

Fit to individual subjects¶

model = build.Model(data=data, fit_on='subjects')

Other "kinds" of models...¶

• Currently only Dependent Process Model (kind='dpm') and Independent Race Model (kind='irace')
• Tell model to include a Dynamic Bias Signal ('xb') by adding an 'x' to the front of model kind
• To implement the Dependent Process Model with dynamic bias:

  model = build.Model(data=data, kind='xdpm')
  

• without:

  model = build.Model(data=data, kind='dpm')
  

• To implement the Independent Race Model with dynamic bias:

  model = build.Model(data=data, kind='xirace')
  

• without:

  model = build.Model(data=data, kind='irace')
  

Other dependencies...¶

• Maybe subjects change their boundary height or go onset time across conditions
• Model with dynamic gain free across conditions:

model = build.Model(data=data, depends_on={'xb':'Cond'})

Optimization parameters...¶

• model.set_basinparams(nsuccess=50, tol=1e-30, ninits=10, nsamples=10000)
• model.set_fitparams(maxfev=5000, tol=1e-35)
• Check out the wts vectors for extreme vals
  • Try re-running the fits with an unweighted model (all wts = 1)
    • m = build.Model(data=data, ... weighted=False)
• Error RTs can be particularly troublesome, sometimes un-shootably so...