---
title: "FEH statistical analysis"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{FEH statistical analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```


## Introduction
This quick guide provides an overview of the core FEH functions that can be used to generate extreme flow estimates from single-site or pooling (gauged & ungauged) approaches. It is assumed that the reader has knowledge of the theory of these approaches; this guide is intended to demonstrate how to use the UKFE package rather than to explain the FEH statistical method.

The following table provides the key functions for implementing the FEH methods.

```{r echo = FALSE, results = 'asis'}
library(knitr)
kable(data.frame(Function = c("GetCDs/CDsXML", "GetAM/AMImport", "QuickResults", "QMED", "DonAdj", "Pool", "H2", "DiagPlots", "Zdists", "PoolEst", "Uncertainty", "GenLog/GEV/GenPareto/Gumbel/Kappa3", "POTextract or POTt", "AnnualStat"), Purpose = c("Get catchment descriptors using the NRFA gauge reference number or import them from XML files downloaded from the FEH Web Service", "Get annual maximum data from sites suitable for pooling or from AM files", "Provides default pooling results (gauged, ungauged & fake ungauged), directly from catchment descriptors", "Estimates QMED from catchment descriptors with donor options (either donor IDs or a number of donors - which are then automatically chosen based on distance)", "Provides a list of donor candidates based on proximity to the study catchment", "Creates a pooling group with the option of excluding one or more sites", "The heterogeneity measure for pooling groups", "Provides a number of diagnostic plots to compare the sites in a pooling group", "The Zdist statistic to determine the distribution with the best fit to the pooling group", "Flood estimates from the pooling group", "Quantifies uncertainty for gauged (enhanced single-site) and ungauged pooling groups", "Four functions for each distribution to fit the distribution and obtain flood frequency estimates and growth factors", "Functions for extracting peaks over threshold data", "Extract annual statistics. The default is maximum, but any statistic can be used such as sum or mean.")))
```

We need to load the package before we can use any of the functions:
```{r setup}
library(UKFE)
```

## Overview and assumptions of the FEH statistical method

This section provides a brief overview of the FEH statistical method, and its core assumptions. Examples of how to apply the method can be found later.

The FEH statistical method comes under the broader term of 'regional frequency analysis' and employs the "index flood" procedure for estimating peak flows at a single location on a river line. The procedure has two main steps, which are multiplied together for a final result:

* Estimate the index flood (in the FEH method, the index flood is QMED - the median annual maximum flow).

* Calculate the growth curve: This provides growth factors as a function of return period and is a standardised model of the AMAX distribution. In this case a growth factor is the extreme flow of interestd (such as the 100-year flow) divided by QMED.

Multiplying the index flood by the growth curve scales the standardised growth factor to produce the flow estimate for the catchment of interest.
When the catchment is gauged, the index flood is estimated as the median of the observed annual maximum sample (assuming suitable data quality). When the catchment is ungauged, the index flood is estimated using a multiple regression equation, using catchment descriptors as predictors.

The growth curve is estimated by pooling gauged data from catchments which are similar to the catchment of interest. A statistical distribution is assumed. Then, weighted average parameters of the fitted distribution (L-moment ratios) from the AMAX samples in the pooling group are used to estimate the growth curve.

There are numerous choices and possible adjustments which can be made for each step. The UKFE package, in general, provides default settings with no adjustments; however, these can be changed by the user.

There are several key underlying assumptions in FEH statistical analysis. As stated in Hosking and Wallis (1997), the index-flood procedure makes the following assumptions:

* Observations at any given site are identically distributed.

* Observations at any given site are serially independent.

* Observations at different sites are independent.

* Frequency distributions at different sites are identical apart from a scale factor.

* The mathematical form of the regional growth curve is correctly specified.

The first of these leads to the assumption of stationarity (i.e. observations of past flood events are assumed to represent the behaviour of future events). Where this is not the case, non-stationary methods can be applied. For example, non-stationary distributions can be fitted using other software such as the `nonstat` R package (Hecht and Zitzmann, 2025).

The QMED regression equation was calibrated using an approach similar to the generalised least squares method. This allows for covariance of the model residuals (errors).
Firstly, the observed QMEDs (on which it was calibrated) were logarithmically transformed. The catchment descriptors used were also transformed. With this in mind, the assumptions of the QMED equation are as follows:

* The relationship between the transformed independent variables (predictors) and the logarithmically transformed dependent variable (QMED) is linear.

* Observations (observed QMEDs used in calibration) are independent of each other.

* Model and sampling errors are independent of each other.

* Model and sampling errors are normally distributed and have a mean of zero.

* Predictor variables are independent.

* The cross-correlation of model errors can be described by the distance between catchment centroids, and the form of the associated correlation matrix is known for the calibration process. 

## Getting catchment descriptors

One of the first steps in any FEH analysis is to obtain catchment descriptors (CDs). See the 'Data input' article for information on functions for importing CDs.

## Quick results

Once you have some catchment descriptors, you can start to use many of the FEH based functions. The full process, covered shortly, is encompassed in the `QuickResults()` function. This provides results from a default pooling group, automatically accounts for urbanisation, and provides a donor transfer from 8 sites by default (this number can be adjusted). Note that if the site is gauged and in the PeakFlowData dataframe (and the descriptors are the same) the QMED will default to the observed QMED. The default distribution is the generalised logistic distribution. Further details and assumptions are provided in the function's help file. 

First, we'll obtain some CDs for site 55002 and then apply the QuickResults function.
The results from the default analysis at gauge 55002 are printed to the console. This is a list with three elements. The first is a data frame with the following columns: return period (RP), peak flow estimates (Q) and growth factor estimates (GF). The second element is the estimated L-CV and L-skew (linear coefficients of variation and skewness). The third element provides distribution parameters (location, scale, and shape) for the chosen distribution (generalised logistic by ddault).  

```{r}
# Get catchment descriptors for NRFA site 55002 and save to an object called CDs.55002
CDs.55002 <- GetCDs(55002)

# Obtain estimates and plots for the gauged scenario using default settings
QuickResults(CDs.55002)
```


## Estimating QMED

There is a `QMED()` function which provides the ungauged QMED (median annual maximum flow) estimate from catchment descriptors. There is an option to adjust this using donor sites. The simplest case is to provide catchment descriptors with no donor transfer:

```{r}
# Estimate QMED from the catchment descriptors for site 55002
QMED(CDs.55002)
```

By default, this gives urban adjusted QMED results. 

To undertake a donor transfer (with the aim of improving the estimate), you can state how many donors you wish to apply with the no.Donors argument (this applied the closest - by catchment centroid - n donors). Or you can specify a vector of IDs. In which case a list of donor site candidates is necessary. For this, there is the `DonAdj()` function, which provides the closest *N* gauges (default of 10) to the site (geographically by catchment centroid) that are suitable for QMED, along with catchment descriptors.  

```{r}
# Identify donor catchments based on the catchment descriptors for site 55002
DonAdj(CDs.55002)
```

You can apply one or multiple donors using the DonorIDs argument within the `QMED()` function. Or, as noted just specify a number of donors. The following, shows both of these. Note that the one specifying no.Donors = 3, will default to the observed QMED because site 55002 is both the catchment of interest and the closest donor and will have all the weighting in the adjustment. Lastly you can also output all the detail of the donors by using the 

```{r}
# Estimate QMED from the catchment descriptors for site 55002 using gauges with NRFA
# IDs 55007 and 55016 as donor sites
QMED(CDs.55002, DonorIDs = c(55007, 55016, 55013))

# Use the closest three donor sites
QMED(CDs.55002, no.Donors = 3)

#Use the specified three donors and return the associated details of the adjustment
QMED(CDs.55002, DonorIDs = c(55007, 55016, 55013), ReturnDetails = TRUE)

```




The user should be aware of the following specifics about the UKFE implementation of QMED estimation. However, note that there is an inherent flexibility in 'R', so any code can be edited.

* Note that the `QMEDDonEq()` function can also be used, and a number of donors could be used in turn - the user would need to choose the donor weights. 

* The `DonAdj()` function only provides candidate donor sites from the NRFA "suitable for QMED and pooling" dataset. However, the `QMEDDonEq()` can be applied, which is entirely flexible regarding inputs.

* By default, the `QMED()` function uses the URBEXT2015 catchment descriptor with an urban expansion (set to current year). There is an argument within the function called `UrbanExpansion`, which is set to `TRUE`. This can be set to `FALSE` if required. There is also the `UEF()` function, which can be applied independently, and for any year. 


## Creating a pooling group

To create a pooling group, the `Pool()` function can be used. If it is for a gauged site, the site will be automatically included in the group, unless it is has an URBEXt2015 above UrbMax (default 0.03). If you wish to include the urban site, the `include` argument can be set be used (just put in the gauge ID). Alternatively UrbMax can be increased (note that this will allow other ubran sites into the group). There is a `DeUrb` argument which is `TRUE` by default. It de-urbanises the LCVs of site in the pooling group. To create and view the pooling group:

```{r}
# Create a pooling group based on the catchment descriptors for site 55002 and 
# save to an object called Pool.55002
Pool.55002 <- Pool(CDs = CDs.55002)

# View the pooling group
Pool.55002
```

You may now wish to assess the group for heterogeneity and see which distribution fits best. The `H2()` function can be used to check for heterogeneity.

```{r, echo = FALSE}
set.seed(5)
```

```{r}
# Check for heterogeneity
H2(Pool.55002)
```

This group appears to be homogeneous, but if a group is heterogeneous, you may wish to use the `DiagPlots()` function to check for any outliers etc. You can also look at the discordancy measures, which can be found in the last two columns of the pooling group created.

```{r fig.alt="Diagnostic plot consisting of 10 panels used to aid review of catchment descriptors. The first seven panels show histograms of all catchments, with each of the ten sites marked by small crosses, for the following variables: Area, SAAR, PROPWET, FARL, FPEXT, BFIHOST19, and URBEXT2000. Two scatter plots display L-skew versus L-CV and L-skew versus L-kurtosis. The next panel shows the catchment locations as points on a UK outline, plotted by Eastings and Northings. The final panel shows the average seasonality of peaks for all the sites in the pooling group. Together, the plots allow comparison of these sites against the distribution of all catchments. The ten sites appear generally central and unremarkable across descriptor space, indicating they are fairly typical relative to the wider catchment set."}
# Create diagnostic plots to aid review of the pooling group
DiagPlots(Pool.55002)
```

Let's pretend that site 8010 and site 76017 are to be removed (because they have the highest discordancy - in reality, a more thorough analysis would be undertaken). To remove the sites, we recreate the pooling group and use the `exclude` option. The function automatically adds the next site(s) with the lowest similarity distance measure if the total number of years of AMAX in the group has dropped below 500 (or whatever we set as the `N` argument in the `Pool()` function, which has a default of 800).

```{r}
# Recreate the pooling group, excluding NRFA sites 8010 and 76017
Pool.55002 <- Pool(CDs = CDs.55002, exclude = c(8010, 76017))

# View the new pooling group
Pool.55002

# Check for heterogeneity
H2(Pool.55002)
```

### Choosing a distribution. 

Once you are happy with the pooling group, you can check for the best distribution. The `Zdists()` function can be used to do this. This compares the Generalised Extreme Value (GEV), Generalised logistic (GenLog), Gumbel, and Kappa3 distributions. The method behind it is outlined in Hosking and Wallis (1997). Note that it is not the same as the method detailed in Science Report SC050050 (Kjeldsen *et al*., 2008). The difference is described in the function's help file and is necessary to encompass the two-parameter Gumbel distribution.

```{r}
# Calculate goodness-of-fit measure
Zdists(Pool.55002)
```

The output is a list, where the first element contains the Z-scores for each distribution, and the second element states which distribution has the best fit.  

The L-CV and L-skew for each site are provided in the pooling group and are derived from the AMAX data. If the user has further information, for example, an extra year of data to provide a new annual maximum, the L-CV and L-skew can be calculated using the `Lcv()` & `LSkew()` functions. The user's data can be pasted in or read into the 'R' environment. Data can be read in using a function such as `read.csv()`, or a column of data can be pasted in from Excel by running the following code and then pressing "Ctrl+V" in the console (assuming you have already copied the required data):

```{r, eval = FALSE}
MyData <- scan()
```

This creates a numeric vector called `MyData` that contains your flow data. The `Lcv()` & `LSkew()` functions take a numeric vector of flows as their only input, so you can then use `Lcv(MyData)` and `LSkew(MyData)` to derive the L-moment ratios.

The relevant results can be updated in the pooling group by applying the `LRatioChange()` function. For example, if you wanted a new pooling group with updated L-moment ratios for site 55002, the following code could be used:

```{r}
# Update the pooling group to use user-supplied L-CV and L-skew values for site 55002
UpdatedPool <- LRatioChange(Pool.55002, SiteID = 55002, lcv = 0.187, lskew = 0.164)

# View the updated pooling group
UpdatedPool
```

## Final estimates (the growth curve and results)

To derive the final estimates, the `PoolEst()` function can be applied to the pooling group. It is necessary to provide an input for the `QMED` argument. Therefore, we will first obtain a QMED estimate (in this case, we will use the data available from the NRFA). To estimate QMED, we can either use the `GetAM()` function and assess the data manually, or we can get it more directly from the `PeakFlowData` data frame embedded within the package (these are calculated as the median of the AMAX data (rejected years excluded) from the NRFA Peak Flow Dataset). The following code provides the latter.

```{r}
# Extract QMED from the QMEDdata data frame within the UKFE package
GetQMED(55002)
```

We can use functions within functions so we can include the above in `PoolEst()` (two things to note; we need to include the CDs and the Gauged = TRUE argument only impacts the uncertainty calculation):

```{r}
# Derive the growth curve and final flow estimates based on the pooling group and QMED
PoolEst(Pool.55002, CDs = CDs.55002, QMED = GetQMED(55002), Gauged = TRUE)
```

The results are printed to the console. This is a list with three elements. The first is a data frame with the following columns: return period (RP), peak flow estimates (Q), growth factor estimates (GF), and factorial standard error. The latter is calculating by bootstrapping the pooling group and sampling from the QMED uncertainty. The second element is the estimated pooled L-CV and L-skew (linear coefficients of variation and skewness). The third element provide parameters (location, scale, and, shape) for the chosen distribution. 

The user should be aware of the following specifics about how UKFE employs the pooled method:

* The default QMED FSE is set to 1.5. This value represents the factorial standard error associated with the combined model and sampling uncertainty of the QMED model. You can update QMED FSE with the fseQMED argument. Generally it should be lower if donor catchments are applied or if the QMED is estimated from a gauged sample. 

* The default distribution applied is the generalised logistic distribution.

* The final FSE comprises the FSE associated with QMED estimation and growth curve estimation. The FSE associated with the growth curve quantifies sampling error (aleatoric uncertainty) only. The FSE associated with QMED includes model uncertainty. In general, the uncertainty associated with the modelling process or the hydrometric data collection is not included in the FSE. Many choices are made during the process, and deriving a range of results from a range of justifiable choices is a way to quantify some of the epistemic uncertainty.

* By default, urban adjustment is applied. The user can change this with the `UrbAdj` argument, by setting it to `FALSE`. 

The `EVPool()` function can be used to plot a growth curve.

```{r fig.alt="Plot showing growth curves from FEH QuickResults calculations. Pooled estimates are shown as a solid line, single-site estimates as a dashed green line, and observed flood data as blue circles. The x-axis shows the logistic reduced variate and the y-axis shows Q/QMED, that is, flow divided by the median annual flood. A scale for the return period is displayed at the bottom right. The plot illustrates how well the pooled and single-site growth curves map onto the observed flood data. Neither curve matches all observations exactly, but both provide a reasonable fit overall. Both curves follow the expected growth-curve shape, with increasing discharge at higher return periods."}
# Plot the growth curve
EVPool(Pool.55002, gauged = TRUE)
```

For the **ungauged catchment**, the same process can be followed. First, we'll exclude the site from our pooling group. We will assume we went through a rigorous pooling group analysis and concluded all was well:

```{r fig.alt="Plot showing how multiple single-site growth curves are pooled to form the default ‘fake ungauged’ growth curve. The x-axis shows the logistic reduced variate and the y-axis shows Q/QMED, that is, flow divided by the median annual flood. Several solid black lines represent growth curves from individual sites, and a red dashed line represents the pooled growth curve with the target site excluded. All curves follow the expected growth-curve shape, with increasing discharge at higher return periods. The pooled curve lies within the spread of the individual site curves and, as expected, closely resembles their average, illustrating how pooling smooths variability to produce a representative ungauged growth curve."}
# Recreate the pooling group, excluding NRFA gauge 55002 to make it an ungauged analysis
PoolUG.55002 <- Pool(CDs.55002, exclude = 55002)

# View the new pooling group
PoolUG.55002

# Estimate QMED from the catchment descriptors with two donors (note that the QMED 
# value needs to be extracted using `$QMEDs.adj` as there are three elements to the 
# output when two donors are used and the PoolEst function expects just the QMED value)
CDsQmed.55002 <- QMED(CDs.55002, DonorIDs = c(55007, 55016))

# Use the pooling group and the ungauged QMED estimate to provide the pooled estimates
Results55002 <- PoolEst(PoolUG.55002, QMED = CDsQmed.55002, CDs = CDs.55002)

# Print the results
Results55002

# Plot a growth curve for the pooling group
EVPool(PoolUG.55002)
```

The same results are provided as for the gauged case.

## Single site analysis

The quickest way to get results for a single sample are the `GEV()`, `GenLog()`, `Kappa3()`, `Gumbel()` and `GenPareto()` functions. For example, the following code can be used to obtain a 75-year return period estimate for an AMAX sample called `MyAMAX` using the generalised logistic distribution:

```{r, eval = FALSE}
GenLogAM(MyAMAX, RP = 75)
```

An annual maximum series can be obtained for sites suitable for pooling and QMED using the `GetAM()`. The `AMImport()` or GetDataNRFA functions can also be used (see the 'Data input' article). We will retrieve and plot the AMAX for our site 55002:

```{r fig.alt="Bar chart of annual maximum river flow. The x-axis shows years, and the y-axis shows peak flow in cubic meters per second. Each bar represents the highest flow in that year. The flows vary from year to year, with several notably high peaks in recent years."}
# Extract the AMAX data for NRFA site 55002
AM.55002 <- GetAM(55002)

# View the head of the AMAX series
head(AM.55002)

# Plot the AMAX data
AMplot(AM.55002)
```

The `AMplot()` function returns a bar plot of the AMAX series.

There are four functions for each of the GEV, GenLog, Gumbel, Kappa3, and GenPareto distributions.
In the function names below, `XXX` should be replaced by the distribution prefix (`GEV`, `GenLog`, `Gumbel`, `Kapp3`, or `GenPareto`).
For the Generalised Pareto distribution, the first function is named `GenParetoPOT()` instead of `GenParetoAM()`.

* `XXXAM()` - estimates return levels or return periods directly from an AMAX (or POT for GenPareto) series.

* `XXXEst()` - estimates return levels or return periods from user-supplied parameters.

* `XXXGF()` - provides the growth factor as a function of L-CV, L-skew, and RP (the Gumbel version only requires the L-CV and RP).

* `XXXPars` - estimates distribution parameters from an AMAX or POT series (using L-moments, or MLE in some cases).

We can look at the fit of the distributions using the `ERPlot()` or `EVPlot()` functions. By default, `ERPlot()` compares the percentage difference between simulated results with the observed for each rank of the data. Another option (see the `ERType` argument) compares the simulated flows for each rank of the sample with the observed flows of the same rank. For both plots, 500 simulated samples are used. The function's help file provides more information, including how to compare the distribution estimated from pooling with the observed data. This is an updated version of that described in Hammond (2019).

The `EVPlot()` function plots the extreme value frequency curve or growth curve with observed sample points. The uncertainty is quantified by parametric bootstrapping.

The following code generates an extreme rank plot and an extreme value plot using the GEV distribution:

```{r fig.alt="Extreme rank frequency curve illustrating uncertainty around the estimated growth curve. The x-axis shows the rank and the y-axis shows the percent difference. Observed data are shown as blue circles, plotting along the 0% line. The central modelled curve, shown as a solid line, tracks closely to the observations, while dashed lines represent the 90% bootstrapped confidence interval. The uncertainty bounds remain narrow through the central portion of the curve and expand slightly toward the lower and upper ranks."}
# Generate an extreme rank plot for the AMAX data for NRFA site 55002 and the GEV 
# distribution and set a user-defined title (using the 'main' argument)
ERPlot(AM.55002$Flow, dist = "GEV", main = "Extreme rank plot for NRFA site 55002: GEV")
```

```{r fig.alt="Extreme rank frequency curve illustrating uncertainty around the estimated growth curve. The x-axis shows the logistic reduced variate and the y-axis shows Q/QMED, that is, flow divided by the median annual flood. Observed data are shown as blue circles, forming a concave-up increasing curve. The central modelled estimate, shown as a solid line, tracks these observations closely. Dashed lines represent the 90% bootstrapped confidence interval, which remains tight across the centre of the distribution and widens toward the extremes."}
# Generate an extreme value plot for the AMAX data for NRFA site 55002 and the GEV
# distribution and set a user-defined title (using the 'Title' argument)
EVPlot(AM.55002$Flow, dist = "GEV", Title = "Extreme value plot for NRFA site 55002: GEV")
```

As an example, we can estimate the 100-year flow. Then, estimate the return period of the maximum observed flow. *Make sure to select the "Flow" column*.

```{r}
# Estimate the 100-year flow directly from the AMAX data for NRFA site 55002 using
# the GEV distribution
GEVAM(AM.55002$Flow, RP = 100)

# Estimate the return period of the maximum observed flow directly from the AMAX 
# data using the GEV distribution
GEVAM(AM.55002$Flow, q = max(AM.55002$Flow))
```

We can also estimate the parameters of the GEV from the data, using the default L-moments, and estimate quantiles as a function of return period (and vice versa) from user input parameters.

```{r}
# Estimate the parameters of the GEV distribution from the AMAX data for NRFA site 55002
GEVPars55002 <- GEVPars(AM.55002$Flow)

# View the parameters
GEVPars55002

# Estimate the 75-year flow using these GEV parameters 
GEVEst(loc = GEVPars55002$Loc, scale = GEVPars55002$Scale, shape = GEVPars55002$Shape, RP = 75)
```

If we derive L-moments for the site, we can then estimate growth factors. The `LMoments()` function does this. We can also calculate L-CV and L-skew separately using the `Lcv()` and `LSkew()` functions.

```{r}
# Estimate the 75-year GEV growth factor from the L-moments for the AMAX data for 
# NRFA site 55002
GEVGF(lcv = Lcv(AM.55002$Flow), lskew= LSkew(AM.55002$Flow), RP = 75)
```

We can estimate the 75-year flow by combining this latter result with QMED (the median of the AMAX record).

```{r}
# Estimate the 75-year flow by multiplying the growth factor by the median AMAX flow
GEVGF(lcv = Lcv(AM.55002$Flow), lskew= LSkew(AM.55002$Flow), RP = 75) * median(AM.55002$Flow)
```

This gives a different 75-year flow estimate from the more direct estimate using the GEV parameters. This is because the former calculates the T-year peak directly, whereas the latter assumes a distribution scaled by the index flood.

## Peaks over threshold

A peaks over threshold (POT) series can also be obtained from a time series using the `POTextract()` function. A daily sample of the Thames at Kingston flow is included in the package for example purposes (the data frame `ThamesPQ`). The `POTextract()` function can be used on a single vector or a data frame, with dates in the first column and the time series of interest in the second. In the `ThamesPQ` data frame, the Thames flow is in the third column and the date in the first. Therefore, the POT series can be extracted as follows:

```{r fig.alt="Time series plot of flow showing extraction of peaks over threshold events. The x-axis shows time and the y-axis shows flow magnitude. Red circles mark exceedances above the blue threshold line. A green line indicates div, which defines independence between peaks; in this case it is set to the mean. Exceedances occur throughout the record, with a lull around 2005, followed by a marked increase in frequency afterwards."}
# Extract a POT series from the Thames at Kingston daily mean flow data, selecting 
# the first and third columns from the ThamesPQ data frame which contain the date 
# and flow data respectively. The threshold is set as the 90th percentile.
POT.Thames <- POTextract(ThamesPQ[, c(1, 3)], thresh = 0.9)

# View the output
POT.Thames
```

This produces a figure showing the daily mean flow data with peaks over the threshold shown as red circles. The blue line is the threshold, and the green line is the mean flow. Independent peaks are chosen above the threshold. If peaks are separated by a drop below the mean flow (green line), they are considered independent. More details can be found in the function's help file (including other options for defining independence and how to rename axis labels and the plot title).

Note that you can also use the `POTt()` function for extracting peaks over threshold. It works significantly faster, but only has the option of separating peaks by the number of timesteps. 

The POT approach is similar to the AMAX approach, except that we need to convert from the POT RP scale to the annual RP scale. To do so, the functions have a peaks per year (ppy) argument. The `POTextract()` function provides the number of peaks per year (in this case, 1.867). To estimate the associated 100-year daily mean flow peak, the GenPareto function can be used as follows:

```{r}
# Estimate the 100-year flow for the Thames at Kingston using the extracted POT series 
GenParetoPOT(POT.Thames$peak, ppy = 1.867, RP = 100)
```

## Uncertainty

Methods for quantifying aleatoric (sampling) uncertainty are provided here for three estimation methods: single-site, gauged pooling, and ungauged pooling. The uncertainty is predicated on forming a sampling distribution for the statistic of interest, and in the three cases (single site, gauged pooling, and ungauged pooling), it is done using bootstrapping. In the pooling approaches (gauged and ungauged), the whole pooling group is bootstrapped, and the weighted estimates provided on every iteration.

### Single-site

Extract an annual maximum sample with the following code.

```{r}
# Extract the AMAX for NRFA site 55002
AM.55002 <- GetAM(55002)
```

The uncertainty for the 100-year flow can be estimated using the `Bootstrap()` function. The 95% confidence intervals are provided by default, but this can be adjusted using the `Conf` argument. The lower limit is provided in the 'Lower' column and the upper limit in the 'Upper' column.

Note that the `Bootstrap()` function can also be applied for any function/statistic, and the sampling distribution developed can be returned. 

```{r}
# Estimate uncertainty for the GEV-estimated 100-year flow
Bootstrap(AM.55002$Flow, Stat = GEVAM, RP = 100)
```

### Ungauged pooling

The `Uncertainty()` function can be used for the pooled estimates to quantify the sampling uncertainty. The estimated QMED and the pooling group are required. Furthermore, the factorial standard error needs to be provided for the ungauged case, using the `fseQMED` argument (the default is 1.5). Note that by default the generalised logistic distribution is selected for the estimates, but this can be changed.

The following example to estimate uncertainty for the ungauged pooled estimates for NRFA site 55002 uses the catchment descriptors object created earlier (`CDs.55002`). This provides a dataframe with return period in the first column and Factorial Standard Error (FSE) in the second. For confidence intervals you can divide and multiply your estimates for each return period by the FSE. For 95% intervals do the same but raise the FSE to the power of 1.96. 


```{r}
# Estimate QMED for site 55002 from catchment descriptors using two donors
QMEDEst <- QMED(CDs.55002, DonorIDs = c(55007, 55016))

# Create an ungauged pooling group for site 55002
PoolUG.55002 <- Pool(CDs.55002, exclude = 55002)

# Estimate the uncertainty of the pooling results 
Uncertainty(PoolUG.55002, QMEDEstimate = QMEDEst)
```

### Gauged pooling

The gauged case is similar to the above, but the `Gauged` argument is set to `TRUE`, and there is no need to set `QMEDEstimate`; it is derived from the observed AMAX if `Gauged = TRUE`. 

It is important to note that when `Gauged` equals `TRUE`, the top site in the pooling group is assumed to be the "gauged" site of interest. This will happen by default because the similarity distance measure should be zero - assuming the that "gauged" site is within the `NRFAData` data frame.

The following example uses the gauged pooling group we previously derived:

```{r}
# Create a gauged pooling group for site 55002
Pool.55002 <- Pool(CDs = CDs.55002, exclude = c(8010, 76017))

# Estimate the uncertainty of the pooling results 
Uncertainty(Pool.55002, Gauged = TRUE)
```

## Exporting results
For use outside of 'R', outputs can be saved as objects and then written to CSV files. For example, to save the peak flow estimates from earlier:

```{r, eval = FALSE}
# Save the peak flow estimates to an object called 'Results.55002'
Results.55002 <- PoolEst(PoolUG.55002, QMED = CDsQmed.55002)[[1]]

# Write to csv
write.csv(Results.55002, "my/file/path/Results55002.csv", row.names = FALSE)
```
