# Valid causal inference with unobserved confounding in high-dimensional
settings

Niloofar Moosavi, Tetiana Gorbach, Xavier de Luna

Department of Statistics, USBE, Umeå University, Umeå, Sweden
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Various methods have recently been proposed to estimate causal effects with
confidence intervals that are uniformly valid over a set of data generating
processes when high-dimensional nuisance models are estimated by post-model-
selection or machine learning estimators. These methods typically require that
all the confounders are observed to ensure identification of the effects. We
contribute by showing how valid semiparametric inference can be obtained in
the presence of unobserved confounders and high-dimensional nuisance models.
We propose uncertainty intervals which allow for unobserved confounding, and
show that the resulting inference is valid when the amount of unobserved
confounding is small relative to the sample size; the latter is formalized in
terms of convergence rates. Simulation experiments illustrate the finite
sample properties of the proposed intervals and investigate an alternative
procedure that improves the empirical coverage of the intervals when the
amount of unobserved confounding is large. Finally, a case study on the effect
of smoking during pregnancy on birth weight is used to illustrate the use of
the methods introduced to perform a sensitivity analysis to unobserved
confounding.

Keywords: Average causal effect; Double robust estimator; Inverse probability
weighting; Sensitivity analysis.

##  1 Introduction

In contrast to randomized experiments, observational studies are prone to the
presence of confounding variables which are not balanced among treated and
control individuals (Fisher, 1958). In such studies, all confounders are often
assumed to be observed for identifiability of the causal parameter of interest
(Rubin, 1974, 1990). Efforts to make this assumption plausible and the use of
flexible nuisance models often yield a high-dimensional setting, where the
number of nuisance parameters to be fitted may be at least of the order of the
sample size. Nonetheless, some confounders may still be unobserved. Therefore,
well-conducted observational studies should investigate the sensitivity of the
inference to the assumption of no unmeasured confounding (Cornfield et al.,
1959).

In this paper, we study sensitivity analysis for semiparametric inference on
the average causal effect of a binary treatment in high-dimensional
observational studies. In this context, under certain conditions, augmented
inverse probability weighting (AIPW) (Robins, Rotnitzky, and Zhao, 1994) and
targeted learning (Van der Laan and Rose, 2011) estimators yield uniformly
valid inferences even when nuisance models are fitted with flexible machine
learning algorithms; see, e.g., Farrell (2015); Van der Laan and Gruber
(2010); Chernozhukov et al. (2018) and Belloni, Chernozhukov, and Hansen
(2014); Moosavi, Häggström, and de Luna (2023) for a recent review. Roughly,
uniformly valid inference means that the finite sample behavior of the
estimator can be well approximated by its asymptotic distribution, even if
preliminary analysis, such as variable selection, has been performed on the
nuisance models prior to final estimation. Farrell (2015) studied uniformly
valid inference associated with the AIPW estimator when flexible estimators of
nuisance models are weakly consistent and fulfill multiplicative rate
conditions (See Assumption 2 below). Here, we build on these results to allow
for unobserved confounding and study the uniform validity of the resulting
inference. For this, we first specify the confounding bias of the AIPW
estimator of the causal effect as a function of unobserved confounding. We
then propose uncertainty intervals for the causal effect which account for
such bias (Vansteelandt et al., 2006; Genbäck and de Luna, 2019). Our
suggested inference on the causal effect using the uncertainty intervals
ignores the finite sample bias and randomness in the estimation of the
confounding bias. It provides uniformly valid inference given assumptions on
the amount of unobserved confounding relative to the sample size; the latter
is formalized in terms of convergence rate. We provide simulations to
illustrate the relevance of the theory for finite samples.

This work contributes to the sensitivity analysis literature that originated
in Cornfield et al. (1959) and has been studied further by many others (to
cite only a few: Rosenbaum, 1987; Ding and VanderWeele, 2016; Franks, D’Amour,
and Feller, 2020; Bonvini and Kennedy, 2022; Zhang and Tchetgen Tchetgen,
2022; Zhao, Small, and Bhattacharya, 2019; Gabriel, Sjölander, and Sachs,
2023). More specifically, we contribute to research which considers the
sensitivity parameter to be the correlation between the potential outcomes and
the treatment assignment given the observed covariates induced by unmeasured
confounders (Copas and Li, 1997; Imai, Keele, and Tingley, 2010; Genbäck and
de Luna, 2019). Our contribution is to consider post-model-selection
inference, including high-dimensional situations, and inference following
machine learning fits of the nuisance models. In this respect, our work is
most related to Scharfstein et al. (2021), which also considers inference on a
causal parameter based on flexible estimation of nuisance models and allows
for unobserved confounders via sensitivity analysis. Their approach differs in
how the confounding bias is parameterized. We contrast the methods using a
case study in Section 4.

The rest of the paper is organised as follows. Section 2 describes the context
and states our result on the uniform validity of the inference for a post-
model-selection estimator of a causal parameter, for a given amount of
unobserved confounding described by a correlation parameter. In practice, this
correlation is unknown and a plausible range of correlation values may be
considered as a sensitivity analysis. In Section 3, simulations are provided
to illustrate the relevance of the theory for finite samples, both in low- and
high-dimensional settings. Section 4 presents a case study of the effect of
maternal smoking on birth weight. This illustrates the use of the herein
proposed methods implemented in the R-package hdim.ui
(https://github.com/stat4reg/hdim.ui). Section 5 concludes the paper.

##  2 Theory and method

We aim to study the causal effect of a binary treatment $T$ on an outcome of
interest $Y.$ Let $Y(t)$, $t\in\\{0,1\\}$, be the potential outcomes that
would have been observed under a corresponding treatment level, and assume the
observed outcome is $Y=TY(1)+(1-T)Y(0)$, for all units in the study. The
average causal effect of the treatment is defined as
$E\left(Y\left(1\right)-Y\left(0\right)\right)$. Without loss of generality,
we focus on the parameter $\tau=E\left(Y\left(1\right)\right)$; see the
Appendix for other parameters of interest. We denote the vector of pre-
treatment covariates and their transformations by
$X=(1,X^{(1)},\ldots,X^{(p-1)})$. For simplicity, we call this the set of
covariates. The dimension $p$ of this vector is considered to increase with
the sample size $n$.

Consider
$O_{i,n}=\left\\{\left(X_{i,n},Y_{i,n},T_{i,n}\right)\right\\}_{i=1}^{n}$ to
be an i.i.d. sample that follows a distribution $P_{n}$ which is allowed to
vary with the sample size $n$. We drop the subscript $n,$ which implies
dependence on $n,$ when it is clear from the context. Moreover, let $n_{t}$
denote the number of individuals with observed treatment $T=1$. We use the
notation $E_{n}[W_{i}]=n^{-1}\Sigma_{i=1}^{n}W_{i}$,
$E_{n}[W_{i}]^{q}=(n^{-1}\Sigma_{i=1}^{n}W_{i})^{q}$ and
$E_{n_{t}}[W_{i}]=n_{t}^{-1}\Sigma_{i=1}^{n}T_{i}W_{i}$.

We consider the augmented inverse propensity weighting estimator (AIPW ; See
Robins et al. (1994) and Scharfstein, Rotnitzky, and Robins (1999)):

$\hat{\tau}_{\textsc{AIPW}}=E_{n}\left[\dfrac{T_{i}\left(Y_{i}-\hat{m}\left(X_{i}\right)\right)}{\hat{e}\left(X_{i}\right)}+\hat{m}\left(X_{i}\right)\right],$
(1)

where $\hat{m}(X)$ and $\hat{e}(X)$ are arbitrary estimators of
$m(X)=E\left(Y\left(1\right)\mid X,T=1\right)$ and ${e}(X)=E(T\mid X)$,
respectively. Let us consider the following assumptions.

###### Assumption 1

$e(X)>0$.

###### Assumption 2

Assume $\hat{e}(X)$ is constant conditional on the observed value of
$\\{X_{i},T_{i}\\}_{i=1}^{n}$ (assumption of ”no additional randomness”,
Farrell, 2018, a revision of Farrell, 2015) and

  1. a.

$E_{n}\left[\left(\hat{e}\left(X_{i}\right)-{e}\left(X_{i}\right)\right)^{2}\right]=o_{P}(1)$
and
$E_{n}\left[\left(\hat{m}\left(X_{i}\right)-{m}\left(X_{i}\right)\right)^{2}\right]=o_{P}(1)$.

  2. b.

$E_{n}\left[\left(\hat{m}\left(X_{i}\right)-{m}\left(X_{i}\right)\right)^{2}\right]^{1/2}E_{n}\left[\left(\hat{e}\left(X_{i}\right)-{e}\left(X_{i}\right)\right)^{2}\right]^{1/2}=o_{P}(n^{-1/2})$.

  3. c.

$E_{n}\left[\left(\hat{m}\left(X_{i}\right)-{m}\left(X_{i}\right)\right)\left(1-T_{i}/{e}\left(X_{i}\right)\right)\right]=o_{P}(n^{-1/2})$.

The following theorem gives the asymptotic behaviour of the AIPW estimator and
its asymptotic target parameter. The theorem generalizes Theorem 3 in Farrell
(2018) by avoiding the assumption $E\left(Y\left(1\right)\mid
X=x,T=1\right)=E\left(Y\left(1\right)\mid X=x\right)$, i.e., allowing for
unobserved confounders.

###### Theorem 1

Let Assumptions 1, 2, and Assumption A.1 in the Appendix hold. The AIPW
estimator (1) is asymptotically linear as follows:

$\displaystyle\sqrt{n}\left(\hat{\tau}_{\textsc{AIPW}}-\tau^{-}\right)=\Sigma_{i=1}^{n}\Psi_{i}/\sqrt{n}+o_{p}(1),$

where $\tau^{-}=E\left({m}\left(X\right)\right)$ and
$\Psi_{i}=\dfrac{T_{i}\left(Y_{i}-{m}\left(X_{i}\right)\right)}{{e}(X_{i})}+{m}(X_{i})-E\left({m}\left(X_{i}\right)\right)$.

Following Genbäck and de Luna (2019), we study the consequences of the
presence of unobserved confounders of the treatment-outcome relationship on
inference about $\tau$ by postulating a sensitivity model:

###### Assumption 3

Let $Y(1)=E\left(Y\left(1\right)\mid X\right)+\xi$, $T^{*}=g(X)+\eta,$
$T=\bm{I}(T^{*}>0),$ where $\eta\sim\mathcal{N}(0,1)$ (i.e., a probit link is
used), $\eta\perp\\!\\!\\!\\!\perp X,$ and $\bm{I}$ represents an indicator
function. Moreover, let $\xi=\rho\sigma\eta+\epsilon$, where
$\rho=\text{corr}(\xi,\eta)$, $\sigma^{2}=\text{var}(\xi)<\infty$ with
$\epsilon$ satisfying $\epsilon\perp\\!\\!\\!\\!\perp(X,\eta)$ and
$E(\epsilon)=0$.

The case $\rho=0$ corresponds to the inclusion of all confounders in $X$.
Assumption 3, therefore, considers a departure from the commonly used
identification assumption $E\big{(}Y\left(1\right)\mid
X=x,T=1\big{)}=E\big{(}Y\left(1\right)\mid X=x\big{)}$ (no unobserved
confounders) in order to perform a sensitivity analysis, and is thus not a
restriction of the semiparametric inference performed under the condition of
no unobserved confounders. Instead, Assumption 3 allows us to consider a
potential amount of unobserved confounding through an interpretable parameter
$\rho$, as well as to compute the size of the confounding bias of the AIPW
estimator with the following result.

###### Theorem 2

Suppose that for some $g^{0}(X)$ and $m^{0}(X)$ we have
$E_{n}\left[|\hat{g}\left(X_{i}\right)-g^{0}\left(X_{i}\right)|\right]=o_{P}(1)$
and
$E_{n}\left[\left(\hat{m}\left(X_{i}\right)-m^{0}\left(X_{i}\right)\right)^{2}\right]=o_{P}(1)$.
Assume 1 and 3. Furthermore, assume A.1 holds for $m^{0}(X)$. If
$g^{0}(X)=g(X)$ or $m^{0}(X)=m(X)$, the (asymptotic) confounding bias of the
AIPW estimator (1), which we denote by $b=\tau^{-}-\tau$, is equal to

$b=\rho\sigma E\left(\lambda\left(g\left(X\right)\right)\right),$

where $\lambda(\cdot)=\phi(\cdot)/\Phi(\cdot),$ and $\phi$ and $\Phi$ are the
probability density function and the cumulative distribution function of the
standard normal distribution, respectively.

Note that under the conditions of Theorem 1, the AIPW estimator (1) is a
semiparametric efficient estimator of $\tau^{-}$(see Hines et al., 2022,
Example 5). However, if confounders of the treatment-outcome relationship are
not observed (not included in $X$), $\rho$ and consequently $b$ would be
nonzero and hence the AIPW estimate would not be a consistent estimate of
$\tau$. It can be shown that for given values of $\rho$ and $\sigma$, the
efficient influence function of the bias term has the form
$\rho\sigma\lambda\left(g\left(X)\right)\right)-b$. This motivates us to
estimate the confounding bias $b$ empirically by plugging in the estimated
$g(X_{i})$ values. For a more realistic scenario, the parameter $\sigma$ must
also be estimated. First, we use
$\hat{b}=\rho\hat{\sigma}E_{n}\left[\lambda\left(\hat{g}\left(X\right)\right)\right]$
as an estimator of the parameter $b$, where
$\hat{\sigma}=E_{n_{t}}\left[\left(Y\left(1\right)_{i}-\hat{m}\left(X_{i}\right)\right)^{2}\right]^{1/2}$.
The estimator of the variance, $\hat{\sigma}^{2}$, is biased due to both high
dimensionality/variable selection and unobserved confounders. Below, we give
the asymptotic properties of $\hat{\tau}_{\textsc{AIPW}}-\hat{b}$ as an
estimator of $\tau$, assuming approximate sparsity and known $\rho.$

###### Theorem 3

Let Assumptions 1–3 hold. Moreover, let Assumptions A.1 and A.3 hold and
assume $\sqrt{n}\rho^{3}=o(1)$, $\sqrt{n}\rho
E_{n}\left[\hat{g}\left(X_{i}\right)-g\left(X_{i}\right)\right]=o_{p}(1)$ and
$\sqrt{n}\rho
E_{n_{t}}\left[\left({m}(X_{i})-\hat{m}\left(X_{i}\right)\right)^{2}\right]^{1/2}=o_{P}(1)$.
We have:

  * a.

$\sqrt{n}\left(\left(\hat{\tau}_{\textsc{AIPW}}-\hat{b}\right)-\tau\right)=\Sigma_{i=1}^{n}\Psi_{i}/\sqrt{n}+o_{p}(1),$

  * b.

$V^{-1/2}\sqrt{n}\left(\left(\hat{\tau}_{\textsc{AIPW}}-\hat{b}\right)-\tau\right)\rightarrow_{d}\mathcal{N}(0,1)$,

  * c.

$\hat{V}-V=o_{P}(1)$,

where $V=E(\Psi_{i}^{2}),$ and
$\hat{V}=E_{n}\left[\dfrac{T_{i}\left(Y_{i}-\hat{m}\left(X_{i}\right)\right)^{2}}{\hat{e}\left(X_{i}\right)^{2}}\right]+E_{n}\left[\left(\hat{m}\left(X_{i}\right)-\hat{\tau}_{\textsc{AIPW}}\right)^{2}\right]$.

As a corollary, the following 95% confidence interval for $\tau$ given $\rho$
is uniformly valid.

###### Corollary 1

For each $n$, let ${\mathcal{P}}_{n}$ be the set of distributions obeying the
assumptions of Theorem 3. Then, we have:

$\displaystyle\sup_{P\in{\mathcal{P}}_{n}}\left|{\text{pr}}_{P}\left(\tau\in\left\\{\left(\hat{\tau}_{\textsc{AIPW}}-\hat{b}\right)\pm
c_{\alpha}\sqrt{\hat{V}/n}\right\\}\right)-\left(1-\alpha\right)\right|\rightarrow
0,$

where $c_{\alpha}=\Phi^{-1}(1-\alpha/2)$.

Proof. The corollary follows from Theorem 3 (see Farrell, 2015, Corollary 2).

The proof of Theorems 1–3 can be found in Appendix. The assumptions on $\rho$
imply that the correlation between the error terms $\xi$ and $\eta$ has to
converge to zero as the sample size grows, even if the nuisance models are
estimated parametrically. In other words, the amount of unobserved confounding
needs to be small relative to the sample size. In practice, $\rho$ is not
known and a sensitivity analysis is obtained by considering the 95%
uncertainty interval constructed as the union of all 95% confidence intervals
obtained by varying $\rho\in(\rho_{min},\rho_{max})$, a plausible interval for
$\rho$. If the latter interval contains the true value of $\rho$, then this
uncertainty interval covers the parameter $\tau$ with at least 95% probability
(Gorbach and de Luna, 2018, Corollary 4.1.1). When no parametric models are
assumed and/or the number of covariates is greater than the number of
observations, lasso regression (Tibshirani, 1996) can be used in the first
step to select low-dimensional sets of variables for fitting linear and probit
nuisance models. In other words, a linear (in the parameters) regression can
be fitted using variables selected by a preliminary lasso regression. A probit
regression for the treatment can be fitted using variables selected by a
preliminary probit-lasso regression. When it comes to the rate conditions in
Assumption 2, one can use a post-selection linear model fit for the outcome
under common sparsity assumptions on the true data generating process (for
post-lasso regression, see Farrell, 2015, Corollary 5). We are unaware of any
similar result for an estimator in a sparse probit model. However, we
investigate the performance of a post-lasso probit regression in the
simulation section.

As noted above, the estimator of $\sigma$ used to obtain $\hat{b}$ in Theorem
3 is biased. We, therefore, propose to use instead a corrected estimator of
$\sigma$. Given the true value of the sensitivity parameter, this corrected
estimator of variance gives us a consistent estimate of the causal parameter.

###### Theorem 4

Let
$E_{n}\left[|\hat{g}\left(X_{i}\right)-g\left(X_{i}\right)|\right]=o_{P}(1)$
and
$E_{n}\left[\left(\hat{m}\left(X_{i}\right)-m\left(X_{i}\right)\right)^{2}\right]=o_{P}(1)$.
Assume, further, that assumptions 1, 3, A.1 and A.2 hold. Let

$\hat{\sigma}^{2}_{\text{c}}=E_{n_{t}}[(Y(1)_{i}-\hat{m}(X_{i}))^{2}]/(1-\rho^{2}E_{n_{t}}[\hat{g}(X_{i})\lambda(\hat{g}(X_{i}))]-\rho^{2}E_{n_{t}}[\lambda^{2}(\hat{g}(X_{i}))]).$

Then, we have

$\displaystyle\hat{\sigma}^{2}_{\text{c}}\xrightarrow[]{p}{\sigma}^{2},$
$\displaystyle\hat{\tau}_{\text{AIPW}}-\hat{b}_{\text{c}}\xrightarrow[]{p}\tau,$

where
$\hat{b}_{\text{c}}=\rho\hat{\sigma}_{\text{c}}E_{n}[\lambda(\hat{g}(X_{i}))]$.

Furthermore, both Theorem 3 and Corollary 1 hold under weaker regularity
conditions when using $\hat{\sigma}_{\text{c}}$ instead of $\sigma$.

###### Theorem 5

Let Assumptions 1–3 hold. Moreover, let Assumptions A.1 and A.4 hold and
assume $\rho=o(1)$, $\sqrt{n}\rho
E_{n}\left[\hat{g}\left(X_{i}\right)-g\left(X_{i}\right)\right]=o_{p}(1)$ and
$\sqrt{n}\rho
E_{n_{t}}\left[\left({m}(X_{i})-\hat{m}\left(X_{i}\right)\right)^{2}\right]^{1/2}=o_{P}(1)$.
Then, results in Theorem 3 and Corollary 1 hold for the corrected estimator
$\hat{\tau}_{\text{AIPW}}-\hat{b}_{\text{c}}$.

Note that the condition on unobserved confounding $\sqrt{n}\rho^{3}=o(1)$ in
Theorem 3 is replaced with the weaker condition $\rho=o(1)$ in the above
theorem. This improvment is illustrated in the simulation study below, where
the finite sample performance of both estimators is studied.

In the following sections, we use the R package hdim.ui
(https://github.com/stat4reg/hdim.ui), where the estimators above have been
implemented and can be used to obtain uncertainty intervals and perform
sensitivity analysis in high-dimensional situations. The package builds on
code from the ui package (Genbäck and de Luna, 2019)
(https://cran.r-project.org/web/packages/ui/index.html).

##  3 Simulation study

The simulation study in this section aims to illustrate the finite sample
behavior of the proposed estimators. Here we consider a high-dimensional
setting, while the results for low-dimensional settings are reported in the
Supplementary Material. Data was generated according to the model in
Assumption 3. All the covariates are generated to be independently and
normally distributed with mean $0$ and variance $1$. The error terms are
generated using $(\eta,\xi)\sim\mathcal{N}(0,\Sigma)$ with
$\Sigma=\begin{bmatrix}1&\rho\\\ \rho&1\end{bmatrix}$. We let $p=n$ and
simulate outcome and treatment as:

$\displaystyle Y(1)$ $\displaystyle=2+X\beta-\rho\lambda(X\gamma)+\xi,\hskip
28.45274ptT^{*}=X\gamma+\eta,$ $\displaystyle\beta$
$\displaystyle=0.6(1,1/2,1/3,1/4,1/5,1,1/2,1/3,1/4,1/5,0,\ldots,0)^{\prime},$
$\displaystyle\gamma$
$\displaystyle=0.3(1,1/2,1/3,1/4,1/5,1,1,1,1,1,0,\ldots,0)^{\prime},$

where both models include covariates that are weakly associated with the
dependent variable. Such covariates are likely to be missed in a variable
selection step.

In order to investigate the performance of the inference under variable
selection, we use
$\hat{\tau}_{\textsc{AIPW}}^{\text{refit}}-\hat{b}^{\text{refit}}$ and
$\hat{\tau}_{\textsc{AIPW}}^{\text{refit}}-\hat{b}^{\text{refit}}_{\text{c}}.$
In these estimators, $\hat{m}(X)$ is found by refitting the linear regression
using variables selected by preliminary linear-lasso regression and
$\hat{g}(X)$ is found by refitting the probit model using variables selected
by preliminary probit-lasso regression. The linear-lasso regression is
implemented using the hdm package in R (R Core Team, 2019; Chernozhukov,
Hansen, and Spindler, 2016). The probit-lasso is implemented using the glmnet
package (Friedman et al., 2021), where the regularization parameter is found
by cross-validation. Then, 95% confidence intervals are constructed for a
range of values of the sensitivity parameter $\rho$ according to Corollary 1
using the influence-curve based standard error estimator. If empirical
coverages are close to nominal, then the corresponding uncertainty intervals
will be conservative by construction, as mentioned above (see Gorbach and de
Luna, 2018; Genbäck and de Luna, 2019). We are particularly interested in
investigating coverages of confidence intervals when varying $\rho$ as sample
size increases ($n=500,1000,1500$), since Theorem 3 requires
$\sqrt{n}\rho^{3}=o(1)$ for valid inference. The results of the simulation
study are based on $500$ Monte-Carlo replications.

Table 1 reports empirical coverages of the 95% confidence intervals for $\tau$
using $\hat{\tau}_{\textsc{AIPW}}^{\text{refit}}-b^{*}$,
$\hat{\tau}_{\textsc{AIPW}}^{\text{refit}}-\hat{b}^{\text{refit}}$ and
$\hat{\tau}_{\textsc{AIPW}}^{\text{refit}}-\hat{b}^{\text{refit}}_{\text{c}}$.
Here, $b^{*}$ is an approximation of the confounding bias $b$ using the true
values of $\sigma$ and $\rho$ and the Monte Carlo estimate of
$E\left(\lambda\left(g\left(X\right)\right)\right).$

Table 1: Empirical coverages of 95% confidence intervals for $\tau$ Estimator | $\hat{\tau}_{{AIPW}}^{{refit}}-b^{*}$ |  | $\hat{\tau}_{{AIPW}}^{{refit}}-\hat{b}^{{refit}}$ |  | $\hat{\tau}_{{AIPW}}^{{refit}}-\hat{b}^{{refit}}_{{c}}$
---|---|---|---|---|---
$n$ \$\rho$ | $0.8$ | $0.6$ | $0.4$ | $0.2$ |  | $0.8$ | $0.6$ | $0.4$ | $0.2$ |  | $0.8$ | $0.6$ | $0.4$ | $0.2$
$500$ | $0.92$ | $0.91$ | $0.91$ | $0.94$ |  | $0.80$ | $0.90$ | $0.92$ | $0.95$ |  | $0.82$ | $0.92$ | $0.93$ | $0.93$
$1000$ | $0.96$ | $0.92$ | $0.94$ | $0.94$ |  | $0.60$ | $0.87$ | $0.93$ | $0.94$ |  | $0.82$ | $0.91$ | $0.94$ | $0.96$
$1500$ | $0.94$ | $0.94$ | $0.91$ | $0.94$ |  | $0.39$ | $0.85$ | $0.91$ | $0.94$ |  | $0.84$ | $0.90$ | $0.94$ | $0.94$

The confidence intervals constructed using
$\hat{\tau}_{\textsc{AIPW}}^{\text{refit}}-b^{*}$ have minimum empirical
coverage of $0.91.$ However, when a plug-in estimator of the confounding bias
$\hat{b}^{\text{refit}}$ is used low empirical coverage may arise when $\rho$
is not small enough (relative to the sample size). In particular, when model
selection is performed, the coverages of the confidence intervals calculated
using $\hat{\tau}_{\textsc{AIPW}}^{\text{refit}}-\hat{b}^{\text{refit}}$ are
low if $\rho$ is large (relative to the sample size). The minimum coverage is
0.39, when $\rho=0.8$. However, as expected from Theorem 3, coverages get
closer to the nominal level at smaller $\rho$ values for any given sample
size. Because of the weaker condition on $\rho$ in Theorem 5, the coverage for
the estimator
$\hat{\tau}_{\textsc{AIPW}}^{\text{refit}}-\hat{b}^{\text{refit}}_{\text{c}}$
is, as expected, closer to nominal level compared to the estimator
$\hat{\tau}_{\textsc{AIPW}}^{\text{refit}}-\hat{b}^{\text{refit}}$, and this
improvement is more pronounced for large sample sizes.

##  4 Case study: Effect of maternal smoking on birth weight

We re-examine a study that aims to assess the effect of smoking during
pregnancy on birth weight (Almond et al., 2005; Cattaneo, 2010; Scharfstein et
al., 2021). The data comes from an open-access sample of 4996 individuals, a
sub-sample of approximately 500 000 singleton births in Pennsylvania between
1989 and 1991 (See Scharfstein et al. 2021 and resources therein for more
details on the data). Weight at birth in grams is the outcome variable. The
treated group consists of pregnant women who smoked during pregnancy, and the
control group consists of women who did not smoke.

We use the same set of covariates as in Scharfstein et al. (2021). As they
argue, sensitivity analysis is necessary to account for potential unobserved
confounders, such as genetic factors not observed. The observed covariates
that we consider are maternal data including numerical (age and the number of
prenatal visits) and categorical variables (education $-$less than high
school, high school, more than high school$-$ and birth order $-$one, two,
larger than two) and binary variables including white, hispanic, married,
foreign and alcohol use. Higher order and interaction terms of the numerical
variables of order up to three are also considered. Finally, all the second-
order interactions with categorical variables are added, which gives a total
of 80 terms.

The analysis in Almond et al. (2005) estimates that smoking has an average
effect of -203.2 grams on birth weight, determined by a regression model
without taking into account the effect of unobserved confounding. In
Scharfstein et al. (2021), a semiparametric approach yields an estimate of
-223 grams (95% confidence interval [-274, -172]). However, after accounting
for unobserved confounding in a sensitivity analysis, the latter study
suggests that the average effect is not more extreme than -200 grams; see (2)
below for details on the clinical assumptions used.

In our analysis, we use the AIPW estimator of average causal effect,
$E(Y(1))-E(Y(0))$. The results presented in this paper are directly applicable
to $E(Y(0))$ as well, and, therefore, to $E(Y(1))-E(Y(0))$ under corresponding
assumptions. In particular, there are now two sensitivity parameters through
the sensitivity model of Assumption 3, for both $Y(1)$ and $Y(0)$, denoted
respectively $\rho_{1}$ and $\rho_{0}$. We use subset = ’refit’ in function
ui.causal to exploit the built-in variable selection function (ui.causal
included in the R-package hdim.ui). The resulting AIPW estimation gives an
estimate of the average causal effect of -221 grams (95% confidence interval:
[-273, -168]).

In their sensitivity analysis, Scharfstein et al. (2021) argued that the
following inequalities make clinical sense:

$\displaystyle E(Y(1)|T=1)$ $\displaystyle<E(Y(0)|T=1)<E(Y(0)|T=0),$ (2)
$\displaystyle E(Y(1)|T=1)$ $\displaystyle<E(Y(1)|T=0)<E(Y(0)|T=0).$

For instance from the second line, under smoking, non-smokers are expected to
have higher average birth weight than smokers, $E(Y(1)|T=0)>E(Y(1)|T=1)$. We
choose to also use these prior clinical assumptions, which yield bounds for
our sensitivity parameters $\rho_{1}$ and $\rho_{0}$. That is because the
terms on the right and left-hand sides of the inequalities in (2) can be
unbiasedly estimated using sample averages since the data are fully observed,
whereas the terms in the middle are functions of the sensitivity parameter.
The Appendix provides a description of how they are estimated. According to
the range of sensitivity parameters obtained from the above restrictions (See
Figure 1), we have derived estimates and the uniformly valid confidence
intervals along with a resulting uncertainty interval for
$E\left(Y\left(0\right)\right)$ and $E\left(Y\left(1\right)\right)$ (See
Figure 2). The resulting uncertainty interval for the average causal effect of
maternal smoking on birth weight is [-333, 49]. Thus, while our AIPW
estimation and 95% confidence interval under the assumption that all
confounders are observed are similar to those of Scharfstein et al. (2021),
the sensitivity analysis is quite different, since they suggest that the
average effect is not more extreme than -200 grams while our analysis does not
discard a potentially much larger adverse effect (up to 333 grams weight loss)
of maternal smoking depending on the strength of confounding.

The difference in conclusions may be due to the fact that Scharfstein et al.
(2021) analysis excludes the unconfoundedness situation as one of possible
scenarios while ours do not. This discrepancy might be due to different
modeling assumptions and thereby different influence functions and
corresponding estimation procedures. For instance, Scharfstein et al. (2021)
sensitivity model establishes a link between the observed and unobserved
potential outcome densities which necessitates the requirement of common
supports; a condition stating that the support of each of the missing
potential outcomes must be a subset of the support of the corresponding
observed potential outcome. This assumption may affect the validity of the
results (See section 7 in Franks et al., 2020). Furthermore, our method only
uses clinical assumptions for bounding sensitivity parameters, whereas
Scharfstein et al. (2021) uses those assumptions for both selecting a valid
tilting function (defining the sensitivity model) and bounding sensitivity
parameters. The effect on their analysis of using alternative tilting
functions is unclear to us.

Figure 1: Left plot: estimated expection of birth weight for smokers under not
smoking, $E(Y(0)|T=1)$. This estimate is constrained by the sample average
estimations of $E(Y(1)|T=1)$ in red and $E(Y(0)|T=0)$ in blue, as(2) suggests.
Right plot: the estimation of $E(Y(1)|T=0)$. See Appendix for a description of
the estimators. From these results, plausible values for the sensitivity
parameters $\rho_{0}$ and $\rho_{1}$ are (-0.25, 0.05) and (-0.2, 0.06)
respectively; i.e., where crossing with the red and blue bounds occur in
respective plots.  Figure 2:  Estimates of average birth weight under smoking
and non-smoking during pregnancy are given in the right and left plots,
respectively; using output from ui.causal function (hdim.ui package). Solid
black lines for bias-corrected AIPW estimates,
$\hat{\tau}_{\textsc{AIPW}}^{\text{refit}}-\hat{b}^{\text{refit}}$,and 95%
confidence intervals (gray area). Red intervals are confidence intervals
assuming no unobserved confounding ($\rho_{t}=0,$ $t=0,1$), and blue intervals
are 95% uncertainty intervals(using $-0.25\leq\rho_{0}\leq 0.05$ and
$-0.20\leq\rho_{1}\leq 0.06$).

##  5 Discussion

Unobserved confounding cannot be discarded nor empirically investigated in
observational studies and therefore sensitivity analysis to the
unconfoundedness assumption should be common practice. Moreover, high-
dimensional settings are typical in observational studies using large data
sets and machine learning for nuisance models. We have presented here a novel
method to conduct sensitivity analysis in such situations using uniformly
valid estimators, which is essential in high-dimensional setting. In
particular, the sensitivity analysis proposed is based on the construction of
an uncertainty interval for the causal effect of interest which we show has
uniformly valid coverage. Finite sample experiments confirm the asymptotic
results.

We use a sensitivity model with a sensitivity parameter which is easy to
interpret and discuss with subject-matter scientists. As all sensitivity
models, ours describes potential departures from the unconfoundedness
assumption. If sensitivity is detected as is the case in the presented
application on the effect of smoking on birth weight, then this is important
information. If no sensitivity is detected then one might argue that this does
not preclude the analysis to be sensitive to other departures from the
unconfoundedness assumption. Lastly, note that our results show that we need
to let $\rho$ tend to zero with increasing sample size unless we can estimate
bias due to unobserved confounding, and hence the propensity score, with a
root-n convergence rate. An alternative to bias estimation is to fix the bias
itself as the sensitivity parameter, at the expense of interpretability.

##  6 Appendix

We use the notations $f(X)=E\left(Y\left(1\right)|X\right)$,
$\tilde{\sigma}^{2}=E_{n_{t}}\left[\left(Y(1)_{i}-{m}(X_{i})\right)^{2}\right]^{1/2}$and
$v\bm{=}E((Y(1)-m(X))^{2}|T=1)$. Further, we use $a\leq_{P}b$ to denote
$a=O_{P}(b)$.

######  Assumption A.1 (Part of Assumption 2 in Farrell (2018))

Let $U=Y(1)-m(X).$ $P_{n}$ obeys the following conditions, with bounds uniform
in $n$.

  * a.

$E\left[|U|^{4}\mid X\right]\leq\mathcal{U}^{4}$.

  * b.

For some $r>0$ : $E\left[\left|m\left(x_{i}\right)\right|^{2+2r}\right]$ and
$E\left[\left|u_{i}\right|^{4+r}\right]$ are bounded.

###### Assumption A.2

We have $E\left(\lambda^{4}\left(\hat{g}\left(X\right)\right)\right)<\infty,$
$E({g}^{4}(X))<\infty$, $E(g^{2}(X)\lambda^{2}(g(X)))<\infty$,
$E(\lambda^{4}({g}(X)))<\infty$, and
$E(\lambda^{2}(\hat{g}(X)\lambda^{2}({g}(X)))<\infty$.

###### Assumption A.3

Assume
$E\left(\left(Y\left(1\right)-f\left(X\right)\right)^{4}\right)<\infty$,
$E\left(\lambda^{2}\left({g}\left(X\right)\right)\right)<\infty$,
$E\left(\left(Y(1)-\hat{m}(X)\right)^{4}\right)<\infty$,
$E\left(g\left(X\right)\lambda\left(g\left(X\right)\right)|T=1\right)<\infty$
and
$E\left(\lambda^{2}\left(g\left(X\right)\right)|T=1)\right)<\infty$.

###### Assumption A.4

Assume
$E\left(\left(Y\left(1\right)-f\left(X\right)\right)^{4}\right)<\infty$,
$E\left(\left(Y(1)-\hat{m}(X)\right)^{4}\right)<~{}\infty$,
$E({g}^{2}(X)\lambda^{2}({g}(X)))=O_{P}(1)$,
$E(\lambda^{4}({g}(X)))=O_{P}(1)$, $\lambda(\hat{g}(X))=O_{P}(1)$ and
${g}(X)=O_{P}(1)$.

We frequently use the following lemma which is a direct result of Bishop,
Fienberg, and Holland (2007, Theorem 14.1-1). The lemma is used to translate
some regularity conditions in the form of order in probability to moment
conditions.

###### Lemma 1

Assume that var$(A_{i})<C$ across $i$ and $n$, then
$E_{n}[A_{i}]-E(A_{i})=O_{P}(n^{-1/2}).$

###  6.1 Proof of Theorem 1

Suppose that Assumptions 1, A.1, and 2 hold. By Farrell (2015, Theorem 3(1))
we have

$\displaystyle\sqrt{n}\left(E_{n}\left[\dfrac{T_{i}\left(Y_{i}-\hat{m}\left(X_{i}\right)\right)}{\hat{e}\left(X_{i}\right)}+\hat{m}\left(X_{i}\right)\right]-\tau^{-}\right)$
$\displaystyle=\dfrac{1}{\sqrt{n}}\sum_{i=1}^{n}\left(\dfrac{T_{i}\left(Y_{i}-{m}\left(X_{i}\right)\right)}{{e}(X_{i})}+{m}(X_{i})-\tau^{-}\right)+o_{P}(1).$

In the above representation, unlike Farrell (2015), the expectation
$m(X)=E\left(Y\left(1\right)\mid X,T=1\right)$ is not necessarily equal to
$E\left(Y\left(1\right)\mid X\right)$ and, therefore, the asymptotic
expectation of the AIPW estimator, denoted by $\tau^{-},$ can be different
from $\tau$.

###  6.2 Proof of Theorem 2

The steps in the proof follow Genbäck and de Luna (2019). However, the
parametric modelling assumptions are dropped here. By Farrell (2015, Theorem
2) we have that under the consistency of one of the nuisance models and other
regularity conditions specified
$\hat{\tau}_{\text{AIPW}}\xrightarrow[]{p}E\left(m\left(X\right)\right)$.
Also,

$\displaystyle E(\xi\mid X,T=1)$ $\displaystyle=E(\rho\sigma\eta+\epsilon\mid
X,T=1)$ $\displaystyle=\rho\sigma E(\eta\mid\eta>-g(X))$
$\displaystyle=\rho\sigma\lambda\left(g\left(X\right)\right).$ (3)

Using (6.2) and Assumption 3, $\tau^{-}=E\left(E\left(Y\left(1\right)\mid
X,T=1\right)\right)=\tau+\rho\sigma
E\left(\lambda\left(g\left(X\right)\right)\right).$

###  6.3 Proof of Theorem 3

If we show that $\sqrt{n}(\hat{b}-b)=o_{P}(1)$, the asymptotic linearity
result is a direct result of Theorem 1. We have

$\displaystyle\sqrt{n}(\hat{b}-b)$ $\displaystyle\begin{aligned}
=&\sqrt{n}\rho\hat{\sigma}E_{n}\left[\lambda\left(\hat{g}\left(X_{i}\right)\right)-\lambda\left({g}\left(X_{i}\right)\right)\right]\\\
&+\sqrt{n}\rho
E_{n}\left[\lambda\left({g}\left(X_{i}\right)\right)\right]\left(\hat{\sigma}-\sigma\right)\\\
&+\sqrt{n}\rho\sigma\left\\{E_{n}\left[\lambda\left({g}\left(X_{i}\right)\right)\right]-E\left(\lambda\left({g}\left(X_{i}\right)\right)\right)\right\\}.\end{aligned}$

For the first term, we have

$\displaystyle\sqrt{n}\rho\hat{\sigma}E_{n}\left[\lambda\left(\hat{g}\left(X_{i}\right)\right)-\lambda\left({g}\left(X_{i}\right)\right)\right]\leq_{P}\sqrt{n}\rho
E_{n}\left[\hat{g}\left(X_{i}\right)-g\left(X_{i}\right)\right]=o_{p}(1),$

where the inequality in probability is derived using Assumption A.3, Lemma 1,
Lipschitz continuity of inverse Mills ratio while the equality follows from
Assumptions on $\rho$ in Theorem 3. For the third term, using Lemma 1 and
Assumption A.3, one can show that
$\sqrt{n}\rho\sigma\big{\\{}E_{n}\left[\lambda\left({g}\left(X_{i}\right)\right)\right]-E\left(\lambda\left({g}\left(X_{i}\right)\right)\right)\big{\\}}=o_{p}(1)$.
This shows that the first and third terms in the above decomposition are
negligible. Therefore, using some moment conditions in Assumption A.3, we have

$\displaystyle\sqrt{n}(\hat{b}-b)$ $\displaystyle\leq_{P}\sqrt{n}\rho
E_{n}\left[\lambda\left({g}\left(X_{i}\right)\right)\right]\left(\hat{\sigma}-\sigma\right)$
$\displaystyle\leq_{P}\sqrt{n}\rho\left(\hat{\sigma}^{2}-\sigma^{2}\right)$
$\displaystyle\begin{aligned}
\leq_{P}&\sqrt{n}\rho(\hat{\sigma}^{2}-\tilde{\sigma}^{2})(=\bm{R_{21}})\\\
&+\sqrt{n}\rho(\tilde{\sigma}^{2}-v)(=\bm{R_{22}})\\\
&+\sqrt{n}\rho\left(v-\sigma^{2}\right)(=\bm{R_{23}}).\end{aligned}$

For $\bm{R_{21}}$, using Assumption A.3 and a condition on $\rho$ stated we
have

$\displaystyle\sqrt{n}\rho(\hat{\sigma}^{2}-\tilde{\sigma}^{2})$
$\displaystyle\leq_{P}\sqrt{n}\rho(\hat{\sigma}-\tilde{\sigma})$
$\displaystyle\begin{aligned} \leq_{P}&\sqrt{n}\rho
E_{n_{t}}\left[\left(Y(1)_{i}-{m}(X_{i})\right)^{2}\right]^{1/2}\\\
&+\sqrt{n}\rho
E_{n_{t}}\left[\left({m}(X_{i})-\hat{m}\left(X_{i}\right)\right)^{2}\right]^{1/2}\\\
&-\sqrt{n}\rho
E_{n_{t}}\left[\left(Y(1)_{i}-{m}(X_{i})\right)^{2}\right]^{1/2}=o_{P}(1),\end{aligned}$

where the second inequality in probability holds by Minkowski inequality.

For $\bm{R_{22}}$, using $\rho=o(1)$, Assumptions 1 and A.3 we have

$\displaystyle\sqrt{n}\rho(\tilde{\sigma}^{2}-v)$
$\displaystyle\begin{aligned}
=&\sqrt{n}\rho\left(\frac{n}{n_{t}}E_{n}[T_{i}(Y(1)_{i}-{m}(X_{i}))^{2}]-v\right)\\\
\leq_{P}&\sqrt{n}\rho\left(\frac{n}{n_{t}}E_{n}[T_{i}(Y(1)_{i}-{m}(X_{i}))^{2}]-P(T=1)^{-1}E_{n}[T_{i}(Y(1)_{i}-{m}(X_{i}))^{2}]\right)\\\
&+\sqrt{n}\rho\left(P(T=1)^{-1}E_{n}[T_{i}(Y(1)_{i}-{m}(X_{i}))^{2}]\right.\\\
&\left.-P(T=1)^{-1}E\left(T\left(Y\left(1\right)-{m}\left(X\right)\right)^{2}\right)\right)\\\
\leq_{P}&\sqrt{n}\rho\left(E_{n}[T_{i}]^{-1}-P(T=1)^{-1}\right)\\\
&+\sqrt{n}\rho\left(E_{n}[T_{i}(Y(1)_{i}-{m}(X_{i}))^{2}]-E\left(T\left(Y\left(1\right)-{m}\left(X\right)\right)^{2}\right)\right)=o_{P}(1).\end{aligned}$

For $\bm{R_{23}}$, we have

$\displaystyle E((Y(1)-m(X))^{2}|T=1)$ $\displaystyle\begin{aligned}
=&E((Y(1)-f(X)-\rho\sigma\lambda(g(X)))^{2}|T=1)\\\
=&E((Y(1)-f(X))^{2}+\rho^{2}\sigma^{2}\lambda^{2}(g(X))-2\rho\sigma(Y(1)-f(X))\lambda(g(X))|T=1)\\\
=&E((Y(1)-f(X))^{2}|T=1)+\rho^{2}\sigma^{2}E(\lambda^{2}(g(X))|T=1)\\\
&-2\rho\sigma E((Y(1)-f(X))\lambda(g(X))|T=1),\end{aligned}$

where

$\displaystyle
E((Y(1)-f(X))^{2}|T=1)=-\sigma^{2}\rho^{2}E(g(X)\lambda(g(X))|T=1)+\sigma^{2},$

using Equation (A.2) in Gorbach and de Luna (2018) and

$\displaystyle 2\rho\sigma E((Y(1)-f(X))\lambda(g(X))|T=1)=2\rho\sigma
E(\xi\lambda(g(X))|T=1)$ $\displaystyle=2\rho\sigma
E((\rho\sigma\eta+\epsilon)\lambda(g(X))|T=1)$
$\displaystyle=2\rho^{2}\sigma^{2}E(\eta\lambda(g(X))|T=1)+2\rho\sigma
E(\epsilon\lambda(g(X))|T=1)$
$\displaystyle=2\rho^{2}\sigma^{2}E(\lambda(g(X))E(\eta|X,T=1)|T=1)+2\rho\sigma
E(\lambda(g(X))E(\epsilon|X,T=1)|T=1)$
$\displaystyle=2\rho^{2}\sigma^{2}E(\lambda^{2}(g(X))|T=1).$

Therefore,

$\displaystyle E((Y(1)-m(X))^{2}|T=1)$
$\displaystyle=\sigma^{2}(1-\rho^{2}E(g(X)\lambda(g(X))|T=1)-\rho^{2}E(\lambda^{2}(g(X))|T=1)).$

Finally, using assumptions on $\rho$ and Assumption A.3, we have the following
for $\bm{R_{23}}$

$\displaystyle\sqrt{n}\rho\left(v-\sigma^{2}\right)$
$\displaystyle=\sqrt{n}\rho^{3}E\left(g\left(X\right)\lambda\left(g\left(X\right)\right)|T=1\right)-\sqrt{n}\rho^{3}\sigma^{2}E\left(\lambda^{2}\left(g\left(X\right)\right)|T=1)\right)=o_{p}(1).$

Theorem 3(b) is a direct result of Theorem 3(a) and the moment condition in
Assumption 2(c).

Theorem 3(c) holds based on the proof of Theorem 3.3 in Farrell (2015),
Assumptions 1–2 and Theorem 3(b).

###  6.4 Proof of Theorem 4

Step 1. First we find the limit of the term in the numerator of the variance
estimator. By the triangle inequality, we have

$\displaystyle|E_{n_{t}}^{1/2}[(Y(1)_{i}-{m}(X_{i}))^{2}]-E_{n_{t}}^{1/2}[(m(X_{i})-\hat{m}(X_{i}))^{2}]|$
$\displaystyle\leq E_{n_{t}}^{1/2}[(Y(1)_{i}-\hat{m}(X_{i}))^{2}]$
$\displaystyle\leq
E_{n_{t}}^{1/2}[(Y(1)_{i}-{m}(X_{i}))^{2}]+E_{n_{t}}^{1/2}[(m(X_{i})-\hat{m}(X_{i}))^{2}],$

where using Assumption 1 (which implies $n/n_{t}=O_{P}(1)$) and consistency of
$\hat{m}(X)$ we have

$\displaystyle E_{n_{t}}^{1/2}[(m(X_{i})-\hat{m}(X_{i}))^{2}]$
$\displaystyle\leq\frac{n}{n_{t}}E_{n}^{1/2}[T_{i}(m(X_{i})-\hat{m}(X_{i}))^{2}]\xrightarrow[]{p}0.$

To bound $E_{n_{t}}^{1/2}[(Y(1)_{i}-\hat{m}(X_{i}))^{2}]$ by the squeeze
theorem, we just need to find the limit of
$E_{n_{t}}^{1/2}[(Y(1)_{i}-{m}(X_{i}))^{2}]$. We have

$\displaystyle E_{n_{t}}^{1/2}[(Y(1)_{i}-{m}(X_{i}))^{2}]$
$\displaystyle=\frac{n}{n_{t}}E_{n}[T_{i}(Y(1)_{i}-{m}(X_{i}))^{2}]$
$\displaystyle\begin{aligned}
\xrightarrow[]{p}&(P(T=1))^{-1}E((Y(1)-m(X))^{2}|T=1)P(T=1)\\\
&=E((Y(1)-m(X))^{2}|T=1),\end{aligned}$

where the convergence is the result of Lemma 1 and Assumption A.2.

Note that based on the proof of Theorem 3 for $\bm{R_{23}}$, the limit found
above has the following relationship with the true parameter $\sigma^{2}$.

$\displaystyle E((Y(1)-m(X))^{2}|T=1)$
$\displaystyle=\sigma^{2}(1-\rho^{2}E(g(X)\lambda(g(X))|T=1)-\rho^{2}E(\lambda^{2}(g(X))|T=1)).$

Step 2. It now remains to show that

$\displaystyle
E_{n_{t}}[\hat{g}(X_{i})\lambda(\hat{g}(X_{i}))]-E_{n_{t}}[\lambda^{2}(\hat{g}(X_{i}))]$
$\displaystyle\xrightarrow[]{p}E(g(X)\lambda(g(X))|T=1)-E(\lambda^{2}(g(X))|T=1).$

We have

$\displaystyle
E_{n_{t}}[\hat{g}(X_{i})\lambda(\hat{g}(X_{i}))]-E(g(X)\lambda(g(X))|T=1)$
$\displaystyle\begin{aligned}
=&\frac{n}{n_{t}}E_{n}[T\hat{g}(X_{i})\lambda(\hat{g}(X_{i}))]-\frac{n}{n_{t}}E_{n}[T{g}(X_{i})\lambda(\hat{g}(X))]\\\
&+\frac{n}{n_{t}}E_{n}[T{g}(X_{i})\lambda(\hat{g}(X_{i}))]-\frac{n}{n_{t}}E_{n}[T{g}(X_{i})\lambda({g}(X_{i}))]\\\
&+\frac{n}{n_{t}}E_{n}[T{g}(X_{i})\lambda({g}(X_{i}))]-E(g(X)\lambda(g(X))|T=1)\end{aligned}$
$\displaystyle\begin{aligned}
\leq&\frac{n}{n_{t}}E_{n}^{1/2}[\lambda^{2}(\hat{g}(X_{i}))]E_{n}^{1/2}[T(\hat{g}(X_{i})-{g}(X_{i}))^{2}]\\\
&+\frac{n}{n_{t}}E_{n}^{1/2}[{g}^{2}(X_{i})]E_{n}^{1/2}[T(\lambda(\hat{g}(X_{i}))-\lambda({g}(X_{i})))^{2}]\\\
&+\frac{n}{n_{t}}E_{n}[T{g}(X_{i})\lambda({g}(X_{i}))]-E(g(X)\lambda(g(X))|T=1)\xrightarrow[]{p}0,\end{aligned}$

where the inequality is due to the Cauchy-Schwarz inequality and the
convergences follows from Assumption 1, consistency assumption for
$\hat{g}(X)$, Lipschitz continuity of the inverse Mills ratio $\lambda(.)$,
Lemma 1 and Assumption A.2. Similarly, we have

$\displaystyle
E_{n_{t}}[\lambda^{2}(\hat{g}(X_{i}))]-E(\lambda^{2}(g(X))|T=1)$
$\displaystyle=\frac{n}{n_{t}}\left(E_{n}[\lambda^{2}(\hat{g}(X_{i}))]-E_{n}[\lambda^{2}({g}(X_{i}))]+E_{n}[\lambda^{2}({g}(X_{i}))]\right)-E(\lambda^{2}(g(X))|T=1)$
$\displaystyle\begin{aligned}
=&\frac{n}{n_{t}}E_{n}[(\lambda(\hat{g}(X_{i}))-\lambda({g}(X_{i})))(\lambda(\hat{g}(X_{i}))+\lambda({g}(X_{i})))]\\\
&+\frac{n}{n_{t}}E_{n}[\lambda^{2}({g}(X_{i}))]-E(\lambda^{2}(g(X))|T=1)\end{aligned}$
$\displaystyle\begin{aligned}
\leq&\frac{n}{n_{t}}E_{n}^{1/2}[(\lambda(\hat{g}(X_{i}))-\lambda({g}(X_{i})))^{2}]E_{n}^{1/2}[(\lambda(\hat{g}(X_{i}))+\lambda({g}(X_{i})))^{2}]\\\
&+\frac{n}{n_{t}}E_{n}[\lambda^{2}({g}(X_{i}))]-E(\lambda^{2}(g(X))|T=1)\xrightarrow[]{p}0.\end{aligned}$

Step 3.  The consistency of the causal parameter estimator can be shown by the
consistency of the variance estimator $\hat{\sigma}^{2}_{\text{c}}$, Slutsky’s
theorem and
$E_{n}[\lambda(\hat{g}(X_{i}))]\xrightarrow[]{p}E(\lambda({g}(X_{i})))$. The
latter holds based on consistency of $\hat{g}(X)$, Lipschitz continuity of
$\lambda(.)$, Lemma 1 and Assumption A.2.

###  6.5 Proof of Theorem 5

We have

$\displaystyle\sqrt{n}(\hat{b}-b)$ $\displaystyle\leq_{P}\sqrt{n}\rho
E_{n}\left[\lambda\left({g}\left(X_{i}\right)\right)\right]\left(\hat{\sigma}_{\text{c}}-\sigma\right)\leq_{P}\sqrt{n}\rho\left(\hat{\sigma}_{\text{c}}^{2}-\sigma^{2}\right)$
$\displaystyle\begin{aligned}
\leq_{P}&\sqrt{n}\rho(\hat{\sigma}^{2}-\tilde{\sigma}^{2})/(1-\rho^{2}E_{n_{t}}[\hat{g}(X_{i})\lambda(\hat{g}(X_{i}))]-\rho^{2}E_{n_{t}}[\lambda^{2}(\hat{g}(X_{i}))])(=\bm{R_{21}})\\\
&+\sqrt{n}\rho(\tilde{\sigma}^{2}-v)/(1-\rho^{2}E_{n_{t}}[\hat{g}(X_{i})\lambda(\hat{g}(X_{i}))]-\rho^{2}E_{n_{t}}[\lambda^{2}(\hat{g}(X_{i}))])(=\bm{R_{22}})\\\
&+\sqrt{n}\rho\left(v/(1-\rho^{2}E_{n_{t}}[\hat{g}(X_{i})\lambda(\hat{g}(X_{i}))]-\rho^{2}E_{n_{t}}[\lambda^{2}(\hat{g}(X_{i}))])-\sigma^{2}\right)(=\bm{R_{23}}).\end{aligned}$

Where the first inequality can be shown from the first lines of the proof of
Theorem 3 in Section 6.3. Moreover, both $\bm{R_{21}}$ and $\bm{R_{22}}$ can
be shown to be $o_{p}(1)$ from the proof of Theorem 4 regarding the terms with
the same name and under some extra moment conditions in Assumption A.4. From
the proof concerning the term named $\bm{R_{23}}$ in Theorem 3 we have the
following. If Assumption 3 holds, we have
$v/(1-\rho^{2}E({g}(X_{i})\lambda({g}(X_{i})))-\rho^{2}E(\lambda^{2}({g}(X_{i}))))=\sigma^{2}.$
To complete the proof we have to show that

$\displaystyle\begin{aligned}
\sqrt{n}\rho\big{(}&v/(1-\rho^{2}E_{n_{t}}[\hat{g}(X_{i})\lambda(\hat{g}(X_{i}))]-\rho^{2}E_{n_{t}}[\lambda^{2}(\hat{g}(X_{i}))])\\\
&-v/(1-\rho^{2}E({g}(X_{i})\lambda({g}(X_{i})))-\rho^{2}E(\lambda^{2}({g}(X_{i})))\big{)}\end{aligned}$
$\displaystyle=o_{p}(1).$

Using Assumption A.4 and assumptions on $\rho$ stated in the theorem we have

$\displaystyle\bm{R_{23}}$ $\displaystyle\begin{aligned}
\leq_{P}&\sqrt{n}\rho^{3}\left(E_{n_{t}}[\hat{g}(X_{i})\lambda(\hat{g}(X_{i}))]+E_{n_{t}}[\lambda^{2}(\hat{g}(X_{i}))]\right)\\\
&-\sqrt{n}\rho^{3}\left(E({g}(X_{i})\lambda({g}(X_{i})))-E(\lambda^{2}({g}(X_{i}))\right)\\\
\leq_{P}&\sqrt{n}\rho^{3}E_{n_{t}}[\hat{g}(X_{i})\lambda(\hat{g}(X_{i}))-{g}(X_{i})\lambda(\hat{g}(X_{i}))]\\\
&+\sqrt{n}\rho^{3}E_{n_{t}}[{g}(X_{i})\lambda(\hat{g}(X_{i}))-{g}(X_{i})\lambda({g}(X_{i}))]\\\
&+\sqrt{n}\rho^{3}E_{n_{t}}[\lambda^{2}(\hat{g}(X_{i}))-\lambda^{2}({g}(X_{i}))]\\\
&+\sqrt{n}\rho^{3}\left(E_{n_{t}}[{g}(X_{i})\lambda({g}(X_{i}))]-E({g}(X_{i})\lambda({g}(X_{i})))\right)\\\
&+\sqrt{n}\rho^{3}\left(E_{n_{t}}[\lambda^{2}({g}(X_{i}))]-E(\lambda^{2}({g}(X_{i}))\right)=o_{P}(1),\end{aligned}$

which completes the proof.

###  6.6 Other parameters of interest

Theorem 1 in the paper only concerns the causal parameter
$E\left(Y\left(1\right)\right)$. The asymptotic linearity of the AIPW
estimators of the causal parameters $E\left(Y\left(0\right)\right)$,
$E\left(Y\left(1\right)|T=0\right)$ and $E\left(Y\left(0\right)|T=1\right)$ is
shown under similar regularity conditions Farrell (2018, Theorem 3-4). Under
asymptotic normality and using Assumption 3, the asymptotic confounding bias
of these estimators can be found similar to the one found for
$E\left(Y\left(1\right)\right)$ in Theorem 2. Here we show the proof for the
parameter $\tau_{10}=E\left(Y\left(1\right)|T=0\right)$. We have

$\displaystyle\begin{aligned}
\dfrac{1}{P(T=0)}\bigg{[}&E\bigg{(}E((1-T)m(X)|X)\bigg{)}\\\
&+E\bigg{(}E\big{(}\dfrac{T(Y-m(X)(1-e(X))}{e(X)}|X\big{)}\bigg{)}\bigg{]}-\tau_{10}\end{aligned}$
$\displaystyle\begin{aligned}
=\dfrac{1}{P(T=0)}\bigg{[}&E\big{(}(1-e(X))m(X)\big{)}\\\
&+E\big{(}\dfrac{1-e(X)}{{e}(X)}E(T\big{(}Y-m(X)\big{)}|X)\big{)}\bigg{]}-\tau_{10}\end{aligned}$
$\displaystyle\begin{aligned} =\dfrac{1}{P(T=0)}\bigg{[}&E((1-e(X))m(X))\\\
&+E\big{(}\dfrac{1-e(X)}{{e}(X)}E(Y(1)-m(X)|X,T=1)Pr(T=1|X)\big{)}\bigg{]}-\tau_{10}\end{aligned}$
$\displaystyle=\dfrac{1}{P(T=0)}\bigg{[}E\bigg{(}(1-Pr(T=1|X))E\big{(}Y(1)|X,T=1\big{)}\bigg{)}\bigg{]}-\tau_{10}$
$\displaystyle\begin{aligned}
=\dfrac{1}{P(T=0)}\bigg{[}&E\bigg{(}E\big{(}Y(1)|X,T=1\big{)}\\\
&-Pr(T=1|X)E\big{(}Y(1)|X,T=1\big{)}\bigg{)}\bigg{]}-\tau_{10}\end{aligned}$
$\displaystyle\begin{aligned}
=\dfrac{1}{P(T=0)}\bigg{[}&E\bigg{(}E\big{(}Y(1)|X\big{)}\\\
&+\rho\sigma\lambda\big{(}g(X)\big{)}-Pr(T=1|X)E\big{(}Y(1)|X,T=1\big{)}\bigg{)}\bigg{]}-\tau_{10}\end{aligned}$
$\displaystyle\begin{aligned} =\dfrac{1}{P(T=0)}\bigg{[}&\rho\sigma
E\bigg{(}\lambda\big{(}g(X)\big{)}\\\
&+Pr(T=0|X)E\big{(}Y(1)|X,T=0\big{)}\bigg{)}\bigg{]}-\tau_{10}\end{aligned}$
$\displaystyle=\dfrac{\rho\sigma
E\bigg{(}\lambda\big{(}g(X)\big{)}\bigg{)}}{P(T=0)}+\dfrac{E\bigg{(}Pr(T=0|X)E\big{(}Y(1)|X,T=0\big{)}\bigg{)}}{P(T=0)}-\tau_{10}$
$\displaystyle=\dfrac{\rho\sigma
E\bigg{(}\lambda\big{(}g(X)\big{)}\bigg{)}}{P(T=0)}+\dfrac{E\bigg{(}E\big{(}(1-T)Y(1)|X\big{)}\bigg{)}}{P(T=0)}-\tau_{10}$
$\displaystyle=\dfrac{\rho\sigma
E\bigg{(}\lambda\big{(}g(X)\big{)}\bigg{)}}{P(T=0)}+\dfrac{E\big{(}(1-T)Y(1)\big{)}}{P(T=0)}-\tau_{10}$
$\displaystyle=\dfrac{\rho\sigma
E\bigg{(}\lambda\big{(}g(X)\big{)}\bigg{)}}{P(T=0)}+\dfrac{P(T=0)E\big{(}Y(1)|T=0\big{)}}{P(T=0)}-\tau_{10}\dfrac{\rho\sigma
E\bigg{(}\lambda\big{(}g(X)\big{)}\bigg{)}}{P(T=0)}.$

For estimating the confounding bias at a given $\rho$ value we use
$\rho\hat{\sigma}_{\text{c}}E_{n}[\lambda(\hat{g}(X_{i}))]/E_{n}[T_{i}]$.

##  7 Supplementary material

Supplementary material contains additional simulation results and is available
online at
https://github.com/stat4reg/hdim.ui#readme.

## Acknowledgements

We are grateful to Minna Genbäck and Mohammad Ghasempour for their helpful
comments. Funding from the Marianne and Marcus Wallenberg Foundation and the
Swedish Research Council for Health, Working Life and Welfare is acknowledged.

## References

  * Almond et al. (2005) Almond, D., K. Y. Chay, and D. S. Lee (2005).  The costs of low birth weight.  The Quarterly Journal of Economics 120(3), 1031–1083.
  * Belloni et al. (2014) Belloni, A., V. Chernozhukov, and C. Hansen (2014).  Inference on treatment effects after selection among high-dimensional controls.  The Review of Economic Studies 81(2), 608–650.
  * Bishop et al. (2007) Bishop, Y. M., S. E. Fienberg, and P. W. Holland (2007).  Discrete multivariate analysis: theory and practice.  Springer Science & Business Media.
  * Bonvini and Kennedy (2022) Bonvini, M. and E. H. Kennedy (2022).  Sensitivity analysis via the proportion of unmeasured confounding.  Journal of the American Statistical Association 117(539), 1540–1550.
  * Cattaneo (2010) Cattaneo, M. D. (2010).  Efficient semiparametric estimation of multi-valued treatment effects under ignorability.  Journal of Econometrics 155(2), 138–154.
  * Chernozhukov et al. (2018) Chernozhukov, V., D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, W. Newey, and J. Robins (2018).  Double/debiased machine learning for treatment and structural parameters.  The Econometrics Journal 21(1), C1–C68.
  * Chernozhukov et al. (2016) Chernozhukov, V., C. Hansen, and M. Spindler (2016).  hdm: High-dimensional metrics.  arXiv preprint arXiv:1608.00354.
  * Copas and Li (1997) Copas, J. B. and H. G. Li (1997).  Inference for non-random samples.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 59(1), 55–95.
  * Cornfield et al. (1959) Cornfield, J., W. Haenszel, E. C. Hammond, A. M. Lilienfeld, M. B. Shimkin, and E. L. Wynder (1959).  Smoking and Lung Cancer: Recent Evidence and a Discussion of Some Questions.  JNCI: Journal of the National Cancer Institute 22(1), 173–203.
  * Ding and VanderWeele (2016) Ding, P. and T. J. VanderWeele (2016).  Sensitivity analysis without assumptions.  Epidemiology (Cambridge, Mass.) 27(3), 368–377.
  * Farrell (2015) Farrell, M. H. (2015).  Robust inference on average treatment effects with possibly more covariates than observations.  Journal of Econometrics 189(1), 1–23.
  * Farrell (2018) Farrell, M. H. (2018).  Robust inference on average treatment effects with possibly more covariates than observations.  arXiv:1309.4686v3.
  * Fisher (1958) Fisher, R. (1958).  Cigarettes, cancer, and statistics.  The Centennial Review of Arts & Science 2, 151–166.
  * Franks et al. (2020) Franks, A., A. D’Amour, and A. Feller (2020).  Flexible sensitivity analysis for observational studies without observable implications.  Journal of the American Statistical Association 115(532), 1730–1746.
  * Friedman et al. (2021) Friedman, J., T. Hastie, R. Tibshirani, and B. Narasimhan (2021).  Package ‘glmnet’.  CRAN R Repositary.
  * Gabriel et al. (2023) Gabriel, E. E., A. Sjölander, and M. C. Sachs (2023).  Nonparametric bounds for causal effects in imperfect randomized experiments.  Journal of the American Statistical Association 118(541), 684–692.
  * Genbäck and de Luna (2019) Genbäck, M. and X. de Luna (2019).  Causal inference accounting for unobserved confounding after outcome regression and doubly robust estimation.  Biometrics 75(2), 506–515.
  * Gorbach and de Luna (2018) Gorbach, T. and X. de Luna (2018).  Inference for partial correlation when data are missing not at random.  Statistics & Probability Letters 141, 82–89.
  * Hines et al. (2022) Hines, O., O. Dukes, K. Diaz-Ordaz, and S. Vansteelandt (2022).  Demystifying statistical learning based on efficient influence functions.  The American Statistician 76(3), 292–304.
  * Imai et al. (2010) Imai, K., L. Keele, and D. Tingley (2010).  A general approach to causal mediation analysis.  Psychological methods 15(4), 309–334.
  * Moosavi et al. (2023) Moosavi, N., J. Häggström, and X. de Luna (2023).  The Costs and Benefits of Uniformly Valid Causal Inference with High-Dimensional Nuisance Parameters.  Statistical Science 38(1), 1 – 12.
  * R Core Team (2019) R Core Team (2019).  R: A Language and Environment for Statistical Computing.  Vienna, Austria: R Foundation for Statistical Computing.
  * Robins et al. (1994) Robins, J. M., A. Rotnitzky, and L. P. Zhao (1994).  Estimation of regression coefficients when some regressors are not always observed.  Journal of the American Statistical Association 89(427), 846–866.
  * Rosenbaum (1987) Rosenbaum, P. R. (1987).  Sensitivity analysis for certain permutation inferences in matched observational studies.  Biometrika 74(1), 13–26.
  * Rubin (1974) Rubin, D. B. (1974).  Estimating causal effects of treatments in randomized and nonrandomized studies.  Journal of educational Psychology 66(5), 688–701.
  * Rubin (1990) Rubin, D. B. (1990).  Formal mode of statistical inference for causal effects.  Journal of Statistical Planning and Inference 25(3), 279–292.
  * Scharfstein et al. (1999) Scharfstein, D., A. Rotnitzky, and J. Robins (1999).  Rejoinder to comments on “adjusting for non-ignorable drop-out using semiparametric non-response models?”.  Journal of the American Statistical Association 94, 1121–1146.
  * Scharfstein et al. (2021) Scharfstein, D. O., R. Nabi, E. H. Kennedy, M.-Y. Huang, M. Bonvini, and M. Smid (2021).  Semiparametric sensitivity analysis: Unmeasured confounding in observational studies.  arXiv preprint arXiv:2104.08300.
  * Tibshirani (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Methodological) 58(1), 267–288.
  * Van der Laan and Gruber (2010) Van der Laan, M. J. and S. Gruber (2010).  Collaborative double robust targeted maximum likelihood estimation.  The International Journal of Biostatistics 6(1), Article 17.
  * Van der Laan and Rose (2011) Van der Laan, M. J. and S. Rose (2011).  Targeted learning: causal inference for observational and experimental data.  Springer Science & Business Media.
  * Vansteelandt et al. (2006) Vansteelandt, S., E. Goetghebeur, M. G. Kenward, and G. Molenberghs (2006).  Ignorance and uncertainty regions as inferential tools in a sensitivity analysis.  Statistica Sinica 16(3), 953–979.
  * Zhang and Tchetgen Tchetgen (2022) Zhang, B. and E. J. Tchetgen Tchetgen (2022).  A semi‐parametric approach to model‐based sensitivity analysis in observational studies.  Journal of the Royal Statistical Society Series A 185(S2), 668–691.
  * Zhao et al. (2019) Zhao, Q., D. S. Small, and B. B. Bhattacharya (2019).  Sensitivity analysis for inverse probability weighting estimators via the percentile bootstrap.  Journal of the Royal Statistical Society Series B 81(4), 735–761.


Avg. Wealth Index | 0.056 | 0.155 | 0.216 | 0.026 | 0.141 | 0.207
Region: Somali | -0.027 | 0.075 | 0.186 | -0.035 | 0.085 | 0.253
Avg. Precipitation | -0.038 | 0.098 | 0.184 | -0.055 | 0.119 | 0.255
% Flush Toilet | 0.024 | 0.060 | 0.179 | 0.032 | 0.069 | 0.256
Avg. HAZ | 0.014 | 0.073 | 0.153 | 0.022 | 0.082 | 0.240
Region: Oromiya | -0.015 | 0.042 | 0.152 | -0.018 | 0.047 | 0.210
Avg. WAZ | 0.023 | 0.104 | 0.136 | 0.036 | 0.117 | 0.234
Urban Cluster | 0.018 | 0.057 | 0.133 | 0.025 | 0.066 | 0.194
% Tertiary Educ. | 0.016 | 0.052 | 0.129 | 0.023 | 0.061 | 0.199
% Married Women | 0.013 | 0.040 | 0.128 | 0.023 | 0.053 | 0.231
Avg. No. Bedrooms | -0.009 | 0.031 | 0.108 | -0.020 | 0.046 | 0.222
% Owns Agric. Land | -0.013 | 0.048 | 0.105 | -0.021 | 0.059 | 0.180
% Owns Livestock | -0.014 | 0.061 | 0.101 | -0.014 | 0.064 | 0.156
Altitude | 0.011 | 0.049 | 0.092 | 0.019 | 0.065 | 0.162
St. Dev. Wealth Index | 0.007 | 0.032 | 0.089 | 0.009 | 0.035 | 0.136
Nightlight Intensity | 0.008 | 0.041 | 0.088 | 0.012 | 0.051 | 0.149
St. Dev HAZ | 0.006 | 0.027 | 0.080 | 0.008 | 0.032 | 0.126
% Finished Walls | 0.007 | 0.039 | 0.078 | 0.014 | 0.055 | 0.151
Avg. Temperature | 0.003 | 0.043 | 0.078 | 0.007 | 0.067 | 0.149
% Surface Water | -0.005 | 0.026 | 0.075 | -0.009 | 0.036 | 0.142
% Health Insured | 0.004 | 0.024 | 0.073 | 0.007 | 0.031 | 0.120
St. Dev. WAZ | 0.005 | 0.024 | 0.072 | 0.010 | 0.035 | 0.147
% Non-Flush Toilet | -0.003 | 0.018 | 0.068 | -0.006 | 0.025 | 0.124
Under-5 Mortality Rate | -0.004 | 0.021 | 0.068 | -0.006 | 0.027 | 0.124
Pop. Density | 0.005 | 0.032 | 0.068 | 0.008 | 0.040 | 0.136
St. Dev. WHZ | -0.004 | 0.022 | 0.065 | -0.009 | 0.032 | 0.138
% Female Children | -0.003 | 0.019 | 0.065 | -0.007 | 0.028 | 0.134
% Piped Water | 0.004 | 0.023 | 0.062 | 0.004 | 0.029 | 0.108
Region: SNNP | 0.001 | 0.019 | 0.059 | 0.000 | 0.028 | 0.098
Avg. WHZ | -0.005 | 0.062 | 0.057 | -0.007 | 0.072 | 0.110
% Well Water | 0.002 | 0.017 | 0.055 | 0.004 | 0.024 | 0.099
% Secondary Educ. | -0.003 | 0.024 | 0.053 | -0.010 | 0.041 | 0.143
Region: Benishangul-Gumuz | 0.001 | 0.016 | 0.052 | 0.002 | 0.023 | 0.100
% Owns Car | 0.002 | 0.020 | 0.050 | 0.005 | 0.033 | 0.125
Region: Gambela | -0.001 | 0.018 | 0.050 | -0.001 | 0.028 | 0.109
% Illiterate | -0.002 | 0.022 | 0.049 | -0.003 | 0.033 | 0.105
Time to Water Source | 0.000 | 0.015 | 0.049 | 0.001 | 0.020 | 0.090
% Rudimentary Walls | -0.001 | 0.018 | 0.049 | -0.001 | 0.026 | 0.097
Region: Amhara | 0.002 | 0.017 | 0.048 | 0.005 | 0.028 | 0.105
St. Dev. Women’s Age | 0.000 | 0.012 | 0.047 | 0.000 | 0.019 | 0.096
% Owns Phone | -0.002 | 0.018 | 0.045 | -0.006 | 0.031 | 0.119
Avg. Women’s Age | 0.000 | 0.015 | 0.045 | 0.000 | 0.022 | 0.092
% Petrol as Fuel | 0.001 | 0.014 | 0.043 | 0.002 | 0.021 | 0.098
% Non-Residents | -0.001 | 0.012 | 0.043 | -0.003 | 0.020 | 0.108
Malaria Prevalence | 0.000 | 0.012 | 0.043 | 0.001 | 0.018 | 0.087
% Primary Educ. | 0.000 | 0.012 | 0.041 | 0.000 | 0.020 | 0.092
$\alpha$ | -0.279 | 0.049 | - | -0.283 | 0.051 | -
$\sigma^{2}$ | 0.171 | 0.077 | - | 0.229 | 0.072 | -
$g$ | 305.000 | 0.000 | - | 88.300 | 56.028 | -
Model Size | 13.569 | 3.900 | - | 17.041 | 4.709 | -

Table A4: Top Five Highest Probability Models Using UIP Prior (Measles Vaccination Data). | Model #1 | Model #2 | Model #3 | Model #4 | Model #5
---|---|---|---|---|---
% Knows Modern FP | x | x |  | x |
Avg. Wealth Index | x |  |  | x | x
% Owns Radio |  |  |  | x |
% Has Electricity |  | x | x |  |
% Coal as Fuel |  | x | x |  |
Avg. Children’s Age | x | x | x | x | x
St. Dev. Children’s Age | x | x | x | x | x
Latitude | x | x | x | x | x
Longitude |  |  | x |  | x
Region: Affar | x | x | x | x | x
Region: Harari | x | x |  | x |
Region: Dire Dawa |  |  | x |  |
Posterior Model Probability | 0.002 | 0.002 | 0.002 | 0.001 | 0.001
Table A5: Top Five Highest Probability Models Using Hyper-$g/n$ Prior (Measles Vaccination Data). | Model #1 | Model #2 | Model #3 | Model #4 | Model #5
---|---|---|---|---|---
St. Dev. WAZ |  |  |  | x |
% Knows Modern FP |  | x |  |  |
% Illiterate |  |  |  | x |
Avg. Wealth Index |  | x | x |  |
% Owns Radio |  |  |  | x |
% Has Electricity | x |  |  |  | x
% Pit Toilet | x |  | x | x |
% Coal as Fuel | x |  |  |  |
% Dung/Crops as Fuel |  |  |  |  | x
% Finished Walls |  |  |  | x |
% Female HH Heads |  |  |  | x |
Avg. No. Births per Woman |  |  |  | x | x
Avg. Children’s Age | x | x | x | x | x
St. Dev. Children’s Age | x | x | x | x | x
Latitude | x | x | x | x | x
Longitude | x |  |  |  | x
Time to Urban Center |  |  |  | x |
Region: Affar | x | x | x | x | x
Region: Somali |  |  |  | x |
Region: Harari |  | x | x | x |
Region: Addis Ababa |  |  |  | x | x
Region: Dire Dawa | x |  |  |  | x
Posterior Model Probability | 0.001 | 0.000 | 0.000 | 0.000 | 0.000
Table A6: Summary Statistics Bilateral Migration Data

Variable | N | Mean | Std. Dev. | Min | Q1 | Q3 | Max | Data Source
---|---|---|---|---|---|---|---|---
Migration Flow | 1406 | 13366.951 | 67942.166 | 0 | 191.25 | 5874 | 1635815 | Abel and Cohen (2019)
GDPPC (Orig.) | 1406 | 10.208 | 0.74 | 8.325 | 9.565 | 10.699 | 11.532 | CEPII Gravity DB
GDPPC (Dest.) | 1406 | 10.208 | 0.74 | 8.325 | 9.565 | 10.699 | 11.532 | CEPII Gravity DB
% Rural Pop. (Orig.) | 1406 | 21.446 | 10.408 | 2.021 | 13.49 | 28.787 | 46.223 | World Bank
% Rural Pop. (Dest.) | 1406 | 21.446 | 10.408 | 2.021 | 13.49 | 28.787 | 46.223 | World Bank
Contiguity | 1406 | 0.053 | 0.223 | 0 | 0 | 0 | 1 | CEPII Gravity DB
Distance | 1406 | 8.154 | 1.153 | 4.007 | 7.224 | 9.166 | 9.896 | CEPII Gravity DB
Common Colonizer | 1406 | 0.004 | 0.065 | 0 | 0 | 0 | 1 | CEPII Gravity DB
Any Colonial Relation | 1406 | 0.003 | 0.053 | 0 | 0 | 0 | 1 | CEPII Gravity DB
Common Official Language | 1406 | 0.08 | 0.271 | 0 | 0 | 0 | 1 | CEPII Gravity DB
Common Popular Language | 1406 | 0.094 | 0.292 | 0 | 0 | 0 | 1 | CEPII Gravity DB
Population (Orig.) | 1406 | 9.5 | 1.493 | 5.802 | 8.599 | 10.784 | 12.679 | CEPII Gravity DB
Population (Dest.) | 1406 | 9.5 | 1.493 | 5.802 | 8.599 | 10.784 | 12.679 | CEPII Gravity DB
Both EU Members | 1406 | 0.329 | 0.47 | 0 | 0 | 1 | 1 | CEPII Gravity DB
Gini Index (Orig.) | 1406 | 34.202 | 6.527 | 24.55 | 30.033 | 36.683 | 51.083 | World Bank
Gini Index (Dest.) | 1406 | 34.202 | 6.527 | 24.55 | 30.033 | 36.683 | 51.083 | World Bank
% Employed (Orig.) | 1406 | 69.006 | 6.997 | 50.315 | 65.394 | 73.933 | 84.433 | World Bank
% Employed (Dest.) | 1406 | 69.006 | 6.997 | 50.315 | 65.394 | 73.933 | 84.433 | World Bank
% Tertiary School Enrollment (Orig.) | 1406 | 73.057 | 21.614 | 19.76 | 62.431 | 84.974 | 136.695 | World Bank
% Tertiary School Enrollment (Dest.) | 1406 | 73.057 | 21.614 | 19.76 | 62.431 | 84.974 | 136.695 | World Bank
% Secondary School Enrollment (Orig.) | 1406 | 112.066 | 17.125 | 66.897 | 101.848 | 117.715 | 158.052 | World Bank
% Secondary School Enrollment (Dest.) | 1406 | 112.066 | 17.125 | 66.897 | 101.848 | 117.715 | 158.052 | World Bank
Existing Migrant Population | 1406 | 6.949 | 3.406 | 0 | 5.19 | 9.335 | 16.27 | UNDESA PD
% Working Age (Orig.) | 1406 | 65.61 | 2.565 | 59.108 | 64.133 | 66.881 | 72.814 | World Bank
% Working Age (Dest.) | 1406 | 65.61 | 2.565 | 59.108 | 64.133 | 66.881 | 72.814 | World Bank
Island Country (Orig.) | 1406 | 0.184 | 0.388 | 0 | 0 | 0 | 1 | Various
Island Country (Dest.) | 1406 | 0.184 | 0.388 | 0 | 0 | 0 | 1 | Various
Health Care % of GDP (Orig.) | 1406 | 8.911 | 2.315 | 4.28 | 7.094 | 10.503 | 17.029 | World Bank
Health Care % of GDP (Dest.) | 1406 | 8.911 | 2.315 | 4.28 | 7.094 | 10.503 | 17.029 | World Bank
GDP Growth (Orig.) | 1406 | 1.62 | 1.693 | -0.893 | 0.642 | 2.307 | 9.323 | World Bank
GDP Growth (Dest.) | 1406 | 1.62 | 1.693 | -0.893 | 0.642 | 2.307 | 9.323 | World Bank
Stability Index (Orig.) | 1406 | 66.528 | 22.383 | 8.995 | 57.74 | 81.757 | 98.58 | World Bank
Stability Index (Dest.) | 1406 | 66.528 | 22.383 | 8.995 | 57.74 | 81.757 | 98.58 | World Bank
Service % of GDP (Orig.) | 1406 | 64.211 | 6.198 | 54.322 | 58.432 | 69.433 | 79.575 | World Bank
Service % of GDP (Dest.) | 1406 | 64.211 | 6.198 | 54.322 | 58.432 | 69.433 | 79.575 | World Bank
Area (Orig.) | 1406 | 12.114 | 1.72 | 7.858 | 10.841 | 13.017 | 16.116 | CEPII
Area (Dest.) | 1406 | 12.114 | 1.72 | 7.858 | 10.841 | 13.017 | 16.116 | CEPII
Landlocked Country (Orig.) | 1406 | 0.158 | 0.365 | 0 | 0 | 0 | 1 | CEPII
Landlocked Country (Dest.) | 1406 | 0.158 | 0.365 | 0 | 0 | 0 | 1 | CEPII
Latitude Capital (Orig.) | 1406 | 38.041 | 26.02 | -44.283 | 37.5 | 52.533 | 64.15 | CEPII
Latitude Capital (Dest.) | 1406 | 38.041 | 26.02 | -44.283 | 37.5 | 52.533 | 64.15 | CEPII
Longitude Capital (Orig.) | 1406 | 9.201 | 60.926 | -99.167 | -6.25 | 24.1 | 174.783 | CEPII
Longitude Capital (Dest.) | 1406 | 9.201 | 60.926 | -99.167 | -6.25 | 24.1 | 174.783 | CEPII
Asian Country (Orig.) | 1406 | 0.079 | 0.27 | 0 | 0 | 0 | 1 | CEPII
European Country (Orig.) | 1406 | 0.684 | 0.465 | 0 | 0 | 1 | 1 | CEPII
Pacific Country (Orig.) | 1406 | 0.053 | 0.223 | 0 | 0 | 0 | 1 | CEPII
Asian Country (Dest.) | 1406 | 0.079 | 0.27 | 0 | 0 | 0 | 1 | CEPII
European Country (Dest.) | 1406 | 0.684 | 0.465 | 0 | 0 | 1 | 1 | CEPII
Pacific Country (Dest.) | 1406 | 0.053 | 0.223 | 0 | 0 | 0 | 1 | CEPII
Agriculture % of GDP (Orig.) | 1406 | 2.483 | 1.617 | 0.22 | 1.15 | 3.607 | 6.513 | World Bank
Agriculture % of GDP (Dest.) | 1406 | 2.483 | 1.617 | 0.22 | 1.15 | 3.607 | 6.513 | World Bank
Battle Deaths (Orig.) | 1406 | 0.02 | 0.097 | 0 | 0 | 0 | 0.603 | World Bank
Battle Deaths (Dest.) | 1406 | 0.02 | 0.097 | 0 | 0 | 0 | 0.603 | World Bank
Homicide Rate (Orig.) | 1406 | 4.526 | 10.911 | 0.264 | 0.722 | 2.561 | 61.21 | World Bank
Homicide Rate (Dest.) | 1406 | 4.526 | 10.911 | 0.264 | 0.722 | 2.561 | 61.21 | World Bank

Table A7: Estimation Results (Bilateral Migration Data).

| Unit Information Prior | Hyper-$g/n$ Prior
---|---|---
Variable | Post. Mean | Post. SD | PIP | Post. Mean | Post. SD | PIP
GDPPC (Orig.) | 0.400 | 0.067 | 1.000 | 0.377 | 0.084 | 0.995
Distance | -0.538 | 0.047 | 1.000 | -0.523 | 0.048 | 1.000
Population (Orig.) | 0.586 | 0.055 | 1.000 | 0.563 | 0.062 | 1.000
% Employed (Dest.) | 0.356 | 0.058 | 1.000 | 0.357 | 0.053 | 1.000
Existing Migrant Population | 1.807 | 0.049 | 1.000 | 1.808 | 0.048 | 1.000
% Working Age (Dest.) | 0.351 | 0.042 | 1.000 | 0.350 | 0.038 | 1.000
GDP Growth (Dest.) | 0.251 | 0.041 | 1.000 | 0.251 | 0.042 | 1.000
Area (Dest.) | 0.512 | 0.088 | 1.000 | 0.501 | 0.088 | 1.000
Latitude Capital (Dest.) | -0.872 | 0.140 | 1.000 | -0.866 | 0.117 | 1.000
Longitude Capital (Orig.) | -0.725 | 0.146 | 1.000 | -0.770 | 0.144 | 1.000
Longitude Capital (Dest.) | 0.777 | 0.151 | 1.000 | 0.811 | 0.154 | 1.000
Pacific Country (Orig.) | 0.788 | 0.139 | 1.000 | 0.841 | 0.143 | 1.000
Pacific Country (Dest.) | -0.947 | 0.173 | 1.000 | -0.971 | 0.166 | 1.000
Island Country (Dest.) | -0.191 | 0.037 | 0.999 | -0.188 | 0.037 | 1.000
Battle Deaths (Dest.) | -0.179 | 0.043 | 0.999 | -0.180 | 0.041 | 0.999
% Rural Pop. (Dest.) | -0.177 | 0.071 | 0.987 | -0.181 | 0.063 | 0.994
Population (Dest.) | 0.232 | 0.082 | 0.974 | 0.225 | 0.073 | 0.984
Landlocked Country (Dest.) | -0.130 | 0.056 | 0.907 | -0.144 | 0.044 | 0.982
Service % of GDP (Dest.) | 0.174 | 0.080 | 0.895 | 0.173 | 0.068 | 0.949
Service % of GDP (Orig.) | 0.126 | 0.060 | 0.881 | 0.105 | 0.065 | 0.817
% Tertiary School Enrollment (Orig.) | 0.104 | 0.051 | 0.880 | 0.110 | 0.050 | 0.914
European Country (Orig.) | 0.279 | 0.155 | 0.868 | 0.313 | 0.141 | 0.933
European Country (Dest.) | 0.260 | 0.140 | 0.817 | 0.223 | 0.165 | 0.698
% Tertiary School Enrollment (Dest.) | -0.106 | 0.068 | 0.778 | -0.122 | 0.058 | 0.898
Battle Deaths (Orig.) | 0.083 | 0.057 | 0.777 | 0.105 | 0.056 | 0.886
Stability Index (Dest.) | -0.140 | 0.103 | 0.727 | -0.160 | 0.098 | 0.826
Asian Country (Orig.) | 0.182 | 0.140 | 0.703 | 0.210 | 0.135 | 0.807
% Employed (Orig.) | -0.062 | 0.063 | 0.576 | -0.080 | 0.063 | 0.729
% Secondary School Enrollment (Dest.) | -0.065 | 0.068 | 0.554 | -0.082 | 0.065 | 0.707
GDP Growth (Orig.) | -0.027 | 0.045 | 0.314 | -0.038 | 0.049 | 0.468
Gini Index (Orig.) | -0.038 | 0.069 | 0.276 | -0.039 | 0.066 | 0.346
Gini Index (Dest.) | 0.048 | 0.092 | 0.273 | 0.034 | 0.072 | 0.262
Asian Country (Dest.) | -0.062 | 0.127 | 0.247 | -0.110 | 0.152 | 0.418
% Secondary School Enrollment (Orig.) | -0.016 | 0.036 | 0.223 | -0.043 | 0.056 | 0.477
Agriculture % of GDP (Orig.) | -0.031 | 0.069 | 0.221 | -0.073 | 0.100 | 0.444
Homicide Rate (Dest.) | -0.024 | 0.056 | 0.200 | -0.056 | 0.075 | 0.440
Common Popular Language | 0.010 | 0.025 | 0.187 | 0.019 | 0.032 | 0.341
GDPPC (Dest.) | -0.018 | 0.055 | 0.150 | -0.025 | 0.062 | 0.232
Island Country (Orig.) | -0.005 | 0.019 | 0.114 | -0.005 | 0.020 | 0.160
Agriculture % of GDP (Dest.) | 0.003 | 0.047 | 0.107 | 0.003 | 0.057 | 0.180
Health Care % of GDP (Orig.) | 0.008 | 0.033 | 0.098 | 0.015 | 0.043 | 0.190
Common Official Language | 0.004 | 0.014 | 0.091 | 0.006 | 0.019 | 0.167
Stability Index (Orig.) | -0.003 | 0.027 | 0.084 | 0.010 | 0.046 | 0.170
Health Care % of GDP (Dest.) | -0.002 | 0.028 | 0.073 | 0.000 | 0.028 | 0.116
% Rural Pop. (Orig.) | -0.004 | 0.023 | 0.072 | -0.022 | 0.050 | 0.254
% Working Age (Orig.) | -0.001 | 0.011 | 0.064 | -0.003 | 0.015 | 0.117
Landlocked Country (Orig.) | 0.002 | 0.012 | 0.062 | 0.002 | 0.014 | 0.114
Contiguity | 0.002 | 0.009 | 0.060 | 0.003 | 0.012 | 0.112
Latitude Capital (Orig.) | -0.003 | 0.024 | 0.056 | -0.002 | 0.037 | 0.123
Area (Orig.) | -0.003 | 0.019 | 0.055 | -0.002 | 0.029 | 0.120
Homicide Rate (Orig.) | -0.002 | 0.015 | 0.052 | -0.005 | 0.023 | 0.124
Common Colonizer | 0.001 | 0.007 | 0.047 | 0.003 | 0.011 | 0.120
Both EU Members | 0.000 | 0.009 | 0.040 | 0.000 | 0.014 | 0.094
Any Colonial Relation | 0.000 | 0.004 | 0.029 | 0.001 | 0.007 | 0.088
$\alpha$ | 6.894 | 0.023 | - | 6.894 | 0.023 | -
$\sigma^{2}$ | 0.720 | 0.029 | - | 0.725 | 0.029 | -
$g$ | 1406.000 | 0.000 | - | 392.997 | 113.159 | -
Model Size | 29.518 | 2.018 | - | 32.793 | 2.653 | -

Table A8: Top Five Highest Probability Models Using UIP Prior (Bilateral
Migration Data).

| Model #1 | Model #2 | Model #3 | Model #4 | Model #5
---|---|---|---|---|---
GDPPC (Orig.) | x | x | x | x | x
% Rural Pop. (Dest.) | x | x | x | x | x
Distance | x | x | x | x | x
Population (Orig.) | x | x | x | x | x
Population (Dest.) | x | x | x | x | x
Gini Index (Orig.) |  |  |  | x |
% Employed (Orig.) |  | x |  | x |
% Employed (Dest.) | x | x | x | x | x
% Tertiary School Enrollment (Orig.) | x | x | x | x | x
% Tertiary School Enrollment (Dest.) | x | x | x | x | x
% Secondary School Enrollment (Dest.) |  |  | x |  | x
Existing Migrant Population | x | x | x | x | x
% Working Age (Dest.) | x | x | x | x | x
Island Country (Dest.) | x | x | x | x | x
GDP Growth (Dest.) | x | x | x | x | x
Stability Index (Dest.) | x | x | x | x | x
Service % of GDP (Orig.) | x | x | x | x | x
Service % of GDP (Dest.) | x | x | x | x | x
Area (Dest.) | x | x | x | x | x
Landlocked Country (Dest.) | x | x | x | x | x
Latitude Capital (Dest.) | x | x | x | x | x
Longitude Capital (Orig.) | x | x | x | x | x
Longitude Capital (Dest.) | x | x | x | x | x
Asian Country (Orig.) | x | x | x |  | x
European Country (Orig.) | x | x | x |  | x
Pacific Country (Orig.) | x | x | x | x | x
Asian Country (Dest.) |  |  |  |  | x
European Country (Dest.) | x | x | x | x |
Pacific Country (Dest.) | x | x | x | x | x
Battle Deaths (Orig.) | x |  | x | x | x
Battle Deaths (Dest.) | x | x | x | x | x
Homicide Rate (Dest.) |  |  |  |  | x
Posterior Model Probability | 0.007 | 0.005 | 0.004 | 0.004 | 0.003

Table A9: Top Five Highest Probability Models Using Hyper-$g/n$ Prior
(Bilateral Migration Data).

| Model #1 | Model #2 | Model #3 | Model #4 | Model #5
---|---|---|---|---|---
GDPPC (Orig.) | x | x | x | x | x
% Rural Pop. (Dest.) | x | x | x | x | x
Distance | x | x | x | x | x
Common Popular Language | x |  |  |  |
Population (Orig.) | x | x | x | x | x
Population (Dest.) | x | x | x | x | x
Gini Index (Orig.) |  |  |  |  | x
% Employed (Orig.) |  |  | x |  | x
% Employed (Dest.) | x | x | x | x | x
% Tertiary School Enrollment (Orig.) | x | x | x | x | x
% Tertiary School Enrollment (Dest.) | x | x | x | x | x
% Secondary School Enrollment (Orig.) |  |  | x |  |
% Secondary School Enrollment (Dest.) | x |  | x | x | x
Existing Migrant Population | x | x | x | x | x
% Working Age (Dest.) | x | x | x | x | x
Island Country (Dest.) | x | x | x | x | x
GDP Growth (Orig.) | x | x |  |  |
GDP Growth (Dest.) | x | x | x | x | x
Stability Index (Dest.) | x | x | x | x | x
Service % of GDP (Orig.) | x | x | x | x | x
Service % of GDP (Dest.) | x | x | x | x | x
Area (Dest.) | x | x | x | x | x
Landlocked Country (Dest.) | x | x | x | x | x
Latitude Capital (Dest.) | x | x | x | x | x
Longitude Capital (Orig.) | x | x | x | x | x
Longitude Capital (Dest.) | x | x | x | x | x
Asian Country (Orig.) | x | x | x | x |
European Country (Orig.) | x | x | x | x |
Pacific Country (Orig.) | x | x | x | x | x
Asian Country (Dest.) | x |  | x |  |
European Country (Dest.) |  | x |  | x | x
Pacific Country (Dest.) | x | x | x | x | x
Agriculture % of GDP (Orig.) | x | x |  |  |
Battle Deaths (Orig.) | x | x |  | x | x
Battle Deaths (Dest.) | x | x | x | x | x
Homicide Rate (Dest.) | x |  | x |  |
Posterior Model Probability | 0.001 | 0.001 | 0.001 | 0.001 | 0.001

###  A12.2 Visualizations of the Data and of the Model Averaging Results

Fig. A6: Map of Survey Clusters in Ethiopia DHS Survey 2019. Cluster
transparency is inversely proportional to local sample size.

(a) Post. Inclusion Probabilities.

(b) Coefficient Post. Means.

Fig. A7: Estimation Results for Measles Vaccination Data. Fig. A8: Bilateral
Migration Flows between OECD countries. Arrow size is proportional to size of
migration flow.

(a) Post. Inclusion Probabilities.

(b) Coefficient Post. Means.

Fig. A9: Estimation Results for Bilateral Migration Flow Data.

(a) Highest Probability Models.

(b) Posterior Distribution of Model Size.

Fig. A10: Estimation Results (Austrian Migration Flow Data Subset). Highest
probability models plot includes variables with estimated $\text{PIP}>0.1$
under unit information prior.

(a) Post. Inclusion Probabilities.

(b) Coefficient Post. Means.

Fig. A11: Estimation Results for Austrian Subset of Bilateral Migration Flow
Data.

(a) Highest Probability Models.

(b) Posterior Distribution of Model Size.

Fig. A12: Estimation Results (Measles Vaccination Data Subset). Highest
probability models plot includes variables with estimated $\text{PIP}>0.1$
under unit information prior.

(a) Post. Inclusion Probabilities.

(b) Coefficient Post. Means.

Fig. A13: Estimation Results for Low Vaccination Rate Subset of Measles
Vaccination Data.

(a) Full Vaccination Data (All Models).

(b) Full Migration Data (All Models).

(c) Full Vaccination Data (ULLGM Models).

(d) Full Migration Data (ULLGM Models).

Fig. A14: Boxplots of log predictive density evaluations across 100 random
training-test partitions for the full data sets. For the migration data, the
results are on a double log scale.

(a) Subset Vaccination Data (All Models).

(b) Subset Migration Data (All Models).

(c) Subset Vaccination Data (ULLGM Models).

(d) Subset Migration Data (ULLGM Models).

Fig. A15: Boxplots of log predictive density evaluations based on leave-one-
out cross-validation in the reduced samples. For the migration data, the
results are on a double log scale.

## Supplementary Material References

  * G. J. Abel and J. E. Cohen (2019) Bilateral international migration flow estimates for 200 countries.  Scientific Data 6 (1),  pp. 82.  Cited by: §A10.
  * J. H. Albert and S. Chib (1993) Bayesian analysis of binary and polychotomous response data.  Journal of the American statistical Association 88 (422),  pp. 669–679.  Cited by: §A2.
  * J. J. Azose and A. E. Raftery (2019) Estimation of emigration, return migration, and transit migration between all pairs of countries.  Proceedings of the National Academy of Sciences 116 (1),  pp. 116–122.  Cited by: §A10.
  * J. Bijak (2010) Forecasting international migration in europe: a bayesian view.  Vol. 24, Springer Science & Business Media.  Cited by: §A10.
  * C. Carrere (2006) Revisiting the effects of regional trade agreements on trade flows with proper specification of the gravity model.  European economic review 50 (2),  pp. 223–247.  Cited by: §A10.
  * A. B. Chan and N. Vasconcelos (2009) Bayesian Poisson regression for crowd counting.  In 2009 IEEE 12th international conference on computer vision,  pp. 545–551.  Cited by: §A11.
  * H. De Haas, S. Castles, and M. J. Miller (2019) The age of migration: international population movements in the modern world.  Bloomsbury Publishing.  Cited by: §A10.
  * C. Holmes and L. Held (2006) Bayesian auxiliary variable models for binary and multinomial regression.  Bayesian Analysis 1,  pp. 145––168.  Cited by: §A2.
  * J. Mitchell, N. Pain, and R. Riley (2011) The drivers of international migration to the uk: a panel-based bayesian model averaging approach.  The Economic Journal 121 (557),  pp. 1398–1444.  Cited by: §A10.
  * C. Morley, J. Rosselló, and M. Santana-Gallego (2014) Gravity models for tourism demand: theory and use.  Annals of tourism research 48,  pp. 1–10.  Cited by: §A10.
  * D. Owen (1980) A table of normal integrals.  Communications in Statistics-Simulation and Computation 9 (4),  pp. 389–419.  Cited by: §A1.
  * J. S. Silva and S. Tenreyro (2006) The log of gravity.  The Review of Economics and statistics 88 (4),  pp. 641–658.  Cited by: §A10.
  * N. G. Welch and A. E. Raftery (2022) Probabilistic forecasts of international bilateral migration flows.  Proceedings of the National Academy of Sciences 119 (35),  pp. e2203822119.  Cited by: §A10.