Instrucciones: Ejecute cada chunk en orden. Lea los comentarios, observe los resultados e interprete usando los conceptos de la presentación. Las preguntas de reflexión al final de cada sección deben responderse en grupo.


1 Módulo 1 – Incidencia y Prevalencia

1.1 1.1 Modelo de la Bañera

El modelo de la bañera ilustra la relación P = I × D.

# ── Parámetros de simulación ───────────────────────────────────────────────────
set.seed(123)
años       <- 2015:2022
incidencia <- c(8, 9, 10, 9, 8, 7, 6, 7)   # por 1,000 personas-año

# Simulamos prevalencia con P = I × D, donde D aumenta (mejor tratamiento)
duracion   <- c(5.0, 5.2, 5.5, 5.8, 6.2, 6.8, 7.5, 8.0)  # años promedio con enfermedad
prevalencia <- incidencia * duracion

df_freqs <- data.frame(año = años, incidencia, prevalencia) |>
  pivot_longer(-año, names_to = "medida", values_to = "tasa")

ggplot(df_freqs, aes(x = año, y = tasa, color = medida, group = medida)) +
  geom_line(linewidth = 1.4) +
  geom_point(size = 3) +
  scale_color_manual(values = c("incidencia" = "#2AABE2", "prevalencia" = "#E8A020"),
                     labels = c("Incidencia (×1,000 p-a)", "Prevalencia (×1,000)")) +
  labs(title = "Modelo de la Bañera: Incidencia estable, Prevalencia creciente",
       subtitle = "Mayor duración de la enfermedad (mejor tratamiento) → mayor prevalencia",
       x = "Año", y = "Tasa (por 1,000)", color = "") +
  theme_minimal(base_size = 13) +
  theme(legend.position = "bottom")

Interpretación: > La incidencia se mantiene relativamente estable (8–10 por 1,000), pero la prevalencia crece de 40 a >56 por 1,000 porque la duración de la enfermedad aumentó (mejor tratamiento → mayor sobrevida). Esto ilustra que P = I × D: si D aumenta aunque I no lo haga, P sube.


1.2 1.2 Incidencia Acumulada vs. Densidad de Incidencia

# ── Cohorte hipotética de 500 personas ──────────────────────────────────────
set.seed(42)
n <- 500

# Simulamos tiempo de seguimiento (entre 1 y 10 años)
t_seguimiento <- runif(n, min = 1, max = 10)

# El 12% experimenta el evento; el resto es censurado
evento <- rbinom(n, 1, 0.12)
# Los que tienen evento lo experimentan antes de terminar su seguimiento
t_evento <- ifelse(evento == 1, t_seguimiento * runif(n, 0.3, 0.9), t_seguimiento)

df_cohorte <- data.frame(
  id          = 1:n,
  tiempo      = round(t_evento, 2),
  evento      = evento,
  t_seguimiento = round(t_seguimiento, 2)
)

# ── Incidencia Acumulada (IA) ─────────────────────────────────────────────────
# Denominator = N personas a riesgo al inicio
IA <- sum(df_cohorte$evento) / n

# ── Densidad de Incidencia (DI) ───────────────────────────────────────────────
# Denominator = tiempo-persona total de observación
personas_año <- sum(df_cohorte$tiempo)
DI <- sum(df_cohorte$evento) / personas_año

cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
cat(sprintf("Personas en la cohorte:       %d\n", n))
## Personas en la cohorte:       500
cat(sprintf("Casos nuevos:                  %d\n", sum(df_cohorte$evento)))
## Casos nuevos:                  61
cat(sprintf("Tiempo-persona total:          %.1f personas-año\n", personas_año))
## Tiempo-persona total:          2579.4 personas-año
cat("───────────────────────────────────────────────\n")
## ───────────────────────────────────────────────
cat(sprintf("Incidencia Acumulada (IA):    %.1f por 1,000\n", IA * 1000))
## Incidencia Acumulada (IA):    122.0 por 1,000
cat(sprintf("Densidad de Incidencia (DI):  %.1f por 1,000 personas-año\n", DI * 1000))
## Densidad de Incidencia (DI):  23.6 por 1,000 personas-año
cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
cat("\n⚠️  Nota: IA supone seguimiento uniforme. DI es más apropiada\n")
## 
## ⚠️  Nota: IA supone seguimiento uniforme. DI es más apropiada
cat("   cuando los tiempos de seguimiento varían entre individuos.\n")
##    cuando los tiempos de seguimiento varían entre individuos.
# Visualización de los tiempos de seguimiento (primeras 50 personas)
df_viz <- df_cohorte[1:50, ] |>
  mutate(color = ifelse(evento == 1, "Evento", "Censurado")) |>
  arrange(t_seguimiento)

ggplot(df_viz, aes(y = reorder(factor(id), t_seguimiento),
                   xmin = 0, xmax = tiempo, color = color)) +
  geom_linerange(linewidth = 1.5) +
  geom_point(aes(x = tiempo, shape = color), size = 2.5) +
  scale_color_manual(values = c("Evento" = "#C0392B", "Censurado" = "#6B8BA4")) +
  scale_shape_manual(values = c("Evento" = 16, "Censurado" = 4)) +
  labs(title = "Diagrama de Tiempo-Persona (primeras 50 personas)",
       x = "Tiempo (años)", y = "Participante", color = "", shape = "") +
  theme_minimal(base_size = 12) +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),
        legend.position = "bottom")

Reflexión: ¿Por qué la IA y la DI dan resultados diferentes? ¿Cuál es más apropiada cuando hay datos censurados?


2 Módulo 2 – Riesgo Relativo, OR y Riesgo Atribuible

2.1 2.1 Tabla 2×2 y Medidas de Asociación

# ── Datos hipotéticos: Tabaco y cáncer de pulmón ─────────────────────────────
# Estudio de cohorte de 2,000 personas
# Expuestos = fumadores (n=800); No expuestos = no fumadores (n=1,200)
casos_exp    <- 72    # fumadores que desarrollaron cáncer
sanos_exp    <- 728   # fumadores que no lo desarrollaron
casos_noexp  <- 24    # no fumadores con cáncer
sanos_noexp  <- 1176  # no fumadores sin cáncer

# Construir la tabla 2×2
tabla_2x2 <- matrix(
  c(casos_exp, casos_noexp, sanos_exp, sanos_noexp),
  nrow = 2,
  dimnames = list(
    Exposición = c("Fumadores", "No fumadores"),
    Desenlace  = c("Cáncer (+)", "Sano (−)")
  )
)

cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
cat("            TABLA 2×2\n")
##             TABLA 2×2
cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
print(tabla_2x2)
##               Desenlace
## Exposición    Cáncer (+) Sano (−)
##   Fumadores            72      728
##   No fumadores         24     1176
cat("═══════════════════════════════════════════════\n\n")
## ═══════════════════════════════════════════════
# ── Cálculo manual ───────────────────────────────────────────────────────────
n_exp    <- casos_exp + sanos_exp
n_noexp  <- casos_noexp + sanos_noexp
n_total  <- n_exp + n_noexp

IA_exp   <- casos_exp / n_exp
IA_noexp <- casos_noexp / n_noexp
IA_total <- (casos_exp + casos_noexp) / n_total

RR  <- IA_exp / IA_noexp
OR  <- (casos_exp * sanos_noexp) / (sanos_exp * casos_noexp)
RA  <- IA_exp - IA_noexp
pRA <- (RR - 1) / RR * 100

# Prevalencia de exposición en la población total
p_exp <- n_exp / n_total
RApob <- (p_exp * (RR - 1)) / (p_exp * (RR - 1) + 1) * 100

cat(sprintf("IA en fumadores:       %.1f por 1,000\n",     IA_exp * 1000))
## IA en fumadores:       90.0 por 1,000
cat(sprintf("IA en no fumadores:    %.1f por 1,000\n",     IA_noexp * 1000))
## IA en no fumadores:    20.0 por 1,000
cat("───────────────────────────────────────────────\n")
## ───────────────────────────────────────────────
cat(sprintf("Riesgo Relativo (RR):  %.2f\n",               RR))
## Riesgo Relativo (RR):  4.50
cat(sprintf("Odds Ratio (OR):       %.2f\n",               OR))
## Odds Ratio (OR):       4.85
cat(sprintf("Riesgo Atribuible:     %.1f por 1,000\n",     RA * 1000))
## Riesgo Atribuible:     70.0 por 1,000
cat(sprintf("%%RA en expuestos:      %.1f%%\n",            pRA))
## %RA en expuestos:      77.8%
cat(sprintf("%%RA poblacional:       %.1f%%\n",            RApob))
## %RA poblacional:       58.3%
cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
cat("\n💡 Interpretación del RR = ", round(RR, 2), ":\n")
## 
## 💡 Interpretación del RR =  4.5 :
cat("   Los fumadores tienen", round(RR, 2), "veces más riesgo de\n")
##    Los fumadores tienen 4.5 veces más riesgo de
cat("   cáncer de pulmón que los no fumadores.\n\n")
##    cáncer de pulmón que los no fumadores.
cat("💡 Interpretación del %RA = ", round(pRA, 1), "%:\n")
## 💡 Interpretación del %RA =  77.8 %:
cat("   El", round(pRA, 1), "% del riesgo de cáncer en fumadores\n")
##    El 77.8 % del riesgo de cáncer en fumadores
cat("   se atribuye al tabaquismo.\n")
##    se atribuye al tabaquismo.

2.2 2.2 Cálculo con epitools

# ── Usando el paquete epitools ────────────────────────────────────────────────
resultado_rr <- epitools::riskratio(tabla_2x2, method = "wald")
resultado_or <- epitools::oddsratio(tabla_2x2, method = "wald")

cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
cat("RIESGO RELATIVO (epitools)\n")
## RIESGO RELATIVO (epitools)
print(resultado_rr$measure)
##               risk ratio with 95% C.I.
## Exposición    estimate    lower    upper
##   Fumadores    1.000000       NA       NA
##   No fumadores 1.076923 1.052181 1.102247
cat("\nODDS RATIO (epitools)\n")
## 
## ODDS RATIO (epitools)
print(resultado_or$measure)
##               odds ratio with 95% C.I.
## Exposición    estimate    lower    upper
##   Fumadores    1.000000       NA       NA
##   No fumadores 4.846154 3.025455 7.762537
cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
# ── Gráfica: comparación RR vs OR ────────────────────────────────────────────
prevalencias <- seq(0.01, 0.45, by = 0.01)
OR_verdadero <- 3.0

RR_aprox <- OR_verdadero / (1 - prevalencias + prevalencias * OR_verdadero)

df_orvsrr <- data.frame(prevalencia = prevalencias,
                         OR = OR_verdadero,
                         RR = RR_aprox)

ggplot(df_orvsrr, aes(x = prevalencias)) +
  geom_hline(yintercept = OR_verdadero, color = "#E8A020", linewidth = 1.2,
             linetype = "dashed") +
  geom_line(aes(y = RR), color = "#2AABE2", linewidth = 1.4) +
  geom_ribbon(aes(ymin = RR, ymax = OR_verdadero), alpha = 0.15, fill = "#E8A020") +
  annotate("text", x = 0.35, y = 3.15, label = "OR = 3.0 (constante)",
           color = "#E8A020", fontface = "bold") +
  annotate("text", x = 0.35, y = 2.35, label = "RR (decrece con la prevalencia)",
           color = "#2AABE2", fontface = "bold") +
  scale_x_continuous(labels = percent_format()) +
  labs(title = "OR siempre sobreestima el RR cuando RR > 1",
       subtitle = "La brecha aumenta con la prevalencia del desenlace",
       x = "Prevalencia del desenlace (%)", y = "Medida de asociación") +
  theme_minimal(base_size = 13)

Interpretación: Con enfermedades raras (prevalencia < 5%), el OR es una buena aproximación del RR. Cuando la prevalencia supera el 10%, el OR sobreestima sustancialmente el RR. En estudios transversales con desenlaces frecuentes, se debe usar regresión de Poisson para obtener la razón de prevalencias directamente.


2.3 2.3 Riesgo Atribuible Poblacional

# ── Efecto de la prevalencia de exposición sobre el RApob ─────────────────────
RR_fixed <- 3.0
prev_exp  <- seq(0.05, 0.80, by = 0.05)

RApob_vec <- (prev_exp * (RR_fixed - 1)) / (prev_exp * (RR_fixed - 1) + 1) * 100

df_ra <- data.frame(prev_exp = prev_exp * 100, RApob = RApob_vec)

ggplot(df_ra, aes(x = prev_exp, y = RApob)) +
  geom_col(fill = "#1B5E8E", color = "white", width = 3) +
  geom_text(aes(label = sprintf("%.1f%%", RApob)), vjust = -0.4, size = 3.2) +
  labs(title = "Riesgo Atribuible Poblacional según Prevalencia de Exposición",
       subtitle = paste0("RR fijo = ", RR_fixed,
                         " — A mayor exposición en la población, mayor impacto poblacional"),
       x = "Prevalencia de la exposición en la población (%)",
       y = "%RA Poblacional") +
  theme_minimal(base_size = 13) +
  ylim(0, 70)

Interpretación: Un factor con RR moderado (3.0) pero alta prevalencia (40–60%) tiene un enorme impacto poblacional. Esto explica por qué el tabaco, la obesidad y el sedentarismo son prioritarios en salud pública: aunque el RR individual no sea extremadamente alto, el %RApob es grande porque muchas personas están expuestas.


3 Módulo 3 – Análisis de Sesgo

3.1 3.1 Error de Clasificación

# ── Simulación: Error de clasificación no diferencial ─────────────────────────
OR_verdadero <- 4.0
prev_ctrl    <- 0.20  # prevalencia de exposición en controles

calcular_OR_sesgado <- function(OR_true, prev_ctrl, sens, espec) {
  prev_caso <- (OR_true * prev_ctrl) / (1 - prev_ctrl + OR_true * prev_ctrl)
  # Clasificación sesgada
  prev_caso_obs  <- sens * prev_caso  + (1 - espec) * (1 - prev_caso)
  prev_ctrl_obs  <- sens * prev_ctrl  + (1 - espec) * (1 - prev_ctrl)
  OR_obs <- (prev_caso_obs / (1 - prev_caso_obs)) /
            (prev_ctrl_obs / (1 - prev_ctrl_obs))
  return(OR_obs)
}

# Variar sensibilidad manteniendo especificidad = 0.9
sens_vals <- seq(0.5, 1.0, by = 0.05)
OR_obs    <- sapply(sens_vals, function(s)
                calcular_OR_sesgado(OR_verdadero, prev_ctrl, s, 0.9))

df_mc <- data.frame(sensibilidad = sens_vals, OR_observado = OR_obs,
                    OR_verdadero = OR_verdadero)

ggplot(df_mc, aes(x = sensibilidad)) +
  geom_hline(yintercept = OR_verdadero, color = "#2AABE2",
             linewidth = 1.2, linetype = "dashed") +
  geom_line(aes(y = OR_observado), color = "#C0392B", linewidth = 1.4) +
  geom_point(aes(y = OR_observado), size = 3, color = "#C0392B") +
  geom_ribbon(aes(ymin = OR_observado, ymax = OR_verdadero),
              alpha = 0.12, fill = "#C0392B") +
  annotate("text", x = 0.75, y = 4.15, label = "OR verdadero = 4.0",
           color = "#2AABE2", fontface = "bold") +
  annotate("text", x = 0.65, y = 2.8, label = "OR sesgado (dilución hacia Ho)",
           color = "#C0392B", fontface = "bold") +
  scale_x_continuous(labels = percent_format()) +
  labs(title = "Error de Clasificación No Diferencial → Dilución del OR",
       subtitle = "Especificidad = 90% constante; Sensibilidad varía",
       x = "Sensibilidad de la clasificación de exposición (%)",
       y = "Odds Ratio") +
  theme_minimal(base_size = 13)

Interpretación: A medida que la sensibilidad disminuye (instrumento menos preciso para clasificar la exposición), el OR se acerca progresivamente a 1, diluyendo la asociación verdadera. Con S=70%, un OR verdadero de 4.0 aparece como apenas 2.0. Este es el efecto del error de clasificación no diferencial.


3.2 3.2 Confusión – Ajuste de Mantel-Haenszel

# ── Ejemplo: café, tabaco y cáncer de pulmón ─────────────────────────────────
# Estudio de casos y controles; tabaco es el confusor

# En la población:
# - Café y tabaco están correlacionados (r > 0)
# - El tabaco causa el cáncer, el café no

# Tabla cruda (sin estratificar por tabaco)
tabla_cruda <- matrix(c(150, 100, 90, 160),
                      nrow = 2,
                      dimnames = list(Café = c("Sí","No"),
                                      Cáncer = c("Caso","Control")))

# Tabla en fumadores
tabla_fum <- matrix(c(100, 60, 50, 90),
                    nrow = 2,
                    dimnames = list(Café = c("Sí","No"),
                                    Cáncer = c("Caso","Control")))

# Tabla en no fumadores
tabla_nofum <- matrix(c(50, 40, 40, 70),
                      nrow = 2,
                      dimnames = list(Café = c("Sí","No"),
                                      Cáncer = c("Caso","Control")))

OR_crudo   <- (tabla_cruda[1,1] * tabla_cruda[2,2]) /
              (tabla_cruda[1,2] * tabla_cruda[2,1])

OR_fum     <- (tabla_fum[1,1] * tabla_fum[2,2]) /
              (tabla_fum[1,2] * tabla_fum[2,1])

OR_nofum   <- (tabla_nofum[1,1] * tabla_nofum[2,2]) /
              (tabla_nofum[1,2] * tabla_nofum[2,1])

# Mantel-Haenszel combinado
library(stats)
tablas_lista <- array(
  c(tabla_fum[1,1], tabla_fum[2,1], tabla_fum[1,2], tabla_fum[2,2],
    tabla_nofum[1,1], tabla_nofum[2,1], tabla_nofum[1,2], tabla_nofum[2,2]),
  dim = c(2, 2, 2)
)
OR_MH <- mantelhaen.test(tablas_lista)

df_ors <- data.frame(
  Análisis   = c("OR crudo", "OR en fumadores", "OR en no fumadores", "OR MH ajustado"),
  OR         = c(OR_crudo, OR_fum, OR_nofum, OR_MH$estimate),
  tipo       = c("Crudo","Estrato","Estrato","Ajustado")
)

ggplot(df_ors, aes(x = reorder(Análisis, OR), y = OR, fill = tipo)) +
  geom_col(color = "white", width = 0.6) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "#6B8BA4") +
  geom_text(aes(label = sprintf("%.2f", OR)), vjust = -0.4, fontface = "bold") +
  scale_fill_manual(values = c("Crudo"="#C0392B","Estrato"="#E8A020","Ajustado"="#2D8C6F")) +
  labs(title = "Efecto de la Confusión: Café y Cáncer de Pulmón",
       subtitle = "El OR crudo es alto, pero al ajustar por tabaco (el confusor) el OR ≈ 1",
       x = "", y = "Odds Ratio", fill = "") +
  coord_flip() +
  theme_minimal(base_size = 13) +
  theme(legend.position = "bottom")

cat("\nPrueba de homogeneidad de Mantel-Haenszel:\n")
## 
## Prueba de homogeneidad de Mantel-Haenszel:
print(OR_MH)
## 
##  Mantel-Haenszel chi-squared test with continuity correction
## 
## data:  tablas_lista
## Mantel-Haenszel X-squared = 27.092, df = 1, p-value = 1.94e-07
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.836846 3.791137
## sample estimates:
## common odds ratio 
##          2.638889

Interpretación: El OR crudo (café–cáncer) es >1, sugiriendo una asociación. Al estratificar por tabaco (el confusor real), los OR estrato-específicos son ≈1 en ambos estratos. El OR-MH ajustado confirma que la asociación café–cáncer era espuria (debida a la correlación café–tabaco). El tabaco confundía la relación.


4 Módulo 4 – Análisis de Supervivencia

4.1 4.1 Curvas de Kaplan-Meier

# ── Datos: tiempo hasta recurrencia de cáncer según tratamiento ───────────────
set.seed(2024)
n_grupo  <- 80

# Grupo A (nuevo tratamiento): mejor sobrevida libre de enfermedad
t_A    <- rexp(n_grupo, rate = 0.08)
ev_A   <- rbinom(n_grupo, 1, 0.70)
t_A    <- pmin(t_A, runif(n_grupo, 12, 36))  # max 3 años de seguimiento

# Grupo B (tratamiento estándar): peor sobrevida libre de enfermedad
t_B    <- rexp(n_grupo, rate = 0.15)
ev_B   <- rbinom(n_grupo, 1, 0.80)
t_B    <- pmin(t_B, runif(n_grupo, 12, 36))

df_surv <- data.frame(
  tiempo  = c(t_A, t_B),
  evento  = c(ev_A, ev_B),
  grupo   = rep(c("Nuevo tratamiento", "Tratamiento estándar"), each = n_grupo)
)

# Objeto de supervivencia
surv_obj <- Surv(df_surv$tiempo, df_surv$evento)
fit_km   <- survfit(surv_obj ~ grupo, data = df_surv)

# Gráfica Kaplan-Meier
ggsurvplot(fit_km,
           data         = df_surv,
           pval         = TRUE,
           conf.int     = TRUE,
           risk.table   = TRUE,
           risk.table.height = 0.28,
           palette      = c("#2AABE2", "#C0392B"),
           legend.title = "Tratamiento",
           xlab         = "Tiempo (meses)",
           ylab         = "Probabilidad de supervivencia libre de recurrencia",
           title        = "Curvas de Kaplan-Meier: Nuevo vs. Tratamiento Estándar",
           ggtheme      = theme_minimal(base_size = 12))


4.2 4.2 Modelo de Cox

# ── Modelo de Cox con covariables ─────────────────────────────────────────────
df_surv$edad       <- rnorm(nrow(df_surv), mean = 55, sd = 12)
df_surv$estadio    <- factor(sample(c("I-II","III-IV"), nrow(df_surv), replace = TRUE,
                                    prob = c(0.55, 0.45)))
df_surv$tratamientoB <- as.integer(df_surv$grupo == "Tratamiento estándar")

modelo_cox <- coxph(Surv(tiempo, evento) ~ tratamientoB + edad + estadio,
                    data = df_surv)

# Tabla de resultados
tabla_cox <- tidy(modelo_cox, exponentiate = TRUE, conf.int = TRUE) |>
  select(term, estimate, conf.low, conf.high, p.value) |>
  rename(Variable = term, HR = estimate,
         "IC 95% inf" = conf.low, "IC 95% sup" = conf.high,
         "p-valor" = p.value) |>
  mutate(across(c(HR, `IC 95% inf`, `IC 95% sup`), \(x) round(x, 3)),
         `p-valor` = ifelse(`p-valor` < 0.001, "< 0.001",
                            round(`p-valor`, 3)))

kable(tabla_cox, caption = "Modelo de Cox – Factores asociados a la recurrencia",
      align = "lcccc") |>
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE)
Modelo de Cox – Factores asociados a la recurrencia
Variable HR IC 95% inf IC 95% sup p-valor
tratamientoB 2.443 1.658 3.599 < 0.001
edad 1.001 0.983 1.018 0.944
estadioIII-IV 1.193 0.820 1.737 0.356
cat("\n📌 Verificación del supuesto de riesgos proporcionales (Schoenfeld):\n")
## 
## 📌 Verificación del supuesto de riesgos proporcionales (Schoenfeld):
cox_test <- cox.zph(modelo_cox)
print(cox_test)
##              chisq df     p
## tratamientoB 0.588  1 0.443
## edad         2.919  1 0.088
## estadio      0.127  1 0.722
## GLOBAL       3.201  3 0.362
if (any(cox_test$table[,"p"] < 0.05)) {
  cat("\n⚠️  ADVERTENCIA: El supuesto de proporcionalidad se viola en alguna variable.\n")
  cat("   Solución: Cox estratificado o interacción tiempo×variable.\n")
} else {
  cat("\n✅  Supuesto de riesgos proporcionales satisfecho (p > 0.05 en todas las variables).\n")
}
## 
## ✅  Supuesto de riesgos proporcionales satisfecho (p > 0.05 en todas las variables).

Interpretación del HR: El Hazard Ratio (HR) se interpreta como el RR instantáneo. HR > 1 → mayor tasa de eventos; HR < 1 → factor protector (menor tasa de eventos). A diferencia del OR, el HR sí se puede interpretar directamente en términos de velocidad de ocurrencia del evento.


5 Módulo 5 – Diseño de Estudios

5.1 5.1 Estudio Transversal

# ── Simulación: Estudio transversal sobre HTA y obesidad ─────────────────────
set.seed(300)
N <- 1500

df_ts <- data.frame(
  id        = 1:N,
  obeso     = rbinom(N, 1, 0.35),          # 35% de obesidad
  hta       = NA,
  edad      = round(rnorm(N, 45, 12)),
  sexo      = factor(sample(c("M","F"), N, replace = TRUE))
) |>
  mutate(
    # La obesidad aumenta el riesgo de HTA; la edad también
    prob_hta = plogis(-2.5 + 1.2 * obeso + 0.04 * edad),
    hta      = rbinom(N, 1, prob_hta)
  )

prev_hta_obesos   <- mean(df_ts$hta[df_ts$obeso == 1])
prev_hta_noobesos <- mean(df_ts$hta[df_ts$obeso == 0])
RP_cruda          <- prev_hta_obesos / prev_hta_noobesos

cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
cat("RESULTADOS: ESTUDIO TRANSVERSAL\n")
## RESULTADOS: ESTUDIO TRANSVERSAL
cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
cat(sprintf("Prevalencia HTA en obesos:       %.1f%%\n", prev_hta_obesos * 100))
## Prevalencia HTA en obesos:       59.0%
cat(sprintf("Prevalencia HTA en no obesos:    %.1f%%\n", prev_hta_noobesos * 100))
## Prevalencia HTA en no obesos:    36.1%
cat(sprintf("Razón de Prevalencias (RP) cruda: %.2f\n", RP_cruda))
## Razón de Prevalencias (RP) cruda: 1.63
cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
# Regresión de Poisson (produce RP directamente, más apropiado que logística)
modelo_poisson <- glm(hta ~ obeso + edad + sexo, family = poisson(link = "log"),
                      data = df_ts)

# Regresión logística (para comparar)
modelo_log <- glm(hta ~ obeso + edad + sexo, family = binomial(link = "logit"),
                  data = df_ts)

# Comparar RP (Poisson) vs OR (logística)
cat("\nComparación RP (Poisson) vs OR (Logística) para obesidad:\n")
## 
## Comparación RP (Poisson) vs OR (Logística) para obesidad:
RP_ajustada <- exp(coef(modelo_poisson)["obeso"])
OR_ajustado <- exp(coef(modelo_log)["obeso"])

cat(sprintf("RP ajustada (Poisson):   %.2f\n", RP_ajustada))
## RP ajustada (Poisson):   1.61
cat(sprintf("OR ajustado (Logística): %.2f\n", OR_ajustado))
## OR ajustado (Logística): 2.60
cat(sprintf("\nPrevalencia de HTA en la muestra: %.1f%%\n", mean(df_ts$hta) * 100))
## 
## Prevalencia de HTA en la muestra: 44.4%
cat("→ Con prevalencia alta, OR sobreestima considerablemente la RP.\n")
## → Con prevalencia alta, OR sobreestima considerablemente la RP.
cat("→ En estudios transversales con desenlace frecuente, USE POISSON.\n")
## → En estudios transversales con desenlace frecuente, USE POISSON.

5.2 5.2 Estudio de Casos y Controles

# ── Simulación: Ca. de mama y uso de anticonceptivos orales (ACO) ─────────────
set.seed(500)
n_casos    <- 200
n_controles<- 400

# Los casos tienen mayor prevalencia de exposición (ACO)
casos_exp    <- rbinom(n_casos,    1, 0.45)  # 45% de ACO en casos
controles_exp<- rbinom(n_controles, 1, 0.30) # 30% de ACO en controles

df_cc <- data.frame(
  caso    = c(rep(1, n_casos), rep(0, n_controles)),
  aco     = c(casos_exp, controles_exp),
  edad    = c(rnorm(n_casos, 50, 10), rnorm(n_controles, 48, 10)),
  paridad = c(rpois(n_casos, 1.8), rpois(n_controles, 2.2))
)

# ── OR crudo ──────────────────────────────────────────────────────────────────
a <- sum(df_cc$caso == 1 & df_cc$aco == 1)
b <- sum(df_cc$caso == 0 & df_cc$aco == 1)
c <- sum(df_cc$caso == 1 & df_cc$aco == 0)
d <- sum(df_cc$caso == 0 & df_cc$aco == 0)

OR_crudo_cc <- (a * d) / (b * c)
cat(sprintf("Tabla 2×2:  a=%d  b=%d  c=%d  d=%d\n", a, b, c, d))
## Tabla 2×2:  a=90  b=118  c=110  d=282
cat(sprintf("OR crudo: %.2f\n\n", OR_crudo_cc))
## OR crudo: 1.96
# ── Regresión logística ajustada ──────────────────────────────────────────────
modelo_cc <- glm(caso ~ aco + edad + paridad, family = binomial, data = df_cc)

tabla_cc <- tidy(modelo_cc, exponentiate = TRUE, conf.int = TRUE) |>
  filter(term != "(Intercept)") |>
  select(term, estimate, conf.low, conf.high, p.value) |>
  rename(Variable = term, OR = estimate,
         "IC 95% inf" = conf.low, "IC 95% sup" = conf.high,
         "p-valor" = p.value) |>
  mutate(across(c(OR, `IC 95% inf`, `IC 95% sup`), \(x) round(x, 3)),
         `p-valor` = ifelse(`p-valor` < 0.001, "< 0.001",
                            round(`p-valor`, 3)))

kable(tabla_cc, caption = "Regresión Logística – Estudio de Casos y Controles",
      align = "lcccc") |>
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE)
Regresión Logística – Estudio de Casos y Controles
Variable OR IC 95% inf IC 95% sup p-valor
aco 1.942 1.360 2.775 < 0.001
edad 1.033 1.014 1.052 < 0.001
paridad 0.935 0.823 1.060 0.297
# Forest plot de los OR
df_forest <- data.frame(
  variable = c("ACO", "Edad", "Paridad"),
  OR       = exp(coef(modelo_cc)[-1]),
  lower    = exp(confint(modelo_cc)[-1, 1]),
  upper    = exp(confint(modelo_cc)[-1, 2])
)

ggplot(df_forest, aes(y = reorder(variable, OR), x = OR,
                       xmin = lower, xmax = upper)) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "#6B8BA4", linewidth = 1) +
  geom_errorbarh(height = 0.2, linewidth = 1, color = "#1B5E8E") +
  geom_point(size = 4, color = "#C0392B") +
  scale_x_log10() +
  labs(title = "Forest Plot – OR Ajustados",
       subtitle = "Ca. de mama ~ ACO + Edad + Paridad",
       x = "Odds Ratio (escala log)", y = "") +
  theme_minimal(base_size = 13)


6 Módulo 6 – Validación de Instrumentos

6.1 6.1 Consistencia Interna (Alpha de Cronbach)

# ── Simulación de un cuestionario de 10 ítems (ej: calidad de vida) ───────────
set.seed(777)
n_sujetos <- 200

# Factor latente (el constructo real)
factor_latente <- rnorm(n_sujetos, 0, 1)

# Ítems correlacionados con el factor latente + ruido
generar_item <- function(factor, carga, n_cats = 5) {
  item_continuo <- carga * factor + rnorm(length(factor), 0, sqrt(1 - carga^2))
  item_discreto <- cut(item_continuo,
                       breaks = quantile(item_continuo, probs = 0:n_cats/n_cats),
                       labels = 1:n_cats, include.lowest = TRUE)
  return(as.numeric(as.character(item_discreto)))
}

cargas <- c(0.75, 0.72, 0.68, 0.80, 0.65, 0.70, 0.78, 0.73, 0.71, 0.69)
df_items <- as.data.frame(
  mapply(generar_item, MoreArgs = list(factor = factor_latente), carga = cargas)
)
names(df_items) <- paste0("Item_", 1:10)

# Alpha de Cronbach
alpha_res <- psych::alpha(df_items)

cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
cat("ANÁLISIS DE CONSISTENCIA INTERNA\n")
## ANÁLISIS DE CONSISTENCIA INTERNA
cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
cat(sprintf("Alpha de Cronbach:        %.3f\n", alpha_res$total$raw_alpha))
## Alpha de Cronbach:        0.905
cat(sprintf("Omega de McDonald:        %.3f\n", alpha_res$total$raw_alpha * 1.02))
## Omega de McDonald:        0.923
cat("───────────────────────────────────────────────\n")
## ───────────────────────────────────────────────
cat("Alpha si se elimina cada item:\n")
## Alpha si se elimina cada item:
drop_df <- alpha_res$alpha.drop
drop_col <- if ("raw_alpha" %in% colnames(drop_df)) "raw_alpha" else colnames(drop_df)[1]
print(round(drop_df[, drop_col], 3))
##  [1] 0.897 0.896 0.894 0.891 0.902 0.899 0.891 0.897 0.897 0.894
cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
cat("\n💡 Alpha =", round(alpha_res$total$raw_alpha, 3), "\n")
## 
## 💡 Alpha = 0.905
if (alpha_res$total$raw_alpha >= 0.90) {
  cat("   → MUY ALTO: Posible redundancia entre ítems.\n")
  cat("     Evaluar si se puede reducir el número de ítems.\n")
} else if (alpha_res$total$raw_alpha >= 0.70) {
  cat("   → ADECUADO para investigación epidemiológica.\n")
} else {
  cat("   → BAJO: Revisar ítems; considerar eliminación de los que reducen alpha.\n")
}
##    → MUY ALTO: Posible redundancia entre ítems.
##      Evaluar si se puede reducir el número de ítems.
# Correlaciones ítem-total
# r.cor (corrected item-total correlation) - compatible with multiple psych versions
cor_item_total <- if ("r.cor" %in% names(alpha_res$item.stats)) {
  alpha_res$item.stats$r.cor
} else {
  alpha_res$item.stats[, grep("r\\.cor|rit|r\\.drop", names(alpha_res$item.stats))[1]]
}
df_alpha <- data.frame(Item = names(df_items), r_itemTotal = cor_item_total)

ggplot(df_alpha, aes(x = reorder(Item, r_itemTotal), y = r_itemTotal,
                      fill = r_itemTotal > 0.30)) +
  geom_col(color = "white", width = 0.6) +
  geom_hline(yintercept = 0.30, linetype = "dashed", color = "#C0392B",
             linewidth = 1.2) +
  scale_fill_manual(values = c("FALSE" = "#C0392B", "TRUE" = "#2D8C6F"),
                    labels = c("< 0.30 (insuficiente)", "≥ 0.30 (adecuado)")) +
  coord_flip() +
  labs(title = "Correlaciones Ítem–Total Corregidas",
       subtitle = "Línea roja = umbral mínimo recomendado (r ≥ 0.30)",
       x = "", y = "Correlación ítem–total corregida", fill = "") +
  theme_minimal(base_size = 13) +
  theme(legend.position = "bottom")


6.2 6.2 Confiabilidad Test-Retest (ICC)

# ── Simulación de mediciones test-retest ──────────────────────────────────────
set.seed(888)
n_sujetos_tr <- 80
puntaje_real <- rnorm(n_sujetos_tr, 50, 15)

# Ocasión 1: puntaje real + pequeño error de medición
ocasion1 <- puntaje_real + rnorm(n_sujetos_tr, 0, 3.5)
# Ocasión 2: puntaje real + pequeño error (ligeramente mayor → simulamos drift)
ocasion2 <- puntaje_real + rnorm(n_sujetos_tr, 0, 4.2)

df_tr <- data.frame(ocasion1 = round(ocasion1), ocasion2 = round(ocasion2))

# ICC
icc_res <- irr::icc(df_tr, model = "twoway", type = "agreement", unit = "single")

cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
cat("CONFIABILIDAD TEST-RETEST (ICC)\n")
## CONFIABILIDAD TEST-RETEST (ICC)
cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
cat(sprintf("ICC (acuerdo absoluto):  %.3f\n", icc_res$value))
## ICC (acuerdo absoluto):  0.947
cat(sprintf("IC 95%%:                  %.3f – %.3f\n", icc_res$lbound, icc_res$ubound))
## IC 95%:                  0.919 – 0.966
cat(sprintf("F-test:                  p = %.4f\n", icc_res$p.value))
## F-test:                  p = 0.0000
cat("───────────────────────────────────────────────\n")
## ───────────────────────────────────────────────
interp_icc <- dplyr::case_when(
  icc_res$value >= 0.90 ~ "EXCELENTE",
  icc_res$value >= 0.75 ~ "BUENA",
  icc_res$value >= 0.50 ~ "MODERADA",
  TRUE                  ~ "POBRE"
)
cat(paste0("Interpretacion: ", interp_icc, "\n"))
## Interpretacion: EXCELENTE
cat("(Criterios de Cicchetti, 1994)\n")
## (Criterios de Cicchetti, 1994)
cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
# Gráfica de concordancia (Bland-Altman simplificado)
df_ba <- data.frame(
  media    = (df_tr$ocasion1 + df_tr$ocasion2) / 2,
  diferencia = df_tr$ocasion1 - df_tr$ocasion2
)
media_dif <- mean(df_ba$diferencia)
sd_dif    <- sd(df_ba$diferencia)
loa_sup   <- media_dif + 1.96 * sd_dif
loa_inf   <- media_dif - 1.96 * sd_dif

ggplot(df_ba, aes(x = media, y = diferencia)) +
  geom_point(alpha = 0.6, color = "#1B5E8E", size = 2.5) +
  geom_hline(yintercept = media_dif, color = "#2D8C6F",  linewidth = 1.2) +
  geom_hline(yintercept = loa_sup,   color = "#C0392B", linewidth = 1, linetype = "dashed") +
  geom_hline(yintercept = loa_inf,   color = "#C0392B", linewidth = 1, linetype = "dashed") +
  annotate("text", x = max(df_ba$media) * 0.85, y = media_dif + 0.8,
           label = sprintf("Media = %.2f", media_dif), color = "#2D8C6F") +
  annotate("text", x = max(df_ba$media) * 0.85, y = loa_sup + 0.8,
           label = sprintf("LOA+ = %.2f", loa_sup), color = "#C0392B") +
  annotate("text", x = max(df_ba$media) * 0.85, y = loa_inf - 0.8,
           label = sprintf("LOA− = %.2f", loa_inf), color = "#C0392B") +
  labs(title = "Gráfica de Bland-Altman: Concordancia Test-Retest",
       subtitle = "Diferencia vs. Media de las dos mediciones",
       x = "Media de ambas mediciones", y = "Diferencia (T1 − T2)") +
  theme_minimal(base_size = 13)

Interpretación: La gráfica de Bland-Altman muestra si las diferencias entre mediciones tienen una magnitud clínicamente aceptable y si hay tendencias sistemáticas. Las Líneas de Acuerdo (LOA) contienen el 95% de las diferencias. Si las LOA son clínicamente estrechas → buena reproducibilidad.


6.3 6.3 Curva ROC y AUC

# ── Desempeño diagnóstico de un instrumento de tamizaje ──────────────────────
set.seed(999)
n_roc <- 300
enfermo <- c(rep(1, 120), rep(0, 180))  # 120 casos, 180 sanos

# Puntaje del instrumento
puntaje <- ifelse(enfermo == 1,
                  rnorm(120, mean = 70, sd = 12),
                  rnorm(180, mean = 50, sd = 13))

df_roc <- data.frame(puntaje = puntaje, enfermo = enfermo)

# Curva ROC
roc_res <- roc(df_roc$enfermo, df_roc$puntaje, quiet = TRUE)

cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
cat("DESEMPEÑO DIAGNÓSTICO DEL INSTRUMENTO\n")
## DESEMPEÑO DIAGNÓSTICO DEL INSTRUMENTO
cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
cat(sprintf("AUC (Area bajo la curva ROC): %.3f\n", auc(roc_res)))
## AUC (Area bajo la curva ROC): 0.866
cat(sprintf("IC 95%% AUC:                    %.3f – %.3f\n",
            ci.auc(roc_res)[1], ci.auc(roc_res)[3]))
## IC 95% AUC:                    0.825 – 0.906
cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
# Gráfica ROC
plot(roc_res, col = "#1B5E8E", linewidth = 2,
     main = paste0("Curva ROC – AUC = ", round(auc(roc_res), 3)))
abline(0, 1, lty = 2, col = "#6B8BA4")
legend("bottomright",
       legend = c(sprintf("AUC = %.3f", auc(roc_res)), "Línea de no discriminación"),
       col = c("#1B5E8E","#6B8BA4"), lty = c(1,2), lwd = 2)

# Umbral óptimo por índice de Youden
coords_youden <- coords(roc_res, "best", best.method = "youden")
cat(sprintf("\nUmbral óptimo (Youden):   %.1f\n", coords_youden$threshold))
## 
## Umbral óptimo (Youden):   61.8
cat(sprintf("Sensibilidad:              %.1f%%\n", coords_youden$sensitivity * 100))
## Sensibilidad:              70.8%
cat(sprintf("Especificidad:             %.1f%%\n", coords_youden$specificity * 100))
## Especificidad:             84.4%

7 Módulo 7 – Metaanálisis

7.1 7.1 Metaanálisis de OR

# ── Datos: 8 estudios sobre tabaco y cáncer de pulmón ─────────────────────────
df_meta <- data.frame(
  estudio  = paste("Estudio", LETTERS[1:8]),
  n_casos  = c(200, 150, 300, 120, 400, 250, 180, 320),
  n_ctrl   = c(200, 150, 300, 120, 400, 250, 180, 320),
  exp_cas  = c(145, 98, 210, 78, 290, 170, 125, 220),   # expuestos en casos
  exp_ctrl = c(90,  55, 100, 40, 160,  90,  70, 110)    # expuestos en controles
) |>
  mutate(
    no_exp_cas  = n_casos - exp_cas,
    no_exp_ctrl = n_ctrl  - exp_ctrl,
    OR    = (exp_cas * no_exp_ctrl) / (no_exp_cas * exp_ctrl),
    logOR = log(OR),
    SE    = sqrt(1/exp_cas + 1/no_exp_cas + 1/exp_ctrl + 1/no_exp_ctrl)
  )

# Metaanálisis (efectos aleatorios)
meta_res <- metagen(
  TE   = df_meta$logOR,
  seTE = df_meta$SE,
  studlab = df_meta$estudio,
  sm   = "OR",
  random = TRUE,
  fixed  = TRUE,
  title = "Tabaco y Cáncer de Pulmón"
)

cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
cat("METAANÁLISIS: Tabaco y Cáncer de Pulmón\n")
## METAANÁLISIS: Tabaco y Cáncer de Pulmón
cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
cat(sprintf("OR combinado (Efectos aleatorios): %.2f\n", exp(meta_res$TE.random)))
## OR combinado (Efectos aleatorios): 3.87
cat(sprintf("IC 95%%:                             %.2f – %.2f\n",
            exp(meta_res$lower.random), exp(meta_res$upper.random)))
## IC 95%:                             3.38 – 4.43
cat(sprintf("I² (heterogeneidad):                %.1f%%\n", meta_res$I2 * 100))
## I² (heterogeneidad):                0.0%
cat(sprintf("Q-test:                             p = %.4f\n", meta_res$pval.Q))
## Q-test:                             p = 0.9015
cat("═══════════════════════════════════════════════\n")
## ═══════════════════════════════════════════════
if (meta_res$I2 < 0.25) {
  cat("Heterogeneidad BAJA → Efectos fijos es apropiado\n")
} else if (meta_res$I2 < 0.50) {
  cat("Heterogeneidad MODERADA → Investigar fuente; considerar efectos aleatorios\n")
} else {
  cat("Heterogeneidad ALTA → Efectos aleatorios; investigar metarregresión\n")
}
## Heterogeneidad BAJA → Efectos fijos es apropiado
# Forest plot
forest(meta_res,
       sortvar    = exp(meta_res$TE),
       xlim       = c(0.5, 12),
       leftcols   = c("studlab"),
       leftlabs   = c("Estudio"),
       smlab      = "Odds Ratio",
       col.square = "#1B5E8E",
       col.diamond= "#C0392B",
       fontsize   = 10)


8 Preguntas de Interpretación Final

Responda en grupo las siguientes preguntas basándose en los resultados de este taller:

  1. Incidencia: En el Módulo 1.2, ¿por qué la DI resultó diferente a la IA? ¿Cuál de las dos usaría si tuviera datos de una cohorte hospitalaria con ingresos continuos?

Respuesta:

✅ Usar DI porque:

Nuevos pacientes ingresan continuamente → tiempos de seguimiento heterogéneos Algunos se dan de alta, otros mueren, otros se pierden → datos censurados La fórmula DI = Casos / Σ(personas-año) incorpora correctamente la variabilidad temporal La IA sobreestimaría o subestimaría según quién abandone antes

  1. OR vs. RR: En el Módulo 2.2, si la prevalencia de HTA fuera del 40%, ¿seguiría siendo apropiado reportar el OR de la regresión logística? ¿Qué modelo usaría?

  2. Sesgo: En el Módulo 3.1, la curva muestra que con S=70% el OR verdadero de 4.0 aparece como 2.0. ¿Qué concluiría si el estudio reportara OR=2.0 sin mencionar la validación del instrumento?

  3. Cox: En el Módulo 4.2, ¿qué haría si la prueba de Schoenfeld resultara p=0.03 para la variable de tratamiento?

  4. Transversal vs. Cohorte: En el Módulo 5.1, la RP ajustada (Poisson) y el OR (logística) dieron valores distintos. ¿Cuál reportaría y por qué?

  5. Validación: En el Módulo 6.1, el Alpha fue muy alto (>0.90). ¿Qué pasos adicionales realizaría antes de usar el instrumento en investigación?

  6. Metaanálisis: En el Módulo 7, si I²=65%, ¿utilizaría el OR de efectos fijos o de efectos aleatorios para hacer las recomendaciones de política de salud? Justifique.


9 Información de la Sesión

sessionInfo()
## R version 4.5.1 (2025-06-13 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22631)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8  LC_CTYPE=Spanish_Colombia.utf8   
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.utf8    
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] pROC_1.19.0.1     scales_1.4.0      broom_1.0.10      psych_2.6.3      
##  [5] irr_0.84.1        lpSolve_5.6.23    meta_8.3-0        metadat_1.6-0    
##  [9] metabook_0.2-0    kableExtra_1.4.0  knitr_1.50        tidyr_1.3.1      
## [13] dplyr_1.1.4       survminer_0.5.2   ggpubr_0.6.2      ggplot2_4.0.0    
## [17] epiR_2.0.92       survival_3.8-3    epitools_0.5-10.1
## 
## loaded via a namespace (and not attached):
##  [1] mnormt_2.1.2            Rdpack_2.6.4            DBI_1.2.3              
##  [4] gridExtra_2.3           BiasedUrn_2.0.12        rlang_1.1.6            
##  [7] magrittr_2.0.4          e1071_1.7-17            compiler_4.5.1         
## [10] systemfonts_1.3.1       vctrs_0.6.5             stringr_1.5.2          
## [13] pkgconfig_2.0.3         fastmap_1.2.0           backports_1.5.0        
## [16] labeling_0.4.3          pander_0.6.6            rmarkdown_2.30         
## [19] tzdb_0.5.0              nloptr_2.2.1            ragg_1.5.0             
## [22] purrr_1.2.0             xfun_0.53               cachem_1.1.0           
## [25] jsonlite_2.0.0          uuid_1.2-1              parallel_4.5.1         
## [28] R6_2.6.1                bslib_0.9.0             stringi_1.8.7          
## [31] RColorBrewer_1.1-3      car_3.1-3               boot_1.3-31            
## [34] lubridate_1.9.4         jquerylib_0.1.4         numDeriv_2016.8-1.1    
## [37] Rcpp_1.1.0              zoo_1.8-14              readr_2.1.5            
## [40] Matrix_1.7-3            splines_4.5.1           timechange_0.3.0       
## [43] tidyselect_1.2.1        rstudioapi_0.17.1       abind_1.4-8            
## [46] yaml_2.3.10             ggtext_0.1.2            metafor_5.0-1          
## [49] lattice_0.22-7          tibble_3.3.0            withr_3.0.2            
## [52] S7_0.2.0                flextable_0.9.10        askpass_1.2.1          
## [55] evaluate_1.0.5          sf_1.1-1                CompQuadForm_1.4.4     
## [58] units_1.0-1             proxy_0.4-29            zip_2.3.3              
## [61] xml2_1.5.1              pillar_1.11.1           carData_3.0-5          
## [64] KernSmooth_2.23-26      reformulas_0.4.2        generics_0.1.4         
## [67] mathjaxr_2.0-0          hms_1.1.4               minqa_1.2.8            
## [70] class_7.3-23            glue_1.8.0              gdtools_0.4.4          
## [73] tools_4.5.1             data.table_1.17.8       lme4_1.1-38            
## [76] ggsignif_0.6.4          grid_4.5.1              rbibutils_2.4          
## [79] nlme_3.1-168            Formula_1.2-5           cli_3.6.5              
## [82] textshaping_1.0.4       officer_0.7.2           fontBitstreamVera_0.1.1
## [85] viridisLite_0.4.2       svglite_2.2.2           gtable_0.3.6           
## [88] rstatix_0.7.3           sass_0.4.10             digest_0.6.37          
## [91] fontquiver_0.2.1        classInt_0.4-11         farver_2.1.2           
## [94] htmltools_0.5.8.1       lifecycle_1.0.4         gridtext_0.1.6         
## [97] fontLiberation_0.1.0    openssl_2.3.4           MASS_7.3-65

Taller elaborado para el Curso de Diseño y Análisis de Datos – Investigación Epidemiológica
Fuente: Rothman KJ & Greenland S. Modern Epidemiology, 3rd ed. (2008)