dgevbhm vignette
In this vignette, we present an example workflow to apply our R code for fitting the duration-dependent Generalized Extreme Value distribution (dGEV) in a Bayesian Hierarchical framework (Rischmuller et al. 2026). We use sub-daily precipitation records from the Global Sub-Daily Rainfall dataset (Lewis et al. 2019) over the United Kingdom.
library(dgevbhm)
#> Package dgevbhm loaded successfully.
#> NOTE: Always update using pak::pkg_install('alexleer/dgevbhm') to avoid long C++ compile times.The data_read() function reads in sub-daily
precipitation data directly from a directory dir_path of
text files. missing_data specifies the minimum percentage
of missing data per year; otherwise, the year is disregarded from the
data. durs are the durations of the precipitation maxima
for IDF/dGEV fitting. cores allows the user to specify the
number of CPU cores for parallel processing, significantly speeding up
the data reading process.
To reproduce this workflow on your own machine, you will first need to obtain the raw sub-daily precipitation text files (e.g., from the Global Sub-Daily Rainfall dataset). Please download the data and place the files into a local directory on your computer.
Once your directory is set up, you can read the data into R by
pointing the dir_path argument to your chosen folder:
# Replace "/path/to/directory/" with the actual path where you saved the data
UK_data_full <- data_read(dir_path = "/path/to/directory/",
missing_data = 5,
durs = c(1, 3, 6, 12, 24, 48),
cores = 4)(Note: For this vignette, the processed UK_data_full object is already bundled and loaded with the package.)
Out of 1903 files, 47 stations were selected based on the requirements mentioned above.
Selecting the rain gauges
dgev_bhm() is the central function of our package. It is
able to fit the Bayesian Hierarchical dGEV models for a chosen sample of
rain gauges. There, the user can choose if the shape parameter of the
dGEV is fixed over durations (resulting in the same shape parameter for
all durations, \(\xi_d\)) or fixed over
space (resulting in the same shape parameter for all rain gauges, \(\xi_j\)).
# mean and median hourly precipitation rate (mm)
mean(UK_data_full$amp)
#> [1] 0.1024251
median(UK_data_full$amp)
#> [1] 0.0882582Below we show an example of selecting similar stations using annual
mean precipitation amp:
med_amp <- median(UK_data_full$amp)
# similar stations: selection of 25 stations that are closest to the median precipitation rate
sim_stations <- which(rank(abs(UK_data_full$amp - med_amp)) <= 25)
# plotting a histogram of the 25 precipitation rates
hist(UK_data_full$amp[sim_stations], main = 'histogram', xlab = 'stations')
# clean the data: select the 25 rain gauges and store them in a new object
UK_data <- list()
UK_data$y <- UK_data_full$y[sim_stations]
UK_data$amp <- UK_data_full$amp[sim_stations]
UK_data$latlon <- UK_data_full$latlon[sim_stations, ]
UK_data$stationID <- UK_data_full$stationID[sim_stations]
# don't forget the durations!
UK_data$durs <- c(1, 3, 6, 12, 24, 48)Fit the Bayesian hierarchical models
# Fit the dGEV-BHM with the shape parameter fixed over durations (shp_d = 'd'). We opt for 2000 iterations of the MCMC sampler, where the first half is used as burn-in.
# We assume that 2 cores are available for computation.
# As you can only use one core per chain, we use 2 chains.
# WARNING: May take some time. You can also load the prepared fitted models.
fit_xi_d <- dgev_bhm(data = UK_data, chains = 2, iter = 2000, cores = 2, shp_d = 'd')
#>
#> Computing initial parameters... Complete!
#>
#> Loading pre-compiled Stan model (durational)... Complete!
#>
#> SAMPLING FOR MODEL 'dgev_bhm_xi_d' NOW (CHAIN 1).
#>
#> SAMPLING FOR MODEL 'dgev_bhm_xi_d' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 0.00265 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 26.5 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.002564 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 25.64 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 115.911 seconds (Warm-up)
#> Chain 1: 84.154 seconds (Sampling)
#> Chain 1: 200.065 seconds (Total)
#> Chain 1:
#> Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 133.551 seconds (Warm-up)
#> Chain 2: 86.945 seconds (Sampling)
#> Chain 2: 220.496 seconds (Total)
#> Chain 2:# Fit the dGEV-BHM with the shape parameter fixed over space (shp_d = 'j').
# WARNING: May take some time. You can also load the prepared fitted models.
fit_xi_j <- dgev_bhm(data = UK_data, chains = 2, iter = 2000, cores = 2, shp_d = 'j')
#>
#> Computing initial parameters... Complete!
#>
#> Loading pre-compiled Stan model (spatial)... Complete!
#>
#> SAMPLING FOR MODEL 'dgev_bhm_xi_j' NOW (CHAIN 1).
#>
#> SAMPLING FOR MODEL 'dgev_bhm_xi_j' NOW (CHAIN 2).
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.003085 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 30.85 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 2:
#> Chain 2: Gradient evaluation took 0.00306 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 30.6 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 114.441 seconds (Warm-up)
#> Chain 2: 83.456 seconds (Sampling)
#> Chain 2: 197.897 seconds (Total)
#> Chain 2:
#> Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 116.834 seconds (Warm-up)
#> Chain 1: 84.748 seconds (Sampling)
#> Chain 1: 201.582 seconds (Total)
#> Chain 1:# structure of the fitted models
names(fit_xi_d)
#> [1] "pars" "data" "information_criteria"
#> [4] "mcmc_diagnostics" "durs" "j"
#> [7] "d" "shp_d" "mcsamp"
#> [10] "stationID"
names(fit_xi_j)
#> [1] "pars" "data" "information_criteria"
#> [4] "mcmc_diagnostics" "durs" "j"
#> [7] "d" "shp_d" "mcsamp"
#> [10] "stationID"There are 10 list items in the fitted model:
pars– posterior parameters samplesdata– original input datainformation_criteria- widely applicable information criteria
- leave-one-out information criteria
mcmc_diagnostics- R-hat
- Effective sample size
durs– duration vector of the modelj– number of localitiesd– number of durationsshp_d– hierarchical dimension with respect to space"j"or durations"d"mcsamp– MCMC sample sizestationID– vector of station IDs
Bayesian summary statistics
From the fitted dgev_bhm() object, we can obtain MCMC
convergence diagnostics.
Here for the parameter \(\tilde{\mu}\):
# MCMC diagnostics for mut
fit_xi_d$mcmc_diagnostics$rhat$mut
#> mut[1] mut[2] mut[3] mut[4] mut[5] mut[6] mut[7] mut[8]
#> 0.9990515 0.9995442 0.9997762 1.0006114 0.9996957 1.0011561 0.9992602 0.9991278
#> mut[9] mut[10] mut[11] mut[12] mut[13] mut[14] mut[15] mut[16]
#> 0.9999091 1.0002548 0.9990023 0.9991194 1.0007376 1.0002650 1.0000417 0.9993420
#> mut[17] mut[18] mut[19] mut[20] mut[21] mut[22] mut[23] mut[24]
#> 1.0001756 0.9992059 0.9995288 0.9994842 0.9993190 0.9999393 0.9995993 0.9991572
#> mut[25]
#> 0.9997953
fit_xi_d$mcmc_diagnostics$ess$mut
#> mut[1] mut[2] mut[3] mut[4] mut[5] mut[6] mut[7] mut[8]
#> 2802.562 3064.385 3543.107 2901.723 3629.802 3441.356 3368.871 4296.135
#> mut[9] mut[10] mut[11] mut[12] mut[13] mut[14] mut[15] mut[16]
#> 3332.189 3474.975 3974.260 2762.584 2584.031 3484.062 3536.488 3055.644
#> mut[17] mut[18] mut[19] mut[20] mut[21] mut[22] mut[23] mut[24]
#> 2962.470 3078.171 3292.360 2817.474 3516.275 2638.827 3279.816 2851.406
#> mut[25]
#> 3808.894The \(\hat{R}\) values for \(\tilde{\mu}\) are near 1.000, indicating perfect MCMC convergence and chain mixing. Furthermore, the large Effective Sample Sizes (5,700–8,700) confirm an abundance of independent draws. Together, these diagnostics guarantee highly stable and reliable parameter estimates for the dGEV model.
Information criteria
Examine out-of-sample predictive accuracy with information criteria (IC).
fit_xi_j$information_criteria
#> $waic
#> Estimate SE
#> elpd_waic -5763.9974 82.836328
#> p_waic 110.5185 4.134955
#> waic 11527.9947 165.672657
#>
#> $looic
#> Estimate SE
#> elpd_loo -5765.4518 82.911794
#> p_loo 111.9729 4.346921
#> looic 11530.9036 165.823588
fit_xi_d$information_criteria
#> $waic
#> Estimate SE
#> elpd_waic -5787.58149 83.58895
#> p_waic 97.76875 4.44916
#> waic 11575.16298 167.17791
#>
#> $looic
#> Estimate SE
#> elpd_loo -5788.68305 83.652387
#> p_loo 98.87032 4.718404
#> looic 11577.36611 167.304773The dGEV-BHM-\(\xi_j\) is the better model, because the IC values are lower, indicating superior out-of-sample predictive accuracy.
Checking the model fits
Now we will proceed with goodness-of-fit metrics.
Probability-probability plots:
Ideally, the dots should line up close to the 1:1 diagonal. Furthermore, no single duration should deviate significantly, as this could indicate that the dGEV model is inappropriate for capturing duration dependence. While both models yield visually similar results in the PP-plots, the information criteria previously established that dGEV-BHM-\(\xi_j\) provides a better overall fit for this dataset.
Posterior predictive metrics
The figures below generated by the post_pred_plot
function are posterior predictive histograms. The histogram generates
the theoretical statistic (e.g., quantiles, median, max, min) of the
posterior parameter samples. In our case, the statistics are computed
using the GEV itself. We compare this with the empirical statistic from
the data sample. With this, we get a sense of the goodness of fit as
well as the uncertainty.
par(mfrow = c(2, 1))
post_pred_plot(fit = fit_xi_d, stat = 'median', j = 1, d = 6) # posterior predictive diagnostic plot for rain gauge #1 and duration 48h for the fit_xi_d model
post_pred_plot(fit = fit_xi_j, stat = 'median', j = 1, d = 6) # posterior predictive diagnostic plot for rain gauge #1 and duration 48h for the fit_xi_j model
par(mfrow = c(2, 1))
post_pred_plot(fit = fit_xi_d, stat = 'min', j = 6, d = 1) # posterior predictive diagnostic plot for rain gauge #6 and duration 1h for the fit_xi_d model
post_pred_plot(fit = fit_xi_j, stat = 'min', j = 6, d = 1) # posterior predictive diagnostic plot for rain gauge #6 and duration 1h for the fit_xi_j model
In the posterior predictive histograms, not all stations are going to
fit perfectly. This is especially, and almost exclusively, true for
'min' and 'max' statistics. This is due to the
asymptotic nature of the GEV, which models the limit of block extremes
(the histogram) and therefore often fails to capture a single, rare
absolute record (the red line) in a finite sample.
par(mfrow = c(2, 1))
post_pred_plot(fit = fit_xi_d, stat = 'max', j = 3, d = 1)
post_pred_plot(fit = fit_xi_j, stat = 'max', j = 3, d = 1)
Instead of using a diagram, we can calculate the posterior predictive
p-value using the post_pred_pval() function. The inputs are
identical to the post_pred_plot() function. An ideal fit
for the model would yield a p-value of 0.5.
post_pred_pval(fit = fit_xi_d, stat = 'median', j = 2, d = 6)
#> [1] 0.248
post_pred_pval(fit = fit_xi_j, stat = 'median', j = 2, d = 6)
#> [1] 0.327While p-values can fluctuate across different gauge samples, the dGEV-BHM-\(\xi_j\) model’s median aligns with the quantitative information criteria results, demonstrating a reliable and superior fit for this specific spatial sample.
Return level plots
The inputs are:
j- spatial localityd- durationmax_rp- maximum plotted return periodalpha- significance level for credible intervals (defaultalpha=0.05)
# let's try a significance level of 0.1
return_level_plot(fit_xi_d, j = 5, d = 1, max_rp = 100, alpha = 0.1)
Intensity-duration-frequency curves
The inputs are:
j- spatial localityrp- return periods (i.e. frequency)alpha- significance level for credible intervals (defaultalpha=0.05)
Here we have a single IDF curve on the plot:
idf_plot(fit = fit_xi_d,j = 5, rp = 100, alpha = 0.05)
#> Computing d-dimensional shape IDF curve...
Multiple IDF curves on the plot:
idf_plot(fit = fit_xi_d,j = 5, rp = c(2,10,100), alpha = 0.05)
#> Computing d-dimensional shape IDF curve...
The IDF curves are the final result, and they can be produced for all 25 rain gauges, where the user can select return period(s) and the confidence interval. Rischmuller et al. (2026) have shown the added value of the dGEV Bayesian Hierarchical estimation compared to standard station-wise GEV or dGEV-MLE estimations in terms of higher accuracy and robustness.