00-network-analysis.Rmd

---
title: "Network analysis"
author: "Villa-Almeyda, Diego"
date: "`r Sys.Date()`"
output: html_document
---

```{r message=FALSE}
library(dplyr)
library(tidygraph)
library(ggraph)
library(gtsummary)
```

# Data reading

```{r}
processed_path <- "data/processed/"
```

```{r}
nodes_path <- fs::path(processed_path, "00-nodes.csv")
nodes <- readr::read_csv(nodes_path, col_types = "ccccdcd")
glimpse(nodes)
```

```{r}
edges_path <- fs::path(processed_path, "01-edges.csv")
edges <- readr::read_csv(edges_path, col_types = "cccd")
glimpse(edges)
```

```{r}
nodes_watershed <- nodes |> 
  mutate(watershed_ana = watershed) |> 
  tidyr::nest(nodes = -watershed_ana) 
```

```{r}
edges_watershed <- edges |> 
  mutate(watershed_ana = watershed) |> 
  tidyr::nest(edges = -watershed_ana) 
```

```{r}
networks <- nodes_watershed |>  
  inner_join(edges_watershed, by = "watershed_ana") |> 
  mutate(
    networks = purrr::map2(
      nodes, edges, 
      ~ tbl_graph(
        directed = FALSE, node_key = "village_id", nodes = .x, edges = .y
      )
    )
  ) |> 
  mutate(
    order = purrr::map(networks, ~ igraph::gorder(.x)),
    size = purrr::map(networks, ~ igraph::gsize(.x)),
    edge_density = purrr::map(networks, ~ igraph::edge_density(.x))
  ) |> 
  tidyr::unnest(c(order, size, edge_density))
```

Strength and eigenvector centrality measures interpret weights as connection strength (higher weights, more connected), whereas closeness and betweenness centrality measures interpret weights as distances (higher weights, less connected).

```{r}
source("R/get_centralities.R")
```

```{r warning=FALSE}
centralities <- networks |> 
  mutate(centralities = purrr::map(networks, get_centralities))
```

```{r}
village_centralities <- centralities |> 
  mutate(
    centrality_tables = purrr::map(centralities, \(x) as_tibble(activate(x, nodes)))
  ) |> 
  select(centrality_tables) |> 
  tidyr::unnest(centrality_tables)
```

```{r include=FALSE}
centralities_filepath <- fs::path(processed_path, "02-village-centralities.csv")
readr::write_csv(village_centralities, centralities_filepath, na = "")
```

# Summary tables

```{r}
centrality_labels <- list(
  dist_stre = "Strength D",
  dist_pop_gravity_stre = "Strength GM-PD",
  dist_forest_gravity_stre = "Strength GM-FD",
  time_stre = "Strength T",
  time_pop_gravity_stre = "Strength GM-PT",
  time_forest_gravity_stre = "Strength GM-FT",
  dist_close = "Closeness D",
  dist_pop_gravity_close = "Closeness GM-PD",
  dist_forest_gravity_close = "Closeness GM-FD",
  time_close = "Closeness T",
  time_pop_gravity_close = "Closeness GM-PT",
  time_forest_gravity_close = "Closeness GM-FT",
  dist_between = "Betweenness D",
  dist_pop_gravity_between = "Betweenness GM-PD",
  dist_forest_gravity_between = "Betweenness GM-FD",
  time_between = "Betweenness T",
  time_pop_gravity_between = "Betweenness GM-PT",
  time_forest_gravity_between = "Betweenness GM-FT",
  dist_eigen = "Eigenvector D",
  dist_pop_gravity_eigen = "Eigenvector GM-PD",
  dist_forest_gravity_eigen = "Eigenvector GM-FD",
  time_eigen = "Eigenvector T",
  time_pop_gravity_eigen = "Eigenvector GM-PT",
  time_forest_gravity_eigen = "Eigenvector GM-FT"
)
```

```{r}
centralities_summary <- village_centralities |> 
  select(watershed, dist_stre:time_forest_gravity_eigen) |> 
  tbl_summary(
    by = watershed,
    statistic = all_continuous() ~ c("{mean} ({sd})"),
    digits = all_continuous() ~ c(2, 2),
    type = everything() ~ "continuous",
    label = centrality_labels
  ) |> 
  add_overall(last = TRUE) |> 
  modify_footnote(all_stat_cols() ~ NA) |> 
  modify_header(
    label = "Centrality index", all_stat_cols() ~ "{level}, N = {n}"
  ) |> 
  as_tibble()
```

```{r}
centralities_summary_path <- fs::path("tables/centralities-summary.csv")
readr::write_csv(centralities_summary, centralities_summary_path)
```

```{r}
weight_labels <- list(
  distance = "Distance (D)",
  dist_pop_gravity = "Gravity model: Population and distance (GM-PD)",
  dist_forest_gravity = "Gravity model: Forest loss and distance (GM-FD)",
  time = "Travel time (T)",
  time_pop_gravity = "Gravity model: Population and travel time (GM-PT)",
  time_forest_gravity = "Gravity model: Forest loss and travel time (GM-FT)"
)
```

```{r}
weights_summary <- edges |> 
  select(
    watershed, distance:dist_forest_gravity, time, 
    time_pop_gravity:time_forest_gravity
  ) |> 
  tbl_summary(
    by = watershed,
    statistic = all_continuous() ~ c("{mean} ({sd})"),
    digits = all_continuous() ~ c(2, 2),
    type = everything() ~ "continuous",
    label = weight_labels
  ) |> 
  add_overall(last = TRUE) |> 
  modify_footnote(all_stat_cols() ~ NA) |> 
  modify_header(
    label = "Edge weight", all_stat_cols() ~ "{level}, N = {n}"
  ) |> 
  as_tibble()
```

```{r}
weights_summary_path <- fs::path("tables/weights-summary.csv")
readr::write_csv(weights_summary, weights_summary_path)
```

# Graph visualization

## Single graphs

```{r}
graphs <- purrr::reduce(centralities$centralities, bind_graphs)
```

```{r}
set.seed(1)
base_palette <- sample(
  c(
    "#F94144", "#F3722C", "#F8961E", "#F9844A", "#F9C74F", "#90BE6D", "#43AA8B",
    "#4D908E", "#577590", "#277DA1", "#F7A072", "#EDDEA4", "#8EA604", "#00CFC1",
    "#485696", "#37323E", "#985277", "#5C374C", "#7678ED", "#CE4257", "#FFBD00"
  )
)
graph_palette <- colorRampPalette(base_palette)
```

```{r}
layout_kk <- create_layout(graphs, layout = "kk")
```

```{r}
source("R/plot_single_graph.R")
```

### Strength

```{r}
plot_dist_stre <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = distance, 
  centrality = dist_stre,
  weight_label = "Euclidean distance",
  centrality_label = "Strength centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/strength/dist_stre.png", 
  plot = plot_dist_stre, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_dist_pop_gravity_stre <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = dist_pop_gravity, 
  centrality = dist_pop_gravity_stre,
  weight_label = "Gravity model with population and distance",
  centrality_label = "Strength centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/strength/dist_pop_gravity_stre.png", 
  plot = plot_dist_pop_gravity_stre, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_dist_forest_gravity_stre <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = dist_forest_gravity, 
  centrality = dist_forest_gravity_stre,
  weight_label = "Gravity model with forest loss and distance",
  centrality_label = "Strength centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/strength/dist_forest_gravity_stre.png", 
  plot = plot_dist_forest_gravity_stre, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_time_stre <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = time, 
  centrality = time_stre,
  weight_label = "Travel time",
  centrality_label = "Strength centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/strength/time_stre.png", 
  plot = plot_time_stre, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_time_pop_gravity_stre <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = time_pop_gravity, 
  centrality = time_pop_gravity_stre,
  weight_label = "Gravity model with population and travel time",
  centrality_label = "Strength centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/strength/time_pop_gravity_stre.png", 
  plot = plot_time_pop_gravity_stre, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_time_forest_gravity_stre <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = time_forest_gravity, 
  centrality = time_forest_gravity_stre,
  weight_label = "Gravity model with forest loss and travel time",
  centrality_label = "Strength centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/strength/time_forest_gravity_stre.png", 
  plot = plot_time_forest_gravity_stre, 
  height = 10, width = 9.5, units = "in"
)
```

### Closeness

```{r}
plot_dist_close <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = distance, 
  centrality = dist_close,
  weight_label = "Euclidean distance",
  centrality_label = "Closeness centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/closeness/dist_close.png", 
  plot = plot_dist_close, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_dist_pop_gravity_close <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = dist_pop_gravity, 
  centrality = dist_pop_gravity_close,
  weight_label = "Gravity model with population and distance",
  centrality_label = "Closeness centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/closeness/dist_pop_gravity_close.png", 
  plot = plot_dist_pop_gravity_close, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_dist_forest_gravity_close <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = dist_forest_gravity, 
  centrality = dist_forest_gravity_close,
  weight_label = "Gravity model with forest loss and distance",
  centrality_label = "Closeness centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/closeness/dist_forest_gravity_close.png", 
  plot = plot_dist_forest_gravity_close, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_time_close <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = time, 
  centrality = time_close,
  weight_label = "Travel time",
  centrality_label = "Closeness centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/closeness/time_close.png", 
  plot = plot_time_close, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_time_pop_gravity_close <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = time_pop_gravity, 
  centrality = time_pop_gravity_close,
  weight_label = "Gravity model with population and travel time",
  centrality_label = "Closeness centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/closeness/time_pop_gravity_close.png", 
  plot = plot_time_pop_gravity_close, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_time_forest_gravity_close <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = time_forest_gravity, 
  centrality = time_forest_gravity_close,
  weight_label = "Gravity model with forest loss and travel time",
  centrality_label = "Closeness centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/closeness/time_forest_gravity_close.png", 
  plot = plot_time_forest_gravity_close, 
  height = 10, width = 9.5, units = "in"
)
```

### Betweenness

```{r}
plot_dist_between <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = distance, 
  centrality = dist_between,
  weight_label = "Euclidean distance",
  centrality_label = "Betweenness centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/betweenness/dist_between.png", 
  plot = plot_dist_between, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_dist_pop_gravity_between <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = dist_pop_gravity, 
  centrality = dist_pop_gravity_between,
  weight_label = "Gravity model with population and distance",
  centrality_label = "Betweenness centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/betweenness/dist_pop_gravity_between.png", 
  plot = plot_dist_pop_gravity_between, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_dist_forest_gravity_between <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = dist_forest_gravity, 
  centrality = dist_forest_gravity_between,
  weight_label = "Gravity model with forest loss and distance",
  centrality_label = "Betweenness centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/betweenness/dist_forest_gravity_between.png", 
  plot = plot_dist_forest_gravity_between, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_time_between <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = time, 
  centrality = time_between,
  weight_label = "Travel time",
  centrality_label = "Betweenness centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/betweenness/time_between.png", 
  plot = plot_time_between, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_time_pop_gravity_between <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = time_pop_gravity, 
  centrality = time_pop_gravity_between,
  weight_label = "Gravity model with population and travel time",
  centrality_label = "Betweenness centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/betweenness/time_pop_gravity_between.png", 
  plot = plot_time_pop_gravity_between, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_time_forest_gravity_between <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = time_forest_gravity, 
  centrality = time_forest_gravity_between,
  weight_label = "Gravity model with forest loss and travel time",
  centrality_label = "Betweenness centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/betweenness/time_forest_gravity_between.png", 
  plot = plot_time_forest_gravity_between, 
  height = 10, width = 9.5, units = "in"
)
```

### Eigenvector

```{r}
plot_dist_eigen <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = distance, 
  centrality = dist_eigen,
  weight_label = "Euclidean distance",
  centrality_label = "Eigenvector centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/eigenvector/dist_eigen.png", 
  plot = plot_dist_eigen, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_dist_pop_gravity_eigen <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = dist_pop_gravity, 
  centrality = dist_pop_gravity_eigen,
  weight_label = "Gravity model with population and distance",
  centrality_label = "Eigenvector centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/eigenvector/dist_pop_gravity_eigen.png", 
  plot = plot_dist_pop_gravity_eigen, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_dist_forest_gravity_eigen <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = dist_forest_gravity, 
  centrality = dist_forest_gravity_eigen,
  weight_label = "Gravity model with forest loss and distance",
  centrality_label = "Eigenvector centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/eigenvector/dist_forest_gravity_eigen.png", 
  plot = plot_dist_forest_gravity_eigen, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_time_eigen <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = time, 
  centrality = time_eigen,
  weight_label = "Travel time",
  centrality_label = "Eigenvector centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/eigenvector/time_eigen.png", 
  plot = plot_time_eigen, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_time_pop_gravity_eigen <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = time_pop_gravity, 
  centrality = time_pop_gravity_eigen,
  weight_label = "Gravity model with population and travel time",
  centrality_label = "Eigenvector centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/eigenvector/time_pop_gravity_eigen.png", 
  plot = plot_time_pop_gravity_eigen, 
  height = 10, width = 9.5, units = "in"
)
```

```{r}
plot_time_forest_gravity_eigen <- plot_single_graph(
  graphs = graphs, layout = layout_kk,  
  weight = time_forest_gravity, 
  centrality = time_forest_gravity_eigen,
  weight_label = "Gravity model with forest loss and travel time",
  centrality_label = "Eigenvector centrality",
  group = watershed
)

ggsave(
  "figures/single-graph/eigenvector/time_forest_gravity_eigen.png", 
  plot = plot_time_forest_gravity_eigen, 
  height = 10, width = 9.5, units = "in"
)
```

## Multiple graphs

```{r}
source("R/plot_multiple_graph.R")
```

```{r}
layout_stress <- create_layout(graphs, layout = "stress")
```

### Strength

```{r}
facet_dist_stre <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = distance, centrality = dist_stre, 
  weight_label = "Euclidean distance",
  centrality_label = "Strength centrality"
)

ggsave(
  "figures/multiple-graph/strength/facet_dist_stre.png", 
  plot = facet_dist_stre, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_dist_pop_gravity_stre <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = dist_pop_gravity, centrality = dist_pop_gravity_stre, 
  weight_label = "Gravity model weight with population and distance",
  centrality_label = "Strength centrality"
)

ggsave(
  "figures/multiple-graph/strength/dist_pop_gravity_stre.png", 
  plot = facet_dist_pop_gravity_stre, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_dist_forest_gravity_stre <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = dist_forest_gravity, centrality = dist_forest_gravity_stre, 
  weight_label = "Gravity model weight with forest loss and distance",
  centrality_label = "Strength centrality"
)

ggsave(
  "figures/multiple-graph/strength/dist_forest_gravity_stre.png", 
  plot = facet_dist_forest_gravity_stre, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_time_stre <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = distance, centrality = time_stre, 
  weight_label = "Travel time",
  centrality_label = "Strength centrality"
)

ggsave(
  "figures/multiple-graph/strength/facet_time_stre.png", 
  plot = facet_time_stre, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_time_pop_gravity_stre <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = time_pop_gravity, centrality = time_pop_gravity_stre, 
  weight_label = "Gravity model weight with population and travel time",
  centrality_label = "Strength centrality"
)

ggsave(
  "figures/multiple-graph/strength/time_pop_gravity_stre.png", 
  plot = facet_time_pop_gravity_stre, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_time_forest_gravity_stre <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = time_forest_gravity, centrality = time_forest_gravity_stre, 
  weight_label = "Gravity model weight with forest loss and travel time",
  centrality_label = "Strength centrality"
)

ggsave(
  "figures/multiple-graph/strength/time_forest_gravity_stre.png", 
  plot = facet_time_forest_gravity_stre, 
  height = 20, width = 14, units = "in"
)
```

### Closeness

```{r}
facet_dist_close <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = distance, centrality = dist_close, 
  weight_label = "Euclidean distance",
  centrality_label = "Closeness centrality"
)

ggsave(
  "figures/multiple-graph/closeness/facet_dist_close.png", 
  plot = facet_dist_close, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_dist_pop_gravity_close <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = dist_pop_gravity, centrality = dist_pop_gravity_close, 
  weight_label = "Gravity model weight with population and distance",
  centrality_label = "Closeness centrality"
)

ggsave(
  "figures/multiple-graph/closeness/dist_pop_gravity_close.png", 
  plot = facet_dist_pop_gravity_close, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_dist_forest_gravity_close <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = dist_forest_gravity, centrality = dist_forest_gravity_close, 
  weight_label = "Gravity model weight with forest loss and distance",
  centrality_label = "Closeness centrality"
)

ggsave(
  "figures/multiple-graph/closeness/dist_forest_gravity_close.png", 
  plot = facet_dist_forest_gravity_close, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_time_close <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = distance, centrality = time_close, 
  weight_label = "Travel time",
  centrality_label = "Closeness centrality"
)

ggsave(
  "figures/multiple-graph/closeness/facet_time_close.png", 
  plot = facet_time_close, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_time_pop_gravity_close <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = time_pop_gravity, centrality = time_pop_gravity_close, 
  weight_label = "Gravity model weight with population and travel time",
  centrality_label = "Closeness centrality"
)

ggsave(
  "figures/multiple-graph/closeness/time_pop_gravity_close.png", 
  plot = facet_time_pop_gravity_close, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_time_forest_gravity_close <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = time_forest_gravity, centrality = time_forest_gravity_close, 
  weight_label = "Gravity model weight with forest loss and travel time",
  centrality_label = "Closeness centrality"
)

ggsave(
  "figures/multiple-graph/closeness/time_forest_gravity_close.png", 
  plot = facet_time_forest_gravity_close, 
  height = 20, width = 14, units = "in"
)
```

### Betweenness

```{r}
facet_dist_between <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = distance, centrality = dist_between, 
  weight_label = "Euclidean distance",
  centrality_label = "Betweenness centrality"
)

ggsave(
  "figures/multiple-graph/betweenness/facet_dist_between.png", 
  plot = facet_dist_between, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_dist_pop_gravity_between <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = dist_pop_gravity, centrality = dist_pop_gravity_between, 
  weight_label = "Gravity model weight with population and distance",
  centrality_label = "Betweenness centrality"
)

ggsave(
  "figures/multiple-graph/betweenness/dist_pop_gravity_between.png", 
  plot = facet_dist_pop_gravity_between, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_dist_forest_gravity_between <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = dist_forest_gravity, centrality = dist_forest_gravity_between, 
  weight_label = "Gravity model weight with forest loss and distance",
  centrality_label = "Betweenness centrality"
)

ggsave(
  "figures/multiple-graph/betweenness/dist_forest_gravity_between.png", 
  plot = facet_dist_forest_gravity_between, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_time_between <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = distance, centrality = time_between, 
  weight_label = "Travel time",
  centrality_label = "Betweenness centrality"
)

ggsave(
  "figures/multiple-graph/betweenness/facet_time_between.png", 
  plot = facet_time_between, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_time_pop_gravity_between <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = time_pop_gravity, centrality = time_pop_gravity_between, 
  weight_label = "Gravity model weight with population and travel time",
  centrality_label = "Betweenness centrality"
)

ggsave(
  "figures/multiple-graph/betweenness/time_pop_gravity_between.png", 
  plot = facet_time_pop_gravity_between, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_time_forest_gravity_between <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = time_forest_gravity, centrality = time_forest_gravity_between, 
  weight_label = "Gravity model weight with forest loss and travel time",
  centrality_label = "Betweenness centrality"
)

ggsave(
  "figures/multiple-graph/betweenness/time_forest_gravity_between.png", 
  plot = facet_time_forest_gravity_between, 
  height = 20, width = 14, units = "in"
)
```

### Eigenvector

```{r}
facet_dist_eigen <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = distance, centrality = dist_eigen, 
  weight_label = "Euclidean distance",
  centrality_label = "Eigenvector centrality"
)

ggsave(
  "figures/multiple-graph/eigenvector/facet_dist_eigen.png", 
  plot = facet_dist_eigen, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_dist_pop_gravity_eigen <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = dist_pop_gravity, centrality = dist_pop_gravity_eigen, 
  weight_label = "Gravity model weight with population and distance",
  centrality_label = "Eigenvector centrality"
)

ggsave(
  "figures/multiple-graph/eigenvector/dist_pop_gravity_eigen.png", 
  plot = facet_dist_pop_gravity_eigen, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_dist_forest_gravity_eigen <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = dist_forest_gravity, centrality = dist_forest_gravity_eigen, 
  weight_label = "Gravity model weight with forest loss and distance",
  centrality_label = "Eigenvector centrality"
)

ggsave(
  "figures/multiple-graph/eigenvector/dist_forest_gravity_eigen.png", 
  plot = facet_dist_forest_gravity_eigen, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_time_eigen <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = distance, centrality = time_eigen, 
  weight_label = "Travel time",
  centrality_label = "Eigenvector centrality"
)

ggsave(
  "figures/multiple-graph/eigenvector/facet_time_eigen.png", 
  plot = facet_time_eigen, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_time_pop_gravity_eigen <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = time_pop_gravity, centrality = time_pop_gravity_eigen, 
  weight_label = "Gravity model weight with population and travel time",
  centrality_label = "Eigenvector centrality"
)

ggsave(
  "figures/multiple-graph/eigenvector/time_pop_gravity_eigen.png", 
  plot = facet_time_pop_gravity_eigen, 
  height = 20, width = 14, units = "in"
)
```

```{r}
facet_time_forest_gravity_eigen <- plot_multiple_graph(
  graphs = graphs, layout = layout_stress, facet = watershed,
  weight = time_forest_gravity, centrality = time_forest_gravity_eigen, 
  weight_label = "Gravity model weight with forest loss and travel time",
  centrality_label = "Eigenvector centrality"
)

ggsave(
  "figures/multiple-graph/eigenvector/time_forest_gravity_eigen.png", 
  plot = facet_time_forest_gravity_eigen, 
  height = 20, width = 14, units = "in"
)
```

## 3D graphs

```{r}
nodes_layout_kk <- igraph::layout_with_kk(
  graphs, dim = 3
)

# layout_3d = igraph::layout_with_fr(
#   graphs, dim = 3, minx = rep(-30, nrow(nodes)), maxx = rep(30, nrow(nodes)), 
#   miny = rep(-30, nrow(nodes)), maxy = rep(30, nrow(nodes)), 
#   minz = rep(-30, nrow(nodes)), maxz = rep(30, nrow(nodes))
# )

nodes_layout_drl <- igraph::layout_with_drl(
  graphs, dim = 3
)

# layout_3d = graphlayouts::layout_with_stress3D(
#   graphs
# )

colnames(nodes_layout_kk) <- c("x", "y", "z")
colnames(nodes_layout_drl) <- c("x", "y", "z")

nodes_layout_kk_attr <- bind_cols(data.frame(nodes_layout_kk), village_centralities)
nodes_layout_drl_attr <- bind_cols(data.frame(nodes_layout_drl), village_centralities)
```

```{r}
edges_layout_kk_x <- rep(0, nrow(edges) * 3)
edges_layout_kk_y <- rep(0, nrow(edges) * 3)
edges_layout_kk_z <- rep(0, nrow(edges) * 3)

edges_layout_id <- rep(0, nrow(edges) * 3)

edges_layout_drl_x <- rep(0, nrow(edges) * 3)
edges_layout_drl_y <- rep(0, nrow(edges) * 3)
edges_layout_drl_z <- rep(0, nrow(edges) * 3)

for (row in 1:nrow(edges)) {
  from = edges$from[row]
  to = edges$to[row]
  x_from = nodes_layout_kk_attr$x[nodes_layout_kk_attr$village_id == from]
  x_to = nodes_layout_kk_attr$x[nodes_layout_kk_attr$village_id == to]
  y_from = nodes_layout_kk_attr$y[nodes_layout_kk_attr$village_id == from]
  y_to = nodes_layout_kk_attr$y[nodes_layout_kk_attr$village_id == to]
  z_from = nodes_layout_kk_attr$z[nodes_layout_kk_attr$village_id == from]
  z_to = nodes_layout_kk_attr$z[nodes_layout_kk_attr$village_id == to]
  edges_layout_kk_x[((row*3)-2):(row*3)] = c(x_from, x_to, NA)
  edges_layout_kk_y[((row*3)-2):(row*3)] = c(y_from, y_to, NA)
  edges_layout_kk_z[((row*3)-2):(row*3)] = c(z_from, z_to, NA)
  edges_layout_id[((row*3)-2):(row*3)] = c(from, to, NA)
}

for (row in 1:nrow(edges)) {
  from = edges$from[row]
  to = edges$to[row]
  x_from = nodes_layout_drl_attr$x[nodes_layout_drl_attr$village_id == from]
  x_to = nodes_layout_drl_attr$x[nodes_layout_drl_attr$village_id == to]
  y_from = nodes_layout_drl_attr$y[nodes_layout_drl_attr$village_id == from]
  y_to = nodes_layout_drl_attr$y[nodes_layout_drl_attr$village_id == to]
  z_from = nodes_layout_drl_attr$z[nodes_layout_drl_attr$village_id == from]
  z_to = nodes_layout_drl_attr$z[nodes_layout_kk_attr$village_id == to]
  edges_layout_drl_x[((row*3)-2):(row*3)] = c(x_from, x_to, NA)
  edges_layout_drl_y[((row*3)-2):(row*3)] = c(y_from, y_to, NA)
  edges_layout_drl_z[((row*3)-2):(row*3)] = c(z_from, z_to, NA)
}

edges_layout_kk <- tibble(
  x = edges_layout_kk_x, y = edges_layout_kk_y, z = edges_layout_kk_z,
  `village_id` = edges_layout_id
)

edges_layout_drl <- tibble(
  x = edges_layout_drl_x, y = edges_layout_drl_y, z = edges_layout_drl_z,
  `village_id` = edges_layout_id
)

edges_layout_kk_attr <- edges_layout_kk |> 
  left_join(nodes, by = "village_id")
edges_layout_drl_attr <- edges_layout_drl |> 
  left_join(nodes, by = "village_id")
```

```{r}
library(plotly)

fig <- plot_ly()

fig_edges_kk <- fig |>
  add_trace(
    data = edges_layout_kk_attr, x = ~x, y = ~y, z = ~z, type = "scatter3d", 
    mode = "lines", color = ~watershed, showlegend = FALSE,
    line = list(width = 0.5, opacity = 0.5), opacity = 0.5, hoverinfo = "none"
  )

fig_network_kk <- fig_edges_kk |>
  add_trace(
  data = nodes_layout_kk_attr, type = "scatter3d", mode = "markers", x = ~x, 
  y = ~y, z = ~z, size = ~population, color = ~watershed,
  marker = list(symbol = "circle", sizemode = "diameter", opacity = 1), 
  sizes = c(5, 20), 
  opacity = 1, showlegend = FALSE, hoverinfo = "text", 
  text = ~paste(
    "Watershed:", watershed, "<br>Village:", village, "<br>Population:", population
  )
  )

axis <- list(
  showbackground = FALSE, showline = FALSE, zeroline = FALSE, showgrid = FALSE,
  showticklabels = FALSE, title = ""
)

fig_3d_network_kk <- fig_network_kk |> 
  layout(
    title = "", scene = list(xaxis = axis, yaxis = axis, zaxis = axis),
    margin = list(b = 40, l = 40, r = 40, t = 40)
  )

htmlwidgets::saveWidget(
  fig_3d_network_kk, file = "figures/3d-network-kk.html", selfcontained = TRUE
)
```

```{r}
fig_edges_drl <- fig |>
  add_trace(
    data = edges_layout_drl_attr, x = ~x, y = ~y, z = ~z, type = "scatter3d", 
    mode = "lines", color = ~watershed, showlegend = FALSE,
    line = list(width = 0.5, opacity = 0.5), opacity = 0.5, hoverinfo = "none"
  )

fig_network_drl <- fig_edges_drl |>
  add_trace(
  data = nodes_layout_drl_attr, type = "scatter3d", mode = "markers", x = ~x, 
  y = ~y, z = ~z, size = ~population, color = ~watershed,
  marker = list(symbol = "circle", sizemode = "diameter", opacity = 1), 
  sizes = c(5, 20), 
  opacity = 1, showlegend = FALSE, hoverinfo = "text", 
  text = ~paste(
    "Watershed:", watershed, "<br>Village:", village, "<br>Population:", population
  )
  )

fig_3d_network_drl <- fig_network_drl |> 
  layout(
    title = "", scene = list(xaxis = axis, yaxis = axis, zaxis = axis),
    margin = list(b = 40, l = 40, r = 40, t = 40)
  )

htmlwidgets::saveWidget(
  fig_3d_network_drl, file = "figures/3d-network-drl.html", selfcontained = TRUE
)
```