% Culster based permutation testing tic %% define some variables clear inclCtrl = 0; % include control data? (0: no, 1: 50Per, 2: MR, 3: DevOnly) omRespI = 1; % 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 % cluster settings lenMin = 5; %% 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']) load(['DD_avData_50Per_',repRate,'_',subStd,'_perm.mat']) dataAv_cell_Ctrl = dataAv_cell; % rename dataAvP_cell_Ctrl = dataAvP_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']) load(['DD_avData_MR_',repRate,'_',subStd,'_perm.mat']) dataAv_cell_Ctrl = dataAv_cell; % rename dataAvP_cell_Ctrl = dataAvP_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']) load(['DD_avData_LowOnly_',repRate,'_',subStd,'_perm.mat']) dataAv_cell_LowOnly = dataAv_cell; % rename dataAvP_cell_LowOnly = dataAvP_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']) load(['DD_avData_HighOnly_',repRate,'_',subStd,'_perm.mat']) dataAv_cell_HighOnly = dataAv_cell; % rename dataAvP_cell_HighOnly = dataAvP_cell; % rename end % load oddball data load(['DD_avData_Oddball_',repRate,'_',subStd,'.mat']) load(['DD_avData_Oddball_',repRate,'_',subStd,'_perm.mat']) case 1 inclCtrl = 0; % no control for OM data % load OM data load(['DD_avData_HighOm_',repRate,'.mat']) load(['DD_avData_HighOm_',repRate,'_perm.mat']) dataAv_cell_HighOm = dataAv_cell; % rename dataAvP_cell_HighOm = dataAvP_cell; % rename load(['DD_avData_LowOm_',repRate,'.mat']) load(['DD_avData_LowOm_',repRate,'_perm.mat']) dataAv_cell_LowOm = dataAv_cell; % rename dataAvP_cell_LowOm = dataAvP_cell; % rename end %% do some additional calulations that are required nSmpl = size(dataAv_cell{1}{1},1); % number of points of each recording nStims = 2; % number of different stimuli per combination nFilt = nFilt_cell{1}; switch inclCtrl case 0 switch omRespI case 0 fileName = 'DD_CbpData_noCtrl'; case 1 fileName = 'DD_CbpData_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}; dataAvP_cell{c}{1} = dataAvP_cell_HighOm{c}{3}; dataAvP_cell{c}{2} = dataAvP_cell_HighOm{c}{4}; dataAvP_cell{c}{3} = dataAvP_cell_LowOm{c}{1}; dataAvP_cell{c}{4} = dataAvP_cell_LowOm{c}{2}; end end nCond = 2; % 2 conditions: dev and std data = dataAv_cell; % only oddball data dataP = dataAvP_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)); dataP = cell(size(dataAvP_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 dataAvP_cell_Ctrl{c}{1} = dataAvP_cell_LowOnly{c}{1}; % assign DevOnly Low data as control of stimulus 1 dataAvP_cell_Ctrl{c}{2} = dataAvP_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) % observed data 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 % permutation data for p = 1:nPerm switch cn % change dataset depending on condition (Oddball vs. control) case {1,2} % if oddball dataP{c}{(s-1)*nCond+cn,p} = dataAvP_cell{c}{(s-1)*2+cn}(:,fileI_Oddb_cell_Ctrl{c},:,p); % concatenate oddball data of respective stimulus & condition (dev/std) case 3 % if control dataP{c}{(s-1)*nCond+cn,p} = dataAvP_cell_Ctrl{c}{s}(:,fileI_Ctrl_cell_Ctrl{c},:,p); % concatenate control data of respective stimulus end end end end end end %% calulate statistical variable (RMS, mean V, etc.) for each sample (points or windows) varS = cell(1,nComb); % preallocate varSP = 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 varSP{c} = zeros(nFiles,nStims,nCond,nFilt,nSmpl,nPerm); % 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) % observed data varS{c}(:,s,cn,ft,sp) = (data{c}{(s-1)*nCond+cn}(sp,:,ft)); % calculate variable per point % permutation data for p = 1:nPerm varSP{c}(:,s,cn,ft,sp,p) = (dataP{c}{(s-1)*nCond+cn}(sp,:,ft,p)); % 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 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 hVP = zeros(nComb,nStims,nFilt,nSmpl,nPerm); % preallocate pVP = zeros(nComb,nStims,nFilt,nSmpl,nPerm); % preallocate ciVP = zeros(nComb,nStims,nFilt,nSmpl,2,nPerm); % preallocate statsVP = cell(nComb,nStims,nFilt,nSmpl,nPerm); % preallocate case {1,2,3} anovaV = cell(nComb,nStims,nFilt,nSmpl); % preallocate multcompV = cell(nComb,nStims,nFilt,nSmpl); % preallocate anovaVP = cell(nComb,nStims,nFilt,nSmpl,nPerm); % preallocate multcompVP = cell(nComb,nStims,nFilt,nSmpl,nPerm); % 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 % observed data 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 % permutation data for p = 1:nPerm [hVP(c,s,ft,sp,p),pVP(c,s,ft,sp,p),ciVP(c,s,ft,sp,:,p),statsVP{c,s,ft,sp,p}] = ttest(varSP{c}(:,s,1,ft,sp,p),varSP{c}(:,s,2,ft,sp,p)); % run t-test end switch inclCtrl case {1,2,3} % if control is included % observed data 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 % permutation data for p = 1:nPerm anovaVP{c,s,ft,sp,p} = ranova(rm); % run ANOVA multcompVP{c,s,ft,sp,p} = multcompare(rm,'stim_prob','ComparisonType','bonferroni'); % run multiple comparison to find differences between responses end end end end end end %% cluster analysis tStat = zeros(1,nSmpl); % preallocate vecCl = zeros(nComb,nStims,nFilt,nSmpl); vecClI = zeros(nComb,nStims,nFilt,nSmpl); nCl = zeros(nComb,nStims,nFilt); sCl = cell(nComb,nStims,nFilt); tStatP = zeros(1,nSmpl,nPerm); % preallocate vecClP = zeros(nComb,nStims,nFilt,nSmpl); vecClIP = zeros(nComb,nStims,nFilt,nSmpl); nClP = zeros(nComb,nStims,nFilt); sClP = cell(nComb,nStims,nFilt); sClPmax = zeros(nComb,nStims,nFilt); clLaI = cell(nComb,nStims,nFilt); pCl = cell(nComb,nStims,nFilt); hCl = cell(nComb,nStims,nFilt); sigCl = zeros(nComb,nStims,nFilt,nSmpl); vecSigCl = zeros(nComb,nStims,nFilt,nSmpl); vecSigClI = zeros(nComb,nStims,nFilt,nSmpl); nSigCl = zeros(nComb,nStims,nFilt); sSigCl = cell(nComb,nStims,nFilt); 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 % observed data tStat(sp) = statsV{c,s,ft,sp}.tstat; % permutation data for p = 1:nPerm tStatP(sp,p) = statsVP{c,s,ft,sp,p}.tstat; end end % observed data [vecCl(c,s,ft,:),vecClI(c,s,ft,:),nCl(c,s,ft),sCl{c,s,ft}] = findClu(reshape(hV(c,s,ft,:),1,nSmpl),tStat,lenMin); % identify clusters in observed data % permutation data for p = 1:nPerm [vecClP(c,s,ft,:,p),vecClIP(c,s,ft,:),nClP(c,s,ft,p),sClP{c,s,ft,p}] = findClu(reshape(hVP(c,s,ft,:,p),1,nSmpl),tStatP(:,p)',lenMin); % identify clusters in observed data if isempty(sClP{c,s,ft,p})==1 sClP{c,s,ft,p} = 0; % replace size of non-existent clusters with 0 end sClPmax(c,s,ft,p) = max(sClP{c,s,ft,p}); % idintify largest cluster for each permutation end % compare observed data with value distribution of % perutation data for cl = 1:nCl(c,s,ft) % once for each cluster in the observed data clLaI{c,s,ft}(cl,:) = sClPmax(c,s,ft,:)>sCl{c,s,ft}(cl); % identify how many of the permutation clusters are larger than the observed clusters pCl{c,s,ft}(cl) = sum(clLaI{c,s,ft}(cl,:),2)/nPerm; % sum up number of permutation clusters that are larger than observed cluster to generate p-value hCl{c,s,ft}(cl) = pCl{c,s,ft}(cl)<0.3; % identify observed clusters that are larger than 95 % of permutation clusters end [vecSigCl(c,s,ft,:),vecSigClI(c,s,ft,:),nSigCl(c,s,ft),sSigCl{c,s,ft}] = findSigClu(reshape(vecCl(c,s,ft,:),1,nSmpl),hCl{c,s,ft}); % identify clusters in observed data 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,'varS','hV','pV','ciV','statsV','cohensD_V','vecSigCl','vecSigClI','nSigCl','sSigCl','lenMin') case {1,2,3} save(fileName,'varS','hV','pV','ciV','statsV','anovaV','multcompV','cohensD_V','vecSigCl','vecSigClI','nSigCl','sSigCl','lenMin') end end end beep toc