1

I'm using the POT package to carry out certain calculations in R. The output of the analysis is stored in a object of class uvtop. Now, I'd like to export the result of the analysis, rather than just plot it within an R window.

Here it goes an example, using sample data from this package.

data(ardieres)
events1 <- clust(ardieres, u = 6, tim.cond = 8/365, clust.max = TRUE)
npy <- length(events1[,"obs"]) / (diff(range(ardieres[,"time"], na.rm
= TRUE)) - diff(ardieres[c(20945,20947),"time"]))
mle <- fitgpd(events1[,"obs"], thresh = 6, est = "mle")
par(mfrow=c(2,2))
plot(mle, npy = npy)

With this, I get the image below: test output

OK, but what I want is to export the necessary data to reproduce the Return Level Plot (bottom right panel) somewhere else, i.e. the values represented by circles, the solid and both dashed lines.

Rui Barradas
  • 70,273
  • 8
  • 34
  • 66
Pythonist
  • 1,937
  • 1
  • 14
  • 25
  • Please post a reproducible example. Read, for instance, [How to make a great R reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). It is almost surely possible to do what you want but we need much more information than just a package name and the object class. (Produced by function `fitgpd`?) – Rui Barradas Nov 02 '18 at 11:57
  • Thanks, I have updated the question accordingly. – Pythonist Nov 02 '18 at 12:15
  • 1
    The third instruction `npy <- length(etc` is incomplete. – Rui Barradas Nov 02 '18 at 12:26
  • What do you mean by *"data to reproduce the Return Level Plot (bottom right panel) somewhere else"*? Does somewhere else mean: in R on someone else's computer, excel, python, SAS, etc.? – emilliman5 Nov 02 '18 at 12:31
  • I have fixed it now. Thanks. – Pythonist Nov 02 '18 at 12:31
  • @emilliman5, what I mean is that I want an ASCII file with the coordinates of the circles, other with the solid line, and other with the dashed lines. This way, I can reproduce the figure with whatever other plotting software I'm more familiar with. – Pythonist Nov 02 '18 at 12:33
  • 1
    What you are asking requires some detective work. First of all run `str(mle)` to know what is in the object to be plotted. Then, run `getAnywhere(plot.uvpot)`. This will tell you that the function that plots the 4th graph is `getAnywhere(retlev.uvpot)`. And if you now have the source code you can know what are the data structures you need. – Rui Barradas Nov 02 '18 at 12:53

2 Answers2

1

To get the data that is plotted for Return Level, we have to dig into the retlev function. Basically, I did my best to strip all of the plotting out and construct a data.frame of the required points.

getRetLevData <- function (fitted, npy) {
  data  <- fitted$exceed
  loc   <- fitted$threshold[1]
  scale <- fitted$param["scale"]
  shape <- fitted$param["shape"]
  n     <- fitted$nat

  pot.fun <- function(T) {
    p <- rp2prob(T, npy)[, "prob"]
    return(qgpd(p, loc, scale, shape))
  }

  eps <- 10^(-3)

  if (!is.null(fitted$noy)){ 
    npy <- n/fitted$noy
    } else if (missing(npy)) {
    warning("Argument ``npy'' is missing. Setting it to 1.")
    npy <- 1
  }

  xlimsup    <- prob2rp((n - 0.35)/n, npy)[, "retper"]
  fittedLine <- pot.fun(seq(1/npy + eps, xlimsup, length.out = n))
  p_emp      <- (1:n - 0.35)/n
  xPoints    <- 1/(npy * (1 - p_emp))
  empPoints  <- sort(data)
  samp       <- rgpd(1000 * n, loc, scale, shape)
  samp       <- matrix(samp, n, 1000)
  samp       <- apply(samp, 2, sort)
  samp       <- apply(samp, 1, sort)
  ci_inf     <- samp[25, ]
  ci_sup     <- samp[975, ]

  rst <- data.frame(xPoints, fittedLine, empPoints, 
                    ci_inf, ci_sup)
}

x <- getRetLevData(mle, npy)
head(x)
#    fittedX fittedLine  xPoints empPoints   ci_inf   ci_sup
#1  1.001000   6.003716 1.011535      6.09 6.001557 6.239971
#2  3.891288  11.678503 1.029810      6.09 6.014536 6.363070
#3  6.781577  14.402517 1.048758      6.09 6.042065 6.470195
#4  9.671865  16.282306 1.068416      6.19 6.074348 6.583290
#5 12.562153  17.740710 1.088825      6.44 6.114193 6.684118
#6 15.452441  18.942354 1.110029      6.45 6.146098 6.812058
write.csv(x, "my_pot_results.csv")

An overlay of the extracted data vs retlev plot. The CI's are a bit different because of sampling.

enter image description here

emilliman5
  • 5,816
  • 3
  • 27
  • 37
  • Thanks, it looks like black magic to be, but it works. What puzzles me most is that '0.35'. Do you mind commenting on this part? Anyway, I'll try now to apply to my particular problem, and then I'll mark the question as solved. – Pythonist Nov 02 '18 at 14:24
  • I have no idea why `0.35` is used in the math or what the `pot.fun` does except call `rp2prob`. – emilliman5 Nov 02 '18 at 14:27
  • Thanks, I it finally works, but it requires a final tweak (that took me several hours of desperation and reverse engineering to figure out). The x-coordinates for the fitted curve are different as for the dots and the dashed lines. I had to edit the function to output the two x-vectors, but it works nicely now. Thanks! – Pythonist Nov 04 '18 at 16:05
  • I noticed that. That's why my function returns `fittedx` and `xpoints`, I should have made that more clear – emilliman5 Nov 04 '18 at 16:14
0

If you don't want to read the file with other applications than R, just save it with:

save(mle, file="myfilename.Rdata")

or

saveRDS(mle, file="myfilename.Rds")

To read it back in, load the library that generated the data and then use 

load("myfilename.Rdata")

to load the object back to your workspace or 

mle <- readRDS("myfilename.Rds")

save saves more environments together with the object than saveRDS depending on how the library is implemented this might make a difference. I'd suggest to use save unless the datasets are getting too large.

snaut
  • 2,261
  • 18
  • 37