Resumen: Este manual guía una clase intensiva de inferencia causal aplicada con R, usando el caso paradigmático “cigüeñas y nacimientos” para ilustrar correlación espuria y métodos formales. Integra teoría (DAGs, ecuaciones estructurales, OLS, IV/2SLS, matching, DiD, RDD), práctica en R y diagnóstico riguroso. Inspirado en Pearl (2009), Hernán & Robins (2020), Imbens & Rubin (2015), Angrist & Pischke (2008), Greenland et al. (1999), y Rosenbaum & Rubin (1983).
¡Bienvenidos! Hoy nos convertiremos en detectives de datos. Nuestra misión: resolver el misterio de la persistente correlación entre el número de cigüeñas y el número de bebés nacidos.
La pregunta clave no es si estas dos variables se mueven juntas, sino si una intervención sobre las cigüeñas (ej. una política de protección que aumente su población) tendría un efecto causal real sobre los nacimientos.
Para cada país i, existen dos realidades potenciales o contrafactuales:
Y_i(1): El número de bebés que nacerían si el país tuviera “muchas” cigüeñas (el resultado potencial bajo “tratamiento”).
Y_i(0): El número de bebés que nacerían si el país tuviera “pocas” cigüeñas (el resultado potencial bajo “control”).
El efecto causal individual es la diferencia: Y_i(1)−Y_i(0). El problema fundamental es que nunca podemos observar ambos escenarios a la vez para el mismo país. La inferencia causal es el conjunto de herramientas que nos permite estimar el efecto promedio de estos contrafactuales a partir de los datos que sí observamos.
Nuestra hipótesis, basada en el conocimiento del mundo, es que el verdadero efecto causal es cero. Veamos si los datos intentan engañarnos.
library(ggdag)
library(dagitty)
# DAG que muestra la estructura de confusión
dag <- dagify(
births ~ storks )
tidy_dagitty(dag) %>%
ggdag_classic(text_col = "black", stylized = TRUE) +
labs(title = "DAG:") +
theme_dag()
La correlación entre el número de cigüeñas y nacimientos es un ejemplo clásico de correlación espuria: dos variables muestran alta asociación, pero no existe relación causal directa. El reto central de la inferencia causal es distinguir efectos causales verdaderos de asociaciones generadas por confusores o por el azar.
Pregunta causal: ¿Cuál sería el resultado \(Y\) si \(X\) tomara otro valor, manteniendo todo lo demás constante? Buscamos \(\mathbb{E}[Y\mid do(X=x)]\).
Supongamos [ Y_i = X_i + _1 C_i + u_i, ] donde \(Y_i\) (nacimientos), \(X_i\) (cigüeñas), \(C_i\) (confusores: tamaño, ruralidad), y \(u_i\) errores no observados.
El Problema Fundamental de la Inferencia Causal: para un individuo no podemos observar a la vez el resultado tratado y no tratado. Trabajamos con contrafactuales y supuestos de identificación.
Para entender cómo surgen las correlaciones engañosas, vamos a simular un mundo donde nosotros definimos las reglas. En nuestro universo, la verdad absoluta es que no hay relación causal entre cigüeñas y bebés. En cambio, una tercera variable, el tamaño del país (área), es una causa común de ambos.
La Verdadera Estructura Causal: El DAG
Podemos visualizar la estructura de nuestro universo con un Diagrama Acíclico Dirigido (DAG).
library(ggdag)
library(dagitty)
# DAG que muestra la estructura de confusión
dag <- dagify(
births ~ area,
storks ~ area,
exposure = "storks",
outcome = "births",
labels = c("storks" = "Cigüeñas",
"births" = "Nacimientos",
"area" = "Área del País")
)
tidy_dagitty(dag) %>%
ggdag_classic(text_col = "black", stylized = TRUE) +
labs(title = "DAG: La Verdadera Estructura Causal de Nuestro Mundo") +
theme_dag()
En nuestro DAG, existe un camino no causal que conecta “Cigüeñas” y “Nacimientos”:
Cigüeñas ← Área del País → Nacimientos
Este es un camino de puerta trasera. Si no lo “bloqueamos” estadísticamente, la asociación que fluye a través de él se mezclará con el efecto causal (que sabemos que es cero). Esto crea una correlación espuria. Nuestro principal objetivo como investigadores es identificar y cerrar estas puertas traseras.
n <- 50 # Número de países en nuestro universo
area_km2 <- runif(n, 500, 300000) # El CONFUSOR
# Las cigüeñas son una función del área + ruido aleatorio
storks <- round(area_km2 / 5000 + rnorm(n, 0, 50), 0)
# Los nacimientos también son una función del área + ruido aleatorio
births <- round(0.1 * area_km2 + rnorm(n, 0, 5000), 0)
# Variable irrelevante
gdp_per_capita <- runif(n, 2000, 50000)
df <- tibble(
country = paste("Country", 1:n),
storks, births, area_km2, gdp_per_capita
)
kable(head(df, 10), caption = "Primeros 10 países de nuestro universo simulado") %>%
kable_styling(bootstrap_options = "striped", full_width = FALSE)| country | storks | births | area_km2 | gdp_per_capita |
|---|---|---|---|---|
| Country 1 | 121 | 22706 | 219030.98 | 2723.956 |
| Country 2 | -16 | 31637 | 270769.63 | 41352.118 |
| Country 3 | -33 | 17279 | 176817.43 | 33871.782 |
| Country 4 | -9 | 6661 | 63093.35 | 14869.084 |
| Country 5 | -1 | 21249 | 210652.87 | 21814.977 |
| Country 6 | 22 | 8786 | 86203.31 | 16141.046 |
| Country 7 | 9 | -5683 | 27877.59 | 48433.998 |
| Country 8 | 32 | 8707 | 40270.88 | 5592.618 |
| Country 9 | 21 | -4844 | 19816.06 | 46834.160 |
| Country 10 | 45 | 29969 | 287945.15 | 17974.929 |
# No puede haber cigüeñas negativas
df<- df %>%
mutate(
storks=if_else(storks<0,0,storks),
births=if_else(births<0,0,births)
)# Visualización de la correlación aparente
ggplot(df, aes(storks)) +
geom_histogram()+theme_classic() 
# Calculamos el promedio para usarlo en el gráfico
mean_births <- mean(df$births)
ggplot(df, aes(x = births)) +
# 1. Histograma con estética mejorada y mapeo a densidad
geom_histogram(
aes(y = after_stat(density)), # Mapea el eje Y a densidad en lugar de conteo
fill = "steelblue", # Color de relleno
color = "white", # Color del borde
alpha = 0.7 # Transparencia
) +
labs(
title = "Distribución del Número de Nacimientos por País",
x = "Número de Nacimientos",
y = "Densidad"
) +
theme_minimal(base_size = 14)
# Visualización de la correlación aparente
ggplot(df, aes(storks, births)) +
geom_point(color = "darkblue", alpha = 0.7) +
geom_smooth(method = "lm", se = FALSE, color = "red", linetype = "dashed") +
labs(
title = "Correlación Aparente: Una Evidencia Engañosa",
subtitle = "Una fuerte relación positiva es visible a simple vista",
x = "Número de Cigüeñas", y = "Nacimientos Anuales"
)
# Regresión lineal simple (Modelo Ingenuo)
lm_naive <- lm(births ~ storks, data = df)
summary(lm_naive)
Call:
lm(formula = births ~ storks, data = df)
Residuals:
Min 1Q Median 3Q Max
-16538.4 -5342.0 -557.7 6137.2 18665.4
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12971.6 1846.3 7.026 6.72e-09 ***
storks 102.3 37.6 2.721 0.00903 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8809 on 48 degrees of freedom
Multiple R-squared: 0.1337, Adjusted R-squared: 0.1156
F-statistic: 7.406 on 1 and 48 DF, p-value: 0.009029# Instala el paquete si no lo tienes
# install.packages("performance")
# Carga la librería
library(performance, see)
# ✅ Ejecuta el diagnóstico completo con una sola función
performance::check_model(lm_naive)
Al ejecutar check_model(lm_naive), obtendrás un panel con varias gráficas. Las más importantes son:
Linearity (Linealidad): Los puntos (residuos) deben distribuirse aleatoriamente alrededor de la línea horizontal en cero. Si siguen un patrón curvo, la relación podría no ser lineal.
Homoscedasticity (Homocedasticidad): La dispersión de los puntos debe ser constante a lo largo del eje X. Si ves una forma de cono o embudo, es un signo de heterocedasticidad (lo que justifica el uso de errores robustos).
Normality of Residuals (Normalidad de los Residuos): Los puntos deben seguir de cerca la línea diagonal. Si se desvían mucho, los residuos no se distribuyen de forma normal.
Influential Observations (Observaciones Influyentes): Destaca los puntos que podrían tener una influencia desproporcionada en los resultados de la regresión (outliers o puntos de alto apalancamiento).
# Ahora creamos los gráficos de diagnóstico
plot(lm_naive,1)
plot(lm_naive,2)
plot(lm_naive,3)
plot(lm_naive,4)
Residuals vs Fitted (Residuos vs. Valores Ajustados)plot(lm_naive, 1)
¿Para qué sirve? Para comprobar dos supuestos a la vez:
Linealidad: ¿La relación entre las variables es realmente lineal?
Homocedasticidad: ¿La varianza de los errores es constante?
✅ ¿Qué buscar? (La buena señal)
Una nube de puntos dispersa aleatoriamente alrededor de la línea horizontal discontinua en y=0. No debe haber ningún patrón claro.
La línea roja (que suaviza la tendencia de los puntos) debe ser casi plana y estar cerca del cero.
❌ ¿Qué es una mala señal?
Patrón curvo: Si la línea roja tiene forma de “U” o de “U” invertida, sugiere que la relación no es lineal. Quizás necesitas incluir un término cuadrático (x^2) en tu modelo.
Forma de cono o embudo: Si los puntos se dispersan más a medida que te mueves hacia la derecha (o hacia la izquierda), es un signo de heterocedasticidad. Esto significa que la varianza de los errores no es constante, y los errores estándar de tu modelo no son fiables (¡lo que justifica usar errores robustos!)
# Instala el paquete si no lo tienes
# install.packages("lmtest")
library(lmtest)
# Ejecuta el test de Breusch-Pagan sobre el modelo
bptest(lm_naive)
studentized Breusch-Pagan test
data: lm_naive
BP = 1.9948, df = 1, p-value = 0.1578Hipótesis:
H_0 (Hipótesis Nula): La varianza de los residuos es constante (homocedasticidad).
H_1 (Hipótesis Alternativa): La varianza de los residuos no es constante (heterocedasticidad).
Normal Q-Q (Gráfico Quantil-Quantil Normal)plot(lm_naive, 2)
¿Para qué sirve? Para comprobar si los residuos siguen una distribución normal.
✅ ¿Qué buscar? (La buena señal)
❌ ¿Qué es una mala señal?
¿Qué significa si hay un problema? Si los residuos no son normales, los intervalos de confianza y los p-valores pueden no ser del todo precisos, sobre todo en muestras pequeñas.
# Extraemos los residuos del modelo
residuos_modelo <- residuals(lm_naive)
# Ejecutamos el test de Shapiro-Wilk sobre los residuos
shapiro.test(residuos_modelo)
Shapiro-Wilk normality test
data: residuos_modelo
W = 0.98218, p-value = 0.6468Hipótesis:
H_0: Los residuos siguen una distribución normal.
H_1: Los residuos no siguen una distribución normal.
Scale-Location (Escala-Ubicación)plot(lm_naive, 3)
¿Para qué sirve? Es otra forma, a veces más clara, de detectar la homocedasticidad (varianza constante de los errores). En el eje Y no están los residuos, sino la raíz cuadrada de los residuos estandarizados, lo que ayuda a visualizar mejor la varianza.
✅ ¿Qué buscar? (La buena señal)
Una nube de puntos dispersa aleatoriamente y sin patrones.
La línea roja debe ser aproximadamente horizontal.
❌ ¿Qué es una mala señal?
# Ejecuta el test RESET de Ramsey sobre el modelo
resettest(lm_naive)
RESET test
data: lm_naive
RESET = 0.25708, df1 = 2, df2 = 46, p-value = 0.7744H_0: El modelo está correctamente especificado (la relación es lineal).
H_1: El modelo está mal especificado (la relación no es lineal).
Residuals vs Leverage (Residuos vs. Apalancamiento)plot(lm_naive, 4)
¿Para qué sirve? Para identificar observaciones influyentes. Un punto influyente es aquel que, si se elimina, cambia drásticamente el resultado de la regresión.
Conceptos clave:
Leverage (Apalancamiento) (eje X): Mide qué tan atípico es un punto en función de sus predictores (la variable x). Un punto con alto apalancamiento está muy lejos de la media de los valores de x.
Standardized Residuals (eje Y): Mide qué tan grande es el error de predicción para ese punto.
Distancia de Cook (líneas rojas discontinuas): Es una métrica que combina el residuo y el apalancamiento para medir la influencia total de un punto.
✅ ¿Qué buscar? (La buena señal)
❌ ¿Qué es una mala señal?
# Calculamos la Distancia de Cook para cada observación
distancias_cook <- cooks.distance(lm_naive)
# Definimos el umbral de influencia
n <- nobs(lm_naive) # nobs() obtiene el número de observaciones
umbral <- 4/n
# Vemos qué observaciones superan el umbral
which(distancias_cook > umbral)2
2 # install.packages("estimatr")
library(estimatr)
# AJUSTAMOS EL MODELO ESPECIFICANDO EL TIPO DE ERROR ESTÁNDAR DE STATA
modelo_stata_replica <- lm_robust(births ~ storks, data = df, se_type = "stata")
# Imprimimos el resultado
modelo_stata_replica Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper
(Intercept) 12971.6050 1934.927 6.703924 2.092307e-08 9081.17321 16862.0367
storks 102.3277 32.598 3.139079 2.897258e-03 36.78501 167.8703
DF
(Intercept) 48
storks 48Este modelo sufre de un grave sesgo de variable omitida. El coeficiente estimado para storks no representa su efecto causal real. En cambio, está contaminado, capturando el efecto de la variable omitida
# Modelo ajustado
lm_adjusted <- lm(births ~ storks + area_km2 + gdp_per_capita, data = df)
summary(lm_adjusted)
Call:
lm(formula = births ~ storks + area_km2 + gdp_per_capita, data = df)
Residuals:
Min 1Q Median 3Q Max
-9842.8 -3038.1 -538.1 2460.8 10596.9
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.914e+03 2.091e+03 1.394 0.170
storks 1.280e+01 2.251e+01 0.569 0.572
area_km2 8.893e-02 8.078e-03 11.009 1.75e-14 ***
gdp_per_capita -2.070e-02 5.864e-02 -0.353 0.726
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4688 on 46 degrees of freedom
Multiple R-squared: 0.7649, Adjusted R-squared: 0.7495
F-statistic: 49.88 on 3 and 46 DF, p-value: 1.677e-14# Después de correr lm_naive y lm_adjusted
library(modelsummary)
# Crea una lista con los modelos que quieres comparar
lista_de_modelos <- list(
"Modelo Ingenuo" = lm_naive,
"Modelo Ajustado" = lm_adjusted
)
# Crea la tabla
modelsummary(lista_de_modelos, stars = TRUE, gof_map = "nobs,r.squared")| Modelo Ingenuo | Modelo Ajustado | |
|---|---|---|
| (Intercept) | 12971.605*** | 2914.497 |
| (1846.323) | (2091.053) | |
| storks | 102.328** | 12.802 |
| (37.600) | (22.513) | |
| area_km2 | 0.089*** | |
| (0.008) | ||
| gdp_per_capita | -0.021 | |
| (0.059) | ||
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
df <- df %>%
mutate(high_storks = if_else(storks > median(storks), "Muchas Cigüeñas", "Pocas Cigüeñas"))
df %>%
group_by(high_storks) %>%
summarise( N_Paises = n(), Area_Promedio_km2 = mean(area_km2), Nacimientos_Promedio = mean(births) ) # A tibble: 2 × 4
high_storks N_Paises Area_Promedio_km2 Nacimientos_Promedio
<chr> <int> <dbl> <dbl>
1 Muchas Cigüeñas 25 177303. 18982.
2 Pocas Cigüeñas 25 134114. 14378.Observación: Vemos un claro desequilibrio. Los países con “Muchas Cigüeñas” son, en promedio, mucho más grandes. Esto confirma que no podemos comparar directamente los grupos. El ajuste por regresión es lo que nos permite hacer una comparación justa, controlando por estas diferencias preexistentes.
df<-df %>%
mutate(
red = rbinom(n, 1, 0.5),
blue = rbinom(n, 1, 0.5),
white = rbinom(n, 1, 0.5),
green = rbinom(n, 1, 0.5),
yellow = rbinom(n, 1, 0.5),
black = rbinom(n, 1, 0.5)
)
lm_bad_adjusted <- lm(births ~ storks + gdp_per_capita +red +blue +white +green+yellow+black , data = df)
summary(lm_bad_adjusted)
Call:
lm(formula = births ~ storks + gdp_per_capita + red + blue +
white + green + yellow + black, data = df)
Residuals:
Min 1Q Median 3Q Max
-16041.7 -5428.8 -703.8 6208.3 16801.0
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13646.5286 4119.5158 3.313 0.00194 **
storks 126.1587 39.8018 3.170 0.00288 **
gdp_per_capita 0.1650 0.1179 1.399 0.16941
red 1090.2778 2668.1967 0.409 0.68495
blue -3113.0047 2593.4925 -1.200 0.23691
white -4369.2541 2638.5407 -1.656 0.10537
green -1318.5649 2757.7260 -0.478 0.63509
yellow -1679.4062 2709.0579 -0.620 0.53874
black -3952.7979 2727.0162 -1.449 0.15481
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8866 on 41 degrees of freedom
Multiple R-squared: 0.2504, Adjusted R-squared: 0.1042
F-statistic: 1.712 on 8 and 41 DF, p-value: 0.1246df<-df %>%
mutate(
stadistical_register = births + storks + rnorm(n, 0, 50),
)
lm_bad_adjusted <- lm(births ~ storks + area_km2 + stadistical_register , data = df)
summary(lm_bad_adjusted)
Call:
lm(formula = births ~ storks + area_km2 + stadistical_register,
data = df)
Residuals:
Min 1Q Median 3Q Max
-98.461 -30.213 -5.229 36.336 73.613
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.572e+00 1.278e+01 0.671 0.506
storks -1.220e+00 1.956e-01 -6.238 1.27e-07 ***
area_km2 -5.674e-05 1.381e-04 -0.411 0.683
stadistical_register 1.000e+00 1.332e-03 750.951 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 42.4 on 46 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 7.974e+05 on 3 and 46 DF, p-value: < 2.2e-16Imagina un mundo donde cuando ves un nido con dos cigueñas, formando un corazon con sus cuellos, decides enviar una foto a tu persona favorita y eso termina en bebe.
df1<-df %>%
mutate(
photo =storks/2 + rnorm(n, 0, 5),
births = round( 2*photo + rnorm(n, 0, 50), 0))
lm_fake <- lm(births ~ storks , data = df1)
summary(lm_fake)
Call:
lm(formula = births ~ storks, data = df1)
Residuals:
Min 1Q Median 3Q Max
-128.967 -32.333 -1.495 35.897 130.033
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.9665 11.3143 0.969 0.33727
storks 0.8039 0.2304 3.489 0.00105 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 53.98 on 48 degrees of freedom
Multiple R-squared: 0.2023, Adjusted R-squared: 0.1857
F-statistic: 12.17 on 1 and 48 DF, p-value: 0.001049lm <- lm(births ~ storks +photo , data = df1)
summary(lm)
Call:
lm(formula = births ~ storks + photo, data = df1)
Residuals:
Min 1Q Median 3Q Max
-127.589 -30.770 -1.309 35.759 130.115
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.5822 11.7676 0.899 0.373
storks 0.7048 0.7567 0.931 0.356
photo 0.2116 1.5372 0.138 0.891
Residual standard error: 54.54 on 47 degrees of freedom
Multiple R-squared: 0.2026, Adjusted R-squared: 0.1687
F-statistic: 5.971 on 2 and 47 DF, p-value: 0.004888lm_fake <- lm(births ~ storks , data = df1)
summary(lm_fake)
Call:
lm(formula = births ~ storks, data = df1)
Residuals:
Min 1Q Median 3Q Max
-128.967 -32.333 -1.495 35.897 130.033
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.9665 11.3143 0.969 0.33727
storks 0.8039 0.2304 3.489 0.00105 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 53.98 on 48 degrees of freedom
Multiple R-squared: 0.2023, Adjusted R-squared: 0.1857
F-statistic: 12.17 on 1 and 48 DF, p-value: 0.001049lm <- lm(births ~ storks +photo , data = df1)
summary(lm)
Call:
lm(formula = births ~ storks + photo, data = df1)
Residuals:
Min 1Q Median 3Q Max
-127.589 -30.770 -1.309 35.759 130.115
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.5822 11.7676 0.899 0.373
storks 0.7048 0.7567 0.931 0.356
photo 0.2116 1.5372 0.138 0.891
Residual standard error: 54.54 on 47 degrees of freedom
Multiple R-squared: 0.2026, Adjusted R-squared: 0.1687
F-statistic: 5.971 on 2 and 47 DF, p-value: 0.004888Hemos resuelto el misterio porque pudimos medir y controlar por area_km2. Pero, ¿y si el confusor fuera algo inobservable como la “cultura rural” o la “calidad del ecosistema”?
En ese caso, el ajuste por regresión falla. Para estos escenarios de confusión no observada, los investigadores recurren a métodos más avanzados.
::: {.callout-info title=“Adelanto de la Próxima Sesión: Variables Instrumentales (IV)”} Un instrumento es una variable Z que nos ofrece una fuente de variación “como si fuera aleatoria” en nuestro tratamiento X (cigüeñas). Para ser válido, un instrumento debe cumplir tres condiciones:
Relevancia: Z tiene un efecto causal sobre X. (Ej: Una enfermedad que solo afecta a las cigüeñas y reduce su población).
Restricción de Exclusión: Z afecta a Y (nacimientos) únicamente a través de su efecto en X. (La enfermedad de las cigüeñas no puede afectar directamente a la fertilidad humana).
Independencia: Z es independiente de cualquier confusor no observado. (La enfermedad no puede estar relacionada con la “cultura rural” del país).
Si encontramos un instrumento válido, podemos estimar el efecto causal incluso sin medir todos los confusores. :::
# --- Creación del Instrumento Fuerte y Válido ---
# Creamos una nueva variable 'proteccion_ciguenas'.
# La hacemos una función directa de 'storks' para asegurar que sea RELEVANTE.
# Añadimos un poco de ruido para que no sea una relación perfecta.
# No la relacionamos con 'births' ni 'area_km2' para asegurar que sea VÁLIDA.
df <- df %>%
mutate(
family_payment_program = round(0.05 * births + rnorm(n, mean = 0, sd = 5), 0),
stork_policy = round(0.1 * storks + rnorm(n, mean = 0, sd = 5), 0))
# Verificamos la correlación (debería ser alta)
cor(df$storks, df$stork_policy)[1] 0.6161176Para que una variable instrumental funcione, debe cumplir dos condiciones clave:
Relevancia (Fuerza): El instrumento debe estar fuertemente correlacionado con la variable problemática (la endógena). En nuestro caso, el instrumento debe predecir bien el número de cigüeñas (storks).
Exclusión (Validez): El instrumento solo puede afectar a la variable resultado (births) a través de la variable problemática (storks). No puede tener un efecto directo sobre los nacimientos ni estar relacionado con la variable de confusión (area_km2).
Tus instrumentos anteriores fallaban principalmente en la primera condición (eran débiles). El instrumento family_payment_program también fallaba en la segunda, ya que lo hiciste afectar directamente a births.
base::summary(lm(births ~ storks))
Call:
lm(formula = births ~ storks)
Residuals:
Min 1Q Median 3Q Max
-20543.8 -4949.7 -410.7 6269.3 19125.5
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13816.42 1707.73 8.091 1.61e-10 ***
storks 81.56 33.54 2.432 0.0188 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9580 on 48 degrees of freedom
Multiple R-squared: 0.1097, Adjusted R-squared: 0.09115
F-statistic: 5.914 on 1 and 48 DF, p-value: 0.0188library(AER)
# Modelo 2SLS usando el instrumento incorrecto
iv_mod_malo <- AER::ivreg(
births ~ storks|
family_payment_program , # 'family_payment_program' como instrumento
data = df
)
# El summary() con diagnostics=TRUE nos dará el F-test de instrumento débil
summary(iv_mod_malo, diagnostics = TRUE)
Call:
AER::ivreg(formula = births ~ storks | family_payment_program,
data = df)
Residuals:
Min 1Q Median 3Q Max
-58615 -17513 3622 17702 42595
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -10958.1 10672.3 -1.027 0.30967
storks 762.6 279.2 2.732 0.00879 **
Diagnostic tests:
df1 df2 statistic p-value
Weak instruments 1 48 7.47 0.00876 **
Wu-Hausman 1 47 543181.85 < 2e-16 ***
Sargan 0 NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 24000 on 48 degrees of freedom
Multiple R-Squared: -5.432, Adjusted R-squared: -5.567
Wald test: 7.462 on 1 and 48 DF, p-value: 0.008793 # --- 1ª ETAPA: Relevancia (Z -> X) ---
# Vemos si la política de multas (Z) realmente afecta a las cigüeñas (X)
primera_etapa <- lm(storks ~ stork_policy, data = df)
print("--- Resultados de la Primera Etapa (Relevancia) ---")[1] "--- Resultados de la Primera Etapa (Relevancia) ---"broom::tidy(primera_etapa)# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 27.0 4.13 6.54 0.0000000376
2 stork_policy 3.55 0.654 5.42 0.00000191 # Modelo 2SLS usando el instrumento incorrecto
iv_mod <- AER::ivreg(
births ~ storks|
stork_policy , # 'stork_policy' como instrumento
data = df
)
# El summary() con diagnostics=TRUE nos dará el F-test de instrumento débil
summary(iv_mod, diagnostics = TRUE)
Call:
AER::ivreg(formula = births ~ storks | stork_policy, data = df)
Residuals:
Min 1Q Median 3Q Max
-17551.3 -6209.7 -650.4 6015.5 19911.9
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11725.08 2560.44 4.579 3.32e-05 ***
storks 136.72 61.56 2.221 0.0311 *
Diagnostic tests:
df1 df2 statistic p-value
Weak instruments 1 48 29.370 1.91e-06 ***
Wu-Hausman 1 47 0.507 0.48
Sargan 0 NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8886 on 48 degrees of freedom
Multiple R-Squared: 0.1186, Adjusted R-squared: 0.1002
Wald test: 4.933 on 1 and 48 DF, p-value: 0.0311 SUTVA es un supuesto con dos partes que a menudo damos por sentado:
No Interferencia: El tratamiento de una unidad no afecta al resultado de otra. Violación de SUTVA: Una política de protección de cigüeñas en Francia causa una migración masiva de estas a España, afectando a los nacimientos (correlacionados) en España.
Consistencia: La definición del tratamiento es la misma y está bien definida para todas las unidades. Violación de SUTVA: Un país cuenta todas las especies de cigüeñas, mientras que otro solo cuenta las cigüeñas blancas. El “tratamiento” no es consistente.
Consejo para tu Investigación 💡: Piensa activamente en cómo SUTVA podría fallar en tu estudio. ¿Existen efectos de contagio (spillovers) entre tus unidades? ¿La intervención que estudias se implementó de la misma forma para todos? Ser explícito sobre estas potenciales limitaciones no debilita, sino que fortalece tu investigación.
df2 <- df %>%
mutate(
stork_fine_policy = rbinom(n, 1, 0.3),
storks = storks + 2* storks * stork_fine_policy+ rnorm(n, 0, 5))%>%
mutate(storks = if_else(storks < 0, 0, storks))
primera_etapa <- lm(births ~ storks, data = df2)
broom::tidy(primera_etapa)# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 13975. 1604. 8.71 1.90e-11
2 storks 43.1 16.1 2.68 9.97e- 3primera_etapa <- lm(storks ~ stork_fine_policy, data = df2)
broom::tidy(primera_etapa)# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 33.8 11.0 3.07 0.00356
2 stork_fine_policy 96.6 20.1 4.80 0.0000158# Modelo 2SLS usando el instrumento incorrecto
iv_mod <- AER::ivreg(
births ~ storks|
stork_fine_policy , # 'stork_fine_policy' como instrumento
data = df2
)
# El summary() con diagnostics=TRUE nos dará el F-test de instrumento débil
summary(iv_mod, diagnostics = TRUE)
Call:
AER::ivreg(formula = births ~ storks | stork_fine_policy, data = df2)
Residuals:
Min 1Q Median 3Q Max
-14871.7 -5764.7 -616.8 5930.5 18706.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12930.80 2189.76 5.905 3.5e-07 ***
storks 59.71 28.50 2.095 0.0415 *
Diagnostic tests:
df1 df2 statistic p-value
Weak instruments 1 48 23.05 1.58e-05 ***
Wu-Hausman 1 47 0.51 0.479
Sargan 0 NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8924 on 48 degrees of freedom
Multiple R-Squared: 0.111, Adjusted R-squared: 0.09246
Wald test: 4.389 on 1 and 48 DF, p-value: 0.04146 # 1ª etapa (relevancia)
fs <- lm(storks ~ stork_policy, data = df)
broom::tidy(fs)# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 27.0 4.13 6.54 0.0000000376
2 stork_policy 3.55 0.654 5.42 0.00000191 # F parcial aprox. del instrumento
fs_anova <- anova(fs)
X1 <- lm(storks ~ family_payment_program, data = df)
rss_r <- sum(resid(X1)^2); rss_u <- sum(resid(fs)^2)
k <- 1; nobs <- nrow(df); p <- length(coef(fs))
F_partial <- ((rss_r - rss_u)/k) / (rss_u/(nobs - p))
tibble(F_partial = F_partial)# A tibble: 1 × 1
F_partial
<dbl>
1 19.0# 2SLS
iv_mod <- AER::ivreg(
births ~ storks |
stork_policy,
data = df
)
summary(iv_mod, diagnostics = TRUE)
Call:
AER::ivreg(formula = births ~ storks | stork_policy, data = df)
Residuals:
Min 1Q Median 3Q Max
-17551.3 -6209.7 -650.4 6015.5 19911.9
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11725.08 2560.44 4.579 3.32e-05 ***
storks 136.72 61.56 2.221 0.0311 *
Diagnostic tests:
df1 df2 statistic p-value
Weak instruments 1 48 29.370 1.91e-06 ***
Wu-Hausman 1 47 0.507 0.48
Sargan 0 NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8886 on 48 degrees of freedom
Multiple R-Squared: 0.1186, Adjusted R-squared: 0.1002
Wald test: 4.933 on 1 and 48 DF, p-value: 0.0311 # Visualización
df$storks_hat <- fitted(fs)
ggplot(df, aes(storks_hat, births)) + geom_point(alpha=.6) +
geom_smooth(method="lm", se=FALSE) +
labs(title="2ª etapa (intuición): births vs storks_hat")
# install.packages("haven") # Descomenta y ejecuta si no lo tienes instalado
library(haven) # Paquete para leer y escribir formatos de Stata, SPSS, etc.
# --- PASO 1: Define la ruta donde quieres guardar el archivo ---
# Reemplaza "C:/Users/TuUsuario/Documents/" con la ruta real de tu ordenador.
# Consejo: Usa la barra inclinada "/" en lugar de "\" en las rutas de R.
ruta <- "C:/Users/Usuario/Desktop/curso us/datos_paises.dta"
# --- PASO 2: Usa la función write_dta() para guardar los datos ---
# El primer argumento es el data.frame que quieres guardar.
# El segundo argumento es la ruta completa con el nombre del archivo.
write_dta(df, ruta)
# Mensaje de confirmación (opcional, para saber que se ejecutó)
print(paste("Archivo guardado exitosamente en:", ruta))[1] "Archivo guardado exitosamente en: C:/Users/Usuario/Desktop/curso us/datos_paises.dta"# --- PASO 1: Usa la función read_dta() para cargar los datos ---
# La función lee el archivo y lo carga en un nuevo objeto en R.
# Lo llamaremos 'df_cargado' para diferenciarlo del original.
df_cargado <- read_dta(ruta)
# --- PASO 2: Verifica que los datos se han cargado correctamente ---
# Mostramos las primeras filas del nuevo data.frame.
head(df_cargado)# A tibble: 6 × 16
country storks births area_km2 gdp_per_capita high_storks red blue white
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
1 Country 1 121 22706 219031. 2724. Muchas Cigü… 0 1 0
2 Country 2 0 31637 270770. 41352. Pocas Cigüe… 0 0 0
3 Country 3 0 17279 176817. 33872. Pocas Cigüe… 1 0 0
4 Country 4 0 6661 63093. 14869. Pocas Cigüe… 0 0 0
5 Country 5 0 21249 210653. 21815. Pocas Cigüe… 1 1 0
6 Country 6 22 8786 86203. 16141. Pocas Cigüe… 1 0 0
# ℹ 7 more variables: green <dbl>, yellow <dbl>, black <dbl>,
# stadistical_register <dbl>, family_payment_program <dbl>,
# stork_policy <dbl>, storks_hat <dbl>Nuestro objetivo es saber si tener un número “alto” de cigüeñas (nuestro “tratamiento”) tiene un efecto causal en el número de nacimientos.
El problema es que los países con muchas cigüeñas son diferentes a los que tienen pocas (por ejemplo, son más grandes). Compararlos directamente sería como comparar manzanas con naranjas y concluir que las manzanas engordan más sin tener en cuenta otros factores. El PSM nos ayuda a encontrar una naranja lo más parecida posible a cada manzana para que la comparación sea justa.
library(MatchIt) # Para hacer el matching
library(cobalt) # Para visualizar el balance
# --- Crear la variable de tratamiento binario ---
# Si ciguenas > mediana, el país es "tratado" (1), si no, es "control" (0).
df <- df %>%
mutate(t_storks = as.integer(storks > median(storks)))El Propensity Score es la probabilidad de que un país reciba el “tratamiento” (tener muchas cigüeñas), dadas sus características (como el área ). Lo calculamos con una regresión logística.
# Modelo logístico para predecir la probabilidad de ser tratado
ps_modelo <- glm(t_storks ~ area_km2+gdp_per_capita,
data = df,
family = binomial())
# Añadimos el propensity score (la probabilidad predicha) a nuestros datos
df$ps <- predict(ps_modelo, type = "response")
# Vemos los primeros propensity scores
head(select(df, country, t_storks, ps),10)# A tibble: 10 × 3
country t_storks ps
<chr> <int> <dbl>
1 Country 1 1 0.786
2 Country 2 0 0.531
3 Country 3 0 0.458
4 Country 4 0 0.466
5 Country 5 0 0.625
6 Country 6 0 0.490
7 Country 7 0 0.160
8 Country 8 1 0.521
9 Country 9 0 0.162
10 Country 10 1 0.757Antes de hacer nada, veamos si los grupos son realmente diferentes. Este gráfico muestra la distribución de los propensity scores. Si las curvas están muy separadas, significa que los grupos son muy distintos.
# Gráfico de densidad para ver el solapamiento de los propensity scores
ggplot(df, aes(x = ps, fill = factor(t_storks))) +
geom_density(alpha = 0.5) +
labs(title = "Distribución de Propensity Scores antes del Matching",
x = "Probabilidad de ser tratado",
fill = "Grupo")
Como vemos, el grupo “tratado” (1) tiende a tener una probabilidad más alta, lo que confirma nuestro sesgo de selección.
Ahora usamos el paquete MatchIt para encontrar un “gemelo” (el país más parecido en su propensity score) para cada país tratado dentro del grupo de control.
# Realizamos el matching 1 a 1 usando el vecino más cercano ("nearest")
match_resultado <- matchit(t_storks ~ area_km2 +gdp_per_capita,
data = df,
method = "nearest",
distance = "logit", # Usa el propensity score de un modelo logit
ratio = 1) # 1 control por cada tratado
# Vemos un resumen del resultado del matching
summary(match_resultado)
Call:
matchit(formula = t_storks ~ area_km2 + gdp_per_capita, data = df,
method = "nearest", distance = "logit", ratio = 1)
Summary of Balance for All Data:
Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance 0.5493 0.4507 0.6547 0.9314 0.1752
area_km2 177303.0940 134113.9084 0.5050 0.8733 0.1288
gdp_per_capita 24143.1902 28630.6130 -0.3773 1.0016 0.0952
eCDF Max
distance 0.32
area_km2 0.28
gdp_per_capita 0.20
Summary of Balance for Matched Data:
Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
distance 0.5493 0.4507 0.6547 0.9314 0.1752
area_km2 177303.0940 134113.9084 0.5050 0.8733 0.1288
gdp_per_capita 24143.1902 28630.6130 -0.3773 1.0016 0.0952
eCDF Max Std. Pair Dist.
distance 0.32 0.6547
area_km2 0.28 1.0408
gdp_per_capita 0.20 1.0378
Sample Sizes:
Control Treated
All 25 25
Matched 25 25
Unmatched 0 0
Discarded 0 0l resumen nos dice cuántas unidades se han emparejado y cuántas se han descartado.
¿Funcionó el matching? ¿Son ahora nuestros grupos “tratado” y “control” comparables? El Love Plot es la mejor herramienta para verificarlo.
# Creamos un Love Plot para visualizar el balance de las variables
love.plot(match_resultado, thresholds = c(m = 0.1))
Interpretación del Love Plot:
Puntos rojos (All): Es el desequilibrio antes del matching. Vemos que area_km2 estaban muy desbalanceados (lejos de la línea del cero).
Puntos azules (Matched): Es el desequilibrio después del matching. Vemos que ahora están muy cerca de la línea del cero.
Conclusión: ¡Éxito! Hemos creado dos grupos que ahora son muy similares en sus características.
Ahora que tenemos grupos comparables, podemos finalmente estimar el efecto del tratamiento de una forma justa. Para ello, extraemos los datos emparejados y comparamos la media de nacimientos entre los grupos.
# 1. Extraemos el dataset con los datos emparejados
datos_match <- match.data(match_resultado)
# 2. Comparamos la media de nacimientos entre los grupos con un t-test
# Esta es la estimación del "Average Treatment Effect on the Treated" (ATT)
t.test(births ~ t_storks, data = datos_match)
Welch Two Sample t-test
data: births by t_storks
t = -1.7756, df = 47.886, p-value = 0.08215
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
-9818.9816 609.7016
sample estimates:
mean in group 0 mean in group 1
14377.64 18982.28 interpretación del resultado final: El t.test te dará la diferencia de medias entre los grupos y un p-valor. Si el p-valor es alto (ej. > 0.05), concluiremos que, una vez que comparamos grupos similares, no hay un efecto estadísticamente significativo de tener muchas cigüeñas sobre el número de nacimientos. ¡Hemos corregido la correlación espuria!
Pregunta de investigación: En 2025 se introduce una nueva ley de basuras (waste_law) que afecta a las poblaciones de cigüeñas. Queremos saber si esta ley tuvo algún efecto secundario en el número de nacimientos (births).
El escenario simulado:
Tratamiento: 30 países aprueban la waste_law en 2025. 20 países no la aprueban.
Efecto sobre las cigüeñas (storks): Las poblaciones de cigüeñas disminuyen un 3% anual de forma natural. En los países con la ley, esta caída se acelera al 10% anual.
Efecto sobre los nacimientos (births): Sabemos que la ley no tiene efecto en los nacimientos. Sin embargo, en todos los países, los nacimientos crecen un 2% anual por una tendencia demográfica general.
El desafío: ¿Puede el método DiD llegar a la conclusión correcta de que el efecto de la ley sobre los nacimientos es cero, a pesar de que las poblaciones de cigüeñas y bebés están cambiando?
Construiremos una base de datos de panel con 50 países durante 3 años (2024-2026), reflejando exactamente las condiciones del problema.
# --- Cargar librerías ---
library(tidyverse)
# --- Configuración inicial ---
set.seed(2025)
N_TOTAL <- 50
N_TRATADOS <- 30
N_CONTROL <- 20
AÑO_LEY <- 2025
# --- 1. Crear los datos base para el año 2024 ---
area_km2 <- runif(N_TOTAL, 500, 300000)
storks <- round(area_km2 / 5000 + rnorm(N_TOTAL, 0, 50), 0)
births <- round(0.1 * area_km2 + rnorm(N_TOTAL, 0, 5000), 0)
gdp_per_capita <- runif(N_TOTAL, 2000, 50000)
df_base_2024 <- tibble(
country_id = 1:N_TOTAL,
storks, births, area_km2, gdp_per_capita
)
# --- 2. Crear la estructura del panel de datos para 3 años ---
# Asignamos el tratamiento y expandimos para tener una fila por país y año
panel_df <- tibble(
country_id = 1:N_TOTAL,
waste_law = sample(c(rep(1, N_TRATADOS), rep(0, N_CONTROL)))
) %>%
crossing(year = 2024:2026) %>%
# Unimos los datos base de 2024
left_join(df_base_2024, by = "country_id") %>%
# Definimos la variable post-tratamiento
mutate(post_law = ifelse(year >= AÑO_LEY, 1, 0))
# --- 3. Aplicar las tendencias para 2025 y 2026 ---
datos_did <- panel_df %>%
group_by(country_id) %>%
mutate(
# Aplicamos las tendencias de crecimiento/decrecimiento
storks = case_when(
year == 2024 ~ storks,
year > 2024 & waste_law == 1 ~ storks * (1 - 0.10)^(year - 2024), # Caída del 10% en tratados
year > 2024 & waste_law == 0 ~ storks * (1 - 0.03)^(year - 2024) # Caída del 3% en control
),
births = case_when(
year == 2024 ~ births,
year > 2024 ~ births * (1 + 0.02)^(year - 2024) # Crecimiento del 2% para TODOS
)
) %>%
ungroup() %>%
mutate(across(c(storks, births), round)) # Redondeamos los resultados
datos_did<- datos_did %>%
mutate(
storks=if_else(storks<0,0,storks),
births=if_else(births<0,0,births)
)
# Vemos las primeras filas de nuestros datos finales
head(datos_did)# A tibble: 6 × 8
country_id waste_law year storks births area_km2 gdp_per_capita post_law
<int> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 0 2024 0 26979 219920. 25138. 0
2 1 0 2025 0 27519 219920. 25138. 1
3 1 0 2026 0 28069 219920. 25138. 1
4 2 1 2024 0 11573 142991. 47585. 0
5 2 1 2025 0 11804 142991. 47585. 1
6 2 1 2026 0 12041 142991. 47585. 1ruta <- "C:/Users/Usuario/Desktop/curso us/datos_did.dta"
write_dta(datos_did, ruta)datos_did <- read_dta(ruta)Un gráfico es la mejor forma de ver si la “suposición de tendencias paralelas” se cumple y de intuir el resultado.
ggplot(datos_did, aes(x = year, y = births, color = factor(waste_law), group = factor(waste_law))) +
# Usamos stat_summary para graficar la media de cada grupo
stat_summary(fun = mean, geom = "line", linewidth = 1.2) +
stat_summary(fun = mean, geom = "point", size = 3) +
# Línea vertical para marcar el inicio del tratamiento
geom_vline(xintercept = AÑO_LEY - 0.5, linetype = "dashed", color = "grey40") +
annotate("text", x = 2024.4, y = max(datos_did$births), label = "Pre-Ley", hjust = 1) +
annotate("text", x = 2025.4, y = max(datos_did$births), label = "Post-Ley", hjust = 0) +
scale_x_continuous(breaks = 2024:2026) +
labs(
title = "Evolución de los Nacimientos: Grupos Tratado vs. Control",
subtitle = "Las líneas permanecen paralelas, sugiriendo un efecto nulo de la ley",
x = "Año",
y = "Media de Nacimientos",
color = "Grupo (1 = Ley, 0 = Sin Ley)"
) +
theme_minimal()
Interpretación del Gráfico: Las líneas son prácticamente paralelas. Aunque los países con la ley (waste_law = 1) tienen un número diferente de nacimientos, la trayectoria de crecimiento es idéntica a la de los países sin la ley. Esto es una fuerte evidencia visual de que la ley no tuvo ningún efecto adicional sobre los nacimientos.
Finalmente, confirmamos nuestra intuición visual con el modelo de regresión. El coeficiente del término de interacción (waste_law:post_law) nos dará la estimación del efecto causal.
# La variable de interés es la interacción entre el grupo y el tiempo
modelo_did <- lm(births ~ waste_law * post_law, data = datos_did)
summary(modelo_did)
Call:
lm(formula = births ~ waste_law * post_law, data = datos_did)
Residuals:
Min 1Q Median 3Q Max
-19281 -7288 1681 6800 18226
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 18716.0 2167.2 8.636 9.25e-15 ***
waste_law -5092.0 2797.8 -1.820 0.0708 .
post_law 565.3 2654.3 0.213 0.8317
waste_law:post_law -153.8 3426.6 -0.045 0.9643
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9692 on 146 degrees of freedom
Multiple R-squared: 0.06663, Adjusted R-squared: 0.04745
F-statistic: 3.474 on 3 and 146 DF, p-value: 0.01772Estos son los tres supuestos fundamentales que has tomado implícitamente al usar Propensity Score Matching:
(Ignorability / Unconfoundedness)
¿Qué significa en teoría?: Una vez que controlamos por las características observables (las variables que usaste para calcular el propensity score), la asignación al tratamiento es “como si fuera aleatoria”. Dicho de otro modo, no hay ninguna variable “oculta” o no medida que afecte tanto a la probabilidad de recibir el tratamiento como al resultado.
¿Qué significa en nuestro ejemplo?: Estás asumiendo que area_km2 y gdp_per_capita son todos los factores importantes que influyen simultáneamente en que un país tenga muchas cigüeñas y en su número de nacimientos. Estás afirmando que no existe un confusor “fantasma” (como, por ejemplo, la “cultura rural” o la “calidad del ecosistema local”) que no hayas medido y que esté afectando a ambas variables.
¿Es un supuesto verificable?: No. Esta es la suposición más fuerte y no se puede probar con los datos. Es una afirmación teórica que debes justificar con tu conocimiento del tema. Es el “talón de Aquiles” de este método.
(Common Support / Overlap)
¿Qué significa en teoría?: Para cualquier combinación de características, existe una probabilidad mayor que cero de ser tanto tratado como no tratado. Nadie tiene una probabilidad de 0 o 1 de recibir el tratamiento.
¿Qué significa en nuestro ejemplo?: Estás asumiendo que existen tanto países con muchas cigüeñas como con pocas en todos los niveles de tamaño y PIB. No puede haber una situación en la que, por ejemplo, todos los países grandes (area_km2 > 200,000) tengan muchas cigüeñas y todos los pequeños tengan pocas. Si fuera así, los países grandes no tendrían con quién compararse (no habría un “gemelo” en el grupo de control) y el matching sería imposible para ellos.
¿Es un supuesto verificable?: Sí. El “Paso 2: Diagnóstico Visual del Solapamiento” que ya hiciste es exactamente la forma de comprobarlo. Tu gráfico de densidad muestra que las curvas de los propensity scores se solapan, lo que indica que el supuesto se cumple. Si las curvas estuvieran completamente separadas, tendrías un problema de soporte común.
Este es un supuesto fundamental en casi toda la inferencia causal y tiene dos partes:
a) No Interferencia: El tratamiento de una unidad (un país) no afecta al resultado de otra unidad.
b) Consistencia: La definición del “tratamiento” es la misma y está bien definida para todas las unidades.
t_storks = 1) es homogénea. No hay diferentes “versiones” de tener muchas cigüeñas que podrían tener efectos distintos sobre los nacimientos.¿Es un supuesto verificable?: No. Al igual que la ignorabilidad, SUTVA se debe justificar teóricamente según el diseño de tu estudio.
Pregunta de investigación: Un gobierno decide otorgar la prestigiosa (y bien financiada) designación de “Santuario de Cigüeñas” a los países que superen un índice de calidad medioambiental de 50 puntos. Queremos saber si recibir esta designación (el “tratamiento”) tiene algún efecto causal en la tasa de natalidad.
El problema es que los países con un índice ambiental alto (ej. 70) son muy diferentes a los que lo tienen bajo (ej. 30). No podemos simplemente compararlos.
La idea de RDD es simple y brillante: en lugar de comparar países muy diferentes, comparamos a los que están justo a ambos lados del punto de corte.
Un país con un índice de 49.9 es prácticamente idéntico a uno con 50.1 en todas sus características, excepto en una: el segundo recibió la designación de “Santuario” y el primero no. La asignación al tratamiento en este punto es casi aleatoria. Por lo tanto, si vemos un “salto” o discontinuidad en la tasa de natalidad justo en el umbral de 50, podemos atribuir ese salto al efecto causal del tratamiento.
Simularemos datos para 100 países. El número de nacimientos aumentará naturalmente con el índice ambiental (países más sanos, más nacimientos), pero sabemos que la designación de “Santuario” no tiene ningún efecto adicional.
# --- Simular datos para RDD ---
set.seed(2025)
datos_rdd <- tibble(
pais_id = 1:100,
# 1. La variable continua que determina el tratamiento ("running variable")
indice_ambiental = runif(100, 0, 100)
) %>%
# 2. El punto de corte estricto para recibir el tratamiento
mutate(
santuario_ciguenas = ifelse(indice_ambiental >= 50, 1, 0),
# 3. El resultado (nacimientos)
# Es una función del índice, pero SIN EFECTO del tratamiento.
# El +1000 solo se aplica si eres tratado, pero es para que se vea la discontinuidad
# En este caso, el efecto es CERO
nacimientos = 5000 + 50 * indice_ambiental + rnorm(100, 0, 500)
) %>%
arrange(indice_ambiental)
head(datos_rdd)# A tibble: 6 × 4
pais_id indice_ambiental santuario_ciguenas nacimientos
<int> <dbl> <dbl> <dbl>
1 14 2.38 0 5806.
2 23 3.96 0 4787.
3 89 6.31 0 4796.
4 75 6.92 0 5471.
5 84 8.57 0 4753.
6 37 9.46 0 6125.El gráfico es la herramienta más importante en un análisis RDD. Trazamos la variable de resultado (nacimientos) contra la variable de asignación (indice_ambiental) y marcamos el punto de corte.
ggplot(datos_rdd, aes(x = indice_ambiental, y = nacimientos, color = factor(santuario_ciguenas))) +
geom_point(alpha = 0.6) +
# Línea vertical en el punto de corte
geom_vline(xintercept = 50, linetype = "dashed", color = "black", linewidth = 1) +
# Añadimos líneas de tendencia a cada lado del corte
geom_smooth(data = filter(datos_rdd, indice_ambiental < 50), method = "lm", se = FALSE) +
geom_smooth(data = filter(datos_rdd, indice_ambiental >= 50), method = "lm", se = FALSE) +
labs(
title = "Análisis de Regresión por Discontinuidad",
subtitle = "No hay un 'salto' visible en el punto de corte de 50",
x = "Índice de Calidad Medioambiental",
y = "Número de Nacimientos",
color = "Recibió 'Santuario' (1=Sí, 0=No)"
) +
theme_minimal()
Interpretación del Gráfico: Observa el comportamiento de los puntos y las líneas de tendencia justo alrededor de la línea discontinua en 50. No hay un “salto” (o discontinuidad). La línea de tendencia de la izquierda parece conectar de forma natural con la de la derecha. Esto es una fuerte evidencia visual de que el tratamiento (la designación de “Santuario”) no tuvo ningún efecto sobre los nacimientos.
Para cuantificar el “salto”, podemos usar un modelo de regresión. La forma más simple es centrarse en los datos que están en una “ventana” o “ancho de banda” (bandwidth) cerca del punto de corte y ajustar un modelo lineal.
El coeficiente de la variable de tratamiento (santuario_ciguenas) nos dará la magnitud del salto.
# Nos centramos en una ventana, por ejemplo, +/- 15 puntos alrededor del corte
datos_ventana <- datos_rdd %>%
filter(indice_ambiental >= 35 & indice_ambiental <= 65)
# Ajustamos un modelo lineal simple en esta ventana
modelo_rdd <- lm(nacimientos ~ indice_ambiental + santuario_ciguenas, data = datos_ventana)
summary(modelo_rdd)
Call:
lm(formula = nacimientos ~ indice_ambiental + santuario_ciguenas,
data = datos_ventana)
Residuals:
Min 1Q Median 3Q Max
-635.9 -314.4 -103.1 270.3 694.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4261.61 705.82 6.038 2.23e-06 ***
indice_ambiental 61.99 16.44 3.771 0.000848 ***
santuario_ciguenas 86.73 278.49 0.311 0.757969
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 418.3 on 26 degrees of freedom
Multiple R-squared: 0.6675, Adjusted R-squared: 0.6419
F-statistic: 26.1 on 2 and 26 DF, p-value: 6.073e-07Interpretación del resultado de la regresión:
(Intercept) y indice_ambiental: Capturan la tendencia de la relación entre el índice y los nacimientos.
santuario_ciguenas: Este es nuestro estimador RDD. El coeficiente será un número muy cercano a cero y su p-valor será muy alto. Esto confirma formalmente lo que vimos en el gráfico: no hay un salto estadísticamente significativo en el punto de corte.
El método RDD nos permitió estimar el efecto causal de la designación de “Santuario” de una manera muy creíble. Al comparar países que eran prácticamente idénticos justo a ambos lados del umbral de elegibilidad, pudimos concluir correctamente que la política no tuvo ningún efecto en los nacimientos. RDD es una herramienta poderosa cuando el tratamiento se asigna en base a un punto de corte estricto.
Aquí están los supuestos que has tomado:
(Continuity of Potential Outcomes)
¿Qué significa en teoría?: La relación entre la variable de asignación (el índice ambiental) y el resultado (los nacimientos) sería una función continua y suave en el punto de corte si el tratamiento no existiera.
¿Qué significa en nuestro ejemplo?: Estás asumiendo que no existe ninguna otra razón, aparte de la designación de “Santuario”, por la que la tasa de natalidad saltaría bruscamente justo en el valor 50 del índice ambiental. La relación entre la calidad ambiental y los nacimientos es suave.
¿Cómo podemos (parcialmente) verificarlo?: Aunque el supuesto en sí no es verificable, podemos ganar confianza en él. Hacemos “pruebas de placebo” verificando que otras variables pre-tratamiento (como el PIB, la población, etc.) no presenten un salto en el punto de corte. Si nada más salta en el umbral, es más plausible que el salto en nuestro resultado se deba al tratamiento.
(No Precise Manipulation of the Running Variable)
¿Qué significa en teoría?: Las unidades (países) no pueden manipular con precisión su posición en la variable de asignación para caer justo por encima o por debajo del umbral y así seleccionar su estado de tratamiento.
¿Qué significa en nuestro ejemplo?: Estás asumiendo que los países no pueden manipular con precisión su indice_ambiental. Un país con un índice real de 49.8 no puede “hacer trampas” fácilmente para reportar un 50.1 y así conseguir la financiación del santuario. Si pudieran hacerlo, los países justo a la derecha del corte serían sistemáticamente diferentes (más astutos, con más recursos para manipular, etc.) a los que están justo a la izquierda, y la comparación dejaría de ser válida.
¿Cómo podemos (parcialmente) verificarlo?: Sí. La forma estándar es hacer un gráfico de densidad de la variable de asignación (indice_ambiental). Si hay una “agrupación” sospechosa de países justo a la derecha del 50 y un “vacío” justo a la izquierda, es una fuerte señal de manipulación. A esto se le conoce como el Test de McCrary.
(Sharp Assignment)
¿Qué significa en teoría?: La probabilidad de recibir el tratamiento salta de 0 a 1 exactamente en el punto de corte.
¿Qué significa en nuestro ejemplo?: Estás asumiendo que todos los países con un índice de 50 o más recibieron la designación de “Santuario” y ninguno por debajo de 50 la recibió. No hay excepciones, incumplimientos ni errores en la asignación. Tu código santuario_ciguenas = ifelse(indice_ambiental >= 50, 1, 0) implementa perfectamente este supuesto.
¿Cómo podemos (parcialmente) verificarlo?: Sí, por definición. Simplemente se comprueba en los datos que la regla de asignación se cumple a la perfección. La alternativa es un RDD “difuso” (Fuzzy RDD), donde cruzar el umbral solo aumenta la probabilidad de ser tratado (ej. te hace “elegible” pero no garantiza el tratamiento).
Has aprendido DiD, que es perfecto para comparar un grupo de tratados con un grupo de controles. Pero, ¿qué pasa si solo un país (“Cigüeñalandia”) implementó la ley? Compararlo con el promedio de otros 49 países puede ser engañoso.
El Método de Control Sintético resuelve esto. En lugar de usar una media simple, crea un “clon” o contrafactual sintético de Cigüeñalandia usando una combinación ponderada de los países de control. Por ejemplo, podría determinar que la mejor réplica de Cigüeñalandia antes de la ley es “30% Normalia + 50% Ruritania + 20% Freedonia”. Luego, compara la trayectoria real de Cigüeñalandia después de la ley con la de su clon sintético. Es la herramienta de elección para estudios de caso causales.
La Correlación NO Implica Causalidad: Es el mantra que debe guiar todo análisis de datos. Siempre busca explicaciones alternativas.
La Teoría Guía al Análisis: Usa DAGs para mapear tus supuestos sobre cómo funciona el mundo. La inferencia causal es pensar antes de computar.
Ajustar es Cerrar Puertas Traseras: La regresión múltiple es una herramienta poderosa para eliminar la confusión, pero su validez depende críticamente del supuesto de ignorabilidad.
Conoce tus Supuestos: La validez de tus conclusiones depende de la plausibilidad de los supuestos de ignorabilidad, solapamiento y SUTVA en el contexto de tu investigación.
El objetivo de la ciencia de datos moderna no es solo predecir, sino entender las palancas causales del mundo. ¡Espero que esta sesión les haya proporcionado las herramientas conceptuales para empezar este emocionante viaje!
Hacer un modelo es solo el principio. La credibilidad de tu investigación depende de la transparencia y el rigor con que reportes tus hallazgos. Un buen análisis causal es una narrativa argumentada. Esta es tu guía para contar esa historia de forma convincente.
La Trama (Hipótesis y DAGs): Empieza siempre con la teoría.
Mecanismo Conceptual: Explica con palabras el mecanismo causal que propones. ¿Por qué crees que X debería causar Y?
DAG Explícito: Dibuja tu Diagrama Acíclico Dirigido (DAG). Este no es un adorno, es el mapa de tus supuestos causales. Obliga a pensar en los confusores, mediadores y colisionadores antes de tocar el código. Es tu contrato con el lector.
La Escena del Crimen (Datos y Limpieza):
El Arma (Estrategia de Identificación):
Declara tu Estrategia: Sé explícito. Di: “Mi estrategia de identificación se basa en un diseño de Regresión por Discontinuidad…” o “Uso Variables Instrumentales para…”.
Defiende tus Supuestos: Cada método tiene un “talón de Aquiles” (un supuesto clave no verificable). Tu trabajo es argumentar por qué ese supuesto es plausible en tu contexto.
Regresión/Matching: ¿Por qué crees que has medido todos los confusores importantes (Ignorabilidad)?
IV: ¿Por qué crees que tu instrumento es válido (Restricción de Exclusión)?
DiD: ¿Por qué crees que las tendencias habrían sido paralelas? Muestra evidencia gráfica con datos pre-tratamiento si los tienes.
RDD: ¿Por qué no habría un salto en el umbral por otra razón (Continuidad)?
El Veredicto (Resultados):
Más Allá del p-valor: No te limites a decir “la relación es significativa”. Reporta el tamaño del efecto en unidades comprensibles, junto con su intervalo de confianza. Ejemplo: “La ley de basuras causó una disminución promedio de 50 cigüeñas (IC 95%: [-70, -30])”. Esto comunica la magnitud y la incertidumbre.
Visualiza: Un buen gráfico (como un plot de DiD o RDD) es a menudo más poderoso que una tabla.
El Contrainterrogatorio (Robustez y Sensibilidad):
Anticipa las Críticas: Un analista robusto se pregunta: “¿Y si mis supuestos están un poco equivocados?”.
Errores Estándar Robustos: Usa errores estándar robustos a la heterocedasticidad (como los de estimatr o sandwich) por defecto. Si tienes datos agrupados (ej. países a lo largo del tiempo), usa errores estándar clúster.
Análisis de Sensibilidad: Demuestra qué tan frágil es tu conclusión. Menciona técnicas como el E-value o el coeficiente delta de Oster para cuantificar cuánto sesgo de una variable omitida se necesitaría para anular tu efecto.
El Legado (Reproducibilidad):
Haz un Regalo a tu Futuro ‘Yo’: Tu código debe ser un manual de instrucciones claro. Coméntalo bien. Fija una semilla (set.seed()) para la aleatoriedad.
Haz un Regalo a la Ciencia: Usa herramientas como renv para crear un entorno de R reproducible. Comparte tu código y datos (si es posible) en un repositorio. Al final de tu documento, incluye siempre la salida de sessionInfo() para registrar las versiones de los paquetes.
Ser “bilingüe” en el software estadístico es una gran ventaja. Aquí tienes una guía de traducción rápida para los métodos que hemos visto.
| Método | R (Paquete::función) | Stata (comando) |
| OLS con Controles | lm(y ~ x + c1, data) |
reg y x c1 |
| OLS (Errores Robustos) | estimatr::lm_robust(y ~ x, data) |
reg y x, robust |
| Variables Instrumentales | AER::ivreg(y ~ x | z, data) |
ivregress 2sls y (x = z) |
| Propensity Score | MatchIt::matchit() |
teffects psmatch |
| Diferencias en Difs. | fixest::feols(y ~ D | id+t) |
reghdfe y D, absorb(id t) |
| Regresión Discontinuidad | rdrobust::rdrobust() |
rdrobust |
Has completado el campamento base de la inferencia causal. La montaña es alta y emocionante. Aquí tienes algunas rutas para continuar la expedición:
Causalidad y Machine Learning (Double/Debiased ML): ¿Cómo usar la potencia predictiva de algoritmos como Random Forest o Lasso para controlar por cientos de confusores de forma robusta, sin caer en el sobreajuste? Esta es una de las fronteras más activas de la econometría. (Paquete en R: DoubleML).
Modelos de Event Study Dinámicos: DiD con múltiples periodos de tiempo y tratamientos escalonados (staggered) tiene problemas que han sido descubiertos recientemente. Los nuevos estimadores (Callaway & Sant’Anna, Sun & Abraham) son ahora el estándar de la industria. (Implementados en fixest y did).
Diseños de IV y RDD Avanzados: Explora el mundo del RDD “difuso” (Fuzzy RDD), donde el umbral solo afecta la probabilidad de ser tratado, y el concepto de LATE, que nos dice a quién se aplica realmente nuestro efecto causal.
Transportabilidad y Validez Externa: ¿Bajo qué condiciones los resultados de mi estudio en España se pueden aplicar a Portugal? La teoría de los DAGs ofrece un marco formal para responder a esta pregunta.
sessionInfo()R version 4.3.0 (2023-04-21 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 26100)
Matrix products: default
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: Europe/Madrid
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] cobalt_4.6.0 MatchIt_4.7.1 haven_2.5.4 AER_1.2-10
[5] survival_3.5-5 sandwich_3.1-1 car_3.1-2 carData_3.0-5
[9] modelsummary_1.4.3 estimatr_1.0.0 lmtest_0.9-40 zoo_1.8-12
[13] performance_0.15.1 dagitty_0.3-4 ggdag_0.2.13 kableExtra_1.3.4
[17] knitr_1.44 lubridate_1.9.2 forcats_1.0.0 stringr_1.5.0
[21] dplyr_1.1.3 purrr_1.0.2 readr_2.1.4 tidyr_1.3.0
[25] tibble_3.2.1 ggplot2_4.0.0 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] gridExtra_2.3 rlang_1.1.1 magrittr_2.0.3 multcomp_1.4-28
[5] compiler_4.3.0 mgcv_1.8-42 systemfonts_1.2.2 vctrs_0.6.3
[9] rvest_1.0.3 crayon_1.5.3 pkgconfig_2.0.3 fastmap_1.1.1
[13] backports_1.4.1 labeling_0.4.3 ggraph_2.2.1 effectsize_1.0.1
[17] utf8_1.2.3 rmarkdown_2.25 tzdb_0.4.0 xfun_0.40
[21] cachem_1.0.8 jsonlite_1.8.7 highr_0.10 chk_0.10.0
[25] tweenr_2.0.3 broom_1.0.5 R6_2.5.1 tables_0.9.17
[29] stringi_1.7.12 RColorBrewer_1.1-3 boot_1.3-28.1 estimability_1.4.1
[33] Rcpp_1.0.11 parameters_0.28.2 Matrix_1.6-4 splines_4.3.0
[37] igraph_2.0.3 timechange_0.2.0 tidyselect_1.2.0 abind_1.4-5
[41] rstudioapi_0.15.0 yaml_2.3.7 viridis_0.6.4 codetools_0.2-19
[45] curl_5.0.2 lattice_0.21-8 withr_2.5.0 bayestestR_0.17.0
[49] S7_0.2.0 evaluate_0.21 polyclip_1.10-7 xml2_1.3.5
[53] texreg_1.39.4 pillar_1.9.0 checkmate_2.2.0 DT_0.33
[57] insight_1.4.2 generics_0.1.3 hms_1.1.3 scales_1.4.0
[61] xtable_1.8-4 glue_1.6.2 emmeans_1.8.9 tools_4.3.0
[65] see_0.12.0 webshot_0.5.5 mvtnorm_1.2-4 graphlayouts_1.2.2
[69] tidygraph_1.3.1 grid_4.3.0 datawizard_1.2.0 nlme_3.1-162
[73] patchwork_1.3.2 ggforce_0.4.2 Formula_1.2-5 cli_3.6.1
[77] fansi_1.0.4 viridisLite_0.4.2 svglite_2.1.3 V8_4.4.2
[81] gtable_0.3.6 digest_0.6.33 ggrepel_0.9.5 TH.data_1.1-3
[85] htmlwidgets_1.6.2 farver_2.1.1 memoise_2.0.1 htmltools_0.5.6
[89] lifecycle_1.0.3 httr_1.4.7 MASS_7.3-58.4