R用ggplot绘制置信带 [英] R Plotting confidence bands with 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 ）* mp $ wav + rnorm（nrow（mp））* mp $ wav ^ 3 
 
 m01 <-gls（wow_poly（wav，3），data = mp，correlation = corARMA ）（）（））
 
 mp $ fit<  -  as.numeric（拟合（m01））
 
p < -  ggplot（mp，aes（year，wow）） + geom_point（）+ geom_line（aes（year，fit））
p 
  

拟合的数值和数据，我想要一些样式为

  p < -  ggplot（mp，aes（一年，哇））+ geom_point（）+ geom_smooth（）
p 
  

生成的gls模型。

谢谢！

解决方案

require（ggplot2）
require（nlme）

set.seed（101）
mp< -data.frame（year = 1990：2010）
N < - nrow（mp）

mp < - 在（mp，
{
wav -norm（N）* cos（2 * pi * year）+ rnorm（N）* sin（2 * pi * year） + 5
wow -norm（N）* wav + rnorm（N）* wav ^ 3
}）

m01 <-gls（wow〜poly（wav ，3），data = mp，correlation = corARMA（p = 1））


与 m01 $ fitted 相同）

$ p $ 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 而不是年：

pre $ （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（） fits - 参见？makepredictcall 。

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!

解决方案

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.

这篇关于R用ggplot绘制置信带的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！