[R] Return level plots

[David Winsemius]

On Sep 21, 2012, at 7:17 AM, MichelleNCSU wrote:

> Hello,
> 
> First of all, let me apologize that my statistics background is modest at
> best.  
> 
> I am doing some extreme value analysis on model output (WRF) which have the
> following dimensions:
> 
> speed(time,lat,lon)

How is this object structured? Are there multiple time layers where speed is measured at lat-lon points are successive times?

> 
> I am trying to fit the GPD (gpd.fit) to each point (time,lat,lon) to get a
> return level plot with values at each grid point.  (Map with return level by
> location.)    

I'm not really sure conducting extreme value analysis is a safe procedure when your stats background is "modest at best".

> 
> Here is some code I tried, following similar structure to languages I'm more
> familiar with, but it isn't working.
> 
> Y = as.matrix(time,lat,lon)

This suggests to me that you come from a different computing universe where the as.matrix() function allows multiple objects (presumably vectors) to be placed "side-by-side" and have a matrix object returned. (In this computing universe that role is played by the cbind function.) Perhaps:

  Y = cbind(time,lat,lon)   # "c" standing for column

> 
> for (t in 1:67)
> + for (j in 1:106)
> + for (i in 1:193)
> + fit(t,j,i)<-gpd.fit(speed(t,j,i), threshold=17,ydat = Y)
> 
> I receive errors at this point, and cant figure out how to get individual
> fits at each grid point.

It's a puzzle to me why you think that passing a single point to a regression function will allow any solution. I did a "??" search and find that there is `gpd.fit` function in the 'ismev' package, but it (like all other regression functions of which I am aware)  appears to take full data objects rather than single points. You would not need to use for-loops to pass the object. Perhaps:

  fit <-gpd.fit(Y, threshold=17,ydat = Y)  
     # Not sure where "speed" entered the picture.
     # as noted before there is ambiguity in the problem statement
     # Or did you do a prior differentiation operation?
     # perhaps a 1/ first difference on time?
   # perhaps ... if gpd.fit follows the usual R conventions
   # return the first point
   i=1,j=1,k=1
   pred.ijk <- predict(fit, data,frame(time=i, lat,=j,lon=k) )

   ? expand.grid  # to cover a range

Your placement of a functional form on the LHS of a formula also suggests recent migration from another statistical universe where assignment is done into functions, i.e. forms using parentheses, rather than the extraction/insertion operators: "[" and  "[<-", to put values into structures with dimensions, like matrices and dataframes. 

I think you really need to do some more self-study of the introductory R material rather than making wild guesses at what "might" work based on experience with Python, Perl, or (less likely in view of that effort to assign into a function) Matlab. 

You should also read the Posting Guide. 

-- 
David Winsemius, MD
Alameda, CA, USA