Obrączki z ggplot
Chciałbym stworzyć opaskę dla modelu wyposażonego w gls takiego:
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
To tylko wykresy dopasowanych wartości i danych, a ja chciałbym coś w stylu
p <- ggplot(mp, aes(year, wow))+ geom_point()+ geom_smooth()
p
Ale z pasmami generowanymi przez model gls.
Dzięki! 45
1 answers
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)
Normalnie moglibyśmy użyć czegoś w rodzaju predict(...,se.fit=TRUE)
, aby uzyskać przedziały ufności w przewidywaniach, ale gls
nie zapewnia tej możliwości. Używamy przepisu podobnego do tego pokazanego na http://glmm.wikidot.com/faq:
V <- vcov(m01)
X <- model.matrix(~poly(wav,3),data=mp)
se.fit <- sqrt(diag(X %*% V %*% t(X)))
Ułóż razem "ramkę przewidywania":
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))
Łatwiej zobaczyć, że mamy właściwą odpowiedź, jeśli spiskujemy przeciwko wav
zamiast 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))
Byłoby miło zrobić prognozy z większą rozdzielczością, ale jest to trochę trudne do zrobienia z wynikami poly()
pasuje -- patrz ?makepredictcall
.
61
Author: Ben Bolker,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2013-03-25 16:43:47
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/doraprojects.net/template/agent.layouts/content.php on line 54
2013-03-25 16:43:47