predict.geolm_cmodMan {gear}

Predict method for geostatistical models

Description

predict calculates the predicted values at specified locations. The method can additionally provide the mean square prediction error (mspe) and perform conditional simulation.

Usage

## S3 method for class 'geolm_cmodMan'
predict(
  object,
  newdata,
  nsim = 0,
  vop,
  vp,
  return_type = "SpatialPointsDataFrame",
  dmethod = "chol",
  compute_mspe = TRUE,
  sp = NULL,
  ...
)

Arguments

object

An object produced by the geolm function.

newdata

An optional data frame in which to look for the coordinates at which to predict. If omitted, the observed data locations are used.

nsim

A non-negative integer indicating the number of realizations to sample at the specified coordinates using conditional simulation.

vop

The cross-covariance matrix between the observed responses and the responses to predict.

vp

The covariance matrix of the responses to predict.

return_type

A character string indicating the type of object that should be returned. The default is "SpatialPointsDataFrame" for easy plotting of results (see Examples). Other options include "data.frame", "geardf", and "sf".

dmethod

The method used to decompose the covariance matrix for conditional simulation. Valid options are "chol", "eigen", and "svd". The default is "chol".

compute_mspe

A logical value indicating whether the mean square prediction error should be calculated. Default is TRUE.

sp

This argument will be deprecated in the future. Please use the return_type argument. A logical value indicating whether to object returned should be of class SpatialPointsDataFrame for easier plotting with the sp package. Default is NULL.

...

Currently unimplemented.

Details

The newdata data frame must include the relevant covariates for the prediction locations, where the covariates are specified on the right side of the ~ in object$formula. newdata must also include the coordinates of the prediction locations, with these columns having the names provided in object$coordnames.

Value

A data.frame, SpatialPointsDataFrame, geardf, or sf object with the kriging predictions pred, kriging variance/mean-square prediction error (mspe), the root mean-square prediction error mspe (rmspe), and the conditional simulations sim.1, sim.2, etc. sim.1, sim.2, etc.

Author(s)

Joshua French

Examples

# generate response
y = rnorm(10)
# generate coordinates
x1 = runif(10); x2 = runif(10)

# data frame for observed data
data = data.frame(y, x1, x2)
coords = cbind(x1, x2)
d = as.matrix(dist(coords))
psill = 2 # partial sill
r = 4 # range parameter
evar = .1 # error variance
fvar = .1 # add finescale variance
# one can't generally distinguish between evar and fvar, but
# this is done for illustration purposes

# manually specify an exponential covariance model
v = psill * exp(-d/r) + (evar + fvar) * diag(10)
mod_man = cmod_man(v = v, evar = evar)

# coordinate names
cnames = c("x1", "x2")

# geolm for universal kriging
gearmod_uk = geolm(y ~ x1 + x2, data = data, mod = mod_man,
                 coordnames = cnames)

# newdata must have columns with prediction coordinates
# add 5 unsampled sites to sampled sites
newdata = data.frame(x1 = c(x1, runif(5)), x2 = c(x2, runif(5)))
newcoords = newdata[, cnames]
# create vop and vp using distances
dop = geodist(as.matrix(coords), as.matrix(newcoords))
dp = geodist(newcoords)

# manually create cross-covariance and covariance for
# prediction locations
vop = psill * exp(-dop/r) + fvar * (dop == 0)
vp = psill * exp(-dp/r) + fvar * diag(nrow(newcoords))

# prediction for universal kriging, with conditional simulation,
# using manual covariance matrices
pred_uk_man = predict(gearmod_uk, newdata, nsim = 2,
                      vop = vop, vp = vp, dmethod = "svd")

# do the same thing, but using cmod_std

# prediction for universal kriging, with conditional simulation
mod_std = cmod_std("exponential", psill = psill, r = r,
                   evar = evar, fvar = fvar)
gearmod_uk2 = geolm(y ~ x1 + x2, data = data, mod = mod_std,
                    coordnames = c("x1", "x2"))
pred_uk_std = predict(gearmod_uk2, newdata, nsim = 2, dmethod = "svd")

# compare results
all.equal(pred_uk_man$pred, pred_uk_std$pred)
all.equal(pred_uk_man$mspe, pred_uk_std$mspe)