69

次のような 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 モデルによって生成されたバンドを使用します。

ありがとう！

4

1 に答える 1

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

wav vs ワウ

より高い解像度で予測を行うのは良いことですが、当てはめの結果でこれを行うのは少し難しいですpoly()- を参照してください?makepredictcall。

于 2012-12-25T20:24:52.873 に答える