inst/doc/POT.R

In martindurocher/floodStat: FloodNet Regional Frequency Analysis

## ------------------------------------------------------------------------
library(floodnetRfa)
data(flowStJohn)

## Extract the year. 
y <- as.numeric(format(flowStJohn$date,'%Y'))

## Verify which year are complete. Print the five years with the less observations.
sort(sapply(split(y,y), length))[1:5]

## Keep complet year
xd <- cbind(flowStJohn[y>=1927,], year = y[y>=1927])

## ------------------------------------------------------------------------
## Fitting POT using MLE
fit <- FitPot(rgpa(2000, 1, 0))
print(fit)

## Show confidence interval for estimated parameters
coef(fit, ci = TRUE)

## Fitting POT using L-moments
coef(FitPot(rgpa(2000, 1, 0), method = 'lmom', varcov = FALSE))

## ----fig.height = 4,fig.width = 6----------------------------------------

## Define a threshold and the minimum separating time
thresh <- 500
r0 <- round(4 + log(14700)) # 14 days

## Extract peaks based on run declustering
peaks1 <- which.clusters(flow~date, xd, u = thresh, r = r0)

## Extract peaks based on WRC recommendations
peaks2 <- which.floodPeaks(flow~date, xd, u = thresh, r = r0, rlow = 0.75)

## Plot the peaks extracted for the year 1927
plot(flow~date, xd[xd$year == 1927,],type = 'l', col = 'grey')
abline(h = thresh, col = 2 ,lwd = 2)
points(flow~date, xd[peaks1,], pch = 16, col = 'blue', cex=2)
points(flow~date, xd[peaks2,], pch = 17, col = 'red')
legend('topleft', col = c('blue','red'), pch = 16:17, 
       legend = c('Cluster','WRC'))


## ------------------------------------------------------------------------
## Fit a POT model after declustering
fit <- FitPot(flow~date, xd, u = 1000, declust = 'wrc', r = r0)
print(fit)

## ---- fig.height=6, fig.width=9------------------------------------------
## Predict flood quantile of return period 10 and 100 years
predict(fit, q = c(.9, .99), se = TRUE, ci = 'profile')

## ---- fig.height=4, fig.width=6------------------------------------------
## List of candidate thresholds
ulst <- seq(500,2000, 10)

## Mean residual life plot
PlotMrl(flow~date, xd, u = ulst, declust = 'wrc', r = r0)

## ---- fig.height=4, fig.width=6------------------------------------------
## Function that extract the p-value of a GOF for a specific threshold.
Fun <- function(u){
  f <- FitPot(flow~date, xd, u = u, declust = 'wrc', r = r0)
  c(u,GofTest(f)$pvalue)
}

## Compute p-values for all candidate thresholds.
candidates <- t(sapply(ulst, Fun))

## Plot the p-values with respect to the candidate thresholds.
plot(candidates[,c(1,2)], type = 'l',
     xlab = 'Threshold', ylab = 'AD test (p-value)')
abline(h = c(.05,.25), lty = 3)

## Print best candidates
u0 <-which(candidates[,2]>.25)[1]
candidates[u0,1]