Sensitivity Analysis with tvbounds

2.1 A small simulated experiment

We simulate a randomized experiment with 500 individuals in 50 villages and treatment assigned at the village level. Follow-up response is markedly higher under treatment and depends on an unobserved ability that also drives the outcome, so that attrition is selective and the comparison of respondents across arms confounds the effect of the program with composition effects.

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
#>         0         1 
#> 0.5333333 0.8192308

By construction the treatment effect equals 0.35, whereas only roughly 53% of control units respond against 82% of treated units, and high-ability units respond more often in both arms.

2.2 Total variation bounds

tvbounds_attrition() requires the column names of the outcome \(Y\), of the binary treatment \(D\), and of the response indicator \(S\), which equals 1 when the outcome is observed and may be paired with an NA outcome otherwise. We compute the bounds on a grid of budgets, with a small number of bootstrap replications for the purposes of this vignette; the seed argument makes the bootstrap reproducible without altering the caller’s random-number stream.

fit_tv <- tvbounds_attrition(rct,
  outcome = "y", treatment = "d", response = "s",
  delta = seq(0, 1, by = 0.05), B = 200, seed = 1)
fit_tv
#> <tvbounds> attrition bounds under a total-variation neighborhood
#>   estimand: treatment effect; baseline point estimate (delta = 0): 0.465; n = 500
#>   budget grid: 21 values of delta in [0, 1]
#>   inference: 95% bootstrap confidence bands (B = 200)
#>   Use summary() for breakdown and price-of-robustness measures; plot() to display.

The estimand is the average treatment effect for the always-observed subpopulation, the units that would respond under either arm: \[ \tau_0 := \mathbb{E}_{P_0}[Y(1) - Y(0) \mid S(0) = 1, S(1) = 1], \] where \(Y(d)\) and \(S(d)\) are the potential outcome and the potential response indicator under treatment \(d\), and \(P_0\) is the distribution of the data. Under random assignment and the monotonicity assumption of Lee (2009) there are no defiers, so that the outcome distribution of the treated respondents, \(P_{\mathsf{T}}\), is a mixture of the outcome distribution of the compliers, \(P_{\mathsf{C}}\), the units that respond only because they were treated, and that of the always-observed units, \(P_{\mathsf{AO}}\): \[ P_{\mathsf{T}} = \pi P_{\mathsf{C}} + (1 - \pi) P_{\mathsf{AO}}, \qquad \pi = 1 - \frac{r_0}{r_1}. \] The mixing weight \(\pi\), the complier share among treated respondents, is identified by comparing the response rates of the two arms, \(r_1\) and \(r_0\); the package reports the complier share \(\pi\) as p_star:

fit_tv$details$p_star
#> [1] 0.3489828

Neither \(P_{\mathsf{C}}\) nor \(P_{\mathsf{AO}}\) is identified, and missingness completely at random requires that the two coincide. A natural way to relax it is thus to let them lie at most \(\delta\) apart in total variation, \(\mathsf{TV}(P_{\mathsf{C}} \,\|\, P_{\mathsf{AO}}) \le \delta\), so that the two outcome distributions may disagree on at most a \(\delta\)-fraction of their mass. Rescaling that restriction onto the identified law \(P_{\mathsf{T}}\) pins down one of the two distributions and lets the candidate complier distribution \(Q\), of which \(P_{\mathsf{C}}\) is one, range over the robustness set \[ \mathcal{Q}_{\mathsf{C}}(\delta) := \{Q \in \Delta(\mathcal{Y}) : \mathsf{TV}(Q \,\|\, P_{\mathsf{T}}) \le (1 - \pi)\delta, \;\; \pi Q \le P_{\mathsf{T}}\}, \] where \(\Delta(\mathcal{Y})\) collects the distributions on the outcome space and the second restriction encodes the mixture structure of the treated arm. The resulting sensitivity bounds \(\underline{\tau}(\delta)\) and \(\overline{\tau}(\delta)\) on \(\tau_0\) admit a closed form, as trimmed means of \(P_{\mathsf{T}}\) net of the control mean \(\mu^{\mathsf{AO}}(0) = \mathbb{E}_{P_0}[Y \mid D = 0, S = 1]\), which is identified because, absent defiers, the control respondents are always-observed units. The bounds data frame carries one row per budget, with bootstrap standard errors and the outer percentile confidence interval:

head(fit_tv$bounds, 4)
#>   delta     lower     upper  lower_se  upper_se   ci_lower  ci_upper
#> 1  0.00 0.4649220 0.4649220 0.1232847 0.1232847 0.24236710 0.7092659
#> 2  0.05 0.3773408 0.5535522 0.1240051 0.1237240 0.14357117 0.7937557
#> 3  0.10 0.3084373 0.6223460 0.1251768 0.1247722 0.07599813 0.8581739
#> 4  0.15 0.2467927 0.6825040 0.1265714 0.1262810 0.01420249 0.9189803

plot(), an alias of tvbounds_plot(), shades the identified region, draws the outer confidence interval as dashed lines, marks the baseline point estimate at \(\delta = 0\), and flags the plug-in breakdown budget at which the lower bound crosses zero:

plot(fit_tv)

Total variation bounds for the simulated experiment.

2.3 The endpoints: baseline at \(\delta = 0\), Lee bounds at \(\delta = 1\)

The two endpoints of the budget grid anchor the interpretation. At \(\delta = 0\) the baseline is the only admissible distribution, the two unobserved outcome distributions coincide, \(P_{\mathsf{C}} = P_{\mathsf{AO}}\), and the bounds collapse to the difference in means among respondents \(\tau_{\mathsf{MCAR}}(P_0)\), the point estimate one would report under missingness completely at random:

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)
#>    point    naive 
#> 0.464922 0.464922

At \(\delta = 1\) the total variation restriction is vacuous, the robustness set imposes nothing beyond the mixture structure of the treated arm, and the bounds coincide with the trimming bounds of Lee (2009), \(\underline{\tau}_{\mathsf{Lee}}\) and \(\overline{\tau}_{\mathsf{Lee}}\), stored in details$lee:

subset(fit_tv$bounds, delta == 1, select = c(delta, lower, upper))
#>    delta      lower    upper
#> 21     1 -0.2010406 1.105157
unlist(fit_tv$details$lee)
#>      lower      upper 
#> -0.2010406  1.1051573

The exercise therefore interpolates between the two poles of the attrition literature, the point estimate that assumes attrition away and the bounds of Lee (2009) that restrict the two unobserved outcome distributions only through the mixture structure, and indexes the continuum between them by a budget that reads as a share of misspecified mass. Both endpoints are computed internally even when the delta grid omits them.

2.4 The contamination neighborhood

neighborhood = "contamination" replaces the total variation restriction with the one-sided mixture restriction \(P_{\mathsf{C}} = (1 - \delta) P_{\mathsf{AO}} + \delta R\) for some distribution \(R\), equivalently \(P_{\mathsf{C}} \ge (1 - \delta) P_{\mathsf{AO}}\) as measures. Since every such mixture lies within total variation distance \(\delta\) of \(P_{\mathsf{AO}}\), the contamination bounds are weakly tighter at every budget, while the two sets of bounds share both endpoints, the baseline at \(\delta = 0\) and the Lee bounds at \(\delta = 1\):

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
#> [1] TRUE
all(fit_ct$bounds$upper <= fit_tv$bounds$upper + 1e-12)
#> [1] TRUE

Because bounds is a plain data frame, the two neighborhoods are readily compared with ggplot2:

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")

Total variation versus contamination bounds.

Which neighborhood to report is a modeling choice. The total variation bounds guard against arbitrary misspecification of a \(\delta\)-fraction of the baseline mass, whereas the contamination bounds impose in addition that the remaining \((1-\delta)\)-share of the compliers is distributed exactly as the always-observed units.

2.5 Cluster bootstrap

Treatment is assigned at the village level, so resampling villages is the appropriate bootstrap. Passing cluster makes each replicate draw whole clusters with replacement. Every replicate also redraws the arm-specific response rates, so that the standard errors reflect the estimation uncertainty in the complier share \(\pi\), which a plug-in formula would omit.

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
#> [1] 50

data.frame(delta      = fit_tv$bounds$delta,
           se_iid     = fit_tv$bounds$lower_se,
           se_cluster = fit_cl$bounds$lower_se)[1:4, ]
#>   delta    se_iid se_cluster
#> 1  0.00 0.1232847  0.1656593
#> 2  0.05 0.1240051  0.1642646
#> 3  0.10 0.1251768  0.1642461
#> 4  0.15 0.1265714  0.1643545

As one would expect with village-level assignment, the clustered standard errors are markedly larger than those obtained by resampling units independently.

2.6 Covariates and pooled bounds

When covariates is supplied, units are stratified on the interaction of the covariate columns, which plays the role of a discrete covariate \(X\) with support \(\mathcal{X}\), and cells with fewer than min_obs observed outcomes in either arm are dropped with a warning. Within a cell the complier share \(\pi(x)\) and the outcome distribution of the treated respondents \(P_{\mathsf{T}}(x)\) are identified, and the cell-level construction is the one of Section 2.2.

The budget can then be allocated across cells in two ways. The within-stratum (“pointwise”) restriction imposes \(\mathsf{TV}(P_{\mathsf{C}}(x) \,\|\, P_{\mathsf{AO}}(x)) \le \delta\) in every cell separately, giving one robustness set \(\mathcal{Q}^{\mathsf{pw}}_{\mathsf{C}}(\delta; x)\) per cell, whereas the pooled restriction caps only the average departure, \[ \int_{\mathcal{X}} \mathsf{TV}(P_{\mathsf{C}}(x) \,\|\, P_{\mathsf{AO}}(x)) \, \mathrm{d}P_{X \mid \mathsf{AO}}(x) \le \delta , \] and so allows heterogeneity across cells inside the single robustness set \(\mathcal{Q}_{\mathsf{C},X}(\delta)\). Under the total variation neighborhood the reported bounds are the pooled (joint) bounds of the paper, \(\underline{\tau}_{X}(\delta)\) and \(\overline{\tau}_{X}(\delta)\), for which a single budget \(\delta\) is allocated optimally across the covariate cells rather than imposed cell by cell: the least favorable distributions concentrate the misspecified mass in the cells in which it moves the aggregate effect the most.

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
#>   stratum    weight    p_star n_treated_obs n_control_obs       mu0
#> 1       0 0.5703125 0.3266925           117            73 0.9145775
#> 2       1 0.4296875 0.3760314            96            55 1.5721558

The per-stratum table records the weight of each retained cell, that is its share of control respondents, which under monotonicity is the covariate distribution of the always-observed population, \(P_{X \mid D = 0, S = 1} = P_{X \mid \mathsf{AO}}\), together with the complier share and the sample sizes. The within-stratum reference curve \(\underline{\tau}^{\mathsf{pw}}_{X}(\delta)\), \(\overline{\tau}^{\mathsf{pw}}_{X}(\delta)\) is kept in details$pooled$pw. The pooled bounds contain it by construction, since a common budget in every cell is one of the ways of spending the average one, and the two coincide at both endpoints:

pw <- fit_x$details$pooled$pw
all(fit_x$bounds$upper >= pw$upper - 1e-12)
#> [1] TRUE
all(fit_x$bounds$lower <= pw$lower + 1e-12)
#> [1] TRUE

Conditioning on covariates tightens the worst case. At \(\delta = 1\) both covariate constructions collapse to the covariate Lee bounds \(\underline{\tau}_{\mathsf{Lee},X}\) and \(\overline{\tau}_{\mathsf{Lee},X}\), which aggregate cell-level trimming bounds and are tighter than the Lee bounds that ignore the covariate, stored in details$lee_nocov for reference:

rbind(with_covariate    = unlist(fit_x$details$lee),
      without_covariate = unlist(fit_x$details$lee_nocov))
#>                        lower    upper
#> with_covariate    -0.1705010 1.039245
#> without_covariate -0.2010406 1.105157

With neighborhood = "contamination" the covariate version imposes the common budget \(\delta\) within every retained cell and then aggregates, and those bounds remain weakly inside the total variation ones at every budget.

2.7 Summary measures

summary(), an alias of tvbounds_summary(), condenses the bounds into the summary measures of the paper. These are stated for generic bound paths \(\underline{\tau}(\delta)\) and \(\overline{\tau}(\delta)\), here the attrition bounds of Section 2.2, and are evaluated on the path adjacent to the reference value \(\tau_{\star}\), which is set through the argument tau_star and defaults to 0; here that is the lower path, since the baseline estimate is positive:

summary(fit_tv)
#> <tvbounds summary> treatment effect (attrition application, total-variation neighborhood)
#>   direction: lower bound path relative to tau_star = 0 (baseline point = 0.465, n = 500)
#>   breakdown budget: plug-in = 0.419; certified = 0.162; normal floor = 0.149
#>   at delta = 0.419: shadow price eta = 0.73; robustness SE varsigma = 4.15 (scale-free 0.186)
#>   certification frontier: n* = 998 at budget 0.162 + jump 0.05 (Delta n = +498, cost per pp = 4980)

We read the output line by line.

• Plug-in breakdown budget (delta_b) — the breakdown budget \(\delta_b(\tau_{\star})\), the smallest budget at which the estimated bounds reach \(\tau_{\star}\), linearly interpolated between grid points. Here the positive-effect conclusion survives until roughly 0.42: the outcome distributions of the compliers and of the always-observed units would have to disagree on more than 42% of their mass to overturn the sign. When the path never reaches \(\tau_{\star}\), the breakdown is censored: in view of the convention of the paper that the infimum over an empty set equals one, it is reported at the right endpoint of the budget grid, which is 1 for the total variation and contamination neighborhoods, with censored = TRUE rather than as “greater than one”.

• Certified breakdown budget (delta_b_ci) — the certified breakdown budget \(\widehat{\delta}_{b}^{\,\mathsf{C}}(\alpha)\), read off as the first crossing of the outer confidence limit, that is of the same band the figures draw, so that tables and figures agree on one number. It is the largest budget at which the conclusion is still certified at the level of the band and, since the band is wider than the estimated bounds, it is never larger than the plug-in breakdown: sampling uncertainty erodes the robustness one can certify.

• Normal floor (delta_b_ci_norm) — a diagnostic that recomputes the certified breakdown from the normal approximation \(\widehat{\underline{\tau}}_{n}(\delta) - z_{1-\alpha/2} \widehat{\sigma}_{n}(\delta) / \sqrt{n}\). It estimates the same population quantity as delta_b_ci and differs from it only when the bootstrap distribution of the bound is asymmetric.

• Shadow price (eta) — the shadow price of robustness \(\underline{\eta}(\delta) = -\underline{\tau}'(\delta)\), minus the slope of the bound path at the evaluation budget, which measures how fast the bound deteriorates per unit of misspecification. A large eta indicates that conclusions erode quickly.

• Robustness standard error (varsigma) — the quantity \(\varsigma_{b} = \sigma(\delta_b) / \underline{\eta}(\delta_b)\), computed as se * sqrt(n) / eta, which converts the sampling uncertainty of the bound from units of the estimand into budget units and thereby makes robustness comparable across studies; varsigma_sc is its finite-sample analogue se / eta, that is \(\varsigma_b / \sqrt{n}\).

• Certification frontier (n_star, delta_n, cost_per_pp) — the frontier \(n^{\star}(\delta;\alpha)\), the sample size at which the normal confidence limit would just certify the conclusion at the certified budget raised by jump (default 0.05). delta_n is the additional sample relative to the realized n, and cost_per_pp prices one percentage point of certified budget at cost_per_unit (default 50 USD per sampled unit, the paper’s benchmark) — a discrete version of the certification elasticity \(\mathrm{d}\log n^{\star}(\delta;\alpha) / \mathrm{d}\delta\), the paper’s exchange rate between observations and misspecified probability mass.

By default the shadow price and the robustness standard error are evaluated at the plug-in breakdown budget. Passing delta evaluates them at a budget of substantive interest instead, and the one-row data frame measures gives access to all of them programmatically:

s01 <- tvbounds_summary(fit_tv, delta = 0.1)
s01$measures[c("delta_eval", "eta", "varsigma", "n_star")]
#>   delta_eval     eta varsigma n_star
#> 1        0.1 1.30548 2.144069    506

3.1 A small formula-instrument design

Many causal questions in economics involve treatments or instruments that combine several sources of variation according to a known formula, shift-share designs being the leading example. Write \(y_i = \beta x_i + \varepsilon_i\) for the structural equation of unit \(i\), with \(\beta\) the parameter of interest, \(x_i\) the endogenous regressor and \(\varepsilon_i\) the unobserved residual; let \(v \in \mathcal{V}\) be the vector of exogenous shocks, \(w\) the predetermined covariates, and \(f_i(\cdot\,;w) : \mathcal{V} \to \mathbb{R}\) the known formula, so that \(z_i = f_i(v;w)\) is the candidate instrument at the realized shocks. Borusyak and Hull (2023) point out that exogeneity of the shocks does not, by itself, deliver a consistent instrumental variables estimator, because units are nonrandomly exposed to the shocks, and show that validity is restored once the instrument is recentered at the expected instrument \(\mu_i(P_{*}) = \mathbb{E}_{P_{*}}[f_i(v;w) \mid w]\), its average across the shocks that might as well have been realized under the assignment process \(P_{*}\) the researcher postulates, which gives the recentered instrument \(\widetilde{z}_i(P_{*}) = z_i - \mu_i(P_{*})\). In a randomized experiment the assignment process is given by the randomization protocol, but in natural experiments it is rarely known or justified by economic theory, and the sensitivity exercise of Palomba (2026) accordingly asks how much the estimate would move when up to a \(\delta\)-fraction of \(P_{*}\) is misspecified.

We simulate 150 regions exposed to \(K = 10\) sectors through the employment shares collected in W, so the formula is the shift-share \(f_i(v;w) = \sum_{k=1}^{K} w_{ik} v_k\) and the realized instrument is \(z_i = \sum_{k} w_{ik} v_k\) at the realized sector shocks g0. The researcher postulates that the shocks are i.i.d. standard normal and recenters using S = 80 counterfactual shock configurations \(v^{(1)}, \dots, v^{(S)}\) drawn from that postulated process, the columns of G. The matrix Fmat collects the formula evaluated at each configuration, so that Fmat[i, s] is \(f_i(v^{(s)};w)\), the instrument of region i under the s-th draw. The realized shocks are drawn with a small positive mean, so that the postulated mean-zero assignment process is itself mildly misspecified, and the outcome loads on a sectoral confounder, so that nonrandom exposure matters.

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

3.2 Bounds under both neighborhoods

tvbounds_riv() takes the outcome \(y_i\), the endogenous regressor \(x_i\), the realized instrument \(z_i\), and the matrix of counterfactual draws. Optionally, p supplies the probabilities that \(P_{*}\) attaches to the configurations, uniform by default, and controls collects the covariates to be partialled out of y, x, z, and every column of Fmat, in view of the Frisch–Waugh–Lovell theorem. As in the paper, all analyses are conducted conditionally on the realized sample, so that the bounds are deterministic and no standard errors or confidence intervals are reported.

riv_tv <- tvbounds_riv(y, x, z, Fmat, delta = seq(0, 1, by = 0.01))
riv_tv
#> <tvbounds> riv bounds under a total-variation neighborhood
#>   estimand: IV coefficient; baseline point estimate (delta = 0): 0.434; n = 150
#>   budget grid: 101 values of delta in [0, 1]
#>   inference: none attached
#>   Use summary() for breakdown and price-of-robustness measures; plot() to display.

riv_ct <- tvbounds_riv(y, x, z, Fmat, delta = seq(0, 1, by = 0.01),
                       neighborhood = "contamination")

The baseline point is the recentered instrumental-variables estimate \(\widehat{\beta}_{*} = \widehat{\beta}(P_{*})\), so that the exercise is centered at the reported estimate. The formula enters only through the two criteria \(g_y(\cdot) = \sum_{i=1}^{n} y_i f_i(\cdot\,;w)\) and \(g_x(\cdot) = \sum_{i=1}^{n} x_i f_i(\cdot\,;w)\), whose recentered values are the reduced form \(G_y(P) = g_y(v) - \mathbb{E}_P[g_y]\) and the first stage \(G_x(P) = g_x(v) - \mathbb{E}_P[g_x]\), so that \(\widehat{\beta}(P) = G_y(P) / G_x(P)\). The bounds trace the range of \(\widehat{\beta}(P)\) as the distribution over the S configurations varies either in the total variation ball \[ \mathcal{P}^{\mathsf{FI}}_{\mathsf{TV}}(\delta) := \{P \in \Delta(\mathcal{V}) : \mathsf{TV}(P, P_{*}) \le \delta\}, \] which delivers \(\underline{\beta}_{\mathsf{TV}}(\delta)\) and \(\overline{\beta}_{\mathsf{TV}}(\delta)\), or in the contamination neighborhood \[ \mathcal{P}^{\mathsf{FI}}_{\mathsf{cont}}(\delta) := \{P \in \Delta(\mathcal{V}) : P = \delta R + (1 - \delta) P_{*}, \; R \in \Delta(\mathcal{V})\}, \] which delivers \(\underline{\beta}_{\mathsf{cont}}(\delta)\) and \(\overline{\beta}_{\mathsf{cont}}(\delta)\). As in the attrition application, the contamination bounds are nested within the total variation bounds at every budget:

all(riv_ct$bounds$lower >= riv_tv$bounds$lower - 1e-10, na.rm = TRUE)
#> [1] TRUE

plot(riv_tv)

Total variation bounds for the recentered IV estimate.

The breakdown budget for the sign of \(\beta\), obtained at the reference value \(\tau_{\star} = 0\) set through the argument tau_star, is interior here, and summarizes in a single number how robust the design is. Since the exercise carries no inference, summary() reports the plug-in measures and returns NA for the certified ones, together with a note to that effect.

riv_tv$details$delta_breakdown
#> [1] 0.8426913
summary(riv_tv)
#> <tvbounds summary> IV coefficient (riv application, total-variation neighborhood)
#>   direction: lower bound path relative to tau_star = 0 (baseline point = 0.434, n = 150)
#>   breakdown budget: plug-in = 0.843; certified = ---; normal floor = ---
#>   at delta = 0.843: shadow price eta = 0.948; robustness SE varsigma = --- (scale-free ---)
#>   certification frontier: n* = --- at budget --- + jump 0.05 (Delta n = ---, cost per pp = ---)
#>   note: no confidence band is attached to the lower bound path, so the certified breakdown and the certification frontier are reported as NA (this application carries no inference by design). 
#>   note: no bootstrap standard errors are attached to the lower bound path, so the normal-floor diagnostic and the robustness standard error are reported as NA.

3.3 The first-stage breakdown budget and its censoring convention

A formula-instrument design can fail in a more fundamental way than by a change of sign. Once the budget is large enough that some distribution in the neighborhood makes the recentered first stage \(G_x(P) = \sum_{i=1}^{n} \widetilde{z}_i(P)\, x_i\) vanish, the estimate \(\widehat{\beta}(P)\) is no longer well defined over the whole neighborhood and the identified set is the entire real line. The smallest such budget is the first-stage breakdown budget, \(\delta^{\mathsf{TV}}_{\mathsf{FS}}\) under the total variation ball and \(\delta^{\mathsf{cont}}_{\mathsf{FS}}\) under the contamination neighborhood, reported, for the neighborhood in use, in details$delta_fs together with a censoring flag:

c(delta_fs = riv_tv$details$delta_fs,
  censored = riv_tv$details$delta_fs_censored)
#> delta_fs censored 
#>        1        1

Following the convention of the paper, the infimum over an empty set equals one, so that when the first stage never breaks down on \([0, 1]\), as is the case here, delta_fs is reported as 1 with delta_fs_censored = TRUE, which records that the first stage never breaks down over the budget range rather than that breakdown occurs at 1; the summary measures adopt the same convention for censored breakdown budgets. The first-stage budgets of both neighborhoods, \(\delta^{\mathsf{TV}}_{\mathsf{FS}}\) and \(\delta^{\mathsf{cont}}_{\mathsf{FS}}\), are always stored, each carrying its own censoring attribute, in details$delta_fs_tv and details$delta_fs_cont.

With a weaker first stage the breakdown becomes interior, and the rows of bounds beyond it carry NA, the bounds being vacuous there:

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)
#> delta_fs censored 
#> 0.604252 0.000000
tail(riv_weak$bounds, 3)                     # vacuous budgets are NA
#>     delta lower upper
#> 99   0.98    NA    NA
#> 100  0.99    NA    NA
#> 101  1.00    NA    NA

plot(riv_weak, breakdown = FALSE)

A weak first stage: the bounds diverge and become vacuous at the first-stage breakdown budget.

The figure makes the mechanics apparent: the bounds fan out as the budget approaches details$delta_fs and are not drawn beyond it. Reporting the first-stage breakdown budget alongside the bounds keeps two rather different fragilities distinct, a conclusion whose sign is overturned and a design whose identifying variation can be undone altogether.

4.1 Requirements and the first-use check

The solver is written in Julia and reached through the ‘JuliaCall’ package. One-time setup:

1. install the R bridge: install.packages("JuliaCall");
2. install Julia (>= 1.9) from https://julialang.org/downloads/ and make sure the julia binary is on the PATH (or set the JULIA_HOME environment variable to its bin/ directory);
3. install KNITRO and obtain a valid license — KNITRO is a commercial solver by Artelys, and free academic trials are available at https://www.artelys.com/solvers/knitro/; make sure KNITRO.jl can locate the installation (e.g. set the KNITRO_DIR environment variable before starting R).

On the first call in each R session the package initializes the embedded Julia, activates and instantiates its own Julia environment, downloading KNITRO.jl, ForwardDiff.jl, Optim.jl, and their dependencies into the standard Julia depot on first use, which requires network access. It then sources the solver module and verifies that KNITRO.jl loads and that a KNITRO solver context can be created, which exercises the license. On success a message reports that KNITRO is available; on failure the call stops and reports installation and license guidance. The check is performed once per R session.

Two session-level caveats are worth recording:

• The package activates its own Julia project inside the embedded session. If JuliaCall is used for other work in the same R session, that Julia code will run with the project environment of tvbounds active from the first tvbounds_counterfactual() call onward.
• Julia’s global random-number generator is re-seeded by the solver, since the multi-start restarts are deterministic by design. The random-number state of R is saved and restored around each call, so R streams are unaffected, but the continuity of the global stream of Julia should not be relied upon across a call.

4.2 The moments function: the Julia route

The model is communicated through a single moments function with the in-place signature moments!(K, G, theta, U, obj): given the current parameter theta and the M x u_dim matrix of latent draws U, it fills the M-vector K with the counterfactual values \(g(U^{(j)};\theta)\) and the M x d matrix G with the moment functions \(m(U^{(j)};\theta)\), one row per draw \(U^{(j)}\), \(j = 1, \dots, M\). The paper writes the counterfactual integrand \(g\) and the moment function \(m\) while the code arrays are named K and G, so read K as \(g\) and G as \(m\) throughout. The package ships a toy example in inst/julia/examples/toy.jl — the latent draw is a scalar \(U \sim \mathrm{Uniform}(0,1)\), the single moment pins its mean to \(\theta\), and the counterfactual is \(U\) itself:

function tvb_toy_moments!(K, G, theta, U, obj)
    M = size(U, 1)
    @inbounds for m in 1:M
        K[m]    = U[m, 1]
        G[m, 1] = U[m, 1] - theta[1]
    end
    return nothing
end

For your own model, put a function with the same signature in a .jl file. The obj.gamma field carries an arbitrary payload passed from R through the gamma argument (an R list arrives in Julia as an ordered dictionary keyed by symbols, so read entries as obj.gamma[:name]), and the scrambled-Halton draws are uniform on the unit cube \((0,1)^{\mathsf{d}_z}\), whose dimension \(\mathsf{d}_z\) is set by u_dim, so the moments function is responsible for mapping them into baseline draws, typically through quantile transforms:

# mymodel.jl -- a template
using Distributions   # available in the package's Julia environment

function my_moments!(K, G, theta, U, obj)
    M = size(U, 1)
    a = obj.gamma[:a]                       # payload passed from R
    @inbounds for m in 1:M
        u1 = quantile(Normal(), U[m, 1])    # map uniforms into N(0,1) draws
        K[m]    = exp(theta[1] + theta[2] * u1)       # counterfactual
        G[m, 1] = u1 * (theta[1] + theta[2] * u1 - a) # moment condition 1
        G[m, 2] = theta[1] + theta[2] * u1 - a        # moment condition 2
    end
    return nothing
end

Write the function generically in the element type of theta (as above: no Float64 annotations on the values it computes): the outer optimization differentiates the moments in \(\theta\) by automatic differentiation (ForwardDiff), which feeds dual numbers through your code. For a non-generic function, pass gradient = "fd" to fall back on finite differences.

4.3 Running the solver

The complete call for the toy model (not run — Julia and a licensed KNITRO are required):

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)

For this toy model the bounds under "TVmix" equal the endpoints of the \(\theta\) box (up to the smoothing described below) whenever the moment condition is satisfiable, which is what the package’s integration tests check.

Beyond the moments specification, the main arguments are the following.

• d is the number of moment conditions \(\mathsf{d}_m\); the common length of theta_lb/theta_ub fixes the dimension \(\mathsf{d}_\theta\) of \(\theta\), and the outer problems optimize over that box. A degenerate box (theta_lb == theta_ub) skips the outer optimization and reports fixed-\(\theta\) bounds at theta_init — useful for pure distribution-sensitivity exercises.
• side selects which problems to solve per budget: "both" (default), "lower", or "upper".
• U supplies the M x u_dim matrix of latent draws \(U^{(1)}, \dots, U^{(M)}\) directly; when U = NULL the package generates M scrambled-Halton uniforms (Owen, 2017) of dimension u_dim (at most 15), seeded by seed. The paper’s default is M = 50000.
• theta_init should be set to the baseline estimate of the model: the reported baseline point is the plug-in counterfactual \(\mathsf{k}(\theta;P_{*}) = \mathbb{E}_{P_{*}}[g(U;\theta)]\) at theta_init and is only returned when theta_init is supplied. The budget-zero baseline is not part of the delta grid (all budgets must be strictly positive); it is exactly this point.
• Failed budgets are reported as NA in bounds (for "TVmix" an NA typically signals an empty robustness set at that budget).

Divergences. The divergence keyword selects the entropy function \(\phi\) whose divergence \(D_\phi(P \,\|\, P_{*})\) defines the neighborhood; the budget must lie in \((0, 1]\) for the total-variation family and may be any positive number for the others:

The kinked total-variation conjugates \(\phi^{*}_{\mathsf{TV}}\) are Huber-smoothed and the per-draw maxima log-sum-exp-smoothed; both smoothings lie above the exact functions, so computed bounds are outward-conservative (wider, never narrower) at order \(10^{-3}\).

4.4 The moments function: the R route

The model may alternatively be supplied as an R function of the form function(theta, U, gamma) returning list(K = <numeric length M>, G = <numeric M x d matrix>):

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)

Two caveats are worth recording. First, the R route is much slower: every objective evaluation crosses the R/Julia boundary. Second, ForwardDiff cannot differentiate through R code, so the outer optimization uses finite differences by default; alternatively, supply gradient as an R function function(theta, U, gamma) returning either the stacked (M*(d+1)) x l Jacobian of c(K, G) (the K rows first, then G in column-major order) or a list with components K (M x l) and G (M x d x l), where l is the dimension \(\mathsf{d}_\theta\) of \(\theta\). For sustained use we recommend porting the moments to a small Julia file; the template above is complete.

4.5 Tuning with tvbounds_control()

tvbounds_control() collects the solver tuning options, and its defaults reproduce the settings of the paper. The options most commonly adjusted are the following.

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)

• maxsolves restarts the outer optimization over \(\theta\) from perturbed starting points (startptrange controls the perturbation size); each budget is warm-started at the previous budget’s optimum, and the reported bound is the inner dual value re-solved at the best candidate, which makes the bound curves monotone in the budget by construction.
• inner_opt / outer_opt point at KNITRO option files for the inner (dual) and outer problems. The package ships four under system.file("julia", "opt", package = "tvbounds"): the defaults inner.opt and outer.opt, plus outer_fast.opt (analytic envelope gradient with looser tolerances, suited to plotting grids) and outer_boot_tv.opt.
• knitro_options overrides individual KNITRO options by name (e.g. list(maxit = 500, outlev = 2)); the overrides are merged into both option files by writing merged copies to tempdir() — the shipped files are never modified.
• use_optim = TRUE replaces the outer KNITRO solve with Optim.jl (projected L-BFGS with the analytic envelope-theorem gradient), which avoids nested KNITRO contexts — these segfault with some KNITRO.jl versions; time_limit, iterations, and outer_iterations govern that fallback.
• The remaining options (eta_min, lower_limit, psi_tv_eps, tvac_tau, tvmix_tau, purekl_acap, tvmix_kappa) expose the smoothing and guard constants of the divergence machinery; see ?tvbounds_control for the full list.

Diagnostics are collected in details$solver (per-budget outer flags, inner KNITRO status codes, and timings for each side), details$theta_lower / details$theta_upper (the optimal structural parameters at each budget), and details$control (the resolved control list, including the option files actually used).

Citation

If you use tvbounds in your work, please cite the paper:

Palomba, F. (2026). “Sensitivity Analysis in Population Shares.” Working paper.

citation("tvbounds") produces a citation entry for the package itself.