Using gumboot

Theory

The sampling uncertainty in the GOF estimates is quantified using a mixture of Jackknife and Bootstrap methods. First, we use the Jackknife and Bootstrap methods to compute the standard error in the GOF estimates. These methods resample from the original data sample using the Non-overlapping Block Bootstrap (NBB) strategy using data blocks of length one year. The use of data blocks of length one year reduces the issues with substantial seasonal nonstationarity in shorter data blocks, while preserving the within-year autocorrelation and seasonal periodicity of streamflow series. Bootstrapping methods are only effective if the blocks used are approximately independent. Second, we use the Bootstrap methods to compute tolerance intervals for the GOF estimates, where the 90% tolerance intervals are defined as the difference between the 95th and 5th percentile of the empirical probability distribution of the GOF estimates. Tolerance intervals differ from confidence intervals, because tolerance intervals are intervals corresponding to a random variable, rather than random confidence intervals around some true value. These bootstrap tolerance intervals are computed using 1000 bootstrap samples. Finally, we use the Jackknife-After-Bootstrap method (Efron, 1992) to estimate the standard error in the Bootstrap tolerance intervals, which enables us to evaluate how sensitive the resulting uncertainty intervals are to individual years (blocks).

References: Efron, B. and Gong, G., 1983. A leisurely look at the bootstrap, the jackknife, and cross-validation. The American Statistician, 37(1), pp.36-48.

Efron, B., 1992. Jackknife‐after‐bootstrap standard errors and influence functions. Journal of the Royal Statistical Society: Series B (Methodological), 54(1), pp.83-111.

How to use gumboot

The most important function is bootjack() which computes the bootstap and jackknife statistics for a single data set, which is a data frame containing dates, observed and simulated flows.

A data set (flows_1030500) is provided.

flows_1030500 <- flows_1030500 
head(flows_1030500)
#>         date      obs      sim
#> 1 1989-10-01 0.753092 1.895452
#> 2 1989-10-02 0.679782 1.738436
#> 3 1989-10-03 0.623800 1.556237
#> 4 1989-10-04 0.582480 1.367636
#> 5 1989-10-05 0.545159 1.185386
#> 6 1989-10-06 0.511836 1.017480

Plotting the values shows fairly good agreement between the simulated and observed values

library(ggplot2)
library(reshape2)
melted <- melt(flows_1030500, id.vars = "date")
ggplot(melted, aes(date, value, colour = variable)) +
  geom_line() +
  xlab("") +
  labs(y = bquote('Daily streamflow'~(m^3/s)), x = "")

To perform the bootstrap and jackknife analyses, the values are passed to bootjack. There are many options. The default is to calculate statistics for both NSE (Nash-Sutcliffe efficiency), and KGE (Kling-Gupta efficiency) values. In this example, we will compute the statistics of the NSE.

Note that the name of the observed variable must be obs and the name of the simulated variable must be sim.

NSE_values <- bootjack(flows_1030500, GOF_stat = "NSE")
NSE_values
#>   GOF_stat     seJack     seBoot       p05       p50       p95     score
#> 1      NSE 0.04850763 0.04719664 0.4746045 0.5547303 0.6290767 0.5541248
#>      biasJack      biasBoot     seJab
#> 1 -0.00196514 -0.0007371399 0.0445264

In this example, the standard error of the NSE statistic, as calculated by jackknifing is 0.0485076; as calculated by bootstrapping it is 0.0471966, and as calculated by jackknifing after the bootstrapping (JAB) is 0.0445264, showing the uncertainty in the statistic.

It is important to note that the bootstrap and jackknife-after-bootstrap standard errors are dependant on random samples of the years, so the values will change with each execution of the function. If you want to get the same values each time, for example to compare your results with another person’s analyses, you have two options.

If you set the option seed to have a value, the R random number generator will always return the same sequence of values, as shown below:

bootjack(flows_1030500, GOF_stat = "NSE", seed = 1)
#>   GOF_stat     seJack     seBoot       p05       p50       p95     score
#> 1      NSE 0.04850763 0.04716093 0.4741728 0.5564785 0.6237025 0.5541248
#>      biasJack     biasBoot      seJab
#> 1 -0.00196514 -0.001495579 0.03160307
bootjack(flows_1030500, GOF_stat = "NSE", seed = 1)
#>   GOF_stat     seJack     seBoot       p05       p50       p95     score
#> 1      NSE 0.04850763 0.04716093 0.4741728 0.5564785 0.6237025 0.5541248
#>      biasJack     biasBoot      seJab
#> 1 -0.00196514 -0.001495579 0.03160307

Note that the value of seJack above is identical to the previously determined value, as the jackknifing always uses the same set of years of data.

If you want to compare the results with code not written in R you can save the randomly selected years to a file. If the specified file does not exist, it will be written to. If the file does exist, then the years will be read from it.

NSE_samples <- bootjack(flows_1030500, GOF_stat = "NSE", returnSamples = TRUE)
names(NSE_samples)
#> [1] "statsBoot" "statsJack"

ggplot(NSE_samples$statsBoot, aes(NSE)) + 
  geom_histogram() +
  ggtitle("Bootstrap samples")

Multiple locations

The function CAMELS_bootjack() applies bootjack() to model runs over the “CAMELS” catchments across the contiguous US (CONUS). The model runs are not supplied, and need to be stored in a NetCDF file.

Newman et al. (2015) and Addor et al. (2017) provide details on the hydrometeorological and physiographical characteristics of the CAMELS catchments. The CAMELS catchments are those with minimal human disturbance (i.e., minimal land use changes or disturbances, minimal water withdrawals), and are hence almost exclusively smaller, headwater-type catchments (median basin size of 336 km²).

References: Addor, N., Newman, A.J., Mizukami, N. and Clark, M.P., 2017. The CAMELS data set: catchment attributes and meteorology for large-sample studies. Hydrology and Earth System Sciences, 21(10), pp.5293-5313.

Newman, A.J., Clark, M.P., Sampson, K., Wood, A., Hay, L.E., Bock, A., Viger, R.J., Blodgett, D., Brekke, L., Arnold, J.R. and Hopson, T., 2015. Development of a large-sample watershed-scale hydrometeorological data set for the contiguous USA: data set characteristics and assessment of regional variability in hydrologic model performance. Hydrology and Earth System Sciences, 19(1), p.209.

In addition to the NetCDF file containing the model runs CAMELS_bootjack() a data frame containing the site numbers and their latitudes and longitudes. This is contained in the supplied data frame CAMELS_sites. Note that you can subset the data frame if you only wish to test some of the runs.

CAMELS_bootjack() has most of the same options as bootjack() and returns the same values. Unless you tell it not to (setting quiet = TRUE) a progress bar will be displayed.

This example does the analyses for all of the CAMELS data, returning statistics for runs which were optimised using NSE and KGE run targets, using both NSE and KGE goodness of fit statistics.

CAMELS_sites <- hcdn_conus_sites
nc_file <- "/home/kevin/data/projects/bootstrappR_test/hess2019/results_hcdn_flow.nc"
CAMELS_stats <- CAMELS_bootjack(CAMELS_sites, nc_file)

Having computed the statistics for all CAMELS basins, the uncertainties in the NSE and KGE values can be plotted using ggplot_estimate_uncertainties(), which returns a ggplot2 plot object.

Because some of the CAMELS basins have data sets which do not meet the default criteria for bootjack(), they will return NA_real_ values for their statistics. It is a good idea to first remove these stations from the data set before calling ggplot_estimate_uncertainties(), by using na.omit().

CAMELS_stats_cleaned <- na.omit(CAMELS_stats)
ggplot_estimate_uncertainties(CAMELS_stats_cleaned, fill_colour = "orange")