predict.georob {georob}

Predict Method for Robustly Fitted Spatial Linear Models

Description

Robust and customary external drift Kriging prediction based on a spatial linear models fitted by georob. The predict method for the class georob computes fitted values, point and block Kriging predictions as well as model terms for display by termplot.

Usage

## S3 method for class 'georob'
predict(object, newdata, type =  c("signal", "response", "trend", "terms"),
    terms = NULL, se.fit = TRUE, signif = 0.95, locations,
    variogram.model = NULL, param = NULL, aniso = NULL, variogram.object = NULL,
    control = control.predict.georob(), verbose = 0, ...)

control.predict.georob(full.covmat = FALSE, extended.output = FALSE,
    mmax = 10000, ncores = pcmp[["max.ncores"]], pwidth = NULL, pheight = NULL,
    napp = 1, pcmp = control.pcmp())

Arguments

object

an object of class "georob" (mandatory argument), see georobObject.

newdata

an optional data frame, SpatialPointsDataFrame, SpatialPixelsDataFrame, SpatialGridDataFrame, SpatialPolygonsDataFrame or an (optional) object of class SpatialPoints, SpatialPixels or SpatialGrid, in which to look for variables with which to compute fitted values or Kriging predictions, see Details.
If newdata is a SpatialPolygonsDataFrame then block Kriging predictions are computed, otherwise point Kriging predictions.

type

a character keyword defining what target quantity should be predicted (computed). Possible values are

• "signal": the “signal” Z(\boldsymbol{s}) = \boldsymbol{x}(\boldsymbol{s})^\mathrm{T} \boldsymbol{\beta} + B(\boldsymbol{s}) of the process (default),

• "response": the observations Y(\boldsymbol{s}) = Z(\boldsymbol{s}) + \varepsilon(\boldsymbol{s}),

• "trend": the external drift \boldsymbol{x}(\boldsymbol{s})^\mathrm{T} \boldsymbol{\beta},

• "terms": the model terms.

terms

If type = "terms", which terms (default is all terms).

se.fit

a logical scalar, only used if type is equal to "terms", see predict.lm.

signif

a positive numeric scalar equal to the tolerance or confidence level for computing respective intervals. If NULL no intervals are computed.

locations

an optional one-sided formula specifying what variables of newdata are the coordinates of the prediction points
(default: object[["locations.objects"]][["locations"]]).

variogram.model

an optional character keyword defining the variogram model to be used for Kriging, see georob and Details.

param

an optional named numeric vector with values of the variogram parameters used for Kriging, see georob and Details.

aniso

an optional named numeric vector with values of anisotropy parameters of a variogram used for Kriging, see georob and Details.

variogram.object

an optional list that defines a possibly nested variogram model used for Kriging, see georob and Details.

control

a list with the components full.covmat, extended.output, mmax, ncores, pwidth, pheight, napp and pcmp or a function such as control.predict.georob that generates such a list.

full.covmat

a logical scalar controlling whether the full covariance matrix of the prediction errors is returned (TRUE) or only the vector with its diagonal elements (FALSE, default), see Value for an explanation of the effect of full.covmat.

extended.output

a logical scalar controlling whether the covariance matrices of the Kriging predictions and of the data should be computed, see Details (default FALSE).

mmax

a positive integer equal to the maximum number (default 10000) of prediction items, computed in a sub-task in parallelized computations, see Details.

ncores

a positive integer controlling how many cores are used for parallelized computations, see Details.

pwidth, pheight, napp

numeric scalars, used to tune numeric integration of semi-variances for block Kriging, see preCKrige.

pcmp

a list of arguments passed to pmm and other functions that carry out parallelized computations or a function such as control.pcmp that generates such a list (see control.pcmp for allowed arguments).

verbose

a positive integer controlling logging of diagnostic messages to the console.
verbose = 0 (default) largely suppresses such messages.

...

arguments passed to control.predict.georob.

Details

If newdata is an object of class SpatialPoints, SpatialPixels or SpatialGrid then the drift model may only use the coordinates as covariates (universal Kriging). Furthermore the names used for the coordinates in newdata must be the same as in data when creating object (argument locations of predict.georob should not be used). Note that the result returned by predict.georob is then an object of class SpatialPointsDataFrame, SpatialPixelsDataFrame or SpatialGridDataFrame.

The predict method for class georob uses the packages parallel and snowfall for parallelized computation of Kriging predictions. If there are m items to predict, the task is split into ceiling(m/mmax) sub-tasks that are then distributed to ncores CPUs. Evidently, ncores = 1 suppresses parallel execution. By default, the function uses all available CPUs as returned by detectCores.
Note that if full.covmat is TRUE mmax must exceed m (and parallel execution is not possible).

The argument extended.output = TRUE is used to compute all quantities required for (approximately) unbiased back-transformation of Kriging predictions of log-transformed data to the original scale of the measurements by lgnpp. In more detail, the following items are computed:

• trend: the fitted values, \boldsymbol{x}(\boldsymbol{s})\mathrm{^T}\widehat{\boldsymbol{\beta}},

• var.pred: the variances of the Kriging predictions, \mathrm{Var}_{\hat{\theta}}[\widehat{Y}(\boldsymbol{s})] or \mathrm{Var}_{\hat{\theta}}[\widehat{Z}(\boldsymbol{s})],

• cov.pred.target: the covariances between the predictions and the prediction targets,
  \mathrm{Cov}_{\hat{\theta}}[\widehat{Y}(\boldsymbol{s}),Y(\boldsymbol{s})] or \mathrm{Cov}_{\hat{\theta}}[\widehat{Z}(\boldsymbol{s}),Z(\boldsymbol{s})],

• var.target: the variances of the prediction targets \mathrm{Var}_{\hat{\theta}}[Y(\boldsymbol{s})] or \mathrm{Var}_{\hat{\theta}}[Z(\boldsymbol{s})].

Note that the component var.pred is also present if type is equal to "trend", irrespective of the choice for extended.output. This component contains then the variances of the fitted values.

Value

The method predict.georob returns, depending on its arguments, the following objects:

If type is equal to "terms" then a vector, a matrix, or a list with prediction results along with bounds and standard errors, see predict.lm. Otherwise, the structure and contents of the output generated by predict.georob are determined by the class of newdata and the logical flags full.covmat and extended.output:

If full.covmat is FALSE then the result is an object of a "similar" class as newdata (data frame, SpatialPointsDataFrame, SpatialPixelsDataFrame SpatialGridDataFrame,
SpatialPolygonsDataFrame).

The data frame or the data slot of the Spatial...DataFrame objects have the following components:

• the coordinates of the prediction points (only present if newdata is a data frame).

• pred: the Kriging predictions (or fitted values).

• se: the root mean squared prediction errors (Kriging standard errors).

• lower, upper: the limits of tolerance/confidence intervals,

• trend, var.pred, cov.pred.target, var.target: only present if extended.output is TRUE, see Details.

If full.covmat is TRUE then predict.georob returns a list with the following components:

• pred: a data frame or a Spatial...DataFrame object as described above for
  full.covmat = FALSE.

• mse.pred: the full covariance matrix of the prediction errors, Y(\boldsymbol{s})-\widehat{Y}(\boldsymbol{s}) or Z(\boldsymbol{s})-\widehat{Z}(\boldsymbol{s}) see Details.

• var.pred: the full covariance matrix of the Kriging predictions, see Details.

• cov.pred.target: the full covariance matrix of the predictions and the prediction targets, see Details.

• var.target: the full covariance matrix of the prediction targets, see Details.

The function control.predict.georob returns a list with control parameters to steer predict.georob, see arguments of the function above for its components.

Author(s)

Andreas Papritz [e-mail].

References

Nussbaum, M., Papritz, A., Baltensweiler, A. and Walthert, L. (2014) Estimating soil organic carbon stocks of Swiss forest soils by robust external-drift kriging. Geoscientific Model Development, 7, 1197–1210. doi:10.5194/gmd-7-1197-2014.

Künsch, H. R., Papritz, A., Schwierz, C. and Stahel, W. A. (2011) Robust estimation of the external drift and the variogram of spatial data. Proceedings of the ISI 58th World Statistics Congress of the International Statistical Institute. doi:10.3929/ethz-a-009900710

See Also

georobPackage for a description of the model and a brief summary of the algorithms;

georob for (robust) fitting of spatial linear models;

georobObject for a description of the class georob;

profilelogLik for computing profiles of Gaussian likelihoods;

plot.georob for display of RE(ML) variogram estimates;

control.georob for controlling the behaviour of georob;

georobModelBuilding for stepwise building models of class georob;

cv.georob for assessing the goodness of a fit by georob;

georobMethods for further methods for the class georob;

lgnpp for unbiased back-transformation of Kriging prediction of log-transformed data;

georobSimulation for simulating realizations of a Gaussian process from model fitted by georob; and finally

sample.variogram and fit.variogram.model for robust estimation and modelling of sample variograms.

Examples

data(meuse)

data(meuse.grid)
coordinates(meuse.grid) <- ~x+y
meuse.grid.pixdf <- meuse.grid
gridded(meuse.grid.pixdf) <- TRUE

data(meuse.blocks, package = "constrainedKriging")

r.logzn.rob <- georob(log(zinc) ~ sqrt(dist), data = meuse, locations = ~ x + y,
    variogram.model = "RMexp", param = c(variance = 0.15, nugget = 0.05, scale = 200),
    tuning.psi = 1., control = control.georob(cov.bhat = TRUE, full.cov.bhat = TRUE))

## point predictions of log(Zn)
r.pred.points.1 <- predict(r.logzn.rob, newdata = meuse.grid.pixdf,
    control = control.predict.georob(extended.output = TRUE))
str(r.pred.points.1, max = 3)

## back-transformation of point predictions
r.backtf.pred.points <- lgnpp(r.pred.points.1)
str(r.backtf.pred.points, max = 3)

spplot(r.backtf.pred.points, zcol = "lgn.pred", main = "Zn content")

## predicting mean Zn content for whole area
if(interactive()){
  ## example is run only in interactive session because cpu times exceeds 5 s
  ## recompute point predictions with argument full.covmat = TRUE
  r.pred.points.2 <- predict(r.logzn.rob, newdata = meuse.grid.pixdf,
      control = control.predict.georob(extended.output = TRUE, full.covmat = TRUE))
  str(r.pred.points.2, max = 3)
  r.block <- lgnpp(r.pred.points.2, is.block = TRUE, all.pred = r.backtf.pred.points@data)
  r.block
}

## block predictions of log(Zn)
if(interactive()){
  ## example is run only in interactive session because cpu times exceeds 5 s
  r.pred.block <- predict(r.logzn.rob, newdata = meuse.blocks,
      control = control.predict.georob(extended.output = TRUE,
          pwidth = 75, pheight = 75, mmax = 50))
  r.backtf.pred.block <- lgnpp(r.pred.block, newdata = meuse.grid)

  spplot(r.backtf.pred.block, zcol = "lgn.pred", main = "block means Zn content")
}