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 , , . 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   .
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 . 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 . Infants diagnosed with ASD, tend to feature late development of linguistic skills . 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 . 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 . Methods for the analysis of multivariate functional data are more sporadic, but tend to follow similar strategies with specific attention payed to data scales  or potential heterogeneity in the functional evaluation domain . The Bayesian literature on the subject is more sparse, with notable exceptions to be found in , who consider nested spatially correlated functional data under assumptions of separability, and the more general work of , 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 . 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 . 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.
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  has been proposed to pool information within anatomical regions with minimal loss of information. Applications to EEG data have been recently proposed by , and . Our approach follows closely the development in , 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 , 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, 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 .
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 .
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 . 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 , 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 , 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 . 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 . 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 , 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 , .
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 https://github.com/Qian-Li.
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 . 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.
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 . 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.
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 , while gamma waves are thought to index conscious perception and tend to correlate with implicit learning processes .
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 Section6.
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.
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  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. , 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.
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.
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 . 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]http://www.e-publications.org/ims/support/dowload/imsart-ims.zip A user-friendly implementation of the proposed method, with examples, is available online as an Rcpp package at: https://github.com/Qian-Li/HFM.
-  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).
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).
-  Bhattacharya, A., and Dunson, D. B. Sparse bayesian infinite factor models. Biometrika (2011), 291–306.
-  Brillinger, D. R. Time Series: Data Analysis and Theory. Holden-Day Series in Time Series Analysis. Holden-Day, San Francisco, 1981.
-  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.
-  Chen, K., and Lynch, B. Weak separablility for two-way functional data: Concept and test. arXiv preprint arXiv:1703.10210 (2017).
-  Chiou, J. M., Chen, Y., and Yang, Y. Multivariate functional principal component analysis: a normalization approach. Statistica Sinica (2014), 1571–1596.
-  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.
-  Fink, A., and Benedek, M. Eeg alpha power and creative ideation. Neuroscience & Biobehavioral Reviews 44 (2014), 111–123.
-  Franks, A. M., and Hoff, P. Shared subspace models for multi-group covariance estimation. arXiv preprint arXiv:1607.03045 (2016).
-  Gelman, A., Hwang, J., and Vehtari, A. Understanding predictive information criteria for Bayesian models. Statistics and Computing 24, 6 (2014), 997–1016.
-  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.
-  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.
-  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.
-  James, G., Hastie, T., and Sugar, C. Principal component models for sparse functional data. Biometrika 87, 3 (2000), 587–602.
-  Kuhl, P. K. Early language acquisition: cracking the speech code. Nature reviews neuroscience 5, 11 (2004), 831–843.
-  Lang, S., and Brezger, A. Bayesian p-splines. Journal of computational and graphical statistics 13, 1 (2004), 183–212.
-  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.
-  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.
-  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.
-  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.
-  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.
-  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.
-  Rojas, D. C., and Wilson, L. B. -band abnormalities as markers of autism spectrum disorders. Biomarkers in medicine 8, 3 (2014), 353–368.
-  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.
-  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).
-  Uhlhaas, P. J., and Singer, W. Abnormal neural oscillations and synchrony in schizophrenia. Nature reviews neuroscience 11, 2 (2010), 100–113.
-  Wakefield, J. Bayesian and Frequentist regression methods. Springer-NY, 2013.
-  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).
-  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.
-  Zhu, H., Li, R., and Kong, L. Multivariate varying coefficient model for functional responses. Annals of Statistics 40, 5 (2012), 2634–2666.