DevianceDetectionHumanABR/Repetition Suppression_Humans/RS_Stats.m

%% define some variables
clear

% define time windows to be analysed (in s)
% low
win1L = [0.000,0.01]; % ABR (whole)
win2L = [0.0055,0.009]; % ABR (peak 5 only)
win3L = [0,0.023]; % ABR (late)
win4L = [0.025,0.035]; % late component (overlaid with ABR in fast data set)
% high
win1H = [0,0.008]; % ABR (whole)
win2H = [0.0045,0.008]; % ABR (peak 5 only)
win3H = [0,0.023]; % ABR (late)
win4H = [0.025,0.035]; % late component (overlaid with ABR in fast data set)


winTxt = ["ABR","ABR early","ABR late","Omission Response"];
%% run stats on suppression / entrainment effects within blocks and correlation between ABRs and slow waves

for daSe = 1:2 % run stats for each SOA dataset
    % load data and define filename
    switch daSe
        case 1
            load('RS_avData_25ms.mat')
            filename = 'RS_statsData_25ms';
            int = [1,25]; % intervals for curve fitting (Rep# and time)
        case 2
            load('RS_avData_50ms.mat')
            filename = 'RS_statsData_50ms';
            int = [1,50]; % intervals for curve fitting (Rep# and time)
    end
    data = dataAv_cell; % rename
    
    winAllSti = cell(2,2);
    varS = cell(2,2); % preallocate
    varAvS = cell(2,2); % preallocate
    varStdDS = cell(2,2); % preallocate
    varSeS = cell(2,2); % preallocate
    meanV = cell(2,2); % preallocate
    meanGrAv = cell(2,2); % preallocate
    maxV = cell(2,2); % preallocate
    maxI = cell(2,2); % preallocate
    timS = cell(2,2); % preallocate
    timAvS = cell(2,2); % preallocate
    timStdDS = cell(2,2); % preallocate
    timSeS = cell(2,2); % preallocate
    for s = 1:2 % run once for each stimulus (low/high)
        for f = 1:2 % once for each stimulation type (block-stim and o1st-stim)
            % do some additional calulations that are required
            c = (f-1)*2+s; % index of dataset to be used
            fs = fs_cell{c};
            nResp = nResp_cell{c};
            pSiResp = pSiResp_cell{c};
            pts2begin = pts2begin_cell{c};
            winAll{1} = round([win1L;win2L;win3L;win4L]*fs)+pts2begin+round(blcEnd(s)*fs+1); % combine time windows, convert them into points and correct for stim delay (pts2begin)
            winAll{2} = round([win1H;win2H;win3H;win4H]*fs)+pts2begin+round(blcEnd(s)*fs+1);
            nWin = size(winAll{1},1);
            winAllSti{c} = zeros(nWin,nResp,2); % preallocate
            for r = 1:nResp
                winAllSti{c}(:,r,:) = [(r-1)*pSiResp+winAll{s}(:,1),(r-1)*pSiResp+winAll{s}(:,2)]; % concatenate windows for each stimulus of the block
            end
            nFilt = nFilt_cell{1};
            % calulate RMS and mean values within defined time windows
            nFiles = nFiles_cell{c};
            varS{s,f} = zeros(nFiles,nFilt,nResp,nWin); % preallocate
            varAvS{s,f} = zeros(nFilt,nResp,nWin); % preallocate
            varStdDS{s,f} = zeros(nFilt,nResp,nWin); % preallocate
            varSeS{s,f} = zeros(nFilt,nResp,nWin); % preallocate
            maxV{s,f} = zeros(nFiles,nFilt,nResp,nWin); % preallocate
            maxI{s,f} = zeros(nFiles,nFilt,nResp,nWin); % preallocate
            timS{s,f} = zeros(nFiles,nFilt,nResp,nWin); % preallocate
            timAvS{s,f} = zeros(nFilt,nResp,nWin); % preallocate
            timStdDS{s,f} = zeros(nFilt,nResp,nWin); % preallocate
            timSeS{s,f} = zeros(nFilt,nResp,nWin); % preallocate
            for ft = 1:nFilt % once for each filter
                for r = 1:nResp
                    for w = 1:nWin % once for each time window
                        % amplitude
                        varS{s,f}(:,ft,r,w) = peak2peak(data{c}(winAllSti{c}(w,r,1):winAllSti{c}(w,r,2),:,ft)); % calculate stats value
                        varAvS{s,f}(ft,r,w) = mean(varS{c}(:,ft,r,w),1); % calculate grand average of the stats value
                        varStdDS{s,f}(ft,r,w) = std(varS{c}(:,ft,r,w),0,1); % calculate standard deviation
                        varSeS{s,f}(ft,r,w) = varStdDS{s,f}(ft,r,w)/sqrt(nFiles); % calculate std error
                        % timing
                        [maxV{s,f}(:,ft,r,w),maxI{s,f}(:,ft,r,w)] = max(data{c}(winAllSti{c}(w,r,1):winAllSti{c}(w,r,2),:,ft)); % identify peak and peak index (=timing) per window
                        timS{s,f}(:,ft,r,w) = (maxI{s,f}(:,ft,r,w)+(winAllSti{c}(w,1,1)-pts2begin-round(blcEnd(s)*fs+1)))/fs*1000; % convert points to timing
                        % timS{s,f}(:,ft,r,w) = (maxI{s,f}(:,ft,r,w))/fs*1000; % convert points to timing
                        timAvS{s,f}(ft,r,w) = mean(timS{c}(:,ft,r,w),1); % calculate grand average of the stats value
                        timStdDS{s,f}(ft,r,w) = std(timS{c}(:,ft,r,w),0,1); % calculate standard deviation
                        timSeS{s,f}(ft,r,w) = timStdDS{s,f}(ft,r,w)/sqrt(nFiles); % calculate std error
                        % other stuff
                        meanV{s,f}(:,ft,r,w) = mean(data{c}(winAllSti{c}(w,r,1):winAllSti{c}(w,r,2),:,ft)); % calculate mean
                        meanGrAv{s,f}(ft,r,w) = mean(meanV{c}(:,ft,r,w),1); % calculate grand average of the mean
                    end
                end
            end
        end
    end
    
    % run t-test on each time window
    cohensD_V = zeros(nComb,nFilt,nResp,nWin); % preallocate
    hV = zeros(nComb,nFilt,nResp,nWin); % preallocate
    pV = zeros(nComb,nFilt,nResp,nWin); % preallocate
    ciV = zeros(nComb,nFilt,nResp,nWin,2); % preallocate
    statsV = cell(nComb,nFilt,nResp,nWin); % preallocate
    stdD_rms = zeros(nComb,nFilt,nResp,nWin); % preallocate
    cohensD_V_Bvs1 = zeros(nComb,nFilt,2,nWin); % preallocate
    hV_Bvs1 = zeros(nComb,nFilt,2,nWin); % preallocate
    pV_Bvs1 = zeros(nComb,nFilt,2,nWin); % preallocate
    ciV_Bvs1 = zeros(nComb,nFilt,2,nWin,2); % preallocate
    statsV_Bvs1 = cell(nComb,nFilt,2,nWin); % preallocate
    stdD_rms_Bvs1 = zeros(nComb,nFilt,2,nWin); % preallocate
    ciV_1st = zeros(2,nFilt,nWin,2,2); % preallocate
    fV = cell(2,2,nFilt,nWin); % preallocate
    gofV = cell(2,2,nFilt,nWin); % preallocate
    mdl = cell(2,2,nFilt,nWin); % preallocate
    for s = 1:2 % run once for each stimulus (low/high)
        for f = 1:2 % once for each stimulation type (block-stim and o1st-stim)
            c = (f-1)*2+s; % index of dataset to be used
            nResp = nResp_cell{c};
            for ft = 1:nFilt % once for each filter
                for r = 1:nResp % once for each stimulus
                    for w = 1:nWin % once for each time window
                        [hV(c,ft,r,w),pV(c,ft,r,w),ciV(c,ft,r,w,:),statsV{c,ft,r,w}] = ttest(varS{c}(:,ft,1,w),varS{c}(:,ft,r,w)); % run t-test; the first one vs. every window
                        cohensD_V(c,ft,r,w) = CohensD(varS{c}(:,ft,1,w),varS{c}(:,ft,r,w)); % calculate Cohens D; the first one vs. every window
                        stdD_rms(c,ft,r,w) = std(varS{c}(:,ft,r,w));
                    end
                end
                % fit curve to RMS values of block-responses
                switch f
                    case 1
                        for w = 1:nWin
                            % fitTy = fittype('V*exp(t/tau)+b','independent','t'); % model for exponential growth (different to bat data)
                            % fitTy = fittype('L/(1+exp(-k*(t-x)))','independent','t'); % model for exponential growth (different to bat data)
                            fitTy = fittype('poly1'); % model
                            x = (1:nResp-1)';
                            for v = 1:2 % once for amplitude, once for timing
                                switch v
                                    case 1
                                        varTemp = varS;
                                    case 2
                                        varTemp = timS;
                                end
                                pd = fitdist(varTemp{s,f}(:,ft,1,w),'Normal'); % Fit a normal distribution object to the RMS data of the first response
                                ciV_1st(s,ft,w,:,:) = paramci(pd); % calculate confidence intervals for the RMS data of the first response
                                yMax = max(reshape(mean(varTemp{s,f}(:,ft,1:end-1,w),1),nResp-1,1)); % caluclate min of y (for normalisation)
                                % y = reshape(mean(varTemp{s,f}(:,ft,1:end-1,w),1),nResp-1,1)/yMax; % normalise y to max of y
                                y = reshape(mean(varTemp{s,f}(:,ft,1:end-1,w),1),nResp-1,1); % y without normalisation
                                [fV{v,s,ft,w},gofV{v,s,ft,w}] = fit(x,y,fitTy); % fit model to data (fexible model)
                                mdl{v,s,ft,w} = fitlm(x,y); % fit model to data (linear regression model)
                            end
                        end
                    case 2
                        % compare o1st response vs. first and last response of block stim
                        for r = 1:2 % once for first and once for last response of block stim
                            for w = 1:nWin % once for each time window
                                switch r
                                    case 1 % the first block response vs. the o1st response
                                        [hV_Bvs1(s,ft,r,w),pV_Bvs1(s,ft,r,w),ciV_Bvs1(s,ft,r,w,:),statsV_Bvs1{s,ft,r,w}] = ttest(varS{s,1}(:,ft,1,w),varS{s,2}(:,ft,1,w)); % run t-test
                                        cohensD_V_Bvs1(s,ft,r,w) = CohensD(varS{s,1}(:,ft,1,w),varS{s,2}(:,ft,1,w)); % calculate Cohens D
                                    case 2 % the last block response vs. the o1st response
                                        [hV_Bvs1(s,ft,r,w),pV_Bvs1(s,ft,r,w),ciV_Bvs1(s,ft,r,w,:),statsV_Bvs1{s,ft,r,w}] = ttest(varS{s,1}(:,ft,end,w),varS{s,2}(:,ft,1,w)); % run t-test
                                        cohensD_V_Bvs1(s,ft,r,w) = CohensD(varS{s,1}(:,ft,end,w),varS{s,2}(:,ft,1,w)); % calculate Cohens D
                                end
                            end
                        end
                end
            end
        end
    end
    
    for c = 1:nComb % once for each stimulus combination
        nResp = nResp_cell{c};
        % determine correlation between ABR-size and average value of the slow wave
        waveAbrSize = zeros(2,nResp,nComb,nWin);
        waveAbrNorm = zeros(2,nResp,nComb,nWin);
        RPea = zeros(nComb,nWin); % preallocate
        pRPea = zeros(nComb,nWin); % preallocate
        RSpea = zeros(nComb,nWin); % preallocate
        pRSpea = zeros(nComb,nWin); % preallocate
        for w = 1:nWin
            % individual averages
            abrSize = zeros(1,size(varS{c},1)*nResp); % preallocate
            waveSize = zeros(1,size(meanV{c},1)*nResp); % preallocate
            abrSize(:) = varS{c}(:,3,:,w); % RMS of the ABR response for the respective stimulus and window
            waveSize(:) = meanV{c}(:,5,:,w); % RMS of the ABR response for the respective stimulus and window
            abrNorm = (abrSize-min(abrSize))/(max(abrSize)-min(abrSize)); % normalise data between 0 and 1
            waveNorm = (waveSize-min(waveSize))/(max(waveSize)-min(waveSize)); % normalise data between 0 and 1
            % grand averages (work much better)
            abrSizeGrAv = varAvS{c}(3,:,w); % RMS of the ABR response for the respective stimulus and window
            abrNormGrAv = (abrSizeGrAv-min(abrSizeGrAv))/(max(abrSizeGrAv)-min(abrSizeGrAv)); % normalise data between 0 and 1
            waveSizeGrAv = meanGrAv{c}(5,:,w); % mean value of the slow wave for the respective stimulus and window
            waveNormGrAv = (waveSizeGrAv-min(waveSizeGrAv))/(max(waveSizeGrAv)-min(waveSizeGrAv)); % normalise data between 0 and 1
            waveAbrSize(:,:,c,w) = [abrSizeGrAv;waveSizeGrAv];
            waveAbrNorm(:,:,c,w) = [abrNormGrAv;waveNormGrAv];
            [RPea(c,w),pRPea(c,w)] = corr(waveAbrNorm(1,:,c,w)',waveAbrNorm(2,:,c,w)','Type','Pearson');
            [RSpea(c,w),pRSpea(c,w)] = corr(waveAbrNorm(1,:,c,w)',waveAbrNorm(2,:,c,w)','Type','Spearman');
        end
    end
    
    % save data
    save(filename,'int','winAllSti','winTxt','varS','varAvS','varSeS','meanV',...
        'timS','timAvS','timSeS','meanGrAv','hV','pV','ciV','statsV',...
        'stdD_rms','cohensD_V','ciV_1st','fV','gofV','hV_Bvs1',...
        'pV_Bvs1','ciV_Bvs1','statsV_Bvs1','stdD_rms_Bvs1','cohensD_V_Bvs1',...
        'waveAbrSize','waveAbrNorm','RPea','pRPea','RSpea','pRSpea','mdl')
end