Motivating Example

Two primary dangers exist in the analysis of incomplete data: ignoring missing values with casewise deletion, and applying the wrong method for imputation. Both of these lead to biased coefficient estimates and misleading results. A study of the Special Migration Statistics from the 1991 British Census by Rees and Duke-Williams illustrates both of these dangers clearly.Footnote 7 The authors are concerned that the deliberate embargo of some values in the data to protect the privacy of respondents harms the accurate estimate of patterns of migration flows. These researchers were able to fill in much, but not all, of the missing data with a method called logical data patching. This approach uses observed structure in the data to deterministically fill in missing values. As a simple example, suppose that the response to a survey question asking whether or not the respondent was college educated is missing, but the respondent answered ‘solicitor’ to the occupation question. The analyst can then with high confidence patch in ‘yes’ for the missing value. Rees and Duke-Williams focus on migration between local government districts within Great Britain, and are missing the destinations for some individual cases.Footnote 8 Yet they were able to deterministically obtain many of these values by using sums of migrations in the margins of tables, since columns and rows must sum to the marginals. The rest of the data was filled in with a technique that lowers the level of measurement for variables with missing values and, therefore, will not appeal broadly to empirical political scientists. After these methods are applied, they find that ignoring the plentiful missing information about the economic position of the family head would have seriously biased coefficient estimates in subsequent model fits.

Logical data patching is an effective means of filling missing data since it is done deterministically andconsequently introduces no uncertainty in the form of variance. The two problems with this approach are that it consumes experts’ time, and that it usually cannot fill in all the missing data. In the case of the migration study, the economic position of the family head is a nominal-measured variable (Self-employed, Other employed, Unemployed, Retired, Students (inactive), Other inactive, Students (active)), and most methods of imputing such missing data assume interval or near-interval measured data. Even highly sophisticated implementations of missing data do not handle such cases well, including elaborate dummy-variable set-ups and approximation methods. Thus, there is a need for some method focused specifically on missing categorical data.

Bias from Missing Values

Until recently, the most common method of addressing missing data in empirical models has been to delete every case with even a single missing value: a technique called casewise deletion (also called list-wise deletion or – optimistically – complete-case analysis). With casewise deletion, rows of the data matrix corresponding to individual cases are removed from the dataset, which reduces the effective sample size and lowers the degrees of freedom for any subsequent modelling exercise. For example, if some survey respondents did not answer all of the items on a questionnaire, then those respondents would be effectively deleted, leaving a dataset of only respondents who answered every single question. Unfortunately, this approach is the default with many statistical packages, and there is typically not even a warning message. The obvious convenience of this approach, however, is that the subsequent process of statistical analysis is then unhampered by the missing values.

Generally, casewise deletion induces bias with political data and cannot, therefore, be recommended when analysing data with missing values. The extent of the bias induced by casewise deletion depends on what proportion of the data are missing, how far the pattern of missing values is from being completely random throughout the data matrix, and how different the missing data are from the observed data.Footnote 9 In the special circumstance where missing values occur completely randomly, casewise deletion is unbiased in expectation regardless of what proportion of the data is missing, although small sample issues may emerge in columns of the data matrix. This bias from casewise deletion in real applications ranges from mild to catastrophic.Footnote 10 Replications of published work in political science have shown that the bias induced by casewise deletion can be severe enough to produce faulty inferences.Footnote 11 It is difficult to imagine a situation in which casewise deletion is an ideal treatment of missing values, because it is unbiased only for a narrow definition of data, and always produces greater inefficiency since some fraction of the data are removed. Therefore, a consensus has developed in the statistical, econometric and methodological literatures that missing values should be filled in (imputed) using some principled method if possible.

There are three basic mechanisms by which data can be missing:Footnote 12 it can be missing completely at random (MCAR), missing at random (MAR) or non-ignorable (NI). Each missing- data generating mechanism is described by a simplification of p(M|D), where D is the data matrix (including all explanatory variables as well as the outcome variable) and M is a dichotomous indicator matrix of the same dimensions as D, with values of 1 if a datum is missing and values of 0 if it is present.Footnote 13 Unfortunately, it is not possible to test directly for any of these mechanisms; one of them must be assumed. So researchers with this problem need to pay attention to the substantive nature of how the data became missing. For example, it is well-known that survey respondents are less likely to answer the personal income question if their income is very high or very low. This not only means that casewise deletion gives an inaccurate view of the subjects’ income, it also means that characteristics of interest related to these two groups will be censored from the data as well.

The best-case scenario for casewise deletion is when observations are missing completely at random (MCAR). MCAR means that there is no underlying associative process that causes the absence of data: the occurrence of missing values is uniformly random throughout the data matrix. When the missing values are orthogonal (mathematically unrelated) to any values of the data, missing or observed, the missing data are MCAR: p(M|D) = p(M). A simple way to think of MCAR missing values is that both the observed and missing data, independently, have the properties of a random sample of the population. More specifically, the data for which we observe responses, D_R, and the data for which we do not have responses, D_NR, have the same underlying distribution. Though unbiased, casewise deletion in this circumstance may not necessarily be ideal since information is lost by discarding incomplete observations and, even if no bias ensues, precision is lost and standard errors from regression-style models will be artificially large.

Unfortunately, data where missing values are completely random are not particularly common in political science: survey respondents are more sensitive to some questions than others, certain types of regimes are more likely to refuse to provide economic or political data, some information is more difficult to recall, and so on. This is why the term ‘independent variables’ is a poor description for the social sciences. It is difficult to think of a situation in political science, other than a computer malfunction, that would result in missing values being entirely unrelated to any attribute or political phenomena, observed or unobserved.

More commonly, missing data is assumed to be missing at random (MAR). If missing values are MAR, then their occurrence is related to the observed data D_R, but not to any unobserved data D_NR. The MAR assumption is expressed formally as p(M|D) = p(M|D_R). MAR is a fairly intuitive concept: it simply means that missing values in one variable, say income, can be related to other variables (education, occupation, neighbourhood and so on) but those other variables must be recorded in the dataset.

When missing values are MAR, casewise deletion will always result in bias. This is because the missing data, D_NR, are not a random subsample of all the data D. The bias induced by casewise deletion under MAR is similar to sampling bias in survey research: the sample of values that are recorded, D_R, is biased away from subjects with characteristics that make them less likely to respond and, consequentially, the recorded data D_R are not representative of the full sample, D. Multiple imputation techniques, conversely, are ideal for dealing with continuously measured MAR data since they produce efficient and unbiased estimators with relatively little effort.

Table 1 Assumptions about Missing Data

Finally, missing data can be non-ignorable (NI). Non-ignorable missing data occur when missing values are related to unknown and/or unobserved parameters. Consider again the example of missing values in an income question on a survey, but now we do not observe the variables of education, occupation and neighbourhood that can help fill in the missing income data. In such a situation, the missing data are NI. Formally, when NI holds, the expression for the missing data mechanism, p(M|D), cannot be simplified and specific information is fully missing from the dataset. With NI missing data, both casewise deletion and imputation tools produce biased statistical results since the missing data is purposefully different from the observed data. Non-ignorable missing data also biases subsequent modelling with multiple imputation since there is no observed information from which to build an imputation process for filling in missing values.

One difficulty with non-ignorable missing values is that the condition of non-ignorability is, by definition, impossible to test for and analysts will never know with certainty whether missing values in their dataset are non-ignorable.Footnote 14 Whether missing values are treated as MCAR, MAR, or NI comes strictly from an assumption researchers make based on their understanding of the data generation and data missing process. The required assumption about the structure of missing values has fundamental implications for the modelling process: assuming MCAR or MAR allows the researcher to handle the missing data in a principled manner, but, if this assumption is incorrect, then subsequent inferences are likely to be misguided. However, recognizing NI missing data means that no standard statistical method exists for repairing the problem. An obvious course of action is to gather further data in an attempt to create conditional information such that the missing values are now MAR, but this is often not possible. If the researcher has a well-developed understanding of the response mechanism, Bayesian model-based imputation through a Markov chain Monte Carlo process may be helpful because assumptions are made about relationships in the data and information may be provided through the prior distribution.Footnote 15

The three patterns of missing data can be further illustrated by looking at the development of a likelihood function for the data and an arbitrary parameter vector $$\[-->$<>{{{\theta }}<$> <!--\]$$ to be estimated. Segment the data matrix, X, by rows into two constituent parts: X = [X_R, X_NR] for responses and non-responses, and write the likelihood:

$${{\bi L}}({{\bi \bitheta }}&#x007C;{\bf {{X}}}) = f({{{\bf {{X}}}}_R},\,{{{\bf {{X}}}}_{NR}}&#x007C;{{\bi \bitheta }}) = f({{{\bf {{X}}}}_R}&#x007C;{{\bi \bitheta }})f({{{\bf {{X}}}}_{NR}}&#x007C;{{{\bf {{X}}}}_R},\,{{\bi \bitheta }})\eqno\$$

The second equality comes from the standard definition of conditional probability. If we use the standard logarithmic transformation of the likelihood function, then the product above becomes a sum:

$$\ell ({{\bi \bitheta }}&#x007C;{\bf {{X}}}) = \ell ({{\bi \bitheta }}&#x007C;{{{\bf {{X}}}}_R},\,{{{\bf {{X}}}}_{NR}}) = \log f({{{\bf {{X}}}}_R}&#x007C;{{\bi \bitheta }}) + \log f({{{\bf {{X}}}}_{NR}}&#x007C;{{{\bf {{X}}}}_R},\,{{\bi \bitheta }}).\eqno\$$

We first notice that casewise deletion is biased whenever logf(X_NR|X_R, θ) is non-zero since then $$\[-->$<>\ell ({{\bi \bitheta }}&#x007C;{\bf {{X}}})\ne \ell ({{\bi \bitheta }}&#x007C;{{{\bf {{X}}}}_R}) <$> <!--\]$$ . If the data are MCAR, then there is no conditional relationship between X_NR and X_R or other missing data. Therefore, the second term is not material and casewise deletion is acceptable. However, if there is a conditional relationship between missing data and observed data (MAR) then log f (X_NR|X_R, θ) is clearly non-zero, but a tool that informs the likelihood using this conditionality can overcome the problem. This is what imputation tools do: impute missing values conditionally based on their relationship with observed values. In the case of NI data, the second term in Equation 1 is actually log f (X_NR|X_ONR, X_R, θ), where X_ONR is ‘other non-responses’. In this case, the term here is non-zero and not estimable with the resources at hand. Therefore, assumptions about how to treat the missing data are actually assumptions about the relationship between X_NR and other data elements, observed or missing.

The Origins of Hot Deck Imputation

Because of the many problems associated with casewise deletion, academic and professional survey researchers began looking for a better way to analyse incomplete data in the 1960s. An early remedy developed roughly in the 1970s by market researchers and census takers was called ‘hot decking’ as a literal reference to taking draws from a deck of computer punch-cards.Footnote 16 In hot deck imputation, if a respondent has a missing value, a set of ‘similar’ respondents are pulled by hand from the data and one of them is randomly selected to provide a fill-in value for the missing datum. So there are two steps to the process: human selection of a sub-sample of cases, and uniform random selection from within this sub-sample. An obvious challenge is that it is not always clear what is meant by ‘similar’, although practitioners of the time became quite accomplished at performing this task over regular iterations. Hot deck imputation methods subsequently grew from simple random sampling to more complicated algorithms in attempts to find respondents as similar as possible to those with missing responses.

Hot deck imputation preserves the integrity of the data to the extent that only observed values of the variable being imputed are candidates for imputation. This, in turn, means that it is impossible for hot deck imputation to impute values outside the observed range of the data, a reasonable limit on how extremes of values can be imputed. Also, within the observed range of the data, only values that have occurred for other cases can be imputed, thus ensuring that nonsensical values cannot be imputed. The technique works best for discrete data because, while it is technically impossible that two respondents will have the exact same value of a continuous variable in the absence of rounding, it is quite likely that several respondents will have the same value for discrete variables (such as gender).

The major problem with traditional hot decking is that it does not reflect uncertainty in the imputed values: since only one value is imputed for each missing datum, all imputations are subsequently treated as factual responses rather than subjective probabilistic imputations. So there is no direct way to account for the uncertainty resulting from the process of sub-sample specification and then case selection from this sub-sample. Therefore, the subsequent model fit produces smaller standard errors and tighter confidence bounds than is appropriate.Footnote 17 In addition, the magnitude of this underestimation increases as the amount of missing data increases and as the quality of the hot deck matches decreases. We address this problem in the section ‘ Multiple Hot Deck Imputation’.

Parametric Multiple Imputation and Its Shortcomings

A crucial improvement in the empirical analysis of incomplete data was precipitated by Rubin when he developed multiple imputation.Footnote 18 Rubin and his collaborators further refined the technique over the next two decades,Footnote 19 and it is now the dominant strategy for handling missing data in the social sciences and elsewhere for those not casewise-deleting.

With multiple imputation each missing datum is replaced by several (M≥2) imputed values from a conditional distribution using other present values for that case. Thus, several imputed datasets are created with the appropriate missing datum being imputed by one of the M imputed values to produce M complete datasets. Then the complete-data statistic, $$\[-->$<>{{{\theta }}<$> <!--\]$$ (for example, a regression coefficient), is estimated M times and the $$\[-->$<>{{{{\bi \hat{\theta }}}}_m} <$> <!--\]$$ , m = 1,…M statistics are averaged to create a single estimate of $$\[-->$<>{{\bar{\theta }}_{{\bi M}}} <$> <!--\]$$ .Footnote 20 A graphical overview of the process is shown in Figure 1.

Fig. 1 Multiple imputation Note: This figure illustrates the method of multiple imputation by which a dataset is copied several times, imputations are drawn for the missing values in each of the datasets, a statistic of interest is computed for each of the datasets, and, finally, those statistics are combined into a single result.

While parametric multiple imputation remains the state-of-the-art technique for filling in missing observations, it can be problematic when used on discrete data. Multiple imputation, as commonly implemented, assumes a continuous metric for the missing data and, therefore, produces imputations that are continuous. The multiple imputation process then transforms categorical definitions into continuous values.Footnote 21 In this way a strictly binary variable for gender might be imputed with a value of 0.6 with parametric multiple imputation. Such problems intensify as the granularity of the discrete data declines. Alternatively, while binary response models will not fail if a nonsensical value is imputed into an explanatory variable, it may produce impossible or nonsensical results. A more basic problem is that such imputations fail to respect the meaning of the measure they are imputing; even if no bias is introduced, they degrade the substantive meaning of a variable of interest.

The most common solution to the categorical measurement problem with parametric multiple imputation has been to round continuous imputations to their nearest discrete value.Footnote 22 This approach allows for the use of multiple imputation and its attractive theoretical properties, while allowing the user to use complete data techniques appropriate for their discrete data. King and his coauthors have implemented this approach in their easy-to-use multiple imputation program called ‘Amelia’,Footnote 23 which has become a popular implementation of multiple imputation in political science.Footnote 24 Since models produced from parametric multiple imputation produce unbiased and consistent coefficient vectors, the acceptance of multiple imputation in the field has been a major improvement over previous approaches and an informal consensus has emerged holding that the process of rounding continuous imputations is not problematic.

Problems with the rounding approach, though subtle, are non-trivial. Chief among the pitfalls of this approach is that it produces biased results in nearly every case. Consider the imputation of missing values in a binary variable using parametric multiple imputation and rounding. If an imputation of 0.3 is made, it is rounded to 0. The imputation model however, did not predict that the value should be zero, it predicted that the value should be 0.3. Therefore, we know that parametric multiple imputation produces unbiased imputations for MAR data and we still deliberately bias the results by rounding with the entirely arbitrary cutpoint of 0.5. As we continue to draw imputations closer to zero than one and round them to zero, we are biasing our imputations further and further from what ought to be the centre of their distribution. The parametric multiple imputation with rounding technique will only be unbiased in the unlikely situation of imputations having a unimodal symmetric distribution centred exactly on 0.5. In all other cases, using parametric multiple imputation with rounding on discrete data will necessarily result in some degree of bias. Furthermore, coefficients from parametric multiple imputation and rounding on discrete data will have inappropriate standard errors because the distance between the imputation and the nearest integer is lost. Further still, in multinomial models when no latent variable may be assumed, an imputed and rounded continuous value can lead to substantially different results because of rounding error.

Affinity Scoring

The selection of donor values is critically important for the validity of the values to be imputed. Donor values should be drawn from respondents (observations) as similar as possible to the respondent with a missing value.Footnote 29 In order to measure the extent to which a respondent with a missing value is similar to other respondents, we create a set of affinity scores bounded by 0 and 1, whose components are denoted α_i,j and measure the degree of similarity recipient i has to each potential donor j. Affinity is defined in terms of the degree to which each potential donor matches the recipient's values across all variables other than the one being imputed.

For each respondent we have the vector (y_i, x_i) where y_i indicates the outcome variable and x_i is a k-length vector of purely discrete explanatory variables (an unrealistic assumption that we will relax later), either of which may contain missing values. If the ith case under consideration has q_i missing values in x_i, then a potential donor vector, x_j, j ≠ i, will have between 0 and k – q_i exact matches with i.

Now define z_i,j as the number of variables for which the potential donor j and the recipient i have different values. Thus, k – q_i – z_i,j is the number of variables on which j and i are perfectly matched. This value, scaled by the highest number of possible matches (k – q_i) is then the affinity score:

$${{\alpha }_{i,j}} = \frac{{k{\rm{ - }}{{q}_i}{\rm{ - }}{{z}_{i,j}}}}{{k{\rm{ - }}{{q}_i}}}.\eqno\$$

Notice that the affinity score as defined in Equation 3 has the desirable properties that $$\[-->$<> {{\alpha }_{i,j}} = 1 <$> <!--\]$$ for $$\[-->$<> i\,{{{\in }\,{{{\bf {{D}}}}_R} <$> <!--\]$$ (data with responses) and $$\[-->$<> {{\alpha }_{i,j}} = 0 <$> <!--\]$$ for $$\[-->$<> i \in {{{\bf {{D}}}}_{NR}} <$> <!--\]$$ (data missing responses). Now write the vector of all α_i,j as α_i. Table 3 shows how the value of α_i,j decreases as the number of matches decrease. Cases where the recipient and the donor are both missing values in the same covariate are deducted from k and q_i prior to the calculation of Equation 3.

Table 3 Affinity Scoring

Note: As the number of matches declines, the affinity score moves towards zero.

As a simple example, consider the illustrative data in Table 4. For the single missing value of x ₁ at observation 6, a vector of affinity scores is computed using the formula in Equation 3. Observations that match observation 6 perfectly on the other x variables get affinity scores of 1. For observations that match the recipient observation (observation 6) on one variable but not the other, α_i,j = (3–2)/(3–1) = 0.5. Lastly, those observations that do not match the recipient on the other x values get an affinity score of 0.

Table 4 Affinity Scoring Example

Note: A vector of affinity score is computed for observation i=6, which has only one missing value denoted by R language notation for missing values: NA. Three observations with perfect matches on the X variables get a sore of 1, three that match on half get 0.5 and those with no matches get 0. The affinity score for i=6 (denoted ‘–’ ) is not a missing element in the α vector, but rather not an element of that vector at all.

When variables are continuous (or discrete with many categories such that their structure begins to resemble that of continuous variables), the probability of a matched value between a potential donor and a recipient goes to 0 due to the measurement. It is wrong to say that any non-matched values should be treated equally as non-matches in the discrete sense. Potential donors, j, for the hth variable in $$\[-->$<>{{{\rm{x}}}_{j{\rm{[}}h{\rm{]}}}} <$> <!--\]$$ that are ‘close’ to the corresponding $$\[-->$<>{{{\rm{x}}}_{\ell {\rm{[}}h{\rm{]}}}} <$> <!--\]$$ value of the recipient should have greater affinity. This means that there needs to be a definition of ‘close’ here. If a value $$\[-->$<>{{{\rm{x}}}_{j{\rm{[}}h{\rm{]}}}} <$> <!--\]$$ is closer to $$\[-->$<>{{{\rm{x}}}_{i{\rm{[}}h{\rm{]}}}} <$> <!--\]$$ than some other case's value $$\[-->$<>{{{\rm{x}}}_{\ell {\rm{[}}h{\rm{]}}}} <$> <!--\]$$ , then in terms of the variable h, j has greater affinity to i than $$\[-->$<>\ell <$> <!--\]$$ . At the same time, if $$\[-->$<>{{{\rm{x}}}_{j{\rm{[}}h{\rm{]}}}} <$> <!--\]$$ and $$\[-->$<>{{{\rm{x}}}_{\ell {\rm{[}}h{\rm{]}}}} <$> <!--\]$$ are very close to one another, but $$\[-->$<>{{{\rm{x}}}_{j{\rm{[}}h{\rm{]}}}} <$> <!--\]$$ is somewhat closer to $$\[-->$<>{{{\rm{x}}}_{i{\rm{[}}h{\rm{]}}}} <$> <!--\]$$ , it does make sense to say that the two values are so disparate that they should carry very different affinity values. To allow our affinity scoring technique to address these circumstances, we extend the technique in the following way: for a continuous variable in x, we count $$\[-->$<>{{{\rm{x}}}_{j{\rm{[}}h{\rm{]}}}} <$> <!--\]$$ as matching $$\[-->$<>{{{\rm{x}}}_{i{\rm{[}}h{\rm{]}}}} <$> <!--\]$$ if $$\[-->$<>{{{\rm{x}}}_{i{\rm{[}}h{\rm{]}}}} <$> <!--\]$$ and $$\[-->$<>{{{\rm{x}}}_{j{\rm{[}}h{\rm{]}}}} <$> <!--\]$$ are in the same concentric standard deviation from the mean of variable h. Values outside of this range are penalized in affinity score by distance. These decisions are the defaults in our R package (hot.deck), but they may be changed by users.

Another option for dealing with continuous variables is to use a metric that measures the distance of observation i and observation j over several x variables. An example of such a technique would be the vector norm || x || = $$\[-->$<> {{{\rm{(\bf x}}^{\prime}{\bf V}{\rm{\bf x)}}}^{{\rm{1/2}}}} <$> <!--\]$$ , with positive definite weight matrix V. One may then define $$\[-->$<>\left\Vert {{{{\rm{x}}}_i}{\rm{ - }}{{{\rm{x}}}_j}} \right\Vert <$> <!--\]$$ to be the distance between the vectors x_i and x_j, where different standard choices for the matrix V are possible. For example, by putting the inverse of the of the sample variances of x (in each cell) on the main diagonal of V or the Mahalanobis distance metric (the inverse of the sample variance-covariance matrix of x (in each cell)). The drawback of such approaches is that, while they provide information about the distance between values of x, whether a given distance would count as a match or not when computing the affinity scores would still depend on arbitrary cutoffs.

Our nonparametric affinity score is based on empirical description, and many other variants exist in various fields. The general affinity score is only one type of similarity measure, and the full class has many uses. Much of the work in this area is done in biology where the concern focuses on similarity and dissimilarity of species types (taxonomy). There is a set of classic alternatives measures in this area including:

• • Czekanowski's Quantitative Index $$\[-->$<> (C{{Z}_{ik}} = 2\sum _{{j = 1}}^{n} \,\min ({{x}_{ij}},{{x}_{kj}})/\sum _{{j = 1}}^{n} ({{x}_{ij}} + {{x}_{ik}})) <$> <!--\]$$ ,

• • Morisita's Index $$\[-->$<> ({{M}_{ik}} = 2\sum _{{j = 1}}^{n} ({{x}_{ij}}{{x}_{kj}})\,/\,(\sum _{{j = 1}}^{n} x_{{ij}}^{2} + \sum _{{j = 1}}^{n} x_{{kj}}^{2} )) <$> <!--\]$$ ,

• • Canberra Metric $$\[-->$<> ({{C}_{ik}} = \sum _{{j = 1}}^{n} [({{x}_{ij}}{\rm{ - }}{{x}_{kj}})\,/\,({{x}_{ij}} + {{x}_{kj}})]\,/\,n) <$> <!--\]$$ ,

matrix-based scores, and more.Footnote 30 The tradeoffs between these alternatives are usually application specific. This is also an ongoing research area in information theory where research is centred on comparing text or electronic messages, as well as image restoration. In recent political science methodological research there has been substantial work on causal inference. In this area propensity scores are a common tool for matching cases in observational data as a means of asserting pseudo-experimental control over the assignment of treatment to cases.Footnote 31 The propensity score is the probability that a case receives the treatment effect, relative to the control, given this case's covariate levels. Thus, propensity scores are analogous to our affinity scores since they are also a single number summary for each case in the data, bounded by [0:1], and reflect unit-wise similarity. They differ in that propensity scores are an n-length vector relative to a single criteria (treatment versus control), whereas affinity scores are an n×n lower-triangular matrix of cross-relations between all cases in the data. Propensity scores also differ substantially from our nonparametric affinity score in that they are model based where a statistical specification (usually a logit/probit model) is applied to produce the assignment probabilities for the cases. A critical advantage to our approach, over such strategies, is that a regression model specification, with all of the concerns about covariate and functional form choices, is not necessary.

Imputation Using Affinity Scores

From the vector of affinity scores, our multiple hot deck imputation method is quite intuitive and possesses attractive theoretical properties. Specifically, within the best imputation cell (the set of donors with the highest affinity scores), the realization of missing values is independent of the data. Unless the missing data is non-ignorable, it must be the case that missing values within this best cell are missing completely at random (MCAR). Therefore, we can draw randomly from the best cell to fill in missing values without introducing bias. Also, we expect random draws from the best cell to be high quality imputations; because all donors in the cell have the highest possible affinity with the recipient, not only must the occurrence of missing values be completely random within this cell, but so must be the values of the variable itself. In Appendix A, we describe a hybrid hot deck-parametric algorithm to handle datasets that have a mixture of discrete and continuous variables with missing values.

Consider a dataset D, with row vectors (y_i, x_i) and missing values in some or all of its variables, where D_R continues to denote observed responses and D_NR continues to denote non-responses. For some variable h, the value of x_{i [h ]} is missing. Intuitively, we could not produce reliable imputations for the missing x_{i [h ]} by randomly drawing from the vector x_{[h ]} in D_R since realizations of x_{[h ]} may be related to the design of the data. But it is possible to subset the observed data (and thus the x_{[h ]} in D_R) in such a way that the realization of x_{[h ]} is unrelated to the design.

Each respondent with a missing x_{i [h ]} can be considered a member of an imputation cell. Imputation cells (sometimes called estimation cells) are subsets of the sample to which respondents are assigned based on characteristics. For instance, if the characteristics of interest are gender and voter turnout, we would have four imputation cells: male voters, male non-voters, female voters, and female non-voters. Assume, for purposes of exposition, that there are an arbitrary number of respondents in D_R who match recipient i perfectly on all x_{i [−h ]} variables.

If the imputation cell, C, contains the set of best possible donors (those with the highest affinity scores), random draws from the observed x_{j [h ]} in the cell will be unbiased and have the least amount of imputation variance theoretically supportable. Formally, the observations composing best imputation cell, C, are those for which α_i,j = max α_i for all j in the cell C. In some cases this set will consist of only perfect matches and in others it will simply consist of the best matches. Returning to the example from Table 4 for intuition, the best imputation cell would consist of only those observations with the highest (in this case, 1) affinity score. So, observations 2, 3, 5 and 10 would be chosen for the imputation cell. Their values on x₁ are respectively: {1, 0, 1, 1}.

By subsetting respondents appropriately into imputation cells, we are able to treat the responses in the best imputation cell C as the realization of independent and identically distributed (iid) random variables having mean μ _C and variance $$\[-->$<>\sigma _{{\rm{C}}}^{2} <$> <!--\]$$ . Now,

$${{{\rm{x}}}_{i[h]}}\,{\rm{\sim }}\,f({{\mu }_{\rm{C}}},\,\sigma _{{\rm{C}}}^{2} ),\,{\rm{for}} \ {\rm{all}} \ i \ {\rm{in}} \ {\bf {{C}}}{\rm{.}}\eqno\$$

So, all values of x_{i [h ]} in the cell will be drawn from the same distribution, whatever that distribution turns out to be.

More importantly for our purposes, dividing the respondents into cells independent of the sampling and response mechanisms implies that the distribution in Equation 4 holds for both observed and unobserved values alike.Footnote 32 Formally,

$$ {{{\rm{x}}}_{i[h]}}\,&#x007C;\,({\bf {{D}}}{\rm{,}}\,{{{\bf {{D}}}}_R},\,{{{\bf {{D}}}}_{NR}})\,{\rm{\sim }}\,f({{\mu }_{\rm{C}}},\,\sigma _{{\rm{C}}}^{2} ), \ {\rm{for}} \ {\rm{all}} \ i \ {\rm{in}} \ {\bf {{C}}}{\rm{.}} \eqno\$$

This shows that the distribution of x_{i [h ]} in the imputation cell is the same regardless of whether it is conditioned on all the data D, the recorded data D_R, or the missing data D_NR. Even though the values of x_{i [h ]} may be related to the data in the population of interest and the sample, the division into imputation cells removes that dependence within the cell.Footnote 33 Therefore, the missing values must be missing completely at random (MCAR) within the imputation cell.

Since the missing values within the best imputation cell are now MCAR and there exist realisations of x_{[h ]} that are independent of the design, we can draw randomly from the (non-missing) values of x_{[h ]} to impute the missing value x_{i [h ]}. We still need multiple imputations in order to capture the imputation variance. The next step is to draw with replacement M ≥ 2 values of x_{[h ]}, then assign the mth draw to the mth duplicate dataset as x_{i [h ]} for D_m, m = 1, …, M. This process is repeated for every missing value in every column of the dataset. The result is M imputed datasets each with no missing values.

Estimation after Multiple Hot Deck Imputation

The final step, as outlined in the original work on multiple imputation by Rubin et al., is to run M copies of the model specification, taking the mean of the coefficient estimates and the weighted mean of the standard errors, where the weighting accounts for between and within variation of the imputations.Footnote 34 Let $$\[-->$<>{{\mathop{\hat \theta }}_m} <$> <!--\]$$ , m = 1,…,M be coefficient estimates computed individually from the M imputed datasets, and $$\[-->$<>{{\sum }_m} <$> <!--\]$$ , m = 1,…,M be the associated variances for $$\[-->$<>{{\mathop{\hat \theta }_m} <$> <!--\]$$ . A single estimate of $$\[-->$<>{{{\theta }}<$> <!--\]$$ ,

$$ {{\bar{\theta }}_{{\bi M}}} = \frac{1}{M}\sum _{{m = 1}}^{M} {{\mathop{\hat\theta }}_m}, \eqno\$$

is created from the mean of $$\[-->$<>{{\mathop{\hat \theta }}_m} <$> <!--\]$$ over the m values.

Calculating the variability of the estimate of $$\[-->$<>{{{\theta }}<$> <!--\]$$ is slightly more involved than the mean used to produce $$\[-->$<>{{\bar{\theta }}_M} <$> <!--\]$$ since the total variance is composed of the variation of the coefficient estimates within each imputed dataset and the variation of the coefficient estimates between the imputed datasets. The within imputation variance is the mean of individual coefficient variances across models:

$$ {{W}_M} = \frac{1}{M}\sum _{{m = 1}}^{M} {{\rSigma }_m}, \eqno\$$

The between imputation variance is the variance of the M coefficient estimates:

$$ {{B}_M} = \frac{1}{{M{\rm{ - 1}}}}\,\sum _{{m = 1}}^{M} {{({{\mathop{\theta }\limits^{\frown}}_m}{\rm{ - }}{{\bar{\theta }}_{{\bi M}}})}^2} . \eqno\$$

The total variance for $$\[-->$<>{{\bar{\theta }}_M} <$> <!--\]$$ is an ANOVA-style weighted sum:

$$ {{T}_M} = {{W}_M} + \left( {1 + \frac{1}{M}} \right){{B}_M}. \eqno\$$

where [1+1/M] is the adjustment for finite M, and the new degrees of freedom are:Footnote 35

$$d{{f}_{\!\!M1}} = (M{\rm{ - 1}})\left[ {1 + \frac{1}{{M + 1}}\,\,\,\frac{{{{W}_M}}}{{{{B}_M}}}} \right].\eqno\$$

These quantities are produced automatically for the user by statistical software implementing multiple imputation.Footnote 36

Monte Carlo Experiments

Researchers in political science with incomplete data generally either casewise delete or use software for parametric multiple imputation. The two most common parametric MI solutions are the R packages ‘Amelia’Footnote 37 and ‘mice’.Footnote 38 Other approaches exist, but are not often considered by researchers in the discipline, including specification of the missing data into a Bayesian model as (nuisance) parameters to be estimated, or application of EM.Footnote 39 Both of these are model-specific approaches, meaning that they must be uniquely specified for each application, and are, therefore, not directly comparable to imputation tools.

To evaluate the performance of multiple hot deck imputation relative to the most commonly used missing data techniques, we provide a series of Monte Carlo experiments. For all of the experiments, we randomly generated datasets of 500 observations, each with five binary variables with fairly high correlations between them produced by a multivariate normal, $$\[-->$<>{\cal N} <$> <!--\]$$ (μ, Σ), with:

$${\bf {{\bimu }}}\,{\rm{ = }}\,{\rm{[0,0,0,0,0],}}\,\,\,\,{\bf {{\brSigma}}}\,{\rm{ = }}\,\left[ {\matrix{ {\rm{1}} &amp; {{\rm{0}}{\rm{.8}}} &amp; {{\rm{0}}{\rm{.8}}} &amp; {{\rm{0}}{\rm{.8}}} &amp; {{\rm{0}}{\rm{.8}}} \cr {{\rm{0}}{\rm{.8}}} &amp; {\rm{1}} &amp; {{\rm{0}}{\rm{.8}}} &amp; {{\rm{0}}{\rm{.8}}} &amp; {{\rm{0}}{\rm{.8}}} \cr {{\rm{0}}{\rm{.8}}} &amp; {{\rm{0}}{\rm{.8}}} &amp; {\rm{1}} &amp; {{\rm{0}}{\rm{.8}}} &amp; {{\rm{0}}{\rm{.8}}} \cr {{\rm{0}}{\rm{.8}}} &amp; {{\rm{0}}{\rm{.8}}} &amp; {{\rm{0}}{\rm{.8}}} &amp; {\rm{1}} &amp; {{\rm{0}}{\rm{.8}}} \cr {{\rm{0}}{\rm{.8}}} &amp; {{\rm{0}}{\rm{.8}}} &amp; {{\rm{0}}{\rm{.8}}} &amp; {{\rm{0}}{\rm{.8}}} &amp; {\rm{1}} }} \right]\eqno\$$

Binary values were created by rounding the draws to 0 or 1 depending on whether they were below or above 0, respectively. This data setup maximally advantages parametric multiple imputation and maximally handicaps our multiple hot deck imputation technique, since it is based originally on a continuous parametric specification. We introduce missing values into two of the five variables where the values are missing at random (MAR), meaning that the occurrence of missing values is related to the observed variables.Footnote 40 We designed the introduction of MAR missing values to cause the greatest possible bias. This was done at different levels: the probability that a value was missing from the two affected variables, conditional on the observed values of the other variables all being zero, was applied at 20 per cent, 50 per cent and 80 per cent, respectively. These experiments are labelled MAR-1, MAR-2, and MAR-3, and are run for 10,000 Monte Carlo iterations each, and the results are summarized in Table 5.

Table 5 Means and Percentages Correctly Imputed for Different Imputation Techniques in the Monte Carlo Experiments

Note: Standard deviations are in parentheses.

The left portion of Table 5 provides the means and standard deviations of the variables with missing values for: the true data (before missing values were introduced), casewise deleted data, multiply hot decked results, and multiply imputed results (as implemented by Amelia). For each of the three MAR studies, the mean of the true data was 0.5 from the data setup. Casewise deletion consistently results in bias, and it is a bias that increases in intensity as the number of missing values increases: a mean of 0.530 for MAR-1 to a mean of 0.656 for MAR-3). Parametric multiple imputation displays a smaller bias than casewise deletion, and this is a bias that also increases as the amount of missing data increases. Notice that the means from the multiply hot decked data display no bias across these three studies, although there is a slight increase in the standard deviation of the means as the proportion of missing values increases, which indicates increasing uncertainty (a desirable result as the volume of missing values goes up). The right portion of Table 5 shows the percentage of imputed values that were imputed correctly for both approaches using the nave criteria. Both multiple hot decking and parametric multiple imputation perform remarkably well in this regard as more missing values are added: hot deck imputation getting about 83 per cent of the missing values correct and parametric multiple imputation getting about 75 per cent correct.

Figure 2 presents the results of logistic regressions conducted at each iteration of the Monte Carlo studies; for each regression, one of the variables with missing values was the outcome variable and two fully observed variables served as explanatory variables. The figure displays the distribution of point estimates over the 10,000 iterations of the Monte Carlo studies for β ₀ (upper-left cell), β ₁ (lower-left cell) and β ₂ (lower-right cell). Within each cell are three sub-plots, from left to right, for the MAR-1, MAR-2 and MAR-3 separate studies. Within each study, looking from left to right, the first boxplot displays the true values of the point estimates, the second boxplot displays the point estimates under casewise deletion, the third boxplot displays those from multiple hot deck imputation, and the fourth boxplot for parametric multiple imputation.

Fig. 2 Comparison of logistic regression results from Monte Carlo experiments

Figure 2 shows again that casewise deletion results in consistently biased point estimates and that this bias becomes more severe as progressively more missing values are added to the data. As with the mean summary in Table 5, parametric multiple imputation typically displays less bias than casewise deletion, but also a bias that gets worse as the volume of missing values increases. Interestingly, for the estimation of β ₁, parametric multiple imputation produces a more severe bias than casewise deletion. The bias of the point estimates produced by parametric multiple imputation is due to the rounding discussed in the section ‘Parametric Multiple Imputation and Its Shortcomings’. Here, multiple hot deck imputation performs quite well: its median point estimates are consistently closer to the median of the true point estimates than the other techniques. The interquartile range of the multiply hot decked point estimates increases as more missing values are added to the data, which usefully describes increased uncertainty as there is less data on which to condition.