## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      fig.width = 7, fig.height = 4.2)

## ----setup, message = FALSE---------------------------------------------------
library(shewhartr)
library(ggplot2)
library(dplyr)

## ----replay-frozen------------------------------------------------------------
replay_frozen <- function(d, base = 10, k = 10, rule = "we_seven_same") {
  n <- nrow(d); s <- 1L; size <- base; starts <- d$date[0]
  repeat {
    end <- s + size - 1L
    if (end >= n) break
    cal <- calibrate(d[s:end, ], chart = "regression",
                     value = new_deaths, index = date,
                     model = "log", limits_scale = "model",
                     phase_changes = integer(0), rules = rule)
    mon <- monitor(d[(end + 1L):n, ], cal)
    hit <- which(mon$augmented[[paste0(".flag_", rule)]])
    if (length(hit) == 0L) break
    s <- end + hit[1] + 1L                 # first day after the run
    if (s > n) break
    starts <- c(starts, d$date[s]); size <- k
  }
  starts
}

## ----replay-expanding---------------------------------------------------------
replay_expanding <- function(d, base = 10, k = 10, len = 7L) {
  n <- nrow(d); s <- 1L; size <- base; starts <- d$date[0]
  side <- numeric(0); t <- s + size
  while (t <= n) {
    cal <- calibrate(d[s:(t - 1L), ], chart = "regression",
                     value = new_deaths, index = date,
                     model = "log", limits_scale = "model",
                     phase_changes = integer(0), rules = "we_seven_same")
    today <- monitor(d[t, ], cal)$augmented
    side <- c(side, sign(today$.model_value - today$.model_center))
    m <- length(side)
    if (m >= len && abs(sum(side[(m - len + 1L):m])) == len) {
      s <- t + 1L                          # first day after the run
      if (s > n) break
      starts <- c(starts, d$date[s]); size <- k
      side <- numeric(0); t <- s + size
    } else {
      t <- t + 1L
    }
  }
  starts
}

## ----replay-run---------------------------------------------------------------
rec <- subset(cvd_recife, date >= as.Date("2020-04-30") &
                date <= as.Date("2020-07-24"))
br  <- subset(cvd_brazil, region == "BR" &
                date >= as.Date("2020-03-16") & date <= as.Date("2020-07-24"))

one_shot <- function(d) {
  fit <- shewhart_regression(d, value = new_deaths, index = date,
                             model = "log", limits_scale = "model",
                             start_base = 12, phase_rule = "we_seven_same")
  fit$augmented$date[!duplicated(fit$augmented$.phase)][-1]
}
article <- list(
  Recife = as.Date(c("2020-05-12", "2020-05-19", "2020-06-05",
                     "2020-06-16", "2020-06-23", "2020-07-05")),
  Brazil = as.Date(c("2020-03-28", "2020-04-05", "2020-04-13",
                     "2020-05-17", "2020-05-28", "2020-06-13"))
)
system.time(runs <- list(
  Recife = list(frozen    = replay_frozen(rec, base = 12),
                expanding = replay_expanding(rec, base = 12),
                one_shot  = one_shot(rec)),
  Brazil = list(frozen    = replay_frozen(br, base = 12),
                expanding = replay_expanding(br, base = 12),
                one_shot  = one_shot(br))
))

## ----replay-table-------------------------------------------------------------
show_dates <- function(x, ref) {
  if (length(x) == 0L) return("none")
  gap <- vapply(x, function(z) {
    dd <- as.numeric(z - ref); dd[which.min(abs(dd))]
  }, numeric(1))
  paste0(format(x, "%m-%d"), " (", sprintf("%+d", as.integer(gap)), ")",
         collapse = ", ")
}
tab <- do.call(rbind, lapply(names(article), function(city) {
  data.frame(
    series = city,
    article = paste(format(article[[city]], "%m-%d"), collapse = ", "),
    frozen = show_dates(runs[[city]]$frozen, article[[city]]),
    expanding = show_dates(runs[[city]]$expanding, article[[city]]),
    one_shot = show_dates(runs[[city]]$one_shot, article[[city]])
  )
}))
knitr::kable(tab, col.names = c("Series", "Article (phase starts)",
  "Replay, frozen limits", "Replay, refit daily", "phase_rule (one shot)"))

## ----article-check------------------------------------------------------------
after_end <- function(d, ends) {
  starts <- c(d$date[1], ends + 1)
  vapply(seq_along(ends), function(j) {
    cal <- calibrate(d[d$date >= starts[j] & d$date <= ends[j], ],
                     chart = "regression", value = new_deaths, index = date,
                     model = "log", limits_scale = "model",
                     phase_changes = integer(0), rules = "we_seven_same")
    nxt <- monitor(head(d[d$date > ends[j], ], 10), cal)$augmented
    paste(ifelse(nxt$.model_value > nxt$.model_center, "+", "-"),
          collapse = "")
  }, character(1))
}
data.frame(
  phase_end = format(article$Recife - 1),
  next_10_days = after_end(rec, article$Recife - 1)
)

## ----replay-plot, fig.height = 4.6--------------------------------------------
fit_replay <- shewhart_regression(
  rec, value = new_deaths, index = date,
  model = "log", limits_scale = "model", lower_bound = 0,
  phase_changes = runs$Recife$frozen,
  rules = c("nelson_1_beyond_3s", "we_seven_same"), locale = "pt"
)
autoplot(fit_replay, phase_dates = TRUE, legend_position = "inside") +
  coord_cartesian(ylim = c(0, 80)) +
  labs(x = "Data", y = "\u00d3bitos di\u00e1rios")

## ----fast-replay--------------------------------------------------------------
first_run <- function(side, len) {
  r <- rle(side)
  ok <- which(r$values != 0 & r$lengths >= len)
  if (length(ok) == 0L) return(NA_integer_)
  j <- ok[1]
  cumsum(r$lengths)[j] - r$lengths[j] + len
}
replay_fast <- function(y, base = 10, k = 10, len = 7L,
                        variant = c("frozen", "expanding")) {
  variant <- match.arg(variant)
  g <- log1p(y); n <- length(g); s <- 1L; size <- base
  starts <- integer(0)
  repeat {
    end <- s + size - 1L
    if (end >= n) break
    if (variant == "frozen") {
      b <- stats::lm.fit(cbind(1, seq_len(size)), g[s:end])$coefficients
      proj <- b[1] + b[2] * (((end + 1L):n) - s + 1L)
    } else {
      # fit on days s..t-1 for every t, by cumulative sums
      gg <- g[s:(n - 1L)]; m <- seq_along(gg)
      sy <- cumsum(gg); sxy <- cumsum(m * gg)
      sx <- m * (m + 1) / 2; sxx <- m * (m + 1) * (2 * m + 1) / 6
      b1 <- (m * sxy - sx * sy) / (m * sxx - sx^2)
      b0 <- (sy - b1 * sx) / m
      proj <- (b0 + b1 * (m + 1))[m >= size]
    }
    h <- first_run(sign(g[(end + 1L):n] - proj), len)
    if (is.na(h)) break
    s <- end + h + 1L
    if (s > n) break
    starts <- c(starts, s); size <- k
  }
  starts
}

stopifnot(
  identical(rec$date[replay_fast(rec$new_deaths, base = 12)],
            runs$Recife$frozen),
  identical(br$date[replay_fast(br$new_deaths, base = 12)],
            runs$Brazil$frozen),
  identical(rec$date[replay_fast(rec$new_deaths, base = 12,
                                 variant = "expanding")],
            runs$Recife$expanding),
  identical(br$date[replay_fast(br$new_deaths, base = 12,
                                variant = "expanding")],
            runs$Brazil$expanding),
  identical(rec$date[replay_fast(rec$new_deaths, base = 12, len = 9)],
            replay_frozen(rec, base = 12, rule = "nelson_2_nine_same"))
)

## ----mc-run-------------------------------------------------------------------
n_days <- 120
tau    <- c(41, 76)
mu <- exp(log(5) + 0.05 * (pmin(seq_len(n_days), 40) - 1) -
            0.03 * pmax(seq_len(n_days) - 75, 0))

set.seed(2020)
n_sim <- 1000
Y <- replicate(n_sim, rpois(n_days, mu))

evaluate <- function(signals) {
  delay <- function(j) vapply(signals, function(x) {
    nxt <- if (j < length(tau)) tau[j + 1] else n_days + 1
    hit <- x[x >= tau[j] & x < nxt]
    if (length(hit)) hit[1] - tau[j] + 1 else NA_real_
  }, numeric(1))
  q <- function(x) sprintf("%g [%g, %g]", stats::median(x, na.rm = TRUE),
                           stats::quantile(x, 0.25, na.rm = TRUE),
                           stats::quantile(x, 0.75, na.rm = TRUE))
  d1 <- delay(1); d2 <- delay(2)
  data.frame(
    spurious = sprintf("%.0f%%", 100 * mean(vapply(
      signals, function(x) any(x < tau[1]), logical(1)))),
    delay_1 = q(d1), missed_1 = sprintf("%.0f%%", 100 * mean(is.na(d1))),
    delay_2 = q(d2), missed_2 = sprintf("%.0f%%", 100 * mean(is.na(d2)))
  )
}
grid <- expand.grid(len = c(7L, 9L), variant = c("frozen", "expanding"),
                    stringsAsFactors = FALSE)
mc <- do.call(rbind, lapply(seq_len(nrow(grid)), function(i) {
  sig <- lapply(seq_len(n_sim), function(j) {
    replay_fast(Y[, j], len = grid$len[i], variant = grid$variant[i]) - 1L
  })
  cbind(rule = c(`7` = "we_seven_same", `9` = "nelson_2_nine_same")[
          as.character(grid$len[i])],
        limits = c(frozen = "frozen", expanding = "refit daily")[
          grid$variant[i]],
        evaluate(sig))
}))
knitr::kable(mc, row.names = FALSE, col.names = c(
  "Rule", "Limits", "Spurious phase before day 41",
  "Delay, change 1: median [IQR]", "Missed 1",
  "Delay, change 2: median [IQR]", "Missed 2"))

## ----theory-------------------------------------------------------------------
p_run <- function(n, k) {
  # P(at least one run of k equal signs in n fair +/- signs)
  p <- c(1, rep(0, k - 2))            # current run length 1..k-1
  for (i in seq_len(n - 1)) {
    p <- c(sum(p) / 2, p[-length(p)] / 2)
  }
  1 - sum(p)
}
p_theory <- c(seven = p_run(30, 7), nine = p_run(30, 9))
round(p_theory, 3)

## ----arl0---------------------------------------------------------------------
set.seed(1)
n_ic <- 3000
Y_ic <- replicate(200, rpois(10 + n_ic, exp(log(20) + 0.001 * (0:(9 + n_ic)))))
arl <- expand.grid(len = c(7L, 9L), variant = c("frozen", "expanding"),
                   base = c(10L, 40L), stringsAsFactors = FALSE)
arl[c("ARL0", "no_signal")] <- vapply(seq_len(nrow(arl)), function(i) {
  rl <- apply(Y_ic, 2, function(y) {
    s <- replay_fast(y, base = arl$base[i], len = arl$len[i],
                     variant = arl$variant[i])
    if (length(s)) s[1] - 1 - arl$base[i] else NA_real_
  })
  c(mean(rl, na.rm = TRUE), sum(is.na(rl)))
}, numeric(2)) |> t() |> round()
arl$variant <- ifelse(arl$variant == "expanding", "refit daily", "frozen")
arl$theory <- 2^arl$len - 1
knitr::kable(arl, col.names = c("Run length", "Limits", "Base (days)",
                                "Simulated ARL0 (days)",
                                "Series without a signal", "2^k - 1"))