DevianceDetectionHumanABR/Deviance Detection_Humans/DD_Cbp.m

% Culster based permutation testing
tic
%% define some variables
clear

inclCtrl = 0; % include control data? (0: no, 1: 50Per, 2: MR, 3: DevOnly)
omRespI = 1; % run stats for omission responses? (0: no, 1: yes)

% define number of std-position datasets
switch omRespI 
    case 0
        nSubStd = 2; % 2 for all datasets
    case 1
        nSubStd = 1; % except Om dataset (which has only std after dev)
end

% cluster settings
lenMin = 5;
%% load datasets
% load control data and rename important variables
for repRateI = 1:3 % once for each rep rate dataset
    for subStdI = 1:nSubStd % once for each std position dataset
        switch inclCtrl 
            case {0,3} % use different rep rate datasets only if no ctrl or dev only control is included (because there are no control data for 25 and 30 ms)
                switch repRateI
                    case 1
                        repRate = '25ms';
                    case 2
                        repRate = '30ms';
                    case 3
                        repRate = '50ms';
                end
            case {1,2} % if 50Per or MR control is included, only use 50 ms dataset
                repRate = '50ms';
        end
        switch subStdI
            case 1
                subStd = 'SaD';
            case 2
                subStd = 'SbD';
        end
        switch omRespI
            case 0
                switch inclCtrl
                    case 1
                        load(['DD_avData_50Per_',repRate,'_',subStd,'.mat'])
                        load(['DD_avData_50Per_',repRate,'_',subStd,'_perm.mat'])
                        dataAv_cell_Ctrl = dataAv_cell; % rename
                        dataAvP_cell_Ctrl = dataAvP_cell; % rename
                        fileI_Oddb_cell_Ctrl = fileI_Oddb_cell; % rename
                        fileI_Ctrl_cell_Ctrl = fileI_Ctrl_cell; % rename
                    case 2
                        load(['DD_avData_MR_',repRate,'_',subStd,'.mat'])
                        load(['DD_avData_MR_',repRate,'_',subStd,'_perm.mat'])
                        dataAv_cell_Ctrl = dataAv_cell; % rename
                        dataAvP_cell_Ctrl = dataAvP_cell; % rename
                        fileI_Oddb_cell_Ctrl = fileI_Oddb_cell; % rename
                        fileI_Ctrl_cell_Ctrl = fileI_Ctrl_cell; % rename
                    case 3
                        load(['DD_avData_LowOnly_',repRate,'_',subStd,'.mat'])
                        load(['DD_avData_LowOnly_',repRate,'_',subStd,'_perm.mat'])
                        dataAv_cell_LowOnly = dataAv_cell; % rename
                        dataAvP_cell_LowOnly = dataAvP_cell; % rename
                        fileI_Oddb_cell_Ctrl = fileI_Oddb_cell; % rename - not necessary to do the same with corresponding "high" file indices, since they should be the same between high and low
                        fileI_Ctrl_cell_Ctrl = fileI_Ctrl_cell; % rename - not necessary to do the same with corresponding "high" file indices, since they should be the same between high and low
                        load(['DD_avData_HighOnly_',repRate,'_',subStd,'.mat'])
                        load(['DD_avData_HighOnly_',repRate,'_',subStd,'_perm.mat'])
                        dataAv_cell_HighOnly = dataAv_cell; % rename
                        dataAvP_cell_HighOnly = dataAvP_cell; % rename
                end
                % load oddball data
                load(['DD_avData_Oddball_',repRate,'_',subStd,'.mat'])
                load(['DD_avData_Oddball_',repRate,'_',subStd,'_perm.mat'])
            case 1
                inclCtrl = 0; % no control for OM data
                % load OM data
                load(['DD_avData_HighOm_',repRate,'.mat'])
                load(['DD_avData_HighOm_',repRate,'_perm.mat'])
                dataAv_cell_HighOm = dataAv_cell; % rename
                dataAvP_cell_HighOm = dataAvP_cell; % rename
                load(['DD_avData_LowOm_',repRate,'.mat'])
                load(['DD_avData_LowOm_',repRate,'_perm.mat'])
                dataAv_cell_LowOm = dataAv_cell; % rename
                dataAvP_cell_LowOm = dataAvP_cell; % rename
        end
        %% do some additional calulations that are required
        
        nSmpl = size(dataAv_cell{1}{1},1); % number of points of each recording
        
        nStims = 2; % number of different stimuli per combination
        nFilt = nFilt_cell{1};
        
        switch inclCtrl
            case 0
                switch omRespI
                    case 0
                        fileName = 'DD_CbpData_noCtrl';
                    case 1
                        fileName = 'DD_CbpData_Om';
                        % replace data in dataAv_cell with omission response data
                        for c = 1:nComb
                            dataAv_cell{c}{1} = dataAv_cell_HighOm{c}{3};
                            dataAv_cell{c}{2} = dataAv_cell_HighOm{c}{4};
                            dataAv_cell{c}{3} = dataAv_cell_LowOm{c}{1};
                            dataAv_cell{c}{4} = dataAv_cell_LowOm{c}{2};

                            dataAvP_cell{c}{1} = dataAvP_cell_HighOm{c}{3};
                            dataAvP_cell{c}{2} = dataAvP_cell_HighOm{c}{4};
                            dataAvP_cell{c}{3} = dataAvP_cell_LowOm{c}{1};
                            dataAvP_cell{c}{4} = dataAvP_cell_LowOm{c}{2};
                        end
                end
                nCond = 2; % 2 conditions: dev and std
                data = dataAv_cell; % only oddball data
                dataP = dataAvP_cell; % only oddball data
            case {1,2,3}
                switch inclCtrl
                    case 1
                        fileName = 'DD_statsData_50Per';
                    case 2
                        fileName = 'DD_statsData_MR';
                    case 3
                        fileName = 'DD_statsData_DevO';
                end
                nCond = 3; % 3 conditions: dev, std and control
                data = cell(size(dataAv_cell));
                dataP = cell(size(dataAvP_cell));
                for c = 1:nComb % once for each stimulus combination
                    if inclCtrl==3 % for DevOnly control
                        dataAv_cell_Ctrl{c}{1} = dataAv_cell_LowOnly{c}{1}; % assign DevOnly Low data as control of stimulus 1
                        dataAv_cell_Ctrl{c}{2} = dataAv_cell_HighOnly{c}{3}; % assign DevOnly High data as control of stimulus 2

                        dataAvP_cell_Ctrl{c}{1} = dataAvP_cell_LowOnly{c}{1}; % assign DevOnly Low data as control of stimulus 1
                        dataAvP_cell_Ctrl{c}{2} = dataAvP_cell_HighOnly{c}{3}; % assign DevOnly High data as control of stimulus 2
                    end
                    for s = 1:nStims % once for each stimulus frequency
                        for cn = 1:nCond % once for each stimulus condition (dev, std, ctrl)
                            % observed data
                            switch cn % change dataset depending on condition (Oddball vs. control)
                                case {1,2} % if oddball
                                    data{c}{(s-1)*nCond+cn} = dataAv_cell{c}{(s-1)*2+cn}(:,fileI_Oddb_cell_Ctrl{c},:); % concatenate oddball data of respective stimulus & condition (dev/std)
                                case 3 % if control
                                    data{c}{(s-1)*nCond+cn} = dataAv_cell_Ctrl{c}{s}(:,fileI_Ctrl_cell_Ctrl{c},:); % concatenate control data of respective stimulus
                            end
                            % permutation data
                            for p = 1:nPerm
                                switch cn % change dataset depending on condition (Oddball vs. control)
                                    case {1,2} % if oddball
                                        dataP{c}{(s-1)*nCond+cn,p} = dataAvP_cell{c}{(s-1)*2+cn}(:,fileI_Oddb_cell_Ctrl{c},:,p); % concatenate oddball data of respective stimulus & condition (dev/std)
                                    case 3 % if control
                                        dataP{c}{(s-1)*nCond+cn,p} = dataAvP_cell_Ctrl{c}{s}(:,fileI_Ctrl_cell_Ctrl{c},:,p); % concatenate control data of respective stimulus
                                end
                            end
                        end
                    end
                end
        end
        %% calulate statistical variable (RMS, mean V, etc.) for each sample (points or windows)
        
        varS = cell(1,nComb); % preallocate
        varSP = cell(1,nComb); % preallocate
        for c = 1:nComb % once for each stimulus combination
            nFiles = nFiles_cell{c};
            varS{c} = zeros(nFiles,nStims,nCond,nFilt,nSmpl); % preallocate
            varSP{c} = zeros(nFiles,nStims,nCond,nFilt,nSmpl,nPerm); % preallocate
            for s = 1:nStims % once for each stimulus frequency
                for cn = 1:nCond % once for each stimulus condition (dev, std, (ctrl))
                    for ft = 1:nFilt % once for each filter
                        for sp = 1:nSmpl % once for each sample (windows or points)
                            % observed data
                            varS{c}(:,s,cn,ft,sp) = (data{c}{(s-1)*nCond+cn}(sp,:,ft)); % calculate variable per point
                            % permutation data
                            for p = 1:nPerm
                                varSP{c}(:,s,cn,ft,sp,p) = (dataP{c}{(s-1)*nCond+cn}(sp,:,ft,p)); % calculate variable per point
                            end
                        end
                    end
                end
            end
        end
    	%% run t-test or ANOVA on each time window
        
        cohensD_V = zeros(nComb,nStims,nFilt,nSmpl,3); % preallocate
        
        switch inclCtrl
            case 0
                hV = zeros(nComb,nStims,nFilt,nSmpl); % preallocate
                pV = zeros(nComb,nStims,nFilt,nSmpl); % preallocate
                ciV = zeros(nComb,nStims,nFilt,nSmpl,2); % preallocate
                statsV = cell(nComb,nStims,nFilt,nSmpl); % preallocate
                hVP = zeros(nComb,nStims,nFilt,nSmpl,nPerm); % preallocate
                pVP = zeros(nComb,nStims,nFilt,nSmpl,nPerm); % preallocate
                ciVP = zeros(nComb,nStims,nFilt,nSmpl,2,nPerm); % preallocate
                statsVP = cell(nComb,nStims,nFilt,nSmpl,nPerm); % preallocate
            case {1,2,3}
                anovaV = cell(nComb,nStims,nFilt,nSmpl); % preallocate
                multcompV = cell(nComb,nStims,nFilt,nSmpl); % preallocate
                anovaVP = cell(nComb,nStims,nFilt,nSmpl,nPerm); % preallocate
                multcompVP = cell(nComb,nStims,nFilt,nSmpl,nPerm); % preallocate
                within = table([1,2,3]','VariableNames',{'stim_prob'}); % define within model
        end
        for c = 1:nComb % once for each stimulus combination
            for s = 1:nStims % once for each stimulus frequency
                for ft = 1:nFilt % once for each filter
                    for sp = 1:nSmpl % once for each time window
                        % observed data
                        cohensD_V(c,s,ft,sp,1) = CohensD(varS{c}(:,s,1,ft,sp),varS{c}(:,s,2,ft,sp)); % calculate Cohens D; dev vs. std
                        [hV(c,s,ft,sp),pV(c,s,ft,sp),ciV(c,s,ft,sp,:),statsV{c,s,ft,sp}] = ttest(varS{c}(:,s,1,ft,sp),varS{c}(:,s,2,ft,sp)); % run t-test
                        % permutation data
                        for p = 1:nPerm
                            [hVP(c,s,ft,sp,p),pVP(c,s,ft,sp,p),ciVP(c,s,ft,sp,:,p),statsVP{c,s,ft,sp,p}] = ttest(varSP{c}(:,s,1,ft,sp,p),varSP{c}(:,s,2,ft,sp,p)); % run t-test
                        end
                        switch inclCtrl
                            case {1,2,3} % if control is included
                                % observed data
                                t = table(varS{c}(:,s,1,ft,sp),varS{c}(:,s,2,ft,sp),varS{c}(:,s,3,ft,sp),'VariableNames',{'dev','std','ctrl'}); % required for ANOVA
                                rm = fitrm(t,'dev-ctrl~1','WithinDesign',within); % required for ANOVA
                                anovaV{c,s,ft,sp} = ranova(rm); % run ANOVA
                                multcompV{c,s,ft,sp} = multcompare(rm,'stim_prob','ComparisonType','bonferroni'); % run multiple comparison to find differences between responses
                                cohensD_V(c,s,ft,sp,2) = CohensD(varS{c}(:,s,1,ft,sp),varS{c}(:,s,3,ft,sp)); % calculate Cohens D; dev vs. ctrl
                                cohensD_V(c,s,ft,sp,3) = CohensD(varS{c}(:,s,3,ft,sp),varS{c}(:,s,2,ft,sp)); % ctrl vs. std
                                % permutation data
                                for p = 1:nPerm
                                    anovaVP{c,s,ft,sp,p} = ranova(rm); % run ANOVA
                                    multcompVP{c,s,ft,sp,p} = multcompare(rm,'stim_prob','ComparisonType','bonferroni'); % run multiple comparison to find differences between responses
                                end
                        end
                    end
                end
            end
        end
        %% cluster analysis

        tStat = zeros(1,nSmpl); % preallocate
        vecCl = zeros(nComb,nStims,nFilt,nSmpl);
        vecClI = zeros(nComb,nStims,nFilt,nSmpl);
        nCl = zeros(nComb,nStims,nFilt);
        sCl = cell(nComb,nStims,nFilt);
        tStatP = zeros(1,nSmpl,nPerm); % preallocate
        vecClP = zeros(nComb,nStims,nFilt,nSmpl);
        vecClIP = zeros(nComb,nStims,nFilt,nSmpl);
        nClP = zeros(nComb,nStims,nFilt);
        sClP = cell(nComb,nStims,nFilt);
        sClPmax = zeros(nComb,nStims,nFilt);
        clLaI = cell(nComb,nStims,nFilt);
        pCl = cell(nComb,nStims,nFilt);
        hCl = cell(nComb,nStims,nFilt);
        sigCl = zeros(nComb,nStims,nFilt,nSmpl);
        vecSigCl = zeros(nComb,nStims,nFilt,nSmpl);
        vecSigClI = zeros(nComb,nStims,nFilt,nSmpl);
        nSigCl = zeros(nComb,nStims,nFilt);
        sSigCl = cell(nComb,nStims,nFilt);

        for c = 1:nComb % once for each stimulus combination
            for s = 1:nStims % once for each stimulus frequency
                for ft = 1:nFilt % once for each filter
                    for sp = 1:nSmpl % once for each time window
                        % observed data
                        tStat(sp) = statsV{c,s,ft,sp}.tstat;
                        % permutation data
                        for p = 1:nPerm
                            tStatP(sp,p) = statsVP{c,s,ft,sp,p}.tstat;
                        end
                    end
                    % observed data
                    [vecCl(c,s,ft,:),vecClI(c,s,ft,:),nCl(c,s,ft),sCl{c,s,ft}] = findClu(reshape(hV(c,s,ft,:),1,nSmpl),tStat,lenMin); % identify clusters in observed data
                    % permutation data
                    for p = 1:nPerm
                        [vecClP(c,s,ft,:,p),vecClIP(c,s,ft,:),nClP(c,s,ft,p),sClP{c,s,ft,p}] = findClu(reshape(hVP(c,s,ft,:,p),1,nSmpl),tStatP(:,p)',lenMin); % identify clusters in observed data
                        if isempty(sClP{c,s,ft,p})==1
                            sClP{c,s,ft,p} = 0; % replace size of non-existent clusters with 0
                        end
                        sClPmax(c,s,ft,p) = max(sClP{c,s,ft,p}); % idintify largest cluster for each permutation
                    end
                    % compare observed data with value distribution of 
                    % perutation data
                    for cl = 1:nCl(c,s,ft) % once for each cluster in the observed data
                        clLaI{c,s,ft}(cl,:) = sClPmax(c,s,ft,:)>sCl{c,s,ft}(cl); % identify how many of the permutation clusters are larger than the observed clusters
                        pCl{c,s,ft}(cl) = sum(clLaI{c,s,ft}(cl,:),2)/nPerm; % sum up number of permutation clusters that are larger than observed cluster to generate p-value
                        hCl{c,s,ft}(cl) = pCl{c,s,ft}(cl)<0.3; % identify observed clusters that are larger than 95 % of permutation clusters
                    end
                    [vecSigCl(c,s,ft,:),vecSigClI(c,s,ft,:),nSigCl(c,s,ft),sSigCl{c,s,ft}] = findSigClu(reshape(vecCl(c,s,ft,:),1,nSmpl),hCl{c,s,ft}); % identify clusters in observed data
                end
            end
        end
        %% save data
        fileName = append(fileName,['_',repRate]); % append rep Rate to file name
        if omRespI==0
            fileName = append(fileName,['_',subStd]); % append Std position to file name
        end

        switch inclCtrl
            case 0
                save(fileName,'varS','hV','pV','ciV','statsV','cohensD_V','vecSigCl','vecSigClI','nSigCl','sSigCl','lenMin')
            case {1,2,3}
                save(fileName,'varS','hV','pV','ciV','statsV','anovaV','multcompV','cohensD_V','vecSigCl','vecSigClI','nSigCl','sSigCl','lenMin')
        end
    end
end
beep
toc