Multi-task learning with summary statistics
Authors: Parker Knight
Rui Duan
Venue: NeurIPS 2023 poster (2023)
Keywords: multi-task learning
genetic risk prediction
summary statistics
TL;DR: We develop a multi-task learning framework for integrating summary statistics from distinct sources to form genetic risk predictions.

Abstract: Multi-task learning has emerged as a powerful machine learning paradigm for integrating data from multiple sources, leveraging similarities between tasks to improve overall model performance. However, the application of multi-task learning to real-world settings is hindered by data-sharing constraints, especially in healthcare settings. To address this challenge, we propose a flexible multi-task learning framework utilizing summary statistics from various sources. Additionally, we present an adaptive parameter selection approach based on a variant of Lepski's method, allowing for data-driven tuning parameter selection when only summary statistics are accessible. Our systematic non-asymptotic analysis characterizes the performance of the proposed methods under various regimes of the source datasets' sample complexity and overlap.  We demonstrate our theoretical findings and the performance of the method  through extensive simulations. This work offers a more flexible tool for training related models across various domains, with  practical implications in genetic risk prediction and many other fields.

Multi-Task Learning with Summary Statistics
Parker Knight Department of Biostatistics Harvard University
Rui Duan
Department of Biostatistics Harvard University Boston, MA aaaaanhhshhhhaavvee edu
Boston, MA pknightēg. hwrvard. edu
Abstract
Multi-task learning has emerged as a powerful machine learning paradigm for integrating data from multiple sources, leveraging similarities between tasks improve overall model performance. However, the application of multi-task learn- ing real-world settings is hindered by data-sharing constraints, especially in healthcare settings. To address this challenge, we propose a flexible multi-task learning framework utilizing summary statistics from various sources. Additionally, we present an adaptive parameter selection approach based on a variant of Lepski's method, allowing for data-driven tuning parameter selection when only summary statistics are available. Our systematic non-asymptotic analysis characterizes the performance of the proposed methods under various regimes of the sample com- plexity and overlap. We demonstrate our theoretical findings and the performance of the method through extensive simulations. This work offers a more flexible tool for training related models across various domains, with practical implications in genetic risk prediction and many other fields.
Introduction
The growing availability of extensive and intricate datasets presents opportunity to integrate data from multiple sources. Multi-task learning has emerged as a promising machine learning approach that enables the simultaneous learning of multiple related models, leveraging shared structure between tasks to enhance the performance on each task individually [36, 35]. In healthcare and biomedical research, the practical application of multi-task learning is often hindered by data-sharing constraints, which stem from concerns about the ownership and privacy of individual-level data [32, 8J. Patient data in these domains typically sensitive and less likely be publicly available or shared across study sites, limiting researchers' access to individual-led data from different domains.
To overcome this limitation, researchers have increasingly integrated summary statistics into analysis pipelines as a substitute for individual-level data [22, 13, 12]. Summary statistics are straightforward, interpretable measures derived from raw data that offer insights into data distribution, variability, and relationships among variables. Furthermore, they can be aggregated across studies to facilitate data integration and reused in various research projects. Recently, the use of summary statistics has garnered interest in healthcare and biomedical research. For example, many genetic risk prediction methods rely on summary-lewe level statistics such as associations from Genome-wide Association Studies (GWAS), Linkage Disequilibrium estimations (LD), and minor allele frequencies (MAFs) [6]. These summary statistics can help predict an individual's likelihood of developing specific diseases based on their genetic profile.
Inspired by a potential use case genetic risk prediction, we propose multi-task learning framework that enables simultaneous learning of multiple genetic risk prediction models using only publicly available summary statistics. Our proposed framework can be used in the context of predicting genetic risks for multiple traits leveraging potentially shared genetic pathways, and can also be used
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

to develop trans-ethnic genetic risk prediction models that account for potential heterogeneity across populations, improving generalizability and real- world applicability. Beyond genetic risk prediction. the ability to learn from summary statisties offers a versatile tool for developing models across a wide range of domains, including healtheare, finance, and marketing-
To summarize, the contributions of this work are threefold: First, we propose a flexible multi-task learning framework which allows training multiple models simultaneously using basic summary statistics characterizing marginal relationship between outcomes and features, which are often publicly available. We allow summary statistics corresponding each task to be generated from distinct or potentially overlapping samples. Secondly, we conducted a systematic characterizes how the performance of the proposed methods are influenced by the characteristics of summary statistics. In particular, we show that there are multiple regimes of performance depending on the sample complexity of the source datasets and their overlap. The theoretical results are supported with extensive simulations. Lastly, We propose an adaptive scheme for tuning parameter selection based on the variant of Lepski's method [14] given in [5]. This allows us to select a data-driven tuning parameter when only summary statistics are available and cross-validation is not feasible. We prove that tuning parameters chosen by this method satisfy an oracle inequality with high probability, and demonstrate the effectiveness of the method via simulations.
-ayympttttiic analysis which
1.1 Related work
The use of summary statistics for regression modeling has been considered in the statistical genetics literature [6]. The lassosum method for polygenic risk prediction was introduced by [22]. which considered fitting a Ln penalized linear regression with summary statistics, and its theoretical properties were studied in depth by [15]. In [4], the authors extend these ideas to polygenic risk prediction with binary traits. The summary statistics used in these methods include the marginal associations between genetic variants and phenotypes, and statistics summarizing the covariance structures among all genetic variants oftentimes derived from a reference genotype dataset. Empirical studies have demonstrated that the efficacy of such models is significantly influenced by the choices of the GWAS summary statistics and the reference dataset [30]. However, there's still limited theoretical understanding regarding how the overlap of samples and the inherent heterogeneity between datasets impact the model performance. Moreover, most current approaches devise models for a single trait within a single ancestral population. Considering shared genetic architectures could potentially enhance performance by employing multi-task learning strategy [23].The authors of [24] take this approach, and describe a multi-task estimator for multi-ancestry PQTL analysis. However, they do not consider the setting when only summary statistics are available for each task.
Our methods build upon classical multi-task learning techniques, and enable fitting models only using basic summary statistics which are often made publicly available. The sparse regularized estimator extends the group-sparse estimators studied in [19, 18], while the nuclear norm estimator expands on the low-rank regression model described in [25]. The nuclear norm approach is closely related to the linear representation learning problem [7, 31], which constrains the regression coefficients to a shared oiwdeimesio subspace.
Another recent line of work studies the multi-task learning problem under data Mata-sharing constraints. In [17], the authors describe a federated multi-task learning linear regression model for privacy- preserving data analysis. Similarly, [2] presents computational framework for multi-task learning under DataSHIELD [1]] constraints. The formulation of these methods is conceptually similar to ours, but they do not provide theoretical guarantees for their estimators, and we consider a more flexible setting where the summary statistics can be derived from different sources.
Finally, our methods are elosely related to the one-shot federated learning paradigm, in which only one round of communication is permitted between the primary local research site and additional sites. [9] presents an algorithm for fitting logistic regression models using summary statistics from different research sites. The works of [20, 34] extend these ideas to linear mixed effects models and generalized mixed effects models, respectively. [21] presents a federated algorithm for fitting the Cox proportional hazards model, and [16] studies federated transfer learning methods for fitting generalized linear models. Nevertheless, the summary statistics addressed in our research are frequently reported in existing studies and can be employed across various models. This is in contrast to the one-shot federated algorithm, where the summary statisties are model specific and the implementation relies on the infrastructure of collaborative environment.

---

2 Problem setup and methods
Consider the setting where we are interested in learning a total Q tasks simultaneously. For each q E [Q], we posit the linear model
Y x g(9) + e(0))
where Y(4) Ree,, X(o) R" and g(9) is mean-zero random noise. Each index q corresponds to the qth task. The dataset D(a) (Y ( contains the individual-level observations of the outcome and features respectively for the Qta task. We consider the generic setting where the features D(g) might be collected from either overlapping or non-overlapping samples across tasks. Our estimand of interest is the matrix B = [B(1) ..... {(2) E R where the qth column of B is B(0),
Furthermore, let e, denote the ith standard basis vector, so that , eg.
If all the individual-level observations D(g) are available, natural estimator of B is the regularized multi-task least-squares estimator
B arg arg min B
1 Y(o) 2mg yEIQ)
| X
x(() Be
|2 XP
(I)
where P is a suitable penalty, chosen to enforce similarity structure between tasks, with tuning parameter À >> O. However, in many applications, we are less likely to observe D(a), Rather, summary statistics which contains information of the feature-outcome and feature-feature relationships may be
more likely to be made publicly available. Motivated by the use case in genetic risk prediction, we assume that only summary statistics ((9) and (4) are observable, where S(4)
1 (()) (4) are ny
derived from (X(9), Y( which we termed as the discovery data, and ((4) =
1 X (4) is
sample covariance matrix computed from the proxy data (4) E RAAXP, which may may not have overlap with X (t)
Our goal is to estimate B" using two sets of summary statistics S(1) and ((4), We note that ( is not necessarily equal to X In practice, the studies which report S(9) may not be the same as the ones reporting (g) Intuitively, we hope that (() ()) is generated from similar population as ( but this may not hold in general. In Section 3, our theoretical analysis reveals how the overlap
between X and (1) and their distributional shift can influence the accuracy of multi-task learning.
To construct an estimator that uses only the information provided by D(4), we notice that the least- squares loss can be written as
((8) = - x - -13, x
By dropping the constant term, we arrive loss function that be computed using only summary- level information, namely the matrices X and ' This motivates our general strategy for constructing an estimator only using summary statistics: we substitute S( and E(0) where ap- propriate in each least-square loss function in Equation and arrive the following optimization problem.
B
arg min B
S Beg -S((1) Beg) XP(B)
2
VEIQI
There are many possible choices of P for enforcing structure similarities across tasks. For instance, the recent works of [29] and [12] study low-rank and angle-based penalties for enforcing a shared orientation among the task-specific parameters. In this work, we study two estimators obtained under the 1,, norm penalty, denoted ||-1l,, and the nuclear norm penalty, denoted ||-I|_. These penalties are chosen for their intuitive interpretation: the f,1 penalty is more likely to effective a common set of variables are active across the tasks. If the task-specific parameters tend to be "correlated",

---

in the sense that they lie in a low-dimensionss subspace, the nuclear norm penalty is preferred. In practice, certain domain knowledge can be incorporated determine the penalty structure, or it can be chosen in a data-driven way in the existence of a validation dataset.
The corresponding estimators expressed as follows:
((xp) arg min
(1)
1/21/22 Beg
-(((1) BBe) + A |
[B
(3)
9E[Q1
B(fr) = arg min
Σ (E(Q)
(4)1
BBg
Beg
(S(i),Beg)) + A B
(4)
2
B
The superscripts (sp) and (Ir) stand for "sparse" and "low-rank" respectively.
3 Theoretical guarantees
Before presenting our theoretical results, we first introduce the relevant notation. Let N denote the total size of discovery observations and proxy observation across Q tasks. Formally,
N
Σ ny + Rgg) 9=1
We note that N may double-count individuals who are part of both the proxy data and the discovery data. Define the subset Ig C [T] as the index for the discovery data points the quh task; in other
words E Ig implies X, E R is a row of ((1), We define I, analogously for the proxy data; E Ig implies X; is row of (a), Let Pa Iy n I/l denote the proportion of proxy samples which are also in the discovery dataset for the qua task. In the results that follow, let
74 1+ 11++|30001
ng "g + 27g) flg
and take 7 = max, Tg- Additionally, let E\t RPxQ be the matrix with its qth column equal to (Σ(9) Σ )(((), where Σ(1) and x are the population-lewe covariance matrices of (() and ( respectively. The quantities y and E play important roles in our results that follow. In particular, is a multiplicative factor in our bounds that represents the cost of using proxy data rather than individual- data. Similarly, E will represent the cost of using a proxy dataset with a distributional shift from the discovery data. Finally, we will let Memi and Amin denote the smallest sample size of discovery and proxy data, respectively. All proofs are given in the supplement.
3.1 Guarantees for f2,1 -nor estimator
In this section, formally state our assumptions and results for the B (sp) estimator. The assumptions are standard for nigh-dimensional regularized estimators, see [26] for a deeper discussion these conditions.
Assumption 3.1 "Sub-gaussian design and noise). The following holds for each [Q]: The rows of X are independent and identically distributed according to sub-Gaussian distribution with covariance matrix x(9) Rpxp, Similarly, the rows of (4) are independent and identically distributed according to Sub-Gaussian distribution with covariance Σ RPX,, The matrices Σ(() and Σ( have bounded eigenvalues. The entries of cla) are independent and identically distributed according to a uub-Gaussian distribation with parameter o². The X(0 and el are independent of one another.

---

Assumption 3.2 (Shared support). There exists a subset S C [p] such that supp( = S" for each q.
Definition 3.1 (Sparse cone). For any C [p], let
Ca( RoxQ 11, [|Asl121]
Assumption 3.3 (Restricted strong convexity). There exists constant and a sequence aN as N + 20 such that the following inequality holds for each A E C2(S*) with probability at least 1 oN.
Σ 2011222
-A
Theorem 3.1. Under assumptions 3.1, 3.2, and 3.3, there exist constants ei and C2 depending only on the d² and the eigenvalues of E1 and ( such that if T'in A Aumin 2 || (Q + log
and O (VQ + log p) mmin 2,00) the following inequality holds with probability at least 1 - e log ON?
(V (1re + log
) / /
B ( Bd
S S02
T'min
In the subsequent discussion, take q" arg gax "y and (n, n, P) (ng', M.: p.g)) that the triplet (n, n, p) corresponds to the same sizes and overlap factor used compute 7-
There are three main quantities in this upper bound that are novel interest: the ratio of discovery data size to proxy data size n the proportion of overlap between the discovery and proxy data p. and the error in specifying the proxy data distribution E\t (I ) The first two of these are captured by the factor 7. Our results show that with fixed n and n. the larger proportion of overlap leads to better estimation accuracy. When the proxy data and discovery data are precisely the same, meaning that I, Ig for all penalized multi-task learning problem established by Theorem 6. of [18]. If the proxy data and the discovery data are disjoint, meaning that I, n nI,, = 9 for all g, the error increased relative to minimax rate by a factor of (1 + nz/n) ||@1|2 This recovers the result of Theorem 2. in [15] up to constant factor, assuming that O. The novelty Theorem 3.1 is that we are able to characterize convergence rate (sp) for any values of n/ū, p, and E. Additionally, we emphasize that form of the term implies that price is paid anytime when (q) is not fully contained in X (4), Indeed, if أ < 1/2 and we take ñ oo we still have that y > 1 as long the signal nonzero. Counter intuite-ituuiivvel.. this indicates that an oracle model which has full access to the population-leve covariance matrix of the covariates will perform worse in terms of estimation error than an estimator which has access to individual-lert estimator will out-perform the oracle estimator that uses the population covariance matrix. These phenomena are validated in our simulation studies in Section 5.
recover the minimax rate of estimation for the &,1
data. Furthermore, if أ > 1/2, our theorem predicts that the
3.2 Guarantees for the nuclear norm estimator
Now we state our results for the low-rank proxy data estimator, when the penalty taken be the nuclear norm. Once again, these assumptions are standard for high- high-dimensional regression problems with the nuclear norm [26, 33].
Assumption 3.4 (Low rank). The matrix B* has rank r <<pA Q. Let lI and V denote the column space and row space of B" respectively. Note that Ut and V each have dimension r.
Definition 3.2 (Subspaces). Let ll denote a dimension lo S p A Q subspace of IR", and let V denote dimension i & p subspace RQ, Define
M M(U, {A RPXQ row A) V, vol()(() 4}
5

---

M- = M = { E RexQ row(4) T V. col(A) 1 2
Furthermore, for any subspace Q of BpxQ, let An denote the projection of A onto 12
We will denote M Mi4", V*).
Definition 3.3 (Low rank cone). For any set M. as defined above, let
(M (A RPxQ :1 ll. |a [|
Assumption 3.5 (Restricted strong convexity). There exists a constant and a sequence by - as N such that the following inequality holds for each A C2(M)) with probability at least by:
Q
Σ
1 |
Theorem 3.2. Under assumptions 3.1, 3.4, and 3.5, there exist constants ei and 02 depending
only o² and the eigenvalues of ((() and ΣXe) such that Imin Dmin 2 ci |B* (Q p)
and À O( (((1(Q + P)/Bmin |3 |llop) the following inequality holds with probability at least 1 P bn:
(/
+ m(( (m((+P)) + Vr F
B B c2
Amin
This theorem recovers precisely the same behavior with respect y and as Theorem 3.1. As + 1, we achieve the minimax rate of estimation for low-rank regression derived in [27] as long as E
Tuning parameter selection with Lepski's method
key challenge applying penalized regression models to summary statistics is that model tuning based on data splitting (e.g., training and validation) is no longer option. Model selection methods based on information criteria require knowing the log squared loss log || Y - (9) 8112 which cannot be recovered from S(e) and (4) [1]. To address this, we propose
tuning scheme based on Lepski's method [14], classical tool of nonparametric statistics for adaptive estimation with unknown tuning parameters. The authors of [5] apply the ideas of Lepski to the LASSO, providing a fast algorithm for model tuning with non-asymptotic guarantees. In this section, we extend the methods in [5] to tune the multi-task estimators described in the present work.
The results and ideas in this section apply to both (xp) and (Ir),, without loss of generality, let (Bx, P) denote generic estimator-regularizer pair with tuning paramter which may refer either
(B(*)) .II-l2,.) ) or (Ir) II-II.). Additionally, let D denote the dual of P, meaning that
(*(X)
sup (X,Y)
Y:P(Y)1
Finally, we let L denote the loss function for both estimators and let VL denote its gradient.
The intuition behind the adaptive tuning procedure that the tuning parameter should be chosen large enough to control fluctuations in the gradient of the loss function, but not too large such that too much bias is incurred. [26] articulates that the performance of regression estimator with a convex penalty is contingent the following event occurring with high probability:
A(A) = p' PVVVL((B)) < 1

---

where we use our problem's 's notation for continuity. The proofs of Theorem 3.1 and 3.2 involve showing that A(x) holds with high probability under our stated conditions. It is straightforward to prove the following proposition, which states that conditional on A, the gradient of the loss function generic estimator B is close to the gradient at the true parameter B*.
Proposition 4.1. Let (Bx, P) denote generic estimator estimator-reggulatzee regularizer pair. Conditional the event A()), there exists a constant C > O such that the following inequality is satisfied almost surely:
p P'((L((B)) VL(())) SA
This proposition motivates the following definition, which we adopt from [5].
Definition 4.1. Let A {A1, Aa,...,AMM ..M denote grid of potential tuning parameters ordered such that 0 < A1 < 22 < < AM < co. Fix 6 E (0, 1). The oracle tuning parameter X3 is defined as
X = arg (P (A())) 2 21-83
The oracle tuming parameter provides the tightest upper bound in Proposition 4.1, but is unknowable in practice, since do not observe B* and hence cannot verify A. The aim of our Lepski-type method is to mimic the performance of λž in an entirely data-driven fashion. Letting denote our ordered grid of potential tuning parameters, we follow [5] and choose the tuning parameter E A that satisfies
À = arg min XEA
P" P'(V[(Bx)) - > C(X + xx1 XAX
max
(5)
where C' is constant chosen by the statistician. The following theorem states that Ả recovers the behavior of Xī with high probability, as long as sufficiently large.
Theorem 4.1. Let C denote the constant in Proposition I. If X is chosen as in Equation 5 with C C, then the following inequalities hold simultaneously with probability least 1 - 6:
1. 1
2. Pd PVV(L(B)) VL(B")) < C*
where C* CC
This theorem is generalization of Theorem in [5], adapted to our setting. The primary advantage of this Lepski-style tuning scheme is that it can be performed using only the gradient of the loss function, which in setting consists only of summary-level statistics. This is marked improvement other summary statistic- based estimators, which typically require an additional set individual-level data for tuning.
The adaptive tuning scheme does hold some disadvantages. First of all, it requires a choice of constant C, which should be taken to be close to the constant in Proposition 4. as possible. Remark 10 in [5] offers some guidance as how to choose C but unfortunately their analysis corresponds only to the LASSO. Deriving the exact constant in Proposition 4.1 may be possible under stronger assumptions on the data-generating process (i.e. Gaussianity), and we view this as an area future work. Furthermore, Theorem 4. offers only a bound on ( - VL(B*). Translating this to bound on B - B" will require strong element-wise conditions each of the matrices Σ(9), which we not explore in the present work. Nevertheless, our simulations in Section indicate that the adaptive tuning method performs well in terms of the MSE B, suggesting that adaptive tuning is good option for model selection when only summary statistics are available.
5 Numerical experiments and real data application
We validate our theory and demonstrate the effectiveness multi-task learning proxy data settings via extensive experiments. In each experiment, take the proxy dataset to be well-specified;

---

in other words, we assume that E-t 0. When E is nonzero. this predictably leads to worse performance, which we demonstrate in the supplement. The code, further implementation details, and additional simulations which explore the use of our adaptive tuning procedure are also available the supplement.
First, we consider the effect of varying proxy data size on empirical MSE per task. We generate synthetic Gaussian data with Témin 100, 100, Amin 77min for {0.5, 2, 5, 10), and P = 0 for each q. The number of tasks was fixed at 8. Furthermore, we generate row-sparse B* matrix with 10 nonzero rows and a B with rank 2 for the sparse and low-rank multi-task estimators, respectively. We then fit the proxy data multi-task learning estimator and compare the prediction MSE per task to the estimator that has access all of the individual level data and the estimator that uses the true covariance matrix Σ. The results of this simulation are given in Figure 1.
Prediction MSE Lask
20
15
Estimator
Prediction MSE
Estimator Proxy Cov True_Cov
Cov Proxy Cov Thue Cov
20
0.5
0.5
10
n_tilde/
n_tilde n
Figure 1: Average prediction MSE per task after 100 repetitions plotted against T n/n. The left hand side corresponds to the sparse estimator, and the right hand side is the low-rank estimator. The red boxes indicate the estimator that uses all of the individual-level data (IL_Cov), the blue boxes indicate the estimator that uses the true covariance matrix of the features (true _ov), and the green boxes correspond to the estimator that uses just the proxy data PProy__Covi
We observe performance gap between the estimators that use the true covariance matrix and the individual level estimators, as predicted by our theory in Section 3. The performance of the proxy data estimators increase with increasing proxy sample size, but are unable to match the performance of the individual level estimator, as expected.
Next we study the effect of varying the proportion of overlapping samples between the discovery and proxy datasets. Simiarly, we generate synthetic data with n the proportion of proxy data points that are also in the discovery dataset. With Q 8, we generate B" the same way as in the previous simulation. These results are given in Figure
100, and vary ý. which indicates
Once again, we observe the expected performance gap between the estimators with the true covariance and the individual-level estimators. As the proportion of overlap between the proxy dataset and the discovery dataset grows, we see that the performance of the proxy data estimator converges to that of the individual level estimator. This is anticipated by Theorems 3.1 and 3.2.
0.35
0.07
MSE per task
0.30
MSE per task
0.08
Estimator
Estimator
Cov Praxy Cov Truc. Truc_Cov
0.25
0.05
Praxy Cov True True_Coo Cav
0.20
0.04
40
0.5...5
0.250.50.75 Proportion overlap
Proportion overlap
Figure 2: Average MSE per task after 100 repetitions plotted against P. The orientation and colors are the same in Figure!

---

Finally, we have applied our method to analyze real genetic data to demonstrate the real-world applicability of our method. We use multi-site data obtained from the electronic Medical Records and Genomics (eMERGE) network [28], which includes individual-level genotype data from multiple research sites in the United States. Our goal is to predict levels of low-density lipoprotein (LDL) across five adult sites, treating the data from each site as separate task. We split the data (with sample sizes n 3813, 122 = 546, n2 = 2666, n4 = 1435, 725 = 525) at each task into a training and test set (with test set data size of 100 for each task) and evaluate the performance of our method using the prediction MSE on the test set. The training data from each site is used to construct the
discovery summary statistics S(4) for each task. For approximating Σ(9) we choose two different approaches: is to the half of the genotype data from each site (this approach is labeled as Proxy_MTLI> MTL1); the other approach is to use X, (genotype data from site 1) to approximate E for all the sites. This approach is labeled as Proxy Proxy_MTL2, MTL2. We use these two approaches to demonstrate potential trade-off in the construction of the reference panel: Proxy MTL1 uses well well-specified reference dataset with a smaller sample size; Trox__MTI_ uses a larger reference dataset that may suffer from distribution shift. For comparison, we also fit multi-task learning estimator that uses all of the individual level training data for each task, which is labeled 'Individual _MTI',, and we fit a ridge regression estimator that models each task separately, for comparison our multi-task learning approach. The ridge estimator uses the proxy sample covariance instead of the individual level covariance matrix for fair comparison with our method. We repeat the train-test split process 10 times, which admits the distribution of prediction MSE values that we report in the figure. We use the nuclear norm penalized multi-tak task learning estimators in this application, because we believe that the genetic effects in this dataset are dense.
Prediction MSE per Task
0.70
Estimator
Individual_M
Proxy_MTL1 MTL1
Proxy_ MTL2
Pidge_Ridge
0.65
Figure 3: Prediction MSE per task after 10 splits of the eMERGE data.
Our results in Figure demonstrate that our method is highly practical when only summary-leve information is available, the prediction MSE of our method is nearly the same as the estimator which uses the individual-level data, despite a slight cost in performance. Furthermore, all multi-task learning estimators outperform the ridge-estimator, confirming that multi-task learning is strong approach when there is shared structure between tasks.
Discussion, Limitations, and Broader Impacts
We have described flexible multi-task framework incorporating summary statistics from distinct sources with general data-driven tuning scheme for selecting tuning parameters. Our theoretical analysis sheds light on the intrinsic price of using summary-level information from distinct sources for statistical analysis, and suggests that more overlap between the sources, less distributional shift, and larger proxy data sample sizes alleviate this cost. Our data-driven tuning scheme allows

---

models to be trained without sample splitting, making it more applicable to real-world settings with only summary statistics available.
The limitations of our work are summarized as follows. First of all, our methods depend on a linear relationship between the covariates and the outcomes. This assumption is often satisfied in target application of genetic risk prediction [3], but it does limit the applicability of our method non-linear models, we may use the second-order Taylor approximation of the loss function as in [16]. However, the summary statistics used by such an algorithm are not found in existing literature or publicly available databases. Additionally, our theoretical results provide only upper bounds on the estimation error of the two estimators that we consider in this work. To fully characterize the cost of using summary statistics for multi-task learning, lower bounds resembling Theorem 2.2 of [15] are needed. We conjecture that our estimators converge at a minimax optimal rate, and we view the proof of this conjecture as an important future direction. Finally, we may also extend the framework of [10] to our summary statistic based setting adjust for potential differences between tasks.
other domains. To extend our framework
Nevertheless, our results have important implications beyond high-dimensional statistical theory. The trade-off between proxy data sample size and discovery-proxy overlap may inform how polygenic risk models are built in real-world applications: Practitioners should prioritise alignment between the sources of summary statistics that they use build these models, rather than optimizing for large sample sizes. This guidance may lead to more accurate polygenic scores, which have emerged as an important predictive tool in the field of precision medicine.
We recognize that the development of polygenic risk scores, if done without care, may worsen existing health disparities [23]. This is potential negative societal impact of our work. We hope that multi-task learning framework may be used to incorporate data from diverse populations to improve generalizability and transportability of genetic risk predictions to overcome these negative impacts.
Acknowledgements
Rui Duan is supported by National Institute of General Medical Sciences (NIGM) R01GM148494,
Parker Knight is supported by NSF Graduate Research Fellowship. We would like to thank Rajarshi Mukherjee for an insightful discussion of Lepski's method, and thank the reviewers for their comments and feedback which greatly improved the paper.
10

---

References
[1] Kenneth Burnham and David R. Anderson. "Multimodel Inference: Understanding AIC and BIC in Model Selection". In: Sociological Methods & Research 33.2 (Nov. 2004), pp. 261-304.
ISSN: 1049-1241, 1552-8294, DOI: 10.
[2] Han Cao et al. "dsMTL: a computational framework for srivacy-preeevvng distributed multi- task machine learning". In: Bioinformatiis 38.21 (Oct. 31, 2022). Ed. by Jonathan Wren, pp. 4919-4926. ISSN: 1367-4803, 1367-4811. DOI: 10. 1et1/099//bbo wooonnoormmticc.o bbtcc/6..........
[3] Nilanjan Chatterjee, Jianxin Shi, and Montserrat Gareía-Closas. "Developing and evaluating polygenic risk prediction models for stratified disease prevention". In: Nature Reviews Genetics 17.7 (July 2016), pp. 392-406. ISSN: 1471-0056, 1471-0064. DOE: 10. 1008/008//// 2016. 27.
[4] Ting-Huei Chen et al. "A Penalized Regression Framework for Building Polygenic Risk Models Based on Summary Statistics From Cenome-Wide Association Studies and Incorporating Ex- ternal Information". In: Journal of the American Statistical Association 116.33 (Jan. 2, 2021), pp. 133-143. ISSN: 0162-1459, 1537-274X. DOI: 10. 000/01622449 2020. 1764849.
[5] Michaël Chichignoud, Johannes Lederer, and Martin Wainwright. A Practical Scheme and Fast Algorithm to Tune the Lasso With Optimality Guarantees. arXiv:1410.0247 type: article. arXiv, Nov. 8, 2016. arviv: 1410. 0247 [math, stat].
[6] Shing Wan Choi, Timothy Shin-Heng Mak, and Paul F. O' O'Reilly. "Tutorial: a guide to perform- ing polygenic risk score analyses". In: Nature Protocols 15.9 (Sept. 1, 2020), pp. 2759-2772. ISSN: 1754-21 1154, 1750-2799. DOT: 10. 1038/ 241594-020----
[7] Simon S. Du et al. Few-Shot Learning via Learning the Representation, Provably. arXiv:2002.09434 type: article. arXiv, Mar. 30, 2021. arXiv: 2002. 09434 [cs, math, stat].
[8] Rui Duan, Yang Ning, and Yong Chen. "Heterogeneity- awaae and communication-efficcenn distributed statistical inference". In: Biometrika (Feb. 12, 2021). asab007. ISSN: 0006-3444, 1464-3510. DOT: 10. 1093/099/bbb amet/asab007
[9] Rui Duan et al. "Learning from electronic health records across multiple sites: A caanuunccation-fffccieeeee and American Medical Informatics Association 27.3 (Mar. 1, 2020), pp. 376-385
distributed algorithm". In: Journal of the
ISSN: 1527-
974X. DOE: 10. 1..11.........1 1093/j ami a/c2199.
[10] Yaqi Duan and Kaizheng Wang. Adaptive and Robust Multi-task Learning. arXiv:2202.05250.
type: article. arXiv, Apr. 1, 2022. arXiv: 2202. 05250 [cs, math, stat].
[11] Amadou Gaye et al. "DataSHIELD taking the analysis to the data, not the data to the analysis". In: International Journal of Epidemiology 43.6 (Dec. 2014), pp. 1929-1944. ISSN: 1464-3685. 0300-5771. DOT: 10. 1 1093/111100 jeee/eeee//// dyu188.
[12] Tian Gu, Yi Han, and Rui Duan. Robust angle-based transfer learning in high dimensions. Apr. 7, 2023, arXiv: 2210. 12759 [stat].
[13] Tian Gu, Phil H. Lee, and Rui Duan. "COMMUTE: Communication-effin transfer learning
for multi-site risk prediction"" In: Journal of Biomedical Informatics 137 (Jan. 2023), p. 104243. ISSN: 15320464. DOI: 10. 1016// jjb. 20i2.022... 104243.
[14] O. V. Lepski and V. G. Spokoiny. "Optimal pointwise adaptive methods in nonparametric estimation". In: The Annals of Statistics 25.6 (Dec. 1, 1997). ISSN: 0090-5364. DOI: 10. 1214/
aos/ aos/103074a 00000041088
[15] Sai Li, T. Tony Cai, and Hongzhe Li. "Estimation and Inference with Proxy Data and its Genetic Applications". In: arXiv:2201.c 0.2027777777 [math, stat] (Jan. 10, 2022). arXiv: 2201. 03727. [16] Sai Li, Tianxi Cai, and Rui Duan. Targeting Underreprsentedd Populations Precision Medicine: A Federated Transfer Learning Approach. Aug. 27, 2021. arXiv: 2108. 12112 stat].
[17] Kunpeng Liu et al. "Privacy-Preservine Multi-task Learning.. In: 2018 JEEE International Conference on Data Mining (ICDM). 2018 IEEE International Conference on Data Mining (ICDM). Singapore: IEEE, Nov. 2018, pp. 1128-1133, ISBN: 978-1-5386--15955 DOI: 10.
1109/ICDD 2018 00147.
[18] Karim Lounici et al. "Oracle inequalities and optimal inference under group sparsity". In:
The Annals of Statistics 39.4 (Aug. 2011). Publisher: Institute of Mathematical Statistics,
pp. 2164 2204. ISSN: 0090-5364, 2168-8966. DOI: 10. 1211//11-100899
[19] Karim Louniei et al. Taking Advantage of Sparsity in Multi-Task Learning. arXiv:0903. 1468.
type: article. arxiv, Mar. 8, 2009. arXiv: 0903. 1468 [math, stat].
11

---

[20] Chongliang Luo et al. "DLMM as lossless one-shot algorithm for collaborative multi-site distributed linear mixed models". In: Nature Communications 13.1 (Mar. 30, 2022), p. 1678.
ISSN: 2041-1723. Dor: 10. 1038/ 467-022-2910--4
[21] Chongliang Luo al. "ODACH: one-shot distributed algorithm for Cox model with heteroge- neous multi-center data". In: Scientific Reports 12.1 (Apr. 22, 2022), p. 6627. ISSN: 2045-2322
DOI: 10. 10111108/4440 s118///1198--1 0122-09069--0
[22] Timothy Shin Heng Mak et al. "Polygenic scores via penalized regression on summary statistics: MAK et al." In: Genetic Epidemiology 41.6 (Sept. 2017), pp. 469-480. ISSN:
07410395. DOI: 10. 1002/gepi. 22050.
[23] Alicia R. Martin et al. "Clinical use of current polygenic risk scores may exacerbate health disparities". In: Nature Genetics 51.4 (Apr. 2019), pp. 584-591. ISSN: 1061-4036, 1546-1718. DOL: 10. 10011008/4444 841588-
0379-X.
[24] Aaron J. Molstad et al. /Heterogeneity-awarr integrative analyses for ancestry-specifie associa- tion studies. June 8, 2023, arXiv: 2306. 05571 [stat].
[25] Sahand Negahban and Martin J. Wainwright. "Estimation of (near) low-rank matrices with noise and high-dimensional scaling". In: The Annals of Statisties 39.2 (Apr. 1, 201 1). ISSN: 0090-5364. DOI: 10. 1214/111//11 1022/1011200 10-A05050
[26] Sahand N. Negahban et al. "A Unified Framework for High-Dimensional Analysis of $MS- Estimators with Decomposable Regularizers". In: Statistical Science 27.4 (Nov. 1, 2012). ISSN: 0883-4237. DOT: 10. 11111/2111111000//2
[27] Angelika Rohde and Alexandre B. Tsybakov. "Estimation of high-dimensional low-rank matrices". In: The Annals of Statistics 39.2 (Apr. 1, 201 ISSN: 0090-5364. DOI: 10. 1214/ 10-A0S860.
[28] the eMERGE Team et al. "The eMERGE Network: A consortium of biorepositories linked to electronic medical records data for conducting genomic studies". In: BMC Medical Genomics 4.1 (Dec. 2011), p. 13. ISSN: 1755-8794. DOr: 10. (186/1755-8794-1
[29] Ye Tian, Yuqi Gu, and Yang Feng. Learning from Similar Linear Representations: Adaptivity, Minimaxity, and Robustness. arXiv:2303. type: article. arXiv, Mar. 30, 2023, arXiv:
17765 [cs, stat].
[30] Ali Torkamani, Nathan E. Wineinger, and Eric J. Topol. "The personal and clinical utility of polygenic risk scores". In: Nature Reviews Geneties 19.9 (Sept. 2018), pp. 581-590. ISSN:
1471-0056, 407144 DOI: 10. 10w11/033116-0000 018-
[31] Nilesh Tripuraneni, Chi Jin, and Michael I. Jordan. Provable Meta-Learning of Linear Repre- sentations. arXiv:2002.1 11684. type: article. arXiv, Dec. 31, 2021, arXiv: 2002. 11684 [ca,
stat].
[32] Joshua R Vest and Larry D Gamm. "Health information exchange: persistent challenges and new strategies: Table 1", In: Journal of the American Medical Informatics Association 17.3 (May 2010), pp. 288-294, ISSN: 1067-5027, 1527-974X. DOI: 10. 0006/113/550 /a136//jama.... a. 003673.
[33] Martin Wainwright. High-dimensional statistics: a non-asymptotic viewpoint. Cambridge series in statistical and probabilistic mathematics 48. Cambridge: New York, NY: Cambridge University Press, 2019. ISBN: 978-I-108-49800-99
[34] Zhiyu Yan et al. "A preserving and computation-efficiai federated algorithm for generalized linear mixed models to analyze correlated electronic health records data". In: PLOS ONE 18.1 (Jan. 17, 2023). Ed. by Sathishkumar V. F.,, e0280192. ISSN: 1932-6203. DOT: 10. 1371/journan pone. 0280192.
[35] Yu Zhang and Qiang Yang. "A Survey on Multi-Tass Learning". In: JEEE Transactions on Knowledge and Data Engineering 34.12 (Dec. 1, 2022), pp. 5586-56666
ISSN: 1041-4347,
1558-2191, 2326-3865. DOI: 10. 1109/TKDE. 2021. Ttt.22222..22.... 3070203.
[36] Yu Zhang and Qiang Yang. "An overview of multi-task learning". In: National Science Review
5.1 (Jan. 1, 2018), pp. 30 43. ISSN: 2095-5138, 2053-714X. DOI: 10. 1093/ nsr///x105..
12


Reviews:
[{"review_id": "hoTO70JPWV", "reviewer": "Reviewer_LHqp", "summary": "Multi-task learning is a powerful machine learning paradigm for integrating data from multiple sources to improve overall model performance. However, data-sharing constraints in healthcare settings hinder its application. To address this challenge, a flexible multi-task learning framework utilizing summary statistics from various sources is proposed, along with an adaptive parameter selection approach based on a variant of Lepski's method. A systematic non-asymptotic analysis characterizes the performance of the proposed methods under various regimes of sample complexity and overlap. Extensive simulations demonstrate the theoretical findings and the performance of the method, offering a flexible tool for training related models across various domains with practical implications in genetic risk prediction and other fields.", "questions": "1. How does the performance of this method compare to existing classical multi-task learning methods?\n2. Can this method be integrated with existing deep learning models, and if so, how?\n3. In what type of tasks is this method most suitable, and what type of tasks may require other methods and models?", "limitations": "N/A", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Multi-task learning is a promising approach to integrating data from multiple sources and improving individual task performance.\nData-sharing constraints in healthcare and biomedical research limit access to individual-level data, making summary statistics a valuable substitute. The proposed multi-task learning framework allows for simultaneous learning of multiple models using only publicly available summary statistics. A systematic non-asymptotic analysis characterizes the performance of the proposed methods under various regimes of sample complexity and overlap. An adaptive scheme for tuning parameter selection based on a variant of Lepski's method is proposed, allowing for data-driven tuning when only summary statistics are available. The framework has practical applications in genetic risk prediction and can be used to develop trans-ethnic prediction models. The ability to learn from summary statistics offers a versatile tool for developing models across various domains.\n\n\nIndeed, a major advantage of this method is that it has theoretical guarantees. Through systematic non-asymptotic analysis, this method provides theoretical guarantees for the performance of the multi-task learning framework based on summary statistics, especially under different regimes of sample complexity and overlap. These theoretical guarantees help us better understand the performance of this method and provide guidance for practical applications. Additionally, this method proposes an adaptive parameter selection approach based on a variant of Lepski's method, which can perform data-driven tuning of parameters when only summary statistics are available. The effectiveness of this method is also theoretically guaranteed, which can help us better select parameters and improve the performance of the model in practical applications. Therefore, this method has great practical value and provides a feasible approach for multi-task learning under data sharing constraints.\n\n\n", "weaknesses": "1. One major drawback of this method is the lack of comparison with classical multi-task learning methods. While this method provides a novel and useful approach for dealing with data sharing constraints, comparing its performance with existing multi-task learning methods would help to further evaluate its effectiveness. By comparing their performance, we can gain a better understanding of the strengths and weaknesses of these methods and provide better guidance for practical applications. Therefore, future research could consider comparing this method with existing multi-task learning methods.\n2. Another potential drawback of this method is that it may not be able to integrate with existing deep learning models. This method is based on training multiple models using basic summary statistics, while deep learning models typically require large amounts of raw data and complex network structures for training. Therefore, integrating this method with existing deep learning models may face many challenges, such as how to convert summary statistics into input data and how to design network structures suitable for summary statistics. While this method is useful in utilizing summary statistics, it may not be applicable in tasks that require processing large amounts of raw data and performing complex computations. Therefore, in practical applications, appropriate methods and models need to be selected based on task requirements.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "**Q1** Our methods build upon classical multi-task learning techniques, and enable fitting models only using basic summary statistics which are often made publicly available. The sparse $\\ell_{2,1}$ regularized estimator extends the group-sparse estimators studied in [4,6], while the nuclear norm estimator expands on the low-rank regression model described in [1]. The nuclear norm approach is closely related to the linear representation learning problem [2,5], constraining regression coefficients to a shared low-dimensional subspace.\n\nThe main motivations behind employing these multi-task learning methods but not others are\nthreefold: (1) In genetic risk prediction modeling, the additive nature of genetic effects makes\nlinear modeling the most effective choice. (2) Across populations and related phenotypes,\nsimilarities in genetic architectures can be characterized by distance measures in model\nparameters, leading to the enforcement of shared structures through penalty terms. (3) The\nwidespread availability of GWAS summary statistics, capturing the marginal correlation between each SNP-phenotype pair, further supports our approach.\n\nIn terms of other multi-task learning methods, our general approach of using proxy variables to compute the second-order structure of the predictor variables may be applied to Multi-Task Kernel Ridge Regression in Reproducing Kernel Hilbert Spaces, following [3]. In particular, we may use the proxy variables to compute the kernel matrix between the features. However, the limited availability of publicly accessible proxy kernel matrices in scientific applications hinders the potential usefulness of this method.\n\n\n\n**Q2** In real-world applications, fitting a deep learning model with only $X^TX$ and $X^TY$ available\nposes significant challenges, if not impossibilities. The summary stats represent highly\ncondensed summaries of linear relationships, whereas the effectiveness of deep learning\nmodels lie in their capacity to capture non-linear relationships. Hence, relying solely on such\ninformation is intuitively insufficient for successful deep learning model training.\n\nNevertheless, we hold the anticipation that training neural models could be feasible if new classes of summary statistics are created and incorporated. This approach would require practitioners to explore and develop innovative summary-level information. \n\nIn a more extreme scenario, the problem could be formulated as a federated learning problem [7], wherein data owners iteratively share gradient-type information for model updating. However, this approach differs from our motivation, which aims to leverage pre-existing summary stats instead of requiring continuous sharing during model training.\n\nAn intriguing avenue worth exploring is the integration of external basic summary statistics to assist in the training of neural models using internal datasets. This direction is currently under investigation, and we are actively exploring its potential implications.\n\n\n**Q3** Our method is particularly well-suited for applications where the dominant signals exhibit\nlinearity. For instance, in genetic risk prediction modeling, linear models have proven to be\neffective. This choice is driven by the constraint of utilizing only publicly available summary\nstatistics. We advocate that in many scenarios, linear models serve as robust and stable\nworking models. With our method, researchers without access to individual-level data sources\ncan still gain insight from summary statistics in terms of classification or risk modeling.\nMoving forward, we recognize the importance of developing novel methods and theories for\nother canonical statistical problems. These endeavors may necessitate summary statistics in\nnew forms, tailored to specific tasks, and represent crucial areas for future research.\n\n[1] S. Negahban and M. J. Wainwright, “Estimation of (near) low-rank matrices with noise and high-dimensional scaling,” Ann. Statist., vol. 39, no. 2, Apr. 2011, doi: 10.1214/10-AOS850.\n[2] S. S. Du, W. Hu, S. M. Kakade, J. D. Lee, and Q. Lei, “Few-Shot Learning via Learning the Representation, Provably.” arXiv, Mar. 30, 2021. Accessed: Feb. 10, 2023. [Online]. Available: http://arxiv.org/abs/2002.09434\n[3] Charles Micchelli and Massimiliano Pontil, “Kernels for Multi-task Learning,” Advances in Neural Information Processing Systems, vol. 17, 2004.\n[4] K. Lounici, M. Pontil, S. van de Geer, and A. B. Tsybakov, “Oracle inequalities and optimal inference under group sparsity,” The Annals of Statistics, vol. 39, no. 4, pp. 2164–2204, Aug. 2011, doi: 10.1214/11-AOS896.\n[5] N. Tripuraneni, C. Jin, and M. I. Jordan, “Provable Meta-Learning of Linear Representations.” arXiv, Dec. 31, 2021. Accessed: Sep. 08, 2022. [Online]. Available: http://arxiv.org/abs/2002.11684\n[6] K. Lounici, M. Pontil, A. B. Tsybakov, and S. van de Geer, “Taking Advantage of Sparsity in Multi-Task Learning.” arXiv, Mar. 08, 2009. Accessed: Sep. 08, 2022. [Online]. Available: http://arxiv.org/abs/0903.1468\n[7] M. I. Jordan, J. D. Lee, and Y. Yang, “Communication-Efficient Distributed Statistical Inference,” Journal of the American Statistical Association, vol. 114, no. 526, pp. 668–681, Apr. 2019, doi: 10.1080/01621459.2018.1429274.\n"}, {"review_id": "ODqqCFp1bm", "reviewer": "Reviewer_PbVj", "summary": "The paper addresses the problem of multi-task learning in settings where only summary statistics (instead of individual-level data) is available, which is a common scenario e.g. in medicine. The paper presents a framework for linear relationships between the covariates and the outcomes which uses summary-level information from distinct sources with a general data-driven tuning scheme for selecting tuning parameters. The results from the theoretical analysis of this framework are confirmed with numerical experiments on synthetically generated data. ", "questions": "- The framework is designed for multi-task learning with linear relationships between the covariates and the outcomes. Can you provide real-world examples that satisfy this constraint? \n- You note that ${\\bf X}^{(q)}$ and $\\tilde{\\bf X}^{(q)}$ are not necessarily the same. Can you give a practical example where they are the same, and one where they are not?\n- To obtain equation (2), you replace ${\\bf X}^T{\\bf X}$ with $\\tilde{\\bf \\Sigma}$, although you before mention that $\\tilde{\\bf \\Sigma}\\propto \\tilde{\\bf X}^T \\tilde{\\bf X}$, and ${\\bf X}$ and $\\tilde{\\bf X}$ are not the same. Why can you do this replacement here then?", "limitations": "authors adequately addressed the limitations and potential negative societal impact of their work", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The paper addresses an important problem for fields where individual level data is not publicly available, which is often the case in health care applications. The numerical experiments confirm the theoretical analysis. ", "weaknesses": "While the paper addresses an important topic, the presented framework comes with strong limitations. Particularly, the framework can only be applied for linear models which drastically reduces the potential applications of this approach. Additionally, the experimnts are limited to rather simple, synthetic examples that indeed show of the expected behaviour based on the theoretical results, but do not provide any insights on the applicability of the framework in real-world scenarios. Finally, the presentation could be improved (see comments below). \n\nComments on presentation:   \n- on related work: first, related work should also be a numbered section (not an unnumbered subsection of introduction). second, there was also quite recently some work by Meija et al. (2022) on estimating causal effects only based on summary statistics reported in medical studies, that seems quite related to the presented problem setup. There the authors show that merging datasets with maximum entropy improves the predictive power compared to using the observed marginal distributions as predictors.   \n- Line 99: public available -> publicly available. \n- After equation (1) you already start using variables with and without tilde ($\\tilde{\\bf X}$ and ${\\bf X}$, $\\tilde{\\mathcal{D}}$ and $\\mathcal{D}$), but you have not really introduced what the difference between them is and only introduce proxy and discovery data (without explaining the terms) in section 3. You should introduce the proxy data (and what that means) much earlier.   \n- Line 124: in a data-driven in the ... -> in a data-driven way in the ...  \n- Line 135 and equation above: In the equation use \\left( and \\right) for the brackets to scale them appropriately. Further, you use $\\varepsilon$ as random noise before (above line 88), and now define $\\Varepsilon$ as cost of using proxy data that is different from the discovery data. That's confusing. Similarly confusing is that you have now $\\tilde{\\Sigma}$, $\\Sigma_1$ and $\\Sigma_2$. Isn't $\\Sigma_2$ the same as $\\tilde{\\Sigma}$? (if not, please clarify what's the difference).  \n- Line 175: what is $\\tilde{p}$? wasn't there only $p$ and $\\tilde{\\rho}$ before?\n- position parameters of figure 2 and 3 should be set to [t]\n- Overall structure could be improved. Suggest to have a section on notation and assumptions used throughout the paper before the current section 2. Then you can introduce the problem setup. In that section you right now already mention how summary statistics get included. One could consider doing that only in the next section that focusses on your approach. This way it should be (hopefully) clearer what you propose and what is new about it. Right now it's difficult to understand what are common assumption, what is common knowledge, and what was added on top by you.  ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for excellent feedback about the presentation of our paper, and for raising several important questions. Here we address each of the comments on the paper’s presentation:\n* (Regarding the related work) Thank you for your feedback on formatting and for pointing out the recent work of Meija et al (2022). We will include this paper in our literature in the final version of our paper. \n* (Regarding the comment on Line 99) Thank you for pointing out this typo, this will be fixed.\n* (Regarding the comment on equation (1)) Thank you for this feedback, we will introduce a notation section in the final version of our paper to better clarify the difference between the discovery data X and the proxy data \\tilde{X}.\n* (Regarding the comment on Line 124) We will fix this typo.\n* (Regarding the comment on Line 135) We will change the notation from $\\mathcal{E}$ to a different character in the accepted version of our manuscript to prevent any confusion with the error term $\\varepsilon$. $\\Sigma_1$ is the population-level covariance of $X$, meaning that $\\Sigma_1 = E[\\frac1n X^TX]$. Similarly $\\Sigma_2 = E[\\frac1n\\tilde{X}^T\\tilde{X}]$ is the population-level covariance of the proxy data. We will clarify this in the notation section in our final version of this paper.\n* (Regarding the comment on Line 175) Thank you for pointing out this typo, this character should be $\\tilde{\\rho}$.\n\n**Q1** In our target application, which is genetic risk prediction, genetic variants contribute to observed traits additively, making linear modeling the most effective choice [3,4,5]. We will include a discussion of this assumption in the final version of our paper.\n\n**Q2** In biomedical applications, $X^TY$ can be obtained from studies that focus on marginal\nassociations between covariates and the outcome of interest. However, it is often the case that\nthese same studies may not report the correlations between the covariates. In this case, such\ninformation has to be derived from other studies or reference datasets. \nFor instance, in genetic studies, many research papers provide GWAS summary statistics\n(https://www.ebi.ac.uk/gwas/), which offer insights into the associations of genetic variants with\nthe outcome. However, the correlations between these genetic variants are often reported in\nonly a few studies, such as the UK Biobank [2] and the 1000 Genomes data [1].\n\n**Q3** The reviewer is correct that $X$ and $\\tilde{X}$ are not the same, and this is due to\nthe availability of summary statistics---studies that report $X^TY$ may not report $X^TX$,\nand we need to find a proxy for $X^TX$ from publicly available data. This replacement is how our estimator is defined.\nThe main contribution of our theory essentially is to quantify the\nfundamental cost of performing this replacement. If both $X^TY$ and $X^TX$ are\navailable from the same study, our theory suggests using the summary stats\nfrom the same study, instead of using a proxy dataset.\n\n\n[1] The 1000 Genomes Project Consortium et al., “A global reference for human genetic variation,” Nature, vol. 526, no. 7571, pp. 68–74, Oct. 2015, doi: 10.1038/nature15393.\n[2] C. Sudlow et al., “UK Biobank: An Open Access Resource for Identifying the Causes of a Wide Range of Complex Diseases of Middle and Old Age,” PLoS Med, vol. 12, no. 3, p. e1001779, Mar. 2015, doi: 10.1371/journal.pmed.1001779.\n[3] Chatterjee, Nilanjan, Jianxin Shi, and Montserrat García-Closas. \"Developing and evaluating polygenic risk prediction models for stratified disease prevention.\" Nature Reviews Genetics 17.7 (2016): 392-406.\n[4] Torkamani, Ali, Nathan E. Wineinger, and Eric J. Topol. \"The personal and clinical utility of polygenic risk scores.\" Nature Reviews Genetics 19.9 (2018): 581-590.\n[5] Choi, Shing Wan, Timothy Shin-Heng Mak, and Paul F. O’Reilly. \"Tutorial: a guide to performing polygenic risk score analyses.\" Nature protocols 15.9 (2020): 2759-2772.\n"}, {"review_id": "Boy7mKxaCg", "reviewer": "Reviewer_TqYc", "summary": "This paper presents an approach to learning predictive models from summary statistics in the setting of multi-task learning.  Linear model is assumed and least square solution was derived with either the \\ell_{2,1}-norm or kernel norm of the parameter matrix to capture relativeness among tasks. Theoretical results for bounding the estimator when proxy data are used to obtain the covariance matrix and hyperparameter tunning are provided. Synthetic data was used in their empirical evaluation.", "questions": "- No test set is mentioned in the description of their experiment. Is there a separate test set? I do not think training MSE is a good measure as it is affected by the sparse driving term in the optimization.  \n-  I found the results shown in Figure 2 is a bit confusing. Given n = \\tilde{n}, \\tilde{p} = 1 means the proxy data is exactly the true data, meaning the proxy covariance is exactly the true covariance. Also, covariance is the only input from the proxy data to obtain the solution according to the estimator. These together mean the solution obtained with proxy_cov should be the same as that with true_cov at \\tilde{p}=1. How the two plots show the former is the same as that with IL_cov?", "limitations": "NA", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "+ The studied problems, both training predictive models from summary statistics and the use of proxy data to estimate covariance are of great importance, given the constraints on data sharing in medical research and the lack of availability of covariance in current practice of GWAS. \n+ Theoretical results can be very useful in assessing the impact of the use of proxy data. \n+ The approach proposed to tune the hyperparameters is important given there is typically no validation set consisting of individual level data for such tunning.", "weaknesses": "- Even though the studied problems are significant, the proposed approach lacks technical innovations, seeming trivial extensions of existing works for both the bound derivation and hyperparameter tunning. The authors may point out the technique challenges in these extensions, highlighting their technical contribution. \n- Lack of evaluation on practical datasets. It is unclear how useful this method is in practical problems.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the valuable feedback, and for recognizing the strengths of our paper and the importance of the problem that we address. \n\nRegarding the technical innovation of our paper: The proofs of the theoretical results present unique challenges in the proxy data setting. The derivation of the error bounds relies on a careful analysis of the contributed variance from each data point $X_i$ and $\\tilde{X}_i$, which to our knowledge is entirely novel for high dimensional linear models. The details of this analysis can be found in the proof of Lemma A.4 in our supplement. The closest related result to ours is Theorem 2.1 of [2], but our results generalize this theorem by accounting for the overlap between the reference data and the discovery data and by allowing for multiple outcomes. Our derivations for the hyperparameter tuning procedure are generalizations of the techniques used in [1]. To our knowledge, this is the first application of Lepski’s method to high-dimensional regression problems beyond the Lasso and is the first time that the use of Lepski’s method has been motivated by data-sharing restrictions. Our proofs for the model tuning procedure rely on the properties of general convex regularizers, analogous to the framework provided in [3]. \n\nRegarding evaluation on practical datasets: Thank you for this suggestion. We have performed a data analysis demonstrating the performance of our method on a real genetic dataset, where the goal is to predict low-density lipoprotein (LDL) levels using genotype data. The results and an in-depth description of the analysis are given in our global response.\n\n**Q1** Our previous simulation results present the mean-squared error of parameter estimation per task; in other words, we present the quantity $\\frac{1}{nQ}\\|\\|\\hat{B} - B^*\\|\\|^2$. We use the entire dataset to compute $\\hat{B}$, and so we did not use a test set. We will clarify this in the final version of our manuscript. To address your concerns, our new real data application uses a standard train-test split to evaluate the prediction performance. We provide more details in the global response. Additionally, we have performed an additional simulation study that analyzes the impact of the proxy data sample size on prediction performance for each of our estimators, under the regime of no-overlap between the proxy and discovery data. The details and results of this simulation are in the additional PDF provided in the global response.\n\n**Q2** When $\\rho = 1$, the proxy data is the same as the individual-level data, meaning that $\\hat{\\Sigma} = \\tilde{\\Sigma}$. In this case, the Proxy_Cov estimator is expected to perform as well as the IL_Cov estimator, which is confirmed by our simulation results in Figure 2. The ‘true’ covariance label refers to the underlying population-level covariance which we do not observe in practice. The True_Cov estimator uses this population-level matrix as the input to our estimator. This is not possible in practice, but we include the results on our figure to demonstrate that overlap between the proxy data and discovery data is more important for good statistical performance than having an increasingly large reference panel.\n\n\n[1] M. Chichignoud, J. Lederer, and M. Wainwright, “A Practical Scheme and Fast Algorithm to Tune the Lasso With Optimality Guarantees.” arXiv, Nov. 08, 2016. Accessed: Apr. 24, 2023. [Online]. Available: http://arxiv.org/abs/1410.0247\n[2] S. Li, T. T. Cai, and H. Li, “Estimation and Inference with Proxy Data and its Genetic Applications,” arXiv:2201.03727 [math, stat], Jan. 2022, Accessed: Mar. 07, 2022. [Online]. Available: http://arxiv.org/abs/2201.03727\n[3] S. N. Negahban, P. Ravikumar, M. J. Wainwright, and B. Yu, “A Unified Framework for High-Dimensional Analysis of $M$-Estimators with Decomposable Regularizers,” Statist. Sci., vol. 27, no. 4, Nov. 2012, doi: 10.1214/12-STS400.\n"}, {"review_id": "3YGQqgLveK", "reviewer": "Reviewer_Z4Yp", "summary": "The paper proposed a multi-task method to learn individual models without having access to the raw data but using summary statistics data for each task. The method is only applicable to linear models. ", "questions": "1. Why is the discovery data estimated from the actual observed data, unlike the sample covariance matrix was estimated from a proxy data?\n2. The proposed method works only for linear models. Can the same analysis be applied using non-linear (deep learning) models?\n3. In Figure 1, can you explain the behavior that the proxy covariance gave lower MSE than the tru covariance?", "limitations": "Yes. ", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1.  Theoretical guarantee for the optimal estimator \n2.  Good experiments on simulated data", "weaknesses": "1. Missing experiments on real data\n2. Applicable only to linear models. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Regarding experiments on real data: We have included additional results demonstrating our method on a prediction task with real genetic data. We have provided a description of the analysis and the results in our global response.\n\n**Q1** In statistical genetics applications, it is often the case that $X^TY$ is reported from studies that analyze the marginal associations between the covariates, which often correspond to genetic markers, and the outcome of interest. However, these studies rarely report the covariance structure between genetic markers, and this information needs to be obtained from additional external datasets, such as the 1000 Genomes data [1]. Our formulation captures this discrepancy between the data used to compute the marginal associations and the covariance matrix.\n\n**Q2** Our results are tailored towards linear models, but our approach may be applied to any estimator whose loss function depends only on the sample covariance matrix and summary association statistics. We anticipate that similar results may be possible for other classes of (non-linear) models such as deep neural networks. However, to fit a neural network using only summary-level information, practitioners need to be able to share new classes of summary statistics that correspond to more complex loss functions. The problem could be formulated as a federated learning problem [2], in which data owners iteratively share gradient-type information for model updating. However, this approach differs from our motivation, which aims to leverage pre-existing summary stats instead of requiring continuous sharing during model training.\n\n**Q3** The behavior in Figure 1 is due to a pre-processing step in the simulation that regularizes the proxy data for better numerical stability across replications. If this pre-processing step is not performed, the behavior of the proxy data estimator converges to that of the true covariance, as predicted by our theory. In the new simulation result described in our global response, we do not perform this pre-processing step, and the estimator behaves as expected according to our theory.\n\n[1] The 1000 Genomes Project Consortium et al., “A global reference for human genetic variation,” Nature, vol. 526, no. 7571, pp. 68–74, Oct. 2015, doi: 10.1038/nature15393.\n[2] M. I. Jordan, J. D. Lee, and Y. Yang, “Communication-Efficient Distributed Statistical Inference,” Journal of the American Statistical Association, vol. 114, no. 526, pp. 668–681, Apr. 2019, doi: 10.1080/01621459.2018.1429274.\n"}, {"review_id": "vhkny2YThk", "reviewer": "Reviewer_wP1d", "summary": "The work considers multi-task learning in settings where for each task only summary statistics X.T@Y and ẍ.T@ẍ are made available. The setting is motivated by healthcare and biomedical related research, where sharing individual level microdata is restricted by regulations due to privacy concerns. The authors assume a linear and sparse or low-rank underlying model. A regularized least-squares type of optimization framework is used, and it is highlighted that only summary statistics are needed to solve it. A main contribution of the paper is theoretical analysis bounding the the quality of the sparse or low-rank estimator taking into account quantities such as overlap and distributional shift between the so called discovery and proxy data used to calculate the summary statistics. Additionally, the problem of hyperparameter tuning is considered as basic (cross-)validation methods are not applicable without access to individual data, and Lepski’s method is proposed as a practical alternative that requires only summary level data as well as some prior knowledge in form of a constant C parameter. Experiments on simulated data show results that are consistent with the theoretical analysis.", "questions": "Can you clarify on the assumption that the discovery and proxy data sets can be separate, when would this in practice be an important consideration?\nFigure 3: am I reading this correctly that the holdout method does not benefit at all from having larger validation sample size? Any explanation for this?\n", "limitations": "-", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "- In general, approaches that allow learning from summary statistics are becoming ever more relevant to the community as regulation on sharing individual level data is tightening \n- The exact formulation of the multi-task summary statistic learning problem considered in the paper, where separate discovery and proxy data sets are assumed appears to be novel, as consequently the theoretical bounds that take this into account provide also novel insights about the effects this can have on learning. \n- Both L1,2 and nuclear norm regularization based variants of the learning problem are analyzed\n- The proposed approach for model selection provides a practical tool for hyperparameter tuning when only summary statistics are available, and might have implications also for other variations of the problem setting?\n- Fairly clear writing and technically rigorous", "weaknesses": "- Real-world relevance: while applications such as polygenic risk prediction are mentioned, there are no convincing examples, real applications of benchmark data sets provided where the assumptions of the learning setting (having related tasks, only X(q).T@Y(q) and ẍ(q).T@ẍ(q) available, perhaps X(q) != ẍ(q)) would hold. Thus the question about practical significance of the exact considered problem setting and proposed solution remain unclear.\n- Related to above point - no experiments on real-world data that would demonstrate the benefits of the approach\n- Somewhat restrictive assumptions such as assumptions about the linearity of the underlying model and relationships between the tasks ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Regarding the real-world application: To address this concern, we have included a real data application to polygenic risk prediction. Details of our analysis and the results are provided in the global response.\n\n**Q1** In statistical genetics applications, many studies that investigate the marginal relationships between covariates and the outcome will report $X^TY$ but do not report $X^TX$. For this reason, the sample covariance must be computed from an external reference data set, such as the 1000 Genomes data [1]. The formulation of our problem and the theory that we develop captures this discrepancy between the datasets. \n\n**Q2** The holdout method does not benefit from a larger sample size because the estimation accuracy will depend only on the sample size of the training dataset, as long as the tuning parameter is chosen to be the right order of magnitude. For the results in Figure 3, the grid of tuning parameters over which the holdout validation was performed was specified to be roughly the appropriate order. \n\n[1] The 1000 Genomes Project Consortium et al., “A global reference for human genetic variation,” Nature, vol. 526, no. 7571, pp. 68–74, Oct. 2015, doi: 10.1038/nature15393.\n\n"}]

Author rebuttal:
We would like to thank the reviewers for their valuable and insightful feedback. In this global response we would like to describe the major changes we made to address the reviewers’ comments, which greatly improved the quality of our work. 

**First, we add a new simulation study evaluating our method based on out-sample prediction accuracy.** The setup of the simulation is the same as what has been considered in the previous submission. We simulated a dataset with $n_{min} = 100$, $p = 100$,  $\tilde{n} = \tau n_{min}$ for $\tau \in (0.5, 1, 5, 10)$ and $\tilde{\rho}_q= 0$ for each q. The number of tasks was fixed at 8. We generated a row-sparse matrix with 10 nonzero rows and a rank 2 $B^*$ for the sparse and low-rank multi-task estimators. After fitting our model and obtained the parameter estimate for $B^*$ , we generate a new test set of size 100 according to the same data generating process and evaluate our estimator's ability to predict the outcome in this test set. The results are given in the first figure in our new additional PDF document. We found that the relative performance of our method compared to the individual level estimator and the true covariance estimator is similar to the estimation task we showed in our previous submission. The performance of the proxy data estimator improves as the size of the proxy data grows, matching that of the oracle estimator that uses the true covariance matrix as its input. We observe the same performance gap between the proxy data estimator and the individual-level data estimator.

**Secondly, we have applied our method to analyze real genetic data to demonstrate the real-world applicability of our method.** We use a multi-site data obtained from the electronic Medical Records and Genomics (eMERGE) network [1], which includes individual-level genotype data from multiple research sites in the United States. Our goal is to predict levels of low-density lipoprotein (LDL) across five adult sites, treating the data from each site as a separate task. We split the data (with sample sizes $n_1 = 3813, n_2 =  546, n_3 = 2666, n_4 = 1435, n_5 = 525$) at each task into a training and test set (with a test set data size of 100 for each task) and evaluate the performance of our method using the prediction MSE on the test set. 
The training data from each site is used to construct the discovery summary statistics $X_q^TY_q$ for each task. For approximating $\Sigma_q$, we choose two different approaches: one is to use the half of the genotype data from each site (this approach is labeled as Proxy_MTL1); the other approach is to use $X_1$ (genotype data from site 1) to approximate $\Sigma_q$ for all the sites. This approach is labeled as Proxy_MTL2. We use these two approaches to demonstrate a potential trade-off in the construction of the reference panel: Proxy_MTL1 uses a well-specified reference dataset with a smaller sample size; Proxy_MTL2 uses a larger reference dataset that may suffer from a distribution shift. For comparison, we also fit a multi-task learning estimator that uses all of the individual level training data for each task, which is labeled ‘Individual_MTL’, and we fit a ridge regression estimator that models each task separately, for comparison to our multi-task learning approach. The ridge estimator uses the proxy sample covariance instead of the individual level covariance matrix for a fair comparison with our method. We repeat the train-test split process 10 times, which admits the distribution of prediction MSE values that we report in the figure. We use the nuclear norm penalized multi-task learning estimators in this application, because we believe that the genetic effects in this dataset are dense. 

Our results demonstrate that our method is highly practical when only summary-level information is available, as the prediction MSE of our method is nearly the same as the estimator which uses the individual-level data, despite a slight cost in performance. Furthermore, all multi-task learning estimators outperform the ridge-estimator, confirming that multi-task learning is a strong approach when there is shared structure between tasks. 

**Finally, we clarified the motivation and provided real-world justifications for considering the discovery data to be different from the proxy data, and we clarified our technical innovation.** Due to the space limit, we refer it to the point-to-point responses to each reviewer.

Once again, we thank the reviewers for their comments, and we look forward to further discussion.

[1] McCarty, C.A., Chisholm, R.L., Chute, C.G. et al. The eMERGE Network: A consortium of biorepositories linked to electronic medical records data for conducting genomic studies. BMC Med Genomics 4, 13 (2011). https://doi.org/10.1186/1755-8794-4-13


Meta-review: This paper addresses the problem of multi-task learning where only summary statistics available. This setting arises in medical or biological applications such as genetic risk prediction. The approach is limited to a linear setting. A theoretical analysis and an evaluation are provided.  


The reviewers have acknowledged that the paper studies an important an interesting problem and the theoretical results are a strong point of the paper. Some reservations were also raised on the limited experimental evaluation, on the limitation to linear models, and on the paper presentation. 

 
During rebuttal, authors have addressed the points raised by the reviewers, including:  
(i) a new experimental study (simulation study on out-sample accuracy, analysis of real genetic data).  
(ii) some modifications to improve the presentation of the paper.  
(iii) the addition of a little discussion about the lack of comparison with other multi-task learning methods, the interest of the method for applications where dominant signal exhibits linearly (reviewer LHqp). 
(iv) a justification of the linear modeling as the most effective choice for the target task genetic risk prediction, as well as with a comment on  the potential use of GWAS summary statistics in biomedical applications (reviewer PbVj).
(v) an answer on the difference between $X$ and $\var{X}$ (reviewers PbVj and Z4Yp).
(vi) a discussion on the perspective of integration with deep learning models (reviewers LHqp and Z4Yp).
(vii) an explanation to the behavior in Figure 1 commented by reviewer Z4Yp. 

Overall, authors have answered mostly to the issues raised by reviewers. The global evaluation is positive.  
I propose then acceptance.  
I encourage the authors to add all the elements provided in the answers of the rebuttal in the final version of the paper
Decision: Accept (poster)