123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364 |
- currentFile = mfilename('fullpath');
- [pathstr,~,~] = fileparts(currentFile);
- cd(fullfile(pathstr,'..'))
- rootpath = pwd;
- pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
- pn.data_erp = fullfile(rootpath, 'data', 'erp');
- pn.data_erf = fullfile(rootpath, 'data', 'erf');
- pn.tools = fullfile(rootpath, 'tools');
- addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
- addpath(fullfile(pn.tools, 'BrewerMap'));
- addpath(fullfile(pn.tools, 'shadedErrorBar'));
- addpath(genpath(fullfile(pn.tools, 'RainCloudPlots')));
- pn.figures = fullfile(rootpath, 'figures');
- %% load erp_bl
- for ind_id = 1:33
- id = sprintf('sub-%03d', ind_id);
- load(fullfile(pn.data_erp, [id,'_erp_bl.mat']));
- for ind_option = 1:numel(conds.scene_category)
- if ind_id == 1
- erpgroup.scene_category.(conds.scene_category{ind_option}) = erp_bl.scene_category{ind_option};
- erpgroup.scene_category.(conds.scene_category{ind_option}) = ...
- rmfield(erpgroup.scene_category.(conds.scene_category{ind_option}), {'avg', 'var', 'dof'});
- erpgroup.scene_category.(conds.scene_category{ind_option}).dimord = 'sub_chan_time';
- end
- erpgroup.scene_category.(conds.scene_category{ind_option}).avg(ind_id,:,:) = erp_bl.scene_category{ind_option}.avg;
- end
- end
- time = erpgroup.scene_category.manmade.time;
- elec = erpgroup.scene_category.manmade.elec;
- channels = erpgroup.scene_category.manmade.label;
- mergeddata = cat(4, erpgroup.scene_category.manmade.avg, ...
- erpgroup.scene_category.natural.avg);
- %% plot the ERPs (for initial inspection)
- % manmade = erpgroup.scene_category.manmade;
- % manmade.avg = squeeze(nanmean(manmade.avg,1));
- % manmade.dimord = 'chan_time';
- %
- % natural = erpgroup.scene_category.natural;
- % natural.avg = squeeze(nanmean(natural.avg,1));
- % natural.dimord = 'chan_time';
- %
- % cfg = [];
- % cfg.layout = 'EEG1010.lay';
- % cfg.interactive = 'yes';
- % cfg.showoutline = 'yes';
- % cfg.xlim = [-.2 .3];
- % ft_multiplotER(cfg, natural, manmade)
- %% plot topography of visual N1
- % set custom colormap
- cBrew = brewermap(500,'RdBu');
- cBrew = flipud(cBrew);
- colormap(cBrew)
- h = figure('units','centimeters','position',[0 0 10 10]);
- set(gcf,'renderer','Painters')
- cfg = [];
- cfg.layout = 'EEG1005.lay';
- cfg.parameter = 'powspctrm';
- cfg.comment = 'no';
- cfg.colormap = cBrew;
- cfg.colorbar = 'EastOutside';
- plotData = [];
- plotData.label = elec.label(1:64); % {1 x N}
- plotData.dimord = 'chan';
- plotData.powspctrm = squeeze(nanmean(nanmean(nanmin(mergeddata(:,:,time>0.04 & time <0.12,1:2),[],3),1),4))';
- [~, sortidx] = sort(plotData.powspctrm, 'ascend');
- idx_chans = sortidx(1);
- cfg.marker = 'off';
- cfg.highlight = 'yes';
- cfg.highlightchannel = plotData.label(idx_chans);
- cfg.highlightcolor = [1 0 0];
- cfg.highlightsymbol = '.';
- cfg.highlightsize = 18;
- cfg.zlim = [-10 10]*10^-4;
- cfg.figure = h;
- ft_topoplotER(cfg,plotData);
- cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Amplitude');
- set(findall(gcf,'-property','FontSize'),'FontSize',12)
- figureName = ['b_topo_minimum'];
- saveas(h, fullfile(pn.figures, figureName), 'epsc');
- saveas(h, fullfile(pn.figures, figureName), 'png');
- %% visualize N1 over negative maximum
- % avg across channels and conditions
- condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
- condAvg1 = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1),2),4));
- condAvg2 = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,2),2),4));
- h = figure('units','centimeters','position',[0 0 10 8]);
- cla; hold on;
- % new value = old value ? subject average + grand average
- curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
- curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
- standError = nanstd(curData,1)./sqrt(size(curData,1));
- l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
- curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
- curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
- standError = nanstd(curData,1)./sqrt(size(curData,1));
- l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
- % ax = gca; ax.YDir = 'reverse';
- xlabel('Time (s) from stim onset')
- xlim([-.025 .16]); %ylim([-.03 .18])
- ylabel({'ERP';'(microVolts)'});
- xlabel({'Time (s)'});
- set(findall(gcf,'-property','FontSize'),'FontSize',12)
- %% ERP components for subjects 1-17 and 18-33 are shifted in time!
- h = figure('units','centimeters','position',[0 0 10 8]);
- cla; hold on;
- curData = squeeze(nanmean(mergeddata(1:17,idx_chans,:,1),2));
- standError = nanstd(curData,1)./sqrt(size(curData,1));
- l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
- curData = squeeze(nanmean(mergeddata(18:end,idx_chans,:,1),2));
- standError = nanstd(curData,1)./sqrt(size(curData,1));
- l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
- % ax = gca; ax.YDir = 'reverse';
- legend({'Initial 17', 'Final 16'}, 'location', 'NorthWest'); legend('boxoff')
- xlabel('Time (s) from stim onset')
- xlim([-.025 .16]); %ylim([-.03 .18])
- ylabel({'ERP';'(microVolts)'});
- xlabel({'Time (s)'});
- set(findall(gcf,'-property','FontSize'),'FontSize',12)
- figureName = ['b_twogroups'];
- saveas(h, fullfile(pn.figures, figureName), 'epsc');
- saveas(h, fullfile(pn.figures, figureName), 'png');
- %% align individual subjects N1 to first negative peak
- % find individual minimum (avg. across conditions) between 40 and 120 ms
- time2search = find(time>0.04 & time <0.12);
- newtime = 0-100*(time(2)-time(1)):(time(2)-time(1)):0+100*(time(2)-time(1));
- for indID = 1:size(condAvg,1)
- tmp = find(islocalmin(condAvg(indID,time2search), ...
- 'FlatSelection', 'center', ...
- 'MinSeparation', 25));
- minVal1(indID) = condAvg1(indID,time2search(tmp(1)));
- curmin = time2search(tmp(1));
- alignedN1_1(indID,:) = condAvg1(indID,curmin-100:curmin+100);
- minVal2(indID) = condAvg2(indID,time2search(tmp(1)));
- curmin = time2search(tmp(1));
- alignedN1_2(indID,:) = condAvg2(indID,curmin-100:curmin+100);
-
- alignedTopo(indID,:,:) = squeeze(nanmean(mergeddata(indID,:,curmin-100:curmin+100,1:2),4));
- end
- % alternatively: consider global minimum
- % for indID = 1:size(condAvg,1)
- % [~, minLoc1(indID)] = min(condAvg(indID,time2search));
- % curmin = time2search(minLoc1(indID));
- % minVal1(indID) = condAvg1(indID,curmin);
- % alignedN1_1(indID,:) = condAvg1(indID,curmin-100:curmin+100);
- % [~, minLoc2(indID)] = min(condAvg(indID,time2search));
- % curmin = time2search(minLoc2(indID));
- % minVal1(indID) = condAvg2(indID,curmin);
- % alignedN1_2(indID,:) = condAvg2(indID,curmin-100:curmin+100);
- % end
- mergeddata_aligned = cat(3, alignedN1_1, alignedN1_2);
- [h, p, ci, stats] = ttest(minVal1, minVal2)
- % avg across channels and conditions
- condAvg_al = squeeze(nanmean(mergeddata_aligned(:,:,1:2),3));
- % check that troughs are aligned
- % figure; imagesc(zscore(condAvg_al,[],2))
- % figure; imagesc(condAvg_al)
- h = figure('units','centimeters','position',[0 0 10 8]);
- cla; hold on;
- % new value = old value ? subject average + grand average
- curData = squeeze(mergeddata_aligned(:,:,1));
- curData = curData-condAvg_al+repmat(nanmean(condAvg_al,1),size(condAvg_al,1),1);
- standError = nanstd(curData,1)./sqrt(size(curData,1));
- l1 = shadedErrorBar(newtime*1000,nanmean(curData,1),standError, 'lineprops', {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
- curData = squeeze(mergeddata_aligned(:,:,2));
- curData = curData-condAvg_al+repmat(nanmean(condAvg_al,1),size(condAvg_al,1),1);
- standError = nanstd(curData,1)./sqrt(size(curData,1));
- l2 = shadedErrorBar(newtime*1000,nanmean(curData,1),standError, 'lineprops', {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
- % ax = gca; ax.YDir = 'reverse';
- legend({'manmade', 'natural'}, 'location', 'NorthWest'); legend('boxoff')
- xlabel('Time (s) from local minimum')
- xlim([-100 100]); %ylim([-.03 .18])
- ylabel({'ERP';'(microVolts)'});
- xlabel({'Time (ms) from local minimum'});
- set(findall(gcf,'-property','FontSize'),'FontSize',12)
- figureName = ['b_scenecat_N1'];
- saveas(h, fullfile(pn.figures, figureName), 'epsc');
- saveas(h, fullfile(pn.figures, figureName), 'png');
- %% plot topography around detected trough
- h = figure('units','centimeters','position',[0 0 10 10]);
- set(gcf,'renderer','Painters')
- cfg = [];
- cfg.layout = 'EEG1010.lay';
- cfg.parameter = 'powspctrm';
- cfg.comment = 'no';
- cfg.colormap = cBrew;
- cfg.colorbar = 'EastOutside';
- plotData = [];
- plotData.label = elec.label(1:64); % {1 x N}
- plotData.dimord = 'chan';
- plotData.powspctrm = squeeze(nanmean(nanmean(alignedTopo(:, :, newtime>-.01 & newtime < .01),3),1))';
- [~, sortidx] = sort(plotData.powspctrm, 'ascend');
- idx_chans = sortidx(1);
- idx_chans_visual = idx_chans;
- cfg.marker = 'off';
- cfg.highlight = 'yes';
- cfg.highlightchannel = plotData.label(idx_chans);
- cfg.highlightcolor = [1 0 0];
- cfg.highlightsymbol = '.';
- cfg.highlightsize = 18;
- cfg.zlim = [-5 5]*10^-4;
- cfg.figure = h;
- ft_topoplotER(cfg,plotData);
- cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Amplitude');
- set(findall(gcf,'-property','FontSize'),'FontSize',12)
- figureName = ['b_trough_topo'];
- saveas(h, fullfile(pn.figures, figureName), 'epsc');
- saveas(h, fullfile(pn.figures, figureName), 'png');
- %% plot raincloudplot of extracted N1 peak amplitudes
- % Note that the stats reported here may vary from above in the case of
- % outliers.
- conds = {'manmade'; 'natural'};
- uData{1} = [minVal1; minVal2];
- cBrew(1,:) = [.6 .6 .6];
- idx_outlier = cell(1); idx_standard = cell(1);
- for indGroup = 1
- dataToPlot = uData{indGroup}';
- % define outlier as lin. modulation that is more than three scaled median absolute deviations (MAD) away from the median
- X = [1 1; 1 2]; b=X\dataToPlot'; IndividualSlopes = b(2,:);
- IndividualSlopes = uData{indGroup}(2,:)-uData{indGroup}(1,:);
- outliers = isoutlier(IndividualSlopes, 'median');
- idx_outlier{indGroup} = find(outliers);
- idx_standard{indGroup} = find(outliers==0);
- end
- h = figure('units','centimeter','position',[0 0 10 8]);
- for indGroup = 1
- dataToPlot = uData{indGroup}';
- % read into cell array of the appropriate dimensions
- data = []; data_ws = [];
- for i = 1:2
- for j = 1:1
- data{i, j} = dataToPlot(:,i);
- % individually demean for within-subject visualization
- data_ws{i, j} = dataToPlot(:,i)-...
- nanmean(dataToPlot(:,:),2)+...
- repmat(nanmean(nanmean(dataToPlot(:,:),2),1),size(dataToPlot(:,:),1),1);
- data_nooutlier{i, j} = data{i, j};
- data_nooutlier{i, j}(idx_outlier{indGroup}) = [];
- data_ws_nooutlier{i, j} = data_ws{i, j};
- data_ws_nooutlier{i, j}(idx_outlier{indGroup}) = [];
- % sort outliers to back in original data for improved plot overlap
- data_ws{i, j} = [data_ws{i, j}(idx_standard{indGroup}); data_ws{i, j}(idx_outlier{indGroup})];
- end
- end
- % IMPORTANT: plot individually centered estimates, stats on uncentered estimates!
- set(gcf,'renderer','Painters')
- cla;
- cl = cBrew(indGroup,:);
- [~, dist] = rm_raincloud_fixedSpacing(data_ws, [.8 .8 .8],1);
- h_rc = rm_raincloud_fixedSpacing(data_ws_nooutlier, cl,1,[],[],[],dist);
- view([90 -90]);
- axis ij
- box(gca,'off')
- set(gca, 'YTickLabels', {'natural'; 'manmade'});
- yticks = get(gca, 'ytick'); ylim([yticks(1)-(yticks(2)-yticks(1))./2, yticks(2)+(yticks(2)-yticks(1))./1.5]);
- minmax = [min(min(cat(2,data_ws{:}))), max(max(cat(2,data_ws{:})))];
- xlim(minmax+[-0.2*diff(minmax), 0.2*diff(minmax)])
- ylabel('scene category'); xlabel({'ERP'; '[Individually centered]'})
- % test linear effect
- curData = [data_nooutlier{1, 1}, data_nooutlier{2, 1}];
- X = [1 1; 1 2]; b=X\curData'; IndividualSlopes = b(2,:);
- [~, p, ci, stats] = ttest(IndividualSlopes);
- title(['M:', num2str(round(mean(IndividualSlopes),6)), '; p=', num2str(round(p,6))])
- end
- set(findall(gcf,'-property','FontSize'),'FontSize',12)
- figureName = ['b_rcp'];
- saveas(h, fullfile(pn.figures, figureName), 'epsc');
- saveas(h, fullfile(pn.figures, figureName), 'png');
- %% plot topography for different timewindows
- %
- % timewins = [0.03 0.075; ...
- % 0.075, 0.1; ...
- % 0.1, 0.15; ...
- % 0.495, 0.535];
- %
- % % set custom colormap
- % cBrew = brewermap(500,'RdBu');
- % cBrew = flipud(cBrew);
- % colormap(cBrew)
- %
- % h = figure('units','centimeters','position',[0 0 10 10]);
- % set(gcf,'renderer','Painters')
- %
- % for indTime = 1:size(timewins,1)
- % subplot(3,2,indTime)
- %
- % cfg = [];
- % cfg.layout = 'EEG1010.lay';
- % cfg.parameter = 'powspctrm';
- % cfg.comment = 'no';
- % cfg.colormap = cBrew;
- % cfg.colorbar = 'EastOutside';
- %
- % plotData = [];
- % plotData.label = elec.label(1:64); % {1 x N}
- % plotData.dimord = 'chan';
- % plotData.powspctrm = squeeze(nanmean(nanmean(nanmean(mergeddata(:,:,time>=timewins(indTime,1) & time <=timewins(indTime,2),1),3),1),4))';
- % [~, sortidx] = sort(plotData.powspctrm, 'ascend');
- % idx_chans = sortidx(1:6);
- %
- % cfg.marker = 'off';
- % cfg.highlight = 'yes';
- % cfg.highlightchannel = plotData.label(idx_chans);
- % cfg.highlightcolor = [1 0 0];
- % cfg.highlightsymbol = '.';
- % cfg.highlightsize = 18;
- % cfg.zlim = [-10 10]*10^-4;
- % cfg.figure = h;
- % ft_topoplotER(cfg,plotData);
- % cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Amplitude');
- %
- % end
|