geoR (version 1.2-5)

krige.bayes: 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, 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"),
              tausq.rel = 0, 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 m

coords

an $n \times 2$ matrix where each row has the 2-D coordinates of the $n$ data locations. By default it takes the component coords of the argument geodata, if provided.

data

a vector with n data values. By default it takes the component data of the argument geodata, if provided.

locations

an $N \times 2$ matrix or data-frame with the 2-D coordinates of the $N$ prediction locations. 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

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.

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 argu

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 docume

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 parame

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

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.

beta.prior

prior distribution for the mean (vector) parameter $\beta$. The options are "flat" (default), "normal" or "fixed" (known mean).

beta

mean hyperparameter for the distribution of the mean (vector) parameter $\beta$. Only used if beta.prior = "normal" or beta.prior = "fixed". For the later beta defines the value of the known mean.

beta.var.std

standardised (co)variance hyperparameter(s) for the prior for the mean (vector) parameter $\beta$. The (co)variance matrix for$\beta$ is given by the multiplication of this matrix by $\sigma^2$. Only used if beta.prior = "normal"

sigmasq.prior

specifies the prior for the parameter $\sigma^2$. If "reciprocal" (the default), the prior $\frac{1}{\sigma^2}$ is used. Otherwise the parameter is regarded as fixed.

sigmasq

fixed value of the sill parameter $\sigma^2$. Only used if sigmasq.prior = FALSE.

df.sigmasq

numerical. Number of degrees of freedom for the prior for the parameter $\sigma^2$. Only used if sigmasq.prior = "sc.inv.chisq".

phi.prior

prior distribution for the range parameter $\phi$. Options are: "uniform", "exponential", "reciprocal" , "squared.reciprocal" and "fixed". Alternativelly, a user defined

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 $\frac{\tau^2}{\sigma^2}$. If tausq.rel.prior = "fixed" the relative nugget is considered known (fixed) with value given by the argument tausq.rel

tausq.rel

fixed value for the relative nugget parameter. Only used if tausq.rel.prior = "fixed".

tausq.rel.discrete

support points of the discrete prior for the relative nugget parameter $\frac{\tau^2}{\sigma^2}$.

obj

an object of the class krige.bayes or posterior.krige.bayes with the output of a call to krige.bayes. The function post2prior takes the posterior distribution computed by one call to kr

Value

• An object with class "krige.bayes" and "kriging". The attribute prediction.locations containing the name of the object with the coordinates of the prediction locations (argument locations) is assigned to the object. Returns a list with the following components:
• posterior
  results on on the posterior distribution of the model parameters. A list with the following possible components:
• • beta
  {summary information on the posterior distribution of the mean parameter $\beta$. }
• sigmasq
  summary information on the posterior distribution of the variance parameter $\sigma^2$ (partial sill).
• phi
  summary information on the posterior distribution of the correlation parameter $\phi$ (range parameter) .
• tausq.rel
  summary information on the posterior distribution of the relative nugget variance parameter $\tau^2_{rel}$.
• joint.phi.tausq.rel
  information on discrete the joint distribution of these parameters.
• sample
  a data.frame with a sample from the posterior distribution. Each column corresponds to one of the basic model parameters.

cr

• 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. }

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 tht 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: [object Object],[object Object],[object Object],[object Object],[object Object]

We remember that 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. 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.

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.

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 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: 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 with samples from the posterior
par(mfrow=c(2,1))
hist(ex.bayes)
par(mfrow=c(1,1))
#
# Plotting empirical variograms and some Bayesian estimates
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)
#
# Plotting some prediction 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)

<testonly>ex.data <- grf(50, cov.pars=c(10, .25))
ex.post <- krige.bayes(ex.data)
ex1 <- krige.bayes(ex.data, prior = list(phi.prior = "fixed", phi = 0.3))
ex1 <- krige.bayes(ex.data, model = list(cov.model="spherical"))
ex1 <- krige.bayes(ex.data, output = list(n.posterior = 100))
ex.grid <- as.matrix(expand.grid(seq(0,1,l=6), seq(0,1,l=6)))
ex.bayes <- krige.bayes(ex.data, loc=ex.grid, prior =
                 prior.control(phi.discrete=seq(0, 2, l=3),
                 tausq.rel.discrete=seq(0, 2, l=3)),
                 output=output.control(n.post=100))
plot(ex.data)
lines(ex.bayes)
plot(ex.bayes)
lines(ex.bayes, summ="median", lty=2)
lines(ex.bayes, summ="mean", lwd=2, lty=2)
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 \npredictivedistribution")
image(ex.bayes, val= "simulation", number.col=2,
      main="another simulation from \nthepredictive distribution")
PC <- prior.control(phi.prior = c(.1, .2, .3, .2, ,.1, .1), phi.disc=seq(0.1, 0.6, l=6), tausq.rel.prior=c(.1, .4, .3, .2), tausq.rel.discrete=c(0,.1,.2,.3))
ex.user <- krige.bayes(ex.data, prior = PC)  

# Simulating data at 2 different "times"
ap1 <- grf(50, cov.pars=c(1, .3))
ap2 <- grf(70, cov.pars=c(1, .3))
## A initial "usual" analysis
ap1.kb <- krige.bayes(ap1)
## using the previous posterior as prior for next call
ap2.kb <- krige.bayes(ap2, prior=post2prior(ap1.kb))
##
## Another example with "user defined" prior
##
PC <- prior.control(phi.prior=c(.2,.3,.2,.1,.1,.1), phi.discrete = seq(0,.5,l=6))
ap3.kb <- krige.bayes(ap1, prior = PC)
##
ap4.kb <- krige.bayes(ap2, prior=post2prior(ap3.kb))
##
## Now include tausq
##
PC <- prior.control(tausq.rel.prior = "uni", tausq.rel.discrete = seq(0, .5, l=6))
ap5.kb <- krige.bayes(ap1)
## using the previous posterior as prior for next call
ap6.kb <- krige.bayes(ap2, prior=post2prior(ap5.kb))
##
##
##
PC <- prior.control(phi.prior=c(.2,.3,.2,.1,.1,.1), phi.discrete = seq(0,.5,l=6), tausq.rel.prior=c(.3,.4,.3), tausq.rel.discrete = c(0, .1, .2)) 
ap7.kb <- krige.bayes(ap1, prior = PC)
#
ap8.kb <- krige.bayes(ap2, prior=post2prior(ap7.kb))

## with trend
data(s100)
prior2.b9 <- prior.control(beta.prior = "normal", beta = c(0,0,0),
beta.var = cbind(c(2,1.5,0),c(1.5,1.8,.2),c(0,0.2,1.5)),
phi.prior = "exponential", phi = 2.5, phi.discrete = c(2.5,3),
sigmasq.prior = "sc.inv.chisq", df.sigmasq = 5, sigmasq = 0.5)
ap <- krige.bayes(s100,prior=prior2.b9, model=model.control(trend.d = "1st") )

prior2.b9 <- prior.control(beta.prior = "normal", beta = c(0,0,0),
beta.var = cbind(c(2,1.5,0),c(1.5,1.8,0.5),c(0,0.5,1.5)),
phi.prior = "fixed", phi = 2.5,sigmasq.prior = "sc.inv.chisq",
df.sigmasq = 5, sigmasq = 0.5)
ap <- krige.bayes(s100,prior=prior2.b9, model=model.control(trend.d="1st") )
prior2.b9 <- prior.control(beta.prior = "normal", beta = 0,beta.var =1,phi.prior = "fixed", phi = 2.5,sigmasq.prior ="sc.inv.chisq",df.sigmasq = 5, sigmasq = 0.5)
ap <- krige.bayes(s100,prior=prior2.b9)
par(op)</testonly>