Calculating probabilistic excursion sets and related quantities using excursions

The R software package excursions contains methods for calculating probabilistic excursion sets, contour credible regions, and simultaneous confidence bands for latent Gaussian stochastic processes and fields. It also contains methods for uncertainty quantification of contour maps and computation of Gaussian integrals. This article describes the theoretical and computational methods used in the package. The main functions of the package are introduced and two examples illustrate how the package can be used.


Introduction
The ability to find regions where a stochastic process exceeds a certain level, or is significantly different from some reference level, is important in several areas of application. Examples in geosciences include studies of air pollution (Cameletti, Lindgren, Simpson, and Rue 2013), temperature (Furrer, Knutti, Sain, Nychka, and A. 2007), precipitation (Sain, Furrer, and Cressie 2011), and vegetation (Eklundh and Olsson 2003;Bolin, Lindström, Eklundh, and Lindgren 2009), and similar problems can be found in a wide range of scientific fields including brain imaging (Marchini and Presanis 2003) and astrophysics (Beaky, Scherrer, and Villumsen 1992). A related problem is uncertainty quantification of contour curves and more generally of contour maps, which are often used to display estimates of continuous surfaces. The number of contours used in a contour map should typically reflect the uncertainty in the estimate, since one should be allowed to draw many contours if the uncertainty of the estimated surface is low and fewer contours if the uncertainty is high. The ability to quantify the uncertainty in the contour map is important if one should be able to choose the number of contours in a rigorous way.
The R (R Core Team 2013) software package excursions (Bolin and Lindgren 2016a) contains functions for solving these problems for latent Gaussian models (LGMs), which is a large model class that is widely used in applications (see e.g., Rue, Martino, and Chopin 2009). The computational methods are based on the theory introduced by Lindgren (2015, 2016b) and are especially well-suited for models where the latent field has Markov properties. Solving the problems involves compting high-dimensional Gaussian integrals, which can be done more efficiently if Markov properties can be utilized. With the ability to efficiently compute Gaussian integrals, one can also compute simultaneous credible bands for latent Gaussian processes, and more generally for mixtures of Gaussian processes. This was investigated by Bolin, Guttorp, Januzzi, Jones, Novak, Podschwit, Richardson, Särkkä, Sowder, and Zimmerman (2015) and excursions contains a slightly more general implementation of the methods from these papers.
The package supports at least three ways of specifying the model that shuold be analysed. The standard method for purely Gaussian models is to specify the model by providing the parameters of the Gaussian process. For more general models, the input can either be given as Monte Carlo simulations of the process or as the result from an analysis using the R-INLA software package (Lindgren and Rue 2015, package available from http://r-inla. org/download/). The package is available to install from the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=excursions. A development version of the package is available via the repository https://bitbucket.org/davidbolin/ excursions. The development version is updated more frequently, and can be easily be installed directly in R as described on the repository homepage.
The following sections describe the theoretical methods used in the package and provide an introduction to the implementation. Section 2 summarizes the theory, Section 3 introduces the main functions in the package, and Section 4 contains two examples that illustrate how the package can be used. Finally, future plans for the package is discussed in Section 5.

Definitions and computational methodology
Hierarchical models are of great importance in many areas of statistics. In the simplest form, a hierarchical model has a likelihood distribution π(Y|X, θ) for observed data Y, which is specified conditionally on a latent process of interest, X, which has a distribution π(X|θ). For Bayesian hierarchical models, one also specifies prior distributions for the model parameters θ. The most important special case of these models are the LGMs, which are obtained by assuming that X|θ has a Gaussian distribution. Numerous applications can be studied using models of this form, and these are therefore the main focus of the methods in excursions.
A statistical analysis using an LGM often concludes with reporting the posterior mean E(X|Y) as a point estimate of the latent field, possibly together with posterior variances as a measure of uncertainty. In many applications, however, reporting posterior means and variances are not enough. As stated in the introduction, one may be interested in computing regions where the latent field exceeds some given threshold, contour curves with their associated uncertainty, or simultaneous confidence bands. In some applications, only a contour map of the posterior mean is reported, where the number of levels in the contour map should represent the uncertainty in the estimate. These are quantities that can be computed with excursions and we now define these in more detail before outlining how they can be computed. For details we refer to Lindgren (2015, 2016b).

Definitions
The main quantities that can be computed using excursions are (1) excursion sets, (2) contour credible regions and level avoiding sets, (3) excursion functions, (4) contour maps and their quality measures, (5) simultaneous confidence bands. This section defines these in more detail.
Throughout the section, X(s) will denote a stochastic process defined on some domain of interest, Ω, which we assume is open with a well-defined area |Ω| < ∞. Since it is not necessary for the definitions, we will not explicitly state the dependency on the data, but simply allow X to have some possible non-stationary distribution. In practice, however, the distribution of X will typically be a posterior distribution conditionally on data, X(s)|Y. For frequentist models, the distribution of X could also be conditionally on for example a maximum likelihood estimate of the parameters, X(s)|Y, θ.

Excursion sets
An excursion set is a set where the process X(s) exceeds or goes below some given level of interest, u. A where X(s) > u is referred to as a positive excursion set, whereas a set where X(s) < u is referred to as a negative excursion set. If X(s) = f (s) is a known function, these sets can be computed directly as A + u (f ) = {s ∈ Ω; f (s) > u} and A − u (f ) = {s ∈ Ω; f (s) < u} respectively. If X(s) is a latent random process, one can only provide a region where it with some (high) probability exceeds the level. More specifically, the positive level u excursion set with probability α, E + u,α (X), is defined as the largest set so that with probability 1 − α the level u is exceeded at all locations in the set, Similarly, the negative u excursion set with probability α, E − u,α (X), is defined as the largest set so that with probability 1 − α the process is below the level u at all locations in the set. This set is obtained by replacing A + u (X) with A − u (X) in (1).

Contour credible regions and level avoiding sets
For a function f , a contour curve (or set in general) of a level u is defined as the set of all level u crossings. Formally, the level curve is defined as B o denotes the interior of the set B and B c denotes the complement. Note that A c u (f ) not only includes the set of locations where f (s) = u, but also all discontinuous level crossings.
For a latent random process X, one can only provide a credible region for the contour curve. A level u contour credibility region, E c u,α (X), is defined as the smallest set such that with probability 1 − α all level u crossings of X are in the set. This set can be seen as the complement of the level u contour avoiding set E u,α (X), which is defined as the largest union where the sets (D + , D − ) are open. The sets M + u,α and M − u,α are denoted as the pair of level u avoiding sets. The notion of level avoiding sets can naturally be extended to multiple levels u 1 < u 2 < · · · < u k , which is needed when studying contour maps. In this case, the multilevel contour avoiding set is denoted C u,α (X) (For a formal definition, see Bolin and Lindgren 2016b). Bolin and Lindgren (2015) introduced excursion functions as a tool for visualizing excursion sets simultaneously for all values of α. For a level u, the positive and negative excursion functions are defined as F + u (s) = 1 − inf{α; s ∈ E + u,α } and F − u (s) = 1 − inf{α; s ∈ E − u,α }, respectively. Similarly, the contour avoidance function, and the contour function are defined as F u (s) = 1 − inf{α; s ∈ E u,α } and F c u (s) = sup{α; s ∈ E c u,α }, respectively. Finally, for levels u 1 < u 2 < · · · < u k , one can define a contour map function as F (s) = sup{1 − α; s ∈ C u,α }.

Excursion functions
These functions take values between zero and one and each set E • u,α can be retrieved as the 1 − α excursion set of the function F • u (s). An example of an excursion set and the corresponding excursion function is shown in Figure 2.

Contour maps and their quality measures
For a function f (s), a contour map C f with contour levels u 1 < u 2 < . . . < u K is defined as the collection of contour curves A c u 1 (f ), . . . , A c u K (f ) and associated level sets G k = {s : u k < f (s) < u k+1 }, for 0 ≤ k ≤ K, where one defines u 0 = −∞ and u K+1 = ∞. In practice, a process X(s) is often visualized using a contour map of the posterior mean E(X(s)|Y). The contour maps is visualized either by just drawing the contour curves labeled by their values, or by also visualising each level set in a specific color. The color for a set G k is typically chosen as the color that corresponds to the level u e k = (u k + u k+1 )/2 in a given color map. An example of this is shown in Figure 2.
In order to choose an appropriate number of contours, one must be able to quantify the uncertainty of contour maps. The uncertainty can be represented using a contour map quality measure P , which is a function that takes values in [0, 1]. Here, P should be chosen in such a way that P ≈ 1 indicates that the contour map, in some sense, is appropriate as a description of the distribution of the random field, whereas P ≈ 0 should indicate that the contour map is inappropriate.
An example of a contour map quality measure is the normalized integral of the contour map function The most useful quality measure is denoted P 2 and is defined as the simultaneous probability for the level crossings of (u e 1 , . . . , u e K ) all falling within their respective level sets (G 1 , . . . , G K ) (for details, see Bolin and Lindgren 2016b).
An intuitively interpretable approach for choosing the number of contours in a contour map is to find the largest K such that P 2 is above some threshold. For a joint credibility of 90%, say, choose the largest number of contours such that P 2 ≥ 0.9. How this can be done using excursions is illustrated in Section 4.2.

Simultaneous confidence bands
Especially for time series applications, the uncertainty in the latent process is often visualized using pointwise confidence bands. A pointwise confidence interval for X at some location s is given by [q α/2 (s), q 1−α/2 (s)], where q α (s) denotes the α-quantile in the marginal distribution of X(s).
A problem with using pointwise confidence bands is that there is not joint interpretation, and one is therefore often of interested in computing simultaneous confidence bands. For a process X(s), s ∈ Ω, we define a simultaneous confidence band as the region {(s, y) : s ∈ Ω, q ρ (s) ≤ y ≤ q 1−ρ (s)}. Here ρ is chosen such that P(q ρ (s) < X(s) < q 1−ρ (s), s ∈ Ω) = 1 − α. Thus α controls the probability that the process is inside the confidence band at all locations in Ω. An example of pointwise and simultaneous confidence bands is given in Figure 3.

Computational methods
If the latent process X(s) is defined on a continuous domain Ω, one has to use a discretization, x, of the process in the statistical analysis. The vector x may be the process evaluated at some locations of interest or more generally the weights in a basis expansion X(s) = i ϕ i (s)x i . Computations in the package are in a first stage performed using the distribution of x. If Ω is continuous, computations in a second stage interpolates the results for x to Ω. In this section, we briefly outline the computational methods used in these two stages. As most quantities of interest can be obtained using some excursion function, we focus on how these are computed. As a representative example, we outline how F + u (s) is computed in the following sections. As before, let Y and θ be a vectors respectively containing observations and model parameters.

Computing an excursion function F
This means that F u can be obtained by first reordering the nodes and then computing a sequential integral. The reordering is in this case obtained by sorting the marginal probabilities P(x i > u) (for other options, see Bolin and Lindgren 2015). After reordering, the i: Using sequential importance sampling as described below, the whole sequence of integrals can be obtained with the same cost as computing only one integral with i = n, making the computation of F + u feasible also for large problems.

Gaussian integrals
The basis for the computational methods in the package is the ability to compute the required integral when the posterior distribution is Gaussian. In this case, one should compute an integral Here µ and Q are the posterior mean and posterior precision matrix respectively. To take advantage of the possible sparsity of Q if a Markovian model is used, the integral is rewritten as where, if the problem has a Markov structure, x i |x i+1:d only depends on a few of the elements in x i+1:d given by the Markov structure. The integral is calculated using a sequential importance sampler by starting with the integral of the last component π(x d ) and then moving backward in the indices, see Bolin and Lindgren (2015) for further details.

Handling non-Gaussian data
Using the sequential importance sampler above, F + u can be computed for Gaussian models with known parameters. For more complicated models, the latent Gaussian structure has to be handled, and this can be done in different ways depending on the accuracy that is needed. excursions currently supports the following five methods described in detail in Bolin and Lindgren (2015): Empirical Bayes (EB), Quantile Correction (QC), Numerical Integration (NI), Numerical integration with Quantile Corrections (NIQC), and improved Numerical integration with Quantile Corrections (iNIQC).
The EB method is the simplest method and is based on using a Gaussian approximation of the posterior, π(x|Y) ≈ π G (x|Y, θ). The QC method uses the same Gaussian approximation but modifies the limits in the integral to improve the accuracy. The three other methods are intended for Bayesian models, where the posterior is obtained by integrating over the parameters, The NI method approximates the integration with respect to the parameters as in the INLA method, using a sum of representative parameter configurations, and the NIQC and iNIQC methods combines this with the QC method to improve the accuracy further.
In general, EB and QC are suitable for frequentist models and for Bayesian models where the posterior distribution conditionally on the parameters is approximately Gaussian. The methods are equivalent if the posterior is Gaussian and in other cases QC is more accurate than EB. For Bayesian models, the NI method is, in general, more accurate than the QC method, and for non-Gaussian likelihoods, the NIQC and iNIQC methods can be used for improved results. In general the accuracy and computational cost of the methods are as follows If the main purpose of the analysis is to construct excursion or contour sets for low values of α, we recommend using the QC method for problems with Gaussian likelihoods and the NIQC method for problems with non-Gaussian likelihoods. The increase in accuracy of the iNIQC method is often small compared to the added computational cost.

Continuous domain interpolations
For a continuous spatial domain, the excursion function F + u (s) can be approximated using F + u computed at discrete point locations. The main idea is to interpolate F + u assuming monotonicity of the random field between the discrete computation locations. Specifically, assume that the values of F + u correspond to the values at the vertices of some triangulated mesh such as the one shown in the left panel of Figure 2. If the excursion set E + u,α (X) should be computed for some specific value of α, one has to find the 1 − α contour for F + u (s). For interpolation to a location s within a specific triangle T with corners in s 1 , s 2 , and s 3 , excursions by default uses log-linear interpolation, are the barycentric coordinates of s within the triangle.
Further technical details of the continuous domain construction are given in Bolin and Lindgren (2016b). Studies of the resulting continuous domain excursion sets in Bolin and Lindgren (2016b) indicate that the log-linear interpolation method results in sets with coverage probability on target or slightly above target for large target probabilities. An example of a continuous domain excursion set for a triangulated mesh is shown in the middle panel of

Implementation
The functions in excursions can be divided into four main categories depending on what they compute: (1) Excursion sets and credible regions for contour curves, (2) Quality measures for contour maps, (3) Simultaneous confidence bands, and (4) Utility such as Gaussian integrals and continuous domain mappings. The main functions come in at least three different versions taking different input: (1) The parameters of a Gaussian process, (2) results from an analysis using the R-INLA software package, and (3) Monte Carlo simulations of the process. These different categories are described in further detail below.
Much of the computations in the package is done in C functions. These functions use methods from a number of C and Fortran libraries, such as BLAS (Dongarra, Du Croz, Hammarling, and Duff 1990), LAPACK (Anderson, Bai, Bischof, Blackford, Demmel, Dongarra, Du Croz, Hammarling, Greenbaum, McKenney, , and Sorensen 1999), and CHOLMOD (Chen, Davis, Hager, and Rajamanickam 2008) for efficient matrix manipulations together with function from the GNU Scientific library (Galassi and Gough 2006) and several different reordering methods. Notably, the CAMD library Duff 1996, 2004) is used for constrained approximate minimum degree orderings.

Excursion sets and contour credible regions
The main function for computing excursion sets and contour credible regions is excursions.
A typical call to the function looks like R> res.exc <-excursions(mu = mu.post, Q = Q.post, alpha = 0.1, type = '>', + u = 0, F.limit = 1) Here, mu and Q are the mean vector and precision matrix for the joint distribution of the Gaussian vector x. The arguments alpha, u, and type are used to specify what type of excursion set that should be computed: alpha is the error probability, u is the excursion or contour level, and type determines what type of region that is considered: '>' for positive excursion regions, '<' for negative excursion regions, ' !=' for contour avoiding regions, and '=' for contour credibility regions. Thus, the call above computes the excursion set E + 0,0.1 . The argument F.limit is used to specify when to stop the computation of the excursion function. In this case with F.limit=1, all values of F + u are computed, but the computation time can be reduced by decreasing the value of F.limit.
The function has the EB method as default strategy for handling the possible latent Gaussian structure. In the simulated example, the likelihood is Gaussian and the parameters are assumed to be known, so the EB method is exact. The QC method can be used instead by specifying method='QC'. In this case, the argument rho should be used to also provide a vector with point-wise marginal probabilities: P (x i > u) for positive excursions and contour regions, and P (x i < u) for negative excursions. In the situation when π(x|Y, θ) is Gaussian but π(x|Y) is not, the marginal probabilities should be calculated under the distribution π(x|Y) and mu and Q should be chosen as the mean and precision for the distribution π(x|Y,θ) wherê θ is the MAP or ML estimate of the parameters.

The INLA interface
The function excursions.inla is used to compute excursion sets and credible regions for contour curves for latent Gaussian models that have been estimated using R-INLA. It takes the same arguments excursions, except that mu and Q are replaced with arguments related to R-INLA. A basic call to the function excursions.inla looks like R> excursions.inla(result.inla, name, alpha, u, type) Here result.inla is the output from a call to the inla function, and name is the name of the component in the output that the computations should be done for. For more complicated models, one typically specifies the model in R-INLA using an inla.stack object. If this is done, the call to excursions.inla will instead look like R> excursions.inla(result.inla, stack, tag, alpha, u, type) Here stack is now the stack object and tag is the name of the component in the stack for which the computations should be done. The typical usage of the tag argument is to have one part of the stack that contains the locations where measurements are taken, and another that contains the locations where the output should be computed.
excursions.inla has support for all strategies described in Section 2 for handling latent Gaussian structures: The argument method can be one of 'EB', 'QC', 'NI', 'NIQC', and 'iNIQC'.

Analysis of Monte Carlo samples
The function excursions.mc can be used to post-process Monte Carlo model simulations in order to compute excursion sets and credible regions. For this function, the model is not specified explicitly. Instead a d × N matrix X containing N Monte Carlo simulations of the d dimensional process of interest is provided. A basic call to the function looks like R> excursions.mc(X, u, type) where u again determines the level of interest and type defines the type of set that should be computed. It is important to note that this function does all computations purely based on the Monte Carlo samples that are provided, and it does not use any of the computational methods based on sequential importance sampling for Gaussian integrals that the is the basis for the previous methods. This means that this function in one sense is more general as X can be samples from any model, not necessarily a latent Gaussian model. The price that has to be payed for this generality is that the only way of increasing the accuracy of the results is to increase the number of Monte Carlo samples that are provided to the function.

Analysis of contour maps
The main function for analysis of contour maps is contourmap. A basic call to the function looks like R> res.con <-contourmap(mu = mu.post, Q = Q.post, n.levels = 4, + alpha = 0.1, compute = list(F = TRUE, measures = c("P0"))) Here, mu is again the mean value and Q is the precision matrix of the model. The parameter n.levels sets the number of contours that should be used in the contour map, and these are spaced equidistant in the range of mu by default. Other types of contour maps can be obtained using the type argument. For manual specification of the levels, the levels argument can be used instead. By default, the function computes the specified contour map but no quality measures and it does not compute the contour map function. If quality measures should be computed, this is specified using the compute argument. This argument is also used to decide whether the contour map function F should be computed.
As for excursions, this function comes in two other versions depending on the form of the input: contourmap.inla for model specification using an R-INLA object, or contourmap.mc for model specification using Monte Carlo simulations of the model. The model specification using these functions is identical to that in the corresponding excursions functions.

Continuous domain interpretations
A common scenario is that the input used in contourmap or excursions represents the value of the model at some discrete locations in a continuous domain of interest. In this case, the function continuous can be used to interpolate the discretely computed values by assuming monotonicity of the random field in between the discrete computation locations, as discussed in Section 2. A typical calls to the function looks like R> sets.exc <-continuous(ex = res.exc, geometry = mesh, alpha = 0.1) Here ex is the result of the call to contourmap or excursions and alpha is the error probability of interest for the excursion set or credible region computation. The argument geometry specifies the geometric configuration of the values in input ex, either as a general triangulation geometry or as a lattice. A lattice can be specified as an object of the form list(x, y) where x and y are vectors, or as list(loc, dims) where loc is a two-column matrix of coordinates, and dims is the lattice size vector. If R-INLA is used, the lattice can also be specified as an inla.mesh.lattice object. In all cases, the input is treated topologically as a lattice with lattice boxes that are assumed convex. A triangulation geometry is specified as an inla.mesh object. Finally, an argument output can be used to specify what type of object should be generated. The options are currently sp which gives a SpatialPolygons (Bivand, Pebesma, and Gómez-Rubio 2013) object, and inla which gives an inla.mesh.segment object.

Simultaneous confidence bands
The function simconf can be used for calculating simultaneous confidence bands for a Gaussian process X(s). A basic call to the function looks like R> simconf (alpha, mu, Q) where alpha is the coverage probability, mu is the mean value vector for the process, and Q is the precision matrix for the process. The function has a few optional arguments similar to those of excursions, all listed in the documentation of the function. The function returns upper and lower limits for both pointwise and simultaneous confidence bands.
As for excursions and contourmap, there is also a version of simconf that can be used to analyze R-INLA models (simconf.inla) and a version that can analyze Monte Carlo samples (simconf.mc). Furthermore, there is a version simconf.mixture which is used to compute simultaneous confidence regions for Gaussian mixture models with a joint distribution on the form This particular function was used to analyze the models in  and Guttorp, Januzzi, Novak, Podschwit, Richardson, Sowder, Zimmerman, Bolin, and Särkkä (2014), but is also used internally by simconf.inla.

Gaussian integrals
Among the utility functions in the package, the function gaussint can be especially useful also in a larger context. It contains the implementation of the sequential importance sampling method for computing Gaussian integrals, described in Section 2. This function has two features that separates it from many other functions for computing Gaussian integrals: Firstly it is based on the precision matrix of the Gaussian distribution, and sparsity of this matrix can be utilized to decrease computation time. Secondly, the integration can be stopped as soon as the value of the integral in the sequential integration goes below some given value 1 − α. If one only is interested in the exact value of the integral given that it is larger than some value 1 − α, this option can save a lot of computation time.
A basic call to the function looks like

R> gaussint(mu, Q, a, b)
where mu is the mean value vector, Q is the precision matrix, a is a vector of the lower limits in the integral, and b contains the upper integration limits. If the Cholesky factor of Q is known beforehand, this can be supplied to the function using the Q.chol argument. An argument alpha is used to set the computational 1 − α limit for the integral. The function returns an estimate of the integral as well as an error estimate. If the error estimate is too high, the precision can be increased by increasing the n.iter argument of the function.

Plotting
The excursions package contains various functions that are useful for visualization. The function tricontourmap can be used for visualization of contour maps computed on triangulated meshes. The following code plots the posterior mean using the contour map we previously computed.

Two applications
In this section, two examples are used to illustrate how excursions can be used.

Time series data: Tokyo rainfall
To illustrate the methods in the package, we use the much analyzed binomial time series from Kitagawa (1987). Each day during the years 1983 and 1984, it was recorded whether there was more than 1 mm rainfall in Tokyo. Of interest is to study the underlying probability of rainfall as a function of day of the year. The data is modelled as y i ∼ Bin(n i , p i ) for calendar day i = 1, ..., 366. Here n i = 2 for all days except for February 29 (i = 60) which only occurred during the leap year of 1984. The probability p i is modeled as a logit-transformed Gaussian process.
R> result <-inla(formula, family = "binomial", data = data, + Ntrials = data$Ntrials, control.predictor = list( + A = inla.stack.A(stack), link = data$link, compute = TRUE), + control.compute = list(config = TRUE)) We now have estimates of the posterior mean and marginal confidence intervals for the probability of rain for each day. However, if we also want joint confidence bands, we can estimate these using simconf.inla as R> res <-simconf.inla(result, stack, tag = "est", alpha = 0.05, link = TRUE) Note the argument link=TRUE which tells the function that the results should be returned in the scale of the data, and not in the scale of the linear predictor. Next, we plot the results, showing the marginal confidence bands with dashed lines and the simultaneous confidence band with dotted lines. The results are shown in Figure 3.

Spatial data: Parana precipitation
We will now use a spatial dataset to illustrate how excursions can be used. In order to keep the parts of the code that are not relevant to excursions brief, we again use data available q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q in R-INLA. The dataset consists of daily rainfall data for each day of the year 2011, at 616 locations in and around the stated of Paraná in Brazil. In the following analysis, we analyse the data from the month of January.
The statistical model used for the data is a latent Gaussian model, where the precipitation measurements are assumed to be Γ-distributed with a spatially varying mean. The mean is modeled as a log-Gaussian Matérn field specified as an SPDE model. Details of the current model model and the following INLA implementation can be found in the SPDE tutorial available on the R-INLA homepage, see Wallin and Bolin (2015) for an analysis of the data using a different non-Gaussian SPDE model.
We start by loading the data and defining the model: R> coords <-as.matrix(PRprec[ind, 1 : 2]) R> b <-inla.nonconvex.hull(coords, -0.03, -0.05, resolution = c(100, 100)) R> prmesh <-inla.mesh.2d(boundary = b, max.edge = c(.45, 1), cutoff = 0.2) R> A <-inla.spde.make.A(prmesh, loc = coords) R> spde <-inla.spde2.matern(prmesh, alpha = 2) R> mesh.index <-inla.spde.make.index(name = "field", n.spde = spde$n.spde) The measurement stations are spatially irregular, but we are interested in making predictions to a regular lattice within the state. In order to do continuous domain interpretations, we define the lattice locations as a lattice object using the submesh function of the excursions package. The function inout from the package splancs (Rowlingson and Diggle 2015) is used to find the locations on the lattice that are within the region of interest. We now define the stack objects and estimate the model using inla. Again note that we have to set the control.compute argument of the results is to be used by excursions.
R> con <-tricontourmap(submesh, z=exc$mean, levels = log (7)) We plot the resulting continuous domain excursion function together with the contour curve using the following commands. The result is shown in the right panel of Figure 4.
To visualize the posterior mean using a contour map using the following commands R> lp <-contourmap.inla(r, stack = stk, tag = "prd", n.levels = 2, + compute = list(F = FALSE)) R> tmap <-tricontourmap(submesh, z = exc$mean, levels = lp$meta$levels) R> plot(tmap$map, col = contourmap.colors(lp, col = cmap)) Here, the contourmap.inla computes the levels of the contour map and tricontourmap computes the contour map on the mesh. Finally, contourmap.colors is used to compute appropriate colours for visualising the contour map. The result is shown in the left panel of Figure 4.
The contour map we computed had two contours, and a relevant question is now if this is an appropriate number. To investigate this, we compute the P 2 quality measure for this contour map and for contour maps with one and three levels.

Discussion
The excursions package was first developed for calculating probabilistic excursion sets and contour credible regions for latent Gaussian stochastic processes and fields. Since the early versions, the scope of the package has grown and it now contains functions for several other related computations, such as simultaneous confidence bands, uncertainty quantification of contour maps, and computations of Gaussian integrals.
Some of the functionaliy in excursions can also be found in other packages. For example, Ex-ceedanceTools (French 2014) provides an alternative method for computing excursion sets and uncertainty regions for contour curves for Gaussian models. The package mvtnorm (Genz, Bretz, Miwa, Mi, Leisch, Scheipl, and Hothorn 2016) contains functions for computing integrals of Gaussian distributions, but does not have support for computations based on sparse precision matrices. Finally, there are numerous packages with the ability to compute simultaneous confidence bands for various models, such as semiparametric regression models using AdaptFitOS (Wiesenfarth, Krivobokova, Klasen, and Sperlich 2012) or nonparametric regression models for functional data using SCBmeanfd (Degras 2016). Comparing the methods in excursions to the corresponding methods in these packages is an intersting topic for future studies.
Future work also includes further development of the package, and new features are added as they are needed. One main focus for the current development is to add functionality for large scale computations, where sparse Cholesky factorisation is computationally infeasible. We also plan to add functionality for computations needed in spatial extreme value theory, where certain Gaussian integrals often are needed during likelihood inference. Finally, the technical aspects of the functions are also being improved, and future releases will introduce improvements to both the continuous domain interpretations and to the algorithms for finding the largest excursion sets.