05_deaths_association.qmd

---
title: "Relación con muertes"
author: "Brian Norman Peña-Calero"
date: "`r Sys.Date()`"
format: 
  html:
    toc: true
    number-sections: true
    anchor-sections: true
    code-line-numbers: true
    code-overflow: scroll
    code-copy: hover
    code-link: true
    code-fold: show
    lang: es
    embed-resources: true
    dpi: 400
lightbox: true
editor_options: 
  chunk_output_type: console
---


# Cargar paquetes

```{r}
library(tidyverse)
library(lme4)
# library(lmerTest)
# library(merDeriv)
library(performance)
library(parameters)
library(glmmTMB)
library(sjPlot)
# library(MASS)
# options(scipen = 9)

# conflicted::conflicts_prefer(dplyr::filter)
# conflicted::conflicts_prefer(dplyr::select)
# conflicted::conflicts_prefer(lmerTest::lmer)

options(parameters_labels = TRUE) 
```

# Importar datos

```{r}
gbmt_fit_final <- readRDS("01_data/processed/gbmt_fit_final.rds")
deaths_by_year <- read_csv(file = "01_data/processed/deaths_by_year.csv")
```

# Unificar datos

```{r}
gbmt_log <- gbmt_fit_final %>% 
  select(ng, gbmt_fit_total) |> 
  unnest(gbmt_fit_total)

gbmt_log_merged <- gbmt_log %>% 
  select(ng, ISO2, pubsalid1, SALID1, 
         YEAR = year, population_imp_norm:group,
         country, City) %>% 
  inner_join(
    deaths_by_year,
    by = join_by(ISO2, SALID1, YEAR),
    relationship = "many-to-many"
  ) %>% 
  mutate(
    deaths_per_100k = (deaths / population_imp) * 100000,
    population_imp_log = log(population_imp),
    country = str_to_title(country)
  )

labelled::var_label(gbmt_log_merged) <- list(
  country = "Country",
  population_imp = "Population",
  bectuareal1ux_imp = "Urban Area",
  deaths_per_100k	 = "Deaths per 100K",
  group = "Clusters",
  GTHGROUP = "Deaths per 100k"
)
```

# Tabla descriptiva

```{r}
library(gtsummary)
```

```{r}
table1_descriptive <- function(data) {  
  if (length(unique(data$group)) == 10) {
    data <- data %>% 
      mutate(
        group = fct_relevel(group, "10",
                            after = Inf)
      )
  }

  tab1 <- data %>% 
    filter(Tier == "Overall") %>% 
    select(country, group) %>% 
    tbl_summary(
      by = group,
      percent = "row"
    ) %>% 
    modify_header(label = "") %>%
    bold_labels()

  tab2 <- data %>% 
    filter(Tier == "Overall") %>% 
    select(group, population_imp, bectuareal1ux_imp) %>% 
    tbl_summary(
      by = group,
      percent = "row"
    ) %>% 
    modify_header(label = "") %>%
    bold_labels() 

  data_tab3 <- data %>% 
    filter(Tier %in% c("Overall", "T1")) %>% 
    mutate(GTHGROUP = fct_relevel(GTHGROUP, "Overall"))

  tab_3_1 <- data_tab3 %>% 
    select(group, GTHGROUP) %>% 
    tbl_summary(
      by = group,
      percent = "row"
    ) %>% 
    bold_labels()

  tab_3_2 <- data_tab3 %>% 
    select(group, GTHGROUP, deaths_per_100k) %>% 
    tbl_strata(
      strata = group,
      .tbl_fun = 
        ~ .x %>% 
        tbl_summary(
          by = GTHGROUP,
          percent = "row"
        ) 
    ) 

  tab_3_2_mod <- tab_3_2$table_body %>%
    pivot_longer(
      cols = starts_with("stat_"),
      names_to = c("stat_type", "group"),
      names_pattern = "stat_([0-9]+)_([0-9]+)"
    ) %>% 
    pivot_wider(
      names_from = group,
      values_from = value,
      names_glue = "stat_{group}"
    ) %>% 
    mutate(
      label = factor(
        stat_type,
        labels = c(
          "Overall",
          "Communicable, maternal, perinatal and nutritional conditions",
          "Injuries",
          "Noncommunicable diseases"
        )
      )
    ) %>% 
    select(label, all_stat_cols(stat_0 = FALSE))

  tab_3_1$table_body <- tab_3_1$table_body %>% 
    select(-all_stat_cols(stat_0 = FALSE)) %>% 
    left_join(
      tab_3_2_mod,
      by = join_by(label)
    )
    
  if (length(unique(data$group)) != 1) {
    tab1 <- tab1 %>% 
      add_overall(last = TRUE)
    tab2 <- tab2 %>% 
      add_overall(last = TRUE)

    tab_3_3 <- data_tab3 %>% 
      select(GTHGROUP, deaths_per_100k) %>% 
      tbl_summary(
        by = GTHGROUP
      ) 

    tab_3_3_mod <- tab_3_3$table_body %>% 
      pivot_longer(
        cols = starts_with("stat_"),
        names_to = c("group"),
        names_pattern = "stat_([0-9]+)"
      ) %>% 
      mutate(
        label = factor(
          group,
          labels = c(
            "Overall",
            "Communicable, maternal, perinatal and nutritional conditions",
            "Injuries",
            "Noncommunicable diseases"
          )
        )
      ) %>% 
      select(label, stat_0 = value)

    tab_3_1$table_body <- tab_3_1$table_body %>% 
      left_join(
        tab_3_3_mod,
        by = join_by(label)
      )

    tab_3_1$table_styling$header <- tab_3_1$table_styling$header %>% 
      bind_rows(
        tab_3_3$table_styling$header %>% 
          filter(modify_stat_level == "Overall") %>% 
          mutate(column = "stat_0")
      )
  }

  tbl_stack(
    list(tab1,
         tab2,
         tab_3_1),
    quiet = TRUE
  )
}
```

```{r}
gbmt_log_tab1 <- gbmt_log_merged %>% 
  group_nest(ng) %>%
  mutate(
    tab1 = map(
      data,
      table1_descriptive
    )
  )

walk2(
  gbmt_log_tab1$tab1,
  gbmt_log_tab1$ng,
  ~ .x %>% 
    as_flex_table() %>%
    flextable::save_as_docx(path = paste0("./02_output/tables/descriptive1",
                                          "_group_", .y, ".docx"))
)
```


# Gráfico inicial de promedio de muertes

```{r}
gbmt_log_merged %>% 
  group_by(Clusters, group, YEAR) %>% 
  summarise(
    deaths_per_100k = mean(deaths_per_100k)
  ) %>% 
  ungroup() %>% 
  filter(Clusters == 4) %>% 
  ggplot(aes(x = YEAR, y = deaths_per_100k, color = as.factor(group))) +
  geom_line(aes(group = group)) +
  geom_point() +
  
  labs(title = "Promedio de Muertes por 100,000 habitantes por Grupo a lo largo del Tiempo",
       x = "Año", y = "Muertes por 100,000 habitantes", color = "Grupo") +
  # facet_wrap(vars(Clusters)) +
  theme_minimal()
```


```{r} 
#| eval: false
ggsave(
  filename = "02_output/plots/promedio_muertes_group_log.png",
  dpi = 300,
  bg = "white"
)
```


# Poisson - Modelo lineal mixto con Efectos Aleatorios 

## Con efectos aleatorios 

### Clusters: 3

```{r}
m1_poisson_w_re_3_glmmTMB <- gbmt_log_merged %>%
  filter(Clusters == 3) %>% 
  glmmTMB(
    deaths ~ YEAR * group + (1|SALID1) + offset(population_imp_log),
    data = .,
    family = poisson,
    na.action = na.omit
  )

performance(m1_poisson_w_re_3_glmmTMB) %>% print_html()
parameters(m1_poisson_w_re_3_glmmTMB) %>% print_html()
parameters(m1_poisson_w_re_3_glmmTMB, exponentiate = TRUE) %>% print_html()


plot_model(m1_poisson_w_re_3_glmmTMB, type = "pred",
           terms = c("YEAR",
                     "group")) +
  geom_line(linewidth = 1.2) +
  labs(
    y = "Deaths",
    x = "Year",
    color = "Cluster"
  )  + 
  scale_x_continuous(
    breaks = seq(1997, 2015, 3)
  ) +
  theme_sjplot2()
```

Diagnóstico de supuestos del modelo:

```{r}
check_autocorrelation(m1_poisson_w_re_3_glmmTMB)
check_collinearity(m1_poisson_w_re_3_glmmTMB)
check_overdispersion(m1_poisson_w_re_3_glmmTMB)
simulate_residuals(m1_poisson_w_re_3_glmmTMB) %>% 
  check_residuals()
```


### Clusters: 4

```{r}
m1_poisson_w_re_4_glmmTMB <- gbmt_log_merged %>%
  filter(Clusters == 4) %>% 
  glmmTMB(
    deaths ~ YEAR * group + (1|SALID1) + offset(population_imp_log),
    data = .,
    family = poisson,
    na.action = na.omit
  )

performance(m1_poisson_w_re_4_glmmTMB) %>% print_html()
parameters(m1_poisson_w_re_4_glmmTMB) %>% print_html()
parameters(m1_poisson_w_re_4_glmmTMB, exponentiate = TRUE) %>% print_html()


plot_model(m1_poisson_w_re_4_glmmTMB, type = "pred",
           terms = c("YEAR",
                     "group")) +
  geom_line(linewidth = 1.2) +
  labs(
    y = "Deaths",
    x = "Year",
    color = "Cluster"
  )  + 
  scale_x_continuous(
    breaks = seq(1997, 2015, 3)
  ) +
  theme_sjplot2()
```

Diagnóstico de supuestos del modelo:

```{r}
check_autocorrelation(m1_poisson_w_re_4_glmmTMB)
check_collinearity(m1_poisson_w_re_4_glmmTMB)
check_overdispersion(m1_poisson_w_re_4_glmmTMB)
simulate_residuals(m1_poisson_w_re_4_glmmTMB) %>% 
  check_residuals()
```


## Con efectos aleatorios (YEAR recodificado) 

### Clusters: 3

```{r}
m2_poisson_w_re_3_glmmTMB <- gbmt_log_merged %>%
  filter(Clusters == 3) %>% 
  mutate(
    YEAR = YEAR - 1996,
    group = factor(group, levels = c(3, 1, 2))
  ) %>% 
  glmmTMB(
    deaths ~ YEAR * group + (1|SALID1) + offset(population_imp_log),
    data = .,
    family = poisson,
    na.action = na.omit
  )

performance(m2_poisson_w_re_3_glmmTMB) %>% print_html()
parameters(m2_poisson_w_re_3_glmmTMB) %>% print_html()
parameters(m2_poisson_w_re_3_glmmTMB, exponentiate = TRUE) %>% print_html()


plot_model(m2_poisson_w_re_3_glmmTMB, type = "pred",
           terms = c("YEAR",
                     "group")) +
  geom_line(linewidth = 1.2) +
  labs(
    y = "Deaths",
    x = "Year",
    color = "Cluster"
  )  + 
  scale_x_continuous(
    breaks = seq(1997, 2015, 3)
  ) +
  theme_sjplot2()
```

Diagnóstico de supuestos del modelo:

```{r}
check_autocorrelation(m2_poisson_w_re_3_glmmTMB)
check_collinearity(m2_poisson_w_re_3_glmmTMB)
check_overdispersion(m2_poisson_w_re_3_glmmTMB)
simulate_residuals(m2_poisson_w_re_3_glmmTMB) %>% 
  check_residuals()
```


### Clusters: 4

```{r}
m2_poisson_w_re_4_glmmTMB <- gbmt_log_merged %>%
  filter(Clusters == 4) %>% 
    mutate(
    YEAR = YEAR - 1996,
    group = factor(group, levels = c(4, 1, 2, 3))
  ) %>% 
  glmmTMB(
    deaths ~ YEAR * group + (1|SALID1) + offset(population_imp_log),
    data = .,
    family = poisson,
    na.action = na.omit
  )

performance(m2_poisson_w_re_4_glmmTMB) %>% print_html()
parameters(m2_poisson_w_re_4_glmmTMB) %>% print_html()
parameters(m2_poisson_w_re_4_glmmTMB, exponentiate = TRUE) %>% print_html()


plot_model(m2_poisson_w_re_4_glmmTMB, type = "pred",
           terms = c("YEAR",
                     "group")) +
  geom_line(linewidth = 1.2) +
  labs(
    y = "Deaths",
    x = "Year",
    color = "Cluster"
  )  + 
  scale_x_continuous(
    breaks = seq(1997, 2015, 3)
  ) +
  theme_sjplot2()
```

Diagnóstico de supuestos del modelo:

```{r}
check_autocorrelation(m2_poisson_w_re_4_glmmTMB)
check_collinearity(m2_poisson_w_re_4_glmmTMB)
check_overdispersion(m2_poisson_w_re_4_glmmTMB)
simulate_residuals(m2_poisson_w_re_4_glmmTMB) %>% 
  check_residuals()
```


# Binomial Negativa - Modelo lineal mixto con Efectos Aleatorios 

## Con efectos aleatorios 

### Clusters: 3

```{r}
m1_bn_w_re_3_glmmTMB <- gbmt_log_merged %>%
  filter(Clusters == 3) %>% 
  glmmTMB(
    deaths ~ YEAR * group + (1|SALID1) + offset(population_imp_log),
    data = .,
    family = nbinom2,
    na.action = na.omit
  )

performance(m1_bn_w_re_3_glmmTMB) %>% print_html()
parameters(m1_bn_w_re_3_glmmTMB) %>% print_html()
parameters(m1_bn_w_re_3_glmmTMB, exponentiate = TRUE) %>% print_html()


plot_model(m1_bn_w_re_3_glmmTMB, type = "pred",
           terms = c("YEAR",
                     "group")) +
  geom_line(linewidth = 1.2) +
  labs(
    y = "Deaths",
    x = "Year",
    color = "Cluster"
  )  + 
  scale_x_continuous(
    breaks = seq(1997, 2015, 3)
  ) +
  theme_sjplot2()
```

Diagnóstico de supuestos del modelo:

```{r}
check_autocorrelation(m1_bn_w_re_3_glmmTMB)
check_collinearity(m1_bn_w_re_3_glmmTMB)
check_overdispersion(m1_bn_w_re_3_glmmTMB)
simulate_residuals(m1_bn_w_re_3_glmmTMB) %>% 
  check_residuals()
```


### Clusters: 4

```{r}
m1_nb_w_re_4_glmmTMB <- gbmt_log_merged %>%
  filter(Clusters == 4) %>% 
  glmmTMB(
    deaths ~ YEAR * group + (1|SALID1) + offset(population_imp_log),
    data = .,
    family = nbinom2,
    na.action = na.omit
  )

performance(m1_nb_w_re_4_glmmTMB) %>% print_html()
parameters(m1_nb_w_re_4_glmmTMB) %>% print_html()
parameters(m1_nb_w_re_4_glmmTMB, exponentiate = TRUE) %>% print_html()


plot_model(m1_nb_w_re_4_glmmTMB, type = "pred",
           terms = c("YEAR",
                     "group")) +
  geom_line(linewidth = 1.2) +
  labs(
    y = "Deaths",
    x = "Year",
    color = "Cluster"
  )  + 
  scale_x_continuous(
    breaks = seq(1997, 2015, 3)
  ) +
  theme_sjplot2()
```

Diagnóstico de supuestos del modelo:

```{r}
check_autocorrelation(m1_nb_w_re_4_glmmTMB)
check_collinearity(m1_nb_w_re_4_glmmTMB)
check_overdispersion(m1_nb_w_re_4_glmmTMB)
simulate_residuals(m1_nb_w_re_4_glmmTMB) %>% 
  check_residuals()
```


## Con efectos aleatorios (YEAR recodificado) 

### Clusters: 3

```{r}
m2_nb_w_re_3_glmmTMB <- gbmt_log_merged %>%
  filter(Clusters == 3) %>% 
  mutate(
    YEAR = YEAR - 1996,
    group = factor(group, levels = c(3, 1, 2))
  ) %>% 
  glmmTMB(
    deaths ~ YEAR * group + (1|SALID1) + offset(population_imp_log),
    data = .,
    family = nbinom2,
    na.action = na.omit
  )

performance(m2_nb_w_re_3_glmmTMB) %>% print_html()
parameters(m2_nb_w_re_3_glmmTMB) %>% print_html()
parameters(m2_nb_w_re_3_glmmTMB, exponentiate = TRUE) %>% print_html()


plot_model(m2_nb_w_re_3_glmmTMB, type = "pred",
           terms = c("YEAR",
                     "group")) +
  geom_line(linewidth = 1.2) +
  labs(
    y = "Deaths",
    x = "Year",
    color = "Cluster"
  )  + 
  scale_x_continuous(
    labels = ~ .x + 1996,
    breaks = seq(1, 19, 3)
  ) +
  theme_sjplot2()
```

Diagnóstico de supuestos del modelo:

```{r}
check_autocorrelation(m2_nb_w_re_3_glmmTMB)
check_collinearity(m2_nb_w_re_3_glmmTMB)
check_overdispersion(m2_nb_w_re_3_glmmTMB)
simulate_residuals(m2_nb_w_re_3_glmmTMB) %>% 
  check_residuals() 
```


### Clusters: 4

```{r}
m2_nb_w_re_4_glmmTMB <- gbmt_log_merged %>%
  filter(Clusters == 4) %>% 
    mutate(
    YEAR = YEAR - 1996,
    group = factor(group, levels = c(3, 1, 2, 4))
  ) %>% 
  glmmTMB(
    deaths ~ YEAR * group + (1|SALID1) + offset(population_imp_log),
    data = .,
    family = nbinom2,
    na.action = na.omit
  )

performance(m2_nb_w_re_4_glmmTMB) %>% print_html()
parameters(m2_nb_w_re_4_glmmTMB) %>% print_html()
parameters(m2_nb_w_re_4_glmmTMB, exponentiate = TRUE) %>% print_html()


m2_nb_w_re_4_glmmTMB_plot1 <- plot_model(m2_nb_w_re_4_glmmTMB, type = "pred",
           terms = c("YEAR",
                     "group")) +
  geom_line(linewidth = 1.2) +
  labs(
    y = "Deaths",
    x = "Year",
    color = "Cluster"
  )  + 
  scale_x_continuous(
    labels = ~ .x + 1996,
    breaks = seq(1, 19, 3)
  ) +
  labs(
    title = "Deaths per 100,000 inhabitants by Cluster",
     x = NULL, y = NULL, color = "Cluster", fill = "Cluster"
  ) +
  theme_sjplot2()

m2_nb_w_re_4_glmmTMB_plot1

plot_model(m2_nb_w_re_4_glmmTMB, type = "pred", 
           terms = c("group")) +
  theme_sjplot2()
```


```{r} 
#| eval: false
ggsave(
  filename = "02_output/plots/muertos_nb_cluster4.png",
  plot = m2_nb_w_re_4_glmmTMB_plot1,
  dpi = 300,
  bg = "white",
  width = 7.5,
  height = 5
)
```

Diagnóstico de supuestos del modelo:

```{r}
check_autocorrelation(m2_nb_w_re_4_glmmTMB)
check_collinearity(m2_nb_w_re_4_glmmTMB)
check_overdispersion(m2_nb_w_re_4_glmmTMB)
simulate_residuals(m2_nb_w_re_4_glmmTMB) %>% 
  check_residuals()
```




Ingreso de los valores predictivos para los modelos de 3 y 4 clústers:

```{r}
library(broom.mixed)

gbmt_log_merged4 <- gbmt_log_merged %>% 
  filter(Clusters == 4) %>% 
  augment(m2_nb_w_re_4_glmmTMB, data = ., type.predict = "response")

gbmt_log_merged4 <- gbmt_log_merged4 %>% 
  mutate(
    predicted_deaths_100k = .fitted / population_imp * 100000,
    # Calcula los intervalos de confianza para cada predicción
    lower_ci = (.fitted - 1.96 * .se.fit) / population_imp * 100000,
    upper_ci = (.fitted + 1.96 * .se.fit) / population_imp * 100000
  )
```


```{r}
#| dpi: 300
#| fig-width: 7
#| fig-height: 4.5

plot_binomialoffset_g <- gbmt_log_merged4 %>% 
  ggplot(aes(x = YEAR, 
           y = predicted_deaths_100k,
           color = group)) + 
  geom_point(aes(y = deaths_per_100k), alpha = 0.5,
             position = position_jitter(h = 0.2)) +  # Datos observados
  geom_line(aes(group = SALID1)) +
  labs(
    title = "Deaths per 100,000 inhabitants by Cluster",
     x = NULL, y = NULL, color = "Cluster", fill = "Cluster"
  ) +
  facet_grid(vars(Clusters), vars(group)) +
  scale_color_brewer(palette = "Set1") +
  # facet_wrap(vars(Clusters)) +
  theme_minimal()

plot_binomialoffset_g
```


# Mapa

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

```{r}
salid1_pubsalid1 <- readr::read_csv("01_data/raw/SALID1_PUBSALID1.csv")
gdb_path <- "01_data/raw/Data Request MS242_09222023/MS242_L1.gdb"
gdb_salurbal <- st_layers(gdb_path)

# cities_l1ad <- st_read(
#   gdb_path,
#   gdb_salurbal$name[[1]]
# ) %>%
#   as_tibble() %>%
#   st_as_sf() %>%
#   left_join(
#     salid1_pubsalid1
#   ) %>%
#   relocate(SALID1, .after = PUBSALID1)

cities_l1ad_noislands <- st_read(
  gdb_path,
  gdb_salurbal$name[[2]]
) %>%
  as_tibble() %>%
  st_as_sf() %>%
  left_join(
    salid1_pubsalid1
  ) %>%
  relocate(SALID1, .after = PUBSALID1)
```


```{r}
get_shp <- function(ISO2) {
  # Obtener shp_1
  shp_0 <- geodata::gadm(
    ISO2,
    path = "01_data/raw/shp/",
    level = 0
  ) %>%
    st_as_sf()
  
  # Obtener shp_1
  shp_1 <- geodata::gadm(
    ISO2,
    path = "01_data/raw/shp/",
    level = 1
  ) %>%
    st_as_sf()

  # Obtener shp_2 y aplicar las transformaciones
  shp_2 <- cities_l1ad_noislands %>%
    filter(COUNTRYAB == ISO2) %>%
    janitor::clean_names()

  # Devolver la lista con shp_1 y shp_2
  list(
    shp_0 = shp_0,
    shp_1 = shp_1,
    shp_2 = shp_2
  )
}


shp_list <- unique(gbmt_fit_final$gbmt_fit_total[[1]]$ISO2) %>%
  append(
    c("GUY", "SUR", "BOL", "PRY",
      "URY", "ECU", "BLZ", "HND",
      "TTO", "VEN")
  ) %>% 
  set_names() %>%
  map(get_shp)
```


```{r}
join_data_shp3 <- function(data, shp_list) {
  # Proceso de unión con shapefiles
  data %>%
    mutate(
      pubsalid1 = as.numeric(as.character(pubsalid1)),
      # City = str_to_lower(City),
      group = factor(group)
    ) %>%
    group_nest(ISO2) %>%
    drop_na(ISO2) %>%
    mutate(
      # country = str_to_title(country),
      shp_1 = map(ISO2, ~ shp_list[[.x]]$shp_1),
      shp_2 = map(ISO2, ~ shp_list[[.x]]$shp_2),
      shp_inner = map2(
        data,
        shp_2,
        ~ inner_join(
          .x,
          .y,
          by = join_by(pubsalid1)
        ) %>%
          st_as_sf()
      )
    )
}
```

```{r}
gbmt_interest_deaths <- gbmt_log_merged %>% 
  filter(Clusters == 4) %>% 
  group_by(
    across(c(ISO2:SALID1, group))
  ) %>% 
  summarise(
    mean_deaths_per_100k = mean(deaths_per_100k, na.rm = TRUE)
  ) %>% 
  ungroup()

gbmt_interest_deaths_shp <- map(
  list(gbmt_interest_deaths),
  ~ join_data_shp3(.x, shp_list)
)

city_data_deaths <- gbmt_interest_deaths_shp[[1]] %>% 
  select(ISO2, shp_inner) %>% 
  unnest(cols = shp_inner) %>%
  st_as_sf() %>%
  st_centroid()
```


```{r}
country_data <- map(shp_list, 1) %>% 
  bind_rows() %>% 
  mutate(
    study = c(rep(TRUE, 11), rep(FALSE, 10))
  )

countries_by_study <- country_data %>% 
  as_tibble() %>% 
  select(GID_0, study)

level1_data <- map(shp_list, 2) %>% 
  bind_rows() %>% 
  left_join(countries_by_study)
```

## Biclass plot

```{r}
#| fig-dpi: 600
#| fig-width: 7
#| fig-height: 7

city_data_deaths_biscale <- bi_class(city_data_deaths,
  x = mean_deaths_per_100k,
  y = group,
  style = "quantile", dim = 4
)

custom_pal <- c(
  "1-1" = "#e6ccf5",  # Morado muy claro
  "2-1" = "#c98bd9",  # Morado claro intensificado
  "3-1" = "#8a5599",  # Morado oscuro
  "4-1" = "#562f61",  # Morado muy oscuro
  "1-2" = "#b9f2a1",  # Verde muy claro
  "2-2" = "#6ed965",  # Verde claro intensificado
  "3-2" = "#348c3f",  # Verde oscuro
  "4-2" = "#1e5630",  # Verde muy oscuro
  "1-3" = "#a3e4fd",  # Azul celeste muy claro
  "2-3" = "#5cb8e6",  # Azul claro intensificado
  "3-3" = "#2d6899",  # Azul oscuro
  "4-3" = "#1a4467",  # Azul marino muy oscuro
  "1-4" = "#ffaba8",  # Coral muy claro
  "2-4" = "#ff5747",  # Rojo claro intensificado
  "3-4" = "#c7322c",  # Rojo oscuro
  "4-4" = "#8b1b15"   # Burdeos muy oscuro
)



ggmap <- level1_data %>%
  ggplot() +
  geom_sf(
    aes(fill = study),
    linewidth = 0.05
  ) +
  scale_fill_manual(
    values = c("#bdbdbd", "#e5e5e5"),
    guide = "none"
  ) + 
  geom_sf(
    data = country_data,
    linewidth = 0.35, 
    color = "#3f3f3f",
    fill = NA
  ) +
  geom_sf(
    data = city_data_deaths_biscale,
    mapping = aes(color = bi_class),
    #color = "white", 
    size = 2.5,
    show.legend = FALSE
  ) +
  bi_scale_color(pal = custom_pal, dim = 4) +
  bi_theme()

legend <- bi_legend(
  pal = custom_pal,
  dim = 4,
  xlab = "Mortality per 100K inhabitants",
  ylab = "Cluster (1-4)",
  size = 5
)

x_range_map2 <- city_data_deaths_biscale %>% 
  separate_wider_delim(
    bi_class,
    delim = "-",
    names = c("X", "Y")
  ) %>% 
  group_by(X) %>% 
  summarise(
    min = min(mean_deaths_per_100k),
    max = max(mean_deaths_per_100k)
  ) %>% 
  mutate(
    across(
      min:max,
      round
    ),
    range = paste0("[", min, "; \n", max, "]")
  ) 

legend <- legend + 
  scale_x_continuous(
    labels = c(NA, x_range_map2$range, NA)
  ) +
  scale_y_continuous(
    labels = c(NA, 1:4, NA)
  ) + 
  theme(
    axis.text = element_text(),
    axis.title.y = element_text(
      margin = margin(r = 4)
    ),
    axis.title.x = element_text(
      margin = margin(t = 4)
    )
  )


library(cowplot)

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

finalPlot
```

```{r}
#| eval: false

ggsave("ggmap2.png",
       finalPlot,
       height = 7,
       width = 7,
       dpi = 600)
```

## Gráfico de barras de distribución

```{r}
city_data_deaths_biscale %>% 
  separate_wider_delim(
    cols = bi_class,
    delim = "-",
    names = c("X", "Y")
  ) %>%
  count(X, Y) %>% 
  mutate(
    n = n/sum(n),
    .by = X
  ) %>% 
  mutate(
    X = paste0("Cluster ", X)
  ) %>% 
  ggplot() +
  geom_col(
    aes(
      y = Y,
      x = n,
      fill = Y
    ),
    show.legend = FALSE
  ) +
  scale_fill_brewer(palette = "Set1") +
  # innovar::scale_fill_innova("jama") +
  scale_x_continuous(labels = scales::percent_format()) +
  geom_label(
    aes(
      y = Y,
      x = n,
      label = scales::percent(n, accuracy = 0.1)
    )
  ) +
  facet_wrap(vars(X)) +
  labs(
    y = "Cuartil de Mortalidad"
  ) +
  theme_minimal()
  

ggsave("02_output/plots/Distribución mortalidad y cluster.png",
       dpi = 300,
       bg = "white",
       height = 5,
       width = 8)
```

## Cluster por colores

<!-- jama -->
<!-- mlobo -->
<!-- btran -->
<!-- nasa -->
<!-- mortality -->

```{r}
ggmap2 <- level1_data %>%
  ggplot() +
  geom_sf(
    aes(fill = study),
    linewidth = 0.05
  ) +
  scale_fill_manual(
    values = c("#bdbdbd", "#e5e5e5"),
    guide = "none"
  ) + 
  geom_sf(
    data = country_data,
    linewidth = 0.35, 
    color = "#3f3f3f",
    fill = NA
  ) +
  geom_sf(
    data = city_data_deaths,
    mapping = aes(size = mean_deaths_per_100k,
                  color = group),
    alpha = 0.75
    #color = "white", 
    # size = 2.5,
    # show.legend = FALSE
  ) +
  scale_color_brewer(palette = "Set1") +
  # innovar::scale_color_innova("jama") +
  labs(
    color = "Clusters",
    size = "Mortality per \n100K inhabitants"
  ) +
  theme_minimal()

ggmap2
```


```{r}
#| eval: false

ggsave("ggmap8.png",
       device = grDevices::png,
       ggmap2,
       height = 7,
       width = 10,
       dpi = 600,
       bg = "white")
```

```{r}
gg_stacked_bar <- function(datos) {
  datos %>%
    count(country, group) %>%
    group_by(country) %>%
    mutate(
      percent = n / sum(n),
      country = str_to_title(country)
    ) %>%
    ungroup() %>%
    ggplot(aes(x = country, y = percent)) +
    geom_col(aes(fill = group)) +
    scale_fill_brewer(palette = "Set1") +
    # innovar::scale_fill_innova("jama") +
    labs(
      fill = "Cluster",
      y = NULL,
      x = NULL
    ) +
    geom_label(
      aes(
        group = group,
        label = scales::percent(percent, accuracy = 1)
      ),
      position = position_stack(vjust = 0.5)
    ) +
    scale_y_continuous(
      labels = scales::percent_format()
    ) +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 20, hjust = 1))
}

stacked_bar <- gbmt_log_4 %>% 
  gg_stacked_bar()
```

```{r}
library(patchwork)

map_stacked <- ggmap2 + stacked_bar
map_stacked
```

```{r}
#| eval: false

ggsave("map_stacked.png",
       device = grDevices::png,
       map_stacked,
       height = 7,
       width = 15,
       dpi = 600,
       bg = "white")
```