%% 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