{'fraction_non_alphanumeric': 0.04068208895949166, 'fraction_numerical': 0.02256254964257347, 'mean_word_length': 4.564779005524862, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modeling Extremal Streamflow using Deep Learning Approximations and a Flexible Spatial Process

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 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 process mixture model 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.


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., 2008Milly et al., , 2015Vogel 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 extremeal 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, i.e., 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 borrowing 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.
The primary objective of our study is to identify regions within the US where changing climate conditions have affected the behaviour of extreme streamflow. We consider locations in the United States Geological Survey's Hydro-Climatic Data Network (HCDN-2009); the HCDN is screened to exclude sites where human activities affect the natural flow of the watercourse. The network locations have a long historical record, making it suitable for analyzing hydrologic variations and trends in a climatic context (Lins, 2012). 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;Buishand et al., 2008;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., Huser and Davison, 2014;Reich et al., 2013).
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 andTawn, 2012, 2014;Wadsworth, 2015a). For general MSPs, Castruccio et al. (2016) stated that full likelihood inference seemed limited to n = 13. 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 , finite-sample bias when using all pairs of observations (Sang and Genton, 2014;Wadsworth, 2015b;Castruccio et al., 2016), and computational challenges posed by computing likelihoods at all O(n 2 ) pairs. More recently,  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. Deep learning has been used to estimate parameters in spatial models (Gerber and Nychka, 2021;Lenzi et al., 2021) 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. However, it is difficult to extend them to problems with many parameters such as the spatio-temporally varying coefficients (STVC) case that are crucial to our application. 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). For lower dimensional problems, approximate Bayesian computation (ABC) can replace likelihood evaluation with repeated simulation from the MSP model (Erhardt and Smith, 2012).
A second challenge posed by MSPs is that they are restrictive in the class of dependence types they can incorporate. Environmental data often has weakening spatial dependence with increasing levels of extremeness; 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 with 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. From a modeling perspective, the mixture of the two spatial processes allows asymptotic dependence or independence for locations separated by distance h, with asymptotic independence as h → ∞. 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 process mixture model to analyze changes in annual maximum streamflow within the US over the past 50 years.

A Process Mixture Model for Spatial Extremes
Let Y(s) be the extreme observation at spatial location s. In this section we consider Y(s) to be defined as a block maximum, but the methods can be extended to peaks over a threshold (scenario 5 in Section 4.2). We assume a potentially different marginal distribution for each spatial location s and denote F s as the marginal cumulative distribution function (CDF) for site s. For example, we assume that F s is the generalized extreme value (GEV) distribution with location µ(s), scale σ(s) and shape ξ(s) so that marginally Y(s) ∼ GEV{µ(s), σ(s), ξ(s)}.

Then the transformed variables
share common uniform marginal distributions across the spatial domain. This transformation separates residual spatial dependence in U(s) from the spatial dependence induced by spatial variation in the GEV parameters, which we model using GP priors over s.
We define our spatial dependence model on U(s) by taking U(s) = G{V(s)}, such that where R(s) is a max-stable process (MSP), W(s) is a Gaussian process (GP), g R and g W are transformations that guarantee g R {R(s)} and g W {W(s)} follow the standard exponential distribution, and δ ∈ [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 Huser and Wadsworth (2019), who assumed a standard Pareto random variable R common to all locations, by replacing it with an MSP. Since (2) mixes two processes, we refer to it as the process mixture model.
By construction, V(s) marginally follows the two-parameter hypoexponential distribution with CDF Without loss of generality, we assume that R(s) has a marginal GEV(1,1,1) distribution and W(s) has a standard normal marginal distribution. In this case, the transformations Although other options are possible, we model the correlation of the GP W(s) using the isotropic powered-exponential correlation function Cor{W(s 1 , s 2 )} = exp{−(h/ρ W ) α W } with distance h = ||s 1 − s 2 ||, smoothness α W ∈ (0, 2), and range ρ W > 0. The MSP R(s) is assumed to have isotropic Brown-Resnick spatial dependence defined by the variogram γ(h) = (h/ρ R ) α R for smoothness α R ∈ (0, 2) and range ρ R > 0.
Extremal spatial dependence of the process at locations s 1 and s 2 is often measured using the conditional exceedance probability where u ∈ (0, 1) is a threshold. The random variables U(s 1 ) and U(s 2 ) are defined as asymptotically dependent if the limit is positive and independent if χ(s 1 , s 2 ) = 0. Since we assume both R(s) and W(s) are  isotropic processes, we simply write χ u (h) and χ(h) as a function of the distance between locations.

Mixture Model
Fitting the process mixture model introduced in Section 2 poses computation challenges, especially for large datasets. The joint distribution for the GP W(s) is available in closed form but is cumbersome for large datasets; the joint distribution of the MSP R(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;Jabot et al., 2014;Wilkinson, 2014;Gutmann and Corander, 2016;Price et al., 2018;Drovandi et al., 2018;Wang and Li, 2018;Acerbi, 2018;Järvenpää et al., 2019Järvenpää et al., , 2021Li 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 s 1 , ..., s n . Partition the parameters into those that affect the marginal distributions, denoted θ 1 , and those that affect the spatial dependence, denoted θ 2 . For the model in Section 2, θ 1 includes the GEV parameters does not depend on θ 1 . 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 ; θ 1 , We approximate the joint likelihood in (6) using a Vecchia approximation (Vecchia, 1988;Stein et al., 2004;Datta et al., 2016;Katzfuss and Guinness, 2021), for u (i) = {u j ; j ∈ N i } and neighboring set N i ⊆ {1, ..., i − 1}, e.g., the k locations in N i that are closest to s i . Here, we use the notation that the collection of variables z i over the neighboring set is denoted z (i) = {z j ; j ∈ N i }. Of course, not all locations that are dependent with location i need be included in N i because distant observations may be approximately independent after conditioning on more local observations.
The conditional distributions for the process mixture model do not have closed-form expressions. We consider two related approximations -the local approximation where the n conditional density functions are fit separately, and a global approximation where a single density estimate is fit and applied to all density functions. Both approaches are described in further detail at the end of this section. In the global approximation, we fit a density regression viewing s (i) , u (i) and θ 2 as features (covariates), denoted as x i . Since the process is assumed to be stationary in space, only the differences s j − s i influence the regression model. To ensure that x i has the same length for all sites, we inflate the feature vector for sites i ∈ {1, . . . , K − 1} with large values of s j − s i and u j ∼ Uniform(0, 1). For the local approximation, u (i) and θ 2 make up the feature set.
For density regression, we use the model introduced in Xu and Reich (2021): where π k (x, W ) ≥ 0 are probability weights with ∑ K k=1 π k (x) = 1 that depend on the parameters W and B k (u) ≥ 0 are M-spline basis functions that, by definition, satisfy B k (u)du = 1 for all k. By increasing the number of basis functions K and appropriately selecting the weights π k (x), 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 L h neurons in hidden layer h and multinomial logistic weights. The model is are the parameters to be estimated and ψ is the activation function. Building on the universal approximation theorem for FFNNs (Hornik et al., 1989), Xu and Reich (2021) argue that (9) with H = 1 and large K and L 1 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 W. Unlike a typical statistical learning problem, observational data are not used to estimate W. Rather, the weights are learned from training data generated from the process mixture model with parameters θ 2 ∼ p * , and then a realization from the process over sites i and N i from the model conditioned on θ 2 . 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 θ 2 . Given these values, we generate U(s) at s ∈ {s i , s (i) }. The feature set x i for modeling u i at location i thus contains the spatial parameters θ 2 , process values at the neighboring locations U(s (i) ), as well as the spatial configuration of the 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.
For the local approximation, the weights in (8) are estimated separately for each location. That is, each component f i (u i |θ 2 , u (i) ) in (7) is modeled using its own FFNN, for i = 2, . . . , n. For the global approximation, a single FFNN is used to model the weights in (8) for all the locations. Both methods have their advantages. Each location has a unique spatial configuration of its neighbors for ungridded data, but the local approximation is not affected by these differences. Note that each local model is for an individual location, and conditional on that, the spatial configuration of the neighbors is constant; this reduces the number of features in the model and also allows us to use shallow networks. The global approximation requires more features and, therefore, benefits from a deeper network. It can however be computationally more attractive, as a single FFNN will take significantly less computational resources than n − 1 shallower FFNNs for most real life examples. 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.
Algorithms 1 and 2 outline the global and local procedures, respectively.

Algorithm 1 Global SPQR approximation
Require: Locations s 1 , . . . , s n with neighbor locations s (1) , . . . , s (n) Require: Design distribution p * , sample size N k ← 1 while k ≤ N do Draw sample location s l k , where l k ∈ {2, . . . , n} Draw values of θ 2k ∼ p * , using (2) Generate U(s) = G{V(s)} at s ∈ {s l k , s (l k ) }, using (2) Define features Given the approximate model in (6) for f y with an SPQR approximation for f u , a

Algorithm 2 Local SPQR approximation
Require: Locations s 1 , . . . , s n with neighbor locations s (1) , . . . , s (n) Require: Design distribution p * , training sample size Bayesian analysis using Markov Chain Monte Carlo (MCMC) methods is straightforward. We use Metropolis updates for both θ 1 and θ 2 . For a spatially-varying coefficient model with local GEV coefficients for location i, we update parameters {µ(s i ), σ(s i ), ξ(s i )} as a block sequentially by site, and exploit the Vecchia approximation to use only terms in the likelihood corresponding to sites j i ∈ 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 available at https://github.com/reetamm/SPQR-for-spatial-extremes.

Gaussian process model
To illustrate the potential of the proposed algorithm, we first use data generated from a GP. Of course, the conditional distributions of a GP are univariate Gaussian and so the SPQR approximation is unnecessary. Starting with this simple case, however, will The data in each case are generated from a GP with mean 0, variance 1, and exponential correlation Cor{Y t (s i ), Y t (s j )} = (1 − r)I(i = j) + r exp(−||s i − s j ||/ρ) for r ∈ (0, 1) and ρ > 0. The mean and variance, θ 1 = (µ, σ 2 ), are handled using the Jacobian terms in (6) and not the estimated f u . The proportion of the variance explained by the spatial error, r, and the spatial range ρ make up θ 2 . After drawing θ 2 ∼ p * from the distributions given in Table 1, the t th realization Y t1 , ..., Y t100 is generated from a multivariate normal distribution with correlation defined by θ 2t = (r t , ρ t ). For each Y tj , t > 1, its m neighbors Epochs 20 20 are identified; u t1 = Φ(Y t1 ) is selected as the response, and θ t2 and the m remaining Y tj constitute X t . In practice, using Y tj as features instead of u tj = Φ(Y tj ) provided a better model fit in our studies. For the local SPQR approximation, the feature set consists of Y tj along with log ρ t and log(r t /(1 − r t )). For the global SPQR model, the feature set consists of Y tj , , and log r t 1−r t , where s j = (s j1 , s j2 ) for j = 2, . . . , 11. Hyperparameter tuning for the global and local SPQR approximations was carried out by comparing fitted models on a validation data set of 10 6 observations. Models were compared on the basis of the log-score and the Kullback-Leibler divergence between the estimated and true densities. The lower section of Table 1 lists the hyperparameter configuration chosen for the global and local SPQR approximations. The global model has more layers and higher complexity than the local models, since it contains more features and more variability in the data. Increasing the complexity further in either model results in diminishing improvements in log-scores and KL divergences, and our hyperparameter choices reflect a balance of computational cost and goodness of fit.
We compute the maximum likelihood estimate (MLE) of the model weights W us-  ing the global and local SPQR approximations described in Algorithms 1 and 2, both of which employ stochastic gradient descent with ADAM (Kingma and Ba, 2014). To improve the stability of the global SPQR, we split the training data into ten training sets each of size 10 7 and obtain R = 10 estimates of W, denotedŴ 1 , ...,Ŵ R . The averaged probabilities π k (X) = ∑ R r=1 π r k (X)/R, where π r k (X) is evaluated usingŴ r , are used to evaluate the approximate densities. Each global SPQR model takes approximately 232 minutes to fit, while each local SPQR with all neighbors takes around 21 minutes of computation time. We were able to reduce total computation time using the doParallel package in R to parallelize SPQR model fits.  Figure 4 plots the probability integral transformation (PIT) scores for the two approximations. The PIT for the true model is F(Y tj ) where F is the true Gaussian conditional distribution of Y tj given its neighbors, and this is plotted against the same measure for the fitted CDFF(Y tj ) obtained from the local (left) and global (right) SPQR fits. The PIT statistics falling on the Y = X line shows that the models fit well.
Finally, Figure 5 plots the variable importances (VI) of the five most important variables across all quantiles of the local SPQR model; see Xu and Reich (2021) for more details. These include the spatial parameters ρ and r, as well as the three nearest neighbors based on the Vecchia approximation. These plots show which features are most important for explaining different aspects of the conditional distribution. For both locations 11 (left) and 45 (right), we note that the nearest neighbors have lower variable importance for the extreme quantiles at either end and higher variable importance in the middle. While they are ordered by their importance (the nearest neighbor has highest importance and so on), their positions relative to each other are different for the two locations in Figure 5 and likely depends on the spatial configuration of the neighbors.
The opposite behavior is seen for the spatial parameters r and ρ, which have highest variable importance for the extreme quantiles. The importance of the spatial parameters suggests that they can capture the spatial behavior in the data.
We conduct a simulation study to explore how density-estimation errors propagate to parameter-estimation errors. Each simulated dataset consists of 5 independent realizations of a GP at the 100 locations shown in Figure 2. The GP has mean µ, variance σ 2 and exponential correlation function with range ρ and variance parameter r.
Therefore, there are four parameters to be estimated: θ 1 = (µ, σ) and θ 2 = (ρ, r). For the priors, we assume that µ, log(σ) iid ∼ Normal(0, 10 2 ), log(ρ) ∼ Normal(−2, 1), and log{r/(1 − r)} ∼ Normal(0, 1.5 2 ). Three scenarios are considered with the true values of θ ≡ (µ, σ, ρ, r) set to: For each scenario, we generate 200 datasets and use 10,000 MCMC samples after a burnin of 1,000 iterations. For local SPQR models, runtimes were approximately 13 seconds  local SPQR also has less sampling variance than the global approximation, and its parameter estimates are comparable to those generated using the exact conditionals. This is especially noticeable in estimates for the range ρ, where the global SPQR has sig-  nificantly higher bias and variablility than the other methods. Since nearest-neighbor configurations vary for each location, local SPQRs where each model is already conditioned on a specific spatial configuration seem to be better at estimating the spatial range. Based on these results and taking into consideration the local SPQR's amenability to parallelization despite its higher overall computational cost, we focus on the local SPQR approximation for density estimation in the process mixture model.

Process mixture model
To evaluate the performance of the method on extreme value data, we simulated data from the process mixture model defined in Section 2. Each dataset consists of n = 50 spatial locations distributed on the unit square and 50 independent (time) replications.
To put the MSPs and GPs on the same scale, we assume common smoothness parameter α R = α W = α = 1 and parameterize ρ W and ρ R to give the same effective range, i.e., the distance at which the correlation of the GP reaches 0.05 (h = ρ W log(20)) and the χ-coefficient of the MSP reaches 0.05 (h = ρ R 4Φ −1 (1 − 0.05/2) 2 where Φ is the standard normal distribution function). This results in ρ = ρ W and ρ R = 0.19ρ. For all simulations, the GEV location and scale are µ = 2 and σ = 1 and the spatial dependence parameters are α = 1 and ρ = 0.15.
The simulations vary based on the GEV shape ξ ∈ {−0.1, 0.1} and asymptotic dependence parameter δ ∈ {0.2, 0.8}. We also add a fifth scenario with observations set to be missing (completely at random over space and time) with probability π M = 0.05 and censored below the threshold T, set to the sample medianq 0.5 (over space and time). The scenarios are: For each of the five scenarios defined by ξ and δ we simulate 100 datasets. The locations are ordered by their distance from their origin, and the neighbor sets are comprised of the m = 15 nearest neighbors. Based on the Gaussian simulation in the previous section, we use the local SPQR model with two hidden layers with 30 and 15 neurons, 15 output nodes, a learning rate 0.001, batch size 100 and 50 epochs. We assume that α = 1 and train the SPQR model with design distribution p * , generating samples uniformly on ρ ∈ (0.0, 0.5) and δ ∈ (0, 1). Figure 7a   For the priors, we select µ, log(σ) ∼ Normal(0, 10 2 ), ξ ∼ Normal(0, 0.25 2 ), δ ∼ Uniform(0, 1) and ρ ∼ Uniform(0.0, 0.5). The posterior distribution is approximated using MCMC with 11,000 (21,000 for scenario 5) iterations and Metropolis candidate distributions tuned to have acceptance probability near 0.4. After discarding the first 1,000 iterations as burn-in, the remaining samples are used to compute the posterior mean and 95% interval for each parameter. For the first four scenarios, runtimes were approximately 1 minute per 1,000 MCMC iterations for the first 4 scenarios and about 6 minutes for scenario 5. posterior mean estimator for the GEV parameters generally has low bias and coverage near the nominal level. While the sampling variance of the posterior mean estimator of δ is high, the method is clearly able to distinguish between the two asymptotic regimes with expected value near 0.30 for the asymptotic independence Scenarios 1 and 3 compared to roughly 0.78 for the asymptotic dependence Scenarios 2 and 4. As expected, the sampling variance increases in Scenario 5 with missing data and censoring, but the method is still able to reliably estimate the model parameters.

Data description and exploratory analysis
We apply the proposed methods to model extreme streamflow from 1972-2021 at 489 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 hu-(a) Conditional exceedance χ u (h) for log annual maximum streamflow computed for different distances.
(b) Sample variogram for log annual maximum streamflow, averaged over 50 years of data. Figure 10: Spatial behaviour of log annual maximum streamflow in terms of the conditional exceedance and the variogram. man activities. 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 cubic feet/second, and for each of the T = 50 years and n = 489 stations, we take as the response Y t (s) 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 s. Figure 9 plots the sample 0.9 quantile of the T observations at each station, which show considerable spatial variation.
In order to study the dependence structure of the process, especially at its extremes, we consider the conditional exceedance probability χ u (h) of maximum streamflow at pairs of stations separated by a distance h = ||s j − s k || in kilometers (km) as defined in (4). Figure 10a plots χ 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 10b 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.
For the marginals at each location, we assume GEV distributions with spatio-temporally varying parameters, where X t = (year t − 1996.5)/10 for 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 µ 1 (s) would suggest an increase in the magnitude of the annual extremal streamflow. Figure 11 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. Figure 12 depicts the goodness-offit for the marginal GEV models fitted using MLE. The gray lines indicate the probability integral transforms (PITs) at each location, whereas the black line represents the mean PIT for 489 locations. The mean falls along the Y = X line (red), which establishes the suitability of the log transform for extremal streamflow and reinforces the suitability of the STVC specification for the marginal GEV parameters.
The slope µ 1 (s), the log scale σ(s), and the shape ξ(s) are modeled similarly using GPs, with the spatial range ρ * common to all four processes. For the residual model, we use the process mixture model in Section 2 for spatial dependence and assume independence across years. Since the mean variogram in Figure 10b suggests the presence of a nugget in the residual model, we also add a nugget r to the process mixture. The GP W(s) and MSP R(s) that form our process mixture are thus constructed as follows: Cor{W(s 1 , s 2 )} = r · exp{−(h/ρ)}, where R 1 (s) is an MSP and R 2 (s) iid ∼ GEV(1, 1, 1) distributed independently of R 1 (s).

Results
The local SPQR model for the log of annual streamflow maxima is fitted using settings similar to the simulation study presented in Section 4.2. The local SPQR approximation is thus based around θ 2 = (δ, ρ, r). Single local SPQR models are fitted at each of the 488 locations with neighbors using 2 × 10 6 synthetic observations. Once the local SPQR models have been fitted, we run two MCMC chains for 20,000 iterations each, with two different starting values of δ. The first 5000 iterations from each chain are discarded as burn-in; see Appendix for further computational details.
The posterior means (standard deviations) of the spatial parameters areδ = 0.47 (0.02), ρ = 1004 (80),ρ * = 17907 (1806) andr = 0.56 (0.07). The SVC parameters are therefore much smoother over space than the year-to-year variation captured by ρ. The posterior of δ puts the process in the asymptotic independence regime with high probability, but    In Figure 14a, we see that the log scale σ(s) are highest around south and central US, and in Figure 14b, we note that the posterior means of the shape parameters ξ(s) are negative at all 489 stations, in some cases very negative, on the log transformed data.
We also note that σ(s) are often low at stations where ξ(s) is high, and vice versa. Our primary interest, however, are in estimates of the location parameters across the USA.   Mexico are also of interest since these have relatively low 0.9 quantile values in Figure 9, suggesting that extreme streamflow is starting to have large impacts in these areas.
To study the behavior of extreme streamflow jointly for multiple locations, we con- The joint exceedance for 2021 has a mean of 0.04 with an SD of 0.012. Note that if we assume independence across locations, the joint exceedance probability for the 7 locations would be approximately 10 −7 . Finally, the probability that the joint exceedance in 2021 is higher than in 1972 is 0.975, providing strong evidence in favor of increased extremal streamflow in the area, possibly due to changing climate.

Discussion
In this paper, we proposed a process mixture model 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 process mixture model 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 use 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 he process mixture model was used to analyze changes in annual maximum streamflow within the US over the past 50 years. For this study, we used the annual Future work will focus on the theoretical properties of this model. While it is straightforward to derive an analytical expression for χ u (h) for the trivial case of a shared MSP (see Appendix), i.e., for R(s) = R, obtaining an analytical expression for the general case is more challenging. This would also enable us to investigate the properties of χ u (h) for the process mixture model for different values of δ and as h → ∞. 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. For example, while 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 process mixture model. 3δ−1 , for δ > 0.5.

Appendix A.2: Computational Details
The local SPQR models for extreme streamflow data was fitted using the SPQR package on R 4.x.x. Each model had 2 hidden layers with 30 and 20 neurons respectively, 15 output knots, a batch size of 1000, learning rate of 0.01, and up to 15 neighbors. Of the 200,000 synthetic observations used to fit the models, 80% was used for training and 20% for validation. The model was run for 50 epochs, and the local SPQR takes approximately 14 minutes for locations with a full feature set (i.e., all 15 neighbors).
For the MCMC, we ran the 2 chains in parallel. Each iteration of the MCMC takes around 6 seconds. Figure 16 overlays the trace plots from the two chains, and we see that they are well mixed. The first 5000 iterates from each chain are discarded as burn-in, and the remaining samples from both chains are used to derive our posterior estimates.