CKrige: Constrained, Covariance-matching Constrained and Universal Kriging

Description

Function computes constrained, covariance-matching constrained and universal (external drift) Kriging predictions for points or blocks (of any shape) in a global neighbourhood and for isotropic covariance models.

Usage

CKrige(formula, data, locations, object, ...)
# S4 method for formula,data.frame,formula,preCKrigePolygons
CKrige(formula,
  data, locations, object, method = 2, ex.out = FALSE, ncores = 1L,
  fork = !identical(.Platform[["OS.type"]], "windows"))
# S4 method for formula,data.frame,formula,preCKrigePoints
CKrige(formula,
  data, locations, object, method = 2, ex.out = FALSE)

Value

If ex.out = FALSE then the function CKrige

returns an object of class “SpatialPointsDataFrame” or “SpatialPolygonsDataFrame”. The class of the argument object of CKrige
(“preCKrigePoints” or “preCKrigePolygons”) determines whether a SpatialPointsDataFrame or SpatialPolygonsDataFrame

is returned.

Irrespective of the chosen Kriging method, the data frame in the data slot of the returned object contains the columns:

prediction: a numeric vector with the Kriging predictions by the chosen method.
prediction.se: a numeric vector with the root mean square prediction errors (Kriging standard errors).

For constrained Kriging (argument method = 2) the data frame has 3 additional columns, see

constrainedKriging-package:

P1: a numeric vector with $P_{1} = (\mathrm{Var}[Y(B_1)] - \mathrm{Var}[\mathbf{x}(B_1)^{\prime}\widehat{\boldsymbol{\boldsymbol{\beta}}}_{\mathrm{GLS}}] )^{0.5}. $
Q1: a numeric vector with $Q_{1} = (\mathrm{Var}[\widehat{Y}_{\mathrm{UK}}(B_1)] - \mathrm{Var}[ \mathbf{x}(B_1)^{\prime}\widehat{\boldsymbol{\boldsymbol{\beta}}}_{\mathrm{GLS}}] )^{0.5}. $
K: a numeric vector with $K = P_{1} / Q_{1}$.

For covariance-matching constrained Kriging (argument method = 3) the data frame has also 3 additional columns, see constrainedKriging-package:

P1.11: a numeric vector with first diagonal elements (which correspond to the target blocks) of the matrices $$ \mathbf{P}_{1} = ( \mathrm{Cov}[\mathbf{Y}, \mathbf{Y}^{\prime}] - \mathrm{Cov}[\mathbf{X}_{m}\widehat{\boldsymbol{\boldsymbol{\beta}}}_{\mathrm{GLS}}, (\mathbf{X}_{m}\widehat{\boldsymbol{\boldsymbol{\beta}}}_{\mathrm{GLS}})^{\prime}] )^{\frac{1}{2}}. $$
Q1.11: a numeric vector with first diagonal elements of the matrices $$ \mathbf{Q}_{1} = (\mathrm{Cov}[\widehat{\mathbf{Y}}_{\mathrm{UK}}, \widehat{\mathbf{Y}}_{\mathrm{UK}}^{\prime}] - \mathrm{Cov}[\mathbf{X}_{m}\widehat{\boldsymbol{\boldsymbol{\beta}}}_{\mathrm{GLS}}, (\mathbf{X}_{m}\widehat{\boldsymbol{\boldsymbol{\beta}}}_{\mathrm{GLS}})^{\prime}])^{\frac{1}{2}}. $$
K.11: a numeric vector with first diagonal elements of the matrices of the matrices $$ \mathbf{K} = \mathbf{Q}_{1}^{-1}\mathbf{P_{1}}. $$

If the argument ex.out = TRUE is used and the argument method is either equal to 1 (universal) or 2

(constrained Kriging) then the function CKrige returns an object either of class

“CKrige.exout.polygons” or “CKrige.exout.points”, which are lists with the following components:

object: either an object of class “SpatialPointsDataFrame” or
“SpatialPolygonsDataFrame” as described above.
krig.method: a numeric scalar, number of the chosen Kriging method (1, 2).
parameter: a list with 2 components. First component beta.coef is the vector with the generalized least squares estimates of the linear regression coefficients, and the second component cov.beta contains the covariance matrix of the generalized least squares estimates.
sk.weights: a matrix with the simple Kriging weights. The ith column contains the simple Kriging weights for the ith prediction target.
inv.Sigma: a matrix with the inverse covariance matrix $\boldsymbol{\Sigma}^{-1}$ of the observations.
residuals: a numeric vector with the generalized least squares residuals of the linear regression.

If the argument ex.out = TRUE is used and the argument method is either equal to 3 (covariance-matching constrained kriging) then the function CKrige returns an object either of class

“CKrige.exout.polygons” or “CKrige.exout.points”, which are lists with the following components:

object

either an object of class “SpatialPointsDataFrame” or
“SpatialPolygonsDataFrame” as described above.

krig.method

a numeric scalar, number of the chosen Kriging method (3).

CMCK.par

a list of 3 lists with the following components:

P1 a list of matrices $\mathbf{P}_1$ for all point or polygon neighbourhood configurations (the first row and columns correspond to the target points or blocks of the configurations).
Q1 a list of matrices $\mathbf{Q}_1$ for all point or polygon neighbourhood configurations (the first row and columns correspond to the target points or blocks of the configurations).
K a list of matrices $\mathbf{K}$ for all polygon neighbourhood configurations (the first row and columns correspond to the target points or blocks of the configurations).

parameter

a list with 2 components. First component beta.coef is the vector with the generalized least squares estimates of the linear regression coefficients, and the second component cov.beta contains the covariance matrix of the generalized least squares estimates.

sk.weights

a list with the simple Kriging weights of all point or polygon neighbourhood configurations (the first columns contain the weights of the target points or blocks of the configurations).

inv.Sigma

a matrix with the inverse covariance matrix $\boldsymbol{\Sigma}^{-1}$ of the observations.

residuals

a numeric vector with the generalized least squares residuals of the linear regression.

Arguments

formula: a formula for the linear regression model of the spatial variable in the form response ~ terms of covariates, for ordinary Kriging use the formula
response ~ 1.
data: a data frame with the values of the covariates for the observations. The names of the covariates used in formula must match the column names of data.
locations: a one-sided formula that describes the coordinates of the observations (e.g. ~ x+y).
object: either an object of class “preCKrigePolygons” for block Kriging or of class “preCKrigePoints” for point Kriging as generated by the preCKrige function.
...: further arguments to control the spatial interpolation method and the output.
method: a numeric scalar to choose the Kriging approach: method = 1: universal (external drift), method = 2: constrained (default) and method = 3: covariance-matching constrained Kriging.
ex.out: a logical scalar. If ex.out is set equal to TRUE CKrige returns an extended output with additional items, see Value for more information, by default ex.out = FALSE.
ncores: a positive integer scalar with the number of CPUs to use for parallel computations.
fork: a logical scalar to control whether parallel computations are done by forking using mclapply (non-windows OSes) or by socket clusters using parLapply (windows OS).

Author

Christoph Hofer, christoph.hofer@alumni.ethz.ch

Details

The CKrige function depends on a preCKrige output object that contains the parameters of the isotropic covariance model, the covariates of the prediction targets and the variance-covariance matrices of the prediction targets.

References

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

Hofer, C. and Papritz, A. (2011). constrainedKriging: an R-package for customary, constrained and covariance-matching constrained point or block Kriging. Computers & Geosciences. 37, 1562--1569, tools:::Rd_expr_doi("10.1016/j.cageo.2011.02.009").

Examples

Run this code

# \donttest{
### load meuse data
data(meuse, package = "sp")
data(meuse.grid, package = "sp")
data(meuse.blocks)

### convert data frame meuse.grid to SpatialPointsDataFrame
coordinates(meuse.grid) <- ~ x+y

### plot blocks along with observation locations
plot(meuse.blocks)
points(y ~ x, meuse, cex = 0.5)


### example 1: lognormal constrained block Kriging
### cf. Hofer & Papritz, 2011, pp. 1567

### compute the approximated block variance of each block in meuse.blocks
### without any neighbouring blocks (default, required for in universal
### and constrained Kriging) for an exponential covariance function without
### a measurement error, a nugget = 0.15 (micro scale white noise process),
### a partial sill variance = 0.15 and a scale parameter = 192.5
### approximation of block variance by pixels of size 75m x 75m
preCK_block <- preCKrige(newdata = meuse.blocks, model = covmodel(modelname =
    "exponential", mev = 0, nugget = 0.05, variance = 0.15,
    scale = 192.5), pwidth = 75, pheight = 75)

### block prediction by constrained Kriging on the log scale
### extended output required for backtransformation
CK_block <- CKrige(formula = log(zinc) ~ sqrt(dist), data = meuse,
  locations = ~ x+y, object = preCK_block, ex.out = TRUE)

### backtransformation of the CK predictions to original scale
### cf Hofer & Papritz, 2011, eq. (14), p. 1567
### note that Var(\hat{Y}{B}) = Var(Y(B))!
beta <- CK_block$parameter$beta.coef
M <- meuse.blocks@data$M
var.blockmeans <- unlist(preCK_block@covmat)
CK_block$object@data$Zn <- exp(CK_block$object@data$prediction
  + 0.5*(0.2 + beta[2]^2 * M - var.blockmeans))

### upper limits U of the relative 95%  prediction intervals for CK
### U multiplied by the predictions CK_block$object@data$Zn gives
### the upper limits of the 95% prediction intervals
CK_block$object@data$U <- exp(CK_block$object@data$prediction
    + 1.96 * CK_block$object@data$prediction.se) / CK_block$object@data$Zn

### plot the CK predictions of Zn and the upper limits of relative 95%
### prediction intervals by function spplot of sp package
### function ck.colors(n) creates a rainbow-like color vector
breaks <- seq(0, 1850, by = 185)
spplot(CK_block$object, zcol = "Zn", at = breaks, col.regions = ck.colors(10),
  colorkey = list(labels = list(at = breaks, labels = breaks)),
  main = "CK predictions of Zn block means")

### plot of upper limits of the relative 95% prediction intervals
breaks <- seq(1, 3.2, by = 0.2)
spplot(CK_block$object, zcol = "U", at = breaks, col.regions = ck.colors(11),
  colorkey = list(labels = list(at = breaks, labels = breaks)),
  main = "Upper limits rel. prediction intervals CK predictions")



### example 2: lognormal covariance-matching constrained block Kriging

### define neighbours by using the poly2nb function of spdep package
if(!requireNamespace("spdep", quietly = TRUE)){
  stop("install package spdep to run example")
}
neighbours_block <- spdep::poly2nb(meuse.blocks)
class(neighbours_block)
### neighbours_block should be an object of class "list"
class(neighbours_block) <- "list"

### compute the approximated block variance-covariance matrices of each
### polygon neighbourhood configuration (= target block plus its
### neighbours)
preCMCK_block <- preCKrige(newdata = meuse.blocks, neighbours = neighbours_block,
  model = covmodel("exponential", mev = 0, nugget = 0.05, variance = 0.15,
    scale = 192.5), pwidth = 75, pheight = 75)

### block prediction by covariance-matching constrained Kriging on log
### scale
CMCK_block <- CKrige(formula = log(zinc) ~ sqrt(dist), data = meuse,
  locations = ~ x+y, object = preCMCK_block, method = 3, ex.out = TRUE)

### backtransformation of the CK predictions to the original scale
### cf Hofer & Papritz, 2011, eq. (14), p. 1567
beta <- CMCK_block$parameter$beta.coef
M <- meuse.blocks@data$M
var.blockmeans <- sapply(preCMCK_block@covmat, function(x) x[1, 1])
CMCK_block$object@data$Zn <- exp(CMCK_block$object@data$prediction
  + 0.5*(0.2 + beta[2]^2 * M - var.blockmeans))

### plot the CMCK predictions of Zn by function spplot of sp package
### function ck.colors(n) creates a rainbow-like color vector
breaks <- seq(0, 1850, by = 185)
spplot(CMCK_block$object, zcol = "Zn", at = breaks, col.regions = ck.colors(10),
  colorkey = list(labels = list(at = breaks, labels = breaks)),
  main = "CMCK predictions of Zn block means")



### example 3: lognormal universal block Kriging

### block prediction by universal Kriging on log scale
UK_block <- CKrige(formula = log(zinc) ~ sqrt(dist), data = meuse,
  locations = ~ x+y, object = preCK_block, method = 1, ex.out = TRUE)

### backtransformation of the CK predictions to the original scale
### cf Hofer & Papritz, 2011, Appendix B, pp. 1568 - 1569
beta <- UK_block$parameter$beta.coef
cov.beta <- UK_block$parameter$cov.beta.coef
M <- meuse.blocks@data$M
var.blockmeans <- unlist(preCK_block@covmat)
SKw <-  UK_block$sk.weights
X <-  model.matrix(~sqrt(dist), meuse)
XB <- model.matrix(~sqrt(dist), meuse.blocks)
c1 <- 0.5 * (0.2 + beta[2]^2*M - var.blockmeans +
  UK_block$object@data$prediction.se^2)
c2 <- rowSums((XB - t(SKw) %*% X) * (XB %*% cov.beta))
UK_block$object@data$Zn <- exp(UK_block$object@data$prediction + c1 - c2)

### upper limits U of the relative 95%  prediction intervals for CK
### U multiplied by the predictions UK_block$object@data$Zn gives
### the upper limits of the 95% prediction intervals
UK_block$object@data$U <- exp(UK_block$object@data$prediction
    + 1.96 * UK_block$object@data$prediction.se) / UK_block$object@data$Zn

### plot the UK predictions of Zn by function spplot of sp package
### function ck.colors(n) creates a rainbow-like color vector
breaks <- seq(0, 1850, by = 185)
spplot(UK_block$object, zcol = "Zn", at = breaks, col.regions = ck.colors(10),
  colorkey = list(labels = list(at = breaks, labels = breaks)),
  main = "UK predictions of Zn block means")

### plot of upper limits of the relative 95% prediction intervals
breaks <- seq(1, 3.2, by = 0.2)
spplot(UK_block$object, zcol = "U", at = breaks, col.regions = ck.colors(11),
  colorkey = list(labels = list(at = breaks, labels = breaks)),
  main = "Upper limits rel. prediction intervals UK predictions")



### example 4: constrained point Kriging
### generate point CK preCKrige object for locations in meuse.grid
preCK_point <- preCKrige(newdata = meuse.grid, model = covmodel(modelname =
    "exponential", mev = 0, nugget = 0.05, variance = 0.15,
    scale = 192.5))

### point prediction by constrained Kriging on the log scale
### no extended output required for backtransformation
CK_point <- CKrige(formula = log(zinc) ~ sqrt(dist), data = meuse,
  locations = ~ x+y, object = preCK_point)

### backtransformation of the CK predictions to the original scale
### cf. Cressie, 2006, eq. (20), p. 421
### note that Var(\hat{Y}{s_0}) = Var(Y(s_0))!
CK_point@data$Zn <- exp(CK_point@data$prediction)

### convert results object to SpatialGridDataFrame
gridded(CK_point) <- TRUE
fullgrid(CK_point) <- TRUE

### plot the CK predictions of Zn by function spplot of sp package
breaks <- seq(0, 2600, by = 185)
spplot(CK_point, zcol = "Zn", at = breaks, col.regions = ck.colors(20),
  colorkey = list(labels = list(at = breaks, labels = breaks)),
  main = "CK predictions of Zn point values")



### example 5: covariance-matching constrained point Kriging
### define 4 nearest neighbours to each grid location by using the
### knearneigh function of the spdep package
if(!requireNamespace("spdep", quietly = TRUE)){
  stop("install package spdep to run example")
}
neighbours_point <- spdep::knearneigh(meuse.grid, k = 4)
### convert matrix with neighbours indices to list
neighbours_point <- apply(neighbours_point$nn, 1, FUN = function(x) x,
  simplify = FALSE)

### generate point CMCK preCKrige object for locations in meuse.grid
preCMCK_point <- preCKrige(newdata = meuse.grid, neighbours = neighbours_point,
  model = covmodel(modelname = "exponential", mev = 0, nugget = 0.05,
    variance = 0.15, scale = 192.5))

### point prediction by covariance-matching constrained Kriging on the log
### scale
### no extended output required for backtransformation
CMCK_point <- CKrige(formula = log(zinc) ~ sqrt(dist), data = meuse,
  locations = ~ x+y, object = preCMCK_point, method = 3)

### backtransformation of the CMCK predictions to the original scale
### cf. Cressie, 2006, eq. (20), p. 421
### note that Var(\hat{Y}{s_0}) = Var(Y(s_0))!
CMCK_point@data$Zn <- exp(CMCK_point@data$prediction)

### convert results object to SpatialGridDataFrame
gridded(CMCK_point) <- TRUE
fullgrid(CMCK_point) <- TRUE

### plot the CMCK predictions of Zn by function spplot of sp package
breaks <- seq(0, 2600, by = 185)
spplot(CMCK_point, zcol = "Zn", at = breaks, col.regions = ck.colors(20),
  colorkey = list(labels = list(at = breaks, labels = breaks)),
  main = "CMCK predictions of Zn point values")



### example 6: universal point Kriging
UK_point <- CKrige(formula = log(zinc) ~ sqrt(dist), data = meuse,
  locations = ~ x+y, object = preCK_point, method = 1, ex.out = TRUE)

### backtransformation of the UK predictions to the original scale cf.
### Cressie, 2006, eq.  (20), p.  421 and Hofer & Papritz, 2011, Appendix
### B, p.  1569
cov.beta <- UK_point$parameter$cov.beta.coef
SKw <- UK_point$sk.weights
X <-  model.matrix(~sqrt(dist), meuse)
Xgrid <- model.matrix(~sqrt(dist), meuse.grid)
### universal kriging variances
c1 <- 0.5 * UK_point$object@data$prediction.se^2
### \psi^' x(s_0)
c2 <- rowSums((Xgrid - t(SKw) %*% X) * (Xgrid %*% cov.beta))
UK_point$object@data$Zn <- exp(UK_point$object@data$prediction + c1 - c2)

### convert results object to SpatialGridDataFrame
gridded(UK_point$object) <- TRUE
fullgrid(UK_point$object) <- TRUE

### plot the CMCK predictions of Zn by function spplot of sp package
breaks <- seq(0, 2600, by = 185)
spplot(UK_point$object, zcol = "Zn", at = breaks, col.regions = ck.colors(20),
  colorkey = list(labels = list(at = breaks, labels = breaks)),
  main = "UK predictions of Zn point values")
# }



### example 7: universal block Kriging of mean over whole domain
### cf chapter 4 of vignette of package georob
### (https://CRAN.R-project.org/package=georob)
if(!requireNamespace("gstat", quietly = TRUE)){
  stop("install package gstat to run example")
}
### load coalash data
data(coalash, package="gstat")

### generate SpatialPointsDataFrame covering entire domain of coalash data
coalash.domain <- rbind(c(0.5,0), c(16.5,0), c(16.5,24), c(0.5,24), c(0.5,0))
coalash.domain <- SpatialPolygonsDataFrame(
  SpatialPolygons(list(Polygons(list(Polygon(coalash.domain)), ID= "domain"))),
  data=data.frame(x=8.5,y=12,row.names="domain"))

### universal block Kriging
preCK_coalash <- preCKrige(newdata = coalash.domain, model = covmodel(modelname =
    "exponential", mev = 1.023, nugget = 0, variance = 0.268,
    scale = 1.907), pwidth = 16, pheight = 24)
UK_coalash <- CKrige(formula = coalash ~ x, data = coalash,
  locations = ~ x+y, object = preCK_coalash, method = 1)
slot( UK_coalash, "data")

Run the code above in your browser using DataLab