% study_stats
% -----------------
% Curso HEP Granada 2018
% Pandelis Perakakis, Luis Ciria
% peraka@ugr.es, lciria@ugr.es
% ------------------------------
%
% start eeglab and load study
% get erpdata
[STUDY, erpdata, times] = std_readerp(STUDY, ALLEEG,...
    'channels', { ALLEEG(1).chanlocs.labels });
%
% to continue the exercise load 'erpdata_HEP'
%
% get one variable for each group*condition 
HEPB1_HF = erpdata{2,1};
HEPB1_LF = erpdata{2,2};
HEPB2_HF = erpdata{4,1};
HEPB2_LF = erpdata{4,2};
%
% if I want to study interactions
% subtract B1-B2
% HEP_HF=HEPB1_HF-HEPB2_HF;
% HEP_LF=HEPB1_LF-HEPB2_LF;
%
% average B1&B2
HEP_HF=HEPB1_HF+HEPB2_HF./2;
HEP_LF=HEPB1_LF+HEPB2_LF./2;
%
% convert data for use in fieldtrip
ft_HEP_HF = erp2ft(HEP_HF, times./1000, chanlocs, 3, fsample);
ft_HEP_LF = erp2ft(HEP_LF, times./1000, chanlocs, 3, fsample);
%
% fieldtrip statistics
%
% add fieldtrip to path without subfolders
% execute : ft_defaults
%
% configure neighbour electrodes using the fieldtrip template
% define configuration parameters
cfg = [];
cfg.method = 'template';
cfg.elecfile = [cd '/biosemi64.lay'];
% electrode layout file originally found at:
% fieldtrip-20161017/template/layout
cfg.template = [cd '/biosemi64_neighb.mat'];
% neighbours template originally found at:
% fieldtrip-20161017/template/neighbours
neighbours = ft_prepare_neighbours(cfg); % prepare neighbours
%
% configure parameters for statistical analysis
% define variables
elecs = 1:64; % select channels
latency = [0.150 0.400]; % select latency of interest
design = [ones(1,20) ones(1,17)*2]; % create the stats design
% configure parameters
cfg = [];
cfg.design = design;
cfg.channel = chanlocs; % chanlocs(elecs);
cfg.latency = latency; % 'all';
cfg.method = 'montecarlo'; % Monte Carlo Method for significance prob.
% cfg.statistic   = 'indepsamplesT';
cfg.statistic = 'ft_statfun_indepsamplesT';
cfg.correctm = 'cluster';
cfg.clusterstatistic = 'maxsum'; 
cfg.minnbchan = 2;
cfg.numrandomization = 4000;
cfg.tail = 0; % -1, 1 or 0 (default = 0); one-sided or two-sided test
cfg.clustertail = 0;
cfg.alpha = 0.025;  
cfg.ivar  = 1;
cfg.neighbours = neighbours;
%
% run analysis
[stat] = ft_timelockstatistics(cfg, ft_HEP_HF{:}, ft_HEP_LF{:});
%
% extract positive clusters
pos_cluster_pvals = [stat.posclusters(:).prob]; 
pos_signif_clust = find(pos_cluster_pvals < stat.cfg.alpha);
pos = ismember(stat.posclusterslabelmat, pos_signif_clust);
[chan,time] = find(pos == 1);
chan_pos = unique(chan);
time_pos = unique(time);
time_range_pos = [stat.time(time_pos(1)) stat.time(time_pos(end))];
% 0.2539 : 0.3438
%
% extract negative clusters
neg_cluster_pvals = [stat.negclusters(:).prob];
neg_signif_clust = find(neg_cluster_pvals < stat.cfg.alpha);
neg = ismember(stat.negclusterslabelmat, neg_signif_clust);
[chan,time] = find(neg == 1);
chan_neg = unique(chan);
time_neg = unique(time);
time_range_neg = [stat.time(time_neg(1)) stat.time(time_neg(end))];
% 0.2578 : 0.3789
%
% auto clusterplot
% average data
cfg = [];
cfg.channel = chanlocs;
latency = [0.150 0.400];
cfg.latency = latency;
GA_HEP_HF = ft_timelockgrandaverage(cfg,ft_HEP_HF{:});  
GA_HEP_LF = ft_timelockgrandaverage(cfg,ft_HEP_LF{:});
%
% subtract HF-LF
cfg = [];
cfg.operation = 'subtract';
cfg.parameter = 'avg';
GA_HF_vs_LF = ft_math(cfg,GA_HEP_HF,GA_HEP_LF);
%
% plot clusters on scalp topography
stat.raw = (GA_HF_vs_LF.avg);
cfg = [];
cfg.subplotsize = [5 5];
%cfg.saveaspng = 'topos';
cfg.alpha  = stat.cfg.alpha;
cfg.parameter = 'stat';
cfg.layout = [cd '/biosemi64.lay'];
ft_clusterplot(cfg, stat);
%
% plot GA heps for significant clusters
% positive cluster
%
% create GAs
cfg = [];
cfg.channel = chanlocs(chan_pos);
GA_HEP_HF = ft_timelockgrandaverage(cfg,ft_HEP_HF{:});  
GA_HEP_LF = ft_timelockgrandaverage(cfg,ft_HEP_LF{:});
%
% average
GA_HF = mean(GA_HEP_HF.avg);
GA_LF = mean(GA_HEP_LF.avg);
% variance
var_HF = mean(GA_HEP_HF.var);
var_LF = mean(GA_HEP_LF.var);
% mean standard error
SEM_HF = sqrt(var_HF)./sqrt(20);
SEM_LF = sqrt(var_LF)./sqrt(17);
%
xx = ft_HEP_LF{1}.time*1000; % time vector
range = 1:205;
% initialise figure
h = figure ('Name', 'Positive cluster');
% plot axes
plot(xx(range),GA_HF(range),'r','LineWidth',2.5); hold on;
plot(xx(range),GA_LF(range),'b--','LineWidth',2.5)
%
% fill HF
xpoints = xx(range);
upper = GA_HF(range)+SEM_HF(range);
lower = GA_HF(range)-SEM_HF(range);
filled = [upper,fliplr(lower)];
xpoints = [xpoints,fliplr(xpoints)];
fillhandle = fill(xpoints,filled,'r');%plot the data
set(fillhandle,'EdgeColor','r','FaceAlpha',0.3,'EdgeAlpha',0.3);
%
% fill LF
xpoints = xx(range);
upper = GA_LF(range)+SEM_LF(range);
lower = GA_LF(range)-SEM_LF(range);
filled = [upper,fliplr(lower)];
xpoints = [xpoints,fliplr(xpoints)];
fillhandle = fill(xpoints,filled,'b');%plot the data
set(fillhandle,'EdgeColor','b','FaceAlpha',0.3,'EdgeAlpha',0.3);

% set figure properties
Fsize = 18; % font size
ylim ([-4 8])
set(get(gca,'YLabel') ,'String','Voltage(uV)',...
    'FontName','Arial',...
    'FontSize',Fsize,'FontWeight','Bold');
set(get(gca,'XLabel'),'String','Time (msec)',...
    'FontName','Arial',...
    'FontSize',Fsize,'FontWeight','Bold');
set(gca,'FontName','Arial','FontSize',Fsize, 'XMinorTick', 'off',...
    'YMinorTick', 'off','box','off','TickLength',[0.02 0.02],'FontWeight','Bold')
%
% shade significant area = 0.2539 : 0.3438
x1 = findnearest (xx,253.9);
x2 = findnearest (xx,343.8);
gr = [0.5 0.5 0.5];
fillhandle = fill ([xx(x1) xx(x2) xx(x2) xx(x1)],[-4 -4 8 8], gr);
set(fillhandle,'EdgeColor','none','FaceAlpha',0.3,'EdgeAlpha',0);