Auditoría crítica del análisis epidemiológico y espacial del accidente ofídico en Eswatini
Reanálisis de Padidar et al. (2023)
Reinaldo Pérez, Eduardo Leiva, John Cleghorn y Randhy Rodríguez
04 de septiembre de 2026
1. Presentación del artículo
El estudio evaluado es el de Padidar et al. (2023), titulado Snakebite epidemiology, outcomes and multi-cluster risk modelling in Eswatini. Los autores describieron 932 mordeduras de serpiente registradas entre 2019 y 2021 y desarrollaron modelos epidemiológicos y mapas de riesgo espacial.
2. Objetivo de la auditoría
El objetivo es evaluar la organización y calidad de las bases públicas, recalcular los principales resultados, identificar posibles errores de captura, análisis o interpretación y proponer mejoras reproducibles.
La revisión diferencia tres tipos de hallazgos:Tres ejes de la auditoría
1
Errores comprobables
Valores o resultados que no coinciden entre la base pública y el artículo.
Ejemplo: porcentajes, frecuencias o rangos diferentes.
2
Información insuficiente
Procedimientos, denominadores o transformaciones que no se describen con suficiente detalle.
Ejemplo: subconjunto espacial no identificado.
3
Limitaciones metodológicas
Decisiones analíticas que reducen la validez o el alcance de las conclusiones.
Ejemplo: ausencia de población expuesta.
3. Importación y descripción de las bases
3.1 Localización y lectura de archivos
buscar_archivo <- function(nombre) {
candidatos <- c(
nombre,
file.path("data", nombre),
file.path("upload", nombre)
)
encontrados <- candidatos[file.exists(candidatos)]
if (length(encontrados) == 0) {
stop(paste("No se encontró el archivo:", nombre))
}
encontrados[1]
}
ruta_clinica <- buscar_archivo("Eswatini Snakebite Harvard.xlsx")
ruta_geografica <- buscar_archivo(
"Eswatini Snakebite Tin Elevation Harvard.xlsx"
)
df_clinical_raw <- read_excel(ruta_clinica, sheet = 1)
df_elevation_raw <- read_excel(ruta_geografica, sheet = 1)
estructura_bases <- tibble(
base = c("Clínica-epidemiológica", "Geográfica-altitudinal"),
filas = c(nrow(df_clinical_raw), nrow(df_elevation_raw)),
columnas = c(ncol(df_clinical_raw), ncol(df_elevation_raw)),
hoja = c(excel_sheets(ruta_clinica)[1], excel_sheets(ruta_geografica)[1])
)
kable(estructura_bases, caption = "Estructura de las bases públicas")| base | filas | columnas | hoja |
|---|---|---|---|
| Clínica-epidemiológica | 932 | 29 | Sheet1 |
| Geográfica-altitudinal | 932 | 4 | Snakebite |
3.2 Estandarización de nombres
La base original se conserva sin modificaciones. La limpieza se realiza sobre copias de trabajo para mantener la trazabilidad.
limpiar_nombres <- function(x) {
x <- trimws(tolower(x))
x <- iconv(x, from = "", to = "ASCII//TRANSLIT")
x <- gsub("[^a-z0-9]+", "_", x)
x <- gsub("^_+|_+$", "", x)
x
}
df_clinical <- df_clinical_raw
df_elevation <- df_elevation_raw
names(df_clinical) <- limpiar_nombres(names(df_clinical))
names(df_elevation) <- limpiar_nombres(names(df_elevation))
columnas_clinicas_esperadas <- c(
"age", "gender", "occupation_type", "date_snakebite",
"time_snakebite", "circumstances_snakebite", "position_snakebite",
"firstaid_any", "clinical_presentation", "av_administered",
"pre_med_adrenaline", "av_vials", "final_outcome"
)
columnas_geo_esperadas <- c("patient", "long", "lat", "elevation")
if (length(setdiff(columnas_clinicas_esperadas, names(df_clinical))) > 0) {
stop("La base clínica no contiene todas las columnas esperadas.")
}
if (length(setdiff(columnas_geo_esperadas, names(df_elevation))) > 0) {
stop("La base geográfica no contiene las cuatro columnas esperadas.")
}3.3 Limpieza de categorías
normalizar_si_no <- function(x) {
x_chr <- trimws(tolower(as.character(x)))
case_when(
is.na(x) | x_chr == "" ~ NA_character_,
x_chr == "yes" ~ "Yes",
x_chr == "no" ~ "No",
TRUE ~ trimws(as.character(x))
)
}
variables_si_no <- intersect(
c(
"firstaid_any", "firstaid_tourniquet", "firstaid_bandage",
"firstaid_herbal_applied", "firstaid_herbal_ingested",
"firstaid_incision", "av_administered", "pre_med_adrenaline",
"surgery", "debridement", "fasciotomy", "skin_graft",
"amputation", "surgery_other"
),
names(df_clinical)
)
df_clinical <- df_clinical %>%
mutate(
record_order = row_number(),
age_raw = trimws(as.character(age)),
age_clean = case_when(
age_raw == ".11-19" ~ "10-19",
TRUE ~ age_raw
),
gender = trimws(as.character(gender)),
position_snakebite = trimws(tolower(as.character(position_snakebite))),
circumstances_snakebite = trimws(as.character(circumstances_snakebite)),
clinical_presentation = trimws(as.character(clinical_presentation)),
final_outcome = trimws(as.character(final_outcome)),
date_snakebite = as.Date(date_snakebite)
) %>%
mutate(across(all_of(variables_si_no), normalizar_si_no))
df_elevation <- df_elevation %>%
mutate(
patient = as.integer(patient),
long = as.numeric(long),
lat = as.numeric(lat),
elevation = as.numeric(elevation)
)
secuencia_geo_correcta <- isTRUE(all(
df_elevation$patient == seq_len(nrow(df_elevation))
))3.4 Vinculación provisional
La base clínica no contiene un identificador de paciente. Por ello, no es posible demostrar que cada fila clínica corresponde al mismo paciente de la base geográfica. Para fines exploratorios se conserva el orden original, pero esta unión no debe presentarse como una vinculación confirmada.
df_combined <- df_clinical %>%
left_join(
df_elevation %>%
transmute(
record_order = patient,
patient_geo = patient,
long,
lat,
elevation
),
by = "record_order"
)
control_union <- tibble(
control = c(
"La base geográfica contiene pacientes consecutivos 1-932",
"La base clínica contiene un identificador explícito",
"Número de filas coincidente"
),
resultado = c(
secuencia_geo_correcta,
"patient" %in% names(df_clinical_raw),
nrow(df_clinical_raw) == nrow(df_elevation_raw)
)
)
kable(control_union, caption = "Controles de vinculación de las bases")| control | resultado |
|---|---|
| La base geográfica contiene pacientes consecutivos 1-932 | TRUE |
| La base clínica contiene un identificador explícito | FALSE |
| Número de filas coincidente | TRUE |
4. Evaluación de la calidad de los datos
4.1 Datos faltantes
es_faltante <- function(x) {
if (is.character(x)) {
x_limpio <- trimws(tolower(x))
is.na(x) | x_limpio %in% c("", "unknown", "not recorded")
} else {
is.na(x)
}
}
tabla_faltantes <- tibble(
variable = names(df_clinical_raw),
faltantes = vapply(
df_clinical_raw,
function(x) sum(es_faltante(x)),
numeric(1)
)
) %>%
mutate(
porcentaje = 100 * faltantes / nrow(df_clinical_raw)
) %>%
arrange(desc(porcentaje))
variables_clave <- c(
"Age", "Gender", "Occupation_type", "Date_snakebite",
"Time_snakebite", "Circumstances_snakebite", "Position_snakebite",
"Clinical_presentation", "AV_administered", "Final_outcome"
)
tabla_faltantes_clave <- tabla_faltantes %>%
filter(variable %in% variables_clave) %>%
mutate(variable = factor(variable, levels = rev(variables_clave)))
kable(
tabla_faltantes,
digits = 1,
caption = "Cantidad y porcentaje de datos faltantes por variable"
)| variable | faltantes | porcentaje |
|---|---|---|
| Other | 896 | 96.1 |
| AV_reaction3 | 324 | 34.8 |
| AV_reaction_other | 321 | 34.4 |
| AV_reaction2 | 303 | 32.5 |
| Time_snakebite | 282 | 30.3 |
| Circumstances_snakebite | 194 | 20.8 |
| Occupation_type | 168 | 18.0 |
| Age | 93 | 10.0 |
| FirstAid_Herbal_ingested | 74 | 7.9 |
| Final_outcome | 43 | 4.6 |
| FirstAid_incision | 41 | 4.4 |
| Clinical_presentation | 41 | 4.4 |
| AV_reaction1 | 40 | 4.3 |
| FirstAid_Any | 31 | 3.3 |
| FirstAid_Tourniquet | 31 | 3.3 |
| FirstAid_Bandage | 31 | 3.3 |
| FirstAid_Herbal_applied | 31 | 3.3 |
| Pre-med_adrenaline | 29 | 3.1 |
| AV_vials | 29 | 3.1 |
| Position_snakebite | 25 | 2.7 |
| Surgery | 24 | 2.6 |
| Debridement | 24 | 2.6 |
| Fasciotomy | 24 | 2.6 |
| Skin_graft | 24 | 2.6 |
| Surgery_Other | 24 | 2.6 |
| AV_administered | 23 | 2.5 |
| Amputation | 18 | 1.9 |
| Gender | 0 | 0.0 |
| Date_snakebite | 0 | 0.0 |
Los espacios en blanco de variables de respuesta múltiple, como tipos de primeros auxilios, reacciones y procedimientos quirúrgicos, pueden significar ausencia del evento y no necesariamente pérdida de información. Por esa razón, la gráfica se limita a variables principales y no interpreta automáticamente todos los espacios en blanco como datos faltantes.
4.2 Duplicados, coordenadas y consistencia lógica
geo_completo <- complete.cases(
df_elevation[, c("long", "lat", "elevation")]
)
columnas_reaccion <- intersect(
c("av_reaction1", "av_reaction2", "av_reaction3", "av_reaction_other"),
names(df_clinical)
)
reaccion_registrada <- Reduce(
`|`,
lapply(
df_clinical[columnas_reaccion],
function(x) !es_faltante(as.character(x))
)
)
controles_calidad <- tibble(
control = c(
"Filas completamente duplicadas en la base clínica",
"Filas completamente duplicadas en la base geográfica",
"Registros geográficos completos",
"Registros sin coordenadas ni elevación",
"Pares de coordenadas diferentes",
"Antiveneno administrado con número de viales faltante",
"Reacción registrada en personas sin antiveneno",
"Resultados clínicos finales faltantes"
),
resultado = c(
sum(duplicated(df_clinical_raw)),
sum(duplicated(df_elevation_raw)),
sum(geo_completo),
sum(!geo_completo),
nrow(distinct(filter(df_elevation, geo_completo), long, lat)),
sum(
df_clinical$av_administered == "Yes" & is.na(df_clinical$av_vials),
na.rm = TRUE
),
sum(reaccion_registrada & df_clinical$av_administered != "Yes", na.rm = TRUE),
sum(es_faltante(df_clinical$final_outcome))
)
)
kable(controles_calidad, caption = "Controles básicos de calidad y coherencia")| control | resultado |
|---|---|
| Filas completamente duplicadas en la base clínica | 0 |
| Filas completamente duplicadas en la base geográfica | 0 |
| Registros geográficos completos | 862 |
| Registros sin coordenadas ni elevación | 70 |
| Pares de coordenadas diferentes | 59 |
| Antiveneno administrado con número de viales faltante | 8 |
| Reacción registrada en personas sin antiveneno | 709 |
| Resultados clínicos finales faltantes | 43 |
Interpretación
La estructura general es rectangular y no presenta duplicados completos. Sin embargo, la calidad se reduce por la ausencia de un identificador compartido, los valores faltantes y las categorías inconsistentes. La repetición de coordenadas demuestra que los registros fueron asignados a localidades o áreas administrativas, por lo que no deben considerarse puntos espaciales independientes.
5. Reproducción de los resultados descriptivos
5.1 Perfil demográfico
n_total <- nrow(df_clinical)
tabla_genero <- df_clinical %>%
mutate(gender_plot = if_else(is.na(gender), "Sin dato", gender)) %>%
count(gender_plot, name = "casos") %>%
mutate(porcentaje = 100 * casos / n_total)
orden_edades <- c(
"0-9", "10-19", "20-29", "30-39", "40-49", "50-59",
"60-69", "70-79", "80-89", "90-99", "100-109", "Sin dato"
)
tabla_edad <- df_clinical %>%
mutate(
age_plot = if_else(is.na(age_clean), "Sin dato", age_clean),
age_plot = factor(age_plot, levels = orden_edades)
) %>%
count(age_plot, name = "casos", .drop = FALSE) %>%
mutate(porcentaje = 100 * casos / n_total)
n_menores_30 <- sum(
df_clinical$age_clean %in% c("0-9", "10-19", "20-29"),
na.rm = TRUE
)
n_edad_conocida <- sum(!is.na(df_clinical$age_clean))
resumen_edad <- tibble(
indicador = c(
"Menores de 30 años sobre todos los registros",
"Menores de 30 años entre las edades conocidas"
),
numerador = c(n_menores_30, n_menores_30),
denominador = c(n_total, n_edad_conocida),
porcentaje = 100 * numerador / denominador
)
kable(tabla_genero, digits = 1, caption = "Distribución por género")| gender_plot | casos | porcentaje |
|---|---|---|
| Female | 421 | 45.2 |
| Male | 511 | 54.8 |
kable(resumen_edad, digits = 1, caption = "Cálculo de pacientes menores de 30 años")| indicador | numerador | denominador | porcentaje |
|---|---|---|---|
| Menores de 30 años sobre todos los registros | 542 | 932 | 58.2 |
| Menores de 30 años entre las edades conocidas | 542 | 839 | 64.6 |
5.2 Presentación clínica
traducir_presentacion <- function(x) {
case_when(
x == "None: Asymptomatic" ~ "Asintomático",
x == "Mild swelling" ~ "Inflamación leve",
x == "Painful progressive swelling (Cytotoxic)" ~
"Edema progresivo doloroso (citotóxico)",
x == "Progressive weakness (Neurotoxic)" ~
"Debilidad progresiva (neurotóxico)",
x == "Opthalmic (Venom in the eye)" ~ "Exposición oftálmica",
x == "Bleeding (Hemotoxic)" ~ "Sangrado (hemotóxico)",
TRUE ~ x
)
}
tabla_presentacion <- df_clinical %>%
filter(!es_faltante(clinical_presentation)) %>%
mutate(presentacion = traducir_presentacion(clinical_presentation)) %>%
count(presentacion, name = "casos") %>%
mutate(porcentaje = 100 * casos / sum(casos)) %>%
arrange(desc(casos))
kable(
tabla_presentacion,
digits = 1,
caption = "Presentaciones clínicas entre los registros con información disponible"
)| presentacion | casos | porcentaje |
|---|---|---|
| Asintomático | 329 | 36.9 |
| Inflamación leve | 277 | 31.1 |
| Edema progresivo doloroso (citotóxico) | 218 | 24.5 |
| Debilidad progresiva (neurotóxico) | 44 | 4.9 |
| Exposición oftálmica | 21 | 2.4 |
| Sangrado (hemotóxico) | 2 | 0.2 |
De los registros con presentación clínica conocida, los cuadros asintomáticos y la inflamación leve fueron los más frecuentes. La gráfica describe categorías registradas y no permite identificar por sí sola la especie responsable ni confirmar la gravedad clínica individual.
5.3 Manejo con antiveneno y desenlaces
tabla_antiveneno <- df_clinical %>%
mutate(
av_plot = case_when(
av_administered == "Yes" ~ "Sí",
av_administered == "No" ~ "No",
TRUE ~ "Sin dato"
)
) %>%
count(av_plot, name = "casos") %>%
mutate(porcentaje = 100 * casos / n_total)
tabla_desenlace <- df_clinical %>%
mutate(
outcome_plot = case_when(
final_outcome == "Full recovery, no permanent damage" ~ "Recuperación completa",
final_outcome == "Death" ~ "Fallecimiento",
final_outcome == "Permanent disability" ~ "Discapacidad permanente",
TRUE ~ "Sin dato"
)
) %>%
count(outcome_plot, name = "casos") %>%
mutate(porcentaje = 100 * casos / n_total)
viales_tratados <- df_clinical %>%
filter(av_administered == "Yes", !is.na(av_vials))
resumen_viales <- tibble(
pacientes_tratados = sum(df_clinical$av_administered == "Yes", na.rm = TRUE),
pacientes_con_viales_registrados = nrow(viales_tratados),
media_viales = mean(viales_tratados$av_vials),
mediana_viales = median(viales_tratados$av_vials),
minimo = min(viales_tratados$av_vials),
maximo = max(viales_tratados$av_vials)
)
kable(tabla_antiveneno, digits = 1, caption = "Administración de antiveneno")| av_plot | casos | porcentaje |
|---|---|---|
| No | 709 | 76.1 |
| Sin dato | 23 | 2.5 |
| Sí | 200 | 21.5 |
kable(tabla_desenlace, digits = 1, caption = "Desenlace clínico final")| outcome_plot | casos | porcentaje |
|---|---|---|
| Discapacidad permanente | 8 | 0.9 |
| Fallecimiento | 10 | 1.1 |
| Recuperación completa | 871 | 93.5 |
| Sin dato | 43 | 4.6 |
kable(resumen_viales, digits = 2, caption = "Número de viales entre pacientes con información disponible")| pacientes_tratados | pacientes_con_viales_registrados | media_viales | mediana_viales | minimo | maximo |
|---|---|---|---|---|---|
| 200 | 192 | 5.19 | 5 | 1 | 25 |
5.4 Manejo según la presentación clínica
tabla_av_presentacion <- df_clinical %>%
filter(
!es_faltante(clinical_presentation),
av_administered %in% c("Yes", "No")
) %>%
mutate(
presentacion = traducir_presentacion(clinical_presentation),
antiveneno = recode(av_administered, "Yes" = "Sí", "No" = "No")
) %>%
count(presentacion, antiveneno, name = "casos") %>%
group_by(presentacion) %>%
mutate(
total = sum(casos),
porcentaje = 100 * casos / total
) %>%
ungroup()
tabla_cirugia_presentacion <- df_clinical %>%
filter(
!es_faltante(clinical_presentation),
surgery %in% c("Yes", "No")
) %>%
mutate(
presentacion = traducir_presentacion(clinical_presentation),
cirugia = recode(surgery, "Yes" = "Sí", "No" = "No")
) %>%
count(presentacion, cirugia, name = "casos") %>%
group_by(presentacion) %>%
mutate(
total = sum(casos),
porcentaje = 100 * casos / total
) %>%
ungroup()
La administración de antiveneno se concentró principalmente en los registros con debilidad progresiva y edema doloroso progresivo. Esta distribución muestra una asociación clínica, pero no demuestra por sí sola que todas las decisiones terapéuticas hayan sido apropiadas.
Las intervenciones quirúrgicas se concentraron en la presentación citotóxica. Los registros faltantes no fueron reclasificados como “No”, evitando reducir artificialmente la proporción de cirugía.
5.5 Primeros auxilios registrados
tabla_primeros_auxilios <- tibble(
tipo = c(
"Sin primeros auxilios",
"Torniquete",
"Vendaje",
"Hierbas aplicadas",
"Hierbas ingeridas",
"Incisión",
"Otro procedimiento"
),
casos = c(
sum(df_clinical$firstaid_any == "No", na.rm = TRUE),
sum(df_clinical$firstaid_tourniquet == "Yes", na.rm = TRUE),
sum(df_clinical$firstaid_bandage == "Yes", na.rm = TRUE),
sum(df_clinical$firstaid_herbal_applied == "Yes", na.rm = TRUE),
sum(df_clinical$firstaid_herbal_ingested == "Yes", na.rm = TRUE),
sum(df_clinical$firstaid_incision == "Yes", na.rm = TRUE),
sum(!es_faltante(as.character(df_clinical$other)))
)
) %>%
mutate(porcentaje = 100 * casos / n_total)
kable(
tabla_primeros_auxilios,
digits = 1,
caption = "Frecuencia de primeros auxilios registrados"
)| tipo | casos | porcentaje |
|---|---|---|
| Sin primeros auxilios | 286 | 30.7 |
| Torniquete | 508 | 54.5 |
| Vendaje | 16 | 1.7 |
| Hierbas aplicadas | 43 | 4.6 |
| Hierbas ingeridas | 138 | 14.8 |
| Incisión | 83 | 8.9 |
| Otro procedimiento | 36 | 3.9 |
Una persona pudo recibir más de un procedimiento, por lo que estas frecuencias no deben sumarse para obtener un total general. Los espacios en blanco de las variables múltiples tampoco se interpretaron automáticamente como ausencia.
5.6 Distribución temporal
tabla_mensual <- df_clinical %>%
filter(!is.na(date_snakebite)) %>%
mutate(mes = as.Date(format(date_snakebite, "%Y-%m-01"))) %>%
count(mes, name = "casos") %>%
arrange(mes)
mes_maximo <- tabla_mensual %>%
slice_max(casos, n = 1, with_ties = FALSE)
kable(mes_maximo, caption = "Mes con el mayor número de registros en la base pública")| mes | casos |
|---|---|
| 2020-12-01 | 103 |
6. Discrepancias entre el artículo y las bases públicas
presentaciones_indicacion <- c(
"Painful progressive swelling (Cytotoxic)",
"Progressive weakness (Neurotoxic)",
"Bleeding (Hemotoxic)"
)
n_sin_av_con_indicacion <- sum(
df_clinical$av_administered == "No" &
df_clinical$clinical_presentation %in% presentaciones_indicacion,
na.rm = TRUE
)
n_retirados <- sum(df_clinical$occupation_type == 14, na.rm = TRUE)
n_primer_auxilio_no <- sum(df_clinical$firstaid_any == "No", na.rm = TRUE)
n_otro_auxilio <- sum(!es_faltante(as.character(df_clinical$other)))
n_premed_si <- sum(
df_clinical$av_administered == "Yes" &
df_clinical$pre_med_adrenaline == "Yes",
na.rm = TRUE
)
n_premed_no <- sum(
df_clinical$av_administered == "Yes" &
df_clinical$pre_med_adrenaline == "No",
na.rm = TRUE
)
n_premed_faltante <- sum(
df_clinical$av_administered == "Yes" &
is.na(df_clinical$pre_med_adrenaline),
na.rm = TRUE
)
n_reaccion_premed_si <- sum(
df_clinical$av_administered == "Yes" &
df_clinical$pre_med_adrenaline == "Yes" & reaccion_registrada,
na.rm = TRUE
)
n_reaccion_premed_no <- sum(
df_clinical$av_administered == "Yes" &
df_clinical$pre_med_adrenaline == "No" & reaccion_registrada,
na.rm = TRUE
)
discapacidad <- df_clinical$final_outcome == "Permanent disability"
n_fasciotomia_discapacidad <- sum(
discapacidad & df_clinical$fasciotomy == "Yes",
na.rm = TRUE
)
elevacion_disponible <- df_elevation$elevation[geo_completo]
tabla_discrepancias <- tibble(
indicador = c(
"Pacientes menores de 30 años",
"Pacientes retirados",
"RR de 30-39 años e IC95%",
"Pacientes sin primeros auxilios",
"Otros primeros auxilios",
"Sin antiveneno pese a indicación clínica",
"Premedicación entre los 200 tratados",
"Reacción al antiveneno según premedicación",
"Fasciotomías entre pacientes con discapacidad",
"Registros con coordenadas y elevación",
"Rango de elevación",
"Mes de mayor frecuencia"
),
publicado = c(
"55%",
"6%",
"1.68 (0.97-1.40)",
"289 en la Tabla S2",
"39 en la Tabla S2",
"105",
"167 sí y 33 no",
"39% con premedicación y 15% sin premedicación",
"3",
"821",
"71-1499 m",
"Enero"
),
recalculado = c(
sprintf("%d/%d = %.1f%%", n_menores_30, n_total, 100 * n_menores_30 / n_total),
sprintf("%d/%d = %.1f%%", n_retirados, n_total, 100 * n_retirados / n_total),
"El estimador puntual queda fuera del intervalo",
as.character(n_primer_auxilio_no),
as.character(n_otro_auxilio),
as.character(n_sin_av_con_indicacion),
sprintf("%d sí, %d no y %d faltantes", n_premed_si, n_premed_no, n_premed_faltante),
sprintf(
"%.1f%% (%d/%d) frente a %.1f%% (%d/%d); %d sin dato de premedicación",
100 * n_reaccion_premed_si / n_premed_si,
n_reaccion_premed_si,
n_premed_si,
100 * n_reaccion_premed_no / n_premed_no,
n_reaccion_premed_no,
n_premed_no,
n_premed_faltante
),
as.character(n_fasciotomia_discapacidad),
as.character(sum(geo_completo)),
sprintf("%.0f-%.0f m", min(elevacion_disponible), max(elevacion_disponible)),
sprintf("%s (%d casos)", format(mes_maximo$mes, "%B de %Y"), mes_maximo$casos)
),
evaluacion = c(
"Porcentaje incorrecto o denominador no declarado",
"Error en la redacción del artículo",
"Error matemático o tipográfico",
"La cifra publicada produce un total imposible de 935",
"No coincide con la base pública",
"El total 105 excluye los dos registros hemotóxicos mencionados en el texto",
"Los datos faltantes fueron tratados aparentemente como ausencia",
"El grupo sin premedicación reproducible es 27, no 33; el 15% no se reproduce",
"Inconsistencia entre la base y el texto",
"No se reproduce el subconjunto espacial publicado",
"No se reproduce con el archivo público",
"No se reproduce con las fechas públicas"
)
)
kable(
tabla_discrepancias,
caption = "Comparación entre los resultados publicados y los recalculados"
)| indicador | publicado | recalculado | evaluacion |
|---|---|---|---|
| Pacientes menores de 30 años | 55% | 542/932 = 58.2% | Porcentaje incorrecto o denominador no declarado |
| Pacientes retirados | 6% | 33/932 = 3.5% | Error en la redacción del artículo |
| RR de 30-39 años e IC95% | 1.68 (0.97-1.40) | El estimador puntual queda fuera del intervalo | Error matemático o tipográfico |
| Pacientes sin primeros auxilios | 289 en la Tabla S2 | 286 | La cifra publicada produce un total imposible de 935 |
| Otros primeros auxilios | 39 en la Tabla S2 | 36 | No coincide con la base pública |
| Sin antiveneno pese a indicación clínica | 105 | 107 | El total 105 excluye los dos registros hemotóxicos mencionados en el texto |
| Premedicación entre los 200 tratados | 167 sí y 33 no | 167 sí, 27 no y 6 faltantes | Los datos faltantes fueron tratados aparentemente como ausencia |
| Reacción al antiveneno según premedicación | 39% con premedicación y 15% sin premedicación | 91.6% (153/167) frente a 88.9% (24/27); 6 sin dato de premedicación | El grupo sin premedicación reproducible es 27, no 33; el 15% no se reproduce |
| Fasciotomías entre pacientes con discapacidad | 3 | 4 | Inconsistencia entre la base y el texto |
| Registros con coordenadas y elevación | 821 | 862 | No se reproduce el subconjunto espacial publicado |
| Rango de elevación | 71-1499 m | 149-1406 m | No se reproduce con el archivo público |
| Mes de mayor frecuencia | Enero | diciembre de 2020 (103 casos) | No se reproduce con las fechas públicas |
La comparación visual conserva separados los seis pacientes sin información de premedicación. “Sin reacción registrada” significa que las columnas de reacción están vacías; no garantiza que el paciente haya sido evaluado y confirmado como libre de reacción.
7. Evaluación crítica de las pruebas estadísticas
7.1 Reconstrucción de las tablas de contingencia
La Figura 3 muestra diez categorías de circunstancias y nueve sitios anatómicos. Una tabla de 10 por 9 debería tener 72 grados de libertad, pero el artículo informa 64. Para edad y sitio anatómico ocurre una discrepancia similar: la figura muestra diez grupos de edad y nueve sitios, mientras el texto informa solamente 24 grados de libertad.
df_chi <- df_clinical %>%
mutate(
circunstancia_agrupada = case_when(
circumstances_snakebite %in% c(
"Indoors_sleeping", "Outside_sleeping", "Sleeping"
) ~ "Sleeping",
circumstances_snakebite %in% c("Indoors_home", "Indoors") ~ "Indoors",
circumstances_snakebite %in% c(
"Indoors_walking", "Outside_walking", "Walking"
) ~ "Walking",
circumstances_snakebite %in% c(
"Outside_object", "Indoors_object", "Object"
) ~ "Object",
circumstances_snakebite %in% c(
"Outside_homestead", "Outside"
) ~ "Outdoors",
circumstances_snakebite %in% c(
"Outside_homestead_playing", "Outside_playing"
) ~ "Playing",
circumstances_snakebite == "Outside_water" ~ "Water",
circumstances_snakebite == "Outside_toilet" ~ "Toilet",
circumstances_snakebite %in% c(
"Farm_work", "Outside_work", "Indoors_work", "Work"
) ~ "Work",
circumstances_snakebite %in% c("Inside_car", "Misc") ~ "Misc",
TRUE ~ NA_character_
),
edad_figura = case_when(
age_clean %in% c("90-99", "100-109") ~ "90+",
TRUE ~ age_clean
)
) %>%
mutate(
circunstancia_es = case_when(
circunstancia_agrupada == "Sleeping" ~ "Dormir",
circunstancia_agrupada == "Indoors" ~ "Dentro de la vivienda",
circunstancia_agrupada == "Walking" ~ "Caminar",
circunstancia_agrupada == "Object" ~ "Manipular objetos",
circunstancia_agrupada == "Outdoors" ~ "Exterior del hogar",
circunstancia_agrupada == "Playing" ~ "Jugar",
circunstancia_agrupada == "Water" ~ "Actividad acuática",
circunstancia_agrupada == "Toilet" ~ "Uso del sanitario",
circunstancia_agrupada == "Work" ~ "Trabajo",
circunstancia_agrupada == "Misc" ~ "Otros",
TRUE ~ NA_character_
),
posicion_agrupada = case_when(
position_snakebite %in% c("foot", "leg") ~ "Miembro inferior",
position_snakebite %in% c("hand", "arm") ~ "Miembro superior",
position_snakebite %in% c("head", "eye") ~ "Cabeza y ojos",
position_snakebite %in% c("buttock", "torso", "multiple") ~
"Tronco u otros",
TRUE ~ NA_character_
)
)
tabla_circunstancia <- table(
df_chi$circunstancia_agrupada,
df_chi$position_snakebite,
useNA = "no"
)
tabla_edad_sitio <- table(
df_chi$edad_figura,
df_chi$position_snakebite,
useNA = "no"
)
chi_circunstancia <- suppressWarnings(
chisq.test(tabla_circunstancia, correct = FALSE)
)
chi_edad <- suppressWarnings(
chisq.test(tabla_edad_sitio, correct = FALSE)
)
set.seed(20260904)
p_mc_circunstancia <- chisq.test(
tabla_circunstancia,
simulate.p.value = TRUE,
B = 9999
)$p.value
set.seed(20260904)
p_mc_edad <- chisq.test(
tabla_edad_sitio,
simulate.p.value = TRUE,
B = 9999
)$p.value
cramer_v <- function(prueba, tabla) {
sqrt(
as.numeric(prueba$statistic) /
(sum(tabla) * min(nrow(tabla) - 1, ncol(tabla) - 1))
)
}
resumen_chi <- tibble(
asociacion = c("Circunstancia × sitio", "Edad × sitio"),
dimensiones = c(
paste(dim(tabla_circunstancia), collapse = " × "),
paste(dim(tabla_edad_sitio), collapse = " × ")
),
chi2_recalculado = c(
as.numeric(chi_circunstancia$statistic),
as.numeric(chi_edad$statistic)
),
gl_recalculados = c(
as.numeric(chi_circunstancia$parameter),
as.numeric(chi_edad$parameter)
),
gl_publicados = c(64, 24),
porcentaje_esperadas_menor_5 = c(
100 * mean(chi_circunstancia$expected < 5),
100 * mean(chi_edad$expected < 5)
),
cramer_v = c(
cramer_v(chi_circunstancia, tabla_circunstancia),
cramer_v(chi_edad, tabla_edad_sitio)
),
p_monte_carlo = c(p_mc_circunstancia, p_mc_edad)
)
kable(
resumen_chi,
digits = 4,
caption = "Reevaluación de las asociaciones presentadas en la Figura 3"
)| asociacion | dimensiones | chi2_recalculado | gl_recalculados | gl_publicados | porcentaje_esperadas_menor_5 | cramer_v | p_monte_carlo |
|---|---|---|---|---|---|---|---|
| Circunstancia × sitio | 10 × 9 | 315.8023 | 72 | 64 | 75.5556 | 0.2324 | 0.0001 |
| Edad × sitio | 10 × 9 | 131.2191 | 72 | 24 | 74.4444 | 0.1414 | 0.0144 |
Interpretación
La gran proporción de frecuencias esperadas menores de cinco incumple el supuesto clásico del chi-cuadrado. Una alternativa es definir agrupaciones clínicamente justificadas antes del análisis y calcular la significancia por simulación de Monte Carlo. El artículo también debería explicar exactamente qué categorías fueron combinadas para obtener los grados de libertad publicados.
El uso de Kruskal-Wallis y comparaciones de Dunn tampoco se encuentra descrito con suficiente detalle. Estas pruebas no deben aplicarse directamente a variables nominales codificadas con números, porque los códigos no representan una escala cuantitativa.
7.2 Modelos espaciales, agrupamiento y lenguaje de riesgo
evaluacion_modelos <- tribble(
~componente, ~debilidad_observada, ~mejora_necesaria,
"Unidad de análisis", "La base pública contiene pacientes, pero los predictores territoriales pertenecen a áreas", "Publicar una base por área con el conteo de casos, la población y el tiempo de observación",
"Denominador", "Los casos observados por sí solos no permiten calcular incidencia", "Incluir población-tiempo como offset y áreas con cero casos",
"Dependencia espacial", "Numerosos pacientes comparten las mismas coordenadas", "Modelar la agrupación por área y evaluar autocorrelación espacial de los residuos",
"Agrupamiento", "No se documentan con suficiente detalle la elección del algoritmo, el número de grupos y su estabilidad", "Comparar soluciones, justificar el número de grupos y evaluar estabilidad por remuestreo",
"Probabilidad", "Un índice normalizado o una pertenencia a conglomerados se interpreta como probabilidad de riesgo", "Reservar el término probabilidad para predicciones calibradas contra desenlaces observados",
"Desenlace adverso", "No se muestra una validación directa del índice contra la variable final_outcome", "Definir el desenlace, separar entrenamiento y validación y reportar discriminación y calibración",
"MaxEnt", "No se publican todos los detalles de sesgo de muestreo, selección de fondo, ajuste y validación espacial", "Reportar esos parámetros y utilizar validación cruzada espacial",
"Incertidumbre", "Los mapas se presentan sin intervalos ni estabilidad de clasificación", "Agregar intervalos, mapas de incertidumbre y análisis de sensibilidad"
)
kable(
evaluacion_modelos,
caption = "Evaluación de reproducibilidad e interpretación de los modelos espaciales"
)| componente | debilidad_observada | mejora_necesaria |
|---|---|---|
| Unidad de análisis | La base pública contiene pacientes, pero los predictores territoriales pertenecen a áreas | Publicar una base por área con el conteo de casos, la población y el tiempo de observación |
| Denominador | Los casos observados por sí solos no permiten calcular incidencia | Incluir población-tiempo como offset y áreas con cero casos |
| Dependencia espacial | Numerosos pacientes comparten las mismas coordenadas | Modelar la agrupación por área y evaluar autocorrelación espacial de los residuos |
| Agrupamiento | No se documentan con suficiente detalle la elección del algoritmo, el número de grupos y su estabilidad | Comparar soluciones, justificar el número de grupos y evaluar estabilidad por remuestreo |
| Probabilidad | Un índice normalizado o una pertenencia a conglomerados se interpreta como probabilidad de riesgo | Reservar el término probabilidad para predicciones calibradas contra desenlaces observados |
| Desenlace adverso | No se muestra una validación directa del índice contra la variable final_outcome | Definir el desenlace, separar entrenamiento y validación y reportar discriminación y calibración |
| MaxEnt | No se publican todos los detalles de sesgo de muestreo, selección de fondo, ajuste y validación espacial | Reportar esos parámetros y utilizar validación cruzada espacial |
| Incertidumbre | Los mapas se presentan sin intervalos ni estabilidad de clasificación | Agregar intervalos, mapas de incertidumbre y análisis de sensibilidad |
Por tanto, los mapas publicados son útiles para generar hipótesis y priorizar áreas, pero no deberían describirse como probabilidades individuales o tasas de incidencia sin una validación adicional.
Se revisó también el panel denominado “Padidar et al. vs. cohorte local” del material complementario del grupo. No se incorporó porque ambas columnas proceden de los mismos 932 registros y varios porcentajes fueron escritos manualmente. Esa comparación sería circular y no constituye una cohorte independiente ni una validación externa.
8. Reanálisis exploratorio de la elevación
8.1 Por qué el modelo del archivo inicial es incorrecto
El modelo glm.nb(patient ~ elevation) utiliza como
respuesta el número consecutivo del paciente. Esto equivale a preguntar
si los identificadores 1, 2, 3, …, 932 cambian con la elevación. El
identificador no es una frecuencia ni una medida de incidencia, por lo
que su coeficiente y valor p carecen de interpretación
epidemiológica.
8.2 Distribución espacial descriptiva
La base geográfica presenta numerosos pacientes con coordenadas idénticas. Para evitar superponer cientos de puntos, se agruparon los registros por coordenada y elevación. El tamaño del punto representa la cantidad de casos asignados a cada localidad y el color representa su elevación.
La figura muestra concentración de registros en un número reducido de coordenadas. Sin los límites administrativos y la población de cada área no es posible convertir esta distribución en tasas de incidencia ni afirmar que los puntos más grandes representan mayor riesgo poblacional.
8.3 Modelo exploratorio con frecuencias por localidad
La corrección mínima consiste en agrupar los registros que comparten coordenadas y elevación y utilizar como respuesta el número de casos. Este modelo es exploratorio: no estima incidencia porque no incorpora población, tiempo-persona ni localidades con cero casos.
modelo_elevacion <- MASS::glm.nb(
casos ~ I(elevation / 100),
data = geo_localidades
)
coef_altitud <- coef(summary(modelo_elevacion))[2, ]
beta_altitud <- unname(coef_altitud["Estimate"])
se_altitud <- unname(coef_altitud["Std. Error"])
resultado_modelo <- tibble(
localidades = nrow(geo_localidades),
casos_incluidos = sum(geo_localidades$casos),
beta_por_100_m = beta_altitud,
razon_frecuencias = exp(beta_altitud),
cambio_porcentual = 100 * (exp(beta_altitud) - 1),
ic95_inferior = exp(beta_altitud - 1.96 * se_altitud),
ic95_superior = exp(beta_altitud + 1.96 * se_altitud),
valor_p = unname(coef_altitud["Pr(>|z|)"]),
theta = modelo_elevacion$theta,
aic = AIC(modelo_elevacion)
)
kable(
resultado_modelo,
digits = 4,
caption = "Modelo binomial negativo exploratorio por localidad geocodificada"
)| localidades | casos_incluidos | beta_por_100_m | razon_frecuencias | cambio_porcentual | ic95_inferior | ic95_superior | valor_p | theta | aic |
|---|---|---|---|---|---|---|---|---|---|
| 59 | 862 | -0.0641 | 0.9379 | -6.2057 | 0.8943 | 0.9837 | 0.0084 | 3.0766 | 418.3094 |
En este subconjunto, el modelo estima un cambio de -6.2% en la frecuencia registrada por cada 100 m de elevación. Es una asociación entre localidades que ya aportaron casos; no debe interpretarse como efecto causal ni como razón de tasas poblacionales.
8.4 Modelo epidemiológico recomendado
Para estimar incidencia se requeriría una base con todas las áreas censales, incluidas las que registraron cero casos, junto con sus poblaciones y tiempos de observación. La estructura apropiada sería:
modelo_ideal <- MASS::glm.nb(
casos ~ I(elevacion / 100) + pobreza + uso_suelo + acceso_salud +
offset(log(poblacion * anos_observacion)),
data = datos_por_area_censal
)El exponencial del coeficiente de elevación se interpretaría como razón de tasas de incidencia por cada 100 metros, manteniendo constantes las demás variables.
9. Propuestas de mejora
tabla_mejoras <- tribble(
~componente, ~problema, ~propuesta,
"Organización", "No existe identificador compartido", "Incluir un Patient_ID anónimo en ambas bases",
"Diccionario", "Los códigos de ocupación no están explicados", "Publicar etiquetas, unidades, categorías y códigos faltantes",
"Calidad", "Se confunden ausencia, no aplicación y dato no registrado", "Diferenciar No, No aplica y Sin dato",
"Descriptivos", "Algunos porcentajes no declaran denominador", "Mostrar n/N y porcentaje en cada resultado",
"Incidencia", "Se analizan casos sin población expuesta", "Usar población-tiempo como offset",
"Dependencia espacial", "Muchas observaciones comparten coordenadas", "Agrupar por área y utilizar modelos espaciales o multinivel",
"Chi-cuadrado", "Existen numerosas frecuencias esperadas pequeñas", "Agrupar categorías a priori o usar simulación de Monte Carlo",
"Selección de variables", "El procedimiento escalonado no se describe", "Definir un modelo completo basado en hipótesis y comparar sensibilidad",
"Mapas", "Los índices normalizados se denominan probabilidades", "Denominarlos perfiles relativos o validar un modelo predictivo",
"Reproducibilidad", "No se publicaron scripts ni productos espaciales intermedios", "Compartir código, semillas, versiones y base por área censal"
)
kable(tabla_mejoras, caption = "Matriz de problemas y propuestas de mejora")| componente | problema | propuesta |
|---|---|---|
| Organización | No existe identificador compartido | Incluir un Patient_ID anónimo en ambas bases |
| Diccionario | Los códigos de ocupación no están explicados | Publicar etiquetas, unidades, categorías y códigos faltantes |
| Calidad | Se confunden ausencia, no aplicación y dato no registrado | Diferenciar No, No aplica y Sin dato |
| Descriptivos | Algunos porcentajes no declaran denominador | Mostrar n/N y porcentaje en cada resultado |
| Incidencia | Se analizan casos sin población expuesta | Usar población-tiempo como offset |
| Dependencia espacial | Muchas observaciones comparten coordenadas | Agrupar por área y utilizar modelos espaciales o multinivel |
| Chi-cuadrado | Existen numerosas frecuencias esperadas pequeñas | Agrupar categorías a priori o usar simulación de Monte Carlo |
| Selección de variables | El procedimiento escalonado no se describe | Definir un modelo completo basado en hipótesis y comparar sensibilidad |
| Mapas | Los índices normalizados se denominan probabilidades | Denominarlos perfiles relativos o validar un modelo predictivo |
| Reproducibilidad | No se publicaron scripts ni productos espaciales intermedios | Compartir código, semillas, versiones y base por área censal |
10. Conclusiones
Los principales conteos demográficos y clínicos pueden reproducirse, pero se identificaron discrepancias numéricas, categorías inconsistentes y porcentajes con denominadores poco claros. Las limitaciones más importantes se concentran en la vinculación de las bases, el manejo de datos faltantes, las pruebas con tablas dispersas y la falta de información necesaria para reproducir el modelo de incidencia publicado.
El agrupamiento espacial realizado por los autores debe interpretarse como una clasificación exploratoria de perfiles territoriales y no como una probabilidad calibrada de mordedura o mal desenlace. La recomendación es reconstruir la base por área censal, incorporar la población expuesta, validar los modelos espacialmente y presentar intervalos de incertidumbre.
La conclusión general no es que todos los resultados del artículo sean incorrectos. La descripción del problema sanitario conserva valor, pero la trazabilidad entre los datos, los análisis y algunas interpretaciones necesita mejorarse.
Referencias principales
- Benchimol, E. I., et al. (2015). The REporting of studies Conducted using Observational Routinely-collected health Data (RECORD) statement. PLOS Medicine, 12(10), e1001885. https://doi.org/10.1371/journal.pmed.1001885
- Padidar, S., et al. (2023). Snakebite epidemiology, outcomes and multi-cluster risk modelling in Eswatini. PLOS Neglected Tropical Diseases, 17(11), e0011732. https://doi.org/10.1371/journal.pntd.0011732
- von Elm, E., et al. (2007). The Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) statement. Annals of Internal Medicine, 147(8), 573-577. https://doi.org/10.7326/0003-4819-147-8-200710160-00010
Información de la sesión de R
sessionInfo()## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=Spanish_Panama.utf8 LC_CTYPE=Spanish_Panama.utf8
## [3] LC_MONETARY=Spanish_Panama.utf8 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Panama.utf8
##
## time zone: America/Panama
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] scales_1.4.0 knitr_1.51 ggplot2_4.0.3 stringr_1.6.0 tidyr_1.3.2
## [6] dplyr_1.2.1 readxl_1.5.0
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 jsonlite_2.0.0 compiler_4.6.1 tidyselect_1.2.1
## [5] dichromat_2.0-1 jquerylib_0.1.4 yaml_2.3.12 fastmap_1.2.0
## [9] R6_2.6.1 labeling_0.4.3 generics_0.1.4 MASS_7.3-66
## [13] tibble_3.3.1 RColorBrewer_1.1-3 bslib_0.12.0 pillar_1.11.1
## [17] rlang_1.3.0 cachem_1.1.0 stringi_1.8.9 xfun_0.60
## [21] S7_0.2.2 sass_0.4.10 otel_0.2.0 cli_3.6.6
## [25] withr_3.0.3 magrittr_2.0.5 digest_0.6.39 grid_4.6.1
## [29] rstudioapi_0.19.0 lifecycle_1.0.5 vctrs_0.7.3 evaluate_1.0.5
## [33] glue_1.8.1 farver_2.1.2 cellranger_1.1.0 stats4_4.6.1
## [37] rmarkdown_2.32 purrr_1.2.2 tools_4.6.1 pkgconfig_2.0.3
## [41] htmltools_0.5.9