An exploratory data analysis of electroencephalograms using the functional boxplots approach

Many model-based methods have been developed over the last several decades for analysis of electroencephalograms (EEGs) in order to understand electrical neural data. In this work, we propose to use the functional boxplot (FBP) to analyze log periodograms of EEG time series data in the spectral domain. The functional bloxplot approach produces a median curve—which is not equivalent to connecting medians obtained from frequency-specific boxplots. In addition, this approach identifies a functional median, summarizes variability, and detects potential outliers. By extending FBPs analysis from one-dimensional curves to surfaces, surface boxplots are also used to explore the variation of the spectral power for the alpha (8–12 Hz) and beta (16–32 Hz) frequency bands across the brain cortical surface. By using rank-based nonparametric tests, we also investigate the stationarity of EEG traces across an exam acquired during resting-state by comparing the spectrum during the early vs. late phases of a single resting-state EEG exam.


Introduction
Electroencephalograms (EEGs) have been used for many decades to study the complex spatiotemporal dynamics of brain processes (Nunez and Srinivasan, 2006). Due to its excellent temporal resolution (sampling rates usually range from 100 to 1000 Hz), EEGs can capture transient changes in brain activity, identify oscillatory behavior and study cross-dependence between EEG components. Since EEGs indirectly measure neuronal electrical activity, they can be used to infer the statistical properties of the underlying brain stochastic process. One such statistical property is the spectrum (or power spectrum) which decomposes the total variability in the EEG according to the contribution of oscillations at different frequencies. Most approaches to analyzing EEGs focus immediately on statistical modeling and spectral estimation. Here, we offer a systematic framework for exploring structures, patterns and features in the signal-prior to formal modeling. We explore the spectral properties only in a single channel using EEG traces from several epochs.
One approach to estimating the spectrum using EEG traces is to fit a parametric time domain model, such as the autoregressive moving average (ARMA) model. Applications of parametric modeling of EEGs have a long history. See (Bohlin, 1973;Isaksson et al., 1981;Krystal et al., 1999;Jain and Deshpande, 2004) among many others. When the spectrum of the EEG evolves over time (e.g., within an epoch), one could still use the ARMA model but allow the coefficients to vary over time. A key element in ARMA models is the order of the autoregressive (AR) and moving average (MA) components. These can be obtained objectively using an information-theoretic criterion such as the Akaike information criterion (AIC) and the Bayesian information criterion (BIC). Using these criteria, we obtain an optimal AR and MA order that jointly gives the best fit with the least complexity (as determined by the orders). BIC puts a heavier penalty for complexity compared to AIC and thus often gives a model with lower orders (lower complexity). From the parametric fit, we derive the estimates of the auto-correlation function and the spectrum. The theoretical background for parametric models are developed in Priestley (1981), Shumway and Stoffer (2000), and Brockwell and Davis (2009).
One could also estimate the spectrum without resorting to a parametric model. Under this approach, the EEGs are considered to be superpositions of sines and cosines (Fourier waveforms) with different frequencies and random amplitudes. These random amplitudes (or coefficients) are computed using the fast Fourier transform (FFT). The squared magnitude of these amplitudes, often called the periodograms, are the dataanalogs of the spectrum defined on discrete frequencies. The theoretical background on the frequency domain approach to time series is developed in Brillinger (1981) and Percival and Walden (1993). This approach to analyzing EEGs continues to be popular in the cognitive and brain sciences. The following papers cover both methods and applications of spectral analysis to EEGs: Pfurtscheller and Aranibar (1979), Bressler and Freeman (1980), Makeig (1993), and Srinivasan and Deng (2012), to name a few.
The common practice prior to spectral estimation is to preprocess EEGs, often to remove artifacts (Makeig et al., 1996). After artifact rejection and segmentation according to epochs, the spectrum is estimated from each EEG trace. As noted, there is a lack a systematic framework for exploring structures, patterns and features in the signal-prior to formal modeling. Due to the complexity of EEG data, exploratory data analysis (EDA) plays an important role, especially when data are recorded from many epochs or trials during an experiment. For example, it is often expected that brain responses to the same stimulus ought to be relatively uniform, with minimal variation across epochs. In contrast, greater variability across epochs may be expected during neuroimaging studies that examine the brain in restingstate, as cognitive processes can vary within and across sessions for individual subjects and across subjects. An appropriate EDA methods can provide insights into features of EEG, including similarities and variability of the brain responses across epochs to facilitate the statistical model. In this paper, we propose to use the functional boxplot (FBP) method originally developed by Sun and Genton (2011) to address these questions.
The methods presented in this paper are motivated by a motor skill acquisition study at the Neuro-rehabilitation laboratory at the University of California, Irvine (Principal Investigator: Steven C. Cramer). In the previous study, EEG was recorded from 17 subjects both during resting-state prior to motor skill training and during motor skill training using dense-array EEG (256 electrodes) as shown in Figure 1. The resting-state EEG exam was 3 min, and during post-processing, was segmented into 1-s non-overlapping epochs. As demonstrated in Wu et al. (2014), the spectral features of the resting-state EEGs when combined with a partial least squares regression analysis, was predictive of an individual's subsequent ability to acquire a novel motor skill. These may be of clinical importance to the field of rehabilitation, as improved methods for stratifying patients may significantly improve response to treatment and assist allotment of limited resources.
We present an exploratory spectral analysis (ESA) of restingstate EEG traces using FBPs for one subject. In spectral analysis, the spectrum is an important stochastic property of the signal. It indicates the amount (or proportion) of variance that is explained by each frequency bin. Thus, the spectrum or the log spectrum of the EEG signal can be used to examine relative amounts of variability explained by slow (delta or theta) waves and fast (alpha or beta) waves. Throughout this analysis, we obtain a sample spectral curve by smoothing the log periodograms of each 1-s EEG epoch, and treat it as one observation unit in the FBP. By using the FBP, we address three primary objectives. The first objective is to identify the median, i.e., the most characteristic spectral curve rather than the pointwise frequencyspecific medians. In addition, outliers are demonstrated by their unusual sample log spectral curve, and can be caused by extrabrain artifacts, including eye blinks, eye movements, and muscle movements in the EEG signal. Subsequently, confirmed outliers will be removed from subsequent analyses. The advantage of the FBP approach, over the usual pointwise boxplot method, is that it identifies epochs that have potential outlying spectral curves.
The second objective is to compare the median curves and the variability of the spectral curves from multiple phases of the resting state period. To test the stationarity of the EEG signal over the entire recording, we compare the spectral curves and the frequency-specific spatial distribution of spectral power during the early phase (first 60 epochs) vs. the late phase (last 60 epochs). Evidence against stationarity must be taken seriously since this would suggest an evolution of brain processes across the recording (Fiecas and Ombao, under revision). Moreover, the FBP approach is able to provide some characterization of the variation of the sample log spectral curves across EEG recording. In experiments comparing more than one group (e.g., healthy controls vs. patients with stroke), it would be also interesting to determine whether groups differ with respect to consistency (uniformity) of the EEG signal over time.
The third objective is to investigate the spatial variability of spectral power across the brain for a given frequency band using the surface boxplot, which is a generalization of the FBP. Using the surface boxplots approach, it is possible to identify cortical regions (or channels) that, relative to the other channels, exhibit a high proportion of beta power. The beta band is particular interest to neuroscientists, as changes in beta activity have a good association with motor function (Roopun et al., 2006;Joundi et al., 2012).
The remainder of the paper is organized as follows. In Section 2, we present a comprehensive exploratory method which consists of the following: a review of the spectra in Section 2.1, a demonstration of automatic bandwidth selector for periodogram smoothing using the gamma generalized crossvalidation criterion in Section 2.2, some remarks on smoothing the periodogram in Section 2.3, a description of the FBPs in Section 2.4, a description of the surface boxplots in Section 2.5, and a demonstration of testing for differences in mean curves between families of curves in Section 2.6. In Section 3, we examine the finite sample performance of the proposed exploratory method. In Section 4, the resting-state EEG data are analyzed. Finally, in Section 5, conclusions and future work are discussed.

Method for Exploratory Spectral Analysis (ESA)
In this section, we review the methods that are needed for ESA of the EEG data. In Section 2.1, we first formally define the spectrum and then discuss a consistent estimator which is obtained by smoothing the periodogram using a bandwidth that is automatically selected by the gamma generalized crossvalidation (Gamma-GCV) method described in Section 2.2. Next, we highlight two remarks on smoothing the periodogram in Section 2.3, then we present the FBPs method in Section 2.4 and surface boxplots method in Section 2.5. Finally, we present a rank sum test which tests for differences in median curves or surfaces between families of curves or surfaces in Section 2.6.

Spectrum
The spectrum of an EEG signal (which is assumed to be stationary) can give the amount of variance contributed by oscillatory components (from delta to beta band activity). Let X(t), t = . . . , −1, 0, 1, . . . be a zero-mean stationary time series with covariance function γ (τ ) = E X(t)X(t + τ ) (τ = . . . − 1, 0, 1, . . .) that is assumed to be absolutely summable, i.e., The starting point for estimating f (ω) is the periodogram. Denote I(ω k ) to be the periodogram computed from a finite sample of the stationary process X(0), X(2), . . . , X(T − 1) at frequency ω k = k/T which is defined to be where [[T/2]] is the quotient of T/2. To characterize the spectra of the EEG signals, we classify the oscillatory patterns of periodograms into four primary frequency bands: delta (0-4 Hz), theta (4-8 Hz), alpha (8-16 Hz), beta (16-32 Hz), and gamma (32-50 Hz) as shown in Figure 2. Since each frequency band is defined by a range, we define S( ) to be the estimated spectral power at the band: It is well-known that the periodogram I(ω k ) is an asymptotically unbiased estimator for f (ω k ), but it is inconsistent because its variance approaches a positive constant when T → ∞. Therefore, to reduce the variance, we smoothed the periodogram. A number of nonparametric smoothing methods have been proposed including the kernel smoother (Lee, 1997;Ombao et al., 2001), wavelet (Gao, 1997), smoothing spline (Wahba, 1980;Pawitan and O'sullivan, 1994), or local polynomial (Fan and Kreutzberger, 1998). For kernel smoothing, Ombao et al. (2001) developed an automatic span selector via the generalized crossvalidation criterion for generalized additive models based on the deviance which is discussed in Section 2.2.

Automatic Span Selector Using the Gamma Generalized Crossvalidation Method
From Brillinger (1981) (Theorem 5.2.6), I(ω k ) follows an asymptotic distribution where I(ω 0 ), . . . , I(ω T/2 ) are independent. As a caveat, we note here that the actual result requires that the number of frequencies is fixed and does not depend on T. However, in most applications, this is often ignored. This result can be equivalently stated as I(ω k )/f (ω k ) ∼ ǫ k where ǫ k∼ χ 2 (1) when k = 0 or T/2 and ǫ k∼ 1 2 χ 2 (2) when k = 1, . . . , T/2 − 1. As noted, we need to smooth the periodogram I(ω k ) to produce a consistent estimator for f (ω k ). Let f p (ω k ) be a smoothed periodogram estimator of f (ω k ) which we define to be f p (ω k ) = p j = −p W p,j I(ω k + j ) k = 0, . . . , T/2, and j = −p, . . . , p where 2p + 1 is the smoothing span and W p,j are nonnegative weights that satisfy the following conditions for any fixed p: W p,j = W p,−j (j = 1, . . . , p), p j=−p W p,j = 1.
Frontiers in Neuroscience | www.frontiersin.org Generally, the weights are chosen so that W p,j is a decreasing function of p, but (Priestley, 1981) shows that the choice of the weights W p,j is of secondary importance to the value of the span or bandwidth. Thus, for simplicity, we use the boxcar smoother with weights defined by W p,j = 1/(2p + 1) for all j = −p, . . . , p. The gamma generalized crossvalidation method selects p to minimize the generalized crossvalidated deviance function and Nelder, 1989). Here, q j = 1 − 0.5I{j = 0, M − 1}, and I is the indicator function. The H p is the smoother matrix with smoothing parameter p, and the term (1 − tr(H p )/M) 2 often referred to as the model degrees of freedom, can be expressed in terms of the weight at the center of the smoothing window: (1 − W p,0 ) 2 . Then, the generalized crossvalidated deviance function can be written as

Remarks
For frequencies over 100 Hz, the periodogram values are almost negligible because the signals underwent low-pass filtering at 100 Hz. , so for simplicity, we will only show the spectrum over the frequency range of 0-100 Hz. In Figure 3, we show the location of channel 197 in right pre-motor region at the resting-state. Figure 4 gives an illustration of smoothing the periodograms for randomly selected epochs 3, 85, and 160 for a fixed channel 197. It can be seen that the power at these periodograms are dominated by low frequencies, and the values of smoothing span minimizing the generalized crossvalidated deviance function are about 3-5. Also, the smoothing lines reasonably approximate the periodograms and the small bandwidths preserve the peaks. Second, since the distribution of I(ω k ) is a multiple of the spectral density, its variance [which depends on f (ω k )] also changes across the frequencies ω k . To stabilize the variance across frequencies and to standardize comparisons of median curves across two phases (early vs. late phases of the resting-state EEG recording) we will use the log transformed periodograms. It is convenient then, that the variance of the log periodograms at each frequency is constant and takes the approximate value of π 2 6 . Moreover, while the periodogram is approximately unbiased for the spectrum, the log periodogram is no longer (approximately) unbiased for the log spectrum due to Jensen's inequality. This is easily fixed by adding the Euler Mascheroni constant 0.57721 to log transformed periodograms to obtain the log bias-corrected periodograms (Wahba, 1980). Let g(ω k ) be the true log spectrum, then Y r (ω k ), the log bias of the corrected periodogram at epoch r, is defined as Y r (ω k ) = g(ω k ) + 0.57721, k = 0, 1, . . . , T/2. Figure 5 gives the log bias-corrected periodograms, Y r (ω k ), corresponding to Figure 4. Throughout this paper, we will apply the gamma crossvalidation method to obtain the optimal smoother of log bias-corrected periodograms.

Functional Boxplots
The FBP is constructed in a similar manner to the classical (pointwise) boxplot. Each observation will be sorted based on decreasing values of some depth measure, and band depth is one notion. A curve is said to be "deeply situated" within a sample of curves if it is covered by many bands from pairs of curves. This idea is an extension of a pointwise boxplot where the median is also located "deep" in a sample because it is situated in the middle of the boxplot and hence covered by many pairs of points. Here, our observation units are curves (or real-valued functions) which are the log bias-corrected periodograms Y r (ω k ), k = 0, . . . , T/2 over many epochs r. The notion of a band depth was introduced in López-Pintado and Romo (2009) through a graph-based approach to order all sample curves which we briefly describe. Suppose that a curve The band in R 2 can be delimited by a number J of curves, and this number is fixed as J = 2 in our study. Now, let Y α , Y β be two continuous functions, Let Y 1 , . . . , Y n be n independent sample curves, then the band depth for a given curve Y i , i = 1, . . . , n is defined as where I(·) is the indicator function. When J = 2, there are n 2 possible bands delimited by two curves. The limit of the band depth BD is that it does not measure the proportion of curve inside the band. Thus, López-Pintado and Romo (2009) also proposed a modified band depth method (MBD), which measures the proportion of a curve Y i that is actually in a band: , and λ is a Lebesgue measure on A. We notice that the MBD computation will be time-consuming when n is large, so we use an exact fast method from Sun et al. (2012) to compute the MBD for the EEG data. Based on the ranks of the depths of the curves, the FBPs can provide the descriptive statistics, such as the 50% central region, the median curve, and the maximum and minimum nonoutlying curves. Moreover, the potential outliers can be detected by the 1.5 times inter-quartile range (IQR) empirical rule, which is commonly used for classical boxplots. The boundary region is defined as 1.5 times the height of the 50% central region. Any curves outside this region are considered potential outliers. In contrast with a constant factor 1.5 in classical boxplot, a factor 1.5 in FBP can be modified due to potential spatio-temporal outliers. This is because the curves from different locations will be spatially correlated, and there can be dependence in time/frequency for each curve (Sun and Genton, 2012a).

Surface Boxplots
Similar to FBPs, one can compute the data depth of all the observations, then order them according to decreasing depth values. Suppose that the observed sample surfaces, z 1 (s), . . . , z n (s), s ∈ S, where S is a region in R 2 . The information unit for such a dataset is the entire surface. To order sample surfaces, we need to generalize univariate order statistics to surfaces. To this end, we generalize the MBD with J = 2 to R 3 through a volume. Genton et al. (2014) define the sample modified volume depth (MVD) to be , if λ is the Lebesgue measure on R 3 . A sample median surface is a surface from the sample with the largest sample MVD value, designed by arg max z∈z 1 ,...,z n MVD n (z). If there are ties, the median will be the average of the surfaces maximizing the sample MVD.
The first step for constructing surface boxplots is the surface ordering. Sample surfaces are ordered from the center outwards based on their MVD values, inducing the order z [1] , z [2] , . . . , z [n] . The sample α central region is naturally defined as the volume delimited by the α proportion (0 < α < 1) of the deepest surfaces. In particular, the sample 50% central region is where [n/2] is the smallest integer not less than n/2. The border of the 50% central region is defined as the inner envelope representing the box in a surface boxplot. This is the surface analog of the first and third quartiles of the classical boxplot. The median surface in the box is the one with the largest depth value. Because the ordering is from the center outwards, the volume of the central region increases as α increases. Hence, the maximum envelope, or the outer envelope, is defined as the border of the maximum non-outlying central region. To determine this region, we propose to identify outlying surfaces by an empirical rule similar to the 1.5 times the 50% central region rule in a FBP. The fences (or the upper and lower surface boundaries for flagging potential outliers) are obtained by inflating the inner envelope (as defined above) by 1.5 times the height of the 50% central region. Any surface crossing the fences are flagged as potential outliers. The factor 1.5 can be also adjusted as in the adjusted FBPs to take into account spatial autocorrelation and possible correlations between surfaces.

Testing for Differences in Median Between Families of Curves or Surfaces
To compare the median curves from two populations of curves, López-Pintado and Romo (2009) proposed the rank sum test. Let µ Y andμ Y ′ be the median curves of two populations Y and Y ′ , respectively. Define the null hypothesis to be Suppose that we observe two sets of curves, namely {y 1 , . . . , y n } and {y ′ 1 , . . . , y ′ m }. Then define the reference sample to be {r 1 , . . . , r k } which is from one of the two observed sets with k ≥ max(n, m). The position of a particular y i for i = 1, . . . , n, or y ′ j for j = 1, . . . , m with respect to the reference sample r, is defined as where MBD is the MBD defined in previous section, and I is the indicator. Then, we can order the values R(y i ) and R(y ′ i ) from the smallest to the largest, and their ranks are between 1 and n+m. The test statistics T = m l=1 rank R(y ′ j ), then under the null hypothesis H 0 , the distribution of T is the distribution of the sum of m numbers that are randomly chosen from 1, 2, . . . , n + m (Sun and Genton, 2012b).

Remarks on the Applications of Functional and Surface Boxplots
In this paper, we use functional and surface boxplots to explore the structure of EEGs. However, these methods are general and can be applied to other types of data such as growth data and climate time series (Sun and Genton, 2012b).

Simulation Study
The purpose of the simulation study is to examine the performance of the exploratory spectral methods under various experimental settings. In Section 3.1, we demonstrate the performance of the FBP on the smoothed log periodograms of a mixture of two first order AR time series, denoted AR(1). In Section 3.2, we illustrate the rank sum test to compare the functional median from two families of curves.

Functional Boxplot Simulation Study
For the r th epoch, let U 1r (t) be an AR(1) process with its spectra dominated by high frequencies and U 2r (t) be another AR(1) with its spectra mostly containing low frequencies. The AR(1) parameters are allowed to vary across epochs. Here, we set t ∈ T = {1, . . . , 1000}. We define X r (t) to be the mixture of U 1r (t) and U 2r (t), such that X r (t) = a 1r U 1r (t) + a 2r U 2r (t) where r = 1, . . . , 220, a 1r and a 2r are weighted coefficients of U 1r (t) and U 2r (t), respectively. Then, the model for high and low frequency AR(1) processes are defined as U ℓr (t) = φ ℓr U ℓr (t − 1) + W rt where ℓ = 1, 2 and W(t) is white noise. In this setting, the high and low frequency AR(1) are distinguished by the value of φ ℓr . For example, for high frequency U 1r (t), we set φ 1r = 0.9 + ξ r , where ξ r are independent and identically distributed from N (0, 0.001). Similarly, for low frequency U 2r (t), we set φ 2r = −0.5 + η r , and η r are also independent and identically distributed from N (0, 0.001). Here, we need the variance of ξ r and η r to be small so that it guarantees causality, i.e., ξ r ∈ (−1, 1) and η r ∈ (−1, 1). Next, we split the 220 subjects into two groups, such that the first group will include both high and low frequency series, U 1r (t) and U 2r (t), while the second group will only have the high frequency series U 1r (t). To split X r (t) into two groups, we set the weight coefficients a 1r and a 2r as following a 1r ∼ N (10, 1) for r = 1, . . . , 220 a 2r ∼ N (5, 1), for r = 1, . . . , 120, and a 2r ∼ N (0, 0.001) for r = 121, . . . , 220.
The two groups of X r (t) are shown in Figure 6. Using the gamma generalized crossvalidation method, Figure 7 displays the log bias-corrected periodograms for each group, and Figure 8 shows the corresponding FBPs. Note that group 1 is dominated by both high (right) and low (left) frequencies while group 2 includes only low frequencies. Thus, the functional median of group 1 should have two peaks, one each in high and low frequency ranges, while the functional median of group 2 has only one peak in the low frequency range. In Figure 8, the black curve is the median curve in the center of the FBP. The two median curves from each group have clearly summarized the typical power distribution for each group. The blue curves in the center form the envelope of the 50% central region. The blue curves outside of the 50% central region are the non-outlying minimum and maximum curves. It is worth remarking that the envelope of group 1 is smaller than the envelope of group 2, and therefore, we demonstrate that group 2 has more dispersion than group 1. Moreover, the envelope of group 1 is in the middle of the non-outlying minimum and maximum curves, while the envelope of group 2 tends to move upwards. This indicates that group 2 shows more skewness than group 1. The red dashed curve in Figure 8 denotes the outliers. We see that the curves from group 1 that are dominated by high frequencies only are detected as outliers while the curves from group 2 that include both high and low frequencies are detected as outliers.
In order to illustrate the usefulness of the FBP compared to the pointwise boxplot, we introduce a simulation study which randomly chooses 10 bias-corrected log periodograms among 160 total periodograms. We simulate an outlying curve by adding additional noise across the 0-100 Hz frequency range, and close to the center for the remaining frequencies. Figure 9A shows the simulation data including the 10 random bias-corrected log periodograms (gray curves) and a simulated outlying curve (red curve). In Figure 9B, the FBP successfully detects the simulated outlying curve and other outliers. However, Figure 9C shows that the pointwise boxplot fails to detect the simulated outlying curve, and provides some disconnected outlying curves across frequencies. We also notice that the non-outlying maximum and minimum curves of pointwise boxplot are actually the outlying curves detected by FBP. Figure 9D compares the two median curves from these two methods, and by visual inspection, there is a slight difference between the two median curves at low frequencies. Thus, FBP can be a non-parametric method to obtain the median curve and the variability around it for EEG data compared to pointwise boxplot.

Rank Sum Test Simulation Study
To investigate the performance of this nonparametric test, we simulated two sets of curves, which are defined as below: where r = 50, ℓ = 1, 2, g(ω k ) = 1 for all ω k , and ω k is defined as ω k = k/100, where k = 1, . . . , 100. In the model, curves, respectively. Let the function f 1 (ω k ) be defined as f 1 (ω k ) = 5 · 1000 · ω k , and consider three different cases: 1. The two means are identical, let f 2 (ω k ) = f 1 (ω k ) for all ω k .
We applied the kernel average smoother with window size 7 to smooth each curve from these two families. Figure 10 illustrates the simulated curves (left panel) and the smoothed curves (right panel). In order to investigate the rank sum test performance in each case, we simulated two families of curves and obtain pvalues of rank sum test; this procedure was repeated 1000 times. Let the type I error α be 5%, we report the percentage of time that the rank sum test rejects H o : f 1 (ω k ) = f 2 (ω k ) for all ω k in Table 1.
Overall, the rank sum test method performed well in each case. When the two families are identical, this method rejected the null hypothesis of equality only 44 times (4.4%) out of 1000 times, which is close to the nominal α. When the two families are nearly identical, this method rejects 605 times (the power is 60.5%), and when the two families are completely different, the power is 100%. Thus, this method demonstrates power and sensitivity to differences.

Data Description
In this paper, we analyze EEG data from one participant in a resting-state EEG study approved by the Institutional Review Board of the University of California, Irvine. The overarching aim of this study was to identify a pattern of EEGderived coherence acquired during rest-state that could predict subsequent response to training on a novel motor skill. During EEG acquisition, subjects sat quietly with both feet flat on the floor, and were instructed to fixate their gaze to the center of a fixation cross. Each recording was 3 min in duration. While the original EEG recording included 256 channels, only 194 were used in subsequent analyses, as extra-brain artifacts, including cheek and neck muscle artifacts, and heart rhythms, are more likely to contaminate EEG signals recorded from electrodes overlying cheek and neck regions. Following data acquisition, pre-processing steps included: 100 Hz low pass filter; EEG segmentation into 1-s consecutive, non-overlapping epochs; mean detrend; and EEG signal re-reference to mean signal across all 194 channels. In addition, a combination of visual inspection and Infomax Independent Component Analysis decomposition were used to remove extra-brain artifacts, including eye blinks, eye movements, muscle artifact, and heart rhythm artifacts. The final dataset consisted of 160 epochs, with each epoch lasting 1 s, and T = 1000 time points for each epoch.
The goals of the present analysis are as follows: In Section 4.2, we closely examined a representative channel in the pre-motor region (specifically channel 197 in this dataset). Since EEGs are not well-localized in space (as opposed to local field potentials), conclusions are constrained to the sensor space. However, electrical activity captured in channel 197 reflects activity roughly around the pre-motor area. Specifically, we estimated the (log) spectrum for each epoch to identify any frequency bin or frequency band that accounts for the majority of the power spectrum. Moreover, using the method of estimating the functional medians, we obtained an estimate of the median curve from the log periodogram curves obtained from several epochs. The median curve is interpreted as a "typical" (log) spectral profile across several epochs. Using this method, we also identified outlier curves which could also be interpreted as epochs with "unusual" EEG activity. In Section 4.3, we investigated the possibility of non-stationarity across the 3 min resting-state EEG recording. Our specific goal was to compare the log spectrum during the early phase (first 60 A B C D E F FIGURE 10 | The two families of simulated curves, Y 1,r and Y 2,r . The gray shaded area represents the first family, and the yellow shaded area is for the second family. The red and blue lines are the first and second mean functions, f 1 and f 2 , respectively. epochs) of the recording with the log spectrum during the late phase (last 60 epochs) of the recording, and identify frequency bands that exhibit any differences between the early vs. late phases. In Section 4.4, we studied the spatial variation of power, at each of the five frequency bands: delta, theta, alpha, beta, and gamma, across all 194 channels, with the goal of identifying regions that exhibit relatively greater proportion of spectral power in each of the five frequency bands of interest. Finally, we compared the spatial variation for each of the five bands during the early vs. late phases of the resting-state EEG recording.

Functional Medians of the Pre-motor Log Spectral Curves
The log of the bias-corrected periodograms at the representative channel (channel 197) that approximately overlies cortex of the pre-motor region recorded for several traces and the FBPs are displayed in Figure 11A. The functional median curve is represented by the black curve, which is located inside the 50% central region, shaded area. The two blue curves outside of the shaded area are the non-outlying maximum and minimum curves. Similar to a FBP, we show in Figure 11B the pointwise boxplot (per frequency point), where the black curve is the median obtained by connecting the medians at each frequency point; the blue curves form the central region (50-th percentile region); the green curves are two non-outlying extreme curves. We compared these two median curves in Figure 11C and noted a slight discrepancy between these median curves derived using a FBP and the pointwise boxplot, with an emphasis on the low frequency range. The main difference between the functional median and the point-wise median curve is in the interpretation. The former is one of the curves from a recorded epoch, whereas the latter may not be an actual curve. Hence the latter cannot really be interpreted as a "typical" curve from a family of curves formed from several epochs. Moreover, the FBPs approach allows us to identify specific epochs that produce "unusual" or outlying log bias-corrected periodogram curves. Note that in the plots, the gray curves are the log bias-corrected periodograms of 160 epochs and the red curves are outliers. Figure 11B also shows that these outlying curves are discontinuous around the frequency bin centered at 100 Hz.

Testing for Stationarity of EEG Epochs Across the Entire Resting-state
In the previous section, the FBP provided descriptive statistics for the log bias-corrected periodograms of 160 epochs from the pre-motor region. Note that there were originally 180 epochs but 20 had to be removed from further analysis due to extra-brain artifact contamination. Our interest now is to test whether resting-state brain activity evolved across the 3 min EEG recording. While there are many ways to characterize such an "evolution" of the underlying brain processes, here we will specifically look into changes on the log spectral curves for early vs. late phases of the resting-state EEG recording. In this case, a change in the log spectral power in early vs. late phases would indicate non-stationarity of the EEG signal across the resting-state recording. The null hypothesis of stationarity here is that the true median curves of the early and last phrases are identical. We test this hypothesis using the rank sum test with the significance level set to 0.05. We defined the early phase to include the first 60 epochs (60 s) of the 3 min recording and the late phase to include the last 60 epochs. In Figure 12, we display the FBPs and the other descriptive statistics for each phase. A visual inspection suggests that the median curves are only slightly different from each other for electrodes that approximately overlie the pre-motor region. More significant differences are noted for electrodes that approximately overlie the prefrontal region (see Figure 12C). Moreover, the rank sum test failed to reject the null hypothesis, as the p-value is 0.56. Therefore, the two median curves are not significantly different and the hypothesis of stationarity in the pre-motor regions is not rejected. This is not entirely unexpected since the whole 3-min recording was purely restingstate. There was no experimental stimulus and the time frame was short.
Next, we use the same testing procedure at this particular channel in the pre-motor region (channel 197) to test the same null hypothesis of non-evolution of the brain process at each of the other channels across the 3 min EEG recording. Among the 194 total channels, 18 channels were identified that demonstrated a significant difference in median curves during the early vs. late phase at a significance level of 0.05. These channels are represented by colored circles in Figure 13. Of these 18, channel 29 (approximately overlying the supplementary motor area) has the lowest p-value at 10 −4 . Since we repeat the same test for 194 channels, we used the Bonferroni correction so that the significance level for each test was set to 0.05/194 = 2 × 10 −4 . Indeed, only channel 29 (anterior supplementary motor area) survived the stringent threshold after the Bonferroni correction. The tests for temporal stationarity at each channel (local spatial tests) revealed several channels having a significant difference between the median curves of the early vs. last phases of the EEG recording. As a next step, we studied stationarity in each of 19 predefined regions of the cortex. In this analysis, the representative EEG signal for each region was obtained by averaging the EEG signal-epochs over all channels within each region.
The plots in Figure 14 suggest that the median curves for the early vs. late phases of the EEG recording are similar for EEG signals recorded from channels that approximately overlie right pre-motor and anterior supplementary motor regions, but different in the right pre-frontal and left parietal regions. Indeed, we conclude from the rank sum test that there is significant difference between the early vs. late phases in cluster of channels that approximately overlie the right pre-frontal (p = 0.01) and the parietal regions (p = 0.029). We found that the right pre-frontal region is significantly non-stationary (i.e., early and late phases differ) at level 0.05 (see Figure 15). This result overlaps with the channel-specific tests, in which several of the channels identified to be non-stationary in the single channel tests are included in the predefined right pre-frontal region. In contrast, while the cluster of electrodes that overlie the left parietal region was found to be non-stationary in the region-by-region tests, none of the 18 channels that were identified to be non-stationary in the single channel tests are part of the left parietal cluster. Therefore, the additional averaging step across group of channels may improve signal-to-noise in this type of analysis. A similar phenomenon was also noted for predefined clusters of electrodes overlying at the left pre-frontal region.

The Variation of Spectral Power at Each Frequency Band Across the Entire Cortex
Our goal here is to test whether the spectral power at each frequency band differed across the cortical surface. We first computed the estimate of the spectral power for each channel at each epoch. Starting with the delta band, for each epoch FIGURE 15 | Testing for difference between the early and late phases of the resting-state for each region. The right pre-frontal regions (blue circles) and the left parietal regions (red circles) exhibit significant non-stationarity at level 0.05.
we construct a 2 − D surface plot of the delta power across the entire cortical surface of 194 channels. These surfaces were then grouped according to the early and late phases of the resting-state. We then applied the surface boxplot method for each frequency band to obtain the median surfaces. In Figure 16, we present the median surface for five frequency bands in the early and late phases. The color blue represents the low spectral power while red is for high power. In Figure 16, it is interesting that even during resting-state there is relatively high spectral power at the beta and gamma bandswhich are both associated with higher cognitive processing (Engel and Fries, 2010). The next step is to test for differences between the early and late phases of the EEG recording for each of the five frequency bands of interest. Using the rank sum test, the delta and alpha bands do not have significant difference between the early and late phases. However, theta, beta and gamma bands show significant differences. In Figure 17, the colored regions indicate significant differences between the first and last phases while the gray color regions indicate no significant differences between these two phases. For the theta band, the rank sum test rejected the null hypothesis at only one region which is the cluster of electrodes overlying anterior supplementary motor. For the beta band, the rank sum test identified differences at the left medial parietal region. For the gamma band, there were 13 regions (out of 19) with significant difference between the early and late phases. Since the gamma band is wider than other bands, an estimated spectrum powers' variation across channel in gamma band is expected to be smaller than the estimated spectrum powers' variation in other bands. In Section 4.3, we tested the stationarity for each region. Figure 15 shows two regions, namely, the right pre-fontal and left parietal, which are significantly non-stationary across all frequencies between the early and late phases. Figure 17 shows that the cluster of electrodes overlying the left parietal region exhibits significant non-stationarity in the beta and gamma bands while the cluster of electrodes overlying the right pre-fontal region is significantly non-stationary only in the gamma band.

Conclusion
This study has extended the use of the classical boxplot to FBP, which is a new visualization tool to analyze functional neuroimaging data, including EEG. The primary findings from the current study demonstrate the FBP is useful for both characterizing the spectral distribution of both simulated and real EEG data and identifying potential outliers in a continuous EEG signal.
In the current implementation of the FBP, ranked sample curves are used to characterize the EEG spectrum by defining a 50% central region, a median curve, and maximum and minimum non-outlying curves. Thus, the shape, size, and length A C B FIGURE 17 | Rank sum test results for regions. Color circles in a bounded curve represents a region, which has a significant difference between the first 60 epochs and the last 60 epochs.
of the FBP can be used to characterize the distribution of the dataset, including the skewedness and degree of variability of the EEG recording. Therefore, potential application of the FBP in this context includes comparing FBPs derived from EEG recordings before and after an experimental intervention (e.g., across a period of motor skill training), comparing mean FBPs derived from EEG recordings in healthy and diseased experimental groups, and comparing mean FBPs derived from EEG recordings during resting-state vs. task.
An additional use of the FBP, as demonstrated by the current results, is to identify potential outliers of the EEG recording. Extra-brain artifacts, including eye blinks, eye movements, heart rhythms measured at pulse points downstream, and muscle movements can cause large deviations in the EEG signal, and represent a significant hurdle in EEG signal processing (Delorme et al., 2007). As a method for identifying outliers in the EEG signal, the FBP could be used to rapidly identify periods of an EEG recording that show high likelihood for contamination by artifacts. In clinical applications, the continuous EEG recording has demonstrated promise as a method for monitoring neural function in patients who have compromised level of consciousness (Fyntanidou et al., 2012) or changes in neural function in patients undergoing neurosurgical interventions (de Vos et al., 2008). The use of FBP to identify outliers in the EEG recording represents a novel method for determining periods of the EEG recording that represent changes in consciousness in patients with a compromised level of consciousness, or for determining changes in neural function across neurosurgical intervention.
The current study also presents an application of the FBP to examine resting-state EEG data acquired from a single individual by comparing EEG signals acquired during early vs. late phases of the 3 min EEG recording. This result has important implications for resting-state studies of neural activity, as many neuroimaging studies that examine resting-state brain function assume resting-state neural activity to be static. However, recent studies that examine dynamic changes in resting-state neural activity suggest momentary change in cognitive processes can cause non-stationarity in resting-state function (Chang and Glover, 2010;Hansen et al., 2015). In contrast, the current results show that the majority of channels demonstrate stationarity across the recording period, and provide support for the assumption that the average EEG signal is static across a 3 min EEG recording. Combined with previous findings, the current results suggest that while momentary changes in cognitive processes result in non-stationary fluctuations of the time series, when averaged across a 60 s subset of the complete 3 min EEG recording, the EEG signal is relatively static. This is supported by the current results that show channels which demonstrate nonstationarity of the EEG signal when comparing early and late phases of the recording include electrodes that overlie the right prefrontal region, which is associated with higher-order cognitive processes (Logue and Gould, 2014). Thus, the assumption of stationarity in resting-state functional neuroimaging studies may be more appropriate for non-cognitive networks, including the motor network. Regardless, further work is needed to determine the minimal time-frame in which EEG signal demonstrate stationarity.
Additional future work is focused on developing a new method for computing confidence bands for the median curve. This method needs to consider the data as a whole. One possible approach is a re-sampling method, in which the notion of band depth is used to construct a 95% confidence band. A potential limitation of the re-sampling method is that there is the potential for multiple curves demonstrating ties with respect to band depth, thus affecting the resultant confidence band. One of the assumptions of the current smoothed periodogram method is that the log bias-corrected periodogram is an unbiased estimator of spectrum. Future work will provide further investigation of this assumption as the current method includes several levels of periodogram manipulation, including smoothing with the gamma generalized crossvalidation, log transformation, and correction by adding Euler Mascheroni constant. In conclusion, the current study presents a novel implementation of the FBP and demonstrates promise as a method for exploratory analysis of complex, high-dimensional neuroimaging datasets, including EEG data.