Auto stash before merge of "master" and "origin/master"

NiiStat.m

@@ -84,7 +84,7 @@ function NiiStat(xlsname, roiIndices, modalityIndices,numPermute, pThresh, minOv
 if exist(xlsname,'file') ~= 2
     error('Unable to find Excel file named %s\n',xlsname);
 end
-[designMat, designUsesNiiImages, minOverlapValFile] = nii_read_design (xlsname);
+[designMat, designUsesNiiImages, minOverlapValFile, nuisanceMat] = nii_read_design (xlsname);
 if ~exist('minOverlap','var')
     if isempty(minOverlapValFile)
         minOverlap = 0;
@@ -132,7 +132,7 @@ function NiiStat(xlsname, roiIndices, modalityIndices,numPermute, pThresh, minOv
         'Minimum overlap (1..numSubj):',...
         ['ROI (0=voxels ' sprintf('%s',kROInumbers) ' negative for correlations [multi OK]'],...
         ['Modality (' sprintf('%s',kModalityNumbers) ') [multiple OK]'],...
-        'Special (1=explicit voxel mask, 2=regress lesion volume, 3=de-skew, 4=include WM/CSF connectivity, 5=customROI, 6=TFCE, 7=reportROImeans, 8=SVM, 9=LowRes, 10=LH only, 11=RH only; 12=interhemispheric) [multi OK]',...
+        'Special (1=explicit voxel mask, 2=control for lesion volume, 3=de-skew, 4=include WM/CSF connectivity, 5=customROI, 6=TFCE, 7=reportROImeans, 8=SVM, 9=LowRes, 10=LH only, 11=RH only; 12=interhemispheric) [multi OK]',...
         'Statistics name [optional]'
         };
     dlg_title = ['Options for analyzing ' xlsname];
@@ -258,7 +258,7 @@ function NiiStat(xlsname, roiIndices, modalityIndices,numPermute, pThresh, minOv
         end
         fprintf('Analyzing roi=%d, modality=%d, permute=%d, %sdesign=%s\n',roiIndex, modalityIndex,numPermute,specialStr, xlsname);
         %Roger added GUI as last argument
-        processExcelSub(designMat, roiIndex, modalityIndex,numPermute, pThresh, minOverlap, regressBehav, maskName, GrayMatterConnectivityOnly, deSkew, customROI, doTFCE, reportROIvalues, xlsname, kROIs, doSVM, doVoxReduce, hemiKey, interhemi, statname,GUI); %%GY
+        processExcelSub(designMat, roiIndex, modalityIndex,numPermute, pThresh, minOverlap, regressBehav, maskName, GrayMatterConnectivityOnly, deSkew, customROI, doTFCE, reportROIvalues, xlsname, kROIs, doSVM, doVoxReduce, hemiKey, interhemi, statname,GUI, nuisanceMat); %%GY
     end
 end
 %end nii_stat_mat()
@@ -333,7 +333,7 @@ function NiiStat(xlsname, roiIndices, modalityIndices,numPermute, pThresh, minOv
 % %end readDesign()
 
 %Roger added GUI as input argument
-function processExcelSub(designMat, roiIndex, modalityIndex,numPermute, pThresh, minOverlap, regressBehav, mask_filename, GrayMatterConnectivityOnly, deSkew, customROI, doTFCE, reportROIvalues, xlsname, kROIs, doSVM, doVoxReduce, hemiKey, interhemi, statname,GUI) %%GY
+function processExcelSub(designMat, roiIndex, modalityIndex,numPermute, pThresh, minOverlap, regressBehav, mask_filename, GrayMatterConnectivityOnly, deSkew, customROI, doTFCE, reportROIvalues, xlsname, kROIs, doSVM, doVoxReduce, hemiKey, interhemi, statname,GUI, nuisanceMat) %%GY
 %GrayMatterConnectivityOnly = true; %if true, dti only analyzes gray matter connections
 %kROIs = strvcat('bro','jhu','fox','tpm','aal','catani'); %#ok<*REMFF1>
 %kModalities = strvcat('lesion','cbf','rest','i3mT1','i3mT2','fa','dti','md'); %#ok<REMFF1> %lesion, 2=CBF, 3=rest
@@ -495,6 +495,9 @@ function processExcelSub(designMat, roiIndex, modalityIndex,numPermute, pThresh,
             end
             data.filename = in_filename;
             data.behav = designMat(i); % <- crucial: we inject behavioral data from Excel file!
+            if ~isempty (nuisanceMat)
+                data.nuisance = nuisanceMat(i);
+            end
             subj_data{idx} = data; %#ok<AGROW>
         else
             %dat = load (in_filename, subfield);
@@ -525,6 +528,9 @@ function processExcelSub(designMat, roiIndex, modalityIndex,numPermute, pThresh,
                 idx = idx + 1;
                 subj_data{idx}.filename = in_filename; %#ok<AGROW>
                 subj_data{idx}.behav = designMat(i); %#ok<AGROW>
+                if ~isempty (nuisanceMat)
+                    subj_data{idx}.nuisance = nuisanceMat(i);%#ok<AGROW>
+                end
                 if isempty(voxMask)
                     %dat.(ROIfield).mean = normSub(dat.(ROIfield).mean, cbfMean, cbfStd);
                     subj_data{idx}.(ROIfield)  = dat.(ROIfield); %#ok<AGROW>
@@ -580,6 +586,7 @@ function processExcelSub(designMat, roiIndex, modalityIndex,numPermute, pThresh,
 % make sure all the subjects have all numeric fields
 beh = zeros(n_subj,n_beh);
 beh(:) = nan;
+nuisance = [];
 for i = 1 : n_subj
     for j = 1:n_beh %length(beh_names)
         if isfield (subj_data{i}.behav, beh_names{j})
@@ -596,7 +603,14 @@ function processExcelSub(designMat, roiIndex, modalityIndex,numPermute, pThresh,
             disp (['Warning! Subject ' subj_data{i}.filename ' does not have a field ' beh_names{j}]);
         end
     end
+
+    if ~isempty (nuisanceMat)
+        nuisance (i, :) = structfun (@(x) x, subj_data{i}.nuisance); %#ok<AGROW>
+    end
+
 end
+
+
 if regressBehav
     vol = zeros(n_subj,1);
     vol(:) = nan;
@@ -610,23 +624,25 @@ function processExcelSub(designMat, roiIndex, modalityIndex,numPermute, pThresh,
         end;
     end;
     if sum(~isnan(vol(:))) > 1
-        for i = 1:n_beh
-            %beh_names1 = deblank(beh_names(i));
-            beh1 = beh(:,i);
-            good_idx = intersect (find(~isnan(beh1)), find(~isnan(vol)));
-            dat = beh1(good_idx)'; %behavior is the data
-            reg = vol(good_idx)'; %lesion volume is our regressor
-            preSD = std(dat);
-            if ~isnan(std(dat)) && (preSD ~= 0) && (std(reg) ~= 0) %both data and regressor have some variability
-                G = ones (2, size(dat,2)); %constant
-                G (2, :) = reg; % linear trend
-                G_pseudoinv = G' / (G * G'); %aka: G_pseudoinv = G' * inv (G * G');
-                Beta = dat * G_pseudoinv;
-                dat = dat - Beta*G; %mean is zero
-                fprintf('Regressing %s with lesion volume reduces standard deviation from %f to %f\n',char(deblank(beh_names(i))),preSD, std(dat) );
-                beh(good_idx,i) = dat;
-            end
-        end
+        nuisance = [nuisance vol];
+
+%         for i = 1:n_beh
+%             %beh_names1 = deblank(beh_names(i));
+%             beh1 = beh(:,i);
+%             good_idx = intersect (find(~isnan(beh1)), find(~isnan(vol)));
+%             dat = beh1(good_idx)'; %behavior is the data
+%             reg = vol(good_idx)'; %lesion volume is our regressor
+%             preSD = std(dat);
+%             if ~isnan(std(dat)) && (preSD ~= 0) && (std(reg) ~= 0) %both data and regressor have some variability
+%                 G = ones (2, size(dat,2)); %constant
+%                 G (2, :) = reg; % linear trend
+%                 G_pseudoinv = G' / (G * G'); %aka: G_pseudoinv = G' * inv (G * G');
+%                 Beta = dat * G_pseudoinv;
+%                 dat = dat - Beta*G; %mean is zero
+%                 fprintf('Regressing %s with lesion volume reduces standard deviation from %f to %f\n',char(deblank(beh_names(i))),preSD, std(dat) );
+%                 beh(good_idx,i) = dat;
+%             end
+%         end
     end
 end %if regressBehav - regress behavioral data using lesion volume
 
@@ -1036,10 +1052,10 @@ function processExcelSub(designMat, roiIndex, modalityIndex,numPermute, pThresh,
         %end
 
         if doSVM
-            nii_stat_svm(les1, beh1, beh_names1,statname, les_names, subj_data, roiName, logicalMask1, hdr, pThresh, numPermute);
+            nii_stat_svm(les1, beh1, beh_names1,statname, les_names, subj_data, roiName, logicalMask1, hdr, pThresh, numPermute, nuisance);
         else
             %nii_stat_core(les1, beh1, beh_names1,hdr, pThresh, numPermute, minOverlap,statname, les_names,hdrTFCE, voxMask);
-            nii_stat_core(les1, beh1, beh_names1,hdr, pThresh, numPermute, logicalMask1,statname, les_names,hdrTFCE);
+            nii_stat_core(les1, beh1, beh_names1,hdr, pThresh, numPermute, logicalMask1,statname, les_names,hdrTFCE, nuisance);
         end
 
 %         diary off
@@ -1054,9 +1070,9 @@ function processExcelSub(designMat, roiIndex, modalityIndex,numPermute, pThresh,
     %les_names(2:2:end)=[]; % Remove even COLUMNS: right in AALCAT: analyze left
     %les(:,2:2:end)=[]; % Remove even COLUMNS: right in AALCAT: analyze left
     if doSVM
-        nii_stat_svm(les, beh, beh_names, statname, les_names, subj_data, roiName, logicalMask, hdr, pThresh, numPermute);
+        nii_stat_svm(les, beh, beh_names, statname, les_names, subj_data, roiName, logicalMask, hdr, pThresh, numPermute, nuisance);
     else
-        nii_stat_core(les, beh, beh_names,hdr, pThresh, numPermute, logicalMask,statname, les_names, hdrTFCE);
+        nii_stat_core(les, beh, beh_names,hdr, pThresh, numPermute, logicalMask,statname, les_names, hdrTFCE, nuisance);
     end
 end
 
@@ -1131,6 +1147,7 @@ function processExcelSub(designMat, roiIndex, modalityIndex,numPermute, pThresh,
 for i = 1:length (hemi_regexp)
     hemi_idx = [hemi_idx; find(cellfun(@length, regexp (labels, hemi_regexp{i})))];
 end
+hemi_idx = unique (hemi_idx);
 %end extract_hemi_idxSub
 
 

nii_read_design.m

@@ -1,4 +1,4 @@
-function [designMat, designUsesNiiImages, CritN] = nii_read_design (xlsname, worksheetname)
+function [designMat, designUsesNiiImages, CritN, nuisanceMat] = nii_read_design (xlsname, worksheetname)
 designUsesNiiImages = false;
 CritN = [];
 if exist(xlsname,'file') ~= 2
@@ -15,9 +15,9 @@
 if strcmpi(x,'.tab') || strcmpi(x,'.txt')  || strcmpi(x,'.val')
     [dMat, CritN] = nii_tab2mat(xlsname);
 else
-    dMat = nii_xls2mat(xlsname , worksheetname,'', true);
+    [dMat, nuisance] = nii_xls2mat(xlsname , worksheetname,'', true);
     if isempty(dMat)
-        dMat = nii_xls2mat(xlsname , 'Data (2)','', true);
+        [dMat, nuisance] = nii_xls2mat(xlsname , 'Data (2)','', true);
     end
 end
 if isempty(dMat)
@@ -33,6 +33,9 @@
 numNII = 0; %number of NIfTI files
 numMat = 0; %number of Mat files
 numOK = 0;
+if isempty (nuisance)
+    nuisanceMat = [];
+end
 %designMat = [];
 for i=1:size(dMat,2)
     matname = deblank( dMat(i).(SNames{1}));
@@ -49,6 +52,15 @@
         fprintf('Warning: no valid behavioral data for %s\n',matname);
         matname = '';
     end
+    % added by GY: 
+    % all nusiance variables must be present for all participants!
+    % exclude participants with missing nuisance variables (i.e. NaNs)
+    if ~isempty (nuisance) && ~isempty (matname)
+        if sum (structfun (@isnan, nuisance(i))) 
+            fprintf ('Warning: excluding %s because some nuisance regressors are missing\n', matname);
+            matname = '';
+        end
+    end
     if ~isempty(matname)
         [matname] = findMatFileSub(matname,imgpath, xlsname);
         [~, ~, ext] = fileparts(matname);
@@ -61,7 +73,9 @@
             dMat(i).(SNames{1}) = matname;
             numOK = numOK + 1;
             designMat(numOK) = dMat(i); %#ok<AGROW>
-
+            if ~isempty (nuisance)
+                nuisanceMat(numOK) = nuisance(i);%#ok<AGROW>
+            end
         end
     end
 end

nii_stat_core.m

@@ -1,5 +1,4 @@
-%function nii_stat_core(les,beh, beh_names,hdr, kPcrit, kNumRandPerm, kOnlyAnalyzeRegionsDamagedInAtleastNSubjects,statname, roi_names, hdrTFCE, voxMask)
-function nii_stat_core(les,beh, beh_names,hdr, kPcrit, kNumRandPerm, logicalMask,statname, roi_names, hdrTFCE)
+function nii_stat_core(les,beh, beh_names,hdr, kPcrit, kNumRandPerm, logicalMask,statname, roi_names, hdrTFCE, nuisance)
 %Generates statistical tests. Called by the wrappers nii_stat_mat and nii_stat_val
 % les    : map of lesions/regions one row per participant, one column per voxel
 % beh    : matrix of behavior, one row per participant, one column per condition
@@ -21,6 +20,7 @@ function nii_stat_core(les,beh, beh_names,hdr, kPcrit, kNumRandPerm, logicalMask
 % statname : (optional) results saved to disk will include this name. If not provided this becomes 'anonymous'
 % roi_names : (optional) if one per column of les, then a labeled table is produced
 % hdrTFCE   : (optional, voxelwise analyses only) - dimensions of volume
+% nuisance: (optional) matrix of nuisance variables
 %Examples:
 %     five0five1 = [0 0 0 0 0 1 1 1 1 1]';
 %     ascend1to10 = [1:10]';
@@ -58,6 +58,9 @@ function nii_stat_core(les,beh, beh_names,hdr, kPcrit, kNumRandPerm, logicalMask
 if ~exist('hdrTFCE','var')
     hdrTFCE = [];
 end
+if ~exist('nuisance','var')
+    nuisance = [];
+end
 
 if size(les,1) ~= size(beh,1)
     error('Error: les and beh must have the same number of rows (one image and one set of behavioral data per participant)');
@@ -118,7 +121,7 @@ function nii_stat_core(les,beh, beh_names,hdr, kPcrit, kNumRandPerm, logicalMask
 end
 if (kNumRandPerm < -1) && (size(beh,2) <= 1) %special case: nuisance regressors which requires multiple behaviors
     error('Error: only one behavioral variable provided: unable to control for nuisance regressors (solution: either provide more behavioral data or specify positive number of permutations');
-elseif (kNumRandPerm < -1) && (size(beh,2) > 1) %special case: nuisance regressors.
+elseif (kNumRandPerm < -1) && (size(beh,2) > 1) %special case: nuisance regressors with permutations.
     X1 = [beh, ones(numObs,1)];%add constant column for intercept
     global global_flContrast
     if isempty(global_flContrast)
@@ -159,9 +162,18 @@ function nii_stat_core(les,beh, beh_names,hdr, kPcrit, kNumRandPerm, logicalMask
 else %behavior and/or lesions is continuous
 	fprintf('Computing glm (pooled-variance t-test, linear regression) for %d regions/voxels with analyzing %d behavioral variables (positive Z when increased image brightness correlates with increased behavioral score).\n',length(good_idx),size(beh,2));
 
-    for i = 1: size(beh,2)
-        [z(:,i), threshMin(i), threshMax(i)] = glm_permSub(les,beh(:,i), kNumRandPerm, kPcrit, good_idx, hdrTFCE);
-    end;
+    % GY, Aug 2019: nuisance regressors
+    if isempty (nuisance)
+        % the good old way
+        for i = 1: size(beh,2)
+            [z(:,i), threshMin(i), threshMax(i)] = glm_permSub(les,beh(:,i), kNumRandPerm, kPcrit, good_idx, hdrTFCE);
+        end;
+    else
+        % the new way, with nuisance regressors
+        for i = 1: size(beh,2)
+            [z(:,i), threshMin(i), threshMax(i)] = glm_permSub_nuisance(les,beh(:,i), nuisance, kNumRandPerm, kPcrit, good_idx, hdrTFCE);
+        end;    
+    end
 end
 %global global_powerMap %option to save parameters for power analysis
 %if ~isempty(global_powerMap) && global_powerMap
@@ -662,6 +674,86 @@ function saveStatMapSub(hdr,statImg, statName,  minThreshold, maxThreshold)
 threshMax = spm_t2z(threshMax,df); %report Z scores so DF not relevant
 %end glm_permSub()
 
+function [uncZ, threshMin, threshMax] = glm_permSub_nuisance(Y, X, nuisance, nPerms, kPcrit, good_idx, hdrTFCE)
+%returns uncorrected z-score for all voxels in Y given single column
+%predictor X and a matrix of nuisance regressors in Z
+% Y: volume data
+% X: single column predictor
+%if either X or Y is binomial, results are pooled variance t-test, else correlation coefficient
+% nPerms: Number of permutations to compute
+% kPcrit: Critical threshold
+%Example
+% glm_t([1 1 0 0 0 1; 0 0 1 1 1 0]',[1 2 3 4 5 6]') %pooled variance t-test
+%
+%inspired by Ged Ridgway's glm_perm_flz
+if ~exist('nPerms','var')
+    nPerms = 0;
+end
+if ~exist('kPcrit','var')
+    kPcrit = 0.05;
+end
+% Basics and reusable components
+[n f] = size(X); %rows=observations, columns=factors
+[nY v] = size(Y); %#ok<NASGU> v is number of tests/voxels
+if (f ~= 1), error('glm_permSub is for one column of X at a time (transpose?)'); end;
+if nY ~= n, error('glm_permSub X and Y data sizes are inconsistent'); end
+X = [X nuisance ones(size(X,1),1)];
+c = [1 zeros(1, size(nuisance, 2)) 0]'; %contrast (incl. nuisance regressors)
+df = n - rank(X);
+pXX = pinv(X)*pinv(X)'; % = pinv(X'*X), which is reusable, because
+pX  = pXX * X';         % pinv(P*X) = pinv(X'*P'*P*X)*X'*P' = pXX * (P*X)'
+% Original design (identity permutation)
+t = glm_quick_t(Y, X, pXX, pX, df, c);
+if any(~isfinite(t(:)))
+    warning('glm_permSub zeroed NaN t-scores'); %CR 3Oct2016
+    t(~isfinite(t))= 0; %CR patch
+end
+uncZ = spm_t2z(t,df); %report Z scores so DF not relevant
+if nPerms < 2
+    threshMin = -Inf;
+    threshMax = Inf;
+    return
+end
+if numel(hdrTFCE) == 3
+   img = zeros(hdrTFCE);
+   img(good_idx) = t;
+   img = tfceMex(img, 100);
+   t = img(good_idx);
+   uncZ = spm_t2z(t,df); %report Z scores so DF not relevant
+end
+% Things to track over permutations
+peak = nan(nPerms,1);
+nadir = nan(nPerms,1);
+peak(1) = max(t);
+nadir(1) = min(t);
+perm5pct = round(nPerms * 0.05);
+startTime = tic;
+for p = 2:nPerms
+    if p == perm5pct
+        fprintf('Expected permutation time is %f seconds\n',toc(startTime)*20)
+    end
+    Xp  = X(randperm(n), :);
+    pXX = pinv(Xp)*pinv(Xp)'; %??? supposedly not require- reusable?
+    pXp = pXX * Xp'; % = pinv(Xp)
+    tp  = glm_quick_t(Y, Xp, pXX, pXp, df, c);
+    if numel(hdrTFCE) == 3
+       img = zeros(hdrTFCE);
+       img(good_idx) = tp;
+       img = tfceMex(img, 100);
+       tp = img(good_idx);
+    end
+    peak(p) = max(tp(:));
+    nadir(p) = min(tp(:));
+end
+threshMin = permThreshLowSub (nadir, kPcrit);
+threshMax = permThreshHighSub (peak, kPcrit);
+threshMin = spm_t2z(threshMin,df); %report Z scores so DF not relevant
+threshMax = spm_t2z(threshMax,df); %report Z scores so DF not relevant
+%end glm_permSub()
+
+
+
+
 function savePowerSub(les, beh ,good_idx, hdr, roi_names, beh_names)
 %save components required for a power analysis
 m.les = les;