ramps (version 0.6-12)

georamps: Bayesian Geostatistical Model Fitting with RAMPS

Description

General function for fitting Bayesian geostatistical models using the reparameterized and marginalized posterior sampling (RAMPS) algorithm of Yan et al. (2007).

Usage

georamps(fixed, random, correlation, data, subset, weights,
            variance = list(fixed = ~ 1, random = ~ 1, spatial = ~ 1),
            aggregate = list(grid = NULL, blockid = ""), kmat = NULL,
            control = ramps.control(...), contrasts = NULL, ...)

Arguments

fixed

two-sided linear "formula" object describing the main effects in the mean structure of the model, with the response on the left of a ~ operator and the terms, separated by + operators, on the right.

random

optional one-sided formula of the form ~ 1 | g, specifying random intercepts for groups defined by the factor g. Several grouping variables may be simultaneously specified, separated by the * operator, as in ~

correlation

'corSpatial' object describing the spatial correlation structure. See the corClasses documentation for a listing of the available structures.

data

optional data frame containing the variables named in fixed, random, correlation, weights, variance, and subset.

subset

optional expression indicating the subset of rows in data that should be used in the fit. This can be a logical vector, or a numerical vector indicating which observation numbers are to be included, or a character vector of the row names to b

weights

optional numerical vector of measurement error variance (inverse) weights to be used in the fitting process. Defaults to a value of 1 for point-source measurements and the number of grid points for areal measurements (see the aggregate argume

variance

optional list of one-sided formulas, each of the form ~ g where g defines a grouping factor for the following elements: fixed for measurement error variances; random for random effects error variances; a

aggregate

optional list of elements: grid a data frame of coordinates to use for Monte Carlo integration over geographic blocks at which areal measurements are available; and blockid a character string specifying the column by which to mer

kmat

optional $n \times s$ design matrix for mapping spatial sites to outcome responses, where $n$ is the number of responses and $s$ the number of unique sites. Unique sites are ordered first according to those supplied to the data argument and

control

list of parameters for controlling the fitting process. See the ramps.control documentation for details.

contrasts

optional list. See the contrasts.arg of model.matrix.

...

further arguments passed to or from other methods.

Value

• An object of class 'ramps' containing the following elements:
• params
  'mcmc' object of monitored model parameters with variable labels in the column names and MCMC iteration numbers in the row names.
• z
  'mcmc' object of monitored latent spatial parameters with variable labels in the column names and MCMC iteration numbers in the row names.
• loglik
  vector of data log-likelihood values at each MCMC iteration.
• evals
  vector of slice sampler evaluations at each MCMC iteration.
• call
  the matched function call to georamps.
• y
  response vector.
• xmat
  design matrix for the main effects.
• terms
  the 'terms' object for xmat.
• xlevels
  list of the factor levels for xmat.
• etype
  grouping factor for the measurement error variances.
• weights
  weights used in the fitting process.
• kmat
  matrix for mapping the spatial parameters to the observed data.
• correlation
  specified 'corSpatial' object for the spatial correlation structure.
• coords
  matrix of unique coordinates for the measurement and grid sites.
• ztype
  grouping factor for the spatial variances.
• wmat
  matrix for mapping the random effects to the observed data.
• retype
  grouping factor for the random effects variances.
• control
  a list of control parameters used in the fitting process.

References

Yan, J., Cowles, M.K., Wang, S., and Armstrong, M. (2007) Parallelizing MCMC for Bayesian Spatiotemporal Geostatistical Models, Statistics and Computing, 17(4), 323-335. Smith, B. J., Yan, J., and Cowles, M. K. (2008) Unified Geostatistical Modeling for Data Fusion and Spatial Heteroskedasticity with R Package ramps, Journal of Statistical Software, 25(10), 1-21.

See Also

corClasses, ramps.control, mcmc, DIC.ramps, plot.ramps, predict.ramps, summary.ramps, window.ramps

Examples

## Load the included uranium datasets for use in this example
data(NURE)

## Geostatistical analysis of areal measurements
NURE.ctrl1 <- ramps.control(
   iter = 25,
   beta = param(0, "flat"),
   sigma2.e = param(1, "invgamma", shape = 2.0, scale = 0.1, tuning = 0.75),
   phi = param(10, "uniform", min = 0, max = 35, tuning = 0.50),
   sigma2.z = param(1, "invgamma", shape = 2.0, scale = 0.1)
)

NURE.fit1 <- georamps(log(ppm) ~ 1,
   correlation = corRExp(form = ~ lon + lat, metric = "haversine"),
   weights = area,
   data = NURE,
   subset = (measurement == 1),
   aggregate = list(grid = NURE.grid, blockid = "id"),
   control = NURE.ctrl1
)
print(NURE.fit1)
summary(NURE.fit1)


## Analysis of point-source measurements
NURE.ctrl2 <- ramps.control(
   iter = 25,
   beta = param(0, "flat"),
   sigma2.e = param(1, "invgamma", shape = 2.0, scale = 0.1, tuning = 0.75),
   phi = param(10, "uniform", min = 0, max = 35, tuning = 0.5),
   sigma2.z = param(1, "invgamma", shape = 2.0, scale = 0.1)
)

NURE.fit2 <- georamps(log(ppm) ~ 1,
   correlation = corRExp(form = ~ lon + lat, metric = "haversine"),
   data = NURE,
   subset = (measurement == 2),
   control = NURE.ctrl2
)
print(NURE.fit2)
summary(NURE.fit2)


## Joint analysis of areal and point-source measurements with
## prediction only at grid sites
NURE.ctrl <- ramps.control(
   iter = 25,
   beta = param(rep(0, 2), "flat"),
   sigma2.e = param(rep(1, 2), "invgamma", shape = 2.0, scale = 0.1, tuning = 0.75),
   phi = param(10, "uniform", min = 0, max = 35, tuning = 0.5),
   sigma2.z = param(1, "invgamma", shape = 2.0, scale = 0.1),
   z.monitor = NURE.grid
)

NURE.fit <- georamps(log(ppm) ~ factor(measurement) - 1,
   correlation = corRExp(form = ~ lon + lat, metric = "haversine"),
   variance = list(fixed = ~ measurement),
   weights = area * (measurement == 1) + (measurement == 2),
   data = NURE,
   aggregate = list(grid = NURE.grid, blockid = "id"),
   control = NURE.ctrl
)
print(NURE.fit)
summary(NURE.fit)


## Discard initial 5 MCMC samples as a burn-in sequence
fit <- window(NURE.fit, iter = 6:25)
print(fit)
summary(fit)

## Deviance Information Criterion
DIC(fit)

## Prediction at unmeasured sites
ct <- map("state", "connecticut", plot = FALSE)
lon <- seq(min(ct$x, na.rm = TRUE), max(ct$x, na.rm = TRUE), length = 20)
lat <- seq(min(ct$y, na.rm = TRUE), max(ct$y, na.rm = TRUE), length = 15)
grid <- expand.grid(lon, lat)

newsites <- data.frame(lon = grid[,1], lat = grid[,2],
                       measurement = 1)
pred <- predict(fit, newsites)

plot(pred, func = function(x) exp(mean(x)),
     database = "state", regions = "connecticut",
     resolution = c(200, 150), bw = 5,
     main = "Posterior Mean",
     legend.args = list(text = "ppm", side = 3, line = 1))

plot(pred, func = function(x) exp(sd(x)),
     database = "state", regions = "connecticut",
     resolution = c(200, 150), bw = 5,
     main = "Posterior Standard Deviation",
     legend.args = list(text = "ppm", side = 3, line = 1))