intkrige {intkrige}

Algorithmic implementation of interval valued kriging.

Description

Function to implement the interval valued extension of ordinary and simple kriging. Includes all necessary input checks and error handling. Essentially acts as a switch function between the R and c++ versions of the algorithm.

Usage

intkrige(
  locations,
  newdata,
  models,
  formulas = list(center ~ 1, radius ~ 1),
  eta = 0.75,
  A = c(1, 1, 0),
  trend = NULL,
  thresh = 100,
  tolq = 0.001,
  maxq = 100,
  tolp = 0.001,
  maxp = 100,
  r = 1,
  useR = TRUE,
  fast = FALSE,
  weights = FALSE,
  cores = 1
)

Arguments

locations

An object of class intsp, specifying the prediction locations with an interval-valued response.

newdata

An object of class SpatialPointsDataFrame or SpatialPixelsDataFrame specifying the locations at which to predict intervals.

models

A list of variogram models of class vgm (see vgm) When specified, the third model represents the center/radius interaction.

formulas

A list of two formulas specifying the centering and scaling of the center and radius respectively. Limitations to these formulas are specified in the details.

eta

A growth/shrink parameter for penalty term. For simple kriging: eta > 1. For ordinary kriging: eta < 1.

A

vector of length three representing the weights of the generalized L2 distance: the vector of three contains the weights for the center, radius, and center/radius respectively. A = c(1, 1, 0) assumes the regular L2 distance calculation for intervals.

trend

If null, use ordinary kriging. When specified, represents the known mean of the stationary process and is an indication to use simple kriging.

thresh

Let n = length(locations). When abs(lam_i) < 1/(n*thresh), this lambda value is set to 0.

tolq

For a set penalty term, convergence is satisfied if max(abs(lamUp-lam)) < tolq.

maxq

For a set penalty term, the max number of iterations allowed for convergence.

tolp

When abs(sum(abs(lam)) - 1) < tolp, consider the constraints satisfied.

maxp

Maximum number of allowed iterations to satisfy equation constraints.

r

The starting value of the penalty term. Must be relatively large to ensure that the initial solution stays within the feasible region.

useR

If TRUE, use the R version of the algorithm. If FALSE, use the rcppArmadillo version.

fast

(Simple kriging only). If TRUE, allows lambdas to converge to 0 and subsets matrices accordingly. When FALSE, runs simple kriging using a barrier penalty at 0. Fast = TRUE is orders of magnitude faster than the full implementation. However, it is not recommended when input measurements are sparse as it is known to have convergence issues in these cases.

weights

If TRUE, return the vector kriging weights for each prediction. If false, simply return the predicted output.

cores

An integer (for parallel computing): specify the number of cores that will be devoted to the computation. Note that the argument 'all' will use all available cores minus one. Parallel processing is only relevant if you are predicting for more than one location. Note there is no parallel option when useR = FALSE.

Details

The formulas argument is current fairly limited in its use. For example the center argument can accept no transformations of the dependent variable. Similarly, the radius argument can accept no variable arguments as independent variables. The idea behind this limited use of formulas is that any transformations should be applied to the entire interval prior to input into interval-valued kriging. This ensures that input into the interval-valued kriging algorithm are well-defined intervals with properly ordered upper and lower endpoints. The transformation that are allowed within this function are linear shifts of the center, and linear scaling of the radius. Note that the scaling term for the radius can contain a division parameter but it must be encapsulated in parenthesis and included on the right hand side of the multiplication parameter.

Value

A matrix with 4 columns where rows correspond to the prediction locations and columns correspond to:

- center prediction

- radius predictions

- kriging prediction variance

- warning column for non-convergent optimization problem (0 - no warning, 1 - warning)

Examples

# First, define the location and elevation of interest.
# (In this case we pick coordinates of Utah State University)
templocs <- data.frame(lat = 41.745, long = -111.810, ELEVATION = 1456)
sp::coordinates(templocs) <- c("long", "lat")
sp::proj4string(templocs) <- "+proj=longlat +ellps=WGS84
+datum=WGS84 +no_defs +towgs84=0,0,0"

# Load the Utah Snow Load Data
data(utsnow)
utsnow.sp <- utsnow

# Convert to an 'intsp' object that inherits a SpatialPointsDataFrame
sp::coordinates(utsnow.sp) <- c("LONGITUDE", "LATITUDE")
sp::proj4string(utsnow.sp) <- sp::proj4string(templocs)
interval(utsnow.sp) <- c("minDL", "maxDL")
# Define the formulas we will use to define the intervals.
temp_formulas <- list(center ~ ELEVATION,
                     radius*(ELEVATION/median(ELEVATION)) ~ 1)

# Define, fit and check the variogram fits.
varios <- intvariogram(utsnow.sp,
                       formulas = temp_formulas)
varioFit <- fit.intvariogram(varios, models = gstat::vgm(c("Sph", "Sph", "Gau")))
preds <- intkrige::intkrige(locations = utsnow.sp,
newdata = templocs,
models = varioFit,
formulas = temp_formulas)