Last data update: 2014.03.03

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.plot``contour`, `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

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 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.

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 )
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
# 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)
title( "MAM min temperatures without elevation")
image.plot( CO.full, col=coltab, zlim=zr)
title( "Including elevation")

## 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
)
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")
image.plot( fhat,zlim=c(373,381))
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.
'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)
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

# 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 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   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
> # 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     title( "MAM min temperatures without elevation")
> ##D   image.plot( CO.full, col=coltab, zlim=zr)
> ##D     title( "Including elevation")
> ## 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   image.plot( fhat,zlim=c(373,381))
> ##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
>

```