Multiple Imputation in Two Stages 

Ofer Harel and Joseph L. Schafer 

Ofer Harel, University of Washington 

F-600 Health Sciences Building, Seattle, WA 98195, ofer@u.washington.edu 

Joseph L. Schafer, Pennsylvania State University

325 Thomas Building, University Park, PA 16802, jls@stat.psu.edu. 

Abstract

Conventional multiple imputation (MI) (Rubin, 1987) replaces the missing values in a dataset by m > 1 sets of simulated

values. We describe a two-stage extension of MI in which the missing values are partitioned into two groups and imputed N = 

mn times in a nested fashion. Two-stage MI divides the missing information into two components of variability, lending 

insight when the missing values are of two qualitatively different types. It also opens new possibilities for making different

assumptions about the mechanisms producing the two kinds of missing values. Point estimates and standard errors from the N

complete-data analyses are consolidated by simple rules derived by analogy to nested analysis of variance. After reviewing 

the theory and practice of two-stage MI, we illustrate the method with a simple analysis of binary variables from a 

longitudinal survey. 

Introduction 

Missing values in a dataset may be of two different types. In response to a sensitive question in an interview survey, for 

example, some participants may simply refuse to answer, whereas others may say, “I don’t remember.” Other examples

include unit versus item nonresponse, planned missingness (e.g. as would arise from double sampling or matrix sampling)

versus unplanned missing values, and mortality versus dropout for other reasons. In this article, we explore a two-stage

version of Rubin’s (1987) multiple imputation in which we impute the first kind of missing value m times; then, for each 

imputation of the first type, we impute the second type n times, treating the imputed values for the first type as if they were 

fixed and known. 

Imputing in stages has several potential advantages. From a computational standpoint, it is sometimes possible to identify a 

small set of missing values which, if known, would simplify the process of imputing the rest. An application of this type is

described by Rubin (2003), who first completed a monotone pattern by computationally intensive iterative procedure, then

imputed the remaining missing values by a less demanding, non-iterative method. Second, two-stage MI allows us to identify 

the amount of uncertainty in population estimates contributed by each type of missing value, which may have important 

implications for interpreting current results and planning future studies. 

Finally, two-stage MI may allow us to posit different assumptions about the probabilistic mechanisms generating the two

types of missing values. Suppose that some values are missing by design but others are missing for reasons beyond the data 

collectors’ control; the missing at random (MAR) condition (e.g., Little and Rubin, 1987) is known to hold for the former but 

not the latter. With two-stage MI, it may be feasible to impute for unplanned missing values under some alternative nonMAR methods, then impute for planned missing values under an MAR assumption. In that case, however, we would need to

clarify what it means for some missing values to be MAR and others to be non-MAR. Making this notion precise requires us 

to extend Rubin’s (1976) concepts of MAR and ignorability to mechanisms that partition the complete data into three parts 

(observed, missing for one reason, missing for the other reason). A full treatment of that topic has been given by Harel (2003)

but is beyond the scope of this paper. 

In many cases, two-stage MI can be carried out by repeatedly applying algorithms and software designed for conventional

MI, such as the missing-values library in S-PLUS or SAS PROC MI. Once the N = mn imputed datasets exist, they are 

analyzed by complete data methods, and the results are combined by simple algebraic rules derived by Shen (2000). In the 

remaining sections, we describe in generic terms the methods for generating the imputations and combining the results; we 

also discuss how to estimate the amount of missing information due to each set of missing values. We conclude with a 

simple example involving a binary variable from two adjacent waves of the National Crime Survey. 

Page 91 
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Imputing in Two Stages 

Let Ycom = (Yobs,Ymis ) denote the complete dataset, part of which is observed (Yobs ) and part of which is missing (Ymis ) . In 

conventional MI, we would typically propose a parametric model for the complete data, say P(Y |θ ), and a Bayesian prior com

(1) (2) (m)

distribution P(θ ) for the unknown model parameters. The m multiple imputations, Y ,Y ,K,Y , are independent draws mis mis mis

from the posterior predictive distribution

P(Y |Y ) = ∫ mis P(Y |Y ,θ ) P(θ |Y ) dθ , obs mis obs obs

which reflects uncertainty about θ as well as Ymis . Computational algorithms for producing these draws under a variety of 

multivariate models are described by Schafer (1997). Software for MI is available from a variety of commercial and noncommercial sources; see Horton and Lipsitz (2001) for a review. 

A B Now suppose that the missing data are divided into two parts, Ymis = (Ymis ,Ymis ) . To carry out two-stage MI, we first draw m

A(1) A(2) A(m) A

independent values Ymis ,Ymis ,K,Ymis the posterior predictive distribution P(Ymis |Yobs ) . Then, for each of these, we draw

B( j,1) B( j,2) B( j,n) B A( j)

n conditionally independent values Ymis ,Ymis ,K,Ymis from P(Ymis |Yobs,Ymis ), j =1,K,m . Care must be taken to

reflect uncertainty in the parameterθ at each stage. The resulting N = mn completed datasets, 

( j,k ) A( j) B( j,k )

Y = (Y ,Y ,Y ), j =1,K,m, k =1,Kn ,

com obs mis mis

B( j,1) B( j,2) B( j,n) A

are not independent, because each block or nest Ymis ,Ymis ,K,Ymis contains identical values for Ymis . Notice, 

(1,1) (2,1) (m,1) however, that gathering the first completed dataset from each nest, Y ,Y ,K,Y , does give us m independent com com com

imputations from P(Ymis |Y ) . Therefore, we can perform two-stage MI by first creating m conventional imputations for obs

B A

Ymis , then discarding Ymis and recreating n −1 additional imputations of it for each imputed version of Ymis . 

In this discussion, we have not made use of the response indicators, the set of binary random variables that indicate for each 

element of Y whether the item is observed or missing. In conventional MI, we would need to condition upon the response com

indicators if the missing data were not MAR (Rubin, 1987). With two-stage MI, we might also need to consider the process

that bifurcates Ymis ; that is, we might need to condition upon the set of three-level indicators that partition Ycom into Yobs , 

A B Ymis and Ymis . Using an extended theory of ignorability presented by Harel (2003), one can show that these three-level 

A B indicators can ignored in two-stage MI if (a) Ymis is MAR, and (b) the process that dividesYmis into Ymis and Ymis does not

depend on any portion of Ymis . Condition (b) allows a data analyst classifies the missing values into groups on the basis of

A

information gleaned from Yobs or from the realized pattern of nonresponse—for example, by takingYmis to be a set of missing 

values that would complete a monotone pattern, as suggested by Rubin (2003). Other conditions under which part or all of

the missingness mechanism can be ignored in one or both stages of two-stage MI are described by Harel (2003). 

Consolidating Results from Post-Imputation Analyses

After imputation, each of the N = mn imputed datasets from two-stage MI should analyzed as if it had no missing values.

Estimates and standard errors from these analyses may then be consolidated using Shen’s (2000) extension of Rubin’s (1987)

ˆ( , ) ( , ) method for conventional MI. Let Q denote a scalar population quantity to be estimated, and let Q

j k and Uj k denote the 

point and variance estimates for Q calculated from the ( j,k) th completed dataset. The overall point estimate for Q is 

simply the grand average

1 m n

ˆ( j,k ) 1 m ( j, ⋅) Q⋅⋅ = ∑ ∑ Q = ∑ Q ,

mn j=1 k=1 m j=1

( j, ⋅) 1 n

ˆ( j,k ) where Q =

k 1

Q . The uncertainty associated with this overall estimate involves three components: the 

n ∑ =

estimated complete-data variance 

Page 92 
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

⋅⋅

⋅⋅

1 m n ( j,k ) U = U ,

⋅⋅ mn∑j=1∑k=1

the between-nest variance 

1 m ( j,⋅)

2

B = ∑ =

(Q − Q⋅⋅) ,

m j 1 −1

and the within-nest variance

1 m 1 n

ˆ( j,k ) ( j,⋅)

2

W = ∑ ∑ (Q −Q ) .

m n j=1 k=1 −1

The total variance is 

⎛

−

1 ⎞ ⎛ 1 ⎞

T =U⋅⋅ + ⎜1 ⎟W + ⎜1+ ⎟B ,

⎝ n ⎠ ⎝ m ⎠

and inferences about Q are based on the Student’s t-approximation (Q − Q ) T ~ tν with degrees of freedom 

2 2

1 ⎛ (1−1 n)W ⎞ 1 ⎛ (1+1 m)B⎞

ν = ⎜ ⎟ + ⎜ ⎟ . 

m(n −1)⎝ T ⎠ m −1⎝ T ⎠

Qˆ( j,k ) These rules arises from analysis-of-variance techniques, regarding the ’s as measurements from a balanced experiment

A with Ymis as a random blocking factor. For the special case of n =1, the term involving W drops out of the total variance, and 

the method becomes identical to that of Rubin (1987) for conventional MI with m imputations. 

Rates of Missing Information and Relative Efficiency

Shen (2000) did not consider rates of missing information, but these rates are easily estimated from the variance components.

The estimated overall rate of missing information for Q is 

ˆ

B + (1−1 n)W rˆ

λ = = ,

U⋅⋅ + B + (1−1 n)W 1+ rˆ

where rˆ = (B + (1−1 n)W) U is the relative increase in variance due to nonresponse. The estimated rate of missing

B A information due to Ymis if Ymis were known is

B|A

ˆB |A W rˆ

λ = = ,

U B|A

⋅⋅ + W 1+ rˆ

B|A ˆA ˆ ˆB|A where rˆ = W U⋅⋅

. The difference λ = λ − λ estimates the amount by which the missing information would drop if

A ˆA ˆ A

Ymis became known, and λ λ is the proportionate reduction in missing information if Ymis became known. It is possible for 

λ

ˆA

to be negative, in which case the estimate should be set to zero. 

Rubin (1987) showed that the relative efficiency of a point estimate based on m conventional MI’s to a fully efficient one (i.e. 

−

based on an infinite number of imputations) is approximately (1+ λ m)

1

, where λ is the overall rate of missing information 

−

for Q . For two-stage MI, one can show that the relative efficiency of Q⋅⋅lies between (1+ λ m)

1

and (1+ λ (mn))−1, 

A

reaching the upper bound when Ymis carries no information about Q . The degrees of freedom for estimating the within- and 

between-nest variance components are approximately equal when n = 2 . In practice, we typically set n = 2 and then choose

m large enough to achieve a reasonable level of efficiency. As with conventional MI, only a few imputations are usually

needed to achieve a highly efficient estimate. Unfortunately, the estimated rates of missing information can be rather noisy

when m is small. Asymptotic results on the variability of these estimated rates are given by Harel (2003). 

A Simple Example 

The data in Table 1, previously analyzed by Kadane (1985) and Schafer (1997), were obtained from the National Crime

Survey conducted by the U.S. Bureau of the Census. Occupants of sampled housing units were asked whether they had been 

Page 93 
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

victimized by crime in the preceding half-year; six months later, occupants of the same units were asked the question again. 

Of the 765 sampled units, 561 (74%) provided responses at both occasions, 42 (6%) responded only the first time, 38 (5%)

responded only the second time, and 115 (15%) did not respond at either time.

Table 1: Victimization status of housing unit occupants in two successive six-month periods

Victimized second period? 

Victimized first period? No Yes Missing 

No 392 55 33 

Yes 76 38 9 

 Missing 31 7 115 

For simplicity of presentation, we will regard this as a simple random sample from the population of interest. Let Yt denote 

victimization status (1=no, 2=yes) during period t =1,2 , and let π ij denote the proportion of units in the population with 

Y1 = i,Y2 = j . Assuming a multinomial model with ignorable nonresponse, the loglikelihood function is

l = 392 log π + 55 log π + 76 log π + 38 log π + 33 log π + 9 log π + 31log π + 7 log π +115 log π , 11 12 21 22 1+ 2+ +1 +2 ++

where ‘+’ denotes summation over a subscript. Maximum-likelihood (ML) estimates computed by the EM algorithm of

Chen and Fienberg (1974) are πˆ

11 = .6971, πˆ12 = .0986 , πˆ21 = .1358 and πˆ22 = .0685 (Schafer, 1997, pp. 42-45). 

Multiple imputations for missing data in cross-classified contingency tables under the multinomial model can be generated 

using a data augmentation (DA) algorithm, a Markov chain Monte Carlo procedure described by Schafer (1997, Chap. 7).

This algorithm, which has been implemented in the missing-data library of S-PLUS (Schimert et al., 2001), can also be used

A B

for imputation in two stages. To illustrate, we defined Ymis to be the missing values of Y1

and Ymis to be the missing values

of Y 2

. Using the default prior density 

−.5 −.5 −.5 −

P(π ,π ,π ) ∝ π π π (1−π −π −π )

.5

, 11 12 21 11 12 21 11 12 21

we generated ten imputations in two stages with m = 5 and n = 2 in the following manner. First, we generated five ordinary

MI's using 100 cycles of DA between imputations. Then, for each of the five datasets, we removed the imputed values for Y2

and re-imputed them once using another 100 cycles of DA. Frequencies for the ten imputed datasets are shown in Table 2.

Table 2: Frequencies for victimization status in two periods after two-stage multiple imputation, with estimated odds ratios 

and differences 

j = 1 j = 2 j = 3 j = 4 j = 5 

k = 1 k = 2 k = 1 k = 2 k = 1 k = 2 k = 1 k = 2 k = 1 k = 2 

No, No 536 531 526 523 528 526 526 530 532 527 

No, Yes 70 75 75 78 71 73 75 71 72 77 

Yes, No 102 92 104 101 108 102 101 107 103 102 

Yes, Yes 48 58 51 54 49 55 54 48 49 50 

αˆ 3.60 4.46 3.44 3.58 3.37 3.89 3.75 3.35 3.52 3.36 

logαˆ 1.28 1.50 1.24 1.28 1.22 1.36 1.32 1.21 1.26 1.21 

SE(logαˆ) 0.22 0.21 0.21 0.21 0.21 0.21 0.21 0.21 0.21 0.21 

δ

ˆ .042 .022 .038 .030 .049 .038 .048 .048 .041 .033 

SE(δ

ˆ) .017 .017 .018 .018 .018 .017 .018 .018 .017 .018 

Page 94 
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Notice that in each of the blocks j =1,K,5 , the imputed values for the marginal Y1frequencies are constant. Two parameters 

−1 −1

of interest are the odds ratio α = π π π π and the change in victimization rate δ = π −π = π −π . Estimates of 11 12 21 22 2+ +2 21 12

α , logα and δ from each imputed dataset are shown in Table 2, along with standard errors calculated as 

−1 −1 −1 −1

SE(logαˆ) = x + x + x + x 11 12 21 22

and 

x ⎛ x ⎞ x ⎛ x ⎞ x x ˆ 12 12 21 21 12 21 SE( ) δ = ⎜1− ⎟ ⎜1− ⎟ 2 ,

2 ⎜ ⎟

+

2 ⎜ ⎟

+

3

x x x ++ ⎝

x++ ⎠ ++ ⎝x++ ⎠ ++

where xij is the sample frequency of Y1 = i,Y2 = j (e.g. Agresti, 1990). Combining these quantities by Shen’s rules, the 

overall estimate of logα becomes Q⋅⋅ =1.29 with a standard error of T = 0.23 , and the estimate of δ becomes Q⋅⋅ = 0.38

with a standard error of T = 0.19 . Although these estimates are based only on a small number of imputations, the large 

values for the degrees of freedom (ν = 345 and ν = 282 , respectively) indicate that the loss of precision relative to estimates 

based on more imputations is slight. For comparison, we computed asymptotic standard errors by evaluating the Hessian of

the actual loglikelihood at the ML estimate; this loglikelihood method gives logαˆ =1.27 with a standard error of 0.26, and 

δ

ˆ = .037 with a standard error of .020. 

Although the estimates and standard errors from such a small number of imputations are reliable, the estimated rates of 

missing information are not. To estimate these rates very precisely, we repeated the procedure to create m = 500 blocks of

n = 2 imputations each. The estimated parameters, standard errors and rates of missing information from this larger

simulation are displayed in Table 3. 

Table 3: Estimated parameters, standard errors and rates of missing information for crime victimization data from multiple 

imputation with m = 500 and n = 2. 

est SE λ

ˆ λˆB|A

λ

ˆA λ

ˆA ˆλ

logα 1.27 0.25 .27 .21 .06 .24 

δ .038 .020 .21 .13 .08 .39 

The rates of missing information are intuitively reasonable. First, consider the log odds ratio, which measures the strength of

association between Y1

and Y2. Because 26% of the units had missing values for one or both of these variables, we would 

expect the overall rate of missing information to be approximately 0.26, and indeed it is ( λ

ˆ = .27 ). Because the odds ratio 

can be regarded as a property of the conditional distribution of Y2 given Y1

, recovery of the missing values of Y1should

λ

ˆA ˆB|A

provide little additional information about this parameter, so we would expect to be much smaller than λ . The

difference δ , however, is a function of the Y1

and Y2 marginal distributions, so we might expect the overall rate of missing

information for this parameter to be something like the proportion of units for which Y1

and Y2are both missing, plus onehalf the proportion of units for which one of the two measurements is missing. Indeed, the estimate ( λ

ˆ = 0.21) is very close 

to .15 + (.05+.06)/2. If Y1 were filled in, Y2 would still be missing for 21% of the units, and indeed λ

ˆ = .13 is close to onehalf of .21. 

Discussion 

The simple example presented here is not intended to represent a typical application of two-stage MI; we provide it merely to 

illustrate the method and to demonstrate the reasonableness of its results. The method was originally proposed by Shen 

A B A (2000) for computational efficiency when imputation of Y mis is costly but Y mis given Y mis is cheap. Situations of that type

do exist (e.g. Rubin, 2003). However, we believe that the most valuable contribution of two-stage MI is that opens up many

new possibilities for interesting and innovative data analysis. Some of these possibilities are listed below. 

Page 95 
 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Questionnaires with planned missingness (Raghunathan and Grizzle, 1995) and longitudinal data-collection schedules with 

intentionally dropped occasions may reduce respondent burden and can be very cost-effective; with two-stage MI, one can 

isolate the effects of planned and unplanned missing values to inform decisions about the design of current and future studies.

Analysis of longitudinal data with nonignorably missing values has received great attention in recent years (e.g. Little, 1995). 

In real applications, it is often reasonable to view the missing values as a mixture, with some missing for reasons closely

related to the phenomena being measured and others missed for unrelated reasons. With two-stage MI, we can potentially 

separate the effects of the two types of missing values. As shown by Harel (2003), however, what it means for some values to 

be ignorably missing and others to be nonignorably missing needs to be made precise.

Another important set of applications pertains to data with both measurement errors and missing values. Viewing observed

items as imperfect measures of unobservable or latent true scores has been fundamental to the social sciences for decades,

and is the central idea of item-response theory, factor analysis and latent-class modeling. Despite these models' long history, 

the role of missing values in the manifest items remains largely unexplored. We believe that two-stage MI will prove useful 

for investigating not only the relative contributions of response error and nonresponse to overall uncertainty, but how the 

nonresponse and measurement processes interact. 

Finally, two-stage MI can also be used to measure the expected increase in information if a single missing datum or a group

of missing values were suddenly observed. If this type of analysis were implemented in a survey during the data-gathering

process, data collectors could identify the nonrespondents who, if converted, would yield the greatest increase in information 

about various population quantities of interest, and then focus their efforts for nonresponse followup on these units. 

References 

Agresti, A. (1990) Categorical Data Analysis. J. Wiley & Sons, New York. 

Chen, T.T. and Fienberg, S.E. (1974) Two-dimensional contingency tables with both completely and partially cross-classified 

data. Biometrics, 32, 133-144. 

Harel, O. (2003) Strategies for data analysis with two types of missing values. Ph.D dissertation, Department of Statistics,

Pennsylvania State University, University Park, PA. 

Horton, N.J. and Lipsitz, S.R. (2001). Multiple imputation in practice: comparison of software packages for regression

models with missing variables. American Statistician, 55, 244-254.

Kadane, J.B. (1985) Is victimization chronic? A Bayesian analysis of multinomial missing data. Journal of Econometrics, 

29, 47-67, correction 35, 393.

Little, R.J.A. (1995) Modeling the drop-out mechanism in repeated-measure studies. Journal of the American Statistical

Association, 90, 1112-1121. 

Little, R.J.A. and Rubin, D.B. (1987) Statistical Analysis with Missing Data. J. Wiley & Sons, New York. 

Raghunathan, E.R. and Grizzle, J.E. (1995). A split questionnaire survey design. Journal of the American Statistical

Association, 90, 54-63. 

Rubin, D.B. (1987) Multiple Imputation for Nonresponse in Surveys. J. Wiley & Sons, New York. 

Rubin, D.B. (2003) Nested multiple imputation of NMES via partially incompatible MCMC. Statistica Neerlandica, 57, 3-

18. 

Schafer, J.L. (1997) Analysis of Incomplete Multivariate Data. Chapman & Hall, London.

Shen, Z.J. (2000). Nested multiple imputation. Ph.D dissertation, Department of Statistics, Harvard University, Cambridge,

MA.

Page 96