autocorrelation for elimination

[atredennick]

Need to add an analysis that implements Eamon's approach for quantifying autocorrelation as outlined in the Distance to epidemic threshold paper. It should serve as the EWS for elimination, not lag-1 autocorrelation as currently presented.

Snippet of code here:

get_fit <- function(y, tstep = 1/52, est_K = FALSE, cutoff = .06) {
  x <- (seq_along(y) - 1) * tstep
  start <- list()
  im <- match(TRUE, abs(y) < cutoff)
  xs <- x[1:(im - 1)]
  ys <- y[1:(im - 1)]
  start$gamma <- try(unname(coef(lm(log(I(abs(ys))) ~ xs))["xs"]))
  if (!inherits(start$gamma, "try-error")){
    spec <- spectrum(y, plot = FALSE, na.action = na.exclude)
    start$omega <- spec$freq[which.max(spec$spec)] / tstep
    start$a <- 0
    fit_osc <- try(minpack.lm::nlsLM(
      y~sqrt(1 + a^2) * exp(x * gamma) * sin(x * omega + atan2(1, a)),
      start = start, na.action = na.exclude,
      control = minpack.lm::nls.lm.control(maxiter = 1000)))
    if (est_K) {
      fit_decay <- try(minpack.lm::nlsLM(y~K * exp(x * gamma),
                                         start = list(gamma = start$gamma, K = y[1]), na.action = na.exclude,
                                         control = minpack.lm::nls.lm.control(maxiter = 1000)))
    } else {
      K <- y[1]
      fit_decay <- try(minpack.lm::nlsLM(y~K * exp(x * gamma),
                                         start = list(gamma = start$gamma),na.action = na.exclude,
                                         control = minpack.lm::nls.lm.control(maxiter = 1000)))
    }
    if (inherits(fit_osc , "try-error")) {
      e_osc <- Inf
    } else {
      e_osc <- fit_osc$m$resid()
    }
    if (inherits(fit_decay, "try-error")) {
      e_decay <- Inf
    } else {
      e_decay <- fit_decay$m$resid()
    }
    nll <- function(resids) {
      n <- length(resids)
      (sum(resids ^ 2))
    }
    aic <- c(constant = nll(y), fit_decay = nll(e_decay) + 2 * (1 + est_K),
             fit_osc = nll(e_osc) + 2 * 3)
    fits <- list(constant = "constant_y=0", fit_decay = fit_decay, fit_osc = fit_osc)
    
    coefests <- try(coef(fits[[which.min(aic)]])[c("omega", "gamma", "a")])
    if (inherits(coefests, "try-error")){
      coefests <- c(NA, NA, NA)
    }
    names(coefests) <- c("omega", "gamma", "a")
    c(list(coef = coefests), fits)
  } else {
    c(list(coef = c("omega" = NA, "gamma" = NA, "a" = NA),
           fits = list(constant = "contant_y=0",
                       fit_decay = NA, fit_osc = NA)))
  }
}

cases <- readRDS(datafile)
cases <- cases %>%
  filter(year > 1994) %>%
  filter(region == "Maradi (City)") %>%
  pull(cases)

y <- acf(cases, lag.max = length(cases)-30, plot = TRUE)
acf_fits <- get_fit(y = as.numeric(y[[1]]))
g <- coef(acf_fits$fit_osc)["gamma"]
w <- coef(acf_fits$fit_osc)["omega"]
distance_to_threshold <- sqrt((g)^2 + (w)^2)

[atredennick]

I've tried implementing Eamon's function above, and the oscillating decay model is never the best model for the ACF. The exponential decay model is always chosen, perhaps because the of the strong signal in the first few lags compared to farther out lags (see example ACF plots below).

@e3bo does this sound right to you? That the oscillating model wouldn't be selected? Just want to make sure I am using your function correctly. Relevant code is in this notebook: docs/notebooks/autocorrelation-across-lags.Rmd (https://github.com/project-aero/measles-ews/blob/master/docs/notebooks/autocorrelation-across-lags.Rmd).

If, in fact, the exponential decay model is the best, should we compare decay rates between intervals near and far from the critical transition? Thus looking at the rate at which autocorrelation drops off with lag rather than just lag-1 autocorrelation.