SubcorticalDevianceDetectionHumans/Deviance Detection_Humans/DD_ConSAnalysisPh.m

%% load ANOVA data

load('DD_ConSData')

%% 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,nConS-1); % preallocate
pV = zeros(nStims,nFilt,(nStPa*nRepRate*(nConS-1))); % preallocate
pVSort = zeros(nStims,nFilt,(nStPa*nRepRate*(nConS-1))); % preallocate
pVOrd = zeros(nStims,nFilt,(nStPa*nRepRate*(nConS-1))); % 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,:,:) = conSP(:,s,:,ft,w,stPaPh(pa));
        end
        pV(s,ft,:) = reshape(pVUc(s,ft,:,:,:),1,(nStPa*nRepRate*(nConS-1)));
        [pVSort(s,ft,:),pVOrd(s,ft,:)] = sort(pV(s,ft,:));
    end
end

%% perform post hoc correction

% preallocate
pVCoTemp = zeros(nStims,nFilt,2,(nStPa*nRepRate*(nConS-1)));
pVCo = zeros(nStims,nFilt,2,nStPa,nRepRate,nConS-1);

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
            pVCo(s,ft,ph,:,:,:) = reshape(pVCoTemp(s,ft,ph,:),nStPa,nRepRate,(nConS-1)); % restore original organisation of uncorrected p-values
        end
    end
end

%% save data

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