Hierarchical Transformed Scale Mixtures for Flexible Modeling of Spatial Extremes on Datasets with Many Locations

Introduction

Modeling the dependence structure in the extremes of spatial processes is of great consequence for risk analysis of extreme events. In this paper, we make a slight alteration to the flexible class of randomly scaled transformed Gaussian process models to sidestep computational bottlenecks normally encountered in the likelihood. Our modification enables high-dimensional inference, while preserving submodels that transition smoothly between extremal dependence classes.

Generally, the probability that two spatially-indexed random variables exceed a high level simultaneously varies by their separation distance, and the particular way in which this joint probability decays must be well-represented in models if one hopes to accurately to assess risks posed by spatial extremal phenomena. Good estimation of how the dependence changes both with distance in space and as one moves farther into the joint tail will enable us to accurately calculate exceedance probabilities of areal quantities, predict at un-observed locations, and, secondarily, get a more realistic picture of marginal quantities.

In classical spatial modeling, Gaussian processes have been widely used due to their mathematical simplicity and tractability for larger datasets. However, the Gaussian density function is very light-tailed, and thus has the potential to underestimate probabilities associated with extreme events; furthermore, Gaussian models stipulate that the dependence among rare events at distinct locations will always diminish such that the probability of observing an extreme at one location, conditional upon an extreme at another location, is zero in the limit. This property is termed asymptotic independence, but models that only exhibit this phenomenon may be too inflexible for applications where the true tail dependence structure is uncertain.

Max-stable processes form an important class of models that exhibit the alternative scenario of asymptotic dependence. They are the natural extension of classical univariate extreme value theory to infinite dimensional settings, and therefore can provide an asymptotically-justified modeling framework for datasets consisting of block-maxima. Counterparts of max-stable processes suitable for threshold exceedances are called generalized Pareto processes (Ferreira and de Haan, 2014; Thibaud and Opitz, 2015). These processes are also asymptotically dependent, but possess the advantage that they bypass many of the computational difficulties of max-stable processes.

Despite the theoretical appeal of limiting max-stable and generalized Pareto processes, there are two main drawbacks to these models: (i) the assumption of asymptotic dependence may be incorrect and (ii) even if the data are asymptotically dependent, they will often not appear to follow such limiting models at sub-asymptotic levels. Both max-stable and generalized Pareto process dependence structures exhibit stability properties, meaning that their dependence structures are invariant to the operations of taking maxima and conditioning upon exceedances of higher thresholds, respectively. If the true data generating process exhibits weakening dependence in the un-observed region of the tail, inference drawn under these models about the far joint tail will over-estimate risk, sometimes substantially.

On account of the limitations of limiting models, it is desirable to find a family of spatial models that can transition between asymptotic dependence and asymptotic independence. In particular, we will be examining a class of marginally transformed Gaussian scale mixture models, which includes those recently proposed by Huser et al. (2017) and Huser and Wadsworth (2019). These models are of great interest due to their appealing theoretical properties and ability to flexibly capture both types of extremal dependence structure. Unfortunately, inference for these models is not feasible for large numbers of observed sites, as calculation of the censored likelihood, the preferred method for fitting joint tail models, entails integration over high-dimensional multivariate Gaussian distribution functions. To increase scalability, we propose an adaptation of the model by adding an independent measurement error term to each component. By adding this nugget effect, the new model circumvents the lengthy computation of the multivariate normal distribution function. Also it can elegantly avoid the integral of the process below the censoring threshold by considering the uncensored process as latent and drawing from it using Gibbs sampling, allowing for truly high-dimensional inference. Furthermore, we show that the modified models retain all the significant asymptotic properties of the original smooth models despite the presence of the measurement errors, which lays a solid theoretical foundation for correctly capturing the sub-asymptotic dependence behavior.

The article is organized as follows. Section 2 provides a brief literature review on the measures of extremal dependence and hybrid spatial extreme models, and further explains the intractability of the existing censored likelihood approaches for inference on these hybrid models. Section 3 describes our new model that alleviates the computational problems, and studies its extremal dependence properties. Section 4 includes a marginal transformation in the hierarchical model and details the inference using Gibbs sampling. Section 5 presents a simulation study that validates the methodology. We apply our model to a dataset of the Fosberg Fire Weather Index (FFWI) on the Great Plains in Section 6. Section 7 concludes with some discussion. Appendix A provides proofs of all the theoretical results. Appendix B includes supplementary diagnostics for the data application.

Spatial Dependence for Extremes

For a stochastic process $\{X(\boldsymbol{s}):\boldsymbol{s}\in\mathcal{S}\}$, we write $X_{j}=X(\boldsymbol{s}_{j})$ and so forth for simplicity, where $\boldsymbol{s}_{j}$ denotes the $j$th spatial location. It is useful to summarize the extremal dependence implied by the observed process concisely.

We restrict the scope to the bivariate case, focusing on stationary and isotropic random fields. One example of a bivariate dependence measure is the upper tail dependence coefficient:

where $X_{j}\sim F_{j}$, $X_{k}\sim F_{k}$, and $h=\|\boldsymbol{s}_{j}-\boldsymbol{s}_{k}\|$. Joe (1993) defined the upper tail dependence parameter as the limit $\chi(h)=\lim_{u\rightarrow 1}\chi_{u}(h)$. Asymptotic dependence is attained if and only if $\chi(h)>0$, while $\chi(h)=0$ defines asymptotic independence.

For max-stable processes, $\chi_{u}(h)=2-V(1,1)+O(1-u),\;u\to 1$, where $V(\cdot,\cdot)=\log F_{jk}(\cdot,\cdot)$. Max-stable distributions can be associated to a generalized Pareto counterpart, for which $\chi_{u}(h)\equiv\chi(h)=2-V(1,1)$ for all $u$ above a certain level (Rootzén et al., 2018). The fact that $\chi_{u}(h)$ does not depend on $u$ is a manifestation of the threshold-stability of generalized Pareto processes. In practice, empirical estimates of (1) from data tend to show $\chi_{u}(h)$ decreasing with both $h$ and $u$, meaning that realistic models should also have this property.

When $\chi(h)=0$, i.e., the case of asymptotic independence, further detail about the behavior of $\chi_{u}(h)$ is obtained by exploiting the joint tail assumption of Ledford and Tawn (1996):

where $\mathcal{L}$ is slowly varying at zero, that is, $\lim_{t\rightarrow 0}\mathcal{L}(tx)/\mathcal{L}(t)=1$ for any $x>0$, and $\eta_{X}(h)\in(0,1]$ is the coefficient of tail dependence of the process $X$. The pair of variables $(X_{j},X_{k})$ are asymptotically dependent when $\eta_{X}(h)=1$ and $\mathcal{L}(\cdot)\nrightarrow 0$. The remaining cases are all asymptotically independent, and the value of $\eta_{X}(h)$ characterizes the strength of extremal dependence in the upper joint tail. In the case of a Gaussian process, $\eta_{X}(h)=\{1+\rho(h)\}/2$, where $\rho(h)$ is the correlation at lag $h$. The variables are called positively associated when $\eta_{X}(h)>1/2$ and negatively associated when $\eta_{X}(h)<1/2$. Near independence corresponds to $\eta_{X}(h)=1/2$.

Gaussian processes are asymptotically independent for all correlations $\rho(h)\neq 1$. They might be considered candidates for modeling the joint tail of asymptotically independent phenomena, but as there is no theory to specifically recommend Gaussian processes in this scenario, it is desirable to consider other models as well. As an alternative, Opitz (2016) captures spatial dependence in asymptotically independent processes by construction of Laplace random fields, defined as mixtures of Gaussian processes with a random variance that is exponentially distributed. Wadsworth and Tawn (2012) proposed the class of inverted max-stable processes, for which the tail decay is specified fully by $\eta_{X}(h)$, although inference is computationally challenging.

In real datasets, it is difficult to conclude definitively whether data exhibit asymptotic independence or asymptotic dependence, and incorrectly assuming an asymptotically independent model can lead to equally severe problems with bias as incorrectly assuming an asymptotically dependent model. Because of this, a recent focus in the literature has been on models that can encompass both scenarios.

Wadsworth and Tawn (2012) were the first to introduce hybrid models that combine max-stable and inverted max-stable processes so that asymptotic dependence prevails at short distances, and asymptotic independence at long distances. However, inference for this model is difficult because there are a fairly large number of parameters involved, and the transition between the dependence classes takes place at the boundary of the parameter space.

Recently, several Gaussian scale mixture models were proposed to allow more flexible transitions between dependence classes. Through multiplying an asymptotically independent Gaussian process by a random effect that governs the extremal dependence, these models can be described by a small number of parameters and have non-trivial asymptotically independent and asymptotically dependent submodels. More precisely, suppose $\{Z(\boldsymbol{s}),\boldsymbol{s}\in\mathcal{S}\}$ is a standard isotropic and stationary Gaussian process with covariance function $C_{{\boldsymbol{\theta}}_{C}}(h)$ indexed by a parameter vector ${\boldsymbol{\theta}}_{C}$, where $h$ is the length of the separation vector, so that $\Sigma_{{\boldsymbol{\theta}}_{C}}$ is the covariance matrix of associated finite-dimensional distributions. The class of Gaussian scale mixture models can be constructed as

where $g(\cdot)$ is a link function, and $R>0$ is a random scaling factor, from distribution $F_{R}$ indexed by ${\boldsymbol{\theta}}_{R}$, that can be interpreted as a constant random process over spatial domain $\mathcal{S}$ with perfect dependence. Impacting simultaneously the whole domain $\mathcal{S}$, heavier tailed $R$ induces asymptotic dependence in $X^{*}$, whereas lighter tailed $R$ induces asymptotic independence. Engelke et al. (2019) provide a fuller description of how extremal dependence of $X^{*}$ relates to the relative marginal tail heaviness of $R$ and $g(Z)$.

Morris et al. (2017a) uses a space-time model based on skew-$t$ process, where $g(\cdot)$ is a identity function, $R^{2}\sim\text{IG}(a/2,b/2)$ is an inverse gamma random variable, and $C_{{\boldsymbol{\theta}}_{C}}$ is a Matérn covariance function. On top of the mixture, they added covariate effects and a skew term. Since the inverse gamma distribution is heavy tailed, the skew-$t$ process is asymptotically dependent for $a<\infty$. Asymptotic independence is achieved only when $a\rightarrow\infty$.

Huser et al. (2017) also used an identity link function, but placed few assumptions on the random scale, and provided more general results on the joint tail decay rates of the mixture processes. They showed that a wide class of Weibull-like tail decay in $R$ yields asymptotic independence, while a Pareto-like tail that is regularly varying at infinity gives asymptotic dependence. They also proposed a parametric model that bridges the two asymptotic regimes and provides a simple transition, in which $R$ is a two-parameter distribution

where $\gamma>0$, and the support is $[1,\infty)$. Since $(r^{\beta}-1)/\beta$ converges to $\log r$ as $\beta$ approaches 0, (4) forms a continuous parametric family on $\beta$. When $\beta>0$, (4) constitutes a class of Weibull-type distributions and thus assures asymptotic independence. When $\beta=0$, the variable $R$ is Pareto distributed and thus gives asymptotic dependence. This shows that the model provides greater flexibility and can transition from asymptotic dependence to independence via adjusting the value of $\beta$.

However, the previous two Gaussian scale mixture models both make the transition between the dependence classes at the limit or the boundary of the parameter space. They are also inflexible in their representation of asymptotic dependence structures because there is dominating preference over one dependence class. It may be more desirable to find a model for which the transition takes place in the interior of the parameter space so one we can quantify the uncertainty about the dependence class in a simpler manner. To overcome this, Huser and Wadsworth (2019) proposed a marginally transformed Gaussian scale mixture model, where $g(\cdot)$ transforms a standard Gaussian variable to standard Pareto, and $R$ itself is Pareto distributed:

Here the type of asymptotic dependence is determined by the value of $\delta$. When $\delta\leq 1/2$, $R$ is lighter tailed or equivalent to standard Pareto, which induces asymptotic independence; when $\delta>1/2$, the converse is true, which induces asymptotic dependence. Specifically, the upper tail dependence parameter $\chi_{X^{*}}=\frac{2\delta-1}{\delta}E\left[\min\{g(Z_{i}),g(Z_{k})\}^{(1-\delta)/\delta}\right]$ when $\delta>1/2$ and 0 otherwise, while the coefficient of tail dependence is

where $\eta_{Z}$ is the coefficient of tail dependence for $(Z_{i},Z_{k})$ (Huser and Wadsworth, 2019).

The model in (5) provides a smooth transition through asymptotically independent and asymptotically dependent submodels. It has many appealing asymptotic properties. However, inference for models of the form (3) is typically made via censored likelihood. This requires computing an integral where the integrand contains the Gaussian distribution function in $|\mathcal{C}|$ dimensions, where $|\mathcal{C}|$ is the number of components below a designated high threshold. Such integrals are computationally prohibitive for even moderately-sized datasets. In Section 3 we introduce a slight alteration to this model to make it tractable while preserving all the desired asymptotic results.

2 The Censored Likelihood

In multivariate and spatial extremes, the preferred approach to fitting the dependence structure is using a censored likelihood, which prevents observations from the bulk of the distribution from affecting the estimation of the extremal dependence structure. It provides a reasonable compromise between bias and variance compared to alternative approaches, although different censoring schemes have been adopted (Thibaud and Opitz, 2015; Huser et al., 2016).

For a process of the form (3) observed at $D$ spatial locations $\boldsymbol{s}_{1},\cdots,\boldsymbol{s}_{D}\in\mathcal{S}$, we obtain the distribution function by conditioning on $R$ as

where $r^{*}=\min(x^{*}_{1},\cdots,x^{*}_{D})$, and $\Phi_{D}$ denotes the $D$-variate Gaussian distribution with zero mean and covariance matrix $\mathbf{\Sigma}_{{\boldsymbol{\theta}}_{C}}$.

Let $\mathcal{C}\subseteq\{1,\ldots,D\}$ be the set of locations with censored observations—that is, the set of locations where the components are below a high threshold; let $\mathcal{U}$ be the set of locations with uncensored observations. For any index set $A,B\subset\{1,\ldots,D\}$, denote $\boldsymbol{x}_{A}=\{\boldsymbol{x}_{i}:i\in A\}$, $\mathbf{\Sigma}_{A;B}$ as the matrix $\mathbf{\Sigma}$ restricted to the rows in $A$ and the columns in $B$, and let $\mathbf{\Sigma}_{A|B}$ be the Schur complement of $B$ in $\mathbf{\Sigma}_{A;B}$. The likelihood is obtained via taking partial derivatives of (6) with respect to $\mathcal{U}$:

Although only one-dimensional integral appears in (7), the integrand includes a $|\mathcal{C}|$-dimensional Gaussian distribution function. When approximating the integral using standard quadrature or Monte Carlo methods, one needs to compute $\Phi_{|\mathcal{C}|}$ for each sample point taken on $(1,r^{*})$. This is only feasible when the number of locations $D$ is moderate. Additionally, this calculation will have to be repeated for each time replicate.

To avoid the integrating the process below the threshold, one could instead think of $X^{*}(\boldsymbol{s})$ as latent and draw from it using Monte Carlo methods. Consequently there is no need to compute the awkward likelihood (7). However, to update the Markov chain each time, it is now necessary to draw $\boldsymbol{x}^{*}_{\mathcal{C}}$ from a high-dimensional truncated distribution, which might again be computationally intensive.

Therefore, we propose to make a slight adjustment to the model in (3). Our new model is markedly more amenable to higher-dimensional inference, yet it keeps hold of the joint tail decay rates attained in the original model (e.g. Huser and Wadsworth, 2019; Huser et al., 2017). Equivalently, our new model has non-trivial asymptotically dependent and asymptotically independent submodels with the transition taking place in the interior of the parameter space in the case of our modified version of (5).

Model

We alter the models in Section 2.1 by adding an independent measurement error term to each component,

where $\epsilon_{i}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\tau^{2}),\;i=1,\ldots,D$, and distribution of $R$ and the link function remain the same. That is, we add a simple nugget effect to the smooth process $X^{*}(\boldsymbol{s}_{i})$. When drawing the latent processes below the threshold, we can condition on the smooth process and simply update the noisy one. Because these error terms are independent of each other, there is only a univariate integral involved in the full conditional likelihood. Also, when we update the smooth process $X^{*}(\boldsymbol{s})$ given the noisy process $X(\boldsymbol{s})$, no truncation or censoring is present and it is much easier to sample from the corresponding likelihood. Section 4 contains more details on the Markov Chain Monte Carlo (MCMC) updating scheme, where we show how this small alteration can hugely facilitate inference.

2 Dependence Properties

We begin with the model (5) from Huser and Wadsworth (2019), modified as in (8). Recall that $g(Z(\boldsymbol{s}))$ is a stationary process with standard Pareto margins possessing asymptotic independence; i.e., $P(g(Z(\boldsymbol{s}))>x)=x^{-1}$ and

where $\mathcal{L}_{Z}(x)$ is slowly varying at infinity, and $\eta_{Z}(h)=(1+\rho(h))/2<1$ for the Gaussian correlation $\rho(h)<1$.

Figure 1 illustrates the estimated coefficient of tail dependence $\eta_{X}$ as a function of $\delta$ for $\eta_{Z}=0.1,\ldots,0.9$. For each combination of $\delta$ and $\eta_{Z}$, we generate 5,000,000 replicates from model (5) (i.e. $\tau^{2}=0$) and model (8) with $\tau^{2}=1$ respectively. For each replicate, we sample $(Z_{i},Z_{k})$ from a Gaussian copula with correlation $2\eta_{Z}-1$. We then numerically approximate the joint survival probability in (2) to obtain an estimate of $\eta_{X}$. The left panel of Figure 1 clearly shows that the smooth transition from asymptotic independence to asymptotic dependence takes place around $\delta=1/2$, confirming the results from Huser and Wadsworth (2019) with a reasonable bias; the right panel shows that adding a measurement error has little effect on the tail dependence because $\eta_{X}$ exhibits similar behavior. This result invites investigation of whether the flexible asymptotic properties in Huser and Wadsworth (2019) are preserved in the altered model. In the following, we generalize the problem from the specific model of Huser and Wadsworth (2019) for the process $X^{*}$ in (8), to any $X^{*}$ with a wide class of marginal tail behaviors.

The impact of the additive Gaussian nugget effect on the extremal dependence of $X$ depends upon the marginal tail heaviness of $X^{*}$: roughly, the heavier the tail of $X^{*}$, the less the impact of the noise. Since we take a copula-like approach and employ $X$ as the spatial dependence model, we assume its margins, and those of $X^{*}$, are identical over space.

Regularly varying tails are defined through the survival function being regularly varying at infinity, i.e., $\mbox{P}(X^{*}>x)\in\mbox{RV}_{-\alpha}$, $\alpha>0$. Weibull-like tails are defined through the survival function

where $u\in\mbox{RV}_{\kappa}$, and $\alpha$ is termed the Weibull index. We also assume that $X^{*}$ has a density satisfying $f_{X^{*}}(x)\sim v(x)\exp(-\theta x^{\alpha})$, with $v(x)=u(x)(\theta\alpha x^{\alpha-1})$. The main results are now summarized in Proposition 3.1.

If $X^{*}$ has a regularly varying tail, or Weibull-like with Weibull index $\alpha<1$, then $\chi_{X}=\chi_{X^{*}}$ and $\eta_{X}=\eta_{X^{*}}$.

If $X^{*}$ has a Weibull-like tail with Weibull index $\alpha=1$ then

and $\eta_{X}=\eta_{X^{*}}$. Note if $\chi_{X^{*}}=0$ then so is $\chi_{X}$.

If $X^{*}$ has a Weibull-like tail with Weibull index $\alpha>1$ then:

If $\alpha\in(1,2)$, then $\eta_{X}=\eta_{X^{*}}$

If $\alpha=2$, then an interval can be given for $\eta_{X}$ (see Expression (27)).

The proof of Proposition 3.1 is given in Appendix A.

For the process of Huser and Wadsworth (2019), described in equation (5). In this case, $X^{*}$ always has a regularly varying tail, so Proposition 3.1 part 1 gives $\eta_{X}(h)=\eta_{X^{*}}(h)$, and $\chi_{X}(h)=\chi_{X^{*}}(h)$, with $\eta_{X^{*}}(h)$ and $\chi_{X^{*}}(h)$ given in Section 2.1. This means the flexible asymptotic properties in Huser and Wadsworth (2019) are preserved in the altered model.

Another popular model for spatial data is the $t$-process (Røislien and Omre, 2006), which is a Gaussian scale mixture for which the mixing variable $R^{2}$ follows an inverse gamma distribution. If $X^{*}$ is a $t$-process, then it is asymptotically dependent with a regularly varying tail, so $\eta_{X}(h)=\eta_{X^{*}}(h)=1$ and $\chi_{X}(h)=\chi_{X^{*}}(h)>0$ by Proposition 3.1. Similarly, the skew-$t$ process (Padoan, 2011; Morris et al., 2017b) is regularly varying and asymptotically dependent, so the same conclusions apply.

For the Gaussian scale mixture of Huser et al. (2017), $X^{*}$ either has regularly varying or Weibull-like tails depending on the distribution of the scaling variable $R$. In particular when $\beta>0$ in distribution (4), $X^{*}$ has a Weibull-like tail with Weibull index $\alpha=2\beta/(\beta+2)<2$. As such, parts 1, 2 and 3a of Proposition 3.1 are relevant. In all cases $\eta_{X^{*}}(h)=\eta_{X}(h)=\{(1+\rho(h))/2\}^{\beta/(\beta+2)}$, and $\chi_{X}(h)=\chi_{X^{*}}(h)=0$ for $\beta>0$. When $\beta=0$ then $X^{*}$ has a regularly varying tail, and is asymptotically dependent with $\eta_{X}(h)=\eta_{X^{*}}(h)=1$, and

where $T_{\nu}$ is the cdf of the Student-$t$ distribution with $\nu$ degrees of freedom.

We note that the process $X^{*}$ in equation (3) is constructed only for its dependence properties, and there is no “natural” scale on which to express it. For example, considering the process of Huser and Wadsworth (2019) we could also write

with $E\,|\,\delta\sim\mbox{Exp}(\delta/(1-\delta))$, $V(s)=\log g(Z(s))$, which has the same dependence structure as defined in (3) and (5), since it is obtained through a monotonic marginal transformation. Taking $X^{*}$ from (11), Proposition 3.1 part 2 gives $\eta_{X}(h)=\eta_{X^{*}}(h)$. Asymptotic dependence of $X^{*}$ implies asymptotic dependence of $X$, but only bounds on $\chi_{X}(h)$ are available.

In practice, if choosing a marginal scale for $X^{*}$, there may be a trade-off between theoretical desires and computational practicality. Supposing that we wish $X$ to inherit the properties of $X^{*}$, a heavy-tailed choice is best.

Bayesian Inference

We define a Bayesian hierarchical model based on the process (8) defined in Section 3.1 and use a MCMC algorithm to fit to the data. For the reasons outlined in Section 2.2, we assume data are censored below a high threshold $u$. In (8), the mixing parameter ${\boldsymbol{\theta}}_{R}$ controls both joint and marginal behavior of the response $X$, which we would prefer to separate. Therefore, motivated by the theory of univariate extremes, we first assume our observations above the same high threshold $u$ are generalized Pareto distributed, and we include a marginal transformation in the hierarchical model.

Let $\{Y(\boldsymbol{s}):\boldsymbol{s}\in\mathcal{S}\}$ denote the observed process. We define a marginal transformation $T(y)$ as follows:

where $F_{X|{\boldsymbol{\theta}}_{R},\tau^{2}}$ is the marginal distribution function for process (8), and

where $p=P(Y(\boldsymbol{s}_{i})\leq u)$, $u$ is a high threshold, and $F_{GPD|u,\sigma,\xi}(y)=1-[1+\xi(y-u)/\sigma]^{-1/\xi}$, with support $\{y\geq u:1+\xi(y-u)/\sigma\geq 0\}$. Conditioning on the smooth process $X^{*}(\boldsymbol{s}_{i})=x^{*}(\boldsymbol{s}_{i})$, which was not truncated, the censored likelihood for an observation $y(\boldsymbol{s}_{i})$ can be derived as

where $\boldsymbol{\theta}_{GPD}=(\sigma,\xi)$. Note that there are only univariate calculations required in (14), compared to (7) for which we have to estimate the $|\mathcal{C}|$-dimensional Gaussian distribution functions. In addition, since $Y_{i}$ and $Y_{k}$ are independent conditioning on the smooth process ($i\neq k$), the joint likelihood of the whole vector $\boldsymbol{y}$ is simply $f_{\boldsymbol{Y}}(\boldsymbol{y})=\prod_{i=1}^{D}\varphi(y(\boldsymbol{s}_{i})\,|\,\boldsymbol{X}^{*},\tau^{2},\boldsymbol{\theta}_{R},\boldsymbol{\theta}_{GPD},p)$. Likelihoods for independent time replicates are simply multiplied together, and the proportion of censored observations $p$ can be treated as a known parameter or an unknown parameter that enters the hierarchical model. See Appendix A.2 for a complete statement of the hierarchical model. The priors for the model parameters are

where halfCauchy(1) refers to the positively truncated standard Cauchy distribution. We implement this methodology for two different Gaussian scale mixture models: that of Huser and Wadsworth (2019), and that of Huser et al. (2017). The priors for $\boldsymbol{\theta}_{R}$ are $\delta\sim U(0,1)$ for the former, and $\beta\sim\text{halfCauchy}(1)$ for the latter. The prior for ${\boldsymbol{\theta}}_{C}$ depends on the choice of the covariance function. In our implementations, we adopt the Matérn covariance function with ${\boldsymbol{\theta}}_{C}=(\rho,\nu)$, where $\rho$ is the range parameter, and $\nu$ is the smoothness parameter. The prior is then set to subject to two independent half Cauchy distributions with scale parameter 1.

2 Gibbs Sampler

To estimate the posterior distribution of the model parameters, we apply random walk Metropolis (RWM) algorithm using Log-Adaptive Proposals (LAP) as our adaptive tuning strategy (Shaby and Wells, 2010). Since conjugate priors are not available, we use random walk Metropolis-Hastings update steps.

At each MCMC iteration, we first update the smooth process $X^{*}$ conditioning on the true values for all non-censored sites, current values for $X^{*}$ and all other model parameters:

where the likelihood function of $X^{*}$ conditioning on the random scaling factor $R$ is calculated in Appendix A.2. We then update $R$ using its conditional posterior distribution

Since time independence is assumed, we can update $\boldsymbol{X}^{*}_{t}$ and $R_{t}$ in a parallel fashion across $t=1,\ldots,T$. The other parameters are updated similarly using adaptively-tuned random walk Metropolis-Hastings updates, with the likelihood (14) multiplied by the corresponding priors in (15).

Simulation Studies

In this section, we present simulation results and conduct coverage analysis to investigate, firstly, whether the MCMC procedure is able to draw accurate inference on model parameters, and secondly, in the case of the modified version of model (5), to check whether our model captures asymptotic dependence characteristics correctly even when the data-generating model is different from the fitted model.

To verify the accuracy of inference made by MCMC sampling, we generate data from model (8) in Section 3.1, in both the special case of the Huser et al. (2017) model (4) and the special case of the Huser and Wadsworth (2019) model (5), with the addition of nugget terms. In both cases, we use $D=200$ sites uniformly drawn from the unit square $^{2}$, with the latent Gaussian processes $Z(\boldsymbol{s})$ are generated using a Matérn covariance with smoothness parameter $\nu=3/2$. The characteristic length scale parameter is set to $\rho=0.05$ in the case of model (4) and $\rho=0.1$ in the case of model (5).

For model (4), we use $T=20$ independent temporal replications, and set $\beta=0.5$, $\gamma=1$ (which is fixed during estimation), and nugget variance $\tau^{2}=0.2^{2}$. This represents a challenging case, with a fairly long tail and nugget that is small relative to the scale of $X^{*}(\boldsymbol{s})$. For model (5), we use $T=40$ independent temporal replications, and set the nugget variance to $\tau^{2}=3^{2}$. Because the latent $Z(\boldsymbol{s})$ is transformed to Pareto, $\tau^{2}=3^{2}$ is still small compared to the scale of smooth process $X^{*}(\boldsymbol{s})$; see the Supplementary Material for a more in-depth discussion on the effects of $\tau^{2}$. We consider two different scenarios for the dependence parameter $\delta$: $\delta=0.3$ and $\delta=0.7$, corresponding to asymptotic independence and asymptotic dependence respectively. Finally, in all cases the processes are marginally transformed to generalized Pareto distribution with $(u,\sigma,\xi)=(11,1,0)$, where $u$ and $p=0.8$, the proportion of censored observations, are treated as known parameters.

The attenuation constants used in the LAP algorithm are $c_{0}=10,c_{1}=0.8$. The prior for $\rho$ is halfCauchy$(1)$, and the priors for the other parameters are specified in (15), where $(\alpha,\beta)=(0.1,0.1)$ so that the prior for $\tau^{2}$ is fairly noninformative. We ran each MCMC chain for 400,000 iterations and thinned the results by a factor of 10. The parallelism of updating $\boldsymbol{X}^{*}_{t}$ and $R_{t}$ is implemented in R via the foreach routine with doParallel package as a backend (Microsoft Corporation and Weston, 2017).

2 Coverage Analysis

We now study the coverage properties of the posterior inference based on the MCMC sampler for the posterior credible intervals with 100 simulated datasets drawn from the Huser et al. (2017) and Huser and Wadsworth (2019) models, under each of the scenarios described in the previous section.

Figures 2 and 3 shows the empirical coverage rates of highest posterior density credible intervals of several sizes, along with standard binomial confidence intervals. In all cases, we can see that the sampler performed well in generating posterior inference that is well calibrated, with close to nominal frequentist coverage. The coverage for larger $\delta$ is may be slightly different than nominal for large $\alpha$, but overall the results are quite good.

3 Simulation with Mis-specified Models

We now fit our model to data generated from other distributions to validate its ability to capture the tail dependence characteristics under mis-specification. We simulate datasets from models referenced in Section 2.1, and use the sampler described in Section 4.2 to fit model (8) with the Huser and Wadsworth (2019) latent process. The data were generated using four different simulation designs:

Skew-$t$ process from Morris et al. (2017a) with $(a,b)=(6,16)$ (asymptotically dependent);

Gaussian scale mixture from Huser et al. (2017) with $\beta=0$ (asymptotically dependent);

Gaussian scale mixture from Huser et al. (2017) with $\beta=1\text{ or }5$ (asymptotically independent).

For each simulation design, we simulate a single dataset using $D=100$ locations uniformly distributed on $^{2}$ and $T=40$ independent time replicates. The Matérn covariance function with $\nu=3/2$ and $\rho=1$ is again specified for the latent Gaussian processes. For the skew-$t$ process, $R_{t}\sim\text{IG}(3,8)$ to give a $t$ distribution with $6$ degrees of freedom, and $\lambda=3$ to simulate moderate skewness. For the last two designs, the $R_{t}$ were generated as described in (4), with $\gamma=1$.

To obtain good starting values for the latent smooth process $\boldsymbol{X}^{*}$ for MCMC, we first marginally transform the simulated data to noisy scale mixture variables $\boldsymbol{X}$ independently at each location using the following procedure. Following the semi-parametric procedure of Coles and Tawn (1991), we estimate each marginal distribution as a blend of the generalized Pareto distribution function above a high marginal threshold, and the empirical distribution function below that threshold. Fixing initial values for $(\delta,\tau^{2})$, we then transform the margins to noisy scale mixtures via $X_{jt}=F_{X|\delta,\tau^{2}}^{-1}\circ\hat{F}_{\boldsymbol{s}_{j}}(Y_{jt})$. The next step is to run a Metropolis algorithm using the full conditional distribution $\varphi(\boldsymbol{X}^{*}\,|\,\cdots)$ (see (16)) 100 times and save the last random walk states as initial values for $\boldsymbol{X}^{*}$. This procedure is also used for data analysis in Section 6. Finally, with initial values in hand, the datasets from each design were fit using a fully Bayesian approach that simultaneously updates marginal and spatial dependence parameters.

Figure 4 displays the nonparametric and model-based estimates of the upper tail dependence $\chi_{u}(h)$ defined in (1). To generate nonparametric estimates of $\chi_{u}(h)$ at distance $h=\|\boldsymbol{s}_{1}-\boldsymbol{s}_{2}\|$, we look at all pairs of points whose locations are $h$ apart (within some small $\varepsilon$ tolerance), and compute the ratio of empirical probabilities $\hat{\chi}_{u}(h)$. This is similar to an empirical variogram estimator. The nonparametric confidence envelopes are obtained by computing pointwise binomial confidence intervals, pretending that each pair of points is independent from each other pair of points. For parametric estimates, we take samples from the converged MCMC chain, and use parameters from each iteration to simulate $10^{8}$ pairs of $(X(\boldsymbol{s}_{1}),X(\boldsymbol{s}_{2}))$ based on our model to generate a smooth $\chi_{u}(h)$ estimate for each MCMC iteration. Then, combining across MCMC iterations, we compute the pointwise average curves and their credible bands. The results in Figure 4 demonstrate that our model provides a sensible approximation to the extremal dependence structure of the mis-specified models. Borrowing strength across locations, the parametric estimators of $\chi_{u}(h)$ are much more reliable than the nonparametric ones in that they are able to discriminate between the two asymptotic classes via estimating a relatively small number of parameters. Especially for the Huser et al. (2017) models, the dependence strength diminishes gradually as $\beta$ becomes greater, eventually resembling the Gaussian copula. With the extra tail flexibility, our model is able to accurately capture these features.

Data Analysis

We consider daily observations of Fosberg Fire Weather Index (FFWI) from 1974 to 2015 at 93 monitoring stations over parts of the Great Plains, mainly from Central Great Plains to South Texas Plains (Dunn et al., 2012). Figure 9 shows the observation locations as black triangles. The FFWI aims to quantify potential wildfire threat. It is a single number summary calculated from temperature, wind speed, and relative humidity; larger index values signify higher flame lengths and more rapid drying (Fosberg, 1978). Due to human activity and changes in the grassland ecosystem, the Great Plains region is becoming an important high-risk region of large wildfires. According to a in-depth study conducted by Donovan et al. (2017), the average total area burned by fires in Great Plains region between 2005 and 2014 was in the millions of hectares per year. In 2017, 809,380 hectares were lost to wildfires in a single week in Texas, Oklahoma, and Kansas alone (Herskovitz, 2017). As a consequence of the prevailing hot and dry air, wildfires in this region are particularly more concentrated in the spring, feasting on grasses made dry by long-term drought. Modeling the tail behavior of FFWI and studying the extremes of the process could have major implications for wildfire planning.

To ensure the independence over time and avoid seasonal effects, we take the maxima of the FFWI values over ten-day intervals during the spring season. Figure 5 shows the 50-year return levels estimated using the block maxima with 10-year sliding windows for 12 randomly-selected stations. There is no clear evidence for a systematic increase or decrease in the return levels. Although treating meteorological data as constant over time is often problematic, particularly for temperature data, the behavior in Figure 5 suggests that an assumption of constant marginal parameters in time is appropriate.

To account for the physical features of the terrain in the Great Plains, we describe the scale parameter in $\boldsymbol{\theta}_{GPD}$ by the trend surface:

where $\text{lon}(\boldsymbol{s})$ and $\text{lat}(\boldsymbol{s})$ are the longitude and latitude of the stations at which the data are observed. We constrain the joint prior of $(\beta_{0},\beta_{1},\beta_{2})^{\prime}$ such that the support of $\sigma(\boldsymbol{s})$ is the positive real line. We model the shape parameter $\xi(\boldsymbol{s})$ as constant over the spatial domain, as suggested by exploratory analysis (see Appendix B.1).

Similar to the procedure in Section 5.3, to obtain starting values for MCMC, we fit generalized Pareto distributions to model events above the 98% quantile, $u_{98}$, and empirical distributions to those below $u_{98}$, of the time series at each station separately, and then use the fitted models to transform the data to have noisy scale mixture distributions. We then ran the MCMC chain for 50,000 iterations thinned by 10 steps and discarded a burn-in period of 25,000 iterations.

First and foremost, we examine whether the estimate $\delta$ falls within $(0,1/2]$ or $(1/2,1)$ in accordance to whether the data-generating process is asymptotically independent or dependent. Table 1 reports the posterior means and 95% credible intervals for the model parameters. Trace plots for $\delta$ and $\tau$ can be found in Appendix B.2. For this dataset, the MCMC results show that the range of the mixing parameter $\delta$ is close to the interface between the two dependence class, while demonstrating asymptotic dependence. Nonetheless, the value of $\delta$ being close to 1/2 means that $\chi_{u}(h)$ will still decrease with $u$ before eventually reaching a positive limit.

To better evaluate the model fit, we randomly hold out 5 stations for validation purposes, and exclude them from the MCMC analysis; see the red points in Figure 9. We then compare the empirical distributions of mean, and maxima for each time point at the 5 held-out stations with those simulated with parameters from MCMC samples. Though we modeled our data as censored observations, the model may still work a bit further into the center of the distribution. Results are displayed in Figure 6 where we only show values exceeding 80% threshold for mean, and 90% threshold for maxima. We can see that the fit displays a good match against the observed mean and maxima in the upper quantiles.

As a comparison, we change $\boldsymbol{X}^{*}$ to be a $t$ process, and re-fit the model using MCMC. We apply proper scoring rules (Gneiting and Raftery, 2007), log scores and continuously-ranked probability scores (CRPS), to compare the quality of probabilistic forecasts. While running MCMC, we interpolate the latent process $\boldsymbol{X}^{*}$ at the held-out locations for each iteration using the full conditional likelihood. Plugging the predictive draws at the held-out observations into the equation (14), we obtain the log score (simply the log-likelihood) as

where $\{\boldsymbol{r}_{j}\,|\,j=1,\ldots,5\}$ are the validation stations. The left panel of Figure 7 compares the log scores between two models, showing that the transformed Gaussian scale mixture process clearly outperforms the $t$ process. Additionally, we calculate the CPRS (Matheson and Winkler, 1976) for both models,

where $F_{Y}$ is the marginal distribution estimated using parameters using one MCMC iteration, and $y$ is the observed value. The right panel of Figure 7 shows the averaged CRPS for the two models, where our model clearly has better results.

Similar to Figure 4, we show the empirical and model-based values of $\chi_{u}(h)$ and $\bar{\chi}_{u}(h)$ for the block maxima of the FFWI in Figure 8, which confirms that our model captures the extremal spatial dependence in the data quite well. The quantity $\bar{\chi}_{u}(h)$ is an alternative dependence measure useful in the situation $\chi(h)=0$, and is defined as

Recall that the posterior mean for $\delta$ is greater than 0.5, which means $\bar{\chi}_{h}(u)\to 1$ as $u\to 1$. Interestingly, the black dashed curve of the right panel of Figure 8 seems to have a limit less than 1. This is because to attain the correct limit in this case, we would need to compute $\bar{\chi}_{h}(u)$ for values of $u$ that are very close to 1, which is very difficult numerically.

2 Results

To get an idea of what a realization of the fitted process looks like, the left panel of Figure 9 shows one realization of the latent $X(s)$ scale mixture process using parameters from one MCMC iteration. The extreme values are mainly concentrated in two small regions. The right panel shows the same realization, now marginally transformed to the scale of the FFWI values. Since we modeled our data as partially censored observations, the map here only displays the areas where threshold exceedances are observed.

A quantity of great interest is areal exceedance probabilities, which represent the amount of territory simultaneously at extreme risk for wildfire. To obtain a Monte Carlo estimate of these joint probabilities, we use parameters from each MCMC sample to simulate 100 processes (on a $15km\times 15km$ grid), and calculate the total area that has FFWI over a designated threshold. Figure 10 shows the results. These curves represent total area at risk for various FFWI thresholds. The curves for the higher thresholds decay faster than those for the lower thresholds, which confirms that extreme events simultaneously occurring across large areas becomes less common when the threshold increases. This also shows that the joint tail of the fire threat index exhibits a weakening dependence structure, with more extreme events being more localized. It is not possible to capture this behavior using limiting extreme value models like max-stable or generalized Pareto processes.

Discussion

In this paper, we have proposed a new modeling approach, based on the class of transformed Gaussian scale mixture models, which includes those recently proposed by Huser et al. (2017) and Huser and Wadsworth (2019). We added a measurement error to the mixture, hence avoiding the need to calculate the onerous $|\mathcal{C}|$-dimensional Gaussian distribution function when dealing with the censored likelihood. We also circumvent the need to draw from a high-dimensional truncated distribution by treating the smooth process as latent and updating repeatedly using MCMC. In addition to its computational advantages, the presence of a measurement error term may also make the model more realistic for data collected by real-world instruments. Indeed, nugget effects are ubiquitous in environmental statistics, not just to represent measurements errors, but also as a result of small scale effects that are not included in the large scale model for spatial dependence.

Even with the presence of the measurement error, the model is still able to capture qualitatively different types of sub-asymptotic dependence behavior of spatial processes. In the case of our modification of the Huser and Wadsworth (2019) model, a smooth transition between both extremal dependence paradigms takes place in the interior of the parameter space, which enables inference of the dependence class in a simple manner. We proved that all the appealing asymptotic properties found in the original smooth process are preserved in the modified model, regardless of the size of the measurement error variance.

The model allows inference on spatial extreme-value datasets with relatively large numbers of locations. The computational limitations are similar to those of conventional spatial Gaussian process models. We have defined the model conditionally as a Bayesian hierarchical model, for which standard MCMC techniques can be used to fit the data. Computation is facilitated greatly by parallelizing over time $t$ and migrating some basic linear algebra to C/C++ via Rcpp. Even so, the lack of closed form marginal transformations creates a significant (though embarrassingly parallel) computational challenge that scales with the total number of exceedances, rather than the usual case of scaling with the number of spatial locations.

Despite easing computational limitations associated with the Huser et al. (2017) and Huser and Wadsworth (2019) models, our modified versions inherit the same theoretical limitations. Neither model is able to account for the possibility of independence between observations as the distance between sites becomes large, nor are they able to transition from asymptotic dependence at short range to asymptotic independence at longer range. Wadsworth and Tawn (2019) presents an alternative approach to modeling spatial extremes that begins to address these issues.

Another interesting possibility to explore would be to include the nugget term inside the link function $g$. This would result in closed-form marginal distributions in some cases, perhaps making computations easier. However, it would change the dependence structure in ways that would not vanish in the limit, which is not the behavior we were aiming for here, but could be useful nonetheless.

Acknowledgements

We gratefully acknowledge support from NSF grant DMS-1752280 and EPSRC grant EP/P002838/1, along with seed grants from the Institute for CyberScience and the Institute for Energy and the Environment at Pennsylvania State University. Computations for this research were performed on the Pennsylvania State University’s Institute for CyberScience Advanced CyberInfrastructure (ICS-ACI).

Appendix A Technical appendix

For the proof of Proposition 3.1, we begin by recalling useful results from the literature.

The first is Breiman’s lemma, see e.g. Breiman (1965) and Cline and Samorodnitsky (1994), and a corollary for sums of a regularly varying and light-tailed random variables.

Suppose that $Q=ST$ where $\mbox{P}(S>s)\in\mbox{RV}_{-\alpha}$, $\alpha\geq 0$ and $\operatorname{E}(T^{\alpha+\delta})<\infty$ for some $\delta>0$. Then

Suppose that $\mbox{P}(X>x)\in\mbox{RV}_{-\alpha}$, $\alpha\in(0,\infty)$, and $\epsilon$ is a random variable such that $\operatorname{E}(e^{\delta\epsilon})<\infty$ for some $\delta>0$. Then

Since $\mbox{P}(X>x)=:\bar{F}_{X}(x)\in\mbox{RV}_{-\alpha}$, for positive finite $\alpha$, $\mbox{P}(e^{X}>x)=\bar{F}_{X}(\log(x))\in\mbox{RV}_{0}$. By assumption $\operatorname{E}(e^{\delta\epsilon})<\infty$, and so applying Breiman’s Lemma to $e^{X}e^{\epsilon}$ yields

The second result relates to sums of Weibull-tailed variables, i.e., those with survival functions satisfying (10) and the associated density

The following lemma can be verified directly from Theorem 4.1 of Asmussen et al. (2018), by identifying that the conditions in Section 4 of that paper hold for this subclass when $\alpha>1$.

Let $Y_{1},Y_{2}$ be variables with density (20), with regularly varying functions $v_{j},u_{j}$ and parameters $\theta_{j}>0,\alpha_{j}>1$, $j=1,2$. Then the density and survival function of the convolution satisfy

where $\psi_{+}(x)=\theta_{1}q_{1}(x)^{\alpha_{1}}+\theta_{2}q_{2}(x)^{\alpha_{2}}$, with $q_{1}(x),q_{2}(x)$ determined by solving

Finally we note also the following two useful inequalities

where the two lower bounds are equal. We can now prove Proposition 3.1.

1) When $X^{*}$ has a regularly varying tail, Corollary A.1.1 provides that $\mbox{P}(X^{*}+\epsilon>x)\sim\mbox{P}(X^{*}>x)$. For the existence of $\eta_{X^{*}}$, the variable $\min(X_{1}^{*},X_{2}^{*})$ also has a regularly varying tail. Further,

so Corollary A.1.1 thus gives $\mbox{P}(X^{*}_{1}+\max(\epsilon_{1},\epsilon_{2})>x,X^{*}_{2}+\max(\epsilon_{1},\epsilon_{2})>x)\sim\mbox{P}(X^{*}_{1}>x,X^{*}_{2}>x)$, and similarly for the lower bound. Hence $\mbox{P}(X^{*}_{1}+\epsilon_{1}>x,X^{*}_{2}+\epsilon_{2}>x)\sim\mbox{P}(X^{*}_{1}>x,X^{*}_{2}>x)$. Therefore it follows that $\chi_{X^{*}}=\chi_{X}$. If $f(x)\sim g(x)$ for $f(x)\to 0$ then $\log f(x)\sim\log g(x)$, so the result for $\eta$ follows as well.

When $X^{*}$ has a Weibull-like tail with Weibull index $\alpha<1$, then $\mbox{P}(e^{X^{*}}>x)\sim u(\log x)\exp\{-\theta(\log x)^{\alpha}\}\in\mbox{RV}_{0}$, and the rest of the argument follows as above.

2) When $\alpha=1$, $\mbox{P}(e^{X^{*}}>x)\sim u(\log x)x^{-\theta}\in\mbox{RV}_{-\theta}$. Breiman’s lemma now yields

As a consequence, inequality (22) does not lead to a precise asymptotic relationship for $\mbox{P}(X^{*}_{1}+\epsilon_{1}>x,X^{*}_{2}+\epsilon_{2}>x)$, but rather that it is asymptotically bounded within the range

Combining (24) and (25), we get the stated bound for $\chi_{X^{*}}$. It follows also that $\log\mbox{P}(X>x)\sim\log\mbox{P}(X^{*}>x)$ and $\log\mbox{P}(X_{1}>x,X_{2}>x)\sim\log\mbox{P}(X^{*}_{1}>x,X^{*}_{2}>x)$, and hence $\eta_{X^{*}}=\eta_{X}$.

3) Here we use Lemma A.2, where different values of $q_{1},q_{2}$ are found for the three cases $\alpha\in(1,2)$, $\alpha=2$ and $\alpha>2$. To make notation more obvious, we replace $q_{1},q_{2}$ with $q_{*},q_{\epsilon}$ for summation of $X^{*},\epsilon$, and $q_{*}^{\wedge},q_{\epsilon}^{\vee}$ etc., for summation of $\min(X_{1}^{*},X_{2}^{*})$ and $\max(\epsilon_{1},\epsilon_{2})$, for example. The three cases are considered for the value of $\alpha_{*}$; we always have $\alpha_{\epsilon}=2$.

which implies $-\log\mbox{P}(X^{*}+\epsilon>x)\sim-\log\mbox{P}(X^{*}>x)$.

To understand the joint behavior, we again use (22). Note that

Consequently, $\min(\epsilon_{1},\epsilon_{2})$ and $\max(\epsilon_{1},\epsilon_{2})$ both have $\alpha=2$, with different $\theta$. As $\theta_{\epsilon}$ does not affect the leading order behavior of $\psi_{+}$ in (26), the tails of $\min(X^{*}_{1},X^{*}_{2})+\min(\epsilon_{1},\epsilon_{2})$ and $\min(X^{*}_{1},X^{*}_{2})+\max(\epsilon_{1},\epsilon_{2})$ both have $\psi_{+}(x)\sim\theta_{*}^{\wedge}x^{\alpha_{*}}$ in the exponent, where

Consequently $-\log\mbox{P}(X^{*}_{1}+\epsilon_{1}>x,X^{*}_{2}+\epsilon_{2}>x)\sim-\log\mbox{P}(X^{*}_{1}>x,X^{*}_{2}>x)$ and so $\eta_{X^{*}}=\eta_{X}$.

b) When $\alpha_{*}=\alpha_{\epsilon}=2$, solving equations (21) provides

For the margins, and maximum $\theta_{\epsilon}=\theta_{\epsilon}^{\vee}=1/(2\tau^{2})$, whilst for the minimum, $\theta_{\epsilon}^{\wedge}=1/\tau^{2}$. We can now use both sets of inequalities (22) and (23), to give that the range of $\eta_{X}$ is

Some simplification arises upon noting that $\theta_{*}^{\vee}=\theta_{*}$, since

As $\tau^{2}\to 0$, i.e., as the nugget effect disappears, the endpoints converge to $\theta_{*}/\theta_{*}^{\wedge}=\eta_{X^{*}}$. As $\tau^{2}\to\infty$, i.e., as the nugget effect dominates, both endpoints converge to $1/2$.

c) When $\alpha_{*}>\alpha_{\epsilon}=2$, solving equations (21) provides

and $\psi_{+}(x)\sim\theta_{\epsilon}x^{2}$. Using inequality (23),

noting that $\theta_{\epsilon}=1/(2\tau^{2})$ and $\theta_{\epsilon}^{\wedge}=1/\tau^{2}$ leads to $\eta_{X}=1/2$.

A.2 Hierarchical model

Under spatiotemporal setting, we observe $\{Y_{t}(\boldsymbol{s})\,|\,t=1,\cdots,T,\boldsymbol{s}\in\mathcal{S}\}$, and the temporal dependence is ignored. Then the hierarchical model can be described as

where $\boldsymbol{\theta}_{R}$ corresponds to $\delta\in$ for the Huser and Wadsworth (2019) model and $\beta\in[0,\infty)$ for the Huser et al. (2017) model, and $Z_{t}$ is a Gaussian process with Matérn covariance $\mathbf{\Sigma}_{\theta_{C}}$. The priors for $\boldsymbol{\theta}_{R}$ are $\delta\sim U(0,1)$ and $\beta\sim\text{halfCauchy}(1)$ respectively.

In Section 4, we have formulated the full conditional likelihood $\varphi(Y_{t}(\boldsymbol{s})\,|\,\boldsymbol{X}_{t}^{*},\tau^{2},\boldsymbol{\theta}_{R},\boldsymbol{\theta}_{GPD},p)$ for a fixed time and location; see Equation (14). Next we are going to work out the conditional joint likelihood $\varphi(\boldsymbol{X}^{*}_{t}\,|\,R_{t},\boldsymbol{\theta}_{C})=\varphi(X^{*}_{t}(\boldsymbol{s}_{1}),\cdots,X^{*}_{t}(\boldsymbol{s}_{D})\,|\,R_{t},\boldsymbol{\theta}_{C})$ for a fixed time $t$. Keeping the notations same as the main text, we denote $g(\cdot)$ as the link function that transforms the scale of the latent Gaussian process. For the model of Huser et al. (2017), $g(\cdot)$ is an identity function, and thus $\boldsymbol{X}_{t}^{*}$ conditioning on $R_{t}$ remains to be a multivariate Gaussian random variable. For the Huser and Wadsworth (2019) model, the conditional density needs more careful scrutiny.

Let $\boldsymbol{X}^{*}=(X^{*}_{1},\cdots,X^{*}_{D})=(R\cdot g(Z_{1}),\cdots,R\cdot g(Z_{D}))$, where $(Z_{1},\cdots,Z_{D})\sim N(0,\mathbf{\Sigma}_{\boldsymbol{\theta}_{C}})$. For the Huser and Wadsworth (2019) model, the density conditional on $R$ is

Fixing $R$ as a constant, we can apply the chain rule to obtain

Then the Jacobian matrix can be written as

Denote $f_{Z}$ as the density function of $N(0,\mathbf{\Sigma}_{\boldsymbol{\theta}_{C}})$. Then

and we will get the desired result. $\Box$

When updating certain variables for a random walk Metropolis-Hastings step, we calculate the corresponding full conditional likelihood for both current values and proposed values. Since $\varphi(Y_{t}(\boldsymbol{s})\,|\,\boldsymbol{X}_{t}^{*},\tau^{2},\boldsymbol{\theta}_{R},\boldsymbol{\theta}_{GPD},p)$ entails most parameters and latent variables, we will evaluate this density over and over again. Looking back at its analytic form in (14), the most compute-intensive part is to calculate the marginal transformation $T(\cdot)$ as defined in (12), which requires the computation of the marginal quantile function of the noisy process, i.e. $F_{X\,|\,{\boldsymbol{\theta}}_{R},\tau^{2}}^{-1}$. Note that the marginal distribution function $F_{X|{\boldsymbol{\theta}}_{R},\tau^{2}}$ can be obtained through a convolution:

where $\phi_{\tau}$ is the density of $N(0,\tau^{2})$, and $F_{X^{*}|{\boldsymbol{\theta}}_{R}}$ is the marginal distribution function of the smooth process, whose forms are derived in Huser et al. (2017) and Huser and Wadsworth (2019) for the two models of interest in this paper. Because (30) cannot be further simplified, we compute the improper integral numerically using the QUADPACK algorithms which are implemented within the gsl_integration library in C++. The computations in C++ and R are interfaced using the package Rcpp.

To compute $p$th quantile of $X$, i.e. $F_{X|{\boldsymbol{\theta}}_{R},\tau^{2}}^{-1}(p)$, which is required for likelihood function evaluations, we first evaluate the distribution function $F_{X|{\boldsymbol{\theta}}_{R},\tau^{2}}$ at a fine grid of $x$ values. We then perform cubic spline interpolation through the control points to yield a continuous quantile function estimate. Due to the smoothness of the quantile function, numerical experiments showed that this technique suffered no measurable reduction in accuracy relative to the much slower technique of computing quantiles using a numerical root finder.

Appendix B Additional Diagnostics

To examine the trend surfaces of the marginal parameter $\boldsymbol{\theta}_{GPD}$, we fit univariate generalized Pareto distribution to the spring observations at each station over a high threshold $u_{0.8}$. The estimated parameters are plotted in Figure 11. We can see that there is no obvious spatial pattern for the shape parameter, while there is significant longitudinal effect for the scale parameter. This give grounds for applying linear trend surface to scale, and constant surface to shape (see Equation (17)).

B.2 MCMC results

Batch means (Flegal et al., 2010) is a convenient way to compute Monte Carlo standard errors for MCMC outputs. If one divides a Markov chain $\{X_{n}\}$ into $k$ batches of size $b$, the Monte Carlo estimate of $\mu=E(g(X))$ based on $i$th batch can be obtained as follows:

and batch means estimate of the Monte Carlo standard error can be defined as

where $\hat{\mu}$ is the overall Monte Carlo estimate. See Figure 13 for batch means standard errors computed periodically for $\delta$ and $\rho$. Other parameters have similar results. We report the stabilized batch means standard errors in Table 2.

For the data analysis, we used 20 cores from the CyberScience Advanced Cyber Infrastructure at Penn State. On average, each iteration takes approximately 1.76 seconds (CPU time). The effective sample size (ESS) per second is also reported in Table 2, in which the marginal parameters has lower values.

Supplementary Material

In this document, we examine the performance of the MCMC algorithm when $\tau^{2}\rightarrow 0$, i.e., our proposed model (8) converges to the underlying smooth process from Huser and Wadsworth. One may suspect that the Markov chains will mix very poorly when the value of $\tau^{2}$ is small (although when $\tau=0$, an alternative MCMC scheme might be possible, wherein components of $\mathbf{X}^{*}$ could be updated from truncated distributions, albeit more complicated ones). We conduct another simulation study for datasets generated from $\delta=0.5$ (transition point) and various $\tau^{2}$ values ($9,\;3,\;0.5,\;0.0025$). Other parameter settings remain the same as in Section 5.1. The sites and the smooth processes are simulated using the same random number generating seed for each dataset. We want to see the critical point of $\tau^{2}$ for the algorithm to fail.

Figure 14 displays the results from running the MCMC algorithm, with each row showing one case. The red dashed line signifies the true parameter values, and the blue lines indicates the 95% posterior credible intervals. Each MCMC chain was run for 400,000 iterations. Thinning the results by a factor of 10, we show the last 50,000 iterations. We can see that, while the performance of $\delta$ stays stable for different true $\tau^{2}$ values, the Markov chain for $\tau^{2}$ converges slower when its true values becomes smaller. In the case where $\tau^{2}=0.0025$, the true value is outside of the 95% credible interval and the chain mixes very poorly. Figure 15 shows trace plots of an $X^{*}$ at one specific site and time which is censored for all four simulations due to the same collection of sites and underlying smoothing processes. We can see that for larger value of $\tau^{2}$, there are more fluctuations for $X^{*}$ over the threshold, which might be the reason why the MCMC algorithm works better in this case. The dissatisfying performance when $\tau^{2}$ is very small might also have something to do with the particular sampler we used (random walk M-H with a symmetric proposal kernel), given the parameter is so close to the boundary, rather than a problem with the overall modeling scheme.

References