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.
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.
# ── 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?
# ── 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.
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.
# ── 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.
# ── 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.
# ── 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.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.
# ── 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))# ── 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.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.
# ── 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.
# ── 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)# ── 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.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")# ── 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.
# ── 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%
# ── 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)Responda en grupo las siguientes preguntas basándose en los resultados de este taller:
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
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?
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?
Cox: En el Módulo 4.2, ¿qué haría si la prueba de Schoenfeld resultara p=0.03 para la variable de tratamiento?
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é?
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?
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.
## 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)