Posterior simulation of residuals and ACF for GAMs

Question

I am using R to smooth time series using Generalised Additive Models (GAMs).

A preceding question concerned uncertain serial autocorrelation in the residuals. I was impressed by the diagnostic plots from mvgam in that:

• residuals.mvgam extracts posterior draws of Dunn-Smyth (randomized quantile) residuals, and;
• plot.mvgam(type="residuals") shows the variability in those draws.

The diagnostics from mgcv are more basic (and traditional) in that residuals.gam() returns point estimates.

My specific concern was uncertainty about serial autocorrelation in the residuals. Examining the mvgam residuals, I found that mean(acf(residuals)) is not equal to acf(mean(residuals)), i.e. the acf() for the point estimates mean(residuals) could be insufficient for checking serial correlation. Here are two relevant plots from the previous post:

I thought to apply simulation to the mgcv::gam() model and generate serial autocorrelation diagnostics similar to those from mvgam. Below are two plots of the results and the code is below those.

The ACF plot looks similar to that for the quantile residuals above in that:

• the pattern in the correlations follows that for the point estimates,
• mean(acf(residuals)) is not equal to acf(mean(residuals)), and;
• the 95% confidence (credibility?) intervals for the correlations are wide.

Setting aside the point estimates, the wide confidence intervals suggest that the serial autocorrelation is uncertain and largely within the tolerance interval.

Is this a reasonable (better?) method for interrogating serial autocorrelation in the residuals?

library("mgcv") # mgcv_1.9-1
library("gratia") # gratia_0.9.2
## data
# some annual counts and populations (100,000s)
# the idea is to model the rate = count / 100,000 population
# n = 29 is a short time series
dat <- data.frame(time = 1996:2024,
  count = c(314, 590, 725, 953, 1218, 1688, 2227, 2751,
  3395, 4189, 4865, 5650, 5829, 6664, 7097, 8078, 8720,
  9455, 9773, 10382, 10519, 10576, 10814, 11904, 13176,
  14631, 15075, 16282, 17371),
  pop = c(23.48816, 23.81456, 24.08718, 24.39848, 24.71095,
  25.08078, 25.35108, 25.59605, 25.78745, 26.0291, 26.2812,
  26.66912, 27.13734, 27.61413, 28.0111, 28.35941, 28.73572,
  29.17764, 29.64432, 30.10587, 30.59446, 31.10672, 31.52141,
  31.9123, 32.20657, 32.11941, 32.41207, 33.22432, 34.03657))
## gam model (with negative binomial error distribution)
mod1 <- gam(count ~ offset(log(pop)) + s(time, k=14),
  data = dat, family = nb())
## simulate residuals
nn <- 29
nsims <- 2000
sims <- simulate(mod1, nsim = nsims, data = dat, unconditional=TRUE) # counts
## compute residuals and acf
# Pearson residuals are easy to compute and 
# comparable to deviance and quantile residuals for this example
lmax <- floor(10*log10(nn)) # acf() lag.max
results1 <- array(dim = c(nsims, nn)) # for residuals
results2 <- array(dim = c(nsims, lmax, 2)) # for acf
# for each simulated y
for (i in 1:nsims){
  # compute Pearson residuals = observed - simulated fit
  vars <- sims[, i] + (sims[, i]^2)/6438.266 # nb(mean, theta)
  r <- (dat$count - sims[, i]) / sqrt(vars)
  # store residual
  results1[i, ] <- r
  # compute acf and pacf
  # discard acf = 1 at lag = 0
  results2[i, , 1] <- acf(r, lag.max=lmax, plot=FALSE)$acf[2:(lmax+1)]
  results2[i, , 2] <- pacf(r, lag.max=lmax, plot=FALSE)$acf
}
## check mean(residuals) agrees with the point estimate
rs <- matrix(nrow=nn, ncol=3)
for (i in 1:nn){
  rs[i, 1] <- mean(results1[, i])
  rs[i, 2:3] <- quantile(results1[, i], probs=c(0.025, 0.975))
}
rs <- cbind(dat$time, rs)
rs <- as.data.frame(rs)
names(rs) <- c("time", "mean", "lower", "upper") 
r1 <- residuals(mod1, type="pearson")
# good match
plot(r1, rs$mean)
abline(0, 1)
## summarise and plot acf
acfs <- matrix(nrow=lmax, ncol=3)
for (i in 1:lmax){
  acfs[i, 1] <- mean(results2[, i, 1])
  acfs[i, 2:3] <- quantile(results2[, i, 1], probs=c(0.025, 0.975))
}
acfs <- cbind(dat$time, acfs, acf(r1, plot=F)$acf[2:(lmax+1)])
acfs <- as.data.frame(acfs)
names(acfs) <- c("lag", "mean", "lower", "upper", "acf")
plot(mean ~ lag, acfs, xlab="Lag", ylab="ACF",
    ylim=c(-0.6,0.6), lab=c(14,5,7), pch=19)
  for (i in 1:lmax){
    lines(rep(acfs[i, 1], 2), acfs[i, 3:4])
  }
  with(acfs, points(lag, acf, pch=19, col="red"))
  abline(h=0)
  abline(h=1.96/sqrt(nn), lty=2)
  abline(h=-1.96/sqrt(nn), lty=2)
  legend("topright", pch=19,
    legend=c("mean(acf p residuals)", "acf(mean(p residuals))"), col = c("black", "red"))
   ## summarise and plot pacf
pacfs <- matrix(nrow=lmax, ncol=3)
for (i in 1:lmax){
  pacfs[i, 1] <- mean(results2[, i, 2])
  pacfs[i, 2:3] <- quantile(results2[, i, 2], probs=c(0.025, 0.975))
}
pacfs <- cbind(dat$time, pacfs, pacf(r1, plot=F)$acf)
pacfs <- as.data.frame(pacfs)
names(pacfs) <- c("lag", "mean", "lower", "upper", "acf")
plot(mean ~ lag, pacfs, xlab="Lag", ylab="PACF",
  ylim=c(-0.6,0.6), lab=c(14,5,7), pch=19)
  for (i in 1:lmax){
    lines(rep(pacfs[i, 1], 2), pacfs[i, 3:4])
  }
  with(pacfs, points(lag, acf, pch=19, col="red"))
  abline(h=0)
  abline(h=1.96/sqrt(nn), lty=2)
  abline(h=-1.96/sqrt(nn), lty=2)
  legend("topright", pch=19, legend=c("mean(pacf p residuals)", "pacf(mean(p residuals))"), col = c("black", "red"))

Answer 1

I replied in the other thread but will do so here as well for completeness. In Bayesian inference with MCMC, each posterior draw represents a plausible model configuration. So we can compute expectation values for any function applied to those draws and use them for inference, including taking summaries of those values. So in {mvgam}, the strategy of mean(acf(residual_draws)) is valid and useful. However, there is no guarantee that taking the summary first will return a valid expectation. I'm fact this can give quite misleading results sometimes, which isn't too difficult to show through simulations. But I don't know if the same logic applies here using {mgcv}. Perhaps you should simulate some stochastic autocorrelated series with known AR effects and see what you get with each of the above strategies

Answer 2

I constructed a simple example to further examine the issue that mean(acf(residuals)) is not equal to acf(mean(residuals)). The data is a linear trend with Gaussian AR1 errors.

library(nlme)
## model
set.seed(1234)
nn <- 100
xx <- 1:nn
yy <- 0.5*xx + 0.5*arima.sim(list(ar=0.6), nn)
mod <- lm(yy ~ xx)
## residuals simulation
nsims <- 1000
lmax <- floor(10*log10(nn)) # acf() lag.max
results1 <- array(dim = c(nsims, nn)) # for residuals
results2 <- array(dim = c(nsims, lmax, 2)) # for acf
# for each simulated ys
for (i in 1:nsims){
  ys <- coef(mod)[1] + coef(mod)[2]*xx + rnorm(nn, 0, sigma(mod))
  rr <- yy - ys
  # store residual
  results1[i, ] <- rr
  # compute acf and pacf
  # discard acf = 1 at lag = 0
  results2[i, , 1] <- acf(rr, lag.max=lmax, plot=FALSE)$acf[2:(lmax+1)]
  results2[i, , 2] <- pacf(rr, lag.max=lmax, plot=FALSE)$acf
}

The rest of the code is at the bottom of this post.

Looking at the plots below:

1. The acf(mean(residuals)) lag 1 autocorrelation = 0.69 is close to the true AR1 = 0.6.
2. The mean(acf(residuals)) lag 1 autocorrelation = 0.33 is biased and the confidence interval doesn't include the true value. The mean(acf(residuals)) estimates have generally smaller magnitudes than the point estimates.

Therefore, I have doubts that posterior simulation is useful for interrogating serial autocorrelation in the residuals of a fitted model.

The rest of the code:

## check mean(residuals) at each xx agrees with the point estimates
rs <- matrix(nrow=nn, ncol=3)
for (i in 1:nn){
  rs[i, 1] <- mean(results1[, i])
  rs[i, 2:3] <- quantile(results1[, i], probs=c(0.025, 0.975))
}
rs <- cbind(xx, rs)
rs <- as.data.frame(rs)
names(rs) <- c("xx", "mean", "lower", "upper") 
rr <- residuals(mod, type="pearson")
# good match
plot(rr, rs$mean)
  abline(0, 1)
## summarise and plot acf
png("lmsims.png", width=2*480)
par(mfrow=c(1,2))
  acfs <- matrix(nrow=lmax, ncol=3)
  for (i in 1:lmax){
    acfs[i, 1] <- mean(results2[, i, 1])
    acfs[i, 2:3] <- quantile(results2[, i, 1], probs=c(0.025, 0.975))
  }
  acfs <- cbind(1:lmax, acfs, acf(rr, plot=F)$acf[2:(lmax+1)])
  acfs <- as.data.frame(acfs)
  names(acfs) <- c("lag", "mean", "lower", "upper", "acf")
  plot(mean ~ lag, acfs, xlab="Lag", ylab="ACF",
    ylim=c(-0.4,0.7), lab=c(20,5,7), pch=19)
    for (i in 1:lmax){
      lines(rep(acfs[i, 1], 2), acfs[i, 3:4])
    }
    with(acfs, points(lag, acf, pch=19, col="red"))
    abline(h=0)
    abline(h=1.96/sqrt(nn), lty=2)
    abline(h=-1.96/sqrt(nn), lty=2)
    legend("topright", pch=19,
      legend=c("mean(acf residuals)", "acf(mean(residuals))"), col = c("black", "red"))
  ## summarise and plot pacf
  pacfs <- matrix(nrow=lmax, ncol=3)
  for (i in 1:lmax){
    pacfs[i, 1] <- mean(results2[, i, 2])
    pacfs[i, 2:3] <- quantile(results2[, i, 2], probs=c(0.025, 0.975))
  }
  pacfs <- cbind(1:lmax, pacfs, pacf(rr, plot=F)$acf)
  pacfs <- as.data.frame(pacfs)
  names(pacfs) <- c("lag", "mean", "lower", "upper", "acf")
  plot(mean ~ lag, pacfs, xlab="Lag", ylab="PACF",
    ylim=c(-0.4,0.7), lab=c(20,5,7), pch=19)
    for (i in 1:lmax){
      lines(rep(pacfs[i, 1], 2), pacfs[i, 3:4])
    }
    with(pacfs, points(lag, acf, pch=19, col="red"))
    abline(h=0)
    abline(h=1.96/sqrt(nn), lty=2)
    abline(h=-1.96/sqrt(nn), lty=2)
    legend("topright", pch=19,
      legend=c("mean(pacf p residuals)", "pacf(mean(p residuals))"), col = c("black", "red")) 
par(mfrow=c(1,1))
dev.off()