lgnpp: Unbiased Back-Transformations for Log-normal Kriging

Description Usage Arguments Details Value Author(s) References See Also Examples

View source: R/georob.lgnpp.R

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

1
2
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 log-normal point Kriging is used:

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

Cov[ U(s_i)-hatU(s_i), U(s_j)-hatU(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 hatU denote the log- and back-transformed predictions of the signal, and

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

The expressions for the required covariance terms can be found in the Appendices of Nussbaum et al. (2012). Instead of the signal Z(s), predictions of the log-transformed response Z(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, Appendix C)

hatU(A) = exp( hatS(A) + 1/2 {Var[Z(s)] + hatβ^T M(A) hatβ - Var[hatS(A)]}),

E[ { U(A)- hatU(A) }^2 ] = μ(A)^2 { \exp(Var[Z(A)]) - 2 \exp(Cov[Z(A),hatS(A)]) + \exp(Var[hatS(A)]) }

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

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

within the block A where

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

and

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

This back-transformation is based on the assumption that both the point data U(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 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 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

hatU(A) = 1/K sum_{s_i in A} hatU(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

E[{U(A) - hatU(A)}^2] = 1/K^2 sum_{s_i in A} sum_{s_j in A} Cov[ U(s_i)-hatU(s_i), U(s_j)-hatU(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., 2012):

Value

If is.block is FALSE 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:

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

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 [email protected].

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
## Not run: 
data(meuse)

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

library(constrainedKriging)
data(meuse.blocks)

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 <- predict(r.logzn.rob, newdata = meuse.grid.pixdf, 
    control = control.predict.georob(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, 
    control = control.predict.georob(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 29, 2017, 10:51 a.m.