Add thresholded dual regression as an option

I.e. the correction proposed in the NAF2 paper.

IO/loadDualRegression.m

-function [ P, A ] = loadDualRegression( dualregDir, params )
+function [ P, A, Pt, At ] = loadDualRegression( dualregDir, params )
 
 % Subject maps / timecourses
-P = cell(params.S, 1);
-A = cell(params.S, 1);
+% And thresholded versions too
+P  = cell(params.S, 1);
+Pt = cell(params.S, 1);
+A  = cell(params.S, 1);
+At = cell(params.S, 1);
 sr = 0;
 for s = 1:params.S
-    P{s} = 0.0;
-    A{s} = cell(params.R(s), 1);
+    P{s}  = 0.0;
+    Pt{s} = 0.0;
+    A{s}  = cell(params.R(s), 1);
+    At{s} = cell(params.R(s), 1);
     for r = 1:params.R(s)
         subj = sprintf('subject%05d', sr);
         
@@ -14,16 +19,22 @@ for s = 1:params.S
         Psr = read_avw(fullfile( ...
             dualregDir, ['dr_stage2_' subj '.nii.gz']));
         P{s} = P{s} + reshape(Psr, params.V, params.iN);
+        Psr = read_avw(fullfile( ...
+            dualregDir, ['dr_stage4_' subj '_thresh.nii.gz']));
+        Pt{s} = Pt{s} + reshape(Psr, params.V, params.iN);
         
         % Timecourses
         A{s}{r} = load(fullfile( ...
             dualregDir, ['dr_stage1_' subj '.txt']))';
+        At{s}{r} = load(fullfile( ...
+            dualregDir, ['dr_stage4_' subj '.txt']))';
         
         sr = sr + 1;
     end
 
     % Take mean map
-    P{s} = P{s} ./ params.R(s);
+    P{s}  = P{s} ./ params.R(s);
+    Pt{s} = Pt{s} ./ params.R(s);
 end
 
 end

Methods/DualRegression.sh

@@ -17,7 +17,7 @@ des_norm=1
 design=-1
 n_perm=0
 dual_regression "${group_maps}" \
-    ${des_norm} ${design} ${n_perm} \
+    ${des_norm} ${design} ${n_perm} --thr \
     "${output_dir}" \
     $(cat "${nifti_dir}/DualRegression_SpecFile.txt") \
     > "${output_dir}/TerminalOutput.txt"

RunTests.m

@@ -371,13 +371,18 @@ for n = 1:params.nRepeats
             'sh Methods/DualRegression.sh %s %s %s', ...
             dualregDir, niftiDir, fullfile(melodicDir, 'melodic_IC.nii.gz')));
         
+        % `t` means thresholded as per Bijsterbosch et al., eLife, 2019
         icadr_Pg = loadMELODIC(melodicDir, params);
-        [icadr_P, icadr_A] = loadDualRegression(dualregDir, params);
+        [icadr_P, icadr_A, icadr_Pt, icadr_At] ...
+            = loadDualRegression(dualregDir, params);
         %[icadr_P, icadr_A] = runDR(D, icadr_Pg, params);
-        icadr_pcA = calculateNetmats(icadr_A, params);
+        icadr_pcA  = calculateNetmats(icadr_A, params);
+        icadr_pcAt = calculateNetmats(icadr_At, params);
         
         scores.ICA_DR(n) = calculateDecompositionAccuracy( ...
                 P, icadr_P, A, icadr_A, pcA, icadr_pcA, params);
+        scores.ICA_DRt(n) = calculateDecompositionAccuracy( ...
+                P, icadr_Pt, A, icadr_At, pcA, icadr_pcAt, params);
     catch
         warning('MELODIC run failed.');
     end