DevianceDetectionHumanABR/Deviance Detection_Humans/DD_Stats.m

%% define some variables
clear

% select dataset
daSeI = 1; % 1: only Low/High; 2: 3/8 and Low/High
anaType = 1; % 1: window-based analysis; 2: point-based analysis

% define time windows to be analysed (in s)
% low
win1L = [0.00,0.01]; % ABR (whole)
win2L = [0.0055,0.009]; % ABR (peak 5 only)
win3L = [0.0045,0.0075]; % ABR (late)
win4L = [0.006,0.010]; % ABR (omitted); original window ([0.007,0.017]) is too long for 25 ms dataset
% high
win1H = [0.00,0.008]; % ABR (whole)
win2H = [0.0045,0.008]; % ABR (peak 5 only)
win3H = [0.0045,0.007]; % ABR (late)
win4H = [0.0065,0.01]; % ABR (omitted)

% win1L = [0.000,0.0075]; % ABR (whole)
% win2L = [0.000,0.0045]; % ABR (early)
% win3L = [0.0045,0.0075]; % ABR (late)
% win4L = [0.000,0.007]; % ABR (omitted); original window ([0.007,0.017]) is too long for 25 ms dataset
% % high
% win1H = [0,0.007]; % ABR (whole)
% win2H = [0,0.0045]; % ABR (early)
% win3H = [0.0045,0.007]; % ABR (late)
% win4H = [0,0.007]; % ABR (omitted)


winTxt = ["ABR","ABR early","ABR late","Omission Response"];

inclCtrl = 0; % include control data? (0: no, 1: 50Per, 2: MR, 3: DevOnly)
omRespI = 0; % 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
%% 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'])
                        dataAv_cell_Ctrl = dataAv_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'])
                        dataAv_cell_Ctrl = dataAv_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'])
                        dataAv_cell_LowOnly = dataAv_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'])
                        dataAv_cell_HighOnly = dataAv_cell; % rename
                end
                % load oddball data
                load(['DD_avData_Oddball_',repRate,'_',subStd,'.mat'])
            case 1
                inclCtrl = 0; % no control for OM data
                % load OM data
                load(['DD_avData_HighOm_',repRate,'.mat'])
                dataAv_cell_HighOm = dataAv_cell; % rename
                load(['DD_avData_LowOm_',repRate,'.mat'])
                dataAv_cell_LowOm = dataAv_cell; % rename
        end
        %% do some additional calulations that are required
        
        fs = fs_cell{1}(1);
        pts2begin = pts2begin_cell{1}(1);
        winAll{1} = round([win1L;win2L;win3L;win4L]*fs)+pts2begin+round(blcEnd(1)*fs+1); % combine time windows, convert them into points and correct for stim delay (pts2begin) and baseline correction window (blcEnd = stimulus peak)
        winAll{2} = round([win1H;win2H;win3H;win4H]*fs)+pts2begin+round(blcEnd(2)*fs+1);
        switch anaType
            case 1 % window-based analysis
                nSmpl = size(winAll{1},1); % number of windows to be analysed
            case 2 % point-based analysis
                nSmpl = size(dataAv_cell{1}{1},1); % number of points of each recording
        end
        
        nStims = 2; % number of different stimuli per combination
        nFilt = nFilt_cell{1};
        
        switch inclCtrl
            case 0
                switch omRespI
                    case 0
                        fileName = 'DD_statsData_noCtrl';
                    case 1
                        fileName = 'DD_statsData_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};
                        end
                end
                nCond = 2; % 2 conditions: dev and std
                data = dataAv_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));
                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
                    end
                    for s = 1:nStims % once for each stimulus frequency
                        for cn = 1:nCond % once for each stimulus condition (dev, std, ctrl)
                            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
                        end
                    end
                end
        end
        %% calulate statistical variable (RMS, mean V, etc.) for each sample (points or windows)
        
        varS = cell(1,nComb); % preallocate
        maxV = cell(1,nComb); % preallocate
        maxI = cell(1,nComb); % preallocate
        timS = cell(1,nComb); % preallocate
        varCsS = cell(1,nComb); % preallocate
        maxCsV = cell(1,nComb); % preallocate
        maxCsI = cell(1,nComb); % preallocate
        timCsS = cell(1,nComb); % preallocate
        varAvCsS = cell(1,nComb); % preallocate
        varSeCsS = cell(1,nComb); % preallocate
        timAvCsS = cell(1,nComb); % preallocate
        timSeCsS = cell(1,nComb); % preallocate
        varGrAvCsS = cell(1,nComb); % preallocate
        maxGrAvCsV = cell(1,nComb); % preallocate
        maxGrAvCsI = cell(1,nComb); % preallocate
        timGrAvCsS = 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
            maxV{c} = zeros(nFiles,nStims,nCond,nFilt,nSmpl); % preallocate
            maxI{c} = zeros(nFiles,nStims,nCond,nFilt,nSmpl); % preallocate
            timS{c} = zeros(nFiles,nStims,nCond,nFilt,nSmpl); % preallocate
            varCsS{c} = zeros(nFiles,nStims,nFilt,nSmpl,nConS); % preallocate
            maxCsV{c} = zeros(nFiles,nStims,nFilt,nSmpl,nConS); % preallocate
            maxCsI{c} = zeros(nFiles,nStims,nFilt,nSmpl,nConS); % preallocate
            timCsS{c} = zeros(nFiles,nStims,nFilt,nSmpl,nConS); % preallocate
            varAvCsS{c} = zeros(nStims,nFilt,nSmpl,nConS); % preallocate
            varSeCsS{c} = zeros(nStims,nFilt,nSmpl,nConS); % preallocate
            timAvCsS{c} = zeros(nStims,nFilt,nSmpl,nConS); % preallocate
            timSeCsS{c} = zeros(nStims,nFilt,nSmpl,nConS); % preallocate
            varGrAvCsS{c} = zeros(nStims,nFilt,nSmpl,nConS); % preallocate
            maxGrAvCsV{c} = zeros(nStims,nFilt,nSmpl,nConS); % preallocate
            maxGrAvCsI{c} = zeros(nStims,nFilt,nSmpl,nConS); % preallocate
            timGrAvCsS{c} = zeros(nStims,nFilt,nSmpl,nConS); % 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)
                            switch anaType
                                case 1
                                    varS{c}(:,s,cn,ft,sp) = peak2peak(data{c}{(s-1)*nCond+cn}((winAll{s}(sp,1):winAll{s}(sp,2)),:,ft)); % calculate variable per window
                                    [maxV{c}(:,s,cn,ft,sp),maxI{c}(:,s,cn,ft,sp)] = max(data{c}{(s-1)*nCond+cn}((winAll{s}(sp,1):winAll{s}(sp,2)),:,ft)); % identify peak and peak index (=timing) per window
                                    timS{c}(:,s,cn,ft,sp) = (maxI{c}(:,s,cn,ft,sp)+(winAll{s}(sp,1)-pts2begin-round(blcEnd(s)*fs+1)))/fs*1000; % convert points to timing
                                    % standard cluster analysis
                                    if cn==1 % only 1 condition exists
                                        for cs = 1:nConS % once for each consecutive std
                                            varCsS{c}(:,s,ft,sp,cs) = peak2peak(conSAv{c}((winAll{s}(sp,1):winAll{s}(sp,2)),s,:,cs,ft)); % calculate variable per window
                                            [maxCsV{c}(:,s,ft,sp,cs),maxCsI{c}(:,s,ft,sp,cs)] = max(conSAv{c}((winAll{s}(sp,1):winAll{s}(sp,2)),s,:,cs,ft)); % identify peak and peak index (=timing) per window
                                            timCsS{c}(:,s,ft,sp,cs) = maxCsI{c}(:,s,ft,sp,cs)/fs*1000; % convert points to timing
                                            varAvCsS{c}(s,ft,sp,cs) = mean(varCsS{c}(:,s,ft,sp,cs),1); % calculate average from values of individual responses
                                            varStdDCsS = std(varCsS{c}(:,s,ft,sp,cs),0,1);
                                            varSeCsS{c}(s,ft,sp,cs) = varStdDCsS/sqrt(nFiles);
                                            timAvCsS{c}(s,ft,sp,cs) = mean(timCsS{c}(:,s,ft,sp,cs),1);
                                            timStdDCsS = std(timCsS{c}(:,s,ft,sp,cs),0,1);
                                            timSeCsS{c}(s,ft,sp,cs) = timStdDCsS/sqrt(nFiles);
                                            % additionally
                                            varGrAvCsS{c}(s,ft,sp,cs) = peak2peak(conSGrAv{c}((winAll{s}(sp,1):winAll{s}(sp,2)),s,cs,ft)); % calculate value from grand averaged response
                                            [maxGrAvCsV{c}(s,ft,sp,cs),maxGrAvCsI{c}(s,ft,sp,cs)] = max(conSGrAv{c}((winAll{s}(sp,1):winAll{s}(sp,2)),s,cs,ft)); % identify peak and peak index (=timing) per window
                                            timGrAvCsS{c}(s,ft,sp,cs) = maxGrAvCsI{c}(s,ft,sp,cs)/fs*1000; % convert points to timing
                                        end
                                    end
                                case 2
                                    varS{c}(:,s,cn,ft,sp) = (data{c}{(s-1)*nCond+cn}(sp,:,ft)); % 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
        cohensDT_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
                hTV = zeros(nComb,nStims,nFilt,nSmpl); % preallocate
                pTV = zeros(nComb,nStims,nFilt,nSmpl); % preallocate
                ciTV = zeros(nComb,nStims,nFilt,nSmpl,2); % preallocate
                statsTV = cell(nComb,nStims,nFilt,nSmpl); % preallocate
            case {1,2,3}
                anovaV = cell(nComb,nStims,nFilt,nSmpl); % preallocate
                multcompV = cell(nComb,nStims,nFilt,nSmpl); % 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
                        % amplitude
                        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
                        % timing
                        cohensDT_V(c,s,ft,sp,1) = CohensD(timS{c}(:,s,2,ft,sp),timS{c}(:,s,1,ft,sp)); % calculate Cohens D; dev vs. std
                        [hTV(c,s,ft,sp),pTV(c,s,ft,sp),ciTV(c,s,ft,sp,:),statsTV{c,s,ft,sp}] = ttest(timS{c}(:,s,2,ft,sp),timS{c}(:,s,1,ft,sp)); % run t-test
                        switch inclCtrl
                            case {1,2,3} % if control is included
                                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
                        end
                    end
                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,'winAll','winTxt','varS','hV','pV','ciV','statsV','cohensD_V','timS','hTV','pTV','ciTV','statsTV','cohensDT_V',...
                    'varCsS','timCsS','varAvCsS','varSeCsS','timAvCsS','timSeCsS','varGrAvCsS','timGrAvCsS')
            case {1,2,3}
                save(fileName,'winAll','winTxt','varS','hV','pV','ciV','statsV','anovaV','multcompV','cohensD_V')
        end
    end
end