# A marked sequential point process for disease surveillance: Modeling and optimization

François D'alayer, Edith Gabriel, Samuel Soubeyrand

► To cite this version:

François D'alayer, Edith Gabriel, Samuel Soubeyrand. A marked sequential point process for disease surveillance: Modeling and optimization. Spatial Statistics, 2025, 68, pp.100913. 10.1016/j.spasta.2025.100913. hal-05490031

HAL Id: hal-05490031

https://hal.inrae.fr/hal-05490031v1

Submitted on 2 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

ELSEVIER

Contents lists available at ScienceDirect

Spatial Statistics

journal homepage: www.elsevier.com/locate/spasta

<figure><img src="image_3.png" /></figure>

# A marked sequential point process for disease surveillance:
Modeling and optimization

François d'Alayer ID*, Edith Gabriel, Samuel Soubeyrand

INRAE, BioSP, Avignon, 84914, France

<figure><img src="image_4.png" /></figure>

Check for
updates

## ARTICLE INFO

**Keywords:**

Marked sequential point process
Disease surveillance for prevalence estimation
Self-interaction
Adaptive design
Sequential sampling

## ABSTRACT

Plant disease surveillance is essential for the management of disease outbreaks that pose significant threats to agricultural sustainability. In this study, we present a novel sequential point process model designed for disease surveillance. The model incorporates self-interaction mechanisms to account for the influence of the process' history. To analyze the dynamics of the model, we propose new sequential summary statistics that extend traditional point process methods to scenarios where sequential interactions are critical. This model serves a dual purpose: it is employed both to propose novel and efficient sampling designs, and to characterize existing sampling schemes, implemented in real-world situations, through parameter inference.

## 1. Introduction

Surveillance plays an important role in various domains such as security, health, transport, and industry. Its perimeter varies with the domain, but surveillance is often viewed as a mechanism or a system that involves monitoring, gathering, and analyzing data, to identify and respond to emerging threats or, more generally, to detect beneficial or detrimental changes within the surveyed sector. In agriculture, surveillance in general and sampling in particular are fundamental for protecting plant health (Perry, 1994; Morris et al., 2022). From this perspective, the outputs of surveillance are used for (i) limiting the spread of pathogens and pests in crops, (ii) reducing the need for widespread harsh control measures, (iii) acquiring knowledge on the ecology and the epidemiology of pathogen and pests as well as the risk they carry, and (iv) timely sharing information about the sanitary situations with diverse stakeholders who must make decisions to mitigate risks and reduce yield losses (Soubeyrand et al., 2024). Overall, plant health surveillance not only helps in minimizing economic losses in the agricultural sector (Fuller et al., 2020), but also contributes to establishing sustainable agricultural practices while ensuring crop productivity and food security.

The operational objectives of plant health surveillance include early detection of diseases or pests, estimation of their prevalence, and delimitation of their spatial distribution (Parnell et al., 2017; Mastin et al., 2020). These objectives depend on the status of the disease or the pest of interest, i.e., whether it is not present yet in the geographic area of concern, emerging or endemic. Moreover, they may be directed toward diverse downstream tasks, such as: ensuring disease-free status for some areas, eradicating quarantine pests, estimating characteristics such as dispersal capacities, and tuning control measures. Surveillance is grounded on observational activities, data collection, and data analysis aimed at assessing emergence, prevalence, spread, or impact of biotic threats affecting crops, where the observational activities can involve very diverse approaches, e.g., visual inspections, molecular diagnostics, remote sensing, textual data mining, and can target diverse observation matrices, e.g., the plants themselves, insects vectorizing pathogens, insects feeding on plants, water, air, soil, social media content (Soubeyrand et al., 2024). In this article, we focus on surveillance

<footnote>* Corresponding author.
E-mail address: francois.dalayer@inrae.fr (F. d'Alayer).</footnote>
approaches that consist in directly observing plants and testing them for the presence of a specific disease: a plant is then identified by its spatial location and its sanitary status.

In epidemiology, as in many other domains, surveillance can be biased towards individuals or situations that are considered to be at risk. In the plant health context, risk-based surveillance was defined by Morris et al. (2022) as a surveillance approach in which the sampling intensity (in a location, time, or both) depends on the local measure of risk (e.g., risk of emergence, risk of presence, risk of outbreak). This approach enables more targeted and efficient allocation of resources and interventions (Cameron, 2012; Parnell et al., 2014). In risk-based surveillance, key risk factors associated with the transmission of the disease under study are generally considered. These risk factors may include population density, demographic variables, travel or dispersal patterns, climatic indicators, land-use information, and socioeconomic data, depending on the specific pathosystem being studied. For instance, for a vector-borne disease like malaria, risk factors might include proximity to stagnant water bodies and high prevalence of mosquito vectors. Regardless the availability of such risk factors, prior observations of infected and uninfected individuals from previous surveillance surveys can be particularly useful for designing risk-based surveillance. Indeed, such prior observations are expected to contain important information on the risk of infection, at least for infectious diseases with propagation speeds that are relatively slow compared to the rate of new samples acquisition.

The utilization of prior observations of infected and uninfected individuals for designing risk-based surveillance constitutes the main topic of this article. In particular, we aim to explore how marked sequential spatial point processes can be leveraged and optimized to conceive adaptive risk-based sampling, with the ultimate objective of estimating disease prevalence. In this context, adaptivity encompasses two complementary aspects: first, new samples may be selected based on previously collected data; second, the surveillance practitioner may retain the flexibility to adjust the survey design at a given time in order to optimize the surveillance process. The former corresponds to the joint distribution of the sampling locations, while the latter refers to the capacity to choose or adapt this distribution dynamically. Our approach consists in using optimal design and point process theory to determine successive locations of plants to be tested for the presence of the disease, thereby achieving efficient plant monitoring for disease prevalence estimation.

Section 2 highlights the principles on which our approach to design optimal surveillance is grounded. Section 3 formalizes the key components used in our approach: surveillance schemes, prevalence maps and information maps. In Section 4, we present the approach for designing optimal surveillance schemes to estimate spatial prevalence. In Section 5, we present a parametric sequential point process model for surveillance inspired by actual surveillance schemes implemented to monitor regulated diseases affecting crops. This model can be fitted to data and analyzed to derive optimal parameters for estimating disease prevalence. Section 6 provides an analysis of a case study on *Xylella fastidiosa*, a bacterium that causes plant diseases.

## 2. Methodological principles for optimizing disease surveillance

Here, we describe how our work relates to optimal design and point processes. The theoretical foundations of optimal design are presented by Pukelsheim (2006), while Atkinson et al. (2007) discusses various applications. We consider a population of plants in a window $W \subset \mathbb{R}^2$. We represent individuals by their spatial locations $y$ and their sanitary status $m_y \in \{0, 1\}$, where $m_y = 1$ indicates the presence of the disease, and $m_y = 0$ indicates its absence. We are interested in a statistical property of this disease indicator, typically the spatial function representing the expectation of $m_y$ for an individual (real or hypothetical) located at $y$, referred to as the spatial prevalence and denoted by $\alpha$. However, due to practical limitations, we do not have access to all individual marks $m_y$ within the population, and must select a subset of individuals $(y_1, \dots, y_n)$ from $W$ to observe their marks $(m_{1}, \dots, m_{n})$, and estimate $\alpha$ based on these observations and locations. In optimal design, a common objective is to determine sampling configurations that yield the most informative data with respect to a specified statistical criterion (Pukelsheim, 2006). These criteria are frequently expressed in terms of the eigenvalues of the so-called information matrix, which quantifies the amount of information that the observed data provide about the model parameters. Among these criteria, A-optimality aims to minimize the trace of the inverse of the information matrix. This corresponds to reducing the average variance of the estimated parameters, thereby enhancing the overall precision of the estimation procedure. Analogously, for our surveillance objective, we want to sample the population such as the uncertainty about the spatial disease prevalence estimate $\alpha_n$ of $\alpha$ is minimized.

To formalize this idea, we quantify the spatial uncertainty about $\alpha_n$ using the Integrated Bernoulli Variance (IBV) criterion, as previously proposed by Anyosa et al. (2023) for benthic habitat mapping and by Ge et al. (2023) for water mass classification. Hence, using optimal design, we can determine optimal surveillance schemes for disease prevalence estimation. However, this alone does not reveal what defines a good surveillance protocol in the sense that the optimal design yields a surveillance strategy whose efficiency is not easily interpretable. Modeling the surveillance process as a marked point process and computing appropriate summary statistics that effectively capture the properties of the surveillance process, such as point aggregation or inter-point distances allow us to identify the attributes of an optimal scheme. Indeed, analyzing the summary statistics of an efficient surveillance process may yield insights into its underlying properties. Conversely, one could propose adjustments to an existing surveillance strategy based on its summary statistics in order to increase its efficiency.

A point pattern is an unordered set of points within a subset of $\mathbb{R}^d$, $d \in \mathbb{N}^*$, usually a bounded observation window $W \subset \mathbb{R}^d$. Hereafter, we focus on the specific case where $d = 2$. Point process models are stochastic generators of random point patterns. A point process $Z$ can be defined either as a random set of points, denoted by $Z = \{y_1, y_2, \dots\}$, or as a random counting measure. In the latter case, for every Borel<fnref n="1" /> set $B \subset \mathbb{R}^2$, the notation $N(B)$ denotes the random number of points from $Z$ that fall within the set

<footnote n="1">¹ A Borel set $B \subset \mathbb{R}^2$ is any set that belongs to the smallest $\sigma$-algebra containing all open subsets of $\mathbb{R}^2$.</footnote>
B. N is then referred to as a random counting measure. Formal and comprehensive presentations of point processes are provided by Daley and Vere-Jones (2003), Illian et al. (2008), and Chiu et al. (2013). If a random mark, say $m_{y_i}$, is associated with each point $y_i \in Z$, then the set of pairs $P = \{(y_i, m_{y_i})\}_{i}$ forms a marked spatial point process, where $\{m_{y_i}\}$ is the set of marks taking values in the mark space,<fnref n="2" /> which in this article is simply $\{0, 1\}$.

As mentioned in the introduction, we are interested in marked sequential spatial point processes to design adaptive risk-based surveillance schemes. A sequential spatial point process (SSPP) is defined as an ordered sequence of spatial locations, denoted $\vec{Z} = (y_1, ..., y_n)$, where $n$ is a positive integer (Van Lieshout, 2006b,a; Penttinen and Ylitalo, 2016). Modeling surveillance with this type of process is intuitive since, in practice, the observation of a disease dynamics is often distributed throughout time, even if it is not adaptive. Typically, when a positive (i.e., infected) individual is detected at $y$, additional individuals are sampled around $y$, following concentric or spiral patterns, with the parameters being radii or polar equation of the spiral (Picard et al., 2019). The starting points of these concentric patterns, from which subsequent plants are sampled, can be determined with the RIsk Based Estimate of System Sensitivity (RiBESS) tool (European Food Safety Authority (EFSA), 2012). This tool calculates the number of samples required in an area to ensure a sufficient level of confidence that the area is disease-free. While spatio-temporal point processes allow the time of each point to be a random variable, SSPPs assign a predetermined temporal order to the points, making them naturally suited for sequential sampling designs where the timing of observations is fixed by the sampling procedure itself.

We extend the SSPP model proposed by Penttinen and Ylitalo (2016) to marked point process models tailored for epidemiological surveillance of plant diseases. Our proposed adaptive sampling model determines new samples based on the locations of previous points and their corresponding disease status. This model can be used to create diverse sampling designs according to the chosen objective of epidemiological surveillance. Moreover, we use this model to characterize surveillance conducted in real scenarios and to compare it with a more efficient strategy, and illustrate its application using presence-absence data of *Xylella fastidiosa* in southern France, a bacterial pathogen for plants.

## 3. Surveillance

As mentioned in the introduction, our objective is to estimate the spatial variation in disease prevalence, defined as the local proportion of infected individuals, disregarding whether they were recently infected or not (by contrast, incidence refers to the proportion of newly infected individuals within a given time step). To achieve this, we assume that, throughout the surveillance campaign, the disease distribution remains stable or, more formally:

(A1) The spatial prevalence function $\alpha$ (as defined in **Definition 3.2**) remains constant over time.

### 3.1. Surveillance schemes

Surveillance requires the specification of a strategy for sampling individuals across a surveillance region $W \subset \mathbb{R}^2$, with each individual associated with a point in $W$. The region $W$ is defined as the set of locations where the population intensity is strictly positive; that is, the probability of observing an individual at any given location within $W$ is non-zero. Once individuals to be sampled are selected, their infection statuses are assessed using a method (e.g. visual inspection or molecular test) that classifies each individual as either healthy or infected. Sampling locations can be represented and analyzed through spatial point process models. Sampling locations and corresponding infection statuses can be represented and analyzed using marked spatial point process models.

**Definition 3.1.** A surveillance scheme is defined as a marked point process $P_n = \{(y_i, m_i) : y_i \in W, m_i \in \{0, 1\}, 1 \le i \le n\}$ of size $n$, where the points are located in $W$ and the marks represent the infection statuses (1: infected, 0: healthy).

A surveillance scheme $P_n$ is said to be sequential if, for all $2 \le i \le n$, $y_i$ is defined conditionally on $(y_1, m_1), \dots, (y_{i-1}, m_{i-1})$. In this case, it is denoted as $\bar{P}_n = ((y_1, m_1), \dots, (y_n, m_n))$. Otherwise, $P_n$ is said to be non-sequential.

#### Example 3.1.

* The Binomial point process of size $n$ on $W$, augmented with marks independently drawn from Bernoulli distributions with success probabilities $\{q_i : i = 1, \dots, n\}$, is a surveillance scheme.

* Given a distance $r > 0$, the point process where points $y_i$ are drawn iteratively for increasing $i$, such that $y_i$ ($i > 1$) lies outside discs of radius $r$ centered at the previously selected points $(y_1, \dots, y_{i-1})$, and where points are marked by independent draws from Bernoulli distributions as described above, is a sequential surveillance scheme.

<footnote n="2">² The mark space is generally assumed to be a Polish space, i.e., a complete, separable metric space, which contains the case of discrete marks 0 or 1.</footnote>
**What is adaptivity.** Adaptivity in surveillance schemes refers to the property whereby the selection of new sampling units depends on those previously selected.

More formally, the adaptivity of a sample $P_n = \{(y_1, m_1), \dots, (y_n, m_n)\}$ can be characterized by its non-exchangeability. This means that there exists at least one permutation $\sigma$ of the indices $\{1, \dots, n\}$ such that the ordered sequence $P_n^\sigma = \{(y_{\sigma(1)}, m_{\sigma(1)}), \dots, (y_{\sigma(n)}, m_{\sigma(n)})\}$ does not follow the same probability distribution as the original sample $P_n$. In other words, the ordering of the samples matters, as each sample may influence subsequent ones, leading to a form of dependence among observations.

The degree of adaptivity can vary across different surveillance schemes. A minimal form of adaptivity is Markovian adaptivity, in which each new sample depends only on the most recently observed sample. A stronger form is full adaptivity, where each new sample is conditioned on all observed samples. A common and important class of adaptive schemes is sequential schemes, in which every new sampling decision is influenced by all prior observations.

In contrast, non-adaptive survey schemes are those where the selection of sampling units is independent of other selected units. Examples of such schemes include uniform sampling, where each unit is selected independently with equal probability, and grid sampling, where samples are placed at predetermined fixed locations. Conversely, adaptive sampling schemes, such as two-phase sampling or cluster sampling, incorporate dependencies between sampled units. For instance, in two-phase sampling, an initial sample is used to guide the selection of a subsequent sample, while in cluster sampling, entire groups of units are sampled together, leading to interdependencies within clusters.

In this paper, we introduce a novel parametric model for designing surveillance schemes, encompassing both adaptive and non-adaptive strategies, as well as configurations exhibiting attraction or repulsion between sampling locations. This framework allows for the inference and creation of diverse surveillance protocols. While the schemes generated by the model are inherently sequential, parameter inference is conducted using a maximum likelihood approach.

## 3.2. Prevalence maps

In this article, surveillance schemes are designed to improve the estimation of disease prevalence and its spatial variation. Thus, here, we define the prevalence map and give an example of its estimation based on data collected from a surveillance scheme, using a smoothing operator.

**Definition 3.2.** The prevalence map is the function $\alpha: W \to [0, 1]$ that assigns to each point in $W$ the probability of an individual at this location to be infected.

By the way, we make the following assumption related to the conditional independence of the marks and their distributions:

(A2) Given $\alpha$, the infection statuses of any individual located at $u \in W$ are independent and follow Bernoulli distributions with success probabilities $\alpha(u)$.

Consider a surveillance scheme $P_n = ((y_1, m_1), \dots, (y_n, m_n))$, the prevalence map $\alpha$ can be estimated using various methods, for example, with the following kernel estimator:

$$ \alpha_{P_n}(u) = \frac{\sum_{i=1}^{n} K(y_i, u)m_i}{\sum_{i=1}^{n} K(y_i, u)} \mathbb{1}_{\sum_{i=1}^{n} K(y_i, u)>0}, \quad \forall u \in W, \qquad (1) $$

where $K$ is the kernel function.<fnref n="3" /> $K$ may be specified to mitigate border effects or to account for spatially varying point densities in $P_n$. In addition, if the sampling of individuals does not follow a uniform distribution across the population, the resulting sampling bias should be corrected using importance sampling weights (Kloek and Van Dijk, 1978). This technique is particularly relevant in the case of sequential surveillance schemes.

## 3.3. Information maps

Each sampled individual, regardless of its infection status, provides spatial information. Here, we summarize the contributions of sampled individuals to the information through a function $\iota$, called the information map.

**Definition 3.3.** An information map is a function $\iota: W \to [0, 1]$ that represents the local amount of the available information. It assigns large values (close to 1) where the vast majority of sampled individuals are infected, low values (close to 0) where the vast majority are healthy, and intermediate values (around 0.5) in other cases or in areas with very few sampled individuals.

For any surveillance scheme $P_n = ((y_1, m_1), \dots, (y_n, m_n))$, whether sequential or non-sequential, we propose the following definition of the information map $\iota_{P_n}$ given $P_n$:

$$ \iota_{P_n}(u) = \frac{1}{2} + \frac{\sum_{i=1}^{n} K'(y_i, u)(2m_i - 1)}{2 \sum_{i=1}^{n} K'(y_i, u)} \mathbb{1}_{\sum_{i=1}^{n} K'(y_i, u)>0} \qquad (2) $$

<footnote n="3">³ In this article, a kernel will be any non-negative function $K: W \times W \to \mathbb{R}_+$.</footnote>
where $K'$ is a kernel. For any standard kernel function whose output value (sharply) decreases with the distance between $y_i$ and $u$, $t_{P_n}$ approaches 1 in the vicinity of infected sampled individuals ($y_i, 1$), whereas it tends to 0 near healthy sampled individuals ($y_i, 0$), and to 1/2 in areas with no sampled individuals or at any point $u$ equidistant from two sampled individuals with marks 0 and 1, respectively.

### 4. Optimal design for prevalence estimation

Here, the objective is to design a surveillance scheme $P_n$ that leads to an accurate estimation of the prevalence map $\alpha$ over $W$. For instance, it could aim to minimize the $L^2$-norm error:

$$ \int_W (\alpha(u) - \alpha_{P_n}(u))^2 du, \qquad (3) $$

where $\alpha_{P_n}$ is the estimate of $\alpha$ derived from $P_n$, regardless of whether the surveillance scheme is sequential or non-sequential.

Directly reducing the error expressed in Eq. (3) is not possible because $\alpha$ is unknown. Therefore, we propose to reduce another quantity that is expected to have the same behavior. Our proposal is based on the following statement: we are primarily looking for a way to collect data (i.e., specifying $P_n$), and, only indirectly, a way to estimate a quantity (i.e., estimating $\alpha$). Heuristically, reducing the error in Eq. (3) corresponds to reducing the local uncertainties resulting from the gaps in the collection of data. These local uncertainties can be quantified using the information map introduced in Section 3.3. Thus, we want to minimize the Integrated Bernoulli Variance (IBV) of the information map:

$$ \text{IBV}(t_{P_n}) = \int_W t_{P_n}(u)(1 - t_{P_n}(u))du, \qquad (4) $$

which may be considered as an approximate way of minimizing the error (3). One might question why we choose to minimize IBV($t_{P_n}$) instead of IBV($\alpha_{P_n}$). The rationale lies in the fact that the information function $t$ is a transformation of the prevalence map that exhibits more desirable properties. Specifically, the function $t_{P_n}(u)(1 - t_{P_n}(u))$ reaches its maximum when $t_{P_n}(u) = 1/2$. In an unsampled region, we have $t_{P_n} = 1/2$, whereas $\alpha_{P_n} = 0$. Consequently, sampling in such a region reduces IBV($t_{P_n}$), which aligns with obtaining information from all regions. In contrast, IBV($\alpha_{P_n}$) may increase when sampling in that region, making it a less suitable criterion for our purposes. Moreover, under uniform random sampling and with an increasing sample size, both $t_{P_n}$ and $\alpha_{P_n}$ converge to the true prevalence map $\alpha$.

As we will see below (with Conjecture 1, which states that the minimization of the $L^2$-norm error consists in making a compromise between sampling in regions with high disease prevalence and sampling in a manner that reduces the IBV), it is necessary to add a constraint to the minimization of the IBV to better approximate the minimization of the $L^2$-norm error expressed in Eq. (3). Let us first heuristically motivate the choice of the IBV with a simple example. Consider the following experiment: Suppose we aim to estimate the unknown parameters $p_1, \dots, p_k \in [0, 1]$ for which independent samples can be drawn from Bernoulli distributions with success probabilities $p_1, \dots, p_k$. The variance of the empirical mean used to estimate of $p_i$ is $(1/N_i)p_i(1 - p_i)$, where $N_i$ is the sample size for the $i$th parameter. This variance increases and then decreases as $p_i$ varies from 0 to 1, reaching its maximum when $p_i = 1/2$. Therefore, to achieve the same variance (or uncertainty) for all the estimates, the sampling size $N_i$ decreases with the absolute difference between $p_i$ and $1/2$.

Beyond this simple example, the IBV criterion is commonly encountered in the optimal design literature (Anyosa et al., 2023; Ge et al., 2023; Zhao et al., 2023). IBV quantifies prediction uncertainty of a model, accounting for both epistemic uncertainty (arising from model approximation) and random uncertainty (stemming from the stochastic nature of the data). IBV computes the variance of the predicted probability distribution for each class over the dataset. A high IBV reflects uncertainty in the model's predictions, making it useful for identifying cases where the model is poorly reliable. Conversely, a low IBV suggests that the model generates stable and reliable predictions. IBV serves as a valuable tool for assessing the confidence in the predictions of a classification model. Optimizing this criterion is similar to the I-optimal design in the optimal design theory. In this context, an I-optimal design seeks to minimize the average prediction variance over the entire design space. In other words, rather than focusing solely on estimating model parameters precisely (as in D-optimal designs), I-optimality aims to ensure that the model yields accurate predictions throughout the region of interest.

We now propose the following conjecture:

**Conjecture 1.** Let $n \ge 2$, and let $S$ denote the set of surveillance schemes of size $n$ that include at least one infected and one non-infected individual. Then, minimizing the IBV across surveillance schemes while maximizing the number of infected samples serves as a surrogate for minimizing the $L^2$-norm error (3) across surveillance schemes.

**Reasoning 1.** For $P_n \in S$, we note $P_n = (y_1, m_1), \dots, (y_n, m_n)$ and $\alpha_{P_n}(u) = \frac{\sum_{i=1}^n K(y_i, u)m_i}{\sum_{i=1}^n K(y_i, u)} = \sum_{i=1}^n \omega_i(u)m_i$ with $\omega_i(u) = \frac{K(y_i, u)}{\sum_{j=1}^n K(y_j, u)}$. We consider the objective of prevalence estimation (3). As the prevalence map is a spatial smoothing of the sanitary statuses across all the individuals, we assume that it only depends on the individuals and can be written as

$$ \alpha(u) = \sum_{y \in D} \tilde{\omega}_y(u)m_y $$

where $D$ is the set of all individuals locations and $m_y$ is the sanitary status of the individual at location $y$. $\omega_i$ and $\tilde{\omega}_{y_i}$ are not necessarily the same. An example is $\tilde{\omega}_y(u) = \frac{K'(y, u)}{\sum_{z \in D} K'(z, u)}$ where $K'$ is a kernel.
We then have to minimize the integral of

$$ 
\begin{aligned}
& \left( \sum_{y \in D} \tilde{\omega}_y(u) m_y - \sum_{i=1}^{n} \omega_i(u) m_i \right)^2 = \\
& \left( \sum_{m_y=1}^{n} \tilde{\omega}_y(u)^2 + \sum_{m_i=1}^{n} \omega_i(u)^2 + \sum_{\substack{i \neq j \\ m_i=m_j}}^{n} \omega_i(u) \omega_j(u) - 2 \sum_{m_y=1}^{n} \sum_{m_i=1}^{n} \tilde{\omega}_y(u) \omega_i(u) \right)
\end{aligned}
$$

This is equivalent to maximizing the integral of:

$$ 
\begin{aligned}
& 2 \sum_{m_y=1}^{n} \sum_{m_i=1}^{n} \tilde{\omega}_y(u) \omega_i(u) + 2 \sum_{m_i=1}^{n} \tilde{\omega}_{y_i}(u) \omega_i(u) - \sum_{m_i=1}^{n} \omega_i(u)^2 + \\
& 2 \sum_{\substack{i \neq j \\ m_i=m_j}} \tilde{\omega}_{y_i}(u) \omega_j(u) - \sum_{\substack{i \neq j \\ m_i=m_j}} \omega_i(u) \omega_j(u)
\end{aligned}
\quad (5)
$$

Maximizing the part $\sum_{y \in D \setminus P_n} \sum_{m_i=1}^{n} \tilde{\omega}_y(u) \omega_i(u)$ of (5) corresponds to sampling in high risk regions, as it is a non-decreasing function of the number of non-sampled infected individuals around each positive sample $(y_i, 1)$, whereas maximizing the remaining terms of (5) corresponds to sampling in a way that maximizes information. The reason is that $\alpha$ is built based on all individuals in D so that it uses a smaller scale than $\alpha_{P_n}$ that considers only the sample $P_n$. For instance, if $\tilde{\omega}_y(u) = \frac{K'(y, u)}{\sum_{z \in D} K'(z, u)}$, the bandwidth of $K'$ would typically be smaller than that of $K$, which explains the two behaviors observed in (5).

Thus, the maximization is a compromise between sampling in high risk regions to increase the number of infected individuals sampled and maximizing the information to minimize the IBV. ☐

From Conjecture 1, minimizing information uncertainty while maximizing positive marks in surveillance schemes is an effective approach to determine samples locations for prevalence estimation. Furthermore, multiple optimal surveillance schemes for prevalence estimation may exist, some being clustered and others promoting exploration. When comparing such schemes, preference is given to those that not only optimize prevalence estimation but also ensure spatial exploration (Marjou et al., 2019).

Thus, we focus on reducing information uncertainty while maximizing the number of infected samples as an equivalent manner to get precise prevalence estimation:

$$ \underset{P_n \in S}{\operatorname{argmin}} \left\{ \operatorname{IBV}(t_{P_n}) \text{ while maximizing } \sum_{i=1}^{n} m_i \right\} \quad (6) $$

where $S$ is the set of surveillance schemes. This minimization ensures both accurate prevalence estimation and effective space exploration.

In this work, we propose to solve this problem by modeling the surveillance scheme as a marked sequential point process. At each iterative step $i$, we propose a new surveillance location $y_{i+1}$ based on all previous samples $\vec{P}_i = (y_1, m_1), \dots, (y_i, m_i)$. Regions identified as high risk in the past are favored to maximize the number of infected samples, while some parameters allow to favor regions with higher information content to reduce information uncertainty.

## 5. Sequential surveillance scheme

A sequential spatial point process is generally defined through successive conditioning (Penttinen and Ylitalo, 2016): in a sequence $\vec{Z}_k = (y_1, \dots, y_k)$ of points in $W \subset \mathbb{R}^2$, each point $y_i$ depends on all preceding points $\vec{Z}_{i-1}$. This dependency can be expressed through conditional densities. The density $f$ of a sequence is defined as:

$$ f(\vec{Z}_k) = f_1(y_1) \prod_{i=1}^{k-1} f_{i+1}(y_{i+1} | \vec{Z}_i) $$

where $f_{i+1}(y_{i+1} | \vec{Z}_i)$ is the transition probability from $\vec{Z}_i$ to $(\vec{Z}_i, y_{i+1})$. In Penttinen and Ylitalo (2016), a sequential spatial point process is used to model the eye movements of a person observing a painting. When using such a process, fundamental characteristics of the context and problem should be incorporated into the conditional densities. For eye movement data, fundamental characteristics include:

* **Spatial heterogeneity:** The intensity of the process governing the position of the eyes on the painting over time is not constant. Typically, the center of the painting usually gets more interest than the border. This characteristic is modeled by a saliency map $\alpha$.

* **Short term dependency:** Observations of jump lengths between successive fixations reveal some dependency. This dependency is captured in the model using a kernel $k$ that describes Markovian properties.

* **Long term dependency:** Since a painting is initially inspected globally before the viewer focus on details, a long term dependency is incorporated into the model. This is achieved through a reweighting function $\pi$. Each point $y_i$ is assigned a weight based on its position relative to the convex hull generated by $\vec{Z}_{i-1}$.
The conditional densities are then created to account for these characteristics. In some ways, epidemiological surveillance is similar to the eye movement on a painting (initially inspecting the surveillance area globally, then focusing locally based on risk factors and findings), we generalize this model to our problem, incorporating additional information available for each point: its sanitary status.

## 5.1. Construction

When modeling surveillance, several key characteristics must be considered, including the spatial heterogeneity of the observation area and the dependence between surveillance points. We propose a marked sequential point process model, hereafter referred to as *Sequential Surveillance Scheme (SSS)*, which incorporates these characteristics through three components that reflect the progression of epidemiological surveillance.

It offers: (i) the ability to intensify surveillance in high-risk areas with a risk map, (ii) Markovian properties via a proposal kernel, and (iii) a re-balancing function to learn from the past by integrating interaction with the process history and favoring areas where the available information exhibits specific characteristics.

**Spatial heterogeneity.** The first key characteristic is spatial heterogeneity, indicating that the importance of certain regions is not uniform. This heterogeneity is captured through a risk map $\alpha_r$, which adjusts the importance assigned to specific areas. It is intended to model the probability of presence of the disease. The risk map can be constructed either by using bioclimatic risk factors identified for the disease (e.g., elevation, proximity to specific features, vegetation type) or by developing an empirical prevalence map derived from previous surveillance efforts. Thus, the density of surveillance points is higher in regions where the risk map indicates elevated risk.

**Dependence between consecutive surveillance points.** The second key characteristic is the dependence between consecutive surveillance points, enabling control over the distance between successive surveillance points. This dependence is modeled by assuming that a point $y_{k+1}$ depends on the preceding point $y_k$. A proposal kernel, $k(y_k, y_{k+1})$, such as a Gaussian kernel, is used to model this relationship. It favors regions close to the last point, thus providing a sequential aspect to the surveillance scheme.

**Long-term dependence/self-interaction.** Beyond Markovian<fnref n="4" /> dependence, the model accounts for a memory effect, meaning that information gathered by the process influences the generation of new points. This is achieved through a rebalancing factor that enables self-interaction within the process. The re-balancing factor adjusts the density of points based on local information, penalizing areas with low values and reducing their density. For instance, areas with high disease intensity can be favored by increasing the density in regions where the information map approaches 1.

As previously discussed, the definition of the density $f$ of a sequential point process relies on its conditional density:

$$g(\vec{P}_n) = g_1(P_1) \prod_{k=1}^{n-1} g_{k+1}(y_{k+1}, m_{k+1} | \vec{P}_k), \quad (7)$$

with $\vec{P}_n = (y_1, m_1), \dots, (y_n, m_n)$. We assume that the marks are independently distributed as follows:

$$g_{n+1}(y_{k+1}, m_{k+1} | \vec{P}_n) = f_{n+1}(y_{n+1} | \vec{P}_n) M(y_{n+1}),$$

where $M(y_{n+1}) = B(\alpha(y_{n+1}))$, $\alpha$ being the true prevalence map.

We take into account the characteristics mentioned above through the conditional densities:

$$f_{n+1}(u | \vec{P}_n) \propto \alpha_r(u) k(y_n, u) \pi(\iota_{\vec{P}_n^*}(u)), \quad (8)$$

where $\alpha_r$ is a risk map, $k$ is a proposal kernel, $\pi$ is a re-weighting factor and $\iota_{\vec{P}_n^*}$ is the information map derived from the first $n$ points. We define these three components as follows.

**Information map** The map $\iota_{\vec{P}_n^*}$ gathers the information from the initial $n$ samples. The formula for $\iota_{\vec{P}_n^*}$ is given by:

$$\iota_{\vec{P}_n^*}(u) = \frac{1}{2} + \frac{\sum_{i=1}^{n} (2m_i - 1) \exp\left(\frac{-\|y_i - u\|^2}{2h^2}\right)}{2 \sum_{i=1}^{n} \exp\left(\frac{-\|y_i - u\|^2}{2h^2}\right)}$$

where $h$ is a smoothing parameter describing how the knowledge of the sanitary status of an individual propagates spatially.

<footnote n="4">⁴ A random sequence $(y_1, \dots, y_n)$ is said to be Markovian if the conditional distribution of any point, given all the preceding points, depends only on the immediately preceding point. In other words, for all $i$, the distribution of $y_i$ given $(y_1, \dots, y_{i-1})$ is the same as the distribution of $y_i$ given $y_{i-1}$.</footnote>
**Re-weighting factor** The re-weighting factor $\pi$ is defined as being proportional to:

$$ \bar{\pi}(p) = p^a(1-p)^b, a > -1, b > -1 $$

More precisely, $\bar{\pi}$ is normalized so that its maximum value is 1. The re-weighting factor allows penalizing or favoring areas based on the information they bring. Its shape allows for a wide range of behaviors depending on the values of parameters a and b:

* When a is significantly larger than b, the function $\pi$ prioritizes regions where the information indicates disease presence.

* When b is significantly larger than a, $\pi$ prioritizes regions where the information suggests an absence of the disease.

* When $a = b$, $\pi$ concentrates on regions where no individuals have been sampled.

* When $-1 < a, b < 0$, $\pi$ balances the attention between regions likely to be infected and those likely to be uninfected.

**Kernel** We choose the following isotropic$^5$ proposal kernel:

$$ k(y_i, u) = k(\|y_i - u\|) = \exp\left(\frac{-\|y_i - u\|^2}{2\sigma^2}\right), $$

where the parameter $\sigma^2$ of this kernel determines the distribution of distances between subsequent points.

By estimating model parameters from observed surveillance data, we can describe the underlying decision-making process governing the sampling strategy (See 6.3). Conversely, by adjusting these parameters, the model can be leveraged to extend usual survey design models and to generate efficient sampling schemes tailored to specific surveillance objectives (See 6.2).

**Remark 5.1.** The shape of the sequential density $f_{n+1}$ can be generalized to other contexts or surveillance problems:

* The first term $\alpha_r$ would remain the intensity or inherent risk of the disease, reflecting a preference for sampling in higher-risk regions.

* The second term $k$ would encapsulate the information summarized by all the previous samples in $\vec{P}_n$, representing the self-interaction component of the density.

* The third term $\pi$ would adjust the density to account for unsampled points, representing the expected information missed due to incomplete sampling.

Such a formulation could serve as a convenient surrogate model, for instance in cases where the marked point process of all individuals is a Gibbs point process.

## 5.2. Simulation method

Let $\alpha_r: W \to [0, 1]$ be a risk map intended to approximate the prevalence map, either computed using risk factors or derived from other surveillance schemes. The risk map $\alpha_r$ can be dynamically updated throughout the surveillance campaign. For instance, it may incorporate prevalence estimates derived from previously collected samples, in which case $\alpha_r(u) \approx \alpha_{\vec{P}_n}(u)$. Our model essentially involves proposing points sequentially as follows:

* The first point $y_1$ is randomly simulated, with its density being proportional to $\alpha_r$, for example. Its mark is revealed from the unknown prevalence map $\alpha: m_1 \sim B(\alpha(y_1))$.

* Given the points $\vec{P}_n = (P_1, \dots, P_n) = ((y_1, m_1), \dots, (y_n, m_n))$, the next point $y_{n+1}$ has the following density:

$$ f_{n+1}(u|\vec{P}_n) \propto \alpha_r(u)k(y_n, u)\pi(\vec{P}_n(u)) $$

This is equivalent to proposing the next point with a density proportional to $\alpha_r(u)k(y_n, u)$ and accepting it with probability $\pi(\vec{P}_n(u))$. Given $y_{n+1}$, the mark of the $(n+1)$th sample is revealed from the unknown prevalence map: $m_{n+1} \sim B(\alpha(y_{n+1}))$.

## 5.3. Parameter estimation

Given a sequential surveillance scheme $\vec{P}_n = (y_i, m_i)_{1 \le i \le n}$ supposedly originating from the SSS with unknown parameters $\theta_0 = (h_0^2, a_0, b_0, \sigma_0^2)$, we want to estimate these parameters. The chosen method for this estimation is the Maximum Likelihood approach. We have the following expression for the log-likelihood:

$$ l(\theta) = \sum_{i=1}^{n-1} (1 + m_{i+1}) \log(\alpha(y_{i+1})) + \sum_{i=1}^{n-1} \log k(y_i, y_{i+1}) + \sum_{i=1}^{n-1} \log \pi(\vec{P}_i(y_{i+1})) \\ + \sum_{i=1}^{n-1} (1 - m_{i+1}) \log(1 - \alpha(y_{i+1})) - \sum_{i=1}^{n-1} \int_W \alpha_r(u)k(y_i, u)\pi(\vec{P}_i(u))du $$

<footnote n="5">5 An isotropic kernel is such that $k(u, v)$ only depends on the distance $\|u - v\|$ for all u and v in W.</footnote>
$$ 
\begin{aligned}
l(\theta) = & \sum_{i=1}^{n-1} - \frac{\|y_{i+1} - y_i\|^2}{2\sigma^2} + \sum_{i=1}^{n-1} a \log(\iota_{\vec{P}_i}(y_{i+1})) + b \log(1 - \iota_{\vec{P}_i}(y_{i+1})) \\
& - \sum_{i=1}^{n-1} \log \int_W \alpha_r(u) \exp(-\frac{\|u - y_i\|^2}{2\sigma^2}) \iota_{\vec{P}_i}(u)^a (1 - \iota_{\vec{P}_i}(u))^b dx
\end{aligned} 
\quad (9) 
$$

The normalizing integral depends on the model parameters and will be evaluated based on Riemann sum method.

## 5.4. Parameter optimization

In this section, we investigate the conditions under which our model minimizes the Integrated Bernoulli Variance (IBV). To this end, we propose specific parameter values (a, b) for the Sequential Surveillance Scheme that lead to an effective reduction of the IBV, thereby improving the accuracy of prevalence estimation. We do not optimize the parameters n, h, and σ, as we adopt a hierarchical approach in which the selection of a and b precedes the determination of the other parameters. Once a and b are fixed, the sample size n and the smoothing parameters h and σ are chosen accordingly. Here, we focus exclusively on selecting the pair (a, b).

Let us recall the conditional densities of our process, given a risk map α_r and a sequence of points $\vec{P}_n$:

$$ f_{n+1}(u|\vec{P}_n) = \alpha_r(u)k(y_n, u)\pi(\iota_{\vec{P}_n}(u)) $$

where k denotes the kernel, π is the re-weighting factor and $\iota_n$ is the information map derived from the first n samples.

Among all possible re-weighting factors π : [0, 1] ↦→ [0, 1], the optimal one π* can be identified as the one that minimizes the IBV.

**Conjecture 2.** The optimal re-weighting factor π : [0, 1] ↦→ [0, 1] for minimizing the IBV is π*(p) = 1/2(p).

**Reasoning 2.** For S the set of surveillance schemes of size n with at least one infected individual and one non-infected, $P_n \in S$ and $\iota_{P_n}(u) = \frac{1}{2} + \frac{\sum_{i=1}^n K(y_i, u)(2m_i - 1)}{2 \sum_{j=1}^n K(y_j, u)} = \frac{1}{2} + \frac{\sum_{i=1}^n \omega_i(u)(2m_i - 1)}{2 \sum_{j=1}^n K(y_j, u)}$, we write the objective of IBV minimization (4), which is minimizing the integral of:

$$ 
\begin{aligned}
\iota_{P_n}(u)(1 - \iota_{P_n}(u)) &= \left(\frac{1}{2} + \frac{1}{2}\sum_{i=1}^n \omega_i(u)(2m_i - 1)\right)\left(\frac{1}{2} - \frac{1}{2}\sum_{i=1}^n \omega_i(u)(2m_i - 1)\right) \\
&= \frac{1}{4} - \frac{1}{4}\left(\sum_{i=1}^n \omega_i(u)(2m_i - 1)\right)^2.
\end{aligned}
$$

Therefore, minimizing the integral $\int_W \iota_{P_n}(u)(1 - \iota_{P_n}(u))du$ is similar as maximizing $\int_W (\sum_{i=1}^n \omega_i(u)(2m_i - 1))^2 du$. Developing, this leads to maximize:

$$ \int_W \left( \sum_{i=1}^{n} \omega_i(u)^2 - 2 \sum_{\substack{i \neq j \\ m_i \neq m_j}} \omega_i(u)\omega_j(u) + \sum_{\substack{i \neq j \\ m_i = m_j}} \omega_i(u)\omega_j(u) \right) du. \quad (10) $$

In the absence of knowledge about the marks, the optimal strategy is to sample surveillance schemes such that $\omega_i(u)\omega_j(u) \approx 0$, meaning that the samples must be sufficiently spaced apart. This approach corresponds to sampling in regions where the information is close to 1/2, yielding the optimal reweighting factor $\pi^*$. ☐

This optimal model is impractical to simulate because, once both infected and non-infected individuals have been sampled, the subset of W where the information equals 1/2 has a Lebesgue measure of 0. Consequently, rejection sampling cannot be employed to simulate this model. Hence, a compromise between numerical feasibility and optimality is necessary when choosing parameters (a, b). Since π is proportional to a beta distribution, setting (a, b) = (m, m) with m ∈ ℝ concentrate the mass of π on 1/2. Larger value of m yield better performance, whereas smaller values of m facilitate simulation.

**Example 5.1.** To illustrate the behavior of our model for this choice of parameters, we consider the following simple example in the unit square:

* $\iota_{i+1}(u|\vec{P}_i) \propto \pi(\iota_{\vec{P}_i}(u))$, where a = b = 20. In this scenario, only the information from previous samples influences future sampling locations.

* $\alpha(u) = u_1$, with $u = (u_1, u_2) \in [0, 1]^2$.

We simulate 50 samples from this model, depicted in Fig. 1. As shown in Fig. 1(b), the space is initially explored completely according to the kernel of the information map. Then, in Fig. 1(c), the information gathered during this exploration is used to reduce the uncertainty of the information map.
<figure><img src="image_5.png" /><figcaption>Fig. 1. Plot of the samples generated from the sequential densities of the Example 5.1. Red points represent infected samples, while green points represent non-infected samples. The background image represents the information map of the samples in the last two plots, while in the first plot, it depicts the true prevalence.</figcaption></figure>

## 5.5. Summary statistics

The abilities of the SSS to adopt a large variety of behaviors can be analyzed through new summary statistics, adapted to the sequential nature of the data. Time-dependent summary statistics can be used to describe the evolution of the sampling process: the Cumulative Lagged Positive Proportion (CLPP) that is similar to a cross K-function between all points and points marked by 1, the Ball Union Coverage of Positive Points (BUC) that shows at what speed the space is explored, the Cumulative Scanpath Length (CSL) that shows the distance between a point and the previous point marked by 1, and the Cumulative Recurrence (CR) that is similar to the K-function of the unmarked point process. Let $\vec{P}_k$ be a surveillance scheme, and $r > 0$.

* The Cumulative Lagged Positive Proportion (CLPP) is the cumulative sum, for $i$ ranging from 2 to $k$, of the number of previous points in the sequence that are positive and located within a distance smaller than $r$ from $y_i$:

$$ \text{CLPP}(\vec{P}_k) = \sum_{i=2}^{k} \sum_{\substack{j \le i-1 \\ m_j=1}} \mathbb{1}_{B(y_i, r)}(y_j). $$

This quantity is expected to vary depending on the parameters of the rebalancing function: a high $\pi$ for low proportions is expected to yield a low CLPP, while a high $\pi$ concentrated on high proportions is expected to yield a high CLPP.

* The Ball Union Coverage of Positive Points (BUC) is the normalized area covered by positive points:

$$ \text{BUC}(\vec{P}_k) = \frac{|\bigcup_{i=1}^{k} m_i \cap W|}{|W|}. $$

Different behaviors are expected depending on whether the focus is on positive individuals or not.

* The Cumulative Scanpath Length (CSL) represents the cumulative sum of distances from each point to its nearest preceding positive point in a temporal sequence, following the broken line created by successive points:

$$ \text{CSL}(\vec{P}_k) = \sum_{i=1}^{k} \sum_{j=0}^{M_i} \|y_{i+1-j} - y_{k-j}\|, $$

where $M_i = \inf\{n \in \mathbb{N}^* | m_{i+1-n} = 1\} - 1$. This quantity accounts for both the distances between successive points, reflecting the proposal kernel, and the positivity rate emphasized by the rebalancing function.

* The Cumulative Recurrence (CR) is the cumulative count of previous samples that are within a distance $r$. For sample $k$ and $r > 0$, it is given by:

$$ \text{CR}(\vec{P}_k) = \sum_{i=2}^{k} \sum_{j=1}^{l-1} \mathbb{1}_{y_j \in B(y_i, r)}. $$

This quantity indicates whether the process favors points near previous positions.

These statistics emphasize different aspects of the model: CR focuses on the proposal kernel, BUC focuses on the re-balancing function $\pi$, and CLPP and CSL on both. For similar statistics defined on sequential point processes, see Penttinen and Ylitalo (2016), Yazigi et al. (2021).

## 6. Applications

To validate and demonstrate the practical utility of our proposed Sequential Surveillance Scheme (SSS), we apply it in both simulated environments and on real-world surveillance data. These applications allow us to assess the model's performance in guiding efficient sampling strategies, particularly for optimizing prevalence estimation and exploring its adaptability to complex epidemiological contexts.
<figure><img src="image_6.png" /><figcaption>Fig. 2. Corresponding summary statistics. Pay attention to the scale when comparing the graphs of the (CLPP).</figcaption></figure>

## 6.1. Simulation study

In the following, we set the parameters as: W is the square [0, 10] × [0, 10], $h^2 = 1$, $\sigma^2 = 1$, and $\alpha_r$ is a mixture of two Gaussian densities with means (3, 3) and (9,9) and covariance matrices $0.1I_2$ and $0.05I_2$, where $I_2$ is the $2 \times 2$ unit diagonal matrix.

We analyze our model with four different parameter pairs $(a, b)$:

- (M1) : (1, 1),

- (M2) : (1, 200),

- (M3) : (0, 0),

- (M4) : (0, 100).

The parameters $(a, b)$ guide the proposal of new points based on the information given by previous points in the space. Model (M3) corresponds to the case where the information of an area has no influence on the acceptance of new points. Models (M2) and (M4) focus on areas where the information suggests low prevalence. We present the summary statistics corresponding to each of these realizations.

The CLPP conveys the positivity of an area: a rapid increase of the CLPP indicates a high number of positives in the sampled area, whereas a slower increase suggests fewer positives. (M1) concludes with the highest CLPP, as illustrated in Fig. 2(a). The BUC for the parameter set (M1) ends up lower than the other BUC values, indicating that process (M1) has explored the space less. The CSL represents the spatial exploration: a steep slope indicates longer distances between each positive point, while a linear slope reflects regularly spaced positive points across space and time. Finally, the CR corresponds to the CLPP without considering the marks, and similar behaviors are observed.

To highlight differences between parameter sets, we simulate each of the processes (M2), (M3), and (M4) 20 times and compare their summary statistics.

Models (M2) and (M4), represented by the blue and black curves in the graphs, respectively, focus on low prevalence areas suggested by the information map, and exhibit greater distances between consecutive positive points (see CSL in Fig. 3).
<figure><img src="image_7.png" /><figcaption>Fig. 3. Summary statistics of 20 realizations for each process, denoted as (M2) for black curves, (M3) for red curves, and (M4) for blue curves.</figcaption></figure>

## 6.2. Comparison between different surveillance models

To gain insights into the advantage of using the proposed SSS and into the influence of parameters *a* and *b*, we compare the SSS with multiple sets of parameters (*a*, *b*) and alternative surveillance methods when observations are drawn from various prevalence maps displayed in Fig. 4.

The competing surveillance methods are listed in Fig. 5. They all consist in drawing 500 points in total in the spatial domain and the corresponding infection marks conditional on the prevalence map. In addition to 9 SSS preliminary initiated with 50 samples independently and uniformly drawn in the spatial domain and corresponding to different parameter vectors (*a*, *b*), we consider 3 risk-based methods consisting in preliminary drawing *q* samples independently and uniformly in the spatial domain and independently drawing the following samples from the density $f_{n+1}(u) \propto \alpha_{P_n}(u)$, with $q = 100, 200$ or $250$ and $n = q$. We also consider the completely naive strategy consisting in drawing the 500 points uniformly and independently in the spatial domain.

Regarding the specifications of the SSS model defined by the sequential densities $f_{n+1}(u|\hat{P}_n) \propto \alpha_{P_n}(u)\pi(\iota_{P_n}(u))$, by choosing $a = -0.5$ and $b = 0.5$, we obtain $f_{n+1}(u|\hat{P}_n) \propto \sqrt{\alpha_{P_n}(u)(1 - \iota_{P_n}(u))}$. Given observations from the first uniform survey, we have $\iota_{P_n}(u) \approx \alpha_{P_n}(u)$, and this choice of density concentrates on areas of uncertainty regarding the disease. Alternatively, setting $a = b = 0.5$ leads to $f_{n+1}(u|\hat{P}_n) \approx \alpha_{P_n}(u)\sqrt{\alpha_{P_n}(u)(1 - \alpha_{P_n}(u))}$, thereby giving more weight to uncertain high-risk regions. Preliminary simulations with high and similar values for *a* and *b* indicate that these configurations do not outperform the parameter sets specified in Fig. 5. While high values of $a = b$ efficiently reduce IBV, as stated in 2, the relation 1 suggests that a trade-off between IBV reduction and sampling in high-risk regions is necessary.

For each prevalence map presented in Fig. 4, we simulated 400 surveillance schemes, each consisting of 500 sampling points, based on the models listed in Fig. 5. We then computed the $L_2$ norm error between the true prevalence map and its estimation at regular intervals in terms of sampling size, specifically every 25 samples, starting from the 50th sample. Results are presented in Fig. 5.

For a fixed level of precision in prevalence estimation, our model needs far less samples to achieve this precision than uniform or risk-based sampling, as seen by the lengths of the horizontal gaps between the curves. Moreover, our model outperforms other models in the low-sample regime; however, all models appear to converge toward similar estimation accuracy as the sample size increases. This observation may be explained by the fact that our model provides a more precise estimation of prevalence in regions where the disease has been detected, but a less accurate estimation in regions where the disease is present but undetected. Consequently, in certain scenarios, the disease may remain undetected in specific regions despite its presence. This suggests that while our model maintains high precision in detected regions, its performance is less reliable in undetected regions, whereas a uniform sampling
<figure><img src="image_8.png" /><figcaption>**Fig. 4.** Prevalence maps on the square (0,1) used in the comparison of various sampling designs. The first three maps are similar except that the prevalence in the bottom left corner is either low, medium or high. The next two maps display one disease focus (hotspot) and three disease foci (multi-Hotspot), respectively. The last map is characterized by a rapid transition between low and high prevalence, following a wavelet.</figcaption></figure>

strategy ensures eventual coverage across all areas. This could explain the observed behavior at larger sample sizes. Introducing a fraction of randomly selected uniform samples at certain intervals in areas estimated as non infected could potentially mitigate this limitation.

The results also indicate that the sharper the transition from low to high prevalence, the better our model performs. In other words, if the spatial boundary of the disease is narrow, our model requires fewer samples than other models to estimate this boundary with the same level of precision.

Regarding the performance across different parameter choices within our model, similar behaviors are observed for the displayed parameter sets.

We now provide a summary of the different sampling models and their appropriate use cases.

* **Uniform sampling** consists in selecting sampling locations independently and identically distributed over the spatial domain. This approach ensures an unbiased estimate of the global prevalence function.

**Advantages:** It is simple to implement and does not require prior knowledge about the spatial distribution of the prevalence.

**Disadvantages:** It is inefficient when the spatial domain is large, as the number of samples required for accurate estimation grows rapidly and often unapplicable in practice. When the prevalence function is equal to zero over a large subset of the domain (i.e., the Lebesgue measure of the non-infected regions is large), uniform sampling allocates unnecessary resources to these areas.

* **Two-Phase Sampling Based on Risk** involves an initial uniform sampling phase followed by a second phase where additional samples are collected with density proportional to the level of prevalence identified in the first phase.

**Advantages:** It is more efficient than pure uniform sampling, as it concentrates effort on regions where the prevalence is positive and thus reduces the number of samples wasted in non-infected regions.

**Disadvantages:** It poorly differentiates between highly infected and marginally infected areas, leading to potential inefficiencies in sample allocation. Also, the transition sample size from the first to the second phase requires a predefined threshold, which may not always be optimal. One possible approach to determining the uniform sample size is to choose a sample size such that all infected regions are detected with high probability under uniform sampling.

* Our proposed method aims to minimize uncertainty in the estimation of the prevalence function. The sampling locations are chosen adaptively based on the estimated prevalence function, with a focus on areas where uncertainty is highest.

**Advantages:** It is more efficient than uniform sampling when the prevalence function exhibits sharp spatial transitions and particularly advantageous when the Lebesgue measure of the transition boundaries of the prevalence function is small, as it
<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Comparison of mean *L*² prevalence error for different prevalence maps and different sampling models. The mean error is computed for increasing sampling sizes based on batches of 25 samples collected starting from the 50th sample.</figcaption></figure>

concentrates samples in these uncertain regions. Also, it avoids excessive sampling in non-infected regions, making it more efficient than uniform and two-phase sampling approaches.

**Disadvantages:** The adaptive nature of the method may introduce implementation complexity compared to simpler approaches.

If an infected area is not detected initially, it is more likely to remain undetected compared to uniform sampling.

In Fig. 6, we show as example one realization of each surveillance method and the corresponding estimation of the prevalence. More examples are shown in Appendix B.
<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Examples of 125 samples surveys and corresponding prevalence estimation for the models uniform, two-phase based on risk and SSS with parameters *a* = 0 and *b* = 0.5 from left to right and top to bottom. The prevalence map is “transition”.</figcaption></figure>

## 6.3. Surveillance of the plant pathogen Xylella fastidiosa

We now estimate the parameters of our model on real surveillance data with the MLE method 5.3 and evaluate the fit using the new summary statistics. We also search for a more efficient parameter set to improve sampling in terms of prevalence estimation.

*Epidemiological context.* *Xylella fastidiosa* (*Xf*) is a phytopathogenic bacterium that infects over 600 plant species, including economically significant crops like olives, almonds, and citrus. It spreads through xylem-feeding insects, such as leafhoppers, and colonizes the plant's xylem vessels, forming biofilms that block sap flow. This blockage causes symptoms like leaf scorch, wilting, and plant death in severe cases.

First identified in Europe in 2013 in Italy, *Xf*, subspecies *pauca*, caused a significant epidemic in Puglia's olive groves. France has experienced multiple introductions: *Xf* was detected in 2015 in Corsica and Provence-Alpes-Côte d'Azur (PACA), with *Xf multiplex* being the primary subspecies and isolated cases of *Xf pauca* were reported in PACA olive trees in Menton in 2019. These cases prompted targeted surveillance and plant removal measures to prevent further spread.

The wide host range of *Xf multiplex* complicates management, particularly due to its impact on ornamental plants like *Polygala myrtifolia* in urban settings. Surveillance strategies in France have evolved to include risk-based approaches, utilizing risk maps to prioritize high-risk areas for intensive monitoring. Recent outbreaks in Occitania in 2020 and 2021 underscore the ongoing threat posed by *Xf* across the Mediterranean basin, despite efforts to manage the disease through plant removal, vector control, and intensified monitoring.

In PACA, initial measures focused on urban and peri-urban zones with systematic plant removal in *infected zones* (100-meter radius around infected plants) and *buffer zones* (4-kilometer radius). Infected zones required the removal of all host plants, while buffer zones involved intensive surveillance and monitoring. The detection of *Xf pauca* in 2019 led to localized actions, including biannual inspections and preservation of heritage olive trees using insect-proof nets. With the introduction of European Regulation, management shifted to a more flexible, risk-based approach, balancing eradication efforts with plant preservation. In the following, we will focus on the surveillance of *Xf multiplex* in this region. More details on the biology, the epidemiology and the European history of *Xf* are provided in Sicard et al. (2018), Jacques et al. (2023).

In Corsica, strategies adapted to the region's ecological and cultural context, prioritizing actions based on risk assessments. High-value and heritage plants are preserved, and decision-making incorporates biological, economic, and cultural criteria. The region's surveillance focuses on sensitive areas hosting diverse *Xf multiplex* populations. Here, we consider a machine learning algorithm to estimate the risk map $\alpha_r$ from these observations and various risk factors (as described in Portes et al., 2024 and available from Portes et al., 2024a).

**Fig. 7** illustrates the two regions of interest and the 40,000 surveillance samples (sampled plants are shown in green if *Xylella fastidiosa* is absent, and in red otherwise) collected between 2015 and 2021.

*Estimation of parameters based on the real sampling.* As mentioned, the surveillance strategy in France initially consisted of uniform sampling across the territory, with intensified sampling in the vicinity of a detected positive case. However, in practice, due to land constraints, obstacles, and logistical considerations, surveillance practitioners had to adapt their strategy locally. Our objective is to infer the actual surveillance strategy implemented in practice by determining the parameters of our model that maximize the likelihood of generating the observed surveillance data. Since our model encompasses a wide range of possible behaviors, we optimistically assume that it can approximate the true surveillance strategy with reasonable accuracy.

We estimate the model parameters ($h^2$, *a*, *b*, $\sigma^2$) using the maximum likelihood (ML) method. The estimation is based on 1,000 sequential samples from the PACA region, chosen such that the total prevalence over these samples is the closest to the total prevalence observed in the entire dataset. Here, we fix $h^2 = \sigma^2$ because we assume that the proposal of new surveillance points should operate on the same scale as the information available. The algorithm BOBYQA (Powell, 2009) from the package NLopt (Johnson, 2007) is used to maximize the log-likelihood. The estimated values for the kernel parameters are $\hat{h} = \hat{\sigma} = 15.04$ kilometers (with standard error equal to 0.34 km), while the estimates for the reweighting factors are $\hat{a} = 0.31$ and $\hat{b} = 0.063$ (with standard errors equal to 0.018 and 0.002, respectively).
<figure><img src="image_11.png" /><figcaption>Fig. 7. Real surveillance scheme in the south of France. Red points indicate infected samples while green points represent non-infected samples.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 8. Risk map ($\alpha_r$) derived from environmental and bioclimatic factors.</figcaption></figure>

Data sampled in Corsica are used to estimate the risk map $\alpha_r$ using risk factors, including land types, soil compositions, and climatic conditions. Fig. 8 illustrates $\alpha_r$ obtained from XGBoost (eXtreme Gradient Boosting).

Assessing the goodness-of-fit of the models is performed by comparing the summary statistics (introduced in Section 5.5), which are calculated from the data and from 99 simulated realizations of the fitted models (see Fig. 10). The BUC and SLC statistics demonstrate a good fit for our model, while the CLPP and CR indicators exhibit a poor fit. Specifically, the CLPP and CR statistics show higher values for the real data compared to the simulated ones. The CLPP indicates the presence of points near previous positive points, while the CR captures the presence of points near previous points in general. Multiple reasons can explain this behavior:

* The window considered in practice is likely smaller than the one used in our model, resulting in faster spatial exploration. Moreover, the surveillance window may not be fixed and, instead, may evolve with time. This proposal can be supported by the result of Soubeyrand et al. (2018), which shows from the same data that the preference in sampling at-risk hosts varied across time. An adaptive window approach could be implemented as follows: starting with an initial window $W_1$ contained within the larger surveillance region $W$, new sample points proposed outside $W_1$ are accepted with a certain probability $p$. Accepted samples then expand the initial window $W_1$ into a larger window $W_2$, and this process can be further iterated.
<figure><img src="image_13.png" /><figcaption>Fig. 9. Re-weighting factor estimated from the data.</figcaption></figure>

* In practice, surveillance is likely conducted at multiple scales. Initially, surveillance may occur at the regional scale, with refinement to a finer scale upon the detection of an infected plant. This multi-scale approach can be incorporated into our model by introducing a mixture of kernels. The transition between kernels could be governed either deterministically or stochastically, for example using a Markov chain.

* The risk map employed in our model, calibrated using surveillance data from Corsica, may be overly conservative. In practice, surveillance practitioners are more likely to begin with a uniform risk map, whereas our model's risk map identifies significantly fewer regions as being at risk compared to a uniform distribution.

This is further supported by two facts: the first is that the Cumulative Scanpath Length exhibits a good fit, meaning that the smoothing parameters $\sigma^2 = h^2$ are well chosen. The second is that the Ball Union Coverage of positive points also exhibits a good fit, which means that the parameters a and b are well selected, ensuring there is neither over-selection nor under-selection of positive points by our fitted model.

*Comparison with a more efficient surveillance scheme.* Now that we have gained insight into the surveillance strategy implemented in practice, our next objective is to explore potential improvements, potentially with only slight modifications to our model's parameters. As we can see in Fig. 9, the estimated re-weighting factor for the real surveillance, $\pi_{\text{data}}$, differs from re-weighting factors we found to be efficient. $\pi_{\text{data}}$ focuses more on infected regions than others that favor exploration.

In practice, when a positive sample is detected, additional samples are often taken in the surrounding area to assess the local prevalence. This suggests that the surveillance process operates on multiple scales: it explores large distances while exploiting information from very localized regions. Moreover, not all surveillance samples might serve the same purpose in practice. Surveillance can simultaneously target different objectives such as prevalence estimation, early detection in a new area, and others.

Fig. 11 displays simulations of 1,000 samples using different parameter sets, to compare $\pi_{\text{data}}$ with other re-weighting factors in terms of sampling distribution. For these simulations, all other parameters are fixed to their estimated values from the real data.

As observed, using $\pi_{\text{data}}$, the model tends to sample extensively within infected regions. As the parameters a = b increase, the model exhibits a greater tendency to explore new areas.

In Appendix A, we compare the efficiency of the inferred surveillance process with that of a model that differs only slightly from the estimated one.

## 7. Discussions

We propose a sequential marked point process model to both describe and optimize surveillance surveys for plant diseases. This model integrates a risk map, $\alpha_r$, as the baseline intensity of the surveillance point process to represent spatial risk, and an information map, $\iota$, as an interaction function to capture the knowledge gained in previously explored areas. Directions for generalizing the model are also discussed, including a formulation of sequential densities (Remark 5.1) for other contexts and the potential to account for multiple surveillance scales (e.g., regions, forests, etc.).

The model parameters are inferred using real surveillance data by maximizing the corresponding likelihood. By employing novel summary statistics for sequential marked point processes, we observe that real surveillance exhibits a stronger tendency than our model to revisit high-risk regions. This suggests that inspectors often oversample certain areas rather than exploring previously unvisited regions. Additionally, we identify model parameters that are efficient for prevalence estimation. An effective strategy is a compromise between maximizing the amount of information, as quantified by the IBV, and sampling in high-risk regions.
<figure><img src="image_14.png" /><figcaption>(a) BUC</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) SLC</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(c) CLPP</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(d) CR</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(a) 1,000 samples generated from the estimated parameters on real data.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b) 1,000 samples generated with $a = b = 2$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(c) 1,000 samples generated with $a = b = 4$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(d) 1,000 samples generated with $a = b = 5$.</figcaption></figure>

<caption>Fig. 11. Sets of 1,000 samples obtained from different set of parameters (a, b). The other parameters are fixed at the values estimated from the data. The two first points of the samples are the same for all simulations. Red points are infected individuals and green points are non infected individuals.</caption>
Further work may focus on extending the model by incorporating disease propagation dynamics (Abboud et al., 2019, 2023) and risk factors (Martinetti and Soubeyrand, 2019; Cendoya et al., 2020; Portes et al., 2024) into the self-interaction functions. One may also introduce into the model temporal variations in the parameter values for improving model flexibility and, therefore, goodness-of-fit. It has to be noted that such model flexibility could be handled in the estimation stage with a data set such as the one we used. Indeed, we only used one thousand observations for estimating the fixed parameters, but the data set includes several tens of thousands of observations.

Using a parametric sequential surveillance model is specifically valuable for characterizing a real sampling process (estimation stage), and evaluate how it could be improved for better assessing prevalence across space (optimization stage; by the way this stage could be strengthened by optimizing the other parameters). However, further work could be carried out to optimize surveillance schemes using a model-free approach, where no specific assumptions are made regarding the prevalence estimate or the density of the surveillance process. As disease surveillance often seeks to optimize multiple objectives in practice (e.g., prevalence estimation, early disease detection), we will explore modeling frameworks for multi-objective surveillance.

Finally, in this paper, the sample size *n* of the surveillance schemes is fixed. Future work will focus on designing the sample size *n*, for example, to ensure sufficient estimation precision.

## Acknowledgments

This work was supported by the BEYOND project (grant ANR-20-PCPA-0002) and the HORIZON BeXyl Project (grant 10106 0593). We thank the ESV Platform (https://plateforme-esv.fr/) for providing *Xylella fastidiosa* data.

## Appendix A. Comparing the estimated and optimal surveillance scheme

In practice, the disease within surveillance regions can be completely absent, completely present, or fluctuate. To illustrate these possible surveillance behaviors, we use the true prevalence map shown in Fig. A.12 within the unit square.

We compare 400 realizations of 500 samples generated from six models based on the sequential densities:

$$f_{n+1}(u|\hat{P}_n) \propto \alpha_{\hat{P}_n}(u)\pi(t_{\hat{P}_n}(u)).$$

For the first model, $M_{data}$, we use the parameters $a = 0.31$ and $b = 0.063$, estimated from the real surveillance data. For the second model, $M_2$, we set $a = b = 0.18$, which corresponds to a very slight modification of the estimated parameters. We also use four other parameter sets of our model, found to be efficient by our simulations in 6.2. The parameters are shown in the legend of Fig. A.13, in which the error $\int (\alpha - \alpha_{\hat{P}_n})^2$ is plotted for varying sample sizes *n*. For the final sample size, fewer than half the samples are required to achieve the same precision with the efficient parameters compared to the estimated ones. The results show that the gap in prevalence error between $M_{data}$ and $M_2$ increases as the sample size grows and $M_2$ exhibits better performances for the spatial error. This means that a very little adjustment to the current surveillance strategy could benefit greatly to the surveillance efficiency.

## Appendix B. More examples of the 3 survey models on other prevalence maps

In the next two figures, we show as examples one realization of each surveillance method and the corresponding estimation of the prevalence (see Figs. B.14 and B.15).

<figure><img src="image_22.png" /><figcaption>Fig. A.12. True prevalence map in the unit square.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. A.13. Mean prevalence error for different sample sizes and different parameter sets of our model.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. B.14. Examples of 125 samples surveys and corresponding prevalence estimation for the models uniform, two-phase based on risk and SSS with parameters $a = 0$ and $b = 0.5$ from left to right. The prevalence map is “hotspot”.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. B.15. Examples of 125 samples surveys and corresponding prevalence estimation for the models uniform, two-phase based on risk and SSS with parameters $a = 0$ and $b = 0.5$ from left to right. The prevalence map is “Multi Hotspot”.</figcaption></figure>

## Data availability

The codes for simulating the Sequential Surveillance Scheme (SSS) and estimating the model parameters from real data are available upon request.
References

Abboud, C., Bonnefon, O., Parent, E., Soubeyrand, S., 2019. Dating and localizing an invasion from post-introduction data and a coupled reaction-diffusion-absorption model. J. Math. Biol. 79, 765-789.

Abboud, C., Parent, E., Bonnefon, O., Soubeyrand, S., 2023. Forecasting pathogen dynamics with Bayesian model-averaging: Application to *Xylella fastidiosa*. Bull. Math. Biol. 85, 67.

Anyosa, S., Eidsvik, J., Pizarro, O., 2023. Adaptive spatial designs minimizing the integrated Bernoulli variance in spatial logistic regression models with an application to benthic habitat mapping. Comput. Statist. Data Anal. 179, 107643.

Atkinson, A., Donev, A., Tobias, R., 2007. Optimum Experimental Designs, with SAS, vol. 34, OUP, Oxford.

Cameron, A., 2012. The consequences of risk-based surveillance: Developing output-based standards for surveillance to demonstrate freedom from disease. Prev. Vet. Med. 105 (4), 280-286.

Cendoya, M., Martínez-Minaya, J., Dalmau, V., Ferrer, A., Saponari, M., Conesa, D., López-Quílez, A., Vicent, A., 2020. Spatial Bayesian modeling applied to the surveys of *Xylella fastidiosa* in Alicante (Spain) and Apulia (Italy). Front. Plant Sci. 11, 1204.

Chiu, S.N., Stoyan, D., Kendall, W.S., Mecke, J., 2013. Stochastic Geometry and Its Applications. John Wiley & Sons, Chichester.

Daley, D.J., Vere-Jones, D., 2003. An Introduction to the Theory of Point Processes: Volume I: Elementary Theory and Methods. Springer, New York.

European Food Safety Authority (EFSA), 2012. A framework to substantiate absence of disease: the risk based estimate of system sensitivity tool (RiBESS) using data collated according to the EFSA standard sample description - An example on *Echinococcus multilocularis*. EFSA Support. Publ. 9, 366E.

Fuller, K.B., McIntosh, C., Zidack, N., 2020. Valuing disease prevention in a vegetatively propagated annual crop: Benefits from the Montana Seed Potato Certification Program. Plant Dis. 104, 2060-2067.

Ge, Y., Eidsvik, J., Mo-Bjørkelund, T., 2023. 3-D adaptive AUV sampling for classification of water masses. IEEE J. Ocean. Eng. 48, 626-639.

Illian, J., Penttinen, A., Stoyan, H., Stoyan, D., 2008. Statistical Analysis and Modelling of Spatial Point Patterns. John Wiley & Sons, Chichester.

Jacques, M.-A., Soubeyrand, S., Dupas, E., 2023. Décrypter l'émergence de *Xylella fastidiosa* en Europe et surveiller son expansion. In: Lannou, C., Rasplus, J.-Y., Soubeyrand, S., Gautier, M., Rossi, J.-P. (Eds.), Crises Sanitaires En Agriculture: Les Espèces Invasives Sous Surveillance. Quæ, Versailles.

Johnson, S.G., 2007. The NLopt nonlinear-optimization package. https://github.com/stevengj/nlopt

Kloek, T., Van Dijk, H.K., 1978. Bayesian estimates of equation system parameters: An application of integration by Monte Carlo. Econometrica 46, 1-19.

Marjou, M., Martinetti, D., Michel, L., Soubeyrand, S., 2019. Risk-based surveillance strategies for early detection of *Xylella fastidiosa* in continental France. In: 2nd European Conference on Xylella Fastidiosa. Ajaccio.

Martinetti, D., Soubeyrand, S., 2019. Identifying lookouts for epidemio-surveillance: application to the emergence of *Xylella fastidiosa* in France. Phytopathology 109 (2), 265-276.

Mastin, A.J., Gottwald, T.R., van Den Bosch, F., Cunniffe, N.J., Parnell, S., 2020. Optimising risk-based surveillance for early detection of invasive plant pathogens. PLoS Biol. 18, e3000863.

Morris, C.E., Géniaux, G., Nédellec, C., Sauvion, N., Soubeyrand, S., 2022. One Health concepts and challenges for surveillance, forecasting, and mitigation of plant disease beyond the traditional scope of crop production. Plant Pathol. 71, 86-97.

Parnell, S., van den Bosch, F., Gottwald, T., Gilligan, C.A., 2017. Surveillance to inform control of emerging plant diseases: An epidemiological perspective. Annu. Rev. Phytopathol. 55, 591-610.

Parnell, S., Gottwald, T., Riley, T., Van Den Bosch, F., 2014. A generic risk-based surveying method for invading plant pathogens. Ecol. Appl. 24, 779-790.

Penttinen, A., Ylitalo, A.-K., 2016. Deducing self-interaction in eye movement data using sequential spatial point processes. Stat. Stat. 17, 1-21.

Perry, J.N., 1994. Sampling and applied statistics for pests and diseases. In: Aspects of Applied Biology, vol. 37.

Picard, C., Soubeyrand, S., Jacquot, E., Thébaud, G., 2019. Analyzing the influence of landscape aggregation on disease spread to improve management strategies. Phytopathology 109, 1198-1207.

Portes, C., Gabriel, E., Ienco, D., 2024a. Environmental and Bio-climatic Data over French Mediterranean areas. Recherche Data Gouv. http://dx.doi.org/10.57745/P7XUII.

Portes, C., Ienco, D., Verdin, E., Gabriel, E., 2024. Environmental and bio-climatic data for epidemiological analysis over french mediterranean areas. Environmental Data Science 3, e31. http://dx.doi.org/10.1017/eds.2024.50.

Powell, M.J.D., 2009. The BOBYQA Algorithm for Bound Constrained Optimization Without Derivatives. Technical Report NA2009/06, Department of Applied Mathematics and Theoretical Physics, Cambridge University, Cambridge, UK.

Pukelsheim, F., 2006. Optimal Design of Experiments. SIAM, Philadelphia.

Sicard, A., Zeilinger, A.R., Vanhove, M., Schartel, T.E., Beal, D.J., Daugherty, M.P., Almeida, R.P., 2018. *Xylella fastidiosa*: insights into an emerging plant pathogen. Annu. Rev. Phytopathol. 56, 181-202.

Soubeyrand, S., de Jerphanion, P., Martin, O., Saussac, M., Manceau, C., Hendrikx, P., Lannou, C., 2018. Inferring pathogen dynamics from temporal count data: The emergence of *Xylella fastidiosa* in France is probably not recent. New Phytol. 219, 824-836.

Soubeyrand, S., Estoup, A., Cruaud, A., Malembic-Maher, S., Meynard, C., Ravigné, V., Barbier, M., Barrès, B., Berthier, K., Boitard, S., et al., 2024. Building integrated plant health surveillance: a proactive research agenda for anticipating and mitigating disease and pest emergence. CABI Agric. Biosci. 5, 72.

Van Lieshout, M.N.M., 2006a. Campbell and moment measures for finite sequential spatial processes. In: Hušková, M., Janžura, M. (Eds.), Proceedings Prague Stochastics 2006. Prague, Czech Republic, August 21-25, 2006, Matfyz Press, pp. 215-224.

Van Lieshout, M.N.M., 2006b. Markovianity in space and time. Lect. Notes- Monogr. Ser. 48, 154-168.

Yazigi, A., Penttinen, A., Ylitalo, A.-K., Maltamo, M., Packalen, P., Mehtätalo, L., 2021. Modeling forest tree data using sequential spatial point processes. J. Agric. Biological Environ. Stat. 27, 88-102.

Zhao, J., Tian, G., Qu, H., 2023. Application of I-Optimal design for modeling and optimizing the operational parameters of ibuprofen granules in continuous twin-screw wet granulation. Biomedicines 11 (7), 2030.