The Algorithm

The computation is described formally in Figure 1. It proceeds in three phases. In Phase 1, temporal trend is removed from each voxel's time series. Then regression factors and correlation coefficients with respect to the ideal time series are computed for each voxel.

Phase 2 iteratively constructs the empirical distribution that will be used to generate probability values. At each step of this phase of the algorithm, the temporal sequence of the acquired images is randomised. Using this randomised sequence, correlations with the ideal time series are computed as in Phase 1. The correlation whose magnitude is the greatest in the entire volume is saved, along with the coordinates of the voxel that produced it. This process is repeated ten thousand times, each time reshuffling the temporal sequence, computing correlations over the entire image, and saving the largest correlation. The resulting set of saved correlations is an empirically derived distribution for the maximum correlation value over the entire image that could be expected to arise under the null hypothesis.

In Phase 3, we repeatedly extract the maximal correlation from the set of actual correlations computed in Phase 1. The coordinates of the voxel that produced this correlation are also retrieved. The correlation is ranked within the empirical distribution computed in Phase 2, and this rank projected onto the interval [0,1] is the adjusted probability level for the voxel in question -- namely, the probability that the maximum correlation produced within an image of unactivated tissue will be less than or equal to this voxel's correlation. If this adjusted probability places the correlation within one of the tails of the distribution, then the voxel is defined as activated.

This activation introduces a problem, though, since we've been tacitly assuming that the empirical distribution contains correlations generated from effectively random data, and the data from an activated voxel are not random. We could continue to use this contaminated distribution, but it would make our statistics less sensitive in the cases of voxels whose correlations are less than that of the current voxel.

To overcome this limitation, we modify Phase 2 of the algorithm to save not only the voxel with the largest correlation, but the next largest, and perhaps the next next largest, and so on, building up a short list of voxels that can be inserted into the distribution as substitutes for voxels that have been defined as activated. In Phase 3, then, we delete from the empirical distribution any and all correlations produced by the activated voxel, and replace each of these deleted correlations with substitute correlations whose voxels of origin have not yet been deleted, if any such substitutes remain available. The total number of substitutes to be held in reserve at each element of the distribution is denoted as `NUMSUBSTS' in the algorithm below; in our experience, 2 usually suffices as a value for NUMSUBSTS.

This entire process is repeated, extracting the greatest correlation, ranking it within the empirical distribution to derive an adjusted probability, and inserting substitute correlations into the empirical distribution, until the adjusted probability of the most recently extracted correlation is no longer significant.

Figure 1. The optimised permutation test algorithm.

Because temporal randomisation occurs at the level of whole images rather than at that of individual voxels, the re-orderings at each voxel are the same (although the particular values that are being re-ordered differ). Thus if the randomised series of some voxel happens to yield a large correlation with respect to the ideal time series, the voxels with which this particular voxel is correlated will also yield large correlations with respect to the ideal time series. Since only one value is saved for each randomisation, only the largest of these correlations will be included in the empirical distribution. In this way, the algorithm automatically accounts for spatial correlations that would otherwise inflate the tails of the empirical distribution and so decrease the computed significance levels [Locascio & al. 1997].

Hochberg and Tamhane [1987] define two forms in which a multiple-comparisons test can control Type I error. In weak control, only the omnibus test need be valid, i.e., the probability of declaring that some voxel somewhere is active when in fact none are active must not exceed . Strong control requires validity not only of the omnibus conclusion but also of the voxelwise tests, again with the specified . Because strong control applies to each voxel considered individually, only strong control permits localisation of significant activations. Thus it is strong control in which we are interested.

To see that the permutation test as described above has strong control, consider any region $U~\subseteq~V$, where V is the entire image.^* (U may be as small as a single voxel.) For all $i, 0{\leq}i<M$, the i^th element in the distribution of maximal correlations over this region must be less than or equal in magnitude to the i^th element in the distribution of maximal correlations over the entire image V: $\vert R_{U}[i]\vert~\leq~\vert R_{V}[i]\vert$. So the rank of any particular correlation r against the distribution R_U must always be at least as far into the tail as the rank of r against R_V. Therefore, the significance levels calculated for voxels within V as a whole can never exceed those for voxels within U by itself.

* The proof here follows that presented in [Holmes & al. 1996].