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

## ----helpers------------------------------------------------------------------
phase_starts <- function(fit) {
  a <- fit$augmented
  a[[fit$metadata$index_name]][!duplicated(a$.phase)][-1]
}
article_axes <- labs(x = "Data", y = "\u00d3bitos di\u00e1rios")

## ----fig2-naive---------------------------------------------------------------
br_naive <- subset(cvd_brazil, region == "BR" &
                     date >= as.Date("2020-03-17") &
                     date <= as.Date("2020-07-04"))
fit_naive <- shewhart_i_mr(br_naive, value = new_deaths, index = date)
fit_naive$augmented |>
  slice(1) |>
  select(.center, .upper, .lower)

## ----fig2-compare, fig.height = 3.6-------------------------------------------
fit_one <- shewhart_regression(br_naive, value = new_deaths, index = date,
                               model = "log", limits_scale = "model",
                               lower_bound = 0,
                               phase_changes = integer(0))
both <- bind_rows(
  fit_naive$augmented |>
    transmute(date, .value, .center, .upper, .lower,
              panel = "Unadapted individuals chart (SBPO Fig. 2)"),
  fit_one$augmented |>
    transmute(date, .value, .center, .upper, .lower,
              panel = "Adapted chart, one phase")
)
both$panel <- factor(both$panel, levels = unique(both$panel))
sig <- shewhart_palette("signal")
ink <- shewhart_palette("neutral")
ggplot(both, aes(date)) +
  geom_ribbon(aes(ymin = .lower, ymax = .upper),
              fill = sig["in_control"], alpha = 0.07) +
  geom_line(aes(y = .upper), colour = sig["out_of_control"],
            linetype = "dashed", linewidth = 0.4) +
  geom_line(aes(y = .lower), colour = sig["out_of_control"],
            linetype = "dashed", linewidth = 0.4) +
  geom_line(aes(y = .center), colour = sig["in_control"], linewidth = 0.7) +
  geom_line(aes(y = .value), colour = ink["text_low"], linewidth = 0.25) +
  geom_point(aes(y = .value), colour = ink["text_high"], size = 0.8) +
  facet_wrap(~panel) +
  coord_cartesian(ylim = c(0, 1800)) +
  labs(x = NULL, y = "Daily deaths",
       title = "Brazil, daily COVID-19 deaths, 17 March to 4 July 2020") +
  shewhart_theme()

## ----fig3-manual, fig.height = 4.6--------------------------------------------
br <- subset(cvd_brazil, region == "BR" &
               date >= as.Date("2020-03-16") & date <= as.Date("2020-07-24"))
br_ends <- as.Date(c("2020-03-27", "2020-04-04", "2020-04-12",
                     "2020-05-16", "2020-05-27", "2020-06-12"))
fit_br <- shewhart_regression(
  br, value = new_deaths, index = date,
  model = "log", limits_scale = "model", lower_bound = 0,
  phase_changes = br_ends + 1,
  rules  = c("nelson_1_beyond_3s", "we_seven_same"),
  locale = "pt"
)
unique(fit_br$augmented$.phase_label)
autoplot(fit_br, phase_dates = TRUE, legend_position = "inside") +
  coord_cartesian(ylim = c(0, 1800)) +
  article_axes

## ----fig3-auto, fig.height = 4.6----------------------------------------------
fit_br_auto <- shewhart_regression(
  br, value = new_deaths, index = date,
  model = "log", limits_scale = "model", lower_bound = 0,
  start_base = 12, phase_rule = "we_seven_same",
  rules  = c("nelson_1_beyond_3s", "we_seven_same"),
  locale = "pt"
)
autoplot(fit_br_auto, legend_position = "inside") +
  coord_cartesian(ylim = c(0, 1800)) +
  article_axes

## ----fig3-table---------------------------------------------------------------
n_max <- max(length(br_ends), length(phase_starts(fit_br_auto)))
knitr::kable(data.frame(
  phase = seq_len(n_max),
  article = format(c(br_ends + 1, rep(NA, n_max - length(br_ends)))),
  automatic = format(c(phase_starts(fit_br_auto),
                       rep(NA, n_max - length(phase_starts(fit_br_auto)))))
), col.names = c("New phase", "Article (first day)", "phase_rule (first day)"))

## ----fig6-auto----------------------------------------------------------------
br_rbe <- subset(cvd_brazil, region == "BR" &
                   date >= as.Date("2020-03-17") &
                   date <= as.Date("2020-06-20"))
fit_rbe <- shewhart_regression(
  br_rbe, value = new_deaths, index = date,
  model = "log", limits_scale = "model", lower_bound = 0,
  start_base = 10, phase_rule = "we_seven_same"
)
phase_starts(fit_rbe)

## ----fig4-manual, fig.height = 4.6--------------------------------------------
rec <- subset(cvd_recife, date >= as.Date("2020-04-30") &
                date <= as.Date("2020-07-24"))
rec_ends <- as.Date(c("2020-05-11", "2020-05-18", "2020-06-04",
                      "2020-06-15", "2020-06-22", "2020-07-04"))
fit_rec <- shewhart_regression(
  rec, value = new_deaths, index = date,
  model = "log", limits_scale = "model", lower_bound = 0,
  phase_changes = rec_ends + 1,
  rules  = c("nelson_1_beyond_3s", "we_seven_same"),
  locale = "pt"
)
unique(fit_rec$augmented$.phase_label)
autoplot(fit_rec, phase_dates = TRUE, legend_position = "inside") +
  coord_cartesian(ylim = c(0, 80)) +
  article_axes

## ----fig4-auto, fig.height = 4.6----------------------------------------------
fit_rec_auto <- shewhart_regression(
  rec, value = new_deaths, index = date,
  model = "log", limits_scale = "model", lower_bound = 0,
  start_base = 12, phase_rule = "we_seven_same",
  rules  = c("nelson_1_beyond_3s", "we_seven_same"),
  locale = "pt"
)
phase_starts(fit_rec_auto)
autoplot(fit_rec_auto, legend_position = "inside") +
  coord_cartesian(ylim = c(0, 80)) +
  article_axes

## ----fig7-pe------------------------------------------------------------------
first_death <- function(d) d[seq(which(d$new_deaths > 0)[1], nrow(d)), ]
pe <- first_death(subset(cvd_brazil, region == "PE" &
                           date <= as.Date("2020-06-20")))
fit_pe <- shewhart_regression(
  pe, value = new_deaths, index = date,
  model = "log", limits_scale = "model", lower_bound = 0,
  start_base = 10, phase_rule = "we_seven_same",
  rules  = c("nelson_1_beyond_3s", "we_seven_same"),
  locale = "pt"
)
autoplot(fit_pe, legend_position = "inside") +
  coord_cartesian(ylim = c(0, 160)) +
  article_axes

## ----fig8-sp------------------------------------------------------------------
sp <- first_death(subset(cvd_brazil, region == "SP" &
                           date <= as.Date("2020-06-20")))
fit_sp <- shewhart_regression(
  sp, value = new_deaths, index = date,
  model = "log", limits_scale = "model", lower_bound = 0,
  start_base = 10, phase_rule = "we_seven_same",
  rules  = c("nelson_1_beyond_3s", "we_seven_same"),
  locale = "pt"
)
autoplot(fit_sp, legend_position = "inside") +
  coord_cartesian(ylim = c(0, 500)) +
  article_axes

## ----phase2-------------------------------------------------------------------
cal <- calibrate(
  subset(rec, date <= as.Date("2020-07-04")),
  chart = "regression", value = new_deaths, index = date,
  model = "log", limits_scale = "model", lower_bound = 0,
  phase_changes = rec_ends[-length(rec_ends)] + 1,
  rules = c("nelson_1_beyond_3s", "we_seven_same")
)
mon <- monitor(subset(rec, date > as.Date("2020-07-04")), cal)
mon$augmented |>
  select(date, .value, .center, .lower, .upper) |>
  head(4)

## ----phase2-plot, fig.height = 4.6--------------------------------------------
viol <- mon$augmented |> filter(.flag_any)
ph <- cal$augmented |>
  group_by(.phase) |>
  mutate(fim = format(max(date))) |>
  ungroup() |>
  mutate(label = if_else(.phase == 0,
                         sprintf("Amostra de base (at\u00e9 %s)", fim),
                         sprintf("Fase %d (at\u00e9 %s)", .phase, fim)))
both <- bind_rows(
  ph |> select(date, .value, .center, .lower, .upper, label),
  mon$augmented |>
    select(date, .value, .center, .lower, .upper) |>
    mutate(label = "Monitoramento")
)
both$label <- factor(both$label, levels = unique(both$label))
pal <- shewhart_palette("phase_seq", n = nlevels(both$label))
ink <- shewhart_palette("neutral")
ggplot(both, aes(date)) +
  geom_ribbon(aes(ymin = .lower, ymax = .upper, fill = label), alpha = 0.07) +
  geom_line(aes(y = .upper, colour = label), linetype = "dashed",
            linewidth = 0.4) +
  geom_line(aes(y = .lower, colour = label), linetype = "dashed",
            linewidth = 0.4) +
  geom_line(aes(y = .center, colour = label), linewidth = 0.7) +
  geom_line(aes(y = .value), colour = ink["text_low"], linewidth = 0.25) +
  geom_point(aes(y = .value, colour = label), size = 1.2) +
  geom_point(data = viol, aes(y = .value), shape = 21, size = 3,
             stroke = 0.8, colour = shewhart_palette("signal")["out_of_control"]) +
  scale_colour_manual(values = pal, name = NULL) +
  scale_fill_manual(values = pal, guide = "none") +
  coord_cartesian(ylim = c(0, 80)) +
  labs(title = "Recife: fases calibradas e monitoramento projetado") +
  article_axes +
  shewhart_theme() +
  theme(legend.position = "inside",
        legend.position.inside = c(0.99, 0.99),
        legend.justification.inside = c(1, 1),
        legend.background = element_rect(fill = ink["bg_panel"], colour = NA))

## ----phase2-signals-----------------------------------------------------------
mon$augmented |>
  filter(.flag_any) |>
  select(date, .value, .center, .lower, .upper,
         .flag_nelson_1_beyond_3s, .flag_we_seven_same)