# Estimating Fold Changes from Partially Observed Outcomes with Applications
in Microbial Metagenomics

David S Clausen
and
Amy D Willis
Department of Biostatistics, University of Washington  Contact:
<EMAIL_ADDRESS>The authors gratefully acknowledge support from the National
Institute of General Medical Sciences (R35 GM133420). The authors are grateful
to Sarah Teichman for thoughtful suggestions on the manuscript and code.

###### Abstract

We consider the problem of estimating fold-changes in the expected value of a
multivariate outcome that is observed subject to unknown sample-specific and
category-specific perturbations. We are motivated by high-throughput
sequencing studies of the abundance of microbial taxa, in which microbes are
systematically over- and under-detected relative to their true abundances. Our
log-linear model admits a partially identifiable estimand, and we establish
full identifiability by imposing interpretable parameter constraints. To
reduce bias and guarantee the existence of parameter estimates in the presence
of sparse observations, we apply an asymptotically negligible and constraint-
invariant penalty to our estimating function. We develop a fast coordinate
descent algorithm for estimation, and an augmented Lagrangian algorithm for
estimation under null hypotheses. We construct a model-robust score test, and
demonstrate valid inference even for small sample sizes and violated
distributional assumptions. The flexibility of the approach and comparisons to
related methods are illustrated via a meta-analysis of microbial associations
with colorectal cancer.

##  1 Introduction

Communities of microorganisms known as microbiomes have become a focus of
study across multiple scientific disciplines, including in many areas of human
health and disease (Young, 2017; Yamashita and Takeshita, 2017; Shreiner et
al., 2015; Moffatt and Cookson, 2017). Modern microbiome research leans
heavily on high-throughput sequencing techniques, which enable researchers to
identify large numbers of microbial taxa on the basis of genetic material
detected in a given environment. In this paper, we consider the problem of
making statistical statements about changes in the “absolute” concentrations
of microbes in an environment using only the “relative” information provided
by marker gene and whole genome shotgun sequencing. Our approach also accounts
for the observation that different microbes are consistently over- or under-
counted (McLaren et al., 2019).

Most generally, this paper introduces a method for inference on ratios of
means of a nonnegative multivariate outcome observed only up to unknown
sample-specific and category-specific scaling terms. To our knowledge no
similar methods currently exist, but a related niche in microbiome science is
occupied by a range of regression models that rely on log-ratio
transformations to target a similar estimand. These include ANCOM and its
variants, which attempt to distinguish sampling, structural, and “outlier”
zero values in high-throughput sequencing data prior to transformation and
modeling (Mandal et al., 2015; Kaul et al., 2017; Lin and Peddada, 2020); and
ALDEx2, which employs a Bayesian approach to impute nonzero values before
transformation (Fernandes et al., 2014).

ANCOM and ALDEx2 fall into the larger category of “differential abundance”
methods, which include LEfSe, DESeq2, edgeR, metagenomeSeq, IFAA, and corncob
(Segata et al., 2011; Love et al., 2014; Robinson et al., 2010; Paulson et
al., 2013; Li et al., 2021; Martin et al., 2020), as well as many others (see
Nearing et al. (2022) for a recent empirical comparison). DESeq2 and edgeR
were originally developed for RNAseq data and attempt to address excess-
Poisson dispersion via a negative binomial likelihood (Love et al., 2014;
Robinson et al., 2010). By contrast, metagenomeSeq was developed specifically
for microbiome analysis and employs a zero-inflated Gaussian model motivated
by the high degree of sparsity in microbiome datasets (Paulson et al., 2013).
LEfSe, also developed for microbiome analysis, attempts to identify
“metagenomic biomarkers” via a Kruskal-Wallis test followed by Wilcoxon tests
on normalized sequencing data (Segata et al., 2011), and IFAA fits a log-
linear zero-inflated Gaussian model by identifying a reference taxon whose
abundance is conditionally independent of covariates (Li et al., 2021).
Corncob (Martin et al., 2020) models both taxon relative abundances and
dispersion of relative abundances via a beta-binomial model with logit link
functions. We note that although the regression models above are all referred
to as differential abundance methods (Nearing et al., 2022), they do not
target a common estimand and hence do not produce strictly comparable
inference. Furthermore, most of these methods target estimands defined in
terms of microbial proportions, while our approach targets an estimand defined
on the total microbial abundance scale (e.g., taxon-specific cell counts, DNA
concentration, or copy number).

Our manuscript is structured as follows: we introduce our model and establish
identifiability of equivalence classes of parameters in Section 2. We discuss
estimation via a generalized estimating equations framework in Section 3, and
inference via model misspecification-robust testing in Section 4. We present
simulation results pertaining to Type 1 and 2 error rates in Section 5 and
consider a data analysis in Section 6. We conclude with a discussion in
Section 7. Open-source software implementing the method is available at
github.com/statdivlab/radEmu, and code to reproduce the analyses are available
at github.com/statdivlab/radEmu_supplementary.

##  2 Model

###  2.1 Motivation

Our work is motivated by settings in which an “ideal” dataset $W_{ij}\geq 0$
(samples $i=1,\dots,n$ and categories $j=1,\dots,J$) is not observed, but a
distorted proxy $Y_{ij}\geq 0$ is available. Observing $Y$ rather than $W$
represents a loss of information in two critical ways. Firstly, the magnitude
of $Y_{i}\in\mathbb{R}^{J}_{\geq 0}$ is uninformative with respect to the
magnitude of $W_{i}\in\mathbb{R}^{J}_{\geq 0}$. In addition, the composition
of $Y_{i}$ may be systematically distorted relative to $W_{i}$ as a result of
differential detection. That is, observations from category $j_{1}$ may be
more readily observed than observations from category $j_{2}$, and the
magnitude of the detection effects are unknown. Our goal is to estimate fold-
differences in $\mathbb{E}W_{\cdot j}$ across covariates, that is,
$\mathbb{E}W_{\\{i:X_{i}=x\\}j}/\mathbb{E}W_{\\{i:X_{i}=x^{\prime}\\}j}$. We
formalize our model and approach to estimating fold-differences in Section
2.2.

We our motivating example of microbial communities, $W_{ij}$ represents the
(unobserved) cell or DNA concentration of microbial taxon $j$ in sample $i$,
which are prohibitively time- and labor-intensive to observe for all taxa
$j=1,\ldots,J$ (Williamson et al., 2022). However, it is relatively
straightforward to perform high-throughput sequencing to produce data
$Y_{ij}$, such as the (observed) read depth of taxon $j$ in sample $i$ in
shotgun sequencing, or the number of times sequence variant $j$ is observed in
sample $i$ based on amplicon sequencing. In the case-control microbiome
studies considered in Section 6, we are interested in estimating the fold-
difference in the mean cell concentration of taxa $j$ when comparing cases to
otherwise similar controls.

###  2.2 Model statement and parameter identifiability

Suppose that for pairs $(X_{i},W_{i})$ with $X_{i}\in\mathbb{R}^{p}$ and
$W_{i}\in\mathbb{R}^{J}_{\geq 0}$,

$\displaystyle\log\mathbb{E}[W_{i}|X_{i},\beta^{\star}]=X_{i}\beta^{\star}.$
(1)

Our target of inference is $\beta^{\star}\in\mathbbm{R}^{p\times J}$.
Unfortunately, however, we do not observe $W_{i}$, and so we cannot directly
estimate $\beta^{\star}$ via (1). Instead, we instead observe a perturbed
outcome $Y_{i}\in\mathbb{R}^{J}_{\geq 0}$ such that

$\displaystyle\log\mathbb{E}[Y_{i}|X_{i},\beta^{\star},z_{i},\delta]=z_{i}\mathbf{1}_{J}^{T}+X_{i}\beta^{\star}+\delta^{T},$
(2)

where sample-specific scalings $z_{i}\in\mathbb{R}$ and category-specific
detection effects $\delta\in\mathbbm{R}^{J}$ are unknown. Hence, the mean of
$Y_{ij}$ is equal to the mean of $W_{ij}$ subject to scaling by the sample-
specific factor $\exp(z_{i})$ and category-specific factor $\exp(\delta_{j})$.

As our target of inference is $\beta^{\star},$ we consider
$\mathbf{z}=[z_{1},\dots,z_{n}]^{T}$ and $\delta$ as nuisance parameters.
$\delta$ is straightforward to address if $\mathbf{1}_{n}$ lies in the column
space of $\mathbf{X}=[X_{1}^{T},\dots,X_{n}^{T}]^{T}$, as is typically the
case in practice. Without loss of generality, we define
$\mathbf{X}=[\mathbf{1}_{n}~{}\mathbf{D}]\in\mathbb{R}^{n\times p}$ and
$\mathbf{Y}=[Y_{1}^{T},\dots,Y_{n}^{T}]^{T}\in\mathbb{R}^{n\times J}$, and
rewrite model (14) as follows:

$\displaystyle\log\mathbb{E}[\mathbf{Y}|\mathbf{X}]$
$\displaystyle=\mathbf{z}\mathbf{1}_{J}^{T}+\mathbf{1}_{n}\beta_{0}+\mathbf{D}\beta^{\star(-0)}=\mathbf{z}\mathbf{1}_{J}^{T}+\mathbf{X}\beta$
(3)

for $\beta_{0}=\beta_{0}^{\star}+\delta^{T}$ and
$\beta=[\beta_{0}^{T},\beta_{1}^{\star T},\dots,\beta_{p-1}^{\star T}]^{T}$,
for $\beta_{k}^{\star}\in\mathbb{R}^{1\times J}$ the $k+1$-th row of
$\beta^{\star}$. In particular, because $\beta_{0}$ depends on the nuisance
detection effects $\delta$, it has no scientifically meaningful
interpretation, but $\beta_{1},\ldots,\beta_{p-1}$ contain useful information
regarding $\beta_{1}^{\star},\ldots,\beta_{p-1}^{\star}$, as discussed in more
detail below.

We now turn our attention to $\mathbf{z}$, which is not identifiable under our
model because for any $\alpha\in\mathbb{R}^{p}$, we have that

$\displaystyle z_{i}\mathbf{1}_{J}^{T}+X_{i}(\beta+\alpha\mathbf{1}_{J}^{T})$
$\displaystyle=z^{\dagger}_{i}\mathbf{1}_{J}^{T}+X_{i}\beta$ (4)

for $z^{\dagger}_{i}:=(z_{i}+X_{i}\alpha)$. Hence if
$\beta^{\dagger}=\beta+\alpha\mathbf{1}_{J}^{T}$ for some
$\alpha\in\mathbb{R}^{p}$, we can always find a $z_{i}^{\dagger}\in\mathbb{R}$
such that

$\displaystyle\log\mathbb{E}[Y_{i}|X_{i},\beta,z_{i}]=\log\mathbb{E}[Y_{i}|X_{i},\beta^{\dagger},z_{i}^{\dagger}],$
(5)

and therefore mean function (14) is not identifiable in $\beta$ and $z_{i}$.
However, mean function (14) is partially identifiable because we can define
equivalence classes of $\beta$’s such that for $\beta$ and $\beta^{\dagger}$
in the same equivalence class, there always exist $z_{i}$ and
$z_{i}^{\dagger}$ such that (5) holds, and for $\beta$ and $\beta^{\dagger}$
in distinct classes, there never exist such $z_{i}$ and $z_{i}^{\dagger}$. We
now make the form of these classes explicit and show that they have the
properties discussed above.

######  Lemma 1.

For arbitrary $b\in\mathbb{R}^{p\times J}$, define

$G_{b}=\\{\beta\in\mathbbm{R}^{p\times
J}:\beta=b+\alpha\mathbf{1}_{J}^{T}\text{ for some
}\alpha\in\mathbb{R}^{p}\\}.$

If $\mathbf{X}$ has full column rank, for any
$\beta,\beta^{\dagger}\in\mathbb{R}^{p\times J}$ there exist
$\mathbf{z}\in\mathbb{R}^{n}$ and $\mathbf{z}^{\dagger}\in\mathbb{R}^{n}$ such
that under model (14),
$\log\mathbb{E}[\mathbf{Y}|\mathbf{X},\beta,\mathbf{z}]=\log\mathbb{E}[\mathbf{Y}|\mathbf{X},\beta^{\dagger},\mathbf{z}^{\dagger}]$
if and only if $\beta^{\dagger}\in G_{\beta}$.

We prove Lemma 1 in Supporting Information Section 1. Lemma 1 demonstrates
that while we cannot uniquely estimate $\beta$ from data generated according
to (14), we can estimate it up to the addition of constant terms to each of
its rows. That is, while we cannot estimate differences across levels of
covariate $X_{\cdot k}$ in the log-means of $\mathbf{W}$ given only data
$\mathbf{Y}$, we can estimate relative differences in the log-mean of
$\mathbf{W}$ in category $j_{1}$ compared to in category $j_{2}$. In
microbiome sequencing studies, this represents a small loss of information
because researchers typically aim to identify microbial taxa that show
markedly different abundances across covariate groups of interest when
compared to the typical magnitude of these differences. Therefore, estimating
the equivalence class of $\beta$ answers important scientific questions about
microbial abundances.

We estimate the equivalence class of $\beta$ by imposing identifiability
constraints of the form $g(\beta_{k})=0$ on the rows $\beta_{k}$ of $\beta$
such that each $G_{\beta}$ contains only a single unique element satisfying
these constraints. Fitted means produced by our model are invariant to choice
of constraint, but this choice does change the interpretation of elements of
$\beta$. For example, under constraint $\beta_{kj}=0$, $\beta_{kj^{\prime}}$
represents $\log\mathbb{E}[W_{\cdot
j^{\prime}}|X_{\cdot}=\mathbf{x}_{x_{k}+1}]-\log\mathbb{E}[W_{\cdot
j^{\prime}}|X_{\cdot}=\mathbf{x}]$ minus the corresponding quantity for taxon
$j$, where $\mathbf{x}_{x_{k}+1}$ is equal to $\mathbf{x}$ in all elements
excepting the $k$th, which is increased by 1 unit. We equivalently interpret
$e^{\beta_{kj^{\prime}}}$ as a fold-difference (multiplicative difference).

Figure 1: Because the mean function is only partially identifiable, true
effect sizes $\beta$ cannot be estimated on an “absolute” scale (left). Full
identifiability is established via a constraint function, allowing us to
estimate the log-fold differences in true abundances across groups relative to
typical differences. The above illustrates the favorable qualities of
centering via the pseudo-Huber smoothed median (right) compared to centering
by the mean (middle). $\beta_{j}$ for $j=2,3,4,5$ are unlikely to be of
scientific interest (relative to $j=1,6$), and the smoothed median constrains
these “uninteresting” taxa closer to zero than mean centering.

If we aim to compare log-mean differences across levels of a covariate to
typical such differences across taxa, we might choose to impose
identifiability by requiring $\text{median}_{j=1,\ldots,J}(\beta_{kj})=0$. The
robustness of the median is attractive: if few taxa exhibit large differences
in mean concentration across levels of a covariate, this will in general
minimally impact a median effect if most effects are small. By contrast,
centering by a non-robust measure of location (e.g., mean) may result in
unwanted behavior, because a single large effect $\beta_{kj}$ will pull the
mean away from the bulk of taxa (Figure 1). However, the median is not a
differentiable function of $\beta_{k}$ (its total derivative is not guaranteed
to exist), which substantially complicates inference (see Section 4).
Accordingly, we propose a smooth approximation to the median that shares its
robustness properties, that is, the pseudo-Huber loss
$g_{p}(\beta_{k})=\arg\min_{c}\sum_{j=1}^{J}h_{\delta}\left(\beta_{kj}-c\right)$
for
$h_{\delta}(x)=\delta^{2}\left(\sqrt{1+\left(x/\delta\right)^{2}}-1\right)$
for prespecified parameter $\delta>0$ (Barron, 2019). The pseudo-Huber loss is
approximately quadratic in a neighborhood of 0 and asymptotically approaches
the function $|x|/\delta$ for large $x$. The parameter $\delta$ controls the
size of the quadratic region around zero: as $\delta\downarrow 0$,
$g_{p}(\beta_{k})$ approaches the median, and as $\delta\rightarrow\infty$, it
approaches the mean. We propose $\delta=0.1$ as a reasonable default for many
applications, though empiricially we find that care is warranted when
covariates are on very different scales. While we propose $g_{p}$ as the
default constrained function, any smooth $g$ satisfying
$g(\beta_{k}+\alpha)=g(\beta_{k})+\alpha$ is reasonable, and the most
appropriate constraint will depend on the scientific setting.

##  3 Estimation

###  3.1 Estimation under $H_{A}$

We estimate $\beta$ (our target of inference) and $\mathbf{z}$ (a nuisance
parameter) using an estimating equations approach motivated by a Poisson
likelihood with mean function (14). To address issues of infinite or undefined
MLEs resulting from data separation, which commonly occurs in high-dimensional
microbiome sequencing data (Bull et al., 2002), we introduce a Firth penalty
on $\beta$ derived from a formal equivalence between our model and the
multinomial logistic model. We develop a fast algorithm to find the maximizer
of the penalized likelihood by efficiently solving a sequence of score
equations that converge to the score of the penalized likelihood within a
neighbourhood of the global maximizer. As we will see in Sections 4 and 5, our
estimator and accompanying hypothesis testing show excellent performance even
when data has excess-Poisson variance, and even in small samples. All
algorithms are provided as pseudocode in Supporting Information, and an R
package implementing the method is available at [GitHub Link].

To motivate our estimator, we first introduce a Poisson log likelihood for our
model. Letting $\beta^{j}$ denote the $j$-th column of $\beta$ and $X_{i}$
denote the $i$-th row of $\mathbf{X}$, the Poisson log likelihood is

$\displaystyle l^{\text{Poisson}}(\beta,z)$
$\displaystyle=\sum_{j=1}^{J}\Big{[}\sum_{i=1}^{n}\Big{(}Y_{ij}(X_{i}\beta^{j}+z_{i})-\exp(X_{i}\beta^{j}+z_{i})\Big{)}\Big{]}.$
(6)

For fixed $\beta$, the maximizing value of $\mathbf{z}$ is available in closed
form

$\displaystyle\hat{z}_{i}$
$\displaystyle=\log\sum_{j=1}^{J}Y_{ij}-\log\sum_{j=1}^{J}\exp(X_{i}\beta^{j}),$
(7)

and therefore we can profile out $\mathbf{z}$ in likelihood (6) to give

$\displaystyle l^{\text{profile}}(\beta)$
$\displaystyle=\sum_{j=1}^{J}\Big{[}\sum_{i=1}^{n}\Big{(}Y_{ij}\log\frac{\exp(X_{i}\beta^{j})}{\sum_{j^{\prime}=1}^{J}\exp(X_{i}\beta^{j^{\prime}})}-Y_{i\cdot}\frac{\exp(X_{i}\beta^{j})}{\sum_{j^{\prime}=1}^{J}\exp(X_{i}\beta^{j^{\prime}})}\Big{)}\Big{]}+C$
(8)

where $C$ is constant (with respect to $\beta$) and
$Y_{i\cdot}:=\sum_{j=1}^{J}Y_{ij}$. Note that the profile likelihood (8) is
equal to a multinomial log likelihood with a logistic link (up to a constant).
This accords with both known results about the marginal Poisson distribution
of multinomial random variables (Birch, 1963), and our observations on
identifiability of $\beta$ (only $p\times(J-1)$ parameters can be identified
in a multinomial logistic regression of a $J$-dimensional outcome on $p$
regressors).

Unfortunately, the maximizer of $l^{\text{profile}}$ will not be finite if
there exists a perfect separating hyperplane (Albert and Anderson, 1984). To
both guarantee the existence of an estimate and reduce its bias, we derive the
Firth penalty associated with likelihood (8), observing the impact of the
nuisance parameters $\mathbf{z}$ on the correct bias-reducing penalty
(Kosmidis and Firth, 2011) . Fix $j^{\dagger}\in\\{1,\ldots,J\\}$ and let
$\tilde{\beta}:=[{\beta^{1}}^{T},\ldots,{\beta^{j^{\dagger}-1}}^{T},{\beta^{j^{\dagger}+1}}^{T},\ldots,{\beta^{J}}^{T}]^{T}$,
that is, a column vector containing stacked columns of $\beta$ but with
$\beta^{j^{\dagger}}$ omitted. Let
$\tilde{\mathbf{X}}_{i}=I_{J}^{[-j^{\dagger}]}\otimes X_{i}$, where
$I_{J}^{[-j^{\dagger}]}$ indicates a $J\times J$ identity matrix whose
$j^{\dagger}$-th column has been removed and $\otimes$ is the Kronecker
product, so $\tilde{\mathbf{X}}_{i}\in\mathbbm{R}^{J\times p(J-1)}$ and
$\tilde{\mathbf{X}}:=[\tilde{\mathbf{X}}_{1}^{T},\dots,\tilde{\mathbf{X}}_{n}^{T}]^{T}\in\mathbb{R}^{nJ\times
p(J-1)}$. Define the fitted values
$p_{i}=\exp(\tilde{\mathbf{X}}_{i}\tilde{\beta})\big{/}\mathbf{1}_{J}^{T}\exp(\tilde{\mathbf{X}}_{i}\tilde{\beta})\in\mathbb{R}^{J}$.
We can then show that the information matrix associated with (8) is
$\tilde{\mathbf{X}}^{T}\mathbf{V}\tilde{\mathbf{X}}$ where $\mathbf{V}$ is a
block diagonal matrix with $i$-th $J\times J$ block equal to
$Y_{i\cdot}(\text{diag}(\mathbf{p}_{i})-\mathbf{p}_{i}\mathbf{p}_{i}^{T})$
with $\text{diag}(\mathbf{p}_{i})$ indicating a diagonal matrix with diagonal
elements equal to $p_{i1},\dots,p_{iJ}$. Finally, we can express the Firth
penalized multinomial log likelihood as follows:

$\displaystyle
l^{\text{penalized}}=l^{\text{profile}}+\frac{1}{2}\log\big{|}\tilde{\mathbf{X}}^{T}\mathbf{V}\tilde{\mathbf{X}}\big{|},$
(9)

where the penalty term depends on $\beta$ through $\mathbf{V}$.

Remarkably, the penalized likelihood (9) is invariant to the choice of
$j^{\dagger}$. This is because the Firth penalty is invariant under smooth
reparametrizations for full exponential family models such as the multinomial
distribution (Jeffreys, 1946; Firth, 1993), and the above parameterization
implies the constraint $g(\beta_{k})=\beta_{kj^{\dagger}}=0$. Thus, we can
estimate $\beta$ under the “convenience constraint” $\beta_{kj^{\dagger}}=0$
and impose the desired identifiability constraint post-hoc by substracting the
relevant $g(\tilde{\beta}_{k})$ from each element of $\tilde{\beta}$. In
practice we find that choosing $j^{\dagger}$ to maximize
$\sum_{i}\mathbf{1}_{[Y_{ij^{\dagger}}>0]}$ improves the speed of
optimization.

To find $\arg\max_{\tilde{\beta}}l^{\text{penalized}}$, we apply the iterative
maximum likelihood algorithm of Kosmidis and Firth (2011), which maximizes (9)
by iteratively solving multinomial score equations in augmented data
$\mathbf{Y}^{(t)}$. Specifically, $Y^{(t)}_{ij}=Y_{ij}+h_{ij}/2$ where
$h_{ij}$ is the $(i-1)\times J+j$-th diagonal element of the Poisson “hat”
matrix
$\mathbf{V}_{p}^{1/2}\tilde{\mathbf{G}}(\tilde{\mathbf{G}}^{T}\mathbf{V}_{p}\tilde{\mathbf{G}})^{-1}\tilde{\mathbf{G}}^{T}\mathbf{V}_{p}^{1/2}$
evaluated at current values of $\beta^{(t)}$ for $\mathbf{V}_{p}$ an $nJ\times
nJ$ diagonal matrix with $(i-1)\times J+j$-th diagonal element equal to
$\exp(X_{i}\beta^{j}+z_{i})$ and
$\tilde{\mathbf{G}}=\Big{[}\tilde{\mathbf{X}}:\mathbf{I}_{n}\otimes\mathbf{1}_{J}\Big{]}\in\mathbb{R}^{nJ\times(p(J-1)+n)}$
(i.e., the expanded design matrix for $\beta$, $\tilde{\mathbf{X}}$, side by
side with a design matrix for $\mathbf{z}$,
$\mathbf{I}_{n}\otimes\mathbf{1}_{J}$). We compute data augmentations under
the restriction on $z$ that $\sum_{j}\exp(X_{i}\beta^{j}+z_{i})=Y_{i\cdot}$ by
applying (7). Each iteration $t$ of this algorithm is separable in
$z,\beta^{1},\ldots,\beta^{J}$, and we therefore solve each iteration using
coordinate descent, applying (7) to update $z$ and updating each $\beta^{j}$
with Fisher scoring followed by a line search. When a desired convergence
threshold is reached, we impose the desired identifiability constraint by
centering by $g(\tilde{\beta}_{k}^{(t_{\text{max}})})$ as described
previously.

###  3.2 Estimation under $H_{0}$

In Section 4, we present a robust score test for
$\beta_{k^{\star}j^{\star}}=0$ under constraint $g(\beta_{k^{\star}})=0$,
which requires fitting a model under the restriction
$\beta_{k^{\star}j^{\star}}=g(\beta_{k^{\star}})$. We therefore develop an
augmented Lagrangian algorithm (Nocedal and Wright, 2006, p. 514) to solve
this constrained optimization problem. The augmented Lagrangian function is

$\displaystyle\mathcal{L}(\beta,z;\rho,u)$
$\displaystyle:=-l^{\text{Poisson}}(\beta,z;\mathbf{Y}^{(t_{\text{max}})})+u\left[g\left(\beta_{k^{\star}}\right)-\beta_{k^{\star}j^{\star}}\right]+\frac{\rho}{2}\left[g\left(\beta_{k^{\star}}\right)-\beta_{k^{\star}j^{\star}}\right]^{2}$
(10)

for $\mathbf{Y}^{(t_{\text{max}})}$ from our model fit under $H_{A}$. For
sequences $\rho^{(t)}$ and $u^{(t)}$ described below, we iteratively find
maximizers $(\hat{\beta}_{aug}^{(t)},\hat{z}_{aug}^{(t)})$ of (10) until
$\hat{\beta}_{aug}^{(t)}$ satisfies the constraint
$\beta_{k^{\star}j^{\star}}=g(\beta_{k^{\star}})$ to within a prespecified
tolerance. While ideally (10) would be a function of
$l^{\text{penalized}}(\beta;\mathbf{Y})$ instead of
$l^{\text{Poisson}}(\beta,z;\mathbf{Y}^{(t_{\text{max}})})$, updating data
augmentations within Lagrangian steps adds an additional inner loop, and
furthermore $l^{\text{Poisson}}(\beta,z;\mathbf{Y}^{(t_{\text{max}})})$ is
both a good approximation to $l^{\text{penalized}}(\beta;\mathbf{Y})$ and
asymptotically equivalent under $H_{0}$.

The presence of $g(\beta_{k^{\star}})$ in (10) complicates optimization, as
this term depends on an entire row of $\beta$. As a result, updating $\beta$
column by column (as when fitting models under $H_{A}$) often results in slow
convergence under $H_{0}$. Hence, we instead minimize (10) via coordinate
descent alternating between block updates to $\beta$ via an approximate Newton
step and updates to $z$ via equation (7).

For fixed $z$, the augmented Lagrangian (10) has second derivative in
$\tilde{\beta}$ (the vectorization of $\beta$ with convenience constraint
$\beta^{j^{\dagger}}=\mathbf{0}$) given by

$\displaystyle\nabla^{2}\mathcal{L}$
$\displaystyle=-\nabla^{2}l^{\text{Poisson}}+[u+\rho(g-\beta_{k^{\star}j^{\star}})]\nabla^{2}g+\rho(\nabla
g-\vec{e}_{[k^{\star}j^{\star}]})(\nabla
g-\vec{e}_{[k^{\star}j^{\star}]})^{T}$ (11)

where $\vec{e}_{[k^{\star}j^{\star}]}$ is a vector equal to 1 in the position
in $\tilde{\beta}$ corresponding to $\beta_{k^{\star}j^{\star}}$ and zero
elsewhere. $\nabla^{2}l^{\text{Poisson}}$ is block diagonal and hence
inexpensive to invert, and $\rho(\nabla
g-\vec{e}_{[k^{\star}j^{\star}]})(\nabla
g-\vec{e}_{[k^{\star}j^{\star}]})^{T}$ is a matrix of rank one. We therefore
approximate the Hessian of (10) as $-\nabla^{2}l+\rho(\nabla
g-\vec{e}_{[k^{\star}j^{\star}]})(\nabla^{T}g-\vec{e}_{[k^{\star}j^{\star}]})$,
and compute its inverse as a rank-one update to $-(\nabla^{2}l)^{-1}$ via the
Sherman-Morrison identity (Sherman and Morrison, 1950). For additional
details, see Supporting Information Section 3.3.

In each iteration of the outer loop of our algorithm, we update $\rho$ and $u$
and obtain an approximate solution to (10). The update to $u$ is given by
$u^{(t+1)}=\rho^{(t)}[g(\beta_{k^{\star}}^{(t)})-\beta^{(t)}_{k^{\star}j^{\star}}]$
We update $\rho$ if we have not progressed sufficiently toward a feasible
solution in the past step: if
$|g(\beta_{k^{\star}}^{(t)})-\beta^{(t)}_{k^{\star}j^{\star}}|>\kappa|g(\beta_{k^{\star}}^{(t-1)})-\beta^{(t-1)}_{k^{\star}j^{\star}}|$,
we set $\rho^{(t+1)}=\tau\rho^{(t)}$ for $\kappa\in(0,1)$ and $\tau>1$;
otherwise, we let $\rho^{(t+1)}=\rho^{(t)}$. We find that default values
$\kappa=0.8$ and $\tau=2$ yield stable optimization in a broad range of
settings. Lower $\kappa$ and higher $\tau$ yield potentially faster
optimization but at the cost of stability. In particular, if $\rho$ grows too
large, the augmented Lagrangian subproblem may become ill-conditioned.

We evaluate convergence with two criteria. Firstly, we require (primal)
feasibility be satisfied to within a chosen tolerance:
$|g(\beta_{k^{\star}})-\beta_{k^{\star}j^{\star}}|<\varepsilon$ for a small
$\varepsilon>0$. We additionally require that
$\max_{k,j}|\beta^{(t)}_{kj}-\beta^{(t-1)}_{kj}|<\Delta$ for a small
$\Delta>0$.

##  4 Inference

We now discuss evaluating evidence for a difference in $\mathbb{E}W_{j}$
across levels of a covariate $X_{k}$, holding other covariates fixed. We
explore the simple but common setting of testing
$\overset{\circ}{H}_{0}:\beta_{kj}=0$, but the approach can be generalized to
test more general hypotheses of the form $\mathbf{A}\beta=\mathbf{C}$ for
fixed matrices $\mathbf{A}\in\mathbb{R}^{h\times p}$ and
$\mathbf{C}\in\mathbb{R}^{h\times J}$. We are not able to test
$\overset{\circ}{H}_{0}$ directly, as $\beta$ is not fully identified.
However, we are able to test $H_{0}^{(g)}:\beta_{kj}=g(\beta_{k})$, or
equivalently

$\displaystyle H_{0}:\text{ }\log\mathbb{E}[$ $\displaystyle
W_{j}|X_{k}=x_{k},X_{-k}=x_{-k}]-\log\mathbb{E}[W_{j}|X_{k}=x_{k}+1,X_{-k}=x_{-k}]$
$\displaystyle=g(\log\mathbb{E}[W|X_{k}=x_{k},X_{-k}=x_{-k}]-\log\mathbb{E}[W|X_{k}=x_{k}+1,X_{-k}=x_{-k}])$

Thus, the interpretation of the test depends on the identifiability
constraint. For example, we may test that a log-fold difference equals the
smoothed median of such differences via $g(\beta_{k})=g_{p}(\beta_{k})$, or
that a log-fold difference equals the log-fold difference of category 1 via
$g(\beta_{k})=\beta_{k1}$.

To construct a test of $H_{0}^{(g)}$, we restate $H_{0}^{(g)}$ as follows:
$H_{0}^{(g)}:h(\beta)=0$ with $h(\beta)=\beta_{kj}-g(\beta_{k})$. We test
nulls of this form via a robust score test, using the test statistic given in
White (1982):

$\displaystyle
T_{RS}=S_{H_{0}}^{T}I_{H_{0}}^{-1}F_{H_{0}}^{T}(F_{H_{0}}I_{H_{0}}^{-1}D_{H_{0}}I_{H_{0}}^{-1}F_{H_{0}}^{T})^{-1}F_{H_{0}}I_{H_{0}}^{-1}S_{F_{H_{0}}}$
(12)

where $S_{H_{0}}$ is the score evaluated at the maximum likelihood estimate
$\hat{\beta}_{0}$ under the null, $I_{H_{0}}$ is a consistent estimate of the
information matrix under the null, $F_{H_{0}}$ is given by $\frac{\partial
h}{\partial\beta^{T}}$ evaluated at $\hat{\beta}_{0}$ and
$D_{H_{0}}=\sum_{i=1}^{n}s_{i}(\hat{\beta}_{0})s_{i}(\hat{\beta}_{0})^{T}$
where $s_{i}(\beta)=\tilde{X}_{i}^{T}(Y_{i}-Y_{i\cdot}p_{i}(\beta))$ is the
contribution to the score of $Y_{i}$. In practice, we use a score statistic
$T^{\prime}_{RS}$ equal to $\frac{n}{n-1}T_{RS}$, where the term in front of
$T_{RS}$ is an adjustment due to Guo et al. (2005) that reduces the small-
sample conservatism of the test. Under $H_{0}$, $T^{\prime}_{RS}$ is
asymptotically $\chi^{2}$-distributed with one degree of freedom; empirically,
we observe that this test provides error rate control in large samples and
conservative inference in small to moderate samples (Section 5). We also
consider a robust Wald approach to testing hypotheses of the form
$H_{0}^{(g)}:h(\beta)=0$. The robust Wald statistic is given by

$\displaystyle
T_{RW}=h(\beta)^{T}(F_{H_{A}}^{T}I_{H_{A}}D_{H_{A}}^{-1}I_{H_{A}}F_{H_{A}})^{-1}h(\beta)$
(13)

where $F_{H_{A}}$, $I_{H_{A}}$, and $D_{H_{A}}$ are as above but evaluated
under the full model. Since robust Wald testing does not require refitting a
model under a null hypothesis, it offers a computational advantage over robust
score testing. However, as explored in Section 5, robust Wald tests fail to
adequately control Type 1 error in small samples.

##  5 Simulations

We conduct a small simulation study to examine the Type 1 error and power of
our proposed tests. We simulate in two distributional settings: a correctly
specified setting in which counts are independently Poisson-distributed, and a
misspecified setting in which the mean model is correctly specified but counts
follow independent zero-inflated negative binomial (ZINB) distributions. To
obtain draws from a ZINB distribution, we first draw from a Negative Binomial
distribution with mean $\mu$ given by our model and dispersion parameter
$\phi=5$ (thus, their variance is $\mu+\mu^{2}/\phi$). We multiply this draw
with an independent Bernoulli draw $\xi$ with $Pr(\xi=1)=0.4$ to obtain our
final ZINB draw.

Figure 2: Q-Q plots comparing empirical quantiles (y-axis) of the robust score
(dark blue) and robust Wald (light blue) p-values to theoretical quantiles
(x-axis). A conservative test will produce a curve above the line $x=y$ for
small p-values and an anti-conservative test will produce a curve below it. To
aid in comparison of lower quantiles, both axes are presented on a square-root
scale. 10,000 simulations were performed.

In all simulations, we simulate observations in two groups of equal size and
test nulls of the form $\beta_{2j}-g_{p}(\beta_{2})=0$, where $\beta_{2}$
corresponds to between-group differences in log mean (up to a constant shift).
The pseudo-Huber smoothing parameter for $g_{p}$ is 0.1 in all simulations.
The first row of $\beta$ is fixed so that its odd elements form an linearly
increasing sequence from -3 to 3 and its even elements form a similarly
decreasing sequence from 3 to -3. We construct the first row of $\beta$ in
this way so that it is not strongly correlated with the second row of $\beta$.
In null simulations, the second row of $\beta$, save its $J/2$\- and
$(J+1)/2)$-th elements, which we set equal to zero, is given by
$5\times\text{sinh}(x)/\text{sinh}(10)$ where $x$ consists of linearly
increasing sequence from $-10$ to $10$, and sinh indicates the hyperbolic sine
function. Under alternatives, the second row of $\beta$ is the same as just
described with the exception of its $J/2$-th element, which is set equal to
$1$ (under the “weak alternative”) or $5$ (under the “strong alternative”).
The values of $\beta$ common to all three hypotheses were chosen to reflect a
setting in which most taxa vary similarly in mean concentration across two
groups, with a few taxa deviating from this typical variation. The values of
$\beta_{k^{\star}j^{\star}}$ (i.e., the element of $\beta$ for which we we
test equality to $g(\beta_{k^{\star}})$) chosen for “weak” and “strong”
alternatives are intended to reflect a modest ($\approx$ 3-fold deviation from
typical ratio of mean concentrations across groups) effect that in some
settings might nonetheless be of scientific interest, as well as a strong
($\approx$ 150-fold deviation from typical) effect. We simulate with total
number of taxa $J\in\\{10,50,25\\}$ and total sample size
$n\in\\{10,50,250\\}$. Note that total sample size indicates count of
independent observations across both groups, so $n=10$ corresponds to $5$
observations in each of the two groups.

Figures 2 and 3 summarize the results of our simulations under the null and
alternatives, respectively, with tabular results in Supporting Information
Section 4. As expected, the robust score test exhibits conservatism in
simulations under the null, with greater conservatism at smaller sample sizes
and with ZINB-distributed counts. By contrast, the robust Wald test is anti-
conservative in the smallest samples, with empirical Type 1 error 0.307 in
ZINB simulations and 0.167 in Poisson simulations with $n=10$, $J=250$. In
these two simulation settings, robust score tests have empirical Type 1 error
0.025 and 0.048, respectively. Notably, while the robust Wald test tends
toward greater anti-conservatism when counts are not Poisson distributed, the
robust score test appears to tend toward greater conservatism. At the highest
sample size we considered, $n=250$, empirical Type 1 error is close to nominal
regardless of distribution, number of taxa, or test.

Figure 3: Empirical power to reject $H_{0}$ at the 5% level, estimated from
500 simulations. The panels at left show power to detect a smaller effect
($\approx$ 3-fold ratio of mean concentrations), and the panels at right show
power to detect a larger effect ($\approx$ 150-fold ratio of mean
concentrations). See Supporting Information Figure 1 for power curves.

In simulations under the alternative, we observe higher power when counts are
drawn from a Poisson than from a ZINB distribution, and when sample size is
larger. While the robust Wald test exhibits higher power than the robust score
test at smaller sample sizes, it also fails to control Type 1 error, so we
nonetheless do not recommend this test for use in small or moderate samples.
We also observe a modest impact of total number of taxa $J$ on power: for ZINB
simulations, we observe nondecreasing power in $J$. However, power of the
score test does appear to decrease in $J$ for Poisson simulations with $n=10$.

##  6 Data analysis: Colorectal cancer microbial metaanalysis

In this section we consider a shotgun metagenomic dataset published by Wirbel
et al. (2019). This dataset comprises a meta-analysis of five studies, each
targeting possible associations between the composition of the human gut
microbiome and colorectal cancer (CRC) risk. In each study, fecal samples were
taken from participants with colorectal cancer as well as from healthy
controls. We are interested in identifying microbial strains unusually
enriched or depleted in study participants with colorectal cancer as compared
to otherwise similar control participants.

Each observation $Y_{ij}\in\mathbb{R}_{\geq 0}$ reflects the observed
abundance of microbial strain $j$ in fecal material of study participant $i$,
measured as the depth of coverage of marker genes calculated by the mOTU
profiler v2 (Milanese et al., 2019). We consider the “ideal” but unobserved
data $W_{ij}$ to be the cell concentration of strains in the stool of study
participants. Participant data was pooled from five studies, yielding 575
study participants among whom 2342 unique strains were detected. In addition
to metagenomic data, Wirbel et al. (2019) published clinical and demographic
metadata including colorectal cancer status, age in years, binary gender, BMI,
and timing of sampling relative to colonoscopy (sample taken prior to or
following colonoscopy). Wirbel et al. (2019) applied stringent filtering
criteria to remove low-abundance strains from their dataset. While not
necessary for our methodology, to enable the comparison of our results with
theirs, we applied the same filtering criterion, excluding all mOTUs that
failed to exceed empirical relative abundances of $10^{-3}$ in at least three
studies. We excluded 9 participants with missing covariate data (8 for BMI and
1 for timing of sampling relative to colonoscopy). With these exclusions, we
analyze $n=566$ study participants and $J=845$ metagenomic strains.

We fit our proposed model to the untransformed mOTU data of Wirbel et al.
(2019), applying Firth penalization as well as a smoothed median
identifiability constraint to the rows of $\beta$. We include, in addition to
an intercept and a term for CRC status (1 for CRC and 0 otherwise), terms in
our model for age (linear spline with a single knot at median age 64 years),
BMI (linear spline with a single knot at median BMI 24.7 kg/m2), gender (1 for
Male and 0 otherwise), study (4 indicator variables denoting the US, CN, DE
and AT studies; we take FR to be the reference study), and whether fecal
samples were taken after colonoscopy (1 if after, 0 if prior). Thus,
coefficients on CRC status estimate the difference in log mean cell
concentration comparing CRC patients to non-CRC controls who are alike in
gender, BMI, age, study, and timing of stool sampling relative to colonoscopy,
relative to typical such differences across taxa. We perform robust score
tests and control the false discovery rate (using the Benjamini and Hochberg
(1995) (BH) procedure, following Wirbel et al. (2019)) at $0.1$.

Figure 4: 95% Wald confidence intervals for taxon-specific CRC status effects
estimated in our multivariable model (y-axis), plotted against the proportion
of participants in all studies in which each taxon was detected (x-axis). Taxa
with Benjamini-Hochberg-adjusted robust score p-values below 0.1 are shown in
dark blue, with all others shown in light blue. Significantly differentially
abundant taxa are identified across levels of observed prevalence.

Robust Wald confidence intervals for the coefficients on CRC status are shown
in Figure 4. We identify 30 taxa for which we reject the null hypothesis of
equal log mean cell concentration for CRC patients and (otherwise alike) non-
CRC controls at a false discovery rate of 0.1 based on the robust score test.
These taxa are denoted by navy-colored confidence intervals. Taxa are ordered
along the x-axis by the proportion of participants in which they are detected,
highlighting that significantly differentially abundant taxa are identified
across the spectrum of detection levels. Taxa our analysis identifies as more
enriched than typical among CRC patients relative to similar controls include
an unknown Dialister species estimated to have ratio of mean cell
concentration 49 times greater than typical (95% CI $12-210$; $p=1.1\times
10^{-6}$), Porphyromonas uenonis (30 times greater; 95% CI $6.2-150$;
$p=0.0025$), Fusobacterium nucleatum s. animalis (24 times greater; 95% CI
$11-52$; $p=2.2\times 10^{-5}$), and Gemella morbillorum (23 times greater;
95% CI $7.6-69$; $p=0.0016$). Confidence intervals for all effect sizes that
are significantly associated with CRC status at FDR cutoff $0.1$ are given in
Figure S7.

Figure 5: Taxa associated with CRC status in adjusted model at FDR level 0.1.
Species designations are given on the x-axis; multiple confidence intervals
are reported for unknown Clostridiales as multiple mOTUs with no known species
designation were identified as significantly associated with CRC status in
this class. 95% Wald confidence intervals for effect estimates (taxon-specific
difference from smoothed median difference in log mean cell concentrations
comparing CRC patients to controls, all other covariates held constant) are
shown on the y-axis, with color indicating proportion of study participants in
whom a given taxon was detected.

We identify a number of taxa associated with CRC status that are detected in
relatively high proportion in the study populations. For example, we estimate
for B. kashiwanohense, which was detected in 35.7% of study participants, that
the ratio of mean stool cell concentrations for this taxon comparing CRC
patients to alike controls is 8.6 times lower (95% CI 3.6 – 20.9; unadjusted
robust score $p=0.0022$) than the typical such ratio among taxa included in
our analysis. This is consistent with a possible protective effect of B.
kashiwanohense, but we cannot rule out association due to unmeasured
confounding (e.g., through diet) or reverse causation (development of CRC
induces relative depletion of B. kashiwanohense). We note that B.
kashiwanohense was not identified as a taxon negatively associated with CRC
status in the analysis of Wirbel et al. (2019).

We compare our results to those of Wirbel et al. (2019), as well as to 5 other
popular differential abundance methods in the microbiome literature, in Figure
6. All panels show robust score p-values obtained from our adjusted model on
the x-axis, and p-values obtained from comparator methods on the y-axis. We
briefly describe details for each method here; please see Code Appendix for a
reproducible workflow. (Top left) corncob (Martin et al., 2020) likelihood
ratio test p-values with design matrix as for the proposed method and timing
of sampling and study as dispersion covariates; (top center) ALDEx2 (Fernandes
et al., 2014) with centered log-ratio transform and design matrix as for the
proposed method and 512 Monte Carlo iterations; (top right) DESeq2 (Love et
al., 2014) likelihood ratio test with design matrix as for the proposed method
and “local” estimation of dispersions; (lower left) IFAA (Li et al., 2021)
with design matrix as for the proposed method; (lower center) ANCOM-BC2 (Lin
and Peddada, 2020) with design matrix as for the proposed method and without
prevalence filtering; and (lower right) the results of Wirbel et al. (2019),
who report a blocked Wilcoxon test performed on proportion-scale mOTU data,
with blocking on study and timing of sampling relative to colonoscopy (five
p-values of zero are altered to $10^{-25}$ for visual comparison). While all
methods target slightly different estimands, our proposed methodology is most
similar to the ALDEx2 and ANCOM-BC2 target estimands, which is consistent with
the relatively high sample Spearman correlation of p-values between our method
when compared with ALDEx2 (0.53) and ANCOM-BC2 (0.65).

Finally, we contrast our results with those of Wirbel et al. (2019), finding
substantial variation in the ordering of p-values between our proposed method
and the blocked Wilcoxon test (Spearman correlation $0.44$). Only 45 of the 94
taxa Wirbel et al. (2019) identified as significantly associated with CRC
status are among the 94 taxa most significantly associated with CRC status in
our analysis, and of the 30 taxa we highlight at FDR level $0.1$, only 10 are
among the 30 most highly significant taxa identified by Wirbel et al. (2019).
We believe that this lack of concordance is likely due to two major factors.
Firstly, our model adjusts for a larger set of potential confounders than do
Wirbel et al. (2019), which may account for our discrepant findings. In
addition, we estimate effects on the scale of log mean cell concentration,
whereas Wirbel et al. (2019) conduct tests on the read proportion scale.
Because a perturbation to the abundance of a small number of taxa on the cell
concentration scale induces a dense perturbation on the proportion scale,
analyses on the proportion scale may highlight taxa that do not vary in a
scientifically interesting way on the cell concentration scale. Similarly,
analyses performed on proportion data may fail to detect scientifically
meaningful effects on the cell concentration scale because of fluctuations in
cell concentration among other taxa.

Figure 6: Negative $\log_{10}$ p-values from distinct approaches to
differential abundance analyses based on the Wirbel et al. (2019) dataset. The
proposed method’s robust score test p-values are shown on the x-axis, and
comparator methods on the y-axis. The proposed method produces smaller
p-values in 45%, 68%, 52%, 49%, 51% and 41% of cases for corncob, ALDEx2,
DESeq2, IFAA, ANCOM-BC2 and blocked Wilcoxon respectively, with rank
correlations of 0.45, 0.53, 0.57, 0.14, 0.65 and 0.44. Wirbel et al. (2019)
report some p-values equal to zero; to allow graphing on the log scale, we
replace these values with $10^{-25}$. In addition to the missing/uncomputable
p-values of corncob, DESeq2 and IFAA, two DESeq2 were so small as to make
plotting infeasible; these observations correspond to
$(p_{\text{radEmu}},p_{\text{DESeq2}})$ tuples of $(10^{-1.8},10^{-31})$ and
$(10^{-2.8},10^{-33})$.

##  7 Discussion

In this paper we present a regression method for inference on means of a
nonnegative multivariate outcome $W$ observed after perturbation by unknown
sample-specific scaling terms $z$ and category-specific detection effects
$\delta$. Our method is motivated by microbiome sequencing experiments, in
which outcomes of interest (microbial cell concentration across many species)
are typically observed indirectly via sequencing, which is known to be subject
to both sample-specific scaling and taxon-specific distortion (McLaren et al.,
2019). However, we anticipate that our method may find application in other
settings; in particular, we expect that it may be appropriate for analysis of
various types of ’omics measurements, where similar scaling and detection
effects are likely to distort measurements.

Our method has a number of advantages over existing methods commonly applied
to perform differential abundance analyses in microbiome studies. Firstly, it
targets a population estimand that is distinguishable from sample- and
category-scaling effects and that exists under realistic conditions
(positivity of population mean abundances $\mathbb{E}W_{j}$). Notably, we do
not require positivity of true sample abundances $W_{ij}$. Secondly, our
method does not require data transformations, such as log-ratio
transformations or relative abundance transformations. This allows us to use
information regarding the precision of observations in estimation. It also
avoids the ad-hoc fixes (e.g., adding a small pseudocount to zero
observations, restricting attention to nonzero counts, or attempting to infer
whether observed zeroes are due to insufficient read depth or reflect true
absence of a taxon) that are necessary to apply log-ratio transformations to
data containing zeroes. Additionally, we employ a model-robust score test that
controls Type 1 error even in small samples with large excess-Poisson
dispersion and high degrees of sparsity. Finally, we incorporate an
identifiability constraint that may be flexibly specified to suit varying
scientific settings. In particular, we do not require a choice of reference
taxon. That said, if such a taxon is available, our method is able to take
advantage of this information.

An additional consideration in practical applications for this method is
measurement error. While our target estimand $\beta$ can be estimated in the
presence of systematic measurement error $\delta$, no such guarantees exist
for other kinds of measurement error, such as cross-sample contamination.
Accordingly, we encourage the use of decontamination techniques (Davis et al.,
2018; Clausen and Willis, 2022) prior to use of our proposed method.

Our method admits several possible extensions. In contexts where it is
scientifically plausible that the parameter of interest $\beta_{k}$ is sparse,
enforcing sparsity via $\ell_{1}$ penalization is an attractive option.
However, we note that if most elements of $\beta_{k}$ are truly equal to zero,
$\beta_{k}$ is correctly specified under pseudo-Huber constraints (up to
negligible differences between medians and pseudo-Huber centers). In addition,
while estimation will be valid in the presence of longitudinal or cluster-
dependent observations, it remains to extend our inferential methods to this
setting. Finally, it may be of interest to extend our method to address
settings in which Neyman-Scott-type problems (Neyman and Scott, 1948) are of
concern due to the presence of nuisance parameters whose dimension grows with
number of observations.

Open-source software implementing the method (including long and short-form
documentation) is available at github.com/statdivlab/radEmu. radEmu
abbreviates radEmuAbPill, which denotes “using relative abundance data to
estimate multiplicative differences in absolute abundances with partially
identified log-linear models”.

## References

  * Albert and Anderson (1984) Albert, A. and Anderson, J.A. (1984).  On the existence of maximum likelihood estimates in logistic regression models.  Biometrika 71(1), 1–10.
  * Barron (2019) Barron, J.T. (2019).  A general and adaptive robust loss function.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 4331–4339.
  * Benjamini and Hochberg (1995) Benjamini, Y. and Hochberg, Y. (1995).  Controlling the false discovery rate: a practical and powerful approach to multiple testing.  Journal of the Royal Statistical Society: Series B 57(1), 289–300.
  * Birch (1963) Birch, M. (1963).  Maximum likelihood in three-way contingency tables.  Journal of the Royal Statistical Society Series B: Statistical Methodology 25(1), 220–233.
  * Bull et al. (2002) Bull, S.B., Mak, C. and Greenwood, C.M. (2002).  A modified score function estimator for multinomial logistic regression in small samples.  Computational Statistics & Data Analysis 39(1), 57–74.
  * Clausen and Willis (2022) Clausen, D.S. and Willis, A.D. (2022).  Modeling complex measurement error in microbiome experiments.  arXiv preprint arXiv:2204.12733.
  * Davis et al. (2018) Davis, N.M. et al (2018).  Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data.  Microbiome 6(1), 1–14.
  * Fernandes et al. (2014) Fernandes, A.D. et al (2014).  Unifying the analysis of high-throughput sequencing datasets: characterizing rna-seq, 16s rrna gene sequencing and selective growth experiments by compositional data analysis.  Microbiome 2(1), 1–13.
  * Firth (1993) Firth, D. (1993).  Bias reduction of maximum likelihood estimates.  Biometrika 80(1), 27–38.
  * Guo et al. (2005) Guo, X. et al (2005).  Small-sample performance of the robust score test and its modifications in generalized estimating equations.  Statistics in Medicine 24(22), 3479–3495.
  * Jeffreys (1946) Jeffreys, H. (1946).  An invariant form for the prior probability in estimation problems.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 186(1007), 453–461.
  * Kaul et al. (2017) Kaul, A. et al (2017).  Analysis of microbiome data in the presence of excess zeros.  Frontiers in Microbiology 8, 2114.
  * Kosmidis and Firth (2011) Kosmidis, I. and Firth, D. (2011).  Multinomial logit bias reduction via the poisson log-linear model.  Biometrika 98(3), 755–759.
  * Li et al. (2021) Li, Z. et al (2021).  Ifaa: robust association identification and inference for absolute abundance in microbiome analyses.  Journal of the American Statistical Association 116(536), 1595–1608.
  * Lin and Peddada (2020) Lin, H. and Peddada, S.D. (2020).  Analysis of compositions of microbiomes with bias correction.  Nature Communications 11(1), 3514.
  * Love et al. (2014) Love, M.I., Huber, W. and Anders, S. (2014).  Moderated estimation of fold change and dispersion for rna-seq data with deseq2.  Genome Biology 15(12), 1–21.
  * Mandal et al. (2015) Mandal, S. et al (2015).  Analysis of composition of microbiomes: a novel method for studying microbial composition.  Microbial ecology in health and disease 26(1), 27663.
  * Martin et al. (2020) Martin, B.D., Witten, D. and Willis, A.D. (2020).  Modeling microbial abundances and dysbiosis with beta-binomial regression.  The Annals of Applied Statistics 14(1), 94.
  * McLaren et al. (2019) McLaren, M.R., Willis, A.D. and Callahan, B.J. (2019).  Consistent and correctable bias in metagenomic sequencing experiments.  eLife 8.
  * Milanese et al. (2019) Milanese, A. et al (2019).  Microbial abundance, activity and population genomic profiling with motus2.  Nature Communications 10(1), 1014.
  * Moffatt and Cookson (2017) Moffatt, M.F. and Cookson, W.O. (2017).  The lung microbiome in health and disease.  Clinical Medicine 17(6), 525.
  * Nearing et al. (2022) Nearing, J.T. et al (2022).  Microbiome differential abundance methods produce different results across 38 datasets.  Nature Communications 13(1), 1–16.
  * Neyman and Scott (1948) Neyman, J. and Scott, E.L. (1948).  Consistent estimates based on partially consistent observations.  Econometrica: Journal of the Econometric Society, 1–32.
  * Nocedal and Wright (1999) Nocedal, J. and Wright, S.J. (1999).  Numerical optimization.  Spinger.
  * Nocedal and Wright (2006) Nocedal, J. and Wright, S.J. (2006).  Numerical optimization.  Spinger.
  * Paulson et al. (2013) Paulson, J.N., Pop, M. and Bravo, H.C. (2013).  metagenomeseq: Statistical analysis for sparse high-throughput sequencing.  Bioconductor package 1(0), 191.
  * Robinson et al. (2010) Robinson, M.D., McCarthy, D.J. and Smyth, G.K. (2010).  edgeR: a bioconductor package for differential expression analysis of digital gene expression data.  Bioinformatics 26(1), 139–140.
  * Segata et al. (2011) Segata, N. et al (2011).  Metagenomic biomarker discovery and explanation.  Genome Biology 12(6), 1–18.
  * Sherman and Morrison (1950) Sherman, J. and Morrison, W.J. (1950).  Adjustment of an inverse matrix corresponding to a change in one element of a given matrix.  The Annals of Mathematical Statistics 21(1), 124–127.
  * Shreiner et al. (2015) Shreiner, A.B., Kao, J.Y. and Young, V.B. (2015).  The gut microbiome in health and in disease.  Current Opinion in Gastroenterology 31(1), 69.
  * White (1982) White, H. (1982).  Maximum likelihood estimation of misspecified models.  Econometrica, 1–25.
  * Williamson et al. (2022) Williamson, B.D., Hughes, J.P. and Willis, A.D. (2022).  A multiview model for relative and absolute microbial abundances.  Biometrics 78(3), 1181–1194.
  * Wirbel et al. (2019) Wirbel, J. et al (2019).  Meta-analysis of fecal metagenomes reveals global microbial signatures that are specific for colorectal cancer.  Nature Medicine 25(4), 679–689.
  * Yamashita and Takeshita (2017) Yamashita, Y. and Takeshita, T. (2017).  The oral microbiome and human health.  Journal of Oral Science 59(2), 201–206.
  * Young (2017) Young, V.B. (2017).  The role of the microbiome in human health and disease: an introduction for clinicians.  BMJ 356.

##  S1 Proof of Lemma 2.1

Lemma:  For arbitrary $\beta^{\dagger}\in\mathbb{R}^{p\times J}$, define
$G_{\beta^{\dagger}}=\\{\beta\in\mathbbm{R}^{p\times
J}:\beta=\beta^{\dagger}+\alpha\mathbf{1}_{j}^{T}\text{ for some
}\alpha\in\mathbb{R}^{p}\\}.$ If $\mathbf{X}$ has full column rank, for any
$\beta,\beta^{\star}\in\mathbb{R}^{p\times J}$ there exist
$\mathbf{z}\in\mathbb{R}^{n}$ and $\mathbf{z}^{\star}\in\mathbb{R}^{n}$ such
that under the model

$\displaystyle\log\mathbb{E}[Y_{i}|X_{i},\beta^{\star},z_{i},\delta]=z_{i}\mathbf{1}_{J}^{T}+X_{i}\beta^{\star}+\delta^{T},$
(14)

$\log\mathbb{E}[\mathbf{Y}|\mathbf{X},\beta,\mathbf{z}]=\log\mathbb{E}[\mathbf{Y}|\mathbf{X},\beta^{\star},\mathbf{z}^{\star}]$
if and only if $\beta^{\star}\in G_{\beta}$.

###### Proof.

Suppose $\beta^{\star}\in G_{\beta}$. Then there exists an
$\alpha\in\mathbb{R}^{p}$ such that
$\beta^{\star}=\beta+\alpha\mathbf{1}_{J}^{T}$ and by ((2), main text), for
any $\mathbf{z}$ and $\mathbf{z}^{\star}:=\mathbf{z}+\mathbf{X}\alpha$,
$\log\mathbb{E}[\mathbf{Y}|\mathbf{X},\beta,\mathbf{z}]=\log\mathbb{E}[\mathbf{Y}|\mathbf{X},\beta^{\star},\mathbf{z}^{\star}]$.
Now suppose $\beta^{\star}\not\in G_{\beta}$. If there exist
$\mathbf{z},\mathbf{z}^{\star}$ such that
$\log\mathbb{E}[\mathbf{Y}|\mathbf{X},\beta,\mathbf{z}]=\log\mathbb{E}[\mathbf{Y}|\mathbf{X},\beta^{\star},\mathbf{z}^{\star}]$,
we have

$\displaystyle\mathbf{z}\mathbf{1}_{J}^{T}+\mathbf{X}\beta=\mathbf{z}^{\star}\mathbf{1}_{J}^{T}+\mathbf{X}\beta^{\star}$
$\displaystyle\Rightarrow~{}$
$\displaystyle\mathbf{X}(\beta-\beta^{\star})=(\mathbf{z}^{\star}-\mathbf{z})\mathbf{1}_{J}^{T}$
$\displaystyle\Rightarrow~{}$
$\displaystyle\mathbf{X}^{T}\mathbf{X}(\beta-\beta^{\star})=\mathbf{X}^{T}(\mathbf{z}^{\star}-\mathbf{z})\mathbf{1}_{J}^{T}$
$\displaystyle\Rightarrow~{}$
$\displaystyle(\beta-\beta^{\star})=(\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}(\mathbf{z}^{\star}-\mathbf{z})\mathbf{1}_{J}^{T}$
$\displaystyle\Rightarrow~{}$
$\displaystyle\beta^{\star}=\beta+\tilde{\alpha}\mathbf{1}_{J}^{T}\text{ for
}\tilde{\alpha}:=(\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}(\mathbf{z}^{\star}-\mathbf{z})\in\mathbb{R}^{p}$
$\displaystyle\Rightarrow~{}$ $\displaystyle\beta^{\star}\in G_{\beta}$

which is a contradiction. Hence $\beta^{\star}\not\in G_{\beta}$ guarantees
that for no $\mathbf{z},\mathbf{z}^{\star}$ do we have
$\log\mathbb{E}[\mathbf{Y}|\mathbf{X},\beta,\mathbf{z}]=\log\mathbb{E}[\mathbf{Y}|\mathbf{X},\beta^{\star},\mathbf{z}^{\star}]$.
∎

##  S2 Pseudo-Huber Centering

For a generic vector $x=[x_{1},\dots,x_{J}]\in\mathbb{R}^{J}$ and smoothing
parameter $\delta>0$, we define the pseudo-Huber smoothed median $c^{\star}$
of $x$ as follows:

$\displaystyle c^{\star}$
$\displaystyle=\text{argmin}_{c}\sum_{j}\delta^{2}\Big{[}1+(\frac{x_{j}-c}{\delta})^{2}\Big{]}^{1/2}$

We now derive the form of $\frac{\partial}{\partial x_{j}}c^{\star}$, which is
needed for optimization and testing.

First, we have

$\displaystyle\frac{\partial}{\partial
c}\Big{[}\sum_{j}\delta^{2}\Big{[}1+(\frac{x_{j}-c}{\delta})^{2}\Big{]}^{1/2}\Big{]}|_{c=c^{\star}}=0$
$\displaystyle\Rightarrow~{}$
$\displaystyle\sum_{j}\frac{\delta^{2}}{2}\Big{[}1+(\frac{x_{j}-c}{\delta})^{2}\Big{]}^{-1/2}2(\frac{x_{j}-c}{\delta})(\frac{-1}{\delta})=0$
$\displaystyle\Rightarrow~{}$ $\displaystyle\sum_{j}w_{j}(x_{j}-c)=0\text{ for
}w_{j}:=\Big{[}1+(\frac{x_{j}-c}{\delta})^{2}\Big{]}^{-1/2}$

We also have

$\displaystyle\frac{\partial}{\partial x_{j^{\dagger}}}w_{j}$
$\displaystyle=-\frac{1}{2}\Big{[}1+(\frac{x_{j}-c}{\delta})^{2}\Big{]}^{-3/2}2(\frac{x_{j}-c}{\delta})(\frac{\mathbf{1}_{[j=j^{\dagger}]}-\frac{\partial}{\partial
x_{j^{\dagger}}}c}{\delta})$
$\displaystyle=-w_{j}^{3}(\frac{x_{j}-c}{\delta^{2}})\Big{[}\mathbf{1}_{[j=j^{\dagger}]}-\frac{\partial}{\partial
x_{j^{\dagger}}}c\Big{]}$

Before proceeding, we will derive a useful identity. Let
$y:=\frac{x_{j}-c}{\delta}$. Then we have

$\displaystyle w_{j}-w_{j}^{3}(\frac{x_{j}-c}{\delta})^{2}$
$\displaystyle=(1+y^{2})^{-1/2}-(1+y^{2})^{-3/2}y^{2}$
$\displaystyle=\frac{1}{(1+y^{2})^{1/2}}-\frac{y^{2}}{(1+y^{2})^{3/2}}$
$\displaystyle=\frac{1+y^{2}}{(1+y^{2})^{3/2}}-\frac{y^{2}}{(1+y^{2})^{3/2}}$
$\displaystyle=\frac{1}{(1+y^{2})^{3/2}}$ $\displaystyle=w_{j}^{3}$

We can now use implicit differentiation to find $\frac{\partial}{\partial
x_{j}}c^{\star}$

$\displaystyle\frac{\partial}{\partial x_{j}}\sum_{j}w_{j}(x_{j}-c^{\star})$
$\displaystyle=$
$\displaystyle\sum_{j}\Bigg{[}w_{j}(\mathbf{1}_{[j=j^{\dagger}]}-\frac{\partial}{\partial
x_{j^{\dagger}}}c^{\star})-w_{j}^{3}(\frac{x_{j}-c^{\star}}{\delta^{2}})\Big{[}\mathbf{1}_{[j=j^{\dagger}]}-\frac{\partial}{\partial
x_{j^{\dagger}}}c^{\star}\Big{]}\Bigg{]}=0$ $\displaystyle\Rightarrow$
$\displaystyle\sum_{j}\Big{[}w_{j}-w_{j}^{3}(\frac{x_{j}-c^{\star}}{\delta})^{2}\Big{]}\frac{\partial}{\partial
x_{j^{\dagger}}}c^{\star}=w_{j^{\dagger}}-w_{j^{\dagger}}^{3}(\frac{x_{j^{\dagger}}-c^{\star}}{\delta})^{2}$
$\displaystyle\Rightarrow$
$\displaystyle\sum_{j}w_{j}^{3}\frac{\partial}{\partial
x_{j^{\dagger}}}c^{\star}=w_{j^{\dagger}}^{3}\text{ by the identity above }$
$\displaystyle\Rightarrow$ $\displaystyle\frac{\partial}{\partial
x_{j^{\dagger}}}c^{\star}=\frac{w_{j^{\dagger}}^{3}}{\sum_{j}w_{j}^{3}}$

For a given $x$, to find the pseudo-Huber center $c^{\star}(x)=\text{arg
min}_{c}\sum_{j}\delta^{2}\Big{[}1+(\frac{x_{j}-c}{\delta})^{2}\Big{]}^{1/2}$
with parameter $\delta$ by minimizing successive quadratic approximations to
$f(c;x,\delta):=\sum_{j}\delta^{2}\Big{[}1+(\frac{x_{j}-c}{\delta})^{2}\Big{]}^{1/2}$.
For current value $c^{(t)}$ of $c$, we compute, for $j=1,\dots,J$,
$w_{j}^{(t)}:=\Big{[}1+(\frac{x_{j}-c^{(t)}}{\delta})^{2}\Big{]}^{-1/2}$.
Letting
$h^{(t)}(c)=\sum_{j}\frac{1}{2}w_{j}^{(t)}(\frac{x_{j}-c}{\delta})^{2}$, we
then construct quadratic approximation
$\tilde{f}^{(t)}:=f(c^{(t)})-h^{(t)}(c^{(t)})+h^{(t)}(c)$ to $f$ at $c^{(t)}$.
We update $c^{(t+1)}$ as the minimizer of $\tilde{f}^{(t)}$, which is
conveniently the weighted mean of $x$ with weights
$\frac{1}{\sum_{j}w_{j}}(w_{1},\dots,w_{J})$.

##  S3 Algorithms

###  S3.1 Unconstrained optimization

Algorithm 1 Coordinate Descent for Maximum Likelihood Estimation

  1. 1.

Initiate with

     * •

Data $(X,Y)$

     * •

Starting value $\beta^{(0)}$ for $\beta$

     * •

Identifiability constraint $g(\beta_{k})=0$ for $k=1,\dots,p$

     * •

Convergence tolerance $\epsilon>0$

     * •

Iteration limit $t_{\text{max}}$

     * •

Choice of working identifiability constraint $\beta^{j^{\dagger}}=0$ (e.g.,
$j^{\dagger}$ such that $\sum_{i}\mathbf{1}_{Y_{ij^{\dagger}}>0}$ is
maximized)

  2. 2.

For $k=1,\dots,p$, impose working constraint by setting
$\beta_{k}=\beta_{k}-\beta_{kj^{\dagger}}$

  3. 3.

For iteration $t=1,\dots,t_{\text{max}}$

    1. (a)

For $j\in\\{1,\dots,J\\}/\\{j^{\dagger}\\}$, update $\beta^{j}$, the $j$-th
column of $\beta$ with single iteration of algorithm (2)

    2. (b)

Update $z_{i}$ for $i=1,\dots,n$ via equation (7).

    3. (c)

If $\text{max}_{k,j}|\beta^{(t)}_{kj}-\beta^{(t-1)}_{kj}|<\epsilon$ or if
$t=t_{\text{max}}$

       * •

For $k=1,\dots,p$ enforce original identifiability constraint by setting
$\beta_{k}=\beta_{k}^{(t)}-g(\beta_{k}^{(t)})$.

       * •

Exit algorithm and return $\beta$.

    4. (d)

Otherwise increment $t$ by $1$ and return to step 3a.

Algorithm 2 Scoring Update with Line Search

  1. 1.

Initiate with

     * •

Data $(X,Y^{j})$

     * •

Current values of $\beta^{j}$ and $\mathbf{z}$

     * •

Maximum step size $s$

  2. 2.

Then

    1. (a)

Compute

       * •

gradient $S^{j}=X^{T}(Y^{j}-\boldsymbol{\mu}^{j}))$ for
$\boldsymbol{\mu^{j}}=[\mu_{1}^{j},\dots,\mu_{n}^{j}]^{T}=[\exp(X_{1}\beta^{j}+z_{1}),\dots,\exp(X_{n}\beta^{j}+z_{n})]^{T}$

       * •

information $I^{j}=X^{T}WX$ for $W$ a diagonal matrix with $i$-th diagonal
entry $\mu^{j}_{i}$

    2. (b)

Compute update direction $d^{j}={I^{j}}^{-1}S^{j}$. If $I^{j}$ is poorly
conditioned, attempt to compute update direction as
$(I^{j}+\sigma\mathbf{I})^{-1}S^{j}$, sequentially increasing $\sigma$ until
$(I^{j}+\sigma\mathbf{I})^{-1}$ is numerically invertible. ($I^{j}$ is
theoretically always full rank if $X$ is full rank; in practice, it may be
poorly conditioned numerically.)

    3. (c)

Perform line search to find update $\beta^{j}+\alpha d$ satisfying
$l(\beta^{j}+\alpha d)\geq l(\beta^{j})+c\alpha d^{T}\nabla l$ (i.e., the
Armijo rule) where $l$ is the log likelihood function evaluated in $\beta^{j}$
with $z$ held fixed, $c>0$ is a constant (by default we use $c=10^{-4}$ as
suggested by Nocedal and Wright (1999, Ch. 17), $\nabla l$ is the partial
derivative of the log likelihood with respect to $\beta^{j}$, and $\alpha$ is
initiated at $0.5$ and halved until the Armijo condition is met. After the
Armijo condition is met, optionally rescale $\alpha d$ so that
$\max_{k}|\alpha d_{k}|\leq\xi$ for some maximum step length $\xi$ before
accepting update $\beta^{j}+\alpha d$ to prevent large steps that may lead to
ill conditioning. We by default use $\xi=1$.

Algorithm 3 Penalized Likelihood Estimation via Iterative Maximum Likelihood

  1. 1.

Initiate with

     * •

Data $(X,Y)$

     * •

Starting value $\beta^{(0)}$ for $\beta$

     * •

Identifiability constraint $g(\beta_{k})=0$ for $k=1,\dots,p$

     * •

Convergence tolerance $\epsilon>0$

     * •

Iteration limit $t_{\text{max}}$

  2. 2.

Set $Y^{(1)}=\mathbf{Y}+\delta$ for some $\delta\geq 0$. (Choosing $\delta>0$
is not necessary but can ameliorate slow convergence issues that otherwise
might arise in the first iteration due to separation.)

  3. 3.

For iteration $t=1,\dots,t_{\text{max}}$

    1. (a)

Update $\beta^{(t)}$ using algorithm 2 applied to data $\mathbf{Y}^{(t)}$ with
other arguments as specified above

    2. (b)

Update $\mathbf{Y}^{(t+1)}$ using $\beta^{(t)}$ and $z$ as laid out in Section
3 of main text

    3. (c)

If $\text{max}_{k,j}|\beta_{kj}^{(t)}-\beta_{kj}^{(t-1)}|<\epsilon$ or if
$t=t_{\text{max}}$, exit and return $\beta$. Otherwise increment $t$ and
return to step 3a.

###  S3.2 Constrained optimization

Algorithm 4 Constrained Estimation

  1. 1.

Initiate with

     * •

Covariate matrix $\mathbf{X}$

     * •

Initial value of $\beta$ and augmented outcomes $\tilde{\mathbf{Y}}$ obtained
via fitting an unconstrained penalized model

     * •

Constraint function $g$

     * •

Row and column indexes $k^{\star},j^{\star}$ for the element of
$\boldsymbol{\beta}$ to be constrained

     * •

Initial value of augmented Lagrangian penalty parameter $\rho$

     * •

Values $\kappa\in(0,1)$ and $\tau>1$ controlling scaling of $\rho$ in each
iteration

     * •

Convergence tolerance $\epsilon$

     * •

Constraint tolerance $\gamma$

     * •

Iteration limits $t_{max}$ for outer loop and $t^{inner}_{max}$ for inner loop

  2. 2.

Set $u^{(1)}=\rho(\beta_{k^{\star}j^{\star}}-g(\beta_{k^{\star}}))$ and
$\rho^{(1)}=\rho$

  3. 3.

Then for iteration $t=1,\dots,t_{max}$

    1. (a)

Define
$\mathcal{L}^{(t)}=-l(\beta,z)+u^{(t)}\big{[}\beta_{k^{\star}j^{\star}}-g(\beta_{k^{\star}})\big{]}+(\rho^{(t)}/2)\big{[}\beta_{k^{\star}j^{\star}}-g(\beta_{k^{\star}})\big{]}^{2}$

    2. (b)

Via block coordinate descent in $\beta$ and $z$, obtain
$(\beta^{(t)},z^{(t)})$ satisfying $||\nabla\mathcal{L}^{(t)}||_{2}<\epsilon$

       * •

Coordinate descent steps in $\beta$ are approximate Newton steps using
$-\nabla^{2}l+\rho\nabla g\nabla^{T}g$ as an approximation to the Hessian of
the augmented Lagrangian

       * •

Coordinate descent steps in $z$ are available in closed form

    3. (c)

If $R^{(t)}:=\beta^{(t)}_{k^{\star}j^{\star}}-g(\beta^{(t)}_{k^{\star}})$
satisfies $|R^{(t)}|<\gamma$ or $t=t_{max}$, exit and return $\beta^{(t)}$

    4. (d)

Otherwise, set $\rho^{(t+1)}=\tau\rho^{(t)}$ and
$u^{(t+1)}=\rho^{(t+1)}R^{(t)}$

###  S3.3 Additional Details: Efficient approximate Newton method

As discussed in the main text (Section 3.2), we approximate the Hessian of the
augmented Lagrangian induced by our constrained optimization problem as
$-\nabla^{2}l+\rho(\nabla
g-\vec{e}_{[k^{\star}j^{\star}]})(\nabla^{T}g-\vec{e}_{[k^{\star}j^{\star}]})$
and compute its inverse as a rank-one update to $-(\nabla^{2}l)^{-1}$ via the
Sherman-Morrison identity.

We may justify this choice on several grounds. When $\rho$ is small in early
iterations of our algorithm, $u$ will also be small, and so provided $g$ is
not too large (which in practice it is not), $-\nabla^{2}l$ will well
approximate $\nabla^{2}\mathcal{L}$. In later iterations, we expect to have
$g$ close to zero and $\rho\gg|u|$, so $(u+\rho g)\nabla^{2}g$ will in general
be small relative to $\rho\nabla g\nabla^{T}g$. Of course, these are somewhat
loose justifications, and we may also observe that, as the sum of a positive
definite and a positive semi-definite matrix, $-\nabla^{2}l+\rho\nabla
g\nabla^{T}g$ is positive definite and hence a Newton update using this
approximation to the Hessian will always take us in a descent direction.

##  S4 Type 1 Error & Power: Additional Results

Distribution | $J$ | $n=10$ | $n=50$ | $n=250$
---|---|---|---|---
Poisson | 10 | 0.061 | 0.051 | 0.046
Poisson | 50 | 0.056 | 0.050 | 0.051
Poisson | 250 | 0.049 | 0.056 | 0.050
ZINB | 10 | 0.028 | 0.026 | 0.032
ZINB | 50 | 0.032 | 0.039 | 0.043
ZINB | 250 | 0.025 | 0.043 | 0.044
Table S1: Empirical type 1 error for the robust score test at the $0.05$ level by simulated distribution, number of taxa $J$, and sample size $n$. Results are shown for 10,000 simulations.  Distribution | $J$ | $n=10$ | $n=50$ | $n=250$
---|---|---|---|---
Poisson | 10 | 0.107 | 0.045 | 0.042
Poisson | 50 | 0.159 | 0.069 | 0.053
Poisson | 250 | 0.167 | 0.080 | 0.057
ZINB | 10 | 0.252 | 0.058 | 0.034
ZINB | 50 | 0.293 | 0.093 | 0.057
ZINB | 250 | 0.307 | 0.098 | 0.059
Table S2: Empirical type 1 error for robust Wald test at the $0.05$ level by simulated distribution, number of taxa $J$, and sample size $n$. Results are shown for 10,000 simulations.  Distribution | $J$ | $n=10$ | $n=50$ | $n=250$
---|---|---|---|---
Poisson | 10 | 0.918 | 1.000 | 1.000
Poisson | 50 | 0.718 | 1.000 | 1.000
Poisson | 250 | 0.640 | 1.000 | 1.000
ZINB | 10 | 0.034 | 0.138 | 0.770
ZINB | 50 | 0.040 | 0.180 | 0.872
ZINB | 250 | 0.028 | 0.194 | 0.908
Table S3: Empirical power for detection of smaller magnitude signal ($\beta_{kj}-g(\beta_{k})=1$) for robust score test at the $0.05$ level by simulated distribution, number of taxa $J$, and sample size $n$. Results are shown for 500 simulations.  Distribution | $J$ | $n=10$ | $n=50$ | $n=250$
---|---|---|---|---
Poisson | 10 | 0.996 | 1.000 | 1.000
Poisson | 50 | 0.996 | 1.000 | 1.000
Poisson | 250 | 0.992 | 1.000 | 1.000
ZINB | 10 | 0.284 | 0.276 | 0.814
ZINB | 50 | 0.366 | 0.380 | 0.898
ZINB | 250 | 0.410 | 0.456 | 0.938
Table S4: Empirical power for detection of smaller magnitude signal ($\beta_{kj}-g(\beta_{k})=1$) for robust Wald test at the $0.05$ level by simulated distribution, number of taxa $J$, and sample size $n$. Results are shown for 500 simulations.  Distribution | $J$ | $n=10$ | $n=50$ | $n=250$
---|---|---|---|---
Poisson | 10 | 0.882 | 1.000 | 1.000
Poisson | 50 | 0.986 | 1.000 | 1.000
Poisson | 250 | 0.920 | 1.000 | 1.000
ZINB | 10 | 0.096 | 0.990 | 1.000
ZINB | 50 | 0.386 | 0.986 | 1.000
ZINB | 250 | 0.276 | 0.922 | 1.000
Table S5: Empirical power for detection of larger magnitude signal ($\beta_{kj}-g(\beta_{k})=5$) for robust score test at the $0.05$ level by simulated distribution, number of taxa $J$, and sample size $n$. Results are shown for 500 simulations.  Distribution | $J$ | $n=10$ | $n=50$ | $n=250$
---|---|---|---|---
Poisson | 10 | 1.000 | 1.000 | 1.000
Poisson | 50 | 1.000 | 1.000 | 1.000
Poisson | 250 | 1.000 | 1.000 | 1.000
ZINB | 10 | 0.920 | 1.000 | 1.000
ZINB | 50 | 0.924 | 1.000 | 1.000
ZINB | 250 | 0.962 | 1.000 | 1.000
Table S6: Empirical power for detection of larger magnitude signal
($\beta_{kj}-g(\beta_{k})=5$) over $500$ simulations for robust Wald test at
the $0.05$ level by simulated distribution, number of taxa $J$, and sample
size $n$. Results are shown for 500 simulations.

##  S5 Additional figures

Figure S7: Empirical power to reject $\beta_{kj}=g(\beta_{k})$ at the $0.05$
level for zero-inflated Negative Binomial draws. Results are shown for 100
simulations. The primary purpose of these results is to provide investigators
with power estimates, that can then be refined further using our simulation
code.