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

A preceding question concerned [uncertain serial autocorrelation in the residuals][1]. 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:

[![mvgam diagnostics][2]][2]
[![custom acf plot][3]][3]

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?**

[![custom gam acf plot][4]][4]
[![custom gam pacf plot][5]][5]

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

  [1]: https://stats.stackexchange.com/questions/657495/uncertain-serial-autocorrelation-in-gam-count-model-residuals
  [2]: https://i.sstatic.net/8OSMnATK.png
  [3]: https://i.sstatic.net/AUWmXn8J.png
  [4]: https://i.sstatic.net/7oqy2GOe.png
  [5]: https://i.sstatic.net/6AQwpVBM.png