Search code examples
rplotgam

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


I have ran a GAM with my data and I am plotting the predicted values from the GAM in a graph together with the data points. There are 15 graphs of the same for different areas, and for some of them there is a autocorrelation problem. For these I have run a GAMM AR1 model, and now I want to plot predicted values for these similar to the other areas. Under you will see two graphs, the one to the left are the predicted values from the GAM with confidence intervals, together with the real data. On the right you will see the line from the GAMM AR1.

enter image description here

As you can see the GAM plot has predicted values, CIs and the "real data" x axis with the data points. The GAMM AR1 has a blue line, with the "GAMM values" on the x axis.

How do I plot predicted values from the GAMM AR1 similar to what I do with the GAM? See data and scripts below.

Data (data frame '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 script and plotting of pred values with CIs and data:

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 script and plotting of line:

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

Solution

  • Predicting the values for gamm is almost the same as for gam so you almost had it. The only difference is that a gamm object consists of a gam and an lme object and the predict method only takes gam as argument. So you have to specify this in the predict call. The code below works for your example:

    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