A Guide to the SpatialGEV Package

Introduction to the GEV-GP Model

The generalized extreme value (GEV) distribution is often used to analyze sequences of maxima within non-overlapping time periods. An example of this type of data is the monthly maximum rainfall levels recorded over years at a weather station. Since there are typically a large number of weather stations within a country or a state, it is more ideal to have a model that can borrow information from nearby weather stations to increase inference and prediction accuracy. Such spatial information is often pooled using the Gaussian process.

The GEV-GP model is a hierarchical model with a data layer and a spatial random effects layer. Let \({\boldsymbol{x}}_1, \ldots, {\boldsymbol{x}}_n \in \mathbb{R}^2\) denote the geographical coordinates of \(n\) locations, and let \(y_{ik}\) denote the extreme value measurement \(k\) at location \(i\), for \(k = 1, \ldots, n_i\). The data layer specifies that each observation \(y_{ik}\) has a generalized extreme value distribution, denoted by \(y \sim \mathrm{GEV}(a, b, s)\), whose CDF is given by \[\begin{equation} F(y\mid a, b, s) = \begin{cases} \exp\left\{-\left(1+s\frac{y-a}{b}\right)^{-\frac{1}{s}}\right\} \ \ &s\neq 0,\\ \exp\left\{-\exp\left(-\frac{y-a}{b}\right)\right\} \ \ &s=0, \end{cases} \label{eqn:gev-distn} \end{equation}\] where \(a\in\mathbb{R}\), \(b>0\), and \(s\in\mathbb{R}\) are location, scale, and shape parameters, respectively. The support of the GEV distribution depends on the parameter values: \(y\) is bounded below by \(a-b/s\) when \(s>0\), bounded above by \(a-b/s\) when \(s<0\), and unbounded when \(s=0\). To capture the spatial dependence in the data, we assume some or all of the GEV parameters in the data layer are spatially varying. Thus they are introduced as random effects in the model.

A Gaussian process \(z({\boldsymbol{x}})\sim \mathcal{GP}(\mu({\boldsymbol{x}}_{cov}), k({\boldsymbol{x}}, {\boldsymbol{x}}'))\) is fully characterized by its mean \(\mu({\boldsymbol{x}}_{cov})\) and kernel function \(k({\boldsymbol{x}}, {\boldsymbol{x}}') = {\mathrm{Cov}}( z({\boldsymbol{x}}), z({\boldsymbol{x}}') )\), which captures the strength of the spatial correlation between locations. The mean is a function of parameters \(\boldsymbol{\beta}\) and covariates \({\boldsymbol{x}}_{cov}=(x_1,\ldots,x_p)'\). We assume that given the locations, the data follow independent GEV distributions each with their own parameters. The complete GEV-GP hierarchical model then becomes \[\begin{equation} \begin{aligned} y_{ik} \mid a({\boldsymbol{x}}_i), b({\boldsymbol{x}}_i), s & {\overset{ind}{\sim}}\mathrm{GEV}\big( a({\boldsymbol{x}}_i), \exp( b({\boldsymbol{x}}_i) ), \exp(s)\big)\\ a({\boldsymbol{x}}) \mid \boldsymbol{\beta}_a, {\boldsymbol{\theta}}_a &\sim \mathcal{GP}\big( {\boldsymbol{X}}_a\boldsymbol{\beta}_a, k({\boldsymbol{x}}, {\boldsymbol{x}}' \mid {\boldsymbol{\theta}}_a) \big)\\ log(b)({\boldsymbol{x}}) \mid \boldsymbol{\beta}_b, {\boldsymbol{\theta}}_b &\sim \mathcal{GP}\big( {\boldsymbol{X}}_b\boldsymbol{\beta}_b, k({\boldsymbol{x}}, {\boldsymbol{x}}' \mid {\boldsymbol{\theta}}_b) \big)\\ s({\boldsymbol{x}}) \mid \boldsymbol{\beta}_s, {\boldsymbol{\theta}}_s &\sim \mathcal{GP}\big( {\boldsymbol{X}}_s\boldsymbol{\beta}_s, k({\boldsymbol{x}}, {\boldsymbol{x}}' \mid {\boldsymbol{\theta}}_a) \big). \end{aligned} \end{equation}\] In this package, a uniform prior \(\pi({\boldsymbol{\theta}}) \propto 1\) is specified on the fixed effect and hyperparameters \[{\boldsymbol{\theta}}=(s, \boldsymbol{\beta}_a, \boldsymbol{\beta}_b, \boldsymbol{\beta}_s, {\boldsymbol{\theta}}_a, {\boldsymbol{\theta}}_b, {\boldsymbol{\theta}}_s).\]

What Does SpatialGEV Do?

The package provides an interface to estimate the approximate joint posterior distribution of the spatial random effects \(a\), \(b\) and \(s\) in the GEV-GP model. The main functionalities of the package are:

Details about the approximate posterior inference can be found in Chen, Ramezan, and Lysy (2022).

Installation

SpatialGEV depends on the package TMB to perform the Laplace approximation. Make sure you have TMB installed following their instruction before installing SpatialGEV. Moreover, SpatialGEV uses several functions from the INLA package for approximating the Matérn covariance with the SPDE representation and for creating meshes on the spatial domain. If the user would like to use the SPDE approximation, please first install INLA. Since INLA is not on CRAN, it needs to be downloaded following their instruction here.

To install the latest version of SpatialGEV, run the following:

devtools::install_github("meixichen/SpatialGEV")

Using the SpatialGEV Package

Exploratory analysis

We now demonstrate how to use this package through a simulation study. The simulated data used in this example comes with the package as a list variable simulatedData2, which contains the following:

a, logb and logs are simulated from Gaussian random fields using the R package SpatialExtremes. Using the simulated GEV parameters, we generate 50 to 70 observed data \({\boldsymbol{y}}_i\) at each location \(i, \ i=1,\ldots,400\).

Spatial variation of \(a\), \(\log(b)\) and \(\log(s)\) can be viewed by plotting them on regular lattices:

plot of chunk sim1-plot-abs

Number of observations at each location is shown in the figure below.

plot of chunk sim1-plot-num-obs-per-loc

Below are histograms of observations at \(8\) randomly sampled locations.

plot of chunk sim-1plot-y

Model fitting

The function spatialGEV_fit() in the SpatialGEV package is used to find the mode and quadrature of the Laplace approximation to the marginal hyperparameter posterior, which is the basis of a Normal approximation to the posterior of both hyperparameters and random effects referred to as Laplace-MQ in our paper. The function spatialGEV_sample() is then used to sample from this distribution via a sparse representation of its precision matrix.

To fit the GEV-GP model to this simulated dataset, the first step is calling the spatialGEV_fit() function, for which several arguments must be provided:

There are other arguments which user might want to specify to override the defaults:

The code below fits a GEV-GP model with Matérn SPDE kernel to the simulated data. The shape parameter is constrained to be positive. No covariates are included in this model, so by default an intercept parameter \(\beta_0\) is estimated for the GP of each spatial random effect. We also specify the return levels of interest at 0.5 and 0.9.

The point estimates and the associated standard errors (calculated via the Delta method) for the first 5 locations can be obtained as follows.

Posterior sampling

To obtain posterior samples of \(a\), \(b\), and \(s\), we pass the fitted model object fit to spatialGEV_sample(), which takes in three arguments:

The following line of code draws 2000 samples from the posterior distribution and the posterior predictive distribution:

Then use summary() to view the summary statistics of the posterior samples.

pos_summary <- summary(sam)
pos_summary$param_summary[c(1:5, 
                            (n_loc+1):(n_loc+5), 
                            (2*n_loc+1):(2*n_loc+5),
                            (3*n_loc):(nrow(pos_summary$param_summary))),]
#>                  2.5%        25%        50%        75%       97.5%       mean
#> a1          59.412539 60.5606787 61.1760198 61.7450187 62.81015397 61.1426639
#> a2          59.539727 60.5454594 61.0884820 61.6491748 62.55328984 61.0895984
#> a3          59.581149 60.5069739 61.0161594 61.4839396 62.34364563 60.9965115
#> a4          59.513730 60.4250367 60.8875294 61.3477945 62.33539179 60.8892916
#> a5          59.296019 60.2661151 60.7286919 61.2282584 62.15330084 60.7418139
#> log_b1       2.877005  2.9294962  2.9567235  2.9825525  3.03232503  2.9561419
#> log_b2       2.893010  2.9397056  2.9636650  2.9864951  3.03059358  2.9627034
#> log_b3       2.915094  2.9525660  2.9739638  2.9940322  3.03259781  2.9734103
#> log_b4       2.924364  2.9631214  2.9819243  3.0019072  3.03939871  2.9822945
#> log_b5       2.939242  2.9758818  2.9943037  3.0133760  3.04701306  2.9940867
#> s1          -3.839856 -3.0841189 -2.7307542 -2.3706166 -1.65738654 -2.7388903
#> s2          -3.573531 -2.9649495 -2.6447043 -2.3158603 -1.72693538 -2.6423994
#> s3          -3.389748 -2.8394872 -2.5343934 -2.2401899 -1.68024250 -2.5359593
#> s4          -3.360504 -2.8003850 -2.5143843 -2.2420109 -1.63099629 -2.5115080
#> s5          -3.344200 -2.8611489 -2.6160538 -2.3815445 -1.84034036 -2.6154145
#> s400        -4.414861 -3.6845410 -3.3234348 -2.9691244 -2.23024853 -3.3197572
#> beta_a      59.254986 59.8170439 60.1140164 60.3904082 60.97734166 60.1139866
#> beta_b       2.930002  2.9900897  3.0248382  3.0565364  3.11982666  3.0233803
#> beta_s      -2.815150 -2.5094988 -2.3405593 -2.1767840 -1.86375833 -2.3421628
#> log_sigma_a -2.168320 -1.3661709 -0.9866477 -0.6044773  0.05501866 -0.9967687
#> log_kappa_a -1.448783 -0.8093071 -0.4702390 -0.1377747  0.55418411 -0.4640518
#> log_sigma_b -4.298673 -3.5310981 -3.0931269 -2.6657391 -1.87501163 -3.0952154
#> log_kappa_b -2.002020 -1.4089644 -1.0952879 -0.7842804 -0.19887461 -1.0943473
#> log_sigma_s -2.305678 -1.8257401 -1.5657913 -1.3120916 -0.85406613 -1.5674691
#> log_kappa_s -0.830086 -0.5131705 -0.3312056 -0.1599720  0.17159183 -0.3330161
pos_summary$y_summary[1:5,]
#>        2.5%      25%      50%      75%    97.5%     mean
#> y1 37.16310 55.05763 68.63354 85.51421 143.4623 73.74704
#> y2 35.63351 54.13532 67.33381 85.98369 137.5352 72.47793
#> y3 37.35344 54.85376 68.40772 85.42123 142.6234 73.70693
#> y4 36.04354 54.88409 67.70335 85.37151 142.4695 73.51183
#> y5 36.61537 54.21058 67.26347 87.65584 147.4782 74.02922

Model checking

Since we know the true values of \(a\), \(b\), and \(s\) in this simulation study, we are able to compare the posterior mean with the true values. The posterior means of \(a\), \(b\) and \(s\) at different locations are plotted against the true values below.

plot of chunk sim1-pos-plots

Posterior prediction

Finally, we demonstrate how to make predictions at new locations. This is done using the spatialGEV_predict() function, which takes the following arguments:

In this dataset, the GEV parameters \(a\) and \(\log(b)\) are generated from the surfaces below, and the shape parameter is a constant \(\exp(-2)\) across space.

plot of chunk sim2-plot-ab

We randomly sample 100 locations from the simulated dataset simulatedData as test locations which are left out. Data from the rest 300 training locations are used for model fitting. We fit a GEV-GP model with random \(a\) and \(b\) and log-transformed \(s\) to the training dataset. The kernel used here is the SPDE approximation to the Matérn covariance function (Lindgren, Rue, and Lindström (2011)).

The fitted model object is passed to spatialGEV_predict() for 500 samples from the posterior predictive distributions. Note that this might take some time.

Then we call summary() on the pred object to obtain summary statistics of the posterior predictive samples at the test locations.

Since we have the true observations at the test locations, we can compare summary statistics of the true observations to those of the posterior predictive distributions. In the figures below, each circle represents a test location.

plot of chunk sim2-pred-plot

Case study: Yearly maximum snowfall data in Ontario, Canada

In this section, we show how to use the SpatialGEV package to analyze a real dataset. The data used here are the 1987-2021 monthly total snowfall data (in cm) obtained from Environment and Natural Resources, Government of Canada. The link to download the raw data is https://climate-change.canada.ca/climate-data/#/monthly-climate-summaries. This dataset is automatically loaded with the package and is named ONsnow.

library(dplyr)
#library(maps)
lon_range <- c(-96, -73)
lat_range <-  c(41.5, 55)
summary(ONsnow)
#>     LATITUDE       LONGITUDE      STATION_NAME       CLIMATE_IDENTIFIER   LOCAL_YEAR    LOCAL_MONTH     TOTAL_SNOWFALL  
#>  Min.   :41.75   Min.   :-94.62   Length:63945       Length:63945       Min.   :1987   Min.   : 1.000   Min.   :  0.00  
#>  1st Qu.:43.58   1st Qu.:-81.50   Class :character   Class :character   1st Qu.:1991   1st Qu.: 3.000   1st Qu.:  0.00  
#>  Median :44.23   Median :-79.93   Mode  :character   Mode  :character   Median :1997   Median : 6.000   Median :  1.00  
#>  Mean   :44.96   Mean   :-80.88                                         Mean   :1999   Mean   : 6.469   Mean   : 15.23  
#>  3rd Qu.:45.50   3rd Qu.:-78.97                                         3rd Qu.:2005   3rd Qu.: 9.000   3rd Qu.: 23.60  
#>  Max.   :54.98   Max.   :-74.47                                         Max.   :2021   Max.   :12.000   Max.   :326.00
maps::map(xlim = lon_range, ylim = lat_range)
points(ONsnow$LONGITUDE, ONsnow$LATITUDE)

plot of chunk snow-display-raw-data

Data preprocessing

We first grid the data using cells of length \(0.5^{\circ}\). By doing this, weather stations that are apart by less than \(0.5^{\circ}\) in longitude/latitude are grouped together in the same grid cell. From now on, we refer to each grid cell as a location.

For each location, we find the maximum snowfall amount each year and only keep locations where there are at least two years of records.

plot of chunk snow-display-grid-data

Now we fit the GEV-GP model to the data using the SPDE kernel. Both \(a\) and \(b\) are treated as spatial random effects. \(s\) is constrained to be a positive constant. Note that here we have specified a \(\mathcal{N}(-5,5)\) prior on the log-transformed shape parameter. This is because we found that the shape parameter is estimated close to 0 and such a prior ensures model fitting procedure is numerically stable.

Next, 5000 samples are drawn from the joint posterior distribution of all parameters.

Instead of using the return_levels argument in spatialGEV_fit(), we can also use the posterior samples to calculate the 5-year, 10-year, and 25-year return level posterior samples, which are the upper \(1/5\)%, \(1/10\)%, and \(1/25\)% quantiles of the extreme value distributions at these locations.

Plotted below are the return levels.

plot of chunk snow-return-plots

plot of chunk snow-return-plots

References

Chen, Meixi, Reza Ramezan, and Martin Lysy. 2022. “Fast Approximate Inference for Spatial Extreme Value Models.” arXiv. https://arxiv.org/abs/2110.07051.

Hrafnkelsson, Birgir, Stefan Siegert, Raphael Huser, Haakon Bakka, and Árni V. Jóhannesson. 2021. “Max-and-Smooth: A Two-Step Approach for Approximate Bayesian Inference in Latent Gaussian Models.” Bayesian Analysis 16 (2): 611–38.

Lindgren, F. K., H. Rue, and J. Lindström. 2011. “An Explicit Link Between Gaussian Fields and Gaussian Markov Random Fields: The Stochastic Partial Differential Equation Approach.” Journal of the Royal Statistical Society, Series B 73: 423–98.

V. Jóhannesson, Stefan Siegert, Raphael Huser, Haakon Bakka, and Birgir Hrafnkelsson. 2022. “Approximate Bayesian Inference for Analysis of Spatiotemporal Flood Frequency Data.” The Annals of Applied Statistics 16 (2): 905–35.