Modern Magnetic Resonance Imaging (MRI) methodology has made it possible to study structural elements of the brain for small infants. Of special interest is to utilize this technique for the study of how the brain grows and develops in the early childhood years, a period of rapid brain and cognitive growth (Lenroot and Giedd, 2006), during which the brain’s structural and functional networks are formed. The maturity of the developing brain is customarily measured in terms of level of myelination, which can be extracted from MRI signals. Myelin serves as an insulating sheath around nerve fibers and higher myelin levels are thought to be associated with more efficient signaling and thus improved brain function. Altered brain connectivity has been hypothesized to underlie a number of intellectual, behavioral, and psychiatric disorders (Stiles and Jernigan, 2010; Lebel and Beaulieu, 2011; Shaw et al., 2008). To evaluate the development of brain connectivity, the evolution of correlations between myelin levels for various white matter brain regions as a function of age in normal children is of interest. This is because brain development can be characterized by the degree of myelination and time-varying correlation matrices that quantify the dynamic correlations of myelination between different brain regions pinpoint the regions that are developing at similar age periods, characterizing the spatio-temporal brain development map.
A general problem for the analysis of MRI studies in children is the sparsity of data in the time domain, which is a consequence of cost and logistic difficulties. Even in studies that were intended to be longitudinal, for many infants only measurements at one point in time are available. This motivates the development of the approach that we present here, where we obtain consistent estimates of the underlying time-varying correlation matrices from sparse observations, where the estimates are guaranteed to be symmetric and nonnegative definite matrices at each time argument. We consider here the sparsest possible case, namely that each subject is observed at only one random age. Our methodology is of interest for any study where one aims at time-dynamic information for the correlation/covariance structure of continuous multivariate observations, when one essentially has only one measurement per subject available.
In our application to brain development, the underlying process of each subject corresponds to the myelination levels for brain regions, which evolve continuously over time. However, observations are measured only at a single time point, so that the data available for each subject are the (random) time point at which the measurement was taken as well as a multivariate -vector of region-specific myelin levels. As the data are vectors from which the smoothly varying correlation matrix is the desired target, we show here that an adaptation of Fréchet regression (Petersen and Müller, 2018) provides a simple and effective approach to recover the correlation matrices. The proposed method produces for each age an estimate of the population correlation matrix with the following properties: the estimates are correlation matrices themselves; and they are consistent, with convergence rates that can be derived from the general theory of Fréchet regression.
In Section 3, we will demonstrate these methods with data obtained as part of an on-going accelerated longitudinal neuroimaging study, the Brown University Assessment of Myelination and Behaviour Across Maturation (BAMBAM) study (Deoni et al., 2015, 2011)
, by estimating time-varying correlation matrix functions. The cross-correlation matrices of myelination level in dependency on age entail a spatio-temporal map of development that indicates the locations of early and late development, as well as the changing connectivity of various sets of regions of interest. For the children in the BAMBAM study, at each age where a MRI was obtained, the Early Learning Composite (ELC) score is also available, which quantifies neurocognitive development. As an application of the estimated covariance structure between the myelination levels of various brain regions across age, we use this estimated covariance matrix, in addition to a similarly constructed cross-correlation between the ELC outcome and the brain region myelin levels, to obtain linear regressions of ELC versus the actual myelination levels for thebrain regions. Since the coefficients in the linear regression are also time-dependent, without the age-dependent covariance and cross-covariance estimates this regression would not be possible as for each age where predictors (brain myelination levels) and response (ELC) are available, we have only one measurement. Our analysis pinpoints the relevance of various brain regions for neurocognitive outcomes in dependence of the age of infants.
Previous work on modeling and quantifying brain development has not addressed the modeling of cross-correlations and covariances, with the exception of Dai et al. (2017a), where ad hoc kernel smoothing techniques were employed to obtain pairwise cross-correlation functions for the sparsely observed BAMBAM data. However, the ensemble of the recovered pairwise correlations does not form a valid covariance/correlation matrix function, which at each age argument is required to be symmetric and non-negative definite (Petersen and Müller, 2016). Of these properties, symmetry can be easily guaranteed, while non-negative definiteness is a different matter and is not guaranteed for pairwise correlation function estimation approaches.
In this paper, we study kernel smoothing methods for the coherent estimation of these time-varying covariance matrices, where “coherent” means that the estimators are symmetric nonnegative definite matrices at all time points. In Section 2.2, we describe existing smoothing techniques which allow one to estimate the individual elements of each covariance matrix and discuss their drawbacks in the current setting. Sections 2.3 and 2.4 reformulate the problem as the smoothing of entire covariance matrices, opening a connection to a recent method for regression for complex objects, the Fréchet regression model. Section 3 gives more detailed background on the BAMBAM study and the Fréchet estimation technique is applied to study dynamic myelin correlations. We also demonstrate that the estimated covariance matrix functions are instrumental to establish the relation between the myelination levels over the 21 brain regions with concurrent ELC scores. Section 4 validates the approach by demonstrating in a simulation study the ability of the Fréchet estimator to consistently recover dynamic covariances from such sparse data using similar sample size and dimension, and its superiority to an alternative kernel estimator. Section 5 demonstrates theoretical advantages of the Fréchet smoother and extends theory for consistent estimation, which, by necessity, must be performed using biased and degenerate covariance matrices as only multivariate random vectors are observed.
2 Regression with Covariance as Response
We begin by defining the data objects and targets for estimation. The observed data consist of i.i.d. pairs , , where represents the age at which the child was examined and is a vector that contains the myelin water fraction (MWF) levels at brain regions measured during the visit (only available at one age). It is natural to think of as a point observation of a latent smooth multivariate stochastic process , where the are i.i.d. mean-square continuous processes (Ash and Gardner, 1975; Hsing and Eubank, 2015), and Thus, the process tracks the MWF levels of subject at all ages, and is the value of this process at a random age that is independent of the process The targets of interest are the cross-sectional mean and covariance of the latent process, which can be expressed in terms of the conditional mean and covariance of via
Notice that the covariance target is a function , where is the space of symmetric nonnegative definite matrices of dimension . Similarly, one is interested in the corresponding pointwise correlation matrices . Estimators of the component functions of can be obtained in a variety of ways, including smoothing splines and local polynomial methods or other scatterplot smoothers, and have been well explored while this is not the case for the estimation of . Therefore we will focus our attention on the covariance estimation.
2.2 Cross-Covariance Estimation
Established nonparametric techniques could be invoked for the estimation of the individual matrix elements by smoothing so-called raw covariances
Let be a univariate density function, a bandwidth and . Local constant (Nadaraya-Watson) and local linear smoothing are two common tools that can be enlisted to smooth these raw covariances, with respective estimators and , where
As solutions to weighted least squares problems, both of these kernel estimators take the form of a weighted average (Fan and Gijbels, 1996),
where correspond to Nadaraya-Watson and local linear smoothing, respectively. Letting
and , the weights in (2.3) are, respectively,
both satisfying This approach is closely related to the method used in Dai et al. (2017a).
The above described estimators are known to be consistent on the interior , while the local linear estimator is preferable due to its improved performance near the boundaries, which can make a substantial difference in practical applications (Fan and Gijbels, 1996). However, in our application, these estimators possess some undesirable properties when combined into matrix estimators , . First, while both lead to symmetric matrices, neither is guaranteed to be nonnegative definite, so that standard analyses involving covariance matrices cannot be applied without some post-hoc alterations.
In the case of the Nadaraya-Watson estimator, this can in fact be resolved by setting a common bandwidth . Since the weights are strictly positive and the space of covariance matrices is convex, this will give a coherent unified estimator . However, this apparent fix does nothing to resolve poor performance near the boundaries. On the other hand, even with a common bandwidth, the local linear estimator can still fail to be nonnegative definite, since
can be (and often are) negative, especially near the boundaries. If one requires a true covariance, a natural procedure would be to project onto the space of covariance matrices by truncating negative eigenvalues to zero. In the next section, we justify this approach by reframing the covariance estimation problem as a geometric problem involving Fréchet means and an adapted version of a recently developed regression method for responses that are random objects in a metric space.
2.3 Time-varying Covariance Matrices as Conditional Fréchet Means
A key step is to characterize the matrix as a conditional Fréchet mean. The Fréchet mean of a random element of an arbitrary metric space is
where existence can be guaranteed by compactness, but uniqueness cannot. If is a convex subset of a Euclidean space, and is the Euclidean metric, then the ordinary mean and Fréchet mean are known to coincide. Let be a random -vector with zero mean and covariance matrix and set Then, if and is the Frobenius distance, it can be easily shown that the unique Fréchet mean of is due to the fact that is a Euclidean norm as it arises from an inner product.
Now, let and be generic copies of the sample elements and the latent multivariate processes , and Then, for the Frobenius metric the matrix in (2.1) can naturally be expressed as
This characterization suggests that, while smoothing is an appropriate method for estimation of , it should not be done on the scalar raw covariances in (2.2), as this merely targets the individual entries of . Rather, one should smooth entire covariance matrices and target the Fréchet mean As our data consist of multivariate observations, the data objects closest to our targets are the raw covariance matrices
Before we develop the estimator, we provide two remarks that highlight the unique difficulties and challenges associated with this estimation problem.
The covariances are biased in the sense that
This is because they serve only as approximations for the unobservable quantities
Another feature of the raw covariances , which we would still face if we actually observed the unknown in (2.9), is that they are degenerate, i.e., they have rank one. The challenge is then that we are attempting to estimate the object that lies in the interior of using a sample of objects that, by definition, reside on the boundary of the space.
In recent years, the smoothing of covariance matrices has been studied in various contexts, with two prominent examples being the mapping of neural pathways using DTI (Yuan et al., 2012; Carmichael et al., 2013) and the spatial modeling of speech recordings (Tavakoli et al., 2016). In the DTI application, the data consist of a random sample of pairs , where and is a symmetric positive definite matrix, and the target is the conditional Fréchet mean under two particular Riemannian metrics for Being based on the sampling of non-degenerate covariance matrices, due to the availability of much richer data in this context, these previous approaches were not affected by the challenges summarized in Remarks 1 and 2. Furthermore, since Fréchet means generally depend on the chosen metric, in these applications specific features of the chosen metric could be exploited. This is in contrast to our situation, where the goal is to gain an understanding of the time-varying covariability of multivariate data. The canonical metric for this purpose is the Frobenius or metric for covariance matrices.
The importance of estimating time-varying covariance matrices is also evident in the spatial analysis of sound objects in Tavakoli et al. (2016), where the covariance also has a spatial component. Here, the smoothing of covariances takes place only over space and, at least conceptually, the objects that are smoothed are not random but fixed. In practice, these objects need to be estimated, and assumptions are included under which the estimation error is negligible. Such assumptions are common in functional data settings, but do not apply to the cross-sectional data that we consider here, where the raw data are decidedly inconsistent. In short, the nature of our data requires overcoming the challenges that are highlighted in the above remarks, whereas these challenges do not arise in covariance estimation settings that have been previously studied.
Tavakoli et al. (2016) also introduced the so-called -covariance for multivariate data. In our setting, if is a metric on , the -covariance is
The motivation for the -covariance is that other metrics besides the Frobenius metric are more suitable for analysis on , which is a nonlinear Riemannian manifold. As Tavakoli et al. (2016) demonstrated, if is the square-root metric , one can readily use linear methods on the square-root space, which can be identified with the linear space of symmetric matrices. The -covariance thus provides a natural class of covariance objects that can be estimated using Fréchet methods. Alternative metrics that would have similar limitations as the -covariance are those in the Box-Cox class, where , with appropriately defined matrix logarithms for the case (Pigoli et al., 2014; Petersen and Müller, 2016).
However, since if and only if is the Frobenius metric, one cannot smooth under an alternative metric without losing the ability to interpret the target as a true covariance of the observed random vector. In the particular application of these methods to the study of myelination that is our primary motivation, we work with the Frobenius metric for two reasons. First, practitioners are more comfortable with assessing ordinary covariance and correlation between myelin levels due to its well-established interpretation. Second, many statistical models that one might wish to apply to the study of myelin levels necessitate an estimate of their covariance properties, which would not be provided when applying a smoothing procedure under a different metric, as this results in estimates that target a different quantity. The motivating problem of regressing ELC scores, which assess cognitive development, on myelin levels as predictors provides one such instance.
While we consider Fréchet methods under the Euclidean or Frobenius metric, it bears emphasizing that our approach is not restricted to this choice. Indeed, the theory in Section 5 is based on results in Petersen and Müller (2018), which hold under very general conditions and can readily be extended to the Fréchet estimation of -covariances for various metrics .
2.4 Fréchet Estimation
The individual components of this estimator coincide with those in (2.3) for , with weights as in (2.4) and all bandwidths being equal. So, while the standard Nadaraya-Watson smoother does not incorporate the constraints of the space , it results in a standard Fréchet estimate of , due to convexity. As mentioned before, this estimator does not work well near boundaries due to bias, and a local linear type estimator based on weights (2.5) is preferable.
Local linear kernel estimation does not immediately generalize to nonlinear spaces, regardless of whether they are convex or not, as it introduces negative weights near the boundaries, which is necessary to control bias. In Yuan et al. (2012), intrinsic local polynomial regression (ILPR) was proposed using manifold features of under metrics different from Since the Frobenius metric is the natural choice given our target, ILPR is not directly applicable and, more importantly, more complex than necessary. A simpler, yet more general, re-characterization of local linear smoothing for arbitrary metric spaces (Local Fréchet Regression or Fréchet smoothing) has recently been developed in Petersen and Müller (2018). Fréchet smoothing can be easily applied to the current problem of smoothing covariance matrices, where it becomes
and is represented as a weighted Fréchet mean, with weights as in (2.5).
In Fréchet smoothing, we are not weighting the covariances directly, but rather the squared distances which define the Fréchet functional to be minimized. Again, these weights can be negative, since they are derived from the local linear weights in a standard Euclidean smoothing problem, so that is not guaranteed to be in . However, similar to the Nadaraya-Watson estimator, the local Fréchet estimator has an analytic expression in the case of covariance matrices and the Frobenius metric, given by
where is the projection operator onto under
. This projection operator is easy to compute, as it corresponds to truncating negative eigenvalues to zero while keeping the eigenspaces the same. It is interesting, albeit unsurprising, that, whileas defined in (2.12) is clearly an intrinsic estimator as a minimizer of a Fréchet functional, it is also equivalent to the projection onto of the ordinary, extrinsic local linear estimator. This phenomenon is explained by the fact that Fréchet means under Euclidean metrics coincide with ordinary means, along with the uniqueness of Euclidean projections onto closed, convex spaces such as
3 Regional Co-evolution of Brain Myelination in the Developing Brain
3.1 Background and the BAMBAM Study
Altered early brain development is hypothesized to underlie many neurological, behavioral, and intellectual disorders. Infancy and early childhood are sensitive periods of brain growth, coinciding with the emergence of nearly all cognitive, behavioral, and social functioning abilities. Beginning in utero, myelination advances rapidly over the first 2 years of life in a carefully choreographed caudal-cranial, posterior-to-anterior pattern (Brody et al., 1987; Paus et al., 2001), and continues throughout childhood and adolescence (Bartzokis et al., 2010). This maturation pattern is driven by the tight regulation of myelination by neural activity (Fields, 2005; Ishibashi et al., 2006; Lang et al., 2003), which likely underpins its spatio-temporal coincidence with emerging neurobehavior (Fields, 2008; Nagy, Westerberg and Klingberg, 2004; Casey, Galvan and Hare, 2005).
The prolonged nature of myelination imparts a high degree of flexibility and plasticity to developing neural systems. However, this protracted timeline of development also places these systems at prolonged risk to injury or deviant development (Rodier, 1995). It is increasingly recognized that alterations in myelination timing and/or extent can significantly affect behavioral and cognitive outcomes (Fields, 2008), with altered white matter microstructure being a consistent finding in many neurological and neuropsychiatric disorders (Wolff et al., 2014; Flynn et al., 2003; Xiao et al., 2014)
. Despite the importance of coordinated neural communication, relatively little is known regarding the development of structural and functional networks in infants or toddlers or the relationships linking the emergence of neural networks to evolving cognition or neurobehavioral outcomes.
The Brown University Assessment of Myelination and Behavior Across Maturation (BAMBAM) study provides voxelwise data on the myelination for neurotypical children that were very sparsely measured across age, and which we analyze in Section 3.2 to demonstrate some of these issues. Only children with measurements before the age of 1250 days were included, resulting in a sample size of children. Additionally, for those children who were assessed at multiple time points during this period, one observation per subject was selected randomly.
Imaging measures include anatomical -weighted (), quantitative relaxometry (, ), and quantitative myelin water fraction (MWF), among other measurements. The children included in the study were devoid of major risk factors for neurologic and psychiatric disorders. Children with in utero alcohol or illicit substance exposure, premature ( weeks gestation) or multiple birth, fetal ultrasound abnormalities, complicated pregnancy, APGAR scores , NICU admission, neurological disorder, psychiatric or developmental disorders in the infant, parents or siblings were excluded (Chlebowski et al., 2013; Bilenberg, 1999).
The BAMBAM study acquired high-quality and artifact-free MRI data that inform on different characteristics of the developing brain (Fornari et al., 2007; Miller et al., 2012; Hagmann et al., 2010; Grydeland et al., 2013; Shafee, Buckner and Fischl, 2015), including morphometric (whole-brain, white and gray matter volume), myelin water fraction (which informs more specifically on myelination, a key process of brain connectivity), and other modalities (Deoni et al., 2008a; Deoni, Rutt and Peters, 2006; Deoni et al., 2015). Such high quality MRI measurements (see Figure 1) are a prerequisite to assess brain development (Knickmeyer et al., 2008; Deoni et al., 2008b, 2011; Deoni, Rutt and Peters, 2003; Deoni, 2007; Deoni et al., 2004; Deoni, 2011; Deoni, Peters and Rutt, 2004; Deoni and Kolind, 2015). It is also of interest to relate myelin level correlation across brain regions to concurrent developmental and cognitive scores (Shaw et al., 2006), where we use the Early Learning Composite (ELC), derived from the verbal and non-verbal development quotient (VDQ and NVDQ, respectively) scores of the Mullen Scales of Early Learning (Mullen et al., 1995); these scores are standardized for age. It is of special interest to identify the brain regions for which myelination levels are most predictive of these scores at certain ages. This will help to pinpoint the relative relevance of brain regions for neurocognitive performance across age, a goal that has been elusive so far.
The implementation of Fréchet regression that we propose here is a one-step approach that leads directly to bona fide correlation/covariance matrix functions and is supported by theoretical consistency and rate of convergence results that are applicable even for the extremely sparse data case that we consider here. We demonstrate that also in practice this method leads to sensible results when applied to the very sparse data case, where the observed data do not correspond to a sample covariance matrix paired with each predictor, but rather just one random vector of observed levels for each predictor level. Our goal is to quantify the evolution of the cross-correlations with age from data that contain only one measurement per subject. For this, we select discrete anatomical brain regions as an example. The methods generalize to other local aggregations of the MWF levels, which are recorded at the voxel level. These regions are listed in Table 1 together with their acronyms used throughout the remainder of the paper.
While modeling cross-correlation and cross-covariance of myelination has not been well explored, a related topic with more substantial work is the modeling of the myelin trajectories as a function of age. Nonlinear approaches have been applied previously to model overall myelination level trajectories by averaging over brain regions and focusing on the dependence of myelin levels as a function of age (Shaw et al., 2008; Knickmeyer et al., 2008; Casey, Galvan and Hare, 2005; Deoni et al., 2014, 2015). Recently, Dai et al. (2017b) proposed a method to obtain whole brain MWF trajectories, modifying approaches of Sentürk and Müller (2010); Şentürk and Nguyen (2011). For the case where repeated measurements are available for each subject, ideally densely spaced, related work on modeling multivariate trajectories by functional data analysis methods includes Zhou, Huang and Carroll (2008); Chiou and Müller (2016). However, these approaches are not applicable if only one or very few observations per subject are available.
|body Corpus Callosum||bCC||left Cingulum||lC|
|genu Corpus Callosum||gCC||right Cingulum||rC|
|splenu Coropus Collosum||sCC||left Corona Radiata||lCR|
|left Frontal||lF||right Corona Radiata||rCR|
|right Frontal||rF||left Internal Capsule||lIC|
|left Parietal||lP||right Internal Capsule||rIC|
|right Parietal||rP||left Optic Radiation||lOR|
|left Occipital||lO||right Optic Radiation||rOR|
|right Occipital||rO||left Superior Longitudinal Fasciculus||lSLF|
|left Temporal||lT||right Superior Longitudinal Fasciculus||rSLF|
3.2 Age-varying correlations of regional myelination
The BAMBAM study cohort includes children who were observed during the first days of life, with the earliest observation at age 65 days. Some children were observed repeatedly within this time period for a total of visits, although nearly half (48%) of the children have just one measurement and 90% have at most 3 measurements. For this reason, in our analysis we assume that only a single observation per subject is available, which is chosen uniformly from the available observations. For comparison purposes, we also conducted this analysis using the full data set, with corresponding figures being relegated to the Appendix.
We investigated the patterns of MWF development in the regions/tracts listed in Table 1. A subset of the observed data are shown in the left panel of Figure 2, where the children were binned into thirty-five age groups of equal width, and one child per bin was chosen randomly for display. The right panel of Figure 2 also shows the preliminary estimates of the mean MWF levels, using local linear smoothing with a bandwidth of days. This bandwidth was chosen by five-fold cross validation combined across all regions. Comparison with Figure 7 in Petersen, Müller and Deoni (2018) shows similar patterns for estimates obtained from the complete data set.
We then computed cross-sectional covariance estimates as in (2.12); the alternative estimators in (2.3) and (2.11) are not shown due to their inferior performance, which is expected theoretically and empirically demonstrated in the simulations of Section 4. For bandwidth choice, a standard cross validation criterion would not be expected to perform well as the raw covariance matrices are known to live on the boundary of the space (see Remark 2). Letting denote a leave-out estimate of using bandwidth , we consider the following criteria,
where denotes the pseudo inverse of a matrix . As expected, tends to undersmooth the data, whereas
oversmooths due to the inversion. We found that using the geometric meanprovides a reasonable choice of 150 days for the covariance smoothing, using five-fold cross validation.
(for Functional Principal Component Analysis and related issues, see, e.g., Wang, Chiou and Müller,2016).
We visualize the distinct covariance curves , for four ages in Figure 3. We see from the diagonals that variability of MWF levels increases with age, while covariability also scales but to a lesser extent. An interesting finding is the weak dependence of the right occipital region with all other regions, despite having a similar variability pattern. We can also examine the dependency patterns while normalizing the scale by considering the distinct correlation curves , , this time using Functional Principal Component Analysis (FPCA). We note here that the correlation curves are dependent as they are constructed from the same individuals so that many of the basic results such as consistency of FPCA are not applicable. We also ignore the constraint that correlation curves are between -1 and 1 when implementing FPCA. However, none of this precludes us to use FPCA for purely descriptive purposes, and this is the premise for the following FPCA application.
Figure 3(a) shows the pointwise average of these correlation curves along with a functional boxplot. The mean curve indicates that correlation between regional myelin development remains roughly constant, when averaging across all region pairs. Comparing these results with the estimates for the full data set, we again find consistency in the dynamics of covariance estimates (Figure 3 and Figure 8 in Petersen, Müller and Deoni (2018)), as well as similar patterns in the corresponding correlation curves (Figure 4 and Figure 9 in Petersen, Müller and Deoni (2018)).
The first eigenfunctions representing the principal deviations from the average correlation curve observed in the estimated correlations are shown in Figure 3(b); the fractions of variance explained (FVE) by these three curves are , and , respectively. The first eigenfunction has an increasing trend up to day 800, followed by a decreasing trend. This indicates that variability between correlation patterns amongst distinct region pairs is mostly explained by an increased separation from the mean correlation up until approximately 800 days of age, followed by a regression to the mean, i.e., moving towards the average correlation. The eigenfunction curves in Figure 3(b) can be further interpreted by examining the corresponding functional principal component scores in Figure 3(c). As previously observed, the occipital regions have markedly lower correlations with the remaining regions as indicated by the negative scores in this component. It can also be seen that the strongest correlations in region pairs tend to include the Frontal and Parietal regions as well as the Corona Radiata.
We now aim at modeling the relation between myelin development and cognitive ability, as quantified by the ELC score, building on the above approaches. Denoting by a latent process that tracks the ELC score of the th child continuously over the age domain, we are able to actually observe only the cognitive score at the random time at which the child is studied. Recall that is the vector latent process containing the MWF levels for the various regions across age, with mean curve Since both cognition and myelin development are dynamic processes, we consider a varying coefficient linear model
where in our specific application we can substitute fixed intercept value , owing to the fact that ELC scores are standardized to have a mean value of 100 and constant variance.
In view of the linear structure of (3.2), the solution for the slope vector function is
where is the dynamic vector of covariances between myelin development and the ELC score., and is the time-varying covariance matrix that we are targeting with the approach described above.
We then propose to estimate the parameter functions in (3.2) by first estimating the elements , For this, we compute and then apply local linear smoothing for pairs with bandwidth days. Finally, we use the estimate obtained by Fréchet smoothing to compute the regularized ridge estimate
The bandwidths for the estimation of and were set to and , as done previously, and the ridge parameter was chosen to minimize prediction error using five-fold cross validation, resulting in . Estimates and are shown in Figures 4(a) and 4(b). While the covariability between each MWF level and ELC score has similar patterns over time, the interdependencies betwen the MWF levels cause the estimates of the components of the vector function to be considerably more complex.
For example, these estimates indicate a stark contrast in the effects of myelination level in the Genu of the Corpus Callosum, which is associated with a lower cognitive score, versus white matter in the right Parietal lobe, which is associated with a higher score, particularly pronounced at 600 days, and also between the right Optic Radiation, with lower scores associated, and left Corona Radiata, with higher scores associated and most pronounced at 800 days. Despite the aforementioned similarities in mean and covariance estimates between the reduced and full data sets, comparison of the functional coefficients shown in Figure 5 above and Figure 10 in Petersen, Müller and Deoni (2018) demonstrate their sensitivity to changes in due to the inversion in (3.3).
Another meaningful measure we can extract is an estimate of the time-varying coefficient of determination, or multiple correlation,
which we estimate by
with the resulting estimates shown in Figure 4(c), where also the corresponding estimates using the complete data set are shown. One finds that while MWF levels only moderately account for variability in ELC score at younger ages, there is a sharp and steady increase in their predictive power between the ages of two and three.
These exploratory results clearly demonstrate the evolving relationships between maturing brain structure and developing cognitive functioning and behaviour, particularly throughout early neurodevelopment. In older children, adolescents, and adults, neuroimaging results generally conform to the hypothesis that increased brain structure, integrity, and myelin are associated with improved cognitive function. Our results suggest that this relationship is less clear and exhibits complex nonlinearities during infancy and early childhood. For example, our finding that reduced myelination within portions of the Corpus Callosum is associated with increased cognitive functioning runs contrary to preconceptions derived from adult studies (Fryer et al., 2008), however, may be a natural consequence of other prior findings suggesting that slower maturation is also associated with improved cognitive outcomes (Deoni et al., 2014; Erus et al., 2014).
With respect to our finding that brain structure and function relationships become increasingly stronger with age, this too may provide insight into the increasing importance of white matter with age, with grey matter structure playing the more important role in earlier life (Smyser et al., 2016). However, it is also important to note that early cognitive measures suffer high variability and poor predictive ability prior to 12-18 months of age (Slater, 1997), so that this finding may also reflect intrinsic measurement variability. Clearly, in addition to the estimation techniques introduced here, it will be useful in future work to develop inferential methods for testing and confidence sets in order to assess the strength of evidence for these findings.
The simulations reported in this section aim, first, to demonstrate the advantages of local Fréchet smoothing over the Nadaraya-Watson smoother and, second, to establish the reliability of estimates obtained from real data in Section 3.2. The simulations were conducted under the setting , with
being independently sampled from a beta distribution with shape parameters 0.5 and 1.8, which roughly approximates the shape of the age distribution in the BAMBAM data. Then,, , were generated according to the model
In order to emulate the BAMBAM data, for a fixed dimension , the mean vector had components , , where and were drawn once and fixed for all simulations. These parametric mean functions were chosen based on the estimated mean myelin trajectories (see Fig. 2 in Section 3.2). The covariance was formed by generating a matrix with independent random variables in each entry, then computing A second matrix was generated with elements drawn independently as , from which was computed. Finally, with Exp denoting matrix exponentiation and the Hadamard product, we formed
The polynomial factor emulates the increasing variability of MWF levels seen in the BAMBAM data, with minor fluctuations induced by the sine component. Repeated simulations were run by sampling the random variate independently of as a standard -dimensional Gaussian random vector. Mean trajectory estimates were obtained using local linear regression with the bandwidth chosen by ordinary leave-one-out cross validation, pooled across , and separately for each simulation run. Covariance estimation was performed using both Nadaraya-Watson (2.11), local Fréchet (2.12), and standard local linear estimators, which are obtained as in (2.3), with and weights as specified in (2.5). For the latter, the various bandwidths were all taken to be equal, i.e.
To compare the various estimates, for each fixed bandwidth over a grid, we computed the integrated squared error (ISE)
for each simulation run, where stands for a generic estimate of using bandwidth . The results for comparing (2.11), (2.12), and (2.3) using this metric are in Table 2 for 100 simulations, dimensions , and sample sizes The settings we investigated include parameter values that approximate the BAMBAM data, and go beyond by increasing both dimension and sample size. One finds that average ISE values are uniformly lower for the local linear and Fréchet methods, and the very small optimal bandwidths used for the Nadaraya-Watson estimator particularly reflect its well-known inherent bias issues, and increasing sample size leads only to very small declines or no declines in ISE, in contrast to the local linear and Fréchet methods, where increasing sample sizes lead to noticeable improvements. While the tabulated integrated MSE results demonstrate only slight improvement using local Fréchet estimation versus ordinary local linear estimation, in fact the pointwise discrepancies are uniformly lower for the local Fréchet method, with larger discrepancies near the boundaries where the ordinary local linear estimate tends to suffer from influential negative eigenvalues.
The simulation results are visualized in Figure 6, demonstrating the ISE values (in log scale) over all simulation runs rather than just the mean, where we leave out the ordinary local linear estimator. The local Fréchet approach clearly outperforms the Nadaraya-Watson estimator uniformly over all settings, which, again due to its bias, is seen not to display the improved performance with increasing sample size that is found for the local Fréchet estimator and is expected for a reasonable estimation approach.
To further validate our analysis in Section 3, we implemented two additional estimates of the dynamic covariance for the setting and The first alternative estimate is the local Fréchet estimator under the square root metric ,
The appeal of this estimate is that it is an intrinsic estimator that does not require projection onto However, as noted in Section 2.2, this estimator no longer targets the ordinary covariance, but rather the -covariance of Tavakoli et al. (2016) in (2.10). This inherent bias resulted in a worse average ISE of (on the log scale, and minimized over a grid of bandwidths) when compared to local Fréchet estimation under the Frobenius metric in Table 2, although the performance was better than that of the standard Nadaraya-Watson estimate.
A second alternative estimate was considered in order to assess the relative efficiency of using just one observation per subject versus the full data set in the BAMBAM analysis. This estimate was obtained by generating repeated observations for the subjects and then computing the local Fréchet estimate using all observations, ignoring dependencies. To mimic the BAMBAM data, we generated independent obervation counts with probabilities , resulting in a total of 497 observations including repeats, and these numbers were used across all simulation runs. For indices with this required an adaptation of the generative model in (4.1) in order to include dependencies. This was done by generating independent timepoints , according to the beta distribution as before, followed by a zero-mean multivariate normal random vector with and , Lastly, the random vector
was computed. Unsurprisingly, the increased information resulted in an improved average ISE of (again in log scale). Hence, when repeats are available, they can indeed be beneficial. However, the methodology does not require repeats and, most importantly, the challenges mentioned in Section 2.2 associated with single observation data are not alleviated for very sparse longitudinal data such as those available in the BAMBAM study.
5 Theoretical Justifications
We analyze the pointwise behavior of the estimate and establish the rate of convergence for the metric If were known, so that in (2.9) can be computed, the estimator
corresponds to the estimator studied in Petersen and Müller (2018). Thus, we can make use of the inequality
Suppose that the marginal densities of and of exist, the latter for all . Let be the standard Euclidean norm on and be a bandwidth sequence. We require the following assumptions.
The kernel function
is a symmetric probability density function with support
The density of is twice continuously differentiable and
There is a constant such that almost surely. Futhermore,
The function is continuous with respect to both arguments and twice continuously differentiable with respect to on the interior.
For any and a null sequence the auxiliary estimates satisfy
Assumption (A1) is common for estimation of a function with bounded support, but can be relaxed by controlling the tail behavior of . Assumption (A2) is also widely used in local polynomial estimation settings, while (A3) is a natural assumption which has the implication that the conditional density of is well-behaved. Assumption (A4) can also be relaxed by controlling the pointwise tail behavior of the process . The required rate in (A5) is needed to control the error incurred by using in place of in the estimation, where the upper bound is known to hold for local linear estimators under certain assumptions, which include that the functions are twice continuously differentiable and conditions on the bandwidth (Mack and Silverman, 1982; Fan and Gijbels, 1996).
Let and Under assumptions (A1)–(A5), for any interior point the local Fréchet estimator satisfies
Let and be the correlation matrices corresponding to and , respectively. Under the assumptions of Theorem 1, for any
From the theorem and corollary, the bandwidth choice which leads to the fastest rate of convergence is and the resulting rate is . The first term corresponds to the optimal rate of pointwise convergence for nonparametric estimation of regression functions that are twice continuously differentiable. A proof of Theorem 1 is provided below, while that of the corollary is straightforward and is omitted.
Proof of Theorem 1.
First, assumptions (A1)–(A4) can be used to verify the conditions necessary to invoke Theorems 3 and 4 of Petersen and Müller (2018), with the exception that the space is not bounded. With minor alterations to the proofs, assumption (A3) implies that these continue to hold, since is a subset of a Euclidean space, implying
By the contraction property of the projection , it follows that
where is the Frobenius norm. Assumption (A5) also implies that, for large enough , we have
Under (A1), it can easily be shown that where the term is uniform over . Hence, using (A5), we have
- Ash and Gardner (1975) [author] Ash, Robert B.R. B. Gardner, Melvin F.M. F. (1975). Topics in Stochastic Processes. Academic Press [Harcourt Brace Jovanovich Publishers], New York.
- Bartzokis et al. (2010) [author] Bartzokis, GG., Lu, PHP., Tingus, KK., Mendez, MFM., A, RichardR. et al., Peters DGP. D. (2010). Lifespan trajectory of myelin integrity and maximum motor speed. Neurobiology of Aging 31 1554–1562.
- Bilenberg (1999) [author] Bilenberg, NielsN. (1999). The Child Behavior Checklist (CBCL) and related material: standardization and validation in Danish population based and clinically based samples. Acta Psychiatrica Scandinavica 100 2–52.
- Brody et al. (1987) [author] Brody, BAB., Kinney, HCH., Kloman, ASA. Gilles, FHF. (1987). Sequence of central nervous system myelination in human infancy. I. An autopsy study of myelination. J. Neuropathol. Exp. Neurol. 46 283–301.
Carmichael et al. (2013)
[author] Carmichael, OwenO., Chen, JunJ., Paul, DebashisD. Peng, JieJ. (2013). Diffusion tensor smoothing through weighted Karcher means. Electronic Journal of Statistics 7 1913.
- Casey, Galvan and Hare (2005) [author] Casey, BJB., Galvan, AdrianaA. Hare, Todd AT. A. (2005). Changes in cerebral functional organization during cognitive development. Current opinion in neurobiology 15 239–244.
- Chiou and Müller (2016) [author] Chiou, Jeng-MinJ.-M. Müller, H. G.H. G. (2016). A pairwise interaction model for multivariate functional and longitudinal data. Biometrika 103 377–396.
- Chlebowski et al. (2013) [author] Chlebowski, ColbyC., Robins, Diana LD. L., Barton, Marianne LM. L. Fein, DeborahD. (2013). Large-scale use of the modified checklist for autism in low-risk toddlers. Pediatrics 131 e1121–e1127.
- Şentürk and Nguyen (2011) [author] Şentürk, DamlaD. Nguyen, Danh VD. V. (2011). Varying coefficient models for sparse noise-contaminated longitudinal data. Statistica Sinica 21 1831.
- Dai et al. (2017a) [author] Dai, XiongtaoX., Müller, Hans-GeorgH.-G., Wang, Jane-LingJ.-L. Deoni, Sean CLS. C. (2017a). Age-Dynamic Networks and Functional Correlation for Early White Matter Myelination. Preprint.
- Dai et al. (2017b) [author] Dai, XiongtaoX., Hadjipantelis, PantelisP., Wang, Jane-LingJ.-L., Deoni, Sean CLS. C. Müller, Hans-GeorgH.-G. (2017b). Longitudinal Associations Between White Matter Maturation and Cognitive Development Across Early Childhood Xiongta. Preprint.
- Deoni (2007) [author] Deoni, Sean CLS. C. (2007). High-resolution T1 mapping of the brain at 3T with driven equilibrium single pulse observation of T1 with high-speed incorporation of RF field inhomogeneities (DESPOT1-HIFI). Journal of Magnetic Resonance Imaging 26 1106–1111.
- Deoni (2011) [author] Deoni, Sean CLS. C. (2011). Correction of main and transmit magnetic field (B0 and B1) inhomogeneity effects in multicomponent-driven equilibrium single-pulse observation of T1 and T2. Magnetic Resonance in Medicine 65 1021–1035.
- Deoni and Kolind (2015) [author] Deoni, Sean CLS. C. Kolind, Shannon HS. H. (2015). Investigating the stability of mcDESPOT myelin water fraction values derived using a stochastic region contraction approach. Magnetic Resonance in Medicine 73 161–169.
- Deoni, Peters and Rutt (2004) [author] Deoni, Sean CLS. C., Peters, Terry MT. M. Rutt, Brian KB. K. (2004). Determination of optimal angles for variable nutation proton magnetic spin-lattice, T1, and spin-spin, T2, relaxation times measurement. Magnetic Resonance in Medicine 51 194–199.
- Deoni, Rutt and Peters (2003) [author] Deoni, Sean CLS. C., Rutt, Brian KB. K. Peters, Terry MT. M. (2003). Rapid combined T1 and T2 mapping using gradient recalled acquisition in the steady state. Magnetic Resonance in Medicine 49 515–526.
- Deoni, Rutt and Peters (2006) [author] Deoni, Sean CLS. C., Rutt, Brian KB. K. Peters, Terry MT. M. (2006). Synthetic T 1-weighted brain image generation with incorporated coil intensity correction using DESPOT1. Magnetic resonance imaging 24 1241–1248.
- Deoni et al. (2004) [author] Deoni, Sean CLS. C., Ward, Heidi AH. A., Peters, Terry MT. M. Rutt, Brian KB. K. (2004). Rapid T2 estimation with phase-cycled variable nutation steady-state free precession. Magnetic resonance in medicine 52 435–439.
- Deoni et al. (2008a) [author] Deoni, Sean CLS. C., Williams, Steven CRS. C., Jezzard, PeterP., Suckling, JohnJ., Murphy, Declan GMD. G. Jones, Derek KD. K. (2008a). Standardized structural magnetic resonance imaging in multicentre studies using quantitative T 1 and T 2 imaging at 1.5 T. Neuroimage 40 662–671.
- Deoni et al. (2008b) [author] Deoni, Sean CLS. C., Rutt, Brian KB. K., Arun, TarunyaT., Pierpaoli, CarloC. Jones, Derek KD. K. (2008b). Gleaning multicomponent T1 and T2 information from steady-state imaging data. Magnetic Resonance in Medicine 60 1372–1387.
- Deoni et al. (2011) [author] Deoni, Sean CLS. C., Mercure, EvelyneE., Blasi, AnnaA., Gasston, DavidD., Thomson, AlexA., Johnson, MarkM., Williams, Steven CRS. C. Murphy, Declan GMD. G. (2011). Mapping infant brain myelination with magnetic resonance imaging. The Journal of Neuroscience 31 784–791.
- Deoni et al. (2014) [author] Deoni, Sean CLS. C., O’Muircheartaigh, JonathanJ., Elison, Jed TJ. T., Walker, LindsayL., Doernberg, EllenE., Waskiewicz, NicoleN., Dirks, HollyH., Piryatinsky, IreneI., Dean III, Doug CD. C. Jumbe, NLN. (2014). White matter maturation profiles through early childhood predict general cognitive ability. Brain Structure and Function 1–15.
- Deoni et al. (2015) [author] Deoni, Sean CLS. C., Dean, Douglas CD. C., Remer, JustinJ., Dirks, HollyH. O’Muircheartaigh, JonathanJ. (2015). Cortical maturation and myelination in healthy toddlers and young children. NeuroImage 115 147–161.
- Erus et al. (2014) [author] Erus, GurayG., Battapady, HarshaH., Satterthwaite, Theodore DT. D., Hakonarson, HakonH., Gur, Raquel ER. E., Davatzikos, ChristosC. Gur, Ruben CR. C. (2014). Imaging patterns of brain development and their relationship to cognition. Cerebral Cortex 25 1676–1684.
- Fan and Gijbels (1996) [author] Fan, J.J. Gijbels, I.I. (1996). Local Polynomial Modelling and its Applications. Chapman & Hall, London. MR1383587 (97f:62063)
- Fields (2005) [author] Fields, RDR. (2005). Myelination: an overlooked mechanism of synaptic plasticity? Neuroscientist. SAGE Publications 11 528–531.
- Fields (2008) [author] Fields, R DouglasR. D. (2008). White matter in learning, cognition and psychiatric disorders. Trends in neurosciences 31 361–370.
- Flynn et al. (2003) [author] Flynn, SWS., Lang, DJD., Mackay, ALA., Goghari, VV., Vavasour, IMI., Whittall, KPK., Smith, GNG., Arango, VV., Mann, JJJ., Dwork, AJA. et al. (2003). Abnormalities of myelination in schizophrenia detected in vivo with MRI, and post-mortem with analysis of oligodendrocyte proteins. Molecular psychiatry 8 811–820.
- Fornari et al. (2007) [author] Fornari, EleonoraE., Knyazeva, Maria GM. G., Meuli, RetoR. Maeder, PhilippeP. (2007). Myelination shapes functional activity in the developing brain. Neuroimage 38 511–518.
- Fryer et al. (2008) [author] Fryer, Susanna LS. L., Frank, Lawrence RL. R., Spadoni, Andrea DA. D., Theilmann, Rebecca JR. J., Nagel, Bonnie JB. J., Schweinsburg, Alecia DA. D. Tapert, Susan FS. F. (2008). Microstructural integrity of the corpus callosum linked with neuropsychological performance in adolescents. Brain and Cognition 67 225–233.
- Grydeland et al. (2013) [author] Grydeland, HåkonH., Walhovd, Kristine BK. B., Tamnes, Christian KC. K., Westlye, Lars TL. T. Fjell, Anders MA. M. (2013). Intracortical myelin links with performance variability across the human lifespan: results from T1-and T2-weighted MRI myelin mapping and diffusion tensor imaging. The Journal of Neuroscience 33 18618–18630.
- Hagmann et al. (2010) [author] Hagmann, PatricP., Sporns, OlafO., Madan, NeelN., Cammoun, LeilaL., Pienaar, RudolphR., Wedeen, Van JayV. J., Meuli, RetoR., Thiran, J-PJ.-P. Grant, PEP. (2010). White matter maturation reshapes structural connectivity in the late developing human brain. Proceedings of the National Academy of Sciences 107 19067–19072.
- Hsing and Eubank (2015) [author] Hsing, TailenT. Eubank, RandallR. (2015). Theoretical Foundations of Functional Data Analysis, with an Introduction to Linear Operators. John Wiley & Sons.
Ishibashi et al. (2006)
[author] Ishibashi, TT., Dakin, KAK., Stevens, BB., Lee, PRP., Kozlov, SVS. Stewart, CL et al.C. e. a. (2006). Astrocytes promote myelination in response to electrical impulses. Neuron 49 823–832.
- Knickmeyer et al. (2008) [author] Knickmeyer, Rebecca CR. C., Gouttard, SylvainS., Kang, ChaeryonC., Evans, DianneD., Wilber, KathyK., Smith, J KeithJ. K., Hamer, Robert MR. M., Lin, WeiliW., Gerig, GuidoG. Gilmore, John HJ. H. (2008). A structural MRI study of human brain development from birth to 2 years. The Journal of Neuroscience 28 12176–12182.
- Lang et al. (2003) [author] Lang, DJMD., Yip, EE., MacKay, ALA., Thornton, AEA., Vila-Rodriguez, FF. MacEwan, GW et al.G. e. a. (2003). Abnormalities of myelination in schizophrenia detected in vivo with MRI, and post-mortem with analysis of oligodendrocyte proteins. Mol Psychiatry 8 811–820.
- Lebel and Beaulieu (2011) [author] Lebel, CatherineC. Beaulieu, ChristianC. (2011). Longitudinal development of human brain wiring continues from childhood into adulthood. The Journal of Neuroscience 31 10937–10947.
- Lenroot and Giedd (2006) [author] Lenroot, Rhoshel KR. K. Giedd, Jay NJ. N. (2006). Brain development in children and adolescents: insights from anatomical magnetic resonance imaging. Neuroscience & Biobehavioral Reviews 30 718–729.
- Mack and Silverman (1982) [author] Mack, Y. P.Y. P. Silverman, B. W.B. W. (1982). Weak and Strong Uniform Consistency of Kernel Regression Estimates. Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete 61 405–415.
- Miller et al. (2012) [author] Miller, Daniel JD. J., Duka, TetyanaT., Stimpson, Cheryl DC. D., Schapiro, Steven JS. J., Baze, Wallace BW. B., McArthur, Mark JM. J., Fobbs, Archibald JA. J., Sousa, André MMA. M., Šestan, NenadN., Wildman, Derek ED. E. et al. (2012). Prolonged myelination in human neocortical evolution. Proceedings of the National Academy of Sciences 109 16480–16485.
- Mullen et al. (1995) [author] Mullen, Eileen ME. M. et al. (1995). Mullen scales of early learning. AGS Circle Pines, MN.
- Nagy, Westerberg and Klingberg (2004) [author] Nagy, ZoltanZ., Westerberg, HelenaH. Klingberg, TorkelT. (2004). Maturation of white matter is associated with the development of cognitive functions during childhood. Journal of cognitive neuroscience 16 1227–1233.
- Paus et al. (2001) [author] Paus, TT., Collins, DLD., Evans, ACA., Leonard, GG., Pike, BB. Zijdenbos, AA. (2001). Maturation of white matter in the human brain: a review of magnetic resonance studies. Brain Research Bulletin 54 255–266.
- Petersen and Müller (2016) [author] Petersen, AlexanderA. Müller, Hans-GeorgH.-G. (2016). Fréchet integration and adaptive metric selection for interpretable covariances of multivariate functional data. Biometrika 103 103–120.
- Petersen and Müller (2018) [author] Petersen, AlexanderA. Müller, Hans-GeorgH.-G. (2018). Fréchet regression for random objects with Euclidean predictors. Annals of Statistics, to appear (arXiv preprint arXiv:1608.03012).
- Petersen, Müller and Deoni (2018) [author] Petersen, AlexanderA., Müller, Hans-GeorgH.-G. Deoni, SeanS. (2018). Supplement to “Fréchet Estimation of Time-Varying Covariance Matrices From Sparse Data, With Application to the Regional Co-Evolution of Myelination in the Developing Brain”.
- Pigoli et al. (2014) [author] Pigoli, DavideD., Aston, John ADJ. A., Dryden, Ian LI. L. Secchi, PiercesareP. (2014). Distances and inference for covariance operators. Biometrika 101 409–422.
- Rodier (1995) [author] Rodier, Patricia MP. M. (1995). Developing brain as a target of toxicity. Environmental Health Perspectives 103 73.
- Sentürk and Müller (2010) [author] Sentürk, DamlaD. Müller, Hans-GeorgH.-G. (2010). Functional varying coefficient models for longitudinal data. Journal of the American Statistical Association 105 1256-1264.
- Shafee, Buckner and Fischl (2015) [author] Shafee, RebeccaR., Buckner, Randy LR. L. Fischl, BruceB. (2015). Gray matter myelination of 1555 human brains using partial volume corrected MRI images. Neuroimage 105 473–485.
- Shaw et al. (2006) [author] Shaw, PhilipP., Greenstein, DeannaD., Lerch, JasonJ., Clasen, LivL., Lenroot, RhoshelR., Gogtay, NN., Evans, AlanA., Rapoport, JJ. Giedd, JJ. (2006). Intellectual ability and cortical development in children and adolescents. Nature 440 676–679.
- Shaw et al. (2008) [author] Shaw, PhilipP., Kabani, Noor JN. J., Lerch, Jason PJ. P., Eckstrand, KristenK., Lenroot, RhoshelR., Gogtay, NitinN., Greenstein, DeannaD., Clasen, LivL., Evans, AlanA., Rapoport, Judith LJ. L. et al. (2008). Neurodevelopmental trajectories of the human cerebral cortex. The Journal of Neuroscience 28 3586–3594.
- Slater (1997) [author] Slater, AlanA. (1997). Can measures of infant habituation predict later intellectual ability? Archives of disease in childhood 77 474–476.
- Smyser et al. (2016) [author] Smyser, Tara AT. A., Smyser, Christopher DC. D., Rogers, Cynthia EC. E., Gillespie, Sarah KS. K., Inder, Terrie ET. E. Neil, Jeffrey JJ. J. (2016). Cortical gray and adjacent white matter demonstrate synchronous maturation in very preterm infants. Cerebral Cortex 26 3370–3378.
- Stiles and Jernigan (2010) [author] Stiles, JJ. Jernigan, TLT. (2010). The basics of brain development. Neuropsychol. Review 20 327–348.
- Tavakoli et al. (2016) [author] Tavakoli, ShahinS., Pigoli, DavideD., Aston, John ADJ. A. Coleman, JohnJ. (2016). Spatial modeling of Object Data: Analysing dialect sound variations across the UK. arXiv preprint arXiv:1610.10040.
- Wang, Chiou and Müller (2016) [author] Wang, Jane-LingJ.-L., Chiou, Jeng-MinJ.-M. Müller, Hans-GeorgH.-G. (2016). Functional Data Analysis. Annual Review of Statistics and its Application 3 257–295.
- Wolff et al. (2014) [author] Wolff, Jason JJ. J., Gu, HongbinH., Gerig, GuidoG., Elison, Jed TJ. T., Styner, MartinM., Gouttard, SylvainS., Botteron, Kelly NK. N., Dager, Stephen RS. R., Dawson, GeraldineG., Estes, Annette MA. M. et al. (2014). Differences in white matter fiber tract development present from 6 to 24 months in infants with autism. American Journal of Psychiatry.
- Xiao et al. (2014) [author] Xiao, ZhouZ., Qiu, TingT., Ke, XiaoyanX., Xiao, XiangX., Xiao, TingT., Liang, FengjingF., Zou, BingB., Huang, HaiqingH., Fang, HuiH., Chu, KangkangK. et al. (2014). Autism spectrum disorder as early neurodevelopmental disorder: evidence from the brain imaging abnormalities in 2–3 years old toddlers. Journal of autism and developmental disorders 44 1633–1640.
- Yuan et al. (2012) [author] Yuan, YingY., Zhu, HongtuH., Lin, WeiliW. Marron, JSJ. (2012). Local polynomial regression for symmetric positive definite matrices. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 74 697–719.
- Zhou, Huang and Carroll (2008) [author] Zhou, L.L., Huang, J. Z.J. Z. Carroll, R. J.R. J. (2008). Joint modelling of paired sparse functional data using principal components. Biometrika 95 601-619.