{'fraction_non_alphanumeric': 0.03484051050091641, 'fraction_numerical': 0.017302550823090245, 'mean_word_length': 4.569582318786289, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spatio-temporal exceedance locations and confidence regions

An exceedance region is the set of locations in a spatial domain where a process exceeds some threshold. Examples of exceedance regions include areas where ozone concentrations exceed safety standards, there is high risk for tornadoes or floods, or heavy-metal levels are dangerously high. Identifying these regions in a spatial or spatio-temporal setting is an important responsibility in environmental monitoring. Exceedance regions are often estimated by finding the areas where predictions from a statistical model exceed some threshold. Even when estimation error is quantifiable at individual locations, the overall estimation error of the estimated exceedance region is still unknown. A method is presented for constructing a confidence region containing the true exceedance region of a spatio-temporal process at a certain time. The underlying latent process and any measurement error are assumed to be Gaussian. Conventional techniques are used to model the spatio-temporal data, and then conditional simulation is combined with hypothesis testing to create the desired confidence region. A simulation study is used to validate the approach for several levels of spatial and temporal dependence. The methodology is used to identify regions of Oregon having high precipitation levels and also used in comparing climate models and assessing climate change using climate models from the North American Regional Climate Change Assessment Program.

1. Introduction. Identifying regions of extreme or unusual response is often a vital concern when analyzing environmental data. Identification of these regions can have important health, societal, and political impacts since these areas may indicate unusual events, outbreak of disease, toxic conditions, regions of extreme risk, etc. [Patil (2010)]. These regions of interest are often described as exceedance regions or hotspots, and are defined as the set of locations where the response process of interest exceeds some specified threshold. On the basis that the broad definition of exceed is "to go beyond the bounds or limits of" something [Dictionary.com (2012)], an exceedance region may refer to the area where the responses are above a threshold or the area where the responses fall below a threshold (though the two settings should be carefully distinguished).
The uncertainty associated with estimating exceedance regions has previously received relatively little attention. Our goal in this paper will be to introduce an approach for identifying exceedance regions with a quantifiable amount of certainty, and then using this methodology to confidently identify exceedance regions in two different contexts. First, in Section 4, this methodology will be used to confidently identify the regions of Oregon having extreme precipitation levels in October of 1998 using data from the previous two years. Being able to predict seasonal climate can help to minimize "climate surprises," reduce impacts on society and ecosystems, assess the chances of events like drought and wildfire, and make important economic decisions in fields such as agriculture, energy, insurance, and public health [NOAA (2011)]. In Section 5, we use the proposed methodology to explore the temperature data of climate models from the North American Regional Climate Change Assessment Program (NARCCAP). Discussion will include assessment of future climate change for several regions of North America based on these models and also highlight similarities and differences between the projections from various climate models.
The methodology proposed in this article allows one to construct a confidence region for the exceedance region of a spatio-temporal process at a certain time. The resulting confidence region will contain the entire exceedance region with known confidence. To our knowledge, no methods have previously been available for constructing confidence regions containing the entire exceedance region of a spatial/spatio-temporal process and having the desired coverage properties. The methodology combines optimal spatiotemporal prediction methods and a hypothesis-testing-like approach to construct the confidence regions. Due to the simultaneous nature of the inference, the size and shape of the confidence regions will be directly related to the domain of interest. French (2012) introduced an approach for identifying the level curves of a spatial process with known confidence. Our methodology generalizes that of French (2012) by handling exceedance regions for spatio-temporal processes and also accounting for additional uncertainty when estimating the mean structure of the spatio-temporal process. The exceedance region of a spatial or spatio-temporal process is often predicted as the region(s) where the predicted response exceeds the desired threshold. This approach results in a biased predicted exceedance region since spatial predictors tend to oversmooth the response surface [Zhang, Cressie and Craigmile (2008)]. Craigmile et al. (2005) accounted for the oversmoothing of traditional spatial predictors by using a loss function that assigned more weight to extreme values in prediction; the predicted exceedance regions resulting from this approach are often larger than those of traditional methods. Patil and Taillie (2004) used the upper level set (ULS) scan statistic, a modification of the popular spatial scan statistic [Kulldorff and Nagarwalla (1995); Kulldorff (1997)] to identify regions or clusters of cells with elevated responses when compared to neighboring cells. The spatial scan statistic is often used in the setting where a region is tessellated into cells and response data are available as counts. Patil and Taillie (2004) used maximum likelihood estimation to identify regions with elevated counts through exhaustive search of the parameter space. The ULS scan statistic reduced the computational complexity of their task by reducing the search space, which allowed Patil and Taillie to construct a type of confidence region for the upper level sets (hotspots) by finding all zones where the ULS statistic was not statistically significant. The Progressive Upper Level Set (PULSE) scan statistic was recently proposed by Patil, Joshi and Koli (2010) as a refinement of the ULS statistic for detecting geospatial hotspots. Neither of the previous two approaches incorporates the typical notion of spatial dependence used to model and predict spatial data as exposited by Cressie (1993), Schabenberger and Gotway (2005), and many others. Zhang, Cressie and Craigmile (2008) predicted the exceedance region by finding the region minimizing the posterior loss of an image-based loss function using simulated annealing. The predicted exceedance region found by the approach of Zhang, Cressie and Craigmile (2008) is statistically optimal, but is not a confidence set in the traditional statistical sense. Recently, Sun et al. (2012) developed a unified theoretical and computational framework for false discovery control in multiple testing when trying to identify locations where the mean response of a process exceeds some threshold.
Other subject areas with connections to exceedance regions include the estimation of spatial cumulative distribution functions (CDFs) and the analysis of fMRI experiments to identify activated voxels. The spatial CDF of a process over a domain is a random distribution function that provides a statistical summary of a random field over a given region [Lahiri (1999)]. Spatial CDFs were introduced by Majure et al. (1996) and Lahiri (1999), andLahiri et al. (1999) used subsampling approaches to predict a spatial CDF and study the asymptotic properties of the associated predictors. Zhu, Lahiri and Cressie (2002) developed an approach to detect a change in the spatial CDF of a region over time using the difference between two empirical CDFs and then quantified any change using the weighted integrated squared distance between the empirical CDFs. While spatial CDFs provide a way of identifying the proportion of a region exceeding some threshold, they do not identify where these exceedances occur. In the medical imaging field, there has long been an interest in identifying the areas of an fMRI map depicting brain activity in response to some stimuli. Marchini and Presanis (2004) provide a helpful overview of the various approaches to detecting active voxels. The most common approach for identifying active voxels is to declare a voxel active if some associated test statistic is above a single common threshold. This threshold is determined in a way that allows for control of some error criterion, such as the standard familywise error rate (FWER) or False Discovery Rate [FDR, Benjamini and Hochberg (1995)]. Often, the approach to determining this threshold utilizes the theory of random fields [Adler (2010), Adler and Taylor (2007)] and the expected Euler characteristic [Taylor, Worsley and Gosselin (2007), Adler (2008)] to determine the probability that the maximum test statistic is greater than some threshold. This approach can naturally apply in a geostatistical context such as the present one, but assumptions regarding the stationarity of the test statistics can be overly strict and often produces conservative results [Marchini and Presanis (2004)].
1.1. Outline. The approach proposed in this paper fills a methodological gap in the literature by allowing researchers to construct confidence regions for exceedance regions using traditional geostatistical tools and notions of spatial dependence while having the confidence level properties typically desired. The underlying latent process and any measurement error are assumed to be Gaussian. This paper will proceed in the following manner. In Section 2 we will describe the proposed methodology in detail. Section 3 presents the results of a simulation study used to validate our approach, as well as several case studies exploring the shape of the confidence region for various underlying mean structures. Sections 4 and 5 describe applications of this methodology in the climatological settings mentioned above. We conclude with a brief discussion in Section 6.

Methodology.
2.1. Framework. Consider a spatio-temporal response process Y (·, ·) ≡ {Y (s, t), (s, t) ∈ D × T ⊂ R 2 × R}, which has spatial index s contained in a bounded two-dimensional region of interest D, and time index t which is contained in a bounded set of possible times T ⊂ R. We assume that Y (·, ·) is the sum of a hidden spatio-temporal process Z(·, ·) and a measurement error process ε(·, ·), so that We will refer to Y (·, ·) as the observable process and Z(·, ·) as the hidden or latent process. The hidden process Z(s, t) is composed of a mean struc-ture µ(s, t) capturing large-scale response behavior and a mean-zero spatiotemporal process W (s, t) capturing small-scale behavior, such that Z(s, t) can be decomposed into The mean structure µ(s, t) is assumed to follow the linear model µ(s, t) = x(s, t) T β, where x(s, t) is a k × 1 vector of known space-time covariates for location s at time t, β is a vector of space-time trend parameters, and W (s, t) is assumed to be a Gaussian process having continuous sample paths [cf. Paciorek (2003), page 57; Abrahamsen (1997), Adler and Taylor (2007)]. The covariance between two responses for the hidden process is denoted by cov{Z(s 1 , t 1 ), Z(s 2 , t 2 )} = C(s 1 , s 2 , t 1 , t 2 ).
The error process ε(s, t) is assumed to be a Gaussian white-noise process with mean 0 and variance σ 2 ε , and the covariance between two responses of the observable process is given by where v ε is a known function of the spatial and temporal indices (often returning 1 when s 1 = s 2 and t 1 = t 2 and zero otherwise).
Suppose we have observed n responses y = [Y (s 1 , t 1 ), . . . , Y (s n , t n )] T of a partial realization of Y (·, ·). Our goal is to use the observed responses to construct a confidence region for the locations where the hidden process Z(·, ·) exceeds some threshold u at some time t p . We define the exceedance region above a threshold u for the hidden process Z as E Z u + = {s ∈ D : Z(s, t p ) ≥ u}, while the exceedance region below a threshold u may be defined as E Z u − = {s ∈ D : Z(s, t p ) ≤ u}. Depending on the exceedance region of interest (above or below the threshold), we would like to find a region S u + that contains E Z u + with probability 1 − α, that is,

2.2.
Methodology. For simplicity, we will only describe methodology for determining the confidence region S u + (the methodology for determining S u − being analagous). Our approach for constructing our confidence region S u + will combine geostatistical techniques along with an approach similar to hypothesis testing. We will state a null and alternative hypothesis, determine a test statistic and critical value, then make our conclusion. However, an important difference is that these tests are related to random variables instead of parameters, violating the classical definition of a hypothesis test. In spite of this difference, we will use the traditional terminology of hypothesis testing since the same concepts apply. Additionally, the procedure will call for a large number of statistical tests so the critical value must be adjusted in order to control the confidence level of our confidence region. We note that though this methodology is presented in the context of spatio-temporal processes, the approach is equally applicable to purely spatial data (with no time varying component) simply by assuming we have only one time.
The confidence region S u + can be constructed by testing for each s ∈ D, H 0 : Z(s, t p ) = u versus H a : Z(s, t p ) < u on the basis of a test statistic Z ′ (s, t p ). The region S u + is made up of all locations s where we fail to conclude that Z(s, t p ) < u. Naturally, the test statistic Z ′ (s, t p ) is based on a predictor of Z(s, t p ),Ẑ(s, t p ), from the statistical model outlined in Section 2.1. When the trend parameter vector β must be estimated, the unbiased linear predictor of Z(s, t) minimizing the mean-squared prediction error is known as the universal kriging predictor. Suppose we wish to make a prediction of Z(·, ·) at location s 0 and time t 0 . Let x 0 denote the k × 1 vector of space-time covariates associated with Z(s 0 , t 0 ), X y denote the n×k matrix of covariates for the observed data y = {Y (s 1 , t 1 ), . . . , Y (s n , t n )}, Σ y denote the n × n covariance matrix for the observed data y, and c 0y = [C(s 0 , s 1 , t 0 , t 1 ), . . . , C(s 0 , s n , t 0 , t n )] T be the n × 1 vector of covariances between the response to be predicted and the observed data responses. The optimal predictor for Z(s 0 , t 0 ) iŝ y y is the generalized least squares estimator of β. The associated mean-squared prediction error (kriging variance) of this predictor is given by [Schabenberger and Gotway (2005), page 242]. Using the assumptions of normality made regarding Z(·, ·) and ε(·, ·), and using the properties ofẐ(·, ·) and σ k (·, ·), the quantitŷ has a standard Gaussian distribution. Assuming the null hypothesis is true, (3) leads us to the natural test statistic for location s 0 and time t 0 .

7
The confidence level of our confidence region is controlled by carefully using the duality that exists between confidence intervals and hypothesis testing. Our confidence region S u + is composed of all locations s where we fail to reject the null hypothesis that Z(s, t p ) = u [or perhaps more clearly, where we fail to reject that Z(s, t p ) ≥ u]. Thus, our confidence region will fail to contain the true exceedance region any time a Type I error is made in our hypothesis tests. In order to maintain a confidence level of 1 − α for our confidence region S u + , we must ensure that the probability of making a Type I error for all tests simultaneously is α. Let C α be the critical value at which we will reject H 0 and conclude H a . By definition, we can only make a Type I error at the locations s that are part of the true exceedance region. Consequently, when seeking to control our error rate, we are only concerned with controlling our error rate in the context of the possible realizations of E Z u + conditional on the observed data. Thus, to control the Type I error rate of our hypothesis at level α and maintain the confidence level of our confidence region at 1 − α, we should only conclude H a when Z ′ (s, 2.3. Estimating the critical value C α and constructing S u + . To properly estimate the critical value C α , we must be able to adequately approximate the distribution of within our domain D. Note that the randomness of (5) is driven by the uncertainty in E Z u + , not the uncertainty in {Z ′ (s, t p )}. Conditional on the observed data, {Z ′ (s, t p )} is completely determined. In contrast, we do not know the exceedance region E Z u + generated in the same realization as the observed data y. Though we are interested in the realization of E Z u + generated along with the actual data, there are infinitely many possible exceedance regions compatible with the observed data. Thus, we need to consider the distribution of inf s∈E Z u + {Z ′ (s, t p )} for all possible exceedance region realizations compatible with the observed data.
In practice, to estimate C α in (4), we must assume that the behavior of the hidden process Z(·, t p ) over the continuous domain D can be adequately approximated by considering its behavior over an appropriately sized finite regular grid [cf. Zhang, Cressie and Craigmile (2008)]. We discretize the domain D into m pixels, letting G = {g 1 , . . . , g m } be the set of pixels and s G = {s * 1 , . . . , s * m } be the set of midpoints of these pixels. Let z G = [Z(s * 1 , t p ), . . . , Z(s * m , t p )] T denote the vector of hidden responses on the grid of locations s G . To approximate possible exceedance region realizations conditional on the observed data, we can simulate realizations of z G |y and determine the (discretized) exceedance regions for each realization. The joint distribution of z G and y is multivariate normal, making it straightforward to determine that the conditional distribution z G |y is also multivariate normal with closed form expressions for the associated mean and covariance matrix [Wasserman (2004), Theorem 2.44]. It is easy to simulate directly from a multivariate normal distribution using a decomposition of its covariance matrix [Givens and Hoeting (2005), page 146], but this does not take into account the fact that our mean function is being estimated. Instead, we will generate realizations of z G |y by conditioning our simulation using kriging. In essence, this approach produces realizations of the conditional distribution by perturbing the kriging estimates by possible realizations of the associated kriging error, as described below.
Suppose that recalling that Σ y is the n × n covariance matrix of the observed data, and letting Σ G denote the m × m covariance matrix of z G and Σ Gy denote the m × n cross-covariance matrix between y and z G . The joint distribution of y c and z c has the same covariance as the joint distribution of y and z G . To obtain B realizations of z G |y, we first use an extension of (1) allowing for simultaneous predictions at multiple locations [Schabenberger and Gotway (2005), pages 242-243] to find the n × m weight matrix Λ such thatẑ G = [Ẑ(s * 1 , t p ), . . . ,Ẑ(s m , t p )] T = Λ T y. Next, we simulate B realizations from the distribution in (6). Denote the first n elements of the ith realization by y (i) c and the next m elements by z (i) c . The ith simulated realization of z G |y is then obtained using the expressioñ T denoting the individual responses of this realization. The first part of expression (7) is the kriging predictor and the second part is a realization of possible kriging error. Additional information about this simulation process can be found in Chilès and Delfiner [(1999), pages 465-468].
To estimate C α , we begin by generating B realizations {z (1) , . . . ,z (B) } of the conditional process z G |y. For realization i, we identify EZ the union of the pixels where the simulated discretized conditional process exceeds the threshold u. For each realization of EZ (i) u + , we next consider the minimum of the test statistic Z ′ (·, t p ) at the gridded locations in the realized exceedance region under the assumption that the null hypothesis is true. Specifically, for each realization of EZ Last,Ĉ α , our estimated value of C α , can be obtained by finding the α quantile of the set of minima from the previous step. The confidence region S u + (or at least its discretized version) is then

Details of inference.
We briefly summarize the inferences that can be made using the proposed methodology, as well as some details related to this inference.
2.4.1. Confidence region above. Using the methodology described above, one may construct a confidence region S u + containing the realized exceedance region E Z u + with probability 1 − α. Specifically, the probability that the confidence region S u + produced by our method will contain the realized ex- Note that the confidence properties of S u + are in the frequentist paradigm. Thus, the probability in (8) may refer to the probability of S u + containing the realized exceedance region before we collect the data or the long-term relative frequency of containing the realized exceedance region when applying this methodology to independent trials. A specific confidence region will either contain or not contain the true but unobserved realization of E Z u + , but in repeated application of this methodology to new data sets, the long-term proportion of trials where the realized exceedance region will be entirely contained in the confidence region is 1 − α.

2.4.2.
Confidence region below. One may obtain a confidence region S u − for the exceedance region below a threshold, E Z u − , by considering the inverse of the problem described in Sections 2.2 and 2.3. By testing for all s ∈ D, H 0 : Z(s, t p ) ≤ u versus H a : Z(s, t p ) > u [and modifying the steps in Section 2.3 so that we approximate the critical value C α such that This type of inference will be useful when one desires a confidence region for the entire area where a process may fall below some level, for example, identifying all locations where drought may occur. The confidence properties of S u − are the same as those described for S u + in Section 2.4.1.

2.4.3.
Properties of the complement. The inference related to S u + and S u − is useful when one desires to find a region containing an entire exceedance region with high confidence. However, it is possible that a sizable portion of the confidence region is not part of the exceedance region and that the confidence region is much larger than the true exceedance region. Instead, we may be interested in the region of our domain D where we can be confident that a response will exceed the threshold in question. This type of confidence region is easily obtained as a byproduct of the construction of S u + and S u − .
Suppose we wish to find a confidence set contained entirely in E Z u + , that is, every point in this confidence set has a response exceeding the threshold u. When constructing the confidence region S u − for the locations having a response below the threshold u, S u − has the property that P In other words, we can be confident that every point in S c u − simultaneously has a response exceeding the threshold u. Similar to the properties of the rejection region S c u − , the rejection region S c u + has the property that every point in it will fall below the threshold u with confidence.
By combining the information from S c u − with S u + , one may obtain a liberal and conservative view of where the realized exceedance region E Z u + may be located. Specifically, we can be confident (in a frequentist sense) that every point in S c u − is a member of the true exceedance set E Z u + , but it is likely that S c u − is smaller than the true exceedance set. On the other hand, we can be confident that every location in E Z u + will be contained within S u + , but S u + will likely be larger than E Z u + . A similar interpretation applies to the regions S c u + and S u − in relation to E Z u − . Table 1 summarizes the methodology and inferences resulting from the proposed methodology.
2.4.4. Additional details. We note once again that the inference and properties of confidence regions discussed above are directly related to the domain in question. Since the procedure described in Section 2.3 depends on discretizing the domain of interest D, determining the realized exceedance regions EZ (i) u + over that domain, estimating C α over that domain, and evaluating the hypothesis tests over that domain, the confidence regions are directly linked to the domain chosen. This must be kept in mind when deciding on a domain of interest and considering the scope of inference.
We also note that the confidence properties above are dependent on knowing the true covariance properties of the random field under consideration. Table 1 Summary of procedure and inferences for proposed methodology

Exceedance above Exceedance below
Region In practice, this is rarely the case. While it is typical to simply plug-in the estimated covariance parameters as if they were the truth [Schabenberger and Gotway (2005), page 254], this will likely result in confidence regions that have confidence levels differing from the desired levels. A Bayesian approach to this problem that naturally incorporates the uncertainty in the covariance parameters is currently under development.
3. Simulation studies. In this section we will use simulation experiments of the proposed methodology to compare empirical confidence levels to intended confidence levels for scenarios involving several combinations of mean structures and spatio-temporal dependence. In the experiments, the goal is to construct a confidence region S u + containing the exceedance region E Z u + of the realization in question. The experiments are similar to ones performed by Zhang, Cressie and Craigmile (2008). Discussion will follow regarding how estimating the covariance parameters affects the empirical confidence levels and how the covariance parameters affect the shape and size of the confidence region for the true exceedance set. The analysis in Sections 3, 4, and 5 was performed using R version 2.15.1 [R Core Team (2012)].
3.1. Overall structure. In each of the subsequent experiments, the hidden process Z(·, ·) was assumed to follow the same general form The mean of Z(s, t) was allowed to vary between four different patterns: trend, cone, cup, and waves. The domain D of the experiment depended on the pattern under consideration. The mean structures of these patterns are shown in Figure 1, while the corresponding domain D of the experiment, the covariate vector x(s, t), and coefficient vector β of each pattern are provided in Table 2. The spatio-temporal process W (s, t) was multivariate normal with mean 0 and isotropic covariance function where · is the Euclidean distance metric. This covariance function is obtained by multiplying an exponential covariance function exp( s 1 − s 2 /φ) by ρ |t 1 −t 2 | , the correlation function of an AR(1) temporal process. The strength of the spatial dependence is determined by φ, and the strength of the temporal dependence is governed by ρ. In all experiments, the scale parameter σ 2 W was fixed to be 1, while the spatial dependence parameter φ varied between the values 0.5 (weak spatial dependence), 1.5 (moderate spatial dependence), and 5 (very strong spatial dependence), and the temporal dependence parameter ρ varied between the values 0.1 (weak temporal dependence), 0.5 (moderate temporal dependence), and 0.9 (strong temporal dependence). Experiments were performed assuming no measurement error (σ 2 ε = 0) and when measurement error was present (σ 2 ε = 0.1 or 0.5 depending on the mean structure).
In each experiment, the temporal window of interest T = {1, 2, 3, 4}. It was assumed that responses with time index t = 1, 2, 3 were observed data and that responses with time index t = 4 corresponded to future responses. For each of the first three times, responses were observed at the same 100 locations s 1 , s 2 , . . . , s 100 . The sites of the 100 locations were irregularly spaced and were obtained by drawing independent observations from Table 2 The domain D, covariates vector x(s, t), and coefficient vector β used for the experiments of each pattern shown in Figure 1. The individual spatial coordinates of a spatial location s are denoted sx and sy Pattern D x(s, t) β a Uniform(0, 1) distribution and combining them to form a set of spatial coordinates. To apply the methodology described in Section 2 and obtain "future" data to test empirical confidence levels, the spatial domain for each experiment was discretized into a 50 × 50 grid of pixels. year, future responses were "observed" at all locations s ∈ s G to provide the future test data. For each combination of covariance parameters φ and ρ, the hidden process Z(·, ·) was randomly generated at each of the 100 locations for the first three times, and then at the pixel center points for the fourth time. After generating the hidden process Z(·, ·), independent error values ε 1 , ε 2 , . . . , ε 300 were generated according to a N (0, σ 2 ε ) distribution and added to the 300 hidden process values observed for the first three times to obtain the observed response values y = {Z(s 1 , 1) + ε 1 , . . . , Z(s 100 , 1) + ε 100 , Z(s 1 , 2)+ε 101 , . . . , Z(s 100 , 2)+ε 200 , Z(s 1 , 3)+ε 201 , . . . , Z(s 100 , 3)+ε 300 }. Note that the observed process recovers the hidden process when σ 2 ε = 0. The threshold level u used to create the exceedance set E Z u + was the 90th percentile of the data generated on the 50 × 50 grid for time t = 4. Confidence regions S u + were constructed at confidence levels of 0.90 and 0.95. For each experiment, the critical value C α was estimated using 2000 realizations of z G |y. Empirical confidence levels were calculated by generating 200 independent realizations of Z(·, ·) for each experimental setting, using the procedure of Section 2 to construct the confidence region S u + for the exceedance region of the hidden process at time t = 4, and then finding the proportion of realizations in which the confidence region S u + contained the exceedance region E Z u + .

Empirical confidence levels.
For the trend pattern, experiments were run for all 9 combinations of the two dependence parameters φ and ρ. Experiments were run both without measurement error (σ 2 ε = 0) and with measurement error (σ 2 = 0.1). In the initial experiments, it was assumed that the covariance parameters were all known. Subsequently, the same experiments were performed assuming that the covariance parameters were unknown. In the second set of experiments, the covariance parameters were estimated using restricted maximum likelihood (REML). The measurement error variance σ 2 ε was estimated in each of these experiments, regardless of whether measurement error was actually present.
3.2.1. Results when covariance parameters known. The results of the simulation studies for the trend mean structure when the covariance parameters were assumed known are given below in Table 3. A comparative boxplot of the empirical confidence levels for these experiments grouped by desired confidence level is shown in Figure 2(a). The empirical confidence levels are all fairly close to the desired confidence level, and the empirical confidence levels behave in a manner consistent with sampling variability.

3.2.2.
Results when covariance parameters unknown. The covariance parameters of our data will not be known in practice, so we must use estimated values when making predictions. The empirical confidence level of the confi- dence regions of Section 2 are likely to be lower when estimated covariance parameter values are used instead of the true covariance parameter values.
To get an idea of the magnitude of this drop, the simulation experiments for the trends mean structure of the previous section were repeated using estimated covariance parameters. Utilizing the same random number seed in each set of experiments, the observed data for the two sets of experiments was the same. The only difference in inference was that the covariance parameters were estimated via restricted maximum likelihood (REML) in the second set of experiments and then used in the construction of the confidence regions. Consequently, the results from the two experiments are paired. The differences between the empirical confidence levels for the two sets of experiments are displayed in the boxplot in Figure 2(b). The difference between the empirical confidence level for the two settings is typically between 0.01-0.07, though the values sometimes deviated from this range.

Computational cost.
The proposed methodology has a computational cost essentially the same as that of the conditional simulation algorithm. As pointed out in the discussion of (7), the predicted responses and prediction weights are reused in the conditional simulation process. The only additional cost is the simulation of the unconditional random field in (6), and this cost can vary greatly depending on the algorithm selected. For the simulations described above, the prediction, conditional simulation, and construction of the confidence regions averaged just under 40 seconds on a MacBook Pro running OS X 10.7.4 with a 2.53 GHz Intel Core 2 Duo CPU and 4 GB of RAM. The same process averaged just over 3 seconds on a desktop PC running Fedora Core 16 with a 3.33 GHz Intel Core i7 CPU and 24 GB of RAM since many of the computations can run in parallel. The memory needed to construct the confidence regions depends directly on the number of observed and predicted data values and depends on the implementation of the proposed methodology. The conditional simulation process can be performed one at a time if necessary to minimize memory usage. However, for the simulation experiments performed in this section, the objects directly related to prediction, conditional simulation, and construction of the confidence regions required just over 38 MB of RAM.

Shapes.
A question of interest is how the shape of the confidence region for the exceedance set is affected by factors such as the mean structure, spatio-temporal dependence, and error variance. We briefly studied the effects of these factors graphically using the cone, cup, and waves patterns discussed in Section 3.1. For these patterns, only the combinations of (φ, ρ) equal to (0.1, 0.1), (0.5, 0.5), (0.9, 0.9) were tested; we refer to these combinations of (φ, ρ) as weak, medium, and strong spatio-temporal dependence, respectively. The covariance parameters were assumed known in these additional experiments. In order to evaluate the the general shape of the confidence and exceedance regions for each mean structure, the 200 confidence regions and realized exceedance regions from each simulation experiment were used to construct a "median" exceedance region and confidence region. The "median" in each case was determined by noting the pixels where the realized exceedance region and/or confidence region appeared in at least half of the 200 realizations. Similarly, to evaluate how the strength of spatiotemporal dependence and measurement error affect the size of a confidence region, the number of pixels used to construct each confidence region for the 200 realizations of each experiment was also noted. The shapes of the confidence and exceedance regions for the different mean structures conformed to intuitive expectations. Specifically, the confidence and exceedance regions of each shape patterned the mean structure. For example, the confidence and exceedance region for the cone mean structure were roughly circular since the exceedance region of the cone would simply be the upper portion of the cone. As an example of this, the median confidence region and exceedance set for the cone mean structure having medium spatial and temporal dependence without measurement error are shown in Figure 3(a). In order to assess the relationship between the strength of spatio-temporal dependence and the size of the resulting confidence region, as well as how measurement error affected these regions, we looked at the number of pixels in the confidence region from each of the 200 realizations of each experiment. Intuitively, stronger spatio-temporal dependence would lead to smaller predictive uncertainty, leading to a better estimate and smaller confidence region for the realized exceedance region. This pattern was consistently seen across all mean structures, and an example of this is shown in the boxplots in Figure 3(b) for the waves mean structure. Last, we would expect that the addition of measurement error to the observed responses would increase the size of the confidence regions due to greater uncertainty in predictions. This pattern was seen across all mean structures and levels of dependence, and an example of this for the cups mean structure is shown in Figure 3(c).

Case study 1: Precipitation in Oregon.
We demonstrate application of this methodology using precipitation data from the state of Oregon. Our goal is to identify with 90% confidence the regions of Oregon where the total monthly precipitation exceeds 250 mm (approximately the 96th percentile of the observed data) in October of 1998 using the precipitation measurements available in October of 1996 and 1997. Analysis was performed using the raw data available at http://www.image.ucar.edu/Data/US.monthly.met/ FullData.shtml#precip. The data used in this case study were created from the data archives of the National Climatic Data Center and have previously been used to generate high-resolution maps of precipitation and other meteorological variables [Daly et al. (2001), Johns et al. (2003)]. The raw data provides the longitude (degrees), latitude (degrees), and elevation (m) for 11,918 unique sites throughout the United States and, when available, includes the total monthly precipitation (mm) for each month between the years 1895-1997. There were 447 total observations provided for the state of Oregon between 1996Oregon between and 1997Oregon between , with 191 observations in 1996Oregon between and 256 observations in 1997. The observed data locations are indicated by solid dots on a map of Oregon in Figure 4  Exploratory analysis of the precipitation measurements revealed them to be positively skewed, so a square-root transformation of the measurements was taken to achieve approximate normality. An image plot of average precipitation for the transformed measurements is shown in Figure 4(b). Because the data were irregularly spaced, the bicubic spline interpolation algorithm of Akima (1996) was used to interpolate the average of the square root of the total monthly precipitation for the months of October 1996 and 1997 onto a regular grid before constructing the image plot.
The mean precipitation level of the responses Z(s, t) was assumed to be the same for each year and follow the linear structure µ(s, t) = x(s, t) T β, with the covariates vector x(s, t) = [1, longitude, latitude, elevation] T . The covariance function of the hidden data was modeled as being stationary and fully separable with respect to space and time, while allowing for a spatiotemporal measurement error effect for the observed data that varied by year. Consequently, the covariance function of the observed data may be written as is a purely spatial covariance function, C T (·) is a purely temporal covariance function, and C ε (·, ·) is the covariance function of the measurement errors. We assumed that the spatial covariance C S (·) could be modeled using an isotropic, Matérn covariance model of the form C S (s 1 , s 2 ) = σ 2 2 1−ν (h/φ) ν K ν (h/φ)/Γ(ν), where h = s 1 − s 2 is the Euclidean distance between two data locations, φ is a parameter related to the strength of spatial dependence, σ 2 measures the variance of the hidden process Z(s, t), ν controls the smoothness of the hidden process, and K ν (·) is the modified Bessel function of the second kind of order ν. The temporal covariance C T (t 1 , t 2 ) was modeled using the correlation function of an AR(1) process so that C T (t 1 , t 2 ) = ρ u , where u = |t 1 − t 2 | is the time lag between t 1 and t 2 and ρ measures the strength of spatial dependence. Last, the covariance function of the measurement errors can be written as C ε (s 1 − s 2 , t 1 − t 2 ) = [σ 2 96 1 (96) (t 1 ) + σ 2 97 1 (97) (t 1 )]1 (0,0) (h, u), where σ 2 96 and σ 2 97 are the variances of the measurement errors in years 1996 and 1997, respectively. Restricted maximum likelihood estimation was used to estimate the covariance parameters θ = [σ 2 , φ, ν, ρ, σ 2 96 , σ 2 97 ] T , with the resulting estimates beinĝ θ = [12.46, 1.54, 0.53, 0.88, 0.52, 1.01] T .
Using the estimated covariance parameters to estimate the covariance matrix of the observed data, the generalized least squares estimates of the trend parameter vector β isβ gls = [−332.02, −2.17, 1.73, 0.0039] T . Examining all of the available total monthly precipitation measurements for the month of October in the state of Oregon from 1895-1997, 250 mm of total rain corresponded to roughly the 96th percentile and was chosen as the exceedance threshold for which exceedance locations would be identified in October 1998. Accordingly, on the transformed scale our goal is to identify at a confidence level of 0.90 the exceedance region E Z √ 250 + where the transformed total monthly precipitation measurements exceed u = √ 250 mm. The next step in our analysis was creating a discrete grid in the region of interest D, where D is the state of Oregon. For computational simplicity, a 100 × 100 rectangular grid was superimposed over D, and the pixels with center points contained in D were used as the discrete grid in subsequent analysis. Next, the procedure outlined in Section 2.3 was used to construct both S √ 250 + and S c √ 250 − at a confidence level of 0.90 using 2000 realizations of (7); the resulting confidence regions are shown in Figure 5. The yellow coloring in Figure 5 indicates the regions where the predicted value (using universal kriging) is greater than the designated threshold, the blue coloring indicates the 90 percent confidence region S √ 250 + , while orange coloring indicates the 90 percent confidence region S c √ 250 − . Naturally, the confidence region S √ 250 + contains the predicted exceedance region, while the predicted exceedance region contains S c √ 250 − . The area of Oregon predicted to have total monthly precipitation greater than 250 mm in October of 1998 is mainly limited to the western area of the state closer to the coast. The associated 90 percent confidence region S √ 250 + where 250 mm of rain could fall in October of 1998 is also found in the western portion of the state, though two small regions in the northeast corner of Oregon were also considered as possible candidates for this event. Only a small region (shown in orange) in the northwest part of the state could confidently be predicted to receive at least 250 mm of total monthly precipitation during this time period. On the other hand, the regions without any shading are the regions where we can be confident that total monthly precipitation will not exceed 250 mm, which comprises most of the eastern part of the state.
5. Case study 2: Regional climate projections. We continue to demonstrate the proposed methodology by using it to explore similarities and differences between projections of climate models from the North American Regional Climate Change Assessment Program [NARCCAP; Mearns et al. (2007Mearns et al. ( ), updated (2012]. NARCCAP is an international, multi-disciplinary program exploring "separate and combined uncertainties in regional projections of future climate change resulting from the use of multiple atmosphereocean general circulation models (AOGCMs) to drive multiple regional climate models (RCMs)" as well as "to provide the climate impacts and adaptation community with high-resolution regional climate change scenarios that can be used for studies of the societal impacts of climate change and possible adaptation strategies" [Mearns et al. (2009);]. Data produced by the program are available for numerous combinations of AOGCMS and RCMs, allowing researchers to investigate and study how various models interact, compare, and contrast with each other. Subsequent analysis will focus on combinations of two AOGCMs (CCSM and CGCM3) and four RCMs (CRCM, MM5I, RCM3, and WRFG) with a total of six models considered (CRCM/CCSM, CRCM/CGCM3, MM5I/CCSM, RCM3/CGCM3, WRFG/CCSM, and WRFG/CGCM3). We consider seasonal averages (Dec-Feb, Mar-May, Jun-Aug, Sep-Nov) of temperature (degrees Celsius) for years between 1971 and 2000 (note that the months are consecutive so that the December of the previous year is included in the average for the current year) on a 50 km grid covering Canada, the United States, and the northern part of Mexico. Projections from each model are also available for years between 2041 and 2070. Potential predictors to capture large-scale spatial trends include longitude, latitude, and elevation. The four main goals of this case study are to: 1. Compare and contrast the different AOGCM/RCM models to determine whether they portray the same type of behavior.
2. Explore the impact of season on climate predictions. 3. Study the effect of changing the threshold level on the associated exceedance region.
4. Assess where temperature increase is likely to occur based on these models. The size of the NARCCAP spatial grid and the wide variety of climatological conditions across the domain necessitated choices for how the data were analyzed. In order to capture the potential temperature increase between current and future model runs, future and current model runs were paired by differences of 70 years (e.g., the average winter temperature for 1971 was paired with the average winter temperature of 2041) and the difference between the future and current temperature values were taken. Depending on how the various RCMs were run, this left between 29 to 30 years of differences for each season. To simplify this part of the analysis, the data were separated into 9 year groupings (1971-1979/2041-2049, 1980-1988/2050-2058, and 1989-1997/2059-2067), and the average temperature difference for each season was calculated. Further, several smaller and climatologically consistent regions of the domain were analyzed individually to proceed with analysis. Examples of these smaller regions, Boreal, Central, Desert, East, MtWest, and South, are overlaid on a plot of North America in Figure 6. Using the averages from the three groupings of temperature, statistical inference for each region proceeded by predicting the average difference between seasonal temperature between the years 1998-2007 and 2068-2077, and constructing confidence regions for the exceedance regions of certain temperature thresholds. For the purposes of this paper, we will focus on results from the Boreal, South, and East regions during the winter and summer seasons.
Each climate model was analyzed in the following manner. The mean structure of the responses Z(s, t) was assumed to be constant across nineyear groupings so that µ(s, t) = x(s, t) T β, where the covariates vector x(s, t) = [1, longitude, latitude, elevation] T . The covariance function of the observed data was modeled as being isotropic and fully separable with respect to space and time, using an exponential covariance structure for the spatial covariance and an AR(1) structure for the temporal covariance. Measurement error was assumed to have the same variance for all of the nine-year groupings. The covariance function of the observed data may be written as where h = s 1 − s 2 is the Euclidean distance between the two data locations and u = |t 1 − t 2 | is the time lag between t 1 and t 2 . Restricted maximum likelihood estimation was used to estimate the covariance parameters θ = [σ 2 , φ, ρ, σ 2 ε ] T assuming a multivariate Gaussian distribution for the observed responses.
The next step in our analysis was creating a discrete grid in the region of interest D, where D is the individual region in question. Over each domain, a 100 × 100 grid of regular pixels was overlaid, and then the center points of the pixels within the convex hull bounding the domain were retained for use as prediction locations. Following the procedure outlined in Section 2.3, confidence regions for the exceedance regions of the temperature change were constructed for levels u = 1, 2, and 3 • C using 10,000 realizations of the conditional random field in (7). For each region and exceedance level, both the region S u + containing the true exceedance region E Z u + and the region S c u − for which all points in the region should be part of the exceedance region were obtained. The region S u + is shown in blue and the region S c u − is shown in orange. Note that because of overlap, all orange locations are covering blue.
We begin by examining the effects of season on temperature predictions and the associated confidence regions for the Boreal region during the winter and summer seasons. The confidence regions associated with temperature change of at least 1 • C in winter are shown in Figure 7(a), and for summer in (b).
Comparing the graphs in Figure 7(a) and (b), the pattern of climate change is not necessarily consistent between winter and summer. Specifically, while the orange area (the area that we are confident will be part of the true exceedance region) for the winter data in Figure 7(a) is consistently in the northeastern area of the Boreal region, we do not see a similar pattern for the summer data in Figure 7(b). The areas that we can confidently identify as experiencing temperature increase in one season may not experience the same change in a different season. The confidence region S 1 + (blue + orange) for the true exceedance set E Z 1 + during both winter and summer comprises the entire Boreal region. Consequently, any part of the Boreal region could experience a temperature increase of at least 1 • C in both winter and summer. For the winter data, the orange region consistently makes up a large percentage of the northeastern Boreal region. Based on the agreement between these climate models, there appears to be high confidence that the northeastern part of the Boreal region will experience a temperature increase of at least 1 • C when comparing the average winter temperature between the years 1989-1997 and 2059-2067. On the other hand, while several orange regions appear for the summer temperature data in Figure 7(b), there does not appear to be consistency between the various climate models. The CRCM/CCSM model results (upper left) in Figure 7 Fig. 7. Confidence regions for locations in the Boreal region where temperature increase will be at least 1 • C for various climate models. The results for winter are shown in (a) and for summer in (b). Blue shading indicates where the temperature increase may occur. Orange shading indicates where the temperature increase is likely to occur. a majority of the region will experience temperature change of a least 1 • C during the specified time period, while the WRFG/CGCM3 model does not indicate that any of the region will confidently experience this same change. While individual climate models may deem a certain region as being likely to experience a certain level of temperature increase, other climate models may paint a different picture.
We next consider temperature change in the South region, which comprises much of the southeastern part of the United States, and focus on how different RCM/AOGCM climate model combinations can make fairly different predictions. The results for temperature change of at least 3 • C during winter and summer are shown in Figure 8(a) and (b), respectively.
In the South region, the potential temperature increase portrayed by the climate models is fairly consistent from winter to summer. Additionally, the results for each of the climate models appear relatively similar with the exception of the WRFG/CGCM3 combination. Specifically, this combination predicts temperature increase in the South region as being smaller than each of the other climate models. This difference cannot be explained by the use of the CGCM3 AOGCM or WRFG RCM alone, since other model combinations having these model components do not behave in the same manner. This specific RCM/AOGCM combination interacts in a way that yields fairly different results than the other combinations considered. We note that the nonshaded regions are the regions confidently identified as not having a temperature increase of 3 • C or more.
Last, we consider the effect of changing the threshold level on the resulting confidence region for the exceedance region by looking at results for the East region along the Eastern seaboard of the United States and Canada. Results for winter temperature increases of 1, 2, and 3 • C are shown in Figure 9(a), (b), and (c), respectively. For all three thresholds, only the northeast part of the East region confidently has temperature increases exceeding the threshold in question (the areas shown in orange). However, as the temperature threshold increases, the size of the orange area decreases. This behavior is sensible since the more extreme a threshold is, the less likely it is that a response will exceed that threshold. Consequently, the orange areas where we are confident a response exceeds a threshold become smaller as the threshold increases. Based on the observation that the orange and blue shaded regions make up nearly the entire East region for all three thresholds, nearly any part of the East region of North America may experience an average temperature increase of at least 3 • C when comparing the winter temperature between the years 1989-1997 and 2059-2067. On the other hand, since the orange area is mostly limited to the northeastern part of the East region for a threshold of 1 • C, this is the only region the climate models can confidently identify as experiencing a temperature increase of at least 1 • C. On the other hand, as the temperature threshold is increased to 3 • C, some areas have no shading. These are the regions that can be confidently identified as not experiencing a temperature increase of 3 • C or more. However, the nonshaded regions are not consistent between the various AOGCM/RCM combinations.
(a) (b) Fig. 8. Confidence regions for locations in the South region where temperature increase will be at least 3 • C for various climate models. The results for winter are shown in (a) and for summer in (b). Blue shading indicates where the temperature increase may occur. Areas without shading are regions where the temperature increase is unlikely to occur. (c) Fig. 9. Confidence regions for locations in the East region where winter temperature increase will be at least 1 • C, 2 • C, and 3 • C are shown in (a), (b), and (c), respectively. Blue shading indicates where the temperature increase may occur. Orange shading indicates where the temperature increase is likely to occur. Areas without shading are regions where the temperature increase is unlikely to occur.
We point out that the large blue regions in the preceding discussion may not be precise indicators of the size of the true exceedance region. The large size of the blue regions may simply be an indication of high uncertainty, since the size of the confidence regions S u + will increase with elevated uncertainty.
6. Discussion. We have presented an approach for constructing confidence regions containing the entire exceedance region of a random field. Additionally, by considering the inverse problem, one may also obtain a region where all locations in the region will confidently be part of the exceedance region of the random field. This allows researchers to compare regions where it is possible for an extreme event to occur to the regions where it is likely that an extreme event will occur. The size of the confidence region naturally depends on factors such as the number of observed responses, the spatial and temporal dependence, and the magnitude of measurement error variance. Simulation experiments in Section 3.2 indicate that the procedure produces confidence regions having the appropriate coverage properties. Though all of these experiments used stationary latent processes for simplicity, stationarity is not required and simulation experiments for nonstationary processes have produced similar results. As further explored in Section 3.3, the shape of a confidence region naturally patterns itself after the underlying mean structure.
Using the proposed procedure, we were able to make inference about the regions of Oregon receiving 250 mm of precipitation in October of 1998. It also allowed us to compare climate models and explore future climate change for several combinations of RCMs and AOGCMs using models from NARCCAP. Though the statistical models used were basic, the results from these models supported the view that temperature increases of several degrees are possible for large parts of North America, and that certain areas seem likely to experience temperature change of several degrees. The analysis also revealed that there can be somewhat large discrepancies between climate models (the South region being a clear example).
It should be mentioned that the assumption that the spatio-temporal covariance functions of the statistical models were separable may be unnecessarily restrictive. Nonseparable and/or nonstationary spatio-temporal covariance models such as the ones proposed by Gneiting (2002) or Fuentes, Chen and Davis (2008) are possible alternatives. Additionally, the confidence regions were often quite large. This often was a consequence of the fact that most of the region in question was predicted to be greater than the exceedance threshold, but was sometimes a result of high predictive uncertainty in the areas. This suggests that similar tools with less stringent error criteria [e.g., Sun et al. (2012)] would be useful for climate model exploration. Related to this is the fact that more data brings more information and, consequently, may reduce the size of the confidence regions. Due to the size and complexity of the NARCCAP model output, specific regions of North America were analyzed in isolation. If the proposed procedure were extended to incorporate reduced rank modeling procedures such as fixed-rank kriging [Cressie and Johannesson (2008)] or fixed-rank filtering [Cressie, Shi and Kang (2010)], this limitation could be ameliorated, and this is the subject of ongoing research efforts. Additionally, it was pointed out in Section 2.4.4 that the confidence levels of the proposed methodology assumed that the covariance function of the random process was known, which is rarely the case. A Bayesian approach to this problem that naturally incorporates the uncertainty of the covariance function is under investigation.