autocorrelation-across-lags.Rmd

---
title: "Beyond lag-1 autocorrelation"
author: "Andrew Tredennick"
date: 2019-02-19
output: html_notebook
---

```{r setup, message=FALSE, warning=FALSE}
library(tidyverse)
library(pomp)
knitr::opts_chunk$set(message=FALSE, warning=FALSE, cache = FALSE)
```

# Motivation
For reasons I don't fully understand yet (ask Eamon), lag-1 autocorrelation is not a good EWS for approaching the transcritical bifurcation from above (elimination).
Eamon suggests using his method from the "distance to an epidemic threshold" paper.

## Fitting the oscillating decay model

Function from Eamon.
```{r}
get_fit <- function(y, tstep, 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)))
  }
}
```

Use it on elimination windows.
```{r}
elim_file <- "../../simulations/elimination-simulations-grid-Maradi-1.5e-05.RDS"
all_sims <- readRDS(elim_file)

decays <- tibble()
for(i in unique(all_sims$sim)){
  cases <- all_sims %>%
    filter(sim == i) %>%
   #  mutate(year = trunc(time)) %>%
   #  group_by(year) %>%
   # summarise(reports = sum(reports)) %>%
    pull(reports)
  
  first <- cases[1:(length(cases)/2)]
  second <- cases[((length(cases)/2)+1):length(cases)]
  y1 <- acf(first, lag.max = length(first), plot = FALSE)
  y2 <- acf(second, lag.max = length(second), plot = FALSE)
  
  acf_fits1 <- get_fit(y = as.numeric(y1[[1]]), tstep = 1)
  acf_fits2 <- get_fit(y = as.numeric(y2[[1]]), tstep = 1)
  tmp1 <- acf_fits1$coef["gamma"]
  tmp2 <- acf_fits2$coef["gamma"]
  
  decays %>%
    bind_rows(
      tibble(
        sim = i,
        set = c("null", "test"),
        decay_rate = c(tmp1, tmp2)
      )
    ) -> decays
}

ggplot(decays, aes(x = set, y = decay_rate)) +
  geom_boxplot()

pROC::auc(response = decays$set, predictor = decays$decay_rate)

# 
# acf_func <- function(x, gamma, omega, a){
#   sqrt(1 + a^2) * exp(x * gamma) * sin(x * omega + atan2(1, a))
# }
# 
# d1 <- sqrt(w1["gamma"]^2 + w1["omega"]^2)
# d2 <- sqrt(w2["gamma"]^2 + w2["omega"]^2)
# print(c(d1,d2))
# null <- acf_func(x = seq(0,2,1/52), gamma = w1["gamma"], omega = w1["omega"], a = w1["a"])
# test <- acf_func(x = seq(0,2,1/52), gamma = w2["gamma"], omega = w2["omega"], a = w2["a"])
# plot(null, type = "l")
# lines(test, col = "red")
```

```{r}

ystar <- acf_func(x = seq(0,20,1/52), gamma = -0.07, omega = 2, a = 0)
ystar_noise <- ystar + rnorm(length(ystar), 0, 0)
plot(ystar_noise)
test <- get_fit(ystar_noise, tstep = 1/52)
test$coef
```

```{r, calc-crit-vacc}
calc_R0 <- function(beta, qis, season, eta = (365/8), 
                    nu = 0.05, gamma = (365/5)){
  B <- as.numeric((1 + exp(season %*% qis)) * beta)
  R0 <- (eta*B) / ((eta + nu)*(gamma + nu))
  return(R0)
}

vacc_thresholds <- {}

for(do_city in c("Agadez", "Maradi", "Niamey", "Zinder")){
  
  # Load fitted parameters and pomp model 
  mle_file <- paste0("../../results/initial-mif-lls-", do_city, ".csv")
  mles <- read.csv(mle_file) %>% 
    slice(2:n()) %>%  # ignore first row of storage NAs
    filter(loglik == max(loglik, na.rm = TRUE)) %>%
    dplyr::select(-do_grid, -loglik, -loglik_se)
  
  pomp_file <- paste0("../../code/measles-pomp-object-", do_city, ".RDS")
  fitted_pomp <- readRDS(pomp_file)
  
  # Calculte R0 
  qis <- mles %>%
    dplyr::select(b1, b2, b3, b4, b5, b6) %>%
    as.numeric()
  
  beta <- mles %>%
    pull(beta_mu)
  
  bases <- as_tibble(fitted_pomp@covar) %>%
    dplyr::select(starts_with("x")) %>%
    dplyr::slice(1:365) %>%
    mutate(
      day = 1:365
    ) %>%
    gather(key = base, value = value, -day)
  
  season <- bases %>%
    spread(key = base, value = value) %>%
    dplyr::select(-day) %>%
    as.matrix()
  
  N <- round(mean(fitted_pomp@covar[, "N"]))
  
  R0 <- calc_R0(beta = beta, qis = qis, season = season)
  crit_vacc_cover <- 1 - (1/max(R0))
  vacc_thresholds <- bind_rows(
    vacc_thresholds, 
    tibble(
      city = do_city,
      threshold = crit_vacc_cover
    )
  )
}
```

```{r avg-spectrum}
elim_file <- "../../simulations/elimination-simulations-grid-Maradi-1.5e-05.RDS"
tmpsims <- readRDS(elim_file) %>%
  left_join(vacc_thresholds) %>%
  filter(vacc_coverage < threshold)

window_size <- tmpsims %>%
  filter(vacc_coverage > 0.7 & sim == 1) %>%
  nrow()  # returns number to times in window leading to tcritical
    
tmptimes <- tmpsims %>%
  pull(time) %>%
  unique()  # returns unique observation times across simulations

times_before_vaccine <- tmptimes[tmptimes < 50]
start_ts <- length(times_before_vaccine) - window_size + 1
times_before_vaccine <- times_before_vaccine[start_ts:length(times_before_vaccine)]
times_after_vaccine <- tmptimes[tmptimes >= 50]

simtimes <- c(times_before_vaccine, times_after_vaccine)

ews_time_ids <- tibble(
  time = simtimes,
  half = ifelse(
    time < 50,
    "first",
    "second"
  )
)

# Bandwidth is the full window for each half
window_bandwidth <- length(simtimes)/2

# Merge in `time` and `half` information
tmp_ews_data <- tmpsims %>%
  dplyr::select(city, time, reports, sim, vacc_speed) %>%
  left_join(ews_time_ids, by = "time") %>%
  filter(is.na(half) == FALSE)

# spec_array_null <- array(data = NA, dim = c(500, 2142, 103))
# spec_array_test <- array(data = NA, dim = c(500, 2142, 103))
# for(i in 1:length(unique(tmp_ews_data$sim))){
#   tmp1 <- tmp_ews_data %>%
#     filter(sim == i & half == "first") %>%
#     pull(reports)
#   tmp2 <- tmp_ews_data %>%
#     filter(sim == i & half == "second") %>%
#     pull(reports)
#   
#   spec1 <- t(biwavelet::wt(data.frame(time = 1:length(tmp1),cases = tmp1))$power)
#   spec2 <- t(biwavelet::wt(data.frame(time = 1:length(tmp2),cases = tmp2))$power)
#   spec_array_null[i, , ] <- spec1
#   spec_array_test[i, , ] <- spec2
# }
# 
# spec_mean_null <- colMeans(spec_array_null, dims = 1, na.rm = TRUE)
# spec_mean_test <- colMeans(spec_array_test, dims = 1, na.rm = TRUE)
# par(mfrow = c(1,2))
# image(spec_mean_null, main = "Null")
# image(spec_mean_test, main = "Test")
```

```{r wavelet-ews}

wavelet_median <- function (x){
  y1 <- x
  x1 = seq_along(y1)
  t1 <- data.frame(time=x1, incid=y1)
  
  ## Continuous bias corrected wavelet transform (liu et al 2007)
  wt.t1 <- biwavelet::wt(t1, do.sig=FALSE)
  
  Power <- wt.t1$power.corr
  probfun <- function(x) cumsum(x) / sum(x)
  ppower <- apply(Power, 2, probfun)
  islower50 <- ppower < 0.5
  get_med_ind <- function(x) match(FALSE, x)
  inds <- apply(islower50, 2, get_med_ind)
  median <- wt.t1$scale[inds]
  mean(median)
}

first_mat <- tmp_ews_data %>%
  filter(sim < 200) %>%
  filter(half == "first") %>%
  dplyr::select(time, sim, reports) %>%
  spread(sim, reports) %>%
  dplyr::select(-time) %>%
  as.matrix()
second_mat <- tmp_ews_data %>%
  filter(sim < 200) %>%
  filter(half == "second") %>%
  dplyr::select(time, sim, reports) %>%
  spread(sim, reports) %>%
  dplyr::select(-time) %>%
  as.matrix()

w1 <- apply(first_mat, MARGIN = 2, FUN = wavelet_median)
w2 <- apply(second_mat, MARGIN = 2, FUN = wavelet_median)

power_plot <- tibble(
  first_power = w1,
  second_power = w2
) %>%
  gather()

ggplot(power_plot, aes(x = key, y = value)) +
  geom_boxplot()

pROC::auc(response = power_plot$key, predictor = power_plot$value)

# wave_df <- data.frame(time = 1:length(cases),
#                       cases = cases)
# wt.t1 <- biwavelet::wt(wave_df[2594:5188,])
# Power=wt.t1$power
#   probfun <- function(x) cumsum(x) / sum(x)
#   ppower <- apply(Power, 2, probfun)
#   islower50 <- ppower < 0.5
#   get_med_ind <- function(x) match(FALSE, x)
#   inds <- apply(islower50, 2, get_med_ind)
#   median <- wt.t1$scale[inds]
#   median
# 
#   spectrum(cases[1:(length(cases)/2)], span = c(15,15,15)) -> ttt1
#   spectrum(cases[(length(cases)/2):length(cases)], span = c(15,15,15)) -> ttt2
```

**The Plan**
Calcualte power spectrum in each of null and test interval.
EWS is the average `median` (from above function) in each interval.
Also look at average power spectrum across all simulations to get idea of the "true" spectrum of the model.