Excursion and contour uncertainty regions for latent Gaussian models

Introduction

In many statistical applications, one is interested in finding areas where the studied process exceeds a certain level or is significantly different from some reference level. A typical example is in studies of air pollution, where one is interested in testing if, and where, the pollution level exceeds some given limit value set by some regulatory agency (Cameletti et al., 2012), and similar examples can be found in a wide range of scientific fields including brain imaging (Marchini and Presanis, 2003) and astrophysics (Beaky et al., 1992). In spatio-temporal applications one might be interested in finding regions that have experienced significant changes over the studied time period. This is a common problem in climate science and the studied quantity can for example be temperature (Furrer et al., 2007), precipitation (Sain et al., 2011), or vegetation (Eklundh and Olsson, 2003; Bolin et al., 2009).

The method derived here is based on using a parametric family for the excursion sets in combination with a sequential importance sampling method for estimating joint probabilities. For a specific choice of the parametric family, the method is equivalent to the thresholding methods mentioned above, with the important difference that the correct joint distribution is used when selecting the threshold. The method is extended using more general parametric families, and the related problem of finding uncertainty regions for contour curves is treated using the same methodology.

The structure of the article is as follows. In Section 2, the problem is formulated and definitions for excursion sets and uncertainty regions for contour curves are given. In Section 3, a method for estimating these sets is proposed. Estimating the sets is the most difficult problem as one easily runs into computational difficulties arising from having to evaluate high-dimensional integrals. In Section 4, the methods are tested on a few simulated examples to test the method’s accuracy. Two applications to real data are covered in Section 5, the first considers air pollution data from the North-Italian region Piemonte, and the second considers estimation of spatially dependent vegetation trends in the African Sahel. Finally, a few remarks and comments are given in Section 6.

Problem formulation

There are a number of different ways one could formulate excursion sets, and not all of them are useful from a practical point of view. Hence, in this section we will formalise the problem and discuss how the results should be interpreted. More precisely we look at two connected problems. The first one is to find areas where a stochastic process exceeds a given level with some probability and the second one is to quantify the uncertainty in contour curves of stochastic fields.

The set of all such points is the complement of the union of the interior sets of the positive and negative excursion sets.

where $A^{o}$ is the interior, relative to $\Omega$, of the set $A$ and $A^{c}$ is the complement.

Taking the interiors of the sets $A_{u}^{+}(f)$ and $A_{u}^{-}(f)$ is important. Consider for example the following function on $\Omega=$

is the negative level $u$ excursion set with probability $1-\alpha$.

It is important to realize how the excursion set $E_{{u,\alpha}}^{{+}}(x)$ should be interpreted: It is the largest set so that the level $u$ is exceeded at all locations in the set with probability $1-\alpha$, and therefore is a smaller set than $D_{m}$ defined in (1), which is the set of points where the marginal probability for exceeding the level is at least $1-\alpha$. Another possible definition of an excursion set would be a set that contains all excursions with probability $1-\alpha$. This is a larger set than $D_{m}$, given by $E_{{u,\alpha}}^{{-}}(x)^{c}$. Which set one is interested in depends on the application, but it can be a good idea to calculate both to get a better understanding of the uncertainties in the problem.

In certain applications, one might be interested in joint positive and negative excursions from some level, for example when doing simultaneous regressions and one is interested in finding regions where the slopes are significantly different from zero (see Section 5.2 for a possible scenario of this kind).

Denote the union of these two sets the level avoiding set $M_{u,\alpha}$:

The probability calculation in Definition 2.4 can now be reformulated as an ordinary excursion probability in $y$:

Similarly to how the contour sets for deterministic functions were defined, the pair of level avoiding sets can now be used to define uncertainty regions for contour curves.

is then an uncertainty region for the contour set of level $u$.

The interpretation of this uncertainty region is important. The set $M^{c}_{u,\alpha}$ is the smallest set such that all level $u$ crossings of $x$ are in the set with probability $1-\alpha$. One should note that this definition of the uncertainty region for level curves is different from some other definitions in the literature. For example, Lindgren and Rychlik (1995) define uncertainty regions as a union of intervals where each interval contains a single level crossing with probability $1-\alpha$.

It is somewhat unsatisfactory that the sets defined here are made unique by finding the largest set satisfying a certain restriction. The set $E_{{u,\alpha}}^{{+}}(x)$ is for example defined as the largest set $D$ satisfying $\mathsf{P}(D\subseteq A_{u}^{+}(x))\geq 1-\alpha$, but there are also many other smaller sets satisfying the requirement, and these are not seen if only $E_{{u,\alpha}}^{{+}}(x)$ is reported. Also, if one wants to know where the field likely exceeds the level $u$, the set $E_{{u,\alpha}}^{{+}}(x)$ might not be sufficient since it does not provide any information about the locations not contained in the set. Therefore, it would be good to have something similar to $p$-values, i.e. the marginal probabilities of exceeding the level, but which can be interpreted simultaneously. To that end we introduce the excursion function, level avoidance function, and contour function as visual tools for answering such questions.

The positive and negative $u$ excursion functions are given by

Similarly, the level avoidance and contour functions are given by

Computations

There are now, in principle, two main problems that have to be solved in order to find the excursion sets, level avoidance sets, or contour uncertainty sets:

Use shape optimization to find largest region $D$ satisfying the required probability constraint.

Hence, given a method to solve each of the two problems, one could simply run the shape optimization algorithm and in each iteration calculate the required probability using the integration method. In theory there are no problems doing this, but in practice the integration method will be computationally demanding and it may not be feasible to use this strategy for applications involving large data sets. Therefore, we instead propose a slightly different strategy that will minimize the number of calls to the integration method by solving the problem sequentially. We first outline the strategy in the simplest possible situation, which will be used as a basis for all other more complicated strategies.

The method is based on using an increasing parametric family for the excursion sets in combination with a sequential integration routine for calculating the probabilities. The advantage with using a sequential integration routine is that if the required probability has been calculated for some set $D_{1}$, then the calculation for a larger set $D_{2}\supset D_{1}$ can be based on the result for $D_{1}$, resulting in large computational savings.

Choose a suitable (sequential) integration method for the problem.

Reorder the nodes to the order they will be added to the excursion set when the parameter $\rho$ is increased.

Before extending this method to more general situations, we go into more detail on how to do the steps in Algorithm 3.1 in practice. In Section 3.1, a few sequential integration methods are presented. In Section 3.2, some different parametric families for the excursion sets and level avoidance sets are introduced and Algorithm 3.1 is extended using two-parameter families. The problem of how to optimally reorder the nodes is also discussed in this section. Finally in Section 3.3, three different methods are proposed for calculating excursion sets under the full posterior distribution (2).

The simplest way of approximating (3) is to use Monte-Carlo (MC) integration. However, estimating the probability with any reasonable accuracy using standard MC integration is often too computationally expensive. Fortunately there are a number of variance reduction techniques that can be used to increase the efficiency.

A key step in many numerical integration techniques is to transform the integral to make it more suitable for integration. Notably, Genz (1992) derived such a transformation for the Gaussian integral (3), though similar transformations have been suggested by other authors as well (see e.g. Geweke, 1991). Besides transforming the integral to the unit hyper cube, the transformation also achieves a separation of the variables so that the full problem can be calculated sequentially. The integral can then efficiently be evaluated using a quasi MC (QMC) method where the uniform random numbers in the ordinary MC integrator are replaced by some deterministic sequence of points chosen to reduce the probabilistic error bound of the crude MC integrator, see Genz and Bretz (2009) for details.

A final variance reduction technique for the general integration problem is to reorder the variables before calculating the integral, as first suggested by Schervish (1984) and later improved by Gibson et al. (1994). These reorderings can reduce the error by an order of magnitude, as shown by Genz and Bretz (2002). However, the technique will not be applicable in our situation since the reordering will be determined by a parametric family for the excursion sets.

A common assumption in spatial statics and image analysis is that the latent field can be modeled, or approximated, using a Gaussian Markov random field (GMRF). See Rue and Held (2005) for an introduction to GMRFs, and note that GMRFs are used also for modeling in continuous space, for example using the SPDE approach by Lindgren et al. (2011). One of the motivating reasons for using GMRFs is that it reduces the computational cost for parameter estimation and spatial prediction, and because of this one would also like to be able to use the Markov property in the calculation of (3).

and note that the integral is the normalizing constant to the truncated density

Proceed like this, simulating from the truncated conditional distributions and in each step updating the importance weights recursively through

2 Parametric families

which are easy to calculate using only the marginal posterior distributions. The simplest one-parameter family based on the marginal quantiles is given in the following definition.

Using this one-parameter family in Algorithm 3.1 is equivalent to finding a threshold value for the marginal excursion probabilities to get the correct simultaneous significance level. It is thus similar to the thresholding algorithms discussed in Marchini and Presanis (2003) but with the important difference that the correct joint, often non-stationary, posterior density is used when finding the threshold.

The simple one-parameter family can be extended in a number of ways, for example by considering other levels in the excursion sets.

The sets $D_{1}^{+}(v,\rho)$ and $D_{1}^{-}(v,\rho)$ are increasing in $\rho$ for a fixed $v$.

Let $p_{i}^{\tau}$ be the smoothed marginal positive $u$ excursion probabilities, using a circular averaging filter with radius $\tau$. A two-parameter family for the positive and negative $u$ excursion sets is then given by

Using the two-parameter families requires a modification to Algorithm 3.5, resulting in a slightly more computationally demanding method.

Choose a suitable (sequential) integration method for the problem.

Select a suitable one-dimensional optimization strategy.

Do optimization of the size of $D(\nu,\bullet)$ over $\nu$:

For the current value of $\nu$, reorder the nodes to the order they will be added to the excursion set when the parameter $\rho$ is increased.

$E_{{u,\alpha}}^{{+}}$ is given by the largest set $D$ found in the optimization over $\nu$.

The optimization can in this case be done using a Golden section search or a similar fast optimization procedure for one-dimensional problems. Algorithm 3.5 can also be used to estimate uncertainty regions for contour curves by using the following two-parameter family for the pair of level avoiding sets.

Let $D_{1}^{+}(\rho_{1})$ and $D_{1}^{-}(\rho_{2})$ be given by Definition 3.2. A two-parameter family for the pair of level avoiding sets is obtained as $(D_{1}^{+}(\rho_{1}),D_{1}^{-}(\rho_{2}))$. A one-parameter family is obtained by requiring that $\rho_{1}=\rho_{2}=\rho$.

The one-parameter family in Definition 3.6 can be used in Algorithm 3.1 to estimate level avoiding sets and uncertainty regions for contour curves without having to use the more computationally expensive Algorithm 3.5.

In the case of a GMRF posterior, it is desirable to make the Cholesky factor of the precision matrix as sparse as possible, because it reduces the number of floating point calculations that have to be done and reduces the error of the estimator. Reordering the nodes according to a parametric family does not guarantee good sparsity of the Cholesky factor, but the reordering can be improved by finding upper and lower bounds for the region.

The simplest upper bound for the region is to use

A simple lower bound for the region is obtained using Boole’s inequality as

where $n$ is the number of points in the discretization of the domain. In terms of multiple hypothesis testing, this lower bound is obtained from the classical Bonferroni correction method and an improved lower bound can be obtained using the Holm-Bonferroni method (Holm, 1979) as

The nodes can now be categorized into three classes, the first class contains the nodes included in the lower bound $L_{2}$, the second class contains the nodes in the set $U_{1}\setminus L_{2}$ and the third class contains all other nodes. Since one knows that all nodes in $L_{2}$ will be included in $D$, these can be reordered to maximize the sparsity of the Cholesky factor, for example using an approximate minimum degree permutation. The nodes in the second class are then added in the order determined by the parametric family. Finally, since the nodes in the third class will not be included in the domain, these can be reordered to maximize the sparsity or integrated out of the posterior distribution. Making the bounds more precise will improve the sparsity of the problem and therefore reduce the Monte-Carlo error and the computational complexity.

3 Probability calculations for the latent Gaussian setting

(Numerical Integration) Numerically approximate the excursion probability by approximating the integral in (2) as

Tests on simulated data

In this section, three examples using simulated data are presented to illustrate the methods and test their accuracy. In the first example, we look at a problem in one dimension with known model parameters, where a latent Gaussian process with an exponential covariance is observed under Gaussian measurement noise. In the second example, we compare the different parametric families for contour uncertainty sets for a model in two dimensions with known parameters, where a latent Gaussian Matérn field is observed under Gaussian measurement noise. In the third example, the same spatial model setup is used, but this time the model parameters are estimated from data and the three methods for handling the full posterior distribution are compared.

We begin with a simple one-dimensional example to illustrate the different sets we have previously defined. Let $x(s)$, $s\in$ be a Gaussian process with an exponential covariance function with scaling parameter $\lambda=1$ and mean

By the definition of $F_{0}^{+}(s)$, the positive -excursion set $E_{{0,\alpha}}^{{+}}(x)$, is obtained by calculating the $1-\alpha$ excursion set of the function $F_{0}^{+}(s)$, and this set is shown for $\alpha=0.05$ in red in Figure 1, Panel (b). The grey set shows the upper bound $U_{1}$, which is the set where $\mathsf{P}(x(s)>0)\geq 1-\alpha$, and the dark red set shows the Holm-Bonferroni lower bound $L_{2}$. The black curve shows the kriging estimate of the process given the data. Note that the grey and red sets are obtained as excursion sets of the grey and red functions in Panel (a), and also note that $L_{2}\subset E_{{0,\alpha}}^{{+}}(x)\subset U_{1}$.

Finally in Figure 2, Panel (b), the -contour uncertainty region $M_{0,0.05}^{c}(x)$ is shown in red and the kriging estimate of $x(s)$ is again shown in black. The set was estimated using the two-parameter family for level avoidance sets from Definition 3.6 and Algorithm 3.5. The complement of this set is the union of the level avoiding sets $(M_{0,0.05}^{-}(x),M_{0,0.05}^{+}(x))$, which is the largest pair of sets $(D^{+},D^{-})$ satisfying $\mathsf{P}(D^{-}\subseteq A_{u}^{-}(x),\,D^{+}\subseteq A_{u}^{+}(x))\geq 0.95$.

2 Example 2: 2d Gaussian data with known parameters

To that end, we generate $50$ data sets using the same setup, and for each data set estimate $M_{0,0.05}^{c}(x)$, first using the one-parameter family $(D_{1}^{+}(\rho),D_{1}^{-}(\rho))$, and then using the more general two-parameter family $(D_{1}^{+}(\rho_{1}),D_{1}^{-}(\rho_{2}))$. Since the one-parameter family is a special case of the two-parameter family where ${\rho_{1}=\rho_{2}=\rho}$, the contour sets estimated with the two-parameter family should always be smaller than the one-parameter sets. However, using the two-parameter family, the estimated sets are on average only $0.2\%$ smaller than if the one-parameter family is used, so in this case it is arguably not worth the extra computational effort to use the two-parameter family, although for other levels $u$, or other latent models, the difference might be larger.

3 Example 3: 2d Gaussian data with unknown parameters

In this example we compare the three methods,described in Section 3.3, for handling the full posterior distribution (2) in the calculations. The same Gaussian Matérn model is used as in Example 2, with the difference that we now also estimate the parameters from the data.

There are three possible sources of errors in this comparison. The first one is the Monte-Carlo error from the estimation of $\hat{p}(\alpha)$, which has nothing to do with the accuracy of the method. The second error is the Monte-Carlo error in the probability estimation when estimating the excursion distribution functions. This error is, however many orders of magnitude smaller in this case. The final error is the approximation error induced by using any of the three methods EB, QC, or NI for handling the full posterior distribution.

Applications

In this section, we will use the techniques described above in two different applications. In the first, we study air pollution data from Piemonte region in northern Italy and estimate regions where the daily limit for PM10 (particulate matter with an aerodynamic diameter of less than 10 $\mu$m) is exceeded. In the second application, we study vegetation index data from the African Sahel and estimate areas that experienced a significant increase in vegetation after the drought period in the early 1980’s. In both of these applications the data sets are large, and the Markov structure of the latent Gaussian models has to be used in the calculations.

High levels of air pollution can be harmful for the ecosystems and the human health. The effects on human health ranges from minor effects to the cardio-respiratory system to premature mortality (Cohen et al., 2009; Cameletti et al., 2012). Because of this, environmental agencies have to assess the air quality in order to take proper actions for improving the situation in polluted areas, and an important tool in this process is the ability to produce continuous maps of air pollution.

A region where the daily limit values fixed by the European Union for human health protection (see EU Council Directive 1999/30/EC) are periodically exceeded is the Piemonte region in northern Italy. Recently, Cameletti et al. (2012) proposed a statistical model to capture the complex spatio-temporal dynamics of PM10 concentration in the region and used it to produce daily maps of PM10. They also produced daily maps of exceedance probabilities of the value $50\mu g/m^{3}$, which is the value fixed by the European directive 2008/50/EC for the daily mean concentration that cannot be exceeded more than 35 days in a year. These probability maps only considered the marginal excursion probabilities, and no attempts of producing maps of simultaneous exceedance probabilities were made. In the following, we will therefore consider the same model and data but also estimate the excursion functions for the $50\mu g/m^{3}$ limit value.

where the $p=9$ covariates $z_{k}$ are used and $\xi$ is a spatio-temporal Gaussian random field. Based on the work of Cameletti et al. (2011) the following covariates were used: 1) Daily mean wind speed; 2) daily maximum mixing height; 3) daily precipitation; 4) daily mean temperature; 5) daily emissions; 6) altitude; 7) longitude; 8) latitude; and 9) intercept. These covariates are provided with hourly temporal resolution on a $4$ km $\times$ $4$ km regular grid by the environmental agency of Piemonte region (Arpa Piemonte), see Finardi et al. (2008). The spatio-temporal process $\xi$ is assumed to follow first order autoregressive dynamics in time with spatially dependent innovations:

where $C(\cdot)$ is a Matérn covariance function given by (7). The model parameters and the posterior distribution for the latent field are then estimated using INLA in combination with the SPDE representation of Lindgren et al. (2011), see Cameletti et al. (2012) for details.

Note that the union of $E_{{50,0.1}}^{{+}}(x)$ and $E_{{50,0.1}}^{{-}}(x)$ covers only a small part of the region, indicating that the uncertainty in the problem is large. See the red and blue sets in the left panel of Figure 9. Also, by taking the complement of the set $E_{{50,0.1}}^{{-}}(x)$, we get the region that contains all exceedances of the level with certainty $0.9$, indicated in grey in the left panel of Figure 9. This set is large, indicating that there are many regions where the level possibly is exceeded. Hence, it is important to note that the positive excursion set $E_{{50,0.1}}^{{+}}(x)$ is small because the uncertainty is large in the problem, and not because the other regions certainly have concentrations below the level.

2 Spatially dependent temporal trends in vegetation data

Trends in vegetation cover are related to changes in climatic drivers, feedback mechanisms between the atmosphere and land surface, and human interaction. A region with rapid recent changes is the African Sahel. This zone has received much attention regarding desertification and climatic variations (Olsson, 1993; Nicholson, 2000; Lamb, 1982). Recently, Eklundh and Olsson (2003) observed a strong increase in seasonal vegetation index over parts of the Sahel using Advanced Very High Resolution Radiometer (AVHRR) data from the NOAA/NASA Pathfinder AVHRR Land (PAL) database (Agbu and James, 1994; James and Kalluri, 1994), for the period 1982-1999. The study was based on ordinary least squares linear regression on individual time series extracted for each pixel in the satellite images. The results of Eklundh and Olsson (2003) were later improved by Bolin et al. (2009) where a spatial model for the vegetation was used in the analysis to capture the spatial dependencies in the trend estimation.

To find regions where changes in the vegetation have occurred over the course of the studied time period, both Eklundh and Olsson (2003) and Bolin et al. (2009) used significance testing for the individual pixels in the field. Thus, pixels that individually had significant changes in vegetation were found, but no attempts were made to find simultaneous excursion regions. Here, we will use a similar model to that of Bolin et al. (2009) but also estimate joint excursion regions for the vegetation trends.

Assume that the vegetation measurements year $t$ are generated as,

Discussion

Estimating excursion sets and uncertainty regions for contour curves for stochastic fields are difficult problems, both because of computational issues but also because it might not be clear how such uncertainty regions should be defined. In this work, we have given precise definitions for these uncertainty regions, introduced the concept of excursion functions as a visual tool for illustrating the uncertainty in these regions, and presented a method for calculating these quantities for latent Gaussian models.

The main idea behind the computational method is to use a parametric family for the excursion sets in combination with a sequential integration method to reduce the computational effort required when estimating the sets in practice. Tests on simulated data showed that the method is accurate, and two applications were presented to show that the method is applicable even to large environmental problems.

There are a number of extensions that could be made to this work. First of all, using the one-parameter family for the excursion sets gives a method that falls into the broad category of $p$-value thresholding methods for estimating simultaneous excursion sets. As previously mentioned, the important advantage with the method proposed here compared with other commonly used thresholding methods is that the correct joint distribution is used when selecting the threshold. The disadvantage is that the method is computationally more expensive than many of the standard thresholding methods. It would, therefore, be interesting to do a comparison with other similar methods with respect to the accuracy and computational complexity. Another interesting comparison would be to compare the uncertainty regions for contour curves produced by these methods to those of Lindgren and Rychlik (1995). One could potentially also combine these methods with the work by Polfeldt (1999) to make statements on the quality of contour maps.

We also presented other parametric families that can be used to obtain more complicated methods for estimating the excursion sets, with the possibility of finding more precise estimates under the cost of higher computational complexity. Initial comparisons showed that there is not much gain in using these more complicated methods, but so far these comparisons have only been made using fairly simple latent models, and the gain is likely higher when the latent models are more complex. Hence, more studies are required to verify if this is the case and to investigate in what situations it is appropriate to use the simple one-parameter families. One possible advantage with the more complicated parametric families is when one has prior knowledge regarding the shape of the excursion sets. For example, if one knows that the excursion sets should be large contiguous regions, such knowledge could be incorporated using a two-parameter smoothing family.

As a final note, the method introduced here is available as a C-package with interfaces to both R and Matlab, see the supplementary material for the online version of the article for details.

Acknowledgements

Data used by the authors in the Sahel vegetation study include data produced through funding from the Earth Observing System Pathfinder Program of NASA’s Mission to Planet Earth in cooperation with National Oceanic and Atmospheric Administration. The data were provided by the Earth Observing System Data and Information System (EOSDIS), Distributed Active Archive Center at Goddard Space Flight Center which archives, manages, and distributes this data set. The data used in the PM10 study was provided by the information system Aria Web Regione Piemonte and Arpa Piemonte. The authors are grateful to Johan Lindström and Daniel Simpson for valuable discussions on the subject of excursions and contour curve uncertainty sets, and to Peter Guttorp for highlighting the need for a thorough treatment of the subject.

Appendix A Notes on the MCMC algorithm used in Example 3

References