Package 'LatticeKrig'

Description

Convert back and forth between lon/lat and direction cosines and also project spherical coordinates on a local tangent plane.

Usage

directionCosines(x)
toSphere( Grid)
projectionSphere(x1, x2)

Arguments

Details

The conversion functions are based on straight forward definitions of spherical coordinates. The projection function is done by two rotations, first around the Z axis and then around the Y axis.

Value

directionCosines A three column matrix of direction cosines.

toSphere A two column matrix of longitudes and latitudes.

projectionSphere A two column matrix of Euclidean coordinates on the tangent plane to x1. The convention is that the origin (0,0) is mapped to x1 and the X- axis are points along the meridian through x1. The Y axis are points on the great circle passing through x1 and perpendicular to the meridian.

Author(s)

Doug Nychka

Examples

# 
# icosahedron:
 phi = (1 + sqrt(5))/2
  V = matrix(c(0, 1, phi, 0, 1, -phi, 0, -1, phi, 0, -1, -phi, 
        1, phi, 0, -1, phi, 0, 1, -phi, 0, -1, -phi, 0, phi, 
        0, 1, -phi, 0, 1, phi, 0, -1, -phi, 0, -1), 12, 3, byrow = TRUE)

# check : library( rgl); plot3d( V, size=10, col="red4" )
# as lon lat:
V2<- toSphere( V)
plot( V2)

# lon lat grid
 lGrid<- make.surface.grid(  list(x= seq( -10,10,, 10), y= seq( -20,20,,10)) )
 
 dGrid<- directionCosines( lGrid)
 pairs( dGrid)
 # also try:   library( rgl); plot3d(  dGrid)

Description

This object is mainly designed to work with methods that take a set of locations organized on a grid. The object is a list where there are as many components as dimensions and each list component is a vector of values being the grid points in that dimension. It is consistent with the older use of the older grid.list format used in the fields package. This form is somewhat redundant because for an equally spaced grid all one needs is the beginning value, spacing and number of points but it makes it simpler to pass the grid information to functions such as image and contour.

Usage

gridListInfo(gridList)

Arguments

Objects from the Class

This object is a list where each component is a vector of grid points in a particular dimension. For example

grid<- structure(list( x= seq( -1,1,,20), y= seq( 0,1,,15)), class= "gridList"

would create this object for a 2d grid with 20 and 15 points over the ranges [-1,1] and [0,1]. The component names ( "x" and "y" in this case) are optional.

The function gridListInfo extracts some summary information that is used to support the summary function for this class.

Methods

LKrigDistance

signature(x1 = "matrix", x2 = "gridList", delta = "numeric"): ...

Author(s)

Doug Nychka

See Also

LKrigDistance and LKrigDistGrid LKrigLatticeCenters

Examples

showClass("gridList")
# a 3-d grid
grid<- structure(
 list( x= seq( -1,1,,20), y= seq( 0,1,,15) ,oneMore = 1:10) ,
 class= "gridList" )

Description

Creates a multi-resolution grid based on subdividing triangles starting with faces of an icosahedron.

Usage

IcosahedronGrid(K)

IcosahedronFaces(K)

Arguments

Details

Creates a nearly regular grid by taking the first level as the 12 points from a regular icosahedron. The subsequent levels generate a finer set of points by subdividing each triangular face into 4 new triangles. The three new mid points from the subdivision are added to the previous nodes to give the new level of resolution. The triangles tend to be roughly equilateral and so the nodes will tend to be roughly equally space but there is some variation in distances among nearest neighbors.

To depict the faces and nodes in a snazzy way use the rgl package and the following code

library( rgl)
# show level 3
Level<- 3
    SGrid <- IcosahedronFaces(4)
    Tri <- SGrid$Face[[Level]]
    L <- dim(Tri)[3]
    plot3d(rbind(c(0, 0, 0)), xlim = c(-1, 1), ylim = c(-1, 1), 
        zlim = c(-1, 1), axes = FALSE, xlab = "", ylab = "", 
        zlab = "", type = "p", box = TRUE)
    for (k in 1:L) {
        U <- Tri[, , k]
        rgl.triangles(U, col = "grey80")
    }
    plot3d(SGrid$nodes[[Level]], col = "green4", type = "s", 
        radius = 0.03, add = TRUE)

Value

IcosahedronGrid A list with K components each component is a three column matrix giving the direction cosines for each grid point.

IcosahedronFaces A list with components:

MultiGrid

The same list returned by IcosahedronGrid.

Faces

A list with K-1 components. Each components are the faces at a given level represented as a three dimensional array (3X3XN with N the number of faces and a given level). The array indices are vertices of triangle, coordinates and faces within a resolution level. e.g. to extract the 10th face (out of 80) for the 2nd level:

look<- IcosahedronFaces(3)$Faces
triangle <- (look[[2]])[,,10]
print (triangle)
           [,1]       [,2]      [,3]
[1,] -0.5257311 -0.8506508  0.000000
[2,] -0.8090170 -0.5000000  0.309017
[3,] -0.8090170 -0.5000000 -0.309017
rowSums( triangle^2)
[1] 1 1 1

triangle will be the 10th face for the second level where the columns are the 3d coordinates of the direction cosines and the rows index the three vertices.

Author(s)

Doug Nychka and Zachary Thomas

See Also

toSphere, directionCosines, projectionSphere

Examples

# second level in lon lat coordinates 
look<- IcosahedronGrid(3)
lonlat<- toSphere( look[[3]])
plot( lonlat, xlab="lon", ylab="lat")

Description

This is a sample of 50 (x,y) ordered pairs: the x values are in sorted order, drawn from the interval [0,1]. The smallest x value is exactly 0, and the largest is exactly 1. The y values are the function y=9*x*(1-x)^3 evaluated at each x value, with a small measurement error (normal with mean 0, standard deviation 0.01) added on.

Usage

data(KrigingExampleData)

Format

The data in KrigingExampleData is listed in two vectors, x and y.

Details

The following code was used to generate the data:

set.seed(123)
x <- c(0,sort(runif(48)),1 )  
y <-  9*x*(1-x)^3  + rnorm(50, sd=0.01)
KrigingExampleData<- list( x=x, y=y)
save( KrigingExampleData, file="KrigingExampleData.rda")

Examples

## Not run: 
data(KrigingExampleData)
x<- KrigingExampleData$x
y<- KrigingExampleData$y
plot(x,y)
kFit <- LatticeKrig(x, y)
xGrid <- seq(0,1,0.001)
lines(xGrid, predict(kFit, xGrid), col='red')

## End(Not run)

Description

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 and also see the vignette https://github.com/NCAR/LatticeKrig/tree/master/Vignette for a gentle introduction.

Despite the simple syntax, the LatticeKrig function 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

Usage

LatticeKrig(x, y, Z = NULL,weights = 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, ...)

Arguments

Details

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)= \sum_{l=1}^L \sum_{j=1}^{m(j)} \phi_{m,l}(x) c_{m,l}$

where $\phi$ 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 tau and the marginal variance of $h(x)$ is sigma2. An important derived parameter of the spatial estimate is lambda = tau^2/ sigma2 the noise to signal ratio. tau and sigma2 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 ("tau") and process variance parameter ("sigma2") 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 == tau^2 / sigma2.

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 = tau^2/ sigma2 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.

Value

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 tau and sigma2 parameters.

Author(s)

Doug Nychka

See Also

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

Examples

# 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:
# here x can be any two column matrix of coordinates this
# function returns a vector of predictions
  out.fhat<- predict( obj, xnew= x)
  
# conveniently predict on a 100X100 grid for plotting
# use the grid.list arugment to give more control over the grid choice.
# output object is the standard list with components x, y, and z
# suitable for contour, persp, image, etc.
 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
# this is the main way to specify  the spatial model components 
## 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)

Description

This function builds a matrix in which each entry on a given diagonal has the same value. The matrix can be any size, not necessarily square, and the user can specify arbitrary diagonals to fill. If ncol is not set, the output will be a square matrix. If diags is not set, the entries will be placed as close as possible to the main diagonal; e.g. passing in one value for entries will produce a diagonal matrix, and passing in a vector of 3 values will produce a tridiagonal matrix. If an even number of entries are provided, they will be placed symmetrically around the main diagonal, leaving the main diagonal as zeroes.

Usage

LKDiag(entries, nrow, diags = NULL, ncol = nrow,
full = FALSE)

Arguments

Value

LKDiag returns a (dense) matrix with the given diagonals.

Author(s)

Matt Iverson

Examples

#produces a 5x5 identity matrix
LKDiag(1, 5)

#produces a 6x6 tridiagonal matrix
LKDiag(c(1, -2, 1), 6)
LKDiag(c(1, -2, 1), 6, diags=(-1:1))

#produces a 5x5 matrix with 3 on the main diagonal and the corners
n <- 5
LKDiag(3, n, diags=c(0, n-1, -(n-1)))

# each of the following produces a 4x6 matrix with a 2 in the diagonal 
# below the main diagonal and 5 in the diagonal  above it
LKDiag(c(2, 5), nrow = 4, ncol = 6)
LKDiag(c(2, 5), nrow = 4, ncol = 6, diags = c(-1, 1))
LKDiag(c(5, 2), nrow = 4, ncol = 6, diags = c(1, -1))

Description

These are the lower level functions to compute the distances among two sets of locations but being limited to distances less than a maximum threshold (see delta below ). These functions are useful for generating a sparse matrix on distances and evaluating a compactly supported function (such as the Wendland). The location - location method supports the distance metrics: Euclidean, spherical, component-wise and Manhattan. LKDistComponent and LKDistComponentGrid return the coordinate-wise distances and are useful for evaluating a tensor product basis functions.

Usage

LKDist(x1, x2, delta, max.points = NULL, mean.neighbor = 50,
             distance.type = "Euclidean")
LKDistComponents(x1, x2, delta, max.points = NULL, mean.neighbor = 50, 
    distance.type = "Euclidean")
    
LKDistGrid(x1, gridList, delta, max.points = NULL, mean.neighbor = NULL,
distance.type = "Euclidean", periodic)

LKDistGridComponents(
x1, gridList, delta,
  max.points = NULL, mean.neighbor = NULL,
  distance.type = "Euclidean") 

LKGridCheck(distance.type, x1, gridList )

LKGridFindNmax(n1, max.points, mean.neighbor, delta, gridList)

Arguments

Value

LKDist and LKDistGrid a list representing a sparse matrix in spind format.

LKDistComponent andLKDistGridComponent a list representing a sparse matrix using the spind format except the ra component is now a matrix. The columns of ra being the individual distances for each coordinate.

directionCosine a matrix where rows index the different points and columns index x,y,z.

LKGridFindNmax returns the maximum number of nonzero elements expected in a pairwise distance matrix.

LKGridFindNmax checks that the calling arguments are compatible with the pairwise distance computation.

Author(s)

Doug Nychka

See Also

spind2spam, spind2full

Examples

set.seed( 123)	
x<- matrix( runif(100*2), 100,2)

DMatrix<- LKDist( x,x, delta=.1)
# coerce to spam matrix format
DMatrix2<- spind2spam( DMatrix)

# a grid
gridL<- list( x1= seq(0,1,.2), x2= seq( 0,2,.2) , x3= seq( -1,1,.2))
class(gridL)<- "gridList"	
x1<- cbind( runif( 100), runif(100)*2, 2*(runif( 100) -.5) )
look<- LKDistGrid( x1, gridL, delta=.45)
# check against rdist.
# look2<- rdist( x1, make.surface.grid(gridL))
# look2[ look2 >= .45] <- 0
# max( abs(look- look2)[look>0] )

# test of periodic option
 gridL<- structure(
            list( x1= seq(0,1,.02),
                  x2= seq( 0,1,.02)),
            class="gridList")
 look1<- LKDistGrid( rbind(c(0,0)), gridL, delta=.35,
                     periodic=c(TRUE,FALSE))
 look2<- spind2full(look1)
 image.plot( as.surface( gridL, look2) )
 
 look1<- LKDistGrid( rbind(c(0,0)), gridL, delta=.35,
                      periodic=c(TRUE,TRUE))
 look2<- spind2full(look1)
 image.plot( as.surface( gridL, look2) )

Description

To implement a spatial model for specific geometry, e.g. a rectangular spatial domain, one needs to specify a few functions that set reasonable default choices, create the lattice, define the Markov random field and define the basis functions. One could also include some specific functions that are particularly efficient for some of the computations for that geometry. This help section gives an overview of what is needed to introduce a new geometry.

A quick start

At the outset if you are an example person (and an impatient one) take a look at the functions created for the 1-d LatticeKrig model (LKInterval). (In the source for the package this is the file ModelLKInterval.R in the R subdirectory.) The specific functions are:

LKrigSetupLattice.LKInterval

Creates a 1-d lattice.

LKrigSAR.LKInterval

Creates a sparse 1-d spatial autoregressive matrix.

LKrigLatticeCenters.LKInterval

Returns the spatial locations of the lattice points.

Nothing else is needed to coax the LatticeKrig computation to fit a 1-d spatial model. There are more details and flexibility in specifying the model, but the ones listed above are generic in that they make sense for any LatticeKrig model and are not specific to the 1-d case. The .LKInterval tag on each of these functions indicates that these are S3 versions of a "method". The actual source code in LatticeKrig just relies on calling the generic method, e.g. LKrigSetupLattice, and the geometry, added as a class indicates, which function is actually used. Checking what each of these functions does and then seeing how they are called in LKrigSetup and LKrig.basis will give a simple introduction to how the code is structured. See LKRectangle for details on the most common 2-d geometry.

Use of the LKInfo object

The basic idea behind this package is to consolidate all the details of the spatial model including the geometry in the LKinfo object. The way to create this object is by the function SetupLKinfo and we believe that this function has enough flexibility to avoid modifying this object outside of this function call. Typically for deliberate spatial modeling the first step will be creating the LKinfo object; however, the quick and more "black box" function LatticeKrig is designed to take a minimal amount of information and then calls SetupLKinfo internally to fill out the LKinfo object. In fact, LatticeKrig( x,y) with just x as locations and y as observations will fit a spatial model using reasonable defaults and the returned object from this fit will include the full LKinfo object that it used.

A subtle and perhaps a potential disadvantage with the coding is that the LKinfo object is built in steps within the LKrigSetup function. In particular the method that setups the lattice is called before information is added about the alpha and a.wght parameters. The source code LKrigSetup has been kept concise and it should be clear how the LKinfo is formed. Conceptually, the sequence of steps are:

The calling arguments

Initialize the LKinfo object with the all the passed arguments.

Call: setDefaultsLKinfo

If necessary add some components that fill in specific default parameters for a given method.

Create basis information

The component LKinfo\$basisInfo holds the specific components for the basis functions. Currently there is not much flexibility here are as it is just derived from the calling arguments.

Call: LKrigSetupLattice

Setup the lattice information and put in the component LKinfo\$latticeInfo

Call: LKrigSetupAlpha

Setup the format for the alpha parameters based on what was passed in and out this in the component LKinfo\$alpha

Call: LKrigSetupAwght

Setup the format for the alpha parameters based on what was passed in and out this in the component LKinfo\$a.wght

Check using: LKinfoCheck

Check for the required components in the object.

The LKinfo object is a list with class LKinfo and comes with a print method to make it easier to read and also has the class of the geometry, e.g LKInterval for the 1-d model. The functions that are required for a geometry all take the LKinfo object as their first argument and the S3 object dispatching then calls the correct version. Put more simply, if LKinfoExample has class LKInterval then the code LKrigLatticeCenters(LKinfoExample, Level= 1) will result in computing the lattice centers at the first level and using the 1-d geometry. In this way much of this package can use "generic" functions for the computational steps without revealing the extra complications and exceptions for particular models and geometries. The structure of the LKinfo object is not completely rigid. Certain components of the list are required and these are checked by LKinfoCheck. Other components of this list can be geometry and problem dependent and this is a useful way to pass information to the computations. In general when faced with what should be included in LKinfo follow the style in LKInterval but also include additional components that are needed to define the model. Some examples are mentioned below.

Required methods

For any model there are three functions required for a new geometry. By a new geometry we mean any variation from the regular grids and beyond a second order SAR. As a first step choose a name for the new geometry, e.g. MyNewIdea, and this should also be the name passed as the argument LKGeometry in the LKrigSetup function. Three new functions need to be supplied and that end in .MyNewIdea and the user will recognize these as S3 methods being added to LatticeKrig.

The specific functions you will have to add to LatticeKrig package are:

LKrigSetupLattice.MyNewIdea(LKinfo) This function will take the information in the LKinfo list and create any indexing and do other computations to create the specific lattice of spatial points. The return value is a list with several basic components that are required. For example m is the total number of basis functions in the model. See help( LKrigSetupLattice) for a description of these required arguments. In addition one can add other arguments that help in coding the following two functions. It possible that this function could be avoided by doing all the setup whenever the SAR and the centers functions are used but we suspect this would result in more complicated code and be less modular. The 1-d implementation in the source file ModelInterval.R is a good example to get an idea of the whole process.

LKrigSAR.MyNewIdea(LKinfo, Level) Returns a sparse spatial autoregressive matrix at a specific level of the multi-resolution and in spind sparse matrix format. See help( LKrigSAR). Note that the lattice information from the previous function is in the component latticeInfo and so this function needs to access any lattice information in this way from the LKinfo object. e.g. LKinfo\$latticeInfo\$m are the total number of basis functions.

LKrigLatticeCenters.MyNewIdea(LKinfo, Level) Returns the spatial locations of the lattice points at a specific level of the multi-resolution. This can either be as a matrix with columns indexing dimension and rows indexing points or another class of object. These locations (or object) are subsequently passed to the second argument of the distance function LKrigDistance to create the basis functions.

For initially creating a new geometry one can just rely on this function returning the lattice locations. Once the whole geometry works one can consider returning a more specialized object and pairing these with more efficient distance finding algorithms. For example if the lattice points are on a regular, equally spaced grid the nearest neighbor locations can be found quickly just by simple arithmetic and this is much faster than a search. This is implemented for the LKInterval, LKBox, and LKRectangle geometries by creating the class gridList for the grid of lattice centers and passing this object directly into the second argument of the distance function LKrigDistance. See also help(LKrigLatticeCenters) and help(LKrigDistance)

Optional Methods

To make the modeling more flexible and take advantage of faster computation some additional methods can be added but are not required. Any new version should add your geometry name to the end. e.g. setDefaultsLKinfo should have the new name: setDefaultsLKinfo.MyNewIdea and have the default arguments listed in the generic function (setDefaultsLKinfo).

setDefaultsLKinfo This function takes the skeleton LKinfo list based on the arguments passed to LKrigSetup and can modify and add to them. The returned value is now taken as the initial LKinfo object. The default method is that no changes are made. But adding versions to this method can be very useful for specific cases (e.g. see setDefaultsLKinfo) without complicating the main function with lots of conditional statements.

LKrigSetupAlpha These methods convert the alpha information passed into LKrigSetup to a list format. A default is provided if you do not supply a specific version. See help(LKrigSetupAlpha).

LKrigSetupAwght These methods convert the a.wght information passed into LKrigSetup to a list format.

LKrigNormalizeBasisFast This is a method to implement any fast methods for normalizing the basis functions. There is no default method because it is assumed this is geometry dependent. See LKrigNormalizeBasisFast.LKRectangle as an example of how this is done. Note that there is already a "slow" method, LKrigNormalizeBasis involving a Cholesky solve of the precision matrix for each point where the basis function is evaluated.

LKinfoCheck This function checks the LKinfo object for necessary components and does some rudimentary comparisons to make sure arrays are the right size. It is always good to add more checks!

Author(s)

Doug Nychka

See Also

LKrigSetup, LKrigSAR, LKrigLatticeCenters

Description

This method performs some simple checks on the LKinfo object before it is returned by LKrigSetup. One benefit is that it checks for the essential components in each of the parts of the LKinfo object. Currently there is only a default method supplied.

Here is an example of how the checks work

	LKinfo<- LKrigSetup( cbind(c(0,1)), LKGeometry="LKInterval",
    a.wght=5, nlevel=2, nu=1.5, NC=4)
    LKinfo$alpha
# corrupt it  
    LKinfo$alpha <- list( 1,.5,.25)
try( LKinfoCheck( LKinfo) )

Usage

LKinfoCheck(object, ...)
## Default S3 method:
LKinfoCheck(object, ...)
## S3 method for class 'LKRectangle'
LKinfoCheck(object, ...)

Arguments

Value

There is no returned value. The benefit is the side effect of an error message and a "stop if the are problems in the LKinfo object.

Author(s)

Doug Nychka

See Also

LKrigSetup

Description

The basic LatticeKrig model is for a 2-d spatial domain and is identified by the geometry class LKRectangle. The lattices used in this case are equally spaced regular grids that are nested and the default distance function is standard Euclidean distance. The number of lattice points in the first level is NC inside the spatial domain and the subsequent levels decrease the lattice spacing by factors of 2. The spatial autoregressive coefficients, referred to as the a.wght parameter(s) have several levels of detail depending on the stationarity and isotropy. In the simplest case a.wght is the central value, and should be greater than 4. Also the four nearest neighbors are take to be -1. In the most complicated a.wght can include all 8 nearest neighbors can be specified differently at every lattice node. To be stable, the sum of the a.wght parameters for a node should be greater than zero.

Details

Here is a simple and small three level example that sets up the LatticeKrig model for spatial estimation and prediction. It assumes a spatial domain with extent [-1,1] in the horizontal and [-1,1] in vertical and beginning with 4 grid points in lowest level and with a default of 5 extra points outside the domain for boundary corrections. The assumed SAR model defaults to a central value of 4.1 and the 4 nearest weights are -1. Defining the extent of the spatial domain explicitly was done to simplify the example. Typically one just defines the domain to the minimum and maximum x and y locations of the spatial data (and the locations can passed as the first argument in the example below instead of finding ranges.) Note that does not mean that the spatial locations fill out a rectangle and they do not need to regularly spaced. The nu parameter is a handy way to specify the relative weights given each level. For a larger numbers of levels this parameter is equivalent to the Matern smoothness parameter.

Setting up the LKinfo object.

 sDomain<- cbind( c(-1,1), c( -1, 1))
LKinfo<- LKrigSetup(sDomain, nlevel=3, NC=4, a.wght=4.1, nu=1.0)

with the result

print(LKinfo)
Classes for this object are:  LKinfo LKRectangle
The second class usually will indicate the geometry
     e.g.  2-d rectangle is  LKRectangle
 
Ranges of locations in raw scale:
     [,1] [,2]
[1,]   -1   -1
[2,]    1    1
 
Number of levels: 3
delta scalings: 0.6666667 0.3333333 0.1666667
with an overlap parameter of  2.5
alpha:  0.7619048 0.1904762 0.04761905
based on smoothness nu =  1
a.wght:  4.1 4.1 4.1
 
Basis  type: Radial using  WendlandFunction  and 
Euclidean  distance.
Basis functions will be normalized
 
Total number of basis functions  1014
 Level Basis size      
     1        196 14 14
     2        289 17 17
     3        529 23 23
 
Lambda value:  NA

About the lattice The number of nodes define the number of basis functions and at first may seem a bit mysterious. However at the first level we get 4 lattice points with 5 extra boundary points added in the x direction amounting to 14 total and similarly 14 in the y direction because it is exactly of the same size. The default is the the parameter NC determines the number of lattice points in the larger dimension and the smaller dimension is divided according to the spacing from the longer dimension. For example if the y extent was [-1,.5] the number of lattice points would be spaced according to

delta scalings: 0.6666667 0.3333333 0.1666667

Thus seq( -1,.5, 0.6666667 ) are the (three) generated points within the spatial domain, and of course 5 would also added beyond each endpoint.

To query the LKinfo object this information is in the latticeInfo component. E.g. the first set of lattice locations.

LKinfo$latticeInfo$grid[[1]]
$x
 [1] -4.3333 -3.6667 -3.0000 -2.3333 -1.6667
 [6] -1.0000 -0.3333  0.3333  1.0000  1.6667
[11]  2.3333  3.0000  3.6667  4.3333

$y
 [1] -4.3333 -3.6667 -3.0000 -2.3333 -1.6667
 [6] -1.0000 -0.3333  0.3333  1.0000  1.6667
[11]  2.3333  3.0000  3.6667  4.3333

attr(,"class")
[1] "gridList"

By default equal spacing of the lattice is assumed in the x and y directions. To change this use the optional V argument to scale the coordinates. For example, to set 4 lattice points in both dimension for the above example:

sDomain2<- cbind( c(-1,1), c( -1, .5))
LKinfo<- LKrigSetup(sDomain2, nlevel=3, NC=4,
                 a.wght=4.1, nu=1.0,
                 V= diag(c(2, 1.5)) )

print(LKinfo$latticeInfo$mx)

     [,1] [,2]
[1,]   14   14
[2,]   17   17
[3,]   23   23

About the basis functions With the lattice points defined the default basis functions are radial Wendland functions centered at each point. The scaling of the basis functions is determined by the overlap and the delta values. In R code, if the lattice point is x0 at level l then the basis function evaluated at location x1 is

basisFunctionValue <- Wendland( rdist( x1,x0)/( delta[l]*overlap))

If V is included then x1 is transformed as x1%*%solve(V) before evaluating in the basis function.

About a.wght

For this geometry the basic form of awght is as a list with as many components as levels. However, the LKrigSetup function will reshape a scalar or vector argument in this this format. a.wght can take the following forms:

Scalar value: In this case the value is used at all lattice points and at all levels in the SAR. Four nearest neighbors are set to -1.

List/vector of length nlevel In this case separate values for the a.wght will be used for the central SAR value at each level.

List of vectors With a list of length nlevel and each component is a 9 element vector, the values in the vector correspond to the central lattice point and 8 nearest neighbors with the indexing:

      1 4 7
      2 5 8
      3 6 9

E.g. the 5 element is the center, 3 is the lower left hand corner etc.

Non-stationary models This is the nuclear option to handle a completely non-stationary correlation structure! In this case one specifies the models by passing prediction objects for one or more of sigma2.object, alphaObject or a.wghtObject.

Author(s)

Doug Nychka

See Also

LatticeKrig LKGeometry LKrig

Examples

# the grid with only 2 extra boundary points
  sDomain<- cbind( c(-1,1), c( -1, 1))
  LKinfo<- LKrigSetup(sDomain, nlevel=3, NC=4, a.wght=4.1,
           NC.buffer=2, alpha=c(1,.5,.125) )
  LKgrid<- LKinfo$latticeInfo$grid
  plot(   make.surface.grid(LKgrid[[1]]),
                          pch=16, cex=1.5)
  points( make.surface.grid(LKgrid[[2]]), 
                          pch=15, cex=.8, col="red" )
  points( make.surface.grid(LKgrid[[3]]),
                          pch="+", col="green" )
  rect(sDomain[1,1],sDomain[1,2],
     sDomain[2,1],sDomain[2,2], lwd=3 )

# basis functions on a grid
# this function actually evaluates all of them on the grid.
  xg<- make.surface.grid(
        list(x=seq( -2,2,,80), y=seq( -2,2,,80)) )
  out<- LKrig.basis( xg, LKinfo)
# basis functions 20, 26, 100  and 200
  plot(   make.surface.grid( LKgrid[[1]] ) , 
                          pch=16, cex=.5)
  rect(sDomain[1,1],sDomain[1,2],
     sDomain[2,1],sDomain[2,2], lwd=3,border="grey" )
  contour( as.surface(xg, out[,20]), col="red1",
                                     add=TRUE)
  contour( as.surface(xg, out[,36]), col="red4", 
                                     add=TRUE)
  contour( as.surface(xg, out[,100]), col="blue1",
                                     add=TRUE)
  contour( as.surface(xg, out[,200]), col="blue4",
                                      add=TRUE)
  title( "basis functions 20, 26, 100, 200")

Description

A variation of Kriging with fixed basis functions that uses a compactly supported covariance to create a regular set of basis functions on a grid. The coefficients of these basis functions are modeled as a Gaussian Markov random field (GMRF). Although the precision matrix of the GMRF will be sparse the model can still exhibit longer ranges of spatial dependence. The multi-resolution feature of this model allows for the approximation of a wide variety of standard covariance functions. There are also some simple extensions where this function can be used to solve a linear inverse type problem (see X and U below).

Usage

LKrig( x, y, weights = NULL, Z = NULL, LKinfo = NULL, iseed =
                 NA, NtrA = 20, use.cholesky = NULL, return.cholesky =
                 TRUE, X = NULL, U = NULL, wX = NULL, wU = NULL,
                 return.wXandwU = TRUE,
                 ..., getVarNames = TRUE, verbose = FALSE)

## S3 method for class 'LKrig'
predict(object, xnew = object$x, Znew = NULL, drop.Z = FALSE,
                 just.fixed = FALSE, return.levels = FALSE,
                 collapseFixedEffect=object$collapseFixedEffect, ...)
                 
## S3 method for class 'LKrig'
predictSE(object, xnew = NULL, Znew = object$Z, verbose = FALSE,
                 ...)

## S3 method for class 'LKrig'
surface( object, ...)

## S3 method for class 'LKrig'
predictSurface(object, grid.list = NULL, extrap = FALSE,
 chull.mask = NA, nx = 80, ny = 80, xy = c(1, 2), verbose = FALSE,
                 ZGrid = NULL, drop.Z = FALSE, ...)

## S3 method for class 'LKrig'
print( x, digits=4, ...)

## S3 method for class 'LKrig'
summary( object, digits=4, stripAwght = TRUE,...)

createLKrigObject(x, y, weights = NULL, Z, X, U, LKinfo,
                 xName = "xVar", ZName = "ZVar", UName = "UVar",
                 verbose = FALSE)

Arguments

Details

This method combines compactly supported basis functions and a Markov random field covariance model to provide spatial analysis for large data sets. The model for the spatial field (or spatial process) is

f(x) = N(x) + Z d + sum Phi.j(x) c.j.

x is a location in two dimensions, N(x) is a low order (linear) polynomial, Z is a matrix of spatial covariates and d a coefficient vector. Phi.j for 1 <= j <= m is a set of fixed basis functions and c.j the coefficients. The variance of f(x) is given by the parameter sigma2 throughout this package. As explained below the process f is a sum of nlevel independent processes that have different scales of spatial dependence. The alpha gives the relative weighting between these processes. Thus, the minimum set of parameters needed to describe the covariance of f are the integer NC, two scalars sigma2 and a.wght, and a vector alpha. alpha has a length of the number of multi-resolution levels but we recommend that it be constrained sum to one. The parameter sigma2 is the marginal variance of the process and this multiples alpha. Thus in total there are a 1 + 2 + (nlevel - 1) parameters for a minimal specification of the covariance. Note that this parsimonious specification results in a covariance that is close to being stationary and isotropic when normalize is TRUE. An additional constraint on alpha is to make the weights alpha[j] proportional to exp( - 2*j*nu) where nu controls decay of the alpha weights. There is some theory to suggest that nu is analogous to the smoothness parameter from the Matern family (e.g. nu=.5 approximates the exponential). In this case the covariance model requires just four parameters, NC, sigma2, a.wght, nu.

The data model is

Y.k = f(x.k) + e.k

Y.k are (scalar) observations made at spatial locations x.k with e.k uncorrelated normal errors with variance tau^2/weights. Thus there is a minimum of one new parameter from the data model: tau. Note that prediction only depends on the ratio lambda = tau^2/ sigma2 and not surprisingly lambda plays a key role in specifying and fitting a spatial model. Also taken with the model for f the minimum parameters needed for a spatial prediction are still four NC, a.wght, nu and lambda. For fixed lambda there are closed form expressions for the MLEs for tau and sigma2. LKrig exploits this feature by depending on lambda and then computing the MLEs for tau and sigma2.

The data model can also be written in vector form as

Y = T d + PHI c + e

Where T is a matrix that combines the low order polynomial with other possible covariates and PHI is the matrix basis functions evaluated at the observations. A simple extension is considering a linear inverse problem form. This is an experimental of LatticeKrig that is intended for solving large linear inverse problems. In this case one supplies to the LKrig function two matrices (in spam sparse format) such that

Y = U d + X c + e.

This features allows for observations that are linear functionals of f and not necessarily the function evaluated at the observation locations. This model is useful for working with demographic or tomographic problems where the observations are expressed as particular integrals of f over different regions or different projections. See the last example on how to use this option. If U and X are not supplied the default is to consider a model with observations made at evaluating f at the observation locations.

About dense matrix computations: The value of this package is to handle larger spatial data sets by exploiting sparse matrix methods using the spam package. However, for small problems, checking, or timing it is useful to do the computations using the usual dense matrix decomposition. If the component dense in the LKinfo object is TRUE then dense (i.e. standard) methods will be used. Set this switch from the LKrigSetup function.

Spatial prediction: The basis functions are assumed to be fixed and the coefficients follow a multivariate Gaussian distribution. Given this spatial model for f, it is possible to compute the conditional expectation of f given Y and also maximize the likelihood for the model parameters, lambda, alpha, and a.wght. This setting is known as fixed rank Kriging and is a common strategy for formulating a spatial model. Typically fixed rank Kriging is used to reduce the dimension of the problem by limiting the number of basis functions. We take a different approach in allowing for models that might even have more basis functions than observations. This provides a spatial model that can come close to interpolating the observations and the spatial process is represented with many degrees of freedom. The key is to make sure the model ingredients result in sparse matrices where linear algebra is required for the computations. By doing so in this package it is possible to compute the estimates and likelihood for large numbers of spatial locations. This model has an advantage over covariance tapering or compactly supported covariance functions (e.g. fastTps from fields), because the implied covariance functions can have longer range correlations.

Radial basis functions ( Phi.j ) : The basis functions are two-dimensional radial basis functions (RBFs) that are derived from scaling and translations of a single covariance function. The default in LatticeKrig is to use the Wendland compactly supported stationary covariance (order 2 for 2 dimensions) that is scaled to be zero beyond a distance of 1. For d the distance between spatial locations, this Wendland function has the standard form:

(1 - d)^6 * (35 * d^2 + 18 * d + 3))/3 for d in [0,1]

0 otherwise.

For a single level the RBFs are centered at a regular grid of points and with radial support delta*overlap where delta is the spacing between grid points. We will also refer to this grid of centers as a lattice and the centers are also referred to as "nodes" in the RBF literature. The overlap for the Wendland has the default value of 2.5 and this represents a compromise between the number of nonzero matrix elements for computation and the shape of the covariance functions.

To create a multi-resolution basis, each subsequent level is based on a grid with delta divided by 2. See the example below and help(LKrig.basis) for more details. For multiple levels the basis functions can be grouped according to the resolution levels and the coefficients can be grouped in a similar manner. There is one important difference in the basis construction – a normalization – and this aspect makes it different from a simple radial basis function specification and is described below.

Markov random field (GMRF) for the coefficients (c.j) : Because the coefficients are identified with locations on a lattice it is easy to formulate a Markov random field for their distribution based on the relationship to neighboring lattice points. The distribution on the basis function coefficients is a multivariate normal, with a mean of zero and the the precision matrix, Q, (inverse of Q is the covariance matrix). Q is partitioned in a block diagonal format corresponding to the organization of the basis functions and coefficients into levels of resolution. Moreover, coefficients at different levels are assumed to be independent and so Q will be block diagonal. If nlevels are specified, the ith block has a precision matrix based on a spatial autoregression with a.wght[i] being related to the spatial autoregressive parameter(s). Schematically in the simplest case the weighting for an interior lattice point with its four neighbors is

The fundamental idea is that these weights applied to each point in the lattice will result in a lattice of random variables that are independent. The specific precision matrix for each block (level), Q.i, is implemented by LKrig.MRF.precision. In the case when alpha is a scalar, let C.i be the vector of basis coefficients at the ith level then we assume that B %*% C.i will be independent N(0,1) random variables. By elementary properties of the covariance it follows that the precision matrix for C.i is Q.i= t(B)%*%B. Thus, given B one can determine the precision matrix and hence the covariance matrix. Each row of B, corresponding to a point in the lattice in the interior, is "a" (a.wght[i]) on the diagonal and -1 at each of the four nearest neighbors of the lattice points. Points on the edges and corners just have less neighbors but get the same weights.

This description is a spatial autoregressive model (SAR). The matrix Q will of course have more nonzero values than B and the entries in Q can be identified as the weights for a conditional autoregressive model (CAR). Moreover, the CAR specification defines the neighborhood such that the Markov property holds. Values for a.wght[i] that are greater than 4 give reasonable covariance models. Moreover setting a.wght[i] to 4 and normalize to FALSE in the call to LKrig will give a thin-plate spline type model that does not have a range parameter. An approximate strategy, however, is to set a.wght close to 4 and get some benefit from the normalization to reduce edge effects.

Multi-resolution process Given basis functions and coefficients at each level we have defined a spatial process g.i that can be evaluated at any location in the domain. These processes are weighted by the parameter vector alpha and then added together to give the full process. It is also assumed that the coefficients at different resolution levels are independent and so the processes at each level are also independent. The block diagonal structure for Q does not appear to limit how well this model can approximate standard spatial models and simplifies the computations. If the each g.i is normalized to have a marginal variance of one then g will have a variance that is the sum of the alpha parameters. Usually it is useful to constrain the alpha parameters to sum to one and then include an additional variance parameter, sigma2, to be the marginal variance for g. So the full model for the spatial process used in LatticeKrig is

g(x) = sqrt(sigma2) * sum.i sqrt(alpha[i]) * g.i(x)

The specification of the basis and GMRF is through the components of the object LKinfo, a required component for many LatticeKrig functions. High level functions such as LKrig only require a minimal amount of information and combined with default choices create the LKinfo list. One direct way to create the complete list is to use LKrigSetup as in the example below. For nlevel==1 one needs to specify a.wght, NC, and also lambda related to the measurement error variance. For a multi-resolution setup, besides these parameters, one should consider different values of alpha keeping in mind that if this vector is not constrained in some way ( e.g. sum(alpha)==1) it will not be identifiable from lambda.

The covariance derived from the GMRF and basis functions: An important part of this method is that the GMRF describes the coefficients of the basis functions rather than the field itself. Thus in order to determine the covariance for the observed data one needs to take into account both the GMRF covariance and the expansion using the basis functions. The reader is referred to the function LKrig.cov for an explicit code that computes the implied covariance function for the process f. Formally, if P is the covariance matrix (the inverse of Q) for the coefficients then the covariance between the field at two locations x1 and x2, will be

sum_ij P_ij Phi.i(x1) Phi.j(x2)

Moreover, under the assumption that coefficients at different levels are independent this sum can be decomposed into sums over the separate levels. The function LKrig.cov evaluates this formula based on the LKrig object (LKinfo) at arbitrary groups of locations returning a cross covariance matrix. LKrig.cov.plot is a handy function for evaluating the covariance in the x and y directions to examine its shape. The function LKrig.cov is also in the form to be used with conventional Kriging codes in the fields package (loaded by LatticeKrig) mKrig or Krig and can be used for checking and examining the implied covariance function.

Normalizing the basis functions The unnormalized basis functions result in a covariance that has some non-stationary artifacts (see example below). For a covariance matrix P and for any location x one can evaluate the marginal variance of the process using unnormalized basis functions for each multi-resolution level. Based this computation there is a weighting function, say w.i(x), so that when the normalized basis w.i(x) Phi.i(x) is used the marginal variance for the multi-resolution process at each level will be one. This makes the basis functions dependent on the choice of Q and results in some extra overhead for computation. But we believe it is useful to avoid obvious artifacts resulting from the use of a finite spatial domain (edge effects) and the discretization used in the basis function model. This is an explicit way to make the model stationary in the marginal variance with the result that the covariance also tends to be closer to a stationary model. In this way the discretization and edges effects associated with the GMRF can be significantly diminished.

The default in LKrig is normalize = TRUE. It is an open question as to whether all levels of the multi-resolution need this kind of explicit normalization. There is the opportunity within the LKrig.basis function to only normalize specific levels with the normalize being extended from a single logical to a vector of logicals for the different levels. To examine what these edge effect artifacts are like the marginal variances for a 6 by 6 basis is included at the end of the Examples Section.

Non-stationary and anisotropic modifications to the covariance The simplest way to build in departures from isotropy is to scale the coordinates. The device to specify this transformation is including the V argument in setting up the LKinfo object or passed in the ... for this function. V should be a matrix that represents the transposed, inverse transformation applied to the raw coordinates. For example if x1 are locations d- dimensions ( i.e. x1 has d columns) and V is a dXd matrix then the transformed coordinates are

x1Transformed <- x1 %*% t(solve(V))

and all subsequent computations will use the transformed coordinates for evaluating the basis functions. This also means that the lattice centers are locations in the transformed scale. The following example might be helpful:

set.seed(123)
x<- matrix( runif( 200),100,2)
V<- cbind( c( 2,1), c( -1,1) )
y<-  x[,1] + x[,2]
LKinfo<- LKrigSetup( x, V=V, nlevel=1, a.wght=5, NC=20)

gridCenters<- LKrigLatticeCenters( LKinfo, Level=1)
# grid in transformed space:
plot( make.surface.grid( gridCenters))
# transformed data:
points( x%*%solve(t(V)), col="red", pch=16) 
# the model is fit to these locations. 

Keep in mind that a practical use of this argument is just to scale a variable(s) that is V, a diagonal matrix with diagonal elements being the values to scale the individual coordinates.

Given that the process at each level has been normalized to have marginal variance of one there are several other points where the variance can be modified. The variance at level i is scaled by the parameters alpha[i] and the marginal variance of the process is scaled by sigma2. Each of these can be extended to have some spatial variation and thus provide a model for non-stationarity.

An option in specifying the marginal variance is to prescribe a spatially varying multiplier. This component is specified by the object sigma2.object. By default this is not included (or assumed to be identically one) but, if used, the full specification for the marginal variance of the spatial process at location x is formally: sigma2 * predict(sigma2.object, x) * sum( alpha) There is then a problem of identifiability between these and a good choice is to constrain sum(alpha) ==1 so that sigma2* predict(sigma2.object, x) is associated with the marginal variance of the full spatial process.

A second option is to allow the alpha variance component parameters to vary across space and at each level. The model in this case could be

g(x) = sqrt(sigma2(x)) * sum.i sqrt(alpha(x).i) * g.i(x)

To specify this case alpha should be a list with nlevel components and each component being a "predict" object that will evaluate the alpha variance at a given level at arbitrary spatial locations. Explicitly alpha should be set up so that predict(alpha[[l]], x1) will return a vector of values for alpha at level l and at locations x1. This happens in LKrig.basis and the user is referred to this function to see how it is handled and possibly to make modifications to the model. Similar to the scalar alpha case we recommend that the alpha's across levels sum to one and it still may make sense to have the sigma2 parameter vary across space as well. Here alpha would control a varying correlation range and sigma2 the marginal variance.

LKrig also has the flexibility to handle more general weights in the GMRF. This is accomplished by a.wght being a list with as many components as levels. If each component is a vector of length nine then these are interpreted as the weights to be applied for the lattice point and its 8 nearest neighbors see the commented source code in LKrigSAR.LKRectangle for details. If each component is a matrix then these are interpreted as the (non-stationary) center weights for each lattice point. Finally if the component is an array with three dimensions this specifies the center and 8 nearest neighbors for every point in the lattice. At this point the choice of these weights beyond a stationary model is experimental and we will defer further documentation of these features to a future version of this package.

The smoothing parameter lambda and effective degrees of freedom Consistent with other fields package functions, the two main parameters in the model, tau^2 and sigma2 are parameterized as lambda = tau^2/sigma2 and sigma2. The MLEs for sigma2 and tau can be written in closed form as a function of lambda and these estimates can be substituted into the full likelihood to give a concentrated version that can numerically be maximized over lambda. The smoothing parameter lambda is best varied on a log scale and is sometimes difficult to interpret independent of the particular set of locations, sample size, and covariance. A more useful interpretation of lambda is through the effective degrees of freedom and an approximate value is found by default using a Monte Carlo technique. The effective degrees of freedom will vary with the dimension of the fixed regression part of the spatial model ( typical 3 = constant + linear terms) and the total number of observations. It can be interpreted as the approximate number of "parameters" needed to represent the spatial prediction surface. For a fixed covariance model the spatial estimate at the observation locations can be represented as f hat = A(lambda) y where A is a matrix that does not depend on observations. The effective number of degrees of freedom is the trace of A(lambda) in analogy to the least squares regression "hat" matrix and has an inverse, monotonic relationship to lambda. The Monte Carlo method uses the fact that if e are iid N(0,1) E( t(e) A(lambda) e) = trace( A(lambda)).

Descriptions of specific functions and objects:

LKrig: Find spatial process estimate for fixed covariance specified by nlevel, alpha, a.wght, NC, and lambda or this information in an LKinfo list.

predict.LKrig, predictSE.LKrig: These functions evaluate the model at the the data locations or at xnew if it is included. Note the use of the drop.Z argument to either include the covariates or just restrict the computation to the spatial drift and the smooth component. If drop.Z is FALSE then then Znew needs to be included for predictions off of the observation locations. The standard errors are computed under the assumption that the covariance is known, that it is the TRUE covariance for the process, and both the process and measurement errors are multivariate normal. The formula that is used involves some recondite shortcuts for efficiency but has been checked against the standard errors found from an alternative formula in the fields Krig function. (See the script Lkrig.se.tests.R in the tests sub-directory for details.)

Value

LKrig: A LKrig class object with components for evaluating the estimate at arbitrary locations, describing the fit and as an option (with Mc.return=TRUE) the Cholesky decomposition to allow for fast updating with new values of lambda, alpha, and a.wght. The "symbolic" first step in the sparse Cholesky decomposition can also be used to compute the sparse Cholesky decomposition for a different positive definite matrix that has the same pattern of zeroes. This option is useful in computing the likelihood under different covariance parameters. For the LKrig covariance the sparsity pattern will be the same if NC, level, overlap and the data locations x are kept the same. The returned component LKinfo has class LKinfo and is a list with the information that describes the layout of the multi-resolution basis functions and the covariance parameters for the GMRF. (See help(LKinfo) and also LK.basis as an example.)

predict.LKrig, predictSE.LKrig: A vector of predictions or standard errors.

predictSurface.LKrig: A list in image format (i.e. having components x,y,z) of the surface evaluated on a regular grid. This surface can then be plotted using several different R base package and fields functions e.g. image, image.plotcontour, persp, drape.plot. The surface method just calls this function and then a combination of the image and contour plotting functions.

Author(s)

Doug Nychka

See Also

LatticeKrig, LKrig.sim.conditional, LKrig.MLE, mKrig, Krig, fastTps, Wendland, LKrig.coef, Lkrig.lnPlike, LKrig.MRF.precision, LKrig.precision

Examples

# NOTE: CRAN requires that the "run" examples  execute within  5 seconds. 
# so to meet this constraint many of these have been 
# severely limited to be quick, typically making NC and nlevel
# very small.
	
# Load ozone data set
  data(ozone2)  
  x<-ozone2$lon.lat
  y<- ozone2$y[16,]
# Find location that are not 'NA'.
# (LKrig is not set up to handle missing observations.)
  good <-  !is.na( y)
  x<- x[good,]
  y<- y[good]

# fairly arbitrary choices for covariance parameters and lambda
# just to show a basic level call
  obj1<- LKrig( x,y, a.wght=5, nlevel=3, nu=1.0, NC=10, lambda=.1)
  print(obj1)
  
# estimates of fixed model parameters with their SEs
  summary( obj1)$coefficients
#  
# this is the same as:
#  LKinfoEX<- LKrigSetup( x, a.wght=5, nlevel=3, nu=1.0, NC=4)
#  obj1<- LKrig( x,y, LKinfo= LKinfoEX, lambda=.1)  

# thin plate spline-like model with the lambda parameter estimated by
# maximum likelihood. Default choices are made for a.wght, nlevel, NC
# and alpha.
## Not run: 
  obj<- LatticeKrig( x, y)
# 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)

## End(Not run)

# breaking down the LatticeKrig function into several steps. 
# also use a different covariance model that has fewer basis functions
# (to make the example run more quickly)

## Not run:  
  LKinfo<- LKrigSetup( x, nlevel=3, nu=1, NC=5, a.wght=5,
                        lambda=.01)
# maximize likelihood over lambda see help( LKrig.MLE) for details
# this assumes the value of 5 for a.wght.  In many cases the fit is not
# very sensitive to the range parameter such as a.wght in this case --
# but very sensitive to lambda when varied on a log scale.

  MLE.fit<- LKrig.MLE(x,y, LKinfo=LKinfo)
  MLE.fit$summary # summary of optimization over lambda.
# fit using MLE for lambda MLE function has added MLE value of lambda to
# the LKinfo object.
  obj<- LKrig( x,y, LKinfo=MLE.fit$LKinfo.MLE)  
  print( obj)  
# find prediction standard errors at locations based on fixing covariance
# at MLE's
  out.se<- predictSE( obj, xnew= x)
# one could evaluate the SE on a grid to get the surface of predicted SE's 
# for large grids it is better to use LKrig.sim.conditional to estimate
#  the variances by Monte Carlo

## End(Not run)

##########################################################################
# Use multi-resolution model that approximates an exponential covariance
# Note that a.wght, related to a range/scale parameter
# is specified at a (seemingly) arbitrary value. 
##########################################################################
  
  LKinfo<- LKrigSetup( x, NC=4, nu=1, nlevel=2, a.wght= 5)
# take a look at the implied covariance function solid=along x
#  and  dashed=along y 
  check.cov<- LKrig.cov.plot( LKinfo)
  matplot( check.cov$d, check.cov$cov, type="l", lty=c(1,2))  

############################################################################
# Search over lambda to find MLE for ozone data with approximate exponential
# covariance
###########################################################################
## Not run: 
  LKinfo.temp<-  LKrigSetup( x, NC=6, nu=1,  nlevel=3, a.wght= 5, lambda=.5)
# starting value for lambda optimization 
  MLE.search<- LKrig.MLE(x,y, LKinfo=LKinfo.temp)
# this function returns an LKinfo object with the MLE for lambda included.
  MLE.ozone.fit<- LKrig( x,y,  LKinfo= MLE.search$LKinfo.MLE)

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

  data(COmonthlyMet)
  y.CO<- CO.tmin.MAM.climate
  good<- !is.na( y.CO)
  y.CO<-y.CO[good]
  x.CO<- as.matrix(CO.loc[good,])
  Z.CO<- CO.elev[good]
# single level with large range parameter -- similar to a thin plate spline
#  lambda specified 

# fit with elevation 
# note how for convenience the LKrig parameters are
# just included here and not passed as a separate LKinfo object. 
  obj.CO.elev<- LKrig( x.CO,y.CO,Z=Z.CO, nlevel=1, NC=50, alpha=1, lambda=.005,
                          a.wght=4.1)
# BTW the coefficient for the linear term for elevation  is obj.CO$d.coef[4]
# fitted surface without the elevation term
## Not run: 
   LKinfo<- LKrigSetup( x.CO, nlevel=1, NC=20,alpha=1, a.wght=4.1, lambda=1.0)
# lambda given here is just the starting value for MLE optimization
  CO.MLE<- LKrig.MLE( x.CO,y.CO,Z=Z.CO, LKinfo=LKinfo)
  obj.CO.elev<- LKrig( x.CO,y.CO,Z=Z.CO, LKinfo= CO.MLE$LKinfo.MLE)
  CO.surface2<- predictSurface( obj.CO.elev, drop.Z=TRUE, nx=50, ny=50)
# pull off CO elevations at same locations on grid as the surface
  data( RMelevation) 
# a superset of elevations at 4km resolution
  elev.surface<- interp.surface.grid( RMelevation, CO.surface2)
   CO.full<- predictSurface( obj.CO.elev, ZGrid= elev.surface, nx=50, ny=50)
   
# for comparison fit without elevation as a linear covariate:
  CO.MLE2<- LKrig.MLE( x.CO,y.CO, LKinfo=LKinfo)
  obj.CO<- LKrig( x.CO,y.CO, LKinfo= CO.MLE2$LKinfo.MLE)
# surface estimate
  CO.surface<- predictSurface( obj.CO, nx=50, ny=50)
  set.panel( 2,1)
  coltab<- topo.colors(256)
  zr<- range( c( CO.full$z), na.rm=TRUE)
  image.plot( CO.surface,  col=coltab, zlim =zr)
    US( add=TRUE,lwd=2)
    title( "MAM min temperatures without elevation")
  image.plot( CO.full, col=coltab, zlim=zr)
    title( "Including elevation")
    US( add=TRUE,lwd=2)

## End(Not run)
## Not run: 
#### a 3-d example using alternative geometry
set.seed( 123)
N<- 500
x<-  matrix( runif(3* N,-1,1), ncol=3, nrow=N)
y<-   10*exp( -rdist( x, rbind( c(.5,.5,.6) ) )/.5)
LKinfo<- LKrigSetup( x,  
                 LKGeometry = "LKBox",
                     nlevel = 1,
                     a.wght = 6.01,
                         NC = 5,
                  NC.buffer = 2,
                  normalize = TRUE,
             choleskyMemory = list(nnzR= 2e6),
                    lambda = .1,
                    mean.neighbor=200
                   )  
# NOTE: memory for the spam routines has been increased to 
# avoid warnings  
  out1<- LKrig( x,y, LKinfo=LKinfo)

## End(Not run)  
  
## Not run: 
# MLE search over lambda and  a bigger problem
# but no normalization
N<- 5e4
x<-  matrix( runif(3* N,-1,1), ncol=3, nrow=N)
y<-   10*exp( -rdist( x, rbind( c(.5,.5,.6) ) )/.5)
LKinfo<- LKrigSetup( x,  
                 LKGeometry = "LKBox",
                     nlevel = 1,
                     a.wght = 6.01,
                         NC = 20,
                  NC.buffer = 2,
                  normalize = FALSE,
             choleskyMemory = list(nnzR= 25e6),
                    lambda = .1,
                    mean.neighbor=200
                   ) 
# add in timing                    
 system.time(  out1<- LKrig( x,y, LKinfo=LKinfo) ) # about 25 seconds
# LKinfo$normalize<- TRUE
# system.time(  out2<- LatticeKrig( x,y, LKinfo=LKinfo) )# ~250 seconds

## End(Not run)  
########################################################################
# for a more elaborate search over  a.wght, alpha and lambda to find
# joint MLEs see help(LKrig.MLE)
########################################################################

########################################################################
# A bigger problem: 26K observations and 4.6K basis functions
# fitting takes about 15 seconds on a laptop for a fixed covariance
#  LKrig.MLE to find the MLE (not included) for lambda takes about
#  8 minutes
#######################################################################
## Not run: 
  data(CO2)
  # the Ring geometry will be periodic in the first dimension and rectangular on 
  # second. A useful approximation for spherical data omitting the polar caps. 
  
  LKinfo.CO2<- LKrigSetup(CO2$lon.lat, NC=100,nlevel=1, lambda=.2,
                       a.wght=5, alpha=1, 
                       LKGeometry="LKRing", choleskyMemory = list(nnzR=2e6) )
  print(LKinfo.CO2)                                          
  obj1<- LKrig( CO2$lon.lat,CO2$y,LKinfo=LKinfo.CO2)
# 5700+ basis functions 101X57 lattice  (number of basis functions
# reduced in y direction because of a rectangular domain
  obj1$trA.est # about 2900+ effective degrees of freedom 
#
  glist<- list( x= seq( -180,180,,200),y=seq( -80,80,,100) )
  fhat<- predictSurface( obj1,grid.list=glist)
#Plot data and gap-filled estimate
  set.panel(2,1)
  quilt.plot(CO2$lon.lat,CO2$y,zlim=c(373,381))
  title("Simulated CO2 satellite observations")
  world(add=TRUE,col="magenta")
  image.plot( fhat,zlim=c(373,381))
  world( add=TRUE, col="magenta")
  title("Gap-filled global predictions")

## End(Not run)  
 set.panel()
#########################################################################
#  Comparing LKrig to ordinary Kriging
########################################################################

# Here is an illustration of using the fields function mKrig with the
# LKrig covariance to reproduce the computations of LKrig. The
# difference is that mKrig can not take advantage of any sparsity in
# the precision matrix because its inverse, the covariance matrix, is
# not sparse.  This example reinforces the concept that LKrig finds the
# the standard geostatistical estimate but just uses a particular
# covariance function defined via basis functions and the precision
# matrix.
# Load ozone data set (AGAIN)
## Not run: 
  data(ozone2)  
  x<-ozone2$lon.lat
  y<- ozone2$y[16,]
# Find location that are not 'NA'.
# (LKrig is not set up to handle missing observations.)
  good <-  !is.na( y)
  x<- x[good,]
  y<- y[good]
  a.wght<- 5
  lambda <-  1.5
  obj1<- LKrig( x,y,NC=16,nlevel=1, alpha=1,  lambda=lambda, a.wght=5,
                NtrA=20,iseed=122)
 
# in both calls iseed is set the same so we can compare 
# Monte Carlo estimates of effective degrees of freedom
  obj1$trA.est
# The covariance "parameters" are all in the list LKinfo
# to create this special list outside of a call to LKrig use
  LKinfo.test <- LKrigSetup( x, NC=16, nlevel=1, alpha=1.0,  a.wght=5)

# this call to mKrig should be identical to the LKrig results
# because it uses the LKrig.cov covariance with all the right parameters.
  obj2<- mKrig( x,y, lambda=lambda, m=2, cov.function="LKrig.cov",
                      cov.args=list( LKinfo=LKinfo.test), NtrA=20, iseed=122)
# compare the two results this is also an
# example of how tests are automated in fields
# set flag to have tests print results
  test.for.zero.flag<- TRUE
  test.for.zero( obj1$fitted.values, obj2$fitted.values,
                  tag="comparing predicted values LKrig and mKrig")
# compare standard errors. 
  se1<- predictSE.LKrig( obj1)
  se2<- predictSE.mKrig(obj2)
  test.for.zero( se1, se2,
                  tag="comparing standard errors for LKrig and mKrig")

## End(Not run)

## Not run: 
##################################################################
#  a linear inverse problem 
##################################################################
# NOTE the settings in the model are small to make for a fast example
#
data(ozone2)
x<- ozone2$lon.lat
y<- ozone2$y[16,]
good<- !is.na(y)	
x<- x[good,]
y<- y[good]

LKinfo<- LKrigSetup(x, a.wght=4.5, nlevel=3, nu=1, NC=4, lambda=.01)

# now observe a linear combination
  NNDist<- LKDist(x,x, delta=1.5) 
  A<- NNDist
  A$ra<- exp(-NNDist$ra)
# A is a weight matrix based on neighbors close by and
# in spind sparse matrix format
# now convert to spam format
  A<- spind2spam(A)
  
TMatrix <- get(LKinfo$fixedFunction)(x = x)
# Tmatrix is a 3 column matrix of constant and the two spatial coordinates
#  i.e. a linear function of the spatial variables 
PHI<- LKrig.basis(x, LKinfo)
X<-  A%*% PHI
U<-  A%*%TMatrix 
yIndirect<- A%*%y
#
# X<-  A

out0<- LatticeKrig(x,y, LKinfo=LKinfo)
out1<- LatticeKrig(x,yIndirect, U=U, X=X, LKinfo=LKinfo)

# the predict function evaluates f in this case -- not the fitted values that
# are related to the 
# observations
# partial agreement but not exact -- this is because the observational models
# assume different errors
#
plot( predict(out0,x), predict( out1,x))

out<- LKrig.sim.conditional( out1,M=5, nx=10, ny=10 )

## End(Not run)

Description

Some internal functions for LKrig that estimate the coefficients of the basis functions and compute the likelihood.

Usage

LKrigMakewU(object, verbose = FALSE)
LKrigMakewX(object, verbose = FALSE)
LKrig.coef( GCholesky, wX, wU, wy, lambda,
collapseFixedEffect = FALSE,  verbose = FALSE)
LKrig.lnPlike( GCholesky, Q, quad.form, nObs, nReps, weights, LKinfo)
LKrig.traceA(GCholesky, wX, wU, lambda, weights, NtrA, iseed = NA)
LKrigUnrollZGrid( grid.list, ZGrid=NULL)
LKDefaultVarNames(A, tag)

Arguments

Details

The LatticeKrig article can be used as a reference for the matrix computations and the G matrix from those formulas figures prominently. The GCholesky object in these functions is the cholesky decomposition of this matrix. For compatibility with older version of this package this object may also be named as Mc ( Cholesky of the M matrix) but the user should not identify this M with that in the article. Ideally all coding using Mc should be changed to GCholesky.

createLKrigObject Based on the arguments passed into LKrig forms the prototype LKrig object. This object is added to as one computes additional steps in the LKrig function. The Names argument in the call is awkward device to pass the names of the original x, Z and U objects when substitued these names are used as the default prefixs for column names of these matrices.

LKrigMakewU and LKrigMakewX construct the weighted U and X matrices from what is passed. In the case of observations that are point locations wU is found the weights and using the fixedFunction and wX is found from the weights and the multi-resolution basis functions. Note that X and wX are assumed to be in spam sparse matrix format.

LKrig.coef and LKrig.lnPlike are two low level functions to find the basis function coefficients and to evaluate the likelihood. The coefficients (c.mKrig) are also found because they provide for shortcut formulas for the standard errors and MLE estimates. These coefficients are identical to the basis coefficients (c.coef) found for usual Kriging in the mKrig function. LKrig.lnPlike also finds the profile MLE of tau and sigma2 given a fixed value for lambda (and alpha and a.wght). See the source for LKrig and also MLE.LKrig to see how these functions are used.

LKrig.traceA finds an estimate of the effective degrees of freedom of the smoothing matrix based a simple Monte Carlo scheme. The smoothing matrix A is the matrix for fixed covariance parameters so that y.hat = A y, where y.hat are the predicted values at the data locations. trace(A) is the effective degrees of freedom. If e are iid N(0,1) then the expected value of t(e)% * % A % * % e is equal to the trace of A. This is the basis for estimating the trace and the standard error for this estimate is based on NtrA independent samples.

dfind2d is a fast FORTRAN subroutine to find nearest neighbors within a fixed distance and is called by Wendland.basis. The function dfind3d is currently not used but is intended for future use to determine chordal distance between points on a sphere or cylinder.

LKrigDefaultFixedFunction Is called to construct the fixed part of the spatial model. The default is a polynomial of degree (m-1).

LKDefaultVarNames A handyfunction to create some simple column names to a matrix if they are missing. This is used so that the tables of parameter estimates will have labels.

Value

LKrig.coef

a list with components d.coef the coefficients of the spatial drift and for covariates (Z) and c.coef the basis function coefficients. The logical vector ind.drift from the LKrig object indicates with components of d.coef are associated with the polynomial spatial drift and which are other fixed spatial covariates.

LKrig.lnPlike

has the components:

lnProfileLike

the log likelihood profiled for lambda, alpha and a.wght

sigma2.MLE

the MLE of sigma2 given lambda, alpha and a.wght

shat.MLE

the MLE of tau given lambda, alpha and a.wght

quad.form

the quadratic form in the exponent of the multivariate normal likelihood

lnDetCov

the log determinant of the covariance matrix in the likelihood

LKrigDefaultFixedFunction

A matrix with dimension nrow(x) and columns of the number of polynomial terms and the number of columns of Z if given.

LKDefaultVarNames

The original column names ( e.g. the result of colnames(A) ) or some simple choices filled in.

Author(s)

Doug Nychka

References

Nychka, D., Bandyopadhyay, S., Hammerling, D., Lindgren, F., & Sain, S. (2015). A multi-resolution Gaussian process model for the analysis of large spatial datasets.Journal of Computational and Graphical Statistics, 24(2), 579-599.

See Also

LKrig, LKrig.basis

Description

Some utility functions used internally by higher level LKrig functions. Currently these are simple functions that perform shifts of a matrix and operations on indices for multidimensional arrays.

Usage

LKrig.shift.matrix( A, shift.row=0, shift.col=0, periodic=c(FALSE, FALSE))
LKrig.rowshift.periodic( A, shift.row)
LKrig.rowshift( A, shift.row, shift.col)
LKArrayShift(A, shift, periodic = FALSE)


expandMatrix0( A, B)
expandMatrix( ...)
expandMList( Mlist, byrow=TRUE)


convertIndexPeriodic(I, nGrid, nPad = NULL)
convertIndexArray(I, nGrid)
grid2Index(I, grid)

Arguments

Details

Shift related: These functions are used to create the nearest neighbor indices for the precision matrices.

Expand related: These functions are useful for creating a sets of covariance parameters that follow a factorial pattern. For example repeating the rows of the "alpha" parameters as the "a.wght" parameters are varied. expandMList is particularly useful for creating a factorial design of parameters to pass to LKrig.MLE for searching the likelihood.

Index related: The function convertIndexPeriodic converts a single index for a multidimensional array, into an index that reflects wrapping into a smaller grid to reflect the padding. This is used in the LKDistGrid function to handle distances when the grid has periodic dimensions. The other two functions are used for checking and debugging of going back and forth between the multidimensional and single indexes.

Value

Shift: A matrix of shifted values. Entries that are not defined due to the shift are set to NA. A column shift is done by a combination of transpose operations and a row shift.

A<- matrix( 1:12,3,4)
A
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12

#shift of 2 for rows:
 LKrig.rowshift( A, 2)
    [,1] [,2] [,3] [,4]
[1,]   NA   NA   NA   NA
[2,]   NA   NA   NA   NA
[3,]    1    4    7   10

#periodic case
LKrig.rowshift.periodic( A, 2)
     [,1] [,2] [,3] [,4]
[1,]    2    5    8   11
[2,]    3    6    9   12
[3,]    1    4    7   10

Expand: ExpandMList Returns a list of matrices where the original matrices are repeated so that combinations of rows are represented. The example below illustrates. byrow=FALSE does the repetition by columns instead of rows.

> A
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> B
     [,1]
[1,]   11
[2,]   12
[3,]   13
> C
[1,] 100
> expandMList( list( A=A, B=B, C=C))
$A
     [,1] [,2]
[1,]    1    3
[2,]    2    4
[3,]    1    3
[4,]    2    4
[5,]    1    3
[6,]    2    4

$B
     [,1]
[1,]   11
[2,]   11
[3,]   12
[4,]   12
[5,]   13
[6,]   13

$C
     [,1]
[1,]  100
[2,]  100
[3,]  100
[4,]  100
[5,]  100
[6,]  100

Author(s)

Doug Nychka

Examples

A<- array( 1:90, c( 4,5,3))
	LKArrayShift( A, c( -1,-1,0))	
	
# welcome to the world of unrolling multiarray indices
A<- array( 1:60, c( 4,3,5))	
I<- rbind( c(1,2,1), c( 3,2,5))
look<- grid2Index( I, c( 4,3,5) )
# A has been filled with the right unrolled index
print( look)
print(A[look])

Description

These functions support the LKrig function. Their main function is to create and evaluate radial basis functions of varying support on a nested set of regular grids. This series of grids forms a multi-resolution basis. The Gaussian process model is an expansion in these basis functions where the basis coefficients follow a Markov random field model for each resolution level. This family of functions generate the basis using sparse matrices, evaluate the covariance function of the process, and also simulate realizations of the process. LKrig.cov.plot is a useful function to get a quick plot of the covariance function implied by a LatticeKrig specification.

Usage

#
LKrig.cov(x1, x2 = NULL, LKinfo, C = NA, marginal = FALSE,
                 aRange = NA)
LKrig.cov.plot( LKinfo, NP=200, center = NULL, xlim = NULL, ylim = NULL)
LKrigCovWeightedObs(x1, wX, LKinfo) 

LKrig.basis(x1, LKinfo, Level = NULL, raw = FALSE, verbose = FALSE)
LKrig.precision(LKinfo, return.B = FALSE, verbose=FALSE)  
LKrig.quadraticform( Q, PHI, choleskyMemory = NULL) 
LKrig.spind2spam(obj, add.zero.rows=TRUE)
LKrigMarginalVariance(x1, LKinfo, verbose = FALSE)
LKFindSigma2VarianceWeights(x1, LKinfo)

Arguments

Details

The basis functions are two-dimensional radial basis functions based on the compactly supported stationary covariance function (Wendland covariance) and centered on regular grid points with the scaling tied to the grid spacing.

For a basis at the coarsest level, the grid centers are generated by expanding the two sequences

seq(grid.info$xmin,grid.info$xmax,grid.info$delta)
seq(grid.info$ymin,grid.info$ymax,grid.info$delta) 

into a regular grid of center points. The same spacing delta is used in both directions. The unnormalized basis functions are evaluated at locations x1 by finding the pairwise, radial distances among centers and x1, scaling by grid.info$delta * overlap and then evaluating with the function name passed as BasisFunction. By default this is the 2-d Wendland covariance of order 2. Perhaps the most important point about the LKrig.basis is that it is designed to return a matrix of all basis functions as a sequence of points. There is no need to have a function that evaluates basis functions individually. In R code for a set of locations x1 and a rectangular spatial domain with ranges xmin, xmax, ymin ,ymax:

 centers<- expand.grid(seq(xmin,xmax,delta),
                       seq(ymin,ymax,delta) )
 bigD<- rdist( x1, centers)/(delta*2.5)
 PHI<- Wendland.function( bigD)

Note that there will be nrow(centers) basis functions generated where the precise number depends on the range of the domain and the choice of delta. The basis functions are indexed by the columns in PHI and this is a convention throughout this package. There will be nrow(x1) rows in PHI as each basis function will be evaluated at each 2-d location.

The basis functions are then normalized by scaling the basis functions at each location so that resulting marginal variance of the process is 1. This is done to coax the covariance model closer to a stationary representation. It is easiest to express this normalization by pseudo R code:

If Q is the precision matrix of the basis coefficients then in R/LatticeKrig code:

Omega<-  solve(Q)
process.variance <- diag(PHI%*% Omega %*%t(PHI) )
PHI.normalized <-  diag(1/sqrt(process.variance)) %*% PHI

where Omega is the unnormalized covariance matrix of the basis function coefficients.

Although correct, the code above is not an efficient algorithm to compute the unnormalized process variance. First the normalization can be done level by level rather than dealing with the entire multi-resolution process at once. Also it is important to work with the precision matrix rather than the covariance. The function LKrigNormalizeBasis takes advantage of the sparsity of the precision matrix for the coefficients. LKrigNormalizeBasisFast is a more efficient version for an isotropic type model when a.wght is constant for a given level and takes advantage of the Kronecker structure on the precision matrix at each level. LKrigNormalizeBasisFFTInterpolate is the fastest method, which provides an accurate but approximate method for the normalization, which takes advantage of situations where the number of locations significantly exceeds the number of basis functions.

The precision matrix for the basis coefficients at each resolution has the form t(B)%*% B. These matrices for the individual levels are assembled by LKrig.precision as the block diagonals of a larger precision matrix for the entire vector of coefficients. Note these matrices are created in a sparse format. The specific entries in B, the object created by LKrig.MRF.precision, are a first order Markov random field: without edge adjustments the diagonal elements have the value a.wght and the first order neighbors have the value -1.

Below we give more details on how the weights are determined. Following the notation in Lindgren and Rue a.wght= 4 + k2 with k2 greater than or equal to 0. Some schematics for filling in the B matrix are given below (values are weights for the SAR on the lattice with a period indicating zero weights).

                                                 __________ 
   .   -1     .         |  -e      .            |  a.wght  -e
                        |                       |
  -1  a.wght  -1        | a.wght   -e           |   -e      .
                        |                       |
   .  -1      .         |  -e      .            |    .      .

Interior point        Left edge              Upper left corner

To adjust for edges and corner lattice points we take two strategies. The first is add extra lattice points around the edges so the actual spatial domain has lattice points with the complete number of neighbors. The default for the LKRectangle geometry is to add a buffer of 5 points ( NC.buffer = 5 ) around the edges. In addition the neighbor weights are inflated sum to a fixed value. For example in the case of LKRectangle and following the schematic above we want the neighbor weights to sum to -4. Thus "e" for the middle figure will be set to -4/3 and for the right figure -2. Empirically these adjustments were found to improve how a.wght parameters chosen close to 4 could generate long range spatial correlations in the lattice coefficients. Note that these simple boundary adjustments happen to the buffer points and so are less likely to introduce artifacts into the spatial domain.

Value

LKrig.basis: A matrix with number of rows equal to the rows of x1 and columns equal to the number of basis functions (LKinfo$m). Attached to the matrix is an info attribute that contains the list passed as LKinfo. Usually this value is in spam sparse matrix format.

LKrig.precision: For return.B ==FALSE a sparse, square matrix with dimensions of the number of basis functions. For return.B == TRUE the "B" SAR matrix is returned. This is useful for checking this function.

LKrig.cov: If C=NA a cross covariance matrix with dimensions nrow(x1) and nrow(x2) is used. If C is passed the result of multiplying the cross covariance matrix times C is used.

LKrigMarginalVariance: Gives the marginal variance of the LatticeKrig process at each level and at the locations in x1. Returned value is a matrix with nlevel columns indexing the levels and the number of rows equal to nrow(x1). If varLevels is a row of this matrix then sum( varLevels* LKinfo$alpha) is the marginal variance for the full process when the different levels are weighted by alpha. This is weighted sum is the marginal variance returned by LKrig.cov and marginal=TRUE ( Also assuming that LKinfo$sigma2.object is NULL, which it usually is.) .

LKFindSigma2VarianceWeights Return either the marginal variance of the process as a scalar or the variance at arbitrary locations x1. This second role is part of a nonstationary specification of the model where the process marginal variance can vary over space.

LKrig.sim: A matrix with dimensions of nrow(x1) by M. Each column are vectors of simulated values at the locations x1.

LKrig.cov.plot: Evaluates the covariance specified in the list LKinfo with respect to the point center along a transects in the x and y directions intersecting this point. Note the rectangular extent of the spatial domain is part of the grid information in LKinfo. Returns components u, d and cov. Each of these are two column matrices with the first column being results in the x direction and second column in the y direction. d are the distances of the points from the center and u are the actual x or y spatial coordinates. cov are the values of the covariance function. If normalize is TRUE these will in fact be the correlation functions. To plot the returned list use

 out<- LKrig.cov.plot(LKinfo)
 matplot( out$d, out$cov, type="l")

LKrig.quadraticform: Returns a vector that is diag(t(PHI)%*% solve(Q) %*% PHI)) closely related to the marginal variance of the process.

LKrigNormalizeBasis, LKrigNormalizeBasisFast, LKrigNormalizeBasisFFTInterpolate: A vector of variances corresponding to the unnormalized process at the locations.

LKrig.spind2spam: This converts a matrix in spind sparse format to spam sparse format. Although useful for checking, LatticeKrig now uses the the spam function directly to do the conversion. If obj is a list in spind format then the

obj2<- LKrig.spind2spam(obj)

and

obj2<- spam(obj[c("ind","ra")], nrow=obj$da[1], ncol=obj$da[2])

give the same result.

Author(s)

Doug Nychka

See Also

LKrig, mKrig, Krig, fastTps, Wendland, LKrigSAR, Radial.basis

Examples

# Load ozone data set
  data(ozone2)  
  x<-ozone2$lon.lat
  y<- ozone2$y[16,]
# Find location that are not 'NA'.
# (LKrig is not set up to handle missing observations.)
  good <-  !is.na( y)
  x<- x[good,]
  y<- y[good]
  LKinfo<- LKrigSetup( x,NC=20,nlevel=1, alpha=1, lambda= .3 , a.wght=5)
# BTW lambda is close to MLE 

# What does the  LatticeKrig covariance function look like?
# set up LKinfo object
# NC=10 sets the grid for the first level of basis functions
# NC^2 = 100 grid points in first level if square domain.
# given four levels the number of basis functions
# = 10^2 + 19^2 +37^2 + 73^2 = 5329
# effective range scales as roughly kappa where a.wght =  4 + kappa^2    
# or exponential decreasing marginal variances for the components.
    NC<- 10 
    nlevel <- 4
    a.wght <-  4 + 1/(.5)^2
    alpha<-  1/2^(0:(nlevel-1)) 
    LKinfo2<- LKrigSetup( cbind( c( -1,1), c(-1,1)), NC=NC,
                   nlevel=nlevel, a.wght=a.wght,alpha=alpha)
# evaluate covariance  along the  horizontal line through
# midpoint of region -- (0,0) in this case. 
    look<- LKrig.cov.plot( LKinfo2)
# a plot of the covariance function in x and y with respect to (0,0)
    set.panel(2,1)  
    plot(look$u[,1], look$cov[,1], type="l")
    title("X transect")
    plot(look$u[,2], look$cov[,2], type="l")
    title("Y transect")
    set.panel(1,1)
#
#
## Not run: 
# full 2-d view of the covariance (this example follows the code
# in LKrig.cov.plot)
 x2<- cbind( 0,0)
 x1<- make.surface.grid( list(x=seq( -1,1,,40),  y=seq( -1,1,,40)))
 look<- LKrig.cov( x1,x2, LKinfo2)
 contour( as.surface( x1, look))
# Note nearly circular contours.
# of  course  plot(look[,80/2]) should look like plot above.
#

## End(Not run)

## Not run: 
#Some correlation functions from different models
set.panel(2,1)
# a selection of ranges:
  hold<- matrix( NA, nrow=150, ncol=4)
  kappa<- seq( .25,1,,4)
  x2<- cbind( 0,0)
  x1<-  cbind( seq(-1,1,,150), rep( 0,150))
  for( k in 1:4){
    LKtemp<-  LKrigSetup( cbind( c( -1,1), c(-1,1)), NC=NC,
                   nlevel=nlevel,
                   a.wght= 4  + 1/(kappa[k]^2),
                   alpha=alpha)
    hold[,k]<-  LKrig.cov( x1,x2, LKinfo=LKtemp)
  }
  matplot( x1[,1], hold, type="l", lty=1, col=rainbow(5), pch=16 )
# a selection of smoothness parameters
  ktemp<- .5 # fix range
  alpha.power<- seq( 1,4,,4)
  LKtemp<- LKinfo2
  for( k in 1:4){
   LKtemp<-  LKrigSetup( cbind( c( -1,1), c(-1,1)), NC=NC,
                   nlevel=nlevel,
                   a.wght= 4  + 1/(ktemp^2),
                   alpha=alpha^alpha.power[k])
    hold[,k]<-  LKrig.cov( x1,x2, LKinfo=LKtemp)
  }
  matplot( x1[,1], hold, type="l", lty=1, col=rainbow(5) )
 set.panel()

## End(Not run)
 
## Not run: 
# generating a basis on the domain [-1,1] by [-1,1] with 1 level
# Default number of buffer points are added to each side. 
  LKinfo<- LKrigSetup(cbind( c(-1,1), c(-1,1)), NC=6,
                                 nlevel=1, a.wght=4.5,alpha=1, NC.buffer=0 )
# evaluate the basis functions on a grid to look at them
  xg<- make.surface.grid( list(x=seq(-1,1,,50), y= seq(-1,1,,50)))
  PHI<- LKrig.basis( xg,LKinfo)
  dim(PHI) # should be  2500=50^2  by  36=6^2
# plot the 9th basis function  as.surface is a handy function to
# reformat the vector as an image object
# using the grid information in an attribute of the grid points
  set.panel(1,3) 
  image.plot(as.surface(xg, PHI[,9]))
  points(  make.surface.grid( LKrigLatticeCenters(LKinfo, 1)) , col="grey", cex=.5)
  title("A radial basis function")
# compare to the tensor product basis type
  LKinfo2<- LKrigSetup(cbind( c(-1,1), c(-1,1)), NC=6,
                                 nlevel=1, a.wght=4.5,alpha=1, NC.buffer=0,
                                 BasisType="Tensor" )
  PHI2<- LKrig.basis( xg,LKinfo2)
  image.plot(as.surface(xg, PHI2[,9]))
  points(  make.surface.grid( LKrigLatticeCenters(LKinfo, 1)), col="grey", cex=.5)
  title("Tensor product basis function")
  
  image.plot(as.surface(xg, PHI[,9] - PHI2[,9]))
  points(  make.surface.grid( LKrigLatticeCenters(LKinfo, 1)), col="grey", cex=.5)
  title(" Radial - Tensor for 9th basis function")                       
set.panel()

## End(Not run)
#
# example of basis function indexing
#
## Not run: 
# generating a basis on the domain [-1,1]X[-1,1] with 3 levels
# note that there are no buffering grid points.
  set.panel(3,2)
  LKinfo<-LKrigSetup(cbind( c(-1,1), c(-1,1)), NC=6,
                    a.wght=5, alpha=c(1,.5,.25), nlevel=3,
                    NC.buffer=0)
# evaluate the basis functions on a grid to look at them
  xtemp<- seq(-1,1,,40)
  xg<- make.surface.grid( list(x=xtemp, y= xtemp) )
  PHI<- LKrig.basis( xg,LKinfo)
# coerce to dense matrix format to make plotting easier.
  PHI<- spam2full(PHI)
# first tenth, and last basis function in each resolution level
# basis functions centers are added
 set.panel(3,3)
    for(  j in 1:3){
      id1<- LKinfo$latticeInfo$offset[j]+ 1
      id2<-  LKinfo$latticeInfo$offset[j]+ 10
      idlast<- LKinfo$latticeInfo$offset[j] +
                  LKinfo$latticeInfo$mx[j,1]*LKinfo$latticeInfo$mx[j,2]
   
      centers<-  make.surface.grid(LKrigLatticeCenters(LKinfo, j) )
      image.plot( as.surface(xg, PHI[,id1]))
      points( centers, cex=.2, col="grey")
      image.plot(as.surface(xg, PHI[,id2]))
      points( centers, cex=.2, col="grey")
      image.plot( as.surface(xg, PHI[,idlast]))
      points( centers, cex=.2, col="grey")
}
  set.panel()

## End(Not run)
## Not run: 
# examining the stationarity of covariance model
  temp.fun<- 
     function( NC.buffer=0, NC=4,  a.wght=4.01){
        LKinfo<- LKrigSetup(cbind( c(-1,1), c(-1,1)),nlevel=1, alpha=1,
                                 a.wght=a.wght, NC=NC,   
                                 NC.buffer=NC.buffer,
                                  choleskyMemory=list(nnzR=2e6))
        cov1y<- cov1x<- cov0x<- cov0y<-  matrix( NA, nrow=200, ncol=20)
        cov1dx<- cov1dy<- cov0dx<- cov0dy<- matrix( NA, nrow=200, ncol=20)
        cgrid<- seq( 0,1,,20)
        for( k in 1:20){
            hold<- LKrig.cov.plot( LKinfo,
                            center=rbind( c(cgrid[k], cgrid[k])), NP=200)
            cov1x[,k] <- hold$cov[,1]
            cov1y[,k] <- hold$cov[,2]
            cov1dx[,k] <- hold$d[,1]
            cov1dy[,k] <- hold$d[,2]
            hold<- LKrig.cov.plot( LKinfo,
                             center=rbind( c(cgrid[k],0) ), NP=200)
            cov0x[,k] <- hold$cov[,1]
            cov0y[,k] <- hold$cov[,2]
            cov0dx[,k] <- hold$d[,1]
            cov0dy[,k] <- hold$d[,2]
                }
         matplot( cov1dx, cov1x, type="l", col= rainbow(20),
                         xlab="", ylab="correlation")
         mtext( side=1, line=-1, text="diagonal X")
         title( paste(  " buffer=",NC.buffer), cex=.5)
         matplot( cov1dy, cov1y, type="l", col= rainbow(20),
                        xlab="", ylab="")
         mtext( side=1, line=-1, text="diagonal Y")
         matplot(cov0dx, cov0x, type="l", col= rainbow(20),
                        xlab="",       ylab="")
         mtext( side=1, line=-1, text="middle X")
         matplot( cov0dy, cov0y, type="l", col= rainbow(20),
                         xlab="",   ylab="")
         mtext( side=1, line=-1, text="middle Y")
         title( paste( NC, a.wght), cex=.5)
}


 set.panel(3,4)
par(mar=c(3,4,1,0), oma=c(1,1,1,1))
temp.fun(  NC.buffer=5, NC=4, a.wght=4.05)
temp.fun(  NC.buffer=5, NC=16, a.wght=4.05)
temp.fun(  NC.buffer=5, NC=64, a.wght=4.05)

set.panel(4,4)
par(mar=c(3,4,1,0), oma=c(1,1,1,1))
temp.fun( NC.buffer=0, NC=8)
temp.fun( NC.buffer=2, NC=8)
temp.fun( NC.buffer=4, NC=8)
# this next one takes a while
temp.fun( NC.buffer=8,  NC=8)
# stationary == curves should all line up!


## End(Not run)

Description

Given a list of different covariance parameters for the Lattice Krig covariance model this function computes the likelihood or a profiled version (over lambda) and approximates a generalized cross-validation function at each of the parameter settings. This is an experimental function that has been productively used with a Latin hypercube design package to efficiently search through the LatticeKrig covariance parameter space.

Usage

LKrigFindLambda(x, y, Z = NULL, U = NULL, X = NULL, ..., LKinfo,
                 use.cholesky = NULL, lambda.profile = TRUE,
                 lowerBoundLogLambda = -16, tol = 0.005, verbose =
                 FALSE)
			    
LKrigFindLambdaAwght(x, y, ..., LKinfo, use.cholesky = NULL,
                 lowerBoundLogLambda = -16, upperBoundLogLambda = 4,
                 lowerBoundOmega = -3, upperBoundOmega = 0.75, factr =
                 1e+07, pgtol=1e-1, maxit = 15, verbose = FALSE)

                                 
LambdaAwghtObjectiveFunction(PARS, LKrigArgs, capture.env, verbose=FALSE )

LKrig.MLE( x,y,..., LKinfo, use.cholesky = NULL,
                    par.grid=NULL,
                    lambda.profile=TRUE,
                    verbose=FALSE,
                    lowerBoundLogLambda = -16,
                    nTasks = 1, taskID = 1,
                    tol = 0.005)
LKrig.make.par.grid(par.grid=NULL, LKinfo = NULL) 
omega2Awght (omega, LKinfo)
Awght2Omega (Awght, LKinfo)

Arguments

par.grid$llambda[k]
par.grid$alpha[k,]
par.grid$a.wght[k,]

Details

LKrigFindLambda: Uses a simple one dimensional optimizer optimize. To maximize the log likelihood for log lambda over the range: llambda.start + [-8,5]. This function is used to determine lambda in LatticeKrig.

LKrigFindLambdaAwght: Uses a simple optimizer optim. To maximize the log likelihood for lambda and a.wght over the range.

LKrig.MLE: This is a simple wrapper function to accomplish repeated calls to the LKrig function to evaluate the profile likelihood and/or to optimize the likelihood over the lambda parameters. The main point is that maximization over the lambda parameter (or equivalently for tau and sigma2) is the most important and should be done before considering variation of other parameters. If lambda is specified then one has closed form expressions for tau, sigma2 that can then be substituted back into the log full likelihood. This operation that is the default throughout LatticeKrig (and fields) can concentrate the likelihood on a reduced set of parameters. The further refinement when lambda.profile==TRUE is to maximize the concentrated likelihood over lambda and report this result. This will be a profile likelihood over the remaining parameters of the covariance.

The covariance/model parameters are alpha, a.wght, and log lambda and are separate matrix or vector components of the par.grid list. The cleanest version of this function would just require the par.grid list, however, to be easier to use there are several options to give partial information and let the function itself create the master parameter list. For example, just a search over lambda should be easy and not require creating par.grid outside the function. To follow this option one can just give an LKinfo object. The value for the lambda component in this object will be the starting value with the default starting value being lambda =1.0.

In the second example below most of the coding is getting the grid of parameters to search in the right form. It is useful to normalize the alpha parameters to sum to one so that the marginal variance of the process is only parameterized by sigma2. To make this easy to implement there is the option to specify the alpha parameters in the form of a mixture model so that the components are positive and add to one (the gamma variable below). If a component gamma is passed as a component of par.grid then this is assumed to be in the mixture model form and the alpha weights are computed from this. Note that gamma will be a matrix with (nlevel - 1) columns while alpha has nlevel columns.

For those readers that use which.max these functions are natural extensions and are handy for looking at interpolated surfaces of the likelihood function.

which.max.matrix: Finds the maximum value in a matrix and returns the row/column index.

which.max.image Finds the maximum value in an image matrix and returns the index and the corresponding grid values.

LKrig.make.par.grid: This is usually used as an internal function that converts the list of parameters in par.grid and the LKinfo object into an more complex data structure used by LKrig.MLE. Its returned value is a "list of lists" to make the search over different parameters combinations simple.

omega2Awght, Awght2omega Converts between the Awght parameter (the diagonal elements of the SAR matrix) and the omega parameter that provides an unconstrained range for optimization. The link is

Awght <- LKinfo$floorAwght + exp(omega) * (xDimension)

.

Note that LKinfo supplied the lower bound on the Awght because this is geometry/problem specific. For the 2-d rectangle this is 4.

Value

LKrigFindLambda

LKrig.MLE

which.max.matrix Returns a 2 column matrix with row and column index of maximum.

which.max.image For an object in image format returns components x,y,z giving the location and the maximum value for the image. Also included is the component ind that is the row and column indices for the maximum in the image matrix.

LKrig.make.par.grid Returns a list with components, alpha, a.wght. Each component is a list where each component of the list is a separate set of parameters. This more general format is useful for the non-stationary case when the parameters alpha might be a list of nlevel matrices.

Author(s)

Douglas Nychka

See Also

LKrig LatticeKrig

Examples

# 
# fitting summer precip for  sub region of North America (Florida)
# (tiny subregion is just to make this run under 5 seconds). 
# total precip in 1/10 mm for JJA 
  data(NorthAmericanRainfall)
# rename for less typing
  x<- cbind( NorthAmericanRainfall$longitude, NorthAmericanRainfall$latitude)
  y<- log10(NorthAmericanRainfall$precip)
# cut down the size of this data set so examples run quickly
  ind<- x[,1] > -90 & x[,2] < 35 #
  x<- x[ind,]
  y<- y[ind]

# This is a single level smoother
 
  LKinfo<- LKrigSetup(x,NC=4, nlevel=1, a.wght=5, alpha=1.0)
  lambdaFit<- LKrigFindLambda( x,y,LKinfo=LKinfo)
  lambdaFit$summary

## Not run: 
# grid search over parameters 
  NG<-15
  par.grid<- list( a.wght= rep( 4.05,NG),alpha= rep(1, NG),
                      llambda=  seq(-8,-2,,NG))
  lambda.search.results<-LKrig.MLE( x,y,LKinfo=LKinfo,
                                    par.grid=par.grid,
                                    lambda.profile=FALSE)
  lambda.search.results$summary
# profile likelihood
  plot( lambda.search.results$summary[,1:2], 
         xlab="effective degrees freedom",
         ylab="ln profile likelihood")
# fit at largest likelihood value:
  lambda.MLE.fit<- LKrig( x,y,
                    LKinfo=lambda.search.results$LKinfo.MLE)

## End(Not run)                    
                    
## Not run:                     
# optimizing  Profile likelihood over lambda using optim
# consider 3 values for a.wght (range parameter)
# in this case the log lambdas passed are the starting values for optim.
  NG<-3
  par.grid<- list( a.wght= c( 4.05,4.1,5) ,alpha= rep(1, NG),
                      llambda= c(-5,NA,NA))
# NOTE: NAs in llambda mean use the previous MLE for llambda as the
# current starting value. 
  LKinfo<- LKrigSetup(x,NC=12,nlevel=1, a.wght=5, alpha=1.0) 
  lambda.search.results<-LKrig.MLE(
                              x,y,LKinfo=LKinfo, par.grid=par.grid,
                              lambda.profile=TRUE)
  print(lambda.search.results$summary)
# note first result a.wght = 4.05 is the optimized result for the grid
# search given above.

## End(Not run)
########################################################################    
# search over two multi-resolution levels varying the  levels of alpha's
########################################################################
## Not run: 
# NOTE: search ranges found largely by trial and error to make this
# example work also the grid is quite coarse ( and NC is small) to
# be quick as a help file example
  data(NorthAmericanRainfall)
# rename for less typing
  x<- cbind( NorthAmericanRainfall$longitude, NorthAmericanRainfall$latitude)
# total precip in 1/10 mm for JJA 
 y<- log10(NorthAmericanRainfall$precip)
# cut down the size of this data set so examples run quickly
# examples also work with  the full data set. Also try NC= 100 for a
# nontrivial model.
  ind<- x[,1] > -90 & x[,2] < 35 #
  x<- x[ind,]
  y<- y[ind]
  
  Ndes<- 10  
# NOTE: this is set to be very small just to make this
#       example run fast
  set.seed(124)
  par.grid<- list()
# create grid of alphas to sum to 1 use a mixture model parameterization
#  alpha1 = (1/(1 + exp(gamma1)) ,
# alpha2 = exp( gamma1) / ( 1 + exp( gamma1))
# 
  par.grid$gamma<- cbind(runif( Ndes, -3,2), runif( Ndes, -3,2))
  par.grid$a.wght<- rep( 4.5, Ndes)
# log lambda grid search values
  par.grid$llambda<- runif( Ndes,-5,-3)  
  LKinfo1<- LKrigSetup( x, NC=5, nlevel=3, a.wght=5, alpha=c(1.0,.5,.25))
# NOTE: a.wght in call is not used. Also a better search is to profile over
#  llambda

 alpha.search.results<- LKrig.MLE( x,y,LKinfo=LKinfo1, par.grid=par.grid,
                                    lambda.profile=FALSE)

########################################################################
# Viewing the search results
########################################################################

# this scatterplot is good for a quick look because  effective degrees
# of freedom is a useful summary of fit. 
  plot( alpha.search.results$summary[,1:2], 
         xlab="effective degrees freedom",
         ylab="ln profile likelihood")
#

## End(Not run)

## Not run: 
# a two level model search 
# with profiling over lambda.
data(NorthAmericanRainfall)
# rename for less typing
  x<- cbind( NorthAmericanRainfall$longitude,
             NorthAmericanRainfall$latitude)
# mean total precip in 1/10 mm for JJA 
  y<- log10(NorthAmericanRainfall$precip)

#  This takes a few minutes
  Ndes<- 40 
  nlevel<-2 
  par.grid<- list()
## create grid of alphas to sum to 1 use a mixture model parameterization:
#    alpha1 = (1/(1 + exp(gamma1)) ,
#   alpha2 = exp( gamma1) / ( 1 + exp( gamma1))
  set.seed(123)
  par.grid$gamma<- runif( Ndes,-3,4)
## values for range (a.wght)
  par.grid$a.wght<- 4 + 1/exp(seq( 0,4,,Ndes))
# log lambda grid search values (these are the starting values)
  par.grid$llambda<- rep(-4, Ndes)

  LKinfo1<- LKrigSetup( x, NC=15, nlevel=nlevel, 
                          a.wght=5, alpha=rep( NA,2) ) 
##
## the search over the parameter list in par.grid  maximizing over lambda 
  search.results<- LKrig.MLE( x,y,LKinfo=LKinfo1, par.grid=par.grid,
                                 lambda.profile=TRUE)
# plotting results of likelihood search
set.panel(1,2)
 plot( search.results$summary[,1:2], 
         xlab="effective degrees freedom",
         ylab="ln profile likelihood")
 xtemp<- matrix(NA, ncol=2, nrow=Ndes)
 for( k in 1:Ndes){
   xtemp[k,] <- c( (search.results$par.grid$alpha[[k]])[1],
                  (search.results$par.grid$a.wght[[k]])[1] )
}
 quilt.plot( xtemp,search.results$summary[,2])
#  fit using Tps
 tps.out<- Tps(  xtemp,search.results$summary[,2], lambda=0)
 contour( predictSurface(tps.out), lwd=3,add=TRUE)
 set.panel()

## End(Not run)
## Not run: 
# searching over nu 
data(ozone2)
x<- ozone2$lon.lat
y<- ozone2$y[16,]
good<- !is.na(y)
y<- y[good]
x<- x[good,]
par.grid<- expand.grid( nu = c(.5,1.0, 1.5), a.wght= list(4.1,4.5,5) )
par.grid$llambda<- rep( NA, length(par.grid$nu))
LKinfo<- LKrigSetup(x,  nlevel=3, nu=.5, NC=5, a.wght=4.5)
out<- LKrig.MLE( x,y, LKinfo=LKinfo, par.grid=par.grid)
# take a look
cbind( par.grid, out$summary[,1:5])

## End(Not run)
## Not run: 
# an MLE fit taking advantage of replicated fields
# check based on simulated data

N<- 200
M<-50 # number of independent replicated fields
tau<- sqrt(.01)
set.seed(123)
x<- matrix( runif(N*2), N,2)
                
LKinfo<- LKrigSetup( x, NC=16, nlevel=1,
                 a.wght=4.5, lambda=.01, 
                fixed.Function=NULL,
                normalize=TRUE)  
                
# the replicate fields
truef<-  LKrig.sim(x,LKinfo=LKinfo, M=M )
set.seed(222)
error<- tau*matrix( rnorm(N*M), N,M)
y<- truef + error 
# with correct lambda
obj<- LKrig( x,y, LKinfo=LKinfo, lambda=.01, )
print( obj$tau.MLE.FULL)
print( obj$sigma2.MLE.FULL)

fitMLE1<- LKrigFindLambda( x,y, LKinfo=LKinfo)
fitMLE1$summary
aWghtGrid<-  c( 4.01, 4.05, 4.1, 4.2, 4.5, 4.6, 4.7, 5, 8, 16)
par.grid<- list( a.wght = aWghtGrid)

fitMLE2<- LKrig.MLE( x,y, LKinfo=LKinfo,
                      par.grid= par.grid )
fitMLE2$summary   

LKinfo1<- LKinfoUpdate( LKinfo, lambda=.1, a.wght= 4.2)                   
fitMLE4<- LKrigFindLambdaAwght( x,y, LKinfo=LKinfo1)
fitMLE4$summary

plot(  log( aWghtGrid -4)/2, fitMLE2$summary[,2], type="b",
  xlab="log( a.wght - 4)/2",
  ylab= "log Profile likelihood" )


points( log(fitMLE4$summary["a.wght.MLE"] -4)/2,
     fitMLE4$summary["lnProfLike"], pch="+", col="red"  )
xline( log(fitMLE4$summary["a.wght.MLE"] -4)/2, col="red", lty=2)

## End(Not run)

Description

The fields are Gaussian and can be either simulated unconditionally or conditional on the field values and a set of irregular locations.

Usage

#
LKrig.sim(x1, LKinfo, M=1,just.coefficients = FALSE)
LKrig.sim.conditional( LKrigObj, M=1, x.grid= NULL, grid.list=NULL,
                           nx=80, ny=80,...,Z.grid=NULL, seed=42, verbose=FALSE)

simConditionalDraw(index = 1, LKrigObj, ghat, x.grid, Z.grid, PHIGrid,
                 seeds = 123, verbose = FALSE)

Arguments

Details

The simulation of the unconditional random field is done by generating a draw from the multi-resolution coefficients using a Cholesky decomposition and then multiplying by the basis functions to evaluate the field at arbitrary points. Currently, there is no provision to exploit the case when one wants to simulate the field on a regular grid. The conditional distribution is a draw from the multivariate normal for the random fields conditioned on the observations and also conditioned on covariance model and covariance parameters. If the nugget/measurement error variance is zero then any draw from the conditional distribution will be equal to the observations at the observation locations. In the past conditional simulation was known to be notoriously compute intensive, but the major numerical problems are finessed here by exploiting sparsity of the coefficient precision matrix.

The conditional field is found using a simple trick based on the linear statistics for the multivariate normal. One generates an unconditional field that includes the field values at the observations. From this realization one forms a synthetic data set and uses LKrig to predict the remaining field based on the synthetic observations. The difference between the predicted field and the realization (i.e. the true field) is a draw from the conditional distribution with the right covariance matrix. Adding the conditional mean to this result one obtains a draw from the full conditional distribution. This algorithm can also be interpreted as a variant on the bootstrap to determine the estimator uncertainty. The fixed part of the model is also handled correctly in this algorithm. See the commented source for LKrig.sim.conditional for the details of this algorithm.

simConditionalDraw is low level function that is called to generate each ensemble member i.e. each draw from the conditional distribution. The large number of arguments is to avoid recomputing many common elements during the loop in generating these draws. In particular passing the basis function matrices avoid having to recompute the normalization at each step, often an intensive computation for a large grid.

Value

LKrig.sim: A matrix with dimensions of nrow(x1) by M of simulated values at the locations x1.

LKrig.sim.conditional: A list with the components.

xgrid

The locations where the simulated field(s) are evaluated.

ghat

The conditional mean at the xgrid locations.

g.draw

A matrix with dimensions of nrow(x.grid) by M with each column being an independent draw from the conditional distribution.

Author(s)

Doug Nychka

See Also

LKrig, mKrig, Krig, fastTps, Wendland

Examples

# Load ozone data set
  data(ozone2)  
  x<-ozone2$lon.lat
  y<- ozone2$y[16,]
# Find location that are not 'NA'.
# (LKrig is not set up to handle missing observations.)
  good <-  !is.na( y)
  x<- x[good,]
  y<- y[good]
  LKinfo<- LKrigSetup( x,NC=20,nlevel=1, alpha=1, lambda= .3 , a.wght=5)
# BTW lambda is close to MLE 
# Simulating this  LKrig process
# simulate 4 realizations of process and plot them
# (these have unit marginal variance)
  xg<- make.surface.grid(list( x=seq( -87,-83,,40), y=seq(36.5, 44.5,,40)))
  out<- LKrig.sim(xg, LKinfo,M=4)
## Not run: 
  set.panel(2,2)
  for( k in 1:4){
    image.plot( as.surface( xg, out[,k]), axes=FALSE) }

## End(Not run)
  obj<- LKrig(x,y,LKinfo=LKinfo)
  O3.cond.sim<- LKrig.sim.conditional( obj, M=3,nx=40,ny=40) 
## Not run: 
  set.panel( 2,2)
  zr<- range( c(  O3.cond.sim$draw,  O3.cond.sim$ghat), na.rm=TRUE)
  coltab<- tim.colors()
  image.plot( as.surface( O3.cond.sim$x.grid, O3.cond.sim$ghat), zlim=zr)
  title("Conditional mean")
  US( add=TRUE)
  for( k in 1:3){
    image( as.surface( O3.cond.sim$x.grid, O3.cond.sim$g.draw[,k]),
              zlim=zr, col=coltab)
    points( obj$x, cex=.5)
    US( add=TRUE)
  }
  set.panel()

## End(Not run)

Description

Creates matrix of low order polynomial in the spatial coordinates and adds any other spatial covariates that are part of the linear model.

Usage

LKrigDefaultFixedFunction(x, Z = NULL, m=2,distance.type="Euclidean") 
 LKrigPeriodicFixedFunction(x, Z = NULL, m = 2, distance.type = "Euclidean")
 
predictLKrigFixedFunction(object, xnew, Znew = NULL, drop.Z = FALSE,
collapseFixedEffect = FALSE)

Arguments

Details

LKrigDefaultFixedFunction This function creates the regression matrix for the fixed part of the spatial model. The default is a low order polynomial regression matrix of degree m-1. To this matrix are bound as columns any covariates passed as Z. Typically one would not need to modify this function. For more exotic fixed part models one can specify create and then specify a different function. See LKrig.setup and LKrig. NOTE: If the argument for this function is passed as NULL then the subsequent computations do not include a fixed part in the model.

LKrigDefaultFixedFunction This is same as LKrigDefaultFixedFunction except the first coordinate is ignored. i.e. it is assumed to be periodic so adding a polynomial does not make sense.

predictLKrigFixedFunction This function is simple, but is introduced to make the code modular and to handle the case for cylindrical geometry where only latitude should have a spatial term (to preserve periodicity in longitude).

Value

A matrix where rows index the locations and columns are the different spatial polynomial and covariates.

Author(s)

Doug Nychka

See Also

LKrig.basis, LKrig

Examples

x<- matrix( runif(100), nrow=50)
# linear polynomial 
T.matrix<- LKrigDefaultFixedFunction(x, m=2)
# quadratic polynomial 
T.matrix<- LKrigDefaultFixedFunction(x, m=3)