Merge branch 'develop2'

FinalCode/BootStrapFunction.m

@@ -1,5 +1,7 @@
 function BootStrap = BootStrapFunction(data,Nboot,FieldNames)
-PointEstResults = WIPsubfnFitModel(data);
+
+PointEstResults  = FitProcessModel(data);
+
 % initialize bootstrap values based on the size of the results structure.
 % FieldNames to bootstrap
 
@@ -38,7 +40,7 @@
         Samp = [Samp1; Samp2];
     end
     temp.data = temp.data(Samp,:);
-    Results = WIPsubfnFitModel(temp);
+    Results = FitProcessModel(temp);
     BootStrap.beta(:,:,i) = Results.beta;
     BootStrap.B(:,:,i) = Results.B;
     BootStrap.Paths{i} = Results.Paths;

UnusedCode/WIPsubfnConvertPathsToMatrix.m → FinalCode/ConvertPathsToMatrix.m

@@ -1,4 +1,4 @@
-function data = WIPsubfnConvertPathsToMatrix(Paths,PathNumber)
+function data = ConvertPathsToMatrix(Paths,PathNumber)
 % number of bootstraps
 Nrepeats = size(Paths,2);
 

FinalCode/CreateBCaCI.m

@@ -0,0 +1,66 @@
+function BCaCI = CreateBCaCI(Results,BootStrap,JackKnife,Thresholds)
+% from the point estimate, the bootstrap and the jacknife results return
+% all of the bias-corrected, accelerated confidence intervals.
+
+% Therefore, the whatever bootstrap estimates are calculated this function
+% creates the BCa confidence intervals for it.
+
+% Find the parameters in the bootstrap
+FieldNames = fieldnames(BootStrap);
+
+Nthresh = length(Thresholds);
+NPaths = size(Results.Paths,1);
+% prepare the structure to hold all of the BCa confidence intervals
+BCaCI = {};
+for i = 1:length(FieldNames)
+    Value = getfield(Results,FieldNames{i});
+    if iscell(Value)
+        BlankValue = zeros(numel(Results.Paths{1}),1,2,NPaths,Nthresh);
+    else
+        BlankValue = zeros([size(Value) 2 Nthresh]);
+    end
+    BCaCI = setfield(BCaCI,FieldNames{i},BlankValue);
+end
+%% put the confidence intervals into this structure
+% cycle over the input fieldnames
+for i = 1:length(FieldNames)
+    if ~isempty(strmatch(FieldNames{i},'Paths'))
+        CurrentBCaCI = zeros(size(getfield(BCaCI,FieldNames{i})));
+        % cycle over all the input thresholds
+        for t = 1:Nthresh
+            CurrentAlpha = Thresholds(t);
+            for j = 1:NPaths
+                % This takes that date from each path and flattens it
+                BootStrapData = ConvertPathsToMatrix(BootStrap.Paths,j);
+                JackKnifeData = ConvertPathsToMatrix(JackKnife.Paths,j);
+
+                PointEstimate = reshape(Results.Paths{j},numel(Results.Paths{j}),1);
+
+                % The flattended data allows it to be used by this
+                % program which calculates the adjusted alpha limits used for
+                % determining the confidence intervals
+                [Alpha1 Alpha2 Z p] = CalculateBCaLimits(JackKnifeData,PointEstimate, BootStrapData,CurrentAlpha);
+
+                % find the confidence intervals for these adjusted
+                % alpha limits
+                CurrentBCaCI(:,:,:,j,t) = CalculateBCaCI(BootStrapData,Alpha1,Alpha2,PointEstimate);
+            end
+        end
+        % put the confidence interval data back into the structure
+        % TO DO: somehow UNFLATTEN the data
+        BCaCI = setfield(BCaCI,FieldNames{i},CurrentBCaCI);
+    else
+        CurrentBCaCI = zeros(size(getfield(BCaCI,FieldNames{i})));
+        for t = 1:Nthresh
+            CurrentAlpha = Thresholds(t);
+            BootStrapData = getfield(BootStrap,FieldNames{i});
+            JackKnifeData = getfield(JackKnife,FieldNames{i});
+            PointEstimate = getfield(Results,FieldNames{i});
+            [Alpha1 Alpha2 Z p] = CalculateBCaLimits(JackKnifeData,PointEstimate, BootStrapData,CurrentAlpha);
+
+            % find the confidence intervals
+            CurrentBCaCI(:,:,:,t) = CalculateBCaCI(BootStrapData,Alpha1,Alpha2,PointEstimate);
+        end
+        BCaCI = setfield(BCaCI,FieldNames{i},CurrentBCaCI);
+    end
+end

FinalCode/CycleOverVoxelsProcessBootstrap.m

@@ -8,6 +8,7 @@
 % Prepare the output structure
 Results = cell(Nvoxels,1);
 for i = 1:Nvoxels
+    fprintf(1,'%d of %d voxels\n',i,Nvoxels);
     % Extract the data for this voxel
     OneVoxelModel = ExtractDataFromVoxel(ModelInfo,i);
     % Perform the analysis for this voxel

FinalCode/ExampleSetup.m

@@ -51,6 +51,8 @@
 BEHAVIOR = randn(Nsub,1);
 AgeGroup = round(rand(Nsub,1));
 
+COVARIATES = randn(Nsub,2);
+
 fprintf(1,'Done preparing data.\n');
 
 %% General Setup
@@ -85,7 +87,7 @@
 % correctsed accelerated confidence intervals determined from bootstrap
 % resampling.
 % Are there any voxel-wise bootstrap resamplings?
-Nboot = 5;
+Nboot = 100;
 
 % An alternative to voxel-wise statistics is a map-wise statistics based
 % off of the maximum statistic approach from a series of permutation
@@ -144,11 +146,20 @@
 DataHeader.dt = [16 0];
 ModelInfo.DataHeader = DataHeader;
 
-%% Model Setup
+%% Model Specific Setup Number One
+% THe first model is to perform a simple mediation analysis with the
+% voxel-wise brain data as the mediation between a dichotomous variable
+% (Age group) and a continuous variable (Behavior). 
+
+%     M
+%    / \
+%   /   \
+%  X     Y
+
 Model1 = ModelInfo;
 
 % Name of the output folder
-Tag = 'ExampleModel1_perm';
+Tag = 'ExampleModel1_boot';
 Model1.Tag = Tag;
 
 % The first model is a basic mediation model testing whether the effect of age
@@ -167,7 +178,7 @@
 % VAR_2 = VAR_1
 % VAR_3 = VAR_1 + VAR_2
 
-Direct = zeros(Nvar);
+Direct = zeros(Model1.Nvar);
 Direct(1,[2 3]) = 1;
 Direct(2,3) = 1;
 
@@ -181,28 +192,137 @@
 
 % For this model there are no interactions so this is kept as a zero
 % matrix.
-Inter = zeros(Nvar);
+Inter = zeros(Model1.Nvar);
 
 % The paths to be tested are included in another matrix the same size of
 % the DIRECT matrix. One difference here is that the PATHS matrix may have
 % a third dimension to test multiple paths within a single model. The steps
 % along a path are specified with integers, i.e. step 1, step 2 ...
-Paths = zeros(Nvar);
+Paths = zeros(Model1.Nvar);
 Paths(1,2) = 1;
 Paths(2,3) = 2;
 
 
 Model1.Direct = Direct;
 Model1.Inter = Inter;
 Model1.Paths = Paths;
-%%
+
+%% Run the analysis
 ResultsFolder = PrepareDataForProcess(Model1);
 % The writing of the images should ideally be perfomed by the cluster once
 % the analyses have completed. Therefore, a check is needed or better yet a
 % wait command for the cluster job: 
 % e.g. qsub -hold_jid job1,job2 -N job3 ./c.sh
 WriteOutResults(ResultsFolder)
 
+
+%% Model Specific Setup Number Two
+% The second model is to perform a simple mediation analysis with the
+% voxel-wise brain data as the mediation between a dichotomous variable
+% (Age group) and a continuous variable (Behavior). This analysis now
+% includes covariates also.
+
+
+Model2 = ModelInfo;
+
+% Add the covariate names
+Model2.Names = [Model2.Names 'Cov1' 'Cov2'];
+
+% Add the covariate data
+Model2.data = [Model2.data COVARIATES(:,1) COVARIATES(:,2)];
+
+% Update the number of variables value 
+Model2.Nvar = length(Model2.Names);
+
+% Name of the output folder
+Tag = 'ExampleModel2_boot';
+Model2.Tag = Tag;
+
+% Specify the simple mediation model as above
+
+Direct = zeros(Model2.Nvar);
+Direct(1,[2 3]) = 1;
+Direct(2,3) = 1;
+
+% Now specify the covariates in the models
+Direct([4,5],2) = 1;
+Direct([4,5],3) = 1;
+
+% For this model there are no interactions so this is kept as a zero
+% matrix.
+Inter = zeros(Model2.Nvar);
+
+% The paths do not change
+Paths = zeros(Model2.Nvar);
+Paths(1,2) = 1;
+Paths(2,3) = 2;
+
+
+Model2.Direct = Direct;
+Model2.Inter = Inter;
+Model2.Paths = Paths;
+
+ResultsFolder = PrepareDataForProcess(Model2);
+
+% The writing of the images should ideally be perfomed by the cluster once
+% the analyses have completed. Therefore, a check is needed or better yet a
+% wait command for the cluster job: 
+% e.g. qsub -hold_jid job1,job2 -N job3 ./c.sh
+WriteOutResults(ResultsFolder)
+
+
+%% Model Specific Setup Number Three
+% The third model is to perform a moderated-mediation analysis with the
+% voxel-wise brain data as the mediation between a dichotomous variable
+% (Age group) and a continuous variable (Behavior). 
+% Now there is an interaction between the brain measure and age group in
+% predicting the behavior variable
+
+%     M
+%    / \
+%   /  /\
+%  /  /  \
+% X --    Y
+
+Model3 = ModelInfo;
+
+% Update the number of variables value 
+Model3.Nvar = length(Model3.Names);
+
+% Name of the output folder
+Tag = 'ExampleModel3_boot';
+Model3.Tag = Tag;
+
+% Specify the simple mediation model as above
+
+Direct = zeros(Model3.Nvar);
+Direct(1,[2 3]) = 1;
+Direct(2,3) = 1;
+
+
+% For this model there is interaction
+Inter = zeros(Model3.Nvar);
+Inter([1 2],3) = 1
+
+
+% The paths do not change
+Paths = zeros(Model3.Nvar);
+Paths(1,2) = 1;
+Paths(2,3) = 2;
+
+
+Model3.Direct = Direct;
+Model3.Inter = Inter;
+Model3.Paths = Paths;
+
+ResultsFolder = PrepareDataForProcess(Model3);
+
+% The writing of the images should ideally be perfomed by the cluster once
+% the analyses have completed. Therefore, a check is needed or better yet a
+% wait command for the cluster job: 
+% e.g. qsub -hold_jid job1,job2 -N job3 ./c.sh
+WriteOutResults(ResultsFolder)
+
 %%
 
 % TODO 
@@ -224,22 +344,3 @@
 % transient times of large amounts of data in memory. I am afraid to
 % overload a single computer even transiently.
 % 
-% 
-% 
-% WIPsubfnPrepareDataForProcessPERMUTE
-% % THis performs the permutation testing approach
-% WIPsubfnVoxelWiseProcessPERMUTE(InDataPath,count,Nperm)
-% % This applies the BCaCI to each voxel
-% Results = WIPsubfnProcessData(Model)
-%     
-%     BootStrap = WIPsubfnBootStrap(Model,Model.Nboot,FieldNames);
-% 
-%     % Perform the jack-knife step
-%     JackKnife = WIPJackKnife(Model,Results,FieldNames);
-%     
-%     % Calculate the BCaci values for each parameter
-%     Results.BCaCI = WIPsubfnCreateBCaCI(Results,BootStrap,JackKnife,Model.Thresh);
-% 
-% % THis does the actual regression model fitting.
-% Results = WIPsubfnFitModel(data)
-% 

FinalCode/JackKnifeFunction.m

@@ -22,7 +22,8 @@
     Include = 1:Model.Nsub;
     Include(i) = 0;
     tempModel.data = Model.data(find(Include),:);
-    tempResults = WIPsubfnFitModel(tempModel);
+    tempResults = FitProcessModel(tempModel);
+
     % store all the parameter estimates in the JackKnife structure
     JackKnife.beta(:,:,i) = tempResults.beta;
     JackKnife.B(:,:,i) = tempResults.B;

FinalCode/OneVoxelProcessBootstrap.m

@@ -22,7 +22,7 @@
     JackKnife = JackKnifeFunction(Model,Results,FieldNames);
 
     % Calculate the BCaci values for each parameter
-    Results.BCaCI = WIPsubfnCreateBCaCI(Results,BootStrap,JackKnife,Model.Thresholds);
+    Results.BCaCI = CreateBCaCI(Results,BootStrap,JackKnife,Model.Thresholds);
 
 end
 

FinalCode/PrepareDataForProcess.m

@@ -71,8 +71,19 @@
     % This analysis is run on the host computer
     switch ModelType
         case 'bootstrap'
-            % Need to split the data and cycle over voxels
-            Results = CycleOverVoxelsProcessBootstrap(ModelInfo);
+            % Save the data
+            InDataPath = fullfile(DataFolder,'ModelInfo');
+            Str = ['save ' InDataPath ' ModelInfo '];
+            eval(Str);
+            % Process the data
+            Parameters = CycleOverVoxelsProcessBootstrap(ModelInfo);
+
+            % Save the results
+            ResultsPath = fullfile(OutFolder,'Results');
+            mkdir(ResultsPath)
+            Str = sprintf('save %s %s',fullfile(ResultsPath,sprintf('Bootstrap_count%04d_%dSamp',1,ModelInfo.Nboot)),'Parameters');
+            eval(Str)
+
         case 'permutation'
             % The permutation analysis requires calculation of the point
             % estimate and then calculation of all of the permutations. The

FinalCode/WriteOutParameterMaps.m

@@ -34,8 +34,8 @@ function WriteOutParameterMaps(Tag,Parameters,ModelInfo)
                 FileName = sprintf('%s%sX',FileName,ModelInfo.Names{InterVar(j)});
             end
             FileName = sprintf('%s_%s.nii',FileName(1:end-1),Tag);
-            I = zeros(ModelInfor.DataHeader.dim);
-            I(ModelInfo.Indices) = squeeze(temp(Nvar+2,i,:));
+            I = zeros(ModelInfo.DataHeader.dim);
+            I(ModelInfo.Indices) = squeeze(temp(ModelInfo.Nvar+2,i,:));
             % Create the header for this image
             Vo = ModelInfo.DataHeader;
             Vo.fname = fullfile(ModelInfo.ResultsPath,FileName);