Package 'geoR'

Adapts nlm for Constraints in the Parameter Values

Description

This function adapts the R function nlm to allow for constraints (upper and/or lower bounds) in the values of the parameters.

Usage

.nlmP(objfunc, params, lower=rep(-Inf, length(params)),
      upper=rep(+Inf, length(params)), ...)

Arguments

objfunc	the function to be minimized.
params	starting values for the parameters.
lower	lower bounds for the variables. Defaults to $-Inf$.
upper	upper bounds for the variables. Defaults to $-Inf$.
...	further arguments to be passed to the function nlm.

Details

Constraints on the parameter values are internally imposed by using exponential, logarithmic, and logit transformation of the parameter values.

Value

The output is the same as for the function nlm.

Author(s)

Patrick E. Brown [email protected].
Adapted and included in geoR by
Paulo Justiniano Ribeiro Jr. [email protected]

References

Further information on the package geoR can be found at:
http://www.leg.ufpr.br/geoR/.

See Also

nlm, optim.

---

Converts an Object to the Class "geodata"

Description

The default method converts a matrix or a data-frame to an object of the class "geodata".
Objects of the class "geodata" are lists with two obligatory components: coords and data. Optional components are allowed and a typical example is a vector or matrix with covariate(s) values.

Usage

as.geodata(obj, ...)

## Default S3 method:
as.geodata(obj, coords.col = 1:2, data.col = 3, data.names = NULL, 
                   covar.col = NULL, covar.names = "obj.names",
                   units.m.col = NULL, realisations = NULL,
                   na.action = c("ifany", "ifdata", "ifcovar", "none"),
                   rep.data.action, rep.covar.action, rep.units.action,
                   ...)

## S3 method for class 'geodata'
as.data.frame(x, ..., borders = TRUE)

## S3 method for class 'geodata.frame'
as.geodata(obj, ...)

## S3 method for class 'SpatialPointsDataFrame'
as.geodata(obj, data.col = 1, ...)

is.geodata(x)

Arguments

obj	a matrix or data-frame where each line corresponds to one spatial location. It should contain values of 2D coordinates, data and, optionally, covariate(s) value(s) at the locations. A method for SpatialPointsDataFrame is also provided. It can also take an output of the function grf, see DETAILS below.
coords.col	a vector with the column numbers corresponding to the spatial coordinates.
data.col	a scalar or vector with column number(s) corresponding to the data.
data.names	optional. A string or vector of strings with names for the data columns. Only valid if there is more than one column of data. By default, takes the names from the original object.
covar.col	optional. A scalar or numeric vector with the column number(s) corresponding to the covariate(s). Alternativelly can be a character vector with the names of the covariates.
covar.names	optional. A string or vector of strings with the name(s) of the covariates. By default take the names from the original object.
units.m.col	optional. A scalar with the column number corresponding to the offset variable. Alternativelly can be a character vector with the name of the offset. This option is particularly relevant when using the package geoRglm. All values must be greater then zero.
realisations	optional. A vector indicating the realisation number or a number indicating a column in obj with the realisation indicator variable. See DETAILS below.
na.action	string defining action to be taken in the presence of NA's. The default option "ifany" excludes all points for which there are NA's in the data or covariates. The option "ifdata" excludes points for which there are NA's in the data. The default option "ifcovar" excludes all points for which there are NA's in the covariates. The option "none" do not exclude points.
rep.data.action	a string or a function. Defines action to be taken when there is more than one data at the same location. The default option "none" keeps the repeated locations, if any. The option "first" retains only the first data recorded at each location. Alternativelly a function can be passed and it will be used. For instance if mean is provided, the function will compute and return the average of the data at coincident locations. The non-default options will eliminate the repeated locations.
rep.covar.action	idem to rep.data.locations, to be applied to the covariates, if any. Defaults to the same option set for rep.data.locations.
x	an object which is tested for the class geodata.
rep.units.action	a string or a function. Defines action to be taken on the element units.m, if present when there is more than one data at the same location. The default option is the same value set for rep.data.action.
borders	logical. If TRUE the element borders in the geodata object is set as an attribute of the data-frame.
...	values to be passed for the methods.

Details

Objects of the class "geodata" contain data for geostatistical analysis using the package geoR. Storing data in this format facilitates the usage of the functions in geoR. However, conversion of objects to this class is not obligatory to carry out the analysis.

NA's are not allowed in the coordinates. By default the respective rows will not be included in the output.

Realisations
Tipically geostatistical data correspond to a unique realisation of the spatial process. However, sometimes different "realisations" are possible. For instance, if data are collected in the same area at different points in time and independence between time points is assumed, each time can be considered a different "replicate" or "realisation" of the same process. The argument realisations takes a vector indication the replication number and can be passed to other geoR functions as, for instance, likfit.

The data format is similar to the usual geodata format in geoR. Suppose there are realisations (times) $1, \ldots, J$ and for each realisations $n_1, ..., n_j$ observations are available. The coordinates for different realisations should be combined in a single $n \times 2$ object, where $n=n_1 + \ldots + n_J$. Similarly for the data vector and covariates (if any).

grf objects
If an object of the class grf is provided the functions just extracts the elements coords and data of this object.

Value

An object of the class "geodata" which is a list with two obligatory components (coords and data) and other optional components:

coords	an $n \times 2$ matrix where $n$ is the number of spatial locations.
data	a vector of length $n$, for the univariate case or, an $n \times v$ matrix or data-frame for the multivariate case, where $v$ is the number of variables.
covariates	a vector of length $n$ or an $n \times p$ matrix with covariate(s) values, where $p$ is the number of covariates. Only returned if covariates are provided.
realisations	a vector on size $n$ with the replication number. Only returned if argument realisations is provided.

Author(s)

Paulo Justiniano Ribeiro Jr. [email protected],
Peter J. Diggle [email protected].

References

Further information on the package geoR can be found at:
http://www.leg.ufpr.br/geoR/.

See Also

read.geodata for reading data from an ASCII file and list for general information on lists.

Examples

## Not run: 
## converting the data-set "topo" from the package MASS (VR's bundle)
## to the geodata format:
if(require(MASS)){
topo
topogeo <- as.geodata(topo)
names(topogeo)
topogeo
}

## End(Not run)

---

The Box-Cox Transformed Normal Distribution

Description

Functions related with the Box-Cox family of transformations. Density and random generation for the Box-Cox transformed normal distribution with mean equal to mean and standard deviation equal to sd, in the normal scale.

Usage

rboxcox(n, lambda, lambda2 = NULL, mean = 0, sd = 1)

dboxcox(x, lambda, lambda2 = NULL, mean = 0, sd = 1)

Arguments

lambda	numerical value(s) for the transformation parameter $\lambda$.
lambda2	logical or numerical value(s) of the additional transformation (see DETAILS below). Defaults to NULL.
n	number of observations to be generated.
x	a vector of quantiles (dboxcox) or an output of boxcoxfit (print, plot, lines).
mean	a vector of mean values at the normal scale.
sd	a vector of standard deviations at the normal scale.

Details

Denote $Y$ the variable at the original scale and $Y'$ the transformed variable. The Box-Cox transformation is defined by:

$Y' = \left\{ \begin{array}{ll} log(Y) \mbox{ , if $\lambda = 0$} \cr \frac{Y^\lambda - 1}{\lambda} \mbox{ , otherwise} \end{array} \right.$

.

An additional shifting parameter $\lambda_2$ can be included in which case the transformation is given by:

$Y' = \left\{ \begin{array}{ll} log(Y + \lambda_2) \mbox{ , $\lambda = 0$ } \cr \frac{(Y + \lambda_2)^\lambda - 1}{\lambda} \mbox{ , otherwise} \end{array} \right.$

.

The function rboxcox samples $Y'$ from the normal distribution using the function rnorm and backtransform the values according to the equations above to obtain values of $Y$. If necessary the back-transformation truncates the values such that $Y' \geq \frac{1}{\lambda}$ results in $Y = 0$ in the original scale. Increasing the value of the mean and/or reducing the variance might help to avoid truncation.

Value

The functions returns the following results:

rboxcox	a vector of random deviates.
dboxcox	a vector of densities.

Author(s)

Paulo Justiniano Ribeiro Jr. [email protected],
Peter J. Diggle [email protected].

References

Box, G.E.P. and Cox, D.R.(1964) An analysis of transformations. JRSS B 26:211–246.

See Also

The parameter estimation function is boxcoxfit. Other packages has BoxCox related functions such as boxcox in the package MASS and the function box.cox in the package ‘⁠car⁠’.

Examples

## Simulating data
simul <- rboxcox(100, lambda=0.5, mean=10, sd=2)
##
## Comparing models with different lambdas,
## zero  means and unit variances
curve(dboxcox(x, lambda=-1), 0, 8)
for(lambda in seq(-.5, 1.5, by=0.5))
  curve(dboxcox(x, lambda), 0, 8, add = TRUE)

---

Box-Cox transformation for geodata objects

Description

Method for Box-Cox transformation for objects of the class geodata assuming the data are independent. Computes and optionally plots profile log-likelihoods for the parameter of the Box-Cox simple power transformation $y^lambda$.

Usage

## S3 method for class 'geodata'
boxcox(object, trend = "cte", ...)

Arguments

object	an object of the class geodata. See as.geodata.
trend	specifies the mean part of the model. See trend.spatial for further details. Defaults to "cte".
...	arguments to be passed for the function boxcox.

Details

This is just a wrapper for the function boxcox facilitating its usage with geodata objects.

Notice this assume independent observations which is typically not the case for geodata objects.

Value

A list of the lambda vector and the computed profile log-likelihood vector, invisibly if the result is plotted.

See Also

boxcox for parameter estimation results for independent data and likfit for parameter estimation within the geostatistical model.

Examples

if(require(MASS)){
boxcox(wolfcamp)

data(ca20)
boxcox(ca20, trend = ~altitude)
}

---

Parameter Estimation for the Box-Cox Transformation

Description

Parameter estimation and plotting of the results for the Box-Cox transformed normal distribution.

Usage

boxcoxfit(object, xmat, lambda, lambda2 = NULL, add.to.data = 0, ...)

## S3 method for class 'boxcoxfit'
print(x, ...)

## S3 method for class 'boxcoxfit'
plot(x, hist = TRUE, data = eval(x$call$object), ...)

## S3 method for class 'boxcoxfit'
lines(x, data = eval(x$call$object), ...)

Arguments

object	a vector with the data.
xmat	a matrix with covariates values. Defaults to rep(1, length(y)).
lambda	numerical value(s) for the transformation parameter $\lambda$. Used as the initial value in the function for parameter estimation. If not provided default values are assumed. If multiple values are passed the one with highest likelihood is used as initial value.
lambda2	logical or numerical value(s) of the additional transformation (see DETAILS below). Defaults to NULL. If TRUE this parameter is also estimated and the initial value is set to the absolute value of the minimum data. A numerical value is provided it is used as the initial value. Multiple values are allowed as for lambda.
add.to.data	a constant value to be added to the data.
x	a list, typically an output of the function boxcoxfit.
hist	logical indicating whether histograms should to be plotted.
data	data values.
...	extra parameters to be passed to the minimization function optim (boxcoxfit), hist (plot) or curve (lines).

Value

The functions returns the following results:

boxcoxfit	a list with estimated parameters and results on the numerical minimization.
print.boxcoxfit	print estimated parameters. No values returned.
plot.boxcoxfit	plots histogram of the data (optional) and the model. No values returned. This function is only valid if covariates are not included in boxcoxfit.
lines.boxcoxfit	adds a line with the fitted model to the current plot. No values returned. This function is only valid if covariates are not included in boxcoxfit.

Author(s)

Paulo Justiniano Ribeiro Jr. [email protected],
Peter J. Diggle [email protected].

References

Box, G.E.P. and Cox, D.R.(1964) An analysis of transformations. JRSS B 26:211–246.

See Also

rboxcox and dboxcox for the expression and more on the Box-Cox transformation. Parameter(s) are estimated using the minimization function optim. Other packages have BoxCox related functions such as boxcox in the package MASS and the function box.cox in the package ‘⁠car⁠’.

Examples

set.seed(384)
## Simulating data
simul <- rboxcox(100, lambda=0.5, mean=10, sd=2)
## Finding the ML estimates
ml <- boxcoxfit(simul)
ml
## Ploting histogram and fitted model
plot(ml)
##
## Comparing models with different lambdas,
## zero  means and unit variances
curve(dboxcox(x, lambda=-1), 0, 8)
for(lambda in seq(-.5, 1.5, by=0.5))
  curve(dboxcox(x, lambda), 0, 8, add = TRUE)
##
## Another example, now estimating lambda2
##
simul <- rboxcox(100, lambda=0.5, mean=10, sd=2)
ml <- boxcoxfit(simul, lambda2 = TRUE)
ml
plot(ml)
##
## An example with a regression model
##
boxcoxfit(object = trees[,3], xmat = trees[,1:2])

---

Calcium content in soil samples

Description

This data set contains the calcium content measured in soil samples taken from the 0-20cm layer at 178 locations within a certain study area divided in three sub-areas. The elevation at each location was also recorded.

The first region is typically flooded during the rain season and not used as an experimental area. The calcium levels would represent the natural content in the region. The second region has received fertilisers a while ago and is typically occupied by rice fields. The third region has received fertilisers recently and is frequently used as an experimental area.

Usage

data(ca20)

Format

The object ca20 belongs to the class geodata and is a list with the following elements:

coords

a matrix with the coordinates of the soil samples.

data

calcium content measured in $mmol_c/dm^3$.

covariate

a data-frame with the covariates

altitude

a vector with the elevation of each sampling location, in meters ($m$).

area

a factor indicating the sub area to which the locations belongs.

borders

a matrix with the coordinates defining the borders of the area.

reg1

a matrix with the coordinates of the limits of the sub-area 1.

reg1

a matrix with the coordinates of the limits of the sub-area 2.

reg1

a matrix with the coordinates of the limits of the sub-area 3.

Source

The data was collected by researchers from PESAGRO and EMBRAPA-Solos, Rio de Janeiro, Brasil and provided by Dra. Maria Cristina Neves de Oliveira.

Capeche, C.L.; Macedo, J.R.; Manzatto, H.R.H.; Silva, E.F. (1997) Caracterização pedológica da fazenda Angra - PESAGRO/RIO - Estação experimental de Campos (RJ). (compact disc). In: Congresso BRASILEIRO de Ciência do Solo. 26., Informação, globalização, uso do solo; Rio de Janeiro, 1997. trabalhos. Rio de Janeiro: Embrapa/SBCS.

References

Oliveira, M.C.N. (2003) Métodos de estimação de parâmetros em modelos geoestatísticos com diferentes estruturas de covariâncias: uma aplicação ao teor de cálcio no solo. Tese de Doutorado, ESALQ/USP/Brasil.

Further information on the package geoR can be found at:
http://www.leg.ufpr.br/geoR/.

---

Calcium and magnesium content in soil samples

Description

This data set contains the calcium content measured in soil samples taken from the 0-20cm layer at 178 locations within a certain study area divided in three sub-areas. The elevation at each location was also recorded.

The first region is tipically flooded during the rain season and not used as an experimental area. The calcium levels would represent the natural content in the region. The second region has received fertilizers a while ago and is tipically occupied by rice fields. The third region has recieved fertilizers recently and is frequently used as an experimental area.

Usage

data(camg)

Format

A data frame with 178 observations on the following 10 variables.

east

east-west coordinates, in meters.

north

north-south coordinates, in meters.

elevation

elevation, in meters

region

a factor where numbers indicate different sub-regions within the area

ca020

calcium content in the 0-20cm soil layer, measured in $mmol_c/dm^3$.

mg020

calcium content in the 0-20cm soil layer, measured in $mmol_c/dm^3$.

ctc020

calcium content in the 0-20cm soil layer.

ca2040

calcium content in the 20-40cm soil layer, measured in $mmol_c/dm^3$.

mg2040

calcium content in the 20-40cm soil layer, measured in $mmol_c/dm^3$.

ctc2040

calcium content in the 20-40cm soil layer.

Details

More details about this data-set, including coordinates of the region and sub-region borders can be found in the data object ca20.

Source

The data was collected by researchers from PESAGRO and EMBRAPA-Solos, Rio de Janeiro, Brasil and provided by Dra. Maria Cristina Neves de Oliveira.

Capeche, C.L.; Macedo, J.R.; Manzatto, H.R.H.; Silva, E.F. (1997) Caracterização pedológica da fazenda Angra - PESAGRO/RIO - Estação experimental de Campos (RJ). (compact disc). In: Congresso BRASILEIRO de Ciência do Solo. 26., Informação, globalização, uso do solo; Rio de Janeiro, 1997. trabalhos. Rio de Janeiro: Embrapa/SBCS.

Examples

plot(camg[-(1:2),])
mg20 <- as.geodata(camg, data.col=6)
plot(mg20)
points(mg20)

---

Geometric Anisotropy Correction

Description

Transforms or back-transforms a set of coordinates according to the geometric anisotropy parameters.

Usage

coords.aniso(coords, aniso.pars, reverse = FALSE)

Arguments

coords	an $n \times 2$ matrix with the coordinates to be transformed.
aniso.pars	a vector with two elements, $\psi_A$ and $\psi_R$, the anisotropy angle and the anisotropy ratio, respectively. Notice that the parameters must be provided in this order. See section DETAILS below for more information on anisotropy parameters.
reverse	logical. Defaults to FALSE. If TRUE the reverse transformation is performed.

Details

Geometric anisotropy is defined by two parameters:

Anisotropy angle

defined here as the azimuth angle of the direction with greater spatial continuity, i.e. the angle between the y-axis and the direction with the maximum range.

Anisotropy ratio

defined here as the ratio between the ranges of the directions with greater and smaller continuity, i.e. the ratio between maximum and minimum ranges. Therefore, its value is always greater or equal to one.

If reverse = FALSE (the default) the coordinates are transformed from the anisotropic space to the isotropic space. The transformation consists in multiplying the original coordinates by a rotation matrix $R$ and a shrinking matrix $T$, as follows:

$X_m = X R T ,$

where $X_m$ is a matrix with the modified coordinates (isotropic space) , $X$ is a matrix with original coordinates (anisotropic space), $R$ rotates coordinates according to the anisotropy angle $\psi_A$ and $T$ shrinks the coordinates according to the anisotropy ratio $\psi_R$.

If reverse = TRUE, the back-transformation is performed, i.e. transforming the coordinates from the isotropic space to the anisotropic space by computing:

$X = X_m (R T)^{-1}$

Value

An $n \times 2$ matrix with the transformed coordinates.

Author(s)

Paulo Justiniano Ribeiro Jr. [email protected]
Peter J. Diggle [email protected].

References

Further information on the package geoR can be found at:
http://www.leg.ufpr.br/geoR/.

Examples

op <- par(no.readonly = TRUE)
par(mfrow=c(3,2))
par(mar=c(2.5,0,0,0))
par(mgp=c(2,.5,0))
par(pty="s")
## Defining a set of coordinates
coords <- expand.grid(seq(-1, 1, l=3), seq(-1, 1, l=5))
plot(c(-1.5, 1.5), c(-1.5, 1.5), xlab="", ylab="", type="n")
text(coords[,1], coords[,2], 1:nrow(coords))
## Transforming coordinates according to some anisotropy parameters
coordsA <- coords.aniso(coords, aniso.pars=c(0, 2))
plot(c(-1.5, 1.5), c(-1.5, 1.5), xlab="", ylab="", type="n")
text(coordsA[,1], coordsA[,2], 1:nrow(coords))
##
coordsB <- coords.aniso(coords, aniso.pars=c(pi/2, 2))
plot(c(-1.5, 1.5), c(-1.5, 1.5), xlab="", ylab="", type="n")
text(coordsB[,1], coordsB[,2], 1:nrow(coords))
##
coordsC <- coords.aniso(coords, aniso.pars=c(pi/4, 2))
plot(c(-1.5, 1.5), c(-1.5, 1.5), xlab="", ylab="", type="n")
text(coordsC[,1], coordsC[,2], 1:nrow(coords))
##
coordsD <- coords.aniso(coords, aniso.pars=c(3*pi/4, 2))
plot(c(-1.5, 1.5), c(-1.5, 1.5), xlab="", ylab="", type="n")
text(coordsD[,1], coordsD[,2], 1:nrow(coords))
##
coordsE <- coords.aniso(coords, aniso.pars=c(0, 5))
plot(c(-1.5, 1.5), c(-1.5, 1.5), xlab="", ylab="", type="n")
text(coordsE[,1], coordsE[,2], 1:nrow(coords))
##
par(op)

---

Operations on Coordinates

Description

Functions for shifting, zooming and envolving rectangle of a set of coordinates.

Usage

coords2coords(coords, xlim, ylim, xlim.ori, ylim.ori)

zoom.coords(x, ...)

## Default S3 method:
zoom.coords(x, xzoom, yzoom, xlim.ori, ylim.ori, xoff=0, yoff=0, ...)

## S3 method for class 'geodata'
zoom.coords(x, ...)

rect.coords(coords, xzoom = 1, yzoom=xzoom, add.to.plot=TRUE,
            quiet = FALSE, ...)

Arguments

coords, x	two column matrix or data-frame with coordinates.
xlim	range of the new x-coordinates.
ylim	range of the new y-coordinates.
xlim.ori	optional. Range of the original x-coordinates, by default the range of the original x-coordinates.
ylim.ori	optional. Range of the original y-coordinates, by default the range of the original y-coordinates.
xzoom	scalar, expanding factor in the x-direction.
yzoom	scalar, expanding factor in the y-direction.
xoff	scalar, shift in the x-direction.
yoff	scalar, shift in the y-direction.
add.to.plot	logical, if TRUE the retangle is added to the current plot.
quiet	logical, none is returned.
...	further arguments to be passed to rect.

Value

coords2coords and zoom.coords	return an object of the same type as given in the argument coords with the transformed coordinates.
rect.coords	returns a matrix with the 4 coordinates of the rectangle defined by the coordinates.

Author(s)

Paulo Justiniano Ribeiro Jr. [email protected],
Peter J. Diggle [email protected].

See Also

subarea, rect

Examples

foo <- matrix(c(4,6,6,4,2,2,4,4), nc=2)
foo1 <- zoom.coords(foo, 2)
foo1
foo2 <- coords2coords(foo, c(6,10), c(6,10))
foo2
plot(1:10, 1:10, type="n")
polygon(foo)
polygon(foo1, lty=2)
polygon(foo2, lwd=2)
arrows(foo[,1], foo[,2],foo1[,1],foo1[,2], lty=2)
arrows(foo[,1], foo[,2],foo2[,1],foo2[,2])
legend("topleft", 
       c("foo", "foo1 (zoom.coords)", "foo2 (coords2coords)"), lty=c(1,2,1), lwd=c(1,1,2))

## "zooming" part of The Gambia map
gb <- gambia.borders/1000
gd <- gambia[,1:2]/1000
plot(gb, ty="l", asp=1, xlab="W-E (kilometres)", ylab="N-S (kilometres)")
points(gd, pch=19, cex=0.5)
r1b <- gb[gb[,1] < 420,]
rc1 <- rect.coords(r1b, lty=2)

r1bn <- zoom.coords(r1b, 1.8, xoff=90, yoff=-90)
rc2 <- rect.coords(r1bn, xz=1.05)
segments(rc1[c(1,3),1],rc1[c(1,3),2],rc2[c(1,3),1],rc2[c(1,3),2], lty=3)

lines(r1bn)
r1d <- gd[gd[,1] < 420,]
r1dn <- zoom.coords(r1d, 1.7, xlim.o=range(r1b[,1],na.rm=TRUE), ylim.o=range(r1b[,2], na.rm=TRUE), 
                    xoff=90, yoff=-90)
points(r1dn, pch=19, cex=0.5)
text(450,1340, "Western Region", cex=1.5)

---

Computes Value of the Covariance Function

Description

Computes the covariances for pairs variables, given the separation distance of their locations. Options for different correlation functions are available. The results can be seen as a change of metric, from the Euclidean distances to covariances.

Usage

cov.spatial(obj, cov.model= "matern",
            cov.pars=stop("no cov.pars argument provided"),
            kappa = 0.5)

Arguments

obj	a numeric object (vector or matrix), typically with values of distances between pairs of spatial locations.
cov.model	string indicating the type of the correlation function. Available choices are: "matern", "exponential", "gaussian", "spherical", "circular", "cubic", "wave", "power", "powered.exponential", "cauchy", "gencauchy", "gneiting", "gneiting.matern", "pure.nugget". See section DETAILS for available options and expressions of the correlation functions.
cov.pars	a vector with 2 elements or an $ns \times 2$ matrix with the covariance parameters. The first element (if a vector) or first column (if a matrix) corresponds to the variance parameter $\sigma^2$. The second element or column corresponds to the range parameter $\phi$ of the correlation function. If a matrix is provided, each row corresponds to the parameters of one spatial structure (see DETAILS below).
kappa	numerical value for the additional smoothness parameter of the correlation function. Only required by the following correlation functions: "matern", "powered.exponential", "cauchy", "gencauchy" and "gneiting.matern".

Details

Covariance functions return the value of the covariance $C(h)$ between a pair variables located at points separated by the distance $h$. The covariance function can be written as a product of a variance parameter $\sigma^2$ times a positive definite correlation function $\rho(h)$:

$C(h) = \sigma^2 \rho(h).$

The expressions of the covariance functions available in geoR are given below. We recommend the LaTeX (and/or the corresponding .dvi, .pdf or .ps) version of this document for better visualization of the formulas.

Denote $\phi$ the basic parameter of the correlation function and name it the range parameter. Some of the correlation functions will have an extra parameter $\kappa$, the smoothness parameter. $K_\kappa(x)$ denotes the modified Bessel function of the third kind of order $\kappa$. See documentation of the function besselK for further details. In the equations below the functions are valid for $\phi>0$ and $\kappa>0$, unless stated otherwise.

cauchy

$\rho(h) = [1+(\frac{h}{\phi})^2]^{-\kappa}$

gencauchy (generalised Cauchy)

$\rho(h) = [1+(\frac{h}{\phi})^{\kappa_{2}}]^{-{\kappa_1}/{\kappa_2}}, \kappa_1 > 0, 0 < \kappa_2 \leq 2$

circular
Let $\theta = \min(\frac{h}{\phi},1)$ and

$g(h)= 2\frac{(\theta\sqrt{1-\theta^2}+ \sin^{-1}\sqrt{\theta})}{\pi}.$

Then, the circular model is given by:

$\rho(h) = \left\{ \begin{array}{ll} 1 - g(h) \mbox{ , if $h < \phi$}\cr 0 \mbox{ , otherwise} \end{array} \right.$

cubic

$\rho(h) = \left\{ \begin{array}{ll} 1 - [7(\frac{h}{\phi})^2 - 8.75(\frac{h}{\phi})^3 + 3.5(\frac{h}{\phi})^5-0.75(\frac{h}{\phi})^7] \mbox{ , if $h<\phi$} \cr 0 \mbox{ , otherwise.} \end{array} \right.$

gaussian

$\rho(h) = \exp[-(\frac{h}{\phi})^2]$

exponential

$\rho(h) = \exp(-\frac{h}{\phi})$

matern

$\rho(h) = \frac{1}{2^{\kappa-1}\Gamma(\kappa)}(\frac{h}{\phi})^\kappa K_{\kappa}(\frac{h}{\phi})$

spherical

$\rho(h) = \left\{ \begin{array}{ll} 1 - 1.5\frac{h}{\phi} + 0.5(\frac{h}{\phi})^3 \mbox{ , if $h$ < $\phi$} \cr 0 \mbox{ , otherwise} \end{array} \right.$

power (and linear)
The parameters of the this model $\sigma^2$ and $\phi$ can not be interpreted as partial sill and range as for the other models. This model implies an unlimited dispersion and, therefore, has no sill and corresponds to a process which is only intrinsically stationary. The variogram function is given by:

$\gamma(h) = \sigma^2 {h}^{\phi} \mbox{ , } 0 < \phi < 2, \sigma^2 > 0$

Since the corresponding process is not second order stationary the covariance and correlation functions are not defined. For internal calculations the geoR functions uses the fact the this model possesses locally stationary representations with covariance functions of the form:

$C_(h) = \sigma^2 (A - h^\phi)$

, where $A$ is a suitable constant as given in Chiles & Delfiner (pag. 511, eq. 7.35).

The linear model corresponds a particular case with $\phi = 1$.

powered.exponential (or stable)

$\rho(h) = \exp[-(\frac{h}{\phi})^\kappa] \mbox{ , } 0 < \kappa \leq 2$

gneiting

$C(h)=\left(1 + 8 sh + 25 (sh)^2 + 32 (sh)^3\right)(1-sh)^8 1_{[0,1]}(sh)$

where $s=0.301187465825$. For further details see documentation of the function CovarianceFct in the package RandomFields from where we extract the following :
It is an alternative to the gaussian model since its graph is visually hardly distinguishable from the graph of the Gaussian model, but possesses neither the mathematical and nor the numerical disadvantages of the Gaussian model.

gneiting.matern
Let $\alpha=\phi\kappa_2$, $\rho_m(\cdot)$ denotes the $\mbox{Mat\'{e}rn}$ model and $\rho_g(\cdot)$ the Gneiting model. Then the $\mbox{Gneiting-Mat\'{e}rn}$ is given by

$\rho(h) = \rho_g(h|\phi=\alpha) \, \rho_m(h|\phi=\phi,\kappa=\kappa_1)$

wave

$\rho(h) = \frac{\phi}{h}\sin(\frac{h}{\phi})$

pure.nugget

$\rho(h) = k$

where k is a constant value. This model corresponds to no spatial correlation.

Nested models Models with several structures usually called nested models in the geostatistical literature are also allowed. In this case the argument cov.pars takes a matrix and cov.model and lambda can either have length equal to the number of rows of this matrix or length 1. For the latter cov.model and/or lambda are recycled, i.e. the same value is used for all structures.

Value

The function returns values of the covariances corresponding to the given distances. The type of output is the same as the type of the object provided in the argument obj, typically a vector, matrix or array.

Author(s)

Paulo J. Ribeiro Jr. [email protected],
Peter J. Diggle [email protected].

References

For a review on correlation functions:
Schlather, M. (1999) An introduction to positive definite functions and to unconditional simulation of random fields. Technical report ST 99-10, Dept. of Maths and Statistics, Lancaster University.

Chilès, J.P. and Delfiner, P. (1999) Geostatistics: Modelling Spatial Uncertainty, Wiley.

Further information on the package geoR can be found at:
http://www.leg.ufpr.br/geoR/.

See Also

matern for computation of the $\mbox{Mat\'{e}rn}$ model, besselK for computation of the Bessel function and varcov.spatial for computations related to the covariance matrix.

Examples

#
# Variogram models with the same "practical" range:
#
v.f <- function(x, ...){1-cov.spatial(x, ...)}
#
curve(v.f(x, cov.pars=c(1, .2)), from = 0, to = 1,
      xlab = "distance", ylab = expression(gamma(h)),
      main = "variograms with equivalent \"practical range\"")
curve(v.f(x, cov.pars = c(1, .6), cov.model = "sph"), 0, 1,
      add = TRUE, lty = 2)
curve(v.f(x, cov.pars = c(1, .6/sqrt(3)), cov.model = "gau"),
      0, 1, add = TRUE, lwd = 2)
legend("topleft", c("exponential", "spherical", "gaussian"),
       lty=c(1,2,1), lwd=c(1,1,2))
#
# Matern models with equivalent "practical range"
# and varying smoothness parameter
#
curve(v.f(x, cov.pars = c(1, 0.25), kappa = 0.5),from = 0, to = 1,
      xlab = "distance", ylab = expression(gamma(h)), lty = 2,
      main = "models with equivalent \"practical\" range")
curve(v.f(x, cov.pars = c(1, 0.188), kappa = 1),from = 0, to = 1,
      add = TRUE)      
curve(v.f(x, cov.pars = c(1, 0.14), kappa = 2),from = 0, to = 1,
      add = TRUE, lwd=2, lty=2)      
curve(v.f(x, cov.pars = c(1, 0.117), kappa = 2),from = 0, to = 1,
      add = TRUE, lwd=2)      
legend("bottomright",
       expression(list(kappa == 0.5, phi == 0.250), 
         list(kappa == 1, phi == 0.188), list(kappa == 2, phi == 0.140),
         list(kappa == 3, phi == 0.117)), lty=c(2,1,2,1), lwd=c(1,1,2,2))
# plotting a nested variogram model
curve(v.f(x, cov.pars = rbind(c(.4, .2), c(.6,.3)),
          cov.model = c("sph","exp")), 0, 1, ylab='nested model')

---

Locates duplicated coordinates

Description

This funtions takes an object with 2-D coordinates and returns the positions of the duplicated coordinates. Also sets a method for duplicated

Usage

dup.coords(x, ...)
## Default S3 method:
dup.coords(x, ...)
## S3 method for class 'geodata'
dup.coords(x, incomparables, ...)
## S3 method for class 'geodata'
duplicated(x, incomparables, ...)

Arguments

x	a two column numeric matrix or data frame.
incomparables	unused. Just for compatibility with the generic function duplicated.
...	arguments passed to sapply. If simplify = TRUE (default) results are returned as an array if possible (when the number of replicates are the same at each replicated location)

Value

Function and methods returns NULL if there are no duplicates locations.

Otherwise, the default method returns a list where each component is a vector with the positions or the rownames, if available, of the duplicates coordinates.

The method for geodata returns a data-frame with rownames equals to the positions of the duplicated coordinates, the first column is a factor indicating duplicates and the remaning are output of as.data.frame.geodata.

Author(s)

Paulo Justiniano Ribeiro Jr. [email protected]
Peter J. Diggle [email protected].

See Also

as.geodata for the definition of geodata class, duplicated for the base function to identify duplicated values and jitterDupCoords for a function which jitters duplicated coordinates.

Examples

## simulating data
dt <- grf(30, cov.p=c(1, .3)) 
## "forcing" some duplicated locations
dt$coords[4,] <- dt$coords[14,] <- dt$coords[24,] <- dt$coords[2,]
dt$coords[17,] <- dt$coords[23,] <- dt$coords[8,]
## output of the method for geodata
dup.coords(dt)
## which is the same as a method for duplicated()
duplicated(dt)
## the default method:
dup.coords(dt$coords)

---

Surface Elevations

Description

Surface elevation data taken from Davis (1972). An onject of the class geodata with elevation values at 52 locations.

Usage

data(elevation)

Format

An object of the class geodata which is a list with the following elements:

coords

x-y coordinates (multiples of 50 feet).

data

elevations (multiples of 10 feet).

Source

Davis, J.C. (1973) Statistics and Data Analysis in Geology. Wiley.

Examples

summary(elevation)
plot(elevation)

---

Interactive Variogram Estimation

Description

Function to fit an empirical variogram "by eye" using an interactive Tcl-Tk interface.

Usage

eyefit(vario, silent = FALSE)

Arguments

vario	An empirical variogram object as returned by the function variog.
silent	logical indicating wheather or not the fitted variogram must be returned.

Value

Returns a list of list with the model parameters for each of the saved fit(s).

Author(s)

Andreas Kiefer [email protected]
Paulo Justiniano Ribeiro Junior [email protected].

See Also

variofit for least squares variogram fit, likfit for likelihood based parameter estimation and krige.bayes to obtain the posterior distribution for the model parameters.

---

Gambia Malaria Data

Description

Malaria prevalence in children recorded at villages in The Gambia, Africa.

Usage

data(gambia)

Format

Two objects are made available:

1. gambia
   A data frame with 2035 observations on the following 8 variables.

   x

   x-coordinate of the village (UTM).

   y

   y-coordinate of the village (UTM).

   pos

   presence (1) or absence (0) of malaria in a blood sample taken from the child.

   age

   age of the child, in days

   netuse

   indicator variable denoting whether (1) or not (0) the child regularly sleeps under a bed-net.

   treated

   indicator variable denoting whether (1) or not (0) the bed-net is treated (coded 0 if netuse=0).

   green

   satellite-derived measure of the green-ness of vegetation in the immediate vicinity of the village (arbitrary units).

   phc

   indicator variable denoting the presence (1) or absence (0) of a health center in the village.

2. gambia.borders
   A data frame with 2 variables:

   x

   x-coordinate of the country borders.

   y

   y-coordinate of the country borders.

References

Thomson, M., Connor, S., D Alessandro, U., Rowlingson, B., Diggle, P., Cresswell, M. & Greenwood, B. (1999). Predicting malaria infection in Gambian children from satellite data and bednet use surveys: the importance of spatial correlation in the interpretation of results. American Journal of Tropical Medicine and Hygiene 61: 2–8.

Diggle, P., Moyeed, R., Rowlingson, B. & Thomson, M. (2002). Childhood malaria in The Gambia: a case-study in model-based geostatistics, Applied Statistics.

Examples

plot(gambia.borders, type="l", asp=1)
points(gambia[,1:2], pch=19)
# a built-in function for a zoomed map
gambia.map()
# Building a "village-level" data frame
ind <- paste("x",gambia[,1], "y", gambia[,2], sep="")
village <- gambia[!duplicated(ind),c(1:2,7:8)]
village$prev <- as.vector(tapply(gambia$pos, ind, mean))
plot(village$green, village$prev)

---

Defunct Functions in the Package geoR

Description

The functions listed here are no longer part of the package geoR as they are no longer needed.

Usage

geoRdefunct()

Details

The following functions are now defunct:

olsfit

functionality incorporated by variofit starting from package version ‘⁠1.0-6⁠’.

wlsfit

functionality incorporated by variofit starting from package version ‘⁠1.0-6⁠’.

likfit.old

functionality incorporated by likfit starting from package version ‘⁠1.0-6⁠’. The related functions were also made defunct:
likfit.nospatial, loglik.spatial, proflik.nug, proflik.phi, proflik.ftau.

distdiag

functionally is redundant with dist.

See Also

variofit

---

Computes global variance

Description

Global variance computation for a set of locations using the covarianve model

Usage

globalvar(geodata, locations, coords = geodata$coords, krige)

Arguments

geodata	an object of the class geodata
locations	n by 2 matrix with a set of locations, typically a prediction grid
coords	data coordinates
krige	a list defining the model components and the type of kriging. It can take an output to a call to krige.control or a list with elements as for the arguments in krige.control.

Value

An scalar with the value of the global variance

Author(s)

Paulo Justiniano Ribeiro Jr. [email protected],
Peter J. Diggle [email protected].

References

Isaaks, E.S and Srivastava, R.M. (1989) An Introduction to Applied Geostatistics, pag. 508, eq. 20.7. Oxford University Press.

See Also

krige.conv for the kriging algorithm.

---

Simulation of Gaussian Random Fields

Description

grf() generates (unconditional) simulations of Gaussian random fields for given covariance parameters.

Usage

grf(n, grid = "irreg", nx, ny, xlims = c(0, 1), ylims = c(0, 1),
    borders, nsim = 1, cov.model = "matern",
    cov.pars = stop("missing covariance parameters sigmasq and phi"), 
    kappa = 0.5, nugget = 0, lambda = 1, aniso.pars,
    mean = 0, method, messages)

Arguments

n	number of points (spatial locations) in each simulations.
grid	optional. An $n \times 2$ matrix with coordinates of the simulated data.
nx	optional. Number of points in the X direction.
ny	optional. Number of points in the Y direction.
xlims	optional. Limits of the area in the X direction. Defaults to $[0,1]$.
ylims	optional. Limits of the area in the Y direction. Defaults to $[0,1]$.
borders	optional. Typically a two coluns matrix especifying a polygon. See DETAILS below.
nsim	Number of simulations. Defaults to 1.
cov.model	correlation function. See cov.spatial for further details. Defaults to the exponential model.
cov.pars	a vector with 2 elements or an $n \times 2$ matrix with values of the covariance parameters $\sigma^2$ (partial sill) and $\phi$ (range parameter). If a vector, the elements are the values of $\sigma^2$ and $\phi$, respectively. If a matrix, corresponding to a model with several structures, the values of $\sigma^2$ are in the first column and the values of $\phi$ are in the second.
kappa	additional smoothness parameter required only for the following correlation functions: "matern", "powered.exponential", "cauchy" and "gneiting.matern". More details on the documentation for the function cov.spatial.
nugget	the value of the nugget effect parameter $\tau^2$.
lambda	value of the Box-Cox transformation parameter. The value $\lambda = 1$ corresponds to no transformation, the default. For any other value of $\lambda$ Gaussian data is simulated and then transformed.
aniso.pars	geometric anisotropy parameters. By default an isotropic field is assumed and this argument is ignored. If a vector with 2 values is provided, with values for the anisotropy angle $\psi_A$ (in radians) and anisotropy ratio $\psi_A$, the coordinates are transformed, the simulation is performed on the isotropic (transformed) space and then the coordinates are back-transformed such that the resulting field is anisotropic. Coordinates transformation is performed by the function coords.aniso.
mean	a numerical vector, scalar or the same length of the data to be simulated. Defaults to zero.
method	simulation method with options for "cholesky", "svd", "eigen". Defaults to the Cholesky decomposition. See section DETAILS below.
messages	logical, indicating whether or not status messages are printed on the screen (or output device) while the function is running. Defaults to TRUE.

Details

For the methods "cholesky", "svd" and "eigen" the simulation consists of multiplying a vector of standardized normal deviates by a square root of the covariance matrix. The square root of a matrix is not uniquely defined. These three methods differs in the way they compute the square root of the (positive definite) covariance matrix.

The argument borders, if provides takes a polygon data set following argument poly for the splancs' function csr, in case of grid="reg" or gridpts, in case of grid="irreg". For the latter the simulation will have approximately “n” points.

Value

grf returns a list with the components:

coords	an $n \times 2$ matrix with the coordinates of the simulated data.
data	a vector (if nsim = 1) or a matrix with the simulated values. For the latter each column corresponds to one simulation.
cov.model	a string with the name of the correlation function.
nugget	the value of the nugget parameter.
cov.pars	a vector with the values of $\sigma^2$ and $\phi$, respectively.
kappa	value of the parameter $\kappa$.
lambda	value of the Box-Cox transformation parameter $\lambda$.
aniso.pars	a vector with values of the anisotropy parameters, if provided in the function call.
method	a string with the name of the simulation method used.
sim.dim	a string "1d" or "2d" indicating the spatial dimension of the simulation.
.Random.seed	the random seed by the time the function was called.
messages	messages produced by the function describing the simulation.
call	the function call.

Author(s)

Paulo Justiniano Ribeiro Jr. [email protected],
Peter J. Diggle [email protected].

References

Wood, A.T.A. and Chan, G. (1994) Simulation of stationary Gaussian process in $[0,1]^d$. Journal of Computatinal and Graphical Statistics, 3, 409–432.

Schlather, M. (1999) Introduction to positive definite functions and to unconditional simulation of random fields. Tech. Report ST–99–10, Dept Maths and Stats, Lancaster University.

Schlather, M. (2001) Simulation and Analysis of Random Fields. R-News 1 (2), p. 18-20.

Further information on the package geoR can be found at:
http://www.leg.ufpr.br/geoR/.

See Also

plot.grf and image.grf for graphical output, coords.aniso for anisotropy coordinates transformation and chol, svd and eigen for methods of matrix decomposition.

Examples

sim1 <- grf(100, cov.pars = c(1, .25))
# a display of simulated locations and values
points(sim1)   
# empirical and theoretical variograms
plot(sim1)
## alternative way
plot(variog(sim1, max.dist=1))
lines.variomodel(sim1)
#
# a "smallish" simulation
sim2 <- grf(441, grid = "reg", cov.pars = c(1, .25)) 
image(sim2)
##
## 1-D simulations using the same seed and different noise/signal ratios
##
set.seed(234)
sim11 <- grf(100, ny=1, cov.pars=c(1, 0.25), nug=0)
set.seed(234)
sim12 <- grf(100, ny=1, cov.pars=c(0.75, 0.25), nug=0.25)
set.seed(234)
sim13 <- grf(100, ny=1, cov.pars=c(0.5, 0.25), nug=0.5)
##
par.ori <- par(no.readonly = TRUE)
par(mfrow=c(3,1), mar=c(3,3,.5,.5))
yl <- range(c(sim11$data, sim12$data, sim13$data))
image(sim11, type="l", ylim=yl)
image(sim12, type="l", ylim=yl)
image(sim13, type="l", ylim=yl)
par(par.ori)

## simulating within borders
data(parana)
pr1 <- grf(100, cov.pars=c(200, 40), borders=parana$borders, mean=500)
points(pr1)
pr1 <- grf(100, grid="reg", cov.pars=c(200, 40), borders=parana$borders)
points(pr1)
pr1 <- grf(100, grid="reg", nx=10, ny=5, cov.pars=c(200, 40), borders=parana$borders)
points(pr1)

---

Head observations in a regional confined aquifer

Description

Measurements of potentiometric head at 29 locations in a regional confined sandstone aquifer. Extract from Kitanidis' book.

Usage

data(head)

Format

An object of the class geodata which is a list with the elements:

coords

coordinates of the data location.

data

the data vector with head measurements (feet).

Source

Kitanidis, P.K. Introduction to geostatistics - applications in hidrology (1997). Cambridge University Press.

Examples

summary(head)
plot(head)

---

Plots Sample from Posterior Distributions

Description

Plots histograms and/or density estimation with samples from the posterior distribution of the model parameters

Usage

## S3 method for class 'krige.bayes'
hist(x, pars, density.est = TRUE, histogram = TRUE, ...)

Arguments

x	an object of the class krige.bayes, with an output of the funtions krige.bayes.
pars	a vector with the names of one or more of the model parameters. Defaults to all model parameters. Setting to -1 excludes the intercept.
density.est	logical indication whether a line with the density estimation should be added to the plot.
histogram	logical indicating whether the histogram is included in the plot.
...	further arguments for the plotting functions and or for the density estimation.

Value

Produces a plot in the currently graphics device.
Returns a invisible list with the components:

histogram	with the output of the function hist for each parameter
density.estimation	with the output of the function density for each parameter

Author(s)

Paulo J. Ribeiro Jr. [email protected],
Peter J. Diggle [email protected].

See Also

krige.bayes, hist, density.

Examples

## See documentation for krige.bayes()

---

Data for spatial analysis of experiments

Description

The hoef data frame has 25 rows and 5 columns.
The data consists of a uniformity trial for which artificial treatment effects were assign to the plots.

Usage

data(hoef)

Format

This data frame contains the following columns:

x1

x-coordinate of the plot.

x2

y-coordinate of the plot.

dat

the artificial data.

trat

the treatment number.

ut

the data from the uniformity trial, without the treatment effect.

Details

The treatment effects assign to the plots are:

• Treatment 1: $\tau_1 = 0$

• Treatment 2: $\tau_2 = -3$

• Treatment 3: $\tau_3 = -5$

• Treatment 4: $\tau_4 = 6$

• Treatment 5: $\tau_5 = 6$

Source

Ver Hoef, J.M. & Cressie, N. (1993) Spatial statistics: analysis of field experiments. In Scheiner S.M. and Gurevitch, J. (Eds) Design and Analysis of Ecological Experiments. Chapman and Hall.

Examples

hoef.geo <- as.geodata(hoef, covar.col=4)
summary(hoef)
summary(hoef.geo)
points(hoef.geo, cex.min=2, cex.max=2, pt.div="quintiles")

---

Image, Contour or Perspective Plot of Simulated Gaussian Random Field

Description

Methods for image, contour or perspective plot of a realisation of a Gaussian random field, simulated using the function grf.

Usage

## S3 method for class 'grf'
image(x, sim.number = 1, borders, x.leg, y.leg, ...)
## S3 method for class 'grf'
contour(x, sim.number = 1, borders, filled = FALSE, ...)
## S3 method for class 'grf'
persp(x, sim.number = 1, borders, ...)

Arguments

x	an object of the class grf, typically an output of the function grf.
sim.number	simulation number. Indicates the number of the simulation top be plotted. Only valid if the object contains more than one simulation. Defaults to 1.
borders	optional. Typically a two coluns matrix especifying a polygon. Points outside the borders will be set no NA
x.leg, y.leg	limits for the legend in the horizontal and vertical directions.
filled	logical. If FALSE the function contour is used otherwise filled.contour. Defaults to FALSE.
...	further arguments to be passed to the functions image, contour or persp.

Value

An image or perspective plot is produced on the current graphics device. No values are returned.

Author(s)

Paulo Justiniano Ribeiro Jr. [email protected],
Peter J. Diggle [email protected].

References

Further information about the package geoR can be found at:
http://www.leg.ufpr.br/geoR/.

See Also

grf for simulation of Gaussian random fields, image and persp for the generic plotting functions.

Examples

# generating 4 simulations of a Gaussian random field
sim <- grf(441, grid="reg", cov.pars=c(1, .25), nsim=4)
op <- par(no.readonly = TRUE)
par(mfrow=c(2,2), mar=c(3,3,1,1), mgp = c(2,1,0))
for (i in 1:4)
  image(sim, sim.n=i)
par(op)

---

Plots Results of the Predictive Distribution

Description

This function produces an image or perspective plot of a selected element of the predictive distribution returned by the function krige.bayes.

Usage

## S3 method for class 'krige.bayes'
image(x, locations, borders,
                  values.to.plot=c("mean", "variance",
                            "mean.simulations", "variance.simulations",
                            "quantiles", "probabilities", "simulation"),
                  number.col, coords.data, x.leg, y.leg, messages, ...)
## S3 method for class 'krige.bayes'
contour(x, locations, borders, 
                  values.to.plot = c("mean", "variance",
                       "mean.simulations", "variance.simulations",
                       "quantiles", "probabilities", "simulation"),
                  filled=FALSE, number.col, coords.data,
                  x.leg, y.leg, messages, ...)
## S3 method for class 'krige.bayes'
persp(x, locations, borders,
                  values.to.plot=c("mean", "variance",
                       "mean.simulations", "variance.simulations",
                       "quantiles", "probabilities", "simulation"),
                  number.col, messages, ...)

Arguments

x	an object of the class krige.bayes, typically an output of the function krige.bayes.
locations	an $n \times 2$ matrix with the coordinates of the prediction locations, which should define a regular grid in order to be plotted by image or persp. By default does not need to be provided and evaluates the attribute "prediction.locations" from the input object.
borders	an $n \times 2$ matrix with the coordinates defining the borders of a region inside the grid defined by locations. Elements in the argument values are assigned to locations internal to the borders and NA's to the external ones.
values.to.plot	select the element of the predictive distribution to be plotted. See DETAILS below.
filled	logical. If FALSE the function contour is used otherwise filled.contour. Defaults to FALSE.
number.col	Specifies the number of the column to be plotted. Only used if previous argument is set to one of "quantiles", "probabilities" or "simulation".
coords.data	optional. If an $n \times 2$ matrix with the data coordinates is provided, points indicating the data locations are included in the plot.
x.leg, y.leg	limits for the legend in the horizontal and vertical directions.
messages	logical, if TRUE status messages are printed while running the function.
...	extra arguments to be passed to the plotting function image or persp.

Details

The function krige.bayes returns summaries and other results about the predictive distributions. The argument values.to.plot specifies which result will be plotted. It can be passed to the function in two different forms:

• a vector with the object containing the values to be plotted, or

• one of the following options: "moments.mean", "moments.variance", "mean.simulations", "variance.simulations", "quantiles", "probability" or "simulation".

For the last three options, if the results are stored in matrices, a column number must be provided using the argument number.col.

The documentation for the function krige.bayes provides further details about these options.

Value

An image or persp plot is produced on the current graphics device. No values are returned.

Author(s)

Paulo J. Ribeiro Jr. [email protected],
Peter J. Diggle [email protected].

References

Further information on the package geoR can be found at:
http://www.leg.ufpr.br/geoR/.

See Also

krige.bayes for Bayesian Kriging computations and, image and persp for the generic plotting functions.

Examples

#See examples in the documentation for the function krige.bayes().

---

Image or Perspective Plot with Kriging Results

Description

Plots image or perspective plots with results of the kriging calculations.

Usage

## S3 method for class 'kriging'
image(x, locations, borders, values = x$predict,
              coords.data, x.leg, y.leg, ...)
## S3 method for class 'kriging'
contour(x, locations, borders, values = x$predict,
              coords.data, filled=FALSE, ...)
## S3 method for class 'kriging'
persp(x, locations, borders, values = x$predict, ...)

Arguments

x	an object of the class kriging, typically with the output of the functions krige.conv or ksline.
locations	an $n \times 2$ matrix with the coordinates of the prediction locations, which should define a regular grid in order to be plotted by image or persp. By default does not need to be provided and evaluates the attribute "prediction.locations" from the input object.
borders	an $n \times 2$ matrix with the coordinates defining the borders of a region inside the grid defined by locations. Elements in the argument values are assigned to locations internal to the borders and NA's to the external ones.
values	a vector with values to be plotted. Defaults to obj$predict.
coords.data	optional. If an $n \times 2$ matrix with the data coordinates is provided, points indicating the data locations are included in the plot.
x.leg, y.leg	limits for the legend in the horizontal and vertical directions.
filled	logical. If FALSE the function contour is used otherwise filled.contour. Defaults to FALSE.
...	further arguments to be passed to the functions image, contour, filled.contour, persp or legend.krige. For instance, the argument zlim can be used to set the the minimum and maximum ‘z’ values for which colors should be plotted. See documentation for those function for possible arguments.

Details

plot1d and prepare.graph.kriging are auxiliary functions called by the others.

Value

An image or perspective plot is produced o the current graphics device. No values are returned.

Author(s)

Paulo J. Ribeiro Jr. [email protected],
Peter J. Diggle [email protected].

References

Further information on the package geoR can be found at:
http://www.leg.ufpr.br/geoR/.

See Also

krige.conv and ksline for kriging calculations. Documentation for image, contour, filled.contour and persp contain basic information on the plotting functions.

Examples

loci <- expand.grid(seq(0,1,l=51), seq(0,1,l=51))
kc <- krige.conv(s100, loc=loci,
                 krige=krige.control(cov.pars=c(1, .25)))
image(kc)
contour(kc)
image(kc)
contour(kc, add=TRUE, nlev=21)
persp(kc, theta=20, phi=20)
contour(kc, filled=TRUE)
contour(kc, filled=TRUE, color=terrain.colors)
contour(kc, filled=TRUE, col=gray(seq(1,0,l=21)))
# adding data locations
image(kc, coords.data=s100$coords)
contour(kc,filled=TRUE,coords.data=s100$coords,color=terrain.colors)
#
# now dealing with borders
#
bor <- matrix(c(.4,.1,.3,.9,.9,.7,.9,.7,.3,.2,.5,.8),
              ncol=2)
# plotting just inside borders
image(kc, borders=bor)
contour(kc, borders=bor)
image(kc, borders=bor)
contour(kc, borders=bor, add=TRUE)
contour(kc, borders=bor, filled=TRUE, color=terrain.colors)
# kriging just inside borders
kc1 <- krige.conv(s100, loc=loci,
                 krige=krige.control(cov.pars=c(1, .25)),
                 borders=bor)
image(kc1)
contour(kc1)
# avoiding the borders 
image(kc1, borders=NULL)
contour(kc1, borders=NULL)

op <- par(no.readonly = TRUE)
par(mfrow=c(1,2), mar=c(3,3,0,0), mgp=c(1.5, .8,0))
image(kc)
image(kc, val=sqrt(kc$krige.var))

# different ways to add the legends and pass arguments:
image(kc, ylim=c(-0.2, 1), x.leg=c(0,1), y.leg=c(-0.2, -0.1))
image(kc, val=kc$krige.var, ylim=c(-0.2, 1))
legend.krige(y.leg=c(-0.2,-0.1), x.leg=c(0,1), val=sqrt(kc$krige.var))

image(kc, ylim=c(-0.2, 1), x.leg=c(0,1), y.leg=c(-0.2, -0.1), cex=1.5)
image(kc, ylim=c(-0.2, 1), x.leg=c(0,1), y.leg=c(-0.2, -0.1), offset.leg=0.5)

image(kc, xlim=c(0, 1.2))
legend.krige(x.leg=c(1.05,1.1), y.leg=c(0,1), kc$pred, vert=TRUE)
image(kc, xlim=c(0, 1.2))
legend.krige(x.leg=c(1.05,1.1), y.leg=c(0,1),kc$pred, vert=TRUE, off=1.5, cex=1.5)

par(op)

---

The (Scaled) Inverse Chi-Squared Distribution

Description

Density and random generation for the scaled inverse chi-squared ($\chi^2_{ScI}$) distribution with df degrees of freedom and optional non-centrality parameter scale.

Usage

dinvchisq(x, df, scale, log = FALSE)
rinvchisq(n, df, scale = 1/df)

Arguments

x	vector of quantiles.
n	number of observations. If length(n) > 1, the length is taken to be the number required.
df	degrees of freedom.
scale	scale parameter.
log	logical; if TRUE, densities d are given as log(d).

Details

The inverse chi-squared distribution with df$= n$ degrees of freedom has density

$f(x) = \frac{1}{{2}^{n/2} \Gamma (n/2)} {(1/x)}^{n/2+1} {e}^{-1/(2x)}$

for $x > 0$. The mean and variance are $\frac{1}{(n-2)}$ and $\frac{2}{(n-4)(n-2)^2}$.

The non-central chi-squared distribution with df$= n$ degrees of freedom and non-centrality parameter scale $= S^2$ has density

$f(x) = \frac{{n/2}^{n/2}}{\Gamma (n/2)} S^n {(1/x)}^{n/2+1} {e}^{-(n S^2)/(2x)}$

, for $x \ge 0$. The first is a particular case of the latter for $\lambda = n/2$.

Value

dinvchisq gives the density and rinvchisq generates random deviates.

See Also

rchisq for the chi-squared distribution which is the basis for this function.

Examples

set.seed(1234); rinvchisq(5, df=2)
set.seed(1234); 1/rchisq(5, df=2)

set.seed(1234); rinvchisq(5, df=2, scale=5)
set.seed(1234); 5*2/rchisq(5, df=2)

## inverse Chi-squared is a particular case
x <- 1:10
all.equal(dinvchisq(x, df=2), dinvchisq(x, df=2, scale=1/2))

---

Data from Isaaks and Srisvastava's book

Description

Toy example used in the book An Introduction to Geostatistics to illustrate the effects of different models and parameters in the kriging results when predicting at a given point.

Usage

data(isaaks)

Format

An object of the class geodata which is a list with the elements:

coords

coordinates of the data location.

data

the data vector.

x0

coordinate of the prediction point.

Source

Isaaks, E.H. & Srisvastava, R.M. (1989) An Introduction to Applied Geostatistics. Oxford University Press.

Examples

isaaks
summary(isaaks)
plot(isaaks$coords, asp=1, type="n")
text(isaaks$coords, as.character(isaaks$data))
points(isaaks$x0, pch="?", cex=2, col=2)

---

Jitters (duplicated) coordinates.

Description

Jitters 2D coordinates uniformily on a region around (duplicated) points.

Usage

jitter2d(coords, max, min = 0.2 * max, fix.one = TRUE,
         which.fix = c("random", "first", "last")) 

jitterDupCoords(x, ...)

## Default S3 method:
jitterDupCoords(x, ...)

## S3 method for class 'geodata'
jitterDupCoords(x, ...)

Arguments

x, coords	a matrix or data frame with 2D coordinates or geodata object.
max	numeric scalar defining maximum jittering distance.
min	numeric scalar defining minimum jittering distance.
fix.one	logical. Whether or not one of the coordinates should not be jittered.
which.fix	single element vector of integer or character, defining which coordinate won't be jittered. Only used if fix.one=TRUE.
...	arguments passed to jitter2d.

Value

jitter2d returns an object of the same type fo the input with jittered values

jitterDupCoords returns an object of the same type fo the input with jittered coordinate values only at the duplicated locations

Author(s)

Paulo Justiniano Ribeiro Jr. [email protected],
Peter J. Diggle [email protected].

See Also

dup.coords, duplicated.geodata for functions identifying duplicated locations.

Examples

## simulating data
dt <- grf(30, cov.p=c(1, .3)) 
dt$coords <- round(dt$coords, dig=2)
## "forcing" some duplicated locations
dt$coords[4,] <- dt$coords[14,] <- dt$coords[24,] <- dt$coords[2,]
dt$coords[17,] <- dt$coords[23,] <- dt$coords[8,]

## jittering a matrix of duplicated coordinates
dt$coords[c(2,4,14,24),]
jitter2d(dt$coords[c(2,4,14,24),], max=0.01)

## jittering only the duplicated locations and comparing with original
cbind(dt$coords, jitterDupCoords(dt$coords, max=0.01))

## creating a now geodata object jittering the duplicated locations of the original one:
dup.coords(dt)
dt1 <- jitterDupCoords(dt, max=0.01)
dup.coords(dt1)

---

Kattegat basin salinity data

Description

Salinity measurements at the Kattegat basin, Denmark.

Usage

data(kattegat)

Format

An object of the class "geodata", which is list with three components:

coords

the coordinates of the data locations. The distance are given in kilometers.

data

values of the piezometric head. The unit is heads to meters.

dk

a list with cooordinates of lines defining borders and islands across the study area.

Source

National Environmental Research Institute, Arhus University, Denmark and the Swedish Meteorological and Hydrological Institute.

References

Diggle, P. J. and Lophaven, S. (2006). Bayesian geostatistical design. Scandinavian Journal of Statistics, 33: 55-64.

Examples

plot(c(550,770),c(6150,6420),type="n",xlab="X UTM",ylab="Y UTM")
points(kattegat, add=TRUE)
lapply(kattegat$dk, lines, lwd=2)

---

Bayesian Analysis for Gaussian Geostatistical Models

Description

The function krige.bayes performs Bayesian analysis of geostatistical data allowing specifications of different levels of uncertainty in the model parameters.
It returns results on the posterior distributions for the model parameters and on the predictive distributions for prediction locations (if provided).

Usage

krige.bayes(geodata, coords = geodata$coords, data = geodata$data,
            locations = "no", borders, model, prior, output)

model.control(trend.d = "cte", trend.l = "cte", cov.model = "matern",
              kappa = 0.5, aniso.pars = NULL, lambda = 1)

prior.control(beta.prior = c("flat", "normal", "fixed"),
              beta = NULL, beta.var.std = NULL,
              sigmasq.prior = c("reciprocal", "uniform",
                                "sc.inv.chisq", "fixed"),
              sigmasq = NULL, df.sigmasq = NULL,
              phi.prior = c("uniform", "exponential","fixed",
                            "squared.reciprocal", "reciprocal"),
              phi = NULL, phi.discrete = NULL,
              tausq.rel.prior = c("fixed", "uniform", "reciprocal"),
              tausq.rel, tausq.rel.discrete = NULL)

post2prior(obj)

Arguments

geodata	a list containing elements coords and data as described next. Typically an object of the class "geodata" - a geoR data-set. If not provided the arguments coords and data must be provided instead.
coords	an $n \times 2$ matrix where each row has the 2-D coordinates of the $n$ data locations. By default it takes the component coords of the argument geodata, if provided.
data	a vector with n data values. By default it takes the component data of the argument geodata, if provided.
locations	an $N \times 2$ matrix or data-frame with the 2-D coordinates of the $N$ prediction locations, or a list for which the first two components are used. Input is internally checked by the function check.locations. Defaults to "no" in which case the function returns only results on the posterior distributions of the model parameters.
borders	optional. If missing, by default reads the element borders from the geodata object, if present. Setting to NULL preents this behavior. If a two column matrix defining a polygon is provided the prediction is performed only at locations inside this polygon.
model	a list defining the fixed components of the model. It can take an output to a call to model.control or a list with elements as for the arguments in model.control. Default values are assumed for arguments not provided. See section DETAILS below.
prior	a list with the specification of priors for the model parameters. It can take an output to a call to prior.control or a list with elements as for the arguments in prior.control. Default values are assumed for arguments not provided. See section DETAILS below.
output	a list specifying output options. It can take an output to a call to output.control or a list with elements as for the arguments in output.control. Default values are assumed for arguments not provided. See documentation for output.control for further details.
trend.d	specifies the trend (covariates) values at the data locations. See documentation of trend.spatial for further details. Defaults to "cte".
trend.l	specifies the trend (covariates) at the prediction locations. Must be of the same type as defined for trend.d. Only used if prediction locations are provided in the argument locations.
cov.model	string indicating the name of the model for the correlation function. Further details in the documentation for cov.spatial.
kappa	additional smoothness parameter. Only used if the correlation function is one of: "matern", "powered.exponential", "cauchy" or "gneiting.matern". In the current implementation this parameter is always regarded as fixed during the Bayesian analysis.
aniso.pars	fixed parameters for geometric anisotropy correction. If aniso.pars = FALSE no correction is made, otherwise a two elements vector with values for the anisotropy parameters must be provided. Anisotropy correction consists of a transformation of the data and prediction coordinates performed by the function coords.aniso.
lambda	numerical value of the Box-Cox transformation parameter. The value $\lambda = 1$ corresponds to no transformation. The Box-Cox parameter $\lambda$ is always regarded as fixed and data transformation is performed before the analysis. Prediction results are back-transformed and returned is the same scale as for the original data. For $\lambda = 0$ the log-transformation is performed. If $\lambda < 0$ the mean predictor doesn't make sense (the resulting distribution has no expectation).
beta.prior	prior distribution for the mean (vector) parameter $\beta$. The options are "flat" (default), "normal" or "fixed" (known mean).
beta	mean hyperparameter for the distribution of the mean (vector) parameter $\beta$. Only used if beta.prior = "normal" or beta.prior = "fixed". For the later beta defines the value of the known mean.
beta.var.std	standardised (co)variance hyperparameter(s) for the prior for the mean (vector) parameter $\beta$. The (co)variance matrix for$\beta$ is given by the multiplication of this matrix by $\sigma^2$. Only used if beta.prior = "normal".
sigmasq.prior	specifies the prior for the parameter $\sigma^2$. If "reciprocal" (the default), the prior $\frac{1}{\sigma^2}$ is used. Otherwise the parameter is regarded as fixed.
sigmasq	fixed value of the sill parameter $\sigma^2$. Only used if sigmasq.prior = FALSE.
df.sigmasq	numerical. Number of degrees of freedom for the prior for the parameter $\sigma^2$. Only used if sigmasq.prior = "sc.inv.chisq".
phi.prior	prior distribution for the range parameter $\phi$. Options are: "uniform", "exponential", "reciprocal" , "squared.reciprocal" and "fixed". Alternativelly, a user defined discrete distribution can be specified. In this case the argument takes a vector of numerical values of probabilities with corresponding support points provided in the argument phi.discrete. If "fixed" the argument $\phi$ must be provided and is regarded as fixed when performing predictions. For the exponential prior the argument phi must provide the value for of hyperparameter $\nu$ which corresponds to the expected value for this distribution.
phi	fixed value of the range parameter $\phi$. Only needed if phi.prior = "fixed" or if phi.prior = "exponential".
phi.discrete	support points of the discrete prior for the range parameter $\phi$. The default is a sequence of 51 values between 0 and 2 times the maximum distance between the data locations.
tausq.rel.prior	specifies a prior distribution for the relative nugget parameter $\frac{\tau^2}{\sigma^2}$. If tausq.rel.prior = "fixed" the relative nugget is considered known (fixed) with value given by the argument tausq.rel. If tausq.rel.prior = "uniform" a discrete uniform prior is used with support points given by the argument tausq.rel.discrete. Alternativelly, a user defined discrete distribution can be specified. In this case the argument takes the a vector of probabilities of a discrete distribution and the support points should be provided in the argument tausq.rel.discrete.
tausq.rel	fixed value for the relative nugget parameter. Only used if tausq.rel.prior = "fixed".
tausq.rel.discrete	support points of the discrete prior for the relative nugget parameter $\frac{\tau^2}{\sigma^2}$.
obj	an object of the class krige.bayes or posterior.krige.bayes with the output of a call to krige.bayes. The function post2prior takes the posterior distribution computed by one call to krige.bayes and prepares it to be used a a prior in a subsequent call. Notice that in this case the function post2prior is used instead of prior.control.

Details

krige.bayes is a generic function for Bayesian geostatistical analysis of (transformed) Gaussian where predictions take into account the parameter uncertainty.

It can be set to run conventional kriging methods which use known parameters or plug-in estimates. However, the functions krige.conv and ksline are preferable for prediction with fixed parameters.

PRIOR SPECIFICATION

The basis of the Bayesian algorithm is the discretisation of the prior distribution for the parameters $\phi$ and $\tau^2_{rel} =\frac{\tau^2}{\sigma^2}$. The Tech. Report (see References below) provides details on the results used in the current implementation.
The expressions of the implemented priors for the parameter $\phi$ are:

"uniform":

$p(\phi) \propto 1$.

"exponential":

$p(\phi) = \frac{1}{\nu} \exp(-\frac{1}{\nu} * \phi)$.

"reciprocal":

$p(\phi) \propto \frac{1}{\phi}$.

"squared.reciprocal":

$p(\phi) \propto \frac{1}{\phi^2}$.

"fixed":

fixed known or estimated value of $\phi$.

The expressions of the implemented priors for the parameter $\tau^2_{rel}$ are:

"fixed":

fixed known or estimated value of $\tau^2_{rel}$. Defaults to zero.

"uniform":

$p(\tau^2_{rel}) \propto 1$.

"reciprocal":

$p(\tau^2_{rel}) \propto \frac{1}{\tau^2_{rel}}$.

Apart from those a user defined prior can be specifyed by entering a vector of probabilities for a discrete distribution with suport points given by the argument phi.discrete and/or tausq.rel.discrete.

CONTROL FUNCTIONS

The function call includes auxiliary control functions which allows the user to specify and/or change the specification of model components (using model.control), prior distributions (using prior.control) and output options (using output.control). Default options are available in most of the cases.

Value

An object with class "krige.bayes" and "kriging". The attribute prediction.locations containing the name of the object with the coordinates of the prediction locations (argument locations) is assigned to the object. Returns a list with the following components:

posterior

results on on the posterior distribution of the model parameters. A list with the following possible components:

beta

summary information on the posterior distribution of the mean parameter $\beta$.

sigmasq

summary information on the posterior distribution of the variance parameter $\sigma^2$ (partial sill).

phi

summary information on the posterior distribution of the correlation parameter $\phi$ (range parameter) .

tausq.rel

summary information on the posterior distribution of the relative nugget variance parameter $\tau^2_{rel}$.

joint.phi.tausq.rel

information on discrete the joint distribution of these parameters.

sample

a data.frame with a sample from the posterior distribution. Each column corresponds to one of the basic model parameters.

predictive

results on the predictive distribution at the prediction locations, if provided. A list with the following possible components:

mean

expected values.

variance

expected variance.

distribution

type of posterior distribution.

mean.simulations

mean of the simulations at each locations.

variance.simulations

variance of the simulations at each locations.

quantiles.simulations

quantiles computed from the the simulations.

probabilities.simulations

probabilities computed from the simulations.

simulations

simulations from the predictive distribution.

prior	a list with information on the prior distribution and hyper-parameters of the model parameters ($\beta, \sigma^2, \phi, \tau^2_{rel}$).
model	model specification as defined by model.control.
.Random.seed	system random seed before running the function. Allows reproduction of results. If the .Random.seed is set to this value and the function is run again, it will produce exactly the same results.
max.dist	maximum distance found between two data locations.
call	the function call.

Author(s)

Paulo J. Ribeiro Jr. [email protected],
Peter J. Diggle [email protected].

References

Diggle, P.J. & Ribeiro Jr, P.J. (2002) Bayesian inference in Gaussian model-based geostatistics. Geographical and Environmental Modelling, Vol. 6, No. 2, 129-146.

The technical details about the implementation of krige.bayes can be found at:
Ribeiro, P.J. Jr. and Diggle, P.J. (1999) Bayesian inference in Gaussian model-based geostatistics. Tech. Report ST-99-08, Dept Maths and Stats, Lancaster University.
Available at: http://www.leg.ufpr.br/geoR/geoRdoc/bayeskrige.pdf

Further information about geoR can be found at:
http://www.leg.ufpr.br/geoR/.

For a extended list of examples of the usage see http://www.leg.ufpr.br/geoR/tutorials/examples.krige.bayes.R and/or the geoR tutorials page at http://www.leg.ufpr.br/geoR/tutorials/.

See Also

lines.variomodel.krige.bayes, plot.krige.bayes for outputs related to the parameters in the model, image.krige.bayes and persp.krige.bayes for graphical output of prediction results. krige.conv and ksline for conventional kriging methods.

Examples

## Not run: 
# generating a simulated data-set
ex.data <- grf(70, cov.pars=c(10, .15), cov.model="matern", kappa=2)
#
# defining the grid of prediction locations:
ex.grid <- as.matrix(expand.grid(seq(0,1,l=21), seq(0,1,l=21)))
#
# computing posterior and predictive distributions
# (warning: the next command can be time demanding)
ex.bayes <- krige.bayes(ex.data, loc=ex.grid,
                 model = model.control(cov.m="matern", kappa=2),
                 prior = prior.control(phi.discrete=seq(0, 0.7, l=51),
                             phi.prior="reciprocal"))
#
# Prior and posterior for the parameter phi
plot(ex.bayes, type="h", tausq.rel = FALSE, col=c("red", "blue"))
#
# Plot histograms with samples from the posterior
par(mfrow=c(3,1))
hist(ex.bayes)
par(mfrow=c(1,1))

# Plotting empirical variograms and some Bayesian estimates:
# Empirical variogram
plot(variog(ex.data, max.dist = 1), ylim=c(0, 15))
# Since ex.data is a simulated data we can plot the line with the "true" model 
lines.variomodel(ex.data, lwd=2)
# adding lines with summaries of the posterior of the binned variogram
lines(ex.bayes, summ = mean, lwd=1, lty=2)
lines(ex.bayes, summ = median, lwd=2, lty=2)
# adding line with summary of the posterior of the parameters
lines(ex.bayes, summary = "mode", post = "parameters")

# Plotting again the empirical variogram
plot(variog(ex.data, max.dist=1), ylim=c(0, 15))
# and adding lines with median and quantiles estimates
my.summary <- function(x){quantile(x, prob = c(0.05, 0.5, 0.95))}
lines(ex.bayes, summ = my.summary, ty="l", lty=c(2,1,2), col=1)

# Plotting some prediction results
op <- par(no.readonly = TRUE)
par(mfrow=c(2,2), mar=c(4,4,2.5,0.5), mgp = c(2,1,0))
image(ex.bayes, main="predicted values")
image(ex.bayes, val="variance", main="prediction variance")
image(ex.bayes, val= "simulation", number.col=1,
      main="a simulation from the \npredictive distribution")
image(ex.bayes, val= "simulation", number.col=2,
      main="another simulation from \nthe predictive distribution")
#
par(op)

## End(Not run)
##
## For a extended list of exemples of the usage of krige.bayes()
## see http://www.leg.ufpr.br/geoR/tutorials/examples.krige.bayes.R
##

---