3.1. Earth observations

We leverage biophysical and biochemical properties such as Leaf Area Index (LAI), Normalized Difference Vegetation Index (NDVI), chlorophyll content, as well as data on evapotranspiration and soil moisture. These factors play a crucial role in monitoring pest population dynamics. The data is derived from the Sentinel-1/2 and Terra/Aqua (MODIS) satellite missions that provide open access to optical multi-spectral and Synthetic Aperture Radar (SAR) images.

3.2. Terrain & soil characteristics

We incorporate data from open-access digital elevation models and information on topsoil physical properties and soil organic carbon content (de Brogniez et al., Reference de Brogniez, Ballabio, Stevens, Jones, Montanarella and vanWesemael2015; Ballabio et al., Reference Ballabio, Panagos and Monatanarella2016). This allows us to include fixed or long-term characteristics specific to the area of interest.

3.3. Numerical weather predictions (NWP) and reanalysis of environmental datasets

Any high spatial resolution weather forecast can be used. We utilize a custom configuration of WRF-ARW (Skamarock et al., Reference Skamarock, Klemp, Dudhia, Gill, Liu, Berner, Wang, Powers, Duda and Barker2019) at a spatial resolution of 2 km. Hourly predictions are made, and for each trap location (i.e., where we have measurements about pest abundance), we obtain daily values for air (2 m) and soil temperature (0 m), relative humidity (RH), accumulated precipitation (AP), dew point (DP), and wind speed (WS). These parameters have been widely used in related work and are extremely valuable for learning from past (reanalysis) and future (NWP) pest states.

3.4. In-field measurements

In-field measurements involve ground observations of pest abundance using pheromone traps specifically designed for monitoring the cotton bollworm, known by the scientific name Helicoverpa armigera (H. armigera). These traps contain the active ingredients Z-11-hexadecen-1-al and Z-9-hexadecenal. The traps are used from the beginning of the first generation until the end of the season, with regular replacement every 4 to 6 weeks. The company Corteva Agriscience Hellas has established a dense (in time and space) trap network (Figure 2) that covers almost all areas in the Greek mainland where cotton is cultivated. The traps are strategically positioned at suitable distances from each other to prevent interference and ensure accurate data collection. An agronomist examines the traps and counts the trapped insects at regular intervals every 3–5 days. Corteva Agriscience Hellas provides historical data consisting of 398 trap sequences and 8202 unique data points from 2019 to 2022 (Table 1). They also provide auxiliary data on pesticide application, potential crop damage from pests, the severity of the damage, trap replacements, and scouter comments.

Figure 2. Traps distribution in the Greek mainland for 2019–2022. Colors indicate the different agroclimatic zones in which traps from the dataset belong. These zones have been identified based on the study conducted by Ceglar et al. (Reference Ceglar, Zampieri, Toreti and Dentener2019).

Table 1. Summary of trap data

4.1. Causal graph for representing domain knowledge

We constructed a causal graph (Figure 3) based on domain knowledge and expertise, denoted as $ G $ , that represents the underlying causal relationships within the pest-farm ecosystem for the H. armigera case. The graph $ G $ comprises vertices $ V $ , which represent the variables in the system, and directed edges $ E $ , which symbolize the cause-and-effect relationships between these variables. Besides helping us articulate domain knowledge, the causal graph $ G $ will benefit the downstream technical analyses in various ways. For instance, $ G $ will be employed for effect identification via graphical tests (Pearl, Reference Pearl2009), where the structure of $ G $ is integral to discerning causal relationships. Conversely, in the case of estimating conditional average treatment effects within the potential outcomes framework, $ G $ will be utilized as a conceptual guide for considering causal structures during the control phase. In invariant causal prediction, the graph will facilitate the construction of an accurate list of invariant features using causal parents of the target outcome. Moreover, the structural knowledge captured in $ G $ could benefit invariant learning methods by guiding the environment $ E $ definition. This diverse and tailored incorporation of $ G $ is aimed at optimizing the utilization of domain knowledge by the specifications and objectives of each analytical technique.

Figure 3. Causal graph of a pest-farm ecosystem for Helicoverpa armigera case.

Specifically, in the current case of the pest-farm ecosystem of H. armigera, various biotic and abiotic factors (Table 2) can influence the population dynamics $ Y $ of H. armigera (Sharma et al., Reference Sharma, Kumar, Vyas, Sharma and Shrivastava2012). Temperature $ T $ plays a crucial role, affecting the insect’s growth, development, fecundity, and survival (Howe, Reference Howe1967). The size $ SG $ of the first generation is related to the size of the second generation, and the Southern Oscillation Index $ SOI $ has a significant correlation with the size of the first spring generation (Maelzer and Zalucki, Reference Maelzer and Zalucki1999, Reference Maelzer and Zalucki2000). Additionally, the life cycle $ LC $ of H. armigera is temperature-dependent, with completion occurring between 17.5°C and 32.5°C (Mironidis and Savopoulou-Soultani, Reference Mironidis and Savopoulou-Soultani2014). The presence of parasitoids and natural enemies in cotton cultivation is crucial to many IPM programs, including the control of H. armigera (Pereira et al., Reference Pereira, Reigada, Diniz and Parra2019). Many egg parasitoids of different families are known for their high parasitism $ P $ rates and their effectiveness in reducing the population of H. armigera (Noor-ul-Ane et al., Reference Noor-ul-Ane, Arif, Gogi and Khan2015). Nevertheless, parasitism rates are influenced by temperature and relative humidity (Kalyebi et al., Reference Kalyebi, Sithanantham, Overholt, Hassan and Mueke2005; Noor-ul-Ane et al., Reference Noor-ul-Ane, Arif, Gogi and Khan2015). Moreover, the efficacy of spray application $ Sp $ also impacts population dynamics (Wardhaugh et al., Reference Wardhaugh, Room and Greenup1980). The efficacy of $ Sp $ is significantly influenced by the plant growth stage $ PGS $ . During the seedling stage, limited leaf surface area reduces spray coverage, while the vegetative stage offers more extensive leaf area, enhancing spray interception. However, dense canopies at later stages may impede spray penetration. Plant physiology also varies, affecting the absorption and translocation of sprayed substances (Fishel and Ferrell, Reference Fishel and Ferrell2010).

Table 2. Pest-farm ecosystem variables

Other environmental factors come into play as well. Precipitation $ \mathit{\Pr} $ affects the population size, with heavy precipitation leading to a decrease in the population (Ge et al., Reference Ge, Liu, Ding, Wang and Zhao2003). It also increases soil water content $ SW $ which affects the emergence rate of H. armigera similar to air relative humidity $ RHa $ (Fajun et al., Reference Fajun, Baoping and Xiaoxi2003). The presence of fruiting organs during the plant growth stage $ PGS $ is important for population dynamics, as it serves as the oviposition site for females (Fitt, Reference Fitt1989). Crop variety $ V $ , such as transgenic Bt cotton, can suppress the second generation of H. armigera, while both different cropping systems $ CS $ and adjacent crops $ AC $ can influence the population structure (Wardhaugh et al., Reference Wardhaugh, Room and Greenup1980; Gao et al., Reference Gao and Zhai2010; Lu et al., Reference Lu, Zalucki, Perkins, Wang and Wu2013). Finally, wind $ W $ and wind direction play a significant role in the emergence of H. armigera, influencing the distance covered during migration from nearby locations. Additionally, wind conditions at the time of spraying $ {W}_s $ can also impact the effectiveness of the intervention. These various factors collectively shape the population dynamics of H. armigera in a complex and interconnected manner as defined through domain knowledge and depicted in the causal graph (Figure 3).

4.2. Invariant & causal learning for robust pest prediction

Our goal is to predict near-future pest populations ( $ {Y}_{t+1} $ ) using Earth observation (EO) and environmental data ( $ {X}_t $ ) along with weather forecasts ( $ {W}_{t+1} $ ) by learning the function $ {y}_{t+1}=f\left({x}_t,{w}_{t+1}\right) $ . Pest management recommendations heavily depend on these predictions. Conventional machine learning methods (Aparecido et al., Reference Aparecido, Rolim, Moraes, Costa and Souza2019; Skawsang et al., Reference Skawsang, Nagai, Tripathi and Soni2019; Xiao et al., Reference Xiao, Li, Kai, Chen, Zhang and Wang2019; Zhang et al., Reference Zhang, Huang, Pu, González-Moreno, Yuan, Wu and Huang2019), which often assume that data points are independent and identically distributed (i.i.d.), struggle to generalize to unseen environments, capture spatiotemporal variability, and adapt to climate change. These methods are prone to learning spurious correlations, limiting their effectiveness in dynamic and non-i.i.d. scenarios.

To address these challenges, we turn to causal learning (Schölkopf and von Kügelgen, Reference Schölkopf and von Kügelgen2022), which leverages domain knowledge and is grounded in the principle of independent causal mechanisms. This principle suggests that joint probabilities can be decomposed into separate mechanisms, each reflecting an underlying causal relationship that remains stable despite environmental changes. By incorporating this principle, our models can improve generalization and robustness across varying conditions.

We achieve this by integrating invariant learning with causality and categorizing dataset units into environments $ E $ as different agroclimatic zones or host crops (Figure 4). While $ E $ influences feature $ {x}_t,{w}_{t+1} $ , it does not directly affect the target $ {Y}_t $ . Utilizing Invariant Causal Prediction (ICP) (Heinze-Deml et al., Reference Heinze-Deml, Peters and Meinshausen2018), Directed Acyclic Graphs (DAGs), and Invariant Risk Minimization (IRM) (Arjovsky et al., Reference Arjovsky, Bottou, Gulrajani and Lopez-Paz2019), we can select causal features, identify potential causal relationships, and capture latent causal structures. These tools allow us to build models that are effective in current conditions and adaptable to future environmental changes.

Figure 4. Invariant learning for robust predictions. Stable and accurate predictions in diverse environments, such as when H. armigera feeds on different crops exhibiting variations in phenotype, agricultural management practices, and spatial distribution. Traditional ML methods risk capturing spurious correlations, such as associating pest abundance with a specific crop (e.g., cotton) due to its higher frequency in the dataset, leading to biased predictions based on the underlying crop rather than true pest presence.

4.3. Explainability & counterfactual reasoning for short-term advice

We define the problem as a binary classification of pest presence or absence at the next time step, using Earth observation (EO) data ( $ {X}_t $ ) and weather forecasts ( $ {W}_{t+1} $ ). The goal is to predict the pest population value at the next time step, $ {Y}_{t+1} $ , by learning the function $ {y}_{t+1}=f\left({x}_t,{w}_{t+1}\right) $ . To enhance the trustworthiness of our predictions, we employ Explainable Boosting Machines (EBM) (Nori et al., Reference Nori, Jenkins, Koch and Caruana2019). This glass-box model achieves high performance while providing inherent explanations at both global and local levels. EBM’s additive nature allows for the sorting and visualization of feature contributions on a local scale for each one of predictions and a global level to summarize the general behavior of the model depending on features (Figure 5), which facilitates a better understanding of the primary drivers of the model and enhances trust in its outputs.

Figure 5. Explainability for trustworthiness enhancement, on the right, with local and global explanations of each prediction and general model behavior, respectively, & Counterfactual explanations as agricultural actionable recommendations on the left.

We propose generating counterfactual examples as recommended interventions to bolster trust further and provide actionable insights. Following the setup of (Mothilal et al., Reference Mothilal, Sharma and Tan2020), we search for minimal perturbations to the feature values $ \left({x}_t,{w}_{t+1}\right) $ that would alter the prediction to the desired class using the same model $ f $ . These counterfactual examples represent proposed actions that could be implemented in natural farm systems, ensuring practicality and feasibility (Wachter et al., Reference Wachter, Mittelstadt and Russell2017; Mothilal et al., Reference Mothilal, Sharma and Tan2020). The approach ensures that the generated counterfactuals are close to the original input but predicted in the desired class, providing feasible and actionable recommendations for IPM (Figure 5).

4.4. Heterogeneous treatment effects for long-term advice

We provide long-term pest prevention and suppression advice by assessing how agricultural practices (e.g., crop rotation, balanced fertilization, sowing dates) impact pest harmfulness and yield indices. Since different agro-environments may respond variably to the same practice, it is crucial to account for this heterogeneity. We estimate the conditional average treatment effect (CATE) following the potential outcomes framework (Rubin, Reference Rubin2005).

The CATE quantifies the difference in potential outcomes, represented as $ \unicode{x1D53C}\left[Y\left(T=1\right)-Y\left(T=0\right)|X\right] $ , where $ Y(T) $ denotes the value of a random variable $ Y $ (e.g., pest harmfulness and yield) if a unit is treated with treatment $ T\in \left\{0,1\right\} $ . By controlling for field characteristics $ X $ —which capture the heterogeneity across different agro-environmental conditions—we can better understand how specific practices affect outcomes in various contexts (Figure 6). This approach allows us to provide tailored and effective long-term IPM advice sensitive to each field’s unique conditions (Giannarakis et al., Reference Giannarakis, Sitokonstantinou, Lorilla and Kontoes2022).

Figure 6. Conditional Average Treatment Effect (CATE) is seen as long-term personalized guidance. By accounting for each land unit’s unique characteristics, we can estimate a distinct treatment effect for each land unit. For example, how differences in land’s characteristics can change the impact of fertilizer application on increasing the risk of pest emergence in the future.

4.5. Causal inference for evaluating advice effectiveness

We employ causal inference techniques to assess the effectiveness of our pest control recommendations, adapting approaches recently introduced in agricultural contexts (Tsoumas et al., Reference Tsoumas, Giannarakis, Sitokonstantinou, Koukos, Loka, Bartsotas, Kontoes and Athanasiadis2023). Specifically, in the case of pest management and with available panel data (Table 1), we utilize causal models such as difference-in-differences (DiDs) (Abadie, Reference Abadie2005), synthetic control (Arkhangelsky et al., Reference Arkhangelsky, Athey, Hirshberg, Imbens and Wager2021) and synthetic DiDs (Abadie, Reference Abadie2021) to quantify the treatment effect of adhering to our framework’s recommendations (treated units) compared to those who did not (control units). Historical intervention data retrospectively annotated based on whether our framework recommended action, will serve as the basis for advice evaluation. Causal inference will be performed per-environment to ensure comparability between treatment and control groups, adhering to the parallel trends assumption (Lechner et al., Reference Lechner2011).

However, digital agriculture requires a two-level evaluation of interventions to disentangle the effectiveness resulting from the accuracy of the recommendation (for intervention) in terms of space–time from the inherent efficacy of the intervention. It is crucial to determine what effect, if any, is attributable to the space and time of application and what is due to the pesticide itself.

In this context, we conducted an initial analysis using the aforementioned panel data to quantify the impact of pesticide application on pest abundance in a real-world setting without expert system guidance, employing staggered DiDs with fixed effects (Eq. 4.1).

Table 3. Results of staggered DiDs with controls for unobserved heterogeneity at the unit and time levels by including fixed effects

Note: It includes point estimates, 95% confidence intervals, and p-value. Numbers represent the increase/decrease of accumulated pest catchments at the trap level after the intervention.

The staggered approach accounts for units receiving treatment at different periods. We include unit-fixed effects to control for each unit’s time-invariant characteristics and time-fixed effects to capture overall time trends that affect all units in each period. The unit of analysis is the plot where the pest trap is located, with periods modeled at the weekly level. Here, $ {Y}_{it} $ represents the outcome variable, accumulated pest abundance, for each unit $ i $ at the time $ t $ , and $ \mathrm{treated}\_{\mathrm{time}}_{it} $ is an indicator of whether the unit $ i $ receives treatment (pesticide application) in a period $ t $ (in a staggered manner across units). Specifically, $ {\beta}_0 $ is the intercept, $ {\beta}_1 $ is the treatment effect coefficient, $ {\alpha}_i $ represents unit fixed effects, $ {\gamma}_t $ captures time fixed effects, and $ {\varepsilon}_{it} $ is the error term. Thus, $ {\beta}_1 $ provides the average causal effect of the treatment (pesticide application) on the outcome (accumulated pest abundance) for treated units (ATT), as presented in Table 3 for each cultivation period from 2019 to 2022.

(4.1) $$ {Y}_{it}={\beta}_0+{\beta}_1\cdot \mathrm{treated}\_{\mathrm{time}}_{it}+{\alpha}_i+{\gamma}_t+{\unicode{x025B}}_{it} $$

For the years 2021 and 2022, we observe a statistically significant reduction in pest abundance, while for 2019 and 2020, we find the opposite effect. At first glance, this contradiction may seem unusual, but several reasonable explanations could account for it. Since the data come from real-world agricultural practice, it likely encapsulates some of the following issues: (i) Some interventions may have been applied incorrectly regarding timing and method, reducing or eliminating their efficacy in the pest-infested plots. This could lead to a biased estimate that the pest population increased after pesticide application (Figure 7). This occurs because the counterfactual is constructed by taking the growth trend from a plot without intervention, which might not experience the same infestation or pest pressure level. So, a mistreated plot that probably follows a steeper population increase, simply due to its higher infestation levels, can lead to this fallacy that pesticide application increases pest population. (ii) After discussions with the data provider (Corteva Agriscience Hellas), noise within the control group labels may be possible. The company is confident in the labels for treated plots, as they receive this information directly from farmers. However, they cannot be as certain about the control group. Some farmers may have applied pest control practices in their plots but chose not to report them for various reasons, such as using less expensive pesticides from competitor companies or participating in eco-schemes prohibiting pesticide use. Consequently, we face a scenario of positively labeled and unlabeled data, a common issue in machine learning. (iii) The assumption of parallel trends may not hold universally, or unobserved confounders may vary over time and between units.

Figure 7. A visual example of DiDS for assessing the real-world impact of pesticide application. It demonstrates how, even when the parallel trends assumption holds in both conditions, applying an intervention (i.e., spray) at a non-recommended time can lead to unexpected effects compared to applying the intervention at the recommended time.

In a more robust causal analysis, we can technically or conceptually address these issues. Technically, we could retrospectively employ a recommendation system or consult experts, as aforementioned, to annotate each time–space slot as favorable or unfavorable for intervention. On the other hand, we can conceptually accept reality and precisely define what causal effect we retrieve. In this case, the ATT in a real-world setting includes different application accuracy levels, farmer’s skills, expert guidance, and proper timing. To address the second issue, we plan to use Positive-Unlabeled (PU) learning methods (Bekker and Davis, Reference Bekker and Davis2020) to train a classifier on covariates, as they are outlined in Section 3. Using the positively labeled (treated) units only as ground truth and PU learning for training, this classifier will help establish a control group consisting only of unlabeled units that are classified there with high confidence. Lastly, a formal investigation with statistical tests is required to retain only cases where the parallel trends assumption holds. Clear assumptions statements should also be made regarding the potential of unobserved confounders that may vary by time and unit. By leveraging these techniques, we aim to rigorously evaluate the impact of our recommendations on pest control outcomes and attribute the effects to the right factors, providing robust evidence for the effectiveness of our framework in diverse agricultural environments.