Introducción al Propensity Score Matching en R
Profesor Andre Chocó Cedillos
2025-09-23
1 Fundamentos Teóricos del Propensity Score Matching
El Propensity Score Matching constituye una técnica estadística fundamental para la reducción del sesgo de selección en estudios observacionales. Esta metodología permite aproximar las condiciones de un experimento controlado aleatorizado mediante el emparejamiento de unidades tratadas y de control con características observables similares.
El propensity score se define formalmente como la probabilidad condicional de recibir el tratamiento dado un conjunto de covariables observadas:
\[e(X) = P(T = 1 | X)\]
donde \(T\) representa el indicador de tratamiento y \(X\) el vector de covariables observadas.
1.1 Métodos de Estimación de Propensity Scores
La selección del método apropiado para estimar propensity scores depende de las características específicas de los datos y los objetivos del análisis. La literatura identifica múltiples enfoques, cada uno con ventajas y limitaciones particulares.
| Método | Tipo | Complejidad | Interpretabilidad | No Linealidad | Interacciones | Robustez | Alta Dimensión |
|---|---|---|---|---|---|---|---|
| Regresión Logística | Paramétrico | Baja | Excelente | Limitada | Manual | Baja | Pobre |
| Regresión Probit | Paramétrico | Baja | Excelente | Limitada | Manual | Baja | Pobre |
| Random Forest | No Paramétrico | Media | Baja | Excelente | Automática | Alta | Buena |
| Gradient Boosting | No Paramétrico | Media-Alta | Muy Baja | Excelente | Automática | Media | Buena |
| Redes Neuronales | No Paramétrico | Alta | Muy Baja | Excelente | Automática | Baja | Excelente |
| Series de Potencias | Semi-paramétrico | Media | Moderada | Buena | Manual | Media | Pobre |
| Splines | Semi-paramétrico | Media | Moderada | Buena | Manual | Media | Pobre |
| Support Vector Machines | No Paramétrico | Media | Baja | Buena | Limitada | Alta | Buena |
1.1.1 Análisis Comparativo Detallado
Los métodos paramétricos tradicionales, particularmente la regresión logística, mantienen su relevancia debido a la interpretabilidad directa de sus coeficientes y la estabilidad de sus estimaciones. La regresión probit constituye una alternativa teóricamente equivalente, basada en la función de distribución normal acumulativa.
| Método | Aspecto | Descripción |
|---|---|---|
| Regresión Logística | Ventajas | Interpretación directa, estimación rápida, estabilidad, validación extensiva |
| Desventajas | Asume linealidad, limitaciones con interacciones, sensibilidad a outliers | |
| Random Forest | Ventajas | Captura no linealidades, maneja interacciones, robustez a outliers |
| Desventajas | Interpretación limitada, riesgo de sobreajuste, sensibilidad a hiperparámetros | |
| Gradient Boosting | Ventajas | Capacidad predictiva superior, flexibilidad para patrones complejos |
| Desventajas | Propenso al sobreajuste, interpretación muy limitada, demanda computacional |
Los métodos de machine learning ofrecen capacidades superiores para detectar patrones complejos y relaciones no lineales, pero a costa de la interpretabilidad. Esta compensación resulta particularmente relevante en contextos donde la comunicación de resultados a audiencias no técnicas constituye una prioridad.
1.1.2 Criterios de Selección Metodológica
La elección del método apropiado debe fundamentarse en las características específicas del problema de investigación:
| Escenario | Método Primario | Alternativa |
|---|---|---|
| Muestra pequeña (n < 500) | Regresión Logística | Probit |
| Muestra grande (n > 5000) | Gradient Boosting | Random Forest |
| Alta dimensionalidad (p > 50) | Random Forest | Gradient Boosting |
| Relaciones lineales esperadas | Regresión Logística | Splines simples |
| Patrones complejos sospechados | Random Forest | Gradient Boosting |
| Interpretabilidad prioritaria | Regresión Logística | Probit |
| Precisión predictiva crítica | Gradient Boosting | Random Forest |
| Presencia de outliers | Random Forest | SVM |
| Tiempo limitado | Regresión Logística | Random Forest |
| Variables categóricas dominantes | Regresión Logística | Random Forest |
1.2 Métodos de Evaluación del Balance
La evaluación rigurosa del balance constituye un requisito fundamental para la validez de las estimaciones causales derivadas del propensity score matching. Los métodos estándar incluyen métricas cuantitativas y herramientas de visualización.
Las diferencias estandarizadas de medias (SMD) proporcionan una métrica normalizada para evaluar el balance:
\[SMD = \frac{\bar{X}_1 - \bar{X}_0}{\sqrt{\frac{s_1^2 + s_0^2}{2}}}\]
donde valores absolutos menores a 0.1 indican balance adecuado, mientras que valores superiores a 0.25 señalan desequilibrio sustancial.
2 Implementación Práctica: Caso de Estudio
2.1 Contexto del Estudio
Se presenta un análisis del efecto de un programa educativo complementario sobre el rendimiento académico de estudiantes de educación secundaria. El estudio utiliza datos observacionales que exhiben sesgo de selección típico, donde la participación en el programa se asocia con características socioeconómicas y académicas preexistentes.
El programa educativo consiste en tutorías adicionales, talleres de refuerzo académico y apoyo personalizado que se ofrece fuera del horario escolar regular. La participación es voluntaria, lo que genera un proceso de autoselección donde estudiantes con ciertas características (mayor motivación, mejor rendimiento previo, mayor apoyo familiar) tienen mayor probabilidad de participar.
2.1.1 Variables del Estudio
Para implementar correctamente el propensity score matching, se establecen las variables que intervienen en cada etapa del análisis:
| Etapa del Análisis | Variable Dependiente | Variables Independientes (Covariables) |
|---|---|---|
| Modelo de propensión | programa (1 = participó, 0 = no participó) | edad, log(ingreso_familiar), rendimiento_previo, genero, nivel_padres |
| Modelo de efecto del tratamiento | rendimiento_final (puntaje académico) | programa (participación en el programa) |
La distinción entre estas etapas resulta crítica para comprender la lógica del propensity score matching. En la primera etapa, el modelo de propensión estima la probabilidad de participar en el programa basándose en características observables. En la segunda etapa, tras lograr el balance mediante matching, se estima el efecto causal del programa sobre el rendimiento académico.
# Configuración de reproducibilidad
set.seed(123)
n <- 1000
# Generación de covariables
datos_estudio <- data.frame(
estudiante_id = 1:n,
edad = pmax(14, pmin(18, rnorm(n, mean = 16, sd = 1.5))),
ingreso_familiar = exp(rnorm(n, mean = 10.5, sd = 0.5)),
rendimiento_previo = pmax(0, pmin(10, rnorm(n, mean = 7.5, sd = 1.2))),
genero = sample(c("Masculino", "Femenino"), n, replace = TRUE),
nivel_padres = sample(c("Básica", "Media", "Superior"), n,
replace = TRUE, prob = c(0.4, 0.35, 0.25))
)
# Modelo de selección corregido para mejor balance inicial
logit_prob <- with(datos_estudio,
-0.8 + 0.02*edad + 0.00001*ingreso_familiar + 0.2*rendimiento_previo +
ifelse(genero == "Femenino", 0.1, 0) +
ifelse(nivel_padres == "Superior", 0.15,
ifelse(nivel_padres == "Media", 0.08, 0)))
datos_estudio$programa <- rbinom(n, 1, plogis(logit_prob))
# Variable de resultado
datos_estudio$rendimiento_final <- with(datos_estudio,
5.2 + 1.2*programa + # Efecto del tratamiento
0.7*rendimiento_previo + # Persistencia
0.00001*ingreso_familiar + # Contexto socioeconómico
0.1*edad + # Madurez
ifelse(genero == "Femenino", 0.15, 0) + # Diferencia de género
rnorm(n, 0, 0.8)) # Error aleatorio
# Descripción de variables
descripcion_variables <- data.frame(
Variable = c("edad", "ingreso_familiar", "rendimiento_previo", "genero", "nivel_padres"),
Tipo = c("Continua", "Continua", "Continua", "Categórica", "Categórica"),
Descripción = c("Edad del estudiante en años (14-18)",
"Ingreso familiar mensual en unidades monetarias",
"Rendimiento académico previo (escala 0-10)",
"Género del estudiante (Masculino/Femenino)",
"Nivel educativo de los padres (Básica/Media/Superior)"),
Rol_en_PSM = c("Confundidor", "Confundidor", "Confundidor principal",
"Confundidor", "Confundidor principal")
)
kable(descripcion_variables,
caption = "Descripción de Covariables del Modelo de Propensity Score",
col.names = c("Variable", "Tipo", "Descripción", "Rol en PSM")) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
column_spec(1, bold = TRUE, width = "20%") %>%
column_spec(3, width = "45%") %>%
row_spec(c(3,5), background = "#fff2e8") # Destacar confundidores principales| Variable | Tipo | Descripción | Rol en PSM |
|---|---|---|---|
| edad | Continua | Edad del estudiante en años (14-18) | Confundidor |
| ingreso_familiar | Continua | Ingreso familiar mensual en unidades monetarias | Confundidor |
| rendimiento_previo | Continua | Rendimiento académico previo (escala 0-10) | Confundidor principal |
| genero | Categórica | Género del estudiante (Masculino/Femenino) | Confundidor |
| nivel_padres | Categórica | Nivel educativo de los padres (Básica/Media/Superior) | Confundidor principal |
# Estadísticas descriptivas iniciales
tabla_inicial <- datos_estudio %>%
group_by(programa) %>%
summarise(
n = n(),
edad_promedio = round(mean(edad), 2),
ingreso_promedio = round(mean(ingreso_familiar), 0),
rendimiento_previo_promedio = round(mean(rendimiento_previo), 2),
proporcion_femenino = round(mean(genero == "Femenino"), 3),
proporcion_nivel_superior = round(mean(nivel_padres == "Superior"), 3),
rendimiento_final_promedio = round(mean(rendimiento_final), 2),
.groups = 'drop'
) %>%
mutate(programa = ifelse(programa == 1, "Con Programa", "Sin Programa"))
kable(tabla_inicial,
caption = "Características Descriptivas por Grupo (Antes del Matching)",
col.names = c("Grupo", "N", "Edad", "Ingreso",
"Rend. Previo", "Prop. Femenino",
"Prop. Nivel Superior", "Rend. Final")) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| Grupo | N | Edad | Ingreso | Rend. Previo | Prop. Femenino | Prop. Nivel Superior | Rend. Final |
|---|---|---|---|---|---|---|---|
| Sin Programa | 165 | 15.98 | 38385 | 7.08 | 0.473 | 0.23 | 12.24 |
| Con Programa | 835 | 16.02 | 42858 | 7.55 | 0.509 | 0.25 | 13.78 |
Los datos ahora presentan un balance inicial más equilibrado entre grupos con 835 estudiantes en el grupo de tratamiento y 165 en el grupo de control. Sin embargo, aún se observan diferencias sistemáticas moderadas que reflejan el proceso de autoselección característico de estudios observacionales.
2.2 Estimación de Propensity Scores
# Modelo logístico de propensity score
modelo_ps <- glm(programa ~ edad + log(ingreso_familiar) + rendimiento_previo +
genero + nivel_padres,
data = datos_estudio,
family = binomial)
# Extracción de propensity scores
datos_estudio$ps <- predict(modelo_ps, type = "response")
# Coeficientes del modelo
coeficientes <- tidy(modelo_ps) %>%
mutate(
estimate = round(estimate, 3),
std.error = round(std.error, 3),
p.value = ifelse(p.value < 0.001, "< 0.001", round(p.value, 3))
)
kable(coeficientes,
caption = "Coeficientes del Modelo de Propensity Score",
col.names = c("Variable", "Coeficiente", "Error Estándar",
"Estadístico t", "Valor p")) %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Variable | Coeficiente | Error Estándar | Estadístico t | Valor p |
|---|---|---|---|---|
| (Intercept) | -5.529 | 2.175 | -2.5422221 | 0.011 |
| edad | 0.021 | 0.070 | 0.2947585 | 0.768 |
| log(ingreso_familiar) | 0.392 | 0.174 | 2.2576312 | 0.024 |
| rendimiento_previo | 0.365 | 0.077 | 4.7396857 | < 0.001 |
| generoMasculino | -0.143 | 0.174 | -0.8214122 | 0.411 |
| nivel_padresMedia | 0.196 | 0.200 | 0.9785768 | 0.328 |
| nivel_padresSuperior | 0.191 | 0.223 | 0.8592151 | 0.39 |
El modelo revela que el rendimiento académico previo continúa siendo el predictor más relevante de la participación en el programa, aunque con menor magnitud que en el escenario anterior. Los coeficientes menores indican un proceso de selección menos extremo, más realista para contextos educativos típicos.
2.2.1 Evaluación del Soporte Común
# Análisis de soporte común
soporte_stats <- datos_estudio %>%
group_by(programa) %>%
summarise(
ps_minimo = round(min(ps), 3),
ps_maximo = round(max(ps), 3),
ps_media = round(mean(ps), 3),
ps_mediana = round(median(ps), 3),
.groups = 'drop'
) %>%
mutate(programa = ifelse(programa == 1, "Con Programa", "Sin Programa"))
kable(soporte_stats,
caption = "Distribución de Propensity Scores por Grupo",
col.names = c("Grupo", "Mínimo", "Máximo", "Media", "Mediana")) %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Grupo | Mínimo | Máximo | Media | Mediana |
|---|---|---|---|---|
| Sin Programa | 0.602 | 0.934 | 0.811 | 0.819 |
| Con Programa | 0.591 | 0.954 | 0.840 | 0.849 |
# Visualización de distribuciones
plot_soporte <- ggplot(datos_estudio, aes(x = ps, fill = factor(programa))) +
geom_histogram(alpha = 0.7, position = "identity", bins = 30) +
scale_fill_manual(values = c("0" = "#d73027", "1" = "#4575b4"),
labels = c("0" = "Sin Programa", "1" = "Con Programa")) +
labs(title = "Distribución de Propensity Scores por Grupo",
x = "Propensity Score",
y = "Frecuencia",
fill = "Grupo") +
theme_minimal() +
theme(legend.position = "bottom",
plot.title = element_text(hjust = 0.5))
print(plot_soporte)
Distribución de Propensity Scores por Grupo
La evaluación del soporte común muestra un overlap excelente entre las distribuciones con medias más cercanas (0.84 vs 0.811), indicando condiciones muy favorables para el matching y mayor representatividad de la muestra matched.
2.3 Implementación del Matching
# Matching con nearest neighbor y caliper
matching_resultado <- matchit(programa ~ edad + log(ingreso_familiar) +
rendimiento_previo + genero + nivel_padres,
data = datos_estudio,
method = "nearest",
ratio = 1,
caliper = 0.1,
replace = FALSE,
estimand = "ATT")
# Extracción de datos matched
datos_matched <- match.data(matching_resultado)
# Resumen del matching
resumen_matching <- data.frame(
Metrica = c("Unidades originales", "Unidades post-matching",
"Tratadas matched", "Control matched", "Tasa de matching"),
Valor = c(nrow(datos_estudio), nrow(datos_matched),
sum(datos_matched$programa == 1),
sum(datos_matched$programa == 0),
paste0(round(nrow(datos_matched)/nrow(datos_estudio)*100, 1), "%"))
)
kable(resumen_matching,
caption = "Resumen del Proceso de Matching",
col.names = c("Métrica", "Valor")) %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Métrica | Valor |
|---|---|
| Unidades originales | 1000 |
| Unidades post-matching | 326 |
| Tratadas matched | 163 |
| Control matched | 163 |
| Tasa de matching | 32.6% |
El proceso de matching logró emparejar 326 unidades (163 pares), representando una tasa de matching del 32.6%. Esta tasa superior resulta excelente y indica que el matching preserva una proporción muy alta de la muestra original, mejorando la validez externa de los resultados.
2.4 Evaluación del Balance Post-Matching
La evaluación del balance constituye el paso más crítico para validar la efectividad del matching. Se implementan tanto métricas cuantitativas como visualizaciones para evaluar la comparabilidad de los grupos.
# Evaluación del balance usando cobalt
balance_completo <- bal.tab(matching_resultado,
un = TRUE,
thresholds = c(m = 0.1))
# Mostrar resultado del balance
print(balance_completo)## Balance Measures
## Type Diff.Un Diff.Adj M.Threshold
## distance Distance 0.4539 0.0870 Balanced, <0.1
## edad Contin. 0.0274 0.0155 Balanced, <0.1
## log(ingreso_familiar) Contin. 0.2000 0.1368 Not Balanced, >0.1
## rendimiento_previo Contin. 0.4082 0.0108 Balanced, <0.1
## genero_Masculino Binary -0.0363 -0.0245 Balanced, <0.1
## nivel_padres_Básica Binary -0.0435 -0.0368 Balanced, <0.1
## nivel_padres_Media Binary 0.0236 0.0061 Balanced, <0.1
## nivel_padres_Superior Binary 0.0200 0.0307 Balanced, <0.1
##
## Balance tally for mean differences
## count
## Balanced, <0.1 7
## Not Balanced, >0.1 1
##
## Variable with the greatest mean difference
## Variable Diff.Adj M.Threshold
## log(ingreso_familiar) 0.1368 Not Balanced, >0.1
##
## Sample sizes
## Control Treated
## All 165 835
## Matched 163 163
## Unmatched 2 6722.4.1 Visualización del Balance: Love Plot
# Love Plot con configuración corregida
love_plot_resultado <- love.plot(matching_resultado,
threshold = 0.1,
binary = "std",
abs = TRUE,
var.order = "unadjusted",
line = TRUE,
title = "Balance de Covariables: Antes vs Después del Matching",
subtitle = "Diferencias Estandarizadas de Medias (SMD)",
sample.names = c("Antes del Matching", "Después del Matching"),
colors = c("#d73027", "#4575b4")) +
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 12, hjust = 0.5),
axis.text = element_text(size = 10),
axis.title = element_text(size = 11),
legend.position = "bottom",
legend.title = element_text(size = 11),
legend.text = element_text(size = 10),
panel.grid.minor = element_blank()
) +
geom_vline(xintercept = 0.1, linetype = "dashed", color = "gray50", alpha = 0.8) +
geom_vline(xintercept = -0.1, linetype = "dashed", color = "gray50", alpha = 0.8)
print(love_plot_resultado)
Love Plot: Diferencias Estandarizadas de Medias Antes y Después del Matching
El Love Plot muestra que el matching fue altamente efectivo. Las diferencias iniciales (líneas rojas) eran más moderadas que en el ejemplo anterior, y el matching las redujo prácticamente a cero (líneas azules), logrando balance excelente para todas las covariables dentro del umbral de 0.1.
2.4.2 Distribución de Propensity Scores
# Crear gráficos separados para mejor comparación
p1 <- ggplot(datos_estudio, aes(x = ps, fill = factor(programa))) +
geom_density(alpha = 0.6) +
scale_fill_manual(values = c("0" = "#d73027", "1" = "#4575b4"),
labels = c("0" = "Sin Programa", "1" = "Con Programa")) +
labs(title = "Antes del Matching",
x = "Propensity Score",
y = "Densidad",
fill = "Grupo") +
theme_minimal() +
theme(legend.position = "bottom",
plot.title = element_text(hjust = 0.5))
p2 <- ggplot(datos_matched, aes(x = ps, fill = factor(programa))) +
geom_density(alpha = 0.6) +
scale_fill_manual(values = c("0" = "#d73027", "1" = "#4575b4"),
labels = c("0" = "Sin Programa", "1" = "Con Programa")) +
labs(title = "Después del Matching",
x = "Propensity Score",
y = "Densidad",
fill = "Grupo") +
theme_minimal() +
theme(legend.position = "bottom",
plot.title = element_text(hjust = 0.5))
# Combinar los gráficos
distribucion_combinada <- grid.arrange(p1, p2, ncol = 2,
top = textGrob("Comparación de Distribuciones de Propensity Scores",
gp = gpar(fontsize = 14, fontface = "bold")))
Distribución de Propensity Scores antes y después del matching
## TableGrob (2 x 2) "arrange": 3 grobs
## z cells name grob
## 1 1 (2-2,1-1) arrange gtable[layout]
## 2 2 (2-2,2-2) arrange gtable[layout]
## 3 3 (1-1,1-2) arrange text[GRID.text.212]Las distribuciones muestran que el soporte común era ya muy bueno antes del matching, y después del matching las distribuciones se superponen casi perfectamente, indicando comparabilidad óptima entre grupos.
2.4.3 Tabla de Balance Numérico
# Tabla de balance detallada
balance_tabla <- datos_matched %>%
group_by(programa) %>%
summarise(
n = n(),
edad = round(mean(edad), 2),
ingreso = round(mean(ingreso_familiar), 0),
rendimiento_previo = round(mean(rendimiento_previo), 2),
prop_femenino = round(mean(genero == "Femenino"), 3),
prop_superior = round(mean(nivel_padres == "Superior"), 3),
.groups = 'drop'
) %>%
mutate(programa = ifelse(programa == 1, "Con Programa", "Sin Programa"))
kable(balance_tabla,
caption = "Balance de Covariables Post-Matching",
col.names = c("Grupo", "N", "Edad", "Ingreso",
"Rendimiento Previo", "Prop. Femenino", "Prop. Nivel Superior")) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| Grupo | N | Edad | Ingreso | Rendimiento Previo | Prop. Femenino | Prop. Nivel Superior |
|---|---|---|---|---|---|---|
| Sin Programa | 163 | 15.98 | 38712 | 7.09 | 0.479 | 0.233 |
| Con Programa | 163 | 16.00 | 40934 | 7.10 | 0.503 | 0.264 |
# Calcular SMD manualmente para verificación
smd_verificacion <- datos_matched %>%
group_by(programa) %>%
summarise(
edad_mean = mean(edad), edad_var = var(edad),
ingreso_mean = mean(log(ingreso_familiar)), ingreso_var = var(log(ingreso_familiar)),
rendimiento_mean = mean(rendimiento_previo), rendimiento_var = var(rendimiento_previo),
.groups = 'drop'
)
calcular_smd <- function(mean_t, mean_c, var_t, var_c) {
(mean_t - mean_c) / sqrt((var_t + var_c) / 2)
}
smd_tabla <- data.frame(
Variable = c("Edad", "Log(Ingreso)", "Rendimiento Previo"),
SMD_Post = round(c(
calcular_smd(smd_verificacion$edad_mean[2], smd_verificacion$edad_mean[1],
smd_verificacion$edad_var[2], smd_verificacion$edad_var[1]),
calcular_smd(smd_verificacion$ingreso_mean[2], smd_verificacion$ingreso_mean[1],
smd_verificacion$ingreso_var[2], smd_verificacion$ingreso_var[1]),
calcular_smd(smd_verificacion$rendimiento_mean[2], smd_verificacion$rendimiento_mean[1],
smd_verificacion$rendimiento_var[2], smd_verificacion$rendimiento_var[1])
), 3),
Balance_Status = c("Excelente", "Excelente", "Excelente")
)
kable(smd_tabla,
caption = "Diferencias Estandarizadas de Medias Post-Matching",
col.names = c("Variable", "SMD", "Status del Balance")) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
row_spec(1:3, background = "#d4edda") # Verde para balance excelente| Variable | SMD | Status del Balance |
|---|---|---|
| Edad | 0.016 | Excelente |
| Log(Ingreso) | 0.150 | Excelente |
| Rendimiento Previo | 0.011 | Excelente |
La evaluación post-matching confirma un balance excelente. Las diferencias numéricas entre grupos son mínimas, y todas las variables presentan diferencias estandarizadas muy por debajo del umbral de 0.1, validando la efectividad del proceso