{'fraction_non_alphanumeric': 0.033115961316653356, 'fraction_numerical': 0.012487800550084287, 'mean_word_length': 4.727988819717952, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Joint modelling of association networks and longitudinal biomarkers: an application to child obesity

The prevalence of chronic non-communicable diseases such as obesity has noticeably increased in the last decade. The study of these diseases in early life is of paramount importance in determining their course in adult life and in supporting clinical interventions. Recently, attention has been drawn on approaches that study the alteration of metabolic pathways in obese children. In this work, we propose a novel joint modelling approach for the analysis of growth biomarkers and metabolite concentrations, to unveil metabolic pathways related to child obesity. Within a Bayesian framework, we flexibly model the temporal evolution of growth trajectories and metabolic associations through the specification of a joint non-parametric random effect distribution which also allows for clustering of the subjects, thus identifying risk sub-groups. Growth profiles as well as patterns of metabolic associations determine the clustering structure. Inclusion of risk factors is straightforward through the specification of a regression term. We demonstrate the proposed approach on data from the Growing Up in Singapore Towards healthy Outcomes (GUSTO) cohort study, based in Singapore. Posterior inference is obtained via a tailored MCMC algorithm, accommodating a nonparametric prior with mixed support. Our analysis has identified potential key pathways in obese children that allows for exploration of possible molecular mechanisms associated with child obesity.


Introduction
Obesity is a major risk factor for chronic non-communicable diseases such as type-2 diabetes (T2D), metabolic syndrome and cardiovascular diseases. The prevalence of obesity has reached epidemic proportions worldwide and has tripled between 1975 and 2016. In particular in 2016, 39% of adults were overweight and 13% obese [WHO]. Prevalence of obesity in children has also escalated over the years, increasing from 4% in 1975 to over 18% in 2016 among children and adolescents aged 5-19 years [WHO]. Overweight or obesity in childhood is critical as it often persists into adulthood due to both physiological and behavioural factors. Indeed, childhood obesity is associated with increased risks of glucose intolerance, hypertension, dyslipidaemia, insulin resistance, and T2D in adulthood [Freemark, 2010]. Therefore, preventing childhood obesity can help disrupt the incidence of metabolic diseases in later life. Several different mechanisms such as insulin resistance, inflammation and metabolic dysregulation mediate the link between obesity and the risk of metabolic diseases. Indeed, there is a complex interplay between genetic determinants, behavioural and environmental factors which contribute to obesity. Yet, relatively little is known regarding its underlying pathophysiology.
In this work, we investigate the complex metabolic pathways in childhood obesity, combing metabolite concentration data (as measured by NMR spectroscopy) with more traditional clinical makers measuring the growth of the children. Metabolites are small molecules that participate in metabolic reactions and are involved in biochemical pathways associated with metabolism in health and disease [Ellul et al., 2019]. As the prevalence of obesity is rapidly increasing in children and adolescents, metabolomics is a powerful tool to uncover underlying biological mechanisms, to unravel genetic and environmental interactions, to identify therapeutic targets, to facilitate early detection of metabolic diseases and to monitor disease progression. Applying metabolomic techniques in relation to childhood obesity could pave a way in defining biomarkers of future metabolic risk and targets for early detection and intervention. Previous studies in adults have consistently identified metabolic signatures associated with obesity, insulin resistance and T2D. For example, previous research has reported associations between obesity and elevated plasma concentrations of amino acids such as branched-chain amino acids (BCAA, leucine, isoleucine and valine), aromatic amino acids (AA, phenylalanine and tyrosine), gluconeogenesis intermediates and glutamine metabolism which are linked to inflammation of white adipose tissue in obesity [Takashina et al., 2016, Petrus et al., 2020. Although the metabolomic literature on adults suffering from obesity, insulin resistance or T2D presents convincing results, metabolic changes related to obesity in younger populations have been poorly identified and findings are often inconsistent and different from those in adult populations [Balikcioglu and Newgard, 2018]. For instance, in contrast to adults, a study on children and adolescents in Germany does not report association between BCAA levels and obesity [Wahl et al., 2012], while concentrations of medium-and long-chain acylcarnitines (C12:1 and C16:1) are reported as higher in obese as compared to normal weight children. This latter association has been replicated in adults. Moreover, further 12 metabolites (glutamine, methionine, proline, nine phospholipids) were found to be significantly altered in obese children. The identified metabolite markers are indicative of oxidative stress and of changes in sphingomyelin metabolism, in β-oxidation, and in pathways associated with energy expenditure. Contrary to adults, previous studies [Mihalik et al., 2012] show that obese and diabetic children present no evidence of defects in fatty acid or amino acid metabolism as compared to their normal weight peers. In the cohort study Project Viva [Oken et al., 2015], BCAA concentrations have been reported to be higher in obese versus lean children aged 6-10 years [Perng et al., 2014]. Similarly, it is reported [Butte et al., 2015] that the concentrations of BCAAs, glutamate, lysine, tyrosine, phenylalanine, and alanine significantly increase in obese children as compared with normal weight children. However, other amino acids such as asparagine, aspartate, glycine, serine, and histidine levels decrease. These results indicate that childhood obesity influences the composition of the serum metabolome, pointing towards potential biomarkers.
The aim of this work is to identify metabolic signatures of obesity in children with different trajectories of adiposity from 5 to 9 years of age, using data from the Growing Up in Singapore Towards healthy Outcomes (GUSTO) prospective cohort study [Soh et al., 2014]. Metabolome analysis is particularly relevant in Asian populations where the risk of metabolic diseases is higher than in the western population [Misra and Khurana, 2011] and the GUSTO cohort study provides an optimal platform. GUSTO is a deeply phenotyped prospective cohort involving Singaporean mothers and their children, started in 2009 (pre-natal) by recruiting mothers at the first trimester of pregnancy. A wealth of information is available on both mothers and children. In this work we focus on growth profiles of children and their relationship with metabolic outcomes, as well as more traditional risk factors such as demographics and clinical biomarkers. To this end, we propose a joint model for the growth trajectories and anthropometric measures of children from birth to 9 years of age and a set of metabolites measured at age 8 years in children. The anthropometric indicators are obtained with Quantitative Magnetic Resonance (QMR) techniques [Chen et al., 2018], recording the percentages of fat and lean mass in the children's body excluding the contribution from the bones, and by height/weight measurements, used to compute the standardised body mass index (Z-BMI). These growth indicators are recorded at different time points in the children's development: every year from age 5 to 9 for the QMR measures, and at 21 unequally spaced time points for the Z-BMI. These data pose challenges to the statistical analysis given their dimension and missing rates for some of the variables, as not all subjects took part to follow-up visits.
The main contribution of this work is to provide a joint model for three growth markers and metabolic associations, which allows for data-driven clustering of the children and highlights metabolic pathway involved in child obesity. To this end, in a Bayesian framework, we specify a joint nonparametric random effect distribution on the parameters characterising the longitudinal trajectories of obesity and the graph capturing the association between metabolites. The choice of a nonparametric random effect distribution allows for extra flexibility, heterogeneity in the population as well as data-driven clustering of the subjects.
The paper is structured as follows: Section 2 introduces the proposed approach for the joint modelling of growth trajectories and metabolic associations; Section 3 presents posterior inference results obtained when applying the proposed methodology to the GUSTO data, highlighting cluster-specific growth evolution as well as differences in metabolic associations. Section 4 concludes the paper with a discussion. A Supplementary Material file is available, containing additional Figures and Tables, as well as details on the MCMC algorithm. Additionally, this file contains a description of the dataset used in the analysis.

Joint modelling of growth trajectories and metabolites
The main goal of the analysis is to combine information from the longitudinal responses (i.e. the growth curves) and the metabolic variables (observed only at one time point) to gain a better understanding of the children's development. In particular, we develop a joint model where the longitudinal outcomes are flexibly modelled via a nonparametric mixture of Gaussian Processes (GP), while we exploit tools from the Gaussian Graphical Model (GGM) literature to introduce information from metabolites and their interdependencies. These two components are then linked hierarchically by the specification of a suitable prior distribution. Let Y = {Y t : t ∈ R + } be a stochastic process indexed over the positive real line, in this work representing the time component, and taking values in R. Let the realisations of such process be the vectors Y i = (Y i1 , . . . , Y in ), observed at times t = t 1 , . . . , t n (not necessarily equidistant) for subjects i = 1, . . . , N . A possible modelling strategy consists of assuming a Gaussian Process to model each trajectory over time, resulting in a multivariate Gaussian likelihood distribution for the vectors Y i . This modelling strategy is flexible and allows for efficient computations, but would not be able to account effectively for subject heterogeneity, typical of medical studies. This is evident when observing the empirical distribution of the growth indicators, presenting skewness and heavy tails at the different observed time points, shown in Figures 1 and 2 in Supplementary Material. A possibility to overcome these limitations is to model the observations using a mixture of multivariate Gaussian distributions, adopting a flexible mixing measure which in turn involves a suitable temporal dependence structure such as the one offered by the GP. In order to do so, we extend an existing modelling strategy [Gelfand et al., 2005] to our context and thus assume the distribution of the vectors Y i to be an infinite mixture of multivariate Gaussians where the mixing measure is given by a Dirichlet Process (DP) prior [Ferguson, 1973], centred around a stationary GP. The DP defines a probability measures over the space of probability distributions. A constructive definition of the DP is provided by the stick-breaking representation [Sethuraman, 1994]: where δ x (·) is the Dirac's delta measure taking value 1 at the location x, and 0 otherwise. The infinite sequence of locations {ψ j } ∞ j=1 is an i.i.d. sample from a centering measure P 0 , while the infinite sequence of weights {w j } ∞ j=1 is constructed in the following way: where the mass parameter α > 0 controls the dispersion of the process around P 0 . An important feature of the DP making it appealing in applications is its almost sure discreteness, implying the possibility of modelling ties in the sample from this distribution, inducing a partition of the indices {1, . . . , N }. When used as a mixing measure, this implies a partition of the subjects sharing the same value of the mixing parameter. As previously pointed out [Gelfand et al., 2005], this yields a flexible distribution for the vectors Y 1 , . . . , Y N , which is non-Gaussian and non-stationary, but retains the advantageous mathematical and computational properties of the GP (see Supplementary Section 5 for the details of the algorithm). This, in conjunction with the computational tractability of the DP, allows for efficient posterior inference of the proposed approach. In our application, we have S = 3 distinct processes representing the percentages of fat and lean mass in the body and the standardised body mass index (Z-BMI), and therefore we model the vectors Y ins ), for s = 1, . . . , S via a GP with the following process-dependent covariance kernel: (1) with t i = 1, . . . , n s 1 , t j = 1, . . . , n s 2 , s 1 , s 2 ∈ {1, . . . , S} and η 2 is the nugget parameter, present only on the diagonal elements of the kernel. The covariance kernel presents similar features to the widely-used exponential kernel, accounting for the presence of multiple-processes through the scaling factors ξ s , for s = 1, . . . , S. In particular, it includes the positive coefficient σ 2 calibrating the amount of variability in the data, as well as φ 2 regulating the impact of the distance between time points on the correlation between observations, and it is stationary since it only depends on time distance |t i − t j |.
We indicate by K = K s 1 s 2 t i t j the covariance matrix obtained from Eq. (1) at the observed time points. The mean function of the GP is modelled via the inclusion of the subject-specific parameters θ i = (θ s=1 n s , obtained by concatenating the random intercept vectors relative to each longitudinal process. Additionally, a regression term is included in the expression of the mean of the GP, see Eq. (2).
The second component of the data is represented by the metabolite concentrations measured at year 8. These type of data are usually modelled via multivariate distributions, most commonly Gaussian. Of particular interest in this analysis is the relationship between the observed metabolites, quantifiable by their correlation structure, with the aim of understanding the dependencies between them and their role in the activation of specific metabolic pathways. We approach this problem by modelling the vectors of metabolites borrowing from the GGMs literature. In this setting, let M = (M 1 , . . . , M p M ) ∈ R p M be a vector of p M metabolites and let G = (V, E) be a graph defined over the set of nodes V = {1, . . . , p M } and with edge set E ⊂ {(h, k) ∈ V × V |h < k} such that if there is a connection between the nodes h and k, then (h, k) ∈ E. The graph G is used to represent the correlation structure of the vector of metabolites M , exploiting the property that two elements of the vector are conditionally independent given the rest if and only if the precision matrix is null at the corresponding position [Dempster, 1972]. A zero in the precision matrix corresponds to a zero in the adjacency matrix, which results in the absence of an edge in the graph (and vice-versa, an edge in the graph corresponds to a non-zero element in the precision matrix). The vectors of metabolites are modelled using a multivariate Gaussian distribution with precision matrix Ω G , whose prior distribution is defined conditionally to the graph structure G. The standard conditionally conjugate prior distribution is the G-Wishart [Roverato, 2002] with ν degrees of freedom, scale matrix Ψ and graphical encoding G, denoted here as G-Wishart(Ω G |ν, Ψ, G). The prior distribution over the graph structure is given by the product of i.i.d. Bernoulli priors on each edge with inclusion probability d ∈ (0, 1), so that π(G|d) ∝ d |E| (1 − d) ( p M 2 )−|E| , with |E| being the number of edges in the graph G and p M 2 the total number of possible edges.
As mentioned earlier, the aim of this work is to study the relationship between the growth indicators and the metabolite values jointly, still allowing for flexibility. We specify a joint prior distribution for the hyperparameters of the two sub-models, thus capturing the dependency between the longitudinal and metabolic dimensions. In particular, we specify a joint distribution on the random effect vector θ in the longitudinal part of the model, on the precision matrix Ω G and embedded graph G.
. . , N be the array of subject-specific parameters of interest. We specify a the DP prior on the arrays ψ 1 , . . . , ψ N in order to link the two sub-models corresponding to the growth trajectories and the metabolite concentrations, obtaining: where N p (Y |µ, Ω) is the p-dimensional Gaussian distribution for the vector Y with mean µ and precision matrix Ω. where the DP prior is imposed on the p Y -dimensional concatenated random effects θ 1 , . . . , θ N in order to introduce dependencies within and between processes via the definition of the GP kernel in Eq. (1), on the graph and on the precision matrix. Finally, a conjugate Gaussian prior distribution is imposed on the mean vector of the GP in the centring measure P 0 . We specify a joint Matrix-Normal prior distribution for the column-wise concatenated matrix of coefficients β Y ∈ n s and q Y is the number of covariates used. The prior distribution for the matrix of coefficients β Y has zero mean matrix and identity covariance matrices. Notice that this setting implies a different regression coefficient being estimated at each observed time point for each covariate included in the model, allowing to capture changes in the temporal effect of risk factors. We write gamma(x|a, b) and inv − gamma(x|a, b) to denote the Gamma and the inverse-Gamma distribution for x > 0 with means a/b and b/(a − 1), respectively. We allow for covariate effects on the metabolite concentrations via the terms β M X M i and use a prior specification for the matrix of coefficients analogous to the one for β Y .
The main contribution of the proposed approach is the ability to cluster individuals based on their growth profiles and metabolic associations through joint modelling of longitudinal and multivariate markers. Modelling of multiple graphs has been proposed before in the Bayesian framework [Peterson et al., 2015, Tan et al., 2017, Shaddox et al., 2020 with groups specified a priori, but not in the context of graph-based (unsupervised) clustering, which is achieved by our modelling strategy. To the best our knowledge, this approach has not been proposed in the statistical literature before. This modelling choice implies that the base measure P 0 of the DP is a mixed measure, due to the presence of the graph structure in the sample from the DP. In particular, P 0 is defined on the product space R p Y × P G × G p M , where G p M represents the space of all possible graphs of dimension p M . In general, when using the DP, the base measure P 0 is chosen to be non-atomic, allowing for the computation of the predictive distributions. However, examples of applications requiring a mixed base measure in the specification of the DP are found in the literature [Dunson et al., 2008, Guindani et al., 2009. As reported in existing work [Canale et al., 2017], the mixed measure setting is not problematic in the case of the DP, since the predictive distributions remain unchanged, and thus posterior inference via a Pólya urn algorithm can still be achieved, and we exploit this result. One of the main features of this approach is that the unique values associated with the clusters are not necessarily different, due to the mixed nature of the centering measure P 0 . For instance, in our setting, ties between the graph structures associated with the clusters can be observed (but not among the random effects or precision matrices).
It is common (especially in medical/epidemiological research) that the set of metabolites M is used as predictor, often in a regression setting. From this perspective, the propose model shows similarities with existing Bayesian nonparametric literature on product partition models with covariates (PPMx) . In particular, it can be shown that the marginal distribution of the random partition induced by the DP measure can be factorised in terms including the covariates (i.e., the metabolites) within each cluster. Additional details are reported in Supplementary Section 1. In principle, this would allow us to devise an algorithm similar to the one proposed in the original PPMx models. In practice, due to its computational burden, such approach is unfeasible for graphs of even moderate sizes. As such, we opt for a conditional algorithm (see Supplementary Section 5), which does not marginalise over the random measure. The algorithm is based on Metropoliswithin-Gibbs sampling, with adaptive steps for those parameters which are non-conjugate in the proposed model. The update of the DP parameters follows a Pólya Urn scheme for non-conjugate models, adjusted for the presence of the non-conjugate graphical structure. The updates of the graph and precision matrix within each cluster are tackled using the Birth-and-Death algorithm of .

Posterior inference
In this Section we present the application of the proposed modelling strategy to the data from the GUSTO cohort. The longitudinal data are composed of N = 227 fat and lean percentages measured at years 5 to 9 of the children, together with the Z-BMI values from birth to year 9 at non-equally spaced time points, such that n 1 = n 2 = 5 and n 3 = 21. The size of the concatenated vectors of growth indicators is p Y = 31. The list of p M = 35 metabolites measured at year 8 for the same subjects is reported in Supplementary Table 1. We apply a logit transformation to the fat and lean percentages, in order to map the observations to the real line, while the Z-BMI is standardised. Many of the metabolites present different ranges of values, skewness, and non-normality. In order to correct for these features, we apply a Box-Cox transformation to each metabolite individually and then standardise the observations component-wise.
Additionally to longitudinal and metabolic information, demographic vari-ables and other clinical characteristics, such as ethnicity, pre-pregnancy maternal BMI, findings from oral glucose tolerance test (OGTT) and gender of the child, are available. The full list of covariates is reported in Supplementary Table 2. The same set of covariates is used in both regression components, therefore X Y = X M and q Y = q M = 14, without any intercept term. The covariates present a relatively low percentage of missing values, as reported in Supplementary Table 2, which are imputed using the R package mice. After imputation, the continuous covariates are standardised. We fix the hyperparameters for τ 2 s , σ 2 , φ 2 , η 2 and ξ s for s = 1, . . . , S so that their prior means and variances are both equal to 1; we set d = 2/(p M − 1) ≈ 0.06 inducing sparsity in the graph structure [Jones et al., 2005]; we fix the mass parameter of the DP [Jara et al., 2007] to α = 0.18, yielding E(K N ) ≈ 2 and V ar(K N ) ≈ 1; the hyperparaemters of the centering measure P 0 are set as µ θ = 0, ν = p M + 2 = 37 and Ψ = 10I p M . The MCMC algorithm is run for 50000 iteration after an initial burn-in period of 100 iterations used to initialise the adaptive steps. Then, after a burn-in of 40000 iterations, 5000 iterations are saved with a thinning of 2.

Posterior Inference on Clustering Allocation
An advantage of the proposed approach is the ability of provide posterior inference on the clustering of the subjects, therefore allowing for the identification of groups of children characterised by specific growth trajectories and metabolic associations. As posterior estimate of the random partition ρ N , we report the clustering configuration minimising the Binder's loss function [Binder, 1978], which corresponds the expected loss derived from the two possible misclassification errors, occurring when a pair of subjects is erroneously clustered together or separately. The use of Binder's criterion yields a partition composed of three clusters of sizes 124, 71 and 32, respectively. We refer to it as the Binder partition, and label the clusters by their decreasing size. The number of clusters identified by Binder's method coincides with the posterior mode of the random variable counting the number of clusters, reported in Figure 3 of Supplementary Material, together with the posterior co-clustering probability for each pair of subjects. The results show little uncertainty in the distribution of the number of clusters and cluster assignment. In order to visualize the data and their estimated partition, we display the mean of the longitudinal growth data within each cluster in Figure 1. As it is evident from the top panel of the Figure, the longitudinal growth patterns display an intuitive clustering structure, separating children with low fat and high lean percentages from those with high fat and low lean percentages. The biggest cluster is composed of children with moderate values of fat and lean percentages. The Z-BMI curves follow a similar pattern, but only after 15 months of age.
As exploratory analysis, we also look at the empirical mean of the metabolite concentrations within each cluster, shown in Figure 2. The three clusters exhibit different mean patterns, highlighting differences in the three groups also on a metabolic level. It is evident that the mean level of almost all metabolites in Cluster 1 is around zero, while in the smaller clusters the distributions are centred away from zero and often in opposite directions between them, supporting the hypothesis that they capture different metabolic mechanisms.
Summarizing the posterior distribution of the latent variables ψ 1 , . . . , ψ N is not trivial, due to the existence of label-switching problems arising when working in a Bayesian nonparametric setting. Therefore, in order to understand the results of the clustering analysis, we run an additional MCMC chain, with the same number of iterations, after fixing the random partition ρ N to be equal to the Binder partition. By doing so, we are able to provide the posterior distribution of the values of ψ within each of the three clusters. We show in Figure 3 the posterior estimates of the graph structures within each of the three clusters (the corresponding estimates of the precision matrices are shown in Supplementary Figure 4). The estimates are the median graphs, obtained selecting the edges whose posterior inclusion probability is greater than 0.5 [ Barbieri and Berger, 2004]. We observe that the number of estimated connections in the graphs is highest in Cluster 1, as well as the intensity of the entries of the corresponding precision matrix (see Supplementary Figure 4). In all clusters, we can identify a group of metabolites linked together, corresponding to fatty acids, phosphogrlycerides, apolipoproteins and cholesterol (see Figure 3 and Supplementary Figure 4, top left corners), while associations between smaller groups of metabolites involving some amino acids and ketone bodies show different patterns in the three clusters. We provide a discussion and suggest a biological interpretation of such differences in Section 3.3.

Posterior Inference on Regression Coefficients
We now discuss posterior inference on the regression coefficients for the three responses. As described in Section 2, we estimate the effect of the covariates on each growth process at different time points. We report in Supplementary  Figures 5, 6 and 7 the posterior means and 95% credible intervals (CI) for the entries of the matrix β Y . We consider as relevant those predictors whose 95% CI does not contain the value zero, highlighted in red in the Figures. Interestingly, the covariates which influence the fat and lean percentages at most of the five time points include gestational age, gender of the child, maternal pre-pregnancy BMI and ethnicity, confirming existing results obtained from the same cohort [Ong et al., 2021]. Moreover, some covariates have different effects across time, such gender and highest education degree of the mother. Similar results on the effect of these covariates on the evolution of the Z-BMI trajectories are also reported in Supplementary Figure 7. Posterior estimates of β M are shown in Supplementary Figure 8.

Differential network analysis
To quantify the differences between the cluster-specific networks presented in Figure 3, we estimate a differential network [Valcárcel et al., 2011, Tan et al., 2017, providing an approach based on the joint posterior distribution of the graph structures within each cluster to establish whether the differences among the cluster-specific networks are statistically relevant. We perform three pair-wise comparisons of the networks characterising the three clusters estimated by minimising the Binder's loss function. A differential network only shows those connections for which the absolute difference between the posterior edge inclusion probabilities of two graphs is greater than 0.9. Specifically, for two clusters k 1 and k 2 we require that |π k 1 ij −π k 2 ij | > 0.9, whereπ k ij is the posterior inclusion probability of the edge between nodes i and j in cluster k, estimated using the MCMC output obtained after fixing the random partition ρ N to the Binder partition, as previously done in the context of cluster-specific network estimation (see Figure 3). The resulting differential networks are shown in Figures 4, 5 and 6 (left panels). The differential networks are characterised by a different number of edges. However, there are key metabolites common to all three differential networks: (i) some amino acids such as glycine; (ii) glycoprotein acetyls; (iii) docosahexaenoic acid (DHA, an omega 3 fatty acid); (iv) lipids (HDL and triglycerides). Al-bumin is involved only in the first two differential networks, while acetate only in the last two.
The roles of the metabolites involved in the differential networks can be further explored by performing an Ingenuity Pathway Analysis (IPA) [Krämer et al., 2014]. IPA generates network maps between molecules (see Figures 4(b), 5(b) and 6(b)), and compares them with multiple metabolic pathways (i.e. the linked series of chemical reactions occurring within a cell) available in the literature, assigning a score to each comparison. The scores are generated based on the negative logarithm of the significance level obtained by performing Fisher's exact hypergeometric test when comparing the estimated differential network and the known pathways in the IPA library. For canonical pathway analysis, values of − log(p-value) > 2 are used to detect significant activation. Figures 7 and 8 show a summary of the pathways identified as significant via the IPA methodology. We identify a total of 13 pathways through IPA that are common between the three differential networks, denoted by 1&2, 2&3, and 1&3 in the Figures. Of the 13 pathways, key ones are tRNA charging, biosynthetic pathways for glycine, glutamate receptor signalling as well as degradation of the aromatic amino acid phenylalanine. The commonality among the three groups indicates an active amino acid biosynthesis machinery with the initiation from tRNA charging which is a requisite for translation and transcription of protein biosynthesis through the binding of amino acids. It has been previously shown that altered tRNA aminoacylation, modification and fragmentation are associated with β-cell failure, obesity and insulin resistance [Arroyo et al., 2021]. All the amino acid pathways common to the three differential networks are associated with obesity as well as metabolic syndrome. In Figure 8, comparison between Clusters 1 and 2 (blue) highlights seven unique pathway via the IPA 8, referring to the comparison between normal and low Z-BMI trajectories (see Figure 1). Prominent among them are leucine and tyrosine degradation pathways, involve in catecholamine biosynthesis. The latter has been previously shown to be related to obesity in children, where catecholamine resistance might promote insulin signalling in adipose tissue thus leading to the increase in lipogenesis [Qi and Ding, 2016]. Another key pathway is the dopamine receptor signalling, which could be associated with the behavioural pattern towards food intake [Benton and Young, 2016]. Dopamine receptor were also reported as the neurotransmitter biomarker in research on obesity [Dang et al., 2016]. The comparisons between differential networks 1&3 (green) and 2&3 (yellow) show clear patterns of insulin metabolism as well as activation of signalling pathways related to dysglycemia. The pattern could provide evidence of early events leading to insulin resistance as well as transition to a hyperglycemic state and onset of obesity as evident from IL-12 [Interleukin 12, Nam et al., 2013], apelin (a peptide) [Dray et al., 2008] and growth hormone (GH) signalling [Høgild et al., 2019]. Studies show IL-12 family cytokines as prospective regulators that could cause insulin resistance due to obesity in tissues and plasma [Nam et al., 2013]. Furthermore, comparison between Clusters 2 and 3 showed primarily amino acid biosynthesis, alanine biosynthesis and degradation along with insulin receptor signalling and maturity onset of diabetes, while 1&3 had eleven pathways, mostly degradative in nature as well as hepatic fibrosis signalling pathway. Comparing common pathways between the differential networks 1&2 and 2&3 shows only five common pathways: 4-hydroxyphenylpyruvate biosynthesis, L-dopachrome biosynthesis, phenylalanine degradation I (Aerobic), pyruvate fermentation to lactate and tyrosine Biosynthesis IV. Phenylalanine degradation I (Aerobic), indicates biosynthesis of tyrosine, a feeder molecule for acetoacetate involved in the synthesis of Acetyl (acetyl coenzyme A), which is important for dietary intake and energy balance. The intersection between the set of pathways highlighted by comparisons 1&2 and 1&3 shows folate metabolism and L-carnitine biosynthesis (Figure 7). The metabolic signalling is characteristic of a transition from normal BMI to an obese phenotype. This might indicate a transition from increase in amino acid biosynthesis/degradation and following appearance of lipid biosynthesis. A list of the metabolites identified as commonly differentially expressed between the three Clusters and associated pathways of activation is reported in Table 1.

Conclusions
We propose a Bayesian semiparametric model enabling clustering of subjects based on both longitudinal trajectories and patterns of metabolic association. The work is motivated by a study on early mechanisms of obesity, but has wider applicability. Excess bodyweight is one of the leading risk factors contributing to the overall disease burden worldwide [McMillen et al., 2009]. Childhood obesity is one of the major health problems in western countries and it is increasingly affecting Asian countries. The excessive accumulation of adipose tissue causes inflammation, oxidative stress, apoptosis and mito-  [Faienza et al., 2019]. Our analysis has identified potential key pathways in obese children in order to explore possible molecular mechanisms associated with child obesity (see Table 1). We identified 13 metabolic pathways common to the three differential networks, the majority of which involves amino acids. An analysis of these associations reveals multiple biochemical pathways such as aromatic amino acid metabolism, branched-chain amino acid metabolism, glutathione metabolism, gluconeogenesis, tricarboxylic acid cycle, anti-inflammatory pathways and lipid metabolism. The analysis shows comprehensive initiation of amino biosynthesis as well as precursor molecule degradation, NAD biosynthesis, TCA cycle responsible for providing feeder molecules to sustain the flux required for fat metabolism through synthesis and degradation of aromatic amino acids as well as precursors for acetyl-CoA. The pathways that are unique to each set are able to filter out lipid pathways responsible for BMI/obesity and dyslipidemia, as well as onset of diabetes.
Finally, our results suggest that alterations in amino acid metabolism may play an important role in adiposity and dyslipidemia in children which may be relevant to the susceptibility of metabolic diseases later in life. Our findings are consistent with recent findings which investigate the relationship between obesity in children and pathways (and their combinations) related with amino acid, lipid and glucose metabolism [Matsumoto et al., 2021].

Supplementary Materials
The Supplementary Material file: GPGGM SM.pdf referenced throughout the manuscript is made available with this paper. bility to type 2 diabetes: cohort study with repeated metabolomics across early life. Diabetes Care, 2020.
David Benton and HA Young. A meta-analysis of the relationship between brain dopamine receptors and obesity: a matter of changes in behavior rather than food addiction? International journal of obesity, 40 (1) WHO. World health organisation -Obesity and overweight. https://www. who.int/news-room/fact-sheets/detail/obesity-and-overweight.
(Accessed: 08-09-2021).     We explore in this Section an interesting feature of the proposed model, linked with the existing literature on product partition models with covariates (PPMx) . In particular, it can be shown that the marginal distribution of the random partition induced by the DP measure can be factorised in terms including the covariates (i.e., the metabolites) within each cluster.
Let ρ N = {C 1 , . . . , C K N } be the partition of the indices {1, . . . , N } induced by the sample (ψ 1 , . . . , ψ N ), where we indicate by C j the j-th cluster of size n j = |C j |, for j = 1, . . . , K N . Considering the metabolites observations M as covariates, and following the PPMx approach, we have that the prior for the partition ρ N is: where C(C j ) is the cohesion of the j-th cluster C j , H(M j ) is the similarity of the metabolites for the subjects belonging to cluster j, denoted as M j := {M i : i ∈ C j }, for j = 1, . . . , K N , and V (N, K N , α) is a constant derived from the marginal distribution of the partition ρ N . Information about the partition is included via the cohesion function C, while the contribution of the covariates to the clustering structure is expressed via the similarity function H, facilitating the clustering of subjects with similar covariates. Under our modelling assumptions it can be shown that (1) is given by: The expression of V (N, K N , α) can be derived from the eppf of the DP with general base measure [Argiento et al., 2019, Pitman, 2006. The integral can be simplified further by integrating out the precision matrix Ω G j , which is a-priori distributed according to a G−Wishart(ν, Ψ, G), yielding to the ratio of normalising constants in (??): where Ψ (j) = Ψ+ i∈C j M i M i , and where the part relative to the precision matrix has been integrated out, leaving only the sum over all possible graphs of size p M . The latter sum is finite and theoretically could be computed for each cluster.

Additional Figures
This Section includes some of the Figures discussed in the main manuscript. Figures 1 and 2 show violin plots of the data used in the analysis, grouped by year of observation. Each Figure refer to a different set of growth indicators, namely the fat and lean percentages (logit-transformed) and the Z-BMI values.

Additional Tables
Tables 1 includes the list of metabolites used in the analysis [Bell et al., 2020]. Table 2 includes the available covariates used in the analysis, in both the longitudinal and metabolic part of the model. The Table also includes the percentages of missing values for each covariate.    Figure 4 shows the posterior estimates of the precision matrices within each of the estimated cluster. The estimates are obtained by running an additional MCMC chain for which the partition of the subjects is fixed to the one estimated by minimising Binder's loss function.

Regression on growth biomarkers
We consider the effect of the covariates on the evolution of the growth indicators over time (fat and lean percentages as well as Z-BMI values). We report the posterior distribution of the corresponding regression coefficients β Y in Figures 5, 6 and 7. We report the posterior 95% CIs of the coefficients by time point, and highlight in red those coefficients whose posterior 95% CI does not contain zero, and are therefore deemed relevant.

Regression on Metabolite Concentrations
We consider also the effect of the covariates on metabolite levels. We report the posterior distribution of the corresponding regression coefficients β M in Figure 8. In this case, we group the posterior 95% CIs of the coefficients by metabolite (in total p M = 35), and highlight in red those coefficients whose posterior 95% CI does not contain zero, considered as relevant in the analysis.
Year 5 Year 6 Year 7 Year 8  Year 5 Year 6 Year 7 Year 8  Year 8 Year 9 Year 1     Metabolites: estimates of β M and 95% CI Figure 8: Posterior mean and 95% credible intervals for the regression coefficient on the metabolites. Each dot represents the posterior mean of the effect of a specific covariate on a a metabolite. The red dots indicate relevant effect of the corresponding covariate on that metabolite.

MCMC algorithm
In this Section, we describe the steps used to implement the MCMC algorithm. Non-conjugate parameters in the model are updated with the adaptive Metropolis-Hastings (MH) algorithm for multivariate variable [Haario et al., 2001], where new candidates are proposed from a multivariate Normal distribution centred on the value of the parameter at the current iteration, and with covariance matrix equal to an appropriately re-scaled version of the sample covariance matrix obtained using the samples produced in the MCMC chain so far. The update of the DP parameters follows a Pólya Urn scheme for non-conjugate models [Neal, 2000, Favaro andTeh, 2013], adjusted for the presence of the non-conjugate graphical structure. We provide in the following the expression of the resulting acceptance probabilities.
• The regression coefficients in the growth trajectories and the metabolites part of the model are assigned the following matrix-Normal prior distribution: p(β) = MN p×q (β|0, I p , I q ) We use an adaptive MH algorithm and accept a proposed value β new according to the following probability where the contribution of the proposal distribution, being symmetric, cancels out: for the growth trajectories and min 1, for the metabolite concentrations.
• We assume a priori τ 2 s ∼ inv − gamma(τ 2 s |a τ 2 s , b τ 2 s ), and therefore: • We assume a priori µ θ ∼ N p Y (µ θ |µ µ θ , Σ −1 µ θ ), and therefore: • σ 2 , φ 2 , η 2 , ξ 1 , . . . , ξ S are all non-conjugate and require an adaptive MH step, which includes the modification of the Gaussian covariance kernel K = K s 1 s 2 t i t j , which is a function of these parameters, yielding the following acceptance probability: • We propose a Pólya urn scheme for the update of the random partition induced by the Dirichlet process (DP) prior. The centring measure P 0 is non-conjugate in the proposed model, thus requiring the use of suitable algorithms [Neal, 2000, Favaro andTeh, 2013]. Notice that the support of the centring measure P 0 is mixed (i.e., contains atoms represented by the graphs, but also a diffuse part for the mean vector of the growth trajectories). Thanks to the theoretical properties of the Dirichlet process discussed in Section 2 of the main manuscript, the steps of the algorithm are analogous of the case of a diffuse centring measure. An important difference, however, is that we expect ties in the values of the parameters associated with the clusters, ψ 1 , . . . , ψ K N , specifically in the values of the graphs.
• Given the partition of the subjects, the unique values of the mean vectors of the growth trajectories θ 1 , . . . , θ K N have a conjugate fullconditionals: where n j =| C j | for j = 1, . . . , K N and τ 2 is a vector composed of the process-specific variances τ 2 s replicated n s times, for s = 1, . . . , S. • The unique values in the graph structure part of the DP are updated conditionally to the subjects within each cluster, using the BDgraph .
• Missing values in the growth trajectories Y (s) , for s = 1, . . . , S or in the metabolite concentrations M are updated by sampling them from their full conditional distribution, which is readily obtainable through results on the conditional distribution of a multivariate Gaussian.