I've read various online blogs and forums about obtaining this information, but I cannot seem to get the code to work. I have annual measurements and want to plot where significant changes in slope occur. I have read somewhere that scaling both the response and explanatory variable between 0-1 and examining where the gradient is >-1 or <1 denotes a significant period of changes, and first-order derivatives should be used for this. Can anyone verify this approach? And if so, suggest code how is best to achieve this? Ultimately, I would like to identify and plot (by highlighting in a different colour) periods of significant change, for two types of datasets like those outlined below. I have seen Gavin Simpson's helpful blog on this (https://www.r-bloggers.com/identifying-periods-of-change-in-time-series-with-gams/ ) but I cannot get it to work for my own data.

Example 1) Using a very basic dataset I can construct the following GAM and (I think) obtain first order derivatives. But I am not sure what to do from here.

##Create data and scale between 0 and 1
DF <- as.data.frame(seq(from = 1950, to = 2000, by = 1))
colnames(DF)[1] <- "YEAR"
DF$AVERAGE <- c(1:20,20,20,20,20,20,20:1,1,1,1,1,2,3)
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
DF$YEAR <- range01(DF$YEAR)
DF$AVERAGE <- range01(DF$AVERAGE)

###Create GAM
library(mgcv)
GAM <- gam(AVERAGE ~ s(YEAR), data=DF)
plot(GAM)
newDF <- with(DF, data.frame(YEAR = unique(YEAR)))
X0 <- as.data.frame(predict(GAM, newDF, type = 'lpmatrix'))

Example 2) I'd like to also identify and plot significant changes in slope for two GAM time series whereby values are nested within a factor (SITE here).

DF_NEW <- as.data.frame(seq(from = 1950, to = 2000, by = 1))
colnames(DF_NEW)[1] <- "YEAR"
DF_NEW$AVERAGE <- DF$AVERAGE * -1

DF_NEW <- rbind(DF,DF_NEW)
DF_NEW$SITE <- as.factor(rep(c("A","B"),each = 51))
DF_NEW$YEAR <- range01(DF$YEAR)
DF_NEW$AVERAGE <- range01(DF$AVERAGE)
 
GAM2 <- gam(AVERAGE ~ s(YEAR, by = SITE), data = DF_NEW)

Thank you in advance for any help / suggestions.

UPDATE: I have now figured out how I can do this for a non-nested GAM (i.e. Example 1) based on Gavin Simpson's blog (see above).

I firstly 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 do not know how to run this for nested data (Example 2), unless I simply run separate GAMs for each of the factor levels (SITE in Example 2).

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)[1] <- "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)[1] <- "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.

UPDATE: I have now figured out how I can do this for a non-nested GAM (i.e. Example 1) based on Gavin Simpson's blog (see above).

I firstly 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 do not know how to run this for nested data (Example 2), unless I simply run separate GAMs for each of the factor levels (SITE in Example 2).