# Exceedances and their uncertainty quantification

$\newcommand{\bm}{\mathbf} \newcommand{\sbm}{\boldsymbol}$

## Motivation

Greater awareness of human impact on the environment (e.g., IPCC, 2013) has led to an increase in the use of statistical modelling for predicting environmentally sensitive scientific processes of interest. For example, climate models help us to understand Earth's atmosphere, particularly those aspects associated with global climate change, such as the carbon cycle. Predictions are often made for the mean or median of a spatial process but, for risk assessment, strategic planning, and policy development, prediction of "extremes," or "exceedances" are also of interest; see Figure 1 for an example of rainfall-related exceedance probabilities over Australia.

Here, our research in environmental informatics is concerned with the problem of simultaneously predicting exceedances for a spatially dependent scientific process at many locations concurrently. This entails two challenges: Predicting the location, extent, and possibly the persistence of the exceedance, and then quantifying the prediction uncertainty; see here for a description of our approach to uncertainty quantification for hierarchical models. Further research in this area is given by Cressie and Suesse (2020).

Figure 1. The chance (percentage probability) of exceeding the median rainfall for October to December, 2015. The probabilities are generated from the Predictive Ocean Atmosphere Model for Australia (POAMA), the Bureau's dynamical climate model. Image from Australian Bureau of Meteorology, retrieved 23 Oct, 2015.

#### Definition of an exceedance

Consider a scientific process $Y(\bm{s};t)$, where $\bm{s}$ represents location and $t$ represents time. An exceedance region (a "hotspot" or a "coldspot") is defined as the set of locations $\bm{s}$ in a region $\mathcal{D}$ at time $t$, for which $Y(\bm{s};t)$ is above or below a given threshold, $u$.

Define a 'minimum' exceedance region ('coldspot'), $e^-(u)$, as the set of locations for which the spatio-temporal process $Y(\bm{s}; t)$ is lessthan $u$: $$e^{-}(u) \equiv \{ \bm{s}\in D : Y(\bm{s}; t) < u \};$$ and define a 'maximum' exceedance region ('hotspot'), $e^+(u)$, as the set of locations for which the spatio-temporal process $Y(\bm{s}; t)$ is greater than $u$: $$e^+(u) \equiv \{ \bm{s}\in D : Y(\bm{s}; t) > u \}.$$

### Exceedances for CO$_2$ fluxes

One goal of climate models is to predict spatio-temporal changes in atmospheric carbon dioxide (CO$_2$) concentration, known as CO$_2$ flux $(\Delta_{CO2}(\bm{s}; t))$, and hence to predict regions of large positive CO$_2$ fluxes (sources), large negative CO$_2$ fluxes (sinks), and the movement of CO$_2$ between these regions. This is a challenging problem because sources and sinks are dynamic, and CO$_2$ flux is expensive and difficult to measure directly.

To predict CO$_2$ flux fields, climate models combine (1) predicted total-column atmospheric CO$_2$ concentrations (denoted $\hat{X}_{CO2}(\bm{s}; t)$), which are frequently obtained from remote sensing data (see Statistical remote sensing for a description of $\hat{X}_{CO2}(\bm{s}; t)$ prediction and our research in this area), with (2) the predicted creation, movement, and transformation of CO$_2$ in the atmosphere, which are obtained from highly complex, often deterministic chemical transport models such as Geos-Chem or LMDZ, and (3) a set of assumptions about the atmospheric state. See Regional flux inversion for a methane-based example of our research on predicting flux fields.

In our research, we use the resulting predicted flux field $\{\hat{\Delta}_{CO2}(\bm{s}; t): \bm{s} \in \mathcal{D}, t=1,2,...\}$ to predict the location, extent, and persistence of CO$_2$ exceedances in the atmosphere. For example, in Figure 2, a sink-type, or 'minimum,' exceedance region is defined as the locations for which the predicted CO$_2$ flux is less than $u=-6.59$ gC/m$^2$/day.

Figure 2. Predicted CO$_2$ flux field (gC/m$^2$/day) obtained from the inversion of a simulated prior flux field and atmospheric measurements (left panel); and the same predicted flux field with $e^-(u=-6.59)$ marked in green (right panel). Original flux predictions courtesy of F. Chevallier.

#### Exceedance predictions and their uncertainty

We wish to predict an exceedance for a spatio-temporal process $Y$, such as $\Delta_{CO2}(\bm{s}; t)$. However, due to measurement error and incomplete sampling, we cannot observe $Y$ directly. Instead we observe a prediction (or estimate) or a proxy for $Y$, such as a matrix of predicted flux fields obtained from satellite remote sensing data and a transport model. Define $\bm{Z} = (\hat{\Delta}_{CO2}(\bm{s}_i; t_j): i=1,\ldots,I; j=1,\ldots,J)'$ to be the predicted CO$_2$ flux field.

For a given level of statistical significance, $\alpha$, and a given exceedance threshold, $u$, we consider a $(1 - \alpha)$100% exceedance set for $e^{+}(u)$ or for $e^{-}(u)$ (henceforth, collectively referred to as $e(u)$), based on $\bm{Z}$, which we write as $e(\bm{Z},u; 1-\alpha)$.

To find the $(1 - \alpha)$100% 'minimum' exceedance set, we assume that all locations are NOT exceedances and, using a test statistic $T(\bm{s})$, we test the hypotheses, $H_0: Y(\bm{s}) \ge u$ versus $H_1: Y(\bm{s}) < u$, for each $\bm{s} \in \mathcal{D}$. The locations, $\bm{s}$, for which we reject $H_0$ are predicted exceedances, whilst other locations are not. In general, the predicted exceedance set $e(\bm{Z},u; 1-\alpha)$ is contained within the true exceedance set $e(u)$ with probability $(1-\alpha)$.

Methods for predicting exceedance regions must consider both the spatial dependence in the data and the multiple-testing problem. Prediction methods that account for spatial dependence include kriging-based approaches (e.g., Cressie, 1993), and hierarchical-model-based approaches (see here for a description of these). Alternatively, methods based on extreme value theory (Gumbel, 1958; Kharin and Zwiers, 2005; Cooley et al., 2007), max stable processes (Schlather and Tawn, 2003; Cooley et al., 2007), and Bayesian quantile regression (Jeffreys, 1961; Reich et al., 2011) are also used (see Davison et al., 2012, for a recent review of these methods).

The control of hypothesis-testing errors for large multivariate datasets is presently an active area of research. Historically, methods to control Type 1 error, such as the Bonferroni adjustment (Dunn,1961), were developed, and applied to spatial data (Adler, 1981). More recently, methods based on the False Discovery Rate (FDR) (Benjamini and Hochberg, 1995; Storey, 2003; Storey and Tibshirani, 2003), or the False Nondiscovery Rate (e.g., Craiu and Sun, 2008) have been developed and extended to correlated data (Benjamini and Yekutieli, 2001), and spatial data (Shen et al., 2002; Benjamini and Heller, 2007; Brown et al., 2014; and Sun et al., 2015). To account for multiple testing, recent approaches to prediction of exceedance regions have used a simulation-based approach (French and Sain, 2013), or a Bayesian approach (Bolin and Lindgren, 2015).

To obtain a $(1-\alpha)100\%$ predicted exceedance set, we use a simulation-based hypothesis-testing approach based on that used by French and Sain (2013). This approach applies fast, desktop computing to calculate the critical value used to reject $H_0$, and multiple testing at many locations is automatically accounted for.

### Exceedance predictions for CO$_2$ fluxes

In our research, we have used predicted global flux fields from an observing system simulation experiment (OSSE), courtesy of Frederic Chevallier (Chevallier et al., 2005; Chevallier et al., 2007; Chevallier et al., 2010; Chevallier et al., 2014). Steps used in the OSSE to produce the predicted fluxes are given by Chevallier et al. (2007) as follows:

• Generate a set of pseudo-observations using a climatology of surface fluxes as boundary conditions to a transport model, such as LMDZ.
• Perturb the pseudo-observations consistently with the assumed error statistics.
• Perturb the surface flux climatology consistently with the assumed error statistics.
• Perform Bayesian inversion using the perturbed pseudo-observations as data and the perturbed climatology as the "prior field" to obtain a predicted flux field.
• Compare the predicted flux field to the known surface flux climatology to obtain error statistics.

From the OSSE, predicted flux fields for a global $97 \times 96$ grid (where the area of each grid cell depends on location), have been summarised for daytime (8am - 8pm local time) and nighttime (8pm - 8am local time), for four days per month during the period between January 2009, and January 2011. We fix a time $t$, and hence we drop the index $t$ from the notation. For the grid cell located at $\bm{s}$, a hierarchical-modelling approach was used to obtain a prediction, $\hat{Y}(\bm{s})$, for the hidden flux field, ${Y}(\bm{s})$, and the hypothesis tests described above were conducted using a test statistic of the form: \begin{equation*} T(\bm{s}) = \frac{\hat{Y}(\bm{s}) - u}{\sqrt{var(\hat{Y}(\bm{s}))}}. \end{equation*} Figure 3 gives an example of a 90% predicted exceedance set, for selected flux fields from the OSSE output, obtained using the simulation-based approach outlined above.

Figure 3. A 90% predicted exceedance set, $e(\bm{Z},u; 1-\alpha)$, based on the CO$_2$ flux field (gC/m$^2$/day) given in Figure 2, left panel. Notice that this is different from $e(u)$ shown in Figure 2, right panel, since it satisfies, $\text{Pr}(e(\bm{Z},u; 1-\alpha) \subset e(u))=(1-\alpha)100\%$. Original flux predictions courtesy of F. Chevallier.

Adler, R. J. (1981). The Geometry of Random Fields. Wiley, Chichester, UK. Reprinted by SIAM, 2010.

Benjamini, Y. and Hochberg, Y. (1995). Controlling the false discovery rate: A practical and powerful approach to multiple testing. Journal of the Royal Statistical Society, Series B57, 289-300.

Benjamini, Y. and Heller, R. (2007). False discovery rates for spatial signals. Journal of the American Statistical Association102, 1272-1281.

Benjamini, Y. and Yekutieli, D. (2001). The control of the false discovery rate in multiple testing under dependency. Annals of Statistics29, 1165-1188.

Bolin, D. and Lindgren, F. (2015). Excursion and contour uncertainty regions for latent Gaussian models. Journal of the Royal Statistical Society, Series B77, 85-106.

Brown, D. A., Lazar, N. A., Datta, G. S., Jang, W., and McDowell, J. E. (2014). Incorporating spatial dependence into Bayesian multiple testing of statistical parametric maps in functional neuroimaging. NeuroImage84, 97-112.

Chevallier, F., Fisher, M., Peylin, P., Serrar, S., Bousquet, P., Bréon, F.-M., Chédin, A., and Ciais,P. (2005). Inferring CO$_2$ sources and sinks from satellite observations: Method and application to TOVS data. Journal of Geophysical Research110, D24309, doi:10.1029/2005JD006390.

Chevallier, F., Bréon, F.-M., and Rayner, P. J. (2007). Contribution of the Orbiting Carbon Observatory to the estimation of CO$_2$ sources and sinks: Theoretical study in a variational data assimilation framework. Journal of Geophysical Research112, D09307, doi:10.1029/2006JD007375.

Chevallier, F., Feng, L., Bösch, H., Palmer, P. I., and Rayner, P. J. (2010). On the impact of transport model errors for the estimation of CO$_2$ surface fluxes from GOSAT observations. Geophysical Research Letters37, L21803, doi:10.1029/2010GL044652.

Chevallier, F., Palmer, P. I., Feng, L., Boesch, H., O'Dell, C. W., and Bousquet, P. (2014). Toward robust and consistent regional CO$_2$ flux estimates from in situ and spaceborne measurements of atmospheric CO$_2$. Geophysical Research Letters41, 1065-1070, doi:10.1002/2013GL058772.

Cooley, D., Nychka, D., and Naveau, P. (2007). Bayesian spatial modeling of extreme precipitation return levels, Journal of the American Statistical Association102, 824-840.

Craiu, R. V. and Sun, L. (2008). Choosing the lesser evil: trade-off between false discovery rate and non-discovery rate. Statistica Sinica18, 861-879.

Cressie, N. (1993). Statistics for Spatial Data (revised edition). Wiley, New York, NY.

Cressie, N. and Suesse, T. (2020). Great expectations and even greater exceedances from spatially referenced data.

Davison, A. C., Padoan, S. A., and Ribatet, M. (2012). Statistical modeling of spatial extremes. Statistical Science27, 161-186.

Dunn, O. J. (1961). Multiple comparisons among means. Journal of the American Statistical Association56, 52-64.

French, J. P. and Sain, S. R. (2013). Spatio-temporal exceedance locations and confidence regions. Annals of Applied Statistics7, 1421-1449.

Gumbel, E.J. (1958). Statistics of Extremes. Columbia University Press, New York, NY. Reprinted by Dover Press, 2004.

IPCC (2013). Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Stocker, T. F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V. and Midgley, P. M. (eds.). Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 1535 pp. Obtained from http://www.ipcc.ch.

Kharin, V. and Zwiers, F. (2005). Estimating extremes in transient climate change simulations. Journal of Climate18, 1156-1173.

Reich, B. J., Fuentes, M., and Dunson, D. B. (2011). Bayesian spatial quantile regression. Journal of the American Statistical Association106, 6-20.

Schlather, M. and Tawn, J. (2003). A dependence measure for multivariate and spatial extremes: Properties and inference. Biometrika90, 139-156

Storey, J. D. (2003). The positive false discovery rate: A Bayesian interpretation and the q-value. Annals of Statistics31, 2013-2035.

Storey, J. D. and Tibshirani, R. (2003). Statistical significance for genome-wide studies. Proceedings of the National Academy of Sciences of the United States of America100, 9440-9445.

Sun, W., Reich, B. J., Cai, T. T., Guindani, M., and Schwartzman, A. (2015). False discovery control in large-scale spatial multiple testing. Journal of the Royal Statistical Society, Series B77, 59-83.

Authored by S. Burden, 2015.