include analysis scripts and data

.gitignore

@@ -0,0 +1,15 @@
+# do not consider anything
+
+*
+
+# except for files and directories
+
+!*.*
+!*/
+
+# with the following exclusions
+
+*.m~
+*.DS_Store
+*._*
+*~*

.gitmodules

@@ -0,0 +1,12 @@
+[submodule "tools/BrewerMap"]
+	path = tools/BrewerMap
+	url = git@github.com:DrosteEffect/BrewerMap.git
+[submodule "tools/RainCloudPlots"]
+	path = tools/RainCloudPlots
+	url = git@github.com:jkosciessa/RainCloudPlots.git
+[submodule "tools/fieldtrip"]
+	path = tools/fieldtrip
+	url = https://github.com/fieldtrip/fieldtrip.git
+[submodule "tools/shadedErrorBar"]
+	path = tools/shadedErrorBar
+	url = https://github.com/raacampbell/shadedErrorBar.git

code/a_eegmp_calculate_erp_tfr.m

@@ -0,0 +1,86 @@
+function a_eegmp_calculate_erp_tfr(rootpath, id)
+
+%% paths
+if ismac
+    currentFile = mfilename('fullpath');
+    [pathstr,~,~] = fileparts(currentFile);
+    cd(fullfile(pathstr,'..'))
+    rootpath = pwd;
+    id = 'sub-011';
+end
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'erp'); mkdir(pn.data_erp);
+pn.data_erf = fullfile(rootpath, 'data', 'erf'); mkdir(pn.data_erf);
+pn.tools = fullfile(rootpath, '..', 'eegmp_preproc', 'tools');
+    addpath(fullfile(pn.tools, 'fieldtrip')); ft_defaults
+
+%% load data_eeg
+
+% load preprocessed eeg data_eeg
+load(fullfile(pn.data_eeg, [id,'_task-xxxx_eeg_art.mat']), 'data_eeg', 'events');
+
+%% further preprocessing
+
+% apply notch filter
+
+cfg = [];
+cfg.dftfilter = 'yes';
+data_eeg = ft_preprocessing(cfg, data_eeg);
+
+%% CSD transform
+    
+csd_cfg = [];
+csd_cfg.elecfile = 'standard_1005.elc';
+csd_cfg.method = 'spline';
+data_eeg = ft_scalpcurrentdensity(csd_cfg, data_eeg);
+
+%% time-frequency transform using wavelets
+
+freq_cfg = [];
+freq_cfg.channel     = 'all';
+freq_cfg.method      = 'wavelet';
+freq_cfg.width       = 5;
+freq_cfg.keeptrials  = 'yes';
+freq_cfg.output      = 'pow';
+freq_cfg.foi         = 1:1:40;
+freq_cfg.toi         = -1:0.05:2;
+tfr = ft_freqanalysis(freq_cfg, data_eeg);
+
+% single-trial log10-transform
+
+tfr.powspctrm = log10(tfr.powspctrm);
+
+%% single-trial baseline-correction for ERPs
+
+time = data_eeg.time{1};
+
+data_eeg_bl = data_eeg;
+
+for indTrial = 1:numel(data_eeg_bl.trial)
+    curbl = squeeze(nanmean(data_eeg_bl.trial{indTrial}(:,time>-.2 & time<0),2));
+    data_eeg_bl.trial{indTrial} = data_eeg_bl.trial{indTrial}-repmat(curbl,1,numel(time));
+end
+
+%% split tfr and time series data by condition, average across trials
+
+parameter = {'scene_category'; 'old'; 'behavior'; 'subsequent_memory'};
+for ind_param = 1:numel(parameter)
+    conds.(parameter{ind_param}) = unique(events.(parameter{ind_param}));
+    for ind_cond = 1:numel(conds.(parameter{ind_param}))
+        cfg = [];
+        cfg.trials = ismember(events.(parameter{ind_param}), ...
+            conds.(parameter{ind_param}){ind_cond});
+        tfravg.(parameter{ind_param}){ind_cond} = ft_freqdescriptives(cfg, tfr);
+        erp.(parameter{ind_param}){ind_cond} = ft_timelockanalysis(cfg, data_eeg);
+        erp_bl.(parameter{ind_param}){ind_cond} = ft_timelockanalysis(cfg, data_eeg_bl);
+    end
+end
+
+erf = tfravg; clear tfravg; 
+
+%% save data
+
+save(fullfile(pn.data_erp, [id, '_erp.mat']), 'erp', 'conds');
+save(fullfile(pn.data_erp, [id, '_erp_bl.mat']), 'erp_bl', 'conds');
+save(fullfile(pn.data_erf, [id, '_erf.mat']), 'erf', 'conds');

code/a_eegmp_calculate_erp_tfr_START.sh

@@ -0,0 +1,27 @@
+#!/bin/bash
+
+#!/bin/bash
+
+fun_name="a_eegmp_calculate_erp_tfr"
+job_name="eegmp_erptfr"
+
+rootpath="$(pwd)/.."
+rootpath=$(builtin cd $rootpath; pwd)
+
+mkdir ${rootpath}/log
+
+for subj in $(seq 1 33); do
+	cursub=$(printf "sub-%03d\n" ${subj})
+	echo ${cursub}
+	# if not been performed yet
+	#if [ ! -f ${rootpath}/data/erf/sub-$(printf "%03d" ${subj})*_erf.mat ]; then
+		sbatch \
+  			--job-name ${job_name}_${cursub} \
+  			--cpus-per-task 4 \
+  			--mem 8gb \
+  			--time 00:30:00 \
+  			--output ${rootpath}/log/${job_name}_${cursub}.out \
+  			--workdir . \
+  			--wrap=". /etc/profile ; module load matlab/R2016b; matlab -nodisplay -r \"${fun_name}('${rootpath}','${cursub}')\""
+  	#fi
+done

code/b1_taskPLS_scenecat_N1.m

@@ -0,0 +1,109 @@
+% Set up EEG PLS for a task PLS using spectral power
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'erp');
+pn.data_erf = fullfile(rootpath, 'data', 'erf');
+pn.data = fullfile(rootpath, 'data', 'stats'); mkdir(pn.data);
+pn.tools = fullfile(rootpath, 'tools');
+    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'shadedErrorBar'));
+
+%% add seed for reproducibility
+
+rng(0, 'twister');
+    
+%% load event info
+
+load(fullfile(pn.data_eeg, ['sub-001_task-xxxx_eeg_art.mat']), 'events');
+
+parameter = {'scene_category'; 'old'; 'behavior'; 'subsequent_memory'};
+for ind_param = 1:numel(parameter)
+    conds.(parameter{ind_param}) = unique(events.(parameter{ind_param}));
+end
+
+%% load erp
+
+for ind_id = 1:33
+    id = sprintf('sub-%03d', ind_id);
+    load(fullfile(pn.data_erp, [id,'_erp_bl.mat']));
+    for ind_option = 1:numel(conds.scene_category)
+        if ind_id == 1
+             erpgroup.scene_category.(conds.scene_category{ind_option}) = erp_bl.scene_category{ind_option};
+             erpgroup.scene_category.(conds.scene_category{ind_option}) = ...
+                 rmfield(erpgroup.scene_category.(conds.scene_category{ind_option}), {'avg', 'var', 'dof'});
+             erpgroup.scene_category.(conds.scene_category{ind_option}).dimord = 'sub_chan_time';
+        end
+        % only include time points from 0 to 300 ms
+        time = erp_bl.scene_category{ind_option}.time;
+        toi = time >0 & time <0.3;        
+        erpgroup.scene_category.(conds.scene_category{ind_option}).avg(ind_id,:,:) = erp_bl.scene_category{ind_option}.avg(:,toi);
+    end
+end
+
+time = erpgroup.scene_category.manmade.time(toi);
+elec = erpgroup.scene_category.manmade.elec;
+channels = erpgroup.scene_category.manmade.label;
+
+mergeddata = cat(4, erpgroup.scene_category.manmade.avg, ...
+    erpgroup.scene_category.natural.avg);
+
+num_chans = numel(channels);
+num_freqs = 1;
+num_time = numel(time);
+
+num_subj_lst = [33];
+num_cond = 2;
+num_grp = 1;
+
+datamat_lst = cell(num_grp); lv_evt_list = [];
+indCount_cont = 1;
+indCount = 1;
+indGroup = 1;
+for indCond = 1:num_cond
+    for indID = 1:num_subj_lst(indGroup)
+        datamat_lst{indGroup}(indCount,:) = reshape(squeeze(mergeddata(indID,:,:,indCond)), [], 1);
+        lv_evt_list(indCount_cont) = indCond;
+        indCount = indCount+1;
+        indCount_cont = indCount_cont+1;
+    end
+end
+datamat_lst{indGroup}(isnan(datamat_lst{indGroup})) = 0;
+
+%% set PLS options and run PLS
+
+option = [];
+option.method = 1; % [1] | 2 | 3 | 4 | 5 | 6
+option.num_perm = 1000; %( single non-negative integer )
+option.num_split = 0; %( single non-negative integer )
+option.num_boot = 1000; % ( single non-negative integer )
+option.cormode = 0; % [0] | 2 | 4 | 6
+option.meancentering_type = 0;% [0] | 1 | 2 | 3
+option.boot_type = 'strat'; %['strat'] | 'nonstrat'
+
+result = pls_analysis(datamat_lst, num_subj_lst, num_cond, option);
+
+%% rearrange into fieldtrip structure
+
+lvdat = reshape(result.boot_result.compare_u(:,1), num_chans, num_freqs, num_time);
+%udat = reshape(result.u, num_chans, num_freqs, num_time);
+
+stat = [];
+stat.prob = lvdat;
+stat.dimord = 'chan_freq_time';
+stat.clusters = [];
+stat.clusters.prob = result.perm_result.sprob; % check for significance of LV
+stat.mask = lvdat > 3 | lvdat < -3;
+stat.cfg = option;
+stat.time = time;
+
+save(fullfile(pn.data, 'b01_taskpls_erp.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')

code/b1_taskPLS_scenecat_N1_avg.m

@@ -0,0 +1,111 @@
+% Set up EEG PLS for a task PLS using spectral power
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'erp');
+pn.data_erf = fullfile(rootpath, 'data', 'erf');
+pn.data = fullfile(rootpath, 'data', 'stats'); mkdir(pn.data);
+pn.tools = fullfile(rootpath, 'tools');
+    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'shadedErrorBar'));
+
+%% add seed for reproducibility
+
+rng(0, 'twister');
+    
+%% load event info
+
+load(fullfile(pn.data_eeg, ['sub-001_task-xxxx_eeg_art.mat']), 'events');
+
+parameter = {'scene_category'; 'old'; 'behavior'; 'subsequent_memory'};
+for ind_param = 1:numel(parameter)
+    conds.(parameter{ind_param}) = unique(events.(parameter{ind_param}));
+end
+
+%% load erp
+
+for ind_id = 1:33
+    id = sprintf('sub-%03d', ind_id);
+    load(fullfile(pn.data_erp, [id,'_erp.mat']));
+    for ind_option = 1:numel(conds.scene_category)
+        if ind_id == 1
+             erpgroup.scene_category.(conds.scene_category{ind_option}) = erp.scene_category{ind_option};
+             erpgroup.scene_category.(conds.scene_category{ind_option}) = ...
+                 rmfield(erpgroup.scene_category.(conds.scene_category{ind_option}), {'avg', 'var', 'dof'});
+             erpgroup.scene_category.(conds.scene_category{ind_option}).dimord = 'sub_chan_time';
+        end
+        % average in N1 time window
+        time = erpgroup.scene_category.(conds.scene_category{ind_option}).time;
+        toi = time >0.08 & time <0.12;
+        erpgroup.scene_category.(conds.scene_category{ind_option}).avg(ind_id,:,:) = nanmean(erp.scene_category{ind_option}.avg(:,toi),2);
+    end
+end
+
+% average in N1 time window
+
+time = 0.1;
+elec = erpgroup.scene_category.manmade.elec;
+channels = erpgroup.scene_category.manmade.label;
+
+mergeddata = cat(4, erpgroup.scene_category.manmade.avg, ...
+    erpgroup.scene_category.natural.avg);
+
+num_chans = numel(channels);
+num_freqs = 1;
+num_time = numel(time);
+
+num_subj_lst = [33];
+num_cond = 2;
+num_grp = 1;
+
+datamat_lst = cell(num_grp); lv_evt_list = [];
+indCount_cont = 1;
+indCount = 1;
+indGroup = 1;
+for indCond = 1:num_cond
+    for indID = 1:num_subj_lst(indGroup)
+        datamat_lst{indGroup}(indCount,:) = reshape(squeeze(mergeddata(indID,:,:,indCond)), [], 1);
+        lv_evt_list(indCount_cont) = indCond;
+        indCount = indCount+1;
+        indCount_cont = indCount_cont+1;
+    end
+end
+datamat_lst{indGroup}(isnan(datamat_lst{indGroup})) = 0;
+
+%% set PLS options and run PLS
+
+option = [];
+option.method = 1; % [1] | 2 | 3 | 4 | 5 | 6
+option.num_perm = 1000; %( single non-negative integer )
+option.num_split = 0; %( single non-negative integer )
+option.num_boot = 1000; % ( single non-negative integer )
+option.cormode = 0; % [0] | 2 | 4 | 6
+option.meancentering_type = 0;% [0] | 1 | 2 | 3
+option.boot_type = 'strat'; %['strat'] | 'nonstrat'
+
+result = pls_analysis(datamat_lst, num_subj_lst, num_cond, option);
+
+%% rearrange into fieldtrip structure
+
+lvdat = reshape(result.boot_result.compare_u(:,1), num_chans, num_freqs, num_time);
+%udat = reshape(result.u, num_chans, num_freqs, num_time);
+
+stat = [];
+stat.prob = lvdat;
+stat.dimord = 'chan_freq_time';
+stat.clusters = [];
+stat.clusters.prob = result.perm_result.sprob; % check for significance of LV
+stat.mask = lvdat > 3 | lvdat < -3;
+stat.cfg = option;
+stat.time = time;
+
+save(fullfile(pn.data, 'b01_taskpls_erp_avg.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')

code/b1_taskPLS_scenecat_N1_avg_plotLV1.m

@@ -0,0 +1,153 @@
+% Create an overview plot featuring the results of the multivariate PLS
+% comparing spectral changes during the stimulus period under load
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data         = fullfile(rootpath, 'data', 'stats');
+pn.figures      = fullfile(rootpath, 'figures');
+pn.tools        = fullfile(rootpath, 'tools');
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'fieldtrip')); ft_defaults;
+    addpath(genpath(fullfile(pn.tools, 'RainCloudPlots')));
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'winsor'));
+
+% set custom colormap
+cBrew = brewermap(500,'RdBu');
+cBrew = flipud(cBrew);
+colormap(cBrew)
+
+load(fullfile(pn.data, 'b01_taskpls_erp_avg.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')
+load(fullfile(rootpath, 'data','erp', ['sub-001_erp.mat']));
+elec = erp.scene_category{1}.elec;
+
+result.perm_result.sprob
+
+indLV = 1;
+
+lvdat = reshape(result.boot_result.compare_u(:,indLV), num_chans, num_freqs, num_time);
+stat.prob = lvdat;
+stat.mask = lvdat > 3 | lvdat < -3;
+
+% maskNaN = double(stat.mask);
+% maskNaN(maskNaN==0) = NaN;
+
+%% plot multivariate topographies
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+cfg.zlim = [-6 6]; 
+cfg.figure = h;
+plotData.powspctrm = squeeze(stat.mask.*stat.prob); 
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Mean BSR');
+figureName = ['b01_lv1'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');
+
+%% plot using raincloud plot
+
+groupsizes=result.num_subj_lst;
+conditions=lv_evt_list;
+conds = {'manmade'; 'natural'};
+condData = []; uData = [];
+for indGroup = 1
+    if indGroup == 1
+        relevantEntries = 1:groupsizes(1)*numel(conds);
+    elseif indGroup == 2
+        relevantEntries = groupsizes(1)*numel(conds)+1:...
+             groupsizes(1)*numel(conds)+groupsizes(2)*numel(conds);
+    end
+    for indCond = 1:numel(conds)
+        targetEntries = relevantEntries(conditions(relevantEntries)==indCond);        
+        condData{indGroup}(indCond,:) = result.vsc(targetEntries,indLV);
+        uData{indGroup}(indCond,:) = result.usc(targetEntries,indLV);
+    end
+end
+
+%% plot RainCloudPlot (within-subject centered)
+
+cBrew(1,:) = 2.*[.3 .1 .1];
+cBrew(2,:) = [.6 .6 .6];
+
+idx_outlier = cell(1); idx_standard = cell(1);
+for indGroup = 1
+    dataToPlot = uData{indGroup}';
+    % define outlier as lin. modulation that is more than three scaled median absolute deviations (MAD) away from the median
+    X = [1 1; 1 2]; b=X\dataToPlot'; IndividualSlopes = b(2,:);
+    outliers = isoutlier(IndividualSlopes, 'median');
+    idx_outlier{indGroup} = find(outliers);
+    idx_standard{indGroup} = find(outliers==0);
+end
+
+h = figure('units','centimeter','position',[0 0 25 10]);
+for indGroup = 1
+    dataToPlot = uData{indGroup}';
+    % read into cell array of the appropriate dimensions
+    data = []; data_ws = [];
+    for i = 1:2
+        for j = 1:1
+            data{i, j} = dataToPlot(:,i);
+            % individually demean for within-subject visualization
+            data_ws{i, j} = dataToPlot(:,i)-...
+                nanmean(dataToPlot(:,:),2)+...
+                repmat(nanmean(nanmean(dataToPlot(:,:),2),1),size(dataToPlot(:,:),1),1);
+            data_nooutlier{i, j} = data{i, j};
+            data_nooutlier{i, j}(idx_outlier{indGroup}) = [];
+            data_ws_nooutlier{i, j} = data_ws{i, j};
+            data_ws_nooutlier{i, j}(idx_outlier{indGroup}) = [];
+            % sort outliers to back in original data for improved plot overlap
+            data_ws{i, j} = [data_ws{i, j}(idx_standard{indGroup}); data_ws{i, j}(idx_outlier{indGroup})];
+        end
+    end
+
+    % IMPORTANT: plot individually centered estimates, stats on uncentered estimates!
+
+    subplot(1,2,indGroup);
+    set(gcf,'renderer','Painters')
+        cla;
+        cl = cBrew(indGroup,:);
+        [~, dist] = rm_raincloud_fixedSpacing(data_ws, [.8 .8 .8],1);
+        h_rc = rm_raincloud_fixedSpacing(data_ws_nooutlier, cl,1,[],[],[],dist);
+        view([90 -90]);
+        axis ij
+    box(gca,'off')
+    set(gca, 'YTickLabels', {'natural'; 'manmade'});
+    yticks = get(gca, 'ytick'); ylim([yticks(1)-(yticks(2)-yticks(1))./2, yticks(2)+(yticks(2)-yticks(1))./1.5]);
+
+    minmax = [min(min(cat(2,data_ws{:}))), max(max(cat(2,data_ws{:})))];
+    xlim(minmax+[-0.2*diff(minmax), 0.2*diff(minmax)])
+    ylabel('scene category'); xlabel({'Brainscore'; '[Individually centered]'})
+
+    % test linear effect
+    curData = [data_nooutlier{1, 1}, data_nooutlier{2, 1}];
+    X = [1 1; 1 2]; b=X\curData'; IndividualSlopes = b(2,:);
+    [~, p, ci, stats] = ttest(IndividualSlopes);
+    title(['M:', num2str(round(mean(IndividualSlopes),3)), '; p=', num2str(round(p,3))])
+end
+figureName = ['b01_rcp'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');
+
+
+%% get channels to visualize
+
+positive = find(squeeze(stat.mask.*stat.prob)>0);
+negative = find(squeeze(stat.mask.*stat.prob)<0);

code/b1_taskPLS_scenecat_N1_plotLV1.m

@@ -0,0 +1,181 @@
+% Create an overview plot featuring the results of the multivariate PLS
+% comparing spectral changes during the stimulus period under load
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data         = fullfile(rootpath, 'data', 'stats');
+pn.figures      = fullfile(rootpath, 'figures');
+pn.tools        = fullfile(rootpath, 'tools');
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'fieldtrip')); ft_defaults;
+    addpath(genpath(fullfile(pn.tools, 'RainCloudPlots')));
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'winsor'));
+
+% set custom colormap
+cBrew = brewermap(500,'RdBu');
+cBrew = flipud(cBrew);
+colormap(cBrew)
+
+load(fullfile(pn.data, 'b01_taskpls_erp.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')
+load(fullfile(rootpath, 'data','erp', ['sub-001_erp.mat']));
+elec = erp.scene_category{1}.elec;
+
+result.perm_result.sprob
+
+indLV = 1;
+
+lvdat = reshape(result.boot_result.compare_u(:,indLV), num_chans, num_freqs, num_time);
+stat.prob = lvdat;
+stat.mask = lvdat > 3 | lvdat < -3;
+
+% maskNaN = double(stat.mask);
+% maskNaN(maskNaN==0) = NaN;
+
+%% invert solution
+
+% stat.mask = stat.mask;
+% stat.prob = stat.prob.*-1;
+% result.vsc = result.vsc.*-1;
+% result.usc = result.usc.*-1;
+
+h = figure('units','centimeter','position',[0 0 15 10]);
+set(gcf,'renderer','Painters')
+statsPlot = [];
+statsPlot = cat(1, statsPlot,squeeze(nanmax(abs(stat.prob(1:64,:,:)),[],1)));
+plot(stat.time,statsPlot, 'k')
+xlabel('Time [s from stim onset]'); ylabel('max abs BSR');
+title({'ERP changes'; ['p = ', num2str(round(result.perm_result.sprob(indLV),4))]})
+set(findall(gcf,'-property','FontSize'),'FontSize',18)
+figureName = ['b01_pls_traces'];
+saveas(h, fullfile(pn.figures, figureName), 'epsc');
+saveas(h, fullfile(pn.figures, figureName), 'png');
+
+
+h = figure('units','centimeter','position',[0 0 15 10]);
+set(gcf,'renderer','Painters')
+statsPlot = [];
+statsPlot = cat(1, statsPlot,squeeze(mean(stat.prob(1:64,:,:),1)));
+plot(stat.time,statsPlot, 'k')
+xlabel('Time [s from stim onset]'); ylabel('mean BSR');
+title({'ERP changes'; ['p = ', num2str(round(result.perm_result.sprob(indLV),4))]})
+set(findall(gcf,'-property','FontSize'),'FontSize',18)
+xlim([0 0.3])
+figureName = ['b01_pls_traces'];
+saveas(h, fullfile(pn.figures, figureName), 'epsc');
+saveas(h, fullfile(pn.figures, figureName), 'png');
+
+%% plot multivariate topographies
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+cfg.zlim = [-6 6]; 
+cfg.figure = h;
+plotData.powspctrm = squeeze(nanmean(stat.mask(:,:,stat.time >0.08 & stat.time <0.12).*...
+    stat.prob(:,:,stat.time >0.08 & stat.time <0.12),3)); 
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Mean BSR');
+figureName = ['b01_lv1'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');
+
+%% plot using raincloud plot
+
+groupsizes=result.num_subj_lst;
+conditions=lv_evt_list;
+conds = {'manmade'; 'natural'};
+condData = []; uData = [];
+for indGroup = 1
+    if indGroup == 1
+        relevantEntries = 1:groupsizes(1)*numel(conds);
+    elseif indGroup == 2
+        relevantEntries = groupsizes(1)*numel(conds)+1:...
+             groupsizes(1)*numel(conds)+groupsizes(2)*numel(conds);
+    end
+    for indCond = 1:numel(conds)
+        targetEntries = relevantEntries(conditions(relevantEntries)==indCond);        
+        condData{indGroup}(indCond,:) = result.vsc(targetEntries,indLV);
+        uData{indGroup}(indCond,:) = result.usc(targetEntries,indLV);
+    end
+end
+
+%% plot RainCloudPlot (within-subject centered)
+
+cBrew(1,:) = 2.*[.3 .1 .1];
+cBrew(2,:) = [.6 .6 .6];
+
+idx_outlier = cell(1); idx_standard = cell(1);
+for indGroup = 1
+    dataToPlot = uData{indGroup}';
+    % define outlier as lin. modulation that is more than three scaled median absolute deviations (MAD) away from the median
+    X = [1 1; 1 2]; b=X\dataToPlot'; IndividualSlopes = b(2,:);
+    outliers = isoutlier(IndividualSlopes, 'median');
+    idx_outlier{indGroup} = find(outliers);
+    idx_standard{indGroup} = find(outliers==0);
+end
+
+h = figure('units','centimeter','position',[0 0 25 10]);
+for indGroup = 1
+    dataToPlot = uData{indGroup}';
+    % read into cell array of the appropriate dimensions
+    data = []; data_ws = [];
+    for i = 1:2
+        for j = 1:1
+            data{i, j} = dataToPlot(:,i);
+            % individually demean for within-subject visualization
+            data_ws{i, j} = dataToPlot(:,i)-...
+                nanmean(dataToPlot(:,:),2)+...
+                repmat(nanmean(nanmean(dataToPlot(:,:),2),1),size(dataToPlot(:,:),1),1);
+            data_nooutlier{i, j} = data{i, j};
+            data_nooutlier{i, j}(idx_outlier{indGroup}) = [];
+            data_ws_nooutlier{i, j} = data_ws{i, j};
+            data_ws_nooutlier{i, j}(idx_outlier{indGroup}) = [];
+            % sort outliers to back in original data for improved plot overlap
+            data_ws{i, j} = [data_ws{i, j}(idx_standard{indGroup}); data_ws{i, j}(idx_outlier{indGroup})];
+        end
+    end
+
+    % IMPORTANT: plot individually centered estimates, stats on uncentered estimates!
+
+    subplot(1,2,indGroup);
+    set(gcf,'renderer','Painters')
+        cla;
+        cl = cBrew(indGroup,:);
+        [~, dist] = rm_raincloud_fixedSpacing(data_ws, [.8 .8 .8],1);
+        h_rc = rm_raincloud_fixedSpacing(data_ws_nooutlier, cl,1,[],[],[],dist);
+        view([90 -90]);
+        axis ij
+    box(gca,'off')
+    set(gca, 'YTickLabels', {'natural'; 'manmade'});
+    yticks = get(gca, 'ytick'); ylim([yticks(1)-(yticks(2)-yticks(1))./2, yticks(2)+(yticks(2)-yticks(1))./1.5]);
+
+    minmax = [min(min(cat(2,data_ws{:}))), max(max(cat(2,data_ws{:})))];
+    xlim(minmax+[-0.2*diff(minmax), 0.2*diff(minmax)])
+    ylabel('scene category'); xlabel({'Brainscore'; '[Individually centered]'})
+
+    % test linear effect
+    curData = [data_nooutlier{1, 1}, data_nooutlier{2, 1}];
+    X = [1 1; 1 2]; b=X\curData'; IndividualSlopes = b(2,:);
+    [~, p, ci, stats] = ttest(IndividualSlopes);
+    title(['M:', num2str(round(mean(IndividualSlopes),3)), '; p=', num2str(round(p,3))])
+end
+figureName = ['b01_rcp'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');

code/b2a_taskPLS_novelty_frontal_erp.m

@@ -0,0 +1,113 @@
+% Set up EEG PLS for a task PLS using spectral power
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'erp');
+pn.data_erf = fullfile(rootpath, 'data', 'erf');
+pn.data = fullfile(rootpath, 'data', 'stats'); mkdir(pn.data);
+pn.tools = fullfile(rootpath, 'tools');
+    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'shadedErrorBar'));
+
+%% add seed for reproducibility
+
+rng(0, 'twister');
+    
+%% load event info
+
+load(fullfile(pn.data_eeg, ['sub-001_task-xxxx_eeg_art.mat']), 'events');
+
+parameter = {'scene_category'; 'old'; 'behavior'; 'subsequent_memory'};
+for ind_param = 1:numel(parameter)
+    conds.(parameter{ind_param}) = unique(events.(parameter{ind_param}));
+end
+
+%% load erp
+
+for ind_id = 1:33
+    id = sprintf('sub-%03d', ind_id); disp(id)
+    load(fullfile(pn.data_erp, [id,'_erp_bl.mat']));
+    for ind_option = 1:numel(conds.old)
+        if ind_id == 1
+             erpgroup.old.(conds.old{ind_option}) = erp_bl.old{ind_option};
+             erpgroup.old.(conds.old{ind_option}) = ...
+                 rmfield(erpgroup.old.(conds.old{ind_option}), {'avg', 'var', 'dof'});
+             erpgroup.old.(conds.old{ind_option}).dimord = 'sub_chan_time';
+        end
+        % only include time points from 300 - 500 ms
+        time = erpgroup.old.(conds.old{ind_option}).time;
+        toi = time >0.3 & time <0.5;
+        erpgroup.old.(conds.old{ind_option}).avg(ind_id,:,:) = erp_bl.old{ind_option}.avg(:,toi);
+    end
+end
+
+time = erpgroup.old.old.time(toi);
+elec = erpgroup.old.old.elec;
+channels = erpgroup.old.old.label;
+
+mergeddata = cat(4, erpgroup.old.old.avg, ...
+    erpgroup.old.new.avg);
+
+% identify frontal channels, and set everything else to zero
+idx_chans = startsWith(channels, 'F') | startsWith(channels, 'A') | startsWith(channels, 'C') & ~startsWith(channels, 'CP');
+mergeddata(:, idx_chans==0,:,:) = 0;
+
+num_chans = numel(channels);
+num_freqs = 1;
+num_time = numel(time);
+
+num_subj_lst = [33];
+num_cond = 2;
+num_grp = 1;
+
+datamat_lst = cell(num_grp); lv_evt_list = [];
+indCount_cont = 1;
+indCount = 1;
+indGroup = 1;
+for indCond = 1:num_cond
+    for indID = 1:num_subj_lst(indGroup)
+        datamat_lst{indGroup}(indCount,:) = reshape(squeeze(mergeddata(indID,:,:,indCond)), [], 1);
+        lv_evt_list(indCount_cont) = indCond;
+        indCount = indCount+1;
+        indCount_cont = indCount_cont+1;
+    end
+end
+datamat_lst{indGroup}(isnan(datamat_lst{indGroup})) = 0;
+
+%% set PLS options and run PLS
+
+option = [];
+option.method = 1; % [1] | 2 | 3 | 4 | 5 | 6
+option.num_perm = 1000; %( single non-negative integer )
+option.num_split = 0; %( single non-negative integer )
+option.num_boot = 1000; % ( single non-negative integer )
+option.cormode = 0; % [0] | 2 | 4 | 6
+option.meancentering_type = 0;% [0] | 1 | 2 | 3
+option.boot_type = 'strat'; %['strat'] | 'nonstrat'
+
+result = pls_analysis(datamat_lst, num_subj_lst, num_cond, option);
+
+%% rearrange into fieldtrip structure
+
+lvdat = reshape(result.boot_result.compare_u(:,1), num_chans, num_freqs, num_time);
+%udat = reshape(result.u, num_chans, num_freqs, num_time);
+
+stat = [];
+stat.prob = lvdat;
+stat.dimord = 'chan_freq_time';
+stat.clusters = [];
+stat.clusters.prob = result.perm_result.sprob; % check for significance of LV
+stat.mask = lvdat > 3 | lvdat < -3;
+stat.cfg = option;
+stat.time = time;
+
+save(fullfile(pn.data, 'c01_taskpls_erp.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')

code/b2a_taskPLS_novelty_frontal_erp_plotLV1.m

@@ -0,0 +1,181 @@
+% Create an overview plot featuring the results of the multivariate PLS
+% comparing spectral changes during the stimulus period under load
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data         = fullfile(rootpath, 'data', 'stats');
+pn.figures      = fullfile(rootpath, 'figures');
+pn.tools        = fullfile(rootpath, 'tools');
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'fieldtrip')); ft_defaults;
+    addpath(genpath(fullfile(pn.tools, 'RainCloudPlots')));
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'winsor'));
+
+% set custom colormap
+cBrew = brewermap(500,'RdBu');
+cBrew = flipud(cBrew);
+colormap(cBrew)
+
+load(fullfile(pn.data, 'c01_taskpls_erp.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')
+load(fullfile(rootpath, 'data','erp', ['sub-001_erp.mat']));
+elec = erp.scene_category{1}.elec;
+
+result.perm_result.sprob
+
+indLV = 1;
+
+lvdat = reshape(result.boot_result.compare_u(:,indLV), num_chans, num_freqs, num_time);
+stat.prob = lvdat;
+stat.mask = lvdat > 3 | lvdat < -3;
+
+% maskNaN = double(stat.mask);
+% maskNaN(maskNaN==0) = NaN;
+
+%% invert solution
+
+% stat.mask = stat.mask;
+% stat.prob = stat.prob.*-1;
+% result.vsc = result.vsc.*-1;
+% result.usc = result.usc.*-1;
+
+h = figure('units','centimeter','position',[0 0 15 10]);
+set(gcf,'renderer','Painters')
+statsPlot = [];
+statsPlot = cat(1, statsPlot,squeeze(nanmax(abs(stat.prob(1:64,:,:)),[],1)));
+plot(stat.time,statsPlot, 'k')
+xlabel('Time [s from stim onset]'); ylabel('max abs BSR');
+title({'ERP changes'; ['p = ', num2str(round(result.perm_result.sprob(indLV),4))]})
+set(findall(gcf,'-property','FontSize'),'FontSize',18)
+% figureName = ['b01_pls_traces'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');
+
+
+h = figure('units','centimeter','position',[0 0 15 10]);
+set(gcf,'renderer','Painters')
+statsPlot = [];
+statsPlot = cat(1, statsPlot,squeeze(mean(stat.prob(1:64,:,:),1)));
+plot(stat.time,statsPlot, 'k')
+xlabel('Time [s from stim onset]'); ylabel('mean BSR');
+title({'ERP changes'; ['p = ', num2str(round(result.perm_result.sprob(indLV),4))]})
+set(findall(gcf,'-property','FontSize'),'FontSize',18)
+%xlim([0 0.3])
+% figureName = ['b01_pls_traces'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');
+
+%% plot multivariate topographies
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+cfg.zlim = [-2 2]; 
+cfg.figure = h;
+plotData.powspctrm = squeeze(nanmean(stat.mask(:,:,:).*...
+    stat.prob(:,:,:),3)); 
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Mean BSR');
+figureName = ['b01_lv1'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');
+
+%% plot using raincloud plot
+
+groupsizes=result.num_subj_lst;
+conditions=lv_evt_list;
+conds = {'old'; 'new'};
+condData = []; uData = [];
+for indGroup = 1
+    if indGroup == 1
+        relevantEntries = 1:groupsizes(1)*numel(conds);
+    elseif indGroup == 2
+        relevantEntries = groupsizes(1)*numel(conds)+1:...
+             groupsizes(1)*numel(conds)+groupsizes(2)*numel(conds);
+    end
+    for indCond = 1:numel(conds)
+        targetEntries = relevantEntries(conditions(relevantEntries)==indCond);        
+        condData{indGroup}(indCond,:) = result.vsc(targetEntries,indLV);
+        uData{indGroup}(indCond,:) = result.usc(targetEntries,indLV);
+    end
+end
+
+%% plot RainCloudPlot (within-subject centered)
+
+cBrew(1,:) = 2.*[.3 .1 .1];
+cBrew(2,:) = [.6 .6 .6];
+
+idx_outlier = cell(1); idx_standard = cell(1);
+for indGroup = 1
+    dataToPlot = uData{indGroup}';
+    % define outlier as lin. modulation that is more than three scaled median absolute deviations (MAD) away from the median
+    X = [1 1; 1 2]; b=X\dataToPlot'; IndividualSlopes = b(2,:);
+    outliers = isoutlier(IndividualSlopes, 'median');
+    idx_outlier{indGroup} = find(outliers);
+    idx_standard{indGroup} = find(outliers==0);
+end
+
+h = figure('units','centimeter','position',[0 0 25 10]);
+for indGroup = 1
+    dataToPlot = uData{indGroup}';
+    % read into cell array of the appropriate dimensions
+    data = []; data_ws = [];
+    for i = 1:2
+        for j = 1:1
+            data{i, j} = dataToPlot(:,i);
+            % individually demean for within-subject visualization
+            data_ws{i, j} = dataToPlot(:,i)-...
+                nanmean(dataToPlot(:,:),2)+...
+                repmat(nanmean(nanmean(dataToPlot(:,:),2),1),size(dataToPlot(:,:),1),1);
+            data_nooutlier{i, j} = data{i, j};
+            data_nooutlier{i, j}(idx_outlier{indGroup}) = [];
+            data_ws_nooutlier{i, j} = data_ws{i, j};
+            data_ws_nooutlier{i, j}(idx_outlier{indGroup}) = [];
+            % sort outliers to back in original data for improved plot overlap
+            data_ws{i, j} = [data_ws{i, j}(idx_standard{indGroup}); data_ws{i, j}(idx_outlier{indGroup})];
+        end
+    end
+
+    % IMPORTANT: plot individually centered estimates, stats on uncentered estimates!
+
+    subplot(1,2,indGroup);
+    set(gcf,'renderer','Painters')
+        cla;
+        cl = cBrew(indGroup,:);
+        [~, dist] = rm_raincloud_fixedSpacing(data_ws, [.8 .8 .8],1);
+        h_rc = rm_raincloud_fixedSpacing(data_ws_nooutlier, cl,1,[],[],[],dist);
+        view([90 -90]);
+        axis ij
+    box(gca,'off')
+    set(gca, 'YTickLabels', {conds{2}; conds{1}});
+    yticks = get(gca, 'ytick'); ylim([yticks(1)-(yticks(2)-yticks(1))./2, yticks(2)+(yticks(2)-yticks(1))./1.5]);
+
+    minmax = [min(min(cat(2,data_ws{:}))), max(max(cat(2,data_ws{:})))];
+    xlim(minmax+[-0.2*diff(minmax), 0.2*diff(minmax)])
+    ylabel('novelty'); xlabel({'Brainscore'; '[Individually centered]'})
+
+    % test linear effect
+    curData = [data_nooutlier{1, 1}, data_nooutlier{2, 1}];
+    X = [1 1; 1 2]; b=X\curData'; IndividualSlopes = b(2,:);
+    [~, p, ci, stats] = ttest(IndividualSlopes);
+    title(['M:', num2str(round(mean(IndividualSlopes),3)), '; p=', num2str(round(p,3))])
+end
+% figureName = ['b01_rcp'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');

code/b2a_taskPLS_novelty_frontal_erp_plot_trace.m

@@ -0,0 +1,156 @@
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'erp');
+pn.data_erf = fullfile(rootpath, 'data', 'erf');
+pn.tools = fullfile(rootpath, 'tools');
+    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'shadedErrorBar'));
+    
+%% load erp
+
+for ind_id = 1:33
+    id = sprintf('sub-%03d', ind_id);
+    load(fullfile(pn.data_erp, [id,'_erp_bl.mat']));
+    for ind_option = 1:numel(conds.old)
+        if ind_id == 1
+             erpgroup.old.(conds.old{ind_option}) = erp_bl.old{ind_option};
+             erpgroup.old.(conds.old{ind_option}) = ...
+                 rmfield(erpgroup.old.(conds.old{ind_option}), {'avg', 'var', 'dof'});
+             erpgroup.old.(conds.old{ind_option}).dimord = 'sub_chan_time';
+        end
+        erpgroup.old.(conds.old{ind_option}).avg(ind_id,:,:) = erp_bl.old{ind_option}.avg;
+    end
+end
+
+time = erpgroup.old.old.time;
+elec = erpgroup.old.old.elec;
+channels = erpgroup.old.old.label;
+
+%idx_chans = find(ismember(channels, {'O1', 'Oz', 'O2'}));
+%idx_chans = find(ismember(channels, {'PO7', 'PO8'}));
+%idx_chans = find(ismember(channels, {'Pz', 'CPz', 'P1'}));
+
+mergeddata = cat(4, erpgroup.old.old.avg, ...
+    erpgroup.old.new.avg);
+
+%% plot ERP topography
+
+% set custom colormap
+cBrew = brewermap(500,'RdBu');
+cBrew = flipud(cBrew);
+colormap(cBrew)
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+plotData.powspctrm = squeeze(nanmean(nanmean(nanmean(mergeddata(:,:,time>=0.3& time <=0.5,1:2),3),1),4))'; 
+
+cfg.zlim = [-14 14]*10^-4; 
+cfg.figure = h;
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Amplitude');
+% figureName = ['xxx'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');
+
+%% visualize frontal ERP
+
+idx_chans = [38];
+
+% avg across channels and conditions
+condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
+
+h = figure('units','centimeters','position',[0 0 10 8]);
+cla; hold on;
+% highlight relevant phase in background
+patches.timeVec = [0.3 0.5];
+patches.colorVec = [1 .95 .8];
+for indP = 1:size(patches.timeVec,2)-1
+    YLim = [-8 2]*10^-4;
+    p = patch([patches.timeVec(indP) patches.timeVec(indP+1) patches.timeVec(indP+1) patches.timeVec(indP)], ...
+                [YLim(1) YLim(1)  YLim(2), YLim(2)], patches.colorVec(indP,:));
+    p.EdgeColor = 'none';
+end
+% new value = old value ? subject average + grand average
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
+% ax = gca; ax.YDir = 'reverse';
+legend([l1.mainLine, l2.mainLine],{'old', 'new'}, ...
+    'location', 'southwest'); legend('boxoff')
+xlabel('Time (s) from stim onset')
+xlim([-.05 .6]); %ylim([-.03 .18])
+ylabel({'ERP';'(microVolts)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',14)
+
+%% plot difference
+
+idx_chans = [38];
+
+h = figure('units','centimeters','position',[0 0 10 8]);
+cla; hold on;
+% highlight relevant phase in background
+patches.timeVec = [0.3 0.5];
+patches.colorVec = [1 .95 .8];
+for indP = 1:size(patches.timeVec,2)-1
+    YLim = [-10 5]*10^-5;
+    p = patch([patches.timeVec(indP) patches.timeVec(indP+1) patches.timeVec(indP+1) patches.timeVec(indP)], ...
+                [YLim(1) YLim(1)  YLim(2), YLim(2)], patches.colorVec(indP,:));
+    p.EdgeColor = 'none';
+end
+% new value = old value ? subject average + grand average
+curData = squeeze(mergeddata(:,idx_chans,:,2)-mergeddata(:,idx_chans,:,1));
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+legend([l1.mainLine],{'new-old'}, ...
+    'location', 'southwest'); legend('boxoff')
+xlabel('Time (s) from stim onset')
+xlim([-.05 .6]); %ylim([-.03 .18])
+ylabel({'ERP';'(microVolts)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',14)
+
+
+%% plot difference between conditions
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+oldminnew = mergeddata(:,:,:,1)-mergeddata(:,:,:,2);
+plotData.powspctrm = squeeze(nanmean(nanmean(oldminnew(:,:,time>=.3 & time <=.5),3),1))'; 
+
+cfg.figure = h;
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Amplitude');

code/b2b_taskPLS_novelty_frontal_theta.m

@@ -0,0 +1,106 @@
+% Set up EEG PLS for a task PLS using spectral power
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'erp');
+pn.data_erf = fullfile(rootpath, 'data', 'erf');
+pn.data = fullfile(rootpath, 'data', 'stats'); mkdir(pn.data);
+pn.tools = fullfile(rootpath, 'tools');
+    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'shadedErrorBar'));
+
+%% add seed for reproducibility
+
+rng(0, 'twister');
+
+%% load erp
+
+for ind_id = 1:33
+    id = sprintf('sub-%03d', ind_id); disp(id)
+    load(fullfile(pn.data_erf, [id,'_erf.mat']));
+    for ind_option = 1:numel(conds.old)
+        if ind_id == 1
+             erpgroup.old.(conds.old{ind_option}) = erf.old{ind_option};
+             erpgroup.old.(conds.old{ind_option}) = ...
+                 rmfield(erpgroup.old.(conds.old{ind_option}), {'powspctrm'});
+             erpgroup.old.(conds.old{ind_option}).dimord = 'sub_chan_freq_time';
+             freq = erpgroup.old.(conds.old{ind_option}).freq;
+             idx_freq = freq >2 & freq <8;
+        end
+        % only include time points from 300 - 500 ms
+        time = erf.old{ind_option}.time;
+        toi = time >0.3 & time <0.5;        
+        erpgroup.old.(conds.old{ind_option}).avg(ind_id,:,:,:) = erf.old{ind_option}.powspctrm(:,idx_freq,toi);
+    end
+end
+
+time = erpgroup.old.old.time(toi);
+freq = erpgroup.old.old.freq(idx_freq);
+channels = erpgroup.old.old.label;
+
+mergeddata = cat(5, erpgroup.old.old.avg, erpgroup.old.new.avg);
+
+% identify frontal channels, and set everything else to zero
+idx_chans = startsWith(channels, 'F') | startsWith(channels, 'A') | startsWith(channels, 'C') & ~startsWith(channels, 'CP');
+mergeddata(:, idx_chans==0,:,:,:) = 0;
+
+num_chans = numel(channels);
+num_freqs = numel(freq);
+num_time = numel(time);
+
+num_subj_lst = [33];
+num_cond = 2;
+num_grp = 1;
+
+datamat_lst = cell(num_grp); lv_evt_list = [];
+indCount_cont = 1;
+indCount = 1;
+indGroup = 1;
+for indCond = 1:num_cond
+    for indID = 1:num_subj_lst(indGroup)
+        datamat_lst{indGroup}(indCount,:) = reshape(squeeze(mergeddata(indID,:,:,:,indCond)), [], 1);
+        lv_evt_list(indCount_cont) = indCond;
+        indCount = indCount+1;
+        indCount_cont = indCount_cont+1;
+    end
+end
+datamat_lst{indGroup}(isnan(datamat_lst{indGroup})) = 0;
+
+%% set PLS options and run PLS
+
+option = [];
+option.method = 1; % [1] | 2 | 3 | 4 | 5 | 6
+option.num_perm = 1000; %( single non-negative integer )
+option.num_split = 0; %( single non-negative integer )
+option.num_boot = 1000; % ( single non-negative integer )
+option.cormode = 0; % [0] | 2 | 4 | 6
+option.meancentering_type = 0;% [0] | 1 | 2 | 3
+option.boot_type = 'strat'; %['strat'] | 'nonstrat'
+
+result = pls_analysis(datamat_lst, num_subj_lst, num_cond, option);
+
+%% rearrange into fieldtrip structure
+
+lvdat = reshape(result.boot_result.compare_u(:,1), num_chans, num_freqs, num_time);
+%udat = reshape(result.u, num_chans, num_freqs, num_time);
+
+stat = [];
+stat.prob = lvdat;
+stat.dimord = 'chan_freq_time';
+stat.clusters = [];
+stat.clusters.prob = result.perm_result.sprob; % check for significance of LV
+stat.mask = lvdat > 3 | lvdat < -3;
+stat.cfg = option;
+stat.freq = freq;
+stat.time = time;
+
+save(fullfile(pn.data, 'b2b_taskpls_theta.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')

code/b2b_taskPLS_novelty_frontal_theta_plotLV1.m

@@ -0,0 +1,189 @@
+% Create an overview plot featuring the results of the multivariate PLS
+% comparing spectral changes during the stimulus period under load
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data         = fullfile(rootpath, 'data', 'stats');
+pn.figures      = fullfile(rootpath, 'figures');
+pn.tools        = fullfile(rootpath, 'tools');
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'fieldtrip')); ft_defaults;
+    addpath(genpath(fullfile(pn.tools, 'RainCloudPlots')));
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'winsor'));
+
+% set custom colormap
+cBrew = brewermap(500,'RdBu');
+cBrew = flipud(cBrew);
+colormap(cBrew)
+
+load(fullfile(pn.data, 'b2b_taskpls_theta.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')
+load(fullfile(rootpath, 'data','erp', ['sub-001_erp.mat']));
+elec = erp.scene_category{1}.elec;
+
+result.perm_result.sprob
+
+indLV = 1;
+
+lvdat = reshape(result.boot_result.compare_u(:,indLV), num_chans, num_freqs, num_time);
+stat.prob = lvdat;
+stat.mask = lvdat > 3 | lvdat < -3;
+
+% maskNaN = double(stat.mask);
+% maskNaN(maskNaN==0) = NaN;
+
+%% invert solution
+
+stat.mask = stat.mask;
+stat.prob = stat.prob.*-1;
+result.vsc = result.vsc.*-1;
+result.usc = result.usc.*-1;
+
+% h = figure('units','centimeter','position',[0 0 15 10]);
+% set(gcf,'renderer','Painters')
+% statsPlot = [];
+% statsPlot = cat(1, statsPlot,squeeze(nanmax(abs(stat.prob(1:64,:,:)),[],1)));
+% plot(stat.time,statsPlot, 'k')
+% xlabel('Time [s from stim onset]'); ylabel('max abs BSR');
+% title({'ERP changes'; ['p = ', num2str(round(result.perm_result.sprob(indLV),4))]})
+% set(findall(gcf,'-property','FontSize'),'FontSize',18)
+% figureName = ['b01_pls_traces'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');
+
+% h = figure('units','centimeter','position',[0 0 15 10]);
+% set(gcf,'renderer','Painters')
+% statsPlot = [];
+% statsPlot = cat(1, statsPlot,squeeze(mean(stat.prob(1:64,:,:),1)));
+% plot(stat.time,statsPlot, 'k')
+% xlabel('Time [s from stim onset]'); ylabel('mean BSR');
+% title({'ERP changes'; ['p = ', num2str(round(result.perm_result.sprob(indLV),4))]})
+% set(findall(gcf,'-property','FontSize'),'FontSize',18)
+% %xlim([0 0.3])
+% % figureName = ['b01_pls_traces'];
+% % saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% % saveas(h, fullfile(pn.figures, figureName), 'png');
+
+%% plot multivariate topographies
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+cfg.zlim = [-6 6]; 
+cfg.figure = h;
+plotData.powspctrm = squeeze(nanmax(nanmax(stat.mask(:,:,:).*...
+    stat.prob(:,:,:),[],3),[],2));
+
+[~, sortind] = sort(plotData.powspctrm, 'descend');
+cfg.marker = 'off'; 
+cfg.highlight = 'yes';
+cfg.highlightchannel = plotData.label(sortind(1:3));
+cfg.highlightcolor = [0 0 0];
+cfg.highlightsymbol = '.';
+cfg.highlightsize = 15;
+
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Max BSR');
+figureName = ['b01_lv1'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');
+
+%% plot using raincloud plot
+
+groupsizes=result.num_subj_lst;
+conditions=lv_evt_list;
+conds = {'old'; 'new'};
+condData = []; uData = [];
+for indGroup = 1
+    if indGroup == 1
+        relevantEntries = 1:groupsizes(1)*numel(conds);
+    elseif indGroup == 2
+        relevantEntries = groupsizes(1)*numel(conds)+1:...
+             groupsizes(1)*numel(conds)+groupsizes(2)*numel(conds);
+    end
+    for indCond = 1:numel(conds)
+        targetEntries = relevantEntries(conditions(relevantEntries)==indCond);        
+        condData{indGroup}(indCond,:) = -1.*result.vsc(targetEntries,indLV);
+        uData{indGroup}(indCond,:) = -1.*result.usc(targetEntries,indLV);
+    end
+end
+
+%% plot RainCloudPlot (within-subject centered)
+
+cBrew(1,:) = 2.*[.3 .1 .1];
+cBrew(2,:) = [.6 .6 .6];
+
+idx_outlier = cell(1); idx_standard = cell(1);
+for indGroup = 1
+    dataToPlot = uData{indGroup}';
+    % define outlier as lin. modulation that is more than three scaled median absolute deviations (MAD) away from the median
+    X = [1 1; 1 2]; b=X\dataToPlot'; IndividualSlopes = b(2,:);
+    outliers = isoutlier(IndividualSlopes, 'median');
+    idx_outlier{indGroup} = find(outliers);
+    idx_standard{indGroup} = find(outliers==0);
+end
+
+h = figure('units','centimeter','position',[0 0 25 10]);
+for indGroup = 1
+    dataToPlot = uData{indGroup}';
+    % read into cell array of the appropriate dimensions
+    data = []; data_ws = [];
+    for i = 1:2
+        for j = 1:1
+            data{i, j} = dataToPlot(:,i);
+            % individually demean for within-subject visualization
+            data_ws{i, j} = dataToPlot(:,i)-...
+                nanmean(dataToPlot(:,:),2)+...
+                repmat(nanmean(nanmean(dataToPlot(:,:),2),1),size(dataToPlot(:,:),1),1);
+            data_nooutlier{i, j} = data{i, j};
+            data_nooutlier{i, j}(idx_outlier{indGroup}) = [];
+            data_ws_nooutlier{i, j} = data_ws{i, j};
+            data_ws_nooutlier{i, j}(idx_outlier{indGroup}) = [];
+            % sort outliers to back in original data for improved plot overlap
+            data_ws{i, j} = [data_ws{i, j}(idx_standard{indGroup}); data_ws{i, j}(idx_outlier{indGroup})];
+        end
+    end
+
+    % IMPORTANT: plot individually centered estimates, stats on uncentered estimates!
+
+    subplot(1,2,indGroup);
+    set(gcf,'renderer','Painters')
+        cla;
+        cl = cBrew(indGroup,:);
+        [~, dist] = rm_raincloud_fixedSpacing(data_ws, [.8 .8 .8],1);
+        h_rc = rm_raincloud_fixedSpacing(data_ws_nooutlier, cl,1,[],[],[],dist);
+        view([90 -90]);
+        axis ij
+    box(gca,'off')
+    set(gca, 'YTickLabels', {conds{2}; conds{1}});
+    yticks = get(gca, 'ytick'); ylim([yticks(1)-(yticks(2)-yticks(1))./2, yticks(2)+(yticks(2)-yticks(1))./1.5]);
+
+    minmax = [min(min(cat(2,data_ws{:}))), max(max(cat(2,data_ws{:})))];
+    xlim(minmax+[-0.2*diff(minmax), 0.2*diff(minmax)])
+    ylabel('novelty'); xlabel({'Brainscore'; '[Individually centered]'})
+
+    % test linear effect
+    curData = [data_nooutlier{1, 1}, data_nooutlier{2, 1}];
+    X = [1 1; 1 2]; b=X\curData'; IndividualSlopes = b(2,:);
+    [~, p, ci, stats] = ttest(IndividualSlopes);
+    title(['M:', num2str(round(mean(IndividualSlopes),3)), '; p=', num2str(round(p,3))])
+end
+% figureName = ['b01_rcp'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');

code/b2c_taskPLS_novelty_posterior_alpha.m

@@ -0,0 +1,116 @@
+% Set up EEG PLS for a task PLS using spectral power
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'erp');
+pn.data_erf = fullfile(rootpath, 'data', 'erf');
+pn.data = fullfile(rootpath, 'data', 'stats'); mkdir(pn.data);
+pn.tools = fullfile(rootpath, 'tools');
+    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'shadedErrorBar'));
+
+%% add seed for reproducibility
+
+rng(0, 'twister');
+    
+%% load event info
+
+load(fullfile(pn.data_eeg, ['sub-001_task-xxxx_eeg_art.mat']), 'events');
+
+parameter = {'scene_category'; 'old'; 'behavior'; 'subsequent_memory'};
+for ind_param = 1:numel(parameter)
+    conds.(parameter{ind_param}) = unique(events.(parameter{ind_param}));
+end
+
+%% load erp
+
+for ind_id = 1:33
+    id = sprintf('sub-%03d', ind_id); disp(id)
+    load(fullfile(pn.data_erf, [id,'_erf.mat']));
+    for ind_option = 1:numel(conds.old)
+        if ind_id == 1
+             erpgroup.old.(conds.old{ind_option}) = erf.old{ind_option};
+             erpgroup.old.(conds.old{ind_option}) = ...
+                 rmfield(erpgroup.old.(conds.old{ind_option}), {'powspctrm'});
+             erpgroup.old.(conds.old{ind_option}).dimord = 'sub_chan_freq_time';
+             freq = erpgroup.old.(conds.old{ind_option}).freq;
+             idx_freq = freq >7 & freq <13;
+        end
+        % only include time points from 300 - 500 ms
+        time = erf.old{ind_option}.time;
+        toi = time >0.3 & time <0.5;        
+        erpgroup.old.(conds.old{ind_option}).avg(ind_id,:,:,:) = erf.old{ind_option}.powspctrm(:,idx_freq,toi);
+    end
+end
+
+time = erpgroup.old.old.time(toi);
+freq = erpgroup.old.old.freq(idx_freq);
+channels = erpgroup.old.old.label;
+
+mergeddata = cat(5, erpgroup.old.old.avg, ...
+    erpgroup.old.new.avg);
+
+% identify parieto-occipital channels, and set everything else to zero
+idx_chans = startsWith(channels, 'P') | startsWith(channels, 'O');
+mergeddata(:, idx_chans==0,:,:,:) = 0;
+
+num_chans = numel(channels);
+num_freqs = numel(freq);
+num_time = numel(time);
+
+num_subj_lst = [33];
+num_cond = 2;
+num_grp = 1;
+
+datamat_lst = cell(num_grp); lv_evt_list = [];
+indCount_cont = 1;
+indCount = 1;
+indGroup = 1;
+for indCond = 1:num_cond
+    for indID = 1:num_subj_lst(indGroup)
+        datamat_lst{indGroup}(indCount,:) = reshape(squeeze(mergeddata(indID,:,:,:,indCond)), [], 1);
+        lv_evt_list(indCount_cont) = indCond;
+        indCount = indCount+1;
+        indCount_cont = indCount_cont+1;
+    end
+end
+datamat_lst{indGroup}(isnan(datamat_lst{indGroup})) = 0;
+
+%% set PLS options and run PLS
+
+option = [];
+option.method = 1; % [1] | 2 | 3 | 4 | 5 | 6
+option.num_perm = 1000; %( single non-negative integer )
+option.num_split = 0; %( single non-negative integer )
+option.num_boot = 1000; % ( single non-negative integer )
+option.cormode = 0; % [0] | 2 | 4 | 6
+option.meancentering_type = 0;% [0] | 1 | 2 | 3
+option.boot_type = 'strat'; %['strat'] | 'nonstrat'
+
+result = pls_analysis(datamat_lst, num_subj_lst, num_cond, option);
+
+%% rearrange into fieldtrip structure
+
+lvdat = reshape(result.boot_result.compare_u(:,1), num_chans, num_freqs, num_time);
+%udat = reshape(result.u, num_chans, num_freqs, num_time);
+
+stat = [];
+stat.prob = lvdat;
+stat.dimord = 'chan_freq_time';
+stat.clusters = [];
+stat.clusters.prob = result.perm_result.sprob; % check for significance of LV
+stat.mask = lvdat > 3 | lvdat < -3;
+stat.cfg = option;
+stat.freq = freq;
+stat.time = time;
+
+save(fullfile(pn.data, 'b2c_taskpls_alpha.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')

code/b2c_taskPLS_novelty_posterior_alpha_LV1.m

@@ -0,0 +1,179 @@
+% Create an overview plot featuring the results of the multivariate PLS
+% comparing spectral changes during the stimulus period under load
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data         = fullfile(rootpath, 'data', 'stats');
+pn.figures      = fullfile(rootpath, 'figures');
+pn.tools        = fullfile(rootpath, 'tools');
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'fieldtrip')); ft_defaults;
+    addpath(genpath(fullfile(pn.tools, 'RainCloudPlots')));
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'winsor'));
+
+% set custom colormap
+cBrew = brewermap(500,'RdBu');
+cBrew = flipud(cBrew);
+colormap(cBrew)
+
+load(fullfile(pn.data, 'b2c_taskpls_alpha.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')
+load(fullfile(rootpath, 'data','erp', ['sub-001_erp.mat']));
+elec = erp.scene_category{1}.elec;
+
+result.perm_result.sprob
+
+indLV = 1;
+
+lvdat = reshape(result.boot_result.compare_u(:,indLV), num_chans, num_freqs, num_time);
+stat.prob = lvdat;
+stat.mask = lvdat > 3 | lvdat < -3;
+
+% maskNaN = double(stat.mask);
+% maskNaN(maskNaN==0) = NaN;
+
+%% invert solution
+
+% stat.mask = stat.mask;
+% stat.prob = stat.prob.*-1;
+% result.vsc = result.vsc.*-1;
+% result.usc = result.usc.*-1;
+
+h = figure('units','centimeter','position',[0 0 15 10]);
+set(gcf,'renderer','Painters')
+statsPlot = [];
+statsPlot = cat(1, statsPlot,squeeze(nanmax(abs(stat.prob(1:64,:,:)),[],1)));
+plot(stat.time,statsPlot, 'k')
+xlabel('Time [s from stim onset]'); ylabel('max abs BSR');
+title({'ERP changes'; ['p = ', num2str(round(result.perm_result.sprob(indLV),4))]})
+set(findall(gcf,'-property','FontSize'),'FontSize',18)
+% figureName = ['b01_pls_traces'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');
+
+
+h = figure('units','centimeter','position',[0 0 15 10]);
+set(gcf,'renderer','Painters')
+statsPlot = [];
+statsPlot = cat(1, statsPlot,squeeze(mean(stat.prob(1:64,:,:),1)));
+plot(stat.time,statsPlot, 'k')
+xlabel('Time [s from stim onset]'); ylabel('mean BSR');
+title({'ERP changes'; ['p = ', num2str(round(result.perm_result.sprob(indLV),4))]})
+set(findall(gcf,'-property','FontSize'),'FontSize',18)
+%xlim([0 0.3])
+% figureName = ['b01_pls_traces'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');
+
+%% plot multivariate topographies
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+cfg.zlim = [-2 2]; 
+cfg.figure = h;
+plotData.powspctrm = squeeze(nanmean(stat.mask(:,:,:).*...
+    stat.prob(:,:,:),3)); 
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Mean BSR');
+figureName = ['b01_lv1'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');
+
+%% plot RainCloudPlot (within-subject centered)
+
+groupsizes=result.num_subj_lst;
+conditions=lv_evt_list;
+conds = {'manmade'; 'natural'};
+condData = []; uData = [];
+for indGroup = 1
+    if indGroup == 1
+        relevantEntries = 1:groupsizes(1)*numel(conds);
+    elseif indGroup == 2
+        relevantEntries = groupsizes(1)*numel(conds)+1:...
+             groupsizes(1)*numel(conds)+groupsizes(2)*numel(conds);
+    end
+    for indCond = 1:numel(conds)
+        targetEntries = relevantEntries(conditions(relevantEntries)==indCond);        
+        condData{indGroup}(indCond,:) = result.vsc(targetEntries,indLV);
+        uData{indGroup}(indCond,:) = result.usc(targetEntries,indLV);
+    end
+end
+
+cBrew(1,:) = 2.*[.3 .1 .1];
+cBrew(2,:) = [.6 .6 .6];
+
+idx_outlier = cell(1); idx_standard = cell(1);
+for indGroup = 1
+    dataToPlot = uData{indGroup}';
+    % define outlier as lin. modulation that is more than three scaled median absolute deviations (MAD) away from the median
+    X = [1 1; 1 2]; b=X\dataToPlot'; IndividualSlopes = b(2,:);
+    outliers = isoutlier(IndividualSlopes, 'median');
+    idx_outlier{indGroup} = find(outliers);
+    idx_standard{indGroup} = find(outliers==0);
+end
+
+h = figure('units','centimeter','position',[0 0 25 10]);
+for indGroup = 1
+    dataToPlot = uData{indGroup}';
+    % read into cell array of the appropriate dimensions
+    data = []; data_ws = [];
+    for i = 1:2
+        for j = 1:1
+            data{i, j} = dataToPlot(:,i);
+            % individually demean for within-subject visualization
+            data_ws{i, j} = dataToPlot(:,i)-...
+                nanmean(dataToPlot(:,:),2)+...
+                repmat(nanmean(nanmean(dataToPlot(:,:),2),1),size(dataToPlot(:,:),1),1);
+            data_nooutlier{i, j} = data{i, j};
+            data_nooutlier{i, j}(idx_outlier{indGroup}) = [];
+            data_ws_nooutlier{i, j} = data_ws{i, j};
+            data_ws_nooutlier{i, j}(idx_outlier{indGroup}) = [];
+            % sort outliers to back in original data for improved plot overlap
+            data_ws{i, j} = [data_ws{i, j}(idx_standard{indGroup}); data_ws{i, j}(idx_outlier{indGroup})];
+        end
+    end
+
+    % IMPORTANT: plot individually centered estimates, stats on uncentered estimates!
+
+    subplot(1,2,indGroup);
+    set(gcf,'renderer','Painters')
+        cla;
+        cl = cBrew(indGroup,:);
+        [~, dist] = rm_raincloud_fixedSpacing(data_ws, [.8 .8 .8],1);
+        h_rc = rm_raincloud_fixedSpacing(data_ws_nooutlier, cl,1,[],[],[],dist);
+        view([90 -90]);
+        axis ij
+    box(gca,'off')
+    set(gca, 'YTickLabels', {'natural'; 'manmade'});
+    yticks = get(gca, 'ytick'); ylim([yticks(1)-(yticks(2)-yticks(1))./2, yticks(2)+(yticks(2)-yticks(1))./1.5]);
+
+    minmax = [min(min(cat(2,data_ws{:}))), max(max(cat(2,data_ws{:})))];
+    xlim(minmax+[-0.2*diff(minmax), 0.2*diff(minmax)])
+    ylabel('scene category'); xlabel({'Brainscore'; '[Individually centered]'})
+
+    % test linear effect
+    curData = [data_nooutlier{1, 1}, data_nooutlier{2, 1}];
+    X = [1 1; 1 2]; b=X\curData'; IndividualSlopes = b(2,:);
+    [~, p, ci, stats] = ttest(IndividualSlopes);
+    title(['M:', num2str(round(mean(IndividualSlopes),3)), '; p=', num2str(round(p,3))])
+end
+% figureName = ['b01_rcp'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');

code/b2x_novelty_plot_theta_alpha.m

@@ -0,0 +1,127 @@
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'erp');
+pn.data_erf = fullfile(rootpath, 'data', 'erf');
+pn.data = fullfile(rootpath, 'data', 'stats'); mkdir(pn.data);
+pn.tools = fullfile(rootpath, 'tools');
+    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'shadedErrorBar'));
+
+%% add seed for reproducibility
+
+rng(0, 'twister');
+
+%% load erp
+
+for ind_id = 1:33
+    id = sprintf('sub-%03d', ind_id); disp(id)
+    load(fullfile(pn.data_erf, [id,'_erf.mat']));
+    for ind_option = 1:numel(conds.old)
+        if ind_id == 1
+             erpgroup.old.(conds.old{ind_option}) = erf.old{ind_option};
+             erpgroup.old.(conds.old{ind_option}) = ...
+                 rmfield(erpgroup.old.(conds.old{ind_option}), {'powspctrm'});
+             erpgroup.old.(conds.old{ind_option}).dimord = 'sub_chan_freq_time';
+             freq = erpgroup.old.(conds.old{ind_option}).freq;
+        end
+        idx_freq_t = freq >2 & freq <8;
+        thetagroup.old.(conds.old{ind_option}).avg(ind_id,:,:,:) = squeeze(nanmean(erf.old{ind_option}.powspctrm(:,idx_freq_t,:),2));
+        idx_freq_a = freq >7 & freq <13;
+        alphagroup.old.(conds.old{ind_option}).avg(ind_id,:,:,:) = squeeze(nanmean(erf.old{ind_option}.powspctrm(:,idx_freq_a,:),2));
+    end
+end
+
+time = erpgroup.old.old.time;
+channels = erpgroup.old.old.label;
+
+% identify frontal channels
+%idx_chans_t = startsWith(channels, 'F') | startsWith(channels, 'A') | startsWith(channels, 'C') & ~startsWith(channels, 'CP');
+
+% identify parieto-occipital channels
+idx_chans_a = startsWith(channels, 'P') | startsWith(channels, 'O');
+
+%% visualize for frontal theta
+
+mergeddata = cat(4, thetagroup.old.old.avg, thetagroup.old.new.avg);
+
+%figure; imagesc(squeeze(nanmean(thetagroup.old.new.avg-thetagroup.old.old.avg,1)))
+
+idx_chans = [33,38,47];
+condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
+
+h = figure('units','centimeters','position',[0 0 10 8]);
+cla; hold on;
+% highlight relevant phase in background
+patches.timeVec = [0.3 0.5];
+patches.colorVec = [1 .95 .8];
+for indP = 1:size(patches.timeVec,2)-1
+    YLim = [-5 -4.8];
+    p = patch([patches.timeVec(indP) patches.timeVec(indP+1) patches.timeVec(indP+1) patches.timeVec(indP)], ...
+                [YLim(1) YLim(1)  YLim(2), YLim(2)], patches.colorVec(indP,:));
+    p.EdgeColor = 'none';
+end
+% new value = old value ? subject average + grand average
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
+    {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
+    {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
+xlabel('Time (s) from stim onset')
+legend([l1.mainLine, l2.mainLine],{'old', 'new'}, ...
+    'location', 'southeast'); legend('boxoff')
+xlim([-.5 1.5]); ylim(YLim)
+%xlim([-.05 .3]);
+ylabel({'theta power';'(log10)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',12)
+
+%% visualize for posterior alpha
+
+mergeddata = cat(4, alphagroup.old.old.avg, alphagroup.old.new.avg);
+
+idx_chans = idx_chans_a;
+condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
+
+h = figure('units','centimeters','position',[0 0 10 8]);
+cla; hold on;
+% highlight relevant phase in background
+patches.timeVec = [0.3 0.5];
+patches.colorVec = [1 .95 .8];
+for indP = 1:size(patches.timeVec,2)-1
+    YLim = [-5.2 -4.5];
+    p = patch([patches.timeVec(indP) patches.timeVec(indP+1) patches.timeVec(indP+1) patches.timeVec(indP)], ...
+                [YLim(1) YLim(1)  YLim(2), YLim(2)], patches.colorVec(indP,:));
+    p.EdgeColor = 'none';
+end
+% new value = old value ? subject average + grand average
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
+    {'color', 'k','linewidth', 3}, 'patchSaturation', .1);
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
+    {'color', 'r','linewidth', 3}, 'patchSaturation', .1);
+xlabel('Time (s) from stim onset')
+legend([l1.mainLine, l2.mainLine],{'old', 'new'}, 'location', 'northeast'); legend('boxoff')
+xlim([-.5 1.5]); ylim(YLim)
+%xlim([-.05 .3]);
+ylabel({'alpha power';'(log10)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',15)
+

code/b3a_taskPLS_recognition_erp.m

@@ -0,0 +1,106 @@
+% Set up EEG PLS for a task PLS using spectral power
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'erp');
+pn.data_erf = fullfile(rootpath, 'data', 'erf');
+pn.data = fullfile(rootpath, 'data', 'stats'); mkdir(pn.data);
+pn.tools = fullfile(rootpath, 'tools');
+    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'shadedErrorBar'));
+
+%% add seed for reproducibility
+
+rng(0, 'twister');
+    
+%% load event info
+
+load(fullfile(pn.data_eeg, ['sub-001_task-xxxx_eeg_art.mat']), 'events');
+
+parameter = {'scene_category'; 'old'; 'behavior'; 'subsequent_memory'};
+for ind_param = 1:numel(parameter)
+    conds.(parameter{ind_param}) = unique(events.(parameter{ind_param}));
+end
+
+%% load erp_bl
+
+for ind_id = 1:33
+    id = sprintf('sub-%03d', ind_id); disp(id)
+    load(fullfile(pn.data_erp, [id,'_erp_bl.mat']));
+    for ind_option = 1:4%numel(conds.behavior)
+        if ind_id == 1
+             erpgroup.behavior.(conds.behavior{ind_option}) = erp_bl.behavior{ind_option};
+             erpgroup.behavior.(conds.behavior{ind_option}) = ...
+                 rmfield(erpgroup.behavior.(conds.behavior{ind_option}), {'avg', 'var', 'dof'});
+             erpgroup.behavior.(conds.behavior{ind_option}).dimord = 'sub_chan_time';
+        end
+        erpgroup.behavior.(conds.behavior{ind_option}).avg(ind_id,:,:) = erp_bl.behavior{ind_option}.avg;
+    end
+end
+
+time = erpgroup.behavior.hit.time;
+elec = erpgroup.behavior.hit.elec;
+channels = erpgroup.behavior.hit.label;
+
+mergeddata = cat(4, erpgroup.behavior.hit.avg, ...
+    erpgroup.behavior.miss.avg);
+
+num_chans = numel(channels);
+num_freqs = 1;
+num_time = numel(time);
+
+num_subj_lst = [33];
+num_cond = 2;
+num_grp = 1;
+
+datamat_lst = cell(num_grp); lv_evt_list = [];
+indCount_cont = 1;
+indCount = 1;
+indGroup = 1;
+for indCond = 1:num_cond
+    for indID = 1:num_subj_lst(indGroup)
+        datamat_lst{indGroup}(indCount,:) = reshape(squeeze(mergeddata(indID,:,:,indCond)), [], 1);
+        lv_evt_list(indCount_cont) = indCond;
+        indCount = indCount+1;
+        indCount_cont = indCount_cont+1;
+    end
+end
+datamat_lst{indGroup}(isnan(datamat_lst{indGroup})) = 0;
+
+%% set PLS options and run PLS
+
+option = [];
+option.method = 1; % [1] | 2 | 3 | 4 | 5 | 6
+option.num_perm = 1000; %( single non-negative integer )
+option.num_split = 0; %( single non-negative integer )
+option.num_boot = 1000; % ( single non-negative integer )
+option.cormode = 0; % [0] | 2 | 4 | 6
+option.meancentering_type = 0;% [0] | 1 | 2 | 3
+option.boot_type = 'strat'; %['strat'] | 'nonstrat'
+
+result = pls_analysis(datamat_lst, num_subj_lst, num_cond, option);
+
+%% rearrange into fieldtrip structure
+
+lvdat = reshape(result.boot_result.compare_u(:,1), num_chans, num_freqs, num_time);
+%udat = reshape(result.u, num_chans, num_freqs, num_time);
+
+stat = [];
+stat.prob = lvdat;
+stat.dimord = 'chan_freq_time';
+stat.clusters = [];
+stat.clusters.prob = result.perm_result.sprob; % check for significance of LV
+stat.mask = lvdat > 3 | lvdat < -3;
+stat.cfg = option;
+stat.time = time;
+
+save(fullfile(pn.data, 'd01_taskpls_erp.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')

code/b3a_taskPLS_recognition_erp_LV1.m

@@ -0,0 +1,195 @@
+% Create an overview plot featuring the results of the multivariate PLS
+% comparing spectral changes during the stimulus period under load
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data         = fullfile(rootpath, 'data', 'stats');
+pn.figures      = fullfile(rootpath, 'figures');
+pn.tools        = fullfile(rootpath, 'tools');
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'fieldtrip')); ft_defaults;
+    addpath(genpath(fullfile(pn.tools, 'RainCloudPlots')));
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'winsor'));
+
+% set custom colormap
+cBrew = brewermap(500,'RdBu');
+cBrew = flipud(cBrew);
+colormap(cBrew)
+
+load(fullfile(pn.data, 'd01_taskpls_erp.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')
+load(fullfile(rootpath, 'data','erp', ['sub-001_erp.mat']));
+elec = erp.scene_category{1}.elec;
+
+result.perm_result.sprob
+
+indLV = 1;
+
+lvdat = reshape(result.boot_result.compare_u(:,indLV), num_chans, num_freqs, num_time);
+stat.prob = lvdat;
+stat.mask = lvdat > 3 | lvdat < -3;
+
+% maskNaN = double(stat.mask);
+% maskNaN(maskNaN==0) = NaN;
+
+%% invert solution
+
+stat.mask = stat.mask;
+stat.prob = stat.prob.*-1;
+result.vsc = result.vsc.*-1;
+result.usc = result.usc.*-1;
+
+h = figure('units','centimeter','position',[0 0 15 10]);
+set(gcf,'renderer','Painters')
+hold on;
+statsPlot = [];
+statsPlot = cat(1, statsPlot,squeeze(nanmean(stat.mask(1:64,:,:).*stat.prob(1:64,:,:),1)));
+plot(stat.time,statsPlot, 'k')
+xlabel('Time [s from stim onset]'); ylabel('max abs BSR');
+title({'ERP changes'; ['p = ', num2str(round(result.perm_result.sprob(indLV),4))]})
+set(findall(gcf,'-property','FontSize'),'FontSize',18)
+% figureName = ['b01_pls_traces'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');
+
+% 
+% h = figure('units','centimeter','position',[0 0 15 10]);
+% set(gcf,'renderer','Painters')
+% statsPlot = [];
+% statsPlot = cat(1, statsPlot,squeeze(mean(stat.prob(1:64,:,:),1)));
+% plot(stat.time,statsPlot, 'k')
+% xlabel('Time [s from stim onset]'); ylabel('mean BSR');
+% title({'ERP changes'; ['p = ', num2str(round(result.perm_result.sprob(indLV),4))]})
+% set(findall(gcf,'-property','FontSize'),'FontSize',18)
+%xlim([0 0.3])
+% figureName = ['b01_pls_traces'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');
+
+%% plot multivariate topographies
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+cfg.zlim = [-5 5]; 
+cfg.figure = h;
+plotData.powspctrm = squeeze(nanmean(stat.mask(:,:,stat.time>.5 & stat.time<.8).*...
+    stat.prob(:,:,stat.time>0.5 & stat.time<.8),3)); 
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Mean BSR');
+
+%%
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+cfg.zlim = [-2 2]; 
+cfg.figure = h;
+plotData.powspctrm = squeeze(nanmean(stat.mask(:,:,stat.time>1.2 & stat.time<1.8).*...
+    stat.prob(:,:,stat.time>1.2 & stat.time<1.8),3)); 
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Mean BSR');
+
+%% plot using raincloud plot
+
+groupsizes=result.num_subj_lst;
+conditions=lv_evt_list;
+conds = {'hit'; 'miss'};
+condData = []; uData = [];
+for indGroup = 1
+    if indGroup == 1
+        relevantEntries = 1:groupsizes(1)*numel(conds);
+    elseif indGroup == 2
+        relevantEntries = groupsizes(1)*numel(conds)+1:...
+             groupsizes(1)*numel(conds)+groupsizes(2)*numel(conds);
+    end
+    for indCond = 1:numel(conds)
+        targetEntries = relevantEntries(conditions(relevantEntries)==indCond);        
+        condData{indGroup}(indCond,:) = -1.*result.vsc(targetEntries,indLV);
+        uData{indGroup}(indCond,:) = -1.*result.usc(targetEntries,indLV);
+    end
+end
+
+cBrew(1,:) = 2.*[.3 .1 .1];
+cBrew(2,:) = [.6 .6 .6];
+
+idx_outlier = cell(1); idx_standard = cell(1);
+for indGroup = 1
+    dataToPlot = uData{indGroup}';
+    % define outlier as lin. modulation that is more than three scaled median absolute deviations (MAD) away from the median
+    X = [1 1; 1 2]; b=X\dataToPlot'; IndividualSlopes = b(2,:);
+    outliers = isoutlier(IndividualSlopes, 'median');
+    idx_outlier{indGroup} = find(outliers);
+    idx_standard{indGroup} = find(outliers==0);
+end
+
+h = figure('units','centimeter','position',[0 0 25 10]);
+for indGroup = 1
+    dataToPlot = uData{indGroup}';
+    % read into cell array of the appropriate dimensions
+    data = []; data_ws = [];
+    for i = 1:2
+        for j = 1:1
+            data{i, j} = dataToPlot(:,i);
+            % individually demean for within-subject visualization
+            data_ws{i, j} = dataToPlot(:,i)-...
+                nanmean(dataToPlot(:,:),2)+...
+                repmat(nanmean(nanmean(dataToPlot(:,:),2),1),size(dataToPlot(:,:),1),1);
+            data_nooutlier{i, j} = data{i, j};
+            data_nooutlier{i, j}(idx_outlier{indGroup}) = [];
+            data_ws_nooutlier{i, j} = data_ws{i, j};
+            data_ws_nooutlier{i, j}(idx_outlier{indGroup}) = [];
+            % sort outliers to back in original data for improved plot overlap
+            data_ws{i, j} = [data_ws{i, j}(idx_standard{indGroup}); data_ws{i, j}(idx_outlier{indGroup})];
+        end
+    end
+
+    % IMPORTANT: plot individually centered estimates, stats on uncentered estimates!
+
+    subplot(1,2,indGroup);
+    set(gcf,'renderer','Painters')
+        cla;
+        cl = cBrew(indGroup,:);
+        [~, dist] = rm_raincloud_fixedSpacing(data_ws, [.8 .8 .8],1);
+        h_rc = rm_raincloud_fixedSpacing(data_ws_nooutlier, cl,1,[],[],[],dist);
+        view([90 -90]);
+        axis ij
+    box(gca,'off')
+    set(gca, 'YTickLabels', {conds{2}; conds{1}});
+    yticks = get(gca, 'ytick'); ylim([yticks(1)-(yticks(2)-yticks(1))./2, yticks(2)+(yticks(2)-yticks(1))./1.5]);
+
+    minmax = [min(min(cat(2,data_ws{:}))), max(max(cat(2,data_ws{:})))];
+    xlim(minmax+[-0.2*diff(minmax), 0.2*diff(minmax)])
+    ylabel('scene category'); xlabel({'Brainscore'; '[Individually centered]'})
+
+    % test linear effect
+    curData = [data_nooutlier{1, 1}, data_nooutlier{2, 1}];
+    X = [1 1; 1 2]; b=X\curData'; IndividualSlopes = b(2,:);
+    [~, p, ci, stats] = ttest(IndividualSlopes);
+    title(['M:', num2str(round(mean(IndividualSlopes),3)), '; p=', num2str(round(p,3))])
+end

code/b3b_taskPLS_recognition_erf.m

@@ -0,0 +1,107 @@
+% Set up EEG PLS for a task PLS using spectral power
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'erp');
+pn.data_erf = fullfile(rootpath, 'data', 'erf');
+pn.data = fullfile(rootpath, 'data', 'stats'); mkdir(pn.data);
+pn.tools = fullfile(rootpath, 'tools');
+    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'shadedErrorBar'));
+
+%% add seed for reproducibility
+
+rng(0, 'twister');
+    
+%% load event info
+
+load(fullfile(pn.data_eeg, ['sub-001_task-xxxx_eeg_art.mat']), 'events');
+
+parameter = {'scene_category'; 'old'; 'behavior'; 'subsequent_memory'};
+for ind_param = 1:numel(parameter)
+    conds.(parameter{ind_param}) = unique(events.(parameter{ind_param}));
+end
+
+%% load erp
+
+for ind_id = 1:33
+    id = sprintf('sub-%03d', ind_id); disp(id)
+    load(fullfile(pn.data_erf, [id,'_erf.mat']));
+    for ind_option = 1:4%numel(conds.behavior)
+        if ind_id == 1
+             erpgroup.behavior.(conds.behavior{ind_option}) = erf.behavior{ind_option};
+             erpgroup.behavior.(conds.behavior{ind_option}) = ...
+                 rmfield(erpgroup.behavior.(conds.behavior{ind_option}), {'powspctrm'});
+             erpgroup.behavior.(conds.behavior{ind_option}).dimord = 'sub_chan_freq_time';
+        end
+        erpgroup.behavior.(conds.behavior{ind_option}).avg(ind_id,:,:,:) = erf.behavior{ind_option}.powspctrm;
+    end
+end
+
+time = erpgroup.behavior.hit.time;
+freq = erpgroup.behavior.hit.freq;
+channels = erpgroup.behavior.hit.label;
+
+mergeddata = cat(5, erpgroup.behavior.hit.avg, ...
+    erpgroup.behavior.miss.avg);
+
+num_chans = numel(channels);
+num_freqs = numel(freq);
+num_time = numel(time);
+
+num_subj_lst = [33];
+num_cond = 2;
+num_grp = 1;
+
+datamat_lst = cell(num_grp); lv_evt_list = [];
+indCount_cont = 1;
+indCount = 1;
+indGroup = 1;
+for indCond = 1:num_cond
+    for indID = 1:num_subj_lst(indGroup)
+        datamat_lst{indGroup}(indCount,:) = reshape(squeeze(mergeddata(indID,:,:,:,indCond)), [], 1);
+        lv_evt_list(indCount_cont) = indCond;
+        indCount = indCount+1;
+        indCount_cont = indCount_cont+1;
+    end
+end
+datamat_lst{indGroup}(isnan(datamat_lst{indGroup})) = 0;
+
+%% set PLS options and run PLS
+
+option = [];
+option.method = 1; % [1] | 2 | 3 | 4 | 5 | 6
+option.num_perm = 1000; %( single non-negative integer )
+option.num_split = 0; %( single non-negative integer )
+option.num_boot = 1000; % ( single non-negative integer )
+option.cormode = 0; % [0] | 2 | 4 | 6
+option.meancentering_type = 0;% [0] | 1 | 2 | 3
+option.boot_type = 'strat'; %['strat'] | 'nonstrat'
+
+result = pls_analysis(datamat_lst, num_subj_lst, num_cond, option);
+
+%% rearrange into fieldtrip structure
+
+lvdat = reshape(result.boot_result.compare_u(:,1), num_chans, num_freqs, num_time);
+%udat = reshape(result.u, num_chans, num_freqs, num_time);
+
+stat = [];
+stat.prob = lvdat;
+stat.dimord = 'chan_freq_time';
+stat.clusters = [];
+stat.clusters.prob = result.perm_result.sprob; % check for significance of LV
+stat.mask = lvdat > 3 | lvdat < -3;
+stat.cfg = option;
+stat.freq = freq;
+stat.time = time;
+
+save(fullfile(pn.data, 'b3b_taskpls_erf.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')

code/b3b_taskPLS_recognition_erf_LV1.m

@@ -0,0 +1,281 @@
+% Create an overview plot featuring the results of the multivariate PLS
+% comparing spectral changes during the stimulus period under load
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data         = fullfile(rootpath, 'data', 'stats');
+pn.figures      = fullfile(rootpath, 'figures');
+pn.tools        = fullfile(rootpath, 'tools');
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'fieldtrip')); ft_defaults;
+    addpath(genpath(fullfile(pn.tools, 'RainCloudPlots')));
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'winsor'));
+
+% set custom colormap
+cBrew = brewermap(500,'RdBu');
+cBrew = flipud(cBrew);
+colormap(cBrew)
+
+load(fullfile(pn.data, 'b3b_taskpls_erf.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')
+load(fullfile(rootpath, 'data','erp', ['sub-001_erp.mat']));
+elec = erp.scene_category{1}.elec;
+
+result.perm_result.sprob
+
+indLV = 1;
+
+lvdat = reshape(result.boot_result.compare_u(:,indLV), num_chans, num_freqs, num_time);
+stat.prob = lvdat;
+stat.mask = lvdat > 3 | lvdat < -3;
+
+% maskNaN = double(stat.mask);
+% maskNaN(maskNaN==0) = NaN;
+
+%% invert solution
+
+% stat.mask = stat.mask;
+% stat.prob = stat.prob.*-1;
+% result.vsc = result.vsc.*-1;
+% result.usc = result.usc.*-1;
+
+h = figure('units','centimeter','position',[0 0 15 10]);
+set(gcf,'renderer','Painters')
+statsPlot = [];
+%statsPlot = cat(1, statsPlot,squeeze(nanmean(stat.prob(1:64,:,:),1)));
+statsPlot = cat(1, statsPlot,squeeze(nanmean(stat.mask(1:64,:,:).*stat.prob(1:64,:,:),1)));
+range_plot = [-4 1.5];
+imagesc(stat.time,[],statsPlot, range_plot);
+set(gca,'Ydir','Normal');
+% set custom colormap
+ncol = 2000; cBrew = brewermap(ncol,'RdBu'); cBrew = flipud(cBrew);
+% adjust for unequal range
+reldev = min(abs(range_plot))./max(abs(range_plot));
+if abs(range_plot(2))<max(abs(range_plot))
+    nadjust = [1:ncol/2, ncol/2+round(linspace(1, ncol/2, reldev*(ncol/2)))];
+elseif abs(range_plot(1))<max(abs(range_plot))
+    nadjust = [round(linspace(1, ncol/2, reldev*(ncol/2))),ncol/2:ncol];
+else
+    nadjust = 1:ncol;
+end
+cBrew = cBrew(nadjust,:);
+colormap(cBrew); colorbar;
+xlabel('Time [s from stim onset]'); ylabel('Frequency (Hz)');
+title({'ERF changes'; ['p = ', num2str(round(result.perm_result.sprob(indLV),4))]})
+set(findall(gcf,'-property','FontSize'),'FontSize',18)
+set(gca, 'YTickLabels', round(stat.freq(get(gca, 'YTick')),0));
+
+%https://de.mathworks.com/matlabcentral/answers/476715-superimposing-two-imagesc-graphs-over-each-other
+
+% h = figure('units','centimeter','position',[0 0 15 10]);
+% set(gcf,'renderer','Painters')
+% hold all; axis square; 
+% ax1 = axes; 
+%     statsPlot = [];
+%     statsPlot = cat(1, statsPlot,squeeze(nanmean(stat.prob(1:64,:,:),1)));
+%     im1 = imagesc(stat.time,[],statsPlot, [-4 4]);
+%     im1.AlphaData = .5;
+%     axis square; 
+% ax2 = axes; 
+%     statsPlot = [];
+%     statsPlot = cat(1, statsPlot,squeeze(nanmean(stat.mask(1:64,:,:).*stat.prob(1:64,:,:),1)));
+%     im2 = imagesc(stat.time,[],statsPlot, [-4 4]);
+%     im2.AlphaData = .8;
+%     axis square; 
+% linkaxes([ax1,ax2])
+% %%Hide the top axes 
+% ax2.Visible = 'off'; 
+% ax2.XTick = []; 
+% ax2.YTick = []; 
+% colormap(ax1,gray) 
+% colormap(ax2,cBrew)
+% set(gca,'Ydir','Normal');
+% xlabel('Time [s from stim onset]'); ylabel('Frequency (Hz)');
+% title({'ERF changes'; ['p = ', num2str(round(result.perm_result.sprob(indLV),4))]})
+% set(findall(gcf,'-property','FontSize'),'FontSize',18)
+% set(gca, 'YTickLabels', round(stat.freq(get(gca, 'YTick')),0));
+% colorbar;
+
+
+% figureName = ['b01_pls_traces'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');
+
+%% plot multivariate topographies
+
+% set custom colormap
+cBrew = brewermap(500,'RdBu');
+cBrew = flipud(cBrew);
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+cfg.zlim = [-4 4]; 
+cfg.figure = h;
+plotData.powspctrm = squeeze(nanmean(nanmean(stat.mask(:,stat.freq>8 & stat.freq<20,stat.time>1 & stat.time<2).*...
+    stat.prob(:,stat.freq>8 & stat.freq<20,stat.time>1 & stat.time<2),3),2)); 
+[~, sortind] = sort(plotData.powspctrm, 'ascend');
+cfg.marker = 'off'; 
+cfg.highlight = 'yes';
+cfg.highlightchannel = plotData.label(sortind(1:6));
+cfg.highlightcolor = [1 1 1];
+cfg.highlightsymbol = '.';
+cfg.highlightsize = 30;
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Mean BSR');
+title('AlphaBeta 1-2s')
+
+%% 'Theta 1-2s'
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+cfg.zlim = [-2 2]; 
+cfg.figure = h;
+plotData.powspctrm = squeeze(nanmean(nanmean(stat.mask(:,stat.freq>2 &stat.freq<7,stat.time>1 & stat.time<2).*...
+    stat.prob(:,stat.freq>2 & stat.freq<7,stat.time>1 & stat.time<2),3),2)); 
+[~, sortind] = sort(plotData.powspctrm, 'descend');
+cfg.marker = 'off'; 
+cfg.highlight = 'yes';
+cfg.highlightchannel = plotData.label(sortind(1:3));
+cfg.highlightcolor = [1 1 1];
+cfg.highlightsymbol = '.';
+cfg.highlightsize = 30;
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Mean BSR');
+title('Theta 1-2s')
+
+%% title('AlphaBeta 0.5-1s')
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+cfg.zlim = [-2 2]; 
+cfg.figure = h;
+plotData.powspctrm = squeeze(nanmean(nanmean(stat.mask(:,stat.freq>8 & stat.freq<20,stat.time>0.5 & stat.time<1).*...
+    stat.prob(:,stat.freq>8 & stat.freq<20,stat.time>0.5 & stat.time<1),3),2)); 
+[~, sortind] = sort(plotData.powspctrm, 'descend');
+cfg.marker = 'off'; 
+cfg.highlight = 'yes';
+cfg.highlightchannel = plotData.label(sortind(1:3));
+cfg.highlightcolor = [1 1 1];
+cfg.highlightsymbol = '.';
+cfg.highlightsize = 30;
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Mean BSR');
+title('AlphaBeta 0.5-1s')
+
+%% plot using raincloud plot
+
+groupsizes=result.num_subj_lst;
+conditions=lv_evt_list;
+conds = {'hit'; 'miss'};
+condData = []; uData = [];
+for indGroup = 1
+    if indGroup == 1
+        relevantEntries = 1:groupsizes(1)*numel(conds);
+    elseif indGroup == 2
+        relevantEntries = groupsizes(1)*numel(conds)+1:...
+             groupsizes(1)*numel(conds)+groupsizes(2)*numel(conds);
+    end
+    for indCond = 1:numel(conds)
+        targetEntries = relevantEntries(conditions(relevantEntries)==indCond);        
+        condData{indGroup}(indCond,:) = result.vsc(targetEntries,indLV);
+        uData{indGroup}(indCond,:) = result.usc(targetEntries,indLV);
+    end
+end
+
+cBrew(1,:) = 2.*[.3 .1 .1];
+cBrew(2,:) = [.6 .6 .6];
+
+idx_outlier = cell(1); idx_standard = cell(1);
+for indGroup = 1
+    dataToPlot = uData{indGroup}';
+    % define outlier as lin. modulation that is more than three scaled median absolute deviations (MAD) away from the median
+    X = [1 1; 1 2]; b=X\dataToPlot'; IndividualSlopes = b(2,:);
+    outliers = isoutlier(IndividualSlopes, 'median');
+    idx_outlier{indGroup} = find(outliers);
+    idx_standard{indGroup} = find(outliers==0);
+end
+
+h = figure('units','centimeter','position',[0 0 8 10]);
+for indGroup = 1
+    dataToPlot = uData{indGroup}';
+    % read into cell array of the appropriate dimensions
+    data = []; data_ws = [];
+    for i = 1:2
+        for j = 1:1
+            data{i, j} = dataToPlot(:,i);
+            % individually demean for within-subject visualization
+            data_ws{i, j} = dataToPlot(:,i)-...
+                nanmean(dataToPlot(:,:),2)+...
+                repmat(nanmean(nanmean(dataToPlot(:,:),2),1),size(dataToPlot(:,:),1),1);
+            data_nooutlier{i, j} = data{i, j};
+            data_nooutlier{i, j}(idx_outlier{indGroup}) = [];
+            data_ws_nooutlier{i, j} = data_ws{i, j};
+            data_ws_nooutlier{i, j}(idx_outlier{indGroup}) = [];
+            % sort outliers to back in original data for improved plot overlap
+            data_ws{i, j} = [data_ws{i, j}(idx_standard{indGroup}); data_ws{i, j}(idx_outlier{indGroup})];
+        end
+    end
+
+    % IMPORTANT: plot individually centered estimates, stats on uncentered estimates!
+
+    %subplot(1,2,indGroup);
+    set(gcf,'renderer','Painters')
+        cla;
+        cl = cBrew(indGroup,:);
+        [~, dist] = rm_raincloud_fixedSpacing(data_ws, [.8 .8 .8],1);
+        h_rc = rm_raincloud_fixedSpacing(data_ws_nooutlier, cl,1,[],[],[],dist);
+        view([90 -90]);
+        axis ij
+    box(gca,'off')
+    set(gca, 'YTickLabels', {conds{2}; conds{1}});
+    yticks = get(gca, 'ytick'); ylim([yticks(1)-(yticks(2)-yticks(1))./2, yticks(2)+(yticks(2)-yticks(1))./1.5]);
+
+    minmax = [min(min(cat(2,data_ws{:}))), max(max(cat(2,data_ws{:})))];
+    xlim(minmax+[-0.2*diff(minmax), 0.2*diff(minmax)])
+    ylabel('scene category'); xlabel({'Brainscore'; '[Individually centered]'})
+
+    % test linear effect
+    curData = [data_nooutlier{1, 1}, data_nooutlier{2, 1}];
+    X = [1 1; 1 2]; b=X\curData'; IndividualSlopes = b(2,:);
+    [~, p, ci, stats] = ttest(IndividualSlopes);
+    title(['M:', num2str(round(mean(IndividualSlopes),3)), '; p=', num2str(round(p,3))])
+end

code/b3x_rcognition_plot_theta_alpha.m

@@ -0,0 +1,146 @@
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'erp');
+pn.data_erf = fullfile(rootpath, 'data', 'erf');
+pn.data = fullfile(rootpath, 'data', 'stats'); mkdir(pn.data);
+pn.tools = fullfile(rootpath, 'tools');
+    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'shadedErrorBar'));
+
+%% add seed for reproducibility
+
+rng(0, 'twister');
+
+%% load erp
+
+for ind_id = 1:33
+    id = sprintf('sub-%03d', ind_id); disp(id)
+    load(fullfile(pn.data_erf, [id,'_erf.mat']));
+    for ind_option = 1:numel(conds.old)
+        if ind_id == 1
+             erpgroup.behavior.(conds.behavior{ind_option}) = erf.behavior{ind_option};
+             erpgroup.behavior.(conds.behavior{ind_option}) = ...
+                 rmfield(erpgroup.behavior.(conds.behavior{ind_option}), {'powspctrm'});
+             erpgroup.behavior.(conds.behavior{ind_option}).dimord = 'sub_chan_freq_time';
+             freq = erpgroup.behavior.(conds.behavior{ind_option}).freq;
+        end
+        idx_freq_t = freq >2 & freq <7;
+        thetagroup.behavior.(conds.behavior{ind_option}).avg(ind_id,:,:,:) = squeeze(nanmean(erf.behavior{ind_option}.powspctrm(:,idx_freq_t,:),2));
+        idx_freq_a = freq >8 & freq <20;
+        alphagroup.behavior.(conds.behavior{ind_option}).avg(ind_id,:,:,:) = squeeze(nanmean(erf.behavior{ind_option}.powspctrm(:,idx_freq_a,:),2));
+    end
+end
+
+time = erpgroup.behavior.correctreject.time;
+channels = erpgroup.behavior.correctreject.label;
+
+%% get max. channels
+
+load(fullfile(pn.data, 'b3b_taskpls_erf.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')
+
+plotData.powspctrm = squeeze(nanmean(nanmean(stat.mask(:,stat.freq>8 & stat.freq<20,stat.time>1 & stat.time<2).*...
+    stat.prob(:,stat.freq>8 & stat.freq<20,stat.time>1 & stat.time<2),3),2));
+[~, alphaChanMin] = sort(plotData.powspctrm, 'ascend');
+
+plotData.powspctrm = squeeze(nanmean(nanmean(stat.mask(:,stat.freq>2 &stat.freq<7,stat.time>1 & stat.time<2).*...
+    stat.prob(:,stat.freq>2 & stat.freq<7,stat.time>1 & stat.time<2),3),2)); 
+[~, thetaChan] = sort(plotData.powspctrm, 'descend');
+
+plotData.powspctrm = squeeze(nanmean(nanmean(stat.mask(:,stat.freq>8 & stat.freq<20,stat.time>0.5 & stat.time<1).*...
+    stat.prob(:,stat.freq>8 & stat.freq<20,stat.time>0.5 & stat.time<1),3),2)); 
+[~, alphaChanMax] = sort(plotData.powspctrm, 'descend');
+
+%% visualize for frontal theta
+
+mergeddata = cat(4, thetagroup.behavior.correctreject.avg, thetagroup.behavior.falsealarm.avg);
+
+%figure; imagesc(squeeze(nanmean(thetagroup.old.new.avg-thetagroup.old.old.avg,1)))
+
+idx_chans = thetaChan(1:3); channels(idx_chans) %[33,38,47];
+condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
+
+h = figure('units','centimeters','position',[0 0 10 8]);
+cla; hold on;
+% new value = old value ? subject average + grand average
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
+    {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
+    {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
+xlabel('Time (s) from stim onset')
+legend([l1.mainLine, l2.mainLine],{'correctreject', 'falsealarm'}, ...
+    'location', 'southwest'); legend('boxoff')
+xlim([-.5 1.7]); %ylim(YLim)
+ylabel({'theta power';'(log10)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',12)
+
+%% visualize for posterior alpha
+
+mergeddata = cat(4, alphagroup.behavior.correctreject.avg, alphagroup.behavior.falsealarm.avg);
+
+idx_chans = alphaChanMin(1:3); channels(idx_chans) %[33,38,47];
+condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
+
+h = figure('units','centimeters','position',[0 0 10 8]);
+cla; hold on;
+% new value = old value ? subject average + grand average
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
+    {'color', 'k','linewidth', 3}, 'patchSaturation', .1);
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
+    {'color', 'r','linewidth', 3}, 'patchSaturation', .1);
+xlabel('Time (s) from stim onset')
+legend([l1.mainLine, l2.mainLine],{'correctreject', 'falsealarm'}, ...
+    'location', 'northeast'); legend('boxoff')
+xlim([-.5 1.75]);
+ylabel({'alpha power';'(log10)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',15)
+
+%% visualize for central alpha
+
+mergeddata = cat(4, alphagroup.behavior.correctreject.avg, alphagroup.behavior.falsealarm.avg);
+
+idx_chans = alphaChanMax(1:3); channels(idx_chans) %[33,38,47];
+condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
+
+h = figure('units','centimeters','position',[0 0 10 8]);
+cla; hold on;
+% new value = old value ? subject average + grand average
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
+    {'color', 'k','linewidth', 3}, 'patchSaturation', .1);
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
+    {'color', 'r','linewidth', 3}, 'patchSaturation', .1);
+xlabel('Time (s) from stim onset')
+legend([l1.mainLine, l2.mainLine],{'correctreject', 'falsealarm'}, ...
+    'location', 'northeast'); legend('boxoff')
+xlim([-.5 1.75]);
+ylabel({'alpha power';'(log10)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',15)

code/b4a_taskPLS_memory_erp.m

@@ -0,0 +1,106 @@
+% Set up EEG PLS for a task PLS using spectral power
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'erp');
+pn.data_erf = fullfile(rootpath, 'data', 'erf');
+pn.data = fullfile(rootpath, 'data', 'stats'); mkdir(pn.data);
+pn.tools = fullfile(rootpath, 'tools');
+    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'shadedErrorBar'));
+
+%% add seed for reproducibility
+
+rng(0, 'twister');
+    
+%% load event info
+
+load(fullfile(pn.data_eeg, ['sub-001_task-xxxx_eeg_art.mat']), 'events');
+
+parameter = {'scene_category'; 'old'; 'behavior'; 'subsequent_memory'};
+for ind_param = 1:numel(parameter)
+    conds.(parameter{ind_param}) = unique(events.(parameter{ind_param}));
+end
+
+%% load erp_bl
+
+for ind_id = 1:33
+    id = sprintf('sub-%03d', ind_id); disp(id)
+    load(fullfile(pn.data_erp, [id,'_erp_bl.mat']));
+    for ind_option = 2:3
+        if ind_id == 1
+             erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}) = erp_bl.subsequent_memory{ind_option};
+             erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}) = ...
+                 rmfield(erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}), {'avg', 'var', 'dof'});
+             erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}).dimord = 'sub_chan_time';
+        end
+        erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}).avg(ind_id,:,:) = erp_bl.subsequent_memory{ind_option}.avg;
+    end
+end
+
+time = erpgroup.subsequent_memory.subsequent_remembered.time;
+elec = erpgroup.subsequent_memory.subsequent_remembered.elec;
+channels = erpgroup.subsequent_memory.subsequent_remembered.label;
+
+mergeddata = cat(4, erpgroup.subsequent_memory.subsequent_forgotten.avg, ...
+    erpgroup.subsequent_memory.subsequent_remembered.avg);
+
+num_chans = numel(channels);
+num_freqs = 1;
+num_time = numel(time);
+
+num_subj_lst = [33];
+num_cond = 2;
+num_grp = 1;
+
+datamat_lst = cell(num_grp); lv_evt_list = [];
+indCount_cont = 1;
+indCount = 1;
+indGroup = 1;
+for indCond = 1:num_cond
+    for indID = 1:num_subj_lst(indGroup)
+        datamat_lst{indGroup}(indCount,:) = reshape(squeeze(mergeddata(indID,:,:,indCond)), [], 1);
+        lv_evt_list(indCount_cont) = indCond;
+        indCount = indCount+1;
+        indCount_cont = indCount_cont+1;
+    end
+end
+datamat_lst{indGroup}(isnan(datamat_lst{indGroup})) = 0;
+
+%% set PLS options and run PLS
+
+option = [];
+option.method = 1; % [1] | 2 | 3 | 4 | 5 | 6
+option.num_perm = 1000; %( single non-negative integer )
+option.num_split = 0; %( single non-negative integer )
+option.num_boot = 1000; % ( single non-negative integer )
+option.cormode = 0; % [0] | 2 | 4 | 6
+option.meancentering_type = 0;% [0] | 1 | 2 | 3
+option.boot_type = 'strat'; %['strat'] | 'nonstrat'
+
+result = pls_analysis(datamat_lst, num_subj_lst, num_cond, option);
+
+%% rearrange into fieldtrip structure
+
+lvdat = reshape(result.boot_result.compare_u(:,1), num_chans, num_freqs, num_time);
+%udat = reshape(result.u, num_chans, num_freqs, num_time);
+
+stat = [];
+stat.prob = lvdat;
+stat.dimord = 'chan_freq_time';
+stat.clusters = [];
+stat.clusters.prob = result.perm_result.sprob; % check for significance of LV
+stat.mask = lvdat > 3 | lvdat < -3;
+stat.cfg = option;
+stat.time = time;
+
+save(fullfile(pn.data, 'e01_taskpls_erp.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')

code/b4a_taskPLS_memory_erp_LV1.m

@@ -0,0 +1,192 @@
+% Create an overview plot featuring the results of the multivariate PLS
+% comparing spectral changes during the stimulus period under load
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data         = fullfile(rootpath, 'data', 'stats');
+pn.figures      = fullfile(rootpath, 'figures');
+pn.tools        = fullfile(rootpath, 'tools');
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'fieldtrip')); ft_defaults;
+    addpath(genpath(fullfile(pn.tools, 'RainCloudPlots')));
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'winsor'));
+
+% set custom colormap
+cBrew = brewermap(500,'RdBu');
+cBrew = flipud(cBrew);
+colormap(cBrew)
+
+load(fullfile(pn.data, 'e01_taskpls_erp.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')
+load(fullfile(rootpath, 'data','erp', ['sub-001_erp.mat']));
+elec = erp.scene_category{1}.elec;
+
+result.perm_result.sprob
+
+indLV = 1;
+
+lvdat = reshape(result.boot_result.compare_u(:,indLV), num_chans, num_freqs, num_time);
+stat.prob = lvdat;
+stat.mask = lvdat > 3 | lvdat < -3;
+
+% maskNaN = double(stat.mask);
+% maskNaN(maskNaN==0) = NaN;
+
+%% invert solution
+
+% stat.mask = stat.mask;
+% stat.prob = stat.prob.*-1;
+% result.vsc = result.vsc.*-1;
+% result.usc = result.usc.*-1;
+
+h = figure('units','centimeter','position',[0 0 15 10]);
+set(gcf,'renderer','Painters')
+hold on;
+statsPlot = [];
+statsPlot = cat(1, statsPlot,squeeze(nanmean(stat.mask(1:64,:,:).*stat.prob(1:64,:,:),1)));
+plot(stat.time,statsPlot, 'k')
+xlabel('Time [s from stim onset]'); ylabel('max abs BSR');
+title({'ERP changes'; ['p = ', num2str(round(result.perm_result.sprob(indLV),4))]})
+set(findall(gcf,'-property','FontSize'),'FontSize',18)
+
+% h = figure('units','centimeter','position',[0 0 15 10]);
+% set(gcf,'renderer','Painters')
+% statsPlot = [];
+% statsPlot = cat(1, statsPlot,squeeze(nanmax(abs(stat.prob(1:64,:,:)),[],1)));
+% plot(stat.time,statsPlot, 'k')
+% xlabel('Time [s from stim onset]'); ylabel('max abs BSR');
+% title({'ERP changes'; ['p = ', num2str(round(result.perm_result.sprob(indLV),4))]})
+% set(findall(gcf,'-property','FontSize'),'FontSize',18)
+% % figureName = ['b01_pls_traces'];
+% % saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% % saveas(h, fullfile(pn.figures, figureName), 'png');
+
+%% plot multivariate topographies
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+cfg.zlim = [-2 2]; 
+cfg.figure = h;
+plotData.powspctrm = squeeze(nanmean(stat.mask(:,:,stat.time>.3 & stat.time<.5).*...
+    stat.prob(:,:,stat.time>.3 & stat.time<.5),3)); 
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Mean BSR');
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+cfg.zlim = [-.5 .5]; 
+cfg.figure = h;
+plotData.powspctrm = squeeze(nanmean(stat.mask(:,:,:).*...
+    stat.prob(:,:,:),3)); 
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Mean BSR');
+
+
+%% plot using raincloud plot
+
+groupsizes=result.num_subj_lst;
+conditions=lv_evt_list;
+conds = {'forgotten'; 'remembered'};
+condData = []; uData = [];
+for indGroup = 1
+    if indGroup == 1
+        relevantEntries = 1:groupsizes(1)*numel(conds);
+    elseif indGroup == 2
+        relevantEntries = groupsizes(1)*numel(conds)+1:...
+             groupsizes(1)*numel(conds)+groupsizes(2)*numel(conds);
+    end
+    for indCond = 1:numel(conds)
+        targetEntries = relevantEntries(conditions(relevantEntries)==indCond);        
+        condData{indGroup}(indCond,:) = result.vsc(targetEntries,indLV);
+        uData{indGroup}(indCond,:) = result.usc(targetEntries,indLV);
+    end
+end
+
+cBrew(1,:) = 2.*[.3 .1 .1];
+cBrew(2,:) = [.6 .6 .6];
+
+idx_outlier = cell(1); idx_standard = cell(1);
+for indGroup = 1
+    dataToPlot = uData{indGroup}';
+    % define outlier as lin. modulation that is more than three scaled median absolute deviations (MAD) away from the median
+    X = [1 1; 1 2]; b=X\dataToPlot'; IndividualSlopes = b(2,:);
+    outliers = isoutlier(IndividualSlopes, 'median');
+    idx_outlier{indGroup} = find(outliers);
+    idx_standard{indGroup} = find(outliers==0);
+end
+
+h = figure('units','centimeter','position',[0 0 8 10]);
+for indGroup = 1
+    dataToPlot = uData{indGroup}';
+    % read into cell array of the appropriate dimensions
+    data = []; data_ws = [];
+    for i = 1:2
+        for j = 1:1
+            data{i, j} = dataToPlot(:,i);
+            % individually demean for within-subject visualization
+            data_ws{i, j} = dataToPlot(:,i)-...
+                nanmean(dataToPlot(:,:),2)+...
+                repmat(nanmean(nanmean(dataToPlot(:,:),2),1),size(dataToPlot(:,:),1),1);
+            data_nooutlier{i, j} = data{i, j};
+            data_nooutlier{i, j}(idx_outlier{indGroup}) = [];
+            data_ws_nooutlier{i, j} = data_ws{i, j};
+            data_ws_nooutlier{i, j}(idx_outlier{indGroup}) = [];
+            % sort outliers to back in original data for improved plot overlap
+            data_ws{i, j} = [data_ws{i, j}(idx_standard{indGroup}); data_ws{i, j}(idx_outlier{indGroup})];
+        end
+    end
+
+    % IMPORTANT: plot individually centered estimates, stats on uncentered estimates!
+
+    set(gcf,'renderer','Painters')
+        cla;
+        cl = cBrew(indGroup,:);
+        [~, dist] = rm_raincloud_fixedSpacing(data_ws, [.8 .8 .8],1);
+        h_rc = rm_raincloud_fixedSpacing(data_ws_nooutlier, cl,1,[],[],[],dist);
+        view([90 -90]);
+        axis ij
+    box(gca,'off')
+    set(gca, 'YTickLabels', {conds{2}; conds{1}});
+    yticks = get(gca, 'ytick'); ylim([yticks(1)-(yticks(2)-yticks(1))./2, yticks(2)+(yticks(2)-yticks(1))./1.5]);
+
+    minmax = [min(min(cat(2,data_ws{:}))), max(max(cat(2,data_ws{:})))];
+    xlim(minmax+[-0.2*diff(minmax), 0.2*diff(minmax)])
+    ylabel('scene category'); xlabel({'Brainscore'; '[Individually centered]'})
+
+    % test linear effect
+    curData = [data_nooutlier{1, 1}, data_nooutlier{2, 1}];
+    X = [1 1; 1 2]; b=X\curData'; IndividualSlopes = b(2,:);
+    [~, p, ci, stats] = ttest(IndividualSlopes);
+    %title(['M:', num2str(round(mean(IndividualSlopes),3)), '; p=', num2str(round(p,3))])
+    title(['M:', num2str(round(mean(IndividualSlopes),3)), '; p=', num2str(round(result.perm_result.sprob(indLV),3))])
+end
+
+set(findall(gcf,'-property','FontSize'),'FontSize',14)

code/b4a_taskPLS_memory_erp_nooutlier.m

@@ -0,0 +1,109 @@
+% Set up EEG PLS for a task PLS using spectral power
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'erp');
+pn.data_erf = fullfile(rootpath, 'data', 'erf');
+pn.data = fullfile(rootpath, 'data', 'stats'); mkdir(pn.data);
+pn.tools = fullfile(rootpath, 'tools');
+    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'shadedErrorBar'));
+
+%% add seed for reproducibility
+
+rng(0, 'twister');
+    
+%% load event info
+
+load(fullfile(pn.data_eeg, ['sub-001_task-xxxx_eeg_art.mat']), 'events');
+
+parameter = {'scene_category'; 'old'; 'behavior'; 'subsequent_memory'};
+for ind_param = 1:numel(parameter)
+    conds.(parameter{ind_param}) = unique(events.(parameter{ind_param}));
+end
+
+%% load erp_bl
+
+for ind_id = 1:33
+    id = sprintf('sub-%03d', ind_id); disp(id)
+    load(fullfile(pn.data_erp, [id,'_erp_bl.mat']));
+    for ind_option = 2:3
+        if ind_id == 1
+             erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}) = erp_bl.subsequent_memory{ind_option};
+             erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}) = ...
+                 rmfield(erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}), {'avg', 'var', 'dof'});
+             erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}).dimord = 'sub_chan_time';
+        end
+        erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}).avg(ind_id,:,:) = erp_bl.subsequent_memory{ind_option}.avg;
+    end
+end
+
+time = erpgroup.subsequent_memory.subsequent_remembered.time;
+elec = erpgroup.subsequent_memory.subsequent_remembered.elec;
+channels = erpgroup.subsequent_memory.subsequent_remembered.label;
+
+mergeddata = cat(4, erpgroup.subsequent_memory.subsequent_forgotten.avg, ...
+    erpgroup.subsequent_memory.subsequent_remembered.avg);
+
+% remove outlier
+mergeddata(1,:,:,:) = [];
+
+num_chans = numel(channels);
+num_freqs = 1;
+num_time = numel(time);
+
+num_subj_lst = [32];
+num_cond = 2;
+num_grp = 1;
+
+datamat_lst = cell(num_grp); lv_evt_list = [];
+indCount_cont = 1;
+indCount = 1;
+indGroup = 1;
+for indCond = 1:num_cond
+    for indID = 1:num_subj_lst(indGroup)
+        datamat_lst{indGroup}(indCount,:) = reshape(squeeze(mergeddata(indID,:,:,indCond)), [], 1);
+        lv_evt_list(indCount_cont) = indCond;
+        indCount = indCount+1;
+        indCount_cont = indCount_cont+1;
+    end
+end
+datamat_lst{indGroup}(isnan(datamat_lst{indGroup})) = 0;
+
+%% set PLS options and run PLS
+
+option = [];
+option.method = 1; % [1] | 2 | 3 | 4 | 5 | 6
+option.num_perm = 1000; %( single non-negative integer )
+option.num_split = 0; %( single non-negative integer )
+option.num_boot = 1000; % ( single non-negative integer )
+option.cormode = 0; % [0] | 2 | 4 | 6
+option.meancentering_type = 0;% [0] | 1 | 2 | 3
+option.boot_type = 'strat'; %['strat'] | 'nonstrat'
+
+result = pls_analysis(datamat_lst, num_subj_lst, num_cond, option);
+
+%% rearrange into fieldtrip structure
+
+lvdat = reshape(result.boot_result.compare_u(:,1), num_chans, num_freqs, num_time);
+%udat = reshape(result.u, num_chans, num_freqs, num_time);
+
+stat = [];
+stat.prob = lvdat;
+stat.dimord = 'chan_freq_time';
+stat.clusters = [];
+stat.clusters.prob = result.perm_result.sprob; % check for significance of LV
+stat.mask = lvdat > 3 | lvdat < -3;
+stat.cfg = option;
+stat.time = time;
+
+save(fullfile(pn.data, 'e01_taskpls_erp_nooutlier.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')

code/b4b_taskPLS_memory_erf.m

@@ -0,0 +1,107 @@
+% Set up EEG PLS for a task PLS using spectral power
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'erp');
+pn.data_erf = fullfile(rootpath, 'data', 'erf');
+pn.data = fullfile(rootpath, 'data', 'stats'); mkdir(pn.data);
+pn.tools = fullfile(rootpath, 'tools');
+    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'shadedErrorBar'));
+
+%% add seed for reproducibility
+
+rng(0, 'twister');
+    
+%% load event info
+
+load(fullfile(pn.data_eeg, ['sub-001_task-xxxx_eeg_art.mat']), 'events');
+
+parameter = {'scene_category'; 'old'; 'behavior'; 'subsequent_memory'};
+for ind_param = 1:numel(parameter)
+    conds.(parameter{ind_param}) = unique(events.(parameter{ind_param}));
+end
+
+%% load erp
+
+for ind_id = 1:33
+    id = sprintf('sub-%03d', ind_id); disp(id)
+    load(fullfile(pn.data_erf, [id,'_erf.mat']));
+    for ind_option = 2:3
+        if ind_id == 1
+             erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}) = erf.subsequent_memory{ind_option};
+             erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}) = ...
+                 rmfield(erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}), {'powspctrm'});
+             erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}).dimord = 'sub_chan_freq_time';
+        end
+        erpgroup.subsequent_memory.(conds.subsequent_memory{ind_option}).avg(ind_id,:,:,:) = erf.subsequent_memory{ind_option}.powspctrm;
+    end
+end
+
+time = erpgroup.subsequent_memory.subsequent_remembered.time;
+freq = erpgroup.subsequent_memory.subsequent_remembered.freq;
+channels = erpgroup.subsequent_memory.subsequent_remembered.label;
+
+mergeddata = cat(5, erpgroup.subsequent_memory.subsequent_forgotten.avg, ...
+    erpgroup.subsequent_memory.subsequent_remembered.avg);
+
+num_chans = numel(channels);
+num_freqs = numel(freq);
+num_time = numel(time);
+
+num_subj_lst = [33];
+num_cond = 2;
+num_grp = 1;
+
+datamat_lst = cell(num_grp); lv_evt_list = [];
+indCount_cont = 1;
+indCount = 1;
+indGroup = 1;
+for indCond = 1:num_cond
+    for indID = 1:num_subj_lst(indGroup)
+        datamat_lst{indGroup}(indCount,:) = reshape(squeeze(mergeddata(indID,:,:,:,indCond)), [], 1);
+        lv_evt_list(indCount_cont) = indCond;
+        indCount = indCount+1;
+        indCount_cont = indCount_cont+1;
+    end
+end
+datamat_lst{indGroup}(isnan(datamat_lst{indGroup})) = 0;
+
+%% set PLS options and run PLS
+
+option = [];
+option.method = 1; % [1] | 2 | 3 | 4 | 5 | 6
+option.num_perm = 1000; %( single non-negative integer )
+option.num_split = 0; %( single non-negative integer )
+option.num_boot = 1000; % ( single non-negative integer )
+option.cormode = 0; % [0] | 2 | 4 | 6
+option.meancentering_type = 0;% [0] | 1 | 2 | 3
+option.boot_type = 'strat'; %['strat'] | 'nonstrat'
+
+result = pls_analysis(datamat_lst, num_subj_lst, num_cond, option);
+
+%% rearrange into fieldtrip structure
+
+lvdat = reshape(result.boot_result.compare_u(:,1), num_chans, num_freqs, num_time);
+%udat = reshape(result.u, num_chans, num_freqs, num_time);
+
+stat = [];
+stat.prob = lvdat;
+stat.dimord = 'chan_freq_time';
+stat.clusters = [];
+stat.clusters.prob = result.perm_result.sprob; % check for significance of LV
+stat.mask = lvdat > 3 | lvdat < -3;
+stat.cfg = option;
+stat.freq = freq;
+stat.time = time;
+
+save(fullfile(pn.data, 'b4b_taskpls_erf.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')

code/b4b_taskPLS_memory_erf_LV1.m

@@ -0,0 +1,215 @@
+% Create an overview plot featuring the results of the multivariate PLS
+% comparing spectral changes during the stimulus period under load
+
+clear all; cla; clc;
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data         = fullfile(rootpath, 'data', 'stats');
+pn.figures      = fullfile(rootpath, 'figures');
+pn.tools        = fullfile(rootpath, 'tools');
+    addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
+    addpath(fullfile(pn.tools, 'fieldtrip')); ft_defaults;
+    addpath(genpath(fullfile(pn.tools, 'RainCloudPlots')));
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'winsor'));
+
+% set custom colormap
+cBrew = brewermap(500,'RdBu');
+cBrew = flipud(cBrew);
+colormap(cBrew)
+
+load(fullfile(pn.data, 'b4b_taskpls_erf.mat'),...
+    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')
+load(fullfile(rootpath, 'data','erp', ['sub-001_erp.mat']));
+elec = erp.scene_category{1}.elec;
+
+result.perm_result.sprob
+
+indLV = 1;
+
+lvdat = reshape(result.boot_result.compare_u(:,indLV), num_chans, num_freqs, num_time);
+stat.prob = lvdat;
+stat.mask = lvdat > 3 | lvdat < -3;
+
+% maskNaN = double(stat.mask);
+% maskNaN(maskNaN==0) = NaN;
+
+%% invert solution
+
+% stat.mask = stat.mask;
+% stat.prob = stat.prob.*-1;
+% result.vsc = result.vsc.*-1;
+% result.usc = result.usc.*-1;
+
+h = figure('units','centimeter','position',[0 0 15 10]);
+set(gcf,'renderer','Painters')
+statsPlot = [];
+%statsPlot = cat(1, statsPlot,squeeze(nanmean(stat.prob(1:64,:,:),1)));
+statsPlot = cat(1, statsPlot,squeeze(nanmean(stat.mask(1:64,:,:).*stat.prob(1:64,:,:),1)));
+range_plot = [-.5 0.25];
+imagesc(stat.time,[],statsPlot, range_plot);
+set(gca,'Ydir','Normal');
+% set custom colormap
+ncol = 2000; cBrew = brewermap(ncol,'RdBu'); cBrew = flipud(cBrew);
+% adjust for unequal range
+reldev = min(abs(range_plot))./max(abs(range_plot));
+if abs(range_plot(2))<max(abs(range_plot))
+    nadjust = [1:ncol/2, ncol/2+round(linspace(1, ncol/2, reldev*(ncol/2)))];
+elseif abs(range_plot(1))<max(abs(range_plot))
+    nadjust = [round(linspace(1, ncol/2, reldev*(ncol/2))),ncol/2:ncol];
+else
+    nadjust = 1:ncol;
+end
+cBrew = cBrew(nadjust,:);
+colormap(cBrew); colorbar;
+xlabel('Time [s from stim onset]'); ylabel('Frequency (Hz)');
+title({'ERF changes'; ['p = ', num2str(round(result.perm_result.sprob(indLV),4))]})
+set(findall(gcf,'-property','FontSize'),'FontSize',18)
+set(gca, 'YTickLabels', round(stat.freq(get(gca, 'YTick')),0));
+
+
+%% plot multivariate topographies
+% set custom colormap
+cBrew = brewermap(500,'RdBu');
+cBrew = flipud(cBrew);
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+cfg.zlim = [-2.5 2.5]; 
+cfg.figure = h;
+plotData.powspctrm = squeeze(nanmean(nanmean(stat.mask(:,stat.freq>8 & stat.freq<20,stat.time>.8 & stat.time<1).*...
+    stat.prob(:,stat.freq>8 & stat.freq<20,stat.time>.8 & stat.time<1),3),2)); 
+[~, sortind] = sort(plotData.powspctrm, 'ascend');
+cfg.marker = 'off'; 
+cfg.highlight = 'yes';
+cfg.highlightchannel = plotData.label(sortind(1:6));
+cfg.highlightcolor = [1 1 1];
+cfg.highlightsymbol = '.';
+cfg.highlightsize = 30;
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Mean BSR');
+title('AlphaBeta')
+
+%% Theta
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+cfg.zlim = [-.5 .5]; 
+cfg.figure = h;
+plotData.powspctrm = squeeze(nanmean(nanmean(stat.mask(:,stat.freq>2 & stat.freq<7,stat.time>0 & stat.time<1).*...
+    stat.prob(:,stat.freq>2 & stat.freq<7,stat.time>0 & stat.time<1),3),2)); 
+[~, sortind] = sort(plotData.powspctrm, 'descend');
+cfg.marker = 'off'; 
+cfg.highlight = 'yes';
+cfg.highlightchannel = plotData.label(sortind(1:6));
+cfg.highlightcolor = [1 1 1];
+cfg.highlightsymbol = '.';
+cfg.highlightsize = 30;
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Mean BSR');
+title('Theta')
+
+%% plot using raincloud plot
+
+groupsizes=result.num_subj_lst;
+conditions=lv_evt_list;
+conds = {'forgotten'; 'remembered'};
+condData = []; uData = [];
+for indGroup = 1
+    if indGroup == 1
+        relevantEntries = 1:groupsizes(1)*numel(conds);
+    elseif indGroup == 2
+        relevantEntries = groupsizes(1)*numel(conds)+1:...
+             groupsizes(1)*numel(conds)+groupsizes(2)*numel(conds);
+    end
+    for indCond = 1:numel(conds)
+        targetEntries = relevantEntries(conditions(relevantEntries)==indCond);        
+        condData{indGroup}(indCond,:) = result.vsc(targetEntries,indLV);
+        uData{indGroup}(indCond,:) = result.usc(targetEntries,indLV);
+    end
+end
+
+cBrew(1,:) = 2.*[.3 .1 .1];
+cBrew(2,:) = [.6 .6 .6];
+
+idx_outlier = cell(1); idx_standard = cell(1);
+for indGroup = 1
+    dataToPlot = uData{indGroup}';
+    % define outlier as lin. modulation that is more than three scaled median absolute deviations (MAD) away from the median
+    X = [1 1; 1 2]; b=X\dataToPlot'; IndividualSlopes = b(2,:);
+    outliers = isoutlier(IndividualSlopes, 'median');
+    idx_outlier{indGroup} = find(outliers);
+    idx_standard{indGroup} = find(outliers==0);
+end
+
+h = figure('units','centimeter','position',[0 0 25 10]);
+for indGroup = 1
+    dataToPlot = uData{indGroup}';
+    % read into cell array of the appropriate dimensions
+    data = []; data_ws = [];
+    for i = 1:2
+        for j = 1:1
+            data{i, j} = dataToPlot(:,i);
+            % individually demean for within-subject visualization
+            data_ws{i, j} = dataToPlot(:,i)-...
+                nanmean(dataToPlot(:,:),2)+...
+                repmat(nanmean(nanmean(dataToPlot(:,:),2),1),size(dataToPlot(:,:),1),1);
+            data_nooutlier{i, j} = data{i, j};
+            data_nooutlier{i, j}(idx_outlier{indGroup}) = [];
+            data_ws_nooutlier{i, j} = data_ws{i, j};
+            data_ws_nooutlier{i, j}(idx_outlier{indGroup}) = [];
+            % sort outliers to back in original data for improved plot overlap
+            data_ws{i, j} = [data_ws{i, j}(idx_standard{indGroup}); data_ws{i, j}(idx_outlier{indGroup})];
+        end
+    end
+
+    % IMPORTANT: plot individually centered estimates, stats on uncentered estimates!
+
+    subplot(1,2,indGroup);
+    set(gcf,'renderer','Painters')
+        cla;
+        cl = cBrew(indGroup,:);
+        [~, dist] = rm_raincloud_fixedSpacing(data_ws, [.8 .8 .8],1);
+        h_rc = rm_raincloud_fixedSpacing(data_ws_nooutlier, cl,1,[],[],[],dist);
+        view([90 -90]);
+        axis ij
+    box(gca,'off')
+    set(gca, 'YTickLabels', {conds{2}; conds{1}});
+    yticks = get(gca, 'ytick'); ylim([yticks(1)-(yticks(2)-yticks(1))./2, yticks(2)+(yticks(2)-yticks(1))./1.5]);
+
+    minmax = [min(min(cat(2,data_ws{:}))), max(max(cat(2,data_ws{:})))];
+    xlim(minmax+[-0.2*diff(minmax), 0.2*diff(minmax)])
+    ylabel('scene category'); xlabel({'Brainscore'; '[Individually centered]'})
+
+    % test linear effect
+    curData = [data_nooutlier{1, 1}, data_nooutlier{2, 1}];
+    X = [1 1; 1 2]; b=X\curData'; IndividualSlopes = b(2,:);
+    [~, p, ci, stats] = ttest(IndividualSlopes);
+    title(['M:', num2str(round(mean(IndividualSlopes),3)), '; p=', num2str(round(p,3))])
+end

code/b_scenecat_n1.m

@@ -0,0 +1,224 @@
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'erp');
+pn.data_erf = fullfile(rootpath, 'data', 'erf');
+pn.tools = fullfile(rootpath, 'tools');
+    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'shadedErrorBar'));
+    
+%% load event info
+
+load(fullfile(pn.data_eeg, ['sub-001_task-xxxx_eeg_art.mat']), 'events');
+
+parameter = {'scene_category'; 'old'; 'behavior'; 'subsequent_memory'};
+for ind_param = 1:numel(parameter)
+    conds.(parameter{ind_param}) = unique(events.(parameter{ind_param}));
+end
+
+%% load erp
+
+for ind_id = 1:33
+    id = sprintf('sub-%03d', ind_id);
+    load(fullfile(pn.data_erp, [id,'_erp.mat']));
+    for ind_option = 1:numel(conds.scene_category)
+        if ind_id == 1
+             erpgroup.scene_category.(conds.scene_category{ind_option}) = erp.scene_category{ind_option};
+             erpgroup.scene_category.(conds.scene_category{ind_option}) = ...
+                 rmfield(erpgroup.scene_category.(conds.scene_category{ind_option}), {'avg', 'var', 'dof'});
+             erpgroup.scene_category.(conds.scene_category{ind_option}).dimord = 'sub_chan_time';
+        end
+        erpgroup.scene_category.(conds.scene_category{ind_option}).avg(ind_id,:,:) = erp.scene_category{ind_option}.avg;
+    end
+end
+
+time = erpgroup.scene_category.manmade.time;
+elec = erpgroup.scene_category.manmade.elec;
+channels = erpgroup.scene_category.manmade.label;
+
+%idx_chans = find(ismember(channels, {'O1', 'Oz', 'O2'}));
+%idx_chans = find(ismember(channels, {'PO7', 'PO8'}));
+%idx_chans = find(ismember(channels, {'Pz', 'CPz', 'P1'}));
+
+mergeddata = cat(4, erpgroup.scene_category.manmade.avg, ...
+    erpgroup.scene_category.natural.avg);
+
+%% plot topography of N1
+
+% set custom colormap
+cBrew = brewermap(500,'RdBu');
+cBrew = flipud(cBrew);
+colormap(cBrew)
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+plotData.powspctrm = squeeze(nanmean(nanmean(nanmean(mergeddata(:,:,time>0.08 & time <0.12,1:2),3),1),4))'; 
+[~, sortidx] = sort(plotData.powspctrm, 'ascend');
+idx_chans = sortidx(1:6);
+
+cfg.marker = 'off'; 
+cfg.highlight = 'yes';
+cfg.highlightchannel = plotData.label(idx_chans);
+cfg.highlightcolor = [1 0 0];
+cfg.highlightsymbol = '.';
+cfg.highlightsize = 18;
+cfg.zlim = [-5 5]*10^-4; 
+cfg.figure = h;
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Amplitude');
+% figureName = ['xxx'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');
+
+%% visualize N1 over negative maximum
+
+% avg across channels and conditions
+condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
+
+h = figure('units','centimeters','position',[0 0 10 8]);
+cla; hold on;
+% new value = old value ? subject average + grand average
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
+xlabel('Time (s) from stim onset')
+xlim([-1 2]); %ylim([-.03 .18])
+ylabel({'ERP';'(microVolts)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',12)
+
+h = figure('units','centimeters','position',[0 0 15 10]);
+cla; hold on;
+% new value = old value ? subject average + grand average
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
+ax = gca; ax.YDir = 'reverse';
+xlabel('Time (s) from stim onset')
+xlim([0 .3]); %ylim([-.03 .18])
+ylabel({'ERP';'(microVolts)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',12)
+
+% across all channels
+condAvg = squeeze(nanmean(nanmean(mergeddata(:,:,:,1:2),2),4));
+
+h = figure('units','centimeters','position',[0 0 15 10]);
+cla; hold on;
+% new value = old value ? subject average + grand average
+curData = squeeze(nanmean(mergeddata(:,:,:,1),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+curData = squeeze(nanmean(mergeddata(:,:,:,2),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
+ax = gca; ax.YDir = 'reverse';
+xlabel('Time (s) from stim onset')
+xlim([0 .3]); %ylim([-.03 .18])
+ylabel({'ERP';'(microVolts)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',12)
+
+%% plot difference between conditions
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+plotData.powspctrm = squeeze(nanmean(nanmean(nanmean(mergeddata(:,:,time>0.08 & time <0.12,2),3),1),4))'...
+    -squeeze(nanmean(nanmean(nanmean(mergeddata(:,:,time>0.08 & time <0.12,1),3),1),4))'; 
+[~, sortidx] = sort(plotData.powspctrm, 'ascend');
+idx_chans = sortidx(1:6);
+
+cfg.marker = 'off'; 
+cfg.highlight = 'yes';
+cfg.highlightchannel = plotData.label(idx_chans);
+cfg.highlightcolor = [1 0 0];
+cfg.highlightsymbol = '.';
+cfg.highlightsize = 18;
+%cfg.zlim = [-5 5]*10^-4; 
+cfg.figure = h;
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Amplitude');
+
+
+%% visualize for negative pls channels
+
+idx_chans = [28:30];
+condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
+
+h = figure('units','centimeters','position',[0 0 10 8]);
+cla; hold on;
+% new value = old value ? subject average + grand average
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
+xlabel('Time (s) from stim onset')
+xlim([-1 2]); %ylim([-.03 .18])
+ylabel({'ERP';'(microVolts)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',12)
+
+%% visualize for positive pls channels
+
+idx_chans = [21:23, 58:62];
+condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
+
+h = figure('units','centimeters','position',[0 0 10 8]);
+cla; hold on;
+% new value = old value ? subject average + grand average
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
+xlabel('Time (s) from stim onset')
+xlim([-1 2]); %ylim([-.03 .18])
+ylabel({'ERP';'(microVolts)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',12)

code/b_scenecat_n1_bl.m

@@ -0,0 +1,248 @@
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'erp');
+pn.data_erf = fullfile(rootpath, 'data', 'erf');
+pn.tools = fullfile(rootpath, 'tools');
+    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'shadedErrorBar'));
+    
+%% load event info
+
+load(fullfile(pn.data_eeg, ['sub-001_task-xxxx_eeg_art.mat']), 'events');
+
+parameter = {'scene_category'; 'old'; 'behavior'; 'subsequent_memory'};
+for ind_param = 1:numel(parameter)
+    conds.(parameter{ind_param}) = unique(events.(parameter{ind_param}));
+end
+
+%% load erp_bl
+
+for ind_id = 1:33
+    id = sprintf('sub-%03d', ind_id);
+    load(fullfile(pn.data_erp, [id,'_erp_bl.mat']));
+    for ind_option = 1:numel(conds.scene_category)
+        if ind_id == 1
+             erpgroup.scene_category.(conds.scene_category{ind_option}) = erp_bl.scene_category{ind_option};
+             erpgroup.scene_category.(conds.scene_category{ind_option}) = ...
+                 rmfield(erpgroup.scene_category.(conds.scene_category{ind_option}), {'avg', 'var', 'dof'});
+             erpgroup.scene_category.(conds.scene_category{ind_option}).dimord = 'sub_chan_time';
+        end
+        erpgroup.scene_category.(conds.scene_category{ind_option}).avg(ind_id,:,:) = erp_bl.scene_category{ind_option}.avg;
+    end
+end
+
+time = erpgroup.scene_category.manmade.time;
+elec = erpgroup.scene_category.manmade.elec;
+channels = erpgroup.scene_category.manmade.label;
+
+%idx_chans = find(ismember(channels, {'O1', 'Oz', 'O2'}));
+%idx_chans = find(ismember(channels, {'PO7', 'PO8'}));
+%idx_chans = find(ismember(channels, {'Pz', 'CPz', 'P1'}));
+
+mergeddata = cat(4, erpgroup.scene_category.manmade.avg, ...
+    erpgroup.scene_category.natural.avg);
+
+%% plot topography of N1
+
+% set custom colormap
+cBrew = brewermap(500,'RdBu');
+cBrew = flipud(cBrew);
+colormap(cBrew)
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+plotData.powspctrm = squeeze(nanmean(nanmean(nanmean(mergeddata(:,:,time>0.08 & time <0.12,1:2),3),1),4))'; 
+[~, sortidx] = sort(plotData.powspctrm, 'ascend');
+idx_chans = sortidx(1:6);
+
+cfg.marker = 'off'; 
+cfg.highlight = 'yes';
+cfg.highlightchannel = plotData.label(idx_chans);
+cfg.highlightcolor = [1 0 0];
+cfg.highlightsymbol = '.';
+cfg.highlightsize = 18;
+cfg.zlim = [-10 10]*10^-4; 
+cfg.figure = h;
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Amplitude');
+% figureName = ['xxx'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');
+
+%% visualize N1 over negative maximum
+
+% avg across channels and conditions
+condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
+
+% h = figure('units','centimeters','position',[0 0 10 8]);
+% cla; hold on;
+% % new value = old value ? subject average + grand average
+% curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
+% curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+% standError = nanstd(curData,1)./sqrt(size(curData,1));
+% l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+% curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
+% curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+% standError = nanstd(curData,1)./sqrt(size(curData,1));
+% l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
+% xlabel('Time (s) from stim onset')
+% xlim([-.2 .5]); %ylim([-.03 .18])
+% ylabel({'ERP';'(microVolts)'});
+% xlabel({'Time (s)'});    
+% set(findall(gcf,'-property','FontSize'),'FontSize',12)
+
+h = figure('units','centimeters','position',[0 0 10 8]);
+cla; hold on;
+% new value = old value ? subject average + grand average
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
+% ax = gca; ax.YDir = 'reverse';
+xlabel('Time (s) from stim onset')
+xlim([-.05 .3]); %ylim([-.03 .18])
+ylabel({'ERP';'(microVolts)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',12)
+% 
+% % across all channels
+% condAvg = squeeze(nanmean(nanmean(mergeddata(:,:,:,1:2),2),4));
+% 
+% h = figure('units','centimeters','position',[0 0 15 10]);
+% cla; hold on;
+% % new value = old value ? subject average + grand average
+% curData = squeeze(nanmean(mergeddata(:,:,:,1),2));
+% curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+% standError = nanstd(curData,1)./sqrt(size(curData,1));
+% l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+% curData = squeeze(nanmean(mergeddata(:,:,:,2),2));
+% curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+% standError = nanstd(curData,1)./sqrt(size(curData,1));
+% l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
+% ax = gca; ax.YDir = 'reverse';
+% xlabel('Time (s) from stim onset')
+% xlim([0 .3]); %ylim([-.03 .18])
+% ylabel({'ERP';'(microVolts)'});
+% xlabel({'Time (s)'});    
+% set(findall(gcf,'-property','FontSize'),'FontSize',12)
+
+%% plot difference between conditions
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+plotData.powspctrm = squeeze(nanmean(nanmean(nanmean(mergeddata(:,:,time>0.08 & time <0.12,2),3),1),4))'...
+    -squeeze(nanmean(nanmean(nanmean(mergeddata(:,:,time>0.08 & time <0.12,1),3),1),4))'; 
+[~, sortidx] = sort(plotData.powspctrm, 'ascend');
+idx_chans = sortidx(1:6);
+
+cfg.marker = 'off'; 
+cfg.highlight = 'yes';
+cfg.highlightchannel = plotData.label(idx_chans);
+cfg.highlightcolor = [1 0 0];
+cfg.highlightsymbol = '.';
+cfg.highlightsize = 18;
+%cfg.zlim = [-5 5]*10^-4; 
+cfg.figure = h;
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Amplitude');
+
+
+%% visualize for negative pls channels
+
+idx_chans = [28:30];
+condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
+
+h = figure('units','centimeters','position',[0 0 10 8]);
+cla; hold on;
+% new value = old value ? subject average + grand average
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
+xlabel('Time (s) from stim onset')
+xlim([-.2 1]); %ylim([-.03 .18])
+xlim([-.05 .3]);
+ylabel({'ERP';'(microVolts)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',12)
+
+%% visualize for positive pls channels
+
+idx_chans = [21:23, 58:62];
+condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
+
+h = figure('units','centimeters','position',[0 0 10 8]);
+cla; hold on;
+% new value = old value ? subject average + grand average
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
+xlabel('Time (s) from stim onset')
+xlim([-.5 1]); %ylim([-.03 .18])
+xlim([-.05 .3]);
+ylabel({'ERP';'(microVolts)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',12)
+
+%% visualize for frontal channels
+
+idx_chans = [37];%[4,38,39];
+condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
+
+h = figure('units','centimeters','position',[0 0 10 8]);
+cla; hold on;
+% new value = old value ? subject average + grand average
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
+xlabel('Time (s) from stim onset')
+xlim([-1 2]); %ylim([-.03 .18])
+ylabel({'ERP';'(microVolts)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',12)

code/z_archive/b2a_taskPLS_novelty_frontal_erp_plot_trace_tmp.m

@@ -0,0 +1,156 @@
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'erp');
+pn.data_erf = fullfile(rootpath, 'data', 'erf');
+pn.tools = fullfile(rootpath, 'tools');
+    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'shadedErrorBar'));
+    
+%% load erp
+
+for ind_id = 1:33
+    id = sprintf('sub-%03d', ind_id);
+    load(fullfile(pn.data_erp, [id,'_erp.mat']));
+    for ind_option = 1:numel(conds.old)
+        if ind_id == 1
+             erpgroup.old.(conds.old{ind_option}) = erp.old{ind_option};
+             erpgroup.old.(conds.old{ind_option}) = ...
+                 rmfield(erpgroup.old.(conds.old{ind_option}), {'avg', 'var', 'dof'});
+             erpgroup.old.(conds.old{ind_option}).dimord = 'sub_chan_time';
+        end
+        erpgroup.old.(conds.old{ind_option}).avg(ind_id,:,:) = erp.old{ind_option}.avg;
+    end
+end
+
+time = erpgroup.old.old.time;
+elec = erpgroup.old.old.elec;
+channels = erpgroup.old.old.label;
+
+%idx_chans = find(ismember(channels, {'O1', 'Oz', 'O2'}));
+%idx_chans = find(ismember(channels, {'PO7', 'PO8'}));
+%idx_chans = find(ismember(channels, {'Pz', 'CPz', 'P1'}));
+
+mergeddata = cat(4, erpgroup.old.old.avg, ...
+    erpgroup.old.new.avg);
+
+%% plot ERP topography
+
+% set custom colormap
+cBrew = brewermap(500,'RdBu');
+cBrew = flipud(cBrew);
+colormap(cBrew)
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+plotData.powspctrm = squeeze(nanmean(nanmean(nanmean(mergeddata(:,:,time>=0.3& time <=0.5,1:2),3),1),4))'; 
+
+cfg.zlim = [-14 14]*10^-4; 
+cfg.figure = h;
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Amplitude');
+% figureName = ['xxx'];
+% saveas(h, fullfile(pn.figures, figureName), 'epsc');
+% saveas(h, fullfile(pn.figures, figureName), 'png');
+
+%% visualize frontal ERP
+
+idx_chans = [38];
+
+% avg across channels and conditions
+condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
+
+h = figure('units','centimeters','position',[0 0 10 8]);
+cla; hold on;
+% highlight relevant phase in background
+patches.timeVec = [0.3 0.5];
+patches.colorVec = [1 .95 .8];
+for indP = 1:size(patches.timeVec,2)-1
+    YLim = [-8 2]*10^-4;
+    p = patch([patches.timeVec(indP) patches.timeVec(indP+1) patches.timeVec(indP+1) patches.timeVec(indP)], ...
+                [YLim(1) YLim(1)  YLim(2), YLim(2)], patches.colorVec(indP,:));
+    p.EdgeColor = 'none';
+end
+% new value = old value ? subject average + grand average
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
+% ax = gca; ax.YDir = 'reverse';
+legend([l1.mainLine, l2.mainLine],{'old', 'new'}, ...
+    'location', 'southwest'); legend('boxoff')
+xlabel('Time (s) from stim onset')
+xlim([-.05 .6]); %ylim([-.03 .18])
+ylabel({'ERP';'(microVolts)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',14)
+
+%% plot difference
+
+idx_chans = [38];
+
+h = figure('units','centimeters','position',[0 0 10 8]);
+cla; hold on;
+% highlight relevant phase in background
+patches.timeVec = [0.3 0.5];
+patches.colorVec = [1 .95 .8];
+for indP = 1:size(patches.timeVec,2)-1
+    YLim = [-10 5]*10^-5;
+    p = patch([patches.timeVec(indP) patches.timeVec(indP+1) patches.timeVec(indP+1) patches.timeVec(indP)], ...
+                [YLim(1) YLim(1)  YLim(2), YLim(2)], patches.colorVec(indP,:));
+    p.EdgeColor = 'none';
+end
+% new value = old value ? subject average + grand average
+curData = squeeze(mergeddata(:,idx_chans,:,2)-mergeddata(:,idx_chans,:,1));
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+legend([l1.mainLine],{'new-old'}, ...
+    'location', 'southwest'); legend('boxoff')
+xlabel('Time (s) from stim onset')
+xlim([-.05 .6]); %ylim([-.03 .18])
+ylabel({'ERP';'(microVolts)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',14)
+
+
+%% plot difference between conditions
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label; % {1 x N}
+plotData.dimord = 'chan';
+oldminnew = mergeddata(:,:,:,1)-mergeddata(:,:,:,2);
+plotData.powspctrm = squeeze(nanmean(nanmean(oldminnew(:,:,time>=.3 & time <=.5),3),1))'; 
+
+cfg.figure = h;
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Amplitude');

code/z_eegmp_old_new_ERP.m

@@ -0,0 +1,59 @@
+function z_eegmp_old_new_ERP(rootpath, id)
+
+%% paths
+if ismac
+    currentFile = mfilename('fullpath');
+    [pathstr,~,~] = fileparts(currentFile);
+    cd(fullfile(pathstr,'..'))
+    rootpath = pwd;
+    id = 'sub-011';
+end
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'test'); mkdir(pn.data_erp);
+pn.tools = fullfile(rootpath, '..', 'eegmp_preproc', 'tools');
+    addpath(fullfile(pn.tools, 'fieldtrip')); ft_defaults
+
+%% load data_eeg
+
+% load preprocessed eeg data_eeg
+load(fullfile(pn.data_eeg, [id,'_task-xxxx_eeg_art.mat']), 'data_eeg', 'events');
+
+%% further preprocessing
+
+% apply notch filter
+
+cfg = [];
+cfg.dftfilter = 'yes';
+data_eeg = ft_preprocessing(cfg, data_eeg);
+
+%% CSD transform
+    
+% csd_cfg = [];
+% csd_cfg.method = 'spline';
+% data_eeg = ft_scalpcurrentdensity(csd_cfg, data_eeg);
+
+%% single-trial baseline-correction for ERPs
+
+time = data_eeg.time{1};
+
+data_eeg_bl = data_eeg;
+
+for indTrial = 1:numel(data_eeg_bl.trial)
+    curbl = squeeze(nanmean(data_eeg_bl.trial{indTrial}(:,time>-.2 & time<0),2));
+    data_eeg_bl.trial{indTrial} = data_eeg_bl.trial{indTrial}-repmat(curbl,1,numel(time));
+end
+
+%% split tfr and time series data by condition, average across trials
+
+cfg = [];
+cfg.trials = ismember(events.old, 'old');
+erp_bl.old = ft_timelockanalysis(cfg, data_eeg_bl);
+
+cfg = [];
+cfg.trials = ismember(events.old, 'new');
+erp_bl.new = ft_timelockanalysis(cfg, data_eeg_bl);
+
+%% save data
+
+save(fullfile(pn.data_erp, [id, '_erp_bl.mat']), 'erp_bl');

code/z_eegmp_old_new_ERP_START.sh

@@ -0,0 +1,27 @@
+#!/bin/bash
+
+#!/bin/bash
+
+fun_name="z_eegmp_old_new_ERP"
+job_name="eegmp_test"
+
+rootpath="$(pwd)/.."
+rootpath=$(builtin cd $rootpath; pwd)
+
+mkdir ${rootpath}/log
+
+for subj in $(seq 1 33); do
+	cursub=$(printf "sub-%03d\n" ${subj})
+	echo ${cursub}
+	# if not been performed yet
+	#if [ ! -f ${rootpath}/data/erf/sub-$(printf "%03d" ${subj})*_erf.mat ]; then
+		sbatch \
+  			--job-name ${job_name}_${cursub} \
+  			--cpus-per-task 4 \
+  			--mem 6gb \
+  			--time 00:30:00 \
+  			--output ${rootpath}/log/${job_name}_${cursub}.out \
+  			--workdir . \
+  			--wrap=". /etc/profile ; module load matlab/R2016b; matlab -nodisplay -r \"${fun_name}('${rootpath}','${cursub}')\""
+  	#fi
+done

code/z_eegmp_old_new_ERP_trace.m

@@ -0,0 +1,125 @@
+
+currentFile = mfilename('fullpath');
+[pathstr,~,~] = fileparts(currentFile);
+cd(fullfile(pathstr,'..'))
+rootpath = pwd;
+
+pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
+pn.data_erp = fullfile(rootpath, 'data', 'test');
+pn.data_erf = fullfile(rootpath, 'data', 'erf');
+pn.tools = fullfile(rootpath, 'tools');
+    addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
+    addpath(fullfile(pn.tools, 'BrewerMap'));
+    addpath(fullfile(pn.tools, 'shadedErrorBar'));
+    
+%% load erp
+
+for ind_id = 1:33
+    id = sprintf('sub-%03d', ind_id);
+    load(fullfile(pn.data_erp, [id,'_erp_bl.mat']));
+    if ind_id == 1
+         erpgroup.old.old = erp_bl.old;
+         erpgroup.old.old = ...
+             rmfield(erpgroup.old.old, {'avg', 'var', 'dof'});
+         erpgroup.old.old.dimord = 'sub_chan_time';
+
+         erpgroup.old.new = erp_bl.new;
+         erpgroup.old.new = ...
+             rmfield(erpgroup.old.new, {'avg', 'var', 'dof'});
+         erpgroup.old.new.dimord = 'sub_chan_time';
+    end
+    erpgroup.old.old.avg(ind_id,:,:) = erp_bl.old.avg;
+    erpgroup.old.new.avg(ind_id,:,:) = erp_bl.new.avg;
+end
+
+time = erpgroup.old.old.time;
+elec = erpgroup.old.old.elec;
+channels = erpgroup.old.old.label;
+
+mergeddata = cat(4, erpgroup.old.old.avg, ...
+    erpgroup.old.new.avg);
+
+%% visualize frontal ERP
+
+idx_chans = [38];
+
+% avg across channels and conditions
+condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
+
+h = figure('units','centimeters','position',[0 0 10 8]);
+cla; hold on;
+% highlight relevant phase in background
+patches.timeVec = [0.3 0.5];
+patches.colorVec = [1 .95 .8];
+for indP = 1:size(patches.timeVec,2)-1
+    YLim = [-8 2]*10^-4;
+    p = patch([patches.timeVec(indP) patches.timeVec(indP+1) patches.timeVec(indP+1) patches.timeVec(indP)], ...
+                [YLim(1) YLim(1)  YLim(2), YLim(2)], patches.colorVec(indP,:));
+    p.EdgeColor = 'none';
+end
+% new value = old value ? subject average + grand average
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
+curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
+% ax = gca; ax.YDir = 'reverse';
+legend([l1.mainLine, l2.mainLine],{'old', 'new'}, ...
+    'location', 'southwest'); legend('boxoff')
+xlabel('Time (s) from stim onset')
+xlim([-.05 .6]); %ylim([-.03 .18])
+ylabel({'ERP';'(microVolts)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',14)
+
+%% plot difference
+
+idx_chans = [38];
+
+h = figure('units','centimeters','position',[0 0 10 8]);
+cla; hold on;
+% highlight relevant phase in background
+patches.timeVec = [0.3 0.5];
+patches.colorVec = [1 .95 .8];
+for indP = 1:size(patches.timeVec,2)-1
+    YLim = [-10 5]*10^-5;
+    p = patch([patches.timeVec(indP) patches.timeVec(indP+1) patches.timeVec(indP+1) patches.timeVec(indP)], ...
+                [YLim(1) YLim(1)  YLim(2), YLim(2)], patches.colorVec(indP,:));
+    p.EdgeColor = 'none';
+end
+% new value = old value ? subject average + grand average
+curData = squeeze(mergeddata(:,idx_chans,:,2)-mergeddata(:,idx_chans,:,1));
+standError = nanstd(curData,1)./sqrt(size(curData,1));
+l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
+legend([l1.mainLine],{'new-old'}, ...
+    'location', 'southwest'); legend('boxoff')
+xlabel('Time (s) from stim onset')
+xlim([-.05 .6]); %ylim([-.03 .18])
+ylabel({'ERP';'(microVolts)'});
+xlabel({'Time (s)'});    
+set(findall(gcf,'-property','FontSize'),'FontSize',14)
+
+%% plot difference between conditions
+
+h = figure('units','centimeters','position',[0 0 10 10]);
+set(gcf,'renderer','Painters')
+
+cfg = [];
+cfg.layout = 'biosemi64.lay';
+cfg.parameter = 'powspctrm';
+cfg.comment = 'no';
+cfg.colormap = cBrew;
+cfg.colorbar = 'EastOutside';
+
+plotData = [];
+plotData.label = elec.label(1:64); % {1 x N}
+plotData.dimord = 'chan';
+oldminnew = mergeddata(:,:,:,1)-mergeddata(:,:,:,2);
+plotData.powspctrm = squeeze(nanmean(nanmean(oldminnew(:,1:64,time>=.3 & time <=.5),3),1))'; 
+
+cfg.figure = h;
+ft_topoplotER(cfg,plotData);
+cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Amplitude');

data/erf/sub-001_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/qv/XM/MD5E-s13011631--618198f320d3695cba8a67d03fab7b5a.mat/MD5E-s13011631--618198f320d3695cba8a67d03fab7b5a.mat

data/erf/sub-002_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/Ff/Mw/MD5E-s11869185--d9c2b10d3c0a5dd9160202cfeaf8ef94.mat/MD5E-s11869185--d9c2b10d3c0a5dd9160202cfeaf8ef94.mat

data/erf/sub-003_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/mZ/PZ/MD5E-s11899882--2e2999433771e61b84f7b38322cff3a8.mat/MD5E-s11899882--2e2999433771e61b84f7b38322cff3a8.mat

data/erf/sub-004_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/GV/fZ/MD5E-s12931033--39021b5851025a6805bd4f1f39f60167.mat/MD5E-s12931033--39021b5851025a6805bd4f1f39f60167.mat

data/erf/sub-005_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/Gv/Mq/MD5E-s12971485--aa8b6e411a586f7ad098955823073cb4.mat/MD5E-s12971485--aa8b6e411a586f7ad098955823073cb4.mat

data/erf/sub-006_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/7F/1v/MD5E-s11851437--ee5c73c68a4fd1e7266b9ee6f7b482aa.mat/MD5E-s11851437--ee5c73c68a4fd1e7266b9ee6f7b482aa.mat

data/erf/sub-007_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/27/KV/MD5E-s12976161--dbfac87116f388a1067b137b9a97079a.mat/MD5E-s12976161--dbfac87116f388a1067b137b9a97079a.mat

data/erf/sub-008_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/7k/m4/MD5E-s11847365--01bdb409bf9a904b16d37aaa25bad5ab.mat/MD5E-s11847365--01bdb409bf9a904b16d37aaa25bad5ab.mat

data/erf/sub-009_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/5X/05/MD5E-s12962734--9c372145748bec57182951ca23d8100c.mat/MD5E-s12962734--9c372145748bec57182951ca23d8100c.mat

data/erf/sub-010_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/vM/4g/MD5E-s12999530--16e890cc703e1d1473b9362019df83cd.mat/MD5E-s12999530--16e890cc703e1d1473b9362019df83cd.mat

data/erf/sub-011_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/kg/Zf/MD5E-s11908825--8cffcf35aa0c30676ad8981dc111b527.mat/MD5E-s11908825--8cffcf35aa0c30676ad8981dc111b527.mat

data/erf/sub-012_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/gW/j3/MD5E-s11887442--bdcc491df2b69c965b680e9694253b38.mat/MD5E-s11887442--bdcc491df2b69c965b680e9694253b38.mat

data/erf/sub-013_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/qX/px/MD5E-s12961623--29065ac232af8402a2422c8b74a9607e.mat/MD5E-s12961623--29065ac232af8402a2422c8b74a9607e.mat

data/erf/sub-014_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/0j/Kg/MD5E-s11913704--f2e99f18643c13cbb46fa99c2dfa4e8b.mat/MD5E-s11913704--f2e99f18643c13cbb46fa99c2dfa4e8b.mat

data/erf/sub-015_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/zv/kM/MD5E-s11885215--be7983d8363359fb1b204cbd16610a91.mat/MD5E-s11885215--be7983d8363359fb1b204cbd16610a91.mat

data/erf/sub-016_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/54/16/MD5E-s12963932--9f6bdf51c2ef40e036a323ce5f855d99.mat/MD5E-s12963932--9f6bdf51c2ef40e036a323ce5f855d99.mat

data/erf/sub-017_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/w1/Gg/MD5E-s11914288--f674d39f6a7ef3cc0352649a857f35b5.mat/MD5E-s11914288--f674d39f6a7ef3cc0352649a857f35b5.mat

data/erf/sub-018_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/4j/p9/MD5E-s12962243--198dc92256655e255cd0aad0a53f8963.mat/MD5E-s12962243--198dc92256655e255cd0aad0a53f8963.mat

data/erf/sub-019_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/1W/f2/MD5E-s12967756--7a8dcb024456c8cf25dab80f54b88074.mat/MD5E-s12967756--7a8dcb024456c8cf25dab80f54b88074.mat

data/erf/sub-020_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/PJ/8K/MD5E-s11915343--d300ec4cded166492669c07260629017.mat/MD5E-s11915343--d300ec4cded166492669c07260629017.mat

data/erf/sub-021_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/q3/4m/MD5E-s12960339--058d843ae534ad32a647e5768c9afc4a.mat/MD5E-s12960339--058d843ae534ad32a647e5768c9afc4a.mat

data/erf/sub-022_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/m0/2Z/MD5E-s11834806--893da53041cbceddc15b23bf4db6f013.mat/MD5E-s11834806--893da53041cbceddc15b23bf4db6f013.mat

data/erf/sub-023_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/wp/j6/MD5E-s12969754--7bf33fd675b5af9bf6f36f38f813b9d3.mat/MD5E-s12969754--7bf33fd675b5af9bf6f36f38f813b9d3.mat

data/erf/sub-024_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/f0/7K/MD5E-s11891030--60f0817f17cae7f9e75ea7fa447d1f1b.mat/MD5E-s11891030--60f0817f17cae7f9e75ea7fa447d1f1b.mat

data/erf/sub-025_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/5k/jw/MD5E-s12986876--a4e76ba0bf4853f2ad30acd189d70e6e.mat/MD5E-s12986876--a4e76ba0bf4853f2ad30acd189d70e6e.mat

data/erf/sub-026_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/0j/FJ/MD5E-s12976047--24299a8c737da252d56b97a6d8480f13.mat/MD5E-s12976047--24299a8c737da252d56b97a6d8480f13.mat

data/erf/sub-027_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/K3/pF/MD5E-s12997497--026e54ce6f1e85b20414613bc9fa88f2.mat/MD5E-s12997497--026e54ce6f1e85b20414613bc9fa88f2.mat

data/erf/sub-028_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/Qp/Qx/MD5E-s12986598--73e61f013b38a3567aa8cc24a6bd2733.mat/MD5E-s12986598--73e61f013b38a3567aa8cc24a6bd2733.mat

data/erf/sub-029_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/jV/Fj/MD5E-s12982199--6c29544a9b39c25ca855df0524ec07af.mat/MD5E-s12982199--6c29544a9b39c25ca855df0524ec07af.mat

data/erf/sub-030_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/89/MX/MD5E-s12962834--84c7ac2d0834f9f371ef980edcb543ae.mat/MD5E-s12962834--84c7ac2d0834f9f371ef980edcb543ae.mat

data/erf/sub-031_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/4w/Mz/MD5E-s12981115--1691c4bddfeba645a9357a13df6d5dd9.mat/MD5E-s12981115--1691c4bddfeba645a9357a13df6d5dd9.mat

data/erf/sub-032_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/j3/pg/MD5E-s11835025--255092095f6c8a00aeee06bc08b7e5ba.mat/MD5E-s11835025--255092095f6c8a00aeee06bc08b7e5ba.mat

data/erf/sub-033_erf.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/G6/kF/MD5E-s11864811--188a6bed6fc0d89a63a459ded09bb149.mat/MD5E-s11864811--188a6bed6fc0d89a63a459ded09bb149.mat

data/erp/sub-001_erp.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/5Z/Q2/MD5E-s18856740--c86ecd68e17ee00a22976fcb19528038.mat/MD5E-s18856740--c86ecd68e17ee00a22976fcb19528038.mat

data/erp/sub-001_erp_bl.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/X4/m8/MD5E-s18850165--a512d0c843e81001e470c4cb21f82e2e.mat/MD5E-s18850165--a512d0c843e81001e470c4cb21f82e2e.mat

data/erp/sub-002_erp.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/w3/59/MD5E-s17943384--02d3d81e912025a9f36836d895d11049.mat/MD5E-s17943384--02d3d81e912025a9f36836d895d11049.mat

data/erp/sub-002_erp_bl.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/QP/f7/MD5E-s17937290--687f3f7e543ce83fbc968d2201903b35.mat/MD5E-s17937290--687f3f7e543ce83fbc968d2201903b35.mat

data/erp/sub-003_erp.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/J6/4j/MD5E-s17977217--0e849452fd478ffe97e8e4e9baf66be4.mat/MD5E-s17977217--0e849452fd478ffe97e8e4e9baf66be4.mat

data/erp/sub-003_erp_bl.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/xv/0Q/MD5E-s17975600--c29e1d75533cede6665ca1b92e4febe8.mat/MD5E-s17975600--c29e1d75533cede6665ca1b92e4febe8.mat

data/erp/sub-004_erp.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/4w/p2/MD5E-s18867198--1b2f1439953d4e0101916c5a278abadd.mat/MD5E-s18867198--1b2f1439953d4e0101916c5a278abadd.mat

data/erp/sub-004_erp_bl.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/7q/8j/MD5E-s18868863--b1e494af08c030b0f9f6bb2ae0ded6bb.mat/MD5E-s18868863--b1e494af08c030b0f9f6bb2ae0ded6bb.mat

data/erp/sub-005_erp.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/8Q/jF/MD5E-s18858655--2bde988ba7774f9ecb3d198fb5819a34.mat/MD5E-s18858655--2bde988ba7774f9ecb3d198fb5819a34.mat

data/erp/sub-005_erp_bl.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/GZ/Jz/MD5E-s18842612--db23f809166f3d959c96ce1b06e29ce6.mat/MD5E-s18842612--db23f809166f3d959c96ce1b06e29ce6.mat

data/erp/sub-006_erp.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/kk/2p/MD5E-s17940517--18f70987afbfe7db807c3e6f89eb2b32.mat/MD5E-s17940517--18f70987afbfe7db807c3e6f89eb2b32.mat

data/erp/sub-006_erp_bl.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/VG/41/MD5E-s17925532--b4aee9072123bf98a3065d1f5ee4ceff.mat/MD5E-s17925532--b4aee9072123bf98a3065d1f5ee4ceff.mat

data/erp/sub-007_erp.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/0X/M2/MD5E-s19566600--b3c9e27f4fe6be145994e900c0f17075.mat/MD5E-s19566600--b3c9e27f4fe6be145994e900c0f17075.mat

data/erp/sub-007_erp_bl.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/50/f1/MD5E-s19571255--19d50cad20d2714eefa98bd3a09090d1.mat/MD5E-s19571255--19d50cad20d2714eefa98bd3a09090d1.mat

data/erp/sub-008_erp.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/j0/G5/MD5E-s17939745--7eddf80c85c043fbda01dacc3c141583.mat/MD5E-s17939745--7eddf80c85c043fbda01dacc3c141583.mat

data/erp/sub-008_erp_bl.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/Mm/Fm/MD5E-s17932101--7c3deebbead27909ac4a3ce026e37e68.mat/MD5E-s17932101--7c3deebbead27909ac4a3ce026e37e68.mat

data/erp/sub-009_erp.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/6g/m1/MD5E-s19586555--66ce978ff4be02b290b906ea2745379d.mat/MD5E-s19586555--66ce978ff4be02b290b906ea2745379d.mat

data/erp/sub-009_erp_bl.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/pV/jV/MD5E-s19571289--059591711ef9242201875fcb76a72fa8.mat/MD5E-s19571289--059591711ef9242201875fcb76a72fa8.mat

data/erp/sub-010_erp.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/6F/P4/MD5E-s18857427--788631d31df8cddccf9c96e01a205934.mat/MD5E-s18857427--788631d31df8cddccf9c96e01a205934.mat

data/erp/sub-010_erp_bl.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/JG/gp/MD5E-s18855910--4a2ad0f3bb399b38c3cd9fca8a54da62.mat/MD5E-s18855910--4a2ad0f3bb399b38c3cd9fca8a54da62.mat

data/erp/sub-011_erp.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/W3/Zk/MD5E-s17944934--8d9ee244f63d42b27f9f34e057a0daae.mat/MD5E-s17944934--8d9ee244f63d42b27f9f34e057a0daae.mat

data/erp/sub-011_erp_bl.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/8g/68/MD5E-s17942065--4346be1a8596203dd163212b9f06aee5.mat/MD5E-s17942065--4346be1a8596203dd163212b9f06aee5.mat

data/erp/sub-012_erp.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/Fj/5x/MD5E-s17971931--82eaf726be0400442e2bf901f4518837.mat/MD5E-s17971931--82eaf726be0400442e2bf901f4518837.mat

data/erp/sub-012_erp_bl.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/Zz/8g/MD5E-s17963130--8b8e9fc9ffc62aede4ec6c6c4b5231e7.mat/MD5E-s17963130--8b8e9fc9ffc62aede4ec6c6c4b5231e7.mat

data/erp/sub-013_erp.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/pm/M1/MD5E-s19580130--1bda614c700b7c9d7e7bdeaffae2e495.mat/MD5E-s19580130--1bda614c700b7c9d7e7bdeaffae2e495.mat

data/erp/sub-013_erp_bl.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/Qq/WW/MD5E-s19583510--57643bfda3091f2b245a2b4cc1cb0516.mat/MD5E-s19583510--57643bfda3091f2b245a2b4cc1cb0516.mat

data/erp/sub-014_erp.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/FG/VQ/MD5E-s17970794--5e106e74f403ca5f163a5d3a28871433.mat/MD5E-s17970794--5e106e74f403ca5f163a5d3a28871433.mat

data/erp/sub-014_erp_bl.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/wM/km/MD5E-s17966722--d18a02249d1324ad9cc1cb5f4d3fb68f.mat/MD5E-s17966722--d18a02249d1324ad9cc1cb5f4d3fb68f.mat

data/erp/sub-015_erp.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/ZJ/gV/MD5E-s17965313--03b16fb2b682a048af6bb40ff10b6a4c.mat/MD5E-s17965313--03b16fb2b682a048af6bb40ff10b6a4c.mat

data/erp/sub-015_erp_bl.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/53/j0/MD5E-s17958084--7f59057057f79bb170a5e07e17d2e4ca.mat/MD5E-s17958084--7f59057057f79bb170a5e07e17d2e4ca.mat

data/erp/sub-016_erp.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/kZ/vg/MD5E-s19575544--5b7580b7304377e2781ee43974a41365.mat/MD5E-s19575544--5b7580b7304377e2781ee43974a41365.mat

data/erp/sub-016_erp_bl.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/W4/5k/MD5E-s19580731--e6e61604a87da32b39b61e3a26f95e06.mat/MD5E-s19580731--e6e61604a87da32b39b61e3a26f95e06.mat

data/erp/sub-017_erp.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/mw/Pq/MD5E-s17937951--e90d30d1b7929a7bdefcdd074750005a.mat/MD5E-s17937951--e90d30d1b7929a7bdefcdd074750005a.mat

data/erp/sub-017_erp_bl.mat

@@ -0,0 +1 @@
+../../.git/annex/objects/FG/vQ/MD5E-s17943499--15f272343aa3e976d081362768058c0d.mat/MD5E-s17943499--15f272343aa3e976d081362768058c0d.mat