SubcorticalDevianceDetectionHumans/Deviance Detection_Humans/DD_2wayANOVAn.m

%% define some variables

clear

% define time windows to be analysed (in s)
% low
win1L = [0.0027,0.0092]; % ABR (whole)
% win2L = [0.0027,0.0092]; % ABR (peak 5 only)
win2L = [0.005,0.0092]; % ABR (peak 5 only)

% high
% win1H = [0,0.0092]; % ABR (whole)
% win2H = [0.004,0.0092]; % ABR (peak 5 only)
win1H = [0.0027,0.0092]; % ABR (whole)
win2H = [0.005,0.0092]; % ABR (peak 5 only)

% name time windows
winTxt = ["ABR","Wave V"];

%% load and concatenate data

% load and rename variables
load("DD_avData_Oddball_25ms_SbD.mat")
dataAv_25 = dataAv_cell{1};
pts2begin_25 = pts2begin_cell{1}(1);
fs_25 = fs_cell{1}(1);
load("DD_avData_Oddball_30ms_SbD.mat")
dataAv_30 = dataAv_cell{1};
pts2begin_30 = pts2begin_cell{1}(1);
fs_30 = fs_cell{1}(1);
load("DD_avData_Oddball_50ms_SbD.mat")
dataAv_50 = dataAv_cell{1};
pts2begin_50 = pts2begin_cell{1}(1);
fs_50 = fs_cell{1}(1);

% concatenate variables
data = cat(1,dataAv_25,dataAv_30,dataAv_50);
pts2begin = cat(1,pts2begin_25,pts2begin_30,pts2begin_50);
fs = cat(1,fs_25,fs_30,fs_50);

% define some variables
nStims = 2;
nRepRate = size(data,1);
nProb = size(data,2)/nStims;
nFilt = size(data{1},3);
nFiles = size(data{1},2);

%% do some additional calculations

% preallocate
winAll = cell(nRepRate,1);

for r = 1:nRepRate
    winAll{r,1} = round([win1L;win2L]*fs(r))+pts2begin(r)+round(blcEnd(1)*fs(r)+1); % combine time windows, convert them into points and correct for stim delay (pts2begin) and baseline correction window (blcEnd = stimulus peak)
    winAll{r,2} = round([win1H;win2H]*fs(r))+pts2begin(r)+round(blcEnd(2)*fs(r)+1);
end
nWin = size(winAll{1},1); % number of windows per stimulus
nStPa = 5; % number of statistical parameters

%% calculate statistical parameters

% preallocate
stPa = zeros(nFiles,nRepRate,nStims,nProb,nFilt,nWin,nStPa);

for r = 1:nRepRate % for each rep rate
    for s = 1:nStims % for each stimulus (A and B)
        for pr = 1:nProb % for each probability (dev and std)
            for ft = 1:nFilt % for each filter
                for w = 1:nWin % for each window
                    for pa = 1:nStPa % for each stat param
                        switch pa
                            case 1 % peak 2 peak
                                stPa(:,r,s,pr,ft,w,pa) = peak2peak(data{r,(s-1)*nProb+pr}(winAll{r,s}(w,1):winAll{r,s}(w,2),:,ft));
                            case 2 % RMS
                                stPa(:,r,s,pr,ft,w,pa) = rms(data{r,(s-1)*nProb+pr}(winAll{r,s}(w,1):winAll{r,s}(w,2),:,ft));
                            case 3 % AUC
                                stPa(:,r,s,pr,ft,w,pa) = trapz(abs(data{r,(s-1)*nProb+pr}(winAll{r,s}(w,1):winAll{r,s}(w,2),:,ft)));
                            case {4,5} % peak latency & peak width
                                pks = zeros(1,nFiles);
                                loc = zeros(1,nFiles);
                                wid = zeros(1,nFiles);
                                for f = 1:nFiles
                                    [pksTemp,locTemp,widTemp,~] = findpeaks(data{r,(s-1)*nProb+pr}(winAll{r,s}(w,1):winAll{r,s}(w,2),f,ft),"WidthReference",'halfprom');
                                    [maxV,maxI] = max(pksTemp);
                                    pks(f) = pksTemp(maxI);
                                    loc(f) = locTemp(maxI);
                                    wid(f) = widTemp(maxI);
                                end
                                switch pa
                                    case 4
                                        stPa(:,r,s,pr,ft,w,pa) = (loc+(winAll{r,s}(w,1)-pts2begin(r)-round(blcEnd(s)*fs(r)+1)))/fs(r)*1000; % convert points to time
                                    case 5
                                        stPa(:,r,s,pr,ft,w,pa) = wid/fs(r)*1000; % convert points to time
                                end
                        end
                    end
                end
            end
        end
    end
end

%% prepare variables for ANOVA

% preallocate
stPaAno = zeros((nFiles*nRepRate*nProb),nStims,nFilt,nWin,nStPa);
faProb = strings((nFiles*nRepRate*nProb),nStims,nFilt,nWin,nStPa);
faRepRate = strings((nFiles*nRepRate*nProb),nStims,nFilt,nWin,nStPa);

count = 0;

for r = 1:nRepRate % for each rep rate
    switch r
        case 1
            RepRateName = "25ms";
        case 2
            RepRateName = "30ms";
        case 3
            RepRateName = "50ms";
    end
    for pr = 1:nProb % for each probability (dev and std)
        switch pr
            case 1
                probName = "Dev";
            case 2
                probName = "Std";
        end
        count = count+1;
        for s = 1:nStims % for each stimulus
            for ft = 1:nFilt % for each filter
                for w = 1:nWin % for each window
                    for pa = 1:nStPa % for each stat param
                        stPaAno((nFiles*(count-1))+1:(nFiles*(count-1))+nFiles,s,ft,w,pa) = stPa(:,r,s,pr,ft,w,pa); % reorganise stPa to have all conditions that should be tested in anova in one dimension
                        % create factors for anova
                        faProb((nFiles*(count-1))+1:(nFiles*(count-1))+nFiles,s,ft,w,pa) = probName;
                        faRepRate((nFiles*(count-1))+1:(nFiles*(count-1))+nFiles,s,ft,w,pa) = RepRateName;
                    end
                end
            end
        end
    end
end

%% run ANOVA & multiple comparison test (without correction)

% preallocate
anoP = cell(nStims,nFilt,nWin,nStPa);
anoTbl = cell(nStims,nFilt,nWin,nStPa);
anoStats = cell(nStims,nFilt,nWin,nStPa);
anoTerms = cell(nStims,nFilt,nWin,nStPa);
mcoC = cell(nStims,nFilt,nWin,nStPa);
mcoM = cell(nStims,nFilt,nWin,nStPa);
mcoH = cell(nStims,nFilt,nWin,nStPa);
mcoGnames = cell(nStims,nFilt,nWin,nStPa);
mcoTbl = cell(nStims,nFilt,nWin,nStPa);

for s = 1:nStims % for each stimulus
    for ft = 1:nFilt % for each filter
        for w = 1:nWin % for each window
            for pa = 1:nStPa % for each stat param
                y = stPaAno(:,s,ft,w,pa); % data to be tested
                group = {faProb(:,s,ft,w,pa),faRepRate(:,s,ft,w,pa)}; % combine factors in group variable
                [anoP{s,ft,w,pa},anoTbl{s,ft,w,pa},anoStats{s,ft,w,pa},anoTerms{s,ft,w,pa}]  = anovan(y,group,'model','interaction','varnames',{'faProb','faRepRate'},'display','off');
                [mcoC{s,ft,w,pa},mcoM{s,ft,w,pa},mcoH{s,ft,w,pa},mcoGnames{s,ft,w,pa}] = multcompare(anoStats{s,ft,w,pa},'Dimension',[1,2],"CriticalValueType","lsd",'display','off');
                mcoTbl{s,ft,w,pa} = array2table(mcoC{s,ft,w,pa},"VariableNames",["Group A","Group B","Lower Limit","AtoB","Upper Limit","Pvalue"]);
                mcoTbl{s,ft,w,pa}.("Group A") = mcoGnames{s,ft,w,pa}(mcoTbl{s,ft,w,pa}.("Group A"));
                mcoTbl{s,ft,w,pa}.("Group B") = mcoGnames{s,ft,w,pa}(mcoTbl{s,ft,w,pa}.("Group B"));
            end
        end
    end
end

%% calculate paird t-tests and effect size

% preallocate
effS = zeros(nRepRate,nStims,nFilt,nWin,nStPa);
h = zeros(nRepRate,nStims,nFilt,nWin,nStPa);
p = zeros(nRepRate,nStims,nFilt,nWin,nStPa);
ci = cell(nRepRate,nStims,nFilt,nWin,nStPa);
stats = cell(nRepRate,nStims,nFilt,nWin,nStPa);

for r = 1:nRepRate
    for s = 1:nStims % for each stimulus
        for ft = 1:nFilt % for each filter
            for w = 1:nWin % for each window
                for pa = 1:nStPa % for each stat param
                    [h(r,s,ft,w,pa),p(r,s,ft,w,pa),ci{r,s,ft,w,pa},stats{r,s,ft,w,pa}] = ttest(stPa(:,r,s,1,ft,w,pa),stPa(:,r,s,2,ft,w,pa)); % run t-test
                    effS(r,s,ft,w,pa) = meanEffectSize(stPa(:,r,s,1,ft,w,pa),stPa(:,r,s,2,ft,w,pa),'Effect','cohen').Effect; % calculate effect size; dev vs. std
                end
            end
        end
    end
end

%% save data

save('DD_2WanovanData','stPa','stPaAno','faProb','faRepRate','anoP',...
    'anoTbl','anoStats','anoTerms','mcoC','mcoM','mcoH','mcoGnames',...
    'mcoTbl','nRepRate','nStims','nProb','nFilt','nWin','nStPa','effS',...
    'h','p','ci','stats','combName','winAll')