Stats (e.g. 3dMVM) on anatomical data (e..g VBM)

Robin · September 14, 2021, 1:36pm

Hi!

If I would have normal fMRI-data and run a factorial 3dMVM analysis I know the steps to get “valid” cluster-corrections using the errts-files to in 3dFWHMx -ACF to generate the acf-paramters followed by simulations in 3dClustSim where you have to pick a “maximum” a per-voxel threshold of 0.002 to be able to claim that the results are false-pos corrected according to AFNI-guidelines.

But this time we have single 3D volume structural maps (from the FSL-VBM pipeline). We now want to do a 3x2 factorial analysis using 3dMVM. The analysis went fine but how do we tackle the cluster-correction question here?

Our approach, since this is not bold-data, is to use the data itself with brain-masks to get smoothness via 3dFWMHx -ShowMeClassicFWHM . We argued that we can use the old -FMHM method since this is not temporal bold-data. I.e. get the x-y-z smoothness estimates and plug them into 3dclustsim -fwhmxyz with the group-test-mask used in 3dMVM.

Would this be correct? If so, is there any “maximum” allowed per voxel p-value to use, like the 0.002 in the fMRI guidelines?

Best,
Robin

dglen · September 15, 2021, 2:30pm

I’m not sure if this is helpful, but these two papers prefer permutation approaches (the second paper offers a cluster based approach too) instead for VBM analysis.

http://dx.doi.org/10.1016/j.neuroimage.2008.01.003
https://pubmed.ncbi.nlm.nih.gov/27785843/

Gang · September 15, 2021, 3:06pm

Robin,

A couple of suggestions. Use the option -resid in 3dMVM to obtain the residuals. Then estimate the spatial smoothness with 3dFWHMx, and proceed with usual ritual. On the other hand, the cluster approach is a Band-Aid solution anyway, so here is my radical opinion:

https://afni.nimh.nih.gov/afni/community/board/read.php?1,166126,166126#msg-166126

pmolfese · September 15, 2021, 3:50pm

I’ll echo Gang’s answer about using the residuals from 3dMVM as likely the best option.

Years ago, 2014 apparently - I did a similar approach to what you proposed of estimating smoothness with 3dFWHMx in this paper: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3926206/

That was of course early in the entire cluster exploration process and would still say is a valid approach, especially if you’re conservative in your cluster sizes and p-value thresholds. I had also evaluated these data using permutation methods and found the results to be rather similar.

Robin · September 16, 2021, 9:57am

Just want to thank you all for your very helpful input!
So very appreciated!

Robin · September 17, 2021, 11:14am

We will try the “randomize” function but if I was to use the approach with 3dMVM -resid:
This gives me a _resid output file with 1800 values per voxel (we have 1800 subjects).

Should I loop through this sub-brik with 3dFWHMx, where the loop index is used to give 3dFWHMx one sub-brick at a time or can I just give it the full 4D-volume?
Should I use any of these option?


-demed      = If the input dataset has more than one sub-brick
                (e.g., has a time axis), then subtract the median
                of each voxel's time series before processing FWHM.
                This will tend to remove intrinsic spatial structure
                and leave behind the noise.
                [Default = don't do this]
  -unif       = If the input dataset has more than one sub-brick,
                then normalize each voxel's time series to have
                the same MAD before processing FWHM.  Implies -demed.
                [Default = don't do this]
  -detrend [q]= Instead of demed (0th order detrending), detrend to

>>>>>
  -geom      }= If the input dataset has more than one sub-brick,
    *OR*     }= compute the final estimate as the geometric mean
>>>>>
  -combine    = combine the final measurements along each axis into
                one result

Should I for structural data still use -ACF or is classical -fwhm sufficient?

I tried this:


3dFWHMx -input resid+tlrc. -ShowMeClassicFWHM

Got this:


++ start ACF calculations out to radius = 8.19 mm
 + ACF done (0.00 CPU s thus far)
 2.73652  2.60893  2.85135     2.73047
 0.773993  1.5801  7.99712    4.18496
++ ACF 1D file [radius ACF mixed_model gaussian_NEWmodel] written to 3dFWHMx.1D
++ 1dplot: AFNI version=AFNI_21.0.12 (Feb 25 2021) [64-bit]
++ Authored by: RWC et al.
 + and 1dplot-ed to file 3dFWHMx.1D.png

Are those values for all 1800 subjects?:
2.73652 2.60893 2.85135 2.73047

Gang · September 21, 2021, 2:42am

Should I loop through this sub-brik with 3dFWHMx, where the loop index is used to give 3dFWHMx one sub-brick at a time or can I just give it the full 4D-volume?

Do it on the 4D data (residuals for the model).

Should I use any of these option?

Those options do not seem to be applicable in this case.

Should I for structural data still use -ACF or is classical -fwhm sufficient?

ACF is preferable.

Robin · September 23, 2021, 9:16am

Thanks Gang, a final question:

I ran 3dMVM using a merged 4D volume with the vbm-gray-matter-modulated volumes of all 1800 subjects:


3dinfo ../GM_mod_merg.nii.gz 
++ 3dinfo: AFNI version=AFNI_21.0.12 (Feb 25 2021) [64-bit]

Dataset File:    ../GM_mod_merg.nii.gz
Identifier Code: NII_FIGEgKOgUVe0sJpMi6fe6A  Creation Date: Thu Sep 23 11:08:02 2021
Template Space:  TLRC
Dataset Type:    Echo Planar (-epan)
Byte Order:      LSB_FIRST {assumed} [this CPU native = LSB_FIRST]
Storage Mode:    NIFTI
Storage Space:   6,498,928,800 (6.5 billion) bytes
Geometry String: "MATRIX(2,0,0,-90,0,-2,0,126,0,0,2,-72):91,109,91"
Data Axes Tilt:  Plumb
Data Axes Orientation:
  first  (x) = Right-to-Left
  second (y) = Posterior-to-Anterior
  third  (z) = Inferior-to-Superior   [-orient RPI]
R-to-L extent:   -90.000 [R] -to-    90.000 [L] -step-     2.000 mm [ 91 voxels]
A-to-P extent:   -90.000 [A] -to-   126.000 [P] -step-     2.000 mm [109 voxels]
I-to-S extent:   -72.000 [I] -to-   108.000 [S] -step-     2.000 mm [ 91 voxels]
Number of time steps = 1800  Time step = 2.00000s  Origin = 0.00000s
  -- At sub-brick #0 '?' datum type is float
  -- At sub-brick #1 '?' datum type is float
  -- At sub-brick #2 '?' datum type is float
** For info on all 1800 sub-bricks, use '3dinfo -verb' **

As you can see it has 1800 sub-briks.

I then refered to each of these sub-briks in the mvm-table (yes, they are in the correct order):


[robka@localhost mvm_test]$ cat table_no_smooth.txt 
Subj PRS status InputFile
1001008 HighPRS status1../GM_mod_merg.nii.gz[0]
1001794 HighPRS status2 ../GM_mod_merg.nii.gz[1]
1005915 HighPRS status2 ../GM_mod_merg.nii.gz[2]
etc..

And I ran 3dMVM:


3dMVM -prefix results_no_smooth -dataTable @table_no_smooth.txt -mask ../GM_mask.nii.gz -resid resid_no_smooth -jobs 4 \
-bsVars "PRS*status"

I just don’t understand the output. Usually I get a “results_no_smooth+tlrc.” file(s). With 3 sub-briks: Factor1(PRS), Factor2(status) Factor3(Interaction) and the Intercept (4 sub-briks).
But when I run the above the resulting “results_no_smooth+tlrc.” file(s) have 1804 sub-briks. I.e. all the subjects plus the stats. This has never happaned before. The only new thing is that I refer to a large 4D volume containing all data, instead of individual files like I usually do with smaller samples.


3dinfo results_no_smooth+tlrc.
++ 3dinfo: AFNI version=AFNI_21.0.12 (Feb 25 2021) [64-bit]

Dataset File:    results_no_smooth+tlrc
Identifier Code: XYZ_DvVyj6f44Rr5K15d5G0OjC  Creation Date: 
Template Space:  TLRC
Dataset Type:    Func-Bucket (-fbuc)
Byte Order:      LSB_FIRST [this CPU native = LSB_FIRST]
Storage Mode:    BRIK
Storage Space:   3,256,685,432 (3.3 billion) bytes
Geometry String: "MATRIX(2,0,0,-90,0,-2,0,126,0,0,2,-72):91,109,91"
Data Axes Tilt:  Plumb
Data Axes Orientation:
  first  (x) = Right-to-Left
  second (y) = Posterior-to-Anterior
  third  (z) = Inferior-to-Superior   [-orient RPI]
R-to-L extent:   -90.000 [R] -to-    90.000 [L] -step-     2.000 mm [ 91 voxels]
A-to-P extent:   -90.000 [A] -to-   126.000 [P] -step-     2.000 mm [109 voxels]
I-to-S extent:   -72.000 [I] -to-   108.000 [S] -step-     2.000 mm [ 91 voxels]
Number of values stored at each pixel = 1804
  -- At sub-brick #0 '(Intercept) F' datum type is short:            0 to         32767 [internal]
                                        [*   0.00305185]             0 to           100 [scaled]
     statcode = fift;  statpar = 1 1794
  -- At sub-brick #1 'PRS F' datum type is short:            0 to         32767 [internal]
                                [*  0.000654147]             0 to       21.4345 [scaled]
     statcode = fift;  statpar = 1 1794
  -- At sub-brick #2 'MDDstatus F' datum type is short:            0 to         32767 [internal]
                                      [*  0.000365401]             0 to       11.9731 [scaled]
     statcode = fift;  statpar = 2 1794
** For info on all 1804 sub-bricks, use '3dinfo -verb' **

Is this normal or something to worry about?
Thanks!

Gang · September 25, 2021, 5:25pm

when I run the above the resulting “results_no_smooth+tlrc.” file(s) have 1804 sub-briks. I.e. all the subjects plus the stats.

Do you mean that you got two separate output files: one with those statistical results while the other contains statistical results plus the residuals? If so, that is strange. I never tested the case with a 4D file as an input. So, if the problem persists, I may have to ask you to share a small test dataset to diagnose the problem.

Robin · September 25, 2021, 7:00pm

Hi!

To clarify:
When running the command I get:
results+tlrc.BRIK and results+tlrc.HEAD (expected)
residual+tlrc.BRIK and residuals+tlrc.HEAD (expected)

The residuals file has 1800 sub-briks (expected since we have 1800 subjects and we get one sub-brik in resduals+tlrc per subject)
The results file has 1804 sub-briks (not expected).
After comparing I now see that the first 4 sub-briks are the stats, the next 1800 sub-briks are the residuals. So it seems the program spits out my residuals+tlrc. output as well as adds them into the results+tlrc file.

Is that expected?

Attaching image of viewer to show the sub-briks and an image showing:
left: Sub-brik[0] of the residuals file (residuals for subject0)
right Sub-brik[4] of the results file (also residuals for subject0)

Only two things are new, since I have never seen this:
a) Reffering to a 4D using brackets [] to refer to the correct sub-brik of the input-data (1800 sub-briks for 1800 subjects)
b) Using -resid to get residuals

a) or b) or a+b maybe causes this?

Gang · September 26, 2021, 1:27am

OK, I suspect that your AFNI version is a little old. This was a problem at some point, but it got fixed later. On the other hand, practically you can just ignore or delete those extra 1,800 sub-bricks in results+tlrc