Statistical methods for identifying anisotropy in the Spodoptera frugiperda spatial distribution

Corn is a very important agricultural product, however, some pests may cause damage to the corn productivity such as Spodoptera frugiperda, which prevents the plant from growing in a regular manner. Since the indiscriminate use of the pesticide may cause an increasing resistance of the insect besides an environmental damage, it is important to estimate the areas and the dominant directions where the insect may propagate. The main aim of this work was to study the spreading of the fall armyworm in a commercial agricultural area in the South of Brazil. For this, we considered a set including the location of each corn plant attacked by the insect. In particular, we assumed that the spatial locations given by the geographic coordinates constitute a spatial point pattern following a stationary Poisson point process. In order to detect the presence of possible dominant directions in the distribution of the fall armyworm infestation we studied the anisotropic features of the data by using some second-order spatial point-pattern analysis techniques such as the K directional test, the wavelet-based test, and the quadrat counting test. All the results showed that spatial distribution of fall armyworm may follow a clustered Poisson point process with the presence of an evident anisotropy mainly due to the shape and the distance between corn plants of the experimental area. These preliminary results could be used for reducing and optimizing the use of pesticides with a consequent decrease of the environmental impact. Additional keywords: cluster processes; corn pests; spatial point pattern; directional K function; wavelet based test. Abbreviations used: CSR (Completely Spatial Random), GPS (Global Positioning System), UTM (Universal Transverse Mercator). Authors’ contributions: Conceptualized the paper, statistical analysis of data, final revision and discussion: DTN, ON, MAUO & FDB. Reviewing the literature and editing the working versions of the manuscript: DTN & ON. Citation: Nava, D. T.; Nicolis, O.; Uribe-Opazo, M. A.; De Bastiani, F. (2018). Statistical methods for identifying anisotropy in the Spodoptera frugiperda spatial distribution. Spanish Journal of Agricultural Research, Volume 16, Issue 1, e1003. https://doi. org/10.5424/sjar/2018161-11916 Received: 20 Jun 2017. Accepted: 13 Mar 2018. Copyright © 2018 INIA. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 International (CC-by 4.0) License. Funding: Fundação Araucária, Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), Post-Graduate Program in Agricultural Engineering (PGEAGRI-Unioeste) and Universidade Tecnológica Federal do Paraná (UTFPR). Competing interests: The authors have declared that no competing interests exist. Correspondence should be addressed to Orietta Nicolis: orietta.nicolis@uv.cl


Introduction
The fall armyworm is one of the pests that cause the greatest damage to the corn (Zea mays L.) crop, leading to a loss of about 34% of the grain production (Cruz et al., 1999). The insects feeds in the whorl of corn from its emergence to the formation of spikes, and if not controlled causes the death of the plant. In the world, corn is one of the most valuable cereal grain crop, it is used for both human and animal consumption. Corn is also used to make numerous non-food products. It can be cultivated in a large part of the world, and Brazil is a big world producer of the grain. Given its usefulness, it is extremely important to take into account its productive potential, resistance to diseases, types of soil besides the climatic conditions concerning to optimize the yield of grain crops.
In general, the pest control is made through the chemical pesticides, which besides polluting the environment also cause the increasing of the resistance of the insect. For this reason control techniques for the insect are researched. For example Sena Jr. et al. (2003) developed an algorithm for identifying damaged corn plants by the fall armyworm using digital colour images. Yu et al. (2003) studied detoxification capability in larval fat bodies and midguts, and adult abdomens of insecticide-resistant and -susceptible fall armyworms. They concluded that the insecticide resistance observed was due to multiple resistance mechanisms. Niu et al. (2013) studied the resistance of the fall armyworm to certain types of protein with regards to propose a hybrid resistant to the caterpillar. In Florida and Lousiana Yang et al. (2013) determined the susceptibility of S. frugiperda to a pyramided Bt corn hybrid containing Agrisure ® Viptera™ 3111 corn.
Statistical methods based on spatial point processes provide useful tools for points randomly distributed on the plane or in space in many fields of knowledge, such as ecology, biology, epidemiology, seismology, archaeology, astronomy, and geography. However, most of the analyses are made assuming isotropy of the process without an adequate confirmation (Guan et al., 2006). This fact is often due to a lack of studies that consider the investigation of the isotropy in this type of data. In ecology, seismology, and biology, Guan et al. (2006), Mateu & Nicolis (2015) and Rajala et al. (2016) presented a few studies for anisotropy. The first one, proposed a formal test for isotropy based on the asymptotic joint normality of the sample second-order intensity function. Mateu & Nicolis (2015) proposed a wavelet-based test, and Rajala et al. (2016) developed a two-stage non-parametric method for quantifying geometric anisotropy for patterns that are compressed or stretched.
The use of spatial point pattern analysis in agriculture was considered by Spósito et al. (2007), who studied the disease of the black spot in an orange crop in Brazil. They evaluated the point pattern of the infected trees by three different methods, the binomial dispersion index, the Ripley K function and Monte Carlo test, and concluded that the disease shows a clustered behavior in the region under study . Ribeiro Jr. et al. (2009) used point processes to detect and model specific patterns in the occurrence of the thrips disease in onion. They used the Mantel test to detect the point pattern and later described the pattern by a mixed spatial Poisson process with a geostatistical random component. They also obtained maps of the prediction of the levels of susceptibility to infestation in the study area.
In this work we study the spreading of the fall armyworm along an agricultural experimental area at Paraná State, Brazil, comparing methods of anisotropy by means of the spatial point patterns. The main purpose is to examine its spatial distribution along the study area, investigating possible dominant direction of infestation of the fall armyworm, and if the pattern presents characteristics of clustering, regular or random spreading. This knowledge provides the farmer a more adequate management of his crop by considering the characteristics of the dispersion of the insect along the area in consortium with the techniques of precision agriculture. The lack of literature that incorporates spatial point processes and agriculture motivated this work.

The data set
The experiment was carried out in an agricultural commercial area of 27 ha located in the city of Cascavel, Paraná state, Brazil (approx. 24.95 o S, 53.57 o W and 650 m asl). Local soil is classified as Rhodic Hapludor (Soil Survey Staff, 2014). The climate of the region is classified as mesothermic and super humid temperate with average annual temperature around 21 o C. A representation of the location of the experimental area is presented in Fig. 1a.
For the data sampling, we considered as a possible event each plant of corn attacked with the fall armyworm in the experimental area during the crop year 2015/2016. For every observed infected plant we registered the geographical coordinates as an event of interest, which were obtained using a GEOEXPLORE 3 GPS positioning system receiver in a Universal Transverse Mercator (UTM) coordinate system.
The study area ( Fig. 1) included 9360 corn plants of which 1303 were infected by the S. frugiperda. This study area was cultivated with the hybrid Pioneer 30F53YH variety, developed by the company DuPont Pioneer, under direct seeding system, with a population of 82,222 plants/ha. Along the season we applied 1095 kg of fertilizers (8% N, 20% P and 18% K). The planting of the seeds began on August 28, 2015 and ended on August 31, 2015 and was performed with line spacing of 0.45 m and between plants of approximately 0.33 m. The sampling was made on October 15, 2015, when the corn crop was in vegetative growth stage.

Statistical analysis of anisotropy
A spatial point pattern is a stochastic mechanism that generates a countable set of events x i in the plane (Diggle, 2003). Events can be scattered in the region of the plane in three forms of relationship, namely completely spatial random (CSR) as represented in Fig.  2a, regular in Fig. 2b, and clustered in Fig. 2c. In spatial point analysis, the main objective is to identify possible properties and forms of relationship of the phenomenon under study taking into account its spatial sites, through geographical coordinates.
A spatial point process is considered to be stationary under translations if the distribution is unchanged when the origin of the index set is translated (Ripley, 1981;Diggle, 2003;Illian et al., 2008), and is said to be isotropic if its distribution is invariant under rotations about the origin. Otherwise the process can be said to be anisotropic. The assumption of stationarity simplifies drastically the statistics of point patterns (Baddeley et al., 2006) because the chance of observing some point configuration at a specific location is independent of the location. Usually, the assumption of isotropy is made by simpler interpretation and ease of analysis (Guan et al., 2006). In other words, stationarity implies that it is possible to estimate some properties of the process from a single realization on the region A, because these properties are the same in different, but geometrically similar, sub regions of A. Isotropy means that they are invariant with respect to directions (Mateu & Nicolis, 2015). In the following, we assume that the process N is stationary, but not necessarily isotropic. The simplest example of a stationary and isotropic point process is the homogeneous Poisson point process.
Figs. 2(a), 2(b) and 2(c) show simulated CSR, regular and cluster spatial point patterns, respectively. The difference of them is the distance between the points. For the regular pattern the distances between neighbor points appear to be constant, while for the cluster pattern there are groups of points over the region. In the CSR pattern the distance between the points is random.
Consider N, a two-dimensional spatial point process, where each event is given by . Let A be the region where the points are taken, denote as |A| the area of this region, and as N(A) the random number of points in A. The intensity function may be described by the first order property, given as where E[.] is the expected value; N(.) is the number of events in the plan region; dx is a small region containing the event x; |dx| is the area of the region dx. For a stationary process λ(x) = λ , where λ is the average number of events per unit area (Daley & Vere-Jones, 2002).
The estimation of the first-order property, or intensity of the point pattern, is usually performed through kernel functions (Baddeley et al., 2006(Baddeley et al., , 2015, where k(.) is the kernel function with ʃ k(u)du = 1, and h>0 is the smoothing parameter. The intensity estimate at the point x depends only on the spatial relationship between x and the elements of the sample, x i (i = 1, …, N(A)), quantified by the metric embedded in the kernel function. In order to estimate λ(x) by using the kernel function as described in Eq.
[2], we used the kernel2d function of the package splancs (Rowlingson & Diggle, 2016) represented by Figs. 2d-2f. The kernel intensity function estimation using quartic smoothing parameter, h=0.2, was performed for each simulated pattern and it is shown in Figs. 2d, 2e and 2f corresponding to the estimation of the first-order property for each type of simulated point pattern. It is worth to note that the intensity seems almost constant for the simulated regular pattern as in Fig. 2e, and for the simulated cluster pattern there are regions with different values of intensity as in Fig. 2f.
The second-order properties of a point process can be observed using the covariance between events of two areas. The second-order intensity provides a way of describing the relationship between pairs of points and is defined as where dx is a small region containing the event x; dy is a small region containing the event y and |dx| and |dy| [3]  (Illian et al., 2008); for a CSR process we have K(r) = πr 2 . For a cluster process the points have more neighbors than expected under CSR, and therefore the estimates for K(r) will be greater than πr 2 . On the other hand, for a regular pattern, the estimates of K(r) will be less than πr 2 (Diggle, 2003). Fig. 3 represents the K function (Eq. [4]) and the empirical envelope of 95% for each pattern simulated in Figs. 2a, 2b and 2c, respectively. A directional counterpart of Ripley K function was proposed by Rajala et al. (2016) as a non-parametric methodology to study anisotropy in spatial point patterns. Instead of counting all data points falling inside a circle of radius r centered at a data point, we could replace the circle by another geometrical shape (Baddeley et al., 2015). For a unit vector u, the conical directional K-function is defined as where λ is the intensity, and C u (r,θ) denotes a double cone in the direction u with an slant height of length r and an apex angle of size 2θ centered in 0 (Redenbach et al., 2009).
The function is obtained by counting how many pairs of events in the pattern have both their vector of difference angle less than θ radians from direction u and difference vector length less than range r. The procedure consists in comparing each of the directions that is defined by the researcher with the function under CSR. In order to do this, it is necessary to assume that are the areas of regions dx and dy, respectively. For a stationary process, λ 2 (x,y) = λ 2 (x-y), and for a stationary, isotropic process, λ 2 (x-y) = λ 2 (r) where r =||x-y|| is the Euclidean distance (Daley & Vere-Jones, 2002). Ripley (1976Ripley ( , 1977 presented an alternative characterization of the second-order properties (Eq. This function has shown its importance as a spatial dependency summary since it allows quantifying spatial dependence between different regions of the process, and is also known as reduced second-order analysis. It has been widely used in two-dimensional point processes because it enables the detection of the point pattern in different distance scales simultaneously and the observed point pattern can be compared to known stochastic process models for different point configurations (Diggle, 2003).
The shape of K(r) relative to that of the CSR provides valuable information on the point process An alternative method to study the isotropic behavior of the data set is the histogram of the k-th nearest neighbor angle. For each event in the point pattern, we calculate the distance to the k-th nearest neighbor and the azimuth angle α between this pair of events where α ϵ [-π, π]. The isotropic behavior of the data set is done by comparison with the theoretical value 1/2π where peaks above this theoretical value indicate a possible dominant angle direction of anisotropy. The histograms can be used to investigate in which nearest neighbor the process reaches isotropy, that is, when there are no more peaks in the histogram. We used the knnangle function of the Kdirectional package (Rajala, 2017) to build the histograms.
To confirm the form of relationship of the point pattern we calculated the quadrat counting test (Illian et al., 2008). It is a test of CSR that uses the χ 2 statistic based on quadrat counts where the study area A is divided into sub regions called quadrats (A 1 , A 2 , …, A q ) of equal area. The test counts the number of points that fall in each quadrat n j = N(A ∩ A j ) for j=1, ..., q. Under the null hypothesis of CSR, the n j are independent and identically distributed Poisson random variables. The test statistic is given by where s 2 is the sample variance. If the distribution is Poisson, the variance is exactly the mean under CSR anisotropy is geometrical, i.e. a result of compression or stretching and modeled by a linear transformation comprised of scaling and rotation (Rajala et al., 2016).
Mateu & Nicolis (2015) developed an anisotropy wavelet-based test by using the empirical logarithm of the directional scalogram, assuming that under the null hypothesis of isotropy a random point pattern is expected having the same value of the directional scalogram for all possible directions. Let f(x), x R 2 , be a function. The continuous directional wavelet transform for scale a and orientation θ is given by where the over line denotes complex conjugate. There are many different directional wavelets Ѱ a,b (x, θ) proposed in literature, see for example Neupauer & Powell (2005). Kumar (1995) proposed the function with regard to identify the behavior of the process in different directions. η(a, θ )characterizes the distribution of the energy at different scales and directions. The directional scalogram is obtained from the variance of wavelet coefficients for each scale and each direction, given by Mateu & Nicolis (2015) applied the fully anisotropic Morlet wavelet transforms given in Neupauer & Powell (2005)    For the data analysis we used the software R version 3.3.3 (R Team, 2017) and the R packages spatstat version 1.47-0 (Baddeley & Turner, 2005) and Kdirectional version 1.01 (Rajala, 2017).

Results
The Fig. 1a represents the location of the corn plants attacked by the S. frugiperda in the study area. The scattering of the points along the region does not seem to present any evident dominant direction for the attack of the fall armyworm.
By applying the kernel method of Eq.
[2] to estimate the intensity function, we observed some directional clustering, corresponding approximately to the angles of 135 and 45 degrees (Fig. 4). The method allowed to see how the concentration changed in the studied area but it was not very useful for detecting anisotropic patterns since it is based on isotropic kernel functions.
A first study of anisotropy is performed by the histogram of the nearest neighbor angle, which shows a descriptive behavior of the isotropy of the point pattern. Fig. 5 represents the histograms for the first and third nearest neighbor angle, respectively. The histogram obtained for the first nearest neighbor angle (Fig. 5a) shows two main directions given by the highest peaks. The first direction corresponds approximately to the angles 135 and 315 degrees, and the second highest peaks due to the angles 45 and 217 degrees. It is worth to note that these directions correspond to the directions used by the planter machine. When considering the histogram of the third nearest neighbor angle (Fig. 5b), there were only two main dominant directions, 135 and 315 degrees, however there was a reduction in the density value. We also noted that the process tended to be below the theoretical line to all the others degrees.
With regard to correct the anisotropy observed in the histogram analysis, we rotated the data set by 45 degrees clockwise. The rotated point pattern is represented by Fig. 6b. And, in order to eliminate the border effect of the design of the experimental area, we selected a number of points in the center of the area which is represented by Fig. 6c.
In the following, we calculated the directional K given by Eq. [5], and the wavelet test of Eq.
[7] to the three data sets for comparing the isotropy/anisotropy features. The second-order property, as defined in Eq.
[3], is estimated by the directional K function.
We used a cone shape with angles 2θ = 0, 45, 90 and 135 degrees with respect to the x-axis, since these directions are normally used in the literature to study anisotropy. Fig. 7 represents the directional K function for each pattern and each of the considered directions with an empirical envelope of 95% for the point pattern under the null hypothesis of CSR. It is evident for the original data strong clustering for each direction, as visualized by the directional K function in Fig. 7a. We observed that the most dominant direction of anisotropy   Fig. 7b, we observed a reduction of anisotropy and a weak clustering. In this case the most dominant direction of anisotropy was 90 degrees, followed by the zero degree direction. By taking the rotated center region of the study area, represented by Fig. 7c, we see a slight clustering and a very weak anisotropy. However, all the considered directional K functions were outside of the envelope, so we could assume that the process was not completely random. Fig. 8 represents the directional scalograms for our data sets by applying the directional Morlet wavelet transform for each degree (from 1 to 180) and scale (from 1 to 30). For the original point pattern (Fig. 8a), we found an evident anisotropy in the direction of approx. 135 degrees in the highest level of resolutions and spreading on a larger angle range (between 80 and 140 degrees) in the lowest levels (close to a=30). For the rotated point pattern (Fig. 8b), the directional scalogram showed a weak anisotropy with a dominant direction at the 80 degrees approximately. Finally, the scalogram evaluated on the central point pattern data (Fig. 8c) did not show any dominant direction.
We then performed the wavelet-based test by using the empirical logarithm of directional scalogram as described in Eq. [6]. The blue line in Fig. 9 represents the values of T(θ i ), i=1, … , 180 for the S. frugiperda data set. In order to test if the process was CSR we simulated 1000 isotropic Poisson point patterns and we evaluated the 2.5% and 97.5% of the total empirical distribution (red lines of Fig. 9). For the original data set (Fig. 9a), the T statistic is out of the 95% empirical confidence interval for all the angles showing a peak around the angle 135. In this case the test clearly rejects the hypothesis that the process was Kd 8 isotropic. When we performed rotation of the data set by 45 degrees, the resulted T statistic (Fig. 9b) did not present significant dominant directions except for the degrees in the range 80-85 which were slightly out from the empirical confidence interval. By observing the T statistic evaluated on the central rotated pattern (Fig. 9c), there were not any dominant directions and the considered point pattern seemed completely random, confirming the results obtained by the K-directional analysis. We think that the anisotropy of the original data set may be due to the shape of the experimental area: it seems that the insects tend to propagate on the closest plants. By correcting the anisotropy, the wavelet analyses showed a second dominant direction close to 80 degrees that was not evident by the K-directional test. Since the centered data set seemed CSR we thought that the anisotropy of Fig. 9b was localized close to the border of the experimental areas. The main advantage of the wavelet test respect other statistical methods is that it is able to detect local anomalies of the data at different levels of resolutions.
In order to analyze the clustering feature of the data we applied the quadrat counting test, as defined by Eq. [8]. Table 1 shows the obtained results. We should reject CSR hypothesis and conclude that our spatial point pattern was clustered. The results showed that for the vegetative growth state of the corn, in which the data sampling was performed, the fall armyworm presented clustering behavior.

Discussion
The knowledge of the population spreading of S. frugiperda, the pest that causes a great deal of damage to the corn crop, is extremely important for the development of efficient techniques of both monitoring and crop management. In this research we applied spatial point processes techniques to study the dispersion of the fall armyworm in an agricultural commercial area.
From the estimation of the intensity function of the fall armyworm data (Fig. 4) we observe a major concentration of this insect in the zone close to the left border of the experimental area. However other small regions seemed characterized by clusters. The clusterization of the armyworm point pattern is also confirmed by the quadrat counting test (Table 1).
Our results also showed a significant anisotropy when considering the original data set. Both the K-directional and the wavelet based methods detected anisotropy along the directions 45 and 135 degrees, probably due to the design given by the planter machine and to the shape of the area. With the view to correct this anisotropy, we rotated the pattern by 45 degrees clockwise. The resulting rotated point pattern did not present any dominant directions by using the K-directional method. Instead, the wavelet method detected a slight anisotropy around the angle 80. The different result obtained by the wavelet method respect the K-directional is probably due to the ability of the   χ wavelet transform to locally analyze the point pattern at each angle and resolution scale.
In order to avoid any border effect, we selected a central subset of the rotated pattern and we estimated the K-directional and the directional scalogram on this last data set. In this case both methods did not detect any dominant directions and showed the completely randomness of the process. The difference between anisotropy studies techniques presented here consists of the number of resolution scales, once the wavelet test was able to identify dominant directions at scales finer than a K directional function. However, when comparing both results, we reached the same conclusion. To our knowledge, there are no results in the literature about the study of anisotropy for spatial point patterns in agriculture, for this reason it is not possible to make comparisons with other results. However, there are some works which analyze the clustering property of the S. frugiperda (e.g., Farias et al., 2008 andRios et al., 2014). We think that our results on clustering are in general coherent with those of Farias et al. (2008) and Rios et al. (2014), although in our study we did not take into account the size of the fall armyworm. In particular, Farias et al. (2008) found that the spreading of the S. frugiperda for small larvae present aggregated behavior, and it becomes increasingly random with larval growth. Hernández-Mendoza et al. (2008) also concluded that the spatial distribution of the fall armyworm was strongly associated with the corn phenological stages. Since the experimental area considered in this work constitutes a small part of the full-size corn crop field, we think that more consistent results on anisotropy could be reached if using a larger area. We also think that better results on the S. frugiperda spatial distribution could be obtained by including the size of this insect. Since our data base did not contain this information we will consider these issues in future works together with the use of highresolution satellite or drones images for having a larger sampling area. Finally, we are planning to improve the estimation of the intensity function by using estimation techniques based on directional kernel and directional wavelet functions.