Bootstrap inference for multiple imputation under uncongeniality and misspecification

Multiple imputation has become one of the most popular approaches for handling missing data in statistical analyses. Part of this success is due to Rubin’s simple combination rules. These give frequentist valid inferences when the imputation and analysis procedures are so-called congenial and the embedding model is correctly specified, but otherwise may not. Roughly speaking, congeniality corresponds to whether the imputation and analysis models make different assumptions about the data. In practice, imputation models and analysis procedures are often not congenial, such that tests may not have the correct size, and confidence interval coverage deviates from the advertised level. We examine a number of recent proposals which combine bootstrapping with multiple imputation and determine which are valid under uncongeniality and model misspecification. Imputation followed by bootstrapping generally does not result in valid variance estimates under uncongeniality or misspecification, whereas certain bootstrap followed by imputation methods do. We recommend a particular computationally efficient variant of bootstrapping followed by imputation.


Introduction
Multiple imputation (MI) has proven to be an extremely versatile and popular tool for handling missing data in statistical analyses. For a recent review, see [13]. Its popularity is due to a number of factors. The imputation and analysis stages are distinct, meaning it is possible for one entity to perform the imputation and another the analysis. It is flexible, in being able to accommodate various constraints and restrictions that the imputer or analyst may want to impose. Auxiliary variables can be used in the imputation process to reduce uncertainty about missing values or make the missing at random (MAR) assumption more plausible, yet need not be included in the analyst's model.
In MI the analysis model of interest is fitted to each imputed dataset. Estimates and standard errors from each of these fits are pooled using 'Rubin's rules' [18]. These give a point estimate as the simple average of the imputed data estimates. Rubin's variance estimator combines the average within imputation variance with the between imputation variance in estimates. This requires an estimator of the complete data variance, which for most estimators is available analytically. In Rubin's original exposition the estimand was characteristic of a fixed finite population of which some units are randomly sampled and data are obtained [18]. Rubin defined conditions for an imputation procedure to be so called 'proper' for a given complete data analysis. If in addition the complete data analysis gives frequentist valid inferences, MI using Rubin's rules yields valid frequentist inferences [18,17,13]. Subsequently Rubin's rules were criticised by some (e.g. [7]) because in certain situations Rubin's variance estimator could be biased relative to the repeated sampling variance of the MI estimator. In response, Meng defined the concept of congeniality between an imputation procedure and an analyst's complete (and incomplete) analysis procedure [11]. If an imputation and analysis procedure are congenial, this implies the imputation is proper for the analysis procedure [14]. Meng showed that for certain types of uncongeniality, Rubin's variance estimator is conservative, ensuring the intervals have at least the advertised coverage level [11]. In other settings however it can be biased downwards, leading to undercoverage of confidence intervals [16].
Rubin's rules have proved fantastically useful since MI's inception, in particular because they facilitate the separation of imputation and analysis into two distinct parts and because they are so simple. Nevertheless, in settings where Rubin's variance estimator is asymptotically biased, if feasible, the analyst may desire sharp frequentist valid inferences. Robins and Wang proposed a variance estimator which is valid without requiring congeniality or correct model specification [16]. Their estimator requires calculation of various quantities depending on the estimating equations corresponding to the particular choice of imputation and analysis models. As such it is arguably harder to apply their approach when the imputer and analyst are separate entities. As far as we are aware, its use has been extremely limited thus far in practice due to these requirements.
Combining bootstrapping with MI was first suggested over 20 years ago [22] and recently a number of papers have investigated a wider variety of approaches to combining them. Schomaker and Heumann investigated four variants which combined bootstrapping with multiple imputation [19]. Their motivation for exploration of using bootstrap with MI was for situations where an analytical complete data variance estimator is not available, or one is concerned that the MI estimator is not normally distributed. On the basis of theoretical and empirical investigation, they recommended three of the four variants for use. They did not explicitly seek to investigate performance under uncongeniality or model misspecification however. von Hippel and Bartlett proposed an alternative combination of bootstrapping with MI in the context of proposing frequentist type (improper) multiple imputation algorithms, and noted that it would be expected to valid under uncongeniality [25]. Lastly, Brand et al investigated six different combinations of MI with bootstrapping in the context of handling skewed data, and recommended using percentile bootstrap confidence intervals with single (stochastic) imputation [3].
In this paper we investigate the properties of the different combinations of MI and bootstrap which have been recommended by these previous papers, giving particular emphasis to their validity under uncongeniality. In Section 2 we review MI and Rubin's combination rules. In Section 3 we describe the various combinations of bootstrapping and MI that have been recently recommended and consider their validity under uncongeniality. Section 4 presents two sets of simulation studies, empirically demonstrating the impacts of uncongeniality on the frequentist performance of the different variants. We conclude in Section 5 with a discussion.

Multiple imputation using Rubin's rules
In this section we briefly review MI and Rubin's combination rules. In MI we first create M imputed datasets. We fit our complete data model to each, obtaining estimatesθ m , m = 1, .., M , and corresponding within imputation variance estimates Var(θ m ). The estimate of θ is then given byθ M = M −1 M m=1θ m , while the variance is estimated by Var(θ m ) For subsequent developments, following von Hippel and Bartlett [25] it will be useful to express each imputation estimate of θ as:θ ∞ , E(a m ) = 0 and Var(a m ) = σ 2 btw . Since the imputation estimates are conditionally independent givenθ ∞ , we have that Assuming the complete data analysis would provide valid frequentist inferences with complete data, if the imputation procedure is proper with respect to the complete data procedure [18], we have that Inference forθ M is then performed assuming that is t-distributed with degrees of freedom given by These results were derived assuming that with complete data the degrees of freedom are infinite and M is finite. In small sample settings the former assumption is questionable, and so Barnard & Rubin subsequently proposed a small sample version of Rubin's rules [1].
Meng subsequently defined an imputation procedure and a complete data analysis to be congenial essentially if there exists a Bayesian joint model for which the posterior distribution of the missing data matches that used by the imputation procedure and for which the complete data posterior mean and variance of the parameters of substantive interest are asymptotically the same as obtained by using the complete data analysis procedure [11]. Meng's congeniality definition in fact incorporated an additional notation of the analyst's incomplete data procedure, but for the present purposes this aspect is not relevant.
When the imputation and complete data analysis procedures models are congenial, this implies the imputation procedure is proper for the complete data analysis, and if in addition the complete data analysis gives frequentist valid inferences, Rubin's variance estimator for finite M is asymptotically unbiased [14]. When this is not case, Rubin's variance estimator can, depending on the configuration, be downwardly or upwardly biased as an estimator of the repeated sampling variance ofθ M [11,16].
Robins and Wang proposed a variance estimator for MI when each dataset is imputed using the maximum likelihood estimate of a parametric imputation model and the imputations are analysed using a non, semi or fully parametric model [16]. Their variance estimator is consistent without requiring the imputation and analysis models to be congenial nor even correctly specified. Hughes et al compared Robins and Wang's proposal to Rubin's rules through a series of simulation studies where the imputation and analysis models were misspecified and/or uncongenial with each other [8]. They demonstrated that Rubin's rules inference could be conservative or anti-conservative, whereas, at least for moderate or large sample sizes, inferences based on Robins and Wang's proposal were valid across their simulation scenarios. Hughes et al noted however that a major practical obstacle to the widespread use of Robins and Wang's method is its implementation is specific to the particular imputation and analysis models, and no software currently implements it.

Combining bootstrapping and multiple imputation
In this section we review the combinations of bootstrapping and MI which have been recommended for use in the recent literature, and consider whether their validity in uncongenial settings.

Imputation followed by bootstrapping
The first collection of methods we consider are where MI is first applied, and then bootstrapping is applied to each imputed dataset.

MI boot Rubin
The first combination considered (and recommended) by Schomaker and Heumann [19] is standard MI using Rubin's rules, but using bootstrapping to estimate the within-imputation complete data variance: This approach is what has often been used when no analytical estimator for the full data variance σ 2 wtn is available, or if one is concerned about whether the analysis model is correctly specified. In the latter case, a sandwich variance estimator has sometimes been used to attempt to provide robustness to misspecification [8].
Provided the analysis model would give valid frequentist frequentist inferences with complete data, we expect MI boot Rubin to give asymptotically unbiased variance estimates ofθ M under congeniality. This is supported by the setting 1 simulation results of Schomaker and Heumann [19]. However, since this approach relies on Rubin's rules, we would not expect it to give unbiased variance estimates in general under uncongeniality. This hypothesis is supported by Schomaker and Heumann's setting 2 with high missingness simulation results, where we believe the imputation and analysis models are uncongenial, and where coverage for one parameter was 91%.

MI boot pooled percentile
The second approach considered and recommended by Schomaker and Heumann [19] is the same as MI boot This approach can be viewed as a route to obtaining a posterior credible interval, and hence assuming the analysis model would give valid inferences with complete data and that the imputation and analysis models are congenial, we expect it to give asymptotically unbiased variance estimates ofθ M . This is because first draws are taken from the posterior of the missing data given observed, and second, conditional on these, bootstrapping and estimating the parameters by their maximum likelihood estimate is in large samples equivalent to taking a draw from the posterior given the imputed missing data and the observed data [9].
To explore this further, we can express the estimate from the mth imputation and bth bootstrap aŝ where a m is as defined earlier and b mb is a term representing the deviation due to bootstrap. The term b mb has mean zero and variance σ 2 wtn , since the between bootstrap variance for a completed dataset corresponds to the complete data variance. For M B large, the sample variance of the pooled sample of M B estimates, which we are effectively treating as a size M B sample from the posterior, is: From standard results for the one-way random intercepts model [21], this is an unbiased estimator of Schomaker and Heumann [19] considered large values of B (e.g. 200) and smaller values of M . For large B the preceding expression is approximately equal to Thus if both M and B are large, this is unbiased for σ 2 btw +σ 2 wtn , which under congeniality is the true posterior variance. If M is not large however, it is biased downwards for the true posterior variance, and so we would expect confidence intervals constructed using the M B sample of estimates, e.g. based on percentiles as suggested by Schomaker and Heumann, to undercover. This concurs with the findings shown in Figure 1 of Schomaker and Heumann, who found that the percentile MI boot pooled confidence intervals undercovered somewhat for small M even under congeniality (Figure 1 of [19]).
Under uncongeniality, there is no reason to expect this approach to result in valid inferences, and the setting 2 with high missingness simulation results of Schomaker and Heumann support this, with coverages between 89% and 92%.

Bootstrap followed by MI
We now consider methods which first bootstrap sample the observed data and then apply MI to each bootstrap sample. This general approach to combining bootstrap with MI was proposed by Shao and Sitter [22] and Little and Rubin [9].

Boot MI percentile
Both Schomaker and Heumann [19] and Brand et al [3] recommended calculating bootstrap percentile in- This approach to direct application of the standard percentile based bootstrap confidence interval to the estimatorθ M and as such, as suggested by Shao and Sitter, we expect it to be asymptotically valid even under uncongeniality [22]. Moreover, provided the point estimator is consistent, asymptotically the resulting confidence intervals should attain norminal coverage even if the imputation and/or analysis models are misspecified. In Schomaker and Heumann's setting 2 simulation results, where we believe the imputation and complete data models are uncongenial, they found coverage rates close to 95%, although for one parameter it was as low as 90%.
Brand et al also found that the Boot MI % approach worked well in simulations [3]. They investigated it using either M = 5 or M = 1, and among the different combinations of bootstrapping and MI recommended using it M = 1. Although we expect the resulting confidence intervals to be valid, even under uncongeniality, we expect the intervals to be unnecessarily wide with M = 1 because only one imputation is used per bootstrap. This is confirmed by the simulation results of Brand et al [3] (Figure 1, panel C), which shows that the bootstrap percentile intervals were wider on average with M = 1 compared with M = 5. Moreover, their results suggested that coverage with M = 1 was slightly above the nominal 95% level, which we investigate further in Section 4.2.

Boot MI von Hippel
Of the various combinations of bootstrapping and imputation described, assuming the MI point estimator is consistent, only Boot MI percentile is expected to give confidence intervals that attain nominal coverage (asymptotically) under uncongeniality or model misspecification. A practical issue however is that the computational burden is high. For standard applications of MI, it is not uncommon now for M to be chosen as 100 or greater, for reasons of statistical efficiency of point estimates and to reduce Monte-Carlo error to an acceptable amount [26,10,24]. For bootstrap confidence intervals, the number of replications B is generally recommended to be at least 200 for variance estimation and at least 1,000 for percentile based intervals [6]. These considerations would imply a potentially very large value of BM , which may be computationally expensive or impractical.
von Hippel and Bartlett proposed an alternative point estimator and confidence interval based on Boot MI [25]. He proposed usingθ BM = B −1 B b=1θ b rather thanθ M , as the point estimator. To construct a confidence interval, von Hippel and Bartlett noted that the estimatesθ b,m can be expressed as: Given this variance components model, we have that This shows that provided B is large,θ BM will have similar efficiency toθ M with M large. The two variance components σ 2 ∞ and σ 2 btw can be estimated by fitting a one-way analysis of variance to the point estimateŝ

Simulations
In this section we report two simulation studies to empirically demonstrate the performance of the previously described combinations of bootstrapping and MI under uncongeniality and/or model misspecification.

Regression models under uncongeniality and misspecification
We first compared the previously described bootstrap and MI combination methods in four scenarios of uncongeniality and/or misspecification of the imputation and analysis models using a simulation study based on one performed by Hughes et al [8].
Briefly, we simulated a hypothetical dataset of one binary variable, sex, and four continuous variables, age, height, weight and natural log of insulin index (hereafter referred to as loginsindex). The data were generated under the following model: sex ∼ Bernoulli(π), age, height|sex ∼ N (α 0 + α 1 sex, Σ), where error W and error L are independent errors and η sex = 1 when sex= 0 and η sex = η when sex= 1.
Parameter values are shown in Table 4.1. Different scenarios were created by setting parameters α 1 , ι 1 , ν and β 1 to their null values. The values of the remaining parameters were fixed. Weight measurements were set to be missing completely at random for 60% of the observations.
The analysis of interest was to estimate θ, the effect of weight on loginsindex after adjustment for age and sex. Both imputation and analysis models were normal linear regression models that assumed homoscedastic errors. Unless otherwise stated, the distributions of error W and error L were normal, weight measurements were missing in men and women, the assumption of homoscedastic errors was true, and the imputation and analysis models were fitted to the entire sample. The following four scenarios were considered: Subgroup analysis scenario The data were simulated such that the continuous variables were identically distributed in men and women, and weights was missing among men only. The imputation and analysis models were uncongenial since the analysis model was fitted to men only whilst the imputation model was fitted to the entire sample ignoring sex (i.e. excluding sex as a predictor).
Heteroscedastic errors The data were simulated such that the variance of weight and loginsindex differed between men and women. The imputation and analysis models were congenial but incorrectly specified because they assumed homoscedastic errors.
Omitted interaction As in all scenarios, the data were simulated such that the effect of weight on loginsindex was the same for men and women. The imputation and analysis models were uncongenial because the analysis model included an interaction term between weight and sex whilst this interaction was, correctly, omitted from the imputation model.
Moderate non-normality Error distributions error W and error L were simulated from the log-normal distribution exp{N (0, 1/4 2 )}. The imputation and analysis models were congenial, but misspecified because they assumed a normal error distribution.
For each scenario we generated 1, 000 independent simulated datasets, where the sample size was 1, 000 observations and the probability of observing weight was 0.4, except for the subgroup analysis scenario where the probability of observing weight was 1 among women and 0.4 among men. We conducted MI Rubin using implying that the estimated coverage probability should lie within the range 0.936 and 0.964 (with 95% probability) [12].
For all methods, the point estimates of θ were either unbiased or the amount of systematic bias was trivial (e.g., at most −0.000289; results available on request). Tables 2 and 3 show the median of the CI widths and CI coverage for the 6 methods under comparison. For the subgroup analysis scenario (Table 2), MI Rubin and both MI then bootstrapping methods resulted in confidence interval overcoverage. Narrower confidence intervals and nominal coverage were achieved with the boot MI percentile method with 10 imputations and boot MI von Hippel. Boot MI percentile with single imputation resulted in wide confidence intervals and overcoverage. This concurs with what was found in the simulations reported by Brand et al. In the Appendix we give a sketch argument for why the Boot MI percentile intervals with M = 1 (or indeed small M more generally) will overcover. Interestingly, this over-coverage does not similarly affect normal based (as opposed to percentile) Boot MI intervals with M = 1.
For the heteroscedastic errors scenario (Table 2), MI Rubin and both MI then bootstrapping methods resulted in confidence interval undercoverage. Again, the boot MI percentile method with 10 imputations and boot MI von Hippel were the best performing methods with close to nominal coverage. The results for the omitted interaction scenario (Table 3) followed a similar pattern noted for the subgroup analysis scenario. For the moderate normality scenario (Table 3), MI boot pooled percentile had slight confidence interval undercoverage and boot MI percentile with single imputation overcovered. The remaining methods had close to nominal coverage with similar median CI widths.

Reference based imputation in clinical trials
Our second simulation study setting is so called control or reference based MI for missing data in randomised trials. Missing data due to study dropout is common on clinical trials, and there is often concern that missing data do not satisfy the missing at random (MAR) assumption. Often dropout in trials coincides with patients' treatments changing. An increasingly popular approach to imputing missing data in trials is using so called reference or control based MI approaches [4]. These involve constructing the imputation distribution for the active treatment arm using a combination of information from the active and control arms, which results in uncongeniality between imputation and analysis models. This uncongeniality results  consistently estimates a sensible variance in the context of MAR sensitivity analyses [5]. We do not here enter this debate, but merely investigate the previously described bootstrap and MI combinations in regards their ability to produce confidence intervals with the correct repeated sampling coverage. In the setting of reference based MI Quan et al applied (we believe) Boot MI to estimate standard errors ofθ M , and found it worked well [15]. We simulated 10,000 datasets of size n = 500 with 250 randomised to control (Z = 0) and 250 (Z = 1) randomised to active treatment. Baseline X and outcome Y were then generated from a bivariate normal model: The analysis model was normal linear regression of Y on X and Z, with the coefficient of treatment Z of primary interest. Values in Y were made missing complete at random with probability 0.5. For each dataset, first the missing values in Y were imputed using a normal linear regression model with X and Z as covariates assuming MAR, such that the imputation and analysis models were congenial. Second they were imputed using the jump to reference method (see [4] for details), such that the two models were uncongenial. The same combinations of bootstrapping and MI were used as in the first simulation study.

Discussion
We have reviewed a number of proposals for combining MI with bootstrapping, in particular in regards their statistical validity when imputation and analysis models are uncongenial or misspecified. Approaches which first impute then bootstrap generally do not give valid inferences under uncongeniality or model misspecification. In contrast, bootstrapping followed by imputation is robust to uncongeniality, and provided the MI point estimator is consistent, the resulting confidence intervals have the correct coverage asymptotically even if the imputation and/or analysis models are misspecified. A drawback of this approach with M large is its high computational cost. Brand et al recommended this method, but using M = 1 imputation per bootstrap, which obviously reduces the computational burden considerably [3]. However, with small M , the MI point estimator is inefficient, and moreover we have shown that percentile based Boot MI intervals over-cover even under congeniality for small M . The alternative boot MI version proposed by von Hippel and Bartlett overcomes these drawbacks, only requiring BM imputations and analyses, and M can be chosen to be two. It does however, like Rubin's rules, assume that the MI estimator is normally distributed. The Boot MI von Hippel approach is implemented in the R package bootImpute, and is available from CRAN [2]. As far as we are aware the only alternative approaches for valid inferences under uncongeniality require complex problem specific calculations which are not conducive to general use [16,23], and in this context the Boot MI von Hippel approach seems very attractive.
As mentioned in the introduction, Rubin originally envisaged the imputer and analyst as distinct entities, with the imputer releasing a single set of multiply imputed datasets to different analysts. A strength of the bootstrap followed by MI approach is that this division of roles is still feasible -the imputer bootstraps and then multiply imputes the observed data, releasing a set of imputations clustered by bootstrap. These can then be analysed by different analysts, and inferences can be obtained using either the boot MI percentile or Boot MI von Hippel approaches.
Combining bootstrapping with MI has some disadvantages compared to inference using Rubin's rules. Compared to regular MI with Rubin's rules, it is considerably more computationally intensive -this is the price paid for being able (in certain situations) to obtain valid inferences under uncongeniality or misspecification. Problems with model (imputation or analysis) convergence are probably more likely due to the large number of bootstraps required. The non-parametric resampling scheme used by bootstrapping relies on an assumption that the data are independent and identically distributed, and further research is warranted to explore the use of other types of bootstrap resampling schemes in conjunction with MI.
Code for the first simulation study (R) and the second simulation study (Stata) are available from https://github.com/jwb133/bootImputePaper.