eegmp_analysis/code/c1_taskPLS_novelty_frontal_erp.m

% 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, 'c1_taskpls_erp.mat'),...
    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')