Robustness to Missing Data: Breakdown Point Analysis

Robustness to Missing Data: Breakdown Point Analysis

Keywords: Missing data, generalized method of moments, robustness, sensitivity analysis.

JEL classification: C01, C14, C18, C21, C25.

\pagenumbering{gobble}

\pagenumbering{arabic}

Introduction

Virtually every economic dataset is plagued by missing and incomplete records. Survey nonresponse is the most visible cause, and appears to be worsening over time. bollinger2019trouble report that the Current Population Survey's Annual Social and Economics Supplement item and whole nonresponse has been increasing, reaching 43 percent in 2015. By linking these data with the Social Security Administration Detailed Earnings Record, the authors show that the distribution of nonresponders differs from that of responders even after conditioning on a large set of covariates.

Samples with missing or incomplete observations fail to identify the population distribution manski2005partial. To make progress, researchers commonly apply standard procedures to the complete observations. This practice is typically justified by assuming the data are “missing completely at random” (MCAR): the assumption that incomplete observations follow the same distribution as that of the complete observations. In many settings such an assumption is implausible. Without it, the conclusions drawn are uncomfortably qualified as being about the distribution of the complete observations, rather than the actual distribution of interest.

This paper proposes a method to investigate the robustness of a conclusion regarding the whole population. Results are more robust when overturning them would require more selection. To make this intuition precise, selection is measured as the divergence from the distribution of complete observations to the distribution of incomplete observations. Many statistical divergences can be used to measure selection. Squared Hellinger is an attractive choice for this purpose, as it can be interpreted as a measure of how well the variables under study would predict an observation being complete. This gives the values of the selection measure context, allowing researchers to gauge how much selection can be expected in a given setting. The breakdown point is the minimum amount of selection needed to overturn a conclusion. Readers who doubt the setting exhibits that much selection will find the conclusion compelling.

The estimator of the breakdown point proposed below can be applied to conclusions drawn from a model identified with the generalized method of moments (GMM) hansen1982large. This includes most models used in applied econometrics, including linear regression, instrumental variable models, binary choice models such as the logit model, and many more. In a model identified with GMM, the breakdown point is the constrained minimum of the value function of a convex optimization problem. Estimators of the breakdown point are constructed from the dual of this convex inner problem, and shown to be $\sqrt{n}$-consistent and asymptotically normal. Lower confidence intervals are simple to construct. Reporting the point estimates and lower confidence intervals of the breakdown point is a simple, concise way to communicate a result's robustness.

As a demonstration of breakdown point analysis, the paper concludes with an investigation of the robustness of results from three randomized controlled trails barham2024experimental, bandiera2020women, giacobino2024schoolgirls. Many randomized controlled trials conducted in developing countries suffer from missing data that results from attrition, often due to study subject migration. Results can change meaningfully when migrants are carefully tracked and included in the sample, suggesting that the data are not missing at random molina2025attrition. The three randomized controlled trials studied here are evaluated through intent-to-treat estimates implemented through linear regression. The breakdown point analyses thus show a range of breakdown point estimates resulting from variation in real world data, rather than variation in methodology. Some claims are notably more robust than others, even when a similar amount of data is missing and the estimates found by dropping incomplete observations are similar.

The breakdown point analysis proposed here has a number of advantages over existing methods for incomplete datasets. Sample selection models consider regressions with samples where the dependent variable is sometimes missing, and obtain point identification by modeling the selection process heckman1979sample, das2003nonparametric. These models often require the data include a variable changing the probability of observation but not the dependent variable. This “exclusion restriction” is difficult to satisfy in many applications. Sample selection models can be identified without such an exclusion restriction provided the researcher makes additional functional form restrictions, as in escanciano2016identification. The breakdown point approach proposed here can be used on most common models identified with GMM, including but not restricted to regressions with missing outcomes. It requires neither additional data nor additional modeling assumptions. The breakdown point can be estimated even if the incomplete observations are in fact completely missing, a distinct possibility when using survey data.

The econometric literature on missing data has also explored bounding the parameter of interest based on the support manski2005partial, horowitz2006identification. If all parameter values within these “worst-case” bounds satisfy the researcher's conclusion, the conclusion is undoubtedly robust. Unfortunately, the bounds may be uninformative in practice. Proponents of this approach are well aware these bounds are conservative, and propose this exercise as a place to begin an investigation rather than end one. Additional identifying assumptions should then be considered, in order to make plain to readers what needs to be assumed to reach a given conclusion manski2013response. The breakdown approach proposed here is a simple version of this exercise, as the assumption that selection is less than the breakdown point leads one to conclude the hypothesis under investigation.

A growing literature advocates for breakdown analysis as a general, tractable method to assess the sensitivity of a result to relaxations of identifying assumptions. The term “identification breakdown point” can be found as early as horowitz1995identification in the context of corrupted data. One well-known example of breakdown analysis is altonji2005selection, which considers linear regressions suffering from omitted variable bias and proposes measuring how strong selection on unobservables would need to be (relative to measured selection on observables) to attribute the entire estimated effect to selection. This idea was developed further in oster2019unobservable, is widely used in empirical economics, and is an area of active research; see, e.g., masten2025effect, diegert2025assessing, and the references therein. Breakdown analyses are also commonly used to evaluate the sensitivity of results to identifying assumptions made for causal inference, with recent examples including masten2020inference, bonvini2022sensitivity, rosenbaum2023sensitivity, and spini2024robustness.

This paper is not the first to notice the appeal of breakdown point analysis in the context of missing data. kline2013sensitivity considers a setting with a missing scalar, propose measuring selection with the maximal Kolmogorov-Smirnov (KS) distance between the conditional distributions of complete and incomplete observations across all values of covariates, and advocates for “reporting the minimal level of selection necessary to undermine a hypothesis,” (p. 233). In their setting, one minus the KS distance can be interpreted as the proportion of the missing population assumed to be missing at random. The methodology proposed here has some notable advantages. First, measuring selection with the maximal KS distance limits researchers to the case where only a scalar is missing, while measuring selection as proposed here allows any number of variables to be missing. Second, in a given setting it is easier to gauge whether the variables under study are likely to be good predictors of missingness than what share of the missing data is missing at random. This makes squared Hellinger a more natural measure of selection than KS distance. Which approach is more tractable will depend on the parameter of interest. kline2013sensitivity derives sharp, closed form bounds to the conditional quantiles of the missing variable, and then frames the conclusion to be investigated in terms of those quantiles. This paper assumes the parameter of interest is identified with GMM and uses the model directly, but gives up closed form solutions. In theory this could lead to computational difficulties, but the simulations in section (ref) and application in section (ref) present no issue.

Estimation of the breakdown point remains tractable due to the use of $f$-divergences to measure selection. These divergences, defined and discussed in section (ref), have seen widespread use in the sensitivity analysis literature. christensen2023counterfactual estimate the identified set of counterfactual predictions from structural models when the distribution of latent variables is allowed to vary within an $f$-divergence neighborhood of a given parametric specification. jin2022sensitivity use $f$-divergences to characterize deviations from unconfoundedness in causal inference. The most famous $f$-divergence may be Kullback-Leibler, which is used to study local model misspecification in bonhomme2022minimizing, sensitivity to the choice of a Bayesian prior in ho2023global, and breakdown of causal inference results in spini2024robustness.

The remainder of this paper is structured as follows. Section (ref) formalizes the setting, the proposed measure of selection, and the breakdown point. The dual problem is presented and discussed in section (ref). Section (ref) defines the estimator and states the main results on estimation and inference, which are proven in the Supplementary Material. Section (ref) presents a simulation study investigating the finite sample performance of the estimators. The methodology is illustrated in section (ref), which contains estimates of the breakdown point of conclusions drawn from several randomized controlled trails suffering from attrition. Section (ref) concludes.

Measuring selection and breakdown analysis

Suppose the available data is the i.i.d. sample $\{(D_i, D_i Y_i, X_i)\}_{i=1}^n$, where $Z_i \equiv (Y_i, X_i)\in \mathbb{R}^{d_y} \times \mathbb{R}^{d_x}$ contains the variables of interest and $D_i \in \{0,1\}$ indicates whether $Y_i$ is observed. Note that $Y_i$ may be a vector, and $X_i$ may be empty. Variables are organized into the vectors $Y_i$ and $X_i$ based on whether they are occasionally missing; there is no need for $Y_i$ to be an outcome in the analysis. Let $p_D \equiv P(D = 1)$ denote the probability of observing $Y$, $P_1$ the distribution of $Z$ conditional on $D = 1$, and $P_0$ the distribution of $Z$ conditional on $D = 0$. $P_1$ and $P_0$ are called the complete case and incomplete case distributions respectively. The distribution of interest is the unconditional distribution of $Z$, given by $p_D P_1 + (1-p_D) P_0$. When $X$ is nonempty, the marginal distribution of $X$ conditional on $D = d$ is denoted $P_{dX}$.

Two assumptions made below are worth highlighting when introducing the setting. First, it is assumed that $P_0$ is absolutely continuous with respect to $P_1$, meaning that for every set $A$ with $P((Y,X) \in A \mid D = 1) = 0$ one has $P((Y, X) \in A \mid D = 0) = 0$ as well. This assumption, denoted $P_0 \ll P_1$, facilitates measuring selection with a statistical divergence. It is natural in some settings and may be restrictive in others; see remark (ref) below for additional discussion. Second, $X$ is assumed to have the same finite support when $D = 0$ as when $D = 1$. This assumption can be relaxed, as discussed in remark (ref).

To fix ideas, consider data collection via survey. $Y$ is a vector of data the survey hopes to collect, which is observed only if the recipient responds ($D = 1$). The survey's response rate, $p_D = P(D = 1)$, is essentially always less than one in practice. It is common for administrative data to provide basic information about a survey recipient (such as age, occupation, etc.), which is collected in $X$.

Analyses based on the complete observations may not convince researchers who worry that $P_0$ differs from $P_1$. Such concerns are common, as few settings plausibly satisfy the missing completely at random assumption. However, it is often similarly implausible that $P_0$ differs greatly from $P_1$. Researchers who convincingly argue that $P_0$ is not too different from $P_1$ can still convince their audience of conclusions drawn from an analysis of $P_1$.\footnote{In some cases, such as correctly specified regression models, it suffices that the conditional distributions $f_{Y \mid X = x, D = 0}(y \mid x)$ are the same as the identified $f_{Y \mid X = x, D = 1}(y \mid x)$. This weaker “missing at random” (MAR) assumption is also rarely plausible in practice, and analyses based on this assumption often rely heavily on the model being correctly specified.}

A quantitative measure of the difference between $P_1$ and $P_0$ is needed to make this argument formal and convincing. The statistics literature provides a natural solution in the form of divergences: functions mapping two probability distributions to the nonnegative real line that take value zero if and only if the two distributions are the same. There are many such functions. To be useful as a measure of selection, a divergence should have a tractable interpretation, so that researchers can gauge whether a given amount of selection is reasonable for their setting.

An interpretable measure of selection

Missing data cause greater concern when researchers expect the variables of interest ($Z$) to be a good predictor of incompleteness ($D$). Consider again the example of data collection via survey. Researchers are rightfully more concerned about survey nonresponse when asking about the respondent's arrest record than when asking for opinions on recent television programming. People with criminal records may be less willing to answer questions about that record.\footnote{For example, brame2012cumulative estimate the cumulative prevalence of arrest from ages 8 to 23 from a survey directly asking about prior arrests. The authors report upper and lower bounds derived by assuming the entire set of nonresponders had or had not been arrested, essentially the worst-case bounds advocated for by manski2005partial.} This suggests that the distribution of responders may look quite different from the distribution of nonresponders, and that criminal records would be a good predictor of nonresponse.

To illustrate this more formally, let $f_1$ and $f_0$ be densities of $P_1$ and $P_0$ with respect to $p_D P_1 + (1-p_D)P_0$ respectively:

An optimist may assume $D$ is independent of $Z$, implying that $P(D = 1 \mid Z = z) = P(D = 1) = p_D$ and $f_0 = f_1 = 1$. This would imply that $P_1$ and $P_0$ are the same distribution; i.e. the data are missing completely at random. In contrast, a pessimist may assume $D$ is close to a deterministic function of $Z$, allowing $Z$ to predict $D$ well. This would imply $P(D = 1 \mid Z = z)$ is close to $1$ or $0$ for many values of $z$, and that $f_1$ differs greatly from $f_0$.

As in the survey example, the setting often makes it clear whether $Z$ would be a good predictor of $D$. This heuristic is useful to identify and discuss selection concerns. The following lemma shows that measuring selection as the squared Hellinger distance between $P_0$ and $P_1$ captures this intuition, with larger values corresponding to $Z$ having greater capability of predicting $D$.\footnote{The Hellinger distance between probability measures $Q$ and $P$ is

where $\lambda$ is any measure dominating both $P$ and $Q$.}

All results are proven in the Supplementary Material. Equation (ref) states that the squared Hellinger distance between $P_0$ and $P_1$ is the expected percent of the standard deviation of $D$ reduced by conditioning on $Z$. In the extreme case where $\text{Var}(D \mid Z) = \text{Var}(D)$, equation (ref) implies $H^2(P_0, P_1) = 0$ and the conditional distributions are the same. As the ability of $Z$ to predict $D$ grows, the variance of $D$ conditional on $Z$ decreases and $H^2(P_0, P_1)$ grows toward one.

Divergences

Squared Hellinger provides an intuitive measure of selection, but there are many other options. A function $d (\cdot \Vert \cdot)$ mapping two probability distributions $P$ and $Q$ to $\mathbb{R}$ is called a divergence if $d(Q \Vert P) \geq 0$, with equality if and only if $P = Q$. Divergences need not be symmetric nor satisfy the triangle inequality. The set of $f$-divergences are particularly well behaved. Given a convex function $f : \mathbb{R} \rightarrow [0,\infty]$ satisfying $f(t) = \infty$ for $t < 0$ and taking a unique minimum of $f(1) = 0$, the corresponding $f$-divergence is given by

where $Q \ll P$ denotes absolute continuity of $Q$ with respect to $P$. Many popular divergences are equal to $f$-divergences when $P$ dominates $Q$.\footnote{By convention, $f$ is specified on its effective domain, the set $\text{dom}(f) \equiv \{t \in \mathbb{R} \; : \; f(t) < \infty\}$. The value of $f$ at $t = 0$ is set equal to $\lim_{t\rightarrow 0^+} f(t)$. Outside of $dom(f)$, $f$ takes the value $+\infty$.}

Although squared Hellinger has intuitive appeal outlined in Section (ref), the breakdown point analysis proposed in this paper remains tractable for any $f$-divergence listed in Table (ref).\footnote{It is worth noting that the family of Cressie-Read divergences nests the other three as special cases. Squared Hellinger corresponds to $\frac{1}{2} f_{1/2}$. L'H\^opital's rule shows that Kullback-Leibler corresponds to $\lim_{\gamma \rightarrow 1} f_\gamma$ and Reverse Kullback-Leibler to $\lim_{\gamma \rightarrow 0} f_\gamma$. See broniatowski2012divergences for additional discussion.} Precise assumptions regarding the $f$-divergence are collected in Assumption (ref) below.

Breakdown analysis in models identified with GMM

Suppose a preliminary analysis supports an alternative hypothesis $H_1$ over a null hypothesis $H_0$. For example, such an analysis may be based on the complete observations assuming MCAR, or using imputation and assuming $Y$ is MAR conditional on $X$. The breakdown point is the minimum amount of selection needed to overturn such a conclusion. When selection is measured in terms of the squared Hellinger distance, the breakdown point translates the claim that $H_0$ is true into a claim about the ability of $Z$ to predict $D$. Specifically, if $H_0$ were true then $1 - \frac{E[\sqrt{\text{Var}(D \mid Z)}]}{\sqrt{\text{Var}(D)}}$ would be weakly larger than the breakdown point. If this is implausible, then $H_0$ is similarly implausible.

This section formalizes this idea for models identified with the generalized method of moments (GMM). Suppose the parameter of interest $\beta \in \textbf{B} \subseteq \mathbb{R}^{d_b}$ is characterized as the unique solution to a finite set of moment conditions,

where the expectation is taken with respect to the unconditional distribution, $p_D P_1 + (1-p_D) P_0$. The conclusion to be investigated is that $\beta$ falls outside a particular set $\textbf{B}_0 \subset \textbf{B}$, motivating the null and alternative hypotheses

Recall that the observed data is $\{(D_i, D_i Y_i, X_i)\}_{i=1}^n$, where $D_i = \mathbbm{1}\{Y_i \text{ is observed}\}$. The sample identifies $P_1$, $p_D$, and $P_{0X}$. A hypothetical distribution of the incomplete observations $Q$ rationalizes the parameter $b$ if it has the identified marginal distribution of $X$, $Q_X = P_{0X}$, and the implied unconditional distribution $p_D P_1 + (1-p_D) Q$ solves the moment conditions for $b$. The set of such distributions implying finite selection is

The breakdown point $\delta^{BP}$ is the minimum selection needed to rationalize the null hypothesis:

where the infimum over the empty set is understood to be infinity. A simple example illustrates the idea.

Breakdown analysis can also be framed as an exercise in partial identification, as in kline2013sensitivity, masten2020inference, and diegert2025assessing. In this framing, the researcher considers assumptions of the form $d_f(P_0, P_1) \leq \delta$ for some $\delta > 0$, which continuously relax the assumption $P_0 = P_1$. The identified set for $\beta$ grows with $\delta$. As long as the identified set is a subset of $\textbf{B} \setminus \textbf{B}_0$, it is clear the researcher's conclusion holds. The breakdown point $\delta^{BP}$ can then be defined as either the largest $\delta$ for which the identified set is contained in $\textbf{B} \setminus \textbf{B}_0$, or the smallest $\delta$ for which the identified set has nontrivial intersection with $\textbf{B}_0$ (the latter of which corresponds to the definition given in display (ref)). For further discussion of this equivalent framing of the breakdown point, see the Supplementary Material.

The remainder of this paper constructs a $\sqrt{n}$-consistent and asymptotically normal estimator of $\delta^{BP}$, and constructs a lower confidence interval for $\delta^{BP}$. Researchers working with partially complete datasets should discuss the plausible amount of selection in their setting, and report the point estimate and the lower confidence interval for $\delta^{BP}$ for each asserted conclusion. This will make plain to readers which conclusions are more sensitive to missing data concerns, and whether crucial results are sufficiently robust.

Preview of results

Estimation of $\delta^{BP}$ proceeds by separating the optimizations in (ref). Define the primal problem

and notice that $\delta^{BP} = \inf_{b \in \textbf{B}_0} \nu(b)$. The first step is to estimate the value function $\nu$ over a set $B \subseteq \textbf{B}$ large enough that $\inf_{b \in \textbf{B}_0} \nu(b) = \inf_{b \in B \cap \textbf{B}_0} \nu(b)$, while the second step estimates $\delta^{BP}$ through a simple plug-in estimator.

The primal problem is an infinite dimensional convex optimization problem over the space of probability distributions, but one that is very well studied in convex analysis. In particular, when $\textbf{P}^b$ defined in display (ref) is characterized by a finite number of moment conditions, the primal problem has a well behaved, finite-dimensional dual problem with the same value function borwein1991duality, borwein1993partially, csiszar1999mem, broniatowski2006minimization. Section (ref) discusses this dual problem and the assumptions needed to make use of it. Under regularity conditions discussed in Section (ref), sample analogue estimators of $\nu$ based on this dual problem are uniformly consistent and asymptotically Gaussian on compact subsets of the parameter space. Differentiability of the infimum then implies convergence in distribution of the plug-in estimator.

To conclude this section, Assumption (ref) collects conditions on the setting, the GMM model, and the $f$-divergence used to measure selection.

The finite support condition, assumption (ref) (ref), simplifies estimation and inference by ensuring that $\textbf{P}^b$ is characterized by a finite number of moments. This assumption is not needed to define the breakdown point, but is used to ensure the duality results of section (ref) hold and can be used to define the estimator. This assumption can also be relaxed. In settings where assumption (ref) (ref) fails, researchers can still perform a breakdown point analysis through a conservative procedure described in remark (ref) below.

Condition (ref) ensures the $f$-divergence used to measure selection is well behaved, and is satisfied by every divergence in Table (ref). In particular, strict convexity of $f$ ensures the primal problem (ref) has a unique solution ($P_1$-almost surely). $f$ is required to be essentially smooth to ensure the dual problem has a unique solution. The requirements that $f(x)$ take a unique minimum of $0$ at $x=1$ and $f(x) = \infty$ for $x < 0$ ensures that $d_f(Q \Vert P)$ is a well-defined $f$-divergence.

Duality

As defined in display (ref), $\nu(b)$ is the value function of an infinite dimensional convex optimization problem. Fortunately, when selection is measured with an $f$-divergence, this minimization becomes a well-studied problem known by various names: maximal entropy csiszar1999mem, partially finite programming borwein1991duality, or $f$-divergence projection broniatowski2006minimization. The convex analysis results in these papers connect the primal problem in display (ref) to a finite dimensional dual problem that is much easier to study and estimate. Under mild conditions, the value function of this dual problem coincides with the value function of the primal.

To state the dual problem, first note that the primal can be viewed as a problem over the set of densities with respect to $P_1$:

where

As shown in borwein1991duality, the dual problem corresponding to (ref) is given by

where $f^*$ is the convex conjugate of $f$, given by $f^*(r) \equiv \sup_{t \in\mathbb{R}}\{rt - f(t)\}$. For convenience, table (ref) summarizes the convex conjugate for several common divergences.

Weak and strong duality

Assumption (ref) suffices to show $V(b) \leq \nu(b)$. This fact is known as weak duality, and implies that

for any $B \subseteq \textbf{B}$. This inequality shows that using the dual problem for estimation of the breakdown point is at worst conservative. If $\inf_{b \in B \cap \textbf{B}_0} V(b)$ is large enough to assuage selection concerns, researchers are assured that the breakdown point can only be larger.

Assuming only slightly more ensures strong duality holds, that is, $V(b) = \nu(b)$. Recall from Assumption (ref) (ref) that the interior of $\text{dom}(f) = \{t \in \mathbb{R} \; : \; f(t) < \infty\}$ is denoted $(\ell, u)$.

That strong duality holds under these conditions is a well-known result.\footnote{To the authors knowledge, the first to show strong duality holds under similar conditions was borwein1991duality. The proof of theorem (ref), found in the Supplementary Material, uses a result due to csiszar1999mem.}

The first order condition of the dual problem (ref) provides intuition. Exchanging expectation and differentiation, the first order condition is

where $\lambda(b) \in \mathbb{R}^{d_g + K}$ solves the dual problem. Consider $(f^*)'(\lambda(b)^\intercal h(y, x, b))$ as a density with respect to $P_1$. Notice that the first $d_g$ equations of the first order condition ensure $p_D E_{P_1}[g(Y, X, b)] + (1-p_D) E_{P_1}[(f^*)'\left(\lambda(b)^\intercal h(Y, X, b)\right)g(Y, X, b)] = 0$, while the remaining $K$ equalities ensure the marginal distribution of $X$ matches $P_{0X}$. In fact, the proof of theorem (ref) shows that under assumptions (ref) and (ref), $(f^*)'\left(\lambda(b)^\intercal h(y, x, b)\right)$ is the $P_1$-density of the solution to the primal problem.

Assumption (ref) ensures the set on which $\nu$ is estimated is large enough to estimate the breakdown point, but not so large as to contain parameter values that cannot be rationalized with a well behaved $P_1$-density. To illustrate, consider again example (ref). $Y$ is a scalar, $\beta = E[Y] = p_D E_{P_1}[Y] + (1-p_D) E_{P_0}[Y]$, and $P_1$ is $\mathcal{U}[0,1]$. For tractability suppose that Kullback-Leibler is used to measure selection. Since $P_0$ takes values on $[0,1]$, the Manski bounds for $\beta$ are $\left[\frac{p_D}{2}, 1 - \frac{p_D}{2}\right]$. The Supplementary Material shows that strong duality is satisfied whenever $b \in \left(\frac{p_D}{2}, 1 - \frac{p_D}{2}\right)$. Thus for this example, $B$ can be any convex, compact set in the interior of the Manski bounds.

Estimation

Assumptions (ref) and (ref) are maintained throughout the remainder of the paper. Accordingly, the notation $\nu$ will be used for the value function of the dual problem as well.

The estimator

The sample analogue of the dual problem provides an estimator of the value function, and suggests a simple plug-in estimator of the breakdown point. The asymptotic properties of these estimators are easier to study if the objective of the dual problem is expressed with a single unconditional expectation, which comes at the cost of additional notation.

First define the matrix $J(D) =

$ where $I_{d_g}$ and $I_K$ are identity matrices. Notice that $E\left[\frac{J(D) h(DY, X, b)}{(1-p_D)}\right] = c(b)$ and

Define

and observe that the dual problem is $\sup_{\lambda \in \mathbb{R}^{d_g + K}} E[\varphi(D, DY, X, b, \lambda, p_D)]$. The estimator of the value function is defined pointwise by

where $\hat{p}_{D,n} \equiv \frac{1}{n}\sum_{i=1}^n D_i$ estimates $p_D$. Finally, $\hat{\delta}_n^{BP} \equiv \inf_{b \in B \cap \textbf{B}_0} \hat{\nu}_n(b)$ estimates the breakdown point.

Asymptotic normality

The following assumption suffices for $\hat{\delta}_n^{BP}$ to be $\sqrt{n}$-consistent and asymptotically normal. First observe that the estimands $\theta_0(b) = (\nu(b), \lambda(b), p_D)$ solve the moment conditions $E[\phi(D, DY, X, b, \theta_0(b))] = 0$, where

Let $\text{Gr}(\theta_0) \equiv \left\{(b, \theta_0(b)) \; : \; b \in B\right\}$ denote the graph of $\theta_0$. For $\eta > 0$, the closed $\eta$-expansion about this graph is $\text{Gr}(\theta_0)^\eta \equiv \left\{(b, \theta) \in B \times \mathbb{R}^{d_g + K + 2} \; : \; \inf_{(b', \theta') \in \text{Gr}(\theta_0)} \lVert (b,\theta) - (b', \theta') \rVert \leq \eta\right\}$.

As previewed in section (ref), $\hat{\delta}_n^{BP}$ is viewed as a two-step estimator where $\hat{\nu}_n$ estimates $\nu$ in the first step, and $\hat{\delta}_n^{BP} = \inf_{b \in B \cap \textbf{B}_0} \hat{\nu}_n(b)$ is a plug-in estimator for $\delta^{BP} = \inf_{b \in B \cap \textbf{B}_0} \nu(b)$. Conditions (ref), (ref), and (ref) imply $\sqrt{n}(\hat{\nu}_n - \nu)$ converges weakly in the space of bounded functions on $B$, to a limiting process that is almost surely continuous. This is shown by linearizing $0 = \frac{1}{n}\sum_{i=1}^n \phi(D_i, D_i Y_i, X_i, b, \hat{\theta}_n(b))$ uniformly over $b \in B$. Conditions (ref) and (ref) ensure minimization over $B \cap \textbf{B}_0$ is a (Hadamard) differentiable map on the set of continuous functions of $B$. The delta method then implies $\sqrt{n}(\hat{\delta}_n^{BP} - \delta^{BP})$ converges in distribution to a normal distribution.

Assumption (ref) (ref) and (ref) are easily verified by inspection of $\textbf{B}_0$ and $g$ respectively. Conditions (ref) and (ref) are similar to conditions required of generalized empirical likelihood estimators (see, e.g., antoine2021robust assumption 1 (v) and assumption 3 (iv), (vii)). Assumption (ref) (ref) deserves additional scrutiny. When $\textbf{B}_0$ is a convex set, condition (ref) holds when $\nu$ is a strictly convex function. The following lemma shows that this is the case when $g(y,x,b)$ describes a linear model with the outcome being the only missing data value.

Lemma (ref) covers instrumental variable models directly, and ordinary least squares as a special case (by setting $X_2 = X_1$). It also covers parameters of the form $\beta = E[\tilde{g}(Y, X)]$, because the OLS regression of $\tilde{g}(Y,X)$ on a constant recovers $E[\tilde{g}(Y, X)]$. Simulation evidence presented in the Supplementary Material suggests data generating processes and models not covered by lemma (ref) also produce convex $\nu$. Remark (ref) below discusses an approach to relaxing assumption (ref) (ref), at the cost of additional complexity.

Theorem (ref) below formally states the convergence in distribution result along with consistency of an estimator of the asymptotic variance. The variance depends on the Jacobian term $\Phi(b) \equiv E[\nabla_\theta \phi(D, DY, X, b, \theta_0(b))]$, which is estimated with

where $\hat{\theta}_n(b) \equiv (\hat{\nu}_n(b), \hat{\lambda}_n(b), \hat{p}_{D,n})$ and $\hat{\lambda}_n(b) \equiv \operatorname*{arg\,max}_{\lambda \in \mathbb{R}^{d_g + K}} \frac{1}{n}\sum_{i=1}^n \varphi(D_i, D_i Y_i, X_i, b, \lambda, \hat{p}_{D,n})$.

Inference

A large breakdown point implies the incomplete distribution $P_0$ would have to differ greatly from $P_1$ to rationalize the null hypothesis. If $\delta^{BP}$ is larger than the plausible amount of selection in the setting, the null hypothesis is similarly implausible. Skeptical readers following this argument may worry the point estimate $\hat{\delta}_n^{BP}$ is larger than $\delta^{BP}$ due to sample noise -- but the force of the argument is only strengthened if $\hat{\delta}_n^{BP}$ falls below $\delta^{BP}$.

To address these concerns, researchers should report lower confidence intervals along with point estimates of the breakdown point. Theorem (ref) implies that under assumptions (ref), (ref), and (ref),

satisfies $\lim_{n \rightarrow \infty} P(\widehat{CI}_{L,n} \leq \delta^{BP}) = 1-\alpha$ when $c_{1-\alpha}$ is the $1-\alpha$ quantile of the standard normal distribution.

Simulations

This section presents simulation results on a variety of different data generating processes. This serves both to illustrate the wide scope of models which can make use of breakdown point analysis and to investigate the finite sample properties of the proposed estimators. In each case, selection is measured using squared Hellinger divergence.

Expectation

Recall example (ref). The parameter of interest is the mean of a scalar random variable $Y$, $\beta = E[Y] = p_D E_{P_1}[Y] + (1-p_D)E_{P_0}[Y]$, and the sample is $\{D_i, D_iY_i\}_{i=1}^n$. The distribution of $Y \mid D =1$ is the uniform distribution on $[0,1]$. The probability of observing $Y$ is $p_D = P(D= 1) = 0.7$. To support the claim $H_1 \; : \; \beta > 0.4$, let $H_0 \; : \; \beta \leq 0.4$. Recall that the true breakdown point, $\delta^{BP}$, of this example is just over $0.2$.

The following table summarizes 1,000 simulations for several different sample sizes.\footnote{Here $\text{CI Length} \equiv \hat{\delta}_n^{BP} - \widehat{CI}_{L,n}$.}

The simulations show little bias. Coverage is slightly above the targeted 95 percent significance level in smaller samples.

Linear model

Linear models are the among the most common tools used by empirical researchers. This subsection uses simulations to investigate linear regression with exogenous regressors.

Consider the model

where $W =

^\intercal$ are the exogenous regressors: $E[W\varepsilon] = 0$. Here $Y_1$ is a continuously distributed dependent variable, $X_1 = \{0,1\}$ is the regressor of interest, $Y_2$ is a continuously distributed regressor, and $X_2 \in \{0, 1, 2\}$ is a discrete regressor. The conclusion to be investigated is that the coefficient on $X_1$ is positive:

The data generating process specification takes inspiration from Mincerian wage equations. For worker $i$, let $Y_{1i}$ be $i$'s log-income, $X_{1i}$ an indicator for $i$ being a college graduate, $Y_{2i}$ be $i$'s work experience, and $X_{2i}$ the number of parents with college degrees ($0$, $1$, or $2$). Specifically, let $X_2$ be multinomial, $X_1 \sim \text{Binomial}\left(\frac{X_2 + 1}{4}\right)$, and $Y_2 \sim \text{Beta}(3-X_1, 3)$.\footnote{The distribution of $X_2$ is $P(X_2 = 0) = 0.4$, $P(X_2 = 1) = 0.25$, and $P(X_2 = 2) = 0.35$.} Let $\tilde{\varepsilon} \sim U[-1,1]$ (independent of all other variables), and $\varepsilon = (X_1+1)\tilde{\varepsilon}$. The coefficients are specified as $\beta_0 = \beta_1 = \beta_2 = 1$ and $\beta_3 = 0.5$. Finally, $Y_1$ is generated according to equation (ref). Notice the support of $(Y_1, Y_2, X_1, X_2)$ is compact, ensuring the moment conditions in assumption (ref) (ref) are satisfied. The Supplementary Material shows simulation evidence that this model and data generating process produces a convex $\nu(\cdot)$, suggesting that assumption (ref) (ref) holds.

For the missing data process, let $D = \mathbbm{1}\{\varepsilon X_1 + 10 X_1 + 5 (X_2 - 1) > \eta\}$, where $\eta \sim N(-5, 15^2)$. The population value of the breakdown point is approximated as the point estimate obtained from a sample with one million observations. This sample reveals $P(D = 1)$ is about $0.71$, and suffers from selection. Specifically, ignoring the incomplete observations is equivalent to solving $\frac{1}{n}\sum_{i=1}^n \frac{D_i}{\hat{p}_{D,n}} (Y_i - W_i^\intercal \hat{\beta}_n^{MCAR})W_i = 0$ for $\hat{\beta}_n^{MCAR}$, which results in $\hat{\beta}_n^{MCAR} = (1.08, 1.34, 1.02, 0.39)$. The squared Hellinger distance between $P_{0X}$ and $P_{1X}$ is about $0.08$. This large sample suggests the breakdown point of the conclusion $\beta_1 > 0$ is about 0.163, which is treated as the truth when evaluating the 1,000 simulations per sample size summarized in the following table:

The simulations again show little bias, with coverage slightly above the targeted 95 percent significance level in smaller samples.

Logit model

The logit model is a popular choice for estimating the conditional probability of an event. Let $Z = (Z_1, Z_{-1}) \in \{0,1\} \times \mathbb{R}^d$ and suppose that $P(Z_1 = 1 \mid Z_{-1}) = \Lambda(Z_{-1}^\intercal \beta)$, where $\Lambda(t) \equiv \frac{\exp(t)}{1 + \exp(t)}$. The log-likelihood is concave, so estimating this model through maximum likelihood is equivalent to solving the first order condition

The model can be viewed as nonlinear GMM, with moment function $g(z,b) = (z_1 - \Lambda(z_{-1}^\intercal b))z_{-1}$. The conclusion to be investigated is that $P(Z_1 = 1 \mid Z_{-1} = \bar{z}) = \Lambda(\bar{z}^\intercal \beta)$ is at least 0.5 for a fixed $\bar{z}$ of interest. The corresponding null and alternative hypotheses are

The data generating process is one where the dependent variable is always observed, and the regressors are sometimes missing. Specifically, $Y = Z_{-1} \in \mathbb{R}^3$ is constructed by drawing $\tilde{Y} \sim N(0, \Omega)$ and setting $Y^{(j)} = 2 \times (\Phi(\tilde{Y}^{(j)}) - 0.5)$ for each $j = 1, 2, 3$; the result is that each $Y^{(j)}$ has uniform marginal distribution on $[-1, 1]$, and together $(Y^{(1)}, Y^{(2)}, Y^{(3)})$ have a nontrivial joint distribution.\footnote{The matrix $\Omega$ is described by $\text{Var}(Y^{(j)}) = 1$ for each $j=1,2,3$, $\text{Cov}(Y^{(1)}, Y^{(2)}) = 0.5$, $\text{Cov}(Y^{(1)}, Y^{(3)}) = -0.1$, and $\text{Cov}(Y^{(2)}, Y^{(3)}) = 0.3$} The outcome is always observed: $X = Z_1$. The true underlying coefficients are $\beta = (1, -1, 0.1)$. Once again, the compact support of $(X, Y)$ ensures the moment conditions in assumption (ref) (ref) are satisfied. Simulation evidence presented in the Supplementary Material suggests that this model and data generating process produces a convex value function. Since $H_0$ is equivalent to $\bar{z}^\intercal \beta \leq \ln(0.5) - \ln(1 - 0.5) = 0$ and therefore defines a convex $\textbf{B}_0$, this suggests that assumption (ref) (ref) holds.

The missing data process is conditionally binomial with $P(D = 1 \mid X = x, Y = y) = \max\{0.8 - X, Y^{(3)}/2 + 0.5 \}$; that is, the probability of an observation being complete is at least $0.8$ when $X = 0$ and grows weakly with $Y^{(3)}$. The resulting samples suffer from selection. A sample with one million observations suggests that $P(D = 1)$ is about $0.65$. Ignoring the incomplete observations is equivalent to solving $\frac{1}{n}\sum_{i=1}^n \frac{D_i}{\hat{p}_{D,n}} g(D_i Y_i, X_i, \hat{\beta}_n^{MCAR}) = 0$, which results in $\hat{\beta}_n^{MCAR} = (1, -1, 0.79)$. The estimated squared Hellinger distance between $P_{0X}$ and $P_{1X}$ is $0.076$. The covariate value of interest is $\bar{y} = (-0.35, -0.25, 0.5)$. The true value for $\Lambda(\bar{y}^\intercal \beta)$ is $0.488$, while the estimate using the complete observations of the large sample above is $\Lambda(\bar{y}^\intercal \hat{\beta}_n^{MCAR}) = 0.573$. The point estimate for the breakdown point of the conclusion described by (ref) using this large sample is $0.108$. This is treated as the truth when evaluating the 1,000 simulations per sample size summarized in the following table:

These simulations show essentially zero bias and correct coverage at relatively small sample sizes.

Application: attrition in randomized controlled trials

This section reports estimates of the breakdown point of conclusions drawn from a number of randomized controlled trials (RCTs) conducted in developing countries. There are several advantages to demonstrating breakdown point analysis on real world data in this way. First, these studies are known to suffer from missing data that is unlikely to be missing at random. Second, missingness in these studies is often a result of study subject migration. This illustrates an important point discussed in section (ref): extra scrutiny should be given to conclusions involving variables that would predict migration, and hence missingness. Finally, these RCTs are evaluated using similar methodologies. The breakdown point estimates below thus show a range of values that might be expected due to variation in real world data, rather than significant variation in methodology. This provides useful context for researchers using breakdown point analysis to investigate the robustness of conclusions drawn from similar studies.

Missing data due to attrition is a prominent concern for studies conducted in developing countries. thomas2012cutting notes that the primary cause of this attrition is researchers being unable to find respondents who moved after the baseline survey. The authors study the Indonesia Family Life Survey, which has a notably low attrition rate despite high mobility of the target population, and provide evidence that migrants differ from non-migrants along dimensions unlikely to be observed at a survey's baseline. Attrition due to study subject migration is also noted in molina2025attrition, which studies randomized controlled trials conducted in developing countries. The authors show through examples that attempting to correct for attrition through inverse propensity weighting with baseline data does not make a notable difference to estimates -- but including individuals found only after intensive (and often costly) tracking does. Both papers suggest that missing data in these contexts is likely due to migration, and not missing at random.

The breakdown point estimates below pertain to results found in barham2024experimental, bandiera2020women, and giacobino2024schoolgirls, which all study randomized controlled trials conducted in developing countries. barham2024experimental studies a conditional cash transfer (CCT) implemented by the Nicaraguan government to address poverty by improving health and education. Study subjects were randomized into early or late treatment groups in the baseline year, 2000, and the authors study the differential effects of receiving the treatment early. The breakdown point analysis below focuses on conclusions regarding the cohort of boys who were aged 9-12 in the year 2000. Those who received the CCT early -- when they were at higher risk of dropping out -- showed higher labor market participation and earnings in 2010. bandiera2020women studies the impact of a program in Uganda designed to increase women's empowerment through training in vocational and life skills. The breakdown point analysis focuses on five conclusions regarding outcomes measured at midline: the treatment increased an index of entrepreneurial ability, increased the probability of being engaged in any income-generating activity, increased the probability of being self-employed, increased the probability of being employed for a wage, and increased expenditures on goods in the last month. Finally, giacobino2024schoolgirls studies a scholarship for adolescent girls in Niger to attend middle school. The intervention was designed to deter child marriage. The breakdown point analysis below focuses on three conclusions: the scholarship reduced the probability of dropping out, reduced the probability of being married by endline, and increased life satisfaction as measured by a standardized 10-point Likert scale.

Table (ref) reports the intent-to-treat (ITT) estimates when incomplete observations are dropped, referred to as missing completely at random (MCAR) estimates. Column (1) reports the total sample size, including subjects that attrited and could not be included in the estimates. Column (2) reports the number of complete observations on which the subsequent estimates are based. Column (3) reports the average of the outcome among untreated subjects. Columns (4) and (5) report coefficient estimates on an indicator for treatment status from a regression of the outcome on a constant, the indicator for treatment, and additional regressors. Column (4) replicates the results from the original papers, including the authors' choice of additional regressors and standard errors. To facilitate comparisons across studies, estimates in column (5) use indicators for the subject's region at baseline as the additional regressors and HC3 standard errors. Consistent with independent randomization of treatment assignment, changing the additional regressors does not meaningfully alter the estimates.

Table (ref) reports breakdown point analyses. In each case, the conclusion is that the ITT parameter takes the sign implied by the point estimate in table (ref). Squared Hellinger is used to measure selection. Every subject's treatment status and region at baseline is observed, and used as the variables in $X$. Also reported is an estimate of $H^2(P_{0X}, P_{1X})$, formed by taking sample analogues of $P(X=x)$ for each possible value in $\mathcal{X}$ and plugging these into the definition of squared Hellinger. This provides an estimated lower bound on the amount of selection in the given setting, as described in remark (ref).

Tables (ref) and (ref) show a number of patterns worth emphasizing. Compared to the other two studies, the larger sample of bandiera2020women resulted in lower standard errors in table (ref) and lower confidence intervals that are closer to the breakdown point estimates in table (ref). However, the larger share of incomplete observations in this study results in generally smaller breakdown points. The magnitude of MCAR estimates and the amount of missing data are clearly determinants of the size of the breakdown point, but superficially similar results can have quite different breakdown points. For example, consider the conclusion that the CCT from barham2024experimental differentially raised the probability of working somewhere other than the recipient's family farm, and the claim that the scholarship studied in giacobino2024schoolgirls reduced the probability the subject is being married at endline. The corresponding MCAR estimates have a similar magnitude ($0.06$ compared to $-0.07$) and have similar standard errors ($0.02$ or $0.03$, depending on the specification). The two studies have a similar share of incomplete data. However, the breakdown point of the latter result appears notably larger than that of the former result. Notice also that the estimated effects of the treatment from bandiera2020women on self-employment and wage employment under MCAR differ considerably ($0.06$ compared to $0.01$), but the claims that the corresponding ITT estimates are positive have similar breakdown points.

Several results appear fragile, while others are quite robust. Consider the claims that the treatment studied in bandiera2020women increased expenditures or the probability of being self- or wage- employed. The breakdown point estimates of these claims are quite close to the estimate of $H^2(P_{0X}, P_{1X})$, implying that it would take only a small amount of selection on these outcomes to rationalize such claims being false. In contrast, the results of giacobino2024schoolgirls appear quite robust. The claim that the scholarship reduced the probability of dropping out could not be rationalized as false in the sample. The other claims from giacobino2024schoolgirls have point estimates for the breakdown point that are close to 0.2. As discussed in remark (ref), these are large values for a breakdown point implying the results are quite robust.

Conclusion

This paper proposes breakdown point analysis as a tractable approach to assessing the sensitivity of a researcher's conclusion to the common assumption that the data are missing at random. When defined with squared Hellinger, the breakdown point $\delta^{BP}$ has a natural interpretation: if the result were false, the variables under study ($Z$) would have to predict an observation being selected into the sample ($D$) at least well enough that $H^2(P_0, P_1) = 1 - E[\sqrt{\text{Var}(D \mid Z)}]/\sqrt{\text{Var}(D)} \geq \delta^{BP}$. Estimators based on the sample analogue of the dual problem are shown $\sqrt{n}$-consistent and asymptotically normal, which facilitates the construction of lower confidence intervals. Researchers working with incomplete datasets should report the breakdown point estimate and lower confidence interval along with standard results, making transparent to their audience how robust the conclusion is to relaxing the assumption that the data are missing at random.

\nocite{bandiera2020womendata} \nocite{barham2024experimentaldata} \nocite{giacobino2024schoolgirlsdata} \nocite{molina2025attritiondata}

\singlespacing