SubcorticalDevianceDetectionHumans/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

% select test and stat parameter
test = 1; % 1: t-test, 2: ranksum
statParam = 1; % 1: P2P, 2: RMS, 3: AUC

% define time windows to be analysed (in s)
% low
win1L = [0.0027,0.0092]; % ABR (whole)
win2L = [0.005,0.0097]; % 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.0027,0.0092]; % ABR (whole)
win2H = [0.004,0.008]; % ABR (peak 5 only)
% win3H = [0.0045,0.007]; % ABR (late)
% win4H = [0.0065,0.01]; % 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 = 2: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]*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]*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 (trapz, mean V, etc.) for each sample (points or windows)
        
        ampS = cell(1,nComb); % preallocate
        ampSNd = cell(1,nComb); % preallocate
        maxV = cell(1,nComb); % preallocate
        maxI = cell(1,nComb); % preallocate
        timS = cell(1,nComb); % preallocate
        widS = cell(1,nComb); % preallocate
        timSNd = cell(1,nComb); % preallocate
        ampCsS = cell(1,nComb); % preallocate
        maxCsV = cell(1,nComb); % preallocate
        maxCsI = cell(1,nComb); % preallocate
        timCsS = cell(1,nComb); % preallocate
        ampAvCsS = cell(1,nComb); % preallocate
        ampSeCsS = cell(1,nComb); % preallocate
        timAvCsS = cell(1,nComb); % preallocate
        timSeCsS = cell(1,nComb); % preallocate
        ampGrAvCsS = 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};
            ampS{c} = zeros(nFiles,nStims,nCond,nFilt,nSmpl); % preallocate
            ampSNd{c} = zeros(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
            widS{c} = zeros(nFiles,nStims,nCond,nFilt,nSmpl); % preallocate
            timSNd{c} = zeros(nStims,nCond,nFilt,nSmpl); % preallocate
            widSNd{c} = zeros(nStims,nCond,nFilt,nSmpl); % preallocate
            ampCsS{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
            ampAvCsS{c} = zeros(nStims,nFilt,nSmpl,nConS); % preallocate
            ampSeCsS{c} = zeros(nStims,nFilt,nSmpl,nConS); % preallocate
            timAvCsS{c} = zeros(nStims,nFilt,nSmpl,nConS); % preallocate
            timSeCsS{c} = zeros(nStims,nFilt,nSmpl,nConS); % preallocate
            ampGrAvCsS{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
            pks = zeros(1,nFiles);
            locs = zeros(1,nFiles);
            w = zeros(1,nFiles);
            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
                                    switch statParam
                                        case 1
                                            ampS{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
                                        case 2
                                            ampS{c}(:,s,cn,ft,sp) = rms(data{c}{(s-1)*nCond+cn}((winAll{s}(sp,1):winAll{s}(sp,2)),:,ft)); % calculate variable per window
                                        case 3
                                            ampS{c}(:,s,cn,ft,sp) = trapz(abs(data{c}{(s-1)*nCond+cn}((winAll{s}(sp,1):winAll{s}(sp,2)),:,ft))); % calculate variable per window
                                    end
                                    % peak latency and width analysis
                                    for f = 1:nFiles
                                        [pksTemp,locTemp,~,~] = findpeaks(data{c}{(s-1)*nCond+cn}((winAll{s}(sp,1):winAll{s}(sp,2)),f,ft));
                                    [maxV,maxI] = max(pksTemp);
                                    pks(f) = pksTemp(maxI);
                                    locs(f) = locTemp(maxI);
                                    end
                                    timS{c}(:,s,cn,ft,sp) = (locs+(winAll{s}(sp,1)-pts2begin-round(blcEnd(s)*fs+1)))/fs*1000; % convert points to timing
                                    % [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

                                    % test for normality
                                    ampSNd{c}(s,cn,ft,sp) = lillietest(ampS{c}(:,s,cn,ft,sp));
                                    timSNd{c}(s,cn,ft,sp) = lillietest(timS{c}(:,s,cn,ft,sp));

                                    % standard cluster analysis
                                    if cn==1 % only 1 condition exists
                                        for cs = 1:nConS % once for each consecutive std
                                            switch statParam
                                                case 1
                                                    ampCsS{c}(:,s,ft,sp,cs) = peak2peak(conSAv{c}((winAll{s}(sp,1):winAll{s}(sp,2)),s,:,cs,ft)); % calculate variable per window
                                                case 2
                                                    ampCsS{c}(:,s,ft,sp,cs) = rms(conSAv{c}((winAll{s}(sp,1):winAll{s}(sp,2)),s,:,cs,ft)); % calculate variable per window
                                                case 3
                                                    ampCsS{c}(:,s,ft,sp,cs) = trapz(abs(conSAv{c}((winAll{s}(sp,1):winAll{s}(sp,2)),s,:,cs,ft))); % calculate variable per window
                                            end

                                            % peak latency and width analysis
                                            for f = 1:nFiles
                                                [pksTemp,locTemp,~,~] = findpeaks(conSAv{c}((winAll{s}(sp,1):winAll{s}(sp,2)),s,f,cs,ft));
                                                [maxV,maxI] = max(pksTemp);
                                                pks(f) = pksTemp(maxI);
                                                locs(f) = locTemp(maxI);
                                            end  
                                            % [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) = (locs+(winAll{s}(sp,1)-pts2begin-round(blcEnd(s)*fs+1)))/fs*1000; % convert points to timing
                                            % timCsS{c}(:,s,ft,sp,cs) = maxCsI{c}(:,s,ft,sp,cs)/fs*1000; % convert points to timing
                                            ampAvCsS{c}(s,ft,sp,cs) = mean(ampCsS{c}(:,s,ft,sp,cs),1); % calculate average from values of individual responses
                                            ampStdDCsS = std(ampCsS{c}(:,s,ft,sp,cs),0,1);
                                            ampSeCsS{c}(s,ft,sp,cs) = ampStdDCsS/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
                                            switch statParam
                                                case 1
                                                    ampGrAvCsS{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
                                                case 2
                                                    ampGrAvCsS{c}(s,ft,sp,cs) = rms(conSGrAv{c}((winAll{s}(sp,1):winAll{s}(sp,2)),s,cs,ft)); % calculate value from grand averaged response
                                                case 3
                                                    ampGrAvCsS{c}(s,ft,sp,cs) = trapz(abs(conSGrAv{c}((winAll{s}(sp,1):winAll{s}(sp,2)),s,cs,ft))); % calculate value from grand averaged response
                                            end
                                            [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
                                    ampS{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
        
        effS = zeros(nComb,nStims,nFilt,nSmpl,3,3); % preallocate
        switch inclCtrl
            case 0
                hV = zeros(nComb,nStims,nFilt,nSmpl,3); % preallocate
                pV = zeros(nComb,nStims,nFilt,nSmpl,3); % preallocate
                ciV = zeros(nComb,nStims,nFilt,nSmpl,3,2); % preallocate
                statsV = cell(nComb,nStims,nFilt,nSmpl,3); % preallocate
            case {1,2,3}
                anovaV = cell(nComb,nStims,nFilt,nSmpl,3); % preallocate
                multcompV = cell(nComb,nStims,nFilt,nSmpl,3); % 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
                        for p = 1:2 % once for each parameter (amplitude, latency, width)
                            switch p
                                case 1
                                    varS = ampS;
                                case 2
                                    varS = timS;
                                case 3
                                    varS = widS;
                            end
                            effS(c,s,ft,sp,p,1) = meanEffectSize(varS{c}(:,s,1,ft,sp),varS{c}(:,s,2,ft,sp),'Effect','cohen').Effect; % calculate effect size; dev vs. std
                            switch test
                                case 1
                                    [hV(c,s,ft,sp,p),pV(c,s,ft,sp,p),ciV(c,s,ft,sp,p,:),statsV{c,s,ft,sp,p}] = ttest(varS{c}(:,s,1,ft,sp),varS{c}(:,s,2,ft,sp)); % run t-test
                                case 2
                                    [pV(c,s,ft,sp,p),hV(c,s,ft,sp,p),statsV{c,s,ft,sp,p}] = ranksum(varS{c}(:,s,1,ft,sp),varS{c}(:,s,2,ft,sp)); % run ranksum test
                            end
                        end
                        switch inclCtrl
                            case {1,2,3} % if control is included
                                t = table(varS{c}(:,s,1,ft,sp,p),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,p} = ranova(rm); % run ANOVA
                                multcompV{c,s,ft,sp,p} = multcompare(rm,'stim_prob','ComparisonType','bonferroni'); % run multiple comparison to find differences between responses
                                effS(c,s,ft,sp,p,2) = meanEffectSize(varS{c}(:,s,1,ft,sp),varS{c}(:,s,3,ft,sp),'Effect','cohen').Effect; % calculate effect size; dev vs. ctrl
                                effS(c,s,ft,sp,p,3) = meanEffectSize(varS{c}(:,s,3,ft,sp),varS{c}(:,s,2,ft,sp),'Effect','cohen').Effect; % 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','ampS','timS','widS','ampSNd','hV','pV','ciV','statsV','effS',...
                    'ampCsS','ampAvCsS','ampSeCsS','timAvCsS','timSeCsS','ampGrAvCsS','timGrAvCsS','timSNd','widSNd')
            case {1,2,3}
                save(fileName,'winAll','winTxt','varS','hV','pV','ciV','statsV','anovaV','multcompV','effS')
        end
    end
end