I have annual measurements from 2 sites and want to plot where significant changes in slope occur in each site (in a nested design). I can achieve this using non-nested data based on Gavin Simpson's helpful blog https://www.r-bloggers.com/identifying-periods-of-change-in-time-series-with-gams/.
Non-nested example -
DF <- as.data.frame(seq(from = 1950, to = 2000, by = 1)) colnames(DF) <- "YEAR" DF$AVERAGE <- c(1:20,20,20,20,20,20,20:1,1,1,1,1,2,3) library(mgcv) GAM <- gam(AVERAGE ~ s(YEAR), data=DF)
Then I copied the functions listed here (https://gist.githubusercontent.com/gavinsimpson/e73f011fdaaab4bb5a30/raw/82118ee30c9ef1254795d2ec6d356a664cc138ab/Deriv.R) into R and ran them so they are stored in the memory. Then the following code runs and provides the output that I want:
want <- seq(1, nrow(DF), length.out = 200) pdat <- with(DF, data.frame(YEAR = YEAR[want])) p2 <- predict(GAM, newdata = pdat, type = "terms", se.fit = TRUE) pdat <- transform(pdat, p2 = p2$fit, se2 = p2$se.fit) colnames(pdat)[2:3] <- c("p2","se2") df.res <- df.residual(GAM) crit.t <- qt(0.025, df.res, lower.tail = FALSE) pdat <- transform(pdat, upper = p2 + (crit.t * se2), lower = p2 - (crit.t * se2)) Term <- "YEAR" m2.d <- Deriv(GAM) m2.dci <- confint(m2.d, term = Term) m2.dsig <- signifD(pdat$p2, d = m2.d[[Term]]$deriv, + m2.dci[[Term]]$upper, m2.dci[[Term]]$lower) plot(p2 ~ YEAR, data = pdat, type = "n") lines(p2 ~ YEAR, data = pdat) lines(upper ~ YEAR, data = pdat, lty = "dashed") lines(lower ~ YEAR, data = pdat, lty = "dashed") lines(unlist(m2.dsig$incr) ~ YEAR, data = pdat, col = "blue", lwd = 3) lines(unlist(m2.dsig$decr) ~ YEAR, data = pdat, col = "red", lwd = 3)
But I can not get this to work for nested data, unless I simply run separate GAMs for each of the factor levels (SITE in Example 2) if this is best?
A nested example -
DF_NEW <- as.data.frame(seq(from = 1950, to = 2000, by = 1)) colnames(DF_NEW) <- "YEAR" DF_NEW$AVERAGE <- DF$AVERAGE * -1.5 DF_NEW <- rbind(DF,DF_NEW) DF_NEW$YEAR <- rep(seq(from = 1950, to = 2000, by = 1),times = 2) DF_NEW$SITE <- as.factor(rep(c("A","B"),each = 51)) GAM2 <- gam(AVERAGE ~ s(YEAR, by = SITE), data = DF_NEW)
Thank you in advance for any help / suggestions.