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):
addpath H:\common\matlab\fieldtrip;
Plot using default FieldTrip plotting functions:
cfg = [];
cfg.channel = {'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:
chaninds = match_str(tl.label, {'MLO12', 'MLO22', 'MLP22', 'MLP31', 'MLP32', 'MLT16'});
% Compute average over channels
erf = mean(tl.avg(chaninds,:), 1);
% Plot
figure();
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.
erf_femtotesla = erf ./ femto;
figure();
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.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():
figure();
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.
% this four-element colour specification [red green blue alpha]
linecolor = [0 0 0 0.02];
figure();
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 <em>every</em> 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.
Dear Eelke, Thank you so much for this post! It really helped me to better visualize my grandaverage ERP results. I just wanted to show my gratitude :)