krige.bayes {geoR}R Documentation

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 = NULL, model = model.control(),
            prior = prior.control(), output = output.control())

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"),
              tausq.rel = 0, tausq.rel.discrete = NULL)

Arguments

geodata a list containing elements coords and data as described next. Typically an object of the class "geodata" - a geoR data-set. If not provided the arguments coords and data must be provided instead.
coords an n x 2 matrix where each row has the 2-D coordinates of the n data locations. By default it takes the component coords of the argument geodata, if provided.
data a vector with n data values. By default it takes the component data of the argument geodata, if provided.
locations an N x 2 matrix or data-frame with the 2-D coordinates of the N prediction locations. Defaults to "no" in which case the function returns only results on the posterior distributions of the model parameters.
borders optional. If a two column matrix defining a polygon is provided the prediction is performed only at locations inside this polygon.
model defines model components. See section DETAILS below.
prior specification of priors for the model parameters. See section DETAILS below.
output Specify output options. See documentation for output.control.
trend.d specifies the trend (covariates) values at the data locations. Possible values are: "cte" - model with constant mean, "1st" - trend is defined as a first order polynomial on the coordinates, "2nd" - trend is defined as a second order polynomial on the coordinates, a formula of the type ~X, where X is a matrix with covariates (external trend) at data locations. Defaults to "cte".
trend.l specifies the trend (covariates) at the prediction locations. Must be of the same type as defined for trend.d. Only used if prediction locations are provided in the argument locations.
cov.model string indicating the name of the model for the correlation function. Further details in the documentation for cov.spatial.
kappa additional smoothness parameter. Only used if the correlation function is one of: "matern", "powered.exponential", "cauchy" or "gneiting.matern". In the current implementation this parameter is always regarded as fixed during the Bayesian analysis.
aniso.pars fixed parameters for geometric anisotropy correction. If aniso.pars = FALSE no correction is made, otherwise a two elements vector with values for the anisotropy parameters must be provided. Anisotropy correction consists of a transformation of the data and prediction coordinates performed by the function coords.aniso.
lambda numerical value of the Box-Cox transformation parameter. The value lambda = 1 corresponds to no transformation. The Box-Cox parameter lambda is always regarded as fixed and data transformation is performed before the analysis. Prediction results are back-transformed and returned is the same scale as for the original data. For lambda = 0 the log-transformation is performed. If lambda < 0 the mean predictor doesn't make sense (the resulting distribution has no expectation).
beta.prior prior distribution for the mean (vector) parameter beta. The options are "flat" (default), "normal" or "fixed" (known mean).
beta mean hyperparameter for the distribution of the mean (vector) parameter beta. Only used if beta.prior = "normal" or beta.prior = "fixed". For the later beta defines the value of the known mean.
beta.var.std standardised (co)variance hyperparameter(s) for the prior for the mean (vector) parameter beta. The (co)variance matrix forbeta is given by the multiplication of this matrix by sigma^2. Only used if beta.prior = "normal".
sigmasq.prior specifies the prior for the parameter sigma^2. If "reciprocal" (the default), the prior 1/sigma^2 is used. Otherwise the parameter is regarded as fixed.
sigmasq fixed value of the sill parameter sigma^2. Only used if sigmasq.prior = FALSE.
df.sigmasq numerical. Number of degrees of freedom for the prior for the parameter sigma^2. Only used if sigmasq.prior = "sc.inv.chisq".
phi.prior prior distribution for the range parameter phi. Options are: \"uniform\", (p(phi) propto 1), \"exponential\" ,((-1/nu) exp(- (1/nu) * phi)) \"reciprocal\" (p(phi) propto 1/phi), \"squared.reciprocal\" (p(phi) propto 1/phi^2), \"fixed\" (known/estimated value of phi). If the latter is defined the parameter phi is fixed when performing prediction. For the exponential prior the argument phi provides the the hyperparameter nu and corresponds to the expected value for this distribution.
phi fixed value of the range parameter phi. Only needed if phi.prior = "fixed".
phi.discrete support points of the discrete prior for the range parameter phi.
tausq.rel.prior specifies a prior distribution for the relative nugget parameter tau^2/sigma^2. If tausq.rel.prior = "fixed" the relative nugget is considered known (fixed) with value given by the argument tausq.rel. If tausq.rel.prior = "uniform" a discrete uniform prior is used with support points given by the argument tausq.rel.discrete.
tausq.rel fixed value for the relative nugget parameter. Only used if tausq.rel.prior = "fixed".
tausq.rel.discrete support points of the discrete prior for the relative nugget parameter tau^2/sigma^2.

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.

The basis for the Bayesian algorithm is to discretize the prior distribution for the parameters phi and tau_rel = tau/sigma. The Tech. Report referenced below provides details on the results used in the current implementation.

CONTROL FUNCTIONS

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

Value

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

posterior results on on the posterior distribution of the model parameters. A list with the following possible components:
beta summary information on the posterior distribution of the mean parameter beta.
sigmasq summary information on the posterior distribution of the variance parameter sigma^2 (partial sill).
phi summary information on the posterior distribution of the correlation parameter phi (range parameter) .
tausq.rel summary information on the posterior distribution of the relative nugget variance parameter tau^2_rel.
joint.phi.tausq.rel information on discrete the joint distribution of these parameters.
sample a data.frame with a sample from the posterior distribution. Each column correspondes to one of the basic model parameters.
predictive results on the predictive distribution at the prediction locations, if provided. A list with the following possible components:
mean expected values.
variance expected variance.
distribution type of posterior distribution.
mean.simulations mean of the simulations at each locations.
variance.simulations variance of the simulations at each locations.
quantiles.simulations quantiles computed from the the simulations.
probabilities.simulations probabilities computed from the simulations.
simulations simulations from the predictive distribution.
prior a list with information on the prior distribution and hyper-parameters of the model parameters (beta, sigma^2, phi, tau^2_rel).
model model specification as defined by model.control.
.Random.seed system random seed before running the function. Allows reproduction of results. If the .Random.seed is set to this value and the function is run again, it will produce exactly the same results.
max.dist maximum distance found between two data locations.
call the function call.

Author(s)

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

References

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.maths.lancs.ac.uk/~ribeiro/publications.html

Further information about geoR can be found at:
http://www.maths.lancs.ac.uk/~ribeiro/geoR.html.

See Also

lines.krige.bayes, image.krige.bayes and persp.krige.bayes for graphical output of the results. krige.conv and ksline for conventional kriging methods.

Examples


# generating a simulated data-set
ex.data <- grf(50, cov.pars=c(10, .25))
#
# defining the prediction grid:
ex.grid <- as.matrix(expand.grid(seq(0,1,l=21), seq(0,1,l=21)))
#
# computing Bayesian posterior and predictive distributions
# (this can take some time to run)
ex.bayes <- krige.bayes(ex.data, loc=ex.grid, prior =
                 prior.control(phi.discrete=seq(0, 2, l=21)))
#
# Prior and posterior for the parameter phi
plot(ex.bayes, type="h", tausq.rel = FALSE, col=c("red", "blue"))
#
# Plot histograms of samples from the posterior
par(mfrow=c(2,1))
hist(ex.bayes)
par(mfrow=c(1,1))
#
# Ploting theoretical and empirical variograms
plot(ex.data)
# adding lines with fitted variograms
lines(ex.bayes)
lines(ex.bayes, summ="median", lty=2)
lines(ex.bayes, summ="mean", lwd=2, lty=2)
#
# Ploting prediction some results
op <- par(no.readonly = TRUE)
par(mfrow=c(2,2))
par(mar=c(3,3,1,1))
par(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)
[Package Contents]