SubcorticalDevianceDetectionHumans/Deviance Detection_Humans/DD_2wayANOVAnPh.m

%% load ANOVA data

load('DD_2WanovanData')

%% define some variables

stPaPh = [1,4]; % define variables that should be included in post hoc analysis (depending on anova results and research question)
rePV = [1,10,15]; % define relevant p-values from multiple comparison

%% concatenate relevant p-values in new array

nStPa = numel(stPaPh); % number of statistical parameters included in post hoc analysis
nRePV = numel(rePV); % number of relevant p-values (= number of repetition rates)

% preallocate
pVUc = zeros(nStims,nFilt,nStPa,nRePV); % preallocate
pV = zeros(nStims,nFilt,(nStPa*nRePV)); % preallocate
pVSort = zeros(nStims,nFilt,(nStPa*nRePV)); % preallocate
pVOrd = zeros(nStims,nFilt,(nStPa*nRePV)); % preallocate

for s = 1:nStims % for each stimulus
    for ft = 1:nFilt % for each filter
        for pa = 1:nStPa % for each statistical parameter
            switch pa
                case 1 % P2P
                    w = 1; % analysis window
                case 2 % latency
                    w = 2;
            end
            pVUc(s,ft,pa,:) = mcoTbl{s,ft,w,stPaPh(pa)}.Pvalue(rePV); % extract relevant p-Values from multcomp table and concatenate them in new matrix
            pV(s,ft,(nRePV*(pa-1))+1:(nRePV*(pa-1))+nRePV) = mcoTbl{s,ft,w,stPaPh(pa)}.Pvalue(rePV); % extract relevant p-Values from multcomp table and concatenate them in new matrix
        end
        [pVSort(s,ft,:),pVOrd(s,ft,:)] = sort(pV(s,ft,:));
    end
end

%% perform post hoc correction

% preallocate
pVCoTemp = zeros(nStims,nFilt,2,(nStPa*nRePV));
pVCo = zeros(nStims,nFilt,2,nStPa,nRePV);

for s = 1:nStims % for each stimulus
    for ft = 1:nFilt % for each filter
        nP = numel(pVSort(s,ft,:));
        for ph = 1:2 % for each post hoc test
            switch ph
                case 1 % Holm-Bonferroni
                    for n = 1:nP
                        cf = nP-(n-1); % correction factor
                        pVCoTemp(s,ft,ph,n) = pVSort(s,ft,n)*cf; % correct p values
                        if n>1
                            if pVCoTemp(s,ft,ph,n)<pVCoTemp(s,ft,ph,n-1)
                                pVCoTemp(s,ft,ph,n) = pVCoTemp(s,ft,ph,n-1); % make sure that p-values never decrease with increasing rank
                            end
                        end
                    end
                case 2 % Benjamini-Hochberg
                    for n = nP:-1:1
                        cf = nP/n; % correction factor
                        pVCoTemp(s,ft,ph,n) = pVSort(s,ft,n)*cf; % correct p values
                        % if n<nP
                        %     if pVCoTemp(s,ft,ph,n)>pVCoTemp(s,ft,ph,n+1)
                        %         pVCoTemp(s,ft,ph,n) = pVCoTemp(s,ft,ph,n+1); % make sure that p-values never increase with decreasing rank
                        %     end
                        % end
                    end
            end
            pVCoTemp(s,ft,ph,pVOrd(s,ft,:)) = pVCoTemp(s,ft,ph,:); % restore original order of p-values
            for pa = 1:nStPa % for each statistical parameter
                pVCo(s,ft,ph,pa,:) = pVCoTemp(s,ft,ph,(nRePV*(pa-1))+1:(nRePV*(pa-1))+nRePV);
            end
        end
    end
end

%% save data

save('DD_2WanovanPhData','pVUc','pVCo','stPaPh','rePV')