Determining significant changes in slope in a nested GAM

Question

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.

Answer 1

In case anyone comes across this and needs a solution, the {marginaleffects} package makes this a breeze:

library(mgcv)
#> Loading required package: nlme
#> This is mgcv 1.8-42. For overview type 'help("mgcv-package")'.
library(dplyr)
#> Warning: package 'dplyr' was built under R version 4.3.2
#> 
#> Attaching package: 'dplyr'
#> The following object is masked from 'package:nlme':
#> 
#>     collapse
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(marginaleffects)
library(ggplot2)
#> Warning: package 'ggplot2' was built under R version 4.3.3
theme_set(theme_bw())

# Origingal (un-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)
GAM <- gam(AVERAGE ~ s(YEAR), data=DF)

# Estimated smooth function (on link scale, intercept included)
plot_predictions(GAM, condition = 'YEAR', type = 'link')

# 1st derivative (slope, again on link scale)
plot_slopes(GAM, variables = 'YEAR', by = 'YEAR', type = 'link') +
  geom_hline(yintercept = 0, linetype = 'dashed')

# Where is the 1st derivative not "statistically different from zero"?
# (use a finer sequence of years for more accurate estimates of slopes)
hypotheses(slopes(GAM, newdata = datagrid(YEAR = seq(from = 1950, 
                                                     to = 2000, 
                                                     by = 0.25)),
                                   variables = 'YEAR', 
                            by = 'YEAR', type = 'link')) %>%
  dplyr::filter(p.value > 0.05)
#> 
#>  Term    Contrast YEAR  Estimate Std. Error       z Pr(>|z|)   S   2.5 % 97.5 %
#>  YEAR mean(dY/dX) 1972  0.033529     0.0239  1.4032    0.161 2.6 -0.0133 0.0804
#>  YEAR mean(dY/dX) 1972 -0.014843     0.0239 -0.6215    0.534 0.9 -0.0617 0.0320
#>  YEAR mean(dY/dX) 1996  0.000772     0.0234  0.0329    0.974 0.0 -0.0452 0.0467
#> 
#> Columns: rowid, term, contrast, YEAR, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted 
#> Type:  link

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

# Estimated smooth functions (on link scale, intercept included)
plot_predictions(GAM2, condition = c('YEAR', 'SITE'), 
                 type = 'link')

# 1st derivatives (slopes, again on link scale)
plot_slopes(GAM2, variables = 'YEAR',
            by = c('YEAR', 'SITE'),
            type = 'link') +
  geom_hline(yintercept = 0, linetype = 'dashed')

# Where are 1st derivatives "statistically different from zero"?
# (use a finer sequence of years for more accurate estimates of slopes)
hypotheses(slopes(GAM2, newdata = datagrid(YEAR = seq(from = 1950, 
                                                     to = 2000, 
                                                     by = 0.25)),
                  variables = 'YEAR', 
                  by = c('YEAR', 'SITE'), type = 'link')) %>%
  dplyr::filter(p.value <= 0.05)
#> 
#>  Term    Contrast YEAR SITE Estimate Std. Error     z Pr(>|z|)   S 2.5 %
#>  YEAR mean(dY/dX) 1982    A   -0.649      0.329 -1.97   0.0484 4.4 -1.29
#>  YEAR mean(dY/dX) 1983    A   -0.658      0.329 -2.00   0.0454 4.5 -1.30
#>  YEAR mean(dY/dX) 1983    A   -0.666      0.329 -2.03   0.0428 4.5 -1.31
#>  YEAR mean(dY/dX) 1983    A   -0.675      0.329 -2.05   0.0406 4.6 -1.32
#>  YEAR mean(dY/dX) 1984    A   -0.683      0.330 -2.07   0.0387 4.7 -1.33
#>    97.5 %
#>  -0.00469
#>  -0.01354
#>  -0.02164
#>  -0.02888
#>  -0.03534
#> --- 21 rows omitted. See ?print.marginaleffects --- 
#>  YEAR mean(dY/dX) 1989    A   -0.778      0.377 -2.06   0.0390 4.7 -1.52
#>  YEAR mean(dY/dX) 1989    A   -0.779      0.381 -2.05   0.0407 4.6 -1.52
#>  YEAR mean(dY/dX) 1990    A   -0.780      0.385 -2.03   0.0427 4.6 -1.53
#>  YEAR mean(dY/dX) 1990    A   -0.781      0.389 -2.01   0.0449 4.5 -1.54
#>  YEAR mean(dY/dX) 1990    A   -0.781      0.394 -1.98   0.0475 4.4 -1.55
#>    97.5 %
#>  -0.03918
#>  -0.03302
#>  -0.02583
#>  -0.01770
#>  -0.00863
#> Columns: rowid, term, contrast, YEAR, SITE, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted 
#> Type:  link

You can read more about how these functions work in the highly-detailed marginaleffects docs

Answer 2

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!