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

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 p.brown@lancaster.ac.uk.
Adapted and included in geoR by
Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br

References

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

See Also

nlm, optim.

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

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

Saturated Hydraulic Conductivity

Description

The data consists of 32 measurements of the saturated hydraulic conductivity of a soil.

Usage

data(Ksat)

Format

The object Ksat is a list of the class geodata with the following elements:

coords

a matrix with the coordinates of the soil samples.

data

measurements of the saturated hidraulic conductivity.

borders

a data-frame with the coordinates of a polygon defining the borders of the area.

Source

Data provided by Dr. Décio Cruciani, ESALQ/USP, Brasil.

Examples

summary(Ksat)
plot(Ksat)

Spatial Interpolation Comparison data

Description

Data from the SIC-97 project: Spatial Interpolation Comparison.

Usage

data(SIC)

Format

Four objects of the class "geodata": sic.all, sic.100, sic.367, sic.some. Each is a list with two components:

coords

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

data

rainfall values. The unit is milimeters.

altitude

elevation values. The unit is milimeters.

Additionally an matrix sic.borders with the borders of the country.

Source

Data from the project Spatial Interpolation Comparison 97; see https://soft.mines-paristech.fr/gstlearn/latest/data/benchmark/SIC97_Readme.pdf.

References

Christensen, O.F., Diggle, P.J. and Ribeiro Jr, P.J. (2001) Analysing positive-valued spatial data: the transformed Gaussian model. In Monestiez, P., Allard, D. and Froidevaux (eds), GeoENV III - Geostatistics for environmental applications. Quantitative Geology and Geostatistics, Kluwer Series, 11, 287–298.

Examples

points(sic.100, borders=sic.borders)
points(sic.all, borders=sic.borders)

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

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:

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

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

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:

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

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

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

Value

The functions returns the following results:

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

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

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. paulojus@leg.ufpr.br
Peter J. Diggle p.diggle@lancaster.ac.uk.

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

Value

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

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

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. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

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

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. paulojus@leg.ufpr.br
Peter J. Diggle p.diggle@lancaster.ac.uk.

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

Value

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

Author(s)

Andreas Kiefer andreas@inf.ufpr.br
Paulo Justiniano Ribeiro Junior paulojus@leg.ufpr.br.

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

geoR internal functions

Description

These are functions internally called by other functions in the package geoR and not meant to be called by the user.

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

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

Value

An scalar with the value of the global variance

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

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

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:

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

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

Value

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

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

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

Value

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

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

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

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. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

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

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. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

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)

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

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. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

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

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:

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.

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.

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

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

Spatial Prediction – Conventional Kriging

Description

This function performs spatial prediction for fixed covariance parameters using global neighbourhood.

Options available implement the following types of kriging: SK (simple kriging), OK (ordinary kriging), KTE (external trend kriging) and UK (universal kriging).

Usage

krige.conv(geodata, coords=geodata$coords, data=geodata$data,
           locations, borders, krige, output)

krige.control(type.krige = "ok", trend.d = "cte", trend.l = "cte",
            obj.model = NULL, beta, cov.model, cov.pars, kappa,
            nugget, micro.scale = 0, dist.epsilon = 1e-10, 
            aniso.pars, lambda)

Arguments

Details

According to the arguments provided, one of the following different types of kriging: SK, OK, UK or KTE is performed. Defaults correspond to ordinary kriging.

Value

An object of the class 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:

Other results can be included depending on the options passed to output.control.

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

output.control sets output options, image.kriging and persp.kriging for graphical output of the results, krige.bayes for Bayesian prediction and ksline for a different implementation of kriging allowing for moving neighborhood. For model fitting see likfit or variofit.

Examples

## Not run: 
# Defining a prediction grid
loci <- expand.grid(seq(0,1,l=21), seq(0,1,l=21))
# predicting by ordinary kriging
kc <- krige.conv(s100, loc=loci,
                 krige=krige.control(cov.pars=c(1, .25)))
# mapping point estimates and variances
par.ori <- par(no.readonly = TRUE)
par(mfrow=c(1,2), mar=c(3.5,3.5,1,0), mgp=c(1.5,.5,0))
image(kc, main="kriging estimates")
image(kc, val=sqrt(kc$krige.var), main="kriging std. errors")
# Now setting the output to simulate from the predictive
# (obtaining conditional simulations),
# and to compute quantile and probability estimators
s.out <- output.control(n.predictive = 1000, quant=0.9, thres=2)
set.seed(123)
kc <- krige.conv(s100, loc = loci,
         krige = krige.control(cov.pars = c(1,.25)),
         output = s.out)
par(mfrow=c(2,2))
image(kc, val=kc$simul[,1], main="a cond. simul.")
image(kc, val=kc$simul[,1], main="another cond. simul.")
image(kc, val=(1 - kc$prob), main="Map of P(Y > 2)")
image(kc, val=kc$quant, main="Map of y s.t. P(Y < y) = 0.9")
par(par.ori)

## End(Not run)

Computes kriging weights

Description

Computes the weights assign for each data point in simple and ordinary krigring

Usage

krweights(coords, locations, krige)

Arguments

Value

A matrix of weights

Examples

## Figure 8.4 in Webster and Oliver (2001), see help(wo)
attach(wo)
par(mfrow=c(2,2))
plot(c(-10,130), c(-10,130), ty="n", asp=1)
points(rbind(coords, x1))
KC1 <- krige.control(cov.pars=c(0.382,90.53))
w1 <- krweights(wo$coords, loc=x1, krige=KC1)
text(coords[,1], 5+coords[,2], round(w1, dig=3))
##
plot(c(-10,130), c(-10,130), ty="n", asp=1)
points(rbind(coords, x1))
KC2 <- krige.control(cov.pars=c(0.282,90.53), nug=0.1)
w2 <- krweights(wo$coords, loc=x1, krige=KC2)
text(coords[,1], 5+coords[,2], round(w2, dig=3))
##
plot(c(-10,130), c(-10,130), ty="n", asp=1)
points(rbind(coords, x1))
KC3 <- krige.control(cov.pars=c(0.082,90.53), nug=0.3)
w3 <- krweights(wo$coords, loc=x1, krige=KC3)
text(coords[,1], 5+coords[,2], round(w3, dig=3))
##
plot(c(-10,130), c(-10,130), ty="n", asp=1)
points(rbind(coords, x1))
KC4 <- krige.control(cov.pars=c(0,90.53), nug=0.382, micro=0.382)
w4 <- krweights(wo$coords, loc=x1, krige=KC4)
text(coords[,1], 5+coords[,2], round(w4, dig=3))
##
## SK vs OK
##
plot(c(-10,130), c(-10,130), ty="n", asp=1)
points(rbind(coords, x1))
KC5 <- krige.control(cov.pars=c(0.382,50))
w5 <- krweights(wo$coords, loc=x1, krige=KC5)
KC6 <- krige.control(type="sk", beta=2, cov.pars=c(0.382,50))
w6 <- krweights(wo$coords, loc=x1, krige=KC6)
text(coords[,1],  5+coords[,2], round(w5, dig=3))
text(coords[,1], -5+coords[,2], round(w6, dig=3))
##
plot(c(-10,130), c(-10,130), ty="n", asp=1)
points(rbind(coords, x1))
KC7 <- krige.control(cov.pars=c(0.382,0))
w7 <- krweights(wo$coords, loc=x1, krige=KC7)
KC8 <- krige.control(type="sk", beta=2, cov.pars=c(0.382,0))
w8 <- krweights(wo$coords, loc=x1, krige=KC8)
text(coords[,1],  5+coords[,2], round(w7, dig=3))
text(coords[,1], -5+coords[,2], round(w8, dig=3))

Spatial Prediction – Conventional Kriging

Description

This function performs spatial prediction for given covariance parameters. Options implement the following kriging types: SK (simple kriging), OK (ordinary kriging), KTE (external trend kriging) and UK (universal kriging).

The function krige.conv should be preferred, unless moving neighborhood is to be used.

Usage

ksline(geodata, coords = geodata$coords, data = geodata$data,
       locations, borders = NULL, 
       cov.model = "matern",
       cov.pars=stop("covariance parameters (sigmasq and phi) needed"),
       kappa = 0.5, nugget = 0, micro.scale = 0,
       lambda = 1, m0 = "ok", nwin = "full",
       n.samples.backtransform = 500, trend = 1, d = 2,
       ktedata = NULL, ktelocations = NULL, aniso.pars = NULL,
       signal = FALSE, dist.epsilon = 1e-10, messages)

Arguments

Value

An object of the class kriging which is a list with the following components:

Note

This is a preliminary and inefficient function implementing kriging methods. For predictions using global neighborhood the function krige.conv should be used instead.

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

krige.conv for a more efficient implementation of conventional kriging methods,
krige.bayes for Bayesian prediction.

Examples

loci <- expand.grid(seq(0,1,l=31), seq(0,1,l=31))
kc <- ksline(s100, loc=loci, cov.pars=c(1, .25))
par(mfrow=c(1,2))
image(kc, main="kriging estimates")
image(kc, val=sqrt(kc$krige.var), main="kriging std. errors")

Data from Landim's book

Description

Artificial or non-specified data from Paulo Landim's book

Usage

data(landim1)

Format

A data frame with 38 observations on the following 4 variables.

EW

a numeric vector with the east-west coordinates.

NS

a numeric vector with the north-south coordinates.

A

a numeric vector with data on a first variable.

B

a numeric vector with data on a second variable.

Source

Landim, P. M. B. (2004) Análise estatística de dados geológicos. Editora Unesp. Data from Table~1, pg.12.

Examples

data(landim)
plot(as.geodata(landim1, data.col=3))
plot(as.geodata(landim1, data.col=4))

Add a legend to a image with kriging results

Description

This function allows adds a legend to an image plot generated by image.kriging or image.krige.bayes. It can be called internally by these functions or directly by the user.

Usage

legend.krige(x.leg, y.leg, values, scale.vals,
             vertical = FALSE, offset.leg = 1, ...)

Arguments

Value

A legend is added to the current plot. No values are returned.

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

image.kriging, image.krige.bayes.

Examples

# See examples in the documentation for image.kriging

Likelihood Based Parameter Estimation for Gaussian Random Fields

Description

Maximum likelihood (ML) or restricted maximum likelihood (REML) parameter estimation for (transformed) Gaussian random fields.

Usage

likfit(geodata, coords = geodata$coords, data = geodata$data,
       trend = "cte", ini.cov.pars, fix.nugget = FALSE, nugget = 0,
       fix.kappa = TRUE, kappa = 0.5, fix.lambda = TRUE, lambda = 1,
       fix.psiA = TRUE, psiA = 0, fix.psiR = TRUE, psiR = 1, 
       cov.model, realisations, lik.method = "ML", components = TRUE,
       nospatial = TRUE, limits = pars.limits(),
       print.pars = FALSE, messages, ...)

## S3 method for class 'likGRF'
fitted(object, spatial = TRUE, ...)

## S3 method for class 'likGRF'
resid(object, spatial = FALSE, ...)

Arguments

Details

This function estimate the parameters of the Gaussian random field model, specified as:

Y(x) = \mu(x) + S(x) + e

where

• x defines a spatial location. Typically Euclidean coordinates on a plane.

• Y is the variable been observed.

• \mu(x) = X\beta is the mean component of the model (trend).

• S(x) is a stationary Gaussian process with variance \sigma^2 (partial sill) and a correlation function parametrized in its simplest form by \phi (the range parameter). Possible extra parameters for the correlation function are the smoothness parameter \kappa and the anisotropy parameters \phi_R and \phi_A (anisotropy ratio and angle, respectively).

• e is the error term with variance parameter \tau^2 (nugget variance).

The additional parameter \lambda allows for the Box-Cox transformation of the response variable. If used (i.e. if \lambda \neq 1) Y(x) above is replaced by g(Y(x)) such that

g(Y(x)) = \frac{Y^\lambda(x) - 1}{\lambda}.

Two particular cases are \lambda = 1 which indicates no transformation and \lambda = 0 indicating the log-transformation.

Numerical minimization

In general parameter estimation is performed numerically using the R function optim to minimise the negative log-likelihood computed by the function negloglik.GRF. If the nugget, anisotropy (\psi_A, \psi_R), smoothness (\kappa) and transformation (\lambda) parameters are held fixed then the numerical minimisation can be reduced to one-dimension and the function optimize is used instead of optim. In this case initial values are irrelevant.

Limits

Lower and upper limits for parameter values can be individually specified using the function link{pars.limits}. For example, including the following in the function call:
limits = pars.limits(phi=c(0, 10), lambda=c(-2.5, 2.5)),
will change the limits for the parameters \phi and \lambda. Default values are used if the argument limits is not provided.

There are internal reparametrisation depending on the options for parameters to be estimated. For instance for the common situation when fix.nugget=FALSE the minimisation is performed in a reduced parameter space using \tau^2_{rel} = \frac{\tau^2}{\sigma^2}. In this case values of \sigma^2 and \beta are then given by analytical expressions which are function of the two parameters remaining parameters and limits for these two parameters will be ignored.

Since parameter values are found by numerical optimization using the function optim, in given circunstances the algorithm may not converge to correct parameter values when called with default options and the user may need to pass extra options for the optimizer. For instance the function optim takes a control argument. The user should try different initial values and if the parameters have different orders of magnitude may need to use options to scale the parameters. Some possible workarounds in case of problems include:

• rescale you data values (dividing by a constant, say)

• rescale your coordinates (subtracting values and/or dividing by constants)

• Use the mechanism to pass control() options for the optimiser internally

Transformation If the fix.lambda = FALSE and nospatial = FALSE the Box-Cox parameter for the model without the spatial component is obtained numerically, with log-likelihood computed by the function boxcox.ns.

Multiple initial values can be specified providing a n \times 2 matrix for the argument ini.cov.pars and/or providing a vector for the values of the remaining model parameters. In this case the log-likelihood is computed for all combinations of the model parameters. The parameter set which maximises the value of the log-likelihood is then used to start the minimisation algorithm.

Alternatively the argument ini.cov.pars can take an object of the class eyefit or variomodel. This allows the usage of an output of the functions eyefit, variofit or likfit be used as initial value.

The argument realisations allows sets of data assumed to be independent replications of the same process. Data on different realisations may or may not be co-located. For instance, data collected at different times can be pooled together in the parameter estimation assuming time independence. The argument realisations takes a vector indicating the replication number (e.g. the times). If realisations = TRUE the code looks for an element named realisations in the geodata object. The log-likelihoods are computed for each replication and added together.

Value

An object of the classes "likGRF" and "variomodel".
The function summary.likGRF is used to print a summary of the fitted model.
The object is a list with the following components:

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

summary.likGRF for summary of the results, plot.variogram, lines.variogram and lines.variomodel for graphical output, proflik for computing profile likelihoods, variofit and for other estimation methods, and optim for the numerical minimisation function.

Examples

## Not run: 
ml <- likfit(s100, ini=c(0.5, 0.5), fix.nug = TRUE)
ml
summary(ml)
reml <- likfit(s100, ini=c(0.5, 0.5), fix.nug = TRUE, lik.met = "REML")
summary(reml)
plot(variog(s100))
lines(ml)
lines(reml, lty = 2)

## End(Not run)

Fits the bivariate Gaussian common component geostatistical model

Description

Computes maximum likelihood estimates of the bivariate Gaussian common component geostatistical model.

Usage

likfitBGCCM(geodata1, geodata2, ini.sigmasq, ini.phi,
            cov0.model="matern", cov1.model="matern", cov2.model="matern",
            kappa0=0.5, kappa1=0.5, kappa2=0.5,
            fc.min = c("optim", "nlminb"), ...)

Arguments

Value

A list with model fitting information to which the class BGCCM is assigned.

Warning

This is a new function and still in draft format and pretty much untested.

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

See Also

optim, nlminb, varcovBGCCM, as.geodata, likfit.

Examples

# see http://www.leg.ufpr.br/geoR/tutorials/CCM.R

Line with a Empirical Variogram

Description

A sample (empirical) variogram computed using the function variog is added to the current plot.

Usage

## S3 method for class 'variogram'
lines(x, max.dist, type = "o",  scaled = FALSE,
         pts.range.cex, ...)

Arguments

Value

A line with the empirical variogram is added to the plot in the current graphics device. No values are returned.

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

variog, lines.variogram, lines.variomodel, variog.model.env, variog.mc.env, plot.grf, lines.

Adds Envelopes Lines to a Variogram Plot

Description

Variogram envelopes computed by variog.model.env or variog.mc.env are added to the current variogram plot.

Usage

## S3 method for class 'variogram.envelope'
lines(x, lty = 3, ...)

Arguments

Value

Lines defining the variogram envelope are added to the plotin the current graphics device.

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

variog for variogram computation, variog.model.env and variog.mc.env for computation of variogram envelopes, and lines for the generic function.

Examples

s100.vario <- variog(s100, max.dist = 1)
s100.ml <- likfit(s100, ini=c(.5, .5))
s100.mod.env <- variog.model.env(s100, obj.variog = s100.vario,
   model = s100.ml) 
s100.mc.env <- variog.mc.env(s100, obj.variog = s100.vario)
plot(s100.vario)
lines(s100.mod.env)
lines(s100.mc.env, lwd=2)

Adds a Line with a Variogram Model to a Variogram Plot

Description

This function adds a line with a variogram model specifyed by the user to a current variogram plot. The variogram is specifyed either by passing a list with values for the variogram elements or using each argument in the function.

Usage

## S3 method for class 'variomodel'
lines(x, ...)
## Default S3 method:
lines.variomodel(x, cov.model, cov.pars, nugget, kappa,
                          max.dist, scaled = FALSE, ...)

Arguments

Details

Adds a line with a variogram model to a plot. In conjuction with plot.variogram can be used for instance to compare sample variograms against fitted models returned by variofit and/or likfit.

Value

A line with a variogram model is added to a plot on the current graphics device. No values are returned.

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

lines.variomodel.krige.bayes, lines.variomodel.grf, lines.variomodel.variofit, lines.variomodel.likGRF, plot.variogram, lines.variogram, variofit, likfit, curve.

Examples

# computing and ploting empirical variogram
vario <- variog(s100, max.dist = 1)
plot(vario)
# estimating parameters by weighted least squares
vario.wls <- variofit(vario, ini = c(1, .3), fix.nugget = TRUE)
# adding fitted model to the plot  
lines(vario.wls)
#
# Ploting different variogram models
plot(0:1, 0:1, type="n")
lines.variomodel(cov.model = "exp", cov.pars = c(.7, .25), nug = 0.3, max.dist = 1) 
# an alternative way to do this is:
my.model <- list(cov.model = "exp", cov.pars = c(.7, .25), nugget = 0.3,
max.dist = 1) 
lines.variomodel(my.model, lwd = 2)
# now adding another model
lines.variomodel(cov.m = "mat", cov.p = c(.7, .25), nug = 0.3,
                 max.dist = 1, kappa = 1, lty = 2)
# adding the so-called "nested" models
# two exponential structures
lines.variomodel(seq(0,1,l=101), cov.model="exp",
                 cov.pars=rbind(c(0.6,0.15),c(0.4,0.25)), nug=0, col=2)
## exponential and spherical structures
lines.variomodel(seq(0,1,l=101), cov.model=c("exp", "sph"),
                 cov.pars=rbind(c(0.6,0.15), c(0.4,0.75)), nug=0, col=3)

Lines with True Variogram for Simulated Data

Description

This functions adds to the graphics device a line with the theoretical (true) variogram used when generating simulations with the function grf.

Usage

## S3 method for class 'grf'
lines.variomodel(x, max.dist, n = 100, ...)

Arguments

Value

A line with the true variogram model is added to the current plot on the graphics device. No values are returned.

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

lines.variomodel, grf, plot.grf, curve.

Examples

sim <- grf(100, cov.pars=c(1, .25)) # simulates data
plot(variog(sim, max.dist=1))       # plot empirical variogram

Adds a Bayesian Estimate of the Variogram to a Plot

Description

Adds a Bayesian estimate of the variogram model to a plot typically with an empirical variogram. The estimate is a chosen summary (mean, mode or mean) of the posterior distribution returned by the function krige.bayes.

Usage

## S3 method for class 'krige.bayes'
lines.variomodel(x, summary.posterior, max.dist, uvec,
                 posterior = c("variogram", "parameters"),  ...)

Arguments

Details

The function krige.bayes returns samples from the posterior distribution of the parameters (\sigma^2, \phi, \tau^{2}_{rel}).

This function allows for two basic options to draw a line with a summary of the variogram function.

"variogram":

for each sample of the parameters the variogram function is computed at the support points defined in the argument uvec. Then a function provided by the user in the argument summary.posterior is used to compute a summary of the values obtained at each support point.

"parameters":

in this case summaries of the posterior distribution of the model parameters as “plugged-in” in the variogram function. One of the options "mode" (default) ,"median" or "mean" can be provided in the argument summary.posterior. The option mode, uses the mode of (\phi, \tau^{2}_{rel}) and the mode of of \sigma^2 conditional on the modes of the former parameters. For the options mean and median these summaries are computed from the samples of the posterior.

Value

A line with the estimated variogram plot is added to the plot in the current graphics device. No values are returned.

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

lines.variomodel, krige.bayes and lines.

Examples

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

Adds a Variogram Line to a Variogram Plot

Description

This function adds a fitted variogram based on the estimates of the model parameters returned by the function likfit to a current variogram plot.

Usage

## S3 method for class 'likGRF'
lines.variomodel(x, max.dist, scaled = FALSE, ...)

Arguments

Details

Adds variogram model(s) to a plot. In conjuction with plot.variogram can be used to compare sample variograms against fitted models returned by likfit.

Value

A line with a variogram model is added to a plot on the current graphics device. No values are returned.

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

lines.variomodel, lines.variomodel.variofit, plot.variogram, lines.variogram, variofit, likfit, curve.

Examples

# compute and plot empirical variogram
vario <- variog(s100, max.dist = 1)
plot(vario)
# estimate parameters
vario.ml <- likfit(s100, ini = c(1, .3), fix.nugget = TRUE)
# adds fitted model to the plot  
lines(vario.ml)

Adds a Line with a Fitted Variogram Model to a Variogram Plot

Description

This function adds a line with the variogram model fitted by the function variofit to a current variogram plot.

Usage

## S3 method for class 'variofit'
lines.variomodel(x, max.dist, scaled = FALSE, ...)

Arguments

Details

Adds fitted variogram model to a plot. In conjuction with plot.variogram can be used to compare empirical variograms against fitted models returned by variofit.

Value

A line with a variogram model is added to a plot on the current graphics device. No values are returned.

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

lines.variomodel, lines.variomodel.likGRF, plot.variogram, lines.variogram, variofit, likfit, curve.

Examples

# compute and plot empirical variogram
vario <- variog(s100, max.dist = 1)
plot(vario)
# estimate parameters
vario.wls <- variofit(vario, ini = c(1, .3), fix.nugget = TRUE)
# adds fitted model to the plot  
lines(vario.wls)

Select prediction locations inside borders

Description

Selects the prediction locations located inside a polygon defining borders of a region where prediction is aimed. Typically internally called by geoR functions krige.bayes, krige.conv, ksline.

Usage

locations.inside(locations, borders, as.is = TRUE, ...)

Arguments

Value

A two columns matrix, data-frame or a list with 2 elements with coordinates of points inside the borders.

See Also

over,coordinates, SpatialPoints.

Examples

gr <- pred_grid(parana$borders, by=20)
plot(gr, asp=1, pch="+")
polygon(parana$borders)
gr.in <- locations.inside(gr, parana$borders)
points(gr.in, col=2, pch="+")

Log-Likelihood for a Gaussian Random Field

Description

This function computes the value of the log-likelihood for a Gaussian random field.

Usage

loglik.GRF(geodata, coords = geodata$coords, data = geodata$data,
           obj.model = NULL, cov.model = "exp", cov.pars, nugget = 0,
           kappa = 0.5, lambda = 1, psiR = 1, psiA = 0,
           trend = "cte", method.lik = "ML", compute.dists = TRUE,
           realisations = NULL)

Arguments

Details

The expression log-likelihood is:

l(\theta) = -\frac{n}{2} \log (2\pi) - \frac{1}{2} \log |\Sigma| - \frac{1}{2} (y - F\beta)' \Sigma^{-1} (y - F\beta),

where n is the size of the data vector y, \beta is the mean (vector) parameter with dimention p, \Sigma is the covariance matrix and F is the matrix with the values of the covariates (a vector of 1's if the mean is constant.

The expression restricted log-likelihood is:

rl(\theta) = -\frac{n-p}{2} \log (2\pi) + \frac{1}{2} \log |F' F| - \frac{1}{2} \log |\Sigma| - \frac{1}{2} \log |F' \Sigma F| - \frac{1}{2} (y - F\beta)' \Sigma^{-1} (y - F\beta).

Value

The numerical value of the log-likelihood.

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

likfit for likelihood-based parameter estimation.

Examples

loglik.GRF(s100, cov.pars=c(0.8, .25), nugget=0.2)
loglik.GRF(s100, cov.pars=c(0.8, .25), nugget=0.2, met="RML")

## Computing the likelihood of a model fitted by ML
s100.ml <- likfit(s100, ini=c(1, .5))
s100.ml
s100.ml$loglik
loglik.GRF(s100, obj=s100.ml)

## Computing the likelihood of a variogram fitted model
s100.v <- variog(s100, max.dist=1)
s100.vf <- variofit(s100.v, ini=c(1, .5))
s100.vf
loglik.GRF(s100, obj=s100.vf)

Computer Values of the Matern Correlation Function

Description

This function computes values of the \mbox{Mat\'{e}rn} correlation function for given distances and parameters.

Usage

matern(u, phi, kappa)

Arguments

Details

The \mbox{Mat\'{e}rn} model is defined as:

\rho(u;\phi,\kappa) =\{2^{\kappa-1} \Gamma(\kappa)\}^{-1} (u/\phi)^\kappa K_\kappa(u/\phi)

where \phi and \kappa are parameters and K_\kappa(\cdot) denotes the modified Bessel function of the third kind of order \kappa. The family is valid for \phi>0 and \kappa>0.

Value

A vector matrix or array, according to the argument u, with the values of the \mbox{Mat\'{e}rn} correlation function for the given distances.

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

See Also

cov.spatial for the correlation functions implemented in geoR, and besselK for computation of the Bessel functions.

Examples

#
# Models with fixed range and varying smoothness parameter
#
curve(matern(x, phi= 0.25, kappa = 0.5),from = 0, to = 1.5,
      xlab = "distance", ylab = expression(rho(h)), lty = 2,
      main=expression(paste("varying  ", kappa, "  and fixed  ", phi)))
curve(matern(x, phi= 0.25, kappa = 1),from = 0, to = 1.5, add = TRUE)
curve(matern(x, phi= 0.25, kappa = 2),from = 0, to = 1.5, add = TRUE,
      lwd = 2, lty=2)
curve(matern(x, phi= 0.25, kappa = 3),from = 0, to = 1.5, add = TRUE,
      lwd = 2)
legend("topright", expression(kappa==0.5, kappa==1.5, kappa==2, kappa==3),
    lty=c(2,1,2,1), lwd=c(1,1,2,2))

#
# Correlations with equivalent "practical range"
# and varying smoothness parameter
#
curve(matern(x, phi = 0.25, kappa = 0.5),from = 0, to = 1,
      xlab = "distance", ylab = expression(gamma(h)), lty = 2,
      main = "models with equivalent \"practical\" range")
curve(matern(x, phi = 0.188, kappa = 1),from = 0, to = 1, add = TRUE)      
curve(matern(x, phi = 0.14, kappa = 2),from = 0, to = 1,
      add = TRUE, lwd=2, lty=2)      
curve(matern(x, phi = 0.117, kappa = 2), from = 0, to = 1,
      add = TRUE, lwd=2)      
legend("topright", 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))

Lists names of the key elements of a geodata object

Description

Produces a list with the names of the main elements of geodata object: coords, data, units.m, covariate and realisation. Can be useful to list names before using {subset.geodata}.

Usage

## S3 method for class 'geodata'
names(x)

Arguments

Value

A list with

See Also

names, subset.geodata, as.geodata.

Examples

names(ca20)

Near location to a point

Description

For a given set of points and locations identified by 2D coordinates this function finds the nearest location of each point

Usage

nearloc(points, locations, positions = FALSE)

Arguments

Value

If positions = FALSE the function returns a matrix, data-frame or list of the same type and size as the object provided in the argument points with the coordinates of the nearest locations.

If positions = FALSE the function returns a vector with the position of the nearest points in the locations object.

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

See Also

loccoords

Examples

set.seed(276)
gr <- expand.grid(seq(0,1, l=11), seq(0,1, l=11))
plot(gr, asp=1)
pts <- matrix(runif(10), nc=2)
points(pts, pch=19)
near <- nearloc(points=pts, locations=gr)
points(near, pch=19, col=2)
rownames(near)
nearloc(points=pts, locations=gr, pos=TRUE)

Defines output options for prediction functions

Description

Auxiliary function defining output options for krige.bayes and krige.conv.

Usage

output.control(n.posterior, n.predictive, moments, n.back.moments,
               simulations.predictive, mean.var, quantile,
               threshold, sim.means, sim.vars, signal, messages)

Arguments

Details

SIGNAL

This function is typically called by the geoR's prediction functions krige.bayes and krige.conv defining the output.

By default, krige.bayes sets signal = TRUE and krige.conv sets signal = FALSE.

The underlying model

Y(x) = \mu + S(x) + \epsilon

assumes that observations Y(x) are noisy versions of a signal S(x) and Var(\epsilon)=\tau^2 is the nugget variance.

If \tau^2 = 0 the Y and S are indistiguishable.

If \tau^2 > 0 and regarded as measurement error, the option signal defines whether the S (signal = TRUE) or the variable Y (signal = FALSE) is to be predicted.
For the latter the predictions will "honor" the data, i.e. predicted values will coincide with the data, at data locations.
For unsampled locations and untransformed data, the predicted values equals data regardless signal = TRUE or FALSE, however predictions variances will differ.

The function krige.conv has an argument micro.scale. If micro.scale > 0 the error term is divided as \epsilon = \epsilon_{ms} + \epsilon_{me} and the nugget variance is divided into two terms: micro-scale variance and measurement error.
If signal = TRUE the term \epsilon_{ms} is regarded as part of the signal and consequently the micro-scale variance is added to the prediction variance.
If signal = FALSE the total error variance \tau^2 is added to the prediction variance.

Value

A list with processed arguments to be passed to the main function.

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

See Also

The prediction functions krige.bayes and krige.conv.

Rainfall Data from Parana State, Brasil

Description

This data-set was used by Diggle and Ribeiro (2001) to illustrate the methods discussed in the paper. The data reported analysis was carried out using the package geoR.

The data refers to average rainfall over different years for the period May-June (dry-season). It was collected at 143 recording stations throughout \mbox{Paran\'{a}} State, Brasil.

Usage

data(parana)

Format

The object parana of the class geodata, which is a list containing the following components:

coords

a matrix with the coordinates of the recording stations.

data

a vector with the average recorded rainfall for the May-June period.

borders

a matrix with the coordinates defining the borders of \mbox{Paran\'{a}} state.

loci.paper

a matrix with the coordinates of the four prediction locations discussed in the paper.

Source

The data were collected at several recording stations at \mbox{Paran\'{a}} State, Brasil, belonging to the following companies: COPEL, IAPAR, DNAEE, SUREHMA and INEMET.

The data base was organized by Laura Regina Bernardes Kiihl (IAPAR, Instituto \mbox{Agron\^{o}mico} do \mbox{Paran\'{a}}, Londrina, Brasil) and the fraction of the data included in this data-set was provided by Jacinta Loudovico Zamboti (Universidade Estadual de Londrina, Brasil). The coordinates of the borders of \mbox{Paran\'{a}} State were provided by \mbox{Jo\~{a}o} Henrique Caviglione (IAPAR).

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.

Examples

summary(parana)
plot(parana)

Set limits for the parameter values

Description

The functions likfit and variofit in the package geoR

Usage

pars.limits(phi = c(lower = 0, upper = +Inf),
            sigmasq = c(lower = 0, upper = +Inf),
            nugget.rel = c(lower = 0, upper = +Inf),
            kappa = c(lower = 0, upper = +Inf),
            kappa2 = c(lower = 0, upper = +Inf), 
            lambda = c(lower = -3, upper = 3),
            psiR = c(lower = 1, upper = +Inf),
            psiA = c(lower = 0, upper = 2 * pi),
            tausq.rel = nugget.rel)

Arguments

Details

Lower and upper limits for parameter values can be individually specified. For example, including the following in the function call in likfit or variofit:
limits = pars.limits(phi=c(0, 10), lambda=c(-2.5, 2.5)),
will change the limits for the parameters \phi and \lambda. Default values are used if the argument limits is not provided.

Value

A list of a 2 elements vector with limits for each parameters

See Also

likfit, variofit

Examples

pars.limits(phi=c(0,10))
pars.limits(phi=c(0,10), sigmasq=c(0, 100))

Exploratory Geostatistical Plots

Description

This function produces a 2 \times 2 display with the following plots: the first indicates the spatial locations assign different colors to data in different quartiles, the next two shows data against the X and Y coordinates and the last is an histogram of the data values or optionally, a 3-D plot with spatial locations and associated data values.

Usage

## S3 method for class 'geodata'
plot(x, coords=x$coords, data = x$data,
             borders, trend="cte", lambda = 1, col.data = 1,
             weights.divide = "units.m", lowess = FALSE, scatter3d = FALSE,
             density = TRUE, rug = TRUE, qt.col, ...)

Arguments

Value

A plot is produced on the graphics device. No values are returned.

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

points.geodata, scatterplot3d, lowess, density, rug.

Examples

require(geoR)
plot(s100)
plot(s100, scatter3d=TRUE)
plot(s100, qt.col=1)

plot(ca20)                       # original data
plot(ca20, trend=~altitude+area) # residuals from an external trend
plot(ca20, trend='1st')          # residuals from a polynomial trend

plot(sic.100, bor=sic.borders)           # original data
plot(sic.100, bor=sic.borders, lambda=0) # logarithm of the data

Plots Variograms for Simulated Data

Description

This function plots variograms for simulated geostatistical data generated by the function grf.

Usage

## S3 method for class 'grf'
plot(x, model.line = TRUE, plot.locations = FALSE, ...)

Arguments

Value

A plot with the empirical variogram(s) is produced on the output device. No values are returned.

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

grf for simulation of Gaussian random fields, plot.variogram for plotting empirical variogram, variog for computation of empirical variograms and plot for the generic plotting function.

Examples

op <- par(no.readonly = TRUE)
par(mfrow=c(2,1))
sim1 <- grf(100, cov.pars=c(10, .25))
# generates simulated data
plot(sim1, plot.locations = TRUE)
#
# plots the locations and the sample true variogram model
#
par(mfrow=c(1,1))
sim2 <- grf(100, cov.pars=c(10, .25), nsim=10)
# generates 10 simulated data
plot(sim1)
# plots sample variograms for all simulations with the true model
par(op)

Plots Prior and/or Posterior Distributions

Description

Produces plots the priors and posteriors distribuitions for the paramters phi and tausq.rel based on results returned by krige.bayes.

Usage

## S3 method for class 'krige.bayes'
plot(x, phi.dist = TRUE, tausq.rel.dist = TRUE, add = FALSE,
                 type=c("bars", "h", "l", "b", "o", "p"), thin, ...)

Arguments

Value

For plot.krige.bayes a plot is produced or added to the current graphics device. No values are returned.

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

See Also

krige.bayes, barplot, matplot.

Examples

## See documentation for krige.bayes

Plots Profile Likelihoods

Description

This function produces plots of the profile likelihoods computed by the function proflik.

Usage

## S3 method for class 'proflik'
plot(x, pages = c("user", "one", "two"), uni.only, bi.only,
             type.bi = c("contour", "persp"), conf.int = c(0.90, 0.95),
             yaxis.lims = c("conf.int", "as.computed"),
             by.col = TRUE, log.scale = FALSE, use.splines = TRUE,
             par.mar.persp = c(0, 0, 0, 0), ask = FALSE, ...)

Arguments

Value

Produces plots with the profile likelihoods on the current graphics device. No values are returned.

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

proflik for computation of the profile likelihoods. For the generic plotting functions see plot, contour, persp. See spline for interpolation.

Examples

# see examples in the documentation for the function proflik()

Plot Directional Variograms

Description

This function plot directional variograms computed by the function variog4. The omnidirectional variogram can be also included in the plot.

Usage

## S3 method for class 'variog4'
plot(x, omnidirectional=FALSE, same.plot=TRUE, legend = TRUE, ...)

Arguments

Value

A plot is produced on the output device. No values returned.

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

variog4 for variogram calculations and matplot for multiple lines plotting.

Examples

s100.v4 <- variog4(s100, max.dist=1)
# Plotting variograms for the four directions
plot(s100.v4)
# changing plot options
plot(s100.v4, lwd=2)
plot(s100.v4, lty=1, col=c("darkorange", "darkblue", "darkgreen","darkviolet"))
plot(s100.v4, lty=1, lwd=2)
# including the omnidirectional variogram
plot(s100.v4, omni=TRUE)
# variograms on different plots
plot(s100.v4, omni=TRUE, same=FALSE)

Plot Empirical Variogram

Description

Plots sample (empirical) variogram computed using the function variog.

Usage

## S3 method for class 'variogram'
plot(x, max.dist, vario.col = "all", scaled = FALSE,
               var.lines = FALSE, envelope.obj = NULL,
               pts.range.cex, bin.cloud = FALSE,  ...)

Arguments

Details

This function plots empirical variograms. Toghether with lines.variogram can be used to compare sample variograms of different variables and to compare variogram models against the empirical variogram.

It uses the function matplot when plotting variograms for more them one variable.

Value

Produces a plot with the sample variogram on the current graphics device. No values are returned.

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk

References

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

See Also

variog for variogram calculations, lines.variogram and lines.variomodel for adding lines to the current plot, variog.model.env and variog.mc.env for variogram envelops computation, matplot for multiple lines plot and plot for generic plot function.

Examples

op <- par(no.readonly = TRUE)
sim <- grf(100, cov.pars=c(1, .2)) # simulates data
vario <- variog(sim, max.dist=1)   # computes sample variogram
par(mfrow=c(2,2))
plot(vario)                     # the sample variogram
plot(vario, scaled = TRUE)      # the scaled sample variogram
plot(vario, max.dist = 1)       # limiting the maximum distance
plot(vario, pts.range = c(1,3)) # points sizes proportional to number of pairs
par(op)

Plot Cross-Validation Results

Description

This function produces ten plots with the results produced by the cross-validation function xvalid.

Usage

## S3 method for class 'xvalid'
plot(x, coords, borders = NULL, ask = TRUE,
            error = TRUE, std.error = TRUE, data.predicted = TRUE,
            pp = TRUE, map = TRUE, histogram = TRUE,
            error.predicted = TRUE, error.data = TRUE, ...)

Arguments

Details

The number of plots to be produced will depend on the input options. If the graphics device is set to just one plot (something equivalent to par(mfcol=c(1,1))) after each graphic being displayed the user will be prompt to press <return> to see the next graphic.

Alternativaly the user can set the graphical parameter to have several plots in one page. With default options for the arguments the maximum number of plots (10) is produced and setting par(mfcol=c(5,2))) will display them in the same page.

The “errors” for the plots are defined as

error = data - predicted

and the plots uses the color blue to indicate positive errors and red to indicate negative erros.

Value

No value returned. Plots are produced on the current graphics device.

See Also

xvalid for the cross-validation computations.

Examples

wls <- variofit(variog(s100, max.dist = 1), ini = c(.5, .5), fix.n = TRUE)
xvl <- xvalid(s100, model = wls)
#
op <- par(no.readonly = TRUE)
par(mfcol = c(3,2))
par(mar = c(3,3,0,1))
par(mgp = c(2,1,0))
plot(xvl, error = FALSE, ask = FALSE)
plot(xvl, std.err = FALSE, ask = FALSE)
par(op)

Plots Spatial Locations and Data Values

Description

This function produces a plot with points indicating the data locations. Arguments can control the points sizes, patterns and colors. These can be set to be proportional to data values, ranks or quantiles. Alternatively, points can be added to the current plot.

Usage

## S3 method for class 'geodata'
points(x, coords=x$coords, data=x$data, data.col = 1, borders,
               pt.divide=c("data.proportional","rank.proportional",
                           "quintiles", "quartiles", "deciles", "equal"),
               lambda = 1, trend = "cte", abs.residuals = FALSE,
               weights.divide = "units.m", cex.min, cex.max, cex.var,
               pch.seq, col.seq, add.to.plot = FALSE,
               x.leg, y.leg = NULL, dig.leg = 2, 
               round.quantiles = FALSE, permute = FALSE, ...)

Arguments

Details

The points can be devided in categories and have different sizes and/or colours according to the argument pt.divide. The options are:

"data.proportional"

sizes proportional to the data values.

"rank.proportional"

sizes proportional to the rank of the data.

"quintiles"

five different sizes according to the quintiles of the data.

"quartiles"

four different sizes according to the quartiles of the data.

"deciles"

ten different sizes according to the deciles of the data.

"equal"

all points with the same size.

a scalar

defines a number of quantiles, the number provided defines the number of different points sizes and colors.

a numerical vector with quantiles and length > 1

the values in the vector will be used by the function cut as break points to divide the data in classes.

For cases where points have different sizes the arguments cex.min and cex.max set the minimum and the maximum point sizes. Additionally, pch.seq can set different patterns for the points and col.seq can be used to define colors. For example, different colors can be used for quartiles, quintiles and deciles while a sequence of gray tones (or a color sequence) can be used for point sizes proportional to the data or their ranks. For more details see the section EXAMPLES.

The argument cex.var allows for displaying 2 variables at once. In this case one variable defines the backgroung colour of the points and the other defines the points size.

The argument permute if set to TRUE randomly realocates the data in the coordinates. This may be used to contrast the spatial pattern of original data against another situation where there is no spatial dependence (when setting permute = TRUE). If a trend is provided the residuals (and not the original data) are permuted.

Value

A plot is created or points are added to the current graphics device.
A list with graphical parameters used to produce the plot is returned invisibily. According to the input options, the list has some or all of the following components:

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

plot.geodata for another display of the data and points and plot for information on the generic R functions. The documentation of par provides details on graphical parameters. For color schemes in R see gray and rainbow.

Examples

op <- par(no.readonly = TRUE)
par(mfrow=c(2,2), mar=c(3,3,1,1), mgp = c(2,1,0))
points(s100, xlab="Coord X", ylab="Coord Y")
points(s100, xlab="Coord X", ylab="Coord Y", pt.divide="rank.prop")
points(s100, xlab="Coord X", ylab="Coord Y", cex.max=1.7,
               col=gray(seq(1, 0.1, l=100)), pt.divide="equal")
points(s100, pt.divide="quintile", xlab="Coord X", ylab="Coord Y")
par(op)

points(ca20, pt.div='quartile', x.leg=4900, y.leg=5850)

par(mfrow=c(1,2), mar=c(3,3,1,1), mgp = c(2,1,0))
points(s100, main="Original data")
points(s100, permute=TRUE, main="Permuting locations")

## Now an example using 2 variable, 1 defining the
## gray scale and the other the points size
points.geodata(coords=camg[,1:2], data=camg[,3], col="gray",
               cex.var=camg[,5])
points.geodata(coords=camg[,1:2], data=camg[,3], col="gray",
               cex.var=camg[,5], pt.div="quint")

Coordinates of Points Inside a Polygon

Description

This function builds a rectangular grid and extracts points which are inside of an internal polygonal region.

Usage

polygrid(xgrid, ygrid, borders, vec.inout = FALSE, ...)

Arguments

Details

The function works as follows: First it creates a grid using the R function expand.grid and then it uses the geoR' internal function .geoR_inout() which wraps usage of SpatialPoints and over from the package sp to extract the points of the grid which are inside the polygon.

Within the package geoR this function is typically used to select points in a non-rectangular region to perform spatial prediction using krige.bayes, krige.conv or ksline. It is also useful to produce image or perspective plots of the prediction results.

Value

A list with components:

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

pred_grid, expand.grid, over, SpatialPoints.

Examples

 poly <- matrix(c(.2, .8, .7, .1, .2, .1, .2, .7, .7, .1), ncol=2)
 plot(0:1, 0:1, type="n")
 lines(poly)
 poly.in <- polygrid(seq(0,1,l=11), seq(0,1,l=11), poly, vec=TRUE)
 points(poly.in$xy)

Pratical range for correlation functions

Description

Computes practical ranges for the correlation functions implemented in the geoR package

Usage

practicalRange(cov.model, phi, kappa = 0.5, correlation = 0.05, ...)

Arguments

Value

A scalar with the value of the practical range.

See Also

cov.spatial

Examples

practicalRange("exp", phi=10)
practicalRange("sph", phi=10)
practicalRange("gaus", phi=10)
practicalRange("matern", phi=10, kappa=0.5)
practicalRange("matern", phi=10, kappa=1.5)
practicalRange("matern", phi=10, kappa=2.5)

Generates a 2D Prediction Grid

Description

This function facilitates the generation of a 2D prediction grid for geostatistical kriging.

Usage

pred_grid(coords, y.coords = NULL, ..., y.by = NULL,
          y.length.out = NULL, y.along.with = NULL)

Arguments

Value

An two column data-frame which is on output of expand.grid.

See Also

See seq and expand.grid which are used internally and locations.inside and polygrid to select points inside a border.

Examples

pred_grid(c(0,1), c(0,1), by=0.25) ## create a grid in a unit square
loc0 <- pred_grid(ca20$borders, by=20)
points(ca20)
points(loc0, pch="+")
points(locations.inside(loc0, ca20$border), pch="+", col=2)

Prediction for the bivariate Gaussian common component geostatistical model

Description

Performs prediction for the bivariate Gaussian common component geostatistical model

Usage

## S3 method for class 'BGCCM'
predict(object, locations, borders,
              variable.to.predict = 1, ...)

Arguments

Value

A list of the class BGCCMpred with components:

Warning

This is a new function and still in draft format and pretty much untested.

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

See Also

likfitBGCCM

Examples

# see http://www.leg.ufpr.br/geoR/tutorials/CCM.R

Prints an summary of of the output from likfitBGCCM.

Description

Prints a short version of an object of the class BGCCM.

Usage

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

Arguments

See Also

format for options to format the output.

Computes Profile Likelihoods

Description

Computes profile likelihoods for model parameters previously estimated using the function likfit.

Usage

proflik(obj.likfit, geodata, coords = geodata$coords,
        data = geodata$data, sill.values, range.values,
        nugget.values, nugget.rel.values, lambda.values, 
        sillrange.values = TRUE, sillnugget.values = TRUE,
        rangenugget.values = TRUE, sillnugget.rel.values = FALSE,
        rangenugget.rel.values = FALSE, silllambda.values = FALSE,
        rangelambda.values = TRUE,  nuggetlambda.values = FALSE,
        nugget.rellambda.values = FALSE,
        uni.only = TRUE, bi.only = FALSE, messages, ...)

Arguments

Details

The functions .proflik.* are auxiliary functions used to compute the profile likelihoods. These functions are internally called by the minimization functions when estimating the model parameters.

Value

An object of the class "proflik" which is a list. Each element contains values of a parameter (or a pair of parameters for 2-D profiles) and the corresponding value of the profile likelihood. The components of the output will vary according to the input options.

Note

1. Profile likelihoods for Gaussian Random Fields are usually uni-modal. Unusual or jagged shapes can be due to the lack of the convergence in the numerical minimization for particular values of the parameter(s). If this is the case it might be necessary to pass control arguments to the minimization functions using the argument .... It's also advisable to try the different options for the minimisation.function argument. See documentation of the functions optim and/or nlm for further details.

2. 2-D profiles can be computed by setting the argument uni.only = FALSE. However, before computing 2-D profiles be sure they are really necessary. Their computation can be time demanding since it is performed on a grid determined by the cross-product of the values defining the 1-D profiles.

3. There is no "default strategy" to find reasonable values for the x-axis. They must be found in a "try-and-error" exercise. It's recommended to use short sequences in the initial attempts. The EXAMPLE section below illustrates this.

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

plot.proflik for graphical output, likfit for the parameter estimation, optim and nlm for further details about the minimization functions.

Examples

op <- par(no.readonly=TRUE)
ml <- likfit(s100, ini=c(.5, .5), fix.nug=TRUE)
## a first atempt to find reasonable values for the x-axis:
prof <- proflik(ml, s100, sill.values=seq(0.5, 1.5, l=4),
                range.val=seq(0.1, .5, l=4))
par(mfrow=c(1,2))
plot(prof)
## a nicer setting 
## Not run: 
prof <- proflik(ml, s100, sill.values=seq(0.45, 2, l=11),
                range.val=seq(0.1, .55, l=11))
plot(prof)
## to include 2-D profiles use:
## (commented because this is time demanding)
#prof <- proflik(ml, s100, sill.values=seq(0.45, 2, l=11),
#                range.val=seq(0.1, .55, l=11), uni.only=FALSE)
#par(mfrow=c(2,2))
#plot(prof, nlevels=16)

## End(Not run)
par(op)

Reads and Converts Data to geoR Format

Description

Reads data from a ASCII file and converts it to an object of the class geodata, the standard data format for the geoR package.

Usage

read.geodata(file, header = FALSE, coords.col = 1:2, data.col = 3,
             data.names = NULL, covar.col = NULL, covar.names = "header",
             units.m.col = NULL, realisations = NULL,
             na.action = c("ifany", "ifdata", "ifcovar", "none"),
             rep.data.action, rep.covar.action, rep.units.action, ...)

Arguments

Details

The function read.table is used to read the data from the ASCII file and then as.geodata is used to convert to an object of the class geodata.

Value

An object of the class geodata. See documentation for the function as.geodata for further details.

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

as.geodata to convert existing R objects, read.table, the basic R function used to read ASCII files, and list for detailed information about lists.

Radionuclide Concentrations on Rongelap Island

Description

This data-set was used by Diggle, Tawn and Moyeed (1998) to illustrate the model-based geostatistical methodology introduced in the paper. discussed in the paper. The radionuclide concentration data set consists of measurements of \gamma-ray counts at 157 locations.

Usage

data(rongelap)

Format

The object is a list with the following components:

coords

the coordinates of data locations.

data

the data.

units.m

n-dimensional vector of observation-times for the data.

borders

a matrix with the coordinates defining the coastline on Rongelap Island.

Source

For further details on the radionuclide concentration data, see Diggle, Harper and Simon (1997), Diggle, Tawn and Moyeed (1998) and Christensen (2004).

References

Christensen, O. F. (2004). Monte Carlo maximum likelihood in model-based geostatistics. Journal of computational and graphical statistics 13 702-718.

Diggle, P. J., Harper, L. and Simon, S. L. (1997). Geostatistical analysis of residual contamination from nuclea testing. In: Statistics for the environment 3: pollution assesment and control (eds. V. Barnet and K. F. Turkmann), Wiley, Chichester, 89-107.

Diggle, P. J., Tawn, J. A. and Moyeed, R. A. (1998). Model-based geostatistics (with Discussion). Applied Statistics, 47, 299–350.

Simulated Data-Sets which Illustrate the Usage of the Package geoR

Description

These two simulated data sets are the ones used in the Technical Report which describes the package geoR (see reference below). These data-sets are used in several examples throughout the package documentation.

Usage

data(s100)

data(s121)

Format

Two objects of the class geodata. Both are lists with the following components:

coords

the coordinates of data locations.

data

the simulated data. Notice that for s121 this a 121 \times 10 matrix with 10 simulations.

cov.model

the correlation model.

nugget

the values of the nugget parameter.

cov.pars

the covariance parameters.

kappa

the value of the parameter kappa.

lambda

the value of the parameter lambda.

References

Ribeiro Jr, P.J. and Diggle, P.J. (1999) geoS: A geostatistical library for S-PLUS. Technical report ST-99-09, Dept of Maths and Stats, Lancaster University.

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

Examples

plot(s100)
plot(s121, type="l")

Simulated Data-Set which Illustrate the Usage of krige.bayes

Description

This is the simulated data-set used in the Technical Report which describes the implementation of the function krige.bayes (see reference below).

Usage

data(s256i)

Format

Two objects of the class geodata. Both are lists with the following components:

coords

the coordinates of data locations.

data

the simulated data.

References

Ribeiro Jr, P.J. and Diggle, P.J. (1999) Bayesian inference in Gaussian model-based geostatistics. Technical report ST-99-08, Dept of Maths and Stats, Lancaster University.

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

Examples

points(s256i, pt.div="quintiles", cex.min=1, cex.max=1)

Sampling from geodata objects

Description

This functions facilitates extracting samples from geodata objects.

Usage

sample.geodata(x, size, replace = FALSE, prob = NULL, coef.logCox,
               external)

Arguments

Details

If prob=NULL and the argument coef.logCox, is provided, sampling follows a log-Cox proccess, i.e. the probability of each point being sampled is proportional to:

exp(b Y(x))

with b given by the value passed to the argument coef.logCox and Y(x) taking values passed to the argument external or, if this is missing, the element data of the geodata object. Therefore, the latter generates a preferential sampling.

Value

a list which is an object of the class geodata.

See Also

as.geodata, sample.

Examples

## Not run: 
par(mfrow=c(1,2))
S1 <- grf(2500,  grid="reg", cov.pars=c(1, .23))
image(S1, col=gray(seq(0.9,0.1,l=100)))
y1 <- sample.geodata(S1, 80)
points(y1$coords, pch=19)
## Now a preferential sampling
y2 <- sample.geodata(S1, 80, coef=1.3)
## which is equivalent topps
## y2 <- sample.geodata(S1, 80, prob=exp(1.3*S1$data))
points(y2$coords, pch=19, col=2)
## and now a clustered (but not preferential)
S2 <- grf(2500,  grid="reg", cov.pars=c(1, .23))
y3 <- sample.geodata(S1, 80, prob=exp(1.3*S2$data))
## which is equivalent to
## points(y3$coords, pch=19, col=4)
image(S2, col=gray(seq(0.9,0.1,l=100)))
points(y3$coords, pch=19, col=4)

## End(Not run)

Samples from the posterior distribution

Description

Sample quadruples (\beta, \sigma^2, \phi, \tau^2_{rel}) from the posterior distribution returned by krige.bayes.

Usage

sample.posterior(n, kb.obj)

Arguments

Value

A n \times 4 data-frame with samples from the posterior distribution of the model parameters.

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

krige.bayes and sample.posterior.

Sample the prior distribution

Description

Sample quadruples (\beta, \sigma^2, \phi, \tau^2_{rel}) from the prior distribution of parameters specifying a Gaussian random field. Typically the prior is specified in the same manner as when calling krige.bayes.

Usage

sample.prior(n, kb.obj=NULL, prior=prior.control())

Arguments

Value

A p+3 \times 4 data-frame with a sample of the prior distribution of model parameters, where p is the length of the mean parameter \beta.

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

See Also

krige.bayes and sample.posterior.

Examples

sample.prior(50, 
             prior=prior.control(beta.prior = "normal", beta = .5, beta.var.std=0.1, 
                                 sigmasq.prior="sc", sigmasq=1.2, df.sigmasq= 2, 
                                 phi.prior="rec", phi.discrete = seq(0,2, l=21)))

Sets Limits to Scale Plots

Description

This is an function typically called by functions in the package geoR to set limits for the axis when plotting spatial data.

Usage

set.coords.lims(coords, borders = coords, xlim, ylim, ...)

Arguments

Value

A 2 \times 2 matrix with limits for the axis.

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

Soil chemistry properties data set

Description

Several soil chemistry properties measured on a regular grid with 10x25 points spaced by 5 meters.

Usage

data(soil250)

Format

A data frame with 250 observations on the following 22 variables.

Linha

x-coordinate

Coluna

y-coordinate

Cota

elevation

AGrossa

a numeric vecto, sand portion of the sample.

Silte

a numeric vector, silt portion of the sample.

Argila

a numeric vector, sand portion of the sample.

pHAgua

a numeric vector, soil pH at water

pHKCl

a numeric vector, soil pH by KCl

Ca

a numeric vector, calcium content

Mg

a numeric vector, magnesium content

K

a numeric vector, potassio content

Al

a numeric vector, aluminium content

H

a numeric vector, hidrogen content

C

a numeric vector, carbon content

N

a numeric vector, nitrogen content

CTC

a numeric vector, catium exchange capability

S

a numeric vector, enxofrar content

V

a numeric vector

M

a numeric vector

NC

a numeric vector

CEC

a numeric vector

CN

a numeric vector, carbon/nitrogen relation

Details

Uniformity trial with 250 undisturbed soil samples collected at 25cm soil depth of spacing of 5 meters, resulting on a regular grid of 25 \times 10 points.

See also the data-set wrc with other variables colected at the same points.

Source

Bassoi thesis

References

Bassoi papers

Examples

data(soil250)
ctc <- as.geodata(soil250, data.col=16)
plot(ctc)

Soya bean production and other variables in a uniformity trial

Description

Data on soya bean production in a uniformity trial measured at plots of size 5x5 meters and other soil properties measured in points given by the data coordinates.

Usage

data(soja98)

Format

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

X

a numeric vector with X-coordinates of the plot centres.

Y

a numeric vector with X-coordinates of the plot centres.

P

a numeric vector, phosphorous content.

PH

a numeric vector, soil pH.

K

a numeric vector, potassium content. (Cmol/dm^3)

MO

a numeric vector, organic matter. (percentage)

SB

a numeric vector, basis saturation.

iCone

a numeric vector, cone index, measuring mechanic resistence of the soil. (kg/cm^2)

Rend

a numeric vector, total soya production within the plot (kg).

PROD

a numeric vector, production converted to productivity by a unit of area - hectare (ton/ha).

Source

Souza, E.G.; Jojann, J. A.; Rocha, J. V.; Ribeiro, S. R. A.; Silva, M. S., Upazo, M. A. U.; Molin, J. P.; Oliveira, E. F.; Nóbrega, L. H. P. (1999) Variabilidade espacial dos atributos químicos do solo em um latossolo roxo distrófico da região de Cascavel-PR. Engenharia Agrícola. Jaboticabal, volume 18, nr. 3, p.80-92.

Examples

data(soja98)
plot(soja98)
require(geoR)
prod98 <- as.geodata(soja98, coords.col=1:2, data.col=)
plot(prod98, low=TRUE)

Summary statistics from predictive distributions

Description

Computes summaries based on simulations of predictive distribution returned by krige.bayes and krige.conv.

Usage

statistics.predictive(simuls, mean.var = TRUE, quantile, threshold,
                      sim.means, sim.vars)

Arguments

Value

A list with one ore more of the following components.

Author(s)

Paulo J. Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

References

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

Selects a Subarea from a Geodata Object

Description

Selects elements of a geodata object wich are within a rectangular (sub)area

Usage

subarea(geodata, xlim, ylim, ...)

Arguments

Details

The function copies the original geodata object and selects values of $coords, $data, $borders, $covariate and $units.m which lies within the selected sub-area. Remaining components of the geodata objects are untouched.

If xlim and/or ylim are not provided the function prompts the user to click 2 points defining an retangle defining the subarea on a existing plot.

Value

Returns an geodata object, subsetting of the original one provided.

Author(s)

Paulo Justiniano Ribeiro Jr. paulojus@leg.ufpr.br,
Peter J. Diggle p.diggle@lancaster.ac.uk.

See Also

zoom.coords, locator

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)

#if(require(geoRglm)){
#data(rongelap)
#points(rongelap)
## zooming the western area
#rongwest <- subarea(rongelap, xlim=c(-6300, -4800))
#points(rongwest)
## now zooming in the same plot
#points(rongelap)
#rongwest.z <- zoom.coords(rongwest, xzoom=3, xoff=2000, yoff=3000)
#points(rongwest.z, add=TRUE)
#rect.coords(rongwest$sub, quiet=TRUE)
#rect.coords(rongwest.z$sub, quiet=TRUE)
#}

Method for subsetting geodata objects

Description

Subsets a object of the class geodata by transforming it to a data-frame, using subset and back transforming to a geodata object.

Usage

## S3 method for class 'geodata'
subset(x, ..., other = TRUE)

Arguments

Value

A list which is an object of the class geodata.

See Also

subset for the generic function and methods and as.geodata for more information on geodata objects.

Examples

subset(ca20, data > 70)
subset(ca20, area == 1)

Summaries for geodata object

Description

Sumarises each of the main elements of an object of the class geodata.

Usage

## S3 method for class 'geodata'
summary(object, lambda =1, add.to.data = 0,
                by.realisations=TRUE, ...)

Arguments

Value

A list with components