library(tidyverse) # ggplot(), %>%, mutate(), and friends
library(ggdag) # Make DAGs
library(scales) # Format numbers with functions like comma(), percent(), and dollar()
library(broom) # Convert models to data frames
library(patchwork) # Combine ggplots into single composite plots
set.seed(1234) # Make all random draws reproducibleExperimentos Aleatorios Controlados (EAC)
Para este ejemplo vamos a usar el archivo (village_randomized.csv)
village_randomized <- read_csv("data/village_randomized.csv")Detalles del programa
En esta situación hipotética, una ONG planea lanzar un programa de capacitación diseñado para aumentar los ingresos. Basándose en su experiencia en programas piloto en otros países, han descubierto que los hombres mayores y con mayor poder adquisitivo tienden a autoseleccionarse para participar en el programa de capacitación. El consultor de evaluación de la ONG (¡tú!) diseñó este modelo causal que explica el efecto del programa en los ingresos de los participantes, considerando la confusión causada por la edad, el sexo y los ingresos previos:
income_dag <- dagify(post_income ~ program + age + sex + pre_income,
program ~ age + sex + pre_income,
exposure = "program",
outcome = "post_income",
labels = c(post_income = "Post income",
program = "Program",
age = "Age",
sex = "Sex",
pre_income = "Pre income"),
coords = list(x = c(program = 1, post_income = 5, age = 2,
sex = 4, pre_income = 3),
y = c(program = 2, post_income = 2, age = 1,
sex = 1, pre_income = 3)))
ggdag_status(income_dag, use_labels = "label", text = FALSE, seed = 1234) +
guides(color = "none") +
theme_dag()
La ONG acaba de recibir financiación para realizar un ensayo controlado aleatorio (ECA) en un pueblo, y usted está entusiasmado porque finalmente puede manipular el acceso al programa: podemos calcular \(E(\text{Post-income} | do(\text{Program})\) siguiendo las reglas de los diagramas causales, podemos eliminar todas las flechas que apuntan al nodo del programa:
income_dag_rct <- dagify(post_income ~ program + age + sex + pre_income,
exposure = "program",
outcome = "post_income",
labels = c(post_income = "Post income",
program = "Program",
age = "Age",
sex = "Sex",
pre_income = "Pre income"),
coords = list(x = c(program = 1, post_income = 5, age = 2,
sex = 4, pre_income = 3),
y = c(program = 2, post_income = 2, age = 1,
sex = 1, pre_income = 3)))
ggdag_status(income_dag_rct, use_labels = "label", text = FALSE, seed = 1234) +
guides(color = "none") +
theme_dag()
1. Checamos balance
Realizaste el estudio con 1.000 participantes a lo largo de 6 meses y acabas de recibir los datos.
Antes de calcular el efecto del programa, primero verificas qué tan bien equilibrada fue la asignación y descubres que la asignación al programa estuvo prácticamente dividida 50/50, lo cual es bueno:
village_randomized %>%
count(program) %>%
mutate(prop = n / sum(n))
## # A tibble: 2 × 3
## program n prop
## <chr> <int> <dbl>
## 1 No program 503 0.503
## 2 Program 497 0.497Luego, verifica qué tan equilibrados estaban los grupos de tratamiento y control en cuanto a las características previas al tratamiento de los participantes:
village_randomized %>%
group_by(program) %>%
summarize(prop_male = mean(sex_num),
avg_age = mean(age),
avg_pre_income = mean(pre_income))
## # A tibble: 2 × 4
## program prop_male avg_age avg_pre_income
## <chr> <dbl> <dbl> <dbl>
## 1 No program 0.584 34.9 803.
## 2 Program 0.604 34.9 801.Estas variables parecen estar bastante bien equilibradas. Para comprobar que no existen diferencias estadísticamente significativas entre los grupos, se crean gráficos (también se podrían realizar pruebas t, pero los gráficos son más fáciles de leer posteriormente para el público no especializado en estadística).
Había más hombres tanto en el grupo de tratamiento como en el de control, pero la proporción es la misma en ambos y no hay una diferencia sustancial en la proporción de sexos:
# Here we save each plot as an object so that we can combine the two plots with
# + (which comes from the patchwork package). If you want to see what an
# individual plot looks like, you can run `plot_diff_sex`, or whatever the plot
# object is named.
#
# stat_summary() here is a little different from the geom_*() layers you've seen
# in the past. stat_summary() takes a function (here mean_se()) and runs it on
# each of the program groups to get the average and standard error. It then
# plots those with geom_pointrange. The fun.args part of this lets us pass an
# argument to mean_se() so that we can multiply the standard error by 1.96,
# giving us the 95% confidence interval.
plot_diff_sex <- ggplot(village_randomized, aes(x = program, y = sex_num, color = program)) +
stat_summary(geom = "pointrange", fun.data = "mean_se", fun.args = list(mult = 1.96)) +
guides(color = "none") +
labs(x = NULL, y = "Proportion male")
# plot_diff_sex # Uncomment this if you want to see this plot by itself
plot_prop_sex <- ggplot(village_randomized, aes(x = program, fill = sex)) +
# Using position = "fill" makes the bars range from 0-1 and show the proportion
geom_bar(position = "fill") +
labs(x = NULL, y = "Proportion", fill = NULL) +
scale_fill_manual(values = c("darkblue", "darkred"))
# Show the plots side-by-side
plot_diff_sex + plot_prop_sex
Los ingresos previos al programa también se distribuyen de la misma manera (y no hay diferencias sustanciales en los promedios) entre los grupos de tratamiento y de control:
plot_diff_age <- ggplot(village_randomized, aes(x = program, y = age, color = program)) +
stat_summary(geom = "pointrange", fun.data = "mean_se", fun.args = list(mult = 1.96)) +
guides(color = "none") +
labs(x = NULL, y = "Age")
plot_hist_age <- ggplot(village_randomized, aes(x = age, fill = program)) +
geom_histogram(binwidth = 1, color = "white") +
guides(fill = "none") +
labs(x = "Age", y = "Count") +
facet_wrap(vars(program), ncol = 1)
plot_diff_age + plot_hist_age
¡Todas la covariates previas al tratamiento se ven bien y equilibradas! Ahora podemos estimar el efecto causal del programa.
2. Estimamos la diferencia
Nos interesa estiumar el efecto causal del programa, o
\[ E[\text{Post income}\ |\ do(\text{Program})] \]
Ese efecto causal lo podemos estimar a través del efecto del tratamiento promedio (ATE)
$$ = E( | = 1) - E( | = 0)
Esto es simplemente el resultado promedio de las personas que participan en el programa menos el resultado promedio de las personas que no participan. Se calculan los promedios de grupo:
village_randomized %>%
group_by(program) %>%
summarize(avg_post = mean(post_income))
## # A tibble: 2 × 2
## program avg_post
## <chr> <dbl>
## 1 No program 1180.
## 2 Program 1279.Eso es 1279 − 1180, o 99, lo que significa que el programa provocó un aumento de $99 en los ingresos, en promedio.
Encontrar esa diferencia requirió cálculos manuales, así que, como atajo, se ejecuta un modelo de regresión con los ingresos posteriores al programa como variable de resultado y la variable indicadora del programa como variable explicativa. El coeficiente para programes el efecto causal (e incluye información sobre los errores estándar y la significancia). Se obtiene el mismo resultado:
model_rct <- lm(post_income ~ program, data = village_randomized)
tidy(model_rct)
## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 1180. 4.27 276. 0
## 2 programProgram 99.2 6.06 16.4 1.26e-53Con base en su RCT, concluyes que el programa genera un aumento promedio de $99,25 en ingresos.
ggplot(village_randomized, aes(x = program, y = post_income, color = program)) +
stat_summary(geom = "pointrange", fun.data = "mean_se", fun.args = list(mult = 1.96)) +
guides(color = "none") +
labs(x = NULL, y = "Post income")