# High-dimensional Log-Error-in-Variable Regression with Applications to
Microbial Compositional Data Analysis111Pixu Shi is Visiting Assistant
Professor (E-mail: pixushi@stat.wisc.edu), Yuchen Zhou is PhD student (E-mail:
yuchenzhou@stat.wisc.edu), and Anru Zhang is Assistant Professor (E-mail:
anruzhang@stat.wisc.edu), Department of Statistics, University of Wisconsin-
Madison, WI 53706.

Pixu Shi1, Yuchen Zhou1, and Anru Zhang

(())

###### Abstract

In microbiome and genomic study, the regression of compositional data has been
a crucial tool for identifying microbial taxa or genes that are associated
with clinical phenotypes. To account for the variation in sequencing depth,
the classic log-contrast model is often used where read counts are normalized
into compositions. However, zero read counts and the randomness in covariates
remain critical issues.

In this article, we introduce a surprisingly simple, interpretable, and
efficient method for the estimation of compositional data regression through
the lens of a novel high-dimensional log-error-in-variable regression model.
The proposed method provides both corrections on sequencing data with possible
overdispersion and simultaneously avoids any subjective imputation of zero
read counts. We provide theoretical justifications with matching upper and
lower bounds for the estimation error. We also consider a general log-error-
in-variable regression model with corresponding estimation method to
accommodate broader situations. The merit of the procedure is illustrated
through real data analysis and simulation studies.

11footnotetext: The first two authors contributed equally and should be
considered as co-first authors.

Keywords: compositional data, error-in-variable, high-dimensional regression,
microbiome study, minimax optimality

##  1 Introduction

Motivated by an array of applications, high-dimensional regression has
attracted enormous attention in contemporary research. The canonical model of
high-dimensional regression can be written as
$y={\mathbf{X}}\beta^{\ast}+\varepsilon$, where
$y=(y_{1},\ldots,y_{n})^{\top}$ is the response vector,
${\mathbf{X}}\in\mathbb{R}^{n\times p}$ is the covariate matrix, and
$\beta^{\ast}\in\mathbb{R}^{p}$ is the unknown coefficient vector of interest.
When sample size $n$ is much smaller than the number of available predictor
$p$, but larger than the number of relevant parameters $s$, one can perform
various well-established methods, such as Lasso (Tibshirani, 1996), SCAD (Fan
and Li, 2001), and Dantzig selector (Candes and Tao, 2005) among others, to
estimate $\beta^{\ast}$.

Much prior attention to high-dimensional regression focused on the “clean
data” case where the covariates are accurately observed. However, in
applications of econometrics (Greene, 2003), genomics (Purdom and Holmes,
2005; Cao et al., 2017), and engineering (Slijepcevic et al., 2002), we also
frequently see covariates corrupted with noise. Previous literature referred
to such scenarios as “error-in-variable” and showed that performing standard
regression methods directly on the corrupted covariates may yield inaccurate
inference results (Hausman, 2001). When the observable covariates are
corrupted by additive Gaussian or sub-Gaussian noises, the methods and
theories for error-in-variable regression have been widely considered
previously. For example, to adjust for the bias of regular least square
estimator in the low-dimensional setting, various methods, e.g., Deming
regression (Deming, 1943) and method of moments (Pal, 1980), were introduced.
In more recent high-dimensional settings, Rosenbaum and Tsybakov
(2010),Rosenbaum and Tsybakov (2013) and Belloni et al. (2016b) introduced
matrix uncertainty selectors and their improved versions, and proposed to
solve by linear or conic programming; Belloni et al. (2017a, b) further
studied the confidence band and pivotal estimation, respectively; Loh and
Wainwright (2012) introduced unbiased surrogates for
${\mathbf{X}}{\mathbf{X}}^{\top}$ and $y^{\top}{\mathbf{X}}$ and cast the
estimation into non-convex optimization problems; Datta and Zou (2017)
proposed convex conditioned Lasso estimator under additive and multiplicative
noise settings; recently, Rudelson et al. (2017) discussed the regression with
correlated error-in-variables.

The focus of high-dimensional error-in-variable regression has so far been
mainly on homoskedastic Gaussian, sub-Gaussian, or bounded corruption setting.
Motivated by applications in high-throughput sequencing in microbiome studies,
we consider a new framework of high-dimensional error-in-variable regression
that adapts to compositional covariates in this paper. Before detailed
explanations, we first discuss the background of microbiome data analysis and
the setup of high-dimensional log-contrast model with sequencing data.

###  1.1 Regression Analysis for Microbial Compositional Data

The human microbiome is the aggregate of all microbes that reside on human
bodies. It has attracted enormous recent attention due to its association with
human health (The Human Microbiome Project Consortium, 2012). Recent studies
found that human microbiome may be closely related with various diseases, such
as cancer (Schwabe and Jobin, 2013), Crohn’s disease (Lewis et al., 2015), and
obesity (Turnbaugh et al., 2006, 2009). Modern next-generation sequencing
technologies, such as 16S ribosomal RNA and shotgun metagenomics sequencing,
provide quantification of the human microbiome by performing direct DNA
sequencing on either whole metagenomes or individual marker genes. By aligning
sequence reads to referential microbial genomes, we can organize the
sequencing data into a count matrix with rows representing samples and columns
representing microbial taxa or genes. Such data can be seen as the random
realization of relative abundance of bacteria in each sample.

To account for the difference in sequencing depth, i.e., total read count,
across samples (see Li (2015) for a survey and the references therein), the
sequencing read counts of bacterial genes are often normalized into
compositions. The resulting data, also called compositional data, pose
statistical challenges due to the collinearity and non-normality that come
from their compositional nature. To address these issues, Aitchison and Bacon-
shone (1984) introduced the log-contrast model:

$y_{i}=\sum_{j=1}^{p-1}\log(Z_{ij}/Z_{ip})\beta_{j}^{\ast}+\varepsilon_{i},\quad
i=1,\dots,n.$ (1)

Here, $W_{ij}$ and $Z_{ij}=W_{ij}/(\sum_{j^{\prime}=1}^{p}W_{ij^{\prime}})$
are the absolute count and the relative abundance of the $j$th component (e.g.
bacterial gene or taxon) in the $i$th sample, respectively;
$y=(y_{1},\ldots,y_{n})^{\top}$ is the regression response;
$\varepsilon=(\varepsilon_{1},\ldots,\varepsilon_{n})^{\top}$ is the noise;
$n$ is the sample size and $p$ is the number of components. The analysis of
log-contrast model (1) is often dependent on the choice of reference component
$Z_{ip}$, especially in high-dimensional settings. To address this issue, Lin
et al. (2014) reformulated (1) by introducing
$\beta_{p}^{\ast}=-\sum_{j=1}^{p-1}\beta_{j}^{\ast}$,

$y_{i}=\sum_{j=1}^{p}\log(Z_{ij})\beta_{j}^{\ast}+\varepsilon_{i},\quad
i=1,\dots,n,\quad\mbox{ subject to }\sum_{j=1}^{p}\beta_{j}^{\ast}=0$ (2)

and estimated $\beta^{\ast}$ through constrained $l_{1}$ regularized estimator

$\hat{\beta}=\mathop{\rm
arg\min}_{\beta:\sum_{j}\beta_{j}=0}\left\\{\frac{1}{2n}\sum_{i=1}^{n}\bigg{(}y_{i}-\sum_{j=1}^{p}\log(Z_{ij})\beta_{j}\bigg{)}^{2}+\lambda||\beta||_{1}\right\\}.$
(3)

More recently, Shi et al. (2016) studied the statistical inference and
confidence intervals for $\beta^{\ast}$, Wang and Zhao (2017) considered the
subcomposition selection in compositional data regression via a tree-guided
regularization method, and Lu et al. (2018) investigated the logistic
regression of compositional data.

The direct application of Models (1) and (2) by normalizing sequencing read
counts, i.e., using $Z_{ij}=W_{ij}/(\sum_{j=1}^{p}W_{ij})$ as covariates, has
several drawbacks. Firstly, it ignores the fact that $Z_{ij}$’s are random
realizations rather than true compositions of the components. In next-
generation sequencing data, $Z_{ij}$ is the proportion of read count of
component $j$ among all components in sample $i$, and is thus a transformation
of discrete random variables that reflect the underlying true compositions
with measurement errors. As mentioned earlier, overlooking the measurement
error in regressors may lead to inaccurate results. By treating $Z_{ij}$ as
the true compositions, it is also overlooking the heteroskedasticity or
overdispersion of $Z_{ij}$ caused by enormous uncontrollable factors of
variation in sequencing, e.g., time, sampling location, or technical
variability (Chen and Li, 2013). Secondly, the procedure requires $Z_{ij}>0$
while in reality, compositional data from next-generation sequencing often
contain a lot of zeros due to rarity of certain components. Strategies to deal
with the zeros include replacing zero counts by a subjectively chosen small
number, such as 0.5, before normalizing counts into compositions (Aitchison,
2003; Martin-Fernandez et al., 2000; Martín-Fernández et al., 2003, 2015), or
imputing the entire composition matrix (Cao et al., 2017) based on low-rank
assumption. However, to the best of our knowledge, there is still no consensus
on the best approach to deal with zero read counts in compositional data
regression.

###  1.2 High-Dimensional Log-Error-in-Variable Regression

To address the aforementioned problems in compositional data regression, we
introduce a _high-dimensional log-error-in-variable regression_ model that
directly handles count covariates without normalization into compositions or
imputation of zeros. Recall $W_{ij}$ is the count of the $j$th component in
the $i$th sample. We assume $W_{i}=(W_{i1},\ldots,W_{ip})^{\top}$ to satisfy
Dirichlet-multinomial distribution (Mosimann, 1962) given the total count
$N_{i}=\sum_{j=1}^{p}W_{ij}$ in the $i$th sample,

$\begin{split}W_{i}|N_{i}&\sim\text{Dirichlet-
Multinomial}(N_{i},\alpha_{i}X_{i1},\ldots,\alpha_{i}X_{in}).\\\ \end{split}$
(4)

where $N_{i}$ follows $\text{Poisson}(\nu_{i})$ to account for the randomness
of sequencing depth, $X_{i}=(X_{i1},\dots,X_{ip})^{\top}$ is the underlying
true composition of the $p$ components, and $\alpha_{i}$ is the overdispersion
parameter of the subject from which the $i$th sample is measured. When
$\alpha_{i}$ goes to infinity, the Dirichlet-multinomial distribution
degenerates to the regular multinomial distribution. More detailed discussions
on Dirichlet-multinomial distribution is given in Section 2.

Since observable counts $W_{i}$ are merely realizations of underlying
compositions $X_{i}$, it is more reasonable to assume association between
$y_{i}$ and $X_{i}$s rather than between $y_{i}$ and $W_{i}$s. We thus assume
the regression response $y_{i}$ to be dependent on $X_{i}$s through the
following log-contrast model,

$y_{i}=\sum_{j=1}^{p}\log(X_{ij})\beta_{j}^{\ast}+\varepsilon_{i},\quad
i=1,\ldots,n,\quad\mbox{subject to}\quad\sum_{j=1}^{p}\beta_{j}^{\ast}=0.$ (5)

We refer to (4) together with (5) as _log-error-in-variable regression model_
, which is well-suited for regression with read count covariates in microbial
studies. Our aim is to estimate $\beta^{\ast}$ based on responses
$y\in\mathbb{R}^{n}$ and error-in-covariates $W\in\mathbb{R}^{n\times p}$.
Most of the results on error-in-variables regression, e.g., Rosenbaum and
Tsybakov (2010, 2013); Belloni et al. (2016a, b); Datta and Zou (2017), deal
with homoscedastic continuous variables and may not be directly applied here
since the $W_{i}$’s are discrete random variables with heteroscedasticity
depending on $X_{i}$ and $\alpha_{i}$. Therefore, new methods are in need for
the estimation of log-error-in-variable model.

In this paper, we propose a surprisingly simple and straightforward estimation
scheme, named _variable correction regularized estimator_ , for the high-
dimensional log-error-in-variable regression. In particular, when the count
observations are without overdispersion, we propose to _add 0.5_ to all counts
$W_{ij}$, then estimate the regression parameters using constrained Lasso; for
overdispersed data, we propose to _add an amount related to the overdispersion
level_ to $W_{ij}$ to alleviate the effect of any overly large or small counts
due to overdispersion. We show that the proposed methods achieves minimax
optimal performance in a general class of settings.

In addition, we consider a more general log-error-in-variable regression model
to accommodate broader types of settings.

$\begin{gathered}y_{i}=\sum_{j=1}^{p}\log(\nu_{ij})\beta_{j}^{\ast}+\varepsilon_{i},\quad
i=1,\ldots,n,\quad\text{ subject to
}\quad{\mathbf{C}}^{\top}\beta^{\ast}=0,\\\
W_{ij}\sim{\mathbb{P}}_{\nu_{ij}}\text{ independently}.\end{gathered}$ (6)

Here, ${\mathbb{P}}_{\nu_{ij}}$ is a general class of distributions with mean
parameter $\nu_{ij}$, $\varepsilon_{i}$ are i.i.d. sub-Gaussian noises with
mean zero and variance $\sigma^{2}$, and $\mathbf{C}$ is a $p\times r$ matrix
with each column representing a linear constraint on the regression
coefficients. We provide the variable correction regularized estimator for
$\beta^{\ast}$ in (6) with theoretical guarantees.

###  1.3 Organization of the Paper

The rest of the article is organized as follows. In Section 2, we introduce
the variable correction regularized estimator under the high-dimensional log-
error-in-variable regression model described by (4) and (5). In Section 3, we
provide theoretical guarantee for the estimation accuracy when overdispersion
is present and absent, respectively. In Section 4, we discuss the more general
error-in-variable model (6), for which we provide both an estimation procedure
and theoretical analysis. A real data application in a longitudinal human gut
microbiome study is presented in Section 5.1. Simulation results are provided
in Section 5.2. Finally, we collect the proofs in Section 7 and the
supplementary materials.

##  2 Methods for Log-Error-in-Variable Regression

We first introduce the notations that will be used in this paper. For any
vector $v$, let $\|v\|_{q}=(\sum_{i}|v_{i}|^{q})^{1/q}$ be the $\ell_{q}$
norm. Especially, $\|v\|_{\infty}=\max_{i}|v_{i}|$. Let $1_{p}$ be the
$p$-dimensional vector with all entries as ones. We use boldface letters, such
as ${\mathbf{W}},{\mathbf{C}}$, to denote matrices. Let ${\mathbf{I}}_{p}$ be
the $p$-by-$p$ identity matrix. For any matrix ${\mathbf{W}}$, define
$\|{\mathbf{W}}\|_{q}=\max_{v}\|{\mathbf{W}}v\|_{q}/\|v\|_{q}$ to be the
$\ell_{q}$ operator norm of ${\mathbf{W}}$. The uppercase $C$, lowercase $c$,
and their variations, e.g., $C_{1},C_{2},c_{0},c^{\prime}$, are used to
represent generic large and small constants, respectively. The actual values
of these symbols for constants may vary across different lines.

To estimate $\beta^{\ast}$ in the log-error-in-variable model (5), one classic
method is to perform simple normalization, i.e., to use $W_{ij}/N_{i}$ as a
surrogate for $X_{ij}$ and implement the classic high-dimensional regularized
estimators with $\log(W_{ij}/N_{i})$ as covariates. As discussed in the
introduction, this idea has two critical issues: first, the zero-valued
$W_{ij}$’s need to be replaced by a small value to make them positive in the
log transformation. The choice of this value is often difficult but critical
to the performance of the final estimates; secondly, even though
${\mathbb{E}}(W_{ij}/N_{i}|N_{i})=X_{ij}$, $\log(W_{ij}/N_{i})$ may be a
biased estimator for $\log(X_{ij})$, which can cause additional inaccuracy to
the regression analysis. To further illustrate the biasness of
$\log(W_{ij}/N_{i})$ and to introduce our fixing plan, we first focus on the
non-overdispersion case, i.e. when $W_{i}|N_{i}$ follows the regular
multinomial distribution:

$\begin{split}&(W_{i1},\ldots,W_{ip})|N_{i}\sim\text{Multinomial}(N_{i},X_{i1},\ldots,X_{ip}),\\\
\end{split}$
$\begin{split}&{\mathbb{P}}\left\\{(W_{i1},\ldots,W_{ip})=(k_{i1},\ldots,k_{ip})\right\\}=\frac{N_{i}!}{k_{i1}!\cdots
k_{ip}!}\prod_{j=1}^{p}X_{ij}^{k_{ij}},\\\
\text{where}\quad&k_{i1}+\cdots+k_{ip}=N_{i},\quad
k_{i1},\ldots,k_{ip}\in\\{0,1,2,\ldots\\}.\end{split}$

In this case, we have $W_{ij}$ following Poisson($\nu_{i}X_{ij}$) and

${\mathbb{E}}(W_{ij})=\nu_{i}X_{ij},\quad{\rm
var}(W_{ij})=\nu_{i}X_{ij}(1-X_{ij})+\nu_{i}X_{ij}^{2}.$

For any $z_{i}\geq 0$, the Taylor expansion of $\log(W_{ij}+z_{i})$ at
$\nu_{i}X_{ij}$ yields the following approximation,

$\begin{split}{\mathbb{E}}\log(W_{ij}+z_{i})&\approx\log(\nu_{i}X_{ij})+\dfrac{{\mathbb{E}}(W_{ij}-\nu_{i}X_{ij}+z_{i})}{\nu_{i}X_{ij}}-\dfrac{{\rm
Var}(W_{ij})+2z_{i}{\mathbb{E}}(W_{ij}-\nu_{i}X_{ij})+z_{i}^{2}}{2\nu_{i}^{2}X_{ij}^{2}}\\\
&=\log(\nu_{i}X_{ij})+\dfrac{z_{i}-1/2}{\nu_{i}X_{ij}}-\dfrac{z_{i}^{2}}{2\nu_{i}^{2}X_{ij}^{2}}.\end{split}$

Since $N_{i}$ is the total number of reads in sample $i$ and is generally
large in practice (e.g., around $10^{4}$ to $10^{5}$ in our real data
example), we can assume $\nu_{i}={\mathbb{E}}N_{i}$ to be large. Then,

$\begin{split}{\mathbb{E}}\log\left((W_{ij}+z_{i})/N_{i}\right)-\log(X_{ij})\approx&\dfrac{z_{i}-1/2}{\nu_{i}X_{ij}}-\dfrac{z_{i}^{2}}{2\nu_{i}^{2}X_{ij}^{2}}+\log(\nu_{i})-{\mathbb{E}}\log(N_{i})\\\
\approx&\dfrac{z_{i}-1/2}{\nu_{i}X_{ij}}-\dfrac{z_{i}^{2}}{2\nu_{i}^{2}X_{ij}^{2}}.\end{split}$

We can see the bias of $\log((W_{ij}+z_{i})/N_{i})$ for estimating
$\log(X_{ij})$ is approximately $-(1/2)\nu_{i}X_{ij}$ and
$-(1/8)\nu_{i}^{2}X_{ij}^{2}$ when $z_{i}=0$ and $z_{i}=1/2$, respectively.
For large $\nu_{i}$, one has
$(1/8)\nu_{i}^{2}X_{ij}^{2}\ll(1/2)\nu_{i}X_{ij}$. Therefore, heuristically
$\log((W_{ij}+1/2)/N_{i})$ is a significantly less biased estimator for
$\log(X_{ij})$ compared to $\log(W_{ij}/N_{i})$. Figure 1 illustrates the bias
of $\log(W\vee 0.5)$ (i.e., replacing zeros by $1/2$) and
$\log(W+c),c=1/4,1/2,3/4,1$ for estimating $\log(\nu)$ when $W$ follows
Poisson($\nu$). The plot suggests that $\log(W+1/2)$ achieves the minimum bias
among these choices. In addition, by adding the positive value $1/2$ to all
$W_{ij}$s, the earlier mentioned zero-replacement issue is simultaneously
solved!

Figure 1: $\log(\nu)$ vs. $\mathbb{E}\log(W\vee 0.5)$ (zero-replace) and
$\mathbb{E}\log(W+c)$ for $c=1/4,1/2,3/4,1$. Here, $W\sim\text{Poisson}(\nu)$.

To account for higher variability in the count data, we also consider the
overdispersed case where $W_{i}|N_{i}$ satisfies Dirichlet-multinomial
distribution (Mosimann, 1962),

$(W_{i1},\ldots,W_{ip})|N_{i}\sim\text{Dirichlet-
Multinomial}(N_{i},\alpha_{i}X_{i1},\ldots,\alpha_{i}X_{ip}),$
$\begin{split}&{\mathbb{P}}\left\\{(W_{i1},\ldots,W_{ip})=(k_{i1},\ldots,k_{ip})\Big{|}(N_{i},q_{i1},\ldots,q_{ip})\right\\}=\frac{N_{i}!}{k_{i1}!\cdots
k_{ip}!}\prod_{j=1}^{p}q_{ij}^{k_{ij}},\\\
&f(q_{i1},\ldots,q_{ip}|N_{i})=\frac{1}{B(\alpha_{i}X_{i1},\ldots,\alpha_{i}X_{ip})}\prod_{j=1}^{p}q_{ij}^{\alpha_{i}X_{ij}-1},\\\
\text{where}\quad&k_{i1}+\cdots+k_{ip}=N_{i},k_{i1},\ldots,k_{ip}\in\\{0,1,2,\ldots\\},\quad
q_{i1}+\cdots+q_{ip}=1,q_{i1},\ldots,q_{ip}\geq 0,\end{split}$

Here,
$B(\alpha_{i}X_{i1},\ldots,\alpha_{i}X_{ip})=\frac{\prod_{j=1}^{p}\Gamma(\alpha_{i}X_{ij})}{\Gamma(\alpha_{i})}$
is the Beta function and $\alpha_{i}$ represents the level of overdispersion.
When $\alpha_{i}<\infty$, we have

${\mathbb{E}}(W_{ij})=\nu_{i}X_{ij},\quad{\rm
var}(W_{ij})=\dfrac{\nu_{i}+\alpha_{i}+1}{\alpha_{i}+1}\cdot\nu_{i}X_{ij}(1-X_{ij})+\nu_{i}X_{ij}^{2}.$

Similarly, by investigating the Taylor expansion of
${\mathbb{E}}\log((W_{ij}+z_{i})/N_{i})$, it can be shown that taking
$z_{i}=\frac{N_{i}+\alpha_{i}+1}{2(\alpha_{i}+1)}$ will make
$\log((W_{ij}+z_{i})/N_{i})$ a better estimator for $\log(X_{ij})$ (more
rigorous argument is postponed to the proof of Lemma 5 in the supplementary
materials). It is noteworthy that $z_{i}$ is an estimate for half of
$\frac{\nu_{i}+\alpha_{i}+1}{\alpha_{i}+1}$, which quantifies the
overdispersion rate of $W_{ij}$ compared with the multinomial distribution.

These heuristic arguments inspire us to the following _variable correction
regularized estimator_ for the log-error-in-variable regression in (4) and
(5):

$\begin{split}&\hat{\beta}=\mathop{\rm
arg\min}_{\beta}\left(\frac{1}{2n}\|y-{\mathbf{B}}_{W}\beta\|_{2}^{2}+\lambda\|\beta\|_{1}\right)\quad\text{subject
to}\quad\sum_{j=1}^{p}\beta_{j}=0,\\\
&\text{where}\quad{\mathbf{B}}_{W}\in\mathbb{R}^{n\times
p}\quad\text{and}\quad({\mathbf{B}}_{W})_{ij}=\log\left(W_{ij}+\frac{N_{i}+\alpha_{i}+1}{2(\alpha_{i}+1)}\right).\end{split}$
(7)

Particularly if $\alpha_{i}=\infty$, i.e., $W_{i\cdot}|N_{i}$ satisfies
multinomial without overdispersion,
$({\mathbf{B}}_{W})_{ij}=\log(W_{ij}+1/2)$.

######  Remark 1 (Interpretation of Variable Correction).

Different from the classic zero-replacement scheme that replaces only the zero
covariates by a fixed value, we propose to add $1/2$ to all covariates in the
non-overdispersion case. For overdispersed sample, say $W_{i}$, we propose to
correct with a larger value:
$z_{i}=\frac{N_{i}+\alpha_{i}+1}{2(\alpha_{i}+1)}$. In particular, with
smaller total count $N_{i}$ or larger degree of overdispersion (i.e, smaller
$\alpha_{i}$), the observable count covariates $W_{i}$ contains much noisier
information about the true underlying composition $X_{i}$. Thus, we add larger
values to alleviate the effect of overly large or small counts due to
overdispersion.

######  Remark 2.

When there are multiple samples $W_{i}$ that share the same $X_{i}$ and
$\alpha_{i}$ and when there is evidence of overdispersion in practice,
$\alpha_{i}$ can be estimated by method of moment estimator (La Rosa et al.,
2012; Mosimann, 1962; Tvedebrink, 2010; Weir and Hill, 2002) or maximum
likelihood estimator (Tvedebrink, 2010). Otherwise, $\alpha_{i}=+\infty$ and
$z_{i}=1/2$ are suggested. The more detailed discussions on the method of
moment estimator of $\alpha_{i}$ is given in Section A of the supplementary
materials.

##  3 Theoretical Analysis

In this section, we investigate the theoretical performance of the proposed
variable correction regularized estimator for the log-error-in-variable
regression model. We first focus on the case of $(W_{i1},\ldots,W_{in})|N_{i}$
following multinomial distribution without overdispersion.

For convenience, denote $\bar{\nu}=\sum_{i=1}^{n}\nu_{i}/n$ and
$\nu=(\nu_{1},\dots,\nu_{n})^{\top}$. We say a matrix
${\mathbf{M}}\in\mathbb{R}^{n\times p}$ satisfies the restricted isometry
property (RIP) (Candes and Tao, 2005; Candes et al., 2007) with constant
$\delta_{s}({\mathbf{M}})\in(0,1)$ if for all $s$-sparse vector
$\beta\in\mathbb{R}^{p}$,

$n(1-\delta_{s}({\mathbf{M}}))\|\beta\|_{2}^{2}\leq\|{\mathbf{M}}\beta\|_{2}^{2}\leq
n(1+\delta_{s}({\mathbf{M}}))\|\beta\|_{2}^{2}.$ (8)

The RIP condition has become one of the most commonly used regularity
condition in high-dimensional regression literature since its introduction.
Denote ${\mathbf{B}}_{W}\in\mathbb{R}^{n\times p}$ with
$({\mathbf{B}}_{W})_{ij}=\log\left(W_{ij}+\frac{N_{i}+\alpha_{i}+1}{2(\alpha_{i}+1)}\right)$
as the corrected design matrix. We assume the centralized and corrected design
matrix
$\bar{{\mathbf{B}}}_{W}={\mathbf{B}}_{W}({\mathbf{I}}_{p}-\frac{1}{p}1_{p}1_{p}^{\top})$
satisfies the RIP condition with high probability.

######  Condition 1 (RIP Condition).

$\bar{{\mathbf{B}}}_{W}$ satisfies RIP condition with constant
$\delta_{2s}(\bar{{\mathbf{B}}}_{W})<1/10$ with probability
$1-\epsilon^{\prime}$ for some small quantity $\epsilon^{\prime}$.

######  Remark 3.

Denote ${\mathbf{V}},\bar{{\mathbf{V}}}\in\mathbb{R}^{n\times p}$ where
$(\mathbf{{\mathbf{V}}})_{ij}=\log(\nu_{i}X_{ij})$ and
$\bar{\mathbf{{\mathbf{V}}}}={\mathbf{V}}({\mathbf{I}}_{p}-\frac{1}{p}1_{p}1_{p}^{\top})$.
Lemma 14 in the supplementary materials shows that Condition 1 holds if the
deterministic matrix $\bar{\mathbf{{\mathbf{V}}}}$ satisfies RIP condition
with constant $\delta_{2s}(\bar{\mathbf{{\mathbf{V}}}})<1/20$ in the case
without overdispersion.

The following sample complexity condition matches the classic results in high-
dimensional regression (see, e.g., Candes et al. (2007); Bickel et al.
(2009)).

######  Condition 2 (Sample Complexity).

Assume a sufficient number of observations can be observed: $n\geq Cs\log p$
for some large constant $C$.

Based on these conditions, we show the proposed variable correction
regularized estimator satisfies the following upper bound with high
probability.

######  Theorem 1 (Upper Bound).

Consider the log-error-in-variable regression model (4) and (5), where
$\alpha_{i}=\infty$, i.e. ${\mathbf{W}}$ has no overdispersion. Suppose
Conditions 1 and 2 hold, and $a\bar{\nu}\leq\nu_{i}\leq b\bar{\nu}$, $a/p\leq
X_{ij}\leq b/p$ for constants $0<a<1<b$. If for some large constant $C>0$,
some $\epsilon>0$, and a constant $C_{\epsilon}$ that only depends on
$\epsilon$, we have
$\bar{\nu}\geq\max\\{p\cdot(sn)^{1/2+\epsilon},Cp\log(np),C_{\epsilon}/n\\}$,
then by choosing $\lambda=C\sqrt{\frac{\log
p}{n}\left(\sigma^{2}+\frac{p}{\bar{\nu}}\|\beta^{*}\|_{2}^{2}\right)}$ for
some large constant $C>0$, the variable correction regularized estimator (7)
satisfies

$\|\hat{\beta}-\beta^{*}\|_{2}^{2}\leq\frac{Cs\log
p}{n}\left(\sigma^{2}+\frac{p}{\bar{\nu}}\|\beta^{*}\|_{2}^{2}\right)$ (9)

with probability at least $1-4p^{-C^{\prime}}-\epsilon^{\prime}$.

######  Remark 4.

Theorem 1 shows that the estimation error upper bound decreases with larger
sample size $n$, smaller dimension $p$, smaller noise variance $\sigma^{2}$,
higher average expectation of count $\bar{\nu}$, smaller signal amplitude
$\|\beta^{\ast}\|_{2}^{2}$, or smaller sparsity level $s$.

######  Remark 5.

We briefly discuss the sketch of the highly non-trivial proof for Theorem 1
here. First, we develop a series of inequalities on the corrected covariates
in Lemmas 1-3 based on the tail probability bounds of Poisson, multinomial,
and sub-Gaussian distributions. Then we develop an upper bound for
$\left\|\bar{{\mathbf{B}}}_{W}^{\top}\left(\bar{{\mathbf{B}}}_{W}\beta^{*}-y\right)\right\|_{\infty}$,
a pivotal term in high-dimensional regression analysis, in Lemma 4. Finally,
we combine these inequalities with Condition 1 and obtain the upper bound for
estimation error. The biggest challenge in the proof is to bound the corrected
covariates $\log(W_{ij}+\frac{1}{2})$. Although the aforementioned Taylor’s
expansion heuristically show $\log\left(W_{ij}+\frac{1}{2}\right)$ are good
estimators for $\log(\nu_{ij})$, more careful analysis is needed to obtain
rigorous upper bounds on their biases and variances. To this end, we perform
truncation on $W_{ij}$ as the direct application of Poisson tail bounds may
not yield sharp enough results.

######  Remark 6.

We introduce the entry-wise upper and lower bounds of $\nu_{i}$ and $X_{ij}$
in Theorem 1 mainly to regularize the tail probability of $W_{ij}$ that
enables us to derive the minimax-optimal error bounds for the proposed
procedure. Actually, similar conditions appeared in a range of literature on
Poisson and multinomial inverse problems, e.g., Poisson matrix completion (Cao
and Xie, 2016, Equation (10)), Poisson sparse regression (Jiang et al., 2015,
Assumption 2.1), composition matrix estimation from sparse count data (Cao et
al., 2017, Remark 3), and Point autoregressive model (Hall et al., 2016,
Definition of $\mathcal{A}_{s}$) among many others. Such assumptions are
mainly for technical purposes and may not be as crucial in practice as
illustrated by the forthcoming real data and simulation studies.

To derive the lower bound for estimation error, we consider the following
class of covariate matrices and parameter vectors,

$\begin{split}\mathcal{F}_{p,n,s}(R,Q)=\bigg{\\{}(\nu,{\mathbf{X}},\beta):&a\bar{\nu}\leq\nu_{i}\leq
b\bar{\nu},a/p\leq X_{ij}\leq b/p\text{ for constants }0<a<1<b;\\\
&\bar{\mathbf{V}}\text{ satisfies RIP conditon with
}\delta_{2s}(\bar{\mathbf{V}})<\frac{1}{20};\\\ &||\beta||_{2}\leq
R,1_{p}^{\top}\beta=0,e^{-\frac{3}{2}}Q\leq\bar{\nu}\leq
e^{\frac{3}{2}}Q\bigg{\\}}.\\\ \end{split}$ (10)

The constraints in $\mathcal{F}_{p,n,s}(R,Q)$ correspond to the regularization
assumptions in Theorem 1. The upper bound in Theorem 1 turns out to match the
minimax lower bound in $\mathcal{F}_{p,n,s}(R,Q)$.

######  Theorem 2 (Lower Bound).

Suppose $\varepsilon_{i}\overset{iid}{\sim}N(0,\sigma^{2})$. If we have $n\geq
Cs\log p$ for some large constant $C>0$,
$R\geq\bar{c}\sqrt{\frac{s\log\left(p/s\right)}{n}\sigma^{2}}$ for some
constant $\bar{c}>0$, $Q\geq p$, and $s\geq 4$, then

$\inf_{\hat{\beta}}\sup_{(\nu,{\mathbf{X}},\beta)\in\mathcal{F}_{p,n,s}(R,Q)}{\mathbb{E}}\left\|\hat{\beta}-\beta\right\|_{2}^{2}\geq\frac{cs\log(p/s)}{n}\left(\sigma^{2}+\frac{p}{Q}R^{2}\right).$
(11)

######  Remark 7.

The lower bound in Theorem 2 consists of two terms:
$\frac{cs\log(p/s)}{n}\sigma^{2}$ and
$\frac{cs\log(p/s)}{n}\cdot\frac{p}{Q}R^{2}$, which originate from the
uncertainty of $\varepsilon$ and ${\mathbf{W}}$, respectively. We thus show
$\inf\sup\mathbb{E}\|\hat{\beta}-\beta\|_{2}^{2}\geq\frac{cs\log(p/s)}{n}\sigma^{2}$
and
$\inf\sup\mathbb{E}\|\hat{\beta}-\beta\|_{2}^{2}\geq\frac{cs\log(p/s)}{n}\cdot\frac{p}{Q}R^{2}$
separately in the proof of this theorem. While the first inequality develops
from the classic high-dimensional regression literature (see, e.g., Rigollet
and Tsybakov (2011)), the proof for the second one is far more complicated. In
particular, we construct a series of instances of
$(\nu^{(i)},{\mathbf{X}}^{(i)},\beta^{(i)})_{i=1}^{N}$ satisfying the
constraints in (10) and ${\rm
diag}(\nu^{(i)}){\mathbf{X}}^{(i)}\beta^{(i)}={\rm
diag}(\nu^{(j)}){\mathbf{X}}^{(j)}\beta^{(j)}$ for all $i,j$. By such the
design, $y^{(i)}={\rm
diag}(\nu^{(i)}){\mathbf{V}}^{(i)}\beta^{(i)}+\varepsilon^{(i)}$ becomes
“nullified” for estimating $\beta^{(i)}$ and we show these instances are non-
distinguishable based a sample of $({\mathbf{W}}^{(i)},y^{(i)})$. Then we
apply the generalized Fano’s lemma and establish the desired error lower
bound.

######  Remark 8.

Theorems 1 and 2 together show that the proposed procedure (7) achieves the
optimal rate of convergence in $\mathcal{F}_{p,n,s}(R)$ as long as $\log(p)$
and $\log(ep/s)$ are in the same order.

Next, we further consider the overdispersed case where $\alpha_{i}<\infty$ in
Model (4) and (5). We have the following upper bound for estimation accuracy
of the proposed variable correction regularized estimator (7) with covariates
${\mathbf{B}}_{W}:({\mathbf{B}}_{W})_{ij}=\log\left(W_{ij}+\frac{N_{i}+\alpha_{i}+1}{2(\alpha_{i}+1)}\right)$.

######  Theorem 3 (Upper Bound for Overdispersed Log-Error-in-Variable
Regression).

Suppose Conditions 1 and 2 hold, and $a\bar{\nu}\leq\nu_{i}\leq b\bar{\nu}$,
$a/p\leq X_{ij}\leq b/p$ for constants $0<a<1<b$. Set
$\zeta_{\max}=\max_{i}\zeta_{i}$, where
$\zeta_{i}=\frac{\nu_{i}+\alpha_{i}+1}{2(\alpha_{i}+1)}$ represents the level
of overdispersion for $i$th sample. If for some $\delta>0$, some large
constant $C$, and a large constant $C(\delta)$ that only depends on $\delta$,
we have

$\nu_{ij}\geq\zeta_{i}^{1+\delta},\quad\frac{\bar{\nu}}{\zeta_{\max}}\geq
p\cdot\max\left\\{C\log(np),\frac{4n}{\log
p}\left(\frac{1+\delta}{\delta}\right)^{2}\log^{2}\left(\frac{n}{\log
p}\left(\frac{1+\delta}{\delta}\right)^{2}\right),C(\delta)\right\\},$

and $n\geq Cs\log p$, then by choosing $\lambda=C\sqrt{\frac{\log
p}{n}}\left(\sigma+\left(\frac{p}{\bar{\nu}}\zeta_{\max}\right)^{\frac{1}{2}}\|\beta^{*}\|_{1}\right)$,
we have

$\|\hat{\beta}-\beta\|_{2}^{2}\leq\frac{Cs\log
p}{n}\left(\sigma^{2}+\frac{p}{\bar{\nu}}\zeta_{\max}\cdot\|\beta^{*}\|_{1}^{2}\right)$
(12)

with probability $1-6p^{-C^{\prime}}-\epsilon^{\prime}$, where $C^{\prime}$ is
a large constant.

##  4 General High-dimensional Log-Error-in-Variable Regression

In this section, we further extend the discussion to general high-dimensional
log-error-in-variable regression that accommodates broader scenarios.
Specifically, let

$\begin{gathered}y={\mathbf{V}}\beta^{\ast}+\varepsilon,\quad\text{or
equivalently}\quad
y_{i}=\sum_{j=1}^{p}\log(\nu_{ij})\beta_{j}^{*}+\varepsilon_{i}\quad\text{
subject to }\quad{\mathbf{C}}^{\top}\beta^{*}=0,\\\
\text{where}\quad{\mathbf{W}}=(W_{ij}),\quad
W_{ij}\sim{\mathbb{P}}_{\nu_{ij}}\text{ independently},\quad
i=1,\ldots,n,j=1,\ldots,p.\end{gathered}$ (13)

Here, ${\mathbf{V}}=(\log(\nu_{ij}))_{\begin{subarray}{c}1\leq i\leq n,1\leq
j\leq p\end{subarray}}$ are unknown underlying covariates,
${\mathbb{P}}_{\nu}$ is some general class of distributions with mean
parameter $\nu$, $\varepsilon_{i}$ are i.i.d. sub-Gaussian noises with mean
zero and variance $\sigma^{2}$, $\beta^{\ast}$ is the sparse parameter of
interest, and $\mathbf{C}$ is a $p\times r$ matrix with each column
representing a linear constraint on the regression coefficients. We aim to
estimate $\beta^{\ast}$ based on $y$ and ${\mathbf{W}}$.

Assume that $W\sim{\mathbb{P}}_{\nu}$. Suppose one can find a good estimate
for $\log(\nu)$, say $\phi(W)$. Based on previous discussions, one can choose
$\phi(W)=\log\left(W+\frac{1}{2}\right)$ if ${\mathbb{P}}_{\nu}$ is Poisson
distribution with parameter $\nu$. Some other choices include:
$\phi(\nu)=\log\left((\nu+\frac{\gamma}{2})\vee 1\right)$ if
${\mathbb{P}}_{\nu}\sim N(\nu,\gamma\nu)$; $\phi(\nu)=\log(\nu+\gamma/2)$ if
${\mathbb{P}}_{\nu}\sim\text{Gamma}(\nu,\gamma)$ (i.e., Gamma distribution
with shape parameter $\nu$, scale parameter $\gamma$, and pdf
$f(x)=\frac{1}{\Gamma(\nu)\gamma^{\nu}}x^{\nu-1}e^{-\frac{x}{\gamma}}$).

Similarly to Section 2, we propose the following variable correction
regularized estimator for $\beta^{*}$ in general high-dimensional log-error-
in-variable regression model (13),

$\begin{split}\hat{\beta}=\mathop{\rm
arg\min}_{\beta}\left(\frac{1}{2n}\|y-{\mathbf{B}}_{W}\beta\|_{2}^{2}+\lambda\|\beta\|_{1}\right)\quad\text{subject
to }\mathbf{C}^{\top}\beta=0.\end{split}$ (14)

Here, ${\mathbf{B}}_{W}\in\mathbb{R}^{n\times p}$ with
$({\mathbf{B}}_{W})_{ij}=\phi(W_{ij})$ and $\lambda$ is some tuning parameter.
Define
${\mathbf{P}}_{\mathbf{C}}={\mathbf{C}}({\mathbf{C}}^{\top}{\mathbf{C}})^{\dagger}{\mathbf{C}}^{\top}$
as the projection matrix onto the column space of ${\mathbf{C}}$ where
$(\cdot)^{\dagger}$ represents matrix pseudoinverse. We introduce the
following regularity conditions for the theoretical analysis of the proposed
estimator.

######  Condition 3.

$\|{\mathbf{I}}_{p}-{\mathbf{P}}_{{\mathbf{C}}}\|_{\infty}\leq k_{0}$ for a
constant $k_{0}$ that is free of $p$. Also, $1_{p}$ belongs to the column
space of ${\mathbf{C}}$, i.e.,
$\left({\mathbf{I}}_{p}-{\mathbf{P}}_{{\mathbf{C}}}\right)1_{p}=0$.

######  Condition 4.

$\bar{{\mathbf{B}}}_{W}={\mathbf{B}}_{W}({\mathbf{I}}_{p}-{\mathbf{P}}_{{\mathbf{C}}})$
satisfies RIP condition (8) with $\delta_{2s}(\bar{{\mathbf{B}}}_{W})<1/10$
with probability at least $1-\epsilon^{\prime}$ for some small quantity
$\epsilon$.

Conditions 3 and 4 are regularization assumptions for the constraint and
covariates, respectively. The requirement of
$\left({\mathbf{I}}_{p}-{\mathbf{P}}_{{\mathbf{C}}}\right)1_{p}=0$ in
Condition 3 can be omitted if we include an intercept in the regression
equation. Condition 4 can be seen as a parallel version of Condition 1 in the
previous section.

######  Condition 5.

Suppose $W\sim\mathbb{P}_{\nu}$. There exists a non-negative increasing
function $\pi(x)$ such that

  1. 1.

there exist a constant $K_{0}$ such that
$\left\|\sqrt{\pi(\nu)}\left(\phi(W_{ij})-{\mathbb{E}}\phi(W_{ij})\right)\right\|_{\psi_{2}}\leq
K_{0}$, where $\|X\|_{\psi_{2}}=\sup_{p\geq
1}p^{-1/2}(\mathbb{E}|X|^{p})^{1/p}$ is the Orlicz norm of random variables;

  2. 2.

$\frac{\pi(tx)}{\pi(x)}\leq C_{t}$ with $C_{t}$ only depending on $t\geq 1$;

  3. 3.

$\lim_{x\to\infty}\pi(x)=\infty$;

  4. 4.

$\max\left\\{\left|{\mathbb{E}}\phi^{2}(W)-\log^{2}\nu\right|,\left(\log\nu\left|{\mathbb{E}}\phi(W)-\log\nu\right|\right)\right\\}=o\left(1/\sqrt{\pi(\nu)}\right)$
as $\nu\to\infty$.

Condition 5 is introduced to regularize the tail property of distribution
$\mathbb{P}_{\nu}$ and the $\pi(\cdot)$ function essentially characterizes the
accuracy of estimator $\phi(W)$ to $\log(\nu)$. In particular, if
${\mathbb{P}}_{\nu}$ are Poisson$(\nu)$, Gamma$(\nu,\gamma)$, or
$N(\nu,\gamma\nu)$, $\pi(t)=t$ satisfies Condition 5.

The following theorem provides an upper bound for the variable correction
regularized estimator in the general high-dimensional log-error-in-variable
regression model.

######  Theorem 4 (General Upper Bound).

Suppose Conditions 2, 3, 4, and 5 hold, and in addition
$\left|\log\nu_{ij}-\log\nu_{kl}\right|\leq a$ for some constant $a>0$ for all
$1\leq i,k\leq n,1\leq j,l\leq p$. Denote
$\bar{\nu}=\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{p}\nu_{ij}$. If for some
uniform constant $C>0$ and constant $C_{n,s}$ that only depends on $s$ and
$n$, we have $\bar{\nu}>pC_{n,s}$, then by choosing $\lambda=C\sqrt{\frac{\log
p}{n}\left(\sigma^{2}+\frac{1}{\pi\left(\bar{\nu}/p\right)}\|\beta^{*}\|_{2}^{2}\right)}$,
we have

$\|\hat{\beta}-\beta^{\ast}\|_{2}^{2}\leq\frac{Cs\log
p}{n}\left(\sigma^{2}+\frac{1}{\pi\left(\bar{\nu}/p\right)}\|\beta^{*}\|_{2}^{2}\right)$
(15)

with probability $1-3p^{-C^{\prime}}-\epsilon$.

The following lower bound result suggests the optimality of the proposed
estimator.

######  Theorem 5 (General Lower Bound).

Suppose $\varepsilon_{i}\overset{i.i.d.}{\sim}N(0,\sigma^{2})$. Assume
$\varphi(x)$ is an increasing function with $\varphi(1)\geq 1$ and the
Kullback-Leibler divergence between $\mathbb{P}_{\nu_{1}}$ and
$\mathbb{P}_{\nu_{2}}$ satisfies
$D_{KL}\left({\mathbb{P}}_{\nu_{1}},{\mathbb{P}}_{\nu_{2}}\right)\leq
C\varphi(\nu_{1})\left(\log(\nu_{1})-\log(\nu_{2})\right)^{2}$ for all
$\nu_{1}$ and $\nu_{2}$ satisfying $|\log(\nu_{1})-\log(\nu_{2})|\leq 1$. If
$n\geq Cs\log p$ for some large constant $C>0$,
$R\geq\bar{c}\sqrt{\frac{s\log\left(\frac{p}{s}\right)}{n}\sigma^{2}}$ for
some constant $\bar{c}>0$ and $Q\geq p,s\geq 4$, then we have

$\inf_{\hat{\beta}}\sup_{(\mathbf{V},{\mathbf{C}},\beta)\in\mathcal{F}_{p,n,s}(R,Q)}{\mathbb{E}}\left\|\hat{\beta}-\beta\right\|_{2}^{2}\geq\frac{cs\log(p/s)}{n}\left(\sigma^{2}+\frac{1}{\varphi\left(CQ/p\right)}R^{2}\right),$
(16)

where

$\begin{split}\mathcal{F}_{p,n,s}(R,Q)=\bigg{\\{}&(\mathbf{V},{\mathbf{C}},\beta):{\mathbf{C}}\text{
satisfies Conditions \ref{con:C1}},\|\beta\|_{2}\leq
R,{\mathbf{C}}^{\top}\beta=0,\\\ &\left|\log\nu_{ij}-\log\nu_{kl}\right|\leq
a\text{ for constant }a>0\text{ for all }1\leq i,k\leq n,1\leq j,l\leq p,\\\
&\bar{{\mathbf{V}}}\text{ satisfies RIP condition with
}\delta_{2s}(\bar{{\mathbf{V}}})<\frac{1}{20},\\\
&e^{-\frac{3}{2}}Q\leq\bar{\nu}\leq e^{\frac{3}{2}}Q\bigg{\\}}.\end{split}$
(17)

By comparing the upper and lower bounds in Theorems 4 and 5, we can see the
proposed method achieves optimal rates when $\phi(\cdot)$ and $\pi(\cdot)$ are
of the same order, which can be guaranteed when $\mathbb{P}_{\nu_{\ast}}$ is a
distribution class with good properties, e.g., the aforementioned
Normal($\nu,\gamma$), Gamma($\nu$, $\gamma$), or Poisson($\nu$).

##  5 Numerical Studies

###  5.1 Regression Analysis for Longitudinal Microbiome Studies

In this section, we apply the proposed procedure to a longitudinal microbiome
study reported by Flores et al. (2014). In this study, 3655 microbiome samples
were taken repeatedly from gut (feces) and several other body sites of 85
college-age adults in a range of three months. Other clinical covariates
including body mass index (BMI), antibiotic disturbance and medication use
were also documented. The microbiome samples were then processed using 16S
ribosomal RNA sequencing, characterized into Operational Taxonomic Units
(OTUs) using the QIIME pipeline, and combined into taxonomic groups at
different taxonomic levels from phyla to species.

We focus on the association between BMI and gut microbiome composition at the
genus level for healthy adults by excluding subjects that have missing BMI,
antibiotic disturbance, or other medication use. For the 352 samples from the
remaining 40 subjects, 92 bacteria genera appear in more than 10% of the
samples and will be used for the analysis forward. For each subject, we
include four samples that are most similar to each other based on Bray-Curtis
distance. We consider these 4 samples to be more representative of their
owners’ true gut microbiome composition since the samples of some subjects
varies significantly across time due to possibility of sample contamination or
other unobserved factors. As a result, we have 160 samples from 40 subjects,
each subject with 4 samples, and 92 bacteria genera for the regression
analysis.

We implement the proposed variable correction regularized estimator (VC) based
on the high-dimensional log-error-in-variable regression model. Specifically,
we assume the four samples of the same subject to share the same unobserved
composition $X_{ij}$ and overdispersion paramter $\alpha_{i}$ and estimate
$\alpha_{i}$ for each subject respectively using the method of moment
estimator ${\alpha}_{i,MOM}$ described in Section A in the supplementary
materials. Then we apply the regression model (18) with $y$ representing BMI,
$W_{ij}$ representing read count of the $i$th sample and $j$th genus:

$\hat{\beta}=\mathop{\rm
arg\min}_{\beta:1_{p}^{\top}\beta=0}\left\\{\dfrac{1}{2n}\sum_{i=1}^{n}\left(y_{i}-\sum_{j=1}^{p}\log\left(W_{ij}+z_{i}\right)\beta_{j}\right)^{2}+\lambda||\beta||_{1}\right\\},$
(18) $\text{where}\quad
z_{i}=\dfrac{N_{i}+\hat{\alpha}_{i,MOM}+1}{2(\hat{\alpha}_{i,MOM}+1)}.$ (19)

For comparison, we also perform the classic zero-replacement method (ZR) in
literature with zero counts changed to $c=0.1$ and $0.5$, respectively:

$\hat{\beta}^{\rm zr}=\mathop{\rm
arg\min}_{\beta:1_{p}^{\top}\beta=0}\left\\{\dfrac{1}{2n}\sum_{i=1}^{n}\left(y_{i}-\sum_{j=1}^{p}\log\left(W_{ij}\vee
c\right)\beta_{j}\right)^{2}+\lambda\|\beta\|_{1}\right\\}.$ (20)

To obtain stable variable selection, we generate 100 bootstrap samples of size
$n/2$, repeat all methods with five-fold cross-validation choosing the tuning
parameter $\lambda$ on each subsample, and record the frequency of each
variable being selected among the 100 bootstrap fittings. We consider a
variable to be selected if its selection frequency is no less than 0.6.

Figure 2 illustrates the selection frequency of 92 genera for each method. The
dashed line corresponds to the selection frequency of 0.6. It can be observed
that the variable correction regularized estimator (VC) selects variables with
either very high or very low frequency, while zero-replacement estimator (ZR)
has much more variables selected with a mid-range frequency. This comparison
indicates that VC has much better stability in variable selection than ZR. The
variables selected by variable correction are: Actinomyces(-), Akkermansia(-),
Bacteroides(-), Coprococcus(+), Desulfovibrio(-), Dialister(+),
Haemophilus(+), Megamonas(+), Porphyromonas(-), Prevotella(+), Sutterella(+),
Veillonella(-), where (+) and (-) are signs of the fitted regression
coefficients. The variables selected by ZR with $c=0.5$ are:
Arcanobacterium(-), Bacteroides(-), Dialister(+), Haemophilus(+),
Prevotella(+), Slackia(+). Here, “(+)” and “(-)” represent positive and
negative effects, respectively. No variable is selected by zero-replacement
with $c=0.1$. These selected genera correspond to the bars exceeding the
dashed line in the left panel (VC) and middle panel (ZR with $c=0.5$) of
Figure 2 from left to right. Among the genera selected by VC but missed by ZR
with $c=0.5$, Akkermansia has been reported to be negatively related with
obesity in extensive literature (Dao et al., 2015; Derrien et al., 2017;
Karlsson et al., 2012; Zhang et al., 2009) and type II diabetes (Qin et al.,
2012). Coprococcus has also been broadly reported to be positively related
with obesity (Kasai et al., 2015) and negatively related with weight loss
induced by diet or gastric bypass surgery, as indicated by Damms-Machado et
al. (2015); Di Luccia et al. (2015); Graessler et al. (2013).

Figure 3 offers a closer look at the selection frequency with regard to the
proportion of $W_{ij}=0$ for each variable. Compared to ZR, the proposed VC
has more preference towards variables with fewer zeros, although the
preference is not decisive. This contrast of behavior between CV and ZR makes
VC more desirable since the bacteria with large proportions of zeros are often
far less reliable for prediction and interpretation purposes since they are
possessed by only a few subjects and are more difficult to generalize to a
larger population. Figure 4 compares the prediction performance of VC and ZR
with $c=0.5$, where the predicted BMI for each sample is obtained using
refitted coefficients of the genera that have selection frequency no less than
0.6. Since each subject has four samples, we also provide the averaged
predicted BMI of each subject in Figure 4. The $R^{2}$ in Figure 4 is computed
using the individual sample points. The $R^{2}$ using averaged predicted BMI
is 0.63 for VC and 0.49 for ZR with $c=0.5$. Here, the proposed VC achieves
much better prediction compared with ZR using both individual predicted BMI
and average predicted BMI.

Figure 2: Selection frequency of 92 genera using different methods Figure 3:
Comparison of selection frequency with regard to proportion of zero counts
Figure 4: Comparison of prediction performance

###  5.2 Simulation Studies

Next, we evaluate the numerical performance of the proposed procedure on
synthetic datasets. To simulate the count matrix ${\mathbf{W}}$ with $n$
samples and $p$ covariates, we first generate $N_{i}$ from negative binomial
distribution with mean $3\times 10^{4}$ and variance $3\times 10^{6}$. Here,
the main purpose of choosing negative binomial instead of the Poisson in
previous theoretical analysis is to show that the Poisson assumption on
$N_{i}$ is not crucial in real practice. Then we set
$X_{ij}=X_{i+n/2,j}=\frac{\exp(\Phi_{ij})}{\sum_{k=1}^{p}\exp(\Phi_{ij})}$ for
$j=1,\dots,p,i=1,\dots,n/2$ with $\Phi_{ij}$s generated independently from
$N(\mu_{j},1.5^{2})$, where $\mu_{1},\dots,\mu_{3}$ are drawn from
Uniform[1,3], $\mu_{4},\dots,\mu_{7}$ are drawn from Uniform[2,4], and
$\mu_{j}$ for $j=8,\dots,p$ are drawn from Uniform[0,2]. With this setting,
the average count of covariates will have a reasonable variation, with causal
covariates slightly more abundant than non-causal ones. Then we generate
$W_{ij}$ from Dirichlet-Multinomial($N_{i},\alpha X_{i1},\dots,\alpha
X_{ip}$), where the overdispersion parameter $\alpha=200,1000,5000$. The $i$th
and $(i+n/2)$th samples are designed to be from the same subject so they share
the same $X_{ij}$ and can be used to estimate their shared overdispersion
parameter. The response $y$ is generated as
$y_{i}=\sum_{j=1}^{p}\log(X_{ij})\beta_{j}+\varepsilon_{i}$, where
$\beta=(1,-0.8,-1.5,0.6,-0.9,1.2,0.4,0,\dots,0)$ is the deterministic
coefficient vector and $\varepsilon_{i}$ are i.i.d. noise generated from
$N(0,0.5^{2})$. We perform simulation study with $n=50,100$ and
$p=100,200,400$.

Using the simulated data, we evaluate the performance of the proposed variable
correction regularized estimator (VC) and the classic zero-replacement with
$c=0.5$ (ZR0.5) in estimation and prediction errors. The results are
aggregated in Figure 5. We can see VC significantly outperforms ZR0.5 in all
parameter configurations.

To evaluate the performance of the proposed method when the response variable
$y$ is shared by samples from the same subject like what we have in the real
data analysis, we repeat the aforementioned simulation with one change: $y$ is
generated with $\varepsilon_{i}=\varepsilon_{i+n/2}$ and
$\varepsilon_{i},i=1,\dots,n/2$ are i.i.d. from $N(0,0.5^{2})$. The results
are summarized in Figure 6. We can see the pattern of performance is similar
to that of Figure 5 and VC still significantly achieves smaller estimation and
prediction errors.

(a) Estimation Error

(b) Prediction Error

Figure 5: Comparison between variable correction regularized estimator (VC)
and zero-replacement estimator with $c=0.5$ (ZR0.5) in simulation analysis.
Noises $\varepsilon_{i}$ are all independent.

(a) Estimation Error

(b) Prediction Error

Figure 6: Comparison between variable correction regularized estimator (VC)
and zero-replacement estimator with $c=0.5$ (ZR0.5) in simulation analysis.
$y$ is shared by samples from the same subject.

##  6 Discussions

Understanding the association between microbial taxa and phenotypes is a
crucial problem in microbiome study. To overcome the difficulties of
randomness in covariates and zero-replacement in classic methods for
compositional data regression, we introduce the novel log-error-in-variable
regression model and the variable correction regularized estimator in this
article. The proposed estimator is surprisingly simple, works on possible
overdispersed data, and avoids any subjective zero replacement. When the count
observations follow multinomial distribution without overdispersion, we
propose to add a half to all counts, then performing the constrained Lasso
estimator; when significant overdispersion exists, we propose to add a larger
and reasonable value to all counts. We show that the proposed methods achieve
minimax optimal performance in a general class of settings.

In addition to the aforementioned microbiome study, the proposed framework can
be used to other applications on regression with count covariates. For
example, in _single-cell RNA-seq data analysis_ , the high-throughput
sequencing technique was performed on each of the single cells and the gene
expressions can be measured as the total number of reads mapped to exonic
regions. This results count matrix with rows and columns representing single
cells and gene expressions, respectively. With the proposed method, we can
perform regression analysis to study the association among the gene
expressions of single cells and clinical phenotypes. Another potential
application is in text mining, where one central task of topic modeling is to
learn the topics of various documents when they share the same vocabulary of
words. By counting the number of words or $n$-grams for each document, one can
obtain count matrix data. Compared to the absolute counts of these words and
$n$-grams, the relative abundances may be more predictive on the topic. Thus,
our proposed method is potentially useful for building classifiers for topics
of documents.

##  7 Proofs

We collect the proofs of the main results in this section. For convenience,
denote $\nu_{\min}=\min_{k,i}\nu_{ki},\nu_{\max}=\max_{k,i}\nu_{ki}$,
$\phi_{1}\left(W\right)=\log\left(W+\frac{1}{2}\right),\phi_{2}\left(W\right)=\log^{2}\left(W+\frac{1}{2}\right)$
and

${\mathbf{A}}_{W}\in\mathbb{R}^{p\times
p},\quad{\mathbf{A}}_{W}={\mathbf{B}}_{W}^{\top}{\mathbf{B}}_{W},\quad{\mathbf{P}}={\mathbf{I}}_{p}-\frac{1}{p}1_{p}1_{p}^{\top}.$

In the log-error-in-variable regression case (Theorems 1, 2, and 3), denote

$\quad\bar{{\mathbf{A}}}_{W}={\mathbf{P}}{\mathbf{A}}_{W}{\mathbf{P}},\quad\bar{{\mathbf{B}}}_{W}={\mathbf{B}}_{W}{\mathbf{P}}.$

In the general case (Theorems 4 and 5), denote

$\bar{{\mathbf{A}}}_{W}=\left({\mathbf{I}}_{p}-{\mathbf{P}}_{{\mathbf{C}}}\right){\mathbf{A}}_{W}\left({\mathbf{I}}_{p}-{\mathbf{P}}_{{\mathbf{C}}}\right),\quad\bar{{\mathbf{B}}}_{W}={\mathbf{B}}_{W}\left({\mathbf{I}}_{p}-{\mathbf{P}}_{{\mathbf{C}}}\right).$

We also denote the Orlicz-$\psi_{1}$ and $-\psi_{2}$ norms as
$\|X\|_{\psi_{1}}=\sup_{p\geq
1}p^{-1}(\mathbb{E}|X|^{p})^{1/p},\|X\|_{\psi_{2}}=\sup_{p\geq
1}p^{-1/2}(\mathbb{E}|X|^{p})^{1/p}$ for any random variable $X$.

###  7.1 Proof of Theorem 1

The model can be summarized as follows,

$y_{i}=\sum_{j=1}^{n}\log(X_{ij})\beta_{j}^{*}+\varepsilon_{i},\quad
i=1,\ldots,n;$ (21) $W_{ij}\sim{\rm Poisson}(\nu_{ij}),\quad 1\leq i\leq
n;1\leq j\leq p.$ (22)

Here $y_{i}$ and $W_{ij}$ are observable and
$\beta_{i}^{\ast},\nu_{ij}=\nu_{i}X_{ij}$ are hidden parameters. The proposed
estimator is as follows,

$\hat{\beta}=\mathop{\rm
arg\min}_{\gamma}\left(\frac{1}{2n}\|y-{\mathbf{B}}_{W}{\mathbf{P}}\gamma\|_{2}^{2}+\lambda\|\gamma\|_{1}\right)\quad\text{subject
to}\quad\quad\sum_{j=1}^{p}\gamma_{j}=0.$ (23)

In the following lemmas, suppose $W\sim\text{Poisson}(v)$ and
$W^{\prime}=W1_{\left\\{\frac{1}{10}v\leq W\leq
10v\right\\}}+v1_{\left\\{W\notin\left[\frac{1}{10}v,10v\right]\right\\}}$.

######  Lemma 1 (Bias of $\log(W^{\prime}+\frac{1}{2})$).

For any $\epsilon>0$, there exists $C_{\epsilon}>0$ that only depends on
$\epsilon$, such that for all $v\geq C_{\epsilon}$, we have
$|\mathbb{E}\log(W^{\prime}+\frac{1}{2})-\log v|\leq
4v^{-\frac{3}{2}+\epsilon}$.

######  Lemma 2 (Bias of $\log^{2}(W^{\prime}+\frac{1}{2})$).

There exists constant $C>0$ such that if $v\geq C$, then
$\left|\mathbb{E}(\log^{2}(W^{\prime}+\frac{1}{2}))-\log^{2}v\right|\leq\frac{4}{v}$.

######  Lemma 3 (Sub-Gaussianity).

There exist positive constants $K_{0}$ and $C$, such that for $v\geq C$, we
have

$\left\|\sqrt{v}\left(\log(W^{\prime}+\frac{1}{2})-{\mathbb{E}}\log(W^{\prime}+\frac{1}{2})\right)\right\|_{\psi_{2}}\leq
K_{0}.$

######  Lemma 4 (Infinity norm bound).

Under the setting of Theorem 1, there exist two constants $C$ and $C^{\prime}$
such that

$\mathbb{P}\left(\left\|\bar{{\mathbf{A}}}_{W}\beta^{*}-\bar{{\mathbf{B}}}_{W}^{\top}y\right\|_{\infty}\leq
C\sqrt{n\log
p\left(\sigma^{2}+\frac{p}{\bar{\nu}}\|\beta^{*}\|_{2}^{2}\right)}\right)\geq
1-4p^{-C^{\prime}}.$

Now let us move to the proof of Theorem 1. Denote $h=\hat{\beta}-\beta^{*}$.
By lemma 4, with probability at least $1-4p^{-C^{\prime}}$,

$\left\|\bar{{\mathbf{A}}}_{W}\beta^{*}-\bar{{\mathbf{B}}}_{W}^{\top}y\right\|_{\infty}\leq
C\sqrt{n\log
p\left(\sigma^{2}+\frac{p}{\bar{\nu}}\|\beta^{*}\|_{2}^{2}\right)}=\frac{n}{2}\lambda.$
(24)

By the definition of the estimator and the fact that
${\mathbf{P}}\beta^{*}=\beta^{*},{\mathbf{P}}\hat{\beta}=\hat{\beta}$, we have

$\frac{1}{2n}\left\|y-\bar{{\mathbf{B}}}_{W}\hat{\beta}\right\|_{2}^{2}+\lambda\left\|\hat{\beta}\right\|_{1}\leq\frac{1}{2n}\left\|y-\bar{{\mathbf{B}}}_{W}\beta^{*}\right\|_{2}^{2}+\lambda\left\|\beta^{*}\right\|_{1},$

which means

$\frac{1}{2n}\left(2h^{\top}\left(\bar{{\mathbf{A}}}_{W}\beta^{*}-\bar{{\mathbf{B}}}_{W}^{\top}y\right)+h^{\top}\bar{{\mathbf{A}}}_{W}h\right)\leq\lambda\left(\|\beta^{*}\|_{1}-\|\hat{\beta}\|_{1}\right).$

Denote $S=\text{supp}\left(\beta^{\ast}\right)$. Noting that

$\begin{split}\|\beta^{*}\|_{1}-\|\hat{\beta}\|_{1}=\|\beta_{S}^{\ast}\|_{1}-\|\hat{\beta}_{S}\|_{1}-\|\hat{\beta}_{S^{c}}\|_{1}\leq\|\beta_{S}^{\ast}-\hat{\beta}_{S}\|_{1}-\|h_{S^{c}}\|_{1}\leq\|h_{S}\|_{1}-\|h_{S^{c}}\|_{1},\end{split}$

we have

$\frac{1}{2n}\left(2h^{\top}\left(\bar{{\mathbf{A}}}_{W}\beta^{*}-\bar{{\mathbf{B}}}_{W}^{\top}y\right)+h^{\top}\bar{{\mathbf{A}}}_{W}h\right)\leq\lambda\left(\|h_{S}\|_{1}-\|h_{S^{c}}\|_{1}\right).$
(25)

In addition,

$\begin{split}&\frac{1}{2n}\left(2h^{\top}\left(\bar{{\mathbf{A}}}_{W}\beta^{*}-\bar{{\mathbf{B}}}_{W}^{\top}y\right)+h^{\top}\bar{{\mathbf{A}}}_{W}h\right)=\frac{1}{2n}\left(2h^{\top}\left(\bar{{\mathbf{A}}}_{W}\beta^{*}-\bar{{\mathbf{B}}}_{W}^{\top}y\right)+\|\bar{{\mathbf{B}}}_{W}h\|_{2}^{2}\right)\\\
\geq&\frac{1}{n}h^{\top}\left(\bar{{\mathbf{A}}}_{W}\beta^{*}-\bar{{\mathbf{B}}}_{W}^{\top}y\right)\geq-\frac{1}{n}\|h\|_{1}\|\bar{{\mathbf{A}}}_{W}\beta^{*}-\bar{{\mathbf{B}}}_{W}^{\top}y\|_{\infty}\\\
=&-\frac{1}{n}\left(\|h_{S}\|_{1}+\|h_{S^{c}}\|_{1}\right)\|\bar{{\mathbf{A}}}_{W}\beta^{*}-\bar{{\mathbf{B}}}_{W}^{\top}y\|_{\infty}.\end{split}$

If (24) holds, then (24), (25) and the previous inequality together imply

$-\frac{1}{2}\left(\|h_{S}\|_{1}+\|h_{S^{c}}\|_{1}\right)\leq\|h_{S}\|_{1}-\|h_{S^{c}}\|_{1},\quad\text{i.e.,}\quad\|h_{S}\|_{1}\geq\frac{1}{3}\|h_{S^{c}}\|_{1}.$

Therefore,

$\|h_{\max(s)}\|_{1}\geq\|h_{S}\|_{1}\geq\frac{1}{3}\|h_{S^{c}}\|_{1}\geq\frac{1}{3}\|h_{-\max(s)}\|_{1},$

where we set $h_{\max(s)}$ as $h$ with all but the largest $s$ entries in
absolute value set to zero, and $h_{-\max(s)}=h-h_{\max(s)}$. By the KKT
condition on (23), we have

$\left\|\bar{{\mathbf{B}}}_{W}^{\top}\left(\bar{{\mathbf{B}}}_{W}\hat{\beta}-y\right)+1_{p}\kappa\right\|_{\infty}\leq
n\lambda$

for some $\kappa\in\mathbb{R}$. Since
$\|{\mathbf{P}}x\|_{\infty}=\|x-\frac{1}{p}\sum_{i=1}^{p}x_{i}\|_{\infty}\leq\|x\|_{\infty}+\left|\frac{1}{p}\sum_{i=1}^{p}x_{i}\right|\leq
2\|x\|_{\infty}$,

$\begin{split}\left\|\bar{{\mathbf{A}}}_{W}\hat{\beta}-\bar{{\mathbf{B}}}_{W}^{\top}y\right\|_{\infty}=\left\|{\mathbf{P}}\left(\bar{{\mathbf{B}}}_{W}^{\top}\left(\bar{{\mathbf{B}}}_{W}\hat{\beta}-y\right)+1_{p}\kappa\right)\right\|_{\infty}\leq
2\left\|\bar{{\mathbf{B}}}_{W}^{\top}\left(\bar{{\mathbf{B}}}_{W}\hat{\beta}-y\right)+1_{p}\kappa\right\|_{\infty}\leq
2n\lambda.\end{split}$

By (24) and the previous inequality, with probability at least
$1-4p^{-C^{\prime}}$,

$\left\|\bar{{\mathbf{A}}}_{W}h\right\|_{\infty}=\left\|\bar{{\mathbf{A}}}_{W}\left(\beta^{\ast}-\hat{\beta}\right)\right\|_{\infty}\leq\|\bar{{\mathbf{A}}}_{W}\hat{\beta}-\bar{{\mathbf{B}}}_{W}^{\top}y\|_{\infty}+\left\|\bar{{\mathbf{A}}}_{W}\beta^{*}-\bar{{\mathbf{B}}}_{W}^{\top}y\right\|_{\infty}\leq\frac{5}{2}n\lambda.$
(26)

Therefore, with probability at least $1-4p^{-C^{\prime}}$,

$h_{\max(s)}\bar{{\mathbf{A}}}_{W}h\leq\left\|h_{\max(s)}\right\|_{1}\left\|\bar{{\mathbf{A}}}_{W}h\right\|_{\infty}\\\
\leq\frac{5}{2}n\lambda\cdot\sqrt{s}\left\|h_{\max(s)}\right\|_{2}.$ (27)

Define $\alpha=\left\|h_{\max(s)}\right\|_{1}/s$, then

$\begin{split}\left\|h_{-max(s)}\right\|_{\infty}&\leq\alpha\leq
3\alpha,\quad\left\|h_{-max(s)}\right\|_{1}\leq
3\left\|h_{\max(s)}\right\|_{1}=3s\alpha.\end{split}$

Apply Lemma 1.1 in Cai and Zhang (2014), $h_{-\max(s)}$ can be expressed as a
convex combinations of sparse vectors:
$h_{-\max(s)}=\sum_{i=1}^{M}\lambda_{i}u_{i}$, where $u_{i}$ is $s$-sparse and

$\left\|u_{i}\right\|_{1}=\left\|h^{(2)}\right\|_{1},\quad\left\|u_{i}\right\|_{\infty}\leq
3\alpha,\quad supp(u_{i})\subseteq supp(h_{-\max(s)}).$

Thus

$\begin{split}\left\|u_{i}\right\|_{2}&\leq\sqrt{\left\|u_{i}\right\|_{0}}\left\|u_{i}\right\|_{\infty}\leq\sqrt{s}\cdot
3\alpha=3\sqrt{s}\alpha.\end{split}$ (28)

Suppose $0\leq\mu\leq 1$ is to be determined. Denote
$\gamma_{i}=h_{\max(s)}+\mu u_{i}$, then

$\begin{split}\sum_{j=1}^{M}\lambda_{j}\gamma_{j}-\frac{1}{2}\gamma_{i}&=h_{\max(s)}+\mu
h_{-\max(s)}-\frac{1}{2}\gamma_{i}=(1-\mu-\frac{1}{2})h_{\max(s)}-\frac{1}{2}\mu
u_{i}+\mu h.\end{split}$

Since $h_{\max(s)}$ and $u_{i}$ are $s$-sparse vectors,
$\gamma_{i}=h_{\max(s)}+\mu
u_{i},\sum_{j=1}^{M}\lambda_{j}\gamma_{j}-\frac{1}{2}\gamma_{i}-\mu
h=(\frac{1}{2}-\mu)h_{\max(s)}-\frac{1}{2}\mu u_{i}$ are all $2s$-sparse
vectors.
Suppose $x=\left\|h_{\max(s)}\right\|_{2}$, by (28),

$\begin{split}\left\|u_{i}\right\|_{2}\leq 3\sqrt{s}\alpha\leq
3\left\|h_{\max(s)}\right\|_{2}=3x.\end{split}$ (29)

Also, we can check that

$\begin{split}&\sum_{i=1}^{M}\lambda_{i}\left\|\bar{{\mathbf{B}}}_{W}\left(\sum_{j=1}^{M}\lambda_{j}\gamma_{j}-\frac{1}{2}\gamma_{i}\right)\right\|_{2}^{2}=\sum_{i=1}^{M}\lambda_{i}\cdot\frac{1}{4}\|\bar{{\mathbf{B}}}_{W}\gamma_{i}\|_{2}^{2}.\end{split}$
$\begin{split}0=&\sum_{i=1}^{M}\left\|\bar{{\mathbf{B}}}_{W}\cdot\left(h_{\max(s)}+\mu
h_{-\max(s)}-\frac{1}{2}(h_{\max(s)}+\mu
u_{i})\right)\right\|_{2}^{2}-\sum_{i=1}^{M}\frac{\lambda_{i}}{4}\|\bar{{\mathbf{B}}}_{W}\gamma_{i}\|_{2}^{2}\\\
=&\sum_{i=1}^{M}\left\|\bar{{\mathbf{B}}}_{W}\left((\frac{1}{2}-\mu)h_{\max(s)}-\frac{\mu}{2}u_{i}+\mu
h\right)\right\|_{2}^{2}-\sum_{i=1}^{M}\frac{\lambda_{i}}{4}\|\bar{{\mathbf{B}}}_{W}\gamma_{i}\|_{2}^{2}\\\
=&\sum_{i=1}^{M}\lambda_{i}\left\|\bar{{\mathbf{B}}}_{W}\left((\frac{1}{2}-\mu)h_{\max(s)}-\frac{\mu}{2}u_{i}\right)\right\|_{2}^{2}+2\mu\left((\frac{1}{2}-\mu)h_{\max(s)}-\frac{\mu}{2}h_{-\max(s)}\right)^{\top}\bar{{\mathbf{A}}}_{W}h\\\
&+\mu^{2}\|\bar{{\mathbf{B}}}_{W}h\|_{2}^{2}-\sum_{i=1}^{M}\frac{\lambda_{i}}{4}\|\bar{{\mathbf{B}}}_{W}\gamma_{i}\|_{2}^{2}\\\
=&\sum_{i=1}^{M}\left\|\bar{{\mathbf{B}}}_{W}\left((\frac{1}{2}-\mu)h_{\max(s)}-\frac{\mu}{2}u_{i}\right)\right\|_{2}^{2}+\mu(1-\mu)h_{\max(s)}^{\top}\bar{{\mathbf{A}}}_{W}h-\sum_{i=1}^{M}\frac{\lambda_{i}}{4}\|\bar{{\mathbf{B}}}_{W}\gamma_{i}\|_{2}^{2}.\end{split}$
(30)

In the third equation above, we used the facts that
$h_{-\max(s)}=\sum_{i=1}^{M}\lambda_{i}u_{i}$ and
$\bar{{\mathbf{A}}}_{W}=\bar{{\mathbf{B}}}_{W}^{\top}\bar{{\mathbf{B}}}_{W}$.
Noting that supp$(h_{\max(s)})$ $\cap$ supp$(u_{i})$ $(i=1,\dots,M)$ are empty
sets and $\gamma_{i},(\frac{1}{2}-\mu)h_{\max(s)}-\frac{\mu}{2}u_{i}$ are all
$2s$-sparse vectors, we can apply Condition 1, (27), and (29) to (30) and
obtain

$\begin{split}0\leq&n(1+\delta_{2s}\left(\bar{{\mathbf{B}}}_{W}\right))\sum_{i=1}^{M}\lambda_{i}\left((\frac{1}{2}-\mu)^{2}\|h_{\max(s)}\|_{2}^{2}+\frac{\mu^{2}}{4}\|u_{i}\|_{2}^{2}\right)+\frac{5}{2}n\mu(1-\mu)\sqrt{s}\lambda\|h_{\max(s)}\|_{2}\\\
&-n(1-\delta_{2s}\left(\bar{{\mathbf{B}}}_{W}\right))\sum_{i=1}^{M}\frac{\lambda_{i}}{4}\left(\|h_{\max(s)}\|_{2}^{2}+\mu^{2}\|u_{i}\|_{2}^{2}\right)\\\
=&n\sum_{i=1}^{M}\lambda_{i}\left\\{\left((1+\delta_{2s}\left(\bar{{\mathbf{B}}}_{W}\right))(\frac{1}{2}-\mu)^{2}-(1-\delta_{2s}\left(\bar{{\mathbf{B}}}_{W}\right))\cdot\frac{1}{4}\right)\|h_{\max(s)}\|_{2}^{2}+\frac{1}{2}\delta_{2s}\left(\bar{{\mathbf{B}}}_{W}\right)\mu^{2}\|u_{i}\|_{2}^{2}\right\\}\\\
&+\frac{5}{2}n\mu(1-\mu)\sqrt{s}\lambda\|h_{\max(s)}\|_{2}\\\
\leq&n\sum_{i=1}^{M}\lambda_{i}\left\\{\left((1+\delta_{2s}\left(\bar{{\mathbf{B}}}_{W}\right))(\frac{1}{2}-\mu)^{2}-(1-\delta_{2s}\left(\bar{{\mathbf{B}}}_{W}\right))\cdot\frac{1}{4}\right)x^{2}+\frac{1}{2}\delta_{2s}\left(\bar{{\mathbf{B}}}_{W}\right)\mu^{2}\cdot
9x^{2}\right\\}\\\ &+\frac{5}{2}n\mu(1-\mu)\sqrt{s}\lambda x\\\
\leq&n\left[(\mu^{2}-\mu)+\frac{1}{2}\delta_{2s}\left(\bar{{\mathbf{B}}}_{W}\right)\left(1-2\mu+11\mu^{2}\right)\right]x^{2}+\frac{5}{2}n\mu(1-\mu)\sqrt{s}\lambda
x\end{split}$

holds with probability at least $1-4p^{-C^{\prime}}-\epsilon^{\prime}$. Set
$\mu=\frac{1}{2}$ and notice that
$\delta_{2s}\left(\bar{{\mathbf{B}}}_{W}\right)\leq\frac{1}{10}$, we know that
with probability at least $1-4p^{-C^{\prime}}-\epsilon^{\prime}$,

$x\leq\left(1-\frac{11}{2}\delta_{2s}\left(\bar{{\mathbf{B}}}_{W}\right)\right)^{-1}\frac{5}{2}\sqrt{s}\lambda\leq\frac{50}{9}\sqrt{s}\lambda.$

Therefore, with probability at least $1-4p^{-C^{\prime}}-\epsilon^{\prime}$,

$\begin{split}\|h\|_{2}=\sqrt{\|h_{\max(s)}\|_{2}^{2}+\|h_{-\max(s)}\|_{2}^{2}}\leq\sqrt{\|h_{\max(s)}\|_{2}^{2}+9\|h_{\max(s)}\|_{2}^{2}}=\sqrt{10}x\leq
18\sqrt{s}\lambda,\end{split}$

which means with probability at least $1-4p^{-C^{\prime}}-\epsilon^{\prime}$,

$\|\hat{\beta}-\beta^{*}\|_{2}^{2}\leq 324s\lambda^{2}=\frac{Cs\log
p}{n}\left(\sigma^{2}+\frac{p}{\bar{\nu}}\|\beta^{*}\|_{2}^{2}\right).$

∎

###  7.2 Proof of Theorem 3

$W_{ij}$ can be seen as $W_{ij}|q_{i}\sim\text{Poisson}(\nu_{i}q_{ij})$ and
$q_{ij}\sim\text{Beta}(\alpha_{i}X_{ij},\alpha_{i}-\alpha_{i}X_{ij})$. We
introduce the truncated version of $W_{ij}$ as

$W^{\prime}_{ij}|W_{ij},q_{i}=\left\\{\begin{array}[]{l}W_{ij},\quad
q_{ij}\in\left[\frac{X_{ij}}{10},10X_{ij}\right]\text{and
}W\in[\frac{\nu_{i}q_{ij}}{10},10\nu_{i}q_{ij}],\\\
\nu_{ij},\quad\text{otherwise}.\end{array}\right.$

Under the setting of Theorem 3, we introduce the following lemmas.

######  Lemma 5.

$\left|{\mathbb{E}}\log\left(W_{ij}^{\prime}+\frac{N_{i}+\lambda_{i}+1}{2(1+\lambda_{i})}\right)-\log\left(\nu_{ij}\right)\right|\leq
10\frac{\zeta_{i}}{\nu_{ij}}$ holds if $\nu_{ij}\geq C(\delta)$.

######  Lemma 6.

$\left|{\mathbb{E}}\log^{2}\left(W_{ij}^{\prime}+\frac{N_{i}+\lambda_{i}+1}{2(1+\lambda_{i})}\right)-\log^{2}\left(\nu_{ij}\right)\right|\leq
C\log(\nu_{ij})\frac{\zeta_{i}}{\nu_{ij}}$ holds if $\nu_{ij}\geq C(\delta)$.

######  Lemma 7 (Sub-Gaussianity).

$\left\|\sqrt{\frac{\nu_{ij}}{\zeta_{\max}}}\left(\log\left(W^{\prime}_{ij}+c_{0}\right)-{\mathbb{E}}\log\left(W^{\prime}_{ij}+c_{0}\right)\right)\right\|_{\psi_{2}}\leq
K_{0}$ for a constant $K_{0}$.

######  Lemma 8 (Infinity norm bound).

Under the setting of Theorem 3, there exists constants $C,C^{\prime}>0$, such
that

$\mathbb{P}\left(\left\|\bar{{\mathbf{A}}}_{W}\beta^{*}-\bar{{\mathbf{B}}}_{W}^{\top}y\right\|_{\infty}\leq
C\sqrt{n\log
p}\left(\sigma+\left(\frac{p}{\bar{\nu}}\zeta_{\max}\right)^{\frac{1}{2}}\|\beta^{*}\|_{1}\right)\right)\geq
1-6p^{-C^{\prime}}.$

Applying these lemmas, we can prove Theorem 3 by essentially the same method
as the proof of Theorem 1 with $\lambda=C\sqrt{\frac{\log
p}{n}}\left(\sigma+\left(\frac{p}{\bar{\nu}}\zeta_{\max}\right)^{\frac{1}{2}}\|\beta^{*}\|_{1}\right)$.∎

## Acknowledgments

We would like to thank Sündüz Keleş for helpful discussions.

## References

  * Aitchison (2003) Aitchison, J. (2003).  The statistical analysis of compositional data.  New Jersey: Blackburn Press.
  * Aitchison and Bacon-shone (1984) Aitchison, J. and Bacon-shone, J. (1984).  Log contrast models for experiments with mixtures.  Biometrika, 71(2):323–330.
  * Belloni et al. (2017a) Belloni, A., Chernozhukov, V., and Kaul, A. (2017a).  Confidence bands for coefficients in high dimensional linear models with error-in-variables.  arXiv preprint arXiv:1703.00469.
  * Belloni et al. (2017b) Belloni, A., Chernozhukov, V., Kaul, A., Rosenbaum, M., and Tsybakov, A. B. (2017b).  Pivotal estimation via self-normalization for high-dimensional linear models with error in variables.  arXiv preprint arXiv:1708.08353.
  * Belloni et al. (2016a) Belloni, A., Rosenbaum, M., and Tsybakov, A. B. (2016a).  An {$\ell_{1},\ell_{2},\ell_{\infty}$}-regularization approach to high-dimensional errors-in-variables models.  Electronic Journal of Statistics, 10:1729–1750.
  * Belloni et al. (2016b) Belloni, A., Rosenbaum, M., and Tsybakov, A. B. (2016b).  Linear and conic programming estimators in high dimensional errors-in-variables models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), pages n/a–n/a.
  * Bickel et al. (2009) Bickel, P. J., Ritov, Y., and Tsybakov, A. B. (2009).  Simultaneous analysis of lasso and dantzig selector.  The Annals of Statistics, 37(4):1705–1732.
  * Cai and Zhang (2014) Cai, T. T. and Zhang, A. (2014).  Sparse representation of a polytope and recovery of sparse signals and low-rank matrices.  IEEE Transactions on Information Theory, 60(1):122–132.
  * Candes et al. (2007) Candes, E., Tao, T., et al. (2007).  The dantzig selector: Statistical estimation when p is much larger than n.  The Annals of Statistics, 35(6):2313–2351.
  * Candes and Tao (2005) Candes, E. J. and Tao, T. (2005).  Decoding by linear programming.  IEEE transactions on information theory, 51(12):4203–4215.
  * Cao and Xie (2016) Cao, Y. and Xie, Y. (2016).  Poisson matrix recovery and completion.  IEEE Trans. Signal Processing, 64(6):1609–1620.
  * Cao et al. (2017) Cao, Y., Zhang, A., and Li, H. (2017).  Multi-sample estimation of bacterial composition matrix in metagenomics data.  arXiv preprint arXiv:1706.02380.
  * Chen and Li (2013) Chen, J. and Li, H. (2013).  Variable selection for sparse dirichlet-multinomial regression with an application to microbiome data analysis.  The annals of applied statistics, 7(1).
  * Damms-Machado et al. (2015) Damms-Machado, A., Mitra, S., Schollenberger, A. E., Kramer, K. M., Meile, T., Königsrainer, A., Huson, D. H., and Bischoff, S. C. (2015).  Effects of surgical and dietary weight loss therapy for obesity on gut microbiota composition and nutrient absorption.  BioMed research international, 2015.
  * Dao et al. (2015) Dao, M. C., Everard, A., Aron-Wisnewsky, J., Sokolovska, N., Prifti, E., Verger, E. O., Kayser, B. D., Levenez, F., Chilloux, J., Hoyles, L., et al. (2015).  Akkermansia muciniphila and improved metabolic health during a dietary intervention in obesity: relationship with gut microbiome richness and ecology.  Gut, pages gutjnl–2014.
  * Datta and Zou (2017) Datta, A. and Zou, H. (2017).  Cocolasso for high-dimensional error-in-variables regression.  The Annals of Statistics, to appear.
  * Deming (1943) Deming, W. E. (1943).  Statistical adjustment of data.  Wiley.
  * Derrien et al. (2017) Derrien, M., Belzer, C., and de Vos, W. M. (2017).  Akkermansia muciniphila and its role in regulating host functions.  Microbial pathogenesis, 106:171–181.
  * Di Luccia et al. (2015) Di Luccia, B., Crescenzo, R., Mazzoli, A., Cigliano, L., Venditti, P., Walser, J.-C., Widmer, A., Baccigalupi, L., Ricca, E., and Iossa, S. (2015).  Rescue of fructose-induced metabolic syndrome by antibiotics or faecal transplantation in a rat model of obesity.  PLoS One, 10(8):e0134893.
  * Fan and Li (2001) Fan, J. and Li, R. (2001).  Variable selection via nonconcave penalized likelihood and its oracle properties.  Journal of the American statistical Association, 96(456):1348–1360.
  * Flores et al. (2014) Flores, G. E., Caporaso, J. G., Henley, J. B., Rideout, J. R., Domogala, D., Chase, J., Leff, J. W., Vázquez-Baeza, Y., Gonzalez, A., Knight, R., et al. (2014).  Temporal variability is a personalized feature of the human microbiome.  Genome biology, 15(12):531.
  * Graessler et al. (2013) Graessler, J., Qin, Y., Zhong, H., Zhang, J., Licinio, J., Wong, M.-L., Xu, A., Chavakis, T., Bornstein, A., Ehrhart-Bornstein, M., et al. (2013).  Metagenomic sequencing of the human gut microbiome before and after bariatric surgery in obese patients with type 2 diabetes: correlation with inflammatory and metabolic parameters.  The pharmacogenomics journal, 13(6):514.
  * Greene (2003) Greene, W. H. (2003).  Econometric analysis.  Pearson Education India.
  * Hall et al. (2016) Hall, E. C., Raskutti, G., and Willett, R. (2016).  Inference of high-dimensional autoregressive generalized linear models.  arXiv preprint arXiv:1605.02693.
  * Hausman (2001) Hausman, J. (2001).  Mismeasured variables in econometric analysis: problems from the right and problems from the left.  Journal of Economic perspectives, 15(4):57–67.
  * Jiang et al. (2015) Jiang, X., Raskutti, G., and Willett, R. (2015).  Minimax optimal rates for poisson inverse problems with physical constraints.  IEEE Transactions on Information Theory, 61(8):4458–4474.
  * Karlsson et al. (2012) Karlsson, C. L., Önnerfält, J., Xu, J., Molin, G., Ahrné, S., and Thorngren-Jerneck, K. (2012).  The microbiota of the gut in preschool children with normal and excessive body weight.  Obesity, 20(11):2257–2261.
  * Kasai et al. (2015) Kasai, C., Sugimoto, K., Moritani, I., Tanaka, J., Oya, Y., Inoue, H., Tameda, M., Shiraki, K., Ito, M., Takei, Y., et al. (2015).  Comparison of the gut microbiota composition between obese and non-obese individuals in a japanese population, as analyzed by terminal restriction fragment length polymorphism and next-generation sequencing.  BMC gastroenterology, 15(1):100.
  * La Rosa et al. (2012) La Rosa, P. S., Brooks, J. P., Deych, E., Boone, E. L., Edwards, D. J., Wang, Q., Sodergren, E., Weinstock, G., and Shannon, W. D. (2012).  Hypothesis testing and power calculations for taxonomic-based human microbiome data.  PloS one, 7(12):e52078.
  * Lewis et al. (2015) Lewis, J. D., Chen, E. Z., Baldassano, R. N., Otley, A. R., Griffiths, A. M., Lee, D., Bittinger, K., Bailey, A., Friedman, E. S., Hoffmann, C., et al. (2015).  Inflammation, antibiotics, and diet as environmental stressors of the gut microbiome in pediatric crohn’s disease.  Cell host & microbe, 18(4):489–500.
  * Li (2015) Li, H. (2015).  Microbiome, metagenomics, and high-dimensional compositional data analysis.  Annual Review of Statistics and Its Application, 2:73–94.
  * Lin et al. (2014) Lin, W., Shi, P., Feng, R., and Li, H. (2014).  Variable selection in regression with compositional covariates.  Biometrika, 101(4):785–797.
  * Loh and Wainwright (2012) Loh, P.-L. and Wainwright, M. (2012).  High-dimensional regression with noisy and missing data: Provable guarantees with nonconvexity.  The Annals of Statistics, 40(3):1637–1664.
  * Lu et al. (2018) Lu, J., Shi, P., and Li, H. (2018).  Generalized linear models with linear constraints for microbiome compositional data.  arXiv preprint arXiv:1801.03238.
  * Martin-Fernandez et al. (2000) Martin-Fernandez, J., Barceló-Vidal, C., and Pawlowsky-Glahn, V. (2000).  Zero replacement in compositional data sets.  In Data analysis, classification, and related methods, pages 155–160. Springer.
  * Martín-Fernández et al. (2003) Martín-Fernández, J. A., Barceló-Vidal, C., and Pawlowsky-Glahn, V. (2003).  Dealing with zeros and missing values in compositional data sets using nonparametric imputation.  Mathematical Geology, 35(3):253–278.
  * Martín-Fernández et al. (2015) Martín-Fernández, J.-A., Hron, K., Templ, M., Filzmoser, P., and Palarea-Albaladejo, J. (2015).  Bayesian-multiplicative treatment of count zeros in compositional data sets.  Statistical Modelling, 15(2):134–158.
  * Mosimann (1962) Mosimann, J. E. (1962).  On the compound multinomial distribution, the multivariate $\beta$-distribution, and correlations among proportions.  Biometrika, 49(1/2):65–82.
  * Pal (1980) Pal, M. (1980).  Consistent moment estimators of regression coefficients in the presence of errors in variables.  Journal of Econometrics, 14(3):349–364.
  * Purdom and Holmes (2005) Purdom, E. and Holmes, S. P. (2005).  Error distribution for gene expression data.  Statistical applications in genetics and molecular biology, 4(1).
  * Qin et al. (2012) Qin, J., Li, Y., Cai, Z., Li, S., Zhu, J., Zhang, F., Liang, S., Zhang, W., Guan, Y., Shen, D., et al. (2012).  A metagenome-wide association study of gut microbiota in type 2 diabetes.  Nature, 490(7418):55.
  * Rigollet and Tsybakov (2011) Rigollet, P. and Tsybakov, A. (2011).  Exponential screening and optimal rates of sparse estimation.  The Annals of Statistics, 39(2):731–771.
  * Rosenbaum and Tsybakov (2010) Rosenbaum, M. and Tsybakov, A. B. (2010).  Sparse recovery under matrix uncertainty.  The Annals of Statistics, 38(5):2620–2651.
  * Rosenbaum and Tsybakov (2013) Rosenbaum, M. and Tsybakov, A. B. (2013).  Improved matrix uncertainty selector.  In From Probability to Statistics and Back: High-Dimensional Models and Processes–A Festschrift in Honor of Jon A. Wellner, pages 276–290. Institute of Mathematical Statistics.
  * Rudelson et al. (2017) Rudelson, M., Zhou, S., et al. (2017).  Errors-in-variables models with dependent measurements.  Electronic Journal of Statistics, 11(1):1699–1797.
  * Schwabe and Jobin (2013) Schwabe, R. F. and Jobin, C. (2013).  The microbiome and cancer.  Nature Reviews Cancer, 13(11):800.
  * Shi et al. (2016) Shi, P., Zhang, A., and Li, H. (2016).  Regression analysis for microbiome compositional data.  The Annals of Applied Statistics, 10(2):1019–1040.
  * Slijepcevic et al. (2002) Slijepcevic, S., Megerian, S., and Potkonjak, M. (2002).  Location errors in wireless embedded sensor networks: sources, models, and effects on applications.  ACM SIGMOBILE Mobile Computing and Communications Review, 6(3):67–78.
  * The Human Microbiome Project Consortium (2012) The Human Microbiome Project Consortium (2012).  A framework for human microbiome research.  Nature, 486(7402):215.
  * Tibshirani (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society. Series B (Methodological), pages 267–288.
  * Turnbaugh et al. (2009) Turnbaugh, P. J., Hamady, M., Yatsunenko, T., Cantarel, B. L., Duncan, A., Ley, R. E., Sogin, M. L., Jones, W. J., Roe, B. A., Affourtit, J. P., et al. (2009).  A core gut microbiome in obese and lean twins.  nature, 457(7228):480–484.
  * Turnbaugh et al. (2006) Turnbaugh, P. J., Ley, R. E., Mahowald, M. A., Magrini, V., Mardis, E. R., and Gordon, J. I. (2006).  An obesity-associated gut microbiome with increased capacity for energy harvest.  nature, 444(7122):1027.
  * Tvedebrink (2010) Tvedebrink, T. (2010).  Overdispersion in allelic counts and $\theta$-correction in forensic genetics.  Theoretical population biology, 78(3):200–210.
  * Wang and Zhao (2017) Wang, T. and Zhao, H. (2017).  Structured subcomposition selection in regression and its application to microbiome data analysis.  The Annals of Applied Statistics, 11(2):771–791.
  * Weir and Hill (2002) Weir, B. S. and Hill, W. G. (2002).  Estimating F-statistics.  Annual review of genetics, 36(1):721–750.
  * Zhang et al. (2009) Zhang, H., DiBaise, J. K., Zuccolo, A., Kudrna, D., Braidotti, M., Yu, Y., Parameswaran, P., Crowell, M. D., Wing, R., Rittmann, B. E., et al. (2009).  Human gut microbiota in obesity and after gastric bypass.  Proceedings of the National Academy of Sciences, 106(7):2365–2370.