SubcorticalDevianceDetectionHumans/Deviance Detection_Humans/DD_2wayANOVAnPhPaired.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)

%% concatenate relevant p-values in new array

nStPa = numel(stPaPh); % number of statistical parameters included in post hoc analysis

% preallocate
pVUc = zeros(nStims,nFilt,nStPa,nRepRate); % preallocate
pV = zeros(nStims,nFilt,(nStPa*nRepRate)); % preallocate
pVSort = zeros(nStims,nFilt,(nStPa*nRepRate)); % preallocate
pVOrd = zeros(nStims,nFilt,(nStPa*nRepRate)); % 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,:) = p(:,s,ft,w,stPaPh(pa)); % extract relevant p-Values and concatenate them in new matrix
            pV(s,ft,(nRepRate*(pa-1))+1:(nRepRate*(pa-1))+nRepRate) = p(:,s,ft,w,stPaPh(pa)); % 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*nRepRate));
pVCo = zeros(nStims,nFilt,2,nStPa,nRepRate);

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,(nRepRate*(pa-1))+1:(nRepRate*(pa-1))+nRepRate);
            end
        end
    end
end

%% save data

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