Análisis Bioestadístico Aplicado a la Toxinología

Integrantes

Descripción de la base de datos

La base Toxinologia_Panama contiene 500,000 casos e integra variables demográficas, clínicas y toxinológicas. Fue construida con una semilla reproducible (set.seed(2026)) y se utiliza exclusivamente con fines académicos.

Preparación de los datos

Todos los análisis parten del mismo objeto generado por BASE.R. La preparación se realizó de forma reproducible antes de cada prueba: selección de las variables necesarias, exclusión únicamente de valores perdidos cuando correspondía, definición explícita de factores y categorías de referencia, y verificación de la estructura de los casos.

Planteamiento

Objetivo: comparar el nivel medio del biomarcador entre dos grupos independientes definidos por sexo.

Hipótesis de investigación: el nivel medio del biomarcador difiere entre pacientes femeninos y masculinos.

H₀: μ_F = μ_M.

H_A: μ_F ≠ μ_M.

Justificación: se compara una variable cuantitativa continua entre dos grupos independientes. La t de Student es el análisis principal cuando la igualdad de varianzas es razonable; Welch se usa como sensibilidad.

Estadística descriptiva

Las medias fueron prácticamente idénticas: 74.79 U/L en mujeres y 74.88 U/L en hombres. Las desviaciones estándar también fueron casi iguales (23.95 y 23.86 U/L).

Figura 1. Media del biomarcador e intervalo de confianza del 95 % según sexo.

Evaluación de supuestos

Ver gráfico Q-Q de los residuos

Diagnóstico complementario de la distribución de los residuos.

Resultados

Interpretación y conclusión

Interpretación estadística. No se obtuvo evidencia estadística suficiente para afirmar una diferencia en el nivel medio del biomarcador entre pacientes femeninos y masculinos. La estimación está muy próxima a 0, el IC95 % incluye 0 y el tamaño del efecto es prácticamente nulo. Welch produjo prácticamente la misma inferencia.

Interpretación toxinológica. En esta comparación no ajustada, el sexo no mostró una asociación relevante con el nivel del biomarcador.

Código R reproducible

# Datos
datos_t <- Toxinologia_Panama |>
  dplyr::select(ID, Sexo, Biomarcador_UL) |>
  dplyr::filter(!is.na(Sexo), !is.na(Biomarcador_UL)) |>
  dplyr::mutate(Sexo = factor(Sexo, levels = c("Femenino", "Masculino")))

# Descriptivos
tabla_t_desc <- datos_t |>
  dplyr::group_by(Sexo) |>
  dplyr::summarise(
    n = dplyr::n(),
    Media = mean(Biomarcador_UL),
    DE = sd(Biomarcador_UL),
    Mediana = median(Biomarcador_UL),
    IQR = IQR(Biomarcador_UL),
    .groups = "drop"
  )

# Residuos para diagnóstico
modelo_t <- lm(Biomarcador_UL ~ Sexo, data = datos_t)

# Brown-Forsythe / Levene centrado en mediana
car::leveneTest(
  Biomarcador_UL ~ Sexo,
  data = datos_t,
  center = median
)

# Student: análisis principal
resultado_t_student <- t.test(
  Biomarcador_UL ~ Sexo,
  data = datos_t,
  var.equal = TRUE,
  alternative = "two.sided",
  conf.level = 0.95
)

# Welch: sensibilidad
resultado_t_welch <- t.test(
  Biomarcador_UL ~ Sexo,
  data = datos_t,
  var.equal = FALSE,
  alternative = "two.sided",
  conf.level = 0.95
)

# Tamaño del efecto
effectsize::cohens_d(
  Biomarcador_UL ~ Sexo,
  data = datos_t,
  pooled_sd = TRUE,
  ci = 0.95
)

Planteamiento

Objetivo: evaluar si la distribución de la severidad clínica difiere según el agente toxinológico.

H₀: Agente y Severidad son independientes.

H_A: Agente y Severidad no son independientes.

Justificación: ambas variables son categóricas y la prueba χ² de independencia permite comparar las frecuencias observadas con las esperadas bajo independencia. La naturaleza ordinal de Severidad se reconoce, aunque χ² no utiliza el orden de las categorías.

Estadística descriptiva y visualización

Evaluación de supuestos

Resultados

Interpretación y conclusión

Interpretación estadística. La prueba χ² no mostró evidencia suficiente de asociación entre agente y severidad, χ²(15) = 19.19, p = 0.205. Además, V de Cramér = 0.0036 indica una asociación de magnitud prácticamente nula.

Interpretación toxinológica. Las proporciones de casos leves, moderados, severos y críticos fueron muy semejantes entre los seis agentes. Las diferencias máximas entre agentes fueron inferiores a 0.7 puntos porcentuales en cada nivel de severidad.

Código R reproducible

# Datos
datos_chi <- Toxinologia_Panama |>
  dplyr::select(ID, Agente, Severidad) |>
  dplyr::filter(!is.na(Agente), !is.na(Severidad)) |>
  dplyr::mutate(
    Agente = factor(Agente),
    Severidad = factor(
      Severidad,
      levels = c("Leve", "Moderada", "Severa", "Crítica"),
      ordered = TRUE
    )
  )

# Tabla y prueba
tabla_chi <- table(datos_chi$Agente, datos_chi$Severidad)
resultado_chi <- chisq.test(tabla_chi, correct = FALSE)

# Tamaño de asociación
chi2 <- unname(resultado_chi$statistic)
n <- sum(tabla_chi)
V_cramer <- sqrt(
  chi2 / (n * min(nrow(tabla_chi)-1, ncol(tabla_chi)-1))
)

# Proporciones por agente
100 * prop.table(tabla_chi, margin = 1)

# Gráfico principal: orden clínico Leve → Crítica
ggplot2::ggplot(datos_chi, ggplot2::aes(x = Agente, fill = Severidad)) +
  ggplot2::geom_bar(
    position = ggplot2::position_fill(reverse = TRUE)
  ) +
  ggplot2::scale_y_continuous(labels = scales::percent_format()) +
  ggplot2::labs(
    title = "Distribución de la severidad según agente toxinológico",
    x = "Agente toxinológico",
    y = "Porcentaje de casos",
    fill = "Severidad"
  ) +
  ggplot2::theme_minimal(base_size = 12)

# Residuos estandarizados: diagnóstico exploratorio
stdres_chi <- resultado_chi$stdres

Planteamiento

Objetivo

Objetivo. Evaluar si el nivel medio del biomarcador toxinológico difiere entre los distintos agentes registrados y, si el contraste global resulta significativo, identificar mediante Tukey qué pares de agentes presentan diferencias.

Hipótesis de investigación: el nivel medio del biomarcador toxinológico varía según el agente involucrado.

H₀: μ₁ = μ₂ = μ₃ = μ₄ = μ₅ = μ₆.

H_A: al menos una de las medias es diferente.

Nivel de significancia: α = 0.05.

Justificación: el ANOVA de un factor es apropiado porque se compara una variable cuantitativa continua entre más de dos grupos independientes definidos por un único factor categórico. Tukey se utiliza después del contraste global para localizar las diferencias controlando el error por comparaciones múltiples.

Estadística descriptiva y visualización

Descripción principal. Bothrops asper presentó la media más alta del biomarcador (80.52 U/L), seguido de Crotalus simus (78.63 U/L). Hymenoptera presentó la media más baja (68.54 U/L). Las desviaciones estándar fueron muy similares entre los seis grupos (23.36–23.55 U/L).

Figura 1. Distribución del biomarcador por agente. La muestra de 3,000 casos por agente se utilizó únicamente para facilitar la visualización; el ANOVA se realizó con los 500,000 casos.

Evaluación de supuestos

Figura 2. Gráfico Q-Q de una muestra reproducible de 5,000 residuos.

Resultado del ANOVA

Comparaciones post hoc de Tukey

El ANOVA global fue significativo, por lo que se aplicó Tukey HSD. Las 15 de 15 comparaciones resultaron estadísticamente significativas (p ajustada < 0.001).

Figura 3. Diferencias de medias e IC95 % de las 15 comparaciones de Tukey. Ningún intervalo cruza 0.

Interpretación y conclusión

Interpretación estadística

Se encontró evidencia estadísticamente significativa de diferencias en el nivel medio del biomarcador entre los seis agentes toxinológicos. Bothrops asper presentó la mayor media (80.52 U/L) y Hymenoptera la menor (68.54 U/L). Tukey confirmó diferencias entre todos los pares de agentes. No obstante, el tamaño del efecto global fue modesto (η² = 0.040), por lo que las diferencias deben valorarse en conjunto con sus magnitudes en U/L y no únicamente mediante los valores p.

Interpretación toxinológica

El nivel del biomarcador no fue uniforme entre los distintos agentes, lo que sugiere que el tipo de exposición está asociado con diferencias en la respuesta biológica representada por este marcador. Sin embargo, el ANOVA de un factor es un análisis no ajustado: no controla otras variables disponibles, como dosis, edad o tratamiento. En consecuencia, las diferencias observadas no deben interpretarse como efectos causales independientes del agente.

Código R reproducible

El código se mantiene dentro del documento para cumplir con el requisito de reproducibilidad, pero permanece plegado para no recargar la lectura del informe.

# ============================================================
# ASIGNACIÓN 3 - GRUPO 1
# ANOVA DE UN FACTOR + TUKEY
# ============================================================

library(dplyr)
library(ggplot2)
library(car)
library(effectsize)
library(knitr)

# BASE.R debe estar en la misma carpeta
source("BASE.R")

# 1. Preparación
datos_anova <- Toxinologia_Panama %>%
  select(ID, Agente, Biomarcador_UL) %>%
  filter(!is.na(Agente), !is.na(Biomarcador_UL)) %>%
  mutate(Agente = droplevels(factor(Agente)))

# 2. Descriptivos
tabla_descriptiva <- datos_anova %>%
  group_by(Agente) %>%
  summarise(
    n = n(),
    Media = mean(Biomarcador_UL),
    DE = sd(Biomarcador_UL),
    Mediana = median(Biomarcador_UL),
    IQR = IQR(Biomarcador_UL),
    .groups = "drop"
  )

# 3. ANOVA
modelo_anova <- aov(
  Biomarcador_UL ~ Agente,
  data = datos_anova
)

summary(modelo_anova)

# 4. Supuestos
set.seed(2026)
residuos_anova <- residuals(modelo_anova)

muestra_residuos <- sample(
  residuos_anova,
  size = min(5000, length(residuos_anova))
)

shapiro.test(muestra_residuos)

car::leveneTest(
  Biomarcador_UL ~ Agente,
  data = datos_anova,
  center = median
)

# 5. Tamaño del efecto
effectsize::eta_squared(
  modelo_anova,
  partial = FALSE,
  ci = 0.95
)

effectsize::omega_squared(
  modelo_anova,
  partial = FALSE,
  ci = 0.95
)

# 6. Tukey
tukey_anova <- TukeyHSD(
  modelo_anova,
  "Agente",
  conf.level = 0.95
)

tukey_anova

# 7. Visualización descriptiva
set.seed(2026)

datos_grafico <- datos_anova %>%
  group_by(Agente) %>%
  slice_sample(n = 3000) %>%
  ungroup()

ggplot(
  datos_grafico,
  aes(x = Agente, y = Biomarcador_UL, fill = Agente)
) +
  geom_violin(alpha = 0.55, trim = FALSE, color = NA) +
  geom_boxplot(width = 0.14, alpha = 0.82, outlier.shape = NA) +
  labs(
    title = "Distribución del biomarcador según agente toxinológico",
    x = "Agente toxinológico",
    y = "Biomarcador (U/L)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "none",
    axis.text.x = element_text(angle = 28, hjust = 1)
  )

# 8. Q-Q de residuos
diagnostico_residuos <- data.frame(Residuo = muestra_residuos)

ggplot(
  diagnostico_residuos,
  aes(sample = Residuo)
) +
  stat_qq(alpha = 0.35) +
  stat_qq_line(linewidth = 0.9) +
  labs(
    title = "Diagnóstico de normalidad de los residuos",
    x = "Cuantiles teóricos",
    y = "Cuantiles observados"
  ) +
  theme_minimal(base_size = 12)

Planteamiento

Objetivo

Evaluar la dirección y magnitud de la relación entre la edad del paciente y el nivel del biomarcador toxinológico.

Hipótesis de investigación: existe una relación entre la edad del paciente y el nivel del biomarcador toxinológico.

H₀: ρ = 0; no existe correlación poblacional.

H_A: ρ ≠ 0; existe correlación poblacional.

Nivel de significancia: α = 0.05.

Justificación: ambas variables son cuantitativas y el propósito es estimar la dirección y fuerza de su asociación. La elección entre Pearson y Spearman se fundamenta principalmente en la forma de la relación, su linealidad y la posible influencia de valores extremos; no se decide únicamente mediante una prueba formal de normalidad.

Estadística descriptiva y visualización

Figura 4. Diagrama de dispersión de edad y biomarcador. Se utilizaron 12,000 casos únicamente para facilitar la visualización; el cálculo de la correlación empleó los 500,000 registros.

El diagrama muestra una pendiente positiva leve, acompañada de una dispersión individual amplia. Visualmente, la relación parece aproximadamente lineal, pero de baja intensidad.

Figura 5. Media del biomarcador por edad. El patrón de las medias es aproximadamente lineal y creciente.

Selección entre Pearson y Spearman

Resultado de la correlación de Pearson

Comprobación de robustez: Spearman

La diferencia absoluta entre ambos coeficientes es de apenas 0.0004. Esto respalda que la conclusión general es estable frente al uso de una medida basada en rangos y que los valores extremos observados no están alterando de forma importante el resultado.

Interpretación y conclusión

Interpretación estadística

Existe evidencia estadística de una relación lineal positiva entre la edad y el biomarcador, pero su magnitud es muy débil (r = 0.0870). El IC95 % es estrecho (0.0843–0.0898) debido al gran tamaño muestral, mientras que r² = 0.0076 indica que un modelo lineal simple con la edad como único predictor explica menos del 1 % de la variabilidad del biomarcador. Spearman produjo un resultado prácticamente idéntico (ρ = 0.0866), reforzando la estabilidad de la conclusión.

Interpretación toxinológica

El patrón observado indica que los valores del biomarcador tienden a aumentar ligeramente a medida que aumenta la edad. Sin embargo, la asociación es muy débil, por lo que la edad por sí sola presenta poca capacidad para describir las diferencias individuales del biomarcador. Este análisis es bivariado y no controla otros factores disponibles, como dosis, agente toxinológico o tratamiento; por ello, no permite establecer una relación causal.

Código R reproducible

El código se mantiene dentro del documento para cumplir con el requisito de reproducibilidad, pero permanece plegado para no recargar la lectura.

# ============================================================
# ANÁLISIS 4 - CORRELACIÓN DE PEARSON / SPEARMAN
# Edad vs biomarcador toxinológico
# ============================================================

library(dplyr)
library(ggplot2)
library(knitr)

source("BASE.R")

datos_cor <- Toxinologia_Panama %>%
  select(ID, Edad, Biomarcador_UL) %>%
  filter(!is.na(Edad), !is.na(Biomarcador_UL))

# Estadística descriptiva
tabla_descriptiva <- data.frame(
  Variable = c("Edad (años)", "Biomarcador (U/L)"),
  n = c(nrow(datos_cor), nrow(datos_cor)),
  Media = c(mean(datos_cor$Edad), mean(datos_cor$Biomarcador_UL)),
  DE = c(sd(datos_cor$Edad), sd(datos_cor$Biomarcador_UL)),
  Mediana = c(median(datos_cor$Edad), median(datos_cor$Biomarcador_UL)),
  IQR = c(IQR(datos_cor$Edad), IQR(datos_cor$Biomarcador_UL)),
  Minimo = c(min(datos_cor$Edad), min(datos_cor$Biomarcador_UL)),
  Maximo = c(max(datos_cor$Edad), max(datos_cor$Biomarcador_UL))
)

# Gráfico de dispersión (submuestra solo para visualización)
set.seed(2026)
datos_grafico <- datos_cor %>%
  slice_sample(n = min(12000, nrow(datos_cor)))

ggplot(datos_grafico, aes(x = Edad, y = Biomarcador_UL)) +
  geom_point(alpha = 0.16) +
  geom_smooth(method = "lm", formula = y ~ x, se = TRUE) +
  labs(
    title = "Relación entre edad y biomarcador toxinológico",
    x = "Edad (años)",
    y = "Biomarcador (U/L)"
  ) +
  theme_minimal(base_size = 12)

# Medias del biomarcador por edad para revisar la forma de la relación
tabla_edad <- datos_cor %>%
  group_by(Edad) %>%
  summarise(
    n = n(),
    Media_biomarcador = mean(Biomarcador_UL),
    .groups = "drop"
  )

ggplot(tabla_edad, aes(x = Edad, y = Media_biomarcador)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ x, se = TRUE) +
  labs(
    title = "Media del biomarcador según edad",
    x = "Edad (años)",
    y = "Media del biomarcador (U/L)"
  ) +
  theme_minimal(base_size = 12)

# Valores extremos del biomarcador
Q1_bio <- quantile(datos_cor$Biomarcador_UL, 0.25)
Q3_bio <- quantile(datos_cor$Biomarcador_UL, 0.75)
IQR_bio <- IQR(datos_cor$Biomarcador_UL)

lim_inf <- Q1_bio - 1.5 * IQR_bio
lim_sup <- Q3_bio + 1.5 * IQR_bio

sum(
  datos_cor$Biomarcador_UL < lim_inf |
  datos_cor$Biomarcador_UL > lim_sup
)

# Pearson: análisis principal
pearson_resultado <- cor.test(
  datos_cor$Edad,
  datos_cor$Biomarcador_UL,
  method = "pearson",
  alternative = "two.sided",
  conf.level = 0.95
)

pearson_resultado

r_pearson <- unname(pearson_resultado$estimate)
r2_pearson <- r_pearson^2
r2_pearson

# Spearman: comprobación complementaria
spearman_resultado <- cor.test(
  datos_cor$Edad,
  datos_cor$Biomarcador_UL,
  method = "spearman",
  alternative = "two.sided",
  exact = FALSE
)

spearman_resultado

Planteamiento

Objetivo

Estimar la asociación ajustada entre la dosis estimada y el nivel del biomarcador toxinológico, considerando simultáneamente la edad, el sexo, el agente toxinológico y el tratamiento administrado.

Hipótesis de investigación: la dosis estimada se relaciona con el nivel del biomarcador aun después de considerar simultáneamente las demás variables incluidas en el modelo.

H₀ principal: β_Dosis = 0.

H_A principal: β_Dosis ≠ 0.

Hipótesis global: H₀: todos los coeficientes de pendiente son 0; H_A: al menos uno es distinto de 0.

Nivel de significancia: α = 0.05.

Justificación: la regresión lineal múltiple es apropiada porque la variable respuesta es cuantitativa continua y se desea estudiar simultáneamente predictores cuantitativos y categóricos. El modelo permite estimar la asociación de la dosis con el biomarcador manteniendo constantes las demás variables y cuantificar la capacidad explicativa global mediante R² y R² ajustado.

Estadística descriptiva y visualización

Figura 6. Relación entre dosis y biomarcador. Se utilizaron 12,000 casos únicamente para facilitar la visualización; el modelo de regresión empleó los 500,000 registros.

El gráfico muestra una asociación positiva clara entre la dosis y el biomarcador. La regresión múltiple permitirá estimar esta relación después de controlar simultáneamente las demás variables incluidas en el modelo.

Evaluación de supuestos

Figura 7. Residuos frente a valores ajustados sobre una muestra reproducible de 6,000 residuos.

Figura 8. Gráfico Q-Q de una muestra reproducible de 6,000 residuos.

Modelo y coeficientes ajustados

Coeficientes

Comprobación complementaria

Como comprobación secundaria se calcularon errores estándar robustos HC3. Los coeficientes y la inferencia fueron prácticamente idénticos a los resultados convencionales, por lo que no se modifican las conclusiones principales.

Interpretación y conclusión

Interpretación estadística

El modelo de regresión lineal múltiple explicó aproximadamente el 60.58 % de la variabilidad observada del biomarcador y fue globalmente significativo. Para la pregunta principal, la dosis mantuvo una asociación positiva y precisa después del ajuste: β = 7.007 U/L por mg/kg, IC95 % 6.989–7.024, p < 0.001. La edad, el agente y el tratamiento también presentaron asociaciones ajustadas, mientras que el sexo no mostró evidencia de una contribución independiente una vez consideradas las demás variables.

Interpretación toxinológica

Después de considerar simultáneamente edad, sexo, agente toxinológico y tratamiento, una mayor dosis se asocia con niveles mayores del biomarcador. La magnitud estimada indica un aumento promedio cercano a 7 U/L por cada incremento de 1 mg/kg, manteniendo constantes las demás variables. Asimismo, el agente y el tratamiento presentan diferencias ajustadas relevantes en el nivel del biomarcador. Estas estimaciones describen asociaciones dentro del modelo y no deben interpretarse automáticamente como efectos causales.

Código R reproducible

El código se mantiene dentro del documento para cumplir con el requisito de reproducibilidad, pero permanece plegado para no recargar la lectura del informe.

# ============================================================
# ANÁLISIS 5 - REGRESIÓN LINEAL MÚLTIPLE
# Biomarcador ~ Dosis + Edad + Sexo + Agente + Tratamiento
# ============================================================

library(dplyr)
library(ggplot2)
library(car)
library(lmtest)
library(sandwich)

source("BASE.R")

datos_reg <- Toxinologia_Panama %>%
  select(
    ID, Biomarcador_UL, Dosis_mgkg, Edad,
    Sexo, Agente, Tratamiento
  ) %>%
  filter(complete.cases(.)) %>%
  mutate(
    Sexo = relevel(factor(Sexo), ref = "Femenino"),
    Agente = relevel(factor(Agente), ref = "Animal marino"),
    Tratamiento = relevel(factor(Tratamiento), ref = "Soporte")
  )

# Modelo principal
modelo_reg <- lm(
  Biomarcador_UL ~
    Dosis_mgkg +
    Edad +
    Sexo +
    Agente +
    Tratamiento,
  data = datos_reg
)

summary(modelo_reg)
confint(modelo_reg, level = 0.95)

# Diagnóstico de homocedasticidad
lmtest::bptest(modelo_reg)

# Multicolinealidad
car::vif(modelo_reg)

# Observaciones influyentes
cook <- cooks.distance(modelo_reg)
max(cook)
sum(cook > 4 / nrow(datos_reg))
100 * mean(cook > 4 / nrow(datos_reg))

# Errores estándar robustos HC3
lmtest::coeftest(
  modelo_reg,
  vcov. = sandwich::vcovHC(
    modelo_reg,
    type = "HC3"
  )
)

# Visualización principal
set.seed(2026)
datos_grafico <- datos_reg %>%
  slice_sample(n = min(12000, nrow(datos_reg)))

ggplot(
  datos_grafico,
  aes(x = Dosis_mgkg, y = Biomarcador_UL)
) +
  geom_point(alpha = 0.16) +
  geom_smooth(method = "lm", formula = y ~ x, se = TRUE) +
  labs(
    title = "Relación entre dosis y biomarcador toxinológico",
    x = "Dosis estimada (mg/kg)",
    y = "Biomarcador (U/L)"
  ) +
  theme_minimal(base_size = 12)

# Diagnóstico gráfico de residuos
set.seed(2026)
residuos_reg <- residuals(modelo_reg)
ajustados_reg <- fitted(modelo_reg)
indices_diag <- sample(
  seq_along(residuos_reg),
  size = min(6000, length(residuos_reg))
)

datos_diag <- data.frame(
  Ajustado = ajustados_reg[indices_diag],
  Residuo = residuos_reg[indices_diag]
)

ggplot(datos_diag, aes(x = Ajustado, y = Residuo)) +
  geom_point(alpha = 0.22) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_smooth(method = "loess", formula = y ~ x, se = FALSE) +
  theme_minimal(base_size = 12)

ggplot(datos_diag, aes(sample = Residuo)) +
  stat_qq(alpha = 0.35) +
  stat_qq_line() +
  theme_minimal(base_size = 12)

Conclusiones generales

Las cinco pruebas responden preguntas distintas y complementarias. Los resultados muestran que la significancia estadística debe interpretarse siempre junto con la magnitud del efecto o de la asociación.