Taller Práctico de Epidemiología en R
Medidas de Frecuencia, Asociación, Sesgos y Diseño de Estudios
José William Martínez | jose01william@utp.edu.co; Estudiantes: Isabel Cristina Peña | Martha Raquel Niño Chacon | Lina Jhojana Roa | Betty Perez | Jenny Santana
mayo 09, 2026
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")## ═══════════════════════════════════════════════## Personas en la cohorte: 500## Casos nuevos: 61## Tiempo-persona total: 2579.4 personas-año## ───────────────────────────────────────────────## Incidencia Acumulada (IA): 122.0 por 1,000## Densidad de Incidencia (DI): 23.6 por 1,000 personas-año## ═══════════════════════════════════════════════##
## ⚠️ Nota: IA supone seguimiento uniforme. DI es más apropiada## 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")## ═══════════════════════════════════════════════## TABLA 2×2## ═══════════════════════════════════════════════## Desenlace
## Exposición Cáncer (+) Sano (−)
## Fumadores 72 728
## No fumadores 24 1176## ═══════════════════════════════════════════════# ── 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## IA en no fumadores: 20.0 por 1,000## ───────────────────────────────────────────────## Riesgo Relativo (RR): 4.50## Odds Ratio (OR): 4.85## Riesgo Atribuible: 70.0 por 1,000## %RA en expuestos: 77.8%## %RA poblacional: 58.3%## ═══════════════════════════════════════════════##
## 💡 Interpretación del RR = 4.5 :## Los fumadores tienen 4.5 veces más riesgo de## cáncer de pulmón que los no fumadores.## 💡 Interpretación del %RA = 77.8 %:## El 77.8 % del riesgo de cáncer en fumadores## 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")## ═══════════════════════════════════════════════## RIESGO RELATIVO (epitools)## risk ratio with 95% C.I.
## Exposición estimate lower upper
## Fumadores 1.000000 NA NA
## No fumadores 1.076923 1.052181 1.102247##
## ODDS RATIO (epitools)## odds ratio with 95% C.I.
## Exposición estimate lower upper
## Fumadores 1.000000 NA NA
## No fumadores 4.846154 3.025455 7.762537## ═══════════════════════════════════════════════# ── 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")
##
## Prueba de homogeneidad de Mantel-Haenszel:##
## 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.638889Interpretació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)| 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 |
##
## 📌 Verificación del supuesto de riesgos proporcionales (Schoenfeld):## 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.362if (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")## ═══════════════════════════════════════════════## RESULTADOS: ESTUDIO TRANSVERSAL## ═══════════════════════════════════════════════## Prevalencia HTA en obesos: 59.0%## Prevalencia HTA en no obesos: 36.1%## Razón de Prevalencias (RP) cruda: 1.63## ═══════════════════════════════════════════════# 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## OR ajustado (Logística): 2.60##
## Prevalencia de HTA en la muestra: 44.4%## → Con prevalencia alta, OR sobreestima considerablemente la RP.## → 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## 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)| 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")## ═══════════════════════════════════════════════## ANÁLISIS DE CONSISTENCIA INTERNA## ═══════════════════════════════════════════════## Alpha de Cronbach: 0.905## Omega de McDonald: 0.923## ───────────────────────────────────────────────## 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## ═══════════════════════════════════════════════##
## 💡 Alpha = 0.905if (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")## ═══════════════════════════════════════════════## CONFIABILIDAD TEST-RETEST (ICC)## ═══════════════════════════════════════════════## ICC (acuerdo absoluto): 0.947## IC 95%: 0.919 – 0.966## F-test: p = 0.0000## ───────────────────────────────────────────────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## (Criterios de Cicchetti, 1994)## ═══════════════════════════════════════════════# 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")## ═══════════════════════════════════════════════## DESEMPEÑO DIAGNÓSTICO DEL INSTRUMENTO## ═══════════════════════════════════════════════## AUC (Area bajo la curva ROC): 0.866## IC 95% AUC: 0.825 – 0.906## ═══════════════════════════════════════════════# 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## Sensibilidad: 70.8%## 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")## ═══════════════════════════════════════════════## METAANÁLISIS: Tabaco y Cáncer de Pulmón## ═══════════════════════════════════════════════## OR combinado (Efectos aleatorios): 3.87## IC 95%: 3.38 – 4.43## I² (heterogeneidad): 0.0%## Q-test: p = 0.9015## ═══════════════════════════════════════════════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:
- 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?
La diferencia se da porque la incidencia acumulada mide el riesgo sin tener en cuenta el momento en que ocurre el evento, mientras que la densidad de incidencia mide el tiempo que la persona estuvo expuesta. Para el caso de una cohorte hospitalaria con ingresos continuos, seria mas apropiado usar la densidad de incidencia porque conocemos el tiempo de exposicion para cada individuo.
- 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?
Cuando la prevalencia es del 40% no seria apropiado usar el OR, debido a que sobreestima sustancialmente el RR. se deberia usar el modelo de regresión de Poisson para obtener la razón de prevalencias directamente.
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? Se concluiria que existe un sesgo relacionado con un error de clasificación no diferencial.
Cox: En el Módulo 4.2, ¿qué haría si la prueba de Schoenfeld resultara p=0.03 para la variable de tratamiento?
Se rechaza de manera estadisticamente significativa ya que se viola el supuesto de riesgos proporcionales para esa variable por ser la p < 0.05, y por tanto se tendria que suspender la exposicion al tratamiento, porque el efecto no seria constante a lo largo del tiempo.
- 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é?
Se reportaria la Razón de Prevalencias (RP) ajustada porque es un estudio transversal y el desenlace tiene una prevalencia alta: 44.4%, por tanto la OR sobreestima la RP
- 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?
Antes de usar el instrumento en investigación se deberia revisar y verificar que no exista redundancia entre los ítems y unificar los que se repitan.
- 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. Un I² de 65% nos da una heterogeneidad alta, por tanto se deberia usar el OR de efectos aleatorios para hacer las recomendaciones de política de salud.
9 Información de la Sesión
## R version 4.6.0 (2026-04-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=Spanish_Spain.utf8 LC_CTYPE=Spanish_Spain.utf8
## [3] LC_MONETARY=Spanish_Spain.utf8 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Spain.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.12 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.51 tidyr_1.3.2
## [13] dplyr_1.2.1 survminer_0.5.2 ggpubr_0.6.3 ggplot2_4.0.3
## [17] epiR_2.0.92 survival_3.8-6 epitools_0.5-10.1
##
## loaded via a namespace (and not attached):
## [1] mnormt_2.1.2 Rdpack_2.6.6 DBI_1.3.0
## [4] gridExtra_2.3 BiasedUrn_2.0.12 rlang_1.2.0
## [7] magrittr_2.0.5 e1071_1.7-17 compiler_4.6.0
## [10] systemfonts_1.3.2 vctrs_0.7.3 stringr_1.6.0
## [13] pkgconfig_2.0.3 fastmap_1.2.0 backports_1.5.1
## [16] labeling_0.4.3 pander_0.6.6 rmarkdown_2.31
## [19] tzdb_0.5.0 nloptr_2.2.1 ragg_1.5.2
## [22] purrr_1.2.2 xfun_0.57 cachem_1.1.0
## [25] jsonlite_2.0.0 uuid_1.2-2 parallel_4.6.0
## [28] R6_2.6.1 bslib_0.10.0 stringi_1.8.7
## [31] RColorBrewer_1.1-3 car_3.1-5 boot_1.3-32
## [34] lubridate_1.9.5 jquerylib_0.1.4 numDeriv_2016.8-1.1
## [37] Rcpp_1.1.1-1.1 zoo_1.8-15 readr_2.2.0
## [40] Matrix_1.7-5 splines_4.6.0 timechange_0.4.0
## [43] tidyselect_1.2.1 rstudioapi_0.18.0 abind_1.4-8
## [46] yaml_2.3.12 ggtext_0.1.2 metafor_5.0-1
## [49] lattice_0.22-9 tibble_3.3.1 withr_3.0.2
## [52] S7_0.2.2 flextable_0.9.11 askpass_1.2.1
## [55] evaluate_1.0.5 sf_1.1-0 CompQuadForm_1.4.4
## [58] units_1.0-1 proxy_0.4-29 zip_2.3.3
## [61] xml2_1.5.2 pillar_1.11.1 carData_3.0-6
## [64] KernSmooth_2.23-26 reformulas_0.4.4 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.1 gdtools_0.5.0
## [73] tools_4.6.0 data.table_1.18.2.1 lme4_2.0-1
## [76] ggsignif_0.6.4 grid_4.6.0 rbibutils_2.4.1
## [79] nlme_3.1-169 Formula_1.2-5 cli_3.6.6
## [82] textshaping_1.0.5 officer_0.7.4 fontBitstreamVera_0.1.1
## [85] viridisLite_0.4.3 svglite_2.2.2 gtable_0.3.6
## [88] rstatix_0.7.3 sass_0.4.10 digest_0.6.39
## [91] fontquiver_0.2.1 classInt_0.4-11 farver_2.1.2
## [94] htmltools_0.5.9 lifecycle_1.0.5 gridtext_0.1.6
## [97] fontLiberation_0.1.0 openssl_2.4.0 MASS_7.3-65Taller 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)