ExceedanceTools.

Description

A package to create confidence or credible regions for the exceedance regions/excursion sets of spatial data.

Colorado precipitation data

Description

Data related to Colorado precipitation in May 1997. Taken from https://www.image.ucar.edu/Data/US.monthly.met/. Data is contained in a list with components odata (containing a transformed precipitation variable) and ocoords containing the longitude and latitude of the associated sites.

Usage

data(colorado)

Format

A list.

Author(s)

Joshua French

Source

National Center for Atmospheric Research

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

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:

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")

Create grid of locations.

Description

create.pgrid creates a grid of locations from the boundaries of domain and other information.

Usage

create.pgrid(
  xmin,
  xmax,
  ymin,
  ymax,
  nx,
  ny,
  midpoints = FALSE,
  poly.coords = NULL
)

Arguments

Details

The key argument in the function midpoints. If this is TRUE, it is assumed that the boundaries of the spatial domain correspond to the midpoints of the cell/pixel in the grid. Otherwise, it is assumed that the boundaries correspond to the actual borders of the region of interest. If poly.coords is supplied, the grid returned is the grid of midpoints contained in the convex hull of poly.coords.

Value

Returns an object of class pgrid with the following components:

Author(s)

Joshua French

Examples

pgrida <- create.pgrid(0, 1, 0, 1, nx = 50, ny = 50, midpoints = FALSE)
pgridb <- create.pgrid(.01, .99, .01, .99, nx = 50, ny = 50, midpoints = TRUE)

Create grid of locations.

Description

create.pgrid2 creates a grid of locations fusing vectors of x and y coordinates.

Usage

create.pgrid2(xgrid, ygrid, midpoints = FALSE, poly.coords = NULL)

Arguments

Details

The key argument in the function midpoints. If this is TRUE, it is assumed that the boundaries of the spatial domain correspond to the midpoints of the cell/pixel in the grid. Otherwise, it is assumed that the boundaries correspond to the actual borders of the region of interest. If poly.coords is supplied, the grid returned is the grid of midpoints contained in the convex hull of poly.coords.

Value

Returns an object of class pgrid with the following components:

Author(s)

Joshua French

Examples

seq1 = seq(0, 1, len = 101)
pgrida <- create.pgrid2(seq1, seq1, midpoint = FALSE)
seq2 = seq(.005, .995, len = 100)
pgridb <- create.pgrid2(seq2, seq2, midpoint = TRUE)
# pgrids produced match
range(pgrida$pgrid - pgridb$pgrid)

Return confidence region

Description

exceedance.ci returns a confidence set for an exceedance region or contour line.

Usage

exceedance.ci(statistic.sim.obj, conf.level = 0.95, type = "null")

Arguments

Value

Returns a numeric vector with the set of pixels comprising the null or rejection region related to statistic.sim.obj.

Author(s)

Joshua French

Examples

library(SpatialTools)

# Example for exceedance regions

set.seed(10)
# Load data
data(sdata)
# Create prediction grid
pgrid <- create.pgrid(0, 1, 0, 1, nx = 26, ny = 26)
pcoords <- pgrid$pgrid
# Create design matrices
coords = cbind(sdata$x1, sdata$x2)
X <- cbind(1, coords)
Xp <- cbind(1, pcoords)

# Generate covariance matrices V, Vp, Vop using appropriate parameters for 
# observed data and responses to be predicted
spcov <- cov.sp(coords = coords, sp.type = "exponential", 
 sp.par = c(1, 1.5), error.var = 1/3, finescale.var = 0, pcoords = pcoords)

# Predict responses at pgrid locations
krige.obj <- krige.uk(y = as.vector(sdata$y), V = spcov$V, Vp = spcov$Vp, 
 Vop = spcov$Vop, X = X, Xp = Xp, nsim = 100, 
 Ve.diag = rep(1/3, length(sdata$y)) , method = "chol")
                
# Simulate distribution of test statistic for different alternatives
statistic.sim.obj.less <- statistic.sim(krige.obj = krige.obj, level = 5, 
 alternative = "less")
statistic.sim.obj.greater <- statistic.sim(krige.obj = krige.obj, level = 5,
 alternative = "greater")
# Construct null and rejection sets for two scenarios
n90 <- exceedance.ci(statistic.sim.obj.less, conf.level = .90, type = "null")
r90 <- exceedance.ci(statistic.sim.obj.greater,conf.level = .90, type = "rejection")       
# Plot results
plot(pgrid, n90, col="blue", add = FALSE, xlab = "x", ylab = "y")
plot(pgrid, r90, col="orange", add = TRUE)
legend("bottomleft", 
 legend = c("contains true exceedance region with 90 percent confidence", 
   "is contained in true exceedance region with 90 percent confidence"),
   col = c("blue", "orange"), lwd = 10)  

Plots pgrid object.

Description

plot.pgrid plots a grid of pixels based on a pgrid object.

Usage

## S3 method for class 'pgrid'
plot(x, set, col = "gray", add = FALSE, type = "confidence", ...)

Arguments

Details

If a vector of pixel indices is supplied to set, then those pixels will be colored col by this function and the type argument has no effect. On the other hand, if the set argument is of class confreg, then the function digs in to display either the confidence or complement set in the confreg object. In that case, type is used to decide which set to display.

Value

This function does not return anything; it only creates a new plot or modifies an existing plot.

Author(s)

Joshua French

Examples

library(SpatialTools)

# Example for exceedance regions

set.seed(10)
# Load data
data(sdata)
# Create prediction grid
pgrid <- create.pgrid(0, 1, 0, 1, nx = 26, ny = 26)
pcoords <- pgrid$pgrid
# Create design matrices
coords = cbind(sdata$x1, sdata$x2)
X <- cbind(1, coords)
Xp <- cbind(1, pcoords)

# Generate covariance matrices V, Vp, Vop using appropriate parameters for 
# observed data and responses to be predicted
spcov <- cov.sp(coords = coords, sp.type = "exponential", 
 sp.par = c(1, 1.5), error.var = 1/3, finescale.var = 0, pcoords = pcoords)

# Predict responses at pgrid locations
krige.obj <- krige.uk(y = as.vector(sdata$y), V = spcov$V, Vp = spcov$Vp, 
 Vop = spcov$Vop, X = X, Xp = Xp, nsim = 100, 
 Ve.diag = rep(1/3, length(sdata$y)) , method = "chol")
                
# Simulate distribution of test statistic for different alternatives
statistic.sim.obj.less <- statistic.sim(krige.obj = krige.obj, level = 5, 
 alternative = "less")
statistic.sim.obj.greater <- statistic.sim(krige.obj = krige.obj, 
 level = 5, alternative = "greater")
# Construct null and rejection sets for two scenarios
n90 <- exceedance.ci(statistic.sim.obj.less, conf.level = .90, 
 type = "null")
r90 <- exceedance.ci(statistic.sim.obj.greater,conf.level = .90, 
 type = "rejection")       
# Plot results
plot(pgrid, n90, col="blue", add = FALSE, xlab = "x", ylab = "y")
plot(pgrid, r90, col="orange", add = TRUE)
legend("bottomleft", 
 legend = c("contains true exceedance region with 90 percent confidence", 
 "is contained in true exceedance region with 90 percent confidence"),
 col = c("blue", "orange"), lwd = 10)  

Synthetic data

Description

A synthetic data set for use in examples. A 100x3 data frame with vectors x1 and x2 (specifying spatial location) and y, the response.

Usage

data(sdata)

Format

A data frame.

Author(s)

Joshua French

Return critical value of distribution.

Description

statistic.cv returns the critical value of the distribution of the test statistics from statistic.sim based on the specified confidence level. However, it is not recommended for general usage. It is recommedned that the exceedance.ci function be used to automatically create confidence regions.

Usage

statistic.cv(statistic.sim.obj, conf.level = 0.95)

Arguments

Value

Returns the desired critical value.

Author(s)

Joshua French

Examples

library(SpatialTools)

# Example for exceedance regions

set.seed(10)
# Load data
data(sdata)
# Create prediction grid
pgrid <- create.pgrid(0, 1, 0, 1, nx = 26, ny = 26)
pcoords <- pgrid$pgrid
# Create design matrices
coords = cbind(sdata$x1, sdata$x2)
X <- cbind(1, coords)
Xp <- cbind(1, pcoords)

# Generate covariance matrices V, Vp, Vop using appropriate parameters for 
# observed data and responses to be predicted
spcov <- cov.sp(coords = coords, sp.type = "exponential", sp.par = c(1, 1.5),
 error.var = 1/3, finescale.var = 0, pcoords = pcoords)

# Predict responses at pgrid locations
krige.obj <- krige.uk(y = as.vector(sdata$y), V = spcov$V, Vp = spcov$Vp, 
 Vop = spcov$Vop, X = X, Xp = Xp, nsim = 100, 
 Ve.diag = rep(1/3, length(sdata$y)) , method = "chol")
                
# Simulate distribution of test statistic for different alternatives
statistic.sim.obj.less <- statistic.sim(krige.obj = krige.obj, level = 5, 
 alternative = "less")
statistic.sim.obj.greater <- statistic.sim(krige.obj = krige.obj, level = 5,
 alternative = "greater")
# Calculate quantiles of distribution of statistic
q90.less <- statistic.cv(statistic.sim.obj.less, conf.level = .90)
q90.greater <- statistic.cv(statistic.sim.obj.greater, conf.level = .90)

Simulates statistics related to exceedance region.

Description

statistic.sim simulates statistics related to the construction of confidence regions for exceedance sets and contour lines.

Usage

statistic.sim(krige.obj, level, alternative = "less", ...)

Arguments

Details

When alternative = "two.sided", the ... argument must include user.cov (a user-specified covariance function), pgrid (the grid of locations to be predicted, produced by create.pgrid or create.pgrid2), X (the matrix of covariates for the observed data), and any other arguments needed by user.cov. Note that user.cov should take cLcoords as its first argument (a matrix containing the coordinates of contour lines under consideration). Additional arguments to user.cov are passed internally using the ... argument. The user.cov function should return a list with values V (the covariance matrix of the observed data), Vop (the cross-covariance matrix between the observed data and the responses with coordinates in cL), Vp (the covariance matrix of the responses with coordinates in cL), and Xp (the matrix of covariates for the coordinates contained in cL). See the Examples section.

Value

Returns a list with components:

Author(s)

Joshua French

Examples

library(SpatialTools)

# Example for exceedance regions

set.seed(10)
# Load data
data(sdata)
# Create prediction grid
pgrid <- create.pgrid(0, 1, 0, 1, nx = 26, ny = 26)
pcoords <- pgrid$pgrid
# Create design matrices
coords = cbind(sdata$x1, sdata$x2)
X <- cbind(1, coords)
Xp <- cbind(1, pcoords)

# Generate covariance matrices V, Vp, Vop using appropriate parameters for
# observed data and responses to be predicted
spcov <- cov.sp(coords = coords, sp.type = "exponential", sp.par = c(1, 1.5),
 error.var = 1/3, finescale.var = 0, pcoords = pcoords)

# Predict responses at pgrid locations
krige.obj <- krige.uk(y = as.vector(sdata$y), V = spcov$V, Vp = spcov$Vp, 
 Vop = spcov$Vop, X = X, Xp = Xp, nsim = 50, 
 Ve.diag = rep(1/3, length(sdata$y)) , method = "chol")
                
# Simulate distribution of test statistic for different alternatives
statistic.sim.obj.less <- statistic.sim(krige.obj = krige.obj, level = 5,
 alternative = "less")
statistic.sim.obj.greater <- statistic.sim(krige.obj = krige.obj, level = 5,
 alternative = "greater")
# Construct null and rejection sets for two scenarios
n90 <- exceedance.ci(statistic.sim.obj.less, conf.level = .90, type = "null")
r90 <- exceedance.ci(statistic.sim.obj.greater,conf.level = .90, 
 type = "rejection")       
# Plot results
plot(pgrid, n90, col="blue", add = FALSE, xlab = "x", ylab = "y")
plot(pgrid, r90, col="orange", add = TRUE)
legend("bottomleft", 
 legend = c("contains true exceedance region with 90 percent confidence", 
 "is contained in true exceedance region with 90 percent confidence"),
 col = c("blue", "orange"), lwd = 10)  

# Example for level curves
data(colorado)
ocoords <- colorado$ocoords
odata <- colorado$odata    

# Set up example
nsim <- 50
u <- log(16)
np <- 26
conf.level <- 0.90   
x.min <- min(ocoords[,1])
x.max <- max(ocoords[,1])
y.min <- min(ocoords[,2])
y.max <- max(ocoords[,2])   

#pixelize the domain
pgrid <- create.pgrid(x.min, x.max, y.min, y.max, nx = np, ny = np)
pcoords <- pgrid$pgrid; upx <- pgrid$upx; upy <- pgrid$upy
names(pcoords) <- c("lon", "lat")    

# Set up covariates matrices
X <- cbind(1, ocoords)
Xp <- cbind(1, pcoords)           

# Estimate covariance parameters
cov.est <- maxlik.cov.sp(X, odata, sp.type = "exponential", range.par = 1.12,
 error.ratio = 0.01, reml = TRUE, coords = ocoords)

# Create covariance matrices
myCov <- cov.sp(coords = ocoords, sp.type = "exponential", 
 sp.par = cov.est$sp.par, error.var = cov.est$error.var, pcoords = pcoords)

# Krige and do conditional simulation 
krige.obj <- krige.uk(y = odata, V = myCov$V, Vp = myCov$Vp, Vop = myCov$Vop, 
 X = X, Xp = Xp, nsim = nsim, Ve.diag = rep(cov.est$error.var, 
 length(odata)))

# Create user covariance function for simulating statistic for confidence 
# regions
user.cov <- function(cLcoords,...)
{
   arglist <- list(...)
   coords <- arglist$coords
   sp.type <- arglist$sp.type
   sp.par <- arglist$sp.par
   V <- arglist$V
   out <- list(V = arglist$V,
               Vp = sp.par[1] * exp(-dist1(cLcoords)/sp.par[2]),
               Vop = sp.par[1] * exp(-dist2(coords, cLcoords)/sp.par[2]))
   out$Xp <- cbind(1, cLcoords)
   return(out)
}

# Simulation statistic for confidence regions
statistic.sim.obj <- statistic.sim(krige.obj = krige.obj, level = u, 
 alternative = "two.sided", user.cov = user.cov, y = odata, pgrid = pgrid, 
 X = X, coords = ocoords, pcoords = pcoords, V = myCov$V, 
 sp.type = "exponential", sp.par = cov.est$sp.par)

# Create 90% confidence region
n90 <- exceedance.ci(statistic.sim.obj, conf.level = conf.level, 
 type = "null")
# Get estimated contour lines
cL <- contourLines(pgrid$upx, pgrid$upy, matrix(krige.obj$pred, nrow = np), 
 level = u)

# Plot results
plot(ocoords, xlab = "longitude", ylab = "latitude", type = "n", 
 cex.lab = 1.5, cex.axis = 1.5)
plot(pgrid, n90, col = "grey", add = TRUE)
plot.contourLines(cL, col="black", lwd=2, lty = 2, add = TRUE)