# Dr.(C) Antonio Villalpando Acuña
# Atlas de Riesgos para la Nutrición de la Niñez en México de Save the Children
# Script de mapeo de riesgos

library(readxl)
library(writexl)
library(caret)
library(xgboost)
library(sf)
library(ggplot2)
library(tmap)
library(tmaptools)
library(leaflet)
library(plyr)
library(dplyr)
library(ggrepel)

# Cargar los datos
datos <- as.data.frame(read_xlsx("dc.xlsx"))
datos03 <- as.data.frame(read_xlsx("dc03.xlsx"))
datos46 <- as.data.frame(read_xlsx("dc46.xlsx"))
datos79 <- as.data.frame(read_xlsx("dc79.xlsx"))

### Procesar y mapear para el grupo de edad 0-3 años

# Selección de variables independientes y dependientes
x_03 <- datos03[, 9:146]
y_03 <- datos03$pesoedad

# Eliminar filas con valores NA en la variable dependiente o en las predictoras
data_complete_03 <- datos03[complete.cases(x_03, y_03), ]
x_complete_03 <- data_complete_03[, 9:143]
y_complete_03 <- data_complete_03$pesoedad

# Datos de entrenamiento y prueba
set.seed(123)  # Para reproducibilidad
trainIndex_03 <- createDataPartition(y_complete_03, p = .8, list = FALSE, times = 1)
x_train_03 <- x_complete_03[trainIndex_03, ]
x_test_03 <- x_complete_03[-trainIndex_03, ]
y_train_03 <- y_complete_03[trainIndex_03]
y_test_03 <- y_complete_03[-trainIndex_03]

# Convertir los datos a la matriz DMatrix de XGBoost
dtrain_03 <- xgb.DMatrix(data = as.matrix(x_train_03), label = y_train_03)
dtest_03 <- xgb.DMatrix(data = as.matrix(x_test_03), label = y_test_03)

# Entrenamiento del modelo XGBoost
params <- list(objective = "reg:squarederror", eval_metric = "rmse")
xgb_model_03 <- xgboost(data = dtrain_03, params = params, nrounds = 100, verbose = 0)

# Predicción con el modelo XGBoost
xgb_predictions_03 <- predict(xgb_model_03, newdata = dtest_03)

# Agregar predicciones al dataframe original
datos03$pesoedad_pred <- NA
datos03$pesoedad_pred[complete.cases(x_03, y_03)] <- predict(xgb_model_03, newdata = xgb.DMatrix(data = as.matrix(x_complete_03)))

# Calcular mediana y desviación estándar de las predicciones
median_pred_03 <- median(datos03$pesoedad_pred, na.rm = TRUE)
sd_pred_03 <- sd(datos03$pesoedad_pred, na.rm = TRUE)

# Calcular riesgo de desnutrición crónica (pesoedad < -2 desviaciones estándar)
threshold_desnutricion_03 <- median_pred_03 - 2 * sd_pred_03
datos03$prob_riesgo_desnutricion <- pnorm(threshold_desnutricion_03, mean = datos03$pesoedad_pred, sd = sd_pred_03)

# Calcular riesgo de obesidad (pesoedad > +2 desviaciones estándar)
threshold_obesidad_03 <- median_pred_03 + 1 * sd_pred_03
datos03$prob_riesgo_obesidad <- 1 - pnorm(threshold_obesidad_03, mean = datos03$pesoedad_pred, sd = sd_pred_03)

# Crear subsets del data frame con las variables de riesgo
r_datos_desnutricion_03 <- datos03[, c("ent", "prob_riesgo_desnutricion")]
r_datos_obesidad_03 <- datos03[, c("ent", "prob_riesgo_obesidad")]

# Guardar los dataframes con las predicciones y los riesgos
write_xlsx(r_datos_desnutricion_03, "datos_de_riesgos_desnutricion_0_3_pesoedad.xlsx")
write_xlsx(r_datos_obesidad_03, "datos_de_riesgos_obesidad_0_3_pesoedad.xlsx")

# Identificar los factores más importantes
importance_matrix_03 <- xgb.importance(feature_names = colnames(x_train_03), model = xgb_model_03)
top_10_factors_03 <- importance_matrix_03[1:10, ]

# Guardar los factores importantes en un archivo Excel
write_xlsx(top_10_factors_03, "factores_importantes_pesoedad_0_3.xlsx")

## MAPAS

mapaest <- sf::st_read("u_territorial_estados_mgn_inegi_2013.shp")
risk_desnutricion_03 <- read_xlsx("datos_de_riesgos_desnutricion_0_3_pesoedad.xlsx")
risk_obesidad_03 <- read_xlsx("datos_de_riesgos_obesidad_0_3_pesoedad.xlsx")

names(mapaest)[names(mapaest) == "cvegeoedo"] <- "ent"

mapaest$ent <- as.integer(mapaest$ent)
risk_desnutricion_03$ent <- as.integer(risk_desnutricion_03$ent)
risk_obesidad_03$ent <- as.integer(risk_obesidad_03$ent)

# Unir los datos
mapa_y_datos_desnutricion_03 <- dplyr::inner_join(mapaest, risk_desnutricion_03, by = "ent")
mapa_y_datos_obesidad_03 <- dplyr::inner_join(mapaest, risk_obesidad_03, by = "ent")

# Crear la gráfica de desnutrición
ggplot(mapa_y_datos_desnutricion_03) +
  geom_sf(aes(fill = prob_riesgo_desnutricion), color = "white", size = 0.2) +
  scale_fill_gradient(low = "#FFE7E7", high = "#C30010", name = "Probabilidad") + 
  labs(title = "Riesgo de desnutrición moderada a grave en niños y niñas de 0 a 3 años",
       subtitle = "Probabilidad de tener un peso para la edad por debajo de -2 desviaciones estándar de la mediana",
       caption = "Resultado del modelo XGBoost para el peso para la edad") +
  geom_text_repel(aes(label = paste0(round(prob_riesgo_desnutricion * 100, 1), "%"), 
                      geometry = st_geometry(mapa_y_datos_desnutricion_03)),
                  stat = "sf_coordinates", size = 4, fontface = "bold", 
                  nudge_x = 0.15, nudge_y = 0.15, max.overlaps = Inf,
                  segment.color = "grey50", segment.size = 0.5) +  
  theme_minimal() +
  theme(
    plot.title = element_text(size = 20, face = "bold"),       
    plot.subtitle = element_text(size = 18),                  
    axis.title = element_blank(),                             
    axis.text = element_blank(),                              
    axis.ticks = element_blank(),                             
    legend.title = element_text(size = 14),                  
    legend.text = element_text(size = 12),                    
    plot.caption = element_text(size = 14)                    
  )

# Crear la gráfica de obesidad
ggplot(mapa_y_datos_obesidad_03) +
  geom_sf(aes(fill = prob_riesgo_obesidad), color = "white", size = 0.2) +
  scale_fill_gradient(low = "#FFE7E7", high = "#C30010", name = "Probabilidad") +
  labs(title = "Riesgo de sobrepeso y obesidad en niños de 0 a 3 años",
       subtitle = "Probabilidad de tener un peso para la edad por encima de +1 desviaciones estándar de la mediana",
       caption = "Resultado del modelo XGBoost para el peso para la edad") +
  geom_text_repel(aes(label = paste0(round(prob_riesgo_obesidad * 100, 1), "%"), 
                      geometry = st_geometry(mapa_y_datos_obesidad_03)),
                  stat = "sf_coordinates", size = 4, fontface = "bold", 
                  nudge_x = 0.15, nudge_y = 0.15, max.overlaps = Inf,
                  segment.color = "grey50", segment.size = 0.5) +  
  theme_minimal() +
  theme(
    plot.title = element_text(size = 20, face = "bold"),       
    plot.subtitle = element_text(size = 18),                  
    axis.title = element_blank(),                             
    axis.text = element_blank(),                              
    axis.ticks = element_blank(),                             
    legend.title = element_text(size = 14),                   
    legend.text = element_text(size = 12),                    
    plot.caption = element_text(size = 14)                    
  )

### Procesar y mapear para el grupo de edad 4-6 años

# Selección de variables independientes y dependientes
x_46 <- datos46[, 9:146]
y_46 <- datos46$pesoedad

# Eliminar filas con valores NA en la variable dependiente o en las predictoras
data_complete_46 <- datos46[complete.cases(x_46, y_46), ]
x_complete_46 <- data_complete_46[, 9:143]
y_complete_46 <- data_complete_46$pesoedad

# Datos de entrenamiento y prueba
set.seed(123)  # Para reproducibilidad
trainIndex_46 <- createDataPartition(y_complete_46, p = .8, list = FALSE, times = 1)
x_train_46 <- x_complete_46[trainIndex_46, ]
x_test_46 <- x_complete_46[-trainIndex_46, ]
y_train_46 <- y_complete_46[trainIndex_46]
y_test_46 <- y_complete_46[-trainIndex_46]

# Convertir los datos a la matriz DMatrix de XGBoost
dtrain_46 <- xgb.DMatrix(data = as.matrix(x_train_46), label = y_train_46)
dtest_46 <- xgb.DMatrix(data = as.matrix(x_test_46), label = y_test_46)

# Entrenamiento del modelo XGBoost
params <- list(objective = "reg:squarederror", eval_metric = "rmse")
xgb_model_46 <- xgboost(data = dtrain_46, params = params, nrounds = 100, verbose = 0)

# Predicción con el modelo XGBoost
xgb_predictions_46 <- predict(xgb_model_46, newdata = dtest_46)

# Agregar predicciones al dataframe original
datos46$pesoedad_pred <- NA
datos46$pesoedad_pred[complete.cases(x_46, y_46)] <- predict(xgb_model_46, newdata = xgb.DMatrix(data = as.matrix(x_complete_46)))

# Calcular mediana y desviación estándar de las predicciones
median_pred_46 <- median(datos46$pesoedad_pred, na.rm = TRUE)
sd_pred_46 <- sd(datos46$pesoedad_pred, na.rm = TRUE)

# Calcular riesgo de desnutrición crónica (pesoedad < -2 desviaciones estándar)
threshold_desnutricion_46 <- median_pred_46 - 2 * sd_pred_46
datos46$prob_riesgo_desnutricion <- pnorm(threshold_desnutricion_46, mean = datos46$pesoedad_pred, sd = sd_pred_46)

# Calcular riesgo de sobrepeso y obesidad (pesoedad > +1 desviaciones estándar)
threshold_obesidad_46 <- median_pred_46 + 1 * sd_pred_46
datos46$prob_riesgo_obesidad <- 1 - pnorm(threshold_obesidad_46, mean = datos46$pesoedad_pred, sd = sd_pred_46)

# Crear subsets del data frame con las variables de riesgo
r_datos_desnutricion_46 <- datos46[, c("ent", "prob_riesgo_desnutricion")]
r_datos_obesidad_46 <- datos46[, c("ent", "prob_riesgo_obesidad")]

# Guardar los dataframes con las predicciones y los riesgos
write_xlsx(r_datos_desnutricion_46, "datos_de_riesgos_desnutricion_4_6_pesoedad.xlsx")
write_xlsx(r_datos_obesidad_46, "datos_de_riesgos_obesidad_4_6_pesoedad.xlsx")

# Identificar los factores más importantes
importance_matrix_46 <- xgb.importance(feature_names = colnames(x_train_46), model = xgb_model_46)
top_10_factors_46 <- importance_matrix_46[1:10, ]

# Guardar los factores importantes en un archivo Excel
write_xlsx(top_10_factors_46, "factores_importantes_pesoedad_4_6.xlsx")

## MAPAS

risk_desnutricion_46 <- read_xlsx("datos_de_riesgos_desnutricion_4_6_pesoedad.xlsx")
risk_obesidad_46 <- read_xlsx("datos_de_riesgos_obesidad_4_6_pesoedad.xlsx")

risk_desnutricion_46$ent <- as.integer(risk_desnutricion_46$ent)
risk_obesidad_46$ent <- as.integer(risk_obesidad_46$ent)

# Unir los datos
mapa_y_datos_desnutricion_46 <- dplyr::inner_join(mapaest, risk_desnutricion_46, by = "ent")
mapa_y_datos_obesidad_46 <- dplyr::inner_join(mapaest, risk_obesidad_46, by = "ent")

# Crear la gráfica de desnutrición
ggplot(mapa_y_datos_desnutricion_46) +
  geom_sf(aes(fill = prob_riesgo_desnutricion), color = "white", size = 0.2) +
  scale_fill_gradient(low = "#FFE7E7", high = "#C30010", name = "Probabilidad") +
  labs(title = "Riesgo de desnutrición moderada a grave en niños y niñas de 4 a 6 años",
       subtitle = "Probabilidad de tener un peso para la edad por debajo de -2 desviaciones estándar de la mediana",
       caption = "Resultado del modelo XGBoost para el peso para la edad") +
  geom_text_repel(aes(label = paste0(round(prob_riesgo_desnutricion * 100, 1), "%"), 
                      geometry = st_geometry(mapa_y_datos_desnutricion_46)),
                  stat = "sf_coordinates", size = 4, fontface = "bold", 
                  nudge_x = 0.15, nudge_y = 0.15, max.overlaps = Inf,
                  segment.color = "grey50", segment.size = 0.5) +  
  theme_minimal() +
  theme(
    plot.title = element_text(size = 20, face = "bold"),       
    plot.subtitle = element_text(size = 18),                 
    axis.title = element_blank(),                             
    axis.text = element_blank(),                              
    axis.ticks = element_blank(),                             
    legend.title = element_text(size = 14),                   
    legend.text = element_text(size = 12),                    
    plot.caption = element_text(size = 14)                    
  )

# Crear la gráfica de obesidad
ggplot(mapa_y_datos_obesidad_46) +
  geom_sf(aes(fill = prob_riesgo_obesidad), color = "white", size = 0.2) +
  scale_fill_gradient(low = "#FFE7E7", high = "#C30010", name = "Probabilidad") +
  labs(title = "Riesgo de sobrepeso y obesidad en niños de 4 a 6 años",
       subtitle = "Probabilidad de tener un peso para la edad por encima de +1 desviaciones estándar de la mediana",
       caption = "Resultado del modelo XGBoost para el peso para la edad") +
  geom_text_repel(aes(label = paste0(round(prob_riesgo_obesidad * 100, 1), "%"), 
                      geometry = st_geometry(mapa_y_datos_obesidad_46)),
                  stat = "sf_coordinates", size = 4, fontface = "bold", 
                  nudge_x = 0.15, nudge_y = 0.15, max.overlaps = Inf,
                  segment.color = "grey50", segment.size = 0.5) +  
  theme_minimal() +
  theme(
    plot.title = element_text(size = 20, face = "bold"),       
    plot.subtitle = element_text(size = 18),                  
    axis.title = element_blank(),                             
    axis.text = element_blank(),                              
    axis.ticks = element_blank(),                             
    legend.title = element_text(size = 14),                   
    legend.text = element_text(size = 12),                    
    plot.caption = element_text(size = 14)                    
  )

### Procesar y mapear para el grupo de edad 7-9 años

# Selección de variables independientes y dependientes
x_79 <- datos79[, 9:146]
y_79 <- datos79$pesoedad

# Eliminar filas con valores NA en la variable dependiente o en las predictoras
data_complete_79 <- datos79[complete.cases(x_79, y_79), ]
x_complete_79 <- data_complete_79[, 9:143]
y_complete_79 <- data_complete_79$pesoedad

# Datos de entrenamiento y prueba
set.seed(123)  # Para reproducibilidad
trainIndex_79 <- createDataPartition(y_complete_79, p = .8, list = FALSE, times = 1)
x_train_79 <- x_complete_79[trainIndex_79, ]
x_test_79 <- x_complete_79[-trainIndex_79, ]
y_train_79 <- y_complete_79[trainIndex_79]
y_test_79 <- y_complete_79[-trainIndex_79]

# Convertir los datos a la matriz DMatrix de XGBoost
dtrain_79 <- xgb.DMatrix(data = as.matrix(x_train_79), label = y_train_79)
dtest_79 <- xgb.DMatrix(data = as.matrix(x_test_79), label = y_test_79)

# Entrenamiento del modelo XGBoost
params <- list(objective = "reg:squarederror", eval_metric = "rmse")
xgb_model_79 <- xgboost(data = dtrain_79, params = params, nrounds = 100, verbose = 0)

# Predicción con el modelo XGBoost
xgb_predictions_79 <- predict(xgb_model_79, newdata = dtest_79)

# Agregar predicciones al dataframe original
datos79$pesoedad_pred <- NA
datos79$pesoedad_pred[complete.cases(x_79, y_79)] <- predict(xgb_model_79, newdata = xgb.DMatrix(data = as.matrix(x_complete_79)))

# Calcular mediana y desviación estándar de las predicciones
median_pred_79 <- median(datos79$pesoedad_pred, na.rm = TRUE)
sd_pred_79 <- sd(datos79$pesoedad_pred, na.rm = TRUE)

# Calcular riesgo de desnutrición crónica (pesoedad < -2 desviaciones estándar)
threshold_desnutricion_79 <- median_pred_79 - 2 * sd_pred_79
datos79$prob_riesgo_desnutricion <- pnorm(threshold_desnutricion_79, mean = datos79$pesoedad_pred, sd = sd_pred_79)

# Calcular riesgo de obesidad (pesoedad > +1 desviaciones estándar)
threshold_obesidad_79 <- median_pred_79 + 1 * sd_pred_79
datos79$prob_riesgo_obesidad <- 1 - pnorm(threshold_obesidad_79, mean = datos79$pesoedad_pred, sd = sd_pred_79)

# Crear subsets del data frame con las variables de riesgo
r_datos_desnutricion_79 <- datos79[, c("ent", "prob_riesgo_desnutricion")]
r_datos_obesidad_79 <- datos79[, c("ent", "prob_riesgo_obesidad")]

# Guardar los dataframes con las predicciones y los riesgos
write_xlsx(r_datos_desnutricion_79, "datos_de_riesgos_desnutricion_7_9_pesoedad.xlsx")
write_xlsx(r_datos_obesidad_79, "datos_de_riesgos_obesidad_7_9_pesoedad.xlsx")

# Identificar los factores más importantes
importance_matrix_79 <- xgb.importance(feature_names = colnames(x_train_79), model = xgb_model_79)
top_10_factors_79 <- importance_matrix_79[1:10, ]

# Guardar los factores importantes en un archivo Excel
write_xlsx(top_10_factors_79, "factores_importantes_pesoedad_7_9.xlsx")

## MAPAS

risk_desnutricion_79 <- read_xlsx("datos_de_riesgos_desnutricion_7_9_pesoedad.xlsx")
risk_obesidad_79 <- read_xlsx("datos_de_riesgos_obesidad_7_9_pesoedad.xlsx")

risk_desnutricion_79$ent <- as.integer(risk_desnutricion_79$ent)
risk_obesidad_79$ent <- as.integer(risk_obesidad_79$ent)

# Unir los datos
mapa_y_datos_desnutricion_79 <- dplyr::inner_join(mapaest, risk_desnutricion_79, by = "ent")
mapa_y_datos_obesidad_79 <- dplyr::inner_join(mapaest, risk_obesidad_79, by = "ent")

# Crear la gráfica de desnutrición
ggplot(mapa_y_datos_desnutricion_79) +
  geom_sf(aes(fill = prob_riesgo_desnutricion), color = "white", size = 0.2) +
  scale_fill_gradient(low = "#FFE7E7", high = "#C30010", name = "Probabilidad") +
  labs(title = "Riesgo de desnutrición moderada a grave en niños y niñas de 7 a 9 años",
       subtitle = "Probabilidad de tener un peso para la edad por debajo de -2 desviaciones estándar de la mediana",
       caption = "Resultado del modelo XGBoost para el peso para la edad") +
  geom_text_repel(aes(label = paste0(round(prob_riesgo_desnutricion * 100, 1), "%"), 
                      geometry = st_geometry(mapa_y_datos_desnutricion_79)),
                  stat = "sf_coordinates", size = 4, fontface = "bold", 
                  nudge_x = 0.15, nudge_y = 0.15, max.overlaps = Inf,
                  segment.color = "grey50", segment.size = 0.5) +  
  theme_minimal() +
  theme(
    plot.title = element_text(size = 20, face = "bold"),       
    plot.subtitle = element_text(size = 18),                  
    axis.title = element_blank(),                             
    axis.text = element_blank(),                              
    axis.ticks = element_blank(),                             
    legend.title = element_text(size = 14),                   
    legend.text = element_text(size = 12),                    
    plot.caption = element_text(size = 14)                    
  )

# Crear la gráfica de obesidad
ggplot(mapa_y_datos_obesidad_79) +
  geom_sf(aes(fill = prob_riesgo_obesidad), color = "white", size = 0.2) +
  scale_fill_gradient(low = "#FFE7E7", high = "#C30010", name = "Probabilidad") +
  labs(title = "Riesgo de obesidad en niños de 7 a 9 años",
       subtitle = "Probabilidad de tener un peso para la edad por encima de +1 desviaciones estándar de la mediana",
       caption = "Resultado del modelo XGBoost para el peso para la edad") +
  geom_text_repel(aes(label = paste0(round(prob_riesgo_obesidad * 100, 1), "%"), 
                      geometry = st_geometry(mapa_y_datos_obesidad_79)),
                  stat = "sf_coordinates", size = 4, fontface = "bold", 
                  nudge_x = 0.15, nudge_y = 0.15, max.overlaps = Inf,
                  segment.color = "grey50", segment.size = 0.5) +  
  theme_minimal() +
  theme(
    plot.title = element_text(size = 20, face = "bold"),       
    plot.subtitle = element_text(size = 18),                  
    axis.title = element_blank(),                             
    axis.text = element_blank(),                              
    axis.ticks = element_blank(),                             
    legend.title = element_text(size = 14),                   
    legend.text = element_text(size = 12),                    
    plot.caption = element_text(size = 14)                    
  )
