{"id":352,"date":"2019-04-02T17:01:22","date_gmt":"2019-04-02T15:01:22","guid":{"rendered":"http:\/\/benediktehinger.de\/blog\/science\/?p=352"},"modified":"2019-04-09T08:11:43","modified_gmt":"2019-04-09T06:11:43","slug":"simulating-data-with-sereega-multiple-comparison-corrections","status":"publish","type":"post","link":"https:\/\/benediktehinger.de\/blog\/science\/simulating-data-with-sereega-multiple-comparison-corrections\/","title":{"rendered":"Simulating EEG\/ERP Data with SEREEGA & multiple comparison corrections"},"content":{"rendered":"\n

For the $i^2 c^2 s^3$ summer school I simulated quite a bit of data and analyzed them with several common multiple comparison methods. I used the SEREEGA<\/a> toolbox for the simulation. All the MatLab-code can be found at the end of this post. In a follow-up blogpost I will extend the toolbox to continuous data that we can analyze with the unfold<\/a> toolbox.<\/p>\n\n\n\n

First, I simulated data based on three effects: Two early dipoles representing the P100, one right lateralized for the N170 and a deep one for the P300.<\/p>\n\n\n\n

\"\"
Two symmetric sources for the P1 in the occipital cortex, one lateralized to the right cortex (third panel right most) for the N170 and a deep one, central for the P300.<\/figcaption><\/figure>\n\n\n\n

I added brown and white noise to the simulated epochs. An exemplary eegplot shows that it kind of looks like EEG data.<\/p>\n\n\n\n

\"\"
5 trials in the default pop_eegplot of eeglab<\/figcaption><\/figure>\n\n\n\n

I generate two conditions with two different condition-differences, one on the N170 and one on the P300 . The same data are depicted in the following three plots:
<\/p>\n\n\n\n

\"\"
Data in butterfly plot & time-aligned points, plotted using the eegvis<\/a> toolbox<\/figcaption><\/figure><\/div>\n\n\n\n

In the first row we see a butterfly plot of activation in all channels (the colored channels are depicted in the next plot). In the first topoplot-row we have the average activation and in the last row we have the condition differences.<\/p>\n\n\n\n

\"\"<\/figure><\/div>\n\n\n\n
\"\"
The three types of components: a short but relatively strong (P100, red), a similar effect, short but strong (N170, condition difference) and a long effect (P300, with weak condition effect). Highlighted gray areas are approximate regions where the components showed condition differences. The noise was a mix between brown and white noise. <\/figcaption><\/figure>\n\n\n\n

The red line is the P100, occipital effect. No difference between conditions. The green one is the N170, temporal effect. Only visible in one condition. The blue one is a P300 like, deep effect. Difference in amplitude between conditions.<\/p>\n\n\n\n

\"\"
A third way to look at the same data. Color indicates activation (blue => negative \u00b5V, yellow=> positive \u00b5V). Each row is one channel, each column one time point. Using such a plot, one looses spatial relations between the channels , but gains an nice overview of all the data. Note that the colorscales are different between the two panels.<\/figcaption><\/figure><\/div>\n\n\n\n

Now we are ready to do some statistical testing and multiple comparison correction. I choose to compare uncorrected, bonferroni-holms, FDR <\/a>(Benjamini-Hochberg), clustermass cluster permutation testing (LIMO toolbox) <\/a>and threshold free cluster enhancement (EPT-TFCE toolbox).<\/a><\/p>\n\n\n\n

<\/p>\n\n\n\n

\"\"
I addapted a display from fMRI depicting t-value vs. significance based on the fmri slice display toolbox<\/a>. Plot can be made using the eegvis toolbox<\/a><\/figcaption><\/figure><\/div>\n\n\n\n

In this instance, uncorrected p-value did not so bad, Bonferroni-Holms is, as expected, quite conservative. FDR hat troubles with the elongated cluster and TFCE\/Cluster-permutation with the short one. <\/p>\n\n\n\n

Note: This is only a single simulation. In order to move from these anectotal findings to proper statements, one would need to repeat the simulation 1000 times and see how often which samples were deemed significant (see Groppe 2011).<\/p>\n\n\n\n

Thanks to Anna Lisa Gert<\/a> for help with writing this post.<\/p>\n\n\n\n

This is the code I used, it is a bit messy because it was programmed in a tutorial-style and I kept my debug-statements in there.<\/p>\n\n\n\n


<\/p>\n\n\n\n

rng(1)\n%% Adding toolboxes\naddpath('lib\\eeglab\\')\nif ~exist('pop_loadset')\n    eeglab\nend\naddpath(genpath('lib\\sereega'))\naddpath(genpath('lib\\eegvis'))\n% Threshold Free cluster enhancement\naddpath(genpath('lib\\unfold\\lib\\ept_TFCE'))\naddpath(genpath('lib\\limo_eeg'))\n%% SEREEGA Data Generation\ncfg = struct('debug',0)\n%% load headmodel \nlf = lf_generate_fromnyhead('montage', 'S64');\n\nif cfg.debug\n    plot_headmodel(lf);\nend\n\n%% set seed\n\n%% noise\nnoise_brown = struct( ...\n    'type', 'noise', ...\n    'color', 'brown', ...\n    'amplitude', 2);\nnoise_brown = utl_check_class(noise_brown);\n\nnoise_white = struct( ...\n    'type', 'noise', ...\n    'color', 'white', ...\n    'amplitude', 2);\nnoise_white = utl_check_class(noise_white);\n\n\n%% visual Cortex N170\nvis = [];\n\nvis.signal{1} = struct();\nvis.signal{1}.peakLatency = 170;      % in ms, starting at the start of the epoch\nvis.signal{1}.peakWidth = 100;        % in ms\nvis.signal{1}.peakAmplitude = -5;      % in microvolt\nvis.signal{1} = utl_check_class(vis.signal{1}, 'type', 'erp');\nvis.signal{2} = noise_brown;\n\n\n% right visual cortex\nvis.source= lf_get_source_nearest(lf, [50 -40 -25]);\nif cfg.debug\n    plot_source_location(vis.source, lf, 'mode', '3d','shrink',0);\nend\n\n% vis.orientation = utl_get_orientation_pseudotangential(vis.source,lf); \nvis.orientationDv = [0 0 0]; % orientation does not change between trials\nvis.orientation = [0.5,-0.46,1]; % I looked at the maps and found that this prientation most closely looked like a N170\n\nif cfg.debug\n    %% This Code I used to generate random orientations and inspect the leadfield projections\n    for k = 1:10\n        vis.orientation = utl_get_orientation_random(1)\n        plot_source_projection(vis.source, lf, 'orientation', vis.orientation,'orientedonly',1);\n        %     title(vis.orientation)\n    end\nend\n%% Generate P100, take from SEREEGA Toolbox\np100 = utl_get_component_fromtemplate('visual_p100_erp',lf);\nif cfg.debug\n    % Visualize the response \n    plot_signal_fromclass(vis.signal{1}, epochs);\nend\n\n%% Generate P300\np3 = [];\n\np3.signal{1} = struct();\np3.signal{1}.peakLatency = 350;      % in ms, starting at the start of the epoch\np3.signal{1}.peakWidth = 400;        % in ms\np3.signal{1}.peakAmplitude = -3;      % in microvolt\np3.signal{1} = utl_check_class(p3.signal{1}, 'type', 'erp');\np3.signal{2} = noise_brown;\n\n\n% somewhere deep in the brain\np3.source= lf_get_source_nearest(lf, [0 -40 -25]);\nif cfg.debug\n    %% visualize it using SEREEGA\n    plot_source_location(p3.source, lf, 'mode', '3d');\nend\n\np3.orientation = utl_get_orientation_pseudotangential(p3.source,lf);\np3.orientationDv = [0 0 0];\n\np3.orientation = [0.03,-0.16,1];\nif cfg.debug\n    %% again find an orientation that somehow matches P300topography\n    for k = 1:10\n        p3.orientation = utl_get_orientation_random(1)\n        plot_source_projection(p3.source, lf, 'orientation', p3.orientation,'orientedonly',1);\n    end\nend\n\n\n%% generate 25 random noise sources and combine all components\nsources = lf_get_source_spaced(lf, 10, 25);\nnoise_brown.amplitude = 5;\ncomps = utl_create_component(sources, noise_brown, lf);\ncomps = utl_add_signal_tocomponent(noise_white,comps);\n[comps1, comps2] = deal(comps);\n\n% add P100 to noise components\ncomps1(end+1:end+2) = p100; %bilateral\ncomps2(end+1:end+2) = p100;\n\n% N170 with condition difference\nvis.signal{1}.peakAmplitude = -3;\ncomps1(end+1) = vis;\nvis.signal{1}.peakAmplitude = -1;\ncomps2(end+1) = vis;\n\n% broad, weak P300, with condition difference\np3.signal{1}.peakAmplitude = -10;\ncomps1(end+1) = p3;\np3.signal{1}.peakAmplitude = -12;\ncomps2(end+1) = p3;\n%% Generate The EEG Data\nepochs = struct();\nepochs.n = 50;             % the number of epochs to simulate\nepochs.srate = 250;        % their sampling rate in Hz\nepochs.length = 500;\n\ndata1 = generate_scalpdata(comps1, lf, epochs);\ndata2 = generate_scalpdata(comps2, lf, epochs);\nEEG1 = utl_create_eeglabdataset(data1, epochs, lf, 'marker', 'event1');\nEEG2 = utl_create_eeglabdataset(data2, epochs, lf, 'marker', 'event2');\nEEG = utl_reorder_eeglabdataset(pop_mergeset(EEG1, EEG2));\n\n%% First inspection, this should match the intended topographies\npop_topoplot(EEG, 1, [0 100 170,300], '', [1 8]);\n%% ERP-Like Figure\nfigure,\nfindtype = @(type)cellfun(@(x)strcmp(type,x),{EEG.event.type});\n\nploterp = @(chan)plot(EEG.times,...\n    [mean(squeeze(EEG.data(chan,:,findtype('event1'))),2),...\n    mean(squeeze(EEG.data(chan,:,findtype('event2'))),2)]);\nfigure,ploterp(44),title(EEG.chanlocs(44).labels)\nhold on,ploterp(63),title(EEG.chanlocs(63).labels)\nhold on,ploterp(30),title(EEG.chanlocs(30).labels)\n\n%% Within-Subject T-Test using LIMO\n\n[m,dfe,ci,sd,n,t,p] = limo_ttest(2,EEG.data(:,:,findtype('event1')),EEG.data(:,:,findtype('event2')),0.05);\n%% simple imagesc plots\nfigure,subplot(2,1,1),imagesc(EEG.times,1:64,mean(EEG.data,3));box off\nsubplot(2,1,2),imagesc(EEG.times,1:64,m),box off\n% Show data\nfigure\nimagesc(m)\n\nfigure\nimagesc(t)\n%% Prepare statistics\n\n\n%%\nept_tfce_nb = ept_ChN2(EEG.chanlocs,0); % the 1 to plot\ntfce_res = ept_TFCE(...\n    permute(EEG.data(:,:,findtype('event1')),[3 1 2]),...\n    permute(EEG.data(:,:,findtype('event2')),[3,1,2]),EEG.chanlocs,'nperm',800,'flag_save',0,'ChN',ept_tfce_nb)\n\n%%\nt_h0 = nan([size(EEG.data,1),size(EEG.data,2),800]);\nfor b = 1:800\n    fprintf('%i\/%i\\n',b,800)\n    r_perm      = randperm(epochs.n*2); % Consider using Shuffle mex here (50-85% faster)...\n    \n    nData       = EEG.data(:,:,r_perm);\n    sData{1}    = nData(:,:,1:epochs.n);\n    sData{2}    = nData(:,:,(epochs.n+1):(epochs.n*2));\n    \n    t_h0(:,:,b) = (mean(sData{1},3)-mean(sData{2},3)).\/sqrt(var(sData{1},[],3)\/epochs.n+var(sData{2},[],3)\/epochs.n);\n    \n    \nend\np_h0 = tpdf(t_h0,199);\n%% Generate Cluster Permutation\nlimo_nb = limo_neighbourdist(EEG,50);\n[limomask,cluster_p,max_th] = limo_clustering(t.^2,p,t_h0.^2,p_h0,struct('data',struct('chanlocs',EEG.chanlocs','neighbouring_matrix',limo_nb)),2,0.05,0);\ncluster_p(isnan(cluster_p(:))) = 1;\n% F_tfce = limo_tfce(2,t.^2,limo_nb);\n\n%% show significance methods\nfigure,\nfor k = 0:4\n    switch k\n        case 0\n            plot_p = p;\n            mask = plot_p<0.05;\n        case 1\n            plot_p = bonf_holm(double(p),0.05);\n            mask = plot_p<0.05;\n        case 2\n            [~,~,~,plot_p] = fdr_bh(p,0.05,'pdep','no'); % see groppe 2011 why pdep is fine\n            mask = plot_p<0.05;\n            \n        case 3\n            plot_p = cluster_p;\n            mask = plot_p<0.05;\n            \n        case 4\n            plot_p = tfce_res.P_Values;\n            mask = plot_p<0.05;\n            \n            \n    end\n    plot_data= plot_p;\n    plot_data = log10(plot_data);\n    subplot(5,1,(k)+1)\n    eegvis_imagesc(m,t,'chanlocs',EEG.chanlocs,'times',EEG.times,...\n        'contour',1,'clustermask',mask,'figure',0,'xlabel',0,'colorbar',k == 4)\n%     caxis([-7,0])\n%     if k ~=4\n%         axis off\n%     end\n%     box off\n    plot_data(~mask) = 1;\n    box off\nend\n\n%%\nhA = plot_topobutter(mean(EEG.data,3),EEG.times,EEG.chanlocs,'colormap',{'div'},'quality',40,'n_topos',12);\n\n%%\nfigure\nhA = plot_topobutter(cat(3,mean(EEG.data,3),m),EEG.times,EEG.chanlocs,'quality',32,'individualcolorscale','row','highlighted_channel',[44,63]);\n\n%% Additional Blog Figures\nplot_source_location([p100.source,vis.source,p3.source], lf, 'mode', '2d');\n\n\n\n\n<\/pre>\n","protected":false},"excerpt":{"rendered":"
For the $i^2 c^2 s^3$ summer school I simulated quite a bit of data and analyzed them with several common multiple comparison methods. I used the SEREEGA toolbox for the simulation. All the MatLab-code can be found at the end of this post. In a follow-up blogpost I will extend the toolbox to continuous data that we can analyze with the unfold toolbox. First, I simulated data based on three effects: Two early dipoles representing the P100, one right lateralized for the N170 and a deep one for the P300. I added brown and white noise to the simulated epochs….<\/p>\n","protected":false},"author":2,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[5],"tags":[],"class_list":["post-352","post","type-post","status-publish","format-standard","hentry","category-blog"],"_links":{"self":[{"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/posts\/352","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/comments?post=352"}],"version-history":[{"count":0,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/posts\/352\/revisions"}],"wp:attachment":[{"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/media?parent=352"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/categories?post=352"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/tags?post=352"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}