{'fraction_non_alphanumeric': 0.04000464630038332, 'fraction_numerical': 0.014124753165292136, 'mean_word_length': 4.542165572293035, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Joint Modeling of An Outcome Variable and Integrated Omic Datasets Using GLM-PO2PLS

In many studies of human diseases, multiple omic datasets are measured. Typically, these omic datasets are studied one by one with the disease, thus the relationship between omics are overlooked. Modeling the joint part of multiple omics and its association to the outcome disease will provide insights into the complex molecular base of the disease. In this article, we extend dimension reduction methods which model the joint part of omics to a novel method that jointly models an outcome variable with omics. We establish the model identifiability and develop EM algorithms to obtain maximum likelihood estimators of the parameters for normally and Bernoulli distributed outcomes. Test statistics are proposed to infer the association between the outcome and omics, and their asymptotic distributions are derived. Extensive simulation studies are conducted to evaluate the proposed model. The model is illustrated by a Down syndrome study where Down syndrome and two omics - methylation and glycomics - are jointly modeled.


Introduction
The biological mechanisms underlying human diseases are often complex. Diverse omic datasets represent various aspects of these mechanisms. Recent advances in highthroughput technologies have made it affordable to measure these omic levels for many studies. Typically, these datasets are studied one-by-one. A good example is the analysis of genomic data in more than 5700 GWAS conducted to identify the genetic risk variants associated with more than 3000 traits and human diseases [49,53]. Other examples include studies of methylation data to pinpoint differentially methylated regions of DNA as indicators of many diseases [32,39], and studies of glycomic data to gain insight into the role of post-translational modification of protein in disease pathways [45,46]. Though these studies on a single omic dataset provided biological insights of diseases from various aspects, they ignored the correlations among the omic levels. Analyzing multiple linked omic datasets jointly can bring further insights into the biological system underlying diseases. In this paper, we propose a new model for two omic datasets and an outcome variable, where the relationship of the omic datasets with the outcome is modeled via the joint parts of the omic datasets.
Our motivating dataset comes from a family-based case-control study of Down syndrome (DS). DS is the most frequent genomic aneuploidy with an incidence of approximately 1 in 700 live-newborn [33], caused by the trisomy of all or part of chromosome 21 (trisomy 21). Studies at the molecular level of DS have reported several alterations in methylation [4,8,15,16] and glycomics [6,9,15]. These alternations are mainly discovered by testing the mean difference of a single CpG site or glycan between the DS subjects and healthy controls. Furthermore, these studies were conducted on each omic level separately, overlooking the influence of methylation on glycosylation [52]. We will use our new model to jointly analyze DS and both omics, aiming to investigate whether the molecules involved in the relationship between methylation and glycomics are related to DS.
Generalized linear models (GLM) are flexible models which link linear predictors with an outcome variable via functions such as identity and logit [31]. However, omic datasets are often high-dimensional (p >> N ) and the features are highly correlated, which leads to (near) collinearity, hence GLM cannot directly be employed. A solution is to use its penalized variants such as ridge regression [19] and elastic net [57], which handle the dimensionality and correlation better. To incorporate information from more than one dataset, stacked sets of omic features can be modeled. However, such approaches do not model the relationship between omics, and hence cannot provide insight on the joint part. Furthermore, when the omic datasets are heterogeneous, regressing on a stacked dataset can lead to inferior performance than using only one of the available omic sources [37,47].
To model the joint part of two omic datasets, several dimension reduction methods have been developed, which map both datasets from the original high-dimensional spaces to low-dimensional spaces including a joint space which retains the relationship. Among these, PLS [10,54] simultaneously decomposes two datasets x and y into joint and residual subspaces. The joint (low-dimensional) subspace of one dataset represents the best approximation of x or y. The joint subspaces of PLS may contain omic-specific variation due to the presence of heterogeneity between the omics. To separate this omic-specific variation from the joint subspaces, two-way orthogonal partial least squares (O2PLS) [12,48] was proposed which decomposes each dataset into joint, data-specific, and residual subspaces. The data-specific subspaces in x and y capture variation unrelated to each other, and improve the estimates of the joint subspaces for the true relationship between x and y. A drawback of PLS and O2PLS is that they are algorithmic, and hence do not model the whole distribution. To enable statistical inference, likelihood-based probabilistic approaches such as supervised integrated factor analysis (SIFA) [25] and probabilistic PLS (PPLS) [13] were proposed. Recently, a general framework probabilistic O2PLS (PO2PLS) [14] was proposed which contains SIFA and PPLS as specific cases. In PO2PLS, the relationship between two omics is inferred by using a Wald-type test-statistic to test the hypothesis that the joint latent variables of the two omics are related. However, PO2PLS is unsupervised, i.e., the the outcome variable is not used in the joint dimension reduction. In this paper, we will propose a new model GLM-PO2PLS which extends PO2PLS by including an outcome variable in the model next to the omic datasets. The relationship between two omic data is modeled by joint and omic-specific latent variables to deal with possible heterogeneity. The joint latent variables are linked to an outcome variable by a generalized linear model. We develop EM algorithms to obtain maximum likelihood estimators of the parameters for normally and Bernoulli distributed outcomes. The relationship between the outcome variable and the omics and that between two omic datasets can be inferred. The code is available on GitHub (github.com/zhujiegu/GLM-PO2PLS).
The rest of the paper is organized as follows. In Section 2, the PO2PLS model is recapped, and the GLM-PO2PLS model is formulated. The EM algorithms to estimate its parameters are proposed. Also, two chi-square tests of the relationship between outcome and both omics are proposed. In Section 3, the performance of GLM-PO2PLS is studied for a range of simulation scenarios where the focus is on parameter estimation and outcome prediction performance. In Section 4, we apply GLM-PO2PLS to the motivating DS datasets. We conclude with a discussion.

Methods
The GLM-PO2PLS was developed based on PO2PLS model which has been described in detail elsewhere [14]. Briefly, let x and y be two random row-vectors of dimensions p and q, respectively. In PO2PLS, x and y are decomposed into joint (t and u of size r), specific (t ⊥ and u ⊥ of size r x resp. r y ) and residual (e and f of size p resp. q) parts. Heterogeneity between the joint parts is represented by an additional random vector h. The PO2PLS model is written as where W (p × r) and C (q × r) are the loading matrices for the joint spaces of x and y respectively and W ⊥ (p × r x ) and C ⊥ (q × r y ) are the loading matrices for the specific parts of x and y respectively. The r × r diagonal matrix B models the relationship between the joint components t and u. With regard to the random vectors, we assume that t, t ⊥ , u ⊥ , h are zero mean multivariate normally distributed, with diagonal covariance matrices Σ t , Σ t⊥ , Σ u⊥ , Σ h , respectively. Since u = tB + h, the covariance matrix of u is Σ u = B Σ t B + Σ h . The residual random vectors e, and f are independent normally distributed, with zero mean and respective diagonal covariance matrices , σ 2 e I p , and σ 2 f I q , where I p and I q are identity matrices of size p and q.

The GLM-PO2PLS model
GLM-PO2PLS jointly models an outcome variable z with two omic datasets x and y, where it is assumed that the effect of x and of y on z is solely through the joint parts of x and y.
Using the same notations as in the PO2PLS model, the GLM-PO2PLS model is given by with a 0 the intercept, a and b both row-vectors of size r and and η the link function which links the outcome z to the linear predictor a 0 + ta + ub . Note that the equations in the first row of (1) are identical to the PO2PLS model. Since the joint latent variables t and u are linked to x, y, and z, GLM-PO2PLS jointly models the outcome and two omics. Now, u is a linear function of t, namely u = tB + h. Hence, the model for z in (1) can equivalently be written in terms of t and the h (the part in u independent of t), i.e.
whereã = a + Bb andb = b. With this equivalent parametrization, instability due to near colinearity in the linear predictor of z is reduced.
In the remainder of the paper, we use the rightmost form in (2) and omit the tildes on a and b.

The GLM-PO2PLS model with a normally distributed outcome
In this subsection, we first consider a continuous outcome z. The details for a binary z is then given in the next subsection. As link function, we use the identity, η(v) = v. We assume that the outcome is centered and since t and h have zero-mean, the intercept a 0 can be omitted. We assume that the residual g (g = z − ta − ub ) is normally distributed, g ∼ N (0, σ 2 g ). Since (x, y, z) is linearly dependent on (t, u, t ⊥ , u ⊥ , e, f, h, g), it follows a multivariate normal distribution N (0, Σ θ ), with a covariance matrix given by e , σ 2 f , Σ h , σ 2 g } is the collection of GLM-PO2PLS model parameters.
Identifiability of GLM-PO2PLS Latent variable models are typically unidentifiable due to rotation indeterminacy of the loading components. In PO2PLS, identifiability up to sign has been shown under mild conditions [14]. Namely, the loading matrices are semi-orthogonal, i.e. W W = C C = I r , W ⊥ W ⊥ = I rx , and C ⊥ C ⊥ = I ry . Additionally, [W W ⊥ ] and [CC ⊥ ] do not have linearly dependent columns. Furthermore, the covariance matrices for the latent variables Σ t , Σ u , Σ t⊥ , Σ u⊥ are diagonal. Finally, the diagonal elements of B are positive and the diagonal elements of Σ t B are strictly decreasing.
We show that these conditions also guarantee the identifiability (up to sign) of the GLM-PO2PLS model. Theorem 2.1. Let (x, y, z) follow the GLM-PO2PLS model where z is normally distributed. Additionally, let the parameters satisfy the PO2PLS conditions as described above. It follows that the GLM-PO2PLS model parameters are identifiable up to a sign.

Maximum likelihood estimation
Since the GLM-PO2PLS model is a latent variable model and the likelihood factorizes in terms which can be maximized separately, we propose an EM algorithm [11] to obtain maximum likelihood estimates of the model parameters.
Suppose we observe the (x, y, z) for N subjects. Since we assume a multivariate normal distribution of (x, y, z) ∼ N (0, Σ θ ), the log-likelihood for one subject is given by Denote the complete data vector by (x, y, z, t, u, t ⊥ , u ⊥ ). For each current estimate θ , the EM algorithm considers the objective function Expectation step The conditional expectation of the complete data log likelihood can be decomposed into different terms, (4) In this equation, the conditioning on x, y, z and θ is dropped, to simplify notation. The individual conditional expectations depend on distinct sets of parameters, yielding separate optimization tasks. Compared to PO2PLS, the extra parameters in GLM- The rest of the terms are identical to the factorized densities in the original PO2PLS EM algorithm, we refer to the PO2PLS paper [14] for the expectation and maximization regarding these terms.
In the expectation step, Q {a,b,σ 2 g } is calculated as Here, the first and second conditional moments of the vector (t, u) given x, y, z and θ are involved. Since (x, y, z, t, u, t ⊥ , u ⊥ ) follows a multivariate normal distribution with zero mean and known covariance matrix, the conditional density f (t, u, t ⊥ , u ⊥ |x, y, z) can be calculated following Lemma 3 in [14]. The conditional moments involved in (5) can then be obtained from the mean and the covariance matrix of (t, u, t ⊥ , u ⊥ |x, y, z) (see the Supplementary material for details).
Maximization step In the maximization (M) step, each conditional expectation in (4) can be optimized separately. Here, we restrict to the description of the term involving the outcome, namely, maximize the Q {a,b,σ 2 g } as given in equation (5). Note that the coefficient vector (a, b) can be separately optimized from the residual parameter σ 2 g , as in the standard linear regressions. We first calculate the derivative with respect to (a, b) and set it to 0, yielding where the conditional moments are calculated in the E step. The maximization with respect to the parameter σ 2 g can then be performed similarly. Details are given in the supplementary material.

Statistical inference
The GLM-PO2PLS method allows for statistical inference on the relationship between the omic data and the outcome. This relationship is captured by the joint parts t and u, and given by the equation η(E[z]) = ta + ub in (1). Here, we propose two tests, one full test for the relationship between z and all the joint components together, and one component-wise test for the relationship between z and each pair of joint components.
For the full test, we consider the null hypothesis, For each component-wise test, we consider the null hypothesis of no relationship between z and the k-th pair of joint components, where a k and b k are the coefficients for t k and h k , respectively. Let α = (a, b) and α k = (a k , b k ). The full test statistic is given by where Π −1 α is the inverse of the covariance matrix ofα. And the pair-wise test statistic is given by: To calculate the (asymptotic) distribution of these test statistics, the asymptotic distribution of all parameters θ needs to be derived.
Asymptotic distribution Under certain regularity conditions, consistency of the estimator θ and its asymptotic distribution N (θ, Π θ ) follows from Shapiro's Proposition 4.2 (Shapiro 1986) applied to the GLM-PO2PLS model. Theorem 2.2. Letθ be the maximum likelihood estimator for θ from the GLM-PO2PLS model. When the sample size N approaches infinity, the distribution ofθ converges to a normal distribution, i.e.
Details and proofs are given in the supplement. In particular,α = (â,b) is asymptotically normally distributed. Therefore, the test statistics T f ull and T comp.wise follow a chi-square distribution with 2r resp. 2 degrees of freedom. An estimate of Π θ is obtained from the inverse observed Fisher information matrix. Let ψ i be an instance of observed data (x, y, z) and ζ i be the latent variables involved. In an EM algorithm, this matrix is given by [27], are the gradient and negative of the second derivative of the log complete likelihood of instance i, respectively, evaluated atθ.
To obtain Πα, the submatrix of I −1 (θ) corresponding toα (denote I −1 (α)) has to be calculated. However, inverting I(θ) is computationally infeasible, even for moderate dimensions. Under additional assumptions thatα andθ/α are asymptotically independent andσ 2 g is non-random, I −1 (α) can be calculated, and be used to approximate Πα. The details are given in supplementary materials.

The GLM-PO2PLS model with a binary outcome
For a binary outcome, we use a Bernoulli distribution for z and the logit link function The model is then given by Here, p(z) = Pr(z = 1|t, h) is the conditional probability of the random variable z being 1, given t and h. Note that the probability p(z) is logit-normally distributed, therefore the linear predictor logit(p(z)) follows a normal distribution N (a 0 , aΣ t a + bΣ h b ). The joint distribution (x, y, logit(p(z))) is multivariate normal with mean vector (0 p+q , a 0 ) and covariance matrix Σ θ in (3) excluding σ 2 g . The collection of parameters in the GLM-PO2PLS model with a binary outcome is Identifiability of GLM-PO2PLS with a binary outcome Theorem 2.1 also appears to hold for a binary z that follows a Bernoulli distribution, under the same conditions. The proof is similar. Specifically, let f (x, y, z|θ) = f (x, y, z|θ) be identical joint distributions under two sets of parameters θ andθ. Then f (x, y|θ) = f (x, y|θ), thus Σ x,y|θ = Σ x,y|θ holds regardless of the distribution of z. The conclusion follows that the parameters involved in PO2PLS model (i.e., {W, C, are identified up to sign. Now consider (x, y, logit(p(z))) which is multivariate normally distributed with mean vector (0 p+q , a 0 ) and covariance matrix Σ θ excluding σ 2 g (denote Σ θ/g 2 ).
Since the mapping f (x, y, z|θ) → f (x, y, logit(p(z))|θ) is one-to-one, it follows that f (x, y, logit(p(z))|θ) = f (x, y, logit(p(z))|θ). Necessarily, the means and covariance matrices of two identical multivariate normal distributions are equivalent, thus (0 p+q , a 0 ) = (0 p+q ,ã 0 ) and Σ θ/g 2 = Σθ /g 2 . It is clear that a 0 =ã 0 from the equivalence of the mean vectors. The identifiablity of a and b can be shown from the equivalence of covariance matrices analogously as in the proof of Theorem 2.1. This shows the identifiability of all the parameters in GLM-PO2PLS with a binary outcome.

EM algorithm for a binary outcome
For a Bernoulli distributed outcome, the log-likelihood of the observed data involves an integral of dimension 2r + r x + r y . Let ν = (t, u) and ξ = (t ⊥ , u ⊥ ), To estimate (8), numerical integration is needed. Note that given ν, the binary outcome z is independent of x, y and ξ, thus the conditional density f (x, y, z|ν, ξ) in (8) can be factorized as f (x, y, z|ν, ξ) = p(z|ν)f (x, y|ν, ξ). The factorization enables to integrate out the specific random vector ξ, hence reducing the dimension of the integral to 2r, Here, the probability mass function p(z|ν) is given by The probability density functions f (x|ν), f (y|ν), and f (ν) follow from the following multivariate normal distributions, where the covariance matrices involved are: Denote the partial complete data vector by (x, y, z, ν). For each current estimate θ , the EM algorithm for a binary outcome considers the objective function Expectation step based on numerical integration Analogously to (4), the conditional expectation in (10) can be decomposed to factors that depend on distinct sets of parameters, Here, the first conditional expectation Q {a0,a,b} has no closed form, To obtain an approximation of the multivariate integral, Gauss-Hermite quadrature can be used. For an integral of form ϕ(ν)p(z|ν)f (x, y|ν)f (ν)dν, where ϕ is any function, we approximate it with . . . is the Cholesky decomposition of Σ ν , and ν * m and w * m are nodes and weights of a M -point standard Gauss-Hermite quadrature rule, which can be found on Page 924 in [1]. The transformation from the standard quadrature nodes ν * m to ν m is to make the sampling range of the integrand in (12) more suitable based on the distribution of ν [26].
The other terms in (11) have explicit expressions in terms of the first and second conditional moments of the vector ν given x, y, z and θ (see for details in the Supplementary materials). Note that the conditional moments of ν are in forms of integrals as follows which can be numerically calculated with (12).
Maximization step based on gradient descent Maximizing Q {a0,a,b} requires iterations as its derivative with respect to β = (a 0 , a, b) has no analytical solutions.
To find an update of β in each EM iteration, we propose a one-step gradient descent strategy. The gradient of Q β is given by To guarantee the increase of Q β in each EM iteration, we search for a step size along the direction of the gradient using the backtracking rule (also known as the Armijo rule) [2]. It is performed by starting with an initial step size of s = 1 for movement along the gradient, and iteratively shrinking the step size (s ← 0.8 * s) until an increase of Q β exceeds the expected increase based on the local gradient. More precisely, we keep shrinking the step size until the following ascent condition is met: The maximization of the other conditional expectation terms in (11) can be found in the supplementary materials.

Simulation
We conduct a simulation study to evaluate the performance of GLM-PO2PLS. Both continuous outcome z c and binary outcome z b are investigated. The datasets are simulated following the GLM-PO2PLS model in (1), with the equations for the continuous and binary outcomes being z c = ta + hb + g, and z b ∼ Bernoulli((1 + exp{−(a 0 + ta + hb )}) −1 ).

Simulation settings
We consider combinations of small and large sample sizes (N = 100, 1000) with low and high dimensionalities (p = 100, 2000; q = 10, 25). The latent variables t, t ⊥ , u ⊥ are simulated from standard normal distribution, and u = tB + h following equation (1). Here, B is the identity matrix and the joint residual h in u that is independent of t determines the level of heterogeneity in the joint parts. To investigate the impact of heterogeneity levels, we vary the variance of h to account for 40% and 80% of the total variance in u. The residual terms e, f are generated from zero-mean normal distributions. In the low noise level scenario, we set the noise proportion in x and y to both 40%. In the high noise level scenario, we investigate the performance of GLM-PO2PLS when integrating a very noisy large dataset and a less noisy small dataset, by increasing the noise in x to 95% and decreasing the noise in y to 5%. The noise term g for the continuous outcome is generated from a zero-mean normal distribution, accounting for 20% of variation in z c . All the loading matrices are generated from standard normal distribution and then semi-orthogonalized. The coefficients a and b are set to 2 and 1, respectively. The number of joint and specific components is set to 1 for simplicity of the model and computational efficiency. For each setting, 500 replications are generated. The settings are summarized in Table 1. The metrics used to assess the performance are listed in Table 2. We first study the estimation accuracy of the coefficients a and b. The errors (â − a) and (b − b) are standardized by a and b to exclude the influence of the parameter scale. The performance of outcome prediction is assessed by root mean square error of prediction (RMSEP), defined as ( E [(ẑ c −z c ) 2 ]) 1 2 for continuous outcome z c , and ( E [(logit(p(z b ))− logit(p(z b ))) 2 ]) 1 2 for binary outcome z b . We compare the performance of GLM-PO2PLS with ridge regression fitted separately on x (denote ridge-x) and on y (denote ridge-y). The shrinkage hyper-parameter in ridge regressions is searched using a 10-fold crossvalidation for each fit. The prediction performance is evaluated on an independent test dataset of size 1000. The accuracy of loading estimation is measured by the inner product between the estimated and the true loading vectors. The performance of feature selection is measured by true positives rate (TPR) calculated as the proportion of true top 25% features among the estimated top 25% in x (i.e., the top 25% of features in x with the largest absolute loading values in GLM-PO2PLS, or with the largest absolute regression coefficients in ridge regression).

Results of simulation study
In Fig 1, results of the coefficient estimation in high-dimensional settings are depicted. Fig 1a shows that for the continuous outcome, overall, the scaled errors of bothâ andb were small. When the sample size was small and the noise was high, the scaled error (â − a)/a was mostly negative, suggesting that a was underestimated. For a large sample size, the estimators appeared to be unbiased. When the heterogeneity between the joint components was increased (from the left panel to the right), the joint residual h had larger variance relative to t and explained a larger proportion of z. Consequently, the estimation of the coefficient b (for h) became more stable, while the estimation of a (for t) became less stable. The results for a binary outcome are shown in Fig 1b. Under a small sample size, the parameter estimation was less stable than the continuous case (note that the scale of y-axis in subplot (a) and (b) are different). The long upper whiskers suggested that the coefficients were overestimated in a some simulation runs. For a large sample size, the scaled errors for all coefficients were close to 0 and stable. Overall, the results for low dimensions were similar, except that the estimation of b was less stable in low dimensions compared to that in high dimensions. Details are given in the supplementary material.
(a) Continuous outcome zc.     Lastly, we briefly present the key results for loading estimation and feature selection. Overall, the loading estimates were accurate for both continuous and binary outcomes, with most inner products between the estimated and the true loadings approaching the optimum. When the sample size was small and the noise level was high, the accuracy of loading estimation for x dropped. This was the same setting in whichâ was biased as is shown in Fig 1a. Regarding feature selection, the lowest median TPR of GLM-PO2PLS was 0.62 in the scenario with a small sample size, large noise proportion, and high heterogeneity. In the other scenarios, the median TPR was above 0.85. The details are given in the supplementary material.

Application to Down syndrome study
We apply the GLM-PO2PLS model to the Down syndrome dataset, aiming to investigate whether the relationship between methylation and glycomics is associated to DS, and select the relevant molecules involved in the relationship. Since Down syndrome is often considered as a model for aging [21], and both methylation and glycomics are associated with biological age [20,24], we expect the DS patients to be more similar to their mothers than siblings.

Data description
The Down syndrome study includes 29 families. Each family consists of one Down syndrome patient (DSP), one non-affected sibling (DSS), and their mother (DSM). The family-based design is used to control for genetic and environmental influences. Two DSS are missing. Thus, the total sample size N is equal to 85 . The ages of the DSPs range from 10 to 43, with a median of 24 years. The ages of the siblings are roughly matched with the DS patients, ranging from 14 to 52 years. The mothers have ages between 41 and 83, with a median of 57 years.
For each individual, the whole blood methylation was measured using Infinium Hu-manMethylation450 BeadChip (Infinium 450k). After quality control following steps described in [4], 450981 CpG sites were retained. Beta value was derived at each CpG site as the ratio of intensities between methylated and unmethylated alleles. White blood cell counts were estimated from the beta values and corrected for using R package 'Meffil' [29]. Age and sex were corrected for using multiple regression. The glycomic dataset consists of 10 plasma N-glycans measured using DNA sequencerassisted fluorophore-assisted carbohydrate electrophoresis (DSA-FACE) [6]. These glycans were logTA normalized [50] and corrected for age and sex.
We will fit a GLM-PO2PLS continuous model and a GLM-PO2PLS binary model to these data. We set methylation as x, glycomics as y, and the DS status as z. The direction from methylation to glycomics (x to y) was chosen based on previous research [52] that suggested the presence of an indirect influence of methylation on glycosylation.

Results of DS data analysis
For the GLM-PO2PLS continuous model, we used 3 joint and 1 specific component for each omic dataset based on the scree plots of the eigenvalues of x y, x x and y y.
We first present the results regarding the relationship between methylation and glycomics, which is represented by the first three equations of the GLM-PO2PLS in (1). The p-value for each pair of methylation and glycomics joint components was 0.0007, 0.03, and 0.20, respectively. Using a threshold of 0.05 for statistical significance, the first (t 1 for methylation and u 1 for glycomics) and second pair (t 2 and u 2 ) of joint components were significantly associated. Fig 3 shows the scores of the first two pairs of joint components. For both t 1 and u 1 , the DSPs were closer to the DSMs, than the DSS group, which was in line with our expectation. No noticeable patterns were observed in the second pair of joint components.  Table 3 shows the results regarding the relationship between the DS status and the omics. The significant test statistic T f ull suggests that DS was associated with the two omics. Component-wise, only the first pair was significant, with a p-value of 6.32 × 10 −5 . Table 3. Results of testing for no relationship between DS and joint components Since t 1 and u 1 were significantly associated with DS, we investigated the CpG sites and glycans in the first component pair. In the first methylation joint component, the 1000 CpG sites with the largest loading values were mapped to their respective target genes, yielding 493 genes. Next, gene ontology (GO) enrichment analysis [3] was performed on this gene set using the GSEA software [30,44]. The top three significant GO terms were listed in Table 4. Among these terms, the cell-cell signaling is a biological function of plasma glycans [23,51]. The cellular component of neuron projection and biological process of neurogenesis were shown to relate to DS [18,22,40,41]. We further searched the mapped geneset in the DisGeNET database [34] for human diseases. The significant diseases found were chronic myeloid leukemia (qvalue 0.0004), common acute lymphoblastic leukemia (q-value 0.045), and glioblastoma multiforme (q-value 0.045). Research has shown that children with Down syndrome have an increased risk for developing acute lymphoblastic leukemia [34]. For chronic myeloid leukemia and glioblastoma multiforme, we did not find evidence linking them with DS. We then checked the 7 genes with the highest GDA score regarding Down syndrome (i.e., with the most evidence of association with DS) in the DisGeNET database, and found the gene RCAN1 which relates to epigenetics was among our top genes mapped from methylation. It has been revealed that RCAN1 plays a critical upstream role in epigenetic regulation of adult neurogenesis [7], hence important in the pathogenesis of Down syndrome [55].
For the first glycomics joint component, the glycan H3N4F1 had the largest absolute loading value. According to the result of a previous study [6] on plasma glycans and DS, H3N4F1 was the top discriminators of DS subjects and siblings.
Next we fitted a GLM-PO2PLS binary model with 1 joint and 1 specific component for each omic dataset. We chose for 1 joint component based on the test results in the continuous model shown in Table 3. The relationship between the two omics was significant with a p-value of 0.022. The top 1000 CpG sites were identified and mapped to genes. The most significant GO terms of the geneset are shown in Table 4. The top two terms were related to membrane organelle, more specifically, Golgi apparatus, which is required for accurate glycosylation [56]. Terms related to DS, such as neurogenesis (q-value 2.33e-6), neuron differentiation (1.11e-5), and synapse (1.33e-5) were also significant. Regarding glycomics, the glycan with the largest absolute loading value was H3N4F1, which was also identified in the GLM-PO2PLS continuous model. The p-value of each annotation was derived by random sampling of the whole genome; the FDR q-value provides the false discovery rate (FDR) analog of the p-value after correcting for multiple hypothesis testing [5,43].
Although for the continuous and binary models the sets of top CpG sites appeared to be relevant to glycosylation and DS, there was little overlap between the top CpG sites. This might be explained by the different number of joint components specified for the two models. Therefore, we performed an additional analysis using a filtered dataset, which was obtained by subtracting the second and third joint components, and the specific components of the GLM-PO2PLS continuous model from the omic data. The estimated parameters of the fitted GLM-PO2PLS binary model were very similar to the ones in the continuous model. More specifically, the inner product of the joint loading vectors from the two models reached 0.99. The correlation between the corresponding joint components of the two models was also high (at 0.99). Note that an inner product of 1 means the loading vectors are the same. Regarding the interpretation of top CpG sites and glycans, we refer to the results in the continuous model, as the loadings were very similar, so as the top features identified.

Discussion
Motivated by the studies on the relationship among Down syndrome, methylation and glycomics, we developed a new statistical model GLM-PO2PLS, which simultaneously models the relationships among an outcome variable and two heterogeneous omic datasets. We studied in detail the models for normally and Benoulli distributed outcome variables. The identifiability of the model was established and EM algo-rithms were developed. For testing, we proposed two chi-square test statistics T f ull and T comp.wise and derived their asymptotic distributions.
Via a simulation study, we have shown that the model parameters were well estimated in various scenarios, and the outcome prediction performance of GLM-PO2PLS was robust against high noise and heterogeneity between omics. GLM-PO2PLS predicted the outcome better than ridge regressions, because it considers all the information in the data jointly, while ridge used each dataset separately. Another advantage of GLM-PO2PLS over ridge regression is that it can provide insights into the relationship between two omic datasets, on top of their relationship with the outcome.
The methylation and glycomics dataset were also analyzed by Bacalini et al. [4] and Borelli et al. [6] using single point approaches for association with Down syndrome, respectively. Since for both omics datasets association with aging have been found and Down syndrome is an aging model, it makes sense to analyze them jointly in this study. Concerning methylation, Bacalini et al. identified four categories of genes. Most of the genes in these categories were also in our obtained geneset: haematopoiesis (RUNX1, DLL1, EBF4, PRDM16), morphogenesis and development (HOXA2, HOXA4, HHIP, NCAM1), neuronal development (NAV1, EBF4, PRDM8, NCAM1), and regulation of chromatin structure (PRDM8, KDM2B). In total four genes mentioned in [4] were not in our gene list, namely, HOXA5, TET1, GABBR1, and HOXA6. It appeared that three out of these four genes rank just below our cut-off point of 1000, namely the CpG site with largest loading value in HOXA5, TET1, and GABBR1 ranked 1059, 1142, and 1535 out of 450K respectively. Concerning the fourth gene HOXA6, we performed univariate logistic regressions of the Down syndrome outcome on each of the 20 CpG sites located in the genetic region, and only identified one significant CpG site (p-value of 0.018). In comparison, the other selected genes from the HOXA family members have more significant CpG sites (such as HOXA2 with 22, HOXA4 with 16, the borderline HOXA5 with 13), and smaller p-values for the most associated CpG sites (HOXA2 0.0003, HOXA5 0.003). Furthermore, there is little evidence linking HOXA6 to the functions of glycans. Therefore, our proposed GLM-PO2PLS seems to better identify the CpG sites relevant to both DS and glycomics.
It is worth mentioning that we expect differences between our approach and the single-omic studies. The single-omic approaches did not consider the presence of correlation between CpG sites and glycomics when modeling the association of CpG sites and Down syndrome. Therefore, some methylation-specific genes that are unrelated with glycomics can rank lower in the joint components in our analysis. Furthermore, in GLM-PO2PLS, we focus on the joint part and the omic-specific parts are not linked to the outcome variable, and hence the top genes mapped from the methylation-specific components are not necessarily associated with the outcome DS. In this regard, an extension of our model which also considers the omic-specific parts in the linear predictor for the outcome variable can provide further insights into the disease from omic-specific aspects.
We have shown evidence for association between the mapped gene set and Down syndrome. Nonetheless, The dedicated "Down syndrome" set in the DisGeNET database was not enriched in our gene set. One reason could be that very few studies have been conducted on DS with methylation data. Furthermore, common diseases and cancers are usually more frequently studied, resulting in possible publication bias in the database. We searched the genes identified by both our study and [4] (namely, RUNX1, DLL1, EBF4, HOXA2, HOXA4, HHIP, NCAM1, NAV1, PRDM8, KDM2B) in the DisGeNET database, and found none of these genes has their highest association score (i.e., amount of evidence) with DS. For example, the RUNX1 gene had the highest association score of 0.8 with acute myeloid leukemia, and a score of only 0.1 with DS.
When estimating a GLM-PO2PLS binary model, we rely on numerical integration. The computational complexity of the numerical estimation is O(M 2r ), with M nodes per dimension. In practice this means that the binary model can only include 1 joint component. A computationally feasible solution is to include only one pair of the joint components in the linear predictor for the binary outcome. Such a model might be suited for our Down syndrome analysis where only one pair of joint components was associated to the outcome. However, the assumption that only one pair of joint components is related to the outcome might not apply to other studies. Therefore, a more efficient numerical integration strategy is needed. One strategy is to use adaptive quadrature. Although for a fixed number of nodes M , the adaptive quadrature is computationally more complex than its non-adaptive counterpart we used, the adaptive variant needs a smaller M to reach an equally precise approximation, thus can be more efficient [35,36]. Another strategy is to decompose the 2r-dimensional integration to r 2-dimensional integrations. This will reduce the computational complexity to O(r × M 2 ).
To calculate the p-values for the tests in (6) and (7), we derived the asymptotic normality of the estimator for the parameters of GLM-PO2PLS with a normally distributed outcome. Asymptotic normality was proved by showing that the mapping (denote τ ) from the parameter vector θ to the moment structure as well as the discrepancy function with respect to the moments satisfy certain regularity conditions [38]. For the GLM-PO2PLS model with a binary outcome, there is not an explicit mapping function τ , and it is difficult to parameterize the likelihood in terms of the moments. Therefore, while the p-values for the binary model can be calculated assuming the asymptotic normality holds, it is unclear whether they are correct. The derivation of asymptotic normality for the binary model is future work.
In this paper our aim was to use one model for all the data, and model the relationships between the omics simultaneously with their relationship with an outcome. While providing a holistic overview, for the binary outcome variable, the approach is computationally intensive. In our data analysis, it appeared that modeling the binary outcome as continuous provided similar information. Alternatively, a two-stage approach might be used. We recently proposed a two-stage PO2PLS approach [17], where we first constructed a few joint latent components that represent the two omics, then linked these latent components to the outcome variable using a linear regression model. In the implementation of two-stage PO2PLS to the DS dataset, the latent variables from the first stage were used as outcomes in several separate regression models in the second stage, thus the interpretation was different from a logistic regression model with DS as outcome. Alternatively, the latent variables can also be used as covariates in the second stage. However, the latent variables contain errors from the dimension reduction process. Ignoring these errors in the covariates can cause attenuated predicted probabilities in the logistic regression [42]. Therefore, to correctly model the outcome, the two-stage approach needs to be augmented with a measurement error model for the latent variables. Here more research is needed.
Several extensions of GLM-PO2PLS might be relevant. For an outcome variable from other members of the exponential family (e.g., Poisson, gamma, etc.), the corresponding EM algorithm can be obtained by modifying the EM algorithm for the binary outcome by replacing p(z|ν) in (9) with the corresponding conditional probability mass/density function. Regarding the relationship between omics and outcome, the omic-specific latent variables are not included in the linear predictor for the out-come variable in GLM-PO2PLS. As discussed above, linking the omic-specific parts to the outcome might provide further insights. Furthermore, since the omic-specific latent variable might also be predictive of the outcome, a model where all the latent variables are linked to the outcome can lead to improved outcome prediction performance in some studies. Extending GLM-PO2PLS to such a model will increase the computational complexity to O(M 2r+rx+ry ) for non-normal outcomes. Another direction is to generalize the model to incorporate more than two omic datasets jointly with an outcome. Such an extension would require to specify the directions of the relationships among more than two sets of variables. A workaround might be to model a common set of latent variables for all sets of variables [28]. For some studies, the directions of the relationships are clear (e.g., among genetics, methylation, and glycomics), and specifying the direction in the model and allowing the joint latent variables for each set of variables to differ can improve model performance. However, the computation will also be intensive for a binary outcome.
To conclude, GLM-PO2PLS is a promising method to model an outcome with two omic datasets and as a base for further extensions.

Data availability statement
The DNA methylation data used in this study are available at the NCBI Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/) under accession number GSE52588.