R: Bayesian Analysis for Gaussian Geostatistical Models

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",
     model = model.control(trend.d = "cte", trend.l = "cte",
                   cov.model = "exponential", kappa = 0.5,
                   aniso.pars = NULL, lambda = 1),
     prior = prior.control(beta.prior = c("flat", "normal",
                                          "fixed"),
                   beta = NULL, beta.var = NULL,
                   sill.prior = c("reciprocal", "fixed"),
                   sill = NULL,
                   range.prior = c("uniform", "exponential",
                                   "fixed", "squared.reciprocal",
                                   "reciprocal"),
                   exponential.prior.par = 1, range = NULL,
                   range.rel.discrete = NULL,
                   nugget.rel.prior = c("fixed", "uniform"),
                   nugget.fixed = 0, nugget.discrete = NULL),
     output = output.control(n.posterior = 1000,
                   n.predictive = NULL, moments = TRUE,
                   simulations.predictive = TRUE,
                   keep.simulations = TRUE, mean.estimator = TRUE,
                   quantile.estimator =  NULL,
                   probability.estimator = NULL,
                   signal = FALSE, messages.screen = TRUE))

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.
`model`	defines model components. See section DETAILS below.
`prior`	specification of priors for the model parameters. See section DETAILS below.
`output`	Defines output options. See section DETAILS below.

Details

krige.bayes is a generic function for Bayesian geostatistical analysis where predictions can 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. The arguments for the control functions are as follows.

ARGUMENTS FOR CONTROL FUNCTIONS

Arguments for model = model.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).

Arguments for prior = prior.control(...) :

beta.prior: prior distribution for the mean (vector) parameter beta. The options are "flat" (default), "normal" or "fixed" (known mean).
beta: 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: (co)variance hyperparameters for the prior of the mean (vector) parameter beta. Only used if beta.prior = "normal".
sill.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.
sill: fixed value of the sill parameter sigma^2. Only used if sill.prior = FALSE.
range.prior: prior distribution for the range parameter phi. Options are: \"uniform\", (propto 1), \"exponential\" ,(exp(- nu * phi)) \"reciprocal\" (1/phi), \"squared.reciprocal\" (1/phi^2), \"fixed\" (known/estimated value of phi). If the latter is defined the parameter phi is fixed when performing prediction.
exponential.prior.par: parameter nu for the exponential prior for the parameter phi. Only used if prior.phi = "exponential".
range: fixed value of the range parameter phi. Only needed if range.prior = "fixed".
range.discrete: support points for the discretisation of the prior for the range parameter phi.
nugget.rel.prior: specifies a prior distribution for the relative nugget parameter tau^2/sigma^2. If nugget.rel.prior = "fixed" the nugget is considered known (fixed) with value given by the argument nugget.fixed. If nugget.rel.prior = "uniform" a discrete uniform prior is used with support points given by the argument nugget.rel.discrete.
nugget.fixed: fixed value for the nugget parameter. Only used if prior.nugget = "fixed".
nugget.rel.discrete: support points for the discretisation of the prior for the relative nugget parameter tau^2/sigma^2.

Arguments for output = output.control(...) :

n.posterior: number of samples to be taken from the posterior distribution.
n.predictive: number of samples to be taken from the predictive distribution. By default equals to n.posterior.
moments: logical. If TRUE moments of the predictive distribution are computed analytically (without sampling). Valid only if lambda = 1 or lambda = 0.
simulations.predictive: logical. Defines whether simulations are drawn from the predictive distribution. Only valid if prediction locations are provided on the argument locations.
keep.simulations: logical. Indicates whether or not the samples of the predictive distributions are returned. Only valid if prediction locations are provided on the argument locations.
mean.estimator: logical. Indicates whether or not the mean and variances of the predictive distributions are computed and returned. If TRUE the objects predict.mean and krige.var are included in the output. Only valid if prediction locations are provided on the argument locations.
quantile.estimator: indicates whether or not quantiles of the predictive distributions are computed and returned. If a vector with numbers in the interval [0,1] is provided the output includes the object quantiles, which contains values of corresponding estimated quantiles. For example, if estimator = c(0.25,0.50,0.75) the function returns the quartiles of the distributions at each of the prediction locations. If quantile.estimator = TRUE, the default, the vector c(0.025, 0.5, 0.975), is assumed. A measure of uncertainty for the predictions, which is analogous to the kriging standard error, can be computed by (quantile0.975 - quantile0.025)/4. Only used if prediction locations are provided in the argument locations.
probability.estimator: the default is FALSE for which case nothing is computed. If some cutoff values are provided instead, an object called probability is included in the output. This object contains, for each prediction location, the probability that the variable is less than or equal to the cutoff value given in the argument.
signal: logical. If TRUE the signal is predicted, otherwise the variable is predicted. If no transformation is performed the expectations are the same in both cases and the kriging variances are different, if the nugget is different of zero.
messages.screen: logical. Indicates whether or not status messages are printed on the screen (or output other device) while the function is running.

Value

An object of the class "krige.bayes" which is a list with the following components:

`posterior`	A list with results for the posterior distribution of the model parameters. The components are:
`beta.summary`	summary for the posterior distribution of the mean parameter beta.
`sigmasq.summary`	summary for the posterior distribution of the variance parameter sigma^2 (partial sill).
`phi.summary`	summary for the posterior distribution of the correlation parameter phi (range parameter) .
`tausq.summary`	summary for the posterior distribution of the nugget variance parameter tau^2.
`beta.samples`	samples from the posterior distribution of the mean parameter beta.
`sigmasq.samples`	samples from the posterior distribution of the variance parameter sigma^2.
`phi.samples`	samples from the posterior distribution of the correlation parameter phi.
`tausq.samples`	samples from the posterior distribution of the nugget variance parameter tau^2.
`phi.marginal`	samples from the marginal posterior distribution of the correlation parameter phi, resulting from averaging the posterior over the distribution of (beta, sigma^2).
`nugget.marginal`	samples from the marginal posterior distribution of the nugget variance parameter tau^2, resulting from averaging the posterior over the distribution of (beta, sigma^2).
`predictive`	A list with results for the predictive distribution of the prediction locations (if provided). The components are:
`moments`	a numerical matrix. The columns contains the estimated moments of the predictive distribution, at each prediction location
`simulations`	a numerical matrix. Each column has a simulation from the predictive distribution. Returned only if `keep.simulations = TRUE`.
`mean.simulations`	a vector with the estimated mean at the prediction locations computed by averaging over the simulations. Returned only if `mean.estimator = TRUE`.
`variance.simulations`	a vector with the estimated variance at the prediction locations, computed using the simulations. Returned only if `mean.estimator = TRUE`.
`quantile`	A matrix or vector with quantile estimators. Returned only if the argument `quantile.estimator` is used.
`probability`	A matrix or vector with probability estimators. Returned only if the argument `probability.estimator` is used.
`type.prediction`	information on the type of prediction performed.
`message.prediction`	information about the parameter uncertainty taken into account. Indicates which parameters has been regarded as random during the analysis.
`.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.
`call`	the function call.

Auxiliary functions

The functions of type krige.bayes.aux and control are auxiliary functions called by krige.bayes.

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.

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=11), seq(0,1,l=11)))
#
# 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(range.discrete=seq(0, 2, l=51)))
#
# Ploting theoretical amd 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.krige.bayes(ex.bayes, loc=ex.grid, main="predicted values")
image.krige.bayes(ex.bayes, val="moments.variance",
                  loc=ex.grid, main="prediction variance")
image.krige.bayes(ex.bayes, val= "simulation", number.col=1,
                  loc=ex.grid,
      main="a simulation from the \npredictive distribution")
image.krige.bayes(ex.bayes, val= "simulation", number.col=2,
                  loc=ex.grid,
     main="another simulation from \nthe predictive distribution")
par(op)