Calculate continuous domain excursion and credible contour sets

Calculates continuous domain excursion and credible contour sets

Usage

continuous(
  ex,
  geometry,
  alpha,
  method = c("log", "linear", "step"),
  output = c("sp", "fm", "inla"),
  subdivisions = 1,
  calc.credible = TRUE
)

Arguments

ex: An excurobj object generated by a call to excursions(), contourmap(), or excursions.regions().
geometry: Specification of the lattice or triangulation geometry of the input. One of list(x, y), list(loc, dims), fm_lattice_2d, inla.mesh.lattice, fm_mesh_2d, or inla.mesh, where x and y are vectors, loc is a two-column matrix of coordinates, and dims is the lattice size vector. The first three versions are all treated topologically as lattices, and the lattice boxes are assumed convex.
alpha: The target error probability. A warning is given if it is detected that the information ex isn't sufficient for the given alpha. Defaults to the value used when calculating ex.
method: The spatial probability interpolation transformation method to use. One of log, linear, or step. For log, the probabilities are interpolated linearly in the transformed scale. For step, a conservative step function is used.
output: Specifies what type of object should be generated. sp gives a SpatialPolygons object, and fm or inla gives a fm_segm object.
subdivisions: The number of mesh triangle subdivisions to perform for the interpolation of the excursions or contour function. 0 is no subdivision. The setting has a small effect on the evaluation of P0 for the log method (higher values giving higher accuracy) but the main effect is on the visual appearance of the interpolation. Default=1.
calc.credible: Logical, if TRUE (default), calculate credible contour region objects in addition to avoidance sets.

Value

A list with elements

M: SpatialPolygons or fm_segm object. The subsets are tagged, so that credible regions are tagged "-1", and regions between levels are tagged as.character(0:nlevels).
F: Interpolated F function.
G: Contour and inter-level set indices for the interpolation.
F.geometry: Mesh geometry for the interpolation.
P0: P0 measure based on interpolated F function (only for contourmap input).

For input from excursions.regions() or excursions.regions.inla(), the continuous domain region of each region is computed from the excursion function of that region, and the regions are tagged by their number, as.character(1:n). F is then a matrix with one column for each region. The excursion functions are only computed for the nodes of each region and its neighbours, so alpha can not be larger than the value used for computing the regions.

References

Bolin, D. and Lindgren, F. (2017) Quantifying the uncertainty of contour maps, Journal of Computational and Graphical Statistics, vol 26, no 3, pp 513-524.

Bolin, D. and Lindgren, F. (2018), Calculating Probabilistic Excursion Sets and Related Quantities Using excursions, Journal of Statistical Software, vol 86, no 1, pp 1-20.

Author

Finn Lindgren finn.lindgren@gmail.com

Examples

if (FALSE) { # \dontrun{
if (require("fmesher") &&
  require("sp")) {
  # Generate mesh and SPDE model
  n.lattice <- 10 # Increase for more interesting, but slower, examples
  x <- seq(from = 0, to = 10, length.out = n.lattice)
  lattice <- fm_lattice_2d(x = x, y = x)
  mesh <- fm_rcdt_2d_inla(lattice = lattice, extend = FALSE, refine = FALSE)

  # Generate an artificial sample
  sigma2.e <- 0.1
  n.obs <- 100
  obs.loc <- cbind(
    runif(n.obs) * diff(range(x)) + min(x),
    runif(n.obs) * diff(range(x)) + min(x)
  )
  Q <- fm_matern_precision(mesh, alpha = 2, rho = 3, sigma = 1)
  x <- fm_sample(n = 1, Q = Q)
  A <- fm_basis(mesh, loc = obs.loc)
  Y <- as.vector(A %*% x + rnorm(n.obs) * sqrt(sigma2.e))

  ## Calculate posterior
  Q.post <- (Q + (t(A) %*% A) / sigma2.e)
  mu.post <- as.vector(solve(Q.post, (t(A) %*% Y) / sigma2.e))
  vars.post <- excursions.variances(chol(Q.post), max.threads = 1)

  ## Calculate contour map with two levels
  map <- contourmap(
    n.levels = 2, mu = mu.post, Q = Q.post,
    alpha = 0.1, F.limit = 0.1, max.threads = 1
  )

  ## Calculate the continuous representation
  sets <- continuous(map, mesh, alpha = 0.1)

  ## Plot the results
  reo <- mesh$idx$lattice
  cols <- contourmap.colors(map,
    col = heat.colors(100, 1, rev = TRUE),
    credible.col = grey(0.5, 1)
  )
  names(cols) <- as.character(-1:2)

  par(mfrow = c(2, 2))
  image(matrix(mu.post[reo], n.lattice, n.lattice),
    main = "mean", axes = FALSE, asp = 1
  )
  image(matrix(sqrt(vars.post[reo]), n.lattice, n.lattice),
    main = "sd", axes = FALSE, asp = 1
  )
  image(matrix(map$M[reo], n.lattice, n.lattice),
    col = cols, axes = FALSE, asp = 1
  )
  idx.M <- setdiff(names(sets$M), "-1")
  plot(sets$M[idx.M], col = cols[idx.M])
}
} # }