| Title: | Confidence/Credible Regions for Exceedance Sets and Contour Lines |
|---|---|
| Description: | Provides methods for constructing confidence or credible regions for exceedance sets and contour lines. |
| Authors: | Joshua French |
| Maintainer: | Joshua French <[email protected]> |
| License: | GPL (>=2) |
| Version: | 1.3.6 |
| Built: | 2024-10-16 02:44:52 UTC |
| Source: | https://github.com/jfrench/exceedancetools |
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.
data(colorado)data(colorado)
A list.
Joshua French
National Center for Atmospheric Research
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.
confreg( obj, level, statistic = NULL, conf.level = 0.95, direction = ">", type = "o", method = "test", greedy = FALSE )confreg( obj, level, statistic = NULL, conf.level = 0.95, direction = ">", type = "o", method = "test", greedy = FALSE )
obj |
An object of the appropriate type ( |
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. |
type |
|
method |
|
greedy |
Only applicable for the direct construction method. Default is |
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.
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. |
Joshua French
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.
# 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")# 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.pgrid creates a grid of locations from the boundaries of domain and other information.
create.pgrid( xmin, xmax, ymin, ymax, nx, ny, midpoints = FALSE, poly.coords = NULL )create.pgrid( xmin, xmax, ymin, ymax, nx, ny, midpoints = FALSE, poly.coords = NULL )
xmin |
The minimum value of the boundary of the x coordinates of the spatial domain. |
xmax |
The maximum value of the boundary of the x coordinates of the spatial domain. |
ymin |
The minimum value of the boundary of the y coordinates of the spatial domain. |
ymax |
The maximum value of the boundary of the y coordinates of the spatial domain. |
nx |
The number of gridpoints/cells/pixels in the x direction. |
ny |
The number of gridpoints/cells/pixels in the y direction. |
midpoints |
A logical value ( |
poly.coords |
An |
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.
Returns an object of class pgrid with the following components:
pgrid |
An |
m |
The number of rows in pgrid. |
p.in.grid |
A vector of 0s and 1s indicating whether the midpoint of each pixel is in the convex hull of |
ubx |
The pixel boundaries in the x direction. |
uby |
The pixel boundaries in the y direction. |
upx |
The pixel midpoints in the x direction. |
upy |
The pixel midpoints in the y direction. |
Joshua French
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)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.pgrid2 creates a grid of locations fusing vectors of x and y coordinates.
create.pgrid2(xgrid, ygrid, midpoints = FALSE, poly.coords = NULL)create.pgrid2(xgrid, ygrid, midpoints = FALSE, poly.coords = NULL)
xgrid |
A vector of locations in the x direction. |
ygrid |
A vector of location in the y direction. |
midpoints |
A logical value ( |
poly.coords |
An |
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.
Returns an object of class pgrid with the following components:
pgrid |
An |
m |
The number of rows in pgrid. |
p.in.grid |
A vector of 0s and 1s indicating whether the midpoint of each pixel is in the convex hull of |
ubx |
The pixel boundaries in the x-direction. |
uby |
The pixel boundaries in the y-direction. |
upx |
The pixel midpoints in the x-direction. |
upy |
The pixel midpoints in the y-direction. |
Joshua French
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)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)
exceedance.ci returns a confidence set for an exceedance region or contour line.
exceedance.ci(statistic.sim.obj, conf.level = 0.95, type = "null")exceedance.ci(statistic.sim.obj, conf.level = 0.95, type = "null")
statistic.sim.obj |
An object returned from the |
conf.level |
The desired confidence level of the confidence region. |
type |
Whether the function should return the null region or rejection region of exceedance confidence region Options are |
Returns a numeric vector with the set of pixels comprising the null or rejection region related to statistic.sim.obj.
Joshua French
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)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)
A package to create confidence or credible regions for the exceedance regions/excursion sets of spatial data.
pgrid object.plot.pgrid plots a grid of pixels based on a pgrid object.
## S3 method for class 'pgrid' plot(x, set, col = "gray", add = FALSE, type = "confidence", ...)## S3 method for class 'pgrid' plot(x, set, col = "gray", add = FALSE, type = "confidence", ...)
x |
An |
set |
A vector which contains the indices of the pixels/cells that should be plotted. OR a |
col |
The color of the plotted pixels. |
add |
A logical value indicating whether the pixels should be added to an existing plot ( |
type |
The type of set of plot if |
... |
Additional arguments that will be passed to the |
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.
This function does not return anything; it only creates a new plot or modifies an existing plot.
Joshua French
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)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)
A synthetic data set for use in examples. A 100x3 data frame with vectors x1 and x2 (specifying spatial location) and y, the response.
data(sdata)data(sdata)
A data frame.
Joshua French
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.
statistic.cv(statistic.sim.obj, conf.level = 0.95)statistic.cv(statistic.sim.obj, conf.level = 0.95)
statistic.sim.obj |
An object returned from the |
conf.level |
The desired confidence level of the confidence interval we want to construct. |
Returns the desired critical value.
Joshua French
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)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)
statistic.sim simulates statistics related to the
construction of confidence regions for exceedance sets
and contour lines.
statistic.sim(krige.obj, level, alternative = "less", ...)statistic.sim(krige.obj, level, alternative = "less", ...)
krige.obj |
An object from the function
|
level |
The threshold/exceedance level under consideration. |
alternative |
Indicates the type of exceedance
region or level curve under consideration. For
exceedances above a threshold, use ( |
... |
Additional arguments when |
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.
Returns a list with components:
statistic |
A vector with the observed values of the test statistic. |
statistic.sim |
A vector with the observed values of the test statistic. |
alternative |
The alternative hypothesis provided
to |
level |
The threshold level under consideration. |
Joshua French
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)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)