# In Nonparametric and High-Dimensional Models, Bayesian Ignorability is an
Informative Prior

Antonio R. Linero  Department of Statistics and Data Sciences, University of
Texas at Austin, email<EMAIL_ADDRESS>
###### Abstract

In problems with large amounts of missing data one must model two distinct
data generating processes: the outcome process which generates the response
and the missing data mechanism which determines the data we observe. Under the
_ignorability_ condition of Rubin, (1976), however, likelihood-based inference
for the outcome process does not depend on the missing data mechanism so that
only the former needs to be estimated; partially because of this
simplification, ignorability is often used as a baseline assumption. We study
the implications of Bayesian ignorability in the presence of high-dimensional
nuisance parameters and argue that ignorability is typically incompatible with
sensible prior beliefs about the amount of selection bias. We show that, for
many problems, ignorability directly implies that the prior on the selection
bias is tightly concentrated around zero. This is demonstrated on several
models of practical interest, and the effect of ignorability on the posterior
distribution is characterized for high-dimensional linear models with a ridge
regression prior. We then show both how to build high-dimensional models which
encode sensible beliefs about the selection bias and also show that under
certain narrow circumstances ignorability is less problematic.

##  1 Introduction

Dealing with missing data is a fundamental problem in data analysis; for
example, missingness complicates inference in clinical trials (National
Research Council,, 2010) and is inherent in the potential outcomes framework
for causal inference (Rubin,, 2005). A common starting point for addressing
missingness is to assume that the mechanism which generated the missingness is
_ignorable_ (Rubin,, 1976). Ignorability allows likelihood-based inference to
proceed without modeling the missing data mechanism, which can greatly
simplify an analysis.

In this paper we consider the Bayesian approach to account for missingness.
For generality, we consider a potential outcome $Y_{i}(a)$ for some exposure
level $a\in\mathscr{A}$ such that we observe both the exposure level $A_{i}$
and its associated potential outcome $Y_{i}(A_{i})$ ($Y_{i}(a)$ is regarded as
missing for all $a\neq A_{i})$. Let $X_{i}$ be a vector of confounders which
are predictive of both $A_{i}$ and $Y_{i}(a)$. By defining $Y_{i}(1)$ as the
outcome of interest, this framework subsumes the standard missing data
problem, where $A_{i}$ is now a missing data indicator such that we observe
the outcome when $A_{i}=1$. We say that the _missing data mechanism_
$f_{\phi}(A_{i}\mid X_{i})$ is Bayesian-ignorable, or simply ignorable, if the
following conditions hold.

  1. IG.1

The potential outcomes $\\{Y_{i}(a):a\in\mathscr{A}\\}$ are conditionally
independent of $A_{i}$ given $X_{i}$.

  2. IG.2

The parameters $\beta$ and $\phi$ are a-priori independent, where $\beta$
parameterizes the model for the potential outcomes and $\phi$ parameterizes
the missing data mechanism. That is, the prior factors as
$\pi(\beta,\phi)=\pi_{\beta}(\beta)\,\pi_{\phi}(\phi)$.

Condition IG.1 constrains the data generating mechanism and is a (type of)
_missing at random_ (MAR) assumption (Seaman et al.,, 2013); in the causal
inference literature, assumptions like IG.1 are sometimes themselves referred
to as strong ignorability assumptions (Rosenbaum and Rubin,, 1983; Imai et
al.,, 2010), and it is sometimes conflated with ignorability in missing data
sense of Rubin, (1976) as well (see Seaman et al.,, 2013, for a through
discussion of MAR and ignorability). Condition IG.2, which constrains the
prior, is also key to ignorability: it guarantees that the posterior
distribution of $\beta$ given the observed data is proportional to
$\pi_{\beta}(\beta)\ \prod_{i}f_{\beta}\\{Y_{i}(A_{i})\mid X_{i}\\}$, which
does not depend on the missing data mechanism. Without IG.2 we are still
obligated to model the missing data mechanism even when the missing data is
MAR, as $A_{i}$ provides information about $\beta$ through $\phi$.

It has been argued before that, from a Frequentist perspective, IG.2 is highly
problematic in high-dimensional problems (Robins and Ritov,, 1997). We
complement this Frequentist view and study the implications of IG.2 from a
Bayesian perspective. In particular we will argue that, while IG.2 is
seemingly innocuous, it actually is highly informative about the selection
bias in high-dimensional problems to the degree that the data has no
reasonable chance at overcoming the prior. We refer to this as _prior
dogmatism_ about the selection bias. We make the following three points.

  1. 1.

Priors which impose IG.2 are dogmatic about the amount of selection bias. This
is particularly true in models which require the use of informative priors,
such as high-dimensional or Bayesian nonparametric models. We conclude that
IG.2 does not reflect substantive prior knowledge in most cases; in the case
of a causal ridge regression model, we are able to quantify these problems
explicitly using random matrix theory (see Dobriban and Wager,, 2018 and
Dicker,, 2016 for related applications of random matrix theory).

  2. 2.

By understanding this induced prior on the selection bias, we are able to
identify several highly effective ways of correcting this problem and unify
several approaches proposed in the Bayesian causal inference literature which
were not motivated by Bayesian considerations. Our remedies take the form of
propensity score adjustments, which have typically been recommended in applied
Bayesian analysis on the grounds of pragmatism and robustness (see, e.g.,
Rubin,, 1985) rather than subjective Bayesian principles.

  3. 3.

We study some relatively narrow settings in which prior dogmatism does not
occur, even in high dimensional problems. For example, strong dependence
structure in $X_{i}$ can act as a shield against dogmatism; in the case of
causal ridge regression, we again use random matrix theory to quantify this
behavior. Despite this, we find little payoff for failing to correct for
dogmatism in these settings.

######  Remark 1.

While we will consider the Frequentist properties of the posterior
distribution, our goal at the outset is not to construct priors specifically
to attain ideal Frequentist properties. It is often quite easy to construct
priors which are doubly robust and attain some semiparametric efficiency bound
if that is our goal from the outset, and various complete class theorems
(Robert,, 2007, Chapter 8) suggest that we can usually construct _some_ Bayes
estimator which is at-least-as-good as any given Frequentist estimator.
Rather, we (i) show that priors of the form IG.2 are inherently dangerous in
purely-Bayesian terms, (ii) explain in which situations the problem is most
acute, and (iii) use dogmatism to show where corrections are needed.

###  1.1 Notation

We let $Y_{i}(\cdot)\in\mathbb{R}$ denote an outcome, $X_{i}\in\mathbb{R}^{P}$
a covariate/confounder, and $A_{i}\in\mathbb{R}$ a treatment/missing data
indicator for $i=1,\ldots,N$. When considering causal inference problems, we
define $Y_{i}=Y_{i}(A_{i})$ to be the observed outcome; when $A_{i}$ is a
missingness indicator, we instead define $Y_{i}=Y_{i}(1)$ so that $Y_{i}$ is
missing when $A_{i}\neq 1$. We set $\bm{Y}=(Y_{1},\ldots,Y_{N})^{\top}$,
$\bm{A}=(A_{1},\ldots,A_{N})$, and let $\bm{X}$ denote an $N\times P$ matrix
obtained by stacking the row vectors $X_{i}^{\top}$. Let $\beta$ parameterize
the distribution of $[Y_{i}(\cdot)\mid X_{i}]$, let $\phi$ parameterize the
distribution of $[A_{i}\mid X_{i}]$, and let $\theta=(\beta,\phi)$. We invoke
IG.1 throughout. We let $\mathbb{E}_{\theta}(\cdot)$ denote the expectation
operator conditional on $\theta$. If the subscript $\theta$ is omitted then
$\mathbb{E}(\cdot)$ is the expectation operator with respect to a prior
distribution on $\theta$, e.g.,
$\mathbb{E}(Y_{i})=\int\mathbb{E}_{\theta}(Y_{i})\,\pi(\theta)\ d\theta$. We
use the Big-O notation $W=O_{p}(V)$ to mean that $|W|/|V|$ is bounded in
probability as $P\to\infty$ (with the dependence of $W$ and $V$ on $P$
suppressed).

###  1.2 Three Illustrative Problems

We consider three problems to illustrate the existence of dogmatism and how to
correct for it. The first two concern causal inference with a continuous
exposure, while the third is a missing data problem. We assume
$X_{i}\sim\operatorname{Normal}(0,\Sigma)$ for some
$\Sigma\in\mathbb{R}^{P\times P}$ to simplify our analysis. All proofs are
deferred to the Supplementary Material.

##### Ridge regression in causal inference

Let $Y_{i}(a)$ denote the outcome observed when individual $i$ receives the
level $a$ of some continuous exposure and let $A_{i}$ denote the value of the
exposure which is actually received. We observe $(A_{i},Y_{i})$ where
$Y_{i}=Y_{i}(A_{i})$. We posit the linear models
$Y_{i}(a)=X_{i}^{\top}\beta+\gamma\,a+\epsilon_{i}(a)$ and
$A_{i}=X_{i}^{\top}\phi+\nu_{i}$ with
$\epsilon_{i}(a)\sim\operatorname{Normal}(0,\sigma^{2}_{y})$ and
$\nu_{i}\sim\operatorname{Normal}(0,\sigma^{2}_{a})$. The Bayesian ridge
regression prior, which satisfies IG.2, takes
$\beta\sim\operatorname{Normal}(0,\tau^{2}_{\beta}\ \mathrm{I})$ and
$\phi\sim\operatorname{Normal}(0,\tau^{2}_{\phi}\ \mathrm{I})$. The parameter
of inference is the mean response at a given exposure
$\mathbb{E}_{\theta}\\{Y_{i}(a)\\}=\gamma\,a$. We assume that $X_{i}$ is high-
dimensional in the sense that $P$ is potentially larger than $N$.

##### Sparsity priors in causal inference

This is the same problem as the ridge regression problem, except that $\beta$
and $\phi$ are sparse. We consider independent spike-and-slab priors for the
coefficients, i.e.,
$\beta_{j}\stackrel{{\scriptstyle\textnormal{iid}}}{{\sim}}(1-p_{\beta})\,\delta_{0}+p_{\beta}\,\operatorname{Normal}(0,\tau^{2}_{\beta})$
and
$\phi_{j}\stackrel{{\scriptstyle\textnormal{iid}}}{{\sim}}(1-p_{\phi})\,\delta_{0}+p_{\phi}\,\operatorname{Normal}(0,\tau^{2}_{\phi})$
where $\delta_{0}$ denotes a point-mass distribution at $0$.

##### Semiparametric regression with missing data

An outcome $Y_{i}=Y_{i}(1)$ is observed if $A_{i}=1$ and missing if $A_{i}=0$,
and our goal is to estimate $\mathbb{E}_{\theta}(Y_{i})$. We consider the
model $Y_{i}=\beta(X_{i})+\epsilon_{i}$ with
$\epsilon_{i}\sim\operatorname{Normal}(0,\sigma^{2}_{y})$ and
$A_{i}\sim\operatorname{Bernoulli}\\{\phi(X_{i})\\}$. We assume that
$\beta(\cdot)$ has a _Gaussian process_ prior (Rasmussen and Williams,, 2006)
with covariance function $\kappa(\cdot,\cdot)$, written
$\operatorname{GP}(0,\kappa)$.

##  2 The Induced Prior on the Selection Bias

The fundamental difficulty with missingness is _selection bias_. When
estimating $\mathbb{E}_{\theta}\\{Y_{i}(a)\\}$ this amounts to the fact that
$\Delta(a)=\mathbb{E}_{\theta}\\{Y_{i}(a)\mid
A_{i}=a\\}-\mathbb{E}_{\theta}\\{Y_{i}(a)\\}\neq 0$. That the selection bias
parameter $\Delta(a)$ is non-zero is the only feature of the problem which
makes estimation of $\mathbb{E}_{\theta}\\{Y_{i}(a)\\}$ non-trivial, as
otherwise we could ignore the covariates $X_{i}$ and directly estimate
$\mathbb{E}_{\theta}\\{Y_{i}(a)\\}$ by estimating
$\mathbb{E}_{\theta}\\{Y_{i}(a)\mid A_{i}=a\\}$ nonparametrically. The
following proposition gives an expression for $\Delta$ in each of our
problems.

######  Proposition 1.

In the ridge and spike-and-slab regression problems, the selection bias is
given by
$\Delta(a)=a\,\dfrac{\phi^{\top}\Sigma\beta}{\sigma^{2}_{a}+\phi^{\top}\Sigma\phi}=a\,\dfrac{\sum_{j}\lambda_{j}\,W_{j}\,Z_{j}}{\sigma^{2}_{a}+\sum_{j}\lambda_{j}\,Z_{j}^{2}}$
where $W=\Gamma^{\top}\beta$, $Z=\Gamma^{\top}\phi$, and
$\Sigma=\Gamma\Lambda\Gamma^{\top}$ is the spectral decomposition of $\Sigma$
with $\Lambda=\operatorname{diag}(\lambda_{1},\ldots,\lambda_{P})$. In the
semiparametric regression problem with missing data, we instead have
$\Delta\equiv\Delta(1)={\operatorname{Cov}_{\theta}\\{\beta(X_{i}),\phi(X_{i})\\}}/{\mathbb{E}_{\theta}\\{\phi(X_{i})\\}}.$

Given the importance $\Delta$ and the working assumption that selection bias
is non-negligible, one would hope that the prior distribution of $\Delta$ is
relatively diffuse. Using Proposition 1 we can gain insight into how the prior
on the selection bias changes as the dimension $P$ increases. For example, for
the ridge regression problem we have the following result.

######  Proposition 2.

Assume the setup of Proposition 1 for the ridge regression problem and suppose
$\beta\sim\operatorname{Normal}(0,\tau^{2}_{\beta}\,\mathrm{I})$ and
$\phi\sim\operatorname{Normal}(0,\tau_{\phi}^{2}\,\mathrm{I})$ independently.
Assume $\frac{1}{P}\sum_{j=1}^{P}\lambda_{j}^{k}$ converges to a positive
constant as $P\to\infty$ for $k=1,2,2+\epsilon$ for some $\epsilon$, and let
$\widetilde{\lambda}$ and $\bar{\lambda}^{2}$ be the limits with $k=1,2$. Then
$\Delta(a)\stackrel{{\scriptstyle{}_{\bullet}}}{{\sim}}\operatorname{Normal}(0,c/P)$
where
$c=a^{2}\,(\tau^{2}_{\beta}/\tau^{2}_{\phi})\,(\bar{\lambda}^{2}/\widetilde{\lambda}^{2})$.

We will return to the conditions on $\Sigma$ (which are moment conditions on
the spectral distribution of $\Sigma$) later and focus on the conclusion
$\Delta(a)\stackrel{{\scriptstyle{}_{\bullet}}}{{\sim}}\operatorname{Normal}(0,c/P)$
for some constant $c$. If selection bias is a-priori of concern for us then it
seems unwise to posit a $\operatorname{Normal}(0,c/P)$ prior for it when $P$
is large. This behavior becomes even more suspect when one considers that the
definition of $\Delta(a)$ is completely free of the $X_{i}$’s, and that
logically the act of measuring additional covariates should not change our
beliefs about $\Delta(a)$. In Section 2.1 we follow up on the inferential
consequences of this.

At a high level, the source of the problem in our three illustrative examples
is the following well-known phenomenon which we refer to as the _orthogonality
principle_.

######  Principle 1 (The Orthogonality Principle).

Let $\widetilde{\beta}$ and $\widetilde{\phi}$ be random unit vectors with
mean $0$ taking values in some high/infinite dimensional Hilbert space
$\mathcal{H}$ with inner product $\langle\cdot,\cdot\rangle$. Then, if
$\widetilde{\phi}$ and $\widetilde{\beta}$ are independent and there is no
other special structure in the problem, with high probability we have
$\langle\widetilde{\beta},\widetilde{\phi}\rangle\approx 0$.

For a concrete example, by the law of large numbers and the central limit
theorem, if
$\beta_{j},\phi_{j}\stackrel{{\scriptstyle\textnormal{iid}}}{{\sim}}\operatorname{Normal}(0,1)$
then $\beta^{\top}\phi/\sqrt{P}\to\operatorname{Normal}(0,1)$ but
$\|\beta\|\,\|\phi\|/P\to 1$ so that
$\langle\beta/\|\beta\|,\phi/\|\phi\|\rangle=O_{p}(P^{-1/2})$ with respect to
the Euclidean inner product. The statement of the orthogonality principle is
intentionally vague as to what it means for the dimension to be “high,” what
$\approx$ means, and what constitutes “special structure.” Nevertheless, it
provides immediate intuition for what to expect: unless one has reason to
believe otherwise, expect high-dimensional unit vectors to be nearly
orthogonal.

The orthogonality principle becomes important when $P$ is large (or in
nonparametric problems) because $\Delta$ is quantifiable in terms of
$\langle\beta,\phi\rangle$ for some suitable inner product (c.f. Proposition
1). If IG.2 holds then the orthogonality principle immediately suggests
$\langle\beta,\phi\rangle\approx 0$ with high probability, implying that our
prior is dogmatic about the selection bias.

###  2.1 Asymptotics for High-Dimensional Ridge Regression

While the dogmatism implied by Proposition 2 is troubling, one might hope that
the informative prior on $\Delta$ is a theoretical curiosity which is
nevertheless swamped by the data. By analogy, a strict prior analysis of the
flat prior $\mu\sim\operatorname{Normal}(0,10^{100})$ in the normal model
$Y_{i}\stackrel{{\scriptstyle\textnormal{iid}}}{{\sim}}\operatorname{Normal}(\mu,1)$
would similarly suggest that we ought to be concerned about the implication of
the prior on the magnitude of $\mu$; instead the data quickly swamps the
diffuse prior to produce a sensible posterior
$\mu\approx\operatorname{Normal}(\bar{Y},1/N)$. That is, in terms of
consequences for inference, what matters is the impact of IG.2 on the
posterior rather than the prior.

We now show that this hopeful scenario is not borne out, and that the prior
concentration on $\Delta$ leads to heavily biased inferences if $P$ grows
sufficiently quickly with $N$. We summarize our main results as follows:

  * •

In the regime $P/N\to r$ for some $r\in(0,\infty)$ (i.e., $P$ grows at the
same rate as $N$), the Bayes estimator which takes a flat prior on $\gamma$
and a Gaussian prior
$\beta\sim\operatorname{Normal}(0,\tau^{2}\,P^{-1}\,\mathrm{I})$ is heavily
biased. Specifically, when selection bias is present through the auxiliary
covariate $\widehat{A}_{i}=X_{i}^{\top}\phi$, the Bayes estimate will have
bias of order $\Delta(1)$.

  * •

In some sense, the setting $\Sigma=\mathrm{I}$ is inherently difficult, and
the problem is generally easier when the components of $X_{i}$ are highly
correlated. We return to this point in Section 5.

We make two sets of assumptions. The first (high-dimensional asymptotics, or
HDA) is used to describe the distribution of the $X_{i}$’s as $N\to\infty$.
The second (random effects model, or REM) describes a particular random
effects model for the regression coefficients. This framework modifies the
framework of Dobriban and Wager, (2018) so that it is suitable for our aims.

  1. HDA.1

The covariates are multivariate normal with
$X_{i}\sim\operatorname{Normal}(0,\Sigma)$.

  2. HDA.2

As $N\to\infty$ we have $P/N\to r$ for some $r\in(0,\infty)$.

  3. HDA.3

The spectral distribution $\sum_{p=1}^{P}\delta_{\lambda_{p}}/P$ associated to
$\Sigma$ converges to some limiting distribution $H$ on $[0,\infty)$, where
$\lambda_{1},\ldots,\lambda_{P}$ are the eigenvalues of $\Sigma$ and
$\delta_{\lambda}$ denotes a point-mass distribution at $\lambda$.

HDA is a standard assumption for understanding the case where $P$ grows like
$N$, though HDA.1 may be replaced with a moment condition on $X_{i}$. HDA.3
allows us to use results from random matrix theory to compute
$\lim_{P\to\infty}\operatorname{tr}\\{(\bm{X}^{\top}\bm{X}+N\lambda\,\mathrm{I})^{-k}\\}$
for $k\in\mathbb{N}$. Under HDA, the empirical distribution of the eigenvalues
of $\underline{S}=\bm{X}\bm{X}^{\top}/N$, namely
$\widehat{F}(dx)=\sum_{i=1}^{N}\delta_{\widehat{\lambda}_{i}}/N$, converges to
a distribution $F(dx)$ called the _empirical spectral distribution_.

Next, we describe a random effects model (REM) for $\beta$ and $\phi$ we will
base our analysis on. Similar models have been used to study both the
prediction risk and minimax-optimality of ridge regression (Dobriban and
Wager,, 2018; Dicker,, 2016). REM is a fruitful assumption for us as it allows
exact formulas for the bias to be derived which are free of the particular
values of $\beta$ and $\phi$.

  1. REM.1

The coefficient vector $\phi$ is randomly sampled as
$\phi\sim\operatorname{Normal}(0,\tau^{2}\,P^{-1}\,\mathrm{I})$.

  2. REM.2

The coefficient vector $\beta$ is randomly sampled as
$\beta\sim\operatorname{Normal}(\omega_{0}\,\phi,\tau^{2}\,P^{-1}\,\mathrm{I})$.

  3. REM.3

Given $\beta$ and $\phi$,
$Y_{i}\sim\operatorname{Normal}(X_{i}^{\top}\beta+A_{i}\,\gamma_{0},1)$ and
$A_{i}\sim\operatorname{Normal}(X_{i}^{\top}\phi,1)$.

To motivate REM.2, note that it is equivalent to setting
$Y_{i}\sim\operatorname{Normal}(X_{i}^{\top}b+\omega_{0}\,\widehat{A}_{i}+\gamma_{0}\,A_{i},1)$
where $\widehat{A}_{i}=X_{i}^{\top}\phi=\mathbb{E}(A_{i}\mid X_{i},\phi)$ and
$b\sim\operatorname{Normal}(0,\tau^{2}\,P^{-1}\,\mathrm{I})$. REM.2 allows for
non-negligible selection bias to enter the model, and priors based on this
parameterization have been used to account for selection bias by other
researchers (Zigler et al.,, 2013; Hahn et al.,, 2018). The parameter
$\omega_{0}$ is intimately connected to the selection bias.

######  Proposition 3.

Suppose that HDA and REM hold and that $\Sigma$ satisfies the conditions of
Proposition 2. Then
$\Delta(1)\to\omega_{0}\frac{\tau^{2}\,\widetilde{\lambda}}{1+\tau^{2}\,\widetilde{\lambda}}$
in probability as $P\to\infty$.

Theorem 1 explicitly computes the bias of the ridge regression estimator under
IG.2 when the prior
$\beta\sim\operatorname{Normal}(0,N^{-1}\lambda^{-1}\mathrm{I})$ is used,
i.e., when we apply the usual ridge regression estimator. We sketch a proof of
Theorem 1 and verify it numerically in the Supplementary Material.

######  Theorem 1.

Suppose HDA and REM hold. Let
$(\widetilde{\gamma},\widetilde{\beta}^{\top})^{\top}$ denote the Bayes
estimate of $(\gamma,\beta^{\top})^{\top}$ under a prior which takes
$\beta\sim\operatorname{Normal}(0,N^{-1}\,\lambda^{-1}\mathrm{I})$ and places
a flat prior on $\gamma$ under IG.2. Then the asymptotic bias of
$\widetilde{\gamma}$ is given by

$\displaystyle\lim_{N,P\to\infty}\mathbb{E}(\widetilde{\gamma}-\gamma_{0})=\frac{\omega_{0}\int
x/(x+\lambda)\ F(dx)}{\int(x+\eta)/(x+\lambda)\
F(dx)}=\omega_{0}\times\frac{1-\lambda\,v(-\lambda)}{1-(\lambda-\eta)\,v(-\lambda)}$
(1)

where $v(z)=\int_{0}^{\infty}\frac{F(dx)}{x-z}$ is the Stieltjes transform of
$F(dx)$ and $\eta=r/\tau^{2}$.

Ideally we would like the bias to be close to $0$ for moderate-to-large values
of $\lambda$ so that we have both small variance and bias; the approach
outlined in Section 4.1 _does_ accomplish this goal for a properly-chosen
$\lambda$. Figure 1 contrasts this alternative method with standard ridge
regression when $\Sigma=\mathrm{I}$ and we see that the bias is quite large
for ridge regression unless $\lambda$ is close to $0$ and $r\leq 1$; this
latter case corresponds to OLS, which (while unbiased) defeats the purpose of
using ridge regression.

Figure 1: Comparison of the bias of naive ridge regression (dashed, blue) to
the direct Z-prior (solid, orange) of Section 4.1 for different values of
$\eta$ and $r$ with $\omega_{0}\equiv 1$. Estimated values of $\lambda$ based
on a single simulated dataset for each combination of $\eta$ and $r$ are given
by the points.

A qualitative observation based on (1) is that smaller bias is obtained when
most of the eigenvalues of $\underline{S}$ are small. For example,
unbiasedness is possible if $F(dx)$ assigns _any_ mass to $0$, since taking
$\lambda\to 0$ will cause $\lambda\,v(-\lambda)\to 0$ (by bounded convergence)
while $\eta\,v(-\lambda)\to\infty$ in the denominator. This occurs naturally
when $N>P$, and incidentally would also occur if duplicate rows of $\bm{X}$
were possible even if $P\gg N$; this latter observation makes intuitive sense,
as we could then identify $\gamma_{0}$ using exact-matching on the $X_{i}$’s.

When $P>N$ the only hope for non-negligible bias is for the eigenvalues of
$\underline{S}$ to be heavily concentrated near $0$. As $\underline{S}$ has
the same non-zero eigenvalues as the sample covariance
$S=\bm{X}^{\top}\bm{X}/N$ this means we should hope for strong colinearities
among the covariates. A particularly unfavorable setting is
$\Sigma=\mathrm{I}$, where the Marchenko-Pastur theorem (see, e.g., Couillet
and Debbah,, 2011, Theorem 2.13) states that if $r\geq 1$ then $F(dx)$ has
density
$q(\lambda)=\frac{\sqrt{(b-\lambda)(\lambda-a)}}{2\,\pi\,\lambda}I(a<\lambda<b)$
where $(a,b)=(1\pm\sqrt{r^{-1}})^{2}$; this places the bulk of the eigenvalues
rather far from $0$. In Section 5 we show that much better results are
obtained when the $X_{i}$’s follow a latent factor model.

###  2.2 Selection Bias Dogmatism for Semiparametric Regression

In the semiparametric regression problem with missing data the selection bias
parameter is given by
$\Delta=\frac{\operatorname{Cov}_{\theta}\\{\beta(X_{i}),\phi(X_{i})\\}}{\mathbb{E}_{\theta}\\{\phi(X_{i})\\}}.$
Figure 2 gives a sense of what to expect for nonparametric priors. In this
figure, $\beta(x)$ and $\Phi^{-1}\\{\phi(x)\\}$ are given independent BART
priors (Chipman et al.,, 2010; Hill,, 2011) where $\Phi(\cdot)$ is a probit
function. We see that as $P$ increases the variance of $\Delta$ decreases
substantially. As in the setting of ridge regression, this is troubling both
because (i) it will typically violate our prior beliefs about $\Delta$ for
large $P$ and (ii) given the definition of $\Delta$ as
$\mathbb{E}_{\theta}(Y_{i}\mid A_{i}=1)-\mathbb{E}_{\theta}(Y_{i})$ there is
no reason for our prior beliefs to be dependent on the number of confounders
we happen to have measured.

Figure 2: Prior distribution of $\Delta$ for the BART model in Section 2.2 for
$P\in\\{1,10,50\\}$.

For convenience we will assume that $\phi(x)$ has a point-mass prior at some
$\phi_{0}$ and that $\beta$ has a Gaussian process prior (Rasmussen and
Williams,, 2006). Recall that $\beta\sim\operatorname{GP}(m,\kappa)$ means
that, for any finite collection $(x_{1},\ldots,x_{M})$, we have
$\big{(}\beta(x_{1}),\ldots,\beta(x_{M})\big{)}^{\top}\sim\operatorname{Normal}(\bm{m},\bm{K})$
where $\bm{m}=\big{(}m(x_{1}),\ldots,m(x_{M})\big{)}^{\top}$ and $\bm{K}$ has
$(j,k)^{\text{th}}$ entry $\kappa(x_{j},x_{k})$. Gaussian processes have been
proposed as priors for causal inference by several authors (Ray and van der
Vaart,, 2020; Ren et al.,, 2021) and they are particularly easy to study
theoretically.

The relevant Hilbert space for applying the orthogonality principle is
$\mathscr{L}_{2}(F_{X})$, the space of square-integrable functions $\\{g:\int
g^{2}\ dF_{X}<\infty\\}$ under the usual inner product
$\langle\beta,\phi\rangle=\int\beta(x)\,\phi(x)\ F_{X}(dx)$, with $F_{X}$
denoting the distribution of $X_{i}$. Let
$\bar{\beta}(x)=\beta(x)-\int\beta(x)\ F_{X}(dx)$ and
$\bar{\phi}(x)=\phi(x)-\int\phi(x)\ F_{X}(dx)$, and define the normalizations
of these functions by $\widetilde{\beta}(x)=\bar{\beta}(x)/\|\bar{\beta}\|$
and $\widetilde{\phi}(x)=\bar{\phi}(x)/\|\bar{\phi}\|$. The following
proposition shows that the selection bias is controlled by
$\langle\widetilde{\beta},\widetilde{\phi}\rangle$, implying that the
orthogonality principle is in effect.

######  Proposition 4.

Suppose $\beta\sim\operatorname{GP}(m,\kappa)$ such that
$\sup_{P}\mathbb{E}\\{\beta(X_{i})^{2}\\}<\infty$ and that there exists
$\delta>0$ such that $\mathbb{E}\\{\phi(X_{i})\\}\geq\delta$ as $P\to\infty$.
Then
$\Delta=\frac{\|\bar{\beta}\|\,\|\bar{\phi}\|}{\mathbb{E}_{\theta}\\{\phi(X_{i})\\}}\
\langle\widetilde{\beta},\widetilde{\phi}\rangle=O_{p}(\langle\widetilde{\beta},\widetilde{\phi}\rangle).$

The question now is how quickly
$\langle\widetilde{\beta},\widetilde{\phi}\rangle$ tends to $0$. In the
nonparametric case, in addition to the dimension ($P$) and distribution of the
covariates ($F_{X}$), this will also depend on the _smoothness_ of
$\widetilde{\beta}$ dictated by the covariance function $\kappa(\cdot,\cdot)$.
Note that $\bar{\beta}=\beta-\int\beta\ dF_{X}$ is also a Gaussian process
with covariance function

$\displaystyle\bar{\kappa}(x,x^{\prime})=\kappa(x,x^{\prime})-\int\kappa(x,z)\
F_{X}(dz)-\int\kappa(x^{\prime},z)\ F_{X}(dz)+\iint\kappa(x,x^{\prime})\
F_{X}(dx)\ F_{X}(dx^{\prime}).$

The following proposition explicitly calculates the prior distribution of
$\Delta$.

######  Proposition 5.

Let $\beta\sim\operatorname{GP}\\{0,\tau^{2}_{\beta}\,\rho(\cdot,\cdot)\\}$
where $\rho(\cdot,\cdot)$ is a correlation function. Then
$\Delta\sim\operatorname{Normal}(0,c)$ where $c$ is

$\displaystyle\frac{\tau^{2}_{\beta}}{\mathbb{E}\\{\phi(X_{i})\\}^{2}}\iint\bar{\phi}(x)\,\bar{\phi}(x^{\prime})\,\bar{\rho}(x,x^{\prime})\
F_{X}(dx)\
F_{X}(dx^{\prime})=\frac{\tau^{2}_{\beta}}{\mathbb{E}\\{\phi(X_{i})\\}^{2}}\sum_{j=1}^{\infty}\lambda_{j}\operatorname{Cov}\\{\phi(X_{i}),v_{j}(X_{i})\\}^{2}$

and
$\rho(x,x^{\prime})=\sum_{j=1}^{\infty}\lambda_{j}\,v_{j}(x)\,v_{j}(x^{\prime})$
is the Karhunen–Loève expansion of $\rho(x,x^{\prime})$ in
$\mathscr{L}_{2}(F_{X})$.

We see that the only way for $c$ to be non-negligible is for $\phi$ to be
highly correlated with some of the leading eigenfunctions of $\rho$; this
suggests, at a minimum, that the kernel $\rho(x,x^{\prime})$ should not be
chosen in a manner which does not reference $\phi(\cdot)$. The next
proposition shows that, should we not incorporate $\phi$ into
$\rho(x,x^{\prime})$, the resulting kernel can be universally poor: the
selection bias can decay exponentially in $P$ irrespective of $\phi$.

######  Proposition 6.

Consider the setup of Proposition 5 where the covariance function
$\rho(x,x^{\prime})$ is given by the Gaussian kernel
$\rho(x,x^{\prime})=\exp\\{-(x-x^{\prime})^{\top}H^{-1}(x-x^{\prime})/2\\}$
for some bandwidth matrix $H$ and suppose (i)
$X_{i}\stackrel{{\scriptstyle\textnormal{iid}}}{{\sim}}\operatorname{Normal}(0,\Sigma)$,
(ii) $\det(\Sigma)^{1/P}/\det(H)^{1/P}$ is bounded away from $0$, and (iii)
there exists a $\delta>0$ such that $\mathbb{E}\\{\phi(X_{i})\\}\geq\delta$
for all $P$. Then $\Delta\sim\operatorname{Normal}(0,c)$ where
$c\leq\exp\\{-CP\\}$ for some constant $C>0$ which is independent of $\phi$.
In particular, this occurs if either:

  1. (a)

$H=k\,\Sigma$ for some $k>0$ with $k^{1/P}$ bounded and $\Sigma$ full-rank; or

  2. (b)

$H=\xi\,\mathrm{I}$ and $\xi^{-1}\prod_{j}\lambda_{j}^{1/P}$ is bounded away
from $0$.

Proposition 6 shows that the Gaussian kernel is dogmatic in a _uniform_ sense:
no matter how favorably $\phi(x)$ is selected, the Gaussian kernel makes the
prior variance on $\Delta$ decrease exponentially in $P$. Moreover, this
exponential decay holds for some potentially-desirable choices of the
bandwidth matrix $H$ (both when $H$ is aligned with $\Sigma$ and when the
kernel is isotropic).

###  2.3 Selection Bias Dogmatism for Spike-and-Slab Priors

A common strategy for dealing with the $N\ll P$ setting in linear regression
is to use a sparsity inducing spike-and-slab prior like the one described in
Section 1.2. Even when we impose sparsity, however, serious problems occur for
the selection bias prior. Suppose that $\Sigma=\sigma^{2}_{x}\mathrm{I}$ and
let $\mathfrak{d}_{j}^{\beta}=I(\beta_{j}\neq 0)$ and
$\mathfrak{d}^{\phi}_{j}=I(\phi_{j}\neq 0)$. Then we can write the selection
bias as
$\Delta(a)=a\frac{\sum_{j:\mathfrak{d}^{\beta}_{j}=\mathfrak{d}^{\phi}_{j}=1}\sigma^{2}_{x}\,\phi_{j}\,\beta_{j}}{\sigma^{2}_{a}+\sum_{j:\mathfrak{d}^{\phi}_{j}=1}\sigma^{2}_{x}\,\phi_{j}^{2}}.$
In this case, the denominator (by the law of large numbers) will be of order
$\sum_{j}\mathfrak{d}_{j}^{\phi}\equiv D_{\phi}$ while the numerator (by the
central limit theorem) will be of order $D_{\phi\cap\beta}^{1/2}$ where
$\sum_{j}\mathfrak{d}_{j}^{\phi}\,\mathfrak{d}_{j}^{\beta}\equiv
D_{\phi\cap\beta}$. If we now use independent spike-and-slab priors for
$\beta$ and $\phi$ which are calibrated to have on average $Q$ variables, we
expect $D_{\phi}\approx Q$ while $D_{\phi\cap\beta}\approx Q/P$, so that the
selection bias will be a-priori negligible in high-dimensional sparse settings
in which spike-and-slab priors are applied. Hence, even if sparsity is
expected (but IG.2 is otherwise in effect) we run into essentially the same
problem as with ridge regression: the prior on $\Delta(a)$ regularizes it
towards zero.

##  3 Specifying Priors Which Violate IG.2: Z-Priors

Part of the appeal of IG.2 is that the treatment $A_{i}$ plays no role in
posterior sampling of the parameter of interest $\beta$. This is
computationally convenient because the updates for $\beta$ and $\phi$ in (say)
a Markov chain Monte Carlo experiment can be carried out independently. It
also prevents the phenomenon of _model feedback_ from occurring, wherein
misspecification of the $A_{i}$-model can result in inconsistent estimation of
$\beta$ even when the $Y_{i}$-model is correctly specified (Robins and Ritov,,
1997; Zigler et al.,, 2013).

Arguably the most natural way to specify a Bayesian model violating IG.2 is to
use the factorization

$\displaystyle\pi(\phi)\,\pi(\beta\mid\phi)\,f_{\phi}(\bm{A}\mid\bm{X})\,f_{\beta}(\bm{Y}\mid\bm{A},\bm{X}).$
(2)

That is, we specify the model in the usual way, but make $\beta$ dependent on
$\phi$. An alternative approach, which is less natural but far more
convenient, is to use the model specification

$\displaystyle\pi(\phi,\beta,\bm{A},\bm{Y}\mid\bm{X})=\pi(\phi)\,f_{\phi}(\bm{A}\mid\bm{X})\,\pi(\beta\mid\bm{A},\bm{X})\,f_{\beta}(\bm{Y}\mid\bm{A},\bm{X}).$
(3)

This approach is argued for by Hahn et al., (2020) who refer to specifications
like (3) as _Zellner priors_ due to the fact that the prior for $\beta$ is
allowed to depend on the design matrix $(\bm{A},\bm{X})$ of $\bm{Y}$ (as is
the case for Zellner’s famous $g$-prior). In order to avoid confusion with
Zellner’s $g$-prior, we refer to priors of this form as Z-priors. Figure 3
shows a schematic comparing these two factorizations with IG.2. This prior
makes $\beta$ and $\phi$ _conditionally independent_ given $(\bm{A},\bm{X})$
so that it retains the principle advantage of IG.2: feedback between the
$Y_{i}$ and $A_{i}$ models is severed, and the updates for $(\beta,\phi)$ are
no longer coupled. Note that (3) induces a dependent prior of the form
$\pi(\phi,\beta)=\pi(\phi)\,\int\pi(\beta\mid\bm{A},\bm{X})\,f_{\phi}(\bm{A}\mid\bm{X})\
d\bm{A}\ dF_{X},$ so that IG.2 is violated.

Figure 3: Directed acyclic graphs showing different conditional independence
structures for model and prior specification; (a) shows the graph implied by
IG.2, (b) shows the graph implied by (2), and (c) shows the graph implied by
(3).

To illustrate the point, in Section 4.1 we will consider a prior for the ridge
regression problem which is of the form
$\beta\sim\operatorname{Normal}(\omega\,\phi,\tau^{2}_{\beta}\,\mathrm{I})$
where $\omega$ is given a diffuse prior. This prior conforms to (2). In our
actual experiments, however, we use the prior
$\beta\sim\operatorname{Normal}(\omega\,\widehat{\phi},\tau^{2}_{\beta}\,\mathrm{I})$
where $\widehat{\phi}=\mathbb{E}(\phi\mid\bm{A},\bm{X})$ is a data-adaptive
ridge estimator. This prior is of the form (3) and is trivial to implement in
a two-stage fashion: fit the model for $A_{i}$, compute $\widehat{\phi}$, and
plug this into the prior for $\beta$ when fitting the $Y_{i}$-model.

In our experience, some Bayesians feel uneasy about using priors like
$\beta\sim\operatorname{Normal}(\omega\widehat{\phi},\tau^{2}_{\beta})$
because it “understates the uncertainty” in $\beta$ due to the fact that it
appears to use a plug-in estimate of $\phi$ rather than $\phi$ itself. While
it is true that using the actual value of $\phi$ rather than an estimate
typically (although not always! see Hirano et al.,, 2003) results in improved
Frequentist performance, the justification of this prior as being of the form
(3) shows that there is no explicit violation of the Bayesian calculus in
using this prior. For example, for the ridge regression Z-prior discussed
above we can explicitly derive the induced prior on
$\pi(\beta,\phi,\omega\mid\bm{X})$ as
$\pi(\beta,\phi,\omega\mid\bm{X})=\pi(\phi,\omega)\,\operatorname{Normal}\big{(}\beta\mid\omega\,\phi,\tau^{2}_{\beta}\,\mathrm{I}+\sigma^{2}_{a}(\bm{X}^{\top}\bm{X})^{-1}\big{)}.$

##  4 Correcting for Dogmatism

###  4.1 Direct Priors for Ridge Regression

A simple approach to addressing dogmatism in the ridge regression setting is
to note that we can make $\beta^{\top}\Sigma\phi$ large by encouraging $\beta$
to align with $\phi$. For example, we might center $\beta$ on $\phi$ by taking
$\beta\sim\operatorname{Normal}(\omega\phi,\tau^{2}_{\beta}\,\mathrm{I}).$
Doing this, we now have $\Delta(a)=a\frac{\phi^{\top}\Sigma
b}{\sigma^{2}_{a}+\phi^{\top}\Sigma\phi}+a\,\omega\frac{\phi^{\top}\Sigma\phi}{\sigma^{2}_{a}+\phi^{\top}\Sigma\phi},$
where $b\sim\operatorname{Normal}(0,\tau^{2}_{\beta}\,\mathrm{I})$. By the
same argument as in Proposition 2, the first term is $O_{p}(P^{-1/2})$; the
second term, however, does not tend to $0$ as $P\to\infty$, preventing prior
dogmatism from taking hold. By centering the prior for $\beta$, we can now
specify a _direct_ prior on $\Delta(a)$ by placing a prior on $\omega$. For
example, we can express prior ignorance about the degree of selection bias by
placing a flat prior on $\omega$.

This approach is related to the targeted maximum likelihood estimation
strategy of introducing a “clever covariate” into the outcome model to account
for selection (see, e.g., van der Laan and Rose,, 2011, Section 4.2.1). The
parameterization $\beta=b+\omega\phi$ gives
$Y_{i}(a)=\beta_{0}+X_{i}^{\top}b+\omega(X_{i}^{\top}\phi)+\gamma\,a+\epsilon_{i}(a),$
which effectively introduces the new covariate
$\widehat{A}_{i}=X_{i}^{\top}\phi$ into the model. A related idea proposed by
(Hahn et al.,, 2018) is to replace $a$ in the outcome model with the residual
$(a-\widehat{A}_{i})$, which is equivalent to setting $\omega=-\gamma$.

#### Bias of the Z-prior Estimate under HDA and REM

In practice, we will not usually use the direct prior described above;
instead, we will use a Z-prior which plugs in a point-estimate of the clever
covariate $\widehat{A}_{i}=X_{i}^{\top}\widehat{\phi}$. Assuming HDA and REM,
the bias of the Bayes estimator under the Z-prior can be shown to be

$\displaystyle\mathbb{E}(\widetilde{\gamma}-\gamma_{0})=\mathbb{E}\left\\{\frac{\bm{A}^{\top}(\mathrm{I}-\widehat{\Psi}\bar{H}\widehat{\Psi}^{\top})\Psi\binom{\omega_{0}}{b}}{\bm{A}^{\top}(\mathrm{I}-\widehat{\Psi}\bar{H}\widehat{\Psi}^{\top})\bm{A}}\right\\}$

where $\widehat{\Psi}=[\bm{X}\widehat{\phi},\bm{X}]$,
$\Psi=[\bm{X}\phi,\bm{X}]$ and
$\bar{H}^{-1}=\widehat{\Psi}^{\top}\widehat{\Psi}+N\lambda\left(\begin{smallmatrix}0&0\\\
0&\mathrm{I}\end{smallmatrix}\right)$. If one happened to know the exact value
of $\phi$ and set $\widehat{\phi}=\phi$ then it is easy to show that the
resulting Bayes estimator $\widetilde{\gamma}$ is unbiased for $\gamma$,
irrespective of the prior for $b$. It is also easy to show that if we take
$\widehat{\phi}=\widetilde{\phi}$ where $\widetilde{\phi}$ is the Bayes
estimate under the correct prior
$\widetilde{\phi}=\mathbb{E}(\phi\mid\bm{A},\bm{X})$ then $\widetilde{\gamma}$
remains unbiased.

The selling point now is that there are moderate values of $\lambda$ for which
ridge regression will have $0$ bias, a situation which was not possible with
the naive prior. In practice, under REM we will know neither $\phi$ nor
$\widetilde{\phi}$ (because we will not know the signal level
$\tau_{\phi}^{2}$). Data is typically quite informative about
$\tau^{2}_{\phi}$, however, and we have had success placing a prior on
$\tau^{2}_{\phi}$ to obtain nearly unbiased estimates. This is seen in Figure
1 ($\Sigma=\mathrm{I}$), where the point on the solid line corresponds to the
bias if we plug in a Bayes estimate of $\tau^{2}_{\phi}$ to construct a ridge
estimator.

It is also possible to show that the asymptotic bias of the Z-prior when
$\widehat{\phi}$ is estimated with ridge regression is

$\displaystyle\left(\omega_{0}\,\lambda\left[\frac{\psi_{11}}{\eta}-\frac{\psi_{22}}{\eta}\frac{(\psi_{21}+\psi_{22}/\eta)}{(\psi_{32}+\psi_{33}/\eta)}\right]\right)/\left({\frac{1-r}{r}+\lambda\left[\psi_{10}+\frac{\psi_{11}}{\eta}-\frac{(\psi_{21}+\psi_{22}/\eta)^{2}}{(\psi_{32}+\psi_{33}/\eta)}\right]}\right)$

where $\psi_{jk}=\int_{0}^{\infty}\frac{x^{k}}{(x+\lambda)^{j}}\ G(dx)$ and
$G(dx)$ is the empirical spectral distribution corresponding to
$\bm{X}^{\top}\bm{X}/N$ (also known as the companion spectral distribution to
$F(dx)$). Each of the $\psi_{jk}$’s can be computed by noting the recursive
identity $\psi_{jk}=\psi_{j-1,k-1}-\lambda\psi_{j,k-1}$ and the fact that
$\psi_{j0}=m^{(j-1)}(-\lambda)/(j-1)!$ where $m(z)$ is the Stieltjes transform
of $G(dx)$.

A plot of the bias is given in Figure 1 when $\Sigma=\mathrm{I}$, and we
numerically our bias formula in the Supplementary Material.Rather curiously,
for $r>1$ we see that the bias of the Z-prior estimate is ill-behaved near
$0$; fortunately, we can estimate $\tau^{2}$ accurately enough to avoid this
region.

#### Evaluation of the Direct Z-Prior via Simulation

In order to determine if there is any benefit to using the Z-prior with the
prior on $\omega$ being $\operatorname{Uniform}(-\infty,\infty)$ relative to
either (i) the naive ridge regression prior or (ii) the approach of (Hahn et
al.,, 2018) which we call the “debiased” approach (equivalent to fixing
$\omega=-\gamma$), we conducted a simulation study. In all cases we set
$N=200$ and $P=1000$ so that $N\ll P$. We consider a dense model with
$\phi=(1,\ldots,1)/\sqrt{P}$, a randomly chosen $\beta$ vector
$\beta_{j}\stackrel{{\scriptstyle\textnormal{iid}}}{{\sim}}\operatorname{Normal}(0,P^{-1})$,
and $\sigma^{2}_{a}=\sigma^{2}_{y}=1$. The methods differ in the treatment
effect size $\gamma$ and the degree to which the coefficients are shifted in
the direction of $\phi$. We considered three simulation settings.

Random

Both $\gamma$ and $\omega$ are $\operatorname{Normal}(0,1)$ random variables
and differ for each replication of the experiment.

Fixed

We set $\gamma=2$ and $\omega=-\gamma/4$ so that $\beta$ is shifted in the
direction of $\phi$, but not by the amount implied by the debiased approach.

Debiased

We set $\gamma=2$ and $\omega=-2$ so that $\beta$ is shifted in the direction
of $\phi$ by exactly the amount implied by the debiased approach.

Naive

We set $\gamma=2$ and $\omega=0$ so that the model corresponds precisely to
the naive ridge model.

The simulation was replicated 200 times for each setting. We evaluated each
procedure according to the following criteria. Coverage: The proportion of
nominal 95% credible intervals which capture the true value of $\gamma$.
Width: The average width of the nominal 95% credible interval. Avg SE: The
average estimated standard error from the model, i.e., the the posterior
standard deviation of $\gamma$ averaged over all replications. RMSE: The root
mean squared error in estimating $\gamma$ with the Bayes estimator
$\widehat{\gamma}$.

Setting | Method | Coverage | Width | Avg SE | RMSE
---|---|---|---|---|---
$\omega\sim\operatorname{Normal}(0,1)$ | Direct | 0.94 | 0.39 | 0.10 | 0.10
Debiased | 0.94 | 0.51 | 0.13 | 0.12
Naive | 0.19 | 0.31 | 0.08 | 0.49
$\omega=\gamma/4$ | Direct | 0.95 | 0.39 | 0.10 | 0.11
Debiased | 0.94 | 0.54 | 0.13 | 0.14
Naive | 0.12 | 0.29 | 0.07 | 0.25
$\omega=-\gamma$ | Direct | 0.96 | 0.39 | 0.10 | 0.11
Debiased | 0.95 | 0.39 | 0.10 | 0.11
Naive | 0.00 | 0.40 | 0.10 | 0.96
Naive ($\omega=0$) | Direct | 0.95 | 0.39 | 0.10 | 0.11
Debiased | 0.93 | 0.63 | 0.16 | 0.16
Naive | 0.90 | 0.28 | 0.07 | 0.08
Table 1: Comparison of different approaches for estimating $\gamma$ under
differing levels of selection bias. Direct denotes the approach which sets
$\omega\sim\operatorname{Flat}$, Debiased denotes the approach of Hahn et al.,
(2018).

Results are compiled in Table 1. The direct and debiased approaches always
attain the nominal coverage level, while the naive approach does not come
close when the selection bias is non-negligible. We also see that the debiased
approach will generally require substantially larger intervals than the direct
approach to cover at the appropriate rate; the only exception is when
$\omega=-\gamma$, which is to be expected as this setting agrees exactly with
the debiased prior. When the naive ridge model actually holds (i.e.,
$\omega=0$) we see that the naive ridge model unsurprisingly performs
substantially better, and is the best in terms of RMSE, with the direct prior
still outperforming the debiased prior.

###  4.2 Variable Sharing for Spike-and-Slab Priors

For the variable selection prior it remains valid to include the “clever
covariate” $X_{i}^{\top}\widehat{\phi}$ in the model to correct for dogmatism.
However, there are other approaches one can take which make specific use of
the variable selection aspect of the model. One possibility is to use _shared
variable selection_ for the two models; in particular, we want to ensure any
variable appearing in the selection model should also appear in the outcome
model. To implement this, we might set
$\phi_{j}\stackrel{{\scriptstyle\textnormal{iid}}}{{\sim}}(1-p_{\phi})\,\delta_{0}+p_{\phi}\,\operatorname{Normal}(0,\tau^{2}_{\phi})$
and conditionally set
$\beta_{j}\stackrel{{\scriptstyle\textnormal{indep}}}{{\sim}}\\{1-p_{\beta}(\phi_{j})\\}\,\delta_{0}+p_{\beta}(\phi_{j})\,\operatorname{Normal}(0,\tau^{2}_{\beta})$.
Setting $p_{\beta}(\phi_{j})=1$ if $\phi_{j}\neq 0$ guarantees that
$\beta_{j}$ will be included whenever $\phi_{j}$ is included.

We conduct a small simulation experiment to justify our claim that shared
variable selection is an effective strategy for combating dogmatism. For our
ground truth we consider $N=P=200$, $\sigma^{2}_{a}=\sigma^{2}_{y}=1$, and
$\gamma=1$. We then sample $\phi_{j}$ from the spike-and-slab prior with
$\tau^{2}_{\phi}=1$ and $p_{\phi}=5/200$. We then consider four different
schemes for sampling $\beta$. Naive: we sample $\beta_{j}$ from the spike-and-
slab prior with $p_{\beta}=5/200$ and $\tau^{2}_{\beta}=1$; Shared: we sample
$\beta_{j}$ from the spike-and-slab prior with $p_{\beta}=5/200$ if
$\phi_{j}=0$ and $p_{\beta}=1$ if $\phi_{j}\neq 0$; Direct: we sample
$\beta_{j}$ according to the Naive prior and then add $-\phi_{j}$ to it; and
Both: we sample $\beta_{j}$ according to the Shared prior then we add
$-\phi_{j}$ to it.

We compare the following prior specifications.

Naive

$\beta_{j}$ has a spike-and-slab prior with $p_{\beta}\equiv 5/200$ and which
is independent of $\phi_{j}$.

Shared

We use a spike-and-slab Z-prior with $p_{\beta,j}=5/200$ if
$\Pr(\phi_{j}=0\mid\bm{A},\bm{X})<0.5$, and $p_{\beta,j}=1$ otherwise.

Direct

We use the Naive prior with the additional covariate
$\widehat{A}_{i}=X_{i}^{\top}\widehat{\phi}$ where
$\widehat{\phi}=\mathbb{E}(\phi\mid\bm{A})$. The variable $\widehat{A}_{i}$ is
included with probability $1$.

Setting | Method | Coverage | Width | Std. Err | RMSE | Bias
---|---|---|---|---|---|---
Shared | Direct | 0.95 | 0.29 | 0.07 | 0.07 | 0.00
Naive | 0.65 | 0.25 | 0.07 | 0.13 | -0.00
Shared | 0.94 | 0.28 | 0.07 | 0.07 | 0.00
Direct | Direct | 0.95 | 0.29 | 0.07 | 0.08 | -0.03
Naive | 0.71 | 0.30 | 0.08 | 0.34 | -0.17
Shared | 0.93 | 0.29 | 0.07 | 0.08 | -0.03
Naive | Direct | 0.96 | 0.28 | 0.07 | 0.07 | -0.00
Naive | 0.94 | 0.14 | 0.04 | 0.04 | -0.00
Shared | 0.97 | 0.28 | 0.07 | 0.07 | -0.00
Both | Direct | 0.94 | 0.29 | 0.07 | 0.07 | 0.00
Naive | 0.77 | 0.27 | 0.07 | 0.12 | -0.01
Shared | 0.94 | 0.28 | 0.07 | 0.07 | 0.00
Table 2: Results for the simulation experiment described in Section 4.2

Results for this simulation are given in Table 4.2. Each simulation setting
was replicated $200$ times. We see that the Direct and Shared methods perform
essentially the same despite correcting for dogmatism in different ways — both
methods have virtually identical coverage, root-mean-squared error, interval
lengths, and bias. The Naive approach, which just applies the spike-and-slab
prior for $\beta$ under IG.2, performs extremely poorly by contrast, unless
the data was generated under the Naive prior.

###  4.3 Semiparametric Regression with Clever Covariates

Mimicking our strategy in Section 4.1 we set
$\beta(x)=\beta^{\star}(x)+g\\{\phi(x)\\}$ for some choice of $g(\cdot)$, with
$\beta^{\star}(x)$ given (say) a Gaussian process prior independent of
$\phi(\cdot)$. The selection bias for this model is given by

$\displaystyle\Delta=\frac{\operatorname{Cov}_{\theta}\\{\beta^{\star}(X_{i}),\phi(X_{i})\\}}{\mathbb{E}_{\theta}\\{\phi(X_{i})\\}}+\frac{\operatorname{Cov}_{\theta}[g\\{\phi(X_{i})\\},\phi(X_{i})]}{\mathbb{E}_{\theta}\\{\phi(X_{i})\\}}\approx\frac{\operatorname{Cov}_{\theta}[g\\{\phi(X_{i})\\},\phi(X_{i})]}{\mathbb{E}_{\theta}\\{\phi(X_{i})\\}}$

by the orthogonality principle. Even if we model $g(\cdot)$ nonparametrically
using (say) a Gaussian process the orthogonality principle will not kick in
because $g(\phi)$ is a function on $\mathbb{R}$ rather than $\mathbb{R}^{P}$.

There are many different choices one can make for the function $g(\cdot)$. If
we are concerned strictly with obtaining good Frequentist properties, an
appropriate choice is to take $g(\phi)=\phi^{-1}$; when $\phi$ is known, this
guarantees $\sqrt{n}$-consistency. Alternatively, we can set
$g\sim\operatorname{GP}(0,\kappa_{g})$ with the covariance function
$\kappa(\phi,\phi^{\prime})=\tau^{2}_{g}\exp\\{-(\phi-\phi^{\prime})^{2}/(2s^{2}_{g})\\}$.
This choice of covariance function was noted by Ren et al., (2021) to induce
matching on the propensity score: individuals with similar propensity scores
have their values of $g(\phi)$ shrunk together. The penalized-spline-of-
propensity approach of Zhou et al., (2019) is similar, except that $g(\phi)$
is chosen to be a spline instead of a Gaussian process. Another approach to
incorporating $\phi(x)$ is to plug it in as a regular covariate, i.e., we
replace $\beta(x)$ with $\beta\\{x,\phi(x)\\}$.

#### Simulation Experiment

We consider the simulation setting of Hahn et al., (2020, Section 6.1) to
evaluate several different approaches to correcting a Gaussian process prior
for dogmatism using the usual evaluation criteria (RMSE, interval width, and
bias, and coverage). We consider the model
$Y_{i}(a)=\mu(X_{i})+a\,\tau(X_{i})+\epsilon_{i},\epsilon_{i}\sim\operatorname{Normal}(0,1)$
with the $X_{ij}$’s being iid $\operatorname{Normal}(0,1)$ random variables
with the exception of $X_{i2}$ which is $\operatorname{Bernoulli}(1/2)$ and
$X_{i4}$ which is uniform on $\\{1,2,3\\}$. We let $\mu(x)$ and $\tau(x)$ be
given

$\displaystyle\tau(x)=\begin{cases}3\quad&\text{homogeneous},\\\\[-10.00002pt]
1+2\,x_{2}\,x_{5}&\text{heterogeneous},\end{cases}\quad\text{and}\quad\mu(x)=\begin{cases}1+g(x_{4})+x_{1}\,x_{3}\quad&\text{linear},\\\\[-10.00002pt]
-6+g(x_{4})+6\,|x_{3}-1|&\text{nonlinear},\end{cases}$ (4)

where $g(1)=2,g(2)=-1$, and $g(3)=-4$. We then set
$A_{i}\sim\operatorname{Bernoulli}\\{\phi(X_{i})\\}$ with
$\phi(x)=0.8\,\Phi\\{3\,\mu(x)/s-0.5\,x_{1}\\}+0.1$, where $s$ is the
empirical standard deviation of the $\mu(X_{i})$’s. In total, we consider 16
possible simulation settings, corresponding to a factorial design with
$N\in\\{250,500\\}$, $P\in\\{5,20\\}$, with four combinations of
linear/nonlinear and homogeneous/heterogeneous.

We consider modeling $\mathbb{E}\\{Y_{i}(a)\mid X_{i}=x\\}=\beta(a,x)$ using a
Gaussian process $\beta\sim\operatorname{GP}(0,\kappa)$ where the kernel
function $\kappa\big{(}(a,x),(a^{\prime},x^{\prime})\big{)}$ is given by the
following choices.

Naive

A kernel which makes no correction for the propensity score:
$\kappa\big{(}(a,x),(a^{\prime},x^{\prime})\big{)}=100(1+a\,a^{\prime})+\lambda\,\exp\\{-b\|(a,x)-(a^{\prime},x^{\prime})\|^{2}_{2}\\}$.

IPW-GP

A kernel which incorporates the inverse propensity score as a “clever
covariate” which enters the model linearly:
$\kappa\big{(}(a,x),(a^{\prime},x^{\prime})\big{)}=100(1+a\,a^{\prime}+w\,w^{\prime}+z\,z^{\prime})+\lambda\,\exp\\{-b\|(a,x)-(a^{\prime},x^{\prime})\|^{2}_{2}\\}$
where $w=a/\phi(x)$ and $z=(1-a)/(1-\phi(x))$.

Spline-of-propensity-GP

A kernel which incorporates the propensity score using a spline basis function
expansion:
$\kappa\big{(}(a,x),(a^{\prime},x^{\prime})\big{)}=100(1+a\,a^{\prime}+\sum_{k}\psi_{k}\,\psi_{k}^{\prime})+\lambda\exp\\{-b\|(a,x)-(a^{\prime},x^{\prime})\|^{2}_{2}\\}$
where $\psi_{k}=\psi_{k}(x)$, $\psi^{\prime}_{k}=\psi(x^{\prime})$, and
$\\{\psi_{1},\ldots,\psi_{K}\\}$ are natural cubic spline basis functions
using $10$ knots.

Spline-of-propensity

Same as Spline-of-propensity-GP but without the Gaussian kernel.

In order to to separate the issue of accurately estimating the propensity
scores from the benefit of using them, we assume that $\phi(x)$ is known
a-priori. For all methods, the kernel hyperparameters and the standard
deviation $\epsilon$ are estimated via empirical Bayes, i.e., by maximizing
the marginal likelihood of the data.

Our main goals are to (i) determine the extent to which the Naive kernel
suffers due to dogmatism, (ii) determine which of the IPW or spline approaches
perform better in this case, and (iii) determine whether the propensity score
alone is sufficient to produce a good estimator. A subset of the results
corresponding to the nonlinear heterogeneous setting with $N=250$ are given in
Figure 4, with the remaining results deferred to the Supplementary Material.
Summarizing these results, we find (i) that the Naive kernel performs well
when $P=5$ where dogmatism is mild, but breaks down completely when $P=20$;
(ii) that the IPW-GP and spline-of-propensity-GP approaches perform comparably
in terms of coverage, but that the spline-of-propensity-GP generally produces
smaller standard errors and RMSEs, suggesting that the spline-of-propensity
approach is more stable while accomplishing the same goals as IPW methods; and
(iii) that the spline-of-propensity-GP produces smaller standard errors and
RMSEs than the spline-of-propensity approach, indicating that there is some
benefit to going beyond simply adjusting for the propensity score.

Figure 4: Results for the simulation study of Section 4.3 in the nonlinear
heterogeneous setting. Naive denotes the naive approach, IPW denotes the IPW-
GP approach, SOP denotes the spline-of-propensity approach, and SOP-GP denotes
the spline-of-propensity-GP approach.

##  5 Factors Mitigating Dogmatism

In demonstrating the issue of dogmatism we made use of the orthogonality
principle, which we noted only holds when “other structures” are not present.
For example, we introduced additional structure into the model by making
$\beta$ and $\phi$ dependent in Section 4, giving us direct control of the
prior on $\Delta(a)$.

Another possible source of structure is dependence structure in $X_{i}$. The
benefit of this is evident in Proposition 2 and Theorem 1, where the spectral
distribution of $\Sigma$ figures prominently. We examine the role of
dependence structure in the ridge and semiparametric regression problems.

###  5.1 Dependence Structure and Ridge Regression

To understand the role of $\Sigma$ in ridge regression, we consider a _latent
factor model_ which takes $X_{i}=\Lambda\eta_{i}+\sigma_{x}\,\nu_{i}$ where
$\Lambda\in\mathbb{R}^{P\times L}$ is a matrix of factor loadings and
$\eta_{i}\in\mathbb{R}^{L}$ is an $L$-dimensional vector of latent factors for
observation $i$. If $\sigma_{x}=0$ in this model then $X_{i}$ is restricted to
be in the $L$-dimensional subspace $\operatorname{span}(\Lambda)$; similarly,
if $\sigma_{x}$ is small, then $X_{i}$ lies very close to
$\operatorname{span}(\Lambda)$.

We first examine the induced prior on the selection bias parameter for such a
$\Sigma$. Assuming $\eta_{i}\sim\operatorname{Normal}(0,\mathrm{I})$ and
$\nu\sim\operatorname{Normal}(0,\mathrm{I})$, we have
$\operatorname{Var}(X_{i})\equiv\Sigma=\Lambda\Lambda^{\top}+\sigma_{x}^{2}\mathrm{I}$.
Letting $\kappa_{1},\ldots,\kappa_{L}$ denote the $L$ non-zero eigenvalues of
$\Lambda\Lambda^{\top}$, Proposition 1 gives

$\displaystyle\Delta(a)$
$\displaystyle=a\frac{\sum_{j=1}^{L}(\kappa_{j}+\sigma^{2}_{x})\,W_{j}\,Z_{j}}{\sigma^{2}_{a}+\sum_{j=1}^{L}(\kappa_{j}+\sigma^{2}_{x})\,Z_{j}^{2}}+a\frac{\sum_{j=L+1}^{P}\sigma^{2}_{x}\,W_{j}\,Z_{j}}{\sigma^{2}_{a}+\sum_{j=L+1}^{P}\sigma_{x}^{2}\,Z_{j}^{2}}\approx
a\frac{\sum_{j=1}^{L}\kappa_{j}\,W_{j}\,Z_{j}}{\sigma^{2}_{a}+\sum_{j=1}^{L}\kappa_{j}\,Z_{j}^{2}}$

where
$W_{j}\stackrel{{\scriptstyle\textnormal{iid}}}{{\sim}}\operatorname{Normal}(0,\tau^{2}_{\beta})$
and
$Z_{j}\stackrel{{\scriptstyle\textnormal{iid}}}{{\sim}}\operatorname{Normal}(0,\tau_{\phi}^{2})$.
The approximation holds when $\sigma_{x}$ is near zero so that $\Sigma$ is
approximately low-rank. Because the left-hand-side depends only on $L$ rather
than $P$, we expect $\Delta(a)$ to be roughly of order $L^{-1/2}$ rather than
$P^{-1/2}$ for the ridge regression prior. Hence, even if $P\gg N$, we may
still avoid dogmatism if $L\ll N$.

Revisiting Theorem 1 we can characterize the effect of $\Sigma$ on the bias.
Recall that Theorem 1 relates the bias of the naive ridge estimator to the
empirical spectral distribution $F(dx)$ of $\bm{X}\bm{X}^{\top}/N$ through the
Stieltjes transform $v(-\lambda)=\int_{0}^{\infty}\frac{F(dx)}{x+\lambda}$.
Smaller values of $\lambda$ such that $v(-\lambda)\approx\lambda^{-1}$ result
in small bias, which occurs when $F(dx)$ places substantial mass near $0$. In
Figure 5 we plot both $v(-\lambda)$ and the bias for the latent class model
with $L=5$ and the entries of $\Lambda$ being iid
$\operatorname{Normal}(0,1)$. We see that as $\sigma_{x}$ decreases we have
substantially less bias. The ridge regression estimator is able to leverage
the fact that $X_{i}$ is approximately in $\operatorname{span}(\Lambda)$
without this being explicitly encoded into the model. In the extreme case
where $\sigma_{x}=0$ this is easy to see, as we can write
$\bm{X}=\bm{E}\bm{D}^{1/2}\Gamma^{\top}$ where $\bm{E}\in\mathbb{R}^{N\times
L}$ has iid $\operatorname{Normal}(0,1)$ entries,
$\Gamma\in\mathbb{R}^{P\times L}$ is semi-orthogonal,
$\bm{D}\in\mathbb{R}^{L\times L}$ is diagonal, and
$\Sigma=\Gamma\bm{D}\Gamma^{\top}$. We can then rewrite
$\bm{X}\beta=\bm{E}\bm{D}^{1/2}\Gamma^{\top}\beta=\bm{E}\zeta$ where
$\zeta\sim\operatorname{Normal}(0,\bm{D}\,\tau^{2}_{\beta})$. Hence performing
ridge regression on the $N\times P$-dimensional $\bm{X}$ is exactly equivalent
to performing ridge regression on the lower-dimensional matrix $\bm{E}$, with
$\zeta_{\ell}$ having variance proportional to the $\ell^{\text{th}}$ non-zero
eigenvalue of $\Sigma$.

Figure 5: Top: $v(-\lambda)$ for the latent factor model for different
$(r,\sigma_{x})$ with a dashed line giving the ideal
$v(-\lambda)=\lambda^{-1}$. Bottom: the associated bias of the ridge
regression estimator.

Figure 6 shows the root mean squared error of the direct approach we proposed
in Section 4.1 and the naive ridge regression estimator. To generate this
figure, we applied the two approaches under the following conditions:
$\sigma^{2}_{a}=\sigma^{2}_{y}=1$; $L=5$; $N=P=200$;
$\Lambda_{p\ell}\stackrel{{\scriptstyle\textnormal{iid}}}{{\sim}}\operatorname{Normal}(0,1)$;
$\gamma=1$; and both $\beta$ and $\phi$ chosen so that
$X_{i}^{\top}\beta=X_{i}^{\top}\phi=\sum_{\ell=1}^{L}\mathbb{E}(\eta_{i\ell}\mid
X_{i})$.

Figure 6:  Root mean squared error in estimating $\gamma$ for direct and naive
methods as a function of $\sigma_{x}$.

###  5.2 Dependence Structure and Semiparametric Regression

We present evidence that the results for ridge regression is part of a much
more general phenomenon which also holds in nonparametric problems: namely,
that low-dimensional structure can shield us from the effects of prior
dogmatism. Rather than assuming that $X_{i}$ is near a hyperplane, we now
assume that $X_{i}$ is concentrated near a smooth manifold of intrinsic
dimension $L$ in $\mathbb{R}^{P}$. Specifically, we take
$\widetilde{X}_{i}=\Lambda(\eta_{i})+\sigma_{x}\,\epsilon_{i},\epsilon_{i}\sim\operatorname{Normal}(0,1)$
where $\Lambda:\mathbb{R}^{L}\to\mathbb{R}^{P}$ is nonlinear; we then scale
$\widetilde{X}_{ij}$ by its standard deviation to get $X_{ij}$, so that
$\sigma_{x}$ indexes how close $X_{i}$ is to
$\mathscr{M}=\\{x:x=\Lambda(\eta),\eta\in\mathbb{R}\\}$. We randomly generated
$\Lambda(\eta)=\big{(}\Lambda_{1}(\eta),\ldots,\Lambda_{P}(\eta)\big{)}^{\top}$
by generating $P$ independent Gaussian processes using the kernel function
$\rho(\eta,\eta^{\prime})=\exp\\{-\|\eta-\eta^{\prime}\|^{2}_{2}\\}$.

After generating the covariates $X_{i}$, we consider a continuous exposure
$A_{i}=r_{a}(X_{i})+\nu_{i}$ and a continuous outcome
$Y_{i}(a)=r_{y}(X_{i})+a\,\gamma+\epsilon_{i}(a)$ with
$\epsilon_{i}(a),\nu_{i}\stackrel{{\scriptstyle\textnormal{iid}}}{{\sim}}\operatorname{Normal}(0,1)$
where $r_{y}(x)=r^{\star}_{y}(x)+r_{a}(x)$. We then generate $r^{\star}_{y}$
and $r_{a}$ as independent Gaussian processes with kernel
$\rho(x,x^{\prime})=\exp\\{-\|x-x^{\prime}\|^{2}_{2}\\}$. Our parameter of
interest is $\gamma$, which represents the causal effect of the exposure on
the outcome. We consider two priors.

Naive

We impose IG.2, but otherwise use the “true” prior for $r^{\star}(x)$ using
the kernel $2\rho(x,x^{\prime})$. We specify a
$\operatorname{Normal}(0,10^{2})$ prior for $\gamma$.

Direct

We use the model
$Y_{i}(a)=r_{y}(X_{i})+\omega\,\widehat{r}_{a}(X_{i})+\gamma\,a$ where
$\widehat{r}_{a}(x)$ is a pilot estimate of $r_{a}(x)$ obtained from fitting a
Gaussian process to the relationship $A_{i}=r_{a}(X_{i})+\nu_{i}$. We specify
a $\operatorname{Normal}(0,10^{2})$ prior for both $\gamma$ and $\omega$.

For the ground truth, we set $\gamma=1$, $L=1$, and consider
$P\in\\{10,200\\}$, $N=300$, and $\sigma_{x}=2^{-j}$ where the $j$’s are
evenly spaced between $-7$ and $-2$. For each $\sigma_{x}$ and $P$ we
generated $200$ simulated datasets and applied the Direct and Naive methods to
estimate $\gamma$ and construct a 95% credible interval.

Figure 7: Results for the semiparametric regression on a manifold problem of
Section 5. Bias denotes the average bias of $\widehat{\gamma}$, coverage
denotes the coverage of nominal 95% intervals, RMSE denotes the root-mean-
squared-error in estimating $\gamma$, and SE denotes the average posterior
standard deviation of $\gamma$.

Results are summarized in Figure 7, which displays the bias, coverage, RMSE,
and average standard error. For $P=10$, we see similar behavior as we did with
the ridge regression problem: the Direct approach performs uniformly well and
the naive approach performs much better as $\sigma_{x}$ is decreased. While
the naive approach never does catch up to the direct approach, we do see that
its deficiencies are attenuated as the $X_{i}$’s are generated closer and
closer to $\mathscr{M}$. For $P=200$, while the naive approach does perform
better as $\sigma_{x}$ is decreased, the problem appears to be too difficult
for methods which do not explicitly account for dogmatism. The behavior of the
direct approach is now very interesting, however. We note a sharp phase
transition around $\sigma_{x}=0.07$ where the problem essentially goes from
infeasible to feasible — the bias, RMSE, and standard error all decrease
dramatically at this point. We also see that the direct approach is much more
honest in terms of its uncertainty: when the problem is infeasible, the model
correctly gives a large posterior standard error, whereas the naive model is
always overconfident about its predictions.

##  6 Discussion

The main concrete recommendation we make in this article is that, for both
causal inference and missing data problems, Bayesian ignorability (and in
particular IG.2) corresponds to an informative prior on the degree of
selection bias. This causes selection bias parameters to be regularized
towards $0$, introducing substantial bias in high-dimensional or nonparametric
problems and should not be imposed in most situations. Instead, Bayesians
should reject IG.2 by default in favor of a prior which allows for more direct
control over the selection bias, and we have illustrated how to do this in
several problems of interest. Of secondary interest, we have noted that
certain features of the design can mitigate prior dogmatism about the
selection bias, and showed that both ridge regression priors and Gaussian
process priors possess some degree of adaptivity towards low-dimensional
structures in $X_{i}$. But this does not change our general recommendation, as
we consistently have observed improved performance of priors which reject IG.2
even when such low-dimensional structures exist.

Dogmatism about other features of the model may also have deleterious effects
on our inferences. In future work, we will extend our results to other
problems in causal inference. Two of potential interest are estimation of the
conditional average treatment effect (CATE) in observational studies and
estimation of the natural direct and indirect effects in mediation analysis.
In the latter setting, one must control for two different selection
mechanisms: the effect of the confounders both on the treatment received and
on the mediating variable.

While we have presented a number of corrections for dogmatism, we have not
presented any coherent framework for _deriving_ corrections. This presents an
important question: are there any objective Bayes principles which
automatically lead to priors which adequately account for dogmatism? Certain
strategies, such as using Jeffreys priors, cannot work because they usually
_imply_ that IG.2 holds. By contrast, other objective principles which are not
parameterization invariant and do not necessarily imply IG.2, such as priors
constructed from decision theoretic principles, entropy maximization, and
reference priors have some chance of working (see Kass and Wasserman,, 1996,
for a review). With the exception of entropy-maximizing priors, the
computational difficulty of implementing these priors makes numeric
experimentation difficult. Interestingly, entropy maximization with respect to
the distribution of the observed data can be used to generate models which
possess very strong Frequentist properties, but these models are (in our
opinion) lacking a satisfying justification.

## References

  * Chipman et al., (2010) Chipman, H. A., George, E. I., and McCulloch, R. E. (2010).  BART: Bayesian additive regression trees.  The Annals of Applied Statistics, 4(1):266–298.
  * Couillet and Debbah, (2011) Couillet, R. and Debbah, M. (2011).  Random Matrix Methods for Wireless Communications.  Cambridge University Press.
  * Dicker, (2016) Dicker, L. H. (2016).  Ridge regression and asymptotic minimax estimation over spheres of growing dimension.  Bernoulli, 22(1):1–37.
  * Dobriban and Wager, (2018) Dobriban, E. and Wager, S. (2018).  High-dimensional asymptotics of prediction: Ridge regression and classification.  The Annals of Statistics, 46(1):247–279.
  * Hahn et al., (2018) Hahn, P. R., Carvalho, C. M., Puelz, D., and He, J. (2018).  Regularization and confounding in linear regression for treatment effect estimation.  Bayesian Analysis, 13(1):163–182.
  * Hahn et al., (2020) Hahn, P. R., Murray, J. S., and Carvalho, C. M. (2020).  Bayesian regression tree models for causal inference: Regularization, confounding, and heterogeneous effects (with discussion).  Bayesian Analysis, 15(3):965–1056.
  * Hill, (2011) Hill, J. L. (2011).  Bayesian nonparametric modeling for causal inference.  Journal of Computational and Graphical Statistics, 20(1):217–240.
  * Hirano et al., (2003) Hirano, K., Imbens, G. W., and Ridder, G. (2003).  Efficient estimation of average treatment effects using the estimated propensity score.  Econometrica, 71(4):1161–1189.
  * Imai et al., (2010) Imai, K., Keele, L., and Tingley, D. (2010).  A general approach to causal mediation analysis.  Psychological Methods, 15(4):309.
  * Kass and Wasserman, (1996) Kass, R. E. and Wasserman, L. (1996).  The selection of prior distributions by formal rules.  Journal of the American statistical Association, 91(435):1343–1370.
  * National Research Council, (2010) National Research Council (2010).  The Prevention and Treatment of Missing Data in Clinical Trials.  The National Academies Press.
  * Rasmussen and Williams, (2006) Rasmussen, C. E. and Williams, C. K. I. (2006).  Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning).  MIT Press, Cambridge.
  * Ray and van der Vaart, (2020) Ray, K. and van der Vaart, A. (2020).  Semiparametric Bayesian causal inference.  The Annals of Statistics, 48(5):2999 – 3020.
  * Ren et al., (2021) Ren, B., Wu, X., Braun, D., Pillai, N., and Dominici, F. (2021).  Bayesian modeling for exposure response curve via Gaussian processes: Causal effects of exposure to air pollution on health outcomes.  arXiv preprint arXiv:2105.03454.
  * Robert, (2007) Robert, C. (2007).  The Bayesian choice: from decision-theoretic foundations to computational implementation.  Springer Science & Business Media.
  * Robins and Ritov, (1997) Robins, J. M. and Ritov, Y. (1997).  Toward a curse of dimensionality appropriate (CODA) asymptotic theory for semi-parametric models.  Statistics in medicine, 16(3):285–319.
  * Rosenbaum and Rubin, (1983) Rosenbaum, P. R. and Rubin, D. B. (1983).  The central role of the propensity score in observational studies for causal effects.  Biometrika, 70(1):41–55.
  * Rubin, (1976) Rubin, D. B. (1976).  Inference and missing data.  Biometrika, 63:581–592.
  * Rubin, (1985) Rubin, D. B. (1985).  The use of propensity scores in applied Bayesian inference.  Bayesian Statistics, 2:463–472.
  * Rubin, (2005) Rubin, D. B. (2005).  Causal inference using potential outcomes.  Journal of the American Statistical Association, 100(469).
  * Seaman et al., (2013) Seaman, S., Galati, J., Jackson, D., and Carlin, J. (2013).  What is meant by “missing at random”?  Statistical Science, 28(2):257–268.
  * van der Laan and Rose, (2011) van der Laan, M. J. and Rose, S. (2011).  Targeted Learning: Causal Inference for Observational and Experimental Data.  Springer Science & Business Media.
  * Zhou et al., (2019) Zhou, T., Elliott, M. R., and Little, R. (2019).  Penalized spline of propensity methods for treatment comparison.  Journal of the American Statistical Association, 114(525):1–19.
  * Zigler et al., (2013) Zigler, C. M., Watts, K., Yeh, R. W., Wang, Y., Coull, B. A., and Dominici, F. (2013).  Model feedback in Bayesian propensity score estimation.  Biometrics, 69(1):263–273.

See pages - of 2021-Informative-Supplement.pdf