*Bounty: 50*

*Bounty: 50*

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