3dClusterize p thresholding

Hi there,

I am trying to understand something about 3dClusterize. This is my code:

3dClusterize
-inset file.nii
-mask mask.nii
-ithr 9
-idat 8
-NN 1
-2sided p=0.001
-clust_nvox K
-pref_map tmp_output.nii

It works as I expect except that when I compare the cluster map to the manual version of looking at the same file (setting the threshold to p < 0.001 with minimum k clusters in the GUI), the clusters do not always match. Through a little testing, it looks like the p threshold to match the two is actually p < 0.002 (double) in the GUI. I think there is something I am not understanding about how this works. Would appreciate some help clarifying. Thanks!

Howdy-

Hmm, that shouldn’t be the case. A few things to consider here:

  1. What version of AFNI are you using (“afni -ver”)? I assume it shouldn’t matter, but let’s be sure.

  2. When you run Clusterize in the GUI and click “Rpt” → “SaveMsk”, it actually makes a call to 3dClusterize using the parameters you have specified (olay and thr bricks, NN value, threshold, sidedness, etc.). For example, when I ran something similar to your command in the AFNI_data6/afni/ Bootcamp directory, namely for the “func_slim*” statistical dset, and then you will see in the terminal what it has run, for example:


3dClusterize -nosum -1Dformat -inset /home/ptaylor/AFNI_data6/afni/func_slim+orig.HEAD -idat 5 -ithr 6 -NN 1 -clust_nvox 40 -bisided -3.3133 3.3133 -mask_from_hdr -pref_map Clust_mask

Checking that out should reduce the ambiguity; note that the internal p-to-stat calculation has been done here—you can check this with the p2dsetstat command:


p2dsetstat -inset "func_slim+orig[6]" -bisided -pval 0.001

… which here outputs:


++ Found input file : func_slim+orig[6]
++ OK stat type     : fitt
++ BRICK_STATAUX    : 0 3 1 430
++        params    : 430
++ Final stat val   : 3.3133

(you could put “-quiet” for scripting purposes). So, this all looks consistent.

  1. Note that 2sided and bisided are very similar but importantly different. In my case above, in the Clusterize panel, I left “bisided” as on here, so that positive-data-voxels and negative-data-voxels don’t cluster together. In general, it is hard not to see that “bisided” would be more appropriate than “2sided”, where in the latter opposite-sided effects can join forces to form a cluster (though proximity in that case doesn’t seem to warrant considering them as part of the same cluster). But, if you really do want 2sided and not bisided clustering, then after selecting “Bisided → no” in the initial Cluster parameter panel, the output SaveMsk would show this:

+ 3dClusterize -nosum -1Dformat -inset /home/ptaylor/AFNI_data6/afni/func_slim+orig.HEAD -idat 5 -ithr 6 -NN 1 -clust_nvox 40 -2sided -0.0001 0.0001 -mask_from_hdr -pref_map Clust_mask

Note that the p2dsetstat calculation for 2sided or bisided is the same, so:


p2dsetstat -inset "func_slim+orig[6]" -2sided -pval 0.001

produces the identical-to-above:


++ Found input file : func_slim+orig[6]
++ OK stat type     : fitt
++ BRICK_STATAUX    : 0 3 1 430
++        params    : 430
++ Final stat val   : 3.3133

  1. A more subtle point in this: I was not loading in a mask dset for 3dClusterize here; there isn’t a way to do that on-the-fly in the GUI. However, if you notice in the SaveMsk output in each case above, there is an option being used: “-mask_from_hdr”. From the 3dClusterize help, we see what this means:

-mask_from_hdr :If 3dClustSim put an internal attribute into the
                 input dataset that describes a mask, 3dClusterize will
                 use this mask to eliminate voxels before clustering,
                 if you give this option (this is how the AFNI
                 Clusterize GUI works by default).  If there is no
                 internal mask in the dataset header, then this
                 doesn't do anything.

The logic for doing this: to do clusterizing, you typically run 3dClustSim or something to determine a cluster threshold; using 3dClustSim, you need to imput a mask to define the volume within which to estimate the spatial smoothness of noise; that is an important mask for further calcs, so let’s keep it close by… even in the header, so it doesn’t get lost.
… and to check if there really is a mask attached to the dset, you can ask 3dClusterize to output it:


-out_mask OM   :specify that you wanted the utilized mask dumped out
                 as a single volume dataset OM.  This is probably only
                 really useful if you are using '-mask_from_hdr'.  If
                 not mask option is specified, there will be no output.

or you can check the attributes like this:


3dAttribute AFNI_CLUSTSIM_NN1_1sided func_slim+orig.

… and indeed, func_slim* here does have one. So, if on a command line test I didn’t include the mask, I might run into troubles. I don’t know that this is your case here, though.

OK, so, after all that (probably a longer reply than you wanted!), to check that things are consistent, can you please use Clusterize → Rpt → SaveMsk to dump out the 3dClusterize command from the GUI, and see what is being output? You can verify with that and p2dsetstat that the p-to-stat calc is consistent with your command line call, too.

Please let me know how that goes, and if things seem consistent. For me, running this in the command line:


3dClusterize -inset func_slim+orig. -ithr 6 -idat 5 -NN 1 -bisided p=0.001 -clust_nvox 40 -pref_map OUT_CMDLINE_map_bisided.nii -out_mask OUT_CMDLINE_mask_bisided.nii -mask_from_hdr -overwrite

… matched with the GUI results after thresholding with p=0.001, opening the Clusterize panel and setting NN=1, clustsize=40, and leaving Bisided->yes, and then dumping things out. But note that I had to include the “-mask_from_hdr” explicitly in my command line call.

–pt

Hi Paul,

Thank you for the detailed response. After investigating this a little more, I see that I have conflated two separate issues. With the example I gave above, the two actually are pretty much identical (though I take your point on the wisdom of using bisided vs. 2sided). Embarrassingly, I think the discrepancy I saw was because I checked the “pos” button and thought I was missing some clusters. Oops!

However, the double p value does come up when I use 3dClusterize with a 1-sided test. So here is what I would request, in this case from an F-statistic.

3dClusterize -inset filename.nii
-mask mask.nii -ithr 5
-idat 5 -NN 1 -1sided RIGHT_TAIL p=0.001
-clust_nvox 11 -pref_map test3.nii

Here is the GUI report for what I think should be the equivalent (it isn’t):
3dClusterize -nosum -1Dformat -inset filename.nii
-idat 5 -ithr 5
-NN 1 -clust_nvox K
-2sided -7.3247 7.3247
-pref_map Clust_mask

  • I can confirm 7.3247 is the threshold for p < 0.001 in the GUI in this example

To actually match the clusters generated through my 3dClusterize command, I need to change the p value to 0.002 in the GUI. I guess this narrows it down to my misunderstanding how the p-value is set for a 1-sided versus 2-sided test. Two follow up questions would be:

  1. Should I be using a bisided test even for this F statistic? (I thought 1-sided was appropriate here).
  2. Since I requested a 1 sided test, is this why the p-value is double in the 3dClusterize output? Do I need to adjust the desired p value to half in the command?

Thanks again for your help!

Howdy-

No worries about the “pos” thing—I have done the same maaaany atime.

Re. thresholding with the F-stat: the issue here is that it is only >=0, so in converting a p-to-stat for thresholding, should one use 1sided or 2sided? (There is no issue with “bisided” here, because 2sided and bisided are exactly the same for this uniformly positive dataset, unlike if we had a t-stat or correlation or something.)

Note that there is a difference between the conversion, whether you choose 1sided or 2sided:


$ p2dsetstat -pval 0.001 -inset func_slim+orig."[0]" -1sided -quiet
6.3053
$ p2dsetstat -pval 0.001 -inset func_slim+orig."[0]" -2sided -quiet
7.01992

Note that when the GUI converts from p-to-stat, it is always using 2sidedness (again, from above, for p-to-stat conversions, 2sided and bisided are the exact same; the difference only comes in with clumping, afterwards). Testing this with func_slim*'s F-stat, when I select the following in the GUI, Thr->“Set p=0.001”, the threshold goes to 7.01992. So, indeed, the GUI is hardwired to always use 2sided conversions.
If I then go through Clusterize->Rpt->SaveMsk, I do see from the command line dump that “-bisided” has been used, so the GUI is doing things consistently—it is just always using 2sided (or bisided), even for an F-stat.
If I instead in the GUI I use Thr->“Set p-value”->0.002, then the threshold goes to 6.3053. So, from above, we see that this “trick” of doubling my p-value P to get the equivalent of a 1sided conversion of P is needed in the GUI. That is the way you would get the “equivalent” of your first 3dClusterize command (the explicitly 1sided one) in the GUI.

Hopefully that makes sense and seems consistent. To me, the command line is clearest. One does just have to know that the GUI’s p-to-stat conversions (at present) are always 2sided, and then work around that for visualization.
And indeed, I think a 1sided test on the F-stat makes sense, yes.

Also, note that you can images of all these clustering things automagically, too, with @chauffeur_afni. Some examples are here:
make https://afni.nimh.nih.gov/pub/dist/doc/htmldoc/tutorials/auto_image/auto_%40chauffeur_afni.html#ex-3-threshold-stats-voxelwise-view-effects
… even using the translucent thresholding, which appeals a lot.

–pt

Thanks Paul, this makes a lot of sense. Luckily it means all my code is actually correct, I just have to adjust when I visualize. A great outcome! :slight_smile:

Great, glad that things are already on a good track code-wise, and that this clarifies the GUI stuff. And indeed, there is a lot going on with these parameters to consider.

–pt