Confidence interval for the change of a GAM over a period

Question

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

Answer 1

I have three approaches.

1. Go back to the original data and test whether the data near the two time points is significantly different.
2. 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.
3. 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)
```