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?

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) return(rlpoints) }

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) return(ci_df) }

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.

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.

LikeLike

Hi, thanks for the code! I have troubles reproducing getrlpoints part. What does npy stand for?

LikeLike

Hi there,

npy (I assume stands for number per year). It represents the rate of exceeding the threshold. In the fevd function it is set by time.units =.

LikeLike

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.

LikeLike