LatticeKrig: User-friendly spatial prediction and inference using a...
In LatticeKrig: Multi-Resolution Kriging Based on Markov Random Fields

Description Usage Arguments Details Value Author(s) See Also Examples

This is a simple and high level function to fit the LatticeKrig spatial model to a data set. In is simplest form for 2-d spatial data: obj<-LatticeKrig(x,y) will fit a LatticeKrig type model to 2d locations x and observation vector y. Several (actually many!) default choices are made for the multi-resolution spatial covariance in this top level function. It uses the defaults that are based on a "thin-plate spline" like model for the spatial estimator and also uses LKrigFindLambda to estimate some covariance parameters (sill and nugget variances) through likelihood maximization (i.e. estimates the measurement and process variances.) For the simplest "black box" use only the observations and their 2-d locations need to be specified. But please see the caveats below in the Details section.

Despite the simple syntax, LatticeKrig still takes advantage of the multi-resolution features of the basic LKrig function and any LKrig parameter can be passed through the function call. See the example below for varying the range parameter. Also, see LKinfo andLKrigSetup for documentation on the complete object that describes the LatticeKrig model and the function to create it easily. See LKGeometry for documentation on extending or adding other spatial models to this package.

The returned value from this function can be used subsequently for prediction, conditional simulation, and other parts of the spatial analysis. See predict.LKrig and LKrig.sim.conditional

LatticeKrig(x, y, Z = NULL, nlevel = 3, findAwght = FALSE, LKinfo = NULL,
 X=NULL, U=NULL, na.rm =
                 TRUE, tol = 0.005, verbose = FALSE, ...)
## S3 method for class 'LatticeKrig'
print( x, digits=4, ...)

`x`	Spatial locations of observations. For the `LatticeKrig` function this should be a matrix where the columns index the spatial dimensions and rows index the observations. For example for 100 2-d locations, `x` would be a 100X2 matrix. Or for the function `print.LatticeKrig` `x` is the returned object from the `LatticeKrig` function.
`y`	Spatial observations. No missing values are allowed.
`Z`	Linear covariates to be included in fixed part of the model that are distinct from the default first order polynomial in `x` (i.e. the spatial drift).
`X`	For linear inverse problems the matrix that maps coefficients of the basis to the predicted values of observations. X must be in spam format. To convert from spind or dense format to spam format see `help(spam)` as an alternative `spind2spam`. See an example for this extension in the `LKrig` help file.
`U`	For linear inverse problems the matrix that maps coefficients of the fixed part of the model to the predicted values of observations. This needs to specified along with `X`
`nlevel`	Number of levels for the multi-resolution basis. Each level increases the number of basis functions by roughly a factor of 4.
`findAwght`	If FALSE the default a.wght parameter (related to correlation range) is set to mimic a thin plate spline. If TRUE this parameter and hence the range is estimated my maximum likelihood.
`LKinfo`	An optional list giving the full specification of the covariance. If this is missing it will be created internally and returned. If passed this parameterization will be used except lambda will be re-estimated by maximum likelihood.
`na.rm`	If TRUE NA's are removed from `y` and `x` is subsetted.
`tol`	Tolerance for the log likelihood used to judge convergence.
`verbose`	If TRUE print out intermediate results.
`...`	Additional arguments to pass to LKrig. The easiest way to pass a full specification is to create an LKinfo object beforehand and then just pass that (see example below.) This gives more control and the setup function will do some error checking on arguments. Also see help(LKrig) for a complete list of arguments to pass. For convenience we note that if you get some pesky memory warnings from spam you can set the storage higher by adding the argument `choleskyMemory`. For example to bump up to 2E6 include: `choleskyMemory=list(nnzR= 2E6)`.
`digits`	Number of significant digits in printed summary.

Keep in mind that overall LatticeKrig is just a specific type of spatial estimate that is designed to handle larger size data sets. It focuses on a specific form of covariance function, but the estimator is still the Kriging/Multivariate Gaussian Conditional Expectation/BLUE that is standard in this field.

The simplest model fit is:

Y_k = p(x_k) + h(x_k) + e_k

Y_k is the k^{th} observation at location x_k with measurement error e_k. Here p(x) is a low order polynomial of degree m-1 with the default m==2. h(x) is a mean zero Gaussian process with the representation:

h(x)= ∑_{l=1}^L ∑_{j=1}^{m(j)} φ_{m,l}(x) c_{m,l}

where φ are multi-resolution basis functions and the coefficients have mean zero and spatial dependence specified by a Markov random field. Keep in mind that this unusual form still implies a specific covariance function for h(x). In fact one can use the Krig or mKrig function from fields to reproduce the LatticeKrig estimate for smaller data sets and check computations. (See LKrig.cov with examples ). Details on the basis functions and the Markov random field are given in the LKrig help function. Throughout this package we assume the standard deviation of e_k is sigma and the marginal variance of h(x) is rho. An important derived parameter of the spatial estimate is lambda = sigma^2/ rho the noise to signal ratio. sigma and rho are estimated my restricted maximum likelihood in LatticeKrig.

This top level function is built on the more basic function LKrig supports a very flexible covariance. LKrig depends on the parameters nlevel, a.wght and alpha specifying all these relevant parameters may be discouraging to a new (or impatient!) user. Thus LatticeKrig is a "wrapper" that generates some simplified, default model choices to call the more general function LKrig. It is useful for users not fully familiar with the LatticeKrig methodology or those that wish to try a default approach to get a quick look at a spatial analysis. You always go back and add some specific non default choices to the LatticeKrig call (e.g. changing a.wght). For the 2-dimensional case the default values are set to give about 4 times as many basis functions as observations, use 5 extra lattice points on the edge to minimize boundary effects, and to use four levels of multi-resolution. An important default is that a linear spatial drift is included in the model so the model will relax to a linear prediction based on the x values in the absence of a spatial component. In other words, the model includes by default a fixed part that is linear in x. The spatial correlation range is nearly stationary and set large to mimic a thin-plate spline. The smoothness mimics the Whittle covariance function ( smoothness = 1 for the Matern). (See LKrig.cov.plot to get a plot of the implied covariance function.) LatticeKrig also provides maximum likelihood estimates of the measurement error standard deviation ("sigma") and process variance parameter ("rho") that are perhaps the parameters that most effect the shape of the estimated spatial field. The ubiquitous parameter lambda throughout LatticeKrig is just the reparameterization lambda == sigma^2 / rho.

This top level function is pitched with all the caveats that statistical model assumptions should always be checked and applying generic methods to a specific problems without checking the appropriateness can lead to misleading results. So plot your data and try several models. Details on the full computations can be found in the LKrig man page. The lambda = sigma2/ rho parameter in LKrig is essential to the Lattice Krig computation and an inappropriate value will result in over or under fitting and incorrect interpolated values. The function LKrigFindLambda is used within LatticeKrig to estimate a lambda value from the data using maximum likelihood.

One interesting feature of this package is the ability to handle spatial processes on different geometries and the form is specified by the LKGeometry argument. The current choices are:

LKRectangle: A 2 dimensional Euclidean spatial domain. The default
LKInterval: A 1 dimensional Euclidean spatial domain.
LKBox: A 3 dimensional Euclidean spatial domain.
LKRing: A 2 dimensional spatial domain where the first coordinate is periodic (on [0,360]) and the second is Euclidean. E.g. a slice around the equator and not having a large latitudinal range.
LKCylinder: A 3 dimension model where an additional coordinate is added to the LKRing geometry. This is useful for representing a small section of the sphere where one also has a height component.
LKSphere: A full 2-d spherical geometry. Coordinates are given in longitude, latitude but the distances and any structures are on the sphere.

One important feature of this package is that the different geometries all use the same computation engine LKrig, following the same computational algorithm. The differences in setting up the problem and in evaluating the function are implemented as S3 methods. The details of this strategy are described in LKGeometry and allow the user to add new geometries.

This function also supports a model where the observations are simply expressed as linear combinations of the basis function coefficients. Equivalently this is a model where the observed data can be expressed as linear functionals applied to the polynomial term and the spatial process. Typically these linear maps represent observing integrals or weighted combinations of the fields and are important for data that is aggregated over by space. See help(LKrig) for an example of how this model is set up at the end of the Examples section.

The main call inside LatticeKrig is to LKrig and so a LKrig object is returned. Thus all of the functionality that comes with LKrig objects such as predict, summary, predictSurface, etc. remain the same as described in LKrig. Also, see the components residuals and fitted.values in the returned object for these parts of the fitted spatial model. The component LKinfo has all the details that specify the basis functions and co variance model. The component MLE gives details of the likelihood evaluations to estimate the sigma and rho parameters.

Doug Nychka

LKrig, LKrig.setup, LKrigFindLambda, LKinfo, LKrig.sim.conditional

# Load ozone data set
  data(ozone2)  
  x<-ozone2$lon.lat
  y<- ozone2$y[16,]

# thin plate spline-like model with the lambda parameter estimated by
# maximum likelihood. Default choices are made for a.wght, nlevel, NC
# and alpha.

  obj<- LatticeKrig( x, y)
## Not run: 
# summary of fit and a plot of fitted surface
  print( obj)
  surface( obj )
  US(add=TRUE)
  points(x)
# prediction standard errors
  out.se<- predictSE( obj, xnew= x)
# predict at observations:
  out.fhat<- predict( obj, xnew= x)
# conveniently predict on a 100X100 grid for plotting
 out.surf<- predictSurface( obj, nx=100, ny=100)
# image.plot( out.surf) 

## End(Not run)
# running an example by first setting up the model object
## Not run: 
# this is just a small model to run quickly
# compare the LKinfo object here  to one created implicitly:  obj$LKinfo
LKinfo1<- LKrigSetup( x, NC=5, nlevel=3, a.wght=4.1, nu=1.0)
obj1<- LatticeKrig( x,y, LKinfo= LKinfo1)

## End(Not run)
#
# In this example lon/lat are treated as just Euclidean coordinates 
# a quick adjustment for small regions is to account for the difference
# in physical distance in N-S verses E_W
# is to just scale the longitude degrees to be comparable to degrees in latitude
# at least in the middle of the domain. The assumption is that for small spatial
# domains this approximation will not be bad for the coordinates at the edges too.
# You accomplish this by adding a scaling, V matrix:
# Here the V argument is rolled into the LKinfo object created within the function
#
## Not run: 
  meanLat<- mean( x[,2])*pi/180
  Vlonlat <- diag(  c( 1/cos(meanLat), 1) )
  obj1<- LatticeKrig( x, y, V = Vlonlat )

## End(Not run)

## Not run: 
# Refit using with just one level of  basis functions
# on a 20X20 grid within the spatial domain ( so about 400) 
# actually number is 720 ( see obj1b$LKinfo) due adding edge nodes
# Add an aspect ratio of spatial domain 
# and find the a.wght parameter along with nugget and process variances.
# this takes a while partly because LatticeKrig model is not optimized for small data sets!
  obj1b<- LatticeKrig( x, y, nlevel=1, NC=20, findAwght=TRUE)
# rudimentary look at how likelihood was optimized
#log lambda and omega =  log(a.wght-4)/2 are useful parameterization ...
  quilt.plot( obj1b$MLE$lnLike.eval[,c("logLambda","omega")],
       obj1b$MLE$lnLike.eval[,"lnProfileLike.FULL"], 
       xlab="loglamda", ylab="omega",
       zlim =c(-640,-612))
  points( obj1b$MLE$lnLike.eval[,c("logLambda","omega")],cex=.25)
      

## End(Not run)
# fitting replicate spatial data sets
# here we use the common observations over days for the ozone
# data set. Whether these are true replicated fields is in question
# but the analysis is still useful

## Not run: 
Y<-  na.omit( t( ozone2$y) ) 
ind<- attr( Y,"na.action")
X<- ozone2$lon.lat[-ind, ]

out1<- LatticeKrig( X, Y, nlevel=1, NC=20, findAwght=TRUE)
out2<- LatticeKrig( X, Y, nlevel=1, NC=20, findAwght=TRUE,
                        collapseFixedEffect=TRUE)
# compare the two models 
# Note second a.wght reflects more spatial correlation when individual 
# fixed effect is not removed ( 4.4 verses 4.07)
# nugget variance is nearly the same!
out1$MLE$summary[1:7]                        
out2$MLE$summary[1:7]
                        



## End(Not run)
## Not run: 
# Refit using the tensor product type of basis functions
# (default is "Radial"). An example how an additional argument that is 
# passed to the LKrigSetup function to create the LKinfo object.
  obj2<- LatticeKrig( x, y, BasisType="Tensor")

## End(Not run)

#
# A 1-d example with 3 levels of basis functions
# See LKrig for an explanation if nlevel, NC,  alpha and a.wght 
# covariance parameters.


## Not run: 
 x<- matrix(rat.diet$t)
 y<- rat.diet$trt
 fitObj<- LatticeKrig( x, y)
# NOTE lots of defaults are set for the model! See print( fitObj)
 plot( x,y)
 xg<- matrix(seq( 0,105,,100))
 lines( xg, predict(fitObj, xg) )

## End(Not run)

## Not run: 
#  a 3D example
set.seed( 123)
N<- 1000
x<-  matrix( runif(3* N,-1,1), ncol=3, nrow=N)
y<-   10*exp( -rdist( x, rbind( c(.5,.5,.6) ) )/.5)

# NOTE setting of memory size for Cholesky. This avoids some warnings and
# extra computation by the spam package
LKinfo<- LKrigSetup( x,  nlevel=1,  a.wght= 6.01, NC=6, NC.buffer=2,
                    LKGeometry="LKBox", normalize=FALSE, mean.neighbor=200,
                    choleskyMemory=list(nnzR= 2E6) )                                      
out1<- LatticeKrig( x,y, LKinfo=LKinfo)

glist<- list( x1=seq( -1,1,,30), x2=seq( -1,1,,30), x3 = 0)
xgrid<- make.surface.grid( glist)

yhat<- predict( out1, xgrid)
# compare yhat to true function created above
image.plot( as.surface( glist, yhat))


## End(Not run)
#
###########################################################################
# Including a covariate (linear fixed part in spatial model)
########################################################################## 
## Not run: 
  data(COmonthlyMet)

  obj  <- LatticeKrig(CO.loc,  CO.tmin.MAM.climate, Z=CO.elev)
  obj2 <- LatticeKrig(CO.loc, CO.tmin.MAM.climate)

# compare with and without linear covariates
  set.panel(1,2)
  surface(obj)
  US(add=TRUE)
  title("With Elevation Covariate")

  surface(obj2)
  US(add=TRUE)
  title("Without Elevation Covariate")


## End(Not run)
## Not run: 
 data(COmonthlyMet)
# Examining a few different "range" parameters
a.wghtGrid<-  4  +  c(.05, .1, .5, 1, 2, 4)^2

#NOTE smallest is "spline like" the largest is essentially independent
# coefficients at each level.  In this case the "independent" end is
# favored but the eff df. of the surface is very similar across models
# indicating about the same separate of the estimates into spatial
# signal and noise
#
for( k in 1:5 ){
obj  <- LatticeKrig(CO.loc,  CO.tmin.MAM.climate, Z=CO.elev, 
                      a.wght=a.wghtGrid[k])
cat( "a.wght:", a.wghtGrid[k], "ln Profile Like:",
            obj$lnProfileLike, "Eff df:", obj$trA.est, fill=TRUE)
}

# MLE
obj0  <- LatticeKrig(CO.loc,  CO.tmin.MAM.climate, Z=CO.elev, 
                     findAwght=TRUE)
print(obj0$MLE$summary)

## End(Not run)

#########################################################################
# Reproducing some of the analysis for the example in the
# JCGS LatticeKrig paper.
#########################################################################

#### Here is an example of dealing with approximate spherical geometry.
## Not run: 
data(NorthAmericanRainfall)
library(mapproj)
x<- cbind(NorthAmericanRainfall$longitude, NorthAmericanRainfall$latitude)
y<- NorthAmericanRainfall$precip
log.y<- log(y)
elev<- NorthAmericanRainfall$elevation
# this is a simple projection as part of this and handled by the mapproj package
x.s<- mapproject( x[,1], x[,2], projection="stereographic")
x.s<- cbind( x.s$x, x.s$y)

# an alternative is to transform coordinates using another projection,
# e.g. a Lambert conformal projection
# with the project function from the rgdal package
# library( rgdal)
# x.s<- project(x,"+proj=lcc +lat_1=22 +lat_2=58 +lon_0=-93 +ellps=WGS84")
# this package has the advantage that the inverse projection is also 
# included ( inv=TRUE) so it is easy to evaluate the surface back on a Mercator grid.
             
obj0<- LatticeKrig(x.s, log.y, Z=elev )

fitSurface<- predictSurface( obj0, drop.Z=TRUE)
fitSurface$z<-  exp(fitSurface$z)/100
colorTable<- designer.colors( 256, c("red4", "orange", "yellow","green1", "green4", "blue"))
image.plot( fitSurface, col=colorTable)
map( "world", add=TRUE, col="grey30", lwd=3, proj="") 


## End(Not run)