Statistical tests on complex
valued time–frequency representations for eventrelated brain dynamics Eduardo MartínezMontes (eduardo@cneuro.edu.cu) Wael ElDeredy (Wael.ElDeredy@manchester.ac.uk) Since our presentation at ICONX in Bodrum in September, a number of colleagues have asked us to make available the Matlab routines from the paper by MartínezMontes et al, 2008 Exploring eventrelated brain dynamics with tests on complex valued timefrequency representations, Statistics in Medicine, 27: 29222947. In the paper, we demonstrated the inadequacy of current measures to distinguish between different types of eventrelated brain dynamics (ERBD). In particular, we showed that current measures of phase locking cannot separate true events of phase reorganisation from changes in mean activity and consequently cannot be considered an accurate measure of phase resetting. We made the case for using separate measures based on the complex valued timefrequency representations to test for each of these possible events and demonstrated four such measures for evoked and induced activity and phase reorganisation. We have packaged these tests into a set of simple Matlab routines and included data to reproduce the figures in the paper. We have tried to add as much help as possible to explain the routines and the syntax of the functions, but all feedback is welcome. Transformation from domain
Time to TimeFrequency (Clouds of complex numbers) The distribution of the energy of the signal in each time and frequency can be calculated with a TimeFrequency transformation using Fourier or equivalent transform: STFT, Morlet Wavelet, Hilbert, Gabor, etc. The result of the transformation are complex coefficients _{}, whose moduli  _{} is a measure of the amplitude of the oscillations (Figure 1) and whose argument arctan( Im(_{})/ Re (_{}) ) is a measure of their phases. Each trial (i) at a given time (t) and frequency (f) is a dot in a complex cloud, forming a net vector with an average amplitude and averaged phase (Figure 2). We shall see later that it is this net vector that confounds measures of amplitude and phase locking. However, we will use the features of this complex cloud to characterise different ERBD scenarios and use it to define new measures to statistically test for each of these scenarios.
Figure 1 From Time to TimeFrequency. The familiar TF plot only displays changes in power (in colour) and phase information is lost . Figure 2 . Complex cloud and net vector Scenarios of
eventrelated brain dynamics Additive event related activity (the classical ERP view)
Figure 3 Additive ERP. (Additional) stimulus evoked activity is superimposed on the ongoing EEG. This type of activity is equivalent to shifting the cloud of complex coefficients from the origin after stimulus presentation; that is, the appearance of a net vector and phase.
Figure 4 Additive ERP shifts the cloud of coefficients such that a mean vector emerges. Phase resetting of
ongoing activity
Figure 5 Partial phase resetting. The stimulus causes a reorganisation of the ongoing EEG, without additional activity. The evoked response obtained by the usual “averaging of trials” procedure is the same as the one obtained when an additive ERP is superimposed to ongoing activity (figure 3). This type of activity causes the cloud of complex coefficients to concentrate around certain phase without change in amplitude and therefore without change in the EEG power spectrum. However, here too a nonzero net vector emerges.
Figure 6 Partial phase resetting (PPR) results in an effective shift of the cloud of coefficients around a certain phase, such that a mean vector emerges. Induced activity (event
related synchronisation and desynchronisation) Here the event related changes in the EEG power spectrum are not due to changes in the mean or to partial phase resetting, but due to changes in the amplitude of the oscillations, which is measured by the variance _{} of the complex coefficients _{}. The variance _{} is related to the EEG power spectrum, so this activity is usually found by comparing power in the pre and post stimulus EEG spectrum. Increases of the variance _{} are termed “synchronization” and decreases “desynchronization”.
Figure 7 Induced activity does not change the orientation of the complex cloud, only its variance Phase locking factor
(or intertrial phase coherence – ITC) cannot be an accurate measure of phase
resetting. Why not?
Intertrial Phase Coherence (ITC) measures the uniformity of the distribution of angles, with respect to the origin  not with respect to the centre of the cloud, such that it cannot truly distinguish between whether the net vector is due to additional evoked activity or phase resetting. Even if we remove the mean, ITC still cannot distinguish between those two cases (Figure 8). _{} For ITC to be a measure of phase resetting we have to assume that phase resetting is the only phenomenon that occurs in neural systems. In fact we have shown that ICT is effectively a nonparametric test for the presence of a nonzero mean or net vector. Figure 8 ITC cannot separate additive ERP from phase
resetting Tests on the complex cloud Necessary conditions For additive ERP: It has to survive a Ttest on the mean (compared to prestimulus baseline) For PPR: It has to specifically test for changes in the shape of the complex cloud, and to survive the subtraction of the mean vector (compared to prestimulus baseline). Statistical tests on
the complex cloud We are proposing to use separate tests for each ERBD scenario. These tests operate on the entire cloud of complex numbers. To do that we have derive some suitable complex statistics similar to the ones we are familiar with in real numbers. From the measure of a particular ERBD scenario, we propose to compute a statistic that accounts for deviations from their values in the baseline period, normalized by the standard deviation in this period. Thus, we will refer to the tests as Tlike statistics. We test for each type of ERBD separately. The argument being the presence of one type of ERBD, say additive ERP, does not rule out the presence of other types, e.g. phase resetting or induced activity as well. These tests are different enough from the familiar Tstat, that we cannot use the Tstats tables to obtain the p values since the measures of the changes are not normal distributed. So we have to use a nonparametric test known as the false discovery rate (or FDR) [Benjamini and Hochberg. 1995, J. Roy Stat Soc B 57:289300) which is commonly used in fMRI data analysis [Genovese et al 2002, NeuroImage 15:870878] and its variant the local false discovery rate (lFDR) [Efron and Tibshirani 2002, Genet. Epidemiol. 23 70–86] to assess the significance of the tests and simultaneously correct for the multiple comparisons in the timefrequency map. The code we included here has the option to use either the global (Benjamini) or the local (Efron) FDR. We are also working on further improvement of this statistic specifically for complex TF data. The detailed description of these tests is in the paper. Here we will focus on the structure of the Matlab code and how to use it. Installing and
running the toolbox
Download the toolbox in .RAR format Unzip the toolbox, it will create a folder ERBD with subfolders. In Matlab, go to the ERBD folder >> cd ~\ERBD Type >> link_me This will add the folders to the Matlab path Main function [measure, Tstat, fdr, threshold] = TestERBD(data, baseline_period, ERBD_type, fdr_type, fdr_cutoff) Inputs data 3D matrix of complex coefficients with the dimensions Trials x Time x Frequency
baseline_period: index to prestim or baseline period to which the test is compared. Must be two numbers indicating the first and the last sample of this period. For example, if the trial is 1100 data points corresponding to 200 ms prestimulus and 2000 ms post stimulus (this means the EEG was sampled at 500Hz), then Baseline_period is [1 100]. Of course you can choose to set the baseline to only a smaller part of the prestim duration say [1 50] or [51 100] corresponding to the first or last 100 ms of the prestim period respectively. Indeed, you can choose your basline_period to be some time after the stimulus say [901 1000] corresponding to 18002000 ms post stimulus in the above example. We highly recommend to avoid using times near the onset and the borders as part of the baseline, due to nuisance effects of timefrequency transformations (such as Wavelets, Hilbert, Fourier). Important thing you MUST specify a baseline against which the possible ERBD are tested. ERBD_type The type of event you want to test for. This is a text variable as follows: ‘evoked’; ‘induced’; ‘ITC’; ‘phase_eig’ or ‘phase_r2’ Each option will invoke the appropriate statistical test: ‘evoked’ tests the mean (equation 3 in MartínezMontes et al, 2008) ‘induced’ tests the variance (equation 4 in the paper) ‘ITC’ gives the Tlike statistic for the intertrial phase coherence or phase locking factor (Makeig et al Science 2002; 295:690–694. TallonBaudry et al J. Neuroscience 1996; 16:4240–4249). (Equation 2 in the paper). ‘phase_eig’ tests for phase reorganisation using the the test for similarity of the eigenvalues of the covariance matrix between real and imaginary part of the complex coeffficients (equation 5 in the paper). This can be seen as a measure of generalizedcorrelation between real and imaginary parts. ‘phase_r2’ tests for phase reorganisation using the second trigonometric moment, which measures uniformity of phases in the bimodal distribution that appears after removing the net vector produced by a phase resetting (equations 6 and 7 in the paper) fdr_type Optional FDR method for finding significant thresholds. 'globalFDR' uses Benjamini’s method, used in fMRI data (Genovese et al 2002). 'localFDR' uses the distributionfree method developed by Efron which is used in MartínezMontes et al, 2008, and is the default. fdr_cutoff Optional Expected fraction of false positives among all points assessed as significant. This is not a pvalue but it plays the role of controlling the Type I error. It is known as qvalue in the FDR terminology and is usually taken below 0.05 for the 'globalFDR' and below 0.2 for the 'localFDR'. These values are taken as default. Outputs Measure Depending on which test you choose (see ERBD_type above) this is the actual value of the measure (i.e the means, variances or phase_r2 if you choose ‘evoked’, ‘induced’ or ‘phase_r2’, respectively). Tstat This is the raw Tlike stats of each test (again depending on the ERBD_type), as in equations 3, 4, 5 or 7. fdr Matrix (frequency x time points) with the fdr probability (if fdr_type='localfdr') or with the pvalues (if fdr_type='globalfdr', assuming a T distribution with df=No. of trials – 1) corresponding to the Tlike statistics computed threshold 2length vector. Thresholds (left and right) corresponding to fdr_cutoff. This is used for plotting the probability map of significant regions on the TF plot.
Paper Data and
Figures In the folder “Paper_Data and Figures” we have included scripts to regenerate the key figures of the paper. If you intend to use these scripts and want to use our data, you’ll need to download it separately. [download here.] There are three files corresponding to the simulated ERP and PPR and the real data. Download the data in the folder “Paper_Data and Figures”. The scripts are named after the figures included: Figure3; Figure4; Figure5; Figure6 and Figure8. If you save your data in the same way as ours, you can easily change the script for figure 8 to load your own data. Examples of plotting
results [EIG, T_Eig, fdr_Eig,
Th_Eig] = TestERBD(WAcomplex, [1 81], 'Phase_eig', 'localFDR', 0.05); %To plot the
thresholded image Thresh_image= zeros(size(Teig)); % Initialises a new
thresholded image Thresh_image (T_Eig <
Th_Eig(1)) = T_Eig(T_Eig< Th_Eig(1)); % This includes the
significant negative values Thresh_image (T_Eig >
Th_Eig(2)) = T_Eig(T_Eig> Th_Eig(2)); % This includes the
significant positive values figure imagesc(time, freq, Thresh_image);
% this assumes you
have kept a vector of the time and frequency values % if not you can also
use imagesc(Thresh_image); axis('xy'); ylabel('Frequency (Hz)');xlabel('Time (ms)'); title('Phase organization:
Teigenvalue') colorbar You can also plot the image of the actual measure or of the un_thresholded T like stats. Figure imagesc(time, freq, EIG);
or figure imagesc(time, freq, T_Eig);

