## ----setup-knitr, include = FALSE--------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 6.5, fig.height = 4, fig.align = "center", out.width = "92%" ) ## ----setup, message = FALSE--------------------------------------------------- library(tvbounds) library(ggplot2) ## ----rct-simulate------------------------------------------------------------- set.seed(20260820) n <- 500 village <- rep(1:50, each = 10) d <- as.integer(rbinom(50, 1, 0.5)[village]) # village-level assignment x <- rbinom(n, 1, 0.4) # baseline covariate ability <- rnorm(n) s <- as.integer(runif(n) < plogis(0.2 + 1.2 * d + 0.5 * ability)) y <- ifelse(s == 1, 1 + 0.35 * d + 0.5 * x + ability + 0.5 * rnorm(n), NA) rct <- data.frame(y = y, d = d, s = s, x = x, village = village) tapply(rct$s, rct$d, mean) # response rates by arm ## ----rct-tv------------------------------------------------------------------- fit_tv <- tvbounds_attrition(rct, outcome = "y", treatment = "d", response = "s", delta = seq(0, 1, by = 0.05), B = 200, seed = 1) fit_tv ## ----rct-pstar---------------------------------------------------------------- fit_tv$details$p_star ## ----rct-bounds-head---------------------------------------------------------- head(fit_tv$bounds, 4) ## ----rct-plot, fig.cap = "Total variation bounds for the simulated experiment."---- plot(fit_tv) ## ----rct-endpoint-0----------------------------------------------------------- naive <- mean(rct$y[rct$d == 1 & rct$s == 1]) - mean(rct$y[rct$d == 0 & rct$s == 1]) c(point = fit_tv$point, naive = naive) ## ----rct-endpoint-1----------------------------------------------------------- subset(fit_tv$bounds, delta == 1, select = c(delta, lower, upper)) unlist(fit_tv$details$lee) ## ----rct-contamination-------------------------------------------------------- fit_ct <- tvbounds_attrition(rct, outcome = "y", treatment = "d", response = "s", delta = seq(0, 1, by = 0.05), neighborhood = "contamination", bootstrap = FALSE) all(fit_ct$bounds$lower >= fit_tv$bounds$lower - 1e-12) # nesting all(fit_ct$bounds$upper <= fit_tv$bounds$upper + 1e-12) ## ----rct-compare-plot, fig.cap = "Total variation versus contamination bounds."---- comp <- rbind( cbind(fit_tv$bounds[c("delta", "lower", "upper")], neighborhood = "total variation"), cbind(fit_ct$bounds[c("delta", "lower", "upper")], neighborhood = "contamination")) ggplot(comp, aes(x = delta)) + geom_line(aes(y = lower, linetype = neighborhood), color = "#1F4E79", linewidth = 0.9) + geom_line(aes(y = upper, linetype = neighborhood), color = "#1F4E79", linewidth = 0.9) + geom_hline(yintercept = 0, linetype = "dotted", color = "gray50") + labs(x = expression(delta), y = "treatment effect") + theme_bw(base_size = 12) + theme(panel.grid = element_blank(), legend.position = "bottom") ## ----rct-cluster-------------------------------------------------------------- fit_cl <- tvbounds_attrition(rct, outcome = "y", treatment = "d", response = "s", delta = seq(0, 1, by = 0.05), B = 200, cluster = "village", seed = 1) fit_cl$details$n_clusters data.frame(delta = fit_tv$bounds$delta, se_iid = fit_tv$bounds$lower_se, se_cluster = fit_cl$bounds$lower_se)[1:4, ] ## ----rct-covariates----------------------------------------------------------- fit_x <- tvbounds_attrition(rct, outcome = "y", treatment = "d", response = "s", covariates = "x", delta = seq(0, 1, by = 0.05), B = 200, seed = 1) fit_x$details$pooled$strata ## ----rct-pooled-vs-pw--------------------------------------------------------- pw <- fit_x$details$pooled$pw all(fit_x$bounds$upper >= pw$upper - 1e-12) all(fit_x$bounds$lower <= pw$lower + 1e-12) ## ----rct-lee-cov-------------------------------------------------------------- rbind(with_covariate = unlist(fit_x$details$lee), without_covariate = unlist(fit_x$details$lee_nocov)) ## ----rct-summary-------------------------------------------------------------- summary(fit_tv) ## ----rct-summary-delta-------------------------------------------------------- s01 <- tvbounds_summary(fit_tv, delta = 0.1) s01$measures[c("delta_eval", "eta", "varsigma", "n_star")] ## ----riv-simulate------------------------------------------------------------- set.seed(1901) n <- 150; K <- 10; S <- 80 W <- matrix(rexp(n * K)^2, n, K) W <- W / rowSums(W) # exposure shares g0 <- rnorm(K, mean = 0.3) # realized sector shocks G <- matrix(rnorm(K * S), K, S) # S draws from the postulated process e_x <- rnorm(n); e_s <- rnorm(K); e_y <- rnorm(n) z <- as.vector(W %*% g0) # realized formula instrument Fmat <- W %*% G # n x S counterfactual draws x <- z + e_x # endogenous regressor (first stage) y <- 0.5 * x + 0.4 * as.vector(W %*% e_s) + 0.5 * e_y ## ----riv-fit------------------------------------------------------------------ riv_tv <- tvbounds_riv(y, x, z, Fmat, delta = seq(0, 1, by = 0.01)) riv_tv riv_ct <- tvbounds_riv(y, x, z, Fmat, delta = seq(0, 1, by = 0.01), neighborhood = "contamination") ## ----riv-nesting-------------------------------------------------------------- all(riv_ct$bounds$lower >= riv_tv$bounds$lower - 1e-10, na.rm = TRUE) ## ----riv-plot, fig.cap = "Total variation bounds for the recentered IV estimate."---- plot(riv_tv) ## ----riv-summary-------------------------------------------------------------- riv_tv$details$delta_breakdown summary(riv_tv) ## ----riv-fs------------------------------------------------------------------- c(delta_fs = riv_tv$details$delta_fs, censored = riv_tv$details$delta_fs_censored) ## ----riv-weak----------------------------------------------------------------- x_weak <- 0.2 * z + e_x # same noise, weaker first stage y_weak <- 0.5 * x_weak + 0.4 * as.vector(W %*% e_s) + 0.5 * e_y riv_weak <- tvbounds_riv(y_weak, x_weak, z, Fmat, delta = seq(0, 1, by = 0.01)) c(delta_fs = riv_weak$details$delta_fs, censored = riv_weak$details$delta_fs_censored) tail(riv_weak$bounds, 3) # vacuous budgets are NA ## ----riv-weak-plot, fig.cap = "A weak first stage: the bounds diverge and become vacuous at the first-stage breakdown budget."---- plot(riv_weak, breakdown = FALSE) ## ----cf-toy, eval = FALSE----------------------------------------------------- # toy <- system.file("julia", "examples", "toy.jl", package = "tvbounds") # # fit_cf <- tvbounds_counterfactual( # moments = c(toy, "tvb_toy_moments!"), # Julia file + function name # d = 1, # number of moment conditions # theta_lb = 0.4, theta_ub = 0.6, # box for the structural parameter # delta = c(0.05, 0.1, 0.25, 0.5, 1), # budgets (strictly positive) # divergence = "TVmix", # side = "both", # lower and upper bound problems # M = 5000, u_dim = 1, # scrambled-Halton draws # theta_init = 0.5, # baseline parameter estimate # seed = 1234, # control = tvbounds_control(maxsolves = 5)) # # fit_cf$bounds # one row per budget # fit_cf$point # plug-in counterfactual at theta_init # fit_cf$details$solver # per-budget KNITRO status codes and timings # plot(fit_cf) ## ----cf-rmoments, eval = FALSE------------------------------------------------ # r_moments <- function(theta, U, gamma) { # u <- qnorm(U[, 1]) # map uniforms into baseline draws # list(K = u + theta[1], # counterfactual values # G = cbind(u - theta[1])) # moment conditions (M x d) # } # # fit_r <- tvbounds_counterfactual( # moments = r_moments, d = 1, # theta_lb = -0.5, theta_ub = 0.5, # delta = c(0.1, 0.5), # divergence = "TV", # M = 2000, u_dim = 1, theta_init = 0, # seed = 1234) ## ----cf-control, eval = FALSE------------------------------------------------- # ctrl <- tvbounds_control( # maxsolves = 3, # multi-start restarts per budget/side # startptrange = 0.05, # spread of the restart perturbations # outer_opt = system.file("julia", "opt", "outer_fast.opt", # package = "tvbounds"), # knitro_options = list(maxit = 300, outlev = 0)) # # fit_fast <- tvbounds_counterfactual( # moments = c(toy, "tvb_toy_moments!"), d = 1, # theta_lb = 0.4, theta_ub = 0.6, # delta = seq(0.1, 1, by = 0.1), divergence = "TVmix", # M = 5000, u_dim = 1, theta_init = 0.5, seed = 1234, # control = ctrl)