Multiple Imputation Using Gaussian Copulas

Abstract

Missing observations are pervasive throughout empirical research, especially in the social sciences. Despite multiple approaches to dealing adequately with missing data, many scholars still fail to address this vital issue. In this article, we present a simple-to-use method for generating multiple imputations (MIs) using a Gaussian copula. The Gaussian copula for MI allows scholars to attain estimation results that have good coverage and small bias. The use of copulas to model the dependence among variables will enable researchers to construct valid joint distributions of the data, even without knowledge of the actual underlying marginal distributions. MIs are then generated by drawing observations from the resulting posterior joint distribution and replacing the missing values. Using simulated and observational data from published social science research, we compare imputation via Gaussian copulas with two other widely used imputation methods: multiple imputation via chained equations and Amelia II. Our results suggest that the Gaussian copula approach has a slightly smaller bias, higher coverage rates, and narrower confidence intervals compared to the other methods. This is especially true when the variables with missing data are not normally distributed. These results, combined with theoretical guarantees and ease of use, suggest that the approach examined provides an attractive alternative for applied researchers undertaking MIs.

Keywords

missing data Bayesian statistics imputation categorical data estimation

Missing data problems are ubiquitous in observational data and common among social science applications. Statistical inference that does not adequately account for the missing data is widely known to lead to biased results and inflated (or deflated) variance estimates (King et al. 2001; Molenberghs et al. 2014; Rubin 1976; White and Carlin 2010). Even though most statistical software platforms provide methods that adequately handle missing data (the most popular of these is multiple imputations [MI]), they are often ignored by applied researchers.¹

In Figure 1, we illustrate the number of articles published in five top Sociology and Political Science journals since 1990 that contain “multiple imputations” in the body of the article.² Our survey of the literature shows the rapid growth of the use of MIs in the social sciences. Nevertheless, as missing data is a feature of almost any observational data set, the annual counts of articles mentioning MIs per year still point to significant underutilization of this method in the social sciences. This may be due to a lack of understanding of the benefits (and assumptions) of common imputation methods.

Figure 1.

Number of references to “multiple imputation” in articles from five top Sociology and Political Science journals since 1990.

This article has two aims. First, we introduce applied researchers in the social sciences to a specific copula method for imputation and discuss its advantages over other methods. The method discussed is easy to implement using the sbgcop package (Holf 2010) in $R$ (R Development Core Team 2004)³ and has theoretical properties that make it attractive. Second, we conduct a systematic evaluation and comparison of the copula method to two commonly used imputation software packages multiple imputation via chained equations (MICEs; Van Buuren and Groothuis-Oudshoorn 2011) and AMELIA II (Honaker, King, and Blackwell 2012) in sociology and political science.

Copulas are often used for the estimation of dependency between variables and are particularly useful in the generation of imputations as they allow for the construction of valid joint distributions of the data, even if the researcher has little knowledge about the actual joint distribution. Given the joint distribution of the data, we can generate imputations by sampling from the conditional distribution of the missing data given the observed data. We highlight a semiparametric Gaussian copula approach to missing data imputation. The Gaussian copula is one particular way of constructing a joint distribution from which missing values can be easily drawn. The method was initially developed by Hoff (2007) to estimate empirical models on multivariate data.

In particular, the Gaussian copula defines the dependence among the distributions of a set of variables which may contain missing values. These variables can include normal, ordinal, and binary variables. Rather than using the distributions themselves, a rank likelihood approximation is used. As a result, the technique does not require the specification of marginal or conditional distributions. This is in stark contrast to other imputation methods using copulas that either require knowledge of the marginals or correlation structure (Käärik 2006, Käärik and Käärik 2009; Robbins, Ghosh, and Habiger 2013) or target different copula parameters via pseudolikelihood methods (Di Lascio, Giannerini, and Reale 2015). The proposed approach allows applied researchers to undertake imputations of their data without relying on prespecification or ad hoc decisions.

The potential use of copulas for MI applications has not been thoroughly discussed within the social sciences. The copula methods we describe are easy to use and are more likely to provide a good representation of the joint distribution of the data than existing methods. Moreover, provided the Markov chain Monte Carlo (MCMC) converges, the output from the copula model represents a valid posterior density. Simply put, this means that we have theoretical guarantees about the posterior distribution from which the imputations are generated that other methods can not provide. Based on an extensive simulation exercise, we show that the method presented here is generally at least as accurate as other commonly used methods—it is often better. It also provides better uncertainty estimates for the imputations. Lastly, as is shown in Bojinov, Pillai, and Rubin (2017), the copula method can also be used to test some of the underlying assumptions about the appropriateness of imputations for a given data set.

Common Approaches to MI

The standard techniques employed to deal with missing data require an assumption regarding the missing data pattern; these were first formalized in Rubin (1976).⁴ To briefly summarize these terms, missing data are missing completely at random (MCAR) when the probability of the observed missing data pattern is unchanged regardless of what values both the observed and missing data take (Marini, Olsen, and Rubin 1980). The missing data are missing at random (MAR) when the probability of observing the missing data pattern is unchanged no matter what values the missing data take. Finally, the missing data are missing not at random (MNAR) when the probability of observing the missing data pattern changes for some values of the missing data.

These definitions are important both from a theoretical and a practical point of view. The most basic methods, such as listwise deletion, generally lead to biased regression coefficients if the missingness process is not MCAR (Graham 2009). To achieve valid inference under the Bayesian and likelihood paradigms, while ignoring the missing data mechanism, we require the weaker MAR assumption.⁵

The most common appropriate approach to dealing with missing data is MI, which refers to any method that replaces the set of missing values with various plausible values, thus obtaining $m$ completed data sets (Rubin 1996). Rubin (1987) initially suggested creating five imputations, but more recently authors recommended using closer to 20 imputations (Van Buuren 2012).⁶ The completed data sets are then separately analyzed using the standard full data techniques, and the resulting quantities of interest from each data set are combined to obtain an overall, average estimate as well as its associated variance. Before moving to introduce the copula method below, we briefly outline two important methods for generating MIs here.

MI with expectation maximization [EM]: This approach uses iterative EM to create complete data sets based on assuming a particular joint distribution. A widely used method for imputation in the social sciences is the Amelia II $R$ package by Honaker and King (2010). In Amelia II, the joint distribution of the data is modeled as a multivariate normal distribution. Amelia II provides an implementation of the EM approach by the use of bootstrapping to derive solutions quickly. One of the disadvantages of imputation via EM is that for large data sets with significant amounts of missing data, it is computationally intensive. This is a trait of EM algorithms in general, as the rate of convergence is proportional to the amount of missing information in the model. Moreover, it is often unclear to what degree modeling the joint distribution of the data as a multivariate normal distribution is appropriate, especially since the data may include binomial and ordinal variables.

Conditional approaches to MI: An alternative method is to model each variable’s imputation via its conditional distribution based on all other variables in the data. One such approach is developed in MICE (Van Buuren 2012), another was developed as the MI package in $R$ (Goodrich et al. 2012). Imputations for fully conditional specification (FCS) methods, such as MICE or MI, are created based on an “appropriate generalized linear model for each variable’s conditional distribution” (Kropko et al. 2014:501). This is done for all variables and iterated until the model converges.

One of the main drawbacks of the FCS is that only under certain conditions, do the individual conditional models define a valid joint distribution. This often leads to pathologies in the convergence of the algorithms (Li, Yu, and Rubin 2012; Chen and Ip 2015). For example, if Y | X is specified to be an exponential random variable with rate $X$ and X | Y is specified to be an exponential random variable with rate $Y$ , it is well known that no joint distribution exists and sequentially sampling from these two distributions generates draws that tend to infinity (Casella and George 1992). More strikingly, example 1 of Li et al. (2012) demonstrates that even when all the conditionals are normal, the order in which the variables are updated in MICE can determine whether the chain will converge to a stationary distribution.

One of the advantages of conditional model specification is that it allows each variable to be modeled based on its specific distribution, which is specified by the researcher. However, this also means the imputation model for each variable in the data has to be correctly specified, which can be “labor-intensive and challenging with even a moderate number of variables” (Murray 2013:41). Moreover, coefficients estimates in the conditional models can suffer significantly when the number of missing observations is large, especially for categorical variables (Murray 2013).

A Copula Approach to Missing Data Imputation

One of the key issues with conditional approaches to imputation, such as MICE, is that they do not necessarily specify a valid joint distribution (such as the example in the previous section).⁷ When a valid joint distribution does not exist, then there are no guarantees that the MI procedure is proper (as defined in Rubin (2004)). A natural approach to overcoming a possibly incompatible conditional specification is by specifying the joint distribution directly. For example, this is done in most EM approaches, such as Amelia II, by simply assuming a multivariate normal distribution. However, while an approximation, most social science data include binary and ordinal variables and thus cannot have a multivariate normal joint distribution. As a result, this misspecification of the joint distribution is problematic. Moreover, specifying the correct joint distribution becomes increasingly complicated as the number of covariates in the model increase.

It is therefore valuable to decouple the specification of the marginal distribution of each covariate from the function that describes the joint behavior of all covariates together. One of the main advantages of using copulas for imputation is that they allow us to do exactly that. Sklar’s (1959) theorem guarantees that every joint distribution can be decomposed in this way:

Theorem 3.1: (Sklar’s theorem) Let $F$ be a $p$ -dimensional joint distribution function with marginals $F_{1}, \dots, F_{p}$ . Then, there exists a copula $C$ with uniform marginals such that

F (x_{1}, \dots, x_{p}) = C (F_{1} (x_{1}), \dots, F_{p} (x_{p})) .

Sklar’s theorem guarantees that the function $C$ is unique if the marginal distributions $F_{1}, \dots, F_{p}$ are continuous. If they are discrete, then it is unique on the cross product of the ranges of the F_j.

Much work has been done studying the class of Gaussian copulas where the multivariate dependence is defined by $C$ via the multivariate normal distribution with a correlation matrix $R$ (Chen, Fan, and Tsyrennikov 2006; Hoff 2007; Klaassen et al. 1997; Pitt, Chan, and Kohn 2006). That is, we define the Gaussian copula function as $C (\cdot | R) = Φ_{p} (Φ^{- 1} (u_{1}), \dots, Φ^{- 1} (u_{p}) | R)$ for $u_{1}, \dots, u_{p} \in {(0, 1)}^{p}$ where $Φ$ is the univariate normal cumulative distribution function (CDF) and $Φ_{p} (\cdot | R)$ is the $p$ -dimensional CDF with correlation matrix $R$ . This means that the joint distribution of the $p$ variables is given by $F (x_{1}, \dots, x_{p}) = Φ_{p} (Φ^{- 1} (F_{1} (x_{1})), \dots, Φ^{- 1} (F_{p} (x_{p})) | R)$ . Simply put, the univariate CDFs $F_{1}, \dots, F_{p}$ of the individual variables are bound together as a multivariate normal CDF where $R$ determines the correlation between the individual variables on the normal scale.

As previously noted, the specification of marginal distributions is difficult in applied settings and so of particular interest is the setting where the researcher does not need to specify the marginal distributions for $F_{1}, \dots, F_{p}$ . In fact, one big advantage of the method discussed here is that we consider a semiparametric approach that does not require parameterizing the $p$ marginal distributions.

In this flexible setting, the estimation procedures described below provide consistent and likely asymptotically efficient estimates of the dependence parameters in the Gaussian copula, that is, $R$ above (Murray et al. 2013; Hoff, Niu, and Wellner 2014). These dependence parameters directly impact the imputation of the missing data, and thus these theoretical results are extremely appealing. The estimation approach we explore below was developed by Hoff (2007) by extending the ideas of the rank likelihood of Pettitt (1982) to the copula setting.

The rank likelihood (Pettitt 1982) is a type of marginal likelihood that bases inference on the ranks of data rather than the full data. In a univariate setting, it is defined as follows: consider $z_{1}, \dots, z_{n} | θ \sim p (z | θ)$ be a sample from some distribution. Instead of observing the actual values $z_{1}, \dots, z_{n}$ , however, consider only observing the ordering of the data $x_{1}, \dots, x_{n}$ (i.e., their rank). Then, the rank likelihood is given by:

L (θ; x_{1}, \dots, x_{n}) = \int_{D} p (z_{1}, \dots, z_{n} | θ) d z_{1}, \dots, d z_{n}

where $D = {z_{α_{1}} < \dots < z_{α_{n}}}$ and $α_{i} = j$ if and only if $z_{j}$ is the $i$ th smallest of $z_{1}, \dots, z_{n}$ .

Hoff (2007) extends the rank likelhood to the multivariate setting by considering the semiparametric Gaussian copula. Let $z_{1}, \dots, z_{n} | R \sim N (0, R)$ , with $z_{i} = (z_{i 1}, \dots, z_{i p})$ , and let $x_{i j} = F_{j}^{- 1} (Φ (z_{i j}))$ . That is, latent data are drawn from a multivariate normal distribution with correlation $R$ and are transformed to the observed scale via an inverse transformation as in the definition of the Gaussian copula above. One can consider the observed data as the ranks of the unobserved latent $Z s$ and define

D = {Z \in R^{n \times p} : max {z_{k j} : x_{k j} < x_{i j}} < z_{i j} < min {z_{k j} : x_{i j} < x_{k j}}} .

It is easy to see that all $Z \in D$ respect the order of the variables on the observed scale. Hoff (2007) shows that $P (Z \in D | R, F_{1}, \dots, F_{p}) = P (Z \in D | R)$ which in turn allows for the decomposition

P (X | R, F_{1}, \dots, F_{p}) = P (Z \in D | R) P (X | Z \in D, F_{1}, \dots, F_{p}) .

The aforementioned results guarantee that inference about $R$ can proceed simply via $P (Z \in D | R)$ . This leverages the ordering of the observed values $x_{1 j}, \dots, x_{n j}$ of each variable to make inference about the parameter $R$ without estimating the CDFs $F_{1}, \dots, F_{p}$ .

This means that regardless of the marginal distributions of the individual variables, all we need is their ordering to facilitate the use of the Gaussian copula model to make inferences about the dependence between these variables, that is, the correlation matrix $R$ . A Bayesian approach to estimating $R$ specifies an inverse Wishart prior for a covariance matrix $V$ such that $R$ is its correlation matrix and a normal prior for the latent $z_{i j}$ . Updates are performed via a Gibbs sampler since full conditional distributions can be derived by conditioning on the ranks of the data alone.⁸

Let us summarize and paraphrase the method in less technical terms. Assume we have two vectors $Z_{1}$ and $Z_{2}$ which come from a bivariate normal distribution with correlation $R$ . We observe $X_{i} = F_{i}^{- 1} (Φ (Z_{i}))$ implying that $X_{i}$ is distributed according to $F_{i}$ . If the $F_{i}$ are continuous and known, we can recreate the vectors $Z_{1}$ and $Z_{2}$ by using the pseudoinverse CDFs on the original data ( $Z_{1} = Φ^{- 1} (F_{1} (X_{1}))$ and $Z_{2} = Φ^{- 1} (F_{2} (X_{2}))$ ). We could then generate a good estimate of $R$ using the transformed vectors $Z_{1}$ and $Z_{2}$ and maximum likelihood estimation, for example $Σ_{i = 1}^{N} \frac{Z_{1}^{(i)} Z_{2}^{(i)}}{N}$ would be a natural estimate for the correlation. However, when either vector is not continuous, the simple pseudoinverse transformation does not allow for correct estimation of the correlation. Now assume $X_{2}$ is a binary or ordinal variable as in many of our cases but that the marginal is not known. Instead of using a plug in value for $Z_{2}$ (say, by estimating the marginal $F_{2}$ ), we contend that the ranks of our latent continuous $Z_{2}$ are the same as those of the observed variable $X_{2}$ . The estimation procedure then iterates the following two steps: Using the ranks of $X_{1}$ and $X_{2}$ and the current estimate of the multivariate correlation $R$ , we can draw values of the latent variables ( $Z_{1}$ and $Z_{2}$ ) that preserve the rank ordering of the observed data. The second step uses the sampled underlying latent variables to sample the correlation $R$ . These steps are iterated until stationarity is reached. Relying on the ranks and latent scale allows us to not specify the marginal distributions of the individual variables and still arrive at a proper solution to estimating $R$ .

When values of $x_{i j}$ are MAR, imputation can be performed first on the latent z_ij scale (since the latent variables are normal, sampling from the conditional distribution of the missing data given the observed data requires a multivariate normal draw) and is then transformed to the observed scale using the empirical cumulative density functions. As this is a Bayesian procedure, we produce a posterior for the missing data. To make our approach comparable to the standard conditional approaches, we only employ a few samples from this posterior and use those as multiply imputed data sets. However, it is natural to consider posterior predictive distributions of parameters of interest or other posterior summaries on a case-by-case basis. For example, the conditional independence graphs of Hoff (2007) succinctly summarize the relationships among many variables.

Comparing Amelia II, sbgcop, and MICE

In this section, we compare the working properties of the copula-based imputation with those of Amelia II and MICE packages. We evaluate each method based on an extensive simulation study as well as an empirical example from the social sciences, discussed in the next section.

Evaluating Imputations

MI procedures are specifically designed to yield valid statistical inference (meaning, asymptotically unbiased with correct standard errors and coverage) for population quantities of interest. Since correct estimation of the coefficients and standard errors is critical for obtaining valid statistical inference, any analysis of MI procedures must focus on studying its frequentist properties. Properties such as empirical coverage, average bias, and average interval length of the estimate of the scientific estimand over repeat samples will be of cardinal interest.

We therefore use the following approach to assess the validity of a MI procedure through simulation:

Define a full data quantity of interest, $θ$ . In our setting, $θ$ is a set of regression coefficients.

Generate a complete data set and apply a prespecified missing data mechanism to remove some observations.

Use the MI procedure to create $m$ completed data sets with the missing values replaced by imputed values.

Use each of the $m$ data sets to obtain an estimate of $θ$ as well as its associated variance and combine them using Rubin’s combining rules (Rubin 2004) to obtain $θ$ and a 95 percent confidence interval (CI).

Report the bias of $θ$ , the CI interval length and whether or not the CI covered the true value (Van Buuren 2012, section 2.5.2).

We repeat steps 2–5 $S$ times to obtain the empirical coverage rate. By varying the full data model and the missing data mechanism, in step 2, we can control the two paths that influence the effectiveness of the MI procedures.

Simulation Study

In regression settings, an outcome $Y$ can depend on many explanatory variables $X = X_{1}, \dots, X_{J}$ , some of which can be costly to measure. As such, it is common that while the outcome Y is measured for all variables, some entries of the design matrix $X$ are missing. In this simulation, we exclusively focus on this situation and restrict the missingness to the explanatory variables. We will further assume that the missingness mechanism does not allow for the missingness to depend on the outcome $Y$ .

In this situation, complete case analysis (or listwise deletion) provides an unbiased estimate of the regression coefficients; however, the reduced sample size often leads to losses in efficiency, through higher standard errors. Another disadvantage of using complete case analysis whenever the number of explanatory variables $J$ is of moderate size is that the probability of having enough complete cases to estimate the regression coefficients is low. In this setting, using an MI procedure is paramount and leads to a significant reduction in the standard errors; however, this can induce a slight bias. White and Carlin (2010) show through an extensive simulation study that the increase in bias often time leads to a decreased empirical coverage rate for both MAR and MNAR data sets.

For our simulation study, we set $J = 40$ , $N = 1, 000$ , and consider $X_{j}$ that include both continuous and discrete variables to demonstrate the versatility of the copula approach without specifying any of the marginal distributions. This is precisely the scenario we described above; the probability of enough complete cases existing to estimate the regression coefficients is effectively 0.⁹

The distributions we consider for the elements of the design matrix are Gaussian, Bernoulli, Poisson, and ordinal. To make imputation feasible, we require the variables to be correlated. To generate correlated variables, we first construct a matrix of correlated Gaussian random variables and then transform the variables to have the appropriate marginals. For example, to generate a pair of correlated Poisson random variables $A$ and $B$ with mean $λ$ , we construct $(Z_{1}, Z_{2}) \sim N (0, Σ)$ where $σ_{11} = σ_{22} = 1$ and $σ_{12} = σ_{21} = ρ$ and set $A = F_{Pois, λ}^{- 1} (F_{N} (Z_{1}))$ and $B = F_{Pois, λ}^{- 1} (F_{N} (Z_{2}))$ . The data generating process thus leads to the following marginal distributions for the entries in $X$ : for $j = 1, \dots,10$ ,

\begin{array}{l} X_{j} \sim N (0, σ_{j}^{2}) & X_{j + 10} \sim Bern (p_{j}) \\ X_{j + 20} \sim Pois (λ_{j}) & X_{j + 30} \sim ordinal (0, 1) \\ X = (X_{1}, \dots, X_{40}) & Y \sim N (\sum_{i = 1}^{40} X_{i},1), \end{array}

where $σ_{j} = 1 + (j - 1) / 9$ , $λ_{j} = 0.2 + 2 (j - 1) / 90$ , and $p_{j} = 2 + 3 (j - 1) / 9$ . Both the amount of missingness (MC) and correlation ( $p$ ) between the different variables are varied according to the specified values given in Table 1.

Table 1.

Simulation Study Configurations.

Correlation (p)	Missingness Coefficient
.2	.3
.35	.4
.5	.5
.65	.6

We consider two missing data mechanisms for $X$ , one that produces MAR data sets and another one that generates MNAR data sets (see Online Appendices A and B for details). The MI procedures we considered are only valid under the MAR assumption; however, it is useful to check how each method performs when this assumption is violated—as is often the case in practice.

Results

We performed 1,000 simulations under each of the possible combinations of the correlation and missingness coefficient, as detailed in Table 1, under both MAR and MNAR missing data mechanisms. For MICE, we specified the correct marginal distributions (e.g., ordered logit model for the ordinal variables). For Amelia II, we used the appropriate variable transformation in accordance with the package help files. For the copula approach, we did not need to specify any distributions/transformations. Using each of the three procedures, we created 20 completed data sets that were used to estimate the regression coefficients and a corresponding 95 percent CI.¹⁰ None of the simulations had enough complete cases to estimate the regression coefficients using listwise deletion.

The most significant source of variation in the simulation was due to the different classes of variables, followed by the correlation and the missingness coefficient. There is only a small difference in the results obtained from the MAR and MNAR data; therefore, our discussion will focus on the former, with the figures for the latter included in Online Appendix C. Figures 2 and 3 illustrate how the bias, coverage, and interval length vary across the interaction of the different variable classes, the correlation, and the missingness coefficient, respectively. Overall, the copula method achieved an empirical coverage rate of 93 percent which was much higher than that of MICE, 87 percent, and Amelia II, 83 percent. Less adversarial regimes were previously studied in White and Carlin (2010), by reducing the number of covariates in our simulation, we can recover similar coverage rates for the MI procedures as are reported there. Both the copula and MICE methods had an absolute average bias of .17. Amelia II performed worse and had a bias of .25. On average, all three methods had approximately the same interval length.

Figure 2.

Simulation study results for the missing at random data as a function of the missingness coefficient, averaging over the correlation. The plot is split by the different variable types (normal, binomial, Poisson, and ordinal) and the three outcomes of interest (bias, coverage, and interval length). The rightmost panel shows the result averaging over the different variable types.

Figure 3.

Simulation study results for the missing at random data as a function of the correlation, averaging over the missingness coefficient. The plot is split by the different variable types (normal, binomial, Poisson, and ordinal) and the three outcomes of interest (bias, coverage, and interval length). The rightmost panel shows the result averaging over the different variable types.

The copula imputations were obtained using 10,000 iterations from Hoff’s (2010) package whose convergence was checked on a subset of simulations. The lag-10 autocorrelation for the thinned chains was less than .18 in absolute value for each of the elements of the latent correlation matrix, and the effective sample size was always above 200 (97.6 percent of the entries were above 500). Since the copula method is sampling from the posterior distribution which requires the MCMC algorithm to converge to the stationary distribution, its computation time depends on the rate of convergence as well as the desired number of imputations. Running multiple MCMC chains in parallel to generate independent imputations can reduce the computation time. This approach is slightly slower than Amelia II but is substantially faster than the standard MICE algorithm where all $J - 1$ variables are used to impute the jth variable. Fortunately, the copula algorithm scales well as the sample size and the number of explanatory variables increases. The copula method had the lowest bias, highest coverage rate, and often the longest interval length. It is noteworthy that even though the semiparametric estimation procedure did not require specification of the marginals, any data transformations, or tuning, it still outperformed the other two procedures.

Since the MICE procedure is iterative, users need to check that the model parameters fully explore the parameter space. Unlike the Bayesian copula method, there are no explicit convergence criteria that can be tracked. We performed a visual check that revealed no abnormalities and also ran each MICE chain for 20 iterations as recommended in Van Buuren and Groothuis-Oudshoorn (2011). The MICE method performed almost as well as the copula method but had a slightly lower coverage rate, meaning the estimated standard errors were too small. MICE also had the smallest average bias for the normal and Poisson variables. Again, however, these results are contingent on specifying the correct conditional distribution which can often be challenging.

Amelia II had the lowest coverage and highest bias both on average and in most scenarios that we considered. It had the smallest average interval length of 1.23, which shows that it was systematically underestimating the variance: leading to the low coverage rates.

Figure 2 shows that the average bias and the interval length increase as a function of the proportion of missing values. This leads to a decrease in the empirical coverage as the bias increases at a faster rate than the interval length. One notable exception was the correct coverage of the copula approach for the regression parameters of the ordinal and binomial variables, both Amelia II and MICE undercovered. Given that these types of variables are frequently encountered in social science applications, these results especially suggest that using a copula approach can lead to better statistical conclusions. Moreover, the overall simulation results indicate that when a normal distribution does not well approximate the data, then the copula approach will consistently outperform both Amelia II and MICE.

Somewhat surprisingly, there seems to be less variation in the bias and the interval length as a function of the correlation, as is shown in Figure 3. Except for the normally distributed variables, the bias decreases as the correlation increases due to the reduction in the relative loss of information from the missing data.

Breaking the MAR assumption did not lead to drastically worse results. We observe a decrease of about 3 percent in the coverage of all three methods and a slight decrease in the average bias. This shows that the methods are somewhat robust to violations of MAR assumption when it is not too severe. Figures C1 and C2 in the Online Appendix C show the results of the simulations when the MAR assumption is violated.

Application Study

In this section, we provide a comparison of the three imputation methods using an application from political science. The empirical example shows how copula methods can be used to generate imputations in a large data set with a variety of variable types.

Inequality and Democratic Support

As we have elaborated above, imputation methods are still underused, especially in the social sciences. There is, however, some visible progress. One example where scholars have taken advantage of one of the imputation methods currently available is “Economic Inequality and Democratic Support” by Krieckhaus et al. (2014) published in the Journal of Politics. Krieckhaus et al. (2014) explore whether the support for democracy within countries is affected by the level of inequality. The authors combine country-level variables (such as inequality) with individual-level survey data from 40 democracies around the world. For multiple countries, several survey waves are included, resulting in 57 country-years and a total of 77,642 observations (Krieckhaus et al. 2014:144). For this replication exercise, we replicate model 1 in Table 1 in Krieckhaus et al. (2014). The dependent variable is a “13-point additive index (ranging from 0 to 12) of democratic support,” which the authors treat as a continuous variable (Krieckhaus et al. 2014:144). The main independent variables of interest are Inequality at the country level and an ordinal Income scale at the individual level (ranging from 1 to 10). Additionally, the authors control for Age, Gender, Institutional Confidence, Interest in Politics, Interpersonal Trust, Education, Prior Regime Evaluation, and Leftist Ideology all drawn from the World Values Survey (2012). As in the original article, all individual-level variables are demeaned “using group-mean centering” after the imputation (Krieckhaus et al. 2014:145). The data are analyzed using a random coefficients model.

Most importantly for this study, the original data suffer from a relatively high number of missing observations. Table 2 shows the share of missing observations for variables included in the replication exercise. We can see that many of the variables have a large share of missing observations. If instead of MIs, the authors used listwise deletion then the number of observations in the regression model would have been approximately halved. Instead, Krieckhaus et al. (2014) use Amelia II to multiple impute five data sets which they analyze. Estimates are then combined using Rubin’s rule.

Table 2.

Share of Missingness in Variables of Interest.

Democracy Support	Inequality	Income	Age
19.9	1.8	12.9	0.2
Gender	Institutional Confidence	Interest in Politics	Interpersonal Trust
0.1	11.7	2.5	3.7
Education	Leftist Ideology	Prior Regime Evaluation
3.9	18.5	21.3

This is an excellent setting for our comparison of MI techniques. The number of missing observations is quite large, and the data set includes different types of variables, continuous, binary, and ordinal. We created 20 multiple imputed data sets using each of the imputation techniques: Amelia II, MICE, and sbgcop. We then reestimate model 1 in table 1 in Krieckhaus et al. (2014:147) and combine the estimation results for each method’s multiple imputed data sets via Rubin’s rule.

For Amelia II, we specify the type of each variable and then generate 20 imputed data sets using the full original data. Similarly, we declare each variable’s type for MICE and estimate the default model for each. We use all variables except the one to be imputed as independent variables in the chained equations. Again, we create 20 multiple imputed data sets and set the maximum number of iterations to 20.

Lastly, we use our preferred method, imputation via the semiparametric Gaussian copula, to generate 20 imputed data sets. We run the MCMC chain for 2,100 iterations and randomly draw 20 data sets from the posterior. Note that, again, we do not have to declare any of the variable types or make any other specification or transformation of the data.

Figure 4 shows the coefficient estimates and 95 percent intervals for the replicated model based on each of the imputation techniques, as well as when listwise deletion is used. First, the results are quite similar for the Inequality, Income, and Age variables. Even for the models based on listwise deletion. For the two main variables of interest, inequality, and income, the results based on different imputation techniques are virtually the same.

Figure 4.

Coefficient estimates and confidence intervals for model 1 in table 1 in Kriechkhaus et al. (2014) based on three imputation techniques and listwise deletion.

On the other hand, there are several significant differences for the other variables included in the model. First, the effect of gender is essentially zero according to the models estimated on the copula imputed data. Based on the data imputed using MICE or Amelia II, females have higher ratings of democracy satisfaction (though the CIs just cover zero). According to the nonimputed data, the effect of gender is quite strong.

Based on the data imputed with the copula method, the estimated association of Institutional Confidence with Democracy Satisfaction is significantly stronger compared to the models based on listwise deletion or other imputation methods. Similarly, the estimated effect of Leftist Ideology is also substantially larger according to the copula imputed data. On the other hand, the association of Education levels with Democracy Satisfaction is significantly smaller. Based on the copula, the relationships of Interest in Politics and Prior Regime Evaluation with the dependent variable of Democracy Satisfaction are all modeled to be weaker, compared to the other methods (and the nonimputed data), though the CIs overlap.

It is interesting to note that, except for one variable (Interpersonal Trust), whenever the estimated coefficient for the copula imputed data differs from the coefficients based on the other imputation methods, it is in the opposite direction of the difference to the listwise deletion coefficient. This is especially easy to see for the Gender and Leftist Ideology variables, where the effect is strongest (weakest) according to the model estimated on the listwise deleted data and weakest (strongest) for the copula-based models.

Based on the simulation results, especially with respect to binary and ordinal variables, and the theoretical properties, we are confident in the accuracy of the copula imputation method. These results suggest that Gender is not associated with people’s satisfaction with democracy, whereas Institutional Confidence and Left ideology both have much stronger effects.

Conclusion

What practical lessons can we learn about how to deal with missing data? In this article, we reemphasize the importance of dealing with missing data and present a copula based approach, developed by Hoff (2007), that is elegant and requires little prespecification of the data. With the rank-based approach introduced by Hoff (2007), the Gaussian copula can be used to impute binary, ordinal, and continuous variables. We discuss the theoretical properties of the copula method and its theoretical attractiveness compared to other commonly employed techniques. In particular, the Gaussian copula introduced here enables researchers to make imputation via draws from a valid posterior of the joint distribution without specifying the distributions of the individual variables. Moreover, we present evidence from simulations that it performs better than either Amelia II or MICE, especially when it comes to nonnormally distributed data.

While the three imputation methods perform relatively similarly, throughout the simulation, the Copula method does have the lowest average bias (tied with MICE) and the highest coverage rate (93 percent). More so, MICE requires specification of the conditional distributions whereas the copula method does not. Recent theoretical results for MICE suggest that good performance heavily relies on being approximately correct in the choice of conditionals (Li et al. 2012). On the other hand, theoretical guarantees for good behavior of copula methods are available. In particular, information bounds for rank-based estimators are the same as the information bounds for estimators based on the full (scale and rank) data (Hoff et al. 2014). Under MAR and MCAR, we inherit all the properties of the full data, and by introducing structure to the imputation, we are likely to have good behavior even under MNAR.

One aspect that we have not addressed herein is the validity and sensitivity to the unassessable assumptions made when analyzing data with missing values (Molenberghs et al. 2014), that is, the type of missingness mechanism. Bojinov et al. (2017) show that the Gaussian copula approach can be used to assess the validity of the missing always at random assumption (a slightly stronger assumption that implies MAR). Their results suggest that by using a Gaussian copula for generating imputations, the analysts can also easily diagnose the assumptions they made and quickly identify variables which are likely to break these assumptions. This adds another benefit to using the method discussed in this article.

Consideration must also be given to the computational cost of any procedure. As indicated by Graham (2009), the disadvantages of EM approaches are especially large when imputing databases with many variables or applications of “big data.” MICE can be computationally less expensive but suffers when the number of variables increases as the correct choice for each of the conditionals becomes increasingly unlikely. The semiparametric copula approach described here relies on MCMC, its speed does not depend on the fraction of missing data and scales nicely in the dimension of the data set. This makes it possible to impute even large database in a relatively timely manner and no prespecification of the data. Moreover, using the copula model to multiply impute missing values provides some of the advantages (such as a proper posterior distribution of the data) but is less burdensome on scholars than imputing values in a fully Bayesian approach (Erler et al. 2016).

Finally, the copula approach is quite flexible and can be employed at different stages of the analysis process. First, it can be used to generate a single estimate of the missing data or the mean of a large number of draws, which is exactly what might be needed in some situations. Second, per the recommendation of Rubin, it can be used to construct multiple databases. As with Amelia II, the copula imputations can be analyzed separately and the results combined using either mitools or Zelig (Imai, King, and Lau 2008) in $R$ . Thus, the copula approach to missing data can be explicitly integrated into the modeling and analysis of observational data in a simplistic, organic fashion.

Supplemental Material

Supplemental Material, Suppl_SMR - Multiple Imputation Using Gaussian Copulas

Supplemental Material, Suppl_SMR for Multiple Imputation Using Gaussian Copulas by Florian M. Hollenbach, Iavor Bojinov, Shahryar Minhas, Nils W. Metternich, Michael D. Ward and Alexander Volfovsky in Sociological Methods & Research

Footnotes

Acknowledgments

For helpful insights, we thank Philippe Loustaunau, among the first of our colleagues to encourage this effort. Stephen Shellman was a strong critic who deserves our thanks too: His criticisms helped us to improve our approach. John Ahlquist, Matt Blackwell, Andreas Beger, Cassy Dorff, Gary King, and Jacob Montgomery provided helpful comments on previous versions of this article.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This project was partially supported by the Office of Naval Research (holding grants to the Lockheed Martin Corporation, Contract N00014-12-C-0066). Nils W. Metternich acknowledges support from the Economic and Social Research Council (ES/L011506/1). The work was completed while Alexander Volfovsky was supported by an NSF MSPRF under DMS-1402235.

ORCID iDs

Florian M. Hollenbach

Nils W. Metternich

Supplemental Material

Supplemental material for this article is available online.

Notes

References

Bojinov

Pillai

Rubin

. 2017. “Diagnosing Missing Always at Random in Multivariate Data.” arXiv preprint arXiv:1710.06891.

Casella

George

E. I.

. 1992. “Explaining the Gibbs Sampler.” The American Statistician 46:167–74.

Chen

S.-H.

E. H.

. 2015. “Behaviour of the Gibbs Sampler When Conditional Distributions Are Potentially Incompatible.” Journal of Statistical Computation and Simulation 85:3266–75.

Chen

Fan

Tsyrennikov

. 2006. “Efficient Estimation of Semiparametric Multivariate Copula Models.” Journal of the American Statistical Association 101:1228–40.

Di Lascio

Giannerini

Reale

. 2015. “Exploring Copulas for the Imputation of Complex Dependent Data.” Statistical Methods & Application 24:159–74.

Erler

N. S.

Rizopoulos

Rosmalen

J. V.

Jaddoe

V. W.

Franco

O. H.

Lesaffre

E. M.

. 2016. “Dealing with Missing Covariates in Epidemiologic Studies: A Comparison between Multiple Imputation and a Full Bayesian Approach.” Statistics in Medicine 35:2955–74.

Goodrich

Kropko

Gelman

Hill

. 2012. “mi: Iterative Multiple Imputation from Conditional Distributions.” R package. Retrieved July 14, 2016. (https://cran.r-project.org/web/packages/mi/).

Graham

J. W.

2009. “Missing Data Analysis: Making It Work in the Real World.” Annual Review of Psychology 60:549–76.

Hoff

2010. “sbgcop: Semiparametric Bayesian Gaussian Copula Estimation and Imputation.” R Package Version 0.975. Retrieved July 14, 2016. (https://CRAN.R-project.org/package=sbgcop).

10.

Hoff

P. D.

2007. “Extending the Rank Likelihood for Semiparametric Copula Estimation.” Annals of Applied Statistics 1:265–83.

11.

Hoff

P. D.

Niu

Wellner

J. A.

. 2014. “Information Bounds for Gaussian Copulas.” Bernoulli: Official Journal of the Bernoulli Society for Mathematical Statistics and Probability 20:604.

12.

Honaker

King

. 2010, April. “What to Do About Missing Values in Time-series Cross-section Data.” American Journal of Political Science 54:561–81.

13.

Honaker

King

Blackwell

. 2012. “AMELIA II: A Program for Missing Data—Documentation.” Retrieved July 14, 2016. (https://r.iq.harvard.edu/docs/amelia/amelia.pdf).

14.

Imai

King

Lau

. 2008. “Toward a Common Framework for Statistical Analysis and Development.” Journal of Computational and Graphical Statistics 17:892–913.

15.

Käärik

2006. “Imputation Algorithm Using Copulas.” Metodoloski zvezki 3:109.

16.

Käärik

. 2009. “Modeling Dropouts by Conditional Distribution, a Copula-based Approach.” Journal of Statistical Planning and Inference 139:3830–35.

17.

King

Honaker

Joseph

Scheve

. 2001. “Analyzing Incomplete Political Science Data: An Alternative Algorithm for Multiple Imputation.” American Political Science Review 95:49–69.

18.

Klaassen

C. A.

Wellner

J. A.

. 1997. “Efficient Estimation in the Bivariate Normal Copula Model: Normal Margins Are Least Favourable.” Bernoulli 3:55–77.

19.

Krieckhaus

Son

Bellinger

Wells

. 2014. “Economic Inequality and Democratic Support.” The Journal of Politics 76:139–51.

20.

Kropko

Goodrich

Gelman

Hill

. 2014. “Multiple Imputation for Continuous and Categorical Data: Comparing Joint Multivariate Normal and Conditional Approaches.” Political Analysis 22:497–519.

21.

Rubin

D. B.

. 2012. “Imputing Missing Data by Fully Conditional Models: Some Cautionary Examples and Guidelines.” Department of Statistical Science, Duke University, Discussion Paper 1124. Retrieved July 2, 2016. (http://ftp.stat.duke.edu/WorkingPapers/11-24.pdf).

22.

Little

R. J.

Rubin

D. B.

. 2002. Statistical Analysis with Missing Data. 2nd ed. New York: Wiley.

23.

Liu

Gelman

Hill

Y.-S.

Kropko

. 2013. “On the Stationary Distribution of Iterative Imputations.” Biometrika 101:155–73.

24.

Marini

M. M.

Olsen

A. R.

Rubin

D. B.

. 1980. “Maximum-likelihood Estimation in Panel Studies with Missing Data.” Sociological Methodology 11:314–57.

25.

Mealli

Rubin

D. B.

. 2015. “Clarifying Missing at Random and Related Definitions, and Implications When Coupled with Exchangeability.” Biometrika 102:995–1000.

26.

Molenberghs

Fitzmaurice

Kenward

M. G.

Tsiatis

Verbeke

. 2014. Handbook of Missing Data Methodology. Boca Raton, FL: Chapman and Hall/CRC.

27.

Murray

J. S

. 2013. “Some Recent Advances in Non- and Semiparametric Bayesian Modeling with Copulas, Mixtures, and Latent Variables.” Dissertation, Department of Statistical Science, Duke University. Retrieved July 2, 2016 (http://dukespace.lib.duke.edu/dspace/handle/10161/8253).

28.

Murray

J. S.

Dunson

D. B.

Carin

Lucas

J. E.

. 2013. “Bayesian Gaussian Copula Factor Models for Mixed Data.” Journal of the American Statistical Association 108:656–65.

29.

Pettitt

1982. “Inference for the Linear Model Using a Likelihood Based on Ranks.” Journal of the Royal Statistical Society. Series B (Methodological) 44:234–43.

30.

Pitt

Chan

Kohn

. 2006. “Efficient Bayesian Inference for Gaussian Copula Regression Models.” Biometrika 93:537–54.

31.

R Development Core Team. 2004. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria (http://www.R-project.org/).

32.

Robbins

M. W.

Ghosh

S. K.

Habiger

J. D.

. 2013. “Imputation in High-dimensional Economic Data as Applied to the Agricultural Resource Management Survey.” Journal of the American Statistical Association 108:81–95.

33.

Rubin

D. B.

1976. “Inference and Missing Data.” Biometrika 63:581–92.

34.

Rubin

D. B

. 1987. Multiple Imputation for Nonresponse in Surveys. Hoboken, NJ: John Wiley.

35.

Rubin

D. B.

1996. “Multiple Imputation after 18+ Years.” Journal of the American Statistical Association 91:473–89.

36.

Rubin

D. B

. 2004. Multiple Imputation for Nonresponse in Surveys. Vol. 81. Hoboken, NJ: John Wiley.

37.

Sklar

1959. “Fonctions de Répartition á N Dimensions et Leur Marges.” Publications de l’Institut Statistique de l’Université Paris 8:229–31.

38.

Van Buuren

2012. Flexible Imputation of Missing Data. Boca Raton, FL: Chapman & Hall/CRC Press.

39.

Van Buuren

Groothuis-Oudshoorn

. 2011. “MICE: Multivariate Imputation by Chained Equations in R.” Journal of Statistical Software 45:1–67.

40.

White

I. R.

Carlin

J. B.

. 2010. “Bias and Efficiency of Multiple Imputation Compared with Complete-case Analysis for Missing Covariate Values.” Statistics in Medicine 29:2920–31.

41.

World Values Survey. 2012. 1981–2008 Integrated Questionnaire. World Values Survey Association (www.worldvaluessurvey.org, accessed January 1, 2012). Aggregate File Producer: JDSystems Data Archive, Madrid, Spain.

42.

Yucel

R. M.

2011. “State of the Multiple Imputation Software.” Journal of Statistical Software 45:1–7.

43.

Zhu

Raghunathan

T. E.

. 2015. “Convergence Properties of a Sequential Regression Multiple Imputation Algorithm.” Journal of the American Statistical Association 110:1112–24.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.42 MB