在具有数据的同一图中绘制 GAM 预测值和 GAMM AR1 的问题

Problems plotting GAM predicted values and GAMM AR1 in the same plot with data

我有 运行 一个 GAM 和我的数据,我正在将 GAM 的预测值与数据点一起绘制在图表中。不同区域有15张相同的图表,其中一些图表存在自相关问题。对于这些,我有 运行 一个 GAMM AR1 模型,现在我想为这些区域绘制类似于其他区域的预测值。下面你会看到两张图,左边的一张是 GAM 的预测值和置信区间,以及真实数据。在右侧,您将看到来自 GAMM AR1 的线。

如您所见,GAM 图具有预测值、CI 和带有数据点的 "real data" x 轴。 GAMM AR1 有一条蓝线,x 轴上有 "GAMM values"。

如何绘制 GAMM AR1 的预测值,类似于我对 GAM 所做的?请参阅下面的数据和脚本。

数据(数据框'eg'):

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 脚本和预测值与 CI 和数据的绘图:

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

预测 gamm 的值与预测 gam 的值几乎相同,所以您几乎成功了。 唯一的区别是 gamm 对象由 gamlme 对象组成,而 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)