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.
55,054 characters
Raking for Estimation and Inference in Panel Models with Nonignorable Attrition and Refreshment
\maketitle
\begin{abstract}
\linespread{1.2}
In panel data subject to nonignorable attrition, auxiliary (refreshment) sampling may restore full identification under weak assumptions on the attrition process.
Despite their generality, these identification strategies have seen limited empirical use, largely because the implied estimation procedure requires solving a functional minimization problem for the target density. We show that this problem can be solved using the iterative proportional fitting (raking) algorithm, which converges rapidly even with continuous and moderately high-dimensional data.
This resulting density estimator is then used as input into a parametric moment condition.
We establish consistency and convergence rates for both the raking-based density estimator and the resulting moment estimator when the distributions of the observed data are parametric.
We also derive a simple recursive procedure for estimating the asymptotic variance.
Finally, we demonstrate the satisfactory performance of our estimator in simulations and provide an empirical illustration using data from the Understanding America Study panel.
\medskip
\noindent \textbf{JEL Classification:} C23
\medskip
\noindent \textbf{Keywords:} panel data, refreshment sample, additively nonignorable attrition, Kullback-Leibler divergence, raking, iterative proportional fitting, semiparametric model
\end{abstract}
\newpage
\section{Introduction}
Longitudinal data have many advantages over (repeated) cross-sections. However, a major drawback is that panel data suffer from panel attrition, in addition to the initial nonresponse common to panel and cross-sectional surveys. Panel attrition reduces the effective sample size. More concerning is that it may bias the estimates obtained from the panel survey. For instance, if we estimate the average change in household income, then households that experienced a decline in income due to the unemployment of the head of the household are more likely to move and be lost to follow-up. This nonrandom attrition that depends on an outcome variable results in an upward bias in the estimate of the average change.
Estimates that depend on outcomes in all panel waves, including the wave in which the subject drops out, are biased. Unbiasedness can be restored if the probability of attrition is restricted to depend on the outcome in the wave that the subject drops out, but not on the outcomes in previous waves \citep{hausman1979attrition}. Estimates can also be unbiased if attrition depends on the outcomes in the waves prior to the wave that the subject drops out, but not on the (unobserved) outcome in the wave that the subject leaves the panel \citep{rubin1976inference,little2019statistical}. The latter case is referred to as Missing At Random (MAR), with MCAR (Missing Completely At Random) being the special case in which the attrition is purely random.
The unselected joint distribution of the outcome variables is, under the assumption of MAR, nonparametrically just identified. Therefore, additional information is needed to relax MAR if we want to allow the attrition to depend on the unobserved outcome in the wave in which the subject drops out. A source of additional information is a supplementary sample drawn to compensate for the loss of subjects due to drop-out. The idea to add such a sample dates back at least to \cite{kish1959replacement}. In panel surveys, these samples are drawn in the second and later waves and are called refreshment samples \citep{ridder1992empirical}. \cite{deng2013handling} and \cite{watson2021refreshment} survey the use of refreshment samples in panel surveys.
A key identification result in panel surveys with selective attrition and refreshment samples was established by \cite{hirano2001combining}. They considered a selective panel with two waves supplemented by a refreshment sample in the second wave. In this setup, the marginal distributions of the first-wave variables and of the second-wave variables are directly identified. They proposed estimating the unselected joint distribution of the variables by the distribution that minimizes the distance to the joint selective observed distribution of the variables in the two waves under the constraints on the marginal distributions.
Because the distance measure is a convex functional and the restrictions are linear, this functional minimization problem has a unique solution. This solution can be used to obtain a semiparametric estimation procedure. Consider an estimator that depends on the outcomes in both waves of the panel. An example is a Fixed Effects (FE) regression where we regress the change in the dependent variable on the change in the independent variables. If there is no selective attrition that depends on the first- and second-wave dependent variables, then the FE estimator of the regression coefficients is unbiased. The estimator is the solution to a sample moment condition. If the attrition does depend on the first and second-wave dependent variables, then the FE estimator is biased.
The first-order condition of the functional minimization problem that defines the estimator of the unselected joint distribution has a unique solution for the probability of observation that is a function of the outcome variables in both waves of the panel. This probability is additively separable in the outcomes for the two waves, which is intuitive because we only observe the selected joint distribution. Under this restriction, the probability is nonparametrically identified.
The probability of observation is used in a weighted generalized moment estimator. \cite{bhattacharya2008inference} studies the properties of this semiparametric estimator, and \cite{hoonhout2019nonignorable} develops the estimator for panel data with three or more waves. \cite{hirano1998combining} and \cite{deng2013handling} choose parametric models for the probability of observation and the binary outcome variables and use Markov Chain Monte Carlo to compute posteriors for the parameters. \cite{deng2013handling} use multiple imputation inference \citep{rubin2004multiple}. A disadvantage of the weighted moment and Bayesian estimators is that they are computationally demanding and require close monitoring by the statistician. \cite{bhattacharya2008inference} reports that the moment estimator does not always converge.\footnote{See \cite{Franguridi2024} for further comments on \cite{bhattacharya2008inference}.} Preferably, we would like a procedure that does not require close monitoring.
In this paper, we follow the generalized moment approach, but instead of estimating weights, we average over the unselected joint distribution that is estimated by minimizing the distance to the observed joint distribution under the constraints on the marginal unselected distributions, a problem that has a unique solution. We exploit the availability of an efficient algorithm for this constrained minimization problem if the distance is the Kullback-Leibler divergence. This algorithm is iterative and involves only arithmetic operations. With discrete data, the algorithm is known as \emph{raking} or \emph{iterative proportional fitting} and was first proposed by \cite{deming1940least}.
Starting from the observed joint distribution of the outcomes, often only a few steps are required for convergence to the unique solution.\footnote{\cite{ruschendorf1995convergence} gave conditions under which the algorithm converges.}
Inputs in the algorithm are the observed joint distribution of the variables in the two waves and the marginal distributions in the first and second waves, with the latter obtained from the refreshment sample. The output is an estimate of the unselected joint distribution of the variables in the two waves. We can estimate the input densities either parametrically or nonparametrically. In this paper, we take the parametric approach. As a by-product, the raking algorithm yields a recursive approximation of the Jacobian and the asymptotic variance of the generalized moment estimator.
The paper is organized as follows.
\Cref{sec:framework} introduces the framework.
\Cref{sec:estimation} considers estimation.
\Cref{sec:inference} discusses inference.
\Cref{sec:numerical} provides details for numerical implementation.
\Cref{sec:mc} reports a simulation study of our estimator.
\Cref{sec:empirical} contains the empirical illustration.
\Cref{sec:conclusion} concludes.
\section{Framework}\label{sec:framework}
We consider a two-period panel data model with attrition and refreshment as in \citet{hirano2001combining}.
Let $Z_{it}$ be a $d$-dimensional vector of all the time-varying variables (both outcomes and covariates) for unit $i$ in period $t \in \{1,2\}$, and let $X_i$ be a vector of all the time-invariant variables for unit $i$.
We often drop the subscript $i$ for convenience.
Our analysis can be made conditional on $X_i$, and hence, throughout the rest of the paper, we drop the variable $X_i$.
In the first period, all the $n$ sampled units are observed: $Z_{i1},$ $i=1,\dots,n$, is a random sample from the marginal distribution of $Z_1$ with density $f_1$ .
In the second period, there is \emph{attrition}: only the units with $W_i=1$ stay in the sample, while the remaining units drop out.
To compensate for this attrition, a \emph{refreshment sample} $Z_{i2}^r$, $i =1,\dots, n_r$, i.e., a random sample from the marginal distribution of $Z_2$ with density $f_2$, is drawn in period 2.
We assume that the refreshment sample is independent of the rest of the data, and its size is of the same order of magnitude as the size of the first-period sample, i.e. $n/n_r \to c\in(0,\infty)$.
Our goal is to estimate and perform inference on the $d_\theta$-dimensional parameter $\theta$ satisfying the moment condition
\begin{align}
\operatorname{\mathbb{E}}_f\left[ \varphi\left( Z_{1},Z_{2},\theta\right) \right] = \iint \varphi\left( z_{1},z_{2},\theta\right) f(z_1,z_2) \, d\mu(z_1,z_2) = 0, \label{eq:moment-theta-unselect}
\end{align}
where $\varphi$ is a $d_{\theta}$-dimensional moment function, and the expectation is taken with respect to the \emph{unconditional} joint density $f$ of $(Z_1,Z_2)$ relative to some dominating probability measure $\mu$ on $\mathbb{R}^{2d}$. The latter allows us to handle both continuous and discrete variables simultaneously.\footnote{We consider the case where $\theta$ is just identified. If the parameter is overidentified, then (\ref{eq:moment-theta-unselect}) is the first-order condition of a quadratic minimization problem.}
Without further restrictions, $f$ is not point identified, and hence neither is $\theta$.
Moreover, the bounds on the identified set are expected to be uninformative: any distribution $f$ with marginals $f_1$ and $f_2$ is consistent with the data.
To understand the intuition, consider the identity for the target distribution with the balanced panel distribution $f^w(z_1,z_2):= f(z_1,z_2|W=1)$ and the probability of observation
\begin{align}
f(z_1,z_2) = \frac{\operatorname{\mathbb{P}}(W=1)}{\operatorname{\mathbb{P}}(W=1|Z_1=z_1,Z_2=z_2)} \cdot f^w(z_1,z_2). \label{eq:key-identity}
\end{align}
Consider the case of discrete data for simplicity.
The only information available to identify the weights $\operatorname{\mathbb{P}}(W=1)/\operatorname{\mathbb{P}}(W=1|Z_1=z_1,Z_2=z_2)$, which generally have $|\operatorname{supp} Z_1|\cdot|\operatorname{supp} Z_2|$ degrees of freedom, is the marginal distributions $f_1$ and $f_2$, which only have $|\operatorname{supp} Z_1|+|\operatorname{supp} Z_2|$ degrees of freedom.
This establishes the lack of identification and suggests imposing restrictions on the weights.
\citet{hirano2001combining} impose the following separability condition on the weights.\footnote{\citet{franguridi2025inference} drop this separability assumption, and \citet{franguridi2025set} further consider the case where refreshment samples are not available. In both cases, the parameters of interest become partially identified.}
\begin{assumption}[additive nonignorability]\label{as:AN}
For a known, increasing function $G$ and unknown functions $k_1,k_2$,
\begin{align*}
\operatorname{\mathbb{P}}(W=1|Z_1=z_1,Z_2=z_2) = G\left( k_1(z_1)+k_2(z_2) \right).
\end{align*}
\end{assumption}
Throughout this paper, we further require the following.
\begin{assumption}[exponential link function]\label{as:AN-exp}
$G(x)=\exp(x)$, $x\in \mathbb{R}$.
\end{assumption}
This link function corresponds to adopting the Kullback-Leibler divergence as the distance mentioned above. This assumption ensures that the functional minimization problem introduced below can be solved by raking.
Since the assumption allows for cross-period interactions in attrition, it does not constitute a significant loss of generality.
Under \cref{as:AN,as:AN-exp}, the key identity \eqref{eq:key-identity} becomes
\begin{align}
f(z_1,z_2) = e^{k_1(z_1)+k_2(z_2)} \cdot f^w(z_1,z_2), \label{eq:key-identity-exp}
\end{align}
where we redefine $k_1,k_2$ appropriately to absorb the constant $\operatorname{\mathbb{P}}(W=1)$ and switch the sign on $k_1,k_2$.
The marginal densities $f_1,f_2$ can now be used to identify $k_1,k_2$, but it is convenient to estimate $f$ using the characterization in \citet{hirano2001combining}.
\section{Estimation}\label{sec:estimation}
\citet{hirano2001combining} show that Assumptions \ref{as:AN} and \ref{as:AN-exp} imply that $f$ is the Kullback-Leibler projection of $f^w$ on the set $\Pi(f_1,f_2)$ of distributions with marginals $f_1$ and $f_2$, i.e.,\footnote{This projection is also known as the \emph{Schr\"{o}dinger bridge} with marginals $f_1$, $f_2$ and the reference distribution $f^w$. It is named after the physicist Erwin Schr\"{o}dinger, who first considered it in \citet{schrodinger1931umkehrung}.}
\begin{align}
f = \arg\min_{\tilde f \in \Pi(f_1,f_2)} \operatorname{KL}(\tilde f, f^w) \label{eq:KL-representation},
\end{align}
where
\[
\operatorname{KL}(\tilde f, f^w) := \iint \tilde f(z_1,z_2) \log \left( \frac{\tilde f(z_1,z_2)}{f^w(z_1,z_2)} \right) \, d\mu(z_1,z_2).
\]
Equivalently, $f$ minimizes the Kullback-Leibler distance to $f^w$ for densities $f$ with given marginals $f_1,f_2$. Because the distance $\operatorname{KL}(\tilde f, f^w)$ viewed as a function of $\tilde f$ is strictly convex in $\tilde f$, the minimum is unique. The density $f$ is a functional of the densities $f_1,f_2$, and $f^w$, which are directly estimable from the data.
The sample analog estimator of $f$ is
\begin{align}
\hat f = \arg\min_{\tilde f \in \Pi(\hat f_1,\hat f_2)} \operatorname{KL}(\tilde f, \hat f^w) \label{eq:KL-representation-1},
\end{align}
where $\hat f_1$, $\hat f_2$, and $\hat f^w$ are density estimators.
The estimator $\hat f$ is then a functional of $\hat f_1, \hat f_2$, and $\hat f^w$.
At first sight, solving the constrained functional optimization problem \eqref{eq:KL-representation} that defines the estimator seems computationally infeasible. However, there exists an algorithm for solving this problem that is guaranteed to converge and requires little beyond basic arithmetic operations.
To describe the algorithm, let $\Pi_1$ and $\Pi_2$ be the operators calculating the first and second marginals, respectively, i.e. $\Pi_1 f(z_1) := \int f(z_1,z_2) \, d \mu_2(z_2) $ and $\Pi_2 f(z_2) := \int f(z_1,z_2) \, d \mu_1(z_1)$.
Given the marginal densities $f_1,f_2$, define the operators $\phi_{1,f_1}$, $\phi_{2,f_2}$, and $\phi_{f_1,f_2}$ by
\begin{align*}
\phi_{1,f_1}(f) := \frac{f_1 \cdot f}{\Pi_1 f},
\quad
\phi_{2,f_2}(f) := \frac{f_2 \cdot f}{\Pi_2 f},
\quad
\phi_{f_1,f_2} := \phi_{2,f_2}(\phi_{1,f_1}(f)).
\end{align*}
where $f$ is a general density, not necessarily with marginal densities $f_1,f_2$. Note that $\phi_{1,f_1}(f)$ is the Kullback-Leibler (KL) projection of $f$ on the set of distributions with the first marginal $f_1$, and similarly for $\phi_{2,f_2}(f)$ \citep[p. 1164]{ruschendorf1995convergence}.
To see this, let us show that the KL projection of $f$ onto the set of distributions with one fixed marginal $f_1$ has to preserve conditional distributions $f(z_2|z_1)$.
Indeed, write an arbitrary distribution $g$ as $g(z_1,z_2)=g(z_1)g(z_2|z_1)$.
Since
\begin{align*}
\frac{g(z_1,z_2)}{f(z_1,z_2)} = \frac{g(z_1) g(z_2|z_1)}{f(z_1)f(z_2|z_1)},
\end{align*}
we have
\begin{align*}
\operatorname{KL}(g\,||\,f) &= \int g(z_1,z_2)\log \frac{g(z_1,z_2)}{f(z_1,z_2)}\,dz_1 dz_2 \\
&= \iint g(z_1,z_2) \log \frac{g(z_1)}{f(z_1)}\, dz_1 dz_2 + \iint g(z_1,z_2) \log \frac{g(z_2|z_1)}{f(z_2|z_1)}\, dz_1 dz_2 \\
&= \int g(z_1)\log \frac{g(z_1)}{f(z_1)}\, dz_1 + \int \left[\int g(z_2|z_1) \log \frac{g(z_2|z_1)}{f(z_2|z_1)} \, dz_2 \right] g(z_1) dz_1 \\
&= \operatorname{KL}(g_1\,||\, f_1) + \int \operatorname{KL}(g(\cdot|z_1) \,||\, f(\cdot|z_1)) g(z_1) \, dz_1.
\end{align*}
Since we are interested in the class of distributions $g$ with fixed first marginal, the first term is constant over this class, and hence is irrelevant.
On the other hand, the second term is nonnegative and is minimized at $g(\cdot|z_1)=f(\cdot|z_1)$ for all $z_1$ in the support of $Z_1$.
Hence, the KL projection is a distribution $g$ with the first marginal $g_1=f_1$ and the conditional distributions equal to those of $f$. There is only one such distribution, and it is exactly $\phi_{1,f_1}(f)$ defined above.
Further discussion can be found in \citet[paragraph 2, p. 155]{csiszar1975divergence}, \citet{ireland1968contingency}, and \citet{kullback1968probability}.
The algorithm requires assumptions on the supports of the involved distributions.
\begin{assumption}\label{as:abs-cont}
\begin{subassumption}
\item $\mu = \bigotimes_{k=1}^{2d} \mu_k$, where $\mu_k$ is either the Lebesgue measure on $\mathbb{R}$ or the counting measure on a countable subset of $\mathbb{R}$, and the joint distribution of $Z_1,Z_2$ is absolutely continuous w.r.t. $\mu$.
\item The support of the joint distribution of $Z_1,Z_2$ in the balanced panel is contained in the product set of the supports of $Z_1$ and $Z_2$ in the balanced panel.
\item The support of the unselected distribution of $Z_1$ is contained in the support of $Z_1$ in the balanced panel and similarly for $Z_2$.
\item There exists $c>0$ such that $\frac{f^w(z_1,z_2) f_2^w(z_2)}{f_2(z_2)} \ge c$ for $f^w$-a.e. $(z_1,z_2)$.
\end{subassumption}
\end{assumption}
\cref{as:abs-cont}(i) states that the variables in the data are either discrete or continuous.
\cref{as:abs-cont}(ii) prevents $f^w$ from concentrating on lower-dimensional subsets of $\mathbb{R}^{2d}$; as long as the latter does not happen, the condition is satisfied.
\cref{as:abs-cont}(iii) posits that the support of the period marginals of the balanced panel distribution cannot be smaller than the support of the corresponding marginals of the target distribution.
Finally, \cref{as:abs-cont}(iv) is satisfied if $f^w$ is bounded away from zero on its support.
Theorem 3.5 in \citet{ruschendorf1995convergence} implies the following result.
\begin{thm}\label{thm:raking}
Under \cref{as:AN,as:AN-exp,as:abs-cont}, the sequence of functions
\begin{align*}
\hat f^{(t)} = \hat \phi_{f_1,f_2}\left (\hat f^{(t-1)}\right ), \quad \hat f^{(0)}= \hat f^w,
\end{align*}
converges to the solution $\hat f$ of \eqref{eq:KL-representation-1} in $L^1$ as $t \to \infty$.
\end{thm}
\begin{proof}
See \Cref{app:thm1}.
\end{proof}
The iterative equation of the theorem is the general data version of \emph{raking}, also known as \emph{iterative proportional fitting} or \emph{Sinkhorn's algorithm}.\footnote{
Raking was introduced by \citet{deming1940least} to adjust contingency tables to given marginals and was generalized to the continuous case by \citet{ireland1968contingency} and \citet{kullback1968probability}.}
It is widely used in survey statistics and, since the seminal contribution of \citet{cuturi2013sinkhorn}, has attracted substantial attention in machine learning due to its application to solving optimal transport problems with entropic regularization.
\cref{thm:raking} suggests the following two-step estimation strategy.
\begin{enumerate}
\item Use (\ref{eq:KL-representation}) with estimated densities $\hat f_1$, $\hat f_2$, and $\hat f^w$ to obtain an estimator $\hat f$ computed using raking as in \cref{thm:raking}:
\begin{align*}
\hat f = \phi_{\hat f_1,\hat f_2}^{(T)}(\hat f^w),
\end{align*}
where the number of iterations $T$ is determined by a stopping rule, e.g., the maximum difference between subsequent iterations being smaller than some tolerance level.
\item Compute $\hat \theta$ as a solution of the moment restriction
\begin{align}
\iint \varphi(z_1,z_2,\hat\theta) \hat f(z_1,z_2) \, d\mu(z_1,z_2) = 0. \label{eq:thetahat}
\end{align}
\end{enumerate}
We now show consistency of $\hat f$ and $\sqrt n$-consistency of $\hat \theta$.
For simplicity, in what follows we assume that $\mu$ is the Lebesgue measure, i.e., all the variables in the data are continuous.
The general case can be handled similarly at the expense of more cumbersome notation.
Denote by $\mathcal{Z} \subset \mathbb{R}^{2d}$ the support of $f_1\times f_2$ and let $\mathcal{K}$ be the set of uniformly bounded, additively separable functions on $\mathcal{Z}$, i.e.,
\[
\mathcal{K} := \left\{ k:\mathcal{Z} \to \mathbb{R} \text{ s.t. } \|k\|_\infty \le K \text{ and } k(z_1,z_2)=k_1(z_1)+k_2(z_2) \text{ for some } k_1,k_2 \right\},
\]
where $K>0$ is a finite constant and $\|k\|_\infty = \sup_{z\in\mathcal{Z}} |k(z)|$ is the uniform norm.
Notice that the functions $k_1$ and $k_2$ are unique only up to addition and subtraction of a constant, but this does not play any role in our arguments.
We make the following assumptions.
\begin{assumption}[additive nonignorability + exponential link function]\label{as:AN-alt}
The function $k_0(z_1,z_2) := \log (f(z_1,z_2)/f^w(z_1,z_2))$ belongs to the set $\mathcal{K}$.
\end{assumption}
\begin{assumption}[identification of $\theta$]\label{as:identification}
\begin{subassumption}
\item The parameter space $\Theta \subset \mathbb{R}^{d_{\theta}}$ is compact.
\item $\operatorname{\mathbb{E}}_f [\varphi(Z_1,Z_2,\theta)]=0$ implies $\theta=\theta_0$.
\item $\varphi(z_1,z_2,\theta)$ is twice continuously differentiable in $\theta$ on $\Theta$ for all $z_1,z_2$.
\item $\operatorname{\mathbb{E}} \left[\frac{\partial}{\partial \theta'}\varphi(Z_1,Z_2,\theta_0)\right]$ is nonsingular.
\item $\varphi(z,\theta)$ and $\frac{\partial \varphi}{\partial\theta'}(z,\theta)$ are bounded on $\mathcal{Z} \times \Theta$.
\end{subassumption}
\end{assumption}
\begin{assumption}[densities]\label{as:densities-alt}
\begin{subassumption}
\item The densities $f^w(z_1,z_2;\gamma_w),f_1(z_1;\gamma_1),f_2(z_2;\gamma_2)$ belong to a parametric family and are continuously differentiable with respect to their parameters $\gamma_1,\gamma_2, \gamma_w$ that are identifiable.
\item The gradients of densities are integrable locally uniformly in the parameters in the following sense: there exists a constant $C$ such that
\begin{align*}
\max\left\{\iint\left\|\frac{\partial f^w}{\partial \gamma_w}(z_1,z_2;\bar\gamma_w)\right\|dz_1 dz_2, \,\, \int\left\|\frac{\partial f_1}{\partial \gamma_1}(z_1;\bar\gamma_1)\right\|dz_1, \,\,
\int\left\|\frac{\partial f_2}{\partial \gamma_2}(z_2;\bar\gamma_2)\right\|dz_2 \right\} \le C
\end{align*}
for all $\bar\gamma_w,\bar\gamma_1,\bar\gamma_2$ in a neighborhood of the population values of $\gamma_w,\gamma_1,\gamma_2$.
\end{subassumption}
\end{assumption}
\Cref{as:AN-alt} is equivalent to imposing both \Cref{as:AN,as:AN-exp}.
\Cref{as:identification} is a standard assumption ensuring identification of $\theta$.
\Cref{as:densities-alt}(i) imposes a parametric model on the observed data.
A sufficient condition for \Cref{as:densities-alt}(ii) is that the information matrices associated with the parametric densities $f^w,f_1,f_2$ are bounded in some neighborhood of the population values of $\gamma_w,\gamma_1,\gamma_2$.
We estimate $\gamma_1,\gamma_2, \gamma_w$ by maximum likelihood and denote the corresponding density estimators $\hat f^w(z_1,z_2)=f^w(z_1,z_2;\hat \gamma_{w})$, $\hat f_1(z_1)=f_1(z_1;\hat \gamma_{1})$, $\hat f_2(z_2)=f_2(z_2;\hat \gamma_{2})$.
We omit the arguments $\gamma_1,\gamma_2,\gamma_w$ when the densities are evaluated at their population parameter values.
\begin{assumption}[maximum likelihood estimators]\label{as:mle}
$\hat\gamma_1, \hat\gamma_2,\hat\gamma_w$ are $\sqrt{n}$-consistent with the expectation bounds
\begin{align*}
\operatorname{\mathbb{E}} \|\hat\gamma_1-\gamma_1\| = O(1/\sqrt{n}), \quad \operatorname{\mathbb{E}} \|\hat\gamma_2-\gamma_2\| = O(1/\sqrt{n}), \quad \operatorname{\mathbb{E}} \|\hat\gamma_w-\gamma_w\| = O(1/\sqrt{n}).
\end{align*}
\end{assumption}
A key idea allowing us to establish the $\sqrt{n}$-convergence of $\hat f$ and $\hat\theta$ is to switch from the primal characterization of $f$ via a \emph{constrained} program \eqref{eq:KL-representation} to the dual characterization \eqref{eq:key-identity-exp} with $k_1,k_2$ solving an \emph{unconstrained} program.
This duality result is not new (see, e.g., Section 4 in \cite{bhattacharya2006iterative} and Theorem 2.3 in \cite{bhattacharya1995general}), but we provide a proof to keep the presentation self-contained.
\begin{lem}[duality]
Denote
\begin{align*}
\mathbb{M}(k) &:= \iint e^{k_1(z_1)+k_2(z_2)} f^w(z_1,z_2)\,dz_1 dz_2 - \int k_1(z_1) f_1(z_1) \, dz_1 - \int k_2(z_2) f_2(z_2) \, dz_2, \\
\mathbb{M}_n(k) &:= \iint e^{k_1(z_1)+k_2(z_2)} \hat f^w(z_1,z_2)\,dz_1 dz_2 - \int k_1(z_1) \hat f_1(z_1) \, dz_1 - \int k_2(z_2) \hat f_2(z_2) \, dz_2.
\end{align*}
Then the population density $f$ (the estimator $\hat f$) is related to the population density $f^w$ (the estimator $\hat f^w$) via
\begin{align}
f(z_1,z_2) &= e^{k_0(z_1,z_2)} f^w(z_1,z_2), \label{eq:f-via-k0} \\
\hat f(z_1,z_2) &= e^{\hat k(z_1,z_2)} \hat f^w(z_1,z_2),
\end{align}
where
\begin{align}
k_0 &= \arg\min_{k \in \mathcal{K}} \mathbb{M}(k), \label{eq:k0-via-dual} \\
\hat k &= \arg\min_{k \in \mathcal{K}} \mathbb{M}_n(k).
\end{align}
\end{lem}
\begin{proof}
See \Cref{app:lem1}.
\end{proof}
Finally, we establish $\sqrt{n}$-consistency of the raking estimator of $f$ and the estimator of $\theta$.
\begin{thm}
Suppose that \cref{as:AN-alt,as:identification,as:abs-cont,as:densities-alt,as:mle} hold.
Then $\|\hat f-f\|_1 = O_p \left (1/\sqrt n \right)$ and $\hat\theta-\theta_0 = O_p\left (1/\sqrt n \right )$, where $||f||_1 =\int |f(z)| dz$.
\end{thm}
\begin{proof}
See \Cref{app:thm2}.
\end{proof}
\section{Inference}\label{sec:inference}
The asymptotic distribution of the moment estimator of $\theta$ depends on the asymptotic distribution of the estimator (\ref{eq:KL-representation-1}) of the pseudo-parameter $f$ defined in (\ref{eq:KL-representation}). That estimator is computed as in \cref{thm:raking}. The steps are analogous to the steps in the derivation of the asymptotic distribution of an MLE that is computed by a numerical optimization algorithm. The algorithm is recursive, and the MLE is the result after a finite number of iterations. The first-order condition is solved for the MLE, and the asymptotic distribution of this solution is the asymptotic distribution of the MLE. The estimator (\ref{eq:KL-representation-1}) is computed by a finite number of iterations as in \cref{thm:raking}. The asymptotic distribution of the estimator (\ref{eq:KL-representation-1}) is therefore equal to $T$-fold recursion in \cref{thm:raking}.
Under \Cref{as:densities-alt}, the marginal densities of $z_1,z_2$ and the joint density of the observed $z_1,z_2$ are in a parametric family. Their parameters can be estimated by MLE. The asymptotic distribution of the estimator of $f$ can therefore be derived using the delta method. Define $\gamma=(\gamma_1',\gamma_2',\gamma_w')'$ and let $D_\gamma$ be the derivative operator with respect to $\gamma$. Define $m(a,b,c) := ab/c$ and let $m_1=b/c$, $m_2=a/c$, and $m_3=-ab/c^2$ be its derivatives with respect to $a,$ $b$, and $c$, respectively. Let $V(\hat\gamma)$ be the asymptotic variance of an estimator $\hat\gamma$ of $\gamma$, such as the inverse information matrix when $\hat\gamma$ is the maximum likelihood estimator.
\begin{thm}\label{thm:asy-var}
The asymptotic variance of $\hat\theta$ is $V(\hat\theta) = A_0 \cdot V(\hat\gamma) \cdot A_0'$, where
\begin{align*}
A_0 = \left( D_\theta \iint \varphi(z_1,z_2;\theta_0) f(z_1,z_2) \, dz_1 dz_2 \right)^{-1} \iint \varphi(z_1,z_2;\theta_0) D_\gamma f(z_1,z_2) \,dz_1 dz_2.
\end{align*}
The matrix $A_0$ can be estimated by
\begin{align*}
\hat A_0 = \left( D_\theta \iint \varphi(z_1,z_2;\hat \theta) \hat f(z_1,z_2) \, dz_1 dz_2 \right)^{-1} \iint \varphi(z_1,z_2;\hat \theta) D_\gamma \hat f(z_1,z_2) \,dz_1 dz_2.
\end{align*}
Here $\hat f(z_1,z_2)=\hat f^{(T)}(z_1,z_2)$ is the $T$-th iterate of
\begin{equation}\label{fhat-iter}
\hat f^{(t)}(z_1,z_2)= \frac{f_2(z_2;\hat \gamma_2)m \left (f_1(z_1;\hat \gamma_1), \hat f^{(t-1)}(z_1,z_2), \Pi_1 \hat f^{(t-1)}(z_1,z_2)\right )}{\Pi_2 m\left ( f_1(z_1; \hat \gamma_1), \hat f^{(t-1)}(z_1,z_2), \Pi_1 \hat f^{(t-1)}(z_1,z_2) \right )}
\end{equation}
with $\hat f^{(0)}(z_1,z_2)=f^w(z_1,z_2; \hat \gamma_w)$,
\[
D_{\gamma} \hat f(z_1,z_2)=D_{\gamma}\hat f^{(T)}(z_1,z_2)
\]
and
\[
D_{\gamma}\hat f^{(t)}(z_1,z_2)= \left ( \begin{array}{c}
D_{\gamma_1}\hat f^{(t)}(z_1,z_2) \\
D_{\gamma_2}\hat f^{(t)}(z_1,z_2)\\
D_{\gamma_w}\hat f^{(t)}(z_1,z_2)
\end{array} \right )
\]
with iteration
\begin{align*}
&D_{\gamma_1}\hat f^{(t)}(z_1,z_2) = \\
&\frac{f_2(z_2;\hat \gamma_2) \left ( m_1 (z_1,z_2) D_{\gamma_1}f_1(z_1;\hat \gamma_1)+m_2 (z_1,z_2 ) D_{\gamma_1}\hat f^{(t-1)}(z_1,z_2)+m_3(z_1,z_2)\Pi_1 D_{\gamma_1}\hat f^{(t-1)}(z_1,z_2)\right )}{\Pi_2 m\left ( f_1(z_1; \hat \gamma_1), \hat f^{(t-1)}(z_1,z_2), \Pi_1 \hat f^{(t-1)}(z_1,z_2) \right )}-\\
&\frac{f_2(z_2;\hat \gamma_2)m \left (f_1(z_1;\hat \gamma_1), \hat f^{(t-1)}(z_1,z_2), \Pi_1 \hat f^{(t-1)}(z_1,z_2)\right )}{\left (\Pi_2 m\left ( f_1(z_1; \hat \gamma_1), \hat f^{(t-1)}(z_1,z_2), \Pi_1 \hat f^{(t-1)}(z_1,z_2) \right )\right)^2} \times \\
&\Pi_2\left ( m_1 (z_1,z_2) D_{\gamma_1}f_1(z_1;\hat \gamma_1)+m_2 (z_1,z_2 ) D_{\gamma_1}\hat f^{(t-1)}(z_1,z_2)+m_3(z_1,z_2)\Pi_1 D_{\gamma_1}\hat f^{(t-1)}(z_1,z_2)\right),
\end{align*}
\begin{align*}
&D_{\gamma_2}\hat f^{(t)}(z_1,z_2) = \\
&\frac{m \left (f_1(z_1;\hat \gamma_1), \hat f^{(t-1)}(z_1,z_2), \Pi_1 \hat f^{(t-1)}(z_1,z_2)\right ) D_{\gamma_2}f_2(z_2;\hat \gamma_2)}{\Pi_2 m\left ( f_1(z_1; \hat \gamma_1), \hat f^{(t-1)}(z_1,z_2), \Pi_1 \hat f^{(t-1)}(z_1,z_2) \right )}+\\
&\frac{f_2(z_2;\hat \gamma_2)\left( m_2(z_1,z_2)D_{\gamma_2}\hat f^{(t-1)}(z_1,z_2)+m_3(z_1,z_2)\Pi_1 D_{\gamma_2}\hat f^{(t-1)}(z_1,z_2)\right)}{\Pi_2 m\left ( f_1(z_1; \hat \gamma_1), \hat f^{(t-1)}(z_1,z_2), \Pi_1 \hat f^{(t-1)}(z_1,z_2) \right )}-\\
&\frac{f_2(z_2;\hat \gamma_2)m \left (f_1(z_1;\hat \gamma_1), \hat f^{(t-1)}(z_1,z_2), \Pi_1 \hat f^{(t-1)}(z_1,z_2)\right )}{\left (\Pi_2 m\left ( f_1(z_1; \hat \gamma_1), \hat f^{(t-1)}(z_1,z_2), \Pi_1 \hat f^{(t-1)}(z_1,z_2) \right )\right)^2} \times \\
&\Pi_2\left (m_2 (z_1,z_2 ) D_{\gamma_2}\hat f^{(t-1)}(z_1,z_2)+m_3(z_1,z_2)\Pi_1 D_{\gamma_2}\hat f^{(t-1)}(z_1,z_2)\right),
\end{align*}
\begin{align*}
&D_{\gamma_w}\hat f^{(t)}(z_1,z_2) = \\
&\frac{f_2(z_2;\hat \gamma_2) \left ( m_2 (z_1,z_2 ) D_{\gamma_w}\hat f^{(t-1)}(z_1,z_2)+m_3(z_1,z_2)\Pi_1 D_{\gamma_w}\hat f^{(t-1)}(z_1,z_2)\right )}{\Pi_2 m\left ( f_1(z_1; \hat \gamma_1), \hat f^{(t-1)}(z_1,z_2), \Pi_1 \hat f^{(t-1)}(z_1,z_2) \right )}-\\
&\frac{f_2(z_2;\hat \gamma_2)m \left (f_1(z_1;\hat \gamma_1), \hat f^{(t-1)}(z_1,z_2), \Pi_1 \hat f^{(t-1)}(z_1,z_2)\right )}{\left (\Pi_2 m\left ( f_1(z_1; \hat \gamma_1), \hat f^{(t-1)}(z_1,z_2), \Pi_1 \hat f^{(t-1)}(z_1,z_2) \right )\right)^2} \times \\
&\Pi_2\left ( m_2 (z_1,z_2 ) D_{\gamma_w}\hat f^{(t-1)}(z_1,z_2)+m_3(z_1,z_2)\Pi_1 D_{\gamma_w}\hat f^{(t-1)}(z_1,z_2)\right).
\end{align*}
In the above equations,
\[
m_k(z_1,z_2)=m_k (f_1(z_1;\hat \gamma_1), \hat f^{(t-1)}(z_1,z_2), \Pi_1 \hat f^{(t-1)}(z_1,z_2)), \ \ \ k=1,2,3,
\]
and the initial conditions are
\[
D_{\gamma}\hat f^{(0)}(z_1,z_2)= \left ( \begin{array}{c}
D_{\gamma_1}\hat f^{(0)}(z_1,z_2) \\
D_{\gamma_2}\hat f^{(0)}(z_1,z_2)\\
D_{\gamma_w}\hat f^{(0)}(z_1,z_2)
\end{array} \right )=
\left ( \begin{array}{c}
0\\
0\\
D_{\gamma_w}f^w(z_1,z_2; \hat \gamma_w)\\
\end{array} \right ).
\]
\end{thm}
\begin{proof}
Differentiation of (\ref{fhat-iter}) with respect to $\gamma$ yields the iteration.
\end{proof}
\section{Numerical implementation}\label{sec:numerical}
Both raking iterations and moment conditions require computing integrals.
To do it quickly and reliably, we suggest using Monte Carlo integration on a fixed rectangular grid in terms of the variables $z_1,z_2$ that we generate randomly from distributions that dominate $Z_1$ and $Z_2$.
To describe our procedure, suppose both $f_1$ and $f_2$ are absolutely continuous with respect to a dominating measure $\mu$. For instance, when the data are one discrete and one continuous variable, $\mu$ is the product of a counting measure and the Lebesgue measure. Notice that integrals in raking iterations and moment conditions are calculated with respect to $\mu$.
Let $\phi_1$ and $\phi_2$ be densities on $\mathbb{R}^d$ such that $f_1$ and $f_2$ are absolutely continuous w.r.t. $\phi_1$ and $\phi_2$, respectively.
Typically, $\phi_1=\phi_2$.
Let $z_{11},\dots, z_{1S_1}$ and $z_{21},\dots,z_{2S_2}$ denote independent draws from distributions $\phi_1$ and $\phi_2$, respectively.
Fix the rectangular grid $\{z_{11},\dots,z_{1S_1}\} \times \{z_{21},\dots,z_{2S_2}\}$ from now on.
We only need to calculate the integrals involved in raking and possibly moment conditions on this grid.
For illustration, suppose $Z=(Y,X)'$, where $X$ is discrete with support $\mathcal{X} = \{x_1,\dots,x_R\}$ and $Y$ is continuous.
Let $\phi_X$ be a probability mass function on $\mathcal{X}$, i.e. $\phi_X(x) >0$ for $x\in\mathcal{X}$, $\phi_X(x)=0$ for $x\notin \mathcal{X}$, and $\phi_X(x_1)+\dots+\phi_X(x_R)=1$.
Let $\phi_Y$ be a probability density on $\mathbb{R}^1$, e.g., a Gaussian density.
Then, for period 1, we can independently sample $S_1$ draws from $\phi_X$ and $S_1$ draws from $\phi_Y$, obtaining a combined sample $\mathcal{Z}_1= \{(x_{11},y_{11}),\dots,(x_{1S_1},y_{1S_1})\}$.
Similarly, for period 2, we obtain a sample $\mathcal{Z}_2= \{(x_{21},y_{21}),\dots,(x_{2S_2},y_{2S_2})\}$ using the same $\phi_X$ and $\phi_Y$. Then $\phi_1=\phi_2$ and the grid is $\mathcal{Z}_1 \times \mathcal{Z}_2$.
The density $\phi=\phi_1=\phi_2$ with respect to the product of the counting measure on $\mathcal{X}$ and the Lebesgue measure on $\mathbb{R}^1$ is
\[
\phi(z)=\phi_X(x)\phi_Y(y) \text{ for } z=(x,y) \in \mathbb{R}^2.
\]
We now describe our integral approximations. The first step of our procedure, raking, requires approximations of integrals with respect to one period variable. For example, an odd iteration of raking is
\begin{align*}
f^{(t+1)}(z_{1i},z_{2j}) = \frac{f_1(z_{1i}) f^{(t)}(z_{1i},z_{2j})}{\int f(z_{1i},z_2) \, d\mu(z_2) }.
\end{align*}
Here, we approximate the integral in the denominator as
\[
\int f(z_{1i},z_2) \, d\mu(z_2) \approx \frac{1}{S_2} \sum_{j=1}^{S_2} \frac{f(z_{1i},z_{2j})}{\phi_2(z_{2j})}.
\]
Similarly, for an even iteration, we approximate the denominator as
\[
\int f(z_1,z_{2j}) \, d\mu(z_1) \approx \frac{1}{S_1} \sum_{i=1}^{S_1} \frac{f(z_{1i},z_{2j})}{\phi_1(z_{1i})}.
\]
The second step of our procedure, calculating the moment conditions, requires the approximation of integrals with respect to the raked distribution.
If the model is implemented reliably in software packages (say, two-way fixed effects regression), we can draw a very large sample from the raking estimator $\hat f$ and treat it as the input to the existing routine.
This will then be numerically close to solving the moment conditions with plugged-in $\hat f$.
On the other hand, if calculating the moment conditions is necessary, we approximate $\operatorname{\mathbb{E}}_f \varphi(Z_1,Z_2)$ as
\begin{align*}
\iint \varphi(z_1,z_2) f(z_1,z_2) \, d\mu(z_1) d\mu(z_2) &\approx \frac{1}{S_1 S_2} \sum_{i=1}^{S_1} \sum_{j=1}^{S_2} \frac{\varphi(z_{1i},z_{2j}) f(z_{1i},z_{2j})}{\phi_1(z_{1i}) \phi_2(z_{2j})}.
\end{align*}
\section{Monte Carlo simulation}\label{sec:mc}
In this section, we illustrate the performance of our estimation procedure in a set of Monte Carlo simulations.
We employ two data-generating processes (DGPs) with discrete and continuous variables.
The discrete DGP is a Markov chain with states $\{0,\dots,d-1\}$, where $Z_1$ is distributed uniformly over the $d$ states ($d=5$ or $d=10$) and $Z_2$ is generated from $Z_1$ using a transition matrix with all the transition probabilities within the interval $[0.05, 0.95]$. From the marginal distribution of $Z_1$ and the conditional distribution of $Z_2|Z_1=z_1$, we calculate the marginal distribution of $Z_2$ that is used to draw the refreshment sample.
The parameter of interest is
\begin{align*}
\beta = \operatorname{\mathbb{E}} [Z_1 Z_2],
\end{align*}
and the selection probability is $e^{k_1(z_1)+k_2(z_2)}$ with
\begin{align*}
k_1(z_1) = -0.1(1+0.5 z_1), \quad k_2(z_2) = -0.2(1+0.5 z_2).
\end{align*}
The overall attrition rate is 40\% for $d=5$ and 53\% for $d=10$.
We assume that the support of the joint distribution of $Z_1, Z_2$ is known.
We use frequency estimators of the marginal probability mass functions of $Z_1$ and $Z_2$ and the selective joint probability mass function of $Z_1, Z_2$ as input to the raking procedure.
The continuous DGP models $(Z_1,Z_2)$ as a two-dimensional Gaussian vector with zero means, unit variances, and $\theta = Cov(Z_1,Z_2)=0.4.$
Again, the attrition process is nonignorable with selection probability $e^{k_1(z_1)+k_2(z_2)}$ and
\begin{align*}
k_1(z_1) = -0.1 |z_1|, \quad k_2(z_2) = -0.3 |z_2|.
\end{align*}
The overall attrition rate is about 30\%. The marginal distributions of $Z_1$ and $Z_2$ and the balanced panel distribution of $Z_1,Z_2$ are specified as normal with means and variances estimated by MLE. These estimated distributions are then an input into the raking estimator.
We report bias, standard deviation (sd), and root-MSE (rmse) of three estimators.
The infeasible estimator $\hat\theta_{\text{infeas}}$ uses the true values $k_1,k_2$ to weight the balanced panel sample. The naive estimator $\hat\theta_{\text{naive}}$ disregards attrition and treats the balanced panel as the full panel, and therefore is inconsistent. Neither of these two estimators uses the refreshment sample.
Finally, $\hat\theta$ is our procedure that uses raking to estimate the balanced panel weights.
\Cref{tab-mc-discr,tab-mc-cont} contain the simulation results.
We denote by $N$ the sum of the sample sizes of the incomplete panel and the refreshment sample, with the latter being $0.4 N$ for all specifications.
Our estimator performs well in terms of both bias and RMSE.
It typically has slightly higher bias than the infeasible estimator, but a lower RMSE, partially driven by its use of the refreshment sample.
\begin{table}[h!]
\centering
\small
\begin{tabular}{c c c c c c c c c c c}
\toprule
& & \multicolumn{3}{c}{$N=1000$} & \multicolumn{3}{c}{$N=5000$} & \multicolumn{3}{c}{$N=10000$}\\
\cmidrule(lr){3-5} \cmidrule(lr){6-8} \cmidrule(lr){9-11}
& & $\hat\theta_{\text{infeas}}$ & $\hat\theta_{\text{naive}}$ & $\hat\theta$
& $\hat\theta_{\text{infeas}}$ & $\hat\theta_{\text{naive}}$ & $\hat\theta$
& $\hat\theta_{\text{infeas}}$ & $\hat\theta_{\text{naive}}$ & $\hat\theta$ \\
\hline
\multirow{3}{*}{$d=5$}
& bias & 0.002 & -0.46 & -0.007 & 0.009 & -0.45 & -0.0008 & -0.0006 & -0.46 & -0.003 \\
& sd & 0.283 & 0.18 & 0.206 & 0.125 & 0.09 & 0.0920 & 0.0829 & 0.05 & 0.063 \\
& rmse & 0.283 & 0.50 & 0.207 & 0.125 & 0.46 & 0.0920 & 0.0829 & 0.46 & 0.064 \\
\midrule
\multirow{3}{*}{$d=10$}
& bias & -0.001 & -3.92 & -0.085 & 0.074 & -3.88 & -0.0095 & -0.0007 & -3.91 & -0.020 \\
& sd & 1.653 & 0.82 & 0.996 & 0.751 & 0.37 & 0.4169 & 0.5044 & 0.25 & 0.303\\
& rmse & 1.653 & 4.01 & 1.000 & 0.755 & 3.90 & 0.4171 & 0.5044 & 3.92 & 0.304\\
\bottomrule
\end{tabular}
\caption{Simulation results for the discrete DGP. Number of simulations $S=1000$.}
\label{tab-mc-discr}
\end{table}
\begin{table}[h!]
\centering
\small
\begin{tabular}{c c c c c c c}
\toprule
& \multicolumn{3}{c}{$N=1000$} & \multicolumn{3}{c}{$N=5000$} \\
\cmidrule(lr){2-4} \cmidrule(lr){5-7}
& $\hat\theta_{\text{infeas}}$ & $\hat\theta_{\text{naive}}$ & $\hat\theta$
& $\hat\theta_{\text{infeas}}$ & $\hat\theta_{\text{naive}}$ & $\hat\theta$ \\
\midrule
bias & -0.012 & -0.109 & -0.018 & -0.009 & -0.106 & 0.002 \\
sd & 0.062 & 0.041 & 0.054 & 0.032 & 0.022 & 0.029 \\
rmse & 0.063 & 0.116 & 0.057 & 0.034 & 0.109 & 0.029 \\
\bottomrule
\end{tabular}
\caption{Simulation results for the continuous DGP. Number of simulations $S=1000$.}
\label{tab-mc-cont}
\end{table}
\section{Empirical illustration}\label{sec:empirical}
\begin{table}[t]
\centering
\small
\begin{tabular}{ll}
\toprule
Population & Frequency \\
\midrule
Balanced Panel (BP) & $3303$ \\
Incomplete Panel (IP) & $904$ \\
Refreshment Sample (RS) & $3180$ \\
\bottomrule
\end{tabular}
\caption{ Sample sizes for the Balanced Panel, the Incomplete Panel, and the Refreshment Sample}
\label{tab:missing_data_pattern}
\end{table}
The empirical infrastructure that spawned the analysis in this paper is the Understanding America Study (UAS).\footnote{See \href{https://uasdata.usc.edu/index.php}{https://uasdata.usc.edu/index.php}}
The UAS is a probability-based Internet panel of about 15000 respondents.
Respondents without prior internet access receive a tablet and a broadband internet connection. Participants complete online surveys once or twice monthly.
Twenty-four core modules fielded biennially (about 400 total minutes) cover physical and mental health, economics, cognition, decision-making, and social determinants. Short monthly health updates capture acute and dynamic events. The system also supports high-frequency and event-triggered assessments (up to 6x/day) for ecological or crisis-response studies.
In addition, the UAS collects genetic information and some digital biomarkers.
Like any longitudinal study, the UAS suffers from attrition (7-8 percent per year), see \citet{kapteyn2024understanding}.
To counter the effect of attrition, but also because the panel is still growing, refreshment samples are drawn regularly (often several times a year). The current paper is a first step in a program to optimally use the refreshment samples for statistical inference.
The illustration in this paper is taken from the first two biennial waves of the UAS (May 20, 2015 - June 1, 2017 and June 1, 2017 - June 18, 2019). For expository reasons, we consider a very simple model of cognition. Most cognition dimensions tend to decrease with age. In addition, the literature has identified a large number of risk factors that increase the chance of dementia and generally may hasten cognitive decline. Here, we consider diabetes and depression as risk factors for cognitive decline. Education is generally found to be protective (\citet{livingston2024dementia}).
We use numeracy as a simple measure of cognition. Appendix \ref{app:variables} contains a description of the variables used in the empirical analysis. The numeracy scores are the result of answers to eight questions, and are calculated using a two-parameter logistic IRT (Item Response Theory) model. Scores are normalized to a mean of 50 and a standard deviation of 10. Higher scores indicate better performance at numeracy problems.
The first five columns of Table \ref{tab:primafacie} display the sample averages of the covariates over the different subsamples.\footnote{The panel and refreshment samples were obtained using stratified sampling. The analyses in this section account for this by using sampling weights.} In these columns, we can see that the Balanced Panel contains $75.8\%$ whites, while the refreshment sample contains $12.1\%$ fewer whites. The refreshment sample contains more depressed respondents and fewer respondents with a college degree, compared to the balanced panel. The dropouts in the incomplete panel suffer less from diabetes. Column (6) shows that dropouts have lower numeracy (as expected), and that numeracy is much lower in the balanced panel than in the refreshment sample. This holds even after correcting for covariates in column (7).
Although these differences are not conclusive, they indicate that attrition may lead to invalid inference when attrition is ignored by using only the balanced panel.
\begin{table}[]
\centering
\small
{
\def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}
\begin{tabular}{l*{7}{c}}
\toprule
&\multicolumn{1}{c}{(1)}&\multicolumn{1}{c}{(2)}&\multicolumn{1}{c}{(3)}&\multicolumn{1}{c}{(4)}&\multicolumn{1}{c}{(5)}&\multicolumn{1}{c}{(6)}&\multicolumn{1}{c}{(7)}\\
&\multicolumn{1}{c}{White}&\multicolumn{1}{c}{Depressed}&\multicolumn{1}{c}{College}&\multicolumn{1}{c}{Diabetes}&\multicolumn{1}{c}{Age}&\multicolumn{1}{c}{Numeracy}&\multicolumn{1}{c}{Numeracy}\\
\hline
IP & -0.009 & 0.034 & -0.016 & -0.034\sym{*} & -0.924 & -0.842\sym{*} & -0.676 \\
& (0.02) & (0.02) & (0.02) & (0.01) & (1.06) & (0.40) & (0.37) \\
RS & -0.121\sym{***}& 0.252\sym{***}& -0.105\sym{***}& -0.034\sym{***}& 0.352 & 1.239\sym{***}& 0.608\sym{*} \\
& (0.01) & (0.01) & (0.01) & (0.01) & (0.56) & (0.24) & (0.31) \\
Diabetes & & & & & & & -1.398\sym{***}\\
& & & & & & & (0.32) \\
College & & & & & & & 5.700\sym{***}\\
& & & & & & & (0.23) \\
Depressed & & & & & & & -2.138\sym{***}\\
& & & & & & & (0.25) \\
White & & & & & & & 4.047\sym{***}\\
& & & & & & & (0.25) \\
Constant & 0.758\sym{***}& 0.260\sym{***}& 0.534\sym{***}& 0.148\sym{***}& 53.995\sym{***}& 49.599\sym{***}& 44.252\sym{***}\\
& (0.01) & (0.01) & (0.01) & (0.01) & (0.36) & (0.14) & (0.28) \\
\bottomrule
\multicolumn{8}{l}{\footnotesize Standard errors in parentheses}\\
\multicolumn{8}{l}{\footnotesize \sym{*} \(p<0.05\), \sym{**} \(p<0.01\), \sym{***} \(p<0.001\)}\\
\end{tabular}
}
\caption{Differences in means between the Balanced Panel (BP, reference group), the Incomplete Panel (IP), and the Refreshment Sample (RS).}
\label{tab:primafacie}
\end{table}
We estimate a random-effects linear regression with random effect $u_i$ and idiosyncratic error $e_{it}$, using the standard random-effects error structure with variances $\sigma^2_u$ and $\sigma^2_e$. We allow the coefficient vector $\beta$ to be different in Period one and Period two.
Table \ref{tab:primafacie} showed that the distribution of the observed variables in the refreshment sample is substantially different from the distribution of these variables in the balanced panel. The additively nonignorable attrition model finds the population distribution that is consistent with the first-period marginal, obtainable from the balanced panel and the incomplete panel, and the second-period marginal, obtainable from the refreshment sample, using the raking weights.
The attrition correction can be expected to affect the conditional mean function of numeracy given the covariates differentially in the two periods.
Table \ref{tab:estimates} shows estimates and standard errors of this model when attrition is ignored (Missing Completely At Random, MCAR) as well as when corrected for possibly nonignorable attrition using the AN model. The estimates are obtained using a weighted GMM procedure using the raking weights. The standard errors are obtained using Theorem \ref{thm:asy-var}.
\begin{table}[]
\centering
\small{
\def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}
\begin{tabular}{@{}lcclcc@{}}
\toprule
& \multicolumn{2}{c}{(1)} & & \multicolumn{2}{c}{(2)} \\ \midrule
& \multicolumn{2}{c}{MCAR (Balanced Panel)} & & \multicolumn{2}{c}{AN Attrition} \\ \cmidrule(lr){2-3} \cmidrule(l){5-6}
& Numeracy, wave 1 & Numeracy, wave 2 & & Numeracy, wave 1 & Numeracy, wave 2 \\ \midrule
Age & -0.01502 & -0.00877 & & -0.0270 & -0.00914 \\
& (0.00829) & (0.00839) & & (0.0365) & (0.0360) \\
Diabetes & -0.74718 & -1.440\sym{***} & & -1.016 & -1.990 \\
& (0.522) & (0.521) & & (1.432) & (1.527) \\
College & 5.381\sym{***} & 5.084\sym{***} & & 5.699\sym{***} & 5.480\sym{***} \\
& (0.362) & (0.371) & & (1.527) & (1.594) \\
Depressed & -2.199\sym{***} & -2.139\sym{***} & & -1.606 & -0.484 \\
& (0.403) & (0.416) & & (1.521) & (1.534) \\
White & 4.118\sym{***} & 3.973\sym{***} & & 4.335\sym{***} & 3.574\sym{***} \\
& (0.417) & (0.427) & & (1.526) & (1.613) \\
Constant & 45.10\sym{***} & 45.04\sym{***} & & 46.41\sym{***} & 46.82\sym{***} \\
& (0.592) & (0.615) & & (2.417) & (2.514) \\ \midrule
$\sigma^{2}_{e}$ & \multicolumn{2}{c}{27.11\sym{***}} & & \multicolumn{2}{c}{27.10\sym{***}} \\
& \multicolumn{2}{c}{(1.092)} & & \multicolumn{2}{c}{(3.347)} \\
$\sigma^{2}_{u}$ & \multicolumn{2}{c}{42.22\sym{***}} & & \multicolumn{2}{c}{43.98\sym{***}} \\
& \multicolumn{2}{c}{(1.544)} & & \multicolumn{2}{c}{(5.27)} \\ \midrule
\multicolumn{6}{l}{\footnotesize Standard errors in parentheses} \\
\multicolumn{6}{l}{\footnotesize \sym{*} \(p<0.05\), \sym{**} \(p<0.01\), \sym{***} \(p<0.001\)} \\ \bottomrule
\end{tabular}
}
\caption{ \small{Estimation results of the model that ignores attrition (Missing Completely At Random, MCAR) and of the Additively Nonignorable (AN) attrition model.}}
\label{tab:estimates}
\end{table}
Note that the standard errors for the AN model are higher because it estimates the population model while taking attrition into account. In the AN model, only attending college and being white are significant covariates. Moreover, the regression coefficients of depressed and white in the two periods in the MCAR model are closer to each other compared to the AN model. The Wald test for the hypothesis that all regression coefficients are the same under MCAR gives an insignificant $\chi^2$ statistic of $3.73$. Under AN, this statistic equals $19.45$, leading to rejection of stability at conventional significance levels. This rejects a standard random effects analysis in this example.
Empirically, the most striking result in Table \ref{tab:estimates} is that under AN, both depressed and diabetes become insignificant, while they have a highly significant effect on numeracy under MCAR.
As Table \ref{tab:primafacie} makes clear, respondents in the balanced panel are less depressed, suffer more from diabetes, and are more likely to have a college degree, relative to respondents in the incomplete panel and refreshment samples. Unlike the MCAR estimates, the AN estimates are aligned with the marginal distributions of all variables in both time periods.
\section{Conclusion}\label{sec:conclusion}
Estimation and inference in panels with attrition and refreshment have long lacked a tractable implementation.
We fill this gap by showing that (a continuous version of) raking can be used to estimate the target distribution.
We derive a convergence rate of our two-step estimation strategy and provide a parametric approximation to the asymptotic variance.
Promising directions for future work include establishing bootstrap validity, deriving the semiparametric efficiency bound and a corresponding variance formula, extending the framework to multi-wave panels, and characterizing the identified set under weaker assumptions on the attrition mechanism.
\section*{Acknowledgements}
We are grateful to Tim Armstrong, Tim Christensen, Vasily Goncharenko, Lidia Kosenkova, Sergey Lototsky, Anna Mikusheva, Francesca Molinari, Elizaveta Rebrova, Azeem Shaikh, Andrei Zeleneev, and seminar participants at USC, UC Irvine, University of Virginia, and Penn State for valuable comments.
All errors and omissions are our own.
\bibliographystyle{ecta}
\bibliography{references}
\newpage
\part*{Appendix}