The unilateral spatial autogressive process for the regular lattice two-dimensional spatial discrete data

SORT 45 (1) January-June 2021, 67-90 DOI: 10.2436/20.8080.02.110 The unilateral spatial autoregressive process for the regular lattice two-dimensional spatial discrete data Azmi Chutoo1, Dimitris Karlis2, Naushad Mamode Khan1 and Vandna Jowaheer3 Abstract This paper proposes a generalized framework to analyze spatial count data under a unilateral regular lattice structure based on thinning type models. We start from the simple spatial integervalued auto-regressive model of order 1. We extend this model in certain directions. First, we consider various distributions as choices for the innovation distribution to allow for additional overdispersion. Second, we allow for use of covariate information, leading to a non-stationary model. Finally, we derive and use other models related to this simple one by considering simplification on the existing model. Inference is based on conditional maximum likelihood approach. We provide simulation results under different scenarios to understand the behaviour of the conditional maximum likelihood. A real data application is also provided. Remarks on how the results extend to other families of models are also given. MSC: 62H11, 62M30. Keywords: Unilateral, Spatial, Regular, Lattice, Thinning. 1 Introduction Problems with spatial count data occur in several disciplines. For example, consider the Human West Nile virus counts spread (Tevie, Bohara and Valdez, 2014), the infant death syndrome for the counties in North Carolina (Cressie and Chan, 1989), the number of vehicle burglary incidents in counties of Texas (Chun, 2014) and most recently, the cases and/or deaths of COVID-19 outbreak. Spatial data usually viewed as an aggregation or average of the events of interest emanates from a lattice structure. There are two main broad ways of representing the spatial observations (Cressie, 1993). The first, and most common way, is when the ob- 1Department of Economics and Statistics, University of Mauritius, Mauritius. 2Department of Statistics, Athens University of Economics and Business, Greece. 3Department of Mathematics, University of Mauritius, Mauritius. Received: December 2020 Accepted: April 2021 68 Spatial models for discrete data servation is in the form of a single indexed variable obtained from an areal unit kwithin some defined boundaries and in which case the spatial event is denoted as Yk, for the kth areal unit. In the second way, the spatial count observation is indexed in two dimensional forms, in terms of the location or site coordinate (or also termed as the latitude and longitude position (i,j)), denoted as Yi j, situated over a regular or irregular grid. Under both representations, the spatial count observation is supplemented by a neighbourhood structure, defined by terms of areal units or sites within the lattice structure. The second form of representation is useful since the two-dimensional representation considers all border cells in the region of interest (see, e.g. Tjøstheim, 1978a; Basu and Reinsel, 1993). Details on these representations can be found in Cressie (1993, Chap 6). In the analysis of spatial data, it is important to investigate the spatial dependence between observations from the different neighbouring areal units or sites. Next, such analysis can also shed light on the possible factors or effects influencing the spatial observations and these can include variables such as the distance metric, elevation, slope, rock type and land use fault types (see, e.g. Tobler, 1969). Several models have been studied in the literature to analyze spatial processes. The majority of the literature treats spatial continuous or discrete processes involving areal units, while only few papers consider the spatial data with a coordinate system, especially for the discrete case. Besides, such spatial observations are seen mostly in the agricultural, disease mapping, environmental and in the field of criminology. Specifically, in agriculture, the plantation field is usually split into small areas or say, square grids or cells with location (i,j)and wherein each cell, the investigator is interested on the number of plants cultivated subject to factors influencing its cultivation (see Kruijer et al. (2007) and references therein) and along with how the plants in the different neighbouring cells impact on the harvest in the (i,j)th position. Similarly in epidemiology, researchers are often concerned on the factors influencing the number of infected or death cases as a result of an outbreak of a virus in a region and how this is affecting the number of cases in the neighbouring regions. Such data has been treated in Cressie and Chan (1989), Wakefield (2007) while some more examples can be found in Lawson et al. (1999) and Lawson and Williams (2001). Moreover, in environmental field, the occurrences of road traffic accidents at different segments also illustrate spatial data analysis. In fact, in a hotspot analysis conducted in Barcelona, it was shown, via the local Moran statistics, that road accidents are concentrated in close neighbouring areas that have a complex road network systems with large roundabouts (Alvarez, 2020). Some other related research include the works by Valverde and Jovanis (2008) and Satria and Castro (2016). Last but not least, Mburu and Bakillah (2016) reported on the number of vehicle burglary incidents in small neighbouring regions of London which were highly spatially autocorrelated. Their study also revealed several influential factors such as unemployment and crimes in these areas of London. Unsurprisingly, there is influx of models for the areal-type spatial data that include mainly the class of conditional autoregressive (CAR) models (Besag, 1974) and its extensions to Intrinsic CAR (ICAR) (Besag and Kooperberg, 1995), the Besag-Yorke- Azmi Chutoo, Dimitris Karlis, Naushad Mamode Khan and Vandna Jowaheer 69 Mollie (BYM) (Besag, Molli´e and York, 1991) models and among several other extended CAR-based models; a review can be found in Obaromi (2019). In fact, for the regular lattice data of discrete nature that are represented in terms of the site coordinates, the only works appearing in the literature so far use the spatial integer-valued auto-regressive model of order 1 (SINAR(1)) by Ghodsi, Shitan and Bakouch (2012) and Ghodsi (2015). The model was constructed by introducing dependence between the observation of interest with its unilateral spatial neighbouring observations via the binomial thinning operator defined by Steutel and Harn (1979). The structure is similar to the observation-driven integer-valued time series models defined in McKenzie (1986). Properties of the model including asymptotic properties of the CML estimators were thereon established and the spatial process was proven to be stationary and ergodic. In the present paper we extend the model in certain directions. Firstly, we introduce different distributions for the innovations to enlarge the model and allow for larger variance, usually observed in spatial data due to clustering effects. Secondly, we allow for further spatial information to be used in the form of covariates that affect the model, leading to a non-stationary model. For this new model we discuss inference based on the CML. Moreover we discuss and apply some models related to the basic one that are parsimonious and easier to interpret, while they allow for easier extension to a broader family of models. Throughout the paper, some computational issues arising are also discussed. The remaining of the paper proceeds as follows: The basic model and its extensions are described in Section 2. Simulations to further support the approach are provided in Section 3. A real data application related to the new models is provided in Section 4. Extensions of the current model and concluding remarks can be found in Section 5. 2 Generalised SINAR(1) model (GSINAR(1)) We consider spatial processes defined on a regular rectangular grid in two dimensions with sites labelled (i,j), with an associated random variable Yi j defined at each site. Examples of such phenomena include data collected on a regular grid of size n1×n2from satellites and from agricultural field trials. The unilateral model (see, e.g. Tjøstheim, 1978b, 1983; Basu and Reinsel, 1993) defines the neighbouring sites that provide information for the site (i,j), namely we denote as Si j the set of indices (k,ℓ)of sites that are considered as neighbours of the site (i,j)and we define this as Si j ={(k,ℓ)∈Z2:k≤i,ℓ ≤j}−{(i,j)}. Tjøstheim (1983) described the model of order (p1,p2)for continuous data as: Yi j = p1 ∑ k=0 p2 ∑ ℓ=0 φkℓYi−k,j−ℓ+ǫi j, 70 Spatial models for discrete data Yi−1,j−1Yi−1,jYi−1,j+1 Yi,j−1Yi,jYi,j+1 Yi+1,j−1Yi+1,jYi+1,j+1 Figure 1:A diagram representing the unilateral model. The arrows indicate which site influence the site under consideration. Only sites located at the left and upper from the site into consideration influence the site. where φ00 =0 and the errors ǫi j follow a N(0,σ2)distribution. The model of order p1=p2=1 is described in Basu and Reinsel (1993). In the present paper we restrict our interest on this spatial structure, i.e. only of order 1. The assumed structure can be seen in Figure 1. One can see that for each point, only three neighbouring points coming from the upper left direction affect the point, implying a restricted spatial structure. In other words the spatial effect is propagated from left to the right and from top to the bottom only. 2.1 The stationary model Ghodsi et al. (2012) extended the model for count data. Hence, they defined that the spatial observation located at site coordinate (i,j)follows an auto-regressive equation of the form: Yi j =α1◦Yi−1,j+α2◦Yi,j−1+α3◦Yi−1,j−1+ǫi j,(1) where i=1,...,n1,j=1,...,n2. The dependence of Yi j on its neighbours as defined by set Si j is handled in equation (1) through the binomial thinning operator “ ◦”. The binomial thinning mechanism emanates from the work of Steutel and Harn (1979) (see also Scotto, Weiß and Gouveia, 2015) for a summary of such operators) and is expressed as: α◦Y= Y ∑ s=1 Bs(α),(2) where α∈[0,1], and Bs(α),s=1,...,Yare identically and independently distributed Bernoulli r.v with P(Bs(α) = 1) = 1−P(Bs(α) = 0)=α. In parsimony, we impose the assumption of independent thinning operator (Du and Li, 1991; Bu, McCabe and Hadri, 2006). In (1), {ǫi j}i=1,...,n1,j=1,...,n2represents the corresponding innovation sequence of independent non-negative integer-valued random variables with finite mean λǫand finite variance τ2 ǫand has a distributional form as Pǫi j (·). Furthermore, at any position (i,j), ǫi j is assumed to be independent of all Yi−k,j−ℓ. In this simple form, the model in (1) is stationary if Azmi Chutoo, Dimitris Karlis, Naushad Mamode Khan and Vandna Jowaheer 71 α1+α2+α3<1.(3) Properties of this SINAR(1) model as well as estimation can be found in Ghodsi et al. (2012) and Ghodsi (2015). The process in equation (1) is proven to be ergodic in Ghodsi (2015) and Markovian in Pickard (1980). 2.2 The non-stationary model Here, we extend the model to the non-stationary case by allowing site specific covariates to influence the mean of the innovation process. We denote the mean and variance of the innovation term as λi j and τ2 ǫi j respectively. We further consider λi j as a function of some position-variant and invariant covariates i.e. λi j =f(xT i jβ β β)with xi j = [xi j1,...,xi jp]Tand regression coefficients β β β= [β β β0,...,β β βp]T. Note that a log function is a standard choice for such cases leading to logλi j =xT i jβ β β. We name this model as Generalized SINAR of order 1 (GSINAR(1)). From the following binomial thinning properties, E(α◦Y) = αE(Y) V(α◦Y) = α(1−α)E(Y) +α2V(Y) Cov(α1◦Y1,α2◦Y2) = α1α2Cov(Y1,Y2),αj∈(0,1),j=1,2 and then we get the unconditional expectation of Yi j to be E(Yi j) = µi j =α1µi−1,j+α2µi,j−1+α3µi−1,j−1+λi j.(4) For the variance we get V(Yi j) = σ2 i j =α1(1−α1)µi−1,j+α2 1σ2 i−1,j +α2(1−α2)µi−1,j+α2 2σ2 i,j−1 +α3(1−α3)µi−1,j−1+α2 3σ2 i−1,j−1 +2α1α2Cov(Yi−1,j,Yi,j−1) +2α1α3Cov(Yi−1,j,Yi−1,j−1) +2α2α3Cov(Yi,j−1,Yi−1,j−1)+τ2 ǫi j . (5) By letting γ(k,ℓ)=Cov(Yi−k,j−ℓ,Yi j), we obtain a difference equation of the form: γ(k,ℓ) = α1γ(k−1,ℓ)+α2γ(k,ℓ−1)+α3γ(k−1,ℓ−1).(6) 72 Spatial models for discrete data Closed form expressions for the above moments are difficult to obtain under non-statio- nary conditions. Whilst, unless assuming weak-stationarity, that is, ǫi j has constant mean λǫand variance τ2 ǫ, we obtain simple expression for the mean and variance with the covariances obtained by solving the difference equation γ(k,ℓ)using the approach in Basu and Reinsel (1993). However, the derivation of the covariance structure in equation (6) is not required in the estimation of the model parameters when using the conditional maximum likelihood equation illustrated as follows. Conditional on the neighbourhood Si j, the probability mass function for the GSINAR(1) model is given by: P(Yi j|Si j) = ∑R1 s1=0∑R2 s2=0∑R3 s3=0pα1(s1|Yi−1,j)pα2(s2|Yi,j−1)pα3(s3|Yi−1,j−1) Pǫi j (Yi j −s1−s2−s3),(7) where R1=min{Yi−1,j,Yi j},R2=min{Yi,j−1,Yi j −s1}and R3=min{Yi−1,j−1,Yi j −s1− −s2}and pα(s|Y) = Y sαs(1−α)Y−s;s=0,1,...,Y, i.e. the probability mass function of a binomial distribution. In the present paper we have considered different choices for Pǫi j (·). The standard assumption of a Poisson distribution limits the variability we expect in the data (see Appendix). A natural improvement is to consider mixed Poisson alternatives. We consider the negative binomial, Poisson-Inverse Gaussian as well as the Poisson-Lindley in order to allow for quite different effects. Also, we consider the COM-Poisson distribution in order to allow for a general model which accounts for underdispersion if we need so. We postpone the details until section 3. Therefore, the log conditional maximum likelihood (CML) equation is then given by: ℓ(θ θ θ) = logL(θ θ θ) = n1 ∑ i=2 n2 ∑ j=2 logP(Yi j|Si j),(8) where θ θ θ=[α1,α2,α3,β β β,ν],νrefers to the dispersion parameter of the innovation distribution if it exists and β β βis the vector of regression coefficients for the mean of the innovation. It can be seen in (Ghodsi, 2015) that ˆ θ−θ∼N(0,I−1(θ)), where I(θ)is the Hessian matrix. The CML equation in (8) is then maximized. Some computational details are the following. We have used the optim function in R. Note that the conditional distribution needs to derive the convolution of three binomials plus the distribution of the innovation term. This can be computationally intensive. We have reduced the computational burden by observing that the probabilities of the binomial distribution are just the coefficients of a polynomial of order Nwhere Nis the number of trials in the binomial. As such computing the convolution of two binomial is equivalent to multiply two polynomials for which there are very fast procedures, like those in the library pracma in R. This reduced the computational effort and improved with respect to the errors produces by huge finite summations. Overall, maximization of (8) was rather simple even for complicated innovation distributions. Azmi Chutoo, Dimitris Karlis, Naushad Mamode Khan and Vandna Jowaheer 73 2.3 Related models The general model in (1) has three parameters to introduce the spatial correlation, namely α1,α2and α3that described the vertical, horizontal and diagonal dependence respectively. One may eliminate some of the effects by setting the corresponding αparameter equal to 0. For example setting α3=0 we assume no-upper diagonal effect, while setting α2=0 we assume no horizontal effect. Such submodels can be very useful in order to examine and interpret the underlying situation for a dataset. For example, we may test and recognize which effect the vertical or horizontal is more important. Another way to simplify the model is by assuming one common effect using the same parameter αfor all directional relationships. Such a model takes the following form due to the properties of thinning operators: Yi j =α◦∑ (k,ℓ)∈Si j Yk,ℓ +ǫi j.(9) The model assumes that all neighbouring sites contribute the same to the structure. Such a model resembles simple INAR(1) time series models. It has the advantage of having less parameters to estimate and explain; all neighbours contribute the same. On the other hand, this may be restrictive since the spatial effects may differ due to direction and thus the model may fail to capture them correctly. Model (9) allows for easy extensions to a general neighbouring structure. It is evident that by considering the set Si j defining the neighbouring sites, this model can be generalized to a large extend including the non-regular lattice case which is more realistic in many applications. The model just assumes that all neighbours contribute to the observation at hand. Properties of such models as well as estimation is straightforward based on the results of the current section. Finally, in the present paper we assume that the covariates enter in the model by the mean of the innovations. One may consider that spatial correlation parameters αj may relate to some covariate information through a logit link function. For example, we can assume that logit(α1i j)=xT i jδ δ δ1, where xi j is some covariate information for the site (i,j)and δ δ δ1some vector of regression coefficients. In this case we assume that each point in space has a different spatial effect α1depending on some characteristics xi j. For example, we may assume that altitude can change the spatial effect, which makes sense if we measure for example something which can be altered due to wind conditions. We believe that such a model, while it has some potential is special cases, it can complicate the model interpretation, especially if we have regression effects in both the mean and the autocorrelation parameters. 74 Spatial models for discrete data 3 Simulation study In this section, we perform simulation experiments, using equation (1), with different innovation distributions namely the Poisson, Negative Binomial (NB), Poisson- Lindley (PL), Conway-Maxwell Poisson (COM-Poisson / CMP) and Poisson-Inverse Gaussian (PIG) with rate or mean parameters commonly indicated by a link predictor λi j = exp(xT i jβ β β)and dispersion index ν. As it is described in more details in the Appendix, using a distribution for the innovations that allows over/under dispersion, we also extend such properties to the observed spatial distribution and hence more realistic and flexible fitting can be achieved. In the present paper we consider the following distributions for the innovations. •For Poisson innovations we assume that P(ǫi j) = e−λi j λǫi j i j ǫi j!;ǫi j =0,1,...;λi j ≥0. •For NB innovations we use the following parametrization: P(ǫi j) = Γ(ν+ǫi j) Γ(ν)ǫi j!λi j λi j +ννν λi j +νǫi j ;ǫi j =0,1,...;λi j ≥0,ν≥0. For this parametrization the mean is λi j and the variance λi j +λ2 i j/ν. •For PL innovations we use P(ǫi j) = λ2(ǫi j +λi j +2) (λi j +1)ǫi j+3;ǫi j =0,1,...;λi j >0. The mean is (λi j +2)/(λi j(λi j +1)) while the variance is λ3 i j +4λ2 i j +6λi j +2 λ2 i j(λi j +1)2 PL can have different shapes than the other Poisson mixtures like the NB and PIG models. •For COM-Poisson innovations we use P(ǫi j) = λǫi j i j (ǫi j)!ν 1 Z(λi j,ν);ǫi j =0,1,...;j=0,1,...;λi j ≥0;ν≥0, where Z(λi j,ν)=∑∞ j=0 λj i j j!ν. Note that λi j is not the mean of the distribution; the mean is hard to be written in closed form, but it is approximated by λ1/ν i j . Azmi Chutoo, Dimitris Karlis, Naushad Mamode Khan and Vandna Jowaheer 75 •For PIG innovations we use the parameterization from the package actuar in R. Namely the pmf is defined with parameters λand dispersion νas P(ǫi j) = 2ν π1/2 expν λi j ai j ν−(x−1 2) Kx−1 2(ai j),x=0,1,..., λi j,ν>0, where a2 i j =2ν1+ν 2λ2 i j and Kx(a)is the modified Bessel function of the third kind. The mean is λi j and the variance is λi j +λ3 i j/ν. The choice of the distributions for the innovation term attempts to cover a range of possible models. So, we have used the Poisson distribution as a starting point, two of the most famous mixed Poisson ones (negative binomial and Poisson-Inverse Gaussian), the COM-Poisson to allow for under-dispersion as well and finally the Poisson-Lindley since this is a tractable mixed Poisson distribution with very different shapes. In order to simulate the grid we followed the following approach: We added an additional row and column with all values equal to 0, i.e. we set Y0j=Yi0=0 for all iand j. Then we simulated the grid Yi j,i=1,...,n1+10 and j=1,...,n2+10 based on the model in (1) using the chosen innovation distribution. Then we rejected the rows and columns from 0 up to 10 so as to keep the grid n1×n2. 3.1 Numerical Results: No covariates For scenario 1, the simulation study assumes the following combinations of (α1,α2,α3, λ,ν)and grids: 1. C1: (0.35, 0.15, 0.2, 5, 0.5) and grid 25 x 25 2. C2: (0.25, 0.25, 0.3, 3, 0.8) and grid 40 x 40 3. C3: (0.6, 0.2, 0.15, 7, 2) and grid 50 x 50 Note that for Poisson-Lindley and COM-Poisson cases parameter λis not the mean while νhas a different interpretation. So, in the simulations this is the value used to simulate the data. For Poisson, negative binomial and Poisson-Inverse Gaussian, λis the mean and νis the dispersion parameters, equal to 1 for the Poisson. Obviously for ν→∞we get the Poisson distribution in such cases. For each scenario 1000 replications were obtained. The simulated mean estimates, their biases, root mean square errors (RMSEs) and standard deviations (SDs) are reported. The results are displayed in Tables 1, 2 and 3. 82 Spatial models for discrete data Table 6:Different models with negative binomial innovation distribution fitted to the Bei data. Models assumes different (or not ta all) spatial dependence. The full model with all three kind of unilateral effects is the chosen one. Model Param Log-lik AIC No restriction 5 −2666.22 5337.442 M1: α1=α2=α33−2686.20 5375.406 M2: α3=0 4 −2674.75 5353.490 M3: α2=0 4 −2765.51 5535.024 M4: α1=0 4 −2717.78 5439.554 M5: α2=α3=0: only verctical 3 −2802.29 5607.586 M6: α1=α3=0: only horizontal 3 −2746.06 5495.128 M7: α1=α2=0: only diagonal 3 −2840.98 5684.954 no spatial effect α1=α2=α3=0 2 −2938.80 5879.606 Table 6 presents for the chosen negative binomial case some more spatial scenarios as mentioned in section 2.3. For example, model M1 assumes that all three α’s are the same, while models M2 to M4 that one of the spatial correlations is not present i.e. we set αj=0 for j=1,2,3 respectively. Actually we remove each time the vertical (M4), horizontal (M3) and the diagonal effects (M2). Finally, models M5 to M7 suggest that only one spatial effect suffices. The full model with all three kind of unilateral effects is the chosen one as judged by AIC, revealing the underlying structure of the data. One can see that the horizontal effect is larger as judged by the change in the LRT when we remove each effect. In addition we use some covariate information available. The Bei data set is accompanied by covariate data giving the elevation (altitude) and slope of elevation in the study region. An important question arises is whether the intensity of Bei trees may be viewed as a spatially varying function of the covariates. We have fitted different models to examine the improvement offered by the covariates. Both covariates were found significant. The selected model can be seen in Table 7. Table 7:Comparison between Poisson innovation and mixed Poisson innovations using Akaike information criterion (AIC): Application to Bei count data with both covariates. Innovations ˆα1ˆα2ˆα3ˆ β0ˆ β1ˆ β2ˆνAIC Poisson Estimates 0.187 0.286 0.102 -6.265 0.037 8.128 s.e (0.012) (0.013) (0.011) (0.671) (0.004) (0.478) 6940 NB Estimates 0.146 0.219 0.055 -6.897 0.044 8.515 0.242 s.e (0.016) (0.017) (0.016) (1.507) (0.010) (1.327) (0.089) 5293 PL Estimates 0.115 0.208 0.000 4.904 -0.029 -6.839 s.e (0.016) (0.017) (0.014) (0.645) (0.004) (0.546) 5716 CMP Estimates 0.111 0.202 0.006 -2.162 0.010 2.089 0.001 s.e (0.016) (0.018) (0.015) (0.116) (0.001) (0.236) (0.031) 5626 PIG Estimates 0.146 0.221 0.056 -15.053 0.097 16.539 0.210 s.e (0.016) (0.017) (0.015) (2.395) (0.016) (3.451) (0.0213) 5294 Azmi Chutoo, Dimitris Karlis, Naushad Mamode Khan and Vandna Jowaheer 83 It is observed from Table 7 that the SINAR(1) model using NB innovation distribution with covariates still yields the lowest AIC value and hence outperforms the other models with different innovations by producing much better fit to the data. PIG is very competitive to the NB. We focus on the selected model with negative binomial innovations. From the results of Table 7 one can see that both covariates are statistically significant, with positive sign, hence increasing the altitude and the slope of elevation we obtain an increased number of trees, having adjusted for the effect of neighbouring areas. From the α’s we see that the larger effect comes from α2that measures the horizontal dependence. All spatial effects are statistically significant at 5%. The model implies a clear spatial dependence. 4.3 Goodness of Fit In order to jude whether the fitted model is satisfactory we have worked a few ideas. To start with, we derived the one step ahead predictions based on the models. Namely, we derived for each data point E(Yi j|Si j) = ˆα1Yi−1,j+ˆα2Yi,j−1+ˆα3Yi−1,j−1+ˆ λi j and log ˆ λi j =ˆ β0+ˆ β1X1i j +ˆ β2X2i j Observed Expected value (Fitted) 0 10 20 30 40 50 60 70 Figure 5:The observed counts in the 40 ×40 grid and the fitted based on the model. Fitted values are expected values based on the conditional expectation. 84 Spatial models for discrete data using the estimated parameters from the selected negative binomial model. The values can be seen in Figure 5 together with the observed counts. We emphasize that the predictions are the expected means that is why they cannot capture the extreme values. From the plots one can see that the model captures in a great extend the pattern. x y 1 10 20 30 40 40 30 20 10 1 Points were predicted value is outside the 95% interval Figure 6:Points in the grid that the true value was outside the 95% prediction interval created. To further exploit the quality of the predictions we have created for each data point (i,j)a 95% confidence interval for the prediction. To do so, we simulated 1000 values from the predictive probability mass function as provided in (7) and then based on them we created the intervals. Out of the 1521 values we predict (we did not predict the first row and column) only 33 (2.1%) values were outside of the interval, implying that the model was quite satisfactory. The values that lay outside the interval are depicted in Figure 6. One can see again that we have failed to predict some of the extreme values as one can see compared to the Figure 5. Another important aspect of the model fitting lies on the ability of the model to capture the spatial dependence structure. To check this aspect we simulated grids of the same size 40 ×40 from the selected negative binomial using the estimated parameters. For each simulation we have estimated the spatial covariance at lags kand ℓby ˆγ(k,ℓ) = 1 n 40 ∑ i=k+1 40 ∑ j=ℓ+1 (Yi j −¯ Y)(Yi+k,j+ℓ−¯ Y) Azmi Chutoo, Dimitris Karlis, Naushad Mamode Khan and Vandna Jowaheer 85 and then we derived the spatial correlation at lag kand ℓas ˆρ(k,ℓ) = ˆγ(k,ℓ) s2 Y , where ¯ Yand s2 Yare the mean and the variance estimated from the data. We have used lags k=0,1,2 and ℓ=0,1,2 and then we compared them with those values observed from the data in order to see whether the observed dependence structure could have been created by the model at hand. We show in Figure 7 the 95% confidence intervals created by 1000 simulations and the dot indicates the observed value. We see a good agreement. The two first values need perhaps improvement with a richer model but overall the model captures the underlying structure in a reasonable way. 0.0 0.1 0.2 0.3 0.4 0.5 spatial correlation (k,m) value (1,0) (0,1) (1,1) (2,0) (0,2) (2,1) (1,2) (2,2) Figure 7:Points in the grid that the true value was outside the 95% prediction interval created. As a summary, we believe that the current model fits satisfactorily the observed spatial structure. 5 Conclusion This paper revisits and extends the simple stationary SINAR(1) model introduced by Ghodsi et al. (2012) and Ghodsi (2015). The SINAR(1) model is the first research in the modelling of the two dimensional unilateral spatial discrete data based on the thinning mechanism that allows to model explicitly the discrete nature of the data. 86 Spatial models for discrete data In the present paper we proposed some novel extensions of the existing SINAR(1) model. These novelties in fact overcome two important limitations of the simple SINAR(1). Firstly, in our model specification, we propose to model the data using overdispersed innovation distributions, while simultaneously allow covariate information to be used leading to a non-stationary model. While not treated in this version, one may also use offset in the regression part, like population size values, if needed. We also discuss parsimonious representations of the model at hand. The model parameters are estimated using the CML approach. We acknowledge some restrictions of the current model which we consider to improve. The proposed model is based on the unilateral regular lattice case. One can extend the model to capture several other cases. For the regular lattice case, define the set of indices (k,ℓ)of the neighbouring observations for the (i,j)observation as Si j. Then, in the general setting, the model can be written as for a general neighbourhood: Yi j =∑ (k,ℓ)∈Si j αkℓ◦Yk,ℓ +ǫi j, where are usual, the ǫi j are the innovations. Defining appropriately the sets Si j one can derive other models at the expense of parsimony. Finally consider the typical case in spatial data where observations are indexed simply as Yjto indicate the value at site jfrom a map with j=1,...,nsites, as for example the different regions of a country. Define as Sjthe indices of its neighbours. In such case the model of order 1 can have the form: Yj=α◦∑ k∈Sj Yk+ǫj or equivalently if we define the n×nadjacency matrix Wwith elements wi j with values equal to 1 if the sites iand jare neighbours and 0 otherwise, then we can write the models as: Yj=α◦∑ k6=j wk jYk+ǫj to mimic typical order 1 models for spatial continuous data. Such generalization will be reported elsewhere. Also note that in this paper we used only spatial model of first order. One may consider SINAR(p)models with higher order effects. Such extension needs special care. It is already known that simple INAR(p)models can have different interpretations /representations, (see the different approaches in Alzaid and Al-Osh (1990) and Jin-Guan and Yuan (1991)). Extending to a SINAR(p)model can have a large number of parameters making inference quite complex. Perhaps more parsimonious models like the one in Section 2.3 are easier to extend to higher order. Azmi Chutoo, Dimitris Karlis, Naushad Mamode Khan and Vandna Jowaheer 87 A Appendix section Based on Ghodsi et al. (2012), for the case of a stationary model, we have for the marginal stationary mean µYand the stationary variance σ2 Ythat µY=µǫ 1−α1−α2−α3 and σ2 Y= µY 3 ∑ i=1 αi(1−αi)+τ2 ǫ 1−(α1+α2α3)λ−(α2+α1α3)η−α2 3 where η=α2+α3λ 1−α1λand λ=(1+α2 1−α2 2−α2 3)−q(1+α2 1−α2 2−α2 3)2−4(α1+α2α3)2 2(α1+α2α3) In the formulas, µǫand τ2 ǫare the mean and the variance of the innovations respectively. Note a misprint in Ghodsi et al. (2012) for the variance. Define the index of dispersion IDY=σ2 Y/µY. Dividing the variance with the mean, we get for the index of dispersion of the spatial data that IDY=σ2 Y µY = 3 ∑ i=1 αi(1−αi)+IDǫ(1−α1−α2−α3) 1−(α1+α2α3)λ−(α2+α1α3)η−α2 3 which relates directly the index of dispersion of the innovation distribution IDǫto that of the marginal, i.e. IDY. Since the denominator is positive and all the quantities in the nominator also are positive, an increase of IDǫwill lead to increase of the IDY. Thus assuming an overdispersed distribution for the innovations we can have much larger overdispersion in the observed spatial data. One can see that even for the Poisson innovations the index of dispersion is larger than 1, however for reasonable values for counts this overdispersion is limited. The introduction of overdispersed innovations increase a lot the overdispersion as one can see in Figure 8. In Figure 8 the two axes depict the marginal mean and variance for a stationary model given above. The different lines correspond to different levels of overdispersion for the innovation distribution. We have used α1=α2=α3=0.2. The diagonal line refers to the case of equidispersion. Therefore, above that line we get overdispersion and below underdispersion. The red line (ID=1) corresponds to Poisson innovations. One can see that in this case we get small overdispersion for the spatial case. Increasing the overdispersion on the innovation, as for example considering a mixture of Poisson we get larger overdispersion. Note that an underdispersed innovation distribution, like the cases of COM-Poisson distribution, can lead to underdispersion. 88 Spatial models for discrete data 0 2 4 6 8 10 12 0 2 4 6 8 10 12 marginal mean marginal variance ID=0.5 ID=1 ID=1.2 ID=1.5 equidispersion line ID=2.5 Figure 8:The marginal mean and variance for a stationary model. The different lines correspond to different levels of overdispersion for the innovation distribution. for the plot α1=α2=α3=0.2. The diagonal line refers to the case of equidispersion. ID implies the index of dispersion of the innovation distribution. One can see that for Poisson (red line) we get small overdispersion. Increasing the overdispersion of the innovations lead to increased overdispersion for the spatial distribution. References Alvarez, M. (2020). Predicting traffic accident hotspots with spatial data science. https://carto.com/blog/ predicting-traffic-accident-hotspots-with-spatial-data-science/ Alzaid, A. and Al-Osh, M. (1990). An integer-valued pth-order autoregressive structure (INAR(p)) process. Journal of Applied Probability, 314-324. Baddeley, A., Rubak, E. and Turner, R. (2015). Spatial Point Patterns: Methodology and Applications with R. London: Chapman; Hall/CRC press. Basu, S. and Reinsel, G. C. (1993). Properties of the spatial unilateral first-order ARMA model. Advances in Applied Probability, 631-648. Besag, J. (1974). Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society. Series B, 36, 192-236. Besag, J. and Kooperberg, C. (1995). On conditional and intrinsic autoregressions. Biometrika, 82, 733-746. Besag, J., Molli´e, A. and York, J. (1991). Bayesian image restoration, with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics, 43, 1-20. Bu, R., McCabe, B. and Hadri, K. (2006). Maximum likelihood estimation of higher-order integer-valud autoregressive processes. Journal of Time Series Analysis, 29, 973-994. Chun, Y. (2014). An application to vehicle burglary. Geographical Analysis, 46, 165-184. Cressie, N. (1993). Statisics for Spatial Data (2nd ed.). New York: Wiley. Cressie, N. and Chan, N. (1989). Spatial modeling of regional variables. Epidemiology and Infection, 84, 393-401. Azmi Chutoo, Dimitris Karlis, Naushad Mamode Khan and Vandna Jowaheer 89 Du, J. and Li, Y. (1991). The integer-valued autoregressive INAR(p) model. Journal of Time Series Analysis, 12, 129-142. Ghodsi, A. (2015). Conditional maximum likelihood estimation of the first-order spatial integer-valued autoregressive SINAR(1,1) model. Journal of The Iranian Statistical Society, 14, 15-36. Ghodsi, A., Shitan, M. and Bakouch, H. (2012). A first-order spatial integer-valued autoregressive SINAR(1,1) model. Communications in Statistics- Theory and Methods, 41, 2773-2787. Jin-Guan, D. and Yuan, L. (1991). The integer-valued autoregressive (INAR(p)) model. Journal of Time Series Analysis, 12, 129-142. Kruijer, W., Stein, A., Schaasfsma, W. and Heijting, S. (2007). Analyzing spatial count data, with an application to weed counts. Environmental and Ecological Statistics, 14, 399-410. Lawson, A., Biggeri, A., Bohning, D., Lesaffre, E., Viel, J. and Bertollini, R. (1999). Disease Mapping and Risk Assessment for Public Health. Chichester, UK: Wiley. Lawson, A. and Williams, F. (2001). An Introductory Guide to Disease Mapping. New York: Wiley Medical Sciences. Mburu, L. and Bakillah, M. (2016). Modeling spatial interactions between areas to assess the burglary risk. International Journal of Geo-Information, 5, 47. McKenzie, E. (1986). Autoregressive moving-average processes with negative binomial and geometric marginal distrbutions. Advanced in Applied Probability, 18, 679-705. Møller, J. and Waagepetersen, R. (2007). Modern spatial point process modelling and inference (with discussion). Scandinavian Journal of Statistics, 34, 643-711. Obaromi, D. (2019). Spatial modelling of some conditional autoregressive priors in a disease mapping model: the Bayesian approach. Biomedical Journal of Scientific and Technical Research, 14. Pickard, D. (1980). Unilateral Markov fields. Advances in Applied Probability, 12, 655-671. Satria, R. and Castro, M. (2016). GIS tools for analyzing accidents and road design: A review. Transportation Research Procedia, 18, 242-247. Scotto, M., Weiß, C. and Gouveia, S. (2015). Thinning-based models in the analysis of integer-valued time series: A review. Statistical Modelling, 15, 590-618. Steutel, F. and Harn, K. V. (1979). Discrete analogues of self-decomposability and statibility. The Annals of Probability, 7, 3893-899. Tevie, J., Bohara, A. and Valdez, R. (2014). Examination of the geographical variation in human west nile virus: a spatial filtering approach. Epidemiology and Infection, 142, 2522-2529. Tjøstheim, D. (1983). Statistical spatial series modelling II: some further results on unilateral lattice processes. Advances in Applied Probability, 15, 562-584. Tjøstheim, D. (1978a). A measure of association for spatial variables. Biometrika, 65, 109-114. Tjøstheim, D. (1978b). Statistical spatial series modelling. Advances in Applied Probability, 10, 130-154. Tobler, W. R. (1969). Geographical filters and their inverses. Geographical Analysis, 1, 234-253. Valverde, J. and Jovanis, P. (2008). Analysis of road crash frequency with spatial models. Transportation Research Record Journal of the Transportation Research Board, 2061, 55-63. Wakefield, J. (2007). Disease mapping and spatial regression with count data. Biostatistics, 8, 158-183.