aguirreg

%
%   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