Statistical tests on complex val

Statistical tests on complex valued time–frequency representations for event-related brain dynamics

Eduardo Martínez-Montes (eduardo@cneuro.edu.cu)

Wael El-Deredy (Wael.El-Deredy@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ínez-Montes et al, 2008 Exploring event-related brain dynamics with tests on complex valued time-frequency representations, Statistics in Medicine, 27: 2922-2947.

In the paper, we demonstrated the inadequacy of current measures to distinguish between different types of event-related 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 time-frequency 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 Time-Frequency (Clouds of complex numbers)

The distribution of the energy of the signal in each time and frequency can be calculated with a Time-Frequency 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 Time-Frequency. 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 event-related 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 non-zero 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 de-synchronisation)

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 non-parametric test for the presence of a nonzero mean or net vector.

$D:\Work\Papers\StatMed\2ndREvision\finalfigures\Figure2rev.bmp$

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 T-test on the mean (compared to pre-stimulus 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 pre-stimulus 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 T-like 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 T-stat, that we cannot use the T-stats tables to obtain the p values since the measures of the changes are not normal distributed. So we have to use a non-parametric test known as the false discovery rate (or FDR) [Benjamini and Hochberg. 1995, J. Roy Stat Soc B 57:289-300) which is commonly used in fMRI data analysis [Genovese et al 2002, NeuroImage 15:870-878] 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 time-frequency 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, T-stat, 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 pre-stim duration say [1 50] or [51 100] corresponding to the first or last 100 ms of the pre-stim period respectively. Indeed, you can choose your basline_period to be some time after the stimulus say [901 1000] corresponding to 1800-2000 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 time-frequency 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ínez-Montes et al, 2008)

‘induced’ tests the variance (equation 4 in the paper)

‘ITC’ gives the T-like statistic for the inter-trial phase coherence or phase locking factor (Makeig et al Science 2002; 295:690–694. Tallon-Baudry 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 distribution-free method developed by Efron which is used in Martínez-Montes 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 p-value but it plays the role of controlling the Type I error. It is known as q-value 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).

T-stat This is the raw T-like 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 p-values (if fdr_type='globalfdr', assuming a T distribution with df=No. of trials – 1) corresponding to the T-like statistics computed

threshold 2-length 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: T-eigenvalue')

colorbar

You can also plot the image of the actual measure or of the un_thresholded T like stats.

Figure

imagesc(time, freq, EIG);

figure

imagesc(time, freq, T_Eig);