Bootstrap based Trans-Gaussian Kriging

Abstract

Trans-Gaussian Kriging is a popular method for predicting unobserved values of a nonstationary and non-Gaussian spatial process. However, the predictions generated by Trans-Gaussian Kriging are often biased, primarily because these are produced by applying nonlinear transformations to unbiased optimal predictors of the underlying Gaussian field. The existing approach to bias correction is based on analytical approximations that only alleviate the bias problem partially. In this paper, we formulate a bootstrap method for bias correction and show that under some conditions, it yields asymptotically unbiased predictions. Results from a moderately large simulation show that the proposed method works well in reducing the bias in finite samples. A real data example is also presented illustrating the methodology.

Keywords

asymptotic unbiasedness bias correction log-normal Kriging non-stationarity spatial prediction

1 Introduction

Nonstationary and non-Gaussian spatial data appear frequently in practice. Although there is well established theory of unbiased prediction for intrinsically stationary spatial processes (see, for example, Cressie, 1991; Stein, 1999), the literature on unbiased prediction is parse for non-Gaussian nonstationary spatial processes. The most common approach to predicting a non-Gaussian spatial process at an unobserved location is Trans-Gaussian Kriging, where the observations are first transformed to a Gaussian process through a (known) link function and then standard theory of unbiased prediction or Kriging is applied to the Gaussian process. The resulting predicted value is then transformed back to the original scale by using the inverse of the link function. This yields a reasonable prediction of the original spatial process. However, the resulting naive predictor is biased, simply due to the fact that in general, E(φ(X)) ≠ φ(E(X)) for a nonlinear function φ.

The reasons for opting for unbiased estimators and predictors are well documented in the literature (cf. Lehmann and Casella, 1998). A statistical procedure with a nontrivial bias produces a distorted picture of the population and may lead to very unreasonable conclusions. To eliminate the bias, Cressie (1991) suggested an additive correction factor based on Taylor’s series expansion and some analytical considerations. Here we formulate a new methodology for bias correction based on the bootstrap that does not require any analytical consideration on the part of the user. The proposed methodology generates bias-corrected predictors essentially using an algorithm that can be implemented quite generally and that is applicable to a large class of nonlinear transformations. However, the method involves computation and provides a way to bypass the required analytical derivations, which are often intractable.

The main results of the paper address the bias properties of the existing as well as new bias-corrected predictors. In Section 4, we show that the bias-corrected predictor based on the Taylor’s series expansion eliminates the bias only partially, except for some very special types of transformations, namely, for low order polynomials. In comparison, we show that under some regularity conditions, the proposed bootstrap based method removes the bias of the naive predictor entirely in the ideal case where some of the unknown process parameters are known, and asymptotically when the unknown parameter are replaced by their consistent estimators. This is an important advancement over the existing solutions in the literature as these do not produce an asymptotically unbiased predictor of the non-Gaussian process in general. Simulation results indicate that the proposed bias corrections do not inflate the mean square prediction errors (MSPEs) by a significant amount. Indeed, the bootstrap based bias-corrected predictors have comparable (and sometimes smaller) MSPEs than the existing predictors.

The layout of the paper is as follows. In Section 2, we briefly describe the Trans-Gaussian Kriging method and present the bias-corrected predictor based on Taylor’s expansions. In Section 3, we describe the bootstrap based bias-corrected predictors—both when the process parameters are known and when these are unknown. Section 4 states the theoretical properties of the bias-corrected estimators. Sections 5 and 6 respectively, present results from a moderately large simulation study and an illustrative example involving real data from Oklahoma on the prediction of Heating Degree Days (HDDs). Proofs of the technical results are given in Section 7.

2 Background

2.1 Basic framework

Let Z(⋅) be a spatial process on a domain D ⊂ ^d and let Z ≡ {Z(s_i): i  1, …, n} denote the observations, where the sites s₁, …, s_n either lie on a regular grid (in which case, we take D ≡ ^d, the d-dimensional integer grid) or are irregularly spaced (in which case D ≡ ^d). Furthermore, let s₀ ∈ D be a site such that Z(s₀) is unobserved. We shall suppose that the observable process Z(⋅) is such that for a known invertible link function φ(⋅), ZφY

with Y(2.1)µβD

where µ(s; β) is a non-random function and where {(s): s ∈ D} is a zero-mean, second order stationary Gaussian process with variogram 2γ_{(⋅; θ), ≡ E((⋅) – (0))² and covariogram C(⋅; θ) ≡ E(⋅) (0). Here we suppose that µ(⋅; β) and 2γ(⋅; θ) (and hence C(⋅; θ) are known, except for the finite dimensional mean parameters β ∈ Band covariance parametersθ∈ Θ. We now give two common examples of such Trans-Gaussian spatial processes.}

Example 1. Log-Normal Kriging:

Suppose that the Z(⋅) variables are positive with probability one and let YZD

be a Gaussian process as in (2.1). Then Z(⋅) is a Log-Normal process.

Example 2. Logit Kriging:

Let Z(s) denote a (0, 1) valued random variable, e.g., a proportion, such as the mortality rate due to a certain disease. Let YZD,

where logit z  log (z/[1–z]), z ∈ (0, 1) and where Y(s) is a Gaussian process as in (2.1). Then Z(⋅) is a Logit-Normal process.

2.2 Trans-Gaussian Kriging

Next we briefly describe the Trans-Gaussian Kriging method (cf. Cressie, 1991). The basic idea is to use the inverse link function to transform the observed spatial process to a Gaussian process and then detrend the data and perform Kriging on the residuals. To fix ideas, first consider the ‘ideal’ case where the parameters β andθare known. In this case, the exact mean of Y(⋅)’s are known and hence, we may work with the Gaussian error variables (s)  Y(s) – µ(s; β) directly. Let $\tilde{Y}$ (S₀) be the Ordinary Kriging predictor (OKP) of Y (s₀) based on Y ≡ (Y(s₁), …, Y(s_n)′, where A′ denotes the transpose of a matrix A. Specifically, if  ≡ ((s₁), …, (s_n))′, then the OKP is

\tilde{Y} (s_{0}) = μ (s_{0}; β) + λ^{'} ϵ = μ (s_{0}; β) + \sum_{i = 1}^{n} λ_{i} ϵ (s_{i}),

(2.2)

where λ′  (λ₁, …, λ_n) are the Ordinary Kriging weights given by (cf. Cressie, 1991 p. 123)

λ^{'} \equiv {(γ + 1^{'} \frac{1 - 1^{'} Σ^{- 1} γ}{1^{'} Σ^{- 1} 1})}^{'} Σ^{- 1} .

(2.3)

where γ ≡ (C(s₀ – s₁; θ), …, C(s₀ – s_n; θ))′ and Σ be the n × n matrix with (i, j)th entry C(s_i – s_j; θ) and 1 is the n ×1 vector of 1 s. The ‘natural’ but biased Trans-Gaussian Kriging predictor of Z(s₀) in the ‘ideal’ case is then given by

\tilde{Z} (s_{0}) = ϕ (\tilde{Y} (s_{0})) .

(2.4)

In practice, the parameters β andθare typically unknown and hence, the error variables  (s_i)’s are not observable. As a result, the OKP in (2.2) cannot be directly used in practice. Although β andθare unknown, we shall suppose that estimators $\hat{β}$ and $\hat{θ}$ of these model parameters are given and instead of  (s), we will work with the residuals

ϵ_{i} (s) = Y (s_{i}) - μ (s_{i}; \hat{β}), i = 1, ..., n .

Also, let $\hat{γ}$ and $\hat{Σ}$ be estimated versions of γ and Σ obtained by replacingθwith $\hat{θ} .$ Define the estimated OKP of Y(s₀) by

\hat{Y} (s_{0}) = μ (s_{0}; \hat{β}) + \sum_{i = 1}^{n} {\hat{λ}}_{i} \hat{ϵ} (s_{i}),

(2.5)

where ${\hat{λ}}^{'} = ({\hat{λ}}_{1}, ..., {\hat{λ}}_{n})$ is given by (2.3) with $γ = \hat{γ}, Σ = \hat{Σ} .$ The naive Trans-Gaussian predictor of Z(s₀) is then given by

\hat{Z} (s_{0}) = ϕ (\hat{Y} (s_{0})) .

(2.6)

2.3 Bias-corrected predictors

As pointed out by a referee, the predictor in (2.4) is the (conditional) ‘maximum likelihood predictor’ of Z(s₀) in the Gaussian setting. However, both the predictors in (2.4) and (2.6) are biased for a nonlinear φ(⋅) such as in Examples 1 and 2. Often, this bias is significant when Y(s₀) is in the tails of its marginal univariate normal distribution. Then, depending on the extremity of the transformation performed by φ(⋅), the bias of $\hat{Z} (s_{0})$ can be pronounced. There are many important applications, particularly in governmental policy formulation and resource allocation, where a nontrivial prediction bias has serious practical implications. In such situations, predictors with a small bias are preferred (cf. Remark 5.1). As a result, often a bias-corrected version of $\hat{Z} (s_{0})$ is constructed. The following bias correction based on the delta method is commonly used (cf. Cressie, 1991, p. 137):

{\hat{Z}}^{C} (s_{0}) = \hat{Z} (s_{0}) + ϕ^{″} ({\hat{μ}}_{0, Y}) (\frac{{\hat{τ}}^{2} (s_{0})}{2} - {\hat{m}}_{Y}),

(2.7)

where ${\hat{μ}}_{0, Y} = μ (s_{0}; \hat{β})$ is an estimate of E(Y(s₀)), ${\hat{τ}}^{2} (s_{0})$ is an estimate of the Kriging variance, or MSPE, of $\tilde{Y} (s_{0}),$ given by

{\hat{τ}}^{2} (s_{0}) = C (o; \hat{θ}) - {\hat{λ}}^{'} \hat{γ} + {\hat{m}}_{Y},

(2.8)

and where ${\hat{m}}_{Y} = m_{Y} (\hat{θ})$ is the estimated Lagrangian multiplier, given by $m_{Y} (θ) = \frac{(1 - 1^{'} Σ^{- 1} γ)}{1 Σ^{- 1} 1}$ . It can be shown (cf. Rister, 2010, Chapter II) that the correction factor in (2.7) is derived from the formal Taylor expansion of φ(⋅), centered at EY(s_o),

ϕ (Y (s_{0})) = \sum_{i} \frac{ϕ^{(i)} (μ (s_{0}; β))}{i!} {[Y (s_{0}) - μ (s_{0}; β)]}^{i},

(2.9)

where φ⁽ⁱ⁾(⋅) denotes the ith derivative of φ(⋅). Notice that the approximate predictor in (2.7) ignores the terms in this summation for all i ≥ 3. While this yields a bias-corrected predictor for low order polynomials, such as φ(y)  y³ (with the higher order derivatives zero), this fails to correct the bias of most nonlinear functions with nonzero higher order derivatives. In the next section, we construct a bias-corrected predictor based on a spatial bootstrap method that produces asymptotically unbiased predictors for a large class of nonlinear functions, thereby improving upon the bias correction method described above.

3 Bootstrap based bias-corrected prediction

There are many ways of adjusting for the bias of the estimated Kriging predictor $\hat{Z} (s_{0})$ in (2.6) using the bootstrap (cf. Hall and Maiti, 2006; Lahiri et al., 2007). Here, we will first use a multiplicative correction factor c such that

{\hat{Z}}^{B C}_{M} (s_{0}) \equiv c \hat{Z} (s_{0})

(3.1)

is approximately unbiased. Here, the subscript M corresponds to a ‘multiplicative’ bias correction. The multiplicative bias correction factor c described in (3.1) works the best for nonnegative Z-variables, such as in the case of log-normal Kriging. For Z-variables with an unrestricted range, such as in power transformation models, an alternative is to use an additive correction factor a based on the bootstrap such that

{\hat{Z}}^{B C}_{A} (s_{0}) \equiv \hat{Z} (s_{0}) + a

(3.2)

is approximately unbiased. To motivate the main steps, we first describe the bootstrap methodology for the simpler ‘ideal’ case where the parameter values β andθare known, in Section 3.1, and then extend it to the unknown parameter case in Section 3.2.

3.1 Known parameters

For now, we will work in the setting where we assume that we know the values of the covariance parameterθand the mean parameter β. Note that in the multiplicative case that the factor that makes $\hat{Z} (s_{0})$ exactly unbiased is given by

c \equiv c (β, θ) = E Z (s_{0}) / E \tilde{Z} (s_{0}),

while in the additive case,

a \equiv a (β, θ) = E Z (s_{0}) - E \hat{Z} (s_{0}) .

The explicit forms of c(β, θ) and a(β, θ) are difficult to derive. In this section, we present a bootstrap method for estimating the factors c(β, θ) and a(β, θ) that does not require the user to do such analytical work; it replaces the tedious and intractable analytical derivation by a simple, albeit computer intensive bootstrap algorithm.

The bootstrap algorithm for the known parameter case for estimating a and c is as follows:

Generate *(s₀), *(s₁), …, *(s_n), from the Gaussian process with covariogram C(⋅; θ). Note that this involves the generation of n + 1 observations.

Compute $\tilde{Y} * (s_{0})$ by replacing { (s_i): i  1, …, n} in (2.2) with the generated variables {*(s_i): i  1, …, n}.

Compute $\tilde{Z} * (s_{0}) = ϕ (\tilde{Y} * (s_{0}))$ and $Z * (s_{0}) = ϕ (Y * (s_{0})) .$

Let ${\tilde{a}}_{n} = E * Z * (s_{0}) - E * \tilde{Z} * (s_{0})$ and ${\tilde{c}}_{n} = E * Z * (s_{0}) / E * \tilde{Z} * (s_{0})$ , where E* denotes the conditional expectation given Υ. Then, the bootstrap-based bias-corrected predictors of Z(s₀) are given by

{\tilde{Z}}_{A}^{B C} (s_{0}) = \tilde{Z} (s_{0}) + {\tilde{a}}_{n}

{\tilde{Z}}_{M}^{B C} (s_{0}) = {\tilde{c}}_{n} \cdot \tilde{Z} (s_{0}) .

In practice, ${\tilde{a}}_{n}$ and ${\tilde{c}}_{n}$ are evaluated by Monte Carlo simulation, where steps 1, 2 and 3 are repeated a large number of times, say B times, resulting in the bootstrap replicates ${\tilde{Z} *^{j} (s_{0}) : j = 1, ..., B}$ and ${Z *^{j} (s_{0}) : j = 1, ..., B},$ for $\tilde{Z} * (s_{0})$ and $Z * (s_{0}),$ respectively. Then, the Monte-Carlo approximations to and ${\tilde{c}}_{n}$ are, respectively, given by

{\tilde{a}}_{n}^{M C} \equiv \frac{1}{B} [\sum_{j = 1}^{B} Z *^{j} (s_{0}) - \sum_{j = 1}^{B} \tilde{Z} *^{j} (s_{0})]

{\tilde{c}}_{n}^{M C} \equiv \frac{\sum_{j = 1}^{B} Z *^{j} (s_{0})}{\sum_{j = 1}^{B} \tilde{Z} *^{j} (s_{0})} .

3.2 Unknown parameters

In real life applications, the true values of the parameters β andθare unknown and hence, the aforementioned predictors must be calculated using their estimates. We suppose that estimators $(\hat{β}, \hat{θ})$ are given, such as the (Gaussian) maximum likelihood estimator or weighted least squares estimators of (β, θ) (cf. Cressie, 1991). The following are the main steps of the bootstrap algorithm in the unknown parameter case:

Generate *(s₀), *(s₁), …, *(s_n) from the Gaussian process with estimated covariogram $C (\cdot; \hat{θ}) .$

For i  0, 1, …, n, define the bootstrap variables $Z * (s_{i}) = ϕ (Y * (s_{i})) w i t h Y * (s_{i}) = μ (s_{i}, \hat{β}) + ϵ * (s_{i}) .$

Use the ‘bootstrap observations’ {Z*(s_i): i  1, …, n} in place of {Z(s_i): i  1, …, n} and the same estimation procedure to find the bootstrap version $(\hat{β} *, \hat{θ} *)$ of the estimators $(\hat{β}, \hat{θ})$ .

Compute $\hat{γ} *$ and $\hat{Σ} *$ , the bootstrap versions of $\hat{γ}$ and $\hat{Σ}$ obtained by replacing $\hat{θ}$ with $\hat{θ} * .$

Next, define the bootstrap version of the estimated OKP $\hat{Y} (s_{0})$ by

{\hat{Y}}^{*} (s_{0}) = μ (s_{0}; \hat{β} *) + \sum_{i = 1}^{n} {\hat{λ}}_{i} * \hat{ϵ} * (s_{i}),

(3.3)

where $\hat{λ} * = ({\hat{λ}}_{1} *, ..., {\hat{λ}}_{n} *)$ is given by (2.3) with $γ = \hat{γ} *, Σ = \hat{Σ} *$ and where $\hat{ϵ} * (s_{i}) = Y * (s_{i}) - μ (s_{i}, \hat{β} *)$ , i  1, …, n.

Set $\hat{Z} * (s_{0}) = ϕ (\hat{Y} * (s_{0}))$ and $Z * (s_{0}) = ϕ (Y * (s_{0})) .$ Define ${\hat{a}}_{n} = E * Z * (s_{0}) - E *$ $\hat{Z} * (s_{0})$ and ${\hat{c}}_{n} = E * Z * (s_{0}) / E * \hat{Z} * (s_{0}),$ where E* denotes the conditional expectation given Υ.

Then, the bootstrap-based bias-corrected predictors of $Z (s_{0})$ are given by

{\hat{Z}}_{A}^{B C} (s_{0}) = \hat{Z} (s_{0}) + {\hat{a}}_{n}

(3.4)

{\hat{Z}}_{M}^{B C} (s_{0}) = {\hat{c}}_{n} \cdot (s_{0}),

(3.5)

respectively. Again, a Monte Carlo simulation is usually used to approximate ${\hat{a}}_{n}$ and ${\hat{c}}_{n} .$ Steps 1–6 are repeated B times resulting in the bootstrap replicates ${\hat{Z} *^{j} (s_{0}) : j = 1, ..., B} a n d {Z *^{j} (s_{0}) : j = 1, ..., B} .$ Then, the Monte Carlo approximation to ${\hat{c}}_{n}$ is given by

{\hat{c}}_{n}^{M C} \equiv \frac{\sum_{j = 1}^{B} Z *^{j} (s_{0})}{\sum_{j = 1}^{B} \hat{Z} *^{j} (s_{0})} .

A similar construction applies to ${\hat{a}}_{n}^{M C}$ in the unknown parameter case.

4 Theoretical results

In this section, as well as in the accompanying proofs, we highlight the dependence of various quantities depending on n. Thus, we shall write $\hat{θ} = {\hat{θ}}_{n},$ $\hat{β} = {\hat{β}}_{n},$ γ  γ_n,_{Σ  Σ_n,
$\tilde{Z} (s_{0}) = {\tilde{Z}}_{n} (s_{0}),$

$\hat{Z} (s_{0}) = {\hat{Z}}_{n} (s_{0}),$
etc. Although we do not impose explicit conditions on the data locations s₁, …,s_n the regularity conditions used in here are more suitable for the increasing domain scenario. See Remark 4.1 below for more details.}

4.1 Bias of the existing predictors

The main results of this section give explicit expressions for the bias of the naive predictor $\hat{Z} (s_{0})$ and its first-order Taylor series based corrected version ${\hat{Z}}^{C} (s_{0}) .$ Note that both of these predictors involve the estimators ${\hat{β}}_{n}$ and ${\hat{θ}}_{n}$ of the unknown parameters. The effects of these on the biases of the predictors $\hat{Z} (s_{0})$ and ${\hat{Z}}^{C} (s_{0})$ are propagated through the mean function µ(⋅; β) and the covariogram C(⋅; θ). When the mean and the covariance functions of the error process { (s)} are smooth and the estimators ${\hat{β}}_{n}$ and ${\hat{θ}}_{n}$ are consistent, the effects of substituting the estimators for the true unknown values on the final bias of the predictors become asymptotically negligible. We shall use the following conditions to study the bias properties of $\hat{Z} (s_{0})$ and ${\hat{Z}}^{C} (s_{0})$ To state the conditions, let Γ(⋅) denote the Gamma function, let $σ_{0}^{2} = C (0; θ)$ and let $σ_{n}^{2} = E \tilde{ϵ} {(s_{0})}^{2} = λ^{'}_{n} Σ_{n} λ^{'}_{n} .$

Conditions:

(C.1)

The function φ(⋅) admits a power series representation given by

ϕ (x) = \sum_{k = 0}^{\infty} d_{k} {(x - μ (s_{0}, β))}^{k}, x \in ℝ,

(4.1)

for some $d_{0}, d_{1}, . . . \in ℝ$ Note that d_k’s depend on µ(s₀, β).

$E {[{\hat{Z}}_{n} (s_{0})]}^{2} = O (1)$ as $n \to \infty$ .

For some K₁ ∈(0, ∞),

\sum_{k = 1}^{\infty} \sum_{j = 1}^{\infty} k j | d_{k} d_{j} | 2^{3 (k + j - 2) / 2} Γ (\frac{k + j - 1}{2}) [σ_{0}^{j + k - 2}] < K_{1.}

(4.2)

(C.2) As n → ∞,

| μ (s_{0}; {\hat{β}}_{n}) - μ (s_{0}; β) | + | \hat{ϵ} (s_{0}) - \tilde{ϵ} (s_{0}) | \to_{p} 0.

(4.3)

Next we comment about the regularity conditions. Condition (4.1) on φ(⋅) implies that φ(⋅) is infinitely diffentiable, although a more complicated expression can be derived for functions φ(⋅) that are only twice differentiable. The power series representation is specifically used for deriving an explicit expression for the bias of the naive predictor in terms of $σ_{0}^{2}$ and $σ_{n}^{2}$ . Note that this condition holds for both Log-normal Kriging and Logit Kriging (Examples 1 and 2). The requirement that $E {[{\hat{Z}}_{n} (s_{0})]}^{2} = O (1)$ can be weakened to $E | {\hat{Z}}_{n} (s_{0}) |^{1 + δ} = O (1)$ for some δ  0. Here, we take δ  1 as it seems to be the most natural condition in this context; the MSPE of ${\hat{Z}}_{n} (s_{0})$ is bounded only when $E {[{\hat{Z}}_{n} (s_{0})]}^{2} = O (1) .$

Condition (4.2) is a sufficient condition for deriving the series representation of Bias $({\hat{Z}}_{n} (s_{0})) .$ Under (4.2), the sum on the right side of (4.4) is absolutely convergent.

Finally, Condition (4.3) is a mild requirement on the estimators ${\hat{β}}_{n}$ and ${\hat{θ}}_{n} .$ In particular, if µ(s₀; ) is continuous at β and ${{\hat{β}}_{n}}_{n \geq 1}$ is consistent for β, then

μ (s_{0}; {\hat{β}}_{n}) - μ (s_{0}; β) \to_{p} 0.

Similarly, if $∥ Σ_{n} {(θ)}^{- 1} ∥ = O (1),$ then it can be shown that consistency of ${{\hat{θ}}_{n}}_{n \geq 1}$ and (equi-)continuity of C(s; ⋅) atθimplies that $\hat{ϵ} (s_{0}) - {\tilde{ϵ}}_{n} (s_{0}) \to_{p} 0.$

The following result gives the leading terms in the expansions for the biases of the naive predictor $\hat{Z} (s_{0})$ and its Taylor series based corrected version. ${\hat{Z}}^{C} (s_{0}) .$

Proposition 4.1 Suppose that Conditions (C.1) and (C.2) hold. Then:

B i a s ({\hat{Z}}_{n} (s_{0}) = \sum_{k = 1}^{\infty} \frac{ϕ^{(2 k)} (μ (s_{0}; β))}{k! 2^{k}} (σ_{n}^{2 k} - σ_{0}^{2 k}) + o (1),

(4.4)

B i a s ({\hat{Z}}^{C}_{n} (s_{0}) = \sum_{k = 2}^{\infty} \frac{ϕ^{(2 k)} (μ (s_{0}; β))}{k! 2^{k}} (σ_{n}^{2 k} - σ_{0}^{2 k}) + o (1) .

(4.5)

Note that under Conditions (C.1) and (C.2), the bias of the naive predictor is the difference between the values of an analytic function of $σ_{n}^{2}$ and $σ_{0}^{2} .$ For pure increasing domain asymptotic structure (cf. Cressie, 1991; Lahiri, 2003) with either regularly spaced or irregularly spaced data points, {σ_n}_n≥1 does not converge to σ₀ and therefore, the bias typically remains bounded away from zero. However, in the presence of an infill component in the spatial sampling design, under mixed or pure infill asymptotics, {σ_n}_n≥1 may converge to σ₀ (e.g., if s₀ is a limit point of s_i’s) and the bias must go to zero in this case. A similar remark applies to the terms for K ≥ 2 for the bias-corrected predictor ${\hat{Z}}_{n}^{C} {(s}_{0}) .$

4.2 Bias of the bootstrap based predictors

We now consider the bias properties of the bootstrap based bias-corrected predictors ${\hat{Z}}_{M, n}^{B C} (s_{0})$ and ${\hat{Z}}_{A, n}^{B C} (s_{0})$ . Let f_o(⋅) denote the probability density function of the N(0, 1) distribution, namely, $f_{0} (x) = \frac{1}{\sqrt{2 π}} e^{- x^{2} / 2},$ for x∈. Let g(b, u) ≡ ∫φ(µ(s₀; b) + $\sqrt{u}$ ⋅x) f₀(x)dx, b ∈ B, u ∈ [0, ∞). Also, let λ  λ_n (θ)  (λ_1n(θ), …, λ_nn(θ))’. In the following, K, K_i, K (⋅) denote constants with values in (0, ∞) that may depend on their arguments (if any), but not on n. For δ  0, write $N (δ) = {(b, t) : ∥ b - β ∥^{2} + ∥ t - θ ∥^{2} \leq δ}$ for the δ neighbourhood of the true parameter value(β, θ).

We shall make use of the following conditions for studying the bias properties of the bootstrap based predictors.

(C.3) There exist δ₀  0 and {α_i}_i≥1 ⊂ [0, ∞) with $α_{\infty} \equiv \sum_{i = 1}^{\infty} α_{i} < \infty$ satisfying the following:

|λ_in(t)| ≤ α_i for all i  1, …, n and t ∈ Θ, and $| λ_{i n} (t) - λ_{i n} (θ) | \leq K_{0} ∥ t - θ ∥ α_{i}$ for 1 ≤ i ≤ n and $∥ t - θ ∥ < δ_{0} .$

{τ_n_{(⋅)}_n≥1 is equicontinuous at the true value θ.}

$| μ (s; b_{1}) - μ (s; b_{2}) | \leq K_{1} ∥ b_{1} - b_{2} ∥$ for all s, b₁, b₂.

φ(⋅) is differentiable and there exists a monotone function Φ⁽¹⁾: [0, ∞) → [0, ∞) such that |φ⁽¹⁾ (x)|≤ Φ⁽¹⁾ (|x|) for all x ∈ , and for n sufficiently large,

E_{b, t} {[Φ^{(1)} (| μ (s_{0}; β) | + 3 K_{0} \sum_{i = 1}^{n} α_{i} | ϵ (s_{i}) | + K_{4} ‖ {\hat{β}}_{n} - β ‖)]}^{2} \leq K_{2}

for all (b, t) ∈ N (δ₀) where K₄  (1 + α_∞) K₁.

There exists a sequence {b_n}_n≥1 with b_n → 0 as n → ∞ such that for all (b, t) ∈ N(δ₀) and n ≥ 1,

E_{b, t} ∥ {\hat{β}}_{n} - β ∥^{4} \leq b_{n}^{4}

P_{b, t} (∥ {\hat{θ}}_{n} - θ ∥ > δ) \leq K (δ) b_{n}

for any δ  0.

C(0, t)  K₃ for all t and C(0, ) is continuous atθwith C(0,θ)  0.

${g^{2} {\hat{β}}_{n}, C (0, {\hat{θ}}_{n}))}_{n \geq 1}$ and ${g^{2} ({\hat{β}}_{n}, τ_{n}^{2} ({\hat{θ}}_{n}))}_{n \geq 1}$ are uniformly integrable.

$E | {\hat{Z}}_{n} (s_{0}) / E * (\hat{Z} * (s_{0})) |^{2} + {[E {\hat{Z}}_{n} (s_{0})]}^{- 1} = O (1)$ as n → ∞.

Now we briefly comment on the conditions. Condition (C.3)(i) is a smoothness condition on the Kriging weights. These can be guaranteed by assuming that the covariance matrix Σ_n(θ) and its inverse are bounded in the spectral norm and that the covariogram decays quickly and is Hölder continuous at the true value θ. See Bandyopadhyay and Lahiri (2010), where a similar condition has been used in the time series context. In the same spirit, Conditions (C.3)(ii) and (iii) are smoothness conditions on the prediction variance τ(⋅) and the mean function µ(⋅, ⋅). Condition (C.3)(iv) is used for proving asymptotic negligibility of second order terms in the Taylor’s expansion in the L¹-metric, uniformly in a neighbourhood of the true parameter value. Condition (C.3)(v) says that the estimators ${\hat{β}}_{n}$ and ${\hat{θ}}_{n}$ converge to the respective true values, locally uniformly. This can be ensured by imposing mild regularity condition on the estimating equations defining the estimators ${\hat{β}}_{n}$ and ${\hat{θ}}_{n} .$ See Pfangagl, 2000 and the references therein. The continuity of the functions in (C.3) and (C.4) at the true value (β, θ) is used for proving consistency of the bootstrap correction factors. Note that the parametric bootstrap resamples that we use here to construct the correction factors are based on the estimated parameter values and hence, we need a form of continuity of the relevant functions to guarantee consistency. Finally, Condition (C.4)(iii) is used for proving asymptotic unbiasedness of the multiplicative bias-corrected predictor ${\hat{Z}}_{M, n}^{B C}$ only—it is not needed for ${\hat{Z}}_{A, n}^{B C}$ .

The main result of this section is the following:

Proposition 4.2. Suppose that Conditions (C.3) and (C.4) hold. Then

B i a s ({\hat{Z}}_{M, n}^{B C} (s_{0})) = o (1) B i a s ({\hat{Z}}_{A, n}^{B C} (s_{0})) = o (1) .

(4.6)

Note that by construction, the ‘ideal’ bootstrap based bias-corrected predictors are exactly unbiased (for any given sample size n). The error term o(1) in (4.6) results from estimating the unknown parameters β andθby ${\hat{β}}_{n}$ and ${\hat{θ}}_{n}$ , respectively. It follows from Proposition 4.2 that Conditions (C.3) and (C.4) are sufficient to guarantee that the bootstrap bias-corrected predictors are asymptotically unbiased. Thus, unlike the naive predictor ${\hat{Z}}_{n} (s_{0})$ or its analytical bias-corrected version ${\hat{Z}}_{n}^{C} (s_{0}),$ which typically have a non-trivial bias for a large class of nonlinear transformations, the bootstrap based bias-corrected predictors’ asymptotic unbiasedness is solely under (C.3) and (C.4).

Remark 4.1: Asymptotic unbiasedness of ${\hat{Z}}_{\cdot, n}^{B C} (s_{0})$ holds for any limiting configuration of the data sites s₁, …, s_n, as long as (C.3) and (C.4) hold. In particular, under pure increasing domain asymptotic structure, $σ_{n} \to σ_{0}$ and in view of (4.4) and (4.5), the naive predictor ${\hat{Z}}_{n} (s_{0})$ or its analytical bias-corrected version ${\hat{Z}}_{n}^{C} (s_{0})$ typically has a nonvanishing bias in the limit, but ${\hat{Z}}_{\cdot, n}^{B C} (s_{0})$ do not suffer from this drawback. On the other hand, the inverse of high dimensional estimated covariance matrix (cf. (2.3) and (2.5)) becomes ill-conditioned under infill asymptotics. As a result, Condition (C.4)(iii) may fail, leading to a worse performance of the proposed version of the bootstrap based bias-corrected predictors.

In the next section, we compare the finite sample performance of the different bias-corrected predictors, both in terms of their biases as well as in terms of their MSPEs.

5 Simulations

5.1 Framework

For the purposes of these simulations, we set D to be a suitable subset of $ℝ^{2}$ with s₀ near its center. Specifically, we set s₀ to be the origin and take D to be one of the following subsets of $ℝ^{2}$ :

$D_{1} : = {(i, j) \in ℤ^{2} : | i |, | j | \leq 3}$

$D_{2} : = {(i, j) \in ℝ^{2} : | i |, | j | \leq 3 and 2 i, 2 j \in ℤ}$

$D_{3} : = {(i, j) \in ℤ^{2} : - 4 \leq i, j \leq 5}$

$D_{4} : = {(i, j) \in ℤ^{2} : - 4 \leq i, j \leq 5 and 2 i, 2 j \in ℤ}$

$D_{5} : = {(i, j) \in ℤ^{2} : - 9 \leq i, j \leq 10} .$

These choices for D allow for the observance of the asymptotic behaviour of the predictor for both expanding and filling in the observed data sites. Notice, for example, that as we move from D₁ to D₃ to D₅, the window expands from a 7 × 7 grid to a 20 × 20 grid. Also, D₂ is a filled in version of D₁, and D₄ is a filled in version of D₂, allowing for closer neighbours to be used for prediction.

We then generated (s₀), …,(s_n) from a zero mean, second order stationary Gaussian process with one of the following covariograms:

Matérn: $C (h) : = \exp (- θ ∥ h ∥)$ with θ 0.5

Exponential: $C (h) : = \exp (- θ_{1} | h_{1} | - θ_{2} | h_{2} |)$ with θ  (0.5, 1).

The first covariogram is a special case of the Matérn (1960) class of covariograms with variance equal to 1 and a smoothness parameter of 0.5, while the second is an exponential covariogram. These two covariograms give examples of both isotropic and anisotropic models.

Next, we let µ(⋅; β) ≡ 1. Then, Yµβ)

Because log-Normal Kriging and the Box-Cox Transformation are frequently used in practice for Trans-Gaussian data, here we consider the following functions to transform to Y(s₀), …,Y(s_n) to Z(s₀), …, Z(s_n):

φ(y):  y³

φ(y):  e^y

φ(y):  e^2y.

Since the mean is constant for all s, we simply used the sample mean $\bar{Y}$ to estimate it. We then calculated the estimated residuals and performed Trans-Gaussian Kriging as described in Sections 2 and 3. All simulations were done for 1000 randomly generated datasets and B  500 bootstrap replicates for each sample.

We compared the four predictors: (i) $\hat{Z} (s_{0}),$ the naive predictor, (ii) ${\hat{Z}}^{C} (s_{0}),$ the bias-corrected predictor with an additive correction factor based on the second derivative of φ(⋅), and two bootstrap predictors, (iii) ${\hat{Z}}_{A}^{B C} (s_{0})$ with an additive correction factor and (iv) ${\hat{Z}}_{M}^{B C} (s_{0})$ with a multiplicative correction factor. The predictors were compared using the following criteria:

Ratio of expected values of an observation and its prediction:

R a t i o = \frac{E Z (s_{0})}{E \hat{Z} (s_{0})} .

Bias of the predictor:

B i a s = E [\hat{Z} (s_{0}) - Z (s_{0})] .

Mean square prediction error (MSPE):

M S P E = E {[\hat{Z} (s_{0}) - Z (s_{0})]}^{2} .

An accurate predictor will have a ratio of expect values close to 1, and a bias close to 0. A smaller MSPE represents less variability in the prediction.

5.2 Simulation results

Overall, the bootstrap bias correction procedure does tend to decrease bias, as does Cressie’s predictor, particularly as the observation grids are filled in (moving from D₁ to D₂ and D₃ to D₄). Which of the two bias correction methods works better depends on the model.

Tables 1 and 2 show the results for φ(y)  y³, with both the isotropic (Table 1) and anisotropic (Table 2) covariograms. Looking at the ratios of expected values and biases, it is easily seen that all three predictors with a correction factor perform better than the naive predictor. Here, the Taylor’s expansion based predictor in (2.7) performs as well as the bootstrap predictors. This is due to the ‘nice’ behaviour of the function φ(y)  y³. The Taylor expansion yields

ϕ (Y (s_{0})) = \sum_{i = 0}^{\infty} \frac{ϕ^{(i)} (s_{0})}{i!} {(Y (s_{0}) - μ (s_{0}; β))}^{i}

= \sum_{i = 0}^{3} \frac{ϕ^{(i)} (s_{0})}{i!} {(Y (s_{0}) - μ (s_{0}; β))}^{i} .

Table 1

Results for φ(y) = y³, C(h) = exp(−θ||h||) withθ= 0.5

		D₁	D₂	D₃	D₄	D₅
Ratios	$\hat{Z}$	1.449	1.194	1.376	1.182	1.393
	${\hat{Z}}^{C}$	1.008	0.995	0.945	1.006	0.993
	${\hat{Z}}_{A}^{B C}$	1.007	0.993	0.944	1.008	0.996
	${\hat{Z}}_{M}^{B C}$	1.066	1.007	0.978	1.003	0.982
Biases	$\hat{Z}$	–1.217	–0.661	–0.976	–0.664	–1.136
	${\hat{Z}}^{C}$	–0.031	0.022	0.208	–0.027	0.029
	${\hat{Z}}_{A}^{B C}$	–0.029	0.027	0.210	–0.032	0.017
	${\hat{Z}}_{M}^{B C}$	–0.244	–0.027	0.079	–0.011	0.075
MSPEs	$\hat{Z}$	29.607	16.622	28.769	18.055	29.679
	${\hat{Z}}^{C}$	26.470	15.511	27.199	17.471	28.310
	${\hat{Z}}_{A}^{B C}$	26.575	15.518	27.285	17.548	28.422
	${\hat{Z}}_{M}^{B C}$	26.979	15.147	27.295	17.669	29.310

Source: Authors’ own

Table 2

Results for φ(y)  y³, C(h)  exp(–θ₁h₁–θ₂h₂) withθ (0.5, 1)′

		D₁	D₂	D₃	D₄	D₅
Ratios	$\hat{Z}$	1.351	1.133	1.331	1.114	1.386
	${\hat{Z}}^{C}$	0.995	1.006	0.962	1.006	0.986
	${\hat{Z}}_{A}^{B C}$	0.998	1.006	0.961	1.003	0.986
	${\hat{Z}}_{M}^{B C}$	1.027	1.002	0.960	1.009	1.005
Biases	$\hat{Z}$	–0.996	–0.464	–0.926	–0.405	–1.011
	${\hat{Z}}^{C}$	0.019	–0.025	0.145	–0.011	0.052
	${\hat{Z}}_{A}^{B C}$	0.008	–0.024	0.151	–0.015	0.050
	${\hat{Z}}_{M}^{B C}$	–0.099	–0.006	0.154	–0.037	–0.017
MSPEs	$\hat{Z}$	26.308	10.236	36.998	10.718	25.210
	${\hat{Z}}^{C}$	24.395	9.822	35.176	10.521	24.094
	${\hat{Z}}_{A}^{B C}$	24.578	9.855	35.303	10.561	24.156
	${\hat{Z}}_{M}^{B C}$	26.157	9.749	32.397	11.109	23.285

Source: Authors’ own

The same also holds for $ϕ (\hat{Y} (s_{0})) .$ When the expected value is taken for a bias calculation, the only terms left are when i  0 and i  2, providing the basis for Taylor’s expansion based predictor.

Finally, we note that the bias and MSPE tend to decrease as the observation grid is filled in, and that the MSPEs are comparable for all four predictors, which was seen in all simulations.

Figure 1 shows the biases for 1000 realizations of all four predictors for the isotropic covariogram. It is easy to see that the bias correction factors reduce the skewness in the distribution of the biases, bringing the outliers closer to the boxes. Recall that the most bias arises when Y(s₀) is far from its mean, and this effect increases as the derivative of φ(⋅) becomes large in absolute value. Figure 2 demonstrates the improvement in bias as a result of the corrected predictors, without a significant increase in the MSPE. With a transformation such as φ(y)  y³, the Taylor’s expansion based bias-corrected predictor is as effective as the bootstrap predictors, as we see in Tables 1 and 2. However, the predictor ${\hat{Z}}^{C}$ begins to fail in correcting the bias compared to the bootstrap predictors for exponential transformation functions, as implied by comparing the biases of the three predictors in Propositions 4.1 and 4.2.

Figure 1

Source: Authors’ own

Figure 2

Source: Authors’ own

This is further illustrated by Table 3 for the Matern covariogram. A similar performance is observed also for the anisotropic case (not reported here to save space; see Rister, 2010, for details). Note that the Taylor’s expansion based predictor is still an improvement over the naive predictor, but it is not comparable to the bootstrap based predictors.

Table 3

Results for φ(y)  e^y, C(h)  exp(–θ||h||) withθ 0.5

		D₁	D₂	D₃	D₄	D₅
Ratios	$\hat{Z}$	1.219	1.095	1.209	1.108	1.254
	${\hat{Z}}^{C}$	1.037	1.012	1.037	1.034	1.096
	${\hat{Z}}_{A}^{B C}$	0.971	0.971	0.972	0.999	1.030
	${\hat{Z}}_{M}^{B C}$	0.986	0.972	0.975	0.993	1.024
Biases	$\hat{Z}$	–0.780	–0.396	–0.730	–0.451	–0.933
	${\hat{Z}}^{C}$	–0.153	–0.052	–0.151	–0.152	–0.404
	${\hat{Z}}_{A}^{B C}$	0.130	0.136	0.122	0.007	–0.133
	${\hat{Z}}_{M}^{B C}$	0.060	0.131	0.107	0.030	–0.107
MSPEs	$\hat{Z}$	15.687	7.684	15.270	9.251	20.070
	${\hat{Z}}^{C}$	14.466	7.484	14.542	8.923	19.205
	${\hat{Z}}_{A}^{B C}$	14.539	7.564	14.694	8.942	19.027
	${\hat{Z}}_{M}^{B C}$	14.265	7.723	14.645	8.636	17.822

Source: Authors’ own

Table 4 shows the results for φ(y)  e^2y with the anisotropic covariogram. Here, we see that as the transformation function becomes more extreme, the predictors become less reliable. The naive predictor performs poorly, and Taylor’s expansion based predictor performs significantly worse than the bootstrap predictors. With this transformation function, we can no longer count on the later terms in Taylor expansion to be small, particularly when φ⁽ⁱ⁾ (y)  2ⁱ e^2y  i!. Compared to the earlier models, the accuracy of the bootstrap predictors has begun to deteriorate, although they both still outperform the existing predictors. Interestingly, the MSPE of the multiplicative predictor appears to be more volatile than the MSPE of the additive predictor. For the isotropic covariogram, similar results are seen (but not reported here). See Rister (2010) for the omitted simulation results.

Table 4

Results for φ(y)  e^2y, C(h)  exp(–θ₁h₁–θ₂h₂) withθ (0.5, 1)′

		D₁	D₂	D₃	D₄	D₅
Ratios	$\hat{Z}$	2.379	1.389	2.171	1.182	1.987
	${\hat{Z}}^{C}$	1.875	1.310	1.733	1.120	1.612
	${\hat{Z}}_{A}^{B C}$	1.010	1.020	0.974	0.896	0.884
	${\hat{Z}}_{M}^{B C}$	1.152	0.981	0.994	0.885	0.927
Biases	$\hat{Z}$	–32.542	–16.496	–27.957	–6.868	–22.382
	${\hat{Z}}^{C}$	–26.199	–13.929	–21.924	–4.767	–17.108
	${\hat{Z}}_{A}^{B C}$	–0.576	–1.141	1.390	5.164	5.940
	${\hat{Z}}_{M}^{B C}$	–7.428	1.166	0.300	5.792	3.526
MSPEs	$\hat{Z}$	174284	66089	34371	14490	29520
	${\hat{Z}}^{C}$	173045	65765	33763	14442	29284
	${\hat{Z}}_{A}^{B C}$	171144	64592	33067	14720	29599
	${\hat{Z}}_{M}^{B C}$	168512	55034	27128	24423	28719

Source: Authors’ own

Remark 5.1: It is evident from Tables 1–4 that the total contribution of the bias part to the MSPE is small (compared to the variance part). This is a typical phenomenon in prediction problems which must be kept in mind while comparing the numerical results here with mean squared errors in estimation problems. The main difference is that here the target is a random variable (cf. Z(s₀)) with a nonvanishing variance which often has the most contribution to the MSPE. However, this does not necessarily make the issue of bias correction unimportant. Indeed, there are many applications, particularly in governmental policy formulation and resource allocation, where a nontrivial prediction bias has serious practical implications. In such situations, the standard squared error loss function is not adequate and more weight must be associated with the bias component, making predictors with a smaller bias preferable. As shown in Tables 1–4, the proposed method provides an effective way of reducing the bias of the naive predictor, particularly in cases where the nonlinear transformation φ(⋅) has a faster growth rate (cf. Table 4). Further, even in terms of the MSPE, the bias-corrected predictors also have (slightly) smaller MSPEs in many (but not all) cases.

6 Oklahoma climatology data

For a real world data example, we look at a problem relating weather and energy consumption. Define Z(s) to be the ‘heating degree days’ (HDDs) at a location s, depending on the average daily outside air temperature T(s). HDD depends on a base temperature, often taken to be 65 Fahrenheit, such that

Z(s) = {\begin{array}{l} 65 ° F - T (s), & for T (s) < 65 ° F, \\ 0, & otherwise . \end{array}

(6.1)

An interpretation of HDD is as follows: if one location s₁ has an HDD value of 20 and another location s₂ has an HDD value of 40, a building at s₂ will require twice as much energy to heat as a similar building at s₁. As implied by the (6.1), HDDs are calculated on a daily basis. For more on HDDs, see Ristinen and Kraushaar (2006). Cooling degree days (CDDs) is a similar measurement regarding the amount of energy required to cool a building. Generally, HDDs are higher in the winter and CDDs in the summer.

Often, HDDs are accumulated over weeks, months or seasons, to study energy usage over an extended period of time, but we will focus on a single day of observations. Figure 3 shows the HDDs for 1 January 2009 at 116 observation sites in Oklahoma, provided by the OCS (2010), a joint venture between the University of Oklahoma and Oklahoma State University that tracks weather throughout the state. Suppose we wish to predict the daily HDD measurement at the points (–98,36), (–96,36.2), (–96,35.2) and (–97.4). These sites are marked by a ‘×’ symbol in Figure 3.

Figure 3

HDD observations in Oklahoma, 1 January 2009, with multiplication symbols (×) representing prediction points

There appears to be a definite relationship between the spatial location and the HDD measurement, with a higher number of HDDs corresponding to the sites in the northeastern part of the state. Further, Figure 4 shows that the distribution of HDDs is highly skewed. The next step is to identify a suitable transformation function φ(⋅).

Figure 4

Source: Authors’ own

We begin by using the transformation suggested by Box and Cox (1964), which has proven to be suitable for positive variables. If negative observations were present, the Box-Cox transformation would fail, and some modified methods, such as the one suggested by Yeo and Johnson (2000), should be tried. The Box-Cox transformation function is given by

g_{λ} (x_{i}) = {\begin{array}{l} \frac{x_{i}^{λ} - 1}{λ {(G M (X))}^{λ - 1}}, & if λ \neq 0, \\ G M (X) \log (x_{i}), & if λ = 0, \end{array}

where x  {x₁, …, x_n} and GM(x) is the geometric mean of x. The geometric mean is a constant that scales the data, and is often omitted. For the HDD data, the optimal value of λ was found to be –3.005236. The omission of the geometric mean for this data set leads to transformed observations with a small variance, causing computational problems later in the analysis. As a result, we included the geometric mean in g_λ (⋅).

Next define x ≡ Z  {Z(s₁), …, Z(s_n)}, Y(s)  g_λ Z(s)) ≡ φ^–1(Z(s)) and k ≡ (GM(Z))^λ–1,^{where λ  –3.005236. Then, it is easy to check that}

Z (s) = ϕ (Y (s)) = g_{λ}^{- 1} (Y (s)) = {(1 + Y (s) λ κ)}^{1 / λ} .

The distribution of transformed HDD observations is shown in Figure 5. It is clear that normality is plausible for Y(s).

Figure 5

Source: Authors’ own

Remark 6.1: As pointed out by a referee, Condition (C.3)(iv) does not hold for the transformation φ(⋅) above. However, by using a direct estimate over I ≡{y: 1 + λκy ≤ δ} (e.g., φ(y)≤ δ^{–1/λ}) and then Taylor’s expansion over I^c  {y: 1 + λκy  δ} as in the proof of Proposition 4.2, one can establish asymptotic unbiasedness of the bias-corrected predictors. Note that for y ∈ I^c, φ⁽¹⁾(y) is bounded by a constant multiple of δ^{–1–1/λ}, which satisfies (C.3)(iv). Now letting n → ∞ first and then, letting δ ↓ 0, one can conclude asymptotic negligibility of the second part. For the first part, one needs to use the key fact that the exponent –1/λ –1 which makes the Lebesgue integral over I to be O(δ^{–1–1/λ})  o(1) as δ ↓ 0.

Next, we calculate $\hat{ϵ} (s_{1}), ..., \hat{ϵ} (s_{n}),$ the estimated residuals, by estimating the mean parameter β. We assume a mean function depending on the latitude and longitude, such that Yµββ₀β₁X₁β₂X₂β₃X₁X₂

where X₁(s) is the longitude of s, X₂(s) is the latitude and µ(s; β)  β₀  β₁ X₁ (s) + β₂ X₂ (s) + β₃ X₁(s) X₂ (s). The estimated coefficients and p-values from a least squares regression model are given in Table 5, with all coefficients significantly different from 0.

Based on some exploratory work on the covariogram, we chose Chθ₁hθ₂

where ||h|| is the Euclidean distance between two locations. θ₁ represents the marginal variance (σ²) of  (s), whereas θ₂ represents the range of correlation. These parameters were estimated using weighted least squares. This covariogram is a member of the Matérn (1960) class, with fixed smoothness parameter equal to 0.5, which is believed to have considerable use in practice (cf. Schabenberger and Gotway, 2005; Stein, 1999).

Table 5

Linear model summary for HDD data

Effect	Coefficient	p-value
Intercept	429250.2	< 2 × 10^–16
Latitude	–11.125688	0.01626
Longitude	37.22961	0.00355
Interaction	0.3514275	0.00728

Source: Authors’ own

WLS procedures give $\hat{θ}$  (9.603985, 2.683276)′. We now have all information needed to calculate the predictors in (2.6), (2.7), (3.1) and (3.2). The predictors for each of the four prediction points are given in Table 6. Each predictor is considered to be an estimate of the amount of energy needed to heat a home at a given location. In this real data example, all four methods give qualitatively similar predicted values, with the bias-corrected prediction methods suggesting slightly lower energy needs. Specifically, ${\hat{Z}}^{C} (s_{0})$ estimates that between 0.14 and 0.25% more energy is required than ${\hat{Z}}^{C} (s_{0})$ does. On the other hand, the bootstrap predictors say that between 0.11 and 0.7% less energy is required than ${\hat{Z}}^{C} (s_{0})$ . While this does not translate into any significant difference in the energy expenditure, the example serves to illustrate the various steps and important aspects associated with implementation of the proposed methodology.

Table 6

Predictors for HDD data

	(–98,36)	(–96,36.2)	(–96,35.2)	(–97.4,34.5)
$\hat{Z} (s_{0})$	25.15465	28.23200	26.25444	25.41312
${\hat{Z}}^{C} (s_{0})$	25.21641	28.27352	26.30458	25.46030
${\hat{Z}}_{A}^{B C} (s_{0})$	25.02503	28.03837	26.19769	25.38501
${\hat{Z}}_{M}^{B C} (s_{0})$	25.02985	28.03462	26.19745	25.38438

Source: Authors’ own

7 Proofs

Here, we once again shall write $\hat{θ} = {\hat{θ}}_{n},$ $\hat{β} = {\hat{β}}_{n},$ $γ = γ_{n},$ $Σ = Σ_{n},$ $\tilde{Z} (s_{0}) = {\tilde{Z}}_{n} (s_{0}),$ $\hat{Z} {(s}_{0}) = {\hat{Z}}_{n} {(s}_{0}),$ etc. Denote the kth derivative of φ(⋅) by φ^(k)(⋅), k ≥ 1. We will also use φ′(⋅), φ″(⋅), φ″′(⋅) for φ⁽¹⁾(⋅), φ⁽²⁾(⋅) and φ⁽³⁾(⋅), respectively. Let 1(⋅) denote the indicator function, with 1(S)  0 or 1 accordingly as a statement S is false or true. Let K, K (⋅) denote generic constants in (0,∞) that may depend on some arguments (if any), but not on n.

7.1 Proof of Proposition 4.1

First, consider (4.4). Let N ∼ N (0, 1). Then, it is easy to check that for any k ≥ 1,

E | N |^{k} = \frac{{(\sqrt{2})}^{k}}{\sqrt{π}} Γ (\frac{k + 1}{2}),

where Γ(⋅) is the Gamma function. Note that by Jensen’s inequality, for any a, b ∈ (0, ∞),

\int_{0}^{\infty} x^{a} e^{- x} d x \leq {(\int_{0}^{\infty} x^{a + b} e^{- x} d x)}^{\frac{a}{a + b}} .

Thus,

Γ (a + 1) \leq {[Γ (a + b + 1)]}^{\frac{a}{a + b}},

for a ≥ 0 and b ≥ 0. Further, for any integer r ≥ 2, it is easy to check that

Γ (r - \frac{1}{2}) \geq \frac{Γ (r) Γ (3 / 2)}{r - 1} .

Hence, by (4.2), it follows that for x ∈ {σ_n, σ₀},

\sum_{k = 1}^{\infty} | d_{k} | 2^{k / 2} Γ (\frac{k + 1}{2}) x^{k}

= {[\sum_{j = 1}^{\infty} \sum_{k = 1}^{\infty} | d_{j} | | d_{k} | 2^{(k + j) / 2} Γ (\frac{k + 1}{2}) Γ (\frac{j + 1}{2}) x^{j + k}]}^{1 / 2}

\leq [\sum_{j = 1}^{\infty} \sum_{k = 1}^{\infty} | d_{j} | | d_{k} | 2^{(k + j) / 2} Γ {(\frac{k - 1}{2} + \frac{j - 1}{2} + 1)}^{\frac{k - 1}{k + j - 2}}

\times Γ {(\frac{k - 1}{2} + \frac{j - 1}{2} + 1)}^{\frac{j - 1}{k + j - 2}} x^{j + k}]^{1 / 2}

\leq {[\sum_{j = 1}^{\infty} \sum_{k = 1}^{\infty} | d_{j} | | d_{k} | 2^{(k + j) / 2} \frac{Γ (\frac{k + j - 1}{2}) (k + j) x^{k + j}}{Γ (\frac{3}{2})}]}^{1 / 2}

O(7.1)n

Note that $ϵ (s_{0}) \sim N (0, σ_{0}^{2})$ and ${\tilde{ϵ}}_{n} (s_{0}) \sim N (0, σ_{n}^{2}) .$ Also, by (7.1), we have

E \sum_{k = 1}^{\infty} | d_{k} | | w_{k} | < \infty,

for $w_{k} \in {ϵ {(s_{0})}^{k}, {\tilde{ϵ}}_{n} {(s_{0})}^{k}},$ k ≥ 1. Also, from (4.1), it follows that φ(⋅) is infinitely differentiable and that d_k  φ^(k) (µ(s₀; β))/k!, k ≥ 1. Hence using Fubini’s theorem, we get

Bias (\hat{Z} (s_{0}))

= E [\tilde{Z} (s_{0}) - Z (s_{0})]

= E [ϕ (μ (s_{0}; β) + {\tilde{ϵ}}_{n} (s_{0})) - ϕ (μ (s_{0}; β) + ϵ (s_{0}))]

= E (\sum_{k \geq 0} \frac{ϕ^{(k)} (μ (s_{0}; β))}{k!} {\tilde{ϵ}}_{n} {(s_{0})}^{k}) - E (\sum_{k \geq 0} \frac{ϕ^{(k)} (μ (s_{0}; β))}{k!} ϵ {(s_{0})}^{k})

= \sum_{k = 1}^{\infty} \frac{ϕ^{(2 k)} (μ (s_{0}; β))}{(2 k)!} {E ({\tilde{ϵ}}_{n} {(s_{0})}^{2 k}) - E (ϵ {(s_{0})}^{2 k})}

= \sum_{k = 1}^{\infty} \frac{ϕ^{(2 k)} (μ (s_{0}; β))}{(2 k)!} [\frac{(2 k)!}{k! 2^{k}} (σ_{n}^{2 k} - σ_{0}^{2 k})]

= \sum_{k = 1}^{\infty} \frac{ϕ^{(2 k)} (μ (s_{0}; β))}{k! 2^{k}} (σ_{n}^{2 k} - σ_{0}^{2 k}) .

(7.2)

Next, fix δ ∈ (0, 1/4). Let

A_{n} = {| μ (s_{0}; {\hat{β}}_{n}) - μ (s_{0}; β) | \leq δ} \cap^{​} {| {\tilde{ϵ}}_{n} (s_{0}) - {\tilde{ϵ}}_{n} (s_{0}) | \leq δ},

for n ≥ 1. Then,

| E [\hat{Z} (s_{0}) - \hat{Z} (s_{0})] 1 (A_{n}) |

= | E {ϕ (μ (s_{0}; {\hat{β}}_{n}) + {\hat{ϵ}}_{n} (s_{0})) - ϕ (μ (s_{0}; β) + {\hat{ϵ}}_{n} (s_{0}))} 1 (A_{n}) |

= | E {[μ (s_{0}; {\hat{β}}_{n}) - μ (s_{0}; β) + {\hat{ϵ}}_{n} (s_{0}) - {\hat{ϵ}}_{n} (s_{0})] \times

\int_{0}^{1} (1 - u) ϕ^{'} [u (μ (s_{0}; {\hat{β}}_{n}) + {\tilde{ϵ}}_{n} (s_{0})) + (1 - u) {μ (s_{0}; β) + {\tilde{ϵ}}_{n} (s_{0})}] d u} 1 (A_{n}) |

\leq \sqrt{2} {[E {(μ (s_{0}; {\hat{β}}_{n}) - μ (s_{0}, β))}^{2} 1 (A_{n}) + E {({\hat{ϵ}}_{n} (s_{0}) - {\hat{ϵ}}_{n} (s_{0}))}^{2} 1 (A_{n})]}^{1 / 2}

\times [\int_{0}^{1} E {ϕ^{'} (μ (s_{0}; β) + {\tilde{ϵ}}_{n} (s_{0}) + u [{μ (s_{0}; {\hat{β}}_{n}) - μ (s_{0}; β)}

+ {{\tilde{ϵ}}_{n} (s_{0}) - {\tilde{ϵ}}_{n} (s_{0})}]) 1 (A_{n})}^{2} d u]^{1 / 2} .

Note that by the Dominated Convergence Theorem,

E {[μ (s_{0}; {\hat{β}}_{n}) - μ (s_{0}; β)]}^{2} 1 (A_{n}) \to 0 a s n \to \infty .

Similarly,

E {({\hat{ϵ}}_{n} (s_{0}) - {\tilde{ϵ}}_{n} (s_{0}))}^{2} 1 (A_{n}) \to 0 a s n \to \infty .

Further, by (4.1), uniformly in u ∈ (0, 1),

E {ϕ^{'} (μ (s_{0}; β) + {\tilde{ϵ}}_{n} (s_{0}) + u [{μ (s_{0}; β) - μ (s_{0}; {\hat{β}}_{n})}

+ {{\hat{ϵ}}_{n} (s_{0}) - {\tilde{ϵ}}_{n} (s_{0})}]) 1 (A_{n})}^{2}

\leq E {[\sum_{k = 1}^{\infty} k | d_{k} | {| {\tilde{ϵ}}_{n} (s_{0}) | + 2 δ}^{k - 1}]}^{2}

\leq \sum_{k = 1}^{\infty} \sum_{j = 1}^{\infty} k j | d_{k} d_{j} | E {(| {\tilde{ϵ}}_{n} (s_{0}) | + 2 δ)}^{k + j - 2}

\leq \sum_{k = 1}^{\infty} \sum_{j = 1}^{\infty} k j | d_{k} d_{j} | 2^{k + j - 2} {E | {\tilde{ϵ}}_{n} (s_{0}) |^{k + j - 2} + {(2 δ)}^{k + j - 2}}

\leq \sum_{k = 1}^{\infty} \sum_{j = 1}^{\infty} k j | d_{k} d_{j} | 2^{3 (k + j - 2) / 2} Γ (\frac{k + j - 1}{2}) σ_{n}^{j + k - 2} + K (δ)

Hence, it follows that

| E [{\hat{Z}}_{n} (s_{0}) - {\tilde{Z}}_{n} (s_{0})] 1 (A_{n}) | = o (1) .

(7.3)

Also, by the Cauchy-Schwarz inequality and Condition 4.3,

| E [{\hat{Z}}_{n} (s_{0}) - {\tilde{Z}}_{n} (s_{0})] 1 (A_{n}^{c}) |

\leq {[E {({\hat{Z}}_{n} (s_{0}))}^{2}]}^{1 / 2} {[P (A_{n}^{c})]}^{1 / 2} + {[E {({\tilde{Z}}_{n} (s_{0}))}^{2}]}^{1 / 2} {[P (A_{n}^{c})]}^{1 / 2}

o(7.4)

From (7.2) to (7.4), (4.4) follows. The proof of (4.5) is similar and hence is omitted.

7.2 Proof of Proposition 4.2

First, consider the multiplicative bias-corrected predictor ${\hat{Z}}_{M, n}^{B C} (s_{0})$ and its ideal version ${\tilde{Z}}_{M, n}^{B C} (s_{0})$ . Note that, with $c_{n} \equiv E {\tilde{Z}}_{n} (s_{0}) / E Z (s_{0})$ ,

E ({\tilde{Z}}_{M, n}^{B C} (s_{0}) - Z (s_{0})) = c E {\tilde{Z}}_{n} (s_{0}) - E Z (s_{0}) = 0,

so that the ideal version ${\tilde{Z}}_{M, n}^{B C} (s_{0})$ is unbiased. It is easy to see that

Bias ({\hat{Z}}_{M, n}^{B C} (s_{0}))

= E ({\hat{Z}}_{M, n}^{B C} (s_{0}) - Z (s_{0}))

= E {\hat{c}}_{n} {\tilde{Z}}_{n} (s_{0}) - E Z (s_{0})

= E {(E * Z * (s_{0}) - E Z (s_{0})) {\hat{Z}}_{n} (s_{0}) / E * {\hat{Z}}_{n}^{*} (s_{0})}

- \frac{E Z (s_{0})}{E {\hat{Z}}_{n} (s_{0})} E [{E * {\hat{Z}}_{n}^{*} (s_{0}) - E^n (s_{0})}^n (s_{0}) / E * {\hat{Z}}_{n}^{*} (s_{0})]

I_1nI_2n

Note that $\hat{ϵ} (s_{0}) = \sum_{i = 1}^{n} λ_{i n} ({\hat{θ}}_{n}) ϵ (s_{i}),$ where λ_n(θ)  (λ_1n(θ), …, λ_nn (θ))′_{is as given by Condition 4.2 Hence}

E * {\hat{Z}}_{n} * (s_{0})

= E * ϕ (μ (s_{0}; {\hat{β}}_{n} *) + \sum_{i = 1}^{n} λ_{i n} ({\hat{θ}}_{n} *) [ϵ * (s_{i}) - {μ (s_{i}; {\hat{β}}_{n} *) - μ (s_{i}; {\hat{β}}_{n})}])

= E * ϕ (μ (s_{0}; {\hat{β}}_{n}) + \sum_{i = 1}^{n} λ_{i n} ({\hat{θ}}_{n}) ϵ * (s_{i}) + \sum_{i = 1}^{n} (λ_{i n} ({\hat{θ}}_{n} *) - λ_{i n} ({\hat{θ}}_{n})) ϵ * (s_{i})

- \sum_{i = 1}^{n} λ_{i n} ({\hat{θ}}_{n} *) {μ (s_{i}; {\hat{β}}_{n}^{} *) - μ (s_{i}; {\hat{β}}_{n})} + [μ (s_{0}; {\hat{β}}_{n} *) - μ (s_{0}; {\hat{β}}_{n})])

= g ({\hat{β}}_{n}, τ_{n} ({\hat{θ}}_{n})) + R_{1 n},

where $τ_{n}^{2} (θ) \equiv E_{θ} {(\sum_{i = 1}^{n} λ_{i n} (θ) ϵ (s_{i}))}^{2} = E_{θ} \tilde{Z} {(s_{0})}^{2},$ and where, on the set ${‖ {\hat{θ}}_{n} - θ ‖ < δ_{0} / 2},$ for any δ ∈ (0, δ₀/2), the remainder term R_1n admits the bound:

| R_{1 n} | \leq E * [(\sum_{i = 1}^{n} α_{i} {K_{0} ‖ θ_{n} * - {\hat{θ}}_{n} ‖ | ϵ * (s_{i}) | + K_{1} ‖ {\hat{β}}_{n} * - {\hat{β}}_{n} ‖} + K_{1} ‖ {\hat{β}}_{n} * - {\hat{β}}_{n} ‖)]

\times Φ^{(1)} (| μ (s_{0}; {\hat{β}}_{n}) | + | \sum_{i = 1}^{n} λ_{i n} ({\hat{θ}}_{n}) ϵ * (s_{i}) | + δ K_{0} \sum_{i = 1}^{n} α_{i} | ϵ * (s_{i}) | + K_{4} ‖ {\hat{β}}_{n} * - {\hat{β}}_{n} ‖)

\times 1 [(‖ θ_{n} * - {\hat{θ}}_{n} ‖ \leq δ)]

+ E * [(2 K_{0} \sum_{i = 1}^{n} α_{i} | ϵ * (s_{i}) | + K_{4} ‖ {\hat{β}}_{n} * - {\hat{β}}_{n} ‖)

\times Φ^{(1)} (| μ (s_{0}; {\hat{β}}_{n}) | + 3 K_{0} \sum_{i = 1}^{n} α_{i} | ϵ * (s_{i}) | + K_{4} ‖ \hat{β} *_{n} - {\hat{β}}_{n} ‖) (‖ θ_{n} * - {\hat{θ}}_{n} ‖ > δ)]

\leq K [δ {E * {(\sum_{i = 1}^{n} α_{i} | ϵ * (s_{i}) |)}^{2}}^{1 / 2} + {(E * ‖ {\hat{β}}_{n} * - {\hat{β}}_{n} ‖^{2})}^{1 / 2}]

\times [E * {Φ^{(1)} (| μ (s_{0}; {\hat{β}}_{n}) | + | \sum_{i = 1}^{n} λ_{i n} ({\hat{θ}}_{n}) ϵ * (s_{i}) |

+ δ K_{0} \sum_{i = 1}^{n} α_{i} | ϵ * (s_{i}) | + K_{4} ‖ {\hat{β}}_{n} * - {\hat{β}}_{n} ‖)}^{2}]^{1 / 2}

+ K {{[E * {(\sum_{i = 1}^{n} α_{i} | ϵ * (s_{i}) |)}^{4}]}^{1 / 4} + {(E ‖ {\hat{β}}_{n} * - {\hat{β}}_{n} ‖^{4})}^{1 / 4}} {[P * (‖ θ_{n} * - {\hat{θ}}_{n} ‖)]}^{1 / 4}

\leq K [{(E * {Φ^{(1)} (| μ (s_{0}; {\hat{β}}_{n}) | + 3 K_{0} \sum_{i = 1}^{n} α_{i} | ϵ * (s_{i}) | + K_{4} ‖ {\hat{β}}_{n}^{*} - {\hat{β}}_{n} ‖)}^{2})}^{1 / 2}]

\times [{(\sum_{i = 1}^{n} α_{i} E * ϵ * {(s_{i})}^{4})}^{1 / 4} {(\sum_{i = 1}^{n} α_{i})}^{3 / 4} {δ + {(P * (‖ θ_{n} * - {\hat{θ}}_{n} ‖ > δ))}^{1 / 4}}

+ {(E * ‖ {\hat{β}}_{n} * - {\hat{β}}_{n} ‖^{4})}^{1 / 4}]

Kδ(7.5)

for all n ≥ n₀ for some n₀  n₀ (K₀, …, K₃, α_∞, δ) and K  K (K₀, …, K₄). By similar arguments,

E {\hat{Z}}_{n} (s_{0}) = g (β, τ_{n} (θ)) + r_{1 n},

(7.6)

where r_1n ≤ Kδ for all n ≥ n₁ for some n₁  n₁ (K₀, …, K₃, α_∞, δ) ≥ 1.

Now, consider I_1n. Note that by (4.2),

| I_{1 n} | = | E [{g ({\hat{β}}_{n}, C (0, {\hat{θ}}_{n})) - g (β, C (0, θ))} {\hat{Z}}_{n} (s_{0}) / E * \hat{Z} * (s_{0})] |

\leq {{(E | g ({\hat{β}}_{n}, C (0, {\hat{θ}}_{n})) - g (β, C (0, θ)) |^{2})}^{1 / 2} {(E | {\hat{Z}}_{n} (s_{0}) / E * \hat{Z} * (s_{0}) |^{2})}^{1 / 2}}

Next, consider I_2n. By (7.5), (7.6) and (4.2), we have

| \frac{E {\hat{Z}}_{n} (s_{0})}{E Z (s_{0})} I_{2 n} | = | E [{E * \hat{Z} *_{n} (s_{0}) - E {\hat{Z}}_{n} (s_{0})} {\hat{Z}}_{n} (s_{0}) / E * {\hat{Z}}_{n}^{*} (s_{0})] |

\leq {{(E {[g ({\hat{β}}_{n}, τ_{n} ({\hat{θ}}_{n})) - g (β, τ_{n} (θ))]}^{2})}^{1 / 2} + [12 (1 + K_{2}^{1 / 2}) K_{3}^{1 / 3} α_{\infty} K_{0}] δ}

\times {[E {({\hat{Z}}_{n} (s_{0}) / E * Z *_{n} (s_{0}))}^{2}]}^{1 / 2}

by letting n → ∞ first and then δ ↓ 0. Hence, the asymptotic unbiasedness of ${\hat{Z}}_{M, n}^{B C}$ now follows.

The proof for ${\hat{Z}}_{A, n}^{B C}$ is similar (and in fact, somewhat simpler) and is thus omitted. This completes the proof of Proposition 4.2.

Footnotes

Acknowledgements

Research partially supported by NSF grant no. DMS–1007703 and NSA grant no. H98230–11–1–0130.

References

Bandyopadhyay

Lahiri

(2010) Resampling-based bias-corrected time series prediction. Journal of Statistical Planning and Inference, 140, 3775–88.

Box

GEP

Cox

(1964) An analysis of transformations. Journal of the Royal Statistical Society B, 26, 211–43.

Cressie

NAC

(1991) Statistics for spatial data, 2nd edn. Wiley: New York.

Hall

Maiti

(2006) On parametric bootstrap methods for small area prediction. Journal of the Royal Statistical Society: Series B, 68, 221–38.

Lahiri

Maiti

Katzoff

Parsons

(2007) Resampling-based empirical prediction: an application to small area estimation. Biometrika, 94, 469–85.

Lahiri

1st edn. (2003) Resampling methods for dependent data. New York: Springer.

Lehmann

Casella

(1998) Theory of point estimation. New York: Springer.

Matérn

(1960) Spatial variation. Meddelanden fran Statens Skogsforskningsinstitut, 49.

OCS (2010) Mesonet. Oklahoma Climatological Survey.

10.

Pfanzagl

(2000) On local uniformity for estimators and confidence limits. Journal of Statistical Planning and Inference, 84, 27–63.

11.

Rister

(2010) Resampling methodology in spatial prediction and repeated measures time series. PhD thesis, Texas A&M University.

12.

Ristinen

Kraushaar

Second edition. (2006) Energy and the environment, Hoboken, NJ: Wiley.

13.

Schabenberger

Gotway

(2005). Statistical methods for spatial data analysis. Boca Raton: Chapman & Hall/CRC.

14.

Stein

(1999). Interpolation of spatial data: some theory for kriging. New York: Springer.

15.

Yeo

I-K

Johnson

(2000) A new family of power transformations to improve normality or symmetry. Biometrika, 87, 954–59.