Modeling Extremal Streamflow using Deep Learning Approximations and a Flexible
Spatial Process
Reetam Majumder111North Carolina State University, Brian J. Reich1 and
Benjamin A. Shaby222Colorado State University

###### Abstract

Quantifying changes in the probability and magnitude of extreme flooding
events is key to mitigating their impacts. While hydrodynamic data are
inherently spatially dependent, traditional spatial models such as Gaussian
processes are poorly suited for modeling extreme events. Spatial extreme value
models with more realistic tail dependence characteristics are under active
development. They are theoretically justified, but give intractable
likelihoods, making computation challenging for small datasets and prohibitive
for continental-scale studies. We propose a process mixture model (PMM) which
specifies spatial dependence in extreme values as a convex combination of a
Gaussian process and a max-stable process, yielding desirable tail dependence
properties but intractable likelihoods. To address this, we employ a unique
computational strategy where a feed-forward neural network is embedded in a
density regression model to approximate the conditional distribution at one
spatial location given a set of neighbors. We then use this univariate density
function to approximate the joint likelihood for all locations by way of a
Vecchia approximation. The PMM is used to analyze changes in annual maximum
streamflow within the US over the last 50 years, and is able to detect areas
which show increases in extreme streamflow over time.

Key words: Gaussian process, Max-stable process, Neural networks, Spatial
extremes, Vecchia approximation.

##  1 Introduction

The Intergovernmental Panel on Climate Change released its Sixth Assessment in
2021 and projected an increased frequency of hydroclimatic extremes. In
addition to changes in the mean of climate variables, the impact of climate
change is more severe with changes in the frequency and magnitude of
hydroclimatic extremes. Floods are responsible for huge economic and human
costs (Hirabayashi et al., 2013; Winsemius et al., 2018), and this cost is
projected to increase due to sea level rise and extreme precipitation events
brought about by our changing climate (Winsemius et al., 2018). Effective
prediction of future flooding events is required for water infrastructure
design, but is challenging due to the complexity of flooding events and
uncertain climate predictions (Merz et al., 2014; Condon et al., 2015;
Kundzewicz et al., 2017; François et al., 2019). Extensive research has been
conducted looking at changing climate signals in historical extreme rainfall
(Knox, 1993; Kunkel et al., 2020) and in flooding (Franks, 2002; Sharma et
al., 2018; Blöschl et al., 2019; Milly et al., 2005; Walter, 2010; Meehl et
al., 2000; Hirsch, 2011; Vogel et al., 2011; Archfield et al., 2016). For
example, Hirsch and Ryberg (2012) found a significant change in annual maximum
streamflow (a key measure of flood risk) at 48 of 200 US Geological Survey
(USGS) gauges and spatial clustering in the direction and magnitude of the
changes. As a result, there is a need to account for spatial and temporal
variability (i.e., nonstationarity) in flood frequency patterns when assessing
current and future risk (Merz et al., 2014; Kundzewicz et al., 2014; Milly et
al., 2008, 2015; Vogel et al., 2011; Salas and Obeysekera, 2014).

One approach to projecting flood risk on the basis of extreme streamflow
involves the statistical extrapolation of the spatiotemporal trends observed
in the historical record. Of particular interest is estimating the joint
probability of extremal streamflow at multiple locations, which is useful for
understanding regional flood impacts and assessment to support federal and
state emergency management agencies. For example, Šraj et al. (2016); Dawdy et
al. (2012); Lima et al. (2016) use extreme value analysis (EVA) methods to
model nonstationarity with regressions or hierarchical models for the
relationship between flooding and watershed characteristics and weather.
Classic non-spatial EVA (Coles et al., 2001) begins by isolating the extreme
events of interest. This is done systematically by either selecting all
exceedances over a threshold or computing the block maximum, e.g., the annual
maximum of daily streamflow. A spatial EVA analyzes exceedances or pointwise
maxima (i.e., computed separately at each spatial location) as a stochastic
process over space. Modeling spatial dependence allows for predictions at
ungauged locations and the estimation of the joint probability of extremes at
multiple locations. It also facilitates the borrowing of information across
locations to estimate the marginal distribution at each location, which is
particularly useful for EVA where data are sparse and low-probability events
are of interest, and gives valid statistical inference for model parameters by
properly accounting for spatial dependence.

In this study, we consider extreme streamflow data from the United States
Geological Survey’s Hydro-Climatic Data Network (HCDN) (Lins, 2012). Our
primary objective is to identify regions within the US where the distribution
of extreme streamflow has changed over time. The HCDN has a long historical
record, and consists of locations that are minimally impacted by anthropogenic
activity while excluding sites where human activities affect the flow of the
watercourse. We focus on the modeling of block maxima of streamflow with the
help of the max-stable process (MSP) (De Haan et al., 2006). MSPs are a
limiting class of models for spatial extremes, featuring strong forms of tail
dependence (Smith, 1990; Tawn, 1990; Schlather, 2002; Kabluchko et al., 2009;
Wadsworth and Tawn, 2012; Reich and Shaby, 2012). They are a natural
asymptotic model for block maxima, but can also be applied to peaks over a
threshold using a censored likelihood (e.g., Reich et al., 2013; Huser and
Davison, 2014).

In practice, MSPs pose two challenges. First, the analytic forms of (censored)
MSP densities are computationally intractable for all but a small number of
spatial locations (Schlather, 2002; Kabluchko et al., 2009; Wadsworth and
Tawn, 2012, 2014; Wadsworth, 2015). For general MSPs, Castruccio et al. (2016)
stated that full likelihood inference seemed limited to $n=13$ locations.
These low dimensional results have led to the use of composite likelihood (CL)
approximations (Padoan et al., 2010). However, CL suffers from statistical
inefficiency for large $n$ (Huser et al., 2016), finite-sample bias when using
all pairs of observations (Sang and Genton, 2014; Wadsworth, 2015; Castruccio
et al., 2016), and computational challenges posed by computing likelihoods at
all $O(n^{2})$ pairs. More recently, Huser et al. (2019) proposed an
expectation-maximization algorithm for full likelihood inference, with
computation time of 19.8 hours with $n=20$. Huser et al. (2022) have also
applied the Vecchia approximation that requires only moderate-dimensional (say
10 or 15) joint distribution functions, which are available for some MSPs.
More recently, deep learning has been used to estimate parameters in spatial
models within a framework of simulation-based inference. They leverage a
likelihood-free approach by simulating datasets with different parameter
values and using deep learning to identify features of the simulated data that
are indicative of particular parameter values. Gerber and Nychka (2021) used
it to estimate covariance parameters for spatial Gaussian process (GP) models,
by training convolutional neural networks (CNNs) to take moderate size spatial
fields or variograms as input and return the range and noise-to-signal
covariance parameters as output. Lenzi et al. (2021) used simulated data as
input and trained CNNs to learn the parameters of an MSP. Finally, Sainsbury-
Dale et al. (2023) have used permutation invariant neural networks for large
spatial extremes datasets in a Bayesian setting for estimating parameters from
independent replicates. However, it is difficult to extend them to problems
with large numbers of parameters; for example, a crucial assumption in our
application is that the marginal distributions have spatio-temporally varying
coefficients (STVC) which substantially expands the parameter space. Bayesian
approaches have also been proposed since they provide stability by
incorporating prior information as available and are often preferred for
uncertainty quantification. But they are restricted to either small $n$
(Ribatet et al., 2012) or very specific models (Reich and Shaby, 2012; Morris
et al., 2019; Bopp et al., 2021). For lower dimensional problems, approximate
Bayesian computation (ABC) can replace likelihood evaluation with repeated
simulation from the MSP model (Erhardt and Smith, 2012). For general
intractable likelihood estimation problems, neural networks can be leveraged
for conditional density estimation. For example, Greenberg et al. (2019) use a
sequential neural posterior estimation method for simulation-based inference.
A related method is normalizing flows (Papamakarios et al., 2021; Kobyzev et
al., 2021), where a simple density is pushed through a series of
transformations, often involving neural networks, to obtain more complex
densities.

A second challenge posed by MSPs is that they are restrictive in the class of
dependence types they can incorporate. Environmental data often have weakening
spatial dependence with increasing levels of extreme quantiles as we go
farther out into the tails of the distributions; however, MSPs are unable to
accommodate this behavior. Wadsworth and Tawn (2012) addressed this with a
max-mixture model that took an MSP and incorporated asymptotic independence at
the boundary point of the parameter space using a mixing parameter. A more
general approach was taken in Huser and Wadsworth (2019) which combined a
Pareto random variable with a GP resulting in a hybrid model which
interpolates between perfect dependence and asymptotic independence, indexed
similarly by a mixing parameter. This flexible model can establish asymptotic
dependence or asymptotic independence from the data without needing a prior
assumption. A limitation of this model is that the Pareto random variable is
shared by the spatial locations, inducing dependence between distant sites.
This might be unrealistic for an analysis over a large spatial domain.
Finally, Hazra et al. (2021) consider a mixture of a GP with a stochastic
scale process; it can capture a range of extremal dependence structures, but
does not employ a mixing parameter and therefore assumes equal contribution
from both its constituent processes.

In this paper, we propose a spatial EVA model and an associated computational
algorithm to address the aforementioned limitations of the MSP and related
approaches. The EVA model is specified as a convex combination of an MSP and a
GP for residual dependency, and has GEV margins with STVC. We refer to it as
the process mixture model (PMM). From a modeling perspective, the mixture of
the two spatial processes allows asymptotic dependence or independence for
locations separated by distance $h$, with independence as $h\to\infty$ (long-
range independence). Furthermore, the STVC can account for temporal
nonstationarity which is key for large-scale climate studies. This flexibility
comes at a computational cost: the model has hundreds of parameters and even
bivariate PDFs do not have a closed form to the best of our knowledge.
Therefore we develop a new computational algorithm that uses a feed-forward
neural network (FFNN) embedded in a density regression model (Xu and Reich,
2021) to approximate the conditional distribution at one spatial location
given a set of neighbors. Following this, the univariate density functions are
used to approximate the joint likelihood for all locations by means of a
Vecchia approximation (Vecchia, 1988). This specification partitions the
parameter space into a low dimensional vector of spatial dependence parameters
and a higher dimensional vector of marginal parameters, and decouples the
likelihood approximation from parameter estimation. The FFNN is trained on
synthetic data generated from a design distribution using different parameter
values; this allows us to avoid data scarcity issues and accommodate a range
of marginal densities. Parameter estimation is carried out using MCMC. This
computational framework is quite general. Unlike many of the approaches
mentioned above, it can be applied to virtually any spatial process (e.g., GP,
MSP, and mixtures), can accommodate high-dimensional STVC margins, as well as
missing and censored data. We use the PMM to analyze changes in annual maximum
streamflow within the US over the past 50 years.

The rest of this paper is organized as follows. Section 2 provides background
on the construction and dependence measures for MSPs. Section 3 introduces the
PMM for block maxima. Section 4 describes inference for the PMM which employs
a deep learning Vecchia approximated density regression approach. Section 5
consists of a detailed simulation study demonstrating the method. Section 6
analyzes annual streamflow maxima data for HCDN stations across the US and
identifies changes in their behavior over the past 50 years. Section 7
concludes with a discussion. Additional theoretical details, simulation
studies, and results from our application are provided in a supplement.

##  2 Background

###  2.1 The max-stable process

A random process $\\{R(\mbox{\bf s}):\mbox{\bf
s}\in\mathcal{S}\subset\mathbb{R}^{d}\\}$ indexed by spatial locations s is
called max-stable if there exists a sequence $\\{X_{i}(\mbox{\bf
s}):i\in\mathbb{N}\\}$ of independent copies of the process $\\{X(\mbox{\bf
s}):\mbox{\bf s}\in\mathcal{S}\\}$, and normalizing functions $a_{n}(\mbox{\bf
s})>0,b_{n}(\mbox{\bf s})\in\mathbb{R}$ such that

$\displaystyle R(\mbox{\bf s})\stackrel{{\scriptstyle
d}}{{=}}\frac{\max_{i=1:n}X_{i}(\mbox{\bf s})-b_{n}(\mbox{\bf
s})}{a_{n}(\mbox{\bf s})}.$

Further, it can be shown that if there exist continuous functions
$c_{n}(\mbox{\bf s})>0,d_{n}(\mbox{\bf s})\in\mathbb{R}$, such that as
$n\to\infty$,

$\displaystyle\frac{\max_{i=1:n}X_{i}(\mbox{\bf s})-d_{n}(\mbox{\bf
s})}{c_{n}(\mbox{\bf s})}\rightarrow R(\mbox{\bf s}),$

then $R(\mbox{\bf s})$ is either degenerate, or an MSP (de Haan, 1984). If
$R(\mbox{\bf s})$ is non-degenerate, the pointwise distributions of
$R(\mbox{\bf s})$ are in the GEV family (de Haan and Ferreira, 2006).

Max-stable processes arise as the pointwise maxima taken over an infinite
number of appropriately rescaled stochastic processes, and are therefore
widely applied for studying spatial extremes, in particular block maxima (with
a block size of $n$). MSPs can be constructed through a spectral
representation (de Haan, 1984; Penrose, 1992). Let
$\\{Z_{i}:i\in\mathbb{N}\\}$ be the points of a Poisson process on
$(0,\infty)$ with intensity $1/z^{2}$. Then there exists a non-negative
stochastic process $W(\mbox{\bf s})$ with continuous sample paths and with
$\mathbb{E}W(\mbox{\bf s})=1$ for all $\mbox{\bf s}\in\mathcal{S}$, such that

$\displaystyle R(\mbox{\bf s})\stackrel{{\scriptstyle d}}{{=}}\max_{i\geq
1}Z_{i}W_{i}(\mbox{\bf s}),$

where $W_{i}(\mbox{\bf s})$ are independent copies of $W(\mbox{\bf s})$.
Common parametric sub-classes of MSPs include mixed moving maxima processes
(Wang and Stoev, 2010), the Schlather processes (Schlather, 2002), and Brown-
Resnick processes (Brown and Resnick, 1977; Kabluchko et al., 2009).

The finite dimensional distribution of an MSP $R(\mbox{\bf s})$ at a set of
locations $(\mbox{\bf s}_{1},\ldots,\mbox{\bf s}_{k})\in\mathcal{S}$ has the
form $Pr\\{R(\mbox{\bf
s}_{j})<r_{j},j=1:k\\}=\exp\\{-V(r_{1},\ldots,r_{k})\\}$, where $V$ is known
as the exponent function and is given by:

$\displaystyle
V(r_{1},\ldots,r_{k})=\mathbb{E}\biggl{[}\max_{j=1:k}\frac{W(\mbox{\bf
s}_{j})}{r_{j}}\biggr{]}.$

###  2.2 Dependence properties

Let $F_{1}$ and $F_{2}$ be the CDFs of $R(\mbox{\bf s}_{1})$ and $R(\mbox{\bf
s}_{2})$, and let $U(\mbox{\bf s}_{i})=F_{i}(R(\mbox{\bf s}_{i}))$ for
$i=1,2$. The joint tail behavior of the two random variables $U(\mbox{\bf
s}_{1})$ and $U(\mbox{\bf s}_{2})$ with uniform marginals can be studied in
terms of the conditional exceedance probability, given by:

$\displaystyle\chi_{u}(\mbox{\bf s}_{1},\mbox{\bf s}_{2}):=Pr[U(\mbox{\bf
s}_{1})>u|U(\mbox{\bf s}_{2})>u]\in(0,1),$

where $u\in(0,1)$ is a threshold. A commonly used measure of extremal
dependence is the upper-tail coefficient (Joe, 1997), defined as:

$\displaystyle\chi(\mbox{\bf s}_{1},\mbox{\bf s}_{2})=\lim_{u\rightarrow
1}\chi_{u}(\mbox{\bf s}_{1},\mbox{\bf s}_{2}).$

The random variables $U(\mbox{\bf s}_{1})$ and $U(\mbox{\bf s}_{2})$ are
considered asymptotically dependent if the upper-tail coefficient is strictly
positive, and asymptotically independent if it is zero.

We note that asymptotic (in)dependence is different from complete
(in)dependence, since asymptotic (in)dependence is determined specifically by
the joint behavior of the tails of the distribution as $u\to 1$. Asymptotic
independence is also different from long-range independence which is
determined by the asymptotic behavior of $\chi_{u}(\mbox{\bf s}_{1},\mbox{\bf
s}_{2})$ as $||\mbox{\bf s}_{1}-\mbox{\bf s}_{2}||\to\infty$.

In the asymptotic independence scenario, the coefficient of tail dependence
proposed by Ledford and Tawn (1996, 1997) is useful to study the joint tail
behavior of the process. Consider $R(\mbox{\bf s})$ with unit Fréchet margins.
The joint survivor function of $R(\mbox{\bf s}_{1})$ and $R(\mbox{\bf s}_{2})$
can be expressed as:

$\displaystyle\bar{F}(r,r):=Pr[R(\mbox{\bf s}_{1})>r,R(\mbox{\bf
s}_{2})>r]\sim\mathcal{L}(r)r^{-1/\eta},\mbox{ as }r\to\infty,$

where $\mathcal{L}$ is a slowly varying function that satisfies
$\mathcal{L}(tr)/\mathcal{L}(r)\to 1$ as $r\to\infty$ for all fixed $t>0$, and
$\eta\in(0,1]$ is a constant that effectively determines the decay rate of
$\bar{F}(r,r)$ for large $r$. The parameter $\eta$ is known as the coefficient
of tail dependence. A value of $\eta=1/2$ indicates independent marginal
variables; values lower and higher than $1/2$ correspond to a negative and a
positive association respectively between the pair of variables.

###  2.3 The Brown-Resnick process

Consider $W(s)=\exp\\{\epsilon(\mbox{\bf s})-\gamma(\mbox{\bf s})\\}$ in the
spectral representation of an MSP, where $\epsilon(\mbox{\bf s})$ is an
intrinsically stationary Gaussian process with semivariogram $\gamma(\cdot)$,
and $\epsilon(0)=0$ almost surely. $W(\mbox{\bf s})$ is continuous and non-
negative. Then $R(\mbox{\bf s})$ is a strictly stationary MSP known as the
Brown-Resnick process, whose distribution depends only on $\gamma(\cdot)$. The
Brown-Resnick process is attractive since a wide range of variograms can be
used with them, and they are relatively easily simulated. The exponent
function that defines the joint distribution for the pair $R(\mbox{\bf
s}_{1})$ and $R(\mbox{\bf s}_{2})$ is

$\displaystyle
V(r_{1},r_{2})=\frac{1}{r_{1}}\Phi\biggl{\\{}\frac{a}{2}-\frac{1}{a}\log\biggl{(}\frac{r_{1}}{r_{2}}\biggr{)}\biggr{\\}}+\frac{1}{r_{2}}\Phi\biggl{\\{}\frac{a}{2}-\frac{1}{a}\log\biggl{(}\frac{r_{2}}{r_{1}}\biggr{)}\biggr{\\}},$

where $a=\\{2\gamma(\mbox{\bf s}_{1}-\mbox{\bf s}_{2})\\}^{1/2}$, and
$\Phi(\cdot)$ denotes the standard normal distribution function.

##  3 A Process Mixture Model for Spatial Extremes

Let $\\{Y(\mbox{\bf s});\mbox{\bf s}\in\mathcal{S}\\}$ be a spatial extremes
process indexed by the set $\mathcal{S}\subset\mathbb{R}^{2}$. In this section
we consider $Y(\mbox{\bf s})$ to be defined as a block maximum, but the
methods can be extended to peaks over a threshold. We assume a potentially
different marginal distribution for each spatial location s and denote
$F_{\mbox{\bf s}}$ as the marginal cumulative distribution function (CDF) for
site s. For example, we assume that $F_{\mbox{\bf s}}$ is the generalized
extreme value (GEV) distribution with location $\mu(\mbox{\bf
s})\in\mathbb{R}$, scale $\sigma(\mbox{\bf s})>0$, and shape $\xi(\mbox{\bf
s})\in\mathbb{R}$, so that marginally

$Y(\mbox{\bf s})\sim\mbox{GEV}\\{\mu(\mbox{\bf s}),\sigma(\mbox{\bf
s}),\xi(\mbox{\bf s})\\}.$

Its CDF $F_{\bf s}(y|\mu(\mbox{\bf s}),\sigma(\mbox{\bf s}),\xi(\mbox{\bf
s})):=Pr[Y(\mbox{\bf s})<y]$ is

$F_{\bf s}(y|\mu(\mbox{\bf s}),\sigma(\mbox{\bf s}),\xi(\mbox{\bf
s}))=\begin{cases}\exp\biggl{[}-\left\\{1+\xi(\bf s)\left(\frac{y-\mu(\bf
s)}{\sigma(\bf s)}\right)\right\\}_{+}^{-1/\xi(\bf s)}\biggr{]},&\xi(\mbox{\bf
s})\neq 0,\\\ \exp\left\\{-\exp\bigl{(}\frac{y-\mu(\bf s)}{\sigma(\bf
s)}\bigl{)}\right\\},&\xi(\mbox{\bf s})=0,\end{cases}$ (1)

with $\\{y\\}_{+}:=\max(0,y)$, and support over the set
$\bigr{\\{}y:1+\xi(\mbox{\bf s})(y-\mu(\mbox{\bf s}))/\sigma(\mbox{\bf
s})>0\bigr{\\}}$ for the CDF. The shape parameter $\xi(\mbox{\bf s})$ controls
the lower and upper bounds of the distribution; the GEV distribution is
bounded above for $\xi(\mbox{\bf s})<0$, and bounded below for $\xi(\mbox{\bf
s})>0$. Therefore, the transformed variables

$U(\mbox{\bf s})=F_{\bf s}\\{Y(\mbox{\bf s})\\}$ (2)

share common uniform marginal distributions across the spatial domain. This
transformation separates residual spatial dependence in $U(\mbox{\bf s})$ from
the spatial dependence induced by spatial variation in the GEV parameters,
which we model using GP priors over s. We note that although we describe the
marginal and residual models separately, we fit a joint hierarchical model to
simultaneously estimate all model parameters.

We define our spatial dependence model on the residual model $U(\mbox{\bf s})$
by taking $U(\mbox{\bf s})=G\\{V(\mbox{\bf s})\\}$, such that

$V(\mbox{\bf s})=\delta R(\mbox{\bf s})+(1-\delta)W(\mbox{\bf s}),$ (3)

where $R(\mbox{\bf s})$ and $W(\mbox{\bf s})$ are an MSP and a GP
respectively, both transformed to have standard exponential margins, and
$\delta\in[0,1]$ is the weight parameter to control relative contribution of
the two spatial processes. Mixing the asymptotically dependent MSP with the
asymptotically independent GP provides a rich model for spatial dependence.
This generalizes the extremal process of Huser and Wadsworth (2019), who
assumed a standard exponential random variable $R$ common to all locations, by
replacing it with an MSP. Since (3) mixes two processes, we refer to it as the
process mixture model (PMM). Further details regarding the transformations
required to obtain standard exponential margins, and the Huser-Wadsworth model
are provided in the Supplementary Material (Majumder et al., 2023, Appendix
A.2–A.3).

Let $G(\cdot)$ denote the CDF of the marginal distribution of $V(\mbox{\bf
s})$. By construction, $V(\mbox{\bf s})$ marginally follows the two-parameter
hypoexponential distribution, and its CDF has the following functional form:

$G(v)=1-\frac{1-\delta}{1-2\delta}\exp\left\\{-\frac{v}{1-\delta}\right\\}+\frac{\delta}{1-2\delta}\exp\left\\{-\frac{v}{\delta}\right\\},$
(4)

where $\delta\in(0,1),\delta\neq 1/2$, and $v>0$. Although other options are
possible, we model the correlation of $W(\mbox{\bf s})$ using the isotropic
powered-exponential correlation function $\mbox{Cor}\\{W(\mbox{\bf
s}_{1},\mbox{\bf s}_{2})\\}=\exp\\{-(h/\rho_{W})^{\alpha_{W}}\\}$ with
distance $h=||\mbox{\bf s}_{1}-\mbox{\bf s}_{2}||$ measured as the $L_{2}$
norm, smoothness $\alpha_{W}\in(0,2]$, and range $\rho_{W}>0$. $R(\mbox{\bf
s})$ is assumed to be an isotropic Brown-Resnick process with the variogram
$\gamma(h)=(h/\rho_{R})^{\alpha_{R}}$, for smoothness $\alpha_{R}\in(0,2]$ and
range $\rho_{R}>0$. We also incorporate a nugget into the process mixture.
Denoting the proportion of the variance explained by the spatial process by
$r$, we construct $W(\mbox{\bf s})$ and $R(\mbox{\bf s})$ to satisfy:

$\displaystyle\mbox{Cor}\bigl{(}W(\mbox{\bf s}_{1}),W(\mbox{\bf
s}_{2})\bigr{)}$ $\displaystyle=r\cdot\exp\\{-(h/\rho_{W})^{\alpha_{W}}\\}$
$\displaystyle R(\mbox{\bf s})$ $\displaystyle=\max\\{r\cdot R_{1}(\mbox{\bf
s}),(1-r)\cdot R_{2}(\mbox{\bf s})\\},$

where $R_{1}(\mbox{\bf s})$ is an MSP, and $R_{2}(\mbox{\bf
s})\stackrel{{\scriptstyle iid}}{{\sim}}\mbox{GEV}(1,1,1)$ distributed
independently of $R_{1}(\mbox{\bf s})$. Since $R(\mbox{\bf s})$ and
$W(\mbox{\bf s})$ are assumed to be isotropic processes, going forward, we can
rewrite $\chi_{u}(\mbox{\bf s}_{1},\mbox{\bf s}_{2})$ and $\chi(\mbox{\bf
s}_{1},\mbox{\bf s}_{2})$ as functions of the distance between locations,
$\chi_{u}(h)$ and $\chi(h)$, respectively.

(a) $\chi_{u}(h)$ as a function of $u$ and $\delta$, at distance $h=0.22$.

(b) $\chi_{u}(h)$ as a function of $h$ and $\delta$, for threshold $u=0.999$.

Figure 1: Behavior of the empirical conditional exceedance: Approximate
$\chi_{u}(h)$ for the PMM plotted as a function of threshold $u$, distance
$h$, and asymptotic dependence parameter $\delta$. Smoothness parameters
$\alpha_{W}=\alpha_{R}=1$, GP range $\rho_{W}=0.2$, and MSP range
$\rho_{R}=0.1$ are fixed for both plots.

Figure 1 plots Monte Carlo approximations of $\chi_{u}(h)$ as a function of
$u$ and $h$ for the PMM. For these plots we fix $\rho_{R}=0.1$,
$\rho_{W}=0.2$, $\alpha_{R}=\alpha_{W}=1$, and
$\delta\in\\{0.2,0.4,0.6,0.8\\}$. Figure 1(a) sets the correlation to $0.4$ by
fixing $h=0.22$, and plots $\chi_{u}(h)$ as a function of the threshold $u$.
The limit is zero for $\delta<0.5$ and positive for $\delta>0.5$. For small
values of $h$, $R(\mbox{\bf s})$ is approximately the same for both sites
(i.e., $R(\mbox{\bf s}_{1})\approx R(\mbox{\bf s}_{2})=R$), and thus the
univariate $R$ result of Huser and Wadsworth (2019) that the process is
asymptotically dependent if and only if $\delta>0.5$ emerges. An analytical
expression for $\chi_{u}(h)$ in this special case is provided in Majumder et
al. (2023, Appendix A.4). From Figure 1(b), we see that as the distance $h$
increases, $\chi_{u}(h)$ converges to zero for all $\delta$ because both
$R(\mbox{\bf s})$ and $W(\mbox{\bf s})$ have diminishing spatial dependence
for long distances. The rate of convergence of $\chi_{u}(h)$ to zero also
depends on the value of $\delta$, with much slower convergence when
$\delta>0.5$. We note that $\chi_{u}(h)$ does not converge to zero for large
$h$ under the common $R$ model of Huser and Wadsworth (2019), which is
unrealistic for studies on a large spatial domain. Additional plots for
different values of $\rho_{R}$ and $\rho_{W}$ are provided in Majumder et al.
(2023, Appendix A.5).

##  4 Deep Learning Vecchia Approximation for the Process Mixture Model

Fitting the PMM introduced in Section 3 poses computational challenges,
especially for large datasets. The joint distribution for $W(\mbox{\bf s})$ is
available in closed form but is cumbersome for large datasets; the joint
distribution of $R(\mbox{\bf s})$ is available only for a moderate number of
spatial locations, and the joint distribution of the mixture model is more
complicated that either of its components. An alternative is to build a
surrogate likelihood for Bayesian computation (e.g., Rasmussen, 2003;
Wilkinson, 2014; Price et al., 2018; Wang and Li, 2018; Järvenpää et al.,
2021; Li et al., 2019). Below we develop a surrogate likelihood based on a
Vecchia decomposition (Vecchia, 1988) and deep learning density regression.

Assume the process is observed at $n$ locations $\mbox{\bf
s}_{1},...,\mbox{\bf s}_{n}$. Partition the parameters into those that affect
the marginal distributions, denoted $\mbox{\boldmath$\theta$}^{MARG}$, and
those that affect the spatial dependence, denoted
$\mbox{\boldmath$\theta$}^{SPAT}$. For the model in Section 3,
$\mbox{\boldmath$\theta$}^{MARG}$ includes the GEV parameters
$\mbox{\boldmath$\theta$}^{MARG}=\\{\mu(\mbox{\bf s}_{i}),\sigma(\mbox{\bf
s}_{i}),\xi(\mbox{\bf s}_{i});i=1,...,n\\}$ and
$\mbox{\boldmath$\theta$}^{SPAT}=\\{\delta,\rho_{R},\alpha_{R},\rho_{W},\alpha_{W}\\}$.
Let $Y(\mbox{\bf s}_{i})\equiv Y_{i}$ and
$U_{i}=F(Y_{i};\mbox{\boldmath$\theta$}^{MARG})$ be the transformation of the
response so that the distribution of $U_{i}\in[0,1]$ does not depend on
$\mbox{\boldmath$\theta$}^{MARG}$. We approximate the spatial model on this
scale and use the standard change of variables formula to define the joint
likelihood on the original scale

$f_{y}(y_{1},...,y_{n};\mbox{\boldmath$\theta$}^{MARG},\mbox{\boldmath$\theta$}^{SPAT})=f_{u}(u_{1},...,u_{n};\mbox{\boldmath$\theta$}^{SPAT})\prod_{i=1}^{n}\left|\frac{dF(y_{i};\mbox{\boldmath$\theta$}^{MARG})}{dy_{i}}\right|,$
(5)

where $f_{y}(\cdot)$ and $f_{u}(\cdot)$ are the joint density functions of
$Y_{1}\ldots,Y_{n}$ and $U_{1},\ldots,U_{n}$ respectively.

We approximate the joint likelihood in (5) using a Vecchia approximation
(Vecchia, 1988; Stein et al., 2004; Datta et al., 2016; Katzfuss and Guinness,
2021),

$f_{u}(u_{1},...,u_{n};\mbox{\boldmath$\theta$}^{SPAT})=\prod_{i=1}^{n}f(u_{i}|\mbox{\boldmath$\theta$}^{SPAT},u_{1},...,u_{i-1})\approx\prod_{i=1}^{n}f_{i}(u_{i}|\mbox{\boldmath$\theta$}^{SPAT},u_{(i)})$
(6)

for $u_{(i)}=\\{u_{j};j\in{\cal N}_{i}\\}$ and ${\cal
N}_{i}\subseteq\\{1,...,i-1\\}$, e.g., the $m$ locations in ${\cal N}_{i}$
that are closest to $\mbox{\bf s}_{i}$. The set of locations $\mbox{\bf
s}_{(i)}$ are analogously defined as $\mbox{\bf s}_{(i)}=\\{\mbox{\bf
s}_{j};j\in{\cal N}_{i}\\}$; the set is referred to as the Vecchia neighboring
set, and its members as the Vecchia neighbors of location $\mbox{\bf s}_{i}$.
Of course, not all locations that are dependent on location $\mbox{\bf s}_{i}$
need be included in ${\cal N}_{i}$ because distant observations may be
approximately independent after conditioning on more local observations. The
approximation therefore entails truncating the dependence that $u_{i}$ has on
all its previous $i-1$ ordered sites to instead consider dependence on only up
to $m$ sites, i.e., $|{\cal N}_{i}|\leq m$. The first term on the right hand
side of (6) is the marginal density $f_{1}(u_{1})$. Different choices are
possible for ordering the locations prior to the Vecchia approximation
(Guinness, 2018). In our work, the spatial locations are scaled to be on the
unit square, and ordered by their distance from the origin.

The conditional distributions for the PMM do not have closed-form expressions.
We consider two related approximations - the local approximation where
individual conditional density functions $f_{i}(\cdot)$ are fit for each
location, and a global approximation where a single conditional density
function $f(\cdot)$ is estimated for all locations $\mbox{\bf
s}_{i},i=2,\ldots,n$. We primarily focus on the local SPQR and present it in
this section. Details of the global SPQR approach can be found in the
Supplementary Material (Majumder et al., 2023, Appendix A.6).

For the local SPQR approximation at location $\mbox{\bf s}_{i}$, we fit a
density regression viewing $u_{(i)}$ and $\mbox{\boldmath$\theta$}^{SPAT}$ as
features (covariates), denoted by $\mbox{\bf x}_{i}$. We approximate the
univariate conditional densities for density regression using the model
introduced in Xu and Reich (2021):

$f_{i}(u_{i}|\mbox{\bf x}_{i},\mathcal{W})=\sum_{k=1}^{K}\pi_{ik}(\mbox{\bf
x}_{i},\mathcal{W}_{i})B_{k}(u_{i}),$ (7)

for $i=2,\ldots,n$, where $\pi_{ik}(\mbox{\bf x}_{i},\mathcal{W}_{i})\geq 0$
are probability weights with $\sum_{k=1}^{K}\pi_{ik}(\mbox{\bf x}_{i})=1$ that
depend on the parameters $\mathcal{W}_{i}$ and $B_{k}(u_{i})\geq 0$ are
M-spline basis functions that, by definition, satisfy $\int B_{k}(u)du=1$ for
all $k$. By increasing the number of basis functions $K$ and appropriately
selecting the weights $\pi_{ik}(\mbox{\bf x}_{i})$, this mixture distribution
can approximate any continuous density function (e.g., Chui et al., 1980;
Abrahamowicz et al., 1992). The weights are modeled using a feed-forward
neural network (FFNN) with $H$ hidden layers with $N_{l}$ neurons in hidden
layer $l$ and multinomial logistic weights, equivalent to the softmax
activation function. The model is

$\displaystyle\pi_{ik}(\mbox{\bf x}_{i},\mathcal{W}_{i})$ $\displaystyle=$
$\displaystyle\frac{\exp\\{\gamma_{Hk}(\mbox{\bf
x}_{i},\mathcal{W}_{i})\\}}{\sum_{j=1}^{K}\exp\\{\gamma_{Hj}(\mbox{\bf
x}_{i},\mathcal{W}_{i})\\}},$ (8) $\displaystyle\gamma_{lk}(\mbox{\bf
x}_{i},\mathcal{W}_{i})$ $\displaystyle=$ $\displaystyle
W_{ilk0}+\sum_{j=1}^{N_{l}}W_{ilkj}\psi\left\\{\gamma_{l-1,j}(\mbox{\bf
x}_{i},\mathcal{W}_{i})\right\\}\mbox{\ \ \ for\ \ \ }l\in\\{1,...,H\\},$
$\displaystyle\gamma_{0k}(\mbox{\bf x}_{i},\mathcal{W}_{i})$ $\displaystyle=$
$\displaystyle W_{i0k0}+\sum_{j=1}^{p}W_{i0kj}x_{ij},$

where $\mbox{\bf x}_{i}=(x_{i1},...,x_{ip})$, $\mathcal{W}_{i}=\\{W_{ilkj}\\}$
are parameters to be estimated and $\psi$ is the activation function.
Activation functions are non-linear transformations applied to each output
element of a layer, and is a key feature of neural networks that allows them
to learn complex, non-linear dependencies in the data. The SPQR methodology of
Xu and Reich (2021) admits most of the commonly used activation functions, and
we use the rectified linear unit (ReLU) (Nair and Hinton, 2010). FFNNs use
optimization to obtain optimum values of $\pi_{ik}(\mbox{\bf
x}_{i},\mathcal{W}_{i})$. In SPQR, the FFNN minimizes the negative log-
likelihood loss associated with the density in (7), using the process values
evaluated at locations $\mbox{\bf s}_{i}$ as the response. Building on the
universal approximation theorem for FFNNs (Hornik et al., 1989), Xu and Reich
(2021) argue that with $H=1$ and large $K$ and $N_{1}$, the model in (8) can
approximate any conditional density function that is smooth in its arguments.

Within this framework, approximating the conditional distributions is
equivalent to estimating the weights $\mathcal{W}$. Unlike a typical
statistical learning problem, observational data are not used to estimate
$\mathcal{W}$. Rather, the weights are learned from training data generated
from the PMM with parameters $\mbox{\boldmath$\theta$}^{SPAT}\sim p^{*}$, and
then a realization from the process over sites $\mbox{\bf s}_{i}$ and
$\mbox{\bf s}_{(i)}$ from the model conditioned on
$\mbox{\boldmath$\theta$}^{SPAT}$. Specifically, we generate data at the
observed spatial location with the same neighbor sets to be used in the
analysis. We select the design distribution $p^{*}$ with support covering the
range of plausible values for $\mbox{\boldmath$\theta$}^{SPAT}$. Given these
values, we generate $U(\mbox{\bf s})$ at $\mbox{\bf s}\in\\{\mbox{\bf
s}_{i},\mbox{\bf s}_{(i)}\\}$. The feature set $\mbox{\bf x}_{i}$ for modeling
$u_{i}$ at location $\mbox{\bf s}_{i}$ thus contains the spatial parameters
$\mbox{\boldmath$\theta$}^{SPAT}$ and the process values $U(\mbox{\bf
s}_{(i)})$ at the neighboring locations.

Therefore, all that is required to build the approximation is the ability to
generate small datasets from the model. The size of the training data is
effectively unlimited, meaning the approximation can be arbitrarily accurate.
Once the weights have been learned, applying the FFNN to the approximate
likelihood is straightforward, and the Vecchia approximation ensures that the
computational burden increases linearly in the number of spatial locations.
The proposed estimation approach is in the same vein as recent simulation
based neural inference methods (Gerber and Nychka, 2021; Lenzi et al., 2021;
Sainsbury-Dale et al., 2023); an important distinction is that instead of
estimating model parameters using a neural network, our approach estimates a
set of conditional densities which approximates the full model likelihood.

The weights in (7) are estimated separately for each location. That is, each
component $f_{i}(u_{i}|\mbox{\boldmath$\theta$}^{SPAT},u_{(i)})$ in (6) is
modeled using its own FFNN, for $i=2,\ldots,n$. The model is fit using the R
package SPQR (Xu and Majumder, 2022) and the fitting process is consequently
referred to as the SPQR approximation. The SPQR package supports hardware
acceleration for systems with a CUDA-compatible NVIDIA graphical processing
unit (GPU), which was used for all SPQR models in this paper and provided
significant speedups for computation times. All computations were carried out
on a mobile workstation with 11th Gen Intel Core i7-11800H processors (8
cores, 16 logical processors), 64 GB of RAM, and an NVIDIA T600 laptop GPU
with CUDA support. Algorithm 1 outlines the local SPQR procedure.

Algorithm 1 Local SPQR approximation

Locations $\mbox{\bf s}_{1},\ldots,\mbox{\bf s}_{n}$ and corresponding sets of
neighboring locations $\mbox{\bf s}_{(1)},\ldots,\mbox{\bf s}_{(n)}$

Design distribution $p^{*}$, training sample size $N$

$i\leftarrow 2$

while $i\leq n$ do

$k\leftarrow 1$

while $k\leq N$ do

Draw values of $\mathbf{\mbox{\boldmath$\theta$}}_{k}^{SPAT}\sim p^{*}$

Generate $U_{k}(\mbox{\bf s})$ at $\mbox{\bf s}\in\\{\mbox{\bf
s}_{i},\mbox{\bf s}_{(i)}\\}$ given $\mathbf{\mbox{\boldmath$\theta$}}_{2k}$
using (3)

Define features $\mbox{\bf
x}_{ik}=(\mathbf{\mbox{\boldmath$\theta$}}_{k}^{SPAT},u_{(i)k})$, where
$u_{(i)k}=\\{U_{k}(\mbox{\bf s});\mbox{\bf s}\in\mbox{\bf s}_{(i)}\\}$

$k\leftarrow k+1$

end while

solve $\hat{{\cal W}_{i}}\leftarrow\underset{{\cal
W}}{\operatorname{argmax}}\prod_{k=1}^{N}f_{i}(u_{ik}|\mbox{\bf x}_{ik},{\cal
W})$ for $f_{i}(u_{i}|\mbox{\bf x}_{i},{\cal W}_{i})$ defined in (7) using
SPQR

$i\leftarrow i+1$

end while

Given the approximate model in (5) for $f_{y}$ with an SPQR approximation for
$f_{u}$, a Bayesian analysis using Markov Chain Monte Carlo (MCMC) methods is
straightforward. We use Metropolis updates for both
$\mbox{\boldmath$\theta$}^{MARG}$ and $\mbox{\boldmath$\theta$}^{SPAT}$. For a
spatially-varying coefficient model with local GEV coefficients for location
$\mbox{\bf s}_{i}$, we update parameters $\\{\mu(\mbox{\bf
s}_{i}),\sigma(\mbox{\bf s}_{i}),\xi(\mbox{\bf s}_{i})\\}$ as a block
sequentially by site, and exploit the Vecchia approximation to use only terms
in the likelihood corresponding to sites $j\mbox{ such that
}j\in\mathcal{N}_{j}$, i.e. sites for which site $i$ is included in the
neighboring set. All Metropolis updates are tuned to give acceptance
probability 0.4, and convergence is diagnosed based on the visual inspection
of the trace plots. Additional computational details are given for specific
analyses below, and MCMC code is provided in a GitHub repository -
https://github.com/reetamm/SPQR-for-spatial-extremes.

##  5 Simulation Study

(a) Variable importance of $\delta$, $\rho$, $r$, and the two nearest
neighbors.

(b) Q-Q plots on the exponential scale for checking goodness-of-fit.

Figure 2: Model diagnostics for local SPQR fit on PMM data: VI plot of the 5
most important variables (left) and Q-Q plot on the exponential scale (right)
for the local SPQR at site 25. The gray lines correspond to Q-Q plots of fits
to different simulated datasets.

We illustrate the potential of the proposed algorithm with a simulation study
to evaluate the performance of the method on extreme value data in terms of
both density estimation and parameter estimation. We consider the PMM as the
underlying spatial process, and an STVC model for the marginal GEV parameters.
The goal of the study is to assess whether our approach can simultaneously
model the marginal GEV distributions and the underlying extremal spatial
process. We simulated data from the PMM defined in Section 3 at $n=50$ spatial
locations distributed randomly on the unit square. To put the MSPs and GPs on
the same scale, we assume common smoothness parameter
$\alpha_{R}=\alpha_{W}=\alpha=1$ and parameterize $\rho_{W}$ and $\rho_{R}$ to
give the same effective range, i.e., the distance at which the correlation of
the GP reaches 0.05 ($h=\rho_{W}\log(20)$) and the $\chi$-coefficient of the
MSP reaches 0.05 ($h=\rho_{R}4\Phi^{-1}(1-0.05/2)^{2}$ where $\Phi$ is the
standard normal distribution function). This results in $\rho=\rho_{W}$ and
$\rho_{R}=0.19\rho$. We also assume the absence of a nugget term, i.e., $r=1$.

The locations are ordered by their distance from their origin. For the
$i^{th}$ location $\mbox{\bf s}_{i},\mbox{ with }i>1$ , the Vecchia neighbor
set $\mbox{\bf s}_{(i)}$ consists of the $m$ nearest neighbors of $\mbox{\bf
s}_{i}$ among the previous $i-1$ locations. Up to 15 conditioning points are
used in the Vecchia neighbor set, i.e., $m=\min(i-1,15)$, where
$m=|\mathcal{N}_{i}|$. We use the local SPQR model outlined in Algorithm 1 to
model the conditional densities at each location, which employs stochastic
gradient descent with the adaptive moment estimation (Adam) optimizer (Kingma
and Ba, 2014). We compared multiple SPQR models as part of the density
estimation process. Models were compared on the basis of the log-score and the
Kullback-Leibler (KL) divergence between the estimated and true densities.
Architectures with the lowest validation loss were chosen in each case. Our NN
architecture for each SPQR model consists of two hidden layers with 30 and 15
neurons, 15 output nodes, a learning rate 0.001, batch size 100, and 50
epochs. We train the SPQR model with design distribution $p^{*}$, generating
samples uniformly on $\rho\in(0.0,0.5)$ and $\delta\in(0,1)$.

We first evaluate the SPQR fits of the PMM full conditional distributions.
Figure 2(a) plots variable importances for the local SPQR model at site 25.
The VI plot identifies the features that are most important for explaining
different aspects of the conditional distribution; the spatial parameters are
found to be more important at extreme quantiles, while the process
realizations at the Vecchia neighbor locations are more important closer to
the median. The plots indicate that the conditional distributions at each
location are sensitive to the spatial parameters; details of the VI metric
used by SPQR are presented in the Supplementary Material (Majumder et al.,
2023, Appendix A.7). To assess goodness of fit, we repeat the process of
fitting the local SPQR model at site 25 for 100 independent datasets simulated
from the PMM. Figure 2(b) is a Q-Q plot based on true and fitted values from
the SPQR models, where each line corresponds to one of the 100 datasets. The
Q-Q plot is presented on the exponential scale (Heffernan and Tawn, 2001) to
verify whether the model can adequately capture tail behavior. The values fall
along the $Y=X$ line, suggesting a good model fit. Computation time for local
SPQR at sites with all 15 neighbors is approximately 22 minutes. The
doParallel package in R was used to parallelize SPQR model fits and improve
computation times.

Table 1: Coverage (in $\%$) for marginal GEV parameters under 2 scenarios based on MCMC simulations over 100 datasets. The 3 values represent the minimum, mean, and maximum coverage across the 50 study locations. | $\mu_{0}$ | $\mu_{1}$ | $\sigma$ | $\xi$
---|---|---|---|---
$\delta=0.2$ | (86, 92, 96) | (91, 96, 99) | (83, 91, 98) | (92, 96, 100)
$\delta=0.8$ | (85, 93, 100) | (90, 95, 100) | (86, 93, 98) | (90, 96, 100)
Figure 3: Sampling distribution of the posterior mean for the asymptotic
dependence parameter $\delta$ for two simulation scenarios. The horizontal
dashed lines are true values and the numbers along the bottom give the
empirical coverage of the 95% intervals.

We conduct a simulation study to explore how density-estimation errors
propagate to parameter-estimation errors for spatially varying GEV parameters.
Two scenarios are considered for the simulation studies, corresponding to
$\delta=\\{0.2,0.8\\}$. For both scenarios,  we assume the true values of
$\rho=0.15$ and $r=0.80$, and simulate 100 datasets at the 50 locations. Each
dataset consists of 50 independent (time) replicates. For spatial coordinates
$\mbox{\bf s}=(s_{1},s_{2})$, the marginal GEV parameters are:

$\displaystyle\mu_{0}(\mbox{\bf s})$ $\displaystyle=\exp\\{2+\cos(2\pi
s_{1})+\cos(2\pi s_{2})\\}$ $\displaystyle\mu_{1}(\mbox{\bf s})$
$\displaystyle=1$ $\displaystyle\sigma(\mbox{\bf s})$
$\displaystyle=\exp\\{\cos(2\pi s_{2})\\}$ $\displaystyle\xi(\mbox{\bf s})$
$\displaystyle=\frac{1}{2}\sin(\frac{\pi}{2}s_{1}).$

We model the marginal GEV parameters using an STVC model. At each location, we
assume the data to be GEV with parameters

$\displaystyle Y_{t}(\mbox{\bf s})\sim GEV\bigl{(}\mu_{0}(\mbox{\bf
s})+\mu_{1}(\mbox{\bf s})X^{*}_{t},\sigma(\mbox{\bf s}),\xi(\mbox{\bf
s})\bigr{)},$

with $X^{*}_{t}=(t-25.5)/10$, $t=1:50$. The variable $X_{t}^{*}$ represents
changes in the location parameter over time.

The intercept process $\mu_{0}(\mbox{\bf s})$ is assigned a GP prior with a
Matérn covariance function and nugget effects allow local heterogeneity:

$\displaystyle\mu_{0}(\mbox{\bf s})$ $\displaystyle=\tilde{\mu}_{0}(\mbox{\bf
s})+e_{0}(\mbox{\bf s})$ $\displaystyle e_{0}(\mbox{\bf s})$
$\displaystyle\stackrel{{\scriptstyle
iid}}{{\sim}}\mbox{Normal}(0,v_{\mu_{0}})$
$\displaystyle\tilde{\mu_{0}}(\mbox{\bf s})$
$\displaystyle\sim\mbox{GP}\bigl{(}\beta_{\mu_{0}},\tau^{2}_{\mu_{0}}K(\mbox{\bf
s},\mbox{\bf s}^{\prime};\rho_{\mu_{0}},\kappa_{\mu_{0}})\bigr{)}$
$\displaystyle\beta_{\mu_{0}}$
$\displaystyle\sim\mbox{Normal}(0,10^{2}),\tau^{2}_{\mu_{0}},v^{2}_{\mu_{0}}\stackrel{{\scriptstyle
iid}}{{\sim}}\mbox{IG}(0.1,0.1)$ $\displaystyle\log\rho_{\mu_{0}}$
$\displaystyle\sim\mbox{Normal}(-1,1),\log\kappa_{\mu_{0}}\sim\mbox{Normal}(-2,1),$

where $K(\mbox{\bf s},\mbox{\bf s}^{\prime};\rho_{\mu_{0}},\kappa_{\mu_{0}})$
is the Matérn correlation function with spatial range $\rho_{\mu_{0}}$ and
smoothness parameter $\kappa_{\mu_{0}}$, and IG$(\cdot,\cdot)$ is the inverse-
Gamma distribution. The slope $\mu_{1}(\mbox{\bf s})$, the log-scale
$\log\sigma(\mbox{\bf s})$, and the shape $\xi(\mbox{\bf s})$ are modeled
similarly using GPs. The spatial parameters have priors
$\delta\sim\mbox{Uniform}(0,1)$, $\rho\sim\mbox{Uniform}(0.0,0.5)$, and
$r\sim\mbox{Uniform(0,1)}$. Runtimes were approximately 1 minute for 1,000
iterations of the MCMC.

Table 1 details coverage of the empirical $95\%$ intervals for the posterior
distribution of the marginal GEV parameters. Mean coverage across locations is
near or at nominal level across different parameters and scenarios. Figure 3
plots the sampling distribution of the posterior mean estimator of $\delta$
for the 2 scenarios, and provides empirical coverage of the $95\%$ posterior
interval. Both scenarios show low bias; coverage is $91\%$ for the asymptotic
independence scenario, and $97\%$ for the asymptotic dependence scenario.
Overall, the SPQR approach is able to distinguish between the two asymptotic
regimes in the presence of spatially varying marginals.

Additional simulation studies are presented in the Supplementary Material
(Majumder et al., 2023, Appendix B) which involve special cases of the PMM.
The first considers a GP as the spatial process, and evaluates both the global
and local SPQR methods. The second study is also for a PMM but in a non-STVC
setting, where were explore a few additional scenarios including negative
shape parameters, and missing and censored data. The final demonstrates the
need of non-linear NN layers by fitting a SPQR model without any hidden NN
layers.

##  6 Analysis of Extreme Streamflow in the US

###  6.1 Data description and exploratory analysis

Figure 4: HCDN annual maxima: Sample 0.9 quantile of the log annual streamflow
maxima $Y_{t}(\mbox{\bf s})$ (in m${}^{3}/$s) at each of the 487 gauges.

We apply the proposed methods to model extreme streamflow from 1972–2021 at
487 stations across the US with complete data. These locations are part of the
USGS Hydro-Climatic Data Network (HCDN) 2009 (Lins, 2012) and are relatively
unaffected by human activities. The data is downloaded using the dataRetrieval
package in R (De Cicco et al., 2022), and the code is made available in our
GitHub repository111https://github.com/reetamm/SPQR-for-spatial-extremes. Our
goal is to identify regions within the US where the distribution of extreme
streamflow is changing over time. The annual maximum of daily streamflow is
measured in m${}^{3}/$s, and for each of the $T=50$ years and $n=487$
stations, the response $Y_{t}(\mbox{\bf s})$ is taken to be the logarithm of
the annual maximum. The log transformation was chosen as a Box-Cox
transformation parameter after comparing parameter values between -2 and 2 on
the basis of goodness of fit, profile likelihood values, and the stability of
initial MLE estimates at the locations. Figure 4 plots the sample 0.9 quantile
of the $T=50$ observations at each station, which show considerable spatial
variation.

(a) Empirical conditional exceedance $\chi_{u}(h)$ for log annual maximum
streamflow at different distances.

(b) Sample variogram for log annual maximum streamflow, averaged over 50 years
of data.

Figure 5: Spatial behavior of log annual maximum streamflow in terms of the
conditional exceedance and the variogram, units of km.

In order to study the dependence structure of the process, especially at its
extremes, we consider the conditional exceedance probability $\chi_{u}(h)$ of
maximum streamflow at pairs of stations separated by a distance $h=||\mbox{\bf
s}_{j}-\mbox{\bf s}_{k}||$ in kilometers (km). Figure 5(a) plots $\chi_{u}(h)$
for rank-standardized streamflow data as a function of $u$ for different
values of $h$. The rank standardization ensures a Uniform$(0,1)$ marginal
distribution at each location. Stations farther away from each other can be
seen to have less extremal dependence, with tail dependence approaching 0 for
stations 1000 km apart. Figure 5(b) plots the mean of the annual variograms of
the streamflow data. It shows a range of over 1500 km, as well as the presence
of a nugget effect. Both plots suggest that extremal streamflow is spatially
dependent at distances of 1000 km or more, even after accounting for the
spatial differences in the marginal distributions.

Figure 6: Sample variogram plots for the MLE of GEV parameters for extreme
streamflow at each location. Parameters are labeled on the top-left of each
panel.

For the marginals at each location, we assume GEV distributions with STVC,

$Y_{t}(\mbox{\bf s})\sim\mbox{GEV}\left[\mu_{0}(\mbox{\bf
s})+\mu_{1}(\mbox{\bf s})X_{t},\sigma(\mbox{\bf s}),\xi(\mbox{\bf s})\right],$
(9)

where $X_{t}=(\text{year}_{t}-1996.5)/10$ for $\text{year}_{t}=1972+t-1$. This
parameterization attempts to capture changes in the location parameter in the
past 50 years due to changing climate; positive values of $\mu_{1}(\mbox{\bf
s})$ would suggest an increase in the magnitude of the annual extremal
streamflow. Figure 6 plots variograms for MLE estimates (estimated separately
by location) of the GEV parameters at each location. All 4 GEV parameters show
spatial dependency, which motivates the STVC specification.

The marginal GEV parameters for each location are assigned GP priors with a
Matérn correlation functions and nugget effects to allow for local
heterogeneity. The prior specification is similar to the one used in Section
5, except with the additional assumption of a common Matérn smoothness
parameter was made for all four GP priors to improve MCMC convergence, i.e.,
$\kappa_{\mu_{0}}=\kappa_{\mu_{1}}=\kappa_{\sigma}=\kappa_{\xi}=\kappa$. For
the residual model, we use the PMM in Section 3 for spatial dependence and
assume independence across years. The specification includes a nugget based on
the mean variogram (Figure 5(b)). As priors for the joint parameters
$\mbox{\boldmath$\theta$}^{SPAT}$, we set $\delta,r\sim\mbox{Uniform}(0,1)$,
and $\rho\sim\mbox{Uniform}(0,6251)$ measured in km.

###  6.2 Results

The local SPQR approximation for the log of annual streamflow maxima is thus
based around $\mbox{\boldmath$\theta$}^{SPAT}=(\delta,\rho,r)$, and models are
fit using 200,000 synthetic observations at each of the 486 locations with
neighbors. Once the local SPQR models have been fit, we run two MCMC chains
for 30,000 iterations each, with two different starting values of $\delta$.
The first 10,000 iterations from each chain are discarded as burn-in;
additional results are provided in the Supplementary Material (Majumder et
al., 2023, Appendix C.2).

####  6.2.1 Parameter estimates:

The posterior means (standard deviations) of the spatial parameters are
${\hat{\delta}}=0.45\,(0.02)$, ${\hat{\rho}}=807\,(45)$ km, and
${\hat{r}}=0.92\,(0.004)$. The posterior of $\delta$ has a $95\%$ interval of
$(0.40,0.49)$ which puts the process in the asymptotic independence regime
with high probability. The GEV Matérn smoothness parameter estimate is
$\hat{\kappa}=0.60\,(0.03)$, and the four range parameters (in km) are
$\hat{\rho}_{\mu_{0}}=12435\,(10645)$, $\hat{\rho}_{\mu_{1}}=27605\,(10689)$,
$\hat{\rho}_{\sigma}=20311\,(11232)$, and $\hat{\rho}_{\xi}=20320\,(11481)$.
The STVC parameters are therefore much smoother over space than the year-to-
year variation captured by $\rho$. The intercept is the most variable GEV
parameter across space and the slope is the least, having the smallest and
largest range parameter estimates, respectively. It is likely that the
intercept varies the most because the magnitude of streamflow at a station is
dependent on very local features. Conversely, the slope parameter may vary
smoothly because the drivers of change are regional rather than local in
nature.

(a) Posterior mean of $\mu_{0}(\mbox{\bf s})$.

(b) Posterior mean of $\mu_{1}(\mbox{\bf s})$.

(c) Posterior mean of $\sigma(\mbox{\bf s})$.

(d) Posterior mean of $\xi(\mbox{\bf s})$.

Figure 7: HCDN GEV parameter estimates: Posterior means for 487 stations in
the USA based on log-transformed data from 1972–2021. Figure 8: Posterior of
$Pr[\mu_{1}(\mbox{\bf s})>0]$ for GEV location parameters based on log-
transformed data from 1972–2021.

Figure 7 plots the spatial distribution of posterior means for the GEV
parameters. The shape and scale parameters are negatively associated for large
parts of the country, possibly a consequence of the constraints on the GEV
parameters. The scale parameter $\sigma(\mbox{\bf s})$ is highest in the
Arkansas-Rio Grande-Texas Gulf and the Missouri basin regions222DOI unified
regions, and the scale parameter is highest in the North Atlantic-Appalachian
and the Columbia-Pacific Northwest regions. Areas with high estimates of
$\mu_{0}(\mbox{\bf s})$ and $\xi(\mbox{\bf s})$ coincide with areas of high
precipitation in the 1991–2020 US Climate Normals333US Climate Normals, US
Climate Atlas, the current official baseline for describing average US
climate. This suggests that an association between precipitation and
streamflow maxima. In Figure 7(d), we note that the posterior means of the
shape parameters $\xi(\mbox{\bf s})$ are negative at all 487 stations on the
log transformed data. The use of the log transform for the streamflow data
leads to negative GEV shape parameter estimates, imposing a finite upper bound
on the distribution even on the original scale. However, we do not expect this
to affect estimation of either the quantiles of the marginal distribution or
the joint exceedance probabilities. To assess the effect of transforming the
data, additional analysis was carried out on PMMs fitted to the original scale
of the data, as well as on the square-root of streamflow. Results from both
these analyses are provided in the Supplementary Material (Majumder et al.,
2023, Appendix C.4).

Our primary interest, however, are in estimates of the location parameters
across the USA. Figure 7(a) plots the posterior mean of the intercept
$\mu_{0}(\mbox{\bf s})$ of the location parameter, and has a spatial
distribution similar to Figure 4. Figure 7(b) plots the posterior mean of the
slope $\mu_{1}(\mbox{\bf s})$ of the location parameters with respect to time
across the USA, and Figure 8 plots the posterior probability of the slope
parameter being positive, $Pr[\mu_{1}(\mbox{\bf s})>0]$. Positive slope
estimates in Figure 7(b) indicate an increase in extreme streamflow over time,
and high probabilities in Figure 8 indicates stronger evidence for the
increase being significant. The majority of the positive slope parameters are
concentrated in the Mississippi and Missouri basins, and the Arkansas-Rio
Grande-Texas Gulf regions. The North Atlantic-Appalachian region in the east
has a large number of catchments with slope estimates near zero, but the
majority of zero and negative slope estimates are concentrated around the
Lower Colorado Basin, Columbia-Pacific Northwest, and the California-Great
Basin regions. An exception is Washington, in the northwest, which has high
estimates of the slope. The upper Colorado basin which includes Wyoming,
Colorado, and New Mexico is of particular interest since HCDN stations in the
region have relatively low 0.9 quantile values in Figure 4 as well as low
$\mu_{0}(\mbox{\bf s})$ estimates in Figure 7(a), suggesting that extreme
streamflow is starting to have large impacts in this region over time.

####  6.2.2 Regional joint exceedance behavior

Figure 9: PIT scores of marginal GEV fits for HCDN locations in Colorado
(left) and New Mexico (right).

To study the behavior of extreme streamflow jointly for multiple locations, we
consider two clusters of HCDN stations in Colorado (CO) and in New Mexico
(NM), comprising of 10 and 11 stations, respectively. Figure 9 plots the
probability integral transform (PIT) scores of marginal GEV fits for locations
within each cluster based on posterior means of GEV parameters, which suggest
adequate marginal fits for both sets of locations. To quantify the effect of
changing climate on extreme streamflow for each cluster, we look at the joint
posterior probability of streamflow maxima exceeding the observed 0.90
quantile values as shown in Figure 4, i.e., Pr$[Y_{t}(\mbox{\bf
s}_{i})>q_{i},i=1,\ldots,n_{i}]$ for the sample 0.90 quantile $q_{i}$ at
location $\mbox{\bf s}_{i}$, where $n_{i}$ is the number of stations within
each cluster. Since our marginal models have STVC, we are able to calculate
posterior probabilities for both 1972 and 2021. The probabilities are
calculated based on 200 post burn-in MCMC samples from the posterior
distribution of the parameters. For each MCMC sample, 20,000 observations are
generated from the fitted model; the 1972 and 2021 probabilities are
calculated based on 10,000 observations each.

In Colorado, the 10 HCDN stations correspond to catchments with drainage
ranging from $15.5\mbox{ km}^{2}$ to $432.9\mbox{ km}^{2}$. The cluster is
well separated from other stations, and have high posterior estimates of
$Pr[\mu_{1}(\mbox{\bf s})>0]$ (minimum = $0.10$, mean = $0.57$). The HCDN
stations are all situated in the Upper Colorado Basin region, and the set of
neighbors are spread across the Upper and Lower Colorado Basin regions. The
joint exceedance probability for 1972 has a mean of 0.075 with an SD of 0.04.
The joint exceedance probability for 2021 has a mean of 0.17 with an SD of
0.046. This corresponds to an increase of around $125\%$ of the mean joint
exceedance probability in the last 50 years. Note that if we assume
independence across locations, the joint exceedance probability for the 10
locations would be approximately $10^{-10}$. Finally, the probability that the
joint exceedance in 2021 is higher than in 1972 is 0.90, providing strong
evidence in favor of increased extremal streamflow in the area, possibly due
to changing climate.

In New Mexico, the 11 HCDN stations correspond to catchments with drainage
ranging from $43.8\mbox{ km}^{2}$ to $4804.9\mbox{ km}^{2}$. The stations are
all situated in the Upper Colorado Basin region, and the set of Vecchia
neighbors are located across the Upper and Lower Colorado Basin as well as the
California-Great Basin region. The catchments have a mix of high and low
posterior estimates of $Pr[\mu_{1}(\mbox{\bf s})>0]$ (minimum = 0.07, mean =
0.48). The joint exceedance probability for 1972 has a mean of 0.045 with an
SD of 0.012. The joint exceedance probability for 2021 has a mean of 0.053
with an SD of 0.017. This corresponds to an increase of around $18\%$ of the
mean joint exceedance probability in the last 50 years. The probability that
the the joint exceedance in 2021 is higher than in 1972 is 0.695, which is
lower than the result for Colorado, but still significantly higher than the
independent scenario. Majumder et al. (2023, Appendix C.2) contains additional
results for both clusters.

####  6.2.3 Model comparison and model fit

Table 2: Estimates and standard errors (in parentheses) from leave-one-out cross validation (LOO-CV) and the Watanabe-Akaike information criterion (WAIC) for comparing the process mixture model (PMM), the Huser-Wadsworth model (HW), the max-stable process model (MSP), and the Gaussian process model (GP). Lower values indicate a better fit. | PMM | HW | MSP | GP
---|---|---|---|---
LOO-CV | 29108 (540) | 29708 (544) | 32058 (583) | 33842 (561)
WAIC | 29559 (549) | 30193 (565) | 33441 (552) | 34440 (585)

We compared the PMM with three spatial processes - a GP, an MSP, and the
Huser-Wadsworth (HW) process. In each case, inference is carried out by using
SPQR for density estimation, and using MCMC afterwards for parameter
estimation. This allows for a comparison of the appropriateness of different
spatial processes for our application. For all three competing models, we used
the same neural network architecture in our local SPQR fits as the PMM.

Table 2 lists estimates and standard errors from leave-one-out cross
validation (LOO-CV) (Vehtari et al., 2016) and the Watanabe-Akaike information
criterion (WAIC) (Watanabe, 2010) for the four models. Based on both metrics,
the PMM has the lowest values and thus the best model fit, followed by the HW
model. The MSP and the GP are noticeably worse than the PMM and HW models, but
their LOO-CV and WAIC estimates are relatively close to each other. This
suggests that models which assume asymptotic (in)dependence (e.g., the MSP and
the GP) are likely to have worse fits when the dependence structure is
complex. An inter-comparison of the posterior means of the slope parameter
from each of these models is presented in the Supplementary Material (Majumder
et al., 2023, Appendix C.3).

Figure 10: Estimates of $\chi_{u}(h)$ for three different regions: Empirical
estimates of the upper tail coefficient in bold, compared against estimates
based on the posterior distribution. The dashed lines and the bands correspond
to the mean and $95\%$ interval, respectively.

Finally, we evaluate the fit when using the PMM as the underlying spatial
process by comparing estimates of $\chi_{u}(h)$ based on the posterior
distribution with empirical estimates obtained from the rank-standardized
observed data. For this, we choose 3 regions within the USA, the first two of
which are the CO and NM clusters studied in Section 6.2.2. Both clusters have
a total of 34 locations (including the Vecchia neighbor locations). The third
area consists of HCDN locations in Oregon (OR) and with its Vecchia neighbors,
accounting for a total of 56 locations in Oregon, Washington, Nevada, and
California. This region provides a contrast to the CO and NM clusters;
locations are farther apart, and have low estimates of the slope and scale but
high estimates of the shape parameter. For each of the 3 clusters, we
generated realizations of the spatial process based 200 samples from the
posterior, and estimated $\chi_{u}(h)$ for each sample.

Figure 10 plots the empirical and posterior estimates of $\chi_{u}(h)$ for the
three regions. The top row plots the upper tail coefficient for high
quantiles, for locations $h=10\mbox{ km}$ apart for CO and NM, and
$h=100\mbox{ km}$ apart for OR. The bottom row plots estimates of
$\chi_{u}(h)$ for $u=0.90$. In each panel, the bold line represents the
empirical estimate, and the dashed line and band represent the mean and $95\%$
interval of the posterior estimates respectively. We note that both the model-
based and empirical estimates of the upper tail coefficient behaves similarly
for the three regions, and are consistent with the empirical estimate that is
obtained from data for the entire USA (Figure 5(a)). This suggests that the
stationarity assumption for the dependence structure in our model is
appropriate. The $95\%$ intervals in the top row includes 0 for all 3 regions,
reflecting the asymptotic independence corresponding to our estimate of
$\delta$. Both the empirical and posterior estimates show behavior similar to
Figure 1(a) for the asymptotic independent case. The tail coefficients can
also been seen to decrease as distance increases, similar to Figure 1(b). We
noted similar behavior for other regions within the country. While the fitted
model tends to overestimate $\chi_{u}(h)$ for several regions, it is able to
capture the behavior of extremal dependence over large distances, and at high
quantile levels.

##  7 Discussion

In this paper, we proposed a process mixture model (PMM) for spatial extremes,
where the marginal distributions at different spatial locations are GEV and
their spatial dependence is captured using a convex combination of a GP and an
MSP. The PMM extends Huser and Wadsworth (2019), and is flexible enough to
accommodate missingness and censoring, as well as STVC for the marginal GEV
distributions. We approximated the joint likelihood for the spatial model
using a Vecchia approximation. We used the density regression model proposed
in Xu and Reich (2021) to approximate this likelihood, whose weights are
modeled using a feed-forward neural network and learned using synthetic data
generated from a design distribution. Parameter estimation for the model is
carried out using MCMC.

We used the PMM to analyze changes in annual maximum streamflow within the US
over the past 50 years. For this study, we used the annual maximum streamflow
measured at 487 stations in the USGS Hydro-Climatic Data Network. The
posterior means of the location parameter have non-zero slope estimates in
several parts of the country. We noted a high concentration of positive slope
estimates in the Mississipi and Missouri basins, and the Arkansas-Rio Grande-
Texas Gulf regions, indicating that extremal streamflow has increased in those
areas over the last 50 years.

Future work will focus on the theoretical properties of this model. While it
is straightforward to derive an analytical expression for $\chi_{u}(h)$ for
the trivial case of a shared MSP and has been provided in the Supplementary
Material (Majumder et al., 2023, Appendix A.4), i.e., for $R(\mbox{\bf s})=R$,
obtaining an analytical expression for the general case is more challenging.
This would also enable us to investigate the properties of $\chi_{u}(h)$ for
the PMM for different values of $\delta$ and as $h\to\infty$. We would like to
further investigate improvements to the computational aspects of this model
and identify reasonable plug-and-play settings for local and global SPQR
approximations. Finally, we would also like to extend this model to provide
climate-informed estimates by regressing the spatiotemporal variability of the
EVA parameters onto large-scale climate drivers from GCM output with
spatially-varying regression coefficients for local calibration. Another area
of future work is to extend the model to accommodate more complex spatial
dependence structures. Recent work on spatial extremes has incorporated
graphical models as additional information for computing distances (Engelke
and Hitz, 2020). While our work considers only spatial coordinates, it is
possible to incorporate additional distance measures in the form of river
network information to reflect the physical structure of watersheds and the
so-called ‘river distance’ between stations (Asadi et al., 2015). While the
stream network information is not readily available for these data,
incorporating this network structure might improve spatial modeling. Spatial
models on stream networks have been developed for both max-stable (Asadi et
al., 2015) and Gaussian (Santos-Fernandez et al., 2022) processes, and so it
should be possible to incorporate these features into the PMM.

## Acknowledgements

This work was supported by grants from the Southeast National Synthesis
Wildfire and the United States Geological Survey’s National Climate Adaptation
Science Center (G21AC10045), the National Science Foundation (CBET2151651,
DMS2152887, and DMS2001433) and the National Institutes of Health
(R01ES031651-01). The authors thank Prof. Sankarasubramanian Arumugam of North
Carolina State University for discussion of the data and scope of the project.

## References

  * Abrahamowicz et al. (1992) Abrahamowicz, M., Clampl, A. and Ramsay, J. O. (1992) Nonparametric density estimation for censored survival data: Regression-spline approach.  Canadian Journal of Statistics, 20, 171–185.
  * Archfield et al. (2016) Archfield, S. A., Hirsch, R. M., Viglione, A. and Blöschl, G. (2016) Fragmented patterns of flood change across the United States.  Geophysical research letters, 43, 10–232.
  * Asadi et al. (2015) Asadi, P., Davison, A. C. and Engelke, S. (2015) Extremes on river networks.  The Annals of Applied Statistics, 9, 2023–2050.
  * Blöschl et al. (2019) Blöschl, G., Hall, J., Viglione, A., Perdigão, R. A., Parajka, J., Merz, B., Lun, D., Arheimer, B., Aronica, G. T., Bilibashi, A., Bohac, M., Bonacci, O., Borga, M., Čanjevac, I., Castellarin, A., Chirico, G., Claps, P., Frolova, N., Ganora, D., Gorbachova, L., Gul, A., Hannaford, J., Harrigan, S., Kireeva, M., Kiss, A., Kjeldsen, T., Kohnova, S., Koskela, J., Ledvinka, O., Macdonald, N., Mavrova-Guirguinova, M., Mediero, L., Merz, R., Molnar, P., Montanari, A., Murphy, C., Osuch, M., Ovcharuk, V., Radevski, I., Salinas, J., Sauquet, E., Sraj, M., Szolgay, J., Volpi, E., Wilson, D., Zaimi, K. and Zivkovic, N. (2019) Changing climate both increases and decreases European river floods.  Nature, 573, 108–111.
  * Bopp et al. (2021) Bopp, G. P., Shaby, B. A. and Huser, R. (2021) A hierarchical max-infinitely divisible spatial model for extreme precipitation.  Journal of the American Statistical Association, 116, 93–106.
  * Brown and Resnick (1977) Brown, B. M. and Resnick, S. I. (1977) Extreme values of independent stochastic processes.  Journal of Applied Probability, 14, 732–739.
  * Castruccio et al. (2016) Castruccio, S., Huser, R. and Genton, M. G. (2016) High-order composite likelihood inference for max-stable distributions and processes.  Journal of Computational and Graphical Statistics, 25, 1212–1229.
  * Chui et al. (1980) Chui, C., Smith, P. and Ward, J. (1980) Degree of $L_{p}$ Approximation by Monotone Splines.  SIAM Journal on Mathematical Analysis, 11, 436–447.
  * Coles et al. (2001) Coles, S., Bawa, J., Trenner, L. and Dorazio, P. (2001) An introduction to statistical modeling of extreme values, vol. 208.  Springer.
  * Condon et al. (2015) Condon, L., Gangopadhyay, S. and Pruitt, T. (2015) Climate change and non-stationary flood risk for the upper Truckee River basin.  Hydrology and Earth System Sciences, 19, 159–175.
  * Datta et al. (2016) Datta, A., Banerjee, S., Finley, A. O. and Gelfand, A. E. (2016) Hierarchical nearest-neighbor Gaussian process models for large geostatistical datasets.  Journal of the American Statistical Association, 111, 800–812.
  * Dawdy et al. (2012) Dawdy, D. R., Griffis, V. W. and Gupta, V. K. (2012) Regional flood-frequency analysis: How we got here and where we are going.  Journal of Hydrologic Engineering, 17, 953–959.
  * De Cicco et al. (2022) De Cicco, L. A., Lorenz, D., Hirsch, R. M., Watkins, W. and Johnson, M. (2022) dataRetrieval: R packages for discovering and retrieving water data available from U.S. federal hydrologic web services.  Reston, VA.  URLhttps://code.usgs.gov/water/dataRetrieval.
  * De Haan et al. (2006) De Haan, L., Ferreira, A. and Ferreira, A. (2006) Extreme value theory: an introduction, vol. 21.  Springer.
  * Engelke and Hitz (2020) Engelke, S. and Hitz, A. S. (2020) Graphical Models for Extremes.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 82, 871–932.
  * Erhardt and Smith (2012) Erhardt, R. J. and Smith, R. L. (2012) Approximate Bayesian computing for spatial extremes.  Computational Statistics & Data Analysis, 56, 1468–1481.
  * François et al. (2019) François, B., Schlef, K., Wi, S. and Brown, C. (2019) Design considerations for riverine floods in a changing climate–a review.  Journal of Hydrology, 574, 557–573.
  * Franks (2002) Franks, S. W. (2002) Identification of a change in climate state using regional flood data.  Hydrology and Earth System Sciences, 6, 11–16.
  * Gerber and Nychka (2021) Gerber, F. and Nychka, D. (2021) Fast covariance parameter estimation of spatial Gaussian process models using neural networks.  Stat, 10, e382.
  * Greenberg et al. (2019) Greenberg, D., Nonnenmacher, M. and Macke, J. (2019) Automatic posterior transformation for likelihood-free inference.  In Proceedings of the 36th International Conference on Machine Learning (eds. K. Chaudhuri and R. Salakhutdinov), vol. 97 of Proceedings of Machine Learning Research, 2404–2414. PMLR.
  * Guinness (2018) Guinness, J. (2018) Permutation and grouping methods for sharpening gaussian process approximations.  Technometrics, 60, 415–429.
  * de Haan (1984) de Haan, L. (1984) A spectral representation for max-stable processes.  The Annals of Probability, 12, 1194–1204.
  * de Haan and Ferreira (2006) de Haan, L. and Ferreira, A. (2006) Extreme value theory: An introduction.  Springer.
  * Hazra et al. (2021) Hazra, A., Huser, R. and Bolin, D. (2021) Realistic and fast modeling of spatial extremes over large geographical domains.  URLhttps://arxiv.org/abs/2112.10248.
  * Heffernan and Tawn (2001) Heffernan, J. E. and Tawn, J. A. (2001) Extreme value analysis of a large designed experiment: A case study in bulk carrier safety.  Extremes, 4, 359–378.
  * Hirabayashi et al. (2013) Hirabayashi, Y., Mahendran, R., Koirala, S., Konoshima, L., Yamazaki, D., Watanabe, S., Kim, H. and Kanae, S. (2013) Global flood risk under climate change.  Nature Climate Change, 3, 816–821.
  * Hirsch (2011) Hirsch, R. M. (2011) A perspective on nonstationarity and water management 1.  JAWRA Journal of the American Water Resources Association, 47, 436–446.
  * Hirsch and Ryberg (2012) Hirsch, R. M. and Ryberg, K. R. (2012) Has the magnitude of floods across the USA changed with global CO2 levels?  Hydrological Sciences Journal, 57, 1–9.
  * Hornik et al. (1989) Hornik, K., Stinchcombe, M. and White, H. (1989) Multilayer feedforward networks are universal approximators.  Neural Networks, 2, 359–366.
  * Huser and Davison (2014) Huser, R. and Davison, A. C. (2014) Space–time modelling of extreme events.  Journal of the Royal Statistical Society, Series B, 76, 439–461.
  * Huser et al. (2016) Huser, R., Davison, A. C. and Genton, M. G. (2016) Likelihood estimators for multivariate extremes.  Extremes, 19, 79–103.
  * Huser et al. (2019) Huser, R., Dombry, C., Ribatet, M. and Genton, M. G. (2019) Full likelihood inference for max-stable data.  Stat, 8, e218.
  * Huser et al. (2022) Huser, R., Stein, M. L. and Zhong, P. (2022) Vecchia likelihood approximation for accurate and fast inference in intractable spatial extremes models.  arXiv preprint arXiv:2203.05626.
  * Huser and Wadsworth (2019) Huser, R. and Wadsworth, J. L. (2019) Modeling spatial processes with unknown extremal dependence class.  Journal of the American Statistical Association, 114, 434–444.
  * Järvenpää et al. (2021) Järvenpää, M., Gutmann, M. U., Vehtari, A. and Marttinen, P. (2021) Parallel Gaussian process surrogate Bayesian inference with noisy likelihood evaluations.  Bayesian Analysis, 16, 147–178.
  * Joe (1997) Joe, H. (1997) Multivariate Models and Multivariate Dependence Concepts.  Chapman and Hall/CRC, 1st edn.
  * Kabluchko et al. (2009) Kabluchko, Z., Schlather, M. and de Haan, L. (2009) Stationary max-stable fields associated to negative definite functions.  The Annals of Probability, 37, 2042–2065.
  * Katzfuss and Guinness (2021) Katzfuss, M. and Guinness, J. (2021) A general framework for Vecchia approximations of Gaussian processes.  Statistical Science, 36, 124–141.
  * Kingma and Ba (2014) Kingma, D. P. and Ba, J. (2014) Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * Knox (1993) Knox, J. C. (1993) Large increases in flood magnitude in response to modest changes in climate.  Nature, 361, 430–432.
  * Kobyzev et al. (2021) Kobyzev, I., Prince, S. J. and Brubaker, M. A. (2021) Normalizing flows: An introduction and review of current methods.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 43, 3964–3979.
  * Kundzewicz et al. (2014) Kundzewicz, Z. W., Kanae, S., Seneviratne, S. I., Handmer, J., Nicholls, N., Peduzzi, P., Mechler, R., Bouwer, L. M., Arnell, N., Mach, K., Muir-Wood, R., Robert Brakenridge, G., Kron, W., Benito, G., Honda, Y., Takahashi, K. and Sherstyukov, B. (2014) Flood risk and climate change: global and regional perspectives.  Hydrological Sciences Journal, 59, 1–28.
  * Kundzewicz et al. (2017) Kundzewicz, Z. W., Krysanova, V., Dankers, R., Hirabayashi, Y., Kanae, S., Hattermann, F. F., Huang, S., Milly, P. C., Stoffel, M., Driessen, P., , Matczak, P., Quevauviller, P. and Schellnhuber, H. (2017) Differences in flood hazard projections in europe–their causes and consequences for decision making.  Hydrological Sciences Journal, 62, 1–14.
  * Kunkel et al. (2020) Kunkel, K. E., Karl, T. R., Squires, M. F., Yin, X., Stegall, S. T. and Easterling, D. R. (2020) Precipitation extremes: Trends and relationships with average precipitation and precipitable water in the contiguous United States.  Journal of Applied Meteorology and Climatology, 59, 125–142.
  * Ledford and Tawn (1996) Ledford, A. W. and Tawn, J. A. (1996) Satistics for near independence in multivariate extreme values.  Biometrika, 83, 169–187.
  * Ledford and Tawn (1997) — (1997) Modelling dependence within joint tail regions.  Journal of the Royal Statistical Society. Series B (Methodological), 59, 475–499.
  * Lenzi et al. (2021) Lenzi, A., Bessac, J., Rudi, J. and Stein, M. L. (2021) Neural networks for parameter estimation in intractable models.  arXiv preprint arXiv:2107.14346.
  * Li et al. (2019) Li, L., Holbrook, A., Shahbaba, B. and Baldi, P. (2019) Neural network gradient Hamiltonian Monte Carlo.  Computational statistics, 34, 281–299.
  * Lima et al. (2016) Lima, C. H., Lall, U., Troy, T. and Devineni, N. (2016) A hierarchical bayesian gev model for improving local and regional flood quantile estimates.  Journal of Hydrology, 541, 816–823.
  * Lins (2012) Lins, H. F. (2012) Usgs hydro-climatic data network 2009 (hcdn-2009).  US Geological Survey Fact Sheet, 3047.
  * Majumder et al. (2023) Majumder, R., Reich, B. J. and Shaby, B. A. (2023) Supplement to “Modeling extremal streamflow using deep learning approximations and a flexible spatial process”.
  * Meehl et al. (2000) Meehl, G. A., Zwiers, F., Evans, J., Knutson, T., Mearns, L. and Whetton, P. (2000) Trends in extreme weather and climate events: issues related to modeling extremes in projections of future climate change.  Bulletin of the American Meteorological Society, 81, 427–436.
  * Merz et al. (2014) Merz, B., Aerts, J., Arnbjerg-Nielsen, K., Baldi, M., Becker, A., Bichet, A., Blöschl, G., Bouwer, L. M., Brauer, A., Cioffi, F., Delgado, J., Gocht, M., Guzzetti, F., Harrigan, S., Hirschboeck, K., Kilsby, C., Kron, W., Kwon, H., Lall, U., Merz, R., Nissen, K., Salvatti, P., Swierczynski, T., Ulbrich, U., Viglione, A., Ward, P., Weiler, M., Wilhelm, B. and Nied, M. (2014) Floods and climate: emerging perspectives for flood risk assessment and management.  Natural Hazards and Earth System Sciences, 14, 1921–1942.
  * Milly et al. (2008) Milly, P., Betancourt, J., Falkenmark, M., Hirsch, R. M., Kundzewicz, Z. W., Lettenmaier, D. P. and Stouffer, R. J. (2008) Stationarity is dead: Whither water management?  Earth, 4.
  * Milly et al. (2015) Milly, P. C., Betancourt, J., Falkenmark, M., Hirsch, R. M., Kundzewicz, Z. W., Lettenmaier, D. P., Stouffer, R. J., Dettinger, M. D. and Krysanova, V. (2015) On critiques of “Stationarity is dead: Whither water management?”.  Water Resources Research, 51, 7785–7789.
  * Milly et al. (2005) Milly, P. C., Dunne, K. A. and Vecchia, A. V. (2005) Global pattern of trends in streamflow and water availability in a changing climate.  Nature, 438, 347–350.
  * Morris et al. (2019) Morris, S. A., Reich, B. J. and Thibaud, E. (2019) Exploration and inference in spatial extremes using empirical basis functions.  Journal of Agricultural, Biological and Environmental Statistics, 24, 555–572.
  * Nair and Hinton (2010) Nair, V. and Hinton, G. E. (2010) Rectified linear units improve restricted boltzmann machines.  In Proceedings of the 27th International Conference on International Conference on Machine Learning, ICML’10, 807–814. Madison, WI, USA: Omnipress.
  * Padoan et al. (2010) Padoan, S., Ribatet, M. and Sisson, S. (2010) Likelihood-based inference for max-stable processes.  Journal of the American Statistical Association, 105, 263–277.
  * Papamakarios et al. (2021) Papamakarios, G., Nalisnick, E., Rezende, D. J., Mohamed, S. and Lakshminarayanan, B. (2021) Normalizing flows for probabilistic modeling and inference.  Journal of Machine Learning Research, 22, 1–64.  URLhttp://jmlr.org/papers/v22/19-1028.html.
  * Penrose (1992) Penrose, M. D. (1992) Semi-Min-Stable Processes.  The Annals of Probability, 20, 1450 – 1463.
  * Price et al. (2018) Price, L. F., Drovandi, C. C., Lee, A. and Nott, D. J. (2018) Bayesian synthetic likelihood.  Journal of Computational and Graphical Statistics, 27, 1–11.
  * Rasmussen (2003) Rasmussen, C. E. (2003) Gaussian processes to speed up hybrid Monte Carlo for expensive Bayesian integrals.  In Seventh Valencia international meeting, dedicated to Dennis V. Lindley, 651–659. Oxford University Press.
  * Reich and Shaby (2012) Reich, B. J. and Shaby, B. A. (2012) A hierarchical max-stable spatial model for extreme precipitation.  The Annals of Applied Statistics, 6, 1430–1451.
  * Reich et al. (2013) Reich, B. J., Shaby, B. A. and Cooley, D. (2013) A hierarchical model for serially-dependent extremes: A study of heat waves in the western US.  Journal of Agricultural, Biological and Environmental Statistics, 19, 119–135.
  * Ribatet et al. (2012) Ribatet, M., Cooley, D. and Davison, A. C. (2012) Bayesian inference from composite likelihoods, with an application to spatial extremes.  Statistica Sinica, 22, 813–845.
  * Sainsbury-Dale et al. (2023) Sainsbury-Dale, M., Zammit-Mangion, A. and Huser, R. (2023) Fast optimal estimation with intractable models using permutation-invariant neural networks.
  * Salas and Obeysekera (2014) Salas, J. D. and Obeysekera, J. (2014) Revisiting the concepts of return period and risk for nonstationary hydrologic extreme events.  Journal of Hydrologic Engineering, 19, 554–568.
  * Sang and Genton (2014) Sang, H. and Genton, M. G. (2014) Tapered composite likelihood for spatial max-stable models.  Spatial Statistics, 8, 86–103.
  * Santos-Fernandez et al. (2022) Santos-Fernandez, E., Ver Hoef, J. M., Peterson, E. E., McGree, J., Isaak, D. J. and Mengersen, K. (2022) Bayesian spatio-temporal models for stream networks.  Computational Statistics & Data Analysis, 170, 107446\.
  * Schlather (2002) Schlather, M. (2002) Models for stationary max-stable random fields.  Extremes, 5, 33–44.
  * Sharma et al. (2018) Sharma, A., Wasko, C. and Lettenmaier, D. P. (2018) If precipitation extremes are increasing, why aren’t floods?  Water Resources Research, 54, 8545–8551.
  * Smith (1990) Smith, R. L. (1990) Max-stable processes and spatial extremes.  Unpublished manuscript.
  * Šraj et al. (2016) Šraj, M., Viglione, A., Parajka, J. and Blöschl, G. (2016) The influence of non-stationarity in extreme hydrological events on flood frequency estimation.  J. Hydrol. Hydromech, 64, 426–437.
  * Stein et al. (2004) Stein, M. L., Chi, Z. and Welty, L. J. (2004) Approximating likelihoods for large spatial data sets.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 66, 275–296.
  * Tawn (1990) Tawn, J. A. (1990) Modelling multivariate extreme value distributions.  Biometrika, 77, 245–253.
  * Vecchia (1988) Vecchia, A. V. (1988) Estimation and model identification for continuous spatial processes.  Journal of the Royal Statistical Society: Series B (Methodological), 50, 297–312.
  * Vehtari et al. (2016) Vehtari, A., Gelman, A. and Gabry, J. (2016) Practical bayesian model evaluation using leave-one-out cross-validation and WAIC.  Statistics and Computing, 27, 1413–1432.
  * Vogel et al. (2011) Vogel, R. M., Yaindl, C. and Walter, M. (2011) Nonstationarity: flood magnification and recurrence reduction factors in the United States.  JAWRA Journal of the American Water Resources Association, 47, 464–474.
  * Wadsworth (2015) Wadsworth, J. L. (2015) On the occurrence times of componentwise maxima and bias in likelihood inference for multivariate max-stable distributions.  Biometrika, 102, 705–711.  URLhttps://doi.org/10.1093/biomet/asv029.
  * Wadsworth and Tawn (2012) Wadsworth, J. L. and Tawn, J. A. (2012) Dependence modelling for spatial extremes.  Biometrika, 99, 253–272.
  * Wadsworth and Tawn (2014) — (2014) Efficient inference for spatial extreme value processes associated to log-Gaussian random functions.  Biometrika, 101, 1–15.  URLhttps://doi.org/10.1093/biomet/ast042.
  * Walter (2010) Walter, M. (2010) Increasing trends in peak flows in the northeastern United States and their impacts on design.  Ph.D. thesis, Tufts University.
  * Wang and Li (2018) Wang, H. and Li, J. (2018) Adaptive Gaussian process approximation for Bayesian inference with expensive likelihood functions.  Neural Computation, 30, 3072–3094.
  * Wang and Stoev (2010) Wang, Y. and Stoev, S. A. (2010) On the structure and representations of max-stable processes.  Advances in Applied Probability, 42, 855–877.
  * Watanabe (2010) Watanabe, S. (2010) Asymptotic equivalence of bayes cross validation and widely applicable information criterion in singular learning theory.  J. Mach. Learn. Res., 11, 3571–3594.
  * Wilkinson (2014) Wilkinson, R. (2014) Accelerating ABC methods using Gaussian processes.  In Artificial Intelligence and Statistics, 1015–1023. PMLR.
  * Winsemius et al. (2018) Winsemius, H. C., Jongman, B., Veldkamp, T. I., Hallegatte, S., Bangalore, M. and Ward, P. J. (2018) Disaster risk, climate change, and poverty: assessing the global exposure of poor people to floods and droughts.  Environment and Development Economics, 23, 328–348.
  * Xu and Majumder (2022) Xu, S. and Majumder, R. (2022) SPQR: Semi-Parametric Quantile Regression.  URLhttps://CRAN.R-project.org/package=SPQR.  R package version 0.1.0.
  * Xu and Reich (2021) Xu, S. G. and Reich, B. J. (2021) Bayesian nonparametric quantile process regression and estimation of marginal quantile effects.  Biometrics, 00, 1–14.  URLhttps://onlinelibrary.wiley.com/doi/abs/10.1111/biom.13576.