Mostrar código
library(simstudy)
library(data.table)
library(ggplot2)
library(future)
library(furrr)
library(parallelly)Métodos doblemente robustos en el contexto de estimación de máxima verosimilitud con validación cruzada y aprendizaje automático: una explicación matemática de TMLE
Una serie en tres partes sobre funciones de influencia, ortogonalidad y el paso de fluctuación en TMLE.
Abstract
Este documento recopila tres entradas del blog ouR data generation (Keith Goldfeld, NYU) en una referencia académica unificada sobre los fundamentos matemáticos de la Estimación por Máxima Verosimilitud Dirigida (TMLE, por sus siglas en inglés). La Parte I introduce cuatro pilares conceptuales: los parámetros como funcionales de distribuciones de probabilidad, el muestreo como perturbación de la distribución generadora de datos, la linealización del error de estimación mediante funciones de influencia y la ortogonalidad respecto al error de estimación de parámetros de perturbación. La Parte II examina la propiedad de ortogonalidad mediante simulación, demostrando que el término de interacción impulsado por los parámetros de perturbación efectivamente se reduce a cero conforme aumenta el tamaño muestral, incluso bajo una mala especificación deliberada del modelo de resultados, mientras permanece centrado en cero a lo largo de todo el proceso. La Parte III deriva el paso de dirección (fluctuación) de TMLE para el efecto promedio del tratamiento, mostrando cómo una actualización unidimensional de la regresión de resultados inicial hace cumplir la ecuación empírica de la función de influencia eficiente (EIF), restaura la condición de centrado alterada por el error de estimación de los parámetros de perturbación y da lugar a la doble robustez. En conjunto, las tres partes ofrecen una explicación accesible pero matemáticamente fundamentada de por qué funciona TMLE, no meramente de cómo ejecutarlo.
Keywords
estimación por máxima verosimilitud dirigida, inferencia causal, funciones de influencia, teoría semiparamétrica, efecto promedio del tratamiento, estimación doblemente robusta, validación cruzada
1 Parte I. Funciones de Influencia y Perturbaciones
1.1 ¿Qué problema intenta resolver TMLE?
La Estimación por Máxima Verosimilitud Dirigida (TMLE), a veces denominada estimación por mínima pérdida dirigida, aborda una tensión que recorre la inferencia causal moderna y la estadística semiparamétrica. Por un lado, los parámetros de interés científico dependen con frecuencia de características complejas del proceso generador de datos: los efectos causales, por ejemplo, requieren modelar conjuntamente tanto los resultados como los mecanismos de asignación al tratamiento. Por otro lado, los componentes necesarios para identificar esos parámetros se estiman cada vez más con métodos flexibles de aprendizaje automático que necesitan muestras grandes para estabilizarse y cuyas propiedades estadísticas resisten el tratamiento analítico clásico.
La inferencia estadística clásica descansa en modelos estructurados en los que el sesgo, la varianza y la incertidumbre se derivan directamente de los supuestos del modelo. Los estimadores flexibles de aprendizaje automático violan muchos de esos supuestos: son muy adaptativos, implican numerosas decisiones de ajuste y, en general, no corresponden a funciones de verosimilitud paramétricas suaves. Aunque tales métodos pueden ofrecer un excelente rendimiento predictivo, sus propiedades inferenciales son difíciles de caracterizar analíticamente.
TMLE resuelve esta tensión desacoplando el objetivo científico de su tratamiento estadístico. En lugar de construir estimadores directamente a partir de modelos de resultados o de tratamiento—que, aunque científicamente importantes, se tratan como componentes de perturbación—TMLE construye estimadores cuyo comportamiento asintótico está gobernado por una función de influencia conocida. Una vez que un estimador posee la función de influencia correcta, su sesgo, varianza e incertidumbre a gran muestra pueden caracterizarse mediante teoría asintótica general, incluso cuando los componentes de perturbación se estiman con métodos flexibles de convergencia más lenta.
El resto de la Parte I desarrolla cuatro ideas que, en conjunto, permiten comprender TMLE:
- ver los parámetros como funcionales de distribuciones de probabilidad;
- interpretar el muestreo como una perturbación de la distribución verdadera;
- linealizar el error de estimación mediante la función de influencia; y
- garantizar la ortogonalidad respecto al error de estimación de los parámetros de perturbación en la estimación de los modelos para el tratamiento y el resultado desde la perspectiva de la teoría de la estimación semiparamétrica.
1.2 Estimadores como funcionales
Un estimador, como la media muestral o un efecto causal, no es simplemente una fórmula aplicada a los datos. Más precisamente, es un funcional: una aplicación que toma una distribución de probabilidad \(P\) de una familia \(\mathcal{P}\) y devuelve un número real \(T(P)\):
\[ T : \mathcal{P} \to \mathbb{R}. \]
Escribir el parámetro de esta manera pone de manifiesto que es una propiedad de la distribución generadora de datos, no de ningún conjunto de datos en particular. Los datos entran únicamente a través de la distribución empírica \(P_n\).
1.2.1 La media muestral como ejemplo canónico
La media muestral ofrece la ilustración más sencilla:
\[ T(P) = \int z\; dP(z). \]
Sea \(P_0\) la distribución generadora de datos verdadera, de modo que el parámetro objetivo es \(T(P_0)\). La distribución empírica es
\[ P_n = \frac{1}{n}\sum_{i=1}^{n}\delta_{Z_i}, \]
donde \(\delta_z\) es la masa puntual en \(z\)—una distribución que concentra toda la probabilidad en el valor único \(z\). Escribir \(P_n\) como un promedio de masas puntuales formaliza el hecho de que los datos observados asignan igual peso a cada punto muestreado; \(P_n\) es por tanto una aproximación discreta a \(P_0\) que coloca masa \(1/n\) en cada observación. El estimador de sustitución es simplemente \(T(P_n)\).
1.2.2 El error de estimación como diferencia funcional
Desde la perspectiva funcional, el error de estimación es la diferencia entre evaluar el mismo funcional en dos distribuciones cercanas:
\[ \hat{\theta} - \theta_0 = T(P_n) - T(P_0). \]
Comprender este error se reduce a entender cómo cambia \(T(P)\) cuando \(P\) pasa de \(P_0\) a su aproximación empírica \(P_n\). Ésta es la pregunta central que motiva las herramientas desarrolladas a continuación.
1.3 El muestreo como perturbación
1.3.1 La medida con signo \(P_n - P_0\)
Cada observación contribuye con una masa puntual de tamaño \(1/n\) a \(P_n\). La diferencia \(P_n - P_0\) no es en sí misma una distribución de probabilidad; es una medida con signo que registra dónde la distribución empírica sobre- o infrarepresenta la verdad. Puede concebirse como un balance: las regiones donde \(P_n\) coloca más masa que \(P_0\) aparecen como entradas positivas, las regiones donde coloca menos masa como entradas negativas, y el total siempre suma cero.
El objeto \(P_n - P_0\) solo cobra significado a través de su acción sobre funciones. Para cualquier función medible \(f\),
\[ (P_n - P_0)f \;\equiv\; \int f(z)\,dP_n(z) - \int f(z)\,dP_0(z), \]
lo que mide cuánto se desvía el promedio empírico de \(f(Z)\) de su media poblacional. Esta notación, aunque compacta, reaparecerá a lo largo del texto; toda cantidad clave de la teoría se expresa como la medida con signo actuando sobre alguna función.
1.3.2 El modelo de contaminación de Hampel
Para caracterizar la sensibilidad de \(T(P)\) a cambios pequeños en \(P\), se necesita una noción precisa de cambio infinitesimal en una distribución, que preserve la masa total de probabilidad y permita la diferenciación. El modelo de contaminación introducido por Hampel (1974) proporciona exactamente esto.
Dada la distribución verdadera \(P_0\) y la masa puntual \(\delta_z\), se define la familia uniparamétrica de distribuciones perturbadas
\[ P_{\epsilon} = (1-\epsilon)P_0 + \epsilon\,\delta_z, \qquad 0 < \epsilon < 1. \]
Esto no equivale a añadir un nuevo punto de datos completo; en cambio, introduce una cantidad infinitesimal de masa de probabilidad en \(z\), trazando un camino suave en el espacio de distribuciones a lo largo del cual se pueden calcular derivadas (pathwise first derivative: derivada direccional).
1.4 Linealización mediante la función de influencia
1.4.1 Definición
La función de influencia de \(T\) en \(P_0\), denotada \(\phi_{P_0}\), se define como la derivada direccional de \(T\) a lo largo de la perturbación de Hampel:
\[ \phi_{P_0}(z) \;\equiv\; \left.\frac{d}{d\epsilon}\,T(P_{\epsilon})\right|_{\epsilon=0} = \left.\frac{d}{d\epsilon}\,T\!\bigl((1-\epsilon)P_0 + \epsilon\,\delta_z\bigr)\right|_{\epsilon=0}. \]
En palabras, la función de influencia mide el efecto de primer orden sobre \(T(P)\) de perturbar infinitesimalmente la distribución en el punto \(z\). Es el análogo funcional de una derivada ordinaria.
1.4.2 Expansión de Taylor funcional
Si \(T\) es suficientemente regular, el error de estimación admite una expansión de primer orden:
\[ T(P_n) - T(P_0) = (P_n - P_0)\phi_{P_0} + R_n, \]
donde el término dominante es lineal en la perturbación \(P_n - P_0\) y \(R_n\) recoge los términos de resto de orden superior. Aplicando la definición de la medida con signo:
\[ (P_n - P_0)\phi_{P_0} = \frac{1}{n}\sum_{i=1}^{n}\phi_{P_0}(Z_i) - \underbrace{E_{P_0}[\phi_{P_0}(Z)]}_{=\;0}, \]
ya que las funciones de influencia se construyen para satisfacer la condición de centrado
\[ E_{P_0}[\phi_{P_0}(Z)] = 0. \]
La expansión debe leerse como una serie de Taylor funcional: \(\phi_{P_0}\) desempeña el papel de la derivada, y \((P_n - P_0)\phi_{P_0}\) es la aproximación lineal al cambio en \(T\) inducido al reemplazar \(P_0\) por \(P_n\).
1.4.3 Normalidad asintótica
Si el resto satisface \(R_n = o_p(n^{-1/2})\), multiplicando ambos lados por \(\sqrt{n}\) se obtiene
\[ \sqrt{n}\bigl(T(P_n) - T(P_0)\bigr) = \frac{1}{\sqrt{n}}\sum_{i=1}^{n}\phi_{P_0}(Z_i) + o_p(1). \]
La escala \(\sqrt{n}\) refleja el hecho conocido de que promediar \(n\) observaciones independientes reduce la desviación típica a una tasa \(n^{-1/2}\), de modo que multiplicar por \(\sqrt{n}\) mantiene visible la fluctuación aleatoria en el límite. Como \(\phi_{P_0}\) tiene media cero, el Teorema Central del Límite implica convergencia en distribución a una variable aleatoria normal centrada, lo que proporciona el marco estándar de inferencia asintótica (intervalos de confianza tipo Wald). TMLE está diseñado explícitamente para que su error de estimación admita esta expansión lineal escalada en \(\sqrt{n}\) incluso cuando los componentes de perturbación se estiman de forma flexible.
1.5 Ortogonalidad bajo estimación de los parámetros de perturbación
1.5.1 El problema
En la práctica, \(\phi_{P_0}\) nunca se conoce. Para parámetros causales, la función de influencia depende de funciones de perturbación —modelos condicionales de resultados y mecanismos de asignación al tratamiento— que deben estimarse a partir de los datos. Al reemplazar la función de influencia verdadera por una estimada \(\phi_{\hat{P}}\), el término dominante se desglosa como
\[ (P_n - P_0)\phi_{\hat{P}} = \underbrace{(P_n - P_0)\phi_{P_0}}_{\text{orden }n^{-1/2}} + (P_n - P_0)\bigl(\phi_{\hat{P}} - \phi_{P_0}\bigr). \]
El segundo término es la interacción impulsada por los parámetros de perturbación. Si entra en primer orden, incluso una mala especificación moderada de los modelos de perturbación puede introducir sesgo persistente.
1.5.2 Análisis de la tasa de convergencia
Supongamos que los estimadores de los parámetros de perturbación convergen a una tasa \(r_n\), es decir,
\[ \|\phi_{\hat{P}} - \phi_{P_0}\| = O_p(r_n). \]
Entonces el término adicional se comporta como
\[ (P_n - P_0)\bigl(\phi_{\hat{P}} - \phi_{P_0}\bigr) = O_p\!\bigl(n^{-1/2}r_n\bigr). \]
Sin una construcción especial, esto tiene el mismo orden que el término dominante, lo que significa que los errores de los parámetros de perturbación contaminan el comportamiento asintótico primario del estimador.
1.5.3 La ortogonalidad como remedio
La ortogonalidad de la función de influencia significa que tiene pendiente cero respecto a perturbaciones infinitesimales de las funciones de perturbación: cambios pequeños en esos modelos auxiliares no desplazan el parámetro objetivo en primer orden. Bajo esta condición, la contribución de los parámetros de perturbación se convierte en
\[ O_p\!\bigl(n^{-1/2}r_n^2\bigr). \]
Si \(r_n \to 0\) más rápido que \(n^{-1/4}\), entonces \(n^{-1/2}r_n^2 = o(n^{-1/2})\) y el error de estimación de los parámetros de perturbación es asintóticamente despreciable.
1.6 Resumen de los cuatro pilares conceptuales
| Concepto | Papel en TMLE |
|---|---|
| Funcionales | Separa el objetivo científico \(T(P)\) del conjunto de datos; el error de estimación es una diferencia funcional |
| Perturbaciones | Formaliza el muestreo como una desviación de medida con signo \(P_n - P_0\); permite la diferenciación en el espacio de distribuciones |
| Función de influencia | Proporciona una aproximación lineal de primer orden al error de estimación; determina la distribución asintótica |
| Ortogonalidad | Garantiza que los errores de estimación de los parámetros de perturbación entran solo en segundo orden; permite el uso de ML flexible para los modelos de perturbación |
2 Parte II. Simulando la Interacción de perturbación (Casi) Desvaneciente
2.1 Motivación
La Parte I argumentó que el término de interacción impulsado por los parámetros de perturbación
\[ (P_n - P_0)\bigl(\phi_{\hat{P}} - \phi_{P_0}\bigr) \]
debería reducirse a cero conforme crece el tamaño muestral, siempre que la función de influencia se haya construido para ser ortogonal a las perturbaciones de perturbación. El propósito de la Parte II es examinar esta afirmación empíricamente mediante datos simulados: ¿desaparece realmente el término de interacción y cómo cambia la tasa de desvanecimiento cuando los modelos de perturbación están deliberadamente mal especificados?
2.2 Recapitulación teórica
El error de estimación puede escribirse como
\[ T(P_n) - T(P_0) \approx (P_n - P_0)\phi_{P_0}. \]
Como la función de influencia verdadera nunca se observa, en la práctica se utiliza la versión estimada \(\phi_{\hat{P}}\). Descomponiendo la cantidad resultante se obtiene
\[ (P_n - P_0)\phi_{\hat{P}} = \underbrace{(P_n - P_0)\phi_{P_0}}_{\text{fluctuación bien comprendida}} + \underbrace{(P_n - P_0)\bigl(\phi_{\hat{P}} - \phi_{P_0}\bigr)}_{\text{interacción de perturbación}}. \]
El primer término es la fluctuación estocástica “buena” que impulsa la normalidad asintótica en \(\sqrt{n}\). El segundo término es la interacción de perturbación: la vía por la que los errores en los modelos auxiliares se filtran al comportamiento dominante del estimador. Cuando la función de influencia está correctamente construida, este término de filtración debería satisfacer
\[ (P_n - P_0)\bigl(\phi_{\hat{P}} - \phi_{P_0}\bigr) \to 0. \]
La simulación que sigue comprueba si esto es observable en muestras finitas, incluso bajo mala especificación del modelo de resultados.
2.3 La función de influencia eficiente para el ATE
2.3.1 Parámetro objetivo
Sea el objetivo el efecto promedio del tratamiento (ATE):
\[ \psi_0 = E_{P_0}[Y(1) - Y(0)]. \]
Bajo las condiciones de identificación estándar —consistencia, intercambiabilidad y positividad— el ATE puede escribirse como un funcional de la distribución de datos observados \(P_0\).
2.3.2 Funciones de perturbación auxiliares (modelo para la varible de resultado y para la variable de tratamiento)
Se definen la regresión de resultados y la puntuación de propensión (propensity score) como
\[ Q_a(X) = E[Y \mid A = a,\, X], \qquad g(X) = P(A = 1 \mid X). \]
2.3.3 Función de influencia eficiente
La función de influencia eficiente (EIF) para el ATE es
\[ \phi_{P_0}(Z) = \bigl(Q_1(X) - Q_0(X) - \psi_0\bigr) + \frac{A}{g(X)}\bigl(Y - Q_1(X)\bigr) - \frac{1-A}{1-g(X)}\bigl(Y - Q_0(X)\bigr), \tag{EIF} \]
donde las funciones de perturbación aparecen en dos roles distintos: directamente a través del término de sustitución \(Q_1(X) - Q_0(X)\), e indirectamente a través de los residuos ponderados por probabilidad inversa \(Y - Q_a(X)\). Esta estructura en capas hace que los errores en los modelos de perturbación entren por ambos canales, pero sus efectos de primer orden se cancelen —el sello distintivo de la ortogonalidad.
2.4 Una nota sobre la validación cruzada
Las simulaciones utilizan validación cruzada para estimar los parámetros de perturbación. Aunque no es estrictamente necesaria para TMLE, la validación cruzada es práctica estándar, especialmente con modelos flexibles. Sin ella, la fluctuación empírica \((P_n - P_0)\) y el error de perturbación \(\phi_{\hat{P}} - \phi_{P_0}\) son funciones de los mismos datos y comparten, por tanto, la misma aleatoriedad, creando un bucle de retroalimentación que puede distorsionar la condición de centrado de la EIF. La validación cruzada separa los roles de entrenamiento y evaluación, permitiendo que su interacción refleje mejor la cantidad teórica de interés.
2.5 Diseño de la simulación
2.5.1 Proceso generador de datos
El proceso generador de datos (PGD) crea una covariable binaria \(X_1\), una covariable continua \(X_2\), un tratamiento binario \(A\) impulsado por ambas covariables, y un resultado continuo \(Y\) que depende del tratamiento, las covariables y su interacción. El parámetro \(\tau\) es el ATE verdadero.
Mostrar código
gen_dgp <- function(n) {
def <-
defData(varname = "x1", formula = .5, dist = "binary") |>
defData(varname = "x2", formula = 0, variance = 1) |>
defData(
varname = "a",
formula = "-0.2 + 0.8 * x1 + 0.6 * x2",
dist = "binary", link = "logit"
) |>
defData(
varname = "y",
formula = "..tau * a + 1.0 * x1 + 1.0 * x2 + 1.5 * x1 * x2",
variance = 1, dist = "normal"
)
genData(n, def)[]
}2.5.2 Ajuste de los modelos de perturbación
Se examinan cuatro escenarios variando si cada modelo de perturbación está correctamente especificado o deliberadamente mal especificado.
Mostrar código
fit_nuisance <- function(dt, scenario) {
# Regresión de resultados Q(a, x)
if (scenario %in% c("both_correct", "g_wrong")) {
Q_fit <- lm(y ~ a + x1 + x2 + x1:x2, data = dt) # correcto
} else {
Q_fit <- lm(y ~ a + x1, data = dt) # mal especificado
}
# Puntuación de propensión g(x)
if (scenario %in% c("both_correct", "Q_wrong")) {
g_fit <- glm(a ~ x1 + x2, data = dt, family = binomial()) # correcto
} else {
g_fit <- glm(a ~ x1, data = dt, family = binomial()) # mal especificado
}
list(Q_fit = Q_fit, g_fit = g_fit)
}2.5.3 Predicciones y funciones de perturbación verdaderas
Mostrar código
predict_Q <- function(Q_fit, dt, a_val) {
nd <- copy(dt); nd[, a := a_val]
as.numeric(predict(Q_fit, newdata = nd))
}
predict_g <- function(g_fit, dt) {
p <- as.numeric(predict(g_fit, newdata = dt, type = "response"))
pmin(pmax(p, 0.01), 0.99) # estabilización de pesos extremos: truncado
}
Q_true <- function(dt, a_val, tau) {
tau * a_val + 1.0 * dt$x1 + 1.0 * dt$x2 + 1.5 * dt$x1 * dt$x2
}
g_true <- function(dt) {
plogis(-0.2 + 0.8 * dt$x1 + 0.6 * dt$x2)
}2.5.4 Construcción de la EIF
Se construyen dos versiones de la EIF: (i) estimada, usando los parámetros de perturbación ajustados; (ii) oracle, usando el PGD conocido. La EIF estimada requiere un estimador de centrado compatible \(\hat\psi\), calculado por validación cruzada a partir de los mismos ajustes de perturbación para garantizar media aproximada cero.
Mostrar código
phi_ate <- function(dt, Q1, Q0, g, psi) {
A <- dt$a; Y <- dt$y
(Q1 - Q0 - psi) + A / g * (Y - Q1) - (1 - A) / (1 - g) * (Y - Q0)
}
psi_hat_from_fits <- function(dt, Q_fit, g_fit) {
Q1 <- predict_Q(Q_fit, dt, 1); Q0 <- predict_Q(Q_fit, dt, 0)
g <- predict_g(g_fit, dt)
A <- dt$a; Y <- dt$y
mean((Q1 - Q0) + A / g * (Y - Q1) - (1 - A) / (1 - g) * (Y - Q0))
}
phi_hat_from_fits <- function(dt, Q_fit, g_fit, psi_hat) {
Q1 <- predict_Q(Q_fit, dt, 1); Q0 <- predict_Q(Q_fit, dt, 0)
g <- predict_g(g_fit, dt)
phi_ate(dt, Q1, Q0, g, psi_hat)
}
phi0_true <- function(dt, tau) {
Q1 <- Q_true(dt, 1, tau); Q0 <- Q_true(dt, 0, tau)
g <- g_true(dt)
phi_ate(dt, Q1, Q0, g, psi = tau)
}2.5.5 Estimación del término de interacción de perturbación
Para cada conjunto de datos se realizan los siguientes pasos: (1) división en dos particiones; (2) ajuste de los modelos de perturbación en cada partición de entrenamiento; (3) cálculo de la EIF con validación cruzada; (4) comparación de la EIF estimada con la EIF oracle tanto en la muestra como en una extracción poblacional independiente, aproximando \((P_n - P_0)(\phi_{\hat{P}} - \phi_{P_0})\).
Mostrar código
est_2T <- function(scenario, dd, tau, dd_pop) {
n <- nrow(dd)
idx <- sample.int(n)
I1 <- idx[seq_len(floor(n / 2))]
I2 <- idx[(floor(n / 2) + 1):n]
fits1 <- fit_nuisance(dd[I1], scenario)
fits2 <- fit_nuisance(dd[I2], scenario)
psi1 <- psi_hat_from_fits(dd[I2], fits1$Q_fit, fits1$g_fit)
psi2 <- psi_hat_from_fits(dd[I1], fits2$Q_fit, fits2$g_fit)
psi_hat_cf <- 0.5 * (psi1 + psi2)
phi_hat_dd <- numeric(n)
phi_hat_dd[I2] <- phi_hat_from_fits(dd[I2], fits1$Q_fit, fits1$g_fit, psi_hat_cf)
phi_hat_dd[I1] <- phi_hat_from_fits(dd[I1], fits2$Q_fit, fits2$g_fit, psi_hat_cf)
dphi_dd <- phi_hat_dd - phi0_true(dd, tau)
dphi_pop <-
0.5 * (phi_hat_from_fits(dd_pop, fits1$Q_fit, fits1$g_fit, psi_hat_cf) +
phi_hat_from_fits(dd_pop, fits2$Q_fit, fits2$g_fit, psi_hat_cf)) -
phi0_true(dd_pop, tau)
data.table(scenario, n, T2 = mean(dphi_dd) - mean(dphi_pop))[]
}2.5.6 Ejecución de la simulación
El bucle externo sobre las 2.000 réplicas se paraleliza con future/furrr, distribuyendo el trabajo entre todos los núcleos disponibles. furrr::future_map es un reemplazo directo de purrr::map que respeta un backend paralelo de future, y furrr_options(seed = TRUE) propaga semillas por trabajador para garantizar la reproducibilidad.
Mostrar código
library(future)
library(furrr)
run_sim <- function(n, tau, dd_pop, scenarios) {
dd <- gen_dgp(n)
rbindlist(lapply(scenarios, function(s) est_2T(s, dd, tau, dd_pop)))
}
# Usar todos los núcleos disponibles menos uno; retrocede a secuencial si solo hay 1
n_cores <- max(1L, parallelly::availableCores() - 1L)
plan(multisession, workers = n_cores)
set.seed(1)
tau <- 5
pop_dd <- gen_dgp(5e5)
n_vec <- rep(c(100, 250, 750, 1000), each = 500)
scenarios <- c("both_correct", "Q_wrong")
res <- rbindlist(
future_map(
n_vec,
function(x) run_sim(x, tau, pop_dd, scenarios),
.options = furrr_options(seed = TRUE)
)
)
# Volver al plan secuencial tras la simulación
plan(sequential)2.6 Resultados
La Figure 1 muestra las estimaciones individuales del término de interacción de perturbación \((P_n - P_0)(\phi_{\hat{P}} - \phi_{P_0})\) a lo largo de los tamaños muestrales y los escenarios de mala especificación.
Mostrar código
scenario_labels <- c(
both_correct = "Ambos modelos correctos",
Q_wrong = "Modelo de resultados mal especificado"
)
res[, scenario_label := factor(
scenario_labels[scenario],
levels = unname(scenario_labels)
)]
res[, n_factor := factor(n, levels = sort(unique(n)))]
ggplot(res, aes(x = n_factor, y = T2)) +
geom_jitter(
width = 0.25, alpha = 0.25, size = 0.8,
colour = "#3B82F6"
) +
geom_boxplot(
aes(group = n_factor),
outlier.shape = NA, fill = NA,
colour = "#1E3A5F", linewidth = 0.5, width = 0.45
) +
geom_hline(
yintercept = 0, linetype = "dashed",
colour = "#EF4444", linewidth = 0.7
) +
facet_wrap(~ scenario_label, scales = "free_y") +
labs(
x = "Tamaño muestral (n)",
y = expression((P[n] - P[0])(phi[hat(P)] - phi[P[0]])),
title = "Término de interacción de perturbación por tamaño muestral",
subtitle = "Los puntos son réplicas individuales; la caja muestra el IQR y la mediana"
) +
theme_bw(base_size = 12) +
theme(
strip.background = element_rect(fill = "#EFF6FF"),
strip.text = element_text(face = "bold"),
panel.grid.minor = element_blank()
)
Modelos correctamente especificados. Las estimaciones están ajustadamente centradas alrededor de cero incluso con \(n = 100\); la variabilidad disminuye rápidamente y es esencialmente despreciable a partir de \(n = 750\). La desviación típica cae de \(0{,}28\) con \(n = 100\) a \(0{,}006\) con \(n = 1000\).
Modelo de resultados mal especificado. La variabilidad es marcadamente mayor en todos los tamaños muestrales (desviación típica \(1{,}34\) con \(n = 100\)), reflejando un error de estimación de perturbación persistente. Sin embargo, las estimaciones permanecen centradas alrededor de cero y la dispersión disminuye claramente con \(n\), cayendo a \(0{,}12\) con \(n = 1000\).
2.7 Limitaciones y transición a la Parte III
El hecho de que la interacción se reduzca no garantiza que sea pequeña relativa al error de muestreo para ningún tamaño muestral dado. La tasa de convergencia depende de la precisión con que mejoran los modelos de perturbación y, sin garantías formales sobre esa tasa, la discrepancia restante puede seguir afectando a la inferencia en muestras finitas. La ortogonalidad garantiza que los errores de perturbación entran únicamente a través de su interacción con la variabilidad muestral, pero no fuerza a que la ecuación empírica de la EIF se cumpla exactamente. Ese es precisamente el papel del paso de dirección de TMLE, que es el objeto de la Parte III.
3 Parte III. Forzando al Objetivo a Comportarse: El Paso de Dirección de TMLE
3.1 Visión general
Las Partes I y II establecieron que los estimadores TMLE deben poseer una función de influencia ortogonal a las perturbaciones de perturbación, y que el término de interacción de perturbación sí se reduce a cero en la simulación. Lo que quedó sin abordar es cómo TMLE hace cumplir activamente la condición de centrado \(P_n\phi_{\hat{P}} \approx 0\) sobre la EIF estimada. Esta parte deriva el paso de dirección (fluctuación) que lo consigue.
La idea central es la siguiente. Si las funciones de perturbación verdaderas fueran conocidas, la media empírica de la función de influencia ya sería (aproximadamente) cero, y el estimador diferiría de la verdad solo a través de la variabilidad muestral. En la práctica, los modelos de perturbación estimados alteran este equilibrio. Incluso errores pequeños pueden desplazar la EIF empírica lejos de cero, introduciendo sesgo impulsado por el modelo que persiste junto con el ruido de muestreo. TMLE elimina este desequilibrio realizando una actualización unidimensional y dirigida de la regresión de resultados inicial, sin reajustar los modelos de perturbación desde cero.
3.2 Fundamento causal
Observamos datos \(Z_i = (X_i, A_i, Y_i)\), \(i = 1,\ldots,n\), que comprenden covariables basales \(X\), un indicador de tratamiento binario \(A \in \{0,1\}\), y un resultado \(Y\). Cada unidad posee dos resultados potenciales \(Y(1),\, Y(0)\), que representan los resultados que se realizarían bajo tratamiento y control respectivamente. El efecto promedio del tratamiento es
\[ \psi^0 = E[Y(1) - Y(0)]. \]
Bajo las condiciones estándar de identificabilidad —consistencia (\(Y = Y(A)\)), intercambiabilidad (\(A \perp\!\!\!\perp Y(a) \mid X\)) y positividad (\(g(X) > 0\) c.s.)— el ATE se identifica a partir de la distribución de datos observados como
\[ \psi(P) = E_P[Q_1(X) - Q_0(X)], \qquad Q_a(x) = E[Y \mid A = a, X = x], \]
con \(g(x) = P(A=1 \mid X=x)\) denotando la puntuación de propensión. El parámetro objetivo es \(\psi_0 = \psi(P_0)\).
3.3 La EIF y la condición de centrado
3.3.1 EIF para el ATE
Siguiendo la construcción de perturbación y diferenciación de la Parte I, la función de influencia eficiente para \(\psi(P)\) es
\[ \phi_P(Z) = \bigl(Q_1(X) - Q_0(X) - \psi(P)\bigr) + \frac{A}{g(X)}\bigl(Y - Q_1(X)\bigr) - \frac{1-A}{1-g(X)}\bigl(Y - Q_0(X)\bigr). \tag{EIF} \]
Cuando se evalúa en los parámetros de perturbación verdaderos, esta función está centrada:
\[ E_{P_0}[\phi_{P_0}(Z)] = 0. \]
3.3.2 Fallo del centrado con parámetros de perturbación estimados
Denotemos las estimaciones iniciales de los parámetros de perturbación como \(\hat{Q}^0_a\) y \(\hat{g}\), y sea \(\hat\psi^0 = P_n(\hat{Q}^0_1 - \hat{Q}^0_0)\). La EIF estimada en el ajuste inicial típicamente tiene media empírica no nula:
\[ P_n\phi_{\hat{P}^0} = \frac{1}{n}\sum_{i=1}^{n}\phi_{\hat{P}^0}(Z_i) \neq 0. \]
Cuando la EIF empírica no está centrada, la expansión de primer orden ideal se ve perturbada: el estimador se aleja de la verdad no solo por el ruido de muestreo, sino también por el desalineamiento introducido por los modelos de perturbación imperfectos.
3.4 Identificación de la fuente del desequilibrio
Por construcción, \(P_n(\hat{Q}^0_1(X) - \hat{Q}^0_0(X) - \hat\psi^0) = 0\), de modo que cualquier desequilibrio en \(P_n\phi_{\hat{P}^0}\) proviene íntegramente de los términos de corrección residual. Introduciendo la covariable inteligente
\[ H_{\hat{g}}(A,X) = \frac{A}{\hat{g}(X)} - \frac{1-A}{1-\hat{g}(X)}, \]
la parte dependiente del resultado de la EIF se convierte en \(H_{\hat{g}}(A,X)\bigl(Y - \hat{Q}^0(A,X)\bigr)\), y el desequilibrio se reduce a
\[ P_n\Bigl[H_{\hat{g}}(A,X)\bigl(Y - \hat{Q}^0(A,X)\bigr)\Bigr] \neq 0. \]
La única palanca disponible es el residuo \(Y - \hat{Q}^0(A,X)\): ajustar ligeramente la regresión de resultados puede llevar este desequilibrio a cero.
3.5 El paso de fluctuación
TMLE introduce una fluctuación unidimensional del ajuste inicial de resultados:
\[ \hat{Q}^{\epsilon}(A,X) = \hat{Q}^0(A,X) + \epsilon\, H_{\hat{g}}(A,X). \]
El escalar \(\epsilon\) se estima regresando \(Y\) sobre \(H_{\hat{g}}(A,X)\) con desplazamiento \(\hat{Q}^0(A,X)\). La ecuación normal de mínimos cuadrados es
\[ P_n\Bigl[H_{\hat{g}}(A,X)\bigl(Y - \hat{Q}^{\hat\epsilon}(A,X)\bigr)\Bigr] = 0. \]
Definiendo \(\hat{Q}^*(A,X) = \hat{Q}^{\hat\epsilon}(A,X)\), la estimación TMLE del ATE es
\[ \hat\psi^* = \frac{1}{n}\sum_{i=1}^{n} \bigl(\hat{Q}^*(1,X_i) - \hat{Q}^*(0,X_i)\bigr). \]
3.6 La EIF actualizada está centrada
La EIF dirigida es
\[ \phi_{\hat{P}^*}(Z) = \underbrace{\bigl(\hat{Q}^*_1(X) - \hat{Q}^*_0(X) - \hat\psi^*\bigr)}_{\text{término de sustitución}} + \underbrace{H_{\hat{g}}(A,X)\bigl(Y - \hat{Q}^*(A,X)\bigr)}_{\text{residuo ponderado}}. \]
Por la construcción de sustitución, la media empírica del primer término es exactamente cero; por la ecuación normal, la media empírica del segundo término también es cero. Por tanto \(P_n\phi_{\hat{P}^*} \approx 0\). El paso de dirección ha restaurado la condición de centrado en la muestra observada.
3.7 Por qué la interacción de perturbación se vuelve despreciable
Partiendo de la descomposición
\[ (P_n - P_0)(\phi_{\hat{P}^*} - \phi_{P_0}) = P_n(\phi_{\hat{P}^*} - \phi_{P_0}) - P_0(\phi_{\hat{P}^*} - \phi_{P_0}), \]
el paso de dirección impone \(P_n\phi_{\hat{P}^*} \approx 0\), por lo que el primer término se reduce a \(-P_n\phi_{P_0}\), que fluctúa al orden \(n^{-1/2}\). El segundo término refleja el error de perturbación a nivel poblacional y, debido a la ortogonalidad, se comporta como el producto \(e_Q \times e_g\) de los errores de estimación respectivos, reduciéndose a lo insignificante a medida que mejoran los parámetros de perturbación. El paso de dirección elimina el desequilibrio dominante; lo que queda está dominado por \((P_n - P_0)\phi_{P_0}\).
3.8 Doble robustez
TMLE hereda la doble robustez de la estructura de la EIF. Los errores de perturbación entran en el estimador de forma multiplicativa: el sesgo dominante se comporta como el producto \(e_Q \times e_g\). En consecuencia, si cualquiera de los modelos de perturbación se estima de forma consistente, \(e_Q \times e_g \to 0\) y el sesgo desaparece. Incluso cuando ambos modelos son imperfectos, su producto puede ser pequeño si los errores no son grandes simultáneamente.
3.9 El papel de la validación cruzada
El paso de dirección puede aplicarse con o sin validación cruzada. Sin embargo, cuando \(\hat{Q}^0\) y \(\hat{g}\) se estiman con métodos flexibles de aprendizaje automático, la validación cruzada ayuda a garantizar que la ecuación empírica de la EIF se comporte como predice la teoría. Sin ella, las mismas observaciones que entrenan los modelos de perturbación también evalúan la función de influencia, creando un bucle de retroalimentación que puede distorsionar la condición de centrado. La validación cruzada separa el entrenamiento de la evaluación, de modo que cada observación se evalúa con estimaciones de perturbación aprendidas a partir de otros datos, permitiendo que la asintótica habitual de la función de influencia aparezca con mayor claridad en muestras finitas.
3.10 Resumen de la serie en tres partes
La Table 2 recoge las ideas principales a lo largo de las tres partes.
| Componente | Qué hace | Por qué importa |
|---|---|---|
| Funcionales | Define el objetivo como \(T(P)\) | Separa el objetivo científico de los datos |
| Perturbaciones | Modela \(P_n - P_0\) como medida con signo | Permite la diferenciación en el espacio de distribuciones |
| Función de influencia | Aproximación de primer orden a \(T(P_n) - T(P_0)\) | Determina la distribución asintótica |
| Ortogonalidad | La EIF tiene pendiente cero respecto a los parámetros de perturbación | Permite estimar con ML; los errores entran en 2.º orden |
| Covariable inteligente \(H_{\hat g}\) | Codifica el residuo ponderado por propensión | Hace visible y corregible el desequilibrio |
| Paso de fluctuación | \(\hat{Q}^\epsilon = \hat{Q}^0 + \epsilon H_{\hat g}\) | Actualización unidimensional que restaura el centrado de la EIF |
| Doble robustez | Sesgo \(\propto e_Q \times e_g\) | Consistente si cualquiera de los modelos de perturbación está correctamente especificado |
| Validación cruzada | Separa entrenamiento de evaluación | Evita retroalimentación; mejora el centrado de la EIF en muestras finitas |
Referencias
Van der Laan, Mark J., and Sherri Rose. Targeted Learning: Causal Inference for Observational and Experimental Data, vol. 4. New York: Springer, 2011.
Hampel, Frank R. “The influence curve and its role in robust estimation.” Journal of the American Statistical Association 69, no. 346 (1974): 383–393. https://doi.org/10.1080/01621459.1974.10482962
Chernozhukov, Victor, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins. “Double/debiased machine learning for treatment and structural parameters.” The Econometrics Journal 21, no. 1 (2018): C1–C68. https://doi.org/10.1111/ectj.12097
Kennedy, Edward H. “Semiparametric doubly robust targeted double machine learning: a review.” arXiv preprint arXiv:2203.06469 (2022).
Luque-Fernandez, Miguel Angel, Michael Schomaker, Bernard Rachet, and Mireille Schnitzer. “Targeted maximum likelihood estimation for a binary treatment: A tutorial.” Statistics in Medicine 37, no. 16 (2018): 2530–2546. https://doi.org/10.1002/sim.7628
Schuler, Megan S., and Sherri Rose. “Targeted maximum likelihood estimation for causal inference in observational studies.” American Journal of Epidemiology 185, no. 1 (2017): 65–73. https://doi.org/10.1093/aje/kww165