# From differential abundance to mtGWAS: accurate and scalable methodology for
metabolomics data with non-ignorable missing observations and latent factors

Shangshu Zhao1, Kedir Turi2, Tina Hartert2, Carole Ober3,
Klaus Bønnelykke4, Bo Chawes4, Hans Bisgaard4, Chris McKennan1,∗
1Department of Statistics, University of Pittsburgh
2Department of Medicine, Vanderbilt University Medical Center
3Department of Human Genetics, University of Chicago
4COPSAC, Copenhagen Prospective Studies on Asthma in Childhood,
Herlev and Gentofte Hospital, University of Copenhagen
***To whom correspondence should be addressed<EMAIL_ADDRESS>
###### Abstract

Metabolomics is the high-throughput study of small molecule metabolites.
Besides offering novel biological insights, these data contain unique
statistical challenges, the most glaring of which is the many non-ignorable
missing metabolite observations. To address this issue, nearly all analysis
pipelines first impute missing observations, and subsequently perform analyses
with methods designed for complete data. While clearly erroneous, these
pipelines provide key practical advantages not present in existing
statistically rigorous methods, including using both observed and missing data
to increase power, fast computation to support phenome- and genome-wide
analyses, and streamlined estimates for factor models. To bridge this gap
between statistical fidelity and practical utility, we developed MS-NIMBLE, a
statistically rigorous and powerful suite of methods that offers all the
practical benefits of imputation pipelines to perform phenome-wide
differential abundance analyses, metabolite genome-wide association studies
(mtGWAS), and factor analysis with non-ignorable missing data. Critically, we
tailor MS-NIMBLE to perform differential abundance and mtGWAS in the presence
of latent factors, which reduces biases and improves power. In addition to
proving its statistical and computational efficiency, we demonstrate its
superior performance using three real metabolomic datasets.

Keywords: Metabolomics; Metabolomic GWAS; Latent factors; Factor analysis;
Confounding; MNAR

##  1 Introduction

Metabolomics is the high-throughput study of small molecule metabolites, and
can help understand human variation and the etiology of disease [1].
Metabolite abundances are typically measured via mass spectrometry, which,
while sensitive, produces a large amount of non-ignorable missing data in
which low abundance metabolites are less likely to be observed [2]. This
precludes the use of the many complete data methods able to perform the three
core metabolomic analyses: differential abundance, metabolome genome-wide
association studies (mtGWAS), and factor analysis [3, 2]. Factor analysis,
while important in its own right, is required in differential abundance
analyses and, as we show in Section 4.4, mtGWAS, as it helps recover latent
factors that plague metabolomic data and confound relationships of interest
[2].

Consequently, nearly all existing analysis pipelines first impute missing
data, which acts as a crude solution to issues of method incompatibility [4]
and offers the following important practical advantages: (i) ensuing analyses
use both observed and missing data to improve power, (ii) downstream
computation is fast enough to perform metabolite phenome- and genome-wide
studies, and (iii) factor models can be estimated. Despite its expedience, it
is well known that imputing non-ignorable missing data can beget biased
estimators and spurious inference [5]. However, to our knowledge, [2] is the
only work to provide a rigorous alternative to imputation while also
considering latent confounding factors. Although a step in the right
direction, their work does not offer the aforementioned advantages of
imputation, as it discards missing data and does not provide methodology to
perform an mtGWAS. And while it does provide a method to perform factor
analysis, its theoretical properties are completely unknown. Therefore, it is
questionable whether the statistical rigor offered by [2] is sufficient to
offset the expediency of imputation.

To bridge the gap between statistical fidelity and practical utility, we
developed MS-NIMBLE (Methods for Non-Ignorable Missing Metabolomic
Observations), a suite of statistically rigorous methods to perform
differential abundance, mtGWAS, and factor analysis in metabolomic data that
offers all of the practical advantages of imputation. Like [2], we estimate
each metabolite’s missingness mechanism once per dataset and store it to
facilitate efficient downstream computation. However, unlike [2], subsequent
estimators use both observed and missing data by leveraging the approximate
conditional normality of metabolite levels. Our method for mtGWAS is able to
partition low rank and idiosyncratic genetic variation, and we prove the
statistical and computational efficiency of our factor analysis-related and
other estimators. We lastly use simulated and three real metabolomic datasets
to show that MS-NIMBLE significantly outperforms the method proposed in [2]
and existing imputation pipelines. An R package and code to reproduce our
simulations are available from https://github.com/chrismckennan/MSNIMBLE.

##  2 Notation, problem setup, and statistical models

Let $[m]=\\{1,\ldots,m\\}$ for $m>0$ and $y_{gi}$ be the possibly missing log-
abundance of metabolite $g\in[p]$ in sample $i\in[n]$. For observed covariates
$\bm{x}_{i}\in\mathbb{R}^{d}$ and latent factors
$\bm{c}_{i}\in\mathbb{R}^{K}$, assume

$\displaystyle
y_{gi}=\bm{\beta}_{g}^{\top}\bm{x}_{i}+\bm{\ell}_{g}^{\top}\bm{c}_{i}+e_{gi},\quad(e_{g1},\ldots,e_{gn})^{\top}\sim
N(0,\sigma_{g}^{2}I_{n}),\quad g\in[p];i\in[n]$ (2.1)

for some unknown and non-random $\bm{\beta}_{g}\in\mathbb{R}^{d}$ and
$\bm{\ell}_{g}\in\mathbb{R}^{K}$. We will assume the number of latent factors
$K$ is known, although we estimate $K$ with parallel analysis [6] in practice.
In differential abundance, $\bm{\beta}_{g}$ is of interest and $\bm{c}_{i}$
confounds the relationship between $\bm{x}_{i}$ and $y_{gi}$. In factor
analysis and mtGWAS, $\bm{\beta}_{g}$ is a nuisance parameter and $\bm{c}_{i}$
and $\bm{\ell}_{g}$ are of interest. Other than assuming the design matrix
with rows $(\bm{x}_{i}^{\top},\bm{c}_{i}^{\top})$, $i\in[n]$, is full rank, we
assume nothing about the relationship between $\bm{x}_{i}$ and $\bm{c}_{i}$,
which facilitates the analysis of data with arbitrarily complex latent
confounding. While our theoretical results require assumptions on the moments
of $\bm{c}_{i}$, our methodology is agnostic to these assumptions, and
therefore postpone their discussion to Section 5. The normality of $e_{gi}$,
which we leverage to design efficient estimators, is a common assumption in
mass spectrometry data [4, 7]. However, we do not require $y_{gi}$ be normal,
as the elements of $\bm{c}_{i}$ are often highly skewed (Figure S1).

It is well known metabolite levels depend on genotype [3]. However, since
genotype does not appear in (2.1), it is possible that its effect is mediated
by $\bm{c}_{i}$ or appears in the idiosyncratic error terms $e_{gi}$, which
belies the canonical factor analysis assumption that $\bm{c}_{i}$ is
independent of $e_{gi}$ [8]. The genetic effects in $e_{gi}$ also imply the
normality of $e_{gi}$ may only be an approximation, and that $e_{gi}$,
$e_{hi}$ may be dependent for $g\neq h$. Our theoretical work in Section 5
accommodate all of these observations.

To describe the missing data model, let $r_{gi}=I(y_{gi}$ is observed$)$. We
follow [2] and assume that for some known cumulative distribution function
$\Psi$ and unknown, metabolite-specific scale and location parameters
$\alpha_{g}\geq 0$ and $\delta_{g}\in\mathbb{R}$,

$\displaystyle\operatorname{Pr}(r_{gi}=1\mid
y_{gi},\bm{x}_{i},\bm{c}_{i})=\operatorname{Pr}(r_{gi}=1\mid
y_{gi})=\Psi\\{\alpha_{g}(y_{gi}-\delta_{g})\\},\quad g\in[p];i\in[n],$ (2.2)

where $\\{r_{gi}\\}_{g\in[p];i\in[n]}$ are independent conditional on
$\\{y_{gi}\\}_{g\in[p];i\in[n]}$. This, along with the assumptions that
$\alpha_{g}\geq 0$ and the distribution of $r_{gi}$ only depends on $y_{gi}$,
is justified because nearly all missing data are due to an artifact of the
mass spectrometer, where analytes with low abundances are less likely to be
observed [2]. [2] contains additional justifications of (2.2).

We assume $\Psi$ in (2.2) is known, which is ostensibly allowed to be any
cumulative distribution function (CDF). While typical choices for $\Psi$
include the CDFs of the logistic and normal distributions [7], our theoretical
work in Section 5 requires the left hand tail of $\Psi$ go to zero no faster
than a polynomial rate. Our default choice for $\Psi$ is therefore the CDF of
the t-distribution with four degrees of freedom, which we show gives excellent
results in real data.

##  3 When do the missing data matter?

Ignoring or incorrectly modeling non-ignorable missing data can bias
estimators [5]. Despite this, differential abundance simulations routinely
suggest that errant imputation techniques have a trivial effect on type I
error [7]. This begs the question when, or if, we have to account for the non-
ignorable missing data in metabolomic analyses. We study this in Proposition
3.1, which analyzes estimates from errantly imputed data.

######  Proposition 3.1.

Let $x_{i}\in\mathbb{R}$. Assume (2.1) satisfies
$y_{gi}=\mu_{g}+x_{i}\beta_{g}+\bm{c}_{i}^{\top}\bm{\ell}_{g}+e_{gi}$,
$(e_{g1},\ldots,e_{gn})^{\top}\sim N(0,\sigma_{g}^{2}I_{n})$, and the
regularity conditions in Section S7 hold. Suppose (2.2) holds and we impute
missing $y_{gi}$’s as $a*\min(\\{y_{gi}\\}_{\\{i:r_{gi}=1\\}})$ for any
constant $a\in\mathbb{R}$. Then for $\hat{\beta}_{g},\hat{s}_{g}$ the
resulting ordinary least squares estimate and standard error for $\beta_{g}$
when $\bm{c}_{i}$ is known, $(\hat{\beta}_{g}-\beta_{g})/\hat{s}_{g}\to
N(0,1)$ as $n\to\infty$ if (i) the null hypothesis $H_{0,g}:\beta_{g}=0$ holds
and (ii) $\bm{\ell}_{g}=0$ or $x_{i}$ is independent of $\bm{c}_{i}$.

######  Remark 3.1.

Minimum imputation from Proposition 3.1 is one of the most common ways to
handle missing metabolomic data [4]. Note $\bm{c}_{i}$ is observed in
Proposition 3.1.

Proposition 3.1 shows errant imputation can beget valid type I error rates
provided (ii) holds, i.e. $\bm{c}_{i}$ does not confound the relationship
between $x_{i}$ and $y_{gi}$. This result explains the abovementioned
befuddling observations that incorrectly modelling simulated non-ignorable
missing metabolomic data has a trivial effect on type I error rates, since
their simulations did not consider confounders.

The proof of the asymptotic normality in Proposition 3.1 relies on $x_{i}$
being independent of $y_{gi}$, which is only true if (i) and (ii) hold. This
suggests properly handling missing $y_{gi}$’s is critical when estimating
intervals for non-zero effects $\beta_{g}$, and when controlling type I error
in the presence of confounding factors $\bm{c}_{i}$, even when $\bm{c}_{i}$ is
observed. We show this using simulated and real data.

##  4 Estimation and inference with MS-NIMBLE

We must overcome several challenging features of (2.1), (2.2), and metabolomic
experiments in general. First, (2.1) is not congruent with existing maximum
likelihood estimators designed for normally distributed data [7], since
$\bm{c}_{i}$’s distribution may be highly non-normal (Figure S1). Second,
leveraging the approximate normality of the errors $e_{gi}$ to improve
estimates requires integrating over missing $y_{gi}$, which can be
prohibitively slow for theoretically valid choices of $\Psi$ discussed in
Section 5. Lastly, our estimators must scale to facilitate phenome- and
genome-wide analyses. Figure 1 gives an overview of the steps in our method.
For simplicity of presentation, we assume in Sections 4 and 5 that all
metabolites have missing data, but provide extensions in supplemental Section
S5 to allow fully observed metabolites. Section 4.1 gives a brief description
of the estimators for $\alpha_{g},\delta_{g}$, as they mirror those from [2].
Sections 4.2-4.4 contain detailed descriptions of our novel methodological
components.

Figure 1: Method overview and how estimators are used to solve different
problems in metabolomics.

###  4.1 Estimating the missingness mechanisms

We follow [2] and estimate $\alpha_{g},\delta_{g}$ from (2.2) using a Bayesian
generalized method of moments estimator. Briefly, for some observed
$\bm{u}_{gi}\in\mathbb{R}^{r}$, we consider the observable sample moment
$\bm{m}_{g}(\tilde{\alpha},\tilde{\delta})=n^{-1/2}\sum_{i=1}^{n}\bm{u}_{gi}[1-r_{gi}/\Psi\\{\tilde{\alpha}(y_{gi}-\tilde{\delta})\\}]$,
which is mean 0 and asymptotically normal when
$(\tilde{\alpha},\tilde{\delta})=(\alpha_{g},\delta_{g})$ and $\bm{u}_{gi}$ is
independent of $r_{gi}$ conditional on $y_{gi}$. Treating $\bm{m}_{g}$ as our
“data”, we estimate $\alpha_{g},\delta_{g}$ as
$(\hat{\alpha}_{g},\hat{\delta}_{g})=\operatorname{\mathbb{E}}\\{(\alpha_{g},\delta_{g})\mid\bm{m}_{g}(\alpha_{g},\delta_{g})\\}$,
where we approximate the posterior
$\operatorname{Pr}\\{\alpha_{g},\delta_{g}\mid\bm{m}_{g}(\alpha_{g},\delta_{g})\\}\propto\operatorname{Pr}\\{\bm{m}_{g}(\alpha_{g},\delta_{g})\mid\alpha_{g},\allowbreak\delta_{g}\\}\operatorname{Pr}(\alpha_{g},\delta_{g})$
assuming $\bm{m}_{g}(\alpha_{g},\delta_{g})$ is normally distributed. Since
$y_{gi}$ must be dependent on $\bm{u}_{gi}$, we let
$\bm{u}_{gi}\in\mathbb{R}^{r}$ be $r$ of the first few principal components of
the data matrix of fully observed metabolites. Sections 3 and 4 of [2] contain
additional details.

Critically, the estimators $\hat{\alpha}_{g},\hat{\delta}_{g}$ only depend on
the dataset $\\{r_{gi}y_{gi}\\}_{g\in[p];i\in[n]}$, and are invariant to the
covariate of interest $\bm{x}_{i}$ and genotype. We therefore only compute
$\hat{\alpha}_{g},\hat{\delta}_{g}$ once per dataset and store the results,
which helps make downstream analyses computationally tractable.

###  4.2 Estimating latent factors

We describe estimates for latent factors $\bm{c}_{i}$ in differential
abundance problems, and show how these can be used to derive estimates in
factor analysis and mtGWAS applications as well in Sections 4.3 and 4.4. Let
$\bm{X}=(\bm{x}_{1}\cdots\bm{x}_{n})^{\top}$ and
$P_{\bm{X}}^{\perp}\in\mathbb{R}^{n\times n}$ be the orthogonal projection
matrix that projects vectors onto the kernel of $\bm{X}^{\top}$. We can
express $\bm{C}=(\bm{c}_{1}\cdots\bm{c}_{n})^{\top}\in\mathbb{R}^{n\times K}$
as $\bm{C}=P_{\bm{X}}^{\perp}\bm{C}+\bm{X}\bm{\Omega}$, where
$\bm{\Omega}=(\bm{X}^{\top}\bm{X})^{-1}\bm{X}^{\top}\bm{C}$. Model (2.1) can
then be re-written as

$\displaystyle
y_{gi}=\bm{b}_{g}^{\top}\bm{x}_{i}+\bm{\ell}_{g}^{\top}[P_{\bm{X}}^{\perp}\bm{C}]_{i*}+e_{gi},\quad\bm{b}_{g}=\bm{\beta}_{g}+\bm{\Omega}\bm{\ell}_{g},\quad(e_{g1},\ldots,e_{gn})\sim
N(0,\sigma_{g}^{2}I_{n}),$ (4.1)

where $[P_{\bm{X}}^{\perp}\bm{C}]_{i*}\in\mathbb{R}^{K}$ is the $i$th row of
$P_{\bm{X}}^{\perp}\bm{C}$. We utilize the paradigm from [8] and sequentially
estimate $P_{\bm{X}}^{\perp}\bm{C}$ and $\bm{\Omega}$, where the latter
estimate adjusts for confounding. It seems natural to use the normality of
$e_{gi}$ in (4.1) to obtain optimal maximum likelihood estimates for
$P_{\bm{X}}^{\perp}\bm{C}$ and $\bm{\Omega}$. However, this would beget
computationally expensive iterative algorithms that require numerically
integrating over all missing $y_{gi}$’s at each iteration. Instead, we use
inverse probability weighting (IPW) to derive computationally efficient
estimators. Remarkably, we prove in Section 5 that the loss of statistical
efficiency that accompanies IPW has an asymptotically negligible effect on
downstream inference.

If $\bm{Y}=[y_{gi}]\in\mathbb{R}^{p\times n}$ were observed, a natural
estimate for $P_{\bm{X}}^{\perp}\bm{C}$ is the first $K$ right singular
vectors of $\bm{Y}P_{\bm{X}}^{\perp}$ [8], which is equivalent to minimizing
$\sum_{g,i}\allowbreak\\{\allowbreak
y_{gi}\allowbreak-\allowbreak(\bm{b}_{g}^{\top}\bm{x}_{i}\allowbreak+\allowbreak\bm{\ell}_{g}^{\top}\bm{C}^{\perp}_{i*})\\}^{2}$
over $\bm{C}^{\perp}$, as well as $\bm{b}_{g}$ and $\bm{\ell}_{g}$, such that
$\bm{X}^{\top}\bm{C}^{\perp}=0$. This motivates estimating
$P_{\bm{X}}^{\perp}\bm{C}$ when $y_{gi}$’s may be missing by solving the
following IPW-inspired optimization problem:

$\displaystyle\begin{aligned}
&\\{P_{\bm{X}}^{\perp}\hat{\bm{C}},\\{\tilde{\bm{b}}_{g},\tilde{\bm{\ell}}_{g}\\}_{g\in[p]}\\}\in\operatorname*{argmin}_{\begin{subarray}{c}\bm{C}^{\perp}\in\mathbb{R}^{n\times
K},\,\bm{b}_{g}\in\mathbb{R}^{d},\,\bm{\ell}_{g}\in\mathbb{R}^{K}\\\
\text{such that
$\bm{X}^{\top}\bm{C}^{\perp}=0$}\end{subarray}}\sum_{g=1}^{p}\sum_{i=1}^{n}\hat{w}_{gi}f_{gi}(\bm{C}^{\perp},\bm{b}_{g},\bm{\ell}_{g})\\\
&f_{gi}(\bm{C}^{\perp},\bm{b}_{g},\bm{\ell}_{g})=\\{y_{gi}-(\bm{b}_{g}^{\top}\bm{x}_{i}+\bm{\ell}_{g}^{\top}\bm{C}^{\perp}_{i*})\\}^{2},\quad\hat{w}_{gi}=r_{gi}/\Psi\\{\hat{\alpha}_{g}(y_{gi}-\hat{\delta}_{g})\\}.\end{aligned}$
(4.2)

Here, $\hat{\alpha}_{g},\hat{\delta}_{g}$ are given in Section 4.1 and the
objective on the first line is observable because $\hat{w}_{gi}=0$ if $y_{gi}$
is missing. If $\hat{\alpha}_{g}=\alpha_{g}$, $\hat{\delta}_{g}=\delta_{g}$,
and we replace $\hat{w}_{gi}$ with its expectation
$\operatorname{\mathbb{E}}(\hat{w}_{gi}\mid y_{gi})=1$, the above discussion
implies (4.2) is equivalent to singular value decomposition. Unlike maximum
likelihood estimators that use the normality of $e_{gi}$ in (4.1), iterative
updates in (4.2) have a closed form and beget fast computation. Since
$P_{\bm{X}}^{\perp}\bm{C}$ is not identifiable in (4.1),
$P_{\bm{X}}^{\perp}\hat{\bm{C}}$ is not unique. While we address this in
factor analysis applications by requiring $P_{\bm{X}}^{\perp}\hat{\bm{C}}$
have orthogonal columns, such an identification criterion is unnecessary in
differential abundance and mtGWAS.

To recover $\bm{\Omega}$, we see (4.2) provides estimates
$\tilde{\bm{b}}_{g},\tilde{\bm{\ell}}_{g}$ for $\bm{b}_{g},\bm{\ell}_{g}$ in
(4.1). Since $P_{\bm{X}}^{\perp}\bm{C}$ is orthogonal to $\bm{X}$, we should
be able to separate variation due to $P_{\bm{X}}^{\perp}\bm{C}$ and $\bm{X}$,
which suggests $\tilde{\bm{b}}_{g},\tilde{\bm{\ell}}_{g}$ are reasonably
accurate. If $\bm{\beta}_{g}=0$ for all $g\in[p]$, then the expression for
$\bm{b}_{g}$ in (4.1) indicates we can estimate $\bm{\Omega}$ by regressing
$(\tilde{\bm{b}}_{1}\cdots\tilde{\bm{b}}_{p})$ onto
$(\tilde{\bm{\ell}}_{1}\cdots\tilde{\bm{\ell}}_{p})$. While not all
$\bm{\beta}_{g}$’s will be 0, covariates of interest encoded by $\bm{X}$
typically correlate with only a few metabolites [1]. We use this to justify
estimating $\bm{\Omega}$ with the aforementioned regression:

$\displaystyle\hat{\bm{\Omega}}=\textstyle\operatorname*{argmin}_{\bm{\Omega}\in\mathbb{R}^{d\times
K}}\textstyle\sum_{g=1}^{p}\lVert\tilde{\bm{b}}_{g}-\bm{\Omega}\tilde{\bm{\ell}}_{g}\rVert_{2}^{2}=(\textstyle\sum_{g=1}^{p}\tilde{\bm{b}}_{g}\tilde{\bm{\ell}}_{g}^{\top})(\sum_{g=1}^{p}\tilde{\bm{\ell}}_{g}\tilde{\bm{\ell}}_{g}^{\top})^{-1}.$
(4.3)

In addition to adjusting for latent confounds, we show $\hat{\bm{\Omega}}$ can
be used to test if latent factors depend on $\bm{X}$ in Sections 5 and 6. We
show in supplemental Section S4 that (4.3) can be further refined by
iteratively removing “outlying” covariate-dependent metabolites from the
regression.

We estimate $\bm{C}$ as
$\hat{\bm{C}}=P_{\bm{X}}^{\perp}\hat{\bm{C}}+\bm{X}\hat{\bm{\Omega}}$ in
differential abundance problems. Since $\bm{X}$ is a nuisance covariate in
factor analysis and mtGWAS, we let $\hat{\bm{C}}$ be the solution to (4.2) in
those applications.

###  4.3 Estimation and inference on coefficients of interest

Here we consider
$\bm{\theta}_{g}=(\bm{\beta}_{g}^{\top},\bm{\ell}_{g}^{\top})^{\top}$, where
$\bm{\beta}_{g}$ is the inferential target in differential abundance and
$\bm{\ell}_{g}$ is important in factor analysis and mtGWAS. Our goal is to
develop statistically efficient estimators that can be computed quickly.
Throughout Section 4.3, we let
$\hat{\bm{z}}_{i}=(\bm{x}_{i}^{\top},\allowbreak\hat{\bm{c}}_{i}^{\top})^{\top}$
for $\hat{\bm{c}}_{i}\in\mathbb{R}^{K}$ the $i$th row of $\hat{\bm{C}}$
defined in Section 4.2 (our estimate for $\bm{c}_{i}$ in (2.1)).

Having estimated
$\\{\alpha_{g},\delta_{g},\bm{z}_{i}=(\bm{x}_{i}^{\top},\bm{c}_{i}^{\top})^{\top}\\}$
as $\\{\hat{\alpha}_{g},\hat{\delta}_{g},\hat{\bm{z}}_{i}\\}$ in Sections 4.1
and 4.2, we consider estimating $\bm{\theta}_{g}$ and $\sigma_{g}$ via the
log-likelihood $h_{g}(\bm{\theta},\sigma)$ of the observed data
$\\{r_{gi}y_{gi}\\}_{i\in[n]}$ implied by (2.1) and (2.2) using the plug-in
estimators $\\{\hat{\alpha}_{g},\hat{\delta}_{g},\hat{\bm{z}}_{i}\\}$:

$\displaystyle\begin{aligned}
h_{g}(\bm{\theta},\sigma)=&\textstyle\sum_{i=1}^{n}-r_{gi}\\{\log(\sigma)+(y_{gi}-\bm{\theta}^{\top}\hat{\bm{z}}_{i})^{2}/(2\sigma^{2})\\}\\\
&+\textstyle\sum_{i=1}^{n}(1-r_{gi})\log[1-\smallint\Psi\\{\hat{\alpha}_{g}(\bm{\theta}^{\top}\hat{\bm{z}}_{i}+\sigma
e-\hat{\delta}_{g})\\}\phi(e)\text{d}e],\end{aligned}$ (4.4)

where $\phi(e)$ is the standard normal density. While Section S8.6 of the
Supplement shows that directly maximizing $h_{g}$ will accurately estimate
$\bm{\theta}_{g}$, this fails to consider the computational cost of
numerically integrating the second line of (4.4). To address this, we design
an appropriately initialized algorithm that only requires a small number of
iterations, and therefore numerical integrations, to accurately approximate
the maximizer of (4.4). Briefly, for $\hat{w}_{gi}$ given in (4.2), let

$\displaystyle\begin{aligned}
\hat{\bm{\theta}}_{g}^{\operatorname{(IPW)}}&=(\textstyle\sum_{i=1}^{n}\hat{w}_{gi}\hat{\bm{z}}_{i}\hat{\bm{z}}_{i}^{\top})^{-1}(\textstyle\sum_{i=1}^{n}\hat{w}_{gi}\hat{\bm{z}}_{i}y_{gi})\\\
\hat{\sigma}_{g}^{\operatorname{(IPW)}}&=[(\textstyle\sum_{i=1}^{n}\hat{w}_{gi})^{-1}\textstyle\sum_{i=1}^{n}\hat{w}_{gi}\\{y_{gi}-\hat{\bm{z}}_{i}^{\top}\hat{\bm{\theta}}_{g}^{\operatorname{(IPW)}}\\}^{2}]^{1/2}\end{aligned}$
(4.5)

be the inverse probability weighted (IPW) estimators of $\bm{\theta}_{g}$ and
$\sigma_{g}$. Since $\hat{w}_{gi}=0$ if $y_{gi}$ is missing, the estimators in
(4.5) only use observed data, and are therefore sub-optimal. However, they are
easy to compute and, as we show in supplemental Section S8.6, consistent,
which make them appropriate starting points. We then iteratively update our
estimates for $\bm{\theta}_{g}$ and $\sigma_{g}$ with Fisher scoring using the
information matrix
$\mathcal{I}_{g}(\bm{\theta},\sigma)=\operatorname{\mathbb{E}}_{\\{\bm{\theta},\sigma\\}}[\nabla^{2}h_{g}(\bm{\theta},\sigma)\mid\\{\hat{\bm{z}}_{i}\\}_{i\in[n]}]$,
where the expectation ignores the uncertainty in $\hat{\bm{z}}_{i}$,
$\hat{\alpha}_{g}$, and $\hat{\delta}_{g}$. While running this algorithm to
completion is potentially computationally expensive, we prove in Section 5.3
that we only require one Fisher scoring step to achieve asymptotically optimal
estimates. In practice, our software default is $\leq 10$ iterations. Letting
$\hat{\bm{\theta}}_{g}=(\hat{\bm{\beta}}_{g}^{\top},\hat{\bm{\ell}}_{g}^{\top})^{\top}$
and $\hat{\sigma}_{g}$ be the resulting estimates, we perform inference on
$\bm{\beta}_{g}$ assuming $\hat{\bm{\beta}}_{g}\approx
N(\bm{\beta}_{g},\hat{\operatorname{\mathbb{V}}}(\hat{\bm{\beta}}_{g}))$ for
$\hat{\operatorname{\mathbb{V}}}(\hat{\bm{\beta}}_{g})$ the first $d\times d$
block of
$\\{-\mathcal{I}_{g}(\hat{\bm{\theta}}_{g},\hat{\sigma}_{g})\\}^{-1}$.

Two features of this procedure cast doubt on its fidelity. The first is the
assumption in (2.1) that $e_{gi}$ is normally distributed, as the existence of
genetic and possibly other non-normal variation in $e_{gi}$ suggest the
likelihood in (4.4) is incorrect. While this is not a concern in fully
observed data, estimates from missing data may be sensitive to distributional
assumptions [5]. The second is $\hat{\bm{\beta}}_{g}$ depends on the estimated
latent factors $\hat{\bm{c}}_{1},\ldots,\hat{\bm{c}}_{n}$ whose theoretical
properties are unknown. We address these concerns in Section 5.3, where we
prove inference with $\hat{\bm{\beta}}_{g}$ is asymptotically equivalent to
knowing both the non-normal genetic effects and latent factors. While the
uncertainty in $\hat{\alpha}_{g},\hat{\delta}_{g}$ ostensibly poses a third
issue, the strong theoretical and simulation results in [2] proving their
accuracy suggest this is trivial.

###  4.4 Metabolite genome-wide association study

We lastly consider performing an mtGWAS. We set $\bm{x}_{i}$ in (2.1) to be 0
for simplicity, but show in supplemental Section S9 how to extend our method
to allow $\bm{x}_{i}\neq 0$. Let $G_{si}\in\\{0,1,2\\}$ be the genotype at
single nucleotide polymorphism (SNP) $s$ in sample $i$. Given (2.1), the
effect of $G_{si}$ on $y_{gi}$ can either appear in the idiosyncratic error
$e_{gi}$, or be mediated by $\bm{c}_{i}$. We therefore assume
$e_{gi}=\gamma^{(e)}_{gs}G_{si}+\Delta_{gi}^{(e)}$ and
$\bm{c}_{i}=\bm{\gamma}_{s}^{(c)}G_{si}+\bm{\Delta}_{i}^{(c)}$, where
$\gamma^{(e)}_{gs}\in\mathbb{R}$, $\bm{\gamma}_{s}^{(c)}\in\mathbb{R}^{K}$
quantify the effect of $G_{si}$ on $e_{gi}$ and $\bm{c}_{i}$, respectively,
and $\Delta_{gi}^{(e)}\in\mathbb{R}$, $\bm{\Delta}_{i}^{(c)}\in\mathbb{R}^{K}$
are mean 0 errors. This implies

$\displaystyle
y_{gi}=\\{\bm{\ell}_{g}^{\top}\bm{\gamma}_{s}^{(c)}+\gamma^{(e)}_{gs}\\}G_{si}+\\{\bm{\ell}_{g}^{\top}\bm{\Delta}_{i}^{(c)}+\Delta_{gi}^{(e)}\\},$
(4.6)

where $\gamma^{(e)}_{gs}$ and $\bm{\ell}_{g}^{\top}\bm{\gamma}_{s}^{(c)}$ are
interpretable as the idiosyncratic and low rank genetic effects. We develop
methodology below to perform inference on $\gamma^{(e)}_{gs}$,
$\bm{\ell}_{g}^{\top}\bm{\gamma}_{s}^{(c)}$, and the total effect
$\bm{\ell}_{g}^{\top}\bm{\gamma}_{s}^{(c)}+\gamma^{(e)}_{gs}$.

Consider testing $H_{0,gs}^{(e)}:\gamma_{gs}^{(e)}=0$. Classic Wald tests
would require optimizing (4.4) for all #metabolites $\times$ #SNPs pairs $g$
and $s$. While this is reasonable for #SNPs $\lesssim 10^{2}$ (i.e. on the
order of a phenome-wide association study), it is infeasible in genome-wide
studies, where #SNPs $\gtrsim 10^{6}$. To circumvent this, we propose a novel
and tractable score test. Briefly, consider the log-likelihood
$h_{gs}(\gamma,\bm{\ell},\sigma)$ for $\\{r_{gi}y_{gi}\\}_{i\in[n]}$ under
(2.1) and (2.2) assuming $e_{gi}\sim N(\gamma G_{si},\sigma^{2})$:

$\displaystyle h_{gs}(\gamma,\bm{\ell},\sigma)=$
$\displaystyle\textstyle\sum_{i=1}^{n}-r_{gi}[\log(\sigma)+\\{y_{gi}-(\bm{\ell}^{\top}\hat{\bm{c}}_{i}+\gamma
G_{si})\\}^{2}/(2\sigma^{2})]$
$\displaystyle+\textstyle\sum_{i=1}^{n}(1-r_{gi})\log[1-\smallint\Psi\\{\hat{\alpha}_{g}(\bm{\ell}^{\top}\hat{\bm{c}}_{i}+\gamma
G_{si}+\sigma e-\hat{\delta}_{g})\\}\phi(e)\text{d}e].$

If $H_{0,gs}^{(e)}:\gamma_{gs}^{(e)}=0$ is true,
$h_{gs}\\{\gamma_{gs}^{(e)},\bm{\ell},\sigma\\}=h_{g}(\bm{\ell},\sigma)$ for
$h_{g}$ as defined in (4.4). Then for $\hat{\bm{\ell}}_{g},\hat{\sigma}_{g}$
the approximate maximizers of $h_{g}$ described in Section 4.3, we define the
score statistic $\eta_{gs}^{(e)}$ to be

$\displaystyle\eta_{gs}^{(e)}=\\{\textstyle\tfrac{\partial}{\partial\gamma}h_{gs}(\gamma,\hat{\bm{\ell}}_{g},\hat{\sigma}_{g})\mid_{\gamma=0}\\}^{2}[\\{-\mathcal{I}_{gs}(0,\hat{\bm{\ell}}_{g},\hat{\sigma}_{g})\\}^{-1}]_{11},$
(4.7)

where $\mathcal{I}_{gs}(\gamma,\bm{\ell},\sigma)$ is the Fisher information
matrix assuming $h_{gs}(\gamma,\bm{\ell},\sigma)$ is the log-likelihood for
$\\{r_{gi}y_{gi}\\}_{i\in[n]}$. A p-value for $H_{0,gs}^{(e)}$ is computed by
comparing $\eta_{gs}^{(e)}$ to the upper quantiles of a $\chi^{2}_{1}$.

Several features of (4.7) make our test computationally and statistically
efficient. First, since $\hat{\bm{\ell}}_{g},\hat{\sigma}_{g}$ are the
approximate maximizers of $h_{g}$ in (4.4), they do not depend on genotype,
and consequently only need to be computed once per metabolite $g$. Therefore,
as we show in supplemental Section S9, (4.7) is a simple function of genotype
and metabolite-specific terms that can be pre-computed. Second, (4.7) uses all
available data and does not errantly impute missing data, which is the
prevailing practice in mtGWAS studies. Lastly, and most importantly, inference
with (4.7) is done conditional on the estimated latent factors
$\hat{\bm{c}}_{i}$, which de-noises the data to substantially improve power by
reducing residual variances. For example, we show that the variance reduction
in our data example is equivalent to increasing the sample size by 67%.

We next consider $\bm{\ell}_{g}^{\top}\bm{\gamma}_{s}^{(c)}$ from (4.6), which
is interpretable as the effect of SNP $s$ on metabolite $g$ that is mediated
through the latent factors $\bm{c}_{i}$. Let $\hat{\bm{\ell}}_{g}$ as defined
above, and let $\hat{\operatorname{\mathbb{V}}}(\hat{\bm{\ell}}_{g})$ be its
its estimated variance obtained using the inferential procedure outlined in
Section 4.3. Since $\bm{\gamma}_{s}^{(c)}$ satisfies
$\operatorname{\mathbb{E}}(\bm{c}_{i}\mid
G_{si})=\bm{\gamma}_{s}^{(c)}G_{si}$, we define $\hat{\bm{\gamma}}_{s}^{(c)}$
and $\hat{\operatorname{\mathbb{V}}}\\{\hat{\bm{\gamma}}_{s}^{(c)}\\}$ to be
the ordinary least squares estimate and its corresponding estimated variance
from the regression of $[\hat{\bm{c}}_{1}\cdots\hat{\bm{c}}_{n}]^{\top}$ onto
$(G_{s1}\cdots G_{sn})^{\top}$, which can be efficiently computed at the
genome-wide scale. If $\hat{\bm{c}}_{i}=\bm{c}_{i}$ and there were no genetic
effects on $e_{gi}$, standard arguments can be used to show
$\hat{\bm{\ell}}_{g}$ is asymptotically independent of
$\hat{\bm{\gamma}}_{s}^{(c)}$. We therefore test
$H_{0,gs}^{(c)}:\bm{\ell}_{g}^{\top}\bm{\gamma}_{s}^{(c)}=0$ by comparing the
following to the upper quantiles of a $\chi^{2}_{1}$:

$\displaystyle\eta_{gs}^{(c)}=\\{\hat{\bm{\ell}}_{g}^{\top}\hat{\bm{\gamma}}_{s}^{(c)}\\}^{2}/[\hat{\bm{\ell}}_{g}^{\top}\widehat{\operatorname{\mathbb{V}}}\\{\hat{\bm{\gamma}}_{s}^{(c)}\\}\hat{\bm{\ell}}_{g}+\\{\hat{\bm{\gamma}}_{s}^{(c)}\\}^{\top}\widehat{\operatorname{\mathbb{V}}}(\hat{\bm{\ell}}_{g})\hat{\bm{\gamma}}_{s}^{(c)}].$
(4.8)

We lastly test whether SNP $s$ has any effect on metabolite $g$’s abundance.
Given (4.6), the classic approach would test the null that
$\bm{\ell}_{g}^{\top}\bm{\gamma}_{s}^{(c)}+\gamma_{gs}^{(e)}=0$. However, as
discussed above, this is not practical because it would require estimating
$\gamma_{gs}^{(e)}$. Instead, since $\bm{c}_{i}$ and $e_{gi}$ are typically
assumed to be independent in metabolomic data [2], we assume their
corresponding genetic effects reflect unrelated variation. This suggests a
metabolite’s abundance is genetically regulated if
$\bm{\ell}_{g}^{\top}\bm{\gamma}_{s}^{(c)}$ or $\gamma_{gs}^{(e)}$ is 0. We
therefore propose testing
$H_{0,gs}^{(c,e)}:\bm{\ell}_{g}^{\top}\bm{\gamma}_{s}^{(c)}=\gamma_{gs}^{(e)}=0$
using $\eta_{gs}^{(c,e)}=\eta_{gs}^{(c)}+\eta_{gs}^{(e)}$, which we show in
Section 5.4 is approximately $\chi^{2}_{2}$ under $H_{0,gs}^{(c,e)}$.

##  5 Theoretical guarantees

Here we justify estimators and inference from Section 4. Since [2] detailed
the theoretical properties of $\hat{\alpha}_{g},\hat{\delta}_{g}$ defined in
Section 4.1, we focus on the properties and impact of the latent factor
estimates $\hat{\bm{c}}_{i}$ from Section 4.2, as their theoretical properties
are unknown but critical to the fidelity of estimators proposed in Sections
4.2-4.4. Given the accuracy of $\hat{\alpha}_{g},\hat{\delta}_{g}$ [2] and the
negligible impact their uncertainty has in real and simulated data (see
Sections 6 and S2), we assume
$\hat{\alpha}_{g}=\alpha_{g},\hat{\delta}_{g}=\delta_{g}$ to make proofs
tractable, which is common in the non-random missing data literature [9].

Section 5.1 details our assumptions, and Sections 5.2-5.4 contain our
theoretical results. In addition to providing the theoretical foundation for
estimators in Section 4, these results help us specify a software default
choice for $\Psi$ defined in (2.2). All proofs are in the supplement.

###  5.1 Assumptions

Let $\bm{X}=[\bm{x}_{1}\cdots\bm{x}_{n}]^{\top}\in\mathbb{R}^{n\times d}$,
$\bm{C}=[\bm{c}_{1}\cdots\bm{c}_{n}]^{\top}\in\mathbb{R}^{n\times K}$, and
$\bm{1}_{n}=(1,\ldots,1)^{\top}\in\mathbb{R}^{n}$. For
$\bm{M}\in\mathbb{R}^{n\times m}$, let
$P_{\bm{M}}^{\perp}\in\mathbb{R}^{n\times n}$ be the orthogonal projection
onto the kernel of $\bm{M}^{\top}$. We first place assumptions on $y_{gi}$.

######  Assumption 5.1.

For $g\in[p]$, $i\in[n]$, and $s\in[S]$, let
$y_{gi}=\bm{\beta}_{g}^{\top}\bm{x}_{i}+\bm{\ell}_{g}^{\top}\bm{c}_{i}+e_{gi}$,
$G_{si}\in\\{0,1,2\\}$, and $\mathcal{G}=\\{G_{si}\\}_{s\in S;i\in[n]}$. Then
the following hold for constants $a_{1}>0$ and $\epsilon\in(0,1/2\wedge
a_{1})$.

  1. (a)

$\bm{X}=[\tilde{\bm{X}},\bm{1}_{n}]$ is non-random,
$n^{-1}\tilde{\bm{X}}^{\top}P_{\bm{1}_{n}}^{\perp}\tilde{\bm{X}}\succeq\epsilon
I_{d-1}$,
$\lVert\tilde{\bm{X}}\rVert_{\infty},\lVert\bm{\beta}_{g}\rVert_{2}\leq
a_{1}$, $\mathcal{G}$’s elements are independent, $\\{G_{si}\\}_{i\in[n]}$ are
identically distributed for each $s\in[S]$, and $\epsilon n\leq p\leq a_{1}n$.

  2. (b)

The eigenvalues $\lambda_{1},\ldots,\lambda_{K}>0$ of
$p^{-1}\sum_{g=1}^{p}\bm{\ell}_{g}\bm{\ell}_{g}^{\top}$ satisfy
$n^{-1/2+\epsilon}\allowbreak\lesssim\allowbreak\lambda_{K}\allowbreak\leq\allowbreak\cdots\leq\allowbreak\lambda_{1}\allowbreak\lesssim
1$, $\lambda_{1}/\lambda_{K}\leq a_{1}$, and
$\left\lVert\bm{\ell}_{g}\right\rVert_{2}\leq a_{1}\lambda_{1}^{1/2}$.
Further,
$\bm{c}_{i}=\bm{f}(\bm{x}_{i})+\sum_{s=1}^{S}\bm{\gamma}_{s}^{(c)}G_{si}+\bm{\Delta}_{i}^{(c)}\in\mathbb{R}^{K}$,
where:

    1. (i)

$\bm{f}:\mathbb{R}^{d}\to\mathbb{R}^{K}$ is a continuous function and
$\\{\bm{\gamma}_{s}^{(c)}\\}_{s\in[S]}$ are non-random and satisfy
$\sum_{s=1}^{S}\lVert\bm{\gamma}_{s}^{(c)}\rVert_{2}\leq a_{1}$,
$\sum_{s=1}^{S}1\\{\bm{\gamma}_{s}^{(c)}\neq 0\\}\leq a_{1}p^{1/2}$, and
$\max_{s\in[S]}\lVert\bm{\gamma}_{s}^{(c)}\rVert_{2}=o(n^{-1/4})$.

    2. (ii)

$\\{\bm{\Delta}_{i}^{(c)}\\}_{i\in[n]}$ are independent, identically
distributed, independent of $\mathcal{G}$,
$\operatorname{\mathbb{V}}\\{\bm{\Delta}_{i}^{(c)}\\}\succeq\epsilon I_{K}$,
and
$\operatorname{\mathbb{E}}\\{\lvert\bm{\Delta}_{i_{k}}^{(c)}\rvert^{m}\\}\leq
b_{m}$ for $k\in[K]$, all $m>0$, and constants $b_{m}>0$.

  3. (c)

For non-random parameters $\\{\gamma_{gs}^{(e)}\\}_{g\in[p];s\in[S]}$,
$e_{gi}=\sum_{s=1}^{S}\gamma_{gs}^{(e)}G_{si}+\Delta_{gi}^{(e)}$ such that:

    1. (i)

$\sum_{s=1}^{S}1\\{\gamma_{gs}^{(e)}\neq 0\\}\leq a_{1}$,
$\max_{g\in[p];s\in[S]}\lvert\gamma_{gs}^{(e)}\rvert=o(n^{-1/4})$,
$\Delta_{gi}^{(e)}\sim N(0,\sigma_{g}^{2})$, $\sigma_{g}^{2}\leq a_{1}$, and
$\\{\Delta_{gi}^{(e)}\\}_{g\in[p];i\in[n]}$ are independent and are
independent of $\\{\mathcal{G},\bm{C}\\}$.

    2. (ii)

Each connected component of the metabolite graph created by placing an edge
between metabolites $g,h\in[p]$ if $\gamma_{gs}^{(e)}\gamma_{hs}^{(e)}\neq 0$
has $\leq a_{1}$ metabolite vertices.

We require $\bm{X}$ contain an intercept in (a). The assumptions on genotype
$G_{si}$ in (a) are akin to assuming each linkage disequilibrium block
contains at most one causal SNP. The eigenvalues in (b) quantify the average
magnitude of $\bm{\ell}_{1},\allowbreak\ldots,\allowbreak\bm{\ell}_{p}$, where
we let eigenvalues be moderate ($\asymp n^{-1/2+\epsilon}$) or large ($\asymp
1$). While some datasets may have eigenvalues even smaller than
$n^{-1/2+\epsilon}$, they likely make a trivial contribution to metabolite
variation and are therefore not considered here.

Since metabolites may be genetically regulated, we allow latent factors
$\bm{c}_{i}$ and errors $e_{gi}$ to be dependent on genotype. This implies
$\bm{c}_{i}$ and $e_{gi}$ may be dependent, which violates the assumptions of
most factor analysis methods [8]. To our knowledge, our theoretical work is
the first to consider genetic dependence between latent factors and errors.

We assume genetic effects $\bm{\gamma}^{(c)}_{s}$ and $\gamma^{(e)}_{gs}$
decay with sample size, which is a common assumption in GWAS [10]. However, we
will have asymptotically perfect power if the genetic effect is $\gtrsim
n^{-1/2+\eta}$ in magnitude for any $\eta>0$ and the number of tested SNPs is
polynomial in $n$. Assumption (c)(c)(ii) assumes metabolites can be
partitioned into pathways where, conditional on latent factors, metabolites in
different pathways are independent, which is a common assumption [1]. We next
place assumptions on the missing data.

######  Assumption 5.2.

Model (2.2) and the following hold for some constants $a_{2}>1$, $m>0$:

  1. (a)

$\\{r_{gi}\\}_{g\in[p];i\in[n]}$ are independent conditional on
$\\{y_{gi}\\}_{g\in[p];i\in[n]}$ and
$\alpha_{g}\in(0,a_{2}),\lvert\delta_{g}\rvert\leq a_{2}$.

  2. (b)

$\Psi$ is a six times continuously differentiable CDF that satisfies (i)
$\Psi(-x)=1-\Psi(x)$, (ii) $\lvert x\rvert^{m}\Psi(x)\geq a_{2}^{-1}$ for all
$x<-a_{2}$, and (iii) $\lvert
x\rvert^{m}\lvert\frac{\text{d}^{(j)}}{\text{d}x^{(j)}}\Psi(x)\rvert\leq
a_{2}$ for all $j\in[6]$ and $\lvert x\rvert>a_{2}$.

######  Remark 5.2.

Assumption (b) is satisfied when $\Psi$ is the CDF of a t-distribution.

Section 2 discusses the conditional independence assumption in (a). Assumption
(b)(ii) requires the left hand tail of $\Psi$ to go to 0 at a polynomial rate,
which ensures the inverse probability weighted estimator in (4.2) is well-
behaved. Remark 5.2 inspires our software-default choice for $\Psi$ to be the
CDF of a t-distribution with four degrees of freedom, which also reduces the
impact of outlying observations on our estimates for $\bm{\beta}_{g}$ (see
supplemental Remark S8.14). Note (b)(ii) excludes the usual assumption that
$\Psi$ is the CDF of a logistic or normal random variable [7], as their left
hand tails go to 0 at exponential and super-exponential rates.

###  5.2 Accuracy of and inference with latent factor estimates

We first consider the accuracy of $P_{\bm{X}}^{\perp}\hat{\bm{C}}$ defined in
(4.2), which is critical to the estimate for $\hat{\bm{C}}$ in differential
abundance and is exactly $\hat{\bm{C}}$ in factor analysis applications.

######  Theorem 5.1.

Suppose Assumptions 5.1 and 5.2 hold, and let
$\hat{\mathcal{P}},\mathcal{P}\in\mathbb{R}^{n\times n}$ be the orthogonal
projections that project vectors onto
$\operatorname{Im}(P_{\bm{X}}^{\perp}\hat{\bm{C}})$ and
$\operatorname{Im}(P_{\bm{X}}^{\perp}\bm{C})$. Then there exists a constant
$\eta>0$ such that if we require
$\lVert\hat{\mathcal{P}}-\mathcal{P}\rVert_{F}\leq\eta$, then
$\lVert\hat{\mathcal{P}}-\mathcal{P}\rVert_{F}^{2}=o_{P}(n^{-1/2})$.

######  Remark 5.3.

The objective in (4.2), which is expressed as a function of the matrix
parameter $\bm{C}^{\perp}$, only depends on $\bm{C}^{\perp}$ through
$\operatorname{Im}(\bm{C}^{\perp})$, and is therefore actually a function of
orthogonal projection matrices. The requirement that
$\lVert\hat{\mathcal{P}}-\mathcal{P}\rVert_{F}\leq\eta$ implies the desired
minimizer of (4.2) may only be a local minima, and we implicitly assume
$\lVert\hat{\mathcal{P}}-\mathcal{P}\rVert_{F}\leq\eta$ in all future
theoretical statements. We show in supplemental Section S5 that, under minor
conditions, we can guarantee
$\lVert\hat{\mathcal{P}}-\mathcal{P}\rVert_{F}\leq\eta$ by initializing (4.2)
using metabolites with fully observed data.

Theorem 5.1 is, to our knowledge, the first result proving the fidelity of
factor analysis in data with non-random missing observations. Remarkably, this
result mirrors the best known factor analysis results when data are observed
[8], and accounts for possible genetic-related dependencies between
$\bm{c}_{i}$ and $e_{gi}$, which are not allowed to exist in most factor
analysis-related theoretical results [11, 8].

We next consider our estimate for $\bm{\Omega}$ from (4.3), which helps ensure
our estimates for $\bm{\beta}_{g}$ are not biased by latent factors
$\bm{c}_{i}$. While its theoretical properties derived in supplemental Section
S8.5 are critical for Sections 5.3 and 5.4, we show in Theorem 5.2 below that
it can also be used to formally test whether $\bm{c}_{i}$ confounds the
relationship between $\bm{x}_{i}$ and $y_{gi}$.

######  Theorem 5.2.

Fix a $j\in[d-1]$. In addition to Assumptions 5.1 and 5.2, suppose (i)
$p^{-1}\sum_{g=1}^{p}\allowbreak 1\\{\bm{\beta}_{g_{j}}\neq
0\\}=o(\lambda_{1}^{1/2}n^{-1/2})$ and (ii)
$\operatorname{\mathbb{E}}(\bm{c}_{i})=\bm{A}^{\top}\bm{x}_{i}$ for some non-
random $\bm{A}\in\mathbb{R}^{d\times K}$. Then if the null hypothesis
$H_{0,j}:\bm{A}_{j*}=0$ is true,
$\hat{\bm{\Omega}}_{j*}^{\top}\hat{\bm{\Omega}}_{j*}/\tilde{x}_{j}^{2}\operatorname{\stackrel{{\scriptstyle\text{d}}}{{\to}}}\chi^{2}_{K}$,
where $\bm{A}_{j*},\hat{\bm{\Omega}}_{j*}\in\mathbb{R}^{K}$ are the $j$th rows
of $\bm{A},\hat{\bm{\Omega}}$ and $\tilde{x}_{j}^{2}$ is the $j$th diagonal of
$(\bm{X}^{\top}\bm{X})^{-1}$.

######  Remark 5.4.

The sparsity assumption in (i) is weaker than the usual assumption
$p^{-1}\sum_{g=1}^{p}\allowbreak 1\\{\bm{\beta}_{g_{j}}\neq
0\\}=o(\lambda_{1}n^{-1/2})$ made by methods that require fully observed data
[8], since $\lambda_{1}<\lambda_{1}^{1/2}$ if $\lambda_{1}<1$. Note (i) is
only required for the $j$th coefficient.

###  5.3 The statistical and computational efficiency of differential
abundance estimates

We next consider our estimate for $\bm{\beta}_{g}$ from Section 4.3. While we
want to ensure its statistical fidelity, we are also interested studying its
computational efficiency, since maximizing the likelihood in (4.4) requires
expensive numerical integrations. We first state a proposition.

######  Proposition 5.2.

Suppose Assumptions 5.1 and 5.2 hold, let
$h^{\operatorname{(known)}}_{g}(\bm{\beta}_{g},\bm{\ell}_{g},\sigma_{g})$ be
the log-likelihood for $\\{r_{gi}y_{gi}\\}_{i\in[n]}$ when $\bm{C}$ and
$\\{\operatorname{\mathbb{E}}(e_{gi}\mid\mathcal{G})\\}_{i\in[n]}$ are known,
and let $\hat{\bm{\beta}}_{g}^{\operatorname{(known)}}$ be $\bm{\beta}_{g}$’s
corresponding consistent maximum likelihood estimate. Then
$\\{\bm{V}_{g}^{\operatorname{(known)}}\\}^{-1/2}\allowbreak\\{\hat{\bm{\beta}}_{g}^{\operatorname{(known)}}\allowbreak-\allowbreak\bm{\beta}_{g}\\}\operatorname{\stackrel{{\scriptstyle\text{d}}}{{\to}}}N(0,I_{d})$
for $\bm{V}_{g}^{\operatorname{(known)}}$ the first $d\times d$ block of
$[-\operatorname{\mathbb{E}}\\{\nabla^{2}h^{\operatorname{(known)}}_{g}(\bm{\beta}_{g},\bm{\ell}_{g},\sigma_{g}^{2})\mid\bm{C},\mathcal{G}\\}]^{-1}$.

Unsurprisingly, estimates are asymptotically normal when we observe the full
covariate matrix $[\bm{X},\bm{C}]$ and know the genetic effects
$\\{\operatorname{\mathbb{E}}(e_{gi}\mid\mathcal{G})\\}_{i\in[n]}$. The latter
is important, since the missing data likelihood is incorrect when the non-
normal genetic effects are unknown, which risks biasing estimates. Remarkably,
Theorem 5.3 shows that our estimator for $\bm{\beta}_{g}$, which replaces
$\bm{C}$ with its estimate from Section 4.2 and ignores genetic effects, is
asymptotically equivalent to $\hat{\bm{\beta}}_{g}^{\operatorname{(known)}}$.

######  Theorem 5.3.

Let $d_{1}\leq d-1$. In addition to Assumptions 5.1 and 5.2, assume (i) in the
statement of Theorem 5.2 holds for all $j\in[d_{1}]$. Suppose we initialize
the optimization to maximize (4.4) at the IPW estimates defined in (4.5), and
let $\hat{\bm{\beta}}_{g}$ be the estimate for $\bm{\beta}_{g}$ after updating
the IPW estimates with one Fisher scoring step. Then for
$\hat{\bm{\beta}}_{g}^{\operatorname{(known)}}$ and
$\bm{V}_{g}^{\operatorname{(known)}}$ defined in Proposition 5.2,

$\displaystyle
n^{1/2}\lVert\hat{\bm{\beta}}_{g(1:d_{1})}-\hat{\bm{\beta}}_{g(1:d_{1})}^{\operatorname{(known)}}\rVert_{2}=o_{P}(1),\quad
n\lVert\hat{\operatorname{\mathbb{V}}}(\hat{\bm{\beta}}_{g})_{(1:d1)}-\bm{V}_{g(1:d_{1})}^{\operatorname{(known)}}\rVert_{2}=o_{P}(1),$
(5.1)

where
$\hat{\bm{\beta}}_{g(1:d_{1})},\hat{\bm{\beta}}_{g(1:d_{1})}^{\operatorname{(known)}}\in\mathbb{R}^{d_{1}}$
are the first $d_{1}$ elements of
$\hat{\bm{\beta}}_{g},\hat{\bm{\beta}}_{g}^{\operatorname{(known)}}$. The
matrices $\bm{V}_{g(1:d_{1})}^{\operatorname{(known)}}$ and the observable
$\hat{\operatorname{\mathbb{V}}}(\hat{\bm{\beta}}_{g})_{(1:d1)}$ are the first
$d_{1}\times d_{1}$ blocks of $\bm{V}_{g}^{\operatorname{(known)}}$ and the
minus inverse Fisher information for the likelihood $h_{g}$ in (4.4) evaluated
at the first Fisher scoring step, respectively.

Result (5.1) indicates both the estimate $\hat{\bm{\beta}}_{g(1:d_{1})}$ and
corresponding inference using
$\hat{\operatorname{\mathbb{V}}}(\hat{\bm{\beta}}_{g})_{(1:d1)}$ is
asymptotically equivalent to that when both $\bm{C}$ and genetic effects are
known. Together with Proposition 5.2, this justifies using standard Wald
intervals and tests to perform inference.

Two features of Theorem 5.3 imply our estimates are computationally efficient.
First, we need only apply a single iteration of Fisher scoring per metabolite.
While we allow more than one iteration in practice, convergence is fast (see
supplemental Section S3). Second, Theorem 5.3 indicates differential abundance
inference incurs no cost when using the computationally efficient, but
statistically sub-optimal, IPW-based estimate for $\bm{C}$ in Section 4.2.
This is critical, since the likelihood-based estimate is prohibitively slow to
compute due to repeated numerical integration.

###  5.4 Fidelity of latent factor-corrected mtGWAS

Here we justify our mtGWAS method from Section 4.4. Recall $\eta_{gs}^{(e)}$,
$\eta_{gs}^{(c)}$, and $\eta_{gs}^{(c,e)}$ are the test statistics that test
whether the genotype at SNP $s$ affects metabolite $g$’s idiosyncratic
variation $e_{gi}$, low-dimensional variation
$\bm{\ell}_{g}^{\top}\bm{c}_{i}$, and total variation
$\bm{\ell}_{g}^{\top}\bm{c}_{i}+e_{gi}$. As we did in Section 4.4, we assume
$\bm{X}=0$ for simplicity, but show in Section S9 that the extension to
general $\bm{X}$ is simple.

######  Theorem 5.4.

Fix a $g\in[p]$, suppose $\bm{X}=0$ and Assumptions 5.1 and 5.2 hold, and let
$\gamma_{gs}^{(e)},\bm{\gamma}_{s}^{(c)}$ be as defined in Assumption 5.1.
Then
$\eta_{gs}^{(e)}\operatorname{\stackrel{{\scriptstyle\text{d}}}{{\to}}}\chi^{2}_{1}$
if $H_{0,gs}^{(e)}:\gamma_{gs}^{(e)}=0$ is true. If
$n^{1/2}\left\lVert\bm{\ell}_{g}\right\rVert_{2}\to\infty$, then
$\eta_{gs}^{(c)}\operatorname{\stackrel{{\scriptstyle\text{d}}}{{\to}}}\chi^{2}_{1}$
if $H_{0,gs}^{(c)}:\bm{\ell}_{g}^{\top}\bm{\gamma}_{s}^{(c)}=0$ is true and
$\eta_{gs}^{(c,e)}\operatorname{\stackrel{{\scriptstyle\text{d}}}{{\to}}}\chi^{2}_{2}$
if
$H_{0,gs}^{(c,e)}:\gamma_{gs}^{(e)}=\bm{\ell}_{g}^{\top}\bm{\gamma}_{s}^{(c)}=0$
is true.

######  Remark 5.5.

The non-trivial effect of latent factors suggests
$n^{1/2}\left\lVert\bm{\ell}_{g}\right\rVert_{2}$ is large for most $g$.

##  6 Real data analysis

We used three metabolomic datasets to evaluate our method MS-NIMBLE. Table 1
describes the data, which were collected from the plasma of children that were
part of the Copenhagen Prospective Study on Asthma in Childhood (COPSAC) [12]
or Infant Susceptibility to Pulmonary Infections and Asthma Following RSV
Exposure Study (INSPIRE) [13] cohorts. We partitioned metabolites into
“observed” metabolites ($<5\%$ missing data) and metabolites with missing data
($\geq 5\%$ but $\leq 50\%$ missing data), and discarded metabolites with
$>50\%$ missing data. We were primarily interested in metabolites with missing
data. Supplemental Section S2 provides simulations further demonstrating MS-
NIMBLE’s superior performance.

Table 1: An overview of the real data analyzed in Section 6, where $m_{g}$ is
the fraction of metabolite $g$’s observations that are missing. The sixth and
seventh columns indicate whether a differential abundance (DA) analysis was
performed using respiratory-related traits and sex. The last column indicates
if the dataset was used to perform the mtGWAS.

###  6.1 Real data differential abundance analyses

Since the COPSAC and INSPIRE studies were designed to investigate respiratory
illness through childhood, we first used MS-NIMBLE to identify respiratory-
related metabolites. Specifically, we considered the phenotypes specific
airway resistance (sRAW), a measure of airway patency in the COPSAC cohort,
and infant wheeze, defined as whether the infant wheezed during the first year
of life in the INSPIRE cohort. Since there was no evidence of sRAW-related
metabolites in infancy, we did not consider the 0.5 year COPSAC dataset in
this analysis.

We compared MS-NIMBLE’s estimators for and inference on $\bm{\beta}_{g}$ from
Section 4.3 to two competing approaches. The first, MetabMiss [2], uses the
estimates for the missingness mechanism parameters from Section 4.1 and takes
a similar approach as that in Section 4.2 to recover latent factors. However,
its estimates for $\bm{\beta}_{g}$ discard missing data, and are therefore
expected to be substantially less powerful than MS-NIMBLE. The second imputes
missing data using one of minimum imputation, singular value decomposition
(SVD), K-nearest neighbors (KNN), or random forest (RF), the four most
commonly used imputation techniques [4], and subsequently estimates
$\bm{\beta}_{g}$ using the latent factor-correction method CATE [11]. While
many methods can adjust for latent factors in imputed data, we found CATE gave
the best simulation results in supplemental Section S2. To facilitate inter-
method comparisons, the number of latent factors was set to be the same for
each method, and, as done previously [2, 1], was estimated via parallel
analysis applied to metabolites with no missing data. Supplemental Section S3
contains additional details, including method-specific software settings.

Figure 2(a) gives the number of respiratory-associated metabolites with
missing data identified by each method at a q-value threshold of 0.2. As
expected, MS-NIMBLE identifies over three times as many metabolites as
MetabMiss, where the three piperine metabolites identified by MetabMiss, whose
relationship with sRAW has previously been explored [2], were also identified
by MS-NIMBLE (Figure 2(b)). Figure 2(c) provides a biological explanation for
the remaining metabolites in Figure 2(b) uniquely identified by MS-NIMBLE,
which helps argue the veracity of MS-NIMBLE’s identifications. The small
p-values in Figure 2(a), which test the null hypothesis from Theorem 5.2,
suggest latent factors confound the relationship between the two respiratory
traits and metabolite levels. As a consequence, Section 3 and supplemental
Section S2 suggest imputation methods are inflating type I error rates,
thereby casting doubt on their results.

Figure 2: Respiratory-related differential abundance results. (a): The number
of metabolites with missing data identified at a q-value threshold of 0.2. MS-
NIMBLE’s p-value is the p-value for the null hypothesis from Theorem 5.2 that
the trait is not related to the latent factors. (b): Named metabolites with
missing data that were identified by MS-NIMBLE and MetabMiss (second column)
and MS-NIMBLE but not MetabMiss (third column). MS-NIMBLE identified two
unnamed wheeze-associated metabolites. (c): Biological plausibility of
metabolites from (b) uniquely identified by MS-NIMBLE. Superscripts are 1:
[14]; 2: [15]; 3: [1]. (d)-(e): Effect estimates and 95% confidence intervals
for selected metabolites. Numbers in parentheses are the fractions of missing
metabolite data.

Having argued hypothesis testing with MS-NIMBLE is sensitive and specific, we
turn our attention to the reliability of MS-NIMBLE’s coefficient estimates for
respiratory-associated metabolites. Since the ground truth is unknown, we
study similar metabolites, as they are likely to have similar effects. Given
our results in Figure 2(b), we consider piperine- and bilirubin-related
metabolites, where we chose the latter because E,Z-bilirubin’s photoisomer
Z,Z-bilirubin was fully observed and shown to be a replicable biomarker for
infant wheeze [1]. Figures 2(d)-(e) provide the results, which illustrate the
consistency of MS-NIMBLE’s estimates. Figure 2(e) is particularly interesting,
as it suggests MS-NIMBLE’s estimates and standard errors for metabolites with
missing data are as reliable as those for fully observed metabolites.

To further explore the fidelity of MS-NIMBLE’s estimates, we compared
estimators for the effect of sex, an important source of metabolite variation,
on metabolite levels in the 0.5 year COPSAC and INSPIRE datasets. Let
$\hat{\beta}_{g}^{(\text{C})},\hat{\beta}_{g}^{(\text{I})}$ be a method’s sex
effect estimates for metabolite $g$ in COPSAC and INSPIRE and
$\hat{\operatorname{\mathbb{V}}}(\cdot)$ their estimated variances. Since
these data were collected from unrelated infants at similar ages, their sex
effects should be the same, meaning the z-score
$\\{\hat{\beta}_{g}^{(\text{C})}\allowbreak-\allowbreak\hat{\beta}_{g}^{(\text{I})}\\}/[\hat{\operatorname{\mathbb{V}}}\\{\hat{\beta}_{g}^{(\text{C})}\\}+\hat{\operatorname{\mathbb{V}}}\\{\hat{\beta}_{g}^{(\text{I})}\\}]^{1/2}$
should be approximately $N(0,1)$. Interestingly, the metabolome-wide z-scores
for imputation-based methods, but not MS-NIMBLE, were significantly inflated
(Table 2), indicating imputation-based estimates and their standard errors are
unreliable. While several factors are likely responsible for this inflation,
we hypothesized errant effect estimates for metabolites with missing data were
partly responsible. Given Section 3 and supplemental Figure S3’s simulation
results showing estimates in trait-associated metabolites are most corrupted
by missing data, we considered z-scores for the 64 sex-associated missing
metabolites, defined as metabolites with missing data and sex q-values $\leq
0.2$ in at least one method, dataset pair. Consistent with our hypothesis,
Table 2 shows these z-scores were inflated in imputation methods, whereas MS-
NIMBLE showed no evidence of inflation. The conclusions were the same even
when we separately examined each method’s sex-associated missing metabolites
(Figure S4), implying differences between MS-NIMBLE and imputation methods
could not be attributed to metabolite selection biases, and indicate MS-
NIMBLE’s estimates and standard errors are accurate.

Table 2: Root mean squared z-score (RMSZ) for all analyzed metabolites (second
row) and the 64 sex-associated missing metabolites (third row), where an RMSZ
$>1$ suggests z-scores are inflated. The p-value is for the null hypothesis
that z-scores are $N(0,1)$.

###  6.2 Metabolite GWAS in the six month COPSAC data

We examined the effect of genotype at 1.4 million SNPs on metabolite levels in
the six month COPSAC data to evaluate the performance of our methodology
proposed in Section 4.4. As far as we are aware, only [3] has considered
controlling for latent sources of variation in mtGWAS studies. However, their
method requires determining a set of latent covariates for each metabolite-SNP
pair, which, as determined by their simulation results, would take 12 CPU
Years if applied to our data. We therefore compared our results to those using
the current state of the art, which involves first imputing missing metabolite
levels and subsequently regressing them onto genotype without considering
latent variation [3]. We present results for minimum imputation, but note
imputation technique did not alter results.

Figure 3(a) contains the results, where the second and third rows imply nearly
all of the genetic effect is idiosyncratic and appears in the error terms
$e_{gi}$, whereas there is no evidence indicating latent factors $\bm{c}_{i}$
mediate genetic effects. This suggests mtGWAS analyses should be performed
conditional on estimated latent factors, as in the second row of Figure 3(a),
which is equivalent to data de-noising. This is recapitulated by Figure 3(b),
which shows such de-noising reduces the residual variance by $\approx 40\%$,
thereby effectively increasing the sample size by 67%.

Figure 3: mtGWAS results for metabolites with missing data. (a): A metabolite
was “significant” if it was associated with at least one SNP at the Bonferroni
p-value threshold $(5\times 10^{-8})/(656+249)$. (b): Reduction in residual
variance after adjusting for latent factors. (c): Named metabolites that were
identified by MS-NIMBLE in the second row of (a), but not minimum imputation.
A metabolite-genomic region association had previous evidence if the region
contained SNPs previously shown to be associated with the metabolite.
Superscripts are m: derived from Metabolon; u: obtained from Uniprot; 1: [16];
2: [17]; 3: [18].

The last row of Figure 3(a) indicates existing approaches are underpowered,
where 11 out of the 13 metabolites identified by minimum imputation were among
the 17 metabolites identified by our proposed method in row two of Figure
3(a). To explore the veracity our method’s results, we sought to evaluate the
biological significance of the four named metabolites uniquely identified by
our method. We did not consider the other two metabolites, since they were
unnamed. Figure 3(c) shows that two out of the four associations have
previously been observed, whereas, to the best of our knowledge, the results
involving 21-hydroxypregnenolone monosulfate and N-linoleoyltaurine are novel.
Critically, their metabolite descriptions and associated gene functions are
congruent, suggesting our method improves power to identify genuine mtGWAS
associations.

##  7 Conclusion

We developed MS-NIMBLE, a rigorous suite of methods to analyze metabolomics
data with non-ignorable missing observations and latent factors that offers
all the practical advantages of missing data imputation. We derived its
theoretical properties and demonstrated its superior performance in
differential abundance and mtGWAS using three real datasets. We believe this
work offers a critical step towards reliable estimation and inference in
metabolomic studies.

## Acknowledgments

This work was supported by the NIH [UL1 TR001857 (C.M.), K01 HL149989 (K.T.),
U19 AI 095227 (T.H.), UG3/UH3 OD023282 (T.H.), UL1 TR002243 (T.H.), R01
HL129735 (C.O.)]. All funding received by COPSAC is listed on www.copsac.com.

## References

  * [1] Kedir N. Turi et al.  “Unconjugated bilirubin is associated with protection from early-life wheeze and childhood asthma”  In _Journal of Allergy and Clinical Immunology_ 148.1, 2021, pp. 128–138
  * [2] 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
  * [3] Apolline Gallois, Joel Mefford, Arthur Ko, Amaury Vaysse, Hanna Julienne, Mika Ala-Korpela, Markku Laakso, Noah Zaitlen, Päivi Pajukanta and Hugues Aschard  “A comprehensive study of metabolite genetics reveals strong pleiotropy and heterogeneity across time and context”  In _Nature Communications_ 10.1, 2019, pp. 4788
  * [4] Runmin Wei, Jingye Wang, Mingming Su, Erik Jia, Shaoqiu Chen, Tianlu Chen and Yan Ni  “Missing Value Imputation Approach for Mass Spectrometry-based Metabolomics Data”  In _Scientific Reports_ 8.1, 2018, pp. 663
  * [5] Jonathan A C Sterne, Ian R White, John B Carlin, Michael Spratt, Patrick Royston, Michael G Kenward, Angela M Wood and James R Carpenter  “Multiple imputation for missing data in epidemiological and clinical research: potential and pitfalls”  In _BMJ (Clinical research ed.)_ 338, 2009, pp. b2393–b2393
  * [6] Andreas Buja and Nermin Eyuboglu  “Remarks on Parallel Analysis”  In _Multivariate Behavioral Research_ 27.4, 1992, pp. 509–540
  * [7] Jasmit Shah, Guy N. Brock and Jeremy Gaskins  “BayesMetab: treatment of missing values in metabolomic studies using a Bayesian modeling approach”  In _BMC Bioinformatics_ 20.S24, 2019  DOI: 10.1186/s12859-019-3250-2
  * [8] 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
  * [9] Jae Kwang Kim and Cindy Long Yu  “A Semiparametric Estimation of Mean Functionals With Nonignorable Missing Data”  In _Journal of the American Statistical Association_ 106.493, 2011, pp. 157–165  DOI: 10.1198/jasa.2011.tm10104
  * [10] Qingyuan Zhao, Jingshu Wang, Gibran Hemani, Jack Bowden and Dylan S. Small  “Statistical inference in two-sample summary-data Mendelian randomization using robust adjusted profile score”  In _The Annals of Statistics_ 48.3 Institute of Mathematical Statistics, 2020  DOI: 10.1214/19-aos1866
  * [11] 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
  * [12] H. Bisgaard et al.  “Deep phenotyping of the unselected COPSAC2010 birth cohort study”  In _Clinical &amp$\operatorname{;}$ Experimental Allergy_ 43.12, 2013, pp. 1384–1394  DOI: 10.1111/cea.12213
  * [13] Emma K Larkin et al.  “Objectives, design and enrollment results from the Infant Susceptibility to Pulmonary Infections and Asthma Following RSV Exposure Study (INSPIRE)”  In _BMC Pulmonary Medicine_ 15.1, 2015
  * [14] Bruce Gillis, Igor M. Gavin, Zarema Arbieva, Stephen T. King, Sundararajan Jayaraman and Bellur S. Prabhakar  “Identification of human cell responses to benzene and benzene metabolites”  In _Genomics_ 90.3 Elsevier BV, 2007, pp. 324–333  DOI: 10.1016/j.ygeno.2007.05.003
  * [15] Rachel Kelly, Michael McGeachie, Kathleen Lee-Sarwar, Priyadarshini Kachroo, Su Chu, Yamini Virkud, Mengna Huang, Augusto Litonjua, Scott Weiss and Jessica Lasky-Su  “Partial Least Squares Discriminant Analysis and Bayesian Networks for Metabolomic Prediction of Childhood Asthma”  In _Metabolites_ 8.4, 2018, pp. 68  DOI: 10.3390/metabo8040068
  * [16] Trisha J. Grevengoed et al.  “N-acyl taurines are endogenous lipid messengers that improve glucose homeostasis”  In _Proceedings of the National Academy of Sciences_ 116.49, 2019, pp. 24770–24778
  * [17] Pirro G. Hysi, Massimo Mangino, Paraskevi Christofidou, Mario Falchi, Edward D. Karoly, NIHR Bioresource Investigators, Robert P. Mohney, Ana M. Valdes, Tim D. Spector and Cristina Menni  “Metabolome Genome-Wide Association Study Identifies 74 Novel Genomic Regions Influencing Plasma Metabolites Levels”  In _Metabolites_ 12.1, 2022  DOI: 10.3390/metabo12010061
  * [18] Natalja Kurbatova et al.  “Urinary metabolic phenotyping for Alzheimer’s disease”  In _Scientific reports_ 10.1, 2020, pp. 21745–21745
  * [19] Christa E Müller and Kenneth A Jacobson  “Xanthines as adenosine receptor antagonists”  In _Handbook of experimental pharmacology_ , 2011, pp. 151–199
  * [20] Raphael Guzman, Darı́o Echeverri, Félix R. Montes, Mariana Cabrera, Angélica Galán and Angélica Prieto  “Caffeine’s Vascular Mechanisms of Action”  In _International Journal of Vascular Medicine_ 2010, 2010, pp. 834060
  * [21] L. Prieto, V. Gutiérrez, J. Liñana and J. Marín  “Bronchoconstriction induced by inhaled adenosine 5’-monophosphate in subjects with allergic rhinitis”  In _European Respiratory Journal_ 17.1 European Respiratory Society, 2001, pp. 64–70
  * [22] Laurent Y Rios, Marie-Paule Gonthier, Christian Rémésy, Isabelle Mila, Catherine Lapierre, Sheryl A Lazarus, Gary Williamson and Augustin Scalbert  “Chocolate intake increases urinary excretion of polyphenol-derived phenolic acids in healthy human subjects”  In _The American Journal of Clinical Nutrition_ 77.4, 2003, pp. 912–918  DOI: 10.1093/ajcn/77.4.912
  * [23] John D. Storey  “A direct approach to false discovery rates”  In _Journal of the Royal Statistical Society: Series B_ 63.3, 2001, pp. 479–498
  * [24] Chris McKennan  “Factor analysis in high dimensional biological data with dependent observations”, 2020  eprint: arXiv:2009.11134
  * [25] Chris McKennan and Dan Nicolae  “Estimating and Accounting for Unobserved Covariates in High-Dimensional Correlated Data”  In _Journal of the American Statistical Association_ , 2020, pp. 1–12
  * [26] Gitta Kutyniok, Yonina C. Eldar and Roman Vershynin  “Introduction to the non-asymptotic analysis of random matrices”  In _Compressed Sensing: Theory and Applications_ Cambridge: Cambridge University Press, 2012, pp. 210–268  DOI: 10.1017/CBO9780511794308.006
  * [27] Rafał Latała  “Some Estimates of Norms of Random Matrices”  In _Proceedings of the American Mathematical Society_ 133.5 American Mathematical Society, 2005, pp. 1273–1282  URL: http://www.jstor.org/stable/4097777

Supplemental material for “From differential abundance to mtGWAS: accurate and
scalable methodology for metabolomics data with non-ignorable missing
observations and latent factors”

##  S1 The normality assumption

We rely on the assumption that $e_{gi}$ in (2.1) is approximately normally
distributed to develop statistically efficient estimators. However, while we
assume $e_{gi}$ is approximately normal, we do not assume $y_{gi}$ is normal.
This is a critical distinction, as Figure S1 indicates the latent factors
$\bm{c}_{i}$ may be highly skewed.

Figure S1: Normality of plasma metabolite levels from [1]. (a) Normal Q-Q
plots for the first three estimated latent factors. (b) Normal Q-Q plots for
four randomly chosen metabolites. (c) Q-Q plots for the same four metabolites,
except after regressing out the $K=19$ latent factors.

##  S2 Simulations

###  S2.1 Simulation setup

We simulated 50 datasets containing $p=1200$ metabolites measured in $n=600$
samples with missing observations and $K=10$ latent factors to best mirror our
real data from Section 6. We partitioned individuals into equal sized
treatment and control groups, where the covariate of interest
$\bm{X}\in\\{0,1\\}^{n}$ denotes treatment status. For some constant
$a\in\mathbb{R}$ controlling the dependence of latent factors on $\bm{X}$,
metabolite levels $y_{gi}$ and missingness indicators $r_{gi}$ were then
simulated according to (S2.1) below.

$\displaystyle\log\left(\alpha_{g}\right)\sim
N_{1}\left(\mu_{\alpha},0.4^{2}\right),\quad\delta_{g}\sim
N_{1}\left(16,1.2^{2}\right),\quad g\in[p]$ (S2.1a)
$\displaystyle\bm{C}=\left(\bm{c}_{1}\cdots\bm{c}_{n}\right)^{\top}\sim
MN_{n\times
K}\left(\left(a\bm{X},\,a\bm{X},\,\bm{0}_{n}\cdots\bm{0}_{n}\right),I_{n},I_{K}\right)$
(S2.1b)
$\displaystyle\bm{\ell}_{g_{k}}\sim\pi_{k}\delta_{0}+\left(1-\pi_{k}\right)N_{1}\left(0,\tau_{k}^{2}\right),\quad
g\in[p];k\in[K]$ (S2.1c) $\displaystyle\mu_{g}\sim
N_{1}\left(18,5^{2}\right),\quad\sigma_{g}^{2}\sim\text{Gamma}\left(0.2^{-2},0.2^{-2}\right),\quad
g\in[p]$ (S2.1d) $\displaystyle\beta_{g}\sim
0.8\delta_{0}+0.2N_{1}\left(0,0.4^{2}\right),\quad g\in[p]$ (S2.1e)
$\displaystyle y_{gi}\sim
N_{1}\left(\mu_{g}+\bm{X}_{i}\beta_{g}+\bm{c}_{i}^{\top}\bm{\ell}_{g},\sigma_{g}^{2}\right),\quad
g\in[p];i\in[n]$ (S2.1f) $\displaystyle
r_{gi}\sim\text{Bernoulli}\left[\tilde{\Psi}\left\\{\alpha_{g}\left(y_{gi}-\delta_{g}\right)\right\\}\right],\quad
g\in[p];i\in[n]$ (S2.1g)

where $\delta_{0}$ is the point mass at 0 and $\mu_{\alpha}$ in (S2.1a) was
set so that if $Z$ has cumulative distribution function
$\tilde{\Psi}\left\\{\exp\left(\mu_{\alpha}\right)x\right\\}$,
$\operatorname{\mathbb{V}}(Z)=1$. To study scenarios where we incorrectly
specify $\Psi$ in (2.2), we let $\tilde{\Psi}$ in (S2.1g) be the cumulative
distribution function (CDF) of a logistic random variable, but analyzed the
data assuming $\Psi$ was the CDF of a t-distribution with four degrees of
freedom. The normal means and variances in (S2.1a) and (S2.1d) were chosen to
match those estimated in the three datasets from Section 6, and the parameters
used to simulate the loadings $\bm{\ell}_{g}$ in (S2.1c) are given in Table
S1. The loadings were chosen so that the eigenvalues
$\lambda_{1},\ldots,\lambda_{K}$ from Assumption 5.1 ranged from
$n^{-0.47}=0.05$ to 0.80 on average, which mirrored the eigenvalues estimated
from the six year COPSAC data (see Table 1). The constant $a$ in (S2.1b) was
chosen so that $\bm{C}$ explained 60% of the variance in $\bm{X}$ on average,
and was chosen to match the substantial correlation between latent factors and
infant wheeze in INSPIRE (see Figure 2(a)). Lastly, we simulated the effects
of interest $\beta_{g}$ in (S2.1e) to violate the sparsity assumption in (i)
of Theorem 5.2, which is also used to prove Theorem 5.3.

Table S1: The $\pi_{k}$ and $\tau_{k}$ values used to simulate
$\bm{\ell}_{1},\ldots,\bm{\ell}_{p}$ ($k=1,\ldots,10$).

Factor number ($k$) | $1$ | $2$ | $3$ | $4$ | $5$ | $6$ | $7$ | $8$ | $9$ | $10$
---|---|---|---|---|---|---|---|---|---|---
$\pi_{k}$ | 0 | 0 | 0.80 | 0.60 | 0.50 | 0.35 | 0.30 | 0.20 | 0.20 | 0.20
$\tau_{k}$ | 0.80 | 0.60 | 0.5 | 0.5 | 0.5 | 0.5 | 0.5 | 0.5 | 0.5 | 0.5

###  S2.2 Simulation results

We compared MS-NIMBLE’s estimates for $\beta_{g}$ to those from MetabMiss [2]
and imputation-based methods, the latter of which first impute missing data
with one of minimum imputation, singular value decomposition (SVD), K-nearest
neighbors (KNN), or random forest (RF), and subsequently use CATE [11] to
estimate latent factors. While other methods are capable of estimating latent
factors in complete data, we found that CATE gave the best results. Imputation
hyperparameters were $K=10$ factors for SVD and the software defaults
recommended in [4] for KNN and RF. The estimates for $\alpha_{g}$ and
$\delta_{g}$, which were used by both MS-NIMBLE and MetabMiss, were obtained
using the method proposed in [2] and outlined in Section 4.1 with 5 potential
instruments. We do not include results when $\bm{C}$ is known or when it is
ignored, as they both performed similarly to and uniformly worse than KNN
imputation, respectively.

On the average, 485 metabolites were fully observed (i.e. missing in $<5\%$ of
samples) and 300 were missing (i.e. $\geq 5\%$ but $\leq 50\%$ missing data).
Metabolites with $>50\%$ missing data were discarded. We first consider each
method’s ability to identify missing metabolites with non-zero $\beta_{g}$.
Figure S2 gives the results, where Figure S2(a) indicates MS-NIMBLE and, to a
lesser extent, MetabMiss are able to control false discovery rates at their
nominal levels. However, Figure S2(b) indicates MS-NIMBLE has 50% greater
power than MetabMiss to identify treatment-related metabolites with missing
data. These results are consistent with the fact that while MetabMiss does use
inverse probability weighting to account for the non-ignorable missing data,
their estimates for $\beta_{g}$ discard missing data, and are therefore less
powerful. On the other hand, imputation-based methods inflate error rates and
have poor power. The former is consistent with our discussion from Section 3,
as their false discovery proportions resembled nominal levels when we
simulated data with latent factors $\bm{C}$ that did not depend on $\bm{X}$.

Figure S2: False discovery proportion (a) and power (b) for metabolites with
missing data at q-value thresholds of 0.05 (left) and 0.2 (right). The dashed
red and blue lines indicate the q-value thresholds and MS-NIMBLE’s median
power, respectively.

We lastly considered each method’s estimates and 95% confidence intervals for
$\beta_{g}$ for metabolites $g$ with missing data, where confidence intervals
were standard Wald intervals assuming estimates for $\beta_{g}$ were
approximately normal. Figure S3(a) contains the results, where only MS-NIMBLE
and MetabMiss return accurate intervals. However, consistent with the above
discussion and results from Section 6.1, MetabMiss’s intervals are on average
over 25% wider than MS-NIMBLE’s (Figure S3(b)). We also see that imputation-
based intervals become less accurate as $\lvert\beta_{g}\rvert$ increases,
which corroborates our discussion in Section 3.

Figure S3: (a): 95% confidence interval coverage for metabolites with missing
data. The dashed grey line indicates 95% coverage. (b): 95% confidence
interval widths for metabolites with missing data. Each point represents a
simulated metabolite with missing data. MS-NIMBLE’s and MetabMiss’s confidence
interval widths did not depend on $\beta_{g}$.

##  S3 Additional real data details and results from Section 6

###  S3.1 Additional real data and analysis details

Raw metabolite intensities were log base 2-transformed. There were no
additional quality control or pre-processing steps.

Missingness mechanism parameters $\alpha_{g},\delta_{g}$, which are used by
both MS-NIMBLE and MetabMiss, were estimated using the procedure outlined in
Section 4.1 and described in detail in [2] with 10 potential instruments.
Missing data were imputed exactly as described in Section S2.2.

Differential abundance regressions in INSPIRE were performed by controlling
for the observed covariates daycare status (yes/no), breast-feeding status
(exclusively breast-fed or not in the first six months of life), age in
months, and sex in the first year wheeze analysis. The sRAW analysis in the
six year COPSAC dataset was done conditional on sex, and we did not include
any observed nuisance covariates in the 0.5 year COPSAC sex regression.

###  S3.2 Additional real data results

We first justify the observed relationship between infant wheeze and
theobromine and 3-(3-hydroxyphenyl)propionate levels in INSPIRE (see Figure
2(c)), where wheezers tended to have higher plasma concentrations of both
metabolites. Theobromine is an alkaloid commonly found in the cacao plant, and
is a notable adenosine receptor antagonist [19]. Higher theobromine
concentrations tend to increase plasma adenosine levels [20], thereby
potentially exacerbating adenosine’s bronchoconstricting properties [20, 21].
The metabolite 3-(3-hydroxyphenyl)propionate is a phenolic degradation product
of proanthocyanidins, the most abundant polyphenols present in chocolate [22],
and therefore may simply correlate with infant wheeze because it correlates
with theobromine levels.

We next consider the sex-related z-scores defined in Section 6.1. To argue
that the inter-method differences in root mean squared z-scores for the 64-sex
related metabolites was not due to metabolite selection bias (i.e. winner’s
curse), we investigated each method’s sex-associated metabolites with missing
data. The results are given in Figure S4, and show that only MS-NIMBLE’s
z-scores show no evidence of inflation. This suggests that differences between
MS-NIMBLE and imputation methods in Table 2 cannot be attributed to metabolite
selection bias.

We lastly consider MS-NIMBLE’s computation time. The most computationally
demanding component in differential abundance analyses is estimating each the
missingness mechanism parameters $\alpha_{g},\delta_{g}$ (see Section 4.1),
which took 40 minutes for the 0.5 year COSAC dataset (the dataset with the
largest sample size). However, this only needed to be computed once, and was
stored for use in all downstream analyses. The subsequent sex analysis in the
0.5 year COSAC dataset took 3.4 minutes.

Figure S4: Q-Q plot for each method’s sex-associated missing metabolites. The
z-score is as defined in Section 6.1 and a metabolite was sex-associated in a
method if it (i) was analyzed in both the six month COPSAC and INSPIRE
datasets, (ii) contained missing data in at least one dataset, and (iii) had a
q-value $\leq 0.2$ using that method. The statistic $\sigma_{z}$ is the
method’s root mean squared z-score for their sex-associated missing
metabolites.

##  S4 Refining our estimator for $\bm{\Omega}$

Here we provide a way to refine our estimator for $\bm{\Omega}$ in (4.3) that
iteratively removes “outlying” metabolites that likely depend on the
covariate(s) of interest. It should be noted that this is our software default
estimator for $\bm{\Omega}$.

Briefly, assume $\bm{X}$ can be written as $\bm{X}=[\bm{X}_{I},\bm{X}_{N}]$,
where $\bm{X}_{I}\in\mathbb{R}^{n\times d_{I}}$ contains the $d_{I}$
covariates of interest and $\bm{X}_{N}$ contains the remaining nuisance
covariates. Let
$\bm{\beta}_{g,I},\allowbreak\hat{\bm{\beta}}_{g,I}^{\operatorname{(IPW)}}\in\mathbb{R}^{d_{I}}$
be the first $d_{I}$ elements of $\bm{\beta}_{g}$ and the inverse probability
weighted estimator $\hat{\bm{\theta}}_{g}^{\operatorname{(IPW)}}$ defined in
(4.5), respectively. Then Lemma S8.20 shows that under the same assumptions
used to prove Theorem 5.3,
$\lVert\hat{\bm{\beta}}_{g,I}^{\operatorname{(IPW)}}-\bm{\beta}_{g,I}\rVert_{2}=O_{P}(n^{-1/2})$
if $d_{I}\leq d_{1}$. Further, it is straightforward to extend Lemma S8.20 to
show that for

$\displaystyle\hat{\operatorname{\mathbb{V}}}\\{\hat{\bm{\beta}}_{g,I}^{\operatorname{(IPW)}}\\}=\left((\textstyle\sum_{i=1}^{n}\hat{w}_{gi}\hat{\bm{z}}_{i}\hat{\bm{z}}_{i}^{\top})^{-1}[\textstyle\sum_{i=1}^{n}\hat{w}_{gi}^{2}\\{y_{gi}-\hat{\bm{z}}_{i}^{\top}\bm{\theta}_{g}^{\operatorname{(IPW)}}\\}^{2}\hat{\bm{z}}_{i}\hat{\bm{z}}_{i}^{\top}](\textstyle\sum_{i=1}^{n}\hat{w}_{gi}\hat{\bm{z}}_{i}\hat{\bm{z}}_{i}^{\top})^{-1}\right)_{1:d_{I},1:d_{I}}$

the sandwich estimator for
$\operatorname{\mathbb{V}}\\{\hat{\bm{\beta}}_{g,I}^{\operatorname{(IPW)}}\\}$,
$x_{g}^{2}=\\{\hat{\bm{\beta}}_{g,I}^{\operatorname{(IPW)}}\\}^{\top}[\hat{\operatorname{\mathbb{V}}}\\{\hat{\bm{\beta}}_{g,I}^{\operatorname{(IPW)}}\\}]^{-1}\hat{\bm{\beta}}_{g,I}^{\operatorname{(IPW)}}\allowbreak\operatorname{\stackrel{{\scriptstyle\text{d}}}{{\to}}}\allowbreak\chi_{d_{I}}^{2}$
under the null hypothesis $H_{0,g}:\bm{\beta}_{g,I}=0$. To refine our estimate
for $\bm{\Omega}$, we compute p-values for $H_{0,g}$ by comparing $x_{g}^{2}$
to the upper quantiles of a $\chi_{d_{I}}^{2}$, use [23] to subsequently
determine q-values, and re-estimate $\bm{\Omega}$ using the regression in
(4.3) after removing metabolites from said regression whose q-values fall
below a user-specified threshold $q$. Our software default is to let $q=0.1$
and iterate this procedure 3 times.

##  S5 Extensions when some metabolites have fully observed data

###  S5.1 Methodological extensions

The factor analysis- and mtGWAS-related estimators are the only estimators
that need to be updated to allow fully observed metabolites. For the former,
we simply let $\hat{w}_{gi}$ in (4.2) be 1 if metabolite $g$ has no missing
data. For the mtGWAS estimators described in Section 4.4, we regress $y_{gi}$
onto genotype $G_{si}$ and estimated latent factors $\hat{\bm{c}}_{i}$ to
estimate $\gamma_{sg}^{(e)}$ and the estimator’s variance. We then use
standard Wald-based inference to test $H_{0,sg}^{(e)}$. Testing
$H_{0,sg}^{(c)}$ remains unchanged. Since the test statistics used to test
$H_{0,sg}^{(e)}$ and $H_{0,sg}^{(c)}$ are asymptotically $\chi_{1}^{2}$ and
independent under Assumptions 5.1 and 5.2, we simply add the test statistics
and compare it to the upper quantiles of a $\chi_{2}^{2}$ to test
$H_{0,sg}^{(c,e)}$.

###  S5.2 Theoretical extensions

The only theoretical extension we must consider is choosing an appropriate
starting point for the estimator $\hat{\mathcal{P}}$ from Theorem 5.1, which
is discussed in Remark 5.3. Let $\mathcal{O}\subset[p]$ be the set of
metabolites with fully observed data,
$\lambda^{(\mathcal{O})}_{1}\geq\cdots\geq\lambda^{(\mathcal{O})}_{K}$ be the
eigenvalues of $\sum_{g\in\mathcal{O}}\bm{\ell}_{g}\bm{\ell}_{g}^{\top}$,
$\hat{\bm{V}}\in\mathbb{R}^{n\times K}$ be the first $K$ right singular
vectors of
$[y_{gi}]_{g\in\mathcal{O};i\in[n]}P_{\bm{X}}^{\perp}\in\mathbb{R}^{\lvert\mathcal{O}\rvert\times
n}$, and define
$\hat{\mathcal{P}}^{(\mathcal{O})}=\hat{\bm{V}}\hat{\bm{V}}^{\top}$. Then
under Assumptions 5.1 and 5.2, the proof of Theorem 4 in [24] can easily be
used to show that
$\lVert\hat{\mathcal{P}}^{(\mathcal{O})}-\mathcal{P}\rVert_{F}^{2}=O_{P}[\\{\lambda^{(\mathcal{O})}_{K}n\\}^{-1}]$.
Therefore, if $\lambda^{(\mathcal{O})}_{K}\gtrsim n^{-1+\epsilon}$ for any
$\epsilon>0$, Corollary S8.3 in Section S8.4 implies $\hat{\mathcal{P}}$ will
satisfy the condition $\lVert\hat{\mathcal{P}}-\mathcal{P}\rVert_{F}\leq\eta$
in Theorem 5.1 when we solve the optimization in (4.2) by initializing
$\bm{C}_{\perp}=\hat{\bm{V}}$ and iteratively updating
$\\{\bm{b}_{g},\bm{\ell}_{g}\\}_{g\in[p]}$ and $\bm{C}_{\perp}$.

##  S6 Outline and notation for the rest of the supplement

###  S6.1 Outline for the remaining supplement

The rest of the supplement is devoted to proving the theoretical statements
made in Sections 3 and 5. Due to its length, we give a compendious outline
below.

  * •

Section S7: we provide the regularity conditions for and prove Proposition 3.1
stated in Section 3.

  * •

Section S8: we prove Theorem 5.1, Theorem 5.2, Proposition 5.2, and Theorem
5.3. The proofs can be found in:

    * –

Theorem 5.1: Corollary S8.4 in Section S8.4.

    * –

Theorem 5.2: Corollary S8.8 in Section S8.5.

    * –

Proposition 5.2: A direct consequence of Lemma S8.22 in Section S8.6. See
Remark S8.10.

    * –

Theorem 5.3: proven in Theorem S8.3 in Section S8.6.

  * •

Section S9: we extend our mtGWAS test statistics to allow $\bm{x}_{i}\neq 0$,
prove an extension of Theorem 5.4 that allows $\bm{x}_{i}\neq 0$, and
illustrate the computational efficiency of our mtGWAS test statistics.

###  S6.2 Notation

For any matrix $\bm{M}\in\mathbb{R}^{m\times n}$, we define
$\bm{M}_{i\mathchoice{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\displaystyle*$}}}}}{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\textstyle*$}}}}}{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\scriptstyle*$}}}}}{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\scriptscriptstyle*$}}}}}}\in\mathbb{R}^{n}$,
$\bm{M}_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\displaystyle*$}}}}}{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\textstyle*$}}}}}{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\scriptstyle*$}}}}}{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\scriptscriptstyle*$}}}}}j}\in\mathbb{R}^{m}$,
and $\bm{M}_{ij}\in\mathbb{R}$ to be the $i$th row, $j$th column, and
$(i,j)$th element of $\bm{M}$, respectively. We also define
$P_{\bm{M}},P_{\bm{M}}^{\perp}\in\mathbb{R}^{n\times n}$ to be the orthogonal
projections matrices that project vectors onto the image of $\bm{M}$ and
kernel of $\bm{M}^{\top}$. Let $\\{\bm{X}_{n}\\}_{n\geq 1}$ be a sequence of
random vectors or matrices. Unless otherwise specified,
$\bm{X}_{n}=O_{p}(a_{n})$ if $\lVert\bm{X}_{n}\rVert_{2}/a_{n}=O_{P}(1)$ and
$\bm{X}_{n}=o_{p}(a_{n})$ if $\lVert\bm{X}_{n}\rVert_{2}/a_{n}=o_{P}(1)$ as
$n\to\infty$. Lastly, for random vector $\bm{e}$, we use the notation
$\bm{e}\sim(\bm{\mu},\bm{V})$ if $\operatorname{\mathbb{E}}(\bm{e})=\bm{\mu}$
and $\operatorname{\mathbb{V}}(\bm{e})=\bm{V}$.

##  S7 Proof of Proposition 3.1

We first state the complete set of sufficient conditions needed to prove
Proposition 3.1.

######  Assumption S7.3 (Proposition 3.1).

In addition to the assumptions in the statement of Proposition 3.1, assume the
following hold:

  1. (a)

The elements of $\\{x_{i},\bm{c}_{i}\\}_{i\in[n]}$ are independent and
identically distributed and are independent of $\\{e_{gi}\\}_{i\in[n]}$.

  2. (b)

$\operatorname{\mathbb{E}}(x_{i}^{4})\leq c$ for some constant $c>0$ and
$\operatorname{\mathbb{E}}(\lvert\bm{c}_{i_{k}}\rvert^{m})\leq c_{m}$ for all
$k\in[K]$, $m>0$, and some constants $c_{m}>0$.

  3. (c)

$\alpha_{g}>0$.

######  Remark S7.6.

The moment assumption on $\bm{c}_{i}$ is the same as that in Assumption 5.1.

###### Proof of Proposition 3.1.

We drop the subscript $g$ to simplify notation. Since
$\\{x_{i},\bm{c}_{i}\\}_{i\in[n]}$ are identically distributed and our design
matrix includes the intercept, it suffices to assume
$\operatorname{\mathbb{E}}(x_{i})$ and $\operatorname{\mathbb{E}}(\bm{c}_{i})$
are 0. Let
$m=a\min_{i:r_{i}=0}y_{i}=a\mu+a\min_{i:r_{i}=0}(\bm{\ell}^{\top}\bm{c}_{i}+e_{i})$.
Since $e_{i}$ is sub-Gaussian and by the moment assumptions on $\bm{c}_{i}$,
$\left\lvert m\right\rvert=O_{P}(n^{\epsilon})$ for any $\epsilon>0$. For any
$M>0$, the Gaussian assumption on $e_{i}$ and the moment assumptions on
$\bm{c}_{i}$ also imply
$\operatorname{Pr}\\{y_{i}\in(-2M,-M)\\}=\delta_{1,M}>0$. Since
$\operatorname{Pr}\\{r_{i}=1\mid
y_{i}\in(-2M,-M)\\}\geq\operatorname{Pr}\\{r_{i}=1\mid
y_{i}=-2M\\}=\delta_{2,M}>0$, this implies the event
$\\{y_{i}\in(-2M,-M),r_{i}=1\\}$ occurs infinitely often as $n\to\infty$,
meaning $m\to-\infty$ as $n\to\infty$.

We consider the cases $\bm{\ell}=0$ and $x_{i}$ is independent of $\bm{c}_{i}$
separately. Suppose first that $\bm{\ell}=0$. Then
$\bm{z}_{i}=(x_{i},\bm{c}_{i}^{\top})^{\top}$ is independent of $y_{i}$ and
the elements of $\\{\bm{z}_{i},y_{i}\\}_{i\in[n]}$ are independent and
identically distributed. Let $\bm{y}=(y_{1},\ldots,y_{n})$,
$\bm{V}=\operatorname{\mathbb{V}}(\bm{z}_{i})$,
$\bm{R}=\operatorname{diag}(r_{1},\ldots,r_{n})$,
$\bm{Z}=(\bm{z}_{1}\cdots\bm{z}_{n})^{\top}$, and
$\bm{y}_{I}=\bm{R}\bm{y}+m(I_{n}-\bm{R})\bm{1}$ be the imputed data. Then for
$\bm{e}_{1}\in\\{0,1\\}^{K+1}$ the first standard basis vector,

$\displaystyle n^{1/2}\hat{\beta}$
$\displaystyle=\bm{e}_{1}^{\top}(\bm{V}^{-1}\hat{\bm{V}})^{-1}\bm{V}^{-1}(n^{-1/2}\bm{Z}^{\top}P_{\bm{1}}^{\perp}\bm{y}_{I}),\quad\hat{\bm{V}}=n^{-1}\bm{Z}^{\top}P_{\bm{1}}^{\perp}\bm{Z}$
$\displaystyle n\hat{s}^{2}$
$\displaystyle=\hat{\sigma}^{2}\bm{e}_{1}^{\top}\hat{\bm{V}}^{-1}\bm{e}_{1},\quad\hat{\sigma}^{2}=(n-K-2)^{-1}\bm{y}_{I}^{\top}P_{[\bm{1},\bm{Z}]}^{\perp}\bm{y}_{I}.$

Since $\lVert\bm{V}^{-1}\hat{\bm{V}})^{-1}-I_{K+1}\rVert_{2}=o_{P}(1)$, we
need only show that for
$\bm{v}=\bm{Z}\bm{V}^{-1}\bm{e}_{1}(\bm{e}_{1}^{\top}\bm{V}^{-1}\bm{e}_{1})^{-1/2}$,

$\displaystyle(\hat{\sigma}^{2})^{-1/2}(n^{-1/2}\bm{v}^{\top}P_{\bm{1}}^{\perp}\bm{y}_{I})\to
N(0,1).$

We start by studying $\hat{\sigma}^{2}$. First, it is easy to see that because
$\lvert m\rvert\to\infty$,
$(n-K-2)^{-1}\bm{y}_{I}^{\top}P_{\bm{1}}^{\perp}\bm{y}_{I}=m^{2}\\{c+o_{P}(1)\\}$
for some $c>0$. Next, for $\tilde{n}=n-K-2$,

$\displaystyle\hat{\sigma}^{2}=\tilde{n}^{-1}\bm{y}_{I}^{\top}P_{\bm{1}}^{\perp}\bm{y}_{I}-\\{\tilde{n}^{-1}\bm{y}_{I}^{\top}(\bm{Z}-\bm{1}\bar{\bm{z}}^{\top})\\}\hat{\bm{V}}^{-1}\\{n^{-1}(\bm{Z}-\bm{1}\bar{\bm{z}}^{\top})^{\top}\bm{y}_{I}\\},\quad\bar{\bm{z}}=n^{-1}\bm{Z}^{\top}\bm{1},$

where because $\bm{Z}$ is independent of $\bm{y}_{I}$ and
$\bar{\bm{z}}=O_{P}(n^{-1/2})$,

$\displaystyle\\{\tilde{n}^{-1}\bm{y}_{I}^{\top}(\bm{Z}-\bm{1}\bar{\bm{z}}^{\top})\\}\hat{\bm{V}}^{-1}\\{n^{-1}(\bm{Z}-\bm{1}\bar{\bm{z}}^{\top})^{\top}\bm{y}_{I}\\}=$
$\displaystyle O_{P}\\{\lVert
n^{-1}\bm{y}_{I}^{\top}\bm{Z}\rVert_{2}^{2}+\lVert
n^{-1}\bm{y}_{I}^{\top}\bm{1}\bar{\bm{z}}^{\top}\rVert_{2}\\}$
$\displaystyle=$ $\displaystyle O_{P}(\lvert m\rvert n^{-1/2})=o_{P}(1).$

Since the entries of $\bm{v}$ are mean 0, variance 1, independent, and
independent of $\bm{y}_{I}$, $\lVert
n^{-1/2}\bm{v}^{\top}P_{\bm{1}}^{\top}\bm{y}_{I}\rVert_{2}=O_{P}(\lvert
m\rvert)$, meaning

$\displaystyle(\hat{\sigma}^{2})^{-1/2}(n^{-1/2}\bm{v}^{\top}P_{\bm{1}}^{\perp}\bm{y}_{I})=\tilde{\sigma}^{-1}(n^{-1/2}\bm{v}^{\top}\tilde{\bm{y}}_{I})+o_{P}(1),\quad\tilde{\bm{y}}_{I}=P_{\bm{1}}^{\perp}\bm{y}_{I},\quad\tilde{\sigma}^{2}=\tilde{n}^{-1}\bm{y}_{I}^{\top}P_{\bm{1}}^{\perp}\bm{y}_{I}.$

We therefore need only show
$\tilde{\sigma}^{-1}(n^{-1/2}\bm{v}^{\top}\tilde{\bm{y}}_{I})\to N(0,1)$. To
prove this, we note that
$\operatorname{\mathbb{E}}\\{\tilde{\sigma}^{-1}(n^{-1/2}\bm{v}^{\top}\tilde{\bm{y}}_{I})\mid\bm{y}\\}=0$,
$\operatorname{\mathbb{V}}\\{\tilde{\sigma}^{-1}(n^{-1/2}\bm{v}^{\top}\tilde{\bm{y}}_{I})\mid\bm{y}\\}=1$,
the elements of $\bm{v}$ are independent conditional on $\bm{y}$, and

$\displaystyle
n^{-2}\sum_{i=1}^{n}\operatorname{\mathbb{E}}(\tilde{\sigma}^{-4}\bm{v}_{i}^{4}\tilde{\bm{y}}_{I_{i}}^{4}\mid\bm{y})\leq
c\tilde{\sigma}^{-4}n^{-1}(\max_{i\in[n]}y_{i}^{4})=o_{P}(1)$

for some constant $c>0$. The asymptotic normality of
$\tilde{\sigma}^{-1}(n^{-1/2}\bm{v}^{\top}\tilde{\bm{y}}_{I})$ follows by the
Lindeberg central limit theorem.

We lastly consider the case when $x_{i}$ is independent of $\bm{c}_{i}$. Here,
$\bm{\ell}$ may not be 0, so $y_{i}$ and $\bm{c}_{i}$ may be dependent. Let
$\bm{x}=(x_{1},\ldots,x_{n})^{\top}$. Let
$v=\operatorname{\mathbb{V}}(x_{i})$. We can express $\hat{\beta}/\hat{s}$ as

$\displaystyle\hat{\beta}/\hat{s}=(n^{-1}\bm{x}^{\top}P_{[\bm{1},\bm{C}]}^{\perp}\bm{x}v^{-1})^{-1/2}\hat{\sigma}^{-1}(n^{-1/2}\tilde{\bm{x}}^{\top}P_{[\bm{1},\bm{C}]}^{\perp}\bm{y}_{I}),\quad\tilde{\bm{x}}=v^{-1/2}\bm{x}.$

Since $n^{-1}\bm{x}^{\top}P_{[\bm{1},\bm{C}]}^{\perp}\bm{x}v^{-1}=1+o_{P}(1)$,
it suffices to show
$\hat{\sigma}^{-1}(n^{-1/2}\tilde{\bm{x}}^{\top}P_{[\bm{1},\bm{C}]}^{\perp}\bm{y}_{I})\operatorname{\stackrel{{\scriptstyle\text{d}}}{{\to}}}N(0,1)$
to complete the proof. The above proof of the asymptotic normality when
$\bm{\ell}=0$ implies this will be true if
$\hat{\sigma}^{2}=m^{2}\\{c+o_{P}(1)\\}$ for some constant $c>0$ and if
$\hat{\sigma}^{2}=\tilde{n}^{-1}\bm{y}_{I}^{\top}P_{[\bm{1},\bm{C}]}^{\perp}\bm{y}_{I}+o_{P}(1)$.
For the former, we see that for $\tilde{\bm{z}}_{i}=\\{1\oplus
v^{-1/2}\oplus\operatorname{\mathbb{V}}(\bm{c}_{i})^{-1/2}\\}(1,\bm{x}_{i},\bm{c}_{i}^{\top})^{\top}$,

$\displaystyle\hat{\sigma}^{2}=O_{P}(m)+m^{2}[\operatorname{\mathbb{E}}(1-r_{1})-\operatorname{\mathbb{E}}\\{(1-r_{1})\tilde{\bm{z}}_{1}\\}^{\top}\operatorname{\mathbb{E}}\\{(1-r_{1})\tilde{\bm{z}}_{1}\\}].$

For any non-random unit vector $\bm{u}\in\mathbb{R}^{K+2}$, Holder’s
inequality implies

$\displaystyle\bm{u}^{\top}[\operatorname{\mathbb{E}}\\{(1-r_{1})\tilde{\bm{z}}_{1}\\}\operatorname{\mathbb{E}}\\{(1-r_{1})\tilde{\bm{z}}_{1}\\}^{\top}]\bm{u}$
$\displaystyle=[\operatorname{\mathbb{E}}\\{(1-r_{1})(\tilde{\bm{z}}_{1}^{\top}\bm{u})\\}]^{2}\leq\operatorname{\mathbb{E}}(1-r_{1})\operatorname{\mathbb{E}}\\{(\tilde{\bm{z}}_{1}^{\top}\bm{u})^{2}\\}$
$\displaystyle=\bm{u}^{\top}\\{\operatorname{\mathbb{E}}(1-r_{1})\operatorname{\mathbb{E}}(\tilde{\bm{z}}_{1}\tilde{\bm{z}}_{1}^{\top})\\}\bm{u},$

where the inequality holds with equality if and only if
$(1-r_{1})\propto\bm{u}^{\top}\tilde{\bm{z}}_{1}$ a.s. Since this does not
hold for any non-random $\bm{u}$, we must have

$\displaystyle\operatorname{\mathbb{E}}\\{(1-r_{1})\tilde{\bm{z}}_{1}\\}\operatorname{\mathbb{E}}\\{(1-r_{1})\tilde{\bm{z}}_{1}\\}^{\top}$
$\displaystyle\prec\operatorname{\mathbb{E}}(1-r_{1})\underbrace{\operatorname{\mathbb{E}}(\tilde{\bm{z}}_{1}\tilde{\bm{z}}_{1}^{\top})}_{=I_{K+2}}$
$\displaystyle\Rightarrow\operatorname{\mathbb{E}}\\{(1-r_{1})\tilde{\bm{z}}_{1}\\}^{\top}\operatorname{\mathbb{E}}\\{(1-r_{1})\tilde{\bm{z}}_{1}\\}$
$\displaystyle=\lVert\operatorname{\mathbb{E}}\\{(1-r_{1})\tilde{\bm{z}}_{1}\\}\operatorname{\mathbb{E}}\\{(1-r_{1})\tilde{\bm{z}}_{1}\\}^{\top}\rVert_{2}<\operatorname{\mathbb{E}}(1-r_{1})\lVert\operatorname{\mathbb{E}}(\tilde{\bm{z}}_{1}\tilde{\bm{z}}_{1}^{\top})\rVert_{2}$
$\displaystyle=\operatorname{\mathbb{E}}(1-r_{1}),$

which implies $\hat{\sigma}^{2}=m^{2}\\{c+o_{P}(1)\\}$ for some constant
$c>0$. Lastly,

$\displaystyle\hat{\sigma}^{2}=$
$\displaystyle\tilde{n}^{-1}\bm{y}_{I}^{\top}P_{[\bm{1},\bm{C}]}^{\perp}\bm{y}_{I}-\\{\tilde{n}^{-1}\bm{y}_{I}^{\top}\bm{x}-\tilde{n}^{-1}\bm{y}_{I}^{\top}[\bm{1},\bm{C}]\hat{\bm{M}}^{-1}(n^{-1}[\bm{1},\bm{C}]^{\top}\bm{x})\\}\hat{v}^{-1}$
$\displaystyle\times\\{n^{-1}\bm{y}_{I}^{\top}\bm{x}-n^{-1}\bm{y}_{I}^{\top}[\bm{1},\bm{C}]\hat{\bm{M}}^{-1}(n^{-1}[\bm{1},\bm{C}]^{\top}\bm{x})\\}$
$\displaystyle\hat{\bm{M}}=$ $\displaystyle
n^{-1}[\bm{1},\bm{C}]^{\top}[\bm{1},\bm{C}],\quad\hat{v}=n^{-1}\bm{x}^{\top}P_{[\bm{1},\bm{C}]}^{\perp}\bm{x}.$

Since $\bm{x}$ is mean 0 and independent of $\\{\bm{y},\bm{C}\\}$, it is easy
to see that

$\displaystyle
n^{-1}\bm{y}_{I}^{\top}\bm{x}-n^{-1}\bm{y}_{I}^{\top}[\bm{1},\bm{C}]\hat{\bm{M}}^{-1}(n^{-1}[\bm{1},\bm{C}]^{\top}\bm{x})=o_{P}(1),$

which completes the proof. ∎

##  S8 Theoretical guarantees for factor analysis and differential abundance

Section S8 proves theoretical statements from Sections 5.2 and 5.3 in the main
text. We prove Theorem 5.4 from Section 5.4 in Section S9.

###  S8.1 Problem statement and preliminaries

We consider the model
$\bm{y}_{g}=\bm{X}\bm{\beta}_{g}+\bm{C}\bm{\ell}_{g}+\bm{e}_{g}$, where
$\bm{e}_{g}\sim(0,\sigma_{g}^{2}I_{n})$ is sub-Gaussian with independent
entries. We also define the diagonal matrix of weights
$\bm{W}_{g}=\operatorname{diag}(w_{g1},\ldots,w_{gn})$ to be
$r_{gi}\\{\pi_{g}(y_{gi})\\}^{-1}$ for
$\pi_{g}(y_{gi})=\Psi\\{\alpha_{g}(y_{gi}-\delta_{g})\\}$. Note that
$\operatorname{\mathbb{E}}(\bm{W}_{g}\mid\bm{y}_{g},\bm{X},\bm{C})=I_{n}$. Let
$\bm{L}=[\bm{\ell}_{1}\cdots\bm{\ell}_{p}]^{\top}\in\mathbb{R}^{p\times K}$,
$\lambda=np^{-1}\operatorname{Tr}(\bm{L}^{\top}\bm{L})$, and
$\bm{B}=[\bm{\beta}_{1}\cdots\bm{\beta}_{p}]^{\top}\in\mathbb{R}^{p\times d}$.
If $\hat{w}_{gi}=w_{gi}$, the optimization problem in (4.2) is equivalent to

$\displaystyle(P_{\bm{X}}^{\perp}\hat{\bm{C}},\hat{\bm{L}},\hat{\bm{B}})\in\operatorname*{argmin}_{\begin{subarray}{c}\bm{C}_{\perp},\bm{L},\bm{B}\\\
\bm{X}^{\top}\bm{C}_{\perp}=\bm{0}\\\
\bm{C}_{\perp}^{\top}\bm{C}_{\perp}=I_{K}\end{subarray}}\frac{1}{2\lambda
p}\sum_{g=1}^{p}(\bm{y}_{g}-\bm{C}_{\perp}\bm{\ell}_{g}-\bm{X}\bm{\beta}_{g})^{\top}\bm{W}_{g}(\bm{y}_{g}-\bm{C}_{\perp}\bm{\ell}_{g}-\bm{X}\bm{\beta}_{g}).$

Define
$\bm{P}_{g}^{\perp}=\bm{W}_{g}-\bm{W}_{g}\bm{X}(\bm{X}^{\top}\bm{W}_{g}\bm{X})^{-1}\bm{X}^{\top}\bm{W}_{g}$.
Solving for $\bm{B}$ and using the fact that
$\hat{\bm{\ell}}_{g}=(\bm{C}_{\perp}^{\top}\bm{P}_{g}^{\perp}\bm{C}_{\perp})^{-1}\bm{C}_{\perp}^{\top}\bm{P}_{g}^{\perp}\bm{y}_{g}$,
the profile likelihood for $\bm{C}_{\perp}$ can be expressed as

$\displaystyle
P_{\bm{X}}^{\perp}\hat{\bm{C}}\in\operatorname*{argmax}_{\begin{subarray}{c}\bm{U}\in\mathbb{R}^{n\times
K}\\\ \bm{X}^{\top}\bm{U}=\bm{0}\\\
\bm{U}^{\top}\bm{U}=I_{K}\end{subarray}}f(\bm{U}),\quad
f(\bm{U})=\frac{1}{2\lambda
p}\sum_{g=1}^{p}\operatorname{Tr}\\{(\bm{U}^{\top}\bm{P}_{g}^{\perp}\bm{U})^{-1}\bm{U}^{\top}\bm{P}_{g}^{\perp}\bm{y}_{g}\bm{y}_{g}^{\top}\bm{P}_{g}^{\perp}\bm{U}\\}$
(S8.1)

Expanding $\bm{y}_{g}\bm{y}_{g}^{\top}$, the objective function can be
expressed as

$\displaystyle\begin{aligned} f(\bm{U})=&\frac{1}{2\lambda
p}\sum_{g=1}^{p}\operatorname{Tr}\\{(\bm{U}^{\top}\bm{P}_{g}^{\perp}\bm{U})^{-1}\bm{U}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{C}}\tilde{\bm{\ell}}_{g}\tilde{\bm{\ell}}_{g}^{\top}\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\bm{U}\\}\\\
&+\frac{1}{\lambda
p}\operatorname{Tr}\\{(\bm{U}^{\top}\bm{P}_{g}^{\perp}\bm{U})^{-1}\bm{U}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{C}}\tilde{\bm{\ell}}_{g}\bm{e}_{g}^{\top}\bm{P}_{g}^{\perp}\bm{U}\\}\\\
&+\frac{1}{2\lambda
p}\operatorname{Tr}\\{(\bm{U}^{\top}\bm{P}_{g}^{\perp}\bm{U})^{-1}\bm{U}^{\top}\bm{P}_{g}^{\perp}\bm{e}_{g}\bm{e}_{g}^{\top}\bm{P}_{g}^{\perp}\bm{U}\\}\\\
\tilde{\bm{C}}=&P_{X}^{\perp}\bm{C}(\bm{C}^{\top}P_{X}^{\perp}\bm{C})^{-1/2},\quad\tilde{\bm{\ell}}_{g}=(\bm{C}^{\top}P_{X}^{\perp}\bm{C})^{1/2}\bm{\ell}_{g}.\end{aligned}$
(S8.2)

We use this expression to prove the consistency of
$P_{\bm{X}}^{\perp}\hat{\bm{C}}$ in Section S8.3 and derive its properties and
rate of convergence in Section S8.4.

###  S8.2 Assumptions

We first re-state the assumptions on $y_{gi}$ and $\Psi$ from Section 5.1 with
a change in the scaling of the eigenvalues $\lambda_{1},\ldots,\lambda_{K}$
defined in Assumption 5.1. Note that the change is without loss of generality.

######  Assumption S8.4.

For $g\in[p]$, $i\in[n]$, and $s\in[S]$, let
$y_{gi}=\bm{\beta}_{g}^{\top}\bm{x}_{i}+\bm{\ell}_{g}^{\top}\bm{c}_{i}+e_{gi}$
and define $G_{si}\in\mathbb{R}$ to be a random variable. Then the following
hold for constant $a>0$ and $\epsilon\in(0,1/2\wedge a)$:

  1. (a)

$\bm{X}=[\bm{x}_{1}\cdots\bm{x}_{n}]^{\top}$ and
$\\{\bm{\beta}_{g}\\}_{g\in[p]}$ are non-random and satisfy
$\lvert\bm{X}_{ij}\rvert,\lVert\bm{\beta}_{g}\rVert_{2}\leq a$,
$\bm{1}_{n}\in\operatorname{Im}(\bm{X})$, and
$n^{-1}\bm{X}^{\top}\bm{X}\succeq\epsilon I_{d}$.

  2. (b)

The matrix $\bm{G}=[G_{si}]\in\mathbb{R}^{S\times n}$ is mean 0, has
independent and uniformly bounded entries, and identically distributed
columns.

  3. (c)

The eigenvalues $\lambda_{1},\ldots,\lambda_{K}>0$ of
$np^{-1}\sum_{g=1}^{p}\bm{\ell}_{g}\bm{\ell}_{g}^{\top}$ and $\bm{\ell}_{g}$
satisfy
$n^{1/2+\epsilon}\allowbreak\lesssim\allowbreak\lambda_{K}\allowbreak\leq\allowbreak\cdots\leq\allowbreak\lambda_{1}\allowbreak\lesssim
n$, $\lambda_{1}/\lambda_{K}\leq a$, and
$\left\lVert\bm{\ell}_{g}\right\rVert_{2}\leq a(\lambda_{1}/n)^{1/2}$.

  4. (d)

$\bm{C}=[\bm{c}_{1}\cdots\bm{c}_{n}]^{\top}=\bm{G}^{\top}\bm{\gamma}^{(c)}+\bm{\Delta}^{(c)}\in\mathbb{R}^{n\times
K}$ for $\bm{\gamma}^{(c)}\in\mathbb{R}^{S\times K}$ and
$\bm{\Delta}^{(c)}_{i*}\in\mathbb{R}^{K}$ such that:

    1. (i)

$\bm{\gamma}^{(c)}$ is non-random and $\lVert\bm{\gamma}^{(c)}\rVert_{1}\leq
a$, $\lVert\bm{\gamma}^{(c)}\rVert_{\infty}=o(n^{-1/4})$,
$\sum_{s=1}^{S}1\\{\bm{\gamma}^{(c)}_{s*}\neq 0\\}\leq ap^{1/2}$.

    2. (ii)

The rows of
$\bm{\Delta}^{(c)}-\operatorname{\mathbb{E}}\\{\bm{\Delta}^{(c)}\\}$ are
independent and identically distributed,
$\operatorname{\mathbb{V}}\\{\bm{\Delta}^{(c)}_{i*}\\}\succeq\epsilon I_{K}$,
and $\operatorname{\mathbb{E}}\\{\lvert\bm{\Delta}^{(c)}_{ik}\rvert^{m}\\}\leq
a_{m}$ for all $i\in[n]$, $k\in[K]$, $m\geq 1$ and constant $a_{m}>0$ that may
depend on $m$.

  5. (e)

$\bm{E}=[e_{gi}]=\\{\bm{\gamma}^{(e)}\\}^{\top}\bm{G}+\bm{\Delta}^{(e)}\in\mathbb{R}^{p\times
n}$, where $\bm{\gamma}^{(e)}\in\mathbb{R}^{S\times p}$,
$\bm{\Delta}^{(e)}\in\mathbb{R}^{p\times n}$ satisfy:

    1. (i)

$\sup_{s\in[S];g\in[p]}\lvert\bm{\gamma}^{(e)}_{sg}\rvert\leq an^{-1/4}$ and
$\sup_{g\in[p]}\sum_{s=1}^{S}1\\{\bm{\gamma}^{(e)}_{sg}\neq 0\\}\leq a$.

    2. (ii)

The columns of $\bm{\gamma}^{(e)}$ can be partitioned into disjoint sets
containing $\leq a$ metabolites, where
$\bm{\gamma}^{(e)}_{sg}\bm{\gamma}^{(e)}_{sh}=0$ if columns $g$ and $h$ lie in
different sets.

    3. (iii)

$\bm{\Delta}^{(e)}\sim MN_{p\times
n}\\{0,\operatorname{diag}(\sigma_{1}^{2},\ldots,\sigma_{p}^{2}),I_{n}\\}$,
where $\sigma_{1},\ldots,\sigma_{p}^{2}\in[\epsilon,a]$.

  6. (f)

$\bm{G}$, $\bm{\Delta}^{(c)}$, and $\bm{\Delta}^{(e)}$ are independent, and
$p,n,S$ satisfy $p\in[\epsilon n,an]$.

  7. (g)

In differential abundance applications, $\bm{X}$ can be written as
$\bm{X}=(\bm{X}_{1},\bm{X}_{2})$ for $\bm{X}_{1}\in\mathbb{R}^{n\times d_{1}}$
and $\bm{X}_{2}\in\mathbb{R}^{n\times d_{2}}$, where $\bm{X}_{1}$ are the
$d_{1}$ covariates of interests and
$p^{-1}\sum_{g=1}^{p}1\\{\bm{\beta}_{g_{j}}\neq
0\\}=o(\lambda_{1}^{1/2}n^{-1})$ for all $j\in[d_{1}]$.

  8. (h)

Assumption 5.2 from the main text holds.

Assumption (a) contains all the regularity conditions on the design matrix
$\bm{X}$ mentioned in Section 5.1. The assumption that
$\bm{1}_{n}\in\operatorname{Im}(\bm{X})$ makes the assumption that
$\operatorname{\mathbb{E}}(\bm{G})=0$ in (b) without loss of generality. Note
$\operatorname{\mathbb{E}}\\{\bm{\Delta}^{(c)}\\}$ in (d) may depend on
$\bm{X}$. The assumptions in (d) are more general than (b) from Assumption 5.1
in the main text, since the latter assumes
$\operatorname{\mathbb{E}}(\bm{c}_{i})-\\{\bm{\gamma}^{(c)}\\}^{\top}\bm{G}_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\displaystyle*$}}}}}{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\textstyle*$}}}}}{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\scriptstyle*$}}}}}{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\scriptscriptstyle*$}}}}}i}$
is continuous in $\bm{x}_{i}$, whereas the former only assumes
$\lVert\operatorname{\mathbb{E}}(\bm{c}_{i})-\\{\bm{\gamma}^{(c)}\\}^{\top}\bm{G}_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\displaystyle*$}}}}}{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\textstyle*$}}}}}{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\scriptstyle*$}}}}}{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\scriptscriptstyle*$}}}}}i}\rVert_{2}$
is bounded from above. The eigenvalues $\lambda_{1},\ldots,\lambda_{K}$ in (c)
have been scaled by a factor of $n$ to make notation in the below theoretical
statements simpler, and to be consistent with [8, 25, 24]. Assumption (g)
gives the sparsity assumption utilized in the statements of Theorem 5.2 and
Theorem 5.3 in the main text. Note this is only needed in differential
abundance applications, and is not needed to prove Theorem 5.1, Propoposition
5.2, or Theorem 5.4. We next prove two useful lemmas about $\bm{C}$ and its
relationship with $\bm{E}$ that we will use in the theoretical results that
follow.

######  Lemma S8.1.

Suppose Assumption S8.4 holds. Then
$\lim_{n\to\infty}\operatorname{Pr}(n^{-1}\bm{C}^{\top}P_{\bm{X}}^{\perp}\bm{C}\succeq\epsilon
I_{K})=1$ and
$\operatorname{\mathbb{E}}(\lvert\bm{C}_{ik}\rvert^{m})\leq\tilde{a}_{m}$ for
all $m\geq 1$ and some constant $\tilde{a}_{m}>0$ that may depend on $m$.

###### Proof.

Define $\bm{\mu}=\operatorname{\mathbb{E}}\\{\bm{\Delta}^{(c)}\\}$. Then

$\displaystyle\bm{C}^{\top}P_{\bm{X}}^{\perp}\bm{C}\succeq$
$\displaystyle\\{\bm{\Delta}^{(c)}-\bm{\mu}\\}^{\top}P_{\bm{X}}^{\perp}\\{\bm{\Delta}^{(c)}-\bm{\mu}\\}+\bm{\mu}^{\top}P_{\bm{X}}^{\perp}\\{\bm{\Delta}^{(c)}-\bm{\mu}\\}+[\bm{\mu}^{\top}P_{\bm{X}}^{\perp}\\{\bm{\Delta}^{(c)}-\bm{\mu}\\}]^{\top}$
$\displaystyle+\\{\bm{\Delta}^{(c)}\\}^{\top}P_{\bm{X}}^{\perp}\bm{G}^{\top}\\{\bm{\gamma}^{(c)}\\}+[\\{\bm{\Delta}^{(c)}\\}^{\top}P_{\bm{X}}^{\perp}\bm{G}^{\top}\\{\bm{\gamma}^{(c)}\\}]^{\top}.$

First, $\bm{G}_{i}^{\top}\\{\bm{\gamma}^{(c)}\\}\leq c$ for some constant
$c>0$ by (b) and (d) in Assumption S8.4, meaning $\lVert
P_{\bm{X}}^{\perp}\bm{G}^{\top}\\{\bm{\gamma}^{(c)}\\}\rVert_{2}^{2}\leq
nc^{2}$. This means
$\lVert\\{\bm{\Delta}^{(c)}\\}^{\top}P_{\bm{X}}^{\perp}\bm{G}^{\top}\\{\bm{\gamma}^{(c)}\\}\rVert_{2}=O_{P}(n^{1/2})$.
Since $\lVert\bm{\mu}\rVert_{2}^{2}=O(n)$ by (d) in Assumption S8.4, we also
have
$\lVert\bm{\mu}^{\top}P_{\bm{X}}^{\perp}\\{\bm{\Delta}^{(c)}-\bm{\mu}\\}\rVert_{2}=O_{P}(n^{1/2})$.
Since the rows of $\bm{\Delta}^{(c)}-\bm{\mu}$ are independent and identically
distributed,
$n^{-1}\\{\bm{\Delta}^{(c)}-\bm{\mu}\\}^{\top}P_{\bm{X}}^{\perp}\\{\bm{\Delta}^{(c)}-\bm{\mu}\\}=\\{\bm{\Delta}^{(c)}-\bm{\mu}\\}^{\top}\\{\bm{\Delta}^{(c)}-\bm{\mu}\\}\succeq\\{c+o_{P}(1)\\}I_{K}$
for some constant $c>0$, which proves
$\lim_{n\to\infty}\operatorname{Pr}(n^{-1}\bm{C}^{\top}P_{\bm{X}}^{\perp}\bm{C}\succeq\epsilon
I_{K})=1$.

Let $\lVert z\rVert_{m}=\\{\operatorname{\mathbb{E}}(\lvert
x\rvert^{m})\\}^{1/m}$ for any random variable $z$. Then for some constant
$c>0$ and $a_{m}$ as defined in Assumption S8.4,
$\lVert\bm{C}_{ik}\rVert_{m}\leq\lVert\bm{G}_{*i}^{\top}\bm{\gamma}_{*k}^{(c)}\rVert_{m}+\lVert\bm{\Delta}_{ik}^{(c)}\rVert_{m}\leq
c+a_{m}$. ∎

######  Lemma S8.2.

Suppose Assumption S8.4 holds, let $\epsilon>0$ be any constant and $m>0$ any
integer, and let $h_{i}:\mathbb{R}^{K}\to\mathbb{R}$, $i\in[n]$, be uniformly
bounded functions with uniformly bounded gradients. Then for
$\bm{L}=[\bm{\ell}_{1}^{\top}\cdots\bm{\ell}_{p}^{\top}]^{\top}$,

$\displaystyle\operatorname{\mathbb{E}}\\{(n^{-1/2}\bm{C}_{*k}^{\top}\bm{E}_{g*})^{2m}\\}\leq
c_{m},\quad k\in[K];g\in[p]$ (S8.3a)
$\displaystyle\textstyle\operatorname{\mathbb{E}}[\\{n^{-1/2}\sum_{i=1^{n}}h_{i}(\bm{C}_{i*})\bm{E}_{gi}\\}^{2m}]\leq
c_{m},\quad g\in[p]$ (S8.3b)
$\displaystyle\lVert(\lambda_{K}p)^{-1/2}\bm{C}^{\top}\bm{E}^{\top}\bm{L}\rVert_{2}=O_{P}(1)$
(S8.3c)
$\displaystyle\lVert(\lambda_{K}p)^{-1/2}\sum_{g=1}^{p}\sum_{i=1}^{n}h_{i}(\bm{C}_{i*})\bm{E}_{gi}\bm{\ell}_{g}\rVert_{2}=O_{P}(1)$
(S8.3d)

for some constant $c_{m}>0$ that may depend on $m$.

###### Proof.

Under Assumption S8.4,

$\displaystyle\bm{C}_{*k}^{\top}\bm{E}_{g*}=\\{\bm{\gamma}^{(c)}_{*k}\\}^{\top}\bm{G}\bm{G}^{\top}\bm{\gamma}^{(e)}_{*g}+\\{\bm{\Delta}^{(c)}_{*k}\\}^{\top}\bm{G}^{\top}\bm{\gamma}^{(e)}_{*g}+\bm{C}_{*k}^{\top}\bm{\Delta}^{(e)}_{*g}$

Let $\epsilon>0$ be any constant. Since $\bm{C}$ is independent of
$\bm{\Delta}^{(e)}_{*g}$,
$\operatorname{\mathbb{E}}\\{(\bm{C}_{*k}^{\top}\bm{\Delta}^{(e)}_{*g})^{2m}\\}\leq
c_{1,m}$ for some constant $c_{1,m}>0$ by Corollary S8.1. Further, since at
most finitely many entries of $\bm{\gamma}^{(e)}_{*g}\in\mathbb{R}^{S}$ are
non-zero, $\bm{G}$ is independent of $\bm{\Delta}^{(c)}$, and the rows of
$\bm{G}$ are mean 0, independent and sub-Gaussian,
$\operatorname{\mathbb{E}}[\\{n^{-1/2}\\{\bm{\Delta}_{*k}^{(c)}\\}^{\top}\bm{G}^{\top}\bm{\gamma}^{(e)}_{*g}\\}^{2m}]\leq
c_{2,m}$ for some constant $c_{2,m}>0$ by Corollary S8.1. Let
$\mathcal{I}_{g}=\\{s\in[S]:\bm{\gamma}^{(e)}_{sg}\neq 0\\}$ and
$\mathcal{C}=\\{s\in[S]:\bm{\gamma}^{(c)}_{sk}\neq 0\\}$. Then

$\displaystyle
n^{-1/2}\\{\bm{\gamma}^{(c)}_{*k}\\}^{\top}\bm{G}\bm{G}^{\top}\bm{\gamma}^{(e)}_{*g}=\sum_{s\in\mathcal{I}_{g}}\tilde{\bm{\gamma}}^{(c)}_{sk}\tilde{\bm{\gamma}}^{(e)}_{sg}(n^{-1}\bm{G}_{s*}^{\top}\bm{G}_{s*})+n^{-1/2}\sum_{r\in\mathcal{I}_{g}^{c}\cap\mathcal{C}}\bm{\gamma}^{(c)}_{rk}\bm{G}_{r*}^{\top}\sum_{s\in\mathcal{I}_{g}}\bm{\gamma}_{sg}^{(e)}\bm{G}_{s*}$
$\displaystyle\tilde{\bm{\gamma}}^{(c)}=n^{-1/4}\bm{\gamma}^{(c)},\quad\tilde{\bm{\gamma}}^{(e)}=n^{-1/4}\bm{\gamma}^{(e)}.$

First, since
$\operatorname{\mathbb{E}}\\{(n^{-1}\bm{G}_{s*}^{\top}\bm{G}_{s*})^{2m}\\}\leq
c_{3,m}$ for all $s\in[S]$ and some constant $c_{3,m}>0$ that may depend on
$m$,

$\displaystyle\operatorname{\mathbb{E}}[\textstyle\\{\sum_{s\in\mathcal{I}_{g}}\tilde{\bm{\gamma}}^{(c)}_{sk}\tilde{\bm{\gamma}}^{(e)}_{sg}(n^{-1}\bm{G}_{s*}^{\top}\bm{G}_{s*})\\}^{2m}]\leq
ac_{3,m}$

for some constant $a>0$ because
$\lvert\tilde{\bm{\gamma}}^{(c)}_{sk}\tilde{\bm{\gamma}}^{(e)}_{sg}\rvert$ and
$\lvert\mathcal{I}_{g}\rvert$ are uniformly bounded from above. Let
$\bm{Z}_{g}=\sum_{s\in\mathcal{I}_{g}}\bm{\gamma}^{(e)}_{sg}\bm{G}_{s*}\in\mathbb{R}^{n}$.
Then $\bm{Z}_{g}$ is mean 0, has independent entries that are bounded above by
$an^{-1/4}$ for some constant $a>0$, and is independent of
$\bm{W}_{g}=\sum_{r\in\mathcal{I}_{g}^{c}\cap\mathcal{C}}\bm{\gamma}^{(c)}_{rk}\bm{G}_{r*}\in\mathbb{R}^{n}$.
Here, $\bm{W}_{g}$ is also mean 0, has independent entries. Further, since
$\\{\bm{G}_{r*}\\}_{r\in\mathcal{I}_{g}^{c}\cap\mathcal{C}}$ are independent
and sub-Gaussian with uniformly bounded sub-Gaussian norm and
$\sum_{r\in[S]}\\{\bm{\gamma}^{(c)}_{rk}\\}^{2}\leq\sum_{r\in[S]}\lvert\bm{\gamma}^{(c)}_{rk}\rvert\leq
a$ for some constant $a>0$ by Assumption S8.4, $\bm{W}_{g}$ is also sub-
Gaussian with uniformly bounded sub-Gaussian norm over $g\in[p]$. All this
implies

$\displaystyle\operatorname{\mathbb{E}}[\textstyle\\{n^{-1/2}\sum_{r\in\mathcal{I}_{g}^{c}\cap\mathcal{C}}\bm{\gamma}^{(c)}_{rk}\bm{G}_{r*}^{\top}\sum_{s\in\mathcal{I}_{g}}\bm{\gamma}_{sg}^{(e)}\bm{G}_{s*}\\}^{2m}]=\operatorname{\mathbb{E}}\\{(n^{-1/2}\bm{W}_{g}^{\top}\bm{Z}_{g})^{2m}\\}\leq
c_{4,m},\quad g\in[p]$

for some constant $c_{4,m}>0$ that may depend on $m$. This completes the proof
of (S8.3a).

Since $h$ is bounded from above, the above work implies we need only show that
$\operatorname{\mathbb{E}}[\\{n^{-1/2}\allowbreak\sum_{i=1}^{n}\allowbreak
h_{i}(\bm{C}_{i*})\bm{G}_{*i}^{\top}\bm{\gamma}_{*g}^{(e)}\\}^{2m}]\allowbreak\leq\allowbreak
c_{m}$ to prove (S8.3b). For simplicity, we assume that
$\mathcal{I}_{g}=\\{s_{g}\\}$, and note that the extension to general
$\mathcal{I}_{g}$ under the Assumption S8.4 is trivial. Then for
$\bm{R}^{(g)}_{i*}=\bm{C}_{i*}-\bm{\gamma}_{s_{g}*}^{(c)}\bm{G}_{s_{g}i}$,

$\displaystyle
n^{-1/2}\sum_{i=1}^{n}h_{i}(\bm{C}_{i*})\bm{G}_{*i}^{\top}\bm{\gamma}_{*g}^{(e)}=$
$\displaystyle
n^{-1/2}\bm{\gamma}_{s_{g}g}^{(e)}\sum_{i=1}^{n}h_{i}\\{\bm{\gamma}_{s_{g}*}^{(c)}\bm{G}_{s_{g}i}+\bm{R}^{(g)}_{i*}\\}\bm{G}_{s_{g}i}$
$\displaystyle=$ $\displaystyle
n^{-1/2}\bm{\gamma}_{s_{g}g}^{(e)}\sum_{i=1}^{n}h_{i}\\{\bm{R}^{(g)}_{i*}\\}\bm{G}_{s_{g}i}+O(1),$

where the second equality follows because $\\{h_{i}\\}_{i\in[n]}$ have
uniformly bounded gradients, $\bm{G}_{s_{g}i}$ is bounded, and
$\lvert\bm{\gamma}_{s_{g}g}^{(e)}\rvert\lVert\bm{\gamma}_{s_{g}*}^{(c)}\rVert_{2}=o(n^{-1/2})$.
An application of Lemma S8.4 then proves the result, which proves (S8.3b).

For (S8.3c),

$\displaystyle(\lambda_{K}p)^{-1/2}\bm{C}^{\top}\bm{E}^{\top}\bm{L}=$
$\displaystyle(\lambda_{K}p)^{-1/2}\bm{C}^{\top}\\{\bm{\Delta}^{(e)}\\}^{\top}\bm{L}+(\lambda_{K}p)^{-1/2}\\{\bm{\Delta}^{(c)}\\}^{\top}\bm{G}^{\top}\bm{\gamma}^{(e)}\bm{L}$
$\displaystyle+(\lambda_{K}p)^{-1/2}\\{\bm{\gamma}^{(c)}\\}^{\top}\bm{G}\bm{G}^{\top}\bm{\gamma}^{(e)}\bm{L}.$

It is straightforward to show that

$\displaystyle\lVert(\lambda_{K}p)^{-1/2}\bm{C}^{\top}\\{\bm{\Delta}^{(e)}\\}^{\top}\bm{L}\rVert_{2},\,\lVert(\lambda_{K}p)^{-1/2}\\{\bm{\Delta}^{(c)}\\}^{\top}\bm{G}^{\top}\bm{\gamma}^{(e)}\bm{L}\rVert_{2}=O_{P}(1).$

For the remaining term, let
$\tilde{\bm{L}}=n^{3/4}\lambda_{K}^{-1/2}\bm{\gamma}^{(e)}\bm{L}\in\mathbb{R}^{S\times
K}$. Then at most $O(p)$ rows of $\tilde{\bm{L}}$ are non-zero and
$\lVert\tilde{\bm{L}}_{s*}\rVert_{2}\leq c$ for some constant $c>0$ by
Assumption S8.4. Then for
$\mathcal{S}=\\{s\in[S]:\bm{\gamma}_{s*},\tilde{\bm{L}}_{s*}\neq 0\\}$ and
$\mathcal{R}=\\{(r,s)\in[S]\times[S]:r\neq
s,\bm{\gamma}_{s*}\odot\tilde{\bm{L}}_{s*}\neq 0\\}$ (where $\odot$ is the
Hadamard product),

$\displaystyle(\lambda_{K}p)^{-1/2}\\{\bm{\gamma}^{(c)}\\}^{\top}\bm{G}\bm{G}^{\top}\bm{\gamma}^{(e)}\bm{L}=$
$\displaystyle
p^{-1/2}\sum_{s\in\mathcal{S}}\tilde{\bm{\gamma}}^{(c)}_{s*}\tilde{\bm{L}}_{s*}^{\top}(n^{-1}\bm{G}_{s*}^{\top}\bm{G}_{s*})$
$\displaystyle+n^{-1}\sum_{(r,s)\in\mathcal{R}}\bm{\gamma}^{(c)}_{s*}\tilde{\bm{L}}_{s*}^{\top}(p^{-1/2}\bm{G}_{r*}^{\top}\bm{G}_{s*}).$

Since
$\lVert\tilde{\bm{\gamma}}^{(c)}_{s*}\rVert_{2},\lVert\tilde{\bm{L}}_{s*}\rVert_{2}\leq
c$ and $\operatorname{\mathbb{E}}(n^{-1}\bm{G}_{s*}^{\top}\bm{G}_{s*})\leq c$
for some constant $c>0$ and all $s\in[S]$,

$\displaystyle\operatorname{\mathbb{E}}\\{\lVert
p^{-1/2}\sum_{s\in\mathcal{S}}\tilde{\bm{\gamma}}^{(c)}_{s*}\tilde{\bm{L}}_{s*}^{\top}(n^{-1}\bm{G}_{s*}^{\top}\bm{G}_{s*})\rVert_{2}\\}\leq
c^{3}p^{-1/2}\lvert\mathcal{S}\rvert=O(1)$

since $\lvert\mathcal{S}\rvert=O(p^{1/2})$ by Assumption S8.4. Lastly, since
$\operatorname{\mathbb{E}}(\bm{G}_{r*}^{\top}\bm{G}_{s*}\bm{G}_{r^{\prime}*}^{\top}\bm{G}_{s^{\prime}*})=0$
for all $(r,s)\neq(r^{\prime},s^{\prime})$,

$\displaystyle\operatorname{\mathbb{E}}[\\{n^{-1}\sum_{(r,s)\in\mathcal{R}}\bm{\gamma}^{(c)}_{sk_{1}}\tilde{\bm{L}}_{sk_{2}}^{\top}(p^{-1/2}\bm{G}_{r*}^{\top}\bm{G}_{s*})\\}^{2}]=O(n^{-2}\lvert\mathcal{R}\rvert)=O(n^{-2}p^{3/2})=O(1),$

which completes the proof of (S8.3c).

Lastly, to prove (S8.3d), we again assume for simplicity that
$\mathcal{I}_{g}=\\{s_{g}\\}$, but note that extending it to general
$\mathcal{I}_{g}$ is trivial under Assumption S8.4. Since $h_{i}$ is bounded,
the proof of (S8.3c) implies we need only consider the behavior of
$\sum_{g=1}^{p}\bm{\gamma}_{s_{g}g}^{(e)}\sum_{i=1}^{n}h_{i}(\bm{C}_{i*})\bm{G}_{s_{g}i}\bm{\ell}_{g}$.
Fix a $k\in[K]$ and let $a_{gk}=n^{1/2}\lambda_{K}^{-1/2}\bm{\ell}_{g_{k}}$,
where $\lvert a_{gi}\rvert$ is uniformly bounded from above by Assumption
S8.4. Let $\mathcal{C}=\\{g\in[p]:\bm{\gamma}_{s_{g}*}^{(c)}\neq 0\\}$. Then
$\lvert\mathcal{C}\rvert=O(p^{1/2})$ by Assumption S8.4 and for any $k\in[K]$,

$\displaystyle(\lambda_{K}p)^{-1/2}\sum_{g=1}^{p}\bm{\gamma}_{s_{g}g}^{(e)}\sum_{i=1}^{n}h_{i}(\bm{C}_{i*})\bm{G}_{s_{g}i}\bm{\ell}_{g_{k}}=(\lambda_{K}p)^{-1/2}\sum_{g\in\mathcal{C}}\bm{\gamma}_{s_{g}g}^{(e)}\sum_{i=1}^{n}h_{i}(\bm{C}_{i*})\bm{G}_{s_{g}i}\bm{\ell}_{g_{k}}$
$\displaystyle+(\lambda_{K}p)^{-1/2}\sum_{g\in\mathcal{C}^{c}}\bm{\gamma}_{s_{g}g}^{(e)}\sum_{i=1}^{n}h_{i}(\bm{C}_{i*})\bm{G}_{s_{g}i}\bm{\ell}_{g_{k}}.$

Since the $\bm{G}_{s_{g}}$ is independent of $\bm{C}$ for
$g\in\mathcal{C}^{c}$, the second term after the equality in the above
expression is $O_{P}(1)$ because $h_{i}$ is uniformly bounded from above. For
the first term, fix a $g\in\mathcal{C}$ and let
$\bm{C}_{i*}=\bm{\gamma}_{s_{g}*}^{(e)}\bm{G}_{s_{g}i}+\bm{R}_{i*}^{(g)}$,
where $\bm{R}_{i*}^{(g)}$ is independent of $\bm{G}_{s_{g}i}$. Then for
$\tilde{\ell}_{gk}=n^{1/2}\lambda_{K}^{-1/2}\bm{\ell}_{gk}$ (which is
uniformly bounded by Assumption S8.4),

$\displaystyle(\lambda_{K}p)^{-1/2}\sum_{g\in\mathcal{C}}\bm{\gamma}_{s_{g}g}^{(e)}\sum_{i=1}^{n}h_{i}(\bm{C}_{i*})\bm{G}_{s_{g}i}\bm{\ell}_{g_{k}}=p^{-1/2}\sum_{g\in\mathcal{C}}\tilde{\ell}_{gk}\bm{\gamma}_{s_{g}g}^{(e)}[n^{-1/2}\sum_{i=1}^{n}h_{i}\\{\bm{R}^{(g)}_{i*}\\}\bm{G}_{s_{g}i}]$
$\displaystyle+p^{-1/2}\sum_{g\in\mathcal{C}}[n^{1/2}\tilde{\ell}_{gk}\bm{\gamma}_{s_{g}g}^{(e)}\\{\bm{\gamma}_{s_{g}*}^{(c)}\\}^{\top}][n^{-1}\sum_{i=1}^{n}\bm{G}_{s_{g}i}^{2}\smallint_{0}^{1}\nabla
h_{i}\\{\bm{R}^{(g)}_{i*}+t\bm{\gamma}_{s_{g}*}^{(c)}\bm{G}_{s_{g}i}\\}dt],$

Since $h_{i}$ and $\nabla h_{i}$ are bounded,

$\displaystyle\operatorname{\mathbb{E}}(\textstyle[n^{-1/2}\sum_{i=1}^{n}h_{i}\\{\bm{R}^{(g)}_{i*}\\}\bm{G}_{s_{g}i}]^{2})\leq
c$ $\displaystyle\operatorname{\mathbb{E}}[\lVert
n^{-1}\sum_{i=1}^{n}\bm{G}_{s_{g}i}^{2}\smallint_{0}^{1}\nabla
h_{i}\\{\bm{R}^{(g)}_{i*}+t\bm{\gamma}_{s_{g}*}^{(c)}\bm{G}_{s_{g}i}\\}dt\rVert_{2}]\leq
c$

for some constant $c>0$. Since $\lvert\mathcal{C}\rvert=O(p^{1/2})$, this
completes the proof. ∎

######  Remark S8.7.

The proof of Lemma S8.2 can easily be extended to show that (S8.3) holds when
we replace $\bm{C}$ with $P_{\bm{X}}^{\perp}\bm{C}$. This will be useful in
Lemma S8.17 below.

###  S8.3 Consistency of $P_{\bm{X}}^{\perp}\hat{\bm{C}}$

Here we use $f(\bm{U})$ from (S8.2) to prove the consistency of
$P_{\bm{X}}^{\perp}\hat{\bm{C}}$. The main results are Lemmas S8.8, S8.12,
S8.13 and Corollary S8.2. For ease of notation, we re-define
$P_{\bm{X}}^{\perp}\hat{\bm{C}}$ to be $\hat{\bm{C}}$ in Sections S8.3 and
S8.4.

######  Lemma S8.3.

Let $\bm{U}\in\mathbb{R}^{n\times K}$ be a matrix with orthonormal columns
that satisfies $\bm{U}^{\top}\bm{X}=\bm{0}$, and define $\delta=\left\lVert
P_{\bm{U}}-P_{\tilde{\bm{C}}}\right\rVert_{2}$. Then there exist
$\bm{v}_{u}\in\mathbb{R}^{K\times K}$, $\bm{z}_{u}\in\mathbb{R}^{(n-K-d)\times
K}$, and some universal constant $c>1$ such that

$\displaystyle\bm{v}_{u}^{\top}\bm{v}_{u}+\bm{z}_{u}^{\top}\bm{z}_{u}=I_{K},\quad\bm{U}=\tilde{\bm{C}}\bm{v}_{u}+\bm{Q}\bm{z}_{u},\quad\left\lVert\bm{v}_{u}-\bm{v}\right\rVert_{2}\in[c^{-1}\delta^{2},c\delta^{2}],\quad\left\lVert\bm{z}_{u}\right\rVert_{2}\in[c^{-1}\delta,c\delta],$

where $\bm{v}=\bm{A}_{u}\bm{B}_{u}^{\top}\in\mathbb{R}^{K\times K}$ for
$\bm{A}_{u}\in\mathbb{R}^{K\times K}$ and $\bm{B}_{u}\in\mathbb{R}^{K\times
K}$ the left and right singular vectors of $\bm{v}_{u}$.

###### Proof.

We can express $P_{\bm{U}}=\bm{U}\bm{U}^{\top}$ and
$P_{\tilde{\bm{C}}}=\tilde{\bm{C}}\tilde{\bm{C}}^{\top}$, and can always
express $\bm{U}=\tilde{\bm{C}}\bm{v}_{u}+\bm{Q}\bm{z}_{u}$ where
$\bm{v}_{u}^{\top}\bm{v}_{u}+\bm{z}_{u}^{\top}\bm{z}_{u}=I_{K}$ by the
Fredholm alternative. Then

$\displaystyle\left\lVert
P_{\bm{U}}-P_{\tilde{\bm{C}}}\right\rVert_{2}^{2}\leq\left\lVert
P_{\bm{U}}-P_{\tilde{\bm{C}}}\right\rVert_{F}^{2}\leq K\left\lVert
P_{\bm{U}}-P_{\tilde{\bm{C}}}\right\rVert_{2}^{2}$ $\displaystyle\left\lVert
P_{\bm{U}}-P_{\tilde{\bm{C}}}\right\rVert_{F}^{2}=2\operatorname{Tr}(I_{K}-\bm{v}_{u}^{\top}\bm{v}_{u})=2\operatorname{Tr}(\bm{z}_{u}^{\top}\bm{z}_{u})=2\left\lVert\bm{z}_{u}\right\rVert_{F}^{2}$
$\displaystyle
K^{-1}\left\lVert\bm{z}_{u}\right\rVert_{F}^{2}\leq\left\lVert\bm{z}_{u}\right\rVert_{2}^{2}\leq\left\lVert\bm{z}_{u}\right\rVert_{F}^{2}.$

The first two lines imply
$\left\lVert\bm{z}_{u}\right\rVert_{F}^{2}\in[\delta^{2}/2,K\delta^{2}/2]$,
which taken with the third line, implies
$\left\lVert\bm{z}_{u}\right\rVert_{2}^{2}\in[\delta^{2}/(2K),K\delta^{2}/2]$.
Note that since $\bm{v}_{u}^{\top}\bm{v}_{u}\preceq I_{K}$, the singular
values of $\bm{v}_{u}$ satisfy $0\leq\sigma_{K}\leq\cdots\leq\sigma_{1}\leq 1$
and

$\displaystyle 1-\sigma_{K}^{2}=\left\lVert
I_{K}-\bm{v}_{u}^{\top}\bm{v}_{u}\right\rVert_{2}=\left\lVert\bm{z}_{u}\right\rVert_{2}^{2}\in[\delta^{2}/(2K),K\delta^{2}/2]$
$\displaystyle\Rightarrow
1-\sigma_{1},\ldots,1-\sigma_{K}\in[\delta^{2}/(4K),K\delta^{2}/2].$

If
$\bm{v}_{u}=\bm{A}\operatorname{diag}(\sigma_{1},\ldots,\sigma_{K})\bm{B}^{\top}$
is the singular value decomposition of $\bm{v}_{u}$, this shows that
$\left\lVert\bm{A}\bm{B}^{\top}-\bm{v}_{u}\right\rVert_{2}\in[\delta^{2}/(4K),K\delta^{2}/2]$,
which completes the proof. ∎

######  Lemma S8.4.

Suppose the random variables $z_{1},\ldots,z_{n}$ satisfy the following for
some integer $m\geq 1$ and constant $c>0$:

  1. 1.

$\operatorname{\mathbb{E}}(z_{i}^{2m})<c$

  2. (ii)

There exists a $\sigma$-algebra $\mathcal{F}$ such that
$\operatorname{\mathbb{E}}(z_{i}\mid\mathcal{F})=0$ for all $i\in[n]$ and
$z_{1},\ldots,z_{n}$ are independent conditional on $\mathcal{F}$.

Then $\operatorname{\mathbb{E}}\\{(\sum_{i=1}^{n}z_{i})^{2m}\\}\leq
cc_{m}n^{m}$, where $c_{m}$ is a constant that only depends on $m$.

###### Proof.

$\displaystyle\operatorname{\mathbb{E}}\left\\{\left(\sum_{i=1}^{n}z_{i}\right)^{2m}\right\\}=$
$\displaystyle\sum_{i_{1},\ldots,i_{2m}\in[n]}\operatorname{\mathbb{E}}(z_{i_{1}}\cdots
z_{i_{2m}})=\sum_{\begin{subarray}{c}i_{1},\ldots,i_{2m}\in[n]:\\\
\text{$\exists j\in[2m]$ such that}\\\
i_{j}\notin\\{i_{s}\\}_{s\in[2m]\setminus\\{j\\}}\end{subarray}}\operatorname{\mathbb{E}}(z_{i_{1}}\cdots
z_{i_{2m}})$
$\displaystyle+\sum_{\begin{subarray}{c}i_{1},\ldots,i_{2m}\in[n]:\\\ \text{
for all $j\in[2m]$, there exists }\\\ \text{$j^{\prime}\neq j$ such that
$i_{j}=i_{j^{\prime}}$}\end{subarray}}\operatorname{\mathbb{E}}(z_{i_{1}}\cdots
z_{i_{2m}}),$

where

$\displaystyle\sum_{\begin{subarray}{c}i_{1},\ldots,i_{2m}\in[n]:\\\
\text{$\exists j\in[2m]$ such that}\\\
i_{j}\notin\\{i_{s}\\}_{s\in[2m]\setminus\\{j\\}}\end{subarray}}\operatorname{\mathbb{E}}(z_{i_{1}}\cdots
z_{i_{2m}}\mid\mathcal{F})=\sum_{\begin{subarray}{c}i_{1},\ldots,i_{2m}\in[n]:\\\
\text{$\exists j\in[2m]$ such that}\\\
i_{j}\notin\\{i_{s}\\}_{s\in[2m]\setminus\\{j\\}}\end{subarray}}\underbrace{\operatorname{\mathbb{E}}(z_{i_{j}}\mid\mathcal{F})}_{=0}\operatorname{\mathbb{E}}\left(\prod_{s\in[2m]\setminus\\{j\\}}z_{i_{s}}\mid\mathcal{F}\right)=0$

and

$\displaystyle\sum_{\begin{subarray}{c}i_{1},\ldots,i_{2m}\in[n]:\\\ \text{
for all $j\in[2m]$, there exists }\\\ \text{$j^{\prime}\neq j$ such that
$i_{j}=i_{j^{\prime}}$}\end{subarray}}\operatorname{\mathbb{E}}(z_{i_{1}}\cdots
z_{i_{2m}})$ $\displaystyle\leq$ $\displaystyle
c\lvert\\{i_{1},\ldots,i_{2m}\in[n]:\text{ for all $j\in[2m]$, there exists
$j^{\prime}\neq j$ such that $i_{j}=i_{j^{\prime}}$}\\}\rvert$
$\displaystyle\leq$ $\displaystyle cc_{m}n^{m}$

for some constant $c_{m}>0$ that only depends on $m$. ∎

######  Corollary S8.1.

Let $r>0$ be an integer, and define the $r$ possibly dependent sets of random
variables $\\{z_{j1},\ldots,z_{jn}\\}$ to be such that $z_{j1},\ldots,z_{jn}$
satisfy the conditions of Lemma S8.4 for each $j\in[r]$. Then for
$S_{j}=n^{-1}\sum_{i=1}^{n}z_{ji}$,
$\tilde{S}_{j}=n^{-1/2}\sum_{i=1}^{n}z_{ji}$, and any $t>0$,
$\operatorname{Pr}(\max_{j\in[r]}\lvert S_{j}\rvert\geq t)\leq
cc_{m}r/(nt^{2})^{m}$ and
$\operatorname{Pr}(\max_{j\in[r]}\lvert\tilde{S}_{j}\rvert\geq t)\leq
cc_{m}r/t^{2m}$ for $c,c_{m}$ defined in the statement of Lemma S8.4.

###### Proof.

Since $\max_{j\in[r]}\lvert
S_{j}\rvert\leq\left(\sum_{j=1}^{r}S_{j}^{2m}\right)^{1/(2m)}$ and
$\max_{j\in[r]}\lvert\tilde{S}_{j}\rvert\leq\left(\sum_{j=1}^{r}\tilde{S}_{j}^{2m}\right)^{1/(2m)}$,
this follows immediately from Lemma S8.4. ∎

######  Remark S8.8.

If $r=n^{a}$ for some $a\in(0,m)$, then Corollary S8.1 implies
$\max_{j\in[r]}\lvert S_{j}\rvert=O_{P}(cn^{-\delta})$ for
$\delta=(1-a/m)/2\in(0,1/2)$.

######  Lemma S8.5.

Let $c>1$ be a constant, and assume $\bm{e}_{g}\sim(\bm{0},\bm{V}_{g})$,
$g\in[p]$, are independent sub-Gaussian random vectors with sub-Gaussian norm
$\left\lVert\bm{e}_{g}\right\rVert_{\Psi_{2}}\leq c$. Then if $p\geq c^{-1}n$
and for $\bm{E}=(\bm{e}_{1}\cdots\bm{e}_{p})$, $\lVert
p^{-1/2}\bm{E}\rVert_{2}=O_{P}(1)$ as $n,p\to\infty$.

###### Proof.

The proof is a simple extension of the proof of Theorem 5.39 in [26], and has
been omitted. ∎

######  Lemma S8.6.

Let $\bm{M}\in\mathbb{R}^{p\times n}$ such that $\bm{M}_{gi}=w_{gi}-1$, and
suppose Assumption S8.4 hold. Then for any constant $\epsilon>0$, $\lVert
p^{-1/2}\bm{M}\rVert_{2}=O_{P}(n^{\epsilon})$ as $n,p\to\infty$.

###### Proof.

Conditional of $\bm{e}_{1},\ldots,\bm{e}_{p}$ and $\bm{C}$, the entries of
$\bm{M}$ are mean 0, independent, and have finite fourth moments, where for
any integer $m>0$ and some constant $c_{m}>0$ that only depends on $m$,

$\displaystyle\operatorname{\mathbb{E}}(w_{gi}^{m}\mid\bm{C}_{i\mathchoice{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\displaystyle*$}}}}}{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\textstyle*$}}}}}{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\scriptstyle*$}}}}}{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\scriptscriptstyle*$}}}}}},\bm{e}_{gi})=$
$\displaystyle\operatorname{\mathbb{E}}([1/\Psi\\{\alpha_{g}(y_{gi}-\delta_{g})\\}]^{(m-1)}\mid\bm{C}_{i\mathchoice{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\displaystyle*$}}}}}{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\textstyle*$}}}}}{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\scriptstyle*$}}}}}{\mathbin{\vbox{\hbox{\scalebox{1.0}{$\scriptscriptstyle*$}}}}}},\bm{e}_{gi})\leq
c_{m}$
$\displaystyle+c_{m}\\{\lvert\bm{e}_{gi}\rvert^{a(m-1)}+\sum_{k=1}^{K}\lvert\bm{C}_{ik}\rvert^{a(m-1)}\\}$

for some constant $a>0$ by (h) in Assumption S8.4. The result then follows by
[27] and the fact that, for any $\epsilon>0$,
$\max_{g\in[p],i\in[n]}\lvert\bm{e}_{gi}\rvert,\max_{i\in[n],k\in[K]}\lvert\bm{C}_{ik}\rvert=O_{P}(n^{\epsilon})$.
∎

######  Lemma S8.7.

Let $\bm{M}\in\mathbb{R}^{p\times n}$ such that
$\bm{M}_{gi}=w_{gi}\bm{e}_{g_{i}}$, and suppose Assumption S8.4 holds. Then
for any fixed constant $\epsilon\in(0,1/2)$, $\left\lVert
p^{-1/2}\bm{M}\right\rVert_{2}=O_{P}(n^{\epsilon})$ as $n,p\to\infty$.

###### Proof.

We can express $\bm{M}$ as $\bm{M}=\bm{M}^{(1)}+\bm{M}^{(2)}$ for
$\bm{M}^{(1)}_{gi}=\bm{e}_{gi}$ and
$\bm{M}^{(2)}_{gi}=(w_{gi}-1)\bm{e}_{g_{i}}$. By Lemma S8.5,
$\lVert\bm{M}^{(1)}\rVert_{2}=O_{P}(1)$, and a simple extension of the proof
of Lemma S8.6 can be used to show $\lVert
p^{-1/2}\bm{M}^{(2)}\rVert_{2}=O_{P}(n^{\epsilon})$. ∎

######  Lemma S8.8.

Let $\bm{U}\in\mathbb{R}^{n\times K}$ be a matrix with orthonormal columns and
define

$\displaystyle f_{1}(\bm{U})=(\lambda
p)^{-1}\sum_{g=1}^{p}\operatorname{Tr}\\{(\bm{U}^{\top}\bm{P}_{g}^{\perp}\bm{U})^{-1}\bm{U}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{C}}\tilde{\bm{\ell}}_{g}\tilde{\bm{\ell}}_{g}^{\top}\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\bm{U}\\}.$

Let
$\delta_{U}=\lVert\tilde{\bm{C}}\tilde{\bm{C}}^{\top}-\bm{U}\bm{U}^{\top}\rVert_{2}$,
$\eta\in(0,1/2)$ be an arbitrarily small constant, and suppose Assumption S8.4
holds. Then there exists a constant $c>1$ that does not depend on $n$ or $p$
such that for all $\bm{U}$ with $\delta_{U}\in(0,c^{-1})$ and any
$\epsilon_{1},\epsilon_{2}>0$, $f_{1}(\tilde{\bm{C}})-f_{1}(\bm{U})\geq
c^{-1}\delta_{U}^{2}\\{1-c\delta_{U}(1+\epsilon_{2}n^{-1/2+\eta})\\}$ with
probability at least $1-\epsilon_{1}$ for all $n,p$ sufficiently large.

###### Proof.

For notational simplicity, we set $\delta=\delta_{U}$. Let
$\bm{U}=\bm{C}\bm{v}_{u}+\bm{Q}\bm{z}_{u}$ for $\bm{Q}$ as defined in Lemma
S8.3, where $\left\lVert\bm{z}_{u}\right\rVert_{2}\in[c^{-1}\delta,c\delta]$
and $\left\lVert\bm{v}_{u}-\bm{v}\right\rVert_{2}\leq c\delta^{2}$ for some
constant $c>1$ and $K\times K$ unitary matrix $\bm{v}$. We let
$\tilde{\bm{z}}_{u}=\bm{Q}\bm{z}_{u}$ for the remainder of the proof, and
without loos of generality, assume $n^{-1}\bm{X}^{\top}\bm{X}=I_{d}$. Provided
$\bm{U}^{\top}\bm{P}_{g}^{\perp}\bm{U}$ is invertible, define

$\displaystyle f_{1g}(\bm{U})=$
$\displaystyle\operatorname{Tr}\\{(\bm{U}^{\top}\bm{P}_{g}^{\perp}\bm{U})^{-1}\bm{U}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{C}}\tilde{\bm{\ell}}_{g}\tilde{\bm{\ell}}_{g}^{\top}\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\bm{U}\\}$
$\displaystyle=$
$\displaystyle\operatorname{Tr}\\{(\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{C}})^{-1/2}\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\bm{U}(\bm{U}^{\top}\bm{P}_{g}^{\perp}\bm{U})^{-1}\bm{U}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{C}}(\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{C}})^{-1/2}\times$
(S8.4)
$\displaystyle\times(\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{C}})^{1/2}\tilde{\bm{\ell}}_{g}\tilde{\bm{\ell}}_{g}^{\top}(\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{C}})^{1/2}\\}\leq\operatorname{Tr}\\{(\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{C}})^{1/2}\tilde{\bm{\ell}}_{g}\tilde{\bm{\ell}}_{g}^{\top}(\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{C}})^{1/2}\\}=f_{1g}(\bm{C}),$

where the inequality follows because the symmetric and positive semi-definite
matrix in the second line has eigenvalues $\leq 1$. We first see that

$\displaystyle\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\bm{U}=\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{C}}\bm{v}_{u}+\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{z}}_{u}.$

Define the $K\times K$ matrix
$\bm{A}_{g}=\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{C}}$. Then the
expression inside the $\operatorname{Tr}$ operator in (S8.4) can be written as

$\displaystyle(\bm{A}_{g}^{1/2}\bm{v}_{u}+\bm{A}_{g}^{-1/2}\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{z}}_{u})(\bm{v}_{u}^{\top}\bm{A}_{g}\bm{v}_{u}+\bm{z}_{u}^{\top}\bm{Q}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{z}}_{u}+\bm{z}_{u}^{\top}\bm{Q}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{C}}\bm{v}_{u}+\bm{v}_{u}^{\top}\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{z}}_{u})^{-1}$
$\displaystyle\times(\bm{A}_{g}^{1/2}\bm{v}_{u}+\bm{A}_{g}^{-1/2}\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{z}}_{u})^{\top}=\bm{B}_{g}\\{\bm{B}_{g}^{\top}\bm{B}_{g}+\bm{z}_{u}^{\top}\bm{Q}^{\top}(\bm{P}_{g}^{\perp}-\bm{P}_{g}^{\perp}\tilde{\bm{C}}\bm{A}_{g}^{-1}\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp})\tilde{\bm{z}}_{u}\\}^{-1}\bm{B}_{g}^{\top}$
$\displaystyle=$ $\displaystyle(I_{K}+\bm{D}_{g})^{-1}$ (S8.5)
$\displaystyle\bm{B}_{g}=\bm{A}_{g}^{1/2}\bm{v}_{u}+\bm{A}_{g}^{-1/2}\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{z}}_{u}$
$\displaystyle\bm{D}_{g}=\bm{B}_{g}^{-\top}\bm{z}_{u}^{\top}\bm{Q}^{\top}(\bm{P}_{g}^{\perp}-\bm{P}_{g}^{\perp}\tilde{\bm{C}}\bm{A}_{g}^{-1}\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp})\tilde{\bm{z}}_{u}\bm{B}_{g}^{-1}.$
(S8.6)

We first prove two lemmas that we will use throughout the proof.

######  Lemma S8.9.

Suppose Assumption S8.4 holds and let
$\tilde{\bm{B}}_{g}=\bm{A}_{g}^{-1/2}\bm{B}_{g}$. Then for all
$\epsilon\in(0,1/2)$ and some constant $c>0$ that does not depend on $n$, $p$,
or $\delta$,

$\displaystyle\max_{g\in[p]}\lVert\bm{A}_{g}-I_{K}\rVert_{2}=O_{P}(n^{-1/2+\epsilon})$
(S8.7)
$\displaystyle\max_{g\in[p]}\lVert\bm{B}_{g}^{\top}\bm{B}_{g}-I_{K}\rVert_{2},\,\max_{g\in[p]}\lVert\tilde{\bm{B}}_{g}^{\top}\tilde{\bm{B}}_{g}-I_{K}\rVert_{2}\leq
c(\delta+\delta^{2})\\{1+O_{P}(n^{-1/2+\epsilon})\\}\text{ as
$n,p\to\infty$.}$ (S8.8)

###### Proof.

Define $\bm{R}=n^{-1}\bm{C}^{\top}P_{X}^{\perp}\bm{C}$ and let $\epsilon>0$ be
an arbitrarily small constant. Then

$\displaystyle\bm{A}_{g}=\tilde{\bm{C}}^{\top}\bm{P}_{g}^{\perp}\tilde{\bm{C}}=$
$\displaystyle\bm{R}^{-1/2}\\{n^{-1}\bm{C}^{\top}P_{X}^{\perp}\bm{W}_{g}P_{X}^{\perp}\bm{C}$
$\displaystyle-$