{"paper_id":"227c27de-6741-4cdb-a274-bab57bddd102","body_text":"1 \n \n \nACMTF-R: supervised multi -omics data integration uncover ing shared 1 \nand distinct outcome-associated variation 2 \nG.R. van der Ploeg 1*, F.T.G. White 1, R.R. Jakobsen 2, J.A. Westerhuis 1, A. Heintz -Buschart1, A.K. 3 \nSmilde1 4 \n1 Biosystems Data Analysis, Swammerdam Institute for Life Sciences, University of Amsterdam, 5 \nAmsterdam, The Netherlands 6 \n2 Department of Food Science, University of Copenhagen, Rolighedsvej 26, DK – 1958 Frederiksberg 7 \nC, Denmark 8 \n 9 \n*Author for correspondence: Geert Roelof van der Ploeg g.r.ploeg@uva.nl 10 \n 11 \nGeert Roelof van der Ploeg:  https://orcid.org/0009-0007-5204-3386 12 \nFred White:    https://orcid.org/0000-0003-0434-6028 13 \nRasmus Riemer Jakobsen  https://orcid.org/0000-0001-7416-1358 14 \nJohan Westerhuis:   https://orcid.org/0000-0002-6747-9779 15 \nAnna Heintz-Buschart:   https://orcid.org/0000-0002-9780-1933 16 \nAge Smilde:    https://orcid.org/0000-0002-3052-4644 17 \n 18 \nAbstract 19 \nThe rapid growth of high -dimensional biological data has necessitated advanced data fusion 20 \ntechniques to integrate and interpret complex multi -omics and longitudinal datasets. Advanced 21 \nCoupled Matrix and Tensor Factorization (ACMTF) has emerged as a powerful framework for 22 \nuncovering common, local, and distinct sources of variation across datasets. However, ACMTF lacks 23 \nthe ability to model variation linked to a dependent variable, limiting its applicability to studies 24 \ninvestigating biological phenotypes. N-way Partial Least Squares (NPLS) is a supervised method that 25 \nidentifies variation in relation to a dependent variable  but lacks the ability to identify common, 26 \nlocal and distinct sources of variation across multiple datasets. To bridge the gap between data 27 \nexploration and prediction, we introduce ACMTF-Regression (ACMTF-R), an extension of ACMTF 28 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n2 \n \n \nthat incorporates a regression term, allowing for the simultaneous decomposition of multi-way data 29 \nwhile explicitly capturing variation associated with an outcome variable. 30 \nWe present a detailed mathematical formulation of ACMTF-R, including its optimisation algorithm 31 \nand implementation. Through extensive simulations, we systematically evaluate its  ability to 32 \nrecover a small 𝒚-related component shared between multiple blocks, its robustness to noise, and 33 \nthe impact of the tuning parameter ( 𝜋) which controls the balance between data exploration and 34 \noutcome prediction. Our results demonstrate that ACMTF -R can robustly identify the 𝒚-related 35 \ncomponent, correctly identifying outcome-associated shared and distinct variation, distinguishing 36 \nit from existing approaches such as N-way Partial Least Squares and ACMTF. 37 \nTo validate its applicability in a real -world setting, we apply ACMTF -R to a multi -omics dataset 38 \nintegrating human milk microbiome, human milk metabolome, and infant faecal microbiome data, 39 \ninvestigating how maternal pre -pregnancy BMI  affects microbial and metabolic signatures. 40 \nACMTF-R successfully identifies novel mother-infant relationships associated with maternal pre-41 \npregnancy BMI, underscoring its utility in multi-omics research. Our findings establish ACMTF-R 42 \nas a versatile tool for multi-way data fusion, offering new insights into complex biological systems 43 \nby integrating common, local, and distinct variation in the context of a dependent variable. 44 \nKey Words 45 \nCoupled matrix and tensor factorization, multi-omics data integration, tensor decomposition 46 \nBackground 47 \nThe advent of high -throughput omics technologies has led to a rapid growth of complex, high -48 \ndimensional biological datasets, including genomics  [1–3], transcriptomics  [4], proteomics  [5], 49 \nmetabolomics [4], and microbiom e [6]. When measured on the same samples, t hese multi-omics 50 \ndatasets give researchers t he opportunity to obtain a systems -level understanding of biological 51 \nprocesses [5]. However, their integration and analysis pose significant challenges due to data 52 \nheterogeneity [5, 7], high dimensionality [4, 6], and complex interactions between different omics 53 \ndatasets [8, 9] . Advanced  data analytical methods are required to extract meaningful biological 54 \ninsights by distinguishing variation that is shared across datasets from variation that is dataset -55 \nspecific [10–13]. 56 \nA critical need in the analysis of multi-omics data is identifying common, local and distinct sources 57 \nof variation (Figure 1; [8–10, 14]). Common variation is shared across all datasets, local variation is 58 \nshared between a subset of  datasets, and d istinct variation is unique to a dataset. Identifying the 59 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n3 \n \n \nCommon, Local and Distinct (CLD) structure of a multi -omics dataset elucidates the complex 60 \nbiological processes that occur within the system. 61 \nAdvanced Coupled Matrix and Tensor Factorization (ACMTF) has emerged as a powerful multi -62 \nway data integration approach capable of uncovering common, local, and distinct sources of 63 \nvariation across datasets [10, 15–18]. However, ACMTF is unsupervised and therefore does not allow 64 \nthe identification of shared and distinct variation that relates to a dependent variable. N-way Partial 65 \nLeast Squares (N-PLS) [19] is a supervised method that identifies variation in relation to a dependent 66 \nvariable, but is limited to analysing one dataset at a time and thus cannot simultaneously identify 67 \ncommon, local, and distinct variation across multiple datasets. 68 \nTo bridge the gap between data exploration and prediction, we introduce Advanced Coupled Matrix 69 \nand Tensor Factorization Regression (ACMTF -R). ACMTF-R extends ACMTF by incorporating a 70 \nregression term into the objective function, enabling the extraction of common, local, and distinct 71 \nvariation while simultaneously predicting an outcome . By introducing a tuning parameter ( 𝜋), 72 \nwhich has been used successfully in other approaches such as Principal Covariates Regression 73 \n(PCovR) [20, 21] and Multiway Covariates Regression ( MCovR) [22], the solution can be steered 74 \ntowards explaining the data or towards predicting 𝒚. 75 \nIn this study, we first provide a detailed overview of the ACMTF framework, and the 76 \nmethodological development represented by ACMTF -R. We then present a simulation -based 77 \napproach to systematically evaluate the ability of ACMTF-R to capture a small, hidden, 𝒚-related 78 \ncomponent under different noise conditions  and tuning parameter values. Finally, we apply 79 \nACMTF-R to a real -world dataset, integrating human milk (HM) microbiome, HM metabolome, 80 \nand infant gut microbiome data to study how maternal pre-pregnancy body mass index (ppBMI) is 81 \nassociated with these data. Our findings highlight the potential of ACMTF-R as a versatile tool for 82 \nmulti-omics data integration, facilitating the discovery of biologically meaningful common, local 83 \nand distinct variation, while accommodating both exploratory and predictive research objectives. 84 \n  85 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n4 \n \n \n 86 \nFigure 1: common, local and distinct variation. Common variation 𝐶123 is shared between all 87 \ndatasets 𝑿(1), 𝑿(2) and  𝑿(3). Local variation is shared between some of the datasets. Hence 𝐿12 is 88 \nshared between 𝑿(1) and 𝑿(2), but not 𝑿(3). Similarly, 𝐿13 is shared between 𝑿(1) and 𝑿(3), but not 89 \n𝑿(2), and 𝐿23 is shared between 𝑿(2) and 𝑿(3), but not 𝑿(1). Distinct variation is unique to one 90 \ndataset: 𝐷1 for block 𝑿(1), 𝐷2 for block 𝑿(2), and 𝐷3 for block 𝑿(3). Together, they make up the 91 \nCommon, Local and Distinct (CLD) structure of the data. 92 \n  93 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n5 \n \n \nMethods 94 \nNotation and definitions 95 \nWe briefly define the mathematical notation that will be used throughout this paper. The notation 96 \nproposed by Kiers [23] is followed with some minor extensions for multi-omics. Scalars are denoted 97 \nby lower-case letters (e.g., 𝑎, 𝑏, 𝑐). Column vectors are denoted by boldface lowercase letters (e.g., 98 \n𝒂, 𝒃, 𝒄). Matrices are denoted by boldface capital letters (e.g., 𝑿, 𝒀, 𝒁). Three-way arrays are denoted 99 \nby underlined boldface capital letters (e.g., 𝑿, 𝒀, 𝒁). While 𝑿 can more generally be used to indicate 100 \nfour-way or higher-order arrays, this paper does not discuss such data. 101 \nThe characters 𝐼, 𝐽, and 𝐾 are reserved to indicate the first, second and third mode of an array. A 102 \ntwo-way array 𝑿 will be assumed to be of size 𝐼 × 𝐽, while a three-way array 𝑿 will be assumed to 103 \nbe of size 𝐼 × 𝐽 × 𝐾. Similarly, lowercase characters 𝑖, 𝑗, and 𝑘 will be used as indices for the first, 104 \nsecond, and third modes, respectively. Columns of a matrix are vectors with a subscript of their 105 \nindex (e.g., 𝒙1, 𝒙2, … , 𝒙𝐽 for some matrix 𝑿). Similarly, the i th-row and jth-column entry is denoted 106 \nas a scalar with subscripts of their indices (e.g., 𝑥𝑖𝑗). For three-way arrays, the entries per mode are 107 \ngiven as subscripts (e.g., 𝑥𝑖𝑗𝑘). 108 \nThe index 𝑝 = 1, … , 𝑃 is used to indicate the block number, which is denoted by a superscript (e.g., 109 \n𝑿(1), 𝑿(2), … , 𝑿(𝑃)) and the size is indicated likewise (e.g., the size of 𝑿(𝑝) is 𝐼 × 𝐽(𝑝) × 𝐾(𝑝)). In 110 \nACMTF and ACMTF-R 𝑿(1), 𝑿(2), … , 𝑿(𝑃) can be any mixture of two-way or three-way arrays, but 111 \nthis study focuses only on the case where all blocks are three-way arrays with a shared subject mode. 112 \nHence the superscript (𝑝) is not used to describe the size of the first mode 𝐼. 113 \nMathematically it is convenient to matricise a three-way array 𝑿, as this transformation facilitates 114 \nthe description of multi -way models using matrix notation [23–25]. There are three different 115 \nvariants of matricisation, depending on which mode of the three-way array is preserved as rows in 116 \nthe resulting matrix. Hence, first mode matrici sation yields a matrix 𝑿𝑎 of size 𝐼 × 𝐽𝐾. Similarly, 117 \nsecond mode matricization produces a matrix 𝑿𝑏 of size 𝐽 × 𝐼𝐾, while third mode matricization 118 \ngenerates a matrix 𝑿𝑐 of size 𝐾 × 𝐼𝐽. 119 \nThe vectorisation of an 𝐼 × 𝐽 matrix 𝑿 is denoted vec(𝑿), such that vec(𝑿) = [𝒙1\n𝑇 𝒙2\n𝑇 … 𝒙𝐽\n𝑇]\n𝑇\n. The 120 \nHadamard product [26] of two equally sized arrays 𝑿 and 𝒀 is denoted 𝑿 ∗ 𝒀, such that (𝑿 ∗ 𝒀)𝑖𝑗 =121 \n𝑥𝑖𝑗𝑦𝑖𝑗. The Kronecker product, denoted ⨂, of two matrices 𝑿 with size 𝐼 × 𝐽 and 𝒀 with size 𝐾 × 𝐿 122 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n6 \n \n \nproduces a block matrix containing all pairwise products of their entries following (𝑿 ⨂ 𝒀)𝑖𝑘,𝑗𝑙 =123 \n𝑥𝑖𝑗𝑦𝑘𝑙. Thus, we have 124 \n𝑿⨂𝒀 =  [\n𝑥11𝒀 ⋯ 𝑥1𝐽𝒀\n⋮ ⋱ ⋮\n𝑥𝐼1𝒀 ⋯ 𝑥𝐼𝐽𝒀\n] (1) 125 \nwhere the outcome of 𝑿⨂𝒀 is size 𝐼𝐾 × 𝐽𝐿. 126 \nThe Khatri-Rao product, denoted ⨀, of two matrices 𝑿 of size 𝐼 × 𝑅 and 𝒀 of size 𝐽 × 𝑅 is defined 127 \nas the column-wise Kronecker product [27, 28] 128 \n𝑿⨀𝒀 = [\n| | ⋯ |\n𝒙1⨂𝒚1  𝒙2⨂𝒚2 ⋯ 𝒙𝑅⨂𝒚𝑅\n| | ⋯ |\n] = [\n| | ⋯ |\nvec(𝒚1𝒙1\n𝑇) vec(𝒚2𝒙2\n𝑇) ⋯ vec(𝒚𝑅𝒙𝑅\n𝑇)\n| | ⋯ |\n] 129 \nwhere the outcome of 𝑿⨀𝒀 has size 𝐼𝐽 × 𝑅. 130 \nThe Frobenius (or Euclidean) norm of a three-way array 𝑿 is defined as 131 \n‖𝑿‖ = √∑ ∑ ∑ 𝑥𝑖𝑗𝑘\n2\n𝐾\n𝑘=1\n𝐽\n𝑗=1\n𝐼\n𝑖=1\n (2) 132 \nand the notation ‖ ∙ ‖ will also be used to refer to the Frobenius norm for matrices and vectors, 133 \nrespectively [24]. The Moore-Penrose inverse [29] of 𝑿 is denoted 𝑿+. 134 \nIn this paper we follow the canonical polyadic form of defining multi-way arrays as the sum of rank-135 \none arrays with normalized factor matrices and a scalar 𝜆 to capture the magnitude per term  [15, 136 \n30–32]. The notation ⟦𝝀𝑝 ; 𝑨, 𝑩(𝑝), 𝑪(𝑝)⟧ is used to describe the model of the three -way array 𝑿(𝑝), 137 \ndefined elementwise as 138 \n𝑥𝑖𝑗𝑘\n(𝑝) = ∑ 𝜆𝑝𝑓𝑎𝑖𝑓𝑏𝑗𝑓\n(𝑝)𝑐𝑘𝑓\n(𝑝) + 𝑒𝑖𝑗𝑘\n(𝑝)\n𝐹\n𝑓=1\n (3) 139 \nwhere columns of the loading matrices 𝑨, 𝑩(𝑝), 𝑪(𝑝) are norm 1 and the column vector 𝝀𝑝 is used to 140 \nmodify the size of each identified component 𝑓 = 1, … , 𝐹 to norm 𝜆𝑝𝑓. For a two -way array, 141 \nequation (3) reduces to ⟦𝑨, 𝑩⟧ = 𝑨𝑩𝑇, assuming that 𝝀 = 1 with size 1 × 𝐹 and that the singular 142 \nvalues are multiplied into 𝑨 to stay in line with common practice. 143 \n 144 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n7 \n \n \nAdvanced Coupled Matrix and Tensor Factorization (ACMTF) 145 \nAdvanced Coupled Matrix and Tensor Factorization is a joint factor isation approach that is 146 \napplicable to mixtures of matrices and tensors [10, 15, 16, 33]. The method seeks a decomposition 147 \nof the datasets 𝑿(1), 𝑿(2), … , 𝑿(𝑃) as 𝑿̂(𝑝) = ⟦𝝀𝑝; 𝑨, 𝑩(𝑝), 𝑪(𝑝)⟧ such that loading matrices 148 \ncorresponding to shared mode s are equal between the blocks. In this paper, we will assume that 149 \nonly the first mode is shared between all blocks. In that case, the matrix containing the first mode 150 \nloadings 𝑨 will be the same for every block 𝑝 and the superscript can be left out. 151 \nThe loss function for ACMTF is defined as 152 \n𝑓(α, β, ε, 𝚲, 𝑨, 𝑩(1), 𝑪(1), … , 𝑩(𝑃), 𝑪(𝑃)) =  ∑‖𝑿(𝑝) − 𝑿̂(𝑝)‖\n2\n𝑃\n𝑝=1\n+ 𝛽 ∑ ∑ √𝜆𝑝𝑓\n2 + 𝜖\n𝐹\n𝑓=1\n𝑃\n𝑝=1\n+ 𝛼 ∑ ∑(‖𝒂𝑓‖− 1)\n2\n𝐹\n𝑓=1\n𝑃\n𝑝=1\n+ 𝛼 ∑ ∑ (‖𝒃𝑓\n(𝑝)‖ − 1)\n2\n𝐹\n𝑓=1\n𝑃\n𝑝=1\n+ 𝛼 ∑ ∑ (‖𝒄𝑓\n(𝑝)‖ − 1)\n2\n𝐹\n𝑓=1\n𝑃\n𝑝=1\n (4)\n 153 \nwhere the 𝑃 × 𝐹 matrix 𝚲 = [𝝀1 𝝀2 … 𝝀𝑃]𝑇 scales the factors per component and per data block, 154 \n𝛼 > 0 is a penalty setting for norm-1 components, 𝛽 > 0 is a sparsity penalty setting for 𝚲, 𝐹 is the 155 \nnumber of components, and 𝑝 is the block index. Here the 1-norm loss term for the elements in 𝚲 156 \nis replaced with a differentiable approximation 𝛽‖𝜆𝑝𝑓‖1 ≈  𝛽√𝜆𝑝𝑓\n2 + 𝜖 for a sufficiently small 𝜖 >157 \n0 [33, 34].  The originally suggested default settings are 𝛼 = 1, 𝛽 = 10−3 and 𝜖 = 10−8 [10]. The 158 \ngradient of the loss function is reported in Supplementary Methods. 159 \nThe ACMTF framework of putting the norm of a component into 𝚲 and constraining the loading 160 \nvectors to become norm one, allows the 𝚲 matrix to encode the Common, Local, and Distinct (CLD) 161 \nstructure of the data (Figure 1). For example, the CLD-structure of a three-block case 𝑿(1), 𝑿(2), 𝑿(3) 162 \ncan be encoded in 𝚲 through 163 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n8 \n \n \n𝚲 =\n𝑿(1)\n𝑿(2)\n𝑿(3)\n𝐶123 𝐿12 𝐿13 𝐿23 𝐷1 𝐷2 𝐷3\n[\n𝜆11 𝜆12 λ13 0 𝜆15 0 0\n𝜆21 𝜆22 0 𝜆24 0 𝜆26 0\n𝜆31 0 𝜆33 𝜆34 0 0 𝜆37\n] 164 \nwhere 𝜆𝑝𝑓 > 0 corresponds to data block 𝑝 contributing to component 𝑓, 0 corresponds to a data 165 \nblock not contributing to a component, 𝐶123 indicates a common component shared between all 166 \nblocks, 𝐿12, 𝐿13 and 𝐿23 indicate local components shared between some blocks, and 𝐷1, 𝐷2 and 𝐷3 167 \nindicate distinct components unique to one block. Sparsity is then achieved by applying an 𝐿1 168 \npenalty on the 𝚲-matrix. While this example uses 7 components to showcase the entire CLD -169 \nstructure, multiple components may be needed to fully describe a specific type of variation in real 170 \ndata (e.g., two components for 𝐶123) and some types of variation may be absent altogether. 171 \nAdvanced Coupled Matrix and Tensor Factorization Regression (ACMTF-R) 172 \nWe present Advanced Coupled Matrix and Tensor Factorization Regression (ACMTF-R) to describe 173 \ncommon, local, and distinct variation of interest in the data blocks in the context of a dependent 174 \nvariable 𝒚 by adding a regression term and a tuning parameter 𝜋 to the loss function.  The 175 \nintroduction of a tuning parameter has been used successfully in other approaches such as Principal 176 \nCovariates Regression (PCovR) [20, 21] and Multiway Covariates Regression (MCovR) [22] to steer 177 \nthe solution towards explaining the data blocks or towards predicting 𝒚. This gives the user more 178 \ninformation about how the addition of 𝒚 affects the model of the data. We only define ACMTF-R 179 \nin the case where the first mode is shared across all data blocks and  𝒚. Hence 𝒚 is a column vector 180 \nof size 𝐼 × 1. The loss function of ACMTF-R is then defined as 181 \n𝑓(α, β, ε, π, 𝚲, 𝑨, 𝑩(1), 𝑪(1), … , 𝑩(𝑃), 𝑪(𝑃)) =  π ∑‖𝑿(𝑝) − 𝑿̂(𝑝)‖\n2\n𝑃\n𝑝=1\n+ (1 − 𝜋)‖𝒚 − 𝑨𝝆‖2\n+ 𝛽 ∑ ∑ √𝜆𝑝𝑓\n2 + 𝜖\n𝐹\n𝑓=1\n𝑃\n𝑝=1\n+ 𝛼 ∑ ∑(‖𝒂𝑓‖− 1)\n2\n𝐹\n𝑓=1\n𝑃\n𝑝=1\n+ 𝛼 ∑ ∑ (‖𝒃𝑓\n(𝑝)‖ − 1)\n2\n𝐹\n𝑓=1\n𝑃\n𝑝=1\n+ 𝛼 ∑ ∑ (‖𝒄𝑓\n(𝑝)‖ − 1)\n2\n𝐹\n𝑓=1\n𝑃\n𝑝=1\n (6)\n                                                                (6) 182 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n9 \n \n \nwhere 𝜋 is a tuning parameter used to focus the model on explaining the data blocks 𝑿(1), 𝑿(2), …, 183 \n𝑿(𝑃) versus predicting 𝒚, and 𝝆 is the 𝐹 × 1 vector of regression coefficients of the factor matrix of 184 \nthe first mode 𝑨 onto 𝒚. Here 𝝆 is given when 𝑨 is defined, through 185 \n𝝆 = 𝑨+𝒚. (7) 186 \nThe gradient of the loss function is reported in Supplementary Methods. 187 \nPrediction of 𝒚 using a new sample 188 \nA fitted ACMTF-R model can be used to predict 𝑦𝑛𝑒𝑤 for a new sample using a similar procedure as 189 \nLatent Root Regression [35, 36]. First, this requires a brief description of obtaining the subject mode 190 \nloadings for a new sample in CANDECOMP/PARAFAC (CP) [30, 32, 37, 38] . Subsequently, this 191 \nresult is extended to ACMTF-R, after which the prediction step is described. 192 \nIn the CP case, finding the subject mode loadings 𝒂𝑛𝑒𝑤 with size 𝐹 × 1 for a sample 𝑿𝑛𝑒𝑤 with size 193 \n𝐽 × 𝐾 requires solving the problem 194 \narg min\n𝒂𝑛𝑒𝑤\n‖𝑿𝑛𝑒𝑤 − 𝑩diag(𝒂𝑛𝑒𝑤)𝑪𝑇 ‖𝐹\n2 = arg min\n𝒂\n‖vec(𝑿𝑛𝑒𝑤) − (𝑪⨀𝑩)𝒂𝑛𝑒𝑤‖𝐹\n2  (9) 195 \nwhere 𝑩 and 𝑪 are the feature and time mode loadings of a previously fitted model  [32]. This 196 \nproblem statement comes down to the least squares solution 197 \n𝒂𝑛𝑒𝑤 = (𝑪⨀𝑩)+vec(𝑿𝑛𝑒𝑤) = 𝒁+vec(𝑿𝑛𝑒𝑤) (10) 198 \nwhere the block matrix 𝒁 = 𝑪⨀𝑩 with size 𝐽𝐾 × 𝐹 is used for convenience and where 𝒁+ is the 199 \nMoore-Penrose inverse of 𝒁. 200 \nIn the case of ACMTF-R, the shared subject mode loadings across the blocks need to be identified. 201 \nThis is done by finding the block matrix 𝒁(𝑝) with size 𝐾(𝑝)𝐽(𝑝) × 𝐹 through 202 \n𝒁(𝑝) = 𝝀𝑝𝑇⨀𝑪(𝑝)⨀𝑩(𝑝) (11) 203 \nwhich can be concatenated across all blocks into one matrix for all data blocks 204 \n1, … , 𝑃 simultaneously as 205 \nZ = (\n𝒁(1)\n𝒁(2)\n⋮\n𝒁(𝑃)\n)  (12) 206 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n10 \n \n \nwhere 𝒁 is of size (𝐾(1)𝐽(1) + 𝐾(2)𝐽(2) + ⋯ + 𝐾(𝑃)𝐽(𝑃)) × 𝐹. The shared subject mode loadings 𝒂𝑛𝑒𝑤 207 \nwith size 𝐹 × 1 for a new sample 𝑿𝑛𝑒𝑤 can then be found by reusing Equation 10. 208 \nFinally, the prediction of 𝑦𝑛𝑒𝑤 is done using the ACMTF-R regression coefficients 𝝆 through 209 \n𝑦𝑛𝑒𝑤 = 𝝆𝑻𝒂𝑛𝑒𝑤. (12) 210 \nWhen multiple new samples are obtained, this procedure is performed per sample. The procedure 211 \nis performed by the npred() function in the supplied R package. 212 \nImplementation and stopping criteria 213 \nThe decomposition of ACMTF and ACMTF -R is achieved through an all -at-once optimi sation 214 \nalgorithm, originally developed in [15] for MATLAB, now available in the CMTFtoolbox package 215 \nfor R on CRAN. The all-at-once optimisation is achieved by defining the loss and gradient function 216 \nand subsequently using the nonlinear conjugate gradient (NCG) method with Hestenes -Stiefel 217 \nupdates and the Moré -Thuente line search as originally suggested [15, 18] . For the provided R 218 \npackage, this was implemented using the mize package (v0.2.4, [39]). The provided package also 219 \nsupports the Limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm (L-BFGS) to speed up 220 \nthe line search as an experimental feature, but this setting is not used for any results in this paper 221 \n[40–42]. 222 \nACMTF-R uses the same default values as ACMTF for the minimum value of the loss function 223 \n(default 10−10) and the minimum change in loss function between two function evaluations (default 224 \n10−10) as stopping criteria [33]. The all-at-once optimization implementation through the mize 225 \npackage also defines three other stopping criteria: the maximum number of iterations  to calculate 226 \n(default 10 000), the maximum number of function evaluation s that are allowed  (default 10 000), 227 \nthe minimum l2-norm of the gradient vector (default 10−10), and the absolute value of the size of 228 \nthe parameter update (default 10−10). Fitting an ACMTF -R model for any of the datasets in this 229 \npaper takes less than a minute of computational time on an average computer (Microsoft Windows 230 \n11 Home v10.0.22631 with an 12 th Gen Intel® Core™ i5-12400F running 2500 MHz with 6 cores 231 \nor 12 logical processors and 16 Gb RAM). Multi -core parallelisation has been implemented as part 232 \nof the provided R package to allow many randomly initiali sed models to be fitted simultaneously.  233 \nThis is needed to efficiently find the appropriate number of components through cross -validation, 234 \nas well as to find the global minimum for a given number of components. 235 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n11 \n \n \nSimulation approach 236 \nWe demonstrate the effect of the ACMTF-R model parameters using various simulations (Figure 2). 237 \nWe follow the guidelines provided in previous work [10]. We generate factor matrices 𝑨 ∈ ℝ𝐼×𝐹, 238 \n𝑩(1) ∈ ℝ𝐽(1)×𝐹, 𝑪(1) ∈ ℝ𝐾(1)×𝐹, 𝑩(2) ∈ ℝ𝐽(2)×𝐹, 𝑪(2) ∈ ℝ𝐾(2)×𝐹, 𝑩(3) ∈ ℝ𝐽(3)×𝐹 and 𝑪(3) ∈ ℝ𝐾(3)×𝐹 239 \nwith random entries drawn from the standard normal distribution 𝑁~(0,1). The columns of the 240 \nfactor matrices are normalized to unit norm. We set 𝐼 = 100, 𝐽(1) = 𝐽(2) = 𝐽(3) = 200, and 𝐾(1) =241 \n𝐾(2) = 𝐾(3) = 5 to resemble a typical pre-processed multi-omics dataset (for example, a mixture of 242 \nmicrobiome and metabolomics data). The factor matrices are used to create three third-order tensors 243 \n𝑿(1) ∈ ℝ100×200×5, 𝑿(2) ∈ ℝ100×200×5, 𝑿(3) ∈ ℝ100×200×5 through 244 \n𝑿(1) = ⟦𝝀1; 𝑨, 𝑩(1), 𝑪(1)⟧, 245 \n𝑿(2) = ⟦𝝀2; 𝑨, 𝑩(2), 𝑪(2)⟧, 246 \nand 247 \n𝑿(3) = ⟦𝝀3; 𝑨, 𝑩(3), 𝑪(3)⟧. 248 \nwhere we encode the required CLD -structure into 𝚲 consisting of eight components total: two 249 \ncommon components shared across all blocks, three local components, and a distinct component for 250 \neach block. This corresponds to the 𝚲-matrix 251 \n𝚲 = [\n1 1 1 𝛿 0 1 0 0\n1 1 1 0 1 0 1 0\n1 1 0 𝛿 1 0 0 1\n] 252 \nwhere all components are norm 1,  except for the second local component whose size 𝛿 varies per 253 \nsimulation. 254 \nNext, randomly distributed noise is added to the third-order tensors through 255 \n𝑿𝑛𝑜𝑖𝑠𝑦\n(𝑝) = 𝑿(𝑝) + 𝜂𝑋\n‖𝑿(𝑝)‖\n‖𝑵(𝑝)‖𝑵(𝑝) (13) 256 \nwhere 𝜂𝑋 indicates the noise level, which is the same across all blocks. 257 \nWe generate 𝒚 equal to the subject mode loadings of the second local component 258 \n𝒚 = 𝒂𝐿12 259 \nand some noise is added through 260 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n12 \n \n \n𝒚𝑛𝑜𝑖𝑠𝑦 = 𝒚 + 𝜂𝑦\n‖𝒚‖\n‖𝒏𝑦‖𝒏𝒚 (14) 261 \nwhere 𝜂𝑦 indicates the noise level on 𝒚. The outcome 𝒚𝑛𝑜𝑖𝑠𝑦 is then normalized to norm 1. 262 \n 263 \nFigure 2: overview of the simulation setup of the study. All simulations contain a CLD-structure of 264 \n2 global components shared between blocks 𝑿(1), 𝑿(2), 𝑿(3), 3 local components , and  a distinct 265 \ncomponent unique to each block. This yields 8 components total for ACMTF -R to identify. All 266 \ncomponents are norm 1, except for the second local component whose norm 𝛿 varies per simulation. 267 \n𝒚 is equal to the subject mode loadings of the second local component . Noise and terms 268 \ncorresponding to components that have no contribution to a block have been removed for visual 269 \nclarity. Subject mode loadings (black lines) are equal between the data blocks per component but 270 \nare different between components. 271 \n  272 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n13 \n \n \nSimulation 1: ACMTF-R can detect a small, hidden component by leveraging 𝒚 273 \nIn the first simulation, the ability of ACMTF-R to detect a small, hidden component by leveraging 274 \n𝒚 was examined. This was done using the three -block simulation with 𝛿 = 0.03, 40% noise on 275 \n𝑿(1), 𝑿(2), 𝑿(3) (𝜂𝑋 = 0.40), and 5% noise on 𝒚 (𝜂𝑦 = 0.05). The data were then block scaled to 276 \nFrobenius norm 1. One hundred randomly initialised eight -component ACMTF models were 277 \ncompared with one hundred randomly initialised eight -component ACMTF -R models using 278 \ndifferent values of 𝜋 ∈ {0.1, 0.2, … ,0.9}, keeping all other parameter and convergence settings at 279 \ntheir default values. Subsequently, factor recovery was assessed using Factor Match Score and 280 \nTucker Congruence Coefficient. The line search algorithm in most ACMTF -R models was 281 \nterminated due to the absolute change for the size of the parameter update stopping criterion, while 282 \nACMTF was terminated only due to the minimum value of the loss function being reached  283 \n(Supplementary Figure 1). 284 \nSimulation 2: the benefit of supervision using 𝒚 versus the size of the hidden component 285 \nIn the second simulation, the improvement of ACMTF -R compared to ACMTF was assessed by 286 \nchanging the size of the small, hidden component. This was done using the three-block simulation 287 \nwith 288 \n𝛿 ∈ {0.01, 0.02, 0.03, 0.04, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5}. 289 \nThe noise on 𝑿(1), 𝑿(2), 𝑿(3) was set at 40% ( 𝜂𝑋 = 0.40), and the noise on 𝒚 was 5% (𝜂𝑦 = 0.05). 290 \nThe data were then block scaled to Frobenius norm 1.  One hundred randomly initialised eight -291 \ncomponent ACMTF models were compared with one hundred randomly initialised eight -292 \ncomponent ACMTF-R models using different values of 𝜋 ∈ {0.1, 0.2, … ,0.9} for each case, keeping 293 \nall other parameter and convergence settings at their default values. Subsequently, factor recovery 294 \nwas assessed using Factor Match Score and Tucker Congruence Coefficient. The line search 295 \nalgorithm in most ACMTF-R models was terminated due to the absolute change for the size of the 296 \nparameter update stopping criterion, while most ACMTF models were terminated only due to the 297 \nminimum value of the loss function being reached (Supplementary Figure 2). 298 \n 299 \n 300 \n 301 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n14 \n \n \nSimulation 3: the benefit of supervision using 𝒚 versus noise on X 302 \nIn the third simulation, the improvement of ACMTF -R compared to ACMTF was assessed by 303 \nchanging the amount of noise on 𝑿(1), 𝑿(2), 𝑿(3). This was done using the three -block simulation 304 \nwith 𝛿 = 0.03 and a variable amount of noise on 𝑿(1), 𝑿(2), 𝑿(3), 305 \n𝜂𝑋 ∈ {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 2, 3, 4, 5, 10, 50, 100}, 306 \nand the noise on 𝒚 was 5% (𝜂𝑦 = 0.05). The data were then block scaled to Frobenius norm 1. One 307 \nhundred randomly initialised eight-component ACMTF models were compared with one hundred 308 \nrandomly initialised eight -component ACMTF -R models using different values of 𝜋 ∈309 \n{0.1, 0.2, … ,0.9} for each case, keeping all other parameter and convergence settings at their default 310 \nvalues. Subsequently, factor recovery was assessed using Factor Match Score and Tucker 311 \nCongruence Coefficient. For most noise levels (𝜂𝑋 < 2), the line search algorithm in most ACMTF-312 \nR models was terminated due to the absolute change for the size of the parameter update stopping 313 \ncriterion, while most ACMTF models were terminated only due to the minimum value of the loss 314 \nfunction being reached ( Supplementary Figure 3 ). At higher noise levels ( 𝜂𝑋 > 2), most ACMTF 315 \nmodels were terminated due to the relative change in the loss function being reached. 316 \nSimulation 4: the benefit of supervision using 𝒚 versus noise on 𝒚 317 \nIn the fourth and final simulation, the improvement of ACMTF-R compared to ACMTF was assessed 318 \nby changing the amount of noise on 𝒚. This is done using the three-block simulation with 𝛿 = 0.03, 319 \nwith 40% noise on 𝑿(1), 𝑿(2), 𝑿(3) (𝜂𝑋 = 0.40), and a variable amount of noise on 𝒚: 320 \n𝜂𝑦 ∈ {0.0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 2, 3, 4, 5, 10, 50, 100}. 321 \nThe data were then block scaled to Frobenius norm 1. One hundred randomly initialised eight -322 \ncomponent ACMTF models were compared with one hundred randomly initialised eight -323 \ncomponent ACMTF-R models using different values of 𝜋 ∈ {0.1, 0.2, … ,0.9} for each case, keeping 324 \nall other parameter and convergence settings at their default values. Subsequently, factor recovery 325 \nwas assessed using Factor Match Score and Tucker Congruence Coefficient. The line search 326 \nalgorithm in most ACMTF-R models was terminated due to the absolute change for the size of the 327 \nparameter update stopping criterion, while ACMTF was terminated only due to the minimum value 328 \nof the loss function being reached (Supplementary Figure 4). 329 \n 330 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n15 \n \n \nSimulation evaluation – Tucker Congruence Coefficient 331 \nIn this paper, the recovered factors are compared against the input factors of the simulation using 332 \nthe absolute Tucker Congruence Coefficient (TCC) [43, 44] through 333 \n|𝜙(𝒙, 𝒚)| = || ∑𝑥𝑖𝑖 𝑦𝑖\n√∑𝑥𝑖\n2𝑖 ∑𝑦𝑖\n2𝑖\n|| (15) 334 \nwhere 𝒙 and 𝒚 are the recovered and true factors, respectively. The TCC can be interpreted as the 335 \ncosine of the angle between 𝒙 and 𝒚. A TCC value in the range of 0.85 -0.94 suggests that the 336 \nrecovered and input factors are fairly similar, whereas a value of 0.95 -1 means that they are 337 \nessentially equal [45]. Although TCC is insensitive to scalar multiplication of the vectors (that is, 338 \n𝜙(𝒙, 𝒚) = 𝜙(𝑎𝒙, 𝑏𝒚) for any non-zero scalars 𝑎 and 𝑏), it is sensitive to sign flips between 𝒙 and 𝒚 339 \n[45]. For this reason, we use the absolute value instead. 340 \nSimulation evaluation – Factor Match Score 341 \nThe Factor Match Score (FMS) is an adaptation of the Tucker Congruence Coefficient and indicates 342 \nthe similarity of two multi-way models across all components. For model evaluation, FMS was used 343 \nto compare the decomposition against the input factors. This depends on the CLD -type of the 344 \nrelevant factor. For the hidden local component 𝐿13 in the simulations, 𝐹𝑀𝑆𝐿13 is defined as 345 \n𝐹𝑀𝑆𝐿13 = |𝒂𝐿13\n𝑇 𝒂̃𝐿13|\n‖𝒂𝐿13‖‖𝒂̃𝐿13‖\n|[𝒃𝐿13\n(1) ]\n𝑇\n𝒃̃𝐿13\n(1) |\n‖𝒃𝐿13\n(1) ‖ ‖𝒃̃𝐿13\n(1) ‖\n|[𝒄𝐿13\n(1) ]\n𝑇\n𝒄̃𝐿13\n(1) |\n‖𝒄𝐿13\n(1) ‖ ‖𝒄̃𝐿13\n(1) ‖\n \n|[𝒃𝐿13\n(3) ]\n𝑇\n𝒃̃𝐿13\n(3) |\n‖𝒃𝐿13\n(3) ‖ ‖𝒃̃𝐿13\n(3) ‖\n|[𝒄𝐿13\n(3) ]\n𝑇\n𝒄̃𝐿13\n(3) |\n‖𝒄𝐿13\n(3) ‖ ‖𝒄̃𝐿13\n(3) ‖\n (16) 346 \nwhere 𝒂𝐿13, 𝒃𝐿13, and 𝒄𝐿13 contain the true subject, feature and time loadings of component 𝐿13, 347 \nand 𝒂̃𝐿13, 𝒃̃𝐿13, and 𝒄̃𝐿13 are the best matching subject, feature and time loadings in the fitted model. 348 \nThe FMS of the other components are defined equivalently. In this paper, the best matching 349 \ncomponent is found using the Hungarian algorithm [46, 47]  for all pairwise combinations of 350 \ncomponents such that the FMS is maximised  and the expected CLD structure is recovered 351 \n(Supplementary Methods). The FMS is a value in the range [0,1], where 𝐹𝑀𝑆 = 1 corresponds to 352 \nthe model correctly recovering the input loadings. 353 \nThe recovery of all input loadings can be quantified through 𝐹𝑀𝑆𝑜𝑣𝑒𝑟𝑎𝑙𝑙, defined as 354 \n𝐹𝑀𝑆𝑜𝑣𝑒𝑟𝑎𝑙𝑙 = 𝐹𝑀𝑆𝐶123 + 𝐹𝑀𝑆𝐶123 + 𝐹𝑀𝑆𝐿12 + 𝐹𝑀𝑆𝐿13 + 𝐹𝑀𝑆𝐿23 + 𝐹𝑀𝑆𝐷1 + 𝐹𝑀𝑆𝐷2 + 𝐹𝑀𝑆𝐷3\n8  (17) 355 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n16 \n \n \nwhere 𝐹𝑀𝑆𝑜𝑣𝑒𝑟𝑎𝑙𝑙 = 1 corresponds to all input components being recovered correctly. Recovery of 356 \nthe CLD structure is assessed using the Lambda Similarity Index, which is defined and reported in 357 \nthe Supplementary Materials. 358 \nExample dataset: Jakobsen2025 359 \nWe include an application study of ACMTF -R on longitudinally measured human milk (HM) 360 \nmicrobiome, HM metabolomics and infant faecal microbiome data. Full details on the preprocessing 361 \nof the data is outlined in the original paper [48].  362 \nBriefly, the zOTU data and taxonomic information of the infant faecal microbiome and mother milk 363 \nmicrobiome were processed in R (v4.4.1, [49]). For the infant faecal microbiome data, features were 364 \nkept if they had ≤75% sparsity across the dataset. This step resulted in 93 out of 565 features being 365 \nselected. For the mother milk microbiome data, features were kept if they had ≤85% sparsity across 366 \nthe dataset. This step resulted in 115 out of 707 features being selected. We performed a centred 367 \nlog-ratio transformation with a pseudo -count of 1 to correct for compositionality [50, 51] . 368 \nSubsequently the data was converted to a three -way tensor, keeping missing samples as a row of 369 \nNAs. This resulted in three-way arrays of size 160 subjects × 93 microbial taxa × 3 time points and 370 \n169 subjects × 115 microbial taxa × 4 time points containing microbial abundances for the infant 371 \nfaecal and mother milk microbiomes, respectively. 372 \nThe human milk (HM) metabolomics data were processed using R (v4.4.1, [49]). Values below the 373 \ndetection limit were imputed with a random value between 0 and the detection limit per metabolite 374 \nto preserve their distribution. Next, the dataset was (natural) log transformed to stabilise the 375 \nvariance. The dataset was then converted to a three-way tensor of size 165 subjects × 70 metabolites 376 \n× 4 time points. 377 \nSince ACMTF -R assumes that the subject mode is shared between all data blocks, we took the 378 \nintersection of the subject identifiers across all three datasets and removed the other subjects. This 379 \nhomogenised the subject mode size 𝐼 across all data blocks to 158 subjects.  The data w ere then 380 \ncentred across the subject mode such that variation per time point is only due to between -subject 381 \nvariation, and scaled within the feature mode to make all features equally important for the 382 \nmodelling procedure [52, 53]. As a result, the models will focus on inter-subject differences and the 383 \nfeatures involved, while the time mode loadings indicate when these differences are largest. 384 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n17 \n \n \nAfter centring and scaling, the datasets were block scaled to Frobenius norm 1 to equali se the 385 \nvariance across the blocks.  The response variable 𝒚 of size 158 × 1 contained the maternal pre -386 \npregnancy BMI for all mother-infant dyads, which was centred and subsequently scaled to norm 1 387 \nprior to modelling.  ACMTF-R models were then created using the appropriate number of 388 \ncomponents and regressing on maternal pre-pregnancy BMI. 389 \nSelecting the appropriate number of ACMTF components 390 \nThis paper presents a two-armed approach for selecting the appropriate number of components for 391 \nACMTF: a K-fold cross-validation scheme, and a random initialisation scheme. The full procedure 392 \nis available as the ACMTF_modelSelection() function in the supplied R package. 393 \nIn the K-fold cross-validation scheme, one fold of the subjects is withheld for all blocks as test set 394 \ndata, while randomly initialised ACMTF models are fitted on the training data using 𝑓 components, 395 \nkeeping all other model parameters ( 𝛼, 𝛽, 𝜖) fixed. To avoid data leakage, the means, standard 396 \ndeviations, and norms that are removed from the training data are also removed from the test set 397 \ndata. Model stability is then reported by pairwise comparing the models with the lowest loss per 398 \nfold and calculating the 𝐹𝑀𝑆𝐶𝑉, defined as  399 \n𝐹𝑀𝑆𝐶𝑉 = 1\n𝐹 ∑\n|𝒃𝑓\n𝑇𝒃̃𝑓|\n‖𝒃𝑓‖‖𝒃̃𝑓‖\n|𝒄𝑓\n𝑇𝒄̃𝑓|\n‖𝒄𝑓‖‖𝒄̃𝑓‖ \n𝐹\n𝑓=1\n(18) 400 \nwhere 𝒃𝑓, and 𝒄𝑓 contain the feature and time loadings for component 𝑓 in one model, 𝒃̃𝑓, and 𝒄̃𝑓 401 \nare the best matching subject, feature and time loadings for component 𝑓 in a second model, and 402 \nthe subject mode term is eliminated from the calculated due to not being comparable between CV 403 \nfolds  [15, 54]. The matching of components is achieved using the Hungarian algorithm [46, 47] for 404 \nall pairwise combinations of components such that the FMS is maximised. Hence, the distribution 405 \nof FMS values across all pairwise comparisons of CV (or jack -knife) folds give an indication of the 406 \nstability of the decomposition for a given number of components 𝐹. Previous work has established 407 \nthat 95% of folds should have 𝐹𝑀𝑆𝐶𝑉 ≥ 0.95 for three-way data blocks and 𝐹𝑀𝑆𝐶𝑉 ≥ 0.9 for two-408 \nway data blocks in the ACMTF case for the decomposition to be considered replicable [54]. The 409 \nappropriate number of components is expected to have a stable 𝐹𝑀𝑆𝐶𝑉 near one across all 410 \ncomparisons. 411 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n18 \n \n \nIn the second arm of the approach, randomly initialised ACMTF models are fitted to the full data 412 \nusing 𝑓 components. Model sensitivity to the starting values and to local minima is then assessed by 413 \npairwise comparing the models and calculating the 𝐹𝑀𝑆𝑟𝑎𝑛𝑑𝑜𝑚 through 414 \n𝐹𝑀𝑆𝑟𝑎𝑛𝑑𝑜𝑚 = 1\n𝐹 ∑\n|𝒂𝑓\n𝑇𝒂̃𝑓|\n‖𝒂𝑓‖‖𝒂̃𝑓‖\n𝐹\n𝑓=1\n|𝒃𝑓\n𝑇𝒃̃𝑓|\n‖𝒃𝑓‖‖𝒃̃𝑓‖\n|𝒄𝑓\n𝑇𝒄̃𝑓|\n‖𝒄𝑓‖‖𝒄̃𝑓‖ (19) 415 \nwhere 𝒂𝑓, 𝒃𝑓, and 𝒄𝑓 contain the subject, feature and time loadings for component 𝑓 in one model, 416 \nand 𝒂̃𝑓, 𝒃̃𝑓, and 𝒄̃𝑓 are the best matching subject, feature and time loadings for component 𝑓 in a 417 \nsecond model [15, 54] . The matching of components is again achieved using the Hungarian 418 \nalgorithm [46, 47] for all pairwise combinations of components such that the FMS is maximised. 419 \nThe distribution of FMS values across all pairwise comparisons of randomly initialised models for a 420 \ngiven number of components 𝐹 gives an indication of the stability of the decomposition.  The 421 \nappropriate number of components is expected to have a stable 𝐹𝑀𝑆𝑟𝑎𝑛𝑑𝑜𝑚 near one across all 422 \nrandomly initialised models. 423 \nAdditionally, it may be possible for ACMTF and ACMTF-R to split up common or local components 424 \ninto distinct components with highly similar subject mode loadings when too many components 425 \nare chosen. This phenomenon can be detected by pairwise comparing all  columns in the  subject 426 \nmode per model and calculating the maximum absolute Tucker Congruence Coefficient 𝜑 [43–45]. 427 \nWe define this as the Degeneracy Score, which is computed through 428 \n𝐷𝑒𝑔𝑒𝑛𝑒𝑟𝑎𝑐𝑦 𝑆𝑐𝑜𝑟𝑒= max\n𝑝≠𝑞\n|𝜑(𝒂𝑝, 𝒂𝑞)| (20) 429 \nwhere 𝑝 and 𝑞 indicate component numbers of the subject mode within the same model,  and 𝑝 ≠430 \n𝑞. The Degeneracy Score is in [0, 1], will be close to one when a factor is split into two highly 431 \ncollinear factors, and will have intermediate values otherwise. This procedure is only performed for 432 \nthe randomly initialised ACMTF models since the complete subject mode is needed. The appropriate 433 \nnumber of components is expected to have a low Degeneracy Score across all randomly initialised 434 \nmodels. 435 \nSelecting the appropriate number of ACMTF-R components 436 \nThis paper presents a  similar two-armed approach for selecting the appropriate number of 437 \ncomponents for ACMTF-R, reporting a few additional metrics related to the prediction and fit of 𝒚. 438 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n19 \n \n \nThe full procedure is available as the ACMTFR_modelSelection() function in the supplied R 439 \npackage. 440 \nIn the cross-validation arm, the root mean squared error of cross-validation (RMSECV) is calculated 441 \nby using the models fitted on the training set to predict 𝒚 in the test set. The appropriate number of 442 \ncomponents is expected to minimise the RMSECV while having a stable 𝐹𝑀𝑆𝐶𝑉 near one across all 443 \ncomparisons. 444 \nIn the random initialisation arm, the root mean squared error (RMSE) of 𝒚 and the variance 445 \nexplained in 𝒚 are reported using the ACMTF -R model with the lowest loss. The appropriate 446 \nnumber of components is expected to minimise the RMSE while having a low Degeneracy Score 447 \nand a stable 𝐹𝑀𝑆𝑟𝑎𝑛𝑑𝑜𝑚 near one across all comparisons.  448 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n20 \n \n \nResults 449 \nSimulation 1: ACMTF-R can detect a small, hidden component by leveraging 𝒚 450 \nWe evaluated how the tuning parameter 𝜋 influences the recovery of a local component that is 451 \nweakly present in  some data blocks and is strongly reflected in 𝒚. This was done by applying 452 \nACMTF-R to a simulated three-tensor dataset with a known structure containing 2 common, 3 local, 453 \nand 3 distinct components and 40% noise. All the components were norm 1, except for the second 454 \nlocal component (𝐿13), which  was present across 𝑿(1) and 𝑿(3) with norm 𝛿 = 0.03. Since this 455 \ncomponent was equal to 𝒚 with 5% noise, we expected that lowering 𝜋 would enhance its detection. 456 \nBy systematically varying the tuning parameter 𝜋 from 1 (fully focused on explaining the data 457 \nblocks) to 0 (fully focused on predicting 𝒚), we examined whether outcome information could reveal 458 \na latent structure that would otherwise remain hidden (Figure 3). One hundred randomly initialised 459 \neight-component ACMTF-R and ACMTF models were created for each setting of 𝜋. 460 \nAs assessed by the Factor Match Score, the lowest -loss ACMTF -R models using 0.6 ≤ 𝜋 ≤ 0.9 461 \nshowed an improvement in factor recovery for all components compared to the lowest-loss ACMTF 462 \nmodel (Figure 3A, Supplementary Figure 5). This was also the case for the input factor recovery of 463 \nthe hidden component 𝐿13 (Figure 3B, Supplementary Figure 6). Meanwhile, the fraction of models 464 \nthat correctly identified the input CLD structure was comparable to ACMTF (Supplementary Figure 465 \n7-8). 466 \nNext, the ACMTF-R models with the lowest loss  for each 𝜋 setting were selected for comparison 467 \nwith the lowest -loss ACMTF model ( Figure 3C-H). The absolute Tucker Congruence Coefficient 468 \n(𝜙) was used to compare the best matching recovered factor loadings with the input loadings of 𝐿13 469 \nper mode (Supplementary Methods). This analysis revealed that the recovered shared subject mode 470 \n(Figure 3C) and time mode loadings of 𝑿(1) (Figure 3E) and 𝑿(3) (Figure 3G) were almost identical 471 \nto the input loadings in the ACMTF-R models for 0.6 ≤ 𝜋 ≤ 0.9, while ACMTF failed to identify 472 \nthem. While the recovered feature mode loadings in 𝑿(1) (Figure 3D) and 𝑿(3) (Figure 3F) were not 473 \nexactly equal to the input factor of the hidden component, their orde r largely matched the input 474 \n(Kendall’s 𝜏 = 0.35 − 0.45 for 0.4 ≤ 𝜋 ≤ 0.9; Supplementary Table 1). 475 \nTaken together, ACMTF-R recovered the hidden component 𝐿13 across intermediate values of 𝜋 by 476 \nleveraging 𝒚, while ACMTF did not.477 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n21 \n \n \n 478 \nFigure 3: recovery of a small, hidden component 𝑳𝟏𝟑 using ACMTF and  ACMTF-R. One hundred 479 \nrandomly initialised ACMTF -R models were fitted per setting of tuning parameter 𝜋. These were 480 \ncompared with one hundred randomly initialised ACMTF models. (A) Overall factor recovery, 481 \nmeasured using 𝐹𝑀𝑆𝑜𝑣𝑒𝑟𝑎𝑙𝑙, was higher in ACMTF -R using 0.6 ≤ 𝜋 ≤ 0.9 compared to ACMTF.  (B) 482 \nLikewise, the 𝐹𝑀𝑆𝐿13 maximum was higher for ACMTF-R than ACMTF. Red data points in (A,B) indicate 483 \nthe model with the lowest loss, which were compared in (C -H). The absolute Tucker Congruence 484 \nCoefficient (𝜙) between the recovered loadings and  the input (C) subject mode loadings, (D) feature 485 \nmode loadings of 𝑿(1), (E) time mode loadings of 𝑿(1), (F) feature mode loadings of 𝑿(3), and (G) time 486 \nmode loadings of 𝑿(3) for hidden component 𝐿13. For all modes, factor recovery using ACMTF-R with 487 \nintermediate values of 𝜋 was higher compared to ACMTF. 488 \n 489 \n 490 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n22 \n \n \nSimulation 2: the benefit of supervision using 𝒚 depends on the size of the hidden component 491 \nWe next evaluated how the size of the hidden component 𝐿13 affected the improvement of ACMTF-492 \nR compared to ACMTF. This was done using the same three-tensor setup as in the first simulation, 493 \nexcept that the size of the  second local component (𝐿13) was varied between 𝛿 = 0.01 to 𝛿 = 0.5. 494 \nWe expected the benefit of supervision to decrease as the size of the hidden component 𝛿 increased. 495 \nAs before, the tuning parameter 𝜋 was varied from 0.1 to 0.9 and one hundred randomly initialised 496 \neight-component ACMTF-R and ACMTF models were created for each setting (Figure 4). 497 \nInput factor recovery across all components improved in ACMTF-R compared to ACMTF at 0.4 ≤498 \n𝜋 ≤ 0.9 when the hidden component was small ( 0.01 ≤ 𝛿 ≤ 0.05), while it deteriorated at lower 499 \nvalues of 𝜋 regardless of hidden component size due to overfitting (Figure 4A). Meanwhile, factor 500 \nrecovery of the hidden component 𝐿13 improved in ACMTF -R at 0.4 ≤ 𝜋 ≤ 0.9 only at 501 \nintermediate component sizes (0.03 ≤ 𝛿 ≤ 0.05) (Figure 4B). At 𝛿 > 0.05, ACMTF and ACMTF-R 502 \nwere both able to identify 𝐿13, while at 𝛿 < 0.03 neither method could recover it (Supplementary 503 \nFigure 9). 504 \nThe recovered shared subject mode loadings of the hidden component showed a consistent recovery 505 \nat small component sizes (0.01 -0.05), regardless of 𝜋 (Figure 4C). Meanwhile, the recovery of the 506 \nfeature and time mode loadings improved at similar component sizes only for intermediate values 507 \nof 𝜋 (0.4-0.9) (Figure 4D-G). Taken together, the simulation revealed that a hidden component that 508 \nis roughly 1/20th the size of the other components can be successfully recovered by supervising with 509 \nACMTF-R using intermediate values of 𝜋. 510 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n23 \n \n \n 511 \nFigure 4. Improvement in f actor recovery for varying sizes of the hidden component 𝑳𝟏𝟑 using 512 \nACMTF-R. One hundred randomly initialised eight -component ACMTF -R and ACMTF models were 513 \nfitted for each combination of hidden component size ( 𝛿) and tuning (𝜋), after which the lowest-loss 514 \nmodels were compared. The improvement in factor recovery is computed by subtracting ACMTF 515 \nperformance from ACMTF-R performance in all panels. Improvement per setting for (A) overall factor 516 \nrecovery, (B) factor recovery of the hidden component, (C) shared subject mode recovery of the hidden 517 \ncomponent, (D) feature mode recovery of 𝑿(1) for the hidden component, (E) time mode recovery of 518 \n𝑿(1) for the hidden component, (F) feature mode recovery of 𝑿(3) for the hidden component, (G) time 519 \nmode recovery of 𝑿(3) for the hidden component.  Improvement compared to ACMTF is shown in 520 \ngreen, while deterioration is shown in purple.  A version of this figure without subtraction of ACMTF 521 \nperformance is shown in Supplementary Figure 9. 522 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n24 \n \n \nSimulation 3: the benefit of supervision using 𝒚 is consistent for biologically realistic noise levels 523 \nWe evaluated how the amount of noise on the independent data blocks 𝑿(1), 𝑿(2), and 𝑿(3) affected 524 \nthe improvement in recovering the hidden component 𝐿13 using ACMTF-R compared to ACMTF. 525 \nThis was done using the same three-tensor setup as in the first simulation, except that the amount 526 \nof noise on the independent data blocks was varied from 𝜂𝑋 = 0.0 to 𝜂𝑋 = 100. We expected the 527 \nbenefit of supervision to increase as the amount of noise increased. As before, the tuning parameter 528 \n𝜋 was varied from 0.1 to 0.9 and one hundred randomly initialised eight-component ACMTF-R and 529 \nACMTF models were created for each setting (Figure 5). 530 \nInput factor recovery across all components improved in ACMTF -R compared to ACMTF at most 531 \n𝜋 values (0.3-0.9) at 0-50% noise (Figure 5A). However, at 10-30% noise, the improvement was low 532 \ndue to both ACMTF and ACMTF-R recovering the hidden component (Supplementary Figure 10). 533 \nLikewise, factor recovery of the hidden component 𝐿13 improved only at intermediate noise levels 534 \n(30-50%) and at intermediate 𝜋 values (0.6 -0.9) ( Figure 5B ). When no noise  was present, the 535 \nimprovement in factor recovery of the hidden component was zero because both methods were 536 \nunable to identify it (Supplementary Figure 10). 537 \nThe recovered shared subject mode loadings of the hidden component showed a consistent recovery 538 \nat all noise levels and values of 𝜋 (Figure 5C). The recovery of the feature mode loadings in block 539 \n𝑿(1) improved only at intermediate noise levels (20-50%) and intermediate 𝜋 values (0.6-0.9), while 540 \nthe recovery of the feature mode loadings in block 𝑿(3) were largely independent of the noise level 541 \n(Figure 5D,F). Surprisingly, the recovery of the time mode loadings in block 𝑿(1) deteriorated at 542 \nhigh noise levels (50+%), while the recovery of the time mode loadings in block 𝑿(3) improved at 543 \nthose levels (Figure 5E,G). Taken together, this simulation revealed that factor recovery by ACMTF-544 \nR improves at biologically realistic noise levels (>30%).  545 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n25 \n \n \n 546 \nFigure 5. The improvement in factor recovery by ACMTF -R under different amounts of noise on X. 547 \nOne hundred randomly initialised eight-component ACMTF-R and ACMTF models were fitted for each 548 \ncombination of noise on X ( 𝜂𝑋) and tuning  parameter (𝜋), after which the lowest -loss models were 549 \ncompared. The improvement in factor recovery is computed by subtracting ACMTF performance from 550 \nACMTF-R performance in all panels. Improvement per setting for (A) overall factor recovery, (B) factor 551 \nrecovery of the hidden component, (C) shared subject mode recovery of the hidden component, (D) 552 \nfeature mode recovery of 𝑿(1) for the hidden component, (E) time mode recovery of 𝑿(1) for the 553 \nhidden component, (F) feature mode recovery of 𝑿(3) for the hidden component, (G) time mode 554 \nrecovery of 𝑿(3) for the hidden component. Improvement compared to ACMTF is shown in green, 555 \nwhile deterioration is shown in purple.  A version of this figure without subtraction of ACMTF 556 \nperformance is shown in Supplementary Figure 10. 557 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n26 \n \n \nSimulation 4: the benefit of supervision using 𝒚 is largely independent of the amount of noise on 𝒚 558 \nWe evaluated how the amount of noise on 𝒚 affects the improvement in recovering the hidden 559 \ncomponent 𝐿13 using ACMTF-R compared to ACMTF. This was done using the same three-tensor 560 \nsetup as in the first simulation, except that the amount of noise on 𝒚 was varied from 𝜂𝑦 = 0.0 to 561 \n𝜂𝑦 = 100. We expected the benefit of supervision to decrease as the amount of noise on 𝒚 increased. 562 \nAs before, the tuning parameter 𝜋 was varied from 0.1 to 0.9 and one hundred randomly initialised 563 \neight-component ACMTF-R and ACMTF models were created for each setting (Figure 6). 564 \nInput factor recovery across all components improved in ACMTF -R compared to ACMTF at 565 \nintermediate 𝜋 values (0.4-0.9) for most noise levels (0-100%) (Figure 6A). Likewise, factor recovery 566 \nof the hidden component 𝐿13 improved at intermediate 𝜋 values (0.5-0.9) for most noise levels (0-567 \n100%) (Figure 6B). At lower 𝜋 values, factor recovery in ACMTF-R deteriorated due to the strong 568 \nfocus on predicting 𝒚, due to ACMTF-R not finding the hidden component (Supplementary Figure 569 \n11). At very high noise levels (100+%), neither ACMTF nor ACMTF -R were able to recover the 570 \nhidden component. 571 \nThe recovered shared subject mode loadings of the hidden component showed a consistent recovery 572 \nat most noise levels (0-100%) and intermediate values of 𝜋 (0.4-0.9) (Figure 6C). Similarly, recovery 573 \nof the feature and time mode loadings was consistent at most noise levels (0-100%) and intermediate 574 \nvalues of 𝜋 (0.4-0.9) (Figure 6D-G). Taken together, the simulation revealed that recovery of the 575 \nhidden component by ACMTF-R is insensitive to biologically realistic noise levels on 𝒚. 576 \n  577 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n27 \n \n \n 578 \nFigure 6. The improvement in factor recovery by ACMTF -R under different amounts of noise on 𝒚. 579 \nOne hundred randomly initialised eight-component ACMTF-R and ACMTF models were fitted for each 580 \ncombination noise on 𝒚 (𝜂𝑦) and tuning  parameter (𝜋), after which the lowest -loss models were 581 \ncompared. The improvement in factor recovery is computed by subtracting ACMTF performance from 582 \nACMTF-R performance in all panels. Improvement per setting for (A) overall factor recovery, (B) factor 583 \nrecovery of the hidden component, (C) shared subject mode recovery of the hidden component, (D) 584 \nfeature mode recovery of 𝑿(1) for the hidden component, (E) time mode recovery of 𝑿(1) for the 585 \nhidden component, (F) feature mode recovery of 𝑿(3) for the hidden component, (G) time mode 586 \nrecovery of 𝑿(3) for the hidden component. Improvement compared to ACMTF is shown in green, 587 \nwhile deterioration is shown in purple.  A version of this figure without subtraction of ACMTF 588 \nperformance is shown in Supplementary Figure 11. 589 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n28 \n \n \nApplication: ACMTF-R enhances multi-omics integration compared to ACMTF and NPLS  590 \nTo demonstrate the advantage of ACMTF -R for integrating multi -omics data, we analysed a 591 \nlongitudinal dataset comprising human milk (HM) metabolomics, HM microbiome and infant faecal 592 \nmicrobiome data, using pre-pregnancy BMI (ppBMI)  as the outcome  (𝒚). To contextuali se 593 \nperformance, we evaluated the results against ACMTF , which focuses only on the CLD structure, 594 \nand N-way Partial Least Squares (NPLS) models per block [48], which focus only on prediction. We 595 \nexpected ACMTF-R to identify more strongly maternal ppBMI associated components compared to 596 \nACMTF. 597 \nA combined random initialisation and cross-validation procedure was used to select the number of 598 \ncomponents for ACMTF and ACMTF-R (Methods). The following tuning parameter values were 599 \ninvestigated: 𝜋 = 0.75, 𝜋 = 0.80, 𝜋 = 0.85, and 𝜋 = 0.90. For each 𝜋 value, 10 randomly initialised 600 \nmodels were fitted and a 10 -fold cross-validation was performed across models with one to five 601 \ncomponents (Supplementary Figure 12). This analysis revealed that one component was optimal for 602 \nall ACMTF-R settings. For ACMTF, a similar procedure indicated that two components were 603 \nappropriate, each explaining a different source of variation  (Supplementary Figure 13). For NPLS, 604 \ncross-validation indicated that one component minimised RMSECV for the infant faecal and HM 605 \nmicrobiome data, while two components were required for the HM metabolomics data 606 \n(Supplementary Figure 14). 607 \nAfter selecting the number of components, 100 randomly initialised one -component ACMTF-R 608 \nmodels were fitted for each 𝜋 value. The model with the lowest overall loss in each case was selected, 609 \nand the variance explained was computed ( Table 1). To contextualise performance, the ACMTF 610 \nmodel explained 14.34%, 9.24%, and 7.03% of the variation in the infant faecal microbiome, HM 611 \nmicrobiome, and HM metabolomics data blocks, respectively, while explaining  9.28% of the 612 \nvariation in 𝒚. Additionally, the NPLS models explained 5.78%, 5.41%, and 8.25% of the variation 613 \nin the infant faecal microbiome, HM microbiome, and HM metabolomics datasets, respectively, 614 \nwhile explaining 11.97%, 11.22%, and 35.95% of the variation in 𝒚. 615 \nThe ACMTF-R model with 𝜋 = 0.90 explained 2.08-8.31% of the variation in the independent data 616 \nblocks, and 24.08% of the variation in 𝒚. At 𝜋 = 0.85, this dropped to between 1.78 -5.38% of the 617 \nvariation in the independent data blocks, and 62.06% in 𝒚. At 𝜋 = 0.80 and 𝜋 = 0.75, the models 618 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n29 \n \n \nwere most heavily tuned towards predicting 𝒚, yielding 1 -2% of explained variation in the 619 \nindependent data blocks and capturing 90+% of the variation in 𝒚 – a clear indication of overfitting. 620 \nTable 1. Variance explained per dataset for the fitted NPLS, ACMTF and ACMTF-R models. 𝐹: the 621 \nnumber of components. 622 \nDataset Model 𝐹 \nVariation explained (%) \nInfant faecal \nmicrobiome \nMother milk \nmicrobiome \nMother milk \nmetabolome \nppBMI \n(𝒚) \nInfant faecal \nmicrobiome NPLS 1 5.78 - - 11.97 \nHM \nmicrobiome NPLS 1 - 5.41 - 11.22 \nHM \nmetabolome NPLS 2 - - 8.25 35.95 \nAll \nACMTF-R \n(𝜋 = 0.75) 1 1.52 1.15 1.24 96.55 \nACMTF-R \n(𝜋 = 0.80) 1 2.45 1.91 1.41 89.97 \nACMTF-R \n(𝜋 = 0.85) 1 5.38 4.35 1.78 62.06 \nACMTF-R \n(𝜋 = 0.90) 1 8.31 6.69 2.08 24.08 \nACMTF 2 14.34 9.24 7.03 9.28 \n 623 \nMLR was performed on the  subject mode loadings to test for associations with subject metadata, 624 \nincluding maternal pre-pregnancy BMI, infant growth at six months, birth mode, maternal secretor 625 \n(𝑆𝑒), and maternal Lewis (𝐿𝑒) blood group status (Table 2). Maternal 𝑆𝑒 and 𝐿𝑒 status determine the 626 \npresence of specific fucosylated human milk oligosaccharides, which may influence the composition 627 \nof the infant faecal microbiome.  628 \nThe subject mode loadings of the NPLS models all contained at least one component that was  629 \nassociated with maternal pre-pregnancy BMI. While both components of the ACMTF model were 630 \nassociated with both ppBMI and maternal secretor status , the first component was most strongly 631 \nassociated with maternal secretor status and the second was most strongly associated with ppBMI . 632 \nThis was not the case for the ACMTF-R models with 𝜋 = 0.90 and 𝜋 = 0.85, which both contained 633 \nsubject mode loadings associated exclusively with ppBMI  after correction . The subject mode 634 \nloadings of the ACMTF-R models with lower values of 𝜋 were strongly associated with ppBMI due 635 \nto overfitting. 636 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n30 \n \n \nDue to the extraction of a component exclusively related to maternal ppBMI, as well as the trade -637 \noff between explaining the independent data and predicting 𝒚, the ACMTF-R model using 𝜋 = 0.90 638 \nwas selected for further interpretation. 639 \nTable 2. MLR results for the subject mode loadings of the created multi -way models . ppBMI: 640 \nmaternal pre-pregnancy body mass index. 𝑆𝑒: maternal secretor status. 𝐿𝑒: maternal Lewis status.  641 \nWHZ: infant weight -of-height z-score at 6 months.  Subject mode loadings were normalised to 642 \nnorm 1 prior to MLR fitting for comparability of the results between the models. Values indicate 643 \nMLR coefficients times one thousand. Stars indicate the original p -values from the MLR model (*: 644 \n𝑝 < 0.05, **: 𝑝 < 0.01, ***: 𝑝 < 0.001), while hats indicate Benjamini-Hochberg corrected p-values 645 \n(^: 𝑞 < 0.05, ^^: 𝑞 < 0.01, ^^^: 𝑞 < 0.001). Coefficients with a significant p-value after correction 646 \nare made bold. The correlation matrix of the metadata is reported in Supplementary Table 2. 647 \nDataset Model F ppBMI Birth \nmode Se Le WHZ \nInfant faecal \nmicrobiome NPLS 1 5.9***, ^^^ 26.2 18.7 -67.4 1.5 \nHM \nmicrobiome NPLS 1 5.3***, ^^^ 22.5 10.1 -25.4 -0.6 \nHM \nmetabolome NPLS \n1 8.7***, ^^^ 11.1 19.7 -78.3 -12.6*, ^ \n2 -0.2 -12.8 -111***, ^^^ 43.0 7.0 \nAll \nACMTF-R \n(𝜋 = 0.75) 1 15.0***, ^^^ 6.7 3.0* -18.5 -0.9 \nACMTF-R \n(𝜋 = 0.80) 1 -14.7***, ^^^ -10.8 -5.3* 31.7 1.4 \nACMTF-R \n(𝜋 = 0.85) 1 -12.5***, ^^^ -18.1 -12.3* 59.6 2.4 \nACMTF-R \n(𝜋 = 0.90) 1 -8.3***, ^^^ -21.5 -21.8 79.4 2.8 \nACMTF \n1 -2.9*, ^ -2.0 -74.8***, ^^^ 65.5 3.2 \n2 -4.4***, ^^ -40.1 -41.3*, ^ 72.3 1.0 \n 648 \nThe 𝚲-matrix of the ACMTF-R model using 𝜋 = 0.90 was examined to interpret the common, local, 649 \nand distinct (CLD) structure of the component. In this model, 𝜆 values of the only available 650 \ncomponent were equal to 0. 34, 0.29, and 0.15, for the infant faecal microbiome, HM microbiome, 651 \nand HM metabolomics, respectively. Based on the relative magnitudes of the 𝜆 values, the 652 \ncomponent appears to be primarily local to the infant faecal and HM microbiome blocks, with a 653 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n31 \n \n \nminor contribution of the HM metabolomics block. This suggests that ppBMI -associated variation 654 \nwas mainly present in the microbiomes of mother and infant. 655 \nIn the infant faecal microbiome, genera associated with higher ppBMI included  Bifidobacterium 656 \nbifidum, Escherichia coli, Staphylococcus epidermidis, and B. breve (Figure 7A). Genera enriched 657 \nin infants of low-ppBMI mothers included Parabacteroides, Enterococcus and Streptococcus. In the 658 \nHM microbiome, higher ppBMI was associated with increased abundance of Staphylococcus and 659 \nStreptococcus spp., while genera associated with lower ppBMI included various taxa such as 660 \nClostridium sensu stricto, Collinsella, and Lactobacillus (Figure 7C). In the HM metabolome, higher 661 \nppBMI was associated with increased levels of many unknown HMOs, caffeine, and LDFT ( Figure 662 \n7E). Metabolites with increased levels in low -ppBMI mothers included several amino acids and 663 \nrelated compounds such as 2-aminobutyrates, aspartate, isoleucine, and methionine. 664 \nDue to the pre -processing of the data, the  largest absolute time loading in the ACMTF -R model 665 \nindicates the time points at which inter -subject differences related to ppBMI are greatest ( Figure 666 \n7B,D,F). Overall, the time loadings were relatively flat, suggesting that ppBMI-associated variation 667 \nremained stable throughout the sampling period. In the HM microbiome, however, these 668 \ndifferences peaked at day 30 and declined thereafter. 669 \nTaken together, ACMTF -R was able to identify variation that was strongly and exclusively 670 \nassociated with maternal ppBMI. Additionally, it revealed that this variation was primarily local to 671 \nthe infant faecal and HM microbiome, with a limited contribution from HM metabolome. This 672 \nhighlights the value of incorporating supervision to isolate outcome -relevant signals in complex 673 \nmulti-omics data.  674 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n32 \n \n \n 675 \nFigure 7. Feature and time mode loadings of the selected ACMTF-R model. Feature mode loadings 676 \nof the selected ACMTF -R model using 𝜋 = 0.90 describing (A) the infant faecal microbiome, (C) 677 \nthe HM microbiome, and (E) the HM metabolomics datasets. Time loadings of the same model are 678 \nshown in the same order in (B,D,F). Signs of the feature mode loadings were adjusted post-hoc such 679 \nthat positive values consistently corresponded to higher ppBMI and negative values to lower ppBMI. 680 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n33 \n \n \nDiscussion 681 \nIn this study, we presented ACMTF-Regression (ACMTF-R) as an extension of Advanced Coupled 682 \nMatrix and Tensor Factorization (ACMTF) that incorporates a regression term to model common, 683 \nlocal and distinct (CLD) variation related to a dependent variable. Through simulations, we 684 \ndemonstrated its ability to capture a hidden component, its robustness to noise, and its flexibility in 685 \ntuning between data exploration and predicting. Our application demonstrated how ACMTF -R 686 \nfinds novel relationships for the effect of maternal ppBMI on the HM metabolome, HM microbiome 687 \nand infant faecal microbiome. 688 \nACMTF-R extends the capabilities of ACMTF by allowing explicit modelling of outcome-associated 689 \nvariation, bridging the gap between exploratory data integration and supervised learning [10]. 690 \nCompared to NPLS, which is designed for predictive modelling but lacks the ability to separate 691 \ncommon, local, and distinct variation across multiple blocks, ACMTF -R provides a more nuanced 692 \nview of multi-way data structures [19]. Unlike classical ACMTF, which only uncovers underlying 693 \ndata patterns, ACMTF -R allows for targeted analysis of biological variation of interest [33]. The 694 \nability to tune between the two makes it particularly useful for cases where a balance between 695 \nexploration and prediction is necessary. 696 \nOur simulations revealed that ACMTF-R reliably recovers the underlying CLD-structure of multi-697 \nway datasets across a range of settings. We found that the tuning parameter (𝜋) governs the trade-698 \noff between reconstructing  the input data and predicting the outcome variable. While 𝜋 = 1 699 \ncorresponds to regular ACMTF, decreasing 𝜋 enhances the identification of factors associated with 700 \nthe dependent variable, even revealing weak structures that would otherwise remain undetected. 701 \nHowever, tuning too far towards 𝒚 (𝜋 < 0.3) leads to overfitting and an inaccurate reconstruction. 702 \nThe observation that intermediate values of 𝜋 (0.6-0.9) lead to recovery of the hidden local 703 \ncomponent is largely in agreement with PCovR and MCovR [20, 22].  704 \nAdditionally, the simulations revealed that the hidden component needed to be sufficiently small 705 \n(𝛿 < 0.05) for ACMTF to fail in its recovery. However, this finding is largely dependent on the 706 \nrelative sizes of the other components, which we have kept at Frobenius norm 1 in our simulations. 707 \nFurther research is needed to investigate the recovery of a hidden component using ACMTF and 708 \nACMTF-R in other scenarios. 709 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n34 \n \n \nFinally, the simulations showed that recovery of the hidden component using ACMTF-R was most 710 \naccurate when a biologically realistic amount of noise was present in the data ( 𝜂𝑋 ≤ 0.5, 𝜂𝑦 < 1). 711 \nThe observation that neither ACMTF nor ACMTF-R can identify the hidden component in the no-712 \nnoise case is surprising. This may occur when the noise has a smoothing effect on the loss space, 713 \nresulting in the global minimum being easier to identify, or reduces the collinearity of the loadings 714 \nin a statistical phenomenon comparable to ridge regression [55]. While the noise level for a real 715 \ndataset is often not known, it is reassuring that ACMTF -R can reliably capture the 𝒚-associated 716 \nsubject mode loadings for most tested noise levels (Supplementary Figure 10-11). 717 \nApplication to a longitudinal multi-omics dataset integrating human milk microbiome, human milk 718 \nmetabolome, and infant gut microbiome data further highlighted the utility of ACMTF -R by 719 \nrecovering a n exclusively ppBMI -associated component [48]. In the infant faecal microbiome, 720 \nelevated abundances of Enterococcus sp., E. coli, Clostridium sp., and Staphylococcus epidermidis 721 \nin infants born to high-ppBMI mothers are consistent with prior findings linking maternal obesity 722 \nto pro -inflammatory microbial profiles [56, 57] . In contrast, the increased abundance of 723 \nBifidobacterium breve and Bifidobacterium bifidum in the same group contradicts earlier studies, 724 \nwhich generally associate these taxa with leaner maternal phenotypes [56, 57].  725 \nInterestingly, the genera enriched in infants from low-ppBMI mothers represent a potential novel 726 \nfinding, as comparable reports are currently lacking and warrant further investigation.  In the HM 727 \nmicrobiome, enrichment of Streptococcus and Staphylococcus spp. in high -ppBMI mothers is in 728 \nagreement with previous research  [58–60]. However, comparable data for taxa enriched in low -729 \nppBMI mothers were not available for direct comparison.  730 \nIn the HM metabolomics data, elevated levels of lactate and multiple unidentified HMOs were 731 \nobserved in high-ppBMI mothers. While the relationship between ppBMI and HMO concentrations 732 \nhas been investigated in several studies, the results remain inconsistent [61, 62]. One study did 733 \nreport higher levels of a lactate derivative in the milk of obese mothers [61], suggesting a possible 734 \nlink to altered maternal energy metabolism. Conversely, the observation of increased 735 \nconcentrations of several amino acids in low-ppBMI mothers is supported by previous findings [61], 736 \nas well as by recent evidence showing enrichment of amino acid biosynthesis pathways in this group 737 \n[63]. These results show that ACMTF-R provides additional biological insights compared to ACMTF 738 \nand N -PLS, effectively disentangling common, local and distinct variation related to maternal 739 \nppBMI. 740 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n35 \n \n \nWhile ACMTF-R demonstrates strong performance in identifying biologically meaningful patterns, 741 \nseveral challenges remain.  The choice of the tuning parameter 𝜋 is dataset-dependent, requiring 742 \ncareful selection through cross -validation or heuristic approaches.  Additionally, in multi -omics 743 \napplications, different data blocks may contain varying amounts of biologically relevant 744 \ninformation. Allowing block -specific 𝜋 values could improve model flexibility, enabling certain 745 \ndatasets to contribute more strongly to the dependent variable while ensuring others remain focused 746 \non data structure. Future research should explore whether adaptive weighting strategies or block -747 \nspecific hyperparameter tuning can further improve performance. 748 \nFor ACMTF and ACMTF-R, the single guiding metric for selecting the number of components was 749 \noften the stability of the decomposition through the factor match score (FMS) [15, 54]. This paper 750 \nhas expanded upon this approach by formalising FMS for a cross -validation and random 751 \ninitialisation scheme, as well as developing the degeneracy score diagnostic to identify 752 \noverfactorisation. However, this approach is computationally intensive for the simulations. Future 753 \nwork should investigate if this approach could be accelerated through the Limited -memory 754 \nBroyden-Fletcher-Goldfarb-Shanno [40–42] and AO-ADMM algorithms [64–66]. 755 \nWhile Factor Match Score (FMS) has been a useful diagnostic for model selection, it suffers from 756 \nthe phenomenon of all blocks contributing to all components in ACMTF outlined above. Due to the 757 \nway the structural model of ACMTF is implemented, 𝐹 trilinear components are fitted per data 758 \nblock, even when the rank of that block is lower than 𝐹 [33] This causes the creation of superfluous 759 \nfeature and time mode loading vectors that contribute to the model due to having non-zero 𝜆 values. 760 \nAs a result, the FMS is negatively impacted by considering these components in the calculation. 761 \nThis problem may be avoided by more strongly penalising the 𝚲-matrix through replacement of the 762 \nL1 penalty term by a smooth L0 penalty, as presented in Relaxed ACMTF (RACMTF) [67, 68]. 763 \nFuture research should investigate if this approach yields more stable ACMTF and ACMTF -R 764 \nmodels for multi-omics data. 765 \nThe multi-omics application presented in this study highlights a critical shortcoming of ACMTF in 766 \nthe identification of common, local, and distinct (CLD) sources of variation by revealing that every 767 \ndata block contributes to all identified components. However, the relative magnitude of the 768 \ncontributing blocks captured by 𝜆 still offer a useful interpretation. For example, a 𝚲-matrix 769 \nthreshold could be applied post-hoc to identify CLD-structures [33]. Another approach could be to 770 \npre-define the expected number of shared and distinct components, as suggested by Common and 771 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n36 \n \n \nDiscriminative subspace Non-negative Tensor Factorization (CDNTF) [69]. Future research should 772 \ninvestigate if either approach or a combination of both might yield more interpretable ACMTF 773 \nmodels for multi-omics data.  774 \nFinally, treating 𝒚 as an additional data block for ACMTF could offer a different approach to 775 \nincorporating outcome variation into ACMTF -R. By embedding 𝒚 within the multi -way 776 \ndecomposition framework, rather than optimizing it separately, the model might be able to better 777 \ncapture latent structures linking 𝒚 to the input datasets. However, this approach would lose the 778 \npredictive aspect of ACMTF-R. Future work should focus on the comparability between ACMTF-R 779 \nand adding 𝒚 as a block to ACMTF. 780 \nConclusions 781 \nOur results demonstrate that ACMTF -R is a powerful framework for multi -way data integration 782 \nthat extends ACMTF by incorporating outcome-associated variation. Through simulations, we 783 \nshowed that ACMTF -R can recover a small hidden component, is  robust to noise, and provides 784 \nflexible tuning between data exploration and outcome prediction. Its application to multi -omics 785 \ndata highlights how ACMTF-R finds novel relationships for the effect of maternal ppBMI on the 786 \nHM metabolome, HM microbiome and infant faecal microbiome. 787 \nEthics approval and consent to participate 788 \nNot applicable. 789 \nConsent for publication 790 \nNot applicable. 791 \nAvailability of the data and materials 792 \nThe package used during the current study is available in the GitHub repository 793 \nhttps://github.com/GRvanderPloeg/CMTFtoolbox. The datasets generated and/or analysed during the 794 \ncurrent study are available in the GitHub repository 795 \nhttps://github.com/GRvanderPloeg/CMTF_toolbox_project. 796 \nCompeting interests 797 \nThe authors declare that they have no competing interests. 798 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n37 \n \n \nFunding 799 \nGRvdP was funded by a grant from the University of Amsterdam, Research Priority Area on 800 \nPersonal Microbiome Health. 801 \nFred White was funded by a grant from the University of Amsterdam, Data Science Centre. 802 \nAuthors’ contributions 803 \nG.R. van der Ploeg: conceptualization, visualization, formal analysis, methodology, programming, 804 \nwriting – original draft; F.T.G. White: conceptualization, methodology, programming, writing – 805 \noriginal draft ; R.R. Jakobsen:  programming, writing – review & editing ; J.A. Westerhuis: 806 \nconceptualization, methodology, supervision, formal analysis, writing – original draft; A. Heintz -807 \nBuschart: conceptualization, formal analysis, supervision, project administration, writing – original 808 \ndraft; A.K. Smilde: conceptualization, methodology, supervision, funding acquisition, writing – 809 \nreview & editing. 810 \nAcknowledgements 811 \nWe want to thank Daniël Rademaker (RU/UvA) for their useful discussions and suggestions. 812 \nReferences 813 \n1. Capobianco E. Ten challenges for systems medicine. Front Genet. 2012;3:193. 814 \n2. Davis-Turak J, Courtney SM, Hazard ES, Glen WB, Da Silveira WA, Wesselman T, et al. Genomics 815 \npipelines and data integration: challenges and opportunities in the research setting. Expert Rev Mol 816 \nDiagn. 2017;17:225–37. https://doi.org/10.1080/14737159.2017.1282822. 817 \n3. Zhao Z, Jin VX, Huang Y, Guda C, Ruan J. Frontiers in integrative genomics and translational 818 \nbioinformatics. BioMed Res Int. 2015;2015. 819 \n4. Bordbar A. Interpreting the deluge of omics data: new approaches offer new possibilities. Blood 820 \nTransfus. 2017;15:189. 821 \n5. Chen C, McGarvey PB, Huang H, Wu CH. Protein Bioinformatics Infrastructure for the Integration 822 \nand Analysis of Multiple High -Throughput “omics” Data. Adv Bioinforma. 2010;2010:1 –19. 823 \nhttps://doi.org/10.1155/2010/423589. 824 \n6. Misra BB, Langefeld C, Olivier M, Cox LA. Integrated omics: tools, advances and future approaches. 825 \nJ Mol Endocrinol. 2019;62:R21–45. 826 \n7. Vitorino R. Transforming Clinical Research: The Power of High -Throughput Omics Integration. 827 \nProteomes. 2024;12:25. https://doi.org/10.3390/proteomes12030025. 828 \n8. van der Kloet FM, Sebastián-León P, Conesa A, Smilde AK, Westerhuis JA. Separating common from 829 \ndistinctive variation. BMC Bioinformatics. 2016;17:S195. https://doi.org/10.1186/s12859-016-1037-2. 830 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n38 \n \n \n9. Måge I, Smilde AK, van der Kloet FM. Performance of methods that separate common and distinct 831 \nvariation in multiple data blocks. J Chemom. 2019;33:e3085. https://doi.org/10.1002/cem.3085. 832 \n10. Acar E, Papalexakis EE, Gürdeniz G, Rasmussen MA, Lawaetz AJ, Nilsson M, et al. Structure -833 \nrevealing data fusion. BMC Bioinformatics. 2014;15:239. https://doi.org/10.1186/1471-2105-15-239. 834 \n11. Tan CS, Salim A, Ploner A, Lehtiö J, Chia KS, Pawitan Y. Correlating gene and protein expression data 835 \nusing Correlated Factor Analysis. BMC Bioinformatics. 2009;10:272. https://doi.org/10.1186/1471 -836 \n2105-10-272. 837 \n12. Xiao X, Moreno -Moral A, Rotival M, Bottolo L, Petretto E. Multi -tissue analysis of co -expression 838 \nnetworks by higher -order generalized singular value decomposition identifies functionally coherent 839 \ntranscriptional modules. PLoS Genet. 2014;10:e1004006. 840 \n13. Berger JA, Hautaniemi S, Mitra SK, Astola J. Jointly analyzing gene expression and copy number 841 \ndata in breast cancer using data reduction models. IEEE/ACM Trans Comput Biol Bioinform. 2006;3:2–842 \n16. 843 \n14. Song Y, Westerhuis JA, Smilde AK. Separating common (global and local) and distinct variation in 844 \nmultiple mixed types data sets. J Chemom. 2020;34:e3197. https://doi.org/10.1002/cem.3197. 845 \n15. Acar E, Kolda TG, Dunlavy DM. All -at-once Optimization for Coupled Matrix and Tensor 846 \nFactorizations. 2011. 847 \n16. Acar E, Bro R, Smilde A. Data Fusion in Metabolomics Using Coupled Matrix and Tensor 848 \nFactorizations. Proc IEEE. 2015;103:1602. https://doi.org/10.1109/JPROC.2015.2438719. 849 \n17. Acar E, Bro R, Smilde AK. Data Fusion in Metabolomics Using Coupled Matrix and Tensor 850 \nFactorizations. Proc IEEE. 2015;103:1602–20. https://doi.org/10.1109/JPROC.2015.2438719. 851 \n18. Acar E, Dunlavy DM, Kolda TG, Mørup M. Scalable tensor factorizations for incomplete data. 852 \nChemom Intell Lab Syst. 2011;106:41–56. https://doi.org/10.1016/j.chemolab.2010.08.004. 853 \n19. Bro R. Multiway calibration. Multilinear PLS. J Chemom. 1996;10:47 –61. 854 \nhttps://doi.org/10.1002/(SICI)1099-128X(199601)10:1<47::AID-CEM400>3.0.CO;2-C. 855 \n20. Vervloet M, Kiers HAL, Noortgate WV den, Ceulemans E. PCovR: An R Package for Principal 856 \nCovariates Regression. J Stat Softw. 2015;65:1–14. https://doi.org/10.18637/jss.v065.i08. 857 \n21. De Jong S, Kiers HA. Principal covariates regression: part I. Theory. Chemom Intell Lab Syst. 858 \n1992;14:155–64. 859 \n22. Smilde AK, Kiers HAL. Multiway covariates regression models. J Chemom. 1999;13:31 –48. 860 \nhttps://doi.org/10.1002/(SICI)1099-128X(199901/02)13:1<31::AID-CEM528>3.0.CO;2-P. 861 \n23. Kiers HAL. Towards a standardized notation and terminology in multiway analysis. J Chemom. 862 \n2000;14:105–22. https://doi.org/10.1002/1099-128X(200005/06)14:3<105::AID-CEM582>3.0.CO;2-I. 863 \n24. Kolda TG. Multilinear operators for higher -order decompositions. Sandia National Laboratories 864 \n(SNL), Albuquerque, NM, and Livermore, CA …; 2006. 865 \n25. De Lathauwer L, De Moor B, Vandewalle J. A Multilinear Singular Value Decomposition. SIAM J 866 \nMatrix Anal Appl. 2000;21:1253–78. https://doi.org/10.1137/S0895479896305696. 867 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n39 \n \n \n26. Styan GP. Hadamard products and multivariate statistical analysis. Linear Algebra Its Appl. 868 \n1973;6:217–40. 869 \n27. McDonald RP. A simple comprehensive model for the analysis of covariance structures: Some 870 \nremarks on applications. Br J Math Stat Psychol. 1980;33:161 –83. https://doi.org/10.1111/j.2044 -871 \n8317.1980.tb00606.x. 872 \n28. Rao CR. GENERALIZED INVERSE OF A MATRIX AND ITS APPLICATIONS. In: Le Cam LM, Neyman J, 873 \nScott EL, editors. Theory of Statistics. University of California Press; 1972. p. 601 –20. 874 \nhttps://doi.org/10.1525/9780520325883-032. 875 \n29. Ben -Israel A, Greville TN. Generalized inverses: theory and applications. Springer Science & 876 \nBusiness Media; 2006. 877 \n30. Carroll JD, Chang J -J. Analysis of individual differences in multidimensional scaling via an N -way 878 \ngeneralization of “Eckart-Young” decomposition. Psychometrika. 1970;35:283–319. 879 \n31. Harshman RA. FOUNDATIONS OF THE PARAFAC PROCEDURE: MODELS AND CONDITIONS FOR AN 880 \n“EXPLANATORY” MULTIMODAL FACTOR ANALYSIS. :84. 881 \n32. Bro R. MULTI-WAY ANALYSIS IN THE FOOD INDUSTRY. 882 \n33. Acar E, Papalexakis EE, Gürdeniz G, Rasmussen MA, Lawaetz AJ, Nilsson M, et al. Structure -883 \nrevealing data fusion. BMC Bioinformatics. 2014;15:239. https://doi.org/10.1186/1471-2105-15-239. 884 \n34. Lee S-I, Lee H, Abbeel P, Ng AY. Efficient l1 regularized logistic regression. In: Aaai. 2006. p. 401–8. 885 \n35. Bertrand D, Qannari EM, Vigneau E. Latent root regression analysis: an alternative method to PLS. 886 \nChemom Intell Lab Syst. 2001;58:227–34. https://doi.org/10.1016/S0169-7439(01)00161-7. 887 \n36. Webster JT, Gunst RF, Mason RL. Latent Root Regression Analysis. Technometrics. 1974;16:513 –888 \n22. https://doi.org/10.1080/00401706.1974.10489232. 889 \n37. Harshman RA. Foundations of the PARAFAC procedure: Models and conditions for an\" explanatory\" 890 \nmultimodal factor analysis. 1970. 891 \n38. Bro R. PARAFAC. Tutorial and applications. Chemom Intell Lab Syst. 1997;:23. 892 \n39. Melville J. mize: Unconstrained Numerical Optimization Algorithms. 2019. 893 \n40. Liu DC, Nocedal J. On the limited memory BFGS method for large scale optimization. Math Program. 894 \n1989;45:503–28. https://doi.org/10.1007/BF01589116. 895 \n41. Malouf R. A comparison of algorithms for maximum entropy parameter estimation. In: proceedings 896 \nof the 6th conference on Natural language learning - Volume 20. USA: Association for Computational 897 \nLinguistics; 2002. p. 1–7. https://doi.org/10.3115/1118853.1118871. 898 \n42. Andrew G, Gao J. Scalable training of L1-regularized log-linear models. In: Proceedings of the 24th 899 \ninternational conference on Machine learning. New York, NY, USA: Association for Computing 900 \nMachinery; 2007. p. 33–40. https://doi.org/10.1145/1273496.1273501. 901 \n43. Burt C. The factorial study of temperamental traits. Br J Psychol. 1948. 902 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n40 \n \n \n44. Tucker LR. A method for synthesis of factor analysis studies. Educational Testing Service Princeton, 903 \nNJ; 1951. 904 \n45. Lorenzo-Seva U, ten Berge JMF. Tucker’s congruence coefficient as a meaningful index of factor 905 \nsimilarity. Methodol Eur J Res Methods Behav Soc Sci. 2006;2:57 –64. https://doi.org/10.1027/1614-906 \n2241.2.2.57. 907 \n46. Kuhn HW. The Hungarian method for the assignment problem. Nav Res Logist Q. 1955;2:83 –97. 908 \nhttps://doi.org/10.1002/nav.3800020109. 909 \n47. Kuhn HW. Variants of the hungarian method for assignment problems. Nav Res Logist Q. 910 \n1956;3:253–8. https://doi.org/10.1002/nav.3800030404. 911 \n48. Jakobsen R, Ploeg GR, Sundekilde U, Astono J, Poulsen K, Fuglsang J, et al. Supervised Modelling of 912 \nLongitudinal Human Milk and Infant Gut Microbiome Reveal Maternal Pre -Pregnancy BMI and Early 913 \nLife Growth Interactions. 2025. https://doi.org/10.21203/rs.3.rs-6244750/v1. 914 \n49. R Core Team R. R: A language and environment for statistical computing. 2013. 915 \n50. Gloor GB, Macklaim JM, Pawlowsky-Glahn V, Egozcue JJ. Microbiome Datasets Are Compositional: 916 \nAnd This Is Not Optional. Front Microbiol. 2017;8:2224. https://doi.org/10.3389/fmicb.2017.02224. 917 \n51. Aitchison J. The Statistical Analysis of Compositional Data. J R Stat Soc Ser B Methodol. 918 \n1982;44:139–60. https://doi.org/10.1111/j.2517-6161.1982.tb01195.x. 919 \n52. Bro R, Smilde AK. Centering and scaling in component analysis. J Chemom. 2003;17:16 –33. 920 \nhttps://doi.org/10.1002/cem.773. 921 \n53. van der Ploeg GR, Westerhuis JA, Heintz-Buschart A, Smilde AK. parafac4microbiome: Exploratory 922 \nanalysis of longitudinal microbiome data using Parallel Factor Analysis. 2024. 923 \nhttps://doi.org/10.1101/2024.05.02.592191. 924 \n54. Li L, Yan S, Horner D, Rasmussen MA, Smilde AK, Acar E. Revealing static and dynamic biomarkers 925 \nfrom postprandial metabolomics data through coupled matrix and tensor factorizations. 2024. 926 \n55. Schott JR. Matrix analysis for statistics. John Wiley & Sons; 2016. 927 \n56. Singh SB, Madan J, Coker M, Hoen A, Baker ER, Karagas MR, et al. Associations of maternal pre -928 \npregnancy BMI and gestational weight gain with the infant gut microbiome differ according to delivery 929 \nmode. Int J Obes 2005. 2019;44:23. 930 \n57. Collado MC, Isolauri E, Laitinen K, Salminen S. Effect of mother’s weight on infant’s microbiota 931 \nacquisition, composition, and activity during early infancy: a prospective follow -up study initiated in 932 \nearly pregnancy. Am J Clin Nutr. 2010;92:1023–30. 933 \n58. Lundgren SN, Madan JC, Karagas MR, Morrison HG, Hoen AG, Christensen BC. Microbial 934 \ncommunities in human milk relate to measures of maternal weight. Front Microbiol. 2019;10:2886. 935 \n59. Cabrera -Rubio R, Collado MC, Laitinen K, Salminen S, Isolauri E, Mira A. The human milk 936 \nmicrobiome changes over lactation and is shaped by maternal weight and mode of delivery1234. Am 937 \nJ Clin Nutr. 2012;96:544–51. https://doi.org/10.3945/ajcn.112.037382. 938 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n41 \n \n \n60. Cortés-Macías E, Selma-Royo M, Rio-Aige K, Bäuerl C, Rodríguez-Lagunas MJ, Martínez-Costa C, et 939 \nal. Distinct breast milk microbiota, cytokine, and adipokine profiles are associated with infant growth 940 \nat 12 months: an in vitro host –microbe interaction mechanistic approach. Food Funct. 2023;14:148 –941 \n59. 942 \n61. Isganaitis E, Venditti S, Matthews TJ, Lerin C, Demerath EW, Fields DA. Maternal obesity and the 943 \nhuman milk metabolome: associations with infant body composition and postnatal weight gain. Am J 944 \nClin Nutr. 2019;110:111–20. 945 \n62. Bardanzellu F, Puddu M, Peroni DG, Fanos V. The Human Breast Milk Metabolome in Overweight 946 \nand Obese Mothers. Front Immunol. 2020;11:1533. https://doi.org/10.3389/fimmu.2020.01533. 947 \n63. Sivalogan K, Liang D, Accardi C, Diaz-Artiga A, Hu X, Mollinedo E, et al. Human Milk Composition Is 948 \nAssociated with Maternal Body Mass Index in a Cross-Sectional, Untargeted Metabolomics Analysis of 949 \nHuman Milk from Guatemalan Mothers. Curr Dev Nutr. 2024;8:102144. 950 \nhttps://doi.org/10.1016/j.cdnut.2024.102144. 951 \n64. Schenker C, Cohen JE, Acar E. A Flexible Optimization Framework for Regularized Matrix -Tensor 952 \nFactorizations With Linear Couplings. IEEE J Sel Top Signal Process. 2021;15:506 –21. 953 \nhttps://doi.org/10.1109/jstsp.2020.3045848. 954 \n65. Schenker C, Cohen JE, Acar E. An Optimization Framework for Regularized Linearly Coupled Matrix-955 \nTensor Factorization. In: 2020 28th European Signal Processing Conference (EUSIPCO). Amsterdam, 956 \nNetherlands: IEEE; 2021. p. 985–9. https://doi.org/10.23919/eusipco47968.2020.9287459. 957 \n66. Roald M. MatCoupLy: Learning coupled matrix factorizations with Python. SoftwareX. 958 \n2023;21:101292. https://doi.org/10.1016/j.softx.2022.101292. 959 \n67. Mohimani H, Babaie-Zadeh M, Jutten C. A fast approach for overcomplete sparse decomposition 960 \nbased on smoothed $\\backslashell^{0} $ norm. IEEE Trans Signal Process. 2008;57:289–301. 961 \n68. Rivet B, Duda M, Guérin-Dugué A, Jutten C, Comon P. Multimodal approach to estimate the ocular 962 \nmovements during EEG recordings: A coupled tensor factorization method. In: 2015 37th Annual 963 \nInternational Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). 2015. p. 964 \n6983–6. https://doi.org/10.1109/EMBC.2015.7319999. 965 \n69. Liu W, Chan J, Bailey J, Leckie C, Ramamohanarao K. Mining Labelled Tensors by Discovering both 966 \ntheir Common and Discriminative Subspaces. In: Proceedings of the 2013 SIAM International 967 \nConference on Data Mining (SDM). Society for Industrial and Applied Mathematics; 2013. p. 614 –22. 968 \nhttps://doi.org/10.1137/1.9781611972832.68. 969 \n 970 \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint \n\n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint","source_license":"CC-BY-4.0","license_restricted":false}