# High-dimensional log contrast models with measurement errors

Wenxi Tan, Lingzhou Xue
Department of Statistics, Penn State University
Songshan Yang
Center for Applied Statistics and Institute of Statistics and Big Data,
Renmin University of China
Xiang Zhan
Department of Biostatistics, Beijing International Center
for Mathematical Research and Center for Statistical Sciences,
Peking University

###### Abstract

High-dimensional compositional data are frequently encountered in many fields
of modern scientific research. In regression analysis of compositional data,
the presence of covariate measurement errors poses grand challenges for
existing statistical error-in-variable regression analysis methods since
measurement error in one component of the composition has an impact on others.
To simultaneously address the compositional nature and measurement errors in
the high-dimensional design matrix of compositional covariates, we propose a
new method named Error-in-composition (Eric) Lasso for regression analysis of
corrupted compositional predictors. Estimation error bounds of Eric Lasso and
its asymptotic sign-consistent selection properties are established. We then
illustrate the finite sample performance of Eric Lasso using simulation
studies and demonstrate its potential usefulness in a real data application
example.

Keywords: Compositional data; Computational biology; Error-in-variable; Log
contrast models; Lasso

##  1 Introduction

Compositional data, representing relative proportions or percentages of
different parts that make up a whole, have a wide range of applications in
many fields, including geology, ecology, social sciences and biology. In
biological and biomedical research, compositional data primarily arise from
high-throughput sequencing technologies-based profiling experiments, which
share a similar measurement process in which the total abundance information
is lost and most sequence counts reflect only the relative abundances (i.e.,
compositional) information of unique sequences of interest (Vandeputte et al.,
2017). Regression analysis with these compositional covariates is essential to
disentangle the relationships between compositions and an outcome of interest.
Due to compositionality, traditional linear regression models fail for
regression analysis with compositional predictors. To address the “curse of
compositionality” in regression analysis, the log contrast model was proposed
in the context of experiments with mixtures (Aitchison and Bacon-Shone, 1984).
Since then, multiple extensions have been developed, especially in the high-
dimensional settings (Hron et al., 2012; Lin et al., 2014; Shi et al., 2016;
Wang and Zhao, 2017; Randolph et al., 2018; Srinivasan et al., 2021; Combettes
and Müller, 2021; Mishra and Müller, 2022; Shi et al., 2022).

Much of the existing work on high-dimensional log contrast regression has
focused on the clean data case. However, measurement errors are ubiquitous in
many scientific endeavors. Taking the compositional sequence count data in
biomedical research that motivate our study as an example, measurement errors
may occur at any stage of the experimental workflow, such as DNA extraction,
PCR amplification, sequencing process and even bioinformatics preprocessing
procedure (McLaren et al., 2019). Moreover, studies have reported that the
problem of data contamination may also be related to genome databases with a
large number of mis-labeled sequences, which could potentially lead to a wrong
sequencing read counts inflated by orders of magnitude (Gihawi et al., 2023).
These measurement errors need to be well accommodated in statistical analysis
in order to avoid potential misleading or invalid scientific findings (Gihawi
et al., 2023).

These measurement errors in sequencing studies are often referred to as the
sequence bias by many authors working in the field of computational biology
(McLaren et al., 2019; Zhao and Satten, 2021). In those sequencing
experiments, when the research of interest is an individual count variable,
the issue of sequencing bias can be appropriately addressed or attenuated by
multilevel modeling techniques, such as the Beta-Binomial regression (Martin
et al., 2020) or the Poisson-Gamma model (Jiang et al., 2023). However, beyond
marginal analysis, it is much more difficult to analyze contaminated
compositional predictors normalized from multiple count variables (Shi et al.,
2022) or cell-type proportions associated with uncertainties (Cai et al.,
2022), since measurement error in one component has a ripple effect on other
components due to the compositional constraint. Similar to the regression
analysis, we are facing the “curse of compositionality” again in measurement
error modelling for compositional data.

The canonical model for high-dimensional regression analysis is expressed as
${\bf y}={\bf X}\bm{\beta}^{*}+\bm{\epsilon}$, where ${\bf
y}=(y_{1},\ldots,y_{n})^{T}$ is the response vector, ${\bf X}\in R^{n\times
p}$ is the design matrix of high-dimensional covariates and $\bm{\beta}^{*}\in
R^{p}$ is the regression coefficient vector of interest. In many applications,
${\bf X}$ may not be accurately measured and a corrupted version ${\bf Z}$ of
${\bf X}$ is often available. In the literature of high-dimensional
statistics, many versatile methods and theories have been developed for
inference of $\bm{\beta}^{*}$ based on ${\bf Z}$ (Rosenbaum and Tsybakov,
2010; Loh and Wainwright, 2012; Belloni et al., 2017; Datta and Zou, 2017). In
terms of compositional data, each row of ${\bf X}$ and ${\bf Z}$ belongs to
simplex
$\mathcal{S}^{p}=\\{(x_{1},\ldots,x_{p}):x_{j}>0,\sum_{j=1}^{p}x_{j}=1\\}$.
Clearly, measurement error in one component has a ripple effect on other
components due to the compositional constraint and thus, existing statistical
methods and theories are not directly applicable to measurement error problems
of compositional data.

A more recent paper considered the measurement error problem in the framework
of log contrast models for regression analysis with compositional predictors
(Shi et al., 2022). In particular, the variable correction regularized
estimator proposed by Shi et al. (2022) requires the knowledge of observed
counts and leverages the Direchlet-Multinomial distribution to correct
observed counts towards unobservable underlying compositions. However, as
pointed out in a recent paper (Gihawi et al., 2023), sequence read counts
could be inflated by many orders of magnitude due to potential contamination
issues of draft reference genomes. Assumptions of the variable correction
regularized estimator no longer hold when measurement errors in counts are
extremely huge. Moreover, it is sometimes less meaningful to analyze counts
from different platforms than compositions (Allali et al., 2017), which limits
the applicability of variable correction regularized estimator to large cohort
studies where samples are typically sequenced at different locations/batches
using different platforms. These potential limitations motivate our
investigation to take a different approach to fill this research gap.

To develop a new method addressing the aforementioned limitations, we utilize
the recently proposed mathematical model to characterize the issue of sequence
bias in next-generation sequencing experiments (McLaren et al., 2019; McLaren
et al., 2022). To further accommodate compositionality in the design matrix,
we build our regression analysis framework upon the Aitchison log contrast
model (Aitchison, 1982; Aitchison and Bacon-Shone, 1984) and then propose a
new method named textbfError-in-compositional (Eric) Lasso to handle corrupted
high-dimensional compositional predictors. The error bound of our Eric Lasso
estimator along with its selection sign consistency property is established.
In summary, we propose a new method to simultaneously address both
compositional nature and measurement errors in regressors, which distinguishes
it from existing ones. The novelty of our work lies in both a new methodology
with desirable statistical properties and interpretations, but also the
particular application which is an important and timely problem for which no
satisfactory analysis methods exists so far.

The rest of this article is organized as follows. In Section 2, we first
introduce some background on regression analysis with compositional covariates
and on existing literature about bias correction for error-in-variable in log
contrast regression. Then, we propose our Error-in-compositional (Eric) Lasso
method to handle regression analysis with contaminated compositional
covariates. The estimation error bounds of our method and its asymptotic sign
consistency selection properties are established in Section 3. We next
demonstrate the superior performance of our method both using simulation
studies and real data analysis application examples in Section 4 and Section
5, respectively. This article concludes with discussion in Section 6. Proofs
of theoretical results are provided in online supplementary materials.

##  2 Methods

###  2.1 Preliminaries

Let $y_{i}$ and $(U_{i1},\ldots,U_{ip})$ denotes response of interest and $p$
compositional covariates (i.e., $\sum_{j=1}^{p}U_{ij}$=1 and $U_{ij}>0,\forall
j$) measured from the $i$th sampling unit. To study relationships between
response and compositional predictors, the following linear log contrast model
(Aitchison and Bacon-Shone, 1984) has been widely used:

$y_{i}=\sum_{j=1}^{p}log(U_{ij})\beta^{*}_{j}+\epsilon_{i},\>s.t.,\>\sum_{j=1}^{p}\beta^{*}_{j}=0,$
(1)

where the intercept term is omitted if both the outcome and predictors are
centered. One most remarkable feature in model (1) is the zero-sum constraint
on regression coefficients, which is essential to guarantee some basic
principles in compositional data analysis: scale invariance, permutation
invariance and subcompositional coherence (Lin et al., 2014; Greenacre et al.,
2023). These principles are often necessary for statistically meaningful
interpretations of compositional data analysis results (Billheimer et al.,
2001). To understand this, one can first see that each individual component
$j$ of the compositional vector only carries relative information and thus
$\beta_{j}^{*}$ itself is less meaningful in compositional data analysis. On
the other hand, contrast between two coefficients
$\beta_{j}^{*}-\beta_{l}^{*}$ is meaningful in that it can measure the
relative importance of component $j$ and $l$. To get rid of potential bias in
selecting a specific reference level $l$, one can use the quantity
$\frac{1}{p}\sum_{l=1}^{p}(\beta^{*}_{j}-\beta^{*}_{l})$ to measure the
relative importance of component $j$ compared to the remaining ones in the
compositional vector. This quantity reduces to $\beta^{*}_{j}$ if and only if
$\sum_{l=1}^{p}\beta^{*}_{l}=0$ holds. In other words, coefficient
$\beta^{*}_{j}$ can measure the relative importance of the $j$th component of
the compositional vector under this zero-sum constraint. Therefore, it is
crucial to develop interpretable regression analysis methods for compositional
data under this zero-sum constraint in practice (Aitchison, 1982; Lin et al.,
2014; Srinivasan et al., 2021; Shi et al., 2022).

In many scenarios, compositions $U_{ij}$ are not observable or measured with
errors. For example, in a study on associations between gut microbial
compositions and body mass index (Wu et al., 2011), the gut microbial
compositions are not measurable and are approximated the microbiota
compositions in stool samples. In a typical sequencing study (e.g., 16S rRNA
microbiome surveys or single-cell RNA-seq studies) that motivates our
research, the observed compositions $O_{ij}$’s are often calculated from
observed sequence counts $W_{ij}$’s. That is, $O_{ij}=W_{ij}/N_{i}$, where
$N_{i}=\sum_{j=1}^{p}W_{ij}$ is the total counts (or sequencing depth) of
sample $i$. These $W_{ij}^{\prime}s$ are often called original scale of
measurements for compositions. By treating these original scale of
measurements as random realizations of a certain distribution with
compositions being its parameters, it has been argued that modelling these
original scale of measurements has advantages over Aitchison’s log-ratio
approaches (Firth and Sammut, 2023). However, this approach (Firth and Sammut,
2023) often treat compositional measurements as response variables, which does
not apply to our setting of regression with corrupted compositional predictors
as explanatory variables in log contrast regression models considered in the
current article. Since true compositions $U_{ij}$’s may not be directly
measurable, a naive idea is to run a surrogate linear log-contrast regression
model with observed measurements $W_{ij}$’s:

$y_{i}=\sum_{j=1}^{p}log(W_{ij})\beta_{j}+\epsilon_{i},\>s.t.,\>\sum_{j=1}^{p}\beta_{j}=0.$
(2)

Clearly, existence of measurement errors will cause departure of estimated
coefficients $\hat{\bm{\beta}}=(\hat{\beta}_{1},\ldots,\hat{\beta}_{p})$ of
model (2) from the true values
$\bm{\beta}^{*}=(\beta^{*}_{1},\ldots,\beta^{*}_{p})$ of interest in model
(1). Using a bias correction approach, Shi et al. (2022) has shown that the
bias term $\hat{\bm{\beta}}-\bm{\beta}^{*}$ can be improved if $W_{ij}$ in
model (2) is replaced by the corrected variables
$W^{c}_{ij}=W_{ij}+\frac{N_{i}+\alpha_{i}+1}{2\alpha_{i}+1}$, where
$\alpha_{i}$ is the over-dispersion parameter associated with sample $i$ in
the Dirichlet-Multinominal distribution (Shi et al., 2022). In other words,
the authors recommend to use the following model (3) to get more accurate
inference on $\bm{\beta}^{*}$:

$y_{i}=\sum_{j=1}^{p}log(W^{c}_{ij})\beta_{j}+\epsilon_{i},\>s.t.,\>\sum_{j=1}^{p}\beta_{j}=0.$
(3)

The estimator solved from (3) requires the knowledge of observed counts and
leverages the Direchlet-Multinomial distribution to correct observed counts
towards unobservable underlying compositions. However, as pointed out in a
recent paper (Gihawi et al., 2023), sequence read counts could be inflated by
many orders of magnitude due to potential contamination issues of draft
reference genomes. Assumptions of the previous method (Shi et al., 2022) may
no longer hold when measurement errors in counts are extremely huge. Moreover,
it is sometimes less meaningful to analyze counts from different platforms
than compositions (Allali et al., 2017), which limits the applicability of
model (3) to large cohort studies where samples are typically sequenced at
different locations/batches using different platforms. Finally, sign
consistency property of the variable correction regularized estimator has not
been established yet (Shi et al., 2022), which further motivates our
investigation to take an alternative approach to fill this research gap.

###  2.2 Log contrast models with measurement errors

A notable nature of measurement errors in compositional covariates is the
ripple effect, that is, measurement error in one component has an impact on at
least one of the other components due to the unit-sum constraint on
compositions. Unfortunately, most existing high-dimensional error-in-variable
regression methods and theories have focused on the unconstrained data
(Rosenbaum and Tsybakov, 2010; Loh and Wainwright, 2012; Datta and Zou, 2017),
which are not directly applicable to compositional data. In a typical sequence
count technologies-based microbiome study, discrepancies between $U_{ij}$ and
$O_{ij}$ (or simply sequencing bias) are thought to approximately act
multiplicatively on the taxon abundances (McLaren et al., 2019). Under this
assumption, the following mathematical model has been proposed by
computational biologists to characterize the sequence bias issue of microbiome
compositional data collected from next generation sequencing experiments
(Clausen and Willis, 2022; McLaren et al., 2022):

$O_{ij}=U_{ij}\cdot\frac{E_{ij}}{\sum_{j=1}^{p}U_{ij}E_{ij}},$ (4)

where $E_{ij}$ denotes the measurement error term for sample $i$ and taxon
$j$. That is, measurement errors are both sample-specific and taxon-specific.
This is because, on the one hand, taxa are not all detected equally well and
measurement error is determined by the interaction between experimental
protocols and the biological/chemical/physical state of that taxon. On the
other hand, samples might come from different sources and be treated by
different specimens or technicians and thus measurement error may also depend
on the sampling process (McLaren et al., 2022). For these reasons, we assume
that measurement errors depend on taxon characteristics and the sampling
process, but not on its relative abundance. That is, we assume that $E_{ij}$
and $U_{ij}$ are independent.

The compositional nature of measurements has been clearly addressed in model
(4) since $\sum_{j=1}^{p}O_{ij}=1$. That is, ripple effects of measurement
errors in one compositional component to another are well accommodated in
model (4). Also, based on model (4), we have
$\sum_{j=1}^{p}\beta_{j}log(O_{ij})=\sum_{j=1}^{p}\beta_{j}log(U_{ij})+\sum_{j=1}^{p}\beta_{j}log(E_{ij})$
holds for any $\beta_{1},\ldots,\beta_{p}$ with $\sum_{j=1}^{p}\beta_{j}=0$.
For ease of presentation, we define
$X_{ij}=log(U_{ij}),Z_{ij}=log(O_{ij}),B_{ij}=log(E_{ij})$. Then, in the scale
of log contrasts, we have ${\bf Z}\bm{\beta}=({\bf X}+{\bf B})\bm{\beta}$,
where $\bm{\beta}$ is an arbitrary zero-sum vector. We emphasize that the
formulation ${\bf Z}={\bf X}+{\bf B}$ is only true under the log contrast
regression model and this formulation will largely facilitate measure error
modelling for regression analysis with contaminated compositional predictors.
In other word, we do have “blessing of compositionality” in measurement error
modelling for log contrast regression with compositional covariates. We
further assume that rows of ${\bf B}$ are independent and identically
distributed with zero mean, finite covariance $\bm{\Sigma}_{B}$ and sub-
Gaussian parameter $\tau^{2}$ for the purpose of developing further
statistical inference.

Let $S=\\{j:\beta^{*}_{j}\neq 0\\}$ denote the support of $\bm{\beta}^{*}$ and
$s=|S|$ is the cardinality of $S$. For any subset $J\subset\\{1,\ldots,p\\}$
and any index $j\in J$, let $J^{c}$ denote the complement of $J$ and define
$J-j=J\backslash\\{j\\}$. For any matrix $K\in\mathbb{R}^{a\times b}$, $K_{J}$
denotes the submatrix of the $j$th column for $j\in J\cap\\{1,\dots,b\\}$ of
$K$, $K_{J_{1}J_{2}}$ is the submatrix formed by $(i,j)$th entries for $i\in
J_{1}\cap\\{1,\dots,a\\},j\in J_{2}\cap\\{1,\dots,b\\}$. $\|\cdot\|_{\infty}$
and $\|\cdot\|_{1}$ denote the $l_{\infty}$ norm and $l_{1}$ norm
respectively. Then, our log contrast model is expressed as: ${\bf y}={\bf
X}\bm{\beta}^{*}+\bm{\epsilon}={\bf
X}_{S}\bm{\beta}^{*}_{S}+\bm{\epsilon},\>s.t.,\>\sum_{j=1}^{p}\beta^{*}_{j}=\sum_{j\in
S}\beta^{*}_{j}=0$, where $\bm{\epsilon}=(\epsilon_{1},\ldots,\epsilon_{n})$
are independent and identically distributed random errors that follow
$N(0,\sigma^{2})$. While predictors ${\bf X}$ may not be observable, we
propose a new method that use its surrogate ${\bf Z}$ to perform inference on
$\bm{\beta}^{*}$. We further introduce some necessary notation before
introducing our new method. Let ${\bf X}^{p}=(log(U_{ij}/U_{ip}))\in
R^{n\times(p-1)}$ and ${\bf Z}^{p}=(log(O_{ij}/O_{ip}))\in R^{n\times(p-1)}$.
Without loss of generality, we pick the last component as the reference level
when defining ${\bf X}^{p}$ and ${\bf Z}^{p}$ in order to facilitate
development and presentation of our method. However, in the subsequent
methodology development, we will guarantee that the proposed method is
invariant on the selection of the reference component. Corresponding to these
log ratios, we can partition the $p$-dimensional regression coefficients as
$\bm{\beta}=(\bm{\beta}_{-p},\beta_{p})$, where
$\bm{\beta}_{-p}=(\beta_{1},\ldots,\beta_{p-1})$ is free of any constraint.
Assuming all columns in ${\bf X}$ are centered and define
$\bm{\Sigma}=\frac{1}{n}{\bf X}^{T}{\bf X}$. Then, one way to fit a sparse log
contrast model in the high-dimensional setting is via the following
$L_{1}$-regularization:

$\displaystyle\hat{\bm{\beta}}$
$\displaystyle=\mathop{\arg\min}_{\bm{\beta}}\frac{1}{2n}\|{\bf y}-{\bf
X}\bm{\beta}\|_{2}^{2}+\lambda\|\bm{\beta}\|_{1}\>,s.t.,\>\sum_{j=1}^{p}\beta_{j}=0$
$\displaystyle=\mathop{\arg\min}_{\bm{\beta}_{-p}}\frac{1}{2n}\|{\bf y}-{\bf
X}^{p}\bm{\beta}_{-p}\|_{2}^{2}+\lambda\|{\bf D}^{p}\bm{\beta}_{-p}\|_{1}$
$\displaystyle=\mathop{\arg\min}_{\bm{\beta}_{-p}}\frac{1}{2}\bm{\beta}_{-p}^{T}\bm{\Sigma}^{p}\bm{\beta}_{-p}-(\bm{\rho}^{p})^{T}\bm{\beta}_{-p}+\lambda\|{\bf
D}^{p}\bm{\beta}_{-p}\|_{1},$ (5)

where $\bm{\Sigma}^{p}=\frac{1}{n}({\bf X}^{p})^{T}{\bf X}^{p}$,
$\bm{\rho}^{p}=\frac{1}{n}({\bf X}^{p})^{T}{\bf y}$ and ${\bf D}^{P}=({\bf
I}_{p-1},-\mathbf{1}_{p-1})^{T}$.

Since ${\bf X}$ is unobserved, we have to choose surrogates
$\tilde{\bm{\Sigma}}^{p}$ and $\tilde{\bm{\rho}}^{p}$ to replace
$\bm{\Sigma}^{p}$ and $\bm{\rho}^{p}$ in (2.2). We begin by first constructing
$\hat{\bm{\Sigma}}$, an unbiased estimator of $\bm{\Sigma}$. As we specify the
observed design matrix ${\bf Z}$ being contaminated by additive measurement
error ${\bf B}$, where the rows of ${\bf B}$ are independent and identically
distributed with zero mean and finite covariance $\bm{\Sigma}_{B}$, then
$\hat{\bm{\Sigma}}=\frac{1}{n}{\bf Z}^{T}{\bf Z}-\bm{\Sigma}_{B}$ is an
unbiased estimator of $\bm{\Sigma}$. Then, following suggestion of CoCoLasso
(Datta and Zou, 2017), The surrogate $\tilde{\bm{\Sigma}}$ and
$\tilde{\bm{\rho}}$ for the unobserved $\bm{\Sigma}=\frac{1}{n}{\bf X}^{T}{\bf
X}$ and $\bm{\rho}=\frac{1}{n}{\bf X}^{T}{\bf y}$ are then defined as:

$\tilde{\bm{\Sigma}}=\mathop{\arg\min}_{{\bf K}\geq 0}\|{\bf
K}-\hat{\bm{\Sigma}}\|_{\max}\>\>,\>\>\tilde{\bm{\rho}}=\frac{1}{n}{\bf
Z}^{T}{\bf y},$

where $\|{\bf K}\|_{\max}$ denotes the element-wise maximum norm of ${\bf K}$
and $\tilde{\bm{\Sigma}}$ is obtained through ADMM algorithm. Through matrix
transformation, we have $\bm{\Sigma}^{p}=\frac{1}{n}({\bf X}^{p})^{T}{\bf
X}^{p}=({\bf D}^{p})^{T}\bm{\Sigma}{\bf D}^{p}$. Then, surrogates
$\tilde{\bm{\Sigma}}^{p}$ and $\tilde{\bm{\rho}}^{p}$ is calculated from

$\tilde{\bm{\Sigma}}^{p}=({\bf D}^{p})^{T}\tilde{\bm{\Sigma}}{\bf
D}^{p}\>\>,\>\>\tilde{\bm{\rho}}^{p}=({\bf D}^{p})^{T}\tilde{\bm{\rho}}.$

Finally, plugging these quantities calculated from observed measurements ${\bf
Z}$ into (2.2), we obtain the following estimator

$\hat{\bm{\beta}}=\mathop{\arg\min}_{\bm{\beta}\,s.t.\bm{\beta}={\bf
D}^{p}\bm{\beta}_{-p}}\frac{1}{2}\bm{\beta}_{-p}^{T}\tilde{\bm{\Sigma}}^{p}\bm{\beta}_{-p}-(\tilde{\bm{\rho}}^{p})^{T}\bm{\beta}_{-p}+\lambda\|{\bf
D}^{p}\bm{\beta}_{-p}\|_{1}.$ (6)

We term this new estimator (6) as the Error-in-composition Lasso estimator or
simply Eric Lasso hereafter in this paper.

_Remark 1_ : The core of CoCoLasso lies in the two-step procedure of
constructing a good estimator for matrix $\bm{\Sigma}$: first getting an
unbiased estimation $\hat{\bm{\Sigma}}$ and then conducting projection to
obtain a positive semi-definite matrix $\tilde{\Sigma}$. In our context of
regression analysis with compositional predictors, the most natural approach
is to directly apply this two-step procedure to log-ratio transformed
measurements ${\bf Z}^{p}$ to obtain surrogates $\tilde{\bm{\Sigma}}^{p}$ and
$\tilde{\bm{\rho}}^{p}$ of ${\bf X}^{p}$, which unfortunately leads to
analysis results that are sensitive to selection of the reference component
and thus lacks valid interpretation for compositional data analysis. In
contrast, the proposed Eric Lasso method is permutation invariant, which gives
the same result under any permutation of the $p$ components. Throughout all
assumptions and proofs in the subsequent parts of this article, we have
intentionally ensured that conditions and conclusions are all independent from
the selection of the reference component to enhance statistical
interpretations of our compositional data analysis.

_Remark 2_ : The error covariance $\bm{\Sigma}_{B}$ might be unknown in
practice, and must be obtained through estimation in such a case. Suppose we
can borrow information from either independent external data or replicated
data to calculate an error matrix $B_{o}\in\mathbb{R}^{n\times p}$ and
correspondingly $\hat{\bm{\Sigma}}_{B}=\frac{1}{n}B_{o}^{T}B_{o}$ as an
estimate of $\bm{\Sigma}_{B}$, which has also been widely assumed in
literature (Loh and Wainwright, 2012; Shi et al., 2022). We can show that our
theoretical analysis still holds when replacing $\bm{\Sigma}_{B}$ with
$\hat{\bm{\Sigma}}_{B}$ in Lemma 1 in Section A of the online supplementary
materials.

##  3 Theoretical analysis

We assume without loss of generality that $p\in S$ such that $S-p$ is well-
defined. Otherwise, we could permute the compositional covariates to ensure
it. This assumption is reasonable as long as we can guarantee all conditions
and proofs are invariant to permutations of covariates indices, which has been
well checked throughout this article. To establish the theoretical results of
Eric Lasso, we need the following two regularity conditions:

######  Condition 1.

The matrix $\Sigma_{SS}$ satisfies

$\Lambda_{\min}(\Sigma_{SS})\geq C_{\min}>0,$ (7)

where $\Lambda_{\min}(K)$ denotes the minimal eigenvalue of matrix $K$.

######  Condition 2.

There exists some $\xi\in(0,1]$ such that

$\|\Sigma^{p}_{S^{c}S}(\Sigma^{p}_{SS})^{-1}\\{\mbox{sgn}(\beta^{*}_{S-p})-\mbox{sgn}(\beta_{p}^{*})1_{s-1}\\}+\mbox{sgn}(\beta_{p}^{*})1_{p-s}\|_{\infty}\leq
1-\xi.$ (8)

Condition 1 is a common assumption in high-dimensional statistics and some
implications to guarantee the permutation invariance of our analysis under
this condition are derived in Lemma 3 and Lemma 4 of the online supplementary
materials. Condition 2 is taken from Lin et al. (2014), which is central to
guaranteed support recovery of $L_{1}$ regularization. An important quantity
associated with Condition 2 is
$\phi=\|D^{p}_{SS}(\Sigma^{p}_{SS})^{-1}(D^{p}_{SS})^{T}\|_{\infty}$.
According to Proposition 2 of Lin et al. (2014), Condition 2 and $\phi$ are
permutation invariant. That is, $\phi$ and the left-hand part of Equation (8)
are both independent from the selection of reference component $p$. Moreover,
(7) is permutation invariant naturally. So far, all our conditions are
permutation invariant. We further assume that $\max_{j}\|X_{j}\|_{2}^{2}\leq
n$. With all these regulatory conditions, we now are ready to establish the
main theoretical results on our Eric Lasso estimator.

######  Theorem 1 (error bound and sign consistency).

Under Condition 1 and Condition 2, let $\zeta=\max(\tau^{4},\sigma^{4},1)$,
for $\lambda\leq\tau^{2}$ and
$\epsilon\leq\min(\epsilon_{1},\lambda/(\lambda\epsilon_{2}+\epsilon_{3}))$
where $\epsilon_{i}$’s are bounded positive constants depending of $\Sigma$,
$\beta^{*}_{S}$, $\theta$ and $\phi$, Then there exists universal constants
$C$ and $c$, with probability at least $1-\delta_{1}$ where
$\delta_{1}=p^{2}C\exp\left(-cn(s-1)^{-4}\epsilon^{2}\zeta^{-1}\right)+pC\exp\left(-cns^{-2}\lambda^{2}\xi^{2}\zeta^{-1}\right)$,
problem(6) has an optimal solution $\hat{\bm{\beta}}$ that satisfies the
following properties: (a) $l_{\infty}$-loss:
$\|\hat{\beta}_{S}-\beta^{*}_{S}\|_{\infty}\leq 9\phi\lambda/2$, (b) sign
consistency: if $\min_{j\in S}|\beta_{j}^{*}|>9\phi\lambda/2$, then
$\mbox{sgn}(\hat{\beta})=\mbox{sgn}(\beta^{*})$.

_Remark 3_ : In CoCoLasso, the parameter $\delta_{1}$ goes to zero when
$s^{2}\log p/n\to 0$ as $n,p\to\infty$ and we need $(s-1)^{4}\log p/n\to 0$ in
Eric Lasso. This discrepancy is due to the fact that the estimation of
$\beta_{S}$ is actually done via estimation of $\beta_{S-p}$ in a
$(s-1)$-dimensional space and then $\beta_{p}=-\sum_{j\neq p}\beta_{j}$ is
obtained as a transformation. As a result, we have the inequality:
$\|\hat{\beta}_{S}-\beta^{*}_{S}\|_{\infty}\leq(s-1)\|\hat{\beta}_{S-p}-\beta^{*}_{S-p}\|_{\infty}$.
In the same framework as CoCoLasso, we can only bound $\beta_{S-p}$, leading
to an error term with an extra factor of $(s-1)$.

_Remark 4_ : To understand the asymptotic implications of Theorem 1, we assume
for simplicity that $\phi$ is constant. From the expression of $\delta_{1}$,
If $(s-1)^{4}\log p/n\to 0\text{ as }n,p\to\infty$ and $\min_{j\in
S}|\beta_{j}^{*}|\gg s(\zeta\log p/n)^{1/2}$, then we can choose a $\lambda$
satisfying both $\lambda\gg s(\zeta\log p/n)^{1/2}$ and $\min_{j\in
S}|\beta_{j}^{*}|>9\phi\lambda/2$ such that $\delta_{1}$ goes to zero and the
sign-consistency of Eric Lasso is achieved.

##  4 Simulation Studies

We have conducted comprehensive numerical studies to evaluate the performance
of the proposed Eric Lasso method. The goals of our simulation studies are: 1)
to provide numerical evidence supporting the theoretical results established
for Eric Lasso estimator in Theorem 1. 2) to assess robustness of the proposed
method using simulated data generated from a variety of different scenarios
including misspecified models.

###  4.1 Simulation setup

We used three different approaches to generate compositional predictors to
evaluate the performance of the proposed method. In the first approach, we
followed the simulation design of a previous paper (Lin et al., 2014) to use
the logistic normal distribution (Aitchison and Shen, 1980) to simulate
unobserved community compositions. In particular, we first simulated a
$n\times p$ latent matrix denoted as $W$ from a multivariate normal
distribution $N(\bm{\theta},\bm{\Sigma}_{W})$ where
$(\bm{\Sigma}_{W})_{ij}=\rho^{|i-j|}$, where $\rho=0.5$ and $\bm{\theta}$ was
set in the following way: $\theta_{j}=\log(0.2p)$ for j = 1 to 5, while for
other indices j, $\theta_{j}$ was set to 0. The unobserved community
compositions were further calculated as
$U_{ij}=\exp(W_{ij})/\sum_{k=1}^{p}\exp(W_{ik})$ and $X_{ij}=log(U_{ij})$.
Compositions generated under this scheme were very heterogeneous in that the
first five components dominate the remaining components. The first a few
coordinates of the true coefficient vector $\bm{\beta}^{*}$ were specified as
[1.2, -0.8, 0.7, 0, 0, -1.5, -1, 1.4] with remaining components being zeros.
The response was then generated as ${\bf y}={\bf X}\bm{\beta}+\bm{\epsilon}$,
where $\epsilon_{1},\ldots,\epsilon_{n}$ are iid errors simulated from
$N(0,0.5^{2})$. Since community compositions are often unobserved, we
simulated a corrupted version by adding errors. In particular, we first
generated an error matrix ${\bf B}$, whose rows were independently simulated
from $N({\bf 0},\tau^{2}{\bf I})$ and then generated the corrupted
compositions $O_{ij}$ according to model (4), which were further transformed
into $Z_{ij}=log(O_{ij})$. In the second simulation scenario, we considered
more homogeneous community compositions by using the Dirichlet distribution to
generate community compositions (Fiksel et al., 2022). In particular, each row
of ${\bf U}$ was independently simulated from the Dirichlet ($\frac{1}{p}{\bf
1}_{p}$), where ${\bf 1}_{p}$ denotes the $p$-dimensional vector of all ones
and we kept all other settings the same as the first scenario. In the third
simulation scenario, we examined method performance under misspecified models
violating our core model assumptions as stated in Equation (4). In particular,
we followed the exact simulation design of variable correction regularized
estimator (Shi et al., 2022) to generate latent counts. We first generated
underlying compositions ($U_{i1},\ldots,U_{ip}$) from the above logistic
normal distribution and total sequence counts $N_{i}$ of sample $i$ from the
negative binomial distribution with mean $3\times 10^{4}$ and variance
$3\times 10^{6}$. Then, observed counts $(W_{i1},\dots,W_{ip})$ were simulated
from DirMult$(N_{i},\alpha U_{i1},\dots,\alpha U_{ip})$ with $\alpha=5000$,
where DirMult denotes Dirichlet-Multinomial distribution. Following suggestion
of Shi et al. (2022) on handling zero counts, we calculated the observed
design matrix ${\bf Z}$ as $Z_{ij}=log((W_{ij}+0.5)/\sum_{j}(W_{ij}+0.5))$.
For ease of presentation, we term these three different data generation
mechanisms as Scenario 1, 2 and 3, respectively, hereafter.

Simulation configuration parameters $(n,p,\tau)$ have not been specified yet.
We followed simulation designs of two major competitors (Lin et al., 2014;
Datta and Zou, 2017) of Eric Lasso to evaluate the performance of our method
in a comprehensive manner. In Simulation I, we fixed $\tau=0.5$ and varied
$(n,p)=(100,200),(250,400),(500,500)$ to evaluate the performance of our
method under different data sizes. This type of setting was considered in the
compositional Lasso paper (Lin et al., 2014). In Simulation II, we fixed
$(n,p)=(100,100)$ to mimic a real data set analyzed later in this paper and
varied $\tau$ to evaluate the robustness of the proposed method with respect
to the noise level of measurement errors. This type of setting was considered
in the CoCoLasso paper (Datta and Zou, 2017). Since Scenario 3 does not
contain $\tau$ as a configuration parameter, it was only considered under
Simulation I.

After each dataset was simulated under a particular scenario and setting. We
applied four different methods to data $({\bf y},{\bf Z})$ to obtain the
regression coefficient estimator $\hat{\bm{\beta}}$ and then compared it to
true regression coefficients
$\bm{\beta}^{*}=(1.2,-0.8,0.7,0,0,-1.5,-1,1.4,0,\ldots,0)$. The four methods
are the proposed Eric Lasso estimator, the compositional data analysis version
of Lasso estimator (Lin et al., 2014), the CoCoLasso estimator (Datta and Zou,
2017) and the vanilla Lasso (Tibshirani, 1996). Among these methods,
compositional Lasso incorporates the compositional nature in data but ignores
measurement errors. CoCoLasso accommodates measurement errors yet fails to
model the compositional constraint in the regression covariates space. The
original Lasso method fails to accommodate neither characteristics of the data
and the proposed Eric Lasso method is the only one to take both aspects into
account. For ease of presentation, we refer to these Lasso methods as Eric,
Coda, CoCo and Vani, respectively, hereafter in this article. The tuning
parameter selection in Eric and CoCo Lasso was done by using the calibrated
cross validation method (Datta and Zou, 2020) and that in Coda and Vani Lasso
was done by using cross validation.

To compare estimation accuracy of different methods, we calculated the
following three metrics based on $\bm{\beta}^{*}-\hat{\bm{\beta}}$ for each
method, including the squared error (SE), the prediction error (PE) and the
$l_{\infty}$ loss:

$SE(\hat{\bm{\beta}})=\|\bm{\beta}^{*}-\hat{\bm{\beta}}\|_{2}^{2},\>\>PE(\hat{\bm{\beta}})=(\bm{\beta}^{*}-\hat{\bm{\beta}})^{T}\bm{\Sigma}(\bm{\beta}^{*}-\hat{\bm{\beta}}),\>\>l_{\infty}(\hat{\bm{\beta}})=\max|\bm{\beta}^{*}-\hat{\bm{\beta}}|.$

To compare the sign consistency on selection of different methods, we
calculated the following two metrics based on
$sgn(\bm{\beta}^{*})-sgn(\hat{\bm{\beta}})$, including the false negative rate
(FNR) and false positive rate (FPR). These quantities are commonly used to
assess the sensitivity and specificity of a method and are defined as:

$FPR=\frac{\\#\\{j:\hat{\beta}_{j}\neq
0\cap\beta^{*}_{j}=0\\}}{\\#\\{j:\beta^{*}_{j}=0\\}},\>\>FNR=\frac{\\#\\{j:\hat{\beta}_{j}=0\cap\beta^{*}_{j}\neq
0\\}}{\\#\\{j:\beta^{*}_{j}\neq 0\\}}.$

For each specific simulation scenario, we repeated 100 times to obtain
multiple values of these metrics and reported their mean values along with the
standard errors of the mean in the following section.

###  4.2 Simulation results

As discussed in Section 2, model interpretation is important in compositional
data analysis. In particular, an important concern in log contrast regression
is the zero-sum constraint on regression coefficients, which can guarantee
some basic principles in compositional data analysis and has been well
recognized in literature (Lin et al., 2014; Shi et al., 2022). Towards this
end, we first examine whether the sum of regression coefficients estimated by
different Lasso methods is close to zero. For each simulation run, we
calculated the value of $\sum_{j=1}^{p}\hat{\beta}_{j}$ and then reported the
average value over 100 replicates. A t-test was also performed to examine
whether it is significantly different from zero. Results under Scenario 1 of
Simulation I are reported in Table 1. As shown in Table 1, both Eric Lasso and
Coda Lasso obtain estimated coefficients with a sum very close to zero, and
thus provide compositional data analysis with meaningful interpretations. On
the other hand, CoCo Lasso and Vanilla Lasso fail to preserve the zero-sum
constraint on regression coefficients, which makes it more difficult to
interpret the corresponding results in the framework of log contrast
regression models for compositional data analysis. We observe a similar
phenomenon under other scenarios (reported in Section B of the online
supplementary materials).

Table 1: The average sum of regression coefficients estimated by different Lasso methods and corresponding t-test p values under Scenario 1 of Simulation I. (n,p) | Model | $\sum\beta_{j}$ | p value
---|---|---|---
(100,200) | Eric | -2.4e-08 | 0.600
| Coda | -3.2e-08 | 0.406
| CoCo | -3.3e-01 | 1.3e-5
| Vani | -4.1e-01 | 1.9e-10
(250,400) | Eric | -2.2e-08 | 0.483
| Coda | -3.0e-08 | 0.218
| CoCo | -4.2e-02 | 0.381
| Vani | -1.3e-01 | 0.002
(500,500) | Eric | -1.6e-08 | 0.628
| Coda | -1.1e-08 | 0.664
| CoCo | 6.9e-02 | 0.026
| Vani | -1.0e-01 | 4.9e-05

We next compare the estimation and selection performance of difference Lasso
methods. Results under Scenario 1 and 2 of Simulation I are reported in Table
2 and those under Scenario 3 of Simulation I (model misspecification) are
reported in Table 3. As shown in these tables, Eric Lasso and CoCoLasso
consistently have better estimation performance (in terms of SE, PE and
$l_{\infty}$ loss) than those of Coda and Vanilla Lasso. Under all three
scenarios, we observe that estimation errors of Eric Lasso tend to decrease as
the sample size $n$ increases, with the only exception being the PE of Eric
Lasso under Scenario 3. As for the $l_{\infty}$-loss of Eric Lasso implicated
in Theorem 1, we do observe that it tends to vanish as sample size get larger
and larger, even when the data are generated from a misspecified model under
Scenario 3. For selection accuracy, the FPR of Coda and Vani Lasso are
significantly higher than those of Eric and CoCo Lasso under Scenario 1. A
similar phenomenon has been observed in the literature that Coda Lasso tends
to select more unnecessary false positives to recover the true model (Susin et
al., 2020; Srinivasan et al., 2021). On the other hand, differences in FNR of
four methods under Scenario 1 are smaller. Patterns under Scenario 2 are
similar to those under Scenario 1, and overall, Eric Lasso has the best
selection performance under these two scenarios. When models are misspecified
under Scenario 3, a remarkable change is that Eric and CoCo Lasso have much
worse FNR compared to Coda and Vani Lasso, which is not surprising given that
the FPR of Coda and Vani are two to five times of those of Eric and CoCo.
There is no methods that are uniformly better than others in terms of both FPR
and FNR under this scenario. In summary, only two methods (Eric and Coda) can
lead to valid statistical interpretations for log contrast models. Between
this two methods, Eric is consistently better than Coda both in terms of
estimation and selection performance across all three scenarios considered in
Simulation I.

Table 2: Comparison of different Lasso estimators under Scenario 1 (top half)
and 2 (bottom half) of Simulation I. Mean values (standard errors) of
different evaluation metrics are calculated based on 100 simulation
replicates.

(n,p) | Model | SE | PE | $l_{\infty}$ | FPR | FNR
---|---|---|---|---|---|---
| Eric | 2.16(0.08) | 0.77(0.03) | 0.82(0.02) | 0.08(0) | 0.17(0.01)
| Coda | 2.91(0.08) | 1.10(0.04) | 0.89(0.02) | 0.11(0.01) | 0.16(0.01)
(100,200) | CoCo | 2.22(0.10) | 0.75(0.03) | 0.82(0.02) | 0.09(0) | 0.16(0.02)
| Vani | 2.96(0.09) | 1.03(0.04) | 0.89(0.02) | 0.13(0.01) | 0.17(0.02)
| Eric | 1.06(0.03) | 0.43(0.01) | 0.59(0.01) | 0.05(0) | 0.03(0.01)
| Coda | 1.95(0.03) | 0.83(0.01) | 0.73(0.01) | 0.09(0) | 0.03(0.01)
(250,400) | CoCo | 1.03(0.04) | 0.42(0.01) | 0.59(0.01) | 0.06(0) | 0.02(0.01)
| Vani | 1.92(0.03) | 0.78(0.02) | 0.72(0.01) | 0.11(0.01) | 0.03(0.01)
| Eric | 0.59(0.02) | 0.26(0.01) | 0.45(0.01) | 0.03(0) | 0(0)
| Coda | 1.54(0.02) | 0.72(0.01) | 0.65(0.01) | 0.07(0) | 0(0)
(500,500) | CoCo | 0.57(0.01) | 0.26(0.01) | 0.46(0.01) | 0.03(0) | 0(0)
| Vani | 1.54(0.02) | 0.70(0.01) | 0.64(0.01) | 0.08(0) | 0(0)
| Eric | 0.07(0) | 0.73(0.02) | 0.13(0.01) | 0.18(0.01) | 0(0)
| Coda | 0.09(0.01) | 0.84(0.02) | 0.16(0.01) | 0.19(0.01) | 0(0)
(100,200) | CoCo | 0.07(0) | 0.76(0.02) | 0.13(0.01) | 0.20(0.01) | 0(0)
| Vani | 0.09(0.01) | 0.85(0.02) | 0.16(0.01) | 0.21(0.01) | 0(0)
| Eric | 0.06(0.01) | 0.32(0.01) | 0.13(0.01) | 0.09(0.01) | 0(0)
| Coda | 0.09(0.01) | 0.48(0.01) | 0.17(0.01) | 0.11(0.01) | 0(0)
(250,400) | CoCo | 0.06(0.01) | 0.34(0.01) | 0.13(0.01) | 0.11(0.01) | 0(0)
| Vani | 0.09(0.01) | 0.49(0.01) | 0.17(0.01) | 0.13(0.01) | 0(0)
| Eric | 0.03(0) | 0.16(0.01) | 0.09(0) | 0.05(0) | 0(0)
| Coda | 0.06(0) | 0.34(0.01) | 0.13(0) | 0.07(0) | 0(0)
(500,500) | CoCo | 0.03(0) | 0.17(0.01) | 0.09(0) | 0.06(0) | 0(0)
| Vani | 0.06(0) | 0.34(0.01) | 0.13(0) | 0.08(0) | 0(0)

Table 3: Comparison of different Lasso estimators under Scenario 3 of Simulation I. Mean values (standard errors) of different evaluation metrics are calculated based on 100 simulation replicates. (n,p) | Model | SE | PE | $l_{\infty}$ | FPR | FNR
---|---|---|---|---|---|---
| Eric | 2.27(0.05) | 0.86(0.02) | 0.81(0.01) | 0.08(0) | 0.17(0.01)
| Coda | 2.44(0.07) | 0.92(0.03) | 0.83(0.01) | 0.14(0.01) | 0.11(0.01)
(100,200) | CoCo | 2.45(0.05) | 0.86(0.02) | 0.83(0.01) | 0.08(0) | 0.22(0.01)
| Vani | 2.43(0.07) | 0.87(0.03) | 0.82(0.01) | 0.16(0.01) | 0.12(0.01)
| Eric | 2.12(0.04) | 0.99(0.02) | 0.78(0.01) | 0.03(0) | 0.15(0.01)
| Coda | 2.23(0.04) | 1.03(0.02) | 0.85(0.01) | 0.13(0.01) | 0.02(0.01)
(250,400) | CoCo | 2.22(0.04) | 0.98(0.02) | 0.78(0.01) | 0.03(0) | 0.18(0.01)
| Vani | 2.15(0.04) | 0.97(0.01) | 0.84(0.01) | 0.15(0.01) | 0.01(0)
| Eric | 1.51(0.04) | 0.73(0.02) | 0.67(0.01) | 0.02(0.01) | 0.04(0.01)
| Coda | 1.90(0.02) | 1.02(0.01) | 0.85(0.01) | 0.12(0) | 0(0)
(500,500) | CoCo | 1.54(0.04) | 0.72(0.02) | 0.64(0.01) | 0.03(0) | 0.03(0.01)
| Vani | 1.88(0.02) | 1.00(0.01) | 0.85(0.01) | 0.13(0.01) | 0(0)

We next evaluate the robustness of the proposed Eric Lasso method against
different noise levels of measurement errors, which is measured by the $\tau$
parameter used in simulation Scenario 1 and 2 of Simulation II. To achieve
this goal, we use ROC curves to compare performance of different Lasso methods
under different $\tau$ values. The ROC curves of four methods are displayed in
Figure 1. As the error level increases, it is more difficult to detect
findings for all methods such that both FPR and TPR reduce. Clearly, the left
panel shows that Eric Lasso and CoCo Lasso have larger areas under curve (AUC)
than Coda and Vanilla Lasso method, which is consistent with our conclusions
found in Simulation I. Furthermore, the right panel of Figure 1 shows that
Eric has a larger AUC than CoCo Lasso. The same is also true when comparing
Coda and Vani. ROC curves under Scenario 2 show the exact same pattern and are
displayed in Section B of the online supplementary materials. Therefore, under
both scenarios of compositional predictors, the special consideration in Eric
Lasso to accommodate compositionality does improve its discriminative power
over the classic CoCo Lasso method developed for high-dimensional Gaussian
data. If the comparison is limited to methods with valid log contrast model
interpretations, the area under ROC curve of Eric is significantly larger than
that of Coda.

|
---|---
Figure 1: The ROC curve with different $\tau$ values under Scenario 1. The
left panel is ROC curves at original FPR and TPR scales and the right panel
zooms in specific regions of FPR and TPR to better distinguish different
methods.

Combining all numerical results in Simulation I and II, the proposed Eric
Lasso method stands out in obtaining a regression model with not only more
accurate coefficients estimation and selection, but also valid model
interpretation for compositional data regression analysis among the four Lasso
methods being evaluated. Between the two methods (Eric and Coda) with valid
interpretation for log contrast regression model, the performance of Eric is
substantially better than that of Coda. It is of note that all aforementioned
results are still true even when data were generated under a misspecified
model, which demonstrates the robustness of Eric Lasso. Therefore, our new
method is the best candidate for high-dimensional regression analysis of
compositional covariates with measurement errors compared to these existing
potential competitors.

##  5 A case study

The human gut microbiota has been shown to play a very important role in
nutrient digestion and absorption (Wu et al., 2011) and most existing analyses
have successfully shown that there exists a significant association between
the overall gut microbiome community and body mass index (BMI) using
permutational multivariate analysis of variance type of approaches (Tang et
al., 2017). However, identifying specific microbial taxa associated with the
outcome is more challenging than detecting an existing overall community-level
association, partially due to the compositional effect of microbiome data
(individual taxa are closely related to or affected by each other). To
illustrate the potential usefulness of Eric Lasso, we applied it to
investigate BMI associated gut microbial taxa using data collected in the
COMBO study (Wu et al., 2011). A total of 3068 non-singleton operational
taxonomic units (OTUs) were detected in the COMBO study. We first aggregated
these OTUs into the genus level and then deleted genera that appeared in less
than 2 samples, ending with up $p=80$ genera kept in subsequential statistical
analysis. After data filtering and quality control, a total of $n=96$ samples
were kept for further analysis and we followed the previous suggestion (Shi et
al., 2022) to transfer these sequence counts into relative abundances for
further analysis. To adjust for potential confounding effects, we first
regressed BMI on total fat and caloric intake and then took residuals as
outcomes for association analysis with gut microbial compositions.

Let ${\bf y}=(y_{1},\ldots,y_{n})$ denote the covariate-adjusted BMI values,
$U_{ij}$ denote the relative abundance of the $j$th genus in the $i$th sample,
$i=1,\ldots,n,j=1,\ldots,p$. To mimic potential measurement error or bias in
the sequencing procedure, we reserved $U_{ij}$ as the true abundances and used
perturbation to generate corrupted abundances. In particular, for each sample
$i$, we first independently simulated each measurement error factor
$e_{ij}\sim Unif(0.1,10)$ and then calculated corrupted abundances as
$U_{ij}e_{ij}$, which were further normalized into a compositional vector
($O_{i1},\ldots,O_{ip}$) for downstream analysis. Let
$X_{ij}=log(U_{ij}),Z_{ij}=log(O_{ij})$ and ${\bf X}=\\{X_{ij}\\},{\bf
Z}=\\{Z_{ij}\\}$ be the corresponding design matrix. Our goal is to use the
noisy version $({\bf y},{\bf Z})$ to infer the relation of $({\bf y},{\bf
X})$. To achieve this goal, we reserved ${\bf X}$ from model fitting and used
it for evaluation purpose only. Following previous analyses (Lin et al., 2014;
Shi et al., 2022) on this dataset, we generated bootstrap samples $({\bf
y}^{b},{\bf X}^{b},{\bf Z}^{b})$ of size $n/2$ from the full dataset $({\bf
y},{\bf X},{\bf Z})$, and then used observations in the bootstrap sample for
model training and the other observations not selected in the bootstrap sample
for prediction evaluation. This whole procedure (including generating ${\bf
Z}$ matrix) was repeated for N=100 times and let $({\bf y}^{b},{\bf
X}^{b},{\bf Z}^{b}),b=1,\ldots,N$ denote the $b$th bootstrap sample.

For each bootstrap sample, we fitted the model using $({\bf y}^{b},{\bf
Z}^{b})$ to obtain $\hat{\bm{\beta}}^{b}$. We keep track of observations used
for prediction evaluation to make sure they are not contained in the bootstrap
samples used for model training. Let $C_{-i}$ denote the indices of bootstrap
samples that do not contain observation $i$. That is,

$C_{-i}=\\{b\in 1,\dots,N|x_{i}\notin{\bf X}^{b}\\}.$

Then, the average (over $N$ bootstrap samples) leave-one-out (LOO) squared
prediction error on observation $i$ is calculated as $\sum_{b\in
C_{-i}}(y_{i}-X_{i}\hat{\beta}^{b})^{2}/|C_{-i}|$, where $|C_{-i}|$ denotes
the cardinality of set $C_{-i}$. And the mean squared error over all
observations is given by:

$MSE_{LOO}=\dfrac{1}{n}\sum_{i=1}^{n}\dfrac{1}{|C_{-i}|}\sum_{b\in
C_{-i}}(y_{i}-X_{i}\hat{\beta}^{b})^{2},$

which is used to compare different Lasso estimators. The other metric is the
mean absolute error, which can be analogously defined as:

$MAE_{LOO}=\dfrac{1}{n}\sum_{i=1}^{n}\dfrac{1}{|C_{-i}|}\sum_{b\in
C_{-i}}|y_{i}-X_{i}\hat{\beta}^{b}|.$

Results on prediction errors are reported in Table 4, where one can see that
the proposed Eric Lasso method has the best prediction performance among all
four Lasso methods being considered. Coda Lasso has the second best
performance in terms of leave-one-out mean squared prediction error and CoCo
Lasso has the second best performance in terms of leave-one-out mean absolute
prediction error.

Table 4: Prediction errors on the gut microbiome data by different Lasso estimators. Model | $MSE_{LOO}$ | $MAE_{LOO}$
---|---|---
Eric | 32.838 | 4.316
Coda | 34.928 | 4.468
CoCo | 35.580 | 4.379
Vani | 35.988 | 4.450

Recall that results established in Theorem 1 guarantee performance accuracy in
both estimation and selection. We thus compare selection results of different
Lasso methods. Out of the $N=100$ bootstrap replicates, the frequency of each
genus taxon being selected by each Lasso method was calculated and presented
in Figure 2. As can be seen in the figure, Coda and Vani Lasso tend to select
more taxa than Eric and CoCo Lasso. For example, there are 1, 6, 1, 4 taxa
that are selected by Eric, Coda, CoCo, Vani Lasso, respectively, over 50 times
out of the 100 bootstrap replicates. The distribution of taxa relative
abundances in this COMBO data is very heterogeneous in that the top 5 most
abundant taxa account for 80% of the total abundances of all $p$ taxa. Hence,
the unbalanced compositions generated under Scenario 1 and 3 of Simulation I
better mimic the distribution of taxa relative abundances in this COMBO
dataset than the more homogeneous case used in Scenario 2 of Simulation I. As
observed in Table 2 and Table 3 presented in the previous section, the FPR of
Coda and Vani under Scenario 1 and 3 are 2-5 times of those of Eric and CoCo
Lasso, which may explain why Coda and Vani Lasso have more findings in this
data. In other words, our previous experience in numerical simulations implies
that many of these additional findings of Coda and Vani in this COMBO dataset
might be false positives. The same phenomenon on spurious findings caused by
ignoring measurement errors in microbiome compositional data in statistical
analysis has also been observed in the literature (Hawinkel et al., 2019;
Gihawi et al., 2023). Finally, the common taxon that is selected more than 50%
of times by all methods is genus Acidaminococcus of the Firmicutes phylum,
which was implicated as important for gut dysbiosis in obese patients (Wu et
al., 2011). In summary, the proposed Eric Lasso method can overall provide the
best model prediction and variable selection performance in this gut
microbiome data analysis collected from the COMBO study.

Figure 2: Selection frequencies of each taxon under N=100 bootstrap
replicates.

##  6 Discussion

Motivated by the recent debate on how data contamination could invalidate
scientific findings in microbiome research (Gihawi et al., 2023; Poore et al.,
2020), we have proposed a novel Eric Lasso approach to combat negative effects
of measurement errors in compositional data analysis in order to get more
accurate and reliable statistical analysis results. While existence of
measurement errors or sequence bias in microbiome compositional data has been
widely observed in the literature (McLaren et al., 2019; McLaren et al., 2022;
Zhao and Satten, 2021), statistical methods to address this issue in the high-
dimensional regression framework are less developed. As far as we are
concerned, the only existing method tailored for this specific problem is the
variable correction regularized estimator (Shi et al., 2022), which treats the
observed counts (contaminated copy) as realizations from an underlying
Dirichlet-Multinomial distribution whose parameters are unobserved
compositions (clean copy). A major limitation of the variable correction
regularized estimator method is requirement on knowledge of total counts of
each sample, which however, can be mistakenly measured by many orders of
magnitude in practice (Gihawi et al., 2023). To gain robustness against such a
huge extent in data contamination, we adopt an alternative approach which
directly targets at compositions rather than counts. As shown in our numerical
studies (i.e., Scenario 3 of Simulation I and the case study in the previous
section), our new Eric Lasso approach has a better overall performance than
the variable correction regularized estimator method and the improvement in
reducing false positive findings is substantial when the underlying true
compositions are heterogeneous. While we have only employed the Eric Lasso
method to microbiome compositional data analysis to illustrate its potential
usefulness and superior performance in this article, the proposed methodology
is very general and can be applied to other types of high-dimensional
compositional data as well. Finally, similar to previous publications in this
research vein (Loh and Wainwright, 2012; Datta and Zou, 2017), we assume the
covariance matrix of measurement errors is known or can be estimated from
independent sources such as external data or multiple replicates of the data.
In the future work, we intend to explore different backup strategies
(e.g.,resampling techniques) when such information is less available.

## References

  * Aitchison (1982) Aitchison, J. (1982).  The statistical analysis of compositional data.  Journal of the Royal Statistical Society: Series B (Methodological) 44(2), 139–160.
  * Aitchison and Bacon-Shone (1984) Aitchison, J. and J. Bacon-Shone (1984).  Log contrast models for experiments with mixtures.  Biometrika 71(2), 323–330.
  * Aitchison and Shen (1980) Aitchison, J. and S. M. Shen (1980).  Logistic-normal distributions: Some properties and uses.  Biometrika 67(2), 261–272.
  * Allali et al. (2017) Allali, I., J. W. Arnold, J. Roach, et al. (2017).  A comparison of sequencing platforms and bioinformatics pipelines for compositional analysis of the gut microbiome.  BMC microbiology 17(1), 1–16.
  * Belloni et al. (2017) Belloni, A., M. Rosenbaum, and A. B. Tsybakov (2017).  Linear and conic programming estimators in high dimensional errors-in-variables models.  Journal of the Royal Statistical Society Series B: Statistical Methodology 79(3), 939–956.
  * Billheimer et al. (2001) Billheimer, D., P. Guttorp, and W. F. Fagan (2001).  Statistical interpretation of species composition.  Journal of the American statistical Association 96(456), 1205–1214.
  * Cai et al. (2022) Cai, B., J. Zhang, H. Li, C. Su, and H. Zhao (2022).  Statistical inference of cell-type proportions estimated from bulk expression data.  arXiv preprint arXiv:2209.04038.
  * Clausen and Willis (2022) Clausen, D. S. and A. D. Willis (2022).  Modeling complex measurement error in microbiome experiments.  arXiv preprint arXiv:2204.12733.
  * Combettes and Müller (2021) Combettes, P. L. and C. L. Müller (2021).  Regression models for compositional data: General log-contrast formulations, proximal optimization, and microbiome data applications.  Statistics in Biosciences 13(2), 217–242.
  * Datta and Zou (2017) Datta, A. and H. Zou (2017).  Cocolasso for high-dimensional error-in-variables regression.  Annals of Statistics 45(6), 2400–2426.
  * Datta and Zou (2020) Datta, A. and H. Zou (2020).  A note on cross-validation for lasso under measurement errors.  Technometrics 62(4), 549–556.
  * Fiksel et al. (2022) Fiksel, J., S. Zeger, and A. Datta (2022).  A transformation-free linear regression for compositional outcomes and predictors.  Biometrics 78(3), 974–987.
  * Firth and Sammut (2023) Firth, D. and F. Sammut (2023).  Analysis of composition on the original scale of measurement.  arXiv preprint arXiv:2312.10548.
  * Gihawi et al. (2023) Gihawi, A., Y. Ge, J. Lu, D. Puiu, A. Xu, C. S. Cooper, D. S. Brewer, M. Pertea, and S. L. Salzberg (2023).  Major data analysis errors invalidate cancer microbiome findings.  Mbio, e01607–23.
  * Greenacre et al. (2023) Greenacre, M., E. Grunsky, J. Bacon-Shone, I. Erb, and T. Quinn (2023).  Aitchison’s compositional data analysis 40 years on: A reappraisal.  Statistical Science 1(1), 1–25.
  * Hawinkel et al. (2019) Hawinkel, S., F. Mattiello, L. Bijnens, and O. Thas (2019).  A broken promise: microbiome differential abundance methods do not control the false discovery rate.  Briefings in bioinformatics 20(1), 210–221.
  * Hron et al. (2012) Hron, K., P. Filzmoser, and K. Thompson (2012).  Linear regression with compositional explanatory variables.  Journal of applied statistics 39(5), 1115–1128.
  * Jiang et al. (2023) Jiang, R., X. Zhan, and T. Wang (2023).  A flexible zero-inflated poisson-gamma model with application to microbiome sequence count data.  Journal of the American Statistical Association 118(542), 792–804.
  * Lin et al. (2014) Lin, W., P. Shi, R. Feng, and H. Li (2014).  Variable selection in regression with compositional covariates.  Biometrika 101(4), 785–797.
  * Loh and Wainwright (2012) Loh, P.-L. and M. J. Wainwright (2012).  High-dimensional regression with noisy and missing data: Provable guarantees with non-convexity.  Annals of Statistics 40(3), 1637–1664.
  * Martin et al. (2020) Martin, B. D., D. Witten, and A. D. Willis (2020).  Modeling microbial abundances and dysbiosis with beta-binomial regression.  Annals of Applied Statistics 14(1), 94.
  * McLaren et al. (2022) McLaren, M. R., J. T. Nearing, A. D. Willis, K. G. Lloyd, and B. J. Callahan (2022).  Implications of taxonomic bias for microbial differential-abundance analysis.  bioRxiv, 2022–08.
  * McLaren et al. (2019) McLaren, M. R., A. D. Willis, and B. J. Callahan (2019).  Consistent and correctable bias in metagenomic sequencing experiments.  Elife 8, e46923.
  * Mishra and Müller (2022) Mishra, A. and C. L. Müller (2022).  Robust regression with compositional covariates.  Computational Statistics & Data Analysis 165, 107315.
  * Poore et al. (2020) Poore, G. D., E. Kopylova, Q. Zhu, C. Carpenter, S. Fraraccio, S. Wandro, T. Kosciolek, S. Janssen, J. Metcalf, S. J. Song, et al. (2020).  Microbiome analyses of blood and tissues suggest cancer diagnostic approach.  Nature 579(7800), 567–574.
  * Randolph et al. (2018) Randolph, T. W., S. Zhao, W. Copeland, M. Hullar, and A. Shojaie (2018).  Kernel-penalized regression for analysis of microbiome data.  The annals of applied statistics 12(1), 540.
  * Rosenbaum and Tsybakov (2010) Rosenbaum, M. and A. B. Tsybakov (2010).  Sparse recovery under matrix uncertainty.  Annals of Statistics 38(5), 2620–2651.
  * Shi et al. (2016) Shi, P., A. Zhang, H. Li, et al. (2016).  Regression analysis for microbiome compositional data.  Annals of Applied Statistics 10(2), 1019–1040.
  * Shi et al. (2022) Shi, P., Y. Zhou, and A. R. Zhang (2022).  High-dimensional log-error-in-variable regression with applications to microbial compositional data analysis.  Biometrika 109(2), 405–420.
  * Srinivasan et al. (2021) Srinivasan, A., L. Xue, and X. Zhan (2021).  Compositional knockoff filter for high-dimensional regression analysis of microbiome data.  Biometrics 77(3), 984–995.
  * Susin et al. (2020) Susin, A., Y. Wang, K.-A. Lê Cao, and M. L. Calle (2020).  Variable selection in microbiome compositional data analysis.  NAR Genomics and Bioinformatics 2(2), lqaa029.
  * Tang et al. (2017) Tang, Z.-Z., G. Chen, A. V. Alekseyenko, and H. Li (2017).  A general framework for association analysis of microbial communities on a taxonomic tree.  Bioinformatics 33(9), 1278–1285.
  * Tibshirani (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society Series B: Statistical Methodology 58(1), 267–288.
  * Vandeputte et al. (2017) Vandeputte, D., G. Kathagen, K. D’hoe, et al. (2017).  Quantitative microbiome profiling links gut community variation to microbial load.  Nature 551(7681), 507–511.
  * Wang and Zhao (2017) Wang, T. and H. Zhao (2017).  Structured subcomposition selection in regression and its application to microbiome data analysis.  Annals of Applied Statistics 11(2), 771–791.
  * Wu et al. (2011) Wu, G. D., J. Chen, C. Hoffmann, et al. (2011).  Linking long-term dietary patterns with gut microbial enterotypes.  Science 334(6052), 105–108.
  * Zhao and Satten (2021) Zhao, N. and G. A. Satten (2021).  A log-linear model for inference on bias in microbiome studies.  Statistical Analysis of Microbiome Data, 221–246.

Supplementary Materials for “High-dimensional log contrast models with
measurement errors”

## A Proofs

For ease of presentation, we first introduce the following notations:

$\begin{array}[]{cc}\tilde{A}=D^{p}_{SS}(\widetilde{\Sigma}^{p}_{SS})^{-1}(D^{p}_{SS})^{T}&A=D^{p}_{SS}(\Sigma^{p}_{SS})^{-1}(D^{p}_{SS})^{T}\\\
E=\widetilde{\Sigma}^{p}-\Sigma^{p}&F_{SS}=(\widetilde{\Sigma}^{p}_{SS})^{-1}-(\Sigma^{p}_{SS})^{-1}\\\
G=\Sigma^{p}_{S^{c}S}(\Sigma^{p}_{SS})^{-1}&\tilde{G}=\widetilde{\Sigma}^{p}_{S^{c}S}(\widetilde{\Sigma}^{p}_{SS})^{-1}\\\
H=\tilde{G}-G&b=\|\beta_{S}^{*}\|_{\infty}\\\
\psi=\|\Sigma_{SS}\|_{\infty}&\nu=\|\Sigma_{S}\|_{\infty}\\\ \end{array}\\\ $
$\Pi=I-\frac{1}{n}X^{p}_{S}(\Sigma^{p}_{SS})^{-1}(X^{p}_{S})^{T}$

We adopt the identical closeness condition as used in CoCoLasso (Datta and
Zou, 2017).

######  Condition 3 (Closeness condition).

Let us assume that the distribution of $\hat{\Sigma}$ and $\tilde{\rho}$ are
identified by a set of parameters $\theta$. Then there exists universal
constants $C$ and $c$, and positive functions $\zeta$ and $\epsilon_{0}$
depending on $\beta^{*}_{S}$, $\theta$ and $\sigma^{2}$ such that for every
$\epsilon\leq\epsilon_{0}$, $\hat{\Sigma}$ and $\tilde{\rho}$ satisfy the
following probability statements:

$\begin{array}[]{c}\rm pr(|\widehat{\Sigma}_{ij}-\Sigma_{ij}|\geq\epsilon)\leq
C\exp\left(-cn\epsilon^{2}\zeta^{-1}\right)\;\;i,j=1,\ldots,p,\\\ \rm
pr(|\tilde{\rho}_{j}-\rho_{j}|\geq\epsilon)\leq
C\exp\left(-cns^{-2}\epsilon^{2}\zeta^{-1}\right)\;\;j=1,\ldots,p.\end{array}$
(S.1)

As indicated in CoCoLasso, the closeness condition holds for
$\widehat{\Sigma}$ with $\zeta=\max(\tau^{4},\sigma^{4},1)$ and
$\epsilon_{0}=\tau^{2}$. In scenarios where $\Sigma_{B}$ is unknown and is
substituted by $\hat{\Sigma}_{B}=\frac{1}{n}B_{o}^{T}B_{o}$, this condition
still satisfies with the same parameter values as guaranteed by the following
Lemma 1:

######  Lemma 1.

$\widehat{\Sigma}_{o}=\frac{1}{n}Z^{T}Z-\widehat{\Sigma}_{B}$ and
$\tilde{\rho}$ satisfy the closeness condition (S.1) with
$\zeta=\max(\tau^{4},\sigma^{4},1)$ and $\epsilon_{0}=\tau^{2}$.

###### Proof.

Since $\widehat{\Sigma}$ satisfies the closeness condition, we have

$\rm pr(|\widehat{\Sigma}_{ij}-\Sigma_{ij}|\geq\epsilon)\leq
C\exp\left(-cn\epsilon^{2}\zeta^{-1}\right).$

Using Lemma B.1 of Datta and Zou (2017),

$\rm pr(|(\widehat{\Sigma}_{B}-\Sigma_{B})_{ij}|\geq\epsilon)\leq
C\exp\left(-cn\epsilon^{2}\tau^{-4}\right)\leq
C\exp\left(-cn\epsilon^{2}\zeta^{-1}\right),$
$\widehat{\Sigma}_{o}-\Sigma=\widehat{\Sigma}_{o}-\widehat{\Sigma}+(\widehat{\Sigma}-\Sigma)=\Sigma_{B}-\widehat{\Sigma}_{B}+(\widehat{\Sigma}-\Sigma).$

Hence by redefining $C$ and $c$, the sub-Gaussianity of $B_{o}$ implies

$\rm pr(|(\widehat{\Sigma}_{o})_{ij}-\Sigma_{ij}|\geq\epsilon)\leq
C\exp\left(-cn\epsilon^{2}\zeta^{-1}\right),$

which completes the proof. ∎

We will now proceed to prove Theorem 1 of the main text under the assumption
that the closeness condition is satisfied. To establish major results
presented in Theorem 1, we begin by proving the following lemmas.

######  Lemma 2.

For any $\epsilon>0$ we have

$\rm pr(\|\widetilde{\Sigma}^{p}-\Sigma^{p}\|_{\max}\geq\epsilon)\leq
p^{2}\max_{k,l}\rm pr(|\hat{\Sigma}_{kl}-\Sigma_{kl}|\geq\epsilon/8).$ (S.2)

###### Proof.

Plugging

$D^{p}=\begin{pmatrix}I_{p-1}\\\ -1_{p-1}\end{pmatrix}$

into $\tilde{\Sigma}^{p}=(D^{p})^{T}\tilde{\Sigma}D^{p}$, we have

$(\widetilde{\Sigma}^{p})_{kl}=\tilde{\Sigma}_{kl}-\tilde{\Sigma}_{kp}-\tilde{\Sigma}_{pl}+\tilde{\Sigma}_{pp}.$

As $\Sigma^{p}=D_{p}^{T}\Sigma D_{p}$, this further implies that

$\|\widetilde{\Sigma}^{p}-\Sigma^{p}\|_{\max}\leq
4\|\tilde{\Sigma}-\Sigma\|_{\max}.$

By definition of $\tilde{\Sigma}$ and the positive semi-definite nature of
$\Sigma$, we have

$\|\tilde{\Sigma}-\Sigma\|_{\max}\leq\|\tilde{\Sigma}-\hat{\Sigma}\|_{\max}+\|\hat{\Sigma}-\Sigma\|_{\max}\leq
2\|\hat{\Sigma}-\Sigma\|_{\max},$ $\rm
pr(\|\tilde{\Sigma}-\Sigma\|_{\max}\geq\epsilon/4)\leq\rm
pr(\|\hat{\Sigma}-\Sigma\|_{\max}\geq\epsilon/8).$

The proof then follows using union bounds over $\rm
pr(|\hat{\Sigma}_{kl}-\Sigma_{kl}|\geq\epsilon/8)$. ∎

The closeness condition requires that $\hat{\Sigma}$ and $\tilde{\rho}$ are
sufficiently close to $\Sigma$ and $\rho$, respectively. A direct implication
of Lemma 2 is

$\|\tilde{\Sigma}^{p}-\Sigma^{p}\|_{\max}\leq
4\|\tilde{\Sigma}-\Sigma\|_{\max}\leq 8\|\hat{\Sigma}-\Sigma\|_{\max},$

which implies that $\tilde{\Sigma}^{p}$ is sufficiently close to $\Sigma^{p}$.
Analogously, since
($\tilde{\rho}^{p}-\rho^{p})=(D^{p})^{T}(\tilde{\rho}-\rho)$, the closeness
condition of $\tilde{\rho}$ on $\rho$ ensures that $\tilde{\rho}^{p}$
approximates $\rho^{p}$ well.

######  Lemma 3.

$\Lambda_{\min}(\Sigma_{SS})\leq\Lambda_{\min}(\Sigma^{p}_{SS})$.

###### Proof.

It can be shown that $\|D_{SS}^{p}v\|^{2}_{2}\geq\|v\|^{2}_{2}$ for any
$v\in\mathbb{R}^{s-1}\backslash\\{0\\}$. Hence we have

$\displaystyle\Lambda_{\min}(\Sigma_{SS})$
$\displaystyle=\frac{1}{n}\min_{v\in\mathbb{R}^{s}\backslash\\{0\\}}\dfrac{\|X_{S}v\|^{2}_{2}}{\|v\|^{2}_{2}}$
$\displaystyle\leq\frac{1}{n}\min_{v\in\mathbb{R}^{s-1}\backslash\\{0\\}}\dfrac{\|X_{S}D_{SS}^{p}v\|^{2}_{2}}{\|D_{SS}^{p}v\|^{2}_{2}}$
$\displaystyle\leq\frac{1}{n}\min_{v\in\mathbb{R}^{s-1}\backslash\\{0\\}}\dfrac{\|X_{S}^{p}v\|^{2}_{2}}{\|v\|^{2}_{2}}$
$\displaystyle=\Lambda_{\min}(\Sigma^{p}_{SS}).$

∎

######  Lemma 4.

$\rm pr(\widetilde{\Sigma}^{p}_{SS}>0)\geq
1-Cp^{2}\exp\left(-cn(s-1)^{-2}\epsilon^{2}\zeta^{-1}\right)$ for all
$\epsilon\leq min(\epsilon_{0},C_{\min}/16)$.

###### Proof.

Using Lemma 3,

$\displaystyle\Lambda_{\min}(\widetilde{\Sigma}^{p}_{SS})\geq$
$\displaystyle\Lambda_{\min}(\Sigma^{p}_{SS})-|\Lambda_{\max}(-E_{SS})|\geq
C_{\min}-\|E_{SS}\|_{2}$ $\displaystyle\geq$ $\displaystyle
C_{\min}-(s-1)\|E_{SS}\|_{\max}\geq C_{\min}-(s-1)\|E\|_{\max}\geq
C_{\min}/2,$

where the last inequality occurs when $\|E\|_{\max}\leq 8(s-1)^{-1}\epsilon$
with probability at least
$1-Cp^{2}\exp\left(-cn(s-1)^{-2}\epsilon^{2}\zeta^{-1}\right)$ for
$\epsilon\leq min(\epsilon_{0},C_{\min}/16)$, according to the closeness
condition and Lemma 2. ∎

######  Lemma 5.

$\|\frac{1}{n}(X_{S})^{T}w\|_{\infty}\leq\lambda/2$ (S.3)

holds with probability at least
$1-s\exp\left(-n\lambda^{2}/(8\sigma^{2})\right)$, and

$\frac{1}{n}\|(X^{p}_{S^{c}})^{T}\Pi w\|_{\infty}\leq\lambda\xi/4$ (S.4)

holds with probability at least
$1-(p-s)\exp\left(-n\lambda^{2}\xi^{2}/(128\sigma^{2})\right)$.

###### Proof.

This lemma is taken from the proof of Lin et al. (2014) and these inequalities
hold as a result of the Gaussian tail bound. ∎

### Proof of Theorem 1

For the optimality, we aim to prove (S.5) and (S.6).

$\displaystyle(\widetilde{\Sigma}^{p}_{\hat{S}\hat{S}},\widetilde{\Sigma}^{p}_{\hat{S}\hat{S}^{c}})\hat{\beta}_{-p}-\widetilde{\rho}^{p}_{\hat{S}}+\lambda\\{\mbox{sgn}(\hat{\beta}_{\hat{S}-p})-\mbox{sgn}(\hat{\beta}_{p})1_{s-1}\\}$
$\displaystyle=0,$ (S.5)
$\displaystyle\|(\widetilde{\Sigma}^{p}_{\hat{S}^{c}\hat{S}},\widetilde{\Sigma}^{p}_{\hat{S}^{c}\hat{S}^{c}})\hat{\beta}_{-p}-\widetilde{\rho}^{p}_{\hat{S}^{c}}-\lambda\mbox{sgn}(\hat{\beta}_{p})1_{p-s}\|_{\infty}$
$\displaystyle\leq\lambda.$ (S.6)

### Part I: proof of (S.5)

Let $u$ denote
$\\{\mbox{sgn}(\widehat{\beta}_{S-p})-\mbox{sgn}(\hat{\beta}_{p})1_{s-1}\\}=(D_{SS}^{p})^{T}\mbox{sgn}(\hat{\beta}_{S})$
for simplicity. Conditioned on Lemma 4, we establish $\hat{\beta}$ such that

$\displaystyle\hat{\beta}_{S-p}$
$\displaystyle=(\widetilde{\Sigma}^{p}_{SS})^{-1}(\widetilde{\rho}^{p}_{S}-\lambda
u),$ (S.7) $\displaystyle\hat{\beta}_{S^{c}}$ $\displaystyle=0,$ (S.8)

where (S.7) is solved from (S.5) by replacing $S$ with $\hat{S}$. We rewrite
it as

$\displaystyle\hat{\beta}_{S}-\beta^{*}_{S}=$ $\displaystyle
D_{SS}^{p}(\widetilde{\Sigma}^{p}_{SS})^{-1}\left(\widetilde{\rho}^{p}_{S}-\rho^{p}_{S}+\Sigma^{p}_{SS}\beta_{S-p}^{*}+\frac{1}{n}(X^{p}_{S})^{T}w-\lambda
u\right)-D_{SS}^{p}\beta_{S-p}^{*}$ $\displaystyle=$ $\displaystyle
D_{SS}^{p}F_{SS}(D_{SS}^{p})^{T}\left(\tilde{\rho}_{S}-\rho_{S}+\Sigma_{SS}\beta_{S}^{*}+\frac{1}{n}(X_{S})^{T}w\right)+A(\tilde{\rho}_{S}-\rho_{S})+\frac{1}{n}A(X_{S})^{T}w-\lambda\tilde{A}\mbox{sgn}(\hat{\beta}_{S}).$

Let
$\eta_{1}=\|D_{SS}^{p}F_{SS}(D_{SS}^{p})^{T}\|_{\infty},\eta_{2}=\|D_{SS}^{p}E_{SS}(D_{SS}^{p})^{T}\|_{\infty}$.
As $(D_{SS}^{p})^{T}D_{SS}^{p}=I_{s-1}$, we have

$\eta_{1}=\|D_{SS}^{p}(\widetilde{\Sigma}^{p}_{SS})^{-1}(D_{SS}^{p})^{T}D_{SS}^{p}E_{SS}(D_{SS}^{p})^{T}A\|_{\infty}\leq(\eta_{1}+\phi)\eta_{2}\phi.$

Therefore, by the Lemma 2 and the closeness condition, for
$\epsilon\leq\min(\epsilon_{0},(2\phi)^{-1})$,

$\eta_{2}\leq 2(s-1)^{2}\|E\|_{\max}\leq\epsilon$ (S.9)

with probability at least
$1-Cp^{2}\exp\left(-\frac{cn\epsilon^{2}}{256(s-1)^{4}\zeta}\right)$. Then

$\eta_{1}\leq\dfrac{\phi^{2}\eta_{2}}{1-\eta_{2}\phi}\leq\dfrac{\phi^{2}\epsilon}{1-\epsilon\phi}\leq\phi,$
(S.10) $\|\tilde{A}\|_{\infty}\leq\phi+\eta_{1}\leq 2\phi.$ (S.11)

The closeness condition and Lemma 5 imply that for $\lambda\leq\epsilon_{0}$

$\|\tilde{\rho}_{S}-\rho_{S}\|_{\infty}\leq\lambda,$ (S.12)
$\|\frac{1}{n}(X_{S})^{T}w\|_{\infty}\leq\lambda/2,$ (S.13)

with probability at least
$1-Cs\exp\left(-cns^{-2}\lambda^{2}\zeta^{-1}\right)-s\exp\left(-n\lambda^{2}\sigma^{-2}/8\right)$.
Combining all these inequalities, for
$\epsilon\leq\min\\{\lambda\phi^{-1}(3\lambda+b\psi)^{-1},\epsilon_{0},\phi^{-1}/2\\}$,
we have

$\displaystyle\|\hat{\beta}_{S}-\beta^{*}_{S}\|_{\infty}\leq$
$\displaystyle\dfrac{\phi^{2}\epsilon}{1-\epsilon\phi}(\lambda+\lambda/2+b\psi)+\phi\|\tilde{\rho}_{S}-\rho_{S}\|_{\infty}+\phi\|\frac{1}{n}(X_{S})^{T}w\|_{\infty}+2\phi\lambda$
$\displaystyle\leq$ $\displaystyle 9\lambda\phi/2.$

Then we can conclude that
$\mbox{sgn}(\hat{\beta}_{S})=\mbox{sgn}(\beta^{*}_{S})$ and (S.5) hold.

### Part II: proof of (S.6)

Using (S.7) and substituting $u$ with
$\mbox{sgn}(\widehat{\beta}_{S-p})-\mbox{sgn}(\hat{\beta}_{p})1_{s-1}$, we
have

$(\widetilde{\Sigma}^{p}_{S^{c}S},\widetilde{\Sigma}^{p}_{S^{c}S^{c}})\hat{\beta}_{-p}-\widetilde{\rho}^{p}_{S^{c}}-\lambda\mbox{sgn}(\hat{\beta}_{p})1_{p-s}=\tilde{G}(\widetilde{\rho}^{p}_{S}-\lambda
u)-\widetilde{\rho}^{p}_{S^{c}}-\lambda\mbox{sgn}(\hat{\beta}_{p})1_{p-s}.$

Taking the absolute values and using triangular inequalities, we have

$\|\tilde{G}(\widetilde{\rho}^{p}_{S}-\lambda
u)-\widetilde{\rho}^{p}_{S^{c}}-\lambda\mbox{sgn}(\hat{\beta}_{p})1_{p-s}\|_{\infty}\\\
\leq\|\tilde{G}\widetilde{\rho}^{p}_{S}-\lambda
Hu-\widetilde{\rho}^{p}_{S^{c}}\|_{\infty}+\lambda\|Gu+\mbox{sgn}(\hat{\beta}_{p})1_{p-s}\|_{\infty}.$

Condition 2 of the main text implies
$\|Gu+\mbox{sgn}(\hat{\beta}_{p})1_{p-s}\|_{\infty}\leq 1-\xi$. We break the
first term as

$\lambda
Hu+\widetilde{\rho}^{p}_{S^{c}}-\tilde{G}\widetilde{\rho}^{p}_{S}=H(\lambda
u-\widetilde{\rho}^{p}_{S})+\\{(\widetilde{\rho}^{p}_{S^{c}}-\rho^{p}_{S^{c}})-G(\widetilde{\rho}^{p}_{S}-\rho^{p}_{S})\\}+(\rho^{p}_{S^{c}}-G\rho^{p}_{S}),$

where
$\|\rho^{p}_{S^{c}}-G\rho^{p}_{S}\|_{\infty}=\frac{1}{n}\|(X^{p}_{S^{c}})^{T}\Pi
w\|_{\infty}\leq\lambda\xi/4$ with probability at least
$1-(p-s)\exp\left(-n\lambda^{2}\xi^{2}/(128\sigma^{2})\right)$, according to
Lemma 5. We further bound $G(D_{SS}^{p})^{T}$ as

$\displaystyle\|G(D_{SS}^{p})^{T}\|_{\infty}$
$\displaystyle=\|\frac{1}{n}(X^{p}_{S^{c}})^{T}X_{S}A\|_{\infty}=\max_{j\in
S^{c}}\|\frac{1}{n}(X_{j}-X_{p})^{T}X_{S}A\|_{1}$ (S.14) $\displaystyle\leq
2\max_{j\in\\{S^{c},p\\}}\|\frac{1}{n}X_{j}^{T}X_{S}A\|_{1}=2\|\Sigma_{\\{S^{c},p\\},S}A\|_{\infty}$
$\displaystyle\leq 2\nu\phi.$

We assume that $\epsilon\leq\frac{\lambda\xi}{4+4\nu\phi}$, hence

$\displaystyle\|(\widetilde{\rho}^{p}_{S^{c}}-\rho^{p}_{S^{c}})-G(\widetilde{\rho}^{p}_{S}-\rho^{p}_{S})\|_{\infty}$
$\displaystyle=\|(\widetilde{\rho}^{p}_{S^{c}}-\rho^{p}_{S^{c}})-G(D_{SS}^{p})^{T}(\tilde{\rho}_{S}-\rho_{S})\|_{\infty}$
$\displaystyle\leq(2+2\nu\phi)\|\tilde{\rho}-\rho\|_{\infty}\leq\lambda\xi/2$

holds with probability at least
$1-Cp\exp\\{-cns^{-2}\epsilon^{2}\zeta^{-1}\\}$. Since the term $H(\lambda
u-\widetilde{\rho}^{p}_{S})=H(D_{SS}^{p})^{T}\left(\lambda\mbox{sgn}(\hat{\beta}_{S})-\tilde{\rho}_{S}\right)$,
we have:

$\displaystyle\|H(D_{SS}^{p})^{T}\|_{\infty}$
$\displaystyle=\|\Sigma^{p}_{S^{c}S}F_{SS}(D_{SS}^{p})^{T}+E_{S^{c}S}(\widetilde{\Sigma}^{p}_{SS})^{-1}(D_{SS}^{p})^{T}\|_{\infty}$
$\displaystyle\leq\|\Sigma^{p}_{S^{c}S}F_{SS}(D_{SS}^{p})^{T}\|_{\infty}+\|E_{S^{c}S}(\widetilde{\Sigma}^{p}_{SS})^{-1}(D_{SS}^{p})^{T}\|_{\infty}$
$\displaystyle=\|\frac{1}{n}(X^{p}_{S^{c}})^{T}X_{S}(A-\tilde{A})\|_{\infty}+\|E_{S^{c}S}(D_{SS}^{p})^{T}\tilde{A}\|_{\infty}$
$\displaystyle\leq\|\frac{1}{n}(X^{p}_{S^{c}})^{T}X_{S}(A-\tilde{A})\|_{\infty}+2\phi\|E_{S^{c}S}(D_{SS}^{p})^{T}\|_{\infty}$
$\displaystyle\leq 2\nu\eta_{1}+4(s-1)\phi\|E\|_{\max},$

where the last inequality follows from the analogous trick in (S.14). As
mentioned in (S.13),we have
$\|\rho_{S}\|_{\infty}\leq\|\frac{1}{n}X_{S}^{T}w\|_{\infty}+\|\Sigma_{SS}\beta^{*}_{S}\|_{\infty}\leq\lambda/2+b\psi$.
Along with (S.12), we have

$\|\lambda\mbox{sgn}(\hat{\beta}_{S})-\tilde{\rho}_{S}\|_{\infty}\leq\|\rho_{S}\|_{\infty}+\|\tilde{\rho}_{S}-\rho_{S}\|_{\infty}+\lambda\leq
5\lambda/2+b\psi.$

Note that (S.9) and (S.10) yield that $\eta_{1}\leq 2\phi^{2}\eta_{2}\leq
4(s-1)^{2}\phi^{2}\|E\|_{\max}\leq 2\phi^{2}\epsilon$. For
$\epsilon\leq\frac{\lambda\xi}{(10\lambda+4\psi)(4\nu\phi^{2}+2\phi)}$,

$\displaystyle\|H(\lambda u-\widetilde{\rho}^{p}_{S})\|_{\infty}$
$\displaystyle\leq(5\lambda/2+b\psi)(2\nu\eta_{1}+4(s-1)\phi\|E\|_{\max})$
$\displaystyle\leq(5\lambda/2+b\psi)[4\nu\phi^{2}\epsilon+2\phi\epsilon/(s-1)]$
$\displaystyle\leq(5\lambda/2+b\psi)(4\nu\phi^{2}+2\phi)\epsilon$
$\displaystyle\leq\lambda\xi/4.$

Gathering them all together, we conclude that for

$\epsilon\leq\min\\{\frac{1}{2\phi},\epsilon_{0},\frac{\lambda}{\phi(3\lambda+b\psi)},\frac{\lambda\xi}{(10\lambda+4\psi)(4\nu\phi^{2}+2\phi)},\frac{\lambda\xi}{4+4\nu\phi},\frac{C_{\min}}{16}\\}$

and $\lambda\leq\epsilon_{0}$, equations (S.5) and (S.6) hold with probability
at least
$1-C_{1}p^{2}\exp\left(-c_{1}n(s-1)^{-4}\epsilon^{2}\zeta^{-1}\right)-C_{2}p\exp\left(-c_{2}ns^{-2}\lambda^{2}\xi^{2}\zeta^{-1}\right).$

## B Additional simulation results

In Table 1 of the main text, we have examined the sum of regression
coefficients $\sum_{j=1}^{p}\hat{\beta}_{j}$ estimated by different Lasso
methods under Scenario 1 of Simulation I. The corresponding results under
Scenario 2 and Scenario 3 of Simulation I are reported in Table S.1 and Table
S.2, respectively. As observed in Table 2 of the main text, the sum of
regression coefficients estimated by Eric Lasso and Coda Lasso is very close
to zero, and the departure from zero for CoCo Lasso or Vani Lasso is more
substantial. It is of interest to observe that test results of CoCo Lasso and
Vani Lasso may not be significant under certain scenarios, especially under
Scenario 2 (i.e., in Table S.1). We take a closer look at the sum of
regression coefficients ($\sum_{j=1}^{p}\hat{\beta}_{j}$) estimated by CoCo
Lasso and Vanilla Lasso and presented histograms of the 100 replicates under
each scenario in Figure S.1. It can be seen that the spread of each empirical
distribution is relatively large compared to its mean, leading to an
insignificant p-value. Like what have been observed in the main text, the
obvious departure of sum of regression coefficients from zero for CoCoLasso
and Vani Lasso makes it difficult to interpret the estimated regression model
coefficients under the framework of compositional data analysis as explained
in Section 2 of the main text.

Table S.1: The average sum of regression coefficients estimated by different Lasso methods and corresponding t-test p values under Scenario 2. (n,p) | Model | $\sum\beta_{j}$ | p value
---|---|---|---
(100,200) | Eric | -2.2e-08 | 0.924
| Coda | 3.5e-08 | 0.847
| CoCo | 1.9e-02 | 0.305
| Vani | 1.8e-02 | 0.343
(250,400) | Eric | -5.0e-07 | 0.269
| Coda | 6.8e-08 | 0.846
| CoCo | -8.0e-03 | 0.650
| Vani | -1.2e-02 | 0.508
(500,500) | Eric | 2.9e-07 | 0.532
| Coda | 1.9e-07 | 0.530
| CoCo | 4.0e-03 | 0.695
| Vani | 5.0e-03 | 0.631
Table S.2: The average sum of regression coefficients estimated by different Lasso methods and corresponding t-test p values under Scenario 3. (n,p) | Model | $\sum\beta_{j}$ | p value
---|---|---|---
(100,200) | Eric | 6.4e-08 | 0.047
| Coda | 6.2e-08 | 0.047
| CoCo | -6.5e-01 | 4.7e-24
| Vani | -1.5e-01 | 0.011
(250,400) | Eric | 1.4e-08 | 0.793
| Coda | -3.0e-08 | 0.138
| CoCo | -4.8e-02 | 9.9e-33
| Vani | 5.0e-02 | 0.238
(500,500) | Eric | -4.1e-08 | 0.549
| Coda | -1.7e-09 | 0.930
| CoCo | -3.8e-02 | 5.0e-29
| Vani | 8.6e-2 | 0.006
Figure S.1: The histograms of the sum of coefficients under Scenario 2 for
model CoCoLasso and Vanilla Lasso.

We next report ROC curves under Scenario 2 of Simulation II in Figure S.2.
Similar to patterns displayed in Figure 1 of the main text, Eric Lasso tends
to have the best performance among the four methods under Scenario 2 of
Simulation II.

|
---|---
Figure S.2: ROC curves with different $\tau$ values under Scenario 2 of
Simulation II. The left panel is ROC curves at original FPR and TPR scales and
the right panel zooms in specific regions of FPR and TPR to better distinguish
different methods.