# A statistical framework for GWAS of high dimensional phenotypes using
summary statistics, with application to metabolite GWAS

Weiqiong Huang1, Emily C. Hector2, Joshua Cape3, Chris McKennan1,
1Department of Statistics, University of Pittsburgh
2Department of Statistics, North Carolina State University
3Department of Statistics, University of Wisconsin
To whom correspondence should be addressed (chm195@pitt.edu).

###### Abstract

The recent explosion of genetic and high dimensional biobank and ‘omic’ data
has provided researchers with the opportunity to investigate the shared
genetic origin (pleiotropy) of hundreds to thousands of related phenotypes.
However, existing methods for multi-phenotype genome-wide association studies
(GWAS) do not model pleiotropy, are only applicable to a small number of
phenotypes, or provide no way to perform inference. To add further
complication, raw genetic and phenotype data are rarely observed, meaning
analyses must be performed on GWAS summary statistics whose statistical
properties in high dimensions are poorly understood. We therefore developed a
novel model, theoretical framework, and set of methods to perform Bayesian
inference in GWAS of high dimensional phenotypes using summary statistics that
explicitly model pleiotropy, beget fast computation, and facilitate the use of
biologically informed priors. We demonstrate the utility of our procedure by
applying it to metabolite GWAS, where we develop new nonparametric priors for
genetic effects on metabolite levels that use known metabolic pathway
information and foster interpretable inference at the pathway level.

Keywords: High dimensional factor models; Bayesian nonparametrics; GWAS;
Summary statistics; Pleiotropy

##  1 Introduction

Genome-wide association studies (GWAS) investigate the relationship between
the genotype at single nucleotide polymorphisms (SNPs) and phenotypes, and
have become an essential tool for deciphering the genetic basis of human
variation [1]. Recently, the explosion of high dimensional biobank and ‘omic’
data have presented biologists with the opportunity to understand the common
genetic origin (pleiotropy) of hundreds to thousands of related phenotypes [2,
3, 4, 5]. Of particular interest are GWAS of metabolites, small molecules that
are by- or end-products of cell metabolism. Metabolites are crucial to
establishing functional links between genotype and disease [6]; thus,
metabolite GWAS (mtGWAS) have the potential to transform our understanding of
the etiology of disease and uncover new biomarkers for targeted therapies [7].

The many unsolved statistical and methodological challenges in mtGWAS
highlight the difficulties of performing GWAS with high dimensional
phenotypes. First, metabolite levels exhibit systemic genetic and non-genetic
correlation [8, 4]. However, existing methods for multi-phenotype GWAS that
incorporate inter-phenotype correlations can only accommodate a small number
of phenotypes [9, 10, 11, 12], can only test whether a SNP is related to at
least one phenotype [13, 14], or provide no inferential guarantees [2, 15].
Consequently, nearly all mtGWAS studies ignore pleiotropy and perform pairwise
marginal regressions with a stringent Bonferroni correction to adjust for the
#SNPs $\times$ #metabolites $\gtrsim 10^{8}$ tests. Second, raw metabolite and
genetic data are rarely available due to privacy concerns. Instead,
practitioners typically only have access to summary statistics from the #SNPs
$\times$ #metabolites marginal regressions of metabolite levels onto SNP
genotypes. Although existing approaches have treated SNPs as analogous to
biological samples and applied standard factor analysis tools to analyze the
matrix of summary statistics [2, 15], we show these approaches lead to
erroneous estimates and inference and that new methods and theoretical results
are required. Lastly, substantial biological information on the grouping of
metabolites into biological pathways should be leveraged to improve power and
interpretability. The latter is essential, since deciphering results from
mtGWAS with thousands of metabolites is exceedingly challenging, whereas
interpreting functional pathways is far simpler.

To address these critical obstacles, we develop a new statistical model,
theoretical framework, and methods for GWAS with high dimensional phenotypes.
Figure 1 gives a graphical depiction of our model, which decomposes the total
genetic effect into an indirect effect mediated by a set of latent factors and
a direct effect. Factors represent potentially genetically determined
biological processes, such as glucose metabolism in metabolomics or gut
architecture in microbiomics [16], that regulate the levels of many phenotypes
and consequently help model pleiotropy and non-genetic correlation. The direct
effect captures remaining idiosyncratic genetic variation.

Figure 1: A graphical model for the effect of genotype on phenotype levels.
Red: indirect effect; blue: direct effect. Parameters above arrows are the
parameters of interest.

Our theoretical framework is essential because our estimates rely on
performing factor analysis with GWAS summary statistics whose statistical
behavior is, to our knowledge, currently unknown. We therefore develop, to our
knowledge, the first set of theoretical guarantees for factor analysis with
summary statistics. These are more than mere technical developments, as we
show that ostensibly trivial dependencies between SNP genotypes cause
estimators’ asymptotic distributions to differ substantially from those
derived from raw data.

We lastly develop two new methods to estimate the indirect and direct genetic
effects in Figure 1 that attempt to solve the abovementioned methodological
challenges. The first is a method to estimate the number of latent factors we
call dBEMA (dependent Bulk Eigenvalue Matching), which accounts for unobserved
dependencies between SNP genotypes that cause existing methods designed for
raw data to drastically misestimate the number of factors when applied to GWAS
summary statistics. The second is a general Bayesian factor analysis method
called HiGSS (High Dimensional GWAS with Summary Statistics), which performs
theoretically justified Bayesian inference on model parameters. Unlike
traditional computationally onerous posterior sampling procedures and sub-
optimal variational approximations, HiGSS leverages our abovementioned
theoretical results to efficiently determine posteriors via empirical Bayes
without sacrificing statistical fidelity. Critically, HiGSS facilitates the
use of biologically-informed priors that would be computationally intractable
in traditional Bayesian factor analysis pipelines, such as phylogenetic and
spatial priors in microbiome and brain image data [17, 18]. While we design
HiGSS to be compatible with general priors, we devote considerable efforts to
constructing appropriate nonparametric priors for the direct and indirect
effects in mtGWAS, which leverage known metabolic pathway information to
facilitate powerful and interpretable inference at the pathway level.

The remainder of this paper is organized as follows. We describe our model for
summary statistics in Section 2. The proposed estimation algorithms and their
theoretical guarantees are presented in Section 3. In Section 4, we construct
nonparametric priors for indirect and direct effects in mtGWAS. Section 5
investigates the finite sample performance of our proposed estimation
approach. We perform a mtGWAS in Section 6 to illustrate the application of
our method. Section 7 concludes.

##  2 A model for the data

###  2.1 Problem Notation

For $\bm{Z}\in\mathbb{R}^{n\times m}$, we let $\bm{Z}_{i*}\in\mathbb{R}^{m}$,
$\bm{Z}_{*j}\in\mathbb{R}^{n}$, and $\bm{Z}_{ij}\in\mathbb{R}$ be the $i$th
row, $j$th column, and $(i,j)$th element of $\bm{Z}$. We use standard matrix
normal notation, i.e. for $\bm{M}\in\mathbb{R}^{n\times m}$,
$\bm{U}\in\mathbb{R}^{n\times n}$, and $\bm{V}\in\mathbb{R}^{m\times m}$,
$\bm{Z}\sim MN(\bm{M},\bm{U},\bm{V})$ if $\operatorname{vec}(\bm{Z})\sim
N(\operatorname{vec}(\bm{M}),\bm{V}\otimes\bm{U})$. Vector
$\bm{x}\sim(\bm{\mu},\bm{A})$ if $\mathbb{E}(\bm{x})=\bm{\mu}$ and
$\operatorname{Var}(\bm{x})=\bm{A}$.

###  2.2 A model for the unobserved genotype and phenotype data

Let $\bm{Y}\in\mathbb{R}^{N\times M}$ be the unobserved levels of $M$
phenotypes (metabolites) measured in $N$ unrelated individuals and
$\bm{X}\in\mathbb{R}^{N\times S}$ be the unobserved genotypes at $S$ SNPs in
the same individuals. Assume

###### Assumption 2.1

$\bm{X}$ is a random matrix with bounded, mean zero, and independent entries
whose rows are identically distributed and
$\operatorname{Var}(\bm{X}_{is})\geq\epsilon$ for some constant $\epsilon>0$
and all $i\in\\{1,\ldots,n\\}=[n]$ and $s\in[S]$.

The mean 0 and independence assumption are equivalent to mean-centering
$\bm{X}$’s columns and pruning SNPs for linkage disequilibrium, where pruning
helps avoid inferring redundant genetic effects. Figure 1 provides a graphical
description of our model for the effect of genotypes on phenotypes, which
decomposes the total effect into a direct effect and an indirect effect
mediated by $K$ latent factors $\bm{C}\in\mathbb{R}^{N\times K}$, and is
equivalent to the following model:

$\displaystyle\begin{aligned} \bm{Y}_{N\times M}=&\bm{C}_{N\times
K}\bm{L}_{M\times K}^{\top}+\bm{X}_{N\times S}\bm{\Delta}_{S\times
M}+\bm{E}_{N\times M},\quad\bm{E}\sim
MN(0,I_{N},\operatorname{diag}(\sigma_{1}^{2},\ldots,\sigma_{M}^{2}))\\\
\bm{C}_{N\times K}=&\bm{X}_{N\times S}\bm{G}_{S\times K}+\bm{\Xi}_{N\times
K},\quad\bm{\Xi}_{1*},\ldots,\bm{\Xi}_{N*}\operatorname{\stackrel{{\scriptstyle
i.i.d}}{{\sim}}}(0,\bm{\Psi}_{K\times K}),\end{aligned}$ (2.1)

where $\bm{L}_{mk}$ is the effect of the $k$th latent factor on the $m$th
phenotype and $\bm{G}_{sk}$ is the effect of the $s$th SNP on the $k$th
factor. Under this model, $\bm{G}$, $\bm{L}$, and $\bm{\Delta}$ are the
parameters of interest and
$\mathbb{E}(\bm{Y}\mid\bm{X},\bm{G},\bm{L},\bm{\Delta})=\bm{X}(\bm{G}\bm{L}^{\top}+\bm{\Delta})$
is the total genetic effect, where $\bm{G}_{s*}^{\top}\bm{L}_{m*}$ and
$\bm{\Delta}_{sm}$ are the indirect and direct effects of the $s$th SNP on the
$m$th phenotype. The non-genetic covariance between phenotypes
$\operatorname{Var}(\bm{Y}_{i*}\mid\bm{X},\bm{G},\bm{L},\bm{\Delta})=\bm{L}\bm{\Psi}\bm{L}^{\top}+\operatorname{diag}(\sigma_{1}^{2},\ldots,\sigma_{M}^{2})$
is low rank plus diagonal, which is a standard assumption in metabolite,
microbiome, and other multi-phenotype data [19, 8, 2]. Since we propose a
Bayesian method to infer $\bm{G}$, $\bm{L}$, and $\bm{\Delta}$, we explicitly
condition on them in the above and below probability statements to avoid
confusion.

Two features of (2.1) help motivate assumptions on $\bm{G}$ and $\bm{\Delta}$.
First, the first term in the total genetic coefficient
$\bm{G}\bm{L}^{\top}+\bm{\Delta}$ captures patterns of inter-phenotype and
inter-SNP effect sharing, suggesting $\bm{\Delta}$ should reflect
idiosyncratic genetic variation. Second, both $\bm{G}$ and $\bm{\Delta}$
represent direct effects of genotype on a set of phenotypes, and are
consequently expected to be sparse with small non-zero entries [20]. We
therefore assume throughout that

###### Assumption 2.2

For some positive integer $A$ and constants $c,r\geq 0$

  1. (a)

$\bm{\Delta}$ is a random matrix with independent entries that satisfy
$\mathbb{P}(\bm{\Delta}_{sm}\neq 0)=\pi_{sm}$ and
$\bm{\Delta}_{sm}\mid(\bm{\Delta}_{sm}\neq
0)\sim\sum_{a=1}^{A}\gamma_{sma}N(0,\tau_{sma}^{2})$, where
$\pi_{sm},\gamma_{sma}\in[0,1]$ and $\sum_{a=1}^{A}\gamma_{sma}=1$. Further,
$\max_{s\in[S];m\in[M];a\in[A]}(N^{1/2}\tau_{sma})\leq c$ and
$\max_{s\in[S];m\in[M]}\pi_{sm}=o(1)$ as $S,M\to\infty$.

  2. (b)

$\max_{s\in[S]}\lVert N^{1/2}\bm{G}_{s*}\rVert_{2}=O_{P}(\log^{r}N)$ and
$S^{-1}\sum_{s=1}^{S}\lVert N^{1/2}\bm{G}_{s*}\rVert_{2}^{2}=o_{P}(1)$ as
$S,N\to\infty$.

$\bm{\Delta}$’s entry-wise independence implies it captures idiosyncratic
genetic variation which, along with its entries’ symmetric spike and slab
distributions, are classic assumptions on genetic effects in GWAS with one and
more than one phenotype [21, 22, 23, 20, 24, 25]. The mixture-normal slab is
quite general and can approximate a large class of symmetric distributions
[26]. The parameters $\tau_{sma}$ determine the magnitude of $\bm{\Delta}$’s
non-zero effects, where the conditions on $N^{1/2}\tau_{sma}$ and
$N^{1/2}\bm{G}_{s*}$ reflect the observation that GWAS z-scores, whose
magnitudes scale with $N^{1/2}$ times the SNPs’ effect sizes, are typically
small to moderate in real data [22, 20]. Sparsity is implied by our condition
on $\pi_{sm}$ in (a) and the sum of squares constraint in (b). It is important
to note that while we assume $\bm{\Delta}$ is random, weak convergence results
for all frequentist estimators in Section 3.4 are stated conditional on
$\bm{\Delta}$. Assuming $\bm{\Delta}$ is random can be thought of as a way to
only consider $\bm{\Delta}$’s likely to arise in genetic data and ignore
unrealistic “corner cases”. We discuss assumptions for (2.1)’s remaining
parameters in Section 3.4.

###  2.3 Observed GWAS summary statistics

We assume the GWAS summary statistics
$\hat{\beta}_{sm}=(\bm{X}_{*s}^{\top}\bm{X}_{*s})^{-1}\bm{X}_{*s}^{\top}\bm{Y}_{*m}$
from the regressions of the $m$th phenotype onto $s$th SNP’s genotype are
observed by the user. We also assume $\bm{X}_{*s}^{\top}\bm{X}_{*s}$ is known,
although this can be estimated from public data if unavailable [27]. We define
the observable matrix of standardized GWAS statistics
$\hat{\bm{B}}\in\mathbb{R}^{S\times M}$ to be
$\hat{\bm{B}}_{sm}=(\bm{X}_{*s}^{\top}\bm{X}_{*s})^{1/2}\hat{\beta}_{sm}$
which, for
$\bm{D}=\operatorname{diag}\\{(\bm{X}_{*1}^{\top}\bm{X}_{*1}),\ldots,(\bm{X}_{*S}^{\top}\bm{X}_{*S})\\}$
and $\bm{R}=\bm{D}^{-1/2}\bm{X}^{\top}\bm{X}\bm{D}^{-1/2}$ (the in-sample
correlation between SNPs), can be expressed as

$\displaystyle\begin{aligned}
&\hat{\bm{B}}=\bm{D}^{-1/2}\bm{X}^{\top}\bm{Y}=\tilde{\bm{G}}\bm{L}^{\top}+\bm{R}(\bm{D}^{1/2}\bm{\Delta})+\tilde{\bm{E}}\\\
&\tilde{\bm{G}}=\bm{R}(\bm{D}^{1/2}\bm{G})+\bm{D}^{-1/2}\bm{X}^{\top}\bm{\Xi},\quad\tilde{\bm{E}}\mid\bm{X}\sim
MN(0,\bm{R},\operatorname{diag}(\sigma_{1}^{2},\ldots,\sigma_{M}^{2})).\end{aligned}$
(2.2)

Standardizing summary statistics by $(\bm{X}_{*s}^{\top}\bm{X}_{*s})^{1/2}$ is
akin to working with z-scores and begets rows (SNPs) with identical
covariances, where $\bm{D}^{1/2}\bm{G}$ and $\bm{D}^{1/2}\bm{\Delta}$ are the
standardized analogues of $\bm{G}$ and $\bm{\Delta}$. Although the diagonal
elements of $\bm{D}$ grow with the sample size, Assumption 2.2 implies
$\bm{D}^{1/2}\bm{G}$ and $\bm{D}^{1/2}\bm{\Delta}$ are moderate. Evidently,
$\hat{\bm{B}}$ depends on $\bm{D}^{1/2}\bm{G}$ and $\bm{D}^{1/2}\bm{\Delta}$
through $\bm{R}(\bm{D}^{1/2}\bm{G})$ and $\bm{R}(\bm{D}^{1/2}\bm{\Delta})$.
However, since the off-diagonal elements of $\bm{R}$ are small under
Assumption 2.1, the sparsity assumptions in Assumption 2.2 imply
$\bm{R}(\bm{D}^{1/2}\bm{G})\approx\bm{D}^{1/2}\bm{G}$ and
$\bm{R}(\bm{D}^{1/2}\bm{\Delta})\approx\bm{D}^{1/2}\bm{\Delta}$.

###  2.4 Parameter identifiability

The parameters of interest $(\bm{G},\bm{L})$ in (2.2) are not identifiable,
since $(\bm{G}\bm{A},\bm{L}\bm{A}^{-\top})$ will give an identical likelihood
for any invertible $K\times K$ matrix $\bm{A}$. We address this with the
following proposition.

###### Proposition 2.1

If Assumptions 2.1, 2.2, and 3.1 in Section 3.4 hold, then for all large
$N,M,S$:

  1. (i)

There exists a $(\bm{G},\bm{L})$ so that
$\mathbb{E}(S^{-1}\tilde{\bm{G}}^{\top}\tilde{\bm{G}}\mid\bm{G})=I_{K}$,
$M^{-1}\bm{L}^{\top}\bm{L}=\operatorname{diag}(\lambda_{1},\ldots,\lambda_{K})$
and $\lambda_{1}>\cdots>\lambda_{K}>0$. Further, the columns of $\bm{G}$ and
$\bm{L}$ are unique up to sign.

  2. (ii)

If $c^{-1}I_{K}\preceq
E(S^{-1}\tilde{\bm{G}}^{\top}\tilde{\bm{G}}\mid\bm{G})\preceq cI_{K}$ for some
constant $c>1$, then for all $\epsilon\in(0,1)$, there exists a
$\delta=\delta({\epsilon})>0$ such that
$\mathbb{P}(\lvert\tilde{\bm{G}}_{sk}\rvert\geq\delta\mid\bm{G})\geq
1-\epsilon$ for all $s\in[S]$, $k\in[K]$.

We use the parametrization in (i) for the remainder of the manuscript, which
is a classic way to identify components of factor models [28]. The
orthogonality of the columns of $\tilde{\bm{G}}$ and $\bm{L}$ also help
motivate Algorithm 1 below, which uses singular value decomposition to derive
estimates for $\bm{G}$, $\bm{L}$, and $\bm{\Delta}$. Part (ii) of Proposition
2.1 suggests that regardless of parametrization, $\tilde{\bm{G}}$ will never
be sparse even if $\bm{G}$ is sparse. This follows because $\tilde{\bm{G}}$ in
(2.2) is a noisy version of $\bm{G}$, whose error
$\bm{D}^{-1/2}\bm{X}^{\top}\bm{\Xi}$ is dense. This is in stark contrast with
existing methods that assume a sparse low rank structure on $\hat{\bm{B}}$
[15], suggesting they enforce non-existent sparsity.

##  3 A general method to estimate model parameters

###  3.1 HiGSS: an algorithm to estimate indirect and direct effects

Here we present a Bayesian method we call HiGSS (High Dimensional GWAS with
Summary Statistics) in Algorithm 1 below to derive the posterior
$\operatorname{pr}(\bm{G},\bm{L},\bm{\Delta}\mid\hat{\bm{B}},\bm{D})$ assuming
the number of latent factors $K$ is known. To circumvent the computational and
statistical issues that arise in Bayesian factor analysis, and as alluded to
by Proposition 2.1, HiGSS first uses singular value decomposition to derive
frequentist estimates for the parameters of interest and subsequently
leverages the asymptotic theory developed Section 3.4 to calculate the
posterior. HiGSS assumes the priors on $\bm{G}$, $\bm{L}$, and $\bm{\Delta}$
are independent.

######  Algorithm 1 (HiGSS)

Input data: Standardized GWAS estimates $\hat{\bm{B}}\in\mathbb{R}^{S\times
M}$,
$\bm{D}=\operatorname{diag}\\{(\bm{X}_{*1}^{\top}\bm{X}_{*1}),\allowbreak\ldots,\allowbreak$
$(\bm{X}_{*S}^{\top}\bm{X}_{*S})\\}$, number of latent factors $K\geq 0$, and
sample size $N$.
Input priors:
$\operatorname{pr}(\bm{G}\mid\bm{D})=\prod_{s=1}^{S}\operatorname{pr}_{s}(\bm{G}_{s*}\mid\bm{D})$,
$\operatorname{pr}(\bm{L}\mid\bm{D})=\prod_{m=1}^{M}\operatorname{pr}_{m}(\bm{L}_{m*})$,
and $\bm{\Delta}$-specific hyperparameters
$\\{\\{\pi_{sm},\\{(\gamma_{sma},\tau_{sma}^{2})\\}_{a\in[A]}\\}\\}_{s\in[S];m\in[M]}$.
Output: Posterior distribution
$\operatorname{pr}(\bm{G},\bm{L},\bm{\Delta}\mid\hat{\bm{B}},\bm{D})$

  1. 1.

Let $\bm{\Gamma}=\operatorname{diag}(\gamma_{1},\ldots,\gamma_{K})$,
$\bm{U}\in\mathbb{R}^{S\times K}$, and $\bm{V}\in\mathbb{R}^{M\times K}$ be
$\hat{\bm{B}}$’s first $K$ singular values and left and right singular
vectors. Define $\hat{\bm{G}}=S^{1/2}\bm{U}$,
$\hat{\bm{L}}=S^{-1/2}\bm{V}\bm{\Gamma}$, and
$\hat{\bm{\Delta}}=\hat{\bm{B}}-\hat{\bm{G}}\hat{\bm{L}}^{\top}$ to be
estimators for $\bm{G}$, $\bm{L}$, and $\bm{\Delta}$.

  2. 2.

Calculate
$\operatorname{pr}(\bm{G}\mid\hat{\bm{G}},\bm{D})\propto\prod_{s=1}^{S}\mathcal{N}(\hat{\bm{G}}_{s*};\bm{D}_{ss}^{1/2}\bm{G}_{s*},I_{K})\operatorname{pr}_{s}(\bm{G}_{s*}\mid\bm{D})$,
where $\mathcal{N}(\bm{x};\bm{\mu},\bm{V})$ is the density at $\bm{x}$ of a
normal with mean $\bm{\mu}$ and variance $\bm{V}$.

  3. 3.

Let $\hat{\sigma}_{m}^{2}=S^{-1}\sum_{s=1}^{S}\hat{\bm{\Delta}}_{sm}^{2}$ be
an estimate for $\sigma_{m}^{2}$. Calculate
$\operatorname{pr}(\bm{L}\mid\hat{\bm{L}},\bm{D})\propto\prod_{m=1}^{M}\mathcal{N}\\{\hat{\bm{L}}_{m*};\bm{L}_{m*},\allowbreak(N^{-1}+S^{-1})\hat{\sigma}_{m}^{2}I_{K}\\}\operatorname{pr}_{m}(\bm{L}_{m*})$.

  4. 4.

Let
$\operatorname{pr}(\bm{\Delta}\mid\hat{\bm{\Delta}},\bm{D})\propto\prod_{m=1}^{M}\prod_{s=1}^{S}\mathcal{N}(\hat{\bm{\Delta}}_{sm};\bm{D}_{ss}^{1/2}\bm{\Delta}_{sm},\hat{\sigma}_{m}^{2})\operatorname{pr}(\bm{\Delta}_{sm}\mid\pi_{sm},\tau_{sm}^{2})$,
where
$\bm{\Delta}_{sm}\mid(\pi_{sm},\tau_{sm}^{2})\sim(1-\pi_{sm})\delta_{0}+\pi_{sm}\sum_{a=1}^{A}\gamma_{sma}N(0,\tau_{sma}^{2})$.

  5. 5.

Return the posterior
$\operatorname{pr}(\bm{G},\bm{L},\bm{\Delta}\mid\hat{\bm{B}},\bm{D})=\operatorname{pr}(\bm{G}\mid\hat{\bm{G}},\bm{D})\operatorname{pr}(\bm{L}\mid\hat{\bm{L}},\bm{D})\operatorname{pr}(\bm{\Delta}\mid\hat{\bm{\Delta}},\bm{D})$.

The intuition for step 1 is if $\bm{\Delta}$ is sparse, (2.2) and the accuracy
of singular value decomposition suggest $\hat{\bm{G}}\approx\tilde{\bm{G}}$,
$\hat{\bm{L}}\approx\bm{L}$, and the residuals $\hat{\bm{\Delta}}$ should be
approximately $\bm{D}^{1/2}\bm{\Delta}$ plus independent noise [28]. The
subscripts “$s$” and “$m$” in $\bm{G}$’s and $\bm{L}$’s prior distributions
indicate these priors may depend on SNP $s$\- and phenotype $m$-specific
parameters, and we justify the likelihoods for $\hat{\bm{G}}$, $\hat{\bm{L}}$,
and $\hat{\bm{\Delta}}$ from steps 2-4 in Section 3.4. Since $\bm{\Delta}$
contains $S\times M\gtrsim 10^{8}$ elements, computing its posterior is
intractable unless an appropriate prior is chosen. The spike and mixture-
normal slab prior in step 4, whose parameters’ definitions match those in
Assumption 2.2, is conjugate and begets tractable computation. Its concordance
with Assumption 2.2 is coincidental and is chosen to match the prevailing
choices of priors in GWAS [20].

Step 5 assumes the total posterior
$\operatorname{pr}(\bm{G},\bm{L},\bm{\Delta}\mid\hat{\bm{B}},\bm{D})$ can be
factored into the product of posteriors defined in the previous three steps.
To see why this is appropriate, note
$\\{\hat{\bm{G}},\hat{\bm{L}},\hat{\bm{\Delta}}\\}$ is a function of
$\hat{\bm{B}}$ and
$\hat{\bm{B}}=\hat{\bm{G}}\hat{\bm{L}}^{\top}+\hat{\bm{\Delta}}$, meaning the
total posterior is equivalent to

$\displaystyle\operatorname{pr}(\bm{G},\bm{L},\bm{\Delta}\mid\hat{\bm{G}},\hat{\bm{L}},\hat{\bm{\Delta}},\bm{D})\propto$
$\displaystyle\operatorname{pr}(\hat{\bm{G}},\hat{\bm{L}},\hat{\bm{\Delta}}\mid\bm{G},\bm{L},\bm{\Delta},\bm{D})\operatorname{pr}(\bm{G}\mid\bm{D})\operatorname{pr}(\bm{L}\mid\bm{D})\operatorname{pr}(\bm{\Delta}\mid\bm{D})$

Step 5 therefore follows provided $\hat{\bm{G}}$, $\hat{\bm{L}}$, and
$\hat{\bm{\Delta}}$ are independent conditional on
$\\{\bm{G},\bm{L},\bm{\Delta},\bm{D}\\}$, which we prove is true
asymptotically in Section 3.4.

Perhaps the most unexpected component of Algorithm 1 is $\hat{\bm{L}}$’s
likelihood in step 3. Since $\hat{\bm{B}}$ has $S$ rows, existing singular
value decomposition results suggest $\hat{\bm{L}}_{m*}$’s asymptotic variance
should be $S^{-1}\hat{\sigma}_{m}^{2}I_{K}$ because $\hat{\bm{L}}_{m*}$ is
ostensibly estimated on $S$ degrees of freedom [28]. This is quite different
from the variance in step 3 which, because $N$ is typically smaller than $S$,
is closer to $N^{-1}\hat{\sigma}_{m}^{2}I_{K}$ than
$S^{-1}\hat{\sigma}_{m}^{2}I_{K}$. This is sensible because $\hat{\bm{B}}$ in
(2.2) is obtained by projecting the $N$ samples up to the space of $S$ SNPs,
implying $\hat{\bm{L}}_{m*}$ is actually estimated on something closer to $N$
degrees of freedom. Theorem 3.2 and Corollary 3.1 in Section 3.4 below, which
prove the veracity of $\hat{\bm{L}}$’s likelihood, are to our knowledge the
first results of their kind and, as we show in Section 5, are critical to the
fidelity of inference on $\bm{L}$.

Algorithm 1 is general in that the likelihoods in steps 2-4 and form for the
posterior in step 5 are the same regardless of the phenotypes being analyzed.
However, the choice of priors will likely depend on the type of phenotypes.
For example, the set of priors chosen in metabolomic applications may be
different that those in microbiome data. We therefore provide general
recommendations for how to choose priors for $\bm{G}$ and $\bm{L}$ and
$\bm{\Delta}$’s hyperparameters in Section 3.2.

###  3.2 Choice of priors

First, recall $\bm{G}$ is the effect of genotype on a small number ($K$) of
latent phenotypes. Since $\hat{\bm{G}}$’s likelihood is normal in Algorithm 1,
we can either treat $\bm{G}$’s $K$ columns as independent and use the many
single-phenotype GWAS priors designed for and normal data (e.g. [22] and
[20]), or model $\bm{G}$’s columns dependencies using recent methods designed
multi-phenotype GWAS [24, 25]. We take the former approach in our metabolomic
application.

The likelihoods for $\hat{\bm{L}}$ and $\hat{\bm{\Delta}}$ in Algorithm 1
suggest the simplest approach would be to treat estimating $\bm{L}$ and
$\bm{\Delta}$ as a normal means problem and use empirical Bayes to derive
column-specific and matrix-wide priors, respectively, which can be done
efficiently using [26]. An alternative would be to incorporate prior knowledge
on the relationships between the $M$ phenotypes into the prior on $\bm{L}$ and
$\bm{\Delta}$’s hyperparameters. For example, it is straightforward to
incorporate phylogenetics in $\bm{L}$’s prior using the methods discussed in
[17], and we show in Section 4 how to leverage metabolic pathways in our prior
on $\bm{L}$ and $\bm{\Delta}$’s hyperparameters.

The likelihoods for $\hat{\bm{G}}$ and $\hat{\bm{L}}$ in Algorithm 1 can be
written as a product of likelihoods across their $K$ columns. Consequently,
provided their priors also factor across their columns, posterior computation
on $\bm{G}$ and $\bm{L}$ can be done in parallel across their columns, which
is not true for traditional Bayesian or recent empirical Bayes-based factor
analysis methods [29, 30].

###  3.3 Determining the number of latent factors

Estimating $K$ is critical, since underestimating it can cause one to miss
important genetic variation and overestimating it by too much can invalidate
statistical guarantees. Given $\hat{\bm{B}}$’s model in (2.2), it is tempting
to reason that because the off-diagonals of $\bm{R}$ (the in-sample
correlation between SNPs) are approximately zero, one can treat the entries of
the error matrix $\tilde{\bm{E}}$ as independent and use one of the many
estimators designed for omic data with independent errors [31, 32, 33, 34].
Unfortunately, as shown in Section 5, ignoring the dependencies produced by
$\bm{R}$’s off-diagonal elements causes these methods to fail.

One way to address the abovementioned issues could be to use [35], which
explicitly models the dependence. However, applying their method in this
setting would require knowing $\bm{R}$, which is unobservable. A second option
might be to use methods whose theoretical results allow $\tilde{\bm{E}}$’s
entries be dependent, such as [36], [37], or [38]. However, the former two
require factor loadings be unrealistically large (i.e. satisfy the pervasive
factor assumption) and the latter patently fails in simulated and real omic
data [39]. To address these deficiencies, note that if we assume for
simplicity that $\bm{\Delta}$ in (2.2) is 0,

$\displaystyle\hat{\bm{B}}=\underbrace{\tilde{\bm{G}}\bm{L}^{\top}}_{\text{rank-$K$
signal}}+\underbrace{\tilde{\bm{E}}}_{\text{noise}}.$

Consequently, any of $\hat{\bm{B}}$’s singular values that exceed those in a
pure noise model must have arisen from the low rank signal. This suggests that
if we knew $\tilde{\bm{E}}$’s maximum singular value $s_{\max}$, a natural
estimator for $K$ is the number of $\hat{\bm{B}}$’s singular values that
exceed $s_{\max}$. We use this reasoning to motivate our novel method dBEMA
(dependent Bulk Eigenvalue Matching), which models $\hat{\bm{B}}$’s bulk
spectrum to predict $\tilde{\bm{E}}$’s maximum singular value. We describe
dBEMA below.

Our method is an extension of [34] to the setting where $\tilde{\bm{E}}$ has
dependent rows. To be consistent with Section 2, we assume throughout that
$\bm{E}$ in (2.1), and therefore $\tilde{\bm{E}}$ in (2.2), is Gaussian.
However, the universality theory of eigenvalues implied by [38] indicates we
can relax this to only require $\bm{E}$’s entries satisfy weaker moment
conditions.

To motivate the method, suppose $S=S(N)$ and $M=M(N)$ such that
$S/N\to\gamma_{S}\in(0,\infty)$ and $M/N\to\gamma_{M}\in(0,\infty)$, and let
$F_{R,N}$ and $F_{\sigma^{2},N}$ be the empirical distributions of $\bm{R}$’s
$S$ eigenvalues and $\\{\sigma_{m}^{2}\\}_{m\in[M]}$, respectively. If the
distributions $\lim_{N\to\infty}F_{R,N}=F_{R}$ and
$\lim_{N\to\infty}F_{\sigma^{2},N}=F_{\sigma^{2}}$ exist and are supported on
compact intervals, Lemmas 1 and 3 in [38] imply $F_{R}$ and $F_{\sigma^{2}}$
determine the behavior of $\tilde{\bm{E}}$’s bulk and largest singular values.
The problem therefore boils down to estimating $F_{R}$ and $F_{\sigma^{2}}$,
since this will allow us to predict $\tilde{\bm{E}}$’s maximum singular value.

First, Assumption 2.1 implies $F_{R}$ is exactly the Marchenko–Pastur
distribution [40], which for $x\geq 0$,
$\lambda_{-}=(1-\sqrt{\gamma_{S}})^{2}$, and
$\lambda_{+}=(1+\sqrt{\gamma_{S}})^{2}$ takes the form

$\displaystyle
F_{R}(x)=(1-\gamma_{S}^{-1})1\\{\gamma_{S}>1\\}+\mathop{\smallint}\limits_{0}^{x}f_{R}(t)dt,\quad
f_{R}(t)=(2\pi\gamma_{S}t)^{-1}\sqrt{(\lambda_{+}-t)(t-\lambda_{-})}1\\{t\in[\lambda_{-},\lambda_{+}]\\}.$

Unlike $F_{R}$, $F_{\sigma^{2}}$ is unknown and must be estimated from the
data. Following [34], we model $F_{\sigma^{2}}$ as a gamma distribution and
use $\hat{\bm{B}}$’s bulk singular values to estimate its parameters. Briefly,
let $\hat{s}_{\min(N,M,S)}\leq\cdots\leq\hat{s}_{1}$ be $\hat{\bm{B}}$’s non-
zero singular values and define
$\mathcal{B}=\\{\hat{s}_{\lfloor(\alpha/2)\min(N,M,S)\rfloor},\ldots,\hat{s}_{\lfloor(1-\alpha/2)\min(N,M,S)\rfloor}\\}$
for some $\alpha\in(0,1)$ to be $\hat{\bm{B}}$’s bulk singular values. Since
$\mathcal{B}$ excludes $\hat{\bm{B}}$’s $K$ largest singular values for
sufficiently large $N$, $\mathcal{B}$’s elements are approximately the inner
$1-\alpha$ fraction of $\tilde{\bm{E}}$’s singular values. Consequently, these
can be predicted given $F_{R}$ and the rate and scale parameters that
parameterize $F_{\sigma^{2}}$. We therefore minimize the difference between
$\mathcal{B}$’s elements and their predicted values as a function of the rate
and scale parameters to estimate $F_{\sigma^{2}}$. Equipped with $F_{R}$ and
an estimate for $F_{\sigma^{2}}$, we lastly determine $\hat{F}_{s_{max}}$, an
estimate for the distribution of $\tilde{\bm{E}}$’s maximum singular value,
and define our estimate for the number of factors to be the number
$\hat{\bm{B}}$’s singular values that exceed the $1-\beta$ quantile of
$\hat{F}_{s_{max}}$ for some user-specified $\beta\in(0,1)$. The algorithm is
given in Section A1.

As mentioned above, the theoretical work used to motivate dBEMA requires
$F_{\sigma^{2}}$ have bounded support, i.e.
$\min\\{x\in\mathbb{R}\cup\\{\infty\\}:F_{\sigma^{2}}(x)=1\\}<\infty$, which
belies our decision to model $F_{\sigma^{2}}$ as a gamma distribution. We use
simulated data in Section 5 to show dBEMA performs well despite this.

###  3.4 Theoretical justifications

Here we provide a theoretical justification for the likelihoods for
$\hat{\bm{G}}$, $\hat{\bm{L}}$, and $\hat{\bm{\Delta}}$ used in Algorithm 1.
Theory for the consistency of our estimate for $K$ proposed in Section 3.3
will appear in a subsequent version of this manuscript. We first present an
assumption on the parameters in (2.1).

###### Assumption 3.1

Let $c>1$ and $\epsilon\in(0,1)$ be constants and define $\lambda_{K+1}=0$.

  1. (a)

$N/S,M/S\in[c^{-1},c]$, $K\leq c$, $\sigma_{m}^{2}\in[c^{-1},c]$ for all
$m\in[M]$, and $\bm{\Xi}_{1k}$ is sub-exponential for all $k\in[K]$.

  2. (b)

$\mathbb{E}(S^{-1}\tilde{\bm{G}}^{\top}\tilde{\bm{G}}\mid\bm{G})=I_{K}$ and
$M^{-1}\bm{L}^{\top}\bm{L}=\operatorname{diag}(\lambda_{1},\ldots,\lambda_{K})\succ
0$ such that $\lambda_{k}/\lambda_{k+1}\geq 1+c^{-1}$ and
$\lambda_{k}\in[c^{-1}N^{\epsilon-1},c]$ for all $k\in[K]$ and
$\lvert\bm{L}_{mk}\rvert\leq c\lambda_{k}^{1/2}$ for all $m\in[M]$ and
$k\in[K]$.

###### Remark 3.1

The conditions in (b) implicitly assume $\bm{G}$ and $\bm{L}$ are non-random,
although we continue to condition on them in all probability statements below
to avoid any confusion. It is straightforward to replace (b) we an equivalent
assumption that allows $\bm{G}$ and $\bm{L}$ to be random.

We assume $N$, $M$, and $S$ are comparable in (a) to be consisted with ‘omic’
and biobank data. The requirement that
$\mathbb{E}(S^{-1}\tilde{\bm{G}}^{\top}\tilde{\bm{G}}\mid\bm{G})=I_{K}$ and
$M^{-1}\bm{L}^{\top}\bm{L}$ be diagonal with decreasing diagonal elements is
without loss of generality by Proposition 2.1, where our condition on the
ratio of successive diagonal elements is a standard eigengap assumption and
ensures we can identify the columns of $\bm{L}$ and $\bm{G}$ up to sign [28].
The parameter $\lambda_{k}$ is the average sum of squares of the $k$th loading
$\bm{L}_{*k}$ and is interpretable as the $k$th latent factor’s signal
strength, where the “pervasive factor assumption” is satisfied if
$\lambda_{1},\ldots,\lambda_{K}\asymp 1$ [41]. While common in the literature,
the pervasive factor assumption is patently violated in ‘omic’ data, where
methods that utilize it often cannot recover non-pervasive, but critical,
sources of omic variation [39]. This is especially important in our data
application whose most genetically regulated factors are weaker. We therefore
allow $\lambda_{k}$ to be large ($\lambda_{k}\asymp 1$, pervasive) or moderate
to small ($\lambda_{k}<<1$, non-pervasive) to ensure we can perform inference
on as many types of factors as possible.

We first justify $\hat{\bm{G}}^{\prime}s$ likelihood in step 2 of Algorithm 1.
We define $\mathscr{C}_{d}$ and $\Phi_{d}$ to be the set of all convex sets
and the distribution function for the standard normal in $\mathbb{R}^{d}$,
respectively, for the remainder of the section.

###### Theorem 3.1

Suppose Assumptions 2.1, 2.2, and 3.1 hold. Then as $N,M,S\to\infty$, there
exist $a_{1},\ldots,a_{K}\in\\{-1,1\\}$ so that for any fixed set
$\mathcal{S}\subset[S]$,

$\displaystyle\sup_{\begin{subarray}{c}R_{s}\in\mathscr{C}_{K}\end{subarray}}\left\lvert\mathbb{P}\left[\mathop{\bigcap}\limits_{s\in\mathcal{S}}\\{\hat{\bm{G}}_{s*}-\bm{D}_{ss}^{1/2}\operatorname{diag}(a_{1},\ldots,a_{K})\bm{G}_{s*}\\}\in
R_{s}\mid\bm{G},\bm{L},\bm{\Delta},\bm{X}\right]-\prod_{s\in\mathcal{S}}\Phi_{K}(R_{s})\right\rvert=o_{P}(1).$

Since the sign multipliers $a_{1},\ldots,a_{K}$ are the same for each SNP $s$,
this shows we can estimate each column of $\bm{G}$ up to sign. Theorem 3.1
also justifies treating $\hat{\bm{G}}$’s rows in Algorithm 1 as independent
standard normals. We next justify step 3 of Algorithm 1.

###### Theorem 3.2

Suppose Assumptions 2.1, 2.2, and 3.1 hold and $a_{1},\ldots,a_{K}$ are as
given in Theorem 3.1. Then the following hold as $N,M,S\to\infty$:

  1. (i)

$\hat{\sigma}_{m}^{2}=\sigma_{m}^{2}+o_{P}(1)$ for all $m\in[M]$.

  2. (ii)

$\mathop{\sup}\limits_{R\in\mathscr{C}_{1}}\left\lvert\mathbb{P}\left\\{\lVert\hat{\bm{L}}_{m*}\rVert_{2}^{-1}(\hat{\bm{G}}_{s*}^{\top}\hat{\bm{L}}_{m*}-\bm{D}_{ss}^{1/2}\bm{G}_{s*}^{\top}\bm{L}_{m*})\in
R\mid\bm{G},\bm{L},\bm{\Delta},\bm{X}\right\\}-\Phi_{1}(R)\right\rvert=o_{P}(1)$
for any $s\in[S]$ and $m\in[M]$ provided
$N^{1/2}\lVert\bm{L}_{m*}\rVert_{2}\to\infty$.

  3. (iii)

$\mathop{\sup}\limits_{\begin{subarray}{c}R_{m}\in\mathscr{C}_{K}\end{subarray}}\left\lvert\mathbb{P}\left[\mathop{\bigcap}\limits_{m\in\mathcal{M}}\\{(N^{-1}+S^{-1})\hat{\sigma}_{m}^{2}\\}^{-1/2}(\hat{\bm{L}}_{m*}-\bm{A}\bm{L}_{m*})\in
R_{m}\mid\bm{G},\bm{L},\bm{\Delta},\bm{X}\right]-\prod\limits_{m\in\mathcal{M}}\Phi_{K}(R_{m})\right\rvert\allowbreak=o_{P}(1)$
for any fixed set $\mathcal{M}\subset[M]$ and some
$\bm{A}\in\mathbb{R}^{K\times K}$ such that
$\lVert\bm{A}-\operatorname{diag}(a_{1},\ldots,a_{K})\rVert_{2}=O_{P}(N^{-1/2})$.

Results (i) and (ii) imply $\hat{\sigma}_{m}^{2}$ in step 3 of Algorithm 1 is
consistent for $\sigma_{m}^{2}$ and we can infer indirect effects
$\bm{G}_{s*}^{\top}\bm{L}_{m*}$. The condition
$N^{1/2}\lVert\bm{L}_{m*}\rVert_{2}\to\infty$ will hold if phenotype $m$ is
non-trivially affected by at least one latent factor, and appears to be
satisfied for nearly all $m$ in our data application. Result (iii) implies
$\hat{\bm{L}}_{m*}$ is asymptotically normal with asymptotic variance as given
in step 3. However, $\hat{\bm{L}}_{m*}$’s asymptotic distribution is centered
around $\bm{A}\bm{L}_{m*}$ rather than
$\operatorname{diag}(a_{1},\ldots,a_{K})\bm{L}_{m*}$, which is a standard
property of singular value decomposition-based estimates [42]. Although their
difference is small, it may not be small enough to guarantee we can infer
$\bm{L}_{m*}$’s sparsity patterns, since $(\bm{A}\bm{L}_{m*})_{k}$ may be
larger than its estimate’s standard deviation even if $\bm{L}_{mk}=0$. This is
especially important when data contain both pervasive and non-pervasive
factors, since $\bm{A}$’s off-diagonal elements risk contaminating sparser,
non-pervasive factor loadings with dense pervasive factor loadings. Corollary
3.1 below shows that no such contamination occurs.

###### Corollary 3.1

In addition to the conditions of Theorem 3.2, suppose $\lambda_{k}\geq\eta$ or
$\lambda_{k}\to 0$ as $M\to\infty$ for all $k\in[K]$ and some constant
$\eta>0$. Define the number of pervasive factors $K_{\operatorname{pf}}$ to be
0 if $\lambda_{1}\to 0$ and $K_{\operatorname{pf}}\in[K]$ if
$\lambda_{K_{\operatorname{pf}}}\geq\eta$ but
$\lambda_{K_{\operatorname{pf}}+1}\to 0$ as $M\to\infty$. Then $\bm{A}$ in
Theorem 3.2 satisfies

$\displaystyle\bm{A}=\begin{cases}\operatorname{diag}(a_{1},\ldots,a_{K})&\text{if
$K_{\operatorname{pf}}=0$}\\\
\bm{A}_{\operatorname{pf}}\oplus\operatorname{diag}(a_{K_{\operatorname{pf}}+1},\ldots,a_{K})&\text{if
$K_{\operatorname{pf}}\in[K-1]$}\\\ \bm{A}&\text{if
$K_{\operatorname{pf}}=K$}\end{cases},\quad\lVert\bm{A}_{\operatorname{pf}}-\operatorname{diag}(a_{1},\ldots,a_{K_{\operatorname{pf}}})\rVert_{2}=O_{P}(N^{-1/2})$

and, for any fixed set $\mathcal{M}\subset[M]$ and non-pervasive factor
$k\in\\{K_{\operatorname{pf}}+1,\ldots,K\\}$,

$\displaystyle\mathop{\sup}\limits_{\begin{subarray}{c}R_{m}\in\mathscr{C}_{1}\end{subarray}}\left\lvert\mathbb{P}\left[\mathop{\bigcap}\limits_{m\in\mathcal{M}}\\{(N^{-1}+S^{-1})\hat{\sigma}_{m}^{2}\\}^{-1/2}(\hat{\bm{L}}_{mk}-a_{k}\bm{L}_{mk})\in
R_{m}\mid\bm{G},\bm{L},\bm{\Delta},\bm{X}\right]-\prod\limits_{m\in\mathcal{M}}\Phi_{1}(R_{m})\right\rvert\allowbreak=o_{P}(1).$

The expression for $\bm{A}$ implies we need not worry about the abovementioned
contamination, which gives rise to the subsequent result stating the
asymptotic distributions for non-pervasive factor loadings are, up to sign,
centered around their true loadings. While this is not necessarily true for
pervasive factor loadings, it has a negligible effect on inference because
contamination will only trivially impact their large non-zero entries and, as
illustrated in our data example, pervasive factors are more likely to reflect
non-genetic variation. We lastly consider steps 4 and 5 of Algorithm 1.

###### Theorem 3.3

Suppose Assumptions 2.1, 2.2, and 3.1 hold and $a_{1},\ldots,a_{K}$, $\bm{A}$
are as given in Theorem 3.1 and 3.2. Then the following hold as
$N,M,S\to\infty$ for any fixed sets $\mathcal{S}\subset[S]$,
$\mathcal{M}\subset[M]$:

  1. (i)

$\mathop{\sup}\limits_{\begin{subarray}{c}R_{sm}\in\mathscr{C}_{1}\end{subarray}}\left\lvert\mathbb{P}\left\\{\mathop{\bigcap}\limits_{\begin{subarray}{c}s\in\mathcal{S}\\\
m\in\mathcal{M}\end{subarray}}\hat{\sigma}_{m}^{-1}(\hat{\bm{\Delta}}_{sm}-\bm{D}_{ss}^{1/2}\bm{\Delta}_{sm})\in
R_{sm}\mid\bm{G},\bm{L},\bm{\Delta},\bm{X}\right\\}-\prod\limits_{\begin{subarray}{c}s\in\mathcal{S}\\\
m\in\mathcal{M}\end{subarray}}\Phi_{1}(R_{sm})\right\rvert\allowbreak=\allowbreak
o_{P}(1)$.

  2. (ii)

Define
$\bm{Z}_{s}^{\operatorname{(G)}}=\hat{\bm{G}}_{s*}-\bm{D}_{ss}^{1/2}\operatorname{diag}(a_{1},\ldots,a_{K})\bm{G}_{s*}$,
$\bm{Z}_{m}^{\operatorname{(L)}}=\\{(N^{-1}+S^{-1})\hat{\sigma}_{m}^{2}\\}^{-1/2}(\hat{\bm{L}}_{m*}-\bm{A}\bm{L}_{m*})$,
and
$z_{sm}^{(\Delta)}=\hat{\sigma}_{m}^{-1}(\hat{\bm{\Delta}}_{ms}-\bm{D}_{ss}^{1/2}\bm{\Delta}_{ms})$.
Then the sets of random variables
$\\{\bm{Z}_{s}^{\operatorname{(G)}}\\}_{s\in\mathcal{S}}$,
$\\{\bm{Z}_{m}^{\operatorname{(L)}}\\}_{m\in\mathcal{M}}$, and
$\\{z_{sm}^{(\Delta)}\\}_{(s,m)\in\mathcal{S}\times\mathcal{M}}$ are
asymptotically independent conditional on
$\\{\bm{G},\bm{L},\bm{\Delta},\bm{X}\\}$.

Part (i) justifies $\hat{\bm{\Delta}}$’s likelihood in Algorithm 1, which
assumes its entries are independent and normally distributed. The asymptotic
independence of the z-statistics in part (ii), which are the standardized
equivalents of $\hat{\bm{G}}$, $\hat{\bm{L}}$, and $\hat{\bm{\Delta}}$,
justifies step 5.

##  4 Pathway-centric priors and inference in metabolite GWAS

Algorithm 1 provides a general method to compute the posterior of our
parameters of interest, where the frequentist estimators in step 1,
likelihoods in steps 2-4, and form for the posterior in step 5 are consisted
across datasets. However, the choice of priors is domain-specific. In this
section we describe a set of priors specific to metabolomics, which also
highlights the power of our methodology to utilize biologically-informed
priors that are intractable in standard Bayesian factor analysis.

To motivate our priors, Figure 2(a) shows that metabolites can be
hierarchically partitioned into a set of super- and sub-pathways, where the
overt pathway-specific variation exhibited in Figure 2(b) suggests metabolites
in the same pathway behave similarly. We therefore follow Section 3.2’s
recommendations and design priors for $\bm{L}$ and $\bm{\Delta}$ that
incorporate metabolic pathway information.

Figure 2: (a): Hierarchical partition of metabolites in our data example.
Colors indicate each super-pathway’s set of sub-pathways and numbers give the
number of metabolites in each super-pathway. Each super-pathway’s two largest
sub-pathways are noted. (b): $\hat{\bm{L}}_{*2}$ in our data example. Red
lines demarcate super-pathways and sub-pathway colors match (a)’s.

###  4.1 Prior on indirect effects

Indirect effects of genotype on metabolite levels are parametrized by $\bm{G}$
and $\bm{L}$. Our prior on $\bm{G}$ follows Section 3.2’s discussion and
assumes each column’s entries are drawn from a spike and slab:

$\displaystyle\bm{G}_{sk}\mid(\pi_{G_{k}},\tau_{G_{k}}^{2},\bm{D})\sim(1-\pi_{G_{k}})\delta_{0}+\pi_{G_{k}}N(0,\tau_{G_{k}}^{2}/\bm{D}_{ss}),\quad
s\in[S];k\in[K],$ (4.1)

which is conjugate under $\hat{\bm{G}}$’s normal likelihood in Algorithm 1.
Standardizing by $\bm{D}_{ss}$ implies standardized effects
$\bm{D}_{ss}^{1/2}\bm{G}_{sk}$ are identically distributed across SNPs $s$ and
reflects the observation that genetic effects are typically inversely
proportional to minor allele frequency [43]. We estimate
$\pi_{G_{k}},\tau_{G_{k}}^{2}$ via empirical Bayes.

Our prior on $\bm{L}$ incorporates metabolic pathway information and is
constructed to reflect the following observations from Figure 2(b) and
inference goals. First, metabolite loadings from the same sub-pathway appear
to be drawn from the same distribution, where sub-pathway distributions within
a super-pathway are more similar than those from different super-pathways.
Second, Figure 2(b) suggests many sub-pathways have identical loading
distributions. For example, the first two orange and green lipid sub-pathways
may have the same distributions. Our prior should therefore facilitate a
clustering on sub-pathways, which couples biological processes to beget more
interpretable inference. Third, it is expected that some metabolites behave as
“outliers” because they were incorrectly partitioned or do not follow their
sub-pathways’ distributions.

To accommodate these observations, we use a hierarchical Dirichlet process
[44] to model sub-pathway distribution parameters. This begets a clustering
over sub-pathways, where two sub-pathways from the same super-pathway are a
priori more likely to lie in the same cluster than sub-pathways from different
super-pathways. Metabolite loadings are then either drawn from their sub-
pathway’s distribution or, to accommodate the above third observation, an
“outlier” distribution. Mathematically, let $p\in[P]$ index super-pathways,
$b\in[B_{p}]$ index constituent sub-pathways, and $m\in[M_{pb}]$ index
constituent metabolites. We assume $\bm{L}$ has independent columns and, to
simplify notation, suppress any dependence on the factor number. As such, we
let $\bm{L}_{pbm}\in\mathbb{R}$ be a factor’s loading for the $m$th metabolite
in the $b$th and $p$th sub- and super-pathway for the remainder of Section
4.1. Then for concentration parameters $\gamma,\alpha_{0}>0$ and probability
measure $H$ on $\mathbb{R}\times\mathbb{R}_{\geq 0}$,

$\displaystyle\begin{aligned}
&F_{0}\mid(\gamma,H)\sim\operatorname{DP}(\gamma,H),\quad\pi_{\operatorname{out}}\sim\text{Beta}(c,d)\\\
&\text{$F_{p}\mid(\alpha_{0},F_{0})\sim\operatorname{DP}(\alpha_{0},F_{0}),\quad
p\in[P]$}\\\ &\text{$(\mu_{pb},\phi^{2}_{pb})\mid F_{p}\sim F_{p},\quad
b\in[B_{p}]$}\\\
&\text{$\bm{L}_{pbm}\mid(\mu_{pb},\phi^{2}_{pb},\pi_{\operatorname{out}},a_{\operatorname{out}},b_{\operatorname{out}})\sim(1-\pi_{\operatorname{out}})N(\mu_{pb},\phi^{2}_{pb})+\pi_{\operatorname{out}}U[a_{\operatorname{out}},b_{\operatorname{out}}],\quad
m\in[M_{pb}]$}.\end{aligned}$ (4.2)

Figure 3(a) provides a graphical description of this prior. Assuming for the
moment that $\pi_{\operatorname{out}}=0$, $\bm{L}_{pbm}$ is drawn from a
hierarchical Dirichlet process mixture model, where mixture components
(clusters) are defined by the unique values of
$\\{(\mu_{pb},\phi^{2}_{pb})\\}_{p\in[P];b\in[B_{p}]}$ [44]. Briefly, the
probability measure $F_{0}$ is first drawn from a Dirichlet process with
baseline measure $H$ and, with probability one, takes the form
$\sum_{j=1}^{\infty}\beta_{0,j}\delta_{(\mu_{j},\phi^{2}_{j})}$ for some
probability weights $\beta_{0,j}$ and
$(\mu_{j},\phi^{2}_{j})\in\mathbb{R}\times\mathbb{R}_{\geq 0}$. Each $F_{p}$
is subsequently drawn from another Dirichlet process with baseline measure
$F_{0}$, meaning $F_{p}$ can also be expressed as
$\sum_{j=1}^{\infty}\beta_{p,j}\delta_{(\mu_{j},\phi^{2}_{j})}$ for some
probability weights $\beta_{p,j}$. The third line in (4.2) implies sub-
pathways within a super-pathway are assigned to clusters based on their
parameters $(\mu_{pb},\phi^{2}_{pb})$, where sub-pathways that share
parameters fall in the same cluster. However, since the $F_{p}$’s share atoms,
sub-pathways from different super-pathways can be assigned to the same
cluster, where larger values of $\alpha_{0}$ increase the probability of this
occurring.

Figure 3: A graphical model for $\bm{L}$’s prior (a) and $\bm{\Delta}$’s prior
(b). Square and circular nodes are fixed and random quantities. Nodes within a
rectangle are replicated.

The parameter $\pi_{\operatorname{out}}$ in (4.2) is the fraction of all
metabolites that are “outliers”, i.e. metabolites whose loadings do not follow
the hierarchical Dirichlet process mixture model. Since we assume their
loadings can take any value with equal probability, we treat them as uniformly
distributed and set $a_{\operatorname{out}},b_{\operatorname{out}}$ so
$[a_{\operatorname{out}},b_{\operatorname{out}}]$ contains all observed
loadings. We let $\pi_{\operatorname{out}}\sim\text{Beta}(c,d)$ so as to allow
the data to help dictate outlier frequency, and set $c=1/5,d=999/5$ so that
out of 1,000 metabolites, we a priori expect one outlier and no more than five
outliers with probability 0.95.

Similar to [44], we place diffuse priors on $\gamma$ and $\alpha_{0}$. We set
the base distribution $H$ so that
$\mathbb{P}\\{(\mu_{pb},\phi_{pb}^{2})=(0,0)\mid H\\}>0$ and
$\mu_{pb}\perp\\!\\!\\!\perp\phi_{pb}^{2}\mid\\{H,(\mu_{pb},\phi_{pb}^{2})\neq(0,0)\\}$.
The former reflects our belief that some sub-pathways will have loadings
exactly equal to 0 and we let
$\mu_{pb}\mid\\{H,(\mu_{pb},\phi_{pb}^{2})\neq(0,0)\\}$ follow a spike and
mixture-normal slab to ensure conjugacy, where the slab is chosen to
approximate the horseshoe prior [45]. We let
$\phi_{pb}^{2}\mid\\{H,(\mu_{pb},\phi_{pb}^{2})\neq(0,0)\\}$ be a densely
discretized gamma distribution, where discretizing facilitates faster
computation by letting us to integrate out $\bm{L}$ in our Gibbs sampler.
Section A2 provides additional details and describes our Gibbs sampler, which
utilizes the Chinese restaurant franchise parametrization of the hierarchical
Dirichlet process to efficiently sample from the posterior.

###  4.2 Prior on direct effects

Similar to our prior on $\bm{L}$, we utilize a hierarchical Dirichlet process
mixture model to place a prior on $\bm{\Delta}$’s hyperparameters in Algorithm
1. However, it differs from the prior presented Section 4.1 in two important
ways. First, it only considers sub-pathways and ignores super-pathway labels.
While this does discard information, we found that incorporating super-
pathways led to trivial changes in estimates and was not worth the extra
computation. Second, our model for $\bm{\Delta}$’s hyperparameters induces a
clustering on metabolites as opposed to sub-pathways, which helps pool
information across related metabolites to improve otherwise uncertain
hyperparameter estimates.

We build $\bm{\Delta}$’s prior to resemble that in [23], which modeled genetic
effects on metabolites as a spike and normal-slab with metabolite-specific
spike probabilities and constant slab variance. We use the same super- and
sub-pathway notation to be consistent with Section 4.1 and re-index
$\bm{\Delta}$ for the purposes of Section 4.2 so that $\bm{\Delta}_{s,pbm}$ is
the direct effect of the $s$th SNP on metabolite $m\in[M_{pb}]$ in sub-pathway
$b\in[B_{p}]$ and super-pathway $p\in[P]$. Then

$\displaystyle\begin{aligned}
&F_{0}^{\operatorname{(\pi)}}\mid\\{\gamma^{\operatorname{(\pi)}},H^{\operatorname{(\pi)}}\\}\sim\operatorname{DP}(\gamma^{\operatorname{(\pi)}},H^{\operatorname{(\pi)}}),\quad\varphi^{2}\mid
H^{\operatorname{(\varphi^{2})}}\sim H^{\operatorname{(\varphi^{2})}}\\\
&\text{$F_{pb}^{\operatorname{(\pi)}}\mid\\{\alpha_{0}^{\operatorname{(\pi)}},F_{0}^{\operatorname{(\pi)}}\\}\sim\operatorname{DP}(\alpha_{0}^{\operatorname{(\pi)}},F_{0}^{\operatorname{(\pi)}}),\quad
p\in[P];b\in[B_{p}]$}\\\ &\text{$\pi_{pbm}\mid
F_{pb}^{\operatorname{(\pi)}}\sim F_{pb}^{\operatorname{(\pi)}},\quad
m\in[M_{pb}]$}\\\
&\text{$\bm{\Delta}_{s,pbm}\mid(\pi_{pbm},\varphi^{2},\bm{D})\sim(1-\pi_{pbm})\delta_{0}+\pi_{pbm}N\left(0,\varphi^{2}\sigma_{pbm}^{2}/\bm{D}_{ss}\right),\quad
s\in[S],$}\end{aligned}$ (4.3)

where $H^{(\pi)}$ and $H^{(\varphi^{2})}$ are probability measures on $[0,1]$
and $(0,\infty)$ and $\sigma_{pbm}^{2}$ is the re-indexed version of
$\sigma_{m}^{2}$ defined in (2.1). We replace it with its consistent estimator
defined in step 3 of Algorithm 1 in practice. Standardizing $\bm{\Delta}$’s
slab variance by $\sigma_{pbm}^{2}$ adjusts for heterogeneity in metabolite
levels, and dividing by $\bm{D}_{ss}$ reflects the fact that genetic effects
are often inversely proportional to minor allele frequency. Figure 3(b) gives
a graphical depiction of (4.3) where, just like we did in Section 4.1, we
place diffuse priors on the concentration parameters
$\gamma^{(\pi)},\alpha_{0}^{(\pi)}$.

The observed data $\hat{\bm{\Delta}}$ will be informative for $\varphi^{2}$
because it is shared across all SNPs and metabolites. However, $\pi_{pbm}$ is
metabolite-specific which, given $\bm{\Delta}$’s sparsity, suggests its
estimates may be uncertain. To address this, our hierarchical Dirichlet
process mixture prior on $\pi_{pbm}$ induces a clustering on metabolites,
where metabolites in the same cluster share $\pi_{pbm}$’s and are more likely
to be grouped with other metabolites in the same sub-pathway. This pools
information across related metabolites, thereby reducing our uncertainty in
$\pi_{pbm}$.

It remains to specify $H^{(\pi)}$ and $H^{(\varphi^{2})}$, where care must be
taken to ensure inference is computationally tractable. Let
$\Theta=\\{\\{\pi_{pbm}\\}_{p\in[P];b\in[B_{p}];m\in[M_{pb}]},\varphi^{2}\\}$.
Then $\bm{\Delta}$’s posterior is

$\displaystyle\operatorname{pr}\\{\bm{\Delta}\mid\hat{\bm{\Delta}},\bm{D},H^{(\pi)},H^{(\varphi^{2})}\\}=$
$\displaystyle\smallint\operatorname{pr}(\bm{\Delta}\mid\hat{\bm{\Delta}},\bm{D},\Theta)\operatorname{pr}\\{\Theta\mid\hat{\bm{\Delta}},\bm{D},H^{(\pi)},H^{(\varphi^{2})}\\}d\Theta.$
(4.4)

The first posterior in the integral matches that in step 4 of Algorithm 1 and
is an entry-wise spike and normal-slab distribution with a closed form.
However, $\Theta$’s posterior is generally intractable, meaning the integral
must be approximated via Markov chain Monte Carlo. Unfortunately, a Gibbs
sampler with standard conjugate inverse-gamma and beta distributions for
$H^{(\varphi^{2})}$ and $H^{(\pi)}$ would require sampling $\bm{\Delta}$ and
would result in intractably slow mixing due $\bm{\Delta}$’s high dimension. We
address this by letting $H^{(\pi)}$ and $H^{(\varphi^{2})}$ be discrete.
Briefly, let $\\{0\\}\cup\mathcal{X}^{(\pi)}\subset[0,1)$ and
$\mathcal{X}^{(\varphi^{2})}\subset(0,\infty)$ be finite sets containing all
possible values of $\pi_{pbm}$ and $\varphi^{2}$. We then set $H^{(\pi)}$ and
$H^{(\varphi^{2})}$ so that

$\displaystyle\mathbb{P}\\{\pi_{pbm}=0\mid
H^{(\pi)}\\}=0.5,\quad\text{$\mathbb{P}\\{\pi_{pbm}=x\mid
H^{(\pi)}\\}=0.5/\lvert\mathcal{X}^{(\pi)}\rvert$ for all
$x\in\mathcal{X}^{(\pi)}$}$
$\displaystyle\text{$\mathbb{P}\\{\varphi^{2}=x\mid
H^{(\varphi^{2})}\\}=1/\lvert\mathcal{X}^{(\varphi^{2})}\rvert$ for all
$x\in\mathcal{X}^{(\varphi^{2})}$},$

where $H^{(\pi)}$’s non-trivial mass at 0 reflects our expectation that some
metabolites will have no direct effects. Our software-default of 0.5 is likely
conservative given that the levels of many metabolites are impacted by at
least one SNP [7]. Critically, the above discrete distribution facilitates
faster mixing in our Gibbs sampler to approximate (4.4) by allowing us to
integrate out $\bm{\Delta}$, and lets us pre-compute the likelihood
$\operatorname{pr}\\{\hat{\bm{\Delta}}\mid\Theta,\bm{D},H^{(\pi)},H^{(\varphi^{2})}\\}=\operatorname{pr}(\hat{\bm{\Delta}}\mid\Theta,\bm{D})$
for all $\Theta$ required to propose new states, which significantly reduces
computation time. Section A3 details how we set $\mathcal{X}^{(\pi)}$ and
$\mathcal{X}^{(\varphi^{2})}$.

###  4.3 Pathway-specific inference on direct effects

Our prior on direct effects in (4.3) lets us compare the “direct” heritability
of metabolic sub-pathways. To see this, let $\Theta$ be as defined in Section
4.2 and $\bm{Y}_{i,pbm}$ be individual $i$’s abundance of metabolite $m$ in
sub- and super-pathways $b$ and $p$. Then under model (2.1) and
$\bm{\Delta}$’s prior in (4.3), the heritability of metabolite $m$ conditional
on latent factors $\bm{C}$ is

$\displaystyle
h_{pbm}^{2}:=\frac{N^{-1}\sum_{i=1}^{N}\operatorname{Var}(\textstyle\sum_{s=1}^{S}\bm{X}_{is}\bm{\Delta}_{s,pbm}\mid\bm{L},\bm{C},\Theta)}{N^{-1}\sum_{i=1}^{N}\operatorname{Var}(\bm{Y}_{i,pbm}\mid\bm{L},\bm{C},\Theta)}=\frac{\varphi^{2}\pi_{pbm}}{\varphi^{2}\pi_{pbm}+1},$

where the numerator and denominator in the first fraction are the average
variance of the metabolite’s direct genetic effect and abundance,
respectively. Since a metabolite’s heritability is an increasing function of
$\pi_{pbm}$, evaluating sub-pathway $b$ in super-pathway $p$’s heritability is
equivalent to studying the distribution of $\\{\pi_{pbm}\\}_{m\in[M_{pb}]}$.
This distribution is exactly $F_{pb}^{(\pi)}$ defined in (4.3), where a sub-
pathway is more heritable when draws from $F_{pb}^{(\pi)}$ tend to be large.
We therefore define sub-pathway $b$ in super-pathway $p$’s direct heritability
score, $DHS_{pb}$, to be

$\displaystyle
DHS_{pb}=\mathbb{P}\\{f_{pb}^{(\pi)}>f_{0}^{(\pi)}\mid\hat{\bm{\Delta}},\bm{D},H^{(\pi)},H^{(\varphi^{2})}\\},\quad
f_{pb}^{(\pi)}\mid F_{pb}^{(\pi)}\sim F_{pb}^{(\pi)},\quad f_{0}^{(\pi)}\mid
F_{0}^{(\pi)}\sim F_{0}^{(\pi)}.$ (4.5)

Since $F_{0}^{(\pi)}=\mathbb{E}\\{F_{pb}^{(\pi)}\mid F_{0}^{(\pi)}\\}$,
$DHS_{pb}$ is interpretable as the probability a randomly chosen metabolite
from sub-pathway $b$ is more heritable than one from the average sub-pathway.
These scores allow us to rank sub-pathways in terms of their heritability,
where sub-pathways with larger $DHS_{pb}$ are more heritable. We demonstrate
their utility in Section 6.

##  5 Simulations

###  5.1 Simulation setup

We simulated metabolomic data from $N=500$ individuals whose $M=257$
metabolites were regulated by $S=5,000$ SNPs to evaluate our estimator for the
number of latent factors and Algorithm 1. Unobserved genotypes $\bm{X}$ and
metabolite levels $\bm{Y}$ were generated assuming $K=10$ latent factors:

$\displaystyle\begin{aligned} &f_{s}\sim
U[0.05,0.5],\quad\bm{X}_{is}\sim\text{Bin}(2,f_{s}),\quad i\in[N];s\in[S]\\\
&\bm{G}_{sk}\sim(1-10^{-3})\delta_{0}+10^{-3}N(0,2/N),\quad\bm{L}_{*k}\sim
HDP(\lambda_{k}),\quad s\in[S];k\in[K]\\\
&\bm{\Delta}_{sm}\sim(1-\pi_{\Delta_{m}})\delta_{0}+\pi_{\Delta_{m}}N(0,\tau_{\Delta
s}^{2}),\quad\sigma_{m}^{2}\sim\text{Gamma}(1,1),\quad s\in[S];m\in[M]\\\
&\bm{C}\sim MN(\bm{X}\bm{G},I_{N},I_{K}),\quad\bm{Y}\sim
MN(\bm{C}\bm{L}^{\top}+\bm{X}\bm{\Delta},\operatorname{diag}(\sigma_{1}^{2},\ldots,\sigma_{M}^{2}),I_{N}).\end{aligned}$
(5.1)

For each dataset we defined GWAS summary statistics $\hat{\beta}_{sm}$ to be
estimates for the slope in the regression of $\bm{Y}_{*m}$ onto $\bm{X}_{*s}$
and the observed matrix $\hat{\bm{B}}$ such that
$\hat{\bm{B}}_{sm}=(\bm{X}_{*s}^{\top}\bm{P}_{\bm{1}_{N}}^{\perp}\bm{X}_{*s})^{1/2}\hat{\beta}_{sm}$
for $\bm{P}_{\bm{1}_{N}}^{\perp}$ the orthogonal projection matrix that
projects vectors onto the kernel of $(1,\ldots,1)\in\mathbb{R}^{1\times N}$.

The parameter $\pi_{\Delta_{m}}$ was simulated as $\pi_{\Delta_{m}}\sim
0.3\delta_{0}+0.7F$ for $F$ a truncated exponential distribution between 0 and
0.01 with mean $0.005$, and $\tau^{2}_{\Delta
s}=5^{2}/(\bm{X}_{*s}^{\top}\bm{P}_{\bm{1}_{N}}^{\perp}\bm{X}_{*s})$. These,
along with $\bm{G}$’s hyperparameters, were all commensurate to the values
estimated in our real data example. We drew $\bm{L}$’s columns independently
from a hierarchical Dirichlet process (HDP) mixture model by partitioning
metabolites into super- and sub-pathways such that the number of constituent
metabolites in each partition mirrored our real data example (see Section A4
for details). The parameters
$\lambda_{k}=M^{-1}\mathbb{E}(\bm{L}_{*k}^{\top}\bm{L}_{*k})$ dictate the
strength of each latent factor where, as we discussed in Section 3.4, factors
with larger $\lambda_{k}$ are stronger and easier to recover. To make
simulations realistic, we let $\lambda_{1}>\cdots>\lambda_{K}$ take a wide
range of values, where $\lambda_{1}=3.35$ corresponded to a large pervasive
factor and $\lambda_{K}=0.25$ was weaker. Notably, these values implied there
was a small eigengap between the signal and noise, since the $K$th singular
value of $\hat{\bm{B}}$ was on the average only 14.5% larger than
$\lVert\hat{\bm{B}}-\mathbb{E}(\hat{\bm{B}}\mid\bm{X},\bm{C},\bm{L},\bm{\Delta})\rVert_{2}$,
the size of the noise. Consequently, as is the case in real data, recovering
all factors was non-trivial.

###  5.2 Simulation results

We first considered estimating the number of latent factors $K$. We compared
our method described in Section 3.3 (dBEMA) to bulk eigenvalue matching (BEMA)
proposed in [34], the eigenvalue difference approach (ED) taken in [38], the
information criterion (PANICr) from [36], parallel analysis (PA) [46], and bi-
cross validation (BCV) [31]. Figure 4(a) contains the results, where our
method clearly outperforms existing methods. Notably, PA and BCV, which are
ubiquitous in omic analysis pipelines, severely overestimate $K$. This is
because the error matrix $\tilde{\bm{E}}$ in the expression for $\hat{\bm{B}}$
in (2.2) has dependent rows (SNPs), which belies their assumption that the
error matrix has independent entries. As a consequence, they mistake the
error’s row-wise dependence for latent structure.

Figure 4: Simulation results for 50 simulated datasets. (a): Estimates for
$K$. dBEMA is our proposed method and BCV’s maximum $K$ was set to 20.
Percents give the fraction of times the estimate matched the true $K=10$. (b):
Inference on $\bm{G}$. Local false sign rates were binned into 10 bins. (c):
Inference on $\bm{\Delta}$. Each point represents a simulated dataset.

Assuming $K$ was known, we next considered estimating $\bm{G}$ and $\bm{L}$,
which parameterize indirect effects, and direct effects $\bm{\Delta}$. We
first used Algorithm 1 and the priors described in Sections 4.1 to estimate
$\bm{G}$ and $\bm{L}$. To address issues of factor identifiability and match
estimated and simulated factors, we rotated our estimators $\hat{\bm{G}}$ and
$\hat{\bm{L}}$ defined in Algorithm 1 by solving Procrustes problem
$\text{argmin}_{\bm{Q}^{\top}\bm{Q}=I_{K}}\lVert\hat{\bm{G}}\bm{Q}-\tilde{\bm{G}}\rVert_{F}$
for $\tilde{\bm{G}}$ defined in (2.2). We used the local false sign rate
($\operatorname{lfsr}$) to perform inference on $\bm{G}$:

$\displaystyle\operatorname{lfsr}(\bm{G}_{sk})=1-\max\\{\mathbb{P}(\bm{G}_{sk}>0\mid\hat{\bm{B}},\bm{D}),\mathbb{P}(\bm{G}_{sk}<0\mid\hat{\bm{B}},\bm{D})\\},\quad
s\in[S];k\in[K].$ (5.2)

The $\operatorname{lfsr}$ is the probability a coefficient’s inferred non-zero
sign is incorrect and is at least as large as the local false discovery rate.
It is used here to avoid issues with the local false discovery rate that can
arise when the null and alternative distributions are hard to distinguish
[26]. Figure 4(b) shows that our estimate for $\bm{G}$’s $\operatorname{lfsr}$
controls the local false discovery rate. Section A4 contains results for
$\bm{L}$, which illustrate the importance of our novel theoretical results
from Section 3.4 giving $\hat{\bm{L}}$’s asymptotic distribution.

We lastly estimated $\bm{\Delta}$. Since the prior for $\bm{\Delta}$
recommended for metabolomic data in Section 4.2 requires non-trivial
computation, we perform inference in these simulated data by controlling the
false discovery rate via the Benjanmini-Hochberg procedure [47] after using
$\hat{\bm{\Delta}}$’s likelihood in Algorithm 1 to compute p-values for the
null hypotheses $H_{0,sm}:\bm{\Delta}_{sm}=0$. Figure 4(c) shows this controls
the false discovery rate, indicating $\hat{\bm{\Delta}}$’s likelihood in
Algorithm 1 is appropriate.

##  6 Real data analysis

We demonstrate the power of our methodology using metabolite summary
statistics derived from $N=6,136$ Finnish adults participating in the METSIM
study [4]. Briefly, we considered SNPs with minor allele frequencies $\geq
5\%$ and subsequently pruned them for linkage disequilibrium. Metabolites
without a name or missing in more than 20% of samples were excluded, which
resulted in $S=70,140$ SNPs and $M=765$ metabolites. Note that while some
metabolites had missing data, the missingness was likely inconsequential, as
more that 80% of metabolites had less than 5% missing data. Figure 2(a)
provides an overview of the hierarchical partition of metabolites into super-
and sub-pathways.

###  6.1 Indirect effects

We first estimated the number of latent factors $K$ and indirect effects.
Table 1(a) gives each method’s estimate for $K$, where our method dBEMA
estimates 36. Compared to dBEMA, the methods ED, PANICr, and BCV behave as
they did in simulations, where ED likely underestimates and PANICr and BCV
likely overestimate $K$. We then used Algorithm 1 with $\hat{K}=36$ latent
factors to determine $\bm{G}$ and $\bm{L}$’s posterior using Section 4.1’s
priors, and determined the SNPs and metabolites involved in indirect effects
using the local false sign rate (lfsr) defined in (5.2), where
$\operatorname{lfsr}(\bm{L}_{mk})$ was defined by replacing $\bm{G}_{sk}$ with
$\bm{L}_{mk}$ in (5.2). We defined a “genetic factor” to be a factor $k$
satisfying $\operatorname{lfsr}(\bm{G}_{sk})\leq 0.05$ for at least one SNP
$s\in[S]$, an “indirect SNP” to be a SNP $s$ with
$\operatorname{lfsr}(\bm{G}_{sk})\leq 0.05$ for at least one $k\in[\hat{K}]$,
and an “indirect metabolite” to be a metabolite $m$ satisfying
$\operatorname{lfsr}(\bm{L}_{mk})\leq 0.05$ for at least one genetic factor
$k$. Table 1(b) provides the results, where all but factors 1, 2, 3, and 6 are
genetic factors. As these factors have the largest loadings, it suggests
large, pervasive factors reflect uninteresting variation and reinforces the
importance of non-pervasive factors.

Table 1: (a): Estimated number of factors. “N/A” means the method’s software
could not be run on these data; the maximum number of factors considered by
BCV was 100. (b): Indirect effect results.

Table 1(b) implies the levels of all metabolites are regulated by a small
number of SNPs, which suggests these SNPs are highly pleiotropic. Our results
are substantiated by [7], who studied 233 serum metabolites measured in
136,016 individuals and found that all metabolites were regulated by 276
genomic regions. Remarkably, these regions contained 306 (73%) of the indirect
SNPs we identified, which is more than four times more than expected by chance
alone (p-value $<10^{-16}$).

Since all metabolites appear to be indirectly genetically regulated, simply
reporting indirect metabolites will do nothing to aid biological
interpretation. Instead, we interpret factors mediating genetic effects by
studying how their loadings behave across metabolic super- and sub-pathways
using $\bm{L}$’s posterior. As discussed in Section 4.1, $\bm{L}$’s posterior
begets a clustering on sub-pathways that couples biological processes to
improve inference. We illustrate the utility of this approach using the first
three genetic factors whose estimated metabolite and SNP loadings are plotted
in the last three rows of Figure 5(a). For each factor $k\in[\hat{K}]$, we
used our Gibbs sampler detailed in Section A2 to sample from
$\operatorname{pr}[\\{(\mu_{pb},\phi_{pb}^{2})\\}_{p\in[P];b\in[B_{p}]}\mid\hat{\bm{B}},\bm{D}]$,
where parameters $(\mu_{pb},\phi_{pb}^{2})$ are defined in (4.2) and
implicitly depend on factor $k$. These parameters characterize the
distribution of factor $k$’s loadings for metabolites in sub-pathway $b$ and
super-pathway $p$, where a “cluster” is a set of sub-pathways that share
parameters. We defined the maximum a posteriori clustering to be our hard
clustering, which we found gave the most parsimonious results compared to
other hard clustering methods. We determined a cluster $\mathcal{C}$’s
importance with enrichment ($\mathscr{E}$) and sign ($\mathscr{S}$) scores:

$\displaystyle\mathscr{E}_{\mathcal{C}}=\lvert\mathcal{C}\rvert^{-1}\textstyle\sum_{(p,b)\in\mathcal{C}}\mathbb{E}(\mu_{pb}^{2}+\phi_{pb}^{2}\mid\hat{\bm{B}},\bm{D}),\quad\mathscr{S}_{\mathcal{C}}=$
$\displaystyle\max\textstyle[2\lvert\mathcal{C}\rvert^{-1}\sum_{(p,b)\in\mathcal{C}}\mathbb{P}\\{N(\mu_{pb},\phi_{pb}^{2})>0\mid\hat{\bm{B}},\bm{D}\\}$
$\displaystyle\textstyle
2\lvert\mathcal{C}\rvert^{-1}\sum_{(p,b)\in\mathcal{C}}\mathbb{P}\\{N(\mu_{pb},\phi_{pb}^{2})<0\mid\hat{\bm{B}},\bm{D}\\}]-1.$

The enrichment score is the average squared magnitude of the cluster’s
loadings. The sign score is bounded between 0 and 1 and will be close to 1 if
the cluster’s loading distribution is far away from zero and 0 if it is
symmetric around zero. Clusters with large scores have large loadings with the
same non-zero signs, and are therefore ideal candidates to help infer a
factor’s function.

Figure 5: (a): Metabolite and SNP loadings for factor 1 and the first three
genetic factors. Dashed red lines demarcate super-pathways whose names are
given in Figure 2(a) and each set of colored metabolites is a sub-pathway;
‘$\triangledown$’ points indicate outlier metabolites in that factor, defined
as metabolites whose posterior probability of being an outlier was $\geq
0.95$. (b): Enrichment and sign scores for clusters identified in the three
genetic factors, where the columns of $\hat{\bm{L}}$ were scaled to all have
the same norm prior to running the Gibbs sampler to make scores comparable
across factors. Scores are negative if more of the cluster’s estimated
loadings are negative than positive.

Figure 5(b) gives a compendious overview of the clustering results for the
first three genetic factors, where the pie charts give memberships for
clusters with the largest enrichment scores. We only display one cluster for
factor 5 because its second highest enrichment and sign scores were small.
Remarkably, our unsupervised clusters mirror biology. Cluster 1 in factors 5
and 7 contain both sphingomyelin sub-pathways and all three classes of sex
hormones, respectively. Grouping sphingomyelins with ceramides in factors 4
and 7 is expected given that a ubiquitous reaction produces sphingomyelins
from ceramides, where the inclusion of diacylglycerols in cluster 14 is likely
because the reaction creates them as a by-product [48].

We use these clusters and their corresponding significant SNPs (the red dots
in Figure 5(b)) to infer indirect mechanisms by which genotype systematically
regulates metabolite levels.

Figure 6: Interpretations of factors 4 and 7.

Factor 4: Glucose homeostasis. Factor 4 contains four significant SNPs that
all map to the gene GCKR, which binds to and inactivates the enzyme
glucokinase in the liver to inhibit hepatic glycolysis and subsequent
lipogenesis [49]. Since the four SNPs’ minor alleles reduce GCKR’s ability to
bind to glucokinase [49], we posit factor 4 relates to glucose homeostasis
characterized by an up-regulation of glycolysis and lipogenesis in the liver
(Figure 6). Excess glycolysis has the effect of increasing the levels of
tyrosine and its derivatives, which include microbiome-derived benzoates [50],
as well as possibly increasing the activity of the enzyme xanthine oxidase
that is responsible for many metabolites in the xanthine sub-pathway [51, 52].
Hepatic lipogenesis raises the levels of sphingomyelin and lactosylceramide
metabolites [53].

Factor 7: Liver function. 43 of the 54 SNPs regulating factor 7 could be
mapped to one of the three genes in Figure 6. SLCO1B1 is uniquely expressed in
the liver and encodes the uptake transport protein responsible for
transporting sex hormones and their related metabolites [54]. The SLCO1B1
polymorphisms we identified tend to increase its activity [55], which likely
explains the decrease in serum sex hormone levels. FADS1 encodes the rate-
limiting enzyme in the synthesis of long chain polyunsaturated fatty acids
(LC-PUFAs), where the polymorphisms we identify reduce its expression to beget
hepatic lipid accumulation and subsequent sphingomyelin synthesis [56, 57].
Less is known about ALMS1, although recent work suggests it contributes to
nonalcoholic fatty liver disease [58].

Factor 5: Hepatic lipid metabolism. Factor 5’s 19 significant SNPs map to
either FADS1/2/3 or SLCO1B1. Unlike factors 4 and 7, factor 5 is characterize
by a substantial number of outlier metabolites with similar large negative
loadings (Figure 5(a)), suggesting these metabolites have related functions
and can be used to help infer factor 5’s function. Remarkably, 24 out of the
25 outliers with loadings $<-0.4$ were omega-3 LC-PUFAs, which is congruent
with the observation that the FADS polymorphisms we identify reduce the
production of these anti-inflammatory LC-PUFAs in the liver [59]. Relatedly,
SLCO1B1’s polymorphisms are potential markers for nonalcoholic fatty liver
disease [60, 61].

###  6.2 Direct effects

We lastly used Algorithm 1 and our pathway-guided prior defined in 4.2 to
infer $\bm{\Delta}$. Figure 7(a) contains an overview of our results. We first
compare our results from HiGSS to those from a standard mtGWAS analysis, which
uses the entries of $\hat{\bm{B}}$ and their standard errors to test
$H_{0,sm}:\bm{B}_{sm}=\bm{G}_{s*}^{\top}\bm{L}_{m*}+\bm{\Delta}_{sm}=0$ and
subsequently controls the family-wise error rate with a Bonferroni correction.
We identify 49% more metabolite-SNP pairs than the standard analysis, which is
quite remarkable given that we are only inferring direct effects whereas the
standard analysis will reject pairs with significant direct or indirect
effects.

Figure 7: Inference on $\bm{\Delta}$. (a): The FSR is the average local false
sign rate. (b): Mean-centered direct heritability scores for the 10 least and
most heritable sub-pathways.

While these results highlight the power of our method, interpreting them
remains a challenge given the number of significant associations. We therefore
considered performing a more interpretable analysis at the pathway level to
determine the sub-pathways that have the largest direct effects. To do so, we
used sub-pathway direct heritability scores defined in (4.5), where a large
score indicates metabolites in that sub-pathway are more heritable. Figure
7(b) contains the results for the 10 least and most heritable sub-pathways.
Unsurprising, xenobiotic sub-pathways, which typically contain exogenous
metabolites, are the least heritable. The most heritable sub-pathway is
“Androgenic Steroids”, which is congruent with existing results from twin
studies suggesting the heritability of circulating androgens may be as high as
50% [62].

##  7 Conclusion

In this work we developed a model, theoretical framework, and set of methods
to perform GWAS of high dimensional phenotypes using summary statistics that
explicitly model pleiotropy, beget fast computation, and facilitate the use of
biologically informed priors. Given the importance of mtGWAS, we spent
considerable effort developing hierarchical Dirichlet priors for indirect and
direct genetic effects that foster interpretable inference at the metabolic
pathway level.

Our real data analysis in Section 6 demonstrates the utility of our approach,
where we were able to use our hierarchical Dirichlet process prior for
indirect effects to couple metabolite sub-pathways by clustering them into
biologically meaningful groups. We showed these were useful for interpreting
factors, thereby providing putative mechanisms by which genotype impacts
metabolite levels. Results using our method also suggested all observed
metabolites were genetically regulated, which recapitulated results from a
study with over 20 times as many samples as our data example. These
conclusions were not possible using the standard analysis, which could only
identify genetic associations for 67% of metabolites and highlights the power
of our approach.

## References

  * [1] Emil Uffelmann, Qin Qin Huang, Nchangwi Syntia Munung, Jantina Vries, Yukinori Okada, Alicia R. Martin, Hilary C. Martin, Tuuli Lappalainen and Danielle Posthuma  “Genome-wide association studies”  In _Nature Reviews Methods Primers_ 1.1, 2021, pp. 59
  * [2] Yosuke Tanigawa et al.  “Components of genetic associations across 2,138 phenotypes in the UK Biobank highlight adipocyte biology”  In _Nature Communications_ 10.1, 2019, pp. 4064
  * [3] Lloyd T. Elliott, Kevin Sharp, Fidel Alfaro-Almagro, Sinan Shi, Karla L. Miller, Gwenaëlle Douaud, Jonathan Marchini and Stephen M. Smith  “Genome-wide association studies of brain imaging phenotypes in UK Biobank”  In _Nature_ 562.7726, 2018, pp. 210–216
  * [4] Xianyong Yin et al.  “Genome-wide association studies of metabolites in Finnish men identify disease-relevant loci”  In _Nature Communications_ 13.1, 2022, pp. 1644
  * [5] Alexander Kurilshikov et al.  “Large-scale association analyses identify host factors influencing human gut microbiome composition”  In _Nature Genetics_ 53.2, 2021, pp. 156–165
  * [6] Xianyong Yin et al.  “Integrating transcriptomics, metabolomics, and GWAS helps reveal molecular mechanisms for metabolite levels and disease risk”  In _The American Journal of Human Genetics_ 109.10 Elsevier BV, 2022, pp. 1727–1741  DOI: 10.1016/j.ajhg.2022.08.007
  * [7] Minna K. Karjalainen et al.  “Genome-wide characterization of circulating metabolic biomarkers reveals substantial pleiotropy and novel disease pathways”  Cold Spring Harbor Laboratory, 2022  DOI: 10.1101/2022.10.20.22281089
  * [8] Chris McKennan, Carole Ober and Dan Nicolae  “Estimation and inference in metabolomics with nonrandom missing data and latent factors”  In _The Annals of Applied Statistics_ 14.2, 2020  DOI: 10.1214/20-aoas1328
  * [9] Guanghao Qi and Nilanjan Chatterjee  “Heritability informed power optimization (HIPO) leads to enhanced detection of genetic associations across multiple traits”  In _PLoS Genetics_ 14.10, 2018, pp. e1007549
  * [10] Mingwei Dai, Xiang Wan, Hao Peng, Yao Wang, Yue Liu, Jin Liu, Zongben Xu and Can Yang  “Joint analysis of individual-level and summary-level GWAS data by leveraging pleiotropy”  In _Bioinformatics_ 35.10, 2019, pp. 1729–1736
  * [11] Jin Liu, Xiang Wan, Chaolong Wang, Chao Yang, Xiaowei Zhou and Can Yang  “LLR: a latent low-rank approach to colocalizing genetic risk variants in multiple GWAS”  In _Bioinformatics_ 33, 2017, pp. 3878–3886
  * [12] Yihao Lu, Meritxell Oliva, Brandon Pierce, Jin Liu and Lin Chen  “Integrative cross-omics and cross-context analysis elucidates molecular links underlying genetic effects on complex traits”, 2022  DOI: 10.21203/rs.3.rs-2137179/v1
  * [13] Chong Wu  “Multi-trait Genome-Wide Analyses of the Brain Imaging Phenotypes in UK Biobank”  In _Genetics_ 215.4 Oxford University Press (OUP), 2020, pp. 947–958  DOI: 10.1534/genetics.120.303242
  * [14] Meida Wang, Xuewei Cao, Shuanglin Zhang and Qiuying Sha  “A clustering linear combination method for multiple phenotype association studies based on GWAS summary statistics”  In _Scientific Reports_ 13.1, 2023, pp. 3389
  * [15] Samantha Petti, Gautam Reddy and Michael M. Desai  “Inferring sparse structure in genotype-phenotype maps”  In _bioRxiv_ Cold Spring Harbor Laboratory, 2022  DOI: 10.1101/2022.09.27.509675
  * [16] Andrew Brantley Hall, Andrew C. Tolonen and Ramnik J. Xavier  “Human genetic variation and the gut microbiome in disease”  In _Nature Reviews Genetics_ 18.11, 2017, pp. 690–699
  * [17] Liangliang Zhang, Yushu Shi, Robert R. Jenq, Kim-Anh Do and Christine B. Peterson  “Bayesian compositional regression with structured priors for microbiome feature selection”  In _Biometrics_ 77.3, 2021, pp. 824–838
  * [18] Per Sidén, Anders Eklund, David Bolin and Mattias Villani  “Fast Bayesian whole-brain fMRI analysis with spatial 3D priors”  In _NeuroImage_ 146, 2017, pp. 211–225
  * [19] Kris Sankaran and Susan P Holmes  “Latent variable modeling for the microbiome”  In _Biostatistics_ 20.4, 2018, pp. 599–614  DOI: 10.1093/biostatistics/kxy018
  * [20] Yan Zhang, Guanghao Qi, Ju-Hyun Park and Nilanjan Chatterjee  “Estimation of complex effect-size distributions using summary-level statistics from genome-wide association studies across 32 complex traits”  In _Nature Genetics_ 50.9, 2018, pp. 1318–1326
  * [21] Xiang Zhou, Peter Carbonetto and Matthew Stephens  “Polygenic Modeling with Bayesian Sparse Linear Mixed Models”  In _PLoS Genetics_ 9.2, 2013, pp. e1003264–
  * [22] Xiang Zhu and Matthew Stephens  “Bayesian large-scale multiple regression with summary statistics from genome-wide association studies”  In _The Annals of Applied Statistics_ 11.3 NIH Public Access, 2017, pp. 1561
  * [23] Helene Ruffieux, Anthony C. Davison, Jorg Hager and Irina Irincheeva  “Efficient inference for genetic association studies with multiple outcomes”  In _Biostatistics_ 18.4, 2017, pp. 618–636  DOI: 10.1093/biostatistics/kxx007
  * [24] Patrick Turley et al.  “Multi-trait analysis of genome-wide association summary statistics using MTAG”  In _Nature Genetics_ 50.2, 2018, pp. 229–237
  * [25] Arunabha Majumdar, Tanushree Haldar, Sourabh Bhattacharya and John S. Witte  “An efficient Bayesian meta-analysis approach for studying cross-phenotype genetic associations”  In _PLoS Genetics_ 14.2, 2018, pp. e1007139–
  * [26] Matthew Stephens  “False discovery rates: a new deal”  In _Biostatistics_ 18.2 Oxford University Press, 2017, pp. 275–294
  * [27] Adam Auton et al.  “A global reference for human genetic variation”  In _Nature_ 526.7571, 2015, pp. 68–74
  * [28] Chris McKennan and Dan Nicolae  “Accounting for unobserved covariates with varying degrees of estimability in high-dimensional biological data”  In _Biometrika_ 106.4, 2019, pp. 823–840  DOI: 10.1093/biomet/asz037
  * [29] Xinyi Zhong, Chang Su and Zhou Fan  “Empirical Bayes PCA in high dimensions”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 84, 2020, pp. 853–878
  * [30] Wei Wang and Matthew Stephens  “Empirical Bayes Matrix Factorization”  In _Journal of Machine Learning Research_ 22.120, 2021, pp. 1–40  URL: http://jmlr.org/papers/v22/20-589.html
  * [31] Art B Owen and Jingshu Wang  “Bi-cross-validation for factor analysis”  In _Statistical Science_ 31.1 Institute of Mathematical Statistics, 2016, pp. 119–139
  * [32] Edgar Dobriban and Art B. Owen  “Deterministic parallel analysis: an improved method for selecting factors and principal components”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 81.1 Oxford University Press (OUP), 2018, pp. 163–183  DOI: 10.1111/rssb.12301
  * [33] Edgar Dobriban  “Permutation methods for factor analysis and PCA”  In _The Annals of Statistics_ 48.5 Institute of Mathematical Statistics, 2020  DOI: 10.1214/19-aos1907
  * [34] Zheng Tracy Ke, Yucong Ma and Xihong Lin  “Estimation of the number of spiked eigenvalues in a covariance matrix by bulk eigenvalue matching analysis”  In _Journal of the American Statistical Association_ Taylor & Francis, 2021, pp. 1–19
  * [35] Chris McKennan and Dan Nicolae  “Estimating and Accounting for Unobserved Covariates in High-Dimensional Correlated Data”  In _Journal of the American Statistical Association_ 117.537 Taylor & Francis, 2022, pp. 225–236
  * [36] Jushan Bai and Serena Ng  “Determining the number of factors in approximate factor models”  In _Econometrica_ 70.1 Wiley Online Library, 2002, pp. 191–221
  * [37] Seung C. Ahn and Alex R. Horenstein  “Eigenvalue Ratio Test for the Number of Factors”  In _Econometrica_ 81.3, 2013, pp. 1203–1227
  * [38] Alexei Onatski  “Determining the number of factors from empirical distribution of eigenvalues”  In _The Review of Economics and Statistics_ 92.4 The MIT Press, 2010, pp. 1004–1016
  * [39] Chris McKennan  “Factor analysis in high dimensional biological data with dependent observations”, 2020  eprint: arXiv:2009.11134
  * [40] V A Marčenko and L A Pastur  “Distribution of eigenvalues for some sets of random matrices”  In _Mathematics of the USSR-Sbornik_ 1.4 IOP Publishing, 1967, pp. 457–483  DOI: 10.1070/sm1967v001n04abeh001994
  * [41] Jianqing Fan, Kaizheng Wang, Yiqiao Zhong and Ziwei Zhu  “Robust High-Dimensional Factor Models with Applications to Statistical Machine Learning”  In _Statistical Science_ 36.2 Institute of Mathematical Statistics, 2021  DOI: 10.1214/20- sts785
  * [42] Jingshu Wang, Qingyuan Zhao, Trevor Hastie and Art B. Owen  “Confounder adjustment in multiple hypothesis testing”  In _The Annals of Statistics_ 45.5, 2017, pp. 1863–1894
  * [43] Ju-Hyun Park, Mitchell H. Gail, Clarice R. Weinberg, Raymond J. Carroll, Charles C. Chung, Zhaoming Wang, Stephen J. Chanock, Joseph F. Fraumeni and Nilanjan Chatterjee  “Distribution of allele frequencies and effect sizes and their interrelationships for common genetic susceptibility variants”  In _Proceedings of the National Academy of Sciences_ 108.44, 2011, pp. 18026–18031
  * [44] Yee Whye Teh, Michael I Jordan, Matthew J Beal and David M Blei  “Hierarchical Dirichlet Processes”  In _Journal of the American Statistical Association_ 101.476, 2006, pp. 1566–1581
  * [45] Carlos M. Carvalho, Nicholas G. Polson and James G. Scott  “Handling sparsity via the horseshoe”  In _Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics_ 5, Proceedings of Machine Learning Research, 2009, pp. 73–80
  * [46] Jeffrey T Leek and John D Storey  “Capturing Heterogeneity in Gene Expression Studies by Surrogate Variable Analysis”  In _PLOS Genetics_ 3.9, 2007, pp. e161–
  * [47] Yoav Benjamini and Yosef Hochberg  “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing”  In _Journal of the Royal Statistical Society: Series B_ 57.1 Wiley, 1995, pp. 289–300  DOI: 10.1111/j.2517-6161.1995.tb02031.x
  * [48] Maristella Villani, Marimuthu Subathra, Yeong-Bin Im, Young Choi, Paola Signorelli, Maurizio Del Poeta and Chiara Luberto  “Sphingomyelin synthases regulate production of diacylglycerol at the Golgi”  In _Biochemical Journal_ 414.1, 2008, pp. 31–41  DOI: 10.1042/bj20071240
  * [49] Martijn C… Brouwers, Chantal Jacobs, Aalt Bast, Coen D.. Stehouwer and Nicolaas C. Schaper  “Modulation of Glucokinase Regulatory Protein: A Double-Edged Sword?”  In _Trends in Molecular Medicine_ 21.10, 2015, pp. 583–594
  * [50] Christopher P. Hopper, Ladie Kimberly De La Cruz, Kristin V. Lyles, Lauren K. Wareham, Jack A. Gilbert, Zehava Eichenbaum, Marcin Magierowski, Robert K. Poole, Jakob Wollborn and Binghe Wang  “Role of Carbon Monoxide in Host–Gut Microbiome Communication”  In _Chemical Reviews_ 120.24, 2020, pp. 13273–13311
  * [51] Christian Hellmuth, Franca Fabiana Kirchberg, Nina Lass, Ulrike Harder, Wolfgang Peissner, Berthold Koletzko and Thomas Reinehr  “Tyrosine Is Associated with Insulin Resistance in Longitudinal Metabolomic Profiling of Obese Children”  In _Journal of Diabetes Research_ 2016, 2016, pp. 1–10  DOI: 10.1155/2016/2108909
  * [52] Maria Elena Hernandez-Hernandez, Enrique Torres-Rasgado, Patricia Pulido-Perez, Leticia Nicolás-Toledo, Margarita Martínez-Gómez, Jorge Rodríguez-Antolín, Ricardo Pérez-Fuentes and Jose R. Romero  “Disordered Glucose Levels Are Associated with Xanthine Oxidase Activity in Overweight Type 2 Diabetic Women”  In _International Journal of Molecular Sciences_ 23.19, 2022, pp. 11177  DOI: 10.3390/ijms231911177
  * [53] Shireesha Sankella, Abhimanyu Garg and Anil K. Agarwal  “Activation of Sphingolipid Pathway in the Livers of Lipodystrophic Agpat2-/- Mice”  In _Journal of the Endocrine Society_ 1.7 The Endocrine Society, 2017, pp. 980–993  DOI: 10.1210/js.2017-00157
  * [54] Eunpi Cho, R. Montgomery and Elahe A. Mostaghel  “Minireview: SLCO and ABC Transporters: A Role for Steroid Transport in Prostate Cancer Progression”  In _Endocrinology_ 155.11, 2014, pp. 4124–4132  DOI: 10.1210/en.2014-1337
  * [55] Xiangyu Wu et al.  “Associations of the SLCO1B1 Polymorphisms With Hepatic Function, Baseline Lipid Levels, and Lipid-lowering Response to Simvastatin in Patients With Hyperlipidemia”  In _Clinical and Applied Thrombosis/Hemostasis_ 24.9_suppl, 2018, pp. 240S–247S  DOI: 10.1177/1076029618805863
  * [56] Anthony D. Gromovsky et al.  “$\Delta$-5 Fatty Acid Desaturase FADS1 Impacts Metabolic Disease by Balancing Proinflammatory and Proresolving Lipid Mediators”  In _Arteriosclerosis, Thrombosis, and Vascular Biology_ 38.1, 2018, pp. 218–231
  * [57] Jorge Simon, Alberto Ouro, Lolia Ala-Ibanibo, Natalia Presa, Teresa Cardoso Delgado and María Luz Martínez-Chantar  “Sphingolipids in Non-Alcoholic Fatty Liver Disease and Hepatocellular Carcinoma: Ceramide Turnover”  In _International Journal of Molecular Sciences_ 21.1 MDPI AG, 2019, pp. 40  DOI: 10.3390/ijms21010040
  * [58] Silvia Bettini et al.  “Liver Fibrosis and Steatosis in Alström Syndrome: A Genetic Model for Metabolic Syndrome”  In _Diagnostics_ 11.5 MDPI AG, 2021, pp. 797  DOI: 10.3390/diagnostics11050797
  * [59] Chaojie Yang et al.  “Impact of Amerind ancestry and FADS genetic variation on omega-3 deficiency and cardiometabolic traits in Hispanic populations”  In _Communications Biology_ 4.1, 2021, pp. 918
  * [60] Xiaoqiang Xiang, Yi Han, Mikko Neuvonen, Marja K. Pasanen, Annikka Kalliokoski, Janne T. Backman, Jouko Laitila, Pertti J. Neuvonen and Mikko Niemi  “Effect of SLCO1B1 polymorphism on the plasma concentrations of bile acids and bile acid synthesis marker in humans”  In _Pharmacogenetics and Genomics_ 19.6 Ovid Technologies (Wolters Kluwer Health), 2009, pp. 447–457  DOI: 10.1097/fpc.0b013e32832bcf7b
  * [61] Puneet Puri et al.  “The presence and severity of nonalcoholic steatohepatitis is associated with specific changes in circulating bile acids”  In _Hepatology_ 67.2 Ovid Technologies (Wolters Kluwer Health), 2017, pp. 534–548  DOI: 10.1002/hep.29359
  * [62] Guangfu Jin et al.  “Genome-wide association study identifies a new locus JMJD1C at 10q21 that may influence serum androgen levels in men”  In _Human Molecular Genetics_ 21.23, 2012, pp. 5222–5228  DOI: 10.1093/hmg/dds361

Appendix for “A statistical framework for GWAS of high dimensional phenotypes
using summary statistics, with application to metabolite GWAS”

##  A1 The dBEMA algorithm

Here we present our algorithm dBEMA, which we described in Section 3.3 and
mirrors Algorithm 2 in [34]. The data matrix input $\hat{\bm{B}}$ is as
defined in (2.2). To be consistent with [34], we parameterize the gamma
distribution as
$\text{Gamma}(\theta,\theta/\phi)\stackrel{{\scriptstyle\text{d}}}{{=}}\phi\times\text{Gamma}(\theta,\theta)$.

Assume $S=S(N)$, $M=M(N)$, $S/N\to\gamma_{S}\in(0,\infty)$, and
$M/N\to\gamma_{M}\in(0,\infty)$, and let $F_{R}$ and $F_{\sigma^{2},N}$ be as
defined in Section 3.3. In what follows, we let $F_{\theta,\phi,N}$ be the
empirical distribution of the non-zero eigenvalues of
$N^{-1}\tilde{\bm{E}}^{\top}\tilde{\bm{E}}$ assuming the $M$ atoms of
$F_{\sigma^{2},N}$ were drawn from a $\text{Gamma}(\theta,\theta/\phi)$. We
note that $F_{\theta,\phi,N}$ will converge to a limiting distribution whose
Stieltjes transform is given in Lemma 1 of [38] if we assume the atoms of
$F_{\sigma^{2},N}$ are drawn from a truncated gamma distribution. This is the
same assumption made in [34].

######  Algorithm A1 (dBEMA)

Input data: Standardized GWAS estimates $\hat{\bm{B}}\in\mathbb{R}^{S\times
M}$, sample size $N$, parameters $\alpha,\beta\in(0,1)$, grid
$0<\theta_{1}<\cdots<\theta_{G}$, and $B\in\mathbb{Z}_{>0}$.
Output: An estimate for $K$, $\hat{K}$.

  1. 1.

Define $\tilde{N}=\min(N,S,M)$ and
$\hat{s}_{1}^{2}\geq\cdots\geq\hat{s}^{2}_{\tilde{N}}$ to be the non-zero
eigenvalues of $N^{-1}\hat{\bm{B}}^{\top}\hat{\bm{B}}$. Estimate $\theta$ and
$\phi$ using the following steps:

    1. (i)

For each
$k\in\\{\lfloor(\alpha/2)\tilde{N}\rfloor,\ldots,\lfloor(1-\alpha/2)\tilde{N}\rfloor\\}$
and $j\in[G]$, determine $q_{kj}$, an estimate of the $(k/\tilde{N})$th
quantile of $F_{\theta_{j},1,N}$.

    2. (ii)

Let $\hat{\phi}_{j}=(\sum_{k}q_{kj}\hat{s}^{2}_{k})/(\sum_{k}q_{kj}^{2})$ for
all $j\in[G]$.

    3. (iii)

Define
$j^{*}=\text{argmin}_{j\in[G]}\sum_{k}(\hat{s}^{2}_{k}-\hat{\phi}_{j}q_{kj})^{2}$
and $\hat{\theta}=\theta_{j^{*}}$, $\hat{\phi}=\hat{\phi}_{j^{*}}$.

  2. 2.

Estimate the distribution of $N^{-1}\tilde{\bm{E}}^{\top}\tilde{\bm{E}}$’s
maximum eigenvalue using the following steps:

    1. (i)

Let $\tilde{\bm{E}}^{(b)}\sim
MN(0,\tilde{\bm{\Sigma}}_{1},\tilde{\bm{\Sigma}}_{2})$, where
$\tilde{\bm{\Sigma}}_{1}\in\mathbb{R}^{S\times S}$ and
$\tilde{\bm{\Sigma}}_{2}\in\mathbb{R}^{M\times M}$ are diagonal matrices whose
elements are sampled from $F_{R}$ and
$\text{Gamma}(\hat{\theta},\hat{\theta}/\hat{\phi})$, respectively.

    2. (ii)

Let $\tilde{\lambda}_{1}^{(b)}$ be the largest eigenvalue of
$N^{-1}\\{\tilde{\bm{E}}^{(b)}\\}^{\top}\tilde{\bm{E}}^{(b)}$.

    3. (iii)

Repeat steps (i)-(ii) and for $b=1,\ldots,B$.

  3. 3.

Let $\hat{K}$ be the number of eigenvalues of
$N^{-1}\hat{\bm{B}}^{\top}\hat{\bm{B}}$ that exceed the $1-\beta$ quantile of
$\\{\tilde{\lambda}_{1}^{(b)}\\}_{b=1}^{B}$.

In all simulations and our data application we let $\alpha=0.4$, $\beta=0.1$,
and $B=500$. Instead of specifying a grid $\theta_{1},\ldots,\theta_{G}$, we
estimated $\theta$ (and subsequently $\phi$) using a line search with lower
and upper end points equal to 0.1 and 5. Lastly, we determined the quantiles
$q_{kj}$ using a modified version of Algorithm 3 in [34], which consisted of
running step 2 but replacing the
$\text{Gamma}(\hat{\theta},\hat{\theta}/\hat{\phi})$ with
$\text{Gamma}(\theta_{j},\theta_{j})$ and $\tilde{\lambda}_{1}^{(b)}$ with
$\tilde{\lambda}_{k}^{(b)}$, the $k$th largest eigenvalue of
$N^{-1}\\{\tilde{\bm{E}}^{(b)}\\}^{\top}\tilde{\bm{E}}^{(b)}$. We defined
$q_{kj}=B^{-1}\sum_{b}\tilde{\lambda}_{k}^{(b)}$.

##  A2 Gibbs sampling using the Chinese restaurant franchise

To perform inference on the parameters in model (4.2), we use a Gibbs sampler
based on the Chinese restaurant franchise (CRF). We start with an analog of
the CRF process. Each pair of parameters $(\mu,\phi^{2})$ corresponds to a
dish, a superpathway $p$ corresponds to a restaurant, a subpathway
$b\in[B_{p}]$ corresponds to a customer. Here, each customer is a group of
metabolites who usually show the same preference for a dish at a table except
for a few outliers. We use $n_{pt}$ to denote the number of customers sitting
at table $t$ in restaurant $p$ and $m_{pk}$ to denote the number of tables
serving dish $k$ in restaurant $p$. Let $n_{p.}=\sum_{t}n_{pt}$,
$m_{.k}=\sum_{p}m_{pk}$ and $m_{..}=\sum_{k}m_{.k}$. We use $t_{pb}$ to denote
the table a subpathway $b$ of superpathway $p$ sits at. $k_{pt}$ denotes the
dish index table $t$ in restaurant $p$ serves.

Our base model described in Section 4.1 is

$H=\pi_{0,0}\delta_{(0,0)}+(1-\pi_{0,0})\operatorname{pr}(\mu)\times\textmd{dIG}(\alpha,\beta)$
(A2.1)

where $\pi_{0,0}$ is the prior probability that a subpathway has loadings
being exactly 0. $\operatorname{pr}(\mu)$ denotes the prior distribution of
mean values, which is a mixture of spike and slab

$\operatorname{pr}(\mu)=\pi_{0}\delta_{0}+\sum_{v=1}^{V}\pi_{v}N(\mu;0,\lambda_{v}^{2})$

where a nonzero $\pi_{0}$ allows for shrinking the means of subpathways to 0.
The slab part is an approximate a Horseshoe prior, which is finite mixture of
multiple normal distributions. The standard deviations $\\{\lambda_{v}\\}$
were chosen as the 0.2, 0.4, 0.6, 0.8 and 0.9 quantiles of a half Cauchy
distribution $C^{+}(0,1)$ and their corresponding weights $\\{\pi_{v}\\}$ are
the standardized and scaled densities of $\\{\lambda_{v}\\}$,

$\pi_{v}^{\prime}=C^{+}(\lambda_{v};0,1),\ \ \
\pi_{v}=\frac{(1-\pi_{0})\pi_{v}^{\prime}}{\sum\pi_{v}^{\prime}}$

$\textmd{dIG}(\alpha,\beta)$ denotes a discritized inverse-gamma prior for
$\phi^{2}$ with shape $\alpha$ and scale $\beta$,

$\textmd{dIG}(\alpha,\beta)=\sum_{w=1}^{W}p_{w}\delta_{\sigma_{w}^{2}}$

where $\\{\sigma_{w}^{2}\\}$ were generated by applying the quantile function
of a inverse-gammma ($\alpha$, $\beta$) distribution to the values
$\\{1/51,2/51,...,50/51\\}$ and their standardized densities were used as
their weights $p_{w}$.

For each factor $k=1,\ldots,K$, after the mean and variance $(\mu,\phi^{2})$
of subpathway $b$ is sampled from the HDP model with base distribution $H$,
concentration parameters $\gamma$ and $\alpha_{0}$, the corresponding entries
of $\bm{L}$ are generated by

$L_{pbm,k}\sim\pi_{out}U[a_{u},b_{u}]+(1-\pi_{out})\mathcal{N}(\mu,\phi^{2}),\
m\in[M_{pb}]$

The estimate $\hat{\bm{L}}$ of $\bm{L}$ is distributed as

$\hat{\bm{L}}_{pbm,k}|\bm{L}_{pbm,k}\sim\mathcal{N}(L_{pbm,k},s_{mk}^{2})$

where $s_{mk}^{2}$ is the variance of the estimator $\hat{\bm{L}}_{pbm,k}$
given $\bm{L}_{pbm,k}$. Here we used a scaling of $\hat{\bm{L}}$ different
from Algorithm 1, which was $\hat{\bm{L}}=\sqrt{M}\bm{V}$ with $M$ being the
number of metabolites and $\bm{V}$ being the first $K$ right singular vectors
as defined in Algorithm 1. With the new scaling, each column of $\hat{\bm{L}}$
has norm $\sqrt{M}$. Scaling $\hat{\bm{L}}$ in this way avoids adjusting the
hyperparameters in the base model according to the strength of each factor and
we can rescale $\hat{\bm{L}}$ and the parameters learned using Gibbs sampler
afterwards. According to Theorem 3.2 (iii), with
$\hat{\bm{L}}=\sqrt{M}\bm{V}$, we have

$s_{mk}^{2}=\frac{M}{\Gamma_{kk}^{2}}(\frac{1}{N}+\frac{1}{S})\hat{\sigma}_{m}^{2}$

We put Beta(1,1) priors on both $\pi_{0,0}$ and $\pi_{0}$, and a
Beta(1/5,999/5) prior on $\pi_{out}$. The latter reflects our prior belief
that we will have on average 1 outlier out of 1,000 metabolites and no more
than 5 outliers with probability 0.95. Following [44], we put Gamma(1,1)
priors on the concentrations parameters $\gamma$ and $\alpha_{0}$. But
different from [44], we discretrized the Gamma prior for $\gamma$ and
$\alpha_{0}$ to avoid time-consuming auxiliary variable sampling.

In our implemention, we first sample the table index $t_{pb}$ for the non-
outlier metabolites of each subpathway, then update the dish indices
$\\{k_{pt}\\}$ for tables and finally update the parameters
$\\{(\mu,\phi^{2})\\}$. The update of other relevant parameters will be
demonstrated in the process. For simplicity of notations, the factor index $k$
is omitted because the Gibbs update will be applied to each column of
$\hat{\bm{L}}$ independently and in parallel. The index $k$ will be used to
indicate parameter $\theta_{k}$ in the content below. We use
$\theta=(\mu,\phi^{2})$ to denote the pair of parameters sampled from the HDP
model (4.2). Thus, for subpathway $b$, the parameter it eventually consumes is
$\theta_{kt_{pb}}=(\mu_{kt_{pb}},\phi^{2}_{kt_{pb}})$.

Updating outliers. For each sub-pathway $b$, we first update the outliers,
because outliers not generated from the HDP model and only the non-outlier
entries will be included in the next round of Gibbs update. The posterior
probability that a metabolite $m$ of subpathway $b$ is an outlier, given
$\hat{\bm{L}}$ and other parameters, can be obtained by integrating $\bm{L}$
out,

$\operatorname{pr}(\textmd{$m$ is an
outlier}|\cdot)\propto\frac{\pi_{out}}{b_{out}-a_{out}}\left[\Phi(\frac{b_{out}-\hat{\bm{L}}_{pbm}}{s_{m}})-\Phi(\frac{a_{out}-\hat{\bm{L}}_{pbm}}{s_{m}})\right]$

where $\Phi$ is the CDF of a standard normal distribution. The parameters of
the uniform distribution, $a_{out}$ and $b_{out}$ should cover the range of
$\hat{\bm{L}}_{k}$. In our implement, $a_{out}$ and $b_{out}$ were set to
centered at zero and cover 1.25 times the range of $\hat{\bm{L}}_{k}$.

The posterior probability that a metabolite is generated from the HDP model
depends on the current Gibbs sample for sub-pathway $b$. If $b$ contains other
non-outlier metabolites, which have been assigned to a table $t$ in the
previous round of Gibbs update, then

$\operatorname{pr}(\textmd{$m$ is not an
outlier}|\cdot)\propto(1-\pi_{out})\operatorname{pr}(\hat{\bm{L}}_{pbm}|\theta_{k_{t_{pb}}},\sigma_{m}^{2})$

where

$\operatorname{pr}(\hat{\bm{L}}_{pbm}|\theta_{k_{t_{pb}}},s_{m})=N(\hat{\bm{L}}_{pbm};\mu_{k_{t_{pb}}},\phi^{2}_{k_{t_{pb}}}+s_{m}^{2}).$

Otherwise, if currently there are no non-outlier metabolites in $b$ except for
$m$, which could happen if $b$ contains only one metabolite or all other
metabolites in $b$ are outliers, we need to integrate $\theta_{k_{t_{pb}}}$
out conditional on the current table and dish assignments,

$\begin{split}\operatorname{pr}(\textmd{$m$ is not an
outlier}|\cdot)\propto&(1-\pi_{out})\left[\sum_{t}\frac{n_{pt}}{n_{p}-1+\alpha_{0}}\operatorname{pr}(\hat{\bm{L}}_{pbm}|\theta_{k_{pt}},s_{m})\right.\\\
+\frac{\alpha_{0}}{n_{p}-1+\alpha_{0}}&\left.\left(\sum_{k}\frac{m_{.k}}{m_{..}+\gamma}\operatorname{pr}(\hat{\bm{L}}_{pbm}|\theta_{k},s_{m})+\frac{\gamma}{m_{..}+\gamma}\operatorname{pr}(\hat{\bm{L}}_{pbm}|H,s_{m})\right)\right]\end{split}$

where

$\operatorname{pr}(\hat{\bm{L}}_{pbm}|H,s_{m})=\pi_{0,0}N(\hat{\bm{L}}_{pbm};0,s_{m}^{2})+(1-\pi_{0,0})\iint
N(\hat{\bm{L}}_{pbm};\mu,\phi^{2}+s_{m}^{2})\operatorname{pr}(\mu)\operatorname{pr}(\phi^{2})d\mu
d\phi^{2}$

The integral on the right hand side of the above equation can be easily
calculated because we discretized $\operatorname{pr}(\mu)$ and
$\operatorname{pr}(\phi^{2})$. Otherwise, the integral is difficult to
evaluate because $N(\hat{\bm{L}}_{pbm};\mu,\phi^{2}+s_{m}^{2})$ is not
conjugate to the inverse-gamma prior due to the existence of $s_{m}^{2}$. This
discretized prior on $\phi^{2}$ helps overcome the non-conjugacy problem and
avoids time-consuming Monte Carlo simulations.

Sampling t. We can treat $t_{pb}$ as the last variable being sampled in
restaurant $p$ due to exchangeability. The likelihood of $t_{pb}=t$ given
other parameters is proportional to the number of subpathways sitting at $t$
if $t$ is not empty, or proportional to $\alpha_{0}$ if $t$ is new. The
likelihood of $\hat{\bm{L}}_{pb}$ given $t_{pb}$ is a product of the
likelihoods of all non-outlier metabolites of subpathway. The likelihood of
$\hat{\bm{L}}_{jg}$ given $t=t_{new}$ can be calculated by integrating out
$k_{pt}$,

$\operatorname{pr}(t_{pb}=t|\cdot)\propto\begin{cases}n_{pt}^{-b}\prod_{m\in[M_{pb}]}\operatorname{pr}(\hat{L}_{pbm}|\theta_{k_{pt}},s_{m})&\textmd{if
$t$ previously used}\\\
\alpha_{0}\prod_{m\in[M_{pb}]}\operatorname{pr}(\hat{L}_{pbm}|t=t_{new},\bm{\theta},\bm{k},s_{m})&\textmd{if
$t=t_{new}$}\end{cases}$ (A2.2)

and

$\operatorname{pr}(\hat{L}_{pbm}|t=t_{new},\bm{\theta},\bm{k},s_{m}^{2})=\sum_{k}\frac{m_{.k}}{m_{..}+\gamma}\prod_{m\in[M_{pb}]}\operatorname{pr}(\hat{\bm{L}}_{pbm}|\theta_{k},s_{m})+\frac{\gamma}{m_{..}+\gamma}\prod_{m\in[M_{pb}]}\operatorname{pr}(\hat{\bm{L}}_{pbm}|H,s_{m})$
(A2.3)

If subpathway $g$ sits at $t_{new}$, we need to assign a dish to $t_{new}$ by
sampling $k_{pt_{new}}$ conditional on data and this will be introduced in the
next step. As a result of updating $t_{pb}$, some table may become unoccupied,
then with probability 0 they will be sampled again in later iterations because
$n_{pt}=0$. As a result, we remove unoccupied tables after updating $t_{pb}$
each time. Since removing tables, as well as updating $k_{pt}$ in the next
step, can potentially cause some dishes unoccupied, we will also delete
unoccupied dishes after updating $\bm{t}$ and $\bm{k}$.

Sampling k. The procedure sampling $k_{pt}$ is similar to sampling $t_{pb}$.
First we assume the exchangeability of dish assignments and treat $k_{pt}$ as
the last variable to be sampled. Since a table can serve multiple subpathways,
the likelihood of $k_{pt}=k$ will contain production across subpathways and
metabolites in those subpathways,

$\operatorname{pr}(k_{pt}=k|\cdot)\propto\begin{cases}m_{.k}\prod_{b:t_{pb}=t}\prod_{m\in[M_{pb}]}\operatorname{pr}(\hat{\bm{L}}_{pbm}|\theta_{k},s_{m})&\text{if
$k$ is previously used}\\\
\gamma\prod_{g:t_{pb}=t}\prod_{m\in[M_{pb}]}\operatorname{pr}(\hat{\bm{L}}_{pbm}|H,s_{m})&\text{if
$k=k_{new}$}\end{cases}$ (A2.4)

If $k_{t}=k_{new}$, we will draw a new pair of parameters
$\theta_{k_{new}}=(\mu_{new},\sigma_{new})$ from $H$ conditional on
$\hat{\bm{L}}_{pbm}$ and $s_{m}$ of all subpathways sitting at table $t$, i.e.
$t_{pb}=t$. We will discretize the inverse-gamma prior on $\phi^{2}$ because
of the non-conjugacy problem arising from $s_{m}^{2}$. See next step for
details.

Updating concentration parameters. Updating $\gamma$ and $\alpha_{0}$ do not
depend on the values of $\bm{\theta}$, it depends on only the number of dishes
and tables, which is fixed after sampling $t$ and $k$. Thus, it was
implemented before updating $\bm{\theta}$.

For updating of concentration parameters, we used the likelihood of $\gamma$
and $\alpha_{0}$ in [44] conditional on $K$, the total number of dishes,
$\\{T_{p}\\}_{1}^{P}$, the number of tables in superpathway $p$ and
$\\{B_{p}\\}_{1}^{P}$, the number of subpathways in superpathway $p$. In our
generative model (4.2) for $\bm{L}$, the concentration parameter $\alpha_{0}$
is shared across all superpathways. Thus, we have

$\operatorname{pr}(T_{1},...T_{P}|\alpha_{0},B_{1},...,B_{P})=\prod_{p=1}^{P}s(B_{p},T_{p})\alpha_{0}^{T_{p}}\frac{\Gamma(\alpha_{0})}{\Gamma(\alpha_{0}+B_{p})}$

and

$\operatorname{pr}(K|\gamma,\sum_{p}T_{p})=s(\sum_{p}T_{p},K)\gamma^{K}\frac{\Gamma(\gamma)}{\Gamma(\gamma+\sum_{p}T_{p})}$

where $s(n,m)$ is unsigned Stirling numbers of the first kind. Instead of
using auxiliary variable sampling or adaptive rejection sampling to
iteratively sample $\gamma$ and $\alpha_{0}$, we discretized the gamma priors
for $\gamma$ and $\alpha_{0}$ and update their values by sampling from their
grids using their posterior probabilities. The grids for both $\gamma$ and
$\alpha_{0}$ were set as $\\{0.1,0.2,...,0.9,1,2,...,10\\}$.

Updating $\bm{\theta}$. Not only sampling new parameters, but also updating
old parameters require sampling the mean values and variances shared by
subpathways conditional on metabolite loading $\hat{\bm{L}}_{pbm}$ and the
variances of error $s_{m}^{2}$. Since parameters $\bm{\theta}$ are shared by
all restaurants(superpathways), we will update $\theta_{k}$ conditional on all
non-outlier metabolites consuming dish $k$ across superpathways and
subpahtways. Let $\\{k\\}$ denote all metabolites consuming dish $k$, we have

$\operatorname{pr}\left(\theta_{k}=(0,0)|\cdot\right)\propto\pi_{0,0}\prod_{m\in\\{k\\}}N(\hat{\bm{L}}_{pbm};0,s_{m}^{2})$
(A2.5)

and

$\operatorname{pr}\left(\theta_{k}\neq(0,0)|\cdot\right)\propto(1-\pi_{0,0})\iint\prod_{m\in\\{k\\}}N(\hat{\bm{L}}_{m};\mu,\phi^{2}+s_{m}^{2})d\mathcal{P}(\mu)d\textmd{IG}(\phi^{2})$
(A2.6)

If $\theta_{k}\neq(0,0)$, we will draw indicator variables

$\operatorname{pr}\left(Z_{v,w}=1|\cdot\right)\propto
p_{v}p_{w}\prod_{m\in\\{k\\}}N(\hat{\bm{L}}_{m};0,\lambda_{v}^{2}+\sigma_{w}^{2}+s_{m}^{2})$
(A2.7)

If $Z_{0,r}=1$, we let $\theta_{k}=(0,\sigma_{w})$ directly. Otherwise, we
draw the value of $\mu_{k}$ given $\lambda_{v}$ and $Z_{v,w}$,

$\sigma_{k}=\sigma_{w},\ \ \mu_{k}\sim
N\left(\frac{\sum_{m\in\\{k\\}}\frac{\hat{\bm{L}}_{m}}{\sigma_{w}^{2}+s_{m}^{2}}+\frac{0}{\lambda_{v}^{2}}}{\sum_{m\in\\{k\\}}\frac{1}{\sigma_{w}^{2}+s_{m}^{2}}+\frac{1}{\lambda_{v}^{2}}},\frac{1}{\sum_{m\in\\{k\\}}\frac{1}{\sigma_{w}^{2}+s_{m}^{2}}+\frac{1}{\lambda_{v}^{2}}}\right)$
(A2.8)

Updating $\pi_{out},\pi_{0,0},\pi_{0}$. The updates of $\pi_{out}$,
$\pi_{0,0}$ and $\pi_{0}$ can be obtained by

$\pi_{out}\sim\textmd{Beta}(\frac{1}{5}+\\#outliers,\frac{999}{5}+M-\\#outliers)$
$\pi_{0,0}\sim\textmd{Beta}(1+\\#spikes,1+K-\\#spikes)$
$\pi_{0}\sim\textmd{Beta}(1+\\#centers,1+K-\\#spikes-\\#centers)$

where $\\#outliers$, $\\#spikes$ and $\\#centers$ denotes the number of
outliers, the number of parameters $\theta=(0,0)$ and the number of parameters
with mean 0 but non-zero variances. $M$ and $K$ denote the total number of
metabolites and number of parameters in the current round of Gibbs update.

##  A3 Estimating direct effects in mtGWAS: Setting grid points

Here we describe our procedure to define the grid points in
$\mathcal{X}^{(\pi)}$ and $\mathcal{X}^{(\varphi^{2})}$ defined in Section
4.2. Briefly, we first estimate $\\{\pi_{pbm}\\}$ and $\varphi^{2}$ via
maximum likelihood as

$\displaystyle\\{\\{\hat{\pi}_{pbm}\\},\hat{\varphi}^{2}\\}=\mathop{\text{argmax}}_{\\{\pi_{pbm}\\},\varphi^{2}}\operatorname{pr}(\hat{\bm{\Delta}}\mid\\{\pi_{pbm}\\},\varphi^{2}),$

assuming

$\displaystyle\hat{\bm{\Delta}}_{s,pbm}\mid\bm{\Delta}_{s,pbm}\sim
N(\bm{D}_{ss}^{1/2}\bm{\Delta}_{s,pbm},\hat{\sigma}_{pbm}^{2})$
$\displaystyle\bm{\Delta}_{s,pbm}\mid\pi_{pbm},\varphi^{2}\sim(1-\pi_{pbm})\delta_{0}+\pi_{pbm}N(0,\varphi^{2}\hat{\sigma}_{pbm}^{2}/\bm{D}_{ss}).$

We then define $\mathcal{X}^{(\pi)}$ and $\mathcal{X}^{(\varphi^{2})}$ to each
be sets with 100 elements, where
$\mathcal{X}^{(\varphi^{2})}=\\{\hat{\varphi}^{2}/2,\ldots,2\hat{\varphi}^{2}\\}$
contains equally spaced elements. We set $\mathcal{X}^{(\pi)}$ so that its
minimum element is $10/S$, its maximum element is 1.5 times the 0.99 quantile
of $\\{\hat{\pi}_{pbm}\\}$, and is equally spaced on a log scale.

##  A4 Simulating L for experiments

To simulate a $\bm{L}$ of $K$ factors, a super-sub pathway structure of
metabolites need to be specify first. In our implement the number of
subpathways of a superpathway and sizes of subpathways were randomly sampled
from Poisson distributions. In summary, we generated a metabolite archives
containing 3 superpathways and 257 metabolites.

Given the pathway structure of metabolites, each column of $\bm{L}$ was
generated independently. The generating procedure of $\bm{L}$ exactly followed
the Chinese restaurant franchise construction of the HDP model in Equation
(A2.1). The global shrinkage parameter $\tau$ of the Horseshoe prior was set
as 0.2, outlier rate $\pi_{out}$ was 0.01, spike rate $\pi_{0,0}$ was 0.1,
$\pi_{0}$ was randomly sampled from a Beta(2,1) distribution and the scale and
rate for the inverse gamma distribution where $\phi^{2}$ will be drawn were
set as 2 and 1. To make sure the outliers were recoverable, the outlier
entries were generated beyond 2 standard deviations away from the mean of the
subpathway, but they cannot exceed (-8,8). The simulation procedure allowed us
to track which subpathways were spikes, which metabolites were outliers and
which subpathways shared the same parameters. Finally, all entries of $\bm{L}$
were multiplied by 0.4 to make the strength of signals comparable to the noise
level.

We lastly show that the Gibbs sampler described in Section A2 can identify
spikes and outliers using the asymptotic variance of $\hat{\bm{L}}$ derived in
Theorem 3.2. To do so, we simulated a phenotype matrix $\bm{Y}$ with 257
metabolites, 5000 SNPs and sample size being 500, making the ratio
$\frac{S}{N}=10$. $\bm{Y}$ contained $K=1$ factor to simplify the simulation
by avoiding identifiability issues that arise with more than one factor. The
$\bm{L}$ used to simulate $\bm{Y}$ is shown in Figure A1 Left, which contains
3 spike subpathways (marked in red) and 3 outliers (marked by inverted
triangles). The Gibbs sampler described in Section A2 was applied to the
estimate $\hat{\bm{L}}$ twice, once assuming
$\operatorname{Var}(\hat{\bm{L}}_{mk}\mid\bm{L})=(1/N+1/S)\sigma_{m}^{2}$ and
another assuming
$\operatorname{Var}(\hat{\bm{L}}_{mk}\mid\bm{L})=(1/S)\sigma_{m}^{2}$. The
second is the variance one would assume if they treat $\bm{R}$ in (2.2) as the
identity [28]. The results are shown in Figure A1 Middle and Right. A
subpathway was concluded as a spike if the posterior probability that it
consumed (0,0) as its parameter exceeded 0.95. Outliers were defined as the
metabolites whose posterior probability of being an outlier was greater than
0.95. It can be observed that the Gibbs sampler was able to identify spikes
and outliers correctly with the right variance, but it failed identifying any
spike subpathways and overestimates the number of outliers.

Figure A1: Left: True $\bm{L}$ simulated using the model in Equation 4.2.
Spike pathways are marked in red. Outliers are marked by inverted triangles.
Middle: Estimates for spikes and outliers with the correct asymptotic variance
of $\hat{\bm{L}}$. Estimated spike pathways are marked in red. Estimated
Outliers are marked by inverted triangles. Right: Estimates for spikes and
outliers with the incorrect variance for $\hat{\bm{L}}$ derived by ignoring
the correlation between SNPs. Estimated spike pathways are marked in red.
Estimated Outliers are marked by inverted triangles.