fields (version 1.7.2)

krig.image: Spatial process estimate for large irregular 2-d dats sets.

Description

Computes the spatial predictions for large numbers of irregularly spaced observations using the standard Kriging equations. The main approximation is that the locations are discretized to a regular grid, but the field need not be observed at all grid boxes.

In Bayesian terms this function computes the posterior mean for the field given the observations under the usual Gaussian assumptions for the fields and observations. The solution is found by the iterative solution of a large linear system using the conjugate gradient algorithm (CGA). Part of the calculations rely on discretizing the spatial locations to a regular grid to make use of the FFT for fast multiplication of a covariance matrix with a vector.

Usage

krig.image(x, Y, cov.function, m=NULL, n=NULL, lambda=0, start=NULL, 
tol=1e-05, kmax=25,  cov.obj=NULL, grid=NULL, 
weights=rep(1, length(Y)), verbose=FALSE, conv.verbose=FALSE, expand=1, ...)

Value

An list object of class krig.image. An explanation of some components:callThe calling sequencecov.functionA copy of the covariance S functionna.indlogical indicating missing values in YxrawPassed spatial locations having removed missing valuesyObservations having omitted missing valuesNLength of yweightspassed weights having omitted missing cases.lambdagridlist with components x an y indicating grid for discretizationcov.objList object to use with cov.functionmNumber of grid point in x axisnNumber of grid point in y axisindexA two column matrix indicating the indices of the closest grid point to each observed location.xObserved locations discretized to nearest grid pointyMObserved values but with a weighted average replacing multiple values associated with the same grid point.xMDiscretized locations associated with yMweightsMWeight vector associated with YM.uniquerowsLogical indicating which rows of x are unique.shat.repPooled standard deviation among observations that fall within the same gird pointsindexMA two column matrix indicating the indices of the closest grid point to each observed location, yM.qr.TQR decomposition of the matrix of constant and linear terms at xMmultAxThe S function that is used for matrix multiplication in the CGA.omega2Parameter vector that describes the spatial process part of the conditional mean.convergeCGA convergence informationbetaConstant, and the two linear parameters for the fixed linear part of the modeldeltaCovariance matrix times delta give the spatial predictions.rhohat, rhoConditional on lambda the MLE for the parameter multiplying the covariance function.sigma2, shat.MLEConditional on lambda the MLE for the parameter dividing the weight function.surfaceA list giving the predicted surface at the grid points.fitted.valuesPredicted values at true locations

References

Large spatial prediction problems and nonstationary fields (1998) Nychka, D., Wikle, C. and Royle, J.A.

FIELDS manual

Details

From a functional point of view krig.image and supporting functions are similar to the class Krig. The main difference is that only 2-dimensional problems are considered and the solution is calculated for a fixed value of lambda. (The Krig function can estimate lambda.) For large data sets a practical way to estimate lambda is by out of sample cross-validation and the FIELDS manual gives a detailed example of this for the precip data set. Also see the manual for an explanation of the computational strategy (Conjugate Gradient) here.

Efficiency for large datasets comes with restrictions on the range of covariance functions and some other features. Currently FIELDS just has two covarince models: exponential/Gaussian and wavelet based. However, it is not difficult to modify these to other models. The default discretization is to a 64X64 grid however even 256X256 is manageable and quite likely to separate irregular locations in most cases. The user should also keep in mind that the estimate is the result of an iterative algorithm and so issues such as good starting values and whether the algorithm converged are present.

The spatial model includes a linear spatial drift and MLE estimates of the nugget variance and sill are found based on the values of lambda. If the weights are all equal to one and the covariance function is actually a correlation function, in the notation of this function, the "sill" is sigma2 + rho and the "nugget" is sigma2. Moreover sigma2 and rho are constrained so sigma2/rho =lambda. This is why lambda is the crucial parameter in this model.

Although the field is only estimated to the resolution of the grid, prediction off of the grid is supported by bilinear interpolation using the FIELDS function interp.surface.

See Also

plot.krig.image, predict.krig.image, exp.image.cov, sim.krig.image

Examples

# 
# fit a monthly precipitation  field over the Rocky Mountains 
# grid is 64X64 
out<- krig.image( x= RMprecip$x, Y = RMprecip$y, m=64,n=64,cov.function= 
exp.image.cov,  
lambda=.5, theta=1, kmax=100)  

#  
# range parameter for exponential here is .5 degree in lon and lat.  
#diagnostic plots.  
plot( out)  

# look at the surface 
image.plot( out$surface)  #or just      surface( out)

#
#simulate 4 realizations from the conditional distribution 
look<- sim.krig.image( out, nreps=4) 
# take a look:  plot( look)
 
# check out another values of lambda reusing some of the objects from the 
# first fit 

out2<- krig.image( RMprecip$x, RMprecip$y, cov.function= exp.image.cov, 
lambda=4, 
start= out$omega2,cov.obj=out$cov.obj)
#
# some of the obsare lumped together into a singel grid box
#
# find residuals among grid box means and predictions
res<- predict( out2, out2$xM) - out2$yM
#compare with sizes of out2$residuals (raw y data) 

#starting values from first fit  in out$omega2 
# covariance and grid information are 
# bundled in the cov.obj 
## 

#
## fitting a thin plate spline. The default here is a linear null space 
## and  second derivative type penalty term.  
##  you will just have to try different values of lambda vary them on  
##  log scale to 
 
out<- krig.image( RMprecip$x, RMprecip$y, cov.function=rad.image.cov, 
lambda=1, m=64, n=64, p=2, kmax=300) 

# take a look: image.plot( out$surface) 

# check out  different values reuse some of the things to make it quicker 
# note addition of kmax argument to increase teh number of iterations

out2<- krig.image( RMprecip$x, RMprecip$y,cov.function=rad.image.cov,  
lambda=.5, start= out$omega2, cov.obj=out$cov.obj, kmax=400) 

# here is something rougher  
out3<- krig.image( RMprecip$x, RMprecip$y,cov.function=rad.image.cov, 
lambda=1e-2, start= out2$omega2, cov.obj=out$cov.obj,kmax=400,
tol=1e-3) 

# here is something close to an interpolation  
out4<- krig.image( RMprecip$x, RMprecip$y,cov.function=rad.image.cov, 
lambda=1e-7, start= out3$omega2, cov.obj=out$cov.obj,kmax=500, tol=1e-3)
 
#compare the the four surfaces: 
# but note the differences in scales ( fix zlim to make them the same) 
# 
# take a look
# set.panel( 2,2) 
# image.plot( out$surface) 
# points( out$x, pch=".") 

# image.plot( out2$surface) 
# image.plot( out3$surface) 
# image.plot( out4$surface) 
 
# some diagnostic plots) 
set.panel( 4,4) 
plot( out, graphics.reset=FALSE) 
plot( out2, graphics.reset=FALSE) 
plot( out3, graphics.reset=FALSE) 
plot( out4, graphics.reset=FALSE) 
set.panel(1,1)