POLI 891: External Validity Replication
Seth L & Marta S
2025-10-22
What does Evalid do?
Evalid’s two primary uses are the T-Pate function and the PCT(Partial Conjunction Test). In this replication we’ll be demonstrating the T-pate function by replicating the study conducted by Broockman & Kalla (2016) which investigated how attitudes to transphobia may be reduced through door-to-door canvassing. Conducting in Miami Dade county Flordia in 2015 among those who completed a pre-experiment survey. Canvasser’s were randomly assigned to either take the perspective of trans individuals, or a placebo conversation. The outcome was formed by four waves of follow-up surveys which indicated that perspective taking conversations improved the outlook toward trans individuals.
#Evalid is the proprietary package developed by Egami & Hartman
library(evalid)
library(grf)
library(foreign)
library(lmtest)
library(sandwich)
library(survey)
library(MASS)
library(parallel)
library(estimatr)
library(bartCause)
library(tidyverse)
library(knitr)
library(kableExtra)
library(estimatr)set.seed(29384)
#recoded results
load("~/R Projects/Poli 819/generated/recoded_broockman_kalla_data.RData")
#previous results
load ("~/R Projects/Poli 819/generated/results_bk_external.RData")
#number of bootstrapped samples.
nboot = 1000## Covariates used for estimation
## Projection models: uses raw age
covariates_proj <- c('vf_age', 'vf_female', 'vf_black', 'vf_white',
names(sample)[str_detect(names(sample), "ideology_t0_factor|religious_t0_factor|pid_t0_factor")][-c(1:3)])
## Weighting models: uses age buckets
covariates <- c('vf_female', 'vf_black', 'vf_white',
names(sample)[str_detect(names(sample), "ideology_t0_factor|religious_t0_factor|pid_t0_factor|vf_age_bucket")][-c(1:4)])
# Pre-specified regression controls used in original analysis
control_vars <- c('miami_trans_law_t0', 'miami_trans_law2_t0', 'therm_trans_t0',
'gender_norms_sexchange_t0', 'gender_norms_moral_t0', 'gender_norms_abnormal_t0',
'ssm_t0', 'therm_obama_t0', 'therm_gay_t0','vf_democrat', 'ideology_t0',
'religious_t0', 'exposure_gay_t0', 'exposure_trans_t0', 'pid_t0', 'sdo_scale',
'gender_norm_daugher_t0', 'gender_norm_looks_t0',
'gender_norm_rights_t0', 'therm_afams_t0', 'vf_female', 'vf_hispanic',
'vf_black', 'vf_age', 'survey_language_es', 'cluster_level_t0_scale_mean')Trans. Tolerance Index Analysis
This section re-analyzes the initial findings from the paper, using the different PATE estimation strategies. The estimators are:
- SATE (Regression) is the bivariate regression with pre-specified controls and HC2 errors, accounting for household-level clusters.
- PATE (IPW) is the weighting based estimator estimated using calibration weights.
- PATE (wLS) is the IPW weighted least squares, controlling for pre-specified controls, using calibration weights.
- PATE (OLS-Outcome) is OLS outcome based estimator (estimated with overlap covariates), estimated separately for Y1 and Y0, projected on to the population. The doubly robust AIPW estimator, using calibration weights, is also included.
- PATE (BART-Outcome) is BART outcome based estimator (estimated with overlap covariates) projected on to the population. The doubly robust AIPW estimator, using calibration weights, is also included.
## Only re-runs analysis if no results are loaded.
## run analysis for each time period
if(!exists("plot_data_tti")) {
plot_data_tti <- NULL
outcome_vars <- c("trans.tolerance.dv.t1" , "trans.tolerance.dv.t2",
"trans.tolerance.dv.t3", "trans.tolerance.dv.t4")
for(outcome in outcome_vars) {
# limit to rows without missing outcome
exp_data <- sample %>% filter(!is.na(sample[, outcome]))
# create formulas for analysis with outcome
formula_outcome <- as.formula(paste0(outcome, " ~ ", paste0(covariates, collapse = " + ")))
formula_outcome_proj <- as.formula(paste0(outcome, " ~ ", paste0(covariates_proj, collapse = " + ")))
formula_outcome_wls <- as.formula(paste0(outcome, " ~ ", paste0(control_vars, collapse = " + ")))
formula_weights <- as.formula(paste0("S ~ ", paste0(covariates, collapse = " + ")))
## Sate Estimate
sate_est <- data.frame(estimator = c("SATE"), type = "SATE",
tpate(formula_outcome = formula_outcome_wls,
formula_weights = formula_weights,
exp_data = exp_data, pop_data = exp_data,
treatment = "treat_ind",
id_cluster = exp_data$block_ind,
est_type = "wls", weights_type = "calibration",
sims = nboot,
compute_sate = FALSE)$tpate[c("est", "se", "ci_lower", "ci_upper")] %>% data.frame())
## PATE Estimates
## Weighting: IPW
ipw_cal <- data.frame(estimator = "ipw-cal", type = "weighting",
tpate(formula_outcome = formula_outcome,
formula_weights = formula_weights,
exp_data = exp_data, pop_data = pop,
treatment = "treat_ind",
id_cluster = exp_data$block_ind,
est_type = "ipw", weights_type = "calibration",
sims = nboot,
weights_max = 10,
compute_sate = FALSE)$tpate[c("est", "se", "ci_lower", "ci_upper")]
%>% data.frame())
## Weighting: Weighted Least Squares
wls_cal <- data.frame(estimator = "wls-cal", type = "weighting",
tpate(formula_outcome = formula_outcome_wls,
formula_weights = formula_weights,
exp_data = exp_data, pop_data = pop,
treatment = "treat_ind",
id_cluster = exp_data$block_ind,
est_type = "wls", weights_type = "calibration",
sims = nboot,
weights_max = 10,
compute_sate = FALSE)$tpate[c("est", "se", "ci_lower", "ci_upper")]
%>% data.frame())
## Outcome: Regression
out_ols <- data.frame(estimator = "OLS-proj", type = "outcome",
tpate(formula_outcome = formula_outcome_proj,
formula_weights = formula_weights,
exp_data = exp_data, pop_data = pop,
treatment = "treat_ind",
id_cluster = exp_data$block_ind,
est_type = "outcome-ols",
sims = nboot,
weights_max = 10,
compute_sate = FALSE)$tpate[c("est", "se", "ci_lower", "ci_upper")]
%>% data.frame())
## Outcome: BART
out_bart <- data.frame(estimator = "BART-proj", type = "outcome",
tpate(formula_outcome = formula_outcome_proj,
formula_weights = formula_weights,
exp_data = exp_data, pop_data = pop,
treatment = "treat_ind",
id_cluster = exp_data$block_ind,
est_type = "outcome-bart",
sims = nboot,
weights_max = 10,
compute_sate = FALSE)$tpate[c("est", "se", "ci_lower", "ci_upper")]
%>% data.frame())
## Doubly Robust Estimators
## Doubly Robust: Augmented-Regression
dr_cal_ols <- data.frame(estimator = "DR-OLS-cal", type = "doubly robust",
tpate(formula_outcome = formula_outcome,
formula_weights = formula_weights,
exp_data = exp_data, pop_data = pop,
treatment = "treat_ind",
id_cluster = exp_data$block_ind,
weights_type = "calibration",
est_type = "dr-ols",
sims = nboot,
weights_max = 10,
compute_sate = FALSE)$tpate[c("est", "se", "ci_lower", "ci_upper")]
%>% data.frame())
## Doubly Robust: Augmented-BART
dr_cal_bart <- data.frame(estimator = "DR-BART-cal", type = "doubly robust",
tpate(formula_outcome = formula_outcome_proj,
formula_weights = formula_weights,
exp_data = exp_data, pop_data = pop,
treatment = "treat_ind",
id_cluster = exp_data$block_ind,
weights_type = "calibration",
est_type = "dr-bart",
sims = nboot,
weights_max = 10,
compute_sate = FALSE)$tpate[c("est", "se", "ci_lower", "ci_upper")]
%>% data.frame())
results = data.frame(time = which(outcome_vars == outcome),
bind_rows(
sate_est,
ipw_cal,
wls_cal,
out_ols,
out_bart,
dr_cal_ols,
dr_cal_bart
)
)
plot_data_tti <- bind_rows(plot_data_tti, results)
}
# Format name columns
plot_data_tti$time_name = factor(plot_data_tti$time, labels = c("+3 Days", "+3 Weeks", "+6 Weeks", "+3 Months"))
plot_data_tti$estimator_type = factor(case_when(str_detect(plot_data_tti$type, "SATE") ~ "SATE",
str_detect(plot_data_tti$type, "outcome") ~ "T-PATE:\nOutcome-based Estimator",
str_detect(plot_data_tti$type, "doubly robust") ~ "T-PATE:\nDoubly Robust Estimator",
TRUE ~ "T-PATE:\nWeighting-based Estimator"),
levels = c("SATE", "T-PATE:\nWeighting-based Estimator",
"T-PATE:\nOutcome-based Estimator",
"T-PATE:\nDoubly Robust Estimator"))
plot_data_tti$estimator_name = factor(case_when(plot_data_tti$estimator == "SATE" ~ "OLS",
plot_data_tti$estimator == "ipw-cal" ~ "IPW",
plot_data_tti$estimator == "wls-cal" ~ "wLS",
plot_data_tti$estimator == "OLS-proj" ~ "OLS",
plot_data_tti$estimator == "BART-proj" ~ "BART",
plot_data_tti$estimator == "DR-OLS-cal" ~ "AIPW\nwith OLS",
plot_data_tti$estimator == "DR-BART-cal" ~ "AIPW\nwith BART"),
levels = c("OLS", "IPW", "wLS", "BART", "AIPW\nwith OLS", "AIPW\nwith BART"))}
## Save results
#save(plot_data_tti, file = paste0("./generated/main_results/results_bk_external.RData")) dw = 0.5
plot_data_tti %>%
ggplot(aes(x = estimator_name, y = est, color = estimator_type)) +
scale_color_manual(values=c("black", "red", "blue", "magenta")) +
geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), position = position_dodge(width = dw)) +
geom_point(aes(y = est), position = position_dodge(width = dw)) +
theme_bw() +
xlab("") +
ylab("Estimated Causal Effects on Transgender Tolerance Scale") +
theme(axis.title.y = element_text(size = 13, margin = margin(t = 0, r = 10, b = 0, l = 0)))+
geom_hline(yintercept = 0, linetype = 2) +
theme(legend.position="none",
axis.text.x = element_text(angle = 0, size = 11, margin = margin(t = 10, r = 0, b = 0, l = 0),
colour = "black"),
axis.text = element_text(size = 11),
axis.title = element_text(size = 11),
legend.text = element_text(size = 11),
panel.grid.major = element_blank(),
strip.text.x = element_text(size = 11),
strip.text.y = element_text(size = 13)) +
facet_grid(time_name ~ estimator_type, scales = "free_x")
#ggsave("./generated/Figures/Figure_7.pdf", height = 8, width = 11)