# A Unified Bayesian Framework for Modeling Measurement Error in Multinomial
Data

Matthew D<EMAIL_ADDRESS>[ Andee
Kaplanlabel=e3]<EMAIL_ADDRESS>[ Mevin
<EMAIL_ADDRESS>[ Department of Statistics,
Colorado State University, Fort Collins, CO, USA.  Department of Statistics,
Colorado State University, Fort Collins, CO, USA.  Department of Statistics
and Data Sciences, The University of Texas, Austin, TX, USA.

(2023+)

###### Abstract

Measurement error in multinomial data is a well-known and well-studied
inferential problem that is encountered in many fields, including engineering,
biomedical and omics research, ecology, finance, and social sciences.
Surprisingly, methods developed to accommodate measurement error in
multinomial data are typically equipped to handle false negatives or false
positives, but not both. We provide a unified framework for accommodating both
forms of measurement error using a Bayesian hierarchical approach. We
demonstrate the proposed method’s performance on simulated data and apply it
to acoustic bat monitoring data.

categorical data,

ecology,

misclassification,

record linkage,

zero-inflation,

###### keywords:

††volume: TBA††issue: TBA

, , and

##  1 Introduction

Measurement error in multinomial data is a well-known and well-studied
inferential problem that is encountered in many fields, including engineering,
biomedical and omics research, ecology, finance, and social sciences (Swartz
et al.,, 2004; Pérez et al.,, 2007; Molinari,, 2008; Datta et al.,, 2021;
Mulick et al.,, 2022). In this work, we consider two types of measurement
error defined as (1) false negatives that occur when a particular category or
class is present in the population but it is not observed in the sample and
(2) false positives that occur when a sampled observation is misclassified
into the wrong category. While both types of measurement error may be present
in the data, existing methods are not designed to accommodate them
simultaneously which may bias inference. Motivated by multispecies wildlife
monitoring data collected in ecological research, we propose a unified
framework for accommodating both forms of measurement error when modeling
multinomial data. Our approach differs from existing methods in that it
explicitly models the probability of misclassification for each observation,
accommodates individual-level covariates associated with the probability of
being a true/false negative, the relative abundances of the multinomial count
data, as well as the probability of misclassification, and is scalable to
high-dimensional classification problems.

Modeling false negatives in multinomial data is closely related to the notion
of handling zero-inflation in univariate and multivariate count data. See
Blasco-Moreno et al., (2019) for an in-depth discussion of zero counts in the
context of ecological research studies and Scharf et al., (2022) for an
overview of hierarchical models for occupancy data. Count data are considered
zero-inflated when the number of zeroes observed in the data set is larger
than expected under the assumptions of the sampling distribution. Zero-
inflated count models are typically constructed as a two-component mixture of
a point mass at zero and a sampling distribution for the count data (e.g.,
Poisson or negative binomial distributions in the univariate setting) (Xu et
al.,, 2015; Zhang and Yi,, 2020; Jiang et al.,, 2021; Shuler et al.,, 2021).
To achieve this, an at-risk indicator is introduced into the model to
differentiate between at-risk zeros (i.e., a zero count is observed, but there
is a positive probability of occurrence) and structural zeros (i.e., a zero
count is observed because there is zero probability of occurrence) (Neelon,,
2019). Mapping these to the measurement error definitions above, at-risk
(structural) zeros are synonymous with false (true) negatives. In multivariate
settings, researchers link zero-inflated univariate count models via latent
parameters that control the dependence structure between counts (Aitchison and
Ho,, 1989; Chiquet et al.,, 2021). This approach models the multivariate
counts unconditionally on the total count for the sample and is therefore not
designed for multinomial classification problems or multivariate compositional
count data when the total number of counts is fixed. Koslovsky, (2023)
introduced a zero-inflated Dirichlet-multinomial (DM) distribution for
handling excess zeros in multivariate compositional count data, which differs
from traditional approaches for modeling zero-inflation (or false negatives)
in count data by assuming a mixture distribution on the count probabilities as
opposed to the sampling distribution. Using a combination of data augmentation
strategies, their approach is scalable to large compositional spaces, can
accommodate covariates associated with zero-inflation and relative abundances,
and has shown promising estimation performance in simulation.

Existing methods designed to model false positives in multinomial data
typically assume the observed classifications follow a multinomial
distribution given the true (latent) classification (Swartz et al.,, 2004;
Wang et al.,, 2020; Pérez et al.,, 2007; Frénay and Verleysen,, 2013). With
this approach, the number of rows for the resulting matrix of
misclassification probabilities equals the number of true categories, and the
number of columns equals the number of observed classes with each row summing
to one. The task of modeling false positives in multinomial data draws
parallels to a popular approach for entity resolution. Entity resolution is
the process of resolving duplicates in many overlapping data sets without the
benefit of a unique identifying attribute. In the hit-miss approach to entity
resolution, observed records are assumed to either represent the true records
associated with an entity (hit) or a distorted version of this truth (miss)
(Tancredi and Liseo,, 2011; Copas and Hilton,, 1990). These potentially noisy
records are then directly modeled using a mixture model in which two records
that are associated with the same latent truth refer to the same entity and
can be duplicated (Steorts et al.,, 2016). This direct approach of modeling
measurement error in the likelihood is an analog to the false positive
misclassification problem we are interested in addressing in that the observed
classifications can either be the true value (hit) or a distorted version of
that truth (miss). Potential benefits of directly modeling the distortion
process in a hit-miss framework include the ability to choose any appropriate
distribution for the misclassification, incorporation of expert knowledge into
the model via the priors on misclassification probabilities, and inference on
the probability of misclassification after obtaining data.

As mentioned previously, this work was motivated by multispecies occupancy
modeling in ecological research. The goal of occupancy modeling in ecological
research studies is to draw inference on true occurrence given a set of
observations that are subject to measurement error (i.e., imperfect
detection). Imperfect detection typically occurs in two different ways; (1) a
species may go undetected (i.e., false negative) and (2) an observed
individual may be misclassified (i.e., false positive). Even with increased
sampling effort, imperfect detection may still occur, resulting in biased
inference if ignored when modeling (Kellner and Swihart,, 2014).

Historically, statistical methods developed to handle imperfect detection have
focused on false negatives (Bayley and Peterson,, 2001; MacKenzie et al.,,
2002; Royle and Nichols,, 2003; MacKenzie et al.,, 2003; Tyre et al.,, 2003;
Broms et al.,, 2015; Dorazio et al.,, 2006, 2011; Devarajan et al.,, 2020).
However, more recently researchers have proposed methods that account for both
false positives and negatives, in part due to the emergence of automated
species detection methods (e.g., unmanned aerial systems and automated
recording units) and volunteer-based surveys (e.g., citizen science) for
monitoring wildlife populations. When developing single- or multispecies
(community) occupancy models, researchers typically take a hierarchical
approach, often referred to as occupancy-detection models, which jointly model
the ecological and observation (or detection) process. This technique allows
researchers to differentiate between latent species occupancy and observed
species detection and effectively account for potential false negatives.
Typically, this is achieved by introducing a site- or location-specific latent
species indicator that models whether or not a species is present. If the
species is present (absent) at that site, there is a positive (zero)
probability of detecting it.

Methods that handle potential species misclassification in occupancy modeling
were initially developed for single species studies (Royle and Link,, 2006;
Miller et al.,, 2011; Chambert et al.,, 2015; Ruiz-Gutierrez et al.,, 2016;
Chambert et al., 2018b, ). Chambert et al., 2018a  introduced a two-species
occupancy model that accounts for both species misidentification and non-
detection. Their approach is based on the premise that false detections for a
given species occur due to the misidentification with a closely related
species. Recently, Wright et al., (2020) developed a multispecies occupancy
model that handles both forms of measurement error for two or more species
aggregately at each site visit. By assuming (1) the true count of each species
follows a Poisson distribution given it is present at the site visit, (2) the
number of detections for each species follows a multinomial distribution given
the true species counts, and (3) the detection counts are independent across
species, the authors demonstrate how the observed/detected counts can be
directly modeled without conditioning on the true count for each species.
Spiers et al., (2022) developed a multispecies occupancy model similar to
Wright et al., (2020) that accommodates individual-level validation data (as
opposed to site-level) which allows for more flexibility when modeling
heterogeneity with covariates and morphospecies. This class of models has
shown promising estimation performance on simulated data that are generated
consistently with the model assumptions. However, it is not clear how these
methods would perform under alternative data generation processes (e.g.,
overdispersed species abundances) commonly encountered in ecological studies
(Lindén and Mäntyniemi,, 2011).

A fundamental issue shared by any model designed to accommodate
misclassification is that the model is not identifiable without additional
information about the zero-inflation or (mis)classification process beyond the
raw data (Swartz et al.,, 2004). Occupancy-detection methods that accommodate
false positives typically deal with identifiability using informative priors
in Bayesian settings, auxiliary/calibration data to estimate the matrix of
misclassification probabilities (typically referred to as a confusion matrix
in ecological research settings) separately from the abundance model, and
validating the true species of a subset of the data (often referred to as
unambiguous/ambiguous detections in coupled classification models from site
confirmation designs) (Wright et al.,, 2020; Chambert et al.,, 2015; Guillera-
Arroita et al.,, 2017). Stratton et al., (2022) explore these strategies
rigorously in simulation. Swartz et al., (2004) provide an extensive
discussion of identifiability issues in multinomial classification models and
propose using constraints to break the symmetry of the model, similar to what
is done to accommodate label switching in Bayesian mixture models (Jasra et
al.,, 2005).

In this work, we propose a novel method for simultaneously modeling false
positives and false negatives in multinomial data. Using a zero-inflated DM
distribution, we accommodate potential false negatives in the underlying true
classification as well as potential overdispersion. We then introduce a latent
hit-miss indicator to model misclassification which allows our approach to
differentiate between true detections and detection by chance. Applied to
wildlife monitoring data, our method can be thought of as a flexible
alternative to existing multispecies occupancy-detection models that forgoes
attempting to correctly model the ecological process in favor of simply
modeling the true classifications. The model is developed at the individual
level which allows for more granular covariate information to model
heterogeneity and accommodate potential morphospecies. We assume a
hierarchical structure for the concentration parameters of the true
classification which allows borrowing of information across sites and/or
visits, an important feature for efficient multispecies occupancy modeling
(Iknayan et al.,, 2014).

##  2 Methods

In this section, we present a general formulation for modeling measurement
error in multinomial data, making connections to relevant occupancy-detection
modeling aspects as necessary. Let the $C$-dimensional vector $\bm{y}_{ijl}$
represent the observed classification for the $i^{th}$ $(i=1,\dots,N$)
observation (or site/location) at the $j^{th}$ $(j=1,\dots,n_{i})$ measurement
(or visit) for the $l^{th}$ $(l=1,\dots,L_{ij})$ individual (or organism),
where $y_{ijlc}=1$ indicates the observed individual was classified into the
$c^{th}$ category (species) and 0 otherwise. In general, the model does not
require $n_{i}>1$, however in various fields, including ecological monitoring,
each site is typically visited multiple times to improve inference (MacKenzie
et al.,, 2002; Lele et al.,, 2012). We let the $T$-dimensional vector
$\bm{z}_{ijl}$ represent the individual’s true classification (or species),
where $z_{ijlt}=1$ indicates the individual truly belongs to the $t^{th}$
category and 0 otherwise. For ease of presentation, we assume $T=C$ and that
the ordering of the elements is the same in $\bm{y}_{ijl}$ and $\bm{z}_{ijl}$
(i.e., $t=c$ corresponds to the same category).

For each observed individual, we introduce a latent hit-miss or
misclassification indicator $\tau_{ijl}\in\\{0,1\\}$, where $0$ indicates that
$\bm{y}_{ijl}=\bm{z}_{ijl}$. To model the observed classifications, we assume

$\bm{y}_{ijl}|\bm{\theta}_{t},\tau_{ijl},\bm{z}_{ijl}\sim\tau_{ijl}\mbox{Multinomial}(1,\bm{\theta}_{t})+(1-\tau_{ijl})\delta_{\bm{z}_{ijl}}(\bm{y}_{ijl}),$
(1)

where $\bm{\theta}_{t}$ is a $C$-dimensional vector of observed classification
probabilities for the ${t}^{th}$ true classification and $\delta_{w}(\cdot)$
is a Dirac delta function at $w$. With this formulation, we assume that if
there is no misclassification (i.e., $\tau_{ijl}=0$), then
$\bm{y}_{ijl}=\bm{z}_{ijl}$, otherwise, the observed individual is considered
misclassified with $\tau_{ijl}=1$. Note that this approach allows for
$\bm{y}_{ijl}$ to be misclassified into the correct category (i.e.,
$\bm{y}_{ijl}=\bm{z}_{ijl}$, but $\tau_{ijl}=1$). As such, our modeling
approach places a positive probability of a “lucky guess” to occur. In section
3, we discuss how to restrict the model to prevent this from occurring if
desired.

We assume the classification probabilities depend on the true classification
of each individual with $\bm{\theta}_{t}\sim\mbox{Dirichlet}(\bm{\nu}_{t})$,
where $\bm{\nu}_{t}$ is a $C$-dimensional vector of concentration
hyperparameters. To allow the classification probabilities to depend on an
observed set of covariates, $\nu_{tc}$ can be replaced with a log-linear
regression model similar to Wadsworth et al., (2017). Next, we let the latent
misclassification indicators,

$\tau_{ijl}|\bm{z}_{ijl},\bm{\beta}_{\psi_{t}},\bm{x}_{ijl}\sim\mbox{Bernoulli}(\psi_{tijl}),$
(2)

where $\text{logit}(\psi_{tijl})=\bm{x}_{ijl}^{\prime}\bm{\beta}_{\psi_{t}}$,
$\bm{x}_{ijl}$ is a $P_{\psi}$-dimensional vector of observed covariates that
are observation-, measurement-, and/or individual-specific (including an
intercept term), and $\bm{\beta}_{\psi_{t}}$ are the corresponding regression
coefficients. We assume
$\beta_{\psi_{tp}}\sim\mbox{Normal}(\mu_{\psi},\sigma_{\psi}^{2})$. Note that
the covariate effects on misclassification are allowed to vary based on the
true classification of the individual.

Next we model the true classification of each observation

$\bm{z}_{ijl}|\bm{\Theta}_{ij}\sim\text{Multinomial}(1,\bm{\Theta}_{ij}),$ (3)

where $\bm{\Theta}_{ij}$ is a $T$-dimensional vector of true classification
probabilities, which we assume follows a Dirichlet$(\bm{\gamma}_{ij})$ with
$\bm{\gamma}_{ij}$ a $T$-dimensional vector of concentration hyperparameters.
With the availability of validation data, some of the $\bm{z}_{ijl}$ will be
known and fixed in the model to inform the estimation of the classification
matrix
$\bm{\theta}=(\bm{\theta}_{1}^{\prime},\dots,\bm{\theta}_{T}^{\prime})^{\prime}$,
similar to Wright et al., (2020) and Spiers et al., (2022). We refer to
$\bm{\theta}$ as a classification matrix to differentiate it from methods that
use a “confusion matrix” to model false positives in multinomial data which,
unlike our model, do not explicitly model misclassfication.

We can equivalently model $\bm{\Theta}_{ij}$ as a set of independent gamma
random variables normalized by their sum (i.e.,
$\bm{z}_{ijl}|\bm{\alpha}_{ij}\sim\text{Multinomial}(1,\frac{\alpha_{ijt}}{\bar{\alpha}_{ij}})$,
where $\alpha_{ijt}\sim\text{Gamma}(\gamma_{ijt},1)$ and
$\bar{\alpha}_{ij}=\sum_{t=1}^{T}\alpha_{ijt}$). This reparameterization
enables us to account for potential zero-inflation (i.e., false negatives or
non-detection) by introducing a latent at-risk (or occupancy) indicator
$\zeta_{ijt}$ for the $t^{th}$ category at the $i^{th}$ observation and
$j^{th}$ measurement. Specifically, we instead let

$\alpha_{ijt}|\zeta_{ijt},\gamma_{ijt}\sim\zeta_{ijt}\mbox{Gamma}(\gamma_{ijt},1)+(1-\zeta_{ijt})\delta_{0}(\alpha_{ijt}),$
(4)

similar to Koslovsky, (2023). To model the latent at-risk indicators, we
assume
$\zeta_{ijt}|\bm{\beta}_{\eta_{t}},\bm{x}_{i}\sim\mbox{Bernoulli}(\eta_{it})$,
where $\text{logit}(\eta_{it})=\bm{x}_{i}^{\prime}\bm{\beta}_{\eta_{t}}$,
$\bm{x}_{i}$ is a $P_{\eta}$-dimensional set of observation-specific
covariates (including an intercept term), and $\bm{\beta}_{\eta_{t}}$
represent the corresponding true classification-specific regression
coefficients. We then let
$\beta_{\eta_{tp}}\sim\mbox{Normal}(\mu_{\eta},\sigma^{2}_{\eta})$. To allow
the relative abundances to depend on a set of covariates, we set
$\log(\gamma_{ijt})=\bm{x}_{ij}^{\prime}\bm{\beta}_{\gamma_{t}}$ with
$\beta_{\gamma_{tp}}\sim\mbox{Normal}(\mu_{\gamma},\sigma^{2}_{\gamma})$ and
$\bm{x}_{ij}$ an observation- and/or measurement-specific set of covariates.

##  3 Posterior Sampling and Inference

For posterior inference, we construct a Metropolis-Hastings within Gibbs
sampler. The full joint distribution is defined as

$\begin{split}&\prod_{i=1}^{N}\prod_{j=1}^{n_{i}}\prod_{l=1}^{L_{ij}}p(\bm{y}_{ijl}|\bm{\theta}_{t},\tau_{ijl},\bm{z}_{ijl})p(\bm{z}_{ijl}|\bm{\Theta}_{ij})p(\tau_{ijl}|\bm{z}_{ijl},\bm{\beta}_{\psi_{t}},\bm{x}_{ijl})p(\omega_{\tau_{ijl}})\\\
&\times\prod_{i=1}^{N}\prod_{j=1}^{n_{i}}\prod_{t=1}^{T}p(\alpha_{ijt}|\zeta_{ijt},\bm{\beta}_{\gamma_{t}},\bm{x}_{ij})p(\zeta_{ijt}|\bm{\beta}_{\eta_{t}},\bm{x}_{i})p(\omega_{\zeta_{tij}})\\\
&\times\prod_{i=1}^{N}\prod_{j=1}^{n_{i}}p(\mu_{ij}|\bar{\alpha}_{ij})\prod_{t=1}^{T}\left[p(\bm{\beta}_{\eta_{t}})p(\bm{\beta}_{\psi_{t}})p(\bm{\beta}_{\gamma_{t}})p(u_{t}|\bar{a}_{t})\prod_{c=1}^{C}p(a_{tc})\right],\end{split}$
(5)

where we introduce an auxiliary parameter
$\mu_{ij}|\bar{\alpha}_{ij}\sim\mbox{Gamma}(1,\bar{\alpha}_{ij})$ for
efficient sampling of $\alpha_{ijt}$. Similarly, we reparameterize
$\theta_{tc}=a_{tc}/\bar{a}_{t}$ and assume
$a_{tc}\sim\mbox{Gamma}(\nu_{tc},1)$ with auxiliary parameter
$u_{t}\sim\mbox{Gamma}(1,\bar{a}_{t})$ and $\bar{a}_{t}=\sum_{c=1}^{C}a_{tc}$.
In addition to enabling efficient sampling of $\bm{\theta}_{t}$, this step
provides the opportunity to easily incorporate covariates and restrict the
model to disallow correct classifications by chance by fixing $a_{tt}=0$.
Additionally, we introduce a latent set of auxiliary parameters
$\omega_{\zeta_{tij}}\sim\mbox{PG}(1,0)$ and
$\omega_{\tau_{ijl}}\sim\mbox{PG}(1,0)$ following Polson et al., (2013). A
graphical representation of the proposed approach for modeling
misclassification in zero-inflated Dirichlet-multinomial models, missZIDM, is
presented in Figure 1. The Markov chain Monte Carlo (MCMC) sampler used to
implement our model is outlined below in Algorithm 1 with more details
provided in the Supplementary Material.

$\bm{z}_{ijl}$$\bm{y}_{ijl}$$\tau_{ijl}$$\bm{\theta}_{tc}$$\alpha_{ijt}$$\bm{x}_{ij}$$\zeta_{ijt}$$\bm{\beta}_{\gamma_{t}}$$\bm{x}_{i}$$\bm{\beta}_{\eta_{t}}$$\bm{\beta}_{\psi_{t}}$$\bm{x}_{ijl}$$N$
$n_{i}$ $L_{ij}$ $T$ $C$ Figure 1: Graphical representation of missZIDM with
covariate dependence. Note that auxiliary parameters and hyperparameters have
been suppressed for clarity. $\bm{N}$ \- total observations; $\bm{n}_{i}$ \-
total measurements per observation; $\bm{L}_{ij}$ \- number of individuals at
each measurement; $\bm{T}$ \- true number of categories; $\bm{C}$ \- observed
number of categories.  Algorithm 1 MCMC Sampler

Input data $\bm{y}_{ijl}$, $\bm{x}_{ijl}$, $\bm{x}_{ij}$, $\bm{x}_{i}$.

Initialize parameters: $\bm{z}_{ijl}$, $\tau_{ijl}$, $\bm{\beta}_{\psi_{t}}$,
$\bm{\beta}_{\eta_{t}}$, $\bm{\beta}_{\gamma_{t}}$, $\omega_{\tau_{ijl}}$,
$\omega_{\zeta_{tij}}$, $\bm{\alpha}_{ij}$, $\bm{\zeta}_{ij}$, $\mu_{ij}$,
$\bm{a}_{t}$, $u_{t}$.

Specify hyperparameters:
$\mu_{\psi},\mu_{\eta},\mu_{\gamma},\sigma^{2}_{\psi},\sigma^{2}_{\eta},\sigma^{2}_{\gamma}$.

for iteration $m=1,\dots,M$ do

for $i=1,\dots,N$ do

for $j=1,\dots,n_{i}$ do

Update $\mu_{ij}\sim\mbox{Gamma}(L_{ij},\bar{\alpha}_{ij})$.

for $l=1,\dots,L_{ij}$ do

if $\bm{y}_{ijl}=\bm{z}_{ijl}$ then Update
$\tau_{ijl}\sim\text{Bernoulli}(\psi_{ijl}^{*})$.

else Update $\tau_{ijl}=1$.

end if

if $\tau_{ijl}$ = 0 then Update
$\bm{z}_{ijl}\sim\delta_{\bm{y}_{ijl}}(\bm{z}_{ijl})$.

else Update
$\bm{z}_{ijl}\sim\text{Multinomial}(1,\bm{\Theta}_{ij}\otimes\bm{\theta}_{c})$.

end if

Update
$\omega_{\tau_{ijl}}\sim\mbox{PG}(1,\bm{x}_{\psi_{ijl}}^{\prime}\bm{\beta}_{\psi_{t}})$.

end for

for $t=1,\dots,T$ do

Jointly update $\alpha_{ijt}$ and $\zeta_{ijt}$ with an Expand/Contract Step
via Koslovsky, (2023).

Update
$\alpha_{ijt}|\zeta_{ijt}=1\sim\mbox{Gamma}\left(\sum_{l=1}^{L_{ij}}I(z_{ijl}=t)+\gamma_{ijt},1+\mu_{ij}\right)$.

Update
$\omega_{\zeta_{tij}}\sim\mbox{PG}(1,\bm{x}_{\eta}^{\prime}\bm{\beta}_{\eta_{t}})$.

end for

end for

end for

for $t=1,\dots,T$ do

Update
$u_{t}\sim\mbox{Gamma}(\sum_{i=1}^{N}\sum_{j=1}^{n_{i}}\sum_{l=1}^{L_{ij}}I(z_{ijl}=t),\bar{a}_{t})$.

Update $\bm{\beta}_{\psi_{t}}\sim
N_{P_{\psi}}(\tilde{\bm{\beta}}_{\psi_{t}},\tilde{\Sigma}_{\psi_{t}})$.

Update $\bm{\beta}_{\eta_{t}}\sim
N_{P_{\eta}}(\tilde{\bm{\beta}}_{\eta_{t}},\tilde{\Sigma}_{\eta_{t}})$.

Update $\bm{\beta}_{\gamma_{t}}$ via a Metropolis-Hastings step.

for $c=1,\dots,C$ do

$a_{tc}\sim\mbox{Gamma}(\sum_{i=1}^{N}\sum_{j=1}^{n_{i}}\sum_{l=1}^{L_{ij}}I(y_{ijl}=c,z_{ijl}=t)+\nu_{tc},u_{t}+1)$.

end for

end for

end for

##  4 Empirical Studies

In this section, we first compare the proposed model in simulation to
alternative methods for handling false positives or false negatives in
multinomial data. In a second scenario, we compare the proposed model to a
multispecies occupancy model in settings designed to mimic the motivating
application.

The first scenario examines the estimation performance of missZIDM with
respect to the at-risk probability, $\bm{\eta}$, the probability of
misclassification, $\bm{\psi}$, the true relative abundance of each category,
$\bm{\Theta}$, and the confusion matrix, $\bm{\theta}^{*}$, at varying
percentages of at-risk observations and misclassification. Estimates for the
confusion matrix, $\bm{\theta}^{*}$, using the proposed model were obtained by
combining the estimated misclassification probabilities, $\bm{\psi}$, and
classification matrix, $\bm{\theta}$, accordingly. In the first scenario, we
compare missZIDM to a similar approach that does not accommodate false
negatives (missDM), an approach that assumes the true and observed
classifications follow Dirichlet-multinomial models which does not explicitly
model misclassification or handle false negatives (DMDM, similar to Swartz et
al., (2004)), and the recently developed zero-inflated Dirichlet-multinomial
model (ZIDM) (Koslovsky,, 2023), which is designed to accommodate false
negatives but not false positives. In this scenario, the models ignore any
potential covariates and therefore only estimate intercept terms for
$\bm{\eta}$, $\bm{\psi}$, and $\bm{\Theta}$. For the DMDM model, we set
$\nu_{tt}=(\mbox{logit}^{-1}(\mu_{\psi})/T+(1-\mbox{logit}^{-1}(\mu_{\psi})))\times
T/\mbox{logit}^{-1}(\mu_{\psi})$, which places a similar prior probability for
correct classification as the methods with misclassification indicators. All
methods were implemented in R using Rcpp (Eddelbuettel and François,, 2011).

In this scenario, we generated $N=50$ observations of $L_{ij}=100$ individuals
to cluster into $C=10$ categories. We assumed that the true number of
categories, $T$, matched the potentially observed number of categories, $C$.
We evaluated the model in four settings with varying percentages of at-risk
observations (1 - % true negatives) and misclassification (false positives).
In these settings, we set $n_{i}=1$ (i.e., no repeated measurements).
Observation-specific at-risk indicators were sampled from a Bernoulli
distribution with the probability of an at-risk observation set to either 0.25
or 0.75. The true classification of each individual was generated from a
Dirichlet-multinomial distribution with concentration parameters set to one
and overdispersion parameter set to 0.01, so that the model assumptions did
not match the true data generation process. Misclassification occurred with
0.25 or 0.75 probability. The observed classifications were generated from a
Dirichlet-multinomial model with a similar overdispersion parameter as above.
We set concentration parameters $\bm{\nu}_{t}$ equal to their index (e.g.,
$\nu_{tc}=c$ ) with the $t^{th}$ element also equal to one, placing the least
probability on a correct classification by chance. No validation data were
used in these settings.

In the second scenario, we investigate how the proposed method performs when
used for inference in multispecies occupancy-detection settings with data
generated to mimic the application data set. We evaluate the performance of
the method with varying percentages of overdispersion in the true counts. We
compare the proposed model to a similar version of the multispecies occupancy-
detection model presented in Wright et al., (2020) that assumes the ecological
processes follow a zero-inflated Poisson distribution and the observation
process follows a Dirichlet-multinomial model, which we refer to as DMZIP.
Code to implement DMZIP was adapted from Stratton, (2022).

Specifically, we generated $N=50$ observations with $n_{i}=5$ measurements per
observation and $C=10$ possible species to observe. No covariates were used in
the baseline simulation setting. The occurrence probability was set using the
observed proportion of zeros in the data, ranging from 23% to 90%. Instead of
assuming the total number of counts was fixed, as in scenario 1, the relative
activity or encounter rates for each species at each site visit was generated
from a negative binomial distribution with mean $\zeta_{it}*\lambda_{t}$ and
variance $(\zeta_{it}*\lambda_{t})^{2}/\sigma$, where $\zeta_{it}$ is the
site-level occupancy indicator for a given species, $\lambda_{t}$ is the
expected number of detections or encounter rate of species $t$ obtained from
the application data which ranged from 2.0 to 28.2, and $\sigma$ is an
overdispersion parameter. Note that as $\sigma$ increases, the variance of the
sampling distribution approaches the mean, and the data are more Poisson-like.
The models were evaluated with $\sigma\in\\{0.1,1,100\\}$ and 25% validation
data.

We then evaluated model performance on data generated similar to scenario 2
but with covariates informing the occupancy probability and the true encounter
rates. We simulated 5 continuous covariates from a standard normal
distribution in both levels of the model. In this setting, we set the
overdispersion parameter for the negative binomial distribution $\sigma=1$.
The intercept terms $\beta_{\eta_{t0}}$ ($\beta_{\gamma_{t0}}$) were randomly
sampled uniformly from $\mbox{logit}(0.25)$ to $\mbox{logit}(0.95)$ ($0$ to
$\log(10)$) with covariate effects set to $\pm 1$ ($\pm 0.2)$ with equal
probability. The off-diagonal elements of the classification matrix
$\bm{\theta}$ were sampled uniformly from 0.01 to 0.2 with diagonal elements
uniformly sampled from 0.5 to 0.95. Thereafter, the rows of $\bm{\theta}$ were
scaled to sum to one, and the individual classifications were sampled from a
Multinomial$(1,\bm{\theta}_{t})$. In this setting, we assumed 25% of the data
were validated. Additionally, we evaluated the models in various other data
generation settings including those with different sample sizes, sampling
efforts, and percent validated data.

In both scenarios, each of the MCMC algorithms were run for 5,000 iterations
treating the first 2,500 as burn-in and thinning to every other iteration,
providing 1,250 iterations for inference. We assumed non- or weakly-
informative priors
$\gamma_{tc}=\nu_{tc}=\sigma_{\eta}^{2}=\sigma_{\psi}^{2}=\sigma_{\gamma}^{2}=1$.
In settings with no covariates in the model, we set $\mu_{\eta}$ and
$\mu_{\psi}$ following the data generation process for all models. In settings
with covariates in the model, the prior mean for the regression coefficients
was set to 0. In the sensitivity analysis presented in the Supplementary
Material, we explore the impact of prior misspecification of these
hyperparameters on inference. To initialize each model, we set the true
classifications $\bm{Z}_{i}$ to the observed classifications $\bm{Y}_{i}$,
with $\bm{\tau}_{i}$ set accordingly. Auxiliary parameters,
$\omega_{\tau_{tijl}}$ and $\omega_{\zeta_{tij}}$, and at-risk indicators,
$\zeta_{tij}$ and $\alpha_{tc}$, were initialized at one. The auxiliary
parameters $u_{t}$ and $\mu_{ij}$ were randomly initialized from a Gamma(1,1).

We evaluated the models in terms of the average absolute value of the
difference between the estimated and true probabilities (ABS), Frobenius norm
(FROB), which is the square root of the sum of the squared difference of the
estimated and true probabilities, and 95% coverage probabilities (COV) for
$\bm{\eta}$, $\bm{\psi}$, $\bm{\theta}^{*}$, and $\bm{\Theta}$. Note that the
proposed method is the only method that provides estimates for all parameters
simultaneously. In settings where covariates were incorporated into the data
generation process (results presented in the Supplementary Material), the
models were compared with respect to the estimation of the occupancy
probability regression coefficients, $\bm{\beta}_{\eta}$, and the confusion
matrix, $\bm{\theta}^{*}$, because these maintained similar interpretation
among all models. Results we report below were obtained by averaging over 50
replicated data sets for each setting.

###  4.1 Results

In the first scenario, the estimation performance for the probability of an
at-risk observation, $\bm{\eta}$, improved as the true percentage of at-risk
observations increased for missZIDM and ZIDM, with missZIDM demonstrating
better performance in settings with more structural zeros. We observed that
the proposed missZIDM always outperformed missDM, which ignores potential
zero-inflation, when estimating the probability of misclassification,
$\bm{\psi}$. All methods obtained relatively similar estimation performance
for the true relative abundances, $\bm{\Theta}$, with the proposed method
demonstrating a slight advantage in the setting with 25% at-risk observations
and misclassification. Estimation accuracy for the confusion matrix,
$\bm{\theta}^{*}$, reduced as the percentage of misclassification increased
for all methods. The proposed method and DMDM both outperformed missDM with
respect to estimating the confusion matrix, $\bm{\theta}^{*}$, however, DMDM
obtained marginally better estimates when data were generated with higher
misclassification percentages. Recall that missZIDM, unlike DMDM, does not
directly provide estimates for $\bm{\theta}^{*}$, but they can be obtained
using the estimated $\bm{\psi}$ and $\bm{\theta}$ values. Both missZIDM and
ZIDM overestimated coverage probabilities for $\bm{\eta}$, while all methods
underestimated coverage probabilities for $\bm{\Theta}$. Coverage for
$\bm{\psi}$ and $\bm{\theta}^{*}$ varied as a function of misclassification
and at-risk observation probabilities, with the proposed model obtaining near
nominal coverage with lower classification probabilities. A similar trend was
observed for missDM when estimating $\bm{\psi}$.

Table 1:  Simulation Results for Scenario 1: Estimation performance for $N=50$ observations, $n_{i}=1$ measurements, $L_{ij}=100$ individuals, and $T=10$ categories at varying percentages of at-risk observations and misclassification with 0% validation data. ABS - absolute value of the difference between the estimated and true probabilities; FROB - Frobenius norm; COV - 95% coverage probabilities.  | 25% at-risk observations and 25% misclassification
---|---
| $\bm{\eta}$ |  | $\bm{\psi}$ |  | $\bm{\Theta}$ |  | $\bm{\theta}^{*}$
| ABS | FROB | COV |  | ABS | FROB | COV |  | ABS | FROB | COV |  | ABS | FROB | COV
missZIDM | 0.05 | 0.25 | 1.00 |  | 0.02 | 0.07 | 0.98 |  | 0.02 | 0.71 | 0.23 |  | 0.01 | 0.14 | 0.95
missDM | - | - | - |  | 0.18 | 0.56 | 0.00 |  | 0.05 | 1.49 | 0.18 |  | 0.04 | 0.62 | 0.34
DMDM | - | - | - |  | - | - | - |  | 0.05 | 1.66 | 0.16 |  | 0.01 | 0.15 | 1.00
ZIDM | 0.13 | 0.43 | 1.00 |  | - | - | - |  | 0.05 | 1.63 | 0.17 |  | - | - | -
| 25% at-risk observations and 75% misclassification
| $\bm{\eta}$ |  | $\bm{\psi}$ |  | $\bm{\Theta}$ |  | $\bm{\theta}^{*}$
| ABS | FROB | COV |  | ABS | FROB | COV |  | ABS | FROB | COV |  | ABS | FROB | COV
missZIDM | 0.09 | 0.33 | 1.00 |  | 0.36 | 1.15 | 0.00 |  | 0.11 | 3.76 | 0.09 |  | 0.08 | 1.29 | 0.66
missDM | - | - | - |  | 0.47 | 1.49 | 0.00 |  | 0.12 | 3.75 | 0.06 |  | 0.09 | 1.61 | 0.60
DMDM | - | - | - |  | - | - | - |  | 0.12 | 3.93 | 0.12 |  | 0.04 | 0.45 | 1.00
ZIDM | 0.14 | 0.45 | 1.00 |  | - | - | - |  | 0.12 | 3.93 | 0.03 |  | - | - | -
| 75% at-risk observations and 25% misclassification
| $\bm{\eta}$ |  | $\bm{\psi}$ |  | $\bm{\Theta}$ |  | $\bm{\theta}^{*}$
| ABS | FROB | COV |  | ABS | FROB | COV |  | ABS | FROB | COV |  | ABS | FROB | COV
missZIDM | 0.03 | 0.12 | 1.00 |  | 0.03 | 0.10 | 0.97 |  | 0.03 | 0.84 | 0.71 |  | 0.01 | 0.17 | 0.99
missDM | - | - | - |  | 0.09 | 0.28 | 0.10 |  | 0.03 | 0.82 | 0.71 |  | 0.02 | 0.35 | 0.89
DMDM | - | - | - |  | - | - | - |  | 0.03 | 0.81 | 0.70 |  | 0.01 | 0.15 | 1.00
ZIDM | 0.02 | 0.07 | 1.00 |  | - | - | - |  | 0.03 | 0.81 | 0.70 |  | - | - | -
| 75% at-risk observations and 75% misclassification
| $\bm{\eta}$ |  | $\bm{\psi}$ |  | $\bm{\Theta}$ |  | $\bm{\theta}^{*}$
| ABS | FROB | COV |  | ABS | FROB | COV |  | ABS | FROB | COV |  | ABS | FROB | COV
missZIDM | 0.03 | 0.11 | 1.00 |  | 0.40 | 1.26 | 0.00 |  | 0.07 | 1.99 | 0.51 |  | 0.08 | 1.40 | 0.64
missDM | - | - | - |  | 0.44 | 1.40 | 0.00 |  | 0.06 | 1.70 | 0.56 |  | 0.09 | 1.53 | 0.64
DMDM | - | - | - |  | - | - | - |  | 0.06 | 1.53 | 0.53 |  | 0.04 | 0.45 | 1.00
ZIDM | 0.01 | 0.04 | 1.00 |  | - | - | - |  | 0.06 | 1.53 | 0.53 |  | - | - | -

In scenario 2, we found the DMZIP model obtained the best estimation
performance for $\bm{\beta}_{\eta}$ and $\bm{\theta}^{*}$ when $\sigma=100$
(Table 2). However in settings with more overdispersion ($\sigma=1$ and
$0.1$), the proposed method outperformed DMZIP with respect to both
parameters. Notably, estimation performance was worse for both methods when
$\sigma=0.1$. Similar trends were observed with 75% validated data
(Supplementary Table S1). These results demonstrate how the proposed method is
preferred in the presence of overdispersion. Additionally, we found that
missZIDM obtained the best estimation performance for $\bm{\beta}_{\eta}$ and
$\bm{\theta}^{*}$ with covariates in the model, more species types, visits,
and sites (Supplementary Table S2).

Table 2:  Simulation Results for Scenario 2: Estimation performance for data generated similar to the application study with 25% validation. ABS - absolute value of the difference between the estimated and true probabilities; FROB - Frobenius norm; COV - 95% coverage probabilities.  | $\sigma=0.1$
---|---
| $\bm{\eta}$ |  | $\bm{\theta}^{*}$
| ABS | FROB | COV |  | ABS | FROB | COV
missZIDM | 0.27 | 0.94 | 0.07 |  | 0.01 | 0.25 | 0.81
DMZIP | 0.46 | 1.54 | 0.00 |  | 0.04 | 1.24 | 0.66
| $\sigma=1$
| $\bm{\eta}$ |  | $\bm{\theta}^{*}$
| ABS | FROB | COV |  | ABS | FROB | COV
missZIDM | 0.15 | 0.60 | 0.24 |  | 0.01 | 0.22 | 0.72
DMZIP | 0.22 | 0.76 | 0.08 |  | 0.03 | 1.04 | 0.75
| $\sigma=100$
| $\bm{\eta}$ |  | $\bm{\theta}^{*}$
| ABS | FROB | COV |  | ABS | FROB | COV
missZIDM | 0.17 | 0.66 | 0.10 |  | 0.02 | 0.28 | 0.73
DMZIP | 0.03 | 0.12 | 0.93 |  | 0.01 | 0.15 | 0.90

One of the major challenges of modeling measurement error in multinomial data
is non-identifiability of the parameters, because there is no information
contained in the raw data to inform zero-inflation or misclassification
probabilities. Our approach is designed to account for non-identifiability
through informative prior specifications and/or incorporating a subset of
validation data to inform parameter estimates. As such, inferential results
obtained by our method, and any method designed to model measurement error in
multinomial data, will be sensitive to the amount of validation data used to
inform the model as well as the specification of the hyperparameters. In the
Supplementary Material, we present an extensive sensitivity analysis of the
proposed model with varying percentages of validated data (Supplementary
Tables S3 and S4) and hyperparameter specification (Supplementary Tables S5
and S6). Based on these results, we recommend taking advantage of validation
data when available. With only 10% validation data, the proposed method was
able to obtain less than 0.05% bias for all probability estimates on average.
In the absence of validation data, the model performed better with lower
concentration parameters for the true relative abundances and the
classification probabilities. Additionally, we found that the model was
relatively robust to misspecification of the misclassification prior for all
other parameters. Similar results were observed for changes in the prior for
the at-risk probability.

##  5 Application Study

We demonstrate the proposed method on data collected in a multispecies bat
acoustic monitoring study conducted in British Columbia, Canada between 2016
and 2020. Details of the study design and data are found in Stratton et al.,
(2022) and Stratton, (2022). Briefly, one to six stationary acoustic recording
devices were placed in 55 sites following the North American Bat Monitoring
Program guidelines and were typically activated for seven nights (Loeb et
al.,, 2015). Similar to Stratton et al., (2022), we analyze detections from
the first and last nights to minimize potential overlap and dependencies.
There were 10 total bat species categories available for analysis, including
an other category for species that were difficult to detect acoustically or
that were not widespread. Each acoustic recording was classified using
Kaleidoscope Pro acoustic classification software for bats
(https://www.wildlifeacoustics.com). A unique attribute of these data is that
each acoustic recording was additionally validated by a bat expert. For this
analysis, we let half of all revisits from every site include validation data,
similar to Stratton et al., (2022). Additionally, we included site-specific
covariates measuring year (categorical with 2016 as the reference), annual
mean elevation (kilometers), precipitation (millimeters), and temperature
(degrees Celcius) for the occupancy (or at-risk) portion of the model. We
included nightly minimum air temperature (degrees Celcius), total
precipitation (millimeters), and percentage of the moon illuminated by the sun
(percent) measured from the centroid of the site at each visit to model the
encounter rates or true classifications depending on the model used. All
continuous covariates were standardized prior to analysis.

In this analysis, the average number of observed individuals per site visit
was 62.6, ranging from 1 to 992. The means, variances, and percent zero counts
for each species across site visits are presented in Table 3. We observed mean
counts ranging from 1.7 to 25.9 with variances ranging from 25.7 to 4339.6 and
percent zeros ranging from 0.26 to 0.83, indicating that the Poisson
assumption was not justified for the number of encounters.

Common Name | Scientific Name | Mean (Variance) | % Zero
---|---|---|---
Big brown bat |  Eptesicus fuscus (EPFU) | 3.5 (256.0) | 66
Hoary bat |  Lasiurus cinereus (LACI) | 3.1 (396.8) | 65
Silver-haired bat |  Lasionycteris noctivagans (LANO) | 12.4 (2703.4) | 34
California myotis |  Myotis californicus (MYCA) | 4.5 (264.3) | 54
Western small-footed myotis |  Myotis ciliolabrum (MYCI) | 2.9 (265.4) | 83
Western long-eared myotis |  Myotis evotis (MYEV) | 1.7 (28.2) | 63
Little brown myotis |  Myotis lucifugus (MYLU) | 25.9 (4339.6) | 26
Long-legged myotis |  Myotis volans (MYVO) | 2.7 (114.3) | 59
Yuma myotis |  Myotis yumanensis (MYYU) | 5.0 (625.5) | 69
Other | - | 0.98 (25.7) | 77
Table 3: Observed mean, variance, and % zero observations for each species in
the case study.

We compared the results of the proposed model and a Bayesian model similar in
construction to Wright et al., (2020) that assumes the true count process
follows a Poisson distribution, which we refer to as DMZIP. Both methods were
run for 10,000 iterations, treating the first 5,000 as burn-in and thinning to
every other iteration. The models were initialized similar to scenario 2 of
the simulation study. We assumed relatively weak or non-informative priors
with $\nu_{tc}=1$, $\gamma_{t}=1$, $\mu_{\beta_{\eta}}=\mu_{\beta_{\psi}}=0$,
and $\sigma_{\eta}^{2}=\sigma_{\psi}^{2}=\sigma_{\gamma}^{2}=1$. Similar
assumptions were made for DMZIP. Convergence and mixing of the models was
visually inspected using traceplots. Additionally, we ran another chain from a
different seed and compared the two chains with the Gelman-Rubin statistic,
which was less than 1.1 for each parameter (Brooks and Gelman,, 1998).

Figure 2 presents the estimated classification probabilities for the proposed
missZIDM model, as well as the differences with the DMZIP model. The models
found relatively similar results. The largest differences in misclassification
estimates were for EPFU, MYCO, and the other category, with DMZIP estimating
more misclassification. Additionally, the proposed method estimated more
misclassification for the LACI and MYYU species compared to DMZIP.

We investigated the estimated covariate associations for occupancy. Since
occupancy is a binary outcome, the estimates plotted in Figure 3 are
interpreted as log odds ratios for occupancy at each site visit. Overall, the
results were quite similar between the models. Neither method found much of an
effect for time on occupancy. However, both models did find a decrease in the
odds of occupancy for EPFU in 2019 compared to 2016, and DMZIP estimated an
increase in the odds of occupancy for LANO in 2020 compared to 2016. We found
that an increase in elevation was associated with an increase in the odds of
occupancy for most species, with the exception of MYVO. We observed mostly
negative associations between precipitation and occupancy, though most of the
95% credible intervals contained 0 for both models. The strongest relation was
for MYCI, where a millimeter increase in precipitation was associated with a
95% decrease in occupancy. EPFU, LACI, LANO, MYCA, MYCI, and MYYU were all
found to have positive associations between temperature and occupancy.
However, temperature was negatively associated with occupancy for MYVO.
Typically the proposed method was more conservative than DMZIP with respect to
parameter uncertainty for all covariate effects.

Figure 4 plots the estimated relations between temperature, precipitation, and
illumination and relative abundance using the proposed method. Exponentiation
of a regression coefficient is interpreted as the multiplicative factor of
change in the proportion of a compositional element with a one unit change in
the corresponding standardized covariate while holding all else constant (Chen
and Li,, 2013). We found that temperature typically had a positive association
with the true relative abundances for all species, with the exception of MYEV.
Again, we found similar negative relations between precipitation and the true
relative abundance. The association between illumination and relative
abundance varied across species. Supplementary Figure S1 presents the
associations between the same covariates and relative activity estimated with
the DMZIP model. It is important to note that the regression coefficients for
missZIDM and DMZIP at this level of the model do not have the same meaning as
the former assumes the true (latent) relative abundances follow a ZIDM
distribution while the latter assumes the latent encounter rate of each
species follows a Poisson distribution. Due to these differences in
interpretation, the results are not directly comparable. However, overall we
found very similar patterns in terms of the directionality of the relations,
with the proposed method typically providing more conservative estimates.

Figure 2: Posterior estimates of classification probabilities for each species
using the proposed missZIDM and DMZIP models.  Figure 3: Posterior estimates
of regression coefficients for occupancy for the proposed missZIDM and DMZIP
models. Dot represents the posterior mean with error bars capturing the 95%
credible intervals.  Figure 4: Posterior estimates of regression coefficients
for relative abundance for the proposed missZIDM model. Dot represents the
posterior mean with error bars capturing the 95% credible intervals.

##  6 Conclusions

In this work, we propose the first method for simultaneously accommodating
false positives and false negatives in multinomial data. Our model can
naturally incorporate existing knowledge of zero-inflation or
misclassification through prior specification and/or easily accommodate
validation data to inform the model. In simulation, we demonstrated that our
approach obtains similar or improved estimation performance for at-risk,
misclassification, true count, and misclassification probabilities compared to
alternative methods that ignore one or both forms of measurement error.
Motivated by data collected in multispecies occupancy-detection studies, we
further showed how the proposed approach can provide more accurate estimation
for occupancy and (mis)classification probabilities than existing approaches
when the observed count data are overdispersed. Existing multispecies
occupancy-detection models could be easily adjusted to accommodate potential
overdispersion by replacing the Poisson distribution for the ecological
process with negative binomial distributions. The implementation of this model
would require modeling the latent counts explicitly, similar to our approach
for modeling misclassification. This approach serves as a potentially viable
alternative for modeling misclassification in overdispersed multivariate count
data, but it would not be appropriate for multinomial or compositional data
settings where the total number of counts is fixed.

Motivated by the task of modeling imperfect classification in acoustic bat
data, the proposed method is applicable to other settings in which multinomial
data with potential misclassification are collected within and outside of
ecological research. Naturally, the method can be applied to other occupancy
monitoring studies in which potential misclassification occurs (e.g., citizen
science, aerial surveys). Another very promising application of the proposed
method in ecological research is to environmental DNA data (Ficetola et al.,,
2016; Schmidt et al.,, 2013; Willoughby et al.,, 2016; Lahoz-Monfort et al.,,
2016). In the context of conservation research and wildlife monitoring
studies, the proposed method can be extended in various directions to answer
pressing research questions. For example, to accommodate occupancy dynamics,
one could assume the probability a site is in a given occupancy state is
governed by a Markov process, similar to Miller et al., (2013). Additionally
in the current formulation, we assume that the total number of species, $T$,
matches the number of species observed. As such, our model does not provide
inference on species richness. To estimate species richness, we could take a
parameter-expanded data augmentation approach, similar to (Royle and Dorazio,,
2012), under the assumption that there were other species in the surveyed
regions that went unobserved, (i.e., $T>>C$).

Supplementary Material The Supplementary Material contains detailed
derivations of the MCMC algorithm, sensitivity analysis, and additional
figures. missZIDM R package This file contains code to generate data similar
to the simulation study, apply the method, and perform inference.

[Acknowledgments] MDK gratefully acknowledges the support of NSF grant
DMS-2245492. The opinions, findings, and conclusions expressed are those of
the authors and do not necessarily reflect the views of the NSF. The authors
would like to thank Dr. Kathi Irvine for providing useful discussions about
modeling bat acoustic data.

## References

  * Aitchison and Ho, (1989) Aitchison, J. and Ho, C. (1989).  The multivariate Poisson-log normal distribution.  Biometrika, 76(4):643–653.
  * Bayley and Peterson, (2001) Bayley, P. B. and Peterson, J. T. (2001).  An approach to estimate probability of presence and richness of fish species.  Transactions of the American Fisheries Society, 130(4):620–633.
  * Blasco-Moreno et al., (2019) Blasco-Moreno, A., Pérez-Casany, M., Puig, P., Morante, M., and Castells, E. (2019).  What does a zero mean? Understanding false, random and structural zeros in ecology.  Methods in Ecology and Evolution, 10(7):949–959.
  * Broms et al., (2015) Broms, K. M., Hooten, M. B., and Fitzpatrick, R. M. (2015).  Accounting for imperfect detection in hill numbers for biodiversity studies.  Methods in Ecology and Evolution, 6(1):99–108.
  * Brooks and Gelman, (1998) Brooks, S. P. and Gelman, A. (1998).  General methods for monitoring convergence of iterative simulations.  Journal of Computational and Graphical Statistics, 7(4):434–455.
  * (6) Chambert, T., Grant, E. H. C., Miller, D. A., Nichols, J. D., Mulder, K. P., and Brand, A. B. (2018a).  Two-species occupancy modelling accounting for species misidentification and non-detection.  Methods in Ecology and Evolution, 9(6):1468–1477.
  * Chambert et al., (2015) Chambert, T., Miller, D. A., and Nichols, J. D. (2015).  Modeling false positive detections in species occurrence data under different study designs.  Ecology, 96(2):332–339.
  * (8) Chambert, T., Waddle, J. H., Miller, D. A., Walls, S. C., and Nichols, J. D. (2018b).  A new framework for analysing automated acoustic species detection data: Occupancy estimation and optimization of recordings post-processing.  Methods in Ecology and Evolution, 9(3):560–570.
  * 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):418–442.
  * Chiquet et al., (2021) Chiquet, J., Mariadassou, M., and Robin, S. (2021).  The Poisson-lognormal model as a versatile framework for the joint analysis of species abundances.  Frontiers in Ecology and Evolution, 9:188.
  * Copas and Hilton, (1990) Copas, J. and Hilton, F. (1990).  Record linkage: Statistical models for matching computer records.  Journal of the Royal Statistical Society: Series A (Statistics in Society), 153(3):287–312.
  * Datta et al., (2021) Datta, A., Fiksel, J., Amouzou, A., and Zeger, S. L. (2021).  Regularized bayesian transfer learning for population-level etiological distributions.  Biostatistics, 22(4):836–857.
  * Devarajan et al., (2020) Devarajan, K., Morelli, T. L., and Tenan, S. (2020).  Multi-species occupancy models: Review, roadmap, and recommendations.  Ecography, 43(11):1612–1624.
  * Dorazio et al., (2011) Dorazio, R. M., Gotelli, N. J., and Ellison, A. M. (2011).  Modern methods of estimating biodiversity from presence-absence surveys.  Biodiversity Loss in a Changing Planet, pages 277–302.
  * Dorazio et al., (2006) Dorazio, R. M., Royle, J. A., Söderström, B., and Glimskär, A. (2006).  Estimating species richness and accumulation by modeling species occurrence and detectability.  Ecology, 87(4):842–854.
  * Eddelbuettel and François, (2011) Eddelbuettel, D. and François, R. (2011).  Rcpp: Seamless R and C++ integration.  Journal of Statistical Software, 40:1–18.
  * Ficetola et al., (2016) Ficetola, G., Taberlet, P., and Coissac, E. (2016).  How to limit false positives in environmental DNA and metabarcoding?  Molecular Ecology Resources, 16(3):604–607.
  * Frénay and Verleysen, (2013) Frénay, B. and Verleysen, M. (2013).  Classification in the presence of label noise: A survey.  IEEE Transactions on Neural Networks and Learning Systems, 25(5):845–869.
  * Guillera-Arroita et al., (2017) Guillera-Arroita, G., Lahoz-Monfort, J. J., van Rooyen, A. R., Weeks, A. R., and Tingley, R. (2017).  Dealing with false-positive and false-negative errors about species occurrence at multiple levels.  Methods in Ecology and Evolution, 8(9):1081–1091.
  * Iknayan et al., (2014) Iknayan, K. J., Tingley, M. W., Furnas, B. J., and Beissinger, S. R. (2014).  Detecting diversity: Emerging methods to estimate species diversity.  Trends in Ecology & Evolution, 29(2):97–106.
  * Jasra et al., (2005) Jasra, A., Holmes, C., and Stephens, D. (2005).  Markov chain Monte Carlo Methods and the label switching problem in Bayesian mixture modeling.  Statistical Science, 20(1):50–67.
  * Jiang et al., (2021) Jiang, S., Xiao, G., Koh, A. Y., Kim, J., Li, Q., and Zhan, X. (2021).  A Bayesian zero-inflated negative binomial regression model for the integrative analysis of microbiome data.  Biostatistics, 22(3):522–540.
  * Kellner and Swihart, (2014) Kellner, K. F. and Swihart, R. K. (2014).  Accounting for imperfect detection in ecology: A quantitative review.  PloS One, 9(10):e111436.
  * Koslovsky, (2023) Koslovsky, M. D. (2023).  A Bayesian zero-inflated Dirichlet-multinomial regression model for multivariate compositional count data.  Biometrics.
  * Lahoz-Monfort et al., (2016) Lahoz-Monfort, J. J., Guillera-Arroita, G., and Tingley, R. (2016).  Statistical approaches to account for false-positive errors in environmental DNA samples.  Molecular Ecology Resources, 16(3):673–685.
  * Lele et al., (2012) Lele, S. R., Moreno, M., and Bayne, E. (2012).  Dealing with detection error in site occupancy surveys: What can we do with a single survey?  Journal of Plant Ecology, 5(1):22–31.
  * Lindén and Mäntyniemi, (2011) Lindén, A. and Mäntyniemi, S. (2011).  Using the negative binomial distribution to model overdispersion in ecological count data.  Ecology, 92(7):1414–1421.
  * Loeb et al., (2015) Loeb, S., Rodhouse, T., Ellison, L., Lausen, C., Reichard, J., Irvine, K., Ingersoll, T., Coleman, J., Thogmartin, W., Sauer, J., et al. (2015).  A plan for the North American Bat Monitoring Program (NABat).  General Technical Report-Southern Research Station, USDA Forest Service.
  * MacKenzie et al., (2003) MacKenzie, D. I., Nichols, J. D., Hines, J. E., Knutson, M. G., and Franklin, A. B. (2003).  Estimating site occupancy, colonization, and local extinction when a species is detected imperfectly.  Ecology, 84(8):2200–2207.
  * MacKenzie et al., (2002) MacKenzie, D. I., Nichols, J. D., Lachman, G. B., Droege, S., Andrew Royle, J., and Langtimm, C. A. (2002).  Estimating site occupancy rates when detection probabilities are less than one.  Ecology, 83(8):2248–2255.
  * Miller et al., (2013) Miller, D. A., Nichols, J. D., Gude, J. A., Rich, L. N., Podruzny, K. M., Hines, J. E., and Mitchell, M. S. (2013).  Determining occurrence dynamics when false positives occur: Estimating the range dynamics of wolves from public survey data.  PLoS One, 8(6):e65808.
  * Miller et al., (2011) Miller, D. A., Nichols, J. D., McClintock, B. T., Grant, E. H. C., Bailey, L. L., and Weir, L. A. (2011).  Improving occupancy estimation when two types of observational error occur: Non-detection and species misidentification.  Ecology, 92(7):1422–1428.
  * Molinari, (2008) Molinari, F. (2008).  Partial identification of probability distributions with misclassified data.  Journal of Econometrics, 144(1):81–117.
  * Mulick et al., (2022) Mulick, A. R., Oza, S., Prieto-Merino, D., Villavicencio, F., Cousens, S., and Perin, J. (2022).  A Bayesian hierarchical model with integrated covariate selection and misclassification matrices to estimate neonatal and child causes of death.  Journal of the Royal Statistical Society Series A: Statistics in Society, 185(4):2097–2120.
  * Neelon, (2019) Neelon, B. (2019).  Bayesian zero-inflated negative binomial regression based on Pólya-Gamma mixtures.  Bayesian Analysis, 14(3):829.
  * Pérez et al., (2007) Pérez, C. J., Girón, F. J., Martín, J., Ruiz, M., and Rojano, C. (2007).  Misclassified multinomial data: A Bayesian approach.  RACSAM, 101(1):71–80.
  * Polson et al., (2013) Polson, N. G., Scott, J. G., and Windle, J. (2013).  Bayesian inference for logistic models using Pólya–Gamma latent variables.  Journal of the American Statistical Association, 108(504):1339–1349.
  * Royle and Dorazio, (2012) Royle, J. A. and Dorazio, R. M. (2012).  Parameter-expanded data augmentation for Bayesian analysis of capture–recapture models.  Journal of Ornithology, 152(2):521–537.
  * Royle and Link, (2006) Royle, J. A. and Link, W. A. (2006).  Generalized site occupancy models allowing for false positive and false negative errors.  Ecology, 87(4):835–841.
  * Royle and Nichols, (2003) Royle, J. A. and Nichols, J. D. (2003).  Estimating abundance from repeated presence–absence data or point counts.  Ecology, 84(3):777–790.
  * Ruiz-Gutierrez et al., (2016) Ruiz-Gutierrez, V., Hooten, M. B., and Campbell Grant, E. H. (2016).  Uncertainty in biological monitoring: A framework for data collection and analysis to account for multiple sources of sampling bias.  Methods in Ecology and Evolution, 7(8):900–909.
  * Scharf et al., (2022) Scharf, H. R., Lu, X., Williams, P. J., and Hooten, M. B. (2022).  Constructing flexible, identifiable and interpretable statistical models for binary data.  International Statistical Review, 90(2):328–345.
  * Schmidt et al., (2013) Schmidt, B. R., Kéry, M., Ursenbacher, S., Hyman, O. J., and Collins, J. P. (2013).  Site occupancy models in the analysis of environmental dna presence/absence surveys: A case study of an emerging amphibian pathogen.  Methods in Ecology and Evolution, 4(7):646–653.
  * Shuler et al., (2021) Shuler, K., Verbanic, S., Chen, I. A., and Lee, J. (2021).  A Bayesian nonparametric analysis for zero-inflated multivariate count data with application to microbiome study.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 70(4):961–979.
  * Spiers et al., (2022) Spiers, A. I., Royle, J. A., Torrens, C. L., and Joseph, M. B. (2022).  Estimating species misclassification with occupancy dynamics and encounter rates: A semi-supervised, individual-level approach.  Methods in Ecology and Evolution, 13(7):1528–1539.
  * Steorts et al., (2016) Steorts, R. C., Hall, R., and Fienberg, S. E. (2016).  A Bayesian approach to graphical record linkage and deduplication.  Journal of the American Statistical Association, 111(516):1660–1672.
  * Stratton, (2022) Stratton, C. (2022).  Strattonch/CoupledUncoupled: Coupling validation effort manuscript release (v1.0.0).  Zenodo. https://doi.org/10.5281/zenodo.6040068.
  * Stratton et al., (2022) Stratton, C., Irvine, K. M., Banner, K. M., Wright, W. J., Lausen, C., and Rae, J. (2022).  Coupling validation effort with in situ bioacoustic data improves estimating relative activity and occupancy for multiple species with cross-species misclassifications.  Methods in Ecology and Evolution, 13(6):1288–1303.
  * Swartz et al., (2004) Swartz, T. B., Haitovsky, Y., Vexler, A., and Yang, T. Y. (2004).  Bayesian identifiability and misclassification in multinomial data.  Canadian Journal of Statistics, 32(3):285–302.
  * Tancredi and Liseo, (2011) Tancredi, A. and Liseo, B. (2011).  A hierarchical Bayesian approach to record linkage and population size problems.  The Annals of Applied Statistics, 5(2B):1553–1585.
  * Tyre et al., (2003) Tyre, A. J., Tenhumberg, B., Field, S. A., Niejalke, D., Parris, K., and Possingham, H. P. (2003).  Improving precision and reducing bias in biological surveys: Estimating false-negative error rates.  Ecological Applications, 13(6):1790–1801.
  * Wadsworth et al., (2017) Wadsworth, W. D., Argiento, R., Guindani, M., Galloway-Pena, J., Shelburne, S. A., and Vannucci, M. (2017).  An integrative Bayesian Dirichlet-multinomial regression model for the analysis of taxonomic abundances in microbiome data.  BMC Bioinformatics, 18(1):94.
  * Wang et al., (2020) Wang, S., Wang, L., and Swartz, T. B. (2020).  Inference for misclassified multinomial data with covariates.  Canadian Journal of Statistics, 48(4):655–669.
  * Willoughby et al., (2016) Willoughby, J. R., Wijayawardena, B. K., Sundaram, M., Swihart, R. K., and DeWoody, J. A. (2016).  The importance of including imperfect detection models in eDNA experimental design.  Molecular Ecology Resources, 4(16):837–844.
  * Wright et al., (2020) Wright, W. J., Irvine, K. M., Almberg, E. S., and Litt, A. R. (2020).  Modelling misclassification in multi-species acoustic data when estimating occupancy and relative activity.  Methods in Ecology and Evolution, 11(1):71–81.
  * Xu et al., (2015) Xu, L., Paterson, A. D., Turpin, W., and Xu, W. (2015).  Assessment and selection of competing models for zero-inflated microbiome data.  PloS One, 10(7):e0129606.
  * Zhang and Yi, (2020) Zhang, X. and Yi, N. (2020).  NBZIMM: Negative binomial and zero-inflated mixed models, with application to microbiome/metagenomics data analysis.  BMC Bioinformatics, 21(1):1–19.