01-analysis_new.qmd

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

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

```{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 = "cccciiccdddddddddddddddddddddd"
)
```

```{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}
junction <- tidyr::expand_grid(
  ubigeo = unique(tb$ubigeo), year = unique(tb$year), 
  type = unique(centralities$type)
)
```

```{r}
tb_central <- tb |> 
  inner_join(junction, by = c("ubigeo", "year")) |> 
  inner_join(centralities, by = c("ubigeo", "type")) |> 
  mutate(
    tb_incidence = 1000 * tb / population,
    tbfp_incidence = 1000 * tbfp / population,
    
  )
```

## Analysis by districts

```{r}
tb_central_nest <- tidyr::nest(tb_central, data = -c(year, type))
```

```{r}
tb_central_corr <- tb_central_nest |> 
  mutate(
    tb_corr_pearson = purrr::map(
      data, 
      \(x) x |> 
        select(c(tb_incidence, starts_with("centrality"))) |> 
        correlate(method = "pearson", quiet = TRUE) |> 
        focus(tb_incidence) |> 
        arrange(desc(tb_incidence))
    ),
    tbfp_corr_pearson = purrr::map(
      data, 
      \(x) x |> 
        select(c(tbfp_incidence, starts_with("centrality"))) |> 
        correlate(method = "pearson", quiet = TRUE) |> 
        focus(tbfp_incidence) |> 
        arrange(desc(tbfp_incidence))
    ),
    tb_corr_spearman = purrr::map(
      data, 
      \(x) x |> 
        select(c(tb_incidence, starts_with("centrality"))) |> 
        correlate(method = "spearman", quiet = TRUE) |> 
        focus(tb_incidence) |> 
        arrange(desc(tb_incidence))
    ),
    tbfp_corr_spearman = purrr::map(
      data, 
      \(x) x |> 
        select(c(tbfp_incidence, starts_with("centrality"))) |> 
        correlate(method = "spearman", quiet = TRUE) |> 
        focus(tbfp_incidence) |> 
        arrange(desc(tbfp_incidence))
    )
  )
```

```{r}
tb_central_corr$tb_corr_pearson[[1]]
```

```{r}
tb_central_corr$data[[1]] |> 
  ggplot(
    aes(centrality_closeness_all, tb_incidence, color = region, size = population)
  ) +
  geom_point()
```


```{r}
tb_central_corr$tb_corr_spearman[[1]]
```

```{r}
tb_central_corr$data[[1]] |> 
  ggplot(
    aes(log(centrality_strength_out), tb_incidence, color = region, size = population)
  ) +
  geom_point()
```


```{r}
tb_central_corr$tbfp_corr_pearson[[1]]
```

```{r}
tb_central_corr$data[[1]] |> 
  ggplot(
    aes(
      log(centrality_strength_all), tbfp_incidence, color = region, size = population
    )
  ) +
  geom_point()
```

```{r}
tb_central_corr$tbfp_corr_spearman[[1]]
```

```{r}
tb_central_corr$data[[1]] |> 
  ggplot(
    aes(log(centrality_strength_out), tbfp_incidence)
  ) +
  geom_point(aes(color = region, size = population)) 
```

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