Using gradients to check sensitivity of MCMC-based analyses to removing data

Tin Nguyen 1 Ryan Giordano 2 Rachael Meager 3 Tamara Broderick 1

Abstract

et al. (2015). More broadly, many researchers analyze datawith Bayesian models and MCMC (Senf et al., 2020; Meager, 2022; Jones et al., 2021; Porter et al., 2022) and wantto know if their conclusions generalize beyond their data.

If the conclusion of a data analysis is sensitive
to dropping very few data points, that conclusion
might hinge on the particular data at hand rather
than representing a more broadly applicable truth.
To check for this sensitivity, one idea is to consider every small data subset, drop it, and re-run
our analysis. But the number of re-runs needed
is combinatorially large. Recent work proposesa differentiable relaxation to find the worst-case
subset, but that work was developed for conclusions based on estimating equations — and does
not directly handle Bayesian posterior approximations using MCMC. We make two principal contributions. We adapt the existing data-dropping
relaxation to estimators computed via MCMC; inparticular, we re-use existing MCMC draws toestimate the necessary derivatives via a covariance relationship. Observing that Monte Carloerrors induce variability in the estimates, we use
a variant of the bootstrap to quantify this uncertainty. Empirically, our method is accurate insimple models, such as linear regression. In models with complex structure, such as hierarchies,
the performance of our method is mixed.

If one is interested in generalization, an intuitive idea ischecking if a conclusion changes after dropping very few
data points. For instance, we show in Section 4.1 that after removing less than 0.1% of the RCT data (Angelucci
et al., 2015), microcredit is estimated to increase, rather
than decrease, profit. So, the conclusion hinges on a small
number of businesses. We do not know if these data points
appear elsewhere. Hence, we worry about generalization.
For more motivation on why small-data sensitivity is related
to generalization, see Appendix A.1.
One might want to know if a similar sensitivity exists forother analyses. A natural idea is to enumerate over every data subset and re-analyze. Unfortunately, the numberof subsets to search over is prohibitively large. See Appendix A.2 for an estimate on brute-force runtime. So weturn to approximations. While there has been previous workon approximating dropping worst-case data subsets (Broderick et al., 2023; Kuschnig et al., 2021; Shiffman et al., 2023;
Moitra and Rohatgi, 2022; Freund and Hopkins, 2023), none
directly apply to MCMC: the focus of past works have been
either estimating equations (Kosorok, 2008)[Chapter 13] orordinary least squares (OLS).

Consider this motivating example. Angelucci et al. (2015)
conducted a randomized controlled trial (RCT) in Mexico
to study whether microcredit loans improve business profits. To analyze this data, one might use a simple Bayesianmodel and Markov chain Monte Carlo (MCMC). Based on
the posterior mean MCMC estimate, one might conclude
that microcredit actually reduces profit in this study. Next, if
a policymaker wants to advocate against microcredit deployment outside of Mexico, they need to know if microcredit
remains detrimental beyond the data gathered in Angelucci

Our work extends Broderick et al. (2023) to conclusions
based on MCMC. In Section 3.1, similar to Broderick et al.
(2023), we use a first-order Taylor series. We observethat the first derivatives can be interpreted as posterior covariances, and we use MCMC to estimate the covariances(Section 3.2). Recognizing that Monte Carlo errors induce
variability in our approximation, in Section 3.3 we use avariant of the bootstrap Efron (1979) to quantify this uncertainty. We provide a longer discussion of related work in
Appendix A.3. Experimentally, in Section 4, while our approximation performs well in simple models such as linear
regression, it is less reliable in complex models.

1
MIT 2 University of California, Berkeley 3 University of New
South Wales. Correspondence to: Tin Nguyen <tdn@mit.edu>.

2. Background

1. Introduction

We introduce notation in two parts. First, we cover thenotation and concepts involved in Bayesian data analysis.

Published at the 2 nd Differentiable Almost Everything Workshop
at the 41 st International Conference on Machine Learning, Vienna, Austria. July 2024. Copyright 2024 by the author(s).

1Small-data sensitivity of MCMC analyses

Second, we extend the notation to dropping data.

and standard likelihoods.
To indicate that an expectation is taken with respect to therandomness β ∼ p(β | w, {d(n) }Nn=1 ), we use Ew . The
value of a posterior functional depends on w. For instance,
the posterior mean under the weighted posterior is Ew g(β).

Bayesian data analysis. Suppose we have a dataset
V{d(n) }Nof interest.
n=1 . Consider a parameter β ∈ R
To estimate β, we take a Bayesian approach. First, we
model the link between β and the data through a likelihood L(d(n) | β). Secondly, we specify a prior distribution
over β, and use p(β) to denote the prior density. Then, thedensity of the posterior distributionof β given the data isQN(n)| β)).
p(β | {d(n) }Nn=1 ) ∝ p(β)
n=1 exp(L(dThe non-robustness concern can be formalized as follows.
For α ∈ (0, 1), let Wα denote the set of all weight vectorsthat correspondthan 100α% of theo datan to dropping no morePN
i.e. Wα := w ∈ {0, 1}N : N1 n=1 (1 − wn ) ≤ α . We
say the analysis is non-robust if there exists a weight w ∈Wα that changes the conclusion.

In practice, an analyst uses a posterior functional to makeconclusions. One example is the posterior mean Eg(β),
where g is a mapping from RV to R. In linear regression,commonly a practitioner will make a decision based on thesign of the posterior mean of a particular regression coefficient. Other decisions are made with uncertainty intervals:
see Appendix B.1 for details.

We focus on decisions satisying the following simplifying
assumption: there exists a posterior functional (ϕ(w)) suchthat ϕ(1N ) < 0 and the conclusion changes if and only if
ϕ(w) > 0. Such a functional will be called a “quantity of interest” (QoI). The decision based on sign fits this framework.
To change this conclusion, if the full-data posterior mean
were positive, we take ϕ(w) = −Ew g(β). Appendix B.1
shows how other decision changes fit this framework.

Computationally, posterior functionals need not have closed
forms. To approximate posterior functionals, practitioners
frequently use MCMC methods. Let (β (1) , . . . , β (S) ) denote the MCMC draws that target the posterior distribution.
We estimate expectations using (β (1) , . . . , β (S) ), and makea decision based on such estimates.

Checking non-robustness is equivalent to a) finding the maximum value of ϕ(w) subject to w ∈ Wα and b) checking its
sign. The outcome of this comparison remains the same if
we maximize the objective function ϕ(w) − c and compare
the optimal value with −c, for any constant c. Out of later
convenience, we set c = ϕ(1N ). As in Broderick et al.
(2023, Section 2), we define the Maximum Influence Perturbation (MIP) to be the largest change induced by dropping
no more than 100α% of the data. In our notation, it is

Dropping data. A Bayesian analyst might be worried
if the substantive decision arising from their data analysis
changed after removing some small fraction α of the data.
For instance, if their decision were based on the sign of theposterior mean, they would be worried if that sign changed.
Other decision changes are given in Appendix B.1.

max (ϕ(w) − ϕ(1N )) .

To describe non-robustness precisely and to develop our
approximation, we need to indicate the dependence of posterior functionals on the presence of data points. We introduce
data weights w = (w1 , w2 , . . . , wN ). This vector defines
the so-called weighted posterior distribution.
Definition 2.1. Let Z(w) be the normalizing constant for
QNp(β) n=1 exp(wn L(d(n) | β)). If Z(w) < ∞, the density of the weighted posterior distribution associated withw is denoted by p(β | w, {d(n) }Nn=1
 ), and is equal toPN1(n)| β) .
n=1 wn L(dZ(w) p(β) exp

w∈Wα

(1)

In general, the brute-force approach to find the MIP takes aprohibitively long time; recall Appendix A.2.

3. Methods
We turn to relaxations of the MIP. We focus on ϕ(w) =Ew g(β); Appendices B.1 and B.3 discusses other QoIs.
3.1. Taylor series
Our first approximation relies on the first-order Taylor seriesof the quantity of interest ϕ(w). We require that the quantityof interest ϕ(w) is differentiable with respect to w. The
combination of Assumption 2.1 and Assumption B.2 from
Appendix B.2 forms a set of mild regularity conditions that
ensure this differentiability.
Theorem 3.1. Assume Assumption 2.1 and Assumption B.2.
For any δ ∈ (0, 1), ϕ(w) is continuously differentiable on{w ∈ [0, 1]N : maxn wn ≥ δ}. Then-th partial derivative
is equal to Covw g(β), L(d(n) | β) .

If wn = 0, the n-th observation is ignored; if wn = 1, the nth observation is fully included. We recover the regular posterior by setting all weightsto 1: w = 1N = (1, 1,. . . , 1).
PN(n)| β) is not
It is possible that p(β) exp
n=1 wn L(dintegrable. For example, the prior p(β) is improper and allweights have been set to zero: w = 0N = (0, 0, . . . , 0). We
assume that any likelihood contribution is enough to define
a proper posterior.
Assumption 2.1. ∀w ∈ [0, 1]N \ {0N }, Z(w) < ∞.

This assumption is immediate in the case of proper prior

See Proof D.1 for the proof, and connections to previous2

Small-data sensitivity of MCMC analyses

an interquantile range of this empirical distribution as ourconfidence interval for ∆(α). For details, see Appendix B.4.

works on sensitivity analysis of posterior expectations.
We define the n-th influence as the partial derivative of
ϕ(w) at w = 1N : ψn := ∂ϕ(w)∂wn w=1N . Then, the firstorder Taylor series approximation of ϕ(w) − ϕ(1N ) isPN
n=1 ψn (wn − 1). We approximately solve Equation (1)by replacing its objective function:
max

w∈Wα

N
X

n=1

(wn − 1)ψn

General MCMC. We now handle the case in which(β (1) , . . . , β (S) ) exhibit dependence. One idea is to use
the previous section’s construction without modification.
Theoretically, it is known that the bootstrap struggles on
non-i.i.d. samples, for even simple estimators. Intuitively,
the bootstrap fails because the bootstrap draws do not haveany dependence, while the original draws do. Appendix B.5
explains this in detail.

(2)

Solving Equation (2) involves a sort.
P For any w ∈ Wα ,the objective function is equal to n:wn =0 (−ψn ). Let
r1 , r2 , . . . , rN sort the ψn in increasing order: ψr1 ≤ ψr2 ≤. . . ≤ ψrN . Then, the optimal value is equal to the negaP⌊N α⌋
tive of m=1 ψrm I{ψrm < 0}, where I{·} is the indicatorfunction. We denote the optimal value by ∆(α), and thedata points to be removed by U (α).

To improve upon the bootstrap, one option is to resample
in a way that respects the original sample’s dependence
structure. In particular, we use the non-overlapping blockbootstrap (Lahiri, 2003; Carlstein, 1986): instead of resampling individual Markov chain states, we first divide theMarkov chain into blocks, and then resample blocks. For
details, see Appendix B.5. The outcome of this resamplingis an interval, denoted by [∆lb (α), ∆ub (α)], which is ourconfidence interval for ∆(α).

3.2. Estimating the influence
To solve Equation (2), we need to compute the influenceψn . Each ψn is a covariance under the full-data posterior.
Therefore, the MCMC draws, which are already used toestimate ϕ(1N ), can be used to estimate ψn : ψn ≈ ψ̂nwhere ψ̂n is the sample covariance between L(d(n) | β (s) )and g(β (s) ) (the sample being (β (1) , . . . , β (S) )).

4. Experiments
Now, we check the quality of our approximations empirically on real data analyses. For a particular MCMC run, we
estimate sensitivity for a range of α values: see Appendix Efor details. For each α, our method proposes an influentialb ) and a change in the quantity of interest, repdata subset (Uresented by a confidence interval ([∆lb (α), ∆ub (α)]). We
plot how the change from re-running minus the proposed
data compares to the confidence interval.

Since ψ̂n is only an approximation of ψn , we are not able
to solve Equation (2) exactly, but only solve an approximation of it. Namely, we replace all instances
PN of ψn with ψ̂nin Equation (2), and solve: maxw∈Wα n=1 (wn − 1)ψ̂n .
The algorithm to solve this is analogous of that for Equation (2). For instance, let the ranks v1 , v2 , . . . , vN be suchthat ψ̂v1 ≤ ψ̂v2 ≤ . . . ≤ ψ̂vN . Then the optimal value isP⌊N α⌋b− m=1 ψ̂vm I{ψ̂vm < 0}. We denote this value by ∆,b.
and the data points to be removed by U

4.1. Linear model
We consider a slight variation of Meager (2019)’s analysisof the Angelucci et al. (2015) RCT. For context and experimental details, see Appendix F.1. In this simple linear
model, our prediction typically contains the rerun.

3.3. Confidence intervals

The most interesting parameter, θ, is the treatment effect: thedifference in profit between treatment and control. Figure 1
plots the histogram of the treatment effect draws as well as
key sample summaries. The sample mean is equal to −4.55.
The sample standard deviation is 5.79. The uncertainty
interval is (−16.10, 6.99). An analyst might conclude thatwhile the posterior mean of the effect of microcredit isnegative, the uncertainty interval covers zero, so they cannot
confidently conclude that microcredit either helps or hurts.

b is a noisy point estimate of ∆(α). One concern regarding∆
b is noise due to sampling uncertainty ofthe quality of ∆(β (1) , . . . , β (S) ). So, we design confidence intervals.
Exact sampling. In certain cases, such as conjugate models (Diaconis and Ylvisaker, 1979), we can generate exact
Monte Carlo draws from the posterior distribution. Then,b is an estimator constructed from an i.i.d.
conceptually, ∆sample. We appeal to the bootstrap (Efron, 1979), a generalpurpose technique to quantify the sampling uncertainty.

In Figure 2, we plot our confidence intervals and the result
after removing the proposed data. For changing sign, ourmethod predicts there exists a data subset at most 0.1% suchthat if we remove it, we change the posterior mean’s sign.
Refitting after removing the proposed data confirms thisprediction. For changing significance, our method predicts

Our confidence interval construction proceeds in three steps.
bFirst, we define the so-called bootstrap distribution of ∆.
Second, we approximate this distribution with an empirical
distribution based on Monte Carlo draws. Finally, we use3

Small-data sensitivity of MCMC analyses

Figure 1. (Linear model) Histogram of treatment effect MCMC
draws. Blue line is sample mean. Dashed red line is zero threshold.
Dotted blue lines are uncertainty interval endpoints.

Figure 3. (Hierarchical model for tree mortality) Histogram ofslope MCMC draws. See the caption of Figure 1 for the meaningof the distinguished vertical lines.

Figure 2. (Linear model) Confidence interval and refit. At maximum, we remove 1% of the data. Each panel corresponds to atarget conclusion change: ‘sign’ is the change in sign, ‘sig’ ischange in significance, and ‘both’ is the change to a significant
effect of the opposite sign. Error bars are confidence interval for
refit after removing the most extreme data subset. Each ‘x’ is therefit after removing the proposed data and re-running MCMC. The
dotted blue line is the fit on the full data.

Figure 4. (Hierarchical model for tree mortality) Confidence interval and refit. See the caption of Figure 2 for the meaning of thepanels and the distinguished lines.

general, our interval is not conservative. The overestimation
is particularly severe for the ‘both’ QoI and the ‘sig’ QoI.
For changing sign, our method predicts there exists a datasubset of relative size at most 0.17% such that if we removeit, we change the posterior mean’s sign; refitting does notconfirm this prediction, however. The smallest α whose
refit’s posterior mean actually changes sign is 0.22%. Forchanging significance, our method predicts there exists adata subset of relative size at most 0.10% such that if weremove it, we change the sign of the right endpoint; refitting
confirms this prediction. For generating a significant resultof the opposite sign, our method predicts there exists a datasubset of relative size at most 0.17% such that if we removeit, we change the sign of the left endpoint; refitting does
not confirm this prediction, however. The smallest α whose
refit’s left endpoint actually changes sign is 1.0%.

there exists a data subset of relative size at most 0.36% suchthat if we remove it, we change the sign of the uncertainty
interval’s right endpoint: refitting confirms this prediction.
Our method is not able to predict whether the result can bechanged to significant effect of the opposite sign for these
α values and this number of draws: we recommend more
MCMC draws, which should decrease interval width.
4.2. Hierarchical model
We consider a slight variation of the analysis of European
tree mortality from Senf et al. (2020). For details, see Appendix F.2. Our approximation struggles: we predict more
extreme changes than realized by re-running.

5. Next Steps

The most interesting parameter, θ, is the association between
water availability and the amount of tree canopy death. Figure 3 plots the histogram of the association effect draws
and sample summaries. The sample mean is equal to −1.88.
The sample standard deviation is 0.48. The uncertainty
interval is (−2.81, −0.94). One might decide that waterbalance has a negative relationship with canopy mortality,
since the posterior mean is negative, and this relationship is
significant, since the uncertainty interval omits zero.

We have provided a fast approximation to what happens toconclusions made with MCMC in Bayesian models whena small percentage of data is removed. In real data experiments, our approximation is accurate in simple models,such as linear regression. In complicated models, such as
hierarchical ones, our methods are less accurate. Naturally,
we would like to understand why the performance in mixed.
In Appendix E, we describe the checks that could be run
to identify error sources. Appendix F.1 and Appendix F.2
show the outcomes of running such checks.

Figure 4 plots our confidence intervals and the reruns. In
4Small-data sensitivity of MCMC analyses

6. Acknowledgments

B. Efron, “Bootstrap Methods: Another Look at the Jackknife,” Annals of Statistics, vol. 7, pp. 1–26, 1979.

The authors are grateful to Hannah Diehl for useful discussions and comments. This work was supported in part by an
ONR Early Career Grant (N000142012532), an NSF CAREER Award (1750286), and SystemsThatLearn@CSAIL
Ignite Grant.

P. Diaconis and D. Ylvisaker, “Conjugate Priors for Exponential Families,” The Annals of Statistics, vol. 7, pp.
269–281, 1979.
S. N. Lahiri, Resampling Methods for Dependent Data.
Springer, 2003.

References

E. Carlstein, “The Use of Subseries Values for Estimating
the Variance of a General Statistic from a Stationary Sequence,” Annals of Statistics, vol. 14, pp. 1171–1179,
1986.

M. Angelucci, D. Karlan, J. Zinman, K. Brennan, E. Degnan,
A. Fishbane, A. Hillis, H. Koizumi, E. Safran, R. Strohm,
B. Torres, A. Troychansky, I. Velez, G. Startz, S. Swamy,
M. White, A. York, and C. Banco, “Microcredit Impacts:
Evidence from a Randomized Microcredit Program Placement Experiment by Compartamos Banco,” American
Economic Journal: Applied Economics, vol. 7, pp. 151–
82, 2015.

R. Meager, “Understanding the Average Impact of Microcredit Expansions: A Bayesian Hierarchical Analysis of
Seven Randomized Experiments,” American Economic
Journal: Applied Economics, vol. 11, pp. 57–91, 2019.
G. Arya, M. Schauer, F. Schäfer, and C. Rackauckas, “Automatic differentiation of programs with discrete randomness,” Advances in Neural Information Processing Systems, vol. 35, pp. 10 435–10 447, 2022.

C. Senf, A. Buras, C. S. Zang, A. Rammig, and R. Seidl,“Excess forest mortality is consistently linked to drought
across europe,” Nature Communications, vol. 11, 12
2020.

G. Arya, R. Seyer, F. Schäfer, K. Chandra, A. K. Lew,
M. Huot, V. K. Mansinghka, J. Ragan-Kelley, C. Rackauckas, and M. Schauer, “Differentiating MetropolisHastings to Optimize Intractable Densities,” 2023.

R. Meager, “Aggregating distributional treatment effects: A
bayesian hierarchical analysis of the microcredit literature,” American Economic Review, vol. 112, no. 6, pp.
1818–47, June 2022.

R. Seyer, “Differentiable Monte Carlo samplers with piecewise deterministic Markov processes,” 2023.

T. C. Jones, G. Biele, B. Mühlemann, T. Veith, J. Schneider, J. Beheim-Schwarzbach, T. Bleicker, J. Tesch, M. L.
Schmidt, L. E. Sander, F. Kurth, P. Menzel, R. Schwarzer,
M. Zuchowski, J. Hofmann, A. Krumbholz, A. Stein,
A. Edelmann, V. M. Corman, and C. Drosten, “Estimating
infectiousness throughout SARS-CoV-2 infection course,”
Science, vol. 373, 7 2021.

J. P. Kleijnen and R. Y. Rubinstein, “Optimization and sensitivity analysis of computer simulation models by thescore function method,” European Journal of Operational
Research, vol. 88, no. 3, pp. 413–427, 1996.
M. C. Fu and J.-Q. Hu, “Conditional Monte Carlo gradientestimation,” in Conditional Monte Carlo: Gradient estimation and optimization applications. Springer Science& Business Media, 2012, vol. 392.

T. Porter, D. C. Molina, A. Cimpian, S. Roberts, A. Fredericks, L. S. Blackwell, and K. Trzesniewski, “GrowthMindset Intervention Delivered by Teachers Boosts
Achievement in Early Adolescence,” Psychological Science, vol. 33, pp. 1086–1096, 7 2022.

B. Heidergott and F. Vázquez-Abad, “Measure-valued differentiation for Markov chains,” Journal of Optimization
Theory and Applications, vol. 136, no. 2, pp. 187–209,
2008.

T. Broderick, R. Giordano, and R. Meager, “An Automatic
Finite-Sample Robustness Metric: Can Dropping a Little
Data Change Conclusions?” 2023.
N. Kuschnig, G. Zens, and J. C. Cuaresma, “Hidden in Plain
Sight: Influential Sets in Linear Models,” 2021.

R. Giordano, T. Broderick, and M. I. Jordan, “Covariances,
Robustness, and Variational Bayes,” Journal of Machine
Learning Research, vol. 19, pp. 1–49, 2018.

M. Shiffman, R. Giordano, and T. Broderick, “Could dropping a few cells change the takeaways from differential
expression?” 2023.

S. Mohamed, M. Rosca, M. Figurnov, and A. Mnih, “Monte
Carlo Gradient Estimation in Machine Learning,” Journalof Machine Learning Research, vol. 21, pp. 1–62, 2020.

A. Moitra and D. Rohatgi, “Provably Auditing Ordinary
Least Squares in Low Dimensions,” 2022.

R. Giordano and T. Broderick, “The Bayesian Infinitesimal
Jackknife for Variance,” 2023.

D. Freund and S. B. Hopkins, “Towards Practical Robustness Auditing for Linear Regression,” 2023.

R. Giordano, R. Liu, M. I. Jordan, and T. Broderick, “Evaluating Sensitivity to the Stick-Breaking Prior in BayesianNonparametrics (with Discussion),” Bayesian Analysis,
vol. 18, pp. 287–366, 2023.

M. R. Kosorok, Introduction to Empirical Processes andSemiparametric Inference. Springer, 2008.
5

Small-data sensitivity of MCMC analyses

P. Diaconis and D. Freedman, “On the Consistency of Bayes
Estimates,” The Annals of Statistics, pp. 1–26, 1986.

Journal of the American Statistical Association, vol. 108,
pp. 1339–1349, 2013.

F. Ruggeri and L. Wasserman, “Infinitesimal sensitivity of
posterior distributions,” The Canadian Journal of Statistic, vol. 21, pp. 195–203, 1993.

W. H. Fleming, Functions of Several Variables, 2nd ed.
Springer, 1977.
A. W. van der Vaart, Asymptotic Statistics. University ofCambridge, 1998.

P. Gustafson, “Local Sensitivity of Posterior Expectations,”
The Annals of Statistics, vol. 24, p. 195, 1996.

C. J. Clopper and E. S. Pearson, “The Use of Confidence orFiducial Limits Illustrated in the Case of the Binomial,”
Biometrika, vol. 26, no. 4, pp. 404–413, 1934.

W. Johnson and S. Geisser, “A Predictive View of the Detection and Characterization of Influential Observations inRegression Analysis,” Journal of the American Statistical
Association, vol. 78, pp. 137–144, 1983.

B. Carpenter, A. Gelman, M. D. Hoffman, D. Lee,
B. Goodrich, M. Betancourt, M. Brubaker, J. Guo, P. Li,
and A. Riddell, “Stan: A Probabilistic Programming Language,” Journal of Statistical Software, vol. 76, pp. 1–32,
2017.

R. E. Mcculloch, “Local Model Influence,” Journal of theAmerican Statistical Association, vol. 84, pp. 473–478,
1989.
M. Lavine, “Local predictive influence in bayesian linear
models with conjugate priors,” Communications in Statistics - Simulation and Computation, vol. 21, pp. 269–283,
1 1992.
B. P. Carlin and N. G. Polson, “An Expected Utility Approach to Influence Diagnostics,” Journal of the American
Statistical Association, vol. 86, pp. 1013–1021, 1991.

A. Reuther, J. Kepner, C. Byun, S. Samsi, W. Arcand, D. Bestor, B. Bergeron, V. Gadepally, M. Houle, M. Hubbell,
M. Jones, A. Klein, L. Milechin, J. Mullen, A. Prout,
A. Rosa, C. Yee, and P. Michaleas, “Interactive supercomputing on 40,000 cores for machine learning anddata analysis,” in 2018 IEEE High Performance extreme
Computing Conference (HPEC). IEEE, 2018, pp. 1–6.

E. C. Marshall and D. J. Spiegelhalter, “Identifying outliers in Bayesian hierarchical models: a simulation-based
approach,” Bayesian Analysis, vol. 2, pp. 409–444, 2007.

S. N. Wood, “Thin plate regression splines,” Journal of theRoyal Statistical Society: Series B (Statistical Methodology), vol. 65, no. 1, pp. 95–114, 2003.

R. B. Millar and W. S. Stewart, “Assessment of Locally
Influential Observations in Bayesian Models,” Bayesian
Analysis, vol. 2, pp. 365–384, 2007.
A. van der Linde, “Local Influence on Posterior Distributions under Multiplicative Modes of Perturbation,”
Bayesian Analysis, vol. 2, pp. 319–332, 2007.
Z. M. Thomas, S. N. MacEachern, and M. Peruggia, “Reconciling Curvature and Importance Sampling Based Procedures for Summarizing Case Influence in BayesianModels,” Journal of the American Statistical Association,
vol. 113, pp. 1669–1683, 10 2018.
M. T. Pratola, E. I. George, and R. E. McCulloch, “Influential Observations in Bayesian Regression Tree Models,”
Journal of Computational and Graphical Statistics, 2023.
J. Huggins, M. Kasprzak, T. Campbell, and T. Broderick,“Validated variational inference via practical posterior
error bounds,” in International Conference on Artificial
Intelligence and Statistics. PMLR, 2020, pp. 1792–1802.
S. Basu, S. R. Jammalamadaka, and W. Liu, “Local posterior robustness with parametric priors: Maximum andaverage sensitivity,” in Maximum Entropy and Bayesian
Methods: Santa Barbara, California, U.S.A., 1993, G. R.
Heidbreder, Ed. Dordrecht: Springer Netherlands, 1996,
pp. 97–106.
N. G. Polson, J. G. Scott, and J. Windle, “Bayesian inference
for logistic models using Pólya-Gamma latent variables,”6Small-data sensitivity of MCMC analyses

A. Other Intro Details
A.1. More motivation of dropping data as generalization check
Standard tools to assess generalization often make assumptions that are not realistic in practice. For instance, an analyst
might use frequentist tools (confidence interval, p-values) to predict whether their inferences hold in the broader population.
The validity of these methods technically depends on the assumption that the gathered data is an independent and identically
distributed (i.i.d.) sample from the broader population. In practice, we have reason to suspect that this assumption is not met;for instance, it might not be reasonable to assume that data collected in Mexico and data collected in a separate country are
i.i.d. from the same distribution.
As pointed out by Shiffman et al. (2023), an analyst might hope that deviations from the i.i.d. assumption are small enoughthat (a) their conclusions remain the same in the broader population and (b) standard tools accurately assess generalization.
On the other hand, the analyst might worry that this hope is misplaced if small, realistic deviations from i.i.d.-ness could
affect the substantive conclusions of an analysis. An often-realistic kind of deviation is the missingness of a small fractionof data; for instance, some percentage of the population might not respond to a survey. So, if it were possible to remove asmall fraction of data and change conclusions, the analyst might worry about generalization.
A.2. Brute-force runtime
We need to enumerate every data subset that drops no more than 100α% of the original data.
 And, for each subset, we
would need to re-run MCMC to re-estimate the quantity of interest. There are more than ⌊NNα⌋ elements in Wα . The RCT
data from Angelucci et al. (2015) has N = 16,560 observations; even for α = 0.001, there are more than 1054 subsets toconsider. Each Markov chain already takes a noticeable amount of time to construct; in this analysis, to generate S = 4,000
draws, we need to run the chain for 1 minute. The total time to find the worst-case data subset is on the order of 1048 years.
A.3. Related work
Our work arguably fits into the intersection of three lines of work.
The first is papers on detecting sensitivity to small-data removal. Broderick et al. (2023) were the first to formulate sensitivity
to dropping a small fraction of data as a check on generalization. Along with the formulation, one contribution of thiswork is a fast approximation to detect sensitivity when the analysis in question is based on estimating equations (Kosorok,2008)[Chapter 13]. Regardless of how estimators are constructed, in general, the brute-force approach to finding an
influential small fraction of data is computationally intractable. One would need to enumerate all possible data subsets
of a given cardinality and re-analyze on each subset: even when the fraction of data removed is small and each analysis
takes little time, there are too many such subsets to consider; see the discussion at the end of Section 3. For estimating
equations, Broderick et al. (2023) approximate the effect of dropping data with a first-order Taylor series approximation;
this approximation can be optimized very efficiently, while the brute-force approach is not at all practical. Neither Broderick
et al. (2023) nor subsequent existing work on small-data removals (Kuschnig et al., 2021; Shiffman et al., 2023; Moitra andRohatgi, 2022; Freund and Hopkins, 2023) can be immediately applied to determine sensitivity in MCMC. In particular,
since MCMC cannot be cast as the root of an estimating equation or the solution to an optimization problem, neither
Broderick et al. (2023) nor Shiffman et al. (2023) apply to our situation. As Kuschnig et al. (2021); Moitra and Rohatgi(2022); Freund and Hopkins (2023) focus on ordinary least squares (OLS), their work does not address our problem, either.
The second line of work estimates the changes that happen to a posterior expectation because of small perturbations to thetotal log likelihood. There are two conceptually distinct approaches to this sensitivity analysis.
• One approach (e.g. Arya et al., 2022, 2023; Seyer, 2023) applies to when the posterior is approximated with a
Metropolis-Hastings algorithm. In particular, this approach computes the gradient of the Metropolis-Hastings sampler
to small perturbations in the total log likelihood. More broadly, there is a literature on estimating gradients for random
processes with discrete components (Kleijnen and Rubinstein, 1996; Fu and Hu, 2012; Heidergott and Vázquez-Abad,2008).
• The other approach does not compute the gradient of the MCMC algorithm or steps within it. Instead, it directlycomputes (and then estimates) the gradient of the posterior expectation. Recent works in this literature include Giordano
et al. (2018); Mohamed et al. (2020); Giordano and Broderick (2023); Giordano et al. (2023), while foundational works
7

Small-data sensitivity of MCMC analyses

include Diaconis and Freedman (1986); Ruggeri and Wasserman (1993); Gustafson (1996).
In our work, we take the second approach. A priori, it is not clear which approach is superior. Two reasons to prefer thesecond approach over the first approach are the following. While the discrete operations in Metropolis-Hastings, e.g. theaccept/reject steps, pose a key challenge in the first approach, they do not cause any issues in the second approach; thesecond approach is “oblivious” to details regarding how the posterior is approximated. In addition, suppose that an analyst
wishes to compute gradients of multiple quantities of interest. If they follow the first approach, for each quantity of interest,
they would need to re-run the sampling algorithm to estimate the gradient. Taking the second approach, the analyst only
needs to run the sampling algorithm once, and use the resulting draws to simultaneously estimate the gradient of multiple
quantities of interest. On the other hand, the first approach might be better than the second approach in the following way.
Our experiments later show that gradient estimates coming from the second approach can be noisy. The first approach,with the promise of variance reduction through a good choice of Markov chain coupling, might produce more accurate
gradient estimates. It is an interesting direction for future work to apply the first approach to our problem and compare theperformance of the two approaches.
While papers taking the second approach have already mentioned how to estimate the effect of dropping an individual
observation, these estimates have not been used to assess whether conclusions based on MCMC are sensitive to the removalof a small data fraction. Some works (e.g. Gustafson, 1996; Giordano et al., 2018, 2023) generate perturbations by varying
prior or likelihood choice. Giordano and Broderick (2023) estimate the frequentist variability of Bayesian procedures, a taskthat can be seen as equivalent to the goal of bootstrap resampling. No existing work aims to find a small fraction of data that,
if dropped, would change conclusions.
The third set of works, in the Bayesian case influence literature, quantifies the importance of individual observations to a
Bayesian analysis. As we will explain, existing works do not tackle our problem. Early works in this area include Johnson
and Geisser (1983); Mcculloch (1989); Lavine (1992); Carlin and Polson (1991), while recent works include Marshall andSpiegelhalter (2007); Millar and Stewart (2007); van der Linde (2007); Thomas et al. (2018); Pratola et al. (2023). Suchpapers focus on the identification of outliers, rather than predictions about whether the conclusion changes after removing
a small amount of data. Generally, this literature defines an observation to be an outlier if the Kullback–Leibler (KL)divergence between the posterior after removing the observation and the original posterior is large. For conclusions basedon posterior functionals, such as the mean, we are not aware of how to systematically connect the KL divergence to thesensitivity of the decision-making process; in fact, recent work (Huggins et al., 2020) has shown that comparing probability
distributions based on the KL divergence can be misleading if an analyst really cared about the comparison between the
distributions’ means or variances.

B. Other Methods Details
B.1. Other decisions and quantities of interest
In Section 2 and Section 3, we focus on decisions based on the sign of the posterior mean. Here, we describe other decisions,which are based on uncertainty intervals, and the corresponding quantities of interest.
An econometrician might declare that an intervention is helping some population if the vast majority of the posterior mass
for a particular coefficient lies above zero. That is, the practitioner checks if the lower bound of an uncertainty interval lies
above zero. This decision might be considered to reflect a Bayesian notion of significance.
A Bayesian analyst might be worried if the substantive decision arising from their data analysis changed after removingsome small fraction α of the data. For instance,• If their decision were based on zero falling outside a credible interval, they would be worried if we can make thecredible interval contain zero.
• If their decision was based on both the sign and the significance, they would be worried if we can both change theposterior mean’s sign and put a majority of the posterior mass on the opposite side of zero.
p
To change the conclusion about significance, if uncertainty interval’s left endpoint1 (E1N g(β) − z0.975 Var1N g(β)) were1
Our uncertainty interval multiplies the posterior standard deviation by z0.975 , which is the 97.5% quantile of the standard normal, but
we can replace this with other scaling without undue effort.

8Small-data sensitivity of MCMC analyses

positive, we takeϕ(w) = −(Ew g(β) − z0.975

p

Varw g(β)).

ϕ(w) > 0 is equivalent to moving the left endpoint below zero, thus changing from a significant result to a non-significant
one. Finally, to change to a significant result of the opposite sign, if the uncertainty interval’s left endpoint were positive, wetakepϕ(w) = −(Ew g(β) + z0.975 Varw g(β)).

On the full data, the right endpoint is above zero. On weight such that ϕ(w) > 0, the right endpoint has been moved below
zero: the conclusion has changed from a positive result to a significant negative result.
B.2. Regularity conditions
To be able to form a Taylor series, we require that the quantity of interest ϕ(w) is differentiable with respect to the weight w.
We are not aware of a complete theory (necessary and sufficient conditions) for this differentiability. However, throughAssumption B.1 and Assumption B.2, we state a set of sufficient conditions.
Assumption B.1. Let g be a function from RV to the real line. ϕ(w) is a linear combination of posterior mean and posterior
standard deviation i.e. there exists constants c1 and c2 , which are independent of w, such thatpϕ(w) = c1 Ew g(β) + c2 Varw g(β).
A typical choice of g is the function that returns the v-th coordinate of a V -dimensional vector.

It might appear that constraining ϕ(w) to be a linear combination of the posterior mean and standard deviation is overly
restrictive. However, this choice encompasses many cases of practical interest: recall from Appendix B.1 that the quantities
of interest for changing sign, changing significance, and producing a significant result of the opposite sign, take the formof Assumption B.1. Furthermore, the choice of constraining ϕ(w) to be a linear combination of the posterior mean andstandard deviation in Assumption B.1 is done out of convenience. Our framework can also handle quantities of interest thatinvolve higher moments of the posterior distribution, and the function that combines these moments need not be linear, but
we omit these cases for brevity. However, we note that posterior quantiles in general do not satisfy Assumption B.1 andleave to future work the question of how to diagnose the sensitivity of such quantities of interest.
Assumption B.2. For any w ∈ [0, 1]N \ {0N }, the following functions have finite expectations under the weighted posterior:|g(β)|, g(β)2 , |L(d(n) | β)| (for all n), |g(β)L(d(n) | β)| (for all n) and |g(β)2 L(d(n) | β)| (for all n).
The assumption is mild. It is satisfied by for instance, linear regression under Gaussian likelihood and g(β) = βv .
Under Assumption 2.1, Assumption B.1, and Assumption B.2, ϕ(w) is continuously differentiable with respect to w.
Theorem B.1. Assume Assumption 2.1, Assumption B.1, and Assumption B.2. For any δ ∈ (0, 1), ϕ(w) is continuously
differentiable with respect to w on {w ∈ [0, 1]N : maxn wn ≥ δ}. The n-th partial derivative2 at w is equal to c1 f + c2 s
where


f = Covw g(β), L(d(n) | β) ,
(3)
and



Covw g(β)2 , L(d(n) | β) − 2Ew g(β) × Covw g(β), L(d(n) | β)p.
s=Varw g(β)(4)

See Proof D.2 for the proof. This theorem is a specific instance of the sensitivity of posterior expectations with respect to log
likelihood perturbations: for further reading, we recommend Diaconis and Freedman (1986); Basu et al. (1996); Gustafson(1996). Theorem 3.1 establishes both the existence of the partial derivatives and their formula. Equation (3) is the partial
derivative of the posterior mean with respect to the weights, while Equation (4) is that for the posterior standard deviation,with the understanding that the derivative is one-sided.
B.3. Estimating the influence for other QoI
Algorithm 1 describes the algorithm to estimate the influence of a quantity of interest ϕ(w) that satisfies Assumption B.1
and Assumption B.2.
2

If wn lies on the boundary, the partial derivative is understood to be one-sided.

9

Small-data sensitivity of MCMC analyses

Algorithm 1 Influence Estimate (EI)
Input: ϕ(w)-defining constants c1 , c2 , Markov chain (β (1) , . . . , β (S) )PS
PS
m = S1 s=1 g(β (s) ), k = S1 s=1 g(β (s) )2
v = k − m2ψ̂ = (0, 0, . . . , 0) {N -dimensional vector}
for n = 1P
to N doS
a = S1 s=1 g(β (s) )L(d(n) | β (s) )PS
1b = S s=1 g(β (s) )2 L(d(n) | β (s) )PS
u = S1 s=1 L(d(n) | β (s) )f = a − mu {Estimate of Equation (3)}g = b − ku√ŝ = (g − 2mf )/( v) {Estimate of Equation (4)}ψ̂n = c1 f + c2 s {Estimate of ψn }end for
Return: ψ̂Algorithm 2 Sum of Sorted Influence Estimate (SoSIE)
Input: ϕ(w)-defining constants c1 , c2 , Markov chain (β (1) , . . . , β (S) ), fraction of data to drop αψ̂ = EI(c1 , c2 , (β (1) , . . . , β (S) ))Find ranks v1 , v2 , . . . , vN such that ψ̂v1 ≤ ψ̂v2 ≤ . . . ≤ ψ̂vN
Find the smallest p such that ψ̂vp+1 ≥ 0. If none exists, set p to N .
w
b is the N -vector where w
bn = 1 for n ∈ {v1 , . . . , vmin(p,⌊N α⌋) } and w
bn = 0 otherwiseb =∅b = {dv , . . . , dv}.
Otherwise,U
If p ≥ 1, U
1min(p,⌊N α⌋)b = − P⌊N α⌋ ψ̂v I{ψ̂v < 0}
∆
m
m
m=1
b U
b
Return: ∆,Algorithm 2 describes the algorithm to estimate the worst-case change in quantity of interest ϕ(w) that satisfies Assumption B.1 and Assumption B.2.
B.4. Confidence interval under exact sampling
For certain prior and likelihoods, we are able to draw exact Monte Carlo samples from the posterior distribution i.e.
(β (1) , . . . , β (S) ) is an i.i.d. sample of size S drawn from the full-data posterior distribution. This happens for conjugate
models (Diaconis and Ylvisaker, 1979), or for models in which convenient augmentation schemes have been discovered,b can besuch as Bayesian logistic regression with Polya-Gamma augmentation (Polson et al., 2013). Conceptually, ∆(n) Nthought of as an estimator constructed from an i.i.d. sample. However, the sample in question is not the data {d }n=1 , butb and (β (1) , . . . , β (S) ), we will use the notation ∆(βb (1) , . . . , β (S) ).
(β (1) , . . . , β (S) ). To highlight the dependence between ∆b is a complex, non-smooth function of the sample: the act of taking the minimum across the estimated
The estimator ∆influences ψ̂n is non-smooth. We do not attempt to prove distributional results for this estimator and use such results toquantify uncertainty. Instead, we appeal to the bootstrap (Efron, 1979), a general-purpose technique to quantify the sampling
uncertainty of estimators.
b
Our confidence interval construction proceeds in three steps. First, we define the so-called bootstrap distribution of ∆.
Second, we approximate this distribution with an empirical distribution based on Monte Carlo draws. Finally, we use therange spanned by quantiles of this empirical distribution as our confidence interval for ∆(α).
To define the bootstrap distribution, consider the empirical distribution of the sample (β (1) , . . . , β (S) ):S

1X
δ (i) (·).
S i=1 {β }We denote one draw from this empirical distribution by β ∗ (s) .
10

A bootstrap sample is a set of S draws:

Small-data sensitivity of MCMC analyses

b is the distribution of ∆(βb ∗ (1) , β ∗ (1) , . . . , β ∗ (S) ), where the(β ∗ (1) , β ∗ (1) , . . . , β ∗ (S) ). The bootstrap distribution of ∆randomness is taken over the bootstrap sample but is conditional on the original sample (β (1) , . . . , β (S) )Clearly, the bootstrap distribution is discrete with finite support. If we chose to, we can enumerate its support and compute its
probability mass function, by enumerating all possible values a bootstrap sample can take. However, this is time-consuming.
b ∗ (1) , β ∗ (1) , . . . , β ∗ (S) ) is
It suffices to approximate the bootstrap distribution with Monte Carlo draws. The draw ∆(βb ∗ : we generate a total number of B such draws. When B increases, the empirical distribution ofabbreviated by ∆b ∗, ∆b ∗, . . . , ∆b ∗ ) becomes a better approximation of the bootstrap distribution. However, the computational cost scales up(∆
1
2
Bwith B. In practice, B in the hundreds are commonplace: our numerical work uses B = 200.
We now define confidence intervals for ∆(α). Each interval is parametrized by η, the nominal coverage level, whichb ∗, . . . , ∆b ∗, ∆b ∗ ), the (1 − η)/2 andis valued in (0, 1). We compute two quantiles of the empirical distribution over (∆2
1B
3(1 + η)/2 quantiles , and define the interval spanned by these two values as our confidence interval. By default, we set
η = 0.95.

One limitation of our current work is that we do not make theoretical claims regarding the actual coverage of such confidence
intervals. Although bootstrap confidence intervals can always be computed, whether the actual coverage matches thenominal coverage η depends on structural properties of the estimator and regularity conditions on the sample. To verify
the quality of these confidence intervals, we turn to numerical simulation. We leave to future work the task of formulating
reasonable assumptions and theoretically analyzing the actual coverage.
B.5. Confidence interval under general MCMC
In the previous section, we made the simplifying assumption that exact sampling were possible. We now lift this assumption
and handle the case in which (β (1) , . . . , β (S) ) truly came from a Markov chain (such as the output of Hamiltonian Monte
Carlo). This case is much more common in practice than the exact sampling case.
To construct confidence intervals, one idea is to use the previous section’s construction without modification. In other words,
apply the bootstrap to a non-i.i.d. sample: recall that the Markov chain states are not independent of each other. Theoretically,
it is known that the bootstrap struggles on non-i.i.d. samples, for even simple estimators. For example, if the estimator inquestion is the sample mean and the draws exhibit positive autocorrelation, under mild regularity conditions, the bootstrap
variance estimate seriously underestimates the true sampling variance, even in the limit of infinite sample size (Lahiri, 2003,
b for
Theorem 2.2). In our case, the bootstrap likely struggles on the sample means that are involved in the definition of ∆:
(1)
(2)
(S)
instance, it is very common for some coordinate v that (βv , βv , . . . , βv ) exhibits positive autocorrelation in practice.
Therefore, we have reason to be pessimistic about the ability of bootstrap confidence intervals to adequately cover ∆(α).
Fundamentally, the bootstrap fails in the non-i.i.d. case because the draws that form the bootstrap sample do not have any
dependence, while the draws that form the original sample do. To improve upon the bootstrap, one option is to resample in a
way that respects the original sample’s dependence structure. We recognize that the sample in question, (β (1) , . . . , β (S) ), isa (multivariate) time series: we focus on methods that perform well under time series dependence. One such scheme is thenon-overlapping block bootstrap (Lahiri, 2003; Carlstein, 1986).4 The sample (β (1) , . . . , β (S) ) is divided up into a numberof blocks: each block is a vector of contiguous draws. Let L be the number of elements in a block, and let M := ⌊S/L⌋
denote the number of blocks. The m-th block is defined as
Bm := β ((m−1)L+1) , . . . , β (mL) .
To generate one sample from the non-overlapping block bootstrap distribution, we first draw with replacement from the set∗of blocks M values: B1∗ , . . . , BM
. Then, we write the elements of these drawn blocks in a contiguous series. For example,
(1)
(S)
(1)
when (β , . . . , β ) = (β , β (2) , β (3) , β (4) ) and L = 2, the two blocks are (β (1) , β (2) ), and (β (3) , β (4) ). The set of
possible samples from resampling include (β (1) , β (2) , β (1) , β (2) ) and (β (3) , β (4) , β (3) , β (4) ) but not (β (1) , β (3) , β (1) , β (3) ).
The name “non-overlapping block bootstrap” comes from the fact that these blocks, viewed as sets, are disjoint from eachother. While the name is needed in Lahiri (2003) to distinguish from other blocking rules, moving forward, as we only3
We use R’s quantile() to compute the sample quantiles. When (1 + η)/2 × B is not an integer, the (1 + η)/2 quantile is defined bylinearly interpolating the order statistics.
4
The original paper, Carlstein (1986), did not use the term “non-overlapping block bootstrap” to describe the technique. The namecomes from Lahiri (2003).

11

Small-data sensitivity of MCMC analyses

consider the above blocking rule, we will refer to the procedure as simply, block bootstrap. Intuitively, the block bootstrap
sample is a good approximation of the original sample if the latter has short-term dependence: in such a case, the original
sample itself can be thought of as the concatenation of smaller, i.i.d. subsamples, and the generation of a block bootstrap
sample mimics that. In well-behaved probabilistic models with well-tuned algorithms, the MCMC draws can be expected toonly have short-term dependence, and the block bootstrap is a good choice.
The block bootstrap has one hyperparameter: the block length L. We would like both L and M to be large: large L captures
time series dependence at larger lags, and large M is close to having many i.i.d. subsamples. However, since their product isconstrained to be S, the choice of L is a trade-off. In numerical studies, we set L = 10.
Our construction of confidence intervals for general MCMC proceeds identically to the previous section’s construction,
except for the step of generating the bootstrap sample: instead of drawing from the vanilla bootstrap, we draw from theblock bootstrap. We will denote the endpoints of such an interval by ∆lb (α) (lower endpoint) and ∆ub (α) (upper endpoint).
Similar to the previous section, we do not make theoretical claims on the actual coverage of our block bootstrap confidence
intervals: we verify the quality of the intervals through later numerical studies.

C. Theory
In this section, we theoretically quantify the approximation errors incurred by our methods. Namely, Appendix C.1 analyzes
the error made by the first-order approximation, while Appendix C.2 analyzes the error made by using MCMC to estimate
influences.
C.1. Accuracy of first-order approximation
In this section, we investigate the error incurred by replacing ϕ(w) − ϕ(1N ) with the Taylor series from Section 3.1. While
the approximation applies to any model that satisfies Assumption 2.1, Assumption B.1, and Assumption B.2, our error
analysis is limited to two models: a normal model and a normal means model. Their salient features are the following.
Both are convenient to analyze and address the same statistical task: derive the population mean based on a finite sample(n){x(n) }N
∈ R. The normal means model is hierarchical: the observations are organized into disjoint
n=1 where xgroups. Each observation d(n) is (x(n) , g (n) ), where g (n) is valued in {1, 2, . . . , G}, with g (n) = g indicating that the n-th
observation belongs to the g-th group. The normal model does not have this structure, as only x(n) is observed and usedin modeling. We show that error in the normal model is qualitatively different from the error in the normal means model.
Roughly speaking, the former depends on the ratio between the number of observations left out, ⌊N α⌋, and the total numberof observations, N . Meanwhile, the later depends on three quantities a) ⌊N α⌋, b) the number of groups G and c) the numberof observations in a group.
Before specializing to different
P models, we pin down the common notion of error. We define error to be the difference
between ϕ(w) − ϕ(1N ) and n (wn − 1)ψn . We mainly care when w encodes the full removal of certain observations
and full inclusion of the remaining ones i.e. w ∈ {0, 1}N . If we let q be the function that returns the zero indices of such aweight (q(w) = {n : wn = 0}), then its inverse q −1 takes a set of observation indices (I ⊂ {1, 2, . . . , N }) and producesa weight valued in {0, 1}N . We reformulate the error as a function of I instead of w by replacing w with q −1 (I) in thedefinition of error. This reformulation reads
Err1st (I) = ϕ(q −1 (I)) − ϕ(1N ) +

X

ψn .

n∈I

C.1.1. N ORMAL MODEL .
We detail the prior and likelihood of the normal model and the associated quantity of interest. The parameter of interest is
the population mean µ. The likelihood of an observation is Gaussiana known standard deviation σ. In other words, the with1n-th log-likelihood evaluated at µ is L(d(n) | µ) = 21 log 2πσ− 2σ1 2 [(x(n) )2 − 2x(n) µ + µ2 ]. We choose the uniform
2
distribution over the real line as the prior for µ. The quantity of interest is the posterior mean of µ.
In this model, expectations under the weighted posterior have closed forms. We can derive an explicit expression for
the error. To display the error, it is convenient
to define the sample average of observations as a function of I: for any
PI ⊂ {1, 2, . . . , N }, let x̄I := (1/|I|) n∈I x(n) . The sample average of the whole dataset will be denote by x̄.
12Small-data sensitivity of MCMC analyses

Lemma C.1. For the normal model, Err1st (I) is qual to|I|2 (x̄ − x̄I )N (N − |I|)
We prove Lemma C.1 in Proof D.3. The error is a function of I through the a) the cardinality of the set |I| and b) the
difference between the whole dataset’s sample mean, x̄, and the sample mean for elements in I. Since the data is fixed, we
can upper bound |x̄ − x̄I | with the constant 2∥x∥∞ , where ∥x∥∞ := maxn |x(n) |. The rate at which the absolute value ofthe error goes to zero is |I|2 : as the ratio |I|/N equals α, this means the error’s absolute value goes to zero like α2 .
C.1.2. N ORMAL MEANS MODEL .
We detail the prior and likelihood of the normal means model and the associated quantity of interest. The parameters ofinterest are the population mean µ and the group means θ = (θ1 , θ2 , . . . , θG ). Observations in group g are modeled as
Gaussian centered at the groupθg with a known standard deviation σ. In other words, the n-th log-likelihood is mean
1
1(n) 2
L(d(n) | µ, θ) = 12 log 2πσ−[(x) − 2x(n) θg(n) + θg2(n) ]. The prior over (µ, θ) is the following. We choose the2
2σ 2
uniform distribution over the real line as the prior for µ. Conditioned on µ, the group means are Gaussian centered at µ,with a known standard deviation τ . The quantity of interest is the posterior mean of µ.
This model, like the normal model, has closed-form posterior expectations. Before displaying the exact formula for the error
Err1st (I), we need to describe the weighted posterior in more detail. For each group g, we define three functions of w:
Ng (w) :=X

wn , Mg (w) :=n:g (n) =g

P

n:g (n) =g wn x(n)Ng (w), Λg (w) :=



σ2
+ τ2
Ng (w)

−1.

While Ng (w) sums up the weights of observations in group g, Mg (w) is the weighted average of observations in this group,and Λg (w) will be used to weigh Mg (w) in forming the posterior mean of µ. Proof D.4 shows that Ew µ is equal toPG

g=1 Λg (w)Mg (w).
PG
g=1 Λg (w)PG
PG
To avoid writing g=1 Λg (w), we define Λ(w) := g=1 Λg (w). To lighten notation, for expectations under the original
posterior, we write µ∗ instead of E1N µ and Ng∗ instead of Ng (1N ). The same shorthand applies to Ng (1N ), Mg (1N ),Λg (1N ) and Λ(1N ). In words, µ∗ is the posterior mean of µ under the full-data posterior, Ng∗ is the number of observations
in group g of the original dataset, and so on. We also utilize the x̄I and x̄ notations defined the normal model section.
The error in the normal means model is given in the following lemma.
Lemma C.2. In the normal means model, let the index set I be such that there exists k ∈ {1, 2, . . . , G} such that g (n) = k
for all n ∈ I. Define|I|2|I| σ 2 Λ∗k ∗F (I) := ∗ ∗(Mk∗ − x̄I ) + ∗(µ − Mk∗ ),Nk [Nk − |I|]
Nk Nk∗|I|
E(I) := ∗ ∗σ 2 Λk (q −1 (I))Λ∗k .
Nk [Nk − |I|]
Then, Err1st (I) is equal to−1

Λk (q (I))
F (I) +
Λ∗

P

∗
∗
−1
(I)))
g̸=k Λg (Mg − Mk (qΛ∗ Λ(q −1 (I))



E(I).

We prove Lemma C.2 in Proof D.4. The constraint where all observations in I belong to the same group k is made out of
convenience: we can derive the error without this constraint, but the formula will be much more complicated.
A corrolary of Lemma C.2 is that the absolute value of the error behaves like |I|2 /(G|Nk∗ |2 ).
13

Small-data sensitivity of MCMC analyses

Corollary C.1. In the normal means model, for all groups g, assume that Ng∗ ≥ σ 2 /τ 2 . Let the index set I be such that
there exists k ∈ {1, 2, . . . , G} such that g (n) = k for all n ∈ I. For this k, assume that Nk∗ − |I| ≥ σ 2 /τ 2 . Then,|Err1st (I)| ≤ C(∥x∥∞ , σ, τ )1 |I|2.
G |Nk∗ |2

where C(∥x∥∞ , σ, τ ) is a constant that only depends on ∥x∥∞ , σ, and τ .
We prove Corollary C.1 in the Proof D.5. In addition to the assumptions Lemma C.2, the corrolary assumes that the numberof observations in each group is not too small, and that after removing I, group k still has enough observations. Thiscondition allows us to approximate Λ∗k and Λg (q −1 (I)) with a constant. The factor ∥x∥∞ in the bound comes from upper(n)bounding |Mg∗ − Mk (q −1 (I))| by 2 maxN|.
n=1 |xFor two reasons, we conjecture that similiar qualitative differences also appear in the comparison between more complicatedhierarchical and non-hierarchical models The fundamental task of estimating the population mean is embedded in manyother statistical tasks, such as regression. In addition, the group structure imposed by the normal means model is also foundin practically relevant hierarchical models.
C.2. Estimator properties
b is the (β (1) , . . . , β (S) )-induced sampling uncerRecall from Section 3.3 that one concern regarding the quality of ∆b is a non-smooth function oftainty. Theoretically analyzing this uncertainty is difficult, with one obstacle being that ∆(1)(S)(β , . . . , β ). In this section, we settle for the easier goal of analyzing the sampling uncertainty of the influence estimatesb but weψ̂n . We expect such theoretical characterizations to play a role in the eventual theoretical characterizations of ∆,leave this step to future work.

In this analysis, we make more restrictive assumptions than those needed for Theorem 3.1 to hold. We assume that the sample(β (1) , . . . , β (S) ) comes from exact sampling: the independence across draws makes it easier to analyze sampling uncertainty.
We focus on the quantity of interest equaling the posterior mean (c1 = 1, c2 = 0 in the sense of Assumption B.1): thescaling c1 = 1 for the posterior mean is made out of convenience, and a smiliar analysis can be conducted when c2 ̸= 0, but
we omit it for brevity. Finally, we need more stringent moment conditions than Assumption B.2.
Assumption C.1. The functions |g(β)2 L(d(i) | β)L(d(j) | β)| (across i, j) have finite expectation under the full-data
posterior.
This moment condition guarantees that the sample covariance of g(β) and L(d(i) | β) has finite variance under the full-data
posterior: it plays the same as role as finite kurtosis in proofs sample variance consistency
With the assumptions in place, we begin by showing that the sampling uncertainty of ψ̂n goes to zero in the limit of S → ∞.

Lemma C.3. Assume Assumption 2.1, Assumption B.1, Assumption B.2, Assumption C.1 holds. Let ψ̂ be output of
Algorithm 1 for c1 = 1, c2 = 0 and (β (1) , . . . , β (S) ) being an i.i.d. sample. Then, there exists a constant C such that for alln, for all S, Var(ψ̂n ) ≤ C/S.
We prove Lemma C.3 in Proof D.7. That the variance of individual ψ̂n goes to zero at the rate of 1/S is not surprising: ψ̂n
is a sample covariance, after all.
We use Lemma C.3 to show consistency of different estimators.

Theorem C.1. Assume Assumption 2.1, Assumption B.1, Assumption B.2, and Assumption C.1 holds Let ψ̂ be outputof Algorithm 1 for c1 = 1, c2 = 0 and (β (1) , . . . , β (S) ) being an i.i.d. sample. Then maxN
n=1 |ψ̂n − ψn | converges inb converges in probability to ∆(α) in the limit S → ∞.
probability to 0 in the limit S → ∞, and ∆b is
We prove Theorem C.1 in Proof D.8. Our theorem states that the vector ψ̂ is a consistent estimator for the vector ψ and ∆a consistent estimator for ∆(α).
Not only is ψ̂ consistent in estimating ψ, it is also asymptotically normal.

Theorem C.2. Assume Assumption 2.1, Assumption B.1, Assumption B.2, and Assumption C.1 holds.
Let ψ̂ be√
output of Algorithm 1 for c1 = 1, c2 = 0 and (β (1) , . . . , β (S) ) being an i.i.d. sample. Then S(ψ̂ − ψ) converges in distribution to N (0N , Σ) where Σ is the N × N matrix whose (i, j) entry, Σi,j , is the covariance between
14

Small-data sensitivity of MCMC analyses(g(β) − E1N g(β)) L(d(i) | β) − E1N L(d(i) | β) and (g(β) − E1N g(β)) L(d(j) | β) − E1N L(d(j) | β) , taken under
the full-data posterior.
We prove Theoremp
C.2 in Proof√ D.9. Heuristically, for each n, the distribution of ψ̂n is the Gaussian centered at ψn , withstandard deviation Σn,n / S.

C.2.1. N ORMAL MODEL WITH UNKNOWN PRECISION .
√p
While Σn,n / S eventually goes to zero, for finite S, this standard deviation can be large, making ψ̂n an imprecise
estimate of ψn . To illustrate this phenomenon, we will derive Σn,n in the context of a simple probabilistic model: a normal
model with unknown precision.
We first introduce the model and the associated quantity of interest. The data is a set of N real values: d(n) = x(n) , wherex(n) ∈ R. The parameters of interest are the mean
precision τ of the population. The log-likelihood of an µ 1and theτ− 2 τ [(x(n) )2 − 2x(n) µ + µ2 ]. The prior is chosen to be the following.
observation based on µ and τ is Gaussian: 12 log 2πµ is distributed from uniform over the real line, and τ is distributed from a gamma distribution. The quantity of interest is
the posterior mean of µ.
For this probabilistic model, the assumptions of Theorem C.2 are satisfied. We show that the variance Σn,n behaves like aquartic function of the observation x(n) .
Lemma C.4. In the normal-gamma model, there exists constants D1 , D2 , and D3 , where D1 > 0, such that for all n, Σn,n
is equal to D1 (x(n) − x̄)4 + D2 (x(n) − x̄)2 + D3 .
We prove Lemma C.4 in Proof D.10. D1 , D2 , D3 are based on the posterior expectations: for instance, the proof shows that
D1 =E1N [τ −1 (τ −E1N τ )2 ]. It is easy to show that for the normal-gamma model,
4N

Cov1N (µ, L(d(n) | µ, τ )) =x(n) − x̄
.
NHence, while the mean of ψ̂n behaves like a linear function of x(n) − x̄, its standard deviation behaves like a quadratic
function of x(n) − x̄. In other words, the more influence an observation has, the harder it is to accurately determine its
influence!

D. Proofs
D.1. Taylor series proofs
Proof D.1 (Proof of Theorem 3.1). Theorem 3.1 is a special case of Theorem B.1. See Proof D.2 for the proof of the latter,
which implies the proof of the former.
Proof D.2 (Proof of Theorem B.1). At a high level, we rely on Fleming (1977, Chapter 5.12, Theorem 5.9) to interchange
integration and differentiation.
Although the theorem statement does not explicitly mention the normalizer, to show that the quantity of interest is continuously
differentiable and compute partial derivatives, it is necessary to show that the normalizer is continuously differentiable andcompute partial derivatives. To do so, we verify the following conditions on the integrand defining Z(w):1. For any β, the mapping w 7→ p(β) exp

P

N
(n)
| β)
n=1 wn L(dis continously differentiable.

2. There exists
function f1 such that for all w ∈ {w ∈ [0, 1]N : maxn wn ≥ δ},P a Lebesgue integrableN(n)| β) ≤ f1 (β).
p(β) exp
n=1 wn L(d3. For each n, there
function f2 such that for all w ∈ {w ∈ [0, 1]N : maxn wn ≥ δ},
Pexists a Lebesgue integrable

N
∂
(n)
| β) ≤ f2 (β).
n=1 wn L(d∂wn p(β) exp

The first condition is clearly satisfied. To construct f1 that satisfies the second condition, we partition the parameter space
RV into a finite number of disjoint sets. To index these sets, we use a subset of {1, 2, . . . , N }. If the indexing subset were15

Small-data sensitivity of MCMC analyses

I = {n1 , n2 , . . . , nM }, the corresponding element of the partition is
BI := {β ∈ RV : ∀n ∈ I, L(d(n) | β) ≥ 0}.

(5)This partition
PNallows us to upper bound the integrand with a function that is independent of w. Suppose β ∈ BI , I ̸= ∅. The
maximum n=1 wn L(d(n) | β) is attained by setting wn = 1 for all n ∈ I and wn = 0 for all n ∈/ I. Suppose β ∈ B∅ . As
PN(n)
L(d | β) < 0 for all 1 ≤ n ≤ N , and we are constrained by maxn wn ≥ δ, the maximum of n=1 wn L(d(n) | β) is
attained by setting wn = δ for arg maxn L(d(n) | β) and wn = 0 for all other n. In short, our envelope function is(Q
if β ∈ BI , I ̸= ∅.
p(β) n∈I exp(L(d(n) | β))f1 (β) :=(n)p(β) maxN
exp(δL(d|β))
ifβ ∈ B∅ .
n=1

The last step is to show f1 is integrable. It suffices to show that the integral of f1 on each BI is finite. On B∅ , integratingp(β) exp(δL(d(n) | β)) over B∅ is clearly finite: by Assumption 2.1, the integral of p(β) exp(δL(d(n) | β)) over RV
is finite, and B∅ is a subset of RV . As f1 (β) is the maximum of a finite number
Q of integrable functions, it is integrable.
Similarly, the integral of f1 over BI where I ̸= ∅ is atmost the integral of p(β) n∈I exp(L(d(n) | β)) over RV , which isfinite by Assumption 2.1. To construct f2 that satisfies the third condition, we use the same partition of RV , and the envelope
function is f2 (β) := L(d(n) | β)f1 (β), since the partial derivative of the weighted log probability is clearly the product of
the n-th log likelihood and the weighted log probability. The integrability of f2 follows from Assumption B.2’s guaranteethat the expectation of |L(d(n) | β)| is finite under different weighted posteriors. In all, we can interchange integration withdifferentiation, and the partial derivatives are
h
i
∂Z(w)
= Z(w) × Ew L(d(n) | β) .
∂wn
We move on to provedifferentiable and find its partial derivatives. The conditions onPthat Ew g(β) is continuously

N
1
(n)
g(β) Z(w) p(β) exp| β) that we wish to check are:n=1 wn L(d1p(β) exp
1. For any β, the mapping w 7→ g(β) Z(w)

P

N
(n)
| β)
n=1 wn L(dis continously differentiable.

f1 such that for all w ∈ {w ∈ [0, 1]N : maxn wn ≥ δ},
2. There exists a LebesgueP integrable functionN1(n)| β) ≤ f3 (β).
g(β) Z(w) p(β) exp
n=1 wn L(d3. For each n, there exists a Lebesgue integrablefunction f4 such that for all w ∈ {w ∈ [0, 1]N : maxn wn ≥ δ},
PN
∂
1
(n)
| β) ≤ f4 (β).
n=1 wn L(d∂wn g(β) Z(w) p(β) exp

We have already proven that Z(w) is continuously differentiable: hence, there is nothing to do for the first condition.
It is straightforward to use Assumption B.2 and check that the second condition is satisfied by the function f3 (β) :=11(n)| β)f1 (β). Hence, we can interchangeZ(w) g(β)f1 (β), and the third condition is satisfied by f4 (β) := Z(w) g(β)L(dintegration with differentiation. The partial derivatives of Ew g(β) is equal to tthe sume of two integrals. The first part is!!
Z
NX∂Z(w)−1(n)wn L(d | β)dβg(β)p(β) exp
∂wn
n=1
!!
N
 h
i Z
X
1
(n)
(n)
= − Ew L(d | β)wn L(d | β)dβg(β)p(β) exp
Z(w)
n=1
h
i
= −Ew L(d(n) | β) × Ew [g(β)] .
The second part isZ

NX1
wn L(d(n) | β)g(β)L(d(n) | β)p(β) exp
Z(w)
n=1

16

!!

h
i
dβ = Ew g(β)L(d(n) | β) .

Small-data sensitivity of MCMC analyses

Putting two parts together, the partial derivative is equal to a covariancehi∂Ew g(β)= Covw g(β), L(d(n) | β) .
∂wn
The proof that Ew g(β)2 is continuously differentiable is similar to that for Ew g(β). The partial derivative ishi∂[Ew g(β)2 ]= Covw g(β)2 , L(d(n) | β) .
∂wn
Since the posterior standard deviation is a differentiably continuous function of the mean and second moment, it is also
differentiably continuous. The partial derivative of the posterior standard deviation is a simple application of the chain rule,and we omit the proof for brevity.
D.2. First-order accuracy proofs
Proof D.3 (Proof of Lemma C.1). Our proof finds exact formulas for the posterior mean and the partial derivatives of theposterior mean with respect to wn . Then, we take the difference between the posterior mean and its Taylor series.
In the normal model, the total log probability at w is equal to



1
1
1
(n) 2
(n)
2
log
− 2 [(x ) − 2x µ + µ ]
wn
2
2πσ 2
2σ
n=1
!2
!
PN
PN
(n)
n=1 wnn=1 wn x=−+ C,µ − PN2σ 2
n=1 wn


N
X

where
PN C is a constant
PN on µ. Hence, the distribution of µ under w is normal with mean
PN that does not depend( n=1 wn x(n) )/( n=1 wn ) and precision ( n=1 wn )/(σ 2 ). The partial derivative of the posterior mean with respect to
wn is
PN
PN
x(n) ( n=1 wn ) − ( n=1 wn x(n) ).
PN( n=1 wn )2
Plugging in w = 1N , we have that ψn is equal to (x(n) − x̄)/N .

After removing the index set I, the actual posterior mean isN x̄ − |I|x̄I,
N − |I|
while the Taylor series approximation isx̄ −X x(n) − x̄
n∈I

N

=

N x̄ + |I|(x̄ − x̄I ).
NThe difference between the actual posterior mean and its approximation is as in the statement of the lemma.
Proof D.4 (Proof of Lemma C.2). Similar to the proof of Lemma C.1, we first find exact formulas for the posterior mean
and its Taylor series.
In the normal means model, the total log probability at w is
G 
X

log

g=1

+



1
2πτ 2



−

1
(θg − µ)2
2τ 2




 
i
1 h (n) 2
1
(n)
2
(x
)
−
2x
θ
+
θ
.
−
wn log
(n)
(n)
g
g
2πσ 2
2σ 2
n=1
N
X

By completing the squares, we know that17

Small-data sensitivity of MCMC analyses• The distribution of µ is normal:
NPG

1
g=1 Λg (w)Mg (w), PG
PG
g=1 Λg (w)g=1 Λg (w)!

• Condition on µ, the group means are independent normals:µ/τ 2 + [Ng (w)Mg (w)]/σ 2
1
θg | µ ∼ N,.
1/τ 2 + Ng (w)/σ 2
1/τ 2 + Ng (w)/σ 2
To express the partial derivative of the posterior mean of µ with respect to wn , it is helpful to define the following“intermediate” value between Ew µ and Ew θg :µ̃g (w) :=Mg (w)Ng (w)/σ 2 + Ew µ/τ 2.
Ng (w)/σ 2 + 1/τ 2

In addition, we need the partial derivatives of the functions Ng , Λg , and Mg .
(0 if g ̸= g (n)
∂Ng
=
,
∂wn
1 if g = g (n)(0if g ̸= g (n)∂Mg,= x(n) −Mg (w)∂wn
if g = g (n)Ng (w)(
0
if g ̸= g (n)
∂Λg
2
=
.
Λ
(w)
g
∂wn
σ 2 N (w)2 if g = g (n)g

If n is in the k-th group, the partial derivative of the posterior mean with respect to wn is
1
1
x(n) − µ̃k (w) .
22
Λ(w) σ + τ Nk (w)After removing only observations from the k-th group, the actual posterior mean isPΛk (q −1 (I))Mk (q −1 (I)) + g̸=k Λg (1N )Mg (1N )P.
Λk (q −1 (I)) + g̸=k Λg (1N )Between w = q −1 (I) and w = 1N , the Ng , Mg , Λg functions do not change for g ̸= k. The Taylor series approximation of
the posterior mean is

 P

P
Λk (1N ) Mk (1N ) + n∈I µ̃k (1N ) − x(n) /Nk (1N ) + g̸=k Λg (1N )Mg (1N )P.
Λk (1N ) + g̸=k Λg (1N )If we denote

A1 :=

X

Λg (1N )Mg (1N ), A2 :=XΛg (1N )

g̸=k

g̸=k

B1 := Λk (q −1 (I))Mk (q −1 (I)), B2 := Λk (q −1 (I))
,
#
"

X
µ̃k (1N ) − x(n) /Nk (1N ) , C2 := Λk (1N )C1 := Λk (1N ) Mk (1N ) +n∈Ithen Err1st (I) is equal to (A1 + B1 )/(A2 + B2 ) − (A1 + C1 )/(A2 + C2 ). The last equation is equal toA2 (B1 − C1 ) + A1 (C2 − B2 ) + (B1 C2 − C1 B2 ).
(A2 + B2 )(A2 + C2 )18

Small-data sensitivity of MCMC analyses

We analyze the differences C2 − B2 , B1 C2 − C1 B2 , and B1 − C1 separately.
C2 − B2 . This difference is

1
1
− 2
.
σ 2 /Nk (1N ) + τ 2
σ /Nk (q −1 (I)) + τ 2

Since we remove |I| from group k, Nk (q −1 (I)) = Nk (1N ) − |I|. Hence, the difference C2 − B2 is|I|,
Nk (1N )(Nk (1N ) − |I|)σ 2 Λk (1N )Λk (q −1 (I))
which is exactly the E(I) mentiond in the lemma statement.
B1 C2 − C1 B2 . The difference isΛk (1N )Λk (q−1((I)) Mk (q−1(I)) − Mk (1N ) −Pn∈I [µ̃k (1N ) − x(n)Nk (1N )]).

We analyze the term in the curly brackets. It is equal to
)
(
P
(n)
X  Mk (1N ) − µ̃k (1N ) 
[M
(1
)
−
x
]
k
N
n∈I
−1
+
Mk (q (I)) − Mk (1N ) −Nk (1N )Nk (1N )n∈I

The left term is equal to|I|2 (Mk (1N ) − x̄I ).
Nk (1N )[Nk (1N ) − |I|]
The right term is equal to|I| σ 2 Λk (1N )(E1N µ − Mk (1N )).
Nk (1N ) Nk (1N )The sum of the two terms is exactly F (I) mentioned in the lemma statement. Overall, the difference B1 C2 − C1 B2 is equalto Λk (1N )Λk (q −1 (I))F (I).
B1 − C1 . If we introduce D := Λk (1N )Mk (q −1 (I)), then the difference B1 − C1 is equal to (B1 − D) + (D − C1 ). The
former term is
Mk (q −1 (I))(B2 − C2 ) = −Mk (q −1 (I))E(I).
The later term is(Λk (1N ) Mk (q−1(I)) − Mk (1N ) −Pn∈I [µ̃k (1N ) − x

Nk (1N )

(n)

]

)

.

We already know that the term in the curly brackets is equal to F (I). Hence B1 − C1 is equal to Λk (1N )F (I) −Mk (q −1 (I))E(I).
With the differences C2 − B2 , B1 C2 − C1 B2 , and B1 − C1 , we can now state the final form of Err1st (I). The final numeratorisXΛk (q −1 (I)) +Λg (1N ) Λk (1N )F (I)+Divide this by the denominator

g̸=k

X

g̸=k

hP

Λg (1N )Mg (1N ) − Mk (q −1 (I))

g Λg (1N )

i hP

X

g̸=k



.

Λg (1N ) E(I)

i
−1
Λ
(q
(I))
, we have proven the lemma.
gg

Proof D.5 (Proof of Corollary C.1). Under the assumption that Ng∗ ≥ σ 2 /τ 2 , we have that Λg (1N ) ∈Mk∗ − |I| ≥ σ 2 /τ 2 , it is also true that Λk (q −1 (I)) ∈ 2τ1 2 , τ12 .
19

 1

1
2τ 2 , τ 2

. SinceSmall-data sensitivity of MCMC analyses

Because of Lemma C.2, an upper bound on Err1st (I) isΛk (q −1 (I))|F (I)| +
Λ∗

P

∗
∗
−1
(I)))
g̸=k Λg (Mg − Mk (qΛ∗ Λ(q −1 (I))



|E(I)| .

The fraction Λk (q −1 (I))/Λ∗ is at most ( τ12 )/ G 2τ12 , which is equal to 2/G. The absolute value |F (I)| is at most2|I|2 ∥x∥∞2|I|∥x∥∞ (σ 2 /τ 2 )2|I|2 ∥x∥∞ (σ 2 /τ 2 + 1)+≤∗∗
(Nk )2(Nk )2(Nk∗ )2

The absolute value

P

∗
∗
−1
(I)))
g̸=k Λg (Mg − Mk (qΛ∗ Λ(q −1 (I))

is at most4∥x∥∞G(1/τ 2 )2∥x∥∞≤.
G2 (1/2τ 2 )G
Finally, the absolute value |E(I)| is at most|I|2 (σ 2 /(4τ 4 ))|I|(σ 2 /(4τ 4 ))≤.
(Nk∗ )2(Nk∗ )2
In all, the constant C(∥x∥∞ , σ, τ ) in the corollary’s statement is∥x∥∞ 4(σ 2 /τ 2 + 1) + σ 2 /τ 4 .
D.3. Consistency and asymptotic normality proofs
The following lemma, on covariance between sample covariances under i.i.d. sampling, will be useful for later proofs.
Lemma D.1. Suppose we have S i.i.d. draws (A(s) , B (s) , C (s) )Ss=1 . Let f1 be the (biased) sample covariance between the
A’s and the B’s. Let f2 be the (biased) sample covariance between the A’s and C’s. In other words,
S

f1 :=

1 X (s) (s)
A B
S s=1
S

f2 :=

1 X (s) (s)
A C
S s=1

!

!

S

−

1 X (s)
A
S s=1
S

−

1 X (s)
A
S s=1

!

!

S

1 X (s)B
S s=1S

1 X (s)
C
S s=1

!

!

,
.

Suppose that the following are finite: E[(A − E[A])2 (B − E[B])(C − E[C])], Cov(B, C), Var(A), Cov(A, B), Cov(A, C).
Then, the covariance of f1 and f2 is equal to(S − 1)2
E[(A − E[A])2 (B − E[B])(C − E[C])]S3
S−1(S − 1)(S − 2)+
Cov(B, C)Var(A) −Cov(A, B)Cov(A, C).
3
S
S3
Proof D.6 (Proof of Lemma D.1). It suffices to prove the lemma in the case where E[A] = E[B] = E[C] = 0. Otherwise,
we can subtract the population mean from the random variable: the value of f1 and f2 would not change (since covarianceis invariant to constant additive changes). In otherwords, we want to show that the covariance between f1 and f2 is equal toS−1(S − 1)(S − 2)
(S − 1)2
E[A2 BC] +E[BC]E[A2 ] −E[AB]E[AC].
3
3
S
S
S3
20

(6)

Small-data sensitivity of MCMC analyses
S−1
Since f1 is the biased sample covariance, Ef1 = S−1
S E[AB]. Similarly, Ef2 = S E[AC]. To compute Cov(f1 , f2 ), weonly need an expression for E[f1 f2 ]. The product f1 f2 is equal to the sum of D1 , D2 , D3 , D4 where:!
!
!
1 X (s) (s)1 X (s)1 X (s)D1 := −
A B
A
C
,
S s
S s
S s
!2
!
!
1 X (s)1 X (s)1 X (s)
A
B
C
,
D2 :=
S s
S s
S s
!
!
!
1 X (s)1 X (s)1 X (s) (s)
,
A C
A
B
D3 := −
S s
S s
S s
!
!
1 X (s) (s)1 X (s) (s)
D4 :=
A B
A C
.
S s
S s

We compute the expectation of each Dj .
D1 . By expanding D1 , we know that ED1 = S13
on the triplet (i, j, k) in the following way:P(k) (k) (i) (j)B A C ]. The value of E[A(k) B (k) A(i) C (j) ] depends
i,j,k E[A



0



E[A2 BC](k) (k) (i) (j)E[A B A C ] = 0E[AB]E[AC]



0

if i = k, j ̸= k
if i = k, j = k
if i ̸= k, j = k
if i ̸= k, j ̸= k, i = j
if i ̸= k, j ̸= k, i ̸= j

We have used independence of (A(s) , B (s) , C (s) )Ss=1 to factorize the expectation E[A(k) B (k) A(i) C (j) ]. For certain triplets,
the factorization reveals that the expectation is zero. By accounting for all triplets, the expectation of D1 is
1 
SE[A2 BC] + S(S − 1)E[AB]E[AC] .
3SP
D2 . By expanding D2 , we know that ED2 = S14 i,j,p,q E[A(i) A(i) B (p) C (q) ]. We can do a similar case-by-case analysis
of how E[A(i) A(i) B (p) C (q) ] depend on the quartet (i, j, p, q). In the end, the expectation of D2 is
1 
E[A2 BC] + (S − 1)E[A2 ]E[BC] + 2(S − 1)E[AB]E[AC] .
S3
D3 . By symmetry betwene D1 and D3 , the expectation of D3 is also
1 SE[A2 BC] + S(S − 1)E[AB]E[AC] .
3SP
D4 . By expanding D4 , we know that ED4 = S12 i,j E[A(i) B (i) A(j) C (j) ].
E[A(i) B (i) A(j) C (j) ] for each (i, j) is simple, and is omitted. The expectation of D4 is

The case-by-case analysis ofS−1
1
E[A2 BC] +E[AB]E[AC].
S
S
Simple algebra reveals that

P4

i=1 E[Di ] −S−1
S−1
S E[AB] S E[AC] is equal to Equation (6).

Proof D.7 (Proof of Lemma C.3). In this proof, we will only consider expectations under the full-data posterior. Hence, toalleviate notation, we shall write E instead of E1N : similarly, covariance and variance evaluations are understood to be at
w = 1N .
21

Small-data sensitivity of MCMC analyses

Applying Lemma D.1, the covariance of ψ̂n and ψ̂n i.e. the variance of ψ̂n is equal to(S − 1)2E{(g(β) − E[g(β)])2 (L(d(n) | β) − E[L(d(n) | β)])2 }S3S−1(S − 1)(S − 2)+
Var(L(d(n) | β))Var(g(β)) −Cov(g(β), L(d(n) | β))2 .
S3
S3Define the constant C to be the maximum over n ofCov(g(β), L(d(n) | β))2 + Var(g(β))Var(L(d(n) | β))+ E{(g(β) − E[g(β)])2 (L(d(n) | β) − E[L(d(n) | β)])2 }.

Simple algebra shows that Var(ψ̂n ) ≤ C
S.

Proof D.8 (Proof of Theorem C.1). Similar to the proof of Lemma C.3, expectations (and variances and covariances) are
understood to be taken under the full-data posterior.

Since ψ̂n is the biased sample variance, we know thatEψ̂n =S−1
ψn .
S

The bias of ψ̂n goes to zero at rate 1/S. Because of Lemma C.3, the variance also goes to zero at rate 1/S. Then, theppapplication of Chebyshev’s inquality shows that ψ̂n −→0→ ψn . Since N is a constant, the pointwise convergence |ψ̂n − ψn | −p→ 0.
implies the uniform convergence maxN
n=1 |ψ̂n − ψn | −pb − ∆(α)| −We now prove that |∆→ 0. We first recall some notation. The ranks r1 , r2 , . . . , rN sort the influences ψr1 ≤ ψr2 ≤P⌊N α⌋. . . ≤ ψrN , and ∆(α) = − m=1 ψrm I{ψrm < 0}. Similarly, v1 , v2 , . . . , vN sort the estimates ψ̂v1 ≤ ψ̂v2 ≤ . . . ≤ ψ̂vN ,b = − P⌊N α⌋ ψ̂v I{ψ̂v < 0}. It suffices to prove the convergence when ⌊N α⌋ ≥ 1: in the case ⌊N α⌋ = 0, both ∆
b
and ∆
m
m
m=1
and ∆(α) are equal to zero, hence the distance between them is identically zero. Denote the T unique values among ψnby u1 < u2 < . . . < uT . If T = 1 i.e. there is only one value, let ω := 1. Otherwise, let ω be the smallest gap between
subsequent values: ω := mint (ut+1 − ut ).
Suppose that maxN
n=1 |ψ̂n − ψn | ≤ ω/3: let A be the indicator for this event. For any n, each ψ̂n is in the interval[ψn − ω/3, ψn + ω/3]. In the case T = 1, clearly all k such that ψ̂k is in [ψn − ω/3, ψn + ω/] satisfy ψk = ψn . In thecase T > 1, since unique values of ψn are at least ω apart, all k such that ψ̂k is in [ψn − ω/3, ψn + ω/] satisfy ψk = ψn .
This means that the ranks v1 , v2 , . . . , vN , which sort the influence estimates, also sort the true influences in ascending order:
ψv1 ≤ ψv2 ≤ . . . ≤ ψvN . Since the ranks r1 , r2 , . . . , rN also sort the true influences, it must be true that ψvm = ψrm forall m. Therefore, we can writeb − ∆(α)| =|∆≤⌊N α⌋ Xm=1

ψvm I{ψvm < 0} − ψ̂vm I{ψ̂vm < 0}

⌊N α⌋

X

m=1

ψvm I{ψvm < 0} − ψ̂vm I{ψ̂vm < 0} .

We control the absolute values ψvm I{ψvm < 0} − ψ̂vm I{ψ̂vm < 0} .

For any index n, by triangle inequality,

ψn I{ψn < 0} − ψ̂n I{ψ̂n < 0} is at mostI{ψ̂n < 0}|ψn − ψ̂n | + |ψn ||I{ψ̂n < 0} − I{ψn < 0}|.
The first term is at most |ψn − ψ̂n |. The second term is at most I{|ψn − ψ̂n | ≥ |ψn |, ψn ̸= 0}. We next prove a boundon ψn I{ψn < 0} − ψ̂n I{ψ̂n < 0} that holds across n. Our analysis proceeds differently based on whether the set{n : ψn ̸= 0} is empty or not.
22

Small-data sensitivity of MCMC analyses• {n : ψn ̸= 0} is empty. This means ψn = 0 for all n. Hence, I{|ψn − ψ̂n | ≥ |ψn |, ψn ̸= 0} is identically zero.
• {n : ψn ̸= 0} is not empty. We then know that minn |ψn | > 0. Hence, I{|ψn − ψ̂n | ≥ |ψn |, ψn ̸= 0} is
upper bounded by I{|ψn − ψ̂n | ≥ minn |ψn |}. Since |ψn − ψ̂n | ≤ maxn |ψn − ψ̂n |, this last indicator is at mostI{maxn |ψn − ψ̂n | ≥ minn |ψn |}.
b − ∆(α)|. When {n : ψn ̸= 0} is empty, on A, |∆b − ∆(α)|
To summarize, we have proven the following upper bounds on |∆is upper bounded by⌊N α⌋ max |ψn − ψ̂n |(7)n=1

b − ∆(α)| is upper bounded by
When {n : ψn ̸= 0} is not empty, on A, |∆⌊N α⌋ max |ψn − ψ̂n | + ⌊N α⌋I{max |ψn − ψ̂n | ≥ min |ψn |}.
n

n=1

n

(8)

b − ∆(α)| > ϵ) converges to zero. For any positive ϵ, we know that
We are ready to show that Pr(|∆b − ∆(α)| > ϵ) ≤ Pr(|∆b − ∆(α)| > ϵ, A) + Pr(Ac ).
Pr(|∆p

The later probability goes to zero because maxN
→ 0.
n=1 |ψ̂n − ψn | −Suppose that {n : ψn ̸= 0} is empty. Using the upper bound Equation (7), we know that event in the former probabilitypN
implies that maxN
→ 0.
n=1 |ψ̂n − ψn | ≥ ϵ/⌊N α⌋: The probability of this event also goes to zero because maxn=1 |ψ̂n − ψn | −Suppose that {n : ψn ̸= 0} is not empty. Using the upper bound Equation (8), we know that event in the former probabilityN
implies that (maxN
n=1 |ψ̂n − ψn | + I{maxn |ψn − ψ̂n | ≥ minn |ψn |}) ≥ ϵ/⌊N α⌋. Since maxn=1 |ψ̂n − ψn | converges
to zero in probability, I{maxn |ψn − ψ̂n | ≥ minn |ψn |} also converges to zero in probability. Hence, the probability that(maxN
n=1 |ψ̂n − ψn | + I{maxn |ψn − ψ̂n | ≥ minn |ψn |}) ≥ ϵ/⌊N α⌋ converges to zero.
b − ∆(α)| > ϵ) goes to zero in both the case where {n : ψn ̸= 0} is empty and the complement case. As the
In all, Pr(|∆pb −choice of ϵ was arbitrary, we have shown ∆→ ∆(α).

Proof D.9 (Proof of Theorem C.2). Similar to the proof of Lemma D.1, we only consider expectations under the full-data
posterior. Hence, we will write E instead of E1N to simplify notation. Variance and covariance operations are also
understood to be taken un der the full-data posteiror. To lighten the dependence of the notation on the parameter β, we will
write g(β) as g and L(d(n) | β) as Ln when talking about the expectation of g(β) and L(d(n) | β).

Define the the following multivariate functionTf (β) := g(β), L(d(1) | β), g(β)L(d(1) | β), . . . , L(d(N ) | β), g(β)L(d(N ) | β) .

As defined, f (·) is a mapping from V -dimensional
space to 2N + 1-dimensional space. Since (β (1) , . . . , β (S) ) is an i.i.d.

(1)
(2)
(S)
sample, f (β ), f (β ), . . . , f (β ) is also an i.i.d. sample. Because of the moment conditions we have assumed, eachf (β) has finite variance. We apply the Lindeberg-Feller multivariate central limit theorem (van der Vaart, 1998, Proposition
2.27), and conclude that
!
√
1X
D
S
f (β (s) ) − Ef (β) −→ N (0, Ξ)S s
where the limit is S → ∞, and Ξ is a symmetric (2N + 1) × (2N + 1) dimensional matrix, which we specify next. It suffices
to write down the formula for (i, j) entry of Ξ where i ≤ j:
Var(g)
if i = j = 1





Cov(g,
L)
if i = 1, j > 1nCov(L , L )if i = 2n, j = 2mnm.
Ξi,j =

Cov(L
,
gL
)
if
i = 2n, j = 2m + 1
n
m





Cov(gLn , Lm )if i = 2n + 1, j = 2m
Cov(gLn , gLm ) if i = 2n + 1, j = 2m + 1
23

Small-data sensitivity of MCMC analyses

To relate the asymptotic distribution of f (β) to that of the vector ψ̂, we now use the delta method. Define the following
function which acts on 2N + 1 dimensional vectors and returns N dimensional vectors:Th([x1 , x2 , . . . , x2N +1 ]T ) := x3 − x1 x2 , x5 − x1 x4 , x7 − x1 x6 , . . . , x2N +1 − x1 x2N .
P
Written this way, clearly h(·) transform the sample mean S1 s f (β (s) ) into the estimated influences: ψ̂ =
P
h S1 s f (β (s) ) . Furthermore, h(·) applied to Ef (β) yields the vector of true influences: ψ = h (Ef (β)). h(·) is
continuously differentiable everywhere: its Jacobian is the N × (2N + 1) matrix
−x2 −x1 1
0
0 ... 0 −x4
0
0 −x1 1 . . . 0
Jh =  .
,
..
.
..
..
 ..
. 0 . . . 0
.
−x2N

0

0

...

0

...

1

which is non-zero. Therefore, we apply the delta method (van der Vaart, 1998, Theorem 3.1) and conclude that



√ 
D
S ψ̂ − ψ −→ N 0, Jh x=Ef (β) Ξ(Jh x=Ef (β) )T .

The (i, j) entry of the asymptotic covariance matrix is the dot product between the i-th row of Jh x=Ef (β) and the j-thcolumn of Ξ(Jh x=Ef (β) )T . The former is[−ELi , 0, 0, . . . , −Eg ,
|{z}

1
|{z}

, . . . , 0].

2i entry (2i+1) entry

The later is(−ELj )Cov(g, g) − (Eg)Cov(g, Lj ) + Cov(g, gLj )..
.
(−ELj )Cov(gLN , g) − (Eg)Cov(gLN , Lj ) + Cov(gLN , gLj )




.

Taking the dot product, we have that the (i, j) entry of the asymptotic covariance matrix is equal toCov(gLi , gLj ) − (Eg) [Cov(gLi , Lj ) + Cov(gLj , Li )]− [(ELj )Cov(g, gLi ) + (ELi )Cov(g, gLj )]+ (ELj )(ELi )Var(g)+ (Eg)2 Cov(Li , Lj )+ (Eg) [(ELj )Cov(g, Li ) + (ELi )Cov(g, Lj )]It is simple to check that the last display is equal to the covariance between (g−E[g])(Lj −E[Lj ]) and (g−E[g])(Li −E[Li ]).
Proof D.10 (Proof of Lemma C.4). We use the (shape, rate) parametrization of the gamma distribution. Let the prior over τ
be Gamma(α, β) where α, β > 0. Conditioned on observations, the posterior distribution of (µ, τ ) is normal-gamma:#!
"NN
N 1 X (n) 22τ ∼ Gamma α + , β +(x ) − x̄
,2
2 N n=1
ϵ ∼ N (0, 1),ϵ.
µ | τ, ϵ = x̄ + √Nτ

In this section, since we only take expectations under the original full-data posterior, we will lighten the notation’s
dependence on w, and write E instead of E1N . Similarly, covariance and variance operators are understood to be under the
full-data posterior.
√
For completeness, we compute Cov(µ, L(d(n) | µ, τ )). We know that µ − Eµ = ϵ/ N τ . The log likelihood, as a functionof τ and ϵ, isτ  11 2 x(n) − x̄ √1
− τ (x(n) − x̄)2 −
log
ϵ + √
ϵ τ.
2
2π
2
2N
N
24

Small-data sensitivity of MCMC analyses√√
The covariance of µ and L(d(n) | µ, τ ) is equal to the covariance between ϵ/ N τ and L(d(n) | µ, τ ). Since ϵ/ N τ iszero mean, the covariance is equal to the expectation of the product. Since ϵ is indedependent of τ , many of the terms thatform the expectatin of the product is zero. After some algebra, the only term that remains is (n)x − x̄ 2
x(n) − x̄
E
ϵ =
.
N
N
2
2
To compute the asymptotic variance of ψ̂n , it suffices to compute the expectation of Nϵ τ L(d(n) | µ, τ ) − EL(d(n) | µ, τ ) .
The calculations are simple, but tedious, and we omit them. We will only state the result. The expectation of
2
ϵ2
(n)
| µ, τ ) − EL(d(n) | µ, τ ) isN τ L(d1E[τ −1 (τ − Eτ )2 ] (x(n) − x̄)44N

3 + E[τ −1 (τ − Eτ )] E[τ −1 (log τ − E log τ )](x(n) − x̄)2
−
+
N2
2N
1
1
1
+
E[τ −1 ] +E[τ −1 (log τ − E log τ )2 ] − 2 E[τ −1 (log τ − E log τ )2 ].
2N 3
2N
N
Since the asymptotic variance is equal to this expectation minus the square of the covariance between L(d(n) | µ, τ ) and µ,
our final expression for the asymptotic variance Σn,n is


1
−1
2
E[τ (τ − Eτ ) ] (x(n) − x̄)44N

2 + E[τ −1 (τ − Eτ )] E[τ −1 (log τ − E log τ )](x(n) − x̄)2
−
+
N2
2N
1
1
1
+
E[τ −1 ] +E[τ −1 (log τ − E log τ )2 ] − 2 E[τ −1 (log τ − E log τ )2 ].
2N 3
2N
N
The constants D1 , D2 , and D3 mentioned in the lemma statement can be read off this last display. It is possible to replace
the posterior functionals of τ with quantities that only depends on the prior (α, β) and the observed data. Such formulas
might be helpful in studying the behavior of Σn,n in the limit where some x(n) becomes very large.

E. Experimental Setup
In this section, we only describe the checks: for the actual results, see Section 4 and Appendices F.1 and F.2.
A practitioner with a particular definition of “small data” can set α to reflect their concern. We consider a number of α
values. We set the maximum value of α to be 0.01. This choice is motivated by Broderick et al. (2023). Many analyses
are non-robust to removing 1% of the data, and we a priori think that α > 1% is a large amount of data to remove. Wevary log10 (α) in an equidistant grid of length 10 from −3 to −2. The ten values are 0.10%, 0.13%, 0.17%, 0.22%, 0.28%,
0.36%, 0.46%, 0.60%, 0.77% and 1.00%.
For the range of dropout fraction specified above and across three common quantities of interest corresponding to sign,
significance, and significant result of opposite sign changes, we walk through what a practitioner would do in practice(although they would choose only one α and one decision). Our method proposes an influential data subset and a change inthe quantity of interest, represented by a confidence interval.
Ideally, we want to check if our interval includes the result of the worst-case data to leave out. We are unable to do so, since
we do not know how to compute the worst-case result in a reasonable amount of time. We settle for the following checks.
In the first check, for a particular MCMC run, we plot how the change from re-running minus the proposed data compares tothe confidence interval. We recommend the user run this check if re-running MCMC a second time is not too computationallyexpensive.
Unfortunately, such refitting does not paint a complete picture of approximation quality. For instance, the MCMC run might
be unlucky since MCMC is random. To be more comprehensive, we run additional checks. We do not expect users to run
25Small-data sensitivity of MCMC analyses

these tests, as their computational costs are high. The central question is how frequently (under MCMC randomness) the
confidence interval includes the result after removing the worst-case data. Since we estimate the worst-case change witha linear approximation, a natural way to answer this question is with two separate checks: while Appendix E.1 checks
how frequently the confidence interval includes the result of the linear approximation i.e. the AMIP, Appendix E.3, checks
whether the linear approximation is good. To understand why we observe the coverage in Appendix E.1, in Appendix E.2
we isolate the impact of the sorting step in the construction of our confidence interval.
E.1. Estimate coverage of confidence interval for AMIP
We estimate how frequently [∆lb (α), ∆ub (α)] covers the AMIP by using another level of Monte Carlo. Recall that[∆lb (α), ∆ub (α)] is intended to be a confidence interval covering ∆(α) a fraction η of the time. If the estimated coverage isfar from η, we have evidence that [∆lb (α), ∆ub (α)] does not achieve the desired nominal coverage.
We draw J Markov chains: we set J = 960. On each chain, we estimate the influences and construct the confidence interval[∆lb (α), ∆ub (α)]. From all the chains, for each n, we have J estimates of ψn . We take the sample mean across chains, anddenote this by ψn∗ : because of variance reduction through averaging, ψn∗ is a much better estimate of ψn than individual ψ̂n .
We denote the indices of the ⌊N α⌋ most negative ψn∗ by U ∗ (α). We sort ψn∗ across n and sum the ⌊N α⌋ most negativeψn∗ . This sum is denoted by ∆∗ (α): we use it in place of the ground truth ∆(α). We use the sample mean of the indicators
I{∆∗ (α) ∈ [∆lb (α), ∆ub (α)]} as the point estimate of the coverage. We also report a 95% confidence interval for thecoverage. This interval is computed using binomial tests designed in Clopper and Pearson (1934) and implemented as R’sbinom.test() function.
E.2. Estimate coverage of confidence intervals for sum-of-influence
It is possible that the estimated coverage of [∆lb (α), ∆ub (α)] is far from the nominal η. We suspect that such a discrepancy
comes from the sorting of ψ̂n to construct ∆(α). To modularize out the sorting, we consider a target of inference
P that is
simpler than ∆(α). At a high level, we fix an index set I, and define the target to be the sum of influences in I: n∈I ψn .
P
b does.
On each sample (β (1) , . . . , β (S) ), our point estimate is n∈I ψ̂n : this estimate does not involve any sorting, while ∆P
We construct the confidence interval, [V lb , V ub ], from the block bootstrap distribution of n∈I ψ̂n . The difference between
b is that the former is not[V lb , V ub ] and [∆lb (α), ∆ub (α)], which is constructed from the block bootstrap distribution of ∆,lb
ub
based on sorting the influence estimates. If the actual coverage of [V , V ] is close to the nominal value, we have evidence
that the miscoverage of [∆lb (α), ∆ub (α)] is due to this sorting.
From Appendix E.1 we use ψn∗ and the associated ∆∗ (α) and U ∗ (α) as replacement for ground truths. We set I to
be U ∗ (α). We run another set of J Markov chains: for each chain, we construct the confidence interval [V lb , V ub ] byP
sampling from the block bootstrap distribution of the estimator n∈I ψ̂n . We report the sample mean of the indicatorsP
I{ n∈I ψn∗ ∈ [V lb , V ub ]} as our point estimate of the coverage. We also report a 95% confidence interval for the coverage.
This interval is computed using binomial tests designed in Clopper and Pearson (1934) and implemented as R’s binom.test()function.
E.3. Re-running MCMC on interpolation path
Ideally, we want to know the difference between the Maximum Influence Perturbation and the AMIP. As we have established,
we do not know how to compute the former
efficiently. We settle for checking the linearity approximation made in Section 3.1
P
i.e. estimating ϕ(w) − ϕ(1N ) with n (wn − 1)ψn . In particular, we expect the first-order Taylor series approximation to
be arbitrarily good for w arbitrarily close to 1N . By necessity, we are interested in some w∗ that has a non-trivial distancefrom 1N . Plotting the quantity of interest ϕ(w) on an interpolation path between 1N and w∗ , we get a sense of how much
we have diverged from linearity by that point.
From Appendix E.1, we have ψn∗ as our replacement for the ground truth ψn . We focus on α = 0.05: 5% is a large
amount of data to remove, and a priori we expect the linear approximation to be poor. Recall that U ∗ (0.05) is the set of⌊0.05N ⌋ observations that are most influential according to sorted ψn∗ . Let w∗ be the N -dimensional weight vector thatis 1 for observations in U ∗ (0.05) and 0 otherwise. For ζ ∈ [0, 1], the linear approximation of ϕ(ζw∗ + (1 − ζ)1N ) isϕ(1N ) + ζ∆∗ (0.05). In the extreme ζ = 0, we do not leave out any data. In the extreme ζ = 1, we leave out the entirety of26

Small-data sensitivity of MCMC analyses

U ∗ (0.05) i.e. 5% of the data. An intermediate value ζ roughly corresponds5 to removing (ζ5)% of the data. We discretize[0, 1] with 15 values: 0, 0.0010, 0.0016, 0.0027, 0.0044, 0.0072, 0.0118, 0.0193, 0.0316, 0.0518, 0.0848, 0.1389, 0.2276,
0.3728, 0.6105, 1. For each value on this grid, we run MCMC to estimate ϕ(ζw∗ + (1 − ζ)1N ), and compare it to the linear
approximation.

F. Other Quality Checks
In our experiments, we find that our approximation works well for a simple linear model. But we find that it can struggle inhierarchical models with more complex structure. While Section 4 has plotted our confidence interval and the result after
removing proposed data, we have not yet discussed the checks in Appendices E.1 to E.3. Furthermore, Section 4 has notdescribed the context to each analysis. This section a) fills in the missing context and b) shows the outcome of the additional
quality checks.
F.1. Linear model
We consider a slight variation of a microcredit analysis from Meager (2019). In Meager (2019), conclusions regarding
microcredit efficacy were based on ordinary least squares (OLS). We refer the reader to Broderick et al. (2023, Section 4.3.2)
for investigations of such conclusions’ non-robustness. Here, we instead consider an analogous Bayesian analysis using
MCMC, and we examine the robustness of conclusions from this analysis. Even for this very simple Bayesian analysis, it ispossible to change substantive conclusions by removing a small fraction of the data.
Our quality checks suggest that our approximation is accurate. Our confidence interval contains the refit after removing
the proposed data. The actual coverage of the confidence interval for AMIP is close to the nominal coverage. The actual
coverage of the confidence interval for sum-of-influence is also close to the nominal coverage. Even for dropping 5% of thedata, the linear approximation is still adequate.
F.1.1. BACKGROUND .
Meager (2019) studies the microcredit data from Angelucci et al. (2015), which was an RCT conducted in Mexico. There
are N = 16,560 households in the RCT. Each observation is d(n) = (x(n) , y (n) ), where x(n) is the treatment status andy (n) is the profit measured. The log-likelihood for the n-th observation is L(d(n) | µ, θ, σ) = − 2σ1 2 (y (n) − θx(n) − µ)2 −1
2
2 log(2πσ ). Here, the model parameters are baseline profit µ, treatment effect θ, and noise scale σ. The most interesting
parameter is θ: as x(n) is binary, θ compares the means in the treatment and control groups. Meager (2019) estimates themodel parameters with OLS.
Our variation of the above analysis is as follows. We put t location-scale distribution priors on the model parameters, with theadditional constraint that the noise scale σ is positive. Recall that the t location-scale distribution has three hyperparameters:
ν, µ, σ. ν is the degrees of freedom, µ is the location, and σ is the scale. The density at y of this distribution is1Γ((ν + 1)/2)√
Γ(ν/2)πνσ 2(y − µ)21+
νσ 2

−(ν+1)/2.

We set the the prior over µ to be t location-scale with degrees of freedom 3, location 0, and scale 1000. We set the theprior over θ to be t location-scale with degrees of freedom 3, location 0, and scale 1000. We set the the prior over σ to be t
location-scale with degrees of freedom 3, location 0, and scale 1000.
We use Hamiltonian Monte Carlo (HMC) as implemented in Stan (Carpenter et al., 2017) to approximate the full-data
posterior. We draw S = 4000 samples.
F.1.2. RUNTIME .
The running of our approximation takes very little time compared to the running of the original analysis. Generating thedraws in Figure 1 took 3 minutes on MIT Supercloud (Reuther et al., 2018). For one α and one quantity of interest, it took
less than 5 seconds to make a confidence interval for what happens if we remove the most extreme data subset. A user might
check approximation quality by dropping a proposed subset and re-running MCMC: each such check took us around 3
5This correspondence is not exact, since for ζ < 1, all observations in U ∗ (0.05) are included in the analysis, only with downplayed
contributions.

27

Small-data sensitivity of MCMC analyses

Figure 5. (Linear model) Monte Carlo estimate of AMIP confidence interval’s coverage. Each panel corresponds to a target conclusionchange. The dashed line is the nominal level η = 0.95. The solid line is the sample mean of the indicator variable for the event thatground truth is contained in the confidence interval. The error bars are confidence intervals for the population mean of these indicators.

Figure 6. (Linear model) Monte Carlo estimate of sum-of-influence confidence interval’s coverage. Each panel corresponds to a target
conclusion change. The dashed line is the nominal level η = 0.95. The solid line is the sample mean of the indicator variable for the
event that ground truth is contained in the confidence interval, and error bars are confidence intervals for the population mean of these
indicators.

minutes, the runtime of the original analysis.
F.1.3. A DDITIONAL QUALITY CHECKS .
Figure 5 shows that the actual coverage of the confidence interval for the AMIP is close to the nominal one, across α. As thehalf-width of each error bar is small (only 0.02), we believe that the difference between the true coverage and our point
estimate of it is small. For either ‘sign’ or ‘both’ QoI, the error bars do not contain the nominal η. However, the difference
between the point estimate and the nominal η is only 0.03 at worst, which is small. For the ‘sig’ QoI, the point estimate iswithin 0.005 of the nominal value, and the error bars contain the nominal η.
Figure 6 shows that the actual coverage of the confidence interval for the sum-of-influence is close to the nominal one acrossα. The absolute errors between our estimate of coverage and the nominal η are similar to those seen in Figure 5. This
success suggests that the default block length, L = 10, is appropriate for this problem.
Figure 7 shows that the linear approximation works very well. It is somewhat remarkable that the linear approximation is
this good even after dropping 5%, which we consider to be a large fraction of data. The horizontal axis (‘scale’) is the sameas ζ in Appendix E.3. For all quantities of interest, the linear approximation and the refit lie mostly on top of each other:
towards the right end of each panel, the approximation slightly underestimates the refit.
28

Small-data sensitivity of MCMC analyses

Figure 7. (Linear model) Quality of the linear approximation. Each panel corresponds to a target conclusion change. The solid blue line is
the full-data fit. The horizontal axis is the distance from the weight that represents the full data. We plot both the refit from rerunning
MCMC and the linear approximation of the refit.

F.2. Hierarchical model on tree mortality data
In the final experiment, we break from microcredit and look at ecological data. In particular, we consider a slight tweak of
the analysis of European tree mortality from Senf et al. (2020). Senf et al. are acutely aware of generalization concerns.
While previous work on tree death had been limited in both time and space, Senf et al. (2020) designs a large study thatstretches across Europe and over 30 years, in hopes of making a broad-scale assessment. Our work shows that, even after an
expansive study with generalization in mind, one might still worry about applying the findings at large, because of small-data
sensitivity.
Our approximation struggles in this case. For the particular MCMC run used to estimate the full-data posterior, ourconfidence interval does not contain the refit after removing the proposed data. As each MCMC run is already highly
time-consuming, we do not run quality checks on the whole dataset. We settle for running quality checks on a subsample ofthe data. On the subsampled data, the confidence interval for AMIP undercovers: the undercoverage is severe for one of thequantities of interest. However, the confidence interval for sum-of-influence is close to achieving the nominal coverage. For
all three quantities of interest, the linear approximation is good up to removing roughly 1.1% of the data. For two of thethree, it breaks down afterwards: for the remaining one, it continues to be good up to 3%, then falters.
As articulated in Appendix E, we think that dropping more than 1% of the data is already removing a large fraction. We are
not worried about the Maximum Influence Perturbation for such α. So, that the linear approximation stops working after
1.1% is not a cause for concern.
F.2.1. BACKGROUND .
Senf et al. (2020) studies the relationship between drought and tree death in Europe. To identify the association, they havecompiled a dataset with N = 87,390 observations. Europe is divided into 2,913 regions, and the data spans 30 years: eachobservation is a set of measurements made in a particular region, which we denote as l(n) , and at a particular year, which wedenote as t(n) . For our purposes, it suffices to know that the measurement of (the opposite of) drought is called climaticwater balance, and we denote it as x(n) : larger values of x(n) indicate that more water is available i.e. there is less drought.
The response of interest, y (n) , is excess death of tree canopy.
In our experiment, we mostly replicate (Senf et al., 2020)’s probabilistic model: we use the same likelihood, and makeonly an immaterial modification in the choice of priors. The likelihood for the n-th observation is exponentially modified
Gaussian with standard deviation σ, scale λ and mean (time)(region)(time)(location)µt(n) + µl(n) + µ + θt(n) + θl(n)+ θ x(n) + f (x(n) ),P10with f (x) := i=1 Bi (x)γi where Bi ’s are fixed thin plate spline basis functions (Wood, 2003) and γi ’s are random:2γi ∼ Normal(0, σ(smooth)). In all, the parameters of interest are• Fixed effects: µ and θ.
29

Small-data sensitivity of MCMC analyses

Figure 8. (Hierarchical model on subsampled tree mortality) Histogram of effect MCMC draws. See Figure 1 for the meaning of the
distinguished lines.
(time)(time)(region)• Random effects: time (µt(n) , θt(n) ) and location (µl(n)(location), θl(n)).

• Degree of smoothing: σ(smooth) .
Since there are many regions (nearly 3,000) and periods of time (30), the number of random effects is large. Senf et al.
(2020) uses brms()’s default priors for all parameters: in this default, the fixed effects are given improper uniform priors over
the real line. To work with proper distributions, we set the priors for the random effects and degree of smoothing in the same
way set by Senf et al. (2020). For fixed effects, we use t location-scale distributions with degrees of freedom 3, location 0,
and scale 1000.
At a high level, both Senf et al. (2020)’s prior and our prior share strength across regions and times by modeling the random
effects as coming from a some common global distributions. However, while Senf et al. (2020) uses an improper prior, we
use a proper one. Numerically, there is no perceptible difference between the two. Theoretically, we prefer working withproper priors to avoid the integrability issue mentioned around Assumption 2.1.
Following Senf et al. (2020), we make conclusions based on posterior functionals of θ. Roughly speaking, θ is the average(across time and space) association effect that water balance has on excess tree death. We use S = 8000 HMC draws toapproximate the posterior.
F.2.2. RUNTIME
The running of our approximation takes very little time compared to the running of the original analysis. Generating thedraws in Figure 3 took 12 hours. For one α and one quantity of interest, it took less than 2 minutes to make a confidence
interval for what happens if we remove the most extreme data subset. A user might check approximation quality by dropping
a proposed subset and re-running MCMC: each such check took us around 12 hours, which is the runtime of the original
analysis.
F.2.3. R ESULTS ON SUBSAMPLED DATA .
Running MCMC on the original dataset of size over 80,000 took 12 hours. In theory, we can spend time (on the order of
thousands of hours) to run our quality checks, but we do not do so. Instead, we subsample 2,000 observations at random
from the original dataset. Each MCMC on this subsample takes only 15 minutes, making it possible to run quality checks ina few hours instead of weeks. We hope that the subsampled data is representative enough of the original data that the quality
checks on the subsampled data are indicative of the quality checks on the original data.
We use the same probabilistic model to analyze the subsampled data. Figure 8 plots the histogram of the association
effect draws and sample summaries. Based on the draws, a forest ecologist might tentatively say that drought is positively
associated with canopy mortality if they relied on the posterior mean, but refrain from conclusively deciding, since theuncertainty interval contains zero.
Figure 9 shows our confidence intervals and the actual refits. Similar to Figure 4, our confidence intervals predict a more
extreme change than realized by the refit. The overestimation is most severe for ‘both’ QoI.
In Figure 10, the confidence interval for AMIP undercovers for all quantities of interest. The actual coverage decreases as30Small-data sensitivity of MCMC analyses

Figure 9. (Hierarchical model on subsampled tree mortality) Confidence interval and refit. See the caption of Figure 2 for the meaning of
the panels and the distinguished lines.

Figure 10. (Hierarchical model on subsampled tree mortality) Monte Carlo estimate of coverage of confidence interval for ∆(α). See
Figure 5 for the meaning of the panels and the distinguished lines.

α increases. The undercoverage is most severe for ‘sig’ QoI: while the nominal level is 0.95, the confidence interval forthe true coverage only contains values less than 0.15. This translates to a relative error of over 84%. In other words, ourconfidence interval for significance change is too narrow, and rarely contains the AMIP. For ‘both’ QoI and ‘sig’ QoI, theworst-case relative error between the nominal and the estimated coverage, which occurs under the largest α, is 15.7%.
In Figure 11, the estimated coverage of the confidence interval for sum-of-influence is close to the nominal coverage. Notethe stark contrast in the vertical scale of the ‘sig’ panel in Figure 10 with that in Figure 11. At worst, our point estimate of
the true coverage is 0.04 less than the nominal level, which is only a 4.2% relative error. This success of the block bootstrapbindicate that the undercoverage observed in Figure 10 can be attributed to the sorting step involved in the definition of ∆.
We leave to future work to investigate why the interference cause by the sorting step is so much more severe for changing
the significance than for changing sign or generating significant result of the opposite sign.
Figure 12 shows that the linear approximation is good for the posterior mean (‘sign’ QoI) and the left credible endpoint(‘both’ QoI) up to ζ = 0.2276: in data percentages, this is roughly 1.1%. For larger ζ, the refit for ‘both’ QoI plateaus while
the linear approximation continues to increase, and the linear approximation for posterior mean slightly underestimates
it. For the left endpoint (‘both’ QoI), the linear approximation is close to the refit up to ζ = 0.6105 (roughly 3% of data);
afterwards, the left endpoint increases while the linear approximation continues to decrease.

31Small-data sensitivity of MCMC analyses

Figure 11. (Hierarchical model on subsampled tree mortality) Monte Carlo estimate of coverage of confidence interval for sum-ofinfluence. See Figure 6 for the meaning of the panels and the distinguished lines.

Figure 12. (Hierarchical model on subsampled tree mortality) Quality of linear approximation. See Figure 7 for the meaning of the panels
and the distinguished lines.

32