Leave-one-out cross-validation of brain-behaviour correlation clusters

AFNI version: AFNI_23.3.01 'Septimius Severus'

Hi everyone! I have a dataset that we used to create brain-behaviour correlation maps for a contrast of interest in our study. Given that we have a sample size of 20, we want to ensure that the significant correlation clusters we found are actually stable. To do so, we want to use a leave-one-out cross-validation or a similar approach. So far, I've tried doing this manually by first generating several aggregate datasets (each with a unique subject left out) and then using 3dTcorr1D to create a leave-one-out correlation map for the retained subjects. My problem after that becomes, how do I see if this new N-1 dataset matches the left-out subject? I've tried assessing this manually by comparing the voxel overlap but I'm getting results that don't make a ton of sense.

I know that leave-one-out cross-validation is typically done by using a predictive model and predicting the test dataset (the subject that was left out), but I can't seem to find a way to do this in AFNI. Maybe I need to use 3dsvm? But that's a rabbit hole I'm not super familiar with...

In essence, I want to determine whether the significant brain-behaviour clusters we find in the full dataset are spatially stable when a given subject is removed. Any help would be appreciated!

Here is the code I've tried so far:


# Set top directory and go to it
set topdir = /mnt/c/Users/Brady\ Roberts/Desktop/fMRIdata/Correlation_Files/LOO_Overlap_Analyses
cd topdir




##################################### Start with Draw - Write brain-behaviour Correlation ########################################################################################################

# Set contrast name in the name of the overall bucket file for each contrast
set contrast_label = DW 
# Sub-brick for this contrast = 19

# Setup output folder for this particular contrast, then go to it
set contrast_folder_name = "LOO_${contrast_label}"

# Check if the folder exists
if (! -e $contrast_folder_name) then
    echo "Creating folder: $contrast_folder_name"
    mkdir $contrast_folder_name
else
    echo "Folder already exists: $contrast_folder_name"
endif

set outputdir = /mnt/c/Users/Brady\ Roberts/Desktop/fMRIdata/Correlation_Files/LOO_Overlap_Analyses/${contrast_folder_name}
cd outputdir



## Subject MM03 left-out #############################################################################################################################################################

# Set subject to leave out 
set left_out = MM03
set row_to_remove = 1
# Sub-brick to leave out of the 3dTcat command below: 0

# Setup output folder for the subject we are leaving out
cd outputdir
set LOO_folder_name = "LOO_${contrast_label}_${left_out}"

# Check if the folder exists
if (! -e $LOO_folder_name) then
    echo "Creating folder: $LOO_folder_name"
    mkdir $LOO_folder_name
else
    echo "Folder already exists: $LOO_folder_name"
endif

cd topdir

# Duplicate the all-subjects behavioural file and then remove subject's corresponding row from the copied file
cp ${contrast_label}_Values_for_Correlations.1D ${contrast_folder_name}/${LOO_folder_name}/${contrast_label}_Behavioural_Values_for_Correlations_except_${left_out}.1D
sed -i "${row_to_remove}d" ${contrast_folder_name}/${LOO_folder_name}/${contrast_label}_Behavioural_Values_for_Correlations_except_${left_out}.1D

# Copy all sub-bricks (in this case, contrast-specific brain data from subjects) from a prior full data set that was created, except for the sub-bricks corresponding to the subject being left out of this iteration
3dTcat -prefix ${contrast_folder_name}/${LOO_folder_name}/all_${contrast_label}_BRIKs_except_${left_out}+tlrc all_DW_BRIKs+tlrc'[1..19]'

# Run correlation between the difference score of two conditions and the activation for the corresponding sub-brik for that contrast in all subjects except the one that was left out
3dTcorr1D -prefix ${contrast_folder_name}/${LOO_folder_name}/BrainBehaviour_Correlations_${contrast_label}_except_${left_out} \
-spearman \
-mask group_mask_.7olap+tlrc \
${contrast_folder_name}/${LOO_folder_name}/all_${contrast_label}_BRIKs_except_${left_out}+tlrc \
${contrast_folder_name}/${LOO_folder_name}/${contrast_label}_Behavioural_Values_for_Correlations_except_${left_out}.1D

# Clusterize to form a mask that only includes the significant voxels at our standard setting of p <.01 and 20 voxel minimum cluster size
3dClusterize                  \
    -nosum \
    -1Dformat \
    -inset ${contrast_folder_name}/${LOO_folder_name}/BrainBehaviour_Correlations_${contrast_label}_except_${left_out}+tlrc.      \
    -idat 0                   \
    -ithr 0                    \
    -mask group_mask_.7olap+tlrc.     \
    -NN 2                      \
    -bisided p=0.01           \
    -clust_nvox 20           \
    -orient LPI \
    -outvol_if_no_clust \
    -pref_dat ${contrast_folder_name}/${LOO_folder_name}/${contrast_label}_Cluster_Mask_except_${left_out}+tlrc.


# This funciton returns the number of non-zero voxels in the input datasets that OVERLAP.
# It does NOT tell us any information about the overlap in correlation intensity etc. 
# Can compare the returned value with 106, which is the number of voxels in the mask with significant clusters found across ALL subjects.
3dOverlap DW_Cluster_Mask_ALL_SUBJECTS+tlrc. ${contrast_folder_name}/${LOO_folder_name}/${contrast_label}_Cluster_Mask_except_${left_out}+tlrc.


# Now imagine that this code repeated for the other 19 subjects in the dataset....
# In the last line of the above code, I get output corresponding to the overlapping voxels from the overall data set relative to the dataset with a particular subject left out. This is a rough metric of cluster stability and gives me some results that I wouldn't expect (like terribly low overlap).

# Another approach I tried is to create a visual conjunction analysis to see how many of the 20 leave-one-out datasets overlap with my full dataset. I followed this tutorial:
# https://www.youtube.com/watch?v=H4b0pzb7N2s


# Now produce a cluster conjunction map that shows the degree of overlap for the clusters visually

set topdir = /mnt/c/Users/Brady\ Roberts/Desktop/fMRIdata/Correlation_Files/LOO_Overlap_Analyses
cd topdir

# Set contrast name in the name of the overall bucket file for each contrast
set contrast_label = DW

# Setup output folder for this particular contrast, then go to it
set contrast_folder_name = "LOO_${contrast_label}"

3dcalc -prefix Conjunction_Map/Cluster_Conjunction_Map \
-a ${contrast_folder_name}/LOO_${contrast_label}_MM03/${contrast_label}_Cluster_Mask_except_MM03+tlrc. \
-b ${contrast_folder_name}/LOO_${contrast_label}_MM04/${contrast_label}_Cluster_Mask_except_MM04+tlrc. \
-c ${contrast_folder_name}/LOO_${contrast_label}_MM05/${contrast_label}_Cluster_Mask_except_MM05+tlrc. \
-d ${contrast_folder_name}/LOO_${contrast_label}_MM06/${contrast_label}_Cluster_Mask_except_MM06+tlrc. \
-e ${contrast_folder_name}/LOO_${contrast_label}_MM07/${contrast_label}_Cluster_Mask_except_MM07+tlrc. \
-f ${contrast_folder_name}/LOO_${contrast_label}_MM08/${contrast_label}_Cluster_Mask_except_MM08+tlrc. \
-g ${contrast_folder_name}/LOO_${contrast_label}_MM09/${contrast_label}_Cluster_Mask_except_MM09+tlrc. \
-h ${contrast_folder_name}/LOO_${contrast_label}_MM10/${contrast_label}_Cluster_Mask_except_MM10+tlrc. \
-i ${contrast_folder_name}/LOO_${contrast_label}_MM11/${contrast_label}_Cluster_Mask_except_MM11+tlrc. \
-j ${contrast_folder_name}/LOO_${contrast_label}_MM12/${contrast_label}_Cluster_Mask_except_MM12+tlrc. \
-k ${contrast_folder_name}/LOO_${contrast_label}_MM13/${contrast_label}_Cluster_Mask_except_MM13+tlrc. \
-l ${contrast_folder_name}/LOO_${contrast_label}_MM14/${contrast_label}_Cluster_Mask_except_MM14+tlrc. \
-m ${contrast_folder_name}/LOO_${contrast_label}_MM15/${contrast_label}_Cluster_Mask_except_MM15+tlrc. \
-n ${contrast_folder_name}/LOO_${contrast_label}_MM16/${contrast_label}_Cluster_Mask_except_MM16+tlrc. \
-o ${contrast_folder_name}/LOO_${contrast_label}_MM17/${contrast_label}_Cluster_Mask_except_MM17+tlrc. \
-p ${contrast_folder_name}/LOO_${contrast_label}_MM18/${contrast_label}_Cluster_Mask_except_MM18+tlrc. \
-q ${contrast_folder_name}/LOO_${contrast_label}_MM19/${contrast_label}_Cluster_Mask_except_MM19+tlrc. \
-r ${contrast_folder_name}/LOO_${contrast_label}_MM20/${contrast_label}_Cluster_Mask_except_MM20+tlrc. \
-s ${contrast_folder_name}/LOO_${contrast_label}_MM21/${contrast_label}_Cluster_Mask_except_MM21+tlrc. \
-t ${contrast_folder_name}/LOO_${contrast_label}_MM101/${contrast_label}_Cluster_Mask_except_MM101+tlrc. \
-expr 'step(a) + \
    2*step(b) + \
    4*step(c) + \
    8*step(d) + \
    16*step(e) + \
    32*step(f) + \
    64*step(g) + \
    128*step(h) + \
    256*step(i) + \
    512*step(j) + \
    1024*step(k) + \
    2048*step(l) + \
    4096*step(m) + \
    8192*step(n) + \
    16384*step(o) + \
    32768*step(p) + \
    65536*step(q) + \
    131072*step(r) + \
    262144*step(s) + \
    524288*step(t)'

# This kind of worked but isn't showing up in the AFNI GUI the way I want. I can't get it to show overlap in the specific clusters of interest from the full dataset without it also showing spurious clusters elsewhere.

Does the outcome of the original analysis involving 20 individual participants align with your expectations? If so, apart from the sample size, what is your primary concern?

In the 3dcalc output generated from the 20 leave-one-out results, what aspect are you particularly interested in: the intersection or the overall overlap?

Gang Chen

Thanks for the quick reply, Gang. Yes the original results align with our expectations. The pursuit of cross-validation came at the request of a reviewer. Essentially, they are worried that our sample size means there are spurious findings and they want us to try and see if our results are robust.

For the 3dcalc part, I believe we are interested in the overall overlap. Essentially we are trying to overlay all 20 of the leave-one-out datasets on top of the original cluster mask formed from the full dataset and then be able to visualize the number of leave-one-out datasets that overlap any given voxel within that mask. For instance I could imagine a figure with an overlay containing 4 colors representing if there are 0-4, 5-9, 10-14, or 15-20 leave-one-out datasets overlapping a voxel in the original brain-behaviour correlation cluster mask.

There is concern that our sample size might yield spurious findings, and they are suggesting we investigate the robustness of our results.

It's essential to acknowledge that a single study, regardless of its sample size, may not conclusively settle matters. Reproducibility is better pursued through a series of individual studies.

The request for cross-validation was prompted by a reviewer's feedback.

If I understand correctly, the reviewer seems to be apprehensive about potential Gaussian distribution violations due to the sample size of 20. However, in the realm of FMRI modeling, this concern is generally less critical compared to other challenges, especially if not explicitly geared towards prediction. For example, see discussion here and here. If you opt for a leave-one-out strategy to address Gaussian distribution concerns, it might be more logical to present the distribution of effect estimates directly, placing less emphasis on statistical evidence like thresholding or cluster size.

Scientific inquiry should prioritize effect estimation and its associated uncertainty over decision-making based solely on statistical evidence. Moreover, the strength of statistical evidence inherently exists on a continuous spectrum. To effectively communicate the strength of statistical evidence, the highlight-but-don't-hide method offers a superior approach compared to traditional methods that solely focus on statistically significant results. Clusterization through arbitrary thresholding, commonly used in neuroimaging inference, lacks substantial neurological meaning. Therefore, the lack of interpretable outcomes from the intersection of those 20 leave-one-out results is not surprising.

Gang Chen

Thanks, Gang. Based on your comment I was thinking of creating histograms to show the distribution of effect sizes in the leave-one-out datasets at the peak voxel coordinate we report in the paper. Is that what you were getting at when you suggested to present the distribution of effect estimates directly?

Yes, that's essentially what I was proposing. You could simply calculate the mean (or median) and standard error of the histogram. The leave-one-out approach suggested by the reviewer is likely akin to those conventional nonparametric methods like bootstrapping and permutations. I suspect these alternatives would have minimal impact on the overall results. In contrast, the highlight-but-don't-hide method in result reporting would provide more informative and transparent insights.

Gang Chen

Thanks for the clarification, Gang. I think to satisfy the reviewer we will go with the histogram route for this particular inquiry. Coincidentally, I had already read and cited the highlight-but-don't-hide article in our new paper anyways and based on it (and conversations I had with Paul Taylor) we include figures with transparent thresholding. Thanks for your help!