Confidence interval for the change of a GAM over a period Similar to Confidence interval for the slope of a GAM, I am fitting a number of gam models to time series data and want to estimate the change (and its uncertainty) over an interval of time (e.g. from x = 15 to 20). I can estimate the instantaneous slope using the derivative() function of the gratis package. However I was wondering about estimating the difference in the predicted model between two points in time (e.g. five years apart). Can I just use a z-test? What about handling correlation?
library(gratia) # https://github.com/gavinsimpson/gratia/
library(mgcv)
library(ggplot2)
library(dplyr)

set.seed(123)

temp3 <- data.frame(z = seq(0, 20, 1/12)) %>%
  mutate(
    x = rlnorm(n(), 1, 1),
    y = 1 + sin(z*2*pi/20) + x/10 + rnorm(n(), 0, 1)
    )

# ggplot() + geom_point(data = temp3, mapping = aes(x = z, y = y))

mod1 <- gam(y ~ s(z), data = temp3, method = "REML") # on time omly
summary(mod1)
# appraise(mod1)
# draw(mod1)
pred1 <- predict(mod1, newdata = temp3, se.fit = TRUE)
resid1 <- mod1$residuals
shap1 <- format(shapiro.test(resid1)$p.value, digits = 2) %>% print()

# mod2 <- gam(y ~ te(z, x), data = temp3, method = "REML") # on time and flow
# mod2 <- gam(y ~ ti(z) + ti(x) + ti(z,x), data = temp3, method = "REML") # on time and flow
mod2 <- gam(y ~ s(z) + s(x), data = temp3, method = "REML") # on time and flow
summary(mod2)
# appraise(mod2)
# draw(mod2)
pred2 <- predict(mod2, newdata = temp3 %>% mutate(x = median(x, na.rm = TRUE)), se = TRUE) # remove flow effect
resid2 <- mod2$residuals
shap2 <- format(shapiro.test(resid2)$p.value, digits = 2) %>% print()

AIC(mod1, mod2)

temp4 <- temp3 %>%
  mutate(
    model1 = pred1$fit,
lower1 = pred1$fit - 1.96 * pred1$se.fit, # 95% C.I.
upper1 = pred1$fit + 1.96 * pred1$se.fit,
model2 = pred2$fit,
    lower2 = pred2$fit - 1.96 * pred2$se.fit, # 95% C.I.
    upper2 = pred2$fit + 1.96 * pred2$se.fit
  )

ggplot(temp4) +
  theme_bw() +
  labs(x = "Year", y = "Concentratoin", colour = "Legend", fill = "") +
  geom_point(mapping = aes(x = z, y = y, colour = "data")) +
  geom_path(mapping = aes(x = z, y = model1, colour = "trend"), size = 2, alpha = 1) +
  geom_ribbon(mapping = aes(x = z, ymin = lower1, ymax = upper1, fill = "trend"), alpha = 0.3) +
  geom_path(mapping = aes(x = z, y = model2, colour = "trendadj"), size = 2, alpha = 1) +
  geom_ribbon(mapping = aes(x = z, ymin = lower2, ymax = upper2, fill = "trendadj"), alpha = 0.3) +
  scale_colour_brewer(palette = "Set1", aesthetics = c("colour", "fill"))


Created on 2021-03-19 by the reprex package (v1.0.0)
 A: I have three approaches.

*

*Go back to the original data and test whether the data near the two time points is significantly different.

*Do what you suggested by using a normal approximation to the predictions.  The downside of this method is that it does not capture the correlation between the predictions at the two time points.

*Bootstrap the experiment so that you can capture the correlation that #2 is missing.  Note that the correlation is high so even small differences are significant

# go back to the sample and test
test_years_apart <- function(yr1, yr2, dat, tol = 0.5)
{
  a <- dat$y[which(dat$z > yr1 - tol & dat$z < yr1 + tol)]
  b <- dat$y[which(dat$z > yr2 - tol & dat$z < yr2 + tol)]
  t.test(a, b)
}
test_years_apart(5, 10, temp3)
test_years_apart(5, 15, temp3)

# assume independent model estimates at two time points
test_years_apart <- function(mod, newdata1, newdata2)
{
  pred1 <- predict(mod, newdata = newdata1, se = TRUE)
  pred2 <- predict(mod, newdata = newdata2, se = TRUE)
  z <- abs(pred1$fit - pred2$fit) / sqrt(pred1$se.fit^2 + pred2$se.fit^2)
  list(interval1 = c(pred1$fit + qnorm(0.025) * pred1$se.fit, pred1$fit + qnorm(0.975) * pred1$se.fit),
       interval2 = c(pred2$fit + qnorm(0.025) * pred2$se.fit, pred2$fit + qnorm(0.975) * pred2$se.fit),
       p.value = 1-pnorm(z))
}
test_years_apart(mod2, 
                 newdata1 = temp3 %>% filter(z == 5) %>% mutate(x = median(x, na.rm = TRUE)),
                 newdata2 = temp3 %>% filter(z == 10) %>% mutate(x = median(x, na.rm = TRUE)))
test_years_apart(mod2, 
                 newdata1 = temp3 %>% filter(z == 5) %>% mutate(x = median(x, na.rm = TRUE)),
                 newdata2 = temp3 %>% filter(z == 15) %>% mutate(x = median(x, na.rm = TRUE)))

# bootstrap to get the correlation right
require(boot)
test_years_apart <- function(yr1, yr2, dat, R = 500)
{
  dat_yr1 <- dat[which(dat$z == yr1)[1],]
  dat_yr2 <- dat[which(dat$z == yr2)[1],]
  b <- boot(dat, statistic = function(d, i){
    mod <- gam(y ~ s(z) + s(x), data = d[i,], method = "REML")
    pred1 <- predict(mod, newdata = dat_yr1)
    pred2 <- predict(mod, newdata = dat_yr2)
    return(c(pred1, pred2, pred1 - pred2))
  }, R = R)
  return(list(interval1 = quantile(b$t[,1], probs = c(0.025, 0.975)),
          interval2 = quantile(b$t[,2], probs = c(0.025, 0.975)),
              p.value = t.test(b$t[,3])$p.value))
}
test_years_apart(5, 10, temp3)
test_years_apart(5, 15, temp3)
```

