%
% Compare the single mechanism (Drift) model with a double mechanism
% (adapt + prototype) model by using a leave-one-out approach.
% The elasticity coefficient for each subject is calculated as the best
% EC from all other subjects' data.
% _______________________________________________
% by Marcelo G Mattar (10/21/2013)
% mattar@sas.upenn.edu
% Parameters
datasets = {'MDS_Multi'}; % Choose between 'MDS_RGI' and 'MDS_Multi'
hemispheres_to_include = {'rh'}; % Choose {'rh'}, {'lh'}, or {'rh','lh'}
% Directories
baseDir = '/jet/aguirre/FaceSpace/';
dataDir = '~/Data/FaceSpace/';
FSpath = '/Applications/freesurfer';
SUBJECTS_DIR = '/Volumes/Blob-iron/freesurfer/subjects';
% Additional parameters
FWHM = 10;
masksize = 800; % Number of vertices to form ROI (within fusiform mask)
SOA=1500;
subjectList_RGI = {...
% [list of subject ID strings redacted]
};
subjectList_Multi = {...
% [list of subject ID strings redacted]
};
% Include auxiliar functions
addpath('./auxFunctions');
filename_endings = {'000' '005' '010' '015' '020' '025' '030' '035' '040' '045' '050' '055' '060' '065' '070' '075' '080' '085' '090' '095' '100'};
coefficient_values = 0:0.05:1;
numCoef = length(filename_endings);
for d = 1:length(datasets)
% Determine which dataset we'll work on
dataset = datasets{d};
switch dataset
case 'MDS_RGI'
subjectList = subjectList_RGI;
dataset_lowercase = '-mds_rgi';
case 'MDS_Multi'
subjectList = subjectList_Multi;
dataset_lowercase = '-mds_multi';
end
numSubjects = length(subjectList);
% Load masks
mask_lh = load_mgh([dataDir 'ROIs/binaryMask_' dataset '_' num2str(masksize) '_lh.mgh']);
mask_rh = load_mgh([dataDir 'ROIs/binaryMask_' dataset '_' num2str(masksize) '_rh.mgh']);
mask_lh = logical(mask_lh);
mask_rh = logical(mask_rh);
% Initialize variables
DriftBetas_lh = zeros(length(mask_lh),numSubjects,numCoef);
DriftBetas_rh = zeros(length(mask_rh),numSubjects,numCoef);
averageBetasROI = zeros(numSubjects,length(filename_endings));
% Load data for all subjects
for k=1:numCoef
DriftBetas_lh(:,:,k) = squeeze(load_mgh([dataDir 'group_glms/' dataset '/Drift_' filename_endings{k} '/lh.Drift_' filename_endings{k} '.fwhm10.mgh']));
DriftBetas_rh(:,:,k) = squeeze(load_mgh([dataDir 'group_glms/' dataset '/Drift_' filename_endings{k} '/rh.Drift_' filename_endings{k} '.fwhm10.mgh']));
% Calculate the average effect within ROI
for s=1:numSubjects
for h=1:length(hemispheres_to_include)
eval(['averageBetasROI(s,k) = averageBetasROI(s,k) + mean(DriftBetas_' hemispheres_to_include{h} '(mask_' hemispheres_to_include{h} ',s,k));']);
end
% Now, take the mean across hemispheres (if looking at both together)
averageBetasROI(s,k) = averageBetasROI(s,k) / length(hemispheres_to_include);
end
end
meanBetas = zeros(numSubjects,21);
bestFit = zeros(numSubjects,1);
%% Leave-one-out approach
for leftout = 1:numSubjects
% Calculate the best elasticity coefficient in the reduced set of
% subjects
meanBetas(leftout,:) = mean(averageBetasROI(1:numSubjects~=leftout,:),1);
[~,bestFit(leftout)] = max(meanBetas(leftout,:));
end
% Pre-allocate variables to output the results of the leave-one-out
betaValues1 = zeros(numSubjects,2);
betaValues2 = zeros(numSubjects,3);
stats1 = zeros(numSubjects,4);
stats2 = zeros(numSubjects,4);
error_series1 = zeros(numSubjects,1);
error_series2 = zeros(numSubjects,1);
adjR2stats1 = zeros(numSubjects,1);
adjR2stats2 = zeros(numSubjects,1);
% Finaly, compare both models for each subject
for s = 1:numSubjects
% First, load the time-series from the ROI
[timeseries_rh, timeseries_lh] = extract_timeseries(dataset, subjectList{s});
% Extract timeseries from selected hemispheres
this_timeseries = 0;
for h=1:length(hemispheres_to_include)
eval(['this_timeseries = this_timeseries + timeseries_' hemispheres_to_include{h} ';']);
end
% Now, take the mean across hemispheres (if looking at both together)
this_timeseries = this_timeseries / length(hemispheres_to_include);
% Extract the best elasticity coefficient from the remaining
% subject's data
this_elasticity = bestFit(s);
% Then, generate the necessary covariates for this subject
drift_covariate = driftCov(dataset, subjectList{s}, str2double(filename_endings{this_elasticity}));
adapt_covariate = driftCov(dataset, subjectList{s}, 100);
prototype_covariate = driftCov(dataset, subjectList{s}, 0);
% Convolve each covariate with the HRF
drift_covariate = hrfconv(drift_covariate,SOA);
adapt_covariate = hrfconv(adapt_covariate,SOA);
prototype_covariate = hrfconv(prototype_covariate,SOA);
% Downsample all covariates to match the length of the timeseries
drift_covariate = downsample(drift_covariate,(length(drift_covariate))/(length(this_timeseries)));
adapt_covariate = downsample(adapt_covariate,(length(adapt_covariate))/(length(this_timeseries)));
prototype_covariate = downsample(prototype_covariate,(length(prototype_covariate))/(length(this_timeseries)));
% Now, run a simple/multiple regression with these covariates
[betaValues1(s,:),~,~,~,stats1(s,:)] = regress(this_timeseries,[ones(size(drift_covariate)) drift_covariate]);
[betaValues2(s,:),~,~,~,stats2(s,:)] = regress(this_timeseries,[ones(size(adapt_covariate)) adapt_covariate prototype_covariate]);
% ------ NOTE -------
% stats is a 1-by-4 vector stats that contains, in order
% R2 statistic
% F statistic
% p value
% an estimate of the error variance.
%{
% Calculate residual series
error_series1(s) = sum((this_timeseries - sum(repmat(betaValues1(s,:),length(this_timeseries),1) .* [ones(size(drift_covariate)) drift_covariate],2)).^2);
error_series2(s) = sum((this_timeseries - sum(repmat(betaValues2(s,:),length(this_timeseries),1) .* [ones(size(adapt_covariate)) adapt_covariate prototype_covariate],2)).^2);
% Calculate all regression statistics (to extract the adjusted r-squared)
regstats1 = regstats(this_timeseries,drift_covariate,'linear',{'beta' 'fstat' 'rsquare' 'adjrsquare'});
regstats2 = regstats(this_timeseries,[adapt_covariate prototype_covariate],'linear',{'beta' 'fstat' 'rsquare' 'adjrsquare'});
adjR2stats1(s) = regstats1.adjrsquare;
adjR2stats2(s) = regstats2.adjrsquare;
%}
end
% Calculate statistics
drift_VS_compound_model_numberOfTimesWin = sum(stats1(:,3)<stats2(:,3));
binomialTest_pvalue = myBinomTest(drift_VS_compound_model_numberOfTimesWin,numSubjects,0.5,'Greater');
% Display summary statistics
display('*****************************************************************************');
display(['Dataset: ' dataset ' - Hemisphere: ' hemispheres_to_include{1}]);
display(sprintf('The Drift model has a smaller p-value than the AdaptPlusID model in %d out of %d subjects.', drift_VS_compound_model_numberOfTimesWin, numSubjects));
display(['One-sided binomial test, B(' num2str(numSubjects) ',0.5), p = ' num2str(binomialTest_pvalue) '.']);
display('*****************************************************************************');
end