# Pattern-Based Prediction of Population Outbreaks

Gabriel R. Palma1,∗ Wesley A.C. Godoy2 Eduardo Engel2 Douglas Lau3 Edgar
Galvan1,4,6 Oliver Mason5 Charles Markham1,4 Rafael A. Moral1,5

1\. Hamilton Institute, Maynooth University, Maynooth, Ireland;

2\. Department of Entomology and Acarology, University of São Paulo,
Piracicaba, Brazil;

3\. Brazilian Agricultural Research Corporation (Embrapa Trigo), Passo Fundo,
Rio Grande do Sul, Brazil;

4\. Department of Computer Science, Maynooth University, Maynooth, Ireland;

5\. Department of Mathematics and Statistics, Maynooth University, Maynooth,
Ireland;

6\. Naturally Inspired Computation Research Group, Maynooth University,
Ireland;

$\ast$ Corresponding author; e-mail<EMAIL_ADDRESS>
## Abstract

  1. 1.

Insect outbreaks are biotic disturbances in forests and agroecosystems that
cause economic and ecological damage. This phenomenon depends on a variety of
biological and physical factors. The complexity and practical importance of
the issue have made the problem of predicting outbreaks a focus of recent
research.

  2. 2.

Here, we propose the Pattern-Based Prediction (PBP) method for predicting
population outbreaks. It is based on the Alert Zone Procedure, combined with
elements from machine learning. It uses information on previous time series
values that precede an outbreak event as predictors of future outbreaks, which
can be useful when monitoring pest species.

  3. 3.

We illustrate the methodology using simulated datasets and real time series
data obtained by monitoring aphids in wheat crops in Southern Brazil. We
obtained an average test accuracy of $84.6\%$ in the simulation studies
implemented with stochastic models, and $95.0\%$ for predicting outbreaks
using the real dataset. This shows the feasibility of the PBP method in
predicting outbreaks in population dynamics.

  4. 4.

We benchmarked our results against established state-of-the-art machine
learning methods, namely Support Vector Machines, Deep Neural Networks, Long
Short Term Memory and Random Forests. The PBP method yielded a competitive
performance, associated with higher true-positive rates in most comparisons,
while being able to provide interpretability rather than being a black-box
method. This is an improvement over current state-of-the-art machine learning
tools, especially when being used by non-specialists, such as ecologists
aiming to use a quantitative approach for pest monitoring.

  5. 5.

We provide open-source code to implement the PBP method in Python, through the
pypbp package, which may be directly downloaded from the Python Package Index
server or accessed through https://pypbp-documentation.readthedocs.io.

Keywords: Alert zone procedure, deep learning, machine learning, population
dynamics, time series.

## 1 Introduction

Automated systems for syndromic surveillance have been reported in many
studies in different contexts, such as public health aiming to predict disease
outbreaks and also agricultural pests (Madden and Wheelis,, 2003; Buckeridge,,
2007; Büntgen et al.,, 2020; Bright et al.,, 2020; Burkom et al.,, 2021).
These studies have demonstrated many possibilities for predicting outbreaks
based on population dynamics, sampling methods, outbreak frequency, and
threshold analysis. Their results have demonstrated potential to help public
health actions, based on the interpretation of results provided by these
algorithms, using different data sources containing simulated and observed
outbreaks (Buckeridge,, 2007; Chan et al.,, 2021).

Historically, the algorithms used to predict outbreaks involved classical time
series methods, such as ARIMA-type models, seasonal models, and partial
differential equations, among others (Buckeridge,, 2007). These tools
positively impacted public health actions by enhancing the possibility of
predicting disease outbreaks, but the applications of these methods are not
restricted to this area. In quantitative ecology, these applications were
expanded so that many authors started representing ecological phenomena with
mathematical, statistical and machine learning methods (Otto and Day,, 2011;
Odum et al.,, 1971; Ross,, 1998). One recent example is the application of
supervised machine learning methods for predicting infestations of pine trees
by a mountain pine beetle (Ramazi et al.,, 2021).

Examples of these applications are the representation of biological systems
and the interactions between the species, such as predator-prey, host-
parasitoid and competition models (Odum et al.,, 1971; Badkundri et al.,,
2019). Among the taxonomic groups used to study these applications, insects
expand the possibility of developing these methods with various models
inspired by different problems. Examples are the LPA (larva-pupa-adult), host-
parasitoid and other models related to herbivory (Otto and Day,, 2011). These
models are commonly used to study biological phenomena such as outbreaks.
Additionally, researchers can implement further statistical and mathematical
modelling studies with insect time series data.

Insects depend on resource availability, as demonstrated by many studies
involving time series (Nair,, 2001, 2007; Santos et al.,, 2017; Lantschner et
al.,, 2019). The co-evolution between pests and plants reinforces this
influence. It shows the complexity of this system, including biochemical
strategies to avoid herbivory and the genetic plasticity of pests, enhancing
their capability to obtain the necessary resources from plants (Wallner,,
1987; Nair,, 2001, 2007). Pest dynamics are also based on physical and
biological conditions, such as temperature, humidity, precipitation or
irrigation, which strongly influence pest density (Wallner,, 1987; Nair,,
2001; Odum et al.,, 1971). Biological factors such as mating system, life
cycle, number of offspring per generation, mortality, and resource
availability are also capable of influencing insect populations (Wallner,,
1987; Godfray and Godfray,, 1994; Nair,, 2001; Hall et al.,, 2017). In
monoculture scenarios, crop phenology is followed by the presence of pests
reinforcing that resource availability has a strong effect on the population
dynamics in agroecosystems (Nair,, 2001, 2007; Santos et al.,, 2017).

Insect outbreaks have frequently been documented in pest populations (Santos
et al.,, 2017; Lynch,, 2009, 2018). They are important biotic disturbances in
forests and agroecosystems (Wallner,, 1987; Nair,, 2001; Lantschner et al.,,
2019), since they may cause economic and/or ecological damage. Studies have
shown possible explanations for outbreaks, such as abundant resources in
monocultures, absence of predators or parasitoids, genetic factors, and
pheromones produced by pests (Hall et al.,, 2017; Tao et al.,, 2012). Biotic
disturbances can be intensified by climate change and human activities, which
makes it essential to study how they influence pest outbreaks (Volney and
Fleming,, 2000; Sharma and Dhillon,, 2020; Phophi et al.,, 2019). One example
of human activities affecting outbreaks is the occurrence of bark beetles in
temperate forests. These insects have devastated a large area of pine trees in
the continental United States (Negron et al.,, 2008).

In Brazil, there are several examples of pest outbreaks in forests and crops,
as for example Eucalyptus with Thyrinteina arnobia and Stenalcidia sp
(Geometridae) (Zanuncio et al.,, 2006); black wattle with Oncideres impluviata
(Cerambycidae) showing annual outbreaks in the state of Rio Grande do Sul (Ono
et al.,, 2014); soybean with Chrysodeixis includes and Anticarsia gemmatalis
(Lepidoptera: Noctuidae) also exhibiting high frequencies (Bueno et al.,,
2010; Santos et al.,, 2017). Studies focused on the species mentioned above
show that outbreaks occur suddenly because of different natural effects that
can result in increased population densities (Nair,, 2001, 2007; Santos et
al.,, 2017; Lantschner et al.,, 2019). However, the exact reason why an insect
species population suddenly increases in number is still an open question
(Ekholm et al.,, 2019). Outbreak forecasting turns out to be an arduous task
requiring a large amount of man-hours, extensive field work and different
types of specialised equipment. Entomologists traditionally have been using
different interventions to reduce economic damage, mainly in agriculture. The
most common method traditionally used for this task is to define an economic
threshold level (Stern et al.,, 1959; Onstad,, 1987).

The economic threshold is based on the fact that pest density must be
considered in the economic context in which it occurs (Mitchell and Onstad,,
2014). One of the first concepts of economic threshold was published by Stern
et al., (1959) which is defined as “the density at which control measures
should be determined to prevent an increasing pest population from reaching
the economic-injury level”. This concept used the term economic-injury level
which was defined in the same paper as “the lowest population density that
will cause economic damage”. Therefore, both concepts help in decision making
to determine the moment of action, which reduces the pest population density.

Nowadays, the possibilities of actions against insect damage in crops can be
found in the integrated pest management domain, which briefly consists of
employing biological, physical, chemical and genetic approaches to reduce the
population densities of a pest (Stern et al.,, 1959; Goodell,, 2009). The
economic threshold concept has been coupled with more complete analyses
involving control functions given by the crop and pest population information
(Mitchell et al.,, 2004; Dun et al.,, 2009; Tinsley et al.,, 2013). When
outbreaks are frequently recorded in insect populations, the probability that
their population size is bigger than the economic-injury level increases,
causing secondary outbreaks (Goodell,, 2009). Given that this biological
disturbance can occur suddenly and pest monitoring can be delayed, this
density can be bigger than the economic threshold in the subsequent monitoring
sample.

Different approaches have been proposed to address problems of this nature.
One example is the Alert Zone Procedure (AZP) (Hilker and Westerhoff,, 2007),
which consists of scanning observations preceding population outbreak events
to obtain profiles associated with these outbreaks. This method can extract
meaningful information about outbreaks because previous densities before this
event allow for the comprehension of ecological patterns. This method can be
used as a basis to improve pest outbreak forecasting. However, it must be
improved to deal with real-world problems. The effectiveness of these
approaches remains an open question that is the focus of this research.

This paper proposes the Pattern-Based Prediction (PBP) method, which is an
extension of the AZP based on statistical machine learning. We begin by
describing the method and then carry out simulation studies to assess the
performance of our method under different conditions. Finally, we illustrate
our proposed method using a dataset obtained from a pest management system
aimed at monitoring aphids, which are important pests for many different
cultures – such as wheat, barley, and mustard (Kranti et al.,, 2021) – and
discuss the feasibility of applying it to the context of pest management.

## 2 Methods

### 2.1 Generating patterns

Let $x_{t}$ represent the population size of a particular species at time
point $t$, $t=1,\ldots,T$. Initially, we set a population size threshold
$x^{*}$ such that when $x_{t}\geq x^{*}$ we have a population outbreak at time
$t$. We then implement the AZP, as proposed by Hilker and Westerhoff, (2007).
This method consists of scanning observations to identify each outbreak event
$i$, $i=1,\ldots,I$, that occurred at time point $t_{i}$, based on the value
of $x^{*}$, and collecting the $m$ observations that precede them, forming a
vector
$\textbf{p}_{i}^{T}=\left\\{p_{i1},p_{i2},...,p_{im}\right\\}=\left\\{x_{t_{i}-1},x_{t_{i}-2},...,x_{t_{i}-m}\right\\}$
per event. If $t_{i}-m<1$, event $i$ is ignored. After that, we group all
population dynamics patterns that precede these events as the matrix

$\textbf{P}=\begin{bmatrix}x_{t_{1}-1}&x_{t_{1}-2}&\cdots&x_{t_{1}-m}\\\
x_{t_{2}-1}&x_{t_{2}-2}&\cdots&x_{t_{2}-m}\\\ \vdots&\vdots&\ddots&\vdots\\\
x_{t_{I}-1}&x_{t_{I}-2}&\cdots&x_{t_{I}-m}\\\
\end{bmatrix}=\begin{bmatrix}\textbf{p}^{T}_{1}\\\ \textbf{p}^{T}_{2}\\\
\vdots\\\ \textbf{p}^{T}_{I}\par\end{bmatrix},$ (1)

where $I$ represents the total number of identified patterns. See Figure 1(a)
for a plot of all rows of a hypothetical $\mathbf{P}$ matrix. Note that time
series pre-processing may be carried out prior to obtaining the pattern matrix
$P$. For instance, in Section 5 we compare the performance of our method using
the raw time series and a pre-processed series using the Empirical Mode
Decomposition method (Kim et al.,, 2012).

Figure 1: a) The representation of patterns $\mathbf{p}_{i}$ within the matrix
$\mathbf{P}$ that precede an outbreak event, using $m=5$. b) The respective
cluster matrices $\mathbf{P}^{\prime}_{c}$ obtained using
$d^{*}_{\mbox{\scriptsize cluster}}=0.4$. These patterns were obtained from
time series data simulated from a Ricker map, with $r=3$ and $K=1000$,
$x_{1}=200$ and $1000$ observations. The population size threshold for the
outbreak event was set as $x^{*}=2224$ representing the 90% percentile of the
data.

Then, if patterns $i$ and $j$ ($i\neq j$) are sufficiently similar, we group
them in the same cluster. We do this based on the association metric

$d(\textbf{p}_{i},\textbf{p}_{j})=\frac{1}{c(\textbf{p}_{i},\textbf{p}_{j})+1},$
(2)

where
$c(\textbf{p}_{i},\textbf{p}_{j})=\sum_{k=1}^{m}\frac{|p_{ik}-p_{jk}|}{|p_{ik}|+|p_{jk}|}>0$
is the Canberra distance (Androutsos et al.,, 1998; Ehsani and Drabløs,, 2020)
between two vectors, where $|\cdot|$ is the Euclidean norm. This distance is
appropriate for non-negative count data (Androutsos et al.,, 1998). Note that
when $c(\textbf{p}_{i},\textbf{p}_{j})\rightarrow\infty$, then
$d(\textbf{p}_{i},\textbf{p}_{j})\rightarrow 0$, and as
$c(\textbf{p}_{i},\textbf{p}_{j})\rightarrow 0$, then
$d(\textbf{p}_{i},\textbf{p}_{j})\rightarrow 1$.

To define the similarity of patterns we set the value
$d^{*}_{\mbox{\scriptsize cluster}}$, representing the minimum association
metric for considering $\textbf{p}_{i}$ similar to $\textbf{p}_{j}$. This
yields the cluster matrices $\mathbf{P}^{\prime}_{c}$, $c~{}=~{}1,\ldots,C$,
that include patterns which are similar to one another. To obtain these, we
start with pattern $\textbf{p}_{1}$, which represents the first row of the
matrix P. We remove $\textbf{p}_{1}$ from P and add it as the first row of
$\textbf{P}^{\prime}_{1}$. After that we compute the association metric
between $\textbf{p}_{1}$ and all subsequent rows of P. If
$d(\textbf{p}_{1},\textbf{p}_{j})\geq d^{*}_{\mbox{\scriptsize cluster}}$, we
add pattern $\textbf{p}_{j}$ as the last row of the cluster matrix
$\textbf{P}^{\prime}_{1}$ and delete it from P. We repeat this process to
obtain the cluster matrices $\textbf{P}^{\prime}_{c},c=1,\ldots,C\leq I,$
until there are no more rows left in P (Algorithm 1 in the supplementary
material). Note that the order of P is important for the clustering procedure.
See Figure 1(b) for an example where the generated patterns were split into
$C=4$ cluster matrices.

After obtaining the $C$ cluster matrices

$\textbf{P}^{\prime}_{c}=\begin{bmatrix}x_{t_{1}-1}&x_{t_{1}-2}&\cdots&x_{t_{1}-m}\\\
x_{t_{2}-1}&x_{t_{2}-2}&\cdots&x_{t_{2}-m}\\\ \vdots&\vdots&\ddots&\vdots\\\
x_{t_{l^{\prime}_{c}}-1}&x_{t_{l^{\prime}_{c}}-2}&\cdots&x_{t_{l^{\prime}_{c}}-m}\\\
\end{bmatrix},$ (3)

where $l^{\prime}_{c}$ is the number of rows of $\mathbf{P}^{\prime}_{c}$, we
compute the vectors of means $\overline{\mathbf{p}}^{\prime}_{c}$, containing
the mean of each column for cluster matrix $\mathbf{P}^{\prime}_{c}$, to form
the rows of the matrix

$\textbf{P}^{\prime}_{{\mbox{\scriptsize
means}}}=\begin{bmatrix}\displaystyle\frac{1}{l^{\prime}_{1}}\sum\limits_{c=1}^{l^{\prime}_{1}}x_{t_{c}-1}&\displaystyle\frac{1}{l^{\prime}_{1}}\sum\limits_{c=1}^{l^{\prime}_{1}}x_{t_{c}-2}&\cdots&\displaystyle\frac{1}{l^{\prime}_{1}}\sum\limits_{c=1}^{l^{\prime}_{1}}x_{t_{c}-m}\\\
\displaystyle\frac{1}{l^{\prime}_{2}}\sum\limits_{c=1}^{l^{\prime}_{2}}x_{t_{c}-1}&\displaystyle\frac{1}{l^{\prime}_{2}}\sum\limits_{c=1}^{l^{\prime}_{2}}x_{t_{c}-2}&\cdots&\displaystyle\frac{1}{l^{\prime}_{2}}\sum\limits_{c=1}^{l^{\prime}_{2}}x_{t_{c}-m}\\\
\vdots&\vdots&\ddots&\vdots\\\
\displaystyle\frac{1}{l^{\prime}_{C}}\sum\limits_{c=1}^{l^{\prime}_{C}}x_{t_{c}-1}&\displaystyle\frac{1}{l^{\prime}_{C}}\sum\limits_{c=1}^{l^{\prime}_{C}}x_{t_{c}-2}&\cdots&\displaystyle\frac{1}{l^{\prime}_{C}}\sum\limits_{c=1}^{l^{\prime}_{C}}x_{t_{c}-m}\\\
\end{bmatrix}=\begin{bmatrix}\overline{\textbf{p}}^{\prime T}_{1}\\\
\overline{\textbf{p}}^{\prime T}_{2}\\\ \vdots\\\
\overline{\textbf{p}}^{\prime T}_{C}\par\end{bmatrix}.$ (4)

The matrix $\textbf{P}^{\prime}_{{\mbox{\scriptsize means}}}$ contains the
information of all cluster matrices $\textbf{P}^{\prime}_{c}$, and is used for
the prediction of a future event. Given a new collection of observations
$\textbf{x}_{\mbox{\scriptsize new}}$, with length $m$, we compute the
association metric between $\textbf{x}_{\mbox{\scriptsize new}}$ and each row
of $\textbf{P}^{\prime}_{{\mbox{\scriptsize means}}}$. If any computed
association is greater or equal to $d^{*}_{\mbox{\scriptsize pred}}$, shown in
Eq. 5, the threshold for prediction, we predict that an event will occur.
Finally $d^{*}_{\mbox{\scriptsize pred}}$ is defined for each row of
$\textbf{P}^{\prime}_{{\mbox{\scriptsize means}}}$ as a function of
$l^{\prime}_{c}$, the number of patterns that generated each vector of means:

$d^{*}_{\mbox{\scriptsize pred}}=f(l^{\prime}_{c})=d^{*}_{\mbox{\scriptsize
base}}+\frac{(1-d^{*}_{\mbox{\scriptsize base}})}{(l^{\prime}_{c})^{\alpha}},$
(5)

where $d^{*}_{\mbox{\scriptsize base}}$ is the baseline value of the
association metric (the smallest it is allowed to be) and $\alpha$ is a
constant that changes the shape of the function $f$ (see Figure 2). When
$l^{\prime}_{c}\rightarrow\infty$, we have that $d^{*}_{\mbox{\scriptsize
pred}}\rightarrow d^{*}_{\mbox{\scriptsize base}}$, and as
$l^{\prime}_{c}\rightarrow 1$, also $d^{*}_{\mbox{\scriptsize
pred}}\rightarrow 1$. This means that for predicting that a new event will
occur, we would need a larger association between
$\mathbf{x}_{\mbox{\scriptsize new}}$ and a particular
$\overline{\mathbf{p}}^{\prime}_{c}$ that was obtained from a small number of
patterns.

Figure 2: The threshold for prediction $d^{*}_{\mbox{\scriptsize pred}}$,
calculated as a function of $l^{\prime}_{c}$ for $\alpha=1$ (red curve),
$\alpha=3$ (green curve) and $\alpha=0.25$ (blue curve), whilst fixing
$d^{*}_{\mbox{\scriptsize base}}=0.6$. The $x$-axis is represented as
$1/l^{\prime}_{c}$ to ease visualisation.

To summarise, the PBP method consists of the following steps:

  1. 1.

Choose the value of the population size threshold $x^{*}$;

  2. 2.

Set the values of $m$ and $d^{*}_{\mbox{\scriptsize cluster}}$;

  3. 3.

Generate the pattern matrix $\mathbf{P}$;

  4. 4.

Obtain the cluster matrices $\mathbf{P}^{\prime}_{c}$ (algorithm presented in
the supplementary material);

  5. 5.

Compute the matrix $\textbf{P}^{\prime}_{{\mbox{\scriptsize means}}}$ from the
column means of each cluster matrix $\mathbf{P}^{\prime}_{c}$;

  6. 6.

Set the values of $d^{*}_{\mbox{\scriptsize base}}$ and $\alpha$ and obtain
$d^{*}_{\mbox{\scriptsize pred}}$ for each row of
$\textbf{P}^{\prime}_{{\mbox{\scriptsize means}}}$;

  7. 7.

Given a new collection of observations $\textbf{x}_{\mbox{\scriptsize new}}$
compute the proposed association metric between $\textbf{x}_{\mbox{\scriptsize
new}}$ and each row of $\textbf{P}^{\prime}_{{\mbox{\scriptsize means}}}$;

  8. 8.

If the computed association coefficient is greater than or equal to the value
of $d^{*}_{\mbox{\scriptsize pred}}$ associated with that row of
$\textbf{P}^{\prime}_{{\mbox{\scriptsize means}}}$, predict that a new event
will occur at the next time step; predict that it will not occur, otherwise.

A schematic diagram of this process is represented in Figure 3.

Figure 3: A schematic representation of the pattern-based method used to
predict an outbreak based on time series data.

### 2.2 Choosing $m$, $d^{*}_{\mbox{\scriptsize cluster}}$,
$d^{*}_{\mbox{\scriptsize base}}$ and $\alpha$ via cross-validation

We propose the use of k-fold cross-validation to choose the values of $m$,
$d^{*}_{\mbox{\scriptsize cluster}}$ and $\alpha$, such that the accuracy of
the method is optimized. Here, the k-fold cross validation briefly consists of
creating $k$ groups of patterns $\mathbf{p}_{i}$ of the pattern matrix
$\mathbf{P}$ and by removing the first group of patterns from $\mathbf{P}$,
obtaining $\textbf{P}^{\prime}_{{\mbox{\scriptsize means}}}$ without using the
information of this group, and carrying out the method to predict the
occurrence of events based on the first group. After obtaining the method
predictions based on this group, we compute

  * •

true positives (TP): the number of times the method accurately predicted an
event;

  * •

true negatives (TN): the number of times the method accurately predicted there
was no event;

  * •

false positives (FP): the number of times the method incorrectly predicted an
event;

  * •

and false negatives (FN): the number of times the method incorrectly predicted
there was no event.

We use these values to obtain the accuracy,
$\mbox{ACC}=\frac{\mbox{TP}+\mbox{TN}}{\mbox{TP}+\mbox{TN}+\mbox{FP}+\mbox{FN}}$,
the true positive rate, $\mbox{TPR}=\frac{\mbox{TP}}{\mbox{TP}+\mbox{FN}}$,
and the false positive rate,
$\mbox{FPR}=\frac{\mbox{FP}}{\mbox{TN}+\mbox{FP}}$. We repeat this process for
each group ending up with $k$ values of these metrics. To measure the overall
performance, we obtain the average of these metrics. Here, we carry out the
analysis using $k=5$.

To optimise the predictive power of the method, firstly we fix the values of
$m$, $d^{*}_{\mbox{\scriptsize cluster}}$ and $\alpha$, and obtain different
TPR and FPR values by varying $d^{*}_{\mbox{\scriptsize base}}$. The TPR and
FPR can be plotted against each other to form a ROC curve (Hastie et al.,,
2004). This curve is bounded between 0 and 1. For a method with good
predictive power, we expect the area under the curve (AUROC) to be close to 1.
Let $g(m,d^{*}_{\mbox{\scriptsize cluster}},\alpha)$ be an objective function
that returns $-\mbox{AUROC}$ based on the described method. Then, we use the
Generalized Simulated Annealing method (Tsallis,, 1988; Tsallis and Stariolo,,
1996; Xiang et al.,, 1997; Xiang and Gong,, 2000; Xiang et al.,, 2013;
Mullen,, 2014) to obtain the values of $m$, $d^{*}_{\mbox{\scriptsize
cluster}}$ and $\alpha$ that minimise $g$. This method speeds up the
computation process compared to a grid search of the variables. Other methods
may be used, such as Differential Evolution (see the supplementary materials).

To construct the ROC curve, we vary $d^{*}_{\mbox{\scriptsize base}}$ from $0$
to $1$ using increments of $0.1$ and calculate the AUROC using the trapezoid
method (Liu and Pierce,, 1994). Finally, we apply one of two forms to choose
$d^{*}_{\mbox{\scriptsize base}}$: the first is based on selecting a minimum
threshold for the true positive rate (e.g. 0.8 or 0.9) and the second on
selecting a maximum threshold for the false positive rate (e.g. 0.1 or 0.2).
In summary, the method consists of the following steps:

  1. 1.

Fix the values of $m$, $d^{*}_{\mbox{\scriptsize cluster}}$ and $\alpha$;

  2. 2.

For different values of $d^{*}_{\mbox{\scriptsize base}}$, carry out $k$-fold
cross validation and obtain the TPR and FPR for each fold, and compute the
AUROC;

  3. 3.

Choose $m$, $d^{*}_{\mbox{\scriptsize cluster}}$ and $\alpha$ such that the
AUROC is the largest;

  4. 4.

Choose $d^{*}_{\mbox{\scriptsize base}}$ based on the minimum TPR or maximum
FPR that would be allowed in the study, based on the ROC curve with the
largest area. (Note that the minimum TPR or maximum FPR allowed depends highly
on the ecological system and objectives of the monitoring programme.)

### 2.3 Sensitivity analysis

We carried out a sensitivity analysis using simulated deterministic population
dynamics, obtained from the Ricker map (Otto and Day,, 2011):

$x_{t+1}=x_{t}\exp{\left[r\left(1-\frac{x_{t}}{K}\right)\right]},$ (6)

where $x_{t}$ denotes the population size of an organism at time $t$ and the
parameters $r>0$ and $K>0$ describe the intrinsic growth rate and carrying
capacity of the environment, respectively. We simulated $100$ generations
using $r=3$, $K=1000$ and an initial value of $x_{1}=200$.

In order to study the influence of different values of $m$ and
$d^{*}_{\mbox{\scriptsize cluster}}$ on $C$ (the total number of cluster
matrices $\mathbf{P}^{\prime}$) and the overall accuracy of the method, we
used the simulated observations from the Ricker map, setting $x^{*}=2224$ as
the threshold for an outbreak event. This value corresponds to the 90%
percentile of the simulated values from the deterministic Ricker map setting
the parameter values as described above. We then employed the methodology
described in the previous sections to obtain the matrices
$\textbf{P}^{\prime}_{{\mbox{\scriptsize means}}}$ for $m$ varying from 2 to
15 in increments of 1, and $d^{*}_{cluster}$ varying from 0 to 1 with
increments of $0.1$.

### 2.4 Method validation under stochastic conditions

To study the accuracy of the proposed method in predicting outbreak and
extinction risk events under stochastic conditions, we simulated from three
different approaches. The first included an additive Gaussian error
$\varepsilon_{t}\sim\mbox{Normal}(0,\sigma^{2})$ in the Ricker map, yielding
the recurrence equation

$x_{t+1}=x_{t}\exp{\left[r\left(1-\frac{x_{t}}{K}\right)\right]}+\varepsilon_{t+1}.$
(7)

Whenever the addition of the random noise term yielded $x_{t+1}<0$, a new
random noise value would be drawn from the normal distribution until
$x_{t+1}>0$, to ensure positive population sizes.

The second approach utilized a state-space formulation using a Poisson
distribution with the mean term $\mu$ given by the Ricker recurrence equation,
i.e.

$\displaystyle X_{1}$ $\displaystyle\sim$
$\displaystyle\mbox{Poisson}(\mu_{1}=x_{1})$ (8) $\displaystyle X_{t+1}|X_{t}$
$\displaystyle\sim$
$\displaystyle\mbox{Poisson}\left(\mu_{t+1}=X_{t}\exp\left[r\left(1-\frac{X_{t}}{K}\right)\right]\right),$
(9)

from which all $x_{t}$ values were drawn recursively. Finally, the third
approach utilizes a state-space formulation based on a negative binomial
distribution to accommodate overdispersion in the simulation study, i.e.
$X_{t+1}|X_{t}\sim\mbox{Negative Binomial}(\mu_{t+1},\phi)$. We estimated $r$
and $K$ based on real time series data of aphid counts in Southern Brazil for
each model formulation, as well as the dispersion parameters $\sigma^{2}$ for
the Gaussian model and $\phi$ for the negative binomial model.

Using the parameter estimates in Table 1, we simulated $20$ samples of size
$400$ for each model. We also simulated $20$ samples of size $400$ using the
negative binomial model with $\phi=3$, to introduce a scenario with stronger
overdispersion. We computed the accuracy, TPR and FPR by training the methods
with the initial $80\%$ observations and testing with $20\%$ of the time
series. Moreover, based on the ROC curve with the largest AUROC, we chose
$d^{*}_{\mbox{\scriptsize base}}$ using four methods:

  1. 1.

‘TPR_08’: choose the $d^{*}_{\mbox{\scriptsize base}}$ value associated with
the smallest TPR value that is equal to or greater than 0.8;

  2. 2.

‘TPR_09’: choose the $d^{*}_{\mbox{\scriptsize base}}$ value associated with
the smallest TPR value that is equal to or greater than 0.9;

  3. 3.

‘FPR_01’: choose the $d^{*}_{\mbox{\scriptsize base}}$ value associated with
the largest FPR value that is equal to or less than 0.1;

  4. 4.

‘FPR_02’: choose the $d^{*}_{\mbox{\scriptsize base}}$ value associated with
the largest FPR value that is equal to or less than 0.2.

Parameter | Model
---|---
Gaussian | Poisson | Negbin
$r$ | 0.15 | 0.28 | 0.57
$K$ | 224 | 310 | 370
$\sigma^{2}$ | $21,866$ | $-$ | $-$
$\phi$ | $-$ | $-$ | 1.2
AIC | $5,226$ | $34,913$ | $4,296$
Table 1: Parameter estimates obtained when fitting the Ricker state-space
model to the aphid data assuming different distributions for the observation
process, namely Gaussian, Poisson and negative binomial, as well as the Akaike
Information Criterion (AIC) for each model fit. Negbin = negative binomial.

### 2.5 Analysis of case-study: Aphid data

To illustrate the predictive performance of our method, we use data obtained
from an aphid monitoring programme implemented in Southern Brazil (State of
Rio Grande do Sul, RS). These insects are considered as important pest species
of many crops. For instance, among the aphid species monitored by this
programme, the species Rhopalosiphum padi and Rhopalosiphum rufiabdominalis
are widely considered important pest species associated with winter cereals,
and are found in the Eurasian region with a cosmopolitan distribution
(Macfadyen and Kriticos,, 2012). Sampling was carried out weekly in an area of
$5500m^{2}$ in a wheat culture region (Coxilha, RS, 710 m altitude,
$28^{\circ}11^{\prime}42.8^{\prime\prime}$ S and
$52^{\circ}19^{\prime}30.6^{\prime\prime}$ W), from 2011 to 2019, totalling
424 observations. The temperature and relative humidity data were monitored at
the Passo Fundo weather station ($28^{\circ}15^{\prime}$ S,
$52^{\circ}24^{\prime}$ W, 684 m), located 10 km from the experimental area.
The field was cultivated under a no-till system.

The species of aphids were monitored using Moericke traps (yellow tray, 45 cm
long x 30 cm wide x 4.5 cm high), filled with a solution (2 L) consisting of
water, 40% formalin (0, 3%) and detergent (0.2%). Each tray had three lateral
holes (5 mm in diameter) close to the edge, protected by a thin screen to
prevent leaks and loss of solid content during rain. Four traps were
distributed at the borders of the crop rotation tests. The traps were levelled
at approximately 20 cm from the floor with bricks. The crop rotation area was
cultivated with cereals (oat, wheat, and triticale), radish and fallow during
the Winter and in the Summer with soybeans, corn, and Brachiaria sp. Every
seven days, the solid content of the trays was separated from the solution
through the sieve and collected. The biological material was preserved in a
glass bottle with 70% alcohol. Aphids and parasitoids were separated,
identified, and counted under a stereomicroscope in the laboratory.

Monitoring is one of the bases for Integrated Pest Management (IPM). For
aphids, their importance stands out mainly due to the ability of these insects
to transmit viruses (Barley Yellow Dwarf Virus) to economically important
crops, such as wheat, oats, barley and rye. The criterion to find a threshold
($x^{*}$) was the total number of aphids observed in the four traps. Usually,
$10\%$ of plants infested by aphids results in an economic threshold (Bell et
al.,, 2015). This percentage corresponds to 50 insects per trap, totalling 200
aphids, which is the threshold value used here to define an outbreak.

We selected 40%, 50%, 60% and 70% of the initial observations of the time
series for training and the complement as test sets to obtain the accuracy,
true-positive rates and false-positive rates. We compared the performance of
the PBP method with Random Forests (RF), Support Vector Machines (SVM), Deep
Neural Networks (DNNs) and Long Short Term Memory (LSTM) algorithms. To obtain
the algorithm performance using these competing methodologies, we created the
matrix P with $m=1$, 4 and 7. Another matrix that did not contain outbreaks
was generated with the same $m$ observations before a threshold lower than the
population size $x^{*}=200$ that defines an outbreak of aphids in the study
area.

For RF, we used 2 splitting predictors per tree (for $m=4$ and 7), and a total
of 1,000 trees. For SVM, we used the linear kernel. For DNN, after
experimenting with different architectures, we used 17 hidden layers with 11,
12, 13, 14, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, and 3 neurons each, all
using the rectified linear unit (ReLU) activation function; for the output
layer we used the sigmoid activation function. We used the LSTM method by
adding a long short term memory block on the DNN’s architecture previously
described. The implementation of the classification methods was carried out in
Python using the libraries TensorFlow and Scikit-learn (Van Rossum and Drake
Jr,, 1995). Finally, we obtained the accuracy, true positive rate and false
positive rate using the aforementioned training-test splits of the time
series.

## 3 Results and discussion

### 3.1 Sensitivity analysis

We found that $d^{*}_{cluster}$ is proportional to the number of cluster
matrices ($C$) created by the proposed method (see Figure 4(a)). However, as
$d^{*}_{cluster}$ reaches values higher than $0.45$, the parameter $m$ did not
influence the accuracy of our methods. It indicates that fixed values of $m$
could be used when we use such values of $d^{*}_{cluster}$ (Figure 4(b)).
These findings highlight the importance of using optimisation procedures to
choose the appropriate value of $d^{*}_{cluster}$ for each study.

Figure 4: The effect of $m$ and $d^{*}_{cluster}$ on (a) $C$ (i.e. number of
cluster-matrices $\textbf{P}^{\prime}_{c}$), and on (b) the accuracy of the
proposed method. These results were obtained from time series data simulated
from a Ricker map, with $r=3$ and $K=1000$, with $x_{1}=200$. The population
size threshold for the outbreak event was set as $x^{*}=2224$ representing the
90% percentile of the data.

The creation of cluster matrices and subsequent sensitivity analysis carried
out here is an increment to the studies conducted by Hilker and Westerhoff,
(2007). By using the AZP as a basis to create $\textbf{P}^{\prime}_{c}$ we
enhance the information which can be extracted from the population dynamics of
any species of interest. This process can extract the different population
states using the dynamics obtained from monitoring programmes. Also, by
grouping these states into different cluster matrices, we can observe the
frequency of each patterns group occurring before the outbreak. The process of
clustering will help us to perform the outbreak classification based on the
different pattern numbers contained in $\textbf{P}^{\prime}_{c}$. This result
reflects the accuracy of $100\%$ obtained in some parameter regions of the
sensitivity study. As shown in the following sections, the results are
dramatically improved when we fully optimise our choice for the parameter
values, even when subject to stochastic effects.

### 3.2 Method validation under stochastic conditions

The accuracy of our method for predicting population outbreaks obtained from
the stochastic simulation scenarios using the raw time series data and pre-
processing the data using Empirical Mode Decomposition (EMD) considering all
models were, respectively, on average $73.8\%$ with a standard deviation of
$23.4\%$ and $73.2\%$ with a standard deviation of $24.1\%$ (see Figure 5).
The average FPR obtained was $25.3\%$ with a standard deviation of $29.0\%$
and $25.1\%$ with a standard deviation of $28.7\%$. Finally, the average TPR
were $53.2\%$ with a standard deviation of $40.9\%$ and $55.6\%$ with a
standard deviation of $41.9\%$. Therefore, we found that there are no
differences in performance when pre-processing the data using EMD.

Figure 5: Accuracy simulation results, TPR (True Positive Rate) and FPR (False
Positive Rate) using the raw simulated time series and pre-processed series
using Empirical Mode Decomposition (EMD). In both scenarios four methods were
used to choose $d^{*}_{\mbox{\scriptsize base}}$: based on a maximum FPR (0.1
and 0.2) or a minimum TPR (0.8 and 0.9).

Also, considering the influence of the model, we found that the rank of models
in which our method produced higher performance is, respectively, starting
with the best one, the negative binomial ($\phi=1.2$), Gaussian, negative
binomial ($\phi=3$) and Poisson stochastic models. On average, we obtain an
accuracy of $84.6\%$ with a standard deviation of $20.5\%$, a false positive
rate of $14.9\%$ with a standard deviation of $24.2\%$ and a true positive
rate of $59.6\%$ with a standard deviation of $42.0\%$ for the negative
binomial ($\phi=1.2$) model. Considering the Gaussian model, we obtained an
accuracy of $75.2\%$ with a standard deviation of $18.0\%$, a false positive
rate of $22.9\%$ with a standard deviation of $22.2\%$ and a true positive
rate of $55.5\%$ with a standard deviation of $37.3\%$. This finding makes our
method a promising prediction tool since we got good results even when using
stochastic approaches to simulate the data.

### 3.3 Analysis of case-study: Aphid data

Figure 6: The time series represents the total aphids collected within the
four traps on the monitoring system on time. The red line represents the
threshold $x^{*}=200$, the green line is the original time series, and the
black line is the result of the empirical mode decomposition method.

To predict the threshold representing an outbreak for the aphid population
dynamics (Figure 6), we select $x^{*}=200$ considering the number of species
collected in the four traps of the monitoring system, which was related to the
$10\%$ of plants infected by aphids resulting in the economic threshold.
Applying PBP using different training data obtained from percentages of the
initial observations of the time series of aphids, the accuracy values were
higher than $70\%$ regardless of the percentage of training using the original
aphid time series.

Table 2:  Prediction accuracy, true-positive rate (TPR) and false positive
rate (FPR) obtained from the Pattern-Based Prediction (PBP) and competing
methods Random Forests (RF), Support Vector Machines (SVM), Deep Neural
Networks (DNNs) and Long Short Term Memory (LSTM) methods with different
values of $m$ (i.e. observations before the event). All methods were carried
out using training sets with $40\%,50\%,60\%,70\%$ and $80\%$ of the initial
observation of the aphid time series.

Metrics Train PBP PBP PBP PBP RF RF RF SVM SVM SVM DNN DNN DNN LSTM LSTM LSTM
percentage ($FPR=0.1$) ($FPR=0.2$) ($TPR=0.8$) ($TPR=0.9$) $(m=1)$ $(m=4)$
$(m=7)$ $(m=1)$ $(m=4)$ $(m=7)$ $(m=1)$ $(m=4)$ $(m=7)$ $(m=1)$ $(m=4)$
$(m=7)$ 0.4 0.81 0.81 0.71 0.49 0.86 0.91 0.90 0.89 0.89 0.89 0.90 0.88 0.93
0.89 0.89 0.94 0.5 0.87 0.87 0.87 0.87 0.87 0.91 0.90 0.87 0.87 0.88 0.89 0.87
0.93 0.87 0.88 0.94 Accuracy 0.6 0.90 0.88 0.83 0.83 0.89 0.92 0.90 0.91 0.89
0.89 0.90 0.89 0.93 0.90 0.89 0.94 0.7 0.94 0.94 0.92 0.87 0.92 0.95 0.95 0.95
0.94 0.94 0.92 0.94 0.98 0.94 0.94 0.96 0.8 0.95 0.90 0.90 0.90 0.95 0.95 0.96
0.96 0.95 0.94 0.94 0.95 0.98 0.95 0.94 0.99 0.4 0.52 0.52 0.78 0.96 0.42 0.60
0.55 0.07 0.04 0.03 0.41 0.00 0.34 0.00 0.11 0.53 0.5 0.18 0.18 0.18 0.18 0.50
0.68 0.60 0.04 0.04 0.00 0.42 0.00 0.36 0.00 0.12 0.64 TPR 0.6 0.12 0.69 0.87
0.87 0.44 0.60 0.55 0.12 0.00 0.00 0.31 0.00 0.31 0.00 0.06 0.56 0.7 0.28 0.28
0.43 0.71 0.29 0.60 0.60 0.14 0.00 0.00 0.14 0.00 0.71 0.00 0.00 0.43 0.8 0.50
0.75 0.75 0.75 0.50 0.50 0.67 0.25 0.00 0.00 0.25 0.00 0.75 0.00 0.00 0.75 0.4
0.15 0.15 0.30 0.57 0.10 0.05 0.04 0.00 0.00 0.00 0.04 0.00 0.00 0.00 0.00
0.01 0.5 0.05 0.05 0.05 0.05 0.07 0.04 0.03 0.00 0.00 0.00 0.04 0.00 0.00 0.00
0.00 0.02 FPR 0.6 0.01 0.09 0.18 0.18 0.07 0.04 0.04 0.00 0.01 0.01 0.03 0.00
0.00 0.00 0.01 0.01 0.7 0.02 0.02 0.05 0.11 0.04 0.02 0.02 0.00 0.00 0.00 0.03
0.00 0.00 0.00 0.00 0.00 0.8 0.03 0.09 0.09 0.09 0.01 0.03 0.01 0.00 0.00 0.00
0.03 0.00 0.00 0.00 0.00 0.00

Table 2 shows that the PBP method is competitive with state-of-the-art machine
learning methods, such as the commonly used Random Forest (RF) algorithm. The
criterion of a false positive rate of at most $0.2$ provides an accuracy of
$91.0\%$, a false positive rate of $8\%$ and a true positive rate of $75\%$.
Only DNN could obtain similar values of true positive rate for this case
study. Moreover, our method using the criteria based on the true positive rate
values of a minimum of $0.8$ and $0.9$ could obtain higher values of $96\%$.
On the other hand, the method got an accuracy of $71\%$ and a false positive
rate of $32\%$ in both cases. It indicates the flexibility of our approach
because users can allow for a trade-off between false positives and false
negatives. In practice, a false positive would result in using pest control
techniques unnecessarily, while a false negative could result in failing to
control the pest, which may cause economic damage due to an outbreak
occurring.

Figure 7: Outbreak patterns obtained from the proposed method using the aphid
time series. Each of the red lines represents a row of the
$\textbf{P}^{\prime}_{{\mbox{\scriptsize means}}}$ matrix. The intervals are
the $25\%$ and $75\%$ percentiles of the patterns that generated each vector
of means. The black line represents an observed series
($\textbf{x}_{\mbox{\scriptsize new}}$), for which the association metric $d$
is calculated between each of the three identified patterns. If
$d>d^{*}_{\mbox{\scriptsize pred}}$, then the PBP method would classify
$\textbf{x}_{\mbox{\scriptsize new}}$ as preceding an outbreak event. The
calculated $d^{*}_{\mbox{\scriptsize pred}}$ values for the three patterns
were $0.37$, $0.61$ and $0.42$, whereas the association metrics between
$\textbf{x}_{\mbox{\scriptsize new}}$ and each pattern were $0.20$, $0.26$ and
$0.15$, respectively. Therefore, $\textbf{x}_{\mbox{\scriptsize new}}$ would
be classified as not preceding an outbreak.

With respect to interpretability, when using RF, it is possible to obtain
variable importance. However in this case, this will tell us which previous
steps were most important when predicting outbreaks, not necessarily how they
relate to its occurrence. When using LSTM and DNN, which are commonly referred
to as ‘black-box’ methods (Liang et al.,, 2021), it is even more challenging
to find explainable frameworks that allow us to study the relationship between
the predictors and the outbreaks. However, each hyperparameter in the PBP
framework provides a clear interpretation, and we are able to create visual
representations of the patterns that occurred before the outbreak (the
$\textbf{P}^{\prime}_{{\mbox{\scriptsize means}}}$ matrix). For instance,
Figure 7 displays the three patterns in
$\textbf{P}^{\prime}_{{\mbox{\scriptsize means}}}$ obtained from employing PBP
using the optimised hyperparameter values for the aphid data using 50% of the
time series for training.

Additionally, the number of patterns encountered in each cluster matrix
presents the importance of each clustered pattern for predicting animal
outbreaks. The parameter $\alpha$ shows how relevant each group of patterns is
for predicting an outbreak. Also, we can assess the importance of clustering
the pattern matrix by looking at the estimate of $d^{*}_{\mbox{\scriptsize
cluster}}$. Larger values typically indicate fewer recognised patterns in
$\textbf{P}^{\prime}_{{\mbox{\scriptsize means}}}$. The $m$ hyperparameter
shows the number of previous observations required to provide a classification
based on our method, so it provides a clear interpretation for ecologists and
farmers in terms of how far in the past to watch for when identifying
patterns. The $d^{*}_{\mbox{\scriptsize base}}$ hyperparameter informs the
minimum degree of similarity that is required to classify an outbreak, based
on previously observed patterns. Therefore, not only is the PBP method
competitive when compared to state-of-the-art machine learning methods, it is
also interpretable, and brings descriptive advantages combined with its
predictive power.

## Acknowledgments

This publication has emanated from research conducted with the financial
support of Fundação de Amparo à Pesquisa do Estado de São Paulo (proc. no.
19/14805-7 and no. 20/06147-7), Agência USP de Inovação and Science Foundation
Ireland under Grant 18/CRT/6049. The opinions, findings and conclusions or
recommendations expressed in this material are those of the authors and do not
necessarily reflect the views of the funding agencies.

## Literature Cited

  * Androutsos et al., (1998) Androutsos, D., Plataniotiss, K., and Venetsanopoulos, A. N. (1998).  Distance measures for color image retrieval.  In Proceedings 1998 International Conference on Image Processing. ICIP98 (Cat. No. 98CB36269), volume 2, pages 770–774. IEEE.
  * Badkundri et al., (2019) Badkundri, R., Valbuena, V., Pinnamareddy, S., Cantrell, B., and Standeven, J. (2019).  Forecasting the 2017-2018 yemen cholera outbreak with machine learning.  arXiv preprint arXiv:1902.06739.
  * Bell et al., (2015) Bell, J. R., Alderson, L., Izera, D., Kruger, T., Parker, S., Pickup, J., Shortall, C. R., Taylor, M. S., Verrier, P., and Harrington, R. (2015).  Long-term phenological trends, species accumulation rates, aphid traits and climate: Five decades of change in migrating aphids.  Journal of Animal Ecology, 84(1):21–34.
  * Bright et al., (2020) Bright, B. C., Hudak, A. T., Meddens, A. J., Egan, J. M., and Jorgensen, C. L. (2020).  Mapping multiple insect outbreaks across large regions annually using landsat time series data.  Remote Sensing, 12(10):1655.
  * Buckeridge, (2007) Buckeridge, D. (2007).  Outbreak detection through automated surveillance: A review of the determinants of detection.  Journal of biomedical informatics, 40:370–9.
  * Bueno et al., (2010) Bueno, R. C., de Freitas Bueno, A., Moscardi, F., Parra, J. R., and Hoffmann-Campo, C. (2010).  Lepidopteran larva consumption of soybean foliage: Basis for developing multiple-species economic thresholds for pest management decisions.  Pest management science, 67:170–4.
  * Büntgen et al., (2020) Büntgen, U., Liebhold, A., Nievergelt, D., Wermelinger, B., Roques, A., Reinig, F., Krusic, P. J., Piermattei, A., Egli, S., Cherubini, P., et al. (2020).  Return of the moth: rethinking the effect of climate on insect outbreaks.  Oecologia, 192(2):543–552.
  * Burkom et al., (2021) Burkom, H., Loschen, W., Wojcik, R., Holtry, R., Punjabi, M., Siwek, M., Lewis, S., et al. (2021).  Electronic surveillance system for the early notification of community-based epidemics (essence): Overview, components, and public health applications.  JMIR public health and surveillance, 7(6):e26303.
  * Chan et al., (2021) Chan, T.-C., Tang, J.-H., Hsieh, C.-Y., Chen, K. J., Yu, T.-H., and Tsai, Y.-T. (2021).  Approaching precision public health by automated syndromic surveillance in communities.  Plos one, 16(8):e0254479.
  * Dun et al., (2009) Dun, Z., Mitchell, P., and Agosti, M. (2009).  Estimating diabrotica virgifera virgifera damage functions with field trial data: applying an unbalanced nested error component model.  Journal of Applied Entomology, 134:409 – 419.
  * Ehsani and Drabløs, (2020) Ehsani, R. and Drabløs, F. (2020).  Robust distance measures for k nn classification of cancer data.  Cancer informatics, 19:1176935120965542.
  * Ekholm et al., (2019) Ekholm, A., Tack, A., Pulkkinen, P., and Roslin, T. (2019).  Host plant phenology, insect outbreaks and herbivore communities – the importance of timing.  Journal of Animal Ecology.
  * Godfray and Godfray, (1994) Godfray, H. C. J. and Godfray, H. (1994).  Parasitoids: behavioral and evolutionary ecology, volume 67.  Princeton University Press.
  * Goodell, (2009) Goodell, P. (2009).  Fifty years of the integrated control concept: The role of landscape ecology in ipm in san joaquin valley cotton.  Pest management science, 65:1293–7.
  * Hall et al., (2017) Hall, A., Johnson, S., Cook, J., and Riegler, M. (2017).  High nymphal host density and mortality negatively impact parasitoid complex during an insect herbivore outbreak.  Insect Science, 26.
  * Hastie et al., (2004) Hastie, T., Tibshirani, R., Friedman, J., and Franklin, J. (2004).  The elements of statistical learning: Data mining, inference, and prediction.  Math. Intell., 27:83–85.
  * Hilker and Westerhoff, (2007) Hilker, F. and Westerhoff, F. (2007).  Preventing extinction and outbreaks in chaotic populations.  The American naturalist, 170:232–41.
  * Kim et al., (2012) Kim, D., Kim, K. O., and Oh, H.-S. (2012).  Extending the scope of empirical mode decomposition by smoothing.  EURASIP Journal on Advances in Signal Processing, 2012(1):1–17.
  * Kranti et al., (2021) Kranti, W., Nivedita, G., and Shindikar, M. (2021).  Understanding the plant aphid interaction: A review.  European Journal of Biology and Biotechnology, 2(6):1–6.
  * Lantschner et al., (2019) Lantschner, V., Aukema, B., and Corley, J. (2019).  Droughts drive outbreak dynamics of an invasive forest insect on an exotic host.  Forest Ecology and Management, 433:762–770.
  * Liang et al., (2021) Liang, Y., Li, S., Yan, C., Li, M., and Jiang, C. (2021).  Explaining the black-box model: A survey of local interpretation methods for deep neural networks.  Neurocomputing, 419:168–182.
  * Liu and Pierce, (1994) Liu, Q. and Pierce, D. A. (1994).  A note on Gauss—Hermite quadrature.  Biometrika, 81(3):624–629.
  * Lynch, (2009) Lynch (2009).  Spruce aphid, elatobium abietinum (walker): Life history and damage to engelmann spruce in the pinaleno mountains, arizona.  The Last Refuge of the Mt. Graham Red Squirrel: Ecology of Endangerment.
  * Lynch, (2018) Lynch, A. (2018).  Socioecological impacts of multiple forest insect outbreaks in the pinaleño spruce–fir forest, arizona.  Journal of Forestry, 117.
  * Macfadyen and Kriticos, (2012) Macfadyen, S. and Kriticos, D. (2012).  Modelling the geographical range of a species with variable life-history.  Public Library of Science One, 7.
  * Madden and Wheelis, (2003) Madden, L. and Wheelis, M. (2003).  The threat of plant pathogens as weapons against u.s. crops.  Annual review of phytopathology, 41:155–76.
  * Mitchell et al., (2004) Mitchell, P., Gray, M., and Steffey, K. (2004).  A composed-error model for estimating pest-damage functions and the impact of the western corn rootworm soybean variant in illinois.  American Journal of Agricultural Economics, 86:332–344.
  * Mitchell and Onstad, (2014) Mitchell, P. D. and Onstad, D. W. (2014).  Valuing pest susceptibility to control.  In Insect Resistance Management, pages 25–53. Elsevier.
  * Mullen, (2014) Mullen, K. M. (2014).  Continuous global optimization in r.  Journal of Statistical Software, 60(1):1–45.
  * Nair, (2001) Nair, K. (2001).  Pest outbreaks in tropical forest plantations: is there a greater risk for exotic tree species?  Center for International Forestry Research, pages 1–82.
  * Nair, (2007) Nair, K. (2007).  Tropical forest insect pests. ecology, impact, and management.  Tropical Forest Insect Pests: Ecology, Impact, and Management, pages 1–404.
  * Negron et al., (2008) Negron, J., Bentz, B., Fettig, C., Gillette, N., Hansen, E., Hayes, J., Kelsey, R., Lundquist, J., Lynch, A., Progar, R., and Seybold, S. (2008).  Us forest service bark beetle research in the western united states: Looking toward the future.  Journal of Forestry, 106:325–331.
  * Odum et al., (1971) Odum, E. P., Barrett, G. W., et al. (1971).  Fundamentals of ecology, volume 3.  Saunders Philadelphia.
  * Ono et al., (2014) Ono, M. A., Ferreira, E. N. L., and Godoy, W. A. C. (2014).  Black wattle insect pest currently in brazil.  Glo Adv Res J Agric Sci, 3(12):409–414.
  * Onstad, (1987) Onstad, D. W. (1987).  Calculation of Economic-injury Levels and Economic Thresholds for Pest Management.  Journal of Economic Entomology, 80(2):297–303.
  * Otto and Day, (2011) Otto, S. P. and Day, T. (2011).  A biologist’s guide to mathematical modeling in ecology and evolution.  In A Biologist’s Guide to Mathematical Modeling in Ecology and Evolution. Princeton University Press.
  * Phophi et al., (2019) Phophi, M., Mafongoya, P., and Lottering, S. (2019).  Perceptions of climate change and drivers of insect pest outbreaks in vegetable crops in limpopo province of south africa.  Climate, 8.
  * Ramazi et al., (2021) Ramazi, P., Kunegel-Lion, M., Greiner, R., and Lewis, M. A. (2021).  Predicting insect outbreaks using machine learning: A mountain pine beetle case study.  Ecology and evolution, 11(19):13014–13028.
  * Ross, (1998) Ross, S. M. (1998).  A First Course in Probability.  Prentice Hall, Upper Saddle River, N.J., fifth edition.
  * Santos et al., (2017) Santos, S., Specht, A., Carneiro, E., Paula-Moraes, S., and Casagrande, M. (2017).  Interseasonal variation of chrysodeixis includens (walker, [1858]) (lepidoptera: Noctuidae) populations in the brazilian savanna.  Revista Brasileira de Entomologia, 61.
  * Sharma and Dhillon, (2020) Sharma, H. C. and Dhillon, M. K. (2020).  Climate change effects on arthropod diversity and its implications for pest management and sustainable crop production.  Agroclimatology: Linking Agriculture to Climate, 60:595–619.
  * Stern et al., (1959) Stern, V., Smith, R., van den Bosch, R., Hagen, K., et al. (1959).  The integration of chemical and biological control of the spotted alfalfa aphid: the integrated control concept.  Hilgardia, 29(2):81–101.
  * Storn and Price, (1997) Storn, R. and Price, K. (1997).  Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces.  Journal of global optimization, 11(4):341–359.
  * Tao et al., (2012) Tao, J., Chen, M., Zong, S., and Luo, Y.-Q. (2012).  Genetic structure in the seabuckthorn carpenter moth (holcocerus hippophaecolus) in china: The role of outbreak events, geographical and host factors.  PloS one, 7:e30544.
  * Tinsley et al., (2013) Tinsley, N., Estes, R., and Gray, M. (2013).  Validation of a nested error component model to estimate damage caused by corn rootworm larvae.  Journal of Applied Entomology, 137.
  * Tsallis, (1988) Tsallis, C. (1988).  Possible generalization of boltzmann-gibbs statistics.  Journal of statistical physics, 52(1):479–487.
  * Tsallis and Stariolo, (1996) Tsallis, C. and Stariolo, D. A. (1996).  Generalized simulated annealing.  Physica A: Statistical Mechanics and its Applications, 233(1-2):395–406.
  * Van Rossum and Drake Jr, (1995) Van Rossum, G. and Drake Jr, F. L. (1995).  Python reference manual.  Centrum voor Wiskunde en Informatica Amsterdam.
  * Volney and Fleming, (2000) Volney, W. and Fleming, R. (2000).  Climate change and impacts of boreal forest insects.  Agriculture, Ecosystems and Environment, 82:283–294.
  * Wallner, (1987) Wallner, W. (1987).  Factors affecting insect population dynamics: Differences between outbreak and non-outbreak species.  Annual Review of Entomology, 32:317–340.
  * Xiang and Gong, (2000) Xiang, Y. and Gong, X. (2000).  Efficiency of generalized simulated annealing.  Physical Review E, 62(3):4473.
  * Xiang et al., (2013) Xiang, Y., Gubian, S., Suomela, B., and Hoeng, J. (2013).  Generalized simulated annealing for global optimization: the gensa package.  R J., 5(1):13.
  * Xiang et al., (1997) Xiang, Y., Sun, D., Fan, W., and Gong, X. (1997).  Generalized simulated annealing algorithm and its application to the thomson model.  Physics Letters A, 233(3):216–220.
  * Zanuncio et al., (2006) Zanuncio, T., Zanuncio, J., Freitas, F., Pratissoli, D., Sediyama, C., and Maffia, V. (2006).  Main lepidopteran pest species from an eucalyptus plantation in minas gerais, brazil.  Revista de biología tropical, 54:553–60.

## Supplementary Materials: Pattern-Based Prediction of Population Outbreaks

We introduce the pypbp package, which is a Python implementation of the
Pattern-Based Prediction (PBP) method, and examples of how to use it in
Section 4. Also, we present the clustering algorithm used by the PBP method in
Section 5. Finally, we describe the methods used for the PBP optimisation
procedure, including a comparison between the performance of Generalised
Simulated Annealing (GSA) and differential evolution in Section 6.

## 4 The pypbp package

To use our package, first install Python 3 using the official website
(https://www.python.org). We recommend the Jupyter Lab environment
(https://jupyter.org) as a GUI for Python, however there are many other
options. Using the pip command in your terminal (cmd in Windows, or terminal
in Mac and Linux operating systems), you may install Jupyter Lab by executing

⬇

1pip install jupyterlab

To install the pypbp package, execute

⬇

1pip install pypbp

You may then open a Jupyter Lab environment using the command

⬇

1jupyter lab \--core-mode

Finally, you may create a new Jupyter notebook file using the Jupyter Lab
environment. Using the cells of the Jupyter notebook file, you can import and
use functions implemented in the pypbp package, as presented below. More
information is available in the package description page (https://pypbp-
documentation.readthedocs.io/en/latest/) and GitHub repository
(https://github.com/GabrielRPalma/PyPBP).

### 4.1 A minimal reproducible example

The following code obtains estimates for the hyperparameters in the PBP method
based on the aphid data used as motivation in the paper.

⬇

1import pypbp as pbp

2

3results, xstar, clustered_patterns, parameters = pbp.pbp_fit(time_series =
pbp.time_series,

4 train_percentage = 0.5,

5 xstar = 200,

6 maxfun = 200)

7pbp.pbp_plot(time_series = pbp.time_series, clustered_patterns =
clustered_patterns,

8 parameters = parameters,

9 xnew = [10, 60, 20, 10, 20, 80, 90])

By default, the optimisation is performed using 5-fold cross-validation by
carrying out the GSA method with the negative area under the ROC curve as the
objective function to be minimised, without pre-processing the time series
using empirical mode decomposition. The function pbp_fit contains other
arguments, such as verbose, which prints the area under the ROC curve at every
iteration of the optimisation process, and maxfun, which sets the number of
evaluations of the objective function by the optimisation algorithm (GSA as
default).

The object clustered_patterns contains the cluster matrices estimated using
the time series presented, the object results contains a data frame with the
metrics obtained from the model: accuracy, fI-Score, precision, recall, true
positive rate (TPR), false positive rate (FPR) and the estimated area under
the ROC curve. It also contains the estimates for $m$,
$d^{*}_{\mbox{\scriptsize cluster}}$ and $\alpha$. These metrics are presented
for each of the four criteria used in our methodology to select
$d^{*}_{\mbox{\scriptsize base}}$ (maximum FPR of 0.1 and 0.2 and minimum TPR
of 0.8 and 0.9). The parameters object contains the estimates for $m$,
$d^{*}_{\mbox{\scriptsize cluster}}$ and $\alpha$ in a dictionary ready for
use. The function pbp_plot can be used to visualise the obtained patterns.
Finally, the details of each function used in the pypbp package are presented
on the website https://pypbp-documentation.readthedocs.io/en/latest/.

## 5 Pattern clustering algorithm

Algorithm 5 starts with pattern $\textbf{p}_{1}$, which represents the first
row of the matrix P. We remove $\textbf{p}_{1}$ from P and add it as the first
row of $\textbf{P}^{\prime}_{1}$. After that we compute the association metric
between $\textbf{p}_{1}$ and all subsequent rows of P. If
$d(\textbf{p}_{1},\textbf{p}_{j})\geq d^{*}_{\mbox{\scriptsize cluster}}$, we
add pattern $\textbf{p}_{j}$ as the last row of the cluster matrix
$\textbf{P}^{\prime}_{1}$ and delete it from P. We repeat this process to
obtain the cluster matrices $\textbf{P}^{\prime}_{c},c=1,\ldots,C\leq I,$
until there are no more rows left in P.

1Input: $d^{*}_{\mbox{\scriptsize cluster}}$ and P

2 set $c=1$

3 set $l=$ number of rows of $\mathbf{P}$

4 for  _$i$ in $\\{1,2,\ldots,l\\}$_ do

5 add $\textbf{p}_{i}$ as the first row of $\mathbf{P}_{c}^{\prime}$ and
delete it from $\mathbf{P}$

6 set $l^{\prime}=$ number of rows of $\textbf{P}^{\prime}_{c}$

7 for  _$j$ in $\\{1,2,\ldots,(I-\sum\limits_{c}l^{\prime}_{c})\\}$_ do

8 if  _$d(\mathbf{p}_{i},\mathbf{p}_{j})\geq d^{*}_{\mbox{\scriptsize
cluster}}$_ then

9 append $\textbf{p}_{j}$ as the last row of $\textbf{P}^{\prime}_{c}$ and
delete $\textbf{p}_{j}$ from P

10 end if

11 update $l^{\prime}=$ number of rows of $\textbf{P}^{\prime}_{c}$

12 end for

13 update $c=c+1$

14 update $l=$ number of remaining rows of $\mathbf{P}$

15 end for

Output: cluster matrices $\textbf{P}^{\prime}_{c},c=1,\ldots,C\leq I$

Algorithm 1 Obtaining cluster matrices $\textbf{P}^{\prime}_{c},c=1,\ldots,C$,
from $\mathbf{P}$.

## 6 Comparison between optimisation methods

Here, we compare the differential evolution (Storn and Price,, 1997) and the
method (Tsallis and Stariolo,, 1996) algorithms used to estimate $\alpha$,
$d^{*}_{\mbox{\scriptsize cluster}}$ and $m$. We use these methods to optimise
the area under the ROC curve (AUROC). We looked at the performance of our
method using the time series data on aphids and parasitoids collected in
Coxilha (Brazil-São Paulo). In Figure 8 we present the accuracy, in Figure 9
the area below the ROC curve (AUROC), in Figures 10 and 11 we present the true
and false positive rates. In Figures 12, 13, and 14 we present the estimated
$\alpha$, $m$ and $d^{*}_{\mbox{\scriptsize cluster}}$ hyperparameters.

Differential evolution is an algorithm based on evolutionary computation,
which briefly consists of setting a population of candidate solutions and
creating new updated candidate solutions based on the existing based on the
score. For more detail see Storn and Price, (1997). On the other hand, the
generalised simulated annealing is a stochastic algorithm used for finding the
global minimum of a given function in a continuous $D$-dimensional space. The
method utilises a generalised entropic form

$S_{q}=k\frac{1-\sum_{i}p^{q}_{i}}{q-1},$

where $q\in\mathbb{R}$, $p_{i}$ are probabilities of the microscopic
configurations and $k$ is a conventional positive constant (Tsallis and
Stariolo,, 1996). Based on this entropy Tsallis and Stariolo, (1996)
generalised its formulation to include the cases of the Boltzmann Cauchy
machines, allowing to find local minima according to an acceptance temperature
set in the generalised metropolis algorithm. Differential evolution was
approximately ten times slower than generalised method. This fact drove our
decision to choose the algorithm as the default algorithm for optimising the
hyperparameters involved in the PBP method. However, users are allowed to
choose which method to use when utilising the pypbp package.

Figure 8: Accuracy obtained on the testing step of the PBP method. In both
algorithms the following methods to choose $d^{*}_{\mbox{\scriptsize base}}$
were used: based on minimum FPR (.1 and .2 as minimum) and TFR maximum (.8 and
.9 as maximum). Figure 9: Area below the ROC curve obtained on the training
step of the PBP method. In both algorithms the following methods to choose
$d^{*}_{\mbox{\scriptsize base}}$ were used: based on minimum FPR (0.1 and 0.2
as minimum) and TFR maximum (0.8 and 0.9 as maximum). Figure 10: True Positive
Rate obtained on the testing step of the PBP method. In both algorithms the
following methods to choose $d^{*}_{\mbox{\scriptsize base}}$ were used: based
on minimum FPR (0.1 and 0.2 as minimum) and TFR maximum (0.8 and 0.9 as
maximum). Figure 11: False Positive Rate obtained on the testing step of the
PBP method. In both algorithms the following methods to choose
$d^{*}_{\mbox{\scriptsize base}}$ were used: based on minimum FPR (0.1 and 0.2
as minimum) and TFR maximum (0.8 and 0.9 as maximum). Figure 12: Estimated
parameter $\alpha$ obtained on the training step of the PBP method. In both
algorithms the following methods to choose $d^{*}_{\mbox{\scriptsize base}}$
were used: based on minimum FPR (0.1 and 0.2 as minimum) and TFR maximum (0.8
and 0.9 as maximum). Figure 13: Estimated parameter $m$ obtained on the
training step of the PBP method. In both algorithms the following methods to
choose $d^{*}_{\mbox{\scriptsize base}}$ were used: based on minimum FPR (0.1
and 0.2 as minimum) and TFR maximum (0.8 and 0.9 as maximum). Figure 14:
Estimated parameter $d^{*}_{\mbox{\scriptsize cluster}}$ obtained on the
training step of the PBP method. In both algorithms the following methods to
choose $d^{*}_{\mbox{\scriptsize base}}$ were used: based on minimum FPR (0.1
and 0.2 as minimum) and TFR maximum (0.8 and 0.9 as maximum).