# Calibrated inference: statistical inference that accounts for both sampling
uncertainty and distributional uncertainty

Yujin Jeong and Dominik Rothenhäusler

###### Abstract

How can we draw trustworthy scientific conclusions? One criterion is that a
study can be replicated by independent teams. While replication is critically
important, it is arguably insufficient. If a study is biased for some reason
and other studies recapitulate the approach then findings might be
consistently incorrect. It has been argued that trustworthy scientific
conclusions require disparate sources of evidence. However, different methods
might have shared biases, making it difficult to judge the trustworthiness of
a result. We formalize this issue by introducing a “distributional uncertainty
model”, which captures biases in the data collection process. Distributional
uncertainty is related to other concepts in statistics, ranging from
correlated data to selection bias and confounding. We show that a stability
analysis on a single data set allows to construct confidence intervals that
account for both sampling uncertainty and distributional uncertainty.

##  1 Introduction

Statistical inferences can be fragile. If we compare two analyses conducted by
different statisticians on different data sets, variation can be due to
sampling, due to distribution shift, or due to a change in methodology. These
issues raise a fundamental question: How can we draw trustworthy scientific
conclusions? A common recommendation is to have independent teams attempt to
replicate the findings of others. While replication is critically important,
it is arguably insufficient. If a study suffers from biases for some reason
and replication studies emulate the study, the findings will be consistently
incorrect [Munafò and Smith, 2018]. To solve this issue, researchers have
advocated investigating independent lines of evidence [Denzin, 1970, Freedman,
1991, Rosenbaum, 2010, Munafò and Smith, 2018]. Ideally, these different lines
of evidence are susceptible to different biases. Intuitively, if results agree
across different methodologies, then a statistical finding is less likely to
be an artifact. However, it might be expensive and impractical to ask several
researchers to run studies independently.

Can we emulate this strategy on a single data set? In fact, stability analyses
have been advocated by many researchers. To be more precise, it has been
recommended to evaluate several reasonable modelling choices for one single
data set [Leamer, 1983, Rosenbaum, 2010, Patel et al., 2015, Steegen et al.,
2016, Yu and Kumbier, 2020]. Common practice includes running differently
specified regressions or considering the perturbations induced by various
forms of data pre-processing. If the estimator-to-estimator variability is
high, then the analyst has reason to distrust the estimates.

This warrants an investigation of the theoretical properties of this practice.
What mathematical problem do we solve with this approach? What distinguishes a
good stability analysis from a bad one? If differently specified regressions
return similar estimates, what criterion tells us whether we should be
concerned or not? Often the decision what qualifies as a “stable result” is
left up to the individual judgements of the analyst. One could analyze the
estimators with a random effects model, but since the estimators are computed
on the same data set, they might share biases. More concretely speaking, since
the structure of the biases is generally unknown, it is not clear how to
define a random effects model that captures the correlation structure of the
estimators’ biases.

In order to address these questions in a formal framework, we take a
distributional perspective. Outliers, sampling bias, confounding, measurement
errors, and unknown correlations between units are omnipresent in data
analysis and make it challenging to gather trustworthy results. We consider a
setting, where the data is drawn from a “perturbed” or “contaminated”
distribution, while the goal is to infer some properties of the uncontaminated
distribution. The classical robust statistics literature [Huber, 1981]
addresses distributional perturbations by investigating the worst-case
behavior of a statistical functional over a fixed neighborhood of the model.
More recently, distributional uncertainty sets based on $f$-divergence have
been linked to distributionally robust optimization [Ben-Tal et al., 2013,
Duchi et al., 2021]. In such models, it is challenging to choose the
appropriate set of distributions, since the size of the perturbation is
generally unknown.

The distributional perturbations we want to model are due to small errors in
the data collection process that potentially induce selection bias,
confounding, batch effects, or measurement errors. We do not believe that
these issues are adversarial in nature. Thus, instead of considering worst-
case distributional perturbations, we define a family of random symmetric
perturbations. We show that modelling distributional perturbations as random
and symmetric has an intriguing consequence: using a stability analysis, it is
possible to estimate the strength of the distributional perturbation. Based on
an estimate of the distributional perturbation strength, we propose confidence
intervals that capture both sampling uncertainty and distributional
uncertainty.

###  1.1 Related work

Considerations of model stability have emerged in Bayesian statistics [Box,
1980, Skene et al., 1986], causal inference [Leamer, 1983, LaLonde, 1986,
Rosenbaum, 1987, Imbens and Rubin, 2015] and in discussions about the data
science lifecycle [Yu, 2013, Steegen et al., 2016, Yu and Kumbier, 2020].
Using different estimation strategies is commonly recommended to corroborate a
causal hypothesis [Freedman, 1991, Rosenbaum, 2010, Karmakar et al., 2019]. In
particular, to evaluate omitted variable bias, it is a common recommendation
to consider the between-estimator variation of several adjusted regressions
[Oster, 2019]. Sensitivity analysis bounds the influence of confounders that
have been omitted in a regression or matching procedure and has played an
influential role in increasing trustworthiness of causal inference from
observational data [Cornfield et al., 1959, Rosenbaum and Rubin, 1983,
VanderWeele and Ding, 2017]. It has been argued that causal mechanisms are
expected to lead to stable associations across settings, if the same mechanism
is shared across settings. Based on this observation, stability principles
have been employed to discover causal relationships based on heterogeneous
data sets [Peters et al., 2016, Rothenhäusler et al., 2015, Bühlmann, 2020,
Pfister et al., 2021]. Stability principles are heavily used in machine
learning, often with the goal of variance reduction. For example, some tree-
based methods employ feature bagging, which can be seen as averaging over
differently specified prediction models [Breiman, 1996, 2001]. Dropout in
neural networks is another form of algorithm perturbation [Srivastava et al.,
2014]. Distributional uncertainty sets based on $f$-divergences have been
linked to distributionally robust optimization [Ben-Tal et al., 2013, Duchi et
al., 2021]. In the context of prediction under distribution shift, stability
or invariance principles have been employed to learn prediction mechanisms
that generalize to new settings [Schölkopf et al., 2012, Zhang et al., 2013,
Rojas-Carulla et al., 2018, Heinze-Deml and Meinshausen, 2021, Rothenhäusler
et al., 2021]. Quasi-likelihoods [Wedderburn, 1974] are a way to allow greater
variability in the data then what is expected from the model. However,
uncertainty quantification in quasi-likelihoods still only deals with sampling
uncertainty, while we aim to quantify uncertainty due to both sampling and
distributional uncertainty. As we show in the appendix, the considered
perturbation model implies equicorrelated data marginally. Thus, there is a
connection between our perturbation model and the problem of inference in
network data [Lee and Ogburn, 2021]. The main difference between our approach
and other approaches is that we do not model the network structure itself;
instead dependencies arise through a symmetric perturbation scheme.

###  1.2 Outline of the paper

In Section 1.3, we will quickly review standard practice for forming
confidence intervals. In Section 2, we introduce the setting of the paper and
discuss why standard statistical practice does not account for all types of
uncertainty in this setting. The setting of our paper arises under a
distributional perturbation model described in Section 2 and sampling
procedures described in the Appendix. We then turn to statistical inference.
In Section 3, we discuss how to form confidence intervals in our setting. This
completes the picture from a statistical viewpoint. In Section 4, we evaluate
the performance of the proposed procedure on a simulated example from causal
inference. In Section 5 we demonstrate that the proposed procedure can
increase the stability of decision-making based on real-world data. We
conclude in Section 6.

###  1.3 Standard approach

Let us consider estimation of the mean $\theta^{0}=\mathbb{E}[D]$ of a square-
integrable real-valued random variable $D\in\mathcal{D}$, $D\sim\mathbb{P}$.
Assume that we are given data
$(D_{i})_{i=1,\ldots,n}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\mathbb{P}$
with $\text{Var}(D_{i})=\sigma^{2}\in(0,\infty)$. We can estimate $\sigma^{2}$
via $\hat{\sigma}^{2}=\frac{1}{n-1}\sum_{i=1}^{n}(D_{i}-\overline{D})^{2}$ to
form asymptotically valid confidence intervals, that means

$\mathrm{P}(\overline{D}-z_{1-\alpha/2}\hat{\sigma}/\sqrt{n}\leq\theta^{0}\leq\overline{D}+z_{1-\alpha/2}\hat{\sigma}/\sqrt{n})\rightarrow
1-\alpha,$ (1)

where $z_{1-\alpha/2}$ is the $1-\alpha/2$ quantile of a standard Gaussian
random variable. This practice is justified by the central limit theorem which
implies

$\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(D_{i}-\mathbb{E}[D])\stackrel{{\scriptstyle
d}}{{=}}\mathcal{N}(0,\text{Var}(D))+o_{p}(1).$ (2)

More generally, for some vector-valued data
$D_{i}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\mathbb{P}$ consider a
parametrized model $\\{p_{\theta},\theta\in\Omega\\}$ of positive probability
densities $p_{\theta}$ with respect to some $\sigma$-finite measure $\mu$.
Assume that the parameter space $\Omega$ is an open subset of
$\mathbb{R}^{d}$. We consider the maximum-likelihood estimator

$\hat{\theta}=\arg\max\sum_{i=1}^{n}\log p_{\theta}(D_{i})$

for some unknown target parameter
$\theta^{0}(\mathbb{P})=\arg\max\mathbb{E}[\log p_{\theta}(D)]$, where
$D\sim\mathbb{P}$. Under regularity assumptions [van der Vaart, 2000, Tsiatis,
2006], for $n\rightarrow\infty$,

$\sqrt{n}(\hat{\theta}-\theta^{0})=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}-\mathbb{E}[\partial_{\theta}^{2}\log
p_{\theta^{0}}(D)]^{-1}\partial_{\theta}\log
p_{\theta^{0}}(D_{i})+o_{p}(1)\stackrel{{\scriptstyle
d}}{{=}}\mathcal{N}(0,\Sigma)+o_{p}(1),$

where $\Sigma=\mathbb{E}[\partial_{\theta}^{2}\log
p_{\theta^{0}}(D)]^{-1}\text{Var}(\partial_{\theta}\log
p_{\theta^{0}}(D))\mathbb{E}[\partial_{\theta}^{2}\log
p_{\theta^{0}}(D)]^{-1}$. Thus, based on a consistent estimator
$\hat{\Sigma}\rightarrow\Sigma$, one can form asymptotically valid confidence
intervals via

$\hat{\theta}_{k}\pm z_{1-\alpha/2}\frac{\sqrt{\hat{\Sigma}_{kk}}}{\sqrt{n}}.$
(3)

A similar approach can be used to construct asymptotically valid confidence
intervals for $M$-estimators. In the following, we discuss situations in which
this approach does not have the desired coverage.

##  2 Distributional uncertainty

There are several reasons why the coverage in equation (1) might not hold. The
main focus of this paper will be violations of (2) due to what we call
distributional uncertainty. Due to batch effects, sampling bias, confounding,
or other distributional shifts, the data scientist might not draw a sample
from the target distribution $\mathbb{P}$ but from some
$\mathbb{P}^{\xi}\neq\mathbb{P}$. The data analyst may try to address this
source of bias by re-weighting, regression adjustment, random effect modeling,
or other statistical techniques. Our viewpoint is that when using such
techniques, it is likely that some residual error remains. Ideally, we would
like to construct confidence intervals that detect residual errors due to
distributional perturbations, and accounts for it, if necessary. We model the
variation due to distributional changes as random. This will allow us to
integrate both distributional uncertainty and sampling uncertainty in a
natural fashion.

Let us make this more concrete by returning to the example of estimating the
mean. Due to unknown sampling bias or other unknown distributional shifts, the
data scientist might not draw a sample from the target distribution
$\mathbb{P}$. Instead, the data $(D_{i})_{i=1,\ldots,n}$ might be drawn i.i.d.
from some perturbed distribution $\mathbb{P}^{\xi}\neq\mathbb{P}$, where $\xi$
is a random variable and $\mathbb{P}^{\bullet}$ is a probability distribution
for each fixed $\bullet\in\text{range}(\xi)$. Then the error of the empirical
mean can be decomposed:

$\frac{1}{n}\sum_{i=1}^{n}D_{i}-\mathbb{E}[D]=\underbrace{\frac{1}{n}\sum_{i=1}^{n}D_{i}-\mathbb{E}^{\xi}[D]}_{\begin{subarray}{c}\text{variation
due to}\\\
\text{sampling}\end{subarray}}+\underbrace{\mathbb{E}^{\xi}[D]-\mathbb{E}[D]}_{\begin{subarray}{c}\text{variation
due to}\\\ \text{distributional perturbation}\end{subarray}}$ (4)

Equation (2) usually does not hold in this setting as distributional
perturbations induce additional variation.

Roughly speaking, in equation (4) we can distinguish three different regimes.
If variation due to sampling is much larger than the variation due to
distributional perturbations, intuitively one can simply ignore the variation
due to the second term and just compute confidence intervals as in the
previous section. Secondly, one can consider the setting where both sources of
variation are of the same order, that means the case where the distributional
variation is of the order $1/\sqrt{n}$. This is the setting we will focus on.
Thirdly, in practice one will encounter cases where the distributional
variation is of higher order than $1/\sqrt{n}$. In this case, one would have
to scale the confidence intervals in equation (3) by a factor that goes to
infinity for $n\rightarrow\infty$. We restrict ourselves to the second setting
mainly for expository reasons: It is desirable to have a methodology that
works in all three settings. To facilitate adoption by practitioners, ideally
such methodology should be a direct extension of current statistical practice.
Thus, it is natural to start with a regime where both types of uncertainty are
of the same order. To summarize, we consider the regime where the overall
variation is of order $1/\sqrt{n}$.

##### A simple model for perturbations.

What is a reasonable model for the distributional perturbation? For an event
$A$ we will think about $\mathbb{P}^{\xi}[A]$ as a random variable that
depends on the (random) perturbation $\xi$. Perhaps the most simplistic form
of distributional perturbation is to randomly up- or down-weight probabilities
of events. For example, for a uniform distribution
$\mathbb{P}(D=k)=\frac{1}{K}$ for $k=1,\ldots,K$ and an exchangeable set of
non-negative weights $\\{\xi_{k}\\}_{k=1,\ldots,K}$ with
$\sum_{k=1}^{K}\xi_{k}=1$ we can define

$\mathbb{P}^{\xi}(D=k)=\xi_{k}.$ (5)

If the data $D_{i}$ is drawn i.i.d. from the perturbed distribution
$\mathbb{P}^{\xi}$, then empirical means are unbiased marginally: For any
function $\psi(\bullet)$

$E[\frac{1}{n}\sum_{i=1}^{n}\psi(D_{i})]=\mathbb{E}[\psi(D)].$

On the left-hand side, we take the expectation over both the randomness due to
drawing from $\mathbb{P}^{\xi}$ and the randomness in $\xi$. Thus, marginally,
empirical means are unbiased for the respective population quantities. Let us
turn to computing the variance of empirical means. Define
$s(k)=\sum_{i=1}^{n}1_{D_{i}=k}$. Note that under the perturbation model
above, the $s(k)$, $k=1,\ldots,K$ are exchangeable marginally. A short
calculation (see Section B.1) shows that

$\displaystyle\begin{split}\text{Var}_{P}(\frac{1}{n}\sum_{i=1}^{n}\psi(D_{i}))&=\frac{K^{2}}{n^{2}(K-1)}\text{Var}_{P}(s(1))\text{Var}_{\mathbb{P}}(\psi(D))\\\
&=\frac{\delta^{2}}{n}\text{Var}_{\mathbb{P}}(\psi(D)),\end{split}$ (6)

where $\delta^{2}=\frac{K^{2}}{n(K-1)}\text{Var}_{P}(s(1))$ and
$D\sim\mathbb{P}$. Interestingly, compared to equation (2), the variance in
(6) contains a scaling factor $\delta$ that does not depend on $\psi$. Since
$\xi$ and $\mathbb{P}$ are unknown in general, also $\delta$ is usually
unknown.

At this point, the reader might wonder whether this property is specific for
the considered example, and whether a slightly different type of perturbation
model for numerical or discrete data would have led to a different change in
the variance. In other words, the reader might wonder, whether this type of
perturbation model is “natural”. To not interrupt the flow of exposition, we
discuss this question in more detail in the Appendix, Section A.1. The short
answer is that our perturbation model arises naturally under a symmetry
assumption. The following result is a direct consequence of Theorem 2 in the
Appendix.

######  Corollary 1 (Uniqueness of perturbation model, informal).

Let $D_{1},\ldots,D_{n}$ be drawn from $\mathbb{P}^{\xi}$. Under a symmetry
assumption and some regularity conditions on $\mathbb{P}^{\xi}$ (stated in
Theorem 2), there exists $\delta>0$, such that for all $\psi\in
L^{2}(\mathbb{P})$ marginally across distributional and sampling uncertainty
we have

$\displaystyle\text{Var}_{P}(\frac{1}{n}\sum_{i=1}^{n}\psi(D_{i}))$
$\displaystyle=\frac{\delta^{2}}{n}\text{Var}_{\mathbb{P}}(\psi(D)).$

The considered symmetry assumptions also directly imply that
$E[\frac{1}{n}\sum_{i=1}^{n}\psi(D_{i})]=\mathbb{E}[\psi(D)]$.

How should we think about the parameter $\delta$? Intuitively, the larger
$\delta$ is, the further away $\mathbb{P}^{\xi}$ is from $\mathbb{P}$. Thus,
assumptions that hold for the unperturbed distribution $\mathbb{P}$ will be
increasingly violated under $\mathbb{P}^{\xi}$. For example, for large
$\delta$ one can think about this model as inducing sampling bias. The
perturbations potentially induce other types of issues such as correlated
data, and confounding. We discuss this in more detail in the Appendix, Section
A.3.

To conduct statistical inference, it is not sufficient to just compute the
variance. We have to characterize the distribution of empirical processes
$\frac{1}{n}\sum_{i=1}^{n}\psi(D_{i})$. In the following, we will see that if
we think about the distributional change as a superposition of small changes,
then marginally across both sampling uncertainty and distributional
uncertainty we get a central limit theorem, albeit with a non-standard
variance formula.

##### CLT for the isotropic perturbation model.

We have already introduced a random perturbation model for the special case of
discrete data and for a fixed perturbation size. Now we will construct a
general isotropic perturbation model for multivariate continuous or discrete
random variables and consider the case where the change in measure is a
superposition of small incremental changes. We will see that under such a
model, we will get a non-standard CLT in the sense that sample means are
asymptotically normal, with a different variance formula compared to the
i.i.d. case.

To recap, in the random perturbation model, the data is not directly drawn
from the target distribution $\mathbb{P}$, but from some random probability
measure $\mathbb{P}^{\xi}$, where $\mathbb{P}^{\xi}$ is close to $\mathbb{P}$.
The idea is that due to small batch effects, measurement error, or convenience
sampling, the actual sampling distribution $\mathbb{P}^{\xi}$ randomly differs
from the target distribution $\mathbb{P}$. Under $\mathbb{P}^{\xi}$,
probabilities of events are slightly up-weighted or down-weighted compared to
$\mathbb{P}$.

We want to construct a random perturbation model that includes many commonly
encountered situations such as distributions on $\mathbb{R}^{d}$ or the
(infinite-dimensional) space of continuous functions on $\mathbb{R}$. A result
from probability theory shows that any random variable $D$ on a finite or
countably infinite dimensional probability space can be written as a
measurable function $D\stackrel{{\scriptstyle d}}{{=}}h(U)$, where $U$ is a
uniform random variable on $[0,1]$.111For any Borel-measurable random variable
$D$ on a Polish (separable and completely metrizable) space $\mathcal{D}$,
there exists a Borel-measurable function $h$ such that
$D\stackrel{{\scriptstyle d}}{{=}}h(U)$ where $U$ follows the uniform
distribution on $[0,1]$ [Dudley, 2018].

Thus, without loss of generality we will construct distributional
perturbations for a uniform distribution on $[0,1]$. With the transformation
$h(\cdot)$ defined above, this construction generalizes to the general cases
by setting

$\mathbb{P}^{\xi}(D\in\bullet)=\mathbb{P}^{\xi}(h(U)\in\bullet).$

Let us now construct the distributional perturbation for a uniform random
variable. Take $m$ bins $I_{k}=[(k-1)/m,k/m]$ for $k=1,\dots,m$. Let
$W_{1},\dots,W_{m}$ be i.i.d. positive random variables with finite variance.
Set $\xi=(W_{1},\ldots,W_{m})$. We define the randomly perturbed distribution
$\mathbb{P}^{\xi}$ by setting

$\mathbb{P}^{\xi}(U\in\bullet)=\sum_{k}\mathbb{P}(U\in
I_{k}\cap\bullet)\cdot\frac{W_{k}}{\sum_{k=1}^{m}W_{k}/m}.$

Let $m=m(n)$ such that $\frac{n}{m(n)}$ converges to some limit
$r\in(0,\infty)$. Note that $\xi$ depends on $m$ and thus also on $n$.
Conditionally on $\xi$, let $(D_{1}^{n},\ldots,D_{n}^{n})$ be i.i.d. draws
from $\mathbb{P}^{\xi}$.

######  Lemma 1 (Asymptotic behaviour of means).

Under the assumptions mentioned above, for any Borel-measurable square-
integrable function $\psi:\mathcal{D}\mapsto\mathbb{R}^{l}$, we have

$\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\psi(D_{i}^{n})-\mathbb{E}[\psi(D)]\stackrel{{\scriptstyle
d}}{{=}}\mathcal{N}(0,\delta^{2}\text{Var}_{\mathbb{P}}(\psi(D)))+o_{p}(1),$
(7)

with

$\delta^{2}=1+\frac{r\text{Var}(W)}{E[W]^{2}}.$ (8)

In other words, the marginal distribution of
$\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\psi(D_{i}^{n})$ is asymptotically Gaussian
with asymptotic variance containing a scaling factor $\delta$.

The proof can be found in the Appendix, Section B.2. The reason we consider a
triangular array of data sets is that we want both the sampling uncertainty
and distributional uncertainty to be of the same order (see discussion above).
As $n$ goes to infinity we thus have to consider distributional perturbations
that go to zero. We believe that the proposed methodology is also useful in
cases where the distributional perturbations have constant order, but a formal
investigation of this setting is beyond the scope of this work.

The asymptotic behaviour shown in equation (7) arises not only under the
distributional perturbation model but other types of sampling procedures that
induce dependence between observations. In Appendix B.3, we discuss other
sampling models that give rise to (7).

Unless explicitly mentioned otherwise, in the following we assume that the
data scientist has access to one such data set $(D_{1}^{n},\ldots,D_{n}^{n})$
for some large $n$. Note that if the data $(D_{1},\ldots,D_{n})$ is
$\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\mathbb{P}$, then equation (7)
holds for $\delta=1$ and $D_{i}^{n}=D_{i}$. Thus, equation (7) is weaker than
assuming that the data is drawn i.i.d. from $\mathbb{P}$.

In the following, we will discuss how estimators behave asymptotically under
equation (7). It turns out that under some regularity assumptions, maximum
likelihood estimators are still consistent and asymptotically normal, but with
the scaling factor $\delta^{2}$ in the variance formula.

##### Asymptotic behaviour of M-estimators.

Here we will consider the asymptotic behaviour of estimators of type
$\hat{\theta}=\arg\min_{\theta\in\Omega}\frac{1}{n}\sum_{i=1}^{n}L(\theta,D_{i}^{n})$
for a target defined via

$\theta^{0}=\arg\min_{\theta\in\Omega}\mathbb{E}[L(\theta,D)],$ (9)

where $L(\theta,\bullet)$ is a Borel-measurable loss function and $\Omega$ is
an open subset of $\mathbb{R}^{d}$. These estimators include maximum
likelihood estimators with $L(\theta,D)=-\log p_{\theta}(D)$. In classical
statistical theory, uncertainty quantification is usually based on showing
that the estimator is asymptotically Gaussian. Since we have a different two-
stage sampling model, one has to verify that a similar approximation – with a
different variance formula – still holds in our setting.

First, we will discuss consistency. Instead of aiming for maximal generality,
we will adapt a simple consistency proof from the literature. In particular,
we will adapt the classical consistency result in van der Vaart [2000],
Section 5.2.1. We expect that other consistency proofs can be adapted
similarly. The main difference in the proof is that since the data is not
i.i.d. from the target distribution we cannot directly rely on the law of
large numbers. The proof can be found in Appendix B.7.

######  Lemma 2 (Consistency of M-estimators).

Consider the M-estimator
$\hat{\theta}=\arg\min_{\theta\in\Omega}\frac{1}{n}\sum_{i=1}^{n}L(\theta,D_{i}^{n})$
and the target $\theta^{0}=\arg\min_{\theta\in\Omega}\mathbb{E}[L(\theta,D)]$,
where $\Omega$ is a compact subset of $\mathbb{R}^{d}$. Furthermore assume
that $\theta\mapsto L(\theta,D)$ is continuous and that
$\inf_{\|\theta-\theta^{\prime}\|_{2}\leq\delta}L(\theta,D)$ is square-
integrable under $\mathbb{P}$ for every $\delta$ and $\theta^{\prime}$ and
that $\inf_{\theta\in\Omega}L(\theta,D)$ is square integrable. We assume that
$\mathbb{E}[L(\theta,D)]$ has a unique minimum. Then,

$\hat{\theta}-\theta^{0}=o_{p}(1).$ (10)

Now let us turn to asymptotic normality. We will modify the proof in van der
Vaart [2000], Section 5.6. Similarly as above, the main difference in the
proof is since the data is not i.i.d. from the target distribution we cannot
directly rely on the law of large numbers or a standard CLT. The proof can be
found in Appendix B.7.

######  Lemma 3 (Asymptotic normality of M-estimators).

For each $\theta$ in an open subset of $\Omega$, let
$\theta\mapsto\partial_{\theta}L(\theta,d)$ be twice continuously
differentiable in $\theta$ for every $d$. Assume that the matrix
$\mathbb{E}[\partial_{\theta}^{2}L(\theta^{0},D)]$ exists and is nonsingular.
Assume that third order partial derivatives of $\theta\mapsto L(\theta,d)$ are
dominated by a fixed function $h(\cdot)$ for every $\theta$ in a neighborhood
of $\theta^{0}$. We assume that $\partial_{\theta}L(\theta^{0},D)$,
$\partial_{\theta}^{2}L(\theta^{0},D)$ and $h(D)$ are square-integrable under
$\mathbb{P}$. Let
$\hat{\theta}=\arg\min\frac{1}{n}\sum_{i=1}^{n}L(\theta,D_{i}^{n})$. Assume
that $\hat{\theta}-\theta^{0}=o_{p}(1)$, where $\theta^{0}$ satisfies the
estimating equation $\mathbb{E}[\partial_{\theta}L(\theta^{0},D)]=0$. Then,

$\sqrt{n}(\hat{\theta}-\theta^{0})=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}[\partial_{\theta}^{2}L(\theta^{0},D)]^{-1}\partial_{\theta}L(\theta^{0},D_{i}^{n})+o_{p}(1).$
(11)

In particular, by Lemma 1 we have that $\sqrt{n}(\hat{\theta}-\theta^{0})$ is
asymptotically normal with mean zero and covariance matrix $\delta^{2}\Sigma$,
where

$\Sigma=\mathbb{E}[\partial_{\theta}^{2}L(\theta^{0},D)]^{-1}\mathbb{E}[\partial_{\theta}L(\theta^{0},D)^{\intercal}\partial_{\theta}L(\theta^{0},D)]\mathbb{E}[\partial_{\theta}^{2}L(\theta^{0},D)]^{-1}.$

The upshot is that $M$-estimators are asymptotically unbiased, marginally
across both sampling uncertainty and distributional uncertainty. However, the
variance formula changes in the sense that there is an (unknown) scaling
factor $\delta^{2}$.

##### The standard mode of inference fails.

Let us quickly sketch why the standard mode of inference fails. Let’s consider
the case of estimating the mean $\theta^{0}=\mathbb{E}[D]$ via
$\hat{\theta}=\frac{1}{n}\sum_{i=1}^{n}D_{i}^{n}$. One may be tempted to use
the standard variance estimate $\hat{\sigma}_{\text{naive}}^{2}/n$, where

$\hat{\sigma}^{2}_{\text{naive}}=\frac{1}{n}\sum_{i=1}^{n}(D_{i}^{n}-\frac{1}{n}\sum_{j=1}^{n}D_{j}^{n})^{2}.$
(12)

However, a short calculation shows that

$\hat{\sigma}^{2}_{\text{naive}}=\frac{1}{n}\sum_{i=1}^{n}(D_{i}^{n}-\frac{1}{n}\sum_{j=1}^{n}D_{j}^{n})^{2}=\frac{1}{n}\sum_{i=1}^{n}(D_{i}^{n})^{2}-(\frac{1}{n}\sum_{j=1}^{n}D_{j}^{n})^{2}=\text{Var}_{\mathbb{P}}(D)+o_{P}(1).$
(13)

Here, we used equation (7) for $\psi(D)=D$ and $\psi(D)=D^{2}$. However, as
shown in Lemma 1, the asymptotic variance of $\hat{\theta}$ is
$\frac{\delta^{2}}{n}\text{Var}_{\mathbb{P}}(D)$. Thus, the standard approach
drastically underestimates variance in our setting. If $\delta$ is known, one
can simply stretch the confidence intervals discussed in Section 1.3
accordingly. However, in general $\delta$ will be unknown and has to be
estimated from data. We will discuss the estimation of $\delta$ in Section 3.

##### Notation.

Let $\mathbb{P}$ denote a probability measure on $\mathcal{D}$. For each fixed
$n$ the random variable $\xi(n)\in\Xi$ encodes the distributional
perturbation. Formally, $\mathbb{P}^{\bullet}$, $\bullet\in\Xi$, is a
stochastic kernel with the target space $\mathcal{D}$. Conditionally on
$\xi(n)$, we draw an i.i.d. sample $(D_{1}^{n},\ldots,D_{n}^{n})$ from
$\mathbb{P}^{\xi(n)}$. $\xi(n)$ might depend on $n$ but we suppress this in
the notation and simply write $\xi$. Similarly, we sometimes suppress the
dependence of $(D_{1}^{n},\ldots,D_{n}^{n})$ on $n$ and simply write
$(D_{1},\ldots,D_{n})$. We write $P$ for the marginal distribution of
$(D_{1},\ldots,D_{n},\xi)$. We write $\text{Var}_{P}$ for the variance under
$P$ and $\text{Var}_{\mathbb{P}}$ for the variance under $\mathbb{P}$.

##  3 Calibrated inference

We will now discuss how to estimate $\delta$ and form asymptotically valid
confidence intervals for $\theta^{0}$. As discussed earlier, data analysts
often have not just one reasonable estimator for a given parameter
$\theta^{0}$, but potentially several reasonable estimators
$\hat{\theta}^{1},\ldots,\hat{\theta}^{K}$. For example, these estimators can
arise from using different specifications in generalized linear models or by
running the analysis for subgroups of the observations.

######  Example 1 (OLS with several specifications).

Let us consider a setting in which the data analyst wants to estimate the
causal effect of some variable $X_{1}$ on a target variable $Y$. On
observational data, this is often done by invoking suitable assumptions and
regressing $Y$ on $X_{1}$ and a suitable set of covariates. Often, the analyst
has several reasonable choices for the set of covariates. Suppose that the
data analyst performs ordinary least-squares on $K$ different subsets of $X$
that include $X_{1}$, denoted by $X^{S_{1}},X^{S_{2}},\dots,X^{S_{K}}$. For
example, $X^{S_{1}}$ can be $(X_{1},X_{2},X_{3})$. Now the data analyst has
$K$ different regression coefficients of $X_{1}$,
$\hat{\theta}^{1},\dots,\hat{\theta}^{K}$ where

$\hat{\theta}^{k}=(\sum_{i=1}^{n}X_{i}^{S_{k}}{(X_{i}^{S_{k}})}^{\intercal})_{1,\bullet}^{-1}\sum_{i=1}^{n}X_{i}^{S_{k}}Y_{i}.$

If the empirical variation between the estimators
$\hat{\theta}^{1},\ldots,\hat{\theta}^{K}$ is low, then the analyst may feel
more confident about conclusions drawn from these estimates than if the
variation between these estimators is very large. As an example, in Chiappori
et al. [2012] the authors write “It is reassuring that the estimates are very
similar in the standard and the augmented specifications”. We will now look at
this practice under the isotropic perturbation model. We will see that in this
setting it is possible to construct a consistent estimator of $\delta$ and
form asymptotically valid confidence intervals that account for both sampling
uncertainty and distributional perturbations.

Let us investigate the behaviour of these estimators in more detail. First, if
the estimators $\hat{\theta}^{k}=\hat{\theta}^{k}(D_{1},\ldots,D_{n})$ are
M-estimators, by Lemma 2 and Lemma 3 the estimators are asymptotically linear
in the sense that

$\hat{\theta}^{k}-\theta^{k}=\frac{1}{n}\sum_{i=1}^{n}\phi^{k}(D_{i})+o_{p}(\frac{1}{\sqrt{n}})$
(14)

for some deterministic $\theta^{k}=\arg\min\mathbb{E}[L^{k}(\theta,D)]$, where
$L^{k}$ is the loss function of estimator $\theta^{k}$. Here the functions
$\phi^{k}$ are assumed to satisfy $\mathbb{E}[\phi^{k}(D)]=0$ and
$\text{Var}_{\mathbb{P}}(\phi^{k}(D))\in(0,\infty)$. Since $\phi^{k}(D)$ is
square integrable, by Lemma 1 the sequence
$\sqrt{n}(\hat{\theta}^{k}-\theta^{k})$ is asymptotically normal with mean
zero and covariance $\delta^{2}\text{Var}_{\mathbb{P}}(\phi^{k}(D))$. We
summarize this behaviour of the estimators as the following assumption for the
convenience of reference later.

######  Assumption 1 (Asymptotic linearity).

The estimators $\hat{\theta}^{k}$, $k=1,\ldots,K$ are asymptotically linear,
that is they satisfy equation (14) for functions $\phi^{k}$ with
$\mathbb{E}[\phi^{k}(D)]=0$ and
$0<\text{Var}_{\mathbb{P}}(\phi^{k}(D))<\infty$.

As discussed above, for the case of $M$-estimators, this assumption can be
justified via Lemma 2 and Lemma 3. We will now formalize the premise that the
data analyst considers each of the $\hat{\theta}^{k}$ a reasonable estimator
for the parameter of interest, $\theta^{0}$.

######  Assumption 2 (Agreement).

We have $\theta^{k}=\theta^{0}$ for $k=1,\ldots,K$.

This assumption must be justified with scientific background knowledge.
Intuitively, the assumption states that if both sampling uncertainty and
distributional uncertainty were negligible, the estimators would agree. In
Section 4, we discuss in a numerical example how the choice of such estimators
can be justified and study the behaviour of our method in cases where the
estimators do not agree asymptotically. If the data scientist does not believe
in asymptotic agreement of the estimators, we present conservative confidence
intervals in the Appendix, Section C.

###  3.1 Confidence intervals

Now let us turn to constructing confidence intervals for $\theta^{0}$. Assume
that the data analyst has access to $K$ different estimators
$\hat{\theta}^{1},\dots,\hat{\theta}^{K}$ that are asymptotically linear for
estimating $\theta^{0}$ with influence functions
$\phi^{1}(D),\dots,\phi^{K}(D)$. For expository simplicity, for now we assume
that their influence functions $\phi^{1}(D),\dots,\phi^{K}(D)$ are
uncorrelated and have the same variance $\sigma^{2}>0$ under $\mathbb{P}$.
Later in the section, we discuss how to construct confidence intervals for
general cases where influence functions are possibly correlated and have
different variances. Since the estimators are asymptotically unbiased,
uncorrelated, and have the same variance, as final estimate we consider the
mean of estimators,
$\hat{\theta}^{\text{pooled}}=\frac{1}{K}\sum_{k}\hat{\theta}^{k}$. In the
following, we will investigate the asymptotic behaviour of this estimator.

By Assumption 1, 2 and Lemma 1, for $k=1,\dots,K$,

$\sqrt{n}(\hat{\theta}^{k}-\theta^{0})\stackrel{{\scriptstyle
d}}{{=}}\delta\sigma Z_{k}+o_{P}(1),$

where $Z_{k}$ are independent standard normal random variables. Thus,

$\sqrt{n}(\hat{\theta}^{\text{pooled}}-\theta^{0})\stackrel{{\scriptstyle
d}}{{=}}\delta\sigma\bar{Z}+o_{P}(1).$ (15)

On the other hand, define the between-estimator variance

$\hat{\sigma}_{\text{bet}}^{2}=\frac{1}{K-1}\sum_{k=1}^{K}(\hat{\theta}^{k}-\hat{\theta}^{\text{pooled}})^{2}.$

Then,

$\hat{\sigma}_{\text{bet}}^{2}\stackrel{{\scriptstyle
d}}{{=}}\frac{\delta^{2}\sigma^{2}}{n}\frac{1}{K-1}\sum(Z_{k}-\bar{Z})^{2}+o_{P}(1/n)\stackrel{{\scriptstyle
d}}{{=}}\frac{\delta^{2}\sigma^{2}}{n}\chi^{2}(K-1)+o_{P}(1/n)$ (16)

where $\chi^{2}(K-1)$ is a chi-square random variable with $K-1$ degrees of
freedom. Let us assume for a moment that $\sigma^{2}$ is known to the data
scientist. In this case, the data scientist may estimate $\hat{\delta}^{2}$
via

$\hat{\delta}^{2}:=\frac{n\hat{\sigma}_{\text{bet}}^{2}}{(K-1)\sigma^{2}}\stackrel{{\scriptstyle
d}}{{=}}\delta^{2}\frac{\chi^{2}(K-1)}{K-1}+o_{P}(1).$

Combining equations (15) and (16), we get

$\frac{\hat{\theta}^{\text{pooled}}-\theta^{0}}{\hat{\sigma}_{\text{bet}}/\sqrt{K-1}}\stackrel{{\scriptstyle
d}}{{=}}t(K-1)+o_{p}(1).$

where $t(K-1)$ is a $t$-distributed random variable with $K-1$ degrees of
freedom. Note that $\delta$, $\sigma$ cancel out.

Without direct estimation of $\delta$ or $\sigma$, we have an $1-\alpha$
confidence interval of $\theta^{0}$:

$\hat{\theta}^{\text{pooled}}\pm
t_{K-1,1-\alpha/2}\frac{\hat{\sigma}_{\text{bet}}}{\sqrt{K-1}},$

where $t_{K-1,1-\alpha/2}$ is the $1-\alpha/2$ quantile of the
$t$-distribution with $K-1$ degrees of freedom. Note that the size of the
confidence intervals goes to zero with rate $1/\sqrt{n}$ as
$\hat{\sigma}_{\text{bet}}=O_{P}(1/\sqrt{n})$. Computing the between-estimator
variance has some similarities with random effect models. We discuss some of
the differences in Remark 3.

Let us make this argument more general. We will now discuss the case where the
estimators $\hat{\theta}^{k}$ have potentially different asymptotic variances
$\text{Var}_{\mathbb{P}}(\phi^{k}(D))$. Instead of using
$\hat{\theta}^{\text{pooled}}=\frac{1}{K}\sum_{k}\hat{\theta}^{k}$ as the
final estimate, we recommend inverse variance weighting. Thus, we first need
to estimate $\text{Var}_{\mathbb{P}}(\phi^{k}(D))$ consistently. While
estimating $\text{Var}_{\mathbb{P}}(\phi^{k}(D))$ is straightforward under
i.i.d. sampling, we also have to verify that this works in our model class. We
estimate $\text{Var}_{\mathbb{P}}(\phi^{k}(D))$ using plug-in estimators of
the influence function $\hat{\phi}^{k}(D)$ as

$\widehat{\text{Var}_{\mathbb{P}}(\phi^{k}(D))}=\frac{1}{n}\sum_{i=1}^{n}\Big{(}\hat{\phi}^{k}(D_{i})-\frac{1}{n}\sum_{i=1}^{n}\hat{\phi}^{k}(D_{i})\Big{)}^{2}.$
(17)

The following proposition shows that
$\widehat{\text{Var}_{\mathbb{P}}(\phi^{k}(D))}$ is a consistent estimator of
$\text{Var}_{\mathbb{P}}(\phi^{k}(D))$.

######  Proposition 1 (Consistency of
$\widehat{\text{Var}_{\mathbb{P}}(\phi^{k}(D))}$).

Suppose that the $\phi^{k}(D)$ has finite fourth moments. Furthermore, suppose
that the estimation of the influence function is consistent in the sense that

$\frac{1}{n}\sum_{i=1}^{n}(\hat{\phi}^{k}(D_{i})-\phi^{k}(D_{i}))^{2}=o_{p}(1).$
(18)

Then, $\widehat{\text{Var}_{\mathbb{P}}(\phi^{k}(D))}$ defined in (17)
satisfies that

$\widehat{\text{Var}_{\mathbb{P}}(\phi^{k}(D))}=\text{Var}_{\mathbb{P}}(\phi^{k}(D))+o_{p}(1).$

######  Remark 1 (OLS).

Note that equation (18) is expected to hold for the plug-in estimators of the
influence function under regularity assumptions. Revisiting Example 1, a plug-
in estimator of the influence function is

$\hat{\phi}^{k}(D_{i})=(\frac{1}{n}\sum_{j=1}^{n}X_{j}^{S_{k}}{(X_{j}^{S_{k}})}^{\intercal})_{1,\bullet}^{-1}X_{i}^{S_{k}}(Y_{i}-(X_{i}^{S_{k}})^{\intercal}\hat{\theta}^{k,OLS}),$

where $\hat{\theta}^{k,OLS}$ is the OLS estimator computed with covariates
$X^{S^{k}}$. One can now justify equation (18) via Lemma 2.

Now we construct asymptotically valid confidence intervals for $\theta^{0}$
using $K$ different estimators $\hat{\theta}^{1},\dots,\hat{\theta}^{K}$ that
are asymptotically linear for estimating $\theta^{0}$. In the following
theorem with Remark 2, influence functions of $K$ different estimators can be
correlated and have different variances.

######  Theorem 1.

(Asymptotic validity of calibrated confidence interval). Suppose Assumption 1
and 2 hold and the influence functions $\phi^{1}(D),\dots,\phi^{K}(D)$ are
uncorrelated. Let
$\hat{\theta}^{W}=\sum_{k=1}^{K}\hat{\alpha}_{k}\hat{\theta}^{k}$ be the
inverse-variance weighted estimator where the weights are

$\hat{\alpha}_{k}=\frac{\frac{1}{\widehat{\text{Var}_{\mathbb{P}}(\phi^{k}(D))}}}{\sum_{j=1}^{K}\frac{1}{\widehat{\text{Var}_{\mathbb{P}}(\phi^{j}(D))}}},$
(19)

with
$\widehat{\text{Var}_{\mathbb{P}}(\phi^{k}(D))}=\text{Var}_{\mathbb{P}}(\phi^{k}(D))+o_{p}(1)$
for $k=1,\dots,K$. Let $\hat{\sigma}_{\text{bet}}$ be the weighted between-
estimator variance defined as

$\hat{\sigma}^{2}_{\text{bet}}=\sum_{k}\hat{\alpha}_{k}(\hat{\theta}^{k}-\hat{\theta}^{W})^{2}.$
(20)

Then for any $\alpha\in(0,1)$, for fixed $K$ and as $n\rightarrow\infty$ we
have

$\mathrm{P}\left(\theta^{0}\in\Big{[}\hat{\theta}^{W}\pm
t_{K-1,1-\alpha/2}\cdot\frac{\hat{\sigma}_{\text{bet}}}{\sqrt{K-1}}\Big{]}\right)\xrightarrow[]{}1-\alpha,$
(21)

where $t_{K-1,1-\alpha/2}$ is the $1-\alpha/2$ quantile of the $t$
distribution with $K-1$ degrees of freedom. To be clear, here we marginalize
over both the randomness due to sampling and the randomness due to the
distributional perturbation.

######  Remark 2 (Correlated estimators).

In practice, the components of $(\phi^{1}(D),\dots,\phi^{K}(D))$ may be
correlated. Then, we can apply a linear transformation to the estimators to
obtain uncorrelated estimators that are asymptotically unbiased for
$\theta^{0}$. We define the transformation matrix
$T_{ij}=\frac{(\hat{\Sigma}^{-1/2})_{ij}}{\sum_{j^{\prime}}(\hat{\Sigma}^{-1/2})_{ij^{\prime}}}$,
where $\hat{\Sigma}$ is an estimate of the covariance matrix of
$\hat{\theta}^{1},\ldots,\hat{\theta}^{K}$. We can then define
$(\hat{\eta}^{1},\ldots,\hat{\eta}^{K})^{\intercal}=T\cdot(\hat{\theta}^{1},\ldots,\hat{\theta}^{K})^{\intercal}$.
If $\|\hat{\Sigma}-\Sigma\|_{2}=o_{p}(1)$ and $\Sigma$ is invertible, then the
estimators $\hat{\eta}^{1},\ldots,\hat{\eta}^{K}$ also satisfy Assumption 1
with influence functions that are pairwise uncorrelated. Furthermore, if the
$\hat{\theta}^{k}$, $k=1,\ldots,K$ satisfy Assumption 2, then also
$\hat{\eta}^{k}$, $k=1,\ldots,K$ satisfy Assumption 2.

######  Remark 3 (Random effect models).

The proposed method shares some similarity with random effect models, since
both can be seen as investigating variability across estimators. In a random
effect model, estimator-to-estimator variability can be modeled as random
additive terms with a correlation structure given by a user-defined model
matrix. This model matrix is usually justified by background knowledge. In our
setup, we capture estimator-to-estimator variability through a multiplicative
term. This multiplicative term is justified by symmetry assumptions on the
sampling and perturbation process. We believe that these procedures will be
used in different scenarios: if the data scientist has sufficient background
knowledge to model the dependency of biases of $\hat{\theta}^{k}$ with a model
matrix, one can use a random effect model. If such knowledge is not available,
using a random effect model is prohibited. In this case, in our view adjusting
confidence intervals with scaling factors is more justifiable than assuming
that the data is i.i.d. from the unperturbed distribution $\mathbb{P}$.

In some cases, the data analyst may trust one of the estimators
$\hat{\theta}^{k}$ more than others. For example, the data analyst may be
convinced that $\theta^{1}=\theta^{0}$ but may not be sure whether
$\theta^{k}=\theta^{0}$ for $k\geq 2$. In this case it is possible to
construct variance estimates that are upwardly biased in the sense that the
resulting confidence intervals are expected to be conservative. In this case,
the data analyst may report the confidence interval for $\theta^{0}$ using
$\hat{\theta}^{1}$ instead of $\hat{\theta}^{W}$ with $\delta$ estimated by
computing the between-estimator variance of the remaining $K-1$ estimators. As
a result, they would lose one degree of freedom in their confidence intervals.
The details can be found in the Appendix, Section C.

In the final variance estimate, there are two effects that are counteracting
each other. Inverse-variance weighting reduces the variance of the final
estimate compared to each of the individual estimators $\hat{\theta}^{k}$. On
the other hand, the new variance formula accounts for distributional
uncertainty and thus potentially inflates the variance.

##### Practical implications for stability analyses.

The proposed model not only leads to a recommendation on how to summarize the
between-estimator uncertainty in confidence intervals but also lets us give
some additional guidance.

First, note that if all estimators $\hat{\theta}^{k}$ have similar influence
functions, the proposed method will be unstable since in Remark 2 we invert
the estimated covariance matrix. This coincides with the following intuition:
Reporting that a large number of extremely similar estimators return similar
results does not automatically increase the trustworthiness of a result. To
corroborate a hypothesis one should have estimators that are susceptible to
different sources of biases. In our model, this corresponds to estimators that
are not highly correlated. Ideally, the estimators are independent. Similar
arguments have appeared in other parts of the literature. For example,
Rosenbaum [2021] writes: “An observational study has two evidence factors if
it provides two comparisons susceptible to different biases that may be
combined as if from independent studies of different data by different
investigators, despite using the same data twice”.

In practice, it may happen that calibrated confidence intervals are very large
compared to traditional sampling-based confidence intervals. Apart from the
estimation error and small $K$, there are two possible explanations.

First, it could be that distributional uncertainty is very large. If
distributional uncertainty is much larger than sampling uncertainty,
conventional (unadjusted) confidence intervals are of limited value. Similar
points have been made in different parts of the literature. For example, Meng
[2018] argues that as the sample size grows, data quality becomes more
important than data quantity and that standard confidence intervals have to be
inflated to account for issues of data quality. In this vein, distributional
confidence intervals can be used as a warning signal that we might be in a
regime where data quality issues are more pressing than sampling uncertainty.

Secondly, it could be that the assumptions are violated (that means
$\theta^{k}\neq\theta^{k^{\prime}}$ for some $k,k^{\prime}$). If the
assumptions are grossly violated, inference will be more conservative. If the
assumptions are correct, inference will be more precise. In other words, the
precision of calibrated inference depends on whether Assumption 2 is satisfied
or not.

If the number of estimators $K$ is very small, then there is an inferential
price to pay for estimating distributional uncertainty in terms of power. This
is reflected in the degrees of freedom of the $t$-distribution.

##  4 Simulation study

In this section, we evaluate the performance of the proposed method via a
simulation study. The estimation accuracy of $\delta$ and the marginal
coverage of calibrated confidence intervals are evaluated on simulated data
sets generated by random perturbation models. In this simulation, we emulate
the situation where a data scientist uses linear regression with an adjustment
set to estimate a causal effect.

Setup. The unperturbed distribution of $D=(X,Y)$ with covariates
$X\in\mathbb{R}^{5}$ and response $Y\in\mathbb{R}$ is generated from the
following structural causal model [Bollen, 1989, Pearl, 2009]:

$\epsilon,\epsilon_{1},\epsilon_{2},X_{3},X_{4},X_{5}\overset{\text{i.i.d}}{\sim}\mathcal{N}(0,1)$,
---
$X_{2}\leftarrow X_{3}+\epsilon_{2}$,
$X_{1}\leftarrow 0.5X_{2}+X_{4}+\epsilon_{1}$,
$Y\,\,\,\leftarrow X_{1}+0.5X_{2}+X_{3}+X_{5}+\epsilon$

The goal is to estimate the direct causal effect of $X_{1}$ on $Y$, which in
this setup corresponds to the regression coefficient of $X_{1}$ in a
regression of $Y$ on the set $S=(X_{1},X_{2})$. Practitioners often conduct
such regressions for different choices of sets $S$ to evaluate the overall
stability of the procedure [Leamer, 1983, Oster, 2019].

In this example, the structural causal model can be used to construct multiple
valid estimators. We look at the case where the data analyst considers $K=8$
different adjustment sets which all include the confounding variable $X_{2}$.
In this case, $K=8$ different regression-adjusted estimators estimate the same
quantity, the direct causal effect of $X_{1}$ on $Y$, under the unperturbed
distribution. We consider following adjustment sets; $\\{X_{1},X_{2}\\}$,
$\\{X_{1},X_{2},X_{3}\\}$, $\\{X_{1},X_{2},X_{4}\\}$,
$\\{X_{1},X_{2},X_{5}\\}$, $\\{X_{1},X_{2},X_{3},X_{4}\\}$,
$\\{X_{1},X_{2},X_{3},X_{5}\\}$, $\\{X_{1},X_{2},X_{4},X_{5}\\}$,
$\\{X_{1},X_{2},X_{3},X_{4},X_{5}\\}$.

We now want to model a random shift between the target and the sampling
distribution due to issues during the sampling procedure. We generate randomly
perturbed data sets in two ways. First, we adopt the random perturbation model
in Lemma 1. We partition the support of the joint distribution of $X$ and $Y$
into $m^{p+1}$ equal probability bins and perturb the probability of each bin
with i.i.d. random weights $Z\cdot W$ where $W\sim\text{Gamma}(1,1)$ and
$Z\sim\text{Ber}(1/m^{p})$. For sufficiently large $m$, this procedure can be
seen as randomly selecting $m$ bins out of $m^{p+1}$ bins and perturbing the
probability of each selected bin with i.i.d random weights
$W\sim\text{Gamma}(1,1)$. In our simulations, we generate $n$ i.i.d. data
points $D_{1},\dots,D_{n}$ from this randomly perturbed distribution. The
strength of the perturbation is given as $\delta^{2}\approx 1+2\cdot n/m$.
Secondly, we employ the random perturbation model described in Example 5 in
the Appendix. We sample $m$ data points from the original distribution and let
randomly perturbed distribution be the empirical distribution of $m$ samples.
The strength of the perturbation is given as $\delta^{2}\approx 1+n/m$.

Our method is carried out for sample sizes $n=200,500,1000$ and for
$m=200,500,1000$ which determines the strength of the perturbation, each with
$N=1000$ replicates. In each replicate, we generate $n$ samples from the
randomly perturbed distribution, obtain $K=8$ different regression-adjusted
estimators from the perturbed data set, and estimate $\hat{\delta}^{2}$. Then
we construct non-calibrated and calibrated $(1-\alpha)$ confidence intervals
for each regression-adjusted estimator. Note that we use
$\max(\hat{\delta},1)$ instead of $\hat{\delta}$ to adjust confidence
intervals. We evaluate the accuracy of $\hat{\delta}$ compared to the ground
truth $\delta$ and the marginal coverage of calibrated confidence intervals
compared to that of non-calibrated confidence intervals for each regression-
adjusted estimator.

Accuracy of $\hat{\delta}$. The estimation accuracy of $\hat{\delta}$ compared
to the ground truth $\delta$ is illustrated in Figure 1. We also look at the
case where the data analyst fails to adjust for the confounding variable
$X_{2}$ in one of the regressions as they consider the adjustment set
$\\{X_{1},X_{3},X_{4}\\}$ instead of $\\{X_{1},X_{2},X_{3},X_{4}\\}$. In this
case, $K=8$ different estimators do not agree asymptotically and therefore
$\hat{\delta}$ is upwardly biased, which leads to overcoverage.

Figure 1: Accuracy of $\hat{\delta}$: The left panel shows results under the
perturbation model in Lemma 1 and the right panel shows results under the
perturbation model in Example 5 from the Appendix. Mean, $2.5\%$ and $97.5\%$
quantiles of the estimated $\hat{\delta}$ for each $m=200,500,1000$ and
$n=200,500,1000$ are provided for two cases when $K$ different estimators
agree asymptotically and when they do not. The dashed lines indicate the true
values of $\delta$.

Marginal coverage of calibrated CI. The marginal coverages of calibrated
confidence intervals and non-calibrated confidence intervals are given in
Figure 2. We see that calibrated confidence intervals have much improved
coverage compared to non-calibrated confidence intervals, especially when $n$
is large and $m$ is small as the variance due to distributional perturbations
dominates the marginal variance.

Figure 2: Marginal coverages of calibrated confidence intervals: The left
panel shows results under the perturbation in Lemma 1 and the right panel
shows results under the perturbation model in Example 5 in the Appendix.
Marginal coverages of calibrated confidence intervals and non-calibrated
confidence intervals are presented for $m=200,500,1000$ and $n=200,500,1000$.
The dashed lines indicate the nominal coverage 0.95.

##  5 Real-world data

Ultimately, the goal of our procedure is to increase stability and
trustworthiness of decision-making. In this section, we demonstrate that our
method can improve stability on a real data set. We will see that even in
situations without distributional perturbations, the proposed method can
increase stability of decision-making. The data set (Cortez and Silva [2008])
was collected by using school reports and questionnaires for the estimation of
final grades of students in secondary education of two Portuguese schools. The
data attributes include student grades, demographic, social and school related
features. It is available at the UCI machine learning repository (Dheeru and
Karra Taniskidou [2017]). We adopt 20 covariates in the data set. The response
$Y$ is the final year grade in Portuguese language. There are 649 students in
total.

The goal is to determine the relative importance of $L=7$ selected binary
covariates: 1) parents’ cohabitation status, 2) whether the student received
extra educational support from the school, 3) whether the student received
family educational support, 4) whether the student is in a relationship, 5)
whether the student had extra paid classes within the course subject, 6)
whether the student’s mother had secondary or higher education, and 7) whether
the student’s father had secondary or higher education. The relative
importance is determined by the rank order of the covariates’ effect sizes in
a linear regression.

In the simulation setup, we want to emulate a situation where as baseline the
analyst has several reasonable choices to conduct a statistical analysis, and
makes these decisions randomly. On the other hand, as a comparison the analyst
aggregates the estimators and conducts uncertainty quantification as proposed
above.

Suppose we are given multiple sets of covariates, all containing the 7
covariates of our interest. We consider the following two methods. In method
1, a statistician randomly chooses one of the sets of covariates, performs a
linear regression, and ranks the effect sizes of 7 covariates. In method 2, a
statistician employs our method. In particular, they perform linear
regressions with multiple sets of covariates and for each covariate, calculate
an inverse-variance weighted estimator and its effect size in consideration of
distributional perturbations as described in Section 3. Then, they rank these
effect sizes.

We evaluate the two methods’ stability in ranking effect sizes. To evaluate
method $i$, we randomly split the data set into two, perform method $i$ on
each split, and compare the rankings resulting from each split. To measure the
stability of the ranking, we compute the set similarity measure between
$S_{1,\ell}=\\{$Top $\ell$ covariates by the effect size on split 1$\\}$ and
$S_{2,\ell}=\\{$Top $\ell$ covariates by the effect size on split 2$\\}$ for
each $\ell=1,\dots,L=7$ as $|S_{1,\ell}\cap S_{2,\ell}|/L$. We repeat this
procedure $N=500$ times and record the average set similarity measure. In each
replicate, we randomly generate $K=10,20$ sets of covariates that include the
7 covariates of our interest. The results can be found in Table 1. Overall, we
see our method (Method 2) improves the stability of the ranking, notably
outperforming Method 1 for $\ell=1,2,3$. Note that the method 1 gives slightly
worse results than random guessing for small $\ell$. One possible explanation
is that sample splitting introduces small negative correlations between
splits: If a regression coefficient is close to zero on the entire data set
and on one split by chance the coefficient is large, then the coefficient is
expected to be small on the other split.

$\ell$ | 1 | 2 | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---|---|---
Method 1 ($K=10$) | 0.102 | 0.203 | 0.407 | 0.648 | 0.817 | 0.898 | 1.000
Method 2 ($K=10$) | 0.210 | 0.296 | 0.449 | 0.658 | 0.828 | 0.912 | 1.000

$\ell$ | 1 | 2 | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---|---|---
Method 1 ($K=20$) | 0.090 | 0.203 | 0.417 | 0.659 | 0.817 | 0.893 | 1.000
Method 2 ($K=20$) | 0.235 | 0.313 | 0.445 | 0.679 | 0.845 | 0.912 | 1.000

Table 1: The stability of the ranking: The table above shows results with
$K=10$ sets of covariates and the table below shows results with $K=20$ sets
of covariates. Mean over $N=500$ iterations of the computed set similarity
measure between $S_{1,\ell}$ and $S_{2,\ell}$ for each $\ell=1,\dots,7$ is
provided for each method.

Additionally, we compare the length of calibrated and non-calibrated
confidence intervals for each selected binary covariate. From the results
provided in Figure 3, one sees that our method is not so conservative given
that we are adjusting confidence intervals with scaling factor
$\max(\hat{\delta},1)$. Moreover, the variance of the length of calibrated
confidence intervals tends to decrease as we increase the number of sets of
covariates from $K=10$ to $K=20$.

Figure 3: The length of calibrated and non-calibrated confidence intervals:
Mean, $2.5\%$ and $97.5\%$ quantiles of the length of calibrated (Method 2)
and non-calibrated (Method 1) confidence intervals for each selected binary
covariate over $N=500$ iterations are provided for $K=10$ and $K=20$.

##  6 Discussion

In practice, data analysts often compute not just one estimator but multiple
estimators for a single target quantity. These estimators are often based on
equally reasonable assumptions. For example, in causal inference from
observational data, it is common to compute regression-adjusted estimators for
different choices of adjustment sets. Often, it is recommended study the
estimator-to-estimator variability between sensible choices of adjustment
sets. If the estimator-to-estimator variability is high, then the analyst
might have reason to not trust the estimates. In these cases, such stability
investigations may be more informative than traditional $p$-values or
confidence regions. This warrants an investigation of the theoretical
properties of this practice. Does this practice have any guarantees and if so,
which? Can we integrate this type of stability analysis into statistical
inference?

We study a variant of this procedure from a distributional perspective. If the
data analyst has access to different estimators that are supposedly estimating
the same quantity, estimator-to-estimator variability can be used to scale
confidence intervals. These scaled confidence intervals can be shown to
account for both sampling uncertainty and distributional uncertainty in an
isotropic perturbation model. This perturbation model assumes that the
distribution shift is a superposition of many small random distributional
changes. Such uncertainty quantification seems desirable, especially in
settings where the sampling uncertainty is of similar or lower order than
other types of uncertainty.

We have shown in the Appendix A.3 that distributional uncertainty is related
to other issues that are common in data analysis, such as violations of the
i.i.d. assumption, confounding, and sampling bias. That being said, the
calibration procedure is not meant to replace existing methods that address
confounding or selection bias via direct modelling (such as regression
adjustment, or weighting procedures, a random effect model, or by bounding
errors). In contrast, the proposed calibration technique allows you to deal
with residual “dense” distributional errors that you are unable to capture
with direct modelling approaches.

The isotropic perturbation model is a strong assumption, but it is a weaker
assumption than assuming that the data is i.i.d. from the target distribution,
which is a commonly made. Thus, the proposed calibration procedure works under
strictly less assumptions on the data generating process than the most common
inferential strategy. Instead of relying on i.i.d. sampling from $\mathbb{P}$,
inference in the proposed model is based on a symmetry assumption and on
scientific background knowledge for finding multiple reasonable estimators.

Of course, in practice perturbations might affect parts of the distribution
differently. In such cases the proposed method can potentially have over-
coverage or under-coverage. Looking forward, it would be desirable to extend
the isotropic perturbation model (which has only one single parameter
$\delta$) to more flexible models that depend on multiple parameters. Such
perturbation models would allow training different uncertainty models for
different parts of the distribution, potentially leading to more realistic and
flexible uncertainty quantification than existing approaches.

A companion R package, calinf, is available at
https://github.com/rothenhaeusler/calinf. Our package allows to draw data
under the distributional uncertainty model and calibrate inference in
generalized linear models. We provide an example of calibrated inference where
the data analyst computes regression-adjusted estimators for different choices
of adjustment sets. If multiple estimators are not available, it is also
possible to estimate $\delta$ using other types of scientific background
knowledge. On the GitHub page, we discuss an example where the data analyst
has background knowledge of population parameters.

##  7 Acknowledgments

We are grateful to several anonymous referees and the associate editor for
valuable comments and suggestions. We thank Bin Yu, Peter Bühlmann, Guido
Imbens, Kevin Guo, Suyash Gupta, Ying Jin, and James Yang for inspiring
discussions.

## References

  * Ben-Tal et al. [2013] A. Ben-Tal, D. Den Hertog, A. De Waegenaere, B. Melenberg, and G. Rennen.  Robust solutions of optimization problems affected by uncertain probabilities.  _Management Science_ , 59(2):341–357, 2013.
  * Bollen [1989] K. Bollen.  _Structural Equations with latent variables_.  John Wiley & Sons, 1989.
  * Box [1980] G. E. Box.  Sampling and bayes’ inference in scientific modelling and robustness.  _Journal of the Royal Statistical Society: Series A_ , 143(4):383–404, 1980.
  * Breiman [1996] L. Breiman.  Bagging predictors.  _Machine learning_ , 24(2):123–140, 1996.
  * Breiman [2001] L. Breiman.  Random forests.  _Machine learning_ , 45(1):5–32, 2001.
  * Bühlmann [2020] P. Bühlmann.  Invariance, causality and robustness.  _Statistical Science_ , 35(3):404–426, 2020.
  * Chiappori et al. [2012] P.-A. Chiappori, S. Oreffice, and C. Quintana-Domeque.  Fatter attraction: anthropometric and socioeconomic matching on the marriage market.  _Journal of Political Economy_ , 120(4):659–695, 2012.
  * Cornfield et al. [1959] J. Cornfield, W. Haenszel, E. Hammond, A. Lilienfeld, M. Shimkin, and E. Wynder.  Smoking and lung cancer: recent evidence and a discussion of some questions.  _Journal of the National Cancer institute_ , 22(1):173–203, 1959.
  * Cortez and Silva [2008] P. Cortez and A. Silva.  Using data mining to predict secondary school student performance.  _Proceedings of 5th Future Business Technology Conference (FUBUTEC 2008)_ , pages 5–12, 2008.
  * Denzin [1970] N. K. Denzin.  _The research act: A theoretical introduction to sociological methods_.  1970\.
  * Dheeru and Karra Taniskidou [2017] D. Dheeru and E. Karra Taniskidou.  UCI machine learning repository, 2017.  URL http://archive.ics.uci.edu/ml.
  * Duchi et al. [2021] J. C. Duchi, P. W. Glynn, and H. Namkoong.  Statistics of robust optimization: A generalized empirical likelihood approach.  _Mathematics of Operations Research_ , 46(3):946–969, 2021.
  * Dudley [2018] R. M. Dudley.  _Real analysis and probability_.  CRC Press, 2018.
  * Freedman [1991] D. Freedman.  Statistical models and shoe leather.  _Sociological methodology_ , pages 291–313, 1991.
  * Heinze-Deml and Meinshausen [2021] C. Heinze-Deml and N. Meinshausen.  Conditional variance penalties and domain shift robustness.  _Machine Learning_ , 110(2):303–348, 2021.
  * Hernán and Robins [2020] M. Hernán and J. Robins.  _Causal Inference: What If_.  CRC, 2020.
  * Huber [1981] P. Huber.  _Robust Statistics_.  Wiley, New York, 1981.
  * Imbens and Rubin [2015] G. Imbens and D. Rubin.  _Causal inference in statistics, social, and biomedical sciences_.  Cambridge University Press, 2015.
  * Karmakar et al. [2019] B. Karmakar, B. French, and D. Small.  Integrating the evidence from evidence factors in observational studies.  _Biometrika_ , 106(2):353–367, 2019.
  * LaLonde [1986] R. LaLonde.  Evaluating the econometric evaluations of training programs with experimental data.  _The American Economic Review_ , pages 604–620, 1986.
  * Leamer [1983] E. Leamer.  Let’s take the con out of econometrics.  _The American Economic Review_ , 73(1):31–43, 1983.
  * Lee and Ogburn [2021] Y. Lee and E. L. Ogburn.  Network dependence can lead to spurious associations and invalid inference.  _Journal of the American Statistical Association_ , 116(535):1060–1074, 2021.
  * Meng [2018] X.-L. Meng.  Statistical paradises and paradoxes in big data (i): Law of large populations, big data paradox, and the 2016 us presidential election.  _The Annals of Applied Statistics_ , 12(2):685–726, 2018.
  * Munafò and Smith [2018] M. Munafò and G. Smith.  Repeating experiments is not enough.  _Nature_ , 553(7689):399–401, 2018.
  * Neyman [1923] J. Neyman.  On the application of probability theory to agricultural experiments. essay on principles. section 9.  _Statistical Science_ , pages 465–472, 1923.
  * Oster [2019] E. Oster.  Unobservable selection and coefficient stability: Theory and evidence.  _Journal of Business & Economic Statistics_, 37(2):187–204, 2019.
  * Patel et al. [2015] C. J. Patel, B. Burford, and J. P. Ioannidis.  Assessment of vibration of effects due to model specification can demonstrate the instability of observational associations.  _Journal of clinical epidemiology_ , 68(9):1046–1058, 2015.
  * Pearl [2009] J. Pearl.  _Causality: Models, reasoning, and inference_.  Cambridge University Press, 2nd edition, 2009.
  * Peters et al. [2016] J. Peters, P. Bühlmann, and N. Meinshausen.  Causal inference by using invariant prediction: Identification and confidence intervals.  _Journal of the Royal Statistical Society, Series B_ , 78(5):947–1012, 2016.
  * Pfister et al. [2021] N. Pfister, E. G. Williams, J. Peters, R. Aebersold, and P. Bühlmann.  Stabilizing variable selection and regression.  _The Annals of Applied Statistics_ , 15(3):1220–1246, 2021.
  * Rojas-Carulla et al. [2018] M. Rojas-Carulla, B. Schölkopf, R. Turner, and J. Peters.  Causal transfer in machine learning.  _Journal of Machine Learning Research_ , 19(36):1–34, 2018.
  * Rosenbaum [1987] P. Rosenbaum.  The role of a second control group in an observational study.  _Statistical Science_ , 2(3):292–306, 1987.
  * Rosenbaum [2010] P. Rosenbaum.  Evidence factors in observational studies.  _Biometrika_ , 97(2):333–345, 2010.
  * Rosenbaum [2021] P. Rosenbaum.  _Replication and Evidence Factors in Observational Studies_.  CRC Press, 2021.
  * Rosenbaum and Rubin [1983] P. Rosenbaum and D. Rubin.  Assessing sensitivity to an unobserved binary covariate in an observational study with binary outcome.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 45(2):212–218, 1983.
  * Rothenhäusler et al. [2015] D. Rothenhäusler, C. Heinze, J. Peters, and N. Meinshausen.  Backshift: Learning causal cyclic graphs from unknown shift interventions.  In _Advances in Neural Information Processing Systems 29 (NIPS)_ , pages 1513–1521, 2015.
  * Rothenhäusler et al. [2021] D. Rothenhäusler, N. Meinshausen, P. Bühlmann, and J. Peters.  Anchor regression: Heterogeneous data meet causality.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 83(2):215–246, 2021.
  * Rubin [1974] D. Rubin.  Estimating causal effects of treatments in randomized and nonrandomized studies.  _Journal of Educational Psychology_ , 66(5):688, 1974.
  * Schölkopf et al. [2012] B. Schölkopf, D. Janzing, J. Peters, E. Sgouritsa, K. Zhang, and J. Mooij.  On causal and anticausal learning.  In _Proceedings of the 29th International Conference on Machine Learning (ICML)_ , pages 1255–1262, 2012.
  * Skene et al. [1986] A. Skene, J. Shaw, and T. Lee.  Bayesian modelling and sensitivity analysis.  _Journal of the Royal Statistical Society_ , 35(2):281–288, 1986.
  * Srivastava et al. [2014] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov.  Dropout: a simple way to prevent neural networks from overfitting.  _The Journal of Machine Learning Research_ , 15(1):1929–1958, 2014.
  * Steegen et al. [2016] S. Steegen, F. Tuerlinckx, A. Gelman, and W. Vanpaemel.  Increasing transparency through a multiverse analysis.  _Perspectives on Psychological Science_ , 11(5):702–712, 2016.
  * Tsiatis [2006] A. Tsiatis.  _Semiparametric theory and missing data_.  Springer, 2006.
  * van der Vaart [2000] A. van der Vaart.  _Asymptotic Statistics_ , volume 3.  Cambridge university press, 2000.
  * VanderWeele and Ding [2017] T. VanderWeele and P. Ding.  Sensitivity analysis in observational research: introducing the e-value.  _Annals of Internal Medicine_ , 167(4):268–274, 2017.
  * Wedderburn [1974] R. W. Wedderburn.  Quasi-likelihood functions, generalized linear models, and the gauss—newton method.  _Biometrika_ , 61(3):439–447, 1974.
  * Yu [2013] B. Yu.  Stability.  _Bernoulli_ , 19(4):1484–1500, 2013.
  * Yu and Kumbier [2020] B. Yu and K. Kumbier.  Veridical data science.  _Proceedings of the National Academy of Sciences_ , 117(8):3920–3929, 2020.
  * Zhang et al. [2013] K. Zhang, B. Schölkopf, K. Muandet, and Z. Wang.  Domain adaptation under target and conditional shift.  In _International Conference on Machine Learning_ , pages 819–827, 2013.

## Appendix

In Section A, we discuss additional properties of the isotropic perturbation
model. Section B contains the proofs. Section C discusses how to form robust
confidence intervals if the data analyst trusts one of the estimators
$\hat{\theta}^{k}$ more than others.

##  Appendix A Properties of the isotropic perturbation model

Recall that conditionally on $\xi$, the data $(D_{i})_{1\leq i\leq n}$ are
drawn i.i.d. from the perturbed distribution $\mathbb{P}^{\xi}(D=\bullet)$,
where $\xi$ is an unobserved random variable. Note that an estimator
$\hat{\theta}=\hat{\theta}(D_{1},\ldots,D_{n})$ for some parameter
$\theta^{0}(\mathbb{P})$ now has two sources of uncertainty: the uncertainty
due to sampling and the uncertainty due to the random perturbation.

$\hat{\theta}-\theta^{0}(\mathbb{P})=\underbrace{\hat{\theta}-\theta(\mathbb{P}^{\xi})}_{\begin{subarray}{c}\text{variation
due to}\\\
\text{sampling}\end{subarray}}+\underbrace{\theta(\mathbb{P}^{\xi})-\theta^{0}(\mathbb{P})}_{\begin{subarray}{c}\text{variation
due to}\\\ \text{random perturbation}\end{subarray}}$ (22)

We refer to the second component as _distributional uncertainty_. In this
section we will study such distributional perturbation models in more detail.
In Section A.1 we will show that under a symmetry assumption, there exists
only one class of perturbation models that is characterized by a one-
dimensional parameter $\delta_{\text{dist}}$. In Section A.2 we will sketch an
extension of the random perturbation model that allows different parts of the
distribution to be affected by different perturbations. In Section A.3 we
discuss how the perturbation model relates to other issues with data, such as
sampling bias and confounding.

###  A.1 Uniqueness of distributional perturbation model

The distributional perturbation models introduced earlier raise the question
whether other distributional perturbation models are reasonable. Let us go
back to the construction of a discrete perturbation model in Section 2. In the
example described in Section 2, $\mathbb{P}$ is the uniform probability
measure on $\\{1,\ldots,K\\}$. In equation (5), we constructed the perturbed
distribution by setting

$\mathbb{P}^{\xi}(D=k)=\xi_{k}.$ (23)

for an exchangeable set of random variables $\xi_{1},\ldots,\xi_{K}\geq 0$
with $\sum\xi_{k}=1$. Analogously as in Section B.1 it can be shown that for
any function $f:\\{1,\ldots,K\\}\rightarrow\mathbb{R}$,

$\text{Var}(\mathbb{E}^{\xi}[f(D)])=\delta_{\text{dist}}^{2}\text{Var}_{\mathbb{P}}(f(D)),$

where $\delta_{\text{dist}}^{2}=\frac{K^{2}}{K-1}\text{Var}(\xi_{1})$. Thus,
under this model, the variance of the perturbation is proportional to the
variance in the unperturbed distribution. This raises the question whether
this property is specific for the considered example or whether it holds in
more generality. We will see that under a symmetry assumption, there exists
only one type of perturbation model, which is equivalent to the one considered
above. Roughly speaking, the symmetry assumption states that two events that
have equal probability under $\mathbb{P}$ are perturbed in a similar fashion.
In the following, we write $\mathbb{Q}$ for the marginal distribution of
$(D,\xi)$, where first the perturbation $\xi$ is drawn and then
$D\sim\mathbb{P}^{\xi}$. The proof of the following result can be found in
Section B.6.

######  Theorem 2 (Characterization of isotropic perturbation models).

Let $(D,\xi)\sim\mathbb{Q}$ and assume that there exists a function
$h(\bullet)$ such that $h(D)$ is uniformly distributed on $[0,1]$. Let
$\mathbb{P}^{\xi}=\mathbb{Q}(\bullet|\xi)$ and let $\mathbb{P}$ denote the
marginal distribution of $D$ under $\mathbb{Q}$. Assume that for any
$D$-measurable events $A$ and $B$ with $\mathbb{P}(A)=\mathbb{P}(B)$,

$\mathrm{Var}({\mathbb{P}^{\xi}(A)})=\mathrm{Var}({\mathbb{P}^{\xi}(B)}).$
(24)

Furthermore, assume that for every sequence of $D$-measurable events $A_{j}$
with $\mathbb{P}(A_{j})\rightarrow 0$,

$\mathrm{Var}({\mathbb{P}^{\xi}(A_{j})})\rightarrow 0.$

Then for any $\phi(D)\in L^{2}(\mathbb{P})$

$\mathrm{Var}(\mathbb{E}^{\xi}[\phi(D)])=\delta_{\text{dist}}^{2}\mathrm{Var}_{\mathbb{P}}(\phi(D)),$
(25)

for some fixed $\delta_{\text{dist}}\geq 0$.

Let us discuss what this result means for the behaviour of empirical means.
Let $D_{1},\ldots,D_{n}$ be i.i.d. drawn from $\mathbb{P}^{\xi}$. Then, for
all square-integrable functions $\psi(D)\in L^{2}(\mathbb{P})$ marginally
across sampling uncertainty and distributional uncertainty we have

$\displaystyle\text{Var}_{P}(\frac{1}{n}\sum_{i=1}^{n}\psi(D_{i}))$
$\displaystyle=\,\,$
$\displaystyle(\frac{1}{n}+\delta_{\text{dist}}^{2}-\frac{\delta_{\text{dist}}^{2}}{n})\text{Var}_{\mathbb{P}}(\psi(D)).$
$\displaystyle=\,\,$
$\displaystyle\frac{\delta^{2}}{n}\text{Var}_{\mathbb{P}}(\psi(D)),$

with $\delta^{2}=1+n\delta_{\text{dist}}^{2}-\delta_{\text{dist}}^{2}$. Since
$\mathbb{P}^{\xi}=\mathbb{Q}(\bullet|\xi)$ we also have
$E[\frac{1}{n}\sum_{i=1}^{n}\psi(D_{i})]=\mathbb{E}[\psi(D)]$.

There are two major assumptions in this theorem. The first assumption says
that two events that have the same probability are perturbed in the same
fashion. This can be seen as a symmetry assumption. The second assumption says
that events that have a small probability are only perturbed by a small
amount. This can be seen as a regularity assumption.

Then, up to a one-dimensional parameter $\delta$, the variance of functions is
uniquely determined. This means that using symmetry we have reduced the
problem of estimating an infinite-dimensional perturbation model to a one-
dimensional quantity $\delta$. Note that the statement in Theorem 2 is
slightly weaker than Lemma 1, since it is only a statement about variances and
not about the asymptotic distribution of empirical processes.

The distributional uncertainty model relates to other concepts in statistics
such as sampling bias, correlated data, and confounding. We discuss these
connections in the Appendix, Section A.3.

In practice, some researchers might object to the symmetry assumption in
equation (24). It turns out that the perturbation model can be generalized. In
the following section, we will give a brief outlook of how perturbation models
can be used to perturb different parts of a distribution differently.

###  A.2 Beyond isotropic distributional perturbations

The discussion in Section A.1 shows that under a symmetry assumption, up to an
unknown scale factor $\delta$, there exists only one type of perturbation
model. However, in practice there might be a situation where one does not
expect a perturbation to affect all parts of the distribution in the same way.
Consider $D=(X,Y)$. For example, one might expect that the distribution of $X$
is perturbed between settings but that the measurement error is invariant.
This may lead one to want to model a situation where $p(x)$ is perturbed but
$p(y|x)$ is not perturbed. Under appropriate regularity conditions on $\psi$
we have

$\displaystyle\mathbb{E}^{\xi}[\psi(X,Y)]-\mathbb{E}[\psi(X,Y)]$
$\displaystyle=\mathbb{E}^{\xi}[\mathbb{E}[\psi(X,Y)|X]]-\mathbb{E}[\mathbb{E}[\psi(X,Y)|X]]$
(26) $\displaystyle\stackrel{{\scriptstyle
d}}{{\approx}}\mathcal{N}(0,\delta_{\text{dist}}^{2}\text{Var}_{\mathbb{P}}(\mathbb{E}[\psi(X,Y)|X])).$
(27)

If $\delta_{\text{dist}}$ is known or can be estimated, this allows us to
adjust variance and confidence intervals to account for uncertainty both due
to sampling and distributional perturbations, similarly as in Section 3.

###  A.3 Relationship to other concepts

In this section we will see that issues that are common in data analysis (such
as sampling bias, confounding, and violations of the i.i.d. assumptions) are
consequences of the isotropic perturbation model. Depending on how strong the
perturbations are, the issues with confounding and sampling bias are more or
less severe. To be clear, the goal in this section is to show that the concept
of distributional uncertainty is connected to sampling bias and confounding,
but neither of these concepts is subsumed by one of the other concepts.

In the following examples, we show that the random perturbation model implies
violations of the independence assumption, sampling bias, and confounding.

######  Example 2 (Violation of i.i.d. assumption).

Draw real-valued random variables
$D_{i}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\mathbb{P}^{\xi}$, where
$\mathbb{P}^{\xi}$ satisfies equation (25). Let $\sigma^{2}<\infty$ denote the
variance of $D_{1}$ under $P$, where $P$ is the marginal distribution of
$(D_{1},D_{2},...,D_{n},\xi)$ where first $\xi$ is drawn and then
$D_{i}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\mathbb{P}^{\xi}$. For all
$i\neq j$, and $\phi(D)\in L^{2}(\mathbb{P})$,

$\displaystyle\mathrm{Cov}(\phi(D_{i}),\phi(D_{j}))$
$\displaystyle=E[\phi(D_{i})\phi(D_{j})]-E[\phi(D_{i})]E[\phi(D_{j})]$ (28)
$\displaystyle=E[\mathbb{E}^{\xi}[\phi(D)]\mathbb{E}^{\xi}[\phi(D)]]-E[\mathbb{E}^{\xi}[\phi(D)]]E[\mathbb{E}^{\xi}[\phi(D)]]$
(29)
$\displaystyle=E[\mathbb{E}^{\xi}[\phi(D)]^{2}]-E[\mathbb{E}^{\xi}[\phi(D)]]^{2}$
(30) $\displaystyle=\mathrm{Var}(\mathbb{E}^{\xi}[\phi(D)])$ (31)
$\displaystyle=\delta_{\text{dist}}^{2}\sigma^{2}.$ (32)

Here, we used that $E[\phi(D_{i})|\xi]=E[\phi(D_{j})|\xi]$ and equation (25).
Thus, marginally, the perturbation model induces equicorrelated data. The
challenge at this point is that $\delta_{\text{dist}}$ is unknown, which means
that the correlation of the covariates is unknown. For a fixed perturbation
$\xi$, even if $n\rightarrow\infty$ it is impossible to distinguish between
i.i.d. data drawn from $\mathbb{P}^{\xi}=\mathbb{P}$ (corresponding to
$\delta_{\text{dist}}=0$) or data drawn from a distributional perturbation
model with $\delta_{\text{dist}}>0$.

######  Example 3 (Sampling bias).

Let $O\in\\{0,1\\}$ denote the indicator whether a unit is observed. We assume
that under the target population there is no sampling bias, that means
$\mathbb{P}(X=\bullet|O=1)=\mathbb{P}(X=\bullet)$. Now assume that
$\mathbb{P}^{\xi}$ satisfies equation (25) for $D=(X,O)$. In the following, we
will see that the variance of
$\mathbb{P}^{\xi}(X=\bullet|O=1)-\mathbb{P}(X=\bullet)$ is nonzero, i.e. the
random perturbation will induce sampling bias with positive probability. For
any $D$-measurable set $A$, using a Taylor expansion,

$\displaystyle\mathbb{P}^{\xi}(X\in A|O=1)-\mathbb{P}(X\in A)$ (33)
$\displaystyle\approx\mathbb{P}(X\in A)+\frac{\mathbb{P}^{\xi}(O=1,X\in
A)-\mathbb{P}(O=1,X\in A)}{\mathbb{P}(O=1)}-\mathbb{P}(X\in
A)\frac{\mathbb{P}^{\xi}(O=1)-\mathbb{P}(O=1)}{\mathbb{P}(O=1)}-\mathbb{P}(X\in
A)$ (34)

It can be shown that this term has variance

$\delta_{\text{dist}}^{2}\frac{\mathbb{P}[X\in A](1-\mathbb{P}[X\in
A])}{\mathbb{P}[O=1]}.$ (35)

Thus, if $\mathbb{P}(X\in A)\not\in\\{0,1\\}$, the variance of
$\mathbb{P}^{\xi}(X\in A|O=1)-\mathbb{P}(X\in A)$ is nonzero. Hence, if we
draw an i.i.d. sample $X=(X_{1},\ldots,X_{n})$ from $\mathbb{P}^{\xi}$ and
only observe units with $O=1$, some members of the sampling population
$\mathbb{P}^{\xi}(\bullet|O=1)$ have lower or higher sampling probability
compared to the target population $\mathbb{P}$. In this sense, the random
perturbation model induces sampling bias.

######  Example 4 (Confounding).

We consider the Neyman-Rubin model [Neyman, 1923, Rubin, 1974] where $Y(t)$
denotes the potential outcome under treatment assignment $T=t\in\\{0,1\\}$ and
$X$ are observed covariates. The statistician observes $Y=TY(1)+(1-T)Y(0)$. In
observational studies, it is common to make the unconfoundedness assumption.
Under unconfoundedness, under $\mathbb{P}$ we have the following conditional
independence:

$\\{Y(1),Y(0)\\}\perp T|X.$ (36)

In the following, for simplicity we assume that $X$ is discrete. If we have
i.i.d. observations drawn from $\mathbb{P}$, under a positivity assumption we
can estimate $\mathbb{E}[Y(1)-Y(0)]$ using inverse probability weighting,
regression adjustment, and other methods [Hernán and Robins, 2020]. In the
distributional uncertainty model, we do not draw samples from $\mathbb{P}$ but
from $\mathbb{P}^{\xi}$, where $\mathbb{P}^{\xi}$ randomly deviates from
$\mathbb{P}$. We will now show that under $\mathbb{P}^{\xi}$, with positive
probability the unconfoundedness assumption is violated.

Ignorability implies that for any function $g$

$\mathbb{E}[(\frac{T}{e(X)}-1)Y(t)g(X)]=0,$ (37)

where $e(x)=\mathbb{P}[T=1|X=x]$. We will show that

$\mathbb{E}^{\xi}[(\frac{T}{e(X)}-1)Y(t)g(X)]\neq 0$ (38)

with positive probability. If equation (38) holds, then with positive
probability under $\mathbb{P}^{\xi}$

$Y(t)\not\perp T|X.$ (39)

Let us now show that with positive probability, equation (38) holds. For
simplicity, set $g(x)=1_{X=x}$. Choose $x$ such that
$\mathbb{P}[X=x]\in(0,1)$, $e(x)\in(0,1)$ and $t$ such that
$\mathbb{E}[Y(t)|X=x]\not\equiv 0$. Then, using equation (25)

$\mathbb{E}^{\xi}[(\frac{T}{e(X)}-1)Y(t)g(X)]-\mathbb{E}[(\frac{T}{e(X)}-1)Y(t)g(X)]$
(40)

has variance

$\delta_{\text{dist}}^{2}\text{Var}_{\mathbb{P}}\left((\frac{T}{e(X)}-1)Y(t)g(X)\right).$
(41)

Under the stated assumptions this variance must be positive. Thus, with
positive probability, equation (38) holds.

Thus, ignorability is violated with positive probability for
$\mathbb{P}^{\xi}$ even if it is satisfied for the unperturbed distribution
$\mathbb{P}$. Note that this model does not imply that there is confounding
for a certain data set (that means, for a certain sampling distribution
$\mathbb{P}^{\xi}$). This model implies that overall, across data sets, some
data sets will be drawn from a distribution $\mathbb{P}^{\xi}$ for which
unconfoundedness is violated. The strength of the confounding is described by
the (a priori unknown) parameter $\delta_{\text{dist}}$. If we want to
construct confidence intervals that cover $\theta^{0}=\mathbb{E}[Y(1)-Y(0)]$
we have to take into account both the uncertainty due to sampling and the
uncertainty due to the distributional perturbation.

##  Appendix B Proofs

###  B.1 Variance calculation for equation (6)

Recall that $D_{1},\ldots,D_{n}$ are drawn i.i.d. from $\mathbb{P}^{\xi}$,
with

$\mathbb{P}^{\xi}[D=k]=\xi_{k},\text{ for $k=1,\ldots,K$, }$

where $\xi_{k}\geq 0$ are exchangeable random variables with $\sum\xi_{k}=1$,
and $s(k)=\sum_{i}1_{D_{i}=k}$. Furthermore, let
$\overline{\psi}=\frac{1}{K}\sum_{k=1}^{K}\psi(k)$. Then,

$\displaystyle\text{Var}_{P}(\frac{1}{n}\sum_{i=1}^{n}\psi(D_{i}))$
$\displaystyle=\text{Var}_{P}(\sum_{k=1}^{K}\frac{s(k)}{n}\psi(k))$
$\displaystyle=\frac{1}{n^{2}}\sum_{k\neq
k^{\prime}}E[s(k)s(k^{\prime})](\psi(k)-\overline{\psi})(\psi(k^{\prime})-\overline{\psi})+\frac{1}{n^{2}}\sum_{k}E[s(k)^{2}](\psi(k)-\overline{\psi})^{2}$
$\displaystyle=\frac{K}{n^{2}}E[s(1)^{2}-s(1)s(2)]\frac{1}{K}\sum_{k}(\psi(k)-\overline{\psi})^{2}$
$\displaystyle=\frac{K}{n^{2}}E[s(1)^{2}-s(1)s(2)]\text{Var}_{\mathbb{P}}(\psi(D))$
$\displaystyle=\frac{K}{n^{2}}E[\frac{K}{K-1}s(1)^{2}-\frac{n^{2}K}{K-1}\frac{1}{K^{2}}]\text{Var}_{\mathbb{P}}(\psi(D))$
$\displaystyle=\frac{K^{2}}{n^{2}(K-1)}\text{Var}_{P}(s(1))\text{Var}_{\mathbb{P}}(\psi(D))$

Here, we used repeatedly that $\sum_{k}s(k)=n$.

###  B.2 Auxiliary results and proof of Lemma 1

Notation: We write $\mathbb{P}$ for the target distribution and
$\mathbb{P}^{\xi}$ for the randomly perturbed distribution from which we draw
$n$ i.i.d. data samples $(D_{i})_{i=1,\ldots,n}$. In both examples $\xi$ can
be seen as a random variable that encodes the perturbations. The expectation
of $f(D_{1},\dots,D_{n})$ over the joint distribution of
$(\xi,D_{1},\dots,D_{n})$ can be written as
$E_{\xi}[\mathbb{E}^{\xi}(f(D_{1},\dots,D_{n})]$ where $E_{\xi}$ means we take
the expectation over $\xi$ and $\mathbb{E}^{\xi}$ means that we take the
expectation over $(D_{1},\ldots,D_{n})$, conditionally on $\xi$.

####  B.2.1 Auxiliary results

Let us first state an auxiliary lemma that will turn out helpful for proving
Lemma 1.

######  Lemma 4.

Let the assumptions of Lemma 1 hold. For the sequence of random variables
$\xi=\xi(n)$, for any bounded $\psi(\bullet)$ we have that

$\mathbb{E}^{\xi}[\psi(D)]-\mathbb{E}[\psi(D)]\stackrel{{\scriptstyle
d}}{{=}}\mathcal{N}(0,\gamma_{n}^{2}\text{Var}(\psi(D)))+o_{p}(\gamma_{n})$
(42)

where $\gamma_{n}^{2}=\frac{\text{Var}(W)}{m(n)E[W]^{2}}$ and where we write
$m(n)$ to make it explicit that $m$ grows with $n$.

###### Proof.

Let $\phi=\psi\circ h$. Without loss of generality, assume that
$\mathbb{E}[\phi(U)]=0.$ Note that

$\sqrt{m}(\mathbb{E}^{\xi}[\phi(U)]-\mathbb{E}[\phi(U)])=\frac{\sqrt{m}\sum_{k=1}^{m}\int_{x\in
I_{k}}\phi(x)dx\cdot(W_{k}-E[W])}{\sum_{k=1}^{m}W_{k}/m}.$

Let

$Y_{m,k}:=\sqrt{m}\int_{x\in I_{k}}\phi(x)dx\cdot(W_{k}-E[W]).$

First, note that

$E[Y_{m,k}]=0$ (43)

for all $k$. As the second step, we want to show that

$\sum_{k=1}^{m}E[Y_{m,k}^{2}]=\text{Var}(W)\cdot
m\sum_{k=1}^{m}\left(\int_{x\in
I_{k}}\phi(x)dx\right)^{2}\xrightarrow[]{}\text{Var}(W)\cdot\text{Var}_{\mathbb{P}}(\phi(U)).$
(44)

For any $f\in L^{2}([0,1])$, define $\Pi_{m}(f)$ as

$\Pi_{m}(f)(x)=\sum_{k=1}^{m}\left(m\int_{x\in I_{k}}f(x)dx\right)\cdot I(x\in
I_{k}).$

Then, we have

$\displaystyle\Bigg{|}m\sum_{k=1}^{m}\left(\int_{x\in
I_{k}}\phi(x)dx\right)^{2}-\text{Var}_{\mathbb{P}}(\phi(U))\Bigg{|}$
$\displaystyle=||\phi-\Pi_{m}(\phi)||_{2}^{2}\xrightarrow[]{}0.$

as $m$ goes to infinity. This is because any bounded function can be
approximated by a sequence of step functions of the form
$\sum_{k=1}^{m}b_{k}\chi_{I_{k}}$. Next we will show that for any
$\epsilon>0$,

$g_{m}(\epsilon)=\sum_{k=1}^{m}E[Y_{m,k}^{2};|Y_{m,k}|\geq\epsilon]\xrightarrow[]{}0.$
(45)

This is implied by the dominated convergence theorem as

$\displaystyle\sum_{k=1}^{m}E[Y_{m,k}^{2};|Y_{m,k}|\geq\epsilon]$
$\displaystyle\leq\sum_{k=1}^{m}\left(\int_{x\in
I_{k}}\phi^{2}(x)dx\right)E[(W_{k}-E[W])^{2}I(||\phi||_{\infty}|W_{k}-E[W]|/\sqrt{m}\geq\epsilon)]$
(46)
$\displaystyle=||\phi||_{2}^{2}E[(W-E[W])^{2}I(||\phi||_{\infty}|W-E[W]|/\sqrt{m}\geq\epsilon)]\xrightarrow[]{}0.$
(47)

Combining equations (43), (44), and (45), we can apply Lindeberg’s CLT. With
Slutsky’s theorem, we have

$\displaystyle\sqrt{m}(\mathbb{E}^{\xi}[\phi(U)]-\mathbb{E}[\phi(U)])=\frac{\sum_{k=1}^{m}Y_{m,k}}{\sum_{k=1}^{m}W_{k}/m}\stackrel{{\scriptstyle
d}}{{=}}\mathcal{N}(0,\text{Var}(W)\text{Var}_{\mathbb{P}}(\phi(U))/E[W]^{2})+o_{p}(1).$

This completes the proof. ∎

######  Lemma 5.

Let the assumptions of Lemma 1 hold. Assume that for a sequence of random
variables $\xi=\xi(n)$ there exists a sequence $\gamma_{n}$ with limit
$\delta^{2}=\lim_{n}(1+n\gamma_{n}^{2})<\infty$ such that for any bounded
$\psi(\bullet)$ we have

$\mathbb{E}^{\xi}[\psi(D)]-\mathbb{E}[\psi(D)]\stackrel{{\scriptstyle
d}}{{=}}\mathcal{N}(0,\gamma_{n}^{2}\text{Var}_{\mathbb{P}}(\psi(D)))+o_{p}(\gamma_{n}).$
(48)

Then, for any bounded $\psi(\bullet)$, it holds that

$\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\psi(D_{i}^{n})-\mathbb{E}[\psi(D)]\stackrel{{\scriptstyle
d}}{{=}}\mathcal{N}(0,\delta^{2}\text{Var}_{\mathbb{P}}(\psi(D)))+o_{p}(1).$
(49)

###### Proof.

In the proof, we suppress the dependence of $\xi$ on $n$. We want to show that
for any $x$,

$E_{\xi}\Bigg{[}\mathbb{P}^{\xi}\left(\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(\psi(D_{i}^{n})-\mathbb{E}[\psi(D)])\leq
x\cdot\sqrt{\delta^{2}\text{Var}_{\mathbb{P}}(\psi(D))}\right)\Bigg{]}=\Phi(x)+o(1),$
(50)

where $\Phi$ is the cdf of a standard Gaussian random variable. Let us define

$Y_{n}=x\cdot\delta\cdot\frac{\sqrt{\text{Var}_{\mathbb{P}}(\psi(D))}}{\sqrt{\text{Var}_{\mathbb{P}^{\xi}}(\psi(D))}}-\frac{\sqrt{n}(\mathbb{E}^{\xi}[\psi(D)]-\mathbb{E}[\psi(D)])}{\sqrt{\text{Var}_{\mathbb{P}^{\xi}}(\psi(D))}},$
(51)

where $\text{Var}_{\mathbb{P}^{\xi}}(\psi(D))$ denotes the variance of
$\psi(D)$ where $D\sim\mathbb{P}^{\xi}$. Then,

$\displaystyle(\ast)$
$\displaystyle=E_{\xi}\Bigg{[}\mathbb{P}^{\xi}\left(\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(\psi(D_{i}^{n})-\mathbb{E}^{\xi}[\psi(D)])+\sqrt{n}(\mathbb{E}^{\xi}[\psi(D)]-\mathbb{E}[\psi(D)])\leq
x\cdot\sqrt{\delta^{2}\text{Var}_{\mathbb{P}}(\psi(D))}\right)\Bigg{]}$ (52)
$\displaystyle=E_{\xi}\Bigg{[}\mathbb{P}^{\xi}\left(\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\frac{\psi(D_{i}^{n})-\mathbb{E}^{\xi}[\psi(D)]}{\sqrt{\text{Var}_{\mathbb{P}^{\xi}}(\psi(D))}}\leq
Y_{n}\right)\Bigg{]}.$ (53)

We define $g_{n}(y;\xi)$ as

$g_{n}(y;\xi)=\mathbb{P}^{\xi}\left(\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\frac{\psi(D_{i}^{n})-\mathbb{E}^{\xi}[\psi(D)]}{\sqrt{\text{Var}_{\mathbb{P}^{\xi}}(\psi(D))}}\leq
y\right).$ (54)

By Berry–Esseen, it holds that

$\sup_{y}\Bigg{|}\mathbb{P}^{\xi}\left(\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\frac{\psi(D_{i})-\mathbb{E}^{\xi}[\psi(D)]}{\sqrt{\text{Var}_{\mathbb{P}^{\xi}}(\psi(D))}}\leq
y\right)-\Phi(y)\Bigg{|}\leq\frac{C\mathbb{E}^{\xi}|{\psi(D)}^{3}|}{(\mathbb{E}^{\xi}|{\psi(D)}^{2}|)^{3/2}\sqrt{n}},$
(55)

for all $n$. Invoking equation (48) for $\psi(D)^{2}$ and $\psi(D)^{3}$, we
have that
$\mathbb{E}^{\xi}|{\psi(D)}^{3}|/(\mathbb{E}^{\xi}|{\psi(D)}^{2}|)^{3/2}$
converges in probability to
$\mathbb{E}|{\psi(D)}^{3}|/(\mathbb{E}|{\psi(D)}^{2}|)^{3/2}<\infty$ as
$n\rightarrow\infty$. Then the right-hand side of the above inequality
converges in probability to 0 as $n\xrightarrow[]{}\infty$, which implies that

$\sup_{y}|g_{n}(y;\xi)-\Phi(y)|\xrightarrow[]{p}0.$ (56)

Using this result,

$\displaystyle(\ast)=E_{\xi}[g_{n}(Y_{n})]$
$\displaystyle=E_{\xi}[g_{n}(Y_{n})]-E_{\xi}[\Phi(Y_{n})]+E_{\xi}[\Phi(Y_{n})]$
(57) $\displaystyle\leq
E_{\xi}[\sup_{y}|g_{n}(y)-\Phi(y)|]+E_{\xi}[\Phi(Y_{n})]$ (58)
$\displaystyle=E_{\xi}[\Phi(Y_{n})]+o(1).$ (59)

by the dominated convergence theorem. Using equation (48),
$\text{Var}_{\mathbb{P}^{\xi}}(\psi(D))\xrightarrow[]{p}\text{Var}_{\mathbb{P}}(\psi(D))$.
Then, we have

$Y_{n}\xrightarrow[]{d}\delta x-\sqrt{\delta^{2}-1}Z,$

where $Z$ is a standard Gaussian random variable. Since $\Phi$ is bounded and
continuous, by Portmanteau Lemma, we get

$\lim_{n\xrightarrow[]{}\infty}E_{\xi}[\Phi(Y_{n})]=E[\Phi(\delta
x-\sqrt{\delta^{2}-1}Z)]=\Phi(x).$ (60)

This completes the proof. ∎

####  B.2.2 Proof of Lemma 1

Now let us show that the Lemma 1 holds.

###### Proof.

Without loss of generality for notational simplicity we restrict ourselves to
the case $l=1$, i.e. $\psi:\mathcal{D}\mapsto\mathbb{R}$. As before we write
$\phi(D)=\psi\circ h(U)$. For any $\psi\in L^{2}(\mathbb{P})$ and for any
$\epsilon>0$, there exits a bounded function $\psi^{B}$ such that
$\mathbb{E}[\psi(D)]=\mathbb{E}[\psi^{B}(D)]$ and
$||\psi-\psi^{B}||_{L^{2}}<\epsilon$. Note that

$\displaystyle\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(\psi(D_{i}^{n})-\mathbb{E}[\psi(D)])$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(\psi(D_{i}^{n})-\psi^{B}(D_{i}))$
(a)
$\displaystyle+\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(\psi^{B}(D_{i}^{n})-\mathbb{E}[\psi^{B}(D)])$
(b)

Without loss of generality, let’s assume that $\mathbb{E}[\psi(D)]=0$. Note
that

$\displaystyle\eqref{eq:eq1}=\underbrace{\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left((\psi-\psi^{B})(D_{i}^{n})-\mathbb{E}^{\xi}[(\psi-\psi^{B})(D)]\right)}_{(a.1)}+\underbrace{\sqrt{n}\mathbb{E}^{\xi}[(\psi-\psi^{B})(D)]}_{(a.2)}.$

The marginal variance of (a.1) is

$\displaystyle
E[\text{Var}(\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\left((\psi-\psi^{B})(D_{i}^{n})-\mathbb{E}^{\xi}[(\psi-\psi^{B})(D)]\right)|\xi)]$
$\displaystyle\leq E_{\xi}[\mathbb{E}^{\xi}[(\psi-\psi^{B})^{2}(D)]]$
$\displaystyle=\mathbb{E}[(\psi-\psi^{B})^{2}(D)].$

Therefore, we can write (a.1) = $\epsilon\cdot s_{n}$ with $E[s_{n}]=0$ and
$\text{Var}(s_{n})\leq 1$. Recall that we write $\phi(D)=\psi\circ h(U)$. Note
that for any $\phi\in L^{2}([0,1])$ such that $\mathbb{E}[\phi(U)]=0$,

$\sqrt{m}(\mathbb{E}^{\xi}[\phi(U)])=\frac{\sqrt{m}\sum_{k=1}^{m}\int_{x\in
I_{k}}\phi(x)dx\cdot(W_{k}-E[W])}{\sum_{k=1}^{m}W_{k}/m}.$

The denominator converges in probability to $E[W]$. The numerator has mean $0$
and its variance is

$\displaystyle\text{Var}(W)\sum_{k=1}^{m}m\left(\int_{x\in
I_{k}}\phi(x)dx\right)^{2}$
$\displaystyle\leq\text{Var}(W)\sum_{k=1}^{m}\int_{x\in I_{k}}\phi^{2}(x)dx$
$\displaystyle=\text{Var}(W)\mathbb{E}[\phi^{2}(U)]$

where the first inequality holds by Jensen’s inequality with $m\int_{x\in
I_{k}}dx=1$. Therefore, using the argument above with
$\phi=(\psi-\psi^{B})\circ h$, we get that

$\displaystyle(a.2)=\frac{\epsilon\cdot s_{n}^{\prime}}{E[W]+o_{p}(1)}.$

with $E[s_{n}^{\prime}]=0$ and $\text{Var}(s_{n}^{\prime})\leq\text{Var}(W)$.
Combining the results, we have (a)$=O_{p}(\epsilon)$. Moreover, with Lemma 4
and Lemma 5, we have that

$\displaystyle\eqref{eq:eq2}$ $\displaystyle\stackrel{{\scriptstyle
d}}{{=}}\mathcal{N}(0,\delta^{2}\text{Var}_{\mathbb{P}}(\psi^{B}(D)))+o_{p}(1)$
$\displaystyle\stackrel{{\scriptstyle
d}}{{=}}\mathcal{N}(0,\delta^{2}\text{Var}_{\mathbb{P}}(\psi(D)))+s_{n}^{\prime\prime}+o_{p}(1).$

where $\delta^{2}=1+r\text{Var}(W)/E[W]^{2}$ and
$s_{n}^{\prime\prime}\sim\mathcal{N}(0,\text{Var}(\psi(D)-\psi^{B}(D))=O_{p}(\epsilon)$.
Note that results hold for arbitrary $\epsilon>0$. Therefore, for any $\psi\in
L^{2}(\mathbb{P})$,

$\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\psi(D_{i}^{n})-\mathbb{E}[\psi(D)]\stackrel{{\scriptstyle
d}}{{=}}\mathcal{N}(0,\delta^{2}\text{Var}_{\mathbb{P}}(\psi(D)))+o_{p}(1).$

∎

###  B.3 Examples of variance inflation induced by non-i.i.d. sampling

##### Non-i.i.d. sampling.

In the following examples we discuss how Assumption 1 with $\delta\neq 1$
arises in non-standard sampling settings. For simplicity, we start with an
artificial example: sampling with replacement from an unknown subpopulation.

######  Example 5 (Sampling with replacement from an unknown subpopulation).

Assume that $D_{1}^{\prime},\dots,D_{m}^{\prime}$ drawn i.i.d. from
$\mathbb{P}$. Set $\xi=(D_{1}^{\prime},\ldots,D_{m}^{\prime})$. We define the
randomly perturbed distribution $\mathbb{P}^{\xi}$ as the empirical measure

$\mathbb{P}^{\xi}(D\in\bullet)=\frac{1}{m}\sum_{i=1}^{m}1_{D_{i}^{\prime}\in\bullet}.$

Let $n\rightarrow\infty$ and assume that $m(n)$ is a sequence of integers such
that $\frac{n}{m(n)}$ converges to some limit $r\in(0,\infty)$. Conditionally
on $\xi$, let $(D_{1}^{n},\ldots,D_{n}^{n})$ be i.i.d. draws from
$\mathbb{P}^{\xi}$. Then equation (7) holds for any $\psi(\bullet)$ with
finite second moment with

$\delta^{2}=1+r.$

###### Proof.

Suppose that $D_{1}^{\prime},\dots,D_{m}^{\prime}$ are drawn from $\mathbb{P}$
for some sequence $m=m(n)$. Let $\mathbb{P}^{\xi}$ denote the empirical
measure of $D_{1}^{\prime},\dots,D_{m}^{\prime}$. Then by the CLT, for any
$\psi(\bullet)$ with finite second moment,

$\displaystyle\mathbb{E}^{\xi}[\psi(D)]-\mathbb{E}[\psi(D)]$
$\displaystyle=\frac{1}{m}\sum_{i=1}^{m}\psi(D_{i}^{\prime})-\mathbb{E}[\psi(D)]$
(61) $\displaystyle\stackrel{{\scriptstyle
d}}{{=}}\mathcal{N}(0,\gamma_{n}^{2}\text{Var}_{\mathbb{P}}(\psi(D)))+o_{p}(\gamma_{n})$
(62)

where $\gamma_{n}^{2}=1/m(n)$. By applying Lemma 5 and following the proof of
Lemma 1, for any $\psi(\bullet)$ with finite second moment, we get

$\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\psi(D_{i}^{n})-\mathbb{E}[\psi(D)]\stackrel{{\scriptstyle
d}}{{=}}\mathcal{N}(0,\delta^{2}\text{Var}_{\mathbb{P}}(\psi(D)))+o_{p}(1)$

where $\delta^{2}=1+r$. ∎

Sampling with replacement from a finite population might seem a bit
artificial. The next example shows that a similar conclusion holds if we
sample clusters, where units in a single cluster are highly correlated, and
units between clusters are independent. If the cluster structure is known, one
can use clustered standard errors. However, in general the dependence
structure might be unknown.

######  Example 6 (Sampling clusters with unobserved membership).

Here, we consider a setting where some observations are associated, but where
the overall dependence structure is unknown. This is similar to the previous
setting, but there are no ties in the data set. Consider $\mathbb{P}$ a
probability distribution with positive density over a compact subset of
$\mathbb{R}^{p}$. Draw i.i.d. observations
$D_{1}^{\prime},\ldots,D_{m}^{\prime}$ from $\mathbb{P}$. Set
$\xi=(D_{1}^{\prime},\ldots,D_{m}^{\prime})$. Conditionally on $\xi$, let
$(D_{1}^{n},\ldots,D_{n}^{n})$ be i.i.d. draws from $\mathbb{P}^{\xi}$, where

$\mathbb{P}^{\xi}(D\in\bullet)=\frac{1}{m}\sum_{i=1}^{m}\mathbb{P}(D\in\bullet|\|D_{i}^{\prime}-D\|_{2}\leq\epsilon_{n}),$

where $\epsilon_{n}>0$ is a deterministic sequence with
$\epsilon_{n}=o(1/\sqrt{n})$. Furthermore, let $n\rightarrow\infty$ and assume
that $m=m(n)$ is a sequence of integers such that $\frac{n}{m(n)}$ converges
to some limit $r\in(0,\infty)$. Then, equation (7) holds for any bounded
Lipschitz continuous $\psi(\bullet)$ with

$\delta^{2}=1+r.$

###### Proof.

Using Lipschitz continuity and $\epsilon_{n}=o(1/\sqrt{n})$, we have

$\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\psi(D_{i}^{n})-\mathbb{E}[\psi(D)]\stackrel{{\scriptstyle
d}}{{=}}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\psi(D_{i}^{\prime\prime})-\mathbb{E}[\psi(D)]+o_{p}(1),$

where the $D_{i}^{\prime\prime}$ are drawn with replacement from
$D_{1}^{\prime},\ldots,D_{m}^{\prime}$. We can now invoke Example 5. ∎

###  B.4 Proof of Proposition 1

###### Proof.

Note that

$\displaystyle\widehat{\text{Var}_{\mathbb{P}}(\phi^{k}(D))}$
$\displaystyle=\underbrace{\frac{1}{n}\sum_{i=1}^{n}\left(\hat{\phi}^{k}(D_{i})-\phi^{k}(D_{i})-\frac{1}{n}\sum_{i=1}^{n}\hat{\phi}^{k}(D_{i})+\frac{1}{n}\sum_{i=1}^{n}\phi^{k}(D_{i})\right)^{2}}_{(i)}$
(63)
$\displaystyle+\underbrace{\frac{2}{n}\sum_{i=1}^{n}\left(\hat{\phi}^{k}(D_{i})-\phi^{k}(D_{i})-\frac{1}{n}\sum_{i=1}^{n}\hat{\phi}^{k}(D_{i})+\frac{1}{n}\sum_{i=1}^{n}\phi^{k}(D_{i})\right)\left(\phi^{k}(D_{i})-\frac{1}{n}\sum_{i=1}^{n}\phi^{k}(D_{i})\right)}_{(ii)}$
(64)
$\displaystyle+\underbrace{\frac{1}{n}\sum_{i=1}^{n}\left(\phi^{k}(D_{i})-\frac{1}{n}\sum_{i=1}^{n}\phi^{k}(D_{i})\right)^{2}}_{(iii)}.$
(65)

As the $\phi^{k}$ has finite fourth moments, we can use Lemma 1 to obtain
(iii) $=\text{Var}_{\mathbb{P}}(\phi^{k}(D))+o_{p}(1)$. Then by Cauchy-
Schwartz inequality and Jensen’s inequality,

(i)
$\displaystyle\leq\frac{2}{n}\sum_{i=1}^{n}\left(\hat{\phi}^{k}(D_{i})-\phi^{k}(D_{i})\right)^{2}+2\left(\frac{1}{n}\sum_{i=1}^{n}\left(\hat{\phi}^{k}(D_{i})-\phi^{k}(D_{i})\right)\right)^{2}$
$\displaystyle\leq\frac{4}{n}\sum_{i=1}^{n}\left(\hat{\phi}^{k}(D_{i})-\phi^{k}(D_{i})\right)^{2}.$

Since our influence function estimators are consistent, (i) $=o_{p}(1)$. Then
again by Cauchy-Schwartz inequality, (ii) $=o_{p}(1)$. Combining results, we
get

$\widehat{\text{Var}_{\mathbb{P}}(\phi^{k}(D))}=\text{Var}_{\mathbb{P}}(\phi^{k}(D))+o_{p}(1).$

This completes the proof. ∎

###  B.5 Proof of Theorem 1

###### Proof.

By Assumption 1, 2 and Lemma 1,

$\begin{pmatrix}\sqrt{n}(\hat{\theta}^{1}-\theta^{0})\\\ \vdots\\\
\sqrt{n}(\hat{\theta}^{K}-\theta^{0})\end{pmatrix}=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\begin{pmatrix}\phi^{1}(D_{i})\\\
\vdots\\\ \phi^{K}(D_{i})\end{pmatrix}+o_{p}(1)\stackrel{{\scriptstyle
d}}{{=}}\delta\textbf{Z}+o_{p}(1)$ (66)

where
$\textbf{Z}=(Z_{1},\dots,Z_{K})^{\intercal}\sim\mathcal{N}(\mathbf{0},\text{diag}(\text{Var}(\phi^{1}),\dots,\text{Var}(\phi^{K})))$.
As $n\xrightarrow[]{}\infty$, using that $\sum\alpha_{k}=1$ and
$\hat{\alpha}_{k}=\alpha_{k}+o_{p}(1)$,

$\sqrt{n}(\hat{\theta}^{W}-\theta^{0})=\sqrt{n}\sum_{k=1}^{K}\hat{\alpha}_{k}(\hat{\theta}^{k}-\theta^{0})+o_{p}(1)\stackrel{{\scriptstyle
d}}{{=}}\delta\sum_{k=1}^{K}\alpha_{k}Z_{k}+o_{p}(1)\stackrel{{\scriptstyle
d}}{{=}}\delta\mathcal{N}(0,\alpha)+o_{p}(1),$ (67)

where
$\alpha=\frac{1}{\sum_{k=1}^{K}\frac{1}{\text{Var}_{\mathbb{P}}(\phi^{k}(D))}}$.
By a similar calculation, we have that

$n\hat{\sigma}^{2}_{bet}=\delta^{2}\sum_{k=1}^{K}\alpha_{k}(Z_{k}-\sum_{j=1}^{K}\alpha_{j}Z_{j})^{2}+o_{p}(1).$
(68)

Thus,

$n\hat{\sigma}^{2}_{bet}\stackrel{{\scriptstyle d}}{{=}}\delta^{2}\alpha
L_{K}+o_{p}(1),$ (69)

where

$L_{K}=\frac{1}{\alpha}\sum_{k=1}^{K}\alpha_{k}(Z_{k}-\sum_{j=1}^{K}\alpha_{j}Z_{j})^{2}.$
(70)

We will now show that

$L_{K}\sim\chi^{2}(K-1)\quad\perp\quad\frac{1}{\sqrt{\alpha}}\sum_{j=1}^{K}\alpha_{j}Z_{j}\sim\mathcal{N}(0,1).$
(71)

First, note that

$L_{K}=\sum_{k=1}^{K}\left(Z_{k}\frac{\sqrt{\alpha_{k}}}{\sqrt{\alpha}}-\frac{\sqrt{\alpha_{k}}}{\sqrt{\alpha}}\sum_{j=1}^{K}\alpha_{j}Z_{j}\right)^{2}$
(72)

With this definition,

$\displaystyle L_{K}$
$\displaystyle=\sum_{k=1}^{K}\left(\frac{\sqrt{\alpha_{k}}}{\sqrt{\alpha}}Z_{k}-\frac{\sqrt{\alpha_{k}}}{\sqrt{\alpha}}\sum_{j=1}^{K}\alpha_{j}Z_{j}\right)^{2}$
$\displaystyle=\sum_{k=1}^{K}\left(\frac{\sqrt{\alpha_{k}}}{\sqrt{\alpha}}Z_{k}-\sqrt{\alpha_{k}}\sum_{j=1}^{K}\sqrt{\alpha_{j}}\frac{\sqrt{\alpha_{j}}}{\sqrt{\alpha}}Z_{j}\right)^{2}$
$\displaystyle=\sum_{k=1}^{K}(\tilde{Z}_{k}-w_{k}\sum_{j=1}^{K}w_{j}\tilde{Z}_{j})^{2}$

where
$\tilde{Z}_{k}:=Z_{k}\sqrt{\alpha_{k}}/\sqrt{\alpha}=Z_{k}/\sqrt{\text{Var}_{\mathbb{P}}(\phi^{k})}$
are i.i.d. standard normal and $w_{k}:=\sqrt{\alpha_{k}}$. Please note that
$\sum_{k}w_{k}^{2}=1$. Thus, we can write this equation

$L_{K}=\|\tilde{Z}-w(w\cdot\tilde{Z})\|_{2}^{2}=\|(\text{Id}-\Pi)\tilde{Z}\|_{2}^{2},$
(73)

where $\Pi$ projects on the one-dimensional subspace spanned by $w$. Let
$b_{1},\ldots,b_{K-1}$ be an orthonormal basis of the span of $\Pi$. Then, by
rotational invariance of the $\ell_{2}$ norm,

$L_{K}=\|(\text{Id}-\Pi)\tilde{Z}\|_{2}^{2}=\sum_{k=1}^{K-1}(b_{k}\cdot\tilde{Z})^{2}.$
(74)

Furthermore, since the $b_{k}$ are orthogonal to each other
$b_{k}\cdot\tilde{Z}$ are independent standard Gaussians. Thus, $L_{K}$
follows a $\chi^{2}(K-1)$ distribution. Furthermore, since the $b_{k}$ are
orthogonal to $w$, $L_{K}$ is independent of

$\sum_{k=1}^{K}w_{k}\tilde{Z}_{k}.$ (75)

Furthermore, by definition

$\frac{1}{\sqrt{\alpha}}\sum_{k=1}^{K}\alpha_{k}Z_{k}=\sum_{k=1}^{K}w_{k}\tilde{Z}_{k}\sim\mathcal{N}(0,1)$
(76)

and thus $L_{K}$ is independent of
$\frac{1}{\sqrt{\alpha}}\sum_{k=1}^{K}\alpha_{k}Z_{k}$. Therefore, (71) holds.
Using (71) with (69) and (67), we get

$\frac{\hat{\theta}^{W}-\theta^{0}}{\hat{\sigma}_{bet}/\sqrt{K-1}}\stackrel{{\scriptstyle
d}}{{=}}\frac{\delta\sum_{k=1}^{K}\alpha_{k}Z_{k}}{\sqrt{\alpha}\delta\sqrt{L_{K}/(K-1)}}+o_{P}(1)=\frac{\frac{\sum_{k}\alpha_{k}Z_{k}}{\sqrt{\alpha}}}{\sqrt{L_{K}/(K-1)}}=T_{K-1}+o_{P}(1).$

where $T_{K-1}$ is a t-distributed random variable with $K-1$ degrees of
freedom. This completes the proof. ∎

###  B.6 Proof of Theorem 2

###### Proof.

In this proof, if not specified otherwise, all variances and covariances are
meant with respect to $\mathbb{Q}$, that is marginally over both the variation
in $D$ and $\xi$. We will directly work with $U=h(D)$. Define

$f(x)=\text{Var}(\mathbb{P}^{\xi}(U\in[0,x)))$

Let $A$ and $B$ be two disjoint subsets of $[0,1]$. Define $a=\mathbb{P}(U\in
A)$ and $b=\mathbb{P}(U\in B)$. Then,

$\displaystyle f(a+b)$ $\displaystyle=\text{Var}(\mathbb{P}^{\xi}(U\in A\cup
B))$ $\displaystyle=\text{Var}(\mathbb{P}^{\xi}(U\in
A))+\text{Var}(\mathbb{P}^{\xi}(U\in B))+2\text{Cov}(\mathbb{P}^{\xi}(U\in
A),\mathbb{P}^{\xi}(U\in B))$
$\displaystyle=f(a)+f(b)+2\text{Cov}(\mathbb{P}^{\xi}(U\in
A),\mathbb{P}^{\xi}(U\in B))$

Thus, for any two disjoint sets $A$ and $B$,

$\text{Cov}(\mathbb{P}^{\xi}(U\in A),\mathbb{P}^{\xi}(U\in
B))=\frac{f(a+b)-f(a)-f(b)}{2}$

Define

$g(a,b)=\frac{f(a+b)-f(a)-f(b)}{2}.$ (77)

Let us first show that $f$ is continuous. Let $a_{n}\rightarrow a$, $a_{n}\geq
a$. Then,

$f(a_{n})-f(a)=f(a)+2g(a_{n}-a,a)+f(a_{n}-a)-f(a)=2g(a_{n}-a,a)+f(a_{n}-a).$

By Cauchy-Schwartz,

$g(a_{n}-a,a)\leq\sqrt{f(a_{n}-a)f(a)}.$

By assumption, $f(a_{n}-a)\rightarrow 0$. Thus, $f(a_{n})\rightarrow f(a)$.
The case $a_{n}\rightarrow a$, $a_{n}\leq a$ can be treated analogously. Thus,
$f(\bullet)$ and $g(\bullet,\bullet)$ are continuous.

Partition the probability space into disjoint $D$-measurable events $A_{i}$,
$i=1,\ldots,n$ with $\mathbb{P}(U\in A_{i})=1/n$. Then,

$0=\text{Var}(\mathbb{P}^{\xi}(\cup A_{i})-\mathbb{P}(\cup
A_{i}))=nf(1/n)+n(n-1)g(1/n,1/n)$ (78)

Thus,

$g(1/n,1/n)=-1/(n-1)f(1/n)$ (79)

We will now show that $f(x)=x(1-x)\delta_{\text{dist}}^{2}$ for $x=1/2^{k}$,
where $\delta_{\text{dist}}^{2}=4f(1/2)$. This will show that up to the
constant $\delta_{\text{dist}}^{2}=4f(1/2)$, $f$ and $g$ are uniquely defined.
First, we will show this equality for $x=1/4$.

$f(1/2)=f(1/4)+f(1/4)+2g(1/4,1/4)=2f(1/4)-2/3f(1/4)$ (80)

Thus,

$f(1/2)=4/3f(1/4)$ (81)

Rearranging,

$f(1/4)=1/4(1-1/4)4f(1/2)=x(1-x)\delta_{\text{dist}}^{2}$ (82)

for $x=1/4$. Induction step: Assume that

$f(x)=x(1-x)\delta_{\text{dist}}^{2}$ (83)

for $x=1/2^{k}$. Now we want to show that

$f(x/2)=x/2(1-x/2)\delta_{\text{dist}}^{2}$ (84)

To this end, using (77) and (79),

$f(x)=f(x/2)+f(x/2)-2/(2/x-1)f(x/2)$ (85)

Thus,

$f(x)=(2-2x/(2-x))f(x/2)=(4-2x-2x)/(2-x)f(x/2)=(4-4x)/(2-x)f(x/2).$ (86)

By induction assumption,

$f(x/2)=(2-x)/(4-4x)x(1-x)\delta_{\text{dist}}^{2}=x/2(1-x/2)\delta_{\text{dist}}^{2}.$
(87)

Thus, by induction for all $x=1/2^{k}$

$f(x)=x(1-x)\delta_{\text{dist}}^{2}.$ (88)

Now we want to show that for any $k$ and $j\leq 2^{k}$ and $x=j/2^{k}$

$f(x)=x(1-x)\delta_{\text{dist}}^{2}.$ (89)

For any $k$ and $j$ with $1\leq j\leq 2^{k}$, using the definition of $f$ and
(79),

$\displaystyle
f(j/2^{k})=jf(1/2^{k})-j(j-1)/(2^{k}-1)f(1/2^{k})=(j2^{k}-j^{2})/(2^{k}-1)f(1/2^{k})$
(90)
$\displaystyle=(j2^{k}-j^{2})/(2^{k}-1)1/2^{k}(1-1/2^{k})\delta_{\text{dist}}^{2}=(j2^{k}-j^{2})1/2^{k}1/2^{k}\delta_{\text{dist}}^{2}=j/2^{k}(1-j/2^{k})\delta_{\text{dist}}^{2}$
(91)

Thus, for all $k$ and $j\leq 2^{k}$, and $x=j/2^{k}$,

$f(x)=x(1-x)\delta_{\text{dist}}^{2}.$ (92)

Using continuity of $f$, for all $x\in[0,1]$

$f(x)=x(1-x)\delta_{\text{dist}}^{2}.$ (93)

We will now derive an explicit formula for $g$. For any $k$ and $j,j^{\prime}$
with $j+j^{\prime}\leq 2^{k}$,

$g(j/2^{k},j^{\prime}/2^{k})=jj^{\prime}g(1/2^{k},1/2^{k})=-jj^{\prime}/(2^{k}-1)f(1/2^{k})=-jj^{\prime}/(2^{k}-1)1/2^{k}(1-1/2^{k})\delta_{\text{dist}}^{2}$
(94)
$=-jj^{\prime}1/2^{k}1/2^{k}\delta_{\text{dist}}^{2}=-j/2^{k}j^{\prime}/2^{k}\delta_{\text{dist}}^{2}.$
(95)

By continuity, for all $x\geq 0$, $y\geq 0$ with $x+y\leq 1$,

$g(x,y)=-xy\delta_{\text{dist}}^{2}.$ (96)

Now assume that for some $D$-measurable disjoint sets $A_{i}$ and some
constants $y_{i}$,

$\psi(D)=\sum 1_{A_{i}}y_{i}$ (97)

Then,

$\text{Var}(\mathbb{E}^{\xi}[\psi(D)]-\mathbb{E}[\psi(D)])=\sum_{i}y_{i}^{2}f(\mathbb{P}(A_{i}))+\sum_{i\neq
j}y_{i}y_{j}g(\mathbb{P}(A_{i}),\mathbb{P}(A_{j}))$ (98)

To simplify, let’s write $p_{i}=\mathbb{P}(A_{i})$. Using explicit formulas
for $f$ and $g$,

$\text{Var}(\mathbb{E}^{\xi}[\psi(D)]-\mathbb{E}[\psi(D)])=\sum_{i}\delta_{\text{dist}}^{2}p_{i}(1-p_{i})-\sum_{i\neq
j}\delta_{\text{dist}}^{2}y_{i}y_{j}p_{i}p_{j}$ (99)

On the other hand,

$\delta_{\text{dist}}^{2}\text{Var}_{\mathbb{P}}(\psi(D))=\delta_{\text{dist}}^{2}(\sum_{i}p_{i}(1-p_{i})y_{i}^{2}+\sum_{i\neq
j}\text{Cov}(1_{A_{i}},1_{A_{j}})y_{i}y_{j})$ (100)

As the sets are disjoint, $\text{Cov}(1_{A_{i}},1_{A_{j}})=-p_{i}p_{j}$. Thus,

$\delta_{\text{dist}}^{2}\text{Var}_{\mathbb{P}}(\psi(D))=\delta_{\text{dist}}^{2}(\sum_{i}p_{i}(1-p_{i})y_{i}^{2}+\sum_{i\neq
j}-p_{i}p_{j}y_{i}y_{j})$ (101)

Combining equation (99) with equation (101),

$\text{Var}(\mathbb{E}^{\xi}[\psi(D)]-\mathbb{E}[\psi(D)])=\delta_{\text{dist}}^{2}\text{Var}_{\mathbb{P}}(\psi(D)).$
(102)

By measure-theoretic induction, this result is extended to any $\psi(D)\in
L^{2}(\mathbb{P})$. ∎

###  B.7 Asymptotic behaviour of _M_ -estimators

####  B.7.1 Proof of Lemma 2

###### Proof.

The proof proceeds as in van der Vaart [2000] with $m_{\theta}=-L(\theta,D)$
and $\mathbb{P}_{n}$ denoting the empirical measure of
$D_{1}^{n},\ldots,D_{n}^{n}$ and $d(\cdot,\cdot)$ the $\ell_{2}$ norm in
$\mathbb{R}^{d}$. The only difference is in the step where van der Vaart
[2000] uses the law of large numbers to show that with probability going to 1,
we have

$\frac{1}{n}\sum_{i=1}^{n}m_{U_{\theta_{j}}}(D_{i}^{n})=\mathbb{E}[m_{U_{\theta_{j}}}(D)]+o_{p}(1)$
(103)

for all $j=1,\ldots,p$. Since the data is not i.i.d. we cannot use the
classical law of large numbers to justify this step. Instead, we apply Lemma 1
to $\psi(D_{i}^{n})=m_{U_{\theta_{j}}}(D_{i}^{n})$ to see that

$\frac{1}{n}\sum_{i=1}^{n}m_{U_{\theta_{j}}}(D_{i}^{n})=\mathbb{E}[m_{U_{\theta_{j}}}(D)]+o_{p}(1).$
(104)

∎

####  B.7.2 Proof of Lemma 3

###### Proof.

The proof follows van der Vaart [2000], Theorem 5.41. with
$\Psi_{n}(\theta)=\frac{1}{n}\sum_{i=1}^{n}\partial_{\theta}L(\theta,D_{i}^{n})$
and
$\dot{\Psi}_{n}=\frac{1}{n}\sum_{i=1}^{n}\partial_{\theta}^{2}L(\theta,D_{i}^{n})$.
There are only three steps where the arguments are slightly different. First,
we cannot justify the bound

$\ddot{\Psi}_{n}(\tilde{\theta}_{n})=O_{P}(1)$ (105)

via the law of large numbers since the data is not i.i.d. Instead, we can use
Lemma 1 for $\psi(D)=h(D)$ to get

$\|\ddot{\Psi}_{n}(\tilde{\theta}_{n})\|\leq\frac{1}{n}\sum_{i=1}^{n}h(D_{i}^{n})=O_{P}(1).$
(106)

Analogously, we cannot apply the law of large numbers to
$\dot{\Psi}_{n}(\theta^{0})=\frac{1}{n}\sum_{i=1}^{n}\partial_{\theta}^{2}L(\theta^{0},D_{i}^{n})$.
Instead, we can use Lemma 1 with
$\psi(D)=\partial_{\theta}^{2}L(\theta^{0},D)$ to get

$\frac{1}{n}\sum_{i=1}^{n}\partial_{\theta}^{2}L(\theta^{0},D_{i}^{n})=\mathbb{E}[\partial_{\theta}^{2}L(\theta^{0},D)]+o_{P}(1).$
(107)

Similarly, we show that
$\Psi_{n}(\theta^{0})=\frac{1}{n}\sum_{i=1}^{n}\partial_{\theta}L(\theta^{0},D_{i}^{n})=O_{P}(1/\sqrt{n})$.
The rest of the proof proceeds as in van der Vaart [2000], Theorem 5.41. ∎

##  Appendix C Robust calibrated inference

In some cases, the data analyst may trust one of the estimators
$\hat{\theta}^{k}$ more than others. For example, the data analyst may be
convinced that $\theta^{1}=\theta^{0}$ but may not be sure whether
$\theta^{k}=\theta^{0}$ for $k\geq 2$. In this case, the data analyst may
report the confidence interval for $\theta^{0}$ using $\hat{\theta}^{1}$
instead of $\hat{\theta}^{W}$ with $\delta$ estimated by looking at the
between-estimator variance of the remaining $K-1$ estimators. Now we present
how to build asymptotic valid confidence intervals in such cases.

######  Theorem 3.

(Asymptotic validity of calibrated confidence interval). Suppose Assumption 1
holds for $k=1,\ldots,K$ and the influence functions
$\phi^{1}(D),\dots,\phi^{K}(D)$ are uncorrelated. Suppose
$\theta^{1}=\theta^{0}$ but $\theta^{k}$ may not be $\theta^{0}$ for $k\geq
2$. Furthermore assume that we have consistent estimates of the variances of
influence functions such that
$\widehat{\text{Var}_{\mathbb{P}}(\phi^{k}(D))}=\text{Var}_{\mathbb{P}}(\phi^{k}(D))+o_{p}(1)$
for $k=1,\dots,K$. Let
$\hat{\theta}^{W}=\sum_{k=2}^{K}\hat{\alpha}_{k}\hat{\theta}^{k}$ be the
inverse-variance weighted estimator of $K-1$ estimators where the weights are

$\hat{\alpha}_{k}=\frac{\frac{1}{\widehat{\text{Var}_{\mathbb{P}}(\phi^{k}(D))}}}{\sum_{j=2}^{K}\frac{1}{\widehat{\text{Var}_{\mathbb{P}}(\phi^{j}(D))}}}.$
(108)

Let $\hat{\sigma}_{\text{bet}}$ be the weighted between-estimator variance of
$K-1$ estimators defined as

$\hat{\sigma}^{2}_{\text{bet}}=\sum_{k=2}^{K}\hat{\alpha}_{k}(\hat{\theta}^{k}-\hat{\theta}^{W})^{2}.$
(109)

Then for any $\alpha\in(0,1)$, it holds that as $n\xrightarrow[]{}\infty$,

$\lim\inf_{n\rightarrow\infty}\mathrm{P}\left(\theta^{0}\in\Big{[}\hat{\theta}^{1}\pm
t_{K-2,1-\alpha/2}\cdot\sqrt{\sum_{j=2}^{K}\frac{\widehat{\text{Var}_{\mathbb{P}}(\phi^{1}(D))}}{\widehat{\text{Var}_{\mathbb{P}}(\phi^{j}(D))}}}\frac{\hat{\sigma}_{\text{bet}}}{\sqrt{K-2}}\Big{]}\right)\geq
1-\alpha,$ (110)

where $t_{K-2,1-\alpha/2}$ is the $1-\alpha/2$ quantile of the $t$
distribution with $K-2$ degrees of freedom. To be clear, here we marginalize
over both the randomness due to sampling and the randomness due to the
distributional perturbation.

The resulting confidence intervals are expected to be conservative. Firstly,
we lose one degree of freedom of the $t$-distribution. Secondly, we get an
overcoverage if $\theta^{k}\neq\theta^{0}$ for $k\geq 2$.

###### Proof.

If $\theta^{k}\neq\theta^{0}$ for some $k\geq 2$, then by asymptotic linearity
$\hat{\sigma}^{2}_{bet}$ converges to some $\tau^{2}>0$. As in the proof of
Theorem 1, we get
$\hat{\theta}^{1}-\theta^{0}=\mathcal{N}(0,\delta^{2}\text{Var}_{\mathbb{P}}(\phi^{1})/n)+o_{P}(1/\sqrt{n})$.
Since the variance estimates are consistent,

$\mathrm{P}\left(\theta^{0}\in\Big{[}\hat{\theta}^{1}\pm
t_{K-2,1-\alpha/2}\cdot\sqrt{\sum_{j=2}^{K}\frac{\widehat{\text{Var}_{\mathbb{P}}(\phi^{1}(D))}}{\widehat{\text{Var}_{\mathbb{P}}(\phi^{j}(D))}}}\frac{\hat{\sigma}_{\text{bet}}}{\sqrt{K-2}}\Big{]}\right)\rightarrow
1.$ (111)

Now let us consider the case $\theta^{0}=\theta^{1}=\ldots=\theta^{K}$. From
the proof of Theorem 1, we know that

$\frac{\sqrt{n}(\hat{\theta}^{1}-\theta^{0})}{\sqrt{\text{Var}_{\mathbb{P}}(\phi^{1}(D))}}\stackrel{{\scriptstyle
d}}{{=}}\delta Z+o_{p}(1)$

where $Z\sim N(0,1)$. Moreover,

$\displaystyle n\hat{\sigma}^{2}_{\text{bet}}\stackrel{{\scriptstyle
d}}{{=}}\delta^{2}\frac{1}{\sum_{j=2}^{K}\frac{1}{\text{Var}_{\mathbb{P}}(\phi^{j}(D))}}\cdot
L_{K-1}+o_{p}(1),$

where $L_{K-1}$ follows the chi-square distribution with $K-2$ degrees of
freedom. Note that $Z$ and $L_{K-1}$ are independent. Then, we get

$\frac{\frac{\hat{\theta}^{1}-\theta^{0}}{\sqrt{\widehat{\text{Var}_{\mathbb{P}(\phi_{1})}}}}}{\hat{\sigma}_{bet}\sqrt{\sum_{j=2}^{K}\frac{1}{\widehat{\text{Var}_{\mathbb{P}(\phi_{k})}}}}/\sqrt{K-2}}\stackrel{{\scriptstyle
d}}{{=}}\frac{Z}{\sqrt{L_{K-1}/(K-2)}}+o_{p}(1).$

Thus, we get

$\displaystyle P\left(\frac{Z}{\sqrt{L_{K-1}/(K-2)}}\leq x\right)$
$\displaystyle=P\left(t_{K-2}\leq x\right)$

where $t_{K-2}$ is a $t$-distributed random variable with $K-2$ degrees of
freedom. This completes the proof. ∎