# Species Distribution Modeling with Expert Elicitation and Bayesian
Calibration

Karel Kaurila  Sanna Kuningas  Antti Lappalainen  Jarno Vanhatalo  Department
of Mathematics and Statistics, University of Helsinki  Natural Resources
Institute Finland  Organismal and Evolutionary Biology Research Programme,
University of Helsinki

(2022)

###### Abstract

Species distribution models (SDMs) are key tools in ecology, conservation and
management of natural resources. They are commonly trained by scientific
survey data but, since surveys are expensive, there is a need for
complementary sources of information to train them. To this end, several
authors have proposed to use expert elicitation since local citizen and
substance area experts can hold valuable information on species distributions.
Expert knowledge has been incorporated within SDMs, for example, through
informative priors. However, existing approaches pose challenges related to
assessment of the reliability of the experts. Since expert knowledge is
inherently subjective and prone to biases, we should optimally calibrate
experts’ assessments and make inference on their reliability. Moreover,
demonstrated examples of improved species distribution predictions using
expert elicitation compared to using only survey data are few as well. In this
work, we propose a novel approach to use expert knowledge on species
distribution within SDMs and demonstrate that it leads to significantly better
predictions. First, we propose expert elicitation process where experts
summarize their belief on a species occurrence proability with maps. Second,
we collect survey data to calibrate the expert assessments. Third, we propose
a hierarchical Bayesian model that combines the two information sources and
can be used to make predictions over the study area. We apply our methods to
study the distribution of spring spawning pikeperch larvae in a coastal area
of the Gulf of Finland. According to our results, the expert information
significantly improves species distribution predictions compared to
predictions conditioned on survey data only. However, experts’ reliability
also varies considerably, and even generally reliable experts had spatially
structured biases in their assessments. This suggests that expert elicitation
can be an efficient tool, for example, in natural resources management and
conservation area planning, but expert information should be analyzed with
care and preferably calibrated.

expert opinion,

Supra Bayes,

hierarchical models,

Gaussian process,

random effects,

bias correction,

fisheries,

62F15, 62P12;,

60G15,

###### keywords:

###### keywords:

[class=MSC]

††volume: TBA††issue: TBA

, , , and

##  1 Introduction

Species distribution models (SDMs) are key tools in ecology, conservation and
management of natural resources. They are used to study, for example, species
habitat preferences (Kallasvuo et al., 2017; Elith and Leathwich, 2009),
interspecific interactions (Vanhatalo et al., 2020; Ovaskainen and Abrego,
2020) and to build species distribution maps by predicting species presence
and abundance over extended regions where species data have not been collected
(Kotta et al., 2019; Mäkinen and Vanhatalo, 2018; Gelfand et al., 2006). SDMs
are commonly trained by controlled survey data but, since organizing
scientific surveys is expensive and human labor intensive, there is also a
need for complementary sources of information on species distributions. One
such source is expert knowledge which has been used to inform species
distribution predictions either independently or together with survey data
(Murray et al., 2009; Pearman-Gillman et al., 2020; Crawford et al., 2020).
However, expert knowledge is inherently subjective posing challenges related
to model calibration (Murray et al., 2009; Di Febbraro et al., 2018).
Moreover, it is not clear what type of expert knowledge should be collected
and who are the experts to trust (Pearce et al., 2001; Di Febbraro et al.,
2018).

In this work, we propose a novel approach to elicit expert knowledge on
species distribution and to calibrate it using survey data within hierarchical
Bayesian framework. Expert elicitation refers to a process during which a
statistician extracts probability distributions, or point estimates, for a
parameter or a variable of interest from one or more people who are believed
to have valuable information about them. Comprehensive, general, treatments of
modern approaches to expert elicitation are given by O’Hagan et al. (2006) and
Dias et al. (2018). Expert elicitation is extensively used in science,
management and other applications (e.g., Burgman et al., 2011; Vanhatalo et
al., 2014; Nevalainen et al., 2018; O’Hagan, 2019; Perälä et al., 2020; LaMere
et al., 2020). In the context of species distribution modeling, expert
elicitation has been used especially in conservation and management
applications. For example, Crawford et al. (2020) used expert opinion to
inform habitat suitability models for concervation planning in US and Di
Febbraro et al. (2018) assessed the feasibility of monitoring habitat quality
for bird communities in Central Italy by using survey data and expert driven
models. Pearman-Gillman et al. (2020) used expert elicitation to model the
distribution of several species in New England. They used a web questionnaire
to elicit species presence probabilities and information on the impact of
covariates on the species occurency. When fitting the SDM, the assessments
were pooled together by weighting them based on the experts’ self-assessment
on their confidence. Murray et al. (2009) elicited informative priors from
experts for a species distribution model, where the presence or absence of a
threatened species was modelled with logistic regression. They provided
experts with an interactive tool to tune the parameters of the beta
distribution to give their assessment of species presence in various
environments.

We will move forward from the above mentioned approaches by explicitly
modeling and correcting for the systematic biases in the expert assessments.
Correcting for possible biases is crucial for reliable use of expert
information in inference since humans are prone to psychological
idiosyncrasies and subjective biases (Dias et al., 2018; Burgman et al., 2011;
Tversky and Kahneman, 1974). We will follow the so called _supra Bayesian_
approach, where expert information is treated as observations which are linked
to model parameters through a conditional probability distribution (likelihood
function; Genest and Schervish, 1985; Lindley and Singpurwalla, 1986; French,
2011; Albert et al., 2012). Our approach is based on three components: i)
expert elicitation process which summarizes experts’ own assessment of their
region of expertize and belief on species occurrence proability within that
area, ii) data from carefully designed survey within the study region that can
be used to calibrate the expert assessments, and iii) a hierarchical Bayesian
model that combines the two information sources.

We motivate our approach with a real wold case study where we map the
distribution of newly hatched pikeperch (_Sander lucioperca_) larvae within
the Porvoo-Sipoo fisheries region (a public corporation whose purpose is to
develop fishery in their region) in the Gulf of Finland (see Figure 1).
Pikeperch is a top predator of Baltic Sea coastal areas and central species in
the coastal ecosystem of the Baltic Sea. The coastal pikeperch population
forms a commercially important fish stock and pikeperch are also highly sought
after by recreational fishers. Hence, knowledge about its spawning areas
(i.e., areas with newly hatched larvae) has both economical and conservation
importance. Pikeperch is of fresh water origin and has specific habitat
requirements for their reproduction, selecting shallow ($<$10 m deep),
vegetated, and sheltered bays that warm up early in the spring (Veneranta et
al., 2011; Kallasvuo et al., 2017). Information on pikeperch spawning areas is
important for implementation of conservation measures to protect the spawning,
the larvae and male who are securing the development of their offspring. Maps
of the spawning areas can be used to guide local fisheries management and
coastal area planning (Kallasvuo et al., 2017) and our aim is to improve these
maps by combining survey data with expert information.

The rest of the paper is organized as follows. We first summarize the data
collection and expert elicitation process used in our motivating case study.
After that we propose a hierarchical Bayesian model for these data and methods
to conduct model comparison. We then present the results from our case study,
discuss them and provide some concluding remarks.

##  2 Data collection and expert elicitation

Figure 1: The study area in Porvoo-Sipoo archipelago, in the Baltic Sea. The
dots show the locations from where the survey data was collected.

###  2.1 Survey data and environmental covariates

The study area consists of coastal environment types ranging from open water
to sheltered bays. In our analysis, we used three environmental covariates to
characterize the environment: _depth_ (m), _distance to deep_ water ($>$10 m)
and _shoreline density_ (km/m2). The latter two covariates are proxies for how
close to shoreline and how sheltered a location is. Each of the covariates
were available throughout the study area as raster maps with 50 m resolution.
More detailed description of the covariates and how they were constructed is
provided by Kallasvuo et al. (2017).

To collect survey data on the distribution of the larvae of pikeperch from the
study area we conducted a field survey of the surface water layer in June 2017
with paired Gulf ichthyoplankton samplers. These are small nets that are
pulled by a small boat over a transect of 500 meters and have been used to
quantitatively monitor the abundance and spatial occurrence of pikeperch
larvae also in earlier studies (see Veneranta et al., 2011; Kallasvuo et al.,
2017; Långnabba et al., 2019, for a more detailed description). We sampled in
total 92 sites (Figure 1) which were dispersed over the entire study area
covering all main habitat types. The sampling was scheduled to the _a priori_
estimated peak larval season (Liu and Vanhatalo, 2020) in June. From each
sample, we counted early-stage larvae (size range of 3–17 mm) and recorded the
effort; that is, the volume of water sampled, which equals the length of the
transect multiplied by the size (area) of the opening of the Gulf
ichthyoplankton sampler. The Gulf ichthyoplankton samplers are accurate method
for larval sampling and the sampling times were set so that the variability in
catchability due to weather conditions was minimized (Långnabba et al., 2019).
Hence, data collected with the Gulf ichthyoplankton samplers will be treated
as the ground truth in this study.

###  2.2 Expert elicitation

We elicited information concerning the pikeperch spawning grounds from ten
active fishermen (to be called experts hereafter) whose fishing areas fall
inside our study area. The experts were suggested by the executive director of
the local fisheries region. Before the elication, the executive director also
confirmed by a phone call to all experts that the experts were well motivated
to participate in the elicitation process. We sent to each expert by post a
printed map of the study area, crayons with three different colors and filling
instructions. First the experts were asked to draw the borders of their
assessment regions to the map; assessment regions were defined to be areas
within which an expert was confident to state his/her beliefs. After this,
each expert was asked to color the areas, within their assessment region,
where they believed pikeperch did or did not spawn. The experts were asked to
describe the strength of their belief using four categories that were
described as follows (translated from the original Finnish and Swedish
versions):

  * •

(Colour 1) A generally known, locally or regionally important pikeperch
spawning area (the probability that pikeperch spawns in the area is over 90%)

  * •

(Colour 2) Another area where pikeperch most likely spawns (the probability
that pikeperch spawns in the area is 50 – 90 %)

  * •

(Colour 3) An area where pikeperch might spawn (the probability that pikeperch
spawns in the area is 10 – 50%)

  * •

(Uncoloured) The areas that are not marked to belong to any of the above three
categories, but are inside an assessment region, are considered areas where
pikeperch does not spawn according to the general knowledge (the probability
that pikeperch spawns in the area is less than 10%)

Sea areas outside an expert’s assessment regions were considered as missing
information from that expert. Figure A.1 shows an example of an expert
elicitation form and a map drawn by an expert. The elicitation process and
elicitation questions were tested before the actual elicition with the
executive director of the local fisheries region. We then revised the
questions and the elicitation process according to her feedback. The written
descriptions of the colour categories were particularly carefully planned
together so that they would be correctly understood by the fishermen. The
expert elicitation was organized in spring 2018 after which we digitized the
expert drawn maps by scanning them and storing them as raster maps. We aligned
the expert assessment raster maps with the raster maps of the environmental
covariates (Section 2.1) using the ArcGIS software so that each expert’s
answers formed one raster map layer whose lattice match that of the
environmental covariates. The grid cells of the expert assessment maps were
labeled by one of the above four categories or NA (see Figure A.2). The latter
denotes an area outside an expert’s assessment region. There are no
significant errors from the digitalization since the questionnaire maps were
printed in the same coordinate system as the raster maps of the environmental
covariates.

The motivation behind asking directly about spawning areas and the reasoning
for using the four above classes to describe experts’ knowledge were the
following. Earlier works have indicated that expert information is most
accurate within areas from where experts have personal experience (Pearce et
al., 2001; Di Febbraro et al., 2018; Crawford et al., 2020). It has also been
shown that experts are typically better in assessing their belief on real
world variables, which they can observe, than on parameters of statistical
models (O’Hagan et al., 2006). The latter are mathematical abstractions whose
interpretation requires expertise in modeling. Hence, we wanted the experts to
directly assess their beliefs concerning a variable that has a clear
undebatable real world meaning instead of asking them to provide prior
information for some parameters of our species distribution model (Section 3).
The definition of pikeperch spawning area is clear and it is understood in
uniform manner by both the fishermen and researchers. Even though the
fishermen cannot directly observe pikeperch spawning they are assumed to be
skilled to assess spawning areas based on their personal experience on where
pikeperch concentrate during spawning season. We also wanted to allow experts
to express the uncertainty in their knowledge since this provides a more
complete picture of their beliefs than hard cut division to spawning and no-
spawning areas (see also O’Hagan et al., 2006). However, since filling in the
questionnaire maps was rather tedious we wanted to restrict the level of
detail. The four categories described above provide a compromise between these
two targets.

##  3 Species distribution models

###  3.1 Model for survey data

####  3.1.1 Model for larval counts

We followed Liu and Vanhatalo (2020) and Kallasvuo et al. (2017) and modeled
the distribution of pikeperch larvae with a log Gaussian Cox process (LGCP)
with intensity function
$\lambda(\operatorname{\mathbf{s}},\operatorname{\mathbf{x}}(\operatorname{\mathbf{s}}))$
where $\operatorname{\mathbf{s}}$ denotes a location inside the study area and
$\operatorname{\mathbf{x}}(\operatorname{\mathbf{s}})$ is a vector of
environmental covariates at that location. We modeled the log intensity with a
linear function of the covariates and a spatial random effect:

$\log\lambda(\textbf{s}_{i},\textbf{x}_{i})=\alpha+\boldsymbol{\beta}^{T}\textbf{x}_{i}+\phi(\textbf{s}_{i}),$
(3.1)

where the intercept $\alpha$ and linear weights, $\boldsymbol{\beta}$ were
given a zero mean Gaussian prior with variance 100. The spatial random effect
$\phi(\textbf{s})$ followed the Barrier model of Bakka et al. (2019), which is
a non-stationary Gaussian process whose covariance does not travel through
physical barriers – islands in our application.

The Barrier model is defined as a Stochastic Partial Differential Equation
(SPDE) which leads to a Gaussian process with a non-stationary covariance
function. For a stationary covariance function the covariance between two
points $\operatorname{\mathbf{s}}$ and $\operatorname{\mathbf{s}}^{\prime}$
depends only on the distance
$d(\operatorname{\mathbf{s}},\operatorname{\mathbf{s}}^{\prime})$. In the
Barrier model the covariance travels only through water, while land areas are
considered as barriers through which the covariance does not travel. This
property is achieved by defining the spatial random effect, $\phi(s)$, as the
solution to the following SPDE:

$\displaystyle\phi(s)-\nabla\frac{r^{2}}{8}\phi(s)$
$\displaystyle=r\sqrt{\frac{\pi}{2}}\sigma_{\phi}\mathcal{W}(s),\text{ for
}s\in\Omega_{w}$ $\displaystyle\phi(s)-\nabla\frac{r^{2}_{b}}{8}\phi(s)$
$\displaystyle=r_{b}\sqrt{\frac{\pi}{2}}\sigma_{\phi}\mathcal{W}(s),\text{ for
}s\in\Omega_{l},$ (3.2)

where $\Omega_{w}$ is water, $\Omega_{l}$ is land, $\mathcal{W}(s)$ is white
noise, $\sigma_{\phi}$ is the marginal standard deviation of the process and
$r$ is the range parameter governing how fast the correlation decreases with
distance through water. The parameter $r_{b}$ is the correlation range on land
and it is set to be a fraction of $r$ to remove the correlation there (Bakka
et al., 2019). Here, we use the default value $r_{b}=0.2r$. The details on
implementing the barrier model are given in Appendix B.1. The parameters of
the barrier model follow a Penalized Complexity (PC) prior (Simpson et al.,
2017; Bakka et al., 2019). These priors shrink the effect towards a base
model, where $\sigma\rightarrow 0$ and $l\rightarrow\infty$. We defined the
prior through probabilities $p(\sigma>1)=0.01$ and $p(l<0.5km)=0.01$.

As detailed in Section 2.1 the transect observations correspond to the number
of pikeperch larvae caught with a net that samples a volume of water over a
transect. We denote by $\operatorname{\mathbf{s}}_{i}$ the coordinates of the
middle point of the $i$th transect, by $y_{i}$ the number of larvae caught
over the transect and by
$\operatorname{\mathbf{x}}_{i}=\operatorname{\mathbf{x}}(\operatorname{\mathbf{s}}_{i})$
the environmental covariates at the middle of the transect. Since survey
transects were short compared to the resolution of the mesh used to implement
the model (see Appendix B.1), we treated the larvae density as a constant over
each transect. Hence, the observation model for the survey data
$\operatorname{\mathbf{y}}=[y_{1},\dots,y_{n}]^{T}$ is
$p(\operatorname{\mathbf{y}}|\boldsymbol{\lambda},\boldsymbol{\epsilon})=\prod_{i=1}^{n}\mathrm{Poisson}(V_{i}\lambda_{i}\epsilon_{i})$
where $V_{i}$ is the volume of water sampled,
$\lambda_{i}=\lambda(\operatorname{\mathbf{s}}_{i},\operatorname{\mathbf{x}}_{i})$
is the intensity of the log Gaussian Cox process and
$\boldsymbol{\epsilon}=[\epsilon_{1},\dots,\epsilon_{n}]^{T}$ is a vector of
independent random effects capturing overdispersion in larval counts arising
from local phenomena and stochasticity in sampling. Note though that, in the
survey data set, each point has the same sampling volume $V_{i}=V$ (28.35 m3)
so we gave Gamma prior for the random effects,
$\epsilon_{i}\sim\mathrm{Gamma}(r,1/r)$, which leads to negative binomial
observation model for the larval counts

$p(\operatorname{\mathbf{y}}|\lambda,r)=\prod_{i=1}^{n}\mathrm{Negative-
Binomial}(V_{i}\lambda_{i},r)$ (3.3)

so that the observed larval counts are actually modeled to arise from an
overdispersed log Gaussian Cox process which approaches the Poisson process as
the overdispersion parameter $r$ increases (see also, Liu and Vanhatalo,
2020). We gave $\text{Gamma}(\sqrt{10},1/\sqrt{10})$ prior for the
overdispersion parameter $r$.

####  3.1.2 Model for larval presence

As an alternative to larval abundance modeling we considered also an
occurrence model where the occurrence of larvae in a survey sample was modeled
with

$\mathbbm{1}(y_{i}>0)\sim\text{Bernoulli}\left(\pi(\textbf{s}_{i},\textbf{x}_{i})\right)$
(3.4)

where $\mathbbm{1}(y_{i}>0)$ is an indicator function returning one if
$y_{i}>0$ and zero otherwise and
$\pi(\textbf{s}_{i},\textbf{x}_{i})=\text{logit}^{-1}\left(\alpha+\boldsymbol{\beta}^{T}\textbf{x}_{i}+\phi(\textbf{s}_{i})\right)$
denotes the probability of presence for pikeperch larvae. The priors for model
parameters $\alpha$, $\beta$ and $\phi$ were the same as in the abundance
model.

###  3.2 Joint model for expert assessments and survey data

####  3.2.1 Model for expert opinions

We constructed the model for expert assessments by first building a model for
experts’ subjective probabilities on the occurrence of larvae. For this, we
denote by $\bar{\pi}_{ji}\in[0,1]$ the $j$th expert’s subjective probability
that larvae are present at location $s_{i}$ and model this probability with
Beta distribution

$\bar{\pi}_{ji}\sim Beta(\bar{\mu}_{ji},\bar{s}_{j}),$ (3.5)

where $\bar{\mu}_{ji}=E[\bar{\pi}_{ji}]$, and $\bar{s}_{j}$ is the (prior
sample size) parameter governing the spread of the Beta distribution. We fixed
$\bar{s}_{j}=2$ to encode vague prior predictive distribution for expert
opinions. We then model the expected value of an expert’s subjective
probability with logistic regression so that it is a function of the log of
the true larval density or the logit of the true probability of presence of
larvae, so that

$\bar{\mu}_{ji}=\mathrm{logit}^{-1}(\bar{\alpha}_{j}+\bar{c}_{j}(\boldsymbol{\beta}^{T}\textbf{x}_{i}+\phi(\textbf{s}_{i}))+\bar{\varphi}_{j}(\textbf{s}_{i})).$
(3.6)

Here, the parameter $\bar{\alpha}_{j}$ is an intercept, $\bar{c}_{j}$ is a
parameter for the expert’s skill and $\bar{\varphi}_{j}(\textbf{s}_{i})$ is a
residual error. Note that, for example,
$\log\lambda_{i}-\alpha=\boldsymbol{\beta}^{T}\textbf{x}_{i}+\phi(\textbf{s}_{i})$,
so $\text{logit}(\mu_{ji})$ is proportional to $\log\lambda_{i}$ (or
$\mathrm{logit}(\pi_{i})$) but we have subtracted $\alpha$ from (3.6) in order
to improve the identifiability of the parameters. The parameter $\bar{c}_{j}$,
thus, describes how strongly an expert’s assessment follows the true
underlying larvae intensity or probability of presence, so that a positive
$\bar{c}_{j}$ indicates that an expert has information about larvae
distribution. The error terms $\bar{\varphi}_{j}(\textbf{s}_{i})$ correct for
spatially correlated local biases in the expert assessment that can result
from actual bias in an expert’s opinion but also from inaccuracies in
expressing them. For example, it is likely that the experts coloured the maps
with smaller resolution than the true resolution in the larvae
presence/absence pattern, so $\bar{\varphi}_{j}(\textbf{s}_{i})$ also accounts
for spatial correlation resulting from this.

In case of a conflict between expert assessment and survey data, we want the
model to explain the expert assessment with $\bar{\varphi}$ instead of
$\bar{c}_{j}(\boldsymbol{\beta}^{T}\textbf{x}_{i}+\phi(\textbf{s}_{i}))$ so
that an expert’s assessment does not provide information on larvae intensity
or occurrence probability. To attain this behaviour we give relatively
stricter prior for $\bar{c}_{j}$ than to $\bar{\varphi}$. Hence, priors for
the parameters were $\bar{\alpha}_{j}\sim N(0,2^{2})$ and $\bar{c}_{j}\sim
N(0,0.5^{2})$. The latter implies 95% prior probability that $\bar{c}_{j}$ is
less than one. We then modeled the expert bias function
$\bar{\varphi}_{j}(\textbf{s}_{i})$, with Besag-York-Mollié (BYM) -model
(Besag et al., 1991), which induces spatial correlation among neighboring
pixels in an expert assessment graph (see Appendix B.1 for details). In the
BYM-model $\bar{\varphi}_{j}\sim N(0,Q^{-1})$, where the precision matrix
$Q=\tau_{u}R+\tau_{v}I$, $I$ is the identity matrix corresponding to the
i.i.d. random effect, and

$R_{kl}=\begin{cases}n_{k},&k=l\\\ -\mathbbm{1}\\{k\sim l\\},&k\neq
l\end{cases}$ (3.7)

where $n_{k}$ is the number of neighbors of vertex $k$, and $k\sim l$
indicates that vertices $k$ and $l$ are neighbors, i.e. they are connected by
an edge. The hyperparameters of the BYM model are precision for the spatially
correlated effect $\tau_{u}$ and precision for the i.i.d. effect $\tau_{v}$.
We gave both of them gamma prior $\Gamma(2,8)$. When marginalizing over the
precision this induces a heavy tailed scaled Student-$t$ distribution with
$\nu=4$ and $s=2$ for $\bar{\varphi}$.

####  3.2.2 Likelihood functions for the expert assessments

As explained in Section 2.2 the expert assessments were coded as raster maps
so that the assessment of the $j$’th expert for grid cell $i$ at location
$\operatorname{\mathbf{s}}_{i}$ is a categorical variable
$z_{ij}\in\\{1,2,3,4\\}$ if the grid cell belongs to the expert’s assessment
region (Figure A.2). If the grid cell is outside expert’s assessment area
$z_{ij}=\text{NA}$, which means that the corresponding grid cell is excluded
from the likelihood function of the expert’s assessments to be detailed below.
We ordered the categories so that $z_{ij}=1$ corresponds to the lowest
($<10\%$) and $z_{ij}=4$ corresponds to the highest ($>90\%$) subjective
probability of presence of an expert.

As the first model for the expert assessments we treated them as binary
statements about presence of pikeperch larvae. In this model, a location
within an expert’s assessment region was considered to be labeled as presence
if the expert had given it over 50% probability for larvae presence (i.e.
$z_{ij}\in\\{3,4\\})$) and otherwise it was considered as an absence. The
observation model for the binary absence assessment is then

$\displaystyle\text{Pr}(z_{ij}\in\\{1,2\\})$
$\displaystyle=\text{Pr}(\bar{\pi}_{ji}\leq
0.5)=F_{\text{Beta}}(0.5|\bar{\mu}_{ji},\bar{s}_{j})$ (3.8)

where $F_{\text{Beta}}(0.5|\mu_{ji},s_{j})$ is the cumulative distribution
function of the Beta distribution. Similarly for the binary presence
assessment we have
$\text{Pr}(z_{ij}\in\\{3,4\\})=1-F_{\text{Beta}}(0.5|\bar{\mu}_{ji},\bar{s}_{j})$
This observation model induces a likelihood function for
$\text{logit}(\bar{\mu}_{ji})$ that resembles closely the logit-function. We
can similarly form the observation model for all expert assessment categories

$\displaystyle\text{Pr}(z_{ij}|\mu_{ji},s_{j})$
$\displaystyle=\begin{cases}F_{\text{Beta}}(0.1|\mu_{ji},s_{j}),&\text{ if
}z_{ij}=1\\\
F_{\text{Beta}}(0.5|\mu_{ji},s_{j})-F_{\text{Beta}}(0.1|\mu_{ji},s_{j}),&\text{
if }z_{ij}=2\\\
F_{\text{Beta}}(0.9|\mu_{ji},s_{j})-F_{\text{Beta}}(0.5|\mu_{ji},s_{j}),&\text{
if }z_{ij}=3\\\ 1-F_{\text{Beta}}(0.9|\mu_{ji},s_{j}),&\text{ if
}z_{ij}=4,\end{cases}$ (3.9)

which again induces likelihood functions for $\text{logit}(\mu_{ji})$ and
through that to the model parameters.

We made a further assumption that, conditionally on the model parameters, the
expert assessments are mutually independent. Hence, the joint distribution for
all expert assessments over the whole study region can be written as

$p\left(\operatorname{\mathbf{z}}_{1},\dots,\operatorname{\mathbf{z}}_{J}|\bar{\alpha},\bar{c},\boldsymbol{\beta},\varphi,\bar{\varphi}\right)=\prod_{j=1}^{J}\prod_{i=1}^{n_{j}}p(z_{ij}|\bar{\alpha}_{j},\bar{c}_{j},\boldsymbol{\beta},\varphi,\bar{\varphi}_{j})$
(3.10)

where $J$ is the number of experts, $n_{j}$ is the number of grid cells inside
the assessment area of the $j$th expert,
$\bar{\alpha}=\\{\bar{\alpha}_{1},\dots,\bar{\alpha}_{J}\\}$,
$\bar{c}=\\{\bar{c}_{1},\dots,\bar{c}_{J}\\}$ and
$\bar{\varphi}=\\{\bar{\varphi}_{1},\dots,\bar{\varphi}_{J}\\}$. When we
combine the two data sets, $\operatorname{\mathbf{y}}$ (the survey
observations) and $\operatorname{\mathbf{z}}_{j},j=1,\dots,J$ (the expert
observations), the observations are independent given the latent function and
model parameters so that, in case of count observations from surveys, the full
observation model is

$p(\operatorname{\mathbf{y}},\operatorname{\mathbf{z}}_{1},\dots,\operatorname{\mathbf{z}}_{J}|\bar{\alpha},\bar{c},\bar{\varphi},\alpha,\boldsymbol{\beta},\varphi,r)=\left[\prod_{i=1}^{n}p(y_{i}|\lambda_{i},r)\right]\prod_{j=1}^{J}\prod_{k=1}^{n_{j}}p(z_{kj}|\bar{\alpha}_{j},\bar{c}_{j},\bar{\varphi}_{j},\boldsymbol{\beta},\varphi).$
(3.11)

The full model with the presence-absence model for the survey observations is
constructed analogously.

##  4 Posterior inference and model comparison

We implemented all the models and conducted the posterior inference using the
Integrated Nested Laplace Approximation (INLA) R package (R-INLA, Rue et al.,
2009). The technical implementation of the barrier model requires a triangular
mesh over the study area on which the model is constructed (Bakka et al.,
2019). The survey observations and expert assessments are then mapped to the
mesh with projection matrices. The details on how we constructed these are
given in the Appendix B.1. R-INLA does not support the exact likelihood
functions derived in Section 3.2, so we searched for the closest match to them
from among the Binomial likelihood functions, which are readily available in
R-INLA. The details of the matching process and the resulting likelihood
functions are given in the Appendix C. After setting up the models, we
constructed the INLA approximation for the model parameters and used the model
to predict the larvae intensity over the whole study area.

Since the survey data were collected with a standardized and extensively
tested sampling method (see Veneranta et al., 2011; Kallasvuo et al., 2017,
and their references) and the sampling time was matched with the peak larval
time (Liu and Vanhatalo, 2020), we considered them to accurately reflect the
larval areas of pikeperch. On the other hand, we had less evidence on the
expert elicitation process and the reliability of the experts. Hence, we
considered the survey data to be more reliable source of information on larvae
distribution than the expert assessments, and compared the alternative models’
based on how well they predict survey data that have been left out from the
training data. The INLA approximation gives access to approximate leave-one-
out (LOO) predictive densities for all observations. These are also called
Conditional Predictive Ordinate (CPO) values and, for the $i$’th survey
observation, it is defined as

$\mathrm{CPO}_{i}=p(y_{i}|\mathbf{x}_{i},\mathbf{s}_{i},D_{\setminus i})$
(4.1)

where $D_{\setminus i}$ includes all the other data except
$\\{y_{i},\mathbf{x}_{i},\mathbf{s}_{i}\\}$. The CPO values can be used to
calculate the widely used average LOO cross-validation log predictive density

$\mathrm{lpd}=\frac{1}{n_{\mathrm{survey}}}\sum_{i=1}^{n_{\mathrm{survey}}}\log\mathrm{CPO}_{i}.$
(4.2)

With models where survey data are modeled as presence-absence data, we can
calculate three additional model comparison metrics. First, we calculated
classification accuracy,
$\mathrm{ACC}=\frac{1}{n_{\mathrm{survey}}}\sum_{i=1}^{n_{\mathrm{survey}}}\mathbbm{1}(\mathrm{CPO}_{i}\geq
0.5)$, i.e. the ratio of correctly classified observations. Because the
observations are skewed towards absences, we also calculated the balanced
classification accuracy $\mathrm{bACC}=\frac{\mathrm{TPR}+\mathrm{TNR}}{2}$,
which is the average of true positive rate (TPR) and true negative rate (TNR)
that are calculated similarly to ACC. Log predictive density is known to be
sensitive to outliers (i.e., observations for which the predictive density is
much lower than in average) for which reason Gneiting et al. (2007) proposed
to use the Continuous Ranked Probability Score (CRPS) as an alternative metric
for comparing probabilistic predictions. The CRPS for a single prediction is
defined as

$\text{CRPS}(F_{i},y_{i})=E_{F_{i}}|Y_{i}-y_{i}|-\frac{1}{2}E_{F_{i}}|Y_{i}-Y_{i}^{\prime}|,$
(4.3)

where $F_{i}$ is the posterior predictive distribution given by the model
being evaluated, $Y_{i}$ and $Y_{i}^{\prime}$ are independent and identical
random variables that follow $F_{i}$ and $y_{i}$ is an observation. In the
case of Bernoulli likelihood for survey data, the LOO posterior predictive
distribution for the $i$’th observation is $Y_{i}\sim
F=\text{Bernoulli}(\hat{\pi}_{i})$ where
$\hat{\pi}_{i}=\text{Pr}(Y_{i}=1|\mathbf{x}_{i},\mathbf{s}_{i},D_{\setminus
i})$. Hence, we can simplify the CRPS for the LOO posterior predictive
distributions of these models to the square of the posterior predictive
probability of the incorrect class:

$\displaystyle\text{CRPS}(F_{i},y_{i})$
$\displaystyle=\hat{\pi}_{i}(1-y_{i})+\left(1-\hat{\pi}_{i}\right)y_{i}-\frac{1}{2}\left(2\hat{\pi}_{i}\left(1-\hat{\pi}_{i}\right)\right)$
$\displaystyle=\begin{cases}\hat{\pi}_{i}^{2}&\text{ if }y_{i}=0\\\
(1-\hat{\pi}_{i})^{2}&\text{ if }y_{i}=1\end{cases}$
$\displaystyle=(1-\mathrm{CPO}_{i})^{2}.$

##  5 Results

We received in total 11 expert assessments out of which ten were used for the
analyses. The assessed areal coverage of experts varied considerably. Some
experts expressed their views for very small regions whereas some of them
covered almost the entire study area (see Figure A.2). One expert assessment
was excluded, as the area assessed by the expert was considered too small for
the analysis. In some cases, the request of first drawing the expert’s own
assessment area in the map was poorly understood and needed to be clarified
later on. Some experts were more uncertain in their statements (they used only
categories 2 or 3) than others (who used all four categories). In general the
digitalization of the expert assessment maps worked reasonably well but in few
cases borders between two colors were hard to distinguish and needed some
interpretation by us.

In general, the models integrating survey data and expert assessments
(hereafter survey+expert models) outperformed the survey only model in their
predictive performance (table 1). In occurrence predictions, that is survey
data were modeled as occurrence data (Bernoulli likelihood (3.4)), the best
model was the one where expert assessments were modeled as presence-absence
statements (3.8). However, the survey+expert models outperformed the survey
only model in terms of classification accuracy (ACC and bACC) and CRPS
statistics but not in LPD statistics (see table 1). Since LPD is more
sensitive to outlying observations than CRPS (Gneiting et al., 2007), this
indicates that survey+expert models made more over-confident false presence-
absence predictions than the survey only models. In count predictions, that is
survey data were modeled as count data (Negative-Binomial likelihood (3.3)),
the survey+expert models had clearly better LPD statistics than the survey
only model (table 1). The differences in LPD statistics between alternative
larval abundance models were rather large and the best model was the
survey+expert model where we used all four expert assessment categories,
observation model (3.9).

Table 1: Posterior predictive model comparison for alternative models. First two columns specify the observation model for survey data and expert elicited maps; ”–” in the latter denotes a survey data only model, ”p/a” denotes the presence absence model ((3.4) for survey data and (3.8) for expert assessments), and ”abu” denotes abundance model (3.3) for survey data and four categories model (3.9) for expert assessments. Other abbreviations are lpd for leave-one out cross validation log predictive density, ACC for classification accuracy, bACC for balanced classification accuracy and CRPS for conditional predictive ordinate. Observation models | lpd | ACC | bACC | CRPS
---|---|---|---|---
Survey data | Expert maps |  |  |  |
p/a | – | -83.6 | 0.718 | 0.683 | 0.538
p/a | p/a | -95.2 | 0.737 | 0.724 | 0.436
p/a | abu | -87.3 | 0.737 | 0.721 | 0.465
abu | – | -452 | - | - | -
abu | p/a | -410 | - | - | -
abu | abu | -294 |  |  |

All the models produced posterior predictive larval distribution maps that
were similar in their overall pattern; that is, the larvae density and
occurrence probability is the highest in the sheltered bays in the northern
parts of the study area and smallest in open sea areas in the south (see
figures 2 and B.4a). However, there are clear differences in the finer scale
patterns of the predictions. For example, the survey+expert models predict
consistently lower densities and occurrence probabilities in southern areas
than the survey only models. The survey+expert models predicted also somewhat
higher larval densities and occurrence probability than the survey only model
in some of the bays in the north whereas in the north-easternmost bay the
survey only models tended to predict higher densities and occurrence
probabilities than the survey+expert model. In general, the survey+expert
models reduce the uncertainty compared to survey only models almost everywhere
in the study area. Only in the south-eastern and south-western areas, the
survey+expert models had larger posterior predictive standard deviation for
the log intensity and logit probability than the survey only models (see
figures 3 and B.4b).

Figure 2: Posterior predictions for larval distribution in the study area with
survey only and survey+experts models and their difference ([survey only] -
[survey+experts]). On top row, the survey data are modeled as occurrence data
and the maps show the posterior predictive mean of the log odds ratio for
larval occurrence and their difference. On bottom row, survey data are modeled
as count data and the maps show the posterior predictive mean of the log
larval density and their difference. In both rows, the expert assessments
included all four categories.

Figure 3: Posterior predictive uncertainty for larval distribution in the
study area with survey only and survey+experts models and their difference
([survey only] - [survey+experts]). On top row, the survey data are modeled as
occurrence data and the maps show the posterior predictive standard deviation
of the log odds ratio for larval occurrence and their difference. On bottom
row, survey data are modeled as count data and the maps show the posterior
predictive standard deviation of the log larval density and their difference.
In both rows, the expert assessments included all four categories.

Experts’ reliability, when measured by the $\bar{c}$ parameters, varied
considerably (see figures 4(b) and B.2). Experts 2-5 were the reliable ones
having significantly positive $\bar{c}$ parameters whereas all the other
experts were unreliable in the sense that the posterior distributions their
$\bar{c}$ parameters did not differ from zero significantly. The third expert
showed some local bias in his/her assessment in a couple of inner bayes in the
north but otherwise the spatial bias terms of the reliable experts were small
(see Figure B.3). The posterior distributions of the parameters of the survey
data model were considerably narrower in the survey+expert models than in the
survey only models (figures 4(a) and B.1). Log larval density and log odds
ratio of larval occurrence responded negatively to depth and openness
(lined3km) and positively to distance to 10m depth curve. The overdispersion
in survey data was significant since the posterior distribution of
overdispersion parameter, $r$, was concentrated in small values. The barrier
model had more significant effect on log larval density and log odds-ratio of
occurrence in survey+expert models than in the survey only models

(a) Survey model parameters

(b) Expert coefficients

Figure 4: Posterior distributions for model parameters in the survey+experts
model with Negative-Binomial observation model for survey data and Beta
observation model for expert assessments.

We also conducted sensitivity analyses for the prior distributions. With too
relaxed barrier model priors, INLA sometimes ran into numerical problems, with
the latent variable $\eta$ evaluating to `NAN`\- or `INF`-values at several
points in the mesh. Numerical problems did not arise when we defined the PC-
priors to correspond the general knowledge on likely scale and spatial
correlation distance in log density and log odds ratio (the priors reported in
Section 2.1). When the BYM-model was left out, or when we restricted its
variance parameters too heavily towards zero by prior distributions, the
survey+expert models effectively ignored the survey data information and
fitted the log abundance and log odds ratio based on experts information only.
However, when expert bias terms (BYM-model components) were _a priori_ allowed
to vary more than the log density or log odds ratio, conflict between survey
data and expert assessment was explained by the bias term, and the expert
assessment did not contribute to the posterior distribution of log density and
log odds ratio in practice. In addition to prior distributions, the resolution
of expert-mesh had an impact on how much weight we gave to expert assessments
so that larger grid cells in the expert-mesh implied reduced the effect of
expert assessments compared to fine grid cells.

##  6 Discussion

Our results show clear differences between the survey only and survey+expert
models both in predictions and parameter inference. In general, the
survey+expert models outperformed the survey only models in their predictive
performance, indicating that adding expert information to SDMs can be
benefitial. The general pattern of the predicted species distributions was
similar across different models, but there were clear differences in some
areas where expert assessments were available (Figure 2). The biggest
differences between the alternative models being in the south. This is
reasonable since all experts, who had marked the southern areas into their
assessment regions, had given the smallest probability for larval occurrence
there (see Figure A.2) whereas there were only a few survey observations from
the southern part of the study region. Hence, the expert assessments have the
largest impact in the south. In addition to changing the predictive mean of
the larvae abundance and occurrence probability, including expert assessments
into the models decreased their predictive uncertainty as well. Similarly, the
posterior distributions of the parameters of the larval density and occurrence
probability models were considerably narrower in the survey+expert models than
in the survey only models (figures 4(a) and B.1). All these results show that
the expert assessments provided significant information on factors affecting
larvae distribution.

These results suggest that expert elicitation can be a viable solution in
studies concerning distribution of species in space and time when resources
for scientific sampling are limited. However, our results demonstrate also
that expert elicitation should be applied with care and optimally combined
with more objective information to reduce the effect of subjective biases.
These findings align well with earlier studies that have shown that the skill
of experts in assessing species distributions varies and depends, for example,
on study area, species and background of an expert (Pearce et al., 2001; Di
Febbraro et al., 2018; Crawford et al., 2020). Bias correction and weighing
expert assessments based on experts’ skills have been emphasized also in many
other expert assessment tasks (O’Hagan et al., 2006; Burgman et al., 2011;
Dias et al., 2018; Perälä et al., 2020).

The expert elicitation process turned out well-functioning. Important first
step was to confirm the commitment of the experts to the process. As a result,
responses were gotten from all candidate experts. For successful
implementation of an expert elicitation by post, it is important that
instructions for filling questionnaire are clear and that the level of details
are kept meaningful. Another, but more consuming option, would be to run a
face-to-face meeting or an expert elicitation workshop to collect the expert’s
views.

The technical implementation of our models was done using the R-INLA software
(Rue et al., 2009) since it allows straightforward implementation of the
barrier model (3.1.1). The barrier model is justified in our application area
since it is scattered by islands, peninsulas and other physical obstacles for
marine species that make traditional stationary spatial random effect an
unrealistic assumption (Bakka et al., 2019). This choice did not come without
price, though, since we were restricted to the built-in likelihood functions
in INLA which affected our practical choice for the expert assessment models
(see Section 3.2 and Appendix C). We believe this choice had only minor effect
on the results though. More fundamental technical consideration is related to
the prior distributions. Firstly, the priors for the hyperparameters of
barrier model had to be chosen carefully to keep the R-INLA calculations
numerically stable. Secondly, to prevent expert assessments from overruling
the survey observations, the mesh resolution and the prior distributions for
the hyperparameters of the expert bias terms (BYM-model parameters) had to be
chosen reasonably. Too small mesh size and too narrow prior for the BYM-model
prevents the model from rejecting biased expert assessments. This is
understandable because we included the expert assessments into our models as
point-wise observations at expert-mesh nodes. Hence, if expert bias terms
(BYM-model components) were restricted to near zero, the combined likelihood
of expert assessments would outweight the likelihood of survey data in
information concerning the log larvae density or log odds ratio of larvae
occurrence probability (there are far more expert-mesh nodes than survey data
observation locations). For this reason, the expert-mesh resolution should be
chosen so that it corresponds to the resolution at which the experts can draw
their maps. Similarly, the variance parameters of the BYM-model should be
given wide priors.

In summary, our results are encouraging for further development of expert
elicitation methods in species distribution modeling. They suggest also that
the methods proposed in this work could also be scaled to larger scale
applications; such as environmental accounting and environmenal management.
Biodiversity loss has become an equally important challenge for humanity and
wellbeing of our planet as climate change. Today, confrontation between nature
and economic wellbeing has started to disappear and general demand for
introducing natural capital into national accounting systems to appear
(Dasgupta, 2021). This demand parallels other calls for biodiversity
preserving measures (e.g., the System of Environmental Economic Accounting
United Nations et al., 2021) in that they all rely on monitoring of species,
biodiversity and ecosystem processes. Current technology does not, however,
allow continuous measuring of biota over space and time so predictive models
are needed to fill in the gaps. Species distribution models are routinely used
for this task so that they are trained with observational data to make
predictions on species occurrence and abundance at spatiotemporal locations
not covered by observations. However, survey data can be expensive and
logistically challenging to collect whereas expert information can often be
collected with relatively inexpensively. For this reason expert elicitation is
a tempting method to collect compelementary information to survey data
(Pearman-Gillman et al., 2020; Murray et al., 2009). The other benefit from
expert elicitation in the context of environmental accounting could be
stakeholder engagement since management and policy decisions are typically
better received by stakeholders and other interest groups if they have been
heard and their knowledge incorporated in the decision making process (La Mere
et al., 2020).

##  7 Conclusions

We have presented a formal Bayesian approach to include expert information
into species distribution modeling. Our approach effectively integrates
spatial expert information with point-wise survey data and allows expert
calibration and assessment of experts’ reliability. Our results show that
expert information can significantly improve species distribution predictions
compared to predictions conditioned on survey data only. This suggests that
expert elicitation can be an efficient tool, for example, in natural resources
management, conservation area planning and in other applications that require
knowledge on distributions of species and other biotic variables. However, our
results also demonstrate that experts’ reliability may considerably vary, and
that even generally reliable experts can have spatially structured biases in
their beliefs. Hence, expert information should be analyzed with care and
preferably calibrated with more objective information sources such as survey
data.

## Acknowledgements

This work has been funded by the Academy of Finland (grant 317255, KK and JV)
and the Finnish Operational Program for the European Maritime and Fisheries
Fund (EMFF, AL and SK).

## References

  * Albert et al. (2012) Albert, I., Donnet, S., Guihenneuc-Jouyaux, C., Low-Choy, S., Mengersen, K., and Rousseau, J. (2012).  “Combining Expert Opinions in Prior Elicitation.”  Bayesian Analysis, 7(3): 503–532.
  * Bakka et al. (2019) Bakka, H., Vanhatalo, J., Illian, J. B., Simpson, D., and Rue, H. (2019).  “Non-stationary Gaussian models with physical barriers.”  Spatial Statistics, 29: 268 – 288.
  * Besag et al. (1991) Besag, J., York, J., and Mollié, A. (1991).  “Bayesian image restoration, with two applications in spatial statistics.”  Annals of the institute of statistical mathematics, 43(1): 1–20.
  * Burgman et al. (2011) Burgman, M., Carr, A., Godden, L., Gregory, R., McBride, M., Flander, L., and Maguire, L. (2011).  “Redefining expertise and improving ecological judgment.”  Conservation Letters, 4(2): 81–87.
  * Crawford et al. (2020) Crawford, B. A., Maerz, J. C., and Moore, C. T. (2020).  “Expert-Informed Habitat Suitability Analysis for At-Risk Species Assessment and Conservation Planning.”  Journal of Fish and Wildlife Management, 11(1): 130–150.
  * Dasgupta (2021) Dasgupta, P. (2021).  “The Economics of Biodiversity: The Dasgupta Review.”  Technical report, London: HM Treasury.
  * Di Febbraro et al. (2018) Di Febbraro, M., Sallustio, L., Vizzarri, M., De Rosa, D., De Lisio, L., Loy, A., Eichelberger, B., and Marchetti, M. (2018).  “Expert-based and correlative models to map habitat quality: Which gives better support to conservation planning?”  Global Ecology and Conservation, 16: e00513.
  * Dias et al. (2018) Dias, L. C., Morton, A., and Quigley, J. (2018).  Elicitation.  Springer International Publishing.
  * Elith and Leathwich (2009) Elith, J. and Leathwich, J. R. (2009).  “Species Distributions Models: Ecological Explanation and Predictions Across Space and Time.”  The Annual Review of Ecology, Evolution and Systematics, 40(677-697).
  * French (2011) French, S. (2011).  “Aggregating expert judgement.”  Revista de la Real Academia de Ciencias Exactas, Fisicas y Naturales. Serie A. Matematicas, 105(1): 181–206.
  * Gelfand et al. (2006) Gelfand, A. E., Jr, J. A. S., Wu, S., Latimer, A., Lewis, P. O., Rebelo, A. G., and Holder, M. (2006).  “Explaining Species Distribution Patterns through Hierarchical Modelling.”  Bayesian Analysis, 1(1): 41–92.
  * Genest and Schervish (1985) Genest, C. and Schervish, M. J. (1985).  “Modeling Expert Judgements for Bayesian Updating.”  The Annals of Statistics, 13(3): 1198–1212.
  * Gneiting et al. (2007) Gneiting, T., Balabdaoui, F., and Raftery, A. E. (2007).  “Probabilistic forecasts, calibration and sharpness.”  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69(2): 243–268.
  * Kallasvuo et al. (2017) Kallasvuo, M., Vanhatalo, J., and Veneranta, L. (2017).  “Modeling the spatial distribution of larval fish abundance provides essential information for management.”  Canadian Journal of Fisheries and Aquatic Sciences, 74(5): 636–649.
  * Kotta et al. (2019) Kotta, J., Vanhatalo, J., Jänes, H., Orav-Kotta, H., Rugiu, L., Jormalainen, V., Bobsien, I., Viitasalo, M., Virtanen, E., Sandman, A. N., Isaeus, M., Leidenberger, S., Jonsson, P. R., and Johannesson, K. (2019).  “Integrating experimental and distribution data to predict future species patterns.”  Scientific Reports, 9(1): 1–14.
  * La Mere et al. (2020) La Mere, K., Mäntyniemi, S., and Haapasaari, P. (2020).  “The effects of climate change on Baltic salmon: Framing the problem in collaboration with expert stakeholders.”  The Science of the Total Environment, 738: 140068.
  * LaMere et al. (2020) LaMere, K., Mäntyniemi, S., Vanhatalo, J., and Haapasaari, P. (2020).  “Making the Most of Mental Models: Advancing the Methodology for Mental Model Elicitation and Documentation with Expert Stakeholders.”  Environmental Modelling & Software, 124: 104589.
  * Långnabba et al. (2019) Långnabba, A., Hyvönen, J., Kuningas, S., Lappalainen, A., Veneranta, L., and Kallasvuo, M. (2019).  “Evaluation of the Gulf sampling method: Report conducted in the VELMU Inventory Programme for the Underwater Marine Environment.”  Technical report, Natural resources and bioeconomy studies 2/2019, Natural Resources Institute Finland.
  * Lindley and Singpurwalla (1986) Lindley, D. V. and Singpurwalla, N. D. (1986).  “Reliability (and fault tree) analysis using expert opinions.”  Journal of the American Statistical Association, 81(393): 87–90.
  * Liu and Vanhatalo (2020) Liu, J. and Vanhatalo, J. (2020).  “Bayesian model based spatiotemporal survey designs and partially observed log Gaussian Cox process.”  Spatial Statistics, 35: 100392.
  * Murray et al. (2009) Murray, J. V., Goldizen, A. W., O’Leary, R. A., McAlpine, C. A., Possingham, H. P., and Choy, S. L. (2009).  “How useful is expert opinion for predicting the distribution of a species within and beyond the region of expertise? A case study using brush-tailed rock-wallabies Petrogale penicillata.”  Journal of Applied Ecology, 46(4): 842–851.
  * Mäkinen and Vanhatalo (2018) Mäkinen, J. and Vanhatalo, J. (2018).  “Hierarchical Bayesian model reveals the distributional shifts of Arctic marine mammals.”  Diversity and Distributions, 24(10): 1381–1394.
  * Nevalainen et al. (2018) Nevalainen, M., Helle, I., and Vanhatalo, J. (2018).  “Estimating the acute impacts of Arctic marine oil spills using expert elicitation.”  Marine Pollution Bulletin, 131: 782 – 792.
  * O’Hagan et al. (2006) O’Hagan, A., Buck, C. E., Daneshkhah, A., Eiser, J. R., Garthwaite, P. H., Jenkinson, D. J., Oakley, J. E., and Rakow, T. (2006).  Uncertain Judgements: Eliciting Experts’ Probabilities.  John Wiley & Sons.
  * Ovaskainen and Abrego (2020) Ovaskainen, O. and Abrego, N. (2020).  Joint Species Distribution Modelling With Applications in R.  Cambridge University Press.
  * O’Hagan (2019) O’Hagan, A. (2019).  “Expert Knowledge Elicitation: Subjective but Scientific.”  The American Statistician, 73(sup1): 69–81.
  * Pearce et al. (2001) Pearce, J. L., Cherry, K., Drielsma, M., Ferrier, S., and Whish, G. (2001).  “Incorporating Expert Opinion and Fine-Scale Vegetation Mapping into Statistical Models of Faunal Distribution.”  Journal of Applied Ecology, 38(2): 412–424.
  * Pearman-Gillman et al. (2020) Pearman-Gillman, S. B., Katz, J. E., Mickey, R. M., Murdoch, J. D., and Donovan, T. M. (2020).  “Predicting wildlife distribution patterns in New England USA with expert elicitation techniques.”  Global Ecology and Conservation, 21: e00853.
  * Perälä et al. (2020) Perälä, T., Vanhatalo, J., and Chrysafi, A. (2020).  “Calibrating Expert Assessments Using Hierarchical Gaussian Process Models.”  Bayesian Anal., 15(4): 1251–1280.
  * Rue et al. (2009) Rue, H., Martino, S., and Chopin, N. (2009).  “Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations (with discussion).”  Journal of the royal statistical society: Series b (statistical methodology), 71(2): 319–392.
  * Simpson et al. (2017) Simpson, D. P., Rue, H., Riebler, A., Martins, T. G., and Sørbye, S. H. (2017).  “Penalising model component complexity: A principled, practical approach to constructing priors.”  Statistical Science, 32(1): 1–28.
  * Tversky and Kahneman (1974) Tversky, A. and Kahneman, D. (1974).  “Judgment under uncertainty: Heuristics and Biases.”  Science, 185(4157): 1124–1131.
  * United Nations et al. (2021) United Nations et al. (2021).  “System of Environmental-Economic Accounting – Ecosystem Accounting (SEEA EA). White cover publication, pre-edited text subject to official editing.”  Technical report, United Nations.
URL https://seea.un.org/ecosystem-accounting

  * Vanhatalo et al. (2020) Vanhatalo, J., Hartmann, M., and Veneranta, L. (2020).  “Additive Multivariate Gaussian Processes for Joint Species Distribution Modeling with Heterogeneous Data.”  Bayesian Anal., 15(2): 415–447.
  * Vanhatalo et al. (2014) Vanhatalo, J., Vetemaa, M., Herrero, A., Aho, T., and Tiilikainen, R. (2014).  “By-Catch of Grey Seals (_Halichoerus grypus_) in Baltic Fisheries – Bayesian Analysis of Interview Survey.”  PLoS ONE, 9: e113836.
  * Veneranta et al. (2011) Veneranta, L., Urho, L., Lappalainen, A., and Kallasvuo, M. (2011).  “Turbidity characterizes the reproduction areas of pikeperch (Sander lucioperca L.) in the northen Baltic Sea.”  Estuarine, Coastal and Shelf Science.


The analogy is further limited by the fact that graph embeddings of transformers must have some (possibly arbitrary) node identifier as input in order to tokenize a graph using the node/edge encoding without losing the ability to associate each node with its incident edges.
However, the juxtaposition of the 1-WL and -based limitations on the abilities of GNNs to solve connectivity-based tasks suggests a fundamental gap in capabilities between models that is apparent in multiple theoretical lenses.

§ EXPERIMENTAL DETAILS

§.§ Datasets

We evaluate our model on the diverse graph reasoning tasks presented in GraphQA <cit.>. We used the public code of the dataset available at <https://github.com/google-research/google-research/tree/master/graphqa>. The code to generate the datasets is licensed under the Apache License, Version 2.0. The tasks in the dataset range in difficulty and encompass the following categories:

* Graph-level:
node counting (counting the number of nodes in a graph), edge counting (counting the number of edges in a graph),
cycle check (determining whether a graph contains a cycle), and
triangle counting (counting the number of triangles in a graph).

* Node-level:
node degree (calculating the degree of a given node in a graph).

* Edge-level:
connectivity (finding if there is a path from one node to another),
edge existence (whether a given edge exists in a graph, and
shortest path (finding the length of the shortest path from one node to another).

The graphs used in the experiments in this paper and the corresponding graph reasoning tasks are taken from <cit.>.
There are $1,000$ graphs in the original train set, $500$ graphs in the dev set, and $500$ graphs in the test set. The graphs are generated randomly using Erdős-Rényi (ER) random graph model <cit.>. Graph size ranges from 5 to 20 nodes.

Train set statistics. Average number of nodes: 11.90; average number of edges: 37.01; average node degree: 5.43.

Test set statistics. Average number of nodes: 12.37; average number of edges: 39.79; average node degree: 5.70.

Histogram of minimum cycle lengths for cycle check instances.

While random instances of graph reasoning tasks provide a valuable assessment of the task complexity on realistic graphs, they do not necessarily reflect the “worst case” graph inputs that convey negative results like <Ref> and <Ref>.
For example, the reduction that establishes that cycle check is “as hard as” graph connectivity and the consequential logarithmic-depth hardness results hinge on the consideration of graph instances with $n$ nodes and polynomial cycle length.
However, as witnessed by <Ref>, the shortest cycles observed in 1000 instances of GraphQA cycle check is almost always of length three, and only 3.2% of instances are larger.
As a consequence, identifying the existence of a cycle on the GraphQA dataset is inherently local, which is reflected by a strong performance by heuristic-based GNN solutions (<Ref>)—despite the fact that efficient GNNs for worst-case cycle check do not exist  (<Ref>).

For our experiments on the effect of the scale of the number of training data points in the final results we obtain, we use the open-source code of GraphQA available to generate a larger training dataset of 100K examples. We follow the original instructions and parameters to create this larger training dataset.

§.§ Implementation Details

Model Hyperparameters.
We fixed the number of iterations as 1,000,000 and train standard decoder-only transformers with $L = 12$ layers, $m = 768$ embedding dimension, $H = 12$ heads, learning rate $5 \cdot 10^{-4}$, and dropout $0.1$.
These models have an approximate parameter count of 60,000,000.

We used random search <cit.> over the following set of hyperparameters to select a universal architecture for all tasks: The range provided for the learning rate and dropout rate are $[10^{-4}, 10^{-1}]$ and $[0, 0.5]$. The number of layers $L$ and embedding dimension $m$ is selected from $L \in \set{4, 6, 8, 10, 12, 14, 16}$ and $m \in \set{192, 384, 576, 768, 960, 1152, 1344, 1536}$. We
employed the GLU <cit.> activation as a non-linearity.

Model Selection.
We implemented our model in JAX <cit.> and used AdamW <cit.> as the optimizer. Optimal hyperparameters for each task and model were determined by training on the GraphQA$_{Train}$ dataset and evaluating performance on the GraphQA$_{Dev}$ dataset. The results presented in the paper are based on the held-out GraphQA$_{Test}$ dataset.

Hardware Acceleration. All experiments were conducted using Google's TPUv3 and TPUv5e accelerators <cit.>.

§.§ Baseline Results

To rigorously evaluate the performance of transformers on graph reasoning tasks, we compare them against three established categories of baselines:

* Prompting-based methods. These methods provide the LLM with a textual descriptions of the graph and question within the prompt. We consider the following variations and copy the results from the original papers:

* zero-shot. In this approach, the model is given a task description and immediately asked to produce the desired output. No additional examples or demonstrations are provided.

* few-shot. This approach provides the model with a few examples of the task and their desired outputs <cit.>. Unlike traditional training, these examples are included directly in the prompt, allowing the model to learn and adapt during the inference.

* CoT. Chain-of-thought (CoT) prompting <cit.> provides examples each showing step-by-step reasoning, teaching the LLM to generate its own thought processes for tackling new tasks.

* zero-cot. Zero-shot CoT <cit.> builds upon Chain-of-Thought (CoT) prompting by eliminating the need for training examples.  The LLM generates its own step-by-step reasoning process using a simple trigger phrase like “Let's think step by step”.

* cot-bag. BAG prompting <cit.> extends cot to improve the performance of LLMs on graph-related tasks by appending “Let's construct a graph with the nodes and edges first” to the prompt.

* Graph-based methods. These models are specifically designed to process graphs as input and are trained task-specific. They leverage the connections between nodes to learn patterns and make predictions, making them ideal for tasks where a graph is involved. We use GCN <cit.>, MPNN <cit.>, and GIN <cit.> from this category. GraphToken <cit.> is a GNN-based model that processes the graph and feed the output of the GNN as soft-tokens to an LLM.

* Transformer models (Ours). The last class of model are task-specific vanilla transformer models <cit.>. The 60M transformer-1K model is the one described above trained on $1,000$ training examples from the GraphQA training set.
To investigate the impact of training data scale, we generated a larger dataset containing $100,000$ examples, ensuring the same distribution as the original training set by using the official GraphQA code and trained 60M transformer-100K on that. The 11B transformer (FT)-1K is a vanilla transformer model that is started with a pre-trained checkpoint of T5 <cit.> and is fine-tuned on the 1K training dataset.
We also include two fine-tuned PaLM 2 <cit.> transformers of size XXS and XS.
Similar to prompting baselines, this model receives a textual description of the graph as input to leverage its textual reasoning capabilities.

The results for zero-shot, zero-cot, few-shot, cot, and cot-bag are taken from <cit.>. Results for soft-prompt and GraphToken are sourced from <cit.>.

We independently evaluated GCN, MPNN, and GIN models on these tasks. We used the original architectures proposed in their respective papers and performed hyperparameter tuning on the GraphQA$_{Dev}$ dataset.

§.§ Further Experimental results

Comparison of train and test scaling on all tasks.

<Ref> presents a comprehensive comparison of the graph reasoning capabilities across various baseline models and our proposed transformer architectures. The results highlight several key findings, which we summarize below:

Transformers Exhibit Strong Performance on Graph-based Reasoning Problems.
While transformers are not explicitly designed for graph reasoning tasks like graph-based models, they demonstrate surprisingly strong performance in this domain. The results of this study indicate that transformers, despite their versatility as a general architecture, can often match or even surpass specialized graph models on a variety of graph reasoning benchmarks.

Transformers Excel at Retrieval Tasks. As proved in <Ref>, retrieval tasks can be solved by transformers. The obtained results confirm that such tasks are relatively easy for transformers as they obtained the full accuracy on most of such tasks. One exception here is the node degree task that GNNs outperform transformers but still transformers perform relatively well. We discuss why GNNs outperform transformers well for this task.

Larger Transformers Excel at Solving Search Tasks. As discussed in <Ref>, transformers are effective for search tasks, albeit requiring a larger number of parameters compared to retrieval tasks. This is empirically evident in the comparison between Transformer-1K and Transformer-1K (pretrained). It's worth noting that the pretrained transformer used here has 11 billion parameters, a significant increase from the 1 million parameters in Transformer-1K.

Transformers Excel at Capturing Global Patterns.
An interesting observation here is the performance gap between Transformers and GNNs across tasks with varying emphasis on local versus global graph structure. Notably:

* Local Structure: The node degree task, which relies heavily on local node information, is best handled by MPNN, a GNN-based model.

* Global Structure: In contrast, tasks like connectivity, triangle counting and shortest path, which require understanding global graph patterns, are dominated by Transformer models.
Notably, vanilla Transformers achieve a remarkable 45% relative improvement over even much larger LLMs augmented with GNN-generated soft prompts (GraphToken) on the shortest path task. This showcases the exceptional ability of Transformers to capture long-range dependencies, a critical factor in understanding global graph structure.

§.§.§ Sample complexity ablations

To develop a more comprehensive understanding of graph reasoning tasks learnability by small transformers, we train a variety of transformers for 1,000,000 steps on each task on a range of sample sizes.
In <Ref>, we demonstrate how the model performance improves as a function of sample complexity.
By doing so, we witness the relative hardness of each task in terms of the marginal benefit of new samples.

* Edge count and node count are “easy” retrieval tasks that can be solved perfectly with as few as 100 samples.

* Edge existence attains near-perfect train and test classification accuracy as the number of training samples approaches 100,000.

* Connectivity, shortest path,
and node degree demonstrate a sharp improvement in evaluation error as a function of the sample size. These models perfectly fit the training set in most sample size regimes, but yield a closer correspondence between training and testing error when trained on 100,000 samples.

* Cycle check and triangle count have persistent gaps between training and testing error and overfit even in the large sample setting.

§ NEURIPS PAPER CHECKLIST

* Claims

Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?

Justification: The empirical and theoretical claims discussed accurately represent the detailed theoretical analysis and the wide-ranging experimentation on graph reasoning benchmarks.

* The answer NA means that the abstract and introduction do not include the claims made in the paper.

* The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.

* The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.

* It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.

* Limitations

Question: Does the paper discuss the limitations of the work performed by the authors?

Justification: Theoretical limitations are discussed in the presentation of the results in <Ref> and detailed more comprehensively in <Ref>. We discuss limitations of experimental results in <Ref>.

* The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.

* The authors are encouraged to create a separate "Limitations" section in their paper.

* The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.

* The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.

* The authors should reflect on the factors that influence the performance of the approach. For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.

* The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.

* If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.

* While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.

* Theory Assumptions and Proofs

Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?

Justification: The theoretical results follow all of the criteria. The proofs are in Appendix.

* The answer NA means that the paper does not include theoretical results.

* All the theorems, formulas, and proofs in the paper should be numbered and cross-referenced.

* All assumptions should be clearly stated or referenced in the statement of any theorems.

* The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition.

* Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.

* Theorems and Lemmas that the proof relies upon should be properly referenced.

* Experimental Result Reproducibility

Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?

Justification: We comprehensively detail the models used and the training hyper-parameters in the appendix. The training dataset is generated by a publicly available source code.

* The answer NA means that the paper does not include experiments.

* If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.

* If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.

* Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.

* While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example

* If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.

* If the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.

* If the contribution is a new model (e.g., a large language model), then there should either be a way to access this model for reproducing the results or a way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).

* We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.

* Open access to data and code

Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?

Justification: While code is not made available, the data can be generated using the GraphQA dataset linked in the appendix. Full training details are made available, which makes the code easily reproducible. The code will be open sourced upon acceptance of the paper.

* The answer NA means that paper does not include experiments requiring code.

* Please see the NeurIPS code and data submission guidelines (<https://nips.cc/public/guides/CodeSubmissionPolicy>) for more details.

* While we encourage the release of code and data, we understand that this might not be possible, so “No” is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).

* The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (<https://nips.cc/public/guides/CodeSubmissionPolicy>) for more details.

* The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.

* The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset of experiments are reproducible, they should state which ones are omitted from the script and why.

* At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).

* Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.

* Experimental Setting/Details

Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?

Justification: All training and evaluation details and model hyper-parameters are included in the paper appendix.

* The answer NA means that the paper does not include experiments.

* The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.

* The full details can be provided either with the code, in appendix, or as supplemental material.

* Experiment Statistical Significance

Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?

Justification: This is due to the computational cost of certain experiments, particularly those involving LLMs.

* The answer NA means that the paper does not include experiments.

* The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.

* The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).

* The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)

* The assumptions made should be given (e.g., Normally distributed errors).

* It should be clear whether the error bar is the standard deviation or the standard error of the mean.

* It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified.

* For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).

* If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.

* Experiments Compute Resources

Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?

Justification: The type of the accelerator used in the experiments is in Appendix.

* The answer NA means that the paper does not include experiments.

* The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.

* The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.

* The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).

* Code Of Ethics

Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics <https://neurips.cc/public/EthicsGuidelines>?

Justification: We reviewed the guidelines and found no relevant violations to our work. Since the experimental datasets are synthetically generated, privacy concerns are moot.

* The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.

* If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.

* The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction).

* Broader Impacts

Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?

Justification: The paper is focused on reasoning capabilities of transformers using theoretical and empirical evaluation. To the best of our knowledge, there is no societal impact for the work.

* The answer NA means that there is no societal impact of the work performed.

* If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.

* Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.

* The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.

* The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.

* If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).

* Safeguards

Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?

Justification: This is because the datasets used are synthetically generated and reflect random graphs, no risks of data leakage are present.

* The answer NA means that the paper poses no such risks.

* Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.

* Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.

* We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.

* Licenses for existing assets

Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?

Justification: The papers are cited, the link to the source codes used is states along with their license.

* The answer NA means that the paper does not use existing assets.

* The authors should cite the original paper that produced the code package or dataset.

* The authors should state which version of the asset is used and, if possible, include a URL.

* The name of the license (e.g., CC-BY 4.0) should be included for each asset.

* For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.

* If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, <paperswithcode.com/datasets> has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset.

* For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.

* If this information is not available online, the authors are encouraged to reach out to the asset's creators.

* New Assets

Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?

Justification: No new assets is introduced in this paper. We have used public data and no code has been released yet.

* The answer NA means that the paper does not release new assets.

* Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.

* The paper should discuss whether and how consent was obtained from people whose asset is used.

* At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.

* Crowdsourcing and Research with Human Subjects

Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?

Justification: No human subjects were involved.

* The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.

* Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper.

* According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.

* Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects

Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?

Justification: No human subjects were involved.

* The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.

* Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.

* We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.

* For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.