## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5) options(digits = 4) # vayr and estimatr are in Suggests, so the chunks that use them are guarded on # both being installed and on vayr being new enough to have # impute_extreme_values(). Everything else runs either way. has_suggests <- requireNamespace("vayr", quietly = TRUE) && requireNamespace("estimatr", quietly = TRUE) && utils::packageVersion("vayr") >= "1.1.0" ## ----setup, message = FALSE, eval = has_suggests------------------------------ library(attrition) library(ggplot2) library(dplyr) library(purrr) library(vayr) library(estimatr) ## ----include = FALSE, eval = !has_suggests------------------------------------ # library(attrition) # library(ggplot2) # library(dplyr) # library(purrr) ## ----------------------------------------------------------------------------- levendusky_replication |> count(Z_condition, R1) ## ----eval = has_suggests------------------------------------------------------ naive_fit <- difference_in_means(formula = Y_polarization_w2 ~ Z, data = filter(levendusky_replication, R1 == 1)) tidy(naive_fit) ## ----------------------------------------------------------------------------- estimator_ev(Y = Y_polarization_w2, Z = Z, R = R1, minY = 0, maxY = 6, data = levendusky_replication) ## ----eval = has_suggests------------------------------------------------------ bounded <- levendusky_replication |> mutate(Y_polarization_w2 = if_else(R1 == 1, Y_polarization_w2, NA_real_)) |> impute_extreme_values(outcome = "Y_polarization_w2", assignment = "Z", range = c(0, 6)) bound_means <- bounded |> group_by(Z_condition, scenario) |> summarise(Y_polarization_w2 = mean(Y_polarization_w2), .groups = "drop") ## ----eval = has_suggests, fig.alt = "Two panels of the same 1,980 subjects. In the lower bound panel the imputed outcomes sit at 0 for the polarized group and at 6 for the moderate group; the upper bound panel reverses them."---- label_df <- bounded |> distinct(scenario, imputed) |> filter(scenario == "Lower bound") |> mutate( Z_condition = "Moderate", Y_polarization_w2 = if_else(imputed == "Outcome imputed", 6.4, 4.9), label = if_else(imputed == "Outcome imputed", "imputed", "reported") ) ggplot(data = bounded, mapping = aes(x = Z_condition, y = Y_polarization_w2)) + geom_point(mapping = aes(colour = imputed, shape = imputed), position = position_jitter(width = 0.25, height = 0), alpha = 0.5, stroke = 0) + geom_point(data = bound_means, size = 3) + geom_text(data = label_df, mapping = aes(label = label, colour = imputed), hjust = 0, nudge_x = 0.15, size = 3.2, show.legend = FALSE) + facet_wrap(facets = ~ scenario) + scale_colour_manual(values = c("#205C8A", "#C67800")) + scale_y_continuous(breaks = 0:6) + labs(x = NULL, y = "Perceived polarization (0 to 6)") + theme_minimal() + theme(legend.position = "none") ## ----eval = has_suggests------------------------------------------------------ bounded |> group_by(scenario) |> reframe(tidy(lm_robust(formula = Y_polarization_w2 ~ Z))) |> filter(term == "Z") |> select(scenario, difference_in_means = estimate) ## ----eval = has_suggests------------------------------------------------------ ev <- estimator_ev(Y = Y_polarization_w2, Z = Z, R = R1, minY = 0, maxY = 6, data = levendusky_replication) bind_rows( bounded |> group_by(scenario) |> reframe(tidy(lm_robust(formula = Y_polarization_w2 ~ Z))) |> filter(term == "Z") |> summarise(interval = "Stacking the two panels", lower = min(conf.low), upper = max(conf.high)), tidy(ev) |> filter(term == "bounds") |> transmute(interval = "Imbens-Manski", lower = conf.low, upper = conf.high) ) |> mutate(width = upper - lower) ## ----------------------------------------------------------------------------- levendusky_replication |> count(Z_condition, Attempt, R2) ## ----------------------------------------------------------------------------- estimator_ds(Y = Y_polarization_w2, Z = Z, R1 = R1, Attempt = Attempt, R2 = R2, minY = 0, maxY = 6, data = levendusky_replication) ## ----------------------------------------------------------------------------- estimator_ds(Y = Y_polarization_w2, Z = Z, R1 = R1, Attempt = Attempt, R2 = R2, strata = X_party_id, minY = 0, maxY = 6, data = levendusky_replication) ## ----------------------------------------------------------------------------- estimator_ds_sens(Y = Y_polarization_w2, Z = Z, R1 = R1, Attempt = Attempt, R2 = R2, delta = 0.5, minY = 0, maxY = 6, data = levendusky_replication) ## ----fig.alt = "Identification regions and confidence intervals as a function of the sensitivity parameter delta"---- sens <- sensitivity_ds(Y = Y_polarization_w2, Z = Z, R1 = R1, Attempt = Attempt, R2 = R2, minY = 0, maxY = 6, alpha = 0.10, data = levendusky_replication) sens sens$sensitivity_plot ## ----------------------------------------------------------------------------- sensitivity_ds(Y = Y_polarization_w2, Z = Z, R1 = R1, Attempt = Attempt, R2 = R2, minY = 0, maxY = 6, alpha = 0.05, data = levendusky_replication)$delta_star ## ----------------------------------------------------------------------------- levendusky_replication |> group_by(Z_condition) |> summarise(response_rate = mean(R1)) ## ----------------------------------------------------------------------------- tidy(estimator_trim(Y = Y_polarization_w2, Z = Z, R = R1, data = levendusky_replication)) ## ----------------------------------------------------------------------------- tidy(estimator_trim(Y = Y_polarization_w2, Z = Z, R = R1, monotonicity = "treatment_decreases_response", data = levendusky_replication)) ## ----------------------------------------------------------------------------- cells <- list( "single sample, monotonicity" = estimator_trim(Y = Y_polarization_w2, Z = Z, R = R1, monotonicity = "treatment_decreases_response", se = "none", data = levendusky_replication), "single sample, no monotonicity" = estimator_trim(Y = Y_polarization_w2, Z = Z, R = R1, monotonicity = "none", se = "none", data = levendusky_replication), "double sampling, monotonicity" = estimator_trim(Y = Y_polarization_w2, Z = Z, R1 = R1, Attempt = Attempt, R2 = R2, monotonicity = "treatment_decreases_response", se = "none", data = levendusky_replication), "double sampling, no monotonicity" = estimator_trim(Y = Y_polarization_w2, Z = Z, R1 = R1, Attempt = Attempt, R2 = R2, monotonicity = "none", se = "none", data = levendusky_replication) ) cells |> map(tidy) |> list_rbind(names_to = "cell") |> filter(term == "bounds") |> transmute(cell, estimate_lower, estimate_upper, width = estimate_upper - estimate_lower) ## ----------------------------------------------------------------------------- set.seed(343) tidy(estimator_trim(Y = Y_polarization_w2, Z = Z, R1 = R1, Attempt = Attempt, R2 = R2, se = "bootstrap", sims = 500, data = levendusky_replication)) ## ----------------------------------------------------------------------------- out <- estimator_ds(Y = Y_polarization_w2, Z = Z, R1 = R1, Attempt = Attempt, R2 = R2, minY = 0, maxY = 6, data = levendusky_replication) out ## ----------------------------------------------------------------------------- tidy(out) ## ----------------------------------------------------------------------------- estimator_ds(Y = Y_polarization_w2 ~ Z, R1 = "R1", Attempt = "Attempt", R2 = "R2", minY = 0, maxY = 6, data = levendusky_replication) ## ----------------------------------------------------------------------------- set.seed(343) fits <- list( "Extreme value" = estimator_ev(Y = Y_polarization_w2, Z = Z, R = R1, minY = 0, maxY = 6, data = levendusky_replication), "Extreme value, double sampled" = estimator_ds(Y = Y_polarization_w2, Z = Z, R1 = R1, Attempt = Attempt, R2 = R2, minY = 0, maxY = 6, data = levendusky_replication), "Extreme value, DS, poststratified" = estimator_ds(Y = Y_polarization_w2, Z = Z, R1 = R1, Attempt = Attempt, R2 = R2, strata = X_party_id, minY = 0, maxY = 6, data = levendusky_replication), "Sensitivity at delta = 0.5" = estimator_ds_sens(Y = Y_polarization_w2, Z = Z, R1 = R1, Attempt = Attempt, R2 = R2, delta = 0.5, minY = 0, maxY = 6, data = levendusky_replication), "Trimming, monotonicity" = estimator_trim(Y = Y_polarization_w2, Z = Z, R = R1, monotonicity = "treatment_decreases_response", se = "bootstrap", sims = 500, data = levendusky_replication), "Trimming, monotonicity, DS" = estimator_trim(Y = Y_polarization_w2, Z = Z, R1 = R1, Attempt = Attempt, R2 = R2, monotonicity = "treatment_decreases_response", se = "bootstrap", sims = 500, data = levendusky_replication), "Trimming, no monotonicity" = estimator_trim(Y = Y_polarization_w2, Z = Z, R = R1, monotonicity = "none", se = "bootstrap", sims = 500, data = levendusky_replication), "Trimming, no monotonicity, DS" = estimator_trim(Y = Y_polarization_w2, Z = Z, R1 = R1, Attempt = Attempt, R2 = R2, monotonicity = "none", se = "bootstrap", sims = 500, data = levendusky_replication) ) estimands <- tribble( ~estimator, ~estimand, "Naive difference in means", "Wave 2\nrespondents", "Extreme value", "All 1,980\nsubjects", "Extreme value, double sampled", "All 1,980\nsubjects", "Extreme value, DS, poststratified", "All 1,980\nsubjects", "Sensitivity at delta = 0.5", "All 1,980\nsubjects", "Trimming, monotonicity", "Always\nreporters", "Trimming, monotonicity, DS", "Always\nreporters", "Trimming, no monotonicity", "Always\nreporters", "Trimming, no monotonicity, DS", "Always\nreporters" ) # The naive estimator claims a point rather than a region, which the plot shows as # an identification region of zero width naive <- lm(Y_polarization_w2 ~ Z, data = filter(levendusky_replication, R1 == 1)) gg_df <- fits |> map(\(fit) filter(tidy(fit), term == "bounds")) |> list_rbind(names_to = "estimator") |> select(estimator, estimate_lower, estimate_upper, conf.low, conf.high) |> add_row(estimator = "Naive difference in means", estimate_lower = coef(naive)[["Z"]], estimate_upper = coef(naive)[["Z"]], conf.low = confint(naive)["Z", 1], conf.high = confint(naive)["Z", 2], .before = 1) |> left_join(estimands, by = "estimator") |> mutate( estimator = factor(estimator, levels = rev(estimands$estimator)), estimand = factor(estimand, levels = unique(estimands$estimand)) ) ## ----fig.height = 5, fig.alt = "Identification regions and confidence intervals for nine estimators, grouped by the population each one is about."---- ggplot(data = gg_df, mapping = aes(y = estimator)) + geom_vline(xintercept = 0, linetype = "dashed", colour = "grey65") + geom_linerange(mapping = aes(xmin = conf.low, xmax = conf.high), linewidth = 0.5, colour = "#9AA5AE") + geom_linerange(mapping = aes(xmin = estimate_lower, xmax = estimate_upper), linewidth = 2.4, colour = "#205C8A") + geom_point(data = filter(gg_df, estimate_lower == estimate_upper), mapping = aes(x = estimate_lower), size = 2.2, colour = "#205C8A") + facet_grid(rows = vars(estimand), scales = "free_y", space = "free_y", switch = "y") + labs(x = "Effect on perceived polarization", y = NULL, subtitle = "Thick: identification region\nThin: 95 percent confidence interval") + theme_minimal(base_size = 10) + theme(strip.placement = "outside", strip.text.y.left = element_text(angle = 0, hjust = 0), panel.grid.major.y = element_blank(), plot.subtitle = element_text(size = 8.5, colour = "grey30"))