packages <- c("tidyverse", "survival", "here", "arm", "tibble", "broom",
              "dagitty", "ggdag", "knitr", "patchwork", "ggridges",
              "gmodels", "xaringanthemer", "latex2exp", "stdReg",
              "epiR", "stats", "epitools",
              "kableExtra", "scales")

invisible(lapply(packages, function(p) {
  suppressMessages(require(p, character.only = TRUE,
                           quietly = TRUE, warn.conflicts = FALSE))
}))

1 Estimación de Efectividad de Vacunas contra COVID-19 en el Reino Unido

Procedo a tomar como referencia el reporte de la Agencia de Seguridad Sanitaria del Reino Unido titulado “COVID-19 Vaccine Surveillance Report”, Semana 41, y la Figura 2(a) del mismo, de esta manera, comienzo por reconstruir dicha gráfica con los datos de la Tabla 2 del reporte.

Esta visualización me permite observar las proporciones de incidencia de COVID-19 por cada 100,000 habitantes, estratificadas por grupo etario y estado vacunal, durante las semanas 37 a 40 de 2021.

datos_casos <- tibble(
  grupo = factor(
    c("< 18", "18-29", "30-39", "40-49", "50-59", "60-69", "70-79", "80+"),
    levels = c("< 18", "18-29", "30-39", "40-49",
               "50-59", "60-69", "70-79", "80+")
  ),
  vacunados    = c(276.5, 402.6, 823.9, 1455.8, 903.1, 589.0, 451.5, 364.6),
  no_vacunados = c(2670.7, 605.0, 709.8, 696.2, 489.3, 314.1, 253.0, 298.5)
)
datos_casos %>%
  pivot_longer(
    cols      = c(vacunados, no_vacunados),
    names_to  = "estado",
    values_to = "tasa"
  ) %>%
  mutate(
    estado = ifelse(estado == "vacunados",
                    "Vacunados (2 dosis o mas)",
                    "No vacunados")
  ) %>%
  ggplot(aes(x = grupo, y = tasa, fill = estado)) +
  geom_col(position = position_dodge(0.7), width = 0.6, alpha = 0.85) +
  geom_text(
    aes(label = format(tasa, big.mark = ",")),
    position = position_dodge(0.7),
    vjust = -0.4, size = 2.5, fontface = "bold"
  ) +
  scale_fill_manual(
    values = c("Vacunados (2 dosis o mas)" = col_vac,
               "No vacunados" = col_nv),
    name = "Estado vacunal"
  ) +
  scale_y_continuous(
    labels = comma_format(),
    expand = expansion(mult = c(0, 0.15))
  ) +
  labs(
    title    = "Proporcion de incidencia de COVID-19 por estado vacunal",
    subtitle = "Semanas 37 a 40 de 2021 - Reporte UKHSA, Semana 41",
    x        = "Grupo etario",
    y        = "Tasa por 100,000 habitantes"
  ) +
  tema_pro
Proporción de incidencia de COVID-19 por estado vacunal y grupo etario, semanas 37 a 40 de 2021. Fuente: COVID-19 Vaccine Surveillance Report, Week 41, UKHSA.

Proporción de incidencia de COVID-19 por estado vacunal y grupo etario, semanas 37 a 40 de 2021. Fuente: COVID-19 Vaccine Surveillance Report, Week 41, UKHSA.

1.1 Pregunta 1a. Medida de frecuencia usada por los autores

La medida de frecuencia que emplean los autores en la Figura 2(a) es una proporción de incidencia, también denominada incidencia acumulada, expresada por cada 100,000 habitantes.

Esta medida se obtiene dividiendo el número de casos nuevos de COVID-19 detectados durante un período determinado entre la población en riesgo al inicio de dicho período, y multiplicando el resultado por 100,000.

En este caso, el período de observación corresponde a cuatro semanas epidemiológicas (semanas 37 a 40 de 2021) y el denominador es la población total de cada grupo, ya sea vacunados con al menos dos dosis o no vacunados (extraída del Sistema Nacional de Gestión de Inmunización del Reino Unido, conocido como NIMS, por sus siglas en inglés).

\[ \text{Proporción de incidencia} = \frac{\text{Casos nuevos en el período}}{\text{Población en riesgo al inicio}} \times 100{,}000 \]

Si bien el reporte utiliza el término “rate” (tasa en inglés), considero que en sentido estricto no se trata de una tasa de incidencia, la cual requeriría persona-tiempo en el denominador, es decir, la sumatoria del tiempo que cada individuo estuvo en riesgo de enfermar; sino que se trata de una proporción de incidencia, dado que el denominador es simplemente el número total de personas en cada grupo al inicio del período de observación.

Ahora bien, cuando el período de seguimiento es corto (como estas cuatro semanas) y la enfermedad es relativamente frecuente, la proporción de incidencia y la tasa de incidencia tienden a ser numéricamente cercanas, de modo que esta distinción, aunque relevante desde el rigor metodológico, considero que no invalida las conclusiones del reporte.

1.2 Pregunta 1b. Qué indican las “tasas” de casos según el estado de inmunización

Al examinar la gráfica que reconstruí, las proporciones de incidencia revelan un patrón aparentemente paradójico, pues en los menores de 18 años, la tasa en no vacunados alcanza 2,670.7 por 100,000 frente a apenas 276.5 en vacunados, lo cual considero que es una diferencia abrumadora que resulta consistente con la expectativa de que la vacunación protege contra la infección.

En el grupo de 18 a 29 años, la relación se mantiene en la misma dirección, aunque con una brecha menor (605.0 en no vacunados frente a 402.6 en vacunados).

Sin embargo, a partir de los 30 años la relación se invierte por completo, pues las tasas de casos son más altas en personas vacunadas que en no vacunadas.

En el grupo de 40 a 49 años, por ejemplo, la tasa en vacunados alcanza 1,455.8 por 100,000 frente a 696.2 en no vacunados, es decir, los vacunados presentan aproximadamente el doble de casos detectados que los no vacunados.

Considero que esta inversión no refleja una falla de la vacuna, sino que es producto de varios fenómenos convergentes que distorsionan la comparabilidad de los grupos, pues en primer lugar, considero que existe un sesgo de testeo diferencial, porque las personas vacunadas en edad laboral activa (entre 30 y 59 años) se sometían a pruebas diagnósticas con mayor frecuencia, ya fuera por requisitos laborales, para viajar o para acceder a eventos masivos.

Esta mayor frecuencia de testeo incrementa la detección de infecciones leves o asintomáticas en el grupo vacunado, inflando artificialmente su proporción de incidencia.

En segundo lugar, la composición de las poblaciones comparadas es fundamentalmente diferente, dado que la cobertura vacunal en mayores de 40 años superaba el 85-90% para la semana 40 de 2021, el grupo de no vacunados en esas edades representaba un residuo pequeño de la población,posiblemente personas menos expuestas socialmente, con menor movilidad, o con menor acceso al sistema de salud; mientras que el grupo vacunado representaba la inmensa mayoría de la población activa.

En tercer lugar, opera la confusión por indicación, porque las personas fueron vacunadas precisamente porque estaban en mayor riesgo por su edad, de modo que al comparar vacunados contra no vacunados sin ajustar por este factor, creo que se genera una asociación espuria entre la vacunación y el desenlace adverso.

Finalmente, la inmunidad natural no contabilizada puede contribuir a este patrón, porque algunos no vacunados podrían haber tenido una infección previa por SARS-CoV-2 que les confirió protección parcial, reduciendo artificialmente su proporción de incidencia sin que esta condición quedara registrada en los datos de vigilancia rutinaria.

Considero que se requiere para poder interpretar correctamente estas tasas, el ajustar por los confusores, especialmente la edad, antes de emitir juicios sobre la efectividad vacunal.

Este patrón considero que es un ejemplo clásico de la Paradoja de Simpson, que consiste en un fenómeno estadístico en el cual una tendencia que aparece en los datos agregados se invierte o desaparece cuando los datos se estratifican por una tercera variable que actúa como confusor.

2 Estimación de Medidas de Asociación

Para esta sección utilizo los datos de la Tabla 5 del reporte, que presenta la asistencia a servicios de emergencia y las muertes de casos Delta secuenciados y genotipificados en Inglaterra, estratificados por estado vacunal, entre el 1 de febrero y el 12 de septiembre de 2021.

En particular, me concentro en la fila resaltada de la tabla: “Cases where presentation to emergency care resulted in overnight inpatient admission (exclusion)”, que corresponde a los casos de variante Delta que acudieron a urgencias y requirieron admisión hospitalaria con pernocta.

De esta tabla extraigo dos columnas fundamentales para construir la tabla de contingencia: la columna de personas vacunadas con esquema completo (con al menos 14 días transcurridos desde la segunda dosis, identificada como “>=14 days post dose 2”) y la columna de personas no vacunadas (identificada como “Un-vaccinated”).

Las filas que utilizo son la de “Delta cases” (el total de casos por variante Delta) y la de “Cases where presentation to emergency care resulted in overnight inpatient admission (exclusion)” (el desenlace de interés).

2.1 Extracción de datos de la Tabla 5

Comienzo por registrar los cuatro valores directamente observables en la Tabla 5, fila “All cases”:

total_vac   <- 157400
total_novac <- 257357

hosp_vac    <- 2361
hosp_novac  <- 3080

cat("Datos extraidos de la Tabla 5, fila 'All cases':\n\n")
## Datos extraidos de la Tabla 5, fila 'All cases':
cat("Total de casos Delta vacunados (2 dosis):  ", format(total_vac, big.mark = ","), "\n")
## Total de casos Delta vacunados (2 dosis):   157,400
cat("Total de casos Delta no vacunados:          ", format(total_novac, big.mark = ","), "\n")
## Total de casos Delta no vacunados:           257,357
cat("Hospitalizados vacunados (2 dosis):          ", format(hosp_vac, big.mark = ","), "\n")
## Hospitalizados vacunados (2 dosis):           2,361
cat("Hospitalizados no vacunados:                 ", format(hosp_novac, big.mark = ","), "\n")
## Hospitalizados no vacunados:                  3,080

A partir de estos cuatro valores puedo derivar, por sustracción, el número de personas que no fueron hospitalizadas en cada grupo, lo cual es este paso necesario porque la Tabla 5 solo reporta los totales y los hospitalizados, pero no explicita los no hospitalizados:

a <- hosp_vac
b <- total_vac - hosp_vac
c <- hosp_novac
d <- total_novac - hosp_novac

cat("Celdas derivadas por sustraccion:\n\n")
## Celdas derivadas por sustraccion:
cat("No hospitalizados vacunados:     ", format(total_vac, big.mark = ","), "-",
    format(hosp_vac, big.mark = ","), "=", format(b, big.mark = ","), "\n")
## No hospitalizados vacunados:      157,400 - 2,361 = 155,039
cat("No hospitalizados no vacunados:  ", format(total_novac, big.mark = ","), "-",
    format(hosp_novac, big.mark = ","), "=", format(d, big.mark = ","), "\n")
## No hospitalizados no vacunados:   257,357 - 3,080 = 254,277

Con las seis cifras completas explicadas (cuatro celdas interiores y dos totales marginales por fila), presento la tabla de contingencia 2x2 con todos sus totales:

tabla_presentar <- data.frame(
  Estado_vacunal = c("Vacunado (2 dosis)", "No vacunado", "Total"),
  Hospitalizado_Si = c(
    format(a, big.mark = ","),
    format(c, big.mark = ","),
    format(a + c, big.mark = ",")
  ),
  Hospitalizado_No = c(
    format(b, big.mark = ","),
    format(d, big.mark = ","),
    format(b + d, big.mark = ",")
  ),
  Total = c(
    format(a + b, big.mark = ","),
    format(c + d, big.mark = ","),
    format(a + b + c + d, big.mark = ",")
  )
)

kable(tabla_presentar,
      col.names = c("Estado vacunal", "Hospitalizado", "No hospitalizado", "Total"),
      align = c("l", "r", "r", "r")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, position = "center") %>%
  row_spec(3, bold = TRUE, background = "#F0F0F0") %>%
  column_spec(1, bold = TRUE)
Estado vacunal Hospitalizado No hospitalizado Total
Vacunado (2 dosis) 2,361 155,039 157,400
No vacunado 3,080 254,277 257,357
Total 5,441 409,316 414,757
tabla_cruda <- matrix(c(a, b, c, d), nrow = 2, byrow = TRUE,
  dimnames = list(
    "Estado_vacunal"  = c("Vacunado", "No_vacunado"),
    "Hospitalizacion" = c("Si", "No")
  ))

2.2 Pregunta 2a. Odds Ratio crudo

La razon de momios, conocida en ingles como Odds Ratio, compara los momios de hospitalización entre los vacunados contra los momios de hospitalización entre los no vacunados.

Los momios representan la razon entre la probabilidad de que un evento ocurra y la probabilidad de que no ocurra. Para una tabla de contingencia 2x2, la formula es:

\[ OR = \frac{a \times d}{b \times c} \]

numerador_or  <- a * d
denominador_or <- b * c
OR_crudo <- numerador_or / denominador_or

cat("Paso 1: Identificar las celdas de la tabla 2x2\n")
## Paso 1: Identificar las celdas de la tabla 2x2
cat("  a (vacunados hospitalizados)       =", a, "\n")
##   a (vacunados hospitalizados)       = 2361
cat("  b (vacunados no hospitalizados)    =", format(b, big.mark = ","), "\n")
##   b (vacunados no hospitalizados)    = 155,039
cat("  c (no vacunados hospitalizados)    =", c, "\n")
##   c (no vacunados hospitalizados)    = 3080
cat("  d (no vacunados no hospitalizados) =", format(d, big.mark = ","), "\n\n")
##   d (no vacunados no hospitalizados) = 254,277
cat("Paso 2: Calcular el numerador (a x d)\n")
## Paso 2: Calcular el numerador (a x d)
cat(" ", a, "x", format(d, big.mark = ","), "=", format(numerador_or, big.mark = ","), "\n\n")
##   2361 x 254,277 = 600,347,997
cat("Paso 3: Calcular el denominador (b x c)\n")
## Paso 3: Calcular el denominador (b x c)
cat(" ", format(b, big.mark = ","), "x", c, "=", format(denominador_or, big.mark = ","), "\n\n")
##   155,039 x 3080 = 477,520,120
cat("Paso 4: Dividir numerador entre denominador\n")
## Paso 4: Dividir numerador entre denominador
cat(" ", format(numerador_or, big.mark = ","), "/", format(denominador_or, big.mark = ","), "=", OR_crudo, "\n")
##   600,347,997 / 477,520,120 = 1.25722
or_check <- oddsratio.wald(tabla_cruda)

cat("Verificacion con la funcion oddsratio.wald del paquete epitools:\n")
## Verificacion con la funcion oddsratio.wald del paquete epitools:
cat("  OR:", or_check$measure[2, 1], "\n")
##   OR: 1.25722
cat("  IC 95%: [", or_check$measure[2, 2], ",", or_check$measure[2, 3], "]\n")
##   IC 95%: [ 1.191151 , 1.326954 ]

El Odds Ratio crudo resulta en 1.2572 (intervalo de confianza del 95%: 1.1912 a 1.3270). Dado que este intervalo no incluye el valor nulo de 1, la asociacion es estadisticamente significativa.

En terminos de interpretacion, este valor me indica que los momios de hospitalización en vacunados son aproximadamente un 25.7% mayores que en no vacunados.

2.3 Pregunta 2b. Riesgo Relativo crudo

El Riesgo Relativo compara directamente la proporcion de incidencia de hospitalización en cada grupo. A diferencia de la razón de momios, esta medida opera sobre probabilidades, lo cual la hace mas intuitiva, y la formula es la siguiente:

\[ RR = \frac{a \,/\, (a + b)}{c \,/\, (c + d)} \]

Calculo cada componente paso a paso:

riesgo_vac   <- a / (a + b)
riesgo_novac <- c / (c + d)
RR_crudo     <- riesgo_vac / riesgo_novac

cat("Paso 1: Calcular el riesgo en vacunados\n")
## Paso 1: Calcular el riesgo en vacunados
cat("  a / (a + b) =", a, "/", format(a + b, big.mark = ","), "=", riesgo_vac, "\n\n")
##   a / (a + b) = 2361 / 157,400 = 0.015
cat("Paso 2: Calcular el riesgo en no vacunados\n")
## Paso 2: Calcular el riesgo en no vacunados
cat("  c / (c + d) =", c, "/", format(c + d, big.mark = ","), "=", riesgo_novac, "\n\n")
##   c / (c + d) = 3080 / 257,357 = 0.01196781
cat("Paso 3: Dividir riesgo en vacunados entre riesgo en no vacunados\n")
## Paso 3: Dividir riesgo en vacunados entre riesgo en no vacunados
cat(" ", riesgo_vac, "/", riesgo_novac, "=", RR_crudo, "\n")
##   0.015 / 0.01196781 = 1.253362
tabla_epi <- matrix(c(a, b, c, d), nrow = 2, byrow = TRUE)

resultado_epi <- epi.2by2(tabla_epi, method = "cohort.count")

cat("Verificacion con epi.2by2 del paquete epiR:\n")
## Verificacion con epi.2by2 del paquete epiR:
cat("  RR:", resultado_epi$massoc.detail$RR.strata.wald$est, "\n")
##   RR: 1.253362
cat("  IC 95%: [", resultado_epi$massoc.detail$RR.strata.wald$lower, ",",
    resultado_epi$massoc.detail$RR.strata.wald$upper, "]\n")
##   IC 95%: [ 1.188373 , 1.321905 ]

El Riesgo Relativo crudo es de 1.253362, con un intervalo de confianza del 95% de [1.188373, 1.321905] que no incluye el valor nulo de 1. Este resultado sugiere que el riesgo de hospitalización en personas vacunadas es aproximadamente un 25% mayor que en personas no vacunadas.

Debo decir que el Odds Ratio (1.25722) y el Riesgo Relativo (1.253362) son numéricamente muy cercanos, dado que la hospitalización es un evento relativamente infrecuente en ambos grupos (la proporción de hospitalizados no supera el 2%), condición bajo la cual ambas medidas tienden a converger.

2.4 Pregunta 2c. Diferencia de riesgo cruda

La diferencia de riesgo cuantifica en términos absolutos cuánto mayor o menor es el riesgo de hospitalización en vacunados con respecto a no vacunados.

A diferencia de la razón de momios y del Riesgo Relativo, que son medidas relativas, la diferencia de riesgo permite estimar el número adicional de eventos por unidad de población:

\[ DR = \frac{a}{a + b} - \frac{c}{c + d} \]

Calculo cada componente paso a paso:

DR_cruda <- riesgo_vac - riesgo_novac

cat("Paso 1: Riesgo en vacunados (ya calculado)\n")
## Paso 1: Riesgo en vacunados (ya calculado)
cat("  a / (a + b) =", a, "/", format(a + b, big.mark = ","), "=", riesgo_vac, "\n\n")
##   a / (a + b) = 2361 / 157,400 = 0.015
cat("Paso 2: Riesgo en no vacunados (ya calculado)\n")
## Paso 2: Riesgo en no vacunados (ya calculado)
cat("  c / (c + d) =", c, "/", format(c + d, big.mark = ","), "=", riesgo_novac, "\n\n")
##   c / (c + d) = 3080 / 257,357 = 0.01196781
cat("Paso 3: Restar ambos riesgos\n")
## Paso 3: Restar ambos riesgos
cat(" ", riesgo_vac, "-", riesgo_novac, "=", DR_cruda, "\n\n")
##   0.015 - 0.01196781 = 0.003032189
cat("Paso 4: Expresar por 100,000 habitantes\n")
## Paso 4: Expresar por 100,000 habitantes
cat(" ", DR_cruda, "x 100,000 =", DR_cruda * 100000, "por 100,000\n")
##   0.003032189 x 100,000 = 303.2189 por 100,000

La diferencia de riesgo cruda es de 0.003032189, lo cual equivale a 303.2189 hospitalizaciones adicionales por cada 100,000 personas vacunadas en comparación con las no vacunadas. En otras palabras, esta medida cruda sugiere que la vacunación se asocia con un exceso absoluto de hospitalización.

2.5 Pregunta 2d. Interpretación global cruda de la “eficacia” vacunal (VE)

A partir del Riesgo Relativo, calculo la efectividad vacunal cruda como su complemento. La efectividad vacunal representa la proporción de reducción del riesgo en los vacunados con respecto a los no vacunados:

\[ VE = (1 - RR) \times 100 \]

VE_cruda <- (1 - RR_crudo) * 100

cat("Paso 1: Tomo el Riesgo Relativo crudo\n")
## Paso 1: Tomo el Riesgo Relativo crudo
cat("  RR =", RR_crudo, "\n\n")
##   RR = 1.253362
cat("Paso 2: Calculo el complemento\n")
## Paso 2: Calculo el complemento
cat("  1 -", RR_crudo, "=", 1 - RR_crudo, "\n\n")
##   1 - 1.253362 = -0.253362
cat("Paso 3: Expreso como porcentaje\n")
## Paso 3: Expreso como porcentaje
cat("  ", 1 - RR_crudo, "x 100 =", VE_cruda, "%\n")
##    -0.253362 x 100 = -25.3362 %

La efectividad vacunal cruda resulta negativa (VE = -25.33621%), lo cual sugeriría de forma absurda, que la vacuna no solo no protege contra la hospitalización, sino que la incrementa.

¿Estoy de acuerdo con este resultado? No. Considero que este resultado constituye un ejemplo clásico de la Paradoja de Simpson y no debe interpretarse como evidencia de que la vacuna incrementa el riesgo de hospitalización.

Considero que la razón central de esta distorsión es la confusión por edad, que actúa como un confusor clásico en esta relación dado que cumple tres condiciones necesarias:

primero, está asociada con la exposición, dado que la campaña de vacunación del Reino Unido priorizó a los adultos mayores, generando coberturas vacunales superiores al 90% en los grupos de mayor edad.

segundo, está asociada con el desenlace, ya que el riesgo de hospitalización por COVID-19 aumenta de forma exponencial con la edad.

y tercero, NO se encuentra en la vía causal entre vacunación y hospitalización, es decir, la vacuna no causa envejecimiento ni la edad es consecuencia de la vacunación.

En consecuencia, al no ajustar por edad, los datos crudos se muestran mezclados con una alta proporción de adultos mayores vacunados (quienes tienen un riesgo basal elevado de hospitalización independientemente de la vacunación) con una alta proporción de adultos jóvenes no vacunados (quienes tienen un riesgo basal bajo).

Por tanto, esta mezcla desproporcionada invierte artificialmente la dirección de la asociación, por lo cual, considero que resulta indispensable estratificar por edad o aplicar métodos de ajuste como el de Mantel-Haenszel, lo cual desarrollaré en la siguiente sección.

3 Estimación de Medidas de Asociación Ajustadas

Ahora utilizo la misma Tabla 5, pero teniendo en cuenta la edad de los participantes. La tabla presenta los datos estratificados en dos grupos etarios, los menores de 50 años (<50) y de 50 años o más (>=50).

Esta estratificación me permite aplicar el método de Mantel-Haenszel para obtener medidas de asociación ajustadas por edad, y de esta manera evaluar si la dirección de la asociación observada en el análisis crudo se mantiene o se invierte una vez que se controla este confusor.

El primer paso consiste en que voy a construir una tabla de contingencia 2x2 para cada estrato de edad, identificando los valores directamente desde la Tabla 5 del reporte.

3.0.1 Medidas estrato-específicas de Menores de 50 años

total_vac_e1   <- 85407
total_novac_e1 <- 248803
hosp_vac_e1    <- 453
hosp_novac_e1  <- 2416

a1 <- hosp_vac_e1
b1 <- total_vac_e1 - hosp_vac_e1
c1 <- hosp_novac_e1
d1 <- total_novac_e1 - hosp_novac_e1

3.1 Estrato 1: Menores de 50 años

tabla_e1_color <- data.frame(
  Estado = c("Vacunado (2 dosis)", "No vacunado", "Total"),
  Hospitalizado = c(
    paste0("a1 = ", format(a1, big.mark = ",")),
    paste0("c1 = ", format(c1, big.mark = ",")),
    format(a1 + c1, big.mark = ",")
  ),
  No_hospitalizado = c(
    paste0("b1 = ", format(b1, big.mark = ",")),
    paste0("d1 = ", format(d1, big.mark = ",")),
    format(b1 + d1, big.mark = ",")
  ),
  Total = c(
    paste0("a1+b1 = ", format(a1 + b1, big.mark = ",")),
    paste0("c1+d1 = ", format(c1 + d1, big.mark = ",")),
    format(a1 + b1 + c1 + d1, big.mark = ",")
  )
)

kable(tabla_e1_color,
      col.names = c("Estado vacunal", "Hospitalizado", "No hospitalizado", "Total"),
      align = c("l", "r", "r", "r")) %>%
  kable_styling(bootstrap_options = c("bordered", "condensed"),
                full_width = FALSE, position = "center") %>%
  row_spec(1, background = "#D6EAF8") %>%
  row_spec(2, background = "#FADBD8") %>%
  row_spec(3, bold = TRUE, background = "#E8E8E8") %>%
  column_spec(1, bold = TRUE) %>%
  column_spec(2, color = "#1A5276") %>%
  column_spec(3, color = "#1A5276")
Estado vacunal Hospitalizado No hospitalizado Total
Vacunado (2 dosis) a1 = 453 b1 = 84,954 a1+b1 = 85,407
No vacunado c1 = 2,416 d1 = 246,387 c1+d1 = 248,803
Total 2,869 331,341 334,210
T1 <- a1 + b1 + c1 + d1

cat("=== ODDS RATIO - Estrato menores de 50 anios ===\n\n")
## === ODDS RATIO - Estrato menores de 50 anios ===
cat("Paso 1: Numerador (a1 x d1)\n")
## Paso 1: Numerador (a1 x d1)
cat(" ", a1, "x", format(d1, big.mark = ","), "=", format(a1 * d1, big.mark = ","), "\n\n")
##   453 x 246,387 = 111,613,311
cat("Paso 2: Denominador (b1 x c1)\n")
## Paso 2: Denominador (b1 x c1)
cat(" ", format(b1, big.mark = ","), "x", format(c1, big.mark = ","), "=", format(b1 * c1, big.mark = ","), "\n\n")
##   84,954 x 2,416 = 205,248,864
OR_e1 <- (a1 * d1) / (b1 * c1)
cat("Paso 3: Dividir\n")
## Paso 3: Dividir
cat("  OR =", format(a1 * d1, big.mark = ","), "/", format(b1 * c1, big.mark = ","), "=", OR_e1, "\n")
##   OR = 111,613,311 / 205,248,864 = 0.543795
cat("=== RIESGO RELATIVO - Estrato menores de 50 anios ===\n\n")
## === RIESGO RELATIVO - Estrato menores de 50 anios ===
riesgo_vac_e1 <- a1 / (a1 + b1)
cat("Paso 1: Riesgo en vacunados\n")
## Paso 1: Riesgo en vacunados
cat("  a1 / (a1+b1) =", a1, "/", format(a1 + b1, big.mark = ","), "=", riesgo_vac_e1, "\n\n")
##   a1 / (a1+b1) = 453 / 85,407 = 0.005304015
riesgo_novac_e1 <- c1 / (c1 + d1)
cat("Paso 2: Riesgo en no vacunados\n")
## Paso 2: Riesgo en no vacunados
cat("  c1 / (c1+d1) =", format(c1, big.mark = ","), "/", format(c1 + d1, big.mark = ","), "=", riesgo_novac_e1, "\n\n")
##   c1 / (c1+d1) = 2,416 / 248,803 = 0.009710494
RR_e1 <- riesgo_vac_e1 / riesgo_novac_e1
cat("Paso 3: Dividir\n")
## Paso 3: Dividir
cat("  RR =", riesgo_vac_e1, "/", riesgo_novac_e1, "=", RR_e1, "\n")
##   RR = 0.005304015 / 0.009710494 = 0.5462147
cat("=== DIFERENCIA DE RIESGO - Estrato menores de 50 anios ===\n\n")
## === DIFERENCIA DE RIESGO - Estrato menores de 50 anios ===
DR_e1 <- riesgo_vac_e1 - riesgo_novac_e1
cat("Paso 1: Restar riesgos\n")
## Paso 1: Restar riesgos
cat("  DR =", riesgo_vac_e1, "-", riesgo_novac_e1, "=", DR_e1, "\n\n")
##   DR = 0.005304015 - 0.009710494 = -0.004406479
cat("Paso 2: Expresar por 100,000\n")
## Paso 2: Expresar por 100,000
cat("  DR =", DR_e1 * 100000, "por 100,000\n")
##   DR = -440.6479 por 100,000
resumen_e1 <- data.frame(
  Medida = c("Razon de momios (Odds Ratio)",
             "Riesgo Relativo",
             "Diferencia de riesgo",
             "Diferencia de riesgo por 100,000"),
  Valor = as.character(c(OR_e1, RR_e1, DR_e1, DR_e1 * 100000)),
  Direccion = c("Menor que 1: proteccion",
                "Menor que 1: proteccion",
                "Negativa: proteccion",
                "Negativa: proteccion")
)

kable(resumen_e1,
      col.names = c("Medida de asociacion", "Valor", "Direccion"),
      align = c("l", "r", "l")) %>%
  kable_styling(bootstrap_options = c("bordered", "condensed"),
                full_width = FALSE, position = "center") %>%
  row_spec(1:4, background = "#D6EAF8") %>%
  column_spec(1, bold = TRUE) %>%
  column_spec(2, bold = TRUE, color = "#1A5276")
Medida de asociacion Valor Direccion
Razon de momios (Odds Ratio) 0.543795024366128 Menor que 1: proteccion
Riesgo Relativo 0.546214742351329 Menor que 1: proteccion
Diferencia de riesgo -0.00440647895113478 Negativa: proteccion
Diferencia de riesgo por 100,000 -440.647895113478 Negativa: proteccion

En el estrato de menores de 50 años, las tres medidas indican de forma consistente que la vacunacion es protectora, pues la razon de momios y el Riesgo Relativo son menores que 1, y la diferencia de riesgo es negativa (-440.65 hospitalizaciones menos por cada 100,000 vacunados), además este resultado contrasta radicalmente con las medidas crudas de la seccion anterior.

3.2 Estrato 2: 50 años o mas

total_vac_e2   <- 71991
total_novac_e2 <- 8551
hosp_vac_e2    <- 1908
hosp_novac_e2  <- 664

a2 <- hosp_vac_e2
b2 <- total_vac_e2 - hosp_vac_e2
c2 <- hosp_novac_e2
d2 <- total_novac_e2 - hosp_novac_e2

cat("Datos de la Tabla 5, fila '>=50', variante Delta:\n\n")
## Datos de la Tabla 5, fila '>=50', variante Delta:
cat("Total de casos vacunados (2 dosis):  ", format(total_vac_e2, big.mark = ","), "\n")
## Total de casos vacunados (2 dosis):   71,991
cat("Total de casos no vacunados:          ", format(total_novac_e2, big.mark = ","), "\n")
## Total de casos no vacunados:           8,551
cat("Hospitalizados vacunados:              ", format(hosp_vac_e2, big.mark = ","), "\n")
## Hospitalizados vacunados:               1,908
cat("Hospitalizados no vacunados:           ", hosp_novac_e2, "\n\n")
## Hospitalizados no vacunados:            664
cat("Celdas derivadas por sustraccion:\n")
## Celdas derivadas por sustraccion:
cat("No hospitalizados vacunados:     ", format(total_vac_e2, big.mark = ","), "-",
    format(hosp_vac_e2, big.mark = ","), "=", format(b2, big.mark = ","), "\n")
## No hospitalizados vacunados:      71,991 - 1,908 = 70,083
cat("No hospitalizados no vacunados:  ", format(total_novac_e2, big.mark = ","), "-",
    hosp_novac_e2, "=", format(d2, big.mark = ","), "\n")
## No hospitalizados no vacunados:   8,551 - 664 = 7,887
tabla_e2_color <- data.frame(
  Estado = c("Vacunado (2 dosis)", "No vacunado", "Total"),
  Hospitalizado = c(
    paste0("a2 = ", format(a2, big.mark = ",")),
    paste0("c2 = ", format(c2, big.mark = ",")),
    format(a2 + c2, big.mark = ",")
  ),
  No_hospitalizado = c(
    paste0("b2 = ", format(b2, big.mark = ",")),
    paste0("d2 = ", format(d2, big.mark = ",")),
    format(b2 + d2, big.mark = ",")
  ),
  Total = c(
    paste0("a2+b2 = ", format(a2 + b2, big.mark = ",")),
    paste0("c2+d2 = ", format(c2 + d2, big.mark = ",")),
    format(a2 + b2 + c2 + d2, big.mark = ",")
  )
)

kable(tabla_e2_color,
      col.names = c("Estado vacunal", "Hospitalizado", "No hospitalizado", "Total"),
      align = c("l", "r", "r", "r")) %>%
  kable_styling(bootstrap_options = c("bordered", "condensed"),
                full_width = FALSE, position = "center") %>%
  row_spec(1, background = "#D5F5E3") %>%
  row_spec(2, background = "#FCF3CF") %>%
  row_spec(3, bold = TRUE, background = "#E8E8E8") %>%
  column_spec(1, bold = TRUE) %>%
  column_spec(2, color = "#1A5276") %>%
  column_spec(3, color = "#1A5276")
Estado vacunal Hospitalizado No hospitalizado Total
Vacunado (2 dosis) a2 = 1,908 b2 = 70,083 a2+b2 = 71,991
No vacunado c2 = 664 d2 = 7,887 c2+d2 = 8,551
Total 2,572 77,970 80,542

3.2.1 Medidas estrato-especificas: 50 años o mas

T2 <- a2 + b2 + c2 + d2

cat("=== ODDS RATIO - Estrato 50 anios o mas ===\n\n")
## === ODDS RATIO - Estrato 50 anios o mas ===
cat("Paso 1: Numerador (a2 x d2)\n")
## Paso 1: Numerador (a2 x d2)
cat(" ", format(a2, big.mark = ","), "x", format(d2, big.mark = ","), "=", format(a2 * d2, big.mark = ","), "\n\n")
##   1,908 x 7,887 = 15,048,396
cat("Paso 2: Denominador (b2 x c2)\n")
## Paso 2: Denominador (b2 x c2)
cat(" ", format(b2, big.mark = ","), "x", c2, "=", format(b2 * c2, big.mark = ","), "\n\n")
##   70,083 x 664 = 46,535,112
OR_e2 <- (a2 * d2) / (b2 * c2)
cat("Paso 3: Dividir\n")
## Paso 3: Dividir
cat("  OR =", format(a2 * d2, big.mark = ","), "/", format(b2 * c2, big.mark = ","), "=", OR_e2, "\n")
##   OR = 15,048,396 / 46,535,112 = 0.3233772
cat("=== RIESGO RELATIVO - Estrato 50 anios o mas ===\n\n")
## === RIESGO RELATIVO - Estrato 50 anios o mas ===
riesgo_vac_e2 <- a2 / (a2 + b2)
cat("Paso 1: Riesgo en vacunados\n")
## Paso 1: Riesgo en vacunados
cat("  a2 / (a2+b2) =", format(a2, big.mark = ","), "/", format(a2 + b2, big.mark = ","), "=", riesgo_vac_e2, "\n\n")
##   a2 / (a2+b2) = 1,908 / 71,991 = 0.02650331
riesgo_novac_e2 <- c2 / (c2 + d2)
cat("Paso 2: Riesgo en no vacunados\n")
## Paso 2: Riesgo en no vacunados
cat("  c2 / (c2+d2) =", c2, "/", format(c2 + d2, big.mark = ","), "=", riesgo_novac_e2, "\n\n")
##   c2 / (c2+d2) = 664 / 8,551 = 0.07765174
RR_e2 <- riesgo_vac_e2 / riesgo_novac_e2
cat("Paso 3: Dividir\n")
## Paso 3: Dividir
cat("  RR =", riesgo_vac_e2, "/", riesgo_novac_e2, "=", RR_e2, "\n")
##   RR = 0.02650331 / 0.07765174 = 0.34131
cat("=== DIFERENCIA DE RIESGO - Estrato 50 o mas ===\n\n")
## === DIFERENCIA DE RIESGO - Estrato 50 o mas ===
DR_e2 <- riesgo_vac_e2 - riesgo_novac_e2

cat("Paso 1: Riesgo en vacunados\n")
## Paso 1: Riesgo en vacunados
cat("  ", a2, "/", format(a2 + b2, big.mark = ","), "=", riesgo_vac_e2, "\n\n")
##    1908 / 71,991 = 0.02650331
cat("Paso 2: Riesgo en no vacunados\n")
## Paso 2: Riesgo en no vacunados
cat("  ", c2, "/", format(c2 + d2, big.mark = ","), "=", riesgo_novac_e2, "\n\n")
##    664 / 8,551 = 0.07765174
cat("Paso 3: Restar (en escala original, sin multiplicar)\n")
## Paso 3: Restar (en escala original, sin multiplicar)
cat("  DR =", riesgo_vac_e2, "-", riesgo_novac_e2, "=", DR_e2, "\n\n")
##   DR = 0.02650331 - 0.07765174 = -0.05114842
cat("Paso 4: Expresar por 100,000\n")
## Paso 4: Expresar por 100,000
cat("  DR =", DR_e2, "x 100,000 =", format(DR_e2 * 100000, big.mark = ","), "por 100,000\n")
##   DR = -0.05114842 x 100,000 = -5,114.842 por 100,000

En el estrato de 50 años o más, las tres medidas confirman nuevamente que la vacunación es protectora,porque la razón de momios es 0.3233772 (menor que 1), el Riesgo Relativo es 0.341309983024987 (menor que 1), y la diferencia de riesgo es -0.05114842 (negativa), equivalente a -5,114.842 hospitalizaciones menos por cada 100,000 personas vacunadas en comparación con las no vacunadas.

Este último valor resulta considerablemente mayor en magnitud absoluta que el del estrato de menores de 50 años (donde la diferencia de riesgo fue de -0.00440647, equivalente a -440.6479 por 100,000).

Esta diferencia se explica porque los adultos mayores tienen un riesgo basal de hospitalización mucho más elevado, pues en el estrato de 50 o más años, el riesgo en no vacunados es de 0.077651736638990, frente a un riesgo de 0.009710494 en el estrato de menores de 50 años, de tal manera que el beneficio absoluto de la vacunación es proporcionalmente mayor en el grupo etario de mayor riesgo.

De esta manera, al estratificar por edad, ambos estratos muestran consistentemente que la vacunación reduce el riesgo de hospitalización, lo cual constituye la evidencia directa de que la Paradoja de Simpson estaba operando en el análisis crudo, donde las medidas de asociación me sugerían falsamente lo contrario.

3.2.2 Estimación combinada por el método de Mantel-Haenszel

Una vez que he calculado las medidas estrato-específicas para ambos grupos de edad, procedo a obtener las medidas de asociación ajustadas por edad mediante el método de Mantel-Haenszel.

Este método pondera la contribución de cada estrato según su tamaño, lo cual me permite obtener una estimación global que controla la confusión por la variable de estratificación (en este caso, la edad).

A continuación presento las fórmulas completas y el cálculo paso a paso para cada medida.

# ============================================================================
# ** CÁLCULOS DEL ODDS RATIO DE MANTEL-HAENSZEL **
# Estos valores alimentan las tablas que presentaré a continuación.
# ============================================================================

T1 <- a1 + b1 + c1 + d1
T2 <- a2 + b2 + c2 + d2

# ** Componentes del numerador **
num_or_e1 <- (a1 * d1) / T1
num_or_e2 <- (a2 * d2) / T2
numerador_OR_MH <- num_or_e1 + num_or_e2

# ** Componentes del denominador **
den_or_e1 <- (b1 * c1) / T1
den_or_e2 <- (b2 * c2) / T2
denominador_OR_MH <- den_or_e1 + den_or_e2

# ** Resultado **
OR_MH <- numerador_OR_MH / denominador_OR_MH
# ============================================================================
# ** TABLA 1: VALORES DE CADA ESTRATO PARA EL ODDS RATIO DE MANTEL-HAENSZEL **
# ============================================================================

tabla_datos_or <- data.frame(
  Estrato = c("Menores de 50 anios", "50 anios o mas"),
  Vacunados_hosp = c(format(a1, big.mark = ","), format(a2, big.mark = ",")),
  Vacunados_no_hosp = c(format(b1, big.mark = ","), format(b2, big.mark = ",")),
  NoVac_hosp = c(format(c1, big.mark = ","), format(c2, big.mark = ",")),
  NoVac_no_hosp = c(format(d1, big.mark = ","), format(d2, big.mark = ",")),
  Total_estrato = c(format(T1, big.mark = ","), format(T2, big.mark = ","))
)

kable(tabla_datos_or,
      col.names = c("Estrato de edad",
                     "Vacunados hospitalizados",
                     "Vacunados no hospitalizados",
                     "No vacunados hospitalizados",
                     "No vacunados no hospitalizados",
                     "Total del estrato"),
      align = c("l", "r", "r", "r", "r", "r"),
      caption = "1° Tabla de Valores de cada estrato extraidos de las tablas de contingencia 2x2") %>%
  kable_styling(
    bootstrap_options = c("bordered", "hover", "condensed"),
    full_width = FALSE, position = "center"
  ) %>%
  row_spec(0, background = "#B3E5FC", color = "#2C3E50", bold = TRUE) %>%
  row_spec(1, background = "#E2F0CB") %>%
  row_spec(2, background = "#FFF1B6") %>%
  column_spec(1, bold = TRUE, width = "10em", color = "#2C3E50") %>%
  column_spec(2, background = "#CDE7FF", color = "#2C3E50") %>%
  column_spec(3, background = "#DCEBFF", color = "#2C3E50") %>%
  column_spec(4, background = "#FFD8BE", color = "#2C3E50") %>%
  column_spec(5, background = "#FFCCBC", color = "#2C3E50") %>%
  column_spec(6, bold = TRUE, background = "#E4C1F9", color = "#2C3E50")
1° Tabla de Valores de cada estrato extraidos de las tablas de contingencia 2x2
Estrato de edad Vacunados hospitalizados Vacunados no hospitalizados No vacunados hospitalizados No vacunados no hospitalizados Total del estrato
Menores de 50 anios 453 84,954 2,416 246,387 334,210
50 anios o mas 1,908 70,083 664 7,887 80,542
# ============================================================================
# ** TABLA 2: CALCULO DEL NUMERADOR DEL ODDS RATIO DE MANTEL-HAENSZEL **
# ============================================================================

tabla_numerador_or <- data.frame(
  Estrato = c("Menores de 50 anios", "50 anios o mas", "SUMA (numerador)"),
  Vac_hosp = c(format(a1, big.mark = ","), format(a2, big.mark = ","), ""),
  Simbolo1 = c("x", "x", ""),
  NoVac_no_hosp = c(format(d1, big.mark = ","), format(d2, big.mark = ","), ""),
  Simbolo2 = c("/", "/", ""),
  Total = c(format(T1, big.mark = ","), format(T2, big.mark = ","), ""),
  Simbolo3 = c("=", "=", "="),
  Resultado = c(num_or_e1, num_or_e2, numerador_OR_MH)
)

kable(tabla_numerador_or,
      col.names = c("Estrato de edad",
                     "Vacunados hospitalizados",
                     "",
                     "No vacunados no hospitalizados",
                     "",
                     "Total del estrato",
                     "",
                     "Resultado"),
      align = c("l", "r", "c", "r", "c", "r", "c", "r"),
      caption = "2° Tabla Numerador: vacunados hospitalizados x no vacunados no hospitalizados / Total del estrato") %>%
  kable_styling(
    bootstrap_options = c("bordered", "striped", "hover", "condensed"),
    full_width = FALSE, position = "center"
  ) %>%
  row_spec(0, background = "#B9F6CA", color = "#2C3E50", bold = TRUE) %>%
  row_spec(1, background = "#E2F0CB") %>%
  row_spec(2, background = "#EAFAF1") %>%
  row_spec(3, background = "#B9F6CA", bold = TRUE, color = "#1B5E20") %>%
  column_spec(1, bold = TRUE, width = "10em", color = "#2C3E50") %>%
  column_spec(2, background = "#CDE7FF", color = "#2C3E50") %>%
  column_spec(3, width = "1.5em", bold = TRUE) %>%
  column_spec(4, background = "#FFCCBC", color = "#2C3E50") %>%
  column_spec(5, width = "1.5em", bold = TRUE) %>%
  column_spec(6, background = "#E4C1F9", color = "#2C3E50") %>%
  column_spec(7, width = "1.5em", bold = TRUE) %>%
  column_spec(8, bold = TRUE, color = "#1B5E20", background = "#E2F0CB")
2° Tabla Numerador: vacunados hospitalizados x no vacunados no hospitalizados / Total del estrato
Estrato de edad Vacunados hospitalizados No vacunados no hospitalizados Total del estrato Resultado
Menores de 50 anios 453 x 246,387 / 334,210 = 333.9616
50 anios o mas 1,908 x 7,887 / 80,542 = 186.8391
SUMA (numerador) = 520.8007
# ============================================================================
# ** TABLA 3: CALCULO DEL DENOMINADOR DEL ODDS RATIO DE MANTEL-HAENSZEL **
# ============================================================================

tabla_denominador_or <- data.frame(
  Estrato = c("Menores de 50 anios", "50 anios o mas", "SUMA (denominador)"),
  Vac_no_hosp = c(format(b1, big.mark = ","), format(b2, big.mark = ","), ""),
  Simbolo1 = c("x", "x", ""),
  NoVac_hosp = c(format(c1, big.mark = ","), format(c2, big.mark = ","), ""),
  Simbolo2 = c("/", "/", ""),
  Total = c(format(T1, big.mark = ","), format(T2, big.mark = ","), ""),
  Simbolo3 = c("=", "=", "="),
  Resultado = c(den_or_e1, den_or_e2, denominador_OR_MH)
)

kable(tabla_denominador_or,
      col.names = c("Estrato de edad",
                     "Vacunados no hospitalizados",
                     "",
                     "No vacunados hospitalizados",
                     "",
                     "Total del estrato",
                     "",
                     "Resultado"),
      align = c("l", "r", "c", "r", "c", "r", "c", "r"),
      caption = "3° Tabla Denominador: vacunados no hospitalizados x no vacunados hospitalizados / Total del estrato") %>%
  kable_styling(
    bootstrap_options = c("bordered", "striped", "hover", "condensed"),
    full_width = FALSE, position = "center"
  ) %>%
  row_spec(0, background = "#FFCCBC", color = "#2C3E50", bold = TRUE) %>%
  row_spec(1, background = "#FFF1B6") %>%
  row_spec(2, background = "#FFE0B2") %>%
  row_spec(3, background = "#FFCCBC", bold = TRUE, color = "#BF360C") %>%
  column_spec(1, bold = TRUE, width = "10em", color = "#2C3E50") %>%
  column_spec(2, background = "#DCEBFF", color = "#2C3E50") %>%
  column_spec(3, width = "1.5em", bold = TRUE) %>%
  column_spec(4, background = "#FFD8BE", color = "#2C3E50") %>%
  column_spec(5, width = "1.5em", bold = TRUE) %>%
  column_spec(6, background = "#E4C1F9", color = "#2C3E50") %>%
  column_spec(7, width = "1.5em", bold = TRUE) %>%
  column_spec(8, bold = TRUE, color = "#BF360C", background = "#FFF1B6")
3° Tabla Denominador: vacunados no hospitalizados x no vacunados hospitalizados / Total del estrato
Estrato de edad Vacunados no hospitalizados No vacunados hospitalizados Total del estrato Resultado
Menores de 50 anios 84,954 x 2,416 / 334,210 = 614.1314
50 anios o mas 70,083 x 664 / 80,542 = 577.7745
SUMA (denominador) = 1191.9059
# ============================================================================
# ** TABLA 4: RESULTADO FINAL DEL ODDS RATIO DE MANTEL-HAENSZEL **
# ============================================================================

tabla_final_or <- data.frame(
  Elemento = c(
    "Numerador (suma de ambos estratos)",
    "Denominador (suma de ambos estratos)",
    "Odds Ratio de Mantel-Haenszel = numerador / denominador",
    "Reduccion porcentual de los momios = (1 - Odds Ratio) x 100"
  ),
  Valor = c(
    numerador_OR_MH,
    denominador_OR_MH,
    OR_MH,
    (1 - OR_MH) * 100
  )
)

kable(tabla_final_or,
      col.names = c("Elemento del calculo", "Valor obtenido"),
      align = c("l", "r"),
      caption = "4°tabla: Resultado del Odds Ratio ajustado por edad (metodo de Mantel-Haenszel)") %>%
  kable_styling(
    bootstrap_options = c("bordered", "striped", "hover", "condensed"),
    full_width = FALSE, position = "center"
  ) %>%
  row_spec(0, background = "#B3E5FC", color = "#2C3E50", bold = TRUE) %>%
  row_spec(1, background = "#E2F0CB", color = "#1B5E20") %>%
  row_spec(2, background = "#FFCCBC", color = "#BF360C") %>%
  row_spec(3, background = "#B3E5FC", bold = TRUE, color = "#1A5276") %>%
  row_spec(4, background = "#FFF1B6", bold = TRUE, color = "#7D6608") %>%
  column_spec(1, bold = TRUE, width = "35em", color = "#2C3E50") %>%
  column_spec(2, bold = TRUE, width = "12em")
4°tabla: Resultado del Odds Ratio ajustado por edad (metodo de Mantel-Haenszel)
Elemento del calculo Valor obtenido
Numerador (suma de ambos estratos) 520.8007290
Denominador (suma de ambos estratos) 1191.9059043
Odds Ratio de Mantel-Haenszel = numerador / denominador 0.4369479
Reduccion porcentual de los momios = (1 - Odds Ratio) x 100 56.3052144

3.2.3 Pregunta 3a. Odds Ratio ajustado por edad (método de Mantel-Haenszel)

Para controlar la confusión por edad que identifiqué en el análisis crudo, apliqué el método de Mantel-Haenszel, pues este método toma las tablas de contingencia 2×2 de cada estrato de edad y las combina en una sola estimación, ponderando la contribución de cada estrato según su tamaño, por tanto la fórmula que empleé es la siguiente:

\[OR_{MH} = \frac{\text{Suma de (vacunados hospitalizados} \times \text{no vacunados no hospitalizados / Total del estrato)}}{\text{Suma de (vacunados no hospitalizados} \times \text{no vacunados hospitalizados / Total del estrato)}}\]

En la primera tabla presenté los valores que alimentan esta fórmula, extraídos directamente de las tablas de contingencia que construí para cada grupo de edad.

En la segunda tabla calculé el numerador, que refleja la proporción ponderada en que la hospitalización se concentra entre los vacunados.

En la tercera tabla calculé el denominador, que refleja la proporción ponderada en que la hospitalización se concentra entre los no vacunados.

El hecho de que el denominador (1,191.9059) sea sustancialmente mayor que el numerador (520.8007) ya me anticipaba el resultado, pues dentro de cada estrato de edad, la hospitalización se concentra proporcionalmente más entre los no vacunados.

En la última tabla, el Odds Ratio ajustado que obtuve es de 0.4369479, lo cual significa que los momios de hospitalización en personas vacunadas con esquema completo son aproximadamente un 56.3% menores que en personas no vacunadas, cuando se comparan personas del mismo grupo etario.

Considero que la magnitud del cambio entre la estimación cruda y la ajustada es lo más revelador de este ejercicio, pues el Odds Ratio crudo fue de 1.25722, lo cual sugería falsamente que la vacunación incrementaba los momios de hospitalización, pero al ajustar por edad, la dirección se inviertió por completo, pues pasó de mayor que 1 (aparente riesgo) a menor que 1 (protección clara).

Esta inversión considero que es la expresión cuantitativa de la Paradoja de Simpson y me confirma que la edad actuaba como un confusor clásico, dado que cumple las tres condiciones, primero que se asocia con la exposición (las personas mayores fueron priorizadas para vacunación),segundo que se asocia con el desenlace (el riesgo de hospitalización aumenta con la edad), y tercero que no se encuentra en la vía causal entre vacunación y hospitalización.

3.2.4 Pregunta 3b. Riesgo Relativo ajustado por edad (metodo de Mantel-Haenszel)

De la misma manera que procedí con el Odds Ratio, ahora aplicaré el método de Mantel-Haenszel para obtener el Riesgo Relativo ajustado por edad, pues a diferencia del Odds Ratio, que compara momios, el Riesgo Relativo compara directamente las proporciones de incidencia de hospitalización entre vacunados y no vacunados, lo cual según la literatura, lo hace una medida más intuitiva en estudios de cohorte.

La fórmula del Riesgo Relativo de Mantel-Haenszel es:

\[RR_{MH} = \frac{\text{Suma de (vacunados hospitalizados} \times \text{total de no vacunados / Total del estrato)}}{\text{Suma de (no vacunados hospitalizados} \times \text{total de vacunados / Total del estrato)}}\]

La estructura es diferente a la del Odds Ratio, porque en el numerador multiplico los vacunados hospitalizados por el total de no vacunados del estrato (no solo los no hospitalizados), y en el denominador multiplico los no vacunados hospitalizados por el total de vacunados del estrato, por lo cual, en esta construcción estoy comparando riesgos (proporciones) y no momios (razones de probabilidades).

A continuación presento el proceso:

# ============================================================================
# ** TABLA 1. VALORES DE CADA ESTRATO PARA EL RIESGO RELATIVO DE MH **
# ============================================================================

tabla_datos_rr <- data.frame(
  Estrato = c("Menores de 50 anios", "50 anios o mas"),
  Vac_hosp = c(format(a1, big.mark = ","), format(a2, big.mark = ",")),
  Total_novac = c(format(c1 + d1, big.mark = ","), format(c2 + d2, big.mark = ",")),
  NoVac_hosp = c(format(c1, big.mark = ","), format(c2, big.mark = ",")),
  Total_vac = c(format(a1 + b1, big.mark = ","), format(a2 + b2, big.mark = ",")),
  Total_estrato = c(format(T1, big.mark = ","), format(T2, big.mark = ","))
)

kable(tabla_datos_rr,
      col.names = c("Estrato de edad",
                     "Vacunados hospitalizados",
                     "Total de no vacunados",
                     "No vacunados hospitalizados",
                     "Total de vacunados",
                     "Total del estrato"),
      align = c("l", "r", "r", "r", "r", "r"),
      caption = "Valores de cada estrato para el calculo del Riesgo Relativo de Mantel-Haenszel") %>%
  kable_styling(
    bootstrap_options = c("bordered", "hover", "condensed"),
    full_width = FALSE, position = "center"
  ) %>%
  row_spec(0, background = "#B3E5FC", color = "#2C3E50", bold = TRUE) %>%
  row_spec(1, background = "#E2F0CB") %>%
  row_spec(2, background = "#FFF1B6") %>%
  column_spec(1, bold = TRUE, width = "10em", color = "#2C3E50") %>%
  column_spec(2, background = "#CDE7FF", color = "#2C3E50") %>%
  column_spec(3, background = "#DCEBFF", color = "#2C3E50") %>%
  column_spec(4, background = "#FFD8BE", color = "#2C3E50") %>%
  column_spec(5, background = "#FFCCBC", color = "#2C3E50") %>%
  column_spec(6, bold = TRUE, background = "#E4C1F9", color = "#2C3E50")
Valores de cada estrato para el calculo del Riesgo Relativo de Mantel-Haenszel
Estrato de edad Vacunados hospitalizados Total de no vacunados No vacunados hospitalizados Total de vacunados Total del estrato
Menores de 50 anios 453 248,803 2,416 85,407 334,210
50 anios o mas 1,908 8,551 664 71,991 80,542

Observo que en el estrato de menores de 50 anios predominan los no vacunados (248,803 frente a 85,407 vacunados), mientras que en el estrato de 50 anios o mas la relacion se invierte drasticamente (71,991 vacunados frente a apenas 8,551 no vacunados).

Esta asimetria es precisamente lo que considero que es la fuente de la confusion que distorsionaba el análisis crudo, y es lo que el método de Mantel-Haenszel corrige al ponderar cada estrato por separado.

# ============================================================================
# ** TABLA 2: CALCULO DEL NUMERADOR DEL RIESGO RELATIVO DE MANTEL-HAENSZEL **
# ============================================================================

tabla_numerador_rr <- data.frame(
  Estrato = c("Menores de 50 anios", "50 anios o mas", "SUMA (numerador)"),
  Vac_hosp = c(format(a1, big.mark = ","), format(a2, big.mark = ","), ""),
  Simbolo1 = c("x", "x", ""),
  Total_novac = c(format(c1 + d1, big.mark = ","), format(c2 + d2, big.mark = ","), ""),
  Simbolo2 = c("/", "/", ""),
  Total = c(format(T1, big.mark = ","), format(T2, big.mark = ","), ""),
  Simbolo3 = c("=", "=", "="),
  Resultado = c(num_rr_e1, num_rr_e2, numerador_RR_MH)
)

kable(tabla_numerador_rr,
      col.names = c("Estrato de edad",
                     "Vacunados hospitalizados",
                     "",
                     "Total de no vacunados",
                     "",
                     "Total del estrato",
                     "",
                     "Resultado"),
      align = c("l", "r", "c", "r", "c", "r", "c", "r"),
      caption = "Numerador: vacunados hospitalizados x total de no vacunados / Total del estrato") %>%
  kable_styling(
    bootstrap_options = c("bordered", "striped", "hover", "condensed"),
    full_width = FALSE, position = "center"
  ) %>%
  row_spec(0, background = "#B9F6CA", color = "#2C3E50", bold = TRUE) %>%
  row_spec(1, background = "#E2F0CB") %>%
  row_spec(2, background = "#EAFAF1") %>%
  row_spec(3, background = "#B9F6CA", bold = TRUE, color = "#1B5E20") %>%
  column_spec(1, bold = TRUE, width = "10em", color = "#2C3E50") %>%
  column_spec(2, background = "#CDE7FF", color = "#2C3E50") %>%
  column_spec(3, width = "1.5em", bold = TRUE) %>%
  column_spec(4, background = "#FFCCBC", color = "#2C3E50") %>%
  column_spec(5, width = "1.5em", bold = TRUE) %>%
  column_spec(6, background = "#E4C1F9", color = "#2C3E50") %>%
  column_spec(7, width = "1.5em", bold = TRUE) %>%
  column_spec(8, bold = TRUE, color = "#1B5E20", background = "#E2F0CB")
Numerador: vacunados hospitalizados x total de no vacunados / Total del estrato
Estrato de edad Vacunados hospitalizados Total de no vacunados Total del estrato Resultado
Menores de 50 anios 453 x 248,803 / 334,210 = 337.2363
50 anios o mas 1,908 x 8,551 / 80,542 = 202.5689
SUMA (numerador) = 539.8053

El numerador pondera, para cada estrato, los vacunados hospitalizados por el tamaño del grupo no vacunado, por ejemplo, el estrato de menores de 50 años aporta 337.2363 y el de 50 años o mas aporta 202.5689, sumando un numerador total de 539.8053.

Observo que el estrato joven contribuye mas al numerador porque tiene un grupo de no vacunados mucho mayor (248,803 frente a 8,551), lo cual le da mayor peso en la ponderación.

# ============================================================================
# ** TABLA 3: CALCULO DEL DENOMINADOR DEL RIESGO RELATIVO DE MANTEL-HAENSZEL **
# ============================================================================

tabla_denominador_rr <- data.frame(
  Estrato = c("Menores de 50 anios", "50 anios o mas", "SUMA (denominador)"),
  NoVac_hosp = c(format(c1, big.mark = ","), format(c2, big.mark = ","), ""),
  Simbolo1 = c("x", "x", ""),
  Total_vac = c(format(a1 + b1, big.mark = ","), format(a2 + b2, big.mark = ","), ""),
  Simbolo2 = c("/", "/", ""),
  Total = c(format(T1, big.mark = ","), format(T2, big.mark = ","), ""),
  Simbolo3 = c("=", "=", "="),
  Resultado = c(den_rr_e1, den_rr_e2, denominador_RR_MH)
)

kable(tabla_denominador_rr,
      col.names = c("Estrato de edad",
                     "No vacunados hospitalizados",
                     "",
                     "Total de vacunados",
                     "",
                     "Total del estrato",
                     "",
                     "Resultado"),
      align = c("l", "r", "c", "r", "c", "r", "c", "r"),
      caption = "Denominador: no vacunados hospitalizados x total de vacunados / Total del estrato") %>%
  kable_styling(
    bootstrap_options = c("bordered", "striped", "hover", "condensed"),
    full_width = FALSE, position = "center"
  ) %>%
  row_spec(0, background = "#FFCCBC", color = "#2C3E50", bold = TRUE) %>%
  row_spec(1, background = "#FFF1B6") %>%
  row_spec(2, background = "#FFE0B2") %>%
  row_spec(3, background = "#FFCCBC", bold = TRUE, color = "#BF360C") %>%
  column_spec(1, bold = TRUE, width = "10em", color = "#2C3E50") %>%
  column_spec(2, background = "#FFD8BE", color = "#2C3E50") %>%
  column_spec(3, width = "1.5em", bold = TRUE) %>%
  column_spec(4, background = "#CDE7FF", color = "#2C3E50") %>%
  column_spec(5, width = "1.5em", bold = TRUE) %>%
  column_spec(6, background = "#E4C1F9", color = "#2C3E50") %>%
  column_spec(7, width = "1.5em", bold = TRUE) %>%
  column_spec(8, bold = TRUE, color = "#BF360C", background = "#FFF1B6")
Denominador: no vacunados hospitalizados x total de vacunados / Total del estrato
Estrato de edad No vacunados hospitalizados Total de vacunados Total del estrato Resultado
Menores de 50 anios 2,416 x 85,407 / 334,210 = 617.4062
50 anios o mas 664 x 71,991 / 80,542 = 593.5043
SUMA (denominador) = 1210.9105

El denominador pondera, para cada estrato, los no vacunados hospitalizados por el tamaño del grupo vacunado.Se observa que el estrato de menores de 50 años aporta 617.4062 y el de 50 años o más aporta 593.5043, sumando un denominador total de 1,210.91.

A diferencia del numerador, comnsidero que aquí ambos estratos contribuyen de forma más equilibrada, lo cual refleja que la hospitalización en no vacunados es proporcionalmente alta tanto en jóvenes como en mayores, pues el hecho de que el denominador (1,210.91) sea más del doble que el numerador (539.8053) ya me anticipa que el Riesgo Relativo ajustado será sustancialmente menor que 1, confirmando el efecto protector de la vacunación.

# ============================================================================
# ** TABLA 4: RESULTADO FINAL DEL RIESGO RELATIVO DE MANTEL-HAENSZEL **
# ============================================================================

tabla_final_rr <- data.frame(
  Elemento = c(
    "Numerador (suma de ambos estratos)",
    "Denominador (suma de ambos estratos)",
    "Riesgo Relativo de Mantel-Haenszel = numerador / denominador",
    "Efectividad vacunal ajustada = (1 - Riesgo Relativo) x 100"
  ),
  Valor = c(
    numerador_RR_MH,
    denominador_RR_MH,
    RR_MH,
    VE_ajustada
  )
)

kable(tabla_final_rr,
      col.names = c("Elemento del calculo", "Valor obtenido"),
      align = c("l", "r"),
      caption = "Resultado del Riesgo Relativo ajustado por edad (metodo de Mantel-Haenszel)") %>%
  kable_styling(
    bootstrap_options = c("bordered", "striped", "hover", "condensed"),
    full_width = FALSE, position = "center"
  ) %>%
  row_spec(0, background = "#B3E5FC", color = "#2C3E50", bold = TRUE) %>%
  row_spec(1, background = "#E2F0CB", color = "#1B5E20") %>%
  row_spec(2, background = "#FFCCBC", color = "#BF360C") %>%
  row_spec(3, background = "#B3E5FC", bold = TRUE, color = "#1A5276") %>%
  row_spec(4, background = "#FFF1B6", bold = TRUE, color = "#7D6608") %>%
  column_spec(1, bold = TRUE, width = "35em", color = "#2C3E50") %>%
  column_spec(2, bold = TRUE, width = "12em")
Resultado del Riesgo Relativo ajustado por edad (metodo de Mantel-Haenszel)
Elemento del calculo Valor obtenido
Numerador (suma de ambos estratos) 539.8052908
Denominador (suma de ambos estratos) 1210.9104661
Riesgo Relativo de Mantel-Haenszel = numerador / denominador 0.4457846
Efectividad vacunal ajustada = (1 - Riesgo Relativo) x 100 55.4215356

3.2.5 Análisis del Riesgo Relativo ajustado por edad

El Riesgo Relativo ajustado por Mantel-Haenszel que obtuve es de 0.4457846, lo cual significa que las personas vacunadas con esquema completo tienen un riesgo de hospitalización aproximadamente un 55.4% menor que las personas no vacunadas, cuando se comparan dentro del mismo grupo etario.

Este valor me permite calcular directamente la efectividad vacunal ajustada, pues al aplicar la fórmula (1 − Riesgo Relativo) × 100 obtengo una efectividad de 55.42%, lo cual indica que la vacunación previene aproximadamente 55 de cada 100 hospitalizaciones que habrían ocurrido en ausencia de vacunación.

El contraste con el Riesgo Relativo crudo de 1.253362 es nuevamente dramático, pues mientras la estimación cruda sugería falsamente que la vacunación incrementaba el riesgo de hospitalización en un 25.3%,la estimación ajustada revela una reducción del riesgo de 55.4%.

Esta inversión confirma, por segunda vía, la Paradoja de Simpson que identifiqué previamente con el Odds Ratio; además observo que el Riesgo Relativo ajustado (0.4457846) y el Odds Ratio ajustado (0.4369479) son numéricamente cercanos, lo cual creo que es esperable dado que la hospitalización es un evento relativamente infrecuente dentro de cada estrato, condición bajo la cual ambas medidas tienden a converger.

Sin embargo, según la literatura, el Riesgo Relativo tiene la ventaja de ser directamente interpretable como una comparación de riesgos, lo cual lo hace más útil para comunicar la magnitud del efecto protector de la vacunación en términos de salud pública.

3.2.6 Pregunta 3c. Diferencia de riesgo ajustada por edad (método de Mantel-Haenszel)

La diferencia de riesgo cuantifica en términos absolutos cuántas hospitalizaciones adicionales (o menos) se observan en vacunados en comparación con no vacunados, y a diferencia del Odds Ratio y del Riesgo Relativo, que son medidas relativas, la diferencia de riesgo me permite estimar directamente el número de eventos evitados por la vacunación.

La fórmula de la diferencia de riesgo de Mantel-Haenszel es la siguiente:

\[DR_{MH} = \frac{\text{Suma de (peso del estrato} \times \text{diferencia de riesgo del estrato)}}{\text{Suma de los pesos de todos los estratos}}\]

por lo cual, el peso de cada estrato se calcula como el producto del total de vacunados por el total de no vacunados, dividido entre el total del estrato, por tanto, este peso refleja el tamaño efectivo de la comparación en cada grupo de edad.

A continuación presento el proceso completo en tablas.

# ============================================================================
# ** TABLA 1: CALCULO DE LOS PESOS DE CADA ESTRATO **
#
# El peso se calcula como:
#   peso = (total de vacunados x total de no vacunados) / Total del estrato
# ============================================================================

tabla_pesos_dr <- data.frame(
  Estrato = c("Menores de 50 anios", "50 anios o mas"),
  Total_vac = c(format(a1 + b1, big.mark = ","), format(a2 + b2, big.mark = ",")),
  Simbolo1 = c("x", "x"),
  Total_novac = c(format(c1 + d1, big.mark = ","), format(c2 + d2, big.mark = ",")),
  Simbolo2 = c("/", "/"),
  Total_estrato = c(format(T1, big.mark = ","), format(T2, big.mark = ",")),
  Simbolo3 = c("=", "="),
  Peso = c(sprintf("%.10f", w1), sprintf("%.10f", w2))
)

kable(tabla_pesos_dr,
      col.names = c("Estrato de edad",
                     "Total de vacunados",
                     "",
                     "Total de no vacunados",
                     "",
                     "Total del estrato",
                     "",
                     "Peso del estrato"),
      align = c("l", "r", "c", "r", "c", "r", "c", "r"),
      caption = "Pesos de cada estrato: total de vacunados x total de no vacunados / Total del estrato") %>%
  kable_styling(
    bootstrap_options = c("bordered", "striped", "hover", "condensed"),
    full_width = FALSE, position = "center"
  ) %>%
  row_spec(0, background = "#D1C4E9", color = "#2C3E50", bold = TRUE) %>%
  row_spec(1, background = "#E2F0CB") %>%
  row_spec(2, background = "#FFF1B6") %>%
  column_spec(1, bold = TRUE, width = "10em", color = "#2C3E50") %>%
  column_spec(2, background = "#CDE7FF", color = "#2C3E50") %>%
  column_spec(3, width = "1.5em", bold = TRUE) %>%
  column_spec(4, background = "#FFCCBC", color = "#2C3E50") %>%
  column_spec(5, width = "1.5em", bold = TRUE) %>%
  column_spec(6, background = "#E4C1F9", color = "#2C3E50") %>%
  column_spec(7, width = "1.5em", bold = TRUE) %>%
  column_spec(8, bold = TRUE, color = "#4A148C", background = "#D1C4E9")
Pesos de cada estrato: total de vacunados x total de no vacunados / Total del estrato
Estrato de edad Total de vacunados Total de no vacunados Total del estrato Peso del estrato
Menores de 50 anios 85,407 x 248,803 / 334,210 = 63581.3345531253
50 anios o mas 71,991 x 8,551 / 80,542 = 7643.1556330858

El peso de cada estrato refleja el tamaño efectivo de la comparación entre vacunados y no vacunados dentro de ese grupo de edad, Y se obtiene multiplicando el total de vacunados por el total de no vacunados y dividiendo entre el total del estrato.

El estrato de menores de 50 años tiene un peso de 63,581.335, mientras que el de 50 años o más tiene un peso de apenas 7,643.156. Esta diferencia de casi nueve veces considero que es porque, el estrato joven tiene un número mucho mayor de personas en ambos grupos de comparación, lo cual le otorga mayor influencia en la estimación ajustada final.

# ============================================================================
# ** TABLA 2a: CALCULO DE LOS RIESGOS ESTRATO-ESPECIFICOS **
#
# El riesgo en cada grupo se obtiene dividiendo los hospitalizados
# entre el total de personas de ese grupo dentro del estrato.
# ============================================================================

tabla_riesgos <- data.frame(
  Estrato = c("Menores de 50 anios", "Menores de 50 anios",
              "50 anios o mas", "50 anios o mas"),
  Grupo = c("Vacunados", "No vacunados",
            "Vacunados", "No vacunados"),
  Hospitalizados = c(format(a1, big.mark = ","), format(c1, big.mark = ","),
                     format(a2, big.mark = ","), format(c2, big.mark = ",")),
  Simbolo1 = c("/", "/", "/", "/"),
  Total_grupo = c(format(a1 + b1, big.mark = ","), format(c1 + d1, big.mark = ","),
                  format(a2 + b2, big.mark = ","), format(c2 + d2, big.mark = ",")),
  Simbolo2 = c("=", "=", "=", "="),
  Riesgo = c(riesgo_vac_e1, riesgo_novac_e1,
             riesgo_vac_e2, riesgo_novac_e2)
)

kable(tabla_riesgos,
      col.names = c("Estrato de edad",
                     "Grupo",
                     "Hospitalizados",
                     "",
                     "Total del grupo",
                     "",
                     "Riesgo (proporcion)"),
      align = c("l", "l", "r", "c", "r", "c", "r"),
      caption = "Riesgos estrato-especificos: hospitalizados / total del grupo") %>%
  kable_styling(
    bootstrap_options = c("bordered", "striped", "hover", "condensed"),
    full_width = FALSE, position = "center"
  ) %>%
  row_spec(0, background = "#B3E5FC", color = "#2C3E50", bold = TRUE) %>%
  row_spec(1, background = "#E2F0CB") %>%
  row_spec(2, background = "#FFCCBC") %>%
  row_spec(3, background = "#E2F0CB") %>%
  row_spec(4, background = "#FFCCBC") %>%
  column_spec(1, bold = TRUE, width = "10em", color = "#2C3E50") %>%
  column_spec(2, bold = TRUE, color = "#2C3E50") %>%
  column_spec(3, background = "#CDE7FF", color = "#2C3E50") %>%
  column_spec(4, width = "1.5em", bold = TRUE) %>%
  column_spec(5, background = "#E4C1F9", color = "#2C3E50") %>%
  column_spec(6, width = "1.5em", bold = TRUE) %>%
  column_spec(7, bold = TRUE, color = "#1A5276", background = "#FFF1B6")
Riesgos estrato-especificos: hospitalizados / total del grupo
Estrato de edad Grupo Hospitalizados Total del grupo Riesgo (proporcion)
Menores de 50 anios Vacunados 453 / 85,407 = 0.0053040
Menores de 50 anios No vacunados 2,416 / 248,803 = 0.0097105
50 anios o mas Vacunados 1,908 / 71,991 = 0.0265033
50 anios o mas No vacunados 664 / 8,551 = 0.0776517

Cada riesgo lo he obtenido dividiendo el número de hospitalizados entre el total de personas de ese grupo dentro del estrato. En el estrato de menores de 50 años, el riesgo en vacunados es 453 / 85,407 = 0.0053040, mientras que el riesgo en no vacunados es 2,416 / 248,803 = 0.0097105.

En el estrato de 50 años o más, el riesgo en vacunados es 1,908 / 71,991 = 0.0265033 y el riesgo en no vacunados es 664 / 8,551 = 0.0776517.

En ambos estratos, el riesgo en no vacunados es mayor que en vacunados, lo cual me confirma que la vacunación es protectora dentro de cada grupo de edad, además, los riesgos del estrato de 50 años o más son considerablemente mayores que los del estrato joven, lo cual refleja que la edad es un factor de riesgo independiente para la hospitalización.

# ============================================================================
# ** TABLA 2: DIFERENCIAS DE RIESGO ESTRATO-ESPECIFICAS **
# ============================================================================

tabla_dr_estratos <- data.frame(
  Estrato = c("Menores de 50 anios", "50 anios o mas"),
  Riesgo_vac = c(riesgo_vac_e1, riesgo_vac_e2),
  Simbolo1 = c("-", "-"),
  Riesgo_novac = c(riesgo_novac_e1, riesgo_novac_e2),
  Simbolo2 = c("=", "="),
  DR = c(DR_e1, DR_e2)
)

kable(tabla_dr_estratos,
      col.names = c("Estrato de edad",
                     "Riesgo en vacunados",
                     "",
                     "Riesgo en no vacunados",
                     "",
                     "Diferencia de riesgo"),
      align = c("l", "r", "c", "r", "c", "r"),
      caption = "Diferencias de riesgo estrato-especificas: riesgo en vacunados - riesgo en no vacunados") %>%
  kable_styling(
    bootstrap_options = c("bordered", "striped", "hover", "condensed"),
    full_width = FALSE, position = "center"
  ) %>%
  row_spec(0, background = "#B3E5FC", color = "#2C3E50", bold = TRUE) %>%
  row_spec(1, background = "#E2F0CB") %>%
  row_spec(2, background = "#FFF1B6") %>%
  column_spec(1, bold = TRUE, width = "10em", color = "#2C3E50") %>%
  column_spec(2, background = "#CDE7FF", color = "#2C3E50") %>%
  column_spec(3, width = "1.5em", bold = TRUE) %>%
  column_spec(4, background = "#FFCCBC", color = "#2C3E50") %>%
  column_spec(5, width = "1.5em", bold = TRUE) %>%
  column_spec(6, bold = TRUE, color = "#1B5E20", background = "#E2F0CB")
Diferencias de riesgo estrato-especificas: riesgo en vacunados - riesgo en no vacunados
Estrato de edad Riesgo en vacunados Riesgo en no vacunados Diferencia de riesgo
Menores de 50 anios 0.0053040
0.0097105 = -0.0044065
50 anios o mas 0.0265033
0.0776517 = -0.0511484
Las diferencias de riesgo estrato-específicas confirman que en ambos grupos de edad la vacunación reduce el riesgo absoluto de hospitalización, dado que ambos valores son negativos.

En el estrato de menores de 50 años, la diferencia es de -0.0044065, lo cual significa que por cada persona vacunada en este grupo etario, el riesgo de hospitalización se reduce en 0.44 puntos porcentuales con respecto a una persona no vacunada de la misma edad.

Esta magnitud relativamente pequeña considero que la razón es porque tanto vacunados como no vacunados jóvenes tienen riesgos basales bajos de hospitalización (0.0053040 y 0.0097105 respectivamente), de modo que la diferencia absoluta entre ambos es modesta.

En el estrato de 50 años o más, la diferencia es de -0.0511484, lo cual representa una reducción absoluta considerablemente mayor de 5.11 puntos porcentuales menos de riesgo de hospitalización en vacunados frente a no vacunados.

Esta magnitud casi doce veces superior a la del estrato joven, considero que una razón es porque los adultos mayores tienen un riesgo basal de hospitalización mucho más elevado (0.0776517 en no vacunados frente a 0.0265033 en vacunados), por lo cual el beneficio absoluto de la vacunación es proporcionalmente mayor.

considero que este hallazgo, tiene una implicación directa para la salud pública, pues aunque la vacunación protege en todas las edades, el impacto absoluto de vacunar a una persona mayor es considerablemente mayor que el de vacunar a una persona joven, lo cual respalda epidemiológicamente la decisión del Reino Unido de haber priorizado a los adultos mayores en las primeras fases de la campaña vacunal.

# ============================================================================
# ** PASO 3: NUMERADOR DE LA DIFERENCIA DE RIESGO DE MANTEL-HAENSZEL **
#
# Lo obtengo multiplicando el peso de cada estrato por su diferencia
# de riesgo correspondiente y sumando ambos productos.
# ============================================================================

prod1 <- w1 * DR_e1
prod2 <- w2 * DR_e2

tabla_num_dr <- data.frame(
  Estrato = c("Menores de 50 anios", "50 anios o mas", "SUMA (numerador)"),
  Peso = c(sprintf("%.10f", w1), sprintf("%.10f", w2), ""),
  Simbolo1 = c("x", "x", ""),
  DR = c(sprintf("%.16f", DR_e1), sprintf("%.16f", DR_e2), ""),
  Simbolo2 = c("=", "=", "="),
  Resultado = c(sprintf("%.10f", prod1), sprintf("%.10f", prod2), sprintf("%.10f", numerador_DR_MH))
)

kable(tabla_num_dr,
      col.names = c("Estrato de edad",
                     "Peso del estrato",
                     "",
                     "Diferencia de riesgo",
                     "",
                     "Resultado"),
      align = c("l", "r", "c", "r", "c", "r"),
      caption = "Numerador: peso del estrato x diferencia de riesgo del estrato") %>%
  kable_styling(
    bootstrap_options = c("bordered", "striped", "hover", "condensed"),
    full_width = FALSE, position = "center"
  ) %>%
  row_spec(0, background = "#B9F6CA", color = "#2C3E50", bold = TRUE) %>%
  row_spec(1, background = "#E2F0CB") %>%
  row_spec(2, background = "#EAFAF1") %>%
  row_spec(3, background = "#B9F6CA", bold = TRUE, color = "#1B5E20") %>%
  column_spec(1, bold = TRUE, width = "10em", color = "#2C3E50") %>%
  column_spec(2, background = "#D1C4E9", color = "#2C3E50") %>%
  column_spec(3, width = "1.5em", bold = TRUE) %>%
  column_spec(4, background = "#FFD8BE", color = "#2C3E50") %>%
  column_spec(5, width = "1.5em", bold = TRUE) %>%
  column_spec(6, bold = TRUE, color = "#1B5E20", background = "#E2F0CB")
Numerador: peso del estrato x diferencia de riesgo del estrato
Estrato de edad Peso del estrato Diferencia de riesgo Resultado
Menores de 50 anios 63581.3345531253 x -0.0044064789511348 = -280.1698123934
50 anios o mas 7643.1556330858 x -0.0511484237248753 = -390.9353629162
SUMA (numerador) = -671.1051753096

El numerador lo he obtenido multiplicando el peso de cada estrato por su diferencia de riesgo y sumando ambos productos, por ejemplo,en el estrato de menores de 50 años, el producto es 63,581.3345531253 × (-0.0044064789511348) = -280.1698123934, y en el de 50 años o más es 7,643.1556330858 × (-0.0511484237248753) = -390.9353629162. La suma de ambos me da un numerador de -671.1051753096.

Los dos productos son negativos porque los pesos son siempre positivos (reflejan tamaños de estrato) mientras que las diferencias de riesgo son negativas en ambos grupos de edad (la vacunación reduce el riesgo), además observo que, a pesar de que el estrato joven tiene un peso casi nueve veces mayor, el estrato de 50 años o más aporta un producto de mayor magnitud absoluta (-390.94 frente a -280.17), lo cual se debe a que la diferencia de riesgo en los adultos mayores es considerablemente más grande (-0.0511 frente a -0.0044).

Considero que esto refleja que el beneficio absoluto de la vacunación es mucho mayor en el grupo de mayor riesgo basal.

# ============================================================================
# ** PASO 4: DENOMINADOR (SUMA DE PESOS) **
# ============================================================================

tabla_den_dr <- data.frame(
  Estrato = c("Menores de 50 anios", "50 anios o mas", "SUMA (denominador)"),
  Peso = c(sprintf("%.10f", w1), sprintf("%.10f", w2), sprintf("%.10f", denominador_DR_MH))
)

kable(tabla_den_dr,
      col.names = c("Estrato de edad",
                     "Peso del estrato"),
      align = c("l", "r"),
      caption = "Denominador: la suma de los pesos de ambos estratos") %>%
  kable_styling(
    bootstrap_options = c("bordered", "striped", "hover", "condensed"),
    full_width = FALSE, position = "center"
  ) %>%
  row_spec(0, background = "#FFCCBC", color = "#2C3E50", bold = TRUE) %>%
  row_spec(1, background = "#FFF1B6") %>%
  row_spec(2, background = "#FFE0B2") %>%
  row_spec(3, background = "#FFCCBC", bold = TRUE, color = "#BF360C") %>%
  column_spec(1, bold = TRUE, width = "10em", color = "#2C3E50") %>%
  column_spec(2, bold = TRUE, color = "#2C3E50")
Denominador: la suma de los pesos de ambos estratos
Estrato de edad Peso del estrato
Menores de 50 anios 63581.3345531253
50 anios o mas 7643.1556330858
SUMA (denominador) 71224.4901862111
# ============================================================================
# ** PASO 5: RESULTADO FINAL DE LA DIFERENCIA DE RIESGO DE MH **


tabla_final_dr <- data.frame(
  Elemento = c(
    "Numerador (suma ponderada)",
    "Denominador (suma de pesos)",
    "Diferencia de riesgo de Mantel-Haenszel = numerador / denominador",
    "Diferencia de riesgo por cada 100 personas vacunadas"
  ),
  Valor = c(
    format(numerador_DR_MH, digits = 22),
    format(denominador_DR_MH, digits = 22),
    format(DR_MH, digits = 22),
    format(DR_MH * 100, digits = 22)
  )
)

kable(tabla_final_dr,
      col.names = c("Elemento del calculo", "Valor obtenido"),
      align = c("l", "r"),
      caption = "Resultado de la diferencia de riesgo ajustada por edad (metodo de Mantel-Haenszel)") %>%
  kable_styling(
    bootstrap_options = c("bordered", "striped", "hover", "condensed"),
    full_width = FALSE, position = "center"
  ) %>%
  row_spec(0, background = "#B3E5FC", color = "#2C3E50", bold = TRUE) %>%
  row_spec(1, background = "#E2F0CB", color = "#1B5E20") %>%
  row_spec(2, background = "#FFCCBC", color = "#BF360C") %>%
  row_spec(3, background = "#B3E5FC", bold = TRUE, color = "#1A5276") %>%
  row_spec(4, background = "#FFF1B6", bold = TRUE, color = "#7D6608") %>%
  column_spec(1, bold = TRUE, width = "35em", color = "#2C3E50") %>%
  column_spec(2, bold = TRUE, width = "18em")
Resultado de la diferencia de riesgo ajustada por edad (metodo de Mantel-Haenszel)
Elemento del calculo Valor obtenido
Numerador (suma ponderada) -671.1051753096478478255
Denominador (suma de pesos) 71224.4901862111146329
Diferencia de riesgo de Mantel-Haenszel = numerador / denominador -0.009422393527213650601948
Diferencia de riesgo por cada 100 personas vacunadas -0.9422393527213650532559

3.2.7 Análisis de la Diferencia de Riesgo ajustada por edad

La diferencia de riesgo ajustada por Mantel-Haenszel que obtuve es de -0.009422393527, lo cual equivale a -0.942239352721 hospitalizaciones menos por cada 100 personas vacunadas en comparación con personas no vacunadas del mismo grupo etario.

Este valor negativo confirma, por tercera vía independiente, que la vacunación reduce el riesgo absoluto de hospitalización una vez que se controla la confusión por edad.

El contraste con la estimación cruda es, una vez más, revelador, pues la diferencia de riesgo cruda fue positiva (0.003032, equivalente a 0.30 hospitalizaciones adicionales por cada 100 personas vacunadas), lo cual sugería falsamente que la vacunación incrementaba el riesgo absoluto de hospitalización, pero cuando ajusté por edad, la dirección se invirtió por completo y la diferencia pasó a ser negativa.

Con las tres medidas ajustadas por el método de Mantel-Haenszel que ya calculé, como fué el Odds Ratio de 0.4369479 (reducción de momios del 56.3%), el Riesgo Relativo de 0.4457846 (efectividad vacunal del 55.4%) y una Diferencia de Riesgo de -0.009422 (0.94 hospitalizaciones menos por cada 100 personas vacunadas), entonces puedo decir que la evidencia es consistente y contundente, en dirección a que la vacunación si confiere protección sustancial contra la hospitalización por variante Delta, y la confusión por edad era la responsable de invertir artificialmente la dirección de la asociación en el análisis crudo.

3.2.8 Verificación cruzada con el paquete epiR

Una vez que ya he completado los cálculos manuales de las tres medidas de asociación ajustadas por Mantel-Haenszel, procedo a verificar los resultados utilizando la función epi.2by2 del paquete epiR, pues haciendo esta verificación puedo confirmar si he aplicado correctamente las fórmulas y asi mismo puedo obtener los intervalos de confianza del 95%, los cuales necesito para evaluar la significancia estadística de las asociaciones.

# ============================================================================
# ** VERIFICACION CON epi.2by2 DEL PAQUETE epiR **
#
# para epi.2by2 le ordeno este formato:
#                    Enfermo(Si)   Enfermo(No)
#   Expuesto(Si)        a              b
#   Expuesto(No)        c              d
#
#  el orden correcto
# de los valores es: a, c, b, d 
# ============================================================================

tabla_estratificada <- array(
  c(a1, c1, b1, d1,
    a2, c2, b2, d2),
  dim = c(2, 2, 2),
  dimnames = list(
    "Vacunacion"      = c("Vacunado", "No_vacunado"),
    "Hospitalizacion" = c("Si", "No"),
    "Edad"            = c("Menor_50", "Mayor_igual_50")
  )
)

resultado_mh <- epi.2by2(tabla_estratificada, method = "cohort.count")

or_mh_est   <- resultado_mh$massoc.detail$OR.mh.wald$est
or_mh_lower <- resultado_mh$massoc.detail$OR.mh.wald$lower
or_mh_upper <- resultado_mh$massoc.detail$OR.mh.wald$upper

rr_mh_est   <- resultado_mh$massoc.detail$RR.mh.wald$est
rr_mh_lower <- resultado_mh$massoc.detail$RR.mh.wald$lower
rr_mh_upper <- resultado_mh$massoc.detail$RR.mh.wald$upper

tabla_verif <- data.frame(
  Medida = c(
    "Odds Ratio de Mantel-Haenszel",
    "Riesgo Relativo de Mantel-Haenszel"
  ),
  Manual = c(
    format(OR_MH, digits = 22),
    format(RR_MH, digits = 22)
  ),
  epiR = c(
    format(or_mh_est, digits = 22),
    format(rr_mh_est, digits = 22)
  ),
  IC_inferior = c(
    format(or_mh_lower, digits = 22),
    format(rr_mh_lower, digits = 22)
  ),
  IC_superior = c(
    format(or_mh_upper, digits = 22),
    format(rr_mh_upper, digits = 22)
  )
)

kable(tabla_verif,
      col.names = c("Medida de asociacion",
                     "Calculo manual",
                     "Verificacion epi.2by2",
                     "Limite inferior IC 95%",
                     "Limite superior IC 95%"),
      align = c("l", "r", "r", "r", "r"),
      caption = "Verificacion cruzada: calculo manual vs. paquete epiR") %>%
  kable_styling(
    bootstrap_options = c("bordered", "striped", "hover", "condensed"),
    full_width = FALSE, position = "center"
  ) %>%
  row_spec(0, background = "#B3E5FC", color = "#2C3E50", bold = TRUE) %>%
  row_spec(1, background = "#E2F0CB") %>%
  row_spec(2, background = "#FFF1B6") %>%
  column_spec(1, bold = TRUE, width = "16em", color = "#2C3E50") %>%
  column_spec(2, background = "#CDE7FF", color = "#2C3E50") %>%
  column_spec(3, background = "#D1C4E9", color = "#2C3E50") %>%
  column_spec(4, background = "#FFCCBC", color = "#2C3E50") %>%
  column_spec(5, background = "#FFCCBC", color = "#2C3E50")
Verificacion cruzada: calculo manual vs. paquete epiR
Medida de asociacion Calculo manual Verificacion epi.2by2 Limite inferior IC 95% Limite superior IC 95%
Odds Ratio de Mantel-Haenszel 0.4369478556062452723197 0.4369478556062452723197 0.4065210469424800354687 0.4696520142188617130685
Riesgo Relativo de Mantel-Haenszel 0.445784643796061197385 0.445784643796061197385 0.4167672919411271048595 0.4768223238412218623061

Los valores obtenidos mediante el cálculo manual coinciden exactamente con los arrojados por la función epi.2by2 del paquete epiR, lo cual me valida la correcta aplicación de las fórmulas del método de Mantel-Haenszel.

El Odds Ratio ajustado es de 0.4369478 con un intervalo de confianza del 95% de [0.4065210, 0.4696520], y el Riesgo Relativo ajustado es de 0.4457846 con un intervalo de confianza del 95% de [0.4167672,0.4768223].

En ambos casos, los intervalos de confianza no incluyen el valor nulo de 1, lo cual confirma que las asociaciones ajustadas son estadísticamente significativas,es decir, puedo afirmar con un 95% de confianza que la vacunación reduce tanto los momios como el riesgo de hospitalización por variante Delta, una vez que se controla la confusión por edad.

3.2.9 Pregunta 3d. Interpretación global ajustada por edad de la “eficacia” vacunal

¿Estoy de acuerdo con los resultados ajustados? Sí, considero que las estimaciones ajustadas por el método de Mantel-Haenszel son confiables y reflejan de manera más veraz la relación entre vacunación y hospitalización por variante Delta.

Encuentro que las tres medidas ajustadas apuntan en la misma dirección de forma consistente, pues el Odds Ratio de 0.4369478 indica una reducción del 56.3% en los momios de hospitalización, el Riesgo Relativo de 0.4457846 indica una reducción del 55.4% en el riesgo de hospitalización (lo cual equivale a una efectividad vacunal ajustada del 55.4%), y la Diferencia de Riesgo de -0.009422 indica que por cada 100 personas vacunadas se evitan aproximadamente 0.94 hospitalizaciones en comparación con personas no vacunadas del mismo grupo etario.

Considero que estos resultados ajustados son creíbles porque, en primer lugar, la dirección de la asociación es biológicamente plausible, dado que las vacunas contra la enfermedad por coronavirus 2019 fueron diseñadas para generar una respuesta inmune que reduzca la probabilidad de enfermedad severa, y los ensayos clínicos controlados demostraron eficacias superiores al 90% contra hospitalización.

En segundo lugar, la consistencia entre las tres medidas refuerza la robustez del hallazgo, pues cuando el Odds Ratio, el Riesgo Relativo y la Diferencia de Riesgo convergen en la misma dirección y con magnitudes coherentes, la probabilidad de que el resultado se deba a un artefacto estadístico se reduce considerablemente.

En tercer lugar, los intervalos de confianza del 95% que obtuve mediante la verificación con el paquete epiR no incluyen el valor nulo en ninguna de las medidas, lo cual me confirma la significancia estadística de las asociaciones.

Finalmente, considero que la comparación entre las estimaciones crudas y las ajustadas constituye la lección epidemiológica central de este ejercicio, pues el análisis crudo me sugería falsamente que la vacunación incrementaba el riesgo de hospitalización (Odds Ratio crudo de 1.2572, Riesgo Relativo crudo de 1.2534, Diferencia de Riesgo cruda positiva de 0.003032).

Sin embargo, esta distorsión era producida por la confusión por edad, que actuaba como un confusor clásico al estar asociada tanto con la exposición (priorización vacunal por edad) como con el desenlace (mayor riesgo de hospitalización en adultos mayores), sin encontrarse en la vía causal.

Cuando ajusté por edad mediante el método de Mantel-Haenszel, la dirección de la asociación se invirtió por completo, revelando el verdadero efecto protector de la vacunación, por lo tanto, esta inversión es la manifestación cuantitativa de la Paradoja de Simpson, y demuestra por qué en epidemiología nunca se deben interpretar medidas de asociación crudas sin antes evaluar la presencia de confusores.

3.2.10 Pregunta 3e. Grafo Acíclico Dirigido (DAG)

A continuación construyo un Grafo Acíclico Dirigido que representa la estructura causal entre las tres variables involucradas en este análisis, como es la vacunación (exposición), la visita al servicio de emergencias con hospitalización (desenlace) y la edad (confusor).

Para ello utilizo los paquetes dagitty y ggdag en R, lo cual me permite definir y visualizar las relaciones causales.

library(dagitty)
library(ggdag)
library(ggplot2)

dag_vacuna <- dagitty('dag {
  Vacuna [exposure, pos="0,0"]
  Urgencias [outcome, pos="6,0"]
  Edad [adjusted, pos="3,4"]
  Edad -> Vacuna
  Edad -> Urgencias
  Vacuna -> Urgencias
}')

cat("Adjustment sets:\n")
## Adjustment sets:
print(adjustmentSets(dag_vacuna, exposure = "Vacuna", outcome = "Urgencias"))
## { Edad }
cat("\nPaths:\n")
## 
## Paths:
print(paths(dag_vacuna, from = "Vacuna", to = "Urgencias"))
## $paths
## [1] "Vacuna -> Urgencias"         "Vacuna <- Edad -> Urgencias"
## 
## $open
## [1] TRUE TRUE
tidy_dag <- tidy_dagitty(dag_vacuna)

tidy_dag$data$rol <- ifelse(tidy_dag$data$name == "Vacuna", "Exposicion",
                     ifelse(tidy_dag$data$name == "Urgencias", "Desenlace",
                     ifelse(tidy_dag$data$name == "Edad", "Confusor (ajustado)", NA)))

ggplot(tidy_dag$data, aes(x = x, y = y, xend = xend, yend = yend)) +
  geom_dag_edges(edge_width = 0.8,
    arrow_directed = grid::arrow(length = grid::unit(10, "pt"), type = "closed")) +
  geom_dag_point(aes(fill = rol), size = 20, shape = 21,
                 color = "black", stroke = 1.3) +
  geom_dag_text(aes(label = name), size = 5, color = "black", fontface = "bold") +
  scale_fill_manual(
    name = "Rol en el DAG",
    values = c("Exposicion" = "#93C5FD",
               "Desenlace"  = "#FDBA74",
               "Confusor (ajustado)" = "#86EFAC")) +
  annotate("text", x = 3, y = -0.55,
           label = "Causal path (efecto directo)",
           size = 3.8, fontface = "bold.italic", color = "#1E40AF") +
  annotate("text", x = 0.7, y = 2.8,
           label = "Backdoor path\nEdad -> Vacuna\n(priorizacion por edad)",
           size = 3.2, fontface = "bold", color = "#DC2626") +
  annotate("text", x = 5.3, y = 2.8,
           label = "Backdoor path\nEdad -> Urgencias\n(riesgo basal por edad)",
           size = 3.2, fontface = "bold", color = "#DC2626") +
  annotate("label", x = 3, y = -2.0,
           label = "Backdoor path:  Vacuna <- Edad -> Urgencias\nSolucion: ajustar por Edad (Mantel-Haenszel)",
           size = 3.5, fill = "grey95", color = "black",
           label.padding = unit(0.5, "lines"), fontface = "bold") +
  coord_cartesian(xlim = c(-1.5, 7.5), ylim = c(-2.2, 5.5)) +
  theme_dag() +
  labs(title = "Grafico Aciclico Dirigido (DAG)",
       subtitle = "Estructura causal: vacunacion, edad y hospitalizacion por variante Delta") +
  theme(
    plot.title    = element_text(face = "bold", size = 16, hjust = 0.5),
    plot.subtitle = element_text(size = 11, hjust = 0.5, color = "grey40"),
    legend.position = "bottom",
    legend.title = element_text(face = "bold", size = 11),
    legend.text  = element_text(size = 10),
    plot.margin = margin(10, 30, 10, 30))

El DAG que obtuve confirma la estructura causal que identifiqué a lo largo de este ejercicio. La función adjustmentSets() de dagitty devuelve { Edad } como el conjunto mínimo suficiente de ajuste, lo cual valida formalmente que la edad es la única variable que necesito condicionar para obtener una estimación no sesgada del efecto causal de la vacunación sobre la hospitalización.

En la grafica se observa como la edad actúa como causa común de la exposición y del desenlace, y se identifican dos rutas de asociación entre vacunación y urgencias:

La primera es el camino causal directo (Vacuna → Urgencias), que representa el efecto biológico de interés, es decir, la capacidad de la vacuna para reducir la probabilidad de hospitalización por variante Delta.

La segunda es el backdoor path (Vacuna ← Edad → Urgencias),que constituye una ruta no causal a través de la cual se genera una asociación espuria.

Esta ruta existe porque la edad influye simultáneamente en la probabilidad de estar vacunado (las personas mayores fueron priorizadas en la campaña vacunal del Reino Unido) y en la probabilidad de requerir hospitalización (el riesgo basal de enfermedad severa aumenta con la edad).

Mientras este backdoor path permanezca abierto, la asociación cruda entre vacunación y hospitalización mezcla el efecto causal verdadero con la confusión generada por la edad, lo cual explica por qué las medidas crudas que calculé en la Sección 2 (OR de 1.2572, RR de 1.2533) sugerían falsamente que la vacunación incrementaba el riesgo.

Al ajustar por edad mediante el método de Mantel-Haenszel, lo que hice fue bloquear el backdoor path, condicionando sobre la variable confusora, por tanto, esto aisló el camino causal directo y me permitió obtener las estimaciones ajustadas (OR de 0.43694, RR de 0.44578), que reflejan el efecto protector real de la vacunación.

4 Otros Sesgos Potenciales

En esta sección analizo la Tabla 6 del reporte, que presenta las tasas de ataque secundario para todas las variantes del SARS-CoV-2 entre el 5 de enero y el 24 de agosto de 2021, con datos de variante actualizados al 13 de septiembre de 2021 y datos de rastreo de contactos actualizados al 14 de septiembre de 2021.

Las tasas de ataque secundario cuantifican la proporción de contactos de un caso índice que desarrollan la infección, y constituyen una medida fundamental para evaluar la transmisibilidad de cada variante.

La tabla distingue entre casos vinculados a viajes internacionales y casos no vinculados a viajes, y entre contactos intradomiciliarios y contactos no domiciliarios.

A continuación presento los datos de la variante Delta extraídos de la Tabla 6, que es la variante de interés.

library(kableExtra)
library(tibble)

delta_tabla6 <- tibble(
  Indicador = c(
    "Casos vinculados a viajes (con contactos)",
    "Casos no vinculados a viajes (con contactos)",
    "Proporcion vinculada a viajes",
    "Tasa de ataque en contactos de casos vinculados a viajes (IC 95%)",
    "Tasa de ataque en contactos intradomiciliarios (no viaje o desconocido) (IC 95%)",
    "Tasa de ataque en contactos no domiciliarios (no viaje o desconocido) (IC 95%)"
  ),
  Valor = c(
    "11,059 (63.2% con contactos)",
    "471,476 (72.6% con contactos)",
    "2.3%",
    "1.8% (1.8% - 1.9%) [2,454/133,378]",
    "10.4% (10.4% - 10.5%) [81,552/781,192]",
    "5.4% (5.3% - 5.5%) [13,956/259,671]"
  )
)

kable(delta_tabla6, 
      caption = "Datos de la variante Delta extraidos de la Tabla 6 del reporte",
      align = c("l", "r")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, font_size = 13) %>%
  row_spec(0, background = "#E8D5F5", bold = TRUE, color = "black") %>%
  row_spec(1:2, background = "#FFF0F5") %>%
  row_spec(3, background = "#F0FFF0") %>%
  row_spec(4:6, background = "#F5F0FF")
Datos de la variante Delta extraidos de la Tabla 6 del reporte
Indicador Valor
Casos vinculados a viajes (con contactos) 11,059 (63.2% con contactos)
Casos no vinculados a viajes (con contactos) 471,476 (72.6% con contactos)
Proporcion vinculada a viajes 2.3%
Tasa de ataque en contactos de casos vinculados a viajes (IC 95%) 1.8% (1.8% - 1.9%) [2,454/133,378]
Tasa de ataque en contactos intradomiciliarios (no viaje o desconocido) (IC 95%) 10.4% (10.4% - 10.5%) [81,552/781,192]
Tasa de ataque en contactos no domiciliarios (no viaje o desconocido) (IC 95%) 5.4% (5.3% - 5.5%) [13,956/259,671]

4.1 Pregunta 4a. Robustez de las tasas de ataque y las proporciones

A partir de la información presentada en la Tabla 6 y de las notas metodológicas al pie de la misma, considero que las tasas de ataque secundario y las proporciones reportadas tienen precisión razonable en los grupos con grandes denominadores, especialmente para la variante Delta, pero no pueden considerarse plenamente robustas como estimaciones absolutas de la transmisión real.

Esto se debe a que el propio reporte señala que las tasas provenientes del sistema NHS Test and Trace deben interpretarse como límites inferiores, dado que la naturaleza del rastreo de contactos y de las pruebas diagnósticas impide captar todos los casos secundarios.

Además, el reporte reconoce que algunos casos vinculados a viajes pueden no ser identificados y quedar clasificados como no relacionados con viajes o como desconocidos, lo que introduce incertidumbre adicional en las proporciones de casos asociados a viajes.

En segundo lugar, la decisión de marcar como “Unavailable” las celdas con menos de 50 contactos o menos de 20 casos protege frente a estimaciones muy inestables derivadas de números pequeños, aunque al mismo tiempo limita la evaluación comparativa de ciertas variantes o subgrupos con baja frecuencia.

Asimismo, el recorte temporal de los datos hasta el 24 de agosto de 2021 mejora la consistencia del seguimiento, porque permite que los contactos tengan tiempo suficiente para convertirse en casos, pero también hace que los conteos sean menores que los de otras fuentes más amplias.

Por último, en la variante Delta se observa que algunas tasas tuvieron intervalos de confianza estrechos, como ocurre en los contactos intradomiciliarios no vinculados a viajes o desconocidos, lo que sugiere buena precisión estadística por el gran tamaño muestral, sin embargo, esa precisión no elimina la posibilidad de que hayan tenido subregistros ni clasificación imperfecta, por lo cual mi conclusión, es que estas estimaciones son útiles y relativamente estables para describir patrones comparativos, pero probablemente hayan subestimado la magnitud real de la transmisión secundaria.

4.2 Pregunta 4b. Sesgos potenciales adicionales

A partir de la revisión de la Tabla 6 y de las notas metodológicas al pie de la misma, identifico la presencia de sesgos potenciales adicionales que pudieron haber afectado las estimaciones de las tasas de ataque secundario.

En primer lugar, reconozco la posibilidad de sesgo de información, especialmente en forma de clasificación errónea, porque la medición tanto del vínculo con viajes como de la ocurrencia real de casos secundarios no fue perfecta, pues esto ocurrió debido a que el propio reporte señaló que algunos casos vinculados a viajes internacionales podían no haber sido identificados por los métodos de vigilancia utilizados y, por tanto, pudieron haber quedado clasificados erróneamente como casos no vinculados a viajes o como casos de origen desconocido.

De esta manera, entiendo que la clasificación de la exposición no fue completamente exacta, además, observo que la identificación de los casos secundarios dependió del funcionamiento del sistema de rastreo y de pruebas diagnósticas, por lo cual no todos los contactos infectados necesariamente fueron detectados como casos secundarios.

Esto pudo haber sucedido, por ejemplo, si algunos contactos permanecieron asintomáticos, si no consultaron, o si no se realizaron una prueba durante el periodo de observación, por lo cual, en ese escenario, interpreto que el desenlace también pudo haber sido medido de manera incompleta, lo que habría llevado a subestimar la transmisión secundaria real.

Así mismo, comprendo que la identificación de contactos dependió del autorreporte y del éxito operativo del rastreo, de tal forma que algunos contactos no declarados o no localizados ni siquiera ingresaron al conjunto observado, lo que también afectó la calidad de la información disponible.

Por otra parte, también identifico la posibilidad de sesgo de selección, porque la inclusión en el análisis no dependió únicamente de que hubiera ocurrido la transmisión, sino de que el caso índice hubiera sido detectado, notificado y rastreado exitosamente por el sistema NHS Test and Trace, además de haber cumplido con los criterios temporales definidos por el estudio.

Considero que esto significó que los casos y contactos finalmente analizados no necesariamente representaron todo el universo real de transmisiones secundarias, sino solamente la fracción que logró ingresar al sistema de vigilancia bajo esas condiciones.

Por eso interpreto que la muestra analizada pudo haber sido selectiva y que las tasas observadas probablemente reflejaron una parte del fenómeno, pero no su magnitud completa.

En conjunto, concluyo que las estimaciones de la Tabla 6 pudieron haber estado afectadas por clasificación errónea del vínculo con viajes, por subdetección de casos secundarios y por selección incompleta de contactos rastreados, lo cual me lleva a interpretar estas tasas con cautela.

Por esa razón considero que los valores reportados eran útiles para describir patrones comparativos de transmisión, pero que debían entenderse más como estimaciones conservadoras o límites inferiores que como una medición exhaustiva de la transmisión secundaria real.

4.3 Pregunta 4c. Efecto del sesgo de información sobre las tasas de ataque y dirección del sesgo

Comienzo diferenciando sesgo de información no diferencial y sesgo de información diferencial, porque cada uno modifica las tasas de ataque en una dirección distinta y, por tanto, cambia también la manera en que interpreto lo que sucedió en la Tabla 6.

En primer lugar, entiendo que ocurrió un mecanismo compatible con sesgo de información no diferencial en la detección de los casos secundarios, por tanto lo planteo así porque muchos contactos que sí pudieron haberse infectado no necesariamente quedaron registrados como casos, ya fuera porque permanecieron asintomáticos, porque no consultaron o porque no se realizaron una prueba diagnóstica durante el periodo de observación.

En ese contexto, lo que sucedió fue que parte de los contactos verdaderamente infectados terminaron clasificados como no infectados, y esa pérdida de información pudo haber afectado de manera relativamente parecida a los distintos grupos comparados.

A partir de ello, las tasas de ataque observadas probablemente quedaron subestimadas frente a la transmisión secundaria real, en consecuencia, dado que este tipo de error no dependió necesariamente del grupo de comparación, considero que su efecto más probable fue hacia el valor nulo, es decir, las diferencias reales entre grupos tendieron a atenuarse.

En otras palabras, la transmisión secundaria pudo haber sido mayor que la reportada, pero que el sistema de observación no logró captarla de forma completa.

En segundo lugar, identifico también un mecanismo plausible de sesgo de información diferencial, especialmente en la clasificación del vínculo con viajes internacionales.

Lo fundamento en que el propio reporte describió que los casos relacionados con viajes fueron investigados mediante múltiples fuentes de vigilancia, incluyendo EpiCell, equipos de protección sanitaria, registros del NHS Test and Trace asociados con viajes internacionales, notificaciones por Regulaciones Sanitarias Internacionales y equipos internacionales de rastreo.

Mientras tanto, los casos no vinculados a viajes dependieron principalmente del sistema rutinario de rastreo, por tanto, desde esa perspectiva, entiendo que no todos los grupos fueron observados con la misma intensidad, y que eso pudo haber generado una probabilidad distinta de clasificar correctamente tanto la exposición como la detección de contactos secundarios.

Lo que sucedió, entonces, fue que los casos vinculados a viajes pudieron haber tenido una búsqueda más exhaustiva, mientras que en otros grupos la captación pudo haber sido menos completa, y por esa razón, considero que las tasas de ataque comparadas entre grupos no necesariamente provinieron de procesos de observación equivalentes.

A partir de esta diferencia en la intensidad del rastreo, el sesgo diferencial pudo haber distorsionado las comparaciones de una manera más compleja y potencialmente alejada del valor nulo,es decir, que parte de las diferencias observadas entre casos vinculados a viajes y casos no vinculados a viajes o desconocidos pudieron no deberse únicamente a diferencias reales en transmisión, sino también a diferencias en la forma en que los casos y contactos fueron identificados, clasificados y seguidos.

En conjunto,en la Tabla 6 sucedieron dos fenómenos simultáneos, por un lado, ocurrió una subdetección general de casos secundarios, que probablemente empujó las tasas de ataque hacia abajo y redujo las diferencias entre grupos.

Por otro lado, pudo haber ocurrido una captación diferencial según el vínculo con viajes, que hizo menos comparables algunos grupos entre sí, es por eso, que las tasas de ataque reportadas fueron útiles para describir patrones generales de transmisión, pero no reflejaron de manera exacta toda la magnitud de la transmisión secundaria real.

4.4 Pregunta 4d. Métodos de corrección del sesgo de información en análisis epidemiológicos

A partir de lo que identifiqué en las preguntas anteriores, considero pertinente mencionar los principales métodos que se han utilizado para corregir o cuantificar el sesgo de información en análisis epidemiológicos.

En primer lugar, los estudios de validación permiten comparar la clasificación obtenida por el sistema habitual con un estándar de referencia de mayor exactitud.

En este caso, si se hubiera seleccionado una muestra de contactos y se les hubiera aplicado una prueba diagnóstica independientemente de la presencia de síntomas, habría sido posible estimar la sensibilidad y la especificidad del sistema NHS Test and Trace, y con ello cuantificar de manera más precisa el subregistro de casos secundarios.

En segundo lugar, el análisis cuantitativo de sesgos permite incorporar de forma explícita los errores de clasificación dentro del cálculo de las tasas o medidas de asociación.

Si se dispone de valores razonables de sensibilidad y especificidad, las estimaciones pueden recalcularse corrigiendo la proporción esperada de falsos negativos y falsos positivos,pues este método puede aplicarse de forma determinística, usando valores fijos, o de forma probabilística, asignando distribuciones a esos parámetros para reflejar la incertidumbre.

En tercer lugar, el análisis de sensibilidad permite explorar cómo cambian las estimaciones bajo distintos supuestos plausibles de error de medición. En un escenario como el de la Tabla 6, habría sido útil recalcular las tasas de ataque suponiendo diferentes niveles de subdetección de casos secundarios, con el fin de valorar qué tan estables eran las conclusiones frente a ese posible sesgo.

Finalmente, la imputación múltiple puede ser útil cuando el sesgo de información se acompaña de datos faltantes, por ejemplo, si parte de la información no estuvo disponible porque algunos contactos no fueron identificados o porque no se registró adecuadamente su desenlace.

Este método permite completar varios conjuntos de datos bajo supuestos explícitos sobre la ausencia de información y combinar luego los resultados en una sola estimación final.

En conjunto, considero que, para un escenario como el de la Tabla 6,los estudios de validación, el análisis cuantitativo de sesgos y el análisis de sensibilidad habrían sido los métodos más pertinentes, porque hubieran podido aproximarse de manera más realista a la magnitud verdadera de la transmisión secundaria, reconociendo que la medición observada había sido incompleta.

4.5 Pregunta 4e. Posibilidad de sesgo de selección

Sí considero que existió posibilidad de sesgo de selección en la Tabla 6, porque la inclusión en el análisis no dependió únicamente de que hubiera ocurrido la transmisión secundaria, sino también de que el caso índice y sus contactos hubieran sido identificados, rastreados e incorporados al sistema NHS Test and Trace.

Por lo tanto, considero que la muestra finalmente analizada no necesariamente representó todo el universo real de transmisiones, sino solo la fracción que logró entrar al sistema de vigilancia.

En primer lugar, considero que operó un mecanismo de captación incompleta de casos índice y de sus contactos, pues el propio reporte reconoció que algunos casos vinculados a viajes internacionales pudieron no haber sido identificados por los métodos empleados.**

De este modo, una parte de los eventos de transmisión nunca ingresó al análisis, por lo que la población observada quedó reducida de manera no aleatoria, y en consecuencia, la muestra analizada pudo diferir sistemáticamente de la población objetivo.

En segundo lugar, considero que existió una intensidad diferencial del rastreo según el tipo de caso, pues los casos relacionados con viajes internacionales fueron investigados mediante múltiples fuentes complementarias, mientras que los casos no vinculados a viajes dependieron principalmente del sistema rutinario.

Por lo que, la probabilidad de que un contacto fuera identificado y quedara incluido en el análisis no fue la misma en todos los grupos, y ello hizo menos comparables las tasas de ataque observadas entre viajeros y no viajeros.

En cuanto a la dirección del sesgo, considero que no fue simple ni completamente predecible, precisamente porque se trató de una selección diferencial, pues si en un grupo se captaron proporcionalmente más contactos, más casos secundarios, o ambos, la estimación pudo desplazarse en una u otra dirección.

Sin embargo, sí considero que este sesgo tuvo la capacidad de alejar la estimación del valor nulo o distorsionar la comparabilidad entre grupos, haciendo que parte de las diferencias observadas no reflejara únicamente la transmisión real, sino también la forma en que los individuos fueron seleccionados para quedar observados.

En conjunto, considero que sí hubo posibilidad de sesgo de selección, porque la entrada de casos y contactos al análisis dependió de procesos de identificación y rastreo que no ocurrieron de manera uniforme entre los grupos comparados.

Por ello, interpreto que las tasas de ataque de la Tabla 6 deben leerse con cautela, pues pudieron estar influidas no solo por diferencias reales de transmisibilidad, sino también por diferencias en la probabilidad de quedar incluidos en el sistema de vigilancia.

Referencias

  1. UK Health Security Agency. COVID-19 vaccine surveillance report, Week 41. London: UKHSA; 14 October 2021. Disponible en: https://assets.publishing.service.gov.uk/government/uploads/system/uploads/attachment_data/file/1025358/Vaccine-surveillance-report-week-41.pdf