# Thread: Graphics of Confidence Interval using R

1. ## Graphics of Confidence Interval using R

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.

2. 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")``````

3. 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?

4. 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?

5. 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)
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)``````

6. Thank you very much.
The article with the g1 graph is on http://www.sendspace.com/file/q4scs2
Lindemberg

 Tweet

#### Posting Permissions

• You may not post new threads
• You may not post replies
• You may not post attachments
• You may not edit your posts