predictSE: Standard errors of predictions for Krig spatial process estimate

Description

Finds the standard error ( or covariance) of prediction based on a linear combination of the observed data. The linear combination is usually the "Best Linear Unbiased Estimate" (BLUE) found from the Kriging equations. This statistical computation is done under the assumption that the covariance function is known.

Usage

predictSE(object, ...) 
# S3 method for Krig
predictSE(object, x = NULL, cov = FALSE, verbose = FALSE,...)
# S3 method for mKrig
predictSE(object, xnew = NULL, XMat = NULL, verbose = FALSE, drop.XMat
                 = FALSE, ...)

Value

A vector of standard errors for the predicted values of the Kriging fit.

Arguments

drop.XMat: If FALSE find standard error without including the additional spatial covariates described by XMat. If TRUE find full standard error with spatial covariates if they are part of the model.
object: A fitted object that can be used to find prediction standard errors. This is usually from fitting a spatial model to data. e.g. a Krig or mKrig object.
xnew: Points to compute the predict standard error or the prediction cross covariance matrix.
x: Same as xnew -- points to compute the predict standard error or the prediction cross covariance matrix.
cov: If TRUE the full covariance matrix for the predicted values is returned. Make sure this will not be big if this option is used. ( e.g. 50X50 grid will return a matrix that is 2500X2500!) If FALSE just the marginal standard deviations of the predicted values are returned. Default is FALSE -- of course.
verbose: If TRUE will print out various information for debugging.
...: These additional arguments passed to the predictSE function.
XMat: Additional matrix of spatial covariates used for prediction. These are used to determine the additional covariance contributed in teh fixed part of the model.

Details

The predictions are represented as a linear combination of the dependent variable, Y. Call this LY. Based on this representation the conditional variance is the same as the expected value of (P(x) + XMat(X) - LY)**2. where P(x)+XMat(x) is the value of the surface at x and LY is the linear combination that estimates this point. Finding this expected value is straight forward given the unbiasedness of LY for P(x) and the covariance for XMat and Y.

In these calculations it is assumed that the covariance parameters are fixed. This is an approximation since in most cases they have been estimated from the data. It should also be noted that if one assumes a Gaussian field and known parameters in the covariance, the usual Kriging estimate is the conditional mean of the field given the data. This function finds the conditional standard deviations (or full covariance matrix) of the fields given the data.

There are two useful extensions supported by this function. Adding the variance to the estimate of the spatial mean if this is a correlation model. (See help file for Krig) and calculating the variances under covariance misspecification. The function predictSE.KrigA uses the smoother matrix ( A(lambda) ) to find the standard errors or covariances directly from the linear combination of the spatial predictor. Currently this is also the calculation in predictSE.Krig although a shortcut is used predictSE.mKrig for mKrig objects.

Examples

Run this code

# 
 data(ozone2)
 sOzone<- ozone2$lon.lat
 yOzone<- ozone2$y[16,]
# Matern smoothness =1.0 sigma2 and tau found by ML
  fit0<- spatialProcess( sOzone,yOzone, aRange=10)
# prediction SEs at  miss data locations
  missingObs<- is.na(yOzone)
  missingSE<- predictSE(fit0, sOzone[missingObs,]) 
# NOTE based on dispatching predictSE.mKrig is called. 
# 95% CI for  ozone surface
   fHat<- predict(fit0, sOzone[missingObs,]) 
   cbind(fHat - 1.96*missingSE, fHat + 1.96*missingSE )
# 95% CI for  ozone _observation_ 
   fHat<- predict(fit0, sOzone[missingObs,])
   tau2<- fit0$summary["tau"]^2
   cbind(fHat - 1.96*sqrt(missingSE + tau2) ,
         fHat + 1.96*sqrt(missingSE + tau2) )

# SE on  a  grid of locations  
  xg<-make.surface.grid( 
  list(x=seq(-91,-83,,30), 
       y=seq( 37, 45,,30))
  )     
  gridSE<- predictSE(fit0,xg)   # std errors of predictions 
  
#at the grid points out is a vector of length 900 
#reshape the grid points into a 30X30 matrix etc.  

   out.p<-as.surface( xg, gridSE) 
   surface( out.p, type="C") 
   points( sOzone[!missingObs,], col="magenta", pch=16)

# this is equivalent to the higher level functions
# (We are lazy and the grid will be a bit different here. See the 
# help file!)
 out<- predictSurfaceSE( fit0, nx=30, ny=30, extra=TRUE) 
 imagePlot( out, col=viridis(256))
 contour( out, add=TRUE, levels=c( 10,20, 20))
 points( sOzone[!missingObs,], col="magenta", pch=16)