R用ggplot绘制置信带

``````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
``````

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

推荐答案

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

``````fit <- predict(m01)
``````

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

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

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

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