Проблемы при построении предсказанных значений GAM и AR1 GAMM на одном графике с данными - PullRequest
Новая
       40

Проблемы при построении предсказанных значений GAM и AR1 GAMM на одном графике с данными

2 голосов
/ 17 июня 2019

Я запустил GAM со своими данными, и я отображаю прогнозные значения из GAM на графике вместе с точками данных. Существует 15 одинаковых графиков для разных областей, и для некоторых из них существует проблема автокорреляции. Для этого я запустил модель GAMM AR1, и теперь я хочу построить прогнозные значения для них, аналогичные другим областям. Ниже вы увидите два графика, один слева - это прогнозируемые значения GAM с доверительными интервалами вместе с реальными данными. Справа вы увидите линию от GAMM AR1.

enter image description here

Как вы можете видеть, график GAM имеет предсказанные значения, CI и ось x «реальных данных» с точками данных. GAMM AR1 имеет синюю линию с «значениями GAMM» на оси x.

Как отобразить прогнозируемые значения из GAMM AR1 аналогично тому, что я делаю с GAM? См. Данные и сценарии ниже.

Данные (фрейм данных «например»): 

structure(list(Year = c(1970, 1971, 1972, 1973, 1974, 1975, 1976, 
1977, 1978, 1979, 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987, 
1988, 1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 
1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 
2010, 2011, 2012, 2013, 2014, 2015, 2016), F = c(0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 16, 14, 24, 10, 15, 26, 20, 
15, 19, 13, 18), M = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 27, 40, 35, 39, 19, 30, 42, 42, 39, 56, 50), U = c(100, 
79, 71, 87, 119, 56, 98, 78, 50, 58, 71, 131, 159, 89, 100, 43, 
28, 89, 108, 95, 110, 131, 114, 45, 49, 56, 52, 51, 69, 81, 85, 
60, 54, 46, 54, 57, 1, 5, 8, 5, 0, 1, 1, 5, 8, 2, 0), Tot = c(100, 
79, 71, 87, 119, 56, 98, 78, 50, 58, 71, 131, 159, 89, 100, 43, 
28, 89, 108, 95, 110, 131, 114, 45, 49, 56, 52, 51, 69, 81, 85, 
60, 54, 46, 54, 57, 44, 59, 67, 54, 34, 57, 63, 62, 66, 71, 68
), ratio = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, 1.6875, 2.85714285714286, 1.45833333333333, 
3.9, 1.26666666666667, 1.15384615384615, 2.1, 2.8, 2.05263157894737, 
4.30769230769231, 2.77777777777778), popsize = c(NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA)), row.names = c(NA, 
-47L), class = c("tbl_df", "tbl", "data.frame"))

GAM-скрипт и построение предустановленных значений с КИ и данными: 

m.eg <- gam(Tot~s(Year),family=poisson,data=eg)

YearP=seq(1970,2016,by=1)
meg.pred=predict(m.eg,newdata=data.frame(Year=YearP),type="response",se.fit=T)

plotCI(x=YearP, y=meg.pred$fit,uiw=2*meg.pred$se.fit, type="l",sfrac=0.003,
       ylim=c(40,140),xlim=c(1970,2016),
       col="red",gap=0,lwd=1.6,cex=1.2,las=1, 
       xlab="", ylab="Number of animals")
points(eg$Year,eg$Tot,pch=19,cex=0.9)

GAMM AR1 скрипт и прорисовка линии: 

m.eg.ar1 <- gamm(Tot~s(Year),family=poisson,correlation = corAR1(form = ~ Year), data=eg)
plot(eg$Year,predict(m.eg),col="white")
lines(eg$Year,predict(m.eg.ar1$gam),col="blue")

1 Ответ

2 голосов
/ 20 июня 2019

Прогнозирование значений для gamm почти такое же, как для gam, так что вы почти получили его.Единственное отличие состоит в том, что объект gamm состоит из gam и lme объекта, а метод predict принимает только gam в качестве аргумента.Поэтому вы должны указать это в вызове predict.Код ниже работает для вашего примера: 

m.eg.ar1 <- gamm(Tot~s(Year),family=poisson,correlation = corAR1(form = ~ Year), data=eg)
megar1.pred=predict(m.eg.ar1$gam,newdata=data.frame(Year=YearP),type="response",se.fit=T)
plotCI(x=YearP, y=megar1.pred$fit,uiw=2*megar1.pred$se.fit, type="l",sfrac=0.003,
       ylim=c(40,140),xlim=c(1970,2016),
       col="red",gap=0,lwd=1.6,cex=1.2,las=1, 
       xlab="", ylab="Number of animals")
points(eg$Year,eg$Tot,pch=19,cex=0.9)

Result of gamm with raw data

...