# Frequentist Prediction Sets for Species Abundance using Indirect Information

Elizabeth Bersson and Peter D. Hoff
Department of Statistical Science
Duke University

###### Abstract

Citizen science databases that consist of volunteer-led sampling efforts of
species communities are relied on as essential sources of data in ecology.
Summarizing such data across counties with frequentist-valid prediction sets
for each county provides an interpretable comparison across counties of
varying size or composition. As citizen science data often feature unequal
sampling efforts across a spatial domain, prediction sets constructed with
indirect methods that share information across counties may be used to improve
precision. In this article, we present a nonparametric framework to obtain
precise prediction sets for a multinomial random sample based on indirect
information that maintain frequentist coverage guarantees for each county. We
detail a simple algorithm to obtain prediction sets for each county using
indirect information where the computation time does not depend on the sample
size and scales nicely with the number of species considered. The indirect
information may be estimated by a proposed empirical Bayes procedure based on
information from auxiliary data. Our approach makes inference for under-
sampled counties more precise, while maintaining area-specific frequentist
validity for each county. Our method is used to provide a useful description
of avian species abundance in North Carolina, USA based on citizen science
data from the eBird database.

Key words: categorical data, conformal prediction, empirical Bayes,
exchangeability, frequentist coverage, nonparametric

##  1 Introduction

Understanding species abundance across heterogeneous spatial areas is an
important task in ecology. Citizen science databases that consist of
observations of species counts gathered by volunteers are increasingly
regarded as one of the richest sources of data for such a task. One of the
largest such data sources is the eBird database in which citizen scientists
throughout the world input counts of bird sightings (Sullivan et al., 2009).
In addition to its use for describing avian species abundance, eBird is a
principal resource for understanding global biodiversity and is widely used in
constructing and implementing conservation action plans (Sullivan et al.,
2017).

More generally, analyses from such databases may be used for informing policy,
conservation efforts, habitat preservation, and more, for which understanding
species prevalence for non-overlapping geographic areas, such as counties
across a state or country, is important. In practice, species abundance from
citizen science data are commonly summarised within areas such as counties by
empirical proportions from a sample, as in, e.g., Arnold et al. (2021);
Camerini and Groppali (2014). Such proportions can be used to construct a
prediction set for each county that provides a description of species
prevalence for that county with guaranteed frequentist coverage.

Given the impact on policy design, corresponding uncertainty quantification is
of particular import (Lele, 2020), and so it is desirable that precise
prediction sets maintain a target coverage rate regardless of the county’s
size or composition. This is challenging as a common feature of citizen
science data is unequal sampling efforts that results in some counties with
large amounts of data information and others with very little. Using direct
procedures that only make use of within-county information, a prediction set
may be imprecise in these counties with low sampling efforts. This suggests
using indirect information such as data from neighboring counties to improve
prediction set precision for a given county.

In this article, we describe species abundance across sampling areas such as
counties with frequentist-valid prediction sets that are constructed to
contain an unobserved bird with $1-\alpha$ probability. That is, a valid
prediction set for a given county is a set of avian species such that an
unobserved bird will belong to one of those species with $1-\alpha$
probability in a frequentist sense. We develop a valid nonparametric
prediction method that allows for information to be shared across counties.
Specifically, our approach results in prediction sets with guaranteed
frequentist coverage for each county that are constructed with the
incorporation of indirect or prior information. We detail and provide code for
an empirical Bayes procedure to estimate such prior information from auxiliary
data such as neighboring counties. If this indirect information used to
construct the prediction sets is accurate, the prediction sets will be smaller
than direct sets that only make use of within-county information.

In Section 4, we detail the usefulness of the proposed approach in summarising
the eBird citizen science data. Sharing information across counties generally
results in smaller prediction sets as compared to direct prediction
approaches, particularly so in counties with low sampling efforts. Moreover,
the prediction sets provide a useful summary of the data that may be used to
compare information across areas and better inform policy.

##  2 Methodology

###  2.1 Background and Notation

For county $j\in\\{1,...,J\\}$, let $\boldsymbol{X}_{j}$ be a vector of length
$K$ where $X_{j,i}=x_{j,i}$ is the observed count of species $i$ over some set
sampling period that may vary across counties. We model $\boldsymbol{X}_{j}$
with a $K$-dimensional multinomial distribution with
$N_{j}=\sum_{i=1}^{K}x_{j,i}$ trials and population proportions vector
${\boldsymbol{\theta}}_{j}$,

$\boldsymbol{X}_{j}\sim MN_{K}({\boldsymbol{\theta}}_{j},N_{j}).$ (1)

We construct a prediction set $A_{\alpha}(\boldsymbol{x}_{j})$ for an
observation of a new bird arising from the same distribution,
$\boldsymbol{Y}_{j}\sim MN_{K}({\boldsymbol{\theta}}_{j},1)$ where
$\boldsymbol{Y}_{j}\in\mathcal{Y}$ for
$\mathcal{Y}=\\{(y_{1},...,y_{K}):\sum_{i=1}^{K}y_{i}=1,y_{i}\in\\{0,1\\}(i=1,....,K)\\}$.
Let $\boldsymbol{y}_{j}^{(k)}\in\mathcal{Y}$ denote a prediction of category
$k$, that is, let $\boldsymbol{y}_{j}^{(k)}$ be a vector of length $K$ with a
one at index $k$ and zeros elsewhere. In particular, we are interested in a
prediction set for $\boldsymbol{Y}_{j}$ that maintains frequentist validity
for some error rate $\alpha$. Formally, we refer to this as an $\alpha$-valid
prediction set:

######  Definition 1 ($\alpha$-Valid Prediction Set)

An $\alpha$-valid prediction set for a predictand
$\boldsymbol{Y}_{j}\in\mathcal{Y}$ is any subset $A_{\alpha}$ of the sample
space $\mathcal{Y}$ that contains $\boldsymbol{Y}_{j}$ with probability
greater than or equal to $1-\alpha$,

$P_{\theta}\left(\boldsymbol{Y}_{j}\in
A_{\alpha}(\boldsymbol{x}_{j})\right)\geq 1-\alpha,\quad\forall\;\theta,$ (2)

where the probability is taken with respect to $\boldsymbol{Y}_{j}$ and
$\boldsymbol{X}_{j}$.

Additionally, small or precise $\alpha$-valid prediction sets are of
particular interest, where prediction set size is measured by expected
cardinality, that is, expected number of the $K$ categories in the sample
space included in the prediction set.

###  2.2 Order-based prediction for a single area

A standard approach to construct $\alpha$-valid prediction sets for each
county or area is with a direct method that only makes use of within-area
information. As such, we first consider construction of a prediction set for a
single area $j$, using only data from county $j$. For ease of notation, we
drop the area-identifying subscript in this subsection.

For multinomial data in general, if the event probability vector
${\boldsymbol{\theta}}$ is known, an $\alpha$-valid prediction set is any
combination of categories such that their event probabilities cumulatively sum
to be greater than or equal to $1-\alpha$. Equivalently stated, an
$\alpha$-valid prediction set may be constructed by excluding categories such
that the cumulative sum of the excluded categories’ event probabilities is
less than $\alpha$. Such a prediction set may be constructed by admitting
categories in some prespecified order into the prediction set until the
cumulative sum of their event probabilities is at least $1-\alpha$. The
resulting prediction set will have $1-\alpha$ coverage regardless of the
ordering used to admit categories. In fact, the class of all $\alpha$-valid
prediction sets may be constructed by following this procedure for non-strict
total orderings of categories.

Perhaps intuitively, constructing such a prediction set by including
categories with the largest event probabilities will result in the smallest
$\alpha$-valid prediction set. In the terminology of ordering, this
corresponds to constructing a prediction set based on an ordering of
categories that matches the ordering of the elements in
${\boldsymbol{\theta}}$. We refer to this optimal ordering as the oracle
ordering:

######  Theorem 1 (Oracle order-based prediction)

Let $\boldsymbol{Y}\sim MN_{K}({\boldsymbol{\theta}},1)$ for
${\boldsymbol{\theta}}$ known. Then,

  1. (a)

the class of all $\alpha$-valid prediction sets for a given
${\boldsymbol{\theta}}$ consists of prediction sets of the form,

$A_{\alpha}^{\theta,o}=\left\\{\boldsymbol{y}^{(k)}\in\mathcal{Y}:\left[\sum_{l=1}^{K}\mathbf{1}\left(o_{k}\geq
o_{l}\right)\theta_{l}\right]>\alpha\right\\},$ (3)

for some vector $\boldsymbol{o}\in\mathbb{R}^{K}$, and

  2. (b)

the oracle ordering is that which corresponds to the increasing order
statistics of ${\boldsymbol{\theta}}$,

$\boldsymbol{o}^{\theta}=\left\\{\boldsymbol{o}:\theta_{m}<\theta_{n}\Rightarrow
o_{m}<o_{n}\;\forall\;m,n\in\\{1,....,K\\},m\neq n\right\\},$

and $A_{\alpha}^{\theta,o^{\theta}}$ has the smallest cardinality among all
orderings.

In practice, ${\boldsymbol{\theta}}$ is unknown, but a prediction set may be
constructed based on an observed sample $\boldsymbol{X}=\boldsymbol{x}$. It
turns out, in fact, that any conditional $\alpha$-valid prediction set can be
written similarly to the previous construction (Equation 3) where the
cumulative sum is computed with respect to the empirical proportions given by
$\boldsymbol{x}$ and $\boldsymbol{y}$. This is a generalization of the
conformal prediction framework, a popular machine learning approach to
construct prediction regions based on measuring conformity (or non-conformity)
of a predictand to an observed sample (Vovk et al., 2005).

######  Theorem 2 ($\alpha$-valid order-based prediction)

Let $\boldsymbol{X}\sim MN_{K}({\boldsymbol{\theta}},N),\boldsymbol{Y}\sim
MN_{K}({\boldsymbol{\theta}},1)$. Then, every conformal $\alpha$-valid
prediction set based on observed data $\boldsymbol{x}$ can be written

$A_{\alpha}(\boldsymbol{x})=\left\\{\boldsymbol{y}^{(k)}\in\mathcal{Y}:\left[\sum_{l=1}^{K}\mathbbm{1}\left(o_{k}\geq
o_{l}\right)\frac{x_{l}+y_{l}^{(k)}}{N+1}\right]>\alpha\right\\},$ (4)

for some vector $\boldsymbol{o}\in\mathbb{R}^{K}$.

Note that the prediction set depends on the vector $\boldsymbol{o}$ only
through the order of its elements.

For any ordering of the $K$ categories, constructing a prediction set
following Theorem 2 results in a prediction set with guaranteed finite-sample
$1-\alpha$ frequentist coverage. The choice of ordering, however, will impact
prediction set precision, that is, the set’s cardinality. For inference for a
single area, a natural approach is to order the categories with respect to
their empirical proportions. The empirical proportions are unbiased for
population proportions, so, if the area has a large sample size, an ordering
based on the empirical proportions will approximate the oracle ordering well.
It turns out this approach is well-motivated by classical prediction
approaches. Specifically, a standard direct prediction method constructs a
prediction set separately for an area based on an area-specific conditional
pivotal quantity (Faulkenberry, 1973; Tian et al., 2022). For a multinomial
population, $\boldsymbol{Y}|\boldsymbol{X}+\boldsymbol{Y}$ is such a quantity
that follows a multivariate hypergeometric distribution which does not depend
on the event probability vector. See Thatcher (1964) for work on prediction
sets of this type for binomial data. A prediction set constructed to contain
species belonging to a highest mass region of this pivotal distribution is
obtained by including species with the largest empirical counts until their
cumulative proportion sum exceeds $1-\alpha$,

$\displaystyle A_{\alpha}^{D}{}$
$\displaystyle(\boldsymbol{x})=\Big{\\{}\boldsymbol{y}^{(k)}\in\mathcal{Y}:$
(5)
$\displaystyle\left[\sum_{l=1}^{K}\mathbbm{1}\left(\left(x_{k}+y_{k}^{(k)}\right)\geq\left(x_{l}+y^{(k)}_{l}\right)\right)\frac{x_{l}+y_{l}^{(k)}}{N+1}\right]>\alpha\Big{\\}}.$

This direct prediction set based on an ordering of the empirical proportions
is appealing as it is easy to interpret and has finite-sample guaranteed
$1-\alpha$ frequentist coverage. For an area with low sampling effort, though,
the empirical proportions will not precisely estimate the true proportions. As
a result, a prediction set may have prohibitively large cardinality such that
it is not practically useful. For such an area, incorporating indirect
information from neighboring counties can improve the estimates of the county
proportions and thereby increase the precision of a prediction set.

###  2.3 Order-based prediction for multiple areas

In general, in analyzing small area data, that is, areal data featuring small
within-area sample sizes in some areas, it is common to utilize indirect
methods that share information across areas (Rao and Molina, 2015). The eBird
database is a rich data source, and inference in any given county may be
improved upon by taking advantage of auxiliary data using an indirect method.
In this subsection, we detail how information from neighboring counties may be
used in estimating an ordering of categories to improve prediction set
precision.

As opposed to a direct prediction set based on an ordering corresponding to
within-county empirical proportions, an indirect prediction set can be
constructed similarly whereby species are admitted into the prediction set
based on an ordering corresponding to empirical posterior proportions
estimated from a hierarchical model. Such an estimate may be obtained based on
a conjugate Dirichlet prior distribution parameterized with a common
concentration hyperparameter for the $J$ areas,

${\boldsymbol{\theta}}_{1},...,{\boldsymbol{\theta}}_{J}\sim{Dirichlet}_{K}(\boldsymbol{\gamma}).$
(6)

Given a hyperparameter $\boldsymbol{\gamma}\in\mathbb{R}^{K}$, the posterior
expectation of the proportions ${\boldsymbol{\theta}}_{j}$ in county $j$ is
$\tilde{\boldsymbol{x}}_{j}/(N_{j}+\sum_{i=1}^{K}\gamma_{i})$ where
$\tilde{\boldsymbol{x}}_{j}=\boldsymbol{x}_{j}+\boldsymbol{\gamma}$. In this
way, $\tilde{\boldsymbol{x}}_{j}$ may be interpreted as a posterior vector of
counts for county $j$. Then, an $\alpha$-valid prediction set based on
$\tilde{\boldsymbol{x}}_{j}$ is,

$\displaystyle A_{\alpha}^{I}{}$
$\displaystyle(\boldsymbol{x}_{j})=\Big{\\{}\boldsymbol{y}^{(k)}\in\mathcal{Y}:$
(7)
$\displaystyle\left[\sum_{l=1}^{K}\mathbbm{1}\left(\left(\tilde{x}_{j,k}+y^{(k)}_{k}\right)\geq\left(\tilde{x}_{j,l}+y^{(k)}_{l}\right)\right)\frac{x_{j,l}+y^{(k)}_{l}}{N_{j}+1}\right]>\alpha\Big{\\}}.$

By Theorem 2, $A^{I}_{\alpha}(\boldsymbol{X}_{j})$ is an $\alpha$-valid
procedure, and it is constructed based on prior information. Specifically, it
differs from the direct set given in Equation 5 in that categories are
admitted into the prediction set based on an ordering determined by posterior
counts that incorporate indirect information $\boldsymbol{\gamma}$, as opposed
to an ordering based on the observed sample. Moreover, it has been shown that
if the indirect information used is accurate,
$A^{I}_{\alpha}(\boldsymbol{X}_{j})$ may be more precise than a direct
prediction set with the same coverage rate (Hoff, 2023; Bersson and Hoff,
2022).

In total, $A^{D}_{\alpha}(\boldsymbol{X}_{j})$ and
$A^{I}_{\alpha}(\boldsymbol{X}_{j})$ are both $\alpha$-valid prediction
procedures. They differ in the order in which species are admitted into the
prediction sets, as species are admitted into the direct set in terms of
decreasing empirical proportions and into the indirect set in terms of
decreasing posterior counts. As a result, for an area with a small sample
size, incorporating accurate prior information can result in an ordering used
to construct a prediction set that more accurately approximates the oracle
ordering as the empirical proportions might be too unstable. Of note, these
two approaches are equivalent for a uniform prior
$\boldsymbol{\gamma}=c\mathbf{1}$, for any constant $c$. This includes, for
example, a standard noninformative prior $c=1$, a standard objective Bayes
Jeffrey’s prior $c=1/2$, and an improper prior $c=0$.

###  2.4 Empirical Bayes estimation of indirect information

To obtain an $\alpha$-valid indirect prediction set for county
$j\in\\{1,...,J\\}$, all that is required is an estimate of the prior
concentration parameter $\boldsymbol{\gamma}$. We propose an empirical
Bayesian approach whereby values of $\boldsymbol{\gamma}$ to be used for
county $j$ are estimated from data collected in neighboring counties.
Specifically, we use the maximum likelihood estimate of the marginal
likelihood based on the conjugate hierarchical model given by Equations 1 and
6,

$\displaystyle\boldsymbol{\gamma}_{j}={}$
$\displaystyle\arg\max_{\boldsymbol{\gamma}}\log p\left(\bigcup_{l\in
L}\boldsymbol{X}_{l}\Big{|}\boldsymbol{\gamma}\right)$ (8) $\displaystyle={}$
$\displaystyle\arg\max_{\boldsymbol{\gamma}}\log\prod_{l\in
L}\left[\frac{\Gamma(\sum_{i=1}^{K}\gamma_{i})}{\Gamma(\sum_{i=1}^{K}x_{l,i}+\gamma_{i})}\times\prod_{i=1}^{K}\frac{\Gamma(x_{l,i}+\gamma_{i})}{\Gamma(\gamma_{i})}\right],$

where $L_{j}\subseteq\\{1,...,K\\}\backslash\\{j\\}$ is a non-empty set
containing the indices of counties neighboring county $j$. Information is
shared across neighboring counties to inform an estimate of the prior for
county $j$, and, when estimated in this way, the prior concentration
represents an across-county pooled prior concentration. This optimization
problem can be solved numerically with a Newton-Raphson algorithm. See
Appendix A for details and derivation of such an algorithm. Code to implement
this procedure in the R Statistical Programming language is available online,
see Section 5.

When $\boldsymbol{\gamma}_{j}$ is estimated using data independent of area $j$
and used to construct $A^{I}_{\alpha}(\boldsymbol{X}_{j})$, the finite sample
coverage guarantee of $A^{I}_{\alpha}(\boldsymbol{X}_{j})$ holds regardless of
the accuracy of the estimated prior hyperparameter. If the estimated vector
$\boldsymbol{\gamma}_{j}$ is accurate, then
$A^{I}_{\alpha}(\boldsymbol{X}_{j})$ may also be more precise than direct
prediction approaches.

##  3 Simulation Study

Figure 1: Monte Carlo approximations (+/- 1 standard deviation) of the
expected cardinality ratio of (red) indirect to direct methods and (blue)
$\alpha$-valid prediction set given the oracle ordering to direct method.

To illustrate how the incorporation of indirect information can affect
precision of prediction sets, we compare expected set cardinality obtained
from the indirect and direct prediction methods for a single simulated area.
In contrast to the eBird data, for example, the analysis of this section
corresponds to that of one county. Because citizen science data such as these
often feature unequal sampling efforts across counties, we are particularly
interested in demonstrating the difference in cardinality between these two
approaches for a range of sample sizes $N=10,100,1000$. Moreover, we compare
results for varying number of categories $K$. Throughout, we consider a low
entropy regime in which $\lceil K/4\rceil$ categories unequally split nearly
all of the probability mass, and the rest of the categories have nearly
probability $0$. While we do not necessarily expect real populations in
practice to have such a distribution, it is chosen to clearly demonstrate the
benefit of including indirect information in the construction of prediction
sets that maintain frequentist coverage.

In one construction of indirect prediction sets, we consider a prior based on
full information with moderate prior precision
$\boldsymbol{\gamma}={\boldsymbol{\theta}}\times 10$. We compare with direct
prediction sets given by Equation 5, or, equivalently, indirect prediction
sets constructed with a uniform prior $\boldsymbol{\gamma}=c\mathbf{1}$.
Finally, we compare the approaches to $\alpha$-valid order-based prediction
sets obtained based on an oracle ordering. Results comparing Monte Carlo
approximations of the expected prediction set cardinality ratios between the
various approaches obtained from $25,000$ replications are displayed in Figure
1.

As all methods considered are $\alpha$-valid procedures, the crucial
difference between them is the incorporation of indirect information.
Utilizing accurate prior information in the construction of prediction sets
generally results in prediction sets distinctly smaller than direct sets,
particularly so if there are a large number of categories relative to the
sample size. This is evidenced by the red dashes in Figure 1 showing the
expected cardinality ratios of the indirect to direct prediction sets are
always at or below a value of 1. An accurate prior may be one that
approximates the true probability mass vector well with large precision
relative to sample size, as seen in the left plot of Figure 1 for sample size
$N=10$. More generally, though, all that is needed is a prior that results in
posterior counts that accurately approximate the oracle ordering of
categories. We discuss the three sample size regimes in detail below.

For a small sample size of $N=10$, the prior $\boldsymbol{\gamma}$ used to
construct the indirect prediction sets is an informative prior with strong
precision in that the scale used is equal to the sample size in this case. As
a result, the posterior distributions contain notably more information than
what is in each simulated dataset. As a result, the ordering of categories
induced by the posterior counts, used to construct the indirect prediction
sets, are accurately approximating the oracle ordering of categories. This is
evidenced by the nearly identical behavior of the two cardinality ratios
explored. In conjunction with the instability of the direct method in the
presence of such a small sample size, this results in notably smaller
cardinality of the indirect set as compared to the direct set, even for
relatively small total numbers of categories. At its best, the indirect
prediction set is about 80% smaller than the direct set.

For a moderate sample size of $N=100$, the prior precision used to construct
the indirect prediction sets is not overwhelming as compared to the sample
size, and hence the posterior counts do not approximate the oracle ordering as
well as in the regime with a smaller sample size. This is evidenced by the
divergence of the red and blue dashes in the middle plot of Figure 1. Still,
particularly as the number of categories increases for fixed $N$, the benefit
of utilizing prior information of this type is highlighted by the decline of
the cardinality ratio of the indirect to direct prediction sets (red lines).
For example, in the case of $N=100$ and $K=150$, the indirect prediction set
constructed with $\boldsymbol{\gamma}$ is about $15\%$ smaller than the direct
prediction set.

A similar but less pronounced pattern is seen in the presence of a larger
sample size of $N=1000$. For this sample size with $K\leq 150$, all methods
considered perform relatively similarly. However, as the number of categories
increases, there is a distinct gain in prediction set precision given the
input of indirect information in prediction set construction.

##  4 Summarizing eBird species abundance data

In this section, we describe avian species abundance in North Carolina, USA
from eBird data obtained from citizen-uploaded complete checklists of species
observations in the first week of May 2023. Across the $99$ counties, $393$
unique species were identified. Some species such as the Northern Cardinal,
Carolina Wren, and American Robin were identified frequently. Many others like
the Northern Saw-whet Owl and the Solitary Sandpiper were rarely seen; in
fact, $50\%$ of species were seen fewer than $100$ times each across the
entire state. Moreover, within-county sample sizes vary drastically (Figure 2)
from approximately $50,000$ individual birds identified in Wake County, one of
the most populous counties in NC that contains the state’s capital, to only
$8$ in Pasquotank County, a small coastal county consisting of about $1/30$th
of the human population of Wake County.

As motivated in the Introduction, describing such data with $\alpha$-valid
prediction sets for each county provides a useful summary with unambiguous
statistical interpretation. That is, with at least probability $1-\alpha$, an
unobserved bird in a given county will belong to a species contained in the
specified prediction set, where the probability is taken with respect to the
random sample and the predictand. Here, we demonstrate the usefulness of this
approach in gaining better understanding of species abundance. Moreover, we
elaborate on the benefit of utilizing indirect information in the construction
of practically useful sets that are precise, particularly for counties with
small within-county sample sizes.

Figure 2: Within-county log sample size. Figure 3: Cardinality ratio of
indirect to direct prediction sets. Prior hyperparameters estimated with an
empirical Bayesian procedure based on five nearest neighbors for each county.
The lowest quantile sample sizes are overlaid on their respective counties.

For each county in NC, we construct an indirect prediction set based on a
prior hyperparameter estimated from data in the five nearest neighboring
counties, following the procedure described in Section 2.4. The eBird data
consist of independent samples collected across the state, so samples are
independent across counties. As a result of this independence, finite-sample
coverage of the indirect prediction approach is guaranteed. We compare the
cardinality of these indirect prediction sets to that of direct prediction
sets, both of which maintain at least $95\%$ coverage for each county. The
cardinality ratios of the indirect to direct prediction sets across the
counties in NC are plotted in Figure 3. To highlight the impact of within-
county sample size, the lower quantile sample sizes are overlaid on their
respective county.

In general, the incorporation of indirect information in the construction of
prediction sets results in notably smaller cardinality of the indirect
prediction sets as compared with that of the direct prediction sets. Of the 99
counties in NC, indirect sets have smaller cardinality in 65, and the two
approaches result in the same cardinality in 20 counties. The improvement in
cardinality is particularly conspicuous in counties with small to moderate
sample sizes, as evidenced by the sample sizes of counties with the brightest
shade of red in Figure 3. Moreover, ten counties have trivial direct sets
consisting of all $K$ species, while only two counties with the smallest
within-county sample sizes, 8 and 14, have trivial indirect prediction sets.
For the county with the third smallest sample size (24), the indirect
prediction set only includes 80 species, or about 20% of all possible species,
while the direct prediction set is the trivial set.

Overall, even in counties with larger sample sizes, it is most common for the
indirect and direct prediction sets to contain a different set of species. In
fact, the indirect and direct prediction sets disagree for nearly every county
in NC. They are equivalent for only six counties where they aren’t both
trivial sets. Commonly, this discrepancy corresponds with smaller indirect
sets, and hence highlights the benefit of inclusion of indirect information in
the construction of prediction sets.

Figure 4: Empirical probability masses of species included in only the
indirect (red text), only the direct (blue text), or both (purple text)
prediction sets, sorted by the posterior proportion. MLE plotted in blue, and
the posterior proportion based on a prior estimated from data in nearest five
neighboring counties.

###  4.1 Order-based prediction in Robeson County

| D. Cormorant | E. Kingbird | Pine Warbler | C. Sparrow
---|---|---|---|---
Robeson | 0.81% | 0.4% | 0.00% | 0.00%
NC-017 | 0.00% | 2.85% | 1.97% | 2.63%
NC-047 | 0.00% | 0.00% | 1.29% | 0.64%
NC-051 | 0.00% | 0.68% | 2.62% | 5.59%
NC-093 | 0.00% | 0.00% | 1.09% | 0.55%
NC-165 | 0.00% | 0.86% | 2.58% | 5.16%
$\boldsymbol{\gamma}$ | 0.00 | 1.33 | 3.42 | 4.69
Table 1: Percentage of all birds observed within each respective county, for
select species included in either the indirect set (red text) or the direct
set (blue text). Estimated prior hyperparameter $\boldsymbol{\gamma}$ for
Robeson County recorded in the last row.

To further compare the two approaches and elucidate the role of the ordering
of the species, we elaborate on the construction of indirect and direct
prediction sets for Robeson County. Robeson is located near the southeastern
border of NC and features a moderately small within-county sample size of
$247$ birds observed, with species-specific observation counts ranging from
zero to ten. The two prediction sets have nearly the same cardinality but
contain differences in species inclusion. Specifically, the indirect
prediction set contains $33$ species, and the direct set contains $32$, with
an overlap of $27$ species.

To illustrate the role of the ordering used in the construction of
$\alpha$-valid prediction sets, the empirical proportions based on the
observed sample (MLE) and posterior proportions (Post.Pred) are plotted in
Figure 4 for the union of included species in the two sets. In the figure, the
species are sorted by increasing posterior proportions. The indirect and
direct sets include species based on the posterior and empirical
distributions, respectively. Discrepancies between the indirect and direct
sets occur when these two distributions disagree. From Figure 4, it is easy to
see the indirect prediction set consists of the species with the 33 largest
posterior predictive proportions. In contrast, the direct set consists of
species with the largest sample probability mass. Naturally, the ordering of
these two estimates agree for species common to the region, and, as such,
there is a fair amount of overlap of species inclusion.

As a result of our estimation procedure for the prior hyperparameter
$\boldsymbol{\gamma}$ for Robeson County, the disparity between inclusion or
exclusion of a species among the two prediction set methods is further
elucidated by examining species presence in neighboring counties. In short,
species with more frequent occurrence in neighboring counties will have a
larger estimated prior count than those seen rarely in neighboring counties.
Species occurrences in neighboring counties are displayed in Table 1 for a
select few species along with the estimated $\boldsymbol{\gamma}$ for Robeson
County, obtained by solving Equation 8 using data in these neighboring
counties.

Intuitively, species that are seen in neighboring counties with some relative
frequency, such as the Chipping Sparrow or Pine Warbler, are probably also
present in Robeson County, and hence should be included in a prediction set.
In practice, these species have a comparatively high estimated prior of about
5 and 4, respectively, and hence are included in the indirect prediction set
even though they weren’t recorded as being observed in Robeson County in the
dataset. Alternatively, consideration of indirect information yields the
conclusion that species like the Eastern Kingbird and Cormorant may be rare in
the area in general, as reflected by small $\boldsymbol{\gamma}$ values, and
thus these species are not included in the indirect prediction set.

###  4.2 Inference among species with tied observed counts in Haywood County

| L. Flycatcher | R. Hawk | C. Yellowthroat | E. Kingbird | Bobolink
---|---|---|---|---|---
Haywood | 0.21% | 0.24% | 0.21% | 0.24% | 0.24%
NC-021 | 0.06% | 0.29% | 0.18% | 0.43% | 0.14%
NC-099 | 0.00% | 0.08% | 0.16% | 0.00% | 0.00%
NC-115 | 0.07% | 0.14% | 0.21% | 0.28% | 0.00%
NC-173 | 0.18% | 0.18% | 0.66% | 0.09% | 1.41%
NC-175 | 0.00% | 0.30% | 1.00% | 0.54% | 0.84%
$\boldsymbol{\gamma}$ | 0.4 | 1.29 | 2.43 | 1.49 | 2.15
Table 2: Percentage of all birds observed within each respective county for
species included in either both prediction sets (purple text) or only the
direct set (blue text). Estimated prior hyperparameter $\boldsymbol{\gamma}$
for Haywood County recorded in the last row. Species are sorted by posterior
proportion.

In species abundance data, particularly for areas or counties with small
sample sizes, it is common for multiple species to have the same observed
count. A feature of the construction of the direct order-based prediction
approach as presented is that species with the same observed counts will
either be jointly included or excluded from the prediction set. As a result, a
direct prediction set constructed from a sample with tied species counts may
have increased cardinality over an indirect prediction set that does not
necessarily jointly admit all species with tied observed counts. If the direct
set has increased cardinality for this reason, the direct set will also have
increased coverage over the indirect set.

When constructing a prediction set based on the empirical proportions without
consideration of indirect information, as in the construction of the direct
set, this may commonly occur, and there is no clear approach to choose among
the species with tied counts without further information than what is provided
in the sample in that county. One could randomly choose to include one of the
species from the set of species with tied counts, for example, but a more
principled manner is to utilize indirect information to determine which
species should be included. This is the mechanism used by the indirect
prediction approach when the prior hyperparameter is a real valued vector
estimated from indirect information. As such, a more nuanced benefit of
utilizing indirect information in the construction of a prediction set is the
capacity to include a select few categories with tied empirical proportions.

To demonstrate, we elaborate on species inclusion in the indirect and direct
prediction sets in Haywood County. Haywood is popular destination in the Blue
Ridge Mountains, located near the western border of North Carolina. It
features a moderately large within-county sample size of roughly $4000$ birds
observed. In Haywood County, the indirect prediction set contains $70$
species, and the larger direct set contains $74$. In the construction of these
prediction sets, the ordering of species with regards to the posterior
proportions and the empirical proportions agree for most species. As a result,
all $70$ species included in the indirect set are also included in the direct
set. The disparity in species inclusion occurs primarily as a result of tied
counts of species occurrence in the sample.

Empirical proportions in Haywood and neighboring counties are reported in
Table 2 for the five species included in Haywood County’s prediction sets with
the smallest posterior proportions. The species with the four smallest
posterior proportions are included only in the direct set, and the other
species, the Bobolink, is included in both the indirect and direct sets. The
Bobolink was observed $9$ times in the sample from Haywood County, or about
$0.24\%$ of the Haywood sample. For an ordering determined by either the
empirical counts or the posterior counts, this species is required to be
included in the order-based prediction set to guarantee $1-\alpha$ coverage.
Two of the other species, the Red-shouldered Hawk and Eastern Kingbird, were
each also observed $9$ times in the sample from Haywood, and, by construction
of the order-based prediction approach, must also be included in the direct
set. When admitting the species into a prediction set by posterior counts
based on the real-valued prior hyperparameter $\boldsymbol{\gamma}$ estimated
from data in neighboring counties, as in the indirect approach considered, the
‘tie’ among these three species is broken, and only one, the Bobolink, is
included in the indirect prediction set.

##  5 Discussion

Species abundance data collected across heterogeneous areas is increasingly
important in understanding biodiversity. Some of the largest sources of such
data are citizen science databases for which volunteers spearhead the data
collection. As a result of the civilian-led scientific effort, such data often
feature unequal sampling across a spatial domain where some areas have large
within-area sample sizes and others have much smaller within-area sample
sizes.

In this article, we propose summarizing species abundance data of this type
with valid prediction sets that are constructed by sharing information across
areas. Utilizing indirect information may result in smaller prediction sets
than otherwise achievable with direct methods. Meanwhile, maintaining validity
of the prediction sets for each area allows for an accessible interpretation
that enables a straightforward comparison across areas. In particular,
maintaining interpretable statistical guarantees on a descriptor of such data
is important as analyses from such data often have far reaching policy
implications. Smaller prediction sets may be attainable based on Bayesian
inference of a spatial hierarchical model such as that presented in Tang et
al. (2023), for example, but these approaches introduce bias and a resulting
prediction set would not retain the nominal frequentist coverage rate
guarantee for each county.

The usefulness of our approach for summarizing citizen science data is
motivated in part to combat the common problem of varying sampling efforts
across areas. We detail how $\alpha$-valid prediction sets can be constructed
with the incorporation of indirect information to improve within-county
prediction set precision and propose an empirical Bayes procedure to do so.
Incorporation of accurate indirect information results in a narrower
prediction set for a given county than a direct prediction set by exploiting
data in nearest neighboring counties. The proposed empirical Bayes procedure
is based on a standard hierarchical model that is straightforward to
understand, and the authors provide code for implementation.

There may, however, be a benefit to utilizing a more structured prior that
incorporates indirect information in a more complex manner such as a prior
that weights data from different parts of the state differently. For example,
a model based on a learned intrinsic distance between counties was shown in
Christensen and Hoff (2022) to fit a subset of the eBird data better than
standard methods based on geographic adjacency structure. In the sample
analyzed in Section 4, we found an indirect prediction set constructed with a
hyperparameter estimated from five nearest neighbors results in overall
narrower prediction sets than a direct approach, but it would be valuable to
explore if this can be further improved upon with a more detailed prior. More
broadly, different applications may warrant an alternative information sharing
prior if, for example, there is no notion of spatial distance across the
different areas. For example, it may be of interest to compare species
abundance variation across different time frames for a given county.

All replication codes for this article, including functions to implement the
empirical Bayes estimation procedure for the prior hyperparameter, are
available at https://https://github.com/betsybersson/FreqPredSets_Indirect.

## References

  * Arnold et al. (2021) Arnold ZJ, Wenger SJ, Hall RJ (2021) Not just trash birds: Quantifying avian diversity at landfills using community science data. PLoS ONE 16(9), DOI 10.1371/journal.pone.0255391
  * Bersson and Hoff (2022) Bersson E, Hoff PD (2022) Optimal conformal prediction for small areas. Tech. rep.
  * Camerini and Groppali (2014) Camerini G, Groppali R (2014) Landfill restoration and biodiversity: A case of study in Northern Italy. Waste Management and Research 32(8):782–790, DOI 10.1177/0734242X14545372
  * Christensen and Hoff (2022) Christensen MF, Hoff PD (2022) A nonstationary spatial covariance model for data on graphs. Tech. rep.
  * Faulkenberry (1973) Faulkenberry GD (1973) A method of obtaining prediction intervals. Journal of the American Statistical Association 68(342):433–435
  * Hoff (2023) Hoff P (2023) Bayes-optimal prediction with frequentist coverage control. Bernoulli 29(2):901–928
  * Lele (2020) Lele SR (2020) How should we quantify uncertainty in statistical inference? Frontiers in Ecology and Evolution 8, DOI 10.3389/fevo.2020.00035
  * Rao and Molina (2015) Rao JNK, Molina I (2015) Small Area Estimation, 2nd edn. John Wiley and Sons, Inc., New York, NY
  * Sullivan et al. (2009) Sullivan B, Wood C, Iliff M, Bonney R, FInk D, Kelling S (2009) eBird: a citizen-based bird observation network in the biological sciences. Biological Conservation 142:2282–2292
  * Sullivan et al. (2017) Sullivan BL, Phillips T, Dayer AA, Wood CL, Farnsworth A, Iliff MJ, Davies IJ, Wiggins A, Fink D, Hochachka WM, Rodewald AD, Rosenberg KV, Bonney R, Kelling S (2017) Using open access observational data for conservation action: A case study for birds. Biological Conservation 208:5–14, DOI 10.1016/j.biocon.2016.04.031
  * Tang et al. (2023) Tang B, Clark JS, Marra PP, Gelfand AE (2023) Modeling Community Dynamics Through Environmental Effects, Species Interactions and Movement. Journal of Agricultural, Biological, and Environmental Statistics 28(1):178–195, DOI 10.1007/s13253-022-00520-3
  * Thatcher (1964) Thatcher AR (1964) Relationships between Bayesian and confidence limits for predictions. Journal of the Royal Statistical Society Series B (Methodological) 26(2):176–210
  * Tian et al. (2022) Tian Q, Nordman DJ, Meeker WQ (2022) Methods to compute prediction intervals: A review and new results. Statistical Science 37(4):580–597, DOI 10.1214/21-STS842
  * Vovk et al. (2005) Vovk V, Gammerman A, Shafer G (2005) Algorithmic Learning in a Random World. Springer US

##  Appendix A Maximization of the marginal multinomial-Dirichlet likelihood

In this section, we detail a Newton-Raphson algorithm to maximize the marginal
log likelihood of a conjugate multinomial-Dirichlet model:

$\displaystyle\boldsymbol{X}_{j}\sim{}$ $\displaystyle
MN_{K}({\boldsymbol{\theta}}_{j},N_{j}),\;\text{independently for }j=1,...,J$
$\displaystyle{\boldsymbol{\theta}}_{1},...,{\boldsymbol{\theta}}_{J}\sim{}$
$\displaystyle Dirichlet_{K}(\boldsymbol{\gamma}).$

The log likelihood of the marginal likelihood is as follows,

$\displaystyle\mathcal{L}(\boldsymbol{\gamma}){}$
$\displaystyle\propto\sum_{j=1}^{J}\Bigg{[}\log\Gamma(\sum_{i=1}^{K}\gamma_{i})-\log\Gamma(N_{j}+\sum_{i=1}^{K}\gamma_{i})+$
$\displaystyle\sum_{i=1}^{K}\log\Gamma(x_{j,i}+\gamma_{i})-\sum_{i=1}^{K}\log\Gamma(\gamma_{i})\Bigg{]}.$

Define

$\Psi(s)=\frac{d}{ds}log\Gamma(s)=-\xi+\sum_{n=0}^{\infty}\left[\frac{1}{n+1}-\frac{1}{n+s}\right],$

where $\xi$ is the Euler-Mascheroni constant. Then, it is straightforward to
obtain the first and second derivatives of the marginal log likelihood,

$\displaystyle\frac{d}{d\gamma_{k}}={}$
$\displaystyle\sum_{j=1}^{J}\Bigg{[}\Psi(\sum_{i=1}^{K}\gamma_{i})-\Psi(N_{j}+\sum_{i=1}^{K}\gamma_{k})+\Psi(x_{j,k}+\gamma_{k})-$
$\displaystyle\hskip 42.67912pt\Psi(\gamma_{k})\Bigg{]}$
$\displaystyle\frac{d}{d\gamma_{k}^{2}}={}$
$\displaystyle\sum_{j=1}^{J}\Bigg{[}\Psi^{\prime}(\sum_{i=1}^{K}\gamma_{i})-\Psi^{\prime}(N_{j}+\sum_{i}\gamma_{i})+\Psi^{\prime}(x_{j,k}+\gamma_{k})-$
$\displaystyle\hskip 42.67912pt\Psi^{\prime}(\gamma_{k})\Bigg{]}$
$\displaystyle\frac{d}{d\gamma_{k}d\gamma_{k^{\prime}}}={}$
$\displaystyle\sum_{j=1}^{J}\left[\Psi^{\prime}(\sum_{i=1}^{K}\gamma_{i})-\Psi^{\prime}(N_{j}+\sum_{i=1}^{K}\gamma_{i})\right],$

where $\Psi^{\prime}$ is the trigamma function. Let $\boldsymbol{g}$ be the
gradient vector of length $K$ and $\boldsymbol{H}$ the Hessian matrix.
Finally, Newton’s method updates $\boldsymbol{\gamma}$ as follows:

$\boldsymbol{\gamma}^{(t+1)}=\boldsymbol{\gamma}^{(t)}-\boldsymbol{H}^{-1}(\boldsymbol{\gamma}^{(t)})\boldsymbol{g}(\boldsymbol{\gamma}^{(t)}),$

where the algorithm is iterated until convergence.


article,A,B,C,D,E,F,G,I,J,K,M,N,Q,R,T,U,V,W,X,Y,Z,A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P,Q,R,S,T,U,V,W,X,Y,Z,a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,A,B,C,D,E,F,G,I,K,L,M,N,O,P,Q,R,S,T,U,V,W,X,Y,Z,A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P,Q,R,T,U,V,W,X,Y,Z
, efinition[theorem] and Definition
Conclusion,Conjecture,Corollary,Fact,Facts,Lemma,Observation,Observations,Proposition,Reminder,Scholium,Summary,Theorem,
efinition[theorem] and Definition
Convention,Counterexample,Counterexamples,Discussion,Example,Examples,
Exercise,Exercises,Explanation,Notation,Project,Projects,Question,Questions,Remark,Remarks,Strategy,Warning
[LeSTre2023-]DifferentiallyLarge

# On ordinary differentially large fields

Omar León Sánchez  Omar León Sánchez, The University of Manchester
Department of Mathematics
Oxford Road
Manchester, M13 9PL, UK<EMAIL_ADDRESS>and  Marcus Tressl
Marcus Tressl, The University of Manchester
Department of Mathematics
Oxford Road
Manchester, M13 9PL, UK
Homepage: http://personalpages.manchester.ac.uk/staff/Marcus.Tressl/
<EMAIL_ADDRESS>
###### Key words and phrases:

differential fields, large fields, formal Laurent series

###### 2010 Mathematics Subject Classification:

Primary: 12H05, 12E99. Secondary: 03C60

_Acknowledgments._ The first author was partially supported by EPSRC grant
EP/V03619X/1