For high-dimensional hierarchical models, consider exchangeability of effects across covariates instead of across datasets
Year: 2022

Abstract:
Hierarchical Bayesian methods enable information sharing across regression problems on multiple groups of data. While standard practice is to model regression parameters (effects) as (1) exchangeable across the groups and (2) correlated to differing degrees across covariates, we show that this approach exhibits poor statistical performance when the number of covariates exceeds the number of groups. For instance, in statistical genetics, we might regress dozens of traits (defining groups) for thousands of individuals (responses) on up to millions of genetic variants (covariates). When an analyst has more covariates than groups, we argue that it is often preferable to instead model effects as (1) exchangeable across covariates and (2) correlated to differing degrees across groups. To this end, we propose a hierarchical model expressing our alternative perspective. We devise an empirical Bayes estimator for learning the degree of correlation between groups. We develop theory that demonstrates that our method outperforms the classic approach when the number of covariates dominates the number of groups, and corroborate this result empirically on several high-dimensional multiple regression and classification problems.

text: Hierarchical Bayesian methods enable information sharing across regression problems on multiple groups of data. While standard practice is to model regression parameters (effects) as (1) exchangeable across the groups and (2) correlated to differing degrees across covariates, we show that this approach exhibits poor statistical performance when the number of covariates exceeds the number of groups. For instance, in statistical genetics, we might regress dozens of traits (defining groups) for thousands of individuals (responses) on up to millions of genetic variants (covariates). When an analyst has more covariates than groups, we argue that it is often preferable to instead model effects as (1) exchangeable across covariates and (2) correlated to differing degrees across groups. To this end, we propose a hierarchical model expressing our alternative perspective. We devise an empirical Bayes estimator for learning the degree of correlation between groups. We develop theory that demonstrates that our method outperforms the classic approach when the number of covariates dominates the number of groups, and corroborate this result empirically on several high-dimensional multiple regression and classification problems.

heading: 1 Introduction
text: Hierarchical modeling is a mainstay of Bayesian inference. For instance, in (generalized) linear models, the unknown parameters are effects, each of which describes the association of a particular covariate with a response of interest. Often covariates are shared across multiple related groups, but the effects are typically allowed to vary both by group and by covariate. A classic methodology, dating back to Lindley and Smith (1972) [44], models the effects as conditionally independent across groups, with a latent (and learnable) degree of relatedness across covariates. From a practical standpoint, the model is motivated by the understanding that it “borrows strength” across different groups [24, Chapter 5.6]. Mathematically, the model is motivated by assuming effects are exchangeable across groups and applying a de Finetti theorem [44, 35]. The methodology of Lindley and Smith is ubiquitous when the number of groups is larger than the number of covariates. It is a standard component of Bayesian pedagogy [[23, Chapter 13.3]; [24, Chapter 15.4]] and software; e.g. it is used in the mixed modeling package lme4 [5], which has over 16 million downloads at the time of writing.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Despite its resounding success when there are more groups than covariates, we show in the present work that this standard methodology performs poorly when there are more covariates than groups. To address the many-covariates case, we turn for inspiration to statistical genetics, where scientists commonly learn linear models relating genetic variants (covariates) to traits (corresponding to different groups) across individuals (which each exhibit a response). These applications may exhibit millions of covariates, thousands of responses, and just a handful of groups. In these cases, [39, 12, 56, 69, 46, 53] use a multivariate Gaussian prior akin to that of Lindley and Smith, but assume conditional independence across covariates and prior parameters that encode correlations across groups, rather than than the other way around.
As we will see, this alternative modeling approach may be motivated from a Bayesian perspective when one begins from an assumption of a priori exchangeability of the effects across covariates (rather than across groups). This exchangeability assumption is reasonable in statistical genetics, where we have little knowledge to distinguish our expectations about the effects of different genetic variants; we argue this modeling approach can be effective other modern high-dimensional analyses of multiple groups of data (beyond statistical genetics) in which large collections of covariates are frequently treated monolithically, e.g. by applying ridge regression. Namely, when there are more covariates than groups, we propose to model the effects as exchangeable across covariates (rather than groups) and learn the degree of relatedness of effects across groups (rather than covariates). In what follows, we refer to this framework as ECov, for exchangeable effects across covariates, and distinguish it from exchangeable effects across groups or EGroup.
While the existing methods in statistical genetics for modeling multiple traits obtain as a special case of ECov, to the best of our knowledge this approach is absent from existing literature on hierarchical Bayesian regression. Brown and Zidek (1980) [10] and Haitovsky (1987) [28] form two exceptions, but these two papers (1) consider only the situation in which a single covariate matrix is shared across all groups (or equivalently, for each data point all responses are observed) and (2) include only theory and no empirics. While Lindley and Smith (and others) discuss a priori exchangeability across covariates in the context of analysis of a single group, to our knowledge no other work has pushed this idea forward to share strength across multiple groups.
We suspect that the historical origins of the methodology in statistical genetics may have hindered earlier expansion of this class of models to a wider audience. In particular, this literature traces back to mixed effects modeling for cattle breeding [57]; here, an even-earlier notion of the genetic contribution of trait correlation (i.e. “genetic correlation;” see Hazel (1943) [29]) informs the covariance structure of random effects. Although genetic correlation is now commonly understood to describe the correlation of effects of DNA sequence changes on different traits [12], its provenance predates even the first identification of DNA as the genetic material in 1944 [3]. As such, this older motivation obviated the need for a more general justification grounded in exchangeability. See Appendix A for further discussion of related work, including more recent works from within the machine learning community on sharing strength across multiple groups of data.
In the present work, we propose ECov as a general framework for hierarchical regression when the number of covariates exceeds the number of groups. We show that the classic model structure from statistical genetics can be seen as an instance of this framework, much as Lindley and Smith give a (complementary) instance of an EGroup framework. To make the ECov approach generally practical, we devise an accurate and efficient algorithm for learning the matrix of correlations between groups. We demonstrate with theory and empirics that ECov is preferred when the number of covariates exceeds the number of groups, while EGroup is preferred when the number of groups exceeds the number of covariates. Our experiments analyze three real, non-genetic groups in regression and classification, including an application to transfer learning with pre-trained neural network embeddings. We provide proofs of theoretical results in the appendix.

heading: 2 Exchangeability and its applications to hierarchical linear modeling
text: We start by establishing the data and model, motivating exchangeability among covariate effects (ECov), and motivating our Bayesian generative model.
Setup and notation. Consider Q groups with D covariates. Let Nq be the number of data points in group q. For the qth group, the Nq ⇥D real design matrix Xq collects the covariates, and Y q is the Nq-vector of responses. The nth datapoint in group q consists of covariate D-vector Xqn and
scalar response Y qn . We let D := {(Xq, Y q)} Q q=1 denote the collection of data from all Q groups. We consider the generalized linear model Y qn |Xqn, q indep⇠ p(·|Xq>n q) with unknown D-vector of real effects q . We collect all effects in a D⇥Q matrix with (d, q) entry qd . The linear form of the likelihood allows interpretation of qd as the association between the dth covariate and the response in group q. In linear regression, the responses are real-valued and the conditional distribution is Gaussian. In logistic regression, the responses are binary, and we use the logit link. The independence assumption conflicts with some models that one might use, for example in some cases when the different groups partially overlap.
Example. As a motivating non-genetics example, consider a study of the efficacy of microcredit. There are seven famous randomized controlled trials of microcredit, each in a different country [48]. We might be interested in the association between various aspects of small businesses (covariates), including whether or not they received microcredit, and their business profit (response). In this case, the dth element of Xqn would be the dth characteristic of the nth small business in the qth country, and Y qn is the profit of this business. See the experiments for additional examples in rates of policing, web analytics, and transfer learning.
Exchangeable effects across groups (EGroup). To fully specify a Bayesian model, we need to choose a prior over the parameters . Lindley and Smith assume the effects are exchangeable across groups. Namely, for every Q-permutation , p( 1, 2, . . . , Q) = p( (1), (2), . . . , (Q)). Assuming exchangeability holds for an imagined growing Q and applying de Finetti’s theorem motivates a conditionally independent prior. Concretely, Lindley and Smith take q i.i.d.⇠ N (⇠, ), for D-vector ⇠ and D⇥D covariance matrix . The (d, d0) entry of captures the degree of relatedness between the effects for covariates d and d0. Both ⇠ and may be learned in an empirical Bayes procedure. However, when D is large relative to Q, learning these parameters can present both computational and inferential challenges, as the O(D2) degrees of freedom in outnumber the O(DQ) effects.
Exchangeable effects across covariates (ECov). We here argue for a complementary approach in settings where D > Q. In the microcredit example, notice that D > Q will arise whenever the experimenter records more characteristics of a small business than there are locations with microcredit experiments; that is, D > 7 in this particular case. Concretely, let d be the Q-vector of effects for covariate d across groups. Then, in the ECov approach, we will assume that effects are exchangeable across covariates instead of across groups. Namely, for every D-permutation , p( 1, 2, . . . , D) = p( (1), (2), . . . , (D)). We will see theoretical and empirical benefits to ECov in later sections, but note that the ECov assumption is often consistent with prior beliefs in high dimensional settings. For instance, regarding microcredit, we may have no prior knowledge about how effects differ for distinct small-business characteristics. And we may a priori believe that different countries could exhibit more similar effects – and wish to learn the degree of relatedness across those countries.
We may apply a similar rationale as Lindley and Smith to motivate a conditionally independent model. Analogous to Lindley and Smith, we propose a Gaussian prior: d
i.i.d.⇠ N (0,⌃). ⌃ is now a Q⇥Q covariance matrix whose (q, q0) entry captures the similarity between the effects in the q and q0 groups. For simplicity, we restrict to E[ d] = 0; see Appendix E.3 for discussion. Another potential benefit to ECov relative to EGroup is that we might expect a statistically easier problem, with O(Q2) rather than O(D2) values to learn in the relatedness matrix. We provide a rigorous theoretical analysis in Sections 4 and 5.

heading: 3 Our method
text: We next describe our inference method for specific instances of the exchangeable covariate effects model of Section 2. We compute the posterior and take an empirical Bayes approach to estimate ⌃. We find that an expectation maximization (EM) algorithm estimates ⌃ effectively; Appendix A.2 compares our approach to existing methods for the related problem of estimating for EGroup.
Notation. We identify estimates of and ⌃ with hats. For instance, ̂LS is the least squares estimate, with ̂qLS := (X
q>Xq) 1Xq>Y q. We will sometimes find it useful to stack the columns of or its estimates into a length DQ vector; we denote such vectors with an arrow; for example,
Algorithm 1 Expectation Maximization for Exchangeability Among Covariate Effects
1: // Initialize covariance 2: ⌃(0) IQ 3: // Run EM algorithm 4: for t = 0, 1, . . . do 5: // Expectation step 6: µ1, . . ., µD, V1, . . ., VD E_Step(⌃(t)) 7: 8: // Maximization step 9: ⌃(t+1) D 1 PD d=1(µdµ > d + Vd)
10: 11: Return ⌃(t+1)
Algorithm 2 E-Step: Linear Regression
1: ~µ, V E[~ |D,⌃],Var[~ |D,⌃] 2: for d = 1, . . ., D do 3: µd (ed ⌦ IQ)>~µ 4: Vd (ed ⌦ IQ)>V (ed ⌦ IQ) 5: Return µ1, . . ., µD, V1, . . ., VD
Algorithm 3 E-Step: Logistic Regression
1: ~µ⇤ argmax~ log p(~ |D,⌃) 2: V [r2 log p(~ |D,⌃)
~ =~µ⇤
] 1
3: for d = 1, . . ., D do 4: µd (ed ⌦ IQ)>~µ⇤ 5: Vd (ed ⌦ IQ)>V (ed ⌦ IQ) 6: Return µ1, . . ., µD, V1, . . ., VD
~ := [ 1>, 2>, . . . , Q>]>. For a natural number N, we use IN ,1N , and eN to denote the N ⇥N identity matrix, N -vector of ones, and N th basis vector, respectively. We use ⌦ to denote the Kronecker product.

heading: 3.1 Posterior inference with a Gaussian likelihood
text: We first consider a Gaussian likelihood: for each group q and observation n, we take Y qn |Xqn, q indep⇠ N (Xq>n q, 2q ) where 2q is a group-specific variance. When the relatedness matrix ⌃ is known, a natural estimate of is its posterior mean. We obtain the full posterior, including its mean, via a standard conjugacy argument; see Appendix B.1:
Proposition 3.1. For each covariate d, let d i.i.d.⇠ N (0,⌃) a priori. For each group q and data point n, let Y qn |Xqn, q indep⇠ N (Xq>n q, 2q ). Then ~ |D,⌃ ⇠ N (~µ, V ) for ~µ = V [ 21 Y 1>X1, . . . , 2Q Y Q>XQ]> and V 1 = ⌃ 1⌦ID+diag( 21 X1>X1, . . . , 2 Q X Q>XQ), where diag( 21 X 1>X1, . . . , 2Q X Q>XQ) denotes a DQ⇥DQ block-diagonal matrix.
At first glance, the posterior mean ~µ for this model might seem to introduce a computational challenge because exact computation of V involves an O(D3Q3)-time matrix inversion. Our experiments (Section 6), however, involve on the order of DQ ⇡ 1,000 parameters, so direct inversion of V demands less than a single second. Moreover, in much larger problems ~µ may still be computed very efficiently using the conjugate gradient algorithm [49, Chapter 5], with convergence in a small number of O(D2Q) time iterations; see Appendix B.2.

heading: 3.2 Empirical Bayes estimation of ⌃ by expectation maximization
text: The posterior mean of in Proposition 3.1 requires ⌃, which is typically unknown. Accordingly, we propose an empirical Bayes approach of estimating ⌃ by maximum marginal likelihood:
̂ECov := E[ | D, ⌃̂] where ⌃̂ := argmax ⌃⌫0 p(D | ⌃). (1)
Equation (1) defines a two step procedure. In the first step, we learn the similarity between groups via estimation of ⌃. In the second step, we use this similarity to compute an estimate, ̂ECov, that correspondingly shares strength. Though we have been unable to identify a general analytic form for ⌃̂, we can compute it with an expectation maximization (EM) algorithm [47, Chapter 1.5]. Algorithm 1 summarizes this procedure; see Appendix B.3 for details.

heading: 3.3 Classification with logistic regression
text: We can extend the approach above to inference for multiple related classification problems. We assume a logistic likelihood; for each q and n, Y qn |Xqn, q
indep⇠ Bern[(1 + exp{ Xq>n q}) 1]. In the classification case, we cannot use Gaussian conjugacy directly, so we apply an approximation. Specifically, we adapt the original E-step in Algorithm 3 by using a Laplace approximation to the posterior [7, Chapter 4.4]. We approximate the posterior mean of by the maximum a posteriori value. We leave extensions to other generalized linear models to future work.

heading: 4 Theoretical comparison of frequentist risk
text: In this section, we prove theory that suggests ECov has better frequentist risk than EGroup when D is large relative to Q. Analyzing ̂ECov directly is challenging due to its non-differentiability as a function of the data, so we take a multipart approach. First, in Theorem 4.3, we show that an ECov estimate based on moment-matching (MM), ̂MMECov, dominates least squares, ̂LS, when D is large relative to Q; ̂LS in turn dominates ̂MMEGroup (a similar estimator for EGroup). Second, in Theorem 4.5, we show that ̂ECov uniformly improves on ̂MMECov.
Setup. Take a fixed value of and an estimator ̂. We use squared error risk, R( , ̂) := E h k̂ k2F | i , as our measure of performance. k · kF is the Frobenius norm of a matrix,
and the expectation is over all observations Y 1, . . . , Y Q jointly. We require the following orthogonal design condition. Condition 4.1. For each group q, 2q Xq>Xq = 2ID for some shared variance 2.
Though restrictive, this condition is useful for theory, as other authors have found; see Appendix C.1. We empirically demonstrate that our theoretical conclusions apply more broadly in Section 6.
ECov vs. EGroup when using moment matching in high dimensions. For ECov, the following estimate for ⌃ is unbiased under correct prior specification: ⌃̂MM := D 1̂>LŜLS D 1diag( 21kX1†k2F , . . . , 2QkXQ†k2F ), where † denotes the Moore-Penrose pseudoinverse of a matrix and ̂LS is the least squares estimate. We define ̂MMECov := E[ |D, ⌃̂MM] to be the resulting parameter estimate, and define ̂MMEGroup analogously for EGroup; see Appendix C.2 for details. While ̂MMECov and ̂ MM EGroup are naturally defined only when D Q and D  Q, respectively, we find it informative to compare how their performances depend on D and Q nonetheless.
Before our theorem, a lemma provides concise expressions for the risks of ̂MMECov and ̂ MM EGroup.
Lemma 4.2. Under Condition 4.1 and when D Q, R( , ̂MMECov) = 2DQ 4D(D 2 2Q)E[k̂†LSk2F | ]. Additionally, when D  Q, R( , ̂MMEGroup) = 2DQ 4Q(Q 2 2D)E[k̂†LSk2F | ].
Lemma 4.2 reveals forms for the risks of ̂MMECov and ̂ MM EGroup that are surprisingly simple. The symmetry between the forms and risks of these estimators, however, is intuitive; under Condition 4.1, ̂MMECov and ̂ MM EGroup can be seen as respectively arising from the same procedure applied to ̂LS and its transpose.
With Lemma 4.2 in hand, we can now compare the risk of ̂MMECov, ̂LS, and ̂ MM EGroup.
Theorem 4.3. Let Condition 4.1 hold. Then (1) if D > 2Q+ 2, ̂MMECov dominates ̂LS with respect to squared error risk. In particular, for any , R( , ̂MMECov) < R( , ̂LS). Additionally, (2) if D > Q/2 1, ̂MMEGroup is dominated by ̂LS.
Since ̂LS is minimax [41, Chapter 5], Theorem 4.3 implies that ̂MMECov has minimax risk in the high-dimensional setting. It follows that, regardless of how well the ECov prior assumptions hold, ̂MMECov will not perform very poorly.
Further improvement with maximum marginal likelihood. The moment based approach analyzed above has a limitation: with positive probability, ⌃̂MM is not positive semi-definite (PSD). Though our expression for ̂MMECov remains well-defined in this case, this non-positive definiteness obscures the interpretation of ̂MMECov as a Bayes estimate. We next show that performance further improves if ⌃ is instead estimated by maximum marginal likelihood (Equation (1)) and is thereby constrained to be PSD.
Our next lemma characterizes the form of the resulting estimator, ̂ECov, and establishes a connection to the positive part James-Stein estimator [4].
Lemma 4.4. Assume D > Q and consider the singular value decomposition ̂LS = V diag( 1 2 )U> where V and U satisfy V >V = U>U = IQ, and is a Q-vector of non-negative reals. Under Condition 4.1, Equation (1) reduces to ⌃̂ = Udiag ⇥ (D 1 21Q)+ ⇤ U> and ̂ECov = V diag h 1 2 (1Q 2D 1)+ i U>, where (·)+ is shorthand for max(·, 0) element-wise, is the
Hadamard (i.e. element-wise) product, and the powers in 1 2 and 1 are applied element-wise.
Lemma 4.4 allows us to see ̂ECov as shrinking ̂LS toward 0 in the direction of each singular vector to an extent proportional to the inverse of the associated singular value. The transition from ̂MMECov to ̂ECov is then analogous to the taking the “positive part” of the James-Stein estimator in vector estimation, which provides a uniform improvement in risk [4]. Though R( , ̂ECov) is not easily available analytically, we nevertheless find that it dominates its moment-based counterpart.
Theorem 4.5. Assume D > Q + 1. Under Condition 4.1 ̂ECov dominates ̂MMECov with respect to squared error loss, achieving strictly lower risk for every value of .
We establish Theorem 4.5 using a proof technique adapted from Baranchik [4]; see also Lehmann and Casella [41][Thm. 5.5.4]. The standard approach we build upon is complicated by the fact that the directions in which we apply shrinkage are themselves random.
Theorem 4.5 provides a strong line of support for using ̂ECov over ̂MMECov that does not rely on any assumption of “correct” prior specification; in particular the risk improves without any subjective assumptions on . We discuss related earlier work in Appendix A.4.

heading: 5 Gains from ECov in the high-dimensional limit
text: The results of Section 4 give a promising endorsement of ECov but face two important limitations. First, the domination results relative to least squares do not directly demonstrate that ̂ECov attains improvements by leveraging similarities across groups in a meaningful way; indeed for a single group (i.e. Q = 1) ̂ECov can be understood as a ridge regression estimate [31], and Theorems 4.3 and 4.5 provide that ̂ECov dominates ̂LS for D > 3. Second, domination results reveal nothing about the size of the improvement or how it depends on any structure of ; intuitively, we should expect better performance when is in some way representative of the assumed prior. To address these limitations, we analyze the size of the gap between the risk of (1) ̂ECov and (2) our method applied to each group independently (ID), which we denote by ̂ID.1 We will characterize the dependence of this gap on .
Reasoning quantitatively about the dependence of the risk on the unknown parameter poses significant analytical challenges. In particular, Lemma 4.2 shows that R( , ̂MMECov) depends on through E[k̂†LSk2F | ]; however, k̂ † LSk2F is the sum of the eigenvalues of a non-central inverse Wishart matrix, a notoriously challenging quantity to work with; see e.g. [42, 30]. To regain tractability, we (1) develop an analysis asymptotic in the number of covariates D and (2) shift to a Bayesian analysis in order to sensibly consider a growing collection of covariate effects. In particular, we consider a sequence of regression problems, with parameters { d}1d=1 distributed as d
i.i.d.⇠ ⇡ for some distribution ⇡. Accordingly, instead of using the frequentist risk as in Section 4, we now use the Bayes risk to measure performance. Specifically, for a group with D covariates and an estimator ̂,
1 Our approach ̂ECov is well defined in the Q = 1 single group case; for each group q, we obtain ̂qID by computing ̂ECov on the group D = {(Xq, Y q)}.
the Bayes risk is RD⇡ (̂) := E⇡[R( , ̂)] where R( , ̂) is the usual frequentist risk. In the following, we describe the results of this analysis with proofs and additional details left to Appendix D.
For a single metric characterizing the benefits of joint modeling, we will define the asymptotic gain as the relative performance between our two estimators of interest here, ̂ECov and ̂ID. Definition 5.1. Consider a sequence of datasets of Q regression problems with an increasing number of covariates D, {DD}1D=1. Assume that for each group Condition 4.1 is satisfied with variance 2 and that each d
i.i.d.⇠ ⇡. The asymptotic gain of joint modeling is Gain(⇡, 2) := limD!1( 2DQ) 1[RD⇡ (̂ID) RD⇡ (̂ECov)].
The factor of 2DQ in Definition 5.1 puts Gain(⇡, 2) on a scale that is roughly invariant to the size and noise level of the problem; for example, ( 2DQ) 1RD⇡ (̂LS) = 1 for any ⇡, D, and Q. In Appendix D.5 we discuss how this asymptotic formulation may allow relaxation of Condition 4.1 if one considers certain random design matrices; for simplicity, the present analysis considers only fixed designs.
Our next lemma gives an analytic expression for Gain(⇡, 2) that provides a starting point for understanding its problem dependence.
Lemma 5.2. Assume ⌃̃ := Var⇡[ 1] is finite and has eigenvalues 1, . . . , Q. If Condition 4.1 satisfied asymptotically, Gain(⇡, 2) = 2Q 1[ PQ q=1( q + 2) 1 PQ q=1(⌃̃q,q + 2) 1].
Lemma 5.2 reveals that the diagonals and eigenvalues and ⌃̃ are key determinants of Gain(⇡, 2), but does not directly provide an interpretation of when ̂ECov offers benefits over ̂ID. Our next theorem demonstrates when an improvement can be achieved from joint modeling.
Theorem 5.3. Gain(⇡, 2) 0, with equality only when ⌃̃ = Var⇡[ 1] is diagonal.
Proof. From Lemma 5.2 we see Gain(⇡, 2) is the difference between a strictly Schur-convex function applied to the eigenvalues of ⌃̃ and to its diagonals (since (x+ 2) 1 is convex on R+). By the Schur-Horn theorem, the eigenvalues of ⌃̃ majorize its diagonals, providing the result.
Theorem 5.3 tells us that ̂ECov succeeds at adaptively learning and leveraging similarities among groups in the high-dimensional limit. In particular, Gain(⇡, 2) reduces to zero only when the eigenvalues of ⌃̃ are arbitrarily close to the entries of its diagonal, which occurs only when the covariate effects are uncorrelated across groups. However, when covariate effects are correlated, we obtain an improvement.
Our next theorem quantifies this relationship through upper and lower bounds.
Theorem 5.4. Let # and `# denote the eigenvalues and diagonals of ⌃̃, respectively, sorted in descending order. Then Gain(⇡, 2)  2 2Q 1k k2k`# #k2/( min + 2)3 and Gain(⇡, 2) 2Q 1k`# #k22/( max + 2)3, where max and min are the largest and smallest, respectively, eigenvalues of ⌃̃.
Theorem 5.4 allows us to see several aspects of when our method will and will not perform well. First, the presence of k`# #k22 in both the upper and lower bounds demonstrates that Gain(⇡, 2) will be small when the eigenvalues are close to the diagonal entries, with Euclidean distance as an informative metric.
As we find in our next corollary, Theorem 5.4 additionally allows us to see that nontrivial gains may be obtained only in an intermediate signal-to-noise regime, where signal is given by the size of the covariate effects and noise is the variance level 2. Notably, under Condition 4.1, 2 relates directly to the variance of ̂LS, and is influenced by both the residual variances and the group sizes; see Appendix C.1. In particular we interpret min as a proxy for signal strength since it captures the magnitude of typical d’s along their direction of least variation. Corollary 5.5. Gain(⇡, 2)  42 min/ 2 and Gain(⇡, 2)  42( min/ 2) 1, where  := max/ min is the condition number of ⌃̃.
Corollary 5.5 formalizes the intuitive result that with enough noise, the little recoverable signal is insufficient to effectively share strength. And furthermore, in the low-noise and high-signal regime ̂ID is very accurate on its own and there is little need for joint modeling. However, when there is a large gap between the largest and smallest eigenvalues of ⌃̃, leading  to be large, the gain could be larger.  will be large, for example, when the covariate effects are very correlated across groups.

heading: 6 Experiments

heading: 6.1 Simulated data
text: We first conduct simulations, where we can directly control the relatedness among groups and where we know the ground truth values of the parameters. We show that ECov is more accurate than EGroup when covariates outnumber groups, whether effects are correlated across groups or not.
In particular, we simulated covariates, parameters, and responses for Q = 10 groups across a range of covariate dimensions. We generated covariate effects as d
i.i.d.⇠ N (0,⌃). We chose ⌃ so that effects were either correlated (Figure 1 Left) or independent (Figure 1 Right) across groups; see Appendix E for details. We compare performance of six estimates on these groups. These are estimates assuming EGroup/ECov using moment matching and maximum marginal likelihood to choose ⌃/ (̂MMEGroup/̂ MM ECov and ̂EGroup/̂ECov, respectively), as well as least squares (̂LS), and
ECov applied to each group independently (̂ID).
Figure 1 reinforces our theoretical conclusions that (1) ̂ECov is more accurate when covariates outnumber groups and (2) ̂EGroup is more accurate when groups outnumber covariates. Our simulated X matrices are somewhat relaxed from a strict orthogonal design (Appendix E), so these experiments suggest that our conclusions hold beyond Condition 4.1. Additionally, ̂ECov and ̂EGroup both outperform their moment based counterparts, ̂MMECov and ̂ MM EGroup.
Even for the simulations with independent effects, Theorem 4.3 suggests ̂ECov should still outperform ̂LS and ̂EGroup in the higher dimensional regime, and we see this behavior in the right panel of Figure 1. Additionally, in agreement with Theorem 5.3, ̂ECov does not improve over ̂ID in the presence of independent effects, and the performances of these two estimators converge as D grows.

heading: 6.2 Real data
text: We find that ECov beats EGroup, as well as least squares and independent estimation, across three real groups. We describe the datasets (with additional details in Appendix E.4) and then our results.
Community level law enforcement in the United States. Policing rates vary dramatically across different communities, mediating disparate impacts of criminal law enforcement across racial and socioeconomic groups [64, 54]. Understanding how demographic and socioeconomic attributes of communities relate to variation in rates of law enforcement is crucial to understanding these impacts. Linear models provide the desired interpretability. We use a dataset [51] consisting of D = 117 community characteristics and their rates of law enforcement (per capita) for different crimes. We consider Q = 4 group subsets corresponding to distinct (region, crime) pairs: (Midwest, Robbery), (South, Assault), (Northeast, Larceny), and (West, Auto-theft). This data setup illustrates a small Q and accords with the independent residuals assumption in the likelihood shared by ECov and EGroup (Section 2). Across q, Nq represents between 400 and 600 communities.
Blog post popularity. We regress reader engagement (responses) on D = 279 characteristics of blog posts (covariates) [13]. We divided the corpus based on an included length attribute into Q = 3 groups, corresponding to (1) long posts, (2) short posts, and (3) posts from an earlier corpus with missing length attribute. We hypothesized that the relationships between the characteristics of posts and engagement would differ across these three groups. We randomly downsampled to Nq = 500 posts in each group to mimic a low sample-size regime, in which sharing strength is crucial.
Multiple binary classifications using pre-trained neural network embeddings on CIFAR10. Modern machine learning methods have proved very successful on large datasets. Translating this success to smaller datasets is one of the most actively pursued algorithmic challenges in machine learning. It has spurred the development of frameworks from transfer learning [65] to one-shot learning [62] to meta-learning [21]. One common and simple strategy starts with a learned representation (or “embedding”) from an expressive neural network fit to a large group. Then one can use this embedding as a covariate vector for classification tasks with few labeled data points.
We take a D = 128 dimensional embedding of the CIFAR10 image group [37, 60]. We create Q = 8 different binary classification tasks using the classes in CIFAR10 (Appendix E.4). We downsampled to Nq varying from 100 to 1000 to mimic a setting in which we hope to share strength from large groups to improve performance on smaller datasets.
Discussion of evaluation and results. In previous sections we have focused on parameter estimation. Here we instead evaluate with prediction error on held-out data since the true parameters are not observed. Specifically we perform 5-fold cross-validation and report the mean squared errors and classification errors on test splits. To reduce variance of out-of-sample error estimates on the applications in which we downsampled, we also evaluate on the additional held-out data. Because the residual variances were unknown, we estimated these for each application and group as ̂2q := kP?XqY qk2/(Nq D), where P?Xq := INq Xq(Xq>Xq) 1Xq > (see e.g. [23, Chapter 18.1]). All methods ran quickly on a 36 CPU machine; computation of ̂ECov, including the EM algorithm, required 2.04 ± 0.64, 6.89 ± 3.19 and 37.14 ± 3.39 seconds (mean ± st-dev across splits) on the law enforcement, blog, and CIFAR10 tasks, respectively.
Our results further reinforce the main aspects of our theory. ̂ECov outperformed ̂EGroup, independent Bayes estimates (̂ID), and least squares (̂LS) in all applications (at > 95% nominal confidence
with a paired t-test).2 Additionally, ̂ECov outperformed the baseline of ignoring heterogeneity, pooling groups together, and using the same effect estimates for every group (“Least Sqrs./MLE Pool”).
Appendix E includes additional results and comparisons. In particular, we provide the performance of the estimators on each component group for each application. Additionally, we report the performances of (1) stable and computationally efficient moment based alternatives to ̂ECov and ̂EGroup and (2) variants of ̂ECov and ̂EGroup that include a learned (rather than zero) prior mean. Appendix E.5 reports the licenses of software we used.

heading: 7 Discussion
text: The Bayesian community has long used hierarchical modeling with priors encoding exchangeability of effects across groups of data (EGroup). In the present work, we have made a case for instead using priors that encode exchangeability across covariates (ECov) – in particular, when the number of covariates exceeds the number of groups. We have presented a corresponding concrete model and inference method. We have shown that ECov outperforms EGroup in theory and practice when the number of covariates exceeds the number of groups.
Our approach is, of course, not a panacea. In some settings, a priori exchangeability among covariate effects will be inconsistent with prior beliefs. For example, imagine in the CIFAR10 application if meta-data covariates (such as geo-location and date) were available, in addition to embeddings. Then we might achieve better performance by treating meta-data covariates as distinct from embedding covariates. Additionally, we focused on a Gaussian prior for convenience. In cases where practitioners have more specific prior beliefs about effects, alternative priors and likelihoods may be warranted, though they may be more computationally challenging. Moreover, while relatively interpretable, linear models have their downsides. The linear assumption can be overly simplistic in many applications. It is common to misinterpret effects as causal rather than associative. Both the linear model and squared error loss lend themselves naturally to reporting means, but in many applications a median or other summary is more appropriate; so using a mean for convenience can be misleading.
Many exciting directions for further investigation remain. For example, the covariance ⌃ may provide an informative measure of task similarity; this similarity measure can be useful in, e.g., meta learning [34] and statistical genetics [12]. Additionally, we here explored two approaches to choosing the covariance matrices in the empirical Bayes step; more sophisticated approaches to covariance estimation may provided improved performance. It also remains to extend our methodology to other generalized linear models.
Acknowledgments and Disclosure of Funding
The authors thank Sameer K. Deshpande, Ryan Giordano, Alex Bloemendal, Lorenzo Masoero, and Diana Cai for insightful discussions on the manuscript, and the anonymous reviewers for their constructive suggestions. This work was supported in part by ONR Award N00014-18-S-F006 and an NSF CAREER Award. BLT is supported by NSF GRFP.

Review summaries:
The authors of this paper study exchangeable prior distribution for regression coefficients in high dimensional regression problems, where the number of covariates is typically much larger than the number of groups.  They argue that when the number of covariates are large, it is preferable to model the covariates as exchangeable rather than the datasets. The authors provide theoretical results by comparing frequentist risk of different empirical Bayes estimators and demonstrate the value of their approach empirically.

This paper studies a setting in which we analyze $Q$ datasets with $D$ covariates. The design matrix for the $q$-th dataset is denoted by $X^q$, a $N_q \times D$-dimensional matrix and the response is a $N_q$ dimensional vector denoted by $Y^q$. The question this paper addresses is the following. How can one share information across datasets and covariates to improve the predictive performance for each individual regression? The paper focuses mostly on the case of linear as well as logistic regression.  The motivation for the paper is that most existing approaches to the above problem posit exchangeability across datasets, for example, that the coefficient vector $\\beta^q \\in \\mathbb R^D$ in the $q$-th dataset, is drawn from the same prior as $\\beta^{q'}$ for another dataset $q'$. In contrast, the main suggestion in this paper is to assume exchangeability across covariates and that $(\\beta^{1}_j, ..., \\beta^{Q}_j) \\in \\mathbb R^Q$ (the vector that collects the coefficient of covariate $j$ across all datasets) is drawn from the same prior distribution for all covariates $j=1,...,D$. It is claimed, that the second approach may be preferable when the number of covariates is larger than the number of datasets.  To implement the abstract proposal above, a multivariate normal prior is assumed for  $(\\beta^{1}_j, ..., \\beta^{Q}_j)$ and an empirical Bayes approach is developed and the prior covariance matrix is estimated by marginal maximum likelihood through the EM algorithm.  There are real data analyses that demonstrate that the approach performs well. There is also theoretical justification. In one part of the theory, James-Stein-type dominance results are established (in a frequentist setting where the true coefficients are deterministic), and in the second part, some results are established in a Bayesian setting (where a prior exists and is well-specified) with asymptotics as $D \to \infty$.

This paper argues that when the number of covariates exceeds the number of datasets, it is more natural to model effects as exchangeable across covariates than across datasets. The paper illustrates the setting using a microcredit example.   It presents a hierarchical model to formalize the perspective and derive an efficient algorithm. The paper evaluates the model using three real-world datasets.