---
title: COVID-19 en población pediátrica menor de 18 años según registros nacionales en Perú, 2020-2023 
authors:
  - name: Brian Norman Peña-Calero
    affiliation: Laboratorio de Innovación en Salud - UPCH
    roles: analyst
    corresponding: true
bibliography: references.bib
date: now
date-format: "DD MMMM, YYYY"
execute: 
  echo: false
  warning: false
  message: false
  dpi: 300
lang: es
editor_options: 
  chunk_output_type: console
---

```{r}
#| cache: false

library(tidyverse)
library(arrow)
library(duckdb)
library(dbplyr)
library(innovar)
library(apyramid)
library(sf)
library(lmerTest)
library(emmeans)
source("03_functions/index_calculator.R")
source("03_functions/compare_levels_per_month.R")
```

# Disponibilidad de los datos

```{r}
#| cache: false

con <- dbConnect(duckdb(), dbdir = "01_data/raw/covid.duckdb", read_only = TRUE)

vacunas <- tbl(con, "vacunas")
vacunas_18 <- tbl(con, "vacunas_18")
fallecidos_total <- tbl(con, "fallecidos_total")
positivos <- tbl(con, "positivos")
positivos_0723 <- tbl(con, "positivos_0723")
hospitalizados <- tbl(con, "hospitalizados")
hosp_final <- tbl(con, "hosp_final")
poblacion <- tbl(con, "poblacion")
peru_mother <- tbl(con, "peru_mother")
# pobreza <- tbl(con, "pobreza")
data(Peru)
```

```{r}
fallecidos_summary <- fallecidos_total %>% 
  summarise(
    n = n(),
    min = min(fecha_fallecimiento),
    max = max(fecha_fallecimiento)
  ) %>% 
  collect() 

positivos_summary <- positivos %>% 
  summarise(
    n = n(),
    min = min(fecha_resultado),
    max = max(fecha_resultado)
  ) %>% 
  collect() 

hospitalizados_summary <- hospitalizados %>% 
  summarise(
    n = n(),
    min = min(fecha_ingreso_hosp),
    max = max(fecha_ingreso_hosp)
  ) %>% 
  collect() 

vacunas_summary <- vacunas %>% 
  summarise(
    n = n(),
    min = min(fecha_vacunacion),
    max = max(fecha_vacunacion)
  ) %>% 
  collect() %>% 
  mutate(
    across(min:max,
           ymd)
  )

max_positivos_0723 <- positivos_0723 %>% 
  summarise(max = max(fecha_resultado)) %>% 
  pull(max)

Peru <- Peru %>%  
  left_join(
    peru_mother %>% 
      as_tibble()
  )

peru_dep_sf <- Peru %>%
  group_by(dep.code, dep) %>%
  summarise() %>% 
  ungroup()

peru_reg_sf <- Peru %>%
  group_by(region) %>%
  summarise() %>% 
  ungroup()

peru_macro_sf <- Peru %>%
  group_by(macrorregion) %>%
  summarise() %>% 
  ungroup()
```

Actualmente los datos disponibles para el análisis tienen las siguientes situaciones:

1.  Fallecidos COVID-19: Ahora mismo, esta data ([Fallecidos por COVID-19 - \[Ministerio de Salud - MINSA\]](https://www.datosabiertos.gob.pe/dataset/fallecidos-por-covid-19-ministerio-de-salud-minsa)) se encuentra actualizada en el portal de datos abiertos. Cuenta con `r scales::number(fallecidos_summary$n)` fallecidos que van desde `r fallecidos_summary$min` al `r fallecidos_summary$max`.

2.  [Casos positivos COVID-19](https://www.datosabiertos.gob.pe/dataset/casos-positivos-por-covid-19-ministerio-de-salud-minsa): Esta base de datos ya fue corregida y actualizada recientemente en la plataforma de datos abiertos. Cuenta con `r scales::number(positivos_summary$n)` casos positivos que van desde `r positivos_summary$min` al `r positivos_summary$max`. Sin embargo, tiene el inconveniente de no tener registrado el `id_persona` que es indispensable para poder calcular los casos nuevos positivos COVID-19. Se usará por ahora una base actualizada a junio del 2023 con algunos casos menos, pero que si tiene esta información. El enlace de descarga de esa data desde el mismo Minsa ya no se encuentra disponible, pero si está en este repositorio Github: <https://github.com/jmcastagnetto/covid-19-peru-limpiar-datos-minsa/tree/main/datos/originales>

3.  [Hospitalizados por COVID-19](https://www.datosabiertos.gob.pe/dataset/hospitalizados-vacunados-y-fallecidos-por-covid-19): Esta data contiene a la vez información de la vacunación y fallecimiento de quienes fueron hospitalizados. Cuenta con `r scales::number(hospitalizados_summary$n)` casos hospitalizados que van desde `r hospitalizados_summary$min` al `r hospitalizados_summary$max`.

4.  [Vacunación COVID-19](https://www.datosabiertos.gob.pe/dataset/vacunaci%C3%B3n-contra-covid-19-ministerio-de-salud-minsa): Esta data contiene información de vacuna. Cuenta con `r scales::number(vacunas_summary$n)` registros de vacunación que van desde `r vacunas_summary$min` al `r vacunas_summary$max`.

5.  [Población Estimada](https://www.minsa.gob.pe/reunis/data/poblacion_estimada.asp): Esta data contiene las estimaciones poblacionales desde el 2020 al 2023 realizada por REUNIS. Presenta subdivisiones por distrito, edad y sexo. Sin embargo, la versión del año 2021, no tiene subdivisión por sexo, por lo que se usa [la estimación realizada por el INEI](https://www.datosabiertos.gob.pe/dataset/poblaci%C3%B3n-peru) que si lo contiene.

Para términos de uniformidad, temporalmente se restringirán todos los datos a junio de 2023.

# Positivos - Prevalencia COVID-19

Para la estimación de la prevalencia, se obtuvo los datos poblacionales por distrito, sexo y grupo de edad para 2020, 2021, 2022 y 2023. A partir de ello, se calculó una prevalencia promedio mensual (@eq-prevalencia-anual) que luego fue promediada para obtener un único valor por distrito, sexo y grupo de edad. Finalmente, se reporta la prevalencia promedio mensual por 1000 habitantes en anualidades (@eq-prevalencia-anual-promedio) u olas covid, a fin de tener una lectura con menor términos decimales.

\small

$$
\text{Prevalencia promedio mensual}_i = \frac{\text{Promedio de Número de casos en el mes del año } i}{\text{Población en el año } i}
$$ {#eq-prevalencia-anual}

\normalsize

Donde:

-   $\text{Número de casos en el año } i$ es el número de positivos COVID-19 en un año en específico (2020-2023)
-   $\text{Población en el año } i$ es la población que puede ser afectada por COVID-19 en un año en específico (2020-2023)

\small

$$
\text{Prevalencia Promedio mensual por 1000 habitantes} = \left( \frac{\sum_{i=2020}^{2023} \text{Prevalencia promedio mensual}_i}{\text{Número de Años}} \right) \times 1000
$$ {#eq-prevalencia-anual-promedio}

\normalsize

## Vista general

```{r}
positivos_0723 <- positivos_0723 %>% 
  filter(edad < 18, 
         sexo != "No registrado") %>% 
  mutate(
    overall = "overall"
  ) 

range_pos_18_0723 <- positivos_0723 %>% 
  summarise(
    min = min(fecha_resultado),
    max = max(fecha_resultado)
  ) %>% 
  collect() %>% 
  mutate(
    min = format(min, "%d-%m-%Y"),
    max = format(max, "%d-%m-%Y"),
    fecha = paste0("(", min, " a ", max, ")")
  )
```

<!-- **Cálculo de prevalencias** -->

```{r}
levels <- list(
  c("overall"),
  c("grupo_edad"),
  c("sexo"),
  c("ola_covid"),
  c("ola_covid", "grupo_edad"),
  c("ola_covid", "sexo"),
  c("ola_covid", "region"),
  c("ola_covid", "macrorregion"),
  c("departamento"),
  c("departamento", "sexo"),
  c("departamento", "grupo_edad"),
  c("departamento", "sexo", "grupo_edad"),
  c("departamento", "ola_covid"),
  c("distrito"),
  c("distrito", "sexo"),
  c("distrito", "grupo_edad"),
  c("distrito", "cuartil_pobreza"),
  c("grupo_edad", "cuartil_pobreza"),
  c("sexo", "cuartil_pobreza"),
  c("ola_covid", "cuartil_pobreza"),
  c("ola_covid", "cuartil_nbi"),
  c("cuartil_pobreza"),
  c("cuartil_nbi"),
  c("region"),
  c("macrorregion")
)

# Calcular la prevalencia
results_prevalencia <- index_calculator(
  data = positivos_0723,
  denom = poblacion %>% mutate(overall = "overall"), 
  levels = levels, 
  type = "month"
)

results_per_month_prevalencia <- index_calculator(
  data = positivos_0723,
  denom = poblacion %>% mutate(overall = "overall"),
  levels = levels, 
  type = "per_month"
)
```

```{r}
positivos_0723_tbl <- positivos_0723 %>% 
  select(edad, overall, grupo_edad, sexo, ola_covid, region, 
         macrorregion, cuartil_pobreza, cuartil_nbi) %>% 
  collect()

labelled::var_label(positivos_0723_tbl) <- list(
  overall = "Overall",
  grupo_edad = "Age Group",
  sexo = "Sex",
  ola_covid = "COVID-19 wave",
  region = "Natural Regions",
  macrorregion = "Macroregions",
  cuartil_pobreza = "Poverty level",
  cuartil_nbi = "Unsatisfied basic needs level"
)
```

```{r}
#| label: tbl-summ-positivos
#| tbl-cap: !expr paste0("Descriptivos de casos positivos menores de 18 años ", range_pos_18_0723$fecha)

tbl_summ_positivos <- positivos_0723_tbl %>% 
  select(-edad) %>% 
  gtsummary::tbl_summary(
    type = list(
      overall ~ "dichotomous"
    ),
    value = list(overall = "overall"),
    missing_text = "No data"
  ) %>% 
  gtsummary::modify_header(label = "**Variables**") %>%
  gtsummary::bold_labels()

tbl_summ_positivos_c2 <- results_prevalencia$overall$index %>% 
  mutate(label = "Overall") %>% 
  rename(variable = overall) %>% 
  bind_rows(
    results_prevalencia$grupo_edad$index %>% 
      mutate(variable = "grupo_edad") %>% 
      rename(label = grupo_edad),
    results_prevalencia$sexo$index %>% 
      mutate(variable = "sexo") %>% 
      rename(label = sexo),
    results_prevalencia$ola_covid$index %>% 
      mutate(variable = "ola_covid") %>% 
      rename(label = ola_covid),
    results_prevalencia$region$index %>% 
      mutate(variable = "region") %>% 
      rename(label = region),
    results_prevalencia$macrorregion$index %>% 
      mutate(variable = "macrorregion") %>% 
      rename(label = macrorregion),
    results_prevalencia$cuartil_pobreza$index %>% 
      mutate(variable = "cuartil_pobreza") %>% 
      rename(label = cuartil_pobreza),
    results_prevalencia$cuartil_nbi$index %>% 
      mutate(variable = "cuartil_nbi") %>% 
      rename(label = cuartil_nbi)
  ) %>% 
  drop_na() %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  )


tbl_summ_positivos_final <- tbl_summ_positivos %>%
  gtsummary::modify_table_body(
    ~ .x %>% 
      left_join(tbl_summ_positivos_c2) %>% 
      mutate(
        label = case_when(
          variable %in% c("cuartil_pobreza", "cuartil_nbi") &
            label == "1" ~ "1 (0-19.9%)",
          variable %in% c("cuartil_pobreza", "cuartil_nbi") & 
            label == "2" ~ "2 (20-39.9%)",
          variable %in% c("cuartil_pobreza", "cuartil_nbi") & 
            label == "3" ~ "3 (40-59.9%)",
          variable %in% c("cuartil_pobreza", "cuartil_nbi") &
            label == "4" ~ "4 (60-100%)",
          .default = label
        )
      )
  ) %>% 
  gtsummary::modify_header(
    rate ~ "**Average monthly prevalence per 1000 inhabitants**"
  )

tbl_summ_positivos_final
```

```{r}
tbl_summ_positivos_final %>% 
  gtsummary::as_hux_xlsx("02_output/tables/positivos_descriptivos.xlsx")

tbl_summ_positivos_final %>%
  gtsummary::as_flex_table() %>%
  flextable::save_as_docx(path = "02_output/tables/positivos_descriptivos.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/positivos_descriptivos.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/positivos_descriptivos.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

<!-- Cálculo de Comparaciones Múltiples -->


```{r}
comparing_levels <- list(
  c("grupo_edad"),
  c("sexo"),
  c("ola_covid"),
  c("cuartil_pobreza"),
  c("cuartil_nbi"),
  c("region"),
  c("macrorregion")
)

prevalence_comparisons <- compare_levels_per_month(
  results_per_month = results_per_month_prevalencia,
  levels = comparing_levels
)

# Extraer las tablas `cld` como tibbles, renombrando y añadiendo una columna "level"
prevalence_comparisons_cld_pre <- lapply(names(prevalence_comparisons), function(level) {
  cld <- prevalence_comparisons[[level]]$cld
  
  if (!is.null(cld)) {
    cld <- cld |> 
      as_tibble() |> 
      rename(group = 1) |>  # Renombra la primera columna a 'group'
      mutate(level = level) |> # Añade una columna 'level' con el nombre de la lista
      relocate(level) 
  }
  
  cld
}) |> 
  bind_rows()

prevalence_comparisons_cld_pre <- 
  prevalence_comparisons_cld_pre |> 
  mutate(
    across(
      where(is.numeric),
      ~ formattable::digits(.x, 3)
    )
  ) 
```

```{r}
#| label: tbl-compar-incidencia
#| tbl-cap: !expr paste0("Comparaciones múltiples de incidencia", range_pos_18_0723$fecha)
#| 
prevalence_comparisons_cld <- prevalence_comparisons_cld_pre |> 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

prevalence_comparisons_cld
```

```{r}
prevalence_comparisons_cld_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:8) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/prevalencia_comparacion_multiple.xlsx",
                       open = FALSE)

prevalence_comparisons_cld %>%
  gt::gtsave("02_output/tables/prevalencia_comparacion_multiple.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_comparacion_multiple.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_comparacion_multiple.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::


<!-- Comparaciones de incidencia con respecto a olas -->

```{r}
#| label: tbl-summ-positivos-edad-ola
#| tbl-cap: !expr paste0("Descriptivos de casos positivos menores de 18 años por edad y olas COVID", range_pos_18_0723$fecha)

tbl_summ_positivos_edad_ola <- positivos_0723_tbl %>% 
  select(-edad) %>% 
  gtsummary::tbl_summary(
    by = ola_covid,
    type = list(
      overall ~ "dichotomous"
    ),
    value = list(overall = "overall"),
    missing_text = "No data"
  ) %>% 
  gtsummary::modify_header(label = "**Variables**") %>%
  gtsummary::bold_labels()


ola_covid_positivos_gtsumm <- results_prevalencia$ola_covid$index %>% 
  rename(stat = ola_covid) %>% 
  mutate(label = "Overall",
         variable = "overall") %>% 
  bind_rows(
    results_prevalencia$ola_covid_grupo_edad$index %>%
      mutate(variable = "grupo_edad") %>% 
      rename(label = grupo_edad,
             stat = ola_covid),
    results_prevalencia$ola_covid_sexo$index %>% 
      mutate(variable = "sexo") %>% 
      rename(label = sexo,
             stat = ola_covid),
    results_prevalencia$ola_covid_region$index %>% 
      mutate(variable = "region") %>% 
      rename(label = region,
             stat = ola_covid),
    results_prevalencia$ola_covid_macrorregion$index %>% 
      mutate(variable = "macrorregion") %>%
      rename(label = macrorregion,
             stat = ola_covid),
    results_prevalencia$ola_covid_cuartil_pobreza$index %>% 
      mutate(variable = "cuartil_pobreza") %>% 
      rename(label = cuartil_pobreza,
             stat = ola_covid),
    results_prevalencia$ola_covid_cuartil_nbi$index %>% 
      mutate(variable = "cuartil_nbi") %>% 
      rename(label = cuartil_nbi,
             stat = ola_covid)
  ) %>% 
  mutate(
    rate = formattable::digits(rate, 3),
    stat = as.character(stat),
    stat = case_match(
      stat,
      "1 (Mar – Oct 2020)" ~ "stat_1",
      "2 (Nov 2020 – Oct 2021)" ~ "stat_2",
      "3 (Nov 2021 – Apr 2022)" ~ "stat_3",
      "4 (May 2022 – Oct 2022)" ~ "stat_4",
      "5 (Nov 2022 – Jun 2023)" ~ "stat_5"
    )
  ) 

tbl_summ_positivos_edad_ola$table_body <- tbl_summ_positivos_edad_ola$table_body %>%
  mutate(
    across(
      gtsummary::all_stat_cols(stat_0 = FALSE),
      ~ str_remove_all(., "\\s*\\(.*?\\)\\s*")
    )
  ) %>% 
  pivot_longer(
    cols = stat_1:stat_5,
    names_to = "stat",
    values_to = "freq"
  ) %>% 
  left_join(
    ola_covid_positivos_gtsumm
  ) %>% 
  mutate(
    value = str_c(freq, ifelse(is.na(rate), "", str_c(" (", rate, ")")))
  ) %>% 
  select(-c(freq, rate)) %>% 
  pivot_wider(
    names_from = stat,
    values_from = value
  )


tbl_summ_positivos_edad_ola <- tbl_summ_positivos_edad_ola %>% 
  gtsummary::modify_footnote(
    gtsummary::all_stat_cols() ~ "n (incidence)"
  )

tbl_summ_positivos_edad_ola
```

```{r}
tbl_summ_positivos_edad_ola %>% 
  gtsummary::as_hux_xlsx("02_output/tables/positivos_olas_descriptivos.xlsx")

tbl_summ_positivos_edad_ola %>%
  gtsummary::as_flex_table() %>%
  flextable::save_as_docx(path = "02_output/tables/positivos_olas_descriptivos.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/positivos_olas_descriptivos.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/positivos_olas_descriptivos.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Gráfico de positivos (serie de tiempo) - Menores de 18 años

```{r}
#| fig-cap: !expr paste0("Positividad COVID Población menor de edad de acuerdo al grupo de edad y ola ", range_pos_18_0723$fecha)
#| label: fig-pos-edad-olas-area
#| fig-width: 10

# 
# range_positivos_und_18 <- positivos %>% 
#   filter(edad < 18, 
#          sexo != "No registrado") %>% 
#   pull(fecha_resultado) %>%
#   range()

Sys.setlocale("LC_TIME", "English")

positive_cases_area <- positivos_0723 %>% 
  mutate(
    fecha_round = case_when(
      fecha_round == as.Date("2024-01-01") ~  as.Date("2023-12-31"),
      .default = fecha_round
    )
  ) %>% 
  count(grupo_edad, fecha_round) %>% 
  collect() %>% 
  ggplot(
    aes(
      x = fecha_round,
      y = n,
      fill = grupo_edad
    )
  ) +
  geom_area(position = 'stack') +
  # innovar::scale_fill_innova("jama") +
  scale_fill_manual(
    values = c("#2C9BB4", "#0B4E60", "#DAE2BC", "#8C3A33")
  ) +
  scale_x_date(
    limits = c(ymd("2020-03-01"), ymd("2023-07-01")),
    breaks = seq(from = ceiling_date(ymd("2020-03-01"), "month"),
                 to = floor_date(ymd("2023-07-01"), "month"),
                 by = "2 months"), 
    # date_breaks = "2 month",
    date_labels = "%b %Y"
  ) +
  labs(
    y = NULL,
    x = NULL,
    title = "COVID-19 cases",
    fill = "Age Group"
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2020-03-01"), 
    xmax = as.Date("2020-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2020-07-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "First wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2020-11-01"), 
    xmax = as.Date("2021-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2021-05-10"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Second wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2021-11-01"), 
    xmax = as.Date("2022-04-30"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
    annotate(
    geom = "text", 
    x = as.Date("2022-02-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Third \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2022-05-01"), 
    xmax = as.Date("2022-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-08-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Fourth \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2022-11-01"), 
    xmax = as.Date("2023-07-01"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2023-03-01"),
    y = Inf, # Posición vertical del texto
    label = "Fifth wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  theme_classic() +
  theme(
    axis.text.x = element_text(
      angle = 75,
      vjust = 0.5
    ),
    axis.ticks.x = element_line()
  ) 

positive_cases_area
```

## Gráfico de Incidencia (serie de tiempo) - Menores de 18 años

```{r}
#| fig-cap: !expr paste0("Incidencia COVID Población menor de edad de acuerdo al grupo de edad y ola ", range_pos_18_0723$fecha)
#| label: fig-inci-edad-olas-area
#| fig-width: 10

prevalence_by_age_area <- results_per_month_prevalencia$grupo_edad$index %>% 
  mutate(
    time = make_date(anio, mes, "1")
  ) %>% 
  ggplot(
    aes(
      x = time,
      y = rate,
      fill = grupo_edad
    )
  ) +
  geom_area(position = 'stack') +
  scale_fill_manual(
    values = c("#2C9BB4", "#0B4E60", "#DAE2BC", "#8C3A33")
  ) +
  scale_x_date(
    limits = c(ymd("2020-03-01"), ymd("2023-07-01")),
    breaks = seq(from = ceiling_date(ymd("2020-03-01"), "month"),
                 to = floor_date(ymd("2023-07-01"), "month"),
                 by = "2 months"), 
    # date_breaks = "2 month",
    date_labels = "%b %Y"
  ) +
  labs(
    y = NULL,
    x = NULL,
    title = "Monthly COVID-19 incidence per 1 000 inhabitants",
    fill = "Age Group"
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2020-03-01"), 
    xmax = as.Date("2020-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2020-07-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "First wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2020-11-01"), 
    xmax = as.Date("2021-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2021-05-10"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Second wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2021-11-01"), 
    xmax = as.Date("2022-04-30"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-02-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Third \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2022-05-01"), 
    xmax = as.Date("2022-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-08-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Fourth \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2022-11-01"), 
    xmax = as.Date("2023-07-01"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2023-03-01"),
    y = Inf, # Posición vertical del texto
    label = "Fifth wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  theme_classic() +
  theme(
    axis.text.x = element_text(
      angle = 75,
      vjust = 0.5
    ),
    axis.ticks.x = element_line()
  ) 

prevalence_by_age_area
```

## Pirámide frecuencia - Edades

```{r}
#| fig-width: 8

positivos_0723_tbl %>% 
  mutate(
    edad = factor(edad)
  ) %>% 
  age_pyramid(
    age_group = "edad",
    split_by = "sexo"
  ) +
  scale_fill_brewer(palette = "Dark2") +
  labs(
    x = "Age",
    fill = "Sex",
    y = "Frequency"
  ) +
  theme_classic()
```

```{r}
#| fig-width: 6

positivos_0723_tbl %>% 
  age_pyramid(
    age_group = "grupo_edad",
    split_by = "sexo"
  ) +
  scale_fill_brewer(palette = "Dark2") +
  labs(
    x = "Age Group",
    fill = "Sex",
    y = "Frequency"
  ) +
  scale_y_continuous(
    labels = ~ scales::number(abs(.x), big.mark = ","),
    n.breaks = 10
  ) +
  theme_classic()
```

## Pirámide proporcional - Edades

```{r}
#| fig-width: 8

positivos_0723_tbl %>% 
  mutate(
    edad = factor(edad)
  ) %>% 
  age_pyramid(
    age_group = "edad",
    split_by = "sexo",
    proportional = TRUE
  ) +
   scale_fill_brewer(palette = "Dark2") +
  labs(
    x = "Age",
    fill = "Sex",
    y = "Percentage"
  ) +
  theme_classic()
```

```{r}
#| fig-width: 6

positivos_0723_tbl %>% 
  age_pyramid(
    age_group = "grupo_edad",
    split_by = "sexo",
    proportional = TRUE
  ) +
   scale_fill_brewer(palette = "Dark2") +
  labs(
    x = "Age Group",
    fill = "Sex",
    y = "Percentage"
  ) +
  theme_classic()
```

## Departamento

```{r}
prev_dep_1k_sf <- peru_dep_sf %>% 
  rename(departamento = dep) %>% 
  inner_join(
    results_prevalencia$departamento$index
  )
```

```{r}
#| fig-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al departamento ", range_pos_18_0723$fecha)
#| label: fig-prev-dep-map
#| fig-width: 5
#| fig-height: 5

prev_dep_1k_map <- prev_dep_1k_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap("Mean monthly COVID-19 incidence", 
               20)
    )
  ) +
  scale_fill_distiller(palette = "Oranges", # old: YlOrBr
                       direction = 1,
                       guide = "legend") + 
  theme_void()

prev_dep_1k_map
```

```{r}
#| tbl-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al departamento ", range_pos_18_0723$fecha)
#| label: tbl-prev-dep

tbl_summ_prev_dep_pre <- results_prevalencia$departamento$index %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  ) %>% 
  mutate(
    departamento = str_to_title(departamento)
  ) 

tbl_summ_prev_dep <- tbl_summ_prev_dep_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

tbl_summ_prev_dep
```

```{r}
tbl_summ_prev_dep_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:2) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/prevalencia_departamento.xlsx",
                       open = FALSE)

tbl_summ_prev_dep %>%
  gt::gtsave("02_output/tables/prevalencia_departamento.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_departamento.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_departamento.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Departamento y sexo

```{r}
#| fig-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al departamento y sexo ", range_pos_18_0723$fecha)
#| label: fig-prev-dep-sex-bar

results_prevalencia$departamento_sexo$index %>% 
  mutate(
    departamento = str_to_title(departamento),
    departamento = fct_reorder2(departamento, sexo, rate)
  ) %>% 
  ggplot(aes(x = departamento, y = rate, fill = sexo)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  theme_bw() +
  labs(title = "Monthly incidence of COVID-19 by department and sex", 
       x = NULL, 
       y = "Mean monthly COVID-19 incidence",
       fill = "Sex",
       caption = paste0("Total cases reported: ", 
                        scales::number(sum(pull(results_prevalencia$departamento_sexo$index_pre,
                                                count))))) +
  innovar::scale_fill_innova("blue_fall") +
  theme(
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1
    )
  )
```

```{r}
prev_dep_sex_1k_sf <- peru_dep_sf %>% 
  rename(departamento = dep) %>% 
  inner_join(
    results_prevalencia$departamento_sexo$index
  ) %>% 
  ungroup()
```

```{r}
#| fig-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al departamento y sexo ", range_pos_18_0723$fecha)
#| label: fig-prev-dep-sex-map

prev_dep_sex_1k_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap(
        "Mean monthly COVID-19 incidence",
        20
      )
    )
  ) +
  scale_fill_distiller(palette = "Oranges", # old: YlOrBr
                       direction = 1) + 
  facet_wrap(vars(sexo)) + 
  theme_void()
```

```{r}
#| tbl-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al departamento y sexo ", range_pos_18_0723$fecha)
#| label: tbl-prev-dep-sex

tbl_summ_prev_dep_sex_pre <- results_prevalencia$departamento_sexo$index %>% 
  mutate(
    rate = formattable::digits(rate, 3),
    departamento = str_to_title(departamento)
  ) %>% 
  pivot_wider(
    names_from = sexo, 
    values_from = rate
  )

tbl_summ_prev_dep_sex <- tbl_summ_prev_dep_sex_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

tbl_summ_prev_dep_sex
```

```{r}
tbl_summ_prev_dep_sex_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:3) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/prevalencia_departamento_sexo.xlsx",
                       open = FALSE)

tbl_summ_prev_dep_sex %>%
  gt::gtsave("02_output/tables/prevalencia_departamento_sexo.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_departamento_sexo.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_departamento_sexo.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Departamento y edad

```{r}
#| fig-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al departamento y edad ", range_pos_18_0723$fecha)
#| label: fig-prev-dep-edad-bar
#| fig-width: 11
#| fig-height: 9

results_prevalencia$departamento_grupo_edad$index %>% 
  mutate(
    departamento = str_to_title(departamento),
    departamento = fct_reorder2(departamento, grupo_edad, rate)
  ) %>% 
  ggplot(aes(x = departamento, y = rate, fill = grupo_edad)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  labs(title = "Monthly incidence of COVID-19 by department and age groups", 
       x = NULL, 
       y = "Mean monthly COVID-19 incidence",
       fill = "Sex",
       caption = paste0("Total cases reported: ", 
                        scales::number(sum(pull(results_prevalencia$departamento_grupo_edad$index_pre,
                                                count))))) +
  innovar::scale_fill_innova("blue_fall") +
  facet_wrap(vars(grupo_edad)) + 
  theme(
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1
    )
  )
```

```{r}
prev_dep_edad_1k_sf <- peru_dep_sf %>% 
  rename(departamento = dep) %>% 
  inner_join(
    results_prevalencia$departamento_grupo_edad$index
  ) %>% 
  ungroup()
```

```{r}
#| fig-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al departamento y edad ", range_pos_18_0723$fecha)
#| label: fig-prev-dep-edad-map

prev_dep_edad_1k_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap(
        "Mean monthly COVID-19 incidence",
        20
      )
    )
  ) +
  scale_fill_distiller(palette = "Oranges", # old: YlOrBr
                       direction = 1) + 
  facet_wrap(vars(grupo_edad), nrow = 1) + 
  theme_void()
```

```{r}
#| tbl-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al departamento y grupo de edad ", range_pos_18_0723$fecha)
#| label: tbl-prev-dep-edad

tbl_summ_prev_dep_edad_pre <- results_prevalencia$departamento_grupo_edad$index %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  ) %>% 
  mutate(
    departamento = str_to_title(departamento)
  ) %>% 
  pivot_wider(
    names_from = grupo_edad, 
    values_from = rate
  )

tbl_summ_prev_dep_edad <- tbl_summ_prev_dep_edad_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

tbl_summ_prev_dep_edad
```

```{r}
tbl_summ_prev_dep_edad_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:5) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/prevalencia_departamento_edad.xlsx",
                       open = FALSE)

tbl_summ_prev_dep_edad %>%
  gt::gtsave("02_output/tables/prevalencia_departamento_edad.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_departamento_edad.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_departamento_edad.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Departamento, sexo y grupo de edad

```{r}
#| tbl-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al departamento, edad y sexo ", range_pos_18_0723$fecha)
#| label: tbl-prev-dep-edad-sex

tbl_summ_prev_dep_edad_sex_pre <- results_prevalencia$departamento_sexo_grupo_edad$index %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  ) %>% 
  pivot_wider(
    names_from = grupo_edad, 
    values_from = rate
  ) %>% 
  arrange(departamento, sexo) %>% 
  collect() %>% 
  mutate(
    departamento = str_to_title(departamento)
  ) 


tbl_summ_prev_dep_edad_sex <- tbl_summ_prev_dep_edad_sex_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )
  # gt::cols_label(
  #   departamento = gt::md("**Departamento**"),
  #   sexo = gt::md("**Sexo**")
  # ) 
  # gt::tab_style(
  #   style = gt::cell_text(weight = "bold"),
  #   locations = gt::cells_column_labels()
  # )

tbl_summ_prev_dep_edad_sex
```

```{r}
tbl_summ_prev_dep_edad_sex_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:6) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/prevalencia_departamento_edad_sexo.xlsx",
                       open = FALSE)

tbl_summ_prev_dep_edad_sex %>%
  gt::gtsave("02_output/tables/prevalencia_departamento_edad_sexo.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_departamento_edad_sexo.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_departamento_edad_sexo.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

```{r}
#| fig-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al departamento, sexo y edad ", range_pos_18_0723$fecha)
#| label: fig-prev-dep-sex-edad-area
#| fig-width: 12
#| fig-height: 7

results_prevalencia$departamento_sexo_grupo_edad$index %>% 
  mutate(
    departamento = str_to_title(departamento)
  ) %>% 
  ggplot(
    aes(
      x = departamento,
      y = rate,
      fill = sexo,
      group = sexo
    )
  ) +
  geom_area() +
  facet_wrap(vars(grupo_edad)) +
  innovar::scale_fill_innova("blue_fall") +
  labs(
    y = "Mean monthly COVID-19 incidence",
    x = NULL,
    title = "Monthly incidence of COVID-19 by age groups and sex",
    fill = "Sex"
  ) + 
  theme_bw() +
  theme(
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1
    )
  )
```

## Departamento y olas

```{r}
#| fig-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al departamento y olas covid ", range_pos_18_0723$fecha)
#| label: fig-prev-dep-olas-bar

results_prevalencia$departamento_ola_covid$index %>% 
  drop_na(ola_covid) %>% 
  mutate(
    departamento = str_to_title(departamento),
    departamento = fct_reorder2(departamento, ola_covid, rate)
  ) %>% 
  ggplot(aes(x = departamento, y = rate, fill = ola_covid)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  labs(title = "Monthly incidence of COVID-19 by department and waves", 
       x = NULL, 
       y = "Mean monthly COVID-19 incidence",
       fill = "Sexo",
       caption = paste0("Total cases reported: ", 
                        scales::number(sum(pull(results_prevalencia$departamento_ola_covid$index_pre,
                                                count))))) +
  innovar::scale_fill_innova("blue_fall") +
  facet_wrap(vars(ola_covid)) + 
  theme(
    axis.text.x = element_text(
      angle = 75,
      hjust = 1,
      vjust = 1
    ),
    legend.position = "bottom"
  )
```

```{r}
prev_dep_olas_100k <- results_prevalencia$departamento_ola_covid$index %>% 
  drop_na(ola_covid) %>% 
  bind_rows(
    results_prevalencia$departamento$index %>% 
      mutate(ola_covid = "Overall")
  ) %>% 
  mutate(
    ola_covid = fct_recode(
      ola_covid,
      "1st wave" = "1 (Mar – Oct 2020)",
      "2nd wave" = "2 (Nov 2020 – Oct 2021)",
      "3rd wave" = "3 (Nov 2021 – Apr 2022)",
      "4th wave" = "4 (May 2022 – Oct 2022)",
      "5th wave" = "5 (Nov 2022 – Jun 2023)"
    ),
    ola_covid = fct_relevel(ola_covid,
                            "Overall")
  ) %>% 
  arrange(ola_covid)


prev_mes_dep_olas_1k_sf <- peru_dep_sf %>%
  rename(departamento = dep) %>% 
  left_join(
    prev_dep_olas_100k
  ) %>%
  ungroup()
```

```{r}
#| fig-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al departamento por Olas ", range_pos_18_0723$fecha)
#| label: fig-prev-dep-olas-map
#| fig-width: 13
#| fig-height: 3

prev_mes_dep_olas_1k_map <- prev_mes_dep_olas_1k_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(str_wrap("Mean monthly COVID-19 incidence", 20))
  ) +
  scale_fill_distiller(palette = "Oranges", # Old: YlOrBr
                       direction = 1) + 
  facet_wrap(vars(ola_covid), nrow = 1) + 
  theme_void()

prev_mes_dep_olas_1k_map
```

```{r}
#| tbl-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al departamento y Olas COVID ", range_pos_18_0723$fecha)
#| label: tbl-prev-dep-olas

tbl_summ_prev_dep_ola_pre <- results_prevalencia$departamento_ola_covid$index %>% 
  drop_na(ola_covid) %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  ) %>% 
  mutate(
    departamento = str_to_title(departamento)
  ) %>% 
  pivot_wider(
    names_from = ola_covid, 
    values_from = rate
  )

tbl_summ_prev_dep_ola <- tbl_summ_prev_dep_ola_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

tbl_summ_prev_dep_ola
```

```{r}
tbl_summ_prev_dep_ola_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:6) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/prevalencia_departamento_ola.xlsx",
                       open = FALSE)

tbl_summ_prev_dep_ola %>%
  gt::gtsave("02_output/tables/prevalencia_departamento_ola.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_departamento_ola.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_departamento_ola.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Regiones

```{r}
#| fig-cap: !expr paste0("Prevalencia COVID-19 de acuerdo a la región ", range_pos_18_0723$fecha)
#| label: fig-prev-reg-bar

results_prevalencia$region$index %>% 
  ggplot(aes(x = region, y = rate, fill = region)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  labs(title = "Monthly incidence of COVID-19 by natural regions", 
       x = NULL, 
       y = "Mean monthly COVID-19 incidence",
       fill = "Sex",
       caption = paste0("Total cases reported: ", 
                        scales::number(sum(pull(results_prevalencia$region$index_pre,
                                                count))))) +
  innovar::scale_fill_innova("blue_fall") +
  theme(
    legend.position = "bottom"
  )
```

```{r}
prev_reg_1k_sf <- peru_reg_sf %>% 
  inner_join(
    results_prevalencia$region$index
  )
```

```{r}
#| fig-cap: !expr paste0("Prevalencia COVID-19 de acuerdo a la región ", range_pos_18_0723$fecha)
#| label: fig-prev-reg-map
#| fig-width: 5
#| fig-height: 5

prev_reg_1k_map <- prev_reg_1k_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap("Mean monthly COVID-19 incidence", 
               20)
    )
  ) +
  ggsflabel::geom_sf_label(aes(label = region)) +
  scale_fill_distiller(palette = "Oranges", # YlOrBr
                       direction = 1) + 
  theme_void()

prev_reg_1k_map
```

```{r}
#| tbl-cap: !expr paste0("Prevalencia COVID-19 de acuerdo a la región ", range_pos_18_0723$fecha)
#| label: tbl-prev-reg

tbl_summ_prev_reg_pre <- results_prevalencia$region$index %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  ) 

tbl_summ_prev_reg <- tbl_summ_prev_reg_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

tbl_summ_prev_reg
```

```{r}
tbl_summ_prev_reg_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:2) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/prevalencia_region.xlsx",
                       open = FALSE)

tbl_summ_prev_reg %>%
  gt::gtsave("02_output/tables/prevalencia_region.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_region.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_region.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Macrorregión

```{r}
#| fig-cap: !expr paste0("Prevalencia COVID-19 de acuerdo a la macrorregión ", range_pos_18_0723$fecha)
#| label: fig-prev-macro-bar

results_prevalencia$macrorregion$index %>% 
  ggplot(aes(x = macrorregion, y = rate, fill = macrorregion)) +
  geom_bar(stat = "identity") +
  theme_bw() +
  labs(title = "Monthly incidence of COVID-19 by macroregions", 
       x = NULL, 
       y = "Mean monthly COVID-19 incidence",
       fill = "Sex",
       caption = paste0("Total cases reported: ", 
                        scales::number(sum(pull(results_prevalencia$macrorregion$index_pre,
                                                count))))) +
  innovar::scale_fill_innova("blue_fall") +
  theme(
    legend.position = "bottom"
  )
```

```{r}
prev_macro_1k_sf <- peru_macro_sf %>% 
  inner_join(
    results_prevalencia$macrorregion$index
  )
```

```{r}
#| fig-cap: !expr paste0("Prevalencia COVID-19 de acuerdo a la macrorregión ", range_pos_18_0723$fecha)
#| label: fig-prev-macro-map
#| fig-width: 5
#| fig-height: 5

prev_macro_map <- prev_macro_1k_sf %>%
  mutate(
    macrorregion = case_match(
      macrorregion,
      "Metropolitan Lima and Callao" ~ "Met Lima & Callao",
      .default = macrorregion
    ),
    macrorregion = str_wrap(macrorregion, width = 15)
  ) %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  ggsflabel::geom_sf_label_repel(aes(label = macrorregion),
                                 force = 100,
                                 seed = 2024) +
  guides(
    fill = guide_legend(
      str_wrap("Mean monthly COVID-19 incidence", 
               20)
    )
  ) +
  scale_fill_distiller(palette = "Oranges", # Old: YlOrBr
                       direction = 1) + 
  theme_void()

prev_macro_map
```

```{r}
#| tbl-cap: !expr paste0("Prevalencia COVID-19 de acuerdo a la macrorregión ", range_pos_18_0723$fecha)
#| label: tbl-prev-macro

tbl_summ_prev_macro_pre <- results_prevalencia$macrorregion$index %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  ) 

tbl_summ_prev_macro <- tbl_summ_prev_macro_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

tbl_summ_prev_macro
```

```{r}
tbl_summ_prev_macro_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:2) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/prevalencia_macrorregion.xlsx",
                       open = FALSE)

tbl_summ_prev_macro %>%
  gt::gtsave("02_output/tables/prevalencia_macrorregion.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_macrorregion.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_macrorregion.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Distrito y sexo

<!-- Agregar un análisis solo de lima -->

```{r}
prev_distr_sex_1k_sf <- Peru %>%
  inner_join(
    results_prevalencia$distrito_sexo$index
  ) %>%
  ungroup()
```

```{r}
#| fig-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al distrito y sexo ", range_pos_18_0723$fecha)
#| label: fig-prev-distr-sex-map

prev_distr_sex_1k_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate), 
          linewidth = 0.01,
          color = "grey75") +
  guides(
    fill = guide_legend(
      str_wrap(
        "Mean monthly COVID-19 incidence",
        20
      )
    )
  ) +
  # scale_fill_distiller(palette = "YlOrBr",
  #                      direction = 1) + 
  scale_fill_gradientn(
    colors = scales::pal_brewer("seq", "Oranges", 1)(6),
    # colors = c("#ffffcc", "#fed976", "#feb24c", "#fc4e2a", "#e31a1c", "#b10026"),
    values = scales::rescale(c(0, 25, 50, 75, 
                               100, max(prev_distr_sex_1k_sf$rate))),
    breaks = c(0, 25, 50, 75, 
               100, max(prev_distr_sex_1k_sf$rate)),
    labels = c("0", "25", "50", "75", "100", "Max"),
    limits = c(0,
               max(prev_distr_sex_1k_sf$rate))
  ) +
  facet_wrap(vars(sexo)) + 
  theme_void()
```

```{r}
#| tbl-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al distrito y sexo ", range_pos_18_0723$fecha, " [Tabla Recortada]")
#| label: tbl-prev-distr-sex

set.seed(2024)

tbl_summ_prev_distr_sex_pre <- results_prevalencia$distrito_sexo$index %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  ) %>% 
  mutate(
    across(departamento:distrito,
           str_to_title)
  ) %>% 
  pivot_wider(
    names_from = sexo, 
    values_from = rate
  ) %>% 
  slice_sample(n = 2, by = c(departamento)) 

tbl_summ_prev_distr_sex <- tbl_summ_prev_distr_sex_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

tbl_summ_prev_distr_sex
```

```{r}
#| eval: false
tbl_summ_prev_distr_sex_pre <- results_prevalencia$distrito_sexo$index %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  ) %>% 
  mutate(
    across(departamento:distrito,
           str_to_title)
  ) %>% 
  pivot_wider(
    names_from = sexo, 
    values_from = rate
  )

tbl_summ_prev_distr_sex <- tbl_summ_prev_distr_sex_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )


tbl_summ_prev_distr_sex_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:6) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/prevalencia_distrito_sexo.xlsx",
                       open = FALSE)

tbl_summ_prev_distr_sex %>%
  gt::gtsave("02_output/tables/prevalencia_distrito_sexo.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_distrito_sexo.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_distrito_sexo.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Distrito y edad

```{r}
prev_distr_edad_mother <- peru_mother %>% 
  select(dep:ubigeo) %>% 
  collect() %>% 
  crossing(
    grupo_edad = levels(results_prevalencia$distrito_grupo_edad$index$grupo_edad)
  ) %>% 
  mutate(
    grupo_edad = fct_relevel(grupo_edad,
                             "Under 1 year",
                             "1 - 5 years",
                             "6 - 11 years")
  )

prev_distr_edad_100k <- prev_distr_edad_mother %>% 
  left_join(results_prevalencia$distrito_grupo_edad$index)

prev_distr_edad_1k_sf <- Peru %>%
  inner_join(
    prev_distr_edad_100k
  ) %>%
  ungroup()
```

```{r}
#| fig-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al distrito y edad ", range_pos_18_0723$fecha)
#| label: fig-prev-distr-edad-map

prev_distr_edad_1k_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate), 
          size = 0.01,
          color = "grey70") +
  guides(
    fill = guide_legend(
      str_wrap("Mean monthly COVID-19 incidence",
               20)
    )
  ) +
  # scale_fill_distiller(palette = "YlOrBr",
  #                      direction = 1) +
  # scale_fill_gradient(
  #   low = "#FEEDDE",
  #   high = "#A63603",
  #   na.value = "white"
  # ) +
  scale_fill_gradientn(
    colors = scales::pal_brewer("seq", "Oranges", 1)(6),
    # colors = c("#ffffcc", "#fed976", "#feb24c", "#fc4e2a", "#e31a1c", "#b10026"),
    values = scales::rescale(c(0, 25, 50, 75,
                               100, max(prev_distr_edad_1k_sf$rate, na.rm = TRUE))),
    breaks = c(0, 25, 50, 75,
               100, max(prev_distr_edad_1k_sf$rate, na.rm = TRUE)),
    labels = c("0", "25", "50", "75", "100", "Max"),
    limits = c(0,
               max(prev_distr_edad_1k_sf$rate, na.rm = TRUE)),
    na.value = "grey70"
  ) +
  facet_wrap(vars(grupo_edad), nrow = 1) + 
  theme_void()
```

```{r}
#| tbl-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al distrito y grupo de edad ", range_pos_18_0723$fecha, " [Tabla Recortada]")
#| label: tbl-prev-distr-edad

set.seed(2024)

tbl_summ_prev_distr_edad_pre <- results_prevalencia$distrito_grupo_edad$index %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  ) %>% 
  mutate(
    across(departamento:distrito,
           str_to_title)
  ) %>% 
  pivot_wider(
    names_from = grupo_edad, 
    values_from = rate
  ) %>% 
  slice_sample(n = 2, by = c(departamento)) 

tbl_summ_prev_distr_edad <- tbl_summ_prev_distr_edad_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

tbl_summ_prev_distr_edad
```

```{r}
#| eval: false
tbl_summ_prev_distr_edad_pre <- results_prevalencia$distrito_grupo_edad$index %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  ) %>% 
  mutate(
    across(departamento:distrito,
           str_to_title)
  ) %>% 
  pivot_wider(
    names_from = grupo_edad, 
    values_from = rate
  ) 

tbl_summ_prev_distr_edad <- tbl_summ_prev_distr_edad_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )


tbl_summ_prev_distr_edad_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:6) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/prevalencia_distrito_edad.xlsx",
                       open = FALSE)

tbl_summ_prev_distr_edad %>%
  gt::gtsave("02_output/tables/prevalencia_distrito_edad.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_distrito_edad.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_distrito_edad.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Distrital y pobreza

<!-- Recalcular las prevalencia/mortalidad/vacunacion/hospitalizacion por quintil, por olas y por grupo de edad (por separado, relacio par a par) -->

```{r}
library(biscale)
```

```{r}
prev_distr_pobr_1k_sf <- Peru %>%
  inner_join(
    results_prevalencia$distrito_cuartil_pobreza$index
  ) %>%
  ungroup()

prev_distr_pobr_biscale <- bi_class(prev_distr_pobr_1k_sf,
  x = rate,
  y = cuartil_pobreza,
  style = "quantile", dim = 4
)
```

```{r}
#| fig-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al distrito y cuartil de pobreza", range_pos_18_0723$fecha)
#| label: fig-prev-distr-pobr-map

custom_pal <- c(
  "1-1" = "#dae7f3", # Azul muy claro - Bajo en COVID, bajo en pobreza
  "2-1" = "#b6cee1", # Azul medio claro - Bajo en COVID, moderado en pobreza
  "3-1" = "#92b5cf", # Azul medio - Bajo en COVID, alto en pobreza
  "4-1" = "#6e9cbd", # Azul oscuro - Bajo en COVID, muy alto en pobreza
  
  
  "1-2" = "#d8f2bc", # Verde muy claro - Alto en COVID, bajo en pobreza
  "2-2" = "#a9d985", # Verde medio claro - Alto en COVID, moderado en pobreza
  "3-2" = "#7bc04e", # Verde medio - Alto en COVID, alto en pobreza
  "4-2" = "#4ca717", # Verde oscuro - Alto en COVID, muy alto en pobreza
  
  
  "1-3" = "#f8d4a0", # Naranja muy claro - Bajo en pobreza, bajo en COVID
  "2-3" = "#f6b86c", # Naranja medio claro - Moderado en pobreza, bajo en COVID
  "3-3" = "#f59b38", # Naranja medio - Alto en pobreza, bajo en COVID
  "4-3" = "#f37e04", # Naranja oscuro - Muy alto en pobreza, bajo en COVID
  

  
  "1-4" = "#f8bfcb", # Rosa muy claro - Muy alto en COVID, bajo en pobreza
  "2-4" = "#f196a0", # Rosa medio claro - Muy alto en COVID, moderado en pobreza
  "3-4" = "#ea6d75", # Rosa medio - Muy alto en COVID, alto en pobreza
  "4-4" = "#e3444a"  # Rosa oscuro - Muy alto en COVID, muy alto en pobreza
)

legend <- bi_legend(
  pal = custom_pal,
  dim = 4,
  xlab = "Incidence",
  ylab = "Poverty Quartiles",
  size = 5
)

ggmap <- prev_distr_pobr_biscale %>%
  ggplot() +
  geom_sf(aes(fill = bi_class), 
          linewidth = 0.01,
          show.legend = FALSE) +
  bi_scale_fill(pal = custom_pal, dim = 4) +
  # bi_scale_color(pal = "GrPink", dim = 3) +
  bi_theme()

library(cowplot)

finalPlot <- ggdraw(ggmap) +
  # draw_plot(ggmap, 0, 0, 1.5, 1.5) +
  draw_plot(legend, 0.1, 0.25, 0.27, 0.27)
finalPlot
```

```{r}
#| label: tbl-summ-positivos-cuartiles
#| tbl-cap: !expr paste0("Cuartiles de pobreza y casos positivos menores de 18 años ", range_pos_18_0723$fecha)

positivos_quartiles <- results_prevalencia$distrito_cuartil_pobreza$index %>% 
  select(prevalencia_promedio = rate, departamento, cuartil_pobreza) %>% 
  gtsummary::tbl_summary(
    by = cuartil_pobreza
  ) %>% 
  gtsummary::add_overall(last = TRUE) %>% 
  gtsummary::modify_header(label = "**Variables**") %>%
  gtsummary::bold_labels()
  

positivos_quartiles
```

```{r}
positivos_quartiles %>% 
  gtsummary::as_hux_xlsx("02_output/tables/prevalencia_distritos_pobreza_descriptivos_by_dep.xlsx")

positivos_quartiles %>%
  gtsummary::as_flex_table() %>%
  flextable::save_as_docx(path = "02_output/tables/prevalencia_distritos_pobreza_descriptivos_by_dep.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_distritos_pobreza_descriptivos_by_dep.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_distritos_pobreza_descriptivos_by_dep.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

```{r}
#| fig-cap: !expr paste0("Prevalencia COVID-19 de acuerdo al cuartil de pobreza", range_pos_18_0723$fecha)
#| label: fig-prev-pobr-area
#| fig-width: 10

prev_pobr_area_plot <- 
  results_per_month_prevalencia$cuartil_pobreza$index %>% 
  drop_na() %>% 
  mutate(
    time = make_date(anio, mes, "1")
  ) %>% 
  ggplot(
    aes(x = time,
        y = rate,
        fill = cuartil_pobreza)
  ) +
  geom_area(position = 'stack') +
  # innovar::scale_fill_innova("jama") +
  scale_fill_manual(
    values = c("#2C9BB4", "#0B4E60", "#DAE2BC", "#8C3A33")
  ) +
  scale_x_date(
    limits = c(ymd("2020-03-01"), ymd("2023-07-01")),
    breaks = seq(from = ceiling_date(ymd("2020-03-01"), "month"),
                 to = floor_date(ymd("2023-07-01"), "month"),
                 by = "2 months"), 
    # date_breaks = "2 month",
    date_labels = "%b %Y"
  ) +
  labs(
    y = NULL,
    x = NULL,
    title = "Monthly incidence of COVID-19",
    fill = "Poverty Quartiles"
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2020-03-01"), 
    xmax = as.Date("2020-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2020-07-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "First wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2020-11-01"), 
    xmax = as.Date("2021-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2021-05-10"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Second wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2021-11-01"), 
    xmax = as.Date("2022-04-30"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-02-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Third \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2022-05-01"), 
    xmax = as.Date("2022-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-08-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Fourth \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2022-11-01"), 
    xmax = as.Date("2023-07-01"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2023-03-01"),
    y = Inf, # Posición vertical del texto
    label = "Fifth wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  theme_classic() +
  theme(
    axis.text.x = element_text(
      angle = 75,
      vjust = 0.5
    ),
    axis.ticks.x = element_line()
  ) 

prev_pobr_area_plot
```

## Distrito Lima

```{r}
#| fig-cap: !expr paste0("Prevalencia COVID-19 de acuerdo a los distritos de lima metropolitana ", range_pos_18_0723$fecha)
#| label: fig-prev-lima-bar
#| fig-height: 9

total_cases_distr_lim_callao <- results_prevalencia$distrito$index_pre %>% 
  filter(provincia %in% c("LIMA", "CALLAO")) %>% 
  pull(count) %>% 
  sum()

results_prevalencia$distrito$index %>% 
  filter(provincia %in% c("LIMA", "CALLAO")) %>% 
  mutate(
    distrito = str_to_title(distrito), 
    distrito = fct_reorder(distrito, rate)
  ) %>% 
  ggplot(aes(y = distrito, x = rate)) +
  geom_bar(stat = "identity",
           color = "black",
           fill = "grey80") +
  theme_bw() +
  labs(title = "Monthly incidence of COVID-19 by districts of metropolitan lima", 
       x = NULL, 
       y = "Mean monthly COVID-19 incidence",
       caption = paste0("Total cases reported: ", 
                        scales::number(total_cases_distr_lim_callao))) +
  innovar::scale_fill_innova("blue_fall") +
  theme(
    legend.position = "bottom"
  )
```

```{r}
prev_lima_distr_1k_sf <- Peru %>% 
  inner_join(
    results_prevalencia$distrito$index %>% 
      filter(provincia %in% c("LIMA", "CALLAO"))
  )
```

```{r}
#| fig-cap: !expr paste0("Prevalencia COVID-19 de acuerdo a los distritos de lima metropolitana ", range_pos_18_0723$fecha)
#| label: fig-prev-lima-map

prev_lima_distr_1k_sf %>%
  mutate(
    label = case_when(
      rate > 5 ~ str_to_title(distrito)
    )
  ) %>% 
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap("Mean monthly COVID-19 incidence", 
               20)
    )
  ) +
  ggsflabel::geom_sf_label_repel(aes(label = label),
                                 force = 10,
                                 size = 2.45,
                                 box.padding = 0.75,
                                 seed = 2024) +
  scale_fill_distiller(palette = "Oranges", # Old: YlOrBr
                       direction = 1) + 
  theme_void()
```

```{r}
#| tbl-cap: !expr paste0("Prevalencia COVID-19 de acuerdo a los distritos de lima metropolitana ", range_pos_18_0723$fecha)
#| label: tbl-prev-lima

tbl_summ_prev_lima_distr_pre <- results_prevalencia$distrito$index %>% 
      filter(provincia %in% c("LIMA", "CALLAO")) %>% 
  mutate(
    across(
      departamento:distrito,
      str_to_title
    ),
    rate = formattable::digits(rate, 3)
  ) %>% 
  select(-departamento)

tbl_summ_prev_lima_distr <- tbl_summ_prev_lima_distr_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

tbl_summ_prev_lima_distr
```

```{r}
tbl_summ_prev_lima_distr_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:4) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/prevalencia_lima_distrito.xlsx",
                       open = FALSE)

tbl_summ_prev_lima_distr %>%
  gt::gtsave("02_output/tables/prevalencia_lima_distrito.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_lima_distrito.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/prevalencia_lima_distrito.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

# Fallecidos - Mortalidad

<!-- En letalidad y hospitalización por 100 casos -->

```{r}
fallecidos_total <- fallecidos_total %>% 
  filter(
    edad_declarada < 18,
    fecha_fallecimiento <= max_positivos_0723
  ) %>% 
  mutate(
    dep.code = str_sub(ubigeo, 1, 2),
    overall = "overall"
  ) 

range_fallecidos <- fallecidos_total %>% 
  summarise(
    min = min(fecha_fallecimiento),
    max = max(fecha_fallecimiento)
  ) %>% 
  collect() %>% 
  mutate(
    min = format(min, "%d-%m-%Y"),
    max = format(max, "%d-%m-%Y"),
    fecha = paste0("(", min, " a ", max, ")")
  )
```

## Vista general

**Cálculo de mortalidad**

```{r}
levels <- list(
  c("overall"),
  c("grupo_edad"),
  c("sexo"),
  c("grupo_edad", "sexo"),
  c("ola_covid"),
  c("ola_covid", "grupo_edad"),
  c("ola_covid", "sexo"),
  c("ola_covid", "region"),
  c("ola_covid", "macrorregion"),
  c("departamento"),
  c("departamento", "sexo"),
  c("departamento", "grupo_edad"),
  c("departamento", "sexo", "grupo_edad"),
  c("departamento", "ola_covid"),
  c("distrito"),
  c("distrito", "sexo"),
  c("distrito", "grupo_edad"),
  c("distrito", "cuartil_pobreza"),
  c("grupo_edad", "cuartil_pobreza"),
  c("sexo", "cuartil_pobreza"),
  c("ola_covid", "cuartil_pobreza"),
  c("ola_covid", "cuartil_nbi"),
  c("cuartil_pobreza"),
  c("cuartil_nbi"),
  c("region"),
  c("macrorregion")
)

# Calcular la mortalidad
results_mortalidad <- index_calculator(
  data = fallecidos_total, 
  denom = poblacion %>% mutate(overall = "overall"), 
  levels = levels, 
  type = "month",
  std = 100000
)

results_per_month_mortalidad <- index_calculator(
  data = fallecidos_total, 
  denom = poblacion %>% mutate(overall = "overall"), 
  levels = levels, 
  type = "per_month",
  std = 100000
)
```

**Cálculo de letalidad**

```{r}
# Calcular la letalidad
results_letalidad_mensual <- index_calculator(
  data = fallecidos_total, 
  denom = positivos_0723 %>% 
    mutate(cantidad = 1), 
  levels = levels, 
  type = "month",
  std = 100
)

results_per_month_letalidad <- index_calculator(
  data = fallecidos_total, 
  denom = positivos_0723 %>% 
    mutate(cantidad = 1), 
  levels = levels, 
  type = "per_month",
  std = 100
)

results_letalidad_overall <- index_calculator(
  data = fallecidos_total, 
  denom = positivos_0723 %>% 
    mutate(cantidad = 1), 
  levels = levels, 
  type = "overall",
  std = 100
)
```

```{r}
fallecidos_total_df <- fallecidos_total %>% 
  select(overall, edad_declarada, grupo_edad, sexo, ola_covid, region, 
         macrorregion, cuartil_pobreza, cuartil_nbi) %>% 
  collect()

labelled::var_label(fallecidos_total_df) <- list(
  overall = "Overall",
  grupo_edad = "Age Group",
  sexo = "Sex",
  ola_covid = "COVID-19 wave",
  region = "Natural Regions",
  macrorregion = "Macroregions",
  cuartil_pobreza = "Poverty level",
  cuartil_nbi = "Unsatisfied basic needs level"
)
```

```{r}
#| label: tbl-summ-fallecidos
#| tbl-cap: !expr paste0("Descriptivos de fallecidos menores de 18 años ", range_fallecidos$fecha)

tbl_summ_fallecidos_g <- fallecidos_total_df %>% 
  select(-edad_declarada) %>% 
  gtsummary::tbl_summary(
    type = list(
      overall ~ "dichotomous"
    ),
    value = list(overall = "overall"),
    missing_text = "No data"
  ) %>% 
  gtsummary::modify_header(label = "**Variables**") %>%
  gtsummary::bold_labels()

tbl_summ_fallecidos_c2 <- results_mortalidad$overall$index %>% 
  mutate(label = "Overall") %>% 
  rename(variable = overall) %>% 
  bind_rows(
    results_mortalidad$grupo_edad$index %>% 
      mutate(variable = "grupo_edad") %>% 
      rename(label = grupo_edad),
    results_mortalidad$sexo$index %>% 
      mutate(variable = "sexo") %>% 
      rename(label = sexo),
    results_mortalidad$ola_covid$index %>% 
      mutate(variable = "ola_covid") %>% 
      rename(label = ola_covid),
    results_mortalidad$region$index %>% 
      mutate(variable = "region") %>% 
      rename(label = region),
    results_mortalidad$macrorregion$index %>% 
      mutate(variable = "macrorregion") %>% 
      rename(label = macrorregion),
    results_mortalidad$cuartil_pobreza$index %>% 
      mutate(variable = "cuartil_pobreza") %>% 
      rename(label = cuartil_pobreza),
    results_mortalidad$cuartil_nbi$index %>% 
      mutate(variable = "cuartil_nbi") %>% 
      rename(label = cuartil_nbi)
  ) %>% 
  drop_na() %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  )

tbl_summ_fallecidos_c3 <- results_letalidad_mensual$overall$index %>% 
  mutate(label = "Overall") %>% 
  rename(variable = overall) %>% 
  bind_rows(
    results_letalidad_mensual$grupo_edad$index %>% 
      mutate(variable = "grupo_edad") %>% 
      rename(label = grupo_edad),
    results_letalidad_mensual$sexo$index %>% 
      mutate(variable = "sexo") %>% 
      rename(label = sexo),
    results_letalidad_mensual$ola_covid$index %>% 
      mutate(variable = "ola_covid") %>% 
      rename(label = ola_covid),
    results_letalidad_mensual$region$index %>% 
      mutate(variable = "region") %>% 
      rename(label = region),
    results_letalidad_mensual$macrorregion$index %>% 
      mutate(variable = "macrorregion") %>% 
      rename(label = macrorregion),
    results_letalidad_mensual$cuartil_pobreza$index %>% 
      mutate(variable = "cuartil_pobreza") %>% 
      rename(label = cuartil_pobreza),
    results_letalidad_mensual$cuartil_nbi$index %>% 
      mutate(variable = "cuartil_nbi") %>% 
      rename(label = cuartil_nbi)
  ) %>% 
  drop_na() %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  ) %>% 
  rename(letalidad_mensual = rate)

tbl_summ_fallecidos_c4 <- results_letalidad_overall$overall$index %>% 
  mutate(label = "Overall") %>% 
  rename(variable = overall) %>% 
  bind_rows(
    results_letalidad_overall$grupo_edad$index %>% 
      mutate(variable = "grupo_edad") %>% 
      rename(label = grupo_edad),
    results_letalidad_overall$sexo$index %>% 
      mutate(variable = "sexo") %>% 
      rename(label = sexo),
    results_letalidad_overall$ola_covid$index %>% 
      mutate(variable = "ola_covid") %>% 
      rename(label = ola_covid),
    results_letalidad_overall$region$index %>% 
      mutate(variable = "region") %>% 
      rename(label = region),
    results_letalidad_overall$macrorregion$index %>% 
      mutate(variable = "macrorregion") %>% 
      rename(label = macrorregion),
    results_letalidad_overall$cuartil_pobreza$index %>% 
      mutate(variable = "cuartil_pobreza") %>% 
      rename(label = cuartil_pobreza),
    results_letalidad_overall$cuartil_nbi$index %>% 
      mutate(variable = "cuartil_nbi") %>% 
      rename(label = cuartil_nbi)
  ) %>% 
  drop_na() %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  ) %>% 
  rename(letalidad_overall = rate)


tbl_summ_fallecidos_final <- tbl_summ_fallecidos_g %>%
  gtsummary::modify_table_body(
    ~ .x %>% 
      left_join(tbl_summ_fallecidos_c2) %>% 
      left_join(tbl_summ_fallecidos_c3) %>% 
      left_join(tbl_summ_fallecidos_c4) %>% 
      mutate(
        label = case_when(
          variable %in% c("cuartil_pobreza", "cuartil_nbi") &
            label == "1" ~ "1 (0-19.9%)",
          variable %in% c("cuartil_pobreza", "cuartil_nbi") & 
            label == "2" ~ "2 (20-39.9%)",
          variable %in% c("cuartil_pobreza", "cuartil_nbi") & 
            label == "3" ~ "3 (40-59.9%)",
          variable %in% c("cuartil_pobreza", "cuartil_nbi") &
            label == "4" ~ "4 (60-100%)",
          .default = label
        )
      )
  ) %>% 
  gtsummary::modify_header(
    rate ~ "**Average monthly mortality per 100K habitants**",
    letalidad_mensual ~ "**Average monthly lethality per 100 cases**",
    letalidad_overall ~ "**Lethality per 100 cases**"
  )

tbl_summ_fallecidos_final
```

```{r}
tbl_summ_fallecidos_final %>% 
  gtsummary::as_hux_xlsx("02_output/tables/fallecidos_descriptivos.xlsx")

tbl_summ_fallecidos_final %>%
  gtsummary::as_flex_table() %>%
  flextable::save_as_docx(path = "02_output/tables/fallecidos_descriptivos.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/fallecidos_descriptivos.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/fallecidos_descriptivos.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

<!-- Cálculo de Comparaciones Múltiples [Mortalidad] -->

```{r}
mortality_comparisons <- compare_levels_per_month(
  results_per_month = results_per_month_mortalidad,
  levels = comparing_levels
)

# Extraer las tablas `cld` como tibbles, renombrando y añadiendo una columna "level"
mortality_comparisons_cld_pre <- lapply(names(mortality_comparisons), function(level) {
  cld <- mortality_comparisons[[level]]$cld
  
  if (!is.null(cld)) {
    cld <- cld |> 
      as_tibble() |> 
      rename(group = 1) |>  # Renombra la primera columna a 'group'
      mutate(level = level) |> # Añade una columna 'level' con el nombre de la lista
      relocate(level) 
  }
  
  cld
}) |> 
  bind_rows()

mortality_comparisons_cld_pre <- 
  mortality_comparisons_cld_pre |> 
  mutate(
    across(
      where(is.numeric),
      ~ formattable::digits(.x, 3)
    )
  ) 
```

```{r}
#| label: tbl-compar-mortalidad
#| tbl-cap: !expr paste0("Comparaciones múltiples de mortalidad", range_pos_18_0723$fecha)
#| 
mortality_comparisons_cld <- mortality_comparisons_cld_pre |> 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

mortality_comparisons_cld
```

```{r}
mortality_comparisons_cld_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:8) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/mortalidad_comparacion_multiple.xlsx",
                       open = FALSE)

mortality_comparisons_cld %>%
  gt::gtsave("02_output/tables/mortalidad_comparacion_multiple.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/mortalidad_comparacion_multiple.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/mortalidad_comparacion_multiple.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

<!-- Cálculo de Comparaciones Múltiples [Letalidad] -->

```{r}
lethality_comparisons <- compare_levels_per_month(
  results_per_month = results_per_month_letalidad,
  levels = comparing_levels
)

# Extraer las tablas `cld` como tibbles, renombrando y añadiendo una columna "level"
lethality_comparisons_cld_pre <- lapply(names(lethality_comparisons), function(level) {
  cld <- lethality_comparisons[[level]]$cld
  
  if (!is.null(cld)) {
    cld <- cld |> 
      as_tibble() |> 
      rename(group = 1) |>  # Renombra la primera columna a 'group'
      mutate(level = level) |> # Añade una columna 'level' con el nombre de la lista
      relocate(level) 
  }
  
  cld
}) |> 
  bind_rows()

lethality_comparisons_cld_pre <- 
  lethality_comparisons_cld_pre |> 
  mutate(
    across(
      where(is.numeric),
      ~ formattable::digits(.x, 3)
    )
  ) 
```

```{r}
#| label: tbl-compar-letalidad
#| tbl-cap: !expr paste0("Comparaciones múltiples de letalidad", range_pos_18_0723$fecha)
#| 
lethality_comparisons_cld <- lethality_comparisons_cld_pre |> 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

lethality_comparisons_cld
```

```{r}
lethality_comparisons_cld_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:8) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/letalidad_comparacion_multiple.xlsx",
                       open = FALSE)

lethality_comparisons_cld %>%
  gt::gtsave("02_output/tables/letalidad_comparacion_multiple.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/letalidad_comparacion_multiple.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/letalidad_comparacion_multiple.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

<!-- Comparaciones de fallecidos con respecto a olas -->

```{r}
#| label: tbl-summ-fallecidos-olas
#| tbl-cap: !expr paste0("Descriptivos de fallecidos menores de 18 años por olas, sexo y grupo de edad ", range_fallecidos$fecha)

tbl_summ_fallecidos_olas <- fallecidos_total_df %>% 
  select(-edad_declarada) %>% 
  gtsummary::tbl_summary(
    by = ola_covid,
    type = list(
      overall ~ "dichotomous"
    ),
    value = list(overall = "overall"),
    missing_text = "No data"
  ) %>% 
  gtsummary::modify_header(label = "**Variables**") %>%
  gtsummary::bold_labels()



ola_covid_fallecidos_gtsumm <- results_mortalidad$ola_covid$index %>% 
  rename(stat = ola_covid) %>% 
  mutate(label = "Overall",
         variable = "overall") %>% 
  bind_rows(
    results_mortalidad$ola_covid_grupo_edad$index %>%
      mutate(variable = "grupo_edad") %>% 
      rename(label = grupo_edad,
             stat = ola_covid),
    results_mortalidad$ola_covid_sexo$index %>% 
      mutate(variable = "sexo") %>% 
      rename(label = sexo,
             stat = ola_covid),
    results_mortalidad$ola_covid_region$index %>% 
      mutate(variable = "region") %>% 
      rename(label = region,
             stat = ola_covid),
    results_mortalidad$ola_covid_macrorregion$index %>% 
      mutate(variable = "macrorregion") %>%
      rename(label = macrorregion,
             stat = ola_covid),
    results_mortalidad$ola_covid_cuartil_pobreza$index %>% 
      mutate(variable = "cuartil_pobreza") %>% 
      rename(label = cuartil_pobreza,
             stat = ola_covid),
    results_mortalidad$ola_covid_cuartil_nbi$index %>% 
      mutate(variable = "cuartil_nbi") %>% 
      rename(label = cuartil_nbi,
             stat = ola_covid)
  ) %>% 
  mutate(
    rate = formattable::digits(rate, 3),
    stat = as.character(stat),
    stat = case_match(
      stat,
      "1 (Mar – Oct 2020)" ~ "stat_1",
      "2 (Nov 2020 – Oct 2021)" ~ "stat_2",
      "3 (Nov 2021 – Apr 2022)" ~ "stat_3",
      "4 (May 2022 – Oct 2022)" ~ "stat_4",
      "5 (Nov 2022 – Jun 2023)" ~ "stat_5"
    )
  ) 

tbl_summ_fallecidos_olas$table_body <- tbl_summ_fallecidos_olas$table_body %>%
  mutate(
    across(
      gtsummary::all_stat_cols(stat_0 = FALSE),
      ~ str_remove_all(., "\\s*\\(.*?\\)\\s*")
    )
  ) %>% 
  pivot_longer(
    cols = stat_1:stat_5,
    names_to = "stat",
    values_to = "freq"
  ) %>% 
  left_join(
    ola_covid_fallecidos_gtsumm
  ) %>% 
  mutate(
    value = str_c(freq, ifelse(is.na(rate), "", str_c(" (", rate, ")")))
  ) %>% 
  select(-c(freq, rate)) %>% 
  pivot_wider(
    names_from = stat,
    values_from = value
  )

tbl_summ_fallecidos_olas <- tbl_summ_fallecidos_olas %>% 
  gtsummary::modify_footnote(
    gtsummary::all_stat_cols() ~ "n (mortality)"
  )

tbl_summ_fallecidos_olas
```

```{r}
tbl_summ_fallecidos_olas %>% 
  gtsummary::as_hux_xlsx("02_output/tables/fallecidos_olas_descriptivos.xlsx")

tbl_summ_fallecidos_olas %>%
  gtsummary::as_flex_table() %>%
  flextable::save_as_docx(path = "02_output/tables/fallecidos_olas_descriptivos.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/fallecidos_olas_descriptivos.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/fallecidos_olas_descriptivos.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Gráfico de muertes (serie de tiempo) - Menores de 18 años

```{r}
#| dpi: 300
#| fig-width: 10

range_fallecidos_und_18 <- fallecidos_total %>%
  summarise(
    min = min(fecha_fallecimiento),
    max = max(fecha_fallecimiento)
  ) %>% 
  collect()

deaths_cases_area <- fallecidos_total %>% 
  count(clasificacion_def, fecha_round) %>% 
  mutate(
    clasificacion_def = case_match(
      clasificacion_def,
      "Criterio serológico" ~ "Serological Criteria",
      "Criterio radiológico" ~ "Radiological Criteria",
      "Criterio SINADEF" ~ "SINADEF Criteria",
      "Criterio virológico" ~ "Virological Criteria",
      "Criterio nexo epidemiológico" ~ "Epidemiological Link Criteria",
      "Criterio clínico" ~ "Clinical Criteria",
      "Criterio investigación Epidemiológica" ~ "Epidemiological Investigation Criteria",
      .default = NA_character_
    )
  ) %>% 
  collect() %>% 
  mutate(
    clasificacion_def = str_wrap(clasificacion_def,
                                 20)
  ) |> 
  ggplot(
    aes(
      x = fecha_round,
      y = n,
      fill = clasificacion_def
    )
  ) +
  geom_area(position = 'stack') +
  innovar::scale_fill_innova("jama") +
  # scale_x_date(
  #   limits = c(ymd("2020-03-01"), ymd("2023-10-31")),
  #   breaks = scales::date_breaks("2 month"), # Cambiar a "6 months" para menos frecuencia
  #   labels = scales::date_format("%b %Y") # Formato de fecha, %b es la abreviatura de mes, %Y es el año
  # ) + 
  scale_x_date(
    limits = c(ymd("2020-03-01"), ymd("2023-07-01")),
    breaks = seq(from = ceiling_date(ymd("2020-03-01"), "month"),
                 to = floor_date(ymd("2023-07-01"), "month"),
                 by = "2 months"), 
    # date_breaks = "2 month",
    date_labels = "%b %Y"
  ) +
  labs(
    y = NULL,
    x = NULL,
    title = "COVID-19 deaths",
    fill = "Classification Criteria"
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2020-03-01"), 
    xmax = as.Date("2020-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2020-07-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "First wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2020-11-01"), 
    xmax = as.Date("2021-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2021-05-10"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Second wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2021-11-01"), 
    xmax = as.Date("2022-04-30"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
    annotate(
    geom = "text", 
    x = as.Date("2022-02-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Third \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2022-05-01"), 
    xmax = as.Date("2022-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-08-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Fourth \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2022-11-01"), 
    xmax = range_fallecidos_und_18$max, 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2023-03-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Fifth wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  theme_classic() +
  theme(
    axis.text.x = element_text(
      angle = 75,
      vjust = 0.5
    ),
    axis.ticks.x = element_line()
  ) 

deaths_cases_area
```

## Gráfico de muertes (serie de tiempo) - Grupo Etáreo

```{r}
#| dpi: 300
#| fig-width: 10

deaths_cases_by_age_area <- fallecidos_total %>%
  count(grupo_edad, fecha_round) %>% 
  collect() %>% 
  ggplot(
    aes(
      x = fecha_round,
      y = n,
      fill = grupo_edad
    )
  ) +
  geom_area(position = 'stack') +
  scale_fill_manual(
    values = c("#2C9BB4", "#0B4E60", "#DAE2BC", "#8C3A33")
  ) +
  # scale_x_date(
  #   limits = c(ymd("2020-03-01"), ymd("2023-10-31")),
  #   breaks = scales::date_breaks("2 month"), # Cambiar a "6 months" para menos frecuencia
  #   labels = scales::date_format("%b %Y") # Formato de fecha, %b es la abreviatura de mes, %Y es el año
  # ) + 
  scale_x_date(
    limits = c(ymd("2020-03-01"), ymd("2023-07-01")),
    breaks = seq(from = ceiling_date(ymd("2020-03-01"), "month"),
                 to = floor_date(ymd("2023-07-01"), "month"),
                 by = "2 months"), 
    # date_breaks = "2 month",
    date_labels = "%b %Y"
  ) +
  labs(
    y = NULL,
    x = NULL,
    title = "COVID-19 deaths",
    fill = "Age Group"
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2020-03-01"), 
    xmax = as.Date("2020-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2020-07-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "First wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2020-11-01"), 
    xmax = as.Date("2021-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2021-05-10"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Second wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2021-11-01"), 
    xmax = as.Date("2022-04-30"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
    annotate(
    geom = "text", 
    x = as.Date("2022-02-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Third \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2022-05-01"), 
    xmax = as.Date("2022-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-08-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Fourth \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2022-11-01"), 
    xmax = as.Date("2023-07-01"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2023-03-01"),
    y = Inf, # Posición vertical del texto
    label = "Fifth wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  theme_classic() +
  theme(
    axis.text.x = element_text(
      angle = 75,
      vjust = 0.5
    ),
    axis.ticks.x = element_line()
  ) 

deaths_cases_by_age_area
```

## Gráfico Mortalidad (serie de tiempo)

```{r}
mortality_cases_by_age_area <- results_per_month_mortalidad$grupo_edad$index %>% 
  mutate(
    time = make_date(anio, mes, "1")
  ) %>% 
  ggplot(
    aes(
      x = time,
      y = rate,
      fill = grupo_edad
    )
  ) +
  geom_area(position = 'stack') +
  scale_fill_manual(
    values = c("#2C9BB4", "#0B4E60", "#DAE2BC", "#8C3A33")
  ) +
  scale_x_date(
    limits = c(ymd("2020-03-01"), ymd("2023-07-01")),
    breaks = seq(from = ceiling_date(ymd("2020-03-01"), "month"),
                 to = floor_date(ymd("2023-07-01"), "month"),
                 by = "2 months"), 
    # date_breaks = "2 month",
    date_labels = "%b %Y"
  ) +
  labs(
    y = NULL,
    x = NULL,
    title = "Monthly COVID-19 mortality per 100 000 inhabitants",
    fill = "Age Group"
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2020-03-01"), 
    xmax = as.Date("2020-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2020-07-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "First wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2020-11-01"), 
    xmax = as.Date("2021-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2021-05-10"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Second wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2021-11-01"), 
    xmax = as.Date("2022-04-30"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
    annotate(
    geom = "text", 
    x = as.Date("2022-02-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Third \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2022-05-01"), 
    xmax = as.Date("2022-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-08-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Fourth \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2022-11-01"), 
    xmax = as.Date("2023-07-01"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2023-03-01"),
    y = Inf, # Posición vertical del texto
    label = "Fifth wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  theme_classic() +
  theme(
    axis.text.x = element_text(
      angle = 75,
      vjust = 0.5
    ),
    axis.ticks.x = element_line()
  ) 

mortality_cases_by_age_area
```

## Gráfico Letalidad (serie de tiempo)

No se está mostrando en la gráfica la información de marzo de 2020, debido a la enorme diferencia entre su valor de letalidad con respecto al resto de letalidades.

```{r}
results_per_month_letalidad$grupo_edad$index |> 
  slice(1:8) |> 
  mutate(rate = round(rate, 3)) |> 
  gt::gt()
```

```{r}
lethality_cases_by_age_area <- results_per_month_letalidad$grupo_edad$index %>% 
  mutate(
    time = make_date(anio, mes, "1")
  ) %>% 
  filter(
    time > ymd("2020-03-01")
  ) |> 
  ggplot(
    aes(
      x = time,
      y = rate,
      fill = grupo_edad
    )
  ) +
  geom_area(position = 'stack') +
  scale_fill_manual(
    values = c("#2C9BB4", "#0B4E60", "#DAE2BC", "#8C3A33")
  ) +
  scale_x_date(
    limits = c(ymd("2020-03-01"), ymd("2023-07-01")),
    breaks = seq(from = ceiling_date(ymd("2020-03-01"), "month"),
                 to = floor_date(ymd("2023-07-01"), "month"),
                 by = "2 months"), 
    # date_breaks = "2 month",
    date_labels = "%b %Y"
  ) +
  labs(
    y = NULL,
    x = NULL,
    title = "Monthly COVID-19 lethality per 100 cases",
    fill = "Age Group"
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2020-03-01"), 
    xmax = as.Date("2020-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2020-07-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "First wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2020-11-01"), 
    xmax = as.Date("2021-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2021-05-10"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Second wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2021-11-01"), 
    xmax = as.Date("2022-04-30"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
    annotate(
    geom = "text", 
    x = as.Date("2022-02-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Third \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2022-05-01"), 
    xmax = as.Date("2022-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-08-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Fourth \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2022-11-01"), 
    xmax = as.Date("2023-07-01"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2023-03-01"),
    y = Inf, # Posición vertical del texto
    label = "Fifth wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  theme_classic() +
  theme(
    axis.text.x = element_text(
      angle = 75,
      vjust = 0.5
    ),
    axis.ticks.x = element_line()
  ) 

lethality_cases_by_age_area
```

## Pirámide frecuencia - Edades

```{r}
#| dpi: 300
#| fig-width: 8

fallecidos_total_df %>% 
  mutate(
    edad_declarada = factor(edad_declarada)
  ) %>% 
  age_pyramid(
    age_group = "edad_declarada",
    split_by = "sexo"
  ) +
  scale_fill_brewer(palette = "Dark2") +
  labs(
    x = "Age",
    fill = "Sex",
    y = "Frequency"
  ) +
  theme_classic()
```

```{r}
#| dpi: 300
#| fig-width: 6

fallecidos_total_df %>% 
  age_pyramid(
    age_group = "grupo_edad",
    split_by = "sexo"
  ) +
  scale_fill_brewer(palette = "Dark2") +
  labs(
    x = "Age Group",
    fill = "Sex",
    y = "Frequency"
  ) +
  theme_classic()
```

## Pirámide proporcional - Edades

```{r}
#| dpi: 300
#| fig-width: 8

fallecidos_total_df %>% 
  mutate(
    edad_declarada = factor(edad_declarada)
  ) %>% 
  age_pyramid(
    age_group = "edad_declarada",
    split_by = "sexo",
    proportional = TRUE
  ) +
   scale_fill_brewer(palette = "Dark2") +
  labs(
    x = "Age",
    fill = "Sex",
    y = "Percentage"
  ) +
  theme_classic()
```

```{r}
#| dpi: 300
#| fig-width: 6

fallecidos_total_df %>% 
  age_pyramid(
    age_group = "grupo_edad",
    split_by = "sexo",
    proportional = TRUE
  ) +
   scale_fill_brewer(palette = "Dark2") +
  labs(
    x = "Age Group",
    fill = "Sex",
    y = "Percentage"
  ) +
  theme_classic()
```

## Departamento

```{r}
fallecidos_dep_100k_sf <- peru_dep_sf %>%
  rename(departamento = dep) %>% 
  left_join(
    results_mortalidad$departamento$index
  ) %>%
  ungroup()
```

```{r}
#| fig-cap: !expr paste0("Mortalidad COVID-19 por 100K habitantes de acuerdo al departamento ", range_fallecidos$fecha)
#| label: fig-mort-dep-map
#| fig-width: 5
#| fig-height: 5

fallecidos_dep_100k_map <- fallecidos_dep_100k_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap("Mean monthly COVID-19 mortality per 100 000 inhabitants",
               20)
    )
  ) +
  scale_fill_distiller(palette = "Oranges", # Old: YlOrBr
                       direction = 1) + 
  theme_void()

fallecidos_dep_100k_map
```

```{r}
#| tbl-cap: !expr paste0("Mortalidad COVID-19 por 100K habitantes de acuerdo al departamento ", range_pos_18_0723$fecha)
#| label: tbl-mort-dep

tbl_summ_mort_dep_pre <- results_mortalidad$departamento$index %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  ) %>% 
  mutate(
    departamento = str_to_title(departamento)
  ) 

tbl_summ_mort_dep <- tbl_summ_mort_dep_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

tbl_summ_mort_dep
```

```{r}
tbl_summ_mort_dep_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:2) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/mortalidad_departamento.xlsx",
                       open = FALSE)

tbl_summ_mort_dep %>%
  gt::gtsave("02_output/tables/mortalidad_departamento.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/mortalidad_departamento.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/mortalidad_departamento.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Departamento y sexo

```{r}
mortalidad_dep_sex_100k_sf <- peru_dep_sf %>% 
  rename(departamento = dep) %>% 
  inner_join(
    results_mortalidad$departamento_sexo$index
  ) 
```

```{r}
#| fig-cap: !expr paste0("Mortalidad COVID-19 por 100K habitantes de acuerdo al departamento y sexo ", range_fallecidos$fecha)
#| label: fig-mort-dep-sex-map

mortalidad_dep_sex_100k_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap("Mean monthly COVID-19 mortality per 100 000 inhabitants",
               20)
    )
  ) +
  scale_fill_distiller(palette = "Oranges", # Old: "YlOrBr"
                       direction = 1) + 
  facet_wrap(vars(sexo)) + 
  theme_void()
```

```{r}
#| tbl-cap: !expr paste0("Mortalidad COVID-19 por 100K habitantes de acuerdo al departamento ", range_fallecidos$fecha)
#| label: tbl-mort-dep-sex

tbl_summ_mort_dep_sex_pre <- results_mortalidad$departamento_sexo$index %>% 
  mutate(
    rate = formattable::digits(rate, 3),
    departamento = str_to_title(departamento)
  ) %>% 
  pivot_wider(
    names_from = sexo, 
    values_from = rate
  )


tbl_summ_mort_dep_sex <- tbl_summ_mort_dep_sex_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

tbl_summ_mort_dep_sex
```

```{r}
tbl_summ_mort_dep_sex_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:3) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/mortalidad_departamento_sexo.xlsx",
                       open = FALSE)

tbl_summ_mort_dep_sex %>%
  gt::gtsave("02_output/tables/mortalidad_departamento_sexo.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/mortalidad_departamento_sexo.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/mortalidad_departamento_sexo.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Departamento y edad

```{r}
mortalidad_dep_edad_100k_sf <- peru_dep_sf %>% 
  rename(departamento = dep) %>% 
  inner_join(
    results_mortalidad$departamento_grupo_edad$index
  ) 
```

```{r}
#| fig-cap: !expr paste0("Mortalidad COVID-19 por 100K habitantes de acuerdo al departamento y grupo de edad ", range_fallecidos$fecha)
#| label: fig-mort-dep-edad-map

mortalidad_dep_edad_100k_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap("Mean monthly COVID-19 mortality per 100 000 inhabitants",
               20)
    )
  ) +
  scale_fill_distiller(palette = "YlOrBr",
                       direction = 1) + 
  facet_wrap(vars(grupo_edad), nrow = 1) + 
  theme_void()
```

```{r}
#| tbl-cap: !expr paste0("Mortalidad COVID-19 por 100K habitantes de acuerdo al departamento y grupo de edad", range_fallecidos$fecha)
#| label: tbl-mort-dep-edad

tbl_summ_mort_dep_edad_pre <- results_mortalidad$departamento_grupo_edad$index %>% 
  mutate(
    rate = formattable::digits(rate, 3),
    departamento = str_to_title(departamento)
  ) %>% 
  pivot_wider(
    names_from = grupo_edad, 
    values_from = rate
  )


tbl_summ_mort_dep_edad <- tbl_summ_mort_dep_edad_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

tbl_summ_mort_dep_edad
```

```{r}
tbl_summ_mort_dep_edad_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:5) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/mortalidad_departamento_edad.xlsx",
                       open = FALSE)

tbl_summ_mort_dep_edad %>%
  gt::gtsave("02_output/tables/mortalidad_departamento_edad.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/mortalidad_departamento_edad.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/mortalidad_departamento_edad.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Departamento y olas

```{r}
mortalidad_dep_olas_mother <- crossing(
  ola_covid = levels(results_mortalidad$departamento_ola_covid$index$ola_covid),
  departamento = unique(results_mortalidad$departamento_ola_covid$index$departamento)
)

mortalidad_dep_olas_100k <- mortalidad_dep_olas_mother %>% 
  left_join(results_mortalidad$departamento_ola_covid$index) %>% 
  bind_rows(
    results_mortalidad$departamento$index %>% 
      mutate(ola_covid = "Overall")
  ) %>% 
  mutate(
    ola_covid = fct_recode(
      ola_covid,
      "1st wave" = "1 (Mar – Oct 2020)",
      "2nd wave" = "2 (Nov 2020 – Oct 2021)",
      "3rd wave" = "3 (Nov 2021 – Apr 2022)",
      "4th wave" = "4 (May 2022 – Oct 2022)",
      "5th wave" = "5 (Nov 2022 – Jun 2023)"
    ),
    ola_covid = fct_relevel(ola_covid,
                            "Overall")
  ) %>% 
  arrange(ola_covid)

mortalidad_dep_olas_100k_sf <- peru_dep_sf %>%
  rename(departamento = dep) %>% 
  left_join(
    mortalidad_dep_olas_100k
  ) %>%
  ungroup()
```

```{r}
#| fig-cap: !expr paste0("Mortalidad COVID-19 por 100K habitantes de acuerdo al departamento y Olas COVID ", range_fallecidos$fecha)
#| label: fig-mort-dep-ola-map
#| fig-width: 13
#| fig-height: 3

mortalidad_dep_olas_100k_map <- 
  mortalidad_dep_olas_100k_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap("Mean monthly COVID-19 mortality per 100 000 inhabitants",
               20)
    )
  ) +
  # scale_fill_distiller(palette = "YlOrBr",
  #                      direction = 1) + 
  scale_fill_gradientn(
    colors = scales::pal_brewer("seq", "Oranges", 1)(4),
    values = scales::rescale(c(0, 1, 2, 3, max(mortalidad_dep_olas_100k_sf$rate, na.rm = TRUE))),
    breaks = c(0, 1, 2, 3, max(mortalidad_dep_olas_100k_sf$rate, na.rm = TRUE)),
    labels = c("0", "1", "2", "3", "Max"),
    limits = c(0,
               max(mortalidad_dep_olas_100k_sf$rate, na.rm = TRUE)),
    na.value = "white"
  ) +
  facet_wrap(vars(ola_covid), nrow = 1) + 
  theme_void()

mortalidad_dep_olas_100k
```

```{r}
#| tbl-cap: !expr paste0("Mortalidad COVID-19 por 100K habitantes de acuerdo al departamento y Olas COVID", range_fallecidos$fecha)
#| label: tbl-mort-dep-ola

tbl_summ_mort_dep_ola_pre <- mortalidad_dep_olas_100k %>% 
  mutate(
    rate = formattable::digits(rate, 3),
    departamento = str_to_title(departamento)
  ) %>% 
  pivot_wider(
    names_from = ola_covid, 
    values_from = rate
  )


tbl_summ_mort_dep_ola <- tbl_summ_mort_dep_ola_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

tbl_summ_mort_dep_ola
```

```{r}
tbl_summ_mort_dep_ola_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:6) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/mortalidad_departamento_ola.xlsx",
                       open = FALSE)

tbl_summ_mort_dep_ola %>%
  gt::gtsave("02_output/tables/mortalidad_departamento_ola.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/mortalidad_departamento_ola.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/mortalidad_departamento_ola.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Regiones

```{r}
mortalidad_reg_100k_sf <- peru_reg_sf %>%
  left_join(
    results_mortalidad$region$index
  ) 
```

```{r}
#| fig-cap: !expr paste0("Mortalidad COVID-19 por 100K habitantes de acuerdo a la región ", range_fallecidos$fecha)
#| label: fig-mort-reg-map
#| fig-width: 5
#| fig-height: 5

mortalidad_reg_100k_map <- mortalidad_reg_100k_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap("Mean monthly COVID-19 mortality per 100 000 inhabitants",
               20)
    )
  ) +
  ggsflabel::geom_sf_label(aes(label = region)) +
  scale_fill_distiller(palette = "Oranges", # Old: YlOrBr
                       direction = 1) + 
  theme_void()

mortalidad_reg_100k_map
```

```{r}
#| tbl-cap: !expr paste0("Mortalidad COVID-19 por 100K habitantes de acuerdo a la región", range_fallecidos$fecha)
#| label: tbl-mort-reg

tbl_summ_mort_reg_pre <- results_mortalidad$region$index %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  )


tbl_summ_mort_reg <- tbl_summ_mort_reg_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

tbl_summ_mort_reg
```

```{r}
tbl_summ_mort_reg_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:2) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/mortalidad_region.xlsx",
                       open = FALSE)

tbl_summ_mort_reg %>%
  gt::gtsave("02_output/tables/mortalidad_region.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/mortalidad_region.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/mortalidad_region.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Macrorregiones

```{r}
mortalidad_macro_100k_sf <- peru_macro_sf %>%
  left_join(
    results_mortalidad$macrorregion$index
  ) 
```

```{r}
#| fig-cap: !expr paste0("Mortalidad COVID-19 por 100K habitantes de acuerdo a la macrorregión ", range_fallecidos$fecha)
#| label: fig-mort-macro-map
#| fig-width: 5
#| fig-height: 5

mort_macro_map <- mortalidad_macro_100k_sf %>%
  mutate(
    macrorregion = str_wrap(macrorregion, width = 15)
  ) %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  ggsflabel::geom_sf_label_repel(aes(label = macrorregion),
                                 force = 100,
                                 seed = 2024) +
  guides(
    fill = guide_legend(
      str_wrap("Mean monthly COVID-19 mortality per 100 000 inhabitants", 
               20)
    )
  ) +
  scale_fill_distiller(palette = "Oranges", # Old: YlOrBr
                       direction = 1) + 
  theme_void()

mort_macro_map
```

```{r}
#| tbl-cap: !expr paste0("Mortalidad COVID-19 por 100K habitantes de acuerdo a la macrorregión", range_fallecidos$fecha)
#| label: tbl-mort-macro

tbl_summ_mort_macro_pre <- results_mortalidad$macrorregion$index %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  )


tbl_summ_mort_macro <- tbl_summ_mort_macro_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

tbl_summ_mort_macro
```

```{r}
tbl_summ_mort_macro_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:2) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/mortalidad_macrorregion.xlsx",
                       open = FALSE)

tbl_summ_mort_macro %>%
  gt::gtsave("02_output/tables/mortalidad_departamento_macrorregion.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/mortalidad_departamento_edad.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/mortalidad_departamento_edad.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Macrorregiones - Letalidad

```{r}
letalidad_macro_100k_sf <- peru_macro_sf %>%
  left_join(
    results_letalidad_overall$macrorregion$index
  ) 
```

```{r}
#| fig-cap: !expr paste0("Letalidad COVID-19 por 100K habitantes de acuerdo a la macrorregión ", range_fallecidos$fecha)
#| label: fig-let-macro-map

let_macro_map <- letalidad_macro_100k_sf %>%
  mutate(
    macrorregion = str_wrap(macrorregion, width = 15)
  ) %>% 
  ggplot() +
  geom_sf(aes(fill = rate)) +
  ggsflabel::geom_sf_label_repel(aes(label = macrorregion),
                                 force = 100,
                                 seed = 2024) +
  guides(
    fill = guide_legend(
      str_wrap("Lethality per 100 cases", 
               20)
    )
  ) +
  scale_fill_distiller(palette = "YlOrBr",
                       direction = 1) + 
  theme_void()

let_macro_map
```

## Distrital y pobreza

```{r}
#| fig-cap: !expr paste0("Mortalidad COVID-19 de acuerdo al cuartil de pobreza", range_fallecidos$fecha)
#| label: fig-mort-pobr-area
#| fig-width: 10

mort_pobr_area_plot <- 
  results_per_month_mortalidad$cuartil_pobreza$index %>% 
  drop_na() %>% 
  mutate(
    time = make_date(anio, mes, "1")
  ) %>% 
  ggplot(
    aes(x = time,
        y = rate,
        fill = cuartil_pobreza)
  ) +
  geom_area(position = 'stack') +
  # innovar::scale_fill_innova("jama") +
  scale_fill_manual(
    values = c("#2C9BB4", "#0B4E60", "#DAE2BC", "#8C3A33")
  ) +
  scale_x_date(
    limits = c(ymd("2020-03-01"), ymd("2023-07-01")),
    breaks = seq(from = ceiling_date(ymd("2020-03-01"), "month"),
                 to = floor_date(ymd("2023-07-01"), "month"),
                 by = "2 months"), 
    # date_breaks = "2 month",
    date_labels = "%b %Y"
  ) +
  labs(
    y = NULL,
    x = NULL,
    title = "Monthly mortality COVID-19 under 18 years old",
    fill = "Poverty Quartiles"
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2020-03-01"), 
    xmax = as.Date("2020-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2020-07-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "First wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2020-11-01"), 
    xmax = as.Date("2021-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2021-05-10"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Second wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2021-11-01"), 
    xmax = as.Date("2022-04-30"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-02-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Third \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2022-05-01"), 
    xmax = as.Date("2022-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-08-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Fourth \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2022-11-01"), 
    xmax = as.Date("2023-07-01"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2023-03-01"),
    y = Inf, # Posición vertical del texto
    label = "Fifth wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  theme_classic() +
  theme(
    axis.text.x = element_text(
      angle = 75,
      vjust = 0.5
    ),
    axis.ticks.x = element_line()
  ) 

mort_pobr_area_plot
```


```{r}
#| fig-cap: !expr paste0("Letalidad COVID-19 de acuerdo al cuartil de pobreza", range_fallecidos$fecha)
#| label: letal-mort-pobr-area
#| fig-width: 10

letal_pobr_area_plot <- 
  results_per_month_letalidad$cuartil_pobreza$index %>% 
  drop_na() %>% 
  mutate(
    time = make_date(anio, mes, "1")
  ) %>% 
  filter(
    time > ymd("2020-03-01")
  ) |> 
  ggplot(
    aes(x = time,
        y = rate,
        fill = cuartil_pobreza)
  ) +
  geom_area(position = 'stack') +
  # innovar::scale_fill_innova("jama") +
  scale_fill_manual(
    values = c("#2C9BB4", "#0B4E60", "#DAE2BC", "#8C3A33")
  ) +
  scale_x_date(
    limits = c(ymd("2020-03-01"), ymd("2023-07-01")),
    breaks = seq(from = ceiling_date(ymd("2020-03-01"), "month"),
                 to = floor_date(ymd("2023-07-01"), "month"),
                 by = "2 months"), 
    # date_breaks = "2 month",
    date_labels = "%b %Y"
  ) +
  labs(
    y = NULL,
    x = NULL,
    title = "Monthly lethality COVID-19 under 18 years old",
    fill = "Poverty Quartiles"
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2020-03-01"), 
    xmax = as.Date("2020-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2020-07-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "First wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2020-11-01"), 
    xmax = as.Date("2021-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2021-05-10"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Second wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2021-11-01"), 
    xmax = as.Date("2022-04-30"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-02-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Third \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2022-05-01"), 
    xmax = as.Date("2022-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-08-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Fourth \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2022-11-01"), 
    xmax = as.Date("2023-07-01"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2023-03-01"),
    y = Inf, # Posición vertical del texto
    label = "Fifth wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  theme_classic() +
  theme(
    axis.text.x = element_text(
      angle = 75,
      vjust = 0.5
    ),
    axis.ticks.x = element_line()
  ) 

letal_pobr_area_plot
```


## Distrito Lima

```{r}
#| fig-cap: !expr paste0("Mortalidad COVID-19 de acuerdo a los distritos de lima metropolitana ", range_fallecidos$fecha)
#| label: fig-mort-lima-bar
#| fig-height: 9

total_mort_distr_lim_callao <- results_mortalidad$distrito$index_pre %>% 
  filter(provincia %in% c("LIMA", "CALLAO")) %>% 
  pull(count) %>% 
  sum()

results_mortalidad$distrito$index %>% 
  filter(provincia %in% c("LIMA", "CALLAO")) %>% 
  mutate(
    distrito = str_to_title(distrito), 
    distrito = fct_reorder(distrito, rate)
  ) %>% 
  ggplot(aes(y = distrito, x = rate)) +
  geom_bar(stat = "identity",
           color = "black",
           fill = "grey80") +
  theme_bw() +
  labs(title = "Mortality COVID-19 in Metropolitan Lima", 
       x = NULL, 
       y = "Mean monthly COVID-19 mortality per 100 000 inhabitants",
       caption = paste0("Total deaths reported: ", 
                        scales::number(total_mort_distr_lim_callao))) +
  innovar::scale_fill_innova("blue_fall") +
  theme(
    legend.position = "bottom"
  )
```

```{r}
mort_lima_distr_1k_sf <- Peru %>% 
  inner_join(
    results_mortalidad$distrito$index %>% 
      filter(provincia %in% c("LIMA", "CALLAO"))
  )
```

```{r}
#| fig-cap: !expr paste0("Mortalidad COVID-19 de acuerdo a los distritos de lima metropolitana ", range_fallecidos$fecha)
#| label: fig-mort-lima-map

mort_lima_distr_1k_sf %>%
  mutate(
    label = case_when(
      rate > 10 ~ str_to_title(distrito)
    )
  ) %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap("Mean monthly COVID-19 mortality per 100 000 inhabitants", 
               20)
    )
  ) +
  ggsflabel::geom_sf_label_repel(aes(label = label),
                                 force = 10,
                                 size = 2.45,
                                 box.padding = 0.75,
                                 seed = 2024) +
  # scale_fill_distiller(palette = "YlOrBr",
  #                      direction = 1) + 
  scale_fill_gradientn(
    colors = scales::pal_brewer("seq", "Oranges", 1)(5),
    values = scales::rescale(c(0, 15, 30, 50, max(mort_lima_distr_1k_sf$rate, na.rm = TRUE))),
    breaks = c(0, 15, 30, 50, max(mort_lima_distr_1k_sf$rate, na.rm = TRUE)),
    labels = c("0", "15", "30", "50", "Max"),
    limits = c(0,
               max(mort_lima_distr_1k_sf$rate, na.rm = TRUE))
  ) +
  theme_void()
```

```{r}
#| tbl-cap: !expr paste0("Mortalidad COVID-19 de acuerdo a los distritos de lima metropolitana ", range_fallecidos$fecha)
#| label: tbl-mort-lima

tbl_summ_mort_lima_distr_pre <- results_mortalidad$distrito$index %>% 
      filter(provincia %in% c("LIMA", "CALLAO")) %>% 
  mutate(
    across(
      departamento:distrito,
      str_to_title
    ),
    rate = formattable::digits(rate, 3)
  ) %>% 
  select(-departamento)

tbl_summ_mort_lima_distr <- tbl_summ_mort_lima_distr_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

tbl_summ_mort_lima_distr
```

```{r}
tbl_summ_mort_lima_distr_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:4) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/mortalidad_lima_distrito.xlsx",
                       open = FALSE)

tbl_summ_mort_lima_distr %>%
  gt::gtsave("02_output/tables/mortalidad_lima_distrito.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/mortalidad_lima_distrito.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/mortalidad_lima_distrito.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

# Incidencia mensual del COVID-19

```{r}
incidencia_positivos <- positivos_0723 %>% 
  mutate(
    #fecha_resultado_floor = floor_date(fecha_resultado, unit = "month"),
    mes = month(fecha_resultado_floor)
  )
```

```{r}
n_incidencia_positivos <- incidencia_positivos %>% 
  count() %>% 
  pull(n)

n_wo_id <- incidencia_positivos %>% 
  filter(is.na(id_persona)) %>% 
  count() %>% 
  pull(n)
```

Los datos utilizados para incidencia positivos tienen datos hasta julio de 2023. Se tiene registro de `r scales::number(n_incidencia_positivos)` casos positivos menores de 18 años. Para poder analizar la incidencia necesitamos saber que si un caso positivos es nuevo o es una reinfección que no se puede considerar aún como caso nuevo. Se tendrá como criterio 3 meses de espacio entre infecciones. Lamentablemente, incluso en estos datos, hay `r scales::number(n_wo_id)` registros sin `id_persona`.

Criterios para casos positivos y de reinfección:

-   Hasta 1 mes (30 días), se considera una misma infección.
-   Espacio de 3 meses a más (90 días), se considera una nueva infección (reinfección)
-   Reportar cantidad de casos que se encuentren entre 1 mes y menos de 3 meses (entre día 31 y 89), con su respectivo método.

```{r}
incidencia_positivos %>% 
  count(caso_intermedio)
```

<!-- ```{r} -->

<!-- incidencia_positivos <- incidencia_positivos %>% -->

<!--   filter(!is.na(id_persona)) %>%  -->

<!--   collect() %>%  -->

<!--   group_by(id_persona) %>% -->

<!--   mutate( -->

<!--     reinfeccion = fecha_resultado - lag(fecha_resultado) > months(3) | is.na(lag(fecha_resultado)), -->

<!--     caso_nuevo = if_else(reinfeccion | row_number() == 1, TRUE, FALSE) -->

<!--   ) %>% -->

<!--   ungroup() -->

<!-- casos_nuevos <- incidencia_positivos %>% -->

<!--   filter(caso_nuevo) -->

<!-- ``` -->

```{r}
# incidencia_positivos %>% 
#   filter(anio == 2020,
#          !is.na(id_persona)) %>% 
#   group_by(fecha_resultado_floor) %>%
#   summarise(
#     casos_sin_reinfecciones = n_distinct(id_persona),
#     casos_con_reinfeccioens = n()
#   ) %>% 
#   ungroup() %>% 
#   arrange(fecha_resultado_floor) %>% 
#   collect()
# 
# incidencia_positivos %>% 
#   filter(anio == 2020, !is.na(id_persona)) %>% 
#   group_by(fecha_resultado_floor) %>%
#   add_count(id_persona) %>% 
#   filter(n > 1) %>% 
#   select(departamento:fecha_resultado, id_persona)

# incidencia_positivos %>% 
#   filter(anio == 2020) %>% 
#   count(fecha_resultado_floor) %>%
#   ungroup() %>% 
#   arrange(fecha_resultado_floor) %>% 
#   collect()
#   
#   count(
#     across(c(fecha_resultado_floor, departamento:distrito, grupo_edad, sexo))
#   )
```

# Tasa de hospitalización

<!-- Que no sea hospitalización promedio mensual, sino, de forma global -->

```{r}
hospitalizados <- hosp_final %>%
  filter(edad < 18,
         fecha_ingreso_hosp <= max_positivos_0723) %>% 
  mutate(
    mes = month(fecha_ingreso_hosp),
    overall = "overall"
  ) 

range_hospi <- hospitalizados %>% 
  summarise(
    min = min(fecha_ingreso_hosp),
    max = max(fecha_ingreso_hosp)
  ) %>% 
  collect() %>% 
  mutate(
    min = format(min, "%d-%m-%Y"),
    max = format(max, "%d-%m-%Y"),
    fecha = paste0("(", min, " a ", max, ")")
  )
```

## Vista general

**Cálculo de Tasa de hospitalización**

```{r}
levels <- list(
  c("overall"),
  c("grupo_edad"),
  c("sexo"),
  c("ola_covid"),
  c("ola_covid", "grupo_edad"),
  c("ola_covid", "sexo"),
  c("ola_covid", "region"),
  c("ola_covid", "macrorregion"),
  c("departamento"),
  c("departamento", "sexo"),
  c("departamento", "grupo_edad"),
  c("departamento", "sexo", "grupo_edad"),
  c("departamento", "ola_covid"),
  c("distrito"),
  c("distrito", "sexo"),
  c("distrito", "grupo_edad"),
  c("distrito", "cuartil_pobreza"),
  c("grupo_edad", "cuartil_pobreza"),
  c("sexo", "cuartil_pobreza"),
  c("ola_covid", "cuartil_pobreza"),
  c("ola_covid", "cuartil_nbi"),
  c("cuartil_pobreza"),
  c("cuartil_nbi"),
  c("region"),
  c("macrorregion")
)

# Calcular la tasa de hospitalizaciones
results_hospi <- index_calculator(
  data = hospitalizados, 
  denom = positivos_0723 %>% 
    mutate(cantidad = 1), 
  levels = levels, 
  type = "overall",
  std = 100
)

results_per_month_hospi <- index_calculator(
  data = hospitalizados, 
  denom = positivos_0723 %>% 
    mutate(cantidad = 1), 
  levels = levels, 
  type = "per_month",
  std = 100
)
```

```{r}
hospi_total_df <- hospitalizados %>% 
  select(edad, overall, grupo_edad, sexo, ola_covid, region, 
         macrorregion, cuartil_pobreza, cuartil_nbi) %>% 
  collect()

labelled::var_label(hospi_total_df) <- list(
  overall = "Overall",
  grupo_edad = "Age Group",
  sexo = "Sex",
  ola_covid = "COVID-19 wave",
  region = "Natural Regions",
  macrorregion = "Macroregions",
  cuartil_pobreza = "Poverty level",
  cuartil_nbi = "Unsatisfied basic needs level"
)
```

```{r}
#| label: tbl-summ-hospi
#| tbl-cap: !expr paste0("Descriptivos de Hospitalizados menores de 18 años ", range_hospi$fecha)

tbl_summ_hospi_g <- hospi_total_df %>% 
  select(overall, grupo_edad, sexo, ola_covid, region, 
         macrorregion, cuartil_pobreza, cuartil_nbi) %>% 
  gtsummary::tbl_summary(
    type = list(
      overall ~ "dichotomous"
    ),
    value = list(overall = "overall"),
    missing_text = "No data"
  ) %>% 
  gtsummary::modify_header(label = "**Variables**") %>%
  gtsummary::bold_labels()

tbl_summ_hospi_c2 <- results_hospi$overall$index %>% 
  mutate(label = "Overall") %>% 
  rename(variable = overall) %>% 
  bind_rows(
    results_hospi$grupo_edad$index %>% 
      mutate(variable = "grupo_edad") %>% 
      rename(label = grupo_edad),
    results_hospi$sexo$index %>% 
      mutate(variable = "sexo") %>% 
      rename(label = sexo),
    results_hospi$ola_covid$index %>% 
      mutate(variable = "ola_covid") %>% 
      rename(label = ola_covid),
    results_hospi$region$index %>% 
      mutate(variable = "region") %>% 
      rename(label = region),
    results_hospi$macrorregion$index %>% 
      mutate(variable = "macrorregion") %>% 
      rename(label = macrorregion),
    results_hospi$cuartil_pobreza$index %>% 
      mutate(variable = "cuartil_pobreza") %>% 
      rename(label = cuartil_pobreza),
    results_hospi$cuartil_nbi$index %>% 
      mutate(variable = "cuartil_nbi") %>% 
      rename(label = cuartil_nbi)
  ) %>% 
  drop_na() %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  )

tbl_summ_hospi_final <- tbl_summ_hospi_g %>%
  gtsummary::modify_table_body(
    ~ .x %>% 
      left_join(tbl_summ_hospi_c2) %>% 
      mutate(
        label = case_when(
          variable %in% c("cuartil_pobreza", "cuartil_nbi") &
            label == "1" ~ "1 (0-19.9%)",
          variable %in% c("cuartil_pobreza", "cuartil_nbi") & 
            label == "2" ~ "2 (20-39.9%)",
          variable %in% c("cuartil_pobreza", "cuartil_nbi") & 
            label == "3" ~ "3 (40-59.9%)",
          variable %in% c("cuartil_pobreza", "cuartil_nbi") &
            label == "4" ~ "4 (60-100%)",
          .default = label
        )
      )
  ) %>% 
  gtsummary::modify_header(
    rate ~ "**Hospital admission rate per 100 cases**"
  ) %>%
  gtsummary::modify_fmt_fun(
    rate ~ scales::number_format(accuracy = 0.001)
  )

tbl_summ_hospi_final
```

```{r}
tbl_summ_hospi_final %>% 
  gtsummary::as_hux_xlsx("02_output/tables/hospitalizacion_descriptivos.xlsx")

tbl_summ_hospi_final %>%
  gtsummary::as_flex_table() %>%
  flextable::save_as_docx(path = "02_output/tables/hospitalizacion_descriptivos.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/hospitalizacion_descriptivos.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/hospitalizacion_descriptivos.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

<!-- Cálculo de Comparaciones Múltiples [Hospitalizaciones] -->

```{r}
hospital_comparisons <- compare_levels_per_month(
  results_per_month = results_per_month_hospi,
  levels = comparing_levels
)

# Extraer las tablas `cld` como tibbles, renombrando y añadiendo una columna "level"
hospital_comparisons_cld_pre <- lapply(names(hospital_comparisons), function(level) {
  cld <- hospital_comparisons[[level]]$cld
  
  if (!is.null(cld)) {
    cld <- cld |> 
      as_tibble() |> 
      rename(group = 1) |>  # Renombra la primera columna a 'group'
      mutate(level = level) |> # Añade una columna 'level' con el nombre de la lista
      relocate(level) 
  }
  
  cld
}) |> 
  bind_rows()

hospital_comparisons_cld_pre <- 
  hospital_comparisons_cld_pre |> 
  mutate(
    across(
      where(is.numeric),
      ~ formattable::digits(.x, 3)
    )
  ) 
```

```{r}
#| label: tbl-compar-hospital
#| tbl-cap: !expr paste0("Comparaciones múltiples de hospitalidad", range_pos_18_0723$fecha)
#| 
hospital_comparisons_cld <- hospital_comparisons_cld_pre |> 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

hospital_comparisons_cld
```

```{r}
hospital_comparisons_cld_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:8) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/hospital_comparacion_multiple.xlsx",
                       open = FALSE)

hospital_comparisons_cld %>%
  gt::gtsave("02_output/tables/hospital_comparacion_multiple.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/hospital_comparacion_multiple.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/hospital_comparacion_multiple.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

<!-- Comparaciones de hospitalizados con respecto a olas -->

```{r}
#| label: tbl-summ-hospi-olas
#| tbl-cap: !expr paste0("Descriptivos de hospitalizados menores de 18 años por olas ", range_hospi$fecha)

tbl_summ_hospi_olas <- hospi_total_df %>% 
  select(overall, grupo_edad, sexo, ola_covid, region, 
         macrorregion, cuartil_pobreza, cuartil_nbi) %>% 
  gtsummary::tbl_summary(
    by = ola_covid,
    type = list(
      overall ~ "dichotomous"
    ),
    value = list(overall = "overall"),
    missing_text = "No data"
  ) %>% 
  gtsummary::modify_header(label = "**Variables**") %>%
  gtsummary::bold_labels()


ola_covid_hospi_gtsumm <- results_hospi$ola_covid$index %>% 
  rename(stat = ola_covid) %>% 
  mutate(label = "Overall",
         variable = "overall") %>% 
  bind_rows(
    results_hospi$ola_covid_grupo_edad$index %>%
      mutate(variable = "grupo_edad") %>% 
      rename(label = grupo_edad,
             stat = ola_covid),
    results_hospi$ola_covid_sexo$index %>% 
      mutate(variable = "sexo") %>% 
      rename(label = sexo,
             stat = ola_covid),
    results_hospi$ola_covid_region$index %>% 
      mutate(variable = "region") %>% 
      rename(label = region,
             stat = ola_covid),
    results_hospi$ola_covid_macrorregion$index %>% 
      mutate(variable = "macrorregion") %>%
      rename(label = macrorregion,
             stat = ola_covid),
    results_hospi$ola_covid_cuartil_pobreza$index %>% 
      mutate(variable = "cuartil_pobreza") %>% 
      rename(label = cuartil_pobreza,
             stat = ola_covid),
    results_hospi$ola_covid_cuartil_nbi$index %>% 
      mutate(variable = "cuartil_nbi") %>% 
      rename(label = cuartil_nbi,
             stat = ola_covid)
  ) %>% 
  mutate(
    rate = formattable::digits(rate, 3),
    stat = as.character(stat),
    stat = case_match(
      stat,
      "1 (Mar – Oct 2020)" ~ "stat_1",
      "2 (Nov 2020 – Oct 2021)" ~ "stat_2",
      "3 (Nov 2021 – Apr 2022)" ~ "stat_3",
      "4 (May 2022 – Oct 2022)" ~ "stat_4",
      "5 (Nov 2022 – Jun 2023)" ~ "stat_5"
    )
  ) 


tbl_summ_hospi_olas$table_body <- tbl_summ_hospi_olas$table_body %>%
  mutate(
    across(
      gtsummary::all_stat_cols(stat_0 = FALSE),
      ~ str_remove_all(., "\\s*\\(.*?\\)\\s*")
    )
  ) %>% 
  pivot_longer(
    cols = stat_1:stat_5,
    names_to = "stat",
    values_to = "freq"
  ) %>% 
  left_join(
    ola_covid_hospi_gtsumm
  ) %>% 
  mutate(
    value = str_c(freq, ifelse(is.na(rate), "", str_c(" (", rate, ")")))
  ) %>% 
  select(-c(freq, rate)) %>% 
  pivot_wider(
    names_from = stat,
    values_from = value
  )


tbl_summ_hospi_olas <- tbl_summ_hospi_olas %>% 
  gtsummary::modify_footnote(
    gtsummary::all_stat_cols() ~ "n (HAR per 100 cases)"
  )

tbl_summ_hospi_olas
```

```{r}
tbl_summ_hospi_olas %>% 
  gtsummary::as_hux_xlsx("02_output/tables/hospitalizados_olas_descriptivos.xlsx")

tbl_summ_hospi_olas %>%
  gtsummary::as_flex_table() %>%
  flextable::save_as_docx(path = "02_output/tables/hospitalizados_olas_descriptivos.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/hospitalizados_olas_descriptivos.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/hospitalizados_olas_descriptivos.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Gráfico de hospitalización (serie de tiempo) - Menores de 18 años

```{r}
#| fig-cap: !expr paste0("Hospitalización COVID-19 de acuerdo a la edad y ola covid ", range_hospi$fecha)
#| label: fig-hospi-edad-ola-area
#| fig-width: 10

hospi_count_by_age_area <- hospitalizados %>% 
  to_arrow() %>% 
  mutate(
    fecha_round = ceiling_date(fecha_ingreso_hosp, unit = "month")
  ) %>%
  compute() %>% 
  to_duckdb() %>% 
  mutate(
    fecha_round = case_when(
      fecha_round == as.Date("2024-01-01") ~  as.Date("2023-12-31"),
      .default = fecha_round
    )
  ) %>%
  count(grupo_edad, fecha_round) %>% 
  collect()

hospi_age_area_plot <- hospi_count_by_age_area %>% 
  bind_rows(
    tibble(
      grupo_edad = "Under 1 year",
      fecha_round = seq(from = as.Date("2020-04-01"), to = as.Date("2023-07-01"), by = "month"),
      n = 0
    )
  ) %>% 
  mutate(
    grupo_edad = factor(
      grupo_edad,
      levels = c("Under 1 year", "1 - 5 years",
                 "6 - 11 years", "12 - 17 years")
    )
  ) %>% 
  ggplot(
    aes(
      x = fecha_round,
      y = n,
      fill = grupo_edad
    )
  ) +
  geom_area(position = 'stack') +
  # innovar::scale_fill_innova("jama") +
  scale_fill_manual(
    values = c("#2C9BB4", "#0B4E60", "#DAE2BC", "#8C3A33")
  ) +
  scale_x_date(
    limits = c(ymd("2020-03-07"), ymd("2023-07-01")),
    breaks = seq(from = ceiling_date(ymd("2020-03-07"), "month"),
                 to = floor_date(ymd("2023-07-01"), "month"),
                 by = "2 months"), 
    # date_breaks = "2 month",
    date_labels = "%b %Y"
  ) +
  labs(
    y = NULL,
    x = NULL,
    title = "COVID-19 hospital admissions",
    fill = "Age Group"
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2020-04-01"), 
    xmax = as.Date("2020-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2020-07-15"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "First wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2020-11-01"), 
    xmax = as.Date("2021-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2021-05-10"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Second wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2021-11-01"), 
    xmax = as.Date("2022-04-30"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
    annotate(
    geom = "text", 
    x = as.Date("2022-02-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Third \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2022-05-01"), 
    xmax = as.Date("2022-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-08-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Fourth \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2022-11-01"), 
    xmax = as.Date("2023-07-01"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2023-03-01"),
    y = Inf, # Posición vertical del texto
    label = "Fifth wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  theme_classic() +
  theme(
    axis.text.x = element_text(
      angle = 75,
      vjust = 0.5
    ),
    axis.ticks.x = element_line()
  ) 

hospi_age_area_plot
```

## Gráfico Tasa de hospitalización (serie de tiempo)

```{r}
per_month_hospi_under_1 <- tibble(
  grupo_edad = "Under 1 year",
  time = seq(from = as.Date("2020-03-01"), to = as.Date("2023-07-01"), by = "month"),
  rate = 0
)

hospi_rate_by_age_area <- results_per_month_hospi$grupo_edad$index %>% 
  mutate(
    time = make_date(anio, mes, "1")
  ) %>% 
  bind_rows(
    per_month_hospi_under_1
  ) %>%  
  mutate(
    grupo_edad = factor(
      grupo_edad,
      levels = c("Under 1 year", "1 - 5 years",
                 "6 - 11 years", "12 - 17 years")
    )
  ) %>% 
  ggplot(
    aes(
      x = time,
      y = rate,
      fill = grupo_edad
    )
  ) +
  geom_area(position = 'stack') +
  scale_fill_manual(
    values = c("#2C9BB4", "#0B4E60", "#DAE2BC", "#8C3A33")
  ) +
  scale_x_date(
    limits = c(ymd("2020-03-01"), ymd("2023-07-01")),
    breaks = seq(from = ceiling_date(ymd("2020-03-01"), "month"),
                 to = floor_date(ymd("2023-07-01"), "month"),
                 by = "2 months"), 
    # date_breaks = "2 month",
    date_labels = "%b %Y"
  ) +
  labs(
    y = NULL,
    x = NULL,
    title = "Monthly COVID-19 hospital admission rate per 100 cases",
    fill = "Age Group"
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2020-03-01"), 
    xmax = as.Date("2020-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2020-07-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "First wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2020-11-01"), 
    xmax = as.Date("2021-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2021-05-10"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Second wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2021-11-01"), 
    xmax = as.Date("2022-04-30"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
    annotate(
    geom = "text", 
    x = as.Date("2022-02-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Third \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2022-05-01"), 
    xmax = as.Date("2022-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-08-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Fourth \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2022-11-01"), 
    xmax = as.Date("2023-07-01"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2023-03-01"),
    y = Inf, # Posición vertical del texto
    label = "Fifth wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  theme_classic() +
  theme(
    axis.text.x = element_text(
      angle = 75,
      vjust = 0.5
    ),
    axis.ticks.x = element_line()
  ) 

hospi_rate_by_age_area
```

## Pirámide frecuencia - Edades

```{r}
#| dpi: 300
#| fig-width: 8

hospi_total_df %>% 
  mutate(
    edad = factor(edad)
  ) %>% 
  age_pyramid(
    age_group = "edad",
    split_by = "sexo"
  ) +
  scale_fill_brewer(palette = "Dark2") +
  labs(
    x = "Edad",
    fill = "Sexo",
    y = "Frecuencia"
  ) +
  theme_classic()
```

```{r}
#| dpi: 300
#| fig-width: 6

hospi_total_df %>% 
  age_pyramid(
    age_group = "grupo_edad",
    split_by = "sexo"
  ) +
  scale_fill_brewer(palette = "Dark2") +
  labs(
    x = "Grupo de Edad",
    fill = "Sexo",
    y = "Frecuencia"
  ) +
  theme_classic()
```

## Pirámide proporcional - Edades

```{r}
#| dpi: 300
#| fig-width: 8

hospi_total_df %>% 
  mutate(
    edad = factor(edad)
  ) %>% 
  age_pyramid(
    age_group = "edad",
    split_by = "sexo",
    proportional = TRUE
  ) +
   scale_fill_brewer(palette = "Dark2") +
  labs(
    x = "Edad",
    fill = "Sexo",
    y = "Porcentaje"
  ) +
  theme_classic()
```

```{r}
#| dpi: 300
#| fig-width: 6

hospi_total_df %>% 
  age_pyramid(
    age_group = "grupo_edad",
    split_by = "sexo",
    proportional = TRUE
  ) +
   scale_fill_brewer(palette = "Dark2") +
  labs(
    x = "Grupo de Edad",
    fill = "Sexo",
    y = "Porcentaje"
  ) +
  theme_classic()
```

## Tasa de Hospitalización

```{r}
#| fig-cap: !expr paste0("Tasa de Hospitalización por 1000 positivos COVID-19 de acuerdo al departamento, sexo y edad ", range_hospi$fecha)
#| label: fig-tasa-hospi-dep-sex-edad-area
#| fig-width: 12
#| fig-height: 7

results_hospi$departamento_sexo_grupo_edad$index %>% 
  mutate(
    departamento = str_to_title(departamento)
  ) %>% 
  ggplot(
    aes(
      x = departamento,
      y = rate,
      fill = sexo,
      group = sexo
    )
  ) +
  geom_area() +
  facet_wrap(vars(grupo_edad),
             nrow = 2) +
  innovar::scale_fill_innova("jama") +
  labs(
    y = NULL,
    x = NULL,
    title = "Hospital admission rate per 100 cases",
    fill = "Sex"
  ) + 
  theme_bw() +
  theme(
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1
    )
  )
```

## Departamento

```{r}
tasa_hospi_dep_sf <- peru_dep_sf %>% 
  rename(departamento = dep) %>% 
  inner_join(
    results_hospi$departamento$index
  ) 
```

```{r}
#| fig-cap: !expr paste0("Tasa de Hospitalización por 100 positivos COVID-19 de acuerdo al departamento ", range_hospi$fecha)
#| label: fig-tasa-hospi-dep-map
#| fig-width: 5
#| fig-height: 5

tasa_hospi_dep_map <- tasa_hospi_dep_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap(
        "Mean COVID-19 hospital admission rate per 100 COVID-19 cases",
        20
      )
    )
  ) +
  scale_fill_distiller(palette = "Oranges", # Old: "YlOrBr",
                       direction = 1) + 
  theme_void()

tasa_hospi_dep_map
```

```{r}
#| tbl-cap: !expr paste0("Tasa de Hospitalización COVID-19 de acuerdo al departamento ", range_pos_18_0723$fecha)
#| label: tbl-tasa-hospi-dep

tbl_summ_hosp_dep_pre <- results_hospi$departamento$index %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  ) %>% 
  mutate(
    departamento = str_to_title(departamento)
  ) 

tbl_summ_hosp_dep <- tbl_summ_hosp_dep_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

tbl_summ_hosp_dep
```

```{r}
tbl_summ_hosp_dep_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:2) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/tasa_hospitalizacion_departamento.xlsx",
                       open = FALSE)

tbl_summ_hosp_dep %>%
  gt::gtsave("02_output/tables/tasa_hospitalizacion_departamento.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/tasa_hospitalizacion_departamento.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/tasa_hospitalizacion_departamento.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Departamento y Sexo

```{r}
tasa_hospi_dep_sex_sf <- peru_dep_sf %>% 
  rename(departamento = dep) %>% 
  inner_join(
    results_hospi$departamento_sexo$index
  ) %>%
  ungroup()
```

```{r}
#| fig-cap: !expr paste0("Tasa de Hospitalización por 100 positivos COVID-19 de acuerdo al departamento y sexo", range_hospi$fecha)
#| label: fig-tasa-hospi-dep-sex-map

tasa_hospi_dep_sex_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap(
        "Mean COVID-19 hospital admission rate per 100 COVID-19 cases",
        20
      )
    )
  ) +
  scale_fill_distiller(palette = "Oranges",  #"YlOrBr",
                       direction = 1) + 
  facet_wrap(vars(sexo)) + 
  theme_void()
```

```{r}
#| fig-cap: !expr paste0("Tasa de Hospitalización por 100 positivos COVID-19 de acuerdo al departamento y sexo ", range_hospi$fecha)
#| label: fig-tasa-hospi-dep-sex-bar

results_hospi$departamento_sexo$index %>% 
  mutate(
    departamento = str_to_title(departamento),
    departamento = fct_reorder(departamento, rate, .desc = TRUE)
    # departamento = fct_reorder2(departamento, sexo, tasa_hospi_mes)
  ) %>%
  ggplot(aes(x = departamento, y = rate, fill = sexo)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  theme_bw() +
  labs(title = "Hospital admissions COVID-19 by department and sex", 
       x = NULL, 
       y = "Mean COVID-19 hospital admission rate per 100 COVID-19 cases",
       fill = "Sex",
       caption = paste0("Total hospital admissions: ", 
                        scales::number(sum(pull(results_hospi$departamento_sexo$index_pre, count))))) +
  innovar::scale_fill_innova("blue_fall") +
  theme(
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1
    )
  )
```

## Departamento y Edad

```{r}
tasa_hospi_dep_edad_sf <- peru_dep_sf %>% 
  rename(departamento = dep) %>% 
  inner_join(
    results_hospi$departamento_grupo_edad$index
  ) %>%
  ungroup()
```

```{r}
#| fig-cap: !expr paste0("Tasa de Hospitalización por 100 positivos COVID-19 de acuerdo al departamento y grupo de edad", range_hospi$fecha)
#| label: fig-tasa-hospi-dep-edad-map

tasa_hospi_dep_edad_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap(
        "Mean COVID-19 hospital admission rate per 100 COVID-19 cases",
        20
      )
    )
  ) +
  scale_fill_distiller(palette = "Oranges", # Old: "YlOrBr",
                       direction = 1) + 
  facet_wrap(vars(grupo_edad),
             nrow = 1) + 
  theme_void()
```

## Departamento y Olas

```{r}
tasa_hospi_dep_olas_mother <- crossing(
  ola_covid = levels(results_hospi$departamento_ola_covid$index$ola_covid),
  departamento = unique(results_hospi$departamento_ola_covid$index$departamento)
)

tasa_hospi_dep_olas <- tasa_hospi_dep_olas_mother %>% 
  left_join(results_hospi$departamento_ola_covid$index) %>% 
  bind_rows(
    results_hospi$departamento$index %>% 
      mutate(ola_covid = "Overall")
  ) %>% 
  mutate(
    ola_covid = fct_recode(
      ola_covid,
      "1st wave" = "1 (Mar – Oct 2020)",
      "2nd wave" = "2 (Nov 2020 – Oct 2021)",
      "3rd wave" = "3 (Nov 2021 – Apr 2022)",
      "4th wave" = "4 (May 2022 – Oct 2022)",
      "5th wave" = "5 (Nov 2022 – Jun 2023)"
    ),
    ola_covid = fct_relevel(ola_covid,
                            "Overall")
  ) %>% 
  arrange(ola_covid)

tasa_hospi_dep_olas_sf <- peru_dep_sf %>%
  rename(departamento = dep) %>% 
  left_join(
    tasa_hospi_dep_olas
  ) %>%
  ungroup()
```

```{r}
#| fig-cap: !expr paste0("Tasa de Hospitalización por 1000 positivos COVID-19 de acuerdo al departamento y ola ", range_hospi$fecha)
#| label: fig-tasa-hospi-dep-ola-map
#| fig-width: 13
#| fig-height: 3

tasa_hospi_dep_olas_map <- tasa_hospi_dep_olas_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap(
        "Mean COVID-19 hospital admission rate per 100 COVID-19 cases",
        20
      )
    )
  ) +
  scale_fill_distiller(palette = "Oranges", #Old: "YlOrBr",
                       direction = 1, na.value = "white") + 
  facet_wrap(vars(ola_covid), nrow = 1) + 
  theme_void()

tasa_hospi_dep_olas_map
```

## Region

```{r}
tasa_hospi_reg_sf <- peru_reg_sf %>% 
  inner_join(
    results_hospi$region$index
  ) 
```

```{r}
#| fig-cap: !expr paste0("Tasa de Hospitalización por 100 positivos COVID-19 de acuerdo a la región ", range_hospi$fecha)
#| label: fig-tasa-hospi-reg-map
#| fig-width: 5
#| fig-height: 5

tasa_hospi_reg_map <- tasa_hospi_reg_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap(
        "Mean COVID-19 hospital admission rate per 100 COVID-19 cases",
        20
      )
    )
  ) +
  ggsflabel::geom_sf_label(aes(label = region)) +
  scale_fill_distiller(palette = "Oranges", #Old: "YlOrBr",
                       direction = 1) + 
  theme_void()

tasa_hospi_reg_map
```

## Macrorregión

```{r}
tasa_hospi_macro_sf <- peru_macro_sf %>% 
  inner_join(
    results_hospi$macrorregion$index
  ) 
```

```{r}
#| fig-cap: !expr paste0("Tasa de Hospitalización por 1000 positivos COVID-19 de acuerdo a la macrorregión ", range_hospi$fecha)
#| label: fig-tasa-hospi-macro-map
#| fig-width: 5
#| fig-height: 5

hospi_macro_map <- tasa_hospi_macro_sf %>%
  mutate(
    macrorregion = str_wrap(macrorregion, width = 15)
  ) %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  ggsflabel::geom_sf_label_repel(aes(label = macrorregion),
                                 force = 100,
                                 seed = 2024) +
  guides(
    fill = guide_legend(
      str_wrap("Mean COVID-19 hospital admission rate per 100 COVID-19 cases", 
               20)
    )
  ) +
  scale_fill_distiller(palette = "Oranges", # Old: YlOrBr
                       direction = 1) + 
  theme_void()

hospi_macro_map
```

## Distrital y pobreza

```{r}
#| fig-cap: !expr paste0("Mortalidad COVID-19 de acuerdo al cuartil de pobreza", range_fallecidos$fecha)
#| label: fig-hospi-pobr-area
#| fig-width: 10

hospi_pobr_area_plot <- 
  results_per_month_hospi$cuartil_pobreza$index %>% 
  drop_na() %>% 
  mutate(
    time = make_date(anio, mes, "1")
  ) %>% 
  ggplot(
    aes(x = time,
        y = rate,
        fill = cuartil_pobreza)
  ) +
  geom_area(position = 'stack') +
  # innovar::scale_fill_innova("jama") +
  scale_fill_manual(
    values = c("#2C9BB4", "#0B4E60", "#DAE2BC", "#8C3A33")
  ) +
  scale_x_date(
    limits = c(ymd("2020-03-01"), ymd("2023-07-01")),
    breaks = seq(from = ceiling_date(ymd("2020-03-01"), "month"),
                 to = floor_date(ymd("2023-07-01"), "month"),
                 by = "2 months"), 
    # date_breaks = "2 month",
    date_labels = "%b %Y"
  ) +
  labs(
    y = NULL,
    x = NULL,
    title = "Monthly COVID-19 hospital admission",
    fill = "Poverty Quartiles"
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2020-03-01"), 
    xmax = as.Date("2020-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2020-07-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "First wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2020-11-01"), 
    xmax = as.Date("2021-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2021-05-10"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Second wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2021-11-01"), 
    xmax = as.Date("2022-04-30"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-02-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Third \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2022-05-01"), 
    xmax = as.Date("2022-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-08-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Fourth \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2022-11-01"), 
    xmax = as.Date("2023-07-01"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2023-03-01"),
    y = Inf, # Posición vertical del texto
    label = "Fifth wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  theme_classic() +
  theme(
    axis.text.x = element_text(
      angle = 75,
      vjust = 0.5
    ),
    axis.ticks.x = element_line()
  ) 

hospi_pobr_area_plot
```

## Distrito Lima

```{r}
#| fig-cap: !expr paste0("Tasa de Hospitalización COVID-19 de acuerdo a los distritos de lima metropolitana ", range_hospi$fecha)
#| label: fig-tasa-hospi-lima-bar
#| fig-height: 9

total_hospi_distr_lim_callao <- results_hospi$distrito$index_pre %>% 
  filter(provincia %in% c("LIMA", "CALLAO")) %>% 
  pull(count) %>% 
  sum()

results_hospi$distrito$index %>% 
  filter(provincia %in% c("LIMA", "CALLAO")) %>% 
  drop_na() %>% 
  mutate(
    distrito = str_to_title(distrito), 
    distrito = fct_reorder(distrito, rate)
  ) %>% 
  ggplot(aes(y = distrito, x = rate)) +
  geom_bar(stat = "identity",
           color = "black",
           fill = "grey80") +
  theme_bw() +
  labs(title = "Hospital admissions COVID-19 in Metropolitan Lima", 
       x = NULL, 
       y = "Mean COVID-19 hospital admission rate per 100 COVID-19 cases",
       caption = paste0("Total hospital admissions: ", 
                        scales::number(total_hospi_distr_lim_callao))) +
  innovar::scale_fill_innova("blue_fall") +
  theme(
    legend.position = "bottom"
  )
```

```{r}
tasa_hospi_lima_distr_1k_sf <- Peru %>% 
  inner_join(
    results_hospi$distrito$index %>% 
      filter(provincia %in% c("LIMA", "CALLAO"))
  )
```

```{r}
#| fig-cap: !expr paste0("Tasa de Hospitalización COVID-19 de acuerdo a los distritos de lima metropolitana ", range_hospi$fecha)
#| label: fig-tasa-hospi-lima-map

tasa_hospi_lima_distr_1k_sf %>%
  mutate(
    label = case_when(
      rate > 25 ~ str_to_title(distrito)
    )
  ) %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap("Mean COVID-19 hospital admission rate per 100 COVID-19 cases", 
               20)
    )
  ) +
  ggsflabel::geom_sf_label_repel(aes(label = label),
                                 force = 10,
                                 size = 2.45,
                                 box.padding = 0.75,
                                 seed = 2024) +
  scale_fill_distiller(palette = "YlOrBr",
                       direction = 1) + 
  theme_void()
```

```{r}
#| tbl-cap: !expr paste0("Tasa de Hospitalización COVID-19 de acuerdo a los distritos de lima metropolitana ", range_hospi$fecha)
#| label: tbl-tasa-hospi-lima

tbl_summ_tasa_hospi_lima_distr_pre <- results_hospi$distrito$index %>% 
      filter(provincia %in% c("LIMA", "CALLAO")) %>% 
  mutate(
    across(
      departamento:distrito,
      str_to_title
    ),
    rate = formattable::digits(rate, 3)
  ) %>% 
  select(-departamento)

tbl_summ_tasa_hospi_lima_distr <- tbl_summ_tasa_hospi_lima_distr_pre %>% 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

tbl_summ_tasa_hospi_lima_distr
```

```{r}
tbl_summ_tasa_hospi_lima_distr_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:4) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/tasa_hospitalizacion_lima_distrito.xlsx",
                       open = FALSE)

tbl_summ_tasa_hospi_lima_distr %>%
  gt::gtsave("02_output/tables/tasa_hospitalizacion_lima_distrito.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/tasa_hospitalizacion_lima_distrito.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/tasa_hospitalizacion_lima_distrito.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

# Tasa de Vacunación

<!-- Agregar la tasa de vacunación global: por edad, sexo. Además, algunos gráficos que sean solo dosis 1 a 3 -->

```{r}
vacunas_18 <- vacunas_18 %>% 
  filter(dosis %in% 1:3) %>% 
  mutate(
    overall = "overall"
  ) 

range_vacu <- vacunas_18 %>% 
  summarise(
    min = min(fecha_vacunacion),
    max = max(fecha_vacunacion)
  ) %>% 
  collect() %>% 
  mutate(
    min = format(min, "%d-%m-%Y"),
    max = format(max, "%d-%m-%Y"),
    fecha = paste0("(", min, " a ", max, ")")
  )
```

## Vista General

<!-- La unidad de análisis debe ser la persona. Trabajar estas tablas en función de la cantidad de vacunados. Y en dosis, las categorías deben ser excluyentes. -->

<!-- Verificar las dosis 8 y 9 en función del uuid y el número de dosis reportado. Alguien con dosis 5, debería tener 4 registros anteriores de sus demás dosis. -->

**Cálculo de tasa de vacunación**

```{r}
levels <- list(
  c("overall"),
  c("dosis"),
  c("grupo_edad"),
  c("sexo"),
  c("ola_covid"),
  c("ola_covid", "grupo_edad"),
  c("departamento"),
  c("departamento", "sexo"),
  c("departamento", "grupo_edad"),
  c("departamento", "sexo", "grupo_edad"),
  c("departamento", "ola_covid"),
  c("distrito"),
  c("distrito", "sexo"),
  c("distrito", "grupo_edad"),
  c("distrito", "cuartil_pobreza"),
  c("grupo_edad", "cuartil_pobreza"),
  c("sexo", "cuartil_pobreza"),
  c("ola_covid", "cuartil_pobreza"),
  c("cuartil_pobreza"),
  c("ola_covid", "cuartil_nbi"),
  c("cuartil_nbi"),
  c("region"),
  c("region", "grupo_edad"),
  c("macrorregion"),
  c("macrorregion", "grupo_edad"),
  c("grupo_edad", "dosis"),
  c("sexo", "dosis"),
  c("ola_covid", "sexo"),
  c("ola_covid", "region"),
  c("ola_covid", "macrorregion"),
  c("ola_covid", "dosis"),
  c("ola_covid", "grupo_edad", "dosis"),
  c("departamento", "dosis"),
  c("departamento", "sexo", "dosis"),
  c("departamento", "grupo_edad", "dosis"),
  c("departamento", "sexo", "grupo_edad", "dosis"),
  c("departamento", "ola_covid", "dosis"),
  c("distrito", "dosis"),
  c("distrito", "sexo", "dosis"),
  c("distrito", "grupo_edad", "dosis"),
  c("distrito", "cuartil_pobreza", "dosis"),
  c("grupo_edad", "cuartil_pobreza", "dosis"),
  c("sexo", "cuartil_pobreza", "dosis"),
  c("ola_covid", "cuartil_pobreza", "dosis"),
  c("cuartil_pobreza", "dosis"),
  c("cuartil_nbi", "dosis"),
  c("region", "dosis"),
  c("region", "grupo_edad", "dosis"),
  c("macrorregion", "dosis"),
  c("macrorregion", "grupo_edad", "dosis")
)

# Calcular la tasa de vacunacion
results_vacunas <- index_calculator(
  data = vacunas_18, 
  denom = poblacion %>% mutate(overall = "overall"), 
  levels = levels, 
  type = "month"
)

# Calcular la tasa de vacunacion
results__per_month_vacunas <- index_calculator(
  data = vacunas_18, 
  denom = poblacion %>% mutate(overall = "overall"), 
  levels = levels, 
  type = "per_month"
)
```

```{r}
vacunas_18_descr_tbl <- vacunas_18 %>% 
  select(edad, overall, grupo_edad, sexo, ola_covid, region, 
         macrorregion, cuartil_pobreza, cuartil_nbi,
         dosis) %>% 
  collect() 

labelled::var_label(vacunas_18_descr_tbl) <- list(
  overall = "Overall",
  grupo_edad = "Age Group",
  sexo = "Sex",
  ola_covid = "COVID-19 wave",
  region = "Natural Regions",
  macrorregion = "Macroregions",
  cuartil_pobreza = "Poverty level",
  cuartil_nbi = "Unsatisfied basic needs level",
  dosis = "Doses"
)
```

```{r}
#| label: tbl-summ-vacunas-g
#| tbl-cap: !expr paste0("Descriptivos de vacunación menores de 18 años ", range_vacu$fecha)

tbl_summ_vacunas_g <- vacunas_18_descr_tbl %>% 
  select(-edad) %>% 
  gtsummary::tbl_summary(
    type = list(
      overall ~ "dichotomous"
    ),
    value = list(overall = "overall"),
    missing_text = "No data"
  ) %>% 
  gtsummary::modify_header(label = "**Variables**") %>%
  gtsummary::bold_labels()


tbl_summ_vacunas_c2 <- results_vacunas$overall$index %>% 
  mutate(label = "Overall") %>% 
  rename(variable = overall) %>% 
  bind_rows(
    results_vacunas$grupo_edad$index %>% 
      mutate(variable = "grupo_edad") %>% 
      rename(label = grupo_edad),
    results_vacunas$sexo$index %>% 
      mutate(variable = "sexo") %>% 
      rename(label = sexo),
    results_vacunas$ola_covid$index %>% 
      mutate(variable = "ola_covid") %>% 
      rename(label = ola_covid),
    results_vacunas$region$index %>% 
      mutate(variable = "region") %>% 
      rename(label = region),
    results_vacunas$macrorregion$index %>% 
      mutate(variable = "macrorregion") %>% 
      rename(label = macrorregion),
    results_vacunas$cuartil_pobreza$index %>% 
      mutate(variable = "cuartil_pobreza") %>% 
      rename(label = cuartil_pobreza),
    results_vacunas$cuartil_nbi$index %>% 
      mutate(variable = "cuartil_nbi") %>% 
      rename(label = cuartil_nbi),
    results_vacunas$dosis$index %>% 
      mutate(dosis = as.character(dosis),
             variable = "dosis") %>% 
      rename(label = dosis)
  ) %>% 
  drop_na() %>% 
  mutate(
    rate = formattable::digits(rate, 3)
  )

tbl_summ_vacunas_final <- tbl_summ_vacunas_g %>%
  gtsummary::modify_table_body(
    ~ .x %>% 
      left_join(tbl_summ_vacunas_c2) %>% 
      mutate(
        label = case_when(
          variable %in% c("cuartil_pobreza", "cuartil_nbi") &
            label == "1" ~ "1 (0-19.9%)",
          variable %in% c("cuartil_pobreza", "cuartil_nbi") & 
            label == "2" ~ "2 (20-39.9%)",
          variable %in% c("cuartil_pobreza", "cuartil_nbi") & 
            label == "3" ~ "3 (40-59.9%)",
          variable %in% c("cuartil_pobreza", "cuartil_nbi") &
            label == "4" ~ "4 (60-100%)",
          .default = label
        )
      )
  ) %>% 
  gtsummary::modify_header(
    rate ~ "**Vaccine dose rate per 1,000 inhabitants**"
  ) %>%
  gtsummary::modify_fmt_fun(
    rate ~ scales::number_format(accuracy = 0.001)
  )

tbl_summ_vacunas_final
```

```{r}
tbl_summ_vacunas_final %>% 
  gtsummary::as_hux_xlsx("02_output/tables/vacunas_descriptivos.xlsx")

tbl_summ_vacunas_final %>%
  gtsummary::as_flex_table() %>%
  flextable::save_as_docx(path = "02_output/tables/vacunas_descriptivos.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/vacunas_descriptivos.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/vacunas_descriptivos.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::



<!-- Cálculo de Comparaciones Múltiples -->

```{r}
comparing_levels <- list(
  c("grupo_edad"),
  c("sexo"),
  c("ola_covid"),
  c("cuartil_pobreza"),
  c("cuartil_nbi"),
  c("region"),
  c("macrorregion")
)

vacunas_comparisons <- compare_levels_per_month(
  results_per_month = results__per_month_vacunas,
  levels = comparing_levels
)

# Extraer las tablas `cld` como tibbles, renombrando y añadiendo una columna "level"
vacunas_comparisons_cld_pre <- lapply(names(vacunas_comparisons), function(level) {
  cld <- vacunas_comparisons[[level]]$cld
  
  if (!is.null(cld)) {
    cld <- cld |> 
      as_tibble() |> 
      rename(group = 1) |>  # Renombra la primera columna a 'group'
      mutate(level = level) |> # Añade una columna 'level' con el nombre de la lista
      relocate(level) 
  }
  
  cld
}) |> 
  bind_rows()

vacunas_comparisons_cld_pre <- 
  vacunas_comparisons_cld_pre |> 
  mutate(
    across(
      where(is.numeric),
      ~ formattable::digits(.x, 3)
    )
  ) 
```

```{r}
#| label: tbl-compar-vacunas
#| tbl-cap: !expr paste0("Comparaciones múltiples de vacunas", range_pos_18_0723$fecha)
#| 
vacunas_comparisons_cld <- vacunas_comparisons_cld_pre |> 
  rename_with(
    .fn = ~ paste0("**", str_to_title(.x), "**")
  ) %>% 
  gt::gt() %>% 
  gt::cols_label(
    .fn = gt::md,
    .process_units = TRUE
  )

vacunas_comparisons_cld
```

```{r}
vacunas_comparisons_cld_pre %>% 
  rename_with(
    .fn = str_to_title
  ) %>% 
  huxtable::as_hux() %>% 
  huxtable::set_bold(1, 1:8) %>% 
  huxtable::quick_xlsx(file = "02_output/tables/vacunas_comparacion_multiple.xlsx",
                       open = FALSE)

vacunas_comparisons_cld %>%
  gt::gtsave("02_output/tables/vacunas_comparacion_multiple.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/vacunas_comparacion_multiple.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/vacunas_comparacion_multiple.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

<!-- Comparaciones de vacunacion con respecto a olas -->


```{r}
#| label: tbl-summ-vacunas-olas
#| tbl-cap: !expr paste0("Descriptivos de vacunación menores de 18 años por olas ", range_vacu$fecha)

tbl_summ_vacunas_olas <- vacunas_18_descr_tbl %>% 
  select(-edad) %>% 
  gtsummary::tbl_summary(
    by = ola_covid,
    type = list(
      overall ~ "dichotomous"
    ),
    value = list(overall = "overall"),
    missing_text = "No data"
  ) %>% 
  gtsummary::modify_header(label = "**Variables**") %>%
  gtsummary::bold_labels()

ola_covid_vacunas_gtsumm <- results_vacunas$ola_covid$index %>% 
  rename(stat = ola_covid) %>% 
  mutate(label = "Overall",
         variable = "overall") %>% 
  bind_rows(
    results_vacunas$ola_covid_grupo_edad$index %>%
      mutate(variable = "grupo_edad") %>% 
      rename(label = grupo_edad,
             stat = ola_covid),
    results_vacunas$ola_covid_sexo$index %>% 
      mutate(variable = "sexo") %>% 
      rename(label = sexo,
             stat = ola_covid),
    results_vacunas$ola_covid_region$index %>% 
      mutate(variable = "region") %>% 
      rename(label = region,
             stat = ola_covid),
    results_vacunas$ola_covid_macrorregion$index %>% 
      mutate(variable = "macrorregion") %>%
      rename(label = macrorregion,
             stat = ola_covid),
    results_vacunas$ola_covid_cuartil_pobreza$index %>% 
      mutate(variable = "cuartil_pobreza") %>% 
      rename(label = cuartil_pobreza,
             stat = ola_covid),
    results_vacunas$ola_covid_cuartil_nbi$index %>% 
      mutate(variable = "cuartil_nbi") %>% 
      rename(label = cuartil_nbi,
             stat = ola_covid),
    results_vacunas$ola_covid_dosis$index %>% 
      mutate(dosis = as.character(dosis),
             variable = "dosis") %>% 
      rename(label = dosis,
             stat = ola_covid)
  ) %>% 
  mutate(
    rate = formattable::digits(rate, 3),
    stat = as.character(stat),
    stat = case_match(
      stat,
      "1 (Mar – Oct 2020)" ~ "stat_1",
      "2 (Nov 2020 – Oct 2021)" ~ "stat_2",
      "3 (Nov 2021 – Apr 2022)" ~ "stat_3",
      "4 (May 2022 – Oct 2022)" ~ "stat_4",
      "5 (Nov 2022 – Jun 2023)" ~ "stat_5"
    )
  ) 

tbl_summ_vacunas_olas$table_body <- tbl_summ_vacunas_olas$table_body %>%
  mutate(
    across(
      gtsummary::all_stat_cols(stat_0 = FALSE),
      ~ str_remove_all(., "\\s*\\(.*?\\)\\s*")
    )
  ) %>% 
  pivot_longer(
    cols = stat_1:stat_5,
    names_to = "stat",
    values_to = "freq"
  ) %>% 
  left_join(
    ola_covid_vacunas_gtsumm
  ) %>% 
  mutate(
    value = str_c(freq, ifelse(is.na(rate), "", str_c(" (", rate, ")")))
  ) %>% 
  select(-c(freq, rate)) %>% 
  pivot_wider(
    names_from = stat,
    values_from = value
  )


tbl_summ_vacunas_olas <- tbl_summ_vacunas_olas %>% 
  gtsummary::modify_footnote(
    gtsummary::all_stat_cols() ~ "n (Vaccine dose rate per 1,000 inhabitants)"
  )

tbl_summ_vacunas_olas
```

```{r}
tbl_summ_vacunas_olas %>% 
  gtsummary::as_hux_xlsx("02_output/tables/vacunas_olas_descriptivos.xlsx")

tbl_summ_vacunas_olas %>%
  gtsummary::as_flex_table() %>%
  flextable::save_as_docx(path = "02_output/tables/vacunas_olas_descriptivos.docx")
```

::: {style="display: flex; justify-content: space-around;"}
```{r}
downloadthis::download_file(
  path = "02_output/tables/fallecidos_olas_descriptivos.xlsx",
  button_label = "Descargar tabla en Excel",
  button_type = "success",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```

```{r}
downloadthis::download_file(
  path = "02_output/tables/fallecidos_olas_descriptivos.docx",
  button_label = "Descargar tabla en Word",
  button_type = "primary",
  has_icon = TRUE,
  icon = "fa fa-save"
)
```
:::

## Gráfico de vacunación (serie de tiempo) - Menores de 18 años

Se restringió la visualización de dosis de 1 a 5.

```{r}
#| fig-cap: !expr paste0("Vacunación COVID Población menor de edad de acuerdo a la dosis y ola ", range_vacu$fecha)
#| label: fig-vac-dosis-olas-area
#| fig-width: 10

summ_vacunas_dosis_fecha <- vacunas_18 %>% 
  to_arrow() %>%
  mutate(
    fecha_round = ceiling_date(fecha_vacunacion, unit = "month"),
    .after = fecha_vacunacion
  ) %>%
  # filter(dosis %in% 1:5) %>% 
  count(dosis, fecha_round) %>% 
  collect()

summ_vacunas_dosis_fecha %>% 
  ggplot(
    aes(
      x = fecha_round,
      y = n,
      fill = as.factor(dosis)
    )
  ) +
  geom_area(position = 'stack') +
  scale_fill_manual(
    values = c("#2C9BB4", "#0B4E60", "#8C3A33")
  ) +
  # scale_x_date(
  #   limits = c(ymd("2020-03-01"), ymd("2023-10-31")),
  #   breaks = scales::date_breaks("2 month"), # Cambiar a "6 months" para menos frecuencia
  #   labels = scales::date_format("%b %Y") # Formato de fecha, %b es la abreviatura de mes, %Y es el año
  # ) + 
  scale_x_date(
    limits = c(ymd("2020-03-01"), ymd("2023-07-01")),
    breaks = seq(from = ceiling_date(ymd("2020-03-01"), "month"),
                 to = floor_date(ymd("2023-07-01"), "month"),
                 by = "2 months"), 
    # date_breaks = "2 month",
    date_labels = "%b %Y"
  ) +
  scale_y_continuous(
    labels = scales::number_format()
  ) +
  labs(
    y = NULL,
    x = NULL,
    title = "Vaccination COVID-19 under 18 years old",
    fill = "Doses"
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2020-03-01"), 
    xmax = as.Date("2020-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2020-07-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "First wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2020-11-01"), 
    xmax = as.Date("2021-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2021-05-10"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Second wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2021-11-01"), 
    xmax = as.Date("2022-04-30"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
    annotate(
    geom = "text", 
    x = as.Date("2022-02-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Third \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2022-05-01"), 
    xmax = as.Date("2022-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-08-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Fourth \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2022-11-01"), 
    xmax = dmy(range_vacu$max), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2023-03-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Fifth wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  theme_classic() +
  theme(
    axis.text.x = element_text(
      angle = 75,
      vjust = 0.5
    ),
    axis.ticks.x = element_line()
  ) 
```

## Gráfico de vacunación (serie de tiempo) - Grupo Etáreo

```{r}
#| fig-cap: !expr paste0("Vacunación COVID Población menor de edad de acuerdo al grupo etáreo y ola ", range_vacu$fecha)
#| label: fig-vac-edad-olas-area
#| fig-width: 10

summ_vacunas_edad_fecha <- vacunas_18 %>% 
  to_arrow() %>%
  mutate(
    fecha_round = ceiling_date(fecha_vacunacion, unit = "month"),
    .after = fecha_vacunacion
  ) %>%
  count(grupo_edad, fecha_round) %>% 
  collect() 

summ_vacunas_edad_fecha %>% 
  ggplot(
    aes(
      x = fecha_round,
      y = n,
      fill = grupo_edad
    )
  ) +
  geom_area(position = 'stack') +
  scale_fill_manual(
    values = c("#2C9BB4", "#0B4E60", "#DAE2BC", "#8C3A33")
  ) +
  # scale_x_date(
  #   limits = c(ymd("2020-03-01"), ymd("2023-10-31")),
  #   breaks = scales::date_breaks("2 month"), # Cambiar a "6 months" para menos frecuencia
  #   labels = scales::date_format("%b %Y") # Formato de fecha, %b es la abreviatura de mes, %Y es el año
  # ) + 
  scale_x_date(
    limits = c(ymd("2020-03-01"), ymd("2023-07-01")),
    breaks = seq(from = ceiling_date(ymd("2020-03-01"), "month"),
                 to = floor_date(ymd("2023-07-01"), "month"),
                 by = "2 months"), 
    # date_breaks = "2 month",
    date_labels = "%b %Y"
  ) +
  scale_y_continuous(
    labels = scales::number_format()
  ) +
  labs(
    y = NULL,
    x = NULL,
    title = "Vaccination COVID-19 under 18 years old",
    fill = "Age Group"
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2020-03-01"), 
    xmax = as.Date("2020-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2020-07-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "First wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2020-11-01"), 
    xmax = as.Date("2021-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2021-05-10"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Second wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2021-11-01"), 
    xmax = as.Date("2022-04-30"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
    annotate(
    geom = "text", 
    x = as.Date("2022-02-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Third \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) + 
  annotate(
    geom = "rect",
    xmin = as.Date("2022-05-01"), 
    xmax = as.Date("2022-10-31"), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2022-08-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Fourth \nwave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 1.25 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  annotate(
    geom = "rect",
    xmin = as.Date("2022-11-01"), 
    xmax = dmy(range_vacu$max), 
    ymin = 0, # El rectángulo se extiende hasta la parte inferior del gráfico
    ymax = Inf, # El rectángulo se extiende hasta la parte superior del gráfico
    color = "black", # Color de la línea
    fill = NA, # Sin relleno
    linetype = "dashed", # Tipo de línea discontinua
    linewidth = 0.5 # Grosor de la línea
  ) + 
  annotate(
    geom = "text", 
    x = as.Date("2023-03-01"), # Punto medio aproximado de la Ola 2
    y = Inf, # Posición vertical del texto
    label = "Fifth wave", 
    color = "black",
    size = 4, # Tamaño del texto
    vjust = 2 # Ajuste vertical para colocar el texto encima del rectángulo
  ) +
  theme_classic() +
  theme(
    axis.text.x = element_text(
      angle = 75,
      vjust = 0.5
    ),
    axis.ticks.x = element_line()
  ) 
```

## Pirámide frecuencia - Edades

```{r}
library(patchwork)

create_age_pyramids <- function(data, age, by, proportional = FALSE, 
                                scale_adjust = 7, angle.axis = NULL) {
  plots <- list()
  
  for (dosis in 1:3) {
    # Filtra por dosis y crea el gráfico
    p <- data %>%
      filter(dosis == dosis) %>%
      age_pyramid(age_group = !!rlang::enquo(age), 
                  split_by = !!rlang::enquo(by),
                  proportional = proportional) +
      scale_fill_brewer(palette = "Dark2") +
      labs(
        subtitle = paste("Dosis", dosis),
        x = NULL,
        fill = "Sexo",
        y = NULL
      ) +
      theme_classic() +
      theme(
        axis.text.x = element_text(
          angle = angle.axis,
          vjust = 0.5,
          hjust = 0.5
        )
      )
    
    if (!proportional) {
      p <- p + 
        scale_y_continuous(
        labels = ~ scales::number(abs(.x)),
        n.breaks = scale_adjust
      ) 
    }
    
    if (dosis %in% 1:2) {
      p <- p +
        theme(
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank()
        )
    }
    
    plots[[dosis]] <- p
  }
  
  # Combina todos los gráficos en uno solo usando patchwork
  combined_plot <- reduce(plots, `+`) +
    plot_layout(guides = 'collect') &
    labs(x = "Edad", y = "Frecuencia") +
    theme(legend.position = 'bottom')

  return(combined_plot)
}
```

```{r}
#| dpi: 300
#| fig-width: 10
#| fig-height: 6
#| warning: false

vacunas_18_descr_tbl %>% 
  mutate(edad = factor(edad)) %>%
  create_age_pyramids(age = edad,
                      by = sexo)
```

```{r}
#| dpi: 300
#| fig-width: 10
#| fig-height: 6
#| warning: false

vacunas_18_descr_tbl %>% 
  create_age_pyramids(age = grupo_edad, 
                      by = sexo, 
                      scale_adjust = 6,
                      angle.axis = 20)
```

## Pirámide proporcional - Edades

```{r}
#| fig-width: 10
#| fig-height: 6
#| warning: false

vacunas_18_descr_tbl %>% 
  mutate(edad = factor(edad)) %>%
  create_age_pyramids(age = edad,
                      by = sexo,
                      proportional = TRUE)
```

```{r}
#| fig-width: 10
#| fig-height: 6
#| warning: false

vacunas_18_descr_tbl %>% 
  create_age_pyramids(age = grupo_edad,
                      by = sexo,
                      proportional = TRUE,
                      angle.axis = 45)
```

## Departamento, edad y Dosis

<!-- Agregar una columna más donde sea "Total por edad" dividido por dosis -->

<!-- Agregar en leyenda que plomo es sin registro -->

```{r}
vac_dep_edad_dosis_1k_mother <- crossing(
  departamento = unique(results_vacunas$departamento_grupo_edad_dosis$index$departamento),
  grupo_edad = unique(results_vacunas$departamento_grupo_edad_dosis$index$grupo_edad),
  dosis = 1:3
)

vac_dep_edad_dosis_1k <- vac_dep_edad_dosis_1k_mother %>% 
  left_join(results_vacunas$departamento_grupo_edad_dosis$index)

vac_dep_dosis_1k_mother <- crossing(
  departamento = unique(results_vacunas$departamento_dosis$index$departamento),
  dosis = 1:3
)

vac_dep_dosis_1k <- vac_dep_dosis_1k_mother %>% 
  left_join(results_vacunas$departamento_dosis$index)


vac_dep_edad_dosis_1k_sf <- peru_dep_sf %>% 
  rename(departamento = dep) %>% 
  inner_join(
    vac_dep_edad_dosis_1k %>% 
      bind_rows(
        vac_dep_dosis_1k %>% 
          mutate(grupo_edad = "Overall")
      ) %>% 
      mutate(
        grupo_edad = fct_relevel(
          grupo_edad,
          "Under 1 year",
          "1 - 5 years",
          "6 - 11 years",
          "12 - 17 years",
          "Overall"
        ),
        dosis = factor(
          dosis,
          labels = c("1st Dose", "2nd Dose",
                     "3rd Dose")
        )
      )
  ) %>% 
  ungroup()
```

```{r}
#| fig-cap: !expr paste0("Tasa de Vacunación COVID-19 de acuerdo al departamento, edad y dosis ", range_vacu$fecha)
#| label: fig-vac-dep-edad-dosis-map

vac_dep_edad_dosis_1k_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap(
        "Vaccination rate per 1000 inhabitants",
        20
      )
    )
  ) +
  scale_fill_distiller(palette = "BuGn",
                       direction = 1,
                       na.value = "grey70") + 
  facet_grid(
    vars(dosis),
    vars(grupo_edad),
    labeller = labeller(grupo_edad = label_wrap_gen(10))
  ) +
  theme_void()
```

## Departamento y edad

<!-- Hacerlo en base a que tengan vacunación completa como mínimo (2 dosis) -->

```{r}
vac_dep_edad_1k_mother <- crossing(
  departamento = unique(results_vacunas$departamento_grupo_edad$index$departamento),
  grupo_edad = unique(results_vacunas$departamento_grupo_edad$index$grupo_edad)
)

vac_dep_edad_1k <- vac_dep_edad_1k_mother %>% 
  left_join(results_vacunas$departamento_grupo_edad$index)

vac_dep_edad_1k_sf <- peru_dep_sf %>% 
  rename(departamento = dep) %>% 
  inner_join(
    vac_dep_edad_1k
  ) %>% 
  ungroup()
```

```{r}
#| fig-cap: !expr paste0("Tasa de Vacunación COVID-19 de acuerdo al departamento y edad ", range_vacu$fecha)
#| label: fig-vac-dep-edad-map

vac_dep_edad_1k_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap(
        "Vaccination rate per 1000 inhabitants",
        20
      )
    )
  ) +
  scale_fill_distiller(palette = "BuGn",
                       direction = 1,
                       na.value = "grey70") + 
  facet_wrap(
    vars(grupo_edad),
    labeller = labeller(grupo_edad = label_wrap_gen(10)),
    nrow = 1
  ) +
  theme_void()
```

## Departamento y dosis

```{r}
vac_dep_dosis_1k_sf <- peru_dep_sf %>% 
  rename(departamento = dep) %>% 
  inner_join(
    vac_dep_dosis_1k %>% 
      mutate(
        dosis = factor(
          dosis,
          labels = c("1st Dose", "2nd Dose",
                     "3rd Dose")
        )
      )
  ) %>% 
  ungroup()
```

```{r}
#| fig-cap: !expr paste0("Tasa de Vacunación COVID-19 de acuerdo al departamento y dosis ", range_vacu$fecha)
#| label: fig-vac-dep-dosis-map

vac_dep_dosis_1k_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap(
        "Vaccination rate per 1000 inhabitants",
        20
      )
    )
  ) +
  scale_fill_distiller(palette = "BuGn",
                       direction = 1,
                       na.value = "grey70") + 
  facet_wrap(
    vars(dosis)
  ) +
  theme_void()
```

## Region, edad vs Dosis

<!-- Agregar en leyenda que plomo es sin registro -->

```{r}
vac_reg_edad_dosis_1k_mother <- crossing(
  region = unique(results_vacunas$region_grupo_edad_dosis$index$region),
  grupo_edad = unique(results_vacunas$region_grupo_edad_dosis$index$grupo_edad),
  dosis = 1:3
)

vac_reg_edad_dosis_1k <- vac_reg_edad_dosis_1k_mother %>% 
  left_join(results_vacunas$region_grupo_edad_dosis$index)


vac_reg_dosis_1k_mother <- crossing(
  region = unique(results_vacunas$region_dosis$index$region),
  dosis = 1:3
)

vac_reg_dosis_1k <- vac_reg_dosis_1k_mother %>% 
  left_join(results_vacunas$region_dosis$index)


vac_reg_edad_dosis_1k_sf <- peru_reg_sf %>% 
  inner_join(
    vac_reg_edad_dosis_1k %>% 
      bind_rows(
        vac_reg_dosis_1k %>% 
          mutate(grupo_edad = "Overall")
      ) %>% 
      mutate(
        grupo_edad = fct_relevel(
          grupo_edad,
          "Under 1 year",
          "1 - 5 years",
          "6 - 11 years",
          "12 - 17 years",
          "Overall"
        ),
        dosis = factor(
          dosis,
          labels = c("1st Dose", "2nd Dose",
                     "3rd Dose")
        )
      )
  ) %>% 
  ungroup()
```

```{r}
#| fig-cap: !expr paste0("Tasa de Vacunación COVID-19 de acuerdo a la región, edad y dosis ", range_vacu$fecha)
#| label: fig-vac-reg-edad-dosis-map

vac_reg_edad_dosis_1k_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap(
        "Vaccination rate per 1000 inhabitants",
        20
      )
    )
  ) +
  scale_fill_distiller(palette = "BuGn",
                       direction = 1,
                       na.value = "grey70") + 
  facet_grid(
    vars(dosis),
    vars(grupo_edad),
    labeller = labeller(grupo_edad = label_wrap_gen(10))
  ) +
  theme_void()
```

## Mapa de Tasa de Vacunación - 1000 habitantes Macrorregión, edad vs Dosis

<!-- Agregar una columna más donde sea "Total por edad" dividido por dosis -->

<!-- Agregar en leyenda que plomo es sin registro -->

```{r}
pob_anio_macro_edad <- poblacion %>% 
  group_by(anio, macrorregion, grupo_edad) %>% 
  summarise(
    pob = sum(cantidad)
  ) %>% 
  window_order(anio, macrorregion, grupo_edad) %>% 
  ungroup() 

vac_mes_macro_edad_dosis <- vacunas_18 %>% 
  filter(dosis %in% 1:3) %>% 
  count(
    anio, mes, macrorregion, grupo_edad, dosis
  )

vac_mes_macro_edad_dosis_adic <- vac_mes_macro_edad_dosis %>% 
  left_join(
    pob_anio_macro_edad
  ) 

vac_macro_edad_dosis_1k <- vac_mes_macro_edad_dosis_adic %>% 
  mutate(
    vacunacion_mes = n/pob
  ) %>% 
  group_by(across(macrorregion:dosis)) %>% 
  summarise(
    vacunacion_mes = mean(vacunacion_mes, na.rm = TRUE)*1000
  ) %>% 
  ungroup() %>% 
  collect() 

vac_macro_edad_dosis_1k_mother <- crossing(
  macrorregion = unique(vac_macro_edad_dosis_1k$macrorregion),
  grupo_edad = unique(vac_macro_edad_dosis_1k$grupo_edad),
  dosis = 1:3
)


vac_macro_edad_dosis_1k <- vac_macro_edad_dosis_1k_mother %>% 
  left_join(vac_macro_edad_dosis_1k)

vac_macro_edad_dosis_1k_sf <- peru_macro_sf %>% 
  inner_join(
    vac_macro_edad_dosis_1k
  ) %>% 
  ungroup()
```

```{r}
vac_macro_edad_dosis_1k_mother <- crossing(
  macrorregion = unique(results_vacunas$macrorregion_grupo_edad_dosis$index$macrorregion),
  grupo_edad = unique(results_vacunas$macrorregion_grupo_edad_dosis$index$grupo_edad),
  dosis = 1:3
)

vac_macro_edad_dosis_1k <- vac_macro_edad_dosis_1k_mother %>% 
  left_join(results_vacunas$macrorregion_grupo_edad_dosis$index)


vac_macro_dosis_1k_mother <- crossing(
  macrorregion = unique(results_vacunas$macrorregion_dosis$index$macrorregion),
  dosis = 1:3
)

vac_macro_dosis_1k <- vac_macro_dosis_1k_mother %>% 
  left_join(results_vacunas$macrorregion_dosis$index)


vac_macro_edad_dosis_1k_sf <- peru_macro_sf %>% 
  inner_join(
    vac_macro_edad_dosis_1k %>% 
      bind_rows(
        vac_macro_dosis_1k %>% 
          mutate(grupo_edad = "Overall")
      ) %>% 
      mutate(
        grupo_edad = fct_relevel(
          grupo_edad,
          "Under 1 year",
          "1 - 5 years",
          "6 - 11 years",
          "12 - 17 years",
          "Overall"
        ),
        dosis = factor(
          dosis,
          labels = c("1st Dose", "2nd Dose",
                     "3rd Dose")
        )
      )
  ) %>% 
  ungroup()
```

```{r}
#| fig-cap: !expr paste0("Tasa de Vacunación COVID-19 de acuerdo a la macrorregión, edad y dosis ", range_vacu$fecha)
#| label: fig-vac-macro-edad-dosis-map

vac_macro_edad_dosis_1k_sf %>%
  ggplot() +
  geom_sf(aes(fill = rate)) +
  guides(
    fill = guide_legend(
      str_wrap(
        "Vaccination rate per 1000 inhabitants",
        20
      )
    )
  ) +
  scale_fill_distiller(palette = "BuGn",
                       direction = 1,
                       na.value = "grey70") + 
  facet_grid(
    vars(dosis),
    vars(grupo_edad),
    labeller = labeller(grupo_edad = label_wrap_gen(10))
  ) +
  theme_void()
```

# Gráficos Unidos (Antes ID Week)

```{r}
library(patchwork)
```

## Area Plot: Casos, muertes, hospitalizados

```{r}
#| fig-width: 14
#| fig-height: 9.5
#| fig-dpi: 400

plots_area_cases <- (positive_cases_area + hospi_age_area_plot) /
(deaths_cases_by_age_area + deaths_cases_area) +
  plot_layout(guides = 'collect') +
    # labs(x = "Edad", y = "Frecuencia") +
    theme(legend.position = 'bottom') +
  plot_annotation(tag_levels = 'A')

plots_area_cases
```


```{r}
#| eval: false
ggsave(
  "plots_area_cases.png",
  plot = plots_area_cases,
  dpi = 400,
  bg = "white",
  width = 14,
  height = 9.5
)
```

## Area Plot: MMI, HAR, MMM y L

```{r}
#| fig-width: 14
#| fig-height: 9.5
#| fig-dpi: 400

plots_area_index <- (prevalence_by_age_area + mortality_cases_by_age_area) /
(lethality_cases_by_age_area + hospi_rate_by_age_area) +
  plot_layout(guides = 'collect') +
  plot_annotation(tag_levels = 'A')

plots_area_index
```

```{r}
#| eval: false
ggsave(
  "plots_area_index.png",
  plot = plots_area_index,
  dpi = 400,
  bg = "white",
  width = 14,
  height = 9.5
)
```

## Area Plot x Pobreza: MMI, HAR, MMM y L


```{r}
#| fig-width: 14
#| fig-height: 9.5
#| fig-dpi: 400

plots_area_pobre_index <- (prev_pobr_area_plot + mort_pobr_area_plot) /
(letal_pobr_area_plot + hospi_pobr_area_plot) +
  plot_layout(guides = 'collect') +
  plot_annotation(tag_levels = 'A')

plots_area_pobre_index
```

```{r}
#| eval: false
ggsave(
  "plots_area_pobre_index.png",
  plot = plots_area_pobre_index,
  dpi = 400,
  bg = "white",
  width = 14,
  height = 9.5
)
```

## Mapa: MMI, HAR, MMM y L

```{r}
#| fig-width: 8.5
#| fig-height: 5
#| fig-dpi: 400

plots_map_index <- (prev_macro_map + mort_macro_map) /
(let_macro_map + hospi_macro_map) +
  plot_annotation(tag_levels = 'A')

plots_map_index
```

```{r}
#| eval: false

ggsave(
  "plots_map_index.png",
  plot = plots_map_index,
  dpi = 400,
  bg = "white",
  width = 8.5,
  height = 8
)
```

## Mapa x Dep x Region x Macro: MMI, HAR, MMM

```{r}
map_prev_joined <- wrap_plots(
  prev_reg_1k_map,
  prev_macro_map,
  prev_dep_1k_map,
  nrow = 1
)

map_prev_joined[[2]] <- map_prev_joined[[2]] + 
  plot_layout(tag_level = 'new')
map_prev_joined[[3]] <- map_prev_joined[[3]] + 
  plot_layout(tag_level = 'new')

map_hospi_joined <- wrap_plots(
  tasa_hospi_reg_map,
  hospi_macro_map,
  tasa_hospi_dep_map,
  nrow = 1
)

map_hospi_joined[[2]] <- map_hospi_joined[[2]] + 
  plot_layout(tag_level = 'new')
map_hospi_joined[[3]] <- map_hospi_joined[[3]] + 
  plot_layout(tag_level = 'new')

map_mort_joined <- wrap_plots(
  mortalidad_reg_100k_map,
  mort_macro_map,
  fallecidos_dep_100k_map,
  nrow = 1
)

map_mort_joined[[2]] <- map_mort_joined[[2]] + 
  plot_layout(tag_level = 'new')
map_mort_joined[[3]] <- map_mort_joined[[3]] + 
  plot_layout(tag_level = 'new')
```

```{r}
#| fig-width: 10
#| fig-height: 9
#| fig-dpi: 400

map_rates_joined <- wrap_plots(
  map_prev_joined,
  map_hospi_joined,
  map_mort_joined,
  nrow = 3
) +
  plot_annotation(tag_levels = 'A')

map_rates_joined
```


```{r}
#| eval: false

ggsave(
  "map_rates_joined.png",
  plot = map_rates_joined,
  dpi = 400,
  bg = "white",
  width = 10,
  height = 9
)
```

## Mapa: Por olas

```{r}
#| fig-width: 13
#| fig-height: 9
#| fig-dpi: 400

plots_map_waves <- wrap_plots(
  prev_mes_dep_olas_1k_map +
    theme(
      strip.text.x = element_text(
        size = 12,
        face = "bold",
        margin = margin(t = -1.25)
      )
    ),
  tasa_hospi_dep_olas_map +
    theme(strip.text.x = element_blank()),
  mortalidad_dep_olas_100k_map +
    theme(strip.text.x = element_blank()),
  nrow = 3
) + 
  plot_annotation(tag_levels = 'A')

plots_map_waves
```

```{r}
#| eval: false

ggsave(
  "plots_map_waves.png",
  plot = plots_map_waves,
  dpi = 400,
  bg = "white",
  width = 13,
  height = 9
)
```




```{r}
dbDisconnect(con, shutdown=TRUE)
```
