SubcorticalDevianceDetectionHumans/Deviance Detection_Humans/DD_2wayANOVA.m

%% define some variables

clear

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

% high
win1H = [0.0027,0.0092]; % ABR (whole)
% win2H = [0.0045,0.008]; % ABR (peak 5 only)
win2H = [0.004,0.008]; % 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:2 % 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 p = 1:nStPa % for each stat param
                        switch p
                            case 1 % peak 2 peak
                                stPa(:,r,s,pr,ft,w,p) = 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,p) = 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,p) = 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));
                                    [maxV,maxI] = max(pksTemp);
                                    pks(f) = pksTemp(maxI);
                                    loc(f) = locTemp(maxI);
                                    wid(f) = widTemp(maxI);
                                end
                                switch p
                                    case 4
                                        stPa(:,r,s,pr,ft,w,p) = (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,p) = 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);

stPaAno = zeros((nFiles*nRepRate),nProb,nStims,nFilt,nWin,nStPa);

count = 0;

for r = 1:nRepRate % for each rep rate
    count = count+1;
    for pr = 1:nProb % for each probability (dev and std)
        for s = 1:nStims % for each stimulus
            for ft = 1:nFilt % for each filter
                for w = 1:nWin % for each window
                    for p = 1:nStPa % for each stat param
                        stPaAno((nFiles*(count-1))+1:(nFiles*(count-1))+nFiles,pr,s,ft,w,p) = stPa(:,r,s,pr,ft,w,p); % reorganise stPa to have all conditions that should be tested in anova in one dimension
                    end
                end
            end
        end
    end
end

%% run ANOVA & post hoc multiple comparison test

% preallocate
anoP = cell(nStims,nFilt,nWin,nStPa);
anoTbl = cell(nStims,nFilt,nWin,nStPa);
anoStats = 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 p = 1:nStPa % for each stat param
                y = stPaAno(:,:,s,ft,w,p); % data to be tested
                reps = nFiles;
                [anoP{s,ft,w,p},anoTbl{s,ft,w,p},anoStats{s,ft,w,p}]  = anova2(y,reps,'off');
                [mcoC{s,ft,w,p},mcoM{s,ft,w,p},mcoH{s,ft,w,p},mcoGnames{s,ft,w,p}] = multcompare(anoStats{s,ft,w,p},"CriticalValueType","bonferroni",'display','off');
                mcoTbl{s,ft,w,p} = array2table(mcoC{s,ft,w,p},"VariableNames",["Group A","Group B","Lower Limit","AtoB","Upper Limit","Pvalue"]);
            end
        end
    end
end

bla = 1;