R Plotting confidence bands with ggplot

Question

I would like to create a confidence band for a model fitted with gls like this:

require(ggplot2) require(nlme)  mp <-data.frame(year=c(1990:2010))  mp$wav <- rnorm(nrow(mp))*cos(2*pi*mp$year)+2*sin(rnorm(nrow(mp)*pi*mp$wav))+5 mp$wow <- rnorm(nrow(mp))*mp$wav+rnorm(nrow(mp))*mp$wav^3  m01 <- gls(wow~poly(wav,3), data=mp, correlation = corARMA(p=1))  mp$fit <- as.numeric(fitted(m01))  p <- ggplot(mp, aes(year, wow))+ geom_point()+ geom_line(aes(year,fit)) p 

This only plots the fitted values and the data, and I would like something in the style of

p <- ggplot(mp, aes(year, wow))+ geom_point()+ geom_smooth() p 

but with the bands generated by the gls model.

Thanks!

Answer 1

require(ggplot2) require(nlme)  set.seed(101) mp <-data.frame(year=1990:2010) N <- nrow(mp)  mp <- within(mp,          {              wav <- rnorm(N)*cos(2*pi*year)+rnorm(N)*sin(2*pi*year)+5              wow <- rnorm(N)*wav+rnorm(N)*wav^3          })  m01 <- gls(wow~poly(wav,3), data=mp, correlation = corARMA(p=1)) 

Get fitted values (the same as m01$fitted)

fit <- predict(m01) 

Normally we could use something like predict(...,se.fit=TRUE) to get the confidence intervals on the prediction, but gls doesn't provide this capability. We use a recipe similar to the one shown at http://glmm.wikidot.com/faq :

V <- vcov(m01) X <- model.matrix(~poly(wav,3),data=mp) se.fit <- sqrt(diag(X %*% V %*% t(X))) 

Put together a "prediction frame":

predframe <- with(mp,data.frame(year,wav,                                 wow=fit,lwr=fit-1.96*se.fit,upr=fit+1.96*se.fit)) 

Now plot with geom_ribbon

(p1 <- ggplot(mp, aes(year, wow))+     geom_point()+     geom_line(data=predframe)+     geom_ribbon(data=predframe,aes(ymin=lwr,ymax=upr),alpha=0.3)) 

It's easier to see that we got the right answer if we plot against wav rather than year:

(p2 <- ggplot(mp, aes(wav, wow))+     geom_point()+     geom_line(data=predframe)+     geom_ribbon(data=predframe,aes(ymin=lwr,ymax=upr),alpha=0.3)) 

It would be nice to do the predictions with more resolution, but it's a little tricky to do this with the results of poly() fits -- see ?makepredictcall.