Assessing spatial dependency under non-standard sampling

Assessing Spatial Dependency under Non-Standard Sampling Raquel Menezes da Mota Leite July 2005 Departamento de Estat´ıstica e I.O. Universidad de Santiago de Compostela Vendo-os assim t˜ao pertinho A Galiza mai’lo Minho. S˜ao como dois namorados Que o rio traz separados Quasi desde o seu nascimento. Deixa-los, pois namorar J´a que os paes para casar Lhes n˜ao d˜ao consentimento Jo˜ao Verde Do livro “Ares da Raya” Acknowledgments I would like to thank my advisors Prof.Pilar Garcia-Soid´an and Prof.Manuel Febrero-Bande for all their support and direction to this thesis. From my visits to Lancaster, I thank the stimulating environment and advice provided by Prof.Peter Diggle and Prof.Jonathan Tawn. I also thank Inˆes, Rose and Susi that, among many others in Lancaster, helped me to find there a comfortable work and leisure environment. My appreciation is also due to Alan, Daniel, Chris and Jamie who struggled to eradicate the worst errors on my non-native english writings. I would like to acknowledge the combined financial support provided by Universidade do Minho, Prodep grant ref.5.3/N/189.015/01 and a Marie-Curie EU grant that together made this work possible. Finally I would to thank my Family and Carlos for their permanent support during this research. Contents 1 Introduction 1 2 Stationary spatial processes 5 2.1 Introduction............................... 5 2.2 Mean, covariance and variogram . . . . . . . . . . . . . . . . . . . . 8 2.3 Scaleofvariation ............................ 11 2.4 Further variogram properties . . . . . . . . . . . . . . . . . . . . . . 12 2.5 Trend and outlier identification . . . . . . . . . . . . . . . . . . . . 16 3 Comparison of valid variograms 21 3.1 Introduction............................... 21 3.2 Traditional three stages . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.2.1 Stage 1 – Empirical variogram estimation . . . . . . . . . . . 25 3.2.2 Stage 2 – Valid model selection . . . . . . . . . . . . . . . . 27 3.2.3 Stage 3 – Model fitting . . . . . . . . . . . . . . . . . . . . . 31 3.2.4 Existing combinations of the previous stages . . . . . . . . . 32 3.3 Simulationstudy ............................ 36 3.3.1 Comparing empirical estimators . . . . . . . . . . . . . . . . 37 3.3.2 Comparing complete approaches . . . . . . . . . . . . . . . . 40 3.3.3 Closingremarks......................... 49 4 Clustered and biased multi-stage sampling 51 4.1 Introduction............................... 51 i ii CONTENTS 4.2 Assessing through simulation . . . . . . . . . . . . . . . . . . . . . . 53 4.2.1 Sample generation algorithm . . . . . . . . . . . . . . . . . . 53 4.2.2 Impact on variogram estimation . . . . . . . . . . . . . . . . 56 4.3 Data exploratory methods . . . . . . . . . . . . . . . . . . . . . . . 60 4.3.1 Detection of dependence . . . . . . . . . . . . . . . . . . . . 61 4.3.2 Detection of sequential dependence . . . . . . . . . . . . . . 61 4.3.3 Impact of sample designs on Eseq,Vseq and V∗ seq ....... 63 4.4 MonteCarlotests............................ 66 4.4.1 Example of a simulated data set . . . . . . . . . . . . . . . . 69 4.4.2 Rongelap island’s data . . . . . . . . . . . . . . . . . . . . . 70 4.4.3 Randomization tests . . . . . . . . . . . . . . . . . . . . . . 72 4.5 Non-standard sampling correctors . . . . . . . . . . . . . . . . . . . 75 4.5.1 Method to adjust for clustering . . . . . . . . . . . . . . . . 76 4.5.2 Sequential biased corrector . . . . . . . . . . . . . . . . . . . 76 4.5.3 Results.............................. 77 4.5.4 Rongelap island’s data . . . . . . . . . . . . . . . . . . . . . 81 5 Properties of bγ(.)robust to clusters 83 5.1 Introduction............................... 83 5.2 Assumptions............................... 85 5.3 Neighbourhood radius selector . . . . . . . . . . . . . . . . . . . . . 86 5.4 Bias of bγ(.)robusttoclusters ..................... 93 5.4.1 Order of a1(u) for u≥Ch ................... 95 5.4.2 Order of a2(u)−a1(u)γ(u) for u≥Ch ............ 98 5.5 Variance of bγ(.) robust to clusters . . . . . . . . . . . . . . . . . . . 99 5.5.1 Order of e1(u) for u≥Ch ...................102 5.5.2 Order of e2(u) for u≥Ch ...................103 5.5.3 Order of e3(u) for u≥Ch ...................105 5.6 Kernel bandwidth selector . . . . . . . . . . . . . . . . . . . . . . . 108 5.6.1 Orderofvariance........................110 CONTENTS iii 5.7 Numericalstudies............................111 5.7.1 Performance of bγ(.) robust to clusters . . . . . . . . . . . . . 111 5.7.2 Analysis of Ed(n, a), Fd(n, a) and Gd(n, a) for large n....114 5.7.3 Estimates of E2(n, a), F2(n, a) and G2(n, a) .........117 6 Assessing preferential sampling 119 6.1 Introduction...............................119 6.2 Class of log-Gaussian Cox processes . . . . . . . . . . . . . . . . . . 121 6.3 Effect on variogram estimation . . . . . . . . . . . . . . . . . . . . . 124 6.3.1 Some simulation details . . . . . . . . . . . . . . . . . . . . 126 6.3.2 Influence of βon bias and variance . . . . . . . . . . . . . . 126 6.3.3 Clustered versus preferential . . . . . . . . . . . . . . . . . . 129 6.4 Impactonprediction ..........................135 6.4.1 Gaussiandata..........................136 6.4.2 Simulationstudy ........................138 7 Moss data and a model-based approach 143 7.1 Application to real data . . . . . . . . . . . . . . . . . . . . . . . . 144 7.1.1 Test if sample is preferential . . . . . . . . . . . . . . . . . . 149 7.1.2 Kriging and cross validation . . . . . . . . . . . . . . . . . . 153 7.2 A model-based approach . . . . . . . . . . . . . . . . . . . . . . . . 157 7.2.1 Relatedwork ..........................158 7.3 Likelihoodinference...........................160 7.3.1 Complete data likelihood . . . . . . . . . . . . . . . . . . . . 161 7.3.2 Likelihood of observed data . . . . . . . . . . . . . . . . . . 164 7.4 Estimation of model parameters . . . . . . . . . . . . . . . . . . . . 166 7.4.1 Profile log-likelihoods . . . . . . . . . . . . . . . . . . . . . . 167 7.4.2 Importance sampling . . . . . . . . . . . . . . . . . . . . . . 169 7.4.3 Direct Monte Carlo approximation . . . . . . . . . . . . . . 170 7.5 Simulationstudy ............................175 Chapter 1 Introduction Nowadays, spatial statistics plays an important role, as current technological development helps the derivation of spatial data. Besides the traditional sources of spatial data, such as maps, census material and aerial photography, more sophisticated and reliable data sources have appeared, such as remote sensors aboard satellites. These technologies resulting from the development of fast computers and specific software, such as image processing software and geographic information systems, are providing better tools and new demands from spatial statistics. Spatial models work with data collected from different spatial locations. These models measure the relationship between observations at various locations. They should reflect the intuitive hypothesis that data items close in space are correlated and that this correlation decreases when distance increases. They should also provide evidence of the existence of spatially correlated errors. The spatial correlation analysis enables us to see how variables, such as pollutant loads, measured at different points in space are related. More important from the practical viewpoint, such relationships can be used for estimating values at sites where no measurements are taken. A sample of data may consist of observations taken in one, two or three dimensions. Measurements of water quality along a river are taken in one-dimension. On the other hand, rainfall or other meteorological variables are measured at par- 1 2CHAPTER 1. INTRODUCTION ticular points, but collectively they constitute a two-dimensional random field. Finally, measurements of a specific mineral in ground is sometimes treated as a three-dimensional problem because it can occur over different levels. The spatial process is a stochastic process. It may be represented as a set of random variables (or vectors) Z(x), indexed by xwhich belongs to a set D⊂IRd, ad-dimensional euclidean space with d= 1,2,3. Its usual notation is {Z(x) : x∈ D}. Let us suppose one has the spatial locations x1, ..., xn, then Z(x1), ..., Z(xn) identifies observed data at those locations. These observations may be obtained from one or more, discrete or continuous, variables. According to Cressie (1993), the nature of Dallows us to identify three major spatial processes, namely lattice processes, point processes and continuous processes. The latter are commonly referred to as geostatistics (Matheron 1963) and, opposite to lattice processes, between any two spatial points associated to existent observations, there is always another point where the random variable could also be observed. The developed work falls within the scope of geostatistics and it applies theory of point processes, presuming that spatial locations have been produced by some form of stochastic mechanism. Under assessing spatial dependency, we consider two distinct issues, namely the estimation of the spatial dependency structure and the subsequent spatial prediction procedure. Within the realm of geostatistics, it is commonly believed that sample locations are equally spread over the observed region. Furthermore, it is assumed that the point process for data locations does not depend on the data process. So, under non-standard sampling, we are considering the failure of one or both of previous assumptions. In Chapter 2, we introduce notation and definitions related to geostatistical data modelling. We review some well-known facts about the convenient assumption of stationarity of the underlying process. The usage of the variogram as a tool to measure spatial dependence between samples is highlighted. In Chapter 3, we start focusing the estimation of spatial dependency under 3 standard sampling. A bibliographic search of current variogram estimators is described. A comparison simulation study is presented, covering different kinds of spatial dependence situations. In Chapter 4, a motivating example is introduced, the Rongelap island, where the data was collected over a two-stage process of uniform and clustered samples, which may have an impact on conclusions. Centered on the multi-stage case, we assess the effect of clustered and biased sampling on spatial dependency estimation. A new variogram estimator for clustered data is proposed. In Chapter 5, we proceed with the theoretical study of the proposed estimator. It is shown to enjoy good properties, such as asymptotic unbiasedness and consistency. In Chapter 6, we introduce the preferential sampling concept, as a formal definition for the dependency of data locations on data values. A flexible class of point processes for preferential sampling, based on log-Gaussian Cox processes, is presented. We then assess the effect of preferential sampling on the classical geostatistical methods typically used for estimation of the spatial dependency and for spatial prediction (kriging). In Chapter 7, a motivating example of pollution data is introduced, reinforcing the importance of doing parametric model analysis when data is suspected to be preferentially sampled. An intuitive approximate model for preferential sampling is proposed. We proceed with likelihood inference to estimate the parameters of this model and a simulation study is performed to show the benefits of this modelbased approach versus the traditional one. We close this Chapter with a short discussion of future work. 4CHAPTER 1. INTRODUCTION Chapter 2 Stationary spatial processes 2.1 Introduction Suppose that {Z(x) : x∈D⊂IRd}is a random spatial process, where Dis a bounded region with positive d-dimensional volume. The mean of Z(x), E[Z(x)], is the mean of all possible realizations of the process at points x. Typically, just one realization of the given process is observed, possibly denoted by z(x). The process is said to act over a random field Ω. Additionally, the difference random process Z(x)−E[Z(x)] represents departures of the original process from the mean at the points considered. The study of such processes is based on the identification of appropriate characteristics of regularity, which is referred to as stationarity in the context of stochastic processes (e.g. Kottegoda and Rosso 1997). The process Z(x) is usually defined through the finite-dimensional distribution Fx1,...,xn(z1, ..., zn) = P{Z(x1)≤z1, ..., Z(xn)≤zn}, n ≥1, which must satisfy the Kolmogorov’s conditions of symmetry, i.e. remain invariant when zjand xjare subject to the same permutation, and the consistency condition, i.e. Fx1,...,xn+m(z1, ..., zn,+∞, ..., +∞) must be equal to Fx1,...,xn(z1, ..., zn). Furthermore, the following hypothesis for the values of the mean and the variance may be considered. 5 6CHAPTER 2. STATIONARY SPATIAL PROCESSES Strict (or strong) stationarity A spatial process {Z(x) : x∈D}is strictly stationary, if for any finite collection {x1, ..., xn}of spatial locations, for any finite collection of real values {z1, ..., zn} and for any vector u∈IRdfor which xi+u∈Dwhenever xi∈D, then: •Fx1,...,xn(z1, ..., zn) = Fx1+u,...,xn+u(z1, ..., zn),∀n, u This means that a stationary process remains invariant when subject to translation transformations of its coordinates. Second-order (or weak or wide-sense) stationarity The spatial process {Z(x) : x∈D}is second-order stationary, if its first moment is a constant and the covariance between two variables is a function of the difference between theirs locations: •E[Z(x)] = µ(x) = µ∀x∈D •Cov[Z(xi), Z(xj)] = c(xi−xj),∀xi,xj∈D The function c(.) is called the stationary covariance function or sometimes covariogram. The function µ(x) is known as the trend of the process. Intrinsic stationarity The spatial process {Z(x) : x∈D}is intrinsically stationary if its first moment is a constant and the variance of the difference between two variables is a function of the difference between theirs locations: •E[Z(x)] = µ(x) = µ∀x∈D •Var[Z(xi)−Z(xj)] = 2γ(xi−xj),∀xi,xj∈D The function 2γ(.) is called the variogram and γ(.) the semivariogram, but the latter is often also referred as to the variogram to simplify terminology. We shall focus on the importance of this function in geostatistics. 2.1. INTRODUCTION 7 One may note that if a process is strictly stationary, then it is also secondorder stationary. Furthermore, if a process is second-order stationary, then it is also intrinsically stationary. Let us confirm this last implication: Var[Z(xi)−Z(xj)] = Var[Z(xi)] + Var[Z(xj)] −2Cov[Z(xi), Z(xj)] = =c(0) + c(0) −2c(xi−xj) = 2γ(xi−xj), being γ(.) = c(0) −c(.) (2.1) Strict and second-order stationarity coincide if the spatial process is Gaussian, i.e. if the joint distribution of any finite collection of variables is Gaussian. Secondorder and intrinsic stationarity coincide if the variance is finite and it does not depend on x, i.e. Var[Z(x)] = c(0) = σ2<∞,∀x. Isotropy The random process {Z(x) : x∈D}is isotropic if it remains invariant when subject to rotations of coordinates, in contrast to the anisotropic process. For example, the intrinsic random process {Z(x) : x∈D}is isotropic if ∀xi,xj∈D then: •E[Z(xi)−Z(xj)] = 0, •Var[Z(xi)−Z(xj)] = 2γ(kxi−xjk) where k.kdenotes the Euclidean norm. The variogram here only depends on the distance between the two locations and not on the direction of the difference vector. The isotropy condition is not so restrictive in practice, since a linear transformation of the coordinates sometimes produces an acceptable approximation to isotropy, or it might be possible to fit a different variogram in different directions in the case of anisotropy. Physically, the former corresponds to a rotation and stretching of the original spatial locations. Algebraically, it means to apply some linear transformation to the space of locations, given some anisotropy angle and anisotropy ratio. 8CHAPTER 2. STATIONARY SPATIAL PROCESSES Ergodicity A subset of the second-order stationary random processes possesses an important property known as ergodicity, which is required for the estimation of the characteristics of a process based on its realizations. This is applicable if the estimates of its moments, taken from the available realizations, converge in probability to the theoretical moments, when the available sample increases. Hence, under ergodicity, one realization will suffice for the estimation of these moments. In practice, this property is normally assumed to hold. 2.2 The mean function, covariance function and variogram A clear difference between the covariance function and the variogram is that the former is a direct function of the association between two variables, whereas the latter measures the disassociation. Variograms are more general than covariance functions, and many important properties have been initially established for covariance functions. Gneiting, Sasv´ari and Schlather (2001) explores the relationship between the two functions and present some analogous results for variograms. For a second-order stationary and isotropic random process, Var[Z(x)] = σ2 and it is useful to write the covariance function as c(u) = σ2ρ(u), where ρ(.) is the correlation function depending on a scalar argument u. As example of an useful correlation functions adopted in geostatistical data modelling, we have the Mat´ern family (Wackernagel 1998) represented in Figure 2.1 and with algebraic form given by ρ(u) = {2κ−1Γ(κ)}−1(u/φ)κKκ(u/φ),(2.2) where κ > 0 and φ > 0 are parameters, and Kκ(.) denotes a Bessel function of order κ. The parameter φdetermines the rate at which the correlation decays to zero with increasing u. The parameter κdetermines the analytic smoothness of Z(.) (see e.g. Ribeiro Jr 2002 for more details). 2.2. MEAN, COVARIANCE AND VARIOGRAM 9 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 distance ρ(u) κ = 0.5 κ = 1.5 κ = 2.5 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 distance ρ(u) φ = 0.05 φ = 0.1 φ = 0.2 Figure 2.1: Examples of the Mat´ern correlation function: φ= 0.15, varying κ(left plot); κ= 1.5, varying φ(right plot). For a process with stationary covariance structure, equation (2.1) shows that the variogram reduces to γ(u) = σ2(1 −ρ(u)).(2.3) Typically the variogram approaches a constant value as the separation distance u increases; this value is known as the sill. In (2.3) the sill is given by σ2. When a variogram has a sill, it means that there is a distance beyond which the correlation between variables is zero; this distance is called the range or radius of influence. For many practical applications, it is useful to consider a Gaussian spatial process with a possibly varying mean function but stationary covariance structure. For such processes, Z(x)−µ(x) is a stationary Gaussian process with zero mean. A possible solution is to specify µ(x) as a regression model, with the aim of fitting a smooth surface to values measured over a sample of points. The regression itself provides a summary of the trend as well as a means of predicting the value at any location within the modelled surface. For example, the analysis of rainfall data in 16 CHAPTER 2. STATIONARY SPATIAL PROCESSES 2.5 Trend and outlier identification A sample variogram can be highly misleading if derived from data with a spatially varying trend. The same happens in the presence of outliers; bear in mind that a single outlying measurement, say Z(xk), might contribute to the variogram estimation with n−1 square differences (Z(xk)−Z(xi))2, i = 1, ...n and i6=k. Hence, it is important to consider some techniques of exploratory data analysis to identify possible non stationarity in the mean or isolated outliers. Returning to the example of rainfall data in Paran´a State in Brazil, we first show the rainfall data against each of the coordinates (two top panels in Figure 2.4). These confirm the trend from east to west and from south to north. We then show similar plots but replace the original data by residuals from a linear trend fitted by ordinary least squares (two bottom panels in Figure 2.4). Assuming that the remaining structure seen in the data can be attributed to the random part of the model, these ordinary least square residuals Z(xi)−bµ(xi) can be used for the subsequent variogram estimations. Moreover, this sample variogram could now be used to make a new estimation of the trend, as more reliable estimations of the deterministic and stochastic components of the spatial process depend on each other. Cressie (1993) summarizes useful methods for exploratory spatial data analysis. If data is located on a regular grid1, one possible solution is to calculate the sample mean or sample median (as a more robust estimator) across rows and down columns. Plots may then be created, summarizing row and column results. This allows us to identify the existence of a linear trend for mean or median along rows or columns. This type of analysis indicates whether and how the spatial location has influence on the variable values. The use of median and mean serves another purpose. The comparison of these two variables has the additional function of highlighting rows or columns that may 1When the data is not located on a regular grid, a low-resolution grouping of observations into a two-way table, still allows these methods to be carried out. 2.5. TREND AND OUTLIER IDENTIFICATION 17 200 400 600 200 300 400 west−east data 100 200 300 400 200 300 400 south−north data 200 400 600 −50 0 50 west−east residuals 100 200 300 400 −50 0 50 south−north residuals Figure 2.4: Paran´a rainfall data shown against each of the coordinates: top two panels display original data and bottom two panels display residuals after fitting linear trend. contain atypical observations. A high value for mean−median indicates a possible outlier. Another solution for stationarity analysis is drawing a bivariate plot of Z(x) and Z(x+u), for a fixed direction u, as xvaries over the data locations. In case of local stationarity, when norm of uis small enough, this plot points should be near the bisectrix. This is also a way of detecting outliers, which correspond to those isolated values found far from the bisectrix. The previous methods prove to be useful in detecting gross trends or isolated 18 CHAPTER 2. STATIONARY SPATIAL PROCESSES outliers. Furthermore, a technique like the one called pocket-plot can allow the identification of localized areas as atypical with respect to a stationary model, i.e. pockets of non stationarity. Once more, this is done by exploiting the spatial nature of the data along rows and columns. These pockets, once discovered, should be removed from variogram estimation; they may be considered afterwards, when modelled and incorporated into final results analysis. Non-stationary spatial processes So far, it has been assumed that the data has come from a stationary process, apart from small pockets of non stationarity. This is a convenient working assumption which can be relaxed in various ways. Sampson and Guttorp (1992) proposes spatial deformation models of the form D(x1,x2) = γ(f(x1), f(x2)) with γ(.) an isotropic variogram and f(.) a smooth nonlinear map from IRdto IRd0. In principle one may permit d06=dthough in most of their work the equality is assumed. The idea is that the map f(.) takes the coordinates from the geographical space into an alternative dispersion space in which stationarity holds. The transformation of the data itself can also help the non stationarity issue. When the responses Z(xi)i= 1, ..., n are continuous but the Gaussian model is clearly inappropriate, some additional flexibility is achieved by introducing an extra parameter λdefining a Box-Cox transformation of the response (see e.g. Ribeiro Jr 2002). Recently a number of methods, including kernel convolutions, deformations and spatially adaptive spectra, have been suggested to allow for non stationarity in the stochastic component of the underlying process. These methods build non stationarity directly into the covariance function (see e.g. Pintore and Holmes 2004 for a brief review). In addition, these authors show how, by working in the spectral domain, one can build non stationary covariance functions which are centred on 2.5. TREND AND OUTLIER IDENTIFICATION 19 popular classes of stationary models such as the Mat´ern or Gaussian. The resulting non stationary models are defined as “localised” versions of their stationary counterparts. Alternative proposals are found in references like Higdon, Swall and Kern (1999), Fuentes (2002) and Stein (2005). Bearing in mind the existence of methods to tackle non stationarity, such as previous ones, and the advantages of increased simplicity coming from more restrictive modelling assumptions, we shall keep the assumption of stationarity of the underlying spatial process in the remaining Chapters. A more restrictive assumption can, indeed, help models remain easily interpretable and, in case of complex models fitted to sparse data, can help to avoid issues of poor identifiability of model parameters. 20 CHAPTER 2. STATIONARY SPATIAL PROCESSES Chapter 3 A comparison of approaches for valid variogram achievement 3.1 Introduction In spatial prediction, a basic question is that, given a set of nobservations of the process Z(.) at points xi,i= 1, ..., n, what is the value taken by the variable at a point, x0, where data are unavailable? The approach differs from regression in that local features can affect the solution. In principle, all measurements should be considered. Having in mind that some measurements in the vicinity of the point investigated, or sometimes elsewhere, are more closely related than others to the true value at point x0, the appropriate procedure would be to adopt a weighted mean: b Z(x0) = Pn i=1 λiZ(xi). This linear combination may be considered to be an optimum estimate if the coefficients, or weights, λiare such that they sum to one and the estimator is unbiased and has minimum variance1. Only data within the radius of influence should be considered. 1More precisely, this estimator is classified as BLUE, i.e. best linear unbiased estimator. 21 22 CHAPTER 3. COMPARISON OF VALID VARIOGRAMS The estimation of a valid variogram plays here a decisive role, as it is commonly used to find the optimal solution to the values of the weights. The method is called kriging, coined this way by Matheron (1963) to honour the mining engineer D.G.Krige. The optimum solution can be found by using the Lagrange multiplier (see e.g. Kottegoda and Rosso 1997). In an isotropic field with estimated variogram values bγ(uij) between points xiand xjat distances uij, the estimated weights, b λj, j = 1, ..., n may be found by solving the following n+1 simultaneous equations:          Pn j=1 λibγ(uij) + λ=bγ(ui0), i = 1, ..., n Pn i=1 λi= 1 where index “0” relates to the unsampled location and λis the Lagrange multiplier. Once the weights are estimated, a prediction value can be easily obtained for the process Z(.) at the point with coordinates x0. Some more details about kriging methods will be discussed in Chapter 6. In this Chapter our main aim is to stress the contribution of the variogram with respect to inference procedures. Moreover, variogram analysis provides a useful tool for summarizing spatial data and it may be used to measure spatial dependence between samples. Commonly used variogram estimators The first proposal, for the presence of stationary processes, for a variogram estimator is due to Matheron (1962). This estimator is based on the method of moments and it is often referred to as the classical estimator: bγ(u) = 1 2|N(u)|X N(u) (Z(xi)−Z(xj))2(3.1) where N(u) = {(xi,xj) : kxi−xjk=u, u ∈IR}and |N(u)|is the total of pairs in N(u). Matheron’s estimator is unbiased2, however it presents some drawbacks 2It is unbiased for γ(.) when Z(.) is intrinsically stationary (Cressie 1993, page 71). 3.2. TRADITIONAL THREE STAGES 23 such as being badly affected by atypical values due to the squared term in the summand of (3.1). In general, its statistical properties are difficult to study. If Z(.) is a Gaussian random process, then bγ(u) is a linear combination of χ2random variables on one freedom degree. According to Journel and Huijbregts (1978), a minimum of 30 pairs is recommended in |N(u)|. When data is not regularly spaced, this estimator can be obtained by considering a tolerance region around u. Note that the squared term is also used to propose a related variogram, called the variogram cloud. If {(Z(xi),xi) : i= 1, ..., n}is the sample data set, then the scatterplot of the points {(uij, vij) : j > i, uij =kxi−xjk, vij = 1 2(Z(xi)−Z(xj))2}identifies the corresponding variogram cloud. As expected, this estimator is also sensitive to outliers. Cressie and Hawkins (1980) has minimized this weakness, by working with square-root absolute differences and, under a Gaussianity assumption, have produced the estimator bγ(u) = n1 2|N(u)|PN(u)|Z(xi)−Z(xj)|1 2o4 0.457 + 0.494 |N(u)| (3.2) where the term 0.457 + 0.494 |N(u)|is used to make it unbiased. Unfortunately, it has been suggested that the estimators in (3.1) and (3.2) should not be used for inference and prediction. The reason for this is that they may fail the conditionally negative-definite property which may lead to absurd negative values for the mean square prediction errors, as proved in Cressie (1993). If this were to occur, then the estimators are deemed to be invalid. 3.2 Traditional three stages A common approach to achieving a valid variogram estimator is to approximate an empirical variogram by some theoretical model which is known to be valid. The 24 CHAPTER 3. COMPARISON OF VALID VARIOGRAMS idea is to select, within the families of valid variograms, a function which captures the underlying spatial dependence of the available data. Traditionally, these type of approaches are accomplished through three distinct stages: 1. Compute an empirical variogram (typically non valid); 2. Choose a theoretical model among the family of valid parametric or nonparametric variograms; 3. Estimate the variogram by fitting the theoretical model to the empirical variogram. To accomplish these tasks, there are several different approaches strongly defended by their authors. In this Chapter, our work’s main purpose was to identify these approaches (Menezes 2002) and compare some of them based on a numerical study, covering different kind of spatial dependence situations. The comparisons are mainly based on the integrated squared errors of the resulting valid estimators. The main contributions of this work appear in Menezes, Garcia-Soid´an and Febrero-Bande (2005a). Note that some authors prefer to group these three stages into two parts, variogram estimation and variogram fitting; the latter part incorporates stages 2 and 3 simultaneously (e.g. Cressie 1993). In contrast, we argue that, when possible, three separate stages allow a better classification of the existing approaches. The output of stage 2 is a vague valid candidate and its complete specification is only obtained from stage 3. Before giving details about the complete approaches that we examined, we make some generic comments on each of the previously listed stages. We shall point out some references, if we think they introduce a relevant idea for the implementation of these tasks. 3.2. TRADITIONAL THREE STAGES 25 3.2.1 Stage 1 – Empirical variogram estimation The word “empirical” means based on observation or experiment. The estimation of the empirical variogram always, unsurprisingly, begins with the observed data, whichever estimator is used. Examples include those estimators introduced in Section 3.1. Robustness to outliers is normally considered an important characteristic for any estimator. In this regard, some other robust empirical estimators have been proposed in addition to (3.2) (see e.g. Dutter 1996, for a general review, or Gunst and Hartfield 1997, for large data sets). They usually avoid one or both of the following issues associated with the classical estimator: •the square term, (Z(xi)−Z(xj))2, because it induces distortion in data values; •the mean, because it is not a robust location estimator. For instance, Armstrong and Delfiner (1980) proposed to use the square of the interquartile range of the differences [UQ{Z(xi)−Z(xj)}−LQ{Z(xi)−Z(xj)}]2 where UQ and LQ stands for upper and lower quartiles, i.e. the 75th and 25th percentiles of the differences Z(xi)−Z(xj). Additionally, they have also considered a sample quantile, the median, of squared differences: med{(Z(xi)−Z(xj))2: (xi,xj)∈N(u)}. Note that these approaches need a correct normalization to make them unbiased. Genton (1998a) proposes a variogram estimator based on the highly robust scale estimator of Rousseeuw and Croux (1992,1993), denoted below by QNu. Considering Z(.) a intrinsically stationary but not necessarily isotropic process, 32 CHAPTER 3. COMPARISON OF VALID VARIOGRAMS Genton (1998b) refuses to accept WLS as the solution for GLS complexity and proposes an explicit formula for the covariance structure V, calling the resulting method GLSE. The basic idea is to obtain a generic covariance structure by using, iteratively, the correlation structure of Matheron’s estimator in the independent case. The main steps of this algorithm for variogram model fitting are: 1. Determine the matrix V=V(θ) such that Vij =Cov(2bγ(ui),2bγ(uj)) = Corr(2bγ(ui),2bγ(uj))γ(ui;θ)γ(uj;θ) pN(ui)N(uj) where Corr(2bγ(ui),2bγ(uj)) can be approximated by the result obtained for the Matheron’s estimator in the independent case; 2. Choose θ(0) randomly or by using OLS or WLS criteria, and let i= 0; 3. Compute matrix V(θ(i)) and determine θ(i+1) which minimizes G(θ) = (2bγ−2γθ)TV(θ(i))−1(2bγ−2γθ) ; 4. Repeat (3) until convergence to obtain b θ. Genton concludes his work by carrying out some simulations to show that the GLSE criterion, combined with a robust variogram estimator, may improve the fit significantly, even in the presence of outliers. 3.2.4 Existing combinations of the previous stages Next, we shall introduce some existing complete approaches to reach our target: a valid variogram estimator. All of them result from distinct combinations of previous stages, summarized in Table 3.1. We shall begin with a mandatory reference, Zimmerman and Zimmerman (1991), where seven different approaches are compared through a Monte Carlo simulation study. This comparative study, in spite of being considerably exhaustive, 3.2. TRADITIONAL THREE STAGES 33 is somehow restricted in scope as it only involves parametric techniques. In fact, these seven approaches are mainly distinguishable by their third stage, four of which are LS-based, two ML-based and one using a modified MIVQU. Their main conclusions may be summarized as follows. They confirm that the performance of each estimator improves, in the sense that its distribution becomes less dispersed as the sample size increases. In general terms, the estimators of parameters perform better as the spatial dependence is weaker. However, the standard 95% prediction intervals perform better when the spatial dependence is stronger. Moreover, they conclude that, in some particular situations, the likelihood-based methods can perform a little better than the least squares methods. These are, however, less computationally demanding and are deemed to be perfectly acceptable in most cases. The paper of Shapiro and Botha (1991) is pioneered towards selecting a valid model in a non-parametric space, as we mentioned in Section 3.2.2. They combine the Matheron’s estimator at stage 1, a broad class of permissible variograms at stage 2 and at the last stage, a WLS fitting criterion where the optimization problem is reduced to a quadratic programming problem. Following Christakos (1984), they define f(u), u ∈IRdas a permissible variogram function, if it is continuous (except possibly at the origin), f(u) = f(−u), f(u)≥0 for all u, and −f(u) is conditionally nonnegative definite. The resulting valid variogram estimator then fulfills equation (3.4). Additionally, requirements such as variogram smoothness, monotonicity or convexity may be incorporated into the fitted variogram γvleading to a better approximation. Suppose the empirical bγ(ui) are scattered, then γvmay change rapidly. In this case, it may be important to impose a smoothness condition, forcing this function’s derivative to be bounded. It may also be important that the γv be monotonically increasing or that it be convex. As expected, the monotonicity condition may be ensured by forcing the derivative to be positive for all u > 0, 34 CHAPTER 3. COMPARISON OF VALID VARIOGRAMS and the convexity condition by forcing the second derivative to be negative. This approach was evaluated by Cherry, Banfield and Quimby (1996), where they conclude that this “non-parametric method is faster, easier to use and more objective than parametric methods”. Gribov, Krivoruchko and Ver Hoef (2000) suggests a new method of computing the empirical variogram of Matheron. The squared differences [Z(xi)− Z(xj)]2are binned into Kdistinct bins, and point estimations of the semivariogram at Kpoints are obtained. The complicating issue is how to best bin the data. They introduce the notion of logarithmic increases in the size of tolerance regions against the traditional fixed size. This new concept allows better results in estimation near the origin. They also propose to use a kernel method to assign, within a given bin, weighted values depending on how close a value is to the center of the bin. This requires fewer elements per bin than the recommended minimum of 30 pairs suggested from the classic guideline of Journel and Huijbregts (1978), as well as the weights’ presence minimizes a possibly existing unequal distribution of lags. With respect to the model fitting stage, they propose a modified WLS3proce- dure. They split this stage into two steps. At step 1, typically with two iterations, they consider logarithmic lag sizes. At step 2, a default lag size obtained from the range estimate in step 1 is used instead. The last approach included in our survey is the one proposed by Garcia- Soid´an et al. (2004). These authors propose the usage of the non-parametric empirical estimator given by equation (3.3), together with the permissible function of Shapiro and Botha (1991). The empirical bγh(u) and the theoretical curve are fitted through a re-iterated WLS criterion. The former is shown to have desirable properties, such as asymptotically unbiasedness and consistency. 3This algorithm is included into the Geostatistical Analyst extension to GIS ArcInfo/ArcView8.1 (Krivoruchko, 1999). 3.2. TRADITIONAL THREE STAGES 35 One should have in mind that an important issue of kernel estimation is the selection of the bandwidth parameter, h. These authors address the problem by asymptotically minimizing the mean square error (MSE) or the mean integrated square error (MISE), in order to derive the local and the global bandwidth, respectively. Both expressions involve the unknown function γ(u). For the purpose of the bandwidth derivation, a simple parametric approach, like the first one presented by Zimmerman and Zimmerman (1991) (see Table 3.1), may be used to estimate γ(u). This isolated parametric estimation can even be improved by being incorporated into an iterated non-parametric procedure. Table 3.1: Taxonomy of existing approaches for valid bγachievement. Bold identifies those approaches selected for the comparative study and, between brackets, main equations references are found. Approaches Stage 1 Stage 2 Stage 3 Zimmerman Matheron(3.1) P model OLS and Cressie-Haw.(3.2) P model WLS Zimmerman Matheron P model WLS (1991) Matheron P model WLS-Delfiner (1976) — P model ML —P model REML Matheron P model OLS+MIVQU Shapiro and Matheron NP function(3.4) WLS Botha (1991) Gribov Matheron-modif. P model WLS-modified et al. (2000) Garcia-Soid´an NW kernel(3.3) NP function(3.4) WLS et al. (2004) 36 CHAPTER 3. COMPARISON OF VALID VARIOGRAMS Outside the boundary, the bias of the Nadaraya-Watson estimator (3.3) is of the order h2; however, the latter order amounts to hfor distances uclose to 0. Then, proceeding as in Kyung-Joon and Shucany (1998), we may denote by ˆγq,h(u) the estimator obtained by substituting a boundary kernel Hqfor the symmetric one Kin (3.3), where q= min {uh−1, C}and Hq(z) = K(z)−rL(z) 1−r, z ∈[−C, q] where Kand Lare symmetric kernel functions, r=c1,K c0,L(c0,Kc1,L)−16= 1 and ci,G =Rq −CziG(z)dz. This particular selection of the boundary kernel Hqproduces a variogram estimator ˆγq,h(u) that makes it negligible the term of order hin the bias and preserves the same convergence orders for all u > 0, as shown in Garcia- Soid´an et al. (2004). 3.3 Simulation study In order to analyze the performance of the previous approaches for valid bγachievement, simulations of spatial data in IR2were carried out for different kinds of dependence situations. We considered the spherical and the exponential variogram models given in (2.5) and (2.6), respectively. Additionally, the wave model was also considered, because of its atypical irregular behaviour: •Wave model: γw(u;θ) = θ0+θ1[1 −θ2sin(u/θ2)/u], u 6= 0. Bear in mind that we want to restrict ourselves, in this Chapter, to the estimation of the spatial dependency under standard sampling. Thus, in all cases, a uniform distribution on [0,1] ×[0,1] was assumed for spatial locations xi= (xi1, xi2), i = 1, .., n, where nrepresents the sample size. Several data sets were generated with Gaussian data, Z(xi), i = 1, .., n, using one of the above variogram models. The parameters of these models were chosen in such a way that the corresponding curves were comparable according to their range. We fix the values for the nugget θ0and θ1to be 0.25 and 5.0, respectively. The third 3.3. SIMULATION STUDY 37 parameter was the one chosen depending on the model: exponential, θ2= 0.167; spherical, θ2= 0.5; and wave, θ2= 0.113. With this selection, the theoretical variograms have a sill of 5.25 and a range (referred to as the minimum value for which the variogram reaches either the sill or 95% of the sill, in case that the range is not finite) of 0.5. More precisely, the wave model oscillates around the sill value and, consequently, the 0.5 value identifies the global maximum of the corresponding variogram function. 3.3.1 Comparing empirical estimators The aim of our first numerical study was to compare the three main empirical estimators used at stage 1 for the approaches included in Table 3.1; these are given in expressions (3.1), (3.2) and (3.3). For data generation, we took a sample size of n= 200 and we started by selecting the exponential model. Unusual estimated values were obtained using estimators (3.1) and (3.2) for the largest lags. Additionally, some of them did not have the recommended minimum of 30 pairs. Therefore, in posterior simulations, we have decided to only consider the first 55% of lags. One may note that this guideline is still less conservative than the one proposed by Journel and Huijbregts (1978), who specifies that the largest used lag ukshould be less than or equal to half of the largest existent lag. As non-parametric estimation requires more lags than those empirically obtained, we have also decided to consider a larger number of lags, equally spaced, within interval [ min(uk),0.55 ∗max(uk) ]. Following these considerations, Figure 3.1 shows the obtained data, as well as two more graphs assuming the spherical and the wave models for data generation. All graphics included in this Chapter use the following notation: lines are used to represent a valid estimator; and isolated symbols, e.g. small squares, are used for empirical estimates. Figure 3.1 demonstrates the behaviour of the estimator when one sample is considered, although it will depend strongly on the sample variability. For this 38 CHAPTER 3. COMPARISON OF VALID VARIOGRAMS a) Exponential b) Spherical c) Wave Figure 3.1: Three empirical estimators and the associated theoretical curve. Data simulated with three distinct models. Sample size equals 200. 3.3. SIMULATION STUDY 39 01234 Matheron Cressie−Haw. NW kernel 01234 01234 a) Exponential b) Spherical c) Wave Figure 3.2: Boxplot of the evaluated ISE from three empirical estimators, using data simulated from three distinct models. The simulation consisted of 100 replications, each with a sample size of 200. reason, we include a second study where 100 independent samples are considered. For each one, the integrated squared error (ISE) between each of the three empirical estimators and the theoretical variogram, given by ISE =Z[bγ(u)−γ(u)]2du, (3.6) was approximated numerically through the trapezoid rule, bγ(u) represents an empirical estimator and γ(u) represents the theoretical curve. This simulation was repeated for the previous models: exponential, spherical and wave. The results are summarized in the boxplot in Figure 3.2. If one compares the median values associated with the three estimators, then the best performance is clearly achieved by the non-parametric estimator, using the Nadaraya-Watson kernel. Another advantage of this non-parametric estimator is that it is a continuous function. In contrast, estimators (3.1) and (3.2) propose point values of the semivariogram for given distances u, making them discontinuous. Most analyses requires knowledge about estimations in a continuous range of γ(u). 40 CHAPTER 3. COMPARISON OF VALID VARIOGRAMS We conclude by bringing attention to the different orders of magnitude of the ISE values for each theoretical model, lower values are associated to the exponential curve and the largest ones to the wave curve. 3.3.2 Comparing complete approaches We highlight three approaches (marked in bold) from Table 3.1, which we consider the most representative of the existing alternatives. For two of the approaches, a valid model is chosen within the space of parametric families. Hereafter, they are identified as the parametric approaches (P), one of which uses WLS as the fitting criterion and the other REML. The third approach, introduced by Garcia-Soid´an et al. (2004), will be referred to as the non-parametric approach (NP). The superior results of the Nadaraya-Watson kernel estimator, when compared to the Matheron’s estimator, led us not to include Shapiro and Botha (1991) in our numerical study. Gribov et al. (2000) was also excluded as, under isotropy, their main contribution is reduced to the usage of weights within a given bin. In this case, the kernel estimator does not differ much from their proposal and may be indeed a better choice. Under the NP approach, we preferred to asymptotically minimize the MSE to derive a local bandwidth parameter. For this purpose, the symmetric Epanechnikov kernel was employed. Additionally, as the bandwidth derivation needs itself an estimation of the variogram, the available WLS parametric estimation was used for this purpose. Near the variogram endpoint 0, a specific asymmetric boundary kernel was constructed from the Epanechnikov kernel and the quartic kernel. For the implementation of the REML fitting criterion, we used the geoR library from R, which provides several functions for geostatistical analysis as explained in Ribeiro Jr and Diggle (2001). Excluding this particular case, we used Fortran to implement our numerical study. Figure 3.3 shows an example of results obtained with the three selected approaches, when using each of our theoretical models for data simulation and a 3.3. SIMULATION STUDY 41 a) Data simulated with an exponential model b) Data simulated with a spherical model c) Data simulated with a wave model Figure 3.3: Approaches to achieving a valid bγ(u): the 2 parametric approaches are on the left and the non-parametric approach is on the right. Data simulated with three distinct models. Sample size equals 50. 48 CHAPTER 3. COMPARISON OF VALID VARIOGRAMS over-estimate the sill approximation when sample size is equal to 200. An aspect also worth mentioning is that, overall, the performance of each estimator improves as sample size increases. Computational costs A final remark about the numerical study is related to the computational cost of the three approaches chosen to achieve a valid estimator. The CPU execution times were recorded for each sample, without considering data simulation but just the time needed to implement all existing stages. The results are summarized next in Table 3.4. n=50 n=200 Pwls 1.3 s 1.5 s Preml 3.0 s ≈30 s NP ≈30 s ≈30 s Table 3.4: Summary of the computational costs of the three approaches chosen for our comparison study. The lowest computational cost was achieved by the Pwls approach, being around 1.3 and 1.5 seconds for n= 50 and n= 200, respectively. The cost for the Preml approach was around 3 seconds for n= 50, this is at least 10 to 15 times greater for n= 200. With respect to the NP approach, we have registered CPU times from 27 to 36 seconds for n= 50, being the lowest values associated to the spherical data and the greatest to the exponential data. These computational costs have only shown a slight increase when we moved to sample sizes of n= 200. Bear in mind that the heavy costs obtained for the NP approach are usually justified by the optimal bandwidth derivation. 3.3. SIMULATION STUDY 49 3.3.3 Closing remarks The problem of estimation of the variogram can be analyzed in practice from several points of view. If the aim is just to obtain an approximation of the dependence structure of the spatial data, then the classical and the Nadaraya-Watson kernel provide good estimators that behave better than the robust estimator proposed by Cressie and Hawkins, using as a term of comparison the values estimated for the median and interquartile range of the ISE; however, the robust estimator reduces the range of variation of the ISE. If we focus on the problem of spatial prediction, we modify the variogram estimators to obtain valid variograms; otherwise, negative mean squared prediction errors may be achieved. From the different alternatives discussed, the valid kernel estimation (referred to as the NP approach) has the best performance for large sample sizes in terms of the values estimated for the ISE, regardless of the parametric model that is considered. In this respect, it is surprising that fitting the correct parametric family does not produce a better fit than the non-parametric method. The misspecification of the parametric family has a second order effect on the kernel estimator, since it affects estimation of values associated to the bandwidth parameter. On the other hand, when considering typical features associated with the variogram (nugget, sill and range), we conclude that the valid kernel estimation provides the lowest MSE values, although the P approaches prove competitive in the estimation of the corresponding median values. In general, the results presented here show that a valid variogram estimator obtained from a NP approach is a good alternative to those valid estimators obtained from the classic parametric approaches. The NP approach has the additional advantage of avoiding problems associated with using the wrong parametric model, which can occur in many conventional approaches. These advantages become even more evident if sample data underlies an atypical spatial dependence, like the one from the wave model. However, one must be prepared to pay an extra computational cost over the cost associated to a simple P approach like the one that fits a 50 CHAPTER 3. COMPARISON OF VALID VARIOGRAMS valid model to some empirical estimations through the WLS criterion. The P approach using REML as fitting criterion is only able to compete with the other methods in the presence of small datasets and, simultaneously, observed data does not follow a wave-type structure. Chapter 4 Assessing the effect of clustered and biased multi-stage sampling 4.1 Introduction The geostatistical methods introduced in Chapter 3 rely on the expected assumption that the sampling design for locations xi,i= 1, ..., n is deterministic or it is stochastic but independent of the data process, and all analyses are carried out conditionally on xi(Diggle, Ribeiro Jr and Christensen 2003). It is then assumed that the sampling points have been chosen independently of the values of the spatial variable. However, dependencies can occur due to the adopted sampling method, such as the favored selection of specific areas that are believed critical (e.g. maximum values search). Schlather, Ribeiro Jr and Diggle (2004) proposes methods to detect the dependence between marks and locations of marked point processes. As described in Mateu and Ribeiro Jr (1999), the random field (that we have been studying) and the marked point process are two type of spatial processes such that: •The former is defined in every point of the observed region, and the sample positions can be determined by the scientist himself (example of deterministic 51 52 CHAPTER 4. CLUSTERED AND BIASED MULTI-STAGE SAMPLING sampling design); •For the latter, the locations are always given by a stochastic point process, and interactions among the locations and the marks are normally expected. Otherwise, one has the so called random field model (marked point process becomes a special class of a random field). If the data are consistent with a random field model, the point pattern and the marks can be analysed separately using standard techniques for point processes (e.g. Ripley 1981 and Diggle 2003) and for geostatistical data (e.g. many references in Chapter 3). Therefore, this analysis is greatly simplified. The examination of second-order characteristics, like the variogram, of a spatial process should consider if data come from a random field or a genuine marked point process. Example of references concerned with this subject are Walder and Stoyan (1996), Mateu and Ribeiro Jr (1999) and Schlather (2002). Schlather et al. (2004) indicates next two likely situations for point and data processes being dependent, and subsequent failure of this important geostatistics assumption. Firstly, if the dependency is an intrinsic property of the data themselves, for example the relative positions of trees impact on their size due to their competition for light and nutrient. This is the case of genuine marked point processes. Alternatively, this dependency can be justified by a prior scientific knowledge of the spatial variable of interest, for example of the expected local level of contamination in air pollution. This can lead to the gathering of samples in areas with atypical values. Our work concerns the problems resulting from the second situation, that we think of major importance in geostatistics because of its high likelihood of occurrence on actual field measurements, and often either ignored or addressed by generic techniques like declustering ones (e.g. Goovaerts 1997 and Isaaks and Srivastava 1989). 4.2. ASSESSING THROUGH SIMULATION 53 In this Chapter, we are motivated by the application example of the radioactivity data from Rongelap island, where a two-stage data collection was used, leading to the presence of clustered data. So that we start by restricting our attention to multi-stage samples, aiming to assess the presence of multi-stage dependence, or also referred to as sequential dependence, where the choice of sampling points is driven by previous measurements. We propose some data exploratory methods which are intended to detect biased multi-stage collection of spatial data. We investigate corrector models that aim to minimize the impact on variogram estimation due to the adoption of the type of non-standard sampling designs just described. Moreover, we assess the effect of these methods on the Rongelap data. 4.2 Assessing through simulation We start by using simulated data to develop and study our diagnostic tools for data analysis. It is well known that simulation allows a level of knowledge and control that leads to more robust and defendable solutions. Using simulated data sets, where the characteristics of the data and the sampling designs are controlled and varied, will help the research of the technique’s potential, and to assess its performance in specific situations. We can gain insight about what happens when assumptions are violated since the true model is known. 4.2.1 Sample generation algorithm Typically, when one carries out some study of geostatistical data, the sample locations are uniformly spread over the observation region. Suppose now that one wishes to proceed with a multi-stage collection of data. If the goal is to better characterize the spatial variability for short distances, then one solution is to include some clusters of locations into later stages. Alternatively, suppose the goal is, as exemplified before, to pursue the maximum values of the spatial variable 54 CHAPTER 4. CLUSTERED AND BIASED MULTI-STAGE SAMPLING of interest, then the complete sample data set is expected to be mainly represented by large data values. All previous situations may condition the sampling design. In our simulation study, we shall then consider four distinct sampling designs: •complete spatial randomness (CSR); •just clustered; •biased but non-clustered; •and, finally, biased and clustered. Furthermore, we start with a two-stage approach for sampling collection, with the second stage potentially influenced by the first. We consider spatial locations x within the unit square [0,1] ×[0,1]. Data sets are generated with Gaussian data, Z(x), using some chosen variogram model. We prefer to propose a more generic algorithm where Kclusters can be generated, each one inside a sub-region Rk. For example, if one wishes to produce a biased sample with just one cluster, this can be done by restricting the sampling points from stage 2 to R1and around1the maximum of measurements from stage 1. The total sample size will be n, with n1from stage 1 and n2from stage 2. The algorithm may be summarized as follows 1. Sample n1points xiat random on [0,1] ×[0,1]; 2. Generate Z= (Z(x1), ..., Z(xn1)) ∼MV N; 3. For k= 1, ..., K (a) If biased=TRUE then select Z(xm,k) = max i{Z(xi)|xi∈Rk} else select Z(xm,k) = random i{Z(xi)|xi∈Rk}; (b) Sample n2,k points at random on [xm,k−δ , xm,k+δ]2; 1A small square with side length 2δwill be considered. Points must be simultaneously inside of the unit square. 4.2. ASSESSING THROUGH SIMULATION 55 4. Consider n2= K P k=1 n2,k; 5. Generate Z∗= (Z(xn1+1), ..., Z(xn1+n2)) where Z∗|Z∼MV N ¡ΣT 12Σ−1 11 Z,Σ22 −ΣT 12Σ−1 11 Σ12 ¢ and Σ22 = var{Z∗}, Σ11 = var{Z}, Σ12 = cov{Z,Z∗}; The conditional distribution from step 5 was derived from the joint distribution using properties of the multivariate Gaussian distribution (Anderson 1984). Under the adopted notation, the sampled points from stage 1 share a common time label t0, while those generated at stage 2 are assigned a time label t1. Bear in mind that a completely random sample can be obtained avoiding stage 2, i.e. n2= 0, or generating the n2points uniformly spread over all unit square. Moreover, the cluster effect tends to disappear for a large K. The two-stage approach reflects more directly the sampling design defined for Rongelap data. The previous algorithm can be easily extended to more than two stages, even though our experience confirmed that similar results are obtained. We also tried a specific multi-stage approach, hereafter termed serial sampling, according to which all points and corresponding data values after stage 1 are generated one at a time. So, we will have {(xi, Z(xi)) : i= 1, ..., n1}from stage 1, (xn1+1, Z(xn1+1)) from stage 2, (xn1+2, Z(xn1+2)) from stage 3, ..., and (xn, Z(xn)) from last stage. Figure 4.1 shows an example of spatial locations derived by our serial sampling algorithm. The latest points, associated to time labels tiwhere i > 0, are conditioned by the maximum of previous measurements. This algorithm also tends to originate a cluster of biased data, within a neighborhood of length approximately equal to 2δ. 56 CHAPTER 4. CLUSTERED AND BIASED MULTI-STAGE SAMPLING 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 X Y time t0 time ti where i>0 Figure 4.1: Serial sampling algorithm. The center point of each square identifies the maximum value at each stage. 4.2.2 Impact on variogram estimation We now want to analyse the impact of clustered and biased multi-stage sampling on variogram estimation. We consider two popular empirical estimators, the classical one from Matheron (1962) and the Nadaraya-Watson kernel estimator, given in (3.1) and (3.3), respectively. Figure 4.2 illustrates an example of degraded behaviour of these estimators under this type of non-standard sample. The resulting estimations are plotted against the theoretical curve of the variogram model chosen for sample generation. In case A, the data set was obtained by random sampling, whereas in case B a serial sampling was considered with 70% of biased clustering2. The results of case B are, at least partially, justified by the sample locations not being sufficiently 2In our serial sampling algorithm, we have specified a total n= 200 and n1= 60 for stage 1. 4.2. ASSESSING THROUGH SIMULATION 57 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 A− spatial locations x y 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 B− spatial locations x y 0.0 0.2 0.4 0.6 0.8 0 2 4 6 8 10 A− empirical estimators distance variance Theoretical Matheron NW kernel 0.0 0.2 0.4 0.6 0.8 0 2 4 6 8 10 B− empirical estimators distance variance Figure 4.2: Behaviour of bγunder random sampling (case A) against clustered and biased sequential sampling (case B). 64 CHAPTER 4. CLUSTERED AND BIASED MULTI-STAGE SAMPLING 0.0 0.2 0.4 0.6 0.8 1.0 1.2 −0.2 0.0 0.2 0.4 0.6 I− Random distance conditional expectation Eseq(u)−E(u) Eseq(u)−E(Z) E(u)−E(Z) CI CI 0.0 0.2 0.4 0.6 0.8 1.0 1.2 −0.2 0.0 0.2 0.4 0.6 II− Just Clustered distance conditional expectation Figure 4.4: Part I - mean values of estimated conditional expectation functions. Total of replicas equals to 1000. 4.3. DATA EXPLORATORY METHODS 65 0.0 0.2 0.4 0.6 0.8 1.0 1.2 −0.2 0.0 0.2 0.4 0.6 III− Just Biased distance conditional expectation Eseq(u)−E(u) Eseq(u)−E(Z) E(u)−E(Z) CI CI 0.0 0.2 0.4 0.6 0.8 1.0 1.2 −1.0 0.0 0.5 1.0 1.5 2.0 IV− Biased and Clustered distance conditional expectation Figure 4.5: Part II - mean values of estimated conditional expectation functions. Total of replicas equals to 1000. 66 CHAPTER 4. CLUSTERED AND BIASED MULTI-STAGE SAMPLING similar behaviour, as the choice of latest data points is randomly affected by previous measurements (not always the maximum). In this case, the CIs are slightly larger than in other cases, maybe because of the random selection of the cluster position. The analysis of the corresponding estimated standard deviation functions4, given by qb Vseq(u), qb V∗ seq(u) and qb V(u), was not so conclusive (see Figure 4.6). According to Schlather et al. (2004), we would expect an approximately constant variance when there is no dependence between data values and data locations, i.e. under random and just clustered sampling designs. However, this only happens when we estimate the second order moment about the theoretical overall expectation, the E(Z) chosen for our simulation study. This means making f(Z(x1), Z(x2)) = (Z(x1)−E(Z))2in the estimator defined in (4.2). Actually, the main pattern found in the variance’s plots is imposed by the clustering issue, responsible for smaller estimates of the variance for smaller lags5. Our function b Vseq(u) seems to under-estimate the theoretical variance, because it is restricted to data values of stage 2, typically with less variability. We decided then to focus on the conditional expectation functions. 4.4 Monte Carlo tests The widely used Monte Carlo significance testing was originally proposed by Barnard (1963) and its basic idea is as follows. Suppose H0is the null hypothesis about the model which generates Y={(xi, Z(xi)) : i= 1, ..., n}, and r1is an observed value of a real valued statistic R=h(Y), which has a distribution function F, possibly mathematically intractable. Moreover, suppose we agree to reject H0 for a large value of r1. Hence, we can use pseudo-random numbers to simulate a random sample 4We decide not to plot the CIs, because they would not add new interpretation hints. 5Note that our simulation study involved a cluster of diameter approximately equal to 0.1. 4.4. MONTE CARLO TESTS 67 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.8 1.2 1.6 I− Random distance conditional std deviation SDseq(u) SDseq*(u) SD(u) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.8 1.2 1.6 II− Clustered distance conditional std deviation 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.8 1.2 1.6 III− Biased distance conditional std deviation 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.0 1.5 2.0 2.5 IV− Biased and Clustered distance conditional std deviation Figure 4.6: Mean values of estimated conditional standard deviation functions (theoretical stddev is displayed in grey). Total of replicas equals to 1000. 68 CHAPTER 4. CLUSTERED AND BIASED MULTI-STAGE SAMPLING r2, ..., rmof m−1 observations from distribution Fand to construct a test by comparing these simulated values with r1. If Fis continuous and k= 1 + #{j: j= 2, ..., m and r1> rj}, then H0will be rejected at the k/m attained significance level, since the rank of r1is uniformly distributed on the integers 1, ..., m when H0is true. See Besag and Diggle (1977) for a general discussion of Monte Carlo tests. Note that the parametric bootstrap techniques work in a similar way to those described in here (see e.g. Gentle 2002 or Hall 1991). In our work, we are interest in a test for the hypothesis that a given data set does not incorporate sequential biasing, so that we shall define H0:Eseq(u)−E(u) = 0.(4.4) Under this hypothesis, the spatial process can be generated by sampling a random field Zat the given locations xi,i= 1, ..., n, with no sequential dependence. In this way, we can simulate m−1 further data sets under H0, and define rjto be a measure of discrepancy between b Ej seq(u) and b Ej(u) over the whole range of u. For example, our test statistic can be given by the integrated squared difference rj=Z{b Ej seq(u)−b Ej(u)}2du. (4.5) We can then proceed to a formal test based on the rank of r1amongst rj, because under H0all ranking of r1are equiprobable. Bear in mind that mis rather smaller than might perhaps be expected, in contrast with the much larger sample which would be needed for accurate estimation of F, the distribution function of R. According to Hope (1968), for a one-sided test at the conventional 5% level of significance, m= 100 is suitable. One may not wish to proceed directly to formal testing. A preliminary rough visual guide to address the problem being investigated can be provided by means of the well-known “simulation envelopes” plot. Testing involves comparing an 4.4. MONTE CARLO TESTS 69 observed test statistic with samples from the model under consideration. Consequently, this visual approach is based directly on the variation in estimates obtained from data generated from the model. The maximum and minimum of the total m−1 independent simulations allow the definition of upper and lower envelopes. See Figure 4.7, for an example. Diggle (2003) emphasises the use of such a plot as a visual aid to interpretation. Comparison of the observed6curve b E1 seq(u)−b E1(u) with that expected from a random arrangement of b Ej seq(u)−b Ej(u), j= 2, ..., m allows an assessment of the overall degree of coverage. If the observed curve lies between the two envelopes, this suggests the acceptance of hypothesis H0given in (4.4). If the observed curve exceeds the envelopes for some distances u, this is an initial and informal indication of the possibility of H0rejection. Anyway, in our case, we prefer to deepen analysis and to proceed with a formal Monte Carlo test. In Sections 4.4.1 and 4.4.2, we take for granted some model assumptions. In Section 4.4.3, we examine an alternative non-parametric approach, applying some randomization test ideas. 4.4.1 Example of a simulated data set We first emphasize the distinction between our proposal and the one described in Schlather et al. (2004), through a simple simulated example. We generate a sample data set {(xi, Z(xi)) : i= 1, ..., n}, incorporating some sequential bias but not forming any obvious cluster of spatial locations. This is our observed data set. We need now to simulate 99 further data sets under H0in (4.4). The idea is to fix the sample locations xiand to generate Gaussian data on them. Bear in mind that this normally requires the estimation of the spatial dependency structure. The accuracy of this estimation may strongly impact the results of Monte Carlo testing, as we shall see in the remainder of this Chapter. In the left panel of Figure 4.7, we plot the observed b Eseq(u)−b E(u) against 6The one derived from the observed data set {(xi, Z(xi)) : i= 1, ..., n}. 70 CHAPTER 4. CLUSTERED AND BIASED MULTI-STAGE SAMPLING 0.0 0.2 0.4 0.6 0.8 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 envelopes for Eseq(u)−E(u) distance expectation 0.0 0.2 0.4 0.6 0.8 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 envelopes for E(u)−E(Z) distance expectation Figure 4.7: Simulation envelopes of Eseq(u)−E(u) and E(u)−E(Z) for a simulated data set, with sequential bias but not clustered: data (solid curve); upper and lower envelopes from 99 simulations of a random field (dashed curves). the corresponding 99 simulations of a random field, which points to a possible rejection of H0. This rejection was confirmed with a formal test based on (4.5) at the 5% level of significance. However, the simulation envelopes of E(u)−E(Z) plotted in the right panel of Figure 4.7 could lead to the acceptance of H0. It is an example of just one simulated data set, but it suggests some caution when using the conditional expectation functions to detect dependency and, whenever applicable, we should take advantage of the possibly available multistage information7. 4.4.2 Rongelap island’s data We now describe the application of our methods to the data set from Rongelap island, whose sampling design has inspired part of the research work described in this Chapter. This island is located in the Pacific Ocean approximately 4000 kilometres south- 7Information about when the collection of the sample data occurred, namely time labels. 4.4. MONTE CARLO TESTS 71 west of Hawaii. The data were collected for the analysis of current levels of radioactivity contamination that resulted from a nuclear weapons testing programme during the 1950s. The scientific problem has been the estimation of the maximum level of radioactivity over the island, as part of a wider investigation to decide whether Rongelap can safely be resettled. See Diggle, Tawn and Moyeed (1998) or Diggle, Harper and Simon (1997) for more detail on these data. −6000 −5000 −4000 −3000 −2000 −1000 0 −3000 −2000 −1000 0 x y time t0 time t1 Figure 4.8: Rongelap’s island: two-stage strategy of uniform and clustered samples. The sampling design defined for data collection is illustrated in Figure 4.8. It started with a coarse grid of 63 locations and ended up with 98 additional measurements within four fine grids. These locations are identified by time label t0and t1, respectively. As this process involved two-stage of uniform and clustered samples, we wonder about the impact on conclusions from a standard analysis that does not account for either of these features. To proceed our analysis, the methods are applied to transformed data with a constant variance as described next. 72 CHAPTER 4. CLUSTERED AND BIASED MULTI-STAGE SAMPLING Our variables of interest from Rongelap data set are the spatial coordinates xi, the counts Yiof radioactive emissions at each location, the length liof time over which the counts are recorded, and the stage in sampling measurements were made. The total sample size is n= 161. Note that the Yiare treated as realizations of mutually independent Poisson random variables with expectations liλ(xi), where λ(x) measures the local radioactivity at location x. We chose the data transformation Zi=pYi/lito make the variability more consistent and more Gaussian8. It was found convenient to start with the maximum likelihood estimation of the spatial dependency structure. So, we consider a variogram estimator, obtained by using the coarse data and, derived through restricted maximum likelihood (REML). In Figure 4.9, we present the results of our Monte Carlo test. We generate 99 simulations of a random field over the total 157 distinct9locations of Rongelap’s island. From this plot, we would not reject the null hypothesis, as confirmed through a formal test. So, we would tend to refuse the existence of sequential bias. However, when one replaces this variogram estimator for one of those proposed in Section 4.5, this tendency is not that clear. 4.4.3 Randomization tests The previous approach requires model assumptions, like Gaussianity of data. If one wishes to avoid it, an alternative Monte Carlo method can be supported by the theory of randomization tests. The basic idea is to calculate a test statistic from the observed data, and then reshuffle the data a large number of times, recalculating the test statistic for each iteration. These statistics are used as before to generate a distribution of values. The observed value can be compared to the distribution to see whether the observed case is a tail value, i.e. an event that is unlikely to occur through chance. 8According to Delta method used to estimate a variance of a transformed parameter, one has Var[G(T)] ≃Var[T]×(G0(µ))2=const, where T=Y/l, E[T] = Var[T] = µand G(T) = √T. 9Four locations were overlapped in fine and coarse grids. 4.4. MONTE CARLO TESTS 73 0 1000 2000 3000 4000 −0.6 −0.4 −0.2 0.0 0.2 0.4 0.6 envelopes for Eseq(u)−E(u) distance expectation Figure 4.9: Simulation envelopes of Eseq(u)−E(u) for Rongelap island’s data, with bγobtained through REML: data (solid curve); upper and lower envelopes from 99 simulations of a random field (dashed curves). Often in comparisons between groups, like in some examples of biology applications (Manly 1991), only the group memberships are randomized while the same set of measurements are maintained. The latter tests are sometimes referred to as permutation ones, because the randomization can be done by reordering the positions of elements in an array. Monte Carlo tests evaluated by permutation are also quite applied to geostatistical data. Examples are the Mantel’s correlation test of two association matrices, that may refer to distances, and the Moran’s I test based on an empirical spatial autocorrelation coefficient. See Cliff and Ord (1981) for more detail on Moran’s I test. 80 CHAPTER 4. CLUSTERED AND BIASED MULTI-STAGE SAMPLING Figure 4.13: Behaviour of bγunder clustered sampling (no sequential biased). Mathe. Kernel RC NewRC Pool. 0.0 1.0 2.0 3.0 h<=0.6 Mathe. Kernel RC NewRC Pool. 0.0 0.2 0.4 0.6 0.8 h<=0.3 Mathe. Kernel RC NewRC Pool. 0.00 0.10 0.20 h<=0.2 Mathe. Kernel RC NewRC Pool. 0.00 0.02 0.04 h<=0.1 Figure 4.14: Behaviour of bγunder sequential biased sampling (no clustered). Mathe. Kernel RC NewRC Pool. 0.0 1.0 2.0 3.0 h<=0.6 Mathe. Kernel RC NewRC Pool. 0.0 0.2 0.4 0.6 0.8 h<=0.3 Mathe. Kernel RC NewRC Pool. 0.00 0.05 0.10 0.15 0.20 h<=0.2 Mathe. Kernel RC NewRC Pool. 0.000 0.010 0.020 0.030 h<=0.1 4.5. NON-STANDARD SAMPLING CORRECTORS 81 4.5.4 Rongelap island’s data We conclude this Chapter by proceeding with the assessment of the proposed variogram estimators on Rongelap data, when testing for the presence of sequential dependency. In Section 4.4.2, the Monte Carlo test suggested for the Rongelap data has employed a variogram estimator derived through maximum likelihood. This estimation was required to model the spatial dependency and to generate 99 simulations of a random field according to the null hypothesis H0in (4.4). We now investigate the influence of adopting the new variogram estimators instead, bearing in mind that we must use their corresponding valid versions. These are obtained by fitting the empirical estimators introduced in Sections 4.5.1 and 4.5.2 to a permissible variogram given by Bochner’s theorem in equation (3.4). Both estimators, with similar outcomes, were applied to all data from coarse and fine grids. In Figure 4.15, we choose to illustrate the Monte Carlo test related to the valid Pooled estimator. The observed curve b Eseq(u)−b E(u) is outside the simulation envelopes for small and large lags. According to these results, the presence of sequential dependency in the samples should not be totally excluded. Actually, the rejection of H0in (4.4) was confirmed with a formal test based on (4.5) at the 5% level of significance. As a last note, we highlight that this data set underlies some characteristics, like a low spatial variance and locations almost forming a straight line due to the island’s layout, requiring a careful estimation. Consequently, corrector methods like the ones proposed are advised. 82 CHAPTER 4. CLUSTERED AND BIASED MULTI-STAGE SAMPLING 0 1000 2000 3000 4000 −0.6 −0.4 −0.2 0.0 0.2 0.4 0.6 envelopes for Eseq(u)−E(u) distance expectation Figure 4.15: Simulation envelopes of Eseq(u)−E(u) for Rongelap island’s data, using the valid Pooled variogram estimator: data (solid curve); upper and lower envelopes from 99 simulations of a random field (dashed curves). Chapter 5 Properties of a kernel variogram estimator for clustered data 5.1 Introduction This Chapter is dedicated to the theoretical study of the new kernel variogram estimator proposed in Section 4.5.1 for clustered data, proving its asymptotic unbiasedness and consistency. Additionally, we shall propose optimal values for its unknown smoothing parameters, two user-adjustable quantities that affect the estimator’s performance. Bearing in mind that {Z(x) : x∈D⊂IRd}is an intrinsic and isotropic random process and denoting by Z(x1), ..., Z(xn) the values of the process observed at spatial locations x1, ..., xn, the suggested variogram estimator is defined as follows: bγ(u) = Pn i=1 Pn j=1 1 √ni×nj×K³u−kxi−xjk h´[Z(xi)−Z(xj)]2 2Pn i=1 Pn j=1 1 √ni×nj×K³u−kxi−xjk h´, u ≥0,(5.1) where ni=PkI{kxi−xkk≤δ}and nj=PkI{kxj−xkk≤δ};hand δrepresent the bandwidth and neighbourhood radius selectors, respectively. 83 84 CHAPTER 5. PROPERTIES OF bγ(.)ROBUST TO CLUSTERS Declustering methods are quite intuitive, and their need is well recognized in the spatial statistics literature to estimate spatially representative mean trends for clustered data (see e.g. Goovaerts (1997) and Isaaks and Srivastava (1989); or Dubois and Saisana (2002) for a comparison of classical declustering methods). In contrast, the corresponding need for the reliable estimation of the second-order spatial structures is not normally considered. The presence of clustered sample data is, however, not negligible at all as shown in Chapter 4. See for example Figure 4.2 or Table 4.1, which exhibit the decay of traditional variogram estimators under unequal samples density. Some of the main reasons for clustering of sample locations are: •External factors, like selection of locations conditional on specific geographic or demographic spots. •The need to better characterize short-range variability, requiring a denser sampling, but sometimes too costly to cover the whole observation region. •Adoption of a denser sampling in areas that are deemed critical. For example, the search of maximum values based on some prior knowledge. The recent paper from Kovitz and Christakos (2004) concerns the clustering issue and the estimation of the second-order structure. These authors suggest a modified form of Matheron’s estimator that also incorporates some declustering weights, but based on zones of proximity. Each zone of any data point is defined by the area of the Voronoi polygon that contains all points closer to that interior data point than to any other data point. The performance of this modified estimator of the variogram is analysed in terms of a numerical application. In our case, we prove that our variogram estimator enjoys good asymptotic properties. A short version of the preliminary theoretical result is found in Menezes, Garcia-Soid´an and Febrero-Bande (2004). As this estimator requires the selection of the bandwidth hand the radius δ, we recommend: the first will be treated via the MSE, i.e. the minimum square error; and the latter will result from the 5.2. ASSUMPTIONS 85 analysis of the density estimation derived on the observation region. The main contributions of this work are described in the extended version Menezes, Garcia- Soid´an and Febrero-Bande (2005b). The remainder of this Chapter is organized as follows. We first introduce additional notation and summarize the main assumptions considered in our asymptotic study. We then include comments about the neighbourhood radius selector. Next, the fundamental properties of the proposed non-parametric variogram estimator are established and corresponding proofs are developed. The results derived for bias and variance are used for the optimal bandwidth selection. We end with some numerical studies and implementation details about the proposed estimator. 5.2 Assumptions To ensure estimation consistency, we follow the strategy proposed by Hall et al. (1994), and recently, adopted by Garcia-Soid´an et al. (2004), according to which the observation region is considered to be increasing. Then, (A1) We start by assuming D=Dn=λD0where λ=λnmay diverge to +∞ and D0⊂IRdis a bounded and fixed region. (A2) Additionally, a random design is assumed for spatial locations xi=λvi, i = 1...n, where viis a realization of a random sample Vifrom f0, the density function defined on D0. (A3) For all v∈D0and for some positive constants d1and d2, one has d1≤f0(v)≤d2.1 (A4) γ(.) admits three continuous derivatives in a neighbourhood of u, for all u > 0. 1This allow us to guarantee that 0 <Rf0(x)idx<+∞, i = 2,3,4. 86 CHAPTER 5. PROPERTIES OF bγ(.)ROBUST TO CLUSTERS (A5) There is a bounded and continuously differentiable function g:IR3d→IR satisfying that Cov [(Z(xi)−Z(xj))2,(Z(xk)−Z(xl))2] = g(xi−xj,xi− xk,xi−xl).We assume lim kx2k≥rWkx3k≥r|g(x1,x2,x3)|= 0,where 0 < r < +∞ (A6) With respect to convergence rates, it is assumed that lim n→∞ {h+λ−1+λdn−1+ (nh)−1}= 0 (A7) Take δ=λa, where δis the neighbourhood radius in Dspace and a > 0 is the equivalent in D0. We assume that ahas an upper bound. Bear in mind that, in the context of a Gaussian process, one has Cov £(Z(xi)−Z(xj))2,(Z(xk)−Z(xl))2¤= = 2 [γ(kxi−xkk) + γ(kxj−xlk)−γ(kxi−xlk)−γ(kxj−xkk)]2 and, afterwards, one may take g(x1,x2,x3) = 2 [γ(kx2k) + γ(kx3−x1k)−γ(kx3k)−γ(kx2−x1k)]2 so that condition (A5) is satisfied provided that the variogram is bounded and has an asymptotic range. Thus a model with no finite range, such as the exponential, is acceptable. We are not considering unbounded variograms, such as the linear, but in real data applications we think it reasonable to restrict to a bounded spatial correlation. 5.3 Neighbourhood radius selector We first apply standard techniques of exploratory data analysis to gain a better understanding on an advisable value for δ. Namely, we used some elementary 5.3. NEIGHBOURHOOD RADIUS SELECTOR 87 theory of spatial point patterns to detect the presence of clusters (Diggle 2003). In this context, two useful functions are the cumulative distribution functions of point-to-point and origin-to-point nearest neighbour distances, Gand Frespectively. Suppose a spatial point pattern dataset with npoints. Let didenote the distance from the ith point to the closest of the other n−1 points. For a grid of ksampling origins, let eidenote the distance from the ith origin to the closest of the npoints. Then, these functions may be derived as b G(u) = n−1X i I{di≤u}and b F(u) = k−1X i I{ei≤u} where I{.}is the indicator function. The estimates of Gor Fcan be used for formal inference purposes about the pattern, when compared to the true value of Gor Ffor a completely random (Poisson) point process, which are G(u) = F(u) = 1 −exp(−λπu2) where λis the intensity (expected number of points per unit area). In Figure 5.1, we exemplify a graphic diagnostic with three distinct spatial models. The first one represents the example of the complete spatial randomness (CSR), where the locations within the unit square were obtained from an uniform distribution. The second and third models were obtained from a mixture of uniform and beta distributions, such that one or two strong clusters were achieved. For each of these models, we plot the estimated Gfunction, as well as, b Gand b Fagainst each other. As expected, the estimates of Gand Fpresent similar values under the first model. In our exploratory analysis, it was found convenient to consider the observation region to be defined in such way that edge effects can safely be ignored. In any case, if one is analysing some clustered area, it is reasonable to presume the cluster itself is not too close to the borders. 88 CHAPTER 5. PROPERTIES OF bγ(.)ROBUST TO CLUSTERS 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.4 0.8 X Y 0.00 0.05 0.10 0.15 0.20 0.0 0.4 0.8 distance estimated G 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.4 0.8 estimated G estimated F 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.4 0.8 X Y 0.00 0.05 0.10 0.15 0.20 0.2 0.6 1.0 distance estimated G 0.2 0.4 0.6 0.8 1.0 0.0 0.4 0.8 estimated G estimated F 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.4 0.8 X Y 0.00 0.05 0.10 0.15 0.20 0.0 0.4 0.8 distance estimated G 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.4 0.8 estimated G estimated F Figure 5.1: b Fand b Gunder three distinct models: CSR, 1 cluster and 2 clusters. 5.3. NEIGHBOURHOOD RADIUS SELECTOR 89 Another useful empirical function to summarise an observed pattern is the K function (variously called “Ripley’s K-function” or the “reduced second moment function”) of a stationary point process. This is defined as the expected number of additional random points within a distance uof a typical random point, divided by the overall intensity of the points. The Kfunction is determined by the second order moment properties of the point process. Under CSR, one expects K(u) = πu2. For any of the previous empirical functions, Monte Carlo simulations can be used to test hypotheses and construct confidence intervals. Testing involves comparing an observed test statistic with samples from the model under consideration, in this case a CSR model. Confidence intervals can be based directly on the variation in estimates observed in data generated from the model. The maximum and minimum of the total independent simulations allow the definition of upper and lower envelopes. The observed test statistic may be compared against these simulation envelopes (Diggle 2003). In Figure 5.2, one finds examples of Kestimates. The dashed curve, from the graphics in the middle side panels, identifies the πu2 curve. The graphics, found in the right side panels, represent the corresponding Monte Carlo tests for the expression qK(u) π−u= 0. The upper and lower envelopes (dashed curves) were derived from 99 CSR simulations. As expected, the second and third models suggest a rejection of the hypothesis “cluster absence”. Note that all these plots were produced using the R package Splancs, presented in Rowlingson and Diggle (1993). Some alternative methodologies for cluster analysis are directly motivated from techniques for density estimation (examples are Wong and Lane 1983, Silverman 1986, Cuevas, Febrero and Fraiman 2001). All of them are based on the natural idea of clusters correspond to modes or peaks in the underlying density function fon IRd. Very often the unknown theoretical function fis replaced by a non- 96 CHAPTER 5. PROPERTIES OF bγ(.)ROBUST TO CLUSTERS Consider the new random variables W2=V1−V2,...,Wn=V1−Vn. Keep also in mind that a realization of W2obeys to limλ→∞ kw2k= 0. This happens since Kis compactly supported, i.e. K(z) = 0 if |z|> C, meaning that λ−1(u−Ch)≤ kw2k ≤ λ−1(u+Ch). Then α=Z... ZK³u−λkw2k h´fn−1(w2, ..., wn) q1 + 2 Pk1≥2I{kwk1k≤ δ λ}+Pk1,k2≥2I{kwk1k≤ δ λ,kwk2−w2k≤ δ λ} dw2...dwn As a marginal distribution, fn−1(w2, ..., wn) can be written as Rfn(w1, ..., wn)dw1 =Rf0(w1)f0(w1−w2)...f0(w1−wn)dw1and, consequently, α=Z... ZK³u−λkw2k h´f0(w1)f0(w1−w2)...f0(w1−wn) s³1+2Pk1≥2I{kwk1k≤a}´+µPk1≥2 k2≥2 I{kwk1k≤a,kwk2−w2k≤a}¶dw1...dwn The expression under the square root of αmay be simplified as follows Ã1+2X k1≥2 I{kwk1k≤a}!+ÃX k1,k2≥2 I{kwk1k≤a,kwk2−w2k≤a}!= =Ã3+2X k1≥3 I{kwk1k≤a}!+   X k1=k2 k1≥2 I{kwk1k≤a}+X k16=k2 k1,k2≥2 I{kwk1k≤a,kwk2k≤a}   = =Ã3+2X k1≥3 I{kwk1k≤a}!+   1+3X k1≥3 I{kwk1k≤a}+X k16=k2 k1,k2≥3 I{kwk1k≤ak}I{kwk2k≤a}   = = 4 + 5 X k1≥3 I{kwk1k≤a}+X k16=k2 k1,k2≥3 I{kwk1k≤ak}I{kwk2k≤a} 5.4. BIAS OF bγ(.)ROBUST TO CLUSTERS 97 Then α=Z... ZK³u−λkw2k h´f0(w1)2f0(w1−w3)...f0(w1−wn) q4+5Pk1≥3I{kwk1k≤a}+Pk1,k2≥3 k16=k2 I{kwk1k≤a}I{kwk2k≤a} dw1...dwn Now, convert w2= (w(1), ..., w(d)) to spherical polar coordinates with the transformation w(i)=rcos θi i−1 Y j=0 sin θj where sin θ0= cos θd= 1, 0 ≤θd−1<2πand 0 ≤θi< π, for i= 1, ..., d −2. The corresponding Jacobian transformation is given by rd−1Jd(θ1, ..., θd−1) = rd−1(sin θ1)d−2(sin θ2)d−3... sin θd−2. Furthermore, suppose Pk1≥3I{kwk1k≤a}=kand apply some basic combinatory rules to obtain α=µZπ 0 ... Zπ 0Z2π 0Zm0 0 rd−1Jd(θ1, ..., θd−1)Kµu−λr h¶drdθ1...dθd−1¶. . n−2 X k=0 ¡n−2 k¢Rf0(w1)2H(a, w1)k(1 −H(a, w1))n−2−kdw1 p4+5k+k(k−1) where m0=sup{kxk:x∈D0}and H(a, w1) = Rkwk≤af0(w1−w)dw. Finally, with the following change of variable t=h−1(u−λr)⇒r=λ−1(u−th)⇒dr =−λ−1h dt, and as Kis compactly supported, the dominant term in αbecomes µZ... ZJd(θ1, ..., θd−1)dθ1...dθd−1¶ÃZu h u−m0λ h¡λ−1(u−th)¢d−1K(t)λ−1h dt!. . n−2 X k=0 ¡n−2 k¢Rf0(w1)2H(a, w1)k(1 −H(a, w1))n−2−kdw1 k+ 2 = 98 CHAPTER 5. PROPERTIES OF bγ(.)ROBUST TO CLUSTERS =ud−1Adλ−dhZf0(w1)2 n−2 X k=0 ¡n−2 k¢H(a, w1)k(1 −H(a, w1))n−2−k k+ 2 dw1 where Ad=Z... ZJd(θ1, ..., θd−1)dθ1...dθd−1(5.4) In the final expression of α, bear in mind that the integral performed over variable w1is bounded due to condition (A3). 5.4.2 Order of a2(u)−a1(u)γ(u)for u≥Ch Similarly, for u≥Ch, one has a2(u)−a1(u)γ(u) = X i6=j K³u−kxi−xjk h´ √ninj (γ(kxi−xjk)−γ(u)) and one may prove that the dominant term in here is given by n2β, where β= E  K³u−λkV1−V2k h´ qPk1,k2I{λkV1−Vk1k≤δ,λkV2−Vk2k≤δ} (γ(λkV1−V2k)−γ(u))  Let us again convert w2= (w(1), ..., w(d)) to spherical polar coordinates and perform the change of variable t=h−1(u−λr), to obtain that β=AdÃZu h u−m0λ h¡λ−1(u−th)¢d−1K(t) (γ(u−th)−γ(u)) λ−1h dt!. . n−2 X k=0 ¡n−2 k¢Rf0(w1)2H(a, w1)k(1 −H(a, w1))n−2−kdw1 k+ 2 Asymptotically, by using condition (A4), the expression (γ(u−th)−γ(u)) may be reduced to the second term of its Taylor expansion, i.e. γ00 (u) 2(−th)2. Then, the dominant term in βbecomes 1 2cKγ00(u)ud−1Adλ−dh3Zf0(w1)2 n−2 X k=0 ¡n−2 k¢H(a, w1)k(1 −H(a, w1))n−2−k k+ 2 dw1 where cK=Rz2K(z)dz. 5.5. VARIANCE OF bγ(.)ROBUST TO CLUSTERS 99 5.5 Variance of bγ(.)robust to clusters For the analysis of the asymptotic efficiency, it is important now to proceed with the derivation of variance for the proposed variogram estimator. A decreasing variance estimate means a growing efficiency of the estimator, as it will tend to be more accurate. Theorem 5.4 Assume the hypotheses required in Theorem 5.1. Additionally, suppose that assumptions (A5) and (A7) are satisfied. Then, for u≥Ch, one has Var [ˆγ(u)] = Bd(u)dK 2ud−1A2 d Ed(n, a)n−2λdh−1+Cd(u) A2 d Fd(n, a)n−1+ +Dd(u) 4A2 d Gd(n, a)λ−d+o(n−2λdh−1+n−1+λ−d+h4) where dK=R(K(z))2dz and Ad,Bd(u),Cd(u),Dd(u),Ed(n, a),Fd(n, a)and Gd(n, a)are as given in (5.4), (5.9), (5.10), (5.11), (5.6), (5.7) and (5.8), respectively. Remark 5.5 In a similar way as in Theorem 5.1, if we assume nsufficiently large then Theorem 5.4 holds for any u > 0. Let us start by considering that Var [ˆγ(u)] = Var [E [ˆγ(u)|V1, ..., Vn]] + E [Var [ˆγ(u)|V1, ..., Vn]] (5.5) By using Theorem 5.1 and Lemma 5.3, it is straightforward to see that for u≥Ch Var [E [ˆγ(u)|V1, ..., Vn]] = o(h4). We need now to check that Var [ˆγ(u)|V1, ..., Vn] = O(n−2λdh−1+n−1+λ−d) and, again by Lemma 5.3, it will lead us to the convergence rate of E [Var [ˆγ(u)|V1, ..., Vn]] . 100 CHAPTER 5. PROPERTIES OF bγ(.)ROBUST TO CLUSTERS Consequently, the convergence rate stated in Theorem 5.4 will be proved to be valid. About the detailed expression obtained for the conditional variance, we have that Var [ˆγ(u)|V1, ..., Vn] = E £(bγ(u)|V1, ..., Vn−E[bγ(u)|V1, ..., Vn])2¤= = E "ÃPi6=jwij(u) ((Z(xi)−Z(xj))2−E[(Z(xi)−Z(xj))2]) 2Pi6=jwij(u)!. .ÃPk6=lwkl(u) ((Z(xk)−Z(xl))2−E[(Z(xk)−Z(xl))2]) 2Pk6=lwkl(u)!#= =Pi6=jwij(u)Pk6=lwkl(u) Cov [(Z(xi)−Z(xj))2,(Z(xk)−Z(xl))2] ³2Pi6=jwij(u)´2= = (2a1(u))−2X i6=j k6=l Kµu−kxi−xjk h¶Kµu−kxk−xlk h¶. .1 √ninj 1 √nknl g(xi−xj,xi−xk,xi−xl) = 2e1(u)+4e2(u) + e3(u) 4(a1(u))2 where e1(u) = X i6=j K³u−kxi−xjk h´2g(xi−xj,0,xi−xj) ninj⇐(i=k∧j=l) e2(u) = X i6=j j6=l K³u−kxi−xjk h´K³u−kxi−xlk h´g(xi−xj,0,xi−xl) √ninj√ninl⇐(i=k) e3(u) = X i6=j,k,l j6=k,l k6=l K³u−kxi−xjk h´K³u−kxk−xlk h´g(xi−xj,xi−xk,xi−xl) √ninj√nknl 5.5. VARIANCE OF bγ(.)ROBUST TO CLUSTERS 101 Then, according to the results from Section 5.4.1 about a1(u), and those from Sections 5.5.1, 5.5.2 and 5.5.3 about e1(u), e2(u) and e3(u), respectively, we obtain: •e1(u) 2(a1(u))2=Bd(u)dK 2ud−1A2 d Ed(n, a)n−2λdh−1+o(n−2λdh−1) a.s. where Ed(n, a) = Rf0(w1)2Pn−2 k=0 (n−2 k)H(a,w1)k(1−H(a,w1))n−2−k (k+2)2dw1 µRf0(w1)2Pn−2 k=0 (n−2 k)H(a,w1)k(1−H(a,w1))n−2−k k+2 dw1¶2(5.6) •e2(u) (a1(u))2=Cd(u) A2 d Fd(n, a)n−1+o(n−1) a.s. where Fd(n, a) = Rf0(w1)3Pn−3 k=0 (n−3 k)H(a,w1)k(1−H(a,w1))n−3−k (k+3)2dw1 µRf0(w1)2Pn−2 k=0 (n−2 k)H(a,w1)k(1−H(a,w1))n−2−k k+2 dw1¶2(5.7) •e3(u) 4(a1(u))2=Dd(u) 4A2 d Gd(n, a)λ−d+o(λ−d) a.s. where Gd(n, a) = Rf0(w1)4Pn−4 k=0 (n−4 k)H(a,w1)k(1−H(a,w1))n−4−k (k+4)2dw1 µRf0(w1)2Pn−2 k=0 (n−2 k)H(a,w1)k(1−H(a,w1))n−2−k k+2 dw1¶2(5.8) The validity of the latter expressions demand that Ed,Fdand Gdare bounded, even for large n, which is proved numerically in Section 5.7.2. For the specific case of IR2, a supplemental simulation study for the dependency analysis of E2,F2and G2on nwas included in Section 5.7.3. 102 CHAPTER 5. PROPERTIES OF bγ(.)ROBUST TO CLUSTERS 5.5.1 Order of e1(u)for u≥Ch The dominant term of e1(u) is given by n2α1, where α1= E    K³u−λkV1−V2k h´2g(λ(V1−V2),0, λ(V1−V2)) Pk1,k2I{kV1−Vk1k≤a,kV2−Vk2k≤a}  = =Z···ZK³u−λkw2k h´2g(λw2,0, λw2)f0(w1)2f0(w1−w3)...f0(w1−wn) 4+5Pk1≥3I{kwk1k≤a}+Pk1,k2≥3 k16=k2 I{kwk1k≤a}I{kwk2k≤a} dw1...dwn As in Section 5.4.1, we may convert w2to spherical polar coordinates and make a change of variable to obtain α1=Zu h u−m0λ hZπ 0 ... Zπ 0Z2π 0 Jd(θ1, ..., θd−1)¡λ−1(u−th)¢d−1K(t)2λ−1h. .g Ã(u−th)(cos θ1, ..., d−1 Y j=0 sin θj),0,(u−th)(cos θ1, ..., d−1 Y j=0 sin θj)!dt dθ1...dθd−1. . n−2 X k=0 ¡n−2 k¢Rf0(w1)2H(a, w1)k(1 −H(a, w1))n−2−kdw1 (k+ 2)2 The dominant term will be given by α1=ud−1Bd(u)dKλ−dhZf0(w1)2 n−2 X k=0 ¡n−2 k¢H(a, w1)k(1 −H(a, w1))n−2−k (k+ 2)2dw1 where dK=R(K(z))2dz and Bd(u) = Zπ 0 ... Zπ 0Z2π 0 Jd(θ1, ..., θd−1). .g Ãu(cos θ1, ..., d−1 Y j=0 sin θj),0, u(cos θ1, ..., d−1 Y j=0 sin θj)!dθ1...dθd−1(5.9) 5.5. VARIANCE OF bγ(.)ROBUST TO CLUSTERS 103 5.5.2 Order of e2(u)for u≥Ch In here, three distinct indices i,j,l, are involved, thus the dominant term of e2(u) will be given by n3α2, where α2= E  K³u−λkV1−V2k h´K³u−λkV1−V3k h´g(λ(V1−V2),0, λ(V1−V3)) qPk1,k2I{kV1−Vk1k≤a,kV2−Vk2k≤a}qPk1,k2I{kV1−Vk1k≤a,kV3−Vk2k≤a}  For random variables Wi=V1−Vi, i = 2,3, as Kis compactly supported, we shall show that the dominant term of the expectation above can be reduced to those values kw2kand kw3ktending to 0. Then, it becomes α2=Z... ZK³u−λkw2k h´K³u−λkw3k h´g(λw2,0, λw3) q1+2Pk1≥2I{kwk1k≤a}+Pk1,k2≥2I{kwk1k≤a,kwk2−w2k≤a} . .f0(w1)f0(w1−w2)...f0(w1−wn) q1+2Pk1≥2I{kwk1k≤a}+Pk1,k2≥2I{kwk1k≤a,kwk2−w3k≤a} dw1...dwn The square root can actually be eliminated, as X k1,k2≥2 I{kwk1k≤a,kwk2−w3k≤a}=X k1,k2≥2 I{kwk1k≤a,kwk2−w2k≤a}=X k1,k2≥2 I{kwk1k≤a,kwk2k≤a} Furthermore 1+2X k1≥2 I{kwk1k≤a}+X k1,k2≥2 I{kwk1k≤a,kwk2k≤a}= =  1 + 2 X k1≥2 I{kwk1k≤a}+X k1=k2 k1≥2 I{kwk1k≤a,kwk1k≤a}  +   X k16=k2 k1,k2≥2 I{kwk1k≤a,kwk2k≤a}   = =Ã7 + 3 X k1≥4 I{kwk1k≤a}!+   2+4X k1≥4 I{kwk1k≤a}+X k1,k2≥4 k16= I{kwk1k≤a}I{kwk2k≤a}    104 CHAPTER 5. PROPERTIES OF bγ(.)ROBUST TO CLUSTERS According to previous results α2=Z... ZK³u−λkw2k h´K³u−λkw3k h´g(λw2,0, λw3) 9+7Pk1≥4I{kwk1k≤a}+Pk1,k2≥4,k16=k2I{kwk1k≤a}I{kwk2k≤a} . .f0(w1)3f0(w1−w4)...f0(w1−wn)dw1...dwn Now, convert w2= (w(1,1), ..., w(d,1)) and w3= (w(1,2), ..., w(d,2)) to spherical polar coordinates with the transformation w(i,1) =r1cos θi,1 i−1 Y j=0 sin θj,1and w(i,2) =r2cos θi,2 i−1 Y j=0 sin θj,2 where, for k= 1,2, sin θ0,k = cos θd,k = 1, 0 ≤θd−1,k <2πand 0 ≤θi,k < π, for i= 1, ..., d −2. The corresponding Jacobian transformations are given by rd−1 kJd(θ1,k, ..., θd−1,2) = rd−1 k(sin θ1,k)d−2(sin θ2,k)d−3... sin θd−2,k. In this way, α2=µZm0 0Zπ 0 ... Zπ 0Z2π 0Zm0 0Zπ 0 ... Zπ 0Z2π 0 rd−1 1Jd(θ1,1, ..., θd−1,1)rd−1 2. .Jd(θ1,2...θd−1,2)gÃλr1(cos θ1,1, ..., d−1 Y j=0 sin θj,1),0, λr2(cos θ1,2, ..., d−1 Y j=0 sin θj,2)!. . K µu−λr1 h¶Kµu−λr2 h¶dr1dr2dθ1,1...dθd−1,1dθ1,2...dθd−1,2¶. . n−3 X k=0 ¡n−3 k¢Rf0(w1)3H(a, w1)k(1 −H(a, w1))n−3−kdw1 9 + 7k+k(k−1) Finally, with the following changes of variable, for k= 1,2 tk=h−1(u−λrk)⇒rk=λ−1(u−tkh)⇒drk=−λ−1h dtk the dominant term becomes 5.5. VARIANCE OF bγ(.)ROBUST TO CLUSTERS 105 α2=λ−2h2Zπ 0 ... Zπ 0Z2π 0Zπ 0 ... Zπ 0Z2π 0 Jd(θ1,1, ..., θd−1,1)Jd(θ1,2...θd−1,2). .g Ãu(cos θ1,1, ..., d−1 Y j=0 sin θj,1),0, u(cos θ1,2, ..., d−1 Y j=0 sin θj,2)!dθ1,1...dθd−1,1dθ1,2...dθd−1,2. .ÃZu h u−λm0 hZu h u−λm0 h (λ−1(u−t1h))d−1(λ−1(u−t2h))d−1K(t1)K(t2)dt1dt2!. . n−3 X k=0 ¡n−3 k¢Rf0(w1)3H(a, w1)k(1 −H(a, w1))n−3−kdw1 9 + 7k+k(k−1) = =u2(d−1)Cd(u)λ−2dh2Zf0(w1)3 n−3 X k=0 ¡n−3 k¢H(a, w1)k(1 −H(a, w1))n−3−k (k+ 3)2dw1 where the integral over variable w1is bounded due to condition (A3), and Cd(u) = Zπ 0 ... Zπ 0Z2π 0Zπ 0 ... Zπ 0Z2π 0 Jd(θ1,1, ..., θd−1,1)Jd(θ1,2...θd−1,2). .g Ãu(cos θ1,1, ..., d−1 Y j=0 sin θj,1),0, u(cos θ1,2, ..., d−1 Y j=0 sin θj,2)!. .dθ1,1...dθd−1,1dθ1,2...dθd−1,2(5.10) 5.5.3 Order of e3(u)for u≥Ch In here, four distinct indices i,j,k,l, are involved, thus the dominant term of e3(u) will be given by n4α3, where α3is equal to E K³u−λkV1−V2k h´K³u−λkV3−V4k h´g(λ(V1−V2), λ(V1−V3), λ(V1−V4)) qPk1,k2I{kV1−Vk1k≤a,kV2−Vk2k≤a}qPk1,k2I{kV3−Vk1k≤a,kV4−Vk2k≤a}  112 CHAPTER 5. PROPERTIES OF bγ(.)ROBUST TO CLUSTERS The region Dis assumed to be equal to λD0, where D0is the bounded and fixed square unit. The new estimator is compared against the estimator of Matheron and the one using the Nadaraya-Watson kernel, given in (3.1) and (3.3), respectively. The symmetric Epanechnikov kernel was employed in the two previous kernel-type estimators. To obtain the optimal local bandwidth, we considered the optimal c0=8 9, and c1= 1. Then, this local bandwidth can be derived as a function of lag u, like hopt(u) = ·B2(u)dKE2(n, a) 2uA2 2c2 Kγ00(u)2¸1/5 n−2/9. The corresponding scale factor, from (5.12), is given by λ= (n8/9)1/2. We considered a sample size n= 100 and a theoretical exponential variogram with a nugget effect of 0.6, a sill of 1.336 and the corresponding range equal to 5.0. About the bandwidth derivation, note that cKidentifies the variance of the Epanechnikov kernel and dKidentifies the integral of this squared kernel. As we are in IR2, the A2expression in (5.4) is reduced to 2π. For the Gaussian case, B2(u) given in (5.9) can be approximated by 8A2γ(u)2. To estimate E2(n, a), given in (5.6), the existing integrals were numerically approximated to a sample average, as Zg(x)dx = E ·g(X) f(X)¸,where f(x) is the density function of X. Section 5.7.3 provides more detail about E2(n, a) estimation. Additionally, as the bandwidth derivation needs itself an estimation of the variogram, a rough parametric estimation was used for this purpose. To proceed with our simulation study, we generate a total of 100 independent data sets and, for each one, derive the integrated square error (ISE) between the estimator and the theoretical variogram. The ISE, defined as Rβ α[bγ(u)−γ(u)]2du, was approximated numerically through the trapezoid rule. In Table 5.1, the mean values of the resulting ISEs are compared for two distinct sampling designs: •A CSR model, where points are uniformly distributed on D; 5.7. NUMERICAL STUDIES 113 u≤0.6λ u ≤0.3λ u ≤0.2λ u ≤0.1λ CSR Matheron 1.270 0.943 0.819 0.763 NW kernel 0.527 0.314 0.276 0.291 RobCluster 0.500 0.307 0.276 0.298 CLUSTER Matheron 1.519 1.141 0.889 0.568 NW kernel 0.582 0.525 0.488 0.400 RobCluster 0.392 0.294 0.243 0.245 Table 5.1: Mean values of the standardized ISEs, from the empirical estimators. The total number of replicas is 100 and in each replica the total sample size is 100. •A clustered model, where 40% of the total points are gathered together into one sub-region of D. As the observation region Ddepends on λ, we decided to group the mean values of the ISEs into four classes of lags: (0,0.6λ), (0,0.3λ), (0,0.2λ) and (0,0.1λ). To easily compare columns, all ISE values were standardized by dividing them by the corresponding integral interval, β−α. According to Table 5.1, the new empirical estimator, named “RobCluster”, offers a better performance in the presence of clustered data. Under a CSR model, the Nadaraya-Watson kernel estimator and the new estimator present similar results, and better than those from Matheron’s proposal. We repeat this same experiment with the corresponding valid versions of the previous three empirical estimators, after fitting them to a permissible variogram defined by Bochner’s theorem in equation (3.4) (see Chapter 3). In Table 5.2, we summarize the mean values of the obtained ISE. Once more, one may confirm the better behaviour of the proposed variogram under clustered data. Alternatively, one might work with a global bandwidth (see Remark 5.6). In this case, the optimal expression for bandwidth hdoes not depend on lag u, as it depends instead on some integrals of u. Bear in mind, a global bandwidth is 114 CHAPTER 5. PROPERTIES OF bγ(.)ROBUST TO CLUSTERS u≤0.6λ u ≤0.3λ u ≤0.2λ u ≤0.1λ CSR Matheron 0.861 0.672 0.609 0.577 NW kernel 0.504 0.285 0.223 0.201 RobCluster 0.479 0.279 0.231 0.214 CLUSTER Matheron 0.925 0.848 0.764 0.519 NW kernel 0.490 0.475 0.473 0.385 RobCluster 0.364 0.278 0.229 0.217 Table 5.2: Mean values of the standardized ISE, from valid estimators. The total number of replicas is 100 and for each replica the total sample size is 100. expected to lead to faster simulations when compared to a local one, as it avoids a specific estimation for each lag u. The natural drawback is that it proposes a less accurate solution. 5.7.2 Analysis of Ed(n, a),Fd(n, a)and Gd(n, a)for large n The goal of the current simulation study is to understand how Ed,Fdand Gd, given in (5.6), (5.7) and (5.8), conduct themselves under a large sample size n. These three expressions share the common denominator ¡Rf0(w1)2S dw1¢2, where S= n−2 X k=0 ¡n−2 k¢H(a, w1)k(1 −H(a, w1))n−2−k k+ 2 (5.14) As H(a, w1) = Rkwk≤af0(w1−w)dw, we shall replace H(a, w1) by Hin S, with 0 < H < 1, and analyse the dependency of Son nand H. The results from this dependency analysis are summarized in Table 5.3. We wish to emphasize that the exact value of Hloses importance with increasing sample size n. In fact, the value derived for the standard deviation decreases when nincreases. The latter conveys that S=O(an), for some bounded sequence an. 5.7. NUMERICAL STUDIES 115 n / H 0.1 0.3 0.5 0.7 Std Dev 100 9.09E-02 3.26E-02 1.98E-02 1.42E-02 3.52E-02 500 1.96E-02 6.64E-03 3.99E-03 2.85E-03 7.74E-03 1000 9.91E-03 3.33E-03 2.00E-03 1.43E-03 3.91E-03 5000 2.00E-03 6.66E-04 4.00E-04 2.86E-04 7.89E-04 10000 9.99E-04 3.33E-04 2.00E-04 1.43E-04 3.95E-04 50000 2.00E-04 6.66E-05 4.00E-05 2.86E-05 7.89E-05 Table 5.3: Values obtained for Sin (5.14), when given nand H. In the last column, are values for the corresponding standard deviation of H, when given n. Let us now consider the following three quotients, for i= 2,3,4: Qi=Pn−i k=0 (n−i k)Hk(1−H)n−i−k (k+i)2 S2(5.15) Table 5.4 presents the values obtained for Q2,Q3and Q4, for the same previous values of H. One notes that these quotients tend to 1 with increasing sample size. This tendency may be observed, for any chosen probability H. Expressions (5.6), (5.7) and (5.8) may now be re-written as Ed(n, a) = Rf0(w1)2Q2S2dw1 ¡Rf0(w1)2S dw1¢2≃Rf0(w1)2O(an)2dw1 ¡Rf0(w1)2O(an)dw1¢2= =OÃRf0(w1)2dw1 ¡Rf0(w1)2dw1¢2! Fd(n, a) = Rf0(w1)3Q3S2dw1 ¡Rf0(w1)2S dw1¢2≃Rf0(w1)2O(an)2dw1 ¡Rf0(w1)3O(an)dw1¢2= =OÃRf0(w1)3dw1 ¡Rf0(w1)2dw1¢2! Gd(n, a) = Rf0(w1)4Q4S2dw1 ¡Rf0(w1)2S dw1¢2≃Rf0(w1)2O(an)2dw1 ¡Rf0(w1)4O(an)dw1¢2= 116 CHAPTER 5. PROPERTIES OF bγ(.)ROBUST TO CLUSTERS Q2n / H 0.1 0.3 0.5 0.7 100 1.08588 1.02326 1.00999 1.00428 500 1.01798 1.00467 1.00200 1.00086 1000 1.00900 1.00233 1.00100 1.00044 5000 1.00180 1.00047 1.00020 1.00009 10000 1.00090 1.00023 1.00010 1.00004 50000 1.00018 1.00005 1.00002 1.00001 Q3n / H 0.1 0.3 0.5 0.7 100 0.90084 0.97543 0.98960 0.99556 500 0.98157 0.99529 0.99798 0.99914 1000 0.99089 0.99765 0.99900 0.99957 5000 0.99820 0.99953 0.99980 0.99991 10000 0.99910 0.99977 0.99990 0.99996 50000 0.99982 0.99995 0.99998 0.99999 Q4n / H 0.1 0.3 0.5 0.7 100 0.76283 0.93098 0.96982 0.98696 500 0.94713 0.98603 0.99399 0.99742 1000 0.97328 0.99301 0.99700 0.99871 5000 0.99461 0.99860 0.99940 0.99974 10000 0.99730 0.99930 0.99970 0.99987 50000 0.99946 0.99986 0.99994 0.99997 Table 5.4: Values obtained for Q2,Q3and Q4in (5.15), when given nand H. 5.7. NUMERICAL STUDIES 117 =OÃRf0(w1)4dw1 ¡Rf0(w1)2dw1¢2! Assumption (A3) allows us to guarantee that 0 <Rf0(w1)idw1<+∞, i = 2,3,4. Consequently, the approximations derived above for Ed,Fdand Gdare of the exact order O(1) and, therefore, they are bounded. 5.7.3 Estimates of E2(n, a),F2(n, a)and G2(n, a) For the specific case of IR2, we now describe a supplemental simulation study for the dependency analysis of E2,F2and G2on n. We also suggest a numeric approximation for the expressions introduced in (5.6), (5.7) and (5.8). Bear in mind that these are defined in the region D0⊂IR2, which must be a bounded and fixed region. We have selected D0to be the unit square, [0,1] ×[0,1]. The density function for the spatial locations on D0is f0. To estimate E2,F2and G2, the corresponding integrals were numerically approximated to a sample average, as Zg(x)dx = E ·g(X) f(X)¸,where f(x) is the density function of X. For instance, b E2may be derived, as follows: b E2(n, a) = 1 NPN i=1 b f0(wi)Pn−2 k=0 (n−2 k)b H(a,wi)k(1−b H(a,wi))n−2−k (k+2)2 µ1 NPN i=1 b f0(wi)Pn−2 k=0 (n−2 k)b H(a,wi)k(1−b H(a,wi))n−2−k k+2 ¶2 where •nis the number of original sampled points; we chose n= 100,200,400; •Nis is the number of extra points generated from density f0and needed for the integral approximation; we chose N= 5000; •b H(a, wi) = ni n, being nithe number of original sampled points within the circle of center wiand radius a; 118 CHAPTER 5. PROPERTIES OF bγ(.)ROBUST TO CLUSTERS CSR n b E2b F2b G2 100 2.142(0.102) 0.927(0.009) 0.497(0.031) 200 1.945(0.056) 0.980(0.006) 0.560(0.019) 400 1.831(0.028) 1.005(0.003) 0.610(0.011) CLUSTER n b E2b F2b G2 100 1.971(1.166) 1.206(0.098) 2.571(1.463) 200 2.096(0.273) 0.892(0.042) 1.010(0.231) 400 2.111(0.097) 0.901(0.012) 0.607(0.054) Table 5.5: Mean values of b E2,b F2and b G2, obtained from a total of 100 independent samples. The corresponding standard deviations are given between brackets. •b f0results from a non parametric density estimation of the spatial locations in D0; we adopted a bivariate kernel-type estimator; The other two estimates, b F2and b G2, may be obtained in a very similar way. We started with a complete spatial randomness (CSR) design. So, we generated nlocations uniformly distributed on D0. This procedure was repeated to obtain 100 independent samples. Table 5.5 presents the average of those 100 replicas and the corresponding standard deviation. The next simulation included one clustered area on D0, where we forced a minimum of 60 points to be restricted to a small square, with area equal to 0.16 × 0.16 instead of the original 1 ×1, and a center randomly chosen. The results were also included in Table 5.5. The main conclusion from both simulations appears to be the absence of an obvious tendency with increasing of sample size. In any case, for any of the three approximations of E2,F2and G2, the standard deviation clearly decreases with increasing of sample size, so that the mean value provides a good estimate of the unknown term. Chapter 6 Assessing the effect of preferential sampling 6.1 Introduction As stated before, in geostatistics, in both prediction and inference contexts, it is commonly assumed that the selection of the sampling locations does not depend on the values of the spatial variable (Diggle et al. 2003). Additionally, most techniques are based on the assumption, possibly tacit, of sampling locations being uniformly distributed over the observed region. In Chapter 4, we assess the effect of the failure of the earlier assumptions concerning the estimation of the correlation structure in the specific case of multistage collection of spatial data. The appraisal of biased data in later stages, conditional on data values from earlier stages, is considered. As the presence of clusters is a natural consequence of non-uniform locations distribution, we propose a kernel estimator robust to clusters. Then, in Chapter 5, we proceed with the theoretical study of the suggested estimator. We now intend to introduce a formal definition directly related to the failure of the independency assumption, and not restricted to multi-stage sample collection. Suppose that, in the nature of the sampling process, involving as it does 119 120 CHAPTER 6. ASSESSING PREFERENTIAL SAMPLING Unobserved field SObserved data (xi, yi) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Figure 6.1: Example of an unobserved field process (highest values represented by lightest colors) and the corresponding observed sample data set (highest values of yirepresented by largest bullets). the search for maximum values, there exists an underlying stochastic relationship between data and locations, then one has preferential sampling. Following the notation used by Diggle et al. (2003), we shall consider that the data for analysis are of the form (xi, yi) : i= 1, ..., n, where x1, ..., xnare locations within an observation region D⊂IR2and y1, ..., ynare measurements associated with these locations. The {xi:i= 1, ..., n}is the sampling design and yiis assumed to be a realization of Yi=Y(xi), where {Y(x) : x∈D}is the measurement process. We also assume the existence of an unobserved field process {S(x) : x∈D}, usually regarded as our goal of prediction. Often, Yican be considered as a noisy version of the underlying random variable S(xi), the value at location xiof process S(.). Figure 6.1, in the left hand panel, illustrates an example of a true field Sand, in the right hand panel, the corresponding sample data set (xi, yi). 6.2. CLASS OF LOG-GAUSSIAN COX PROCESSES 121 Remark 6.1 Under preferential sampling, the sampling design process is assumed to be stochastically dependent of the field process S(.). Consequently, the corresponding geostatistical model (specified by the joint distribution of the processes involved) must take into account the conditional distribution of the sampling design. Considering the presence of preferential sampling under Gaussian assumptions, we shall propose a model-based approach for spatial prediction. This new parametric model will be founded on a flexible class of log-Gaussian Cox processes to be introduced in the next Section. The remainder of this Chapter is devoted primarily to the analysis of the consequences of ignoring preferential sampling and adopting the classic geostatistical methods. In Chapter 7, we proceed with the likelihood inference to estimate our model parameters. To terminate this introductory Section, we want to emphasize the distinction between clustered and preferential sampling. As already mentioned, the clustering of locations may be due to the existence of specific geographic or demographic spots, or they may even be used to describe short-range variability better. These are good examples showing that clustered sampling may not imply preferential sampling. On the contrary, the opposite implication tends to occur, as preferred sample locations normally occur in concentrated areas. For example, some prior scientific knowledge about S(.), such as the expected local ore grade in mine exploration, may cause the concentration of samples in areas with atypically large values. 6.2 A class of log-Gaussian Cox processes The sample locations xare nothing more than realizations of a point process P. Under complete spatial randomness, the point process modelling is typically based on some homogeneous Poisson process. In our case, however, we are not 128 CHAPTER 6. ASSESSING PREFERENTIAL SAMPLING 0.0 0.2 0.4 0.6 0.00 0.10 0.20 0.30 beta=0 distance std dev / bias std dev std error bias 0.0 0.2 0.4 0.6 0.00 0.10 0.20 0.30 beta=1 distance std dev / bias 0.0 0.2 0.4 0.6 0.00 0.10 0.20 0.30 beta=2 distance std dev / bias Figure 6.3: Influence of βon variogram estimation. Comparison of the estimation bias and the corresponding approximation variability, given by the standard deviation and the standard error. The simulation consisted of 500 replications, each with a sample size of 100. 6.3. EFFECT ON VARIOGRAM ESTIMATION 129 σ2→0), then γ(.) and eγM(.) become similar. Consequently, the degree of preferability βseems to lose importance under too much white noise, even though the clustering effect persists. 6.3.3 Clustered versus preferential We now aim to carry out a numerical study to measure the impact of clustered sampling on variogram estimation, both with and without, being also preferential. Three distinct sampling designs are to be considered: completely spatial randomness (CSR), preferential and just clustered. As expected, the CSR sample is obtained for βequal to 0. A pairwise sample generation is adopted, otherwise: we first generate a preferential sampling data set for βequal to 2; we then keep previous locations and generate new multivariate Gaussian data on them to obtain a clustered but non-preferential data set. Remember that eγM(.) is derived from a Monte Carlo average of several independent bγM(.). For all results included in the foregoing Section, we choose as bγM(.) the classic estimator from Matheron, given in (3.1). Here, we intend to compare again the performance of three next variogram estimators: the classic, the NW kernel and the robust to clusters (the two latter are given in (3.3) and (5.1), respectively). In the following simulations, the variance parameters were changed to more significant values, becoming τ2= 0.25 and σ2= 2.25. The simulation consisted of 500 independent replicas, each with a sample size of 100. For each independent data set, we first derive the integrated squared error (ISE) between each empirical estimator and the theoretical variogram, as defined in (3.6). The results are summarized in the boxplots in Figure 6.4, through the quartiles of the ISE values found, when taking all lags smaller than 0.6. These boxplots confirm the results already observed in Chapter 4 (see e.g. Table 4.1), exhibiting the positive contribution of the proposed estimator when sampling is clustered, independently of whether preferential or not. 130 CHAPTER 6. ASSESSING PREFERENTIAL SAMPLING Classic Kernel RobCl 0.00 0.10 0.20 beta=0 ISE 0.069 0.052 0.049 Classic Kernel RobCl 0.00 0.10 0.20 beta=2 ISE 0.086 0.066 0.051 Classic Kernel RobCl 0.00 0.10 0.20 beta=2 just clustered ISE 0.099 0.064 0.058 Figure 6.4: Boxplot of the evaluated ISE from three empirical estimators: Classic, NW kernel and Robust to Clusters. Three sampling designs are considered: CSR (β= 0), preferential (β= 2) and clustered (non-preferential). To gain more information concerning the behaviour of these estimators, we then proceed with some efficiency assessment, by comparing their variances and their mean squared errors (MSE). Recall that the latter is defined as MSE [ˆγM(u)] = (Bias [ˆγM(u)])2+ Var [ˆγM(u)] . In Figure 6.5, in the left hand panels, we plot the square roots of variances, including biases, for our three estimators. In the right hand panels the corresponding square roots of the MSE’s are likewise plotted. 6.3. EFFECT ON VARIOGRAM ESTIMATION 131 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.4 0.8 beta=0 distance STDEV and BIAS Classic NW kernel RobClusters 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.4 0.8 beta=0 distance sqrt(MSE) Classic NW kernel RobClusters 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.4 0.8 beta=2 distance STDEV and BIAS Classic NW kernel RobClusters 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.4 0.8 beta=2 distance sqrt(MSE) Classic NW kernel RobClusters 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.4 0.8 beta=2 just clustered distance STDEV and BIAS Classic NW kernel RobClusters 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.4 0.8 beta=2 just clustered distance sqrt(MSE) Classic NW kernel RobClusters Figure 6.5: Comparison of variogram estimators through their efficiency (measure in terms of √V ar and √Mse) when sampling is clustered, both with and without, being also preferential. The estimators bias is also plotted (dashed-lines). It is being considered τ2= 0.25 and σ2= 2.25. 132 CHAPTER 6. ASSESSING PREFERENTIAL SAMPLING The analysis of these graphics leads us to the following results. Under CSR, the performance of the two kernel estimators is confirmed to be only slightly better than under the classic one. Actually, the latter presents a little more variance for nugget estimation, and for γ(.) estimation regarding those lags larger than the range value, i.e. 0.4. The top-left panel in Figure 6.5 also reveals a minor bias, nearly meaningless in the case of the classic estimator, that it is sometimes referred to as the smoothing bias. Here it is of interest to recall that, in classical geostatistics3, the bias term goes to zero as nincreases, as kernel estimators are asymptotically unbiased. Under clustered sampling, whether or not also preferential, the two kernel estimators always perform significantly better than the classic one. Note that the smaller variance of these estimators, when sampling is also preferential, is an indirect consequence of the smaller variance of the data values under this sampling design. Comparing the two kernel estimators through the MSE and, through the ISE, that robust to clusters always performs at least slightly better. This occurs even when this estimator presents a larger variance, as it is always associated with a smaller bias. We conclude that, if the sample locations are not homogeneously scattered over the observed region, then corrector methods such as that proposed by downweight clustering, are advisable. Remember that the simulation studies included in Chapter 4 pointed to the same conclusion. Where efficiency is concerned, we have now proved that the benefit should not be disregarded. Bear in mind, however, that the proposed estimator (like the other estimators) does not consider the preferability issue. Therefore, under preferential sampling, all of them present a non-negligible bias and the efficiency issue consequently becomes less relevant. 3That corresponds to first and third rows of panels in Figure 6.5. 6.3. EFFECT ON VARIOGRAM ESTIMATION 133 Effect of white noise The last simulation study included in this Section is related to the analysis of the effect of noise. As has already been stated, the Y(.) process can be considered the noisy version of S(.). Suppose we have more noise variance τ2, but a fixed value for the total variance τ2+σ2. We intend to observe how the variogram estimation for the Y(.) process would be affected. To simplify notation, let us represent this estimator by b VY, and that associated with S(.) by b VS. The two are related through the expression b VS≡b VY−bτ2. An estimate of τ2is, however, normally difficult to obtain. We changed τ2from 0.25 to 0.81. The value for σ2was chosen in such a way that the total variance, i.e. sill, remains equal to 2.5. In Figure 6.6, the efficiency results are plotted. The same three sampling designs and the same three empirical estimators are taken into account. Besides efficiency, other performance indicators, such as the ISE, were derived. The main conclusions can be summarized as follows. Under CSR, the classic estimator is no longer a reasonable alternative to either kernel estimator. The presence of strong noise really degrades the performance of this estimator, in terms of variance and bias, for small lags. The two kernel estimators seem considerably resistant to the presence of noise. This may be justified by the adoption of a specific asymmetric boundary kernel, near the variogram endpoint 0, as briefly explained in Chapter 3 (more detail to be found in Garcia-Soid´an et al. 2004). Under non-CSR, the white noise continues to affect the classical estimator but not the kernel ones. Curiously, these two last estimators present a little less bias than in the previous case study. Finally, if the two kernel estimators are compared in terms of efficiency, one may say that the gain from using the proposed kernel estimator is almost negligible compared with that from the NW kernel. We highlight the fact that the main contribution of the former continues to be some bias reduction. 134 CHAPTER 6. ASSESSING PREFERENTIAL SAMPLING 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.4 0.8 beta=0 distance STDEV and BIAS Classic NW kernel RobClusters 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.4 0.8 beta=0 distance sqrt(MSE) Classic NW kernel RobClusters 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.4 0.8 beta=2 distance STDEV and BIAS Classic NW kernel RobClusters 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.4 0.8 beta=2 distance sqrt(MSE) Classic NW kernel RobClusters 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.4 0.8 beta=2 just clustered distance STDEV and BIAS Classic NW kernel RobClusters 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.4 0.8 beta=2 just clustered distance sqrt(MSE) Classic NW kernel RobClusters Figure 6.6: Comparison of variogram estimators through their efficiency (measure in terms of √V ar and √Mse) when sampling is clustered, in both the preferential and non-preferential sub-cases. The estimation bias is also plotted (dashed-lines). It is being considered τ2= 0.81 and σ2= 1.69. 6.4. IMPACT ON PREDICTION 135 6.4 Impact on prediction Optimal prediction, i.e kriging, typically depends on knowledge of spatial dependency, as it is supposed that some measurements in the vicinity of the point investigated should be more closely related than others to the true unknown value. It is then supported by the second-order properties of the spatial process. One point of interest to be addressed is the analysis of the influence of earlier variogram estimators on the consequences for prediction Kriging is indeed considered to provide the best linear unbiased estimator (abbreviated as BLUE) of the unknown characteristic studied. Linear because its estimates result from a weighted linear combination of sample data; unbiased because the mean of prediction errors (deviation between true value and predicted value) is expected to be null; best because the variance of the prediction errors (prediction variance) is at its minimum. For simplicity, suppose now that interest lies in the process Y(.) and not in its noiseless version S(.). Our target for prediction becomes Y(x0) the value of process Y(.) at a generic location x0, given sample data (xi, yi), i= 1,2, ..., n. The prediction problem may be formalized by invoking the conditional distribution of Y(.) given the observed data y. E[Y(x0)|y] = b Y(x0) specifies the predicted value and Var[Y(x0)|y] = E[(Y(x0)−b Y(x0))2] specifies the prediction variance (or prediction mean square error). The optimal predictor b Y(x0) minimises MSE[b Y(x0)] = E[(Y(x0)−b Y(x0))2], where the expectation is specified by the joint distribution of Y(x0) and Y. As discussed in Cressie (1993), this best predictor is not always linear in y(under non-Gaussianity). Therefore, one should additionally require b Y(x0) = n X i=1 λiyi+λ. We now aim to minimize (over coefficients λ1, ..., λn, λ) E ÃY(x0)− n X i=1 λiyi−λ!2 = Var "Y(x0)− n X i=1 λiyi#+õ(x0)− n X i=1 λiµ(xi)−λ!2 , 136 CHAPTER 6. ASSESSING PREFERENTIAL SAMPLING where µ(x) = E[Y(x)],x∈D. The derived solution is λ=µ(x0)−Pn i=1 λiµ(xi) and (λ1, ..., λn)t=ctΣ−1, where c= (C(x0,x1), ..., C(x0,xn))t,Σis a n×nmatrix with (i, j)th element C(xi,xj) and C(.) is the covariance function. The optimal linear predictor becomes b Y(x0) = ctΣ−1(y−µ) + µ(x0), where µ= (µ(x1), ..., µ(xn))t. The minimized mean-squared prediction error is Var[Y(x0)|y] = C(x0,x0)−ctΣ−1c. Such a sampling prediction technique, assuming a known mean function µ(x), was named simple kriging by Matheron (1971). If µ(x) is equal to an unknown constant µ, one is in the presence of ordinary kriging. Alternatively, universal kriging is used when µ(x) is linear in a fixed number of unknown parameters. All these kriging techniques produce optimal linear predictors. 6.4.1 Gaussian data Suppose that we disregard the possibility of data being preferentially sampled and we adopt previous standard kriging techniques. Under a Gaussian model like the one proposed in Section 6.2, we would have Y(x) = S(x) + N(0, τ2), where Y(.) is the measurement process and S(.) is the target of prediction. So, interest is now in process S(.), given the observed data y. [S(x0), Y ] is supposed to be multivariate Gaussian with mean vector (µ, µ1)t and covariance matrix  σ2σ2rt σ2rτ2I+σ2R  6.4. IMPACT ON PREDICTION 137 where ris a vector with elements ri=ρ(kx0−xik) : i= 1, ..., n,Ris a n×n matrix with (i, j)th element ρ(kxi−xjk). The 1is a n−length vector of ones and Iis the n×nidentity matrix. As described in Diggle et al. (2003), the minimum mean square error predictor becomes b S(x0) = σ2rt(τ2I+σ2R)−1(y−µ1) + µ(6.1) and with prediction variance Var[S(x0)|y] = σ2−σ2rt(τ2I+σ2R)−1σ2r(6.2) Consequently, the prediction variance depends on the correlation model, on the spatial configuration of the data and on the prediction location, but it does not directly4depend on the actual data. An extension to this Gaussian model may be taken into account by dealing with a non-constant mean value surface x. Typically, it may be useful to consider µ(x) = Σp j=1βjfj(x), where f1(x), ..., fp(x) are observed functions of location x, or functions of observed covariates, leading to universal kriging or kriging with a trend model (Wackernagel 1998). An estimator for the unknown β= (β1, ..., βp)tmay be derived by maximum likelihood, yielding b β= (FtV−1F)−1FtV−1y,(6.3) where Fis a n×pmatrix with (i, j)th element fj(xi) and V=R+τ2/σ2I. In this case, the expression (6.1) for minimum mean-squared predictor would change slightly to give b S(x0) = σ2rt(τ2I+σ2R)−1(y−Fb β) + F0b β where F0is a 1 ×pmatrix with (1, j)th element fj(x0). 4Note that one may say that the prediction variance indirectly depends on data y, through the estimation of the parameters.