intkrige: Algorithmic implementation of interval valued kriging.
In intkrige: A Numerical Implementation of Interval-Valued Kriging

Description Usage Arguments Details Value Examples

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.

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
)

`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.

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.

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)

# 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)