Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
55,054 characters · 9 sections · 38 citation commands
Raking for Estimation and Inference in Panel Models with Nonignorable Attrition and Refreshment
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 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 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 kish1959replacement. In panel surveys, these samples are drawn in the second and later waves and are called refreshment samples ridder1992empirical. deng2013handling and 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 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. bhattacharya2008inference studies the properties of this semiparametric estimator, and hoonhout2019nonignorable develops the estimator for panel data with three or more waves. hirano1998combining and 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. deng2013handling use multiple imputation inference rubin2004multiple. A disadvantage of the weighted moment and Bayesian estimators is that they are computationally demanding and require close monitoring by the statistician. bhattacharya2008inference reports that the moment estimator does not always converge.\footnote{See Franguridi2024 for further comments on 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 raking or iterative proportional fitting and was first proposed by deming1940least.
Starting from the observed joint distribution of the outcomes, often only a few steps are required for convergence to the unique solution.\footnote{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. (ref) introduces the framework. (ref) considers estimation. (ref) discusses inference. (ref) provides details for numerical implementation. (ref) reports a simulation study of our estimator. (ref) contains the empirical illustration. (ref) concludes.
We consider a two-period panel data model with attrition and refreshment as in 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 attrition: only the units with $W_i=1$ stay in the sample, while the remaining units drop out. To compensate for this attrition, a 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
where $\varphi$ is a $d_{\theta}$-dimensional moment function, and the expectation is taken with respect to the 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)) 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
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.
hirano2001combining impose the following separability condition on the weights.\footnote{franguridi2025inference drop this separability assumption, and franguridi2025set further consider the case where refreshment samples are not available. In both cases, the parameters of interest become partially identified.}
Throughout this paper, we further require the following.
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 (ref), the key identity (ref) becomes
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 hirano2001combining.
hirano2001combining show that Assumptions (ref) and (ref) 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 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 schrodinger1931umkehrung.}
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
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 (ref) 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
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)$ 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
we have
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 csiszar1975divergence, ireland1968contingency, and kullback1968probability.
The algorithm requires assumptions on the supports of the involved distributions.
(ref)(i) states that the variables in the data are either discrete or continuous. (ref)(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. (ref)(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, (ref)(iv) is satisfied if $f^w$ is bounded away from zero on its support.
Theorem 3.5 in ruschendorf1995convergence implies the following result.
The iterative equation of the theorem is the general data version of raking, also known as iterative proportional fitting or Sinkhorn's algorithm.\footnote{ Raking was introduced by deming1940least to adjust contingency tables to given marginals and was generalized to the continuous case by ireland1968contingency and kullback1968probability.} It is widely used in survey statistics and, since the seminal contribution of cuturi2013sinkhorn, has attracted substantial attention in machine learning due to its application to solving optimal transport problems with entropic regularization.
(ref) suggests the following two-step estimation strategy.
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.
(ref) is equivalent to imposing both (ref). (ref) is a standard assumption ensuring identification of $\theta$. (ref)(i) imposes a parametric model on the observed data. A sufficient condition for (ref)(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.
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 constrained program (ref) to the dual characterization (ref) with $k_1,k_2$ solving an unconstrained program. This duality result is not new (see, e.g., Section 4 in bhattacharya2006iterative and Theorem 2.3 in bhattacharya1995general), but we provide a proof to keep the presentation self-contained.
Finally, we establish $\sqrt{n}$-consistency of the raking estimator of $f$ and the estimator of $\theta$.
The asymptotic distribution of the moment estimator of $\theta$ depends on the asymptotic distribution of the estimator ((ref)) of the pseudo-parameter $f$ defined in ((ref)). That estimator is computed as in (ref). 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)) is computed by a finite number of iterations as in (ref). The asymptotic distribution of the estimator ((ref)) is therefore equal to $T$-fold recursion in (ref).
Under (ref), 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.
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
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
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
and the selection probability is $e^{k_1(z_1)+k_2(z_2)}$ with
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
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.
(ref) 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.
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 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 (livingston2024dementia).
We use numeracy as a simple measure of cognition. Appendix (ref) 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) 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.
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) 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) 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).
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) is that under AN, both depressed and diabetes become insignificant, while they have a highly significant effect on numeracy under MCAR. As Table (ref) 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.
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.
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.
\part*{Appendix}