使用ggplot绘制置信区间带的R图表

Question

我想为使用gls拟合的模型创建置信带，就像这样:

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

但是使用gls模型生成的频带。谢谢！

- Leosar

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

获取拟合值（与 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))

如果我们将图表绘制在WoW而不是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。

- Ben Bolker

---

你如何使最后一个图形的多项式平滑，并真正看起来像一个多项式？ - Pertinax

就像我说的那样，这会很棘手。你必须确保以这样的方式构建模型矩阵，以便它使用原始模型的基础 - 参见?makepredictcall。 - Ben Bolker