我想为配备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模型产生的谱带.

谢谢

推荐答案

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

年份vs wow

如果我们以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))

wav vs wow

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

R相关问答推荐

根据矢量化/非矢量化调用,使用 stringr 的字符串操作不起作用

如何对所有 NA 为 0/NA 的行求和

如何做每列报告的最新值的新数据框?

用逗号分隔的引号连接字符串

在R中按组查找最长的值序列

将具有不同名称的嵌套列表转换为 data.frame 填充 NA 并添加列

R Shiny:滑块锚在末端重叠

使用 mutate_at 和 contains 将函数应用于多列

模型未能在 r (lme4) 中收敛或不收敛

空格后删除

如何计算R中每一行的每一列的周期总和

应用 ifelse 而不使用 R 应用 for 循环

如何将 R 汇总与多个数字和基于文本的条件子集一起使用

如何为每个条件添加中位数和标准差

R Plotly Bar Chart - 添加水平线标记

使用 dplyr 进行动态编程 - 使用动态输入改变多个动态列

R pivot_longer 带有存根名称和最后一个下划线

根据 dplyr 中的开始和停止日期生成新变量

查找 R 中两个时间戳之间的重叠以分配班次

查找 R 包中使用的 Fortran 文件