confreg {ExceedanceTools}

Construct confidence regions for exceedance (excursion) sets.

Description

confreg constructs confidence regions for the exceedance (excursions) sets of geostatistical processes. These will actually be credible regions if obj contains samples from the joint posterior predictive distribution in a Bayesian setting.

Usage

confreg(
  obj,
  level,
  statistic = NULL,
  conf.level = 0.95,
  direction = ">",
  type = "o",
  method = "test",
  greedy = FALSE
)

Arguments

obj

An object of the appropriate type (matrix, krigeConditionalSample, or jointPredictiveSample. See Details.

level

The threshold level for the exceedance region.

statistic

The statistic used in constructing the confidence region. Should be a vector containing a value for each location

conf.level

The confidence level of the confidence region. Default is 0.95.

direction

The direction of the exceedance region. ">" indicates the exceedance region is values above a threshold, while "<" indicates values below a threshold.

type

"o" indicates on outer confidence region while "i" indicates in inner confidence region.

method

"test" indicates a testing-based method, while "direct" indicates a direct method using joint probabilities.

greedy

Only applicable for the direct construction method. Default is FALSE. If TRUE, then grid cells are added to the confidence region using a greedy algorithm based on joint probability.

Details

obj can be an object of class matrix, krigeConditionalSample, or jointPredictiveSample. If obj is a matrix, then it should have m rows and nsim columns. In that case, each row of obj corresponds to a sample from the conditional distribution of the response conditional on the observed data. Each row represents a different location. Generally, these locations are assumed to be on a grid spanning the spatial domain of interest. A krigeConditionalSample object can be obtained using the krige.sk, krige.ok, or krige.uk functions in the SpatialTools package. In these functions, the nsim argument must be greater than 0, and indicates the number of samples used to construct the confidence region. A jointPredictiveSample object can be obtained using the spLMPredictJoint function in the SpatialTools package. Since this is in the context of Bayesian statistics, the function actually produces credible region.

If statistic is supplied for the direct construction procedure, then the locations are ordered by marginal probability and then the statistic. statistic should be a vector of length m, where m is the number of prediction locations at which samples were drawn for in obj.

If type == "o", then an outer credible region is constructed. The outer credible region should entirely contain the true exceedanace region with the specified posterior probability. If type == "i", then an inner credible region is constructed. The inner confidence region should be entirely contained within the true exceedanace region with specified posterior probability.

Value

Returns an object of class confreg with the following components:

confidence

The sites included in the confidence region.

complement

The complement of the confidence region.

Author(s)

Joshua French

References

Joshua P. French and Stephan R. Sain (2013). Spatio-temporal exceedance locations and confidence regions. Annals of Applied Statistics. 7(3):1421-1449.

French, J. P. (2014), Confidence regions for the level curves of spatial data, Environmetrics, 25, pages 498–512, DOI: 10.1002/env.2295

French, J. P., and Hoeting, J. A. (2016) Credible regions for exceedance sets of geostatistical data. Environmetrics, 27: 4–14. doi: 10.1002/env.2371.

Examples

# Set parameters
n <- 100
mygrid = create.pgrid(0, 1, 0, 1, nx = 5, ny = 4)
n.samples <- 10
burnin.start <- 1
sigmasq <- 1
tausq <- 0.0
phi <- 1
cov.model <- "exponential"
n.report <- 5

# Generate coordinates
coords <- matrix(runif(2 * n), ncol = 2) 
pcoords <- mygrid$pgrid
# Construct design matrices
X <- as.matrix(cbind(1, coords))
Xp <- cbind(1, pcoords)

# Specify priors
starting <- list("phi" = phi, "sigma.sq"= sigmasq, "tau.sq" = tausq)
tuning <- list("phi"=0.1, "sigma.sq"=0.1, "tau.sq"=0.1)
priors.1 <- list("beta.Norm"=list(c(1, 2, 1), diag(100, 3)), "phi.Unif"=c(0.00001, 10), 
 "sigma.sq.IG"=c(1, 1))

# Generate data
library(SpatialTools)
B <- rnorm(3, c(1, 2, 1), sd = 10)
phi <- runif(1, 0, 10)
sigmasq <- 1/rgamma(1, 1, 1)
V <- simple.cov.sp(D = dist1(coords), cov.model, c(sigmasq, 1/phi), error.var = tausq, 
 smoothness = nu, finescale.var = 0)
y <- X %*% B + rmvnorm(1, rep(0, n), V) + rnorm(n, 0, sqrt(tausq))

# Create spLM object
library(spBayes)
m1 <- spBayes::spLM(y ~ X - 1, coords = coords, starting = starting, tuning = tuning, 
 priors = priors.1, cov.model = cov.model, n.samples = n.samples, verbose = FALSE,
 n.report = n.report)

# Sample from joint posterior predictive distribution
y1 <- spLMPredictJoint(m1, pred.coords = pcoords, pred.covars = Xp, 
 start = burnin.start, verbose = FALSE, method = "chol")
 u = quantile(y, .5)
myfun = function(x)
{
    (mean(x) - u)/sd(x)
}

myfun2 = function(x)
{
 mean(x > u)
}

stat1 = apply(y1, 1, myfun)
stat2 = apply(y1, 1, myfun2)

myconf = confreg(y1, level = u, statistic = NULL, direction = ">", type = "o", method = "direct")
myconf2 = confreg(y1, level = u, statistic = stat1, direction = ">", type = "o")
myconf3 = confreg(y1, level = u, statistic = stat2, direction = ">", type = "o")