[R] Saving fExtremes estimates and k-block return level with confidence intervals.

[Peter Maclean]

Hi:
I am trying to compare the results of "evir" and "fExtreme" packages. I could 
not figure out how to save the "evir" package results. Also, how to pass the 
results to "fExtreme" function "gevrlevelPlot" and "evir" function "rlevel.gev" 
to get the return levels. I just need the values and not graphs. 

 
#Example
library(fExtremes)
library(evir)
require(plyr)
y<- data.frame(rgev(300, xi = 0.1, mu = .5, sigma = 1.6))
colnames(y) <- c("y")
n <- data.frame(n = rep(c(1,2,3), each = 100))
z <- cbind(n,y)
y <- z$y

# Model for grouped data
##fExtremes package
z1 <- split(z,z$n)
rgf <- lapply(z1, function(x){
              m <- as.numeric(x$y)
              gevFit(m, block = 2, type = "mle")
 }) 
#Save results
resf<- ldply(rgf, function(x) x at fit$par.ests)

#Qs: How to transfer rge object to "gevrlevelPlot" function to get the values? 
 
##evir package
rge<- by(z,z[,"n"], function(x) gev(x$y, 2, method = "BFGS", control =list(maxit 
= 10000)))  

 
#Qs:How to save par.ests and par.ses? 
#Qs:How to transfer rge object to "rlevel.gev" function to get the 
values?           

 
#Model for single vector
rlf <- gevFit(y, block = 2, type = "mle") 
rlfp <- gevrlevelPlot(rlf, kBlocks = 2)
rlfp
rle <- gev(y, 2, method = "BFGS", control =list(maxit = 10000))
rlep<- rlevel.gev(rle, k.blocks = 2, add = FALSE)
rlep




----- Original Message ----
From: Dennis Murphy <djmuser at gmail.com>
To: Peter Maclean <pmaclean2011 at yahoo.com>
Sent: Thu, June 30, 2011 3:03:05 PM
Subject: Re: [R] Saving fExtremes estimates and k-block return level with 
confidence intervals.

Hi:

The plyr package can help you out as well. As far as the estimates go,

library(plyr)
> ldply(res2, function(x) x at fit$par.ests)
  .id        xi        mu      beta
1  1 0.1033614 2.5389580 0.9092611
2  2 0.3401922 0.5192882 1.5290615
3  3 0.5130798 0.5668308 1.2105666

You could also look into the l_ply() function, which takes a list as
input and outputs nothing; however, it is usually called if one want
to make a series of similar plots. You need to write a function that
takes a generic list component and outputs a plot. Here's a fairly
trivial example:

testdf <- data.frame(gp = rep(LETTERS[1:3], each = 10), x = 1:10,
                      y = rnorm(30))
l_ply(split(testdf, testdf$gp), function(df) plot(y ~ x, data = df))

Alternatively,

# Observe that a data frame is used as input
pfun <- function(df) plot(y ~ x, data = df)
l_ply(split(testdf, testdf$gp), pfun)

In this case, l_ply() plots y vs. x in each of the three subgroups of
testdf separately. It appears you want to do a similar thing, but with
the list of model outputs instead; in your case, the function you
write should take a list as its sole argument. Any slots or components
that are accessed are then relative to the input list object - see the
anonymous function I wrote for the output data frame above as an
illustration.

Since it appears that gevrlevelPlot() also outputs a vector, you may
want to write a function that returns the output vector as a data
frame and call ldply() instead. I don't know enough about the
fExtremes package to write this myself, but it should be possible to
write a function for a generic model that generates both a plot and an
output data frame. The plyr package is pretty good at this sort of
thing.

HTH,
Dennis

On Wed, Jun 29, 2011 at 9:16 PM, Peter Maclean <pmaclean2011 at yahoo.com> wrote:
> I am estimating a large model by groups. How do you save the results 
>and returns
> the associated quantiles?
> For this example I need a data frame
> n    xi        mu        beta
> 1   0.1033614  2.5389580 0.9092611
> 2   0.3401922  0.5192882 1.5290615
> 3   0.5130798  0.5668308 1.2105666
> I also want to apply gevrlevelPlot() for each "n" or group.
>
> #Example
> n <- c(1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3,3)
> y <- c(2,3,2,3,4,5,6,1,0,0,0,6, 2, 1, 0, 0,9,3)
> z <- as.data.frame(cbind(n,y))
> colnames(z) <- c("n","y")
> library(fExtremes)
> z <- split(z, z$n)
> res2 <-lapply(z, function(x){
>                m <- as.numeric(x$y)
>                gevFit(m, block = 1, type = c("pwm"))
>                 })
>> res2
> $`1`
> Title:
>  GEV Parameter Estimation
> Call:
>  gevFit(x = m, block = 1, type = c("pwm"))
> Estimation Type:
>   gev pwm
> Estimated Parameters:
>        xi        mu      beta
> 0.1033614 2.5389580 0.9092611
> Description
>   Wed Jun 29 23:07:48 2011
>
> $`2`
> Title:
>  GEV Parameter Estimation
> Call:
>  gevFit(x = m, block = 1, type = c("pwm"))
> Estimation Type:
>   gev pwm
> Estimated Parameters:
>        xi        mu      beta
> 0.3401922 0.5192882 1.5290615
> Description
>   Wed Jun 29 23:07:48 2011
>
> $`3`
> Title:
>  GEV Parameter Estimation
> Call:
>  gevFit(x = m, block = 1, type = c("pwm"))
> Estimation Type:
>   gev pwm
> Estimated Parameters:
>        xi        mu      beta
> 0.5130798 0.5668308 1.2105666
> Description
>   Wed Jun 29 23:07:48 2011
>
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>