knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE,
  warning = FALSE,
  fig.width = 8,
  fig.height = 6,
  dpi = 120
)

paquetes <- c(
  "dplyr",
  "ggplot2",
  "mgcv",
  "sf",
  "gstat",
  "sp"
)

faltan <- paquetes[!paquetes %in% rownames(installed.packages())]

if(length(faltan) > 0){
  install.packages(faltan, repos = "https://cloud.r-project.org")
}

library(dplyr)
library(ggplot2)
library(mgcv)
library(sf)
library(gstat)
library(sp)

archivo_csv <- "pntd.0004813.s004 (3)(2).csv"

if(!file.exists(archivo_csv)){
  stop(
    paste0(
      "No se encontró el archivo '", archivo_csv, "'. ",
      "Coloque el CSV en la misma carpeta que este archivo .Rmd antes de usar Knit."
    )
  )
}

datos <- read.csv(
  archivo_csv,
  stringsAsFactors = FALSE,
  check.names = FALSE
)

variables_necesarias <- c(
  "province",
  "district",
  "cluster_number",
  "bite_status",
  "envenoming_status",
  "cluster_centroid_gps_x",
  "cluster_centroid_gps_y"
)

variables_faltantes <- setdiff(variables_necesarias, names(datos))

if(length(variables_faltantes) > 0){
  stop(
    paste(
      "Faltan variables necesarias:",
      paste(variables_faltantes, collapse = ", ")
    )
  )
}

datos <- datos %>%
  mutate(
    bite_status = as.numeric(bite_status),
    envenoming_status = as.numeric(envenoming_status),
    cluster_centroid_gps_x = as.numeric(cluster_centroid_gps_x),
    cluster_centroid_gps_y = as.numeric(cluster_centroid_gps_y),
    province = factor(province),
    district = factor(district),
    cluster_id = paste(
      as.character(province),
      as.character(district),
      cluster_number,
      sep = "_"
    )
  )

cat("Filas:", nrow(datos), "\n")
cat("Columnas:", ncol(datos), "\n")
cat("Mordeduras:", sum(datos$bite_status, na.rm = TRUE), "\n")
cat("Envenenamientos:", sum(datos$envenoming_status, na.rm = TRUE), "\n")

primer_no_na <- function(x){
  y <- x[!is.na(x)]
  if(length(y) == 0) return(NA_real_)
  y[1]
}

datos_cluster <- datos %>%
  group_by(
    province,
    district,
    cluster_number,
    cluster_id
  ) %>%
  summarise(
    n_cluster = n(),
    mordeduras = sum(bite_status, na.rm = TRUE),
    envenenamientos = sum(envenoming_status, na.rm = TRUE),
    longitud = primer_no_na(cluster_centroid_gps_x),
    latitud = primer_no_na(cluster_centroid_gps_y),
    .groups = "drop"
  ) %>%
  mutate(
    no_mordeduras = n_cluster - mordeduras,
    no_envenenamientos = n_cluster - envenenamientos,
    proporcion_mordedura = mordeduras / n_cluster,
    proporcion_envenenamiento = envenenamientos / n_cluster,
    incidencia_mordedura_100k = proporcion_mordedura * 100000,
    incidencia_envenenamiento_100k = proporcion_envenenamiento * 100000,
    province = factor(province),
    district = factor(district)
  )

cat("Clusters:", nrow(datos_cluster), "\n")
cat("Personas:", sum(datos_cluster$n_cluster), "\n")
cat("Mordeduras:", sum(datos_cluster$mordeduras), "\n")
cat("Envenenamientos:", sum(datos_cluster$envenenamientos), "\n")

datos_espaciales <- datos_cluster %>%
  filter(
    !is.na(longitud),
    !is.na(latitud),
    longitud >= -180,
    longitud <= 180,
    latitud >= -90,
    latitud <= 90
  )

datos_sf <- st_as_sf(
  datos_espaciales,
  coords = c("longitud", "latitud"),
  crs = 4326,
  remove = FALSE
)

datos_utm <- st_transform(datos_sf, 32644)
datos_sp <- as(datos_utm, "Spatial")

cat("Clusters totales:", nrow(datos_cluster), "\n")
cat("Clusters con coordenadas válidas:", nrow(datos_espaciales), "\n")
cat("Clusters excluidos por falta/coordenadas inválidas:",
    nrow(datos_cluster) - nrow(datos_espaciales), "\n")

modelo_glm_mordedura <- glm(
  cbind(mordeduras, no_mordeduras) ~ province,
  family = binomial(link = "logit"),
  data = datos_cluster
)

summary(modelo_glm_mordedura)

prueba_glm_mordedura <- anova(
  modelo_glm_mordedura,
  test = "Chisq"
)

prueba_glm_mordedura

modelo_glm_envenenamiento <- glm(
  cbind(envenenamientos, no_envenenamientos) ~ province,
  family = binomial(link = "logit"),
  data = datos_cluster
)

summary(modelo_glm_envenenamiento)

prueba_glm_envenenamiento <- anova(
  modelo_glm_envenenamiento,
  test = "Chisq"
)

prueba_glm_envenenamiento

modelo_gam_mordedura <- gam(
  cbind(mordeduras, no_mordeduras) ~
    s(longitud, latitud, k = 50),
  family = binomial(link = "logit"),
  method = "REML",
  data = datos_espaciales
)

summary(modelo_gam_mordedura)
gam.check(modelo_gam_mordedura)

plot(
  modelo_gam_mordedura,
  select = 1,
  scheme = 2,
  too.far = 0.1,
  main = "GAM espacial - Mordeduras"
)

modelo_gam_envenenamiento <- gam(
  cbind(envenenamientos, no_envenenamientos) ~
    s(longitud, latitud, k = 50),
  family = binomial(link = "logit"),
  method = "REML",
  data = datos_espaciales
)

summary(modelo_gam_envenenamiento)
gam.check(modelo_gam_envenenamiento)

plot(
  modelo_gam_envenenamiento,
  select = 1,
  scheme = 2,
  too.far = 0.1,
  main = "GAM espacial - Envenenamiento"
)

datos_provincia <- datos %>%
  group_by(province) %>%
  summarise(
    n = n(),
    mordeduras = sum(bite_status, na.rm = TRUE),
    envenenamientos = sum(envenoming_status, na.rm = TRUE),
    incidencia_mordedura_100k = mordeduras / n * 100000,
    incidencia_envenenamiento_100k = envenenamientos / n * 100000,
    .groups = "drop"
  )

ggplot(
  datos_provincia,
  aes(
    x = reorder(province, incidencia_mordedura_100k),
    y = incidencia_mordedura_100k
  )
) +
  geom_col() +
  coord_flip() +
  labs(
    title = "Incidencia cruda de mordeduras por provincia",
    x = "Provincia",
    y = "Mordeduras por 100,000"
  ) +
  theme_minimal()

ggplot(
  datos_provincia,
  aes(
    x = reorder(province, incidencia_envenenamiento_100k),
    y = incidencia_envenenamiento_100k
  )
) +
  geom_col() +
  coord_flip() +
  labs(
    title = "Incidencia cruda de envenenamiento por provincia",
    x = "Provincia",
    y = "Envenenamientos por 100,000"
  ) +
  theme_minimal()

p_glm_mord <- tail(na.omit(prueba_glm_mordedura$`Pr(>Chi)`), 1)
p_glm_env  <- tail(na.omit(prueba_glm_envenenamiento$`Pr(>Chi)`), 1)

cat("GLM mordeduras - p global de provincia:", p_glm_mord, "\n")
cat("GLM envenenamiento - p global de provincia:", p_glm_env, "\n")

cat("\nGAM mordeduras:\n")
print(summary(modelo_gam_mordedura)$s.table)

cat("\nGAM envenenamiento:\n")
print(summary(modelo_gam_envenenamiento)$s.table)

ggplot(datos_sf) +
  geom_sf(
    aes(size = incidencia_mordedura_100k),
    alpha = 0.5
  ) +
  labs(
    title = "Distribución espacial de la incidencia de mordedura",
    subtitle = "Incidencia cruda observada por cluster",
    size = "Mordeduras\npor 100,000",
    x = "Longitud",
    y = "Latitud"
  ) +
  theme_minimal()

ggplot(datos_sf) +
  geom_sf(
    aes(size = incidencia_envenenamiento_100k),
    alpha = 0.5
  ) +
  labs(
    title = "Distribución espacial de la incidencia de envenenamiento",
    subtitle = "Incidencia cruda observada por cluster",
    size = "Envenenamientos\npor 100,000",
    x = "Longitud",
    y = "Latitud"
  ) +
  theme_minimal()

variograma_mordedura <- variogram(
  proporcion_mordedura ~ 1,
  data = datos_sp
)

variograma_mordedura

plot(
  variograma_mordedura,
  main = "Variograma empírico - Mordeduras",
  xlab = "Distancia entre clusters (m)",
  ylab = "Semivarianza"
)

variograma_envenenamiento <- variogram(
  proporcion_envenenamiento ~ 1,
  data = datos_sp
)

variograma_envenenamiento

plot(
  variograma_envenenamiento,
  main = "Variograma empírico - Envenenamiento",
  xlab = "Distancia entre clusters (m)",
  ylab = "Semivarianza"
)

cat("Número de intervalos del variograma de mordedura:",
    nrow(variograma_mordedura), "\n")

cat("Número de intervalos del variograma de envenenamiento:",
    nrow(variograma_envenenamiento), "\n")

cat("\nPrimeras filas - mordedura:\n")
print(head(variograma_mordedura))

cat("\nPrimeras filas - envenenamiento:\n")
print(head(variograma_envenenamiento))

var_mord <- var(
  datos_espaciales$proporcion_mordedura,
  na.rm = TRUE
)

dist_max_mord <- max(
  variograma_mordedura$dist,
  na.rm = TRUE
)

modelo_inicial_mordedura <- vgm(
  psill = max(var_mord * 0.8, .Machine$double.eps),
  model = "Exp",
  range = max(dist_max_mord / 3, 1),
  nugget = max(var_mord * 0.2, 0)
)

modelo_mordedura <- fit.variogram(
  variograma_mordedura,
  modelo_inicial_mordedura
)

modelo_mordedura

plot(
  variograma_mordedura,
  modelo_mordedura,
  main = "Variograma y modelo exponencial - Mordeduras",
  xlab = "Distancia entre clusters (m)",
  ylab = "Semivarianza"
)

var_env <- var(
  datos_espaciales$proporcion_envenenamiento,
  na.rm = TRUE
)

dist_max_env <- max(
  variograma_envenenamiento$dist,
  na.rm = TRUE
)

modelo_inicial_envenenamiento <- vgm(
  psill = max(var_env * 0.8, .Machine$double.eps),
  model = "Exp",
  range = max(dist_max_env / 3, 1),
  nugget = max(var_env * 0.2, 0)
)

modelo_envenenamiento <- fit.variogram(
  variograma_envenenamiento,
  modelo_inicial_envenenamiento
)

modelo_envenenamiento

plot(
  variograma_envenenamiento,
  modelo_envenenamiento,
  main = "Variograma y modelo exponencial - Envenenamiento",
  xlab = "Distancia entre clusters (m)",
  ylab = "Semivarianza"
)

paso_malla <- 10000

malla_sf <- st_make_grid(
  datos_utm,
  cellsize = paso_malla,
  what = "centers"
)

malla_sf <- st_sf(geometry = malla_sf)

envolvente <- st_convex_hull(st_union(datos_utm))

malla_sf <- malla_sf[
  lengths(st_intersects(malla_sf, envolvente)) > 0,
]

malla_sp <- as(malla_sf, "Spatial")

cat("Puntos de predicción:", nrow(malla_sf), "\n")

kriging_mordedura <- krige(
  proporcion_mordedura ~ 1,
  locations = datos_sp,
  newdata = malla_sp,
  model = modelo_mordedura
)

kriging_mordedura_sf <- st_as_sf(kriging_mordedura)

kriging_mordedura_sf$prediccion <- pmin(
  pmax(kriging_mordedura_sf$var1.pred, 0),
  1
)

kriging_mordedura_sf$incidencia_100k <-
  kriging_mordedura_sf$prediccion * 100000

summary(kriging_mordedura_sf$incidencia_100k)

ggplot() +
  geom_sf(
    data = kriging_mordedura_sf,
    aes(color = incidencia_100k),
    size = 3
  ) +
  geom_sf(
    data = datos_utm,
    size = 0.3,
    alpha = 0.25
  ) +
  labs(
    title = "Predicción espacial parcial de mordeduras",
    subtitle = "Kriging ordinario sobre la proporción observada por cluster",
    color = "Predicción\npor 100,000"
  ) +
  theme_minimal()

kriging_envenenamiento <- krige(
  proporcion_envenenamiento ~ 1,
  locations = datos_sp,
  newdata = malla_sp,
  model = modelo_envenenamiento
)

kriging_envenenamiento_sf <- st_as_sf(kriging_envenenamiento)

kriging_envenenamiento_sf$prediccion <- pmin(
  pmax(kriging_envenenamiento_sf$var1.pred, 0),
  1
)

kriging_envenenamiento_sf$incidencia_100k <-
  kriging_envenenamiento_sf$prediccion * 100000

summary(kriging_envenenamiento_sf$incidencia_100k)

ggplot() +
  geom_sf(
    data = kriging_envenenamiento_sf,
    aes(color = incidencia_100k),
    size = 3
  ) +
  geom_sf(
    data = datos_utm,
    size = 0.3,
    alpha = 0.25
  ) +
  labs(
    title = "Predicción espacial parcial de envenenamiento",
    subtitle = "Kriging ordinario sobre la proporción observada por cluster",
    color = "Predicción\npor 100,000"
  ) +
  theme_minimal()

cat("MODELO DE MORDERURAS\n")
print(modelo_mordedura)

cat("\nMODELO DE ENVENENAMIENTO\n")
print(modelo_envenenamiento)

extraer_parametros <- function(modelo){
  nugget <- modelo$psill[modelo$model == "Nug"]
  psill_espacial <- sum(modelo$psill[modelo$model != "Nug"])
  sill_total <- sum(modelo$psill)
  range_modelo <- max(modelo$range[modelo$model != "Nug"])

  data.frame(
    nugget = ifelse(length(nugget) == 0, 0, nugget),
    psill_espacial = psill_espacial,
    sill_total = sill_total,
    range_modelo_m = range_modelo,
    rango_practico_aprox_m = 3 * range_modelo
  )
}

parametros_mordedura <- extraer_parametros(modelo_mordedura)
parametros_envenenamiento <- extraer_parametros(modelo_envenenamiento)

cat("\nParámetros de mordeduras:\n")
print(parametros_mordedura)

cat("\nParámetros de envenenamiento:\n")
print(parametros_envenenamiento)

cat(
"
ANÁLISIS 1 — GLM/GAM
Pregunta:
¿Existen diferencias geográficas y patrones no lineales
en mordedura y envenenamiento?

Métodos:
GLM binomial + GAM binomial.

Limitación:
No están disponibles las covariables originales.


ANÁLISIS 2 — VARIOGRAMAS
Pregunta:
¿Existe estructura espacial dependiente de la distancia?

Método:
Variograma empírico de la proporción observada por cluster.

Limitación:
No se reproducen exactamente los residuos del GLM
multivariable original.


ANÁLISIS 3 — MODELO GEOESTADÍSTICO
Pregunta:
¿Puede modelarse y predecirse parcialmente la distribución
espacial utilizando la estructura de dependencia observada?

Métodos:
Modelo exponencial de variograma + kriging ordinario.

Limitación:
Es una aproximación espacial de las proporciones observadas,
no el modelo geoestadístico binomial multivariable original.
"
)

write.csv(
  datos_cluster,
  "datos_cluster_analisis_1_2_3.csv",
  row.names = FALSE
)

write.csv(
  variograma_mordedura,
  "variograma_mordedura.csv",
  row.names = FALSE
)

write.csv(
  variograma_envenenamiento,
  "variograma_envenenamiento.csv",
  row.names = FALSE
)

write.csv(
  parametros_mordedura,
  "parametros_geo_mordedura.csv",
  row.names = FALSE
)

write.csv(
  parametros_envenenamiento,
  "parametros_geo_envenenamiento.csv",
  row.names = FALSE
)

cat("Archivos de resultados guardados correctamente.\n")

sessionInfo()