Edoardo Costantini Kyle M. Lang Tim Reeskens Klaas Sijtsma Tilburg University,
Department of Methodology and Statistics Utrecht University, Department of
Methodology and Statistics Tilburg University, Department of Sociology Tilburg
University, Department of Methodology and Statistics Corresponding author’s
email address:
<EMAIL_ADDRESS>
This article has been submitted to a journal but it was published on
https://arxiv.org before formal peer review.

# High-dimensional imputation for the social sciences: a comparison of state-
of-the-art methods

###### Abstract

Including a large number of predictors in the imputation model underlying a
multiple imputation (MI) procedure is one of the most challenging tasks
imputers face. A variety of high-dimensional MI techniques can help, but there
has been limited research on their relative performance. In this study, we
investigated a wide range of extant high-dimensional MI techniques that can
handle a large number of predictors in the imputation models and general
missing data patterns. We assessed the relative performance of seven high-
dimensional MI methods with a Monte Carlo simulation study and a resampling
study based on real survey data. The performance of the methods was defined by
the degree to which they facilitate unbiased and confidence-valid estimates of
the parameters of complete-data analysis models. We found that using
regularized regression to select the predictors used in the MI model and using
principal component analysis to reduce the dimensionality of auxiliary data
produce the best results.

##  1 Introduction

Today’s social, behavioral, and medical scientists have access to large
multidimensional data sets that can be used to investigate the complex roles
that social, psychological, and biological factors play in shaping individual
and societal outcomes. Large social scientific data sets—such as the World
Values Survey and the European Values Study (EVS)—are easily accessible to
researchers, but making use of the full potential of these data requires
dealing with the crucial problem of multivariate missing data.

###  1.1 The state of imputation in Sociology

Sociologists working with social surveys are usually interested in drawing
inferential conclusions based on a substantively interesting analysis model.
Generally, these analysis models require complete-data, so the researcher must
address any missing values before moving on to their substantive analysis.
There are many possible missing data treatments from which to choose, and
their relative strengths and weaknesses are covered elsewhere (e.g., Little
Rubin, 2002; Enders, 2010; Van Buuren, 2018). In this paper, we will focus on
Rubin (1987)’s multiple imputation (MI), which is one of the most effective
ways of addressing missing values in survey data.

MI is a three-step procedure that entails imputation, analysis, and pooling
phases. The fundamental idea of the _imputation phase_ is to replace each
missing data point with $d$ plausible values sampled from the posterior
predictive distribution of the missing data, given the observed data. This
phase generates $d$ completed versions of the original data set that are each
analyzed separately during the _analysis phase_ , using any standard complete-
data analysis model. Finally, in the _pooling phase_ , the $d$ sets of
estimates from the analysis models are pooled following Rubin’s rules (Rubin,
1987) to create a single set MI parameter estimates and standard errors.

Missing values are one of the main factors impacting the quality of data
gathered with surveys (Meyer ., 2015), and nonresponse rates in large social
survey have risen drastically over the last two decades (Brick Williams, 2013;
Massey Tourangeau, 2013; Williams Brick, 2018). To explore how sociologists
are addressing the issue of nonresponse in their research, we reviewed how
missing data have been discussed in the articles published over the last five
years in two leading sociological journals: American Journal of Sociology
(AJS) and American Sociological Review (ASR). We found that of the 148 AJS
research articles that mentioned using a survey, or some form of sample, for
inferential analysis, 24 addressed the presence of missing values, and 17
conducted some form of imputation. Of these 17, only 13 performed MI, and
among these 13, only 3 articles gave information on which predictors were used
in the imputation models. Turning to ASR, the picture was similar. Of the 191
research articles published between January 2017 and January 2022 that met the
inclusion criteria described above, 20 reported performing MI. Of these 20
articles, only 6 gave information regarding which predictors were used in the
imputation models. Across the two journals, in the 9 papers we found that
described which predictors were used in the imputation models, the predominant
choice was to use only the analysis model variables in the imputation model.

In general, it seems that even when sociologists pay attention to the problem
of missing values, little attention is given to which variables should be used
in the imputation models. Similar conclusions were drawn in other literature
reviews Mustillo (2012); Mustillo Kwon (2015). However, which variables to
include in the imputation models is a crucial decision in MI. Leaving out
important predictors of missingness can induce _missing not at random_ (MNAR)
data (Collins ., 2001), while including good predictors can both correct for
nonresponse bias and improve the efficiency of the parameters estimates
(Collins ., 2001; von Hippel Lynch, 2013).

###  1.2 The challenge of specifying good imputation models

Specifying the imputation model is one of the most challenging steps in
dealing with missing values. As described by Van Buuren . (1999), the task
involves defining two aspects of the model: the model form (e.g., linear,
logistic) and the predictor matrix (i.e., the set of predictors that enter the
imputation model). The first choice is straightforward in virtually any
imputation task, as it depends primarily on the measurement level of the
variables under imputation. The second choice requires a careful selection
process aimed at identifying the subset of variables that will be most useful
in a given imputation model.

Generally, the variables that will be part of the analysis model should also
be included in the imputation model. When some analyzed variables (including
transformations such as polynomials or interactions) are excluded from the
imputation, the analysis and imputation models are said to be uncongenial
(Meng, 1994). Such uncongeniality can lead to biased parameter estimates and
invalid inferences. When designing an imputation model, the range of analysis
models for which the resulting imputations will be congenial is an important
consideration. In the methodological literature, this concept is known as the
scope of the imputation model. Van Buuren (2018, 46) distinguishes three
typical imputation model scopes:

  * •

_Narrow Scope:_ Narrowly scoped imputation models are matched to individual
analyses. In such a scenario, the imputation is a bespoke pre-processing step
intended to facilitate only a single analysis model. When imputing with a
narrow scope, the primary objective is to ensure that all the variables in the
analysis models (including relevant transofmrations) appear in the imputation
model. An analyst who imputes their own data and plans to estimate only one
model (or a single series of nested models) may wish to specify a narrow
scoped imputation model.

  * •

_Intermediate Scope:_ Imputation models with an intermediate scope are
designed to support several different analysis models. The imputer will
generally know approximately which analyses are intended but may not have an
exhaustive list of all variables that will be analyzed. The objective is to
design an imputation model that will be congenial with all planned and
unplanned analysis models. Such analytic contexts frequently arise within
research teams wherein several different analyses contribute to a larger
research program. The evaluation of the Dating Matters intervention (Tharp,
2012) is an example of one such research program. Due to the size and
complexity of the data and the diversity of the intended analyses, treating
the missing data in the Dating Matters evaluation took several months of
dedicated work (Niolon ., 2019). The resulting imputations were then used to
support the substantive analyses by which various dimensions of the
intervention were evaluated (e.g., Vivolo-Kantor ., 2021; Estefan ., 2021).

  * •

_Broad Scope:_ Imputation models with a broad scope are designed to create
imputations that will be congenial to the most general set of analysis models
feasible. The imputer cannot know beforehand which variables will be part of
the analysis models, so the imputation models are designed to be general
enough to accommodate a wide range of potential analyses. Practically
speaking, the objective is to recreate the moments of the hypothetically fully
observed data as closely as possible. Rubin (1987, 3) originally envisioned MI
as a method using broadly scoped imputation models to treat publicly released
data and argues that well-implemented MI can accomodate models that were not
contemplated by the imputer (Little Rubin, 2002, 218). Any data curation
institution imputing data that are intended for public release will need
imputation models with a broad scope. The Federal Reserve Board’s Survey of
Consumer Finances (Kennickell, 1998) and the Luxembourg Wealth Study (LWS,
2020) are two examples of surveys released after performing MI, and used by
sociologists publishing in AJS and ASR.

Despite its importance, congeniality should not be considered as the sole
guiding principle when defining imputation models. Uncongeniality can
theoretically improve standard complete-data procedures by producing standard
errors that are smaller than those that would have been produced by fully
observed data. This phenomenon is known as superefficiency (Meng, 1994,
544–546; Rubin, 1996, 481; Little Rubin, 2002, 217–218). Furthermore, an
imputation model that is congenial to a given analysis model may nevertheless
fail to produce proper imputations. MI relies on the crucial _missing at
random_ (MAR) assumption, and, when this assumption is violated, MI-based
inferences are invalid. Meeting this assumption requires specifying imputation
models that include all variables that correlate with the missingness and the
analysis model variables. Omitting such a variable from the imputation model
results in imputation under MNAR (Collins ., 2001, 339). MI under MNAR can
lead to bias in the parameter estimates and can invalidate inferences
involving the imputed variables (Collins ., 2001, 341–343), even if the
imputation model is congenial to the analysis model. Therefore, including as
many good predictors of the variables under imputation as possible in the
imputation model is generally advisable.

We refer to all variables that are not targets of imputation, as potential
auxiliary variables. This set of potential auxiliaries may include important
predictors of missingness, variables that correlate with the imputation
targets, and variables that are not useful for imputation. Discerning which of
the potential auxiliary variables may be useful predictors in the imputation
model can be a daunting task. Following an inclusive approach (i.e., including
numerous auxiliary variables in the imputation model) reduces the chances of
omitting important correlates of missingness, thereby making the MAR
assumption more plausible (Rubin ., 1995, 826–827; Schafer, 1997, 23; White .,
2011; Van Buuren, 2018, 167). Furthermore, Collins . (2001) showed that the
inclusive strategy reduces estimation bias and increases efficiency. When
designing broad and intermediate imputation models, the inclusive strategy can
also grant congeniality with a wider range of analysis models.

Although following the inclusive strategy may be beneficial for the imputation
procedure, it is often infeasible to use all potential auxiliaries as
predictors with standard imputation methods. Standard imputation methods, such
as imputation under the normal linear model (Van Buuren, 2018, 68), face
computational limitations in the presence of many predictors. For example,
using traditional (unpenalized) regression models for the imputation model
requires the number of predictors ($p$) in the imputation models to be smaller
than the number of observed cases ($n$) to avoid mathematical singularity of
the underlying system of equations (James ., 2013, 203). As a result, imputers
need to balance the benefits of the inclusive strategy with its computational
limit actions. The large number of variables available in modern social
scientific data sets makes the difficult step of deciding which predictors to
include in the imputation models even more arduous.

In addition to their size, other aspects of social surveys and other social
scientific data can further complicate the task of specifying good imputation
models. Sociologists and researchers working with large social surveys often
want to estimate analysis models that use composite scores (i.e., aggregates
of multi-item scales). When working with multi-item scales, the imputer needs
to decide if variables should be imputed at the item level or at the scale
level. When all a scale’s items are usually missing or observed together,
scale-level imputation can be effective (Mainzer ., 2021). When item-level
missing predominates, however, the literature generally suggests imputing
multi-item scales at the item level (Van Buuren, 2010; Gottschall ., 2012;
Eekhout ., 2014), but pursuing such a strategy can lead to increased
dimensionality of the imputation models (Eekhout ., 2018).

Furthermore, social surveys are often longitudinal, and it is usually most
convenient to impute such a data structure in wide format (Van Buuren, 2018,
312). A wide dataset has a single record for each unit, with observations made
at subsequent time points coded as additional columns in the dataset. As a
result, long-running panel studies might easily induce large pools of
potential auxiliary variables with which the imputer must contend.

###  1.3 High-dimensional imputation

The factors discussed above—or combinations thereof—may result in _high-
dimensional_ imputation problems wherein the pool of potential auxiliary
variables is larger than the available sample size. Such high-dimensional
problems preclude a straightforward application of MI and force researchers to
choose which variables to include in the imputation model or otherwise
regularize the imputation model. One possible solution to this problem is
using high-dimensional prediction models as the imputation model. When we say
“high-dimensional prediction”, we are referring to the branch of statistical
prediction concerned with improving prediction in situations where the number
of predictors is larger than the number of observed cases (the so-called $p>n$
problem). Recent developments in high-dimensional imputation techniques
leverage high-dimensional prediction methodology to offer opportunities for
embracing an inclusive strategy while substantially diminishing its downsides.

MI has been combined with high-dimensional prediction models in algorithms
that use shrinkage methods (Zhao Long, 2016; Deng ., 2016) and dimensionality
reduction (Song Belin, 2004; Howard ., 2015) to avoid the obstacles of an
inclusive strategy. Tree-based imputation strategies (Burgette Reiter, 2010;
Doove ., 2014) also have the potential to overcome the computational
limitations of the inclusive strategy. The nonparametric nature of decision
trees bypasses the identification issues most parametric methods face in high-
dimensional contexts. To the best of our knowledge, no study to date has
directly compared the performance of the various high-dimensional MI methods
recommended in the literature.

###  1.4 Scope of the current project

The goal of this project was to compare how different high-dimensional MI (HD-
MI) methods fare when imputing data sets with many variables. In particular,
we were interested in the types of imputation problems that may arise in large
social scientific data sets. Such data sets do not need to be strictly high-
dimensional to be too large for standard MI routines. Even in low-dimensional
settings (i.e., $n<p$), including too many auxiliary variables in the
imputation model can bias analysis model estimates and lead to convergence
problems and other computational issues (Hardt ., 2012). The high-dimensional
imputation approaches we compared in this project can be used to simplify the
process of specifying a good imputation model in both high- and low-
dimensional problems.

We compared seven state-of-the-art HD-MI algorithms in terms of their ability
to support statistically valid analyses. We chose these techniques because
they stood out as the most promising candidates in our review of the HD-MI
literature. The comparison was based on two numerical experiments: a Monte
Carlo simulation study and a resampling study using Wave 7 of EVS. The
simulation study allows us to compare the imputation methods in an artificial
scenario with maximum experimental control. In a simulation study, we are able
to precisely manipulate data features to match our experimental goals because
we define the population model. However, the variables in a simulation study
are usually sampled from simple multivariate distributions with regular,
unrealistic mean and covariance structures. The resampling study allows us to
shed the artifice of the simulation study and compare the methods using real
social scientific data. EVS is a large-scale, cross-national survey on human
values administered in around 50 countries across Europe. The EVS data contain
both numerical and categorical variables associated via a complicated,
heterogeneous covariance structure. Performing a resampling study on this data
set allows us to estimate bias and coverage in a more ecologically
valid—albeit still somewhat artificial—scenario than is possible with a Monte
Carlo simulation study.111We chose to use the EVS data for a resampling study
rather than an applied example because there is no ground truth in an applied
example. The resampling study offers much of the ecological validity of an
applied example with the added benefit of supporting the same types of
generalizable conclusions provided by a simulation study.

The imputation techniques we compared are best suited to data-driven
imputation with an intermediate or broad scope. The potential benefits of HD-
MI methods lie in the automatic imputation model specification that these
techniques offer. Therefore, we focused on data-driven imputation tasks where
the objective is accommodating a wide range of analysis models. However, the
techniques we compared do not exclude the possibility of specifying more
narrowly scoped imputation models. With little tweaking, one can always force
specific variables into the imputation model.

In what follows, we first introduce the missing data treatments that we
compared in our study. Then, we present the methodology and results of the two
numerical experiments, we discuss the implications of the results for applied
researchers, and we provide recommendations. We conclude by discussing the
limitations of the study and suggesting future research directions.

##  2 Imputation methods and algorithms

We use the following notation: scalars, vectors, and matrices are denoted by
italic lowercase, bold lowercase, and bold uppercase letters, respectively. A
scalar belonging to an interval is indicated by $s_{1}\in[s_{2},s_{3}]$, while
a scalar taking the values in a set is represented as
$s_{1}\in\\{s_{2},s_{3}\\}$. We use the scope resolution operator, ::, to
designate a function provided by a specific software package. So, for example,
_mice::quickpred()_ represents the _quickpred()_ function provided by the
_mice_ package.

Consider an $n\times p$ data set, $\bm{Z}$, comprising variables $\bm{z}_{1}$,
$\bm{z}_{2}$, …, $\bm{z}_{p}$. Assume that the first $t$ variables of $\bm{Z}$
have missing values and that these $t$ variables are the targets of
imputation. Denote the columns of $\bm{Z}$ containing $\bm{z}_{1}$ to
$\bm{z}_{t}$ as the $n\times t$ matrix, $\bm{T}$. The remaining $(p-t)$
columns of $\bm{Z}$ contain variables that are not targets of imputation.
These variables constitute a pool of _potential_ auxiliary variables that
could be used to improve the imputation procedure. Let $\bm{A}$ be a
$n\times(p-t)$ matrix denoting this set of potential auxiliary variables and
write $\bm{Z}$ as $\bm{Z}=(\bm{T},\bm{A})$. For a given $\bm{z}_{j}$, with
$j=(1,\dots,p)$, denote its observed and missing components as
$\bm{z}_{j,obs}$ and $\bm{z}_{j,mis}$, respectively. Let
$\bm{Z}_{-j}=(\bm{z}_{1},\dots,\bm{z}_{j-1},\bm{z}_{j+1},\dots,\bm{z}_{p})$ be
the collection of $p-1$ variables in $\bm{Z}$ excluding $\bm{z}_{j}$. Denote
by $\bm{Z}_{-j,obs}$ and $\bm{Z}_{-j,mis}$ the components of $\bm{Z}_{-j}$
corresponding to the data units in $\bm{z}_{j,obs}$ and $\bm{z}_{j,mis}$,
respectively.

###  2.1 Multiple imputation by chained equations

Assume that $\bm{Z}$ is the result of $n$ random samples from a multivariate
distribution defined by an unknown set of parameters $\bm{\theta}$. The
multiple imputation by chained equations (MICE) approach obtains the posterior
distribution of $\bm{\theta}$ by sampling iteratively from conditional
distributions of the form
$P(\bm{z}_{1}|\bm{Z}_{-1},\bm{\theta}_{1}),\dots,P(\bm{z}_{t}|\bm{Z}_{-t},\bm{\theta}_{t})$,
where $\bm{\theta}_{1},\dots,\bm{\theta}_{t}$ are imputation model parameters
specific to the conditional distributions of each variable with missing
values.

More precisely, the MICE algorithm takes the form of a Gibbs
sampler222Technically, the MICE algorithm is only a true Gibbs sampler when
there exists a valid joint distribution for all targets of imputation. When
such a joint distribution does not exist, the MICE algorithm is still valid,
but it is not a true Gibbs sampler.. At the $m$th iteration $(m=1,\dots,M)$,
samples are drawn for the $j$th target variable ($j=1,\dots,t$) from the
following distributions:

$\displaystyle\hat{\bm{\theta}}_{j}^{(m)}$ $\displaystyle\sim
p(\bm{\theta}_{j}|\bm{z}_{j,obs},\dot{\bm{Z}}_{-j,obs}^{(m)})$ (1)
$\displaystyle\bm{z}_{j,mis}^{(m)}$ $\displaystyle\sim
p(\bm{z}_{j,mis}|\dot{\bm{Z}}_{-j,mis}^{(m)},\hat{\bm{\theta}}_{j}^{(m)}),$
(2)

where $\hat{\bm{\theta}}_{j}^{(m)}$ and $\bm{z}_{j,mis}^{(m)}$ are draws from
the parameter’s full conditional posterior distribution (1) and the missing
data posterior predictive distribution (2), respectively.
$\dot{\bm{Z}}_{-j,obs}^{(m)}$ and $\dot{\bm{Z}}_{-j,mis}^{(m)}$ are subsets of
the variables in $\bm{Z}_{-j}^{(m)}$ (potentially every variable in
$\bm{Z}_{-j}^{(m)}$). These subsets are chosen by the imputer to act as
predictors in the elementary imputation model for $\bm{z}_{j}$.

After convergence, $d$ sets of values are sampled from (2) and used as
imputations. Any analysis model can then be estimated on each of the $d$
completed data sets, and the parameter estimates can be pooled using Rubin’s
rules (Rubin, 1987). The high-dimensional imputation methods we considered in
this paper all follow the general MICE framework, but they differ in which
elementary imputation methods they use to define equations (1) and (2).

####  2.1.1 MICE with a fixed ridge penalty

A linear regression model is the standard univariate imputation model for
MICE. However, ordinary linear regression faces computational limitations when
applied to datasets with many predictors. If $n$ is not much larger than $p$,
the regression estimates will have a large variance, and, if $p>n$ there is no
unique solution for the regression coefficients. Ridge regression (Hastie .,
2009, 61–68) is a common _shrinkage method_ that imposes a penalty during
model estimation to shrink the regression slopes toward zero and allow a large
number of predictors to be included in the model, while still controlling the
variance of the estimates. When applied to the imputation model in MICE, a
ridge penalty allows a more inclusive auxiliary variable strategy.

MICE with a fixed ridge penalty uses the Bayesian normal linear model
described by Van Buuren (2018, 68, algorithm 3.1) as the univariate imputation
method. We refer to this method as Bayesian Ridge (BRidge). In this approach,
the sampling of each $\hat{\bm{\theta}}_{j}^{(m)}$ in equation (1) relies on
inverting the cross-products matrix of $\dot{\bm{Z}}_{-j,obs}^{(m)}$333When
estimating ridge regression coefficients, predictors are centered and scaled
to have unit variance. Adding a positive constant (the ridge penalty,
$\kappa$) to the diagonal of the cross-product matrix stabilizes this
inversion. Indeed, if $p>n$, sufficiently large values of $\kappa$ will
facilitate inversion of the cross-products matrix and induce a unique (albeit
biased) solution for the regression coefficients.

In BRidge, every variable in $\bm{Z}_{-j}$ is used as a predictor in the
imputation model, and the ridge penalty is the only precaution taken to
address a large number of predictors. The value of $\kappa$ is usually chosen
to be close to zero (e.g. $\kappa=0.0001$), because values larger than $0.1$
may introduce excessive systematic bias (Van Buuren, 2018, 68). However,
larger values of $\kappa$ may be necessary to adequately stabilize the
estimation in certain scenarios. In the present work, we chose the value of
$\kappa$ by means of cross-validation.

####  2.1.2 Direct use of regularized regression

444The DURR approach is now implemented in the R package _mice_ via the
mice::mice.impute.lasso.norm() and mice::mice.impute.lasso.logreg() elementary
imputation methods.

Lasso regression (_least absolute shrinkage and selection operator_ ;
Tibshirani, 1996) is another popular shrinkage method. Unlike ridge
regression, the lasso penalty achieves both shrinkage and automatic variable
selection (whereas ridge does not exclude any variables). The extent of the
lasso penalization depends on a tuning parameter, $\lambda$, which is selected
from a set of possible values by means of cross-validation. For sufficiently
large values of $\lambda$, lasso will force some coefficient estimates to be
exactly zero thereby excluding the associated predictors from the fitted
model. When applied to an imputation model, lasso will automatically select
which predictors enter the imputation model. Zhao Long (2016) and Deng .
(2016) used lasso regression as the univariate imputation model in a MICE
algorithm to impute high-dimensional data and referred to this approach as
_direct use of regularized regression_ (DURR).

At iteration $m$, for a target variable $\bm{z}_{j}$, DURR replaces equations
(1) and (2) with the following two steps:

  1. 1.

Generate a bootstrap sample $\bm{Z}^{*(m)}$ by sampling with replacement from
$\bm{Z}^{(m)}$, and train a regularized linear regression model (such as lasso
regression) with $\bm{z}_{j,obs}^{*(m)}$ as outcome and
$\bm{Z}_{-j,obs}^{*(m)}$ as predictors555As with ridge regression, predictors
are centered and scaled to have unit variance.. This produces a set of
parameter estimates (regression coefficients and error variance),
$\hat{\bm{\theta}}_{j}^{(m)}$, that can be viewed as a sample from equation
(1).

  2. 2.

Use $\bm{Z}^{(m)}_{-j,mis}$ and $\hat{\bm{\theta}}_{j}^{(m)}$ to predict
$z_{j,mis}$, and obtain draws from the posterior predictive distribution of
the missing data as in equation (2).

Hence, at every iteration, each elementary imputation model is estimated as a
lasso regression, and uncertainty regarding the parameter values is included
by bootstrapping.

In high-dimensional cases, lasso selects at most $n$ predictors (Zou Hastie,
2005). So, when using lasso for imputation, no elementary imputation model
will contain more predictors than the number of observed cases on the
corresponding outcome. Deng . (2016) compared lasso with the elastic net—which
does not have this restriction—for high-dimensional MI, but they did not find
evidence to favor the elastic net over lasso. Lasso is also computationally
simpler than the elastic net because lasso only has one tuning parameter to
estimate whereas the elastic net has two. Therefore, we chose to implement
DURR with lasso as the regularization method.

####  2.1.3 Indirect use of regularized regression

666The IURR approach is now implemented in the R package _mice_ via the
mice::mice.impute.lasso.select.norm() and
mice::mice.impute.lasso.select.logreg() elementary imputation methods.

While DURR simultaneously performs model regularization and parameter
estimation in equation (1), the _indirect use of regularized regression_
(IURR; Zhao Long, 2016; Deng ., 2016) algorithm uses regularized regression
exclusively for variable selection. The selected predictors are then used in a
standard MI procedure.

At iteration $m$, the IURR algorithm performs the following steps for each
target variable, $\bm{z}_{j}$:

  1. 1.

Fit a linear regression model using a regularized method that does variable
selection (e.g., lasso). Take $\bm{z}_{j,obs}$ as the dependent variable and
$\bm{Z}_{-j,obs}^{(m)}$ as the predictors (unlike DURR, IURR uses the original
data, not a bootstrap sample). The regression coefficients that are _not_
shrunk to 0 define the active set of variables that will be used as predictors
in the actual imputation model (i.e., the variables in
$\dot{\bm{Z}}_{-j}^{(m)}$).

  2. 2.

Obtain the maximum likelihood estimates of the regression coefficients and the
error variance from the linear regression of $\bm{z}_{j,obs}$ onto the active
set of predictors defined in step 1. Then, sample new values of these
parameters from a multivariate normal distribution parameterized by the
MLEs777The sampling notation is the same used by Deng . (2016).:

$(\hat{\bm{\theta}}_{j}^{(m)},\hat{\sigma}_{j}^{(m)})\sim
N(\hat{\bm{\theta}}_{MLE}^{(m)},\hat{\bm{\Sigma}}_{MLE}^{(m)})$ (3)

so that equation (3) corresponds to equation (1) in the general MICE
framework.

  3. 3.

Impute $\bm{z}_{j,mis}$ by sampling from the posterior predictive distribution
based on $\dot{\bm{Z}}_{-j,mis}^{(m)}$ and the parameters’ posterior draws,
$(\hat{\bm{\theta}}_{j}^{(m)},\hat{\sigma}_{j}^{(m)})$.

####  2.1.4 MICE with Bayesian lasso

Zhao Long (2016) proposed the MICE with Bayesian Lasso imputation algorithm
(BLasso), an MI procedure that uses the Bayesian lasso as its elementary
imputation method: MICE with Bayesian lasso (BLasso). A Bayesian lasso model
is a regular Bayesian multiple regression model with informative priors on the
slope coefficients that allow interpreting the mode of the slopes’ posterior
distribution as lasso estimates (Park Casella, 2008; Hans, 2009). Following
Zhao Long (2016), we used the Bayesian lasso specification given by Hans
(20102). Given data with a sample size, $n$, a dependent variable, $\bm{y}$,
and a set of predictors, $\bm{X}$, the Bayesian lasso model has the following
form.

$\displaystyle p(\bm{y}|\bm{\beta},\sigma^{2},\tau)$
$\displaystyle=\textrm{N}(\bm{y}|\bm{X\beta},\sigma^{2}\bm{I}_{n})$ (4)
$\displaystyle p(\beta_{j}|\tau,\sigma^{2},\rho)$
$\displaystyle=(1-\rho)\delta_{0}(\beta_{j})+\rho\left(\frac{\tau}{2\sigma}\right)\times\textrm{exp}\left(\frac{-\tau\left\lVert\beta_{j}\right\rVert_{1}}{\sigma}\right)$
(5) $\displaystyle\sigma^{2}$ $\displaystyle\sim\textrm{Inverse-Gamma}(a,b)$
(6) $\displaystyle\tau$ $\displaystyle\sim\textrm{Gamma}(r,s)$ (7)
$\displaystyle\rho$ $\displaystyle\sim\textrm{Beta}(g,h)$ (8)

Equation (4) represents the density function of a multivariate normal random
variable with mean $\bm{X\beta}$ and covariance matrix $\sigma^{2}\bm{I}_{n}$
evaluated at $\bm{y}$. Equation (5) is the mixture prior distribution for the
regression coefficients $\beta_{j}$ proposed by Hans (20102). This formulation
differs from the classical Bayesian lasso prior proposed by Park Casella
(2008) because of the presence of the sparsity parameter, $\rho$ (Ley Steel,
2009, 655–656; Scott Berger, 2010, 2592), and the point mass at zero,
$\delta_{0}(\beta_{j})$. Finally, equations (6) to (8) represent hyper priors
for the residual variance, $\sigma^{2}$, the penalty parameter, $\tau$, and
the sparsity parameter, $\rho$. Our implementation of BLasso imputation
replaced equation (1) with the BLasso model defined by equations (4) to (8)
with $\bm{y}=\bm{z}_{j,obs}$ and $\bm{X}=\bm{Z}_{-j,obs}$.

The R code used to perform the BLasso imputation was based on the R Package
_blasso_ (Hans, 20101) and can be found on the main author’s GitHub page. For
a detailed description of the Bayesian lasso MI algorithm in a univariate
missing data context see Zhao Long (2016).

####  2.1.5 MICE with PCA

By extracting principal components (PCs) from the set of potential auxiliary
variables, $\bm{A}$, the _MICE with PCA_ (MI-PCA) method summarizes the
information contained in $\bm{A}$ with just a few components. These PCs can
then be used as predictors in a standard, low-dimensional application of MICE.
The MI-PCA procedure can be summarized as follows:

  1. 1.

Extract the first PCs that cumulatively explain the desired proportion of the
variance in the set of potential auxiliary variables, $\bm{A}$888The columns
of $\bm{A}$ are centered and scaled to have unit variance., and collect these
components in a new matrix, $\bm{A}^{\prime}$;

  2. 2.

Replace $\bm{A}$ in $\bm{Z}$ with $\bm{A}^{\prime}$ to obtain
$\bm{Z}^{\prime}=(\bm{T},\bm{A}^{\prime})$;

  3. 3.

Use the standard MICE algorithm with a Bayesian normal linear model and no
ridge penalty to obtain multiply imputed datasets from $\bm{Z}^{\prime}$.

The MI-PCA method was inspired by Howard . (2015) and the _PcAux_ R package
(Lang ., 2018). For this study, we used the R function _stats::prcomp()_ to
perform the PCA estimation via truncated singular value decomposition. Hence,
$p>n$ data are not a problem. When $\bm{A}$ has more columns than rows,
_prcomp()_ will simply extract a maximum of $n$ components.

####  2.1.6 MICE with classification and regression trees

MICE with classification and regression trees (MI-CART; Burgette Reiter, 2010)
is a MICE algorithm that uses classification and regression trees (CART) as
the elementary imputation method. Given an outcome variable $\bm{y}$ and a set
of predictors $\bm{X}$, CART is a nonparametric recursive partitioning
technique that models the relationship between $\bm{y}$ and $\bm{X}$ by
sequentially splitting observations into subsets of units with relatively more
homogeneous $\bm{y}$ values. At every splitting stage, the CART algorithm
searches through all variables in $\bm{X}$ to find the best binary
partitioning rule to predict $\bm{y}$. The resulting collection of binary
splits can be visually represented by a decision tree structure where each
terminal node (or _leaf_) represents the conditional distribution of $\bm{y}$
for units that satisfy the splitting rules.

For each $\bm{z}_{j}$, the $m$th iteration of MI-CART proceeds as follows:

  1. 1.

Train a CART model to predict $\bm{z}_{j,obs}$ from the corresponding
$\bm{Z}_{-j,obs}^{(m)}$.

  2. 2.

Assign each element of $\bm{z}_{j,mis}$ to a terminal node by applying the
splitting rules from the fitted CART model to $\bm{Z}^{(m)}_{-j,mis}$.

  3. 3.

Create imputations for each element of $\bm{z}_{j,mis}$ by sampling from the
pool of $\bm{z}_{j,obs}$ in the terminal node containing $\bm{z}_{j,mis}$.
This procedure corresponds to sampling from the missing data posterior
predictive distribution in Equation (2).

This approach does not consider uncertainty in the imputation model parameters
since the tree structure is not perturbed between iterations. Therefore, MI-
CART cannot produce proper imputations in the sense of Rubin (1986). The
implementation of MI-CART used in this paper corresponds to the one presented
by Doove . (2014, 95, algorithm 1) and the _impute.mice.cart()_ function from
the _mice_ package.

CART searches for the best splitting criterion one variable at a time. As a
result, $p>n$ does not pose the same computational limitations that plague
methods based on linear regression. More variables can increase estimation
times but will not result in computational obstructions.

####  2.1.7 MICE with random forests

MICE with random forests (MI-RF) is a MICE algorithm that uses random forests
as the elementary imputation method. The random forest algorithm (e.g., Hastie
., 2009, 588) entails fitting many decision trees (e.g., CART models) to
subsamples of the original data. These subsamples are derived by resampling
rows with replacement and sampling subsets of columns without replacement. The
random forest algorithm results in an ensemble of fitted decision trees that
generate a sample of predictions for each outcome value. Random forests are
generally considered to improve on the use of individual trees by reducing the
variance of an estimated prediction function.

For each $\bm{z}_{j}$, the $m$th iteration of MI-RF proceeds as follows:

  1. 1.

Generate $k$ bootstrap samples from $\bm{Z}_{-j,obs}$.

  2. 2.

Use these bootstrap samples to fit $k$ single trees predicting
$\bm{z}_{j,obs}$ from a random subset of the variables in $\bm{Z}_{-j,obs}$.

  3. 3.

Generate a pool of $k$ terminal nodes for each element of $\bm{z}_{j,mis}$ by
applying the splitting rules from each of the $k$ fitted trees to the
appropriate columns of $\bm{Z}_{-j,mis}$.

  4. 4.

Create imputations for each element of $\bm{z}_{j,mis}$ by sampling from the
$\bm{z}_{j,obs}$ contained in the pool of terminal nodes defined above.

Bootstrapping and random input selection introduce uncertainty regarding the
imputation model parameters (i.e., the tree structure), as required by a
proper MI procedure. For more details on the MI-RF algorithm, see Doove .
(2014, 103). To perform MICE with random forests we used the R function
_mice::impute.mice.rf()_. As with CART, the random forests algorithm is not
subject to computational limitations in high-dimensional problems because
random forests simply aggregate a collection of univariate decision trees.

####  2.1.8 MICE with quickpred

One popular way to select predictors for the imputation model is to use only
those variables whose associations with the variables under imputation, or
their response indicator, exceed some threshold. To implement this approach,
we used the R function _mice::quickpred()_ with the default value of $r=0.1$
as the screening threshold. We refer to this approach as MI-QP.

The MI-QP approach has two main drawbacks. First, it can still select
collinear predictors. If one predictor is highly correlated with another and
with a variable under imputation, the collinear predictors will both be
selected. Second, when applied to $p>n$ scenarios, MI-QP is not guaranteed to
select fewer predictors than observations available for a given imputation
model. It can still select too many predictors. As a result, MI-QP often needs
to be accompanied by other techniques to address collinearity and linear
dependencies in the data.

####  2.1.9 MICE with analysis model variables as predictors

Another common approach is to use only the analysis model variables in the
imputation model. We refer to this approach as MI-AM. We implemented the MI-AM
method by applying the _mice::mice()_ function to only the analysis model
variables. We left all arguments of the _mice()_ function set to their default
values. As long as the analysis model does not include more variables than the
number of observed cases, MI-AM does not need to address the large
dimensionality of the data.

####  2.1.10 Oracle low-dimensional MICE

When dealing with a large set of possible predictors for the imputation model,
a common recommendation in the MI literature is to decide which predictors to
include according to three criteria (Van Buuren, 2018, 168):

  1. 1.

Include all variables that are part of the analysis model(s).

  2. 2.

Include all variables that are related to the nonresponse.

  3. 3.

Include all variables that are correlated with the targets of imputation.

In practice, the first criterion can be met only if the analysis model is
known before imputation, which is not always true. Furthermore, researchers
can never be sure that the second criterion is entirely met, as there is no
way to know exactly which variables are responsible for missingness. However,
with simulated data, we know which variables define the response model. The
oracle low-dimensional MICE approach (MI-OR) is an ideal specification of the
MICE algorithm that uses this knowledge to include only the relevant
predictors in the imputation models. As such, this method cannot be reliably
used in practice, but it provides a useful reference point for the desirable
performance of an MI procedure. The MI-OR imputations were generated using the
Bayesian normal linear model as the univariate imputation method.

###  2.2 Non-MI strategies

####  2.2.1 Complete case analysis

By default, most data analysis software either fails in the presence of
missing values or defaults to list-wise deletion wherein only complete cases
are used for the analysis (R Core Team, 2020; pandas development team, 2020).
As the default behavior of most analysis tools, complete cases analysis (CC)
remains a popular missing data treatment in the social sciences, despite its
well-known flaws (Rubin, 1987, 8; Van Buuren, 2018, 9 Baraldi Enders, 2010).
Therefore, we included CC as a practical benchmark.

####  2.2.2 Gold standard

We estimated the analysis models directly on the fully observed data before
imposing any missing values. In the following, we refer to the results
obtained in this fashion as the gold standard (GS). These results represent
the counterfactual analysis that would have been performed if there had been
no missing data.

##  3 Simulation study

We investigated the performance of the methods described above with a Monte
Carlo simulation study. Following a similar procedure to that employed by
Collins . (2001), we generated $S=1000$ samples of $n=200$ units while varying
two experimental factors: the number of variables in the dataset,
$p\in\\{50,500\\}$, and the proportion of missing cases on each of $t$
variables, $pm\in\\{0.1,0.3\\}$. Table 1 summarizes the four resulting crossed
conditions.

We chose the values of $n$ and $p$ to reflect extreme dimensionality
situations that would tease apart the relative strengths and weaknesses of the
imputation methods considered here. Nonetheless, we selected these values to
be somewhat plausible for real-world social scientific studies. Consider, for
example, that a typical EVS wave has around 55,000 observations and 250 items
in its questionnaire. Therefore, data structures similar to those in both our
low- and high-dimensional conditions could arise by taking reasonable subsets
of EVS data (potentially over several waves). As for the levels of $pm$, we
chose the lower level to match the 10% of missing cases that is typical of
variables in EVS data. We also included a more extreme level to create more
challenging—but still realistic—conditions for the imputation methods. For
every iteration, we imposed missing values on six target items, and then we
used all missing data treatment methods described above to obtain estimates of
the means, variances, and covariances of these incomplete variables.

condition | label | n | p | pm
---|---|---|---|---
1 | low-dim-low-pm | 200 | 50 | 0.1
2 | high-dim-low-pm | 200 | 500 | 0.1
3 | low-dim-high-pm | 200 | 50 | 0.3
4 | high-dim-high-pm | 200 | 500 | 0.3
|  |  |  |
Table 1: Summary of conditions for Experiment 1. Low-dim (high-dim) represent
conditions where the number of predictors is smaller (larger) than the number
of observations available. Low-pm (high-pm) represent conditions where the
proportion of missing values is low (high)

###  3.1 Simulation study procedure

####  3.1.1 Data generation

At every replication, a data matrix $\bm{Z}_{n\times p}$ was generated
according to a multivariate normal model with means of five and variances of
one. The distribution was centered around five as typical 10-point numerical
items in the EVS dataset have means around five. After sampling the data, all
variables were rescaled to have a variance of approximately 5, which reflects
the expected size of the variance of 10-point items in the EVS data. For the
correlation structure, we defined three blocks of variables based on three
levels of strengths of associations: strong, weak, and no correlation. The
first five variables were highly correlated among themselves ($\rho=0.6$);
variables 6 to 10 were weakly correlated with variables in block 1 and among
themselves ($\rho=0.3$), and all the remaining $p-10$ variables were
uncorrelated. Survey data have more complex correlation structures. However,
when performing imputation of survey data, the challenge is often finding a
few important auxiliary predictors in a large collection of possible
predictors. The block structure and the three correlation levels we chose
should be sufficient to replicate this characteristic.

####  3.1.2 Missing data imposition

Missing values were imposed on six items in $\bm{Z}$: three variables in the
block of highly correlated variables $\\{\bm{z}_{1},\bm{z}_{2},\bm{z}_{3}\\}$
and three in the block of lowly correlated variables
$\\{\bm{z}_{6},\bm{z}_{7},\bm{z}_{8}\\}$. Item nonresponse was imposed by
sampling from a Bernoulli distribution with individual probabilities of
nonresponse defined by

$p_{miss}=p(z_{i,j}=miss|\tilde{\bm{Z}})=\frac{exp(\gamma_{0}+\tilde{\bm{z}}_{i}\bm{\gamma})}{1+exp(\gamma_{0}+\tilde{\bm{z}}_{i}\bm{\gamma})},$
(9)

where $z_{i,j}$ is the $i$th subject’s response on $\bm{z}_{j}$,
$\tilde{\bm{z}}_{i}$ is a vector of responses to the set of missing data
predictors for the $i$th individual, $\gamma_{0}$ is an intercept parameter,
and $\bm{\gamma}$ is a vector of slope parameters. $\tilde{\bm{Z}}$ was
specified to include two fully observed variables from the strongly correlated
set and two from the weakly correlated set
$\\{\bm{z}_{4},\bm{z}_{5},\bm{z}_{9},\bm{z}_{10}\\}$. Therefore, the
probability of nonresponse for a variable depended on variables present in the
data, but never on the variable itself. As a result, when the elements of
$\tilde{\bm{Z}}$ are included as predictors in the MI procedures, the MAR
assumption is satisfied. All slopes in $\bm{\gamma}$ were fixed to 1, while
the value of $\gamma_{0}$ was chosen with an optimization algorithm that
minimized the difference between the actual and desired proportion of missing
values 999The pseudo R-squared for the logistic regression of the missing
value indicator on the predictors of missingness was around 15%. The AUC for
the logistic regression was around 0.72. .

####  3.1.3 Imputation

Missing values were treated with all methods described in the second section.
To evaluate the convergence of the imputation models, we ran ten replications
of the _high-dim-high-pm_ condition and generated trace plots of the imputed
values’ means. These plots suggested that the imputation algorithms converged
within 50 iterations. The only exception was BLasso, which required
approximately 2,000 iterations for convergence.

The ridge penalty used in the BRidge algorithm was fixed across iterations. We
chose the value of the penalty term employing a cross-validation procedure,
wherein penalty values of $10^{-1},10^{-2},\dots,10^{-8}$ were used to impute
data with BRidge. We then selected the value that resulted in the smallest
average fraction of missing information (FMI; Rubin, 1987, eq. 3.1.10) across
the analysis model parameters. Both IURR and DURR could have been implemented
with a variety of penalties (e.g., lasso, Tibshirani, 1996; elastic net, Zou
Hastie, 2005; adaptive lasso, Zou, 2006). In this study, we used lasso as it
is computationally efficient, and it performed well for imputation in Zhao
Long (2016) and Deng . (2016). A 10-fold cross-validation procedure was used
at every iteration of DURR and IURR to choose the penalty parameter. To
maintain consistency with previous research, we specified the BLasso hyper-
parameters in equations (6), (7), and (8) as in Zhao Long (2016):
$(a,b)=(0.1,0.1)$, $(r,s)=(0.01,0.01)$, and $(g,h)=(1,1)$. For the MI-PCA
algorithm, the set of possible auxiliary variables in $\bm{A}$ was defined by
all the fully observed variables. Another important decision when using PCA is
the number of components to keep. Howard . (2015) used only the first
component in their simulations. Since this component explained, on average,
40% of the variance in the auxiliary data, they recommend using enough
components to explain 40% of the variance. For our study, we generated more
complex data for which a single component was not likely to suffice. We,
therefore, applied the intuitively appealing—albeit arbitrary—heuristic of
using enough components to explain 50% of the total variance in the data.

Running MI-QP in the high-dimensional procedure led to frequent convergence
failures. A more common use of the method includes accompanying the quickpred
approach with a ridge penalty and data-driven checks that exclude collinear
variables. We decided to run MI-QP in this more favorable manner by applying
the mice package’s usual data-screening procedures. Accordingly, the _mice()_
call for MI-QP was specified with default values for the following arguments:
ridge = $1e-5$, eps = $1e-4$, threshold = $0.999$.

####  3.1.4 Analysis and comparison criteria

After imputation, we estimated the item means, variances, and covariances of
the six variables with missing values. These six variables were defined as the
analysis model variables. After pooling, we compared the performances of the
methods by computing the bias and confidence interval coverage for the six
mean, six variance, and 15 covariance estimates.

For a given parameter of interest $\theta$, we used the absolute percent
relative bias (PRB) to quantify the estimation bias introduced by the
imputation procedure:

$\textit{PRB}=\displaystyle\left\lvert\
\frac{\bar{\hat{\theta}}-\theta}{\theta}\right\rvert\times 100,$ (10)

where $\theta$ is the true value of the focal parameter defined as
$\sum_{s=1}^{S}\hat{\theta}_{s}^{GS}/S$, with $\hat{\theta}_{s}^{GS}$ being
the Gold Standard parameter estimate for the $s$th repetition. The averaged
focal parameter estimate under a given missing data treatment was computed as
$\bar{\hat{\theta}}=\sum_{s=1}^{S}\hat{\theta}_{s}/S$, with $\hat{\theta}_{s}$
being the estimate obtained from the treated incomplete data in the $s$th
repetition. Following Muthén . (1987), we considered $\text{PRB}>10$ as
indicative of problematic estimation bias.

To assess the performance in hypothesis testing and interval estimation, we
evaluated the confidence interval coverage (CIC) of the true parameter value:

$\textit{CIC}=\frac{\sum_{s=1}^{S}I(\theta\in\widehat{\textit{CI}}_{s})}{S},$
(11)

where $\widehat{\textit{CI}}_{s}$ is the confidence interval of the parameter
estimate $\hat{\theta}_{s}$ in the $s$th repetition, and $I(.)$ is the
indicator function that returns 1 if the argument is true and 0 otherwise.

CICs below 0.9 are usually considered problematic for 95% CIs (Van Buuren,
2018, 52; Collins ., 2001, 340) as they imply inflated Type I error rates.
High CICs (e.g., above 0.99) indicate CIs that are too wide, implying inflated
Type II error rates. Therefore, we considered CIs to show severe under-
coverage (over-coverage) if CIC $<0.9$ (CIC $>0.99$). From a testing
perspective, a CIC can be considered as significantly different from the
nominal coverage rate if the magnitude of its difference from the nominal
coverage proportion ($p_{0}$) is more than two times the standard error of
$p_{0}$, $\textit{SE}(p_{0})=\sqrt{p_{0}(1-p_{0})/S}$ (Burton ., 2006). In our
simulation study, the nominal coverage probability was 95%. Therefore, we
considered 95% CI coverages outside the interval $[0.94,0.96]$ to be
significantly different from the nominal coverage rate. We assumed normal
sampling distributions for variances and covariances when computing and
pooling their CIs. This assumption is plausible under large sample conditions.

###  3.2 Results

We computed both PRB and CIC for each of the 27 parameters in the analysis
model (six means, six variances, and 15 covariances). To summarize the
results, we focus on the expected and extreme values of these measures. In
Figures 1 and 2, we report the average, minimum, and maximum PRB and CIC
obtained with the different missing data treatments, for each parameter type.
As the GS estimates were used to define the “true” values of the parameters,
the bias for this method was by definition 0. So, we do not include bias of
the GS estimates in the figure. For ease of presentation, we report the
results only for the large proportion of missing cases ($pm=0.3$) condition.
While the relative performances were independent of the missing data rate, the
performance patterns were clearer with a larger proportion of missing values.
In the supplementary material, we included the same figures for the low
proportion of missing cases ($pm=0.1$) condition.

####  3.2.1 Means

The largest PRB for the means was below 10 for all imputation methods. Only CC
produced problematic degrees of bias. Looking at the relative performances,
IURR, BRidge, MI-PCA, and MI-OR resulted in smaller biases than the other
methods. In terms of CIC, only MI-PCA and MI-OR showed a consistently strong
performance. Neither method demonstrated any extreme under-/over-coverage
(i.e., all CICs $\in$ [0.9, 0.99]), and both methods resulted in only the
highest coverage being significantly different from nominal coverage (max CIC
$>0.96$).

IURR resulted in significant under-coverage of the true means (CICs $<0.94$)
in both the high-dimensional ($p=500$) and low-dimensional ($p=50$)
conditions, although under-coverage was never severe (with CICs
$\in[0.90,0.94]$). DURR and BLasso demonstrated some significant differences
from nominal coverage in the low-dimensional condition, and both led to
extreme under-coverage in the high-dimensional condition. The tree-based
methods and CC performed most poorly. These methods led to CICs significantly
different from nominal coverage rates in all conditions, and they demonstrated
extreme under-coverage even in the low-dimensional condition. MI-QP resulted
in close to nominal coverage in the low-dimensional condition, performing
about as well as MI-OR and MI-PCA. However, in virtually all replications of
the high-dimensional condition, the CIs contained the true parameter values,
thereby producing severe over-coverage. Finally, MI-AM resulted in significant
to extreme under-coverage. This method was not influenced by the
dimensionality of the data as it used the same six variables as predictors in
both conditions.

####  3.2.2 Variances

IURR, BLasso, and the tree-based MI methods resulted in low biases (i.e., PRB
$<10$) in both the high and low dimensional conditions. For BLasso, these low
biases were paired with low deviations from nominal coverage rates. IURR only
demonstrated problematic CICs for the high-dimensional condition, where it
produced extreme under-coverage (largest CIC $<0.9$). MI-CART and MI-RF did
not produce reasonable coverage in either the low- or the high-dimensional
condition, with the largest coverage being significantly different from
nominal (CICs $<0.94$) and the smallest being severely below the nominal level
(CICs $<0.9$).

MI-PCA showed acceptable biases and reasonable coverage rates in the low-
dimensional condition, but it showed large biases and under-coverage in the
high-dimensional condition. In the low-dimensional condition, DURR produced
low bias and reasonable coverage (i.e., only the lowest coverage being
significantly different from nominal), but it resulted in PRBs $>10$ for all
variances in the high-dimensional condition, where it also produced extreme CI
under-coverage. BRidge and CC performed poorly in nearly all conditions. These
methods tended to demonstrate substantial biases and extreme under-coverage.
Although MI-QP performed well in the low-dimensional condition, in the high-
dimensional condition, it resulted in PRBs larger than 50 and CICs close to 1
for all six item variances. MI-AM maintained low bias and acceptable coverage
for all item variances.

####  3.2.3 Covariances

MI-PCA was the only method that showed consistently strong performance when
estimating covariances. MI-PCA showed negligible bias and minimal deviations
from nominal coverage in both low- and high-dimensional conditions. In
particular, the PRB was smaller than 10 for all covariances in both conditions
and was almost as low as the PRB obtained by MI-OR. MI-PCA never produced
extreme under-/over-coverage, and when the CIC was significantly different
from the nominal rate, the CIs showed mild _over_ -coverage (i.e., CICs
greater than 0.96 but smaller than 0.99). After MI-PCA, IURR demonstrated the
second-strongest performance, with negligible bias and acceptable coverage in
the low-dimensional condition. However, in the high-dimensional condition,
IURR produced large biases and extreme under-coverage with the average bias
being above the 10% threshold and even the largest coverage being just around
the 90% threshold.

Both MI-QP and BRidge displayed low bias and acceptable coverage in the low-
dimensional condition, but they resulted in unacceptable biases in the high-
dimensional condition. In the high-dimensional condition, MI-QP led to 100%
coverage of the true values, while BRidge led to extreme under-coverage of the
true values. MI-AM, DURR, BLasso, and the tree-based MI methods tended to
result in PRBs larger than 10, accompanied by under-coverage of the true
covariance values, even in the low-dimensional conditions.

Figure 1:  Minimum, average, and maximum absolute percent relative bias (PRB)
for the 6 item means, 6 variances, and 15 covariances in the simulation study.
If no data points are reported for a method in a panel, all of its PRBs were
larger than 50. The methods reported on the Y-axis are: Direct Use of
Regularized Regression (DURR), Indirect Use of Regularized Regression (IURR),
MICE with Bayesian Lasso (BLasso), MICE with Bayesian Ridge (BRidge), MICE
with Principal Component Analysis (MI-PCA), MICE with CART (MI-CART), MICE
with Random Forests (MI-RF), Oracle low-dimensional MICE (MI-OR), MICE with
quickpred (MI-QP), MICE with Analysis Model (MI-AM), and Complete Case
Analysis (CC).  Figure 2:  Minimum, average, and maximum confidence interval
coverage (CIC) for the 6 item means, 6 variances, and 15 covariances in the
simulation study. If no data points are reported for a method in a panel, all
of its CICs were smaller than 0.80. The methods reported on the Y-axis are:
Direct Use of Regularized Regression (DURR), Indirect Use of Regularized
Regression (IURR), MICE with Bayesian Lasso (BLasso), MICE with Bayesian Ridge
(BRidge), MICE with Principal Component Analysis (MI-PCA), MICE with CART (MI-
CART), MICE with Random Forests (MI-RF), Oracle low-dimensional MICE (MI-OR),
MICE with quickpred (MI-QP), MICE with Analysis Model (MI-AM), Complete Case
Analysis (CC), and Gold Standard Analysis (GS).

##  4 EVS resampling study

We performed a resampling study based on the EVS data to assess whether the
results of our simulation study would replicate in more realistic data. EVS is
a high-quality survey widely used by sociologists for comparative studies
between European countries (EVS, 20202). Furthermore, it is freely available
and represents the type of data social scientists regularly analyze. Variables
in the EVS data are discrete numerical and categorical items following a
variety of distributions.

To perform the resampling study, we treated the original EVS data as a
population. We then resampled $S=1000$ data sets of $n$ units from this
population, and we used these replicates as we used the multivariate normal
samples in the simulation study. For each replicate, we imposed missing
values, and we treated these missing values with the same methods explored in
the simulation study. This procedure was repeated for low-dimensional and
high-dimensional conditions. As the number of predictors in the data was fixed
at $p=243$, we controlled the dimensionality of the data by varying the sample
size ($n\in\\{1000,300\\}$). When the sample size was 300, after dummy coding
categorical predictors, even a small proportion of missing values ($pm=0.1$)
led to a high-dimensional ($p>n$) situation. Although $n=300$ might be too low
to represent the typical use of EVS data, we do not see this as a limitation
on the following results, for two reasons. First, our purpose in conducting
this resampling study was primarily to see if our simulation results would
carry over into data generated from a more realistic population model, not
necessarily to see if those results would hold in a typical social science
dataset. Increasing the sample size to match ranges typically seen in analyzes
of EVS data would remove the high-dimensional condition where we saw the most
interesting results in the simulation, thereby greatly reducing the utility of
the resampling study. Second, many social science studies analyze data with
around 300 observations, so our samples are not unrealistic in a general
sense.

###  4.1 Resampling study procedure

####  4.1.1 Data preparation and sampling

We used the third prerelease of the 2017 wave of EVS data (EVS, 20201) to
create a population data set with _no missing values_. The original data set
contained 55,000 observations from 34 countries. We selected only the four
founding countries of the European Union included in the data set (France,
Germany, Italy, and the Netherlands) because keeping all countries would have
entailed either including a set of 33 dummy codes in the imputation models or
imputing under some form of a multilevel model. Since both of these options
fall outside the scope of the current study, we opted to subset the data as
described. We excluded all columns that contained duplicated information
(e.g., recoded versions of other variables), or metadata (e.g., time of the
interview, mode of data collection).

The original EVS data set contained missing values. We needed to treat these
missing data before we could use the EVS data in the resampling study. We used
the _mice_ package to fill the missing values with a single round of
predictive mean matching (PMM). We used the _quickpred_ function, to select
the predictors for the imputation models. We implemented the variable
selection by setting the minimum correlation threshold in _quickpred_ to 0.3.
The number of iterations in the _mice()_ run was set to 200. We used a single
imputation, and not MI, because this imputation procedure was used only to
obtain a set of pseudo-fully observed data to act as the population in our
resampling study and not for statistical modeling, estimation, or inference
with respect to the true population from which the EVS data were sampled. For
the same reason, the relatively poor performance that we observed for MI-QP in
the simulation study is not relevant here. At the end of the data cleaning
process, we obtained a pseudo-fully observed dataset of 8,045 observations
across four countries with $p=243$ variables. For every replicate in the
resampling study, we generated a bootstrap sample by sampling $n$ observations
with replacement from this dataset.

####  4.1.2 Analysis models

To define plausible analysis models, we reviewed the models reported in the
repository of publications using EVS data that is available on the EVS website
(EVS, 20202). As a result, we defined two linear regression models. Model 1
was inspired by Köneke (2014). The dependent variable was a 10-point item
measuring euthanasia acceptance (‘Can [euthanasia] always be justified, never
be justified, or something in between?’). The predictors included an item
measuring the self-reported importance of religion in one’s life, trust in the
health care system, trust in the state, trust in the press, country, sex, age,
education, and religious denomination. A researcher might estimate this model
to test a hypothesis regarding the effect of religiosity on the acceptance of
end-of-life treatments.

Model 2 was inspired by Immerzeel . (2015). The dependent variable was a
harmonized variable that quantifies the respondents’ tendencies to vote for
left- or right-wing parties, expressed on a 10-point left-to-right continuum.
The predictors included a scale measuring respondents’ attitudes toward
immigrants and immigration (‘nativist attitudes scale’). The scale was
obtained by taking the average of respondents’ agreement, on a scale from 1 to
10, with three statements: ‘immigrants take jobs away from natives’,
‘immigrants increase crime problems’, and ‘immigrants are a strain on welfare
system’. The remaining predictors were: attitudes toward law and order,
attitudes toward authoritarianism, interest in politics, level of political
activity, country, sex, age, education, employment status, socioeconomic
status, importance of religion in life, religious denomination, and the size
of the town where the interview was conducted. A researcher might estimate
this model to test a hypothesis regarding the effect of xenophobia on voting
tendencies.

####  4.1.3 Missing data imposition

We imposed missing data on six variables using the same strategy as in the
simulation study. The targets of missing data imposition were the two
dependent variables in models 1 and 2 (i.e., euthanasia acceptance, and left-
to-right voting tendency), religiosity, and the three items making up the
“nativist attitudes” scale. The response model was the same as in Equation
(9), and three variables were included in $\tilde{\bm{Z}}$: age, education,
and an item measuring trust in new people. 101010 The pseudo R-squared for the
logistic regression of the missing value indicator on the predictors of
missingness was around 14%. The AUC for the logistic regression was around
0.75.  We chose these predictors because older people tend to have higher item
nonresponse rates than younger people, and lower educated people tend to have
higher item non-response rates than higher educated people (Guadagnoli Cleary,
1992; De Leeuw ., 2003). We also assumed that people with less trust in
strangers would have a higher nonresponse tendency as they are likely to
withhold more information from the interviewer (a stranger).

####  4.1.4 Imputation

We treated the missing values with the same methods used in the simulation
study. MI-AM used all the variables present in either of the analysis models
as predictors for the imputation models. MI-PCA was performed considering all
the fully observed variables as possible auxiliary variables. In other words,
the six variables with missing values were used in their raw form, while the
remaining 237 were used to extract PCs. The other imputation methods were
parameterized in the same way as in the simulation study, and convergence
checks were performed in the same way. These convergence checks suggested that
the imputation models had converged after 60 iterations.

###  4.2 Results

When estimating linear regression models, all partial regression coefficients
can be influenced by missing values on a subset of the variables included in
the model. Therefore, it is important to evaluate the estimation bias and CIC
rates for all model parameters. Figure 3 reports the absolute PRBs for the
intercept and all partial slopes from Model 2 obtained after using each
imputation method, for both the low- and high-dimensional conditions. Model 2
has an intercept and 13 regression coefficients. Every horizontal line in the
figure represents the PRB for the estimation of one of these 14 parameters.
Figure 4 reports CIC results in the same way. For ease of presentation,
results for Model 1 are reported in the supplementary materials.

As seen in Figure 3, in both the high- and low-dimensional conditions, DURR,
IURR, BLasso, and MI-CART showed only slightly larger PRBs than MI-OR.
However, even MI-OR did not provide entirely unbiased parameter estimates.
After imputing with MI-OR, almost half of the parameters in Model 2 were
estimated with large bias (PRB $>10\%$). MI-PCA, MI-RF, and CC showed similar
trends but produced larger PRBs (particularly CC). BRidge demonstrated the
same results described in the simulation studies. It was competitive in the
low-dimensional scenario, but it was inadequate with high-dimensional data
(all PRBs $>10$.) In the low-dimensional condition, MI-QP resulted in only
three parameter estimates with acceptable bias and only one in the high-
dimensional condition. MI-AM resulted in six parameter estimates with
acceptable bias in the low-dimensional condition but only one in the high-
dimensional condition.

As shown in Figure 4, MI-OR and DURR resulted in the lowest deviations from
nominal coverage, with only one or two coverages differing significantly from
the nominal level. IURR showed a similar trend but four coverages were
significantly different from nominal in the low-dimensional condition.

BLasso, MI-PCA, MI-CART, MI-RF, and MI-AM all showed similar performance in
the low-dimensional condition. These methods all significantly over-covered
most of the parameters but did not produce any extreme under-/over-coverage,
except for one parameter for MI-RF. BLasso, MI-PCA, and MI-RF maintained
similar performance in the high-dimensional condition, but MI-CART improved to
match the performance of MI-OR, and MI-AM produced extreme over-coverage for
most of the parameters. BRidge performed well in the low-dimensional
condition—around the level of IURR—but produced very poor coverages in the
high-dimensional condition. MI-QP performed poorly in both the low- and high-
dimensional conditions, producing only two non-significant coverages in the
low-dimensional condition and none in the high-dimensional condition. CC
performed quite well, but it had a much more pronounced tendency toward
_under_ -coverage than the MI methods. Notably, very few of the CICs fell into
the range of extreme under-/over-coverage. Only the high-dimensional estimates
from BRidge and MI-AM consistently exhibited extreme under-/over-coverage.

Figure 3:  PRB for all the model parameters in model 2. For each method, the
PRBs are ordered by increasing absolute value. The methods reported on the
Y-axis are: Direct Use of Regularized Regression (DURR), Indirect Use of
Regularized Regression (IURR), MICE with Bayesian Lasso (BLasso), MICE with
Bayesian Ridge (BRidge), MICE with Principal Component Analysis (MI-PCA), MICE
with CART (MI-CART), MICE with Random Forests (MI-RF), Oracle low-dimensional
MICE (MI-OR), MICE with quickpred (MI-QP), MICE with analysis model (MI-AM),
and Complete Case analysis (CC).  Figure 4:  CIC for all model parameters in
model 2. For each method, the CICs are ordered by increasing value. The
methods reported on the Y-axis are: Direct Use of Regularized Regression
(DURR), Indirect Use of Regularized Regression (IURR), MICE with Bayesian
Lasso (BLasso), MICE with Bayesian Ridge (BRidge), MICE with Principal
Component Analysis (MI-PCA), MICE with CART (MI-CART), MICE with Random
Forests (MI-RF), MICE with quickpred (MI-QP), MICE with analysis model (MI-
AM), Oracle low-dimensional MICE (MI-OR), Complete Case analysis (CC), and
Gold Standard analysis (GS).

####  4.2.1 Imputation time

Table 2 reports the average imputation time for the different methods. IURR
and DURR were the most time-consuming methods, with imputation times above one
hour in the low-dimensional condition. In the high-dimensional condition, IURR
and DURR were not as time-intensive due to the smaller sample size but still
took more than ten times longer than MI-PCA and BLasso. MI-PCA was the fastest
method, with imputation times of under a minute in both the high- and low-
dimensional conditions. BLasso, MI-OR, and MI-AM were close seconds, with
imputation times of two minutes or less in both conditions. BRidge, MI-CART,
MI-RF, and MI-QP fell in the middle, with imputations times ranging from
1.24(MI-CART) to 12.68(MI-QP) minutes in the high-dimensional condition and
from 3.46(MI-CART) to 9.40(MI-RF) minutes in the low-dimensional condition.

condition DURR IURR BLasso BRidge MI-PCA MI-CART MI-RF MI-OR MI-QP MI-AM low-
dim (n = 1000) 66.929 69.083 1.250 7.068 0.532 3.461 9.400 1.878 9.036 2.082
high-dim (n = 300) 5.437 8.714 0.444 2.767 0.364 1.235 3.968 1.634 12.677
1.899

Table 2: Average imputation time in minutes for the different MI methods in
the resampling study.

##  5 Discussion

###  5.1 Methods that work well

On balance, IURR and MI-PCA were the two strongest performers across the
simulation study and the resampling study. In the simulation study, IURR
produced trivial estimation bias for all parameters in the low-dimensional
condition and means and variances in the high-dimensional condition.
Furthermore, the covariance estimation bias introduced by IURR in the high-
dimensional condition only slightly exceeded the PRB $=10$ threshold, while
most of the other MI methods resulted in covariance PRBs larger than 20 (with
MI-PCA being the most salient exception). IURR produced good coverages in the
low-dimensional condition, but tended to under-cover in the high-dimensional
condition, especially for variances and covariances. In the resampling study,
IURR was also among the strongest performers. Although it did not demonstrate
the best performance, there were no conditions in which IURR produced
unacceptable results.

From the end-user’s perspective, IURR is an appealing method. IURR does not
require the imputer to make choices regarding which variables are relevant for
the imputation procedure. The only additional decision required of the imputer
is selecting the number of folds to use when cross-validating the penalty
parameter. As a result, an IURR imputation run is easy to specify, which makes
IURR an appealing method for the imputation of large social scientific
datasets. However, IURR is relatively computationally intensive. If the number
of variables with missing values is large, IURR might result in prohibitive
imputation time. In such a scenario, the MI-PCA method might be preferable.

In the simulation study, MI-PCA showed small bias and good coverage for both
item means and covariances. Although it exhibited a large bias of the item
variances, the—arguably more interesting—covariance relations between
variables with missing values were always correctly estimated. Notably, MI-PCA
was the only method resulting in small bias and close-to-nominal CIC for the
covariances, even in the high-dimensional condition. Finally, when the CICs
obtained with MI-PCA deviated significantly from nominal rates, they over-
covered. In most situations, over-coverage is less worrisome than under-
coverage as it leads to conservative, rather than liberal, inferential
conclusions. In the resampling study, MI-PCA demonstrated middle-of-the-pack
performance: somewhat worse than IURR, but still within acceptable levels. We
believe the strong performance demonstrated by MI-PCA in the simulation study
outweighs the mediocre performance shown in the resampling study.
Consequently, we view MI-PCA as a promising approach for data analysts
interested in testing theories on large social scientific datasets with
missing values. We are currently conducting research to explore—and extend—the
capabilities of the MI-PCA approach.

###  5.2 Methods with mixed results

In both the simulation study and the resampling study, BRidge manifested the
same mixed performance. This method worked well when the imputation task was
low-dimensional but led to extreme bias and unacceptable CI coverage in nearly
all the high-dimensional conditions. Our results suggest that BRidge is
effective only for low-dimensional imputation problems. The same can be said
of MI-QP, at least as implemented here. In our studies, its performance was
acceptable only in low-dimensional conditions.

DURR performed very well in the resampling study and quite poorly in the
simulation study. In the resampling study, DURR was probably the best overall
method in terms of bias and coverage, but it performed very badly in the high-
dimensional condition of the simulation study. In the simulation study’s low-
dimensional condition, DURR produced a small bias and good CI coverage for
item means and variances. However, compared to IURR, it suffered from greater
deterioration in performance when applied to high-dimensional data, especially
in terms of coverage. Our results suggest that DURR may have some unique
benefits when treating the types of more discrete data seen in the resampling
study. On balance, though, DURR probably should not be preferred to IURR.

There was little difference in performance between the use of CART and random
forests as elementary imputation methods within the MICE algorithm. In line
with what Doove . (2014) found, when a difference was noticeable, the simpler
CART generally outperformed the more complex random forests. Both MI-CART and
MI-RF produced large covariance bias in the simulation study. Although the
bias for means and variances was acceptable, it was usually larger than that
obtained by other MI methods. Furthermore, in terms of CI coverage, both
methods showed a large under-coverage of the true values in the high-
dimensional condition. In the resampling study, MI-CART and MI-RF both showed
somewhat better performance than in the simulation study but not enough better
to outweigh the mediocre simulation study performance. Although the
nonparametric nature of these approaches elegantly avoids over-
parameterization of imputation models, these methods were still outperformed
by IURR and MI-PCA.

In the simulation study, Blasso resulted in small biases for item means and
variances, even in the high-dimensional conditions, but it produced
unacceptably biased covariance in both the low- and high-dimensional
conditions. On the other hand, BLasso seemed to recover the relationships
between variables in the resampling study well, where the overall bias levels
for the regression coefficients were similar to those of MI-OR. However, in
terms of CI coverage, BLasso showed poor performance in both studies resulting
in either under-coverage or over-coverage for most parameters in the high-
dimensional conditions.

The mixed performance of BLasso is also accompanied by a few obstacles to its
application for social scientific research. Using Hans (20102)’s Bayesian
Lasso requires the specification of six hyper-parameters, which introduces
more researcher degrees of freedom and demands a strong grasp of Bayesian
statistics. Furthermore, the method has not currently been developed for
multi-categorical data imputation, a common task in the social sciences. As a
result, we do not recommend BLasso for the imputation of large social science
data sets.

Finally, we do not recommend the use of MI-AM to impute large social science
data sets. MI-AM bypasses the need to select which of the many potential
auxiliary variables should be included in the imputation models. However, it
can lead to the omission of crucial MAR predictors. In our simulation study,
smaller biases and better coverages could always be achieved by using at least
one of the alternative methods we evaluated.

##  6 Limitations and future directions

The present study aimed to compare current implementations of existing
imputation methods. As a result, the scope of the simulation and resampling
studies was limited by the current development state of the different methods.
For example, DURR, IURR, and MI-PCA allow imputation of any type of data: DURR
and IURR have been developed for categorical data imputation (Deng ., 2016),
and MI-PCA can be performed with any standard imputation model for categorical
data. However, BLasso has not been formally developed for imputing multi-
categorical variables yet. This limitation of BLasso forced us to work with
missing values on variables that are either continuous or usually considered
as such in practice (e.g., Likert-type scales). To maintain a fair comparison
with BLasso, all methods were implemented with the assumption that the imputed
variables are continuous and normally distributed. However, IURR, DURR, and
MI-PCA could have performed differently in the resampling study if we had used
their ordinal data implementations.

More generally, the results reported in this article only apply to the
specific implementations of the algorithms we used. Many of the methods
discussed could have been implemented differently. Zhao Long (2016) proposed
versions of IURR and DURR using the elastic net penalty (Zou Hastie, 2005) and
the adaptive lasso (Zou, 2006) instead of the lasso penalty. Although no
substantial performance differences between penalty specifications emerged
from the work of Zhao Long (2016) or Deng . (2016), we must acknowledge that
we did not investigate the impact of different types of regularization in the
present study. Similarly, we have not investigated the sensitivity of BLasso
to different hyper-parameters choices. Finally, the use of random forests
within the MICE algorithm followed Doove . (2014), the version supported in
the popular R package _mice_. However, Shah . (2014) independently developed
another implementation of random forests within the MICE algorithm, which was
available in the now archived R package CALIBERrfimpute (Shah, 2018). We are
not aware of any evidence or theoretical reason to expect differences between
the two implementations, but we did not verify this empirically.

MI-PCA requires choosing the number of components to extract from the
auxiliary variables. In this study, we decided to retain the first components
that explained 50% of the total variance in the auxiliary variables. However,
this decision was arbitrary and its impact was not explored in this study.
Additionally, the good performance of MI-PCA may have been driven by the fact
that, while imputing the $j$th variable, all other variables under imputation
were used directly as predictors. If the other variables under imputation had
been included in the imputation models through the PCs extraction step, and
not used as separate, individual predictors, the performance of MI-PCA might
have been less favorable. We are currently developing a follow-up study to
assess the effects of these two factors on the MI-PCA method.

##  7 Conclusions

Our objective in this project was to find a good data-driven way to select the
predictors that go into an imputation model. A wide range of methods have been
proposed to address this issue, but little research has been done to compare
their performance. With this article, we start to fill this gap and provide
initial insights into applying such methods in social science research. Both
MI-PCA and IURR showed promising performance when compared to other high-
dimensional imputation approaches. While both of these methods represent good
options for automatically defining the imputation models of an MI procedure,
MI-PCA is the more practically appealing option due to its much greater speed.
The use of Bayesian ridge regression is a good alternative when the imputer
wants to have an automatic way of defining the imputation models in a low-
dimensional setting ($n\gg p$).

##  8 Availability of data and materials

The EVS data that support the findings of this study are openly available in
GESIS Data Archive at https://doi.org/10.4232/1.13511, reference number
ZA7500.

##  9 Code availability

The code used for the study is available at the corresponding author’s GitHub
page: [https://github.com/EdoardoCostantini/imputeHD-comp/releases/tag/v1.0].
Please read the README.md files for instructions on how to replicate the
results.

## References

  * Baraldi Enders (2010) baraldiEnders:2010Baraldi, AN. Enders, CK.  2010\.  An introduction to modern missing data analyses An introduction to modern missing data analyses.  Journal of School Psychology4815–37.
  * Brick Williams (2013) brickWilliams:2013Brick, JM. Williams, D.  2013\.  Explaining rising nonresponse rates in cross-sectional surveys Explaining rising nonresponse rates in cross-sectional surveys.  The ANNALS of the American academy of political and social science645136–59.
  * Burgette Reiter (2010) burgetteReiter:2010Burgette, LF. Reiter, JP.  2010\.  Multiple Imputation for Missing Data via Sequential Regression Trees Multiple imputation for missing data via sequential regression trees.  American Journal of Epidemiology17291070–1076.  10.1093/aje/kwq260
  * Burton . (2006) burtonEtAl:2006Burton, A., Altman, DG., Royston, P. Holder, RL.  2006\.  The design of simulation studies in medical statistics The design of simulation studies in medical statistics.  Statistics in medicine25244279–4292.
  * Collins . (2001) collinsEtAl:2001Collins, LM., Schafer, JL. Kam, CM.  2001\.  A comparison of inclusive and restrictive strategies in modern missing data procedures A comparison of inclusive and restrictive strategies in modern missing data procedures.  Psychological Methods64330–351.  10.1037//1082-989X.6.4.330
  * De Leeuw . (2003) leeuwEtAl:2003De Leeuw, ED., Hox, JJ. Huisman, M.  2003\.  Prevention and treatment of item nonresponse Prevention and treatment of item nonresponse.  Journal of Official Statistics19153–176.
  * Deng . (2016) dengEtAl:2016Deng, Y., Chang, C., Ido, MS. Long, Q.  2016\.  Multiple imputation for general missing data patterns in the presence of high-dimensional data Multiple imputation for general missing data patterns in the presence of high-dimensional data.  Scientific reports621689.  10.1038/srep21689
  * Doove . (2014) dooveEtAl:2014Doove, LL., Van Buuren, S. Dusseldorp, E.  2014\.  Recursive partitioning for missing data imputation in the presence of interaction effects Recursive partitioning for missing data imputation in the presence of interaction effects.  Computational Statistics & Data Analysis7292–104.
  * Eekhout . (2018) eekhoutEtAl:2018Eekhout, I., de Vet, HC., de Boer, MR., Twisk, JW. Heymans, MW.  2018\.  Passive imputation and parcel summaries are both valid to handle missing items in studies with many multi-item scales Passive imputation and parcel summaries are both valid to handle missing items in studies with many multi-item scales.  Statistical Methods in Medical Research2741128–1140.
  * Eekhout . (2014) eekhoutEtAl:2014Eekhout, I., de Vet, HC., Twisk, JW., Brand, JP., de Boer, MR. Heymans, MW.  2014\.  Missing data in a multi-item instrument were best handled by multiple imputation at the item score level Missing data in a multi-item instrument were best handled by multiple imputation at the item score level.  Journal of clinical epidemiology673335–342.
  * Enders (2010) enders:2010Enders, CK.  2010\.  Applied missing data analysis Applied missing data analysis.  New York, NYThe Guilford Press.
  * Estefan . (2021) estefanEtAl:2021Estefan, LF., Vivolo-Kantor, AM., Niolon, PH., Le, VD., Tracy, AJ., Little, TD.others  2021\.  Effects of the dating matters® comprehensive prevention model on health-and delinquency-related risk behaviors in middle school youth: A cluster-randomized controlled trial Effects of the dating matters® comprehensive prevention model on health-and delinquency-related risk behaviors in middle school youth: A cluster-randomized controlled trial.  Prevention Science222163–174.
  * EVS (20201) EVS:2017EVS.  20201\.  European Values Study 2017: Integrated Dataset (EVS 2017). European values study 2017: Integrated dataset (evs 2017).  GESIS Data Archive, Cologne. ZA7500 Data file Version 3.0.0.  10.4232/1.13511
  * EVS (20202) EVSbibEVS.  20202\.  EVS Bibliography. Evs bibliography.  Retrieved September 30, 2020 (https://europeanvaluesstudy.eu/education-dissemination-publications/evs-publications/publications/)
  * Gottschall . (2012) gottschallEtAl:2012Gottschall, AC., West, SG. Enders, CK.  2012\.  A comparison of item-level and scale-level multiple imputation for questionnaire batteries A comparison of item-level and scale-level multiple imputation for questionnaire batteries.  Multivariate Behavioral Research4711–25.
  * Guadagnoli Cleary (1992) guadagnoliCleary:1992Guadagnoli, E. Cleary, PD.  1992\.  Age-related item nonresponse in surveys of recently discharged patients Age-related item nonresponse in surveys of recently discharged patients.  Journal of Gerontology473P206–P212.
  * Hans (2009) hans:2009Hans, C.  2009\.  Bayesian lasso regression Bayesian lasso regression.  Biometrika964835–845.
  * Hans (20101) blassoHans, C.  20101\.  blasso: MCMC for Bayesian Lasso Regression Model blasso: Mcmc for bayesian lasso regression model [].  http://www.stat.osu.edu/~hans/ R package version 0.3
  * Hans (20102) hans:2010Hans, C.  20102\.  Model uncertainty and variable selection in Bayesian lasso regression Model uncertainty and variable selection in bayesian lasso regression.  Statistics and Computing202221–229.
  * Hardt . (2012) hardtRainer:2012Hardt, J., Herke, M. Leonhart, R.  2012\.  Auxiliary variables in multiple imputation in regression with missing X: a warning against including too many in small sample research Auxiliary variables in multiple imputation in regression with missing x: a warning against including too many in small sample research.  BMC Medical Research Methodology1211–13.
  * Hastie . (2009) hastieEtAl:2009Hastie, T., Tibshirani, R. Friedman, JH.  2009\.  The elements of statistical learning: data mining, inference, and prediction The elements of statistical learning: data mining, inference, and prediction ( 2).  Springer.
  * Howard . (2015) howardEtAl:2015Howard, WJ., Rhemtulla, M. Little, TD.  2015\.  Using principal components as auxiliary variables in missing data estimation Using principal components as auxiliary variables in missing data estimation.  Multivariate Behavioral Research503285–299.  10.1080/00273171.2014.999267
  * Immerzeel . (2015) immerzeel:2015Immerzeel, T., Coffé, H. Van der Lippe, T.  2015\.  Explaining the gender gap in radical right voting: A cross-national investigation in 12 Western European countries Explaining the gender gap in radical right voting: A cross-national investigation in 12 western european countries.  Comparative European Politics132263–286.
  * James . (2013) jamesEtAl:2013James, G., Witten, D., Hastie, T. Tibshirani, R.  2013\.  An introduction to statistical learning An introduction to statistical learning ( 112).  Springer.
  * Kennickell (1998) kennickell:1998Kennickell, AB.  1998\.  Multiple imputation in the Survey of Consumer Finances Multiple imputation in the survey of consumer finances.  Proceedings of the Section on Survey Research Methods. Proceedings of the section on survey research methods.
  * Köneke (2014) koneke:2014Köneke, V.  2014\.  Trust increases euthanasia acceptance: a multilevel analysis using the European Values Study Trust increases euthanasia acceptance: a multilevel analysis using the european values study.  BMC Medical Ethics15186.
  * Lang . (2018) PcAuxLang, KM., Little, TD. PcAux Development Team.  2018\.  PcAux: Automatically extract auxiliary features for simple, principled missing data analysis Pcaux: Automatically extract auxiliary features for simple, principled missing data analysis [].  https://github.com/PcAux-Package/PcAux R package version 0.0.0.9013
  * Ley Steel (2009) leySteel:2009Ley, E. Steel, MF.  2009\.  On the effect of prior assumptions in Bayesian model averaging with applications to growth regression On the effect of prior assumptions in bayesian model averaging with applications to growth regression.  Journal of applied econometrics244651–674.
  * Little Rubin (2002) littleRubin:2002Little, RJA. Rubin, DB.  2002\.  Statistical analysis with missing data Statistical analysis with missing data (2 ).  Hoboken, NJWiley-Interscience.
  * LWS (2020) LWS:2020LWS.  2020\.  Luxembourg Wealth Study Database. Luxembourg wealth study database.  Luxembourg: LIS.  https://www.lisdatacenter.org/
  * Mainzer . (2021) mainzerEtAl:2021Mainzer, R., Apajee, J., Nguyen, CD., Carlin, JB. Lee, KJ.  2021\.  A comparison of multiple imputation strategies for handling missing data in multi-item scales: Guidance for longitudinal studies A comparison of multiple imputation strategies for handling missing data in multi-item scales: Guidance for longitudinal studies.  Statistics in Medicine40214660–4674.
  * Massey Tourangeau (2013) masseyTourangeau:2013Massey, DS. Tourangeau, R.  2013\.  Where do we go from here? Nonresponse and social measurement Where do we go from here? nonresponse and social measurement.  The ANNALS of the American Academy of Political and Social Science6451222–236.
  * Meng (1994) meng:1994Meng, XL.  1994\.  Multiple-imputation inferences with uncongenial sources of input Multiple-imputation inferences with uncongenial sources of input.  Statistical Science538–558.
  * Meyer . (2015) meyerEtAl:2015Meyer, BD., Mok, WK. Sullivan, JX.  2015\.  Household surveys in crisis Household surveys in crisis.  Journal of Economic Perspectives294199–226.
  * Mustillo (2012) mustillo:2012Mustillo, S.  2012\.  The effects of auxiliary variables on coefficient bias and efficiency in multiple imputation The effects of auxiliary variables on coefficient bias and efficiency in multiple imputation.  Sociological Methods & Research412335–361.
  * Mustillo Kwon (2015) mustilloSoyoung:2015Mustillo, S. Kwon, S.  2015\.  Auxiliary variables in multiple imputation when data are missing not at random Auxiliary variables in multiple imputation when data are missing not at random.  The Journal of Mathematical Sociology39273–91.
  * Muthén . (1987) muthenEtAl:1987Muthén, B., Kaplan, D. Hollis, M.  1987\.  On structural equation modeling with data that are not missing completely at random On structural equation modeling with data that are not missing completely at random.  Psychometrika523431–462.
  * Niolon . (2019) niolonEtAl:2019Niolon, PH., Vivolo-Kantor, AM., Tracy, AJ., Latzman, NE., Little, TD., DeGue, S.others  2019\.  An RCT of dating matters: Effects on teen dating violence and relationship behaviors An rct of dating matters: Effects on teen dating violence and relationship behaviors.  American Journal of Preventive Medicine57113–23.
  * pandas development team (2020) pandas:2020pandas development team, T.  202002\.  pandas-dev/pandas: Pandas. pandas-dev/pandas: Pandas.  Zenodo.  https://doi.org/10.5281/zenodo.3509134 10.5281/zenodo.3509134
  * Park Casella (2008) parkCasella:2008Park, T. Casella, G.  2008\.  The bayesian lasso The bayesian lasso.  Journal of the American Statistical Association103482681–686.
  * R Core Team (2020) R:2020R Core Team.  2020\.  R: A Language and Environment for Statistical Computing R: A language and environment for statistical computing [].  Vienna, Austria.  https://www.R-project.org/
  * Rubin (1986) rubin:1986Rubin, DB.  1986\.  Statistical matching using file concatenation with adjusted weights and multiple imputations Statistical matching using file concatenation with adjusted weights and multiple imputations.  Journal of Business & Economic Statistics4187–94.
  * Rubin (1987) rubin:1987Rubin, DB.  1987\.  Multiple imputation for nonresponse in surveys Multiple imputation for nonresponse in surveys ( 519).  New York, NYJohn Wiley & Sons.
  * Rubin (1996) rubin:1996Rubin, DB.  1996\.  Multiple imputation after 18+ years Multiple imputation after 18+ years.  Journal of the American Statistical Association91434473–489.
  * Rubin . (1995) rubinEtAl:1995Rubin, DB., Stern, HS. Vehovar, V.  1995\.  Handling “don’t know” survey responses: the case of the Slovenian plebiscite Handling “don’t know” survey responses: the case of the slovenian plebiscite.  Journal of the American Statistical Association90431822–828.
  * Schafer (1997) schafer:1997Schafer, JL.  1997\.  Analysis of incomplete multivariate data Analysis of incomplete multivariate data ( 72).  Boca Raton, FLChapman & Hall/CRC.
  * Scott Berger (2010) scottBerger:2010Scott, JG. Berger, JO.  2010\.  Bayes and empirical-Bayes multiplicity adjustment in the variable-selection problem Bayes and empirical-bayes multiplicity adjustment in the variable-selection problem.  The Annals of Statistics2587–2619.
  * Shah (2018) CALIBERrfimputeShah, AD.  2018\.  CALIBERrfimpute: Imputation in MICE using Random Forest Caliberrfimpute: Imputation in mice using random forest [].  R package version 1.0-1
  * Shah . (2014) shahEtAl:2014Shah, AD., Bartlett, JW., Carpenter, J., Nicholas, O. Hemingway, H.  2014\.  Comparison of random forest and parametric imputation models for imputing missing data using MICE: a CALIBER study Comparison of random forest and parametric imputation models for imputing missing data using mice: a caliber study.  American Journal of Epidemiology1796764–774.
  * Song Belin (2004) songBelin:2004Song, J. Belin, TR.  2004\.  Imputation for incomplete high-dimensional multivariate normal data using a common factor model Imputation for incomplete high-dimensional multivariate normal data using a common factor model.  Statistics in Medicine23182827–2843.  10.1002/sim.1867
  * Tharp (2012) tharp:2012Tharp, AT.  2012\.  Dating matters™: The next generation of teen dating violence prevention Dating matters™: The next generation of teen dating violence prevention.  Prevention Science134398–401.
  * Tibshirani (1996) tibshirani:1996Tibshirani, R.  1996\.  Regression shrinkage and selection via the lasso Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Methodological)581267–288.
  * Van Buuren (2010) vanBuuren:2010Van Buuren, S.  2010\.  Item imputation without specifying scale structure Item imputation without specifying scale structure.  Methodology.
  * Van Buuren (2018) vanBuuren:2018Van Buuren, S.  2018\.  Flexible imputation of missing data Flexible imputation of missing data.  CRC press.
  * Van Buuren . (1999) vanBuurenEtAl:1999Van Buuren, S., Boshuizen, HC. Knook, DL.  1999\.  Multiple imputation of missing blood pressure covariates in survival analysis Multiple imputation of missing blood pressure covariates in survival analysis.  Statistics in Medicine186681–694.
  * Vivolo-Kantor . (2021) vivoloEtAl:2021Vivolo-Kantor, AM., Niolon, PH., Estefan, LF., Le, VD., Tracy, AJ., Latzman, NE.Tharp, AT.  2021\.  Middle school effects of the Dating Matters® comprehensive teen dating violence prevention model on physical violence, bullying, and cyberbullying: A cluster-randomized controlled trial Middle school effects of the dating matters® comprehensive teen dating violence prevention model on physical violence, bullying, and cyberbullying: A cluster-randomized controlled trial.  Prevention Science222151–161.
  * von Hippel Lynch (2013) vonHippel:2013von Hippel, P. Lynch, J.  2013\.  Efficiency gains from using auxiliary variables in imputation Efficiency gains from using auxiliary variables in imputation.  arXiv preprint arXiv:1311.5249.
  * White . (2011) whiteEtAl:2011White, IR., Royston, P. Wood, AM.  2011\.  Multiple imputation using chained equations: Issues and guidance for practice Multiple imputation using chained equations: Issues and guidance for practice.  Statistics in Medicine304377–399.
  * Williams Brick (2018) williamsBrick:2018Williams, D. Brick, JM.  2018\.  Trends in US face-to-face household survey nonresponse and level of effort Trends in us face-to-face household survey nonresponse and level of effort.  Journal of Survey Statistics and Methodology.
  * Zhao Long (2016) zhaoLong:2016Zhao, Y. Long, Q.  2016\.  Multiple imputation in the presence of high-dimensional data Multiple imputation in the presence of high-dimensional data.  Statistical Methods in Medical Research2552021–2035.  10.1177/0962280213511027
  * Zou (2006) zou:2006Zou, H.  2006\.  The adaptive lasso and its oracle properties The adaptive lasso and its oracle properties.  Journal of the American statistical association1014761418–1429.
  * Zou Hastie (2005) zouHastie:2005Zou, H. Hastie, T.  2005\.  Regularization and variable selection via the elastic net Regularization and variable selection via the elastic net.  Journal of the royal statistical society: series B (statistical methodology)672301–320.