# Modelling heterogeneity in the classification process in multi-species
distribution models can improve predictive performance.

Kwaku Peprah Adjei1,200 Robert B. O’Hara1,20 Wouter Koch2,40 and Anders
Finstad 2,30

1Department of Mathematical Sciences, Norwegian University of Science and
Technology, Trondheim, Norway
2Center for Biodiversity Dynamics, Norwegian University of Science and
Technology, Norway
3Department of Natural History , Norwegian University of Science and
Technology, Norway
4Norwegian Biodiversity Information Center , Trondheim , Norway Corresponding
author.E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>wouter.koch<EMAIL_ADDRESS>
###### Abstract

  1. 1.

Species distribution models and maps from large-scale biodiversity data are
necessary for conservation management. One current issue is that biodiversity
data are prone to taxonomic misclassifications. Methods to account for these
misclassifications in multispecies distribution models have assumed that the
classification probabilities are constant throughout the study. In reality,
classification probabilities are likely to vary with several covariates.
Failure to account for such heterogeneity can lead to bias in parameter
estimates.

  2. 2.

Here we present a general multispecies distribution model that accounts for
heterogeneity in the classification process. The proposed model assumes a
multinomial generalised linear model for the classification confusion matrix.
We compare the performance of the heterogeneous classification model to that
of the homogeneous classification model by assessing how well they estimate
the parameters in the model and their predictive performance on hold-out
samples. We applied the model to gull data from Norway, Denmark and Finland,
obtained from Global Biodiversity Information Facility.

  3. 3.

Our simulation study showed that accounting for heterogeneity in the
classification process increased precision by 30 % and reduced accuracy and
recall by 6%. Applying the model framework to the gull dataset did not improve
the predictive performance between the homogeneous and heterogeneous models
due to the smaller misclassified sample sizes. However, when machine learning
predictive scores are used as weights to inform the species distribution
models about the classification process, the precision increases by 70%.

  4. 4.

We recommend multiple multinomial regression to be used to model the variation
in the classification process when the data contains relatively larger
misclassified samples. Machine prediction scores should be used when the data
contains relatively smaller misclassified samples.

Key words: Bayesian models, Citizen Science, False positives, Machine
learning, Misclassification, Multispecies distribution models

##  1 Introduction

Species distribution models are essential ecology and conservation management
tools that predict how natural and human factors affect biodiversity (Elith
and Leathwick, 2009; Vermeiren et al., 2020). With an increasing amount of
biodiversity data from multi-species surveys available to scientists, multi-
species distribution models (hereafter mSDMs) have become widely used in
analysing these data. These mSDMs model data at the community level, identify
the important variables that drive species co-occurrences and predict the
distribution of species in a community (Ovaskainen and Soininen, 2011; Pollock
et al., 2014; Hui et al., 2015).

However, the biodiversity data obtained from these surveys can be subject to
observation errors. These observation errors include false positives
(recording absent species as present and misidentifying one species for
another) and false negatives (failing to record present species). These false
positives are usually a species that has been misclassified. It would be a
great advantage if the observation was correctly classified, rather than being
discarded once it has been identified as a false positive. The false positives
may arise from imperfect classifiers (Wright et al., 2020; Spiers et al.,
2022), observer error, and many other sources. In this study, we collectively
describe all the false positives and misidentifications as misclassification.
We use the term true states in this study to describe the correct or actual
observation identity we are interested in modelling. False negatives are
mostly accounted for in occupancy models by modelling the imperfect detection
in the observation model. Failure to account for or correct these errors leads
to biases in inferences about state variables such as occupancy probabilities,
covariate effects and relative activity (Ferguson et al., 2015; Wright et al.,
2020; Clare et al., 2021), leading to an impairment in decision making
(Hoekman, 2021).

The methods to deal with misclassification from biodiversity data can be
grouped into data review methods and model-based methods (Clare et al., 2021).
Data review methods require complete and proper data collection and processing
methods. This process can be very demanding as it is challenging to control
for misclassification. This makes the model-based methods more popular when
working with large-scale datasets from large-scale biodiversity data vendors
like the Global Biodiversity Information Facility (GBIF hereafter; GBIF.Org,
2022). Model-based methods estimate classification probabilities jointly with
the true state variables of interest. Model-based methods attempting to
account for misclassification in multispecies occupancy models currently
include modelling misclassification with detection heterogeneity (Ferguson et
al., 2015; Louvrier et al., 2018; Clement et al., 2022), integrating multiple
observers records with other methods such as distance sampling and N-Mixture
models (Hoekman, 2021), supervised methods with extra information from
observation confirmation or verification (Ferguson et al., 2015; Guillera-
Arroita et al., 2017), site confirmation (Clare et al., 2021) and other
calibrated methods. These methods need extra data from the verification
process, which helps in estimating the misclassification probabilities in a
semi-supervised setting (Spiers et al., 2022) and makes the parameters in the
model identifiable (Guillera-Arroita et al., 2017). The above-mentioned
studies have either used verified data collected on the site-level (where the
occupancy state of a species is known at a site and not at the individual
sample level; Chambert et al., 2018b), on aggregated individual sample level
using a multinomial model with site-covariates (Wright et al., 2020), or on
individual sample-level validation data which helps in modelling non-species
identities (morphospecies) to species identities (Spiers et al., 2022). It is
also worth stating that some studies have explored misclassification in
abundance (Conn et al., 2013) and capture-recapture (Augustine et al., 2020)
models.

Furthermore, these previous studies assumed that the misclassification
probabilities are homogeneous (constant) across the study. In reality, the
classification probabilities may vary with environmental covariates (such as
field conditions; Conn et al., 2013) or observer experience (especially when
ascertaining how well each observer classifies a report in citizen science
projects will be informative; Arazy and Malkinson, 2021; Johnston et al.,
2022), distance from a transect when using transect data (Conn et al., 2013),
picture quality, etc. An attempt at modelling the heterogeneity in the
classification process would be to model this effect as part of the true state
model (such as abundance or occupancy). However, this approach may not solve
the problem of heterogeneity in the classification process since the estimates
of the true state observation process parameters only serve as informed priors
to the classification process (Spiers et al., 2022).

A more correct approach to model this heterogeneity is adding the covariate
effect to the classification process. Some studies on dynamic false positive
single-species occupancy models have modelled temporal changes in false
positives using year as a covariate (Sutherland et al., 2013; Miller et al.,
2013; Kéry and Royle, 2020), showing the possibility to model
misclassification trends over time. Our study attempts to model variation in
classification probabilities in mSDMs by modelling the probability of
classifying a true state identity with a multinomial generalised linear model.
To our knowledge, no previous work has been done on this. Failure to account
for the heterogeneity in the misclassification probabilities will lead to
biased estimates in the parameters of the true state observation process (such
as species abundance, richness, and occupancy probabilities) and reduce the
model’s predictive performance.

Moreover, recent efforts to correctly classify observations from biodiversity
surveys have relied on machine learning (hereafter ML) algorithms (Suzuki-Ohno
et al., 2022; Saoud et al., 2020; Lotfian et al., 2021; Willi et al., 2019;
Keshavan et al., 2019; Borowiec et al., 2022; Koch et al., 2022). These ML
algorithms use sounds and/or images of observations to predict the true
identity of the individual observations, and they can be trained to mimic
expert verification of observations (Keshavan et al., 2019; Ponti and Seredko,
2022; Langenkämper et al., 2019). These ML algorithms use a prediction score
(a value that shows the weight of predicting the observations as something
else) to predict the possible list of the true identities of the individual
reported observations. These prediction scores and a list of possible true
identities provide information about the classification process of each
observation. They can be used to model heterogeneity in the classification
process. This study is the first to model the heterogeneity in the
classification process by using the prediction scores to weigh the
distribution of the reported observations and predict the distribution of the
actual observation identities.

Here, we present a joint model that simultaneously models the true state
variables of interest (such as occupancy and abundance) and the heterogeneity
in the classification process. Our model set-up extends the work done by
Wright et al. (2020) and Spiers et al. (2022) by a) allowing the
classification probabilities to vary with a covariate, b) allowing additional
morphostates designations in imperfect classifications (for example, having
"unidentified woodpecker" as part of the classified observations), c) using ML
prediction scores as weights to account for heterogeneity in the
classification process and d) performing variable selection on the
classification process covariate to check for potential mSDM overfit. We
compared the classification performance of our model with models that assume a
homogeneous classification probability done by Wright et al. (2020) and
(Spiers et al., 2022) through simulation studies. We parameterise our model
with citizen science data on gulls in Norway, Finland and Denmark from
iNaturalist (Matheson, 2014) downloaded from the GBIF (GBIF.Org, 2022).

##  2 Methodology

###  2.1 Model framework

The proposed framework assumes a model for where the true state’s identity is
"present" (we use "state(s)" in this work to refer to taxon identity as well
as any other identification category or morpho-states, such as an unidentified
group). This true state can be on any taxonomic level. However, the true
states are often unknown and treated as latent. An individual is sampled and
classified as another state, known as the "reported" or "classified" state,
with a probability. These reported states can also be on any taxonomic level
or any unidentified group (such as in studies that include morphospecies;
Spiers et al., 2022). We further assume that the individuals are verifiable
(and that we have information on the verification process) and that the
verified state approximates the true state identity.

####  2.1.1 True state observation process model

We specify a relative abundance model (specifically a multinomial logit model)
for the true state observation process. Our objective is to show how to model
heterogeneity in the classification process and not to make inferences about
the true state’s abundance, so we chose a model that was easier to fit and
understand to explain the distribution of the true state.

Let $\lambda_{is}$ be the average number of $s=1,2,\ldots,S$ true individuals
at location $i=1,2,\ldots,R$, which describes the abundance of the true states
over space. This means that intensity can be modelled as an inhomogeneous
point process (PPM) which assumes that the data are dependent on the
environment covariate, or as a log-Cox Gaussian Point process (LGCP) where we
assume a spatial dependency in the observed data (Renner et al., 2015). Here,
the mean intensity is modelled using the inhomogeneous point process and
defined as:

$\ln(\lambda_{is})=\beta_{0s}+\sum_{j=1}^{n_{e}}x_{ij}\beta_{js},$ (1)

where $\mathbf{\beta_{0s}}$ is the intercept of state $s$, $\beta_{js}$ is the
effect of covariate $j$ on the intensity of true state $s$, $x_{ij}$ is the
$j^{th}$ covariate that affects the observation of the state at location $i$
and $n_{e}$ is the number of covariates in the observation process model. Note
that we assume there are no species interactions in our relative abundance
model, and this could have been added as a random effect in the true state
intensity definition (equation (1)).

Let $p_{is}$ be the probability of true state $s$ being at location $i$. We
estimate this probability from the mean intensities as follows:

$p_{is}=\frac{\lambda_{is}}{\sum_{s}\lambda_{is}},$ (2)

where $\lambda_{is}$ is given by (1).

The true state $s=1,2,\ldots,S$ at location $i=1,2,\ldots,R$ is then a
realisation from a categorical distribution with probability $p_{is}$.

This model specification for the ecological process used here is similar to
the occupancy dynamics and encounter rate model used by Spiers et al. (2022)
by assuming that occupancy probability equals one; and similar to the model
used by Wright et al. (2020) by assuming Poisson counts with intensity
$\lambda_{is}$ (refer to Table 3 for the link between our model framework and
that of Spiers et al. (2022) and Wright et al. (2020)).

####  2.1.2 Classification process model

For each observation, we assume there was imperfect classification and that
the true state identities can be classified into $k=1,2,\ldots,K$ states.
Assuming the true state distribution represents the true point pattern for the
species under study, each reported information can be seen as a draw from any
of the $K$ states under consideration with a given probability. As mentioned
above, these states could be on any taxonomic level or include any
unidentified group. For example, one could have four true states: common gull,
herring gull, Audouin’s gull and Sooty gull. These species can be reported in
three states: large while-headed gulls, large black-headed gulls and other
gulls. An example of the classification probability (confusion matrix) is
shown in Table 1.

|  | Reported states |
---|---|---|---
True states | Large white-headed gulls | Large black-headed gulls | Other gulls
Common gull | 0.8 | 0.1 | 0.1
Herring gull | 0.9 | 0 | 0.1
Audouin’s gull | 0 | 0.9 | 0.1
Sooty gull | 0 | 1 | 0
Table 1: Example of confusion matrix that applies to our model.

Let $\Omega_{sk}$ be the probability that a single true state
$s\in\\{1,\ldots,S\\}$ is classified as state $k\in\\{1,\ldots,K\\}$. The
probabilities across all the possible $k$ states sum to $1$. In studies with
homogeneous classification probabilities, the confusion matrix for the
classification can be expressed as:

$\mathbf{\Omega}=\begin{bmatrix}\mathbf{\Omega}_{1}\\\ \mathbf{\Omega}_{2}\\\
\vdots\\\ \mathbf{\Omega}_{S}\\\
\end{bmatrix}=\begin{bmatrix}\Omega_{11}&\Omega_{11}&\cdots&\Omega_{1k}&\cdots&\Omega_{1K}\\\
\Omega_{21}&\Omega_{22}&\cdots&\Omega_{2k}&\cdots&\Omega_{2K}\\\
\vdots&\vdots&\cdots&\vdots&\cdots&\vdots\\\
\Omega_{S1}&\Omega_{S2}&\cdots&\Omega_{Sk}&\cdots&\Omega_{SK}\\\
\end{bmatrix}$ (3)

where the rows correspond to the true states $s$ and the columns correspond to
the reported states $k$.

We model the heterogeneity in the classification probabilities by fitting a
multinomial generalised linear model (Fahrmeir et al., 2013) to each of the
rows of $\mathbf{\Omega}$ defined in equation (3). We refer to this approach
as the multiple multinomial generalised linear model (MMGLM, hereafter). We
define the linear predictor of the MMGLM as:

$\begin{split}\zeta_{isk}&=\omega_{0sk}+\sum_{j=1}^{n_{c}}z_{ij}\omega_{jsk}\\\
&=\begin{bmatrix}{\omega_{0}}_{11}&{\omega_{0}}_{12}&\cdots&{\omega_{0}}_{1K}\\\
{\omega_{0}}_{21}&{\omega_{0}}_{22}&\cdots&{\omega_{0}}_{2K}\\\
\vdots&\vdots&\cdots&\vdots\\\
{\omega_{0}}_{S1}&{\omega_{0}}_{S2}&\cdots&{\omega_{0}}_{SK}\\\
\end{bmatrix}+\ldots+\begin{bmatrix}{\omega_{n}}_{11}&{\omega_{n}}_{12}&\cdots&{\omega_{n}}_{1K}\\\
{\omega_{n}}_{21}&{\omega_{n}}_{22}&\cdots&{\omega_{n}}_{2K}\\\
\vdots&\vdots&\cdots&\vdots\\\
{\omega_{n}}_{S1}&{\omega_{n}}_{S2}&\cdots&{\omega_{n}}_{SK}\\\
\end{bmatrix}*z_{in_{c}}\\\ \text{and in matrix notation:}\\\
\boldsymbol{\zeta}_{i}&=\boldsymbol{\omega_{0}}+\sum_{j=1}^{n_{c}}\boldsymbol{\omega_{j}}z_{ij};\end{split}$
(4)

where $\boldsymbol{\omega_{0}}$ is an $S\times K$ matrix of intercepts and
$\boldsymbol{\omega_{j}}$ is an $S\times K$ matrix of covariate $j$ effects,
and $z_{ij}$ is the $j^{th}$ covariate that drives the classification process
at location $i$. Using equation (4) as the definition for the linear
predictor, our estimates of the parameters $\omega_{0sk}$ and $\omega_{jsk}$
are identifiable with reference to one reported state. Therefore for each true
state identity $s$, the classification probabilities for each reported state
$k=1,2,\ldots,K-1$ with reference to state $K$ is modelled as:

$\ln\bigg{(}\frac{\Omega_{sk}}{\Omega_{sK}}\bigg{)}=({\omega_{0}}_{sk}-{\omega_{0}}_{sK})+({\omega_{1}}_{sk}-{\omega_{1}}_{sK})*\textbf{z}_{1}+\ldots+({\omega_{n}}_{sk}-{\omega_{n}}_{sK})*\textbf{z}_{n};$
(5)

where $\textbf{z}_{1}$, $\ldots$, $\textbf{z}_{n}$ are covariate vectors.

This general framework has $S\times(K-1)\times(n_{c}+1)$ parameters to be
estimated, where $S$ is the number of true states, $K$ is the number of
reported states, and $n_{c}$ is the number of covariates in the classification
process. Estimating these parameters can be very computationally expensive as
the number of true states, reported states and covariates increases and would
require significant numbers of misclassified state samples to estimate them.
Therefore, we explored simplified forms of the generalised model in equation
(5).

A simplified case of equation (5) assumes that the classification covariate z
only affects the probability of correctly classifying the states. In this
instance, the matrix $\boldsymbol{\omega_{j}}$ has its off-principal diagonal
elements $\omega_{jsk}=0$ for $s\neq k$ for covariate $j$, and these
parameters are not estimated (This is our study scenario "fixed-covariate" in
Table 2). This simplification reduces the number of parameters estimated for
the classification process by $n_{c}\times(K-S-1)$, where $S$ is the number of
true states, $K$ is the number of reported states and $n_{c}$ is the number of
classification covariates. A further simplification would also assume that the
average probability of correctly classifying the true state does not change
with species, that is, $\omega_{0sk}$ is the same for all $s=k$ (This is our
study scenario "fixed-intercov" in Table 2). The classification process
heterogeneity is captured by the covariate effect $\omega_{1sk}$ for all
$s=k$. The latter further reduces the number of parameters estimated by $S-1$.
This last simplification is useful especially when the states are very
similar, and one would expect that their average classification probabilities
would be the same.

Given which true state ($s$) was present at location $i$, the reported state
at that location is a draw from $K$ states with probability $\mathbf{\Omega}$.

$Y_{i}|V_{i}=s\sim\text{Categorical}(\Omega_{is\cdot}),$ (6)

where $\Omega_{is.}$ refers to the probability of being classified as any of
the $K$ states given that the true state identity at location $i$ is state
$s$.

####  2.1.3 Variable selection

Fitting a complex model with these many parameters can result in an overfitted
model. An overfitted model captures the pattern and noise in the training
data, but performs poorly on validation or test data (Montesinos López et al.,
2022). In some cases, the true state observation and classification process
covariates may be highly correlated. These correlated covariates can inflate
standard errors (reduce the precision) of the estimated parameters (Yu et al.,
2015; Roberts et al., 2017; Caradima et al., 2019). To avoid overfitting the
model, there is a need to perform variable selection. In this study, we
performed Bayesian variable selection, specifically the spike and slab prior
to the classification process covariates (for review of Bayesian variable
selection see O’Hara and Sillanpää, 2009). For each of the classification
process covariates, we re-define the linear predictor in equation (4) as:

$\begin{split}\zeta_{isk}&=\omega_{0sk}+\sum_{j=1}^{n_{c}}\psi_{j}z_{ij}\omega_{jsk};\\\
\text{with}\quad\psi_{j}&\sim Bernoulli(q_{j})\end{split}$ (7)

where $\psi_{j}$ is the indicator that variable $j$ is selected with the
expected probability $q_{j}$. The variable selection indicator specified in
equation (7) jointly selects the variable that affects all the true states in
the model but can also be made state-specific (Ovaskainen and Abrego, 2020).
Probabilities $q_{j}$ closer to $1$ indicate that the variable contributes
much to the model and should be selected, and those closer to $0$ indicate
that the variable contributes less and can be discarded.

[sansbold]  Classification type Model True state observation process
Classification process Heterogeneous Variable / covariate
$\ln(\lambda_{is})=\beta_{0s}+x_{1i}\beta_{1s}+x_{2i}\beta_{2s}$
$\ln(\frac{\zeta_{isk}}{\zeta_{isK}})=\omega_{0sk}+\psi_{i}z_{1i}\omega_{1sk}$
fixed-covariate
$\ln(\lambda_{is})=\beta_{0s}+x_{1i}\beta_{1s}+x_{2i}\beta_{2s}$
$\ln(\frac{\zeta_{isk}}{\zeta_{isK}})=\omega_{0sk}+\psi_{i}z_{1i}\omega_{1ss}$,
where $\omega_{1sk}=0$ for $s\neq k$ fixed-intercov
$\ln(\lambda_{is})=\beta_{0s}+x_{1i}\beta_{1s}+x_{2i}\beta_{2s}$
$\ln(\frac{\zeta_{isk}}{\zeta_{isK}})=\omega_{0sk}+\psi_{i}z_{1i}\omega_{1ss}$,
where $\omega_{1sk}=0$ for $s\neq k$ and $\omega_{0sk}$ is the same for $s=k$
Homogeneous intercept
$\ln(\lambda_{is})=\beta_{0s}+x_{1i}\beta_{1s}+x_{2i}\beta_{2s}$
$\ln(\frac{\zeta_{isk}}{\zeta_{isK}})=\omega_{0sk}$ constant
$\ln(\lambda_{is})=\beta_{0s}+x_{1i}\beta_{1s}+x_{2i}\beta_{2s}$
$\Omega_{sk}\sim Dirichlet(\alpha_{sk})$, where $\alpha_{sk}\sim Exp(1)$ main
$\ln(\lambda_{is})=\beta_{0s}+x_{1i}\beta_{1s}+x_{2i}\beta_{2s}+\psi_{i}z_{1i}\beta_{3s}$
$\Omega_{sk}\sim Dirichlet(\alpha_{sk})$, where $\alpha_{sk}\sim Exp(1)$

Table 2: Variations in the MMGLM for the classification process defined by
equation (7) and observation process defined by equation (1) used as our study
scenarios. There are three heterogeneous models: covariate, fixed-covariate
and fixed-intercov, and three homogeneous models: intercept, constant, and
main.

###  2.2 Modelling heterogeneity using ML prediction scores

Some studies give weight to the true identities of the reported observations
or state, for example, because they use machine learning to classify the
observation. These prediction scores (such as $F_{1}$ score, mean square
error, and logarithmic loss) are not classification probabilities but are
values that indicate how well the ML algorithm classifies data in the test
sample. We can use this information to model the heterogeneity in the
classification process and predict the true state identity as follows:

$\begin{split}V_{i}|Y_{i}&\sim\text{Categorical}(p_{iks});\\\ \text{with}\quad
p_{iks}&=\frac{\lambda_{ik}w_{iks}}{\sum_{k}\lambda_{ik}w_{iks}}\end{split}$
(8)

where $\lambda_{ik}$ is the average abundance of reported state $k$ at site
$i$ and $w_{iks}$ is the predictive score of the $k^{th}$ reported state to
true state $s$ at sampling unit $i$.

It must be noted that this approach is a non-parametric approach to account
for heterogeneity in the classification whereas the MMGLM is a parametric
approach. Moreover, modelling the heterogeneity in the classification process
using ML prediction scores models the covariate effects on the expected
abundance of the reported states and corrects them using the prediction scores
as weights to obtain the relative abundance of the true states. However, the
MMGLM models the covariate effect on the expected abundance of true states and
estimates the heterogeneity in the classification process using a parametric
model. The prediction of the true state identity is done by weighing the
expected intensity of true states with the estimated classification covariate.

###  2.3 Generalisation of model framework

The proposed framework generalises the existing mSDMs that account for
misclassification in occupancy models. Wright et al. (2020) provided a
framework to account for the homogeneity in the classification process, and
our model is connected to this by using the relationship between the
multinomial and the Poisson distribution (Steel et al., 1953) for the
observation process as well as using a species by species constant model for
the classification process (Table 3). Wright et al. (2020) further provided
arguments that their proposed models were generalised forms of models for the
binary detection of two species (Chambert et al., 2018a), single species with
count detections (Chambert et al., 2015) and single species with binary
detections (Chambert et al., 2018b). Since our proposed framework can be seen
as a heterogeneous version of Wright et al. (2020), our framework is also a
generalisable form of the models in Chambert et al. (2015, 2018b, 2018a).
Spiers et al. (2022) provided an individual-level semi-supervised approach
that estimates species misclassification with occupancy dynamics and encounter
rates, and our model is connected to this if we assume a homogeneous
classification process (that is $\boldsymbol{\omega_{1}}=\textbf{0}$), and
assume that the occupancy probability is $1$ (Table 3).

[sansbold]  Author Model Framework Link to our model Wright et al. (2020) and
the models their proposed framework generalises such as Chambert et al. (2015,
2018b, 2018a) True state observation process: absolute counts
$z_{is}\sim\text{Bernoulli}(\psi_{is})$
$[V_{ijs}|z_{is}=1]\sim\text{Poisson}(\lambda_{ijs})$ True state observation
process: relative abundance Choose $\psi_{is}=1$, then $z_{is}=1$ for all
locations $i$ and true states $s$.
$[V_{ijs}|z_{is}=1]\sim\text{Poisson}(\lambda_{ijs})$ using the relationship
between multinomial and Poisson distribution (Steel et al., 1953).
Classification process:
$[Y_{ijs}|V_{ijs}=v_{ijs},z_{is}=1]\sim\text{Multinomial}(v_{ijs},\Omega_{sk})$,
where $\Omega_{sk}\sim\text{Dirichlet}(\alpha)$ Classification process:
$[Y_{ijs}|V_{ijs}=1,z_{is}=1]\sim\text{Categorical}(\Omega_{sk})$, where
$\Omega_{sk}$ can be chosen as any of the homogeneous models described in
Table 2 Spiers et al. (2022) True state observation process: occupancy
dynamics and encounter rates $z_{ist}\sim\text{Bernoulli}(\psi_{ist})$
$[V_{ijst}]\sim
Categorical\bigg{(}\frac{\lambda_{ijst}z_{ist}}{\sum_{s}\lambda_{ijst}z_{ist}}\bigg{)}$;
True state observation process: relative abundance Choose $t=1,\psi_{ist}=1$,
then $z_{ist}=1$ for all locations $i$ and true states $s$. $[V_{ijs}]\sim
Categorical\bigg{(}\frac{\lambda_{ijs}}{\sum_{s}\lambda_{ijs}}\bigg{)}$;
Classification process:
$[Y_{ijs}|V_{ijs}]\sim\text{Categorical}(\Omega_{sk})$, where
$\Omega_{sk}\sim\text{Dirichlet}(\alpha)$ Classification process:
$[Y_{ijs}|V_{ijs}=1]\sim\text{Categorical}(\Omega_{sk})$, where $\Omega_{sk}$
can be chosen as any of the homogeneous models described in Table 2.

Table 3: Extensions of our proposed models from homogeneous classification
process studies done by Wright et al. (2020) and Spiers et al. (2022). The
table specifies the true state observation process model for Wright et al.
(absolute abundance model), Spiers et al. (occupancy dynamics and encounter
rate model) and ours (relative abundance model); and also the classification
process model for Wright et al. and Spiers et al. (homogeneous classification
process with classification probabilities simulated from Dirichlet
distribution) and ours from heterogeneous models described in Table 2. Since
our framework extends the work done by (Wright et al., 2020), it is safe to
say that they are also generalised forms of Chambert et al. (2015, 2018b,
2018a).

###  2.4 Simulation Study

To demonstrate how our proposed model works and its use in prediction, we
performed a simulation study for $S=2$ true states, and $K=3$ reported states
over $R=1000$ sites. We simulated two covariates for the true state
observation process and one for the classification, all from a Normal
distribution with a mean of $0$ and variance of $1$. The true intensity was
simulated from equation (1). The intercepts of the model for the two true
states were chosen as $\beta_{01}=-1$, $\beta_{02}=0$ and covariate effect for
each true state was chosen as $\beta_{11}=4$, $\beta_{12}=-2$, $\beta_{21}=0$
and $\beta_{22}=0$ (that is, state $2$ is used as a reference category). The
intercept and covariate effect for the classification process was chosen as
follows:

$\begin{split}\omega_{0}&=\begin{bmatrix}2&0.5&0\\\ 1&1&0\\\
\end{bmatrix};\omega_{1}=\begin{bmatrix}3&-1&0\\\ -1&1&0\\\
\end{bmatrix}.\end{split}$

These values were chosen to obtain significant sample sizes of misclassified
states. We simulated $200$ datasets with a heterogeneous classification
process using the variable model in Table 2 (referred to as the "full model"),
$200$ with a homogeneous classification process by assuming
$\boldsymbol{\omega_{1}}=\textbf{0}$ (a matrix of zeros) in equation (5)
(referred to as "reduced model") and another $200$ with the covariate effect
for the classification process (using the variable model in Table 2) that is
correlated to the true state observation process covariate (referred to as
"correlation model"). The first simulated dataset explored modelling
heterogeneity’s effect on the classification process, whereas the latter
explored the effect of having correlated covariates for the classification and
true state observation process. The second assessed the effect of overfitting
the classification process model (adding heterogeneity to the classification
process when it should be homogeneous). Moreover, we assessed the impact of
the number of misclassified samples on the mSDMs predictive performance. We
increased the principal diagonal components of $\boldsymbol{\omega_{0}}$ by
$6$ to obtain a reduction in the number of misclassified samples simulated.

We randomly withheld two hundred of the true state identities for each dataset
simulated as our validation sample. The number of validation samples was not
varied since Spiers et al. (2022) found that the number of validation samples
had a modest effect on the model’s predictive ability. We fitted the model
under the various scenarios described in Table 2 to the data and evaluated the
model’s predictive performance on the validation sample.

###  2.5 Case study: Gulls dataset

The proposed model was used to analyse gull dataset downloaded from GBIF
(GBIF.Org, 2022). The database hosts over two billion occurrence observations
with well over a million observers (website visited on 17th February 2023). We
were interested in the iNaturalist records since they have community
verifications (Matheson, 2014). The observers collected these occurrence
records and uploaded their observations with images and/or sounds that allowed
for verification. The reported observations go through iNaturalist community
verification and are accepted as research grade when two-thirds of the
community agreed to the taxon identification (Ueda, 2020), at which point they
are published on GBIF. We assumed that the community-accepted taxon name is
the true state V. We checked the iNaturalist platform to track the
identification process of the observations and use the first reported
identification as the reported state Y.

We obtained observations for some species of gulls in Denmark, Finland and
Norway from 2015 to 2022. Specifically, we selected great black-backed gulls
(Larus marinus), herring (Larus argentatus), common gulls (Larus canus) and
lesser black-backed gulls (Larus fuscus) because the iNaturalist website
reported that these species are commonly misclassified as the other. Any other
species reported apart from the above-mentioned species were labelled as
"others". We used annual precipitation (accessed from the raster package;
Hijmans et al., 2015) as true state observation process covariates in the
model since it has been noted in some literature to affect the distribution of
sea birds such as gulls (Jongbloed, 2016; Algimantas and Rasa, 2010).

The data obtained were presence-only records. Exploratory analysis revealed
that there were no multiple observations at the same location for our selected
species. We, therefore, assumed that our locations were discrete and treated
the data as a marked point process, where the species that was reported at a
location is given a value of $1$ and $0$ for the other species in this study.
If we had a species list, then we could have modelled it as a repeated marked
point process at the same location and treated the sites as a random effect.

Citizen science data are known to be affected by several sources of bias. Some
common biases are spatial bias (Johnston et al., 2022; Tang et al., 2021), and
misclassifications (Johnston et al., 2022; Tulloch et al., 2013). We only
accounted for the misclassifications in this study, as we are interested in
explaining the classification process and not making inferences about the
abundance of the gulls. Citizen scientists have been reportedly known to
correctly classify species as they gain experience reporting the species
(Vohland et al., 2021). We, therefore, modelled the variation in the
classification process by using the number of reports made by each observer as
a covariate in the classification process. If an observer has more than ten
observations, the extra number of observations was calibrated at $10$. We used
this number of observations for an observer as a measure of experience
(Johnston et al., 2018; Kelling et al., 2015), although there are other
indices for measuring effort or experience of the citizen scientist (Santos-
Fernandez and Mengersen, 2021; Vohland et al., 2021).

We also used ML algorithm to obtain prediction scores (specifically $F_{1}$
scores) for the possible true state identity of our downloaded data. The ML
algorithm was a Convolutional Neural Network (a modified form in Koch et al.,
2022) trained with data from all citizen science observations of any species
in Norway. Since the ML algorithm is trained with all birds data from GBIF in
Norway, we trained all our six study scenarios models summarised in Table 2
with all data for our selected gull species in Norway and all data reported
before 2022 in Finland and Denmark and used all data reported in 2022 in
Finland and Denmark as our validation sample. The summary of the
classifications (true and false positives) in the training and validation
sample is presented in Supplementary information 2 (S2-2, S2-3).

###  2.6 Fitting and evaluating the model

We ran all the analyses using the Bayesian framework using the Markov chain
Monte Carlo approach in the NIMBLE package (de Valpine et al., 2017) using the
R software (R Core Team, 2022). We chose the priors for all true state
observation parameters from a normal distribution with a mean of $0$ and
standard deviation of $10$, and we chose the priors for the classification
process parameters from Normal distribution with a mean of $0$ and standard
deviation of $1$. For the scenarios: constant and main model, we chose the
priors of the confusion matrix ($\Omega$) from the Dirichlet distribution with
parameter alpha ($\alpha$), which has a prior of an exponential distribution
with mean $1$.

We ran 3 chains, each with 10000 iterations; the first 5000 iterations were
chosen as the burn-in. We check the convergence of the models by visually
inspecting the trace plots and ignoring models with a Gelman-Rubin statistic
(Brooks and Gelman, 1998) value greater than $1.1$. We kept a fifth of the
remaining samples in each of the chains.

We used accuracy (the proportion of predicted true states identities from all
the predictions of the validation samples), precision (the proportion of
mismatched true states in the validation sample that were correctly classified
from the predictions) and recall (the proportion of correct true states in the
validation sample retrieved from the predictions) as performance metrics. We
used a Bayesian approach and got the posterior distributions for the
parameters. Higher values of the validation metrics indicated the preferred
model. We also checked how well the model estimated the true state observation
process parameters ($\boldsymbol{\beta_{0}}$, $\boldsymbol{\beta_{1}}$ and
$\boldsymbol{\beta_{2}}$) and the classification process parameters
($\boldsymbol{\omega_{0}}$ and $\boldsymbol{\omega_{1}}$) by estimating the
bias (difference between the true value and the estimated value) and precision
of the parameters.

##  3 Results

###  3.1 Simulation study

####  3.1.1 Predictive performance

We illustrated the gain in model performance by using the accuracy, recall and
precision of our model’s predictions (Figure 1 A and B). When data was
simulated from the full and correlated models, there was a strong indication
that the predictive performance of mSDMs improved when the variability in the
classification process was included. That is the "variable" model performed
best for the full and correlation models with the highest accuracy, recall and
precision values (Figure 1 A (i - ii), B (i - ii)). The simplified
heterogeneous models (fixed-intercov and fixed-covariate) however did not
perform any better than the homogeneous models (Figure 1 A (i - ii), B (i -
ii)). This suggested that simplifying the heterogeneous classification model
did not yield any improvement in predictive performance, and the heterogeneous
model that captures the entire variability (in this case variable model; Table
2) would be the best predictive model. When the classification process
covariate was modelled as part of the observation process (main model), the
model’s predictive performance also performed similarly to the homogeneous
models (Figure 1 A and B).

Figure 1: Boxplot of validation metrics (accuracy, precision, recall) from the
six study models defined in Table 2 on the two hundred (200) withheld samples
out of the thousand (1000) samples simulated in each dataset. Accuracy is the
proportion of withheld samples that were correctly classified, recall is the
proportion of correctly classified samples that were retrieved from the
withheld samples, and precision is the proportion of the misclassified samples
that were correctly classified. Each boxplot shows the median and the
interquartile range (25 - 75% quartiles). Each column shows the type of model
used to simulate the dataset: "full" refers to using the variable/covariate
model in Table 2, "reduced" refers to using the intercept model in Table 2 and
"correlation" refers to using the variable model in Table 2, but with
correlated true state observation and classification process covariates. The
rows correspond to changes made to the number of misclassified samples in the
simulated dataset: "Baseline" refers to using the values defined in section
2.4 and "Decrease" refers to reducing the number of misclassified samples by
diagonal elements of $\mathbf{\omega}$ by 6 as described in section 2.4

Overfitting a homogeneous classification process with a heterogenous one did
not have any effect on the mSDM’s predictive performance (Figure 1 A iii) and
B iii)). We expected the overfitted heterogeneous models to have poor
predictive performance (Montesinos López et al., 2022), but the heterogeneous
and homogeneous performed similarly (with equal recall, accuracy and precision
across all six study models). The Bayesian variable selection probability
indicated that the homogeneous classification model was better (with the
probability of including classification covariates in heterogeneous models
being $0.359\pm 0.012$; Supplementary information 2; S2-1). Although the
simplified heterogeneous models did not yield improvement in predictive
performance, they performed similarly to the variable model in the variable
selection process.

####  3.1.2 Effect of number of misclassified samples

As we increased the number of misclassified samples in our simulated data, the
precision increased by on average $30\%$ and accuracy and recall reduced by
$6\%$ (Figure 1 A (i - ii), B (i - ii)). This decrease in accuracy and recall
could be attributed to the reduction in the number of correct classifications
in the simulated data as the number of misclassified samples increased.
Moreover, the classification parameters were estimated better when the number
of misclassified samples was higher, leading to the high precision of
predictions (Supplementary information 2; S2-1). This suggested that our
proposed model will be beneficial when one has a substantial number of
misclassified samples.

####  3.1.3 Bias in observation and classification process parameters

Although failure to account for misclassification in mSDMs can result in
biased true state observation process parameters (Wright et al., 2020; Spiers
et al., 2022), any method used to account for misclassification in mSDMs has a
small effect on the accuracy and precision of the true state observation
process parameters. The bias of the observation process parameters was
consistently low for all six models, and the coverage was higher for all the
scenarios under the full and reduced model (Supplementary information 2;
S2-1). All the scenarios studied accounted for misclassification of some sort,
thereby correcting for the bias in the observation parameters estimates
(Wright et al., 2020; Spiers et al., 2022). The classification process
parameters were estimated more accurately and precisely for the variable model
than the other models (Supplementary information; S2-1). This was only
possible in the case where we had enough misclassified samples. This suggests
that if the objective of a study is to predict true state identity with mSDMs,
then modelling the full heterogeneity can improve predictive performance; if
the aim is inference on true state distribution, then heterogeneous models may
not provide any advantage over homogeneous models.

[sansbold]  Method Accuracy Precision Recall Variable selection probability
Variable / Covariate 0.97 0.1 0.99 0.71 Constant 0.97 0.1 0.99 - Intercept
0.97 0.1 0.99 - Main 0.97 0.1 0.99 0.29 Fixed intercov 0.97 0.1 0.99 0.70
Fixed covariate 0.97 0.1 0.99 0.71 Machine Learning 0.89 0.8 0.90 -

Table 4: Validation metrics of the models under study on the withheld gull
dataset. The accuracy is the proportion of correctly classified validated
data, the precision is the proportion of mismatched identities that were
correctly matched and recall is the proportion of correctly matched identities
that were recovered. The number of validated samples was 384 out of which 10
were mismatched.

###  3.2 Case study: Gull dataset

All six study scenario models performed equally well regarding their
predictive performance with high accuracy and recall but smaller precision
(Table 4). The poor precision value could not be attributed to the
insignificance of the classification covariate (observer experience) in
explaining the heterogeneity in the classification process since the variable
selection probabilities are closer to $1$ (Table 4) but to the small
misclassification sample sizes (Supplementary Information 1; S2-2, S2-3).
However, the precision increased from 10% to 80% when the heterogeneity in the
classification process was accounted for by using the prediction scores from
the Machine learning algorithm (Table 4). The ML algorithm’s prediction scores
were individual observation-specific which provided direct information to the
observation process model. However, the six classification models depended on
the misclassified sample size to capture the heterogeneity in the
classification process. This suggests that one remedy to improve mSDM’s
predictive performance for data with very small misclassified samples is to
use ML weights to account for heterogeneity in mSDMs.

Figure 2: Summary of results from the model fit to gull dataset showing A)
Probability of correct classification for the common (Larus canus), herring
(Larus argentatus), great black-backed (Larus marinus) and lesser black-backed
gulls (Larus fuscus) and B) the distribution of the experiences of the
observers used in the modelling. The ribbon around the correct classification
probability estimates represents the 95 % credible interval of the estimates.

Although the study scenario models had smaller precision, it was observed that
the probability of correctly classifying the gull species in Denmark, Finland
and Norway increased with the experience of the observer (Figure 2). The
pattern showed that observers have a high chance of making mistakes on their
first few reports, and they get better as the number of reports increased
(Vohland et al., 2021).

##  4 Discussion

The main objective of this paper was to propose a general framework to account
for misclassifications from imperfect classifications (such as those from
surveys) and uncertain classifications (from automated classifiers) in mSDMs.
This work builds on previous work by Spiers et al. (2022); Wright et al.
(2020) by accounting for the heterogeneity in the classification probabilities
while allowing the classified categories to be more than the verified species
(such as unknown species, morphospecies etc.). Moreover, we assessed the
effect of overfitting a homogeneous classification process on the predictive
performance of mSDMs and provided ways of checking the overfitting of the
classification process model.

Our study bridges the knowledge gap in the literature on accounting for
misclassification in mSDMs by modelling the heterogeneity in the
classification process. False positives are inevitable in biodiversity data
(Miller et al., 2013; Kéry and Royle, 2020). To model these false positives in
this study, we presented the true state observation process as one model and
the classification process as another model in a hierarchical form. We model
the classification probabilities for each true state identity as a multinomial
generalised linear model. This specification generalises the modelling of the
classification process to model effects of covariates as fixed or random
effects or both. For example, one can estimate the classification
probabilities of each observer in volunteer-collected data by assuming a
random observer effect. This formulation for the classification process also
mitigates the modelling problems of using the Dirichlet distribution as the
prior for the classification probabilities (Spiers et al., 2022).

Furthermore, the specification of a separate state-space model for the true
state observation process in the proposed framework allows the use of various
multi-species observation models (such as joint species distribution models
(Tobler et al., 2019; Ovaskainen and Abrego, 2020), Royle-Nichols model for
abundance (Royle and Nichols, 2003), amongst many others) to model the
distribution of the true states. Our simulation study showed that the proposed
model framework could estimate the true state relative abundance parameters
(with the bias of estimated parameters close to zero; Supplementary
information 2; S2-1), an observation noted in previous studies that use
observation confirmation design to model misclassification in mSDMs (Spiers et
al., 2022; Wright et al., 2020; Kéry and Royle, 2020). We have shown that the
true state observation process model is a simplified form of occupancy and
abundance models that account for misclassification (Table 3), so we believe
our proposed framework can be extended to any design used to collect and
verify data on the true states (for example, point processes, distance
sampling, site-confirmation (multi-method) design, etc.), and any model used
to fit the data (for example, multi-state occupancy model (Kéry and Royle,
2020), joint species distribution models (Tobler et al., 2019; Ovaskainen and
Abrego, 2020)).

Modelling the true state observation process with more complex models than the
relative abundance models used in the study would add another level of
hierarchical structure to the true state observation process (for example
modelling detection probability or true occupancy state). This complexity
could introduce confounding of the classification process and observation
process parameters and with frequentist estimation approaches make the
likelihood multi-modal (Kéry and Royle, 2020). This study did not explore such
issues and further work can be done on this. A possible solution in the
Bayesian framework to avoid such confounding issues would be to model the
different processes with separate covariates, choose a good prior for the mSDM
parameters, and use repeated survey visit data to model the true state
observation process (Kéry and Royle, 2020). Moreover, the identifiability or
confounding issues could also be tackled by using data with much information
on detection and false positive detections such as those derived from acoustic
surveys (Clement et al., 2022) and integrating occupancy data that are not
susceptible to misclassifications such as those from camera traps to those
with misclassifications (Kéry and Royle, 2020).

Accounting for the heterogeneity in the classification process increases the
predictive performance of mSDMs. The homogeneous classification models may
sometimes be unable to explain the variation in the classification process
(Conn et al., 2013), leading to poor model predictive performance as a result
of overfitting (Montesinos López et al., 2022). The simulation study showed a
30% increase in precision and a 6% decrease in accuracy and recall when the
heterogeneity in the classification process was accounted for in the mSDMs
(Figure 1 A, B). However, there was no change in predictive performance when a
heterogeneous classification model overfitted a homogeneous classification
process (Figure 1 A, B) due to the small classification covariate effect size,
observed from the bias of parameter estimates and low Bayesian variable
selection probability (Figure 1 A, B). Since the predicted posterior
probability for the true state’s identity heavily relies on the weights from
the misclassification probability (Supplementary Information 1), failure to
account for that would mean our posterior probability would be incorrectly
estimated. The incorrectly predicted probability would lead to the
underestimation of the prediction of the ranges of coverage and possibly
abundance in the true states (Molinari-Jobin et al., 2012).

Fitting a more complex true state observation model with the covariate that
explains the heterogeneity of the classification process does not provide
enough information to improve the mSDM’s predictive performance. Previous
studies have shown that the estimates of the observation process model inform
the estimation of the classification probabilities (Spiers et al., 2022), but
the variability in the classification process cannot be inferred from
variability in the observation process model (Figure 1 A (i - iii),
Supplementary Information 1). Ecologists should therefore model the
variability in the classification in its process model to gain the advantage
in the mSDMs predictive performance.

Our model was parameterised with volunteer-collected gull data. These
volunteer-collected data have several sources of bias in their generation,
such as spatial bias, and misidentification of species, among many others. We
acknowledge that all these sources of biases may be present in the data, but
we only modelled the misidentification of species by using the number of
previously collected data as a proxy measure for the observer’s experience in
the classification process model. The predictive performance of the
homogeneous and heterogeneous models was approximately the same due to small
misclassified samples (19 misclassified out of 1382 samples in training data
(Supplementary Information S2-2) and 10 misclassified out of 378 samples in
validation data (Supplementary Information S2-3)). However, the estimated
covariate effect gives an idea of how the experience affects the probability
of classifying a new observation. Specifically, the probability of correctly
identifying the correct species increases with the observer’s experience, as
is noted in some literature (Santos-Fernandez and Mengersen, 2021; Vohland et
al., 2021; Johnston et al., 2018; Kelling et al., 2015). Therefore, there is a
trade-off between the model’s ability to correctly classify mismatched data
(precision) and understanding the covariate’s effect driving the
classification process when there are relatively small misclassified samples.

The inclusion of ML prediction scores in the mSDMs to account for the
heterogeneity in the classification process increased the precision of our
predictions by 70% (Table 4). These ML prediction scores are observation
specific and provide much information about the classification process to
increase the precision of the model. The information from the ML does not
depend on the misclassified sample sizes but on the quality of the images
(Koch et al., 2022), making them advantageous to use in accounting for
heterogeneity in the classification process when misclassified sample sizes
are small (like we have in our gull data).

The proposed model framework in this study is flexible and can be generalised
into any species distribution model and integrated distribution model. The
framework proposed fits into the frameworks provided by Wright et al. (2020)
and Spiers et al. (2022) and any framework their study generalises. Our
proposed classification process model, MMGLM, improved the predictive
performance of mSDMs, but it heavily relies on the misclassified samples size.
Furthermore, the confusion matrix defined in the model framework allows for
the classification of different taxonomic groups, as opposed to just the
species by species confusion matrix in Wright et al. (2020) and including
morphospecies in the classification categories (Spiers et al., 2022). This
will make it possible for citizen science data analysts to account for the
misclassification of data at any level in the data collection process. We
recommend that variable or model selection is performed during the analysis to
check for overfitting. Moreover, ecologists should explore using ML prediction
scores (where the prediction scores are available) as weights in mSDMs that
aim at predicting true state distributions especially when the data has a
small misclassified sample size.

##  5 Data Availability

The Gulls data used for this paper can be downloaded from GBIF
(https://doi.org/10.15468/dl.h24bp5). All codes used for this paper are
available at data dyrad with https://doi.org/10.5061/dryad.0rxwdbs51.

##  6 Conflict of Interest

The authors declare no conflict of interest.

##  7 Author Contribution

All authors were involved in the idea conception and writing of the
manuscript, KPA and RBO were involved with the methodology development, WK ran
the ML algorithm and KPA ran the entire analysis used in this paper and led
the writing of the manuscript.

## References

  * Algimantas and Rasa (2010) Algimantas, P. and Rasa, U. (2010) Abundance and seasonal migration of gulls(laridae) on the lithuanian baltic sea coast.  Acta Biologica Universitatis Daugavpiliensis, 10, 147–164.
  * Arazy and Malkinson (2021) Arazy, O. and Malkinson, D. (2021) A framework of observer-based biases in citizen science biodiversity monitoring: Semi-structuring unstructured biodiversity monitoring protocols.  Frontiers in Ecology and Evolution, 435.
  * Augustine et al. (2020) Augustine, B. C., Royle, J. A., Linden, D. W. and Fuller, A. K. (2020) Spatial proximity moderates genotype uncertainty in genetic tagging studies.  Proceedings of the National Academy of Sciences, 117, 17903–17912.
  * Borowiec et al. (2022) Borowiec, M. L., Dikow, R. B., Frandsen, P. B., McKeeken, A., Valentini, G. and White, A. E. (2022) Deep learning as a tool for ecology and evolution.  Methods in Ecology and Evolution, 13, 1640–1660.
  * 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, 434–455.
  * Caradima et al. (2019) Caradima, B., Schuwirth, N. and Reichert, P. (2019) From individual to joint species distribution models: A comparison of model complexity and predictive performance.  Journal of Biogeography, 46, 2260–2274.
  * Chambert et al. (2018a) 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, 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, 332–339.
  * Chambert et al. (2018b) 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, 560–570.
  * Clare et al. (2021) Clare, J. D., Townsend, P. A. and Zuckerberg, B. (2021) Generalized model-based solutions to false-positive error in species detection/nondetection data.  Ecology, 102, e03241.
  * Clement et al. (2022) Clement, M. J., Royle, J. A. and Mixan, R. J. (2022) Estimating occupancy from autonomous recording unit data in the presence of misclassifications and detection heterogeneity.  Methods in Ecology and Evolution.
  * Conn et al. (2013) Conn, P. B., McClintock, B. T., Cameron, M. F., Johnson, D. S., Moreland, E. E. and Boveng, P. L. (2013) Accommodating species identification errors in transect surveys.  Ecology, 94, 2607–2618.
  * Elith and Leathwick (2009) Elith, J. and Leathwick, J. R. (2009) Species distribution models: ecological explanation and prediction across space and time.  Annual Review of Ecology, Evolution and Systematics, 40, 677–697.
  * Fahrmeir et al. (2013) Fahrmeir, L., Kneib, T., Lang, S. and Marx, B. (2013) Categorical Regression Models, 325–347.  Berlin, Heidelberg: Springer Berlin Heidelberg.  URL: https://doi.org/10.1007/978-3-642-34333-9_6.
  * Ferguson et al. (2015) Ferguson, P. F., Conroy, M. J. and Hepinstall-Cymerman, J. (2015) Occupancy models for data with false positive and false negative errors and heterogeneity across sites and surveys.  Methods in Ecology and Evolution, 6, 1395–1406.
  * GBIF.Org (2022) GBIF.Org (2022) Gbif: The global biodiversity information facility (year) what is gbif?  URL: https://www.gbif.org/what-is-gbif.
  * 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, 1081–1091.
  * Hijmans et al. (2015) Hijmans, R. J., Van Etten, J., Cheng, J., Mattiuzzi, M., Sumner, M., Greenberg, J. A., Lamigueiro, O. P., Bevan, A., Racine, E. B., Shortridge, A. et al. (2015) Package ‘raster’.  R package, 734.
  * Hoekman (2021) Hoekman, S. T. (2021) Multi-observer methods for estimating uncertain species identification.  Ecosphere, 12, e03648.
  * Hui et al. (2015) Hui, F. K., Warton, D. I. and Foster, S. D. (2015) Multi-species distribution modeling using penalized mixture of regressions.  The Annals of Applied Statistics, 9, 866–882.
  * Johnston et al. (2018) Johnston, A., Fink, D., Hochachka, W. M. and Kelling, S. (2018) Estimates of observer expertise improve species distributions from citizen science data.  Methods in Ecology and Evolution, 9, 88–97.
  * Johnston et al. (2022) Johnston, A., Matechou, E. and Dennis, E. B. (2022) Outstanding challenges and future directions for biodiversity monitoring using citizen science data.  Methods in Ecology and Evolution.
  * Jongbloed (2016) Jongbloed, R. H. (2016) Flight height of seabirds: a literature study.
  * Kelling et al. (2015) Kelling, S., Johnston, A., Hochachka, W. M., Iliff, M., Fink, D., Gerbracht, J., Lagoze, C., La Sorte, F. A., Moore, T., Wiggins, A. et al. (2015) Can observation skills of citizen scientists be estimated using species accumulation curves?  PloS one, 10, e0139600.
  * Kéry and Royle (2020) Kéry, M. and Royle, J. A. (2020) Applied Hierarchical Modeling in Ecology: Analysis of distribution, abundance and species richness in R and BUGS: Volume 2: Dynamic and Advanced Models.  Academic Press.
  * Keshavan et al. (2019) Keshavan, A., Yeatman, J. D. and Rokem, A. (2019) Combining citizen science and deep learning to amplify expertise in neuroimaging.  Frontiers in neuroinformatics, 13, 29.
  * Koch et al. (2022) Koch, W., Hogeweg, L., Nilsen, E. B. and Finstad, A. G. (2022) Maximizing citizen scientists’ contribution to automated species recognition.  Scientific reports, 12, 1–10.
  * Langenkämper et al. (2019) Langenkämper, D., Simon-Lledo, E., Hosking, B., Jones, D. O. and Nattkemper, T. W. (2019) On the impact of citizen science-derived data quality on deep learning based classification in marine images.  PloS One, 14, e0218086.
  * Lotfian et al. (2021) Lotfian, M., Ingensand, J. and Brovelli, M. A. (2021) The partnership of citizen science and machine learning: Benefits, risks, and future challenges for engagement, data collection, and data quality.  Sustainability, 13, 8087.
  * Louvrier et al. (2018) Louvrier, J., Chambert, T., Marboutin, E. and Gimenez, O. (2018) Accounting for misidentification and heterogeneity in occupancy studies using hidden markov models.  Ecological modelling, 387, 61–69.
  * Matheson (2014) Matheson, C. A. (2014) inaturalist.  Reference Reviews.
  * 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, e65808.
  * Molinari-Jobin et al. (2012) Molinari-Jobin, A., Kéry, M., Marboutin, E., Molinari, P., Koren, I., Fuxjäger, C., Breitenmoser-Würsten, C., Wölfl, S., Fasel, M., Kos, I. et al. (2012) Monitoring in the presence of species misidentification: the case of the e urasian lynx in the a lps.  Animal Conservation, 15, 266–273.
  * Montesinos López et al. (2022) Montesinos López, O. A., Montesinos López, A. and Crossa, J. (2022) Overfitting, model tuning, and evaluation of prediction performance.  In Multivariate statistical machine learning methods for genomic prediction, 109–139. Springer.
  * O’Hara and Sillanpää (2009) O’Hara, R. B. and Sillanpää, M. J. (2009) A review of bayesian variable selection methods: what, how and which.  Bayesian analysis, 4, 85–117.
  * Ovaskainen and Abrego (2020) Ovaskainen, O. and Abrego, N. (2020) Evaluating Model Fit and Selecting among Multiple Models, 217–252.  Ecology, Biodiversity and Conservation. Cambridge University Press.
  * Ovaskainen and Soininen (2011) Ovaskainen, O. and Soininen, J. (2011) Making more out of sparse data: hierarchical modeling of species communities.  Ecology, 92, 289–295.
  * Pollock et al. (2014) Pollock, L. J., Tingley, R., Morris, W. K., Golding, N., O’Hara, R. B., Parris, K. M., Vesk, P. A. and McCarthy, M. A. (2014) Understanding co-occurrence by modelling species simultaneously with a joint species distribution model (jsdm).  Methods in Ecology and Evolution, 5, 397–406.
  * Ponti and Seredko (2022) Ponti, M. and Seredko, A. (2022) Human-machine-learning integration and task allocation in citizen science.  Humanities and Social Sciences Communications, 9, 1–15.
  * R Core Team (2022) R Core Team (2022) R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.  URL: https://www.R-project.org/.
  * Renner et al. (2015) Renner, I. W., Elith, J., Baddeley, A., Fithian, W., Hastie, T., Phillips, S. J., Popovic, G. and Warton, D. I. (2015) Point process models for presence-only analysis.  Methods in Ecology and Evolution, 6, 366–379.
  * Roberts et al. (2017) Roberts, D. R., Bahn, V., Ciuti, S., Boyce, M. S., Elith, J., Guillera-Arroita, G., Hauenstein, S., Lahoz-Monfort, J. J., Schröder, B., Thuiller, W. et al. (2017) Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure.  Ecography, 40, 913–929.
  * Royle and Nichols (2003) Royle, J. A. and Nichols, J. D. (2003) Estimating abundance from repeated presence–absence data or point counts.  Ecology, 84, 777–790.
  * Santos-Fernandez and Mengersen (2021) Santos-Fernandez, E. and Mengersen, K. (2021) Understanding the reliability of citizen science observational data using item response models.  Methods in Ecology and Evolution, 12, 1533–1548.
  * Saoud et al. (2020) Saoud, Z., Fontaine, C., Loïs, G., Julliard, R. and Rakotoniaina, I. (2020) Miss-identification detection in citizen science platform for biodiversity monitoring using machine learning.  Ecological Informatics, 60, 101176.
  * 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.
  * Steel et al. (1953) Steel, G. et al. (1953) Relation between poisson and multinomial distributions.
  * Sutherland et al. (2013) Sutherland, C., Elston, D. and Lambin, X. (2013) Accounting for false positive detection error induced by transient individuals.  Wildlife Research, 40, 490–498.
  * Suzuki-Ohno et al. (2022) Suzuki-Ohno, Y., Westfechtel, T., Yokoyama, J., Ohno, K., Nakashizuka, T., Kawata, M. and Okatani, T. (2022) Deep learning increases the availability of organism photographs taken by citizens in citizen science programs.  Scientific Reports, 12, 1–10.
  * Tang et al. (2021) Tang, B., Clark, J. S. and Gelfand, A. E. (2021) Modeling spatially biased citizen science effort through the ebird database.  Environmental and Ecological Statistics, 28, 609–630.
  * Tobler et al. (2019) Tobler, M. W., Kéry, M., Hui, F. K., Guillera-Arroita, G., Knaus, P. and Sattler, T. (2019) Joint species distribution models with species correlations and imperfect detection.  Ecology, 100, e02754.
  * Tulloch et al. (2013) Tulloch, A. I., Possingham, H. P., Joseph, L. N., Szabo, J. and Martin, T. G. (2013) Realising the full potential of citizen science monitoring programs.  Biological Conservation, 165, 128–138.
  * Ueda (2020) Ueda, K. (2020) inaturalist research-grade observations. inaturalist. org.  Occurrence dataset https://doi. org/10.15468/ab3s5x [accessed via GBIF. org on 26 May 2020, https://www. gbif. org/occurrence/2294559370].
  * de Valpine et al. (2017) de Valpine, P., Turek, D., Paciorek, C. J., Anderson-Bergman, C., Lang, D. T. and Bodik, R. (2017) Programming with models: writing statistical algorithms for general model structures with nimble.  Journal of Computational and Graphical Statistics, 26, 403–413.
  * Vermeiren et al. (2020) Vermeiren, P., Reichert, P. and Schuwirth, N. (2020) Integrating uncertain prior knowledge regarding ecological preferences into multi-species distribution models: Effects of model complexity on predictive performance.  Ecological Modelling, 420, 108956.
  * Vohland et al. (2021) Vohland, K., Land-Zandstra, A., Ceccaroni, L., Lemmens, R., Perelló, J., Ponti, M., Samson, R. and Wagenknecht, K. (2021) The science of citizen science.  Springer Nature.
  * Willi et al. (2019) Willi, M., Pitman, R. T., Cardoso, A. W., Locke, C., Swanson, A., Boyer, A., Veldthuis, M. and Fortson, L. (2019) Identifying animal species in camera trap images using deep learning and citizen science.  Methods in Ecology and Evolution, 10, 80–91.
  * 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, 71–81.
  * Yu et al. (2015) Yu, H., Jiang, S. and Land, K. C. (2015) Multicollinearity in hierarchical linear models.  Social science research, 53, 118–136.