« StockwellDs | Software | 3dpc »

randpermute.py

randpermute.py is a Python program that does non-parametric random-permutation analyses of SAM volume data using AFNI. randpermute.py can handle F, t, and z statistical volumes produced by 3dANOVA*, 3dttest, and 3dWilcoxon, for example.

In a nutshell, random permutation analysis performs a test over and over with shuffled data to build up a surrogate distribution of statistical values for each voxel. These are used to calculate corrected p values that more accurately reflect the per-voxel statistics, and are suitable for publication. Note that you need a large enough N to reach the desired level of significance. See the Theory section below for a more detailed description.

Download and Installation

Version 1.6 20070525 Download randpermute.py, make it executable with chmod +x randpermute.py, and put it somewhere in your $PATH.

If you have a version of Python older than 2.3, you'll also need to install the tempfile.py script. You can either install it somewhere along the standard Python system path (run python and type "import sys; sys.path" to find out what that is), or install it anywhere and set the environment variable $PYTHONPATH to include that directory.

An older randpermute program is also still available as source (randpermute.tgz). It is about twice as fast as randpermute.py for the test it does, because it takes advantage of a symmetry.

Usage -— randpermute.py

randpermute.py [-n npermute] [-v] [-d] template

randpermute.py runs a series of AFNI commands to generate the corrected statistic. The AFNI commands to run are built using a template. The template file is structured as a series of sections, each preceded by a section label. The section labels, and the contents of each section, are defined below. The -d option is used for debugging, and does nothing except print the original, unpermuted AFNI command, followed by npermute permuted commands.

cmd

The AFNI command to run, containing the special variable $permute in place of the -dset arguments.

output

The name of the output brik from the AFNI command, as name+view. Do not use -bucket datasets. Subbrik 0 of the output brik must contain the 'intensity' image, and subbrik 1 must contain the statistical threshold.

permute

The list of -dset arguments, where each line contains the entire set of elements (condition labels) to be permuted, with $n positional parameters specifying the values to be permuted. See below for an example.

values

The list of values to be substituted into the $n variables of the permute section.

Comment lines begin with '#', and long permute lines may be continued by ending them with backslashes. Special characters, such as subbrik selectors, must be properly quoted for the shell.

Example:

cmd:
3dANOVA -levels 3 $permute -ftr out

output:
out+tlrc

permute:
-dset $1 S1condA+tlrc -dset $2 S1condB+tlrc -dset $3 S1condC+tlrc
-dset $1 S2condA+tlrc -dset $2 S2condB+tlrc -dset $3 S2condC+tlrc
...
-dset $1 S9condA+tlrc -dset $2 S9condB+tlrc -dset $3 S9condC+tlrc

values:
1
2
3

In this example, the original, uncorrected 3dANOVA is run using conditions A, B, and C for the three levels 1, 2, and 3 in the order given, to generate ftr+tlrc, the F-test for the treatment effect. This brik will be renamed uncorr+tlrc. Then, the dataset level arguments are permuted npermute times to generate a surrogate distribution. npermute defaults to the smaller of 1000 or nval ^np, where nval is the number of values and np is the number of lines in the permute section. Larger npermute values allow smaller p values, but take longer to run. The final output brik will be called corr+tlrc, and contains the same functional image as the uncorr brik, but with the threshold subbrik set according to the tail probability of the surrogate distribution.

Here's another example that shows how to do a 3dttest in a manner equivalent to what the old randpermute program does:

cmd:
3dttest -prefix out -base1 0 -set2 $permute 2>&1

output:
out+tlrc

permute:
"3dcalc( -a S1AvsB+tlrc -expr $1 )"
"3dcalc( -a S2AvsB+tlrc -expr $1 )"
...
"3dcalc( -a S9AvsB+tlrc -expr $1 )"

values:
a
-a

Note the quoting of the permute lines, to protect the special characters. The "2>&1" phrase in the 3dttest command means "stderr is copied from stdout," or in other words, all error output will be sent to the normal output stream; this makes the execution quieter because stdout is normally turned off (use -v to turn it back on), and 3dttest sends a lot of output to stderr.

Also note that value list entries can be any strings, and that there are two of them here, but only $1 is used in the permute section. Here, the choice of -a versus a flips active versus control in the SAM briks, which are assumed to be D3 images. The original, unpermuted order just uses a for all the datasets.

Here's another example showing a more complicated ANOVA:

cmd:
3dANOVA3 -type 4 -alevels 2 -blevels 2 -clevels 8 -fa out $permute

output:
out+tlrc

permute:
-dset $1 1 1 s1-A1B1+tlrc -dset $1 2 1 s1-A1B2+tlrc \
-dset $2 1 1 s1-A2B1+tlrc -dset $2 2 1 s1-A2B2+tlrc
-dset $1 1 2 s2-A1B1+tlrc -dset $1 2 2 s2-A1B2+tlrc \
-dset $2 1 2 s2-A2B1+tlrc -dset $2 2 2 s2-A2B2+tlrc
...
-dset $1 1 8 s8-A1B1+tlrc -dset $1 2 8 s8-A1B2+tlrc \
-dset $2 1 8 s8-A2B1+tlrc -dset $2 2 8 s8-A2B2+tlrc

values:
1
2

In this case, it's a 2x2 design with subjects as a random effect. An important thing to notice here is that you can only permute the labels of one level at a time (here, A). The output BRIK specifies the A level treatment effect, and that is what is being varied to build the surrogate distribution. Notice also that you can break long lines with backslashes.

Usage -— randpermute

randpermute takes a set of SAM volumes (typically, one for each subject) and outputs a bash script that does multiple 3dttest commands.

randpermute subj1+tlrc subj2+tlrc … > moo
sh moo

The result will be an uncorr+tlrc BRIK that is the output of

3dttest -prefix uncorr+tlrc -base1 0 -set2 subj1+tlrc subj2+tlrc …

and a corr+tlrc BRIK that contains the corrected p values (along with some other junk files — you probably will want to run randpermute in a separate subdirectory).

Theory

Since the power estimates in each SAM volume are calculated from a single covariance matrix, the voxels are not all independent. This creates a problem for hypothesis testing, where we want to know if some of the voxels are different from the others. To deal with this, we need to have independent statistics for each voxel.

Random permutation analysis builds a distribution for each voxel by randomly flipping the condition labels (such as active vs. control) on each subject's SAM volume. If there are nval ways to do this for each of n subjects, there are nval ⁿ possible permutations, which are selected at random without replacement.

The corrected p value is based on this random sampling. For each voxel, the number of times that the test statistic (t, z, or F) is at least as large as the original one results in a tail probability when normalized by the total number of permuations. Thus, the final p value estimates the probability that the original value occurred by chance. The calculation is done voxel by voxel, so the p values are independent. The calculation also takes into account that the t- and z-tests are two-tailed and the F-test is one-tailed.