SpatialTools (version 0.4.1)

krige.uk: Performs universal kriging

Description

Performs Universal Kriging using X, the $n \times k$ design matrix for the regression coefficients of the observed data, y, the $n \times 1$ matrix of observed responses, V, the (positive definite) covariance matrix of the observed responses, Xp, the $np \times k$ design matrix of the responses to be predicted, Vp, the $np \times np$ covariance matrix of the responses to be predicted, and Vop, the $n \times np$ matrix of covariances between the observed responses and the responses to be predicted. Uses Armadillo C++ template via RcppArmadillo to perform most of the operations.

Usage

krige.uk(y, V, Vp, Vop, X, Xp, return.w = FALSE, nsim = 0, 
	Ve.diag = NULL, method = "eigen")

Arguments

y

The vector of observed responses. Should be a matrix of size $n \times 1$ or a vector of length $n$.

V

The covariance matrix of the observed responses. The size is $n \times n$.

Vp

The covariance matrix of the responses to be predicted. The size is $np \times np$

Vop

The cross-covariance between the observed responses and the responses to be predicted. The size is $n \times np$

X

The design matrix of the observed data. The size is $n \times k$

Xp

The design matrix of the responses to be predicted. The size is $np \times k$

return.w

A boolean value (TRUE or FALSE) indicating whether the matrix of prediction weights should be returned. By default, this is FALSE (and results in faster computation since there is less to return).

nsim

A non-negative value indicating the number of conditional simulations that should be returned. If this is less than 1, then no conditional simulation is done.

Ve.diag

A vector of length n specifying the measurement error variances of the observed data.

method

A character vector specifying the method used to decompose V. Options are "eigen", "chol", or "svd" (Eigen decomposition, Cholesky decomposition, or Singular value decomposition, respectively).

Value

• The function a list containing the following objects:
• pred
  A vector of length $np$ containing the predicted responses.
• mspe
  A vector of length $np$ containing the mean-square prediction error of the predicted responses.
• coeff
  A vector of length $k$ containing the estimated regression coefficients.
• vcov.coeff
  A $k \times k$ matrix containing the (estimated) covariance matrix of estimated the regression coefficients.
• w
  A $np \times n$ matrix containing the kriging weights used to calculate pred.
• simulations
  An $n \times nsim$ matrix containing the nsim realizations of the conditional realizations. Each column of the matrix represents a realization of the conditional normal distribution.

Details

It is assumed that there are $n$ observed data values and that we wish to make predictions at $np$ locations. We assume that there are $k$ regression coefficients (including the intercept). Both X and Xp should contain a column of 1's if an intercept is desired.

If doing conditional simulation, the Cholesky decomposition should not work when there are coincident locations between the observed data locations and the predicted data locations. Both the Eigen and Singular Value Decompositions should work.

References

Statistical Methods for Spatial Data Analysis, Schabenberger and Gotway (2003). See p. 241-243.

Examples

# create observed and predicted coordinates
ocoords <- matrix(runif(100), ncol = 2)
pcoords <- matrix(runif(200), ncol = 2)

# include some observed locations in the predicted coordinates
acoords <- rbind(ocoords, pcoords)

# create design matrices
X <- as.matrix(cbind(1, ocoords))
Xa <- as.matrix(cbind(1, acoords))

# create covariance matrix
C3 <- cov.sp(coords = ocoords, sp.type = "matern", sp.par = c(2, 1), smoothness = 1, 
	finescale = 0, error = 0.5, pcoords = acoords)

# set values of regression coefficients
coeff <- matrix(c(1, 2, 3), nrow = 1)

# generate data with error
y <- rmvnorm(nsim = 1, mu = tcrossprod(X, coeff), V = C3$V) + rnorm(50, sd = sqrt(.5))

# use universal kriging to make predictions.  Do not return weights matrix or do 
# conditional simulation
krige.obj <- krige.uk(as.vector(y), V = C3$V, Vp = C3$Vp, Vop = C3$Vop, 
	X = X, Xp = Xa, return.w = FALSE, nsim = 0)

#Do return weights matrix and do conditional simulation
krige.obj2 <- krige.uk(as.vector(y), V = C3$V, Vp = C3$Vp, Vop = C3$Vop, 
	X = X, Xp = Xa, return.w = TRUE, nsim = 100, 
	Ve.diag = rep(.5, 50), method = "eigen")