Customizing common M/EEG plots part 1: the event-related potential/field (ERP/F)

I recently (well, 25 June 2018) taught a data visualization workshop as part of the Language in Interaction consortium at the Donders Institute in Nijmegen. As part of my workshop, I thought I’d go over some fairly simple ways in which people could customize (or dare I say: improve?) some commonly used plots when analyzing MEG or EEG data. Since it would be a waste to have all that preparation for the workshop go to waste on just the small audience that was there that day, I’ll post my code with some explanations here on my blog as well. Today: how to better visualize plots like the event-related potential or field. Or, more generally: how to visualize 1D data over time.

We ran timelocked analysis on an experiment that involves a visual stimulation with onset at time point t = 0. The details are unimportant, we’re going to focus on the visualization. We want to display an evoked field over time, for some sensors which we’ve identified as showing an interesting deflection.

First, load some data (available to download from here; 113 MB) and setup the path (adapt to wherever your own FieldTrip toolbox installation is):

load tl;
addpath H:\common\matlab\fieldtrip;

Plot using default FieldTrip plotting functions:

cfg = []; = {'MLO12', 'MLO22', 'MLP22', 'MLP31', 'MLP32', 'MLT16'};
ft_singleplotER(cfg, tl);

This neatly shows the mean ERF, but many things are lacking in this plot. To improve it from the ground up, we’re going to start by using standard Matlab plotting functions, rather than relying on FieldTrip to do these things for us. In general, I like FieldTrip’s plotting functions a lot for exploring data, but for publication-quality plots ‘manual’ plotting code can be better. (It’s actually on my to-do list to incorporate better plots directly into FieldTrip, but hey, one can only work on so much at the same time…) Let’s start by reproducing the mean ERF with standard Matlab code:

% Find the index of the channels of interest in the data
chaninds = match_str(tl.label, {'MLO12', 'MLO22', 'MLP22', 'MLP31', 'MLP32', 'MLT16'});

% Compute average over channels
erf = mean(tl.avg(chaninds,:), 1);

% Plot
plot(tl.time, erf);

Let’s plot again but add some axis labels. In addition, let’s convert the units to fT (femtoTesla), so we don’t have to deal with the awkward 10-14 scaling factor on the axis.

femto = 1e-15;
erf_femtotesla = erf ./ femto;

plot(tl.time, erf_femtotesla);
xlabel('Time (s)');
ylabel('Magnetic gradient (fT)');

Much better. However, the axes are in a slightly weird position, let’s move the axes such that they cross the origin.

ax = gca(); % this Gets the Current Axis so we can set properties
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
ax.TickDir = 'out';

% Remove the box around the plot, while we're at it:
box off;

% And move the x-axis label to underneath the axis:
ax.XLabel.Position(2) = -30;

Starting to look nice. However, an important aspect is still missing: showing variation! Let’s add error bounds. For this we’ll make use of the excellent Matlab package ‘boundedline’, which is available from GitHub or the Matlab FileExchange:


% First, get the individual trial data, averaged over our channels of
% interest:
trialdata = squeeze(mean(tl.trial(:,chaninds,:), 2)) ./ femto;

% Use the standard deviation over trials as error bounds:
bounds = std(trialdata, [], 1);

% boundedline will replace the call to plot():
boundedline(tl.time, erf_femtotesla, bounds, 'alpha'); % alpha makes bounds transparent

% Add all our previous improvements:
xlabel('Time (s)');
ylabel('Magnetic gradient (fT)');
ax = gca();
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
ax.TickDir = 'out';
box off;
ax.XLabel.Position(2) = -60;

Another option to show variation is to simply display all the trials overlaid. This really highlights the amount of variation present in a dataset, which is something I would in general recommend. (Though at the same time: keep in mind what message you want to convey; data visualization always requires a careful balance between these two principles.) Right now we are working with individual trials, but in a typical ERP/F study you will probably be plotting e.g. individual participants along with the grand mean.

% We'll want the indiviual trials to be plotted very translucently, so use
% this four-element colour specification [red green blue alpha]
linecolor = [0 0 0 0.02];

plot(tl.time, trialdata', 'Color', linecolor);

% Now overlay the mean ERF with standard deviation bounds:
hold on;
boundedline(tl.time, erf_femtotesla, bounds, 'alpha', 'transparency', 0.1);

% Maybe scale back the axes a bit; we don't need to see *every* outlying
% individual trial curve:
ylim([-400 800]);

% Add all our previous improvements:
xlabel('Time (s)');
ylabel('Magnetic gradient (fT)');
ax = gca();
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
ax.TickDir = 'out';
box off;
ax.XLabel.Position(2) = -80;

I think it looks quite a bit better than the original; hope you agree. And even if not, hopefully you learned some stuff about Matlab plotting.

Concatenate subfields from Matlab structure into one data array

Don’t you just hate when you have an MxN struct array in Matlab, where each element contains some PxK data field, and you just really want to have all that data in a single PxKxMxN array? OK, maybe this doesn’t happen to you very frequently, but it has happened to me. In fact, it has happened so often that I decided to write a (very simple) general solution to this problem.

Let me motivate the problem with an example. I was dealing with Granger causality spectra that are stored in FieldTrip structures as follows:

>> allconns(1)

ans =

  struct with fields:

            label: {7×1 cell}
           dimord: 'chan_chan_freq'
    grangerspctrm: [7×7×101 double]
             freq: [1×101 double]
              cfg: [1×1 struct]

and I had one of these structures per participant, per condition of my experiment, and per time-direction (forward or reverse; details don’t really matter here). The structures were all tucked into one structure variable allconns such that

>> size(allconns)

ans =

    36     7     2

To get out the data into one big 6-dimensional array, I first concatenated along the fourth dimension:

catdat = cat(4, allconns.grangerspctrm);

which unfortunately flattens all those 36x7x2 dimensions such that

>> size(catdat)

ans =

     7     7   101   504

However, it’s easy to just get those back by one final reshape:

catdat = reshape(catdat, [7 7 101 36 7 2]);

The above works because Matlab stores data in column-major order.

It’s very straightforward to embed this into a simple function (also available on GitHub):

function x = structcat(y, param)
% STRUCTCAT takes a structure y of arbitrary dimensionality Mx...xN, each
% element of which contains a field named of dimensionality
% Px...xK, and returns the data in those fields concatenated into a
% Px...XKxMx...N array.
% Author: Eelke Spaak, 2018.

siz1 = size(y);
siz2 = size(y(1).(param));

x = cat(numel(siz2)+1, y.(param));
x = reshape(x, [siz2 siz1]);


I hope this is useful to someone who’s encountered the same!