Reproduction of the Analytical Steps of WEFSNexus Index

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

plot(cars)

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.

# --- 1. LOAD LIBRARIES ---
library(sf)
library(tidyverse)
library(spdep)
library(stringr)
library(broom)

# --- 2. PREPARE NUMERIC DATA (Excel Results) ---
composite_trends_ols <- composite_index_score %>%
  mutate(Year = as.numeric(as.character(Year))) %>%
  group_by(District) %>%
  do(tidy(lm(`Composite Index Score` ~ Year, data = .))) %>%
  filter(term == "Year") %>%
  mutate(District = str_replace_all(District, "\\|", "i"), # Fixes Pan|pat
         District = str_trim(District),
         CleanName = toupper(District)) %>%
  mutate(CleanName = case_when(
    CleanName == "GURGAON" ~ "GURUGRAM",
    CleanName == "MEWAT"   ~ "NUH",
    CleanName == "CHARKHIDADRI" ~ "CHARKHI DADRI", # Fixes the no-space version
    TRUE ~ CleanName
  )) %>%
  dplyr::select(CleanName, beta_ols = estimate)

# --- 3. PREPARE SPATIAL DATA (KML Map) ---
haryana_map_clean <- st_read("~/Desktop/ec8e264b-aadd-4adf-b349-5d8a5753a22d.kml") %>%
  st_zm() %>% 
  st_make_valid() %>%
  mutate(Name = str_replace_all(Name, "\\|", "i"),
         Name = str_trim(Name),
         CleanName = toupper(Name)) %>%
  mutate(CleanName = case_when(
    CleanName == "GURGAON" ~ "GURUGRAM",
    CleanName == "MEWAT"   ~ "NUH",
    CleanName == "CHARKI DADRI" ~ "CHARKHI DADRI",
    TRUE ~ CleanName
  ))

Reading layer `Haryana Districts' from data source 
  `/Users/zhangchuanhong/Desktop/ec8e264b-aadd-4adf-b349-5d8a5753a22d.kml' 
  using driver `KML'
Simple feature collection with 22 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY, XYZ
Bounding box:  xmin: 74.47326 ymin: 27.65227 xmax: 77.60459 ymax: 30.92854
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84

# --- 4. JOIN (Target: 22 Districts) ---
haryana_spatial <- haryana_map_clean %>%
  inner_join(composite_trends_ols, by = "CleanName")

print(paste("District Count for Analysis:", nrow(haryana_spatial))) 

[1] "District Count for Analysis: 22"

# --- 5. GLOBAL MORAN'S I (Power Test) ---
nb <- poly2nb(haryana_spatial, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

set.seed(123)
moran_mc <- moran.mc(haryana_spatial$beta_ols, lw, nsim = 999, zero.policy = TRUE)
print(moran_mc)

    Monte-Carlo simulation of Moran I

data:  haryana_spatial$beta_ols 
weights: lw  
number of simulations + 1: 1000 

statistic = 0.26042, observed rank = 973, p-value = 0.027
alternative hypothesis: greater

# LISA Cluster of Index
library(sf)
library(tidyverse)
library(spdep)
library(gridExtra)
# 1. Prepare Data: Filter for the 3 key years and ensure naming matches
target_years <- c(2018, 2022, 2025)
map_data_years <- composite_index_score %>%
filter(Year %in% target_years) %>%
mutate(CleanName = toupper(str_trim(District))) %>%
mutate(CleanName = case_when(
CleanName == "GURGAON" ~ "GURUGRAM",
CleanName == "MEWAT"   ~ "NUH",
CleanName == "CHARKHIDADRI" ~ "CHARKHI DADRI",
TRUE ~ CleanName
))
# 2. Join with the Haryana SF map object
haryana_lisa <- haryana_map %>%
left_join(map_data_years, by = "CleanName")
# 3. Create a function to calculate LISA for a specific year
get_lisa_for_year <- function(yr) {
# Filter for the year
temp_sf <- haryana_lisa %>% filter(Year == yr)
# Weights matrix
nb <- poly2nb(temp_sf, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
# Local Moran
local_m <- localmoran(temp_sf$`Composite Index Score`, lw, zero.policy = TRUE)
# Manual Quadrant calculation
scale_val <- scale(temp_sf$`Composite Index Score`)
lag_val <- lag.listw(lw, scale_val)
p_vals <- local_m[, 5]
temp_sf$Category <- "Not Significant"
# Using p <= 0.1 for small N robustness
temp_sf$Category[scale_val > 0 & lag_val > 0 & p_vals <= 0.1] <- "High-High (Hotspot)"
temp_sf$Category[scale_val < 0 & lag_val < 0 & p_vals <= 0.1] <- "Low-Low (Coldspot)"
temp_sf$Category[scale_val < 0 & lag_val > 0 & p_vals <= 0.1] <- "Low-High (Outlier)"
temp_sf$Category[scale_val > 0 & lag_val < 0 & p_vals <= 0.1] <- "High-Low (Outlier)"
return(temp_sf)
}
# 4. Run the function for all 3 years and combine
lisa_2018 <- get_lisa_for_year(2018)
lisa_2022 <- get_lisa_for_year(2022)
lisa_2025 <- get_lisa_for_year(2025)
final_lisa_data <- rbind(lisa_2018, lisa_2022, lisa_2025)
# 5. Create the Side-by-Side Plot
ggplot(data = final_lisa_data) +
geom_sf(aes(fill = Category), color = "black", size = 0.2) +
facet_wrap(~Year, ncol = 3) +
scale_fill_manual(values = c(
"High-High (Hotspot)" = "#d73027", # Red
"Low-Low (Coldspot)"  = "#4575b4", # Blue
"Low-High (Outlier)"  = "#74add1", # Light Blue
"High-Low (Outlier)"  = "#fdae61", # Orange
"Not Significant"     = "#f0f0f0"  # Grey
)) +
theme_minimal() +
theme(
legend.position = "bottom",
axis.text = element_blank(),
panel.grid = element_blank(),
strip.text = element_text(face = "bold", size = 14)
) +
labs(title = "Spatial Dynamic Evolution of Nexus Sustainability",
subtitle = "Local Indicators of Spatial Association (LISA) for 2018, 2022, and 2025",
fill = "LISA Cluster")
# 6. Save for Thesis
ggsave("Triple_LISA_Map.png", width = 15, height = 7, dpi = 300)

# LISA Cluster Evolution of Specific Pillars
# 1. Load Libraries
library(sf)
library(tidyverse)
library(spdep)
library(stringr)

# 2. Prepare the Data
pillar_data_long <- four_pillar %>%
  pivot_longer(cols = c(`Water Pillar`, `Energy Pillar`, `Food Pillar`, `Social Pillar`), 
               names_to = "Pillar", values_to = "Score") %>%
  mutate(CleanName = toupper(str_trim(District)))

# 3. Create Spatial Dataframe
final_lisa_grid <- data.frame()
pillars <- unique(pillar_data_long$Pillar)
years <- sort(unique(pillar_data_long$Year))

for (p in pillars) {
  for (y in years) {
    
    # Subset data for specific pillar and year
    subset_data <- pillar_data_long %>% filter(Pillar == p, Year == y)
    
    # Join with your cleaned haryana_map object
    temp_sf <- haryana_map %>%
      inner_join(subset_data, by = "CleanName")
    
    # Calculate Spatial Weights (Queen contiguity)
    nb <- poly2nb(temp_sf, queen = TRUE)
    lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
    
    # Calculate Local Moran's I
    local_m <- localmoran(temp_sf$Score, lw, zero.policy = TRUE)
    
    # Manual Quadrant Calculation
    scale_val <- as.numeric(scale(temp_sf$Score))
    lag_val <- as.numeric(lag.listw(lw, scale_val, zero.policy = TRUE))
    p_vals <- local_m[, 5]
    
    # Classify the clusters (using p <= 0.1 for small N robustness)
    temp_sf$LISA_Cat <- "Not Significant"
    temp_sf$LISA_Cat[scale_val > 0 & lag_val > 0 & p_vals <= 0.1] <- "High-High (Hotspot)"
    temp_sf$LISA_Cat[scale_val < 0 & lag_val < 0 & p_vals <= 0.1] <- "Low-Low (Coldspot)"
    temp_sf$LISA_Cat[scale_val < 0 & lag_val > 0 & p_vals <= 0.1] <- "Low-High (Outlier)"
    temp_sf$LISA_Cat[scale_val > 0 & lag_val < 0 & p_vals <= 0.1] <- "High-Low (Outlier)"
    
    # Collect results
    final_lisa_grid <- rbind(final_lisa_grid, temp_sf)
  }
}

# 4. Visualization
ggplot(data = final_lisa_grid) +
  geom_sf(aes(fill = LISA_Cat), color = "black", size = 0.05) +
  facet_grid(Pillar ~ Year) +
  scale_fill_manual(values = c(
    "High-High (Hotspot)" = "#d73027", # Red
    "Low-Low (Coldspot)"  = "#4575b4", # Blue
    "Low-High (Outlier)"  = "#74add1", # Light Blue
    "High-Low (Outlier)"  = "#fdae61", # Orange
    "Not Significant"     = "#f0f0f0"  # Grey
  )) +
  theme_minimal() +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank(),
    strip.text.x = element_text(face = "bold", size = 10),
    strip.text.y = element_text(face = "bold", size = 10, angle = 0),
    legend.position = "bottom",
    legend.title = element_text(face = "bold")
  ) +
  labs(title = "LISA Cluster Evolution of Specific Pillars (2018-2025)",
       subtitle = "4x8 Matrix of Local Indicators of Spatial Association (LISA)",
       fill = "LISA Cluster Type")

# 5. Save
ggsave("Nexus_LISA_Matrix_4x8.png", width = 18, height = 12, dpi = 300)

# Transition Velocity LISA Clusters of Specific Pillarslibrary(sf)
library(tidyverse)
library(spdep)
library(broom)

# 1. Calculate the OLS Betas for EACH Pillar separately
pillar_betas <- four_pillar_long %>%
mutate(Year = as.numeric(as.character(Year))) %>%
group_by(District, Pillar) %>%
do(tidy(lm(Score ~ Year, data = .))) %>%
filter(term == "Year") %>%
dplyr::select(District, Pillar, beta = estimate) %>%
mutate(CleanName = toupper(str_trim(District))) %>%
mutate(CleanName = case_when(
CleanName == "GURGAON" ~ "GURUGRAM",
CleanName == "MEWAT"   ~ "NUH",
CleanName == "CHARKHIDADRI" ~ "CHARKHI DADRI",
TRUE ~ CleanName
))

# 2. JOIN and Run Moran's I Loop
pillars <- unique(pillar_betas$Pillar)
moran_summary <- data.frame()
for (p in pillars) {
# Subset for current pillar and join with map
p_data <- pillar_betas %>% filter(Pillar == p)
p_spatial <- haryana_map %>% inner_join(p_data, by = "CleanName")
# Weights Matrix
nb <- poly2nb(p_spatial, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
# Global Moran's I (Monte Carlo)
set.seed(123)
mc_res <- moran.mc(p_spatial$beta, lw, nsim = 999, zero.policy = TRUE)
# Store Result
moran_summary <- rbind(moran_summary, data.frame(
Pillar = p,
Moran_I = mc_res$statistic,
p_value = mc_res$p.value
))
}
print(moran_summary)

                  Pillar     Moran_I p_value
statistic  Energy Pillar  0.59051342   0.001
statistic1   Food Pillar -0.08374114   0.621
statistic2 Social Pillar  0.49570507   0.001
statistic3  Water Pillar  0.37901555   0.004

# 3. Cleanup the Beta Data (Regression Results)
pillar_betas_clean <- pillar_betas %>%
mutate(CleanName = District %>%
str_replace_all("\\|", "i") %>%
str_trim() %>%
toupper()) %>%
mutate(CleanName = case_when(
CleanName == "GURGAON" ~ "GURUGRAM",
CleanName == "MEWAT"   ~ "NUH",
CleanName == "CHARKHIDADRI" ~ "CHARKHI DADRI",
TRUE ~ CleanName
))

# 5. Initialize the loop for spatial analysis
pillars_to_plot <- unique(pillar_betas_clean$Pillar)
lisa_output_list <- list()
for (p in pillars_to_plot) {
# Filter for specific pillar and join with the spatial haryana_map
p_data <- pillar_betas_clean %>% filter(Pillar == p)
# Ensure haryana_map (your KML) is available and has 'CleanName'
p_spatial <- haryana_map %>%
mutate(CleanName = Name %>% str_replace_all("\\|", "i") %>% str_trim() %>% toupper()) %>%
mutate(CleanName = case_when(
CleanName == "GURGAON" ~ "GURUGRAM",
CleanName == "MEWAT"   ~ "NUH",
CleanName == "CHARKI DADRI" ~ "CHARKHI DADRI",
TRUE ~ CleanName
)) %>%
inner_join(p_data, by = "CleanName")
# Calculate Neighbors and Weights
nb <- poly2nb(p_spatial, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
# Local Moran's I math
local_m <- localmoran(p_spatial$beta, lw, zero.policy = TRUE)
# Center the data around the mean
diff_beta <- p_spatial$beta - mean(p_spatial$beta)
# Calculate the spatial lag (what the neighbors are doing)
lag_beta <- lag.listw(lw, diff_beta, zero.policy = TRUE)
# Get p-values
p_vals <- local_m[, 5]
p_spatial$LISA_Cat <- "Not Significant"
# p <= 0.1 for small sample robustness
p_spatial$LISA_Cat[diff_beta > 0 & lag_beta > 0 & p_vals <= 0.1] <- "High-High (Hotspot)"
p_spatial$LISA_Cat[diff_beta < 0 & lag_beta < 0 & p_vals <= 0.1] <- "Low-Low (Coldspot)"
p_spatial$LISA_Cat[diff_beta < 0 & lag_beta > 0 & p_vals <= 0.1] <- "Low-High (Outlier)"
p_spatial$LISA_Cat[diff_beta > 0 & lag_beta < 0 & p_vals <= 0.1] <- "High-Low (Outlier)"
lisa_output_list[[p]] <- p_spatial
}
# Combine all pillars for the final plot
all_pillars_lisa <- do.call(rbind, lisa_output_list)

# 6. Visualization
ggplot(data = all_pillars_lisa) +
geom_sf(aes(fill = LISA_Cat), color = "black", size = 0.2) +
facet_wrap(~Pillar, ncol = 2) +
scale_fill_manual(values = c(
"High-High (Hotspot)" = "#d73027", # Red
"Low-Low (Coldspot)"  = "#4575b4", # Blue
"Low-High (Outlier)"  = "#74add1", # Light Blue
"High-Low (Outlier)"  = "#fdae61", # Orange
"Not Significant"     = "#f0f0f0"  # Grey
), name = "Momentum Type") +
theme_minimal() +
theme(
strip.text = element_text(face = "bold", size = 12),
legend.position = "bottom",
axis.text = element_blank(),
panel.grid = element_blank()
) +
labs(title = "Regional Momentum of the Agroecological Transition",
subtitle = "Local Moran's I Analysis of Pillar Growth Rates (2018-2025)")

# 7. Save
ggsave("Transition_Momentum_LISA_Final.png", width = 12, height = 10, dpi = 300)

# Link to other Indicators
library(readxl)
View(SDG)
library(tidyverse)
library(ggpubr) # For easy correlation plots

# 1. Prepare your data for the year 2024
# Filtering and selecting relevant columns from your dataframes
my_comp_2024 <- composite_index_score %>%
filter(Year == 2024) %>%
select(District, `Composite Index Score`)
my_pillars_2024 <- four_pillar %>%
filter(Year == 2024) %>%
select(District, `Water Pillar`, `Energy Pillar`, `Food Pillar`, `Social Pillar`)

# 2. Join with the SDG dataframe
comparison_df <- my_comp_2024 %>%
left_join(my_pillars_2024, by = "District") %>%
left_join(SDG, by = "District")

# --- VISUAL 1: Composite Comparison ---
# Added na.rm = TRUE inside ggscatter so geom_text_repel inherits it
plot1 <- ggscatter(comparison_df, x = "Composite Index Score", y = "SDG_Composite",
          add = "reg.line", conf.int = TRUE,
          cor.coef = TRUE, cor.method = "spearman",
          xlab = "XXXX Composite Index (2024)", ylab = "Haryana SDG Composite Score (2024)",
          title = "Validation: XXXX Index vs. Official SDG Index",
          label = "District", repel = TRUE, font.label = c(8, "plain"),
          na.rm = TRUE) + 
  theme_minimal()

# --- VISUAL 2: Pillar-Specific Validation (Faceted) ---
# Added na.rm = TRUE to geom_point, geom_smooth, AND stat_cor separately
plot2 <- ggplot(pillar_comparison, aes(x = My_Score, y = SDG_Score)) +
  geom_point(color = "steelblue", na.rm = TRUE) +
  geom_smooth(method = "lm", color = "darkred", fill = "lightgray", na.rm = TRUE) +
  stat_cor(method = "spearman", label.x.npc = "left", label.y.npc = "top", na.rm = TRUE) + 
  facet_wrap(~Pillar, scales = "free") +
  labs(title = "Pillar-Specific Validation against Corresponding SDGs",
       subtitle = "Social-SDG1, Food-SDG2, Water-SDG6, Energy-SDG7",
       x = "XXXX Pillar Scores", y = "Official SDG Scores") +
  theme_bw()

# Display Visuals without warnings popping up
print(plot1)

print(plot2)