georamps(fixed, random, correlation, data, subset, weights,
variance = list(fixed = ~ 1, random = ~ 1, spatial = ~ 1),
aggregate = list(grid = NULL, blockid = ""),
control = ramps.control(...), contrasts = NULL, ...)"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.~ 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 ~ 'corSpatial' object describing the spatial correlation structure. See the corClasses documentation for a listing of the available structures.fixed, random, correlation, weights, variance, and subset.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 baggregate argument below).~ g where g defines a grouping factor for the following elements: fixed for measurement error variances; random for random effects error variances; agrid 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 merramps.control documentation for details.contrasts.arg of model.matrix.'ramps' containing the following elements:'mcmc'"} object of monitored model parameters with variable labels in the column names and MCMC iteration numbers in the row names.}
item{z}{code{'mcmc'} object of monitored latent spatial parameters with variable labels in the column names and MCMC iteration numbers in the row names.}
item{loglik}{vector of data log-likelihood values at each MCMC iteration.}
item{evals}{vector of slice sampler evaluations at each MCMC iteration.}
item{call}{the matched function call to code{georamps}.}
item{y}{response vector.}
item{xmat}{design matrix for the main effects.}
item{terms}{the code{'terms'} object for code{xmat}.}
item{xlevels}{list of the factor levels for code{xmat}.}
item{etype}{grouping factor for the measurement error variances.}
item{weights}{weights used in the fitting process.}
item{kmat}{matrix for mapping the spatial parameters to the observed data.}
item{correlation}{specified code{'corSpatial'} object for the spatial correlation structure.}
item{coords}{matrix of unique coordinates for the measurement and grid sites.}
item{ztype}{grouping factor for the spatial variances.}
item{wmat}{matrix for mapping the random effects to the observed data.}
item{retype}{grouping factor for the random effects variances.}
item{control}{a list of control parameters used in the fitting process.}
}
references{
Yan, J., Cowles, M.K., Wang, S., and Armstrong, M. (2007) dQuote{Parallelizing MCMC for Bayesian Spatiotemporal Geostatistical Models}, emph{Statistics and Computing}, 17(4), 323-335.
Smith, B. J., Yan, J., and Cowles, M. K. (2008) dQuote{Unified Geostatistical Modeling for Data Fusion and Spatial Heteroskedasticity with R Package ramps}, emph{Journal of Statistical Software}, in press.
}
author{
Brian Smith email{brian-j-smith@uiowa.edu}, Jun Yan email{jun.yan@uconn.edu}, and Kate Cowles email{kate-cowles@uiowa.edu}}
seealso{
code{corClasses},
code{ramps.control},
code{mcmc},
code{DIC.ramps},
code{plot.ramps},
code{predict.ramps},
code{summary.ramps},
code{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 = 0.1,
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 = 0.1,
main = "Posterior Standard Deviation",
legend.args = list(text = "ppm", side = 3, line = 1))
}
keyword{models}