aguirreg

%
%   Generate Figure 2C (ROI analysis - bar graphs) for FaceSpace paper
%   _______________________________________________
%   by Marcelo G Mattar (10/21/2013)
%   mattar@sas.upenn.edu

% Parameters
datasets = {'MDS_RGI'}; % Choose between 'MDS_RGI' and 'MDS_Multi'
hemispheres_to_include = {'rh'}; % Choose {'rh'}, {'lh'}, or {'rh','lh'}
masksize = 800; % Number of vertices to form ROI (within fusiform mask)
outputDir = './Figures/'; % Directory where output images will be saved
saveBool = false;

% Data location
dataDir = '~/Data/FaceSpace/';
% Subjects to include
subjectList_RGI = {...
    % [list of subject ID strings redacted]
    };
subjectList_Multi = {...
    % [list of subject ID strings redacted]
    };

% Include auxiliar functions
addpath('./auxFunctions');

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
    prototypeOnly = zeros(1,numSubjects);
    adaptOnly = zeros(1,numSubjects);
    Drift_000_100_prototype = zeros(1,numSubjects);
    Drift_000_100_adapt = zeros(1,numSubjects);
    
    % Load individual subjects' data
    Drift_000_lh = squeeze(load_mgh([dataDir 'group_glms/' dataset '/Drift_000/lh.Drift_000.fwhm10.mgh']));
    Drift_000_rh = squeeze(load_mgh([dataDir 'group_glms/' dataset '/Drift_000/rh.Drift_000.fwhm10.mgh']));
    Drift_100_lh = squeeze(load_mgh([dataDir 'group_glms/' dataset '/Drift_100/lh.Drift_100.fwhm10.mgh']));
    Drift_100_rh = squeeze(load_mgh([dataDir 'group_glms/' dataset '/Drift_100/rh.Drift_100.fwhm10.mgh']));
    Drift_000_100_prototype_lh = squeeze(load_mgh([dataDir 'group_glms/' dataset '/DriftPrototype/lh.DriftPrototype.fwhm10.mgh']));
    Drift_000_100_prototype_rh = squeeze(load_mgh([dataDir 'group_glms/' dataset '/DriftPrototype/rh.DriftPrototype.fwhm10.mgh']));
    Drift_000_100_adapt_lh = squeeze(load_mgh([dataDir 'group_glms/' dataset '/DriftAdapt/lh.DriftAdapt.fwhm10.mgh']));
    Drift_000_100_adapt_rh = squeeze(load_mgh([dataDir 'group_glms/' dataset '/DriftAdapt/rh.DriftAdapt.fwhm10.mgh']));
    
    % Load group data
    Drift_000_betas_rh = load_mgh(['~/Data/FaceSpace/group_glms/' dataset '/Drift_000/rh.Drift_000.glmdir/beta.mgh']);
    Drift_000_betas_lh = load_mgh(['~/Data/FaceSpace/group_glms/' dataset '/Drift_000/lh.Drift_000.glmdir/beta.mgh']);
    Drift_100_betas_rh = load_mgh(['~/Data/FaceSpace/group_glms/' dataset '/Drift_100/rh.Drift_100.glmdir/beta.mgh']);
    Drift_100_betas_lh = load_mgh(['~/Data/FaceSpace/group_glms/' dataset '/Drift_100/lh.Drift_100.glmdir/beta.mgh']);
    Drift_000_100_adapt_betas_rh = load_mgh(['~/Data/FaceSpace/group_glms/' dataset '/DriftAdapt/rh.DriftAdapt.glmdir/beta.mgh']);
    Drift_000_100_adapt_betas_lh = load_mgh(['~/Data/FaceSpace/group_glms/' dataset '/DriftAdapt/lh.DriftAdapt.glmdir/beta.mgh']);
    Drift_000_100_prototype_betas_rh = load_mgh(['~/Data/FaceSpace/group_glms/' dataset '/DriftPrototype/rh.DriftPrototype.glmdir/beta.mgh']);
    Drift_000_100_prototype_betas_lh = load_mgh(['~/Data/FaceSpace/group_glms/' dataset '/DriftPrototype/lh.DriftPrototype.glmdir/beta.mgh']);
    
    for s=1:numSubjects
        % Calculate the average effect within ROI
        for h=1:length(hemispheres_to_include)
            eval(['prototypeOnly(s) = prototypeOnly(s) + mean(Drift_000_' hemispheres_to_include{h} '(mask_' hemispheres_to_include{h} ',s));']);
            eval(['adaptOnly(s) = adaptOnly(s) + mean(Drift_100_' hemispheres_to_include{h} '(mask_' hemispheres_to_include{h} ',s));']);
            eval(['Drift_000_100_prototype(s) = Drift_000_100_prototype(s) + mean(Drift_000_100_prototype_' hemispheres_to_include{h} '(mask_' hemispheres_to_include{h} ',s));']);
            eval(['Drift_000_100_adapt(s) = Drift_000_100_adapt(s) + mean(Drift_000_100_adapt_' hemispheres_to_include{h} '(mask_' hemispheres_to_include{h} ',s));']);
        end
        % Now, take the mean across hemispheres (if looking at both
        % together)
        prototypeOnly(s) = prototypeOnly(s) / length(hemispheres_to_include);
        adaptOnly(s) = adaptOnly(s) / length(hemispheres_to_include);
        Drift_000_100_prototype(s) = Drift_000_100_prototype(s) / length(hemispheres_to_include);
        Drift_000_100_adapt(s) = Drift_000_100_adapt(s) / length(hemispheres_to_include);
    end
    
    % Calculate means and standard errors of the mean
    mean_prototypeOnly = mean(prototypeOnly);
    mean_adaptOnly = mean(adaptOnly);
    mean_prototype = mean(Drift_000_100_prototype);
    mean_adapt = mean(Drift_000_100_adapt);
    sem_prototypeOnly = std(prototypeOnly)/sqrt(numSubjects);
    sem_adaptOnly = std(adaptOnly)/sqrt(numSubjects);
    sem_prototype = std(Drift_000_100_prototype)/sqrt(numSubjects);
    sem_adapt = std(Drift_000_100_adapt)/sqrt(numSubjects);
    
    % Calculate statistics
    [~,p_adaptOnly,ci_adaptOnly,stats_adaptOnly] = ttest(adaptOnly,0,'tail','both');
    [~,p_prototypeOnly,ci_prototypeOnly,stats_prototypeOnly] = ttest(prototypeOnly,0,'tail','both');
    % Display summary statistics
    display('*****************************************************************************');
    display('Discrete Model statistics:');
    display('     Adapt-Only:');
    display(['M = ' num2str(mean_adaptOnly) ', SD: = ' num2str(std(adaptOnly)) ', SEM = ' num2str(sem_adaptOnly) '.']);
    display(['One-sample t(' num2str(stats_adaptOnly.df) ') = ' num2str(stats_adaptOnly.tstat) ', p = ' num2str(p_adaptOnly) ' (two-tailed), 95% CI [' num2str(ci_adaptOnly(1)) ', ' num2str(ci_adaptOnly(2)) '].']);
    display('     Prototype-Only:');
    display(['M = ' num2str(mean_prototypeOnly) ', SD: = ' num2str(std(prototypeOnly)) ', SEM = ' num2str(sem_prototypeOnly) '.']);
    display(['one-sample t(' num2str(stats_prototypeOnly.df) ') = ' num2str(stats_prototypeOnly.tstat) ', p = ' num2str(p_prototypeOnly) ' (two-tailed), 95% CI [' num2str(ci_prototypeOnly(1)) ', ' num2str(ci_prototypeOnly(2)) '].']);
    display('*****************************************************************************');
    
    % Plot the results
    figure(1);
    h = barwitherr([sem_prototypeOnly sem_adaptOnly], [mean_prototypeOnly mean_adaptOnly]);% Plot with errorbars
    ylim([0 0.4*10^(-3)])
    set(gca,'XTickLabel',{'Prototype only','Adapt only'})
    ylabel('Beta values')
    title('Discrete Models');
    if saveBool; print('-depsc',[outputDir 'Figure2_discretemodels_' dataset '.eps']); end;
    if saveBool; print('-dpng',[outputDir 'Figure2_discretemodels_' dataset '.png']); end;
    
    figure(2);
    h = barwitherr([sem_prototype sem_adapt], [mean_prototype mean_adapt]);% Plot with errorbars
    ylim([0 0.4*10^(-3)])
    set(gca,'XTickLabel',{'Adapt&ID Prototype','Adapt&ID Adapt'})
    ylabel('Beta values')
    title('Compound Models');
    if saveBool; print('-depsc',[outputDir 'Figure2_compoundmodels_' dataset '.eps']); end;
    if saveBool; print('-dpng',[outputDir 'Figure2_compoundmodels_' dataset '.png']); end;
end