Description Usage Arguments Details Value Author(s) References See Also Examples
The function lgnpp
backtransforms point or block Kriging
predictions of a logtransformed 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.
1 2 
object 
an object with Kriging predictions of a logtransformed
response variable as obtained by

newdata 
an optional object as passed as argument 
locations 
an optional onesided formula specifying what variables
of 
is.block 
an optional logical (default 
all.pred 
an optional positive integer or an object as obtained by

extended.output 
logical controlling whether the covariance matrix of the errors of the backtransformed point predictions is added as an attribute to the result, see Details. 
The function lgnpp
performs three tasks:
The usual, marginally unbiased backtransformation for lognormal 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 backtransformed 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
logtransformed response Z(s)
or the estimated trend
x(s)^T
hatβ of the logtransformed data can be backtransformed (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.
Block Kriging predictions of a logtransformed response variable are backtransformed 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 backtransformed 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 backtransformation is based on the assumption that both the point data U(s) and the block means U(A) 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(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
.
lgnpp
allows as a further option to backtransform 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):
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 U(A).
The mean squared error is then approximated by
E[{U(A)  hatU(A)}^2] = 1/K^2 sum_{s_i in A} E[ { U(s_i)  hatU(s_i)}^2]
+ (K1)/(K k (k1)) sum_{s_i in sample} sum_{s_j in sample, s_j != s_i} Cov[ U(s_i)hatU(s_i), U(s_j)hatU(s_j)]
The first term of the RHS (and hatU(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
hatU(A) = sum_{s_i in A} hatU(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[ { U(s_i)  hatU(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.
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:
lgn.pred
: the backtransformed Kriging predictions of a
logtransformed response.
lgn.se
: the standard errors of the
backtransformed predictions.
lgn.lower
, lgn.upper
: the bounds of the
backtransformed prediction intervals.
If is.block
is TRUE
or all.pred
not equal to
NULL
a named numeric vector with two elements:
mean
: the backtransformed 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 backtransformed point prediction errors.
Andreas Papritz [email protected].
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://ecollection.library.ethz.ch/eserv/eth:6027/eth602701.pdf
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.
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)
## backtransformation 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)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.