% function [] = RunRegression(i,m)
%i = subject, m = model number

%Single trial regression scripts provided by AGF for single-trial regression workshop at MCC 2019. 
%COMPASS ICA analyses by JRW. Thanks to Theo Gründler (fconv_tg) & Gerhard Jocham.
%If you use the scripts, please quote at least Fischer & Ullsperger, Neuron, 2013 & Fischer et al., Scientific Reports, 2016. 
%AGF, 2019

%% 0. Get started
clear; clc; close all

CF = matlab.desktop.editor.getActive; cd(fileparts(CF.Filename)); %just change to folder of this file

path('STATB/',path);path('STATB/helpers/',path);  %add path of STA-TB to matlab

%define data and output paths
pathEEG = 'data/';                  % path to EEG
pathBeh = 'Part_Behav/';            % path to behavioural data
pathAna = 'AnalysisInfo/';          % path to information about analysis steps
pathIn  = 'output/for_across/';     % calculated 'Simple_Error_Model' for 16 participants
pathPic = 'plots/across/';          % output of plots
pathPart= 'PartInfo/';              % Path to information about participant


%% 1. Average and t-test
s1                = [];

ModelName         = 'Simple_Error_Model';
s1.plotelect      = 'Cz'; % plot electrode

s1.plotsteps      = [-200 -100 0 50 100 150 200 300 400 500]; % define time points of topo plots
s1.ERP_Im_plot    = 1;    % plot an ERP image
s1.PlotPval       = 1;    % plot a time course for the p-values
s1.UseValues      = 'b';  % use b or t values?
s1.GroupStats     = 0;    % plot group statistics or average within participant values
s1.AddString      = [' ' s1.UseValues]; %we add the value used to the output file
s1.exclude        = [];

STA_Plot_Regression( [pathIn ModelName '/'], [pathPic ModelName], s1 )

%% 1.1 Model on ICs
ModelName         = 'Simple_Error_Model_IC_Based';
s1.plotelect      = 'ErrC1'; % plot electrode

STA_Plot_Regression( [pathIn ModelName '/'], [pathPic ModelName], s1 )

%% 1.2 Following trial model on ICs
ModelName         = 'Following_Trial_Model_IC_Based';                                                                                                                                                                                                                                                          
s1.plotelect      = 'ErrC1'; % plot electrode
s1.PlotR2         = 0; % plot electrode
s1.UseValues      = 't'; 
s1.MaskPval       = 0.05; 
s1.GroupStats     = 1;

STA_Plot_Regression( [pathIn ModelName '/'], [pathPic ModelName], s1 )


%% 2. Group stats
ModelName         = 'Simple_Error_Model';
s1.plotelect      = 'Cz'; % plot electrode
s1.plotsteps      = [-200 -100 0 50 100 150 200 300 400 500]; % define time points of topo plots
s1.ERP_Im_plot    = 1;    % plot an ERP image
s1.PlotPval       = 1;    % plot a time course for the p-values
s1.UseValues      = 'b';  % use b or t values?
s1.GroupStats     = 1;    % plot group statistics or average within participant values
s1.AddString      = [' ' s1.UseValues]; %we add the value used to the output file
s1.exclude        = [];
s1.GroupStats     = 1;    % plot group statistics or average within participant values
s1.AddString      = [' ' s1.UseValues ' group t']; %we add the value used to the output file

STA_Plot_Regression( [pathIn ModelName '/'], [pathPic ModelName], s1 )


%% 3. FWER correction: Bonferroni
%Bonferroni overpowers everything
s1.GroupStats     = 1;          % plot group statistics or average within participant values
s1.MaskPval       = 'bonf3';    % Bonferroni correction over every data point, electrode, and regressors
s1.AddString      = [' ' s1.UseValues ' group t bonferroni over everything'];
s1.PlotReg        = [];

STA_Plot_Regression( [pathIn ModelName '/'], [pathPic ModelName], s1 )

%% 3.1 Bonferroni accounting for only one regressor
s1.PlotReg        = [1];        % Plots only the regressors selected here
s1.AddString      = [' ' s1.UseValues ' group t bonferroni for accuracy'];

STA_Plot_Regression( [pathIn ModelName '/'], [pathPic ModelName], s1 )


%% 4. FDR
s1.PlotReg        = []; 
s1.MaskPval       = 'fdr';
s1.fdr_q          = 0.05;
s1.AddString      = [' ' s1.UseValues ' group t fdr'];

STA_Plot_Regression( [pathIn ModelName '/'], [pathPic ModelName], s1 )

%% 4.1 FDR with lower q
s1.fdr_q          = 0.05 / 3;
s1.AddString      = [' ' s1.UseValues ' group t fdr q three regressors'];

STA_Plot_Regression( [pathIn ModelName '/'], [pathPic ModelName], s1 )

%% 4.2 FDR with even lower q
s1.fdr_q          = 0.05 / (3 * 7);
s1.AddString      = [' ' s1.UseValues ' group t fdr q three regressors and electrodes'];

STA_Plot_Regression( [pathIn ModelName '/'], [pathPic ModelName], s1 )


%% 5. Comparing groups
%load information
F         = dir([pathPart '*.mat']);                    % get files in participant information folder
age       = nan(length(F),1);sex = nan(length(F),1);    % pre-allocate
for c = 1 : length(F) %loop through all participants
    load([pathPart F(c).name]); %load participant information
    age(c) = VP.Age;            %assign age
    sex(c) = VP.Gender;         %assign sex
end

%add groups
s2.plotsteps      = [-200 -100 0 50 100 150 200 300 400 500]; % define time points of topo plots
s2.ERP_Im_plot    = 1;    % plot an ERP image
s2.PlotPval       = 1;    % plot a time course for the p-values
s2.PlotR2         = 1;
s2.MaskPval       = 'fdr';
s2.fdr_q          = 0.05 / 7;
s2.GroupStats     = 0;
s2.UseValues      = 'b';
s2.PlotReg        = [1];
s2.Groups         = {[find(sex==1)]';[find(sex==2)]'};  %use sex to define male and female groups
s2.GroupNames     = {'male' 'female'};                  %provide a name for the groups
s2.AddString      = [' ' s2.UseValues ' group t fdr electrodes']; %just a string that is added to the filename of the figure

STA_Plot_Regression( [pathIn ModelName '/'], [pathPic ModelName], s2 )

%% 5.1 Use the same threshold
s2.MaskPval       = 0.0024; %threshold from female dataset
s2.AddString      = [' ' s2.UseValues ' group t female threshold three regressors and electrodes'];

STA_Plot_Regression( [pathIn ModelName '/'], [pathPic ModelName], s2 )

%% 5.2 Contrast effects across groups
ModelName = 'Simple_Error_Model';

DM_2ndL = table(sex, age)                   %make a table from the variables

s2.UseValues    = 't';                      %use b values
s2.ModelName    = 'workshop_model';         %name for the 2nd level model (will be created within the first level plot folder by default)
s2.RegLables    = [{'male' 'female'}; {'younger' 'older'}]; %lables
s2.UseReg       = 1;                        %reduce to effects of accuracy regressor alone
s2.ERPimage     = 1;                        %return an ERP image (which at the 2nd level reflects individual participant averages)

STA_2nd_Level_Regression([pathIn ModelName '/'], DM_2ndL, s2) %run second level model

%% 5.3 Plot results
infold             = [pathIn ModelName '/2ndlevel/workshop_model/']; %the output path from before
F                  = dir([infold '*.mat']);  %find the output files -->  'workshop_model Accuracy.mat'

sp2.plotsteps      = [-200 -100 0 50 100 150 200 300 400 500];
sp2.plotelect      = 'Cz'; 
sp2.FileName       = F(1).name;     %set this as the input file name    
sp2.MaskPval       = 0.3;           %we want to see the topographies well, so we use a low threshold
sp2.level          = 2;             %define that this is a 2nd level statistic plot
sp2.ERP_Im_plot    = 1;             %plot an ERP image (these are so to say within subject averages!)
sp2.ExtrapolS      = length(age);   %overwrite default extrapolation setting to number of subjects for ERP image 
sp2.ERPimageSm     = 5;             %because ERP image is average subject data, we do not want to smooth this too much
sp2.AddString      = [' 2nd Level Regression Model'];

STA_Plot_Regression(infold, [pathPic ModelName], sp2 )
























