W_{j}}\frac{\sigma_{Y}}{\sigma_{W_{j}}}\,\beta_{W_{j}\sim D}.$

By applying the definition of the $R^{2}$-value, we derive

$\beta_{Y\sim D|X}-\beta_{Y\sim D|X,W}=-\beta_{Y\sim D|X,W}\,R^{2}_{D\sim
W}+\sum_{j=1}^{l}R_{Y\sim
W_{j}}\frac{\sigma_{Y}}{\sigma_{W_{j}}}\,\beta_{W_{j}\sim D}.$

Next, we use rule [ii] of the $R^{2}$-calculus – independent additivity – on
$R^{2}_{D\sim W}$ and rewrite $\beta_{Y\sim D|X,W}$ in terms of $R$-values and
$\sigma$-values, i.e. standard deviations:

$\begin{multlined}\beta_{Y\sim D|X}-\beta_{Y\sim D|X,W}\\\ \begin{aligned}
&\stackrel{{\scriptstyle\text{[ii]}}}{{=}}-R_{Y\sim D|W}\frac{\sigma_{Y\sim
W}}{\sigma_{D\sim W}}\sum_{j=1}^{l}R^{2}_{D\sim W_{j}}+\sum_{j=1}^{l}R_{Y\sim
W_{j}}\frac{\sigma_{Y}}{\sigma_{W_{j}}}\,\beta_{W_{j}\sim D}\\\
&=\sum_{j=1}^{l}\beta_{W_{j}\sim D}\bigg{[}R_{Y\sim
W_{j}}\frac{\sigma_{Y}}{\sigma_{W_{j}}}-\frac{\sigma_{D}}{R_{D\sim
W_{j}}\sigma_{W_{j}}}R_{Y\sim D|W}\frac{\sigma_{Y\sim W}}{\sigma_{D\sim
W}}R^{2}_{D\sim W_{j}}\bigg{]}\end{aligned}\end{multlined}\beta_{Y\sim
D|X}-\beta_{Y\sim D|X,W}\\\ \begin{aligned}
&\stackrel{{\scriptstyle\text{[ii]}}}{{=}}-R_{Y\sim D|W}\frac{\sigma_{Y\sim
W}}{\sigma_{D\sim W}}\sum_{j=1}^{l}R^{2}_{D\sim W_{j}}+\sum_{j=1}^{l}R_{Y\sim
W_{j}}\frac{\sigma_{Y}}{\sigma_{W_{j}}}\,\beta_{W_{j}\sim D}\\\
&=\sum_{j=1}^{l}\beta_{W_{j}\sim D}\bigg{[}R_{Y\sim
W_{j}}\frac{\sigma_{Y}}{\sigma_{W_{j}}}-\frac{\sigma_{D}}{R_{D\sim
W_{j}}\sigma_{W_{j}}}R_{Y\sim D|W}\frac{\sigma_{Y\sim W}}{\sigma_{D\sim
W}}R^{2}_{D\sim W_{j}}\bigg{]}\end{aligned}$

In order to extract the factor $\sigma_{Y\sim D+W_{-j}}/\sigma_{W_{j}\sim
D+W_{-j}}$, we apply rule [iii] – decomposition of unexplained variance – six
times and arrive at

$\begin{multlined}\beta_{Y\sim D|X}-\beta_{Y\sim D|X,W}\\\
\stackrel{{\scriptstyle\text{[iii]}}}{{=}}\sum_{j=1}^{l}\beta_{W_{j}\sim
D}\,\frac{\sigma_{Y\sim D+W_{-j}}}{\sigma_{W_{j}\sim
D+W_{-j}}}\Bigg{[}R_{Y\sim W_{j}}\sqrt{\frac{1-R^{2}_{W_{j}\sim
D+W_{-j}}}{1-R^{2}_{Y\sim D+W_{-j}}}}\\\ \qquad\qquad\qquad-R_{Y\sim
D|W}\,R_{D\sim W_{j}}\sqrt{\frac{(1-R^{2}_{W_{j}\sim D+W_{-j}})(1-R^{2}_{Y\sim
W_{j}|W_{-j}})}{(1-R^{2}_{D\sim W})(1-R^{2}_{Y\sim
D|W_{-j}})}}\,\Bigg{]}.\end{multlined}\beta_{Y\sim D|X}-\beta_{Y\sim D|X,W}\\\
\stackrel{{\scriptstyle\text{[iii]}}}{{=}}\sum_{j=1}^{l}\beta_{W_{j}\sim
D}\,\frac{\sigma_{Y\sim D+W_{-j}}}{\sigma_{W_{j}\sim
D+W_{-j}}}\Bigg{[}R_{Y\sim W_{j}}\sqrt{\frac{1-R^{2}_{W_{j}\sim
D+W_{-j}}}{1-R^{2}_{Y\sim D+W_{-j}}}}\\\ \qquad\qquad\qquad-R_{Y\sim
D|W}\,R_{D\sim W_{j}}\sqrt{\frac{(1-R^{2}_{W_{j}\sim D+W_{-j}})(1-R^{2}_{Y\sim
W_{j}|W_{-j}})}{(1-R^{2}_{D\sim W})(1-R^{2}_{Y\sim D|W_{-j}})}}\,\Bigg{]}.$
(32)

We concentrate on the term in brackets, denoted by $T_{j}$. Invoking rule [v]
– reduction of partial correlation – and the (conditional) independence
assumption, we infer

$\displaystyle R_{Y\sim W_{j}}\stackrel{{\scriptstyle\text{[v]}}}{{=}}R_{Y\sim
W_{j}|W_{-j}}\sqrt{1-R^{2}_{Y\sim W_{-j}}},\quad R_{D\sim
W_{j}}\stackrel{{\scriptstyle\text{[v]}}}{{=}}R_{D\sim
W_{j}|W_{-j}}\sqrt{1-R^{2}_{D\sim W_{-j}}},$ $\displaystyle R^{2}_{W_{j}\sim
D+W_{-j}}=R^{2}_{W_{j}\sim D|W_{-j}}.$

We insert these relationships into the expression of $T_{j}$ and simplify it
via rule [iii]. Then, we apply rule [iv] – recursion of partial correlation –
on $R_{Y\sim D|W}$ and simplify the resulting expression

$\displaystyle T_{j}$ $\displaystyle=\begin{multlined}R_{Y\sim
W_{j}|W_{-j}}\sqrt{1-R^{2}_{W_{j}\sim D|W_{-j}}}\sqrt{\frac{1-R^{2}_{Y\sim
W_{-j}}}{1-R^{2}_{Y\sim D+W_{-j}}}}\\\ -R_{Y\sim D|W}\,R_{D\sim
W_{j}|W_{-j}}\sqrt{1-R^{2}_{D\sim W_{-j}}}\sqrt{\frac{(1-R^{2}_{W_{j}\sim
D|W_{-j}})(1-R^{2}_{Y\sim W_{j}|W_{-j}})}{(1-R^{2}_{D\sim W})(1-R^{2}_{Y\sim
D|W_{-j}})}}\end{multlined}R_{Y\sim W_{j}|W_{-j}}\sqrt{1-R^{2}_{W_{j}\sim
D|W_{-j}}}\sqrt{\frac{1-R^{2}_{Y\sim W_{-j}}}{1-R^{2}_{Y\sim D+W_{-j}}}}\\\
-R_{Y\sim D|W}\,R_{D\sim W_{j}|W_{-j}}\sqrt{1-R^{2}_{D\sim
W_{-j}}}\sqrt{\frac{(1-R^{2}_{W_{j}\sim D|W_{-j}})(1-R^{2}_{Y\sim
W_{j}|W_{-j}})}{(1-R^{2}_{D\sim W})(1-R^{2}_{Y\sim D|W_{-j}})}}$
$\displaystyle\stackrel{{\scriptstyle\text{[iii]}}}{{=}}R_{Y\sim
W_{j}|W_{-j}}\sqrt{\frac{1-R^{2}_{D\sim W_{j}|W_{-j}}}{1-R^{2}_{Y\sim
D|W_{-j}}}}-R_{Y\sim D|W}R_{D\sim W_{j}|W_{-j}}\sqrt{\frac{1-R^{2}_{Y\sim
W_{j}|W_{-j}}}{1-R^{2}_{Y\sim D|W_{-j}}}}$
$\displaystyle\stackrel{{\scriptstyle\text{[iv]}}}{{=}}\\!R_{Y\sim
W_{j}|W_{-j}}\sqrt{\frac{1-R^{2}_{D\sim W_{j}|W_{-j}}}{1-R^{2}_{Y\sim
D|W_{-j}}}}\\!-\\!R_{D\sim W_{j}|W_{-j}}\frac{R_{Y\sim
D|W_{-j}}\\!\\!-\\!R_{Y\sim W_{j}|W_{-j}}R_{D\sim
W_{j}|W_{-j}}}{\sqrt{1-R^{2}_{Y\sim D|W_{-j}}}\sqrt{1-R^{2}_{D\sim
W_{j}|W_{-j}}}}$ $\displaystyle=\frac{R_{Y\sim W_{j}|W_{-j}}(1-R^{2}_{D\sim
W_{j}|W_{-j}})-R_{Y\sim D|W_{-j}}R_{D\sim W_{j}|W_{-j}}-R_{Y\sim
W_{j}|W_{-j}}R^{2}_{D\sim W_{j}|W_{-j}}}{\sqrt{1-R^{2}_{Y\sim
D|W_{-j}}}\sqrt{1-R^{2}_{D\sim W_{j}|W_{-j}}}}$ $\displaystyle=\frac{R_{Y\sim
W_{j}|W_{-j}}-R_{Y\sim D|W_{-j}}\,R_{W_{j}\sim D|W_{-j}}}{\sqrt{1-R^{2}_{Y\sim
D|W_{-j}}}\sqrt{1-R^{2}_{W_{j}\sim D|W_{-j}}}}=R_{Y\sim W_{j}|D,W_{-j}}.$

Returning to equation (32), we plug in $T_{j}=R_{Y\sim W_{j}|D,W_{-j}}$ and
thus finish the proof

$\displaystyle\beta_{Y\sim D|X}-\beta_{Y\sim D|X,W}$
$\displaystyle=\sum_{j=1}^{l}\beta_{W_{j}\sim D}\,\frac{\sigma_{Y\sim
D+W_{-j}}}{\sigma_{W_{j}\sim D+W_{-j}}}R_{Y\sim W_{j}|D,W_{-j}}$
$\displaystyle=\sum_{j=1}^{l}\beta_{Y\sim W_{j}|D,W_{-j}}\,\beta_{W_{j}\sim
D}.$

∎

Lemma 6 helps us express the bias of OLS and $k$-class estimands in terms of
partial $R$-values which serve as sensitivity parameters. Whether these are
intuitive, depends on the causal structure of the underlying DAG. In the case
of two independent unmeasured variables $U_{1}$ and $U_{2}$ which confound or
mediate $\beta$ – the direct effect of $D$ on $Y$ –, the sensitivity
parameters $(R_{D\sim U_{1}},R_{D\sim U_{2}})$ and $(R_{Y\sim
U_{1}|D,U_{2}},R_{Y\sim U_{2}|D,U_{1}})$ are indeed intuitive. The former
tuple targets the dependence between $D$ and $U$, the latter tuple focuses on
the direct effects of $U$ on $Y$ regressing out the remaining variables.

The following theorem demonstrates how Lemma 6 can be applied. We identify the
bias of the $k$-class estimand in terms of the intuitive sensitivity
parameters.

######  Theorem 2.

Assume the setting of Theorem 1 and let $U=(U_{1},U_{2})$ be a two-dimensional
random vector. Further, suppose $U_{1}\perp\kern-5.38193pt\perp U_{2}\,|\,X,Z$
holds. Then,

$\beta_{k}-\beta=\Bigg{[}\frac{f_{Y\sim Z|X,D}\,R_{D\sim
Z|X}}{1-k+k\,R^{2}_{D\sim Z|X}}\\\
+\sum_{j=1}^{2}\frac{R_{j}\,f_{j}}{\sqrt{1-f^{2}_{j}\,f^{2}_{-j}+\left(R_{-j}\sqrt{\frac{1-R^{2}_{j}}{1-R^{2}_{-j}}}-R_{j}\,f_{j}\,f_{-j}\right)^{2}}}\Bigg{]}\frac{\sigma_{Y\sim
X+Z+D}}{\sigma_{D\sim X+Z}},$

where $R_{j}$ and $f_{j}$ abbreviate $R_{Y\sim U_{j}|X,Z,D,U_{-j}}\\!$ and
$f_{D\sim U_{j}|X,Z}$, respectively, for ${j\\!\in\\!\\{1,2\\}}$.

###### Proof.

Similarly to the proof of Theorem 1, we expand the difference
$\beta_{k}-\beta$ as a telescoping sum

$\beta_{k}-\beta=(\beta_{k}-\beta_{Y\sim D|X,Z})+(\beta_{Y\sim
D|X,Z}-\beta_{Y\sim D|X,Z,U}),$

which allows us to deal with the two summands separately. Following from the
same arguments, the first summand equals

$\frac{1}{1-k(1-R^{2}_{D\sim Z|X})}(\beta_{Y\sim D|X}-\beta_{Y\sim
D|X,Z})=\frac{f_{Y\sim Z|X,D}\,R_{D\sim Z|X}}{1-k+k\,R^{2}_{D\sim
Z|X}}\,\frac{\sigma_{Y\sim X+Z+D}}{\sigma_{D\sim X+Z}}.$

From here onwards, partialing out $X$ and $Z$ is only indicated in the
estimands. In order to rewrite and simplify the second summand, we invoke
Lemma 6 and the rule on decomposition of unexplained variance

$\displaystyle\beta_{Y\sim D|X,Z}-\beta_{Y\sim D|X,Z,U}$
$\displaystyle=\sum_{j=1}^{2}R_{Y\sim U_{j}|D,U_{-j}}\frac{\sigma_{Y\sim
D+U_{-j}}}{\sigma_{U_{j}\sim D+U_{-j}}}R_{D\sim
U_{j}}\frac{\sigma_{U_{j}}}{\sigma_{D}}$
$\displaystyle\stackrel{{\scriptstyle\text{[iii]}}}{{=}}\sum_{j=1}^{2}R_{Y\sim
U_{j}|D,U_{-j}}R_{D\sim U_{j}}\sqrt{\frac{1-R^{2}_{Y\sim
U_{-j}|D}}{1-R^{2}_{U_{j}\sim D+U_{-j}}}}\frac{\sigma_{Y\sim D}}{\sigma_{D}}.$
(33)

Due to rule [i] and the conditional independence assumption, $R_{U_{1}\sim
U_{2}}=0$ holds. This result can be used to rewrite $R_{U_{j}\sim U_{-j}|D}$
via the recursive partial correlation formula [iv]; moreover, we use the
decomposition of unexplained variance [iii] on $1-R^{2}_{U_{i}\sim D+U_{-i}}$
as follows

$\displaystyle R_{U_{j}\sim
U_{-j}|D}\stackrel{{\scriptstyle\text{[iv]}}}{{=}}\frac{R_{U_{j}\sim
U_{-j}}-R_{U_{j}\sim D}\,R_{U_{-j}\sim D}}{\sqrt{1-R^{2}_{U_{j}\sim
D}}\sqrt{1-R^{2}_{U_{-j}\sim D}}}=-f_{D\sim U_{j}}\,f_{D\sim U_{-j}},$ (34)
$\displaystyle 1-R^{2}_{U_{i}\sim
D+U_{-i}}\stackrel{{\scriptstyle\text{[iii]}}}{{=}}(1-R^{2}_{D\sim
U_{i}})(1-R^{2}_{U_{i}\sim U_{-i}|D}).$

Inserting these relationships into (33), we find

$\displaystyle\beta_{Y\sim D|X,Z}-\beta_{Y\sim D|X,Z,U}$
$\displaystyle=\sum_{j=1}^{2}R_{Y\sim U_{j}|D,U_{-j}}\,f_{D\sim
U_{j}}\sqrt{\frac{1-R^{2}_{Y\sim U_{-j}|D}}{1-f^{2}_{D\sim U_{1}}f^{2}_{D\sim
U_{2}}}}\frac{\sigma_{Y\sim D}}{\sigma_{D}}$ (35)
$\displaystyle=\sum_{j=1}^{2}R_{j}\,f_{j}\sqrt{\frac{1-R^{2}_{Y\sim
U_{-j}|D}}{1-f^{2}_{1}f^{2}_{2}}}\frac{\sigma_{Y\sim D}}{\sigma_{D}}.$

Lastly, we aim to express $\sqrt{1-R^{2}_{Y\sim U_{-j}|D}}$ in terms of the
other sensitivity parameters. To this end, we use the three-variable identity
[vi] with $Y\equiv Y$, $X\equiv U_{-j}$, $W\equiv U_{j}$ and $Z\equiv D$,
where we replace $R_{U_{j}\sim U_{-j}|D}$ according to (34).

$\frac{R_{-j}}{\sqrt{1-R^{2}_{-j}}}\sqrt{1-R^{2}_{j}}-f_{1}f_{2}\,R_{j}\stackrel{{\scriptstyle\text{[vi]}}}{{=}}f_{Y\sim
U_{-j}|D}\sqrt{1-f^{2}_{1}f^{2}_{2}}.$

By definition, the identity $\sqrt{1-R^{2}}=1/\sqrt{1+f^{2}}$ holds true for
any (partial) $R^{2}$ and its corresponding $f^{2}$. Thus, we get

$\sqrt{1-R^{2}_{Y\sim
U_{-j}|D}}=\left[1+\frac{\left(\frac{R_{-j}}{\sqrt{1-R^{2}_{-j}}}\sqrt{1-R^{2}_{j}}-f_{1}f_{2}\,R_{j}\right)^{2}}{1-f_{1}^{2}f_{2}^{2}}\right]^{-1/2}.$

Substituting the $\sqrt{1-R^{2}_{Y\sim U_{-j}|D}}$ term in (35) for the
expression above proves the form of the second summand that was required. ∎

##  Appendix C Derivation of Constraints in Section 4

###  C.1 Ordinary Least Squares

Specifying a comparative bound on $U\rightarrow Y$ that partials out $D$
involves two additional sensitivity parameters, $R_{Y\sim U|X,Z}$ and
$R_{Y\sim U|\tilde{X},\dot{X}_{I},Z,D}$. The former is related to $R_{Y\sim
U|X,Z,D}$ via equation (15); hence, it remains to find a relationship that
connects $R_{Y\sim U|\tilde{X},\dot{X}_{I},Z,D}$ to the other sensitivity
parameters. To this end, we employ rule [v] – reduction of partial correlation
– and the recursive partial correlation formula [iv] for $R_{Y\sim
U|\tilde{X},\dot{X}_{I},Z,D}$ and infer

$\displaystyle R_{Y\sim U|X,Z}$
$\displaystyle\stackrel{{\scriptstyle\text{[v]}}}{{=}}\frac{R_{Y\sim
U|\tilde{X},\dot{X}_{I},Z}}{\sqrt{1-R^{2}_{Y\sim\dot{X}_{I^{c}}|\tilde{X},\dot{X}_{I},Z}}}$
$\displaystyle\stackrel{{\scriptstyle\text{[iv]}}}{{=}}\begin{multlined}\frac{1}{\sqrt{1-R^{2}_{Y\sim\dot{X}_{I^{c}}|\tilde{X},\dot{X}_{I},Z}}}\Bigg{[}R_{Y\sim
D|\tilde{X},\dot{X}_{I},Z}R_{D\sim U|\tilde{X},\dot{X}_{I},Z}\\\ +R_{Y\sim
U|\tilde{X},\dot{X}_{I},Z,D}\sqrt{1-R^{2}_{Y\sim
D|\tilde{X},\dot{X}_{I},Z}}\sqrt{1-R^{2}_{D\sim
U|\tilde{X},\dot{X}_{I},Z}}\Bigg{]}.\end{multlined}\frac{1}{\sqrt{1-R^{2}_{Y\sim\dot{X}_{I^{c}}|\tilde{X},\dot{X}_{I},Z}}}\Bigg{[}R_{Y\sim
D|\tilde{X},\dot{X}_{I},Z}R_{D\sim U|\tilde{X},\dot{X}_{I},Z}\\\ +R_{Y\sim
U|\tilde{X},\dot{X}_{I},Z,D}\sqrt{1-R^{2}_{Y\sim
D|\tilde{X},\dot{X}_{I},Z}}\sqrt{1-R^{2}_{D\sim
U|\tilde{X},\dot{X}_{I},Z}}\Bigg{]}.$

This equation contains the unknown quantity $R_{D\sim
U|\tilde{X},\dot{X}_{I},Z}$ which can be expressed in terms of $R_{D\sim
U|X,Z}$ via rule [v]

$R_{D\sim U|X,Z}\stackrel{{\scriptstyle\text{[v]}}}{{=}}\frac{R_{D\sim
U|\tilde{X},\dot{X}_{I},Z}}{\sqrt{1-R^{2}_{D\sim\dot{X}_{I^{c}}|\tilde{X},\dot{X}_{I},Z}}}.$

Plugging this relationship into the equation above, we arrive at the
constraint

$R_{Y\sim
U|X,Z}=\frac{1}{\sqrt{1-R^{2}_{Y\sim\dot{X}_{I^{c}}|\tilde{X},\dot{X}_{I},Z}}}\Bigg{[}R_{Y\sim
D|\tilde{X},\dot{X}_{I},Z}\sqrt{1-R^{2}_{D\sim\dot{X}_{I^{c}}|\tilde{X},\dot{X}_{I},Z}}R_{D\sim
U|X,Z}\\\ +R_{Y\sim U|\tilde{X},\dot{X}_{I},Z,D}\sqrt{1-R^{2}_{Y\sim
D|\tilde{X},\dot{X}_{I},Z}}\sqrt{1-R^{2}_{D\sim
U|X,Z}\big{(}1-R^{2}_{D\sim\dot{X}_{I^{c}}|\tilde{X},\dot{X}_{I},Z}\big{)}}\Bigg{]}.$

###  C.2 Two-stage Least Squares

The comparative bound on $U\leftrightarrow Z$ is in fact equivalent to a bound
on the sensitivity parameter $R_{Z\sim U|X}$. First, we relate $R_{Z\sim
U|\tilde{X},\dot{X}_{-j}}$ to $R_{Z\sim U|X}$ via the conditional independence
assumption. Recursion of partial correlation [iv] yields

$\displaystyle
0=R_{U\sim\dot{X}_{j}|\tilde{X}_{I},\dot{X}_{-j},Z}\stackrel{{\scriptstyle\text{[iv]}}}{{=}}\frac{R_{U\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j}}-R_{Z\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j}}R_{Z\sim
U|\tilde{X},\dot{X}_{-j}}}{\sqrt{1-R_{Z\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j}}^{2}}\sqrt{1-R_{Z\sim
U|\tilde{X},\dot{X}_{-j}}^{2}}}$ $\displaystyle\Leftrightarrow\quad
R_{U\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j}}=R_{Z\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j}}R_{Z\sim
U|\tilde{X},\dot{X}_{-j}}.$

Employing this relationship and rule [iv] again, we find

$\displaystyle R_{Z\sim U|X}$
$\displaystyle\stackrel{{\scriptstyle\text{[iv]}}}{{=}}\frac{R_{Z\sim
U|\tilde{X},\dot{X}_{-j}}-R_{Z\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j}}R_{U\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j}}}{\sqrt{1-R_{Z\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j}}^{2}}\sqrt{1-R_{U\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j}}^{2}}}$
$\displaystyle=R_{Z\sim
U|\tilde{X},\dot{X}_{-j}}\sqrt{\frac{1-R^{2}_{Z\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j}}}{1-R^{2}_{Z\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j}}R_{Z\sim
U|\tilde{X},\dot{X}_{-j}}^{2}}}.$

As the right-hand side above is monotone in $R_{Z\sim
U|\tilde{X},\dot{X}_{j}}$, we conclude

$R^{2}_{Z\sim U|X}\leq
b_{U\\!Z}\,R^{2}_{Z\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j}}\frac{1-R^{2}_{Z\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j}}}{1-b_{U\\!Z}\,R^{4}_{Z\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j}}}.$

If a practitioner specifies a comparative bound on $Z\rightarrow Y$, we need
to connect $R_{Y\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},U,D}$ to $R_{D\sim
U|X,Z}$ and $R_{Y\sim U|X,Z,D}$. To this end, we employ rule [vi] – the three-
variable identity – with $Y\equiv Y$, $X\equiv\dot{X}_{j}$, $W\equiv U$ and
$Z\equiv(\tilde{X},\tilde{X}_{-j},D)$ which yields

$f_{Y\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z,U,D}\sqrt{1-R^{2}_{Y\sim
U|X,Z,D}}\\\
\stackrel{{\scriptstyle\text{[vi]}}}{{=}}f_{Y\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z,D}\sqrt{1-R^{2}_{U\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z,D}}-R_{Y\sim
U|X,Z,D}\,R_{U\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z,D}.$

Furthermore, we use the conditional independence
$U\perp\kern-5.38193pt\perp\dot{X}_{j}\,|\,\tilde{X},\dot{X}_{-j},Z$ both to
simplify the following recursive partial correlation formula [iv] and to apply
the reduction of partial correlation formula [v] on $R_{D\sim U|X,Z}$

$\displaystyle R_{U\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z,D}$
$\displaystyle\stackrel{{\scriptstyle\text{[iv]}}}{{=}}\frac{R_{U\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z}-R_{D\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z}\,R_{D\sim
U|\tilde{X},\dot{X}_{-j},Z}}{\sqrt{1-R_{D\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z}^{2}}\sqrt{1-R_{D\sim
U|\tilde{X},\dot{X}_{-j},Z}^{2}}}$
$\displaystyle=-f_{D\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z}\,f_{D\sim
U|\tilde{X},\dot{X}_{-j},Z},$ $\displaystyle R_{D\sim
U|\tilde{X},\dot{X}_{-j},Z}$
$\displaystyle\stackrel{{\scriptstyle\text{[v]}}}{{=}}R_{D\sim
U|X,Z}\sqrt{1-R^{2}_{D\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z}}.$

Inserting these two relationships in the three-variable identity above and
cancelling some terms, we arrive at

$\begin{multlined}f_{Y\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z,U,D}\sqrt{1-R^{2}_{Y\sim
U|X,Z,D}}=\bigg{[}f_{Y\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z,D}\sqrt{1-R^{2}_{D\sim
U|X,Z}}\\\ +R_{Y\sim
U|X,Z,D}\,R_{D\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z}\,R_{D\sim
U|X,Z}\bigg{]}\Big{/}\sqrt{1-R^{2}_{D\sim
U|X,Z}(1-R^{2}_{D\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z})}.\end{multlined}f_{Y\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z,U,D}\sqrt{1-R^{2}_{Y\sim
U|X,Z,D}}=\bigg{[}f_{Y\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z,D}\sqrt{1-R^{2}_{D\sim
U|X,Z}}\\\ +R_{Y\sim
U|X,Z,D}\,R_{D\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z}\,R_{D\sim
U|X,Z}\bigg{]}\Big{/}\sqrt{1-R^{2}_{D\sim
U|X,Z}(1-R^{2}_{D\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z})}.$

##  Appendix D Solving the Optimization Problem

Since users can specify any number and kind of bounds on the sensitivity
parameters, the resulting constraint set $\Psi(\hat{\theta})$ is potentially
very complex. It may be non-convex and can contain multiple non-linear
equality- and inequality constraints. This only leaves few standard
optimization algorithms to compute a global solution for (2). These, however,
often require careful choice of hyper-parameters and sometimes fail to solve
the problem. For this reason, we propose an adapted grid search algorithm that
is more robust and tailored to our specific optimization problem by exploiting
the structure of $\beta$. First, we characterize the set of potential
minimizers and maximizers; then, we explain how we can use monotonicity of
equality constraints to reduce the number of dimensions of the grid search
algorithm; finally, we give the pseudocode of the algorithm and discuss its
computational complexity.

###  D.1 Characterization of the Solution

According to Theorem 1, the objective $\beta$ is identified in terms of the
sensitivity parameters $(\psi_{1},\psi_{2})=(R_{D\sim U|X,Z},R_{Y\sim
U|X,Z,D})$. Due to its monotonicity in $\psi_{2}$, the objective $\beta$
attains its optimal values on a subset of the boundary of
$\Psi(\hat{\theta})$. In order to show this, we characterize the feasible set
as

$\Psi(\hat{\theta})=\bigcup_{\psi_{1}\colon
P_{\psi_{1}}\neq\emptyset}P_{\psi_{1}},\text{ where}\quad
P_{\psi_{1}}=\\{\psi_{2}\colon(\psi_{1},\psi_{2})\in\Psi(\hat{\theta})\\}.$

For every fixed $\psi_{1}$ such that $P_{\psi_{1}}\neq\emptyset$, the
objective $\beta$ is a linear function in $\psi_{2}$. This implies that, for
any $\psi_{2}\in P_{\psi_{1}}$, we obtain

$\beta(\hat{\theta},\psi_{1},\min
P_{\psi_{1}})\,\lesseqgtr\,\beta(\hat{\theta},\psi_{1},\psi_{2})\,\lesseqgtr\,\beta(\hat{\theta},\psi_{1},\max
P_{\psi_{1}}),$

where the direction of the inequalities depends on the sign of $\psi_{1}$.
Therefore, $\psi$-values that minimize/maximize $\beta$ are contained in

$\Psi^{*}(\hat{\theta}):=\bigcup_{\psi_{1}\colon
P_{\psi_{1}}\neq\emptyset}\\{\min P_{\psi_{1}},\max P_{\psi_{1}}\\},$ (36)

which is a subset of the boundary of $\Psi(\hat{\theta})$.

Therefore, it suffices to discretize the set $\Psi^{*}(\hat{\theta})$ instead
of $\Psi(\hat{\theta})$ to find an approximate solution to the optimization
problem.

###  D.2 Transfering Bounds via Monotonicity

Regular grid search algorithms are highly computationally expensive as their
complexity grows exponentially in the number of unknown parameters. Yet, the
high computational costs can be significantly reduced by leveraging the
monotonicity of many equality constraints. We illustrate this with an example.

######  Example 1.

Suppose a practitioner specifies the following direct constraint on
${U\rightarrow D}$ and comparative constraint on $U\rightarrow Y$:

$\begin{gathered}R_{D\sim U|X,Z}\in[-0.5,0.5],\\\ R^{2}_{Y\sim
U|\tilde{X},\dot{X}_{-j},Z}\leq
2R^{2}_{Y\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z}\,\,\Leftrightarrow\,\,R^{2}_{Y\sim
U|X,Z}\leq
2\,f^{2}_{Y\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z},\end{gathered}$ (37)

where the latter equivalence is due to (14). In addition, the unknown
parameters $R_{D\sim U|X,Z}$, $R_{Y\sim U|X,Z,D}$ and $R_{Y\sim U|X,Z}$ are
constrained by the recursive partial correlation formula

$R_{Y\sim U|X,Z,D}\stackrel{{\scriptstyle\text{[iv]}}}{{=}}\frac{R_{Y\sim
U|X,Z}-R_{Y\sim D|X,Z}\,R_{D\sim U|X,Z}}{\sqrt{1-R^{2}_{Y\sim
D|X,Z}}\sqrt{1-R^{2}_{D\sim U|X,Z}}}.$ (38)

Note that, for any fixed $R_{D\sim U|X,Z}$ value, $R_{Y\sim U|X,Z,D}$ is a
linear, and hence, monotone function of $R_{Y\sim U|X,Z}$.

In this setting, brute-force grid search creates a three-dimensional grid of
points – one dimension per unknown partial $R$-value – and only keeps those
that (approximately) conform with (37) and (38). (Partial $R$\- and $f$-values
that only depend on $V$ are estimated.) The remaining points are projected
onto the $(R_{D\sim U|X,Z},R_{Y\sim U|X,Z,D})$-plane and, for every fixed
$R_{D\sim U|X,Z}$, we can find the smallest/largest value of $R_{Y\sim
U|X,Z,D}$ to approximate $\Psi^{*}(\hat{\theta})$. Hence, in this example, the
complexity of brute-force grid search is cubic in the number of points per
dimension.

Our algorithm, on the other hand, only needs to create a one-dimensional grid
of $R_{D\sim U|X,Z}$ values, i.e. discretize $[-0.5,0.5]$. For every such
value, we can compute the smallest value for $R_{Y\sim U|X,Z,D}$ by plugging
$R_{Y\sim
U|X,Z}=-\sqrt{2}\,\lvert\hat{f}_{Y\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z}\rvert$
into (38) directly and likewise for the largest value. Therefore, the
complexity only grows linearly in the number of points per dimension.

This principle of using monotonicity of the equality constraints can reduce
the dimension of the grid and applies beyond the above example. In fact, when
only bounds on $U\rightarrow D$ and $U\rightarrow Y$ are given, we solely
require a one-dimensional grid. Hence, the computational complexity of
generating equally spaced points in $\Psi^{*}(\hat{\theta})$ grows linearly in
the number of grid points. In the most general case, when any (finite) number
and kind of bound can be specified, only a three-dimensional grid is needed.
Hence, the worst complexity of the point-generation algorithm is cubic in the
number of points per dimension. Evaluating the objective over
$\Psi^{*}(\hat{\theta})$ has linear complexity in any case.

###  D.3 Adapted Grid-Search Algorithm

Input: lower and upper bounds given by $A^{l}$, $A^{u}$, $B^{l}$, $B^{u}$,
$D^{l}$, $D^{u}$, $E^{u}$, $E^{l}$, $M^{l}$, $M^{u}$, $O^{l}$, $O^{u}$,
$b_{ZY}$

Output: vectors $A$, $L$ and $U$

1 $a^{l}$ $\leftarrow\max\\{A^{l}\\};\qquad$ $a^{u}$
$\leftarrow\min\\{A^{u}\\}$

2 $m^{l}$ $\leftarrow\max\\{M^{l}\\};\qquad$ $m^{u}$
$\leftarrow\min\\{M^{u}\\}$

3 Initialize $A,L,U\in\mathbb{R}^{N_{a}}$

4 for  _$i\in[N_{a}]$_ do

5

6 $A_{i}\leftarrow a^{l}+(i-1)\,(a^{u}-a^{l})/(N_{a}-1)$

7

$d^{l}$ $\leftarrow\max\\{h_{d}(A_{i},E^{l},c_{2},c_{3},c_{4}),\,D^{l}\\}$

/* Pushing bounds onto b */

8 $d^{u}$ $\leftarrow\min\\{h_{d}(A_{i},E^{u},c_{2},c_{3},c_{4}),\,D^{u}\\}$

9 $b^{l}$ $\leftarrow\max\\{h_{b}(A_{i},d^{l}),\,B^{l}\\}$

10 $b^{u}$ $\leftarrow\min\\{h_{b}(A_{i},d^{u}),\,B^{u}\\}$

11

12 if  _$b^{l} >b^{u}$_ then

13 $A_{i}\leftarrow$ NA

14 $L_{i}\leftarrow$ NA

15 $U_{i}\leftarrow$ NA

16

17 else

$f_{g^{l}}$ $\leftarrow h_{f_{g}}(A_{i},M^{l})$

/* Pushing bounds onto g */

18 $f_{g^{u}}$ $\leftarrow h_{f_{g}}(A_{i},M^{u})$

19 $g^{l}$ $\leftarrow f_{g^{l}}/\sqrt{1+f^{2}_{g^{l}}}$

20 $g^{u}$ $\leftarrow f_{g^{u}}/\sqrt{1+f^{2}_{g^{u}}}$

21

$\mathrm{found}$ $\leftarrow$ False

/* Finding $L_{i}$ */

22 for  _$j\in[N_{b}]$ and not $\mathrm{found}$_ do

23 $B_{ij}\leftarrow b^{l}+(j-1)\,(b^{u}-b^{l})/(N_{b}-1)$

$f_{q}\leftarrow h_{f_{q}}(A_{i},B_{ij},c_{7})$

/* Computing bounds on o */

24 $o^{l}\leftarrow\max\\{-\sqrt{b_{ZY}\cdot
f_{q}^{2}/(1+f_{q}^{2})},O^{l}\\}$

25 $o^{u}\leftarrow\min\\{\sqrt{b_{ZY}\cdot f_{q}^{2}/(1+f_{q}^{2})},O^{u}\\}$

26 if  _$o^{l} <=o^{u}$_ then

27 for  _$k\in[N_{g}]$ and not $\mathrm{found}$_ do

28 $G_{ik}\leftarrow g^{l}+(j-1)\,(g^{u}-g^{l})/(N_{g}-1)$

29 $f_{o}\leftarrow h_{f_{o}}(B_{ij},G_{ik})$

30 $o\leftarrow f_{o}/\sqrt{1+f_{o}^{2}}$

31 $\mathrm{found}$ $\leftarrow$ $\mathrm{found}$ $\lor(o^{l}\leq o\land o\leq
o^{u})$

32

33 if  _$\mathrm{found}$_ then

34 $L_{i}$ $\leftarrow$ $B_{ij}$

35

36

37 if  _not $\mathrm{found}$_ then

38 $L_{i}$ $\leftarrow$ NA

39

$\mathrm{found}$ $\leftarrow$ False

/* Finding $U_{i}$ */

$\ldots$

/* Analogously to $L_{i}$ but the $B_{ij}$ decrease */

40

41

return _$A,L,U$_

Algorithm 1 Grid approximation of $\Psi^{*}(\hat{\theta})$

Our proposed algorithm first constructs a set of equally spaced points that
are (approximately) contained in $\Psi^{*}(\hat{\theta})$; then, it evaluates
$\beta$ over this set and takes the minimum/maximum of the obtained
$\beta$-values. The latter step is straightforward whereas the former is
complex when multiple interlocking constraints are present.

In order to keep the the notation short, we introduce some abbreviations:

$\displaystyle a=R_{D\sim U|X,Z},\quad b=R_{Y\sim U|X,Z,D},\quad d=R_{Y\sim
U|X,Z},\quad e=R_{Y\sim U|\tilde{X},\dot{X}_{I},Z,D},$ $\displaystyle
g=R_{Z\sim U|X,D},\quad m=R_{Z\sim U|X},\quad o=R_{Y\sim Z|X,U,D},\quad
q=R_{Y\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z,U,D},$ $\displaystyle
c_{1}=R_{Y\sim D|X,Z},\quad c_{2}=R_{Y\sim D|\tilde{X},\dot{X}_{I},Z},\quad
c_{3}=R_{D\sim\dot{X}_{I^{c}}|\tilde{X},\dot{X}_{I},Z},\quad
c_{4}=R_{D\sim\dot{X}_{I^{c}}|\tilde{X},\dot{X}_{I},Z},$ $\displaystyle
c_{5}=R_{D\sim Z|X},\quad c_{6}=R_{Y\sim Z|X,D},\quad
c_{7}=R_{Y\sim\dot{X}_{j}|\tilde{X},\dot{X}_{-j},Z,D},$ $\displaystyle
h_{b}(a,d)=\frac{d-c_{3}\,a}{\sqrt{1-c_{3}^{2}}\,\sqrt{1-a^{2}}},$
$\displaystyle
h_{d}(a,e,c_{2},c_{3},c_{4})=\frac{1}{\sqrt{1-c_{4}^{2}}}\left[e\,\sqrt{1-c_{2}^{2}}\,\sqrt{1-a^{2}(1-c_{3}^{2})}+c_{2}\,\sqrt{1-c_{3}^{2}}\,a^{2}\right],$
$\displaystyle
h_{f_{g}}(a,m)=\frac{1}{\sqrt{1-a^{2}}}\left[\sqrt{1-c_{5}^{2}}\cdot
f_{m}-c_{5}\,a\right],$ $\displaystyle
h_{f_{o}}(b,g)=\frac{1}{\sqrt{1-b^{2}}}\left[\sqrt{1-g^{2}}\cdot
f_{c_{6}}-b\,g\right],$ $\displaystyle
h_{f_{q}}(a,b,c_{7})=\frac{\sqrt{1-a^{2}}\cdot
f_{c_{7}}+c_{7}\,a\,b}{\sqrt{1-b^{2}}\sqrt{1-a^{2}(1-c_{7}^{2})}}.$

With a slight abuse of notation, the parameter $e$ and the associated
constants $c_{2}$, $c_{3}$ and $c_{4}$ as well as $q$ and the associated
constant $c_{7}$ may be scalars or vectors depending on the number of (13)-
and (22)-constraints, respectively. The notation $f_{s}$ is a shorthand for
the $f$-transformation of some scalar or vector, that is
$f_{s}=s/\sqrt{1-s^{2}}$. The functions $h_{b}$ and $h_{d}$ are abbreviations
for the right-hand sides of the equations (15) and (16), $h_{f_{o}}$ and
$h_{f_{g}}$ stem from (17) and (18) and $h_{f_{q}}$ states (23) in the new
notation. Inserting vectors instead of scalars into the functions is
interpreted as componentwise evaluation.

In order to compute a set of points that is approximately contained in
$\Psi^{*}(\hat{\theta})$, we first discretize the interval of all possible
$a$-values and construct the vector $A\in\mathbb{R}^{N_{a}}$ which contains
$N_{a}$ equally spaced points. (This corresponds to discretizing the interval
$[\min\\{\psi_{1}\colon
P_{\psi_{1}}\neq\emptyset\\},\max\\{\psi_{1}\colon{P_{\psi_{1}}\neq\emptyset}\\}]$.)
Second, we construct the vectors $L,U\in\mathbb{R}^{N_{a}}$ which approximate
the corresponding minima and maxima of $\beta$ at the respective $a$-value.
Thus, we can create the points

$\\{(A_{i},L_{i})\colon i\in[N_{a}]\\}\cup\\{(A_{i},U_{i})\colon
i\in[N_{a}]\\},$

which are (approximately) contained and equally spaced in
$\Psi^{*}(\hat{\theta})$. Evaluating the objective $\beta$ over this set has
complexity $\mathcal{O}(N_{a})$.

In case that only bounds on $U\rightarrow D$ and $U\rightarrow Y$ are
specified, the computational complexity of generating $A$, $L$ and $U$ grows
linearly in $N_{a}$ and the computed points are actually elements of
$\Psi^{*}(\hat{\theta})$ instead of merely approximating it. The two types of
bounds on $U\rightarrow D$ specify direct constraints on $a$. Hence, denoting
the vectors of upper and lower bounds on $a$ stemming from (9) and (10)
$A^{l}$ and $A^{u}$, we can construct $A$ by equally spacing $N_{a}$ points in
the interval $[\max\\{A^{l}\\},\min\\{A^{u}\\}]$. Bounds on $U\rightarrow Y$
directly constrain $b$ (11), $d$ (14) and $e$ (13). Crucially, for every fixed
$a$-value $A_{i}$, the functions $h_{d}$ and $h_{b}$ are linear in $e$ and
$d$, respectively. Hence, we can transfer bounds on $e$ onto $d$ and, thus,
update bounds on $d$; likewise, we can then push forward the bounds on $d$
onto $b$ and compute $L_{i}$ and $U_{i}$.

In case that at least one bound on $U\leftrightarrow Z$ or $Z\rightarrow Y$ is
specified, $A$ can be constructed in the same way as before whereas $L$ and
$U$ are more computationally involved. We again use the observation that many
$h$-functions are monotone in one argument in order to ”push forward” bounds.
For fixed $a$-value, we can transfer bounds on $m$ onto $g$; for fixed $a$\-
and $b$\- value, we can compute bounds on $o$; for fixed $a$\- and $g$-value,
we can compute the corresponding $o$-value. We construct $L_{i}$ and $U_{i}$
by discretizing the range of possible $b$-values (after successively pushing
bounds on $e$ and $d$ onto $b$) into $N_{b}$ points and searching for the
smallest/largest feasible value. To test whether a given $b$-value is
feasible, we construct a sequence of $N_{g}$ values of $g$ and check whether
there is at least one value such that the bounds on $o$ are satisfied.
Therefore, the computational complexity of constructing $A$, $L$ and $U$ is
$\mathcal{O}(N_{a}\cdot N_{b}\cdot N_{g})$.

Algorithm 1 contains the pseudocode of the algorithm to generate $A$, $L$ and
$U$. It concerns the case where at least one bound on $U\leftrightarrow Z$ or
$Z\rightarrow Y$ is specified. Otherwise, we could directly set
$L_{i}\leftarrow b^{l}$ and $U_{i}\leftarrow b^{u}$ in line 15.

A full implementation of the algorithm will be made available in a public
Github repository soon.

##  Appendix E Simulation Study

We investigate the empirical coverage of sensitivity intervals computed with
the bootstrap in two scenarios: a regression model with one additional
covariate and an instrumental variable model. In both set-ups, we set the
nominal level to 90 %.

###  E.1 Linear Regression Simulation

We generate a sample of $n$ i.i.d. random vectors
$(\varepsilon_{U},\varepsilon_{X},\varepsilon_{X},\varepsilon_{Y})^{T}\sim\mathcal{N}(0,\mathrm{Id})$
and compute the variables in the model using the following linear structural
equations:

$\displaystyle U$ $\displaystyle:=\varepsilon_{U},$ $\displaystyle X$
$\displaystyle:=\varepsilon_{X},$ $\displaystyle D$
$\displaystyle:=X+U+\varepsilon_{D},$ $\displaystyle Y$
$\displaystyle:=D+2X+U+\varepsilon_{Y}.$

Based on these structural equations, we derive the covariance matrix of the
involved random variables

$\operatorname*{var}\bigg{(}\begin{pmatrix}U\\\ X\\\ D\\\
Y\end{pmatrix}\bigg{)}=\begin{pmatrix}1&0&1&2\\\ 0&1&1&3\\\ 1&1&3&6\\\
2&3&6&15\end{pmatrix}.$

It can be used to compute (partial) $R$-values as well as the bias
$\beta_{Y\sim D|X}-\beta=1/2$. If the comparative constraints

$R^{2}_{D\sim U}\leq R^{2}_{D\sim X},\qquad R^{2}_{Y\sim
U}\leq\frac{4}{9}\,R^{2}_{Y\sim X}$

are specified, the partially identified region is $[1,(3+\sqrt{3})/2]$. Hence,
the true value $\beta=1$ equals the lower end of the PIR. The bounds above are
sharp in the sense that the lower end of the partially identified range can
only be reached when both inequalities are active, i.e. hold with equality.

In order to construct sensitivity intervals, we generate bootstrap samples of
the observed data and solve the corresponding optimization problems. Then, we
use percentile, basic and BCa bootstrap (Davison and Hinkley, 1997, chap. 5)
to compute the lower and upper end of the sensitivity interval. This approach
is compared to the heuristic sensitivity intervals described in Section 5 as
well as the oracle 90% confidence interval, which could be computed if $U$ was
observed.

We simulate data for different sample sizes $n$ and repeat each such
experiment 1000 times to compute the empirical coverage and length of the
sensitivity/confidence intervals. (For $n=2000$, we do not compute BCa
bootstrap intervals due to the larger number of required samples which
significantly increases computational effort.) More specifically, we evaluate
the empirical coverage of $\beta$ and the PIR for different sensitivity
intervals and adapt the notion of length. In order to account for the fact
that the length of typical confidence intervals approaches zero as
$n\to\infty$ whereas the length of valid sensitivity intervals is lower
bounded by the length of the PIR, we use the distance between the lower end of
an interval and 1, when it covers 1, as length.

Table 2:  Simulation results of the linear regression example.

$n$ Method Coverage Length $\beta$ PIR Mean Median 200 Percentile bootstrap
95.9% 93.1% 2.729 0.733 Basic bootstrap 85.7% 78.6% 0.535 0.312 BCa bootstrap
94.6% 91.4% 2.662 0.708 Heuristic 73.7% 47.6% 0.339 0.213 Oracle 90.2% - 0.127
0.124 500 Percentile bootstrap 96.2% 92.5% 0.426 0.319 Basic bootstrap 88.4%
80.8% 0.247 0.194 BCa bootstrap 95.6% 91.7% 0.415 0.309 Heuristic 71.7% 44.0%
0.160 0.126 Oracle 88.8% - 0.079 0.077 1000 Percentile bootstrap 94.3% 90.9%
0.238 0.203 Basic bootstrap 86.8% 82.0 % 0.173 0.144 BCa bootstrap 93.6% 90.5%
0.234 0.199 Heuristic 69.5% 45.7% 0.111 0.089 Oracle 90.2% - 0.054 0.053 2000
Percentile bootstrap 95.1% 91.6% 0.152 0.139 Basic bootstrap 89.7% 84.4% 0.123
0.111 BCa bootstrap - - - - Heuristic 71.3% 45.6% 0.077 0.066 Oracle 89.5% -
0.040 0.039

The results of this simulation study are summarized in Table 2. Both
percentile and BCa bootstrap exhibit coverage of the PIR close to the
envisaged level of 90%; its coverage of $\beta$ is close to 95%. The latter is
expected as the true value of $\beta$ is the lower end of the PIR. By
comparison, the empirical coverage of sensitivity intervals constructed via
basic bootstrap is 5 to 10 percentage points below the required level. Hence,
we use BCa bootstrap to construct sensitivity intervals in the data example in
the main text. Moreover, this simulation study illustrates that heuristic
sensitivity intervals do not possess frequentist coverage guarantees: the
empirical coverage of the PIR is consistently below 50%. Finally, we see that
sensitivity intervals are substantially longer than the oracle confidence
interval. We attribute the increased length to the uncertainty stemming from
estimating the constraints. In this simulation study, we did not encounter
cases where the estimated constraint set was empty on a bootstrap sample.

Figure 7:  Empirical distribution of the lower and upper end of the PIR as
well as the corresponding bootstrap distributions.

In order to investigate the coverage of bootstrap sensitivity intervals more
closely, we consider the distribution of the estimated upper and lower end of
the PIR as well as the corresponding bootstrap distributions. Figure 7 depicts
the estimates of these distributions based on 1000 repitions of the
experiment. For small sample sizes $n$, we notice that the bootstrap
distribution is both biased and skewed. Both phenomena diminish as $n$ grows
so that the bootstrap distribution approximates the target distribution more
closely. This is in line with the observation that coverage moderately
improves for larger sample sizes, especially for basic bootstrap.

###  E.2 Linear Instrumental Variable Simulation

We generate 100 i.i.d. samples from the distribution
$(\varepsilon_{U},\varepsilon_{Z},\varepsilon_{D},\varepsilon_{Y})^{T}\sim\mathcal{N}(0,\mathrm{Id})$
and compute the variables of the model as follows

$\displaystyle U$ $\displaystyle:=\varepsilon_{U},$ $\displaystyle Z$
$\displaystyle:=\varepsilon_{Z},$ $\displaystyle D$
$\displaystyle:=Z+U+\varepsilon_{D},$ $\displaystyle Y$
$\displaystyle:=D+U+\varepsilon_{Y}.$

This data-generating process fulfills the instrumental variable assumptions
which renders $\beta=1$ point identified. Hence, a sensitivity interval where
the IV-related sensitivity parameters are set to zero ought to be comparable
with the confidence interval that is based on the asymptotic normality of the
TSLS estimator. In order to use Algorithm 1, we slightly relax the IV
assumptions requiring $R_{Z\sim U|X},R_{Y\sim Z|X,U,D}\\!\in[-0.002,0.002]$
and further set $R_{D\sim U|X,Z}\in[-0.999,0.999]$ to bound it away from $-1$
and $1$.

We compute the empirical coverage and length of sensitivity intervals
constructed via percentile, basic and BCa bootstrap, the heuristic sensitivity
intervals and the oracle confidence intervals over 500 repitions of the
experiment. Due to the high computational costs, we conduct this simulation
study only for sample size $n=100$.

Table 3:  Simulation results of the instrumental variable example.

Method Coverage Length Mean Median Percentile bootstrap 91.2% 0.327 0.300
Basic bootstrap 93.8% 0.260 0.235 BCa bootstrap 91.2% 0.329 0.300 Heuristic
98.6% 1.070 0.567 Oracle 93.0% 0.276 0.252

The results of this experiment are stated in Table 3. We notice that the
bootstrap sensitivity intervals are on par with the oracle confidence
interval, both in terms of coverage and length. By contrast, the heuristic
sensitivity intervals exhibit very high coverage but their length is too long
to be informative in practice. In this simulation study, 12 of the $500\cdot
500=250,000$ constructed bootstrap samples led to an empty constraint set. In
these cases, we set the solution of the optimization problem to $-\infty$ and
$\infty$, respectively.

##  Appendix F Derivation of Comparison Points

In Section 7.2, we introduce comparison points to provide context for the
strength of the unmeasured confounder $U$ in relation to an observed covariate
$\dot{X}_{j}$. To shorten the notation in the following, we introduce the
abbreviation $W:=(\tilde{X},\dot{X}_{-j},Z)$. In the main text, we consider
the case where $U$ explains $b$ times as much variance in $D$ and $Y$ as
$\dot{X}_{j}$ does – conditional on $W$ or $(W,D)$. Here, we consider the more
general case, where the multiplicative constant $b$ may be different for the
relationship between $U\rightarrow D$ and $U\rightarrow Y$. We indicate this
with a corresponding subscript.

In the main text, we deduced the relationship

$R_{D\sim U|X,Z}=R_{Y\sim
U|W,\dot{X}_{j}}=\sqrt{b_{D}}\,f_{D\sim\dot{X}_{j}|W};$ (39)

in the following, we derive a similar formula for $R_{Y\sim U|X,Z,D}$ both for
comparisons unconditional on $D$ and conditional on $D$. In either case, we
make use of the following equation which follows from applying rule [vi] of
the $R^{2}$-calculus with $Y\equiv Y$, $X\equiv\dot{X}_{j}$, $W\equiv D$ and
$Z=W$

$R_{Y\sim
D|W,\dot{X}_{j}}=\frac{1}{R_{D\sim\dot{X}_{j}|W}}\\!\left[f_{Y\sim\dot{X}_{j}|W}\,\sqrt{1-R^{2}_{D\sim\dot{X}_{j}|W}}-f_{Y\sim\dot{X}_{j}|W,D}\,\sqrt{1-R^{2}_{Y\sim
D|W,\dot{X}_{j}}}\right]\\!.$ (40)

###  F.1 Comparison unconditional on $D$

Assuming that $U$ explains $b_{Y}$ times as much variance in $Y$ as
$\dot{X}_{j}$ does – conditional on $W$ – and that the respective correlations
have the same sign, (14) implies

$R_{Y\sim
U|W,\dot{X}_{j}}=\sqrt{b_{Y}}\,\frac{R_{Y\sim\dot{X}_{j}|W}}{\sqrt{1-R^{2}_{Y\sim\dot{X}_{j}|W}}}=\sqrt{b_{Y}}f_{Y\sim\dot{X}_{j}|W}.$

Inserting this relationship as well as (39) into (15) yields

$R_{Y\sim U|X,Z,D}=R_{Y\sim
U|W,\dot{X}_{j},D}=\frac{\sqrt{b_{Y}}f_{Y\sim\dot{X}_{j}|W}-R_{Y\sim
D|W,\dot{X}_{j}}\,\sqrt{b_{D}}f_{D\sim\dot{X}_{j}|W}}{\sqrt{1-R^{2}_{Y\sim
D|W,\dot{X}_{j}}}\sqrt{1-b_{D}f_{D\sim\dot{X}_{j}|W}^{2}}}.$

We reformulate this expression by invoking (40), multiplying numerator and
denominator with $\sqrt{1-R^{2}_{Y\sim\dot{X}_{j}|W}}$ and applying rule [iii]
of the $R^{2}$-calculus

$R_{Y\sim
U|X,Z,D}=\frac{\left(\sqrt{b_{Y}}-\sqrt{b_{D}}\right)\,R_{Y\sim\dot{X}_{j}|W}+\sqrt{b_{D}}\,\sqrt{\frac{1-R^{2}_{Y\sim
D+\dot{X}_{j}|W}}{1-R^{2}_{D\sim\dot{X}_{j}|W}}}\,f_{Y\sim\dot{X}_{j}|W,D}}{\sqrt{1-R^{2}_{Y\sim
D+\dot{X}_{j}|W}}\sqrt{1-b_{D}f_{D\sim\dot{X}_{j}|W}^{2}}}.$

The first suggested comparison point (27) follows from setting
$b_{Y}=b_{D}=b$, simplifying the resulting expression and replacing the
population values with their estimates.

###  F.2 Comparison conditional on $D$

Assuming that $U$ explains $b_{Y}$ times as much variance in $Y$ as
$\dot{X}_{j}$ does – conditional on $W$ and $D$ – and that the respective
correlations have the same sign, equation (13) yields

$R_{Y\sim U|W,D}=\sqrt{b_{Y}}R_{Y\sim\dot{X}_{j}|W,D}.$

Plugging this relationship as well as (39) into (16) leads to

$R_{Y\sim
U|W,\dot{X}_{j}}=\frac{1}{\sqrt{1-R^{2}_{Y\sim\dot{X}_{j}|W}}}\bigg{[}R_{Y\sim
D|W}\,\sqrt{b_{D}}\,R_{D\sim\dot{X}_{j}|W}\\\
+\sqrt{b_{Y}}\,R_{Y\sim\dot{X}_{j}|W,D}\,\sqrt{1-R^{2}_{Y\sim
D|W}}\sqrt{1-b_{D}R^{2}_{D\sim\dot{X}_{j}|W}}\bigg{]}.$

To get a formula for $R_{Y\sim U|X,Z,D}$, we insert the equation above into
(15), multiply numerator and denominator by
$\sqrt{1-R^{2}_{Y\sim\dot{X}_{j}|W}}$, simplify the resulting expression via
rule [iii] and obtain

$R_{Y\sim U|X,Z,D}=-\frac{R_{Y\sim
D|W,\dot{X}_{j}}\,\sqrt{b_{D}}f_{D\sim\dot{X}_{j}|W}\sqrt{1-R^{2}_{Y\sim\dot{X}_{j}|W}}}{\sqrt{1-R^{2}_{Y\sim
D+\dot{X}_{j}|W}}\sqrt{1-b_{D}f_{D\sim\dot{X}_{j}|W}^{2}}}\\\\[6.45831pt]
+\frac{R_{Y\sim
D|W}\,\sqrt{b_{D}}\,R_{D\sim\dot{X}_{j}|W}+\sqrt{b_{Y}}\,R_{Y\sim\dot{X}_{j}|W,D}\,\sqrt{1-R^{2}_{Y\sim
D|W}}\sqrt{1-b_{D}R^{2}_{D\sim\dot{X}_{j}|W}}}{\sqrt{1-R^{2}_{Y\sim
D+\dot{X}_{j}|W}}\sqrt{1-b_{D}f_{D\sim\dot{X}_{j}|W}^{2}}}.$

Invoking (40) in the first equality below and applying rule [iv] with $Y\equiv
Y$, $X\equiv\dot{X}_{j}$, $W\equiv D$ and $Z\equiv W$ in the second, we arrive
at

$\displaystyle R_{Y\sim U|X,Z,D}$
$\displaystyle=\begin{multlined}\frac{\sqrt{b_{D}}(R_{Y\sim
D|W}\,R_{D\sim\dot{X}_{j}|W}-R_{Y\sim\dot{X}_{j}|W})+\sqrt{b_{D}}\,\sqrt{\frac{1-R^{2}_{Y\sim
D+\dot{X}_{j}|W}}{1-R^{2}_{D\sim
X_{j}|W}}}\,f_{Y\sim\dot{X}_{j}|W,D}}{\sqrt{1-R^{2}_{Y\sim
D+\dot{X}_{j}|W}}\sqrt{1-b_{D}f_{D\sim\dot{X}_{j}|W}^{2}}}\\\\[4.30554pt]
+\frac{\sqrt{b_{Y}}\,R_{Y\sim\dot{X}_{j}|W,D}\,\sqrt{1-R^{2}_{Y\sim
D|W}}\sqrt{1-b_{D}R^{2}_{D\sim\dot{X}_{j}|W}}}{\sqrt{1-R^{2}_{Y\sim
D+\dot{X}_{j}|W}}\sqrt{1-b_{D}f_{D\sim\dot{X}_{j}|W}^{2}}}\end{multlined}\frac{\sqrt{b_{D}}(R_{Y\sim
D|W}\,R_{D\sim\dot{X}_{j}|W}-R_{Y\sim\dot{X}_{j}|W})+\sqrt{b_{D}}\,\sqrt{\frac{1-R^{2}_{Y\sim
D+\dot{X}_{j}|W}}{1-R^{2}_{D\sim
X_{j}|W}}}\,f_{Y\sim\dot{X}_{j}|W,D}}{\sqrt{1-R^{2}_{Y\sim
D+\dot{X}_{j}|W}}\sqrt{1-b_{D}f_{D\sim\dot{X}_{j}|W}^{2}}}\\\\[4.30554pt]
+\frac{\sqrt{b_{Y}}\,R_{Y\sim\dot{X}_{j}|W,D}\,\sqrt{1-R^{2}_{Y\sim
D|W}}\sqrt{1-b_{D}R^{2}_{D\sim\dot{X}_{j}|W}}}{\sqrt{1-R^{2}_{Y\sim
D+\dot{X}_{j}|W}}\sqrt{1-b_{D}f_{D\sim\dot{X}_{j}|W}^{2}}}$
$\displaystyle=\begin{multlined}\frac{\left(\sqrt{b_{Y}}\sqrt{1-b_{D}R^{2}_{D\sim\dot{X}_{j}|W}}-\sqrt{b_{D}}\sqrt{1-R^{2}_{D\sim\dot{X}_{j}|W}}\right)R_{Y\sim\dot{X}_{j}|W,D}\,\sqrt{1-R^{2}_{Y\sim
D|W}}}{\sqrt{1-R^{2}_{Y\sim
D+\dot{X}_{j}|W}}\sqrt{1-b_{D}f_{D\sim\dot{X}_{j}|W}^{2}}}\\\\[4.30554pt]
+\frac{\sqrt{b_{D}}\,\sqrt{\frac{1-R^{2}_{Y\sim
D+\dot{X}_{j}|W}}{1-R^{2}_{D\sim
X_{j}|W}}}\,f_{Y\sim\dot{X}_{j}|W,D}}{\sqrt{1-R^{2}_{Y\sim
D+\dot{X}_{j}|W}}\sqrt{1-b_{D}f_{D\sim\dot{X}_{j}|W}^{2}}}.\end{multlined}\frac{\left(\sqrt{b_{Y}}\sqrt{1-b_{D}R^{2}_{D\sim\dot{X}_{j}|W}}-\sqrt{b_{D}}\sqrt{1-R^{2}_{D\sim\dot{X}_{j}|W}}\right)R_{Y\sim\dot{X}_{j}|W,D}\,\sqrt{1-R^{2}_{Y\sim
D|W}}}{\sqrt{1-R^{2}_{Y\sim
D+\dot{X}_{j}|W}}\sqrt{1-b_{D}f_{D\sim\dot{X}_{j}|W}^{2}}}\\\\[4.30554pt]
+\frac{\sqrt{b_{D}}\,\sqrt{\frac{1-R^{2}_{Y\sim
D+\dot{X}_{j}|W}}{1-R^{2}_{D\sim
X_{j}|W}}}\,f_{Y\sim\dot{X}_{j}|W,D}}{\sqrt{1-R^{2}_{Y\sim
D+\dot{X}_{j}|W}}\sqrt{1-b_{D}f_{D\sim\dot{X}_{j}|W}^{2}}}.$

The second suggested comparison point (28) follows from setting
$b_{Y}=b_{D}=b$, simplifying the resulting expression and replacing the
population values with their estimates.

##  Appendix G Choice of Hyper-parameters

In Table 4, we list the hyper-parameters of Algorithm 1 that were used for
different data analyses. The mesh size of the grid is the same in every
dimension, that is the numbers of points considered per dimension $N_{a}$,
$N_{b}$, and $N_{g}$ are equal. We define
$N_{\text{grid}}:=N_{a}=N_{b}=N_{g}$. The number of points per dimension for
$b$-contour plots and the number of bootstrap samples are denoted by
$N_{\text{b-contour}}$ and $N_{\text{boot}}$, respectively.

In the simulation study and data example in this work, we found that the PIR
estimates change only marginally for values of $N_{\text{grid}}$ larger than
200. We recommend to consider at least 100 points per grid dimension, i.e.
$N_{\text{grid}}=100$. The rough structure of the $b$-contours often becomes
apparent for $N_{\text{b-contour}}$ as low as $10$. Due to the computational
costs of the optimization problem, we choose a relatively low number of
bootstrap samples. The simulation studies empirically confirm that percentile
and BCa bootstrap sensitivity intervals achieve good coverage nonetheless.

Table 4:  Hyper-parameters for different plots and simulation examples.

$N_{\text{grid}}$ $N_{\text{b-contour}}$ $N_{\text{boot}}$ Figure 2 200 - 3500
Figure 4 200 30 - Figure 4 150 30 - Figure 5 400 - - Figure 6 300 - - Table 2
200 - 1500 Table 3 100 - 500