MEEG time-frequency searchlight

This example shows MVPA analyses performed on MEEG data, using a searchlight across the time, frequency and channel dimensions

The input dataset involved a paradigm with electrical median nerve stimulation for durations of 2s at 20Hz.

The code presented here can be adapted for other MEEG analyses, but there are a few potential caveats: * assignment of targets (labels of conditions) is based here on stimulation periods versus pre-stimulation periods. In typical analyses the targets should be based on different trial conditions, for example as set a FieldTrip .trialinfo field. * assignment of chunks (parts of the data that are assumed to be independent) is based on a trial-by-trial basis. For cross-validation, the number of chunks is reduced to two to speed up the analysis. * the time window used for analyses is rather small. This means that in particular for time-freq analysis a lot of data is missing, especially for early and late timepoints in the lower frequency bands. For typical analyses it may be preferred to use a wider time window. * the current examples do not perform baseline corrections or signal normalizations, which may reduce discriminatory power.

Note: running this code requires FieldTrip.

  1. For CoSMoMVPA's copyright information and license terms, #
  2. see the COPYING file distributed with CoSMoMVPA. #

Contents

get timelock data in CoSMoMVPA format

% set configuration
config=cosmo_config();
data_path=fullfile(config.tutorial_data_path,'meg_20hz');

% show dataset information
readme_fn=fullfile(data_path,'README');
cosmo_type(readme_fn);

% reset citation list
cosmo_check_external('-tic');

% load data
data_fn=fullfile(data_path,'subj102_B01_20Hz_timefreq.mat');
data_tf=load(data_fn);

% convert to cosmomvpa struct
ds=cosmo_meeg_dataset(data_tf);

% set the target (trial condition)
% here, 1=pre-stimulus, 2=peri-stimulus
ds.sa.targets=ds.sa.trialinfo(:,1);

% set the chunks (independent measurements)
% in this dataset, the first half of the samples (in order)
% are the post-trials;
% the second half the pre-trials
ds.sa.chunks=[(1:145) (1:145)]';

% just to check everything is ok
cosmo_check_dataset(ds);

% get rid of features with at least one NaN value across samples
fa_nan_mask=sum(isnan(ds.samples),1)>0;
fprintf('%d / %d features have NaN\n', ...
            sum(fa_nan_mask), numel(fa_nan_mask));
ds=cosmo_slice(ds, ~fa_nan_mask, 2);



% Define a channel neighborhood uses meg_combined_from_planar, which means
% that input are planar channels but the output has
% combined-planar channels. Note that with EEG data there is no need to set
% the chan_type, as there is only a single type.
chan_type='meg_combined_from_planar';
chan_count=10;        % use 10 channel locations (relative to the combined
                      % planar channels)
                      % as we use meg_combined_from_planar there are
                      % 20 channels in each searchlight because
                      % gradiometers are paired

% Use cosmo_meeg_chan_neighborhood to define a channel neighborhood,
% with 'count' set to chan_count, and 'chantype' set to chan_type.
% Assign the result to chan_nbrhood. How many searchlight centers are in
% the neighborhood?
% >@@>
chan_nbrhood=cosmo_meeg_chan_neighborhood(ds, 'count', chan_count, ...
                                              'chantype', chan_type);
fprintf('There are %d channel neighborhood centers\n',...
            numel(chan_nbrhood.neighbors));
% <@@<

% Define the frequency neighborhood using 9 time bins
freq_radius=0; % 4*2+1=9 freq bins

% >@@>
% Use cosmo_interval_neighborhood to define the frequency neighborhood.
% Assign the result to a variable 'freq_nbrhood'
% How many searchlight centers are in the neighborhood?
freq_nbrhood=cosmo_interval_neighborhood(ds,'freq',...
                                            'radius',freq_radius);
fprintf('There are %d frequency neighborhood centers\n',...
            numel(freq_nbrhood.neighbors));
% <@@<

% Define the temporal neighborhood using 5 time bins

time_radius=2; % 2*2+1=5 time bines
% Use cosmo_interval_neighborhood to define the temporal neighborhood.
% Assign the result to a variable 'time_nbrhood'
% How many searchlight centers are in the neighborhood?
% >@@>
time_nbrhood=cosmo_interval_neighborhood(ds,'time',...
                                            'radius',time_radius);
fprintf('There are %d temporal neighborhood centers\n',...
            numel(time_nbrhood.neighbors));

% <@@<


% cross the three neighborhoods to get a chan-freq-time neighborhood,
% using cosmo_cross_neighborhood
% How many searchlight centers are in the crossed neighborhood?

% >@@>
nbrhood=cosmo_cross_neighborhood(ds,{chan_nbrhood,...
                                        freq_nbrhood,...
                                        time_nbrhood});
% <@@<

% Count how many features there are, on average, in each neighborhood,
% and store the results in a variable 'nbrhood_nfeatures'. Then print its
% average and standard deviation

% >@@>
nbrhood_nfeatures=cellfun(@numel,nbrhood.neighbors);
fprintf('Features have on average %.1f +/- %.1f neighbors\n', ...
            mean(nbrhood_nfeatures), std(nbrhood_nfeatures));
% <@@<

% divide .sa.chunks in ds into four chunks using cosmo_chunkize

% >@@>
nfolds=4;
ds.sa.chunks=cosmo_chunkize(ds,nfolds);
% <@@<

% Define partitions using cosmo_nfold_partitioner, and use
% cosmo_balance_partitions afterwards. Assign the result to a variable
% 'partitions'
% >@@>
partitions=cosmo_nfold_partitioner(ds);
partitions=cosmo_balance_partitions(partitions,ds);
% <@@<

% Set the measure arguments
measure_args=struct();
measure_args.partitions=partitions;
Summary
-------
MEG data in raw, time-locked and time-frequency formats.

Contents
--------
- subj102_B01_20Hz.fif            raw MEG recording
- preproc.m                       Matlab preprocessing script (based on 
                                  FieldTrip)
- subj102_B01_20Hz_timefreq.mat   Time-locked data (generated by preproc.m)
- subj102_B01_20Hz_timelock.mat   Time-frequency data  (generated by preproc.m)
- LICENSE                         License file 
- README                          This file

Methods
-------
The dataset involved a paradigm with electrical median nerve stimulation with a 
human participant for durations of 2s at 20Hz. Data was acquired at 1khz using 
a neuromag306 system. 
Trial info in the .mat files: 1=pre-stimulus, 2=peri/post-stimulus

License
-------
The contents are made available by Nathan Weisz <nathanweisz |at| me.com> and 
Gianpaolo Demarchi <gianpaolo.demarchi |at| unitn.it> under the Creative 
Commons CC0 1.0 Universal Public Domain Dedication ("CC0"). See the LICENSE 
file for details, or visit 
http://creativecommons.org/publicdomain/zero/1.0/deed.en.

Acknowledgements
----------------
Thanks to Nathan Weisz and Gianpaolo Demarchi for providing this dataset, and 
an anonymous participant for volunteering during the recordings.

Contact
-------
Nikolaas N. Oosterhof <nikolaas.oosterhof |at| unitn.it>
16968 / 92718 features have NaN
There are 102 channel neighborhood centers
There are 18 frequency neighborhood centers
There are 17 temporal neighborhood centers
+00:00:01 [####################] -00:00:00  crossing neighborhoods
Features have on average 77.2 +/- 32.1 neighbors

run searchlight

% only keep features with at least 10 neighbors
% (some have zero neighbors - in particular, those with low frequencies
% early or late in time)
center_ids=find(nbrhood_nfeatures>10);

% instead of using the cosmo_searchlight function, use the faster
% cosmo_naive_bayes_classifier_searchlight.
% As input it takes a dataset, neighborhood, and measure_args.
% Also provide it with the 'center_ids',center_ids argument as well
% to avoid computing results for features without neighbors.

% It may take a while before results for the first fold are computed
result_ds=cosmo_naive_bayes_classifier_searchlight(ds,...
                                        nbrhood,measure_args,...
                                        'center_ids',center_ids);
+00:00:17 [####################] -00:00:00  

visualize results

% deduce layout from output
layout=cosmo_meeg_find_layout(result_ds);
fprintf('The output uses layout %s\n', layout.name);

% map to FT struct for visualization
tf_map=cosmo_map2meeg(result_ds);

% show figure
figure()
cfg = [];
if cosmo_wtf('is_octave')
    % GNU Octave does not show data when labels are shown
    cfg.interactive='no';
    cfg.showlabels='no';
else
    % Matlab supports interactive viewing and labels
    cfg.interactive = 'yes';
    cfg.showlabels = 'yes';
end
cfg.zlim=[0 1];
cfg.layout       = layout;
ft_multiplotTFR(cfg, tf_map);

% Show citation information
cosmo_check_external('-cite');
The output uses layout neuromag306cmb.lay
the call to "ft_selectdata" took 0 seconds and required the additional allocation of an estimated 5 MB
the call to "ft_prepare_layout" took 0 seconds and required the additional allocation of an estimated 0 MB
the call to "ft_multiplotTFR" took 3 seconds and required the additional allocation of an estimated 16 MB
Warning: A value of class "char" was indexed with no subscripts
specified. Currently the result of this operation is the indexed value
itself, but in a future release, it will be an error. 
If you use CoSMoMVPA and/or some other toolboxes for a publication, please cite:

R. Oostenveld, P. Fries, E. Maris, J.-M. Schoffelen (2011). FieldTrip: Open Source Software for Advanced Analysis of MEG, EEG, and Invasive Electrophysiological Data, Computational Intelligence and Neuroscience, vol. 2011, Article ID 156869, 9 pages.doi:10.1155/2011/156869. FieldTrip toolbox available online from http://fieldtrip.fcdonders.nl

N. N. Oosterhof, A. C. Connolly, J. V. Haxby (2016). CoSMoMVPA: multi-modal multivariate pattern analysis of neuroimaging data in Matlab / GNU Octave. Frontiers in Neuroinformatics, doi:10.3389/fninf.2016.00027.. CoSMoMVPA toolbox available online from http://cosmomvpa.org

The Mathworks, Natick, MA, United States. Matlab 9.1.0.441655 (R2016b) (September 7, 2016). available online from http://www.mathworks.com