The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
86,812 characters
Robustness to Missing Data: Breakdown Point Analysis
\maketitle
\begin{abstract}
Missing data is pervasive in econometric applications, and rarely is it plausible that the data are missing (completely) at random. This paper proposes a methodology for studying the robustness of results drawn from incomplete datasets. Selection is measured as the divergence from the distribution of complete observations to the distribution of incomplete observations. The \textit{breakdown point} is defined as the minimal amount of selection needed to overturn a given result. Reporting point estimates and lower confidence intervals of the breakdown point is a simple, concise way to communicate the robustness of a result. An estimator of the breakdown point is proposed and shown $\sqrt{n}$-consistent and asymptotically normal. This estimator can be applied directly to conclusions drawn from any model identified with the generalized method of moments (GMM) that satisfies mild assumptions. Simulations demonstrate the finite sample performance of the breakdown point estimator on averages, linear regression, and logistic regression. The methodology is illustrated by estimating the breakdown point of conclusions drawn from several randomized controlled trails suffering from missing data due to attrition.
\end{abstract}
\bigskip
\noindent \textbf{Keywords:} Missing data, generalized method of moments, robustness, sensitivity analysis.
\noindent \textbf{JEL classification:} C01, C14, C18, C21, C25.
\pagenumbering{gobble}
\clearpage
\pagenumbering{arabic}
\section{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.
\cite{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 \citep{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 \textit{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) \citep{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 \citep{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 \citep{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 \citep{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 \cite{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 \citep{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 \citep{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 \cite{horowitz1995identification} in the context of corrupted data.
One well-known example of breakdown analysis is \cite{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 \cite{oster2019unobservable}, is widely used in empirical economics, and is an area of active research; see, e.g., \cite{masten2025effect}, \cite{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 \cite{masten2020inference}, \cite{bonvini2022sensitivity}, \cite{rosenbaum2023sensitivity}, and \cite{spini2024robustness}.
This paper is not the first to notice the appeal of breakdown point analysis in the context of missing data.
\cite{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.
\cite{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{Section: simulations} and application in section \ref{Section: application} 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{Section: measuring selection and breakdown analysis, divergences}, have seen widespread use in the sensitivity analysis literature.
\cite{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.
\cite{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 \cite{bonhomme2022minimizing}, sensitivity to the choice of a Bayesian prior in \cite{ho2023global}, and breakdown of causal inference results in \cite{spini2024robustness}.
The remainder of this paper is structured as follows.
Section \ref{Section: measuring selection and breakdown analysis} formalizes the setting, the proposed measure of selection, and the breakdown point.
The dual problem is presented and discussed in section \ref{Section: duality}.
Section \ref{Section: estimators and asymptotics} defines the estimator and states the main results on estimation and inference, which are proven in the Supplementary Material.
Section \ref{Section: simulations} presents a simulation study investigating the finite sample performance of the estimators.
The methodology is illustrated in section \ref{Section: application}, which contains estimates of the breakdown point of conclusions drawn from several randomized controlled trails suffering from attrition.
Section \ref{Section: conclusion} concludes.
\section{Measuring selection and breakdown analysis}
\label{Section: 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 \textit{complete case} and \textit{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{Remark: absolute continuity of P_0 wrt P_1} 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{Remark: always observed variables have finite support}.
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 \textit{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.
\subsection{An interpretable measure of selection}
\label{Section: measuring selection and breakdown analysis, squared hellinger}
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, \cite{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 \cite{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:
\begin{align*}
&f_1(z) = \frac{P(D = 1 \mid Z = z)}{p_D}, &&f_0(z) = \frac{1-P(D=1 \mid Z = z)}{1-p_D}
\end{align*}
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
\begin{equation*}
H(Q, P) \equiv \left(\frac{1}{2} \int \left(\sqrt{\frac{dQ}{d\lambda}(z)} - \sqrt{\frac{dP}{d\lambda}(z)}\right)^2 d\lambda(z) \right)^{1/2}
\end{equation*}
where $\lambda$ is any measure dominating both $P$ and $Q$.}
\begin{restatable}{lemma}{lemmaSquaredHellingerInterpretation}
\label{Lemma: squared hellinger interpretation}
\singlespacing
Let $(Z, D) \in \mathbb{R}^{d_z} \times \{0,1\}$ be random variables with $p_D = P(D=1) \in (0,1)$. Let $Z \mid D = 1 \sim P_1$ and $Z \mid D = 0 \sim P_0$. Then
\begin{equation}
H^2(P_0, P_1) = 1 - \frac{E\left[\sqrt{\text{Var}(D \mid Z)}\right]}{\sqrt{\text{Var}(D)}} \label{Display: squared Hellinger interpretation}
\end{equation}
where the expectation is taken with respect to $p_D P_1 + (1-p_D) P_0$, the marginal distribution of $Z$.
\end{restatable}
All results are proven in the Supplementary Material. Equation \eqref{Display: squared Hellinger interpretation} 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 \eqref{Display: squared Hellinger interpretation} 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.
\begin{remark}
\label{Remark: squared hellinger lower bound}
It is shown in the Supplementary Material that
\begin{equation*}
H^2(P_0, P_1) = 1 - \frac{E\left[\sqrt{\text{Var}(D \mid Y, X)}\right]}{\sqrt{\text{Var}(D)}} \geq 1 - \frac{E[\sqrt{\text{Var}(D \mid X)}]}{\sqrt{\text{Var}(D)}} = H^2(P_{0X}, P_{1X})
\end{equation*}
where $P_{0X}$, $P_{1X}$ are the marginal distributions of $X$ conditional on $D = 0$ and $D = 1$ respectively. This lower bound on selection is identified from the sample, and motivates the common practice of comparing the distribution of $X$ conditional on $D=0$ with that of $X$ conditional on $D=1$; the distributions $P_0$ and $P_1$ can only be ``further'' apart.
\end{remark}
\subsection{Divergences}
\label{Section: measuring selection and breakdown analysis, 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 \textit{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 \textit{$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
\begin{align}
d_f(Q \Vert P) \equiv \begin{cases}
\int f\left(\frac{dQ}{dP}\right)dP & \text{ if } Q \ll P \\
\infty & \text{ otherwise }
\end{cases}, \label{Definition: f divergence}
\end{align}
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 \textit{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$.}
\begin{table}[H]
\begin{small}
\caption{Common $f$-divergences}
\label{Table: Common f divergences}
\begin{center}
\begin{tabular}{|l || l | l | l |} \hline
Name & Common formula & $f(t)$ & $\text{dom}(f)$ \\ \hline \hline
Squared Hellinger & $H^2(Q, P) = \frac{1}{2}\int \left(\sqrt{\frac{dQ}{dP}(z)} - 1\right)^2 dP(z)$ & $f(t) = \frac{1}{2}(\sqrt{t} - 1)^2$ & $t \in [0, \infty)$ \\ \hline
Kullback-Leibler (KL) & $KL(Q \Vert P) = \int \log\left(\frac{dQ}{dP}(z)\right)dQ(z)$ & $f(t) = t \log(t)-t+1$ & $t \in [0,\infty)$ \\ \hline
``Reverse'' KL & $ KL(P \Vert Q) = \int \log\left(\frac{dP}{dQ}(z)\right)dP(z)$ & $f(t) = -\log(t)+t-1$ & $t \in (0, \infty)$ \\\hline
Cressie-Read & -- & $f_\gamma(t) = \frac{t^\gamma - \gamma t + \gamma -1}{\gamma(\gamma - 1)}$ & -- \\ \hline
\end{tabular}
\end{center}
\end{small}
\end{table}
Although squared Hellinger has intuitive appeal outlined in Section \ref{Section: measuring selection and breakdown analysis, squared hellinger}, the breakdown point analysis proposed in this paper remains tractable for any $f$-divergence listed in Table \ref{Table: Common f divergences}.\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 \cite{broniatowski2012divergences} for additional discussion.} Precise assumptions regarding the $f$-divergence are collected in Assumption \ref{Assumption: setting} below.
\begin{remark}
\label{Remark: absolute continuity of P_0 wrt P_1}
Measuring selection with an $f$-divergence facilitates estimation and inference, as the space of distributions $Q$ with $d_f(Q \Vert P_1) < \infty$ corresponds to the set of densities with respect to $P_1$. In essence, measuring selection with an $f$-divergence assumes that $P_0$ is absolutely continuous with respect to $P_1$, denoted $P_0 \ll P_1$, as all distributions failing this requirement have infinite divergence from $P_1$.
Absolute continuity is a natural assumption in some settings, but restrictive in others. For an example where $P_0 \ll P_1$ is natural, suppose $Y$ is a measure of time worked in a week measured in hours, obtained through a survey. Suppose the distribution $P_1$ displays positive mass at $Y = 0$ hours and $Y = 40$ hours, and is otherwise continuous. The assumption $P_0 \ll P_1$ here is natural, as it allows $P_0$ to have an atom of any size at $Y = 0$ and $Y = 40$ while ruling out distributions with positive mass at other points. For an example where absolute continuity fails, suppose $Y$ represents wage data where top coded observations are treated as missing. In this example, $P_1$ puts mass one below the top code while $P_0$ puts mass one above the top code, and $P_0 \not \ll P_1$.
\end{remark}
\subsection{Breakdown analysis in models identified with GMM}
\label{Section: measuring selection and breakdown analysis, breakdown analysis}
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,
\begin{equation*}
E[g(Z, \beta)] = 0 \in \mathbb{R}^{d_g}
\end{equation*}
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
\begin{align*}
&H_0 \; : \; \beta \in \textbf{B}_0, &&H_1 \; : \; \beta \in \textbf{B} \setminus \textbf{B}_0
\end{align*}
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$ \textit{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
\begin{equation}
\textbf{P}^b \equiv \left\{Q \; : \; Q \ll P_1, \; Q_X = P_{0X}, \; p_D E_{P_1}[g(Z, b)] + (1-p_D) E_Q[g(Z, b)] = 0\right\}. \label{Definition: Set of distributions rationalizing parameter}
\end{equation}
The \textit{breakdown point} $\delta^{BP}$ is the minimum selection needed to rationalize the null hypothesis:
\begin{equation}
\delta^{BP} \equiv \inf_{b \in \textbf{B}_0} \inf_{Q \in \textbf{P}^b} d_f(Q \Vert P_1), \label{Definition: breakdown point}
\end{equation}
where the infimum over the empty set is understood to be infinity. A simple example illustrates the idea.
\begin{example}
\label{Example: Squared Hellinger Expectation}
\singlespacing
Let $Y \in \mathbb{R}$ and $\beta = E[Y] = p_D E_{P_1}[Y] + (1-p_D)E_{P_0}[Y]$. Let $p_D = 0.7$ and $P_1$
be $\mathcal{U}[0,1]$. The claim to support is $H_1 \; : \; \beta > 0.4$, and selection is measured with squared Hellinger. $\textbf{P}^b$ is the set of continuous distributions on $[0,1]$ with expectation $\frac{b - p_D/2}{1-p_D}$, so that $Q \in \textbf{P}^b$ implies
\begin{equation*}
p_D E_{P_1}[Y] + (1-p_D)E_Q[Y] = \frac{p_D}{2} + (1-p_D) \frac{b - p_D/2}{1-p_D} = b
\end{equation*}
The inner minimization in display \eqref{Definition: breakdown point} chooses the distribution that minimizes selection while rationalizing $b$. The outer minimization chooses the parameter that minimizes selection while rationalizing $H_0 \; : \; \beta \leq 0.4$. Unsurprisingly, the outer minimization is solved by $b = 0.4$. The breakdown point $\delta^{BP}$ is slightly above $0.2$. A researcher convinced $H^2(P_0, P_1)$ is less than $0.2$ should conclude $\beta > 0.4$.
\begin{figure}[H]
\caption{$\inf_{Q \in \textbf{P}^b} d_f(Q \Vert P_1)$ and $p_D P_1 + (1-pD)Q^*$, where $Q^* \in \textbf{P}^{0.4}$ minimizes selection.}
\label{Figure: numerical example}
\vspace{0.5 cm}
\begin{center}
\includegraphics[scale=0.61, trim={0 0 0 83}, clip]{Fig_Uniform_SqHellinger_population_nub_bdp.png}
\end{center}
\end{figure}
\end{example}
\begin{remark}
\label{Remark: 0.2 is a big number}
Example \ref{Example: Squared Hellinger Expectation}, while simple, can help anchor expectations for the values of $\delta^{BP}$ when squared Hellinger is used to measure selection. It is clear from observation of the left panel of Figure \ref{Figure: numerical example} that the value of $\delta^{BP}$ would rise quite quickly as the $\beta$ in $H_0 : \beta \leq 0.4$ shrinks toward the lower Manksi bound of $0.35$. $\beta$ would need to be quite close to the theoretical bound to see breakdown point values of $0.5$ or higher. The right-hand side panel of Figure \ref{Figure: numerical example} shows the distribution $Q^*$ closest to $P_1$ in terms of squared Hellinger that rationalizes an unconditional mean of $0.4$. $Q^*$ is quite different from the uniform distribution $P_1$, and in many settings where the complete data follows a uniform distribution it would be implausible that the missing observations follow such a distinct distribution. This suggests a breakdown point of $0.2$ should be treated as quite large for a squared Hellinger breakdown point. Depending on the context, even smaller values could provide reassurance to many researchers.
\end{remark}
Breakdown analysis can also be framed as an exercise in partial identification, as in \cite{kline2013sensitivity}, \cite{masten2020inference}, and \cite{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 \eqref{Definition: breakdown point}). 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.
\subsection{Preview of results}
\label{Section: measuring selection and breakdown analysis, preview of results}
Estimation of $\delta^{BP}$ proceeds by separating the optimizations in \eqref{Definition: breakdown point}. Define the \textit{primal problem}
\begin{equation}
\nu(b) \equiv \inf_{Q\in \textbf{P}^b} d_f(Q \Vert P_1), \label{Definition: primal problem}
\end{equation}
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 \eqref{Definition: Set of distributions rationalizing parameter} 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 \citep{borwein1991duality, borwein1993partially, csiszar1999mem, broniatowski2006minimization}. Section \ref{Section: duality} discusses this dual problem and the assumptions needed to make use of it. Under regularity conditions discussed in Section \ref{Section: estimators and asymptotics}, 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{Assumption: setting} collects conditions on the setting, the GMM model, and the $f$-divergence used to measure selection.
\begin{restatable}[Setting]{assumption}{assumptionSetting}
\label{Assumption: setting}
\singlespacing
$\{(D_i, D_i Y_i, X_i)\}_{i=1}^n$ is an i.i.d. sample from a distribution satisfying
\begin{enumerate}[label=(\roman*)]
\item $p_D = P(D=1) \in (0,1)$, \label{Assumption: setting, missing data}
\item $X \mid D = 1$ and $X \mid D = 0$ have the same finite support $\{x_1, \ldots, x_K\}$, \label{Assumption: setting, support of always observed variables}
\item $E\left[\sup_{b \in \boldsymbol{B}} \lVert g(Z, b) \rVert \mid D = 1\right] < \infty$, where $Z = (Y,X)$, \label{Assumption: setting, GMM}
\item $P_0 \ll P_1$, and \label{Assumption: setting, absolute continuity}
\item $f : \mathbb{R} \rightarrow [0,\infty]$ is closed, proper, strictly convex, essentially smooth, takes its unique minimum of $f(t) = 0$ at $t=1$, and satisfies $f(t) = \infty$ for all $t < 0$. The interior of $\text{dom}(f) \equiv \{t \in \mathbb{R} \; : \; f(t) < \infty\}$, denoted $(\ell, u)$, satisfies $\ell < 1 < u$, and $f$ is twice continuously differentiable on $(\ell, u)$. \label{Assumption: setting, divergence}
\end{enumerate}
\end{restatable}
The finite support condition, assumption \ref{Assumption: setting} \ref{Assumption: setting, support of always observed variables}, 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{Section: duality} hold and can be used to define the estimator. This assumption can also be relaxed. In settings where assumption \ref{Assumption: setting} \ref{Assumption: setting, support of always observed variables} fails, researchers can still perform a breakdown point analysis through a conservative procedure described in remark \ref{Remark: always observed variables have finite support} below.
Condition \ref{Assumption: setting, divergence} ensures the $f$-divergence used to measure selection is well behaved, and is satisfied by every divergence in Table \ref{Table: Common f divergences}. In particular, strict convexity of $f$ ensures the primal problem \eqref{Definition: primal problem} 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.
\section{Duality}
\label{Section: duality}
As defined in display \eqref{Definition: primal problem}, $\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 \citep{csiszar1999mem}, partially finite programming \citep{borwein1991duality}, or $f$-divergence projection \citep{broniatowski2006minimization}. The convex analysis results in these papers connect the primal problem in display \eqref{Definition: primal problem} 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$:
\begin{align*}
\nu(b) = &\inf_q E\left[f(q(Y, X)) \mid D = 1\right] \\
&\text{s.t. } E[h(Y, X, b) q(Y,X) \mid D =1] = c(b)
\end{align*}
where
\begin{align}
&h(y, x, b) \equiv
\begin{pmatrix}
g(y, x, b) \\
\mathbbm{1}\{x = x_1\} \\
\vdots \\
\mathbbm{1}\{x = x_K\}
\end{pmatrix},
&&c(b) \equiv
\begin{pmatrix}
\frac{-p_D}{1-p_D} E[g(Y, X, b) \mid D = 1] \\
P(X = x_1 \mid D = 0) \\
\vdots \\
P(X = x_K \mid D = 0)
\end{pmatrix}, \label{Definition: constraint notation, h(z,b) and c(b)}
\end{align}
As shown in \cite{borwein1991duality}, the dual problem corresponding to \eqref{Definition: primal problem} is given by
\begin{equation}
V(b) \equiv \sup_{\lambda \in \mathbb{R}^{d_g + K}} \lambda^\intercal c(b) - E\left[f^*\left(\lambda^\intercal h(Y, X, b)\right) \mid D = 1\right] \label{Definition: dual problem}
\end{equation}
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{Table: Common f divergence conjugates} summarizes the convex conjugate for several common divergences.
\begin{table}[H]
\caption{Common $f$-divergence conjugates and effective domains}
\label{Table: Common f divergence conjugates}
\begin{center}
\begin{tabular}{|l || l | l | l | l |} \hline
Name & $f(t)$ & $\ell$, $u$ & $f^*(r)$ & $\ell^*$, $u^*$ \\ \hline \hline
Squared Hellinger & $\frac{1}{2}(\sqrt{t} - 1)^2$ & $\ell = 0$, $u = \infty$ & $\frac{1}{2}\left(\frac{1}{1-2r} - 1\right)$ & $\ell^* = -\infty$, $u^* = \frac{1}{2}$ \\ \hline
Kullback-Leibler (KL) & $t\log(t) - t + 1$ & $\ell = 0$, $u = \infty$ & $\exp(r) - 1$ & $\ell^* = -\infty$, $u^* = \infty$\\ \hline
``Reverse'' KL & $-\log(t) + t - 1$ & $\ell = 0$, $u = \infty$ & $-\log(1-r)$ & $\ell^* = -\infty$, $u^* = 1$ \\\hline
Cressie-Read & $f_\gamma(t) = \frac{t^\gamma - \gamma t + \gamma -1}{\gamma(\gamma - 1)}$ & -- & $\frac{1}{\gamma}(\gamma r - r + 1)^{\frac{\gamma}{\gamma - 1}} - \frac{1}{\gamma}$ & -- \\ \hline
\end{tabular}
\end{center}
\end{table}
\begin{remark}
\label{Remark: ensuring a probability density}
To ensure $q$ corresponds to a probability density, the constraints must enforce $\int q(z) dP(z) = 1$. This is implied by the constraints ensuring $Q_X = P_{0X}$ when $X$ is present. If $X$ is empty, set $h(z, b) = \begin{pmatrix} g(z, b)^\intercal & 1 \end{pmatrix}^\intercal \in \mathbb{R}^{d_g + 1}$ and $c(b) = \begin{pmatrix} \frac{-p_D}{(1-p_D)} E[g(Y, X, b) \mid D = 1]^\intercal & 1 \end{pmatrix}^\intercal \in \mathbb{R}^{d_g + 1}$ to ensure $q$ integrates to $1$.
\end{remark}
\subsection{Weak and strong duality}
\label{Subsection: weak and strong duality}
Assumption \ref{Assumption: setting} suffices to show $V(b) \leq \nu(b)$. This fact is known as \textit{weak duality}, and implies that
\begin{equation}
\inf_{b \in B \cap \textbf{B}_0} V(b) \leq \inf_{b \in B \cap \textbf{B}_0} \nu(b) = \delta^{BP}
\label{Display: weak duality}
\end{equation}
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 \textit{strong duality} holds, that is, $V(b) = \nu(b)$. Recall from Assumption \ref{Assumption: setting} \ref{Assumption: setting, divergence} that the interior of $\text{dom}(f) = \{t \in \mathbb{R} \; : \; f(t) < \infty\}$ is denoted $(\ell, u)$.
\begin{restatable}[Strong duality]{assumption}{assumptionStrongDuality}
\label{Assumption: strong duality}
\singlespacing
$B \subseteq \textbf{B}$ is convex, compact, and satisfies $\inf_{b \in \textbf{B}_0} \nu(b) = \inf_{b \in B \cap \textbf{B}_0} \nu(b)$. Furthermore, for each $b \in B$,
\begin{enumerate}[label=(\roman*)]
\item there exists $Q^b \in \textbf{P}^b$ such that $\ell < \frac{\partial Q^b}{\partial P_1}(z) < u$, almost surely $P_1$, and \label{Assumption: strong duality, constraint qualification}
\item $\lambda(b)$ solving \eqref{Definition: dual problem} is in the interior of $\{\lambda \; : \; E[\lvert f^*(\lambda^\intercal h(Z, b)) \rvert \mid D = 1] < \infty\}$. \label{Assumption: strong duality, interior dual solution}
\end{enumerate}
\end{restatable}
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 \cite{borwein1991duality}. The proof of theorem \ref{Theorem: strong duality}, found in the Supplementary Material, uses a result due to \cite{csiszar1999mem}.}
\begin{restatable}[Strong duality]{theorem}{theoremStrongDuality}
\label{Theorem: strong duality}
\singlespacing
Suppose assumptions \ref{Assumption: setting} and \ref{Assumption: strong duality} hold. Then for each $b \in B$, $\nu(b) = V(b)$, with dual attainment.
\end{restatable}
The first order condition of the dual problem \eqref{Definition: dual problem} provides intuition. Exchanging expectation and differentiation, the first order condition is
\begin{align*}
\begin{pmatrix}
\frac{-p_D}{1-p_D} E_{P_1}[g(Y, X, b)] \\
P(X = x_1 \mid D = 0) \\
\vdots \\
P(X = x_K \mid D = 0)
\end{pmatrix} =
E_{P_1}\left[(f^*)'\left(\lambda(b)^\intercal h(Y, X, b)\right)
\begin{pmatrix}
g(Y, X, b) \\
\mathbbm{1}\{X = x_1\} \\
\vdots \\
\mathbbm{1}\{X = x_K\}
\end{pmatrix}
\right]
\end{align*}
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{Theorem: strong duality} shows that under assumptions \ref{Assumption: setting} and \ref{Assumption: strong duality}, $(f^*)'\left(\lambda(b)^\intercal h(y, x, b)\right)$ is the $P_1$-density of the solution to the primal problem.
Assumption \ref{Assumption: strong duality} 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{Example: Squared Hellinger Expectation}. $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.
\begin{remark}
\label{Remark: always observed variables have finite support}
Assumption \ref{Assumption: setting} \ref{Assumption: setting, support of always observed variables} asks that $X$ be finitely supported, ensuring that $\textbf{P}^b$ is characterized by a finite number of moments and hence that the dual problem \eqref{Definition: dual problem} is finite dimensional. In settings where assumption \ref{Assumption: setting} \ref{Assumption: setting, support of always observed variables} fails and $X$ is not finitely valued, one can still conduct a conservative breakdown point analysis. Specifically, requiring $Q_X$ match a finite number of moments of $P_{0X}$ will estimate a value no larger than $\delta^{BP}$. If this value is large enough to assuage missing data concerns, the researcher is assured the breakdown point is weakly larger.
To illustrate, consider requiring that $Q_X$ match the first moment of $P_{0X}$. Define
\begin{align*}
&\tilde{h}(z, b) = \tilde{h}(y, x, b) =
\begin{pmatrix}
g(y, x, b) \\
x \\
1
\end{pmatrix},
&&\tilde{c}(b) =
\begin{pmatrix}
\frac{-p_D}{1-p_D} E[g(Y,X,b) \mid D = 1] \\
E[X \mid D = 0] \\
1
\end{pmatrix}
\end{align*}
and consider the value of the problem
\begin{equation*}
\tilde{V}(b) \equiv \sup_{\lambda \in \mathbb{R}^{d_g + d_x}} \lambda^\intercal \tilde{c}(b) - E\left[f^*\left(\lambda^\intercal \tilde{h}(Y, X, b)\right) \mid D = 1 \right].
\end{equation*}
This is the dual problem corresponding to the primal minimization problem $\inf_{Q \in \tilde{\textbf{P}}^b} d_f(Q \Vert P_1)$, where $\tilde{\textbf{P}}^b$ is the set of distributions that zero the moment conditions and match the first moment of $X$:
\begin{equation*}
\tilde{\textbf{P}}^b \equiv \left\{Q \; : \; Q \ll P_1, \; E_Q[X] = E_{P_{0X}}[X], \; p_D E_{P_1}[g(Y, X, b)] + (1-p_D) E_{Q}[g(Y, X, b)] = 0\right\}
\end{equation*}
The minimization problem $\inf_{Q \in \tilde{\textbf{P}}^b} d_f(Q \Vert P_1)$ is less constrained than the problem defining $\nu(b)$ in \eqref{Definition: primal problem}, and so has a lower value function. Minimizing $\tilde{V}(b)$ over $b \in B \cap \textbf{B}_0$ will therefore attain a lower value than $\delta^{BP}$. If $\inf_{b \in B \cap \textbf{B}_0} \tilde{V}(b)$ is large enough to assuage selection concerns, the reader is assured that $\delta^{BP}$ could only be larger. Moreover, the statistical properties of the estimator studied in section \ref{Section: estimators and asymptotics} are essentially unchanged when replacing $h$ and $c$ with $\tilde{h}$ and $\tilde{c}$ respectively.
As \cite{borwein1993failure} shows by counterexample, infinite dimensional analogues of \eqref{Definition: primal problem} and \eqref{Definition: dual problem} can fail to satisfy strong duality. \cite{borwein1993failure} also shows that duality can be restored by relaxing or penalizing the primal problem, and taking appropriate limits. This suggests another promising approach to relaxing assumption \ref{Assumption: setting} \ref{Assumption: setting, support of always observed variables}, but would considerably complicate estimation and so is left for future research.
\end{remark}
\section{Estimation}
\label{Section: estimators and asymptotics}
Assumptions \ref{Assumption: setting} and \ref{Assumption: strong duality} 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.
\subsection{The estimator}
\label{Section: estimators and asymptotics, 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) = \begin{bmatrix} -D I_{d_g} & 0 \\ 0 & (1-D) I_K \end{bmatrix}$ 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
\begin{equation}
\nu(b) = \sup_{\lambda \in \mathbb{R}^{d_g + K}} E\left[\frac{\lambda^\intercal J(D) h(DY, X, b)}{1-p_D} - \frac{D f^*\left(\lambda^\intercal h(DY, X, b)\right)}{p_D} \right]. \label{Definition: dual problem unconditional expectation}
\end{equation}
Define
\begin{equation}
\varphi(D, DY, X, b, \lambda, p) \equiv \frac{\lambda^\intercal J(D) h(DY, X, b)}{1-p} - \frac{D}{p} f^*(\lambda^\intercal h(DY, X, b)), \label{Definition: dual objective integrand}
\end{equation}
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
\begin{equation}
\hat{\nu}_n(b) \equiv \sup_{\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}), \label{Definition: estimator of value function}
\end{equation}
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.
\subsection{Asymptotic normality}
\label{Section: estimators and asymptotics, 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
\begin{equation}
\phi(D, DY, X, b, \theta) = \phi(D, DY, X, b, v, \lambda, p) =
\begin{pmatrix}
\varphi(D, DY, X, b, \lambda, p) - v \\
\nabla_\lambda \varphi(D, DY, X, b, \lambda, p) \\
D - p
\end{pmatrix}, \label{Definition: Z-estimator integrand}
\end{equation}
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\}$.
\begin{restatable}[Estimation]{assumption}{assumptionEstimation}
\label{Assumption: estimation}
\singlespacing
Suppose that
\begin{enumerate}[label=(\roman*)]
\item $\textbf{B}_0$ is closed, \label{Assumption: estimation, closed null hypothesis}
\item $\min_{b \in B \cap \textbf{B}_0} \nu(b)$ has a unique solution, \label{Assumption: estimation, unique solution}
\item the matrix $E[h(Y, X, b)h(Y, X, b)^\intercal \mid D = 1]$ is nonsingular for each $b \in B$, \label{Assumption: estimation, nonsingular second moments}
\item $g(y, x, b)$ is continuously differentiable with respect to $b$ for each $(y,x)$, and \label{Assumption: estimation, continuously differentiable moment functions}
\item there exists a convex, compact set $\Theta^B$ containing $\text{Gr}(\theta_0)^\eta$ for some $\eta > 0$ satisfying
\begin{align*}
&E\left[\sup_{(b, \theta) \in \Theta^B} \lVert \phi(D, DY, X, b, \theta) \rVert^2\right] < \infty &&\text{ and } &&E\left[\left(\sup_{(b,\theta) \in \Theta^B} \lVert \nabla_{(b,\theta)} \phi(D, DY, X, b, \theta) \rVert\right)^2\right] < \infty.
\end{align*}
\label{Assumption: estimation, moment conditions}
\end{enumerate}
\end{restatable}
As previewed in section \ref{Section: measuring selection and breakdown analysis, preview of results}, $\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{Assumption: estimation, nonsingular second moments}, \ref{Assumption: estimation, continuously differentiable moment functions}, and \ref{Assumption: estimation, moment conditions} 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{Assumption: estimation, closed null hypothesis} and \ref{Assumption: estimation, unique solution} 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{Assumption: estimation} \ref{Assumption: estimation, closed null hypothesis} and \ref{Assumption: estimation, continuously differentiable moment functions} are easily verified by inspection of $\textbf{B}_0$ and $g$ respectively. Conditions \ref{Assumption: estimation, nonsingular second moments} and \ref{Assumption: estimation, moment conditions} are similar to conditions required of generalized empirical likelihood estimators (see, e.g., \cite{antoine2021robust} assumption 1 (v) and assumption 3 (iv), (vii)). Assumption \ref{Assumption: estimation} \ref{Assumption: estimation, unique solution} deserves additional scrutiny. When $\textbf{B}_0$ is a convex set, condition \ref{Assumption: estimation, unique solution} 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.
\begin{restatable}[Convex value function, linear models]{lemma}{lemmaConvexDualValueFunctionLinearModels}
\label{Lemma: linear models imply a convex value function}
\singlespacing
Suppose assumptions \ref{Assumption: setting} and \ref{Assumption: strong duality} hold, the sample is $\{D_i, D_i Y_i, X_{i1}, X_{i2}\}_{i=1}^n$ where $Y_i \in \mathbb{R}$, $X_{i1} \in \mathbb{R}^{d_{x1}}$, and $X_{i2} \in \mathbb{R}^{d_{x2}}$, and the parameter $\beta$ is identified by
\begin{equation*}
E[(Y - X_1^\intercal \beta) X_2] = 0
\end{equation*}
Then $\hat{\nu}_n$ and $\nu$ are convex. If in addition $E[X_2X_1^\intercal]$ has full column rank, then $\nu$ is strictly convex.
\end{restatable}
\noindent Lemma \ref{Lemma: linear models imply a convex value function} 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{Lemma: linear models imply a convex value function} also produce convex $\nu$. Remark \ref{Remark: allowing for multiple minimizers} below discusses an approach to relaxing assumption \ref{Assumption: estimation} \ref{Assumption: estimation, unique solution}, at the cost of additional complexity.
Theorem \ref{Theorem: asymptotic normality} 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
\begin{equation}
\hat{\Phi}_n(b) \equiv \frac{1}{n}\sum_{i=1}^n \nabla_{\theta} \phi(D, DY, X, b, \hat{\theta}_n(b)), \label{Definition: estimator of asymptotic variance}
\end{equation}
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})$.
\begin{restatable}[Asymptotic normality]{theorem}{theoremAsymptoticNormality}
\label{Theorem: asymptotic normality}
\singlespacing
Suppose assumptions \ref{Assumption: setting}, \ref{Assumption: strong duality}, and \ref{Assumption: estimation} hold. Let $\hat{b}_n \equiv $ \\ $\operatorname*{arg\,min}_{b \in B \cap \textbf{B}_0} \hat{\nu}_n(b)$ and
\begin{equation*}
\hat{\sigma}_n^2 \equiv \frac{1}{n}\sum_{i=1}^n \left((\hat{\Phi}_n(\hat{b}_n)^{-1})^{(1)} \phi(D, DY, X, \hat{b}_n, \hat{\theta}_n(\hat{b}_n))\right)^2
\end{equation*}
where $(\hat{\Phi}_n(\hat{b}_n)^{-1})^{(1)} $ is the first row of the matrix $\hat{\Phi}_n(\hat{b}_n)^{-1}$. Then $\sqrt{n}(\hat{\delta}_n^{BP} - \delta^{BP})/\hat{\sigma}_n \overset{d}{\rightarrow} N(0,1)$.
\end{restatable}
\subsection{Inference}
\label{Section: estimators and asymptotics, 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{Theorem: asymptotic normality} implies that under assumptions \ref{Assumption: setting}, \ref{Assumption: strong duality}, and \ref{Assumption: estimation},
\begin{equation}
\widehat{CI}_{L,n} \equiv \hat{\delta}_n - \frac{\hat{\sigma}_n}{\sqrt{n}} c_{1-\alpha} \label{Definition: lower confidence interval}
\end{equation}
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.
\begin{comment}
$z_{1-\alpha}$ satisfies $P(Z \leq z_{1-\alpha}) = 1-\alpha$.
\begin{align*}
\lim_{n \rightarrow \infty} P\left(\hat{\delta}_n^{BP} - \frac{\hat{\sigma}_n}{\sqrt{n}} z_{1-\alpha} \leq \delta^{BP}\right) &= \lim_{n \rightarrow \infty} P\left(\hat{\delta}_n^{BP} - \delta^{BP} \leq \frac{\hat{\sigma}_n}{\sqrt{n}} z_{1-\alpha}\right) \\
&= \lim_{n \rightarrow \infty} P\left(\frac{\sqrt{n}(\hat{\delta}_n^{BP} - \delta^{BP})}{\hat{\sigma}_n} \leq z_{1-\alpha}\right) \\
&= P(Z \leq z_{1-\alpha})
\end{align*}
where the last equality follows from the definition of convergence in distribution and the fact that the standard normal CDF is continuous at all points.
\end{comment}
\begin{remark}
\label{Remark: allowing for multiple minimizers}
Assumption \ref{Assumption: estimation} \ref{Assumption: estimation, unique solution} can be relaxed at the cost of additional complexity. Without assumption \ref{Assumption: estimation} \ref{Assumption: estimation, unique solution}, $\sqrt{n}(\hat{\nu}_n - \nu)$ still converges in $\ell^\infty(B)$ to $\mathbb{G}_\nu$, a tight Gaussian process on $B$, and minimization of a function over $B \cap \textbf{B}_0$ remains a (Hadamard) \textit{directionally} differentiable map on the set of continuous functions of $B$. The delta method continues to imply $\sqrt{n}(\hat{\delta}_n^{BP} - \delta^{BP})$ converges in distribution to $\inf_{b \in \textbf{m}(\nu)} \mathbb{G}_\nu(b)$, where $\textbf{m}(\nu)$ is the set of minimizers of $\nu$.
Given a bootstrap $\hat{\nu}_n^*$ such that $\sqrt{n}(\hat{\nu}_n^* - \hat{\nu}_n)$ converges weakly in probability conditional on $\{D_i, D_i Y_i, X_i\}_{i=1}^n$ to $\mathbb{G}_\nu$, confidence intervals can still be constructed by utilizing the tools developed in \cite{fang2019inference}. One approach is to estimate the set $\textbf{m}(\nu)$ through ``near maximizers'' of $\hat{\nu}_n$ and use this estimated set to form an estimator of the map $h \mapsto \inf_{b \in \textbf{m}(\nu)} h(b)$. The confidence interval for $\delta^{BP}$ is formed by replacing $\hat{\sigma}_n c_{1-\alpha}$ in display \eqref{Definition: lower confidence interval} with the $1-\alpha$ quantile of this estimated function applied to the bootstrap sample; see \cite{fang2019inference} theorem 3.2 and appendix lemma S.4.8.. As most cases of interest appear to satisfy assumption \ref{Assumption: estimation} \ref{Assumption: estimation, unique solution}, this extension is left for future research.
\end{remark}
\begin{comment}
\subsection{Alternative for directional differentiability}
I'm not sure it makes sense to work out or code the estimator for the directionally differentiable case. I'll leave my reasoning commented out here in case I need to revisit it.
There are two possibilities:
\begin{enumerate}
\item There's a unique minimizer. Then we have asymptotic normality and can simply estimate the variance.
\begin{itemize}
\item This feels like the most common case anyway.
\end{itemize}
\item There is not a unique minimizer. Two more sub-cases:
\begin{enumerate}[label=(\roman*)]
\item $\nu$ is not convex
\begin{itemize}
\item \cite{fang2019inference} lemma S.4.8 suggests that a valid confidence interval could be constructed as follows: find a valid bootstrap $\sqrt{n}(\hat{\nu}_{n,b}^* - \hat{\nu}_n)$ estimating $\mathbb{G}_\nu$. Choose $\kappa_n \uparrow \infty$ with $\frac{\kappa_n}{\sqrt{n}} \downarrow 0$, and define
\begin{align*}
&\hat{m}_n(\nu) \equiv \left\{b \in B \cap \textbf{B}_0, \; ; \; \hat{\nu}_n(b) \leq \inf_{\tilde{b} \in B \cap \textbf{B}_0} \hat{\nu}_n(\tilde{b}) + \kappa_n\right\}, &&\hat{\iota}_n'(h) \equiv \inf_{b \in \hat{m}_n(\nu)} h(b)
\end{align*}
Use the $1-\alpha$ quantile of $\{\hat{\iota}_n(\sqrt{n}(\hat{\nu}_{n,b}^* - \hat{\nu}_n))\}_{b=1}^B$ in place of $z_{1-\alpha}$ in equation \eqref{Definition: lower confidence interval}.
\item However, with $\nu$ not being convex, $\hat{m}_n(\nu)$ will generally be a non-convex set. $\sqrt{n}(\hat{\nu}_n^* - \hat{\nu}_n)$ will also generally be a non-convex function. So this procedure involves solving $B$ minimization problems (usually \textit{thousands}) which each have non-convex objectives, over non-convex feasible sets... its not computationally feasible.
\end{itemize}
\item $\nu$ is convex
\begin{itemize}
\item Are you sure there isn't a unique minimizer? With $\textbf{B}_0$ convex, all we need is strict convexity of $\nu$ to gaurantee a unique minimizer...
\item If $\nu(b)$ really does ``flatten out'' at its minimum, \cite{fang2019inference} may be helpful and computationally feasible here. But do we have an example falling into this case?
\end{itemize}
\end{enumerate}
\end{enumerate}
\end{comment}
\section{Simulations}
\label{Section: 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.
\subsection{Expectation}
\label{Section: simulations, simple mean}
Recall example \ref{Example: Squared Hellinger Expectation}. 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}$.}
\begin{table}[H]
\caption{Simulations, expectation}
\begin{center}
\begin{tabular}{|c || c|c|c|c |} \hline
$n$ & Bias & St. Dev. & Coverage & Ave. CI Length \\ \hline \hline
1,000 & 0.005 & 0.056 & 98.5 & 0.090 \\ \hline
3,000 & 0.002 & 0.032 & 96.3 & 0.051 \\ \hline
5,000 & 0.001 & 0.025 & 95.8 & 0.039 \\ \hline
10,000 & 0.001 & 0.017 & 95.8 & 0.028 \\ \hline
\end{tabular}
\end{center}
\end{table}
\noindent The simulations show little bias. Coverage is slightly above the targeted 95 percent significance level in smaller samples.
\subsection{Linear model}
\label{Section: simulations, linear models}
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
\begin{equation}
Y_1 = \beta_0 + \beta_1 X_1 + \beta_2 Y_2 + \beta_3 X_2 + \varepsilon = W^\intercal \beta + \varepsilon, \label{Simulation: MLR}
\end{equation}
where $W = \begin{pmatrix} 1 & X_1 & Y_2 & X_2 \end{pmatrix}^\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:
\begin{align}
&H_0 \; : \; \beta_1 \leq 0, &&H_1 \; : \; \beta_1 > 0 \label{Display: linear simulations breakdown point conclusion}
\end{align}
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 \eqref{Simulation: MLR}. Notice the support of $(Y_1, Y_2, X_1, X_2)$ is compact, ensuring the moment conditions in assumption \ref{Assumption: estimation} \ref{Assumption: estimation, moment conditions} are satisfied. The Supplementary Material shows simulation evidence that this model and data generating process produces a convex $\nu(\cdot)$, suggesting that assumption \ref{Assumption: estimation} \ref{Assumption: estimation, unique solution} 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:
\begin{table}[H]
\caption{Simulations, linear model}
\begin{center}
\begin{tabular}{|c || c|c|c|c |} \hline
$n$ & Bias & St. Dev. & Coverage & Ave. CI Length \\ \hline \hline
1,000 & 0.016 & 0.048 & 98.9 & 0.076 \\ \hline
3,000 & 0.007 & 0.026 & 95.8 & 0.041 \\ \hline
5,000 & 0.004 & 0.019 & 95.4 & 0.031 \\ \hline
10,000 & 0.003 & 0.013 & 94.5 & 0.022 \\ \hline
\end{tabular}
\end{center}
\end{table}
\noindent The simulations again show little bias, with coverage slightly above the targeted 95 percent significance level in smaller samples.
\subsection{Logit model}
\label{Section: simulations, logistic regression}
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
\begin{equation*}
E[(Z_1 - \Lambda(Z_{-1}^\intercal \beta)) Z_{-1}] = 0.
\end{equation*}
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
\begin{align}
&H_0 : \Lambda(\bar{z}^\intercal \beta) \leq 0.5, &&H_1 : \Lambda(\bar{z}^\intercal \beta) > 0.5. \label{Display: logistic simulations breakdown point conclusion}
\end{align}
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{Assumption: estimation} \ref{Assumption: estimation, moment conditions} 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{Assumption: estimation} \ref{Assumption: estimation, unique solution} 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 \eqref{Display: logistic simulations breakdown point conclusion} 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:
\begin{table}[H]
\caption{Simulations, logit model}
\begin{center}
\begin{tabular}{| c || c|c|c|c |} \hline
$n$ & Bias & St. Dev. & Coverage & Ave. CI Length \\ \hline \hline
1,000 & 0.003 & 0.018 & 94.5 & 0.029 \\ \hline
3,000 & -0.000 & 0.010 & 96.1 & 0.017 \\ \hline
5,000 & 0.001 & 0.008 & 94.8 & 0.013 \\ \hline
10,000 & -0.000 & 0.005 & 95.9 & 0.009 \\ \hline
\end{tabular}
\end{center}
\end{table}
\noindent These simulations show essentially zero bias and correct coverage at relatively small sample sizes.
\section{Application: attrition in randomized controlled trials}
\label{Section: application}
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{Section: measuring selection and breakdown analysis, squared hellinger}: 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. \cite{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 \cite{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 \cite{barham2024experimental}, \cite{bandiera2020women}, and \cite{giacobino2024schoolgirls}, which all study randomized controlled trials conducted in developing countries. \cite{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. \cite{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, \cite{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{Table: Application, MCAR estimates} 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.
\begin{table}[H]
\caption{Missing completely at random estimates}
\label{Table: Application, MCAR estimates}
\begin{center}
\begin{tabular}{|l | c || c | c | c | c | c |}
\hline
\multirow{3}{*}{Paper} & \multirow{3}{*}{Outcome} & \multicolumn{2}{c|}{Observations} & Untreated & \multicolumn{2}{c|}{MCAR ITT estimates} \\
& & Total $(n)$ & Complete & mean & Replication & Region ind. \\
& & (1) & (2) & (3) & (4) & (5) \\
\hline \hline
\multirow[c]{9}{*}{\begin{tabular}{l} Barham \\ et al. (2024)\end{tabular}} & & & & & & \\
& Off-farm & 1,138 & 1,006 & 0.83 & 0.06 & 0.06 \\
& employment & & & & (0.02) & (0.03) \\
& & & & & & \\
& Rank of earnings & 1,138 & 1,006 & 497.15 & 41.78 & 42.87 \\
& per mo. worked & & & & (19.50) & (21.58) \\
& & & & & & \\
& Read and write & 1,138 & 1,007 & 0.87 & 0.05 & 0.07 \\
& & & & & (0.02) & (0.02) \\
\cline{1-7}
\multirow[c]{15}{*}{\begin{tabular}{l} Bandiera \\ et al. (2020)\end{tabular}} & & & & & & \\
& Entrepreneurial & 5,966 & 4,765 & 71.77 & 5.63 & 5.46 \\
& ability index & & & & (0.98) & (0.71) \\
& & & & & & \\
& Income-generating & 5,966 & 4,831 & 0.10 & 0.07 & 0.07 \\
& activity & & & & (0.02) & (0.01) \\
& & & & & & \\
& Self-employed & 5,966 & 4,831 & 0.06 & 0.06 & 0.06 \\
& & & & & (0.01) & (0.01) \\
& & & & & & \\
& Wage employed & 5,966 & 4,831 & 0.04 & 0.01 & 0.01 \\
& & & & & (0.01) & (0.01) \\
& & & & & & \\
& Expenditures & 5,966 & 4,752 & 11.92 & 4.68 & 4.52 \\
& in UGX, 1,000s & & & & (0.95) & (0.73) \\
\cline{1-7}
\multirow[c]{9}{*}{\begin{tabular}{l} Giacobino \\ et al. (2024) \end{tabular}} & & & & & & \\
& Dropped out & 1,501 & 1,344 & 0.40 & -0.21 & -0.21 \\
& & & & & (0.05) & (0.02) \\
& & & & & & \\
& Married & 1,501 & 1,344 & 0.14 & -0.07 & -0.07 \\
& & & & & (0.03) & (0.02) \\
& & & & & & \\
& Life satisfaction & 1,501 & 1,344 & 0.00 & 0.25 & 0.25 \\
& & & & & (0.11) & (0.05) \\
\cline{1-7}
\end{tabular}
\end{center}
\begin{flushleft}
\vspace{-0.3 cm}
{\scriptsize \textit{Notes:} Column (4) replicates each paper's results, including the authors' choice of additional regressors and standard errors. Column (5) uses region indicators as the additional regressors and reports HC3 standard errors in parentheses. For outcomes from Barham et al. (2024), column (3) reports the mean from subjects receiving the CCT late and estimates in column (5) make use of sampling weights from Molina-Mill\'an \& Macours (2025) rather than inverse propensity weights from Barham et al. (2024). For Bandiera et al. (2020), column (3) reports average outcomes measured at baseline.}
\end{flushleft}
\end{table}
Table \ref{Table: Application, BDP estimates} 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{Table: Application, MCAR estimates}. 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{Remark: squared hellinger lower bound}.
\begin{table}[H]
\caption{Breakdown point analyses}
\label{Table: Application, BDP estimates}
\begin{center}
\begin{tabular}{|l | c || c | c | c | c |}
\hline
\multirow{2}{*}{Paper} & \multirow{2}{*}{Outcome} & \multirow{2}{*}{$\hat{p}_{D,n}$} & \multirow{2}{*}{$\hat{\delta}_n^{BP}$} & \multirow{2}{*}{$\widehat{CI}_{L,n}$} & Estimate of \\
& & & & & $H^2(P_{0X}, P_{1X})$ \\
\hline \hline
\multirow[c]{3}{*}{Barham et al. (2024)} & Off-farm employment & 0.88 & 0.08 & 0.03 & 0.05 \\
& Rank of earnings per mo. worked & 0.88 & 0.08 & 0.03 & 0.05 \\
& Read and write & 0.88 & 0.07 & 0.04 & 0.05 \\
\hline
\multirow[c]{5}{*}{Bandiera et al. (2020)} & Entrepreneurial ability index & 0.81 & 0.08 & 0.06 & 0.03 \\
& Income-generating activity & 0.82 & 0.05 & 0.04 & 0.03 \\
& Self-employed & 0.82 & 0.04 & 0.03 & 0.03 \\
& Wage employed & 0.82 & 0.03 & 0.03 & 0.03 \\
& Expenditures in UGX, 1,000s & 0.81 & 0.04 & 0.04 & 0.03 \\
\hline
\multirow[c]{3}{*}{Giacobino et al. (2024)} & Dropped out & 0.90 & $\infty$ & $\infty$ & 0.05 \\
& Married & 0.90 & 0.17 & 0.06 & 0.05 \\
& Life satisfaction & 0.90 & 0.26 & 0.09 & 0.05 \\
\hline
\end{tabular}
\end{center}
\end{table}
Tables \ref{Table: Application, MCAR estimates} and \ref{Table: Application, BDP estimates} show a number of patterns worth emphasizing. Compared to the other two studies, the larger sample of \cite{bandiera2020women} resulted in lower standard errors in table \ref{Table: Application, MCAR estimates} and lower confidence intervals that are closer to the breakdown point estimates in table \ref{Table: Application, BDP estimates}. 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 \cite{barham2024experimental} differentially raised the probability of working somewhere other than the recipient's family farm, and the claim that the scholarship studied in \cite{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 \cite{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 \cite{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 \cite{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 \cite{giacobino2024schoolgirls} have point estimates for the breakdown point that are close to 0.2. As discussed in remark \ref{Remark: 0.2 is a big number}, these are large values for a breakdown point implying the results are quite robust.
\section{Conclusion}
\label{Section: 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}
\bibliography{./MissingDataBDP_bibliography.bib}
\clearpage
\singlespacing