{'fraction_non_alphanumeric': 0.029743117152401796, 'fraction_numerical': 0.013362515839189033, 'mean_word_length': 4.673986928104575, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Merging versus Ensembling in Multi-Study Prediction: Theoretical Insight from Random Effects

A critical decision point when training predictors using multiple studies is whether these studies should be combined or treated separately. We compare two multi-study learning approaches in the presence of potential heterogeneity in predictor-outcome relationships across datasets. We consider 1) merging all of the datasets and training a single learner, and 2) multi-study ensembling, which involves training a separate learner on each dataset and combining the predictions resulting from each learner. In a linear regression setting, we show analytically and confirm via simulation that merging yields lower prediction error than ensembling when the predictor-outcome relationships are relatively homogeneous across studies. However, as cross-study heterogeneity increases, there exists a transition point beyond which ensembling outperforms merging. We provide analytic expressions for the transition point in various scenarios, study asymptotic properties, and illustrate how transition point theory can be used for deciding when studies should be combined with an application from metabolomics.


INTRODUCTION
Prediction models trained on a single study often perform considerably worse in external validation than in cross-validation (Castaldi et al., 2011;Bernau et al., 2014). Their generalizability is compromised by overfitting, but also by various sources of study heterogeneity, including differences in study design, data collection and measurement methods, unmeasured confounders, and study-specific sample characteristics (Zhang et al., 2018). Using multiple training studies can potentially address these challenges and lead to more replicable prediction models. In many settings, such as precision medicine, multi-study prediction is motivated by systematic data sharing and data curation initiatives. For example, the establishment of gene expression databases such as Gene Expression Omnibus (Edgar et al., 2002) and ArrayExpress (Parkinson et al., 2010) and neuro-imaging databases such as OpenNeuro (Gorgolewski et al., 2017) has facilitated access to sets of studies that provide comparable measurements of the same outcome and predictors. Comparability can often be improved through preprocessing (Lazar et al., 2012;Benito et al., 2004). For problems where such a set of studies is available, it is important to systematically integrate information across datasets when developing prediction models.
One common approach is to merge all of the datasets and treat the observations as if they are all from the same study (Xu et al., 2008;Jiang et al., 2004;Gupta et al., 2020). The resulting increase in sample size can lead to better performance when the datasets are relatively homogeneous. Also, the merged dataset is often representative of a broader reference population than any of the individual datasets. Xu et al. (2008) showed that a prognostic test for breast cancer metastases developed from merged data performed better than the prognostic tests developed using individual studies. Zhou et al. (2017) proposed hypothesis tests for determining when it is beneficial to pool data across multiple sites for linear regression, compared to using data from a single site.
Another approach is to combine results from separately trained models. Meta-analysis and ensembling both fall under this approach. Meta-analysis combines summary measures from multiple studies to increase statistical power (Riester et al., 2012;Tseng et al., 2012;Rashid et al., 2019). A common combination strategy is to take a weighted average of the study-specific summary measures. In fixed effects meta-analysis, the weights are based on the assumption that there is a single true parameter underlying all of the studies, while in random effects meta-analysis, the weights are based on a model where the true parameter varies across studies according to a probability distribution. When learners are indexed by a finite number of common parameters, meta-analysis applied to these parameters can be used for multi-study prediction, with useful results (Riester et al., 2012). Various studies have compared meta-analysis to merging. For effect size estimation, Bravata & Olkin (2001) showed that merging heterogeneous datasets can lead to spurious results while meta-analysis protects against such problematic effects. Taminau et al. (2014) and Kosch & Jung (2018) found that merging with batch effect removal had higher sensitivity than meta-analysis in gene expression analysis, while Lagani et al. (2016) found that the two approaches performed comparably in reconstruction of gene interaction networks. Ensemble learning methods (Dietterich, 2000), which combine predictions from multiple models, can also be used to leverage information from multiple studies. Ensembling can lead to lower variance and higher accuracy, and is applicable to more general classes of learners than meta-analysis. Patil & Parmigiani (2018) proposed multi-study ensembles, which are weighted averages of prediction models trained on different studies, as an alternative to merging. They showed empirically that when the datasets are heterogeneous, multi-study ensembling can lead to improved generalizability and replicability compared to merging and meta-analysis.
In this paper, we provide explicit conditions under which merging outperforms multi-study ensembling and vice versa. We study merged and ensemble learners based on ridge regression with potential basis expansion. We hypothesize a flexible mixed effects model for heterogeneity and show that merging has lower prediction error than ensembling when heterogeneity is low, but as heterogeneity increases, there exists a transition point beyond which ensembling outperforms merging. We characterize this transition point analytically under a fixed weighting scheme for the ensemble. We also discuss optimal ensemble weights for ridge regression. We study the transition point via simulations and compare merging and ensembling in practice, using microbiome data.

PROBLEM DEFINITION
We will use the following matrix notation: is a diagonal matrix with u along its diagonal, and (A) ij is the entry in row i and column j of matrix A. Other notation introduced throughout the paper is summarized in a notation table in the supplementary materials (Table A.1).
Suppose we have K studies that measure the same outcome and the same P predictors, and the datasets have been harmonized so that measurements across studies are on the same scale. For individual j in study k, let Y jk be the outcome and X jk ∈ R p the vector of predictors. Let X ijk be the value of the ith predictor for individual j in study k. Let N k be the number of individuals in study k, Assume the data are generated from the mixed effects model The random effects are independent of the residual errors. For j = 1, . . . , Q, if σ 2 j > 0, then the effect of the corresponding predictor differs across studies, and if σ 2 j = 0, then the predictor has the same effect in each study. Model 1 is nonparametric in the fixed effects but linear in the random effects. Similar models have been considered by Wang (1998) and Gu et al. (2005). The linear mixed effects model where β ∈ R P is the vector of fixed effects is a special case of Model 1. The relationship between the predictors and the outcome in a given study can be seen as a perturbation of the population-level relationship described by f . The degree of heterogeneity in predictor-outcome relationships across studies can be summarized by the sum of the variances of the random effects divided by the number of fixed effects: σ 2 = tr(G)/P . We are interested in comparing the performance of two multi-study prediction approaches as σ 2 varies: 1) merging all of the studies and fitting a single linear regression model, and 2) fitting a linear regression model on each study and forming a multi-study ensemble by taking a weighted average of the predictions. These multi-study strategies can be applied to any type of study-specific model, including non-parametric machine learning models, but our theoretical analysis will be restricted to linear models, as they have a closed form solution. Comparisons involving machine learning models without a closed form solution are also interesting and will be explored in the simulations.
When there is a single predictor, the fixed effect estimator for the linear mixed effects model is a weighted average of the study-specific least squares estimators, which is a form of ensembling. In general, ensembling linear models with appropriate weights may provide a reasonable approximation of a mixed effects model, but ensembling is applicable to a broader range of scenarios than linear models. There is also a connection between linear mixed effects models and ridge regression: the best linear unbiased predictions from a linear mixed effects model are equivalent to the predictions from a ridge regression model where the predictors with fixed effects are not penalized while the predictors with random effects are penalized (de Vlaming & Groenen, 2015).
Learners. We evaluate the performance of merged and ensemble versions of ridge regression with basis expansion. This approach encompasses many commonly used linear models as special cases, including ordinary least squares and penalized spline regression. Though we focus on a class of linear prediction models, with an appropriate choice of basis these models can in theory approximate any sufficiently smooth function to any level of precision (Stone, 1948).
Consider a pre-specified set of basis functions for each predictor: suppose there are M p one-dimensional basis functions for predictor p, given by {h 1p , . . . , h Mpp }. Let h(X jk ) = (h 11 (X 1jk ), . . . , h M 1 1 (X 1jk ), . . . , h 1P (X P jk ), . . . , h M P P (X P jk )). Define the basis-expanded design matrixX Let the first column ofX k be a column of 1's for the intercept if the prediction model will include an intercept. For the merged learner, we fit the model where λ ≥ 0 is a pre-specified regularization parameter and I − M is defined as diag(0, 1, . . . , 1) if there is an intercept and I M otherwise.
For the ensemble learner, we first fit the models separately in studies k = 1, . . . , K to obtain estimatorŝ Then the ensemble estimator, based on pre-specified weights w k satisfying K k=1 w k = 1, iŝ Examples of basis functions for a given predictor include h(x) = x, which gives a model linear in the original predictor, h m (x) = x m for m = 0, . . . , S, which gives a degree S polynomial regression model, and basis splines of degree S with knots ξ 1 , . . . , ξ T : Setting λ = 0 is equivalent to applying ordinary least squares with design matrixX. To make progress analytically, we assume that the weights w k and the regularization parameters λ and λ k are predetermined -for example, using independent data (Tsybakov, 2014). When Model 1 is a linear mixed effects model and there is a single predictor, the fixed effect estimator for this model is a weighted average of the study-specific least squares estimators, so ensembling has a similar form as estimates of the linear mixed effects model. In general, ensembling linear models provides a reasonable approximation of the mixed effects model, but ensembling is applicable to a broader range of scenarios than linear models. For any number of predictors, the best linear unbiased predictions from a linear mixed effects model are equivalent to the predictions from a ridge regression model where the predictors with fixed effects are not penalized while the predictors with random effects are penalized de Vlaming & Groenen (2015). Though the generating model is a mixed effects model, simply fitting a linear mixed effects model to the data may not be ideal. Since Model 1 is nonparametric in the fixed effects, using linear mixed effects methods may result in model misspecification. Also, mixed effects approaches can be computationally intensive when the number of random effects is large.
For linear regression, averaging predictions across study-specific learners is equivalent to averaging the estimated coefficient vectors across study-specific learners and then computing predictions. Multi-study ensembling based on linear regression on scaled variables has points in common with meta-analysis of effect sizes. In the univariate case, a standard meta-analysis is also a weighted average of the study-specific estimates, though typically weights reflect precision of estimates rather than cross-study features of predictions. When P > 1,β ens weights each dimension of the coefficient vector equally in a given study while meta-analytic approaches, which involve either performing separate univariate meta-analyses for each predictor or performing a multivariate meta-analysis (for example, see Jackson et al. 2010Jackson et al. , 2011, do not impose this constraint.
Performance Comparison. For a new test dataset, not used for training, with fixed design matrix X 0 and unknown outcome vector Y 0 , the goal is to identify conditions under which multistudy ensembling has lower mean squared prediction error, conditional on X 0 , than merging, i.e.
where the expectations are taken with respect to Y 0 |X 0 and · 2 denotes the L2 norm. Given a vector (u 1 , . . . , u m ) of length m, the L2 norm (u 1 , . . . , u m ) 2 is m i=1 u 2 i . Predictions are utilizing X 0 only, while Y 0 is held as the gold standard.

Overview
We consider two cases for the structure of G: equal variances and unequal variances. Let σ 2 (1) , . . . , σ 2 (D) be the distinct values on the diagonal of G and let J d be the number of random effects with variance σ 2 (d) . In the equal variances case where r = 1 and σ 2 j = σ 2 for j = 1, . . . , Q, Theorem 1 provides a necessary and sufficient condition for the ensemble learner to outperform the merged learner. In the unequal variances case, Theorem 2 provides sufficient conditions under which the ensemble learner outperforms the merged learner and vice versa. These conditions allow us to characterize a transition point in terms of the heterogeneity measure σ 2 between a regime that favors merging and a regime that favors ensembling. We also provide optimal ensemble weights in Proposition 1, state special cases of the results for prediction models trained using ordinary least squares in Corollaries 1-3, provide an asymptotic version of the transition point for least squares in Corollary 4, and discuss how to interpret the results. LetR the bias of the merged predictions in the test set.

Main Results
Theorem 1 is a special case of Theorem 2, which is proved in Appendix A in the supplementary materials. By Theorem 1, for any fixed weighting scheme that satisfies Condition 9, τ represents a transition point from a regime where merging outperforms ensembling to a regime where ensembling outperforms merging. The transition point depends on the true populationlevel relationship f between the predictors and the outcome since f is part of the bias terms b ens and b merge , so in general, we need to estimate f in order to estimate the transition point. However, if the models provide unbiased predictions, then the dependency on f vanishes. For example, this occurs when the generating model is linear and the prediction models are correctly specified least squares models, which are discussed in Section 3.3. THEOREM 2.

Suppose
and define Then Then Theorem 2 is proved in Appendix A in the supplementary materials. It shows that in the more general scenario where the random effects do not necessarily have the same variance, there is a transition interval such that the merged learner outperforms the ensemble learner when σ 2 is smaller than the lower bound of the interval and the ensemble learner outperforms the merged learner when σ 2 is greater than the upper bound of the interval. PROPOSITION 1. Let Then the optimal weights for the ridge regression ensemble are In general, the optimal weights depend on the true population-level relationship between the predictors and the outcome, f , via the bias terms b k . The optimal weights also depend on G and σ 2 . Using a first-order Taylor approximation, the optimal w k is approximately proportional to where the first two terms depend on the inverse mean squared prediction error and the third term depends on the covariances between the bias terms from study k and the bias terms from other studies, as well as the magnitudes of the studies' prediction errors. If the magnitudes of the prediction error and bias term for study k are held fixed and the other studies are held fixed, then the optimal weight for study k increases as it becomes less correlated with the other studies. This is similar to a result from Krogh & Vedelsby (1995) that decomposes the prediction error of an ensemble into a term that depends on the prediction errors of the individual learners and a correlation term that quantifies the disagreement across the individual learners; the decomposition implies that if the individual prediction errors are held fixed, then the performance of the ensemble improves as the correlation term decreases. The transition point under any fixed weighting scheme provides an upper bound for the optimal weights transition point, which can be calculated numerically. In practice, f , G, and σ 2 are not known, so the optimal weights need to be estimated. Weight estimation increases the variance of the ensemble learner and may result in a higher transition point than when the optimal weights are known. However, as seen in Figure A.9 in the supplementary materials, our simulations suggest that if f , G, and σ 2 can be reasonably estimated -for example, via a mixed effects model -then the estimation will have little impact on the transition point.

Special Case: Least Squares
Least squares is a special case of ridge regression where the regularization parameter is set to 0. The lack of regularization allows for a simplification of the results from Section 3.2. In particular, when the generating model is Model 2, which is linear in X k , a least squares model based on the original predictors X k is correctly specified and unbiased. Here, we present Corollaries 1-3, which are simplified versions of Theorem 1, Theorem 2, and Proposition 1 in this setting, as well as Corollary 4, which provides an asymptotic version of Corollary 2 as the number of training studies tends to infinity. In this setting, assumeX k = X k . Define Assuming R k is not identical for all k, if equal weights w k = 1/K are used, then Condition 16 is satisfied due to Jensen's operator inequality (Hansen & Pedersen, 2003) and τ LS > 0, so the transition point always exists.
COROLLARY 2. Suppose f is a linear function of X k and λ = λ k = 0.

Suppose
and define Then and define COROLLARY 3. Suppose f is a linear function of X k and λ k = 0. Then the optimal weights for the least squares ensemble are In this setting, the optimal weight for study k is proportional to the inverse mean squared prediction error of the least squares learner trained on that study.
COROLLARY 4. Suppose f is a linear function of X k and λ = λ k = 0, and there exist positive definite matrices A 1 , A 2 , A 3 ∈ R p×p such that as K → ∞, where → denotes almost sure convergence. Let w k = 1/K.

If
2. If Corollary 4 presents an asymptotic version of Corollary 2 as the number of studies goes to infinity. If all study sizes are equal to n and the predictors are independent and identically distributed within and across studies, then In the special case where P = Q = 1 and the predictor follows N (0, v), τ LS,1 = τ LS,2 converges to so asymptotically the transition point is controlled simply by the variance of the residuals, the variance of the predictor, and the study sample size.

Interpretation
The covariance matrices of linear regression coefficient estimators can be written as a sum of two components, one driven by between-study variability and one driven by within-study variability. For example, when all predictors have random effects, the covariance of the least squares ensemble estimator for the coefficients is and the covariance of the merged least squares estimator is Since the merged learner ignores between-study heterogeneity, the trace of its first component is generally larger than that of the ensemble learner. However, since the merged learner is trained on a larger sample, the trace of its second component is generally smaller than that of the ensemble learner. When merged and ensemble learners are trained using correctly specified least squares models, they are unbiased, so the transition point depends on the trade-off between these two components. When P = 1, Expression 27 shows that having a higher-variance predictor favors ensembling over merging, since increasing the variance of the predictor amplifies the impact of the random effect. Unlike least squares estimators, ridge regression estimators are biased as a result of regularization even when the models are correctly specified. The transition point for ridge regression depends on the regularization parameters used on the merged and individual datasets. It also depends on the true relationship f between the predictors and the outcome through the squared bias terms in the mean squared prediction errors of the merged and ensemble learners, so an estimate of f is needed to compute the expressions in Theorems 1 and 2. These expressions can vary considerably for different choices of regularization parameters and different forms of f . We did not provide the asymptotic results for ridge regression as K → ∞ with the study sizes held constant because this scenario is not entirely fair to the ensemble learner. For P > N k and sufficiently large K, the merged learner will be in a low-dimensional setting where the number of samples exceeds the number of predictors, while the ensemble learner will remain in the highdimensional setting. As K → ∞, the bias term approaches 0 for the merged learner, assuming λ/K → 0, but not for the ensemble learner, which suggests that when K is sufficiently large, merging will always yield lower mean squared prediction error than ensembling.
In general, the transition points for least squares and ridge regression depend on the design matrix of the test set. However, the test design matrix drops out when it is a scalar multiple of an orthogonal matrix. For example, this occurs when P = 1.

SIMULATIONS
We conducted simulations to verify the theoretical results for least squares and ridge regression and to compare them to the empirical transition points for three methods for which we cannot derive a closed-form solution: lasso, single hidden layer neural network, and random forest. We also ran a linear mixed effects model. We used the R packages glmnet, nnet, randomForest and nlme for ridge regression/lasso, neural networks, random forests, and linear mixed effects models.
The simulations were based on gut microbiome data from studies available through the curatedMetagenomicData R package, which was also the data source for our data example. In particular, we used marker abundance measurements from curatedMetagenomicData studies as predictors in the simulations. We used a generating model with P = 10 predictors, out of which Q = 5 have random effects. The outcome for individual j in study k is generated by the model where h 1,0 , . . . , h 1,3 are cubic basis splines as defined in Equation 8 with a single knot at 0, the coefficients were randomly generated from N (0, 0.3), and Z j k = (X 3jk , X 4jk , X 5jk , X 6jk , X 7jk ).
We used K = 5 datasets of size N k = 50 to train the models and evaluated the mean squared prediction error in a test dataset, k = 0, of size 100. For each study k, the predictors X k were sampled a randomly selected dataset in the curatedMegatenomicData R package.
Our main simulation analysis was based on the setting of Theorem 1, where the random effects have equal variances σ 2 1 = · · · = σ 2 q = σ. We considered choices of G corresponding to different levels of heterogeneity ranging from σ 2 = σ 2 /2 = 0 to σ 2 = σ 2 /2 = 3τ R , where τ R is the theoretical transition point for ridge regression from Theorem 1. We performed another analysis based on the setting of Theorem 2, where the random effects have unequal variances. For this analysis, we considered choices of G corresponding to levels of heterogeneity ranging from σ 2 = 0 to σ 2 = 2τ R,2 , where τ R,2 is the theoretical upper bound of the transition interval from Theorem 2. For each choice of G and each study k, including the test study k=0, the outcomes were generated from X k using the following steps: 1) sample the q random effects γ k from M V N (0, G), 2) sample the residual errors k independently from N (0, σ 2 ) with σ 2 = 1, and 3) generate Y k using Model 1 and the functional form given in Equation 29.
We trained and tested the following approaches: a correctly specified linear mixed effects model and merged and ensemble learners based on least squares, ridge regression, lasso, neural networks, and random forests. For least squares, ridge regression, and lasso, we considered fitting correctly specified models, i.e. using the correct basis expansion for the predictors, in the main analysis, but also conducted a sensitivity analysis with misspecified models. For ridge regression, lasso, random forest, and neural network, the following hyperparameters were chosen using 5-fold cross-validation in a tuning dataset of size 250 that was not used for training: the regularization parameter for ridge regression and lasso, the number of predictors sampled at each split for random forest, and the number of nodes and weight decay parameter for neural network. We restricted to neural networks with a single hidden layer to limit the computational burden. The outcomes in the tuning dataset were generated under σ 2 = 0.The hyperparameters for the merged versions of the models were selected using the entire tuning dataset while the hyperparameters for the study-specific models were selected using a subset of size 50.
In addition to the simulations with K = 5 studies and diagonal G, we performed sensitivity analyses considering scenarios with misspecified least squares/ridge regression/lasso models, i.e. using the original predictors without basis expansion, K = 10 and K = 15 instead of K = 5 training studies, studies of size N k = 100 instead of N k = 50, and scenarios with correlated random effects.
The results from the main simulation analysis based on Theorem 1 are shown in Figure 1 and Figure A.1 while the results from the other scenarios are provided in Figures A.3-A.8. The performance of the models relative to the data generating model in the main simulation analysis is summarized in Figure A.1 for three values of σ 2 . When σ 2 = 0, the merged regression learners perform as well as or slightly better than the generating model and outperform the ensembles. In particular, the merged neural network performs better than the generating model. At the ridge regression transition point, all models perform similarly. Beyond the transition point, the models continue to perform similarly, with the ensembles slightly outperforming the merged learners. In the scenarios with uncorrelated random effects, i.e. the scenarios for which we derived theoretical results (Figure 1 and Figures A.2-A.6), the empirical transition points for least squares and ridge regression agree with the theoretical results regardless of whether the models were correctly specified. Also, in each scenario, the transition points of the different methods did not vary substantially except that the random forest either had no transition point because ensembling always outperformed merging (Figure 1 and Figure A.3) or had a much earlier transition point than the other methods ( Figures A.4-A.6). This may be related to the fact that a single random forest is already an ensemble of decision trees. In the scenarios with correlated random effects ( Figures A.7 and A.8), the least squares transition point occurred earlier than when the random effects were uncorrelated. In the positive correlation scenario, the ridge regression transition point also occurred earlier than when the random effects were uncorrelated, but the opposite was true in the negative correlation scenario.

METAGENOMICS APPLICATION
Growing research on associations between the gut microbiome and health-related outcomes has motivated the development of prediction models based on metagenomic sequencing data (Fu et al., 2015;Gupta et al., 2020). The ability to predict health status from microbiome samples can help guide the design of controlled intervention studies that target the microbiome, for example using diet or probiotics. In our data example, we focus on cholesterol, a strong risk factor for cardiovascular disease. Studies have shown that the microbiome is associated with blood cholesterol levels (Fu et al., 2015;Le Roy et al., 2019;Kenny et al., 2020). In particular, the metabolism of intestinal cholesterol by gut bacteria may reduce the amount of cholesterol absorbed from the intestine, resulting in lower blood cholesterol (Kenny et al., 2020). Recent work has identified bacterial genes involved in cholesterol metabolism (Kenny et al., 2020).
To illustrate in a practical example, we compared the performance of merging and ensembling in developing models to predict cholesterol levels from metagenomic sequencing data. We used datasets from the curatedMetagenomicData R package (Pasolli et al., 2017), which contains a collection of curated, uniformly processed whole-metagenome sequencing data from multiple studies. In addition to metagenomic data such as gene marker abundance, it also includes demographic and clinical data. Cholesterol measurements were available for three of the studies that sequenced gut microbial DNA: 1) a study of Chinese type 2 diabetes patients Fig. 1: Relative performance of multi-study ensembling and merging as a function of heterogeneity in the main simulation scenario where K = 5, N k = 50, P = 10, Q = 5, and the random effects have equal variances. MSPE: mean squared prediction error. The vertical dashed lines correspond to the theoretical transition points calculated using Theorem 1. The empirical transition point occurs at the value of σ 2 where the log ratio of the prediction errors for ensembling and merging is equal to 0. and non-diabetic controls, including N 1 = 151 female patients (Qin et al., 2012), 2) a study of middle-aged European women with normal, impaired or diabetic glucose control (Karlsson et al., 2013), and 3) a study of patients with a family history of type 1 diabetes (Heintz-Buschart et al., 2017). We used samples from female patients in these studies. The sample sizes were 151, 145, and 32 for the three studies.
We used merging and ensembling to train least squares, ridge regression, lasso, neural network, and random forest models to predict cholesterol, calculated the theoretical transition intervals for least squares and ridge regression, and evaluated the performance of the approaches. We considered two scenarios: 1) training on K = 4 different subsets of the same study and testing on a held out subset, and 2) training on K = 2 different studies and testing on an independent study. In the first scenario, we randomly split the samples from Qin et al. (2012) into five datasets of approximately equal size, using K = 4 for training and the remaining one for testing. In second scenario, we used the K = 2 datasets from Qin et al. (2012) and Karlsson et al. (2013) for training and dataset from Heintz-Buschart et al. (2017) for testing. We used age and gene marker abundance for a selected set of gene markers as the predictors. Due to the limited sample sizes, we restricted to a small number of gene markers when fitting the prediction models. In the first scenario, we selected the top five gene markers that were most highly correlated with cholesterol in the training set and in the second scenario, we selected the top twenty gene markers that were most highly correlated with cholesterol in the training set. See Appendix C in the supplementary materials for a list of the selected gene markers. In each scenario, we fit merged and ensemble versions of least squares, ridge regression, lasso, neural network, and random forest. The linear models were specified as Y jk = β T X jk + jk where X jk contained 7 elements in the first scenario -intercept, age, and 5 gene markers -and 22 elements in the second scenario -intercept, age, and 20 gene markers. To estimate the transition intervals, we estimated G by fitting a linear mixed effects model using residual maximum likelihood, allowing each predictor to have a random effect: Y jk = β T X jk + γ k T X jk + jk . We estimated the optimal weights for least squares and ridge regression using Proposition 1 and estimates of G, σ 2 , and β from the mixed effects model. We calculated the theoretical transition bounds from Theorem 2 assuming a linear generating model and compared them to the estimate of σ 2 . We used the optimal weights for ridge regression to ensemble the study-specific lasso, neural network, and random forest models. We evaluated the performance of merging and ensembling by calculating the prediction error on the test set.
In the first scenario, σ 2 was estimated to be 0.005 2 , the lower bound of the transition interval for least squares was τ LS,1 = 0.266 2 , and the lower bound of the transition interval for ridge regression was τ R,1 = 0.050 2 , so merging was expected outperform ensembling. In the test set, the merged versions of all of the models except lasso had lower prediction error than the corresponding ensembles (Figure 2). In the second scenario, σ 2 was estimated to be 18.322 2 , the upper bound of the transition interval for least squares was τ LS,2 = 15.247 2 , and the upper bound of the transition interval for ridge regression was τ R,2 = 4.423 2 , so ensembling was expected to outperform merging. In the test set, the ensemble versions of least squares, ridge regression, and neural network had lower prediction error than the corresponding merged learners, while merging outperformed ensembling for lasso and merging performed similarly to ensembling for random forest (Figure 3). Lasso and random forest both had considerably higher prediction error than the other models.

DISCUSSION
The availability of large though somewhat heterogeneous collections of data for training classifiers is challenging traditional approaches for training and validating prediction and classification algorithms. At the same time, it is creating opportunities for new and more general paradigms. One of these is multi-study ensemble learning, motivated by variation in the relation between predictors and outcomes across collections of similar studies. A natural benchmark for these methods is to combine all training studies to exploit the power of larger training sample sizes. In previous work (Patil & Parmigiani, 2018), merged learners were shown to perform better than ensemble learners in simulations under low-heterogeneity settings. As heterogeneity increased, however, there was a transition point in the heterogeneity scale beyond which acknowledging cross-study heterogeneity became preferable, and the ensemble learners outperformed the merged learners.
In this paper, we provide the first theoretical investigation of multi-study ensembling. Using a flexible data generating framework based on a nonparametric mixed effects model, we have been able to prove that when prediction is done using linear models, multi-study ensembling can, under certain conditions, achieve better optimality properties than merging all studies prior to training, a question that was still open. The class of linear models we have considered includes regularized and basis-expanded models, which are able to provide reasonable approximations of potentially nonlinear data generating models. We characterized the relative performance of merging and ensembling as a function of inter-study heterogeneity, explicitly identifying the transition point beyond which ensembling outperforms merging for least squares and ridge regression. These results are applicable to both correctly specified and misspecified linear prediction models. We confirmed the analytic results in simulations and demonstrated that when the data are generated by a mixed effects model, the least squares and ridge regression solutions can potentially serve as proxies for the transition point under other learning strategies (lasso, neural network) for which closed-form derivation is difficult. Finally, we estimated the transition point in cases of low and high cross-study heterogeneity in microbiome data and showed how it can be used as a guide for deciding when and when not to merge studies together in the course of learning a prediction rule.
We focused on deriving analytic results for ridge regression and its variants, including least squares and penalized spline regression, because of the opportunity to pursue closed-form solutions. Other widely used methods such as lasso, neural networks, and random forests are not as easily amenable to a closed-form solution, so we used simulations to investigate whether the analytic results we derived could potentially serve as an approximation for these methods. In our simulation settings, the empirical transition points of least squares, ridge regression, lasso, and neural networks fell within a similar range, but the random forests showed very different behavior. In one setting, the ensemble version of the random forest always outperformed the merged version and there was no transition point, while in other settings, the transition for the random forest occurred much earlier than for the other models. A key difference between random forests and the other methods considered is that a single random forest is already an ensemble of decision trees, which potentially complicates the comparison of merging to ensembling. Further exploration of this issue is provided by Ramchandran et al. (2020). The simulation results suggest that in some scenarios it may be reasonable to use the transition points for ridge regression or least squares as a proxy when considering machine learning models such as lasso or neural networks. However, the results are not generalizable to random forests and our data example shows that lasso can sometimes behave very differently from ridge regression and least squares, so caution is needed in using such approximations. Also, it is important to consider how the reliability of such approximations could be affected by the choice of model hyperparameters.
In practice, the analytic transition point and transition interval expressions could be used to help guide decisions about whether to merge data from multiple studies when there is potential heterogeneity in predictor-outcome relationships across the study populations. σ 2 can be estimated from the training data and compared to the theoretical transition points or bounds for least squares and/or ridge regression. Various methods can be used to estimate σ 2 , including maximum likelihood and method of moments-based approaches used in meta-analysis (for example, see Jackson et al. 2016), with the caveat that estimates will be imprecise when the number of studies is small.
Under Model 1, fitting a correctly specified mixed effects model will generally be more efficient than both the merged and ensemble versions of least squares. However, more flexible machine learning algorithms can potentially yield better prediction accuracy than the true model. For example, in the main simulation analysis, the mixed effects model was outperformed by a neural network. Moreover, fitting a mixed effects model can be computationally difficult when the number of predictors is large and standard mixed effects models are not appropriate for highdimensional data, though there are methods for penalized mixed effects models (Bondell et al., 2010;Ibrahim et al., 2011;Schelldorfer et al., 2011;Fan & Li, 2012).
A limitation of our derivations is that they treat the following quantities as known: the subset of predictors with random effects, the ensemble weights, and the regularization parameters for ridge regression. In practice, these are usually selected using statistical procedures that introduce additional variability. Furthermore, we obtained closed-form transition point expressions for cases where the ensemble weighting scheme does not depend on the variances of the random effects. Such weighting schemes are generally sub-optimal, as the optimal weights given by Proposition 1 depend on G. Therefore, the closed-form results are based on a conservative estimate of the maximal performance of multi-study ensembling. Another limitation is the assumption that the random effects are uncorrelated, which is often not true in practice.
In summary, although this work is predicated upon the assumption that cross-study heterogeneity manifests as random effects and assumes that weights and regularization parameters are known, we believe it provides a theoretical rationale for multi-study ensembling and a strong foundation for developing practical rules and guidelines to implement it.

APPENDIX A: CALCULATIONS AND PROOFS
A.1. Covariance Matrices and MSPEs The covariances and expectations are all conditional on the predictors.
A.2. Proof of Theorem 2 Theorems 1 is a special case of Theorem 2. Let a d = i:σ 2 i =σ 2 Using Lagrange multipliers, we get the system of equations The solution to the system is where the entries of C are (C) kk = (v k + b k T b k ) and (C) jk = b j T b k for j = k and j = (1, . . . , 1) T ∈ R K .
Least Squares Under a Linear Generating Model. For least squares, we want to minimize K k=1 w 2 k tr(GZ 0 subject to the constraint k w k = 1. Using Lagrange multipliers, we get the system of equations 2w k tr(GZ 0 T Z 0 ) + σ 2 tr(R k −1 R 0 ) = α (for k = 1, . . . , K; α is the Lagrange multiplier) k w k = 1 which is solved by Plugging the optimal weights into the MSPE of the ensemble learner, tr(Cov(X 0βLS,ens )), we get tr(Cov(X 0βLS,ens )) = 1 k tr(GZ 0 In the equal variances setting, tr(Cov(X 0βLS,ens )) = 1 k σ 2 tr(Z 0 T Z 0 ) + σ 2 tr(R k −1 R 0 ) −1 so tr(Cov(X 0βLS,ens )) is approximately linear in σ 2 when tr(R k −1 R 0 ) is similar across studies or when σ 2 is large.       Relative performance of multi-study ensembling and merging as a function of heterogeneity when K = 5, N k = 50, and the random effects have equal variances σ 2 and are positively correlated with correlation 0.25σ 2 . MSPE: mean squared prediction error. The vertical dashed lines correspond to the transition points calculated using Theorem 1 under uncorrelated random effects. Fig. A.8: Relative performance of multi-study ensembling and merging as a function of heterogeneity when K = 5, N k = 50, and the random effects have equal variances σ 2 and are negatively correlated with correlation −0.25σ 2 . MSPE: mean squared prediction error. The vertical dashed lines correspond to the transition points calculated using Theorem 1 under uncorrelated random effects.