2
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

You could try using the group_by() and do() functions in the dplyr package, as follows:

1) Run the group_by() function to define the groups of interest (in your case the groups given by the values taken by the SITE variable)

2) Run the do() function on the object generated by group_by(). The argument of the do() function should reference the process to call, and the process will be run on each group defined by the previous group_by() call.

Here is a simple example on your data to produce a plot by each value of SITE, with the points shown with different col:

library(dplyr)
group_by_SITE = group_by(DF_NEW, SITE)
colors = c("red", "blue")
group_by_SITE %>% plot(AVERAGE ~ YEAR, type="p", pch=21, col=colors[SITE], bg=colors[SITE], data=.)

In your case, the group_by() call would be the same as in my example. The argument of the do() should be a function you define that performs the calculations you do for the non-nested GAM.

Note that the . notation (used in the data= argument of the do() call in my example) refers to the data object that is piped (via the >%> operator) to the do() call from the left (in my example this object is group_by_SITE --which is of classes "grouped_df" "tbl_df" "tbl" "data.frame" as revealed by running class(group_by_SITE); so the object is ALSO a data frame).

Hope this helps!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.