b1_taskPLS_scenecat_N1_plotLV1.m 6.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181
  1. % Create an overview plot featuring the results of the multivariate PLS
  2. % comparing spectral changes during the stimulus period under load
  3. clear all; cla; clc;
  4. currentFile = mfilename('fullpath');
  5. [pathstr,~,~] = fileparts(currentFile);
  6. cd(fullfile(pathstr,'..'))
  7. rootpath = pwd;
  8. pn.data = fullfile(rootpath, 'data', 'stats');
  9. pn.figures = fullfile(rootpath, 'figures');
  10. pn.tools = fullfile(rootpath, 'tools');
  11. addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
  12. addpath(fullfile(pn.tools, 'fieldtrip')); ft_defaults;
  13. addpath(genpath(fullfile(pn.tools, 'RainCloudPlots')));
  14. addpath(fullfile(pn.tools, 'BrewerMap'));
  15. addpath(fullfile(pn.tools, 'winsor'));
  16. % set custom colormap
  17. cBrew = brewermap(500,'RdBu');
  18. cBrew = flipud(cBrew);
  19. colormap(cBrew)
  20. load(fullfile(pn.data, 'b01_taskpls_erp.mat'),...
  21. 'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')
  22. load(fullfile(rootpath, 'data','erp', ['sub-001_erp.mat']));
  23. elec = erp.scene_category{1}.elec;
  24. result.perm_result.sprob
  25. indLV = 1;
  26. lvdat = reshape(result.boot_result.compare_u(:,indLV), num_chans, num_freqs, num_time);
  27. stat.prob = lvdat;
  28. stat.mask = lvdat > 3 | lvdat < -3;
  29. % maskNaN = double(stat.mask);
  30. % maskNaN(maskNaN==0) = NaN;
  31. %% invert solution
  32. % stat.mask = stat.mask;
  33. % stat.prob = stat.prob.*-1;
  34. % result.vsc = result.vsc.*-1;
  35. % result.usc = result.usc.*-1;
  36. h = figure('units','centimeter','position',[0 0 15 10]);
  37. set(gcf,'renderer','Painters')
  38. statsPlot = [];
  39. statsPlot = cat(1, statsPlot,squeeze(nanmax(abs(stat.prob(1:64,:,:)),[],1)));
  40. plot(stat.time,statsPlot, 'k')
  41. xlabel('Time [s from stim onset]'); ylabel('max abs BSR');
  42. title({'ERP changes'; ['p = ', num2str(round(result.perm_result.sprob(indLV),4))]})
  43. set(findall(gcf,'-property','FontSize'),'FontSize',18)
  44. figureName = ['b01_pls_traces'];
  45. saveas(h, fullfile(pn.figures, figureName), 'epsc');
  46. saveas(h, fullfile(pn.figures, figureName), 'png');
  47. h = figure('units','centimeter','position',[0 0 15 10]);
  48. set(gcf,'renderer','Painters')
  49. statsPlot = [];
  50. statsPlot = cat(1, statsPlot,squeeze(mean(stat.prob(1:64,:,:),1)));
  51. plot(stat.time,statsPlot, 'k')
  52. xlabel('Time [s from stim onset]'); ylabel('mean BSR');
  53. title({'ERP changes'; ['p = ', num2str(round(result.perm_result.sprob(indLV),4))]})
  54. set(findall(gcf,'-property','FontSize'),'FontSize',18)
  55. xlim([0 0.3])
  56. figureName = ['b01_pls_traces'];
  57. saveas(h, fullfile(pn.figures, figureName), 'epsc');
  58. saveas(h, fullfile(pn.figures, figureName), 'png');
  59. %% plot multivariate topographies
  60. h = figure('units','centimeters','position',[0 0 10 10]);
  61. set(gcf,'renderer','Painters')
  62. cfg = [];
  63. cfg.layout = 'biosemi64.lay';
  64. cfg.parameter = 'powspctrm';
  65. cfg.comment = 'no';
  66. cfg.colormap = cBrew;
  67. cfg.colorbar = 'EastOutside';
  68. plotData = [];
  69. plotData.label = elec.label; % {1 x N}
  70. plotData.dimord = 'chan';
  71. cfg.zlim = [-6 6];
  72. cfg.figure = h;
  73. plotData.powspctrm = squeeze(nanmean(stat.mask(:,:,stat.time >0.08 & stat.time <0.12).*...
  74. stat.prob(:,:,stat.time >0.08 & stat.time <0.12),3));
  75. ft_topoplotER(cfg,plotData);
  76. cb = colorbar('location', 'EastOutside'); set(get(cb,'ylabel'),'string','Mean BSR');
  77. figureName = ['b01_lv1'];
  78. % saveas(h, fullfile(pn.figures, figureName), 'epsc');
  79. % saveas(h, fullfile(pn.figures, figureName), 'png');
  80. %% plot using raincloud plot
  81. groupsizes=result.num_subj_lst;
  82. conditions=lv_evt_list;
  83. conds = {'manmade'; 'natural'};
  84. condData = []; uData = [];
  85. for indGroup = 1
  86. if indGroup == 1
  87. relevantEntries = 1:groupsizes(1)*numel(conds);
  88. elseif indGroup == 2
  89. relevantEntries = groupsizes(1)*numel(conds)+1:...
  90. groupsizes(1)*numel(conds)+groupsizes(2)*numel(conds);
  91. end
  92. for indCond = 1:numel(conds)
  93. targetEntries = relevantEntries(conditions(relevantEntries)==indCond);
  94. condData{indGroup}(indCond,:) = result.vsc(targetEntries,indLV);
  95. uData{indGroup}(indCond,:) = result.usc(targetEntries,indLV);
  96. end
  97. end
  98. %% plot RainCloudPlot (within-subject centered)
  99. cBrew(1,:) = 2.*[.3 .1 .1];
  100. cBrew(2,:) = [.6 .6 .6];
  101. idx_outlier = cell(1); idx_standard = cell(1);
  102. for indGroup = 1
  103. dataToPlot = uData{indGroup}';
  104. % define outlier as lin. modulation that is more than three scaled median absolute deviations (MAD) away from the median
  105. X = [1 1; 1 2]; b=X\dataToPlot'; IndividualSlopes = b(2,:);
  106. outliers = isoutlier(IndividualSlopes, 'median');
  107. idx_outlier{indGroup} = find(outliers);
  108. idx_standard{indGroup} = find(outliers==0);
  109. end
  110. h = figure('units','centimeter','position',[0 0 25 10]);
  111. for indGroup = 1
  112. dataToPlot = uData{indGroup}';
  113. % read into cell array of the appropriate dimensions
  114. data = []; data_ws = [];
  115. for i = 1:2
  116. for j = 1:1
  117. data{i, j} = dataToPlot(:,i);
  118. % individually demean for within-subject visualization
  119. data_ws{i, j} = dataToPlot(:,i)-...
  120. nanmean(dataToPlot(:,:),2)+...
  121. repmat(nanmean(nanmean(dataToPlot(:,:),2),1),size(dataToPlot(:,:),1),1);
  122. data_nooutlier{i, j} = data{i, j};
  123. data_nooutlier{i, j}(idx_outlier{indGroup}) = [];
  124. data_ws_nooutlier{i, j} = data_ws{i, j};
  125. data_ws_nooutlier{i, j}(idx_outlier{indGroup}) = [];
  126. % sort outliers to back in original data for improved plot overlap
  127. data_ws{i, j} = [data_ws{i, j}(idx_standard{indGroup}); data_ws{i, j}(idx_outlier{indGroup})];
  128. end
  129. end
  130. % IMPORTANT: plot individually centered estimates, stats on uncentered estimates!
  131. subplot(1,2,indGroup);
  132. set(gcf,'renderer','Painters')
  133. cla;
  134. cl = cBrew(indGroup,:);
  135. [~, dist] = rm_raincloud_fixedSpacing(data_ws, [.8 .8 .8],1);
  136. h_rc = rm_raincloud_fixedSpacing(data_ws_nooutlier, cl,1,[],[],[],dist);
  137. view([90 -90]);
  138. axis ij
  139. box(gca,'off')
  140. set(gca, 'YTickLabels', {'natural'; 'manmade'});
  141. yticks = get(gca, 'ytick'); ylim([yticks(1)-(yticks(2)-yticks(1))./2, yticks(2)+(yticks(2)-yticks(1))./1.5]);
  142. minmax = [min(min(cat(2,data_ws{:}))), max(max(cat(2,data_ws{:})))];
  143. xlim(minmax+[-0.2*diff(minmax), 0.2*diff(minmax)])
  144. ylabel('scene category'); xlabel({'Brainscore'; '[Individually centered]'})
  145. % test linear effect
  146. curData = [data_nooutlier{1, 1}, data_nooutlier{2, 1}];
  147. X = [1 1; 1 2]; b=X\curData'; IndividualSlopes = b(2,:);
  148. [~, p, ci, stats] = ttest(IndividualSlopes);
  149. title(['M:', num2str(round(mean(IndividualSlopes),3)), '; p=', num2str(round(p,3))])
  150. end
  151. figureName = ['b01_rcp'];
  152. % saveas(h, fullfile(pn.figures, figureName), 'epsc');
  153. % saveas(h, fullfile(pn.figures, figureName), 'png');