Diseño experimental. Flujo general del estudio: muestreo por conglomerados → recolección de datos → GLM/GAM → evaluación de dependencia espacial → modelo geoestadístico y predicción.
GLM/GAM · estructura espacial · variogramas · modelamiento geoestadístico parcial
Incluye Análisis 1, 2 y 3 · Análisis 4 excluido
Cómo está organizado: cada análisis mantiene el mismo orden lógico: pregunta → hipótesis → selección del método → código R → gráfica → resultados → discusión.
Este documento reproduce, hasta donde lo permiten los datos disponibles en el archivo CSV suplementario, tres etapas del análisis estadístico de Ediriweera et al. (2016):
Limitación fundamental: el CSV utilizado contiene los desenlaces de mordedura y envenenamiento, identificadores geográficos y coordenadas de los clusters, pero no contiene todas las covariables ambientales y sociodemográficas del modelo original (elevación, densidad poblacional, zona climática, porcentaje de trabajadores agrícolas, NDVI, uso/cobertura del suelo, etc.). Por tanto, los modelos multivariables originales no pueden reproducirse exactamente con este archivo. Las partes que sí se ejecutan aquí utilizan exclusivamente variables presentes en el CSV.
El siguiente esquema resume visualmente el diseño del estudio y el orden de los análisis: desde la encuesta nacional y la recolección de datos, pasando por el análisis exploratorio con GLM/GAM y el análisis espacial con residuos y variograma, hasta el modelamiento geoestadístico y la predicción de riesgo.
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)El objeto principal se llamará siempre
datos. Ningún análisis vuelve a importar o
renombrar la base.
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")Los análisis espaciales requieren trabajar a nivel de cluster. Este
objeto se llamará siempre
datos_cluster.
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")¿La probabilidad de mordedura y de envenenamiento varía geográficamente entre las provincias y existe además un patrón espacial no lineal asociado con la localización de los clusters?
Esta pregunta constituye una replicación parcial del enfoque GLM/GAM, porque las covariables ambientales originales no están disponibles en el CSV.
Para el GLM:
[ H_0:_1=_2==_k=0 ]
frente a:
[ H_1:_j ]
Para el término suavizado del GAM:
[ H_0:f(longitud,latitud)=0 ]
frente a:
[ H_1:f(longitud,latitud) ]
Se utiliza un GLM binomial con enlace logit porque cada cluster contiene un número de eventos y no eventos. Se utiliza un GAM binomial para explorar variación geográfica no lineal usando las coordenadas disponibles.
Estos modelos no son una reproducción exacta de los GLM/GAM multivariables originales del artículo.
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_mordeduramodelo_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_envenenamientomodelo_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)El GLM permite determinar si la probabilidad observada de mordedura o envenenamiento difiere entre provincias. El GAM permite explorar si las coordenadas geográficas presentan una asociación suave y no lineal con cada desenlace.
La interpretación debe hacerse separadamente para mordedura y envenenamiento. Si el valor p del GLM es menor que 0.05, existe evidencia de diferencias provinciales. En el GAM, un término suavizado con valor p menor que 0.05 aporta evidencia de un patrón geográfico no lineal.
Limitación: estos resultados no pueden atribuirse a elevación, clima, densidad poblacional o agricultura porque dichas covariables no están en el CSV.
¿Existe una estructura espacial dependiente de la distancia en la incidencia observada de mordeduras y envenenamientos entre clusters?
[ H_0:Cov[Z(s_i),Z(s_j)]=0 ]
frente a:
[ H_1:Cov[Z(s_i),Z(s_j)] ]
para determinadas distancias entre localizaciones.
Se utiliza un variograma empírico porque las observaciones están georreferenciadas y se desea estudiar cómo cambia su similitud con la distancia.
El variograma es fundamentalmente una herramienta de diagnóstico y
modelamiento espacial; por sí solo no se interpreta mediante una regla
automática de p < 0.05.
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_mordeduraplot(
variograma_mordedura,
main = "Variograma empírico - Mordeduras",
xlab = "Distancia entre clusters (m)",
ylab = "Semivarianza"
)
variograma_envenenamiento <- variogram(
proporcion_envenenamiento ~ 1,
data = datos_sp
)
variograma_envenenamientoplot(
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))Si la semivarianza es relativamente baja a distancias cortas y aumenta conforme crece la distancia, existe evidencia descriptiva de que clusters cercanos presentan valores más similares. Esto apoya la existencia de estructura espacial y justifica el modelamiento geoestadístico posterior.
Si el variograma es aproximadamente plano desde las distancias más pequeñas, la evidencia de dependencia espacial será débil.
Limitación: aquí se analiza el variograma de la incidencia observada. No son los residuos estandarizados del GLM multivariable original, porque las covariables necesarias para reconstruir dicho modelo no se encuentran en el CSV.
¿Puede modelarse la distribución espacial de la incidencia observada de mordeduras y envenenamientos y utilizarse dicha estructura para realizar predicciones espaciales entre los clusters muestreados?
[ Y(s)=+S(s)+(s) ]
con:
[ H_0:_S^2=0 ]
frente a:
[ H_1:_S^2>0 ]
Nota sobre las gráficas integradas: el código
mostrado es el código R del R Markdown original. Las gráficas se
calcularon con el mismo CSV para que el HTML pueda abrirse ya ilustrado.
Como este entorno no dispone de R, las superficies suaves y predictivas
se renderizaron con equivalentes numéricos para visualización. Al
ejecutar el .Rmd en RStudio, las figuras exactas serán las
producidas por mgcv, sf y
gstat.
Con las variables disponibles se utiliza:
Esto representa una replicación geoestadística parcial. No es el modelo binomial-logístico geoestadístico multivariable original del artículo.
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_mordeduraplot(
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_envenenamientoplot(
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)Mordeduras y envenenamientos se modelan por separado porque son dos desenlaces diferentes. Una mordedura no implica necesariamente envenenamiento, y cada fenómeno puede presentar distinta intensidad y alcance de autocorrelación espacial.
Los parámetros del variograma permiten describir el
nugget, el componente de varianza espacial, el
sill y el range. En un modelo
exponencial, el rango práctico puede aproximarse como
3 × range.
El kriging utiliza esa estructura para interpolar valores entre los sitios observados.
Limitación importante: el kriging ordinario aplicado a proporciones no incorpora directamente la distribución binomial de los casos ni el tamaño diferente de los clusters como lo haría un modelo geoestadístico binomial. Por ello, estos mapas son exploratorios y constituyen una réplica parcial, no los mapas predictivos originales del artículo.
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()“En conclusión, los autores utilizaron una estrategia estadística coherente y progresiva. Primero aplicaron GLM y GAM para identificar asociaciones y relaciones no lineales. Posteriormente evaluaron la existencia de patrones espaciales mediante mapas de residuos y variogramas. Finalmente incorporaron esa dependencia espacial mediante un modelo geoestadístico, con el cual generaron mapas de riesgo de mordedura y envenenamiento. Por lo tanto, no consideramos que las pruebas utilizadas fueran incorrectas; más bien, proponemos como mejoras una validación espacial más explícita y pruebas complementarias de autocorrelación. Nuestra replicación fue parcial debido a que la base suplementaria no contiene todas las covariables del estudio original.”