SubcorticalDevianceDetectionHumans/Deviance Detection_Humans/DD_ConSAnalysis.m

%% define some variables

clear

% general parameters
fitI = 0; % should curve be fitted to datapoints?

% define time windows to be analysed (in s)
% low
win1L = [0.0027,0.0092]; % ABR (whole)
win2L = [0.005,0.0092]; % ABR (peak 5 only)

% high
win1H = [0.0027,0.0092]; % ABR (whole)
win2H = [0.004,0.008]; % ABR (peak 5 only)

% name time windows
winTxt = ["ABR","Wave V"];

% curve fitting parameters
% fitTy = fittype('poly2'); % define model
fitTy = fittype('V*exp(-t/tau)+b','independent','t'); % model for exponential decay

%% load and concatenate data

% load and rename variables
load("DD_avData_Oddball_25ms_SbD.mat")
dataAv_25 = conSAv;
se_25 = conSSe;
pts2begin_25 = pts2begin_cell{1}(1);
fs_25 = fs_cell{1}(1);
pSmpl_25 = pSmpl;
stimDel_25 = stimDelay_cell{1}(1);
load("DD_avData_Oddball_30ms_SbD.mat")
dataAv_30 = conSAv;
se_30 = conSSe;
pts2begin_30 = pts2begin_cell{1}(1);
fs_30 = fs_cell{1}(1);
pSmpl_30 = pSmpl;
stimDel_30 = stimDelay_cell{1}(1);
load("DD_avData_Oddball_50ms_SbD.mat")
dataAv_50 = conSAv;
se_50 = conSSe;
pts2begin_50 = pts2begin_cell{1}(1);
fs_50 = fs_cell{1}(1);
pSmpl_50 = pSmpl;
stimDel_50 = stimDelay_cell{1}(1);

% concatenate variables
data = cat(1,dataAv_25,dataAv_30,dataAv_50);
se = cat(1,se_25,se_30,se_50);
pts2begin = cat(1,pts2begin_25,pts2begin_30,pts2begin_50);
fs = cat(1,fs_25,fs_30,fs_50);
pSmpl = cat(1,pSmpl_25,pSmpl_30,pSmpl_50);
stimDel = cat(1,stimDel_25,stimDel_30,stimDel_50);

% calculate some variables
nStims = 2;
nRepRate = size(data,1);
nFilt = size(data{1},5);
nFiles = size(data{1},3);
if devConSI==1
    nConS = nConS+1;
end

%% do some additional calculations

% preallocate
winAll = cell(nRepRate,1);

for r = 1:nRepRate
    winAll{r,1} = round([win1L;win2L]*fs(r))+pts2begin(r)+round(blcEnd(1)*fs(r)+1); % combine time windows, convert them into points and correct for stim delay (pts2begin) and baseline correction window (blcEnd = stimulus peak)
    winAll{r,2} = round([win1H;win2H]*fs(r))+pts2begin(r)+round(blcEnd(2)*fs(r)+1);
end
nWin = size(winAll{1},1); % number of windows per stimulus
nStPa = 5; % number of statistical parameters

%% calculate statistical parameters

% preallocate
stPa = zeros(nFiles,nRepRate,nStims,nConS,nFilt,nWin,nStPa);

for r = 1:nRepRate % for each rep rate
    for s = 1:nStims % for each stimulus (A and B)
        for cs = 1:nConS
            for ft = 1:nFilt % for each filter
                for w = 1:nWin % for each window
                    for p = 1:nStPa % for each stat param
                        switch p
                            case 1 % peak 2 peak
                                stPa(:,r,s,cs,ft,w,p) = peak2peak(data{r}(winAll{r,s}(w,1):winAll{r,s}(w,2),s,:,cs,ft));
                            case 2 % RMS
                                stPa(:,r,s,cs,ft,w,p) = rms(data{r}(winAll{r,s}(w,1):winAll{r,s}(w,2),s,:,cs,ft));
                            case 3 % AUC
                                stPa(:,r,s,cs,ft,w,p) = trapz(abs(data{r}(winAll{r,s}(w,1):winAll{r,s}(w,2),s,:,cs,ft)));
                            case {4,5} % peak latency & peak width
                                pks = zeros(1,nFiles);
                                loc = zeros(1,nFiles);
                                wid = zeros(1,nFiles);
                                for f = 1:nFiles
                                    [pksTemp,locTemp,widTemp,~] = findpeaks(data{r}(winAll{r,s}(w,1):winAll{r,s}(w,2),s,f,cs,ft),"WidthReference",'halfprom');
                                    % if r==1&&s==1&&cs==4&&ft==2&&w==1&&p==4
                                    %     figure; findpeaks(data{r}(winAll{r,s}(w,1):winAll{r,s}(w,2),s,f,cs,ft),"WidthReference",'halfprom')
                                    % end
                                    [maxV,maxI] = max(pksTemp);
                                    pks(f) = pksTemp(maxI);
                                    loc(f) = locTemp(maxI);
                                    wid(f) = widTemp(maxI);
                                end
                                switch p
                                    case 4
                                        stPa(:,r,s,cs,ft,w,p) = (loc+(winAll{r,s}(w,1)-pts2begin(r)-round(blcEnd(s)*fs(r)+1)))/fs(r)*1000; % convert points to time
                                    case 5
                                        stPa(:,r,s,cs,ft,w,p) = wid/fs(r)*1000; % convert points to time
                                end
                        end
                    end
                end
            end
        end
    end
end
%% test plotting

% A = mean(stPa,1);
% B = reshape(A(1,3,1,:,3,2,4),1,nConS);
% figure; plot(B,'o','Linewidth',2,'Color','k')

%% fit curve to data

% preallocate
ci1st = cell(nRepRate,nStims,nFilt,nWin,nStPa);
flexFit = cell(nRepRate,nStims,nFilt,nWin,nStPa);
flexGof = cell(nRepRate,nStims,nFilt,nWin,nStPa);
linFit = cell(nRepRate,nStims,nFilt,nWin,nStPa);

if fitI==1
    x = (1:nConS)';
    for r = 1:nRepRate % for each rep rate
        for s = 1:nStims % for each stimulus (A and B)
            for cs = 1:nConS % for each standard in the cluster
                for ft = 1:nFilt % for each filter
                    for w = 1:nWin % for each window
                        for p = 1:nStPa % for each stat param
                            pd = fitdist(stPa(:,r,s,1,ft,w,p),'Normal'); % Fit a normal distribution object to the data of the first response
                            ci1st{r,s,ft,w,p} = paramci(pd); % calculate confidence intervals for the data of the first response
                            y = reshape(mean(stPa(:,r,s,:,ft,w,p),1),nConS,1); % prepare data for curve fitting
                            [flexFit{r,s,ft,w,p},flexGof{r,s,ft,w,p}] = fit(x,y,fitTy,'StartPoint',[1,1,1],'Upper',[10000,10000,10000],'Lower',[0,0,0]); % fit model to data (fexible model)
                            linFit{r,s,ft,w,p} = fitlm(x,y); % fit model to data (linear regression model)
                        end
                    end
                end
            end
        end
    end
end

%% prepare variables for ANOVA

% preallocate
stPaAno = zeros((nFiles*nRepRate*nConS),nStims,nFilt,nWin,nStPa);
faRole = strings((nFiles*nRepRate*nConS),nStims,nFilt,nWin,nStPa);
faRepRate = strings((nFiles*nRepRate*nConS),nStims,nFilt,nWin,nStPa);

count = 0;

for r = 1:nRepRate % for each rep rate
    switch r
        case 1
            RepRateName = "25ms";
        case 2
            RepRateName = "30ms";
        case 3
            RepRateName = "50ms";
    end
    for ro = 1:nConS % for each role (dev, std1, std2, std3)
        switch ro
            case 1
                roleName = "Dev";
            case 2
                roleName = "Std1";
            case 3
                roleName = "Std2";
            case 4
                roleName = "Std3";
        end
        count = count+1;
        for s = 1:nStims % for each stimulus
            for ft = 1:nFilt % for each filter
                for w = 1:nWin % for each window
                    for p = 1:nStPa % for each stat param
                        stPaAno((nFiles*(count-1))+1:(nFiles*(count-1))+nFiles,s,ft,w,p) = stPa(:,r,s,ro,ft,w,p); % reorganise stPa to have all conditions that should be tested in anova in one dimension
                        % create factors for anova
                        faRole((nFiles*(count-1))+1:(nFiles*(count-1))+nFiles,s,ft,w,p) = roleName;
                        faRepRate((nFiles*(count-1))+1:(nFiles*(count-1))+nFiles,s,ft,w,p) = RepRateName;
                    end
                end
            end
        end
    end
end

%% run ANOVA & multiple comparison test (without correction)

% preallocate
anoP = cell(nStims,nFilt,nWin,nStPa);
anoTbl = cell(nStims,nFilt,nWin,nStPa);
anoStats = cell(nStims,nFilt,nWin,nStPa);
anoTerms = cell(nStims,nFilt,nWin,nStPa);
mcoC = cell(nStims,nFilt,nWin,nStPa);
mcoM = cell(nStims,nFilt,nWin,nStPa);
mcoH = cell(nStims,nFilt,nWin,nStPa);
mcoGnames = cell(nStims,nFilt,nWin,nStPa);
mcoTbl = cell(nStims,nFilt,nWin,nStPa);

for s = 1:nStims % for each stimulus
    for ft = 1:nFilt % for each filter
        for w = 1:nWin % for each window
            for p = 1:nStPa % for each stat param
                y = stPaAno(:,s,ft,w,p); % data to be tested
                group = {faRole(:,s,ft,w,p),faRepRate(:,s,ft,w,p)}; % combine factors in group variable
                [anoP{s,ft,w,p},anoTbl{s,ft,w,p},anoStats{s,ft,w,p},anoTerms{s,ft,w,p}]  = anovan(y,group,'model','interaction','varnames',{'faRole','faRepRate'},'display','off');
                [mcoC{s,ft,w,p},mcoM{s,ft,w,p},mcoH{s,ft,w,p},mcoGnames{s,ft,w,p}] = multcompare(anoStats{s,ft,w,p},'Dimension',[1,2],"CriticalValueType","lsd",'display','off');
                mcoTbl{s,ft,w,p} = array2table(mcoC{s,ft,w,p},"VariableNames",["Group A","Group B","Lower Limit","AtoB","Upper Limit","Pvalue"]); % not used. Paird t-tests are used instead (because data is paired for relevant post hoc comparisons)
                mcoTbl{s,ft,w,p}.("Group A") = mcoGnames{s,ft,w,p}(mcoTbl{s,ft,w,p}.("Group A"));
                mcoTbl{s,ft,w,p}.("Group B") = mcoGnames{s,ft,w,p}(mcoTbl{s,ft,w,p}.("Group B"));
            end
        end
    end
end

%% compare all responses to the first one with t-test

% preallocate
conSH = zeros(nRepRate,nStims,nConS-1,nFilt,nWin,nStPa);
conSP = zeros(nRepRate,nStims,nConS-1,nFilt,nWin,nStPa);
conSCi = cell(nRepRate,nStims,nConS-1,nFilt,nWin,nStPa);
conSStats = cell(nRepRate,nStims,nConS-1,nFilt,nWin,nStPa);
conSEffS = zeros(nRepRate,nStims,nConS-1,nFilt,nWin,nStPa);

for r = 1:nRepRate % for each rep rate
    for s = 1:nStims % for each stimulus (A and B)
        for cs = 2:nConS % for each standard in the cluster (except for the first)
            for ft = 1:nFilt % for each filter
                for w = 1:nWin % for each window
                    for p = 1:nStPa % for each stat param
                        [conSH(r,s,cs-1,ft,w,p),conSP(r,s,cs-1,ft,w,p),conSCi{r,s,cs-1,ft,w,p},conSStats{r,s,cs-1,ft,w,p}] = ttest(stPa(:,r,s,1,ft,w,p),stPa(:,r,s,cs,ft,w,p)); % run t-test
                        conSEffS(r,s,cs-1,ft,w,p) = meanEffectSize(stPa(:,r,s,1,ft,w,p),stPa(:,r,s,cs,ft,w,p),'Effect','cohen').Effect; % calculate effect size; dev vs. std
                    end
                end
            end
        end
    end
end


%% save data

save('DD_ConSData','fitTy','stPa','ci1st','flexFit','flexGof','linFit',...
    'pts2begin','fs','pSmpl','stimDel','nStims','nRepRate','nFilt','nFiles','nConS',...
    'nWin','nStPa','winAll','conSH','conSP','conSCi','conSStats',...
    'conSEffS','combName','data','se','blcEnd','anoTbl','mcoTbl')