R: Spatial prediction and inference using a compactly supported...
LKrig
R Documentation
Spatial prediction and inference using a compactly supported
multi-resolution basis and a lattice model for the basis coefficients.
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, ..., verbose = FALSE)
## S3 method for class 'LKrig'
predict(object, xnew = object$x, Znew = NULL, drop.Z = FALSE,
just.fixed = FALSE, return.levels = FALSE, ...)
## 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, ...)
Arguments
x
Spatial locations of observations. Or the LKrig
object for printing.
y
Spatial observations.
weights
A vector that is proportional to the reciprocal
variances of the errors. I.e. errors are assumed to be uncorrelated
with variances sigma^2/weights.
Z
Linear covariates to be included in fixed part of the model
that are distinct from the default first order polynomial in
x (i.e. the spatial drift).
chull.mask
An additional constraint for evaluating predictions see
help(in.poly). Usually this is not needed.
digits
Number of digits in printed output.
drop.Z
If true the fixed part will only be evaluated at the
spatial drift polynomial part of the fixed model. The contribution
from the other, Z, covariates in the fixed component will be omitted.
extrap
If TRUE will extrapolate predictions beyond convex hull of locations.
grid.list
A list with components x and y specifying the grid ( see help( grid.list) ).
iseed
Random seed used in the Monte Carlo technique for
approximating the effective degrees of freedom (trace of the
smoothing matrix). If NA, no seed is set.
just.fixed
If TRUE just the fixed part of the model is evaluated.
LKinfo
An object whose components specify the LatticeKrig
spatial model. This is usually created by the function
LKrigSetup. If NULL, this object is created and returned as
a component of the LKrig object.
nx
Number of grids in x for prediction grid.
ny
Number of grids in y for prediction grid.
NtrA
Number of random samples used in Monte Carlo method for
determining effective degrees of freedom.
object
The LKrig object.
return.cholesky
If TRUE the Cholesky decomposition is included
in the output list (with the name Mc). This is needed for
some of the subsequent computations such as finding prediction
standard errors. Set this argument to FALSE to avoid much larger
objects when the decomposition is not needed. This option is often
paired with a subsequent call to LKrig with
use.cholesky. See the MLE.LKrig source code for details.
return.levels
If TRUE the predicted values for each level are
returned as columns in a matrix.
return.wXandwU
If TRUE the matrices wX and xU are included in the LKrig object.
U
For linear inverse problems the matrix that maps coefficients of the fixed part of model to the predicted
values of observations.
verbose
If TRUE print out intermediate results and
messages.
use.cholesky
Use the symbolic part of the Cholesky
decomposition passed in this argument.
X
For linear inverse problems the matrix that maps coefficients of the basis to the
predicted values of observations. X must be in spam format.
wU
The matrix U multiplied by the weights.
wX
The matrix X multiplied by the weights.
xnew
Matrix of locations for prediction.
xy
Order of evaluating surface coordinates. This would be used if for example
a lon-lat surface needed to be transposed as a "lat-lon" object. Usually not needed.
Znew
Values of covariates, distinct from the spatial drift
for predictions of data locations.
ZGrid
An array or list form of covariates to use for prediction. This must match the
grid.list argument. e.g. ZGrid and grid.list describe same grid.
If ZGrid is an array then the first two indices are the x and y locations in the
grid. The third index, if present, indexes the covariates. e.g. For evaluation on
a 10X15 grid and with 2 covariates. dim( ZGrid) == c(10,15, 2).
If ZGrid is a list then the components x and y should match those of grid.list and
the z component follows the shape described above for the no list case.
...
For the methods additional arguments to pass to generic methods.
For LKrig these extra arguments are interpreted as appropriate for the
LKrigSetup function and in fact they are just passed through to this function to
create an LKinfo object.
Of particular use is resetting the size of the memory allocation for the
sparse decompositions in spam. This is done by the argument
choleskyMemory. This is a list in the format of the spam memory argument for the
sparse cholesky function.
For example if a warning is printed that nnzR needs to be at least 2e5 in size, passing
choleskyMemory= list( nnzR= 2e5) will avoid this warning.
See the help on LKrigSetup for explanation of the options. Typically one would setup the LKinfo object outside of this call and
just pass a few arguments that are being varied. In particular to
use LKrig with different lambda values lambda = 1e-3 in the call
would reset the lambda value to 1e-3.
The minimal set of arguments for a 2-d rectangular
domain, however, are alpha, nlevel, NC, and a.wght.
alpha The sequence of positive weights for each level of the multiresolution process. Usually these add to one and are interpreted as the variance of the process at each level.
a.wght The weight given to the central lattice point in the
spatial autoregression (see details below).
nlevel Number of levels for the multiresolution basis. Each
level increases the number of basis functions by roughly a factor of
4.
NC Number of lattice points in first level and along the largest
dimension. e.g. if NC=5 and the domain is square the domain will
contain 25 lattice points at the first level. Note that there may be additional lattice points added as buffers outside the spatial domain
(default is 5 on each side).
Some additional arguments are V that defines a linear scaling of the
coordinates before the LatticeModel is applied. See details below.
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
rho 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 rho and a.wght, and a vector
alpha. alpha has a length of the number of multiresolution
levels but we recommend that it be constrained sum to one. The parameter rho 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, rho, 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 sigma^2/weights. Thus there
is a minimum of one new parameter from the data model: sigma. Note
that prediction only depends on the ratio lambda = sigma^2/ rho
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 sigma and rho. LKrig exploits this feature
by depending on lambda and then computing the MLEs for sigma and rho.
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.
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
.
.
.
.
.
.
.
-1
.
.
.
-1
a.wght
-1
.
.
.
-1
.
.
.
.
.
.
.
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.
Multiresolution 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, rho, to be the marginal
variance for g. So the full model for the spatial process used in
LatticeKrig is
g(x) = sqrt(rho) * 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 multiresolution
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 multiresolution 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 multiresolution 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.
Nonstationary and anisotropic modifications to the covariance
The simplest way to build in departures from isotopy 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
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 rho. Each of these can be extended to have some
spatial variation and thus provide a model for nonstationarity.
An option in specifying the marginal variance is to prescribe a
spatially varying multiplier. This component is specified by the
object rho.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: rho * predict(rho.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 rho*
predict(rho.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 the lattices at each level. For this case alpha is a list
with nlevel components and each component being a matrix with
the same dimensions as the lattice at that level. The SAR weight
matrix is taken to be the usual weights but each row is scaled by the
corresponding value in the alpha weight matrix. To be precise the
weight matrix is given in pseudo R code by
diag( 1/sqrt( c(alpha[[i]]))%*%B
leading to the precision matrix at level i of
Q.i = t(B) %*% diag( c(1/alpha[[i]]))%*%B
In this code, note that the matrix of weights, alpha[[i]] is
being stacked as a larger vector to match the implicit indexing of the
basis coefficients.
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 help(
LKrig.MRF.precision) and the commented source code for details). If
each component is a matrix then these are interpreted as the
(nonstationary) center weights for each lattice point. Finally if the
component is an array with three dimensions this specifies the center
and 8 nearest neighhors 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.
The smoothing parameter lambda and effective degrees of
freedom Consistent with other fields package functions, the two main
parameters in the model, sigma^2 and rho are parameterized as lambda =
sigma^2/rho and rho. The MLEs for rho and sigma 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)).
An Extension to linear inverse problems
Descriptions of specific functions and objects:
LKrig: Find spatial process estimate for fixed
covariance specificed 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
subdirectory 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 multiresolution 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.
# NOTE: CRAN requires that the "run" examples to execute within 5 seconds.
# so to meet this constraint many of these have been
# severely limited to be quick, typically making NC and nlevel
# unreasonably 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)
#
# this is the same as:
# LKinfoEX<- LKrigSetup( x,y, a.wght5, 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 multiresolution 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 abou
# 8 minutes
#######################################################################
## Not run:
data(CO2)
# the Ring geometry will be periodic in the first dimension and rectagular 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:
# an example of 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,verbose=TRUE )
## End(Not run)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(LatticeKrig)
Loading required package: spam
Loading required package: grid
Spam version 1.3-0 (2015-10-24) is loaded.
Type 'help( Spam)' or 'demo( spam)' for a short introduction
and overview of this package.
Help for individual functions is also obtained by adding the
suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
Attaching package: 'spam'
The following objects are masked from 'package:base':
backsolve, forwardsolve
Loading required package: fields
Loading required package: maps
# maps v3.1: updated 'world': all lakes moved to separate new #
# 'lakes' database. Type '?world' or 'news(package="maps")'. #
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/LatticeKrig/LKrig.Rd_%03d_medium.png", width=480, height=480)
> ### Name: LKrig
> ### Title: Spatial prediction and inference using a compactly supported
> ### multi-resolution basis and a lattice model for the basis
> ### coefficients.
> ### Aliases: LKrig choleskyMemory predict.LKrig predictSE.LKrig print.LKrig
> ### print.LKinfo surface.LKrig predictSurface.LKrig
> ### Keywords: spatial
>
> ### ** Examples
>
> # NOTE: CRAN requires that the "run" examples to execute within 5 seconds.
> # so to meet this constraint many of these have been
> # severely limited to be quick, typically making NC and nlevel
> # unreasonably 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)
Call:
LKrig(x = x, y = y, a.wght = 5, nlevel = 3, nu = 1, NC = 10,
lambda = 0.1)
Number of Observations: 147
Number of parameters in the fixed component 3
Effective degrees of freedom (EDF) 48.39
Standard Error of EDF estimate: 2.267
Smoothing parameter (lambda) 0.1
MLE sigma 10.53
MLE rho 1109
Nonzero entries in Ridge regression matrix 145845
Fixed part of model is a polynomial of degree 1 (m-1)
Basis function type: Radial
Basis function used: WendlandFunction
3 Levels 2699 basis functions with overlap of 2.5 (lattice units)
Level Lattice points Spacing
1 340 1.1791111
2 667 0.5895556
3 1692 0.2947778
Type of distance metric used: Euclidean
Value(s) for weighting (alpha): 0.7619048 0.1904762 0.04761905
Value(s) for lattice dependence (a): 5 5 5
Basis functions normalized so marginal process variance is stationary
> #
> # this is the same as:
> # LKinfoEX<- LKrigSetup( x,y, a.wght5, 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:
> ##D obj<- LatticeKrig( x, y)
> ##D # summary of fit and a plot of fitted surface
> ##D print( obj)
> ##D surface( obj )
> ##D US(add=TRUE)
> ##D points(x)
> ##D # prediction standard errors
> ##D 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:
> ##D
> ##D LKinfo<- LKrigSetup( x, nlevel=3, nu=1, NC=5, a.wght=5,
> ##D lambda=.01)
> ##D # maximize likelihood over lambda see help( LKrig.MLE) for details
> ##D # this assumes the value of 5 for a.wght. In many cases the fit is not
> ##D # very sensitive to the range parameter such as a.wght in this case --
> ##D # but very sensitive to lambda when varied on a log scale.
> ##D
> ##D MLE.fit<- LKrig.MLE(x,y, LKinfo=LKinfo)
> ##D MLE.fit$summary # summary of optimization over lambda.
> ##D # fit using MLE for lambda MLE function has added MLE value of lambda to
> ##D # the LKinfo object.
> ##D obj<- LKrig( x,y, LKinfo=MLE.fit$LKinfo.MLE)
> ##D print( obj)
> ##D # find prediction standard errors at locations based on fixing covariance
> ##D # at MLE's
> ##D out.se<- predictSE( obj, xnew= x)
> ##D # one could evaluate the SE on a grid to get the surface of predicted SE's
> ##D # for large grids it is better to use LKrig.sim.conditional to estimate
> ##D # the variances by Monte Carlo
> ## End(Not run)
>
> ##########################################################################
> # Use multiresolution 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:
> ##D LKinfo.temp<- LKrigSetup( x, NC=6, nu=1, nlevel=3, a.wght= 5, lambda=.5)
> ##D # starting value for lambda optimization
> ##D MLE.search<- LKrig.MLE(x,y, LKinfo=LKinfo.temp)
> ##D # this function returns an LKinfo object with the MLE for lambda included.
> ##D 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:
> ##D LKinfo<- LKrigSetup( x.CO, nlevel=1, NC=20,alpha=1, a.wght=4.1, lambda=1.0)
> ##D # lambda given here is just the starting value for MLE optimization
> ##D CO.MLE<- LKrig.MLE( x.CO,y.CO,Z=Z.CO, LKinfo=LKinfo)
> ##D obj.CO.elev<- LKrig( x.CO,y.CO,Z=Z.CO, LKinfo= CO.MLE$LKinfo.MLE)
> ##D CO.surface2<- predictSurface( obj.CO.elev, drop.Z=TRUE, nx=50, ny=50)
> ##D # pull off CO elevations at same locations on grid as the surface
> ##D data( RMelevation)
> ##D # a superset of elevations at 4km resolution
> ##D elev.surface<- interp.surface.grid( RMelevation, CO.surface2)
> ##D CO.full<- predictSurface( obj.CO.elev, ZGrid= elev.surface, nx=50, ny=50)
> ##D
> ##D # for comparison fit without elevation as a linear covariate:
> ##D CO.MLE2<- LKrig.MLE( x.CO,y.CO, LKinfo=LKinfo)
> ##D obj.CO<- LKrig( x.CO,y.CO, LKinfo= CO.MLE2$LKinfo.MLE)
> ##D # surface estimate
> ##D CO.surface<- predictSurface( obj.CO, nx=50, ny=50)
> ##D set.panel( 2,1)
> ##D coltab<- topo.colors(256)
> ##D zr<- range( c( CO.full$z), na.rm=TRUE)
> ##D image.plot( CO.surface, col=coltab, zlim =zr)
> ##D US( add=TRUE,lwd=2)
> ##D title( "MAM min temperatures without elevation")
> ##D image.plot( CO.full, col=coltab, zlim=zr)
> ##D title( "Including elevation")
> ##D US( add=TRUE,lwd=2)
> ## End(Not run)
> ## Not run:
> ##D #### a 3-d example using alternative geometry
> ##D set.seed( 123)
> ##D N<- 500
> ##D x<- matrix( runif(3* N,-1,1), ncol=3, nrow=N)
> ##D y<- 10*exp( -rdist( x, rbind( c(.5,.5,.6) ) )/.5)
> ##D LKinfo<- LKrigSetup( x,
> ##D LKGeometry = "LKBox",
> ##D nlevel = 1,
> ##D a.wght = 6.01,
> ##D NC = 5,
> ##D NC.buffer = 2,
> ##D normalize = TRUE,
> ##D choleskyMemory = list(nnzR= 2e6),
> ##D lambda = .1,
> ##D mean.neighbor=200
> ##D )
> ##D # NOTE: memory for the spam routines has been increased to
> ##D # avoid warnings
> ##D out1<- LKrig( x,y, LKinfo=LKinfo)
> ## End(Not run)
>
> ## Not run:
> ##D # MLE search over lambda and a bigger problem
> ##D # but no normalization
> ##D N<- 5e4
> ##D x<- matrix( runif(3* N,-1,1), ncol=3, nrow=N)
> ##D y<- 10*exp( -rdist( x, rbind( c(.5,.5,.6) ) )/.5)
> ##D LKinfo<- LKrigSetup( x,
> ##D LKGeometry = "LKBox",
> ##D nlevel = 1,
> ##D a.wght = 6.01,
> ##D NC = 20,
> ##D NC.buffer = 2,
> ##D normalize = FALSE,
> ##D choleskyMemory = list(nnzR= 25e6),
> ##D lambda = .1,
> ##D mean.neighbor=200
> ##D )
> ##D # add in timing
> ##D system.time( out1<- LKrig( x,y, LKinfo=LKinfo) ) # about 25 seconds
> ##D # LKinfo$normalize<- TRUE
> ##D # 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 abou
> # 8 minutes
> #######################################################################
> ## Not run:
> ##D data(CO2)
> ##D # the Ring geometry will be periodic in the first dimension and rectagular on
> ##D # second. A useful approximation for spherical data omitting the polar caps.
> ##D
> ##D LKinfo.CO2<- LKrigSetup(CO2$lon.lat, NC=100,nlevel=1, lambda=.2,
> ##D a.wght=5, alpha=1,
> ##D LKGeometry="LKRing", choleskyMemory = list(nnzR=2e6) )
> ##D print(LKinfo.CO2)
> ##D obj1<- LKrig( CO2$lon.lat,CO2$y,LKinfo=LKinfo.CO2)
> ##D # 5700+ basis functions 101X57 lattice (number of basis functions
> ##D # reduced in y direction because of a rectangular domain
> ##D obj1$trA.est # about 2900+ effective degrees of freedom
> ##D #
> ##D glist<- list( x= seq( -180,180,,200),y=seq( -80,80,,100) )
> ##D fhat<- predictSurface( obj1,grid.list=glist)
> ##D #Plot data and gap-filled estimate
> ##D set.panel(2,1)
> ##D quilt.plot(CO2$lon.lat,CO2$y,zlim=c(373,381))
> ##D title("Simulated CO2 satellite observations")
> ##D world(add=TRUE,col="magenta")
> ##D image.plot( fhat,zlim=c(373,381))
> ##D world( add=TRUE, col="magenta")
> ##D title("Gap-filled global predictions")
> ## End(Not run)
> set.panel()
plot window will lay out plots in a 1 by 1 matrix
> #########################################################################
> # 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:
> ##D data(ozone2)
> ##D x<-ozone2$lon.lat
> ##D y<- ozone2$y[16,]
> ##D # Find location that are not 'NA'.
> ##D # (LKrig is not set up to handle missing observations.)
> ##D good <- !is.na( y)
> ##D x<- x[good,]
> ##D y<- y[good]
> ##D a.wght<- 5
> ##D lambda <- 1.5
> ##D obj1<- LKrig( x,y,NC=16,nlevel=1, alpha=1, lambda=lambda, a.wght=5,
> ##D NtrA=20,iseed=122)
> ##D
> ##D # in both calls iseed is set the same so we can compare
> ##D # Monte Carlo estimates of effective degrees of freedom
> ##D obj1$trA.est
> ##D # The covariance "parameters" are all in the list LKinfo
> ##D # to create this special list outside of a call to LKrig use
> ##D LKinfo.test <- LKrigSetup( x, NC=16, nlevel=1, alpha=1.0, a.wght=5)
> ##D
> ##D # this call to mKrig should be identical to the LKrig results
> ##D # because it uses the LKrig.cov covariance with all the right parameters.
> ##D obj2<- mKrig( x,y, lambda=lambda, m=2, cov.function="LKrig.cov",
> ##D cov.args=list( LKinfo=LKinfo.test), NtrA=20, iseed=122)
> ##D # compare the two results this is also an
> ##D # example of how tests are automated in fields
> ##D # set flag to have tests print results
> ##D test.for.zero.flag<- TRUE
> ##D test.for.zero( obj1$fitted.values, obj2$fitted.values,
> ##D tag="comparing predicted values LKrig and mKrig")
> ##D # compare standard errors.
> ##D se1<- predictSE.LKrig( obj1)
> ##D se2<- predictSE.mKrig(obj2)
> ##D test.for.zero( se1, se2,
> ##D tag="comparing standard errors for LKrig and mKrig")
> ## End(Not run)
>
> ## Not run:
> ##D # an example of a linear inverse problem
> ##D # NOTE the settings in the model are small to make for a fast example
> ##D #
> ##D data(ozone2)
> ##D x<- ozone2$lon.lat
> ##D y<- ozone2$y[16,]
> ##D good<- !is.na(y)
> ##D x<- x[good,]
> ##D y<- y[good]
> ##D
> ##D LKinfo<- LKrigSetup(x, a.wght=4.5, nlevel=3, nu=1, NC=4, lambda=.01)
> ##D
> ##D # now observe a linear combination
> ##D NNDist<- LKDist(x,x, delta=1.5)
> ##D A<- NNDist
> ##D A$ra<- exp(-NNDist$ra)
> ##D # A is a weight matrix based on neighbors close by and
> ##D # in spind sparse matrix format
> ##D # now convert to spam format
> ##D A<- spind2spam(A)
> ##D
> ##D TMatrix <- get(LKinfo$fixedFunction)(x = x)
> ##D # Tmatrix is a 3 column matrix of constant and the two spatial coordinates
> ##D # i.e. a linear function of the spatial variables
> ##D PHI<- LKrig.basis(x, LKinfo)
> ##D X<- A%*% PHI
> ##D U<- A%*%TMatrix
> ##D yIndirect<- A%*%y
> ##D #
> ##D # X<- A##D
> ##D
> ##D out0<- LatticeKrig(x,y, LKinfo=LKinfo)
> ##D out1<- LatticeKrig(x,yIndirect, U=U, X=X, LKinfo=LKinfo)
> ##D
> ##D # the predict function evaluates f in this case -- not the fitted values that
> ##D # are related to the
> ##D # observations
> ##D # partial agreement but not exact -- this is because the observational models
> ##D # assume different errors
> ##D #
> ##D plot( predict(out0,x), predict( out1,x))
> ##D
> ##D out<- LKrig.sim.conditional( out1,M=5, nx=10, ny=10,verbose=TRUE )
> ## End(Not run)
>
>
>
>
>
>
> dev.off()
null device
1
>