+ Reply to Thread
Results 1 to 6 of 6

Thread: Graphics of Confidence Interval using R

  1. #1
    Points: 1,525, Level: 22
    Level completed: 25%, Points required for next Level: 75

    Posts
    3
    Thanks
    0
    Thanked 0 Times in 0 Posts

    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. #2
    TS Contributor
    Points: 5,565, Level: 48
    Level completed: 8%, Points required for next Level: 185

    Location
    St Albans, UK
    Posts
    257
    Thanks
    0
    Thanked 7 Times in 5 Posts
    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. #3
    Points: 1,525, Level: 22
    Level completed: 25%, Points required for next Level: 75

    Posts
    3
    Thanks
    0
    Thanked 0 Times in 0 Posts
    Thank you for helping me.
    That's an advance!
    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. #4
    TS Contributor
    Points: 5,565, Level: 48
    Level completed: 8%, Points required for next Level: 185

    Location
    St Albans, UK
    Posts
    257
    Thanks
    0
    Thanked 7 Times in 5 Posts
    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. #5
    TS Contributor
    Points: 5,565, Level: 48
    Level completed: 8%, Points required for next Level: 185

    Location
    St Albans, UK
    Posts
    257
    Thanks
    0
    Thanked 7 Times in 5 Posts
    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)
    rlplot( fit, ci=0.1, add.ci=TRUE)
    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)

  6. #6
    Points: 1,525, Level: 22
    Level completed: 25%, Points required for next Level: 75

    Posts
    3
    Thanks
    0
    Thanked 0 Times in 0 Posts

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

+ Reply to Thread

           




Similar Threads

  1. Replies: 1
    Last Post: 01-20-2011, 11:57 PM
  2. 95% confidence interval
    By obscene in forum Statistics
    Replies: 2
    Last Post: 11-01-2010, 06:48 PM
  3. SPSS graphics problem!!
    By Jka83 in forum Statistics
    Replies: 0
    Last Post: 01-20-2010, 06:58 AM
  4. Confidence Interval
    By wchao_iris in forum Statistics
    Replies: 7
    Last Post: 08-27-2008, 09:23 AM
  5. confidence interval
    By bethsher in forum Statistics
    Replies: 2
    Last Post: 07-14-2008, 05:58 PM

Posting Permissions

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






Advertise on Talk Stats