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


Something a little different to my usual posts and targeted to people interested in Extreme Value Analysis (EVA).

I have been working on an analysis of extreme sea levels as part of my PhD on flooding risk. If you are familiar with EVA in R you will have likely come across the extRemes package, written by Eric Gilleland. While I am not (yet!) an expert in EVA, I have found this package to be quite accessible.

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?

Example of default plots from the extRemes package (from Gilleland and Waltz, 2016)

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!

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 (here called pot.mle) you will get the information you need to draw a Return Level plot with the following two functions (other types of plots will use a different section from plot.fevd.mle):

The first returns a dataframe of return level point coordinates:

getrlpoints <- function(pot.mle){

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


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

getcidf <- function(pot.mle){

rperiods = c(2, 5, 10, 20, 50, 80, 100, 120, 200, 250, 300, 500, 800)
bds <- ci(pot.mle, 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) 


For example:

rlpoints <- getrlpoints(pot.mle)
ci_df <- getcidf(pot.mle)

And that’s all you need to produce a good looking return level plot!

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 = 2) +
 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), limits = c(1,1000))+
 theme_bw() +
 theme(axis.text = element_text(size=12), axis.title = element_text(size = 14, face = "bold"))



Gilleland, E. and Katz, R. W. (2016) extRemes 2.0: An Extreme Value Analysis Package in R. Journal of Statistical Software, 72:8.


Author: Ulysse Pasquier

Ulysse Pasquier is a writer and PhD Candidate at the University of East Anglia researching how to adapt to climate change and manage future flood risk

5 thoughts on “[R] Plotting results from extRemes package – Extreme Value Analysis”

  1. What’s up,I log on to your blog named “[R] Plotting results from extRemes package – Extreme Value Analysis – Ulysse Pasquier” regularly.Your humoristic style is witty, keep it up! And you can look our website about love spell.


  2. Hi, thanks for the useful codes. However, as the “getrlpoints” function uses “npy” which is only available in the POT model, it will not work for the GEV model. Could you please revise the function? Thank you.


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s