---
title: "Sensitivity Analysis with tvbounds"
author: "Filippo Palomba"
output:
rmarkdown::html_vignette:
toc: true
vignette: >
%\VignetteIndexEntry{Sensitivity Analysis with tvbounds}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r setup-knitr, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.width = 6.5,
fig.height = 4,
fig.align = "center",
out.width = "92%"
)
```
```{r setup, message = FALSE}
library(tvbounds)
library(ggplot2)
```
## 1. The framework in a nutshell
Empirical conclusions often rest on untestable distributional assumptions:
that attrition in a randomized experiment is as good as random, that the
exogenous shocks in a formula instrument were assigned by the mechanism
the design postulates, or that the preference shocks in a discrete choice model
follow an extreme-value distribution. Each of them amounts to a choice of
distribution for an object the data does not identify, and a stakeholder may
reasonably ask by how much the reported conclusion would change were that
choice wrong. **tvbounds** implements the sensitivity analysis of Palomba
(2026), which answers this question and computes sensitivity bounds, for
randomized experiments with attrition, counterfactuals in structural models,
and recentered (formula) instrumental variables.
The exercise rests on three ingredients:
* an **estimand** $\mathbb{E}_P[g(Z;\theta)]$, the expectation of a known
function $g$ of a random vector $Z$ under a distribution $P$, indexed by
a structural parameter $\theta \in \Theta$ in the applications that carry
one;
* a **baseline distribution** $P_{*}$, the distribution the econometrician
would use absent any concern for misspecification;
* a **budget** $\delta \in [0, 1]$, the sensitivity parameter, which caps
the discrepancy from the baseline we are willing to entertain.
The candidate distributions the package optimizes over form the
**robustness set**
$$
\mathcal{P}_\phi(\theta; \rho, P_{*}, \delta)
= \{P : \mathbb{E}_P[m(Z;\theta)] \in \mathcal{M}(\rho), \;\;
D_\phi(P \,\|\, P_{*}) \le \delta\},
$$
which collects the distributions that remain compatible with the model,
where $m$ is the moment function and $\mathcal{M}(\rho)$ the moment
constraint set at the target $\rho$, and that lie within $\delta$ of the
baseline in the divergence $D_\phi$ selected by an entropy function $\phi$.
We abbreviate it $\mathcal{P}_\phi(\theta;\delta)$ whenever the target and
the baseline are clear from the context. At $\delta = 0$ the baseline is the
only admissible distribution.
Our preferred divergence is the one generated by the total variation entropy
$\phi_{\mathsf{TV}}(s) = |s - 1| / 2$, for which
$D_{\phi_{\mathsf{TV}}}(P \,\|\, P_{*}) = \mathsf{TV}(P, P_{*})$. A
distribution within total variation distance $\delta$ of the baseline agrees
with it except on at most a $\delta$-fraction of the baseline probability
mass, so that $\delta$ reads as the fraction of the mass in $P_{*}$ that may
be misspecified. The second neighborhood available in the package is the
*contamination* neighborhood in the tradition of Huber (1964), which
collects the distributions of the mixture form
$P = (1-\delta) P_{*} + \delta R$ for an arbitrary distribution $R$,
equivalently $P \ge (1-\delta) P_{*}$ as measures. It follows that $(1-\delta)$-share of
the population is distributed exactly as the baseline and only the remaining
$\delta$-share is left arbitrary. Since every such mixture lies within total variation
distance $\delta$ of the baseline, the contamination bounds are weakly
tighter at every budget than the total variation one.
Given the estimand, the baseline, and the budget, the package reports
**sensitivity bounds**
$$
\Big[\;\inf_{\theta \in \Theta} \; \inf_{P \in \mathcal{P}_\phi(\theta;\delta)} \mathbb{E}_P[g(Z;\theta)],\;\;
\sup_{\theta \in \Theta} \; \sup_{P \in \mathcal{P}_\phi(\theta;\delta)} \mathbb{E}_P[g(Z;\theta)]\;\Big],
$$
the extreme values the estimand can take over the robustness set. They are
nested extrema of the same sign rather than a minimax, and they trace the
exact range of the estimand: every value inside the bounds is attained by
some distribution in the robustness set, and no value outside is. The
exercise reports the whole path of the bounds as the budget varies over the
range of the chosen divergence. At $\delta = 0$ the interval collapses to
the baseline value of the estimand and, as $\delta$ grows, it widens
monotonically; the budget at which the bounds first reach a
reference value $\tau_{\star}$ (zero, say, at which the sign of the effect
is no longer identified) is the **breakdown budget**
$\delta_b(\tau_{\star})$, the direct answer to the question of how much
misspecification would overturn the conclusion.
The package provides three interfaces, one per empirical application of the
paper.
* **Randomized experiments with attrition** —
`tvbounds_attrition()` bounds the average treatment effect when outcomes
are recorded only for respondents and attrition may be selective. At
$\delta = 0$ the bounds collapse to the baseline difference in means among
respondents and at $\delta = 1$ they coincide with the trimming bounds of
Lee (2009). Bootstrap standard errors and confidence intervals, optionally
with cluster resampling, and covariate-pooled bounds are also available.
* **Recentered instrumental variables** — `tvbounds_riv()` bounds a
recentered (formula) instrumental-variables estimate in the design of
Borusyak and Hull (2023), where identification rests on a distribution for
the shocks that the researcher postulates. The bounds trace how far the
estimate can move when up to a $\delta$-fraction of the postulated
distribution is misspecified, and a first-stage breakdown budget records
the budget at which the design itself collapses. The exercise is
conducted conditionally on the realized sample and carries no inference by
design.
* **Counterfactual predictions in structural models** —
`tvbounds_counterfactual()` bounds a counterfactual of a structural model
when the distribution of the latent variables ranges over a divergence
neighborhood of the simulated baseline $P_{*}$, as in
Christensen and Connault (2023). It is the only function that supports
divergences beyond total variation and contamination, and it relies on
Julia and on the commercial Artelys KNITRO solver (Section 4).
All three return a common `tvbounds` object with `print()`, `plot()`, and
`summary()` methods. `tvbounds_summary()` computes the summary measures of
the paper, namely the breakdown budgets, the shadow price of robustness, the
robustness standard error, and the certification frontier, and
`tvbounds_plot()` draws the bounds against the budget.
## 2. Randomized experiments with attrition
Consider a randomized experiment that grants a treatment to half of the
applicants and then runs a follow-up survey to record their outcomes.
Outcomes are observed only for the applicants who respond, and respondents
may be selected on dimensions that also drive the outcome, confounding the
effects of the program. Practitioners typically assume that attrition is
ignorable or, equivalently, that nonresponse is as good as random, an
assumption known as missingness completely at random. That assumption is
untestable, and the corroborating evidence usually offered, that attrition
is not differential across arms, is hardly conclusive. A stakeholder might
accordingly be interested in quantifying by how much the treatment effect
can change when it fails, which is what `tvbounds_attrition()` does.
### 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.
```{r 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
```
By construction the treatment effect equals 0.35, whereas only roughly
`r round(100 * mean(rct$s[rct$d == 0]))`%
of control units respond against `r round(100 * mean(rct$s[rct$d == 1]))`% 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.
```{r 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
```
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`:
```{r rct-pstar}
fit_tv$details$p_star
```
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:
```{r rct-bounds-head}
head(fit_tv$bounds, 4)
```
`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:
```{r rct-plot, fig.cap = "Total variation bounds for the simulated experiment."}
plot(fit_tv)
```
### 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:
```{r 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)
```
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`:
```{r rct-endpoint-1}
subset(fit_tv$bounds, delta == 1, select = c(delta, lower, upper))
unlist(fit_tv$details$lee)
```
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$:
```{r 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)
```
Because `bounds` is a plain data frame, the two neighborhoods are readily
compared with `ggplot2`:
```{r 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")
```
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.
```{r 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, ]
```
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.
```{r 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
```
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:
```{r 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)
```
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:
```{r rct-lee-cov}
rbind(with_covariate = unlist(fit_x$details$lee),
without_covariate = unlist(fit_x$details$lee_nocov))
```
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:
```{r rct-summary}
summary(fit_tv)
```
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
`r round(summary(fit_tv)$measures$delta_b, 2)`: the outcome distributions
of the compliers and of the always-observed units would have to disagree
on more than
`r round(100 * summary(fit_tv)$measures$delta_b)`% 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:
```{r rct-summary-delta}
s01 <- tvbounds_summary(fit_tv, delta = 0.1)
s01$measures[c("delta_eval", "eta", "varsigma", "n_star")]
```
## 3. Recentered instrumental variables
### 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.
```{r 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
```
### 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.
```{r 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")
```
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:
```{r riv-nesting}
all(riv_ct$bounds$lower >= riv_tv$bounds$lower - 1e-10, na.rm = TRUE)
```
```{r riv-plot, fig.cap = "Total variation bounds for the recentered IV estimate."}
plot(riv_tv)
```
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.
```{r riv-summary}
riv_tv$details$delta_breakdown
summary(riv_tv)
```
### 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:
```{r riv-fs}
c(delta_fs = riv_tv$details$delta_fs,
censored = riv_tv$details$delta_fs_censored)
```
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:
```{r 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
```
```{r 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)
```
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. Counterfactual predictions in structural models
Counterfactual predictions from structural models depend on the distribution
of the latent variables, which is typically not nonparametrically
identifiable and is often chosen for tractability alone; taste shocks,
unobserved heterogeneity, and productivity shocks are common instances. A
stakeholder may thus be interested in knowing how sensitive such predictions
are to that untestable assumption. The third application bounds a
counterfactual
$\mathbb{E}_P[g(U;\theta)]$ of a structural model when the distribution $P$
of the latent variables $U$ ranges over a divergence neighborhood of the
simulated baseline $P_{*}$ and the structural parameter $\theta \in \Theta$
ranges over the values compatible with the moment conditions
$\mathbb{E}_P[m(U;\theta)] \in \mathcal{M}(\rho)$, following Christensen and
Connault (2023) and Palomba (2026). The reported bounds are the nested
extrema
$$
\underline{\mathsf{k}}(\delta)
= \inf_{\theta \in \Theta} \; \inf_{P \in \mathcal{P}_\phi(\theta;\delta)}
\mathbb{E}_P[g(U;\theta)],
\qquad
\overline{\mathsf{k}}(\delta)
= \sup_{\theta \in \Theta} \; \sup_{P \in \mathcal{P}_\phi(\theta;\delta)}
\mathbb{E}_P[g(U;\theta)] .
$$
`tvbounds_counterfactual()` is the only function in the package that
supports divergences beyond total variation and contamination.
**None of the code chunks in this section is executed when the vignette is
built**, since the function relies on Julia (>= 1.9) *and on a licensed copy
of the commercial Artelys KNITRO solver*, neither of which can be assumed on
the machine building the vignette. The code is nevertheless complete and
runs as it stands on a machine with a working installation.
### 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 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
; 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:
```julia
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:
```julia
# 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):
```{r 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)
```
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:
| keyword | neighborhood |
|------------|-------------------------------------------------------------------------------|
| `"KL_chi2"`| hybrid Kullback–Leibler/chi-square divergence of Christensen and Connault (2023); the default |
| `"KL"` | pure Kullback–Leibler divergence, with entropy $\phi_{\mathsf{KL}}(s) = s \log s - s + 1$ |
| `"chi2"` | Pearson chi-square divergence |
| `"TV"` | total variation, with entropy $\phi_{\mathsf{TV}}(s) = \lvert s - 1 \rvert / 2$, so the budget constraint is $\mathsf{TV}(P, P_{*}) \le \delta$ |
| `"TVmix"` | total variation intersected with the mixture (contamination) constraint $P \ge \kappa P_{*}$ at mixing weight $\kappa = 1 - \delta$, solved in its exact reduced form |
| `"TVmixC"` | the literal dual of the same program — slower, kept as a cross-check, and the only mode supporting a mixing weight $\kappa < 1 - \delta$ (set `tvmix_kappa` in the control list) |
| `"TVac"` | total variation restricted to distributions absolutely continuous with respect to the baseline, $P \ll P_{*}$ |
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 = , G = )`:
```{r 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)
```
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.
```{r 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)
```
* `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).
## 5. From the paper to the package
The table below maps each empirical application of Palomba (2026) to its
package interface.
| Application in the paper | Exemplar | Package call |
|---|---|---|
| Randomized experiments with attrition: bounds on the treatment effect for the always-observed, bootstrap inference, covariate-pooled bounds | Lee (2009) | `tvbounds_attrition()` |
| Recentered instrumental variables / formula instruments: bounds on the IV coefficient, first-stage breakdown budget | Borusyak and Hull (2023) | `tvbounds_riv()` |
| Counterfactual predictions in structural models: divergence bounds through Julia/KNITRO | Christensen and Connault (2023) | `tvbounds_counterfactual()` |
| Summary measures: breakdown budgets, shadow price of robustness, robustness standard error, certification frontier | — | `tvbounds_summary()` / `summary()` |
| Sensitivity figures: bounds against the budget | — | `tvbounds_plot()` / `plot()` |
### 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.
## References
Borusyak, K. and Hull, P. (2023). "Nonrandom Exposure to Exogenous
Shocks." *Econometrica*, 91(6), 2155–2185.
Christensen, T. and Connault, B. (2023). "Counterfactual Sensitivity and
Robustness." *Econometrica*, 91(1), 263–298.
Huber, P. J. (1964). "Robust Estimation of a Location Parameter." *Annals
of Mathematical Statistics*, 35(1), 73–101.
Lee, D. S. (2009). "Training, Wages, and Sample Selection: Estimating
Sharp Bounds on Treatment Effects." *Review of Economic Studies*, 76(3),
1071–1102.
Owen, A. B. (2017). "A randomized Halton algorithm in R."
arXiv:1706.02808.
Palomba, F. (2026). "Sensitivity Analysis in Population Shares." Working
paper.