The universal kriging (UK) block predictor of Z(A) is well known to be,

where S = var(Z); c(A) = (C(s₁, A),…,C(s_n, A))' and C(s, A) º (∫_A C_θ(s-u)du)/|A|; C_θ(s-u) º cov(Z(s), Z(u));  X º (x_j(s_i)) is an n ´ p matrix of explanatory variables; x(A) º ∫_A x(u)du/|A|; and  is the best linear unbiased estimator of , namely，

Now consider prediction of g(Z(A)), where g is a smooth nonlinear function and |A| > 0 or A is countable. Cressie (1993a) proposed the constrained kriging (CK) predictor , which is an optimal linear predictor that is unbiased for Z(A) (as is ) and satisfies the extra constraint, . Cressie (1993a) showed that .

In the spirit of wanting the predictor's statistical properties to match those of the target quantity, Aldworth and Cressie (2003) extended this methodology to include further constraints, where (some) covariances are matched in addition to the variance(s). This covariance-matching constrained kriging (CMCK) predictor  includes  as a special case.

Consider the problem of predicting the M-dimensional vector, Z(A) º (Z(A₁),…, Z(A_M))', where A_i Ì D; i = 1, …, M and A º {A₁,…, A_M}. The CMCK predictor, given by Aldworth and Cressie (2003), is:

where X_M º (x_j(A_i)) is an M ´ p matrix, C_M º (c(A₁),…, c(A_M)) is an n ´ M matrix, and K is an M ´ M matrix defined as follows. Write  and  Under the assumption that both of these M ´ M matrices are positive-definite, there exist nonsingular matrices P₁ and Q₁ such that P=P₁'P₁ and Q=Q₁'Q₁. Then K º Q₁^-1P₁. Notice that the two constraints,  and are satisfied; the latter constraint involves covariances as well as variances.

When the CMCK predictor is not defined (i.e., P or Q is not positive-definite), the covariance constraints can be relaxed until it is defined. Cressie and Johannesson (2001) took the approach of including the (M – 1) nearest data locations to prediction region A, and then predicting Z(A) using the implied M(M+1)/2 variance-covariance constraints.  This is what we shall do in the next section on TCDD contamination. Also, in the next section, we shall detrend the data first and apply CMCK to the residual process R(×) with constant trend, namely E(R(×)) = β₀. This amounts to putting p = 1 and x(s) º 1 in the formulas above.

 

Background

In environmental-remediation problems, soil contamination at one location often leads to contamination at other locations because of the conductivity properties of soil. This has been demonstrated by many case studies such as for the TCDD (tetrachlorodibenzo-p-dioxin) data, which were analyzed by Zirschky and Harris (1986) and Waller and Gotway (2004). The analysis we give below has been taken from Cressie, Zhang, and Craigmile (2005).

In 1971, a truck transporting dioxin-contaminated residues dumped an unknown quantity of waste in a rural area of Missouri in order to prevent citations for being overweight. Although the highest concentration occurred where the waste was dumped, contamination had spread to the shoulders of an adjacent state highway. In November 1983, The US Environmental Protection Agency (EPA) collected soil samples along transects and measured the TCDD concentration (in μg/kg) in each sample. Samples were composited along the transects. In our analysis, we consider only the data close to the dumping location, where the TCDD concentrations tend to be larger. (Recall that our goal is to estimate high extrema and their extent.) These data are based on 50-foot transects and there are no non-detects in the region.

Thus, we consider measurements of 31 samples of TCDD within a region of D, a 458 ´ 78 square-foot rectangle. The direction parallel to the highway and the transects is defined to be the x-direction, and the y-direction is then perpendicular to the highway. The x-coordinates of the data are the x-values of the start of each transect. A plot of the data and their spatial locations is shown in Figure 1(a), where the decimal point represents the (x,y) coordinate of the data.

Spatial Analysis

The TCDD data appear to be lognormally distributed; see Figure 1(b).  Let {Z(s_i) : i =1,...,31} denote the observed log concentrations. Based on some exploratory plots and regression analysis, we considered the model,  for the log-concentration process, where μ(s) = x(s)'β, β = (β₀, β₁)', and x(s) º (1, y²)' denotes a quadratic trend in the y-direction for s = (x, y). The ordinary least squares estimates of the regression coefficients are  and. Recall that the x-axis runs along the center of the highway; then we conclude that the quadratic surface of log concentrations in the y-direction is probably due to the drainage system along the highway, which is designed to let water run off the road quickly.  See Figure 1(c).

We then removed the trend component and analyzed the spatial process generated by the residuals:

                         .

Histograms of the residuals show that normality is a fairly good approximation (given the small sample size); see Figure 2(a). 

 

Fig 1. The top panel displays the spatial location and value of the log TCDD concentrations; the decimal points represent the (x, y) coordinates of the data.  The bottom left panel shows a histogram of the log TCDD concentrations.  The bottom right panel displays the estimated trend for the log TCDD concentrations. 

 

Directional empirical semivariograms of the residuals demonstrated that the residual process is severely anisotropic. To make the empirical semivariograms look isotropic, we multiplied the y-values by a factor of 7 (and the x-values were left unchanged). With this transformation, the data almost lie on a square coordinate system (the transformed y-values range between [-266, 280] and the x-values range between [-230, 228]).  We then considered various semivariogram models for the residual process. Using the weighted-least-squares method (Cressie, 1993b, p. 99), we deemed that the spherical model fitted best.

 

The general form for the covariance function derived from the isotropic spherical model is,

where h denotes the Euclidean distance on the transformed coordinate system and the vector of spatial parameters is θ = (c₀, s², r)'. The semivariogram is given by,

.

The estimate of the nugget effect, c₀, was 0; the estimate of the range, r, was 91.61; and the estimate of the partial sill, s², was 0.36. Figure 2(b) shows the estimated semivariogram . Notice that the estimated value  implies a stronger spatial dependence than is apparent from the empirical semivariogram values, also shown in Figure 2(b).  The estimated covariance function is

 

Fig 2.  The left panel displays a histogram of the residuals.  The right panel shows the empirical semivariogram for the residuals; the solid line is the estimated spherical semivariogram .

Comparison of prediction methods

Based on the spatial analysis we carried out in the previous section, we now consider prediction of the log-concentration process. Our aim is to estimate the value and extent of high extremes of the log TCDD concentration process in a region B. For this analysis, we let B º D, the rectangular region that encloses all the observed sites. We shall predict on a 2 ´ 2 ft. grid of points throughout D. We denote this discrete approximation to D by {s_j*: j =1,..., m = 9200}.

For the CMCK predictor, we first predicted the residual process R(×) using CMCK and then we added the estimated trend component, , to obtain predictions of the log-concentration process, Z(×). The exceedance set for the process Z(×) in the region B associated with an absolute threshold t Î IR is defined by , which we estimate by  for some generic predictor ; see Craigmile et al. (2005). The two panels of Figure 3 display the exceedance sets obtained when we use ordinary kriging plus trend and CMCK plus trend with M = 2 (i.e., one neighbor and constraining the variance and the covariance between nearest neighbors). In each case, we let the threshold be the 90^th percentile of the inverse averaged Gaussian cumulative distribution function (averaged over B). The spatial patterns of exceedances show two regions of exceedance at around x » 100 feet and CMCK shows another region of exceedance at around x » -100 feet. Ordinary kriging estimates 3.16% of the prediction locations are in the exceedance set and CMCK estimates 17.09% of the prediction locations are in the exceedance set.

 

Fig 3.  The gray regions denote those locations with predicted values that exceed the 90th percentile of the inverse ACDF for ordinary kriging plus trend and CMCK plus trend.