Al finalizar este capítulo, serás capaz de:
En estudios observacionales, nunca podemos estar seguros de haber controlado todos los confusores. El análisis de sensibilidad evalúa qué tan fuertes tendrían que ser los confusores no medidos para explicar nuestros resultados.
library(ggdag)
library(ggplot2)
sens_dag <- dagify(
Y ~ X + C + U,
X ~ C + U,
coords = list(
x = c(X = 0, Y = 2, C = 1, U = 1),
y = c(X = 0, Y = 0, C = 0.7, U = -0.7)
),
labels = c(
X = "Exposición",
Y = "Resultado",
C = "Confusor\n(medido)",
U = "Confusor\n(no medido)"
)
)
ggdag(sens_dag, text = FALSE, use_labels = "label") +
theme_dag() +
labs(title = "U podría explicar la asociación observada")
El E-value (VanderWeele y Ding 2017) es la asociación mínima que un confusor no medido tendría que tener con tanto la exposición como el resultado para explicar completamente la asociación observada.
\[E\text{-value} = RR + \sqrt{RR \times (RR - 1)}\]
library(EValue)
# Ejemplo: RR observado = 2.5
rr_observado <- 2.5
ic_inferior <- 1.8 # Límite inferior del IC 95%
# Calcular E-value
evalues <- evalues.RR(est = rr_observado, lo = ic_inferior)
evalues point lower upper
RR 2.500000 1.8 NA
E-values 4.436492 3.0 NA# Crear datos para el gráfico
rr_seq <- seq(1.1, 5, 0.1)
evalue_seq <- rr_seq + sqrt(rr_seq * (rr_seq - 1))
plot_data <- data.frame(
RR = rr_seq,
EValue = evalue_seq
)
# Valores del ejemplo anterior
rr_ejemplo <- 2.5
evalue_ejemplo <- 2.5 + sqrt(2.5 * (2.5 - 1))
ggplot(plot_data, aes(x = RR, y = EValue)) +
geom_line(size = 1.2, color = "steelblue") +
geom_vline(xintercept = rr_ejemplo, linetype = "dashed", color = "red") +
geom_hline(yintercept = evalue_ejemplo, linetype = "dashed", color = "red") +
annotate("point", x = rr_ejemplo, y = evalue_ejemplo, size = 4, color = "red") +
labs(x = "Riesgo Relativo observado",
y = "E-value",
title = "E-value en función del RR observado") +
theme_minimal()
Para un RR de 2.5, el E-value es aproximadamente 4.4. Esto significa que un confusor no medido tendría que tener un RR de al menos 4.4 con tanto la exposición como el resultado para explicar completamente la asociación observada.
El paquete sensemakr (Cinelli y Hazlett 2020) proporciona herramientas más sofisticadas.
library(sensemakr)
# Ejemplo con datos de discriminación laboral
data("darfur")
# Modelo ajustado
modelo <- lm(peacefactor ~ directlyharmed + age + farmer_dar +
herder_dar + pastvoted + female, data = darfur)
# Análisis de sensibilidad
sens <- sensemakr(model = modelo,
treatment = "directlyharmed",
benchmark_covariates = "female",
kd = 1:3,
q = 1)
summary(sens)Sensitivity Analysis to Unobserved Confounding
Model Formula: peacefactor ~ directlyharmed + age + farmer_dar + herder_dar +
pastvoted + female
Null hypothesis: q = 1 and reduce = TRUE
-- This means we are considering biases that reduce the absolute value of the current estimate.
-- The null hypothesis deemed problematic is H0:tau = 0
Unadjusted Estimates of 'directlyharmed':
Coef. estimate: 0.0489
Standard Error: 0.0184
t-value (H0:tau = 0): 2.6648
Sensitivity Statistics:
Partial R2 of treatment with outcome: 0.0056
Robustness Value, q = 1: 0.0721
Robustness Value, q = 1, alpha = 0.05: 0.0195
Verbal interpretation of sensitivity statistics:
-- Partial R2 of the treatment with the outcome: an extreme confounder (orthogonal to the covariates) that explains 100% of the residual variance of the outcome, would need to explain at least 0.56% of the residual variance of the treatment to fully account for the observed estimated effect.
-- Robustness Value, q = 1: unobserved confounders (orthogonal to the covariates) that explain more than 7.21% of the residual variance of both the treatment and the outcome are strong enough to bring the point estimate to 0 (a bias of 100% of the original estimate). Conversely, unobserved confounders that do not explain more than 7.21% of the residual variance of both the treatment and the outcome are not strong enough to bring the point estimate to 0.
-- Robustness Value, q = 1, alpha = 0.05: unobserved confounders (orthogonal to the covariates) that explain more than 1.95% of the residual variance of both the treatment and the outcome are strong enough to bring the estimate to a range where it is no longer 'statistically different' from 0 (a bias of 100% of the original estimate), at the significance level of alpha = 0.05. Conversely, unobserved confounders that do not explain more than 1.95% of the residual variance of both the treatment and the outcome are not strong enough to bring the estimate to a range where it is no longer 'statistically different' from 0, at the significance level of alpha = 0.05.
Bounds on omitted variable bias:
--The table below shows the maximum strength of unobserved confounders with association with the treatment and the outcome bounded by a multiple of the observed explanatory power of the chosen benchmark covariate(s).
Bound Label R2dz.x R2yz.dx Treatment Adjusted Estimate Adjusted Se
1x female 0.0027 0.1314 directlyharmed 0.0366 0.0171
2x female 0.0054 0.2627 directlyharmed 0.0243 0.0158
3x female 0.0081 0.3941 directlyharmed 0.0119 0.0144
Adjusted T Adjusted Lower CI Adjusted Upper CI
2.1357 0.0030 0.0703
1.5341 -0.0068 0.0553
0.8273 -0.0163 0.0401plot(sens)
# El Robustness Value indica qué tan fuerte tendría que ser
# un confusor para reducir el efecto a cero
cat("Robustness Value (q=1):",
round(sens$bounds$r2yz.dx[1], 3), "\n")Robustness Value (q=1): 0.131 cat("\nInterpretación: Un confusor tendría que explicar al menos",
round(sens$bounds$r2yz.dx[1] * 100, 1), "% de la varianza\n",
"residual tanto de X como de Y para eliminar el efecto.\n")
Interpretación: Un confusor tendría que explicar al menos 13.1 % de la varianza
residual tanto de X como de Y para eliminar el efecto.La fórmula de sesgo permite calcular cuánto cambiaría la estimación dado un confusor específico:
\[\text{Sesgo} = \frac{(\text{RR}_{UY} - 1) \times (\text{RR}_{UX} - 1)}{\text{RR}_{UY} + (\text{RR}_{UX} - 1) \times P(U)}\]
# Función para calcular sesgo
calcular_sesgo <- function(rr_uy, rr_ux, p_u = 0.5) {
# Fórmula simplificada
bias <- (rr_uy - 1) * (rr_ux - 1) /
(rr_uy + (rr_ux - 1) * p_u)
return(bias)
}
# Ejemplo: ¿Cuánto sesgo introduciría un confusor con RR=1.5
# con exposición y resultado?
sesgo <- calcular_sesgo(rr_uy = 1.5, rr_ux = 1.5)
cat("Sesgo aproximado:", round(sesgo, 3), "\n")Sesgo aproximado: 0.143 # Simulación con confusor no medido
set.seed(606)
n <- 2000
# Confusor medido
C_medido <- rnorm(n)
# Confusor NO medido
U <- rnorm(n)
# Exposición
prob_X <- plogis(-1 + 0.5 * C_medido + 0.8 * U)
X <- rbinom(n, 1, prob_X)
# Resultado
Y <- 2 + 3 * X + 1.5 * C_medido + 2 * U + rnorm(n)
datos_sens <- data.frame(C_medido, U, X, Y)
# Efecto verdadero (controlando U)
efecto_verdadero <- coef(lm(Y ~ X + C_medido + U, data = datos_sens))["X"]
# Efecto sin U (sesgado)
efecto_sesgado <- coef(lm(Y ~ X + C_medido, data = datos_sens))["X"]
# IPTW sin U
library(WeightIt)
pesos <- weightit(X ~ C_medido, data = datos_sens, method = "ps")
efecto_iptw <- coef(lm(Y ~ X, data = datos_sens, weights = pesos$weights))["X"]
cat("Efecto verdadero (con U):", round(efecto_verdadero, 3), "\n")Efecto verdadero (con U): 2.958 cat("Efecto sin U (OLS):", round(efecto_sesgado, 3), "\n")Efecto sin U (OLS): 4.468 cat("Efecto IPTW (sin U):", round(efecto_iptw, 3), "\n")Efecto IPTW (sin U): 4.432 cat("Sesgo por U:", round(efecto_iptw - efecto_verdadero, 3), "\n")Sesgo por U: 1.473 # Supongamos un RR observado
rr_obs <- 1.8
rr_lo <- 1.3
rr_hi <- 2.4
# E-values
ev <- evalues.RR(est = rr_obs, lo = rr_lo, hi = rr_hi)
cat("REPORTE DE SENSIBILIDAD\n")REPORTE DE SENSIBILIDADcat("=======================\n\n")=======================cat("RR observado: ", rr_obs, " (IC 95%: ", rr_lo, "-", rr_hi, ")\n\n", sep = "")RR observado: 1.8 (IC 95%: 1.3-2.4)cat("E-value para el estimado puntual:", round(ev["E-values", "point"], 2), "\n")E-value para el estimado puntual: 3 cat("E-value para el límite inferior del IC:", round(ev["E-values", "lower"], 2), "\n\n")E-value para el límite inferior del IC: 1.92 cat("Interpretación: Para explicar completamente la asociación observada,\n")Interpretación: Para explicar completamente la asociación observada,cat("un confusor no medido tendría que estar asociado con tanto la\n")un confusor no medido tendría que estar asociado con tanto lacat("exposición como el resultado con un RR de al menos",
round(ev["E-values", "point"], 2), ".\n\n")exposición como el resultado con un RR de al menos 3 .cat("Para mover el IC inferior a 1.0, la asociación del confusor\n")Para mover el IC inferior a 1.0, la asociación del confusorcat("tendría que ser de al menos RR =", round(ev["E-values", "lower"], 2), ".\n")tendría que ser de al menos RR = 1.92 .Un estudio observacional encuentra que el consumo de vegetales está asociado con menor mortalidad (HR = 0.75, IC 95%: 0.65-0.85).
Usando sensemakr, realiza un análisis de sensibilidad completo para un modelo de tu elección. Interpreta los resultados y determina si el efecto es robusto.
sensemakr proporciona visualizaciones y métricas avanzadas