Leveraging Random Assignment to Impute Missing Covariates in Causal Studies

TL;DR

This paper examines whether leveraging randomized treatment assignment improves the accuracy of imputing missing covariates in experimental studies, finding only small gains in typical scenarios.

Contribution

It investigates the benefits of using randomization information in imputation methods for missing covariates in randomized experiments, through comprehensive simulation studies.

Findings

Respecting randomization yields small accuracy improvements.

Imputation methods using covariates or outcomes show similar performance.

Guidance on imputing missing covariates for heterogeneous treatment effect estimation.

Abstract

Baseline covariates in randomized experiments are often used in the estimation of treatment effects, for example, when estimating treatment effects within covariate-defined subgroups. In practice, however, covariate values may be missing for some data subjects. To handle missing values, analysts can use imputation methods to create completed datasets, from which they can estimate treatment effects. Common imputation methods include mean imputation, single imputation via regression, and multiple imputation. For each of these methods, we investigate the benefits of leveraging randomized treatment assignment in the imputation routines, that is, making use of the fact that the true covariate distributions are the same across treatment arms. We do so using simulation studies that compare the quality of inferences when we respect or disregard the randomization. We consider this question for imputation routines implemented using covariates only, and imputation routines implemented using the outcome variable. In either case, accounting for randomization offers only small gains in accuracy for our simulation scenarios. Our results also shed light on the performances of these different procedures for imputing missing covariates in randomized experiments when one seeks to estimate heterogeneous treatment effects.

Tables

Table 1. Table 1: Counts for the observed data generated as per ICIN using the loglinear model approach, for a fixed random seed. NA represents a missing value.

	$X_2=0$	$X_2=1$	$X_2=\text{NA}$
$X_1=0$	103	55	65
$X_1=1$	97	133	197
$X_1=\text{NA}$	68	155	127

Table 2. Table 2: Absolute biases, MC standard deviations, estimated standard errors, coverage probabilities, and average CI lengths for β t subscript 𝛽 𝑡 \beta_{t} and β t ​ x 2 subscript 𝛽 𝑡 subscript 𝑥 2 \beta_{tx_{2}} estimates under the high association setting in scenario 1. True values of β t = 0.3 subscript 𝛽 𝑡 0.3 \beta_{t}=0.3 and β t ​ x 2 = 0.5 subscript 𝛽 𝑡 subscript 𝑥 2 0.5 \beta_{tx_{2}}=0.5 .

		${\beta }_t$						${\beta }_{t{x}_2}$
Method		Absolute bias	MC-SD	SE	Coverage	Average CI length		Absolute bias	MC-SD	SE	Coverage	Average CI length
MI-R		0.039	0.192	0.225	0.974	0.881		0.192	0.255	0.370	0.978	1.449
MI-NR		0.040	0.192	0.225	0.975	0.884		0.194	0.255	0.371	0.983	1.453
MI-RY		0.000	0.232	0.228	0.950	0.890		0.021	0.433	0.427	0.946	1.668
MI-NRY		0.000	0.233	0.230	0.957	0.903		0.022	0.436	0.433	0.947	1.696
Mean-R		0.031	0.241	0.239	0.948	0.938		0.018	0.408	0.408	0.949	1.601
Mean-NR		0.030	0.249	0.246	0.949	0.940		0.017	0.415	0.412	0.947	1.606
Mean-NRY		0.078	0.311	0.234	0.850	0.919		0.310	0.616	0.438	0.790	1.715
Reg-R		0.040	0.210	0.207	0.947	0.813		0.199	0.316	0.314	0.902	1.231
Reg-NR		0.038	0.214	0.210	0.941	0.814		0.188	0.322	0.315	0.900	1.234
Reg-NRY		0.001	0.252	0.203	0.884	0.798		0.023	0.489	0.335	0.814	1.313
CCA		0.006	0.310	0.296	0.940	1.164		0.062	0.549	0.556	0.948	2.188

Table 3. Table 3: Absolute biases, MC standard deviations, estimated standard errors, coverage probabilities, and average confidence lengths for β t subscript 𝛽 𝑡 \beta_{t} and β t ​ x 2 subscript 𝛽 𝑡 subscript 𝑥 2 \beta_{tx_{2}} estimates under the low association setting in scenario 1. True values of β t = 0.3 subscript 𝛽 𝑡 0.3 \beta_{t}=0.3 and β t ​ x 2 = 0.015 subscript 𝛽 𝑡 subscript 𝑥 2 0.015 \beta_{tx_{2}}=0.015 .

		${\beta }_t$						${\beta }_{t{x}_2}$
Method		Absolute bias	MC-SD	SE	Coverage	Average CI length		Absolute bias	MC-SD	SE	Coverage	Average CI length
MI-R		0.002	0.168	0.208	0.989	0.816		0.002	0.208	0.302	0.995	1.184
MI-NR		0.002	0.169	0.209	0.988	0.818		0.002	0.209	0.303	0.995	1.187
MI-RY		0.005	0.222	0.220	0.948	0.858		0.012	0.334	0.327	0.948	1.279
MI-NRY		0.005	0.225	0.224	0.945	0.871		0.010	0.340	0.331	0.944	1.309
Mean-R		0.007	0.224	0.219	0.945	0.861		0.011	0.340	0.332	0.942	1.303
Mean-NR		0.007	0.233	0.225	0.945	0.861		0.011	0.345	0.339	0.943	1.303
Mean-NRY		0.014	0.332	0.220	0.809	0.864		0.028	0.569	0.332	0.782	1.307
Reg-R		0.000	0.191	0.189	0.950	0.743		0.003	0.268	0.257	0.943	1.009
Reg-NR		0.003	0.197	0.194	0.946	0.745		0.003	0.276	0.267	0.945	1.012
Reg-NRY		0.007	0.251	0.190	0.866	0.744		0.012	0.399	0.258	0.797	1.013
CCA		0.001	0.286	0.279	0.947	1.094		0.019	0.419	0.416	0.940	1.631

Table 4. Table 4: Point estimates, standard errors, and 95 % confidence intervals for the log-odds scale average treatment effect in the randomized experiment from Foos and Gilardi ( 2019 ) .

Method	Estimate	SE	95% CI
MI-R (ICIN)	-0.462	0.678	(-1.790,0.866)
MI-NR (ICIN)	-0.462	0.678	(-1.791,0.867)
MI-RY (ICIN)	-0.466	0.669	(-1.794,0.831)
MI-NRY (ICIN)	-0.469	0.678	(-1.798,0.860)
MI-R (MAR)	-0.464	0.678	(-1.792,0.864)
MI-NR (MAR)	-0.464	0.678	(-1.792,0.865)
MI-RY (MAR)	-0.451	0.672	(-1.792,0.827)
MI-NRY (MAR)	-0.467	0.678	(-1.795,0.861)
Mean-R	-0.463	0.677	(-1.805,0.944)
Mean-NR	-0.463	0.679	(-1.805,0.944)
Mean-NRY	-0.465	0.677	(-1.806,0.943)
Reg-R	-0.464	0.678	(-1.806,0.944)
Reg-NR	-0.468	0.679	(-1.810,0.940)
Reg-NRY	-0.462	0.678	(-1.804,0.945)
CCA	-0.692	0.715	(-2.090,0.709)

Equations

$$\text{logit }(Pr(Y=1|X_1,X_2,T))={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_{t}T+{\beta }_{t{x}_2}T{X}_2$$

$${\theta }_1$$

$${\theta }_2$$

$${\theta }_3$$

$${\theta }_4$$

$${\theta }_5$$

$${\theta }_6$$

$${\theta }_7$$

$${\theta }_8$$

$${\theta }_9$$

$$f(X_{mis},X_{obs},D_1,D_2,T,Y)$$

$$f(X_{mis}|X_{obs},D_1,D_2,T,Y)$$

$$\propto f(X_{mis},X_{obs},D_1,D_2)f(Y|X_{mis},X_{obs},D_1,D_2,T).$$

$$f(X_{mis}|X_{obs},D_1,D_2,T,Y,\gamma )\propto f(X_{mis},X_{obs},D_1,D_2,\theta )f(Y|X_{obs},X_{mis},D_1,D_2,T,\beta ).$$

$$\text{log }{m}_{(x_1,x_2).(d_1,d_2)}$$

$$+0.5x_1{x}_2+0.75x_1{d}_2+1x_2{d}_1+0.25d_1{d}_2.$$

$${\pi }_{(x_1,x_2).(d_1,d_2)}=\frac{m_{(x_1,x_2).(d_1,d_2)}}{{\sum }_{x_1}{\sum }_{x_2}{m}_{(x_1,x_2).(d_1,d_2)}}.$$

$$\text{logit }(Pr(Y=1|X_1,X_2,T))={\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_{t}T+{\beta }_{t{x}_2}T{X}_2.$$

$$\text{logit }(Pr(Y=1|X_1,X_2,D_1,D_2,T))={\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_{t}T+{\beta }_{t{x}_2}T{X}_2-0.6D_1-0.4D_2,$$

Full text

Keywords: experiment; imputation; non-ignorable; randomization.

1 Introduction

Randomized experiments are widely considered to be the gold standard for causal inference. Their appeal, in large part, is attributed to randomized treatment assignment, which ensures baseline comparability of all treatment groups. Alternatively stated, randomization balances all covariates on average (Rubin, 2008), facilitating simple causal comparisons. Randomized studies can be broadly separated into two stages, namely (1) the stage where the experiment is planned, covariates measured, and treatments assigned, which Rubin (2007) calls the “design stage”, and (2) the stage where the outcome variable is measured and compared across treatment groups, which we henceforth call the “outcome stage”. Since the design stage is executed without any outcome data in view, results do not systematically favor particular treatment groups (Rubin, 2008).

As with all studies, the data from randomized experiments can suffer from missing values, both in the outcomes and the covariates. While much research focuses on missing outcomes (e.g., Frangakis and Rubin, 1999; Chen et al., 2009; Imai, 2009), we focus on missing covariates here, assuming for didactic reasons that no outcome values are missing. In this situation, analysts who estimate average treatment effects with simple comparisons of outcome means can disregard the missingness in the covariates. However, analysts generally cannot do so when they estimate subgroup treatment effects or use regression adjustment. In such cases, a complete case analysis (CCA) sacrifices information and can result in distorted estimates. Moreover, CCA violates the intention to treat (ITT) principle in randomized studies (White and Thompson, 2005).

An alternative is to impute values for the missing observations. For example, in mean imputation, the analyst replaces each missing value in the covariate $X$ with the mean of the observed values of $X$. In stochastic regression imputation, the analyst replaces each missing value in $X$ with a draw from a regression model relating $X$ to other covariates and possibly the treatment or the outcome. These single imputation techniques for missing covariates have been described in detail by Schemper and Smith (1990) and White and Thompson (2005). In multiple imputation (MI) (Rubin, 1987, 1996; Little and Rubin, 2002), the analyst creates multiple completed datasets by filling in missing values with draws from predictive distributions estimated with the observed data. The analyst performs complete-data causal inference on each completed dataset, and combines the results using the inferential methods of Rubin (1987).

When implementing any of these techniques in randomized experiments, analysts can utilize the randomization in the imputation routine. In MI, for example, they can require the imputed data to be samples from the same covariate distribution, regardless of treatment assignment. This may allow analysts to estimate parameters of the MI model more accurately than estimating separate imputation models in each treatment arm. In turn, this could reduce variances in the imputations and increase accuracy in estimating causal effects from the completed datasets. As another example, in mean imputation, the analyst could pool the observed data across treatment arms to estimate the covariate means, as opposed to estimating covariate means separately in each treatment arm.

In this article, we investigate the magnitude of such improvements in moderately-sized randomized experiments for MI, regression imputation, and mean imputation. For all methods, we consider imputation routines that adhere to the tenets of Rubin (2007, 2008) by keeping the design and outcome stages separate; that is, we do not include the outcome in the imputation routine. We also consider methods that use the outcome in the imputation routine, as recommended by other authors (e.g., Vach, 1994; Moons et al., 2006). We assess the benefits of respecting or not respecting the randomization, as well as the benefits of using or not using the outcome variable, via simulation studies based on an ignorable and a non-ignorable missing data scenario. We also use the simulation results to discuss the performances of the single imputation and the MI approaches more generally.

Our results contribute to previous work on handling missing covariates in randomized experiments (e.g., White and Thompson, 2005; Groenwold et al., 2012; Jolani and Safarkhani, 2017; Sullivan et al., 2018; Kayembe et al., 2020). In particular, as part of broader simulation studies, White and Thompson (2005) compare properties of adjusted treatment effects when performing mean imputation within and across treatment arms. Sullivan et al. (2018) compare properties of causal estimates when performing MI with a treatment indicator included in the imputation model, and separately within treatment arms. These works use methods and run simulations for a single covariate, whereas we use more than one covariate. They focus on estimating outcome models assuming constant treatment effects, whereas we estimate outcome models where treatment effects are modified by covariate values. Finally, they use simulations with smaller sample sizes than we do. As a result of these differences, we are able to reach additional conclusions about the benefits of leveraging randomization, and we contribute new insights into the performances of the different imputation methods.

As part of the simulations, we develop and illustrate an MI approach for covariate values that ensures non-parametric identification (Vansteelandt et al., 2006; Robins, 1997). This approach is particularly suited for covariates that are not missing at random. We describe this approach in Section 3, as well as in detail in the supplementary materials.

The remainder of this article is organized as follows. In Section 2, we provide background on randomized experiments and on missing data modeling. In Section 3, we describe how we respect random assignment in the context of MI. In Section 4, we describe the simulation design and results. In Section 5, we implement all the methods on data from Foos and Gilardi (2019). Finally, in Section 6, we conclude with a discussion.

2 Background

We begin with notation and key assumptions in randomized experiments. Here and throughout, we assume a parallel design with an active treatment and a control.

2.1 Randomized experiments

We consider a randomized experiment with $n$ units. For $i=1,\dots,n$, let $T_{i}=1$ when unit $i$ is assigned to the treatment, and $T_{i}=0$ when unit $i$ is assigned to the control. For $i=1,\dots,n$, let $X_{ij}$ denote the value of covariate $j$ measured for unit $i$, and let $X_{i}=(X_{i1},\dots,X_{ip})$ represent measurements on $p$ covariates of interest. Let $Y_{i}(1)$ and $Y_{i}(0)$ be the potential outcomes associated with the treatment and the control, respectively. We assume that all $n$ units’ potential outcomes, $(Y_{i}(1),Y_{i}(0))$, are independent. For unit $i$, let $Y_{i}$ be the outcome value observed at the end of the experiment.

We derive results under the stable unit treatment value assumption (SUTVA) (Rubin, 1978). Specifically, we assume that there is only one version of the treatment and no interference between units. With SUTVA, we can express $Y_{i}$ as a deterministic function of $Y_{i}(1)$, $Y_{i}(0)$, and $T_{i}$, namely $Y_{i}=T_{i}Y_{i}(1)+(1-T_{i})Y_{i}(0)$. For a completely randomized experiment under SUTVA, the treatment status $T_{i}$ is independent of the potential outcomes $(Y_{i}(1),Y_{i}(0))$ and the set of covariates $X_{i}$. We write this as T_{i}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}\{Y_{i}(1),Y_{i}(0),X_{i}\}, where $i=1,\dots,n$.

2.2 Missing data modeling

When covariate values are non-ignorably missing (Rubin, 1976), the analyst needs to model the covariates and their missingness indicators. This joint distribution cannot be identified without recourse to generally untestable assumptions (Imbens and Pizer, 2000; Ding and Geng, 2014), also known as identifying restrictions. In this section, we describe the approach that we use to make such assumptions in MI, following extant work (e.g., Linero and Daniels, 2018).

For ease of exposition, we momentarily ignore the treatment and the outcome, and consider $p$ generic variables measured on $n$ subjects, i.e., $X_{i}=(X_{i1},\dots,X_{ip})$ where $i=1,\dots,n$. We let $D_{i}=(D_{i1},\dots,D_{ip})$ denote their missingness indicators, such that $D_{ij}=1$ when $X_{ij}$ is missing, and $D_{ij}=0$ otherwise, for $i=1,\dots,n$ and $j=1,\dots,p$. In what follows, we forgo the subscript $i$ indexing subjects to simplify notation. For example, $X=(X_{1},\dots,X_{p})$ represents a generic vector of $p$ variables, and $D=(D_{1},\dots,D_{p})$ represents a generic missingness pattern, i.e., the vector of values for the $p$ missingness indicators.

Based on $D$, we can divide $X$ into $X_{obs}$, the observed part of $X$, and $X_{mis}$, the missing part of $X$. The full data distribution, $f(X_{mis},X_{obs},D)$, can be factored into the product of the observed data distribution $f(X_{obs},D)$, and the extrapolation distribution $f(X_{mis}|X_{obs},D)$ (Daniels and Hogan, 2008). Here, $f(X_{obs},D)$ is identifiable from the observed data, whereas $f(X_{mis}|X_{obs},D)$ is not. Thus, we need to construct $f(X_{mis}|X_{obs},D)$, and consequentially $f(X_{mis},X_{obs},D)$, by imposing identifying restrictions. One reasonable desideratum for specifying $f(X_{mis},X_{obs},D)$, when feasible, is to ensure non-parametric identification (Vansteelandt et al., 2006), defined below.

Definition 1** (Non-parametric Identification).**

Let $\mathscr{F}$ denote the class of observed data distributions, and let $\mathscr{F_{A}}$ denote a family of full data distributions under restriction(s) $\mathscr{A}$. $\mathscr{F_{A}}$ is said to be non-parametric identified if given $f\in$ $\mathscr{F}$, there exists a unique $f_{A}\in$ $\mathscr{F_{A}}$ that marginalizes to $f$.

Several restrictions proposed in the literature fulfill this desideratum (see Linero and Daniels, 2018, and references therein). In this article, we confine our attention to two, namely the itemwise conditionally independent non-response (ICIN) assumption (Sadinle and Reiter, 2017) and the missing at random (MAR) assumption (Gill et al., 1997).

Definition 2** (Itemwise Conditionally Independent Non-response (ICIN)).**

$(X_{1},\dots,X_{p})$ are missing according to the itemwise conditionally independent non-response assumption when X_{j}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}D_{j}|X_{-j},D_{-j}, for $j=1,\dots,p$. Here, $X_{-j}=(X_{1},\dots,X_{j-1},X_{j+1},\dots,X_{p})$, and $D_{-j}$ is defined likewise.

In other words, with ICIN we presume that controlling for all other variables and their missingness indicators, the missingness of $X_{j}$ does not predict its value and vice versa.

Definition 3** (Missing at Random (MAR)).**

If $(X_{1},\dots,X_{p})$ are missing at random, then for each missingness pattern $D$, $f(D|X_{mis},X_{obs})=f(D|X_{obs})$. This can be re-expressed as $f(X_{mis}|X_{obs},D)=f(X_{mis}|X_{obs})$.

With MAR, for every $D$, we assume that the missingness is independent of the missing data conditional on the observed data. In the supplementary materials, we illustrate how ICIN and MAR result in identifiable joint distributions for $p=2$ variables.

For MI, we fill in missing values $M>1$ times, by drawing independently from the extrapolation distribution $f(X_{mis}|X_{obs},D)$. This creates $M$ completed data sets, each of which is analyzed separately. Point and variance estimates are combined using the rules in Rubin (1987).

3 MI methods for handling missing covariates

For purposes of illustration, we consider a randomized experiment with two binary covariates, $(X_{1},X_{2})\in\{0,1\}^{2}$, having associated missingness indicators $(D_{1},D_{2})$ $\in\{0,1\}^{2}$. The treatment status $T\in\{0,1\}$, and the binary outcome $Y\in\{0,1\}$, are fully observed. For estimating treatment effects, we focus on analysis models with effect modification. In particular, we assume the logistic regression,

[TABLE]

is of interest. The parameters of importance are $\beta_{t}$ and $\beta_{tx_{2}}$.

3.1 Design stage modeling

Covariates are design stage quantities. To follow the principle that the design stage should remain free of any influence from the outcome, the imputation model for $(X_{1},X_{2})$ should only include quantities observable at the time of treatment assignment. Hence, in our example, the data used for imputation modeling will be $(X_{1},X_{2},D_{1},D_{2},T)$. This can be expressed as a $2^{5}$ contingency table.

Covariates should be measured before randomizing units to treatment or control. When this happens, it is generally the case that the reasons for missingness in the covariates do not depend on treatment group membership. As a result, it is reasonable to regard $(D_{1},D_{2})$ as pretreatment variables, implying that (D_{1},D_{2})\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}T. Since we have (X_{1},X_{2})\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}T, we can collapse the $2^{5}$ table over $T$, forming a $2^{4}$ marginal table for $(X_{1},X_{2},D_{1},D_{2})$. Thus, in the design stage, respecting randomization is equivalent to treating $(X_{1},X_{2},D_{1},D_{2})$ as full data. MI then proceeds by identifying $f(X_{1},X_{2},D_{1},D_{2})$, and drawing imputations from $f(X_{mis}|X_{obs},D_{1},D_{2})$.

Under this method, let $\theta$ denote the vector of probabilities offered by the observed data. We have $\theta=(\theta_{1},\dots,\theta_{9})$, where,

[TABLE]

For MI inference, we treat the corresponding counts as a multinomial sample, and place a Dirichlet$(\frac{1}{9},\dots,\frac{1}{9})$ prior on $\theta$. We further generate an imputation in three steps. We sample a value for $\theta$ from its posterior distribution, which is also Dirichlet. Using this value, we obtain extrapolation distributions for all missingness patterns under the identifying restriction of choice. We derive these distributions for the ICIN and MAR assumptions in the supplementary materials. For each missing value, we then obtain an imputation from the pattern-specific extrapolation distribution. We repeat these steps $M$ independent times, creating $M$ completed datasets.

Alternatively, not utilizing randomization means disregarding the independence of $(X_{1},X_{2})$ and $(D_{1},D_{2})$ with respect to $T$. We hence use the full $2^{5}$ contingency table, and identify $f(X_{1},X_{2},D_{1},D_{2},T)$. Imputations are generated from $f(X_{mis}|X_{obs},D_{1},D_{2},T)$, analogous to separately imputing within the $T=0$ and $T=1$ groups. Under this method, $\theta$ is a vector of eighteen observed probabilities, given by (2) - (10) in each treatment arm. Multiple imputation then proceeds in the aforementioned manner, except that we use a Dirichlet$(\frac{1}{18},\dots,\frac{1}{18})$ prior for $\theta$.

Intuitively, accounting for randomization can offer the potential for improved accuracy in estimating treatment effects. By collapsing over treatment groups, we estimate imputation model parameters using the full study sample. In contrast, by imputing separately within treatment groups, we estimate imputation model parameters in each group, using a smaller sample size. This decreased sample size can result in larger parameter uncertainty, which in turn can result in greater variability in the imputations, and hence the MI inferences.

3.2 Outcome stage modeling

For missing covariates in regression models, it has been generally recommended that the outcome $Y$ be used in MI (Rubin and Schenker, 1991; Vach and Blettner, 1991; Greenland and Finkle, 1995; Barnard and Meng, 1999; Little and Rubin, 2002; Moons et al., 2006; Sterne et al., 2009). In randomized experiments, this has the disadvantage of allowing $Y$ to directly influence the design stage, which can bring the face validity of the final conclusions into question. At the same time, not controlling for $Y$ in imputations can lead to distorted estimates, particularly when regression-adjusted estimators are of interest. Little (1992) elucidates this issue: if a partly missing covariate $X$ is highly predictive of $Y$, then $Y$ will carry information about $X$ that may not be captured by other variables in the imputation model. If $X$ is imputed without using $Y$, then the imputed part of $X$ will have no (conditional) association with $Y$. This could falsely attenuate the overall covariate-outcome association.

In our example, an imputation model in the outcome stage uses $(X_{1},X_{2},D_{1},D_{2},T,Y)$ as data. This forms a contingency table with $2^{6}$ cells. One way to allow randomization to play a role here is to collapse this $2^{6}$ table across $T$. However, this makes a strong assumption that Y\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}T, which ultimately could underestimate the treatment effect (see Lu and Ashmead, 2018, for an illustrative simulation). A more principled way is to factorize the joint distribution $f(X_{mis},X_{obs},D_{1},D_{2},T,Y)$ to naturally represent the design and outcome stages as

[TABLE]

Here, $f(X_{mis},X_{obs},D_{1},D_{2},T)$ represents the joint distribution of the design stage quantities, and $f(Y|X_{mis},X_{obs},D_{1},D_{2},T)$ is the entire outcome response surface. Under random assignment, we have (X_{mis},X_{obs})\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}T, and assume (D_{1},D_{2})\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}T, so that we have

[TABLE]

We can generate imputations from (12) using the data augmentation strategy introduced by Tanner and Wong (1987). For a given draw of the parameter vector $\gamma=(\theta,\beta)$, where $\theta=(\theta_{1},\dots,\theta_{9})$ and $\beta=(\beta_{0},\beta_{1},\beta_{2},\beta_{t},\beta_{tx_{2}})$, we generate an imputation using

[TABLE]

Given the imputed dataset, we update $\theta$ and $\beta$ from their full conditional posterior distributions. For $\theta$, this is a Dirichlet distribution. For $\beta$, the conditional posterior is non-standard; we sample from this distribution using the Polya-Gamma latent variable technique outlined in Polson et al. (2013).

In practice, we specify $f(Y|X_{mis},X_{obs},D_{1},D_{2},T)$ using the outcome model posited for analysis. Often, however, analysis models of interest do not adjust for the missingness indicators, in which case, the conditional independence assumption Y\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}(D_{1},D_{2})|X_{1},X_{2},T is implicitly made.

If we ignore randomization, we essentially use the treatment as well as the outcome in the imputation model. This amounts to identifying $f(X_{1},X_{2},D_{1},D_{2},T,Y)$, and generating imputations from $f(X_{mis}|X_{obs},D_{1},D_{2},T,Y)$. Under this method, $\theta$ is a vector of thirty-six observed probabilities, given by (2) - (10) in each category of $(Y,T)$. We use a Dirichlet$(\frac{1}{36},\dots,\frac{1}{36})$ prior for $\theta$, and carry out MI as described in Section 3.1.

4 Simulation study

We now conduct repeated sampling studies to evaluate the performance of the imputation methods. We use twelve simulation settings, comprising all combinations of three missingness scenarios, two identifying restrictions, and two covariate-outcome associations. We begin by describing the data generation process.

4.1 Data generation

We simulate the randomized experiment from our running example for $n=1000$ units. Each unit is randomized to the treatment or the control arm with equal probabilities, i.e., $T$ is generated from a Bernoulli distribution with mean 0.5. We generate $(X_{1},X_{2})$ and $(D_{1},D_{2})$ as per the MAR and ICIN assumptions, under three missingness scenarios.

4.1.1 Scenario 1

In scenario 1, we regard $(D_{1},D_{2})$ as pre-treatment variables, and generate them independently of $T$.

To simulate an MAR situation under this scenario, we first draw $X_{1}\sim\text{Bernoulli }(0.7)$, followed by $X_{2}|X_{1}=1\sim\text{Bernoulli }(0.6)$ and $X_{2}|X_{1}=0\sim\text{Bernoulli }(0.45)$. Subsequently, we delete at random 35% of the observations from $X_{1}$, and 40% of the observations from $X_{2}$. The true missingness mechanism is thus missing completely at random (MCAR), which is a special case of MAR.

To create missingness as per ICIN, we make use of the relationship between the ICIN assumption and hierarchical loglinear models. A discussion of this relationship and the data generation approach that follows is deferred to the supplementary materials. We begin by generating $D_{1}\sim\text{Bernoulli }(0.35)$ and $D_{2}\sim\text{Bernoulli }(0.40)$. Next, we split the $2^{4}$ contingency table for $(X_{1},X_{2},D_{1},D_{2})$ into four, partial $2^{2}$ tables, controlling for $(D_{1},D_{2})$. Let $m_{(x_{1},x_{2}).(d_{1},d_{2})}$ represent the expected count for $(X_{1},X_{2})=(x_{1},x_{2})$ in the partial table with $(D_{1},D_{2})=(d_{1},d_{2})$. We simulate $m_{(x_{1},x_{2}).(d_{1},d_{2})}$ from the loglinear model

[TABLE]

Further, we obtain multinomial probabilities $\pi_{(x_{1},x_{2}).(d_{1},d_{2})}$ using (Agresti, 2012)

[TABLE]

For each missingness pattern $(D_{1},D_{2})=(d_{1},d_{2})$, we jointly draw $(X_{1},X_{2})=(x_{1},x_{2})$ with probability $\pi_{(x_{1},x_{2}).(d_{1},d_{2})}$, and set $X_{j}$ to missing wherever $D_{j}=1$, where $j=1,2$. Similar to the MAR setting, this approach (asymptotically) gives $P(X_{1}=1)\approx 0.7$, $P(X_{2}=1|X_{1}=0)\approx 0.45$, and $P(X_{2}=1|X_{1}=1)\approx 0.6$. To illustrate, we display counts for one random draw from the ICIN data generation procedure in Table 1.

We generate $Y$ from a Bernoulli distribution with probabilities defined by

[TABLE]

We set $\beta_{t}=0.3$. We vary the strength of the association between $(X_{1},X_{2})$ and $Y$ in two settings. In the high association setting, we set $\beta_{1}=0.8,\beta_{2}=0.9$, and $\beta_{tx_{2}}=0.5$. In the low association setting, we set $\beta_{1}=0.02$, $\beta_{2}=0.05$, and $\beta_{tx_{2}}=0.015$.

4.1.2 Scenario 2

In scenario 2, we regard $(D_{1},D_{2})$ to be post-treatment variables. Here, data are generated using the models for $(X_{1},X_{2},T,Y)$ from scenario 1, except that we draw $D_{1}$ and $D_{2}$ from the conditional distributions $D_{1}|T=1\sim\text{Bernoulli }(0.35)$, $D_{1}|T=0\sim\text{Bernoulli }(0.1)$, and $D_{2}|T=1\sim\text{Bernoulli }(0.40)$, $D_{2}|T=0\sim\text{Bernoulli }(0.1)$. Thus, missingness rates differ by treatment arms, but do not depend on the covariate values.

4.1.3 Scenario 3

In scenario 3, we allow $(D_{1},D_{2})$ to be predictive of $Y$. We generate $(X_{1},X_{2})$ and $(D_{1},D_{2})$ as in scenario 1. We linearly adjust for $(D_{1},D_{2})$ in the outcome generation model, i.e., we use

[TABLE]

where $(\beta_{1},\beta_{2},\beta_{t},\beta_{tx_{2}})$ are defined as in scenario 1.

4.2 Methods

For MI, we create $M=100$ datasets as per the four approaches described in Section 3, namely respecting randomization in the design stage (MI-R), not respecting randomization in the design stage (MI-NR), respecting randomization in the outcome stage (MI-RY), and not respecting randomization in the outcome stage (MI-NRY). For each imputed dataset, we obtain point and variance estimates for $\beta_{t}$ and $\beta_{tx_{2}}$, and combine them for MI inferences.

We also examine the performance of mean imputation, stochastic regression imputation, and complete case analysis (CCA). For mean imputation, we replace the missing values for each covariate by the mean of the observed values of that covariate. Although not typically used for categorical data, mean imputation for dichotomous covariates has been validated in Schemper and Smith (1990) and Sullivan et al. (2018) for estimation of average treatment effects. Specifically, we consider mean imputation across treatment arms (Mean-R), in which we use the covariate means from the combined sample of the treated and control cases; mean imputation within treatment arms (Mean-NR), in which we use the covariate means for the observed treated cases in the treatment arm and the observed control cases in the control arm; and, as a way to use the outcome values, mean imputation within each of the four cells in the cross-classification of the treatment and the binary outcome (Mean-NRY). For regression imputation, we replace missing values for a covariate by draws from the posterior predictive distribution obtained from a logistic regression on the other covariate, where each model is estimated using only the complete cases. In particular, we evaluate regression imputation across treatment arms (Reg-R), in which the model is estimated on the combined sample and independently of the treatment; regression imputation within treatment arms (Reg-NR), in which separate regression models are estimated within the treatment and control arms; and, regression imputation within within each of the four cells in the cross-classification of the treatment and the binary outcome (Reg-NRY).

In each scenario, we generate 1000 independently sampled, simulated datasets. For each method, we compute the Monte Carlo estimate of the absolute bias in the treatment effect estimates, the Monte Carlo estimate of the standard deviation (MC-SD) of the treatment effect estimates, the square root of the average estimated variance (i.e., the estimated standard error, SE), the coverage rate of the 95% confidence intervals (CI), and the average length of the 95% CI.

4.3 Results

Here, we present results under the ICIN assumption. Results under the MAR assumption are qualitatively similar and are included in the supplementary materials.

4.3.1 Scenario 1: High association

Figure 1 displays the coefficient estimates for $\beta_{t}$ and $\beta_{tx_{2}}$ under scenario 1, when the covariate-outcome association is high. Table 2 displays the absolute biases, MC-SDs, estimated SEs, coverage rates, and average CI lengths for the same scenario. The averages of the $\beta_{t}$ and $\beta_{tx_{2}}$ estimates before deletion of data are 0.299 and 0.520, respectively. The MC-SD and 95% CI coverage rate for the $\beta_{t}$ estimates before deletion are 0.204 and 0.951, respectively. The MC-SD and 95% CI coverage rate for the $\beta_{tx_{2}}$ estimates before deletion are 0.327 and 0.953, respectively.

We first consider the MI methods. MI-R and MI-NR produce almost equivalent distributions of point estimates for both coefficients. We also see little difference in the MC-SDs for the two approaches. We ascribe these findings to the relatively large sample size. With 500 units in each treatment arm, we are able to estimate the parameters of the MI model fairly accurately, regardless of whether we combine the data or impute separately in each treatment arm. With binary covariate data, small differences in the precision of the parameters do not substantially alter the posterior predictive distributions for the missing values. The accuracy gain is further diluted, since only a modest fraction of cases are missing. Finally, we note that the analysis model adjusts for both $X_{1}$ and $X_{2}$. This can dampen the effect of the accuracy gain even further, as adjustment for covariates is known to increase precision against any residual imbalances that exist in spite of randomization. We see a similarly negligible effect of leveraging randomization when comparing MI-RY and MI-NRY.

For $\beta_{t}$, the outcome stage MI methods produce approximately unbiased estimates (simulated absolute bias $\leq$ 0.001). The simulated absolute biases for their design stage counterparts are higher ($\approx 0.04$). Differences are more pronounced for $\beta_{tx_{2}}$, with biases for methods MI-R and MI-NR escalating to 0.19 in the negative direction. As $\beta_{tx_{2}}$ measures how the conditional association between $Y$ and $X_{2}$ is modified by $T$, this coefficient is attenuated as a result of imputing $X_{2}$ without using $Y$, explaining the large negative bias.

The design stage MI methods exhibit efficiency gains over the outcome stage MI methods, i.e., they have smaller MC-SDs, especially for $\beta_{tx_{2}}$. The models that do not control for $Y$ eliminate the uncertainty from estimating coefficients relating the covariates and outcomes in the imputation model. However, this efficiency gain comes at the cost of bias that potentially could be substantial. Interestingly, in these simulations, one could argue that the efficiency gains outweigh the added bias; for example, the simulated mean-squared error is around 0.10 for MI-R and 0.19 for MI-RY.

For MI-R and MI-NR, the estimated standard errors for both coefficients are large compared to the corresponding MC-SDs. This accords with findings that the MI variance estimator can be positively biased (Wang and Robins, 1998; Robins and Wang, 2000; Reiter and Raghunathan, 2007), especially when the imputation and analysis models are not congenial (Xie and Meng, 2017). Such large standard errors often lead to wide confidence intervals containing the true parameter values, despite the biases in the point estimates of the coefficients; in fact, both MI-R and MI-NR provide coverage rates over 97% because of the over-estimation of variances. In contrast, MI-RY and MI-NRY appear not to suffer from over-estimation of variances, and as a result exhibit close to nominal coverage rates.

We next evaluate the performance of the mean imputation techniques. Mean-R and Mean-NR perform quite similarly in this scenario, generating nearly identical simulated biases and close to nominal confidence interval coverage rates for both coefficients. The MC-SDs and estimated SEs under Mean-NR are slightly larger than those under Mean-R, suggesting only minor benefits of leveraging the randomization when using mean imputation for this scenario.

Mean-NRY is a clear exception among the mean imputation methods. It generates significantly biased point estimates with the largest MC-SDs, and also has underestimated standard errors. As a result, it has abysmal confidence interval coverage rates.

Turning to the regression imputation approaches, we see that Reg-R and Reg-NR yield similarly biased estimates. The differences in efficiency are modest, with MC-SDs under Reg-NR about 2% larger than those under Reg-R. Thus, once again, leveraging randomization seems not to improve inferences meaningfully in this scenario. The biases in the coefficients under Reg-R and Reg-NR are akin to those for MI-R and MI-NR; however, the confidence interval coverage rates are noticeably lower. This reflects the impact of the bias, which is not compensated for by over-estimation of the variances.

The outcome version of regression imputation, Reg-NRY, offers mixed results. On the one hand, unlike Reg-R and Reg-NR, it results in coefficient estimates with only small simulated biases. However, the estimated standard errors for $\beta_{t}$ and $\beta_{tx_{2}}$ are too small compared to the MC-SDs, resulting in low coverage rates. This is an example of single imputation resulting in under-estimation of variance, emphasized by Rubin (1987) as motivation for MI.

Finally, we consider the performance of CCA. In this scenario, CCA avoids serious biases in $\beta_{t}$ and $\beta_{tx_{2}}$, but it results in considerably higher MC-SDs than those under competitive imputation methods. Unlike other methods, CCA does not take advantage of the partially observed information. The relatively large variance generates the widest confidence intervals across all methods, although the intervals have close to nominal coverage rates.

Looking across imputation methods, in this scenario, it is evident that Mean-NRY and Reg-NRY are not competitive procedures. Additionally, if one demands confidence interval coverage rates to be at least near-to-nominal, Reg-R and Reg-NR are dominated by the other imputation procedures. Finally, if one requires low bias in coefficient estimators, Mean-R (or Mean-NR) and MI-RY (or MI-NRY) are most effective. Arguably, neither method obviously outperforms the other in these simulations.

4.3.2 Scenario 1: Low association

Figure 2 displays the coefficient estimates under scenario 1 when the covariate-outcome association is low. Table 3 provides the corresponding standard errors, coverage rates, and mean CI lengths. The averages of the $\beta_{t}$ and $\beta_{tx_{2}}$ estimates before deletion of data are 0.294 and 0.017 respectively. The MC-SD and 95% CI coverage rate for the $\beta_{t}$ estimates before deletion are 0.191 and 0.955, respectively. The MC-SD and 95% CI coverage rate for the $\beta_{tx_{2}}$ estimates before deletion are 0.256 and 0.946, respectively.

When comparing the design stage methods that leverage random assignment to the corresponding methods that do not, we see at most minor differences in this simulation scenario. MI-R and MI-NR have nearly identical values across all performance metrics. MI-RY is around 1% to 2% more efficient than MI-NRY; Mean-R is around 2% to 4% more efficient than Mean-NR; and, Reg-R is around 3% more efficient than Reg-NR. As in the high association scenario, the sample size of 500 in each treatment arm evidently is large enough that the gains in precision in the imputation model parameters do not translate to substantial gains at the treatment effect estimation stage.

Looking across methods, we continue to see some of the patterns apparent in the high association setting. First, for MI-R and MI-NR, the MI variance estimator continues to be positively biased due to the uncongeniality of the imputation and analysis models. This results in confidence interval coverage rates around 99%, even though MI-R and MI-NR have the smallest average confidence interval lengths. Second, MI-RY and MI-NRY have substantially larger MC-SDs than MI-R and MI-NR, although they tend not to have the large positive bias in the MI variance estimator. We note that the design stage MI methods produce approximately unbiased coefficient estimates. Owing to the low association, the outcome does not carry much information about the partly missing covariates, which reduces the bias due to its exclusion in the covariate imputation model. Third, Mean-R and Mean-NR continue to perform similarly to MI-RY and MI-NRY. In this scenario, there is a hint of under-estimation in the variances associated with Mean-R and Mean-NR, which translates to slightly lower than nominal CI coverage rates. Fourth, Reg-R and Reg-NR continue to offer the smallest average CI lengths. In this scenario, however, the coverage rates for $\beta_{tx_{2}}$ are closer to nominal than in Table 2. As with the MI methods, ignoring the outcome in the regression imputations does not reduce the quality of the imputations meaningfully when the associations between the covariates and outcome are weak. Finally, Mean-NRY and Reg-NRY continue to be ineffective as imputation procedures.

4.3.3 Scenarios 2 and 3

Under scenarios 2 and 3, conclusions about the missing data methods do not fundamentally change. At $n=1000$, MI methods respecting and not respecting randomization produce similar point estimates, standard errors, and coverage rates. When the association between the covariates and the outcome is high, the design stage MI methods continue to be biased—with biases that are more marked than before—with lower standard errors, while the outcome stage MI methods have low bias with higher standard errors. When this association is low, all MI methods produce comparable point estimates. The single imputation methods and CCA continue to perform as in scenario 1. We present the related graphical and tabular displays in the supplementary materials.

5 Application

In this section, we present an application of the missing data methods to a randomized experiment analyzed in Foos and Gilardi (2019), available at https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/BSIFTF. The data comprise $n=612$ women, randomized to receive ($T=1$) or not to receive ($T=0$) an invitation to a career workshop in politics. Approximately two-thirds of the participants receive the treatment.

The primary behavioral outcome $Y$ is binary, with $Y=1$ if the participant applies to a political office mentoring program, and $Y=0$ if she does not. Several pre-treatment covariates are measured in the original experiment. We consider two that are highly associated with the outcome, namely an indicator of whether the participant takes an active interest in planning her career ($X_{1}$), and an indicator of whether the participant wishes to have children in the future ($X_{2}$). $X_{1}$ is binary, and exhibits low levels of missingness (2.3%). $X_{2}$ contains a “perhaps/don’t know” category, which we regard as missing data, as in, for example, Rubin et al. (1995) and Sadinle and Reiter (2017). This results in 26% missingness in $X_{2}$.

The distribution of the outcome variable is highly imbalanced, with only 1.5% of the participants applying to the mentoring program. This leads to analysis models with interaction terms exhibiting perfect prediction issues. We hence focus on estimating the regression-adjusted average treatment effect measured on the log-odds scale, and use the model logit ($\pi$) = $\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\beta_{t}T$, where $\pi=Pr(Y=1|X_{1},X_{2},T)$. The coefficient of interest is $\beta_{t}$. We note that there is non-compliance in the experiment, since not everyone in the treatment group attended the career workshop. Accordingly, we carry out an ITT analysis. We estimate $\beta_{t}$ using the same methods as in the simulations in Section 4. For MI, we use $M=100$ imputations for methods MI-R, MI-NR, and MI-NRY, and iterate to convergence for MI-RY.

Table 4 displays the resultant point estimates, standard errors and 95% confidence intervals for $\beta_{t}$. In line with our simulation results, we see that estimates and standard errors are nearly identical for MI-R and MI-NR and also for MI-RY and MI-NRY. The differences between design and outcome stage MI are not remarkable, likely due to the sparsity of $Y=1$ cases in the sample and low missingness levels in the covariates. Results also seem fairly insensitive to the choice of the identifying assumption in MI. Akin to the simulations, mean imputation as well as the regression imputation approaches perform with efficiency comparable to MI-NRY. All in all, for these data, the choice of imputation method seems not to impact results meaningfully.

We note that CC analysis produces higher standard errors and different point estimates in comparison to the imputation methods. This difference results mainly from the additional variability in CCA that goes along with deleting the partially observed cases. Because the outcome has low frequency–which leads to relatively large standard errors–removing the observations with missing covariates has a noticeable effect on the estimates.

6 Discussion

Our simulations show that, when units have been properly randomized in the design stage and sample sizes are moderately large, results are practically the same whether one respects or ignores randomization in common methods for imputing missing covariates. This is in slight contrast to the findings by White and Thompson (2005), who found that leveraging the randomization can improve inferences when using mean imputation. One key difference in the simulations underpinning those and our findings is the functional form of the outcome model. We estimate non-constant treatment effects for a binary outcome, whereas White and Thompson (2005) estimate homogeneous treatment effects for a continuous outcome. For these functional forms, it is evidently difficult to realize substantial efficiency gains for the sample sizes we consider. For randomized experiments where at least one treatment arm has a small size, say in the tens, the gains in precision when accounting for randomization can be more substantial than what we see here.

The results produced by the design and outcome stage imputation methods can differ notably, especially when the covariates are highly prognostic of the outcome and when one seeks to estimate heterogeneous treatment effects. For the MI methods, we observe a trade-off between bias and efficiency. Estimates produced by MI-R and MI-NR can have large bias with relatively low variances, whereas the opposite is true for the outcome stage MI methods. When the covariate-outcome association is low, all four MI methods yield estimates with low simulated bias, although MI-R and MI-NR continue to be more efficient.

Interestingly, in these scenarios, MI-R and MI-NR offer higher coverage rates and lower average interval length than MI-RY and MI-NRY. Of course, this pattern need not hold in other simulation scenarios. As an extreme illustration of this point, suppose we have the simulation set-up from Section 4.3.1 with very large $n$, say in the 100,000s in each treatment arm. In this case, the biases in MI-R and MI-NR for $\beta_{tx_{2}}$, which will persist at the levels in Table 2, will swamp the MI standard error, even with the variance inflation due to uncongeniality. Confidence interval coverage rates for MI-R and MI-NR will be well below nominal rates. In contrast, MI-RY and MI-NRY will have biases and coverage rates akin to those seen in Table 2. All in all, this suggests that, when one is interested in heterogeneous treatment effects, MI using the outcome potentially has some advantages over MI in the design stage. Of course, this assumes that analysts use a reasonable imputation model, which is hard to be certain of in practice.

Turning to the single imputation methods, it is evident that mean imputation in the design stage performs well in these simulations, offering close-to-nominal coverage rates with modest biases and reasonable average confidence interval lengths. These results accord with findings in Sullivan et al. (2018), who showed that mean imputation generates small biases in simulations with a single binary covariate and constant treatment effect. We note that, in large samples, even small biases can overwhelm standard errors, leading to below nominal coverage rates. Finally, one clear finding from the simulations is that the single imputation methods that account for the outcome perform worse than those that do not, suggesting that the former not be used for imputing missing covariates.

Our simulation scenarios are not exhaustive, and results may vary in alternative situations. By design, our findings are specific to analysis models that include an interaction between the treatment and a covariate. They may not apply when constant treatment effects are being estimated. For example, White and Thompson (2005) show that for a single covariate and no interaction effects, mean imputation in the design stage does not exhibit biases. When covariates are much more strongly associated with one another than in our simulation design, properly specified MI models that account for this association could be notably more efficient than mean imputation.

In practice, it is common to adjust for continuous, or a mix of continuous and categorical covariates in the causal analysis. The MI method in Section 2.2 can be adapted to such situations. For example, when $X$ is continuous, we can define the distribution $f(X_{obs},D)$ non-parametrically, using kernel density estimators, as in Titterington and Mill (1983). Alternatively, for each $D$, $f(X_{obs},D)$ can be given a parametric form (see, for e.g., Little (1993), who used normal densities). Sadinle and Reiter (2017) discuss these methods in detail by way of illustrations.

It is also possible that the outcome as well as the covariates contain missing values. In the framework of MI under non-parametric identification, we can separate the covariates and the outcome into blocks and place different identifying restrictions within these blocks. One such block-based method has been presented in Sadinle and Reiter (2018). We also note that non-parametric identification and subsequent imputation entails breaking down a data set by the observed missingness patterns. For this procedure to work well, sufficient numbers of data points per pattern are required. In practice, this amounts to having a large sample size, especially when the data have more than two covariates with missing values. When sample sizes are not adequate, analysts generally have to sacrifice non-parametric identification.

**SUPPLEMENTARY MATERIALS

** We provide two supplements as supporting materials for this article.