# Improved methods for empirical Bayes multivariate multiple testing and
effect size estimation

Yunqi Yang  Peter Carbonetto  David Gerard  Matthew
<EMAIL_ADDRESS>[ Committee on Genetics, Genomics and
System Biology, University of Chicago  Department of Human Genetics,
University of Chicago  Department of Mathematics and Statistics, American
University  Department of Statistics, University of Chicagopresep=, ]e1

###### Abstract

Estimating the sharing of genetic effects across different conditions is
important to many statistical analyses of genomic data. The patterns of
sharing arising from these data are often highly heterogeneous. To flexibly
model these heterogeneous sharing patterns, Urbut et al. (2019) proposed the
multivariate adaptive shrinkage (MASH) method to jointly analyze genetic
effects across multiple conditions. However, multivariate analyses using MASH
(as well as other multivariate analyses) require good estimates of the sharing
patterns, and estimating these patterns efficiently and accurately remains
challenging. Here we describe new empirical Bayes methods that provide
improvements in speed and accuracy over existing methods. The two key ideas
are: (1) adaptive regularization to improve accuracy in settings with many
conditions; (2) improving the speed of the model fitting algorithms by
exploiting analytical results on covariance estimation. In simulations, we
show that the new methods provide better model fits, better out-of-sample
performance, and improved power and accuracy in detecting the true underlying
signals. In an analysis of eQTLs in 49 human tissues, our new analysis
pipeline achieves better model fits and better out-of-sample performance than
the existing MASH analysis pipeline. We have implemented the new methods,
which we call “Ultimate Deconvolution”, in an R package, udr, available on
GitHub.

multivariate analysis,

empirical Bayes,

normal means,

multiple testing,

adaptive shrinkage,

expectation maximization,

###### keywords:

, and

##  1 Introduction

The problem of testing and estimating effect sizes for many units in multiple
conditions (or on multiple outcomes) arises frequently in genomics
applications. Examples include assessing the effects of many expression
quantitative trait loci (eQTLs) in multiple tissues (GTEx Consortium et al.,
2015) and assessing the effects of many genetic variants on multiple traits
(Zhou and Stephens, 2014; Pickrell et al., 2016; Turchin and Stephens, 2019;
Udler et al., 2018; Zou et al., 2024). The simplest approach to assessing
effects in multiple conditions is to analyze each condition separately.
However, this fails to exploit sharing or similarity of effects among
conditions. For example, a genetic variant that increases expression of a
particular gene in the heart may similarly increase expression in other
tissues, particularly in tissues that are related to heart. Such sharing of
effects could be exploited to improve power and estimation accuracy by
borrowing information across conditions. This is, in essence, the motivation
for meta-analysis methods (Willer, Li and Abecasis, 2010; Han and Eskin, 2011;
Wen and Stephens, 2014), and more generally it motivates consideration of
multivariate approaches to multiple testing and effect size estimation (Urbut
et al., 2019).

While borrowing information across conditions may be a natural idea, getting
it to work well in practice requires confronting some thorny issues. One
challenge is that the extent to which effects are shared among conditions will
vary among data sets; for example, data sets involving very similar conditions
may have very high levels of sharing, whereas data sets involving very
different conditions may exhibit little to no sharing. Furthermore, some data
sets may include some conditions that are very similar and others that are
very dissimilar, and these differences may be difficult to specify in advance.
Assessing the sharing and similarity of effects among conditions may also be
an important goal in itself.

Empirical Bayes (EB) approaches (e.g., Flutre et al. 2013; Urbut et al. 2019)
provide an attractive way to confront these challenges. EB methods estimate a
prior distribution that captures the sharing or similarity of effects among
the conditions, then, using Bayes theorem, they combine the prior with the
observed data to improve the effect estimates. The methods proposed in Urbut
et al. (2019) and implemented in the R package mashr assume a mixture of
multivariate normal distributions for the prior, which has the twin advantages
of being both both flexible and computationally convenient. Indeed, mashr has
been used to analyze very large data sets involving many conditions (Zou et
al., 2024; Lin et al., 2024; Urbut et al., 2021; Barbeira et al., 2020; Araujo
et al., 2023; Li et al., 2022; Soliai et al., 2021; Natri et al., 2024; Bonder
et al., 2021).

Despite this, these methods still have considerable limitations; in
particular, fitting a mixture of multivariate normals prior with unknown
covariance matrices raises statistical and computational challenges. Urbut et
al. (2019) used a two-stage procedure to deal with (or sidestep) these
challenges: in the first stage, the covariance matrices in the prior were
estimated by maximum-likelihood on a subset of the data (using the “Extreme
Deconvolution” algorithm of Bovy, Hogg and Roweis 2011); next, given the
covariances estimated in the first stage, the second stage estimated the
mixture proportions in the prior by maximzing the likelihood from all the data
(using the fast optimization algorithms of Kim et al. 2020). The second stage
is a convex optimization problem, and can be solved efficiently and reliably
for very large data sets. The first stage, however, presents several
challenges, including that the Extreme Deconvolution (ED) algorithm can be
slow to converge, the results of running ED are often sensitive to
initialization, and the estimated covariance matrices can be quite unstable,
particularly when the number of conditions, $R$, is large relative to sample
size, $n$. These challenges motivated this work, which is a closer examination
of these challenges from both a statistical perspective and a computational
one. One of the contributions of this study is a new algorithm, which we call
“truncated eigenvalue decomposition” (TED). TED often converges much faster
than ED. (Noting that ED applies to some settings where TED does not.) We also
explore the use of simple regularization schemes that can improve accuracy
compared with maximum-likelihood estimation, particularly when the sample size
$n$ is small and the number of conditions $R$ is large (that is, the ratio
$R/n$ is large). And we highlight some problems that arise from using low-rank
covariance matrices, which was a strategy previously suggested in Urbut et al.
(2019) to reduce the number of estimated parameters. We provide an R package,
udr (“Ultimate Deconvolution in R”), available at
https://github.com/stephenslab/udr, which implements all these methods
described here within a convenient, user-friendly interface, and that
interacts well with our previous R package, mashr (Urbut et al., 2019). The
version of the R package used to perform the analyses in this manuscript is
also available for download as a supplementary file (Yang, Carbonetto and
Stephens, 2024a).

The paper is structured as follows. Section 2 summarizes notation used in the
mathematical expressions throughout the paper. Section 3 formally introduces
the models, priors and regularization schemes considered, and Section 4
describes the model-fitting algorithms and procedures. Section 5 evaluates the
performance of the different methods and algorithms on data sets simulated
under a variety of scenarios. Section 6 applies the improved methods to an
analysis of a large, multi-tissue eQTL data set from the GTEx Project.
Finally, Section 7 discusses the promise and limitations of our methods, and
future directions.

##  2 Notation used in the mathematical expressions

For the mathematical expressions below, we use bold, capital letters
($\bm{A}$) to denote matrices; bold, lowercase letters ($\bm{a}$) to denote
column vectors; and plain, lowercase letters ($a$) to denote scalars. For a
matrix ${\bm{A}}$, $\|\bm{A}\|_{\ast}$ denotes its nuclear norm, $\bm{A}^{T}$
denotes its transpose, $\bm{A}^{-1}$ denotes its inverse, $\bm{A}^{-T}$
denotes the inverse of $\bm{A}^{T}$, $\text{tr}(\bm{A})$ denotes the trace,
and $|\bm{A}|$ denotes the matrix determinant. We use $\bm{0}$ and $\bm{1}$ to
denote the vectors whose elements are all zeros and all ones respectively;
$\bm{e}_{r}$ for the standard basis vector
$\bm{e}_{r}=(0,\dots,0,1,0,\dots,0)$ with the 1 appearing in the $r$th
position; $\bm{I}_{R}$ is the $R\times R$ identity matrix; and
$\mbox{diag}(\bm{a})$ denotes the diagonal matrix in which the diagonal
entries are given by the entries of vector $\bm{a}$. We use $\mathbb{R}$ for
the set of real numbers, $\mathbb{R}^{R}$ for the set of real-valued vectors
of length $R$, and $\mathbb{R}^{n\times m}$ for the set of real-valued
$n\times m$ matrices. We use $N(\mu,\sigma^{2})$ to denote the normal
distribution on $\mathbb{R}$ with mean $\mu$ and variance $\sigma^{2}$, and
$N(\,\cdot\,;\mu,\sigma^{2})$ denotes its density. And we use
$N_{R}(\bm{\mu},\bm{\Sigma})$ for the multivariate normal distribution on
$\mathbb{R}^{R}$ with mean $\bm{\mu}\in\mathbb{R}^{R}$ and $R\times R$
covariance matrix $\bm{\Sigma}$, and $N_{R}(\,\cdot\,;\bm{\mu},\bm{\Sigma})$
denotes its density. We use $P^{+}_{R}\subseteq\mathbb{R}^{R\times R}$ to
denote the set of real-valued $R\times R$ (symmetric) positive semi-definite
matrices. $\mathcal{S}_{R}\subseteq\mathbb{R}^{R}$ denotes the $R$-dimensional
simplex.

##  3 The empirical Bayes multivariate normal means model

In this section, we define the “empirical Bayes multivariate normal means”
(EBMNM) model. We describe several variations of this model that involve
constraints or penalties on the model parameters.

The EBMNM model assumes that we observe vectors $\bm{x}_{j}\in\mathbb{R}^{R}$,
that are independent, noisy, normally-distributed measurements of underlying
true values $\bm{\theta}_{j}\in\mathbb{R}^{R}$:

$\bm{x}_{j}\mid\bm{\theta}_{j}\sim N_{R}(\bm{\theta}_{j},\bm{V}_{j}),\quad
j=1,\ldots,n,$ (1)

in which the covariances $\bm{V}_{j}\in P^{+}_{R}$ are assumed to be known and
invertible. Our ultimate goal is to perform inference for the unknown means
$\bm{\theta}_{j}$ from the observed data $\bm{x}_{j}$. This model is a natural
generalization of the well-studied (univariate) normal means model (Robbins,
1951; Efron and Morris, 1972; Stephens, 2017; Bhadra et al., 2019; Johnstone,
2019; Sun, 2020), and so we call it the “multivariate normal means model”. An
important special case is when the measurement error distribution is the same
for all observations; i.e., $\bm{V}_{j}=\bm{V}$, $j=1,\ldots,n$. We refer to
this as the “homoskedastic” case. Some of our computational methods are
designed specifically for the homoskedastic case, which is easier to solve
than the “heteroskedastic” case.

The EBMNM model (1) further assumes that the unknown means are independent and
identically distributed draws from some distribution, which we refer to as the
“prior distribution”. While other choices are possible, here we assume the
prior distribution is a mixture of zero-mean multivariate normals:

$p(\bm{\theta}_{j}\mid\bm{\pi},\mathcal{U})=\sum_{k=1}^{K}\pi_{k}N_{R}(\bm{\theta}_{j};\bm{0},\bm{U}_{k}),$
(2)

where $\bm{\pi}\in\mathcal{S}_{K}$ is the set of mixture proportions,
$\bm{U}_{k}\in P_{R}^{+}$ are covariance matrices, and
$\mathcal{U}\colonequals\\{\bm{U}_{1},\dots,\bm{U}_{K}\\}$ denotes the full
collection of covariance matrices.

Combining (1) and (2) yields the marginal distribution

$p(\bm{x}_{j}\mid\bm{\pi},\mathcal{U})=\sum_{k=1}^{K}\pi_{k}N_{R}(\bm{x}_{j};\bm{0},\bm{U}_{k}+\bm{V}_{j}),$
(3)

and the marginal log-likelihood,

$\displaystyle\ell(\bm{\pi},\mathcal{U})\colonequals$
$\displaystyle\;\sum_{j=1}^{n}\log p(\bm{x}_{j}\mid\bm{\pi},\mathcal{U})$
$\displaystyle=$
$\displaystyle\;\sum_{j=1}^{n}\log\sum_{k=1}^{K}\pi_{k}N_{R}(\bm{x}_{j};\bm{0},\bm{U}_{k}+\bm{V}_{j})$
(4)

The EB approach to fitting the model (1–2) proceeds in two stages:

  1. 1.

Estimate the prior parameters $\bm{\pi},\mathcal{U}$ by maximizing a penalized
log-likelihood,

$(\hat{\bm{\pi}},\hat{\mathcal{U}})\colonequals\operatorname*{argmax}_{\bm{\pi}\,\in\,\mathcal{S}_{K},\,\bm{U}_{k}\,\in\,P^{+,k}_{R}}\ell(\bm{\pi},\mathcal{U})-\sum_{k=1}^{K}\tilde{\rho}(\bm{U}_{k}),$
(5)

where $\tilde{\rho}$ denotes a penalty function introduced to regularize
$\bm{U}_{k}$, and $P^{+,k}_{R}\subseteq P^{+}_{R}$ allows for constraints on
$\bm{U}_{k}$ (which may be different for each component of the mixture). The
exact penalties and constraints considered are described below. When
$P^{+,k}_{R}=P^{+}_{R}$ and $\tilde{\rho}(\bm{U}_{k})=0$ for all $\bm{U}_{k}$,
solving (5) corresponds to maximum-likelihood estimation of
$\bm{\pi},\mathcal{U}$.

  2. 2.

Compute the posterior distribution for each $\bm{\theta}_{j}$ given
$\bm{\pi},\bm{U}$ estimated in the first stage:

$\displaystyle p_{\mathsf{post}}({\bm{\theta}}_{j})\colonequals$
$\displaystyle\;p(\bm{\theta}_{j}\mid\bm{x}_{j},\hat{\bm{\pi}},\hat{\mathcal{U}})$
$\displaystyle\propto$
$\displaystyle\;p(\bm{x}_{j}\mid\bm{\theta}_{j})\,p(\bm{\theta}_{j}\mid\hat{\bm{\pi}},\hat{\mathcal{U}}).$
(6)

The posterior distributions (2) have an analytic form, and are mixtures of
multivariate normal distributions (see for example Urbut et al. 2019; Bovy,
Hogg and Roweis 2011) and Section LABEL:appendix:posterior_theta of the
Supplementary Materials; Yang, Carbonetto and Stephens 2024b). Since these
posterior distributions have a closed form, it is relatively straightforward
to compute posterior summaries such as the posterior mean and posterior
standard deviation, and measures for significance testing such as the local
false sign rate (lfsr) (Stephens, 2017). Therefore, we focus on methods for
accomplishing the first step, solving (5).

###  3.1 A reformulation to ensure scale invariance

Intuitively, one might hope that changing the units of measurement of all the
observed $\bm{x}_{j}$ would simply result in corresponding changes to the
units of the estimated $\bm{\theta}_{j}$. This idea can be formalized as
requiring that solutions of the EBMNM model should obey a “scale invariance”
property. This consideration motivates us to reformulate (5). Let
$\hat{\bm{\theta}}_{j}(\bm{x}_{1},\ldots,\bm{x}_{n},\bm{V}_{1},\ldots,\bm{V}_{n})$
denote the posterior mean for $\bm{\theta}_{j}$ computed by solving the EBMNM
problem with data $\bm{x}_{1},\ldots,\bm{x}_{n},\bm{V}_{1},\ldots,\bm{V}_{n}$.
We say that the solution is “scale invariant” if, for any $s>0$, the following
holds:

$\hat{\bm{\theta}}_{j}(s\bm{x}_{1},\ldots,s\bm{x}_{n},s^{2}\bm{V}_{1},\ldots,s^{2}\bm{V}_{n})=s\hat{\bm{\theta}}_{j}(\bm{x}_{1},\ldots,\bm{x}_{n},\bm{V}_{1},\ldots,\bm{V}_{n}).$
(7)

That is, multiplying all the observed data points by $s$ (which multiplies the
corresponding error variance by $s^{2}$) has the effect of multiplying the
estimated means by $s$. (Note: we state scale invariance in terms of the
posterior means only for simplicity; the concept is easily generalized to
require that the whole posterior distribution for $\bm{\theta}_{j}$ scales
similarly.)

Without the penalty $\tilde{\rho}$ in (5), it is easy to show that the scale
invariance property (7) holds provided that the constraints satisfy
$\bm{U}_{k}\in P^{+,k}_{R}\implies s_{k}\bm{U}_{k}\in P^{+,k}_{R},\forall
s_{k}>0$. With penalty, scale invariance holds provided that
$\tilde{\rho}(\bm{U})=\tilde{\rho}(s\bm{U}),\forall s,\bm{U}$; that is,
provided that the penalty function depends only on the “shape” of $\bm{U}$ and
not on its “scale”. To ensure scale invariance, we therefore consider
penalties of the form

$\tilde{\rho}(\bm{U})=\min_{s\,>\,0}\,\rho(\bm{U}/s),$ (8)

where $\rho$ is a penalty that may depend on both the shape and scale of
$\bm{U}$. To give some intuition, suppose that $\rho$ encourages all the
eigenvalues of $\bm{U}$ to be close to 1. Then $\tilde{\rho}$ will encourage
the eigenvalues to be close to each other, without requiring that they
specifically be close to 1. Therefore, plugging (8) into (5), we have

$(\hat{\bm{\pi}},\hat{\mathcal{U}})\colonequals\operatorname*{argmax}_{\bm{\pi}\,\in\,\mathcal{S}_{K},\,\bm{U}_{k}\,\in\,P^{+,k}_{R},\,\bm{s}\,>\,\bm{0}}\ell(\bm{\pi},\mathcal{U})-\sum_{k=1}^{K}\rho(\bm{U}_{k}/s_{k}),$
(9)

where $\bm{s}\colonequals(s_{1},\ldots,s_{K})$. We use (9) for the remainder
of the paper.

###  3.2 Constraints and penalties

Estimating covariance matrices in high-dimensional settings (i.e., large $R$)
is known to be a challenging problem (e.g., Johnstone and Paul 2018; Fan, Liao
and Liu 2016; Ledoit and Wolf 2022). Even in the simpler setting of
independent and identically distributed observations from a single
multivariate normal distribution, the maximum-likelihood estimate of the
covariance matrix can be unstable, and so various covariance regularization
approaches have been proposed to address this issue (Ledoit and Wolf, 2004;
Friedman, Hastie and Tibshirani, 2008; Cai and Liu, 2011; Won et al., 2013;
Chi and Lange, 2014). Interestingly, in the context of using EBMNM for
significance testing, adding penalties have additional benefits (see the
numerical experiments below).

We consider two different penalties that have been previously used for
covariance regularization:

  1. 1.

The “inverse Wishart” (IW) penalty

$\displaystyle\rho_{\lambda}^{\mathrm{IW}}(\bm{U})$
$\displaystyle\colonequals\frac{\lambda}{2}\\{\log|\bm{U}|+\mathrm{tr}({\bm{U}}^{-1})\\}$
(10) $\displaystyle=\frac{\lambda}{2}\sum_{r=1}^{R}\big{(}\log
e_{r}+1/e_{r}\big{)}.$ (11)

  2. 2.

The “nuclear norm” (NN) penalty:

$\displaystyle\rho_{\lambda}^{\mathrm{NN}}(\bm{U})$
$\displaystyle\colonequals\frac{\lambda}{2}\\{0.5\|\bm{U}\|_{\ast}+0.5\|\bm{U}^{-1}\|_{\ast}\\}$
(12)
$\displaystyle=\frac{\lambda}{2}\sum_{r=1}^{R}\big{(}0.5e_{r}+0.5/e_{r}\big{)}.$
(13)

Here, $e_{1},\dots,e_{R}$ denote the eigenvalues of $\bm{U}$, and $\lambda>0$
controls the strength of the penalty. We chose $\lambda=R$ in our simulations,
but one could also use cross-validation to select $\lambda$.

The IW penalty on $\bm{U}_{k}$ corresponds to maximum a posteriori (MAP)
estimation of $\bm{U}_{k}$ under an inverse-Wishart prior, with prior mode
$\bm{I}_{R}$ and $\lambda-R-1$ degrees of freedom (Fraley and Raftery, 2007).
This penalty was also mentioned (but not used or evaluated) in Bovy, Hogg and
Roweis (2011).

The nuclear norm penalty was studied as an alternative to the IW penalty for
estimation of covariance matrices in Chi and Lange (2014). To our knowledge,
this penalty has not been studied in the EBMNM setting. The nuclear norm
penalty in Chi and Lange (2014) includes a hyperparameter $\alpha\in(0,1)$
that controls the balance between $\|\bm{U}\|_{\ast}$ and
$\|\bm{U}^{-1}\|_{\ast}$, but our approach to ensuring scale-invariance has
the consequence that changing $\alpha$ is equivalent to changing $\lambda$
(see the Supplementary Materials; Yang, Carbonetto and Stephens 2024b), so we
set $\alpha=0.5$.

As can be seen from (11) and (13), both penalties depend only on the
eigenvalues of $\bm{U}$, and decompose into additive functions of the $R$
eigenvalues. Both penalties are minimized when $\bm{U}=\bm{I}_{R}$, and more
generally encourage $\bm{U}$ to be well-conditioned by making it closer to the
identity matrix (by pushing the eigenvalues closer to 1).

As an alternative to penalized estimation of $\bm{U}$, we also consider
estimating $\bm{U}$ under different constraints:

  1. 1.

A scaling constraint, $\bm{U}_{k}=c_{k}\bm{U}_{0k}$, for some chosen
$\bm{U}_{0k}\in P^{+}_{R}$.

  2. 2.

A rank-1 constraint, $\bm{U}_{k}=\bm{u}_{k}\bm{u}_{k}^{T}$, for some
$\bm{u}_{k}\in\mathbb{R}^{R}$.

The scaling constraint could be useful for situations in which the
$\bm{\theta}_{j}$ may obey an expected sharing structure, or for sharing
structures that are easier to interpret. For example, $\bm{U}_{0k}=\bm{I}_{R}$
captures the situation in which all the effects
$\theta_{j1},\dots,\theta_{jR}$ are independent, and
$\bm{U}_{0k}=\bm{1}\bm{1}^{T}$ captures the situation in which all the effects
$\theta_{j1},\dots,\theta_{jR}$ are equal. Such covariances are referred to as
“canonical” covariance matrices in Urbut et al. (2019).

The rank-1 constraint also leads to potentially more interpretable covariance
matrices, and can be thought of as a form of regularization because low-rank
matrices have fewer parameters to be estimated. It may also allow for faster
computations. Urbut et al. (2019) in fact make extensive use of the rank-1
constraint. However, our results will show that this constraint can cause
problems for significance testing and so may be better avoided in practice.

##  4 Fitting the EBMNM model

We now describe three algorithms we have implemented for fitting variations of
the EBMNM model described above: the ED algorithm from Bovy, Hogg and Roweis
2011; an algorithm based on methods commonly used in factor analysis (FA); and
another based on the truncated eigenvalue decomposition (TED). Each of these
algorithms applies to a subset of EBMNM models (Table 1). In some situations,
only one algorithm can be applied; for example, only ED can handle
heteroskedastic variances with no constraints on $\bm{U}$. However, in other
settings all three algorithms can be applied (e.g., homoskedastic errors, no
constraints on $\bm{U}$). Below, we empirically assess the relative merits of
the different algorithms in simulations that model these different settings.
See also Supplementary Tables LABEL:tab:comp_complex1 and
LABEL:tab:comp_complex2 (Yang, Carbonetto and Stephens, 2024b) for a
comparison of the algorithms’ computational properties in the different
settings.

Table 1: Summary of the EBMNM algorithms and the situations in which they apply. In this table we consider 4 variations of EBMNM: no constraint on $\bm{U}$ or a rank-1 constraint; and homoskedastic (hom.) errors (all $\bm{V}_{j}$ are the same) or heteroskedastic (het.) errors (one or more $\bm{V}_{j}$ differ). Abbreviations used in this table are: ED = Extreme Deconvolution), FA = factor analysis, TED = truncated eigenvalue decomposition. A checkmark ($\checkmark$) indicates that the algorithm (ED, FA, TED) can be applied to the particular variation. Note that fitting $\bm{U}$ with a scaling constraint involves only a 1-d optimization and is treated separately. | constraints on $\bm{U}$
---|---
| none | rank-1
algorithm | hom. | het. | hom. | het.
ED | $\checkmark$ | $\checkmark$ |  |
FA | $\checkmark$ |  | $\checkmark$ | $\checkmark$
TED | $\checkmark$ |  | $\checkmark$ |

###  4.1 Algorithms for the single-component EBMNM model with no penalty

While these algorithms are implemented for the mixture prior (2) with the
penalties described above, the algorithms are much easier to describe in the
special case of one mixture component ($K=1$), and without penalty. So
initially we focus on this simpler case, and later we extend to the general
form with penalties and $K\geq 1$.

With $K=1$ and no penalty, the prior is $\bm{\theta}_{j}\sim
N(\bm{0},\bm{U})$, the model is

$\bm{x}_{j}\mid\bm{U}\sim N_{R}(\bm{0},\bm{U}+\bm{V}_{j}),\quad j=1,\dots,n.$
(14)

and the goal is to compute the maximum-likelihood estimate of $\bm{U}$:

$\hat{\bm{U}}\colonequals\operatorname*{argmax}_{\bm{U}\,\in\,P_{R}^{+}}\sum_{j=1}^{n}\log
N_{R}(\bm{x}_{j};\bm{0},\bm{U}+\bm{V}_{j})$ (15)

The three algorithms for solving (15) are as follows.

##### Truncated Eigenvalue Decomposition (TED)

This algorithm, which to the best of our knowledge is new in this context,
exploits the fact that, in the special case where $\bm{V}_{j}=\bm{I}_{R}$,
$j=1,\ldots,n$, the maximum-likelihood estimate (15) can be computed exactly.
At first glance, one might try to solve for $\hat{\bm{U}}$ by setting
$\hat{\bm{U}}+\bm{I}_{R}$ to the sample covariance matrix,
$\bm{S}\colonequals\sum_{j=1}^{n}\bm{x}_{j}\bm{x}_{j}^{T}/n$, then recovering
$\hat{\bm{U}}$ as $\hat{\bm{U}}=\bm{S}-\bm{I}_{R}$. However,
$\bm{S}-\bm{I}_{R}$ is not necessarily a positive semi-definite matrix; that
is, it may have one or more eigenvalues that are negative. One could deal this
problem by setting the negative eigenvalues to zero, and indeed this approach
is correct; that is, letting $(\bm{S})_{+}$ be the matrix obtained from
$\bm{S}$ by truncating its negative eigenvalues to 0,
$\hat{\bm{U}}=(\bm{S}-\bm{I}_{R})_{+}$ is the maximum-likelihood estimate of
$\bm{U}$ (Tipping and Bishop, 1999). This idea can also be used to solve the
more general case, $\bm{V}_{j}=\bm{V}$, essentially by transforming the data
to $\bm{R}^{-1}\bm{x}_{j}$ where $\bm{V}=\bm{R}\bm{R}^{T}$, estimating
$\bm{R}^{-1}\bm{U}\bm{R}^{-T}$ from this transformed data, then reversing this
transformation, see Yang, Carbonetto and Stephens 2024b.

##### Extreme Deconvolution (ED)

This is an EM algorithm (Dempster, Laird and Rubin, 1977), an iterative
approach to solving (15); the name comes from Bovy, Hogg and Roweis (2011). ED
uses the natural “data augmentation” representation of (14):

$\displaystyle\bm{\theta}_{j}$ $\displaystyle\sim N_{R}(\bm{0},\bm{U})$ (16)
$\displaystyle\bm{x}_{j}\mid\bm{\theta}_{j}$ $\displaystyle\sim
N_{R}(\bm{\theta}_{j},\bm{V}_{j}).$

Following the usual EM derivation, the updates can be derived as

$\bm{U}^{\mathrm{new}}=\frac{1}{n}\sum_{j=1}^{n}\bm{B}_{j}+\bm{b}_{j}\bm{b}_{j}^{T},$
(17)

where $\bm{b}_{j}$ and $\bm{B}_{j}$ are, respectively, the posterior mean and
covariance of $\bm{\theta}_{j}$ given $\bm{U}$:

$\displaystyle\bm{b}_{j}$
$\displaystyle\colonequals\bm{U}(\bm{U}+\bm{V}_{j})^{-1}\bm{x}_{j}$ (18)
$\displaystyle\bm{B}_{j}$
$\displaystyle\colonequals\bm{U}-\bm{U}(\bm{U}+\bm{V}_{j})^{-1}\bm{U}.$ (19)

The update (17) is guaranteed to increase (or not decrease) the objective
function in (15), and repeated application of (17–19) will converge to a
stationary point of the objective.

##### Factor analysis (FA)

This is also an EM algorithm, but based on a different data augmentation than
ED; the name comes from its close connection to EM algorithms for factor
analysis models (Ghahramani and Hinton, 1996; Rubin and Thayer, 1982; Zhao,
Philip and Jiang, 2008; Liu and Rubin, 1998; McLachlan and Peel, 2000). In its
simplest form, the FA approach imposes a rank-1 constraint on $\bm{U}$,
$\bm{U}=\bm{u}\bm{u}^{T}$, where $\bm{u}\in\mathbb{R}^{R}$ is to be estimated.
The model (14) then admits the following data augmentation representation:

$\displaystyle a_{j}$ $\displaystyle\sim N(0,1)$ (20)
$\displaystyle\bm{x}_{j}\mid\bm{u},\bm{V}_{j},a_{j}$ $\displaystyle\sim
N(a_{j}\bm{u},\bm{V}_{j}).$

The usual EM derivation gives the update

$\bm{u}^{\mathrm{new}}=\bigg{(}\sum_{j=1}^{n}(\mu_{j}^{2}+\sigma_{j}^{2})\bm{V}_{j}^{-1}\bigg{)}^{-1}\bigg{(}\sum_{j=1}^{n}\mu_{j}\bm{V}_{j}^{-1}\bm{x}_{j}\bigg{)},$
(21)

in which $\mu_{j}$ and $\sigma_{j}^{2}$ denote, respectively, the posterior
mean and posterior covariance of $a_{j}$ given $\bm{u}$,

$\displaystyle\mu_{j}$
$\displaystyle\colonequals\sigma_{j}^{2}\bm{u}^{T}\bm{V}_{j}^{-1}\bm{x}_{j}$
(22) $\displaystyle\sigma_{j}^{2}$ $\displaystyle\colonequals
1/(1+\bm{u}^{T}\bm{V}_{j}^{-1}\bm{u}).$ (23)

The update (21) is guaranteed to increase (or not decrease) the objective
function in (15), and repeated application of (21–23) will converge to a
stationary point of the objective. These updates can be extended to higher-
rank covariances where the goal is to find the maximum-likelihood estimate
subject of $\bm{U}$ subject to $\bm{U}$ having rank at most $R^{\prime}$,
where $R^{\prime}\leq R$. However, when $R^{\prime}>1$, the updates are have
closed-form expressions only for homoskedastic errors, $\bm{V}_{j}=\bm{V}$.

The three algorithms have different strengths and weaknesses, and settings to
which they apply (Table 1). The TED algorithm has the advantage of directly
computing the maximum-likelihood estimate, which seems preferable to an
iterative approach. However, TED only applies in the case of homoskedastic
errors. The ED approach is more general, applying to both heteroskedastic and
homoskedastic errors, although it cannot fit rank-1 matrices. The FA approach
is attractive for low-rank matrices, particularly for fitting rank-1 matrices.

Our claim that ED cannot fit rank-1 covariance matrices deserves discussion,
especially since Urbut et al. (2019) used ED to fit such matrices. As pointed
out by Urbut et al. (2019), if ED is initialized to a low-rank matrix with
rank $R^{\prime}$, then the updated estimates (17) are also rank (at most)
$R^{\prime}$. Thus, if ED is initalized to a rank-1 matrix, the final estimate
is also rank-1. However, the ED estimates are not only low rank, but also span
the same subspace as the initial estimates, a property we refer to as
“subspace-preserving”. Thus, if ED is initialized to
$\bm{U}=\bm{u}\bm{u}^{T}$, the updated estimates will be of the form
$a\bm{u}\bm{u}^{T}$ for some scalar $a$. In other words, the ED update does
not change rank-1 matrices, except by a scaling factor, and so the final
estimate will simply be proportional to the initial estimate. This flaw
motivated us to implement the FA method. However, as our numerical comparisons
will illustrate, the rank-1 matrices turn out to have other drawbacks that
lead us not to recommend their use anyhow.

###  4.2 Extending the algorithms to a mixture prior

Input: Data vectors $\bm{x}_{j}\in\mathbb{R}^{R}$ and correspoding covariance
matrices $\bm{V}_{j}\in P_{R}^{+}$, $j=1,\ldots,n$; $K$, the number of mixture
components; initial estimates of the prior covariance matrices
$\mathcal{U}=\\{\bm{U}_{1},\ldots,\bm{U}_{K}\\}$, $\bm{U}_{k}\in P_{R}^{+,k}$,
$k=1,\ldots,K$; initial estimates of the scaling parameters
$\bm{s}=\\{s_{1},\ldots,s_{K}\\}\in\mathbb{R}^{K}$; initial estimates of the
mixture weights $\bm{\pi}\in\mathcal{S}_{K}$.

Output: $\mathcal{U}$, $\bm{\pi}$.

repeat

for  _$j\leftarrow 1$ to $n$_ do

for  _$k\leftarrow 1$ to $K$_ do

Update $w_{jk}$ using (24).

end for

end for

for  _$k\leftarrow 1$ to $K$_ do

$\pi_{k}\leftarrow\sum_{j=1}^{n}w_{jk}/n$

$\bm{U}_{k}\leftarrow\operatorname*{argmax}_{\bm{U}\,\in\,P_{R}^{+,k}}\phi(\bm{U};\bm{w}_{k})-\rho(\bm{U}/s_{k})$

$\triangleright$ Note that some algorithms compute this argmax inexactly.

$s_{k}\leftarrow\operatorname*{argmin}_{s\,>\,0}\rho(\bm{U}_{k}/s)$

end for

until  _convergence criterion is met_ ;

Algorithm 1 EM for fitting the EBMNM model.

All of the algorithms—TED, ED and FA—can be generalized from the $K=1$ case to
the $K\geq 1$ case using the standard EM approach to dealing with mixtures.
The resulting algorithms have a simple common structure, summarized in
Algorithm 1. (This algorithm also allows for an inclusion of a penalty, which
is treated in the next section. See also the Supplementary Materials for a
derivation of this algorithm; Yang, Carbonetto and Stephens 2024b.) This EM
algorithm involves iterating the following steps: (i) compute the weights,
$w_{jk}$ (sometimes called the “responsibilities”), each which represent the
conditional probability that mixture component $k$ gave rise to observation
$j$ given the current estimates of $\bm{\pi},\mathcal{U}$,

$w_{jk}=\frac{\pi_{k}N_{R}(\bm{x}_{j};\bm{0},\bm{U}_{k}+\bm{V}_{j})}{\sum_{k^{\prime}=1}^{K}\pi_{k^{\prime}}N_{R}(\bm{x}_{j};\bm{0},\bm{U}_{k^{\prime}}+\bm{V}_{j})};$
(24)

(ii) update $\bm{\pi}$ by averaging the weights (this is the standard EM
update for estimating mixture proportions, and is the same for all the
algorithms); (iii) update the covariance matrices $\mathcal{U}$ (this is the
step where the algorithms differ); and (iv) update the scaling parameters,
$\bm{s}$. (The update of the scaling parameters is the same for all
algorithms, and depends on the chosen penalty. For details, see the
Supplementary Materials; Yang, Carbonetto and Stephens 2024b.)

With this data augmentation, the updates for the covariance matrices
$\bm{U}_{k}$ have the following form:

$\bm{U}^{\mathrm{new}}=\operatorname*{argmax}_{\bm{U}\,\in\,P^{+,k}_{R}}\phi(\bm{U};\bm{w}_{k}),$
(25)

where

$\phi(\bm{U};\bm{w})\colonequals\sum_{j=1}^{n}w_{j}\log
N_{R}(\bm{x}_{j};\bm{0},\bm{U}+\bm{V}_{j}).$ (26)

The function $\phi$ can be viewed as a weighted version of the log-likelihood
with normal prior (14), and the updates for this weighted problem are very
similar to the updates for a normal prior. For example, the TED update, which
solves the weighted problem exactly in the case $\bm{V}_{j}=\bm{I}_{R}$,
involves truncating the eigenvalues of $\hat{\bm{S}}-\bm{I}_{R}$ where
$\hat{\bm{S}}\colonequals\sum_{j=1}^{n}w_{j}\bm{x}_{j}\bm{x}_{j}^{T}/(\sum_{j=1}^{n}w_{j})$
is the weighted sample covariance matrix. Details of the TED, ED and FA
updates for weighted log-likelihoods are given in the Supplementary Text
(Yang, Carbonetto and Stephens, 2024b).

###  4.3 Modifications to the algorithms to incorporate the penalties

With a penalty, the updates (25) are instead

$\bm{U}^{\mathrm{new}}=\operatorname*{argmax}_{\bm{U}\,\in\,P^{+,k}_{R},\,s_{k}\,>\,0}\phi(\bm{U};\bm{w}_{k})-\rho(\bm{U}/s_{k}).$
(27)

We have adapted the TED approach, in the case $\bm{V}_{j}=\bm{I}_{R}$, to
incorporate either the IW or NN penalty. These extensions replace the simple
truncation of the eigenvalues with solving a (1-d) optimization problem for
each eigenvalue $e_{r}$; while these 1-d optimization problems do not have
closed form solutions, they are easily solved using off-the-shelf numerical
methods. This approach can also be applied, via the data transformation
approach, to the general homoskedastic case $\bm{V}_{j}=\bm{V}$; however, this
transformation approach implicitly changes the penalty so that it encourages
$\bm{U}_{k}/s_{k}$ to be close to $\bm{V}$ rather than being close to
$\bm{I}_{R}$. This change in the penalty appears to be necessary to make the
TED approach work when $\bm{V}\neq\bm{I}_{R}$. See the Supplementary Materials
(Yang, Carbonetto and Stephens, 2024b) for details.

Incorporating an IW penalty into the ED approach is also straightforward
(Bovy, Hogg and Roweis, 2011) and results in a simple change to the closed-
form updates (17). The NN penalty does not result in closed-form updates for
ED and so we have not implemented it.

Incorporating penalties into the FA updates may be possible but we have not
done so.

###  4.4 Significance testing

In EBMNM, inferences are based on the posterior distribution,
$p_{\mathsf{post}}(\bm{\theta}_{j})=p(\bm{\theta}_{j}\mid\bm{x}_{j},\hat{\mathcal{U}},\hat{\bm{\pi}})$,
in which $\hat{\mathcal{U}},\hat{\bm{\pi}}$ denote the estimates returned by
Algorithm 1. To test for significance, we use the local false sign rate
(lfsr), which has been used in both univariate (Stephens, 2017; Xie and
Stephens, 2022) and multivariate (Urbut et al., 2019; Liu et al., 2023; Zou et
al., 2024) settings. The lfsr is defined as

$\mbox{\em lfsr}_{jr}\colonequals\min\\{p_{\mathsf{post}}(\theta_{jr}\geq
0),p_{\mathsf{post}}(\theta_{jr}\leq 0)\\}.$ (28)

In particular, a small lfsr indicates high confidence in the sign of
$\theta_{jr}$. The lfsr is robust to modeling assumptions, which is helpful
for reducing sensitivity to the choice of prior (Stephens, 2017).

##  5 Numerical comparisons

We ran simulations to (i) compare the performance of the different approaches
to updating the covariance matrices ${\bm{U}}_{k}$; (ii) assess the benefits
of the penalties and constraints; and (iii) assess the sensitivity of the
results to the choice of $K$, the number of mixture components in the prior.

We used the Dynamic Statistical Comparisons software
(https://github.com/stephenslab/dsc) to perform the simulations. Code
implementing the simulations is provided in the supplementary ZIP file (Yang,
Carbonetto and Stephens, 2024a), and includes a workflowr website (Blischak,
Carbonetto and Stephens, 2019) for browsing the results. The code and
workflowr website is also available online at
https://github.com/yunqiyang0215/udr-paper.

###  5.1 Data generation

We simulated all data sets from the EBMNM model (1–2); that is, for each data
set, we simulated the “true” means
$\bm{\theta}_{1},\ldots,\bm{\theta}_{n}\in\mathbb{R}^{R}$ from the mixture
prior (2) then we simulated observed data, the “noisy” means
$\bm{x}_{1},\ldots,\bm{x}_{n}\in\mathbb{R}^{R}$, independently given the true
means. Note that model fitting and inferences were performed using only
$\bm{x}_{1},\ldots,\bm{x}_{n}$; the true means
$\bm{\theta}_{1},\ldots,\bm{\theta}_{n}$ were only used to evaluate accuracy
of the inferences. Further, to evaluate the ability of the model to generalize
to other data, we also simulated test sets with true means
$\bm{\theta}_{1}^{\mathrm{test}},\ldots,\bm{\theta}_{n_{\mathrm{test}}}^{\mathrm{test}}$
and observed vectors
$\bm{x}_{1}^{\mathrm{test}},\ldots,\bm{x}_{n_{\mathrm{test}}}^{\mathrm{test}}$.
These test sets were simulated in the same way as the training sets.

In all cases, we set the number of mixture components $K$ to 10, with uniform
mixture weights, $\pi_{1},\ldots,\pi_{10}=1/10$. We generated the $K=10$
covariance matrices $\bm{U}_{1},\ldots,\bm{U}_{10}$ in two different ways,
which we refer to as the “hybrid” and “rank-1” scenarios:

  1. 1.

Hybrid scenario. We used 3 canonical covariance matrices, and randomly
generated an additional 7 covariance matrices randomly from an inverse-Wishart
distribution with scale matrix $\bm{S}=5\bm{I}_{R}$ and $\nu=R+2$ degrees of
freedom. The 3 canonical matrices were as follows:
$\bm{U}_{1}=5\bm{e}_{1}\bm{e}_{1}^{T}$, a matrix of all zeros except for a 5
in the top-left position, which generates “singleton” mean vectors with a
single non-zero element, $\bm{\theta}_{j}=(\theta_{j1},0,\ldots,0)$;
$\bm{U}_{2}=5\bm{1}\bm{1}^{T}$, which generates equal means
$\bm{\theta}_{j}=(\alpha_{j},\ldots,\alpha_{j})$ for some scalar $\alpha_{j}$;
and $\bm{U}_{3}=5\bm{I}_{R}$, which generates mean vectors $\bm{\theta}_{j}$
that are independent in each dimension.

  2. 2.

Rank-1 scenario. We used 5 covariance matrices of the form
$\bm{U}_{k}=5\bm{e}_{k}\bm{e}_{k}^{T}$, $k=1,\ldots 5$, which generate mean
vectors of length $R$ with zeros everywhere except at the $k$th position. The
remaining 5 covariances were random rank-1 matrices of the form
$\bm{U}_{k}=\bm{u}_{k}\bm{u}_{k}^{T}$, $\bm{u}_{k}\sim
N_{R}(\bm{0},\bm{I}_{R})$, $k=6,\dots,10$.

In both scenarios, we simulated large, low-dimension data sets
($n=\mbox{10,000}$, $R=5$), and smaller, high-dimension data sets
($n=\mbox{1,000}$, $R=50$). We refer to these data sets as “large $n/R$” and
“small $n/R$”, respectively. To allow comparisons between TED, ED and FA, in
all cases we simulated data sets with homoskedastic errors; that is,
$\bm{V}_{j}=\bm{I}_{R}$, $j=1,\ldots,n$ and
$\bm{V}_{j}^{\mathrm{test}}=\bm{I}_{R}$, $j=1,\ldots,n_{\mathrm{test}}$.

###  5.2 Results

####  5.2.1 Comparison of convergence

We first focused on comparing the convergence of the different updates (TED,
ED and FA). For brevity, we write “TED” as shorthand for “the algorithm with
TED updates”, and similarly for ED and FA. To compare the updates under the
same conditions, we used FA here to fit full-rank covariance matrices, not
rank-1 matrices. We ran TED and ED with and without the IW penalty.

Figure 1: Illustrative examples comparing convergence of TED, ED and FA. Each
plot shows the algorithms’ progress over iterations on a single simulated data
set. A and C show the results for the same data set, and similarly for B and
D. In all cases, we ran 2,000 iterations, although in some cases, TED updates
stopped early because the updates converged to a stationary point of the
objective (the likelihood or the penalized likelihood). For examples A and B,
no penalty was used; for examples C and D, the inverse Wishart (IW) penalty
was used with penalty strength $\lambda=R$. Log-likelihood and penalized log-
likelihood differences are plotted with respect to the (penalized) log-
likelihood near the initial estimate. An initial round of 20 ED iterations
common to all the runs (the “warm start”) is not shown in these plots.

We ran all methods on 100 “large $n/R$” data sets and 100 “small $n/R$” data
sets simulated in the hybrid scenario, setting $K=10$ to match the simulated
truth. To reduce the likelihood that the different updates converge to
different local solutions, we performed a prefitting stage in which we ran 20
iterations of ED from a random initial starting point (specifically, the
initialization was $\pi_{k}=1/10$, $s_{k}=1$, with randomly generated
$\bm{U}_{k}$, $k=1,\ldots,10$). We call this a “warm start”. We then ran each
algorithm for at most 2,000 iterations after this warm start. Figure 1 shows
illustrative results for two data sets (one large $n/R$ and one small $n/R$),
and Figure 2 summarizes the results across all simulations.

The results in Figure 1 illustrate typical behavior. Among the unpenalized
updates, TED and FA converged much more quickly than ED. For the penalized
updates, TED still converged more quickly than ED, but the difference is less
striking than without the penalty. In the small $n/R$ example the three
methods appear to have converged to different solutions (despite the use of a
warm start). This is not unexpected due to the non-convexity of the objective
function, and illustrates an important general point to keep in mind:
improvements in the quality of solution obtained may be due to either faster
convergence to the solution, convergence to a local solution with higher
objective, or both.

Figure 2: Comparison of convergence of TED, ED and FA. Each plot summarizes
the results from 100 simulations. In each simulation, log-likelihood achieved
after (at most) 2,000 iterations was recorded. Panels E and F show differences
in the penalized log-likelihoods (IW penalty, $\lambda=R$).

The results in Figure 2 confirm that some of the patterns observed in the
illustrative example are true more generally. In the unpenalized case (Panels
A–D), TED often arrived at a better solution than FA or ED, although in the
large $n/R$ setting FA was comparable to TED (Panel A). In the penalized case,
the TED and ED solutions were much more similar for both large and small $n/R$
(Panels E, F).

Figure 3: Results demonstrating the ability of TED to “rescue” ED. The log-
likelihood or penalized log-likelihood obtained by performing at most 2,000
TED updates is compared against performing 2,000 iterations of ED, followed by
another round of (at most) 2,000 iterations of TED updates (“ED+TED”). Each
plot summarizes the results from 100 simulations, the same simulations as in
Figure 2.

The improved performance of (unpenalized) TED vs. ED could be due either to
faster convergence (e.g., Figure 1A) or due to convergence to a better local
optimum (e.g., Figure 1B). To investigate this, we assessed whether TED can
“rescue” ED by running TED initialized to the ED solution (“ED+TED”). If ED
converges to a poorer local optimum, then TED will not rescue it, and ED+TED
will be similar to ED; on the other hand, if ED is simply slow to converge,
then ED+TED will be similar to TED. The results in Figure 3 show that TED
usually rescues ED; the ED+TED estimates were consistently very similar to the
TED estimates, regardless of whether a penalty was used or not. This suggests
that the improved performance of TED is generally due to faster convergence.

To assess how additional iterations improve ED’s performance, we reran ED for
100,000 iterations instead of only 2,000. Since these runs take a long time,
we examined only a few simulated data sets randomly selected from the large
$n/R$ and small $n/R$ scenarios (Supplementary Figure
LABEL:supple:fig1_convergence). Even after 100,000 iterations, ED fell
measurably short of 1,000 updates of TED (average difference of 40.6 log-
likelihood units).

In summary, our experiments confirm that, for homoskedastic errors (which is
the setting where all three methods apply), TED exhibited the best
performance. This is not unexpected since TED solves the subproblem (eq. 25 or
27) exactly whereas ED and FA do not. Our experiments also show that including
a penalty can help improve convergence behavior, especially for ED. Thus,
including the penalty has computational benefits in addition to statistical
benefits demonstrated below.

####  5.2.2 Comparison of the penalties and constraints

Next we evaluated the benefits of different penalties and constraints for
estimation and significance testing of the underlying means $\bm{\theta}_{j}$.
We focused on TED and ED with and without penalties, and on fitting rank-1
matrices using TED. (Since TED solves the subproblem exactly, rather than
iteratively, TED should be better than FA in this setting with homoskedastic
variances. The main benefit of FA is that it can fit rank-1 matrices with
heteroskedastic variances.)

We compared methods using the following evaluation measures:

  * •

Plots of power vs. false sign rate (FSR). These plots are similar to the more
commonly-used plots of “power vs. false discovery rate,” but improve
robustness and generality by requiring “true discoveries” to have the correct
sign. The better methods are those that achieve higher power at a given FSR.
See the Supplementary Text (Yang, Carbonetto and Stephens, 2024b) for
definitions.

  * •

Empirical False Sign Rate (FSR). We report the empirical FSR among tests that
were significant at $\text{\em lfsr}<0.05$. A well-behaved method should have
a small empirical FSR, certainly smaller than 0.05. We consider an FSR
exceeding 0.05 to be indicative of a poorly behaved method.

  * •

Accuracy of predictive distribution. To assess generalizability of the
estimates of $\mathcal{U}$ and $\bm{\pi}$, we compared the marginal predictive
density (3) in test samples,
$p(\bm{x}_{j}^{\mathrm{test}}\mid\hat{\mathcal{U}},\hat{\bm{\pi}},\bm{V}_{j}^{\mathrm{test}})$,
against the “ground-truth” marginal predictive density
$p(\bm{x}_{j}^{\mathrm{test}}\mid\mathcal{U}^{\mathrm{true}},{\bm{\pi}}^{\mathrm{true}},\bm{V}_{j}^{\mathrm{test}})$,
where $\mathcal{U}^{\mathrm{true}},{\bm{\pi}}^{\mathrm{true}}$ denote the
parameters used to simulate the data. We summarized the relative accuracy in
the predictions as

$\frac{1}{n_{\mathrm{test}}}\sum_{j=1}^{n_{\mathrm{test}}}\log\bigg{\\{}\frac{p(\bm{x}_{j}^{\mathrm{test}}\mid\mathcal{U}^{\mathrm{true}},\bm{\pi}^{\mathrm{true}})}{p(\bm{x}_{j}^{\mathrm{test}}\mid\hat{\mathcal{U}},\hat{\bm{\pi}})}\bigg{\\}}.$
(29)

This measure can be interpreted as (an approximation of) the Kullback-Leibler
(K-L) divergence from the true predictive distribution to the estimated
predictive distribution. Smaller K-L divergences are better.

Figure 4: Comparison of penalties, constraints and updates (ED vs. TED) in the
“hybrid” simulated data sets. For the IW and NN penalties, the penalty
strength was set to $\lambda=R$. In C and D, the target FSR is shown as a
dotted horizontal line at 0.05. In A, most of the methods are not visible
because the lines overlap at the top near the oracle result. Note that the
oracle model always achieves a K-L divergence of zero. Figure 5: Comparison of
penalties, constraints and updates (ED vs. TED) in the “rank-1” simulated data
sets. For the IW and NN penalties, the penalty strength was set to
$\lambda=R$. In C and D, the target FSR is shown as a dotted horizontal line
at 0.05. Note that the oracle model always achieves a K-L divergence of zero.

We simulated 20 data sets with large $n/R$ and another 20 data sets with small
$n/R$ under both the hybrid and rank-1 scenarios (80 data sets in total). In
all cases, model fitting was performed as above, again with $K=10$. The
results are summarized in Figures 4 and 5. In all comparisons, we included
results from the “oracle” EBMNM model—that is, the model used to simulate the
data—as a point of reference.

Results for the hybrid setting are shown in Figure 4. The results show a clear
benefit of using a penalty in the small $n/R$ setting: both IW and NN
penalties improved the power vs. FSR and the accuracy of predictive
distributions. For large $n/R$ data sets, the penalties do not provide a clear
benefit, but also do not hurt performance. In both cases TED and ED perform
similarly, suggesting that the poor convergence of ED observed in previous
comparisons may have less impact on performance than might have been expected.
The rank-1 constraints performed very poorly in all tasks, which is perhaps
unsurprising since the true covariances were (mostly) not rank-1.

Results for the rank-1 scenario are shown in Figure 5. In this case, imposing
rank-1 constraints on the covariance matrices improved predictive
performance—which makes sense because the true covariances were indeed
rank-1—but produced worse performance in other metrics. In particular, the
lfsr values from the rank-1 constraint are very poorly calibrated. This is
because, as noted in Liu et al. (2023), the rank-1 constraint leads to lfsr
values that do not differ across conditions; see the Discussion (Section 7)
for more on this. Penalized estimation of the covariance matrices (using
either a IW or NN penalty) consistently achieved the best power at a given FSR
in both the large $n/R$ and small $n/R$ settings. Interestingly, (unpenalized)
TED performed much worse than (unpenalized) ED in the power vs. FSR for large
$n/R$. We speculate that this was due to the slow convergence of ED providing
sort of “implicit regularization”. However, explicit regularization via a
penalty seems preferable to implicit regularization via a poorly converging
algorithm, and overall penalized estimation was the winning (or equally
winning) strategy across a variety of settings.

####  5.2.3 Robustness to mis-specifying the number of mixture components

In the above experiments, we fit all models with a value of $K$ that matched
the model used to simulate the data. In practice, however, $K$ is unknown, and
so we must also consider situations in which $K$ is mis-specified.
Intuitively, one might expect that overstating $K$ may lead to overfitting and
worse performance; Urbut et al. (2019) argued however that the use of the
mixture components centered at zero in the prior (2) makes it robust to
overfitting because the mean-zero constraint limits their flexibility. Here we
assess this claim.

Figure 6: Assessment of robustness to mis-specifying $K$ in “hybrid” simulated
data sets. All results are averages over 20 data sets, each simulated with
$K=10$ mixture components. In A and C, most of the methods are not visible
because they overlap with the “TED-IW” result. Figure 7: Assessment of
robustness to mis-specifying $K$ in “rank-1” simulated data sets. All the
results shown in the plots are averages over the 20 data sets. All data sets
were simulated with $K=10$ mixture components.

In these experiments, we analyzed the same 80 data sets as in the previous
section (simulated with $K=10$). We fit models with different penalties,
constraints and algorithms, with $K$ varying from 2 to 128. We compared
results in both the accuracy of the predictive distribution (K-L divergence)
and the average FSR at an lfsr threshold of 0.05.

Results are shown in Figures 6 and 7. For large $n/R$, all model fits except
those with the rank-1 constraints were robust to overstating $K$, with
similarly good performance even at $K=128$. This is generally consistent with
the claim in Urbut et al. (2019) that the mixture prior should be robust to
overstating $K$. However, for small $n/R$ the story is quite different: all of
the unpenalized algorithms eventually showed a decline in performance when $K$
was too large, presumably due to “overfitting”. In comparison, all the
penalized methods were more robust to overstating $K$; the performance did not
substantially decline as $K$ increased.

The improved robustness of the penalized methods could be achieved in at least
two different ways: they could be using a smaller number of components by
estimating some of the mixture weights $\pi_{k}$ to be very small; or by
estimating some of the components $k$ to have very similar covariances
$\bm{U}_{k}$ (or both). To investigate these explanations, we looked at the
models with $K=100$ components and recorded the number of “important”
components, defined as components $k$ with $\pi_{k}>0.01$. We found that the
penalized methods tended to produce much fewer “important” components than the
unpenalized methods (Supplementary Figure LABEL:supple:compare_weights).
Essentially, the penalties have the effect of shrinking each $\bm{U}_{k}$
toward the identity matrix, so a component is assigned a small weight whenever
the identity matrix does not match the truth.

In summary, our results support the use of penalized methods with a large
value of $K$ as a simple and robust way to achieve good performance in
different settings.

##  6 Analysis of genetic effects on gene expression in 49 human tissues

To illustrate our methods on real data, we used the EBMNM model to analyze the
effects of genetic variants on gene expression (“cis-eQTLs”) in multiple
tissues. The use of the EBMNM model for this purpose was first demonstrated in
Urbut et al. (2019). They used a two-stage procedure to fit the EBMNM models:
(i) fit the EBMNM model to a subset of “strong” eQTLs to estimate the prior
covariance matrices; and (ii) fit a (modified) EBMNM model to all eQTLs—or a
random subset of eQTLs—using the covariance matrices from (i). The model from
(ii) was then used to perform inferences and to test for eQTLs in each tissue.
Our new methods are relevant to (i), and so we focus on stage (i) here. For
(i), Urbut et al. (2019) used the ED algorithm without a penalty. Recognizing
that the results are sensitive to initialization, they described a detailed
initialization procedure.111An updated version of the initialization procedure
of Urbut et al. (2019) is described in the “flash_mash” vignette included in
the mashr R package. See also
https://stephenslab.github.io/mashr/articles/flash_mash.html. We coded an
initialization procedure similar to this, which we refer to as the
“specialized initialization”. We note that the specialized initialization adds
substantially to both the complexity and computation time of the overall
fitting procedure. Our goal was to assess the benefits of different EBMNM
analysis pipelines which consisted of combinations of: TED vs. ED updates;
specialized initialization vs. a simple random initialization; and penalty vs.
no penalty (maximum-likelihood). All these combinations resulted in 8
different analyses of multi-tissue cis-eQTL data.

We analyzed z-scores from tests for association between gene expression in
dozen human tissues and genotypes at thousands of genetic variants. The
z-scores came from running the “Matrix eQTL” software (Shabalin, 2012) on
genotype and gene expression data in release 8 of the Genotype-Tissue
Expression (GTEx) Project (GTEx Consortium, 2020). Following Urbut et al.
(2019), we selected the genetic variants with the largest z-score (in
magnitude) across tissues for each gene. After data filtering steps, we ended
up with a data set of z-scores for $n=\mbox{15,636}$ genes and $R=49$ tissues.
We fit the EBMNM model to these data, with all the $\bm{V}_{j}$ set to a
common correlation matrix; that is, $\bm{V}_{j}={\bm{C}}$, where ${\bm{C}}$ is
a correlation matrix of non-genetic effects on expression. This correlation
matrix ${\bm{C}}$ was estimated from the association test z-scores following
the approach described in Urbut et al. (2019). In all runs, we set $K=40$ to
match the number of covariance matrices produced by our specialized
initialization. And, in all cases, we initialized the mixture weights to
$\pi_{k}=1/40$ and the scaling factors (when needed) to $s_{k}=1$. The
supplementary ZIP file (Yang, Carbonetto and Stephens, 2024a) contains the
code and data implementing these analyses as well as the results we generated.

Figure 8: Plots showing improvement in model fit over time for the GTEx data,
using different initialization schemes, different prior covariance matrix
updates, and penalty vs. no penalty (maximum-likelihood). Log-likelihood
differences and penalized log-likelihood differences are with respect to the
(penalized) log-likelihood near the initial estimate. All models were fit with
$K=40$ mixture components. In B, the inverse wishart (IW) penalty was used
with penalty parameter $\lambda=R$. In all cases, the model fitting was halted
when the difference in the log-likelihood between two successive updates was
less than 0.01, or when 5,000 updates were performed, whichever came first.

Figure 8 shows the improvement of the model fits over time in the 8 different
analyses. In the analyses without a penalty (A), the TED updates with a
specialized initialization produced the best fit, while the ED and TED updates
with random initialization were somewhat worse (e.g., TED with random
initialization produced a fit that was 2,497.78 log-likelihood units worse, or
0.16 log-likelihood units per gene). Strikingly, the ED updates with
specialized initialization resulted in a much worse fit. We attribute this to
the fact that the specialized initialization includes many rank-1 matrices,
and the “subspace preserving” property of the ED updates means that these
matrices are fixed at their initialization (they changed only by a scaling
factor), which substantially limits their ability to adapt to the data. In the
penalized case, consistent with our simulation results, there was less
difference between ED vs. TED. In both cases, the specialized initialization
improved the fit relative to random initialization (e.g., TED increased the
penalized log-likelihood by 3,176, or about 0.2 per gene). Note that adding a
penalty function makes the subspace preserving property of the ED updates
irrelevant by forcing the matrices to be full rank.

To compare the quality of the fits obtained by each method, we computed log-
likelihoods on held-out (“test set”) data, using a 5-fold cross-validation
(CV) design. Following the usual CV setup, in each CV fold 80% of the genes
were in the training set, and the remaining 20% were in the test set. Then we
fit an EBMNM model to the training set following each of the 8 approaches
described above, and measured the quality of the fit by computing the log-
likelihood in the test set. We also recorded the number of iterations. Within
a given fold, the number of components $K$ was the same across all the
analyses, and was set depending on the number of covariance matrices produced
by the specialized initialization. ($K$ was at least $32$ and at most $39$.)

Consistent with the simulations, the inclusion of a penalty consistently
improved the test set log-likelihood (Table 2). With a penalty, the
specialized initialization also improved the test set log-likelihood compared
with a random initialization. Of the 8 approaches tried, the one used by Urbut
et al. (2019)—ED with no penalty and specialized initialization—resulted in
the worst test-set likelihood.

Table 2: Cross-validation results on the GTEx data. The “mean relative log-likelihood” column gives the increase in the test-set log-likelihood over the worst log-likelihood among the 8 approaches compared, divided by total number of genes in each test set. The “average number of iterations” column gives the number of iterations performed until the stopping criterion is met (log-likelihood between two successive updates less than 0.01, up to a maximum of 5,000 iterations), averaged over the 5 CV folds. |  |  | mean relative | average number
---|---|---|---|---
initialization | algorithm | penalty | log-likelihood | of iterations
specialized | ED | none | 0.00 | 1,101
specialized | ED | IW | 1.21 | 1,083
specialized | TED | none | 0.88 | 1,054
specialized | TED | IW | 1.19 | 412
random | ED | none | 0.25 | 5,000
random | ED | IW | 0.86 | 1,377
random | TED | none | 0.20 | 450
random | TED | IW | 0.94 | 584
Figure 9: Comparison of the prior mixture weights $\bm{\pi}$ from the previous
pipeline vs. the new pipeline. The histogram shows the distributions of the
$K40$ prior mixture weights resulting from both pipelines. Figure 10: Top
effect-sharing patterns in the EBMNM model fit to the GTEx data using the
previous pipeline. The “top” patterns are the mixture components with the
largest weights. Each heatmap shows the 49 $\times$ 49 the scaled covariance
matrix $\bm{U}_{k}/\sigma_{k}^{2}$, where $\sigma_{k}^{2}$ is the largest
diagonal element of $\bm{U}_{k}$, so that all elements of the scaled
covariance lie between $-1$ and 1. (The vast majority of the covariances are
positive; negative correlations are unexpected.) The scaled covariance
matrices are arranged in decreasing order by mixture weight. The “scale” above
each heatmap is $\sigma_{k}$. Note that the top three covariance matrices
capture broadly similar effect-sharing patterns, but different effect scales.
Figure 11: Top effect-sharing patterns in the GTEx data generated by the new
analysis pipeline. See the caption to Fig. 10 for more details.

Although the analyses with a specialized initialization resulted in better
fits than the analyses with a random initialization, the specialized
initialization also had substantial computational overhead; running the
initialization procedures on these data took more than 1 hour (by comparison,
running a single iteration of the EBMNM algorithm typically took about 1
second). Therefore, on balance, one might prefer to dispense with the
specialized initialization. Based on these results and considerations, we
subsequently examined in more detail the analysis with ED, no penalty and
specialized initialization—which was the approach used in Urbut et al.
(2019)—and the analysis with TED, IW penalty and random initialization. For
brevity, we refer to these analyses as the “previous pipeline” and “new
pipeline”, respectively.

A notable outcome of the new pipeline is that it produced a prior with weights
that were more evenly distributed across the mixture components (Figure 9).
For example, the new pipeline produced 22 covariances with weights greater
than 1%, whereas the previous pipeline produced only 10 covariances with
weights more than 1%. Also, the top 16 covariances accounted for 85% of the
total weight in the new pipeline, whereas only 6 covariances were needed to
equal 85% total weight in the previous pipeline. Inspecting the individual
covariances generated from the two analysis pipelines, there many strong
similarities in the estimated sharing patterns (Figures 10 and 11). But the
new analysis pipeline learned a greater variety of tissue-specific patterns
(e.g., whole blood, testis, thyroid) and tissue-sharing patterns, many of
which appear to reflect underlying tissue biology. For example, sharing
pattern 6 (see Figure 11) captures sharing of brain-specific effects
(including the pituitary gland, which is found at the base of the brain near
the hypothalamus). Sharing pattern 8 may reflect the fact that skin and the
mucosa layer of the esophagus wall both contain squamous epithelial cells.
Additionally, we looked at more closely at the 8 genes with very strong
posterior weights (¿98%) on component 6, the component capturing eQTL sharing
among brain tissues: HOMER1, JPT1, SRSF2, ABCA1, GPATCH8, SYNGAP1, SPI1 and
LGMN. Several of these have been linked to neurological and neuropsychiatric
conditions, including major depression, schizophrenia, attention dementia and
Alzheimer’s disease. For example, SYNGAP1 is linked to neuronal functions and
psychiatric diseases based on results from the GWAS Catalog (Sollis et al.,
2022) and from functional studies (Jeyabalan and Clement, 2016; Llamosas et
al., 2020). In summary, the improvements to the EBMNM analysis pipeline should
result in the discovery of a greater variety of cross-tissue genetic effects
on gene expression.

##  7 Discussion

The growing interest in deciphering shared underlying biological mechanisms
has led to a surge in multivariate analyses in genomics; among the many recent
examples are multi-trait analyses (Wu et al., 2020; Luo et al., 2020) and
multi-ancestry analyses to improve polygenic risk scores (Zhang et al., 2023).
The EBMNM approach described here and in Urbut et al. (2019) provide a
versatile and robust multivariate approach to multivariate analysis. We have
implemented and compared several algorithms for this problem. These algorithms
not only enhance accuracy but also provide a level of flexibility not achieved
by other methods.

One of our important findings is that using low-rank covariance matrices in
this setting, as was done in Urbut et al. (2019), is not recommended. In
particular, while low-rank matrices may be relatively easy to interpret, they
lead to poorly calibrated lfsr values (e.g., Figure 5). For intuition,
consider fitting a EBMNM model with $K=1$ and a rank-1 covariance matrix,
$\bm{U}_{1}=\bm{u}\bm{u}^{T}$. (We credit Dongyue Xie for this example.) Under
this model, the mean is $\bm{\theta}_{j}=\bm{u}a_{j}$ for some $a_{j}$, and
therefore, given a $\bm{u}$, the signs of all the elements of
$\bm{\theta}_{j}$ are fully determined by $a_{j}$. As a result, all elements
of $\bm{\theta}_{j}$ will have the same lfsr, and so this model cannot capture
situations where one is confident in the sign of some elements of
$\bm{\theta}_{j}$ but not others. This can cause problems even if the model is
correct, that is, when the true covariances are rank-1, as in Figure 5. There
are some possible ways to address these issues, say, by imposing sparsity on
estimates of $\bm{u}$, and this could be an area for future work.

Even when the pitfalls of low-rank matrices are avoided, it is still the case
that EB methods tend to understate uncertainty compared to “fully Bayesian”
methods (e.g., Morris 1983; Wang and Titterington 2005). As a result, one
should expect that the estimated lfsr values may be anti-conservative; that
is, the lfsr values are smaller than they should be. Indeed, we saw this anti-
conservative behavior across many of our simulations and methods
(Supplementary Figures LABEL:supple:calibration and
LABEL:supple:calibration_rank1). For this reason estimated lfsr rates should
be treated with caution, and it would be prudent to use more stringent
significance thresholds than are actually desired (e.g., an lfsr threshold of
0.01 rather than 0.05). In the special case where $\bm{V}=\bm{I}$, it should
be possible to improve calibration of significance tests by using ideas from
Lei and Fithian (2018). Improving calibration in the general case of dependent
multivariate tests seems to be an important area for future research.

[Acknowledgments] We thank Gao Wang and Sarah Kim-Hellmuth for their help in
preparing the GTEx data. We thank Yuxin Zou, Yusha Liu and Dongyue Xie for
helpful discussions on MASH and the mashr software. And we thank the staff at
the Research Computing Center at the University of Chicago for maintaining the
computing resources used in our numerical experiments.

This work was supported by NIH grant R01HG002585 to MS.

Supplementary materials

The supplementary PDF includes supporting mathematical results and additional
tables and figures.

R package and code reproducing the analyses

ZIP file containing the source code for the R package as well as the source
code reproducing the results of the simulations and analysis of the GTEx data.

## References

  * Araujo et al. (2023) [author] Araujo, Daniel SD. S., Nguyen, ChrisC., Hu, XiaoweiX., Mikhaylova, Anna VA. V., Gignoux, ChrisC., Ardlie, KristinK. et al. (2023). Multivariate adaptive shrinkage improves cross-population transcriptome prediction and association studies in underrepresented populations. Human Genetics and Genomics Advances 4 100216\.
  * Barbeira et al. (2020) [author] Barbeira, Alvaro NA. N., Melia, Owen JO. J., Liang, YanyuY., Bonazzola, RodrigoR., Wang, GaoG., Wheeler, Heather EH. E., Aguet, FrançoisF., Ardlie, Kristin GK. G., Wen, XiaoquanX. and Im, Hae KH. K. (2020). Fine-mapping and QTL tissue-sharing information improves the reliability of causal gene identification. Genetic Epidemiology 44 854–867.
  * Bhadra et al. (2019) [author] Bhadra, AnindyaA., Datta, JyotishkaJ., Polson, Nicholas G.N. G. and Willard, BrandonB. (2019). Lasso meets horseshoe: a survey. Statistical Science 34 405–427.
  * Blischak, Carbonetto and Stephens (2019) [author] Blischak, John DJ. D., Carbonetto, PeterP. and Stephens, MatthewM. (2019). Creating and sharing reproducible research code the workflowr way [version 1; peer review: 3 approved]. F1000Research 8\.
  * Bonder et al. (2021) [author] Bonder, Marc JanM. J., Smail, CraigC., Gloudemans, Michael JM. J., Frésard, LaureL., Jakubosky, DavidD., D’Antonio, MatteoM. et al. (2021). Identification of rare and common regulatory variants in pluripotent cells using population-scale transcriptomics. Nature Genetics 53 313–321.
  * Bovy, Hogg and Roweis (2011) [author] Bovy, JoJ., Hogg, David W.D. W. and Roweis, Sam T.S. T. (2011). Extreme Deconvolution: Inferring complete distribution functions from noisy, heterogeneous and incomplete observations. Annals of Applied Statistics 5 1657–1677.
  * Cai and Liu (2011) [author] Cai, TonyT. and Liu, WeidongW. (2011). Adaptive thresholding for sparse covariance matrix estimation. Journal of the American Statistical Association 106 672–684.
  * Chi and Lange (2014) [author] Chi, Eric CE. C. and Lange, KennethK. (2014). Stable estimation of a covariance matrix guided by nuclear norm penalties. Computational statistics and Data Analysis 80 117–128.
  * GTEx Consortium (2020) [author] GTEx Consortium (2020). The GTEx Consortium atlas of genetic regulatory effects across human tissues. Science 369 1318–1330.
  * GTEx Consortium et al. (2015) [author] GTEx Consortium, Ardlie, Kristin GK. G., Deluca, David SD. S., Segrè, Ayellet VA. V., Sullivan, Timothy JT. J., Young, Taylor RT. R., Gelfand, Ellen TE. T., Trowbridge, Casandra AC. A., Maller, Julian BJ. B., Tukiainen, TaruT. et al. (2015). The Genotype-Tissue Expression (GTEx) pilot analysis: multitissue gene regulation in humans. Science 348 648–660.
  * Dempster, Laird and Rubin (1977) [author] Dempster, Arthur PA. P., Laird, Nan MN. M. and Rubin, Donald BD. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B 39 1–22.
  * Efron and Morris (1972) [author] Efron, BradleyB. and Morris, CarlC. (1972). Limiting the risk of Bayes and empirical Bayes estimators—Part II: the empirical Bayes case. Journal of the American Statistical Association 67 130–139.
  * Fan, Liao and Liu (2016) [author] Fan, JianqingJ., Liao, YuanY. and Liu, HanH. (2016). An overview of the estimation of large covariance and precision matrices. Econometrics Journal 19 C1–C32.
  * Flutre et al. (2013) [author] Flutre, TimothéeT., Wen, XiaoquanX., Pritchard, JonathanJ. and Stephens, MatthewM. (2013). A statistical framework for joint eQTL analysis in multiple tissues. PLoS Genetics 9 e1003486.
  * Fraley and Raftery (2007) [author] Fraley, ChrisC. and Raftery, Adrian EA. E. (2007). Bayesian regularization for normal mixture estimation and model-based clustering. Journal of Classification 24 155–181.
  * Friedman, Hastie and Tibshirani (2008) [author] Friedman, JeromeJ., Hastie, TrevorT. and Tibshirani, RobertR. (2008). Sparse inverse covariance estimation with the graphical lasso. Biostatistics 9 432–441.
  * Ghahramani and Hinton (1996) [author] Ghahramani, ZoubinZ. and Hinton, Geoffrey EG. E. (1996). The EM algorithm for mixtures of factor analyzers. Technical Report, University of Toronto.
  * Han and Eskin (2011) [author] Han, BuhmB. and Eskin, EleazarE. (2011). Random-effects model aimed at discovering associations in meta-analysis of genome-wide association studies. American Journal of Human Genetics 88 586–598.
  * Jeyabalan and Clement (2016) [author] Jeyabalan, NallathambiN. and Clement, James PJ. P. (2016). SYNGAP1: mind the gap. Frontiers in Cellular Neuroscience 10 32\.
  * Johnstone (2019) [author] Johnstone, IainI. (2019). Gaussian estimation: sequence and wavelet models. Available at https://imjohnstone.su.domains/GE_08_09_17.pdf.
  * Johnstone and Paul (2018) [author] Johnstone, Iain M.I. M. and Paul, DebashisD. (2018). PCA in high dimensions: an orientation. Proceedings of the IEEE 106 1277–1292.
  * Kim et al. (2020) [author] Kim, YoungseokY., Carbonetto, PeterP., Stephens, MatthewM. and Anitescu, MihaiM. (2020). A fast algorithm for maximum likelihood estimation of mixture proportions using sequential quadratic programming. Journal of Computational and Graphical Statistics 29 261–273.
  * Ledoit and Wolf (2004) [author] Ledoit, OlivierO. and Wolf, MichaelM. (2004). A well-conditioned estimator for large-dimensional covariance matrices. Journal of Multivariate Analysis 88 365–411.
  * Ledoit and Wolf (2022) [author] Ledoit, OlivierO. and Wolf, MichaelM. (2022). The power of (non-)linear shrinking: a review and guide to covariance matrix estimation. Journal of Financial Econometrics 20 187–218.
  * Lei and Fithian (2018) [author] Lei, LihuaL. and Fithian, WilliamW. (2018). AdaPT: an interactive procedure for multiple testing with side information. Journal of the Royal Statistical Society, Series B 80 649–679.
  * Li et al. (2022) [author] Li, QinQ., Gloudemans, Michael JM. J., Geisinger, Jonathan MJ. M., Fan, BomingB., Aguet, FrançoisF., Sun, TaoT., Ramaswami, GokulG., Li, Yang IY. I., Ma, Jin-BiaoJ.-B., Pritchard, Jonathan KJ. K. et al. (2022). RNA editing underlies genetic risk of common inflammatory diseases. Nature 608 569–577.
  * Lin et al. (2024) [author] Lin, WenheW., Wall, Jeffrey D.J. D., Li, GeG., Newman, DeborahD., Yang, YunqiY., Abney, MarkM., VandeBerg, John L.J. L., Olivier, MichaelM., Gilad, YoavY. and Cox, Laura A.L. A. (2024). Genetic regulatory effects in response to a high-cholesterol, high-fat diet in baboons. Cell Genomics 4 100509\.
  * Liu and Rubin (1998) [author] Liu, ChuanhaiC. and Rubin, Donald BD. B. (1998). Maximum likelihood estimation of factor analysis using the ECME algorithm with complete and incomplete data. Statistica Sinica 8 729–747.
  * Liu et al. (2023) [author] Liu, YushaY., Carbonetto, PeterP., Takahama, MichihiroM., Gruenbaum, AdamA., Xie, DongyueD., Chevrier, NicolasN. and Stephens, MatthewM. (2023). A flexible model for correlated count data, with application to analysis of gene expression differences in multi-condition experiments. arXiv 2210.00697.
  * Llamosas et al. (2020) [author] Llamosas, NereaN., Arora, VineetV., Vij, RidhimaR., Kilinc, MuratM., Bijoch, LukaszL., Rojas, CamiloC. et al. (2020). SYNGAP1 controls the maturation of dendrites, synaptic function, and network activity in developing human neurons. Journal of Neuroscience 40 7980–7994.
  * Luo et al. (2020) [author] Luo, LanL., Shen, JudongJ., Zhang, HongH., Chhibber, AparnaA., Mehrotra, Devan VD. V. and Tang, Zheng-ZhengZ.-Z. (2020). Multi-trait analysis of rare-variant association summary statistics using MTAR. Nature Communications 11 2850\.
  * McLachlan and Peel (2000) [author] McLachlan, GeoffreyG. and Peel, DavidD. (2000). Finite mixture models. John Wiley & Sons, Inc., New York, NY.
  * Morris (1983) [author] Morris, Carl N.C. N. (1983). Parametric empirical Bayes inference: theory and applications. Journal of the American Statistical Association 78 47–55.
  * Natri et al. (2024) [author] Natri, Heini MH. M., Del Azodi, Christina BC. B., Peter, LanceL., Taylor, Chase JC. J., Chugh, SagrikaS., Kendle, RobertR. et al. (2024). Cell-type-specific and disease-associated expression quantitative trait loci in the human lung. Nature Genetics 56 595–604.
  * Pickrell et al. (2016) [author] Pickrell, Joseph KJ. K., Berisa, TomazT., Liu, Jimmy ZJ. Z., Ségurel, LaureL., Tung, Joyce YJ. Y. and Hinds, David AD. A. (2016). Detection and interpretation of shared genetic influences on 42 human traits. Nature Genetics 48 709–717.
  * Robbins (1951) [author] Robbins, HerbertH. (1951). Asymptotically subminimax solutions of compound statistical decision problems. In Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability, 1951, vol. II 131–149. University of California Press, Berkeley and Los Angeles, CA.
  * Rubin and Thayer (1982) [author] Rubin, Donald BD. B. and Thayer, Dorothy TD. T. (1982). EM algorithms for ML factor analysis. Psychometrika 47 69–76.
  * Shabalin (2012) [author] Shabalin, Andrey A.A. A. (2012). Matrix eQTL: ultra fast eQTL analysis via large matrix operations. Bioinformatics 28 1353–1358.
  * Soliai et al. (2021) [author] Soliai, Marcus MM. M., Kato, AtsushiA., Helling, Britney AB. A., Stanhope, Catherine TC. T., Norton, James EJ. E., Naughton, Katherine AK. A. et al. (2021). Multi-omics colocalization with genome-wide association studies reveals a context-specific genetic mechanism at a childhood onset asthma risk locus. Genome Medicine 13 157\.
  * Sollis et al. (2022) [author] Sollis, ElliotE., Mosaku, AbayomiA., Abid, AlaA., Buniello, AnnalisaA., Cerezo, MariaM., Gil, LaurentL. et al. (2022). The NHGRI-EBI GWAS Catalog: knowledgebase and deposition resource. Nucleic Acids Research 51 D977–D985.
  * Stephens (2017) [author] Stephens, MatthewM. (2017). False discovery rates: a new deal. Biostatistics 18 275–294.
  * Sun (2020) [author] Sun, LeiL. (2020). Topics on Empirical Bayes Normal Means, PhD thesis, University of Chicago, Chicago, IL.
  * Tipping and Bishop (1999) [author] Tipping, Michael EM. E. and Bishop, Christopher MC. M. (1999). Probabilistic principal component analysis. Journal of the Royal Statistical Society, Series B 61 611–622.
  * Turchin and Stephens (2019) [author] Turchin, Michael CM. C. and Stephens, MatthewM. (2019). Bayesian multivariate reanalysis of large genetic studies identifies many new associations. PLoS Genetics 15 e1008431.
  * Udler et al. (2018) [author] Udler, Miriam SM. S., Kim, JaegilJ., von Grotthuss, MarcinM., Bonàs-Guarch, SílviaS., Cole, Joanne BJ. B., Chiou, JoshuaJ., Christopher D. Anderson on behalf of METASTROKE and the ISGC, Boehnke, MichaelM. et al. (2018). Type 2 diabetes genetic loci informed by multi-trait associations point to disease mechanisms and subtypes: a soft clustering analysis. PLoS Medicine 15 e1002654.
  * Urbut et al. (2019) [author] Urbut, Sarah MS. M., Wang, GaoG., Carbonetto, PeterP. and Stephens, MatthewM. (2019). Flexible statistical methods for estimating and testing effects in genomic studies with multiple conditions. Nature Genetics 51 187–195.
  * Urbut et al. (2021) [author] Urbut, Sarah MS. M., Selvaraj, Margaret SunithaM. S., Pampana, AkhilA., Dattilo, LillianL., Neale, BenjaminB., ODonnell, Christopher JC. J., Peloso, GinaG. and Natarajan, PradeepP. (2021). Bayesian multivariate genetic analysis of lipids improves translational insight. Circulation 144 A9855–A9855.
  * Wang and Titterington (2005) [author] Wang, BoB. and Titterington, D. M.D. M. (2005). Inadequacy of interval estimates corresponding to variational Bayesian approximations. In Proceedings of the 10th International Workshop on Artificial Intelligence and Statistics 373–380.
  * Wen and Stephens (2014) [author] Wen, XiaoquanX. and Stephens, MatthewM. (2014). Bayesian methods for genetic association analysis with heterogeneous subgroups: from meta-analyses to gene-environment interactions. The Annals of Applied Statistics 8 176–203.
  * Willer, Li and Abecasis (2010) [author] Willer, Cristen JC. J., Li, YunY. and Abecasis, Gonçalo RG. R. (2010). METAL: fast and efficient meta-analysis of genomewide association scans. Bioinformatics 26 2190–2191.
  * Won et al. (2013) [author] Won, Joong-HoJ.-H., Lim, JohanJ., Kim, Seung-JeanS.-J. and Rajaratnam, BalaB. (2013). Condition-number-regularized covariance estimation. Journal of the Royal Statistical Society, Series B 75 427–450.
  * Wu et al. (2020) [author] Wu, YuluY., Cao, HongbaoH., Baranova, AnchaA., Huang, HailiangH., Li, ShengS., Cai, LeiL. et al. (2020). Multi-trait analysis for genome-wide association study of five psychiatric disorders. Translational Psychiatry 10 209\.
  * Xie and Stephens (2022) [author] Xie, DongyueD. and Stephens, MatthewM. (2022). Discussion of ”Confidence Intervals for Nonparametric Empirical Bayes Analysis”. Journal of the American Statistical Association 117 1186–1191.
  * Yang, Carbonetto and Stephens (2024a) [author] Yang, YunqiY., Carbonetto, PeterP. and Stephens, MatthewM. (2024a). R package and source code reproducing the analyses for “Improved methods for empirical Bayes multivariate multiple testing and effect size estimation”. 10.1214/provided-by-typesetter
  * Yang, Carbonetto and Stephens (2024b) [author] Yang, YunqiY., Carbonetto, PeterP. and Stephens, MatthewM. (2024b). Supplementary materials for “Improved methods for empirical Bayes multivariate multiple testing and effect size estimation”. 10.1214/provided-by-typesetter
  * Zhang et al. (2023) [author] Zhang, HaoyuH., Zhan, JiananJ., Jin, JinJ., Zhang, JingningJ., Lu, WenxuanW., Zhao, RuzhangR. et al. (2023). A new method for multiancestry polygenic prediction improves performance across diverse populations. Nature Genetics 55 1757–1768.
  * Zhao, Philip and Jiang (2008) [author] Zhao, J-HJ.-H., Philip, LHL. and Jiang, QibaoQ. (2008). ML estimation for factor analysis: EM or non-EM? Statistics and Computing 18 109–123.
  * Zhou and Stephens (2014) [author] Zhou, XiangX. and Stephens, MatthewM. (2014). Efficient multivariate linear mixed model algorithms for genome-wide association studies. Nature Methods 11 407–409.
  * Zou et al. (2024) [author] Zou, YuxinY., Carbonetto, PeterP., Xie, DongyueD., Wang, GaoG. and Stephens, MatthewM. (2024). Fast and flexible joint fine-mapping of multiple traits via the Sum of Single Effects model. bioRxiv doi:10.1101/2023.04.14.536893.