--- title: "GAM" author: "P Wright" date: "14 June 2019" output: html_document --- ```{r} library (mgcv) MyFreqsrd<-read.csv("Timeseries_hhrdkill.csv") MyFreqsrd$date<-as.Date(MyFreqsrd$date, format="%d/%m/%Y") MyFreqsrd$newdate<-format(MyFreqsrd$date, "%Y-%b-%d") m <- gamm(Freq ~ s(Month, bs = "cc", k = 12) + s(Months012000), family=poisson, data = MyFreqsrd) summary(m$gam) layout(matrix(1:2, ncol = 2)) plot(m$gam, scale = 0) layout(1) layout(matrix(1:2, ncol = 2)) acf(resid(m$lme), lag.max = 36, main = "ACF") pacf(resid(m$lme), lag.max = 36, main = "pACF") layout(1) ``` ```{r} ctrl <- list(niterEM = 0, msVerbose = TRUE, optimMethod="L-BFGS-B") ## AR(1) m1 <- gamm(Freq ~ s(Month, bs = "cc", k = 12) + s(Months012000, k = 20), data = MyFreqsrd, correlation = corARMA(form = ~ 1|Year, p = 1), family=poisson, control = ctrl) ## AR(2) m2 <- gamm(Freq ~ s(Month, bs = "cc", k = 12) + s(Months012000, k = 20), data = MyFreqsrd, correlation = corARMA(form = ~ 1|Year, p = 2),family=poisson, control = ctrl) ## AR(3) m3 <- gamm(Freq ~ s(Month, bs = "cc", k = 12) + s(Months012000, k = 20), data = MyFreqsrd, correlation = corARMA(form = ~ 1|Year, p = 3),family=poisson, control = ctrl) ## AR(3) m4 <- gamm(Freq ~ s(Month, bs = "cc", k = 12) + s(Months012000, k = 20), data = MyFreqsrd, correlation = corARMA(form = ~ 1|Year, p = 4),family=poisson, control = ctrl) anova(m$lme, m1$lme, m2$lme, m3$lme,m4$lme) ``` ```{r} layout(matrix(1:2, ncol = 2)) plot(m3$gam, scale = 0) layout(1) summary(m2$gam) layout(matrix(1:2, ncol = 2)) res <- resid(m3$lme, type = "normalized") acf(res, lag.max = 36, main = "ACF - AR(3) errors") pacf(res, lag.max = 36, main = "pACF- AR(3) errors") layout(1) ``` ```{r} want <- seq(1, nrow(MyFreqsrd), length.out = 200) pdat <- with(MyFreqsrd, data.frame(Months012000 = Months012000[want], date = date[want], Month = Month[want])) ## predict trend contributions p <- predict(m$gam, newdata = pdat, type = "terms", se.fit = TRUE) p1 <- predict(m1$gam, newdata = pdat, type = "terms", se.fit = TRUE) p2 <- predict(m2$gam, newdata = pdat, type = "terms", se.fit = TRUE) p3 <- predict(m3$gam, newdata = pdat, type = "terms", se.fit = TRUE) p4 <- predict(m4$gam, newdata = pdat, type = "terms", se.fit = TRUE) ## combine with the predictions data, including fitted and SEs pdat <- transform(pdat, p = p$fit[,2], se = p$se.fit[,2], p1 = p1$fit[,2], se1 = p1$se.fit[,2], p2 = p2$fit[,2], se2 = p2$se.fit[,2], p3 = p3$fit[,2], se3 = p3$se.fit[,2], p4 = p4$fit[,2], se4 = p4$se.fit[,2]) ``` ```{r} op <- par(mar = c(5,4,2,2) + 0.1) ylim <- with(pdat, range(p, p1, p2, p3, p4)) ylim[1] <- floor(ylim[1]) ylim[2] <- ceiling(ylim[2]) plot(Freq - mean(Freq) ~ date, data = MyFreqsrd, type = "n", ylab = "Number of roadkill", ylim = ylim, xlab="Date") lines(p ~ date, data = pdat, col = "black") lines(p1 ~ date, data = pdat, col = "red") lines(p2 ~ date, data = pdat, col = "blue") lines(p3 ~ date, data = pdat, col = "forestgreen", lwd = 1) lines(p4 ~ date, data = pdat, col = "yellow") legend("topleft", legend = c("Uncorrelated Errors", paste0("AR(", 1:4, ") Errors")), bty = "n", col = c("black","red","blue","forestgreen", "yellow"), lty = 1, lwd = c(1,1,1)) par(op) ``` ```{r} summary(m2$gam) layout(matrix(1:2, ncol = 2)) plot(m2$gam, scale = 0) layout(1) library(mgcViz) b <- getViz(m2$gam) o <- plot( sm(b, 1) ) o + l_fitLine(colour = "black") + #l_rug(mapping = aes(x=x, y=y), alpha = 0.8) + l_ciLine(mul = 5, colour = "black", linetype = 2) + scale_x_continuous(name="Month", breaks=c(1,2,3,4,5,6,7,8,9,10,11,12), limits=c(1, 12))+ geom_hline(yintercept = 0, lty=3)+ #l_points(shape = 19, size = 1, alpha = 0.1) + theme_classic() ```