Spatial pattern detection modeling of onion thrips (Thrips tabaci) on onion fields
Paulo Justiniano Ribeiro Jr1; Clarice Garcia. Borges Demétrio2*; Denise Nunes Viola 3; Bryan F. Manly2, Odair Aparecido Fernandes4

1 UFPR – Lab. de Estatística e Geoinformação - C.P. 19.081 - 81531-990, Curitiba, PR - Brasil.
2 USP/ESALQ - Depto de Ciências Exatas, Caixa Postal 9 - 13418-900, Piracicaba, SP – Brasil.
3 UFBA/DES - Depto de Estatística, Av. Adhemar de Barros, s/n - Campus de Ondina - 40170-110 - Salvador, BA – Brasil.
4 UNESP/FCAV - Depto de Fitossanidade, Via de Acesso Prof. Paulo Donato Castellane s/n, 14884-900, Jaboticabal, SP- Brasil
* Corresponding author <clarice@esalq.usp.br>

Spatial pattern detection modeling of onion thrips (Thrips tabaci) on onion fields

ABSTRACT: Onion (Allium caepa) is one of the most cultivated and consumed vegetables in Brazil and its importance is due to the large workforce involved. One of the main pest that affect this crop is the onion thrips (Thrips tabaci), but the spatial distribution of the insect, although important, has not been considered in crop management recommendations, experiment planning or sampling plans. Our purpose here is to consider statistical tools to detect and model spatial patterns in the occurrence of onion thrips. In order to characterize the spatial distribution pattern of the onion thrips a survey was carried out to record the number of insects in each development phase on onion plant leaves, on different dates and sample locations, in four rural properties with neighboring farms with different infestation levels and planting methods. The Mantel randomization test proved to be a useful tool to test for spatial correlation which when detected was described by a mixed spatial Poisson model with a geostatistical random component and parameters allowing for a characterization of the spatial pattern as well as the production of prediction maps of susceptibility to levels of infestation throughout the area.
Key words: spatial statistics, randomization tests, geostatistics, Poisson distribution

Detecção de padrões espaciais na ocorrência do tripes
do prateamento Thrips tabaci na cultura da cebola

RESUMO: A cebola é uma das hortaliças mais cultivadas e consumidas no Brasil e sua importância social se deve à grande demanda por mão-de-obra. Uma das principais pragas que afeta essa cultura é o tripes do prateamento (Thrips tabaci) e sua distribuição espacial, embora importante, não tem sido considerada nas recomendações de manejo da cultura, planejamento de experimentos ou estudos amostrais. O objetivo desse artigo foi considerar métodos estatísticos para detectar e modelar padrões espaciais na ocorrência do tripes do prateamento da cebola. Para caracterizar o padrão espacial da dispersão do tripes do prateamento da cebola foi feito um levantamento anotando-se o número de insetos por fase de desenvolvimento em folhas de plantas de cebola, em diferentes datas e pontos amostrais dentro de quatro propriedades com fazendas vizinhas apresentando diferentes níveis de infestação e métodos de plantio. O teste de aleatorização de Mantel mostrou-se útil para testar a presença de padrão espacial, que quando detectado foi descrito por um modelo de Poisson misto espacial com componente aleatório geoestatístico com parâmetros que possibilititam a caracterização do padrão espacial bem como a obtenção de mapas de predição dos níveis de susceptibilidade à infestação na área.
Palavras-chave: estatística espacial, testes de aleatorização, geoestatística, distribuição de Poisson 

INTRODUCTION
Onion (Allium caepa) is one of the most cultivated and consumed vegetables in Brazil. The social importance of the crop is due to the large workforce involved. It is estimated that 70% of the production is small scale, because it is typically grown on small and medium sized properties. It is an annual plant for bulb production, biannual for seed production, and propagated by direct sowing, bulbs or seedlings planted in beds and transplanted to the field.
One of the main pest that affects onion crops is the onion thrips (Thrips tabaci), which in high infestation levels can damage the harvest (Workman & Martin, 2002), with reduction in the production of up to 80% during hot and dry periods (Sato, 1989).
	The insect is typically found at the base of leaves. It feeds from the sap and the leaves parenchyma causing gray spots, which gradually change to silver as a result of the external tissue damage of the leaves. Massive attacks on the aerial part of the plant cause loss in bulb production, which reduces the size and quality, damaging the commercial value and creating obstacles to exports.
When an attack is very intense, the leaves get yellowish, dry and with wrenched tips, causing the wilting and the death of the plant (Sato, 1989), and also allowing for the entrance of water to the bulb, which gets rotten. The insect is also considered a vector of a phytopathological agent with the capacity to transmit a virus to the plant.
	The insect development occurs in the four phases of egg, nymph, pupa and adult, with the nymph and adult stages damaging the production, because the pupa phase is restricted to the soil. The nymph has low mobility, while the adult, although winged, has restricted movement. The development cycle varies typically from 14 to 30 days, changing to 10 and 11 days when the temperature is over 30oC.
	The spatial distribution of thrips in commercial fields is important for the efficient application of insecticides. However, this has not been considered in crop management recommendations, experiment planning and sampling plans. Considering the low mobility of nymphs and adults it is reasonable to assume that the wind is the main dispersion factor for the thrips that potentially determines the spatial pattern.
	A spatial pattern can be classified as random, aggregate or uniform. The random pattern occurs when there is a constant and independent probability of infestation for all the plants, while the aggregate pattern is associated with low insect mobility. The uniform pattern rarely occurs naturally, but can be induced, for instance, by alternated planting of resistant and susceptible plants. In order to study whether infant leukemia cases tend to be close in space and time, Mantel (1967) proposed a randomization test, based on matrices of time and space distances between observations. This test can be used to test for spatial correlation in an insect distribution, but its usage has not being considered in practical applications, and in particular, in studies of the spatial distribution of the onion thrips.
	It is common in insect distribution studies, to find the use of indices based on the relationship between the variance and the mean, such as the David & More index, the Taylor power law, and the aggregate indices of Lloyd and Iwao, among others (Ruiz et al., 2003). However, these indices ignore the spatial location of the samples, have limited capacity to describe spatial patterns, and strongly depend on the size of the sample unit.
	Geostatistical methods (Isaaks & Srisvastava, 1989; Goovaerts, 1997) have been used to describe insect spatial patterns as, for instance, in Grego et. al. (2006). Such methods were originally developed for continuous response variables, with several computational implementations available for data analysis. The insect counts are discrete and typically distributed in clusters, with many zero counts.  Therefore, the data cannot have a covariance structure of the type assumed by traditional methods of geostatistical analysis, with a stationary spatial covariance structure in the study area (Ruiz, 2002).  For this reason it is appropriate to use models that incorporate explicitly a data generating mechanism such as the Poisson distribution, combined with structures that describe the spatial pattern of the counts. These kinds of models have been proposed in the statistical literature (Diggle et al., 1998) but have had few practical applications.
	This paper describes a study of the spatial distribution of onion thrips with data from surveys of four different properties with different infestation levels and planting methods. We aimed to detect spatial patterns in the occurrence of onion thrips at different production fields and propose an statistical model for such patterns. We adopt the Mantel randomization test (Manly, 2006) to decide for the presence of spatial autocorrelation which when detected was modeled by a mixed spatial Poisson model with a random term given by a geostatistical component. This model allows the characterization of the spatial pattern as well as the production of maps of levels of susceptibility to infestation in different areas.

MATERIAL AND METHODS
Data description
	This work is motivated by a set of data originated from a study involving sampling onion thrips in onion crop in four different farms, located in the municipality of São José do Rio Pardo, São Paulo State, Brazil (21°36’S, 43°53’W; altitude 705m), from June to September, 1996. The aim is to study the spatial and temporal distribution of thrips. The four chosen properties used the onion hybrid Granex 33 and the seedling planting method. The trial areas were chosen with neighbors who adopted different kinds of planting and had different infestation levels.
	Details referring to the kind of planting in the neighborhood and collection dates and numbers of samples collected in the different farms are shown in the Table 1. The São Paulo farm is located at a high elevation of the region and the nearest neighboring onion crop is situated over one kilometer away. The neighborhood of Estância Bela Vista had already had some crops attacked by onion thrips pest.

Table 1 - Characteristics about the data precedence, types of neighbors, sample times and number of samples. 
Farm	Neighborhood	Sample times	Number of samples	
Fazenda 
São Paulo	Isolated from 
other plantings	10/7, 24/7, 31/7, 7/8,
14/8, 21/8, 28/8, 04/9	100, 100, 100, 98 100, 100, 100, 100	
Estância 
Bela Vista	Bulbs	11/7, 1/8, 8/8,
14/8, 9/9	100, 100, 84,
99, 99	
Sítio 
Rosário	Seedlings	21/6, 29/6, 7/7, 14/7, 21/7,
28/7, 4/8, 11/8, 18/8, 25/8, 3/9	50, 50, 48, 50, 50,
50, 50, 50, 50, 50, 50	
Sítio 
Novo II	Seedlings	4/6, 19/6, 27/6, 28/6,
4/7, 11/7, 24/7, 31/7, 7/8	100, 100, 100, 100,
100, 100, 100, 100, 100	

	The sampling unit was a 1m radius circle with a center stake. One plant was then randomly selected from within the circle. The position of the stakes in the four studied farms, in general with a 10x10m grid, but with some variations at Fazenda São Paulo is shown in Figure 1. The measured variables were the stake location on the coordinate axes, the number of nymphs, the number of adult insects and the number of leaves per plant. The number of samples and sampling times varied from farm to farm as shown in Table 1. The response variables are discrete because of result of counting. In some cases, the counts are multiples of 5 or 10 and some values over 100 were truncated to 100.

Figure 1- Localization of the stakes in each farm.
	Figure 2 shows box-plots for the average number of insects per leaf, at the four farms, for the various sample times. There is great variability in the counts and also some outliers, not all of them being influential on the model fitting. At the São Paulo farm the average number of insects and the variability increased with time, while at the other farms, the average increased and then decreased. In all cases the observations above the median are more variable showing positive asymmetry, with some extreme values.

Figure 2 - Box-plots for the average number of insects per leaf.
	At the São Paulo farm the lowest average number of insects per leaf and also the lowest variance wer found at 31/07, with one insect per leaf as the maximum value. In contrast on 04/09 this farm had a much larger average number of insects per leaf and much greater variability. The percentage of infested plants ranged from 35% to 100%. For the Estância Bela Vista, the lowest average number of insects per leaf occurred on 11/07 and 14/08 with 89% to 100% plants infested. The Rosário farm had only 50 plants sampled and the highest average number of insects per leaf on 11/08 and 18/08. Sítio Novo II had the least average for the number of insects per leaf with low variability except for one outlier count of 30.

Mantel’s test for the detection of spatial pattern
The non existence of spatial pattern in the dispersion of insects may be considered a randomization hypothesis, and the existence of a spatial pattern can be tested through the randomization of the order of the observed values (Manly, 2006).
Randomization tests are based on the fact that, if the null hypothesis is true, then all of the possible orders of the data have the same chance of occurrence. Therefore, the value eo of a statistic E is calculated for a set of observations, and then a large number of randomizations are made. For spatial data these randomizations are made by randomly reordering the data. For each randomization a value ea is calculated and the set of the ea values generate an approximation of the randomization distribution of E. Just as for classic statistical tests, the decision is guided by a p-value, which in the case of randomized tests is given by the proportion of the ea values that are larger than or equal to eo, for a one-sided test. For instance, if p < 0.05, it’s concluded that there is evidence that the null hypothesis is not true (Manly, 2006).
Randomized tests have some advantage in comparison to classic statistical tests. For example, the statistics are usually easy to calculate, relatively to the classic statistical tests. They are based on non standard statistics and they do not need previous information about the population from which the samples were taken. Also, they can be applied with non-random samples which can consist only of the data that need to be analyzed (Manly, 2006). However, the randomization tests are easier to justify when the analyzed samples are random or the experimental design suggests a randomization test. 
Usually, when considering spatial data, it is desired to test the null hypothesis of a random spatial pattern versus the alternative of a non-random spatial pattern. A test for this hypothesis was proposed by Mantel (1967). The test is implemented as follows. Let a variable be observed in n locations. Two symmetric matrices A and B are obtained, each with n x n dimensions. The elements represent distances between the observations. These matrices can be denoted as

The matrix A is the Euclidian distances between the stakes with locations given by  and  i.e., with elements of the form  and B is the matrix with elements , where Z is the mean of the number of insects per leaf. The test statistic is given by the Pearson correlation coefficient between the correspondent elements of A and B, i.e.,
                                                            (1)
which produces the  value when calculated for the observed values. For the randomization test the rows and columns of one of the matrices are permutated a large number (N) of times, and the values   are obtained, for k = 1, 2,…, N. The proportion p of values  is then compared with a pre-fixed significance level a (for example, 0.05) and the null hypothesis is rejected if p < a (Manly, 2006). 
As the matrices A and B are symmetric, the correlation amongst all the elements outside the main diagonal is the same as the correlation of the  elements in the upper or lower triangular part of the matrix. Note that the only term of (1) that is altered by changing the order of the elements in one of the two matrices is the sum of products .
Other possible metrics used for the calculation of the distances are Euclidian with standardized data, Euclidian squared, Euclidian squared with standardized data, proportional distance and sample difference. The alternative is given by Snäll et. al. (2003) who built a randomized test using flexible forms for the relation between the distance measurements, given by the structure of additive generalized models.
When the Mantel test rejects the null hypothesis there may be interest in knowing the kind of association amongst the variables. This can be shown by the graph of  versus  One of the possible models of association is the simple linear regression, in which the elements of the A matrix give an explanatory variable and the elements of the B matrix a response variable, so that,

where  e  are parameters to be estimated and  is the error associated with the response assumed to be Gaussian, independently and identically distributed. This assumption is a pragmatic approach avoiding more complex structures for the error term which would require further modeling assumptions we wish to avoid at this exploratory stage. Also, more complex forms of spatial dependence than given by the linear relation can, in principle, also occur. Our approach is to rely on simple assumptions for the randomization tests and leaving more complex structures to be considered by the model discussed in the next Section.
In this work, the randomization test for spatial pattern was carried out on the observations for each sampling date. The test can be extended for the detection of time patterns. However, this raises the question of how to combine the information from several units of observation. Although such alternative has been studied for the data on thrips occurrences, it was decided not to include the results here because of the small number of observations in time and the lack of a specific interest in testing for time patterns.

Modeling the spatial pattern
Having detected a spatial pattern, it may be of interest to describe the pattern by means of a stochastic model. Modeling allows not only the characterization of the dependence pattern but also for the prediction of quantities of interest such as a map of expected levels of infestation over the area, the proportion of the area with infestation above or below a certain threshold, and areas with high and low infestation, among others possible quantities of interest.
One possible way of modeling the spatial distribution is by adopting the geostatistical framework, which associates the level of spatial dependency with distances between sampled plots. Usually the description of the spatial dependence assumes that the closest sampled plots are more alike than those farthest apart (Montagna, 2001). Diggle et al. (2003) uses the term geostatistics to identify a part of the spatial statistical methods in which the used model describes a continuous variation of the observations over the space.
The basic geostatistical data format is  i=1, 2, ..., n, in which  identifies the spatial location, generally in the two dimensions and  is the measure of interest at the  position of the ith observation. The response variable can be potentially measured at any point within the studied region (Diggle & Ribeiro Jr., 2007).
The geostatistical model is specified by assuming two processes over the study region (Diggle et al., 1998; Diggle & Ribeiro Jr., 2007) described as follows.  is a measure process within the study region A which is observed at  a set of locations  to obtain the s, the observed data. This first process is related to a underlying Gaussian process  with mean  variance  and correlation function , where  is the distance between pairs of observations. The values of  are usually not directly observed. Conditional independence is assumed in the sense that the  are independent, conditionally on the values of  meaning all the spatial dependency is modeled through . The exact form of the relation between the two processes may vary according to the type of variable being measured. For instance, when  follows the Gaussian distribution, the model can be written as , in which the  values are mutually independent and follow the normal distribution, with mean 0 and variance  . In this case the observations yi can be seen as a noisy version of  at the location  and, for a finite set of plots, the random vector  follows a multivariate Gaussian distribution. More generally,  may follow other distributions and Diggle et al., (1998) specify a model within the class of the generalized linear model (McCullagh &  Nelder, 1989) in which the S process defines random effects with spatial dependence structure. Diggle and Ribeiro Jr. (2007) call this a generalized linear geostatistical model (GLGM). This model allows the explicit specification of a Poisson distribution for the observations, which is compatible with the insect counting structure of the data considered here. 
The GLGM is a special case of a mixed generalized linear model, in which the , i=1, 2, …, n are conditionally independent given  with expected values given by  and linear predictor  i=1, 2, …, n with a known link function , which, for the Poison model considered here is typically given by the logarithm function. The model can extended allowing for covariates considering , in which  is the observed covariate values and  is the regression parameter vector. (Diggle et al., 1998, Diggle et al., 2003). 
Let  be the observed total number of insects with a Poisson distribution with mean , i=1, 2, …, n in which  represents the number of leaves. Then the probability function is given by

	 The likelihood function is often considered for inference about the model parameters within the context of generalised linear models. However, in this case the likelihood function does not have a closed form and is given by 

where  is the correlation matrix for  and with dimension equal to the number of observations which cannot be solved by analytical or numerical methods. Each element of  is given by the corresponding value of the correlation function of the  process and therefore having model parameters within non-linear functions which explains the lack of such solutions. A possible solution is to use Monte Carlo Markov Chain (MCMC) methods and a computational implementation is available through the package geoRglm (Christensen & Ribeiro Jr., 2002) for the R statistical environment (R Development Core Team, 2007). 
For discrete random variables, the variogram is not a natural summary of the data, but it may be useful as a diagnosis tool, after fitting the mixed generalized linear model (Diggle & Ribeiro Jr., 2007). In this case, the variogram obtained from the estimated parameters can be compared to the experimental variogram, obtained through the residuals from a GLM model fit. The variogram if given, respectively, by

which can be written as 

However, this approach must be used with caution because the variogram is even more erratic then the one usually obtained for data with a symmetric and continuous distribution, because of the asymmetric data.
After the choice of a specific model, a map that describes the behaviour of the study variable over the region can be obtained. Supposing that the parameters are known and that the interest is in the expected insects number given by , for the location  from the  marginal distribution and the  conditional distribution, it is possible to simulate the conditional distribution of , using the MCMC method. The predicted surface is given (Diggle et al., 1998) by 

where  is the process mean in this case because there are no explanatory variables or trend, and  is the linear kriging predictor and Var(x) is the prediction variance.

Results and discussion
Spatial pattern detection through Mantel’s randomization test
	The Mantel’s test was applied separately for each sampling date and each farm and the obtained p-values contrasted with the adopted 5% significance level. For the Fazenda São Paulo, there was evidence of spatial pattern in the number of insects per leaf for the first three data collections on 10th , 24th and 31th of July. These patterns can be observed in the dispersion plots shown in Figure 3 where symbols sizes are proportional to the number of insects per leaf. In general, considering all the farms and dates, the distribution of the mean number of insects per leaf is asymmetric and, does not show trends against the spatial coordinates. Also, the linear regression between the number of insects by leaf distances and the stakes location distances shows that, for the above mentioned dates, there is evidence of positive association in conformity with Table 2 which shows, as well, analogous results for the dates that detected spatial pattern at the other farms.