Genetic Algorithm

Read me!

By Tor Wager, 7/2/02

Wager, T. D. & Nichols, T. E. (2003) Optimization of Experimental Design in fMRI: A General Framework Using a Genetic Algorithm. Neuroimage, 18, 293-309.

The idea behind the genetic algorithm is that it is an efficient way to search through a large and complicated space. In this case, the complicated space is the space of possible fMRI designs. In order to use a GA, you have to 1) be able to parameterize a design as a list of numbers, and 2) evaluate the fitness of a design using some kind of metric testable in computer simulation. This set of Matlab routines will help you to run a GA to optimize experimental designs (hopefully) a little easier than if you didn’t have them. In the best case, it’s “plug and play” – just create (or modify an example script) a script with the experiment parameters you want and run it.

Text script vs. GUI

This set of functions is called from a single script, which is unique to each experiment. I chose to implement it this way instead of creating a graphical user interface because the text version has several advantages. First, in the text version, you can see and save all the input parameters, and you can easily save each script you run for each experiment as a record. Second, text is easier to use (since there are a number of input parameters that must be specified), and easier to transfer across platforms.

Resources

The principal resources you can use to help you are

1)  This readme file

2)  The paper (currently under review) describing the GA

3)  The code

When you should use the GA (and when you shouldn’t)

Using the GA is good if you have to optimize a randomized event-related design, and you want to detect one or more contrasts across different event types. (A contrast is a linear combination of signal estimates for each event type; for example, 2 different “visual” events versus two different “auditory” events.) Or, it could be good if you want to recover hemodynamic response estimates (HRF shapes) for each event type. Or, it could be good if you want to counterbalance the order of event types up to any number of time steps back (i.e., each trial follows each other one equally often).

If you really want to maximize detection power for a single contrast, consider using a blocked design. That will give you the best detection power. If you can’t use a blocked design because you want to isolate particular events (and avoid confounds), the task can’t be done in a blocked fashion, or for any other reason, the GA may be able to help. Some block designs are better than others, especially with more than two conditions, but the problem can be solved analytically—the strong suit of the GA is in cases where the solution is too hard to predict a priori. ER designs with multiple contrasts generally fall under this category.

If you optimize for contrast detection, what the GA will try to do with your event-related design is to create mini-blocks, or clumps of similar event types. This may make the design somewhat predictable to subjects—after all, it’s not truly random once you start selecting specific trial orders. You can counteract this by specifying maximum numbers of the same trial types that can occur in a row (hard constraint) or optimize for HRF shape estimation and/or counterbalancing in addition to contrast detection.

The original GA was designed to work with single events as models, not sustained epochs. If you have an epoch-related design, there is an alternative GA function that has worked for me in the past, but hasn’t been extensively tested. The normal event-related GA function (that does most of the work) is called optimizeGA.m. The epoch-design GA function is called optimizeGA_epochs.m. Using this second function might be frustrating at this point.

What to do to get started

Here are some simple steps to get you going.

1)  Download the zipped tar file from the website. Unzip it using:

tar zxvf OptimizeDesign10.tar.gz

2)  The resulting folders should be included in your matlab path. Some ways to do this are with the addpath command (try help addpath), pathtool, or specifying the GA folders in the matlab path in your .cshrc file.

3)  Choose one of the example scripts stored in the example_scripts subfolder (ga_example_script.m) and open it in a text editor. The script has commenting with each line you should edit for your experiment that explains some of the options. Type the script name to run it. Here are some of the example scripts:

ga_example_script.m A basic 2-condition design with probe periods, nonlinearity modeling, high-pass filtering, and autocorrelation modeling.

ga_example_script_basic.m Stripped-down 2-condition difference detection

ga_example_script_factorial.m

A 2-factor design optimizing main effects and the interaction

ga_text_script_epochs.m An uncommented example script for an epoch design with 3 epochs per trial (long trials). Use of 3 epochs/trial is mandatory with this version.

ga_text_script_hox.m An uncommented example script for a GA that uses ISI as a variable parameter, and attempts to optimize ISI as well as trial order.

4)  If your script terminates successfully, a variable called M (for “Model”) will be created in the workspace. M is a structure that contains information about your experiment. Refer to the fields of M with the . operator; e.g., typing M.stimlist will show the optimized stimulus presentation sequence. You can use this to program your experiment. Try typing:

modelDiagnostics2(M,'plot')

M.ga is a nested structure contains all the input parameters to the GA, and some measures of performance across generations. More explanation on the fields can be obtained by typing

help modelDiagnostics2.

modelDiagnostics2.m is the function that creates the M structure given a stimulus list and a set of filtering and other experiment parameters.

a-optimality vs. d-optimality

A-optimality is currently the default option in the GA. It minimizes the standard error of each contrast or predictor you’re trying to estimate. This is good for testing an a prior hypothesis – e.g., maximizing power to look for an A – B effect. However, it doesn’t take into account the correlations between estimated coefficients, and so isn’t as good if you want to know WHICH of several trial types are producing a change in a particular brain area. For these types of questions, d-optimality may be more appropriate. D-optimality maximizes the determinant of the Fisher Information Matrix X’X, and imposes a more severe penalty for multicolinearity than a-optimality. I recommend using this option for estimating HRFs.

To use it, make sure GA.dflag = 1 in your script.

Optimizing the spacing of events within multi-part trials

The GA optimizes the order of discrete events or trial types in a sequence. However, many experiments have trials with multiple parts (e.g., encoding, delay, probe, and feedback periods), and we’d like to separate responses to each part. The way to do that is to vary the intervals between each part of the trial. A way to optimize the design by selecting the exact placement of events relative to one another is provided in ga2.m, in the ga2 subfolder. Rather than optimizing the order, GA2 optimizes the timing parameters, including presentation length and jitter. It may be desirable to optimize both the order and timing in sequences with dependencies, but that’s not quite implemented yet. GA2 has its own help files, located in the GA2 subfolder.

“Jitter” in designs

There are two easy ways to introduce random jitter into GA designs. The first is by setting the field GA.freqConditions to be a vector that sums to less than 1 (this is set in the GA script that you run, e.g., ga_example_script.m. freqConditions specifies the proportion of the total number of events for each trial type. For example, if you have two conditions, then [.5 .5] means 50% event type 1 and 50% event type 2. [.4 .6] signifies 40% event type 1 and 60% type 2. [.4 .4] signifies 40% for type 1 and 40% for type 2; the remaining 20% will be blank rest intervals.

The second way is to create an extra condition that will be the “baseline” condition. If you have two event types plus rest, you might set GA.conditions = [1 2 3] to specify three trial types: type 1, type 2, and rest. GA.freqConditions = [.4 .4 .2] signifies 40% for type 1, 40% for type 2, and 20% rest (condition 3), and is equivalent to the example above.

How designs are coded

Designs are coded as a column vector of integers. Each integer represents a different condition, or “event type,” and the elements represent the position in time. Zeros represent rests. So the vector represents the stimulus presentation state at each time interval throughout the run, where the intervals are defined at an arbitrary resolution (stored in GA.ISI). ISI stands for inter-stimulus interval; in practice, that’s a pretty good time resolution to use, as design vectors too fine-grained take a long time to process. Another way of saying this is that you must specify some stimulus-presentation state at some regular interval (ISI). Zeros indicate rest intervals, and they are not modeled in the model matrix or considered in counterbalancing or other fitness measures.

Much of the GA script and the main GA function (optimizeGA.m) are for the purpose of constructing a number of randomized lists of event codes according to your specifications. These vectors, concatenated into a matrix of column vectors stored in the variable listMatrix, forms the starting point for the GA. The rest of optimizeGA.m creates model matrices (the X matrix in a statistical GLM analysis) from each vector, tests them with fitness (or “objective”) functions, and performs crossover on the best half of each generation. This is illustrated graphically above.

GA timeline: Sequence of operation

Here’s a slightly more detailed look at what the GA does at what point during its operation.

Your_script.m creates GA structure variable with all necessary input values

For 1 to nmodels (all these names refer to variables defined in the script), runs the function optimizeGA.m

OptimizeGA.m the main GA function; can also substitute optimizeGA_epochs.m (discretion is advised). This function does the following:

1)  SETUP: set up contrasts, nonlinear saturation, periodic rest intervals, filtering and intrinsic autocorrelation matrices, task switching (trial order effects) and mixed block/ER setup (if specified), create initial population of randomized designs (the “dna”), hard-constraint criteria setup, approximately in that order. Prints output so you can check it.

2)  Loop through generations; for each generation, and each list in each generation runs these functions in order (if necessary given your input specifications):

  1. dna2rna.m: Create intermediate design vector (“rna”) from the original event list (“dna”)
  2. getCounterBal.m: Counterbalancing fitness calculation
  3. calcFreqDev.m: Calculate deviation from input frequencies for each event-type (fitness measure)
  4. getMaxInARow.m: Calculate maximum repeats for the list (hard-constraint measure)
  5. designvector2model.m: Contrast detection model matrix builder.
  6. calcEfficiency.m: Compute contrast detection efficiency.
  7. tor_make_deconv_mtx2.m: Construct model matrix for HRF shape estimation efficiency.
  8. CalcEfficiency.m: Run again for HRF shape.

3)  After each generation

  1. GetOverallFitness.m: calculates overall fitness score given results of all fitness tests.
  2. Save the best of this generation.
  3. Breedorgs.m: Take the best half (or roughly so, depending on value of “alpha” input parameter) of designs (“dna”) and crossover, exchanging upper halves. Returns a new population of design vectors.
  4. Print output for the generation.
  5. Every 10 generations, save intermediate results in listMat_working. If the GA quits partway through, you can use the listMatrix (the population list) variable in this file to start again where the GA left off.

4)  After all generations

  1. Create M structure using the best overall design, with all diagnostic measures.
  2. modelDiagnostics2.m: adds resampled and filtered design matrices to M structure, as well as efficiency, variance inflation, condition number, predictor correlation matrix, average time between repeats, and other design fitness/evaluation measures. Also plots predictors and contrasts graphically, if specified.
  3. M structure is saved in model#_<date>.mat, where # is the #’th iteration of the GA (out of nmodels).

Additional functions not further explained here

A number of parameters coded in the scripts and functions are “add-ons” designed to perform specific functions. If you generally leave these at their default values when getting started, you should be OK. See the simple ga_example_script.m for defaults. Some of these functions (and their defaults) are:

GA.doepochs = 0; Build epochs instead of events from design vector

Use with optimizeGA_epochs.m

GA.numhox = 0; “Hox” genes (by analogy) allow you to optimize things like ISI directly with the GA. This is experimental. It’s worked for me, but it hasn’t been widely tested, and probably only works under specific conditions. The idea is simple; include ISI, or any other input parameter, in the GA with a random value between two limits. Then optimize that parameter value with the GA. In practice, some input parameters interact with how model matrices are constructed, and some are very computationally intensive to implement. The current ones are hard-coded to represent [ISI TR cuelen cuerest stimlen stimrest resplen resprest (in s)] in that order. ISI should work, TR does not work (specify and leave at whatever your TR is), and the rest code for epoch lengths for a long single-trial epoch design with multiple periods (cue, stimulation, response) within the trial.

GA.hoxrange = []; Starting ranges for each Hox parameter to optimize.

GA.alph = 2.1; "selection pressure": higher is more extreme selection.

GA.lowerLimit = []; From an older version, not used; don’t worry about it.

GA.restevery = [];

GA.restlength = []; These two parameters insert rest intervals (0’s) of length restlength (in units of the ISI) every restevery ISIs. For example, restevery = 5, restlength = 3 would insert 3 zeros into each design vector (3 ISIs) every 5th element. This is a convenient way to model instruction and probe periods, which occur periodically and regularly, but are not of interest in the study. These periods do affect the design efficiency, so including them makes the simulation a bit more realistic.

GA.trans2switch = 0; This parameter is set to 1 when you have a history-dependent design. For example, if your predictors of interest are switches between two event types [1 2], and you set trans2switch to 1, conditions 1 and 2 will represent events A and B, where A is preceded by A and B is preceded by B (“no switch”). Conditions 3 and 4 will be created automatically, and they will represent A and B, where A is preceded by B and B by A (“switch”) condition. Contrasts will then have 4 conditions, potentially representing main effects of A vs B, main effects of switching, and interactions. This special feature works only when GA.conditions = [1 2].