statistic.sim: Simulates statistics related to exceedance region.

In jpfrench81/ExceedanceTools: Confidence/Credible Regions for Exceedance Sets and Contour Lines

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

krige.obj	An object from the function krige.uk in the SpatialTools package.
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 (alternative = "less"). For exceedances below a threshold, use (alternative = "greater"). For contour lines, use (alternative = "two.sided"). Defaults to "less".
...	Additional arguments when alternative = "two.sided". See Details.

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:

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 statistic.sim.
level	The threshold level under consideration.

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)

jpfrench81/ExceedanceTools documentation built on Aug. 23, 2023, 5:57 a.m.