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.
29,778 characters
Testing selection on observables in parametric models with refreshment samples
\maketitle
\onehalfspacing
\begin{abstract}
\linespread{1.2}
In panels with sample selection (that may occur due to attrition, nonresponse, etc.), the assumption of selection on observables (missing at random, MAR) is commonly imposed despite often being implausible.
However, this assumption becomes testable when a refreshment sample is available.
We develop a statistical test of MAR based on a distance between two estimated distributions: one obtained using the standard inverse probability weighting (IPW) that is valid under MAR and the other obtained using an alternative weighting that is valid under a weaker assumption of additive nonignorability of \citet{hirano2001combining}.
This test implicitly compares the distribution of the IPW-weighted sample in the attrition period with the distribution of the refreshment sample, which coincide if the MAR assumption holds.
We establish that, when the input distributions are parametric, our test statistic converges to the generalized chi-squared distribution under the null of MAR.
This limit distribution can be estimated using the recursive formulas derived by \citet{franguridi2025raking}.
We illustrate the performance of our test in Monte Carlo simulations.
Finally, we apply our test to an empirical example using a subsample of the Understanding America Study (UAS) dataset.
\medskip
\noindent \textbf{JEL Classification:} C23
\medskip
\noindent \textbf{Keywords:} refreshment sample, attrition, sample selection, missing at random, additive nonignorability
\end{abstract}
\newpage
\section{Introduction}
In the last 60 years or so, both empirical applications of panel data models and associated methodological developments have seen tremendous growth \citep[e.g.,][]{sarafidis2021celebrating}. One of the most intricate problems in using panel data is attrition, i.e., units dropping out of the sample in a nonrandom way.
This reduces the effective sample size and can bias the resulting estimates.
These biases can be completely removed when units are \textit{missing at random} (MAR; also called \textit{selection on observables}), i.e., when the probability of attrition depends only on observed variables. In that case, inverse probability weighting of the sample of stayers leads to unbiased estimates.
However, the MAR assumption can be restrictive, and panel data alone is not sufficient to determine its plausibility.
To compensate for attrition, researchers often collect so-called \emph{refreshment samples}, i.e., samples that are designed to ``replace'' the missing units \citep{kish1959variances}.
Refreshment samples are used in many empirical settings, including survey panels \citep{deng2013handling,si2015semi,franguridi2024closed,franguridi2025inference}, public transportation data \citep{ridder1992empirical,zheng2025semiparametric}, and retail scanner data \citep{chen2017retail}.
With a refreshment sample, the MAR assumption underlying the inverse probability weighting becomes testable.
We develop a statistical test of MAR that relies on a modeling framework with \textit{additively nonignorable} (AN) attrition introduced by \citet{hirano2001combining} and further used by \citet{nevo2003using,bhattacharya2008inference,hoonhout2019nonignorable,franguridi2025raking,franguridi2025robust} and many others.
The AN assumption allows the attrition to depend on yet unobserved variables and includes MAR as a special case.
We propose parameterizing the three directly estimable distributions (the first-period distribution, the balanced panel distribution, and the refreshment sample distribution) and compute the Hellinger distance between the joint distribution of data estimated under the MAR assumption (the inverse probability weighted distribution) and that estimated under the AN assumption as a test statistic.
To compute the AN distribution, we use the raking algorithm in \citet{franguridi2025raking}.
We show that our statistic has an asymptotically generalized chi-squared distribution, which can be estimated using the recursive formulas also derived in \citet{franguridi2025raking}.
The remainder of the paper is organized as follows.
\cref{sec:framework} introduces the framework.
\cref{sec:testing} describes the test statistic and derives its asymptotic distribution.
\cref{sec:mc-sim} illustrates the size and power of our test in a Monte Carlo simulation.
\cref{sec:empirical} uses the test to assess the MAR assumption in a subset of the Understanding America Study (UAS) panel survey.
\cref{sec:conclusion} concludes.
Finally, the Appendix specializes our results to the cases of Gaussian data and discrete data and describes the construction of the cognition variable.
\section{Framework}\label{sec:framework}
Consider a two-period panel in which units may drop out of the sample in period 2.
Let $Z_{it}=(Y_{it},X_{it})$ be stacked outcomes and covariates for unit $i$ in period $t=1,2$.
In period 1, data $Z_{i1}$ for all the units are observed, while in period 2, data $Z_{i2}$ are observed only for units who stay in the sample, which we denote by $S_i=1$.
Using only data of the selected sample $\{i:\,S_i=1\}$ may lead to significant bias. To reduce this bias, researchers often assume that the selection probability
\begin{align*}
p(z_1,z_2)=\operatorname{\mathbb{P}}(S_i=1|Z_{1i}=z_1,Z_{2i}=z_2)
\end{align*}
only depends on the observables $z_1$, not the potentially missing $z_2$.
This allows estimating this probability by running a (possibly nonparametric) regression of $S_i$ on $Z_{1i}$ and then using this estimated probability as the inverse weight for the selected sample.
The ``selection on observables'' (MAR) assumption may be implausible in many empirical settings.
Consider, for example, a panel survey.
It is likely that survey respondents may decide whether to continue participating in the survey depending on data yet unobserved to the researcher.
For example, if $Z_{it}$ is respondent $i$'s income at time $t$, then the respondent $i$ may drop out of the sample after experiencing a negative income shock, i.e., when $Z_{i2}$ is significantly smaller than $Z_{i1}$ (before reporting $Z_{i2}$ to the researcher).
Without auxiliary information, the MAR assumption is untestable.
However, in survey panels, researchers often collect so-called \emph{refreshment samples}, i.e., samples that are designed to ``replace'' the missing units.
We denote a refreshment sample in period 2 by $Z_{i2}^r$, $i=1,\dots,n_r$.
What makes it a ``refreshment'' sample is being an independent sample from the target period-2 marginal distribution, i.e., $Z_{i2}^r \overset{d}{=} Z_{j2}$.
Our key observation is that, with a refreshment sample, the MAR assumption becomes testable.
This is because the selection probability $p(z_1,z_2)$ is identifiable with a refreshment sample under the assumption of \emph{additive nonignorability} (AN) introduced by \citet{hirano2001combining}:
\begin{align}
p(z_1,z_2) = \exp(k_1(z_1)+k_2(z_2)) \text{ for some unknown functions } k_1,k_2. \label{eq:AN}
\end{align}
This assumption is significantly weaker than MAR and can be used to develop a computationally feasible test of the null hypothesis of MAR against the alternative hypothesis of AN.\footnote{We focus on the exponential link function $G(x)=\exp(x)$ because it leads to a computationally convenient procedure based on raking, as shown by \citet{franguridi2025raking}. In principle, all the theoretical results of this paper can be generalized to any smooth link function $G$.}
Under AN, \citet{hirano2001combining} showed that the joint density $f_\text{AN}(z_1,z_2)$ of $(Z_1,Z_2)$ is identified and can be recovered from the following three identified objects: the selected density $f^s(z_1,z_2) = f(z_1,z_2|S=1)$, the density $f_1(z_1)$ of $Z_1$ (identified with the first-period data where there is no attrition), and the density $f_2(z_2)$ of $Z_2$ (identified with the refreshment sample).
Specifically, $f_\text{AN}$ is the solution to the functional projection problem
\begin{align}
\min_{\tilde f\in\Pi(f_1,f_2)} \operatorname{KL}(\tilde f,f^s), \label{eq:KL-projection}
\end{align}
where $\Pi(f_1,f_2)$ is the set of joint densities with (multivariate) marginals $f_1$ and $f_2$ and
\begin{align*}
\operatorname{KL}(f,f^s) = \operatorname{\mathbb{E}}_{f} \left[ \log \frac{f(Z_1,Z_2)}{f^s(Z_1,Z_2)} \right]
\end{align*}
is the Kullback-Leibler (KL) divergence between $f$ and $f^s$.\footnote{This result also holds when the data contain a mix of discrete and continuous variables, in which case the densities involved should be understood as densities with respect to an appropriate product of Lebesgue and counting measures. For simplicity, we formulate all our results in the absolutely continuous case.}
In other words, $f_\text{AN}$ is the closest (in the KL sense) density to $f^s$ among densities with marginals $f_1$ and $f_2$.
Notice how the functions $k_1,k_2$ do not appear explicitly in this characterization of $f_\text{AN}$.
They do appear in the dual formulation of this problem, but this formulation is not required for this paper.
\section{Testing procedure}\label{sec:testing}
\subsection{Test statistic}
It is natural to seek evidence against MAR in a statistical distance between distributions $f_\text{AN}$ and $f_\text{MAR}$.
We consider the squared Hellinger distance due to its analytical tractability in our framework,\footnote{In principle, we could consider a test statistic based on any other statistical distance that depends smoothly on the densities of the input distributions, such as the Kullback-Leibler divergence, but the resulting formulas for the asymptotic distribution of such a statistic would be more complicated.}
\begin{align*}
\hat T_n=\frac{1}{2} \int \left( \sqrt{\hat f_\text{AN}(z_1,z_2)}-\sqrt{\hat f_\text{MAR}(z_1,z_2)}\right)^2 \, d z_1 dz_2,
\end{align*}
where $\hat f_\text{AN}$ and $\hat f_\text{MAR}$ are estimates of $f_\text{AN}$ and $f_\text{MAR}$, respectively.
Estimation of $f_\text{MAR}$ is easy. Indeed, the Bayes formula and the MAR assumption imply
\begin{align*}
f_\text{MAR}(z_1,z_2) = \frac{\operatorname{\mathbb{P}}(S=1) f^s(z_1,z_2)}{\operatorname{\mathbb{P}}(S=1|Z_1=z_1,Z_2=z_2)} = \frac{\operatorname{\mathbb{P}}(S=1) f^s(z_1,z_2)}{\operatorname{\mathbb{P}}(S=1|Z_1=z_1)} .
\end{align*}
Integrating with respect to $z_2$ and taking into account that the first marginal of $f_\text{MAR}$ is $f_1$, we obtain
\begin{align*}
f_\text{MAR}(z_1,z_2) = \frac{f_1(z_1) f^s(z_1,z_2)}{\int f^s(z_1,z_2)\, dz_2}.
\end{align*}
Hence, we can use the plug-in estimator
\[
\hat f_\text{MAR}(z_1,z_2) = \frac{\hat f_1(z_1) \hat f^s(z_1,z_2)}{\int \hat f^s(z_1,z_2)\, dz_2},
\]
where $\hat f_1, \hat f_2, \hat f^s$ are chosen estimators of $f_1,f_2,f^s$, respectively.
Notice that this formula forces $\hat f_\text{MAR}$ to have the first marginal $\hat f_1$.
Incidentally, $\hat f_\text{MAR}$ can be obtained from $\hat f^s$ after just one iteration of the raking algorithm described in \Cref{sec:raking}.
Estimating $f_\text{AN}$ is a more difficult problem, and we defer its exposition to \Cref{sec:raking}.
\begin{remark}
While the MAR assumption corresponds to the lack of dependence of the selection probability on $Z_2$, we may want to test its symmetric analog $\operatorname{\mathbb{P}}(S=1|Z_1, Z_2)=\operatorname{\mathbb{P}}(S=1|Z_2)$ instead. This can be handled easily because the target distribution in this case can be written as
\begin{align*}
f(z_1,z_2)= \frac{f_2(z_2) f^s(z_1,z_2)}{\int f^s(z_1,z_2) \, dz_1},
\end{align*}
where $f_2$ is identified using the refreshment sample. Besides, identification under this assumption is possible even without refreshment samples under additional parametric restrictions; see \citet{hausman1979attrition}.
\end{remark}
\subsection{Asymptotic distribution}
We now derive the asymptotic distribution of our test statistic under the assumption that the densities $f_1,f_2,f^s$ are parametric, with the parameters estimable at the standard $\sqrt{n}$ rate. Hence, we impose the following assumption.
\begin{assumption}\label{a:gamma}
\begin{subassumption}
\item The densities $f_1,f_2,$ and $f^s$ belong to parametric families $f_1(\cdot,\gamma_1)$, $f_2(\cdot,\gamma_2)$, and $f^s(\cdot,\gamma^s)$, where $\gamma = (\gamma_1',\gamma_2',(\gamma^s)')'$ is a finite-dimensional parameter.
\item The densities $f_1(\cdot,\gamma_1)$, $f_2(\cdot,\gamma_2)$, and $f^s(\cdot,\gamma^s)$ are twice continuously differentiable with respect to $\gamma$ and are positive on their respective supports at the true value $\gamma_0$ of $\gamma$.
\item For an estimator $\hat\gamma$ of $\gamma$, $\sqrt{n}(\hat\gamma - \gamma_0) \rightsquigarrow N(0,\Omega(\gamma_0))$, where $\Omega(\gamma_0)$ is a positive semidefinite matrix.
\end{subassumption}
\end{assumption}
Here $n=n_1+n_r$ is the combined sample size of the first period and the refreshment sample. Normally, $n_1$ and $n_r$ would be of the same order of magnitude, but this is not required for our results. Our main theorem is as follows.
\begin{theorem}
Suppose \Cref{a:gamma} holds and let $g_0 \sim N(0,\Omega(\gamma_0))$. Then, under $\gamma_0$ satisfying the MAR restriction,
\begin{align*}
n \, \hat T_{n} \rightsquigarrow \frac12 g_0' H(\gamma_0) g_0,
\end{align*}
where
\begin{align*}
H(\gamma_0) = \frac14 \operatorname{\mathbb{E}}_{f_{\gamma_0}} \left( s_{\text{AN}}(Z,\gamma_0) - s_{\text{MAR}}(Z,\gamma_0) \right) \left( s_{\text{AN}}(Z,\gamma_0) - s_{\text{MAR}}(Z,\gamma_0) \right)'
\end{align*}
is the Hessian of the squared Hellinger distance and
\begin{align*}
s_{\text{AN}}(z,\gamma) &= \nabla_\gamma \log f_\text{AN}(z,\gamma), \\
s_{\text{MAR}}(z,\gamma) &= \nabla_\gamma \log f_\text{MAR}(z,\gamma).
\end{align*}
are the score functions of $f_\text{AN}$ and $f_\text{MAR}$, respectively.
\end{theorem}
\begin{proof}
Define a function $T$ by
\begin{align*}
\hat T_n = T(\hat\gamma) = \frac{1}{2} \int \left( \sqrt{\hat f_\text{AN}(z_1,z_2,\hat\gamma)}-\sqrt{\hat f_\text{MAR}(z_1,z_2,\hat\gamma)}\right)^2 \, d z_1 dz_2,
\end{align*}
For brevity, denote $f_\gamma(\cdot) = f_\text{AN}(\cdot,\gamma)$, $g_\gamma(\cdot)=f_\text{MAR}(\cdot,\gamma)$, and $z=(z_1',z_2')'$.
Suppose $\gamma_0$ is such that $f_{\gamma_0}=g_{\gamma_0}$. We have
\begin{align*}
\nabla_\gamma T(\gamma_0)= \int \Big(\sqrt{f_{\gamma_0}(z)}-\sqrt{g_{\gamma_0}(z)}\Big) \Big(\nabla_\gamma \sqrt{f_{\gamma_0}(z)}-\nabla_\gamma \sqrt{g_{\gamma_0}(z)}\Big)\,dz=0.
\end{align*}
Now write $h_\gamma(z)=\sqrt{f_\gamma(z)}-\sqrt{g_\gamma(z)}$ so that
\begin{align*}
T(\gamma)=\frac12\int h_\gamma(z)^2\,dz.
\end{align*}
The $(i,j)$ entry of the Hessian of $T$ at $\gamma_0$ is
\begin{align*}
\partial_{ij} T(\gamma_0) &= \int \partial_i h_{\gamma_0}(z)\,\partial_j h_{\gamma_0}(z)\,dz + \int h_{\gamma_0}(z)\,\partial_{ij} h_{\gamma_0}(z)\,dz \\
&= \int \partial_i h_{\gamma_0}(z)\,\partial_j h_{\gamma_0}(z)\,dz \\
&= \frac14 \int
\left(\frac{\partial_i f_{\gamma_0}(z)}{\sqrt{f_{\gamma_0}(z)}}-\frac{\partial_i g_{\gamma_0}(z)}{\sqrt{g_{\gamma_0}(z)}}\right)
\left(\frac{\partial_j f_{\gamma_0}(z)}{\sqrt{f_{\gamma_0}(z)}}-\frac{\partial_j g_{\gamma_0}(z)}{\sqrt{g_{\gamma_0}(z)}}\right)\,dz \\
&= \frac14 \int
\left(\partial_i \log f_{\gamma_0}(z) - \partial_i \log g_{\gamma_0}(z)\right)
\left(\partial_j \log f_{\gamma_0}(z) - \partial_j \log g_{\gamma_0}(z)\right)
f_{\gamma_0}(z)\,dz.
\end{align*}
Applying the second-order delta method to the function $\gamma \mapsto T(\gamma)$ completes the proof.
\end{proof}
One convenient feature of the asymptotic distribution of our test statistic is that it does not depend on the second derivatives of $f_\text{AN}$ and $f_\text{MAR}$ with respect to $\gamma$, which would be the case for other distances such as the KL divergence.
\subsection{Estimator of the AN density}\label{sec:raking}
Here we describe a raking-based estimator of $f_\text{AN}$ introduced in \citet{franguridi2025raking}.
To this end, we introduce some notation.
First, let $\operatorname{\mathbb{I}}_1$ and $\operatorname{\mathbb{I}}_2$ be integration operators defined for any integrable function $f=f(z_1,z_2)$ by
\begin{align*}
[\operatorname{\mathbb{I}}_1 f](z_2) = \int f(z_1,z_2)\, dz_1, \qquad [\operatorname{\mathbb{I}}_2 f](z_1) = \int f(z_1,z_2)\, dz_2.
\end{align*}
Second, let $\Pi_1$ and $\Pi_2$ be the KL projection operators on the sets of distributions with the first marginal $f_1$ and the second marginal $f_2$, respectively,\footnote{See, e.g., \citet{franguridi2025raking} for a proof that these operators solve the KL projection problem.} i.e., for any density $f$,
\begin{align*}
[\Pi_1 f](z_1,z_2) = \frac{f_1(z_1) f(z_1,z_2)}{[\operatorname{\mathbb{I}}_2 f](z_1)}, \\
[\Pi_2 f](z_1,z_2) = \frac{f_2(z_2) f(z_1,z_2)}{[\operatorname{\mathbb{I}}_1 f](z_2)}.
\end{align*}
Third, define the composition operator $\Pi = \Pi_2 \circ \Pi_1$, i.e.,
\begin{align}
[\Pi f](z_1,z_2)= m & \left[ f_2(z_2,\gamma_2), \right. \notag \\
&m \left(f_1(z_1,\gamma_1), f(z_1,z_2), \operatorname{\mathbb{I}}_2 f(z_1,\cdot)\right), \notag \\
&\left. \operatorname{\mathbb{I}}_1 m\left( f_1(\cdot,\gamma_1), f(\cdot,z_2), \operatorname{\mathbb{I}}_2 f(\cdot,\cdot) \right) \right],
\end{align}
where $m(a,b,c)=ab/c$.
The raking procedure (also called iterative proportional fitting or Sinkhorn's algorithm) is an iteration of the operator $\Pi$ starting from a given density function.
It turns out that under weak conditions, the raking procedure initialized at $f^s$ converges to the solution of the KL projection problem \eqref{eq:KL-projection}, i.e.,
\begin{align*}
\lim_{T \to \infty} \Pi^{(T)} f^s = f_\text{AN},
\end{align*}
where $\Pi^{(T)}$ is the T-fold iteration of $\Pi$ and the convergence is understood in the $L^1$ sense, see \citet{ruschendorf1995convergence,franguridi2025raking}.
In other words, raking starts with $f^s$ and alternates between projection onto the set of distributions with the first marginal $f_1$ and projection onto the set of distributions with the second marginal $f_2$.
Dropping the variables $z_1,z_2$ and reflecting the dependence of the previous iteration on $\gamma$, we can write the $(t+1)$-th raking iteration as
\begin{align}
f^{(t+1)}(\gamma) = m & \left[ f_2(\gamma_2), \right. \\
&m \left(f_1(\gamma_1), f^{(t)}(\gamma), \operatorname{\mathbb{I}}_2 f^{(t)}(\gamma) \right), \\
&\left. \operatorname{\mathbb{I}}_1 m\left( f_1(\gamma_1), f^{(t)}(\gamma), \operatorname{\mathbb{I}}_2 f^{(t)}(\gamma) \right) \right], \label{eq:raking-recursion}
\end{align}
with the initial condition
\begin{align*}
f^{(0)}(\gamma) = f^s(\gamma).
\end{align*}
The raking estimator of \citet{franguridi2025raking} is then a sample analog of this procedure,
\begin{align*}
\hat f_\text{AN} = \hat\Pi^{(T)} \hat f^s,
\end{align*}
where $\hat\Pi$ is the sample raking operator that uses estimators $\hat f_1,\hat f_2$ in place of $f_1,f_2$; $\hat f^s$ is an estimator of $f^s$; and $T$ is the number of iterations chosen in advance.\footnote{The speed of convergence of the raking iterations was studied by \citet{ruschendorf1995convergence}.}
\subsection{Critical values}
The key object needed to find the critical values of our test is the Hessian $H(\gamma_0)$.
For its plug-in estimation, we need the estimators of the derivatives of $f_\text{AN}$ and $f_\text{MAR}$ with respect to $\gamma$.
Fortunately, these estimators were developed in \citet{franguridi2025raking}.
In this section, we will formulate their methodology using the operator notation that is more concise and suitable for programming.
Denote the first derivatives of $m(a,b,c)$ with respect to its arguments by
\begin{align*}
m_1(a,b,c) = b/c, \qquad m_2(a,b,c) = a/c, \qquad m_3(a,b,c) = -ab/c^2.
\end{align*}
Differentiating the raking recursion formula \eqref{eq:raking-recursion} with respect to $\gamma$ yields
\begin{align}
\nabla_\gamma f^{(t+1)}
=
m_1(f_2,h_t,\operatorname{\mathbb{I}}_1 h_t)\,\nabla_\gamma f_2
+
m_2(f_2,h_t,\operatorname{\mathbb{I}}_1 h_t)\,\nabla_\gamma h_t
+
m_3(f_2,h_t,\operatorname{\mathbb{I}}_1 h_t)\operatorname{\mathbb{I}}_1 \nabla_\gamma h_t, \label{eq:Df-recursion}
\end{align}
where
\begin{align*}
h_t &=m(f_1,f^{(t)}, \operatorname{\mathbb{I}}_2 f^{(t)}), \\
\nabla_\gamma h_t &= m_1(f_1,f^{(t)},\operatorname{\mathbb{I}}_2 f^{(t)})\,\nabla_\gamma f_1
+
m_2(f_1,f^{(t)},\operatorname{\mathbb{I}}_2 f^{(t)})\,\nabla_\gamma f^{(t)}
+
m_3(f_1,f^{(t)},\operatorname{\mathbb{I}}_2 f^{(t)})\operatorname{\mathbb{I}}_2 \nabla_\gamma f^{(t)}.
\end{align*}
The right-hand side of \eqref{eq:Df-recursion} depends on the fixed functions $f_1,f_2,\nabla_\gamma f_1,\nabla_\gamma f_2$ and the objects from the previous iteration $f^{(t)}$ and $\nabla_\gamma f^{(t)}$.
Therefore, we can write
\begin{align*}
(f^{(t+1)}, \, \nabla_\gamma f^{(t+1)}) = \mathcal{R}_{\gamma_1,\gamma_2} (f^{(t)}, \, \nabla_\gamma f^{(t)}),
\end{align*}
where $\mathcal{R}_{\gamma_1,\gamma_2}$ is the operator mapping a pair of functions of $z_1,z_2$ to another such pair of functions and which only depends on the fixed functions $f_1,f_2,\nabla_\gamma f_1,\nabla_\gamma f_2$.
The subscript $\gamma_1,\gamma_2$ emphasizes that the operator is applied at specific values $\gamma_1,\gamma_2$ where the functions and their derivatives are being calculated.
Taking into account the initial condition $f^{(0)}=f^s$, we have
\begin{align*}
(f^{(T)}, \, \nabla_\gamma f^{(T)}) = \mathcal{R}_{\gamma_1,\gamma_2}^{(T)} (f^s, \, \nabla_\gamma f^s).
\end{align*}
In particular, since $\hat f = f^{(T)}(\hat\gamma)$,
\begin{align*}
(\hat f, \, \nabla_\gamma {\hat f}) = \mathcal{R}_{\hat \gamma_1,\hat \gamma_2}^{(T)} (f^s(\hat\gamma^s), \, \nabla_\gamma f^s(\hat\gamma^s)).
\end{align*}
Finally, since $f_\text{MAR} = m(f_1,f^s, \operatorname{\mathbb{I}}_2 f^s)$, we have
\begin{align*}
\nabla_\gamma f_\text{MAR} = m_1(f_1,f^s,\operatorname{\mathbb{I}}_2 f^s)\,\nabla_\gamma f_1
+
m_2(f_1,f^s,\operatorname{\mathbb{I}}_2 f^s)\,\nabla_\gamma f^s
+
m_3(f_1,f^s,\operatorname{\mathbb{I}}_2 f^s)\operatorname{\mathbb{I}}_2 \nabla_\gamma f^s.
\end{align*}
\section{Monte Carlo simulation}\label{sec:mc-sim}
In this section, we evaluate the finite-sample performance of our test in a set of Monte Carlo simulations for a Gaussian data-generating process.
We model the latent (target) distribution $f$ of $Z=(Z_1',Z_2')'$ as
\[
Z \sim N(\mu,\Sigma),
\]
where the mean \(\mu=0\) and the covariance matrix
\[
\Sigma=
\begin{pmatrix}
I_d & \rho I_d\\
\rho I_d & I_d
\end{pmatrix},
\qquad \rho=0.5.
\]
We then model the selection probability as
\[
\operatorname{\mathbb{P}}(S=1\mid Z_1,Z_2)=\exp\!\left(\alpha-Z_1' B_1 Z_1 - Z_2' B_2' Z_2\right),
\]
with \(B_1=0\).
These choices imply that the selected density $f^s$ is also Gaussian.
The refreshment sample is drawn from the second marginal of $f$ and always has the same number of observations $n$ as the first-period sample.
We consider two DGPs, one for evaluating the size and the other for evaluating the power of our test.
The first DGP $H_0$ sets \(B_2=0\), leading to a constant selection probability, the simplest case satisfying the null hypothesis of MAR.
The second DGP $H_1$ sets
\[
B_2=0.05\,I_d,
\]
so that
\[
\operatorname{\mathbb{P}}(S=1\mid Z_1,Z_2)= \exp\!\left( \alpha-0.05\,\|Z_2\|^2 \right).
\]
In this case, selection depends on the unobserved component \(Z_2\), generating a mechanism that violates MAR but satisfies AN.
Under both $H_0$ and $H_1$, the intercept \(\alpha\) is chosen to attain the target attrition rate \(\operatorname{\mathbb{P}}(S=0)=0.2\).
To obtain the critical values for our test, we employ the following procedure that involves the nonparametric bootstrap for the parameter estimator
\[
\hat\gamma = (\hat\mu_1',\operatorname{vech}(\hat\Sigma_1),\hat\mu_2',\operatorname{vech}(\hat\Sigma_2),(\hat\mu^{\,s})',\operatorname{vech}(\hat\Sigma^{\,s}))'.
\]
First, we compute the estimate $\hat H(\hat\gamma)$ of the Hessian of the test statistic using the formulas in \Cref{sec:gaussian}.
Second, we generate bootstrap draws \(\hat \gamma^{*,b}\), $b=1,\dots, B$ and compute the bootstrap values of the test statistic
\begin{align*}
\hat T^{*,b} = \frac{1}{2} (\hat\gamma^{*,b}-\hat\gamma)' \hat H(\hat\gamma) (\hat\gamma^{*,b}-\hat\gamma), \quad b=1,\dots, B.
\end{align*}
Finally, we reject the null hypothesis of MAR if the test statistic $\hat T$ exceeds the $1-\alpha$ sample quantile of $\hat T^{*,1},\dots,\hat T^{*,B}$. We use $B=300$ bootstrap samples.
\Cref{tab:mc} reports the empirical rejection rates of our test based on $1000$ Monte Carlo replications across a range of values of $n$ and $d$.
Under the null hypothesis $H_0$, the rejection frequencies are close to the nominal 10\% level across all sample sizes and dimensions, indicating that the test provides excellent size control.
Under the alternative hypothesis $H_1$, the test exhibits substantial power even at a relatively small sample size $n=2000$.
Power then increases rapidly with sample size, reaching $100\%$ when $n=10000$.
Overall, the results suggest that the test maintains accurate size while achieving high power against nonignorable attrition alternatives.
\begin{table}[ht]
\centering
\begin{tabular}{c|cc|cc}
& \multicolumn{2}{c|}{$H_0$} & \multicolumn{2}{c}{$H_1$} \\
$n$ & $d=2$ & $d=4$ & $d=2$ & $d=4$ \\
\midrule
2000 & 0.114 & 0.093 & 0.583 & 0.676 \\
5000 & 0.101 & 0.116 & 0.937 & 0.994 \\
10000 & 0.101 & 0.090 & 1.000 & 1.000 \\
\bottomrule
\end{tabular}
\caption{Empirical rejection rates (size and power) based on $1000$ Monte Carlo replications.}
\label{tab:mc}
\end{table}
\section{Empirical illustration}\label{sec:empirical}
We apply the framework developed here to data from the Understanding America Study (UAS), a panel of approximately 15,000 U.S. respondents. The UAS started in 2014 and is still growing. This means that new batches of respondents are regularly added, both to expand the panel and to replace respondents lost to attrition. The panel collects a large amount of information in a two-year cycle, including health, health behavior and health insurance, financial status and labor market behavior, cognition, and personality. The motivation for collecting a large amount of information across different domains is that it allows much richer analysis than when a dataset concentrates on only one dimension, e.g., economics or cognition.
The example we consider is how wealth holdings are affected by income, education, and cognition. The latter variable is an aggregate obtained from factor analysis of the outcomes of five different cognitive tests. Appendix \ref{sec:app-cognition} describes the construction of the cognition variable. As noted, core information is collected in a two-year cycle. UAS respondents can join at different times, which means that respondents may answer core surveys at different points in a two-year interval.
We consider answers of respondents in the wave that spans 4 July 2021 to 4 July 2023, which, for simplicity, we refer to as Wave 1, and the wave that spans 5 July 2023 to 6 July 2025, referred to as Wave 2.
The wave that is currently in the field will only be complete by July 2027 and hence is not suitable for the analysis in this paper.
To implement the test in this empirical example, the four explanatory variables are all discretized into three categories. Education is categorized as High School or less, Some College, or College or more. Income is categorized into three terciles, and so are cognition and financial wealth.
Wave 1 has 9,725 observations, of whom 7,432 also appear in Wave 2. In addition, Wave 2 has 2,264 fresh respondents (the refreshment sample).
Our test strongly rejects the null of no selection on unobserved Wave 2 variables at any reasonable level, with a p-value of $5\cdot 10^{-6}$. This result suggests that Wave 2 attrition is related to the variables in Wave 2. The test cannot tell us which of these unobserved Wave 2 variables has affected the probability of dropping out of the sample.
Responding to surveys is a cognitively demanding activity that takes increasing effort when cognition declines. Thus, one possible explanation is that respondents with declining cognition are more likely to have dropped out of the sample.
\section{Conclusion}\label{sec:conclusion}
The commonly imposed assumption of selection on observables (MAR) may be implausible in many empirical settings and is generally untestable.
However, this assumption becomes testable when a refreshment sample is available, which is the case for many panel surveys used in empirical research.
We propose to test the MAR assumption against a weaker alternative of additive nonignorability introduced in \citet{hirano2001combining}.
Our test statistic implicitly measures the distributional deviation of the refreshment sample from the inverse probability weighted selected sample in the second period.
When the input distributions are parameterized, the computation of the test statistic and the critical values can be performed using the raking-based recursive algorithm developed in \citet{franguridi2025raking}.
\bibliographystyle{ecta}
\bibliography{references}
\onehalfspacing
\frenchspacing
\newpage
\part*{Appendix}