The management of agricultural fields (i.e. farming) in the era of data science has evolved to use spatial mapping, remote sensing, soil and terrain measurements, weather measurements and other data sources to improve the quantity and quality of crops [21, 9]. The management of agricultural fields using advanced data analytics is referred to, collectively, as “precision agriculture” [5, 8] and is quickly becoming industry standard. Broadly, precision agriculture attempts to use the spatial variability within a field to manage individual crop locations rather than treat a field as spatially and temporally homogeneous. The spatial variability within a field can be influenced by differences in elevation, aspect, or other topographical features of the field that influence both irrigation and crop yield. Given that agriculture fields often occupy multiple acres of space, precision agriculture has been shown to outperform basic farming techniques by increasing crop yield [30, 32].

Crop yield is heavily driven by the volumetric water content (VWC; the ratio of the volume of water currently in the soil to the unit volume of soil, or in other words the quantity of water contained in the soil) and, in arid regions, irrigation is a practiced method for controlling and adjusting the VWC. As such, variable rate irrigation (VRI) is a practice within precision agriculture focused on using data to adjust the amount of water applied throughout the field according to spatial and temporal variations [20, 39, 40, 25, 60]. Such management zones can be temporally dynamic to help inform day-to-day watering decisions [13] or, as is the case in this application, static to inform season-long watering decisions.

One approach to VRI is to partition the field into management zones (or irrigation zones) wherein irrigation rates are homogeneous within each management zone but different across zones rather than utilizing a constant irrigation rate throughout the entire field [18, 61, 57, 38, 37]. Benefits of using VRI with management zones include reductions in water and energy input, increased crop productivity, decreased runoff (i.e. less waste) and reductions in chemical inputs and soil pollution [42, 3, 25].

Effective irrigation zones should partition the agricultural field based on VWC [31]. However, obtaining VWC is a labor intensive process and is often only sparsely measured for the entire field [27]. As an example, consider Figure 1 which shows 66 VWC measurements averaged over four time periods between April and September 2019 for an agricultural field of winter wheat in Rexburg, Idaho along with the associated irrigation zone for each location. Although the VWC was recorded at 66 different locations, the available data do not adequately cover the field to the point where precise irrigation decisions and zones can be precisely defined. That is, based on such sparse measurement the boundaries between zones are unclear.

Rather than using the VWC directly, agricultural scientists have begun to use alternative, more easily obtained, data to inform irrigation zones [50, 56, 26, 24]. Examples of such possible data for the Rexburg field considered in this research are shown in Figure 2 and include historical yields, normalized difference vegetation index (NDVI), elevation, topographical wetness index (TWI), aspect, and slope (more details on these covariates are given in the application Section 3.1 below). These alternative covariates to VWC were recorded for over 5000 unique locations in the field using simple remote sensing mounted on drones or tractors making these data readily available to farmers.

The primary issue of using alternative covariates to VWC such as that in Figure 2 to define VRI zones is that such covariates are non-linearly related to VWC and, hence, may result in less water efficient zones than ones defined directly by VWC. Figure 3 shows scatterplots of a few of the variables against the VWC for the Rexburg field demonstrating clear non-linear relationships. In fact, the scatterplots in Figure 3 show a highly complex relationship between these alternative covariates and VWC suggesting that zones defined directly on the covariates may vary considerably from the more water-efficient zones defined on VWC.

In this research, we seek to implement a method to delineate irrigation zones based on VWC using the data displayed in Figure 2 as covariate information. Specifically, we define a model relating the easily obtained and spatially dense covariate data to the sparse VWC data. To do so, we integrate deep learning into a statistical modeling framework for irrigation zones to capture the complex relationship between the covariate data and the response variable.

Deep learning is a subfield of machine learning focused on using highly flexible algorithms called neural networks to learn complex relationships between covariates and response variables [53]. Deep learning algorithms have been successfully implemented across a wide range of applications leading to their explosion in popularity in recent years [1, 41, 49, 12]. More recently, deep learning (and machine learning more generally) have successfully incorporated spatial correlation [48, 45] into these algorithms.

From a statistical perspective, deep learning algorithms do not inherently account for uncertainty in the predictions and are often so complex that they are uninterpretable (so-called “black box” methods). Yet, recent efforts in the statistics community have started to fold deep learning into statistical models [34, 54, 46, 11]. While we acknowledge that some methods (such as dropout) can give uncertainty estimates [see 15], we consider embedding a deep neural network into a statistical model and follow [52] by estimating associated parameters using a fully Bayesian paradigm. Further, we implement Bayesian versions of partial dependence plots to give some interpretability to the deep neural network between the covariates and the response variable.

With a deep neural network accounting for the complex relationships between the covariates and response, our model also creates smooth irrigation zones by incorporating spatial correlation. Specifically, we use spatial basis functions distinctly from the neural network portion of the model thereby successfully merging deep learning into a spatial modeling framework.

The primary quantity of interest in this research is the irrigation zone of each location on the field – not necessarily the VWC. Admittedly, we could predict the VWC first and then create irrigation zones by splitting based on predicted VWC quantiles. However, here we perform modeling for the associated zone as the response of interest given by the right panel of Figure 1. We focus on the zone rather than VMC for a few reasons. First, our data is unique in that we have the exact VWC measurements. However, given the cost, obtaining exact VWC data is exceptionally rare in agricultural practice (our data being an exception rather than the rule). Predicting continuous VWC would be applicable only to this field in Rexburg and would not be able to be tested on other fields. Rather, most agriculture fields are able to only collect ambiguous VWC levels of “low” or “high.” By focusing on irrigation zone, our statistical model will be more aligned with typical data available in other crop fields. Second, VRI technology is typically done by applying “less than average” water to some locations and “more than average” to others. Hence, what is required to implement VRI is knowledge of which locations are “low” and which are “high” rather than the exact VWC.

Given the focus on irrigation zones, the response variable is thus an ordered multinomial response. Under this discretization, we build our model following the latent variable approach of [2], [22] and [4] to facilitate Bayesian computation. While others have developed statistical models for neural networks [see 16, 34], to our knowledge, this is the first attempt at building a deep spatial statistical model for an ordered categorical response.

We view this work as a contribution to the growing body of literature of spatial deep learning methods [see 58]. Alternative approaches include [7] use spatial bases as inputs into neural networks. However, here we use spatial bases separately from the neural network because using such bases as inputs may sacrifice some spatial smoothness which is crucial to our application. Likewise, [59] use generalized squared error loss functions with a spatial covariance to fit spatial graphical neural networks. Our approach is inherently different from this because our response is non-Gaussian so squared error loss functions are not directly applicable.

The remainder of this paper is outlined as follows. Section 2 describes our spatial neural network model along with our Bayesian strategy to estimate model parameters. Section 3 describes the result of the model fit to the Rexburg field including model tuning, predicted zone delineation and effect interpretation. Finally, Section 4 concludes by pinpointing strengths and weaknesses of our approach along with areas for future research.

Let Y ( s ) ∈ { 1 , … , R } denote the irrigation zone for location s = ( s 1 , s 2 ) ′ ∈ R 2 where R is the number of desired irrigation zones for a field. Given the current state of variable rate irrigation systems, usually R will only be as high as 4 or 5. Further, because each irrigation zone is given a different amount of water (i.e. dry irrigation zones receive more water), we can treat Y ( s ) as an ordered, discrete spatial random variable. Due to the difficulty of performing spatial analysis on a discrete scale, we follow [2], [22] and [4] by augmenting Y ( s ) with a latent variable Z ( s ) ∈ R such that (2.1) ( Y ( s ) ∣ { c r } , Z ( s ) ) = ∑ r = 1 R r 𝟙 ( c r − 1 ≤ Z ( s ) < c r ) where 𝟙 ( · ) is an indicator function and c 0 = − ∞ < c 1 = 0 < c 2 < ⋯ < c R = ∞ are cut points used to determine the probability that location s belongs to a specific zone.

Under Equation (2.1) and given the cut points c 0 , … , c R , Y ( s ) is completely determined by Z ( s ) so that modeling Z ( s ) will induce a statistical model for Y ( s ). Using this relationship between Y ( s ) and Z ( s ), we assume (2.2) ( Z ( s ) ∣ w ( s ) ) ∼ N ( f L ( X ( s ) ) + w ( s ) , 1 ) where f L ( X ( s ) ) is the univariate output function of an L-layer feed forward neural network (FFNN) or, alternatively referred to as a multilayer perceptron (MLP) with input covariates X ( s ) = ( X 1 ( s ) , … , X P ( s ) ) ′ and w ( s ) is a spatial random effect. Under this model structure, we follow [2] by fixing the variance at 1 to identify the scale of the weights of the neural networks (coefficients) and we set c 1 = 0 to identify a bias term. This model, importantly, accounts for the various challenges of determining irrigation zones via statistical modeling discussed in Section 1.2. First, the FFNN in (2.2) utilizes a nonlinear relationship between the set of covariates X ( s ) and the associated irrigation zone (see Figure 3). Second, the spatial random effect w ( s ) in (2.2) ensures that the fitted irrigation zones change smoothly over space thereby creating zones that can be reasonably implemented via VRI. Consider each piece of the model in the following paragraphs.

FFNNs have become known for the strong performance in prediction, especially where non-linear relationships exist between covariates and the response variable as is the case here [12]. A FFNN with L layers, consists of three different types of layers: an input layer, hidden layers and an output layer each consisting of a different number of units. More succinctly, let f l ( X ( s ) ) = ( f ℓ 1 ( X ( s ) ) , … , f ℓ P ℓ ( X ( s ) ) ) ′ be the vector of P ℓ units at layer ℓ of the FFNN. The transformation from layer to layer occurs via (2.3) f ℓ ( X ( s ) ) = a ℓ ( λ ( ℓ − 1 ) 0 + Λ ℓ − 1 f ℓ − 1 ( X ( s ) ) ) where a ℓ ( · ) is a element-wise nonlinear activation function used at layer ℓ, λ ( ℓ − 1 ) 0 = ( λ ( ℓ − 1 ) 01 , … , λ ( ℓ − 1 ) 0 P ℓ ) ′ is the vector of intercepts (biases) applied to layer ℓ − 1 and Λ ℓ − 1 = { λ ( ℓ − 1 ) i j } i , j is the P ℓ × P ℓ − 1 matrix of coefficients (weights) used to transition from layer ℓ − 1 to layer ℓ. For the input layer (i.e. ℓ = 1), we take f 1 ( X ( s ) ) = X ( s ) such that P 1 = P is the number of covariates included in the model. From a statistical point of view, the weights and biases (intercepts and coefficients) { λ ℓ 0 , Λ ℓ } ℓ = 1 L − 1 are the model parameters to be estimated from the data while the number of layers L, the layer dimensions { P ℓ } and the activation functions are known as “tuning parameters” and are fixed via cross-validation to optimize prediction performance.

The flexibility of FFNNs to capture complex relationships comes from appropriate choice of L, { P ℓ } and { a ℓ ( · ) }. First, deeper networks (larger L) may be required to allow for more complex transformations of X ( s ) between the input and output layers. Conversely, in the simplest case, when L = 2, the relationship between X ( s ) and Z ( s ) would be linear provided a L ( x ) = x is the identity function.

Beyond the number of layers, the number of units per layer (denoted by P ℓ ) can also be modified to capture complex relationships between the covariates and the response variable. For example, if P ℓ > P 1 for any ℓ > 1 then the FFNN uses a dimension expansion to model the relationship between covariates and the response [55, 6, 47]. As recently shown by [28] and [33] but pioneered by [35], an infinite dimension expansion would mimic a Gaussian process regression model between the covariates and response. Conversely, if P ℓ < P 1 for any ℓ > 1 then the transition in (2.3) acts similar to a principal components regression model by reducing the dimension of the input space [see 14, Chapter 5].

The activation functions { a ℓ ( · ) } l = 2 L in (2.3) are typically non-linear to ensure that the relationship between X ( s ) and Z ( s ) is also non-linear. Typical choices of a ℓ ( · ) include rectified linear units, hyperbolic tangents and identity (see [36] and [44] for a discussion of these activations). However, simple algebraic manipulations show that if a ℓ ( x ) = x (the identity activation) for all ℓ then the relationship between Z ( s ) and X ( s ) would be linear. Notably, to increase flexibility, different layers can utilize different activation functions but, traditionally, the same activation function is used for all layers with the exception of the output layer where the domain of a L ( · ) must match the support of Z ( s ) (which, in our cases, a L ( · ) is the identity function).

While f L ( X ( s ) ) captures the non-linear relationship between the response and covariates, the spatial random effect w ( s ) serves to smooth the predicted zones over space. To achieve this smoothing, we opt to use the Moran basis function expansion for w ( s ) advocated by, among others, [23]. To calculate the Moran basis functions, we first determine an adjacency matrix A = { a i j } between all observed locations and desired prediction locations according to the inverse weighting distance (2.4) a i j = 0 if i = j 1 ‖ s i − s j ‖ otherwise where s i and s j are two locations on the field (either observed or predicted). Let P = I − 11 ′ / n where 1 is a vector of ones and n is the total data size (number of observed locations plus the number of desired prediction locations). Let B be the n × K matrix of the K eigenvectors associated with the K largest positive eigenvalues of P A P. Then, we set (2.5) w ( s ) = b ′ ( s ) β where b ′ ( s ) the row of B associated with location s and β is a vector of coefficients to be estimated from the data. Generally, as K increases then the associated fitted spatial surface from the Moran basis is more variable. Hence, for purposes of this research, we treat K as an additional tuning parameter that we will choose via cross-validation.

Even though neural networks are typically fit via loss minimization, in an effort to merge machine learning and statistical modeling, we adopt a Bayesian approach for parameter estimation to ensure that our predicted irrigation zone surface accounts for associated parameter uncertainty. Accounting for uncertainty here is important so that the predicted irrigation zones can be potentially altered to more easily be incorporated into a VRI system. For example, if the zone assigned to a certain location is highly uncertain, then that location can be manually assigned to a zone by the farmers to increase the overall efficiency of the VRI.

Under the Bayesian approach, prior distributions are required for all the neural network parameters { λ ℓ 0 , λ ℓ }, the Moran basis coefficients ( β ) and the cut points c 2 , … , c R and parameter estimation is done via posterior inference. Because the β parameter vector is simply coefficients in an ordered probit regression model, we assume a vague N ( 0 , 100 I ) prior distribution because the data can generally estimate these parameters well [22]. Because the cut points are ordered so that 0 < c 2 < ⋯ < c R = ∞, each cutpoint was transformed according to c 2 ⋆ = log ( c 2 ) and (2.6) c r ⋆ = log ( c r − c r − 1 ) for r = 3 , … , R − 1 and a N ( 0 , 10 ) prior was used for each c r ⋆ . The corresponding back transformation c r = ∑ i = 2 r exp { c i ⋆ } ensures the ordering constraint.

The priors for the FFNN weights { Λ ℓ } and biases { λ 0 ℓ } need to be chosen with care to avoid overfitting. As documented by, among others, [10], neural networks can easily overfit training data resulting in poor predictive performance. One common approach for restricting neural networks is via penalization (regularization). In a Bayesian setting, regularization is enforced via informative prior constraints [see 43, for a review]. As such, we assume a priori independent N ( 0 , 0.01 ) priors for all biases and weights. The highly informative 0.01 variance constrains the biases and weights to be near zero; thus preventing overfitting similar to ridge and LASSO regression models. For our application below, we tried a few different values for the prior variance here but ultimately settled on 0.01 as the best performing prior in terms of predictive accuracy (see Section 3.1 below). Admittedly, a Laplace prior could be used for the weights because it corresponds to the LASSO penalty [51]. For our purposes, the above Gaussian prior worked adequately to prevent overfitting as demonstrated by our cross-validation study below.

Posterior inference for our model parameters was accomplished by using Markov chain Monte Carlo (MCMC) sampling. Conditional on all other parameters, the complete conditional distribution for Z ( s ) is a truncated Gaussian distribution with mean f L ( X ( s ) ) + b ′ ( s ) β, variance 1 and endpoints c Y ( s ) − 1 and c Y ( s ) . Because each Z ( s ) is conditionally independent, this sampling can be done independently and efficiently. Next, because Z ( s ) ∈ R, we use assume a L ( · ) is the identity activation and P L = 1 so that (2.7) Z ( s ) ∼ i i d N ( λ 0 ( L − 1 ) + f L − 1 ′ ( X ( s ) ) Λ L − 1 + b ′ ( s ) β , 1 ) which can be rewritten simply as Z ( s ) ∼ N ( X ⋆ ′ ( s ) β ⋆ , 1 ) where β ⋆ = ( λ 0 ( L − 1 ) , Λ L − 1 ′ , β ′ ) ′ and X ⋆ ( s ) = ( 1 , f L − 1 ′ ( X ( s ) ) , b ′ ( s ) ) ′ . Under the above Gaussian priors, the complete conditional for β ⋆ is Gaussian and can be sampled directly.

While { Z ( s ) } and { λ 0 ( L − 1 ) , Λ L − 1 , β } can be sampled directly from the corresponding complete conditional distributions, the cut points c 2 , … , c R − 1 and other FFNN parameters { λ 0 l , Λ l } l = 1 L − 2 need to be sampled indirectly via Metropolis-Hastings or other algorithms. In early phases of this research, the weights and biases were sampled individually or a single column at a time. However, this proved to be computationally expensive and there was high correlation between the different weights within each layer. As a result, our final MCMC algorithms sampled the neural net weight matrix Λ ℓ and the biases λ 0 ℓ were sampled jointly for each layer (with each layer being sampled separately). For our MCMC algorithm, we used the adaptive Metropolis algorithm from [17] to update the proposal variance and achieve better mixing for the FFNN parameters. Finally, this adaptive Metropolis algorithm was again used to sample the transformed cutpoints c 2 ⋆ , … , c R − 1 ⋆ after marginalizing out Z ( s ) as suggested by [22].

For our application, we used ten covariates: elevation, yield, NDVI index for 2018 and 2019, two different measures of slope at a location, the x and y aspect of a location and a topographical wetness index (TWI). Each of the covariates were observed at 5062 10 m 2 areas in the field (see Figure 2). In this application, NDVI (a measure of vegetation greenness) was calculated using an aerial drone at four distinct times throughout a growing season and averaged as a single NDVI measurement for the year. Yield was measured per 10 m 2 unit of land using a grain yield monitor installed on the harvester. TWI, is a static measure, calculated from a digital elevation models, and indicates where water will accumulate in an area with elevation variability. Slope as well as the x and y aspects were likewise calculated from a digital elevation model.

The data of irrigation zone consisted of only the 66 observations shown in the right panel of Figure 1. For our implementation, we set the number of irrigation zones, R, to 3. As previously mentioned, R will rarely if ever be higher than 4 or 5 and, in discussion with the farm owner, 3 zones was determined to be reasonable based on the field’s VRI capacities. We also used rectified linear unit activation functions for all the activation functions with the exception of the output layer which was an identity activation to match the support of Z ( s ). Certainly, other activation functions could be used but the rectified linear units is one of the most common.

Beyond the above model settings, implementation of our spatial neural network model requires tuning the number of layers (L), the number of units per layer ( { P ℓ }), the number of Moran basis functions (K) and the penalization for the weights in the prior (the prior variance). Computationally, fitting the neural network model at a single setting took about 3.5 hours on a 2.60 GHz CPU. For each parameter setting in Table 1 and for each of 0.1, 0.01 and 0.001 for the prior variance of the weights, we implemented a 6-fold cross validation and calculated the average adjusted Rand index (ARI). The ARI is given by A R I = ∑ i j n i j 2 − [ ∑ i a i 2 ∑ j b j 2 ] n 2 1 2 [ ∑ i a i 2 + ∑ j b j 2 ] − [ ∑ i a i 2 ∑ j b j 2 ] n 2 where a r is the number locations predicted to belong to irrigation zone r and b r is the number of true locations belonging to irrigation zone r for r = 1 , … , R and n i j is number of location classified as zone i in the predictions and zone j in the truth. ARI, intuitively, is a measure of similarity between two data clusterings (in our case, the similarity between the true zone and the predicted zone) and ranges from −1 to 1 where 1 indicates perfect agreement between the two clusterings, 0 indicates a random agreement and −1 indicates that the two clusterings are completely different. Prior to fitting the models in Table 1, a larger grid was first used to get a general idea of the reasonable values for the tuning parameters. Given our data set consisted of only 66 observations, an effort was made to keep the total number of the parameters low. Deeper neural networks, generally, require big data to be effective. Hence, the grid search only examines one or two layer neural networks (in addition to the input and output layers). The maximum number of Moran basis functions considered was 10 in order to make sure that the zones were being determined by the neural network predictions rather than being overly driven by the spatial aspect of our model.

Table 1 shows the cross validation results using the 0.01 prior variance because this value was uniformly better in terms of predictive accuracy. Cross-validation finds that a single hidden layer with 10 neurons and 5 Moran basis functions were the ideal parameters for this data. The results displayed in the following subsections are from this model setting. Note that, generally, from Table 1, adding spatial basis functions improved the model’s predictive ability suggesting that merging spatial modeling techniques with deep learning is of value in this particular setting.

When using Markov chain Monte Carlo sampling in a Bayesian framework as is the case here, it is important to make sure that the algorithm provides samples of the parameters from posterior distribution via convergence diagnostics. The supplementary material to this article includes trace plots of f L ( X ( s ) ) + w ( s ) for four different locations along with a trace plot of an example cut point. We focus on these trace plots as f L ( X ( s ) ) + w ( s ) and the cut point because these are the main quantities from which we derive our irrigation zone delineation. These trace plots show that these parameters have converged and they can be used to for posterior inference.

As a means to validate the use of the spatial neural network model proposed in this analysis, three additional alternative models were fit to the data and the ARI was again calculated for the 66 locations across the field in Rexburg. The four models compared are ordinal logistic regression with linear effects (Logistic), ordinal logistic regression with natural splines (NSLogistic), ordinal logistic regression with natural splines and spatial bases (Spatial NSLogistic), the neural network model with no spatial basis functions (NN) and our full neural network model with spatial basis functions (Spatial NN). Ordinal logistic regression is a common approach to modeling discrete ordinal multinomial response variables. One of the drawbacks is the assumptions that the relationship between predictors and the response is linear. Splines are a common approach to modeling non-linear relationships for statisticians, while neural networks are often a machine learning approach to account for these non-linear relationships. These two models can provide a comparison between these two different approaches to model a non-linear relationship. While neural networks and splines are able to account for the non-linear relationships, they don’t incorporate any spatial correlation. For irrigation zone delineation, the spatial correlation is used to smooth out the irrigation zone predictions across the field. As a result, we also included the spatial basis functions beyond the natural splines and/or neural network to help improve the predictions of the zones.

Table 2 shows the adjusted Rand index for the four models using all 66 locations in the field. The results show that the ordinal logistic regression model is least able to delineate zones. This is, perhaps, not surprising given this model doesn’t account for the non-linear relationships between predictors and the response variable or the spatial correlation. The spline model and neural network without spatial basis functions had a very similar performance and both seem to be effective approaches to modeling a non-linear relationship. However, this effectiveness was enhanced in both models by the addition of spatial basis functions. The best model was the neural network with the spatial basis functions, which validates the use of the model proposed in this analysis since it accounts for both the non-linear relationships and spatial correlation.

The primary goal of the analysis was to obtain predictions of the irrigation zones for the whole field for use in the VRI system. The fitted spatial neural network model was used to make predictions of which irrigation zone a location should be classified in across the whole field. Due to Monte Carlo sampling, the uncertainty in the zone delineation was also able to be calculated. Figure 4 shows the highest probability zone by location (top left) along with the probabilities of each location being assigned to each of the three zones. The predictions appear to be fairly spatially contiguous zones (as desired) with the exception of middle of the field in areas between zone 1 and 2.

While zone prediction map in Figure 4 produces zones that are spatially smoothed by the spatial basis function, the covariates produce some non-contiguous regions that would be difficult to implement using variable rate irrigation given the current state of precision agriculture. Rather, for VRI implementation purposes, we desire to further smooth the predicted zones into purely contiguous zones. For purposes of VRI implementation, we use the spatial clustering algorithm of [19] based on finite differences to achieve more contiguous zones. Specifically, we spatially cluster using ward dissimilarity [see 19, for details] based on the expected zone (3.1) E ( Y ( s ) ) = 1 × [ Φ ( 0 − ( f ˆ L ( X ( s ) ) + w ˆ ( s ) ) ) ] + 2 × [ Φ ( c ˆ 2 − ( f ˆ L ( X ( s ) ) + w ˆ ( s ) ) ) (3.2) − Φ ( 0 − ( f ˆ L ( X ( s ) ) + w ˆ ( s ) ) ) ] + 3 × [ 1 − Φ ( c 2 − ( f ˆ L ( X ( s ) ) + w ˆ ( s ) ) ) ] where Φ ( · ) is the standard normal cumulative density function. The result from this clustering process is shown in Figure 5 and compared with the original predictions. While there are still some discrepancies between the clustered zones and the original zones, they are relatively similar and the clustering provides contiguous zones that would be viable to be implemented with VRI. We note that clustering the model results (using any spatial clustering method) is not required but is simply beneficial for our application and the implementation of the estimated zones for the growing season.

While prediction was the primary goal of the analysis, it may also be of interest to examine the effect of the covariates on zone classification. However, one of the drawbacks of using neural network models as we have here is that the parameters have no intuitive interpretation. Indeed, interpretability of neural networks, generally, is an open area of research [29]. One possible solution is to use partial dependence plots (PDPs) and feature importance to intuitively understand the effects of covariates on the response. First, partial dependence plots are used to show the marginal effect of a covariate on the predicted outcome. Mathematically, the partial dependence plot for a covariate, say X p , is calculated as (3.3) PDP ( X p ) = 1 n ∑ i = 1 n ( f ˆ L ( X p , X − p ( s i ) ) + w ˆ ( s i ) ) where X − p ( s i ) is the vector of covariates with X p ( s i ) removed and f ˆ L ( · ) and w ˆ ( · ) is our fitted spatial neural network model. The above is calculated for a grid of X p in the domain of X p ( s ) to produce a curve representing the marginal effect of X p ( s ). Intuitively, the partial dependence measure is calculated by replacing the variable of interest ( X p ( s )) with a fixed singleton value for all observations in the data and averaging the associated prediction across observations. Because we adopted the Bayesian approach, the uncertainty for these partial dependence plots are accounted for as well via Monte Carlo sampling.

Figure 6 shows the partial dependence plots for some of the covariates that were used in the analysis (covariates with no effect were omitted for brevity) where the black line represents the posterior mean of the marginal effect for the covariate and the shaded region is the 95% credible interval for the marginal effect. Elevation, yield, and the NDVI for 2018 and 2019 seemed to heavily influence the probability of belonging to an irrigation zone. For example, as elevation increases, the probability of belonging to a lower-water zone increases (Zone 1 corresponds to the zone with the smallest VWC which, in turn, would get the most irrigation to compensate). The other covariates appear to have a minimal effect on the predicted zone.

Another common measure for covariate influence on machine learning models is feature importance. For this analysis, the permutation feature importance was used as the feature importance measure. The feature importance was calculated by permuting values of the feature (covariate) and calculating the ARI for the model with the permuted feature. The purpose of permuting the covariate is to break the relationship between the feature and the response. Under permutation, the prediction accuracy should decrease for features that are highly predictive of the response. After calculating the ARI for the model with the permuted feature, the feature importance for the p th covariate ( X p ( s )) is calculated as: (3.4) F I p = A R I orig − A R I perm where A R I orig is the ARI for the model before permutation and A R I perm is the ARI for the model after feature X p ( s ) had been permuted. Intuitively, higher F I p equates to more important features in determining the zone delineation.

The feature importance plot in Figure 7 shows the feature importance for the 10 covariates that were included in the model. In correspondence with the partial dependence plots, the results from the feature importance plot also indicates that elevation, yield, and the NDVI index from 2018 and 2019 are important covariates in determining the irrigation zone. Hence, for other agricultural fields where irrigation zones are desired, collecting the highest importance features will provide the most effective zones for variable rate irrigation.