3dClustSim

Hi AFNI team,

I was using 3dclustsim to estimate the probability of false positive (noise-only) clusters. I estimated the smoothness values using 3dFWHMx (-acf). I obtained -acf values of 0.5128 3.8451 11.3839 which I inputted to the 3dclustsim.

The issue I face is that every time I run 3dclustim , I get different cluster size. Sometimes the values had a difference of more than 10. I was wondering if there is any issues with the way I am running 3dclustsim.

I have attached the screenshot of the output (N2 bi-sided) for reference.

Code used:

3dClustSim
-mask $temp_path/MNI152_T1_2mm_brain_GM_02182017.nii.gz
-acf 0.5128 3.8451 11.3839
-nodec
-pthr 0.001
-athr 0.05
-prefix 3dClustSimOutput_censor_MR_6mm_ACF.txt

AFNI version=AFNI_20.2.16

I have also tried running 3dclustsim on a different data set with different smoothness values and yet I find cluster size differing every time I run 3dclustsim .

There is some stochasticity in these calculations, so they will be likely to vary a litttle bit across runs. 3dClustSim calcs are also generally intended to be run in more spherical, whole-brain masklike masks, as well. Having a mask that is composed of set of linear-ish structures with complicated boundaries can affect some of convergence, I would guess. I mention this because from your filenames, I think you are running this in a gray matter only mask, which will be relatively filamentary rather than blobby/spherical-ish.

That being said, I generated a practice GM-only mask on the same dimensions of template you are using, using the GM mask from the MNI152 template resampled onto the 91x109x91, 2mm iso MNI_avg152T1+tlrc brain (see attached image):


3dresample -prefix TMP_RES.nii.gz -master MNI_avg152T1+tlrc. -input MNI152_2009_template_SSW.nii.gz'[4]' -rmode NN

This gave me a GM mask that is approximately the same size as yours:


3dROIstats -nomeanout -quiet -nzvoxels -mask TMP_RES.nii.gz TMP_RES.nii.gz

… shows that this mask is 178964 voxels, and yours was 160,742 (just a 10% diff). NB: this might not be the best mask for your template specifically- I did grab this one from a different MNI152 just for testing purposes, so you should verify that whatever you are using is a good GM map. The MNI_avg152T1+tlrc. template itself is relatively low-resolution and low-contrast, so distinct GM-WM boundaries are difficult to tell, anyways.

Then I ran 3 runs of your 3dClustSim command, using the same ACF parameters (note that I let the whole output tables be produced, rather than just one pvalue and one FWE value):


3dClustSim -mask TMP_RES.nii.gz -acf 0.5128 3.8451 11.3839 -nodec -prefix TMP_RES_CLUST
3dClustSim -mask TMP_RES.nii.gz -acf 0.5128 3.8451 11.3839 -nodec -prefix TMP_RES_CLUST2
3dClustSim -mask TMP_RES.nii.gz -acf 0.5128 3.8451 11.3839 -nodec -prefix TMP_RES_CLUST3

… and the pthr=0.001 with alpha=0.05 are quite similar in all cases: 46, 46, 47. The full output tables are below, just for reference.

So, I am not seeing a huge amount of variability. Would you mind posting an image of your GM mask, perhaps that might help show some relative difference? The GM one used here doesn’t have a lot of sharp structure, it should be noted. It is possible that with a mask that has a complicated structure and boundary, more iterations would be necessary for convergence. The default is 10000 iterations, and you could try increasing it: “-iter 20000” or “-iter 30000”, perhaps. Those would be my main suggestions (checking the GM mask for complicatedness, as a reason why convergence might be more complicated, and then increasing the number of iterations to see if that suggests better converge.)

–pt


# 3dClustSim -mask TMP_RES.nii.gz -acf 0.5128 3.8451 11.3839 -nodec -prefix TMP_RES_CLUST
# bi-sided thresholding
# Grid: 91x109x91 2.00x2.00x2.00 mm^3 (178964 voxels in mask)
#
# CLUSTER SIZE THRESHOLD(pthr,alpha) in Voxels
# -NN 2  | alpha = Prob(Cluster >= given size)
#  pthr  | .10000 .05000 .02000 .01000
# ------ | ------ ------ ------ ------
 0.050000     757    910   1135   1325
 0.020000     294    362    452    539
 0.010000     165    202    256    296
 0.005000     100    121    154    185
 0.002000      56     68     88    103
 0.001000      38     46     61     71
 0.000500      26     32     42     50
 0.000200      16     21     27     33
 0.000100      11     15     20     24


# 3dClustSim -mask TMP_RES.nii.gz -acf 0.5128 3.8451 11.3839 -nodec -prefix TMP_RES_CLUST
# bi-sided thresholding
# Grid: 91x109x91 2.00x2.00x2.00 mm^3 (178964 voxels in mask)
#
# CLUSTER SIZE THRESHOLD(pthr,alpha) in Voxels
# -NN 2  | alpha = Prob(Cluster >= given size)
#  pthr  | .10000 .05000 .02000 .01000
# ------ | ------ ------ ------ ------
 0.050000     757    910   1135   1325
 0.020000     294    362    452    539
 0.010000     165    202    256    296
 0.005000     100    121    154    185
 0.002000      56     68     88    103
 0.001000      38     46     61     71
 0.000500      26     32     42     50
 0.000200      16     21     27     33
 0.000100      11     15     20     24


# 3dClustSim -mask TMP_RES.nii.gz -acf 0.5128 3.8451 11.3839 -nodec -prefix TMP_RES_CLUST3
# bi-sided thresholding
# Grid: 91x109x91 2.00x2.00x2.00 mm^3 (178964 voxels in mask)
#
# CLUSTER SIZE THRESHOLD(pthr,alpha) in Voxels
# -NN 2  | alpha = Prob(Cluster >= given size)
#  pthr  | .10000 .05000 .02000 .01000
# ------ | ------ ------ ------ ------
 0.050000     763    922   1140   1300
 0.020000     291    355    449    513
 0.010000     166    201    248    290
 0.005000     100    123    156    182
 0.002000      56     70     89    104
 0.001000      38     47     61     71
 0.000500      26     33     43     50
 0.000200      16     21     28     34
 0.000100      11     15     20     25

Hi taylor !

Thank you for your prompt response. Yes, I am using a gray matter mask. I have upload the GM mask that we usually use in our lab. This mask was obtained from FSL MNI 152 avg and several changes were done on it.

Link to GM mask: https://indianinstituteofscience-my.sharepoint.com/:u:/g/personal/sahithyans_iisc_ac_in/ETzywpowH-lKsxWaLRZSDBQBo87nuXFVC8maHS0egqoyDQ?e=3eLnDL

Can you please let me know if there is any issue with the GM mask attached ?

I also recently started using the GM mask provided in the SSwarper. So, I resampled it using :

3dresample -dxyz 3.0 3.0 3.0 -prefix MNI152_2009_template_SSW_3mm_brain_GM.nii.gz -input MNI152_2009_template_SSW.nii.gz’[4]’

This resampled GM mask shows that a grid of 64x76x64 3.00x3.00x3.00 mm^3 (53104 voxels in mask)
I also used -iter 30000 and ran thrice using the same ACF parameters ( 0.5128 3.8451 11.3839):

3dClustSim
-mask $temp_path/MNI152_2009_template_SSW_3mm_brain_GM.nii.gz
-acf 0.5128 3.8451 11.3839
-nodec
-iter 30000
-prefix 3dClustSimOutput_censor_MR_ACF.txt

and the pthr=0.001 with alpha=0.05 I get cluster sizes of : 15, 17, 16 . (The full output tables are below, just for reference).

3dClustSim -mask /media/cogemolab/home2/templates/MNI152_2009_template_SSW_3mm_brain_GM.nii.gz -acf 0.5128 3.8451 11.3839 -nodec -iter 30000 -prefix 3dClustSimOutput_censor_MR_ACF.txt

bi-sided thresholding

Grid: 64x76x64 3.00x3.00x3.00 mm^3 (53104 voxels in mask)

CLUSTER SIZE THRESHOLD(pthr,alpha) in Voxels

-NN 2 | alpha = Prob(Cluster >= given size)

pthr | .10000 .05000 .02000 .01000

------ | ------ ------ ------ ------

0.050000 221 267 329 375
0.020000 87 104 132 152
0.010000 49 59 75 87
0.005000 30 37 47 54
0.002000 17 21 27 32
0.001000 12 15 19 22
0.000500 9 11 14 16
0.000200 6 7 9 11
0.000100 4 5 7 8

3dClustSim -mask /media/cogemolab/home2/templates/MNI152_2009_template_SSW_3mm_brain_GM.nii.gz -acf 0.5128 3.8451 11.3839 -nodec -iter 30000 -prefix 3dClustSimOutput_censor_MR_ACF.txt

bi-sided thresholding

Grid: 64x76x64 3.00x3.00x3.00 mm^3 (53104 voxels in mask)

CLUSTER SIZE THRESHOLD(pthr,alpha) in Voxels

-NN 2 | alpha = Prob(Cluster >= given size)

pthr | .10000 .05000 .02000 .01000

------ | ------ ------ ------ ------

0.050000 238 285 353 408
0.020000 96 115 144 167
0.010000 55 67 83 97
0.005000 34 41 52 60
0.002000 20 24 30 36
0.001000 14 17 21 25
0.000500 10 12 15 18
0.000200 7 8 10 12
0.000100 5 6 8 9

3dClustSim -mask /media/cogemolab/home2/templates/MNI152_2009_template_SSW_3mm_brain_GM.nii.gz -acf 0.5128 3.8451 11.3839 -nodec -iter 30000 -prefix 3dClustSimOutput_censor_MR_ACF.txt

bi-sided thresholding

Grid: 64x76x64 3.00x3.00x3.00 mm^3 (53104 voxels in mask)

CLUSTER SIZE THRESHOLD(pthr,alpha) in Voxels

-NN 2 | alpha = Prob(Cluster >= given size)

pthr | .10000 .05000 .02000 .01000

------ | ------ ------ ------ ------

0.050000 228 275 342 395
0.020000 91 110 139 161
0.010000 51 63 79 92
0.005000 32 39 49 57
0.002000 18 23 28 33
0.001000 13 16 20 23
0.000500 9 11 14 17
0.000200 6 8 10 11
0.000100 5 6 7 9

And I also tried running 30000 iteration using the old GM mask (attached file;Grid: 91x109x91 2.00x2.00x2.00 mm^3 (160742 voxels in mask) ), for pthr=0.001 with alpha=0.05 I get cluster sizes of : 49, 46 (The full output tables are below, just for reference)

So, is it an issue with the old GM mask used ?
Also another question, I understand that there will be some stochasticity while calculating the cluster size, but how much variation between runs would be considered a good estimate ?

3dClustSim -mask /media/emocoglab/home2/templates/MNI152_T1_2mm_brain_GM_02182017.nii.gz -acf 0.5128 3.8451 11.3839 -nodec -iter 30000 -prefix 3dClustSimOutput_censor_MR_6mm_ACF.txt

bi-sided thresholding

Grid: 91x109x91 2.00x2.00x2.00 mm^3 (160742 voxels in mask)

CLUSTER SIZE THRESHOLD(pthr,alpha) in Voxels

-NN 2 | alpha = Prob(Cluster >= given size)

pthr | .10000 .05000 .02000 .01000

------ | ------ ------ ------ ------

0.050000 748 897 1099 1279
0.020000 297 362 448 521
0.010000 168 206 261 306
0.005000 103 126 161 191
0.002000 58 72 93 111
0.001000 39 49 65 78
0.000500 27 35 45 54
0.000200 17 22 29 35
0.000100 12 16 21 26

3dClustSim -mask /media/emocoglab/home2/templates/MNI152_T1_2mm_brain_GM_02182017.nii.gz -acf 0.5128 3.8451 11.3839 -nodec -iter 30000 -prefix 3dClustSimOutput_censor_MR_6mm_ACF_run2.txt

bi-sided thresholding

Grid: 91x109x91 2.00x2.00x2.00 mm^3 (160742 voxels in mask)

CLUSTER SIZE THRESHOLD(pthr,alpha) in Voxels

-NN 2 | alpha = Prob(Cluster >= given size)

pthr | .10000 .05000 .02000 .01000

------ | ------ ------ ------ ------

0.050000 716 859 1072 1231
0.020000 284 345 428 490
0.010000 160 193 247 289
0.005000 97 119 154 180
0.002000 54 67 87 102
0.001000 36 46 59 71
0.000500 25 32 42 50
0.000200 16 20 27 33
0.000100 11 14 19 24

Hi-

Thanks for sending the mask. I show it overlaid on the one I had been testing—it is indeed a bit narrower, esp around the ventricles. I don’t know that that is an inherent problem, but those narrow regions are kind of non-ideal for the 3D cluster calcs and smoothing. It might just be part of life when using a GM mask. When I ran with just the default 10,000 iterations, the cluster sizes were 45, 47 and 49 (see below for full output).

I am not sure why you are switching to a different voxel size in that other mask? Note that regridding does introduce blur/smoothing into data, and that may change things slighty in your final analysis.

It looks like the cluster size has a reasonable convergence around 47+/a couple. Note that in all these things, it is good to keep an eye on the main goal: presenting reasonable results. If results appear/disappear because of having 47 vs 48 voxel thresholding, that is a sign of deeper instability in analyses. Reducing the sensitivity to these kinds of thresholds is discussed in this draft, with transparent thresholding:
https://www.biorxiv.org/content/10.1101/2022.10.26.513929v2
and here with ETAC:
https://pubmed.ncbi.nlm.nih.gov/31115252/

–pt


# 3dClustSim -mask MNI152_T1_2mm_brain_GM_02182017.nii.gz -acf 0.5128 3.8451 11.3839 -nodec -prefix TMP_CSIM_TEST2a
# bi-sided thresholding
# Grid: 91x109x91 2.00x2.00x2.00 mm^3 (160742 voxels in mask)
#
# CLUSTER SIZE THRESHOLD(pthr,alpha) in Voxels
# -NN 2  | alpha = Prob(Cluster >= given size)
#  pthr  | .10000 .05000 .02000 .01000
# ------ | ------ ------ ------ ------
 0.050000     694    828   1024   1203
 0.020000     268    332    415    477
 0.010000     155    190    242    277
 0.005000      93    116    147    173
 0.002000      54     66     84    100
 0.001000      36     45     59     70
 0.000500      25     32     42     52
 0.000200      15     20     27     34
 0.000100      11     14     19     24


# 3dClustSim -mask MNI152_T1_2mm_brain_GM_02182017.nii.gz -acf 0.5128 3.8451 11.3839 -nodec -prefix TMP_CSIM_TEST2b
# bi-sided thresholding
# Grid: 91x109x91 2.00x2.00x2.00 mm^3 (160742 voxels in mask)
#
# CLUSTER SIZE THRESHOLD(pthr,alpha) in Voxels
# -NN 2  | alpha = Prob(Cluster >= given size)
#  pthr  | .10000 .05000 .02000 .01000
# ------ | ------ ------ ------ ------
 0.050000     728    881   1088   1262
 0.020000     287    351    441    505
 0.010000     164    202    251    287
 0.005000     101    124    157    183
 0.002000      58     72     91    108
 0.001000      39     49     62     73
 0.000500      27     34     44     52
 0.000200      17     22     28     34
 0.000100      12     16     21     26


# 3dClustSim -mask MNI152_T1_2mm_brain_GM_02182017.nii.gz -acf 0.5128 3.8451 11.3839 -nodec -prefix TMP_CSIM_TEST2c
# bi-sided thresholding
# Grid: 91x109x91 2.00x2.00x2.00 mm^3 (160742 voxels in mask)
#
# CLUSTER SIZE THRESHOLD(pthr,alpha) in Voxels
# -NN 2  | alpha = Prob(Cluster >= given size)
#  pthr  | .10000 .05000 .02000 .01000
# ------ | ------ ------ ------ ------
 0.050000     698    851   1058   1223
 0.020000     277    335    434    505
 0.010000     156    194    247    292
 0.005000      96    120    153    177
 0.002000      55     68     88    105
 0.001000      37     47     61     72
 0.000500      26     33     43     51
 0.000200      16     21     28     34
 0.000100      11     15     20     25

Thanks for taking a look at our GM mask.

I was using the 3mm (SSW mask) grid for a different study, so I quickly ran 3dclustsim on that just to rule out the if the variability was caused solely due to our GM mask.

But looks like our old GM mask is more vulnerable to fluctuations compared to the GM mask obtained from SSW.

and I completely agree if results appear/disappear due to changing the cluster size by +/- few voxels, there might be a bigger issue at stake.

Thanks for taking time and providing a detailed explanation. And also thanks for pointing towards the ETAC paper.

I think the main aspect is: while it seems reasonable to restrict the voxel analysis to GM (fewer voxels, about 1/3 of the brain volume → less adjustment for multiple comparisons), the practical reality of the geometry of GM is that the analyzed region becomes less of a simple volumetric shape. Both these kinds of cluster simulations and random field theory-based corrections depend on the size and shape of the masked region; they generally probably all work better in “simple”, quasi-spherical (= masked whole brain) masks, rather than in masks that have long filamentary and/or planar structures. The reason is that smoothness is estimated and assumed to be spherically symmetric in a 3D volume, and those structures don’t really follow that assumption.

So, it may be that more iterations are needed for convergence. Here, it seemed to converge relatively well. But any “edge case” results are always, well, borderline, and it is good to check for stability of results reporting in any case. All of these adjustment methods are kind of a separate fix for the issue of implementing a massively univariate analysis over a large number of elements, and they should probably be viewed as somewhat approximate. If your threshold were 47 voxels, is a cluster with 46 voxels really worth ignoring totally and considering the brain “inactive” there while one with 47 is a major result to discuss? The “Highlight, Don’t Hide” article above goes into this aspect more, but translucent thresholding would really be recommended, and then many of these instabilities, differences and information loss reduce greatly.

–pt