R: Unbiased Back-Transformations for Log-normal Kriging

lgnpp {georob}

R Documentation

Unbiased Back-Transformations for Log-normal Kriging

Description

The function lgnpp back-transforms point or block Kriging predictions of a log-transformed response variable computed by predict.georob. Alternatively, the function averages log-normal point Kriging predictions for a block and approximates the mean squared prediction error of the block mean.

Usage

lgnpp(object, newdata, locations, is.block = FALSE, all.pred = NULL,
    extended.output = FALSE)

Arguments

`object`	an object with Kriging predictions of a log-transformed response variable as obtained by `predict(georob-object, ...)`.
`newdata`	an optional object as passed as argument `newdata` to `predict.georob`, see Details.
`locations`	an optional one-sided formula specifying what variables of `newdata` are the coordinates of the prediction points, see `predict.georob`.
`is.block`	an optional logical scalar (default `FALSE`) specifying whether point predictions contained in `object` are considered to belong to a single block and should be averaged after back-transformation. Ignored if `object` contains block Kriging predictions, see Details.
`all.pred`	an optional positive integer or an object as obtained by `lgnpp(predict(georob-object, ...))`, see Details.
`extended.output`	a logical scalar controlling whether the covariance matrix of the errors of the back-transformed point predictions is added as an attribute to the result, see Details.

Details

The function lgnpp performs three tasks:

1. Back-transformation of point Kriging predictions of a log-transformed response

The usual, marginally unbiased back-transformation for log-normal point Kriging is used:

\widehat{U}(\boldsymbol{s}) = \exp( \widehat{Z}(\boldsymbol{s}) + 1/2 ( \mathrm{Var}_{\hat{\theta}}[ Z(\boldsymbol{s})] - \mathrm{Var}_{\hat{\theta}}[\widehat{Z}(\boldsymbol{s})])),

\mathrm{Cov}_{\hat{\theta}}[ U(\boldsymbol{s}_i) - \widehat{U}(\boldsymbol{s}_i), U(\boldsymbol{s}_j) - \widehat{U}(\boldsymbol{s}_j) ] = \mu_{\hat{\theta}}(\boldsymbol{s}_i) \mu_{\hat{\theta}}(\boldsymbol{s}_j)

\times \{ \exp(\mathrm{Cov}_{\hat{\theta}}[Z(\boldsymbol{s}_i),Z(\boldsymbol{s}_j)]) -2\exp(\mathrm{Cov}_{\hat{\theta}}[\widehat{Z}(\boldsymbol{s}_i),Z(\boldsymbol{s}_j)]) +\exp(\mathrm{Cov}_{\hat{\theta}}[\widehat{Z}(\boldsymbol{s}_i),\widehat{Z}(\boldsymbol{s}_j)]) \},

where \widehat{Z} and \widehat{U} denote the log- and back-transformed predictions of the signal, and

\mu_{\hat{\theta}}(\boldsymbol{s}) \approx \exp(\boldsymbol{x}(\boldsymbol{s})\mathrm{^T}\widehat{\boldsymbol{\beta}} + 1/2 \mathrm{Var}_{\hat{\theta}}[Z(\boldsymbol{s})]).

The expressions for the required covariance terms can be found in the Appendices of Nussbaum et al. (2014). Instead of the signal Z(\boldsymbol{s}), predictions of the log-transformed response Y(\boldsymbol{s}) or the estimated trend \boldsymbol{x}(\boldsymbol{s})^\mathrm{T}\widehat{\boldsymbol{\beta}} of the log-transformed data can be back-transformed (see georobPackage). The above transformations are used if object contains point Kriging predictions (see predict.georob, Value) and if is.block = FALSE and all.pred is missing.

2. Back-transformation of block Kriging predictions of a log-transformed response

Block Kriging predictions of a log-transformed response variable are back-transformed by the approximately unbiased transformation proposed by Cressie (2006, Appendix C)

\widehat{U}(A) = \exp( \widehat{Z}(A) + 1/2 \{ \mathrm{Var}_{\hat{\theta}}[Z(\boldsymbol{s})] + \widehat{\boldsymbol{\beta}}\mathrm{^T} \boldsymbol{M}(A) \widehat{\boldsymbol{\beta}} - \mathrm{Var}_{\hat{\theta}}[\widehat{Z}(A)] \}),

\mathrm{E}_{\hat{\theta}}[\{U(A) - \widehat{U}(A))^2] = \mu_{\hat{\theta}}(A)^2 \{ \exp(\mathrm{Var}_{\hat{\theta}}[Z(A)]) - 2 \exp(\mathrm{Cov}_{\hat{\theta}}[\widehat{Z}(A),Z(A)]) + \exp(\mathrm{Var}_{\hat{\theta}}[\widehat{Z}(A)]) \}

where \widehat{Z}(A) and \widehat{U}(A) are the log- and back-transformed predictions of the block mean U(A), respectively, \boldsymbol{M}(A) is the spatial covariance matrix of the covariates

\boldsymbol{M}(A) = 1/|A| \int_A ( \boldsymbol{x}(\boldsymbol{s}) - \boldsymbol{x}(A) ) ( \boldsymbol{x}(\boldsymbol{s}) - \boldsymbol{x}(A) )\mathrm{^T} \,d\boldsymbol{s}

within the block A where

\boldsymbol{x}(A) = 1/|A| \int_A \boldsymbol{x}(\boldsymbol{s}) \,d\boldsymbol{s}

and

\mu_{\hat{\theta}}(A) \approx \exp(\boldsymbol{x}(A)\mathrm{^T} \widehat{\boldsymbol{\beta}} + 1/2 \mathrm{Var}_{\hat{\theta}}[Z(A)]).

This back-transformation is based on the assumption that both the point data U(\boldsymbol{s}) and the block means U(A) follow log-normal laws, which strictly cannot hold. But for small blocks the assumption works well as the bias and the loss of efficiency caused by this assumption are small (Cressie, 2006; Hofer et al., 2013).

The above formulae are used by lgnpp if object contains block Kriging predictions in the form of a SpatialPolygonsDataFrame. To approximate \boldsymbol{M}(A), one needs the covariates on a fine grid for the whole study domain in which the blocks lie. The covariates are passed lgnpp as argument newdata, where newdata can be any spatial data frame accepted by predict.georob. For evaluating \boldsymbol{M}(A) the geometry of the blocks is taken from the polygons slot of the
SpatialPolygonsDataFrame passed as object to lgnpp.

3. Back-transformation and averaging of point Kriging predictions of a log-transformed response

lgnpp allows as a further option to back-transform and average point Kriging predictions passed as object to the function. One then assumes that the predictions in object refer to points that lie in a single block. Hence, one uses the approximation

\widehat{U}(A) \approx \frac{1}{K} \sum_{s_i \in A} \widehat{U}(\boldsymbol{s}_i)

to predict the block mean U(A), where K is the number of points in A. The mean squared prediction error can be approximated by

\mathrm{E}_{\hat{\theta}}[\{U(A) - \widehat{U}(A)\}^2] \approx \frac{1}{K^2} \sum_{s_i \in A} \sum_{s_j \in A} \mathrm{Cov}_{\hat{\theta}}[ U(\boldsymbol{s}_i) - \widehat{U}(\boldsymbol{s}_i), U(\boldsymbol{s}_j) - \widehat{U}(\boldsymbol{s}_j) ].

In most instances, the evaluation of the above double sum is not feasible because a large number of points is used to discretize the block A. lgnpp then uses the following approximations to compute the mean squared error (see also Appendix E of Nussbaum et al., 2014):

Point prediction results are passed as object to lgnpp only for a random sample of points in A (of size k), for which the evaluation of the above double sum is feasible.
The prediction results for the complete set of points within the block are passed as argument all.pred to lgnpp. These results are used to compute \widehat{U}(A).
The mean squared error is then approximated by

\mathrm{E}_{\hat{\theta}}[\{U(A) - \widehat{U}(A)\}^2] \approx \frac{1}{K^2} \sum_{s_i \in A} \mathrm{E}_{\hat{\theta}}[ \{ U(\boldsymbol{s}_i) - \widehat{U}(\boldsymbol{s}_i)\}^2]

+ \frac{K-1}{K k (k-1)} \sum_{s_i \in \mathrm{sample}}\sum_{s_j \in \mathrm{sample}, s_j \neq s_i} \mathrm{Cov}_{\hat{\theta}}[ U(\boldsymbol{s}_i) - \widehat{U}(\boldsymbol{s}_i), U(\boldsymbol{s}_j) - \widehat{U}(\boldsymbol{s}_j) ].

The first term of the RHS (and \widehat{U}(A)) can be computed from the point Kriging results contained in all.pred, and the double sum is evaluated from the full covariance matrices of the predictions and the respective targets, passed to lgnpp as object (one has to use the arguments control=control.predict.georob(full.covmat=TRUE) for predict.georob when computing the point Kriging predictions stored in object).
If the prediction results are not available for the complete set of points in A then all.pred may be equal to K. The block mean is then approximated by

\widehat{U}(A) \approx \frac{1}{k} \sum_{s_i \in \mathrm{sample}} \widehat{U}(\boldsymbol{s}_i)

and the first term of the RHS of the expression for the mean squared error by

\frac{1}{kK} \sum_{s_i \in \mathrm{sample}} \mathrm{E}_{\hat{\theta}}[ \{ U(\boldsymbol{s}_i) - \widehat{U}(\boldsymbol{s}_i)\}^2].
By drawing samples repeatedly and passing the related Kriging results as object to lgnpp, one can reduce the error of the approximation of the mean squared error.

Value

If is.block is FALSE and all.pred is equal to NULL lgnpp returns an updated object of the same class as object (see section Value of predict.georob). The data frame with the point or block Kriging predictions is complemented by lgnpp with the following new components:

lgn.pred: the back-transformed Kriging predictions of a log-transformed response.
lgn.se: the standard errors of the back-transformed predictions.
lgn.lower, lgn.upper: the bounds of the back-transformed prediction intervals.

If is.block is TRUE or all.pred not equal to NULL lgnpp returns a named numeric vector with two elements:

mean: the back-transformed block Kriging estimate, see Details.
se: the (approximated) block Kriging standard error, see Details.

If extended.output is TRUE then the vector is supplemented with the attribute mse.lgn.pred that contains the full covariance matrix of the back-transformed point prediction errors.

Author(s)

Andreas Papritz papritz@retired.ethz.ch.

References

Cressie, N. (2006) Block Kriging for Lognormal Spatial Processes. Mathematical Geology, 38, 413–443, doi:10.1007/s11004-005-9022-8.

Hofer, C., Borer, F., Bono, R., Kayser, A. and Papritz, A. 2013. Predicting topsoil heavy metal content of parcels of land: An empirical validation of customary and constrained lognormal block Kriging and conditional simulations. Geoderma, 193–194, 200–212, doi:10.1016/j.geoderma.2012.08.034.

Nussbaum, M., Papritz, A., Baltensweiler, A. and Walthert, L. (2014) Estimating soil organic carbon stocks of Swiss forest soils by robust external-drift kriging. Geoscientific Model Development, 7, 1197–1210. doi:10.5194/gmd-7-1197-2014.

Examples

data(meuse)

data(meuse.grid)
coordinates(meuse.grid) <- ~x+y
meuse.grid.pixdf <- meuse.grid
gridded(meuse.grid.pixdf) <- TRUE

data(meuse.blocks, package = "constrainedKriging")

r.logzn.rob <- georob(log(zinc) ~ sqrt(dist), data = meuse, locations = ~ x + y,
    variogram.model = "RMexp", param = c(variance = 0.15, nugget = 0.05, scale = 200),
    tuning.psi = 1., control = control.georob(cov.bhat = TRUE, full.cov.bhat = TRUE))

## point predictions of log(Zn)
r.pred.points.1 <- predict(r.logzn.rob, newdata = meuse.grid.pixdf,
    control = control.predict.georob(extended.output = TRUE))
str(r.pred.points.1, max = 3)

## back-transformation of point predictions
r.backtf.pred.points <- lgnpp(r.pred.points.1)
str(r.backtf.pred.points, max = 3)

spplot(r.backtf.pred.points, zcol = "lgn.pred", main = "Zn content")

## predicting mean Zn content for whole area
if(interactive()){
  ## example is run only in interactive session because cpu times exceeds 5 s
  ## recompute point predictions with argument full.covmat = TRUE
  r.pred.points.2 <- predict(r.logzn.rob, newdata = meuse.grid.pixdf,
      control = control.predict.georob(extended.output = TRUE, full.covmat = TRUE))
  str(r.pred.points.2, max = 3)
  r.block <- lgnpp(r.pred.points.2, is.block = TRUE, all.pred = r.backtf.pred.points@data)
  r.block
}

## block predictions of log(Zn)
if(interactive()){
  ## example is run only in interactive session because cpu times exceeds 5 s
  r.pred.block <- predict(r.logzn.rob, newdata = meuse.blocks,
      control = control.predict.georob(extended.output = TRUE,
          pwidth = 75, pheight = 75, mmax = 50))
  r.backtf.pred.block <- lgnpp(r.pred.block, newdata = meuse.grid)

  spplot(r.backtf.pred.block, zcol = "lgn.pred", main = "block means Zn content")
}

[Package georob version 0.3-19 Index]