01-analysis.qmd

---
title: "Analysis"
author: "Diego Villa"
format: html
editor: source
---

```{r}
library(dplyr)
```

```{r}
# library(tidyr)
# library(ggplot2)
# library(gtsummary)
# library(tidygraph)
# library(ggraph)
# library(igraph)
# library(corrr)
# library(biscale)
# library(performance)
```

```{r}
centralities_raw <- readr::read_csv(
  "data/processed/centralities.csv",
  col_types = sprintf("cccciic%s", paste0(rep("d", 63), collapse = ""))
)
```

```{r}
centralities <- centralities_raw |> 
  mutate(
    region = factor(
      region, 
      levels = c(
        "NORTHERN LIMA", "EASTERN LIMA", "CENTRAL LIMA", "SOUTHERN LIMA", "CALLAO"
      )
    )
  )
```

```{r}
tb <- readr::read_csv("data/processed/tb.csv", col_types = "ciii")
```

```{r}
tb_central <- tb |> 
  inner_join(centralities, by = c("ubigeo")) 
```

## Analysis by subregions 

```{r}
cir_subregions_table <- districts_analysis |>
  gtsummary::tbl_summary(
    by = subregion,
    missing = "no",
    type = cir ~ "continuous2",
    statistic = cir ~ c(
      "{mean} ({sd})", "{median} ({p25}, {p75})", "{min}, {max}"
    ),
    digits = cir ~ 2,
    include = c(subregion, cir)
  )
cir_subregions_table
```

```{r}
districts_analysis |> 
  group_by(subregion) |> 
  summarise(new_cases = sum(new_cases), .groups = "drop") |> 
  mutate(prop_cases = 100 * new_cases / sum(new_cases))
```

```{r}
subregion_cir <- districts_analysis |> 
  group_by(subregion) |> 
  summarise(
    population = sum(population_permanent),
    new_cases = sum(new_cases),
    cir_mean = mean(cir),
    .groups = "drop"
  ) |>
  mutate(
    prop_cases = 100 * new_cases / sum(new_cases),
    cir = 100000 * new_cases / population
  )
subregion_cir
```

Lima Este es la subregión con mayor incidencia acumulada en el 2017, con 131.91 nuevos casos por 100000 habitantes, seguida de Lima Sur con 114.09 nuevos casos por 100000 habitantes, Lima Norte con 103.64 nuevos casos por 100000 habitantes, y por último Lima Centro con 88.53 nuevos casos por 100000 habitantes. 

```{r}
#| fig.height: 3
#| fig.width: 3
nejm_pal <- ggsci::pal_nejm(alpha = 0.9)(8)
# scales::show_col(nejm_pal)
```

```{r}
# subregions_colors = nejm_pal[1:4]
subregions_colors <- c("#0cc0aa", "#154e56", "#b4dbe7", "#539ab3")
names(subregions_colors) <- levels(districts_analysis$subregion)
```

```{r}
# districts_analysis |> 
#   ggplot(aes(subregion, cir, fill = subregion)) +
#   geom_boxplot() +
#   scale_fill_manual(values = subregions_colors, guide = "none") +
#   theme_classic()
```

```{r}
subregion_cir_plot <- subregion_cir |> 
  ggplot(aes(reorder(subregion, cir), cir, fill = subregion)) +
  geom_col() +
  coord_flip() +
  scale_fill_manual(values = subregions_colors) +
  labs(
    x = "Subregiones", y = "Incidencia",
    title = "A) Incidencia anual por 100,000 habitantes en 2017 por subregion en Lima"
  ) +
  theme_classic() +
  theme(legend.position = "none")

subregion_cir_plot
```


```{r}
#| fig.height: 3
#| fig.width: 3
# npg_pal = ggsci::pal_npg(alpha = 1)(9)
# scales::show_col(npg_pal)
```

```{r}
library(sf)
data("Peru", package = "innovar")
peru <- select(Peru, ubigeo, dep, prov, distr)
```

```{r}
districts_analysis_polygons <- districts_analysis |> 
  left_join(peru, by = "ubigeo") |> 
  st_as_sf()
```

```{r}
subregions_polygons <- districts_analysis_polygons |> 
  group_by(subregion) |> 
  summarise() |> 
  ungroup() |> 
  st_as_sf()
```

```{r}
subregion_map <- ggplot(districts_analysis_polygons, aes(fill = subregion)) +
  geom_sf(color = "white", size = 0.05) +
  geom_sf(data = subregions_polygons, color = "black", fill = "transparent") +
  scale_fill_manual(values = subregions_colors) +
  theme_void() +
  labs(title = "A", fill = "Subregion")
subregion_map
```

```{r}
districts_analysis_centroids <- districts_analysis |> 
  left_join(peru, by = "ubigeo") |> 
  st_as_sf() |> 
  st_centroid() 

coordinate_layout <- districts_analysis_centroids |> 
  st_coordinates() |> 
  as_tibble() |> 
  rename_with(tolower)
```

## Commuting flows 

```{r}
edges <- readr::read_csv("data/processed/edges.csv", col_types = "cccd")
```

```{r}
edges |> 
  group_by(type) |> 
  summarise(weight = sum(weight))
```

En cifras absolutas, 2462577 personas se movilizaron para estudiar, 3714535 para trabajar, y 6027237 ya sea para estudiar o trabajar o ambos. 

```{r}
graphs <- edges |>
  rename(weight = weight) |> 
  tidyr::nest(edges = -type) |>
  mutate(
    graph_raw = purrr::map(
      edges, ~ tbl_graph(
        nodes = districts_analysis, edges = .x, directed = TRUE, node_key = "ubigeo"
      )
    ),
    graph = purrr::map(
      graph_raw, ~ as_tbl_graph(simplify(.x, remove.multiple = FALSE))
    ),
    order = purrr::map(graph, ~ gorder(.x)),
    size = purrr::map(graph, ~ gsize(.x)),
    density = purrr::map(graph, ~ edge_density(.x)),
    density_weighted = purrr::map(edges, ~ mean(.x$weight)),
    transitivity = purrr::map(graph, ~ transitivity(.x)),
    diameter = purrr::map(graph, ~ diameter(.x)),
    radius = purrr::map(graph, ~ radius(.x)),
    mean_dist = purrr::map(graph, ~ mean_distance(.x)),
    adhesion = purrr::map(graph, ~ edge_connectivity(.x)),
    cohesion = purrr::map(graph, ~ vertex_connectivity(.x)),
    reciprocity = purrr::map(graph, ~ reciprocity(.x)),
    assortativity_degree = purrr::map(graph, ~ assortativity_degree(.x)),
    assortativity_cir = purrr::map(graph, ~ assortativity(.x, types1 = V(.x)$cir))
  ) 

graphs |> 
  unnest(order:assortativity_cir) |> 
  select(-c(edges:graph))
```

```{r}
total_flow_lima <- ggraph(graphs$graph[[1]], coordinate_layout) +
  geom_sf(
    data = subregions_polygons, aes(fill = subregion), color = "black",
    alpha = 0.5
  ) +
  geom_sf(
    data = districts_analysis_polygons, fill = "transparent", color = "white", 
    size = 0.05
  ) +
  geom_sf(data = subregions_polygons, color = "black", fill = "transparent", linewidth = 0.8) +
  geom_edge_arc(alpha = 0.08, strength = 0.1) +
  geom_node_point() +
  scale_fill_manual(values = subregions_colors) +
  labs(
    title = "Flujo distrito a distrito en la provincia de Lima representado como un grafo", 
    fill = "Subregión"
  ) +
  theme_void()

total_flow_lima
```

```{r}
#| eval: false
ggsave(
  "figures/total_flow_lima.png", 
  total_flow_lima, height = 13, 
  width = 13.5
)
```

```{r}
study_flow_lima <- ggraph(graphs$graph[[2]], coordinate_layout) +
  geom_sf(
    data = subregions_polygons, aes(fill = subregion), color = "black",
    alpha = 0.5
  ) +
  geom_sf(
    data = districts_analysis_polygons, fill = "transparent", color = "white", 
    size = 0.05
  ) +
  geom_sf(data = subregions_polygons, color = "black", fill = "transparent", linewidth = 0.8) +
  geom_edge_arc(alpha = 0.08, strength = 0.1) +
  geom_node_point() +
  scale_fill_manual(values = subregions_colors) +
  labs(
    title = "Figura 1. Flujo de movilización de estudiantes entre distritos de Lima\nrepresentado como un grafo", 
    fill = "Subregión"
  ) +
  theme_void()

study_flow_lima
```

```{r}
#| eval: false
ggsave(
  "figures/study_flow_lima.png", 
  study_flow_lima, height = 11, 
  width = 11.5
)
```

```{r}
work_flow_lima <- ggraph(graphs$graph[[1]], coordinate_layout) +
  geom_sf(
    data = subregions_polygons, aes(fill = subregion), color = "black",
    alpha = 0.5
  ) +
  geom_sf(
    data = districts_analysis_polygons, fill = "transparent", color = "white", 
    size = 0.05
  ) +
  geom_sf(data = subregions_polygons, color = "black", fill = "transparent", linewidth = 0.8) +
  geom_edge_arc(alpha = 0.08, strength = 0.1) +
  geom_node_point() +
  scale_fill_manual(values = subregions_colors) +
  labs(
    title = "Flujo distrito a distrito en la provincia de Lima representado como un grafo", 
    fill = "Subregión"
  ) +
  theme_void()

work_flow_lima
```

```{r}
#| eval: false
ggsave(
  "figures/work_flow_lima.png", 
  work_flow_lima, height = 13, 
  width = 13.5
)
```

### Commuting flows by subregions

```{r}
subregions <- select(districts_analysis, ubigeo, subregion) 
```

```{r}
edges_subregions_raw <- edges |>
  left_join(subregions, by = c("origin" = "ubigeo")) |>
  left_join(
    subregions, by = c("destination" = "ubigeo"),
    suffix = c("_origin", "_destination")
  )
```

```{r}
edges_subregions <- edges_subregions_raw |> 
  group_by(type, subregion_origin, subregion_destination) |>
  summarise(weight = sum(weight), .groups = "drop")
```

```{r}
edges_subregions_total <- edges_subregions |> 
  filter(type == "total") |> 
  select(-type)
```

```{r}
nodes_subregions <- districts_analysis |> 
  group_by(name = subregion) |> 
  summarise(
    n_districts = n(),
    population = sum(population_permanent)
  )
```

```{r}
graphs_subregions <- edges_subregions %>%
  tidyr::nest(edges = -type) %>%
  mutate(
    graph_raw = purrr::map(
      edges, ~ tbl_graph(
        nodes = nodes_subregions, edges = .x, directed = TRUE, 
        node_key = "name"
      )
    ),
    graph = purrr::map(
      graph_raw, ~ as_tbl_graph(simplify(.x, remove.multiple = FALSE))
    )
  )
```

#### Study

```{r}
study_subregions_adjacency <- igraph::as_adjacency_matrix(
   graphs_subregions$graph_raw[[1]], attr = "weight", sparse = FALSE
  )

study_subregions_adjacency
```

```{r}
rowSums(study_subregions_adjacency)
```

```{r}
round(100 * study_subregions_adjacency / rowSums(study_subregions_adjacency), 2)
```

```{r}
round(1000000 * study_subregions_adjacency / nodes_subregions$population, 2)
```


- De la subregión de Lima Norte, un total de 724144 personas se movilizaron para estudiar. De este total, el 85.68% se movilizó hacia la misma subregión y el restante 14.32% hacia otras subregiones. 

- De la subregión de Lima Este, un total de 760302 personas se movilizaron para estudiar. De este total, el 82.05% se movilizó hacia la misma subregión y el restante 17.95% hacia otras subregiones. 

- De la subregión de Lima Centro, un total de 472074 personas se movilizaron para estudiar. De este total, el 88.55% se movilizó hacia la misma subregión y el restante 11.45% hacia otras subregiones.

- De la subregión de Lima Sur, un total de 506057 personas se movilizaron para estudiar. De este total, el 79.08% se movilizó hacia la misma subregión y el restante 20.92% hacia otras subregiones.

```{r}
colSums(study_subregions_adjacency)
```

```{r}
round(100 * study_subregions_adjacency / colSums(study_subregions_adjacency), 2)
```

- Un total de 648166 personas se movilizaron hacia Lima Norte para estudiar. De este total, el 95.73% provenía de la misma subregión y el restante 4.27% de otras subregiones. 

- Un total de 680497 personas se movilizaron hacia Lima Este para estudiar. De este total, el 91.67% provenía de la misma subregión y el restante 8.33% de otras subregiones. 

- Un total de 711062 personas se movilizaron hacia Lima Centro para estudiar. De este total, el 58.79% provenía de la misma subregión y el restante 41.21% de otras subregiones. 

- Un total de 422852 personas se movilizaron hacia Lima Sur para estudiar. De este total, el 94.65% provenía de la misma subregión y el restante 5.35% de otras subregiones. 

```{r}
#| fig.height: 3
#| fig.width: 3
study_flow_chordplot <- ggplotify::as.ggplot(
    ~circlize::chordDiagram(
      filter(
        graphs_subregions$edges[[1]], subregion_origin != subregion_destination
      ),
      directional = 1, 
      annotationTrack = c("name", "grid", "axis"),
      grid.col = subregions_colors, 
      scale = TRUE,
      link.sort = TRUE, link.decreasing = TRUE
    )
  ) + 
  theme_void() +
  labs(title = "B) Flujo de movilización de estudiantes entre subregiones")

study_flow_chordplot
```

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

```{r}
subregion_cir_study_chord <- subregion_cir_plot / study_flow_chordplot +
  plot_layout(heights = c(1, 2.5)) +
  plot_annotation(title = "Figura 2. Incidencia de TB y flujo de movilización de estudiantes\na nivel subregional en Lima")
subregion_cir_study_chord
```

```{r}
ggsave(
  "figures/subregion_cir_study_chord.png", subregion_cir_study_chord,
  width = 7.3, height = 10, bg = "transparent"
)
```

#### Work 

```{r}
work_subregions_adjacency <- igraph::as_adjacency_matrix(
   graphs_subregions$graph_raw[[3]], attr = "weight", sparse = FALSE
  )

work_subregions_adjacency
```

```{r}
rowSums(work_subregions_adjacency)
```

```{r}
round(100 * work_subregions_adjacency / rowSums(work_subregions_adjacency), 2)
```

- De la subregión de Lima Norte, un total de 1008269 personas se movilizaron para trabajar. De este total, el 67.47% se movilizó hacia la misma subregión y el restante 32.53% hacia otras subregiones. 

- De la subregión de Lima Este, un total de 1128303 personas se movilizaron para trabajar. De este total, el 69.81% se movilizó hacia la misma subregión y el restante 30.19% hacia otras subregiones. 

- De la subregión de Lima Centro, un total de 813049 personas se movilizaron para trabajar. De este total, el 87.20% se movilizó hacia la misma subregión y el restante 12.80% hacia otras subregiones.

- De la subregión de Lima Sur, un total de 764914 personas se movilizaron para estudiar trabajar. De este total, el 63.87% se movilizó hacia la misma subregión y el restante 36.13% hacia otras subregiones.

```{r}
colSums(work_subregions_adjacency)
```

```{r}
round(100 * work_subregions_adjacency / colSums(work_subregions_adjacency), 2)
```

- Un total de 728850 personas se movilizaron hacia Lima Norte para trabajar. De este total, el 93.33% provenía de la misma subregión y el restante 6.67% de otras subregiones. 

- Un total de 935641 personas se movilizaron hacia Lima Este para estudiar o trabajar o ambos. De este total, el 91.67% provenía de la misma subregión y el restante 8.33% de otras subregiones. 

- Un total de 1484885 personas se movilizaron hacia Lima Centro para estudiar o trabajar o ambos. De este total, el 47.75% provenía de la misma subregión y el restante 52.25% de otras subregiones. 

- Un total de 565159 personas se movilizaron hacia Lima Sur para estudiar o trabajar o ambos. De este total, el 86.45% provenía de la misma subregión y el restante 13.55% de otras subregiones. 

```{r}
#| fig.height: 3
#| fig.width: 3
work_flow_chordplot <- ggplotify::as.ggplot(
    ~circlize::chordDiagram(
      filter(
        graphs_subregions$edges[[3]], subregion_origin != subregion_destination
      ),
      directional = 1, 
      annotationTrack = c("name", "grid", "axis"),
      grid.col = subregions_colors, 
      scale = TRUE,
      link.sort = TRUE, link.decreasing = TRUE
    )
  ) + 
  theme_void() +
  labs(title = "B) Flujo de movilización hacia centros de trabajo de subregión a subregión")
work_flow_chordplot
```

```{r}
subregion_cir_work_chord <- subregion_cir_plot / work_flow_chordplot +
  plot_layout(heights = c(1, 2.5))
subregion_cir_work_chord
```

```{r}
ggsave(
  "figures/subregion_cir_work_chord.png", subregion_cir_work_chord,
  width = 11.3, height = 15, bg = "transparent"
)
```

#### Study or work or both

```{r}
total_subregions_adjacency <- igraph::as_adjacency_matrix(
   graphs_subregions$graph_raw[[2]], attr = "weight", sparse = FALSE
  )

total_subregions_adjacency
```

```{r}
rowSums(total_subregions_adjacency)
```

```{r}
round(100 * total_subregions_adjacency / rowSums(total_subregions_adjacency), 2)
```

- De la subregión de Lima Norte, un total de 1693254 personas se movilizaron para estudiar o trabajar o para ambas actividades. De este total, el 74.93% se movilizó hacia la misma subregión y el restante 25.07% hacia otras subregiones. 

- De la subregión de Lima Este, un total de 1836660 personas se movilizaron para estudiar o trabajar o para ambas actividades De este total, el 74.47% se movilizó hacia la misma subregión y el restante 25.53% hacia otras subregiones. 

- De la subregión de Lima Centro, un total de 1253601 personas se movilizaron para trabajar. De este total, el 87.56% se movilizó hacia la misma subregión y el restante 12.44% hacia otras subregiones.

- De la subregión de Lima Sur, un total de 1243722 personas se movilizaron para estudiar trabajar. De este total, el 69.75% se movilizó hacia la misma subregión y el restante 30.25% hacia otras subregiones.

```{r}
colSums(total_subregions_adjacency)
```

```{r}
round(100 * total_subregions_adjacency / colSums(total_subregions_adjacency), 2)
```

- Un total de 1343828 personas se movilizaron hacia Lima Norte para trabajar. De este total, el 94.4% provenía de la misma subregión y el restante 5.6% de otras subregiones.

- Un total de 1570226 personas se movilizaron hacia Lima Este para estudiar o trabajar o ambos. De este total, el 87.11% provenía de la misma subregión y el restante 12.89% de otras subregiones. 

- Un total de 2147292 personas se movilizaron hacia Lima Centro para estudiar o trabajar o ambos. De este total, el 51.12% provenía de la misma subregión y el restante 48.88% de otras subregiones. 

- Un total de 965891 personas se movilizaron hacia Lima Sur para estudiar o trabajar o ambos. De este total, el 89.82% provenía de la misma subregión y el restante 10.18% de otras subregiones. 

```{r}
chord_plot <- function() {
  par(cex = 3, mar = c(0, 0, 0, 0))
  circlize::chordDiagram(
    filter(
      graphs_subregions$edges[[2]], subregion_origin != subregion_destination
    ),
    directional = 1, 
    annotationTrack = c("name", "grid", "axis"),
    grid.col = subregions_colors, 
    scale = TRUE,
    link.sort = TRUE, link.decreasing = TRUE
  )
}
total_chord_plot <- chord_plot()
```


```{r}
#| fig.height: 3
#| fig.width: 3
total_flow_chordplot <- ggplotify::as.ggplot(
    ~circlize::chordDiagram(
      filter(
        graphs_subregions$edges[[2]], subregion_origin != subregion_destination
      ),
      directional = 1, 
      annotationTrack = c("name", "grid", "axis"),
      grid.col = subregions_colors, 
      scale = TRUE,
      link.sort = TRUE, link.decreasing = TRUE
    )
  ) + 
  theme_void() +
  labs(title = "B) Flujo de movilización hacia instituciones educativas o centros de trabajo de subregión a subregión")

total_flow_chordplot
```

```{r}
# ggsave(
#   "figures/subregion_flow_chord_map.png", subregion_flow_chord_map, 
#   width = 10.5, height = 6, bg = "transparent"
# )
```

```{r}
subregion_cir_total_chord <- subregion_cir_plot / total_flow_chordplot +
  plot_layout(heights = c(1, 2.5))
subregion_cir_total_chord
```

```{r}
# ggsave(
#   "figures/subregion_cir_total_chord.png", subregion_cir_total_chord,
#   width = 12, height = 6.5, bg = "transparent"
# )
```

```{r}
ggsave(
  "figures/subregion_cir_total_chord.png", subregion_cir_total_chord,
  width = 11.3, height = 15, bg = "transparent"
)
```

## Analysis by districts

```{r}
sum(districts_analysis$new_cases)
```

```{r}
summary(districts_analysis$cir)
sd(districts_analysis$cir)
```

En promedio, hubieron 90 nuevos casos de TB por 100000 habitantes (SD = 48.23) a nivel distrital en el 2017.

```{r}
districts_analysis |> 
  arrange(cir) |> 
  select(district, ubigeo, name_code, cir, subregion)
```

El distrito con mayor incidencia acumulada fue La Victoria en la subregion de Lima Centro, con 212 nuevos casos por 100000 habitantes. El distrito con menos incidencia fue Santa Maria del Mar en la subregión de Lima Sur, con 0 nuevos casos reportados. 

## Commuting flow by districts

### Properties

```{r}
districts_analysis_long <- districts_analysis |> 
  pivot_longer(
    cols = total_centrality_strength_intra:work_structural_constraint_weighted,
    names_to = c("type", "property", "measure", "mode"),
    names_sep = "_"
  ) |> 
  mutate(type = factor(type, levels = c("study", "work", "total")))
```

```{r}
degree_table <- districts_analysis_long |> 
  filter(measure == "degree") |> 
  pivot_wider(
    id_cols = district:property, names_from = mode, values_from = value
  ) |> 
  tbl_summary(
    by = type,
    type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ c(
      "{mean} ({sd})", "{median} ({p25}, {p75})", "{min}, {max}"
    ),
    digits = all_continuous() ~ 2,
    include = c(out:`in`)
  )

degree_table
```

```{r}
strength_table <- districts_analysis_long |> 
  filter(measure == "strength") |> 
  pivot_wider(
    id_cols = district:property, names_from = mode, values_from = value
  ) |> 
  tbl_summary(
    by = type,
    type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ c(
      "{mean} ({sd})", "{median} ({p25}, {p75})", "{min}, {max}"
    ),
    digits = all_continuous() ~ 2,
    include = c(intra, out, `in`),
    missing = "no"
  )

strength_table
```

```{r}
closeness_table <- districts_analysis_long |> 
  filter(measure == "closeness") |> 
  pivot_wider(
    id_cols = district:property, names_from = mode, values_from = value
  ) |> 
  tbl_summary(
    by = type,
    type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ c(
      "{mean} ({sd})", "{median} ({p25}, {p75})", "{min}, {max}"
    ),
    digits = all_continuous() ~ 2,
    include = c(out:`in`)
  )

closeness_table
```

```{r}
other_centralities_table <- districts_analysis_long |> 
  filter(measure %in% c("betweenness", "eigenvector", "pagerank")) |> 
  pivot_wider(
    id_cols = district:property, names_from = measure, values_from = value
  ) |> 
  tbl_summary(
    by = type,
    type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ c(
      "{mean} ({sd})", "{median} ({p25}, {p75})", "{min}, {max}"
    ),
    digits = all_continuous() ~ 4,
    include = c(betweenness, eigenvector, pagerank)
  )

other_centralities_table
```

```{r}
eccentricity_table <- districts_analysis_long |> 
  filter(measure == c("eccentricity")) |> 
  pivot_wider(
    id_cols = district:property, names_from = mode, values_from = value
  ) |> 
  tbl_summary(
    by = type,
    type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ c(
      "{mean} ({sd})", "{median} ({p25}, {p75})", "{min}, {max}"
    ),
    digits = all_continuous() ~ 2,
    include = c(out:`in`)
  )

eccentricity_table
```

```{r}
knn_table <- districts_analysis_long |> 
  filter(measure == c("knn")) |> 
  pivot_wider(
    id_cols = district:property, names_from = mode, values_from = value
  ) |> 
  tbl_summary(
    by = type,
    type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ c(
      "{mean} ({sd})", "{median} ({p25}, {p75})", "{min}, {max}"
    ),
    digits = all_continuous() ~ 2,
    include = c(outout:inin)
  )

knn_table
```

```{r}
constraint_table <- districts_analysis_long |> 
  filter(measure == c("constraint")) |> 
  pivot_wider(
    id_cols = district:property, names_from = mode, values_from = value
  ) |> 
  tbl_summary(
    by = type,
    type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ c(
      "{mean} ({sd})", "{median} ({p25}, {p75})", "{min}, {max}"
    ),
    digits = all_continuous() ~ 2,
    include = c(weighted)
  )

constraint_table
```

### Correlations

#### Study

```{r}
districts_analysis |> 
  select(c(starts_with("study"))) |> 
  select(-c(contains("eccentricity"))) |> 
  select(-c(ends_with("all"))) |> 
  correlate(method = "pearson", quiet = TRUE) |> 
  stretch()
```


```{r}
districts_analysis |> 
  select(c(starts_with("study"))) |> 
  select(-c(contains("eccentricity"))) |> 
  select(-c(ends_with("all"))) |> 
  correlate(method = "pearson", diagonal = 1, quiet = TRUE) |> 
  rearrange() |> 
  stretch() |> 
  ggplot(aes(x, y, fill = r)) +
  geom_tile() +
  theme(axis.text.x=element_text(angle = 90, vjust = 0.5, hjust = 1))
```

```{r}
districts_analysis |> 
  select(c(cir, starts_with("study"))) |> 
  select(-c(contains("eccentricity"))) |> 
  select(-c(ends_with("all"))) |> 
  correlate(method = "pearson", quiet = TRUE) |> 
  focus(cir) |> 
  arrange(desc(cir))
```

```{r}
confintr::ci_cor(
  districts_analysis$cir, 
  districts_analysis$study_centrality_strength_intra,
  type = "bootstrap"
)
```

```{r}
confintr::ci_cor(
  districts_analysis$cir, 
  districts_analysis$study_centrality_strength_intra
)
```


```{r}
districts_analysis |> 
  select(c(cir, starts_with("study"))) |> 
  select(-c(contains("eccentricity"))) |> 
  select(-c(ends_with("all"))) |> 
  correlate(method = "spearman", quiet = TRUE) |> 
  focus(cir) |> 
  arrange(desc(cir))
```

```{r}
plot_scatter <- function(data, x, y, x_label = "x", y_label = "x") {
  ggplot(data, aes({{x}}, {{y}})) +
    geom_point(size = 6) +
    geom_smooth(method = "loess", formula = "y ~ x") +
    labs(x = x_label, y = y_label) +
    theme_classic()
}
```

```{r}
plot_scatter(
  districts_analysis, study_structural_knn_inout, cir,
  "Intra-district flow rate", "Incidence rate"
)
```

```{r}
plot_scatter(
  districts_analysis, study_centrality_degree_out, cir,
  "Intra-district flow rate", "Incidence rate"
)
```

```{r}
plot_scatter(
  districts_analysis, study_centrality_strength_intra, cir,
  "Intra-district flow rate", "Incidence rate"
)
```

```{r}
districts_analysis |> 
  arrange(desc(study_centrality_strength_intra))
```

```{r}
study_centrality_strength_intra_biscale <- bi_class(
  districts_analysis_polygons, x = study_centrality_strength_intra, y = cir, 
  style = "jenks", dim = 3
)
```

```{r}
study_centrality_strength_intra_biscale_map <- study_centrality_strength_intra_biscale |> 
  ggplot() +
  geom_sf(
    aes(fill = bi_class), color = "white", 
    size = 0.05, show.legend = FALSE
  ) +
  geom_sf(data = subregions_polygons, color = "black", fill = "transparent", linewidth = 0.8) +
  bi_scale_fill(pal = "GrPink", dim = 3, flip_axes = TRUE) +
  theme_void() +
  annotate("text", x = -77.13, y = -12.145, label = "Lima\nCentro", size = 5) +
  annotate("text", x = -76.9, y = -12.36, label = "Lima\nSur", size = 5) +
  annotate("text", x = -76.95, y = -11.67, label = "Lima\nNorte", size = 5) +
  annotate("text", x = -76.8, y = -11.85, label = "Lima\nEste", size = 5) +
  labs(title = "Figura 3. Mapa de incidencia de TB vs. porcentaje de estudiantes que\nse movilizan dentro del mismo distrito en el que habitan")

study_centrality_strength_intra_biscale_map
```

```{r}
biscale_legend <- bi_legend(
    pal = "GrPink", dim = 3, xlab = "Mayor centralidad ", ylab = "Mayor incidencia ", 
    size = 10, pad_width = 0.1, flip_axes = TRUE
  ) +
  theme(
    panel.background = element_rect(fill = "transparent"),
    plot.background = element_rect(fill = "transparent", color = NA),
    axis.ticks.x = element_line(color = "transparent"),
    axis.ticks.y = element_line(color = "transparent"),
    axis.line.x = element_line(color = "transparent"),
    axis.line.y = element_line(color = "transparent"),
    text = element_text(size = 16)
  )
```

```{r}
study_centrality_strength_intra_biscale_plot <- cowplot::ggdraw() +
  cowplot::draw_plot(study_centrality_strength_intra_biscale_map, 0, 0, 1, 1) +
  cowplot::draw_plot(biscale_legend, 0.2, 0.08, 0.22, 0.22) +
  # cowplot::draw_label("Title", y = 0.9, x = 0.3) +
  theme(text = element_text(size = 20)) 

study_centrality_strength_intra_biscale_plot
```

```{r}
#| eval: false
ggsave(
  "figures/study_centrality_strength_intra_biscale_plot.png", 
  study_centrality_strength_intra_biscale_plot, height = 11, 
  width = 11
)
```


```{r}
plot_scatter(
  districts_analysis, study_centrality_degree_out, cir,
  "Out-degree centrality", "Incidence rate"
)
```

```{r}
districts_analysis |> 
  arrange(desc(study_centrality_degree_out))
```

```{r}
plot_scatter(
  districts_analysis, study_centrality_degree_in, cir,
  "In-degree centrality", "Incidence rate"
)
```

```{r}
cir_out_strength_scatterplot <- plot_scatter(
  districts_analysis, study_centrality_strength_out, cir,
  "Out-strength centrality", "Incidence rate"
) +
  annotate(
    "text", x = 170000, y = 195, label = "r = -0.58, 95% CI [-0.75, -0.34]",
    size = 9
  ) +
  theme(text = element_text(size = 25))

cir_out_strength_scatterplot
```

```{r}
districts_analysis |> 
  arrange(desc(study_centrality_strength_out))
```

#### Work

```{r}
districts_analysis |> 
  select(c(cir, starts_with("work"))) |> 
  select(-c(contains("eccentricity"))) |> 
  select(-c(ends_with("all"))) |> 
  correlate(method = "pearson", quiet = TRUE) |> 
  focus(cir) |> 
  arrange(desc(cir))
```

```{r}
districts_analysis |> 
  select(c(cir, starts_with("work"))) |> 
  select(-c(contains("eccentricity"))) |> 
  select(-c(ends_with("all"))) |> 
  correlate(method = "spearman", quiet = TRUE) |> 
  focus(cir) |> 
  arrange(desc(cir))
```

```{r}
plot_scatter(
  districts_analysis, work_centrality_closeness_out, cir,
  "Out-closeness centrality", "Incidence rate"
)
```

```{r}
districts_analysis |> 
  arrange(desc(work_centrality_closeness_out))
```

```{r}
plot_scatter(
  districts_analysis, work_structural_knn_inout, cir,
  "Average nearest in-neighbor out-strength", "Incidence rate"
)
```

```{r}
districts_analysis |> 
  arrange(desc(work_structural_knn_inout))
```

```{r}
plot_scatter(
  districts_analysis, work_centrality_strength_in, cir,
  "In-strenght centrality ", "Incidence rate"
)
```

```{r}
districts_analysis |> 
  arrange(desc(work_centrality_strength_in))
```

```{r}
plot_scatter(
  districts_analysis, work_structural_knn_outin, cir,
  "Average nearest out-neighbor in-strength ", "Incidence rate"
)
```

```{r}
districts_analysis |> 
  arrange(desc(work_structural_knn_outin))
```

#### Total

```{r}
districts_analysis |> 
  select(c(cir, starts_with("total"))) |> 
  select(-c(contains("eccentricity"))) |> 
  select(-c(ends_with("all"))) |> 
  correlate(method = "pearson", quiet = TRUE) |> 
  focus(cir) |> 
  arrange(desc(cir))
```

```{r}
districts_analysis |> 
  select(c(cir, starts_with("total"))) |> 
  select(-c(contains("eccentricity"))) |> 
  select(-c(ends_with("all"))) |> 
  correlate(method = "spearman", quiet = TRUE) |> 
  focus(cir) |> 
  arrange(desc(cir))
```

```{r}
plot_scatter(
  districts_analysis, total_centrality_strength_intra, cir,
  "Intra-district flow rate", "Incidence rate"
)
```

```{r}
districts_analysis |> 
  arrange(desc(total_centrality_strength_intra))
```

```{r}
plot_scatter(
  districts_analysis, total_centrality_strength_out, cir,
  "Out-strength centrality", "Incidence rate"
)
```

```{r}
districts_analysis |> 
  arrange(desc(total_centrality_strength_out))
```

```{r}
total_centrality_strength_out_biscale <- bi_class(
  districts_analysis_polygons, x = total_centrality_strength_out, y = cir, 
  style = "jenks", dim = 3
)
```

```{r}
total_centrality_strength_out_biscale_map <- total_centrality_strength_out_biscale |> 
  ggplot() +
  geom_sf(
    aes(fill = bi_class), color = "white", 
    size = 0.05, show.legend = FALSE
  ) +
  geom_sf(data = subregions_polygons, color = "black", fill = "transparent", linewidth = 0.8) +
  bi_scale_fill(pal = "GrPink", dim = 3, flip_axes = TRUE) +
  theme_void() +
  annotate("text", x = -77.13, y = -12.145, label = "Lima\nCentro", size = 5) +
  annotate("text", x = -76.9, y = -12.36, label = "Lima\nSur", size = 5) +
  annotate("text", x = -76.95, y = -11.67, label = "Lima\nNorte", size = 5) +
  annotate("text", x = -76.8, y = -11.85, label = "Lima\nEste", size = 5)

total_centrality_strength_out_biscale_map
```

```{r}
biscale_legend <- bi_legend(
    pal = "GrPink", dim = 3, xlab = "Mayor centralidad ", ylab = "Mayor incidencia ", 
    size = 10, pad_width = 0.1, flip_axes = TRUE
  ) +
  theme(
    panel.background = element_rect(fill = "transparent"),
    plot.background = element_rect(fill = "transparent", color = NA),
    axis.ticks.x = element_line(color = "transparent"),
    axis.ticks.y = element_line(color = "transparent"),
    axis.line.x = element_line(color = "transparent"),
    axis.line.y = element_line(color = "transparent"),
    text = element_text(size = 16)
  )
```

```{r}
total_centrality_strength_out_biscale_plot <- cowplot::ggdraw() +
  cowplot::draw_plot(total_centrality_strength_out_biscale_map, 0, 0, 1, 1) +
  cowplot::draw_plot(biscale_legend, 0.13, 0.1, 0.25, 0.25) +
  theme(text = element_text(size = 20)) 

total_centrality_strength_out_biscale_plot
```

```{r}
#| eval: false
ggsave(
  "figures/total_centrality_strength_out_biscale_plot.png", 
  total_centrality_strength_out_biscale_plot, height = 13, 
  width = 13
)
```

## Model building

```{r}
districts_model <- districts_analysis |> 
  mutate(
    across(c(co, no2, pm25, o3), \(x) as.vector(scale(x)))
  )
```


```{r}
glm_base_formulas_list <- list(
  "poverty-co" = 
    new_cases ~ offset(log(population_permanent)) + monetary_poverty + co,
  "poverty-no2" = 
    new_cases ~ offset(log(population_permanent)) + monetary_poverty + no2,
  "poverty-pm25" = 
    new_cases ~ offset(log(population_permanent)) + monetary_poverty + pm25,
  "poverty-o3" = 
    new_cases ~ offset(log(population_permanent)) + monetary_poverty + o3,
  "nonmonetary-co" = 
    new_cases ~ offset(log(population_permanent)) + nonmonetary_poverty + co,
  "nonmonetary-no2" = 
    new_cases ~ offset(log(population_permanent)) + nonmonetary_poverty + no2,
  "nonmonetary-pm25" = 
    new_cases ~ offset(log(population_permanent)) + nonmonetary_poverty + pm25,
  "nonmonetary-o3" = 
    new_cases ~ offset(log(population_permanent)) + nonmonetary_poverty + o3
)
```

```{r}
glm_base_formulas <- tibble(
  base = names(glm_base_formulas_list), 
  formula = glm_base_formulas_list
)
```

### GLM

#### Study

```{r}
study_mobility <- 
  colnames(districts_model)[grepl("study_centrality", colnames(districts_model))]
```

```{r}
glm_study_mobility_design <- expand.grid(
    base = names(glm_base_formulas_list), mobility = study_mobility
  ) |> 
  arrange(base)
```

```{r}
glm_study_mobility_formulas_raw <- glm_study_mobility_design |> 
  inner_join(glm_base_formulas, by = c("base"))
```

```{r}
glm_study_mobility_formulas <- glm_study_mobility_formulas_raw |> 
  mutate(
    formula = purrr::map2(
      formula, as.list(mobility),
      ~update.formula(.x, as.formula(sprintf("~ . + %s", .y)))
    )
  ) 
```

```{r}
glm_study_formulas <- bind_rows(glm_base_formulas, glm_study_mobility_formulas)
```

```{r}
glm_study_fits <- glm_study_formulas |> 
  mutate(
    glm_po_fit = purrr::map(
      formula, ~glm(.x, data = districts_model, family = poisson)
    ),
    glm_nb_fit = purrr::map(
      formula, ~MASS::glm.nb(.x, data = districts_model)
    ),
  ) |> 
  tidyr::pivot_longer(
    -c(base, mobility, formula),
    names_to = c("model", "family"),
    names_pattern = "(.*)_(.*)_fit",
    values_to = "fit"
  )
```

```{r}
glm_study_performance <- glm_study_fits |> 
  mutate(aic = purrr::map(fit, ~AIC(.x))) |> 
  tidyr::unnest(aic) |> 
  arrange(aic)

glm_study_performance
```

```{r}
glm_study_performance_base <- glm_study_performance |> 
  filter(is.na(mobility)) |> 
  select(base, family, aic_base = aic) |> 
  arrange(aic_base)

glm_study_performance_base
```

```{r}
glm_study_performance |> 
  drop_na(mobility) |> 
  left_join(glm_study_performance_base, by = c("base", "family")) |> 
  mutate(aic_change = 100 * (aic - aic_base) / aic_base) |> 
  arrange(aic_change) |> 
  filter(family == "nb", base == "poverty-no2")
```

```{r}
summary(glm_study_performance$fit[[1]])
```

```{r}
gratia::appraise(glm_study_performance$fit[[1]])
```

```{r}
check_overdispersion(glm_study_performance$fit[[1]])
```

```{r}
check_collinearity(glm_study_performance$fit[[1]])
```

```{r}
check_model(glm_study_performance$fit[[1]])
```


#### Work

```{r}
work_mobility <- 
  colnames(districts_model)[grepl("work_centrality", colnames(districts_model))]
```

```{r}
glm_work_mobility_design <- expand.grid(
    base = names(glm_base_formulas_list), mobility = work_mobility
  ) |> 
  arrange(base)
```

```{r}
glm_work_mobility_formulas_raw <- glm_work_mobility_design |> 
  inner_join(glm_base_formulas, by = c("base"))
```

```{r}
glm_work_mobility_formulas <- glm_work_mobility_formulas_raw |> 
  mutate(
    formula = purrr::map2(
      formula, as.list(mobility),
      ~update.formula(.x, as.formula(sprintf("~ . + %s", .y)))
    )
  ) 
```

```{r}
glm_work_formulas <- bind_rows(glm_base_formulas, glm_work_mobility_formulas)
```

```{r}
glm_work_fits <- glm_work_formulas |> 
  mutate(
    glm_po_fit = purrr::map(
      formula, ~glm(.x, data = districts_model, family = poisson)
    ),
    glm_nb_fit = purrr::map(
      formula, ~MASS::glm.nb(.x, data = districts_model)
    ),
  ) |> 
  tidyr::pivot_longer(
    -c(base, mobility, formula),
    names_to = c("model", "family"),
    names_pattern = "(.*)_(.*)_fit",
    values_to = "fit"
  )
```

```{r}
glm_work_performance <- glm_work_fits |> 
  mutate(aic = purrr::map(fit, ~AIC(.x))) |> 
  tidyr::unnest(aic) |> 
  arrange(aic)

glm_work_performance
```

```{r}
glm_work_performance_base <- glm_work_performance |> 
  filter(is.na(mobility)) |> 
  select(base, family, aic_base = aic) |> 
  arrange(aic_base)

glm_work_performance_base
```

```{r}
glm_work_performance |> 
  drop_na(mobility) |> 
  left_join(glm_work_performance_base, by = c("base", "family")) |> 
  mutate(aic_change = 100 * (aic - aic_base) / aic_base) |> 
  arrange(aic_change) |> 
  filter(family == "nb")
```

```{r}
summary(glm_work_performance$fit[[1]])
```

#### Total

```{r}
total_mobility <- 
  colnames(districts_model)[grepl("total_centrality", colnames(districts_model))]
```

```{r}
glm_total_mobility_design <- expand.grid(
    base = names(glm_base_formulas_list), mobility = total_mobility
  ) |> 
  arrange(base)
```

```{r}
glm_total_mobility_formulas_raw <- glm_total_mobility_design |> 
  inner_join(glm_base_formulas, by = c("base"))
```

```{r}
glm_total_mobility_formulas <- glm_total_mobility_formulas_raw |> 
  mutate(
    formula = purrr::map2(
      formula, as.list(mobility),
      ~update.formula(.x, as.formula(sprintf("~ . + %s", .y)))
    )
  ) 
```

```{r}
glm_total_formulas <- bind_rows(glm_base_formulas, glm_total_mobility_formulas)
```

```{r}
glm_total_fits = glm_total_formulas |> 
  mutate(
    glm_po_fit = purrr::map(
      formula, ~glm(.x, data = districts_model, family = poisson)
    ),
    glm_nb_fit = purrr::map(
      formula, ~MASS::glm.nb(.x, data = districts_model)
    ),
  ) |> 
  tidyr::pivot_longer(
    -c(base, mobility, formula),
    names_to = c("model", "family"),
    names_pattern = "(.*)_(.*)_fit",
    values_to = "fit"
  )
```

```{r}
glm_total_performance <- glm_total_fits |> 
  mutate(aic = purrr::map(fit, ~AIC(.x))) |> 
  tidyr::unnest(aic) |> 
  arrange(aic)

glm_total_performance
```

```{r}
glm_total_performance_base <- glm_total_performance |> 
  filter(is.na(mobility)) |> 
  select(base, family, aic_base = aic) |> 
  arrange(aic_base)

glm_total_performance_base
```

```{r}
glm_total_performance |> 
  drop_na(mobility) |> 
  left_join(glm_total_performance_base, by = c("base", "family")) |> 
  mutate(aic_change = 100 * (aic - aic_base) / aic_base) |> 
  arrange(aic_change) |> 
  filter(family == "nb")
```

```{r}
summary(glm_total_performance$fit[[1]])
```

## GAM

```{r}
gam_base_formulas_list <- list(
  "poverty-co" = 
    new_cases ~ offset(log(population_permanent)) + s(monetary_poverty) + s(co),
  "poverty-no2" = 
    new_cases ~ offset(log(population_permanent)) + s(monetary_poverty) + 
    s(no2),
  "poverty-pm25" = 
    new_cases ~ offset(log(population_permanent)) + s(monetary_poverty) + 
    s(pm25),
  "poverty-o3" = 
    new_cases ~ offset(log(population_permanent)) + s(monetary_poverty) + 
    s(o3),
  "nonmonetary-co" = 
    new_cases ~ offset(log(population_permanent)) + s(nonmonetary_poverty) + s(co),
  "nonmonetary-no2" = 
    new_cases ~ offset(log(population_permanent)) + s(nonmonetary_poverty) + s(no2),
  "nonmonetary-pm25" = 
    new_cases ~ offset(log(population_permanent)) + s(nonmonetary_poverty) + s(pm25),
  "nonmonetary-o3" = 
    new_cases ~ offset(log(population_permanent)) + s(nonmonetary_poverty) + s(o3)
)
```

```{r}
gam_base_formulas <- tibble(
  base = names(gam_base_formulas_list), 
  formula = gam_base_formulas_list
)
```

### Study

```{r}
gam_study_mobility_design <- expand.grid(
    base = names(gam_base_formulas_list), mobility = study_mobility
  ) |> 
  arrange(base)
```

```{r}
gam_study_mobility_formulas_raw <- gam_study_mobility_design |> 
  inner_join(gam_base_formulas, by = c("base"))
```

```{r}
gam_study_mobility_formulas <- gam_study_mobility_formulas_raw |> 
  mutate(
    formula = purrr::map2(
      formula, as.list(mobility),
      ~update.formula(.x, as.formula(sprintf("~ . + s(%s, k = 5)", .y)))
    )
  ) 
```

```{r}
gam_study_formulas <- bind_rows(gam_base_formulas, gam_study_mobility_formulas)
```

```{r}
gam_study_fits <- gam_study_formulas |> 
  mutate(
    gam_po_fit = purrr::map(
      formula, ~mgcv::gam(
        .x, data = districts_model, family = poisson, method = "ML"
      )
    ),
    gam_nb_fit = purrr::map(
      formula, ~mgcv::gam(
        .x, data = districts_model, family = mgcv::nb(), method = "ML"
      )
    ),
  ) |> 
  tidyr::pivot_longer(
    -c(base, mobility, formula),
    names_to = c("model", "family"),
    names_pattern = "(.*)_(.*)_fit",
    values_to = "fit"
  )
```

```{r}
gam_study_performance <- gam_study_fits |> 
  mutate(aic = purrr::map(fit, ~AIC(.x))) |> 
  tidyr::unnest(aic) |> 
  arrange(aic)

gam_study_performance
```

```{r}
gam_study_performance_base <- gam_study_performance |> 
  filter(is.na(mobility)) |> 
  select(base, family, aic_base = aic) |> 
  arrange(aic_base)

gam_study_performance_base
```

```{r}
gam_study_performance |> 
  drop_na(mobility) |> 
  left_join(gam_study_performance_base, by = c("base", "family")) |> 
  mutate(aic_change = aic - aic_base) |> 
  arrange(aic_change) 
  # filter(family == "nb", base == "nonmonetary-co")
```

```{r}
summary(gam_study_performance$fit[[1]])
```

```{r}
gratia::appraise(gam_study_performance$fit[[1]])
```


```{r}
gratia::draw(gam_study_performance$fit[[1]])
```

```{r}
gam_study_performance_fit_est <- gam_study_performance$fit[[1]] |> 
  gratia::smooth_estimates() |> 
  gratia::add_confint()
```

```{r}
study_poverty_plot <- gam_study_performance_fit_est |> 
  mutate(across(c(est, lower_ci, upper_ci), exp)) |> 
  filter(smooth == "s(ubn)") |> 
  ggplot() +
  geom_line(aes(ubn, est), linewidth = 1) +
  geom_ribbon(
    aes(ymin = lower_ci, ymax = upper_ci, x = ubn), alpha = 0.4,
    fill = "steelblue3"
  ) +
  geom_rug(
    aes(ubn), data = districts_model, sides = "b", length = grid::unit(0.02, "npc")
  ) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  labs(y = "Efecto parcial", x = "Porcentaje de hogares con al menos una necesidad básica insatisfecha") +
  theme_classic()
study_poverty_plot
```

```{r}
study_pollulant_plot <- gam_study_performance_fit_est |> 
  mutate(across(c(est, lower_ci, upper_ci), exp)) |> 
  filter(smooth == "s(no2)") |> 
  ggplot() +
  geom_line(aes(no2, est), linewidth = 1) +
  geom_ribbon(
    aes(ymin = lower_ci, ymax = upper_ci, x = no2), alpha = 0.4,
    fill = "steelblue3"
  ) +
  geom_rug(
    aes(no2), data = districts_model, sides = "b", length = grid::unit(0.02, "npc")
  ) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  labs(y = "Efecto parcial", x = "Concentración de NO2 (estandarizado)") +
  theme_classic()
study_pollulant_plot
```

```{r}
study_centrality_plot <- gam_study_performance_fit_est |> 
  mutate(across(c(est, lower_ci, upper_ci), exp)) |> 
  filter(smooth == "s(study_centrality_strength_intra)") |> 
  ggplot() +
  geom_line(aes(study_centrality_strength_intra, est), linewidth = 1) +
  geom_ribbon(
    aes(ymin = lower_ci, ymax = upper_ci, x = study_centrality_strength_intra), alpha = 0.4,
    fill = "steelblue3"
  ) +
  geom_rug(
    aes(study_centrality_strength_intra), data = districts_model, sides = "b", length = grid::unit(0.02, "npc")
  ) +
  scale_x_continuous(labels = scales::comma) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  labs(y = "Efecto parcial", x = "Porcentaje de estudiantes que se movilizan dentro del distrito donde habitan") +
  theme_classic()
study_centrality_plot
```

```{r}
study_smooth_plots <- study_poverty_plot + study_pollulant_plot + 
  study_centrality_plot + plot_layout(ncol = 2) +
  plot_annotation(title = "Figura 4. Efectos parciales del mejor modelo para la incidencia de TB con una medida de centralidad del grafo de movilización de estudiantes como predictor")
study_smooth_plots
```

```{r}
ggsave(
  "figures/study_smooth_plots.png", study_smooth_plots,
  width = 13, height = 9, bg = "transparent"
)
```

```{r}
gratia::draw(gam_study_performance$fit[[1]])
```


### Work

```{r}
gam_work_mobility_design <- expand.grid(
    base = names(gam_base_formulas_list), mobility = work_mobility
  ) |> 
  arrange(base)
```

```{r}
gam_work_mobility_formulas_raw <- gam_work_mobility_design |> 
  inner_join(gam_base_formulas, by = c("base"))
```

```{r}
gam_work_mobility_formulas <- gam_work_mobility_formulas_raw |> 
  mutate(
    formula = purrr::map2(
      formula, as.list(mobility),
      ~update.formula(.x, as.formula(sprintf("~ . + s(%s, k = 5)", .y)))
    )
  ) 
```

```{r}
gam_work_formulas <- bind_rows(gam_base_formulas, gam_work_mobility_formulas)
```

```{r}
gam_work_fits <- gam_work_formulas |> 
  mutate(
    gam_po_fit = purrr::map(
      formula, ~mgcv::gam(
        .x, data = districts_model, family = poisson, method = "ML"
      )
    ),
    gam_nb_fit = purrr::map(
      formula, ~mgcv::gam(
        .x, data = districts_model, family = mgcv::nb(), method = "ML"
      )
    ),
  ) |> 
  tidyr::pivot_longer(
    -c(base, mobility, formula),
    names_to = c("model", "family"),
    names_pattern = "(.*)_(.*)_fit",
    values_to = "fit"
  )
```

```{r}
gam_work_performance <- gam_work_fits |> 
  mutate(aic = purrr::map(fit, ~AIC(.x))) |> 
  tidyr::unnest(aic) |> 
  arrange(aic)

gam_work_performance
```

```{r}
gam_work_performance |> 
  drop_na(mobility) |> 
  left_join(gam_study_performance_base, by = c("base", "family")) |> 
  mutate(aic_change = aic - aic_base) |> 
  arrange(aic_change) 
  # filter(family == "nb", base == "nonmonetary-co")
```

```{r}
gratia::draw(gam_work_performance$fit[[1]])
```

### Total

```{r}
gam_total_mobility_design <- expand.grid(
    base = names(gam_base_formulas_list), mobility = total_mobility
  ) |> 
  arrange(base)
```

```{r}
gam_total_mobility_formulas_raw <- gam_total_mobility_design |> 
  inner_join(gam_base_formulas, by = c("base"))
```

```{r}
gam_total_mobility_formulas <- gam_total_mobility_formulas_raw |> 
  mutate(
    formula = purrr::map2(
      formula, as.list(mobility),
      ~update.formula(.x, as.formula(sprintf("~ . + s(%s, k = 5)", .y)))
    )
  ) 
```

```{r}
gam_total_formulas <- bind_rows(gam_base_formulas, gam_total_mobility_formulas)
```

```{r}
gam_total_fits <- gam_total_formulas |> 
  mutate(
    gam_po_fit = purrr::map(
      formula, ~mgcv::gam(
        .x, data = districts_model, family = poisson, method = "ML"
      )
    ),
    gam_nb_fit = purrr::map(
      formula, ~mgcv::gam(
        .x, data = districts_model, family = mgcv::nb(), method = "ML"
      )
    ),
  ) |> 
  tidyr::pivot_longer(
    -c(base, mobility, formula),
    names_to = c("model", "family"),
    names_pattern = "(.*)_(.*)_fit",
    values_to = "fit"
  )
```

```{r}
gam_total_performance <- gam_total_fits |> 
  mutate(aic = purrr::map(fit, ~AIC(.x))) |> 
  tidyr::unnest(aic) |> 
  arrange(aic)

gam_total_performance
```

```{r}
gam_total_performance_base <- gam_total_performance |> 
  filter(is.na(mobility)) |> 
  select(base, family, aic_base = aic) |> 
  arrange(aic_base)

gam_total_performance_base
```

```{r}
gam_total_performance |> 
  drop_na(mobility) |> 
  left_join(gam_study_performance_base, by = c("base", "family")) |> 
  mutate(aic_change = aic - aic_base) |> 
  arrange(aic_change) 
  # filter(family == "nb", base == "nonmonetary-co")
```

```{r}
summary(gam_total_performance$fit[[1]])
```

```{r}
gratia::draw(gam_total_performance$fit[[1]])
```

```{r}
gam_total_performance_fit_est <- gam_total_performance$fit[[1]] |> 
  gratia::smooth_estimates() |> 
  gratia::add_confint()
```

```{r}
ubn_plot <- gam_total_performance_fit_est |> 
  mutate(across(c(est, lower_ci, upper_ci), exp)) |> 
  filter(smooth == "s(ubn)") |> 
  ggplot() +
  geom_line(aes(ubn, est), linewidth = 1) +
  geom_ribbon(
    aes(ymin = lower_ci, ymax = upper_ci, x = ubn), alpha = 0.4,
    fill = "steelblue3"
  ) +
  geom_rug(
    aes(ubn), data = districts_model, sides = "b", length = grid::unit(0.02, "npc")
  ) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  labs(y = "Efecto parcial", x = "Hogares con al menos una necesidad básica insatisfecha (%)") +
  theme_classic()
ubn_plot
```

```{r}
co_plot <- gam_total_performance_fit_est |> 
  mutate(across(c(est, lower_ci, upper_ci), exp)) |> 
  filter(smooth == "s(co)") |> 
  ggplot() +
  geom_line(aes(co, est), linewidth = 1) +
  geom_ribbon(
    aes(ymin = lower_ci, ymax = upper_ci, x = co), alpha = 0.4,
    fill = "steelblue3"
  ) +
  geom_rug(
    aes(co), data = districts_model, sides = "b", length = grid::unit(0.02, "npc")
  ) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  labs(y = "Efecto parcial", x = "Concentración de CO (estandarizado)") +
  theme_classic()
co_plot
```


```{r}
strength_plot <- gam_total_performance_fit_est |> 
  mutate(across(c(est, lower_ci, upper_ci), exp)) |> 
  filter(smooth == "s(total_centrality_strength_intra)") |> 
  ggplot() +
  geom_line(aes(total_centrality_strength_intra, est), linewidth = 1) +
  geom_ribbon(
    aes(ymin = lower_ci, ymax = upper_ci, x = total_centrality_strength_intra), alpha = 0.4,
    fill = "steelblue3"
  ) +
  geom_rug(
    aes(total_centrality_strength_intra), data = districts_model, sides = "b", length = grid::unit(0.02, "npc")
  ) +
  scale_x_continuous(labels = scales::comma) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  labs(y = "Efecto parcial", x = "Personas que se movilizan dentro del mismo distrito por cada 100,000 habitantes") +
  theme_classic()
strength_plot
```

```{r}
smooth_plots <- ubn_plot + co_plot + strength_plot + plot_layout(ncol = 2)
smooth_plots
```

```{r}
ggsave(
  "figures/smooth_plots.png", smooth_plots,
  width = 15, height = 12, bg = "transparent"
)
```



```{r}
gratia::draw(gam_total_performance$fit[[1]]) &
  theme_classic()
```