Region-Referenced Spectral Power Dynamics of EEG Signals: A Hierarchical Modeling Approach

11/13/2018 ∙ by Qian Li, et al. ∙ 0

Functional brain imaging through electroencephalography (EEG) relies upon the analysis and interpretation of high-dimensional, spatially organized time series. We propose to represent time-localized frequency domain characterizations of EEG data as region-referenced functional data. This representation is coupled with a hierarchical modeling approach to multivariate functional observations. Within this familiar setting, we discuss how several prior models relate to structural assumptions about multivariate covariance operators. An overarching modeling framework, based on infinite factorial decompositions, is finally proposed to balance flexibility and efficiency in estimation. The motivating application stems from a study of implicit auditory learning, in which typically developing (TD) children, and children with autism spectrum disorder (ASD) were exposed to a continuous speech stream. Using the proposed model, we examine differential band power dynamics as brain function is interrogated throughout the duration of a computer-controlled experiment. Our work offers a novel look at previous findings in psychiatry, and provides further insights into the understanding of ASD. Our approach to inference is fully Bayesian and implemented in a highly optimized Rcpp package.



There are no comments yet.


page 3

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The human brain and its functional relation to biobehavioral processes like motor coordination, memory formation and perception, as well as pathological conditions like Parkinson’s disease, epilepsy, and autism have been a subject of intense scientific scrutiny [5], [27], [22]. An important and highly prevalent imaging paradigm aims to study macroscopic neural oscillations projected onto the scalp in the form of electrophysiological signals, and record them by means of electroencephalography (EEG).

Currently, a typical multi-channel EEG imaging study is carried out with a geodesic EEG-net, composed of up to 256 electrodes. After precise placement on the scalp, electrodes collect electrophysiological signals at high time-resolution, through event-related potentials (ERP) or event-related oscillations (ERO). In this setting, endogenous or exogenous events can result in frequency-specific changes to ongoing EEG oscillations. The spectral features of the resulting time-series are often used to quantify such changes. Specifically, in a time-frequency analysis, the spatiotemporal dynamics of frequency-specific power are examined as they relate to sensory, motor and/or cognitive processes [26] [12] [19].

Our work is motivated by a study of language acquisition in young children with autism spectrum disorder (ASD), a developmental condition that affects an individual’s communication and social interactions [18]. It is thought that typically developing (TD) infants, as young as 6 months old, start to parse continuous streams of speech to actively discern word patterns [16]. Infants diagnosed with ASD, tend to feature late development of linguistic skills [8]. Because, neither verbal instructions nor behavioral evaluations can be performed on toddlers, EEG has been recognized as an effective imaging tool aimed at the interrogations of the neural substates underlying these complex cognitive processes.

In an idealized setting, for each study unit we aim to characterize variability in the relative log-power of a specific frequency band, recorded at a scalp location , at time . This aim is, however, complexed by several estimation- and data-related challenges. Raw EEG recordings, in fact, suffer from exogenous contamination, poor spatial resolution, as well as non-stationarity. To mitigate these issues, several pre-processing and regularization strategies have been proposed in the literature [26]. In this article, we build on our previous work and couple the characterization of EEG spectral dynamics as region-referenced functional data with a Bayesian model for multivariate functional observations. Fig. 1 illustrates the fundamental nature of the data structures considered in our motivating case study, by depicting specific band power trajectories through time-in-trial and across 3 of 11 brain regions in a study of language acquisition. Given this data structure, we are interested in region-specific, time-varying differences in band power dynamics between diagnostic groups. Crucially, valid inference depends on the correct specification of large multivariate covariance functions.

Many approaches to functional data analysis rely on functional principal components analysis as a fundamental tool for advanced modeling

[29] [15]. Methods for the analysis of multivariate functional data are more sporadic, but tend to follow similar strategies with specific attention payed to data scales [7] or potential heterogeneity in the functional evaluation domain [14]. The Bayesian literature on the subject is more sparse, with notable exceptions to be found in [1], who consider nested spatially correlated functional data under assumptions of separability, and the more general work of [30], who consider similar data structures under more flexible covariance models.

This article introduces a probabilistic framework for multivariate functional data which hinges on two ideas, (1) the constructive definition of Gaussian Processes through bases expansions, and (2) the modeling of covariance operators through functional factor analysis [20]. While neither approach is new, we highlight the application of these concepts to functional brain imaging though EEG. In the process, we discuss important consequences associated with specific implementation choices as they relate to transparent assumptions on the structure of multivariate covariance functions. Our proposal seeks to develop a general framework for automatic shrinkage and regularization through infinite factor models [3]. The resulting methodology is suitable for nonparametric estimation with any projection methods, which is amenable for use in a broad range of applications.

In our work statistical inference is based on posterior simulation through Markov chain Monte Carlo. All proposed techniques are implemented within a highly optimized Rcpp package.

The article is organized as follows. In Section 2, we outline a conceptual framework for region-referenced band powers dynamics. A hierarchical functional model for region-referenced functional data is introduced in Section 3. We assess the performance of the proposed model on engineered datasets in Section 4, and apply our work to a neurocognitive study on auditory implicit learning in Section 5. We conclude with a critical discussion of our work in Section 6.

Figure 1: Power Dynamics in Language Acquisition. Log power dynamics both raw and smoothed, in three of eleven anatomical regions: (LT) left temporal, (RT) right temporal, and (RF) right frontal. For each region, we plot the principal component power for two frequency bands (Alpha, Gamma) during time on trial. We select three subjects from each of the three diagnostic groups: (TD) typically developing, (vASD) verbal ASD, and (mvASD) minimally verbal ASD.

2 Time-dependent Region-referenced Band-Power

EEG signals, measured through high resolution geodesic networks of electrodes, allow for highly localized interrogation of cortical functions throughout the duration of a computer-controlled experiment. In large studies, the high-dimensional nature of the resulting time series poses potential modeling and computational challenges.

Due to spatial proximity, neighboring electrodes collect signals likely to be highly multicollinear. In this context, the analysis of EEG data is often preceded by a dimension reduction exercise, aimed at discarding redundant information and aid interpretation through the definition of region-level summaries. Specifically, electrodes are partitioned according to anatomical regions on the scalp. Subsequently, region-level summaries are defined by averaging power spectra or by selecting the power spectrum of a single electrode within the region. Alternatively, given that multicollinearity of signals defines similarity in their spectral features, spectral PCA [4] has been proposed to pool information within anatomical regions with minimal loss of information. Applications to EEG data have been recently proposed by [23], and [26]. Our approach follows closely the development in [26], who perform spectral PCA on non-overlapping EEG segments. Our target summaries however focus on specific frequency bands as opposed to the entire power spectrum.

We outline and define basic concepts used in our work and defer details to Appendix B of the supplementary document available online. Let denote the raw EEG recording, conceptualized as a locally stationary, zero-mean -dimensional random process observed on subject , , within anatomical region , , for segment , , at a sampling rate across discretized time , . For each observation , FFT is performed to obtain Fourier coefficients , at frequency . The raw periodogram matrix is computed, where is the transpose of the complex conjugate of . Following [23], we define kernel smoothed spectral density matrices , by smoothing across , using a Daniell kernel with bandwidths selected through generalized cross validation (GCV). Crucially, smoothing bandwidth are held fixed within region, resulting in being Hermitian and non-negative definite.

A one-dimensional region-level summary may be obtained by defining the principal power

, as the normalized leading eigenvalue of

, with

. Within region, the principal power summarizes common frequency-level variation across electrodes along the direction of the leading eigenvector. In contrast to

[26], who consider region-referenced longitudinal-functional representations, we adopt a simplified view and focus on specific frequency bands. Specifically, for a frequency range , this article emphasizes modeling the time-varying principal band power .

In adults, typical frequency bands and their spectral boundaries are delta (4 Hz), theta (4-8 Hz), alpha (8-15 Hz), beta (15-32 Hz), and gamma (32-50 Hz). This view is motivated by traditional approaches in neurocognitive science, which differentiate functionally distinct frequency bands as they are thought to target distinct neurocognitive and biobehavioral precesses [25].

3 Hierarchical Models for Region-referenced Functional Data

3.1 Data Projections

Given a frequency band of interest, let be the logarithm of the principal band power for subject , anatomical region , evaluated at time . In practice, for each subject, only a finite number of observations are collected discretely at . However, for ease of notation, and without loss of generality, we assume . Collecting all region-level log band-power measurements into a

-dimensional vector

, we characterize subject-level observations as possibly contaminated multivariate functional data. Specifically, let be a set of -dimensional random functions, with each , , in , a Hilbert space of square integrable functions with respect to the Lebesque measure on . We assume,


with residual covariance . Given , at time ,

is assumed to arise from independent realizations of a heteroschedastic multivariate Gaussian distribution.

We characterize the random functions through finite-dimensional projections onto cubic -spline bases . Specifically, let be a matrix of random basis coefficients, we express 3.1 as follows:


The representation above is readily available for hierarchical modeling. Given priors on and , standard posterior inference and computation applies to the seemingly challenging problem of modeling multivariate functional data. Crucially, if the elements of are assumed Gaussian, given , the simple construction in 3.2 implies that follows a -dimensional Gaussian process, with second-order properties fully determined by covariance assumptions on the random basis coefficients. This observation invites some care in the specification of a prior for , as possibly unwarranted restrictions are easily introduced by naïve constructions.

In the following sections we develop an overarching structure for prior specification, including possible dependence on covariate information, and discuss modeling strategies and their implications for flexibility and efficiency in estimation.

3.2 Hierarchical Priors and Dependence on Covariates

A typical study of functional brain imaging aims to relate subject-level time-stable covariate information , with a region-referenced dynamic outcome . We focus on the conditional expectation as the principal measure of association between covariates and brain function outcomes. Let ; a familiar modeling framework for relies on the following matrix mixed effects construction,


where captures subject specific regional and functional variation. Assuming is organized as a regression design vector, let , , be regression coefficients; a matrix linear model defines . In this setting, the regression functions , are -dimensional varying coefficients, to be interpreted in relation to the design structure encoded in [31].

Given a basis projection, the second order behavior indexing functional and region-level dependence for the the multivariate Gaussian process is fully determined by the covariance matrix . While finite-dimensional, this matrix is likely high-dimensional as its spans both the number of regions , and a potentially large number of basis functions . We discuss three modeling approaches, aimed at regularization through shrinkage and simplifying structural assumptions, namely: (NB) a Naïve-Bayes prior, (SS) a Separable Sandwich prior, and (NS) a regularized Non-Separable prior.

(NB) Naïve-Bayes Prior. This construction exploits two simplifying assumptions: (1) the random coefficients matrix is assumed to follow a Matrix Gaussian distribution, and (2) covariance along time is parametrized through Bayesian P-splines smoothing penalties [17]. Let be a deterministic penalty matrix, and be a covariance matrix indexing dependence between anatomical regions. The NB prior assumes


The prior above structures , implying separability of covariation between regions and time points. The matrix serves both the purpose of modeling dependence between brain regions through its off diagonal elements, and the purpose of establishing region-level adaptive smoothing through its diagonal elements as they multiply

. More details about this matrix and the choice of hyperparameters are discussed in Web Appendix C. While greatly simplified, when compared to a completely unstructured

, this construction is potentially problematic as it hinges on a relatively rigid parametrization for the time-covariance, while enforcing only mild Wishart regularization of .

(SS) Separable Sandwich Prior. A more balanced approach to the prior in 3.4 retains the assumption of separability, but implements adaptive regularized estimation of both the time and region covariance. In particular, we extend the approach of [21], and propose a two-way Bayesian latent factor model for . Let , and be two loading matrices. Also let

, be a random matrix with

, , , we write:


where and , both are diagonal with anisotropic components. This construction, implies

Typically, as in latent factor models, a truncation and is selected to define a low-rank representation of the region and time covariances. Rather than selecting the number of latent factors, we consider the multiplicative shrinkage prior of [3], so that:

When , this prior defines stochastically increasing precisions as more columns are added to and . Specific details about this shrinkage strategy and choice of hyperparameters are reported in Web Appendix ??. A related, but more general construction leading to weakly separable covariance structures is reported in [6]. In what follows, we build on latent factor representations and forego the assumption of separability altogether.

(NS) Non-Separable Prior. Assumptions of separability of the region-by-time covariance are often hard to justify in the setting of functional neuroimaging and should always be tested for adequacy [6]. A simple approach, which foregoes this assumption, while still defining a regularized representation of involves dealing with directly. Specifically, let , be a loading matrix, and , we write


with , implying . This formulation leads to a probabilistic version of the multivariate FPCA of [7], with normalization handled through explicit assumptions of heteroskedasticity through . Echoing the approach used in the sandwich prior, the model is completed with a shrinkage prior on the loading matrix , so that

where indicates the -th -dimensional block, , of the -th column of , .

It is easy to show that the vectorized model in 3.6 encompasses both the separable sandwich prior in 3.5 and the naïve-Bayes prior in 3.4. We examine the finite sample properties of these modeling strategies under several simulation scenarios in the section 4.

3.3 Posterior Inference

Posterior inference is based on Markov Chain Monte Carlo simulations from the target measure. All prior models discussed in section 3.2 are amenable to simple Gibbs Sampling implementations. A highly optimized R package for data manipulation and inference is freely available from github at

We note that the loading matrix in 3.6, as well as and in 3.5 are not likelihood-identified due to invariance to orthogonal rotations. Crucially, however: , and are all uniquely identified, leading to likelihood-identifiability of . Given posterior samples from the model parameters, Monte Carlo estimates of all quantities of interest, including simultaneous credible bands are obtained in a relatively straightforward fashion [2]. Detailed calculations, including full conditional distributions are reported in Web Appendix C.

4 Experiments on Engineered Data

To evaluate the finite sample performance of the hierarchical model and priors described in Section 3 we carried out an extensive Monte Carlo study considering data generated under several dependence structure (separable, non-separable), dependence strength, sample size, and signal-to-noise ratio scenarios. The general goal of these experiments aims to assess how well the mean and covariance are recovered under different priors and simulation truths. Further details regarding the data generation scheme are available in Web Appendix D.

Let and denote posterior mean estimates for their parametric counterparts, we consider the following performance metrics:

  • aiMSE: average integrated mean square error

  • fCOV: Frobenius norm between the estimated and true covariance of : .

In relation to these metrics we observe the following behavior.

Separable (SNR )
NB 0.330 16.600 0.236 16.067 0.154 15.914
SS 0.384 16.597 0.272 15.622 0.175 14.590
NS 0.386 17.577 0.272 16.206 0.174 15.080
(SNR )
NB 0.264 26.247 0.186 26.141 0.119 26.419
SS 0.273 10.333 0.192 9.447 0.122 8.719
NS 0.273 13.492 0.192 11.931 0.122 9.956
Non-Separable (SNR )
NB 0.573 147.594 0.414 146.602 0.275 147.927
SS 0.702 128.988 0.492 125.901 0.315 123.170
NS 0.705 121.671 0.493 113.188 0.314 102.681
(SNR )
NB 0.509 231.030 0.354 213.846 0.230 215.735
SS 0.527 121.137 0.364 119.169 0.235 117.887
NS 0.524 96.578 0.362 86.380 0.234 74.582
Table 1: Simulation Study. Mean and Covariance recovery under Naïve-Bayes (NB), Separable-Sandwich (SS), and Non-Separable (NS) priors. Results are assessed under separable and non-separable simulation truths, sample size escalation () and two signal-to-noise ratio scenarios ().

Sample size and SNR. For each of 100 datasets, we simulate individuals in each of 3 diagnostic groups, with functional observations collected over 6 anatomical regions. Each data-set is considered under two signal-to-noise ratio scenarios, SNR. We assume observations are defined on a common time-grid , with random missing patterns, where between and of the time-points are discarded to mimic observed data. Results are summarized in Table 1.

All three priors recover the mean structure equally well under both SNR scenarios. Independently of the generating dependence structure, aiMSE shows improvement as more subjects are included for analysis. Important differences are observed when we consider recovery of the covariance structure. When data are generated under separable covariance, SS priors exhibit the best performance, with improved accuracy for larger sample sizes. Some loss in efficiency is observed when considering NS priors. However, even with only 10 subjects the loss in efficiency is estimated to be only about 30%, reducing to about 14% when sample size escalates to . When data are generated from non-separably covarying processes, NS priors exhibit the highest efficiency, and are the only model improving meaningfully in accuracy as sample size increases. In all simulation settings, NB priors seem to perform poorly in relation to other alternatives.

In summary, our findings agree with well known results in multivariate and longitudinal data analysis [28]. Consistent recovery of the mean function seems to be independent of the covariance. Inference about the mean, however, requires correct recovery of the dependence structure encoded in . While meaningful gains in efficiency can be achieved if separability assumptions are warranted in applications, the encompassing regularized prior (NS) leads to a generally more flexible and efficient modeling framework, without the need for stringent structural assumptions.

Extended results, including sensitivity to model specification are reported in Web Appendix D. While theoretical large sample results are not examined in the manuscript, all simulation results can be reproduced and extended with our R package and supplementary documentation.

5 EEG and Language Acquisition in TD and ASD Infants

5.1 Study Background

We consider a functional brain imaging study carried out by our collaborators in the Jeste laboratory at UCLA. The study aims to characterize differential functional features associated with language acquisition in TD and ASD children. EEG data were recorded for 144s using an 128 electrode HydroCel Geodesic Sensor Net for 9 TD, 32 ASD children ranging between 4 and 12 years of age. The EEG data is divided into non-overlapping segments of 1.024 seconds, producing a maximum of 140 observable segments for each subject at each electrode. Details about data pre-processing are deferred to Web Appendix A. Individual sensors were partitioned between 11 anatomical regions made up of 4-7 electrodes; left and right for the temporal region (LT and RT) and left, right, and middle for the frontal, central, and posterior regions (LF, RF, MF, LC, RC, MC, LP, RP, and MP, respectively). Region-level power dynamics were estimated as outlined in Section 2. Figure 1 illustrates a sample of smoothed power trajectories within the alpha and gamma

frequency bands for both TD and ASD children. ASD children were further classified as verbal ASD (vASD - 14 children) and minimally verbal ASD (mvASD - 19 children) in relation to their verbal Developmental Quotients (vDQ).

The broad scientific goal of this study aims to understand the functional underpinnings associated with language acquisition through the human ability to parse otherwise continuous speech streams into meaningful words. In order to replicate this natural phenomenon, children were exposed to continuous synthetic speech constructed through a collection of 12 phonemes. By defining phoneme triplets (e.g. pa-bi-ku, da-ro-pi) as deterministic pseudo-words and exposing children to random continuous permutations of pseudo-words, study subjects were given the opportunity for implicit statistical learning of word segmentation. During this process, we seek to detect differential neurophysiological response across anatomical regions and oscillation frequencies.

Figure 2: Alpha Power Dynamics. Average log power dynamics stratified by diagnostic group: (TD) typically developing, (vASD) verbal ASD, and (mvASD) minimally verbal ASD. Four of eleven anatomical regions are represented, namely: (LF) left frontal, (RF) right frontal, (LT) left temporal, and (RT) right temporal. In each region, we display the posterior mean, accompanied by simultaneous credible bands, shaded to include 0.2, 0.6, and 0.8 of all posterior samples.

5.2 Group Mean Trajectory Analysis

Our analysis focuses on time-varying region-referenced differences among diagnostic groups (TD, vASD, mvASD) in relation to two frequency bands, namely alpha (8-15 Hz), and gamma (32-50 Hz). Specifically, alpha waves are thought to play an active role in network coordination and communication [9], while gamma waves are thought to index conscious perception and tend to correlate with implicit learning processes [13].

Our analyses are based on the non-separable model in (3.6), with 10 latent factors, and data projected on 12 B-splines basis functions. Main findings across 4 of 11 regions are summarized graphically in Fig 2 for alpha waves and in Fig 3 for gamma waves. In both models, the varying coefficients representing the mean structure in 3.3

, take as input a simple factorial covariate indexing diagnostic group membership, and yield group-level means across anatomical regions. For each time-varying mean, we also represent simultaneous credible bands by probability shading, to include a gradient of 0.2, 0.6 and 0.9 posterior coverage. A brief discussion about sensitivity to the prior model, number of latent factors, and projection methods is deferred to Section


Starting with a group-by-region analysis of the alpha frequency band in Fig. 2, we note that for all three groups, across all 11 brain regions, there is no discernible time-varying pattern in the average relative alpha frequency, which is maintained relatively flat across time-in-trial. The TD and vASD groups are essentially indistinguishable, with minor differences likely due to sampling variability. Crucially, meaningfully lower levels in the average relative alpha frequency are observed for the mvASD group. While we fall short of considering formal notions of statistical significance, we highlight potentially meaningful findings in the RT region, where 90% simultaneous bands fail to overlap for some of the time-on-trial. More importantly, the lower level in relative alpha frequency is consistent across most brain regions, suggesting a potential role of alpha waves in distinguishing brain functional characteristics in the minimally verbal ASD group. Finally, within group, for all diagnostic groups, we confirm a well known phenomenon known as left-right alpha asymmetry, which is is significant over time at frontal regions, diminished at temporal regions, and negligible at central and posterior regions.

Figure 3: Gamma Power Dynamics. Average log power dynamics stratified by diagnostic group: (TD) typically developing, (vASD) verbal ASD, and (mvASD) minimally verbal ASD. Four of eleven anatomical regions are represented, namely: (LF) left frontal, (RF) right frontal, (LT) left temporal, and (RT) right temporal. In each region, we display the posterior mean, accompanied by simultaneous credible bands, shaded to include 0.2, 0.6, and 0.8 of all posterior samples.

Replicating the same analysis for gamma frequencies in Fig.2 we note that, while alpha waves seem to isolate the mvASD group, gamma waves tend to separate TD children from the vASD and mvASD groups. More precisely, the TD group exhibits significantly higher average gamma band power than vASD and mvASD, especially at the early stages of the experiment. This pattern is most evident at temporal regions, both left and right (LT, RT). Furthermore, the average gamma band power for TD children exhibits a sharp increase at the beginning of time-on-trial, followed by a decreasing trend through the end of the experiment. This finding suggests, that differently from the ASD groups, TD children consciously perceive the beginning of the speech stream as a stimulus. Compatible findings by [13] report a decreased response when stimuli recur, which they projected to be linked to a “neural savings” mechanism within a cell assembly representing an object, i.e. a word in our case. As observed for alpha band power, we note significant frontal asymmetry throughout the study and across diagnostic groups. While previous studies, e.g. [24], pointed out left-dominant asymmetry as a discriminative feature between TD and ASD cohorts, our findings do not replicate this observation within the context of a language acquisition experiment.

Figure 4: Varying Coefficients for Age and vDQ. Varying coefficients associated with age and verbal DQ across all anatomical regions for the Gamma and Alpha frequency bands. Each varying coefficient is accompanied by simultaneous credible bands, shaded to include 0.2, 0.6, and 0.9 of all posterior samples.

5.3 Effects of Age and Verbal-DQ

Rather than considering a coarse classification of subjects into TD, vASD and mvASD, we examine how alpha and gamma band power trajectories change as a function of age and verbal DQ. Due to the relatively small-sized sample, the distribution of subjects demographics among three groups is somewhat unbalanced. To be specific, the TD cohort is significantly older (, in months) with hider vDQ (), the v-ASD cohort are younger () with medium vDQ () and mv-ASD has a wide age-range (), however significantly lower by vDQ (). Our analysis is therefore not perfect, but still warranted under the assumption of generalized additive effects. In particular, the mean structure in 3.3 is represented as the additive combination of time-varying coefficients, including an intercept function, a coefficient function for the main effect of age and a coefficient function for the main effect of vDQ.

Our results are summarized graphically in Fig. 4. We display the varying coefficients for age and vDQ for both the alpha and gamma band power, across all 11 brain regions. For each curve we include probability shading for simultaneous 0.2, 0.6, and 0.9 credible bands. The effects of Age and vDQ on alpha band power is consistent across all brain regions. While age does not seem to be associated with the outcome, higher vDQ levels are consistently and significantly associated with higher alpha band power levels. This finding is consistent with the group-level analysis in Section 5.2, where mvASD children were found to exhibit significantly lower alpha levels, when compared to TD and vASD children. On the other hand, the relationship between vDQ and gamma band power is highly heterogenous across brain regions, with vDQ positively associated with higher gamma band power in temporal regions (LT and RT) and negatively associated with the outcome in the right frontal regions (RF). The effect of Age on gamma power is more consisted across brain regions, with older children exhibiting higher gamma in some regions, and significantly so in LF.

6 Discussion

We introduced a hierarchical modeling framework for the analysis of region-referenced functional data in the context of functional brain imaging through EEG. Our work hinges on two classical ideas, namely: the constructive definition of Gaussian processes through basis functions, and a flexible representation of regularized dependence within and across anatomical regions through latent factor expansions. Our proposal is implemented within a high-performance computation R package, supporting three prior models, which include both separable and non-separable covariance structures for region-referenced functional observations. We showed that the proposed approach has satisfactory operating characteristics under extensive numerical experiments, as well as appealing inferential properties, due to straightforward handling of posterior functionals.

The application of our method relies on projections into the space spanned by basis functions. This feature is engineered into our proposal in order to insure expert knowledge can be included in the analysis by not limiting the applicability of our method to B-spline projections, but potentially including functional spaces spanned, for example, by wavelets or periodic functions, which may be more appropriate in some applications of multivariate functional data analysis. Specific modeling choices, like the number and placement of spline knots, can be based on straightforward calculation of information criteria [11]. Similar considerations apply to the choice of prior model, including restrictions on the structure of large covariance oparators. Our implementation depends on choosing the number of latent factors encoding regularized estimation of covariance operators. From our simulation experiments and data analysis, we conclude that regularization through product Gamma priors yields results which are robust to this specific choice. In practice, the number of latent factors included in the analysis can be compared to the number of principal component functions included in a standard FPCA analysis. A crucial difference in our modeling approach is that we interpret regularized estimation through continuous penalization, rather than discrete truncation to a specific number of principal components. Therefore, we recommend relative overparametrization, by using 5 or more latent factors in default analyses.

Our case study, considering region-referenced EEG data, shows how complex data structures may be analyzed within the familiar hierarchical modeling framework, coupled with varying coefficient models. Our work focused on inference for the mean structure of region-referenced functions. However, important strides can be made in a formal characterizations of mean and covariance dependence on predictors. This is notably relevant in the field of functional brain imaging, where large levels of subject heterogeneity are often observed. We note that, for the covariance structure, simple group comparisons are indeed possible within the proposed framework, by simply running separate analyses. Several questions involving regional and functional correlation are, therefore, easily answered from a straightforward examination of posterior summaries. A promising, perhaps more principled, approach would, for example, include a formal representation of covariance heterogeneity through differential subspace structures


Finally, our work discusses statistical significance only informally, through pairwise comparisons of simultaneous credible bands within anatomical regions. Within the context of multivariate functional data analysis, more work on the meaning and construction of uncertainty bounds is, however, needed, with formal considerations of multiplicity within structured dependence settings.


This work was supported by the grant R01 GM111378-01A1 (DS, DT, CS) from the National Institute of General Medical Sciences.

Supplement Web-Based Supplementary Material [url] A user-friendly implementation of the proposed method, with examples, is available online as an Rcpp package at:


  • [1] Baladandayuthapani, V., Mallick, B., Turner, N., Hong, M., Chapkin, R., Lupton, J., and Carrol, R. Bayesian hierarchical spatially correlated functional data analysis with applications to colon carcinogenesis. Biometrics 64, 64-73 (2008).
  • [2] Baladandayuthapani, V., Mallick, B. K., and Carrol, R. J.

    Spatially adaptive bayesian penalized splines with heteroscedastic errors.

    Journal of Computational and Graphical Statistics 14, 378-394 (2005).
  • [3] Bhattacharya, A., and Dunson, D. B. Sparse bayesian infinite factor models. Biometrika (2011), 291–306.
  • [4] Brillinger, D. R. Time Series: Data Analysis and Theory. Holden-Day Series in Time Series Analysis. Holden-Day, San Francisco, 1981.
  • [5] Broyd, S. J., Demanuele, C., Debener, S., Helps, S. K., James, C. J., and Sonuga-Barke, E. J. Default-mode brain dysfunction in mental disorders: a systematic review. Neuroscience & biobehavioral reviews 33, 3 (2009), 279–296.
  • [6] Chen, K., and Lynch, B. Weak separablility for two-way functional data: Concept and test. arXiv preprint arXiv:1703.10210 (2017).
  • [7] Chiou, J. M., Chen, Y., and Yang, Y. Multivariate functional principal component analysis: a normalization approach. Statistica Sinica (2014), 1571–1596.
  • [8] Eigsti, I.-M., de Marchena, A. B., Schuh, J. M., and Kelley, E. Language acquisition in autism spectrum disorders: A developmental review. Research in Autism Spectrum Disorders 5, 2 (2011), 681–691.
  • [9] Fink, A., and Benedek, M. Eeg alpha power and creative ideation. Neuroscience & Biobehavioral Reviews 44 (2014), 111–123.
  • [10] Franks, A. M., and Hoff, P. Shared subspace models for multi-group covariance estimation. arXiv preprint arXiv:1607.03045 (2016).
  • [11] Gelman, A., Hwang, J., and Vehtari, A. Understanding predictive information criteria for Bayesian models. Statistics and Computing 24, 6 (2014), 997–1016.
  • [12] Gou, Z., Choudhury, N., and Benasich, A. A. Resting frontal gamma power at 16, 24 and 36 months predicts individual differences in language and cognition at 4 and 5 years. Behavioural brain research 220, 2 (2011), 263–270.
  • [13] Gruber, T., and Müller, M. M. Effects of picture repetition on induced gamma band responses, evoked potentials, and phase synchrony in the human eeg. Cognitive Brain Research 13, 3 (2002), 377–392.
  • [14] Happ, C., and Greven, S. Multivariate functional principal component analysis for data observed on different (dimensional) domains. Journal of the American Statistical Association 113, 522 (04 2018), 649–659.
  • [15] James, G., Hastie, T., and Sugar, C. Principal component models for sparse functional data. Biometrika 87, 3 (2000), 587–602.
  • [16] Kuhl, P. K. Early language acquisition: cracking the speech code. Nature reviews neuroscience 5, 11 (2004), 831–843.
  • [17] Lang, S., and Brezger, A. Bayesian p-splines. Journal of computational and graphical statistics 13, 1 (2004), 183–212.
  • [18] Lord, C., Risi, S., Lambrecht, L., Cook, E. H., Leventhal, B. L., DiLavore, P. C., Pickles, A., and Rutter, M. The autism diagnostic observation schedule—generic: A standard measure of social and communication deficits associated with the spectrum of autism. Journal of autism and developmental disorders 30, 3 (2000), 205–223.
  • [19] Mills, D. L., Coffey-Corina, S. A., and Neville, H. J. Language acquisition and cerebral specialization in 20-month-old infants. Journal of Cognitive Neuroscience 5, 3 (1993), 317–334.
  • [20] Montagna, S., Tokdar, S. T., Neelon, B., and Dunson, D. B. Bayesian latent factor regression for functional and longitudinal data. Biometrics 68, 4 (2012), 1064–1073.
  • [21] Montagna, S., Tokdar, S. T., Neelon, B., and Dunson, D. B. Bayesian latent factor regression for functional and longitudinal data. Biometrics 68, 4 (2012), 1064–1073.
  • [22] Murias, M., Webb, S. J., Greenson, J., and Dawson, G. Resting state cortical connectivity reflected in eeg coherence in individuals with autism. Biological psychiatry 62, 3 (2007), 270–273.
  • [23] Ombao, H., and Ho, M.-h. R. Time-dependent frequency domain principal components analysis of multichannel non-stationary signals. Computational statistics & data analysis 50, 9 (2006), 2339–2360.
  • [24] Rojas, D. C., and Wilson, L. B. -band abnormalities as markers of autism spectrum disorders. Biomarkers in medicine 8, 3 (2014), 353–368.
  • [25] Saby, J., and Marshall, P. The utitlity of EEG band power analysis in the study of infancy and early childhood. Developmental neuropsychology 37, 3 (2012), 253–273.
  • [26] Scheffler, A., Telesca, D., Li, Q., Sugar, C., Distefano, C., Jeste, S., and Senturk, D. Hybrid principal components analysis for region-referenced longitudinal functional eeg data. Biostatistics In press (2018).
  • [27] Uhlhaas, P. J., and Singer, W. Abnormal neural oscillations and synchrony in schizophrenia. Nature reviews neuroscience 11, 2 (2010), 100–113.
  • [28] Wakefield, J. Bayesian and Frequentist regression methods. Springer-NY, 2013.
  • [29] Yao, F., Müller, H. G., and Wang, J. L. Functional data analysis for sparse longitudinal data. Journal of the American Statistical Association 100, 577-590 (2005).
  • [30] Zhou, L., Huang, J. Z., Martinez, J. G., Maity, A., Baladandayuthapani, V., and Carroll, R. J. Reduced rank mixed effects models for spatially correlated hierarchical functional data. Journal of the American Statistical Association 105, 489 (2010), 390–400.
  • [31] Zhu, H., Li, R., and Kong, L. Multivariate varying coefficient model for functional responses. Annals of Statistics 40, 5 (2012), 2634–2666.