Artifact Rejection
Of course you as an experimenter will try your very best in getting the cleanest data possible. Nevertheless some artefacts are inevitable and they come in different flavours, e.g.:
- Physiological artefacts such as blinks or ECG
- Environmental noise that will mostly affect your magnetometers
- Channel noise such as channel jumps (usually affecting a few trials) or "bad channels" (that are basically noisy over the entire experiment)
With regards to physiological artefacts some are easier to control (e.g. eye blinks) whereas others are inevitable (e.g. you can't ask the participant to stop the heart beating). While the most efficient way to remove occasional artefacts (such as blinks) is simply to remove those trials from the data, you can consider using e.g. ICA to remove other constant artefacts such as ECG. Whether you do so depends on to what extent you think that this factor is a problem in your analysis (e.g. do you think that there will be a systematic ECG difference between conditions / groups; in particular when such artefacts are time-locked to a specific event they are more problematic).
TODO: FOR NOW WE HAVE NOT TESTED ICA ON MIXED SENSOR TYPES AND WE WILL PROCEED SIMPLY BY REMOVING BAD TRIALS
Our Neuromag system comes with different challenges with regards to this preprocessing step:
- Mixture of different sensor types
- Instability of the system, meaning that it is quite variable on a day to day basis which channels will be bad.
With regards to the different sensor types you can of course choose to only read data from one sensor type in the first place (by setting cfg.channel to 'MEGGRAD' or 'MEGMAG' in ft_preprocessing), which is not uncommon (many papers only publish gradiometer data). But at the moment the MEG-group tries to keep the information and data together as long as possible. The instability issue poses particularly big problems when you are interested in sensor level statistics, since the Fieldtrip statistic functions will discard a missing channel in one participant from all other participants. This can lead to a significant loss of information as well as statistical power.
While automatic thresholds or template matching procedures can be used to detect trials with certain artefacts, we think that the best way at the end of the day is still by visual inspection. Since this really is the only step in which you directly interact with the data, you should be careful during this step.
Warning
Some users reported that at this stage they're experiencing 'Out of Memory' issues. This is mainly due to the fact that ft_rejectvisual is making a temporary copy of the data. On decently modern computers (i.e. at least 4 Gb of RAM), this shouldn't be a problem, if configured properly.
As a general reminder, for doing MEG data analysis, you should have available at least 4 Gb of RAM, in the sense that your (Operating System + Matlab) setup should be able to access that amount of RAM.
For that you need to have everything with 64 bits support, i.e. the hardware architecture should be 64 bit capable (all modern machines), the OS as well (i.e. all Linuxes, all MacOsX, Windows XP 64-bit, Windows 7), and you should have also the 64 bit version of Matlab installed. If one piece of the 64 bit chain is missing, you will be unable to user effectively more than 2 GB of RAM. If you see that you have only the 32 bit version of Matlab installed, either install the 64 bit as well, or open a ticket and ask the IT's to do it for you.
For more info, please refer to
cfg = [];
cfg.showlabel='yes';
cfg.channel= {'MEG'};
cfg.latency = [-.5 2]; %home in on relevant period ... choose accordingly. leaving out this option then uses entire epoch
cfg.method = 'summary';
cfg.gradscale = 0.04; %factor to bring Grad and Mag data onto same scale
cfg.layout = 'neuromag306mag.lay';
data_raw = ft_rejectvisual(cfg,data_raw);
This is what you should see (useful to have a large monitor):
In the 'summary' view (check out documentation for other viewing options) you can see three figures. The outer right displays the relevant metric (default is variance, but you can choose within the window) for each channel across the entire epoch (see comment on choosing relevant periods above). In the lower panel the same metric displays the metric for all trials. The upper left figure is then the channel x trial depiction in color coding.
Basically a good strategy is first to remove outlying trials, simply by holding the mouse button down and drawing a window over the relevant trials. In this depiction there are no shockingly outlying trials, but I will remove the top 4 trials in any case. You can see these trials removed in the color coded image immediately. Bad channels can be detected either in the upper right window or in the color coded window as horizontal lines with extremely 'hotter' colors. There may be 1-2, but nothing shocking in this dataset from these depiction. However one can easily identify a 'dead' channel, i.e. the one sticking at 0 in the upper right image. Remove it accordingly as you did for the trials. Trials and channels removed will be displayed on the right hand side. You can toggle them in and out of the data by simply typing the according channel / trial in the relevant fields.
Now that we have done a very course artefact rejection, removing the most prominent stuff it may be useful to look at some trials in detail and remove remaining bad trials by hand. In principle you could of course look at all trials, but in practice it is usually sufficient to focus on the ones that have high values on the relevant metric. For doing the in detail look at the data it is useful to switch on the Figure Toolbar (under View) and using the zooming tools. One can then e.g. zoom in the color coded field into epochs and channels with intensive red color, e.g. trial number 72. Typing 72 into the relevant field will give you a layout of the time series of all channels. Note that for this purpose we only use the magnetometers because here (i.e. before combining the gradients) they are easier to interpret. You should be able to see see something that spreads from frontal to central sensors reaching its highest amplitude there (~2 pT). This layout shows no clear polarity reversal. It's unlikely to be an eye-movement, but about its origin can only be speculated; @GIANPAOLO: COULD IT BE A STIMULATOR ARTEFACT?). You can in this layout depiction also choose channels to display and within this display select time periods in order to see 2D-topographies (broken on my Mac but working under Linux). Since it is highly unlikely to have such high amplitude events with such topographies being genuine brain data I would suggest deleting the trial. So in the toggle trial window type 72 and hit return. Proceed with other trials that look suspicious. Beware that the color coding adapts to min/max so there will always be deeply red trials (similarly when plotting trials, then look at the scale ... MEG / EEG raw data always looks noisy). When you are done click Quit (do not just close the window).
Visual artefact rejection requires practice (seeing lots of data) and you will get experienced with time. Ask some MEG-experienced person around you to look at trials. The worst thing that can happen is redoing it.
At the end it is sometimes worth looking at the other metrics, since they are sensitive to different aspects of artefacts.
In the end I removed ~20 trials and 3 channels (the dead one and another noisy one detected using the 'max' and 'range' metric).
Beware: The computation of the metric sometimes does take time. Don't be impatient and start clicking around wildly.
You can now check data_raw and you will see that trials have been indeed removed. Importantly the 'label' field now only has 303 instead of 306 cells.
After, having removed bad trials and channels we know interpolate rejected channels using a home-made wrapper function to ft_channelrepair. This function currently is in a beta stage so use with care. Normally, the best practice in MEG is simply to discard the bad sensors. However the Elekta system (contrary to CTF or 4d) is unfortunately quite unstable, meaning that bad channels can vary from day to day or sometimes even from block to block. This will introduce issues if you intend to do sensor level statistics, which is now common practice in MEG, since a channel lost in one participant is basically lost for the entire statistics. For these purposes (and also because ft_combineplanar has troubles if gradiometer pairs are incomplete) we decided to add channel interpolation, however keeping track of the original channels before interpolations, so that only these are then used for source analysis.
The call to the wrapper function is:
%at the moments runs only with defaults; will add flexibility later cfg=[]; data_raw = cimec_fixchannels(cfg,data_raw);
The function basically automatically detects the labels of the removed channels and replaces it in the default version by the average of its neighbours (of the respective sensor type; Bravo Gianpaolo!) taking into account the different sensor types. The "repaired" data will now include a 'label' field with all 306 channels and an 'origLabel' field that includes all channels without the discarded ones. This field will become important for source analysis later.
Update: Apart of this automatic checking the function can now also take a hand-written specification of bad channels.
TODO: Make function work when channel was removed during data acquisition.
You can now again check data_raw and you will see that the 'label' field contains 306 cells again.
As stated above: Apart from assuring the best possible data during data acquisition, this is the only time you directly intervene with the data. Do this step with great care!
Now we are ready to continue with sensor level analysis.
;){kind=link}