1
$\begingroup$

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)

$\endgroup$
1
+50
$\begingroup$

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)
```
$\endgroup$
7
  • $\begingroup$ Thank you so much for taking time to look at my question! I think subsetting the data will change the uncertainty, and give the wrong result. I think it might be possible to do it with the "lpmatrix" option of pred.gam() as described in Simon Wood's book but I haven't understood the details yet. $\endgroup$ – Simon Woodward Mar 25 at 21:54
  • $\begingroup$ I agree. I think the bootstrap option (#3) is the right approach. I looked at the lpmatrix option, but don't think it will get you to the same place as the bootstrap for your specific question. $\endgroup$ – R Carnell Mar 26 at 1:55
  • $\begingroup$ Thanks R, isn't the bootstrap method also subsampling the data though? $\endgroup$ – Simon Woodward Mar 26 at 2:04
  • $\begingroup$ Yes. The bootstrap samples with replacement, but does it multiple times (500 times in my example) so that the whole dataset is covered many times over. The main advantage of the bootstrap is that the predictions at the two time points you want to compare are correlated in a joint distribution. $\endgroup$ – R Carnell Mar 26 at 12:17
  • 1
    $\begingroup$ The bootstrap samples the same amount of data as the original set, with replacement, as you guessed. The boostrap does not not construct a pdf for the difference, but instead think of it as an empirical pdf of the point estimates from each sub-sampled model. $\endgroup$ – R Carnell Mar 29 at 12:29

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.