lgnpp: Unbiased Back-Transformations for Lognormal Kriging

In georob: Robust Geostatistical Analysis of Spatial Data

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 lognormal 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 (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	logical 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 lognormal point kriging is used:

hatZ(s)=exp(hatS(s) + 1/2(hatσ^2_n + hatσ^2 - Var[hatS(s)])),

Cov[ Z(s_i)-hatZ(s_i), Z(s_j)-hatZ(s_j)] = μ(s_i) μ(s_j)

*{ exp(Cov[Y(s_i), Y(s_j)]) - 2 exp(Cov[hatY(s_i), Y(s_j)]) + exp(Cov[hatY(s_i), hatY(s_j)]) },

where hatY and hatZ denote the log- and back-transformed predictions of the signal, and

μ(s) = exp( x(s)^T hatβ + 1/2 Var[S(s)]).

The expressions for the required covariance terms can be found in the Appendices of Nussbaum et al. (2012). Instead of the signal S(s), predictions of the log-transformed response S(s) or the estimated trend x(s)^T hatβ of the log-transformed data can be back-transformed (see georobIntro). 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)

hatZ(B) = exp( hatS(B) + 1/2 {Var[S(s)] + hatβ^T M(B) hatβ - Var[hatS(B)]}),

E[ { Z(B)- hatZ(B) }^2 ] = μ(B)^2 { \exp(Var[S(B)]) - 2 \exp(Cov[S(B),hatS(B)]) + \exp(Var[hatS(B)]) }

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

M(B) = 1/|B| int_B (x(s)-x(B) (x(s)-x(B))^T) ds

with

x(B) = 1/|B| int_B x(s) ds

and

μ(B) = exp( x(B)^T hatβ + 1/2 Var[S(B)]).

These equations are based on the assumption that both the point data Z(s) and the block means Z(B) follow lognormal 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 M(B), one needs the covariates on a fine grid within the block B. The covariates are passed to lgnpp as argument newdata, where newdata can be any spatial data frame accepted by predict.georob. For evaluating M(B) the geometry of the blocks is taken from the polygons slot of the SpatialPolygonsDataFrame passed as object to lgnpp.

3. Backtransformation 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 refer to points that lie in a single block. Hence, one uses the approximation

hatZ(B) = 1/K sum_{s_i in B} hatZ(s_i)

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

E[{Z(B) - hatZ(B)}^2] = 1/K^2 sum_{s_i in B} sum_{s_j in B} Cov[ Z(s_i)-hatZ(s_i), Z(s_j)-hatZ(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 B. lgnpp then uses the following approximation for the mean squared error (see also Appendix E of Nussbaum et al., 2012):

• Prediction results are passed as object to lgnpp only for a random sample of points in B (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 Z(B).

• The mean squared error is then approximated by

  E[{Z(B) - hatZ(B)}^2] = 1/K^2 sum_{s_i in B} E[ { Z(s_i) - hatZ(s_i)}^2]

  + (K-1)/(K k (k-1)) sum_{s_i in sample} sum_{s_j in sample, s_j != s_i} Cov[ Z(s_i)-hatZ(s_i), Z(s_j)-hatZ(s_j)]

  The first term of the RHS 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.

• If the prediction results are not available for the complete set of points in B then all.pred may be equal to K. The block mean is then approximated by

  hatZ(B) = sum_{s_i in B} hatZ(s_i)

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

  1(k K) sum_{s_i in sample} E[ { Z(s_i) - hatZ(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 isFALSE and all.pred is equal to NULL 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 a named numeric vector with two elements:

• mean: the back-transformed block kriging estimate, see Details.

• mse: the (approximated) block kriging variance, 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 [e-mail].

References

Cressie, N. (2006) Block Kriging for Lognormal Spatial Processes. Mathematical Geology, 38, 413–443.

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.

Nussbaum, M., Papritz, A., Baltensweiler, A. and Walthert, L. (2012) Organic Carbon Stocks of Swiss Forest Soils, Institute of Terrestrial Ecosystems, ETH Zurich and Swiss Federal Institute for Forest, Snow and Landscape Research (WSL), pp. 51. http://e-collection.library.ethz.ch/eserv/eth:6027/eth-6027-01.pdf

See Also

georobIntro for a description of the model and a brief summary of the algorithms; georob for (robust) fitting of spatial linear models; predict.georob for computing robust kriging predictions.

Examples

## Not run: 
data(meuse )

data(meuse.grid)
meuse.grid.pixdf <- SpatialPixelsDataFrame(points = meuse.grid[, 1:2],
    data = meuse.grid[, -(1:2)])

library(constrainedKriging)
data(meuse.blocks)

r.logzn.rob <- georob(log(zinc) ~ sqrt(dist), data = meuse, locations = ~ x + y,
    variogram.model = "exponential", param = c( variance = 0.15, nugget = 0.05, scale = 200 ),
    tuning.psi = 1., control = georob.control(cov.bhat = TRUE, full.cov.bhat = TRUE,
        cov.bhat.betahat = TRUE, aux.cov.pred.target = TRUE))
        
## point predictions of log(Zn)
r.pred.points <- predict(r.logzn.rob, newdata = meuse.grid.pixdf, extended.output = TRUE,
    full.covmat = TRUE)
str(r.pred.points$pred@data)

## back-transformation of point predictions
r.backtf.pred.points <- lgnpp(r.pred.points)
str(r.backtf.pred.points$pred@data)

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

## predicting mean Zn content for whole area
r.block <- lgnpp(r.pred.points, is.block = TRUE, all.pred = r.backtf.pred.points[["pred"]])
r.block

## block predictions of log(Zn)
r.pred.block <- predict(r.logzn.rob, newdata = meuse.blocks, extended.output = TRUE,
    pwidth = 75, pheight = 75)
r.backtf.pred.block <- lgnpp(r.pred.block, newdata = meuse.grid)

spplot(r.backtf.pred.block, zcol = "lgn.pred", main = "block means Zn content")
## End(Not run)

georob documentation built on May 2, 2019, 6:53 p.m.