R用ggplot绘制置信带

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

获取拟合值(与m01$fitted相同)

fit <- predict(m01)

通常我们可以使用predict(...,se.fit=TRUE)之类的东西来获得预测的置信区间，但gls不提供这种功能.我们使用的配方与http://glmm.wikidot.com/faq处显示的配方相似:

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

现在以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))

如果我们以wav而不是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))

用更高的分辨率进行预测是很好的，但用poly()次拟合的结果进行预测有点棘手——参见?makepredictcall.