用数据在同一图中绘制 GAM 预测值和 GAMM AR1 时出现问题答案

【问题标题】：Problems plotting GAM predicted values and GAMM AR1 in the same plot with data用数据在同一图中绘制 GAM 预测值和 GAMM AR1 时出现问题
【发布时间】：2019-06-17 15:50:08
【问题描述】：

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

如您所见，GAM 图具有预测值、CI 和带有数据点的“真实数据”x 轴。 GAMM AR1 有一条蓝线，“GAMM 值”位于 x 轴上。

如何绘制来自 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 和数据绘制 pred 值：

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

【问题讨论】：

标签： r plot gam

【解决方案1】：

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

【讨论】：

这行得通，谢谢。关系的AIC是否有所改善还有待观察。如何检查新对象的 AIC？使用 AIC(m.eg) 检查原始 GAM 对象“m.eg”的 AIC，但这不适用于“m.eg.ar1”，不适用于“m.eg.ar1$gam”。
您可以使用AIC(m.eg.ar1$lme) 获得 gamm 的 AIC，但是您无法将此 AIC 与为 gam 模型获得的 AIC 进行比较，因为 gam 是由 REML 拟合的，而 gamm 是由惩罚准似然拟合的……但是，如果自相关是第一个模型中的一个问题，检查您是否使用第二个模型解决了自相关问题以支持使用它而不是第一个模型就足够了。这可以使用以下命令完成acf(resid(m.eg.ar1$gam))。
我必须以这种方式处理 5 个类似的数据系列，以校正自相关。该系列中的四个工作正常，这一切都有意义，但是，但是，当在第五个数据系列上运行 gamm 程序时，我收到此错误消息：recalc.corAR1 中的迭代 1 迭代 2 错误（对象 [[i]]， conLin) : 外部函数调用中的 NA/NaN/Inf (arg 1)
这可能是 nlme 无法拟合相关项。如果不查看数据，我无法判断这是如何发生的。也许你应该问一个新问题来解决这个问题。将第五个数据系列的数据与之前的数据系列之一进行比较可能有助于发现问题。