# krige.bayes {geoR}

Bayesian Analysis for Gaussian Geostatistical Models
Package:
geoR
Version:
1.7-5.2

### Description

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

### Usage

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

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

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

post2prior(obj)
```

### Arguments

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

### Details

`krige.bayes` is a generic function for Bayesian geostatistical analysis of (transformed) Gaussian where predictions take into account the parameter uncertainty. It can be set to run conventional kriging methods which use known parameters or plug-in estimates. However, the functions `krige.conv` and `ksline` are preferable for prediction with fixed parameters.

PRIOR SPECIFICATION The basis of the Bayesian algorithm is the discretisation of the prior distribution for the parameters phi and tau_rel = tau/sigma. 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":
(-1/nu) exp(- (1/nu) * phi).
"reciprocal":
p(φ) propto 1/phi.
"squared.reciprocal":
p(φ) propto 1/phi^2.
"fixed":
fixed known or estimated value of phi.

The expressions of the implemented priors for the parameter tausq.rel are:

"fixed":
fixed known or estimated value of tausq.rel. Defaults to zero.
"uniform":
p(tausq.rel) propto 1.
"reciprocal":
p(tausq.rel) propto 1/tausq.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.

### Values

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.relinformation 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.
posterior
results on on the posterior distribution of the model parameters. A list with the following possible components:
predictive
results on the predictive distribution at the prediction locations, if provided. A list with the following possible components:
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.

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

`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
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
predictive distribution")
image(ex.bayes, val= "simulation", number.col=2,
main="another simulation from
the 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
##```

### Author(s)

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

Documentation reproduced from package geoR, version 1.7-5.2. License: GPL (>= 2)