FW

13. Correction for multiple comparisons: Monte Carlo simulation

One major consideration when conducting voxelwise statistical tests is the sheer number of tests that are being performed. In such cases, Type I error inflation is a problem. A common approach to correct for this is to use a clustering method that tries to determine how big a contiguous cluster of voxels, each one significant at an uncorrected threshold of P_u, has to be to in order to be significant at a threshold P_c – that is corrected for multiple comparisons. To determine how big such a cluster needs to be, there are 3 common methods:

1.      Guassian Random Field Theory: Uses the mathematical properties of random fields to derive a corrected cluster size.

2.      Permutation tests: Uses the data itself to create a null distribution of cluster sizes, from which the cluster size corresopnding to a desired corrected P can be read off.

3.      Monte Carlo simulation: Creates multiple simulated null datasets, and from them creates a distribution of cluster sizes, from which the cluster size corresponding to a desired corrected P can be read off.

We could go into a lot of detail about the relative merits of each method here, but that would be somewhat irrelevant for AFNI analysis, since the Monte Carlo method is the only one currently implemented in AFNI (other techniques might become available in the near future).
In order to carry out the Monte Carlo simulation, we need some basic information about the data itself, as well as the uncorrected P_u that we plan on using for creating clusters, and a few other parameters that control the clustering operation.
First we need an estimate of how spatially smooth the noise in the dataset is. To estimate this, we use the program 3dFWHM:
3dFWHM –dset err_ts+orig[0] –mask ‘3dcalc( –a glm_out+orig[0] –expr step(a-500))’
Here we have used the residual time series from an individual subject GLM to estimate the spatial smoothness. We have masked the time series with the thresholded baseline estimate so as to exclude voxels from outside the brain. This command will return output like the following:

Gaussian filter widths:
sigmax = 1.72 FWHMx = 4.05
sigmay = 1.56 FWHMy = 3.67
sigmaz = 1.70 FWHMz = 4.00

So for this dataset, our estimate would be that the full-width-at-half-maximum (FWHM) estimate of spatial smoothness is approximately 4mm. We could then do this for all our subjects and average the results. However, if we used spatial blurring on our statistical percent signal change estimates of more that 4mm, we should use the FWHM of the spatial blur that we applied. In the example above, we used a spatial blur of 5mm, so our estimate of spatial smoothness that we use in Monte Carlo simulations should be 5mm.
The next step is to decide upon the uncorrected P_u that we plan on using for creating clusters. This choice is fairly arbitrary. If we choose a fairly liberal voxelwise P_u, say P_u = 0.01, we will be able to detect large clusters of voxels that are activated at a fairly liberal threshold, but we will miss small clusters that are made up of highly significantly activated voxels. If, on the other hand, we use a conservative P_u, say P_u = 0.0001, we will be able to detect small clusters of highly activated voxels, but not larger clusters of less activated voxels. The final choice on a P_u depends on the size of clusters that we are looking for, as well as considerations of what our overall statistical power is likely to be.
Let’s say for the current example that we choose a P_u = 0.01.
The other parameter that we need to determine is what we want our cluster connection radius to be. This specifies how close two voxels need to be to one another in order to be considered contiguous. Usually we would use the voxel size of the datasets being analysed, in this case 2mm.
So now we can run the simulation:
AlphaSim –mask glm_out_perc+tlrc[0] –fwhm 5 –rmm 2 –pthr 0.01 –iter 1000
The –iter 1000 specifies that we want to run 1000 simulations. This is usually enough. The program will take quite a few minutes to run, and then produce output like the following:

Cl Size Frequency Cum Prop p/Voxel Max Freq Alpha
 1 229840 0.328718 0.01012002 0 1.000000
 2 108491 0.483882 0.00947785 0 1.000000
 3 66742 0.579336 0.00887161 0 1.000000
 4 51187 0.652544 0.00831218 0 1.000000
 5 39186 0.708588 0.00774011 0 1.000000

and so on. What you now need to do is read down the right-most column, which gives the corrected P_c for each minimum cluster size (left hand column).

84 4 0.999920 0.00001630 4 0.058000
 85 5 0.999927 0.00001536 5 0.054000
 86 1 0.999929 0.00001417 1 0.049000
 87 5 0.999936 0.00001393 4 0.048000

So from this we can see that according to Monte Carlo simulations, thresholding the statistical images with an uncorrected P_u = 0.01, and clustering with a cluster connection radius of 2mm, the resulting clusters need to be at least 86 voxels in size in order to achieve a corrected P_c < 0.05. We can now use this information to threshold and cluster our group analysis statistical maps, which is explained in the next section.

14. Correction for multiple comparisons: Clustering

We can now use the information from Monte Carlo simulations to threshold and cluster our group analysis statistical maps, so as to get thresholded statistical images corrected for multiple comparisons. This we must do separately for each contrast or F-test that we are interested in. For the preceding example, imagine that we are interested in the F-test for the interaction of group by face. First we need to find out where this information is stored in the anova+tlrc bucket file:
3dinfo anova+tlrc
This produces output such as:

-- At sub-brick #0 'group:Inten' datum type is short: 0 to 32767 [internal]
 [* 0.000644032] 0 to 21.103 [scaled]
 -- At sub-brick #1 'group:F-stat' datum type is short: 0 to 2562 [internal]
 [* 0.01] 0 to 25.62 [scaled]
 statcode = fift; statpar = 1 14
 -- At sub-brick #2 'face:Inten' datum type is short: 0 to 32767 [internal]
 [* 0.000571107] 0 to 18.7135 [scaled]
 -- At sub-brick #3 'face:F-stat' datum type is short: 0 to 3689 [internal]
 statcode = fift; statpar = 2 28
 -- At sub-brick #4 'groupByFace:Inten' datum type is short: 0 to 32767 [internal]
 [* 0.000571007] 0 to 19.7125 [scaled]
 -- At sub-brick #5 'groupByFace:F-stat' datum type is short: 0 to 3489 [internal]
 statcode = fift; statpar = 2 28

This tells us that the sub-brick we want to threshold with is sub-brick 5. We can also keep the intensity sub-brick for the interaction, although this would make more sense for a t-test of a specific contrast, since the intensity sub-brick in that case contains the value of the contrast itself (e.g. percent signal change contrast between two conditions).
Before we run the cluster command, we need to determine the F-statistic that corresponds to a P_u of 0.01, for the degrees of freedom that we have for the interaction F-test. This can be read off from tables, calculated using Excel or an online calculator, or read off from the threshold slider in AFNI if viewing the interaction statistical map. In this case, the correct F-statistic is 5.45.
Now we can use the 3dmerge command to perform the clustering:
3dmerge –session . –prefix group_by_face_clust –1thresh 5.45 -1clust 2 –86 anova+tlrc[4,5]
3dmerge –session . –prefix group_by_face_clust_order –1thresh 5.45 -1clust_order 2 –86 anova+tlrc[4,5]
The first command will identify clusters of voxels that exceed an F-value of 5.45, using a cluster connectivity radius of 2mm, and minimum cluster size of 86 voxels, using sub-brick 5 as the thresholding brick, saving the thresholded, clustered data in a new bucket dataset called group_by_face+tlrc. The second command will do the same thing, except all voxel intensities within a cluster will be replaced by the cluster size index (largest cluster=1, next largest=2, ...). This will be useful for extracting data from specific clusters later on. The reason why we specify the minimum cluster size as a negative number is that denotes cluster size in voxels. If we were to put a positive 86 there, that would denote 86mm³.
Some of the other options in the command 3dmerge might be useful. In particular, you can play around with the –1erode and –1dilate options, which are useful in separating clusters that are joined by a narrow “neck”.
We would typically want to carry out the same basic operation on all the other F-tests and contrasts of interest. When this is done, it is possible to combine all the resulting clustered datasets into one big buckets file using the command 3dbucket.
This basically completes the voxelwise group analysis of data.

15. Extracting cluster means from individual subjects

One remaining step that you might want to take is to extract each subject’s individual percent signal change parameter estimates from the cluster regions identified in the preceding step, for entry into a statistics package or for graphing results. To do this, you can use the 3dmaskave command. Say, for example, that when looking at the face by group clusters identified in the preceding step, we were interested in a cluster in prefrontal cortex, which happens to be the 5^th largest cluster in the group by face interaction dataset. Thus voxels in this cluster would have a value of 5 in the group_by_face_clust_order+tlrc dataset. We can thus use this to extract mean voxel percent signal change for individual subjects as follows:
3dmaskave –mask group_by_face_clust_order+tlrc –mrange 5 5 control001/glm_out_5mm+tlrc[4]
This command would print on the screen the mean percent signal change for this cluster for control subject 1 for the happy face condition. We could repeat this command for the other face conditions, as well as for all the other subjects. The mean percent signal change values could then be used to create graphs with error bars etc.
Alternatively, if want to output the mean percent signal change for each cluster, we could use the following command:
3dROIstats -–mask group_by_face_clust_order+tlrc control001/glm_out_5mm+tlrc[4]
This would cycle through each cluster in group_by_face_clust_order+orig and output the mean for each one for control subject 1, face condition happy.
Because these commands simply print their output to the screen, you’ll probably want to use them with the Linux redirection operators > and >> to save the output to a text file. Obviously these commands are probably best used in a script that cycles through subjects and conditions. For example, the following simple Python script will cycle through 6 subjects and extract the cluster means from 4 designated clusters and two conditions.

---

#!/usr/bin/python

# filename: extract_clusters.py
# simple script to cycle through 6 subjects and extract the cluster means from 4 designated clusters and two conditions
# means are first appended to a temporary text file, with subject number, cluster number, and condition brick
# appended to a separate temporary text file. The two files are then combined using 1dcat to give a properly indexed
# file which could be imported into a stats package or spreadsheet program

import os,sys

# set up study-specific directories and file names
top_dir = '/study/my_study/analysis/'
cluster_mask_file = top_dir+'group_by_face_clust_order+tlrc'
subject_stats_filename = 'glm_out_5mm+tlrc'
output_file = top_dir+'clust_means.txt'


# specify the subjects, clusters and bricks that you want to extract
subjects = ['001','002','004','008','009','010'] # these are the subjects (should be same as subject directory names)
clusters = [2,5,6,8] # these are the cluster numbers that we want to extract
condition_bricks = [4,6] # these are the bricks in the individual subject data files that contain the conditions of interest

# check to see for existence of temporary files
if os.isfile('tmp.txt') or os.isfile('tmp2.txt'):
        print 'You first need to delete or rename the file(s) tmp.txt and tmp2.txt'
        sys.exit()

# check to see for existence of output file
if os.isfile(output_file):
        print 'You first need to delete or rename the output file '+output_file
        sys.exit()


# now loop through subjects, clusters and bricks to get data
for subject in subjects:
        for cluster in clusters:
                for brick in condition_bricks:
                        subject_dir = top_dir+subject+'/'
                        data_file = subject_dir+subject_stats_filename+'['+str(brick)+']'
                        command = '3dmaskave –mask '+cluster_mask_file+' -mrange '+str(cluster)+' '+str(cluster)+' '+data_file+' >> tmp.txt'
                        print command
                        os.system(command)
                        command = 'echo '+subject+' '+str(cluster)+' '+str(brick)+' >> tmp2.txt'
                        print command
                        os.system(command)

# create the output file with appropriate column headings
command  = 'echo subject cluster brick contrast > '+output_file
print command
os.system(command)

# now combine the two temporary file and then delete them
command  = '1dcat tmp2.txt tmp.txt[0] >> '+output_file
print command
os.system(command)

os.remove('tmp.txt')
os.remove('tmp2.txt')

---

Assuming the script is saved as extract_clusters.py, this script can easily be run by typing:

python extract_clusters.py
You can learn more about Python scripting here.