eegmp_analysis/code/c2_taskPLS_novelty_frontal_theta.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 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, 'c2_taskpls_theta.mat'),...
    'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')