{'fraction_non_alphanumeric': 0.03542405603196695, 'fraction_numerical': 0.010057019709962735, 'mean_word_length': 4.922361074182398, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Neural Processes Mixed-Effect Models for Deep Normative Modeling of Clinical Neuroimaging Data

Normative modeling has recently been introduced as a promising approach for modeling variation of neuroimaging measures across individuals in order to derive biomarkers of psychiatric disorders. Current implementations rely on Gaussian process regression, which provides coherent estimates of uncertainty needed for the method but also suffers from drawbacks including poor scaling to large datasets and a reliance on fixed parametric kernels. In this paper, we propose a deep normative modeling framework based on neural processes (NPs) to solve these problems. To achieve this, we define a stochastic process formulation for mixed-effect models and show how NPs can be adopted for spatially structured mixed-effect modeling of neuroimaging data. This enables us to learn optimal feature representations and covariance structure for the random-effect and noise via global latent variables. In this scheme, predictive uncertainty can be approximated by sampling from the distribution of these global latent variables. On a publicly available clinical fMRI dataset, we compare the novelty detection performance of multivariate normative models estimated by the proposed NP approach to a baseline multi-task Gaussian process regression approach and show substantial improvements for certain diagnostic problems.


Introduction
Recently, there has been great interest in applying machine learning methods to neuroimaging data in order to find structural or functional biomarkers for brain disorders. Such biomarkers can potentially be used for diagnosis or predicting treatment outcome in the spirit of precision medicine (Mirnezami et al., 2012). In psychiatry, this is very challenging because clinical groups are highly heterogeneous in terms of symptoms and underlying biology (Kapur et al., 2012). However, most common analysis approaches ignore such heterogeneity and instead consider groups as distinct entities (Foulkes and Blakemore, 2018), e.g., in a case-control setting where subjects are simply labeled as 'patients' or 'controls'. Supervised machine learning methods have been widely used in such settings but their accuracy is limited by the heterogeneity within each disorder (Wolfers et al., 2015).
Normative modeling (Marquand et al., 2016) is an emerging approach to address this challenge that has shown significant promise in multiple clinical settings Zabihi et al., 2018). Intuitively, normative modeling involves estimating variation across the population in terms of mappings between clinically relevant covariates (e.g., age, cognitive scores) and biology (e.g., fMRI activity patterns). This is analogous to the use of 'growth charts' in pediatric medicine to map variation in height or weight as a function of age. Currently, this is implemented using probabilistic regression methods that provide estimates of predictive uncertainty which map variation across the population. Deviations from the resulting normative model can then be interpreted as subject-specific biomarkers for brain disorders. For example, these can be used in a novelty detection setting for predicting diagnosis in an unsupervised fashion .
Accurate quantification of uncertainty is crucial for normative modeling. In the original framework (Marquand et al., 2016), Gaussian process regression (Williams and Rasmussen, 1996) (GPR) was the central tool used to regress neuroimaging measures from clinical covariates. GPR is appealing because it estimates a distribution over functions, providing coherent estimates of uncertainty to map population variation. However GPR also has limitations: it is computationally prohibitive for large datasets and relies on predefined kernels with restricted functional form. Moreover, in the original implementation, brain measures were regressed independently (i.e., in a mass-univariate manner), which does not capitalize on the rich spatial structure of neuroimaging data. This last problem can be addressed by using multi-task GPR (MT-GPR) (Bonilla et al., 2008) to jointly predict multiple brain measurements. However, applying MT-GPR to neuroimaging data is very computationally demanding because of the need to invert large covariance matrices across both space and subjects. Recently, a combination of a low-rank approximations and Kronecker algebra was proposed to scale MT-GPR to whole brain neuroimaging data , which reduces the computational complexity by one order of magnitude. However, this comes with restrictive assumptions that the spatial structures of the signal and noise can be expressed by sets of orthogonal basis functions.
Neural processes (NP) (Garnelo et al., 2018a,b) are latent variable models that bring all the advantages of deep learning (e.g., representation learning and computationally efficient training and prediction) to a stochastic process framework and can address the problems described above. In the NP framework, a distribution over functions is modeled by learning an approximation to a stochastic process. Here, we present an application of NP to multivariate normative modeling of clinical neuroimaging data. This provides three advantages: i) like GPR, NP provides the necessary estimates of predictive uncertainty at test time; ii) similar to MT-GPR, it provides the possibility of learning structured variation; and iii) unlike alternatives, it is computationally scalable without restrictive assumptions on the orthogonality of lower dimensional representations of data. To this end, we make four contributions: i) in a tensor Gaussian predictive process (TGPP) framework (Kia et al., 2018), we define mixed-effect models of neuroimaging data (Friston et al., 1999) as stochastic processes; ii) we show how NPs can be employed for mixed-effect modeling; iii) we use the resulting NP-based mixed-effect model to estimate a normative model of a clinical functional magnetic resonance imaging (fMRI) dataset; iv) we provide an example application of the proposed deep normative modeling for detecting psychiatric disorders in a novelty detection setting. To this end, we apply the proposed NP-based mixed-effect model to a clinical fMRI dataset (Poldrack et al., 2016) in order to jointly predict task-related fMRI brain activity from a set of impulsivity measures. Our experimental results show that the proposed method more accurately identifies ADHD patients from healthy individuals compared to the GP-based normative modeling.

Methods
In this text, we use respectively calligraphic capital letters, A, boldface capital letters, A, and capital letters, A, to denote tensors, matrices, and scalar numbers. We use × n to denote n-mode tensor products. We denote the vertical vector which results from collapsing the entries of a tensor A into a vector with vec(A). Notation |.| is accordingly used to represents the determinant of a matrix or the size of a set.

Mixed-Effect Modeling of MRIs in the TGPP Framework
Consider a neuroimaging study with N subjects and let X ∈ R N ×D to denote the design matrix of D covariates of interest for N subjects. Let Y ∈ R N ×T 1 ×T 2 ×T 3 to represent a 4-order tensor of MRI data for corresponding N subjects with respectively T 1 , T 2 , and T 3 voxels in x, y, and z axes. In the normative modeling setting, we are interested in finding the function f : X → Y. Adopting the tensor Gaussian predictive process (TGPP)  for structured multi-way mixed-effect modeling of MRI data, we have: where A ∈ R D×T 1 ×T 2 ×T 3 represents the fixed-effect across subjects that contains regression coefficients estimated by solving the following linear equations: In Equation (1), Z ∈ R N ×T 1 ×T 2 ×T 3 is the random-effect that characterizes the spatially structured joint variations from the fixed-effect across individuals in different dimensions of MRIs; and E ∈ R N ×T 1 ×T 2 ×T 3 is heteroscedastic noise. Assuming a tensor-variate normal distribution for Y and a zero-mean tensor-variate normal distribution for Z + E, we have: Intuitively, the distribution of the mixed-effect in the joint hypercubic space of clinical covariates and neuroimaging measures can be described as a multi-dimensional Gaussian distribution with vec(Ŷ) and S respectively serving as its mean and covariance.

Mixed-Effect Models of MRI Data as Stochastic Processes
Let (Ω, Φ, ρ) represent a complete probability space (see Oksendal (2003) for definitions) where Ω is a set of neuroimaging measures for N → +∞ subjects (i.e., |Ω| = N ) and Φ is a σ-algebra on Ω that contains all possible subsets of Ω. In other words, Φ represents a set of all possible random neuroimaging samples drawn from the whole human population. Here, ρ : Φ → [0, 1] represents a probability measure that quantifies the probability of occurrence for any entry in Φ. In this setting, each mixed-effect function f i estimated on the ith entry of Φ is a random variable, i.e., a Φ-measurable function from Ω to a Borel set in R N ×T 1 ×T 2 ×T 3 . Therefore, parametrizing f i on different subsets of Ω and X; and considering the exchangeability and consistency properties of mixed-effect models (McCullagh, 2005;Nie and Yang, 2005), can be defined as stochastic processes (Garnelo et al., 2018b). As a corollary, for each entry in Φ, its joint distribution in Equation (3) can be considered as a marginal for a higher-dimensional joint distribution when N → +∞. We exploit this property to frame the problem of mixed-effect modeling in the neural processes framework (Garnelo et al., 2018b). To this end, given a particular realization of the mixedeffect stochastic process f i , the joint distribution in Equation (3) can be rewritten as: In a NP paradigm, we parametrize the integration over all f i (X) on a lower dimensional , resulting the following generative model: where g(X, Z) is a deep neural network that learns the behavior of the mixed-effect model in an amortized variational inference regime. To this end, following the procedure proposed by Garnelo et al. (2018b) the first challenge is to induce stochasticity, for which we need to define 'context' and 'target' functions. The idea is to reduce the difference between the distribution of random context functions from the target function by minimizing their Kullback-Leibler (KL) divergence in the latent space. In our application in order to learn the distribution of the mixed-effect model in Equation (1), we propose to use the estimatedŶ C ∈ R N ×M ×T 1 ×T 2 ×T 3 (using Equation (2)) on M randomly drawn subsets of the training set as context functions. Then, using the actual corresponding neuroimaging training samples as target functions, the following evidence lower-bound should be optimized: where q(Z | X, Y) is the variational posterior of the global latent variable that is parametrized on an encoder h(X,Ŷ C ). In fact in this setting, each context function is a linear component of the target function that roughly approximates a stochastic process f i . Having enough samples of context functions, large enough M , we expect the distribution of context functions to get rich enough to explain non-linear characteristics of the target function (i.e., the mixed-effect f i ). Figure 1 shows a simplified illustration of this scenario in a 2D space where fitting enough linear models on subsets of noisy observations provides an estimation of the distribution of a non-linear target function. Furthermore, by minimizing the KL term in Equation (6), it is expected that the global latent variable Z will learn characteristics of the variance structure of the random-effect and noise terms (the diagonal elements of S) from the difference between the context and target functions (recall that Y −Ŷ = Z + E).

Deep Normative Modeling using Neural Processes
Using NPs in the TGPP framework brings all the advantages of deep learning methods (e.g., representation learning from structured data and computational efficiency) for modeling the multi-way structured variation in neuroimaging data. It has been shown that modeling such structured variation provides the possibility of accurate unsupervised stratification of psychiatric patients in the normative modeling paradigm . To this end, here we introduce deep normative modeling, which utilizes an NP-based mixed-effect modeling and involves following three steps: Figure 1: A schematic illustration on approximating the distribution of a non-linear target function (red curve), e.g., a mixed-effect, from the distribution of linear context functions (blue lines), e.g., fixed-effects, which are fitted on M random subsets of noisy observations (circles).
1. Encoding phase: where an encoder h(X,Ŷ C ) is learned to transfer the covariates, X, and the estimated fixed-effects on M randomly drawn samples from the training set,Ŷ C , to the parameters of the global latent variable Z. Here, to preserve the 3D MRIs structure in the TGPP framework, we propose to use 3D convolutional neural network (3D-CNN) layers to first transfer theŶ C to a lower dimensional representation of neuroimages RŶ ∈ R N ×T . Then, fully connected (FC) layers can be used to derive a latent representation in the joint space of clinical covariates (X) and neuroimages, R ∈ R N ×T . It is worthwhile to emphasize that in this architecture, the aggregation across M context functions is implicitly done by the 3D-CNN layers as they are considered as M input channels to the CNN. Finally, two separate FC layers are used to transfer R to the means (µ Z ∈ R N ×Q ) and standard deviations (σ Z ∈ R N ×Q ) of Z.

2.
Decoding phase: where a decoder g(X, Z) is learned to transfer back the joint covariateslatent space to the neuroimaging data Y. Fully connected and 3D inverse CNN (3D-ICNN) layers can be accordingly used to reconstruct MRIs in the original space.
3. Normative modeling: let Y * ∈ R N * ×T 1 ×T 2 ×T 3 to represent the reconstructed neuroimaging data by the decoder g(X * , Z) on N * test samples. Following Marquand et al.
(2016), we then compute statistical maps describing the deviation for each individual subject from the normative model, referred to as normative probability maps (NPMs), denoted by N ∈ R N * ×T 1 ×T 2 ×T 3 where N = (Y − Y * )/ √ S. Here, S represents the sum of epistemic and aleatoric uncertainties, which respectively describe uncertainty about the true model parameters and inherent variation in the data (Kendall and Gal, 2017). To be able to calculate the epistemic uncertainty in our NP model, we keep the dropout layers active at test time (Gal and Ghahramani, 2016). In the context of mixed-effect modeling of neuroimaging data (in Equation (1)), the aleatoric uncertainty is the byproduct of two factors: i) the across-subject variability which is captured via the covariance of the random-effect Z; and ii) noise in the data which is captured via covariance of E. In the proposed NP framework, these two sources of uncertainties are learned from data and are summarized in the distribution of the global latent variable Z. Therefore, given a test example of clinical covariates x * ∈ X * , we calculate the associated aleatoric uncertainty by sampling from the distribution of Z.

Experimental Material and Setups
In our experiments, we use the response inhibition (i.e., 'stop signal') task from the UCLA Consortium for Neuropsychiatric Phenomics dataset (Poldrack et al., 2016). Specifically,  We use three layers of 3D-CNNs followed by an FC layer to projectŶ C to RŶ . In each CNN layer, we alternate a 3D-convolutional layer, a batch normalization layer (Ioffe and Szegedy, 2015), an average pooling layer, and a leaky ReLU activation function (Xu et al., 2015) (with negative slope of 0.01). Then, two FC layers are used to transfer the merged RŶ and X to the middle joint representation R. A similar reverse architecture is used for the decoder g(X, Z) to transfer back the Z to Y space. Figure 2 depicts a schematic of the employed NP architecture with detailed hyperparameter descriptions. Due to the small sample size and illustrative purpose of our experiments, we did not optimize the architecture and its hyperparameters (e.g., number of layers, number and the size of filters, number of neurons, etc.). The ADAM optimizer (Kingma and Ba, 2014) with decreasing learning rate (from 10 −2 to 10 −5 ) is used for optimization in 100 epochs.
We compare the normative models derived by NP and scalable multi-task Gaussian process tensor regression (sMT-GPTR) , in terms of their accuracy in detecting healthy subjects from patients. In the sMT-GPTR case, we set the number of basis functions across xyz dimensions of data 5, 10 and 3, 5 for the signal and noise, respectively (as they produced the best results in the original study). We evaluate normative modeling accuracy in a novelty detection scenario where we first train a model on a random subset of subjects (75 healthy, 5 SCHZ, 5 ADHD, and 5 BIPL) and then calculate NPMs on a test set of remaining healthy subjects and patients. We emphasize that the model has no access to the diagnostic labels during the training phase and thus our novelty detection approach is completely unsupervised. As in Marquand et al. (2016), we use extreme value statistics to provide a statistical model for the deviations. Specifically, we use a block-maximum approach on the top 1% values in NPMs and fit these to a generalized extreme value distribution (GEVD) (Davison and Huser, 2015). Then for a given test sample and given the shape parameter of GEVD, we compute the value of the cumulative distribution function of GEVD as the probability of that sample being an abnormal sample (Roberts, 2000). Given these probabilities and actual labels, we evaluate the area under the ROC curve (AUC) to measure the performance of the model. All steps (random sampling, modeling, and evaluation) are repeated 10 times in order to estimate the fluctuations of models trained on different training sets. In all these experiments, ordinary least squares are used to estimate the fixed-effect (Equation (2)) on bootstrapped subsets of the training set. 2 Figure 3 compares the AUC of normative models derived by sMT-GPTR and NP. While sMT-GPTR shows slightly better performance in detecting SCHZ patients, NP provides substantially higher accuracy for ADHD cases. The methods perform similarly for BIPL. Considering the fact that these differences in performance are consistent across different model parameters and repetitions, it can be concluded that sMT-GPTR and NP are capturing different characteristics of the underlying biology of impulsivity. Furthermore, the above chance-level detection rates of NP models in ADHD and BIPL confirm a successful application of the proposed NP-based mixed-effect modeling in unsupervised diagnostic prediction. Another important observation in NP models is the ascending trend of the detection performance as the number of samples from the fixed-effect (i.e., M ) increases. This is compatible with the consistency property of mixed-effects as stochastic processes. Figure 4(a) depicts the average difference in NPMs of patient groups from the healthy population for NP(20) model (see Appendix A for supplementary results). Different patterns of deviations from one diagnosis to another shed light on their different underlying biological causes. For example, the sign of deviations changes from SCHZ to ADHD patients in many regions. To further explore the link between these deviations and the level of impulsivity, we computed the coefficient of determination (R 2 ) between the average NPMs (b) R 2 between the impulsivity and deviation from the normative model across different anatomical brain areas ( * * p < 0.01 and * p < 0.1).

Results
in 9 anatomical brain areas and the first principal component of covariates across different diagnostic groups (see Figure 4(b)). The results show significantly (Bonferroni corrected F-test p-values) greater association between impulsivity and deviations in temporal lobes in ADHD and SCHZ patients compared to healthy individuals. This observation is compatible with previous research on the structural and functional engagement of temporal lobes in SCHZ and ADHD (Suddath et al., 1989;Kobel et al., 2010).

Related Work
Rad et al. (2018) used a convolutional autoencoder for unimodal deep normative modeling of human movements recorded by wearable sensors. They used dropout technique in order to evaluate the parameter uncertainty of the model. Our proposed NP-based approach extends their effort in applying deep architectures to normative modeling from two perspectives: i) it provides the possibility of bimodal normative modeling. This is more appropriate for clinical usages where we are generally interested in interpreting the association between clinical covariates and biological measures (Marquand et al., 2016); ii) using the fully probabilistic NP regime, we are also capable to evaluate aleatoric uncertainties resulting from individual differences and noise in addition to the epistemic parameter uncertainty.

Summary and Discussion
In this paper, we proposed a principled approach for estimating spatially structured mixedeffects in neuroimaging data using neural processes. We demonstrated normative modeling as a possible target application for NP-based mixed-effect modeling. Even though the main focus in this study was on neuroimaging data, our contribution in framing the popular mixed-effect modeling as stochastic processes is quite general and opens the door for a wide range of NP applications in different research areas. Moreover, the presented application of NP for deep normative modeling of clinical neuroimaging data brings the advantages of deep neural networks in representation learning to the applications in precision psychiatry. Finally, the computational efficiency of NP in the training and evaluation phases (provided by its reliance on the variational inference) overcomes the lack of computational tractability of the GP-based normative modeling approaches especially when applied to large cohorts of high-dimensional neuroimaging data. For a possible future direction, we consider applying the proposed deep normative modeling approach to a large clinical neuroimaging cohort.