[R] Plotting results from extRemes package - Extreme Value Analysis

If you are familiar with Extreme Value Analysis (EVA) in R you will have likely come across the extRemes package, written by Eric Gilleland.

The package’s main function, ‘fevd()’, fits an extreme value distribution to data. Methods for plotting are available for objects of class ‘fevd’ to produce diagnostics and result plots, such as a return level plot.

My issue is the restricted control you seem to have over the plot’s layout and general look. You can choose which plots to show, or to add a custom title, but not much else. What if you want to customise the plot to your liking, plot the results in another software/package or display two EVAs on the same figure as shown above?

I could not find an easy way to extract the information the package uses to generate the plots, nor could I find any mention of this problem online. In the off-chance someone is in the same situation, you’ll find a solution below. It seems a little awkward to me and if there is a better way that I am missing, please do let me know!

Let’s use the data from the “extRemes” package, with the example for the Peaks Over Threshold approach:

library(extRemes)
data("damage", package = "extRemes")
fitD <- fevd(Dam, damage, threshold = 6, type = "GP", time.units = "2.06/year")
plot(fitD, type = "rl")

You have to play with the original code for the package’s plot functions. When using the Maximum Likelihood Estimation (MLE) method, the function to change is ‘plot.fevd.mle’. If you are satisfied with the overall default aspect of the plots and just need to tweak it you can customise the code and reload the function. However, I generally use ggplot2 to produce my figures and like to have a consistent template.

Once you have your ‘fevd’ object you will get the information you need to draw a Return Level plot with the following two functions (other types of plots and distributions will use a different section from plot.fevd.mle).

The first returns a dataframe of return level point coordinates:

getrlpoints <- function(fit){

xp2 <- ppoints(fit$n, a = 0)
ytmp <- datagrabber(fit)
y <- c(ytmp[, 1])
sdat <- sort(y)
npy <- fit$npy
u <- fit$threshold
rlpoints.x <- -1/log(xp2)[sdat > u]/npy
rlpoints.y <- sdat[sdat > u]
rlpoints <- data.frame(rlpoints.x, rlpoints.y)

return(rlpoints)
}

The second returns a dataframe of estimated returned levels with confidence intervals:

getcidf <- function(fit){

rperiods = c(2, 5, 10, 20, 50, 80, 100, 120, 200, 250, 300, 500, 800)
bds <- ci(fit, return.period = rperiods)
c1 <- as.numeric(bds[,1])
c2 <- as.numeric(bds[,2])
c3 <- as.numeric(bds[,3])
ci_df <- data.frame(c1, c2, c3, rperiods) 

return(ci_df)
}

That is all you need to produce the plot in ggplot2:

library(ggplot2)

rlpoints <- getrlpoints(fitD)
ci_df <- getcidf(fitD)

ggplot() +
  geom_line(data = ci_df, aes(x = rperiods, y = c2), color = "red") +
  geom_line(data = ci_df, aes(x = rperiods, y = c1), color = "grey") +
  geom_line(data = ci_df, aes(x = rperiods, y = c3), color = "grey") +
  geom_point(data = rlpoints, aes(x = rlpoints.x, y = rlpoints.y), size = 1) +
  ylab("Return Level (m)") +
  xlab("Return Period (Years)") +
  scale_x_log10(expand = c(0, 0), breaks = c(2,5,10,20,50,100,200,500)) +
  theme_bw() +
  theme(axis.text = element_text(size=12), axis.title = element_text(size = 14, face = "bold"))