Multivariate Statistical Analysis in Neuroscience

Contents

Preface

Acknowledgements

1. INTRODUTION TO ELECTROENCEPHALOGRAPHY (GENERAL ASPECTS)
1.1 FUNDAMENTALS OF EEG MEASUREMENT
1.1.1 Activity of the brain
1.1.2 Electroencephalography (EEG)
1.1.3 Quantitative electroencephalography (qEEG)
1.1.4 Frequency and amplitude of the signal
1.1.5 International 10-20 system

2. DATA MINING AND STATISTICAL ANALYSIS
2.1 PRE-PROCESSING PROCEDURES
2.1.1 EEGLAB: statistical software for electro-physiological data analysis
2.1.2 Importing channel location: information about the electrodes placement
2.1.3 Filtering data to minimizing the introduction of artifacts
2.1.4 Extracting data epochs and removing baseline values
2.2 CHANNEL DATA ANALYSIS
2.2.1 Channel data scroll: visualization, normalization and channel rejection procedure
2.2.2 Channel spectra and associated topographical maps
2.2.3 ERP and associated topographical maps
2.2.4 Time/frequency decomposition
2.2.5 Cross-coherences computation
2.2.6 Channel summary
2.2.7 Rejecting artifacts in continuous and epoch data
2.3 COMPONENT DATA ANALYSIS
2.3.1 Independent Component Analysis
2.3.2 ICA Algorithms
2.3.3 Component data scroll
2.3.4 Component spectra and associated topographical maps
2.3.5 Time/frequency decomposition
2.3.6 Computing cross-coherences
2.3.7 Component summary
2.3.8 Rejecting based on independent data components
2.4 MULTIPLE SUBJECT DATA PROCESSING
2.4.1 Channel statistics
2.4.2 Component statistics
2.4.3 Clustering procedure
2.4.4 Preparing to cluster with PCA method
2.4.5 Clustering
2.4.6 Component clusters visualization
2.5 STATISTICAL PROCEDURES
2.5.1 Parametric and non-parametric statistics
2.5.2 Paired/unpaired samples
2.5.3 Re-sampling methods
2.5.4 Multivariate methods (PCA vs ICA)
2.5.5 Correcting for multiple comparisons

3. EXPERIMENTAL STUDY DESIGNS
3.1 ELECTROCORTICAL ACTIVITY IN DANCERS AND NON-DANCERS LISTENING TO DIFFERENT MUSIC GENRE AND DURING IMAGINATIVE DANCE MOTOR ACTIVITY
3.1.1 Abstract
3.1.2 Introduction
3.1.3 Materials and Methods
3.1.4 Statistical analysis
3.1.5 Results
3.1.6 Discussions and conclusions
3.2 ELECTROCORTICAL ACTIVITY DURING MONOSYNAPTIC REFLEX IN ATHLETES
3.2.1 Abstract
3.2.2 Introdution
3.2.3 Materials and methods
3.2.4 Statistical analysis
3.2.5 Results
3.2.6 Discussion and conclusions
3.3 MONITORING OF ELECTROCORTICAL ACTIVITY FOR EVALUATION OF SEASICKNESS
3.3.1 Abstract
3.3.2 Introdution
3.3.3 Materials and methods
3.3.4 Statistical analysis
3.3.5 Results
3.3.6 Discussions and conclusions
3.4 ELECTROCORTICAL ACTIVITY IN DIFFERENT BODY POSITIONS
3.4.1 Abstract
3.4.2 Introdution
3.4.3 Materials and methods
3.4.4 Statistical analysis
3.4.5 Results
3.4.6 Discussions and conclusions
3.5 ELECTROCORTICAL ACTIVITY IN ATHLETES AND NON-ATHLETES DURING BODY BALANCE TASKS
3.5.1 Abstract
3.5.2 Introdution
3.5.3 Materials and methods
3.5.4 Statistical analysis
3.5.5 Results
3.6 ELECTROCORTICAL RESPONSES IN VOLUNTEERS WITH AND WITHOUT SPECIFIC EXPERIENCE WATCHING MOVIES INCLUDING THE EXECUTION OF COMPLEX MOTOR GESTURES
3.6.1 Abstract
3.6.2 Introduction
3.6.3 Materials and methods
3.6.4 Statistical Analysis
3.6.5 Results
3.6.6 Discussions and conclusions

4. OTHER INTERESTING THINGS
4.1 COMPARISON BETWEEN CLINICAL DIAGNOSTIC CRITERIA OF SLEEP BRUXISM AND THOSE PROVIDED BY A VALIDATED PORTABLE HOLTER
4.1.1 Abstract
4.1.2 Introduction
4.1.3 Materials and Methods
4.1.4 Statistical Analysis
4.1.5 Results
4.1.6 Discussion and conclusions

REFERENCES

Preface

Electroencephalography, commonly called 'EEG', estimates through the application of electrodes, the electrical activity of the brain (which is the sum of the electrical activity of each neuron). In recent years, with the goal of making more reliable the EEG, many researchers have turned their interest in the development of tools, methods and software. This thesis describes some best procedures for the experimental design, data visualization and descriptive or inferential statistical analysis. The application of statistical models to single or multiple subjects study-design are also described, including parametric and non-parametric approaches. Methods for processing multivariate data (PCA, ICA, clustering) were described. Re-sampling methods (bootstrap) using many randomly software-generated samples were also described. The aim of this work is to provide, with statistical concepts and examples, information on the qualitative and quantitative approaches related to the electroencephalographic signals. The work consists into three parts: INTRODUTION TO ELECTROENCEPHALOGRAPHY (GENERAL CHARACTERISTICS); DATA MINING AND STATISTICAL ANALYSIS; EXPERIMENTAL STUDY DESIGNS. The six works included in the section called “EXPERIMENTAL STUDY DESIGNS” analyze EEG alterations in the protocols: Electrocortical activity in dancers and non-dancers listening to different music genre and during imaginative dance motor activity; Electrocortical activity during monosynaptic reflex in athletes; Monitoring of electrocortical activity for evaluation of seasickness; Electrocortical activity in different body positions; Electrocortical activity in athletes and non-athletes during body balance tasks; Electrocortical responses in volunteers with and without specific experience watching movies including the execution of complex motor gestures. In the section called “OTHER INTERESTING THINGS” were included one work that analyze EMG (electromyography) alterations in pathological and healthy subjects in the protocol: Comparison between clinical diagnostic criteria of sleep bruxism and those provided by a validated portable holter. The described procedures can be used for clinical trials, although the studies proposed in this work do not refer to samples from pathological subjects. With its multi-specialist approach, through many theoretical and practical feedback, this work will be useful for specializing in neuroscience, statistics, engineering or physiology.

Acknowledgements

I like to express my thanks to: Department of Brain and Behavioral Sciences, Unit of Medical and Genomic Statistics, University of Pavia, Strada Nuova, 65, 27100, Pavia, Italy; Department of Medical Sciences Motor Science Research Center SUISM University of Turin P.za Bernini 12 10143 Torino, Italy; Department of Surgical Sciences, Specialization School of Orthodontics, Dental School, University of Torino, Via Nizza 230, 10126, Torino, Italy; Department of Surgical Sciences, Gnathology Unit, Dental School, University of Torino, Via Nizza 230, 10126, Torino, Italy.

I like to express my thanks to: OT BioLab (version 1.8, OT Bioelettronica, Turin, Italy) for raw signal recording; MATLAB and Statistics Toolbox Release 2012b (the MathWorks, Inc., Natick, Massachusetts, United States) for data processing; EEGLAB software (Swartz Center for Computational Neuroscience, University of San Diego, California) for data visualization and statistical analysis; The R statistical package (version 3.0.1, R Core Team, Foundation for Statistical Computing, Vienna, Austria) for statistical analysis.

I would like to express my deepest gratitude to my advisor, Prof. Alberto Rainoldi, for his excellent guidance and providing me with an excellent atmosphere for doing research. He has inspired my scientific research interest.

I extend my genuine thanks with gratitude to my advisor, Prof. Luisa Bernardinelli. She has inspired my scientific research interest.

I would like to express my special appreciation and thanks to my advisor Dr. Marco Ivaldi, for his motivation, enthusiasm, knowledge and for encouraging my research, allowing me to grow as a research scientist. “You have been a mentor for me”.

Thanks also to all the trainees Eleonora Fiorenti, Francesca Pretari, Valentina Frison, Sara Peracchione, Valentina Verzoletto and Michela Carlucci for providing a good atmosphere in our department.

I would like to thank affectionately the friends of this adventure Simona De Summa, Serena Martire, Alessandra Dentamaro and all teachers.

I would also like to thank my parents, and two brothers. They were always supporting me and encouraging me with their best wishes.

Last but not the least, I would like to thank my girlfriend, Dr. Daniela Donatiello. She was always there cheering me up and stood by me through the good times and bad. “Your support has been fantastic”.

Giovanni Cugliari

1. INTRODUTION TO ELECTROENCEPHALOGRAPHY (GENERAL ASPECTS)

1.1 FUNDAMENTALS OF EEG MEASUREMENT

The aim of this section is not to explain all of the features referred to electroencephalographic signals, but to provide the needed inputs to address to mining and analysis procedures.

However, in order to proceed with the second section, you must have some information about the characteristics of the EEG signals (frequency, amplitude) and importing procedures (electrodes placement on the scalp).

1.1.1 Activity of the brain

Bioelectrical phenomena, that occurs in the cerebral cortex, determines the generation of the electric potential that can be recorded via electrodes placed to the scalp (electroencephalography, EEG) or directly on the cortical surface (electrocorticography, ECOG). The same phenomena also give rise to very weak magnetic fields registered by sensors (magnetoencephalography, MEG).

1.1.2 Electroencephalography (EEG)

Electroencephalography is the recording of brain electrical activity by sensors. The electrodes are arranged on the surface of the head using a suitable amplifying equipment. The EEG signal reflects the electrical events of the skull. These events include cerebral post-synaptic potentials, action potentials, electrical signals of the skin, muscles, blood vessels and eyes.

Some EEG waves are associated with certain states of consciousness and brain-specific pathological conditions (like epilepsy). In some cases, researchers are more interested in the EEG analysis related to certain psychological events. These EEG associated with events (external or internal) are called “event related potentials” (ERPs). Since the EEG amplitude of signal decreases as it spreads from its point of origin, a comparison of the signals recorded from different points of the scalp can sometimes indicate the origin of any particular waves. This explains how the EEG activity is recorded in many points of the scalp. One of the most common potential event related is the sensory evoked potential, the modification of the EEG caused by momentary presence of a sensorial stimulus.

The EEG signal has two components: the response to stimulus and the contemporary background activity. Obviously, the signal is an important part of each trace recorded, while the background activity isn’t so important. The problem in sensory evoked potentials recording of is that the background is usually great as to obscure the event related signal. The analysis of evoked potentials takes into account the peaks or waves present in the average EEG. Each wave is characterized by direction, positive or negative, and latency. These waves are called stem-encephalic potential (far-field potentials) because, although recorded on the scalp, they originate in the sensory nuclei of the brain stem.

Even if the electroencephalography has an excellent temporal resolution, in a first time showed very modest results as to the spatial resolution. With the electroencephalographic conventional procedures, we can estimate the signal source only approximately. The latest techniques, that using sophisticated software, can locate precisely the signals source.

1.1.3 Quantitative electroencephalography (qEEG)

The quantitative electroencephalography (qEEG) is different from clinical EEG due to the mathematical analysis of the brain waves also in not pathological conditions. This procedure is not invasive and allows better performance compared to other neuroimaging techniques, such as magnetoencephalography (MEG), as regards of the time resolution of the recorded signal. Another advantage of this device is that it reduced the required spatial, no larger than a personal computer, and much less expensive than other instruments.

QEEG is accurate, functional, fast and reliable (2 minutes are reliable to 96%). The normal brain activity includes electrical activity that, although attenuated, is still measurable on surface of the scalp and his magnitude is some dozens of microvolt. The continuous oscillation of these electric currents creates the phenomenon of brain waves. The track recorded is the representation of different signals recorded from electrodes and analyzed in a differential way with reference medially of the skull. This type of representation is defined monopolar, since all the signals refer to the signal of a single electrode (therefore defined reference).

The EEG is the representation of the postsynaptic potentials that born in the cortex of the brain by synchronous activity of about 105 neurons; the number is so high because the signal must pass through several layers of non-neutral tissue including the meninges, the intermediate liquids, the bones of the skull and skin before being picked up by the electrode.

The reference electrode should be positioned at a certain distance from active electrodes, it can be placed on the scalp, cephalic reference, or in other body regions (mastoid, earlobes, backs of hands) electrically inert or anyway without electrical activity named non-cephalic reference. The electrical activity of the brain is the sum of several waves at different frequencies and generated by specific signals related to particular tasks (sensorial, movement related or cognitive) in which the subject is involved during the recording.

From the comparison between the spontaneous activity and its variation during the induced activity is possible to identify in real time the areas of change electrical activity. The strength of this type of analysis is the time resolution, in fact is possible to obtain the information every millisecond of electrical activity at the intra-extra cortical level of recordings; its major limitation is the spatial resolution, the electrodes pick up only the current that reaches the surface of the skull and just due reconstruction signal algorithms is possible to locate the source of signal.

The elementary functional units of the cerebral cortex are composed of columnar clusters of neurons organized with perpendicular orientation to the surface of the cerebral cortex. The potentials measured by EEG are associated to the flow of electrical current from the brain to the scalp.

1.1.4 Frequency and amplitude of the signal

Alpha waves (8-12 Hz) are characteristics in the wakefulness and rest condition, but absent in sleep (except for the REM stadium). When subjects show a brain activity increased, the presence of beta rhythm is possible to recognize.

Beta (16-32 Hz) presents an average voltage of 19 uV (8-30 uV). Beta waves are dominant in subjects with open eyes, but also in states of alertness and REM sleep.

The theta rhythm is dominant in newborns and in people in emotional stress (4-8 Hz), has an average voltage of 100 uV.

Finally, the delta waves have a frequency between 0 and 4 Hz and an average voltage of about 150 UV; are not present in physiological conditions in the wakefulness in adulthood, although they are also prevalent in childhood and generally appear in anesthesia and in some brain or general metabolic diseases, such as azotemia. The delta waves are characteristic of REM sleep (slow wave sleep). In different stages of sleep are mainly theta and delta waves, glimpses of alpha activities and, rarely, beta activity.

The resulting signal is a mixture of all these waves, with percentages greater for beta waves and alpha (90-95%) and percentages smaller for theta and delta waves (the first 3-4% and the second 0,5-1%).

illustration not visible in this excerpt

Figure 1 – EEG rhythms show an example of signal filtered with a passband filter for the indicated frequencies (Kevin Roebuck, 2012)

illustration not visible in this excerpt

Figure 2 – Example of power spectrum, after independent component analysis (ICA, )with specifics frequency bands:. Frequency bands subdivision has been highlighted

illustration not visible in this excerpt

Figure 3 – Subdivision of the EEG frequency bands with the description of the related frequencies and several literature references associated to the band (Kevin Roebuck, 2012)

1.1.5 International 10-20 system

The electrodes are applied on the scalp according to the standard International 10-20 system, 10% or 20% refer to the percentage of the distance that separates the electrodes, this distance usually varies from 30 to 36 cm with great interpersonal variability between two landmarks points of cranial "inion" (prominence at the base of the occipital bone) and "nasion" (hanging above the nose). Electrodes (16+3) and mass are placed, along five lines. For each electrode on the scalp there is a code reference. The acronyms identify the position of an electrode and they are formed by one/two letters identifying the region of the cortex (Fp: frontopolar; F: frontal, C central, P: parietal, T: temporal, O: occipital) and a number (or z) that identifies the hemisphere (odd numbers: left, even numbers: right; z: midline).

illustration not visible in this excerpt

Figure 4 – Electrode placement in International 10-20 System. Azure highlighted shows the typical electrodes sequence of used in EEG recording (EEGLAB)

2. DATA MINING AND STATISTICAL ANALYSIS

2.1 PRE-PROCESSING PROCEDURES

This section has the aim to explicate data mining procedures first of all; secondly, to show the signal analysis procedures; finally, to provide the statistical methods to analyze multi-subjects dataset.

This section is divided into four parts: pre-processing procedures; channel data analysis; component data analysis; multiple subject dataset processing; statistical procedures. Procedures show a linear order list of execution. Each procedure is followed by the graph and the associated legend

2.1.1 EEGLAB: statistical software for electro-physiological data analysis

EEGLAB is a MATLAB toolbox for processing event-related and continuous electroencephalography (EEG), magnetoencephalography (MEG) and other electrophysiological data. EEGLAB provides a programming environment for accessing, visualizing, measuring, manipulating, and storing electrophysiological data. EEGLAB primary allows several modes of visualization of the single-trial and averaged data, event-related statistics, independent component analysis, time/frequency analysis and artifact rejection. EEGLAB include several plug-in such as NFT (3-D head and source location modeling), SIFT (3-D source information flow modeling), MPT (3-D source measure projection analysis), BCILAB (Brain-computer interface design & analysis), MoBILAB (Mobile brain/body imaging) and PACT (epileptic spike detection). The chief are Arnaud Delorme and Scott Makeig and now the development is possible due to Swartz Center for Computational Neuroscience (SCCN) of the Institute for Neural Computation at the University of California San Diego (UCSD) supported by the US National Institute of Neurological Disorders and Stroke (NINDS).

2.1.2 Importing channel location: information about the electrodes placement

To plot EEG scalp maps, dataset (loaded manually) must contain information about the locations of the recording electrodes. Usually all the channels use the same reference electrode for recording EEG data. Typical references (Figure 5) are one mastoid (for example, TP10 in the International 10-20 System), or vertex electrode (CZ in the International 10-20 System). There is no best common reference site. Some researches shows that non-scalp references (mastoid, earlobes, nose) introduce more noise than a scalp channel reference. If the data have been recorded with a reference, they can usually be re-referenced to any other reference channel.

illustration not visible in this excerpt

Figure 5 – Electrode placement in International 10-20 System. Red highlighted show the typical references in EEG recording (EEGLAB)

2.1.3 Filtering data to minimizing the introduction of artifacts

Filtering the continuous data minimizes the introduction of artifacts (linear trends) at level of epoch end-lines. High-pass or low-pass filter may be used (applied in different call). Another common use for band-pass filtering is to remove 50-Hz noise. The filtering step uses the linear finite impulse response (FIR), forward and backward, to ensure that delays phase introduced by the filter are nullified.

2.1.4 Extracting data epochs and removing baseline values

Procedures to study the event-related of continuously recorded EEG data:

specifying the data epochs time and the baseline value in each epoch;

extracting data epochs time to events of interest;

removing a mean baseline value from each epoch

If the objective of the analysis is to estimate changes that occur in the data following the time-locking events, the mean value in the pre-stimulus period is effective.

2.2 CHANNEL DATA ANALYSIS

2.2.1 Channel data scroll: visualization, normalization and channel rejection procedure

Data scrolling is useful to reject epochs of data which contains artifact (if the channel acquisition data not is completely compromise) before the data processing (Figure 6). Data visualization, voltage scale (microvolts), time (seconds), and rejection of the portion of continuous data were show in the following Figures.

illustration not visible in this excerpt

Figure 6 – On the left of the plot window there is a channel list ordered by name. On the right there is a vertical scale value (microvolts), which indicates the height of the vertical scale bar. In this case, that value is 0.1 microvolts. In abscissa there is the horizontal scale value (seconds), which indicates the time scale of the horizontal scale bar. In this case, that value is 5 seconds

Normalization procedure (Figure 7) of the voltage allows a better visualization of data, in this way the differences that may be highlighted between the channels are reduced.

illustration not visible in this excerpt

Figure 7 – In this case, the normalized value is 5 microvolts. The horizontal scale value (seconds), which indicates the time scale of the horizontal scale bar after normalizing procedure, because not is dependent on the activation level of the channels

Channel rejection procedure (Figure 8) allows to exclude channels with outliers, after kurtosis or probability tests.

illustration not visible in this excerpt

Figure 8 – Channel rejection by kurtosis and probability tests with threshold limits (max) set at 5%. In this case one channel (P4) was excluded by dataset used for the analysis

2.2.2 Channel spectra and associated topographical maps

Spectral analysis is one method used for EEG quantification. A signal can be analyzed with power spectrum (figure 9), which provides information on the signal power at each frequency. The EEG spectrum includes frequencies from 0.1 Hz to 100 Hz. The Fourier transform decomposes the EEG time series into power spectrum, in which the power is the square of the EEG amplitude, and the amplitude is the integral average of the EEG signal during the epoch sampled. This type of analysis is a mathematical approach to quantify the EEG. The frequency resolution is given by the inverse of the time value of the epoch. For example, with one epoch of 5 s the frequency resolution is 0.20 Hz.

illustration not visible in this excerpt

Figure 9 – Each colored trace represents the spectrum of the activity of one data channel

2.2.3 ERP and associated topographical maps

In ERP plot (Figure 10), EEG data epochs (trials) are before ordered and then smoothed with others trials, and finally color-coded. As opposed to the average ERP, which exists in only one form, the number of possible ERP plots of a set of single trials is nearly infinite. Not all of the sorting orders will give equal insights into the brain dynamics. By default, trials are sorted in the order of temporal occurrence.

illustration not visible in this excerpt

Figure 10 – Each colored trace represents the ERP of the activity of one data channel

2.2.4 Time/frequency decomposition

Time/frequency plot (Figure 11) characterizes changes in the spectral of the data, considered as a sum of sinusoidal functions. A significant ITC indicates that the EEG activity at a given time and frequency in single trials becomes phase-locked. The time/frequency points showing significant ITC and ERSP are not necessarily identical.

illustration not visible in this excerpt

Figure 11 – The top image shows mean event-related changes in spectral power from pre-stimulus (baseline =0 ms) at each time during the epoch and at each frequency (<50 Hz). The upper left panel shows the baseline mean power spectrum, and the lower part of the upper panel shows the ERSP envelope (low and high mean dB values, relative to baseline, at each time in the epoch). The lower image shows the inter-trial coherence (ITC) at all frequency.

2.2.5 Cross-coherences computation

We may plot event-related cross-coherence (Figure 12), to determine the degree of synchronization between the activations of two channels. Channels are maximally independent over the whole time range of the training data and may become transiently partially synchronized in specific frequency bands.

illustration not visible in this excerpt

Figure 12 – The upper panel shows the coherence amplitude between 0 and 1 (1 representing the perfectly synchronizing among the two signals). The lower panel indicates the phase difference among the two signals at time/frequency points where cross-coherence amplitude is significant.

illustration not visible in this excerpt

2.2.6 Channel summary

Figure 13 – Estimated variables (mean, standard deviation, skewness, kurtosis, median and Kolmogorov-Smirnov test at a significance level of P=0.05) are printed as text in the lower panel to facilitate graphical analysis and interpretation. The upper left panel shows the data histogram (blue bars), the data mean (vertical red line), a fitted normal distribution (light blue curve) and a normal distribution fitted to the trimmed data (yellow curve). The P and 1-P percentile points are marked with yellow ticks on the X axis. A horizontal notched-box plot with whiskers is showed below the histogram. The whiskers are lines extending from each edge of the box to show the extent of the rest of the data. The upper right panel shows the empirical quantile-quantile plot (QQ-plot). The QQ-plot helps to determine whether the data sample is drawn from a Normal distribution

2.2.7 Rejecting artifacts in continuous and epoch data

The approach for artifact rejection uses statistical thresholding to suggest rejectable epochs to reject from analysis. It is a semi-automated rejection coupled with visual inspection, as detailed below:

Rejecting data by visual inspection;

Rejecting data channels based on channel statistics;

Rejecting epochs by visual inspection;

Rejecting extreme values;

Rejecting abnormal trends;

Rejecting improbable data;

Rejecting abnormally distributed data;

Rejecting abnormal spectra;

Inspecting current versus previously rejections;

Inspecting results of all rejection measures.

2.3 COMPONENT DATA ANALYSIS

2.3.1 Independent Component Analysis

In the original channel data, each row of the recording matrix represents the time course of the differences between source projections to one data channel and the reference channel. After ICA, each row of the data activation matrix gives the time course of the activity of one component process spatially filtered. Independent component filters produce the maximally temporally independent signals available in the channel data. The mixing process is passive, linear, and adds no information to the data. In contrast to PCA, in which the first component may account 50% of the data, the second 25%, etc., ICA component contributions are much more homogeneous, ranging from roughly 5% down to 0%. PCA makes each successive component from the remaining activity by previously determined components, while ICA makes maximally independent sources of activity. PCA components are temporally and spatially orthogonal, while ICA components of EEG data are maximally temporally independent, but spatially unconstrained and therefore they are able to find maps representing the projection of a partially synchronized region of cortex. The ordering, scalp topography and activity time courses of best-matching components may appear slightly different due to ICA decomposition starts with a random weight matrix. Selecting specific channel types is possible to use for ICA decomposition. Studying with both EEG and EMG channels, any relationship between EEG and EMG signals should involve propagation delays and ICA assumes an instantaneous outcome of the relationship.

2.3.2 ICA Algorithms

All three algorithms return near-equivalent components:

1. Runica Infomax algorithm gives stable decompositions with up to hundreds of channels and components with a supergaussian activity distribution. The component order returned by Infomax is in decreasing order of the EEG variance accounted by each component.
2. Jader algorithm uses 4th-order moments (whereas Infomax uses a combination of higher-order moments) but the storage required for all the 4th-order moments become impractical for datasets with more than 50 channels;
3. Fastica algorithm may be less stable than Infomax for high-dimensional data sets. As a general rule, finding N stable components (from N-channel data) typically requires more than kN^2 data sample points, where N^2 is the number of weights in the unmixing matrix to learn and k is a multiplier.

2.3.3 Component data scroll

Scrolling through the ICA activations (Figure 14), easily spot components accounting for characteristic artifacts.

illustration not visible in this excerpt

Figure 14 – To the left of the plot there is a component list ordered by name. On the right there is a vertical scale value (microvolts), which indicates the height of the vertical scale bar. In this case, that value is 0.1 microvolts. In the abscissa there is the horizontal scale value (seconds), which indicates the time scale of the horizontal scale bar. In this case, that value is 5 seconds

2.3.4 Component spectra and associated topographical maps

Plotting component spectra and topographical maps (Figure 16) are interest to see which components strongly contribute. Is useful, plot the data signal (less the component activity) and estimate the power decrease in comparison to the original signal at one channel. The scale of the component activity time uses arbitrary units. However, the component scalp map values for the component activity time course is in the same unit as the data. The main criteria to determine if a component is cognitively related, a muscle artifact or some other type of artifact are: scalp maps of all component at the set frequency range (Figure 15); power spectrum at specific frequency value (Figure 16). power spectrum and associate scalp map for any singular component (Figure 17);

illustration not visible in this excerpt

Figure 15 – Scalp maps ordered by decreasing level of influence on the activation for the range of frequency considered

illustration not visible in this excerpt

Figure 16 – Scalp maps shows components contribute at 10 Hz. Each colored trace represents the spectrum of the activity of the data component. The second (from left to right) scalp map shows the scalp power distribution at the selected frequency. Other scalp maps shows the contribution of the first 5 ICA (in order of influence)

illustration not visible in this excerpt

Figure 17 – Single scalp map shows the scalp distribution of the component at 3-35 Hz. In this case IC 13 was selected

2.3.5 Time/frequency decomposition

It is more interesting to look at time-frequency (Figure 18) decompositions of component activations than separate channel activities, since independent components may directly index the activity of one brain EEG source, whereas channel activities from different parts of the brain.

illustration not visible in this excerpt

Figure 18 – The top image shows mean event-related changes in spectral power from pre-stimulus (baseline =0 ms) at each time during the epoch and at each frequency (<50 Hz). The upper left panel shows the baseline mean power spectrum, and the lower part of the upper panel shows the ERSP envelope (low and high mean dB values, relative to baseline, at each time in the epoch). The lower image shows the inter-trial coherence (ITC) at all frequency.

2.3.6 Computing cross-coherences

Determining the degree of synchronization between two components, it’s possible to plot event-related cross-coherence (Figure 19). Independent components are (maximally) independent over the whole time range of the training data, they may become synchronized in specific frequency bands.

illustration not visible in this excerpt

Figure 19 – The upper panel shows the coherence amplitude between 0 and 1 (1 representing the perfectly synchronizing among the two signals). The lower panel indicates the phase difference among the two signals at time/frequency points where cross-coherence amplitude is significant.

2.3.7 Component summary

Figure 20 – Estimated variables (mean, standard deviation, skewness, kurtosis, median and Kolmogorov-Smirnov test at a significance level of P=0.05.) are printed as text in the lower panel to facilitate graphical analysis and interpretation. The upper left panel shows the data histogram (blue bars), the data mean (vertical red line), a fitted normal distribution (light blue curve) and a normal distribution fitted to the trimmed data (yellow curve). The P and 1-P percentile points are marked with yellow ticks on the X axis. A horizontal notched-box plot with whiskers is showed below the histogram. The whiskers are lines extending from each edge of the box to show the extent of the rest of the data. The upper right panel shows the empirical quantile-quantile plot (QQ-plot). The QQ-plot helps to determine whether the data sample is drawn from a Normal distribution

2.3.8 Rejecting based on independent data components

The functions below described work when applied to data components:

Visually reject unsuitable portions of the continuous data;

Separate the data into suitable short data epochs;

Perform ICA on these epochs to extract their independent components;

Perform semi-automated and visual-inspection based rejection of data epochs on the derived components;

Visually inspect and select data epochs for rejection;

Reject the selected data epochs;

Perform ICA a second time on the pruned collection of short data epochs. This may improve the quality of the ICA decomposition;

Inspect and reject the components. Components should not be rejected before the second ICA, but after.

2.4 MULTIPLE SUBJECT DATA PROCESSING

2.4.1 Channel statistics

To compare electrophysiological results across subjects, the usual practice identifies scalp channels. The spatial relationship of any electrode site to the underlying cortical areas that generate the activities may be rather different in different subjects, depending on the physical locations, and the orientation of the cortical source areas. That is, data recorded from equivalent channel locations in different subjects may sum activity of different mixtures of underlying cortical EEG sources.

In the below example, a dataset compares two different postures was used. Analyzing each frequency band, statistically significant differences were found at the mu band in the P3 channel (Figure 21). Then, the power spectrum for P3 channel is displayed, first for the entire frequency range, 0-35 Hz (Figure 22), secondly for the mu band, 8-12 Hz (Figure 23); It is noted that statistically significant differences are between 10 and 11 Hz; At this point, we investigate if the difference among the two conditions from 10 to 11 Hz is also present in other channels (Figure 24). Finally, we separately investigate the frequency range to 10 (Figure 25) and 11 (Figure 26) Hz and highlight the channels statistically significant among the two conditions (seated and standing) for specific frequencies, in the first case channels C4, F8, T8, P3, P7 and P8, in the second case channel O1.

illustration not visible in this excerpt

Figure 21 – Scalp map shows the difference in channel activation from 8 to 12 Hz (Mu band) between two different posture conditions (seated vs standing). Left scalp map shows the channel activation during seated condition and the central scalp map shows the channel activation during standing condition. Right scalp map shows (marking in red) statistical difference among the two conditions with the level of p-value set at P<0.05. P3 (left parietal channel) shows a major activation during standing condition

illustration not visible in this excerpt

Figure 22 – Power spectrum shows the difference in P3 from 3 to 35 Hz among two different posture conditions (seated vs standing). The vertical grey line shows statistical difference among the two conditions with the level of p-value set at P<0.05. The two conditions are statistical different from 10 to 11 Hz range frequency band

illustration not visible in this excerpt

Figure 23 – Power spectrum shows the difference in P3 from Mu frequency band (3-35 Hz) among two different posture conditions (seated vs standing). The vertical grey marking shows statistical difference among the two conditions with the level of p-value set at P<0.05. The two conditions are statistical different from 10 to 11 Hz range frequency band

illustration not visible in this excerpt

Figure 24 – Channels F4, F8, T8, P3, P7, P8, O1 shows statistical differences among seated (red line) and standing (green line) conditions. The vertical grey line shows statistical difference among the two conditions with the level of p-value set at P<0.05. Standing condition shows higher power spectrum activation levels, statistically significant differences in all channel

illustration not visible in this excerpt

Figure 25 – Scalp map shows the difference in channel activation from 10 Hz (Mu band) among two different posture conditions (seated vs standing). Left scalp map shows the channel activation during seated condition and the central scalp map shows the channel activation during standing condition. Right scalp map shows (marking in red) statistical difference among the two conditions with the level of p-value set at P<0.05. F4 (middle-right frontal channel), F8 (right frontal channel), T8 (right temporal channel), P3 (middle-left parietal channel), P7 (left parietal channel) and O1 (left occipital channel) shows a major activation during standing condition

illustration not visible in this excerpt

Figure 26 – Scalp map shows the difference in channel activation from 11 Hz (which belongs to Mu band) among two different posture condition (seated vs standing). Left scalp map shows the channel activation during seated condition and the central scalp map shows the channel activation during standing condition. Right scalp map shows (marking in red) statistical difference among the two condition with the level of p-value set at P<0.05. P8 (right parietal channel), shows a major activation during orthostatic condition

2.4.2 Component statistics

This section describes how to manage and process data recorded from multiple subjects, groups, sessions, and conditions of an experimental study. It allows the definition and visualization of clusters with equivalent or similar ICA data components. Missing data channels may be replaced using spherical interpolation. Parametric or bootstrap statistics may be used with correction for multiple comparisons to compare a given measure in any experimental design. Currently, two clustering methods are available: 'kmeans' and 'neural network' clustering.

2.4.3 Clustering procedure

There are several steps in the independent component clustering procedure:

Identifying a set of EEG datasets containing ICA weights;

Specifying the subjects and group, task condition, and session for each dataset;

Identifying the component in each dataset to cluster;

Specifying and computing measures to use in clustering;

Performing component clustering using specific measures;

Viewing the scalp maps, dipole models, and activity measures of the component clusters;

Performing signal processing and statistical estimation on the clusters;

Studying the consistency and properties of the generated component clusters to determine which method of clustering produces clusters adequate for research hypothesis.

There are two types of clustering methods available. One is the PCA method and the other one is the new Measure Product (MP) method. Both methods produce reasonable clusters, but the new MP method has only one parameter and it is therefore recommended. The PCA method is older and its implementation is more stable.

2.4.4 Preparing to cluster with PCA method

We can use several measures to construct the cluster: ERP, power spectrum, ERSP, ITC, component scalp maps and their equivalent dipole model locations. The goal of the pre-clustering function is to compute an N-dimensional cluster position vector for each component. These 'cluster position' vectors will be used to measure the 'distance' of components from each other in the N-dimensional cluster space. The value of N is arbitrary but, for numeric reasons pertaining to the clustering algorithms, should be kept relatively low (<10). At this phase the normalization procedure are required, this involves dividing the measure data of all principal components by the standard deviation of the first PCA component.

2.4.5 Clustering

If anatomical component information (scalp maps and associate dipoles) are using, the experimental design does not impact in the clustering procedure. Once ICA components are clustered, it is possible to compute differences between conditions. Comparing ICA components among conditions is like comparing activities in different channels. Comparing the activities of a cluster of components between conditions could be seen as similar to comparing the activity of an individually channel for each subject. ICA components and electrode channels are both spatial filters. Each electrode channel gives the arithmetic difference between the potential reaching some scalp electrode and the potential reaching a reference electrode (or the mean of the potentials reaching the set of reference electrodes). Each ICA component gives the arithmetic weighted sum/difference of the signals reaching each of the electrodes. Another editing option is to reject “outlier” (a component whose cluster location is more than a given number of standard deviations from the location of the cluster centroid components from a cluster after clustering). The standard threshold for outliers rejection is 3 standard deviations.

2.4.6 Component clusters visualization

In computing the mean cluster scalp maps was first adjusted so as to correlate positively with the cluster mean. Then the map variances were equated. Finally, the normalized means were computed.

illustration not visible in this excerpt

Figure 27 – The top left scalp map cluster shows the power spectra activity referred to outliers distribution and the others shows the cluster for the dataset loaded

illustration not visible in this excerpt

Figure 28 – To see individual component scalp maps for components in the cluster, select the cluster of interest, for example, cluster 2 as in the figure above. Cluster 4 scalp map may differ after you have recomputed the clusters for the sample dataset. Channels missing from any of the datasets do not affect clustering or visualization of cluster scalp maps

In the below example, the same dataset compares two different postures. Analyzing each frequency band, considering all independent components no statistically significant differences were found (Figure 29). Then, the power spectrum for IC 13 is displayed, first for the entire frequency range, 0-35 Hz (Figure 30), secondly for the mu band, 8-12 Hz (Figure 31); Statistically significant differences are between 10 and 11 Hz. There isn’t always a homogeneous behavior among channel and independent component analysis, all the analysis show statistically significant differences at the mu frequency band. In conclusion, standing condition shows a greater activation levels respect to seated condition at mu band.

illustration not visible in this excerpt

Figure 29 – Power spectrum difference considering total ICs from 3 to 35 Hz among two different posture condition (seated vs standing). White bar shows no statistical difference among the two condition with the level of p-value set at P<0.05. The two conditions aren’t statistical different from 3 to 35 Hz frequency band

illustration not visible in this excerpt

Figure 30 – Power spectrum difference considering component 13 from 3 to 35 Hz among two different posture condition (seated vs standing). Black bar shows statistical difference among the two condition with the level of p-value set at P<0.05. The two conditions are statistical different from 3 to 4 Hz (Delta frequency band) and from 10 to 11 Hz (mu frequency band)

illustration not visible in this excerpt

Figure 30 – Power spectrum difference considering component 13 from 3 to 35 Hz among two different posture condition (seated vs standing). Black bar shows statistical difference among the two condition with the level of p-value set at P<0.05. The two conditions are statistical different from 3 to 4 Hz (Delta frequency band) and from 10 to 11 Hz (mu frequency band)

illustration not visible in this excerpt

Figure 31 – Power spectrum difference considering component 13 from 8 to 11 Hz among two different posture condition (seated vs standing). Black bar shows statistical difference among the two condition with the level of p-value set at P<0.05. The two conditions are statistical different from 10 to 11 Hz (mu frequency band)

2.5 STATISTICAL PROCEDURES

2.5.1 Parametric and non-parametric statistics

The same methods for statistical comparison are applied both to component clusters and to groups of data channels. It’s possible performing parametric tests (paired t-test, unpaired t-test, ANOVA) on ERPs, power spectra, ERSPs, and ITCs. The p-values are computed at every frequency, for power spectra analysis and they are computed at every time/frequency point, for ERSP and ITC time/frequency analysis.

We can also compute non-parametric statistics. Statistical test will not compute the mean condition difference, but its t-value (the mean difference divided by the standard deviation of the difference and multiplied by the square root of the number of observations, less one). The result is equivalent to using the mean difference. The advantage is that when there are more conditions, the ANOVA measure can be performed. In EEGLAB, computing the probability density distribution of the t-test or ANOVA is only a "expedient" to be able to obtain a difference measure across all subjects and conditions.

2.5.2 Paired/unpaired samples

There is no direct correspondence between values for unpaired data samples. In contrast, for paired samples, each value in one sample corresponds to a value in the other sample. Note that paired groups must necessarily be the same size. Matched/unmatched data samples are an extension of paired/unpaired data samples when there are more than two samples.

2.5.3 Re-sampling methods

Resampling methods help provide confidence intervals for parameters in situations where these are hard or impossible to derive analytically. Resampling methods also help to perform statistical inference without assuming a known probability distribution for the data.

The bootstrap method consists of drawing random sub-samples and the randomization method consists of shuffling data samples; It’s the most recently developed method to estimate errors and other statistics.

Moreover, bootstrap statistics use the implicit assumption that data samples are representative of the underlying population. Therefore it is not possible to draw direct conclusions about the underlying population. If it lies in the lower 2.5% or upper 2.5% tails, then the bootstrap test may be considered to be significant at the 5% level of significance.

2.5.4 Multivariate methods (PCA vs ICA)

Many multivariate analytical methods involve inference for the parameters (means, variances, and correlation coefficients) based on multivariate normal distribution. Principal component analysis (PCA) would specifically makes each successive component account for as much as possible of the remaining variability uncorrelated with previously determined components.

Recent progresses in signal processing and information theory have seen the development of blind source separation methods, which attempt to find a coordinate frame onto which the data projections have minimal overlap.

Independent component analysis (ICA) is a family of linear blind source separation methods. The core mathematical concept of ICA is to minimize the mutual information among the data projections.

PCA components are orthogonal, usually aren’t a realistic assumption for bio-physical data. To find biologically plausible sources, PCA must be followed by an axis rotation procedure, and ICA can be viewed as a powerful rotation method. ICA seeks to find axes for which the projection of data is maximally non-normal (i.e., contains the maximum amount of information). ICA is free to adapt to the actual projection patterns of source generators, if their activity time courses are (near) independent of one another. Performing ICA decompositions is most appropriate when sources are linearly mixed in the recorded signal, without differential time delays.

ICA is being applied to various biomedical signal processing problems that include performing speech and noise separation, decomposing functional resonance imaging data, and separating brain area activities and artifacts mixed in electro-encephalographic scalp sensors.

2.5.5 Correcting for multiple comparisons

When performing a large number of statistical inferences, it is necessary to correct for multiple comparisons. For example, with a statistical threshold at p<0.05, by definition about 5% of the inferred significant values will be false positives:

Bonferroni: The most conservative method, simply divides the p-value by the number of comparisons. For example, when computing ERSP time-frequency images of 100 frequencies by 200 time points, the number of inferences is 20,000. To correct for multiple comparison at p<0.05, a statistical threshold of 0.05/20000 = 0.0000025 should be applied. This method is quite conservative as, essentially, it assumes (erroneously) that each of the time/frequency point values is independent of the others.

Holms method: The Holms method, also called Holms-Bonferroni method is not as conservative. Actual uncorrected p-values are sorted and, to assess whether a given p-value reaches the corrected threshold for multiple comparison, the lowest uncorrected p-value is compared to the Bonferroni statistical threshold of 0.05/20000. The next lowest is compared to statistical threshold of 0.05/(20000-1). The highest uncorrected p-value is compared to the uncorrected threshold 0.05/(20000-19999)=0.05.

False Discovery Rate: The False Discovery Rate (FDR) method is based on Holms correction. Under FDR, a common corrected threshold is applied to all p-values. This threshold is the last significant threshold calculated using Holms correction. For example, if the 423th strongest p-value is the smallest (and thus the last one in Holms correction) that is significant, this p-value is used as the corrected threshold for all uncorrected p-values. Note that under Holms correction, statistical assessment is performed independently for each p-value. Therefore, after sorting by p-value the 48th p-value might not be significant while the 49th might be.

Max method: the max method is only available when using non-parametric test. At each iteration in computing a surrogate distribution, the maximum statistics across all time-frequency points is considered.

Cluster method: The cluster method is also only available when using non-parametric statistics.

3. EXPERIMENTAL STUDY DESIGNS

The aim of this section is to show the analysis procedures discussed in the previous section applied on research studies conducted to the Motor Science Research Center of Turin:

3.1 ELECTROCORTICAL ACTIVITY IN DANCERS AND NON-DANCERS LISTENING TO DIFFERENT MUSIC GENRE AND DURING IMAGINATIVE DANCE MOTOR ACTIVITY

3.1.1 Abstract

Introdution: The protocol was designed to highlight differences in electrocortical activity between a group of dancers and a group of controls completely unrelated to dance, listening to three different kinds of music: classical, rock, and waltz musics. Two were the experimental conditions: the pure listening and the listening associated to the imagination of an impromptu dance related to the music.

Materials and methods: Twenty one volunteers, 11 belonging to the expert group and 10 to the control one took part in the study; The EEG was recorded using a prototype of a wireless amplifier. For statistical analysis the ANOVA test to intra-group analysis and independent t-test for inter-group analysis was performed.

Results: No statistically significant differences were found between groups at baseline. The analysis of the scalp maps differences (considering channel data) among the three listening conditions (open eyes, closed eyes and imagination) showed statistically significant differences only for the dancers group at delta (Fp1, T3 and O2) and theta (Fp1) frequency bands. Statistically significant differences (p<0,05) emerges in EEG power spectrum, comparing the two groups, in particular during the imaginative dancing motor activity: delta, theta and alpha bands listening to classical music; all bands listening to rock music; delta, beta and gamma bands listening to waltz music. Alpha band power spectrum shows statistically significant differences during imaginative dancing motor activity listening classical and rock music.

Discussions and conclusions: The statistically significant differences between expert dancers and controls during imaginative dancing motor activity could indicate a difference in the attentional effort that could be explained in a different brain activity during the development of motor programs related to the dance. The remarks partially confirm existing findings on the relationship among EEG activity, creative thinking and imaginative motor activity. Findings also extend the knowledge on the electrocortical response to auditory stimuli and during imaginative motor activity, emphasizing the difference between experienced and inexperienced subjects, relative to the field analyzed.

3.1.2 Introduction

Dance is the way through which the dancer expresses her/his being intimately. Dancing is an action involving creativity, technique and skill. Thanks to the imagination, the knowledge gained with study and experience are shaped to live new forms of expression. The imagination is the most powerful tool that human-beings have to go beyond the limits of the achievable. Creativity and imagination are the foundation of what in dance is called "improvisation", a spontaneous representation of the inner world. Dance is movement and rhythm at the same time then and, besides being an artistic expression, it is also and especially motor activity intimately tied to the music. The most interesting neurophysiological aspects relating to dance are the electro-cortical effects induced listening to music, the effect of imagination resulting in CNS (central nervous system) during the improvisation of dance, and the differences between experts and non-experts.

Analyzing the literature it is possible to highlight a number of studies investigating the electroencephalographic activity listening music. Fink, Graif and Neubauer, in 2009, investigated creativity through EEG in 15 professional dancers who have attained a high level of expertise in this domain. The group was compared with 17 novices who have only basic experience in dancing and completed no comprehensive training in this field. The EEG was recorded during performance of two different dancing imagery tasks which differed with respect to creative demands. In the first task participants were instructed to mentally perform a dance which should be as unique and original as possible (improvisation dance). In the waltz task they were asked to imagine dancing the waltz, a standard dance which involves a sequence of monotonous steps (lower creative demands). They observed that professional dancers show different alpha activity in posterior parietal brain regions than novice dancers: during improvisation dance, professional dancers exhibited different alpha activity in right-hemispheric than the group of novices; during imagining dancing the waltz no significant group differences emerged.

Bhattacharya, Petsche, Pereda in 2001 analyzed the electrical activity of different brain regions during the performance of higher cognitive functions. Their goal was to check both task-related differences (e.g. listening to music vs. rest) and training-related differences (musicians vs. non-musicians). The results demonstrate the occurrence of task-related differences in both groups of subjects. Furthermore, subjects with musical training possessed significantly differences in electro-cortical activity than such without musical training while listening to music but not to text.

In a further study (Bruyna and Severtesen, 1984) a greater activation of the beta rhythm while listening to the songs of Chopin or Halpem with open eyes and an increase in the alpha rhythm while listening with closed eyes were observed.

Steinberg, Gunther and Stiltz in 1992 investigated four independent studies of EEG. Findings during music perception were compared and discussed with respect to common results and to current psycho-physiological knowledge about music perception. The influence of music training and gender of the subjects appeared to be important as did the nature of the music stimulation. With musically non-educated subjects, listening to monotonously structured, simple sounds, an alpha power augmentation over the whole skull and a beta increase resulted, involving consistently primary acoustical areas. Theta diminution was seen over left rolandic and occipito-parietal sides of both hemispheres. Perception of classical and light music led to a considerable decrease of theta and alpha power, primarily in left fronto-temporal electrodes, whereas beta decrease in pre-frontal areas seemed to be restricted to males. Also, an increase of higher beta frequency over temporal and parietal-occipital parts of the skull was revealed. Musically trained subjects were examined separately by one group of investigators. They showed more pronounced results, primarily in an ample left-sided alpha suppression when compared to resting conditions. Retesting after several weeks was performed by two research groups, leading to a discussion of the importance of directing attention to careful experimental conditions.

Jaušovec and Habe in 2003 studied the event-related responses of 18 individuals while they were listening to 3 music clips of 6 seconds duration. The music clips differed in the level of their complex structure. They were taken from Mozart's sonata (K. 448), and Brahms' Hungarian dance (no. 5). The third clip was a simplified version of the theme taken from Haydn's symphony (no. 94) played by a computer synthesizer. Significant differences among the 3 music clips were only observed in the lower-1 alpha band which is related to attentional processes. While respondents listened to the Mozart clip, in the lower alpha bands increased more, whereas in the gamma band a less pronounced increase was observed as compared with the Brahms and Haydn clips. The clustering of the three clips based on EEG measures distinguished among the Mozart clip on the one hand, and the Haydn and Brahms clips on the other, even though the Haydn and Brahms clips were at the opposite extremes with regard to the mood they induced in listeners, musical time, and complexity of structure. This would suggest that Mozart's music influences the level of arousal.

[...]