Graphics of Confidence Interval using R

lindemoliveira

New Member
The idea is to make on R a graph like this one (g1) http://img87.imageshack.us/img87/613/59617892.jpg for any vector of data.
So, I have some difficulties on the statistics and on programming in R.
Just to the moment the nearest graph (g2) that I've found is this one:
...............................................................................................
require(car)
x<-seq(1:10)
y<-c(16,16,48.4,49.9,56,57,57,51,54,48.9)
plot(y~x)
a<-lm(y~x)
qq.plot(a,simulate=TRUE,dist="norm",envelope=0.9)
................................................................................................
On the code, y are annual maximum speeds of winds on a meteorological station in Brazil.
The problem is that on the graph of the link indicated "g1" it exists a straight on full line, and teh legend shows that it corresponds to the fitting on GEV-I or Gumbel.
It comes the following doubt: the straight on "g2" is of a normal distribution (?) 'cause the term dist="norm", so how can I change this straight on code to the correct fitting,
and how to make the vaules on ordinates of the "g2" graph assuming the y values of speed (like "g1") and not a reduced variable?
Thanks for the help.

Mike White

TS Contributor
Is the predict function what you are looking for, as in the following example:
Code:
x<-seq(1:10)
y<-c(46,46,48.4,49.9,56,57,57,51,54,48.9)
a<-lm(y~x)
xx<-data.frame(x=seq(1,10,0.2))
cl<-predict(a, newdata=xx, interval="confidence", level=0.9)
cl<-as.data.frame(cl)
matplot(xx,cl, lty=c(1,2,2), type="l", col=c(1,2,2), ylab="predicted y")
points(x, y, pch=22, bg="white")
legend("bottomright", legend=c("Data", "Fit","90% Confidence limit"), pch=c(22,-1,-1),lty=c(-1, 1,2), col=c(1,1,2), bg="white")

lindemoliveira

New Member
Thank you for helping me.
I've changed the data, to an artificial one made by funcion rgumbel(10,35,8) using the package "evd", for testing:
........................................................................................
x<-seq(1:10)
y<-c(70.68355,69.08129,61.11749,40.69553,29.17408,53.89139,59.96864,32.74423,22.94650,38.24350)
a<-lm(y~x)
xx<-data.frame(x=seq(1,10,0.2))
cl<-predict(a, newdata=xx, interval="confidence", level=0.9)
cl<-as.data.frame(cl)
matplot(xx,cl, lty=c(1,2,2), type="l", col=c(1,2,2), ylab="predicted y")
points(x, y, pch=22, bg="white")
legend("bottomright", legend=c("Data", "Fit","90% Confidence limit"), pch=c(22,-1,-1),lty=c(-1, 1,2), col=c(1,1,2), bg="white")
........................................................................................
In this case, i see that a normal distibution delimits the confidence limit, that was on the objective, although the black line on the midle isn't touched by the points of the data, like on the link i've sent on my first message (g1).
My objective in that the line called "Fit" be the ideal ajust of Gumbel distribution, which rule of distribution with parameters loc = a and scale = b is

G(x) = exp{-exp[-(z-a)/b]}

for all real z, where b > 0.

If you undestood my english and if it's possible, could you help me?

Mike White

TS Contributor
I am not sure that I help you much further with this. The vglm and guplot functions in the VGAM package and the gum.fit function in the ismev package might be useful. It might be helpful if you could provide a reference for the source of the image showing the wind speed versus reduced variate?

Mike White

TS Contributor
The best option I have found is to use the the rlplot function in the extRemes package. This uses a GEV distribution. The gumbel distribution is a special case with shape = 0 and the gev fit results with the gumbel data give a shape value that is close to zero.
Code:
rm(list=ls())
library(evd)
library(extRemes)
set.seed(1)
N<-100
y<-rgumbel(N,35,8)
x<-1:N
datx<-data.frame(x,y)
fit <- gev.fit(y)
head(fit$vals) # shape is non-zero An other option is the gum.diag function in the ismev package but this is not as useful because it produces 4 plots and you cannot change the confidence intervals. It is possible to hack the code to produce just the return level plot with the gum.rl function (as shown below) but changing the confidence intervals would be more complex. Code: rm(list=ls()) library(evd) library(ismev) set.seed(1) N<-100 y<-rgumbel(N,35,8) fit1<-gum.fit(y) gum.diag(fit1) gum.rl(c(fit1$mle,0), fit1$cov, fit1$data)