Edoardo Costantini Kyle M. Lang Klaas Sijtsma Tilburg University, Department
of Methodology and Statistics Utrecht University, Department of Methodology
and Statistics Tilburg University, Department of Methodology and Statistics
Corresponding author’s email address:
<EMAIL_ADDRESS>

# Supervised dimensionality reduction for multiple imputation by chained
equations

###### Abstract

Multivariate imputation by chained equations (MICE) is one of the most popular
approaches to address missing values in a data set. This approach requires
specifying a univariate imputation model for every variable under imputation.
The specification of which predictors should be included in these univariate
imputation models can be a daunting task. Principal component analysis (PCA)
can simplify this process by replacing all of the potential imputation model
predictors with a few components summarizing their variance. In this article,
we extend the use of PCA with MICE to include a supervised aspect whereby
information from the variables under imputation is incorporated into the
principal component estimation. We conducted an extensive simulation study to
assess the statistical properties of MICE with different versions of
supervised dimensionality reduction and we compared them with the use of
classical unsupervised PCA as a simpler dimensionality reduction technique.

##  1 Introduction

Multiple Imputation (MI) is a state-of-the-art missing data treatment in
today’s data analysis (Schafer Graham, 2002; Van Buuren, 2018, p. 30).
Multiple imputations is often implemented through the multivariate imputation
by chained equations approach (MICE, Van Buuren Oudshoorn, 2000), which can
accommodate a wide range of data measurement levels. This flexibility comes
from the possibility of modeling the multivariate joint density of the
variables with missing values through a collection of conditional densities
for every variable under imputation.

MICE requires specifying a different univariate imputation model for each
variable under imputation, which entails deciding on the imputation model form
and which predictors to use. The first decision is usually guided by the
measurement level of the variables under imputation. For example, continuous
variables can be imputed using a linear regression model, while binary
variables can be imputed using logistic regression. The second decision
concerns which and how many predictors to include in the imputation models and
therefore it is more difficult. The general recommendation has been to follow
an inclusive strategy (Collins ., 2001), meaning that as many predictors as
possible should be included in the imputation models. Using as much
information as possible from the data leads to multiple imputations that have
minimal bias and maximal efficiency (Meng, 1994; Collins ., 2001).
Furthermore, including more predictors in the imputation models makes the
missing at random assumption (MAR) more plausible (Collins ., 2001, p. 339).
Finally, if the imputation model omits variables that are part of the analysis
model fitted to the data after imputation, the parameter estimates might be
biased (Enders, 2010, p. 229) and estimated confidence intervals might be too
wide (Little Rubin, 2002, p. 218). As a result, including more predictors in
the imputation models increases the range of analysis models that can be
estimated with a given set of imputations (Meng, 1994).

Despite its advantages, the inclusive strategy easily results in singularity
issues (Hastie ., 2009, p. 46) when estimating imputation models.
Consequently, researchers performing MI often face difficult choices on how
many and which variables to use as predictors. High-dimensional prediction
models offer an opportunity to specify the imputation models automatically and
to include more predictors than traditionally possible. For example, ridge
regression (Hoerl Kennard, 1970) can estimate regression models with hundreds
of predictors; lasso regression (Tibshirani, 1996) can perform data-driven
variable selection; decision trees (e.g., Breiman, 2001) can consider hundreds
of variables for their splitting rules; and principal component analysis (PCA,
Jolliffe, 2002) can summarize a large set of predictors using a few
independent linear combinations.

All of these modeling strategies have been implemented in combination with
MICE (Zhao Long, 2016; Deng ., 2016; Burgette Reiter, 2010; Doove ., 2014;
Shah ., 2014; Howard ., 2015). Costantini . (2022) compared their performance
in terms of estimation bias and confidence interval coverage when applied to
data with missing values. They found that using PCA to create summaries of the
many possible imputation model predictors performs particularly well. In a
follow-up study, Costantini . (2023) explored different ways of using PCA with
the MICE algorithm and found that updating the principal components (PCs) for
the imputation of every variable at every iteration provided the lowest bias
and the lowest deviation from nominal confidence interval coverage. However,
these results relied heavily on the number of components computed. With their
simulation study, Costantini . (2023) showed that to achieve small bias and
satisfactory coverage, a researcher imputing the data using PCA to aid
imputation model specification should retain at least as many PCs as the
number of latent variables in the data generating model. This is undesirable
as researchers usually do not know the true number of latent variables.

PCA is an unsupervised dimensionality reduction technique that summarizes the
variability of a set of $P$ variables {$\mathbf{x}_{1},\dots,\mathbf{x}_{p}$}
measured on $N$ observations with a set of $Q$ PCs, with $Q<P$. The PCs can be
used in any regression model as a replacement for the original predictors, an
approach known as Principal Component Regression (PCR; Jolliffe, 2002, pp.
168-173). PCR addresses possible multicollinearity issues afflicting the
model. However, the PCs obtained by PCR cannot take into account variables
that are not part of the set {$\mathbf{x}_{1},\dots,\mathbf{x}_{p}$}, a
feature that can result in PCs that are unrelated or only weakly related to
the dependent variable, which by definition is not included in the set of
predictors. Contrary to PCA, supervised dimensionality reduction (SDR)
techniques use the outcome variable to guide the computation so that the
resulting PCs are both good representations of the predictor variables and
strongly associated with the dependent variable (e.g., Wold, 1975; De Jong
Kiers, 1992; Bair ., 2006). Using SDR within MICE might relax the need to know
the number of latent variables in the data-generating model described by
Costantini . (2023) for PCA. The purpose of this study is to evaluate how SDR
techniques can improve upon unsupervised PCR as a univariate imputation model
in MICE.

In this study, we considered two questions. First, what are the statistical
properties (bias, coverage, confidence interval width) of parameters estimated
from data treated with the MICE algorithm using different versions of SDR as
the univariate imputation models? Second, can using SDR in MICE relax the PCA
requirement of using at least as many PCs as the number of latent variables in
the data-generating model? We used a Monte Carlo simulation study to explore
the performance of different versions of SDR with the MICE algorithm.

The article is structured as follows. In Section 2, we describe the MICE
algorithm, unsupervised and supervised dimensionality reduction, different
versions of SDR, and we propose uses of SDR as a univariate imputation method
for MICE. In Section 4, we describe the Monte Carlo simulation study. Next, we
discuss the main findings (Section 5), our ideas for future research
directions (Section 6), and we provide concluding remarks (Section 7).

##  2 Imputation methods and algorithms

We use the following notation. Indices and scalars are denoted by lowercase
and uppercase letters. For example, $i$ is an index enumerating iterations out
of $I$ total iterations $(i\in\\{1,\dots,I\\})$. Vectors are written in bold
lowercase while matrices are denoted by bold uppercase letters. The
superscript ′ defines the transpose of a matrix. We use the subscript $obs$
and $mis$ to refer to the observed and missing elements in a vector or matrix.

###  2.1 Multivariate imputation by chained equations

Consider data set $\mathbf{Z}$ with $N$ rows and $P$ columns
$\mathbf{z}_{1},\dots,\mathbf{z}_{P}$. We assume that $\mathbf{Z}$ is a random
draw from a multivariate distribution $f(\mathbf{Z}|\mathbf{\theta})$, where
$\mathbf{\theta}$ is a vector of unknown parameters that completely specifies
its multivariate distribution. Let the first $J$ columns of $\mathbf{Z}$ have
missing values. MICE is an iterative algorithm for imputing multivariate
missing data on a variable-by-variable basis. It obtains multiple imputations
for the missing values by drawing from the variable-specific conditional
distributions of the form:

$f(\mathbf{z}_{j}|\mathbf{Z}_{-j},\mathbf{\theta}_{j}),$ (1)

where $\mathbf{z}_{j}$ is a partially observed variable, $\mathbf{Z}_{-j}$ is
the collection of variables in $\mathbf{Z}$ excluding $\mathbf{z}_{j}$, and
$\mathbf{\theta}_{j}$ is a vector of model parameters. The parameters
$\mathbf{\theta}_{j}$ are specific to the respective conditional distributions
and might not determine the unique _true_ joint distribution
$f(\mathbf{Z}|\mathbf{\theta})$. We refer to the conditional distributions in
Equation 1 as (univariate) imputation models.

The MICE algorithm starts by replacing the missing values in each
$\mathbf{z}_{j}$ with initial guesses (e.g., random draws from the observed
values). At iteration $i$, the MICE algorithm imputes successively variables
$\mathbf{z}_{1}$ to $\mathbf{z}_{J}$ by taking draws from the following
distributions:

$\displaystyle\mathbf{\theta}_{j}^{(i)}$ $\displaystyle\sim
f(\mathbf{\theta}_{j}|\mathbf{Z}_{j,\text{obs}},\mathbf{Z}_{-j}^{(i)}),$ (2)
$\displaystyle\mathbf{z}^{(i)}_{j,\text{mis}}$ $\displaystyle\sim
f(\mathbf{z}_{j,\text{mis}}|\mathbf{Z}_{-j}^{(i)},\mathbf{\theta}_{j}^{(i)})$
(3)

Equation 2 is the fully conditional posterior distribution defined by the
product of an uninformative prior distribution for $\mathbf{\theta}_{j}$ and
the likelihood of observing $\mathbf{z}_{j,\text{obs}}$ under the imputation
model for $\mathbf{z}_{j}$. Equation 3 is the posterior predictive
distribution from which updates of the imputations are drawn. In both
equations, $\mathbf{Z}_{-j}^{(i)}$ is
$(\mathbf{z}_{1}^{(i)},\dots,\mathbf{z}_{j-1}^{(i)},\mathbf{z}_{j+1}^{(i-1)},\dots,\mathbf{z}_{J}^{(i-1)},\mathbf{z}_{J+1},\dots,\mathbf{z}_{P})$,
meaning that at all times the most recently imputed values of all variables
are used to impute other variables.

After repeating the sampling steps described by Equation 2 and 3 for every
variable under imputation, the algorithm moves to the next iteration and
repeats the same sampling steps for all variables under imputation. The
convergence of the algorithm is usually assessed by plotting the trends of the
average imputations across iterations for different starting values. After
convergence, the imputations are assumed to be samples from the target
multivariate distribution. With this process, one can generate as many imputed
data sets as desired. Finally, the analysis model used to answer a substantive
researcher question is estimated on each imputed data set, and the parameter
estimates are pooled using Rubin’s rules (Rubin, 1987).

For small values of $P$, the researcher imputing the data can use all of the
columns in $\mathbf{Z}_{-j}$ as predictors in the univariate imputation model
for $\mathbf{z}_{j}$. As $P$ grows larger, the imputer needs to decide which
predictors to include and which to leave out, a task that can require a
considerable amount of expertise in both statistical modeling techniques and
the field of the substantive research question. By summarizing the information
in all of the possible predictors with a few linear combinations of the
columns of $\mathbf{Z}_{-j}$, PCA and other dimensionality reduction
techniques provide an accessible, data-driven way of specifying the imputation
models.

###  2.2 Principal component analysis

PCA is a dimensionality reduction technique that finds a low-dimensional
representation of the variables contained in an $N\times P$ data matrix
$\mathbf{X}$ with minimal loss of information. It does so by finding a
$P\times Q$ matrix of weights $\mathbf{W}$ that defines $Q$ independent linear
combinations of the columns of $\mathbf{X}$111We follow the common practice of
assuming that the columns of $\mathbf{X}$ are mean-centered and scaled to have
a variance of 1. with maximum variance (with $Q\leq P$). These weights project
the columns of $\mathbf{X}$ onto a lower-dimensional subspace to produce the
$N\times Q$ matrix $\mathbf{T}=\mathbf{X}\mathbf{W}$ that summarizes the
information in $\mathbf{X}$.

The $Q$ columns of $\mathbf{T}$ are called the principal components (PCs) of
$\mathbf{X}$. The first PC of $\mathbf{X}$ is the linear combination of the
columns of $\mathbf{X}$ with the largest variance:

$\mathbf{t}_{1}=\mathbf{x}_{1}w_{11}+\mathbf{x}_{2}w_{12}+\dots+\mathbf{x}_{P}w_{1P}=\mathbf{X}\mathbf{w}_{1},$
(4)

with $\mathbf{w}_{1}$ being the $P\times 1$ vector of weights that comprises
the first column of $\mathbf{W}$. The second principal component
($\mathbf{t}_{2}$) is defined by finding the vector of weights
$\mathbf{w}_{2}$ giving the linear combination of
$\mathbf{x}_{1},\dots,\mathbf{x}_{P}$ with maximal variance out of all the
linear combinations that are uncorrelated with $\mathbf{t}_{1}$. Every
subsequent column of $\mathbf{T}$ can be understood in the same way: for
example, $\mathbf{t}_{3}$ is the linear combinations of
$\mathbf{x}_{1},\dots,\mathbf{x}_{P}$ that has maximal variance out of all the
linear combinations that are uncorrelated with $\mathbf{t}_{1}$ and
$\mathbf{t}_{2}$. As a result, all PCs are uncorrelated by definition and
every subsequent PC has a lower variance than the preceding one.

We can also think of PCA as the process of projecting the original data from a
set of oblique vectors in a $P$-dimensional space to a set of orthogonal
vectors in a Q-dimensional subspace. The weight vectors
$\mathbf{w}_{1},\dots,\mathbf{w}_{Q}$ define the directions in which the $N$
observations in $\mathbf{x}_{1},\dots,\mathbf{x}_{P}$ are projected. The
projected values are the principal component scores $\mathbf{T}$.

Estimating PCA is the task of minimizing the criterion:

$\displaystyle\begin{split}(\mathbf{W},\mathbf{P})&=\underset{\mathbf{W},\mathbf{P}}{\operatorname{argmin}}\;\lVert\mathbf{X}-\mathbf{XWP}^{\prime}\rVert^{2}\\\
&=\underset{\mathbf{W},\mathbf{P}}{\operatorname{argmin}}\;\lVert\mathbf{X}-\mathbf{TP}^{\prime}\rVert^{2}\end{split}$
(5)

subject to the constraint $\mathbf{P}^{\prime}\mathbf{P}=\mathbf{I}$, where
$\mathbf{P}^{\prime}$ provides the weights for estimating $\mathbf{X}$ from
$\mathbf{T}$. In other words, PCA finds the matrices $\mathbf{W}$ and
$\mathbf{P}$ that minimize the reconstruction error
$\mathbf{E}_{\mathbf{X}}=\mathbf{X}-\mathbf{XWP}^{\prime}$.

####  2.2.1 Principal component regression

PCR replaces the $P$ predictors of a regression model with $Q$ PCs extracted
from those predictors. Given an outcome variable $\mathbf{y}$ and a set of $P$
predictors $\mathbf{X}$, consider a standard regression model:

$\mathbf{y}=\mathbf{X}\mathbf{\beta}+\epsilon,$ (6)

where $\mathbf{\beta}$ is a $P\times 1$ vector of regression coefficients, and
$\epsilon$ is a $N\times 1$ vector of independent normally distributed errors.
With PCR we use $Q$ PCs of $\mathbf{X}$ in its place so that Equation 6 can be
rewritten as:

$\mathbf{y}=\mathbf{T}\mathbf{\gamma}+\epsilon,$ (7)

where $\mathbf{\gamma}$ is a $Q\times 1$ vector of regression coefficients.
The lower dimensionality of $\mathbf{T}$ compared to $\mathbf{X}$ and the
independence of its columns allow Equation 7 to address the computational
limitations of Equation 6 when $P$ is large or when the variables in
$\mathbf{X}$ are highly correlated.

The optimal number of components $Q$ is never certain. A crucial feature of
PCA is that the first PC explains the maximum amount of variance in all $P$
predictors, the second PC explains the maximum variance of all $P$ residuals,
and so on. This means that the explained variance decreases as fast as the
data allows as more PCs are retained. As a result, $Q$ is usually taken to be
much smaller than $P$; if $Q=P$, the variance in the $P$ predictors would be
redistributed across $P$ new predictors and this would defy the goal of PCA to
reduce the number $P$ considerably while retaining as much variance as
possible. In practice, when estimating PCR, researchers rely on cross-
validation (e.g., Vervloet ., 2016) to guide this decision.

##  3 Multiple imputation with principal component regression

Costantini . (2023) found that the best way to incorporate PCR into MICE is to
extract PCs at every iteration. When imputing $\mathbf{z}_{j}$ in the $i$th
iteration of MICE, the PCs can be estimated from $\mathbf{Z}^{(i)}_{-j}$ and
used as predictors in the univariate imputation model. Each univariate
imputation model can then be defined as:

$f(\mathbf{z}_{j}|\mathbf{T}^{(i)}_{-j},\mathbf{\theta}_{j}),$ (8)

where $\mathbf{T}^{(i)}_{-j}$ is the matrix storing the PC scores estimated on
$\mathbf{Z}^{(i)}_{-j}$. The steps described in Algorithm 1 are followed to
impute $\mathbf{z}_{j}$ with PCR at every iteration. We refer to this use of
PCR within MICE as MI-PCR. This MI-PCR incorporates uncertainty around the
imputation model parameters using bootstrapping following the same principle
as the ‘imputation under the normal linear model with bootstrap’ algorithm
described by Van Buuren (2018, p. 69).

Algorithm 1 Imputation under the PCR model with bootstrap

1:Define $Q$, the number of PCs to compute.

2:For a given $\mathbf{z}_{j}$ variable under imputation, draw a bootstrap
version $\mathbf{z}^{*}_{j,\text{obs}}$ with replacement from the observed
cases $\mathbf{z}_{j,\text{obs}}$, and store as
$\mathbf{Z}^{*}_{-j,\text{obs}}$ the corresponding rows of
$\mathbf{Z}^{(i)}_{-j}$.

3:Center and scale $\mathbf{Z}^{*}_{-j,\text{obs}}$ and store the result as
$\tilde{\mathbf{Z}}^{*}_{-j,\text{obs}}$

4:Center and scale $\mathbf{Z}^{(i)}_{-j,\text{mis}}$ based on the means and
standard deviations of $\mathbf{Z}^{*}_{-j,\text{obs}}$ and store the result
as $\tilde{\mathbf{Z}}_{-j,\text{mis}}$

5:Center $\mathbf{z}^{*}_{j,\text{obs}}$ on its mean value
$\bar{\mathbf{z}}^{*}_{j,\text{obs}}$ and store it in
$\tilde{\mathbf{z}}^{*}_{j,\text{obs}}$.

6:Estimate $\mathbf{W}$ and $\mathbf{P}$ by the eigendecomposition of the
cross-product matrix of $\tilde{\mathbf{Z}}^{*}_{-j,\text{obs}}$

7:Compute the first $Q$ PCs as
$\mathbf{T}^{(i)}_{-j,\text{obs}}=\tilde{\mathbf{Z}}^{*}_{-j,\text{obs}}\mathbf{W}$

8:Regress the mean-centered $\tilde{\mathbf{z}}^{*}_{j,\text{obs}}$ on
$\mathbf{T}^{(i)}_{-j,\text{obs}}$ and store the regression coefficients
$\mathbf{\beta}$

9:Estimate the residual error variance $\sigma^{2}$ as the ratio between the
residual sum of square ($RSS$) and the degrees of freedom ($df$):
$\displaystyle\sigma^{2}$ $\displaystyle=RSS/df$
$\displaystyle=\frac{\Sigma\,(\tilde{\mathbf{z}}^{*}_{j,\text{obs}}-\mathbf{T}^{(i)}_{-j,\text{obs}}\mathbf{\beta})^{2}}{N-Q}$

10:Obtain the predicted values for $\mathbf{z}_{j,miss}$ by
$\hat{z}_{j,miss}=\tilde{\mathbf{Z}}_{-j,\text{mis}}\mathbf{W}\mathbf{\beta}$

11:Obtain imputations by adding normally distributed errors scaled by
$\sigma^{2}$ to these predictions and by adding
$\bar{\mathbf{z}}^{*}_{j,\text{obs}}$ to center them appropriately.

MI-PCR allows the researcher imputating the data to include all predictors in
the imputation model for every variable, bypassing the difficult model
selection step, while preserving the advantages of an inclusive strategy.
However, the performance of MI-PCR is highly sensitive to the number of PCs
computed. Not using enough PCs to adequately represent the latent structure
(i.e., using fewer PCs than the true number of latent variables) will produce
poor imputations (Costantini ., 2023). Furthermore, there is no guarantee the
number of PCs that optimally represent $\mathbf{Z}^{(i)}_{-j}$ will be good
imputation model predictors as the PCs retained might be summarizing
information that is unrelated to the variable under imputation
$\mathbf{z}_{j}$. Finally, performing PCA on large datasets involves demanding
matrix operations. MI-PCR requires repeating these intensive manipulations for
every variable under imputation and every iteration of the MICE algorithm.

###  3.1 Supervised dimensionality reduction and MICE

SDR techniques represent an alternative to PCA that could obviate some of the
limitations of MI-PCR outlined above. In particular, the _superivsion_ of SDR
methods should help computing PCs that are better predictors of the variables
under imputation than the ones produced by PCA, and, as a corollary, it could
also allow to retain fewer PCs than the number of latent variables in the
data-generating model. In what follows we describe three alternative
approaches to finding linear combinations of the predictors that do a good job
of both summarizing the predictors and predicting the outcome variable. For
each approach, we first describe how it works, and then we describe its
implementation as a univariate imputation method in the MICE algorithm.

####  3.1.1 Supervised principal component regression

Bair . (2006) proposed computing the PCs only on the subset of variables that
are associated with the dependent variable. Their approach is straightforward:

  1. 1.

Regress $\mathbf{y}$ onto each column of $\mathbf{X}$ via $P$ separate simple
linear regressions. Because the data are standardized, the regression
coefficients of these simple linear regression are equivalent to correlation
coefficients. The strength of the association is what matters in the
predictive task, so we consider only the absolute value of the correlation and
refer to it as $\hat{\rho}$.

  2. 2.

Define the subset $\mathbf{X}_{s}\in\mathbf{X}$ by discarding all variables
whose correlation $\hat{\rho}$ is less than a selected threshold $\rho_{s}$.

  3. 3.

Use $\mathbf{X}_{s}$ to estimate the PCs.

  4. 4.

Use these PCs as independent variables in the PCR model.

A key aspect of the method is that both the number of PCs and the threshold
value $\rho_{s}$ can be determined by cross-validation. We refer to this
approach as supervised principal component regression (SPCR).

In SPCR, the component weights are estimated by minimizing the same criterion
as in Equation 5, but only a subset of relevant variables from $\mathbf{X}$ is
used for the computation. By doing so, SPCR effectively sets to 0 component
weights for variables that are not relevant predictors of $\mathbf{y}$. As a
result, SPCR produces PCs that are better predictors of $\mathbf{y}$ and
improves the predictive performance of PCR. We refer to the approach of
excluding variables that are uninteresting for the prediction of the dependent
variable as _discrete_ supervision.

A similar approach, known as Sparse PCA (Zou ., 2006), reduces the number of
variables explicitly used in the PC computation by combining a lasso penalty
with the PCA optimization criterion. Similarly to SPCR, Sparse PCA sets
certain loadings to 0, but it does so to increase the interpretability of the
resulting PCs, not to improve their predictive performance. This key
difference makes SPCR a more suitable tool than Sparse PCA for aiding
automatic imputation model specification. Therefore, we considered SPCR and
not Sparse PCA as a means to reduce the dimensionality of the imputation
models.

In the context of imputation, SPCR can be used as a univariate imputation
model in a similar way to PCR. For each partially observed $\mathbf{z}_{j}$,
with $j\in\\{1,\dots,J\\}$, the imputation model can be defined as:

$f(\mathbf{z}_{j}|\mathbf{T}^{(i)}_{s},\mathbf{\theta}_{j}),$ (9)

where $\mathbf{T}^{(i)}_{s}$ is the matrix of PCs computed on
$\mathbf{Z}^{(i)}_{s}$, the subset of variables with $\hat{\rho}>\rho_{s}$, at
the $i$th iteration of the MICE algorithm. The steps described in Algorithm 2
are followed to impute $\mathbf{z}_{j}$ at every iteration. We refer to this
use of SPCR within MICE as MI-SPCR.

Algorithm 2 Imputation under the SPCR model with bootstrap

1:Define $\mathbf{\rho}$, the vector of possible threshold values $\rho_{h}$
(with $h\in\\{1,\dots,H\\}$). Define $Q$, the number of PCs to retain.

2:For a given $\mathbf{z}_{j}$ variable under imputation, draw a bootstrap
version $\mathbf{z}^{*}_{j,\text{obs}}$ with replacement from the observed
cases $\mathbf{z}_{j,\text{obs}}$, and store as
$\mathbf{Z}^{*}_{-j,\text{obs}}$ the corresponding values on
$\mathbf{Z}^{(i)}_{-j,\text{obs}}$.

3:Compute the absolute correlation between $\mathbf{z}^{*}_{j,\text{obs}}$ and
every potential predictor in $\mathbf{Z}^{*}_{-j,\text{obs}}$.

4:For every $\rho_{h}$, create a set of predictors with absolute correlation
higher than $\rho_{h}$.

5:Use _K_ -fold cross-validation to select the value of $\rho_{h}$ and the
associated set of predictors that return the PCR model with the smallest
prediction error. Define $\rho_{s}$ as the selected value.

6:Drop from $\mathbf{Z}^{*}_{-j,\text{obs}}$ and
$\mathbf{Z}^{(i)}_{-j,\text{mis}}$ all variables with an absolute correlation
smaller than $\rho_{s}$ and create $\mathbf{Z}^{*}_{s,obs}$ and
$\mathbf{Z}_{s,mis}$.

7:Center and scale $\mathbf{Z}^{*}_{s,obs}$ and store the result as
$\tilde{\mathbf{Z}}^{*}_{s,obs}$

8:Center and scale $\mathbf{Z}_{s,mis}$ based on based on the means and
standard deviations of $\mathbf{Z}^{*}_{s,obs}$ and store the result as
$\tilde{\mathbf{Z}}_{s,mis}$

9:Center $\mathbf{z}^{*}_{j,\text{obs}}$ on its mean value
$\bar{\mathbf{z}}^{*}_{j,\text{obs}}$ and store it in
$\tilde{\mathbf{z}}^{*}_{j,\text{obs}}$.

10:Estimate $\mathbf{W}$ and $\mathbf{P}$ by the eigendecomposition of the
cross-product matrix of $\tilde{\mathbf{Z}}^{*}_{s,obs}$

11:Compute the first $Q$ PCs as
$\mathbf{T}^{(i)}_{s,obs}=\tilde{\mathbf{Z}}^{*}_{s,obs}\mathbf{W}$

12:Regress the mean centered $\tilde{\mathbf{z}}^{*}_{j,\text{obs}}$ on
$\mathbf{T}^{(i)}_{s,obs}$ and store the regression coefficients
$\mathbf{\beta}$.

13:Estimate the residual error variance $\sigma^{2}$ as the ratio between the
residual sum of square ($RSS$) and the degrees of freedom ($df$):
$\displaystyle\sigma^{2}$ $\displaystyle=RSS/df$
$\displaystyle=\frac{\Sigma\,(\tilde{\mathbf{z}}^{*}_{j,\text{obs}}-\mathbf{T}^{(i)}_{s,obs}\mathbf{\beta})^{2}}{N-Q}$

14:Obtain the predicted values for $\mathbf{z}_{j,miss}$ by
$\hat{z}_{j,miss}=\tilde{\mathbf{Z}}_{s,mis}\mathbf{W}\mathbf{\beta}$

15:Obtain imputations by adding normally distributed errors scaled by
$\sigma^{2}$ to these predictions and by adding
$\bar{\mathbf{z}}^{*}_{j,\text{obs}}$ to center them appropriately.

In our implementation of MI-SPCR, _K_ -fold cross-validation (Hastie ., 2009,
pp. 241–245) is used to select $\rho_{s}$, from a user-defined vector of
possible values. For every threshold value in the interval $[0,1]$, all
predictors of $\mathbf{z}^{*}_{j,\text{obs}}$ in
$\mathbf{Z}^{*}_{-j,\text{obs}}$ with a correlation larger than the threshold
form an active set of predictors. Then, $Q$ PCs are extracted from each active
set and used to predict $\mathbf{z}^{*}_{j,\text{obs}}$ in a _K_ -fold cross-
validation procedure. The active set giving the lowest cross-validated
prediction error is kept. As with MI-PCR, the number of components $Q$ is
considered fixed, but it can be selected by the same cross-validation
procedure. Note that for a given number of components, only certain threshold
values are allowed. We can compute $Q$ components only if the data have at
least $Q$ columns. Therefore, the more components we want to estimate, the
less restrictive $\rho_{s}$ can be. If we ask for as many components as there
are columns in the data, then $\rho_{s}$ must be large enough to keep all
columns of the data, making MI-SPCR equivalent to MI-PCR.

####  3.1.2 Principal covariates regression

Principal covariates regression (PCovR, De Jong Kiers, 1992) is an SDR
approach that modifies the optimization criteria behind PCA to include
information from the outcome variable in the optimization problem. PCovR looks
for a low-dimensional representation of $\mathbf{X}$ that accounts for the
maximum amount of variation in both $\mathbf{X}$ and $\mathbf{y}$. To
understand how PCovR differs from PCR consider the following decomposition of
the data:

$\displaystyle\mathbf{X}$
$\displaystyle=\mathbf{T}\mathbf{P}^{\prime}_{\mathbf{X}}+\mathbf{E}_{\mathbf{X}}$
(10) $\displaystyle\mathbf{y}$
$\displaystyle=\mathbf{T}\mathbf{P}^{\prime}_{\mathbf{y}}+\mathbf{e}_{\mathbf{y}}$
(11) $\displaystyle\mathbf{T}$ $\displaystyle=\mathbf{X}\mathbf{W}$ (12)

where $\mathbf{T}$ and $\mathbf{W}$ are defined as in 5,
$\mathbf{P}_{\mathbf{X}}$ is $\mathbf{P}$ from 5, and
$\mathbf{P}_{\mathbf{y}}$ is the $Q\times 1$ vector of weights relating
$\mathbf{y}$ to the component scores in $\mathbf{T}$.
$\mathbf{E}_{\mathbf{X}}$ and $\mathbf{e}_{\mathbf{y}}$ are reconstruction
errors. They represent the information lost by using $\mathbf{T}$ as a summary
of $\mathbf{X}$ and the errors in the linear regression model, respectively.
PCovR can be formulated as the task of minimizing a weighted combination of
both $\mathbf{E}_{\mathbf{X}}$ and $\mathbf{e}_{\mathbf{y}}$:

$(\mathbf{W},\mathbf{P}_{\mathbf{X}},\mathbf{P}_{\mathbf{y}})=\underset{\mathbf{W},\mathbf{P}_{\mathbf{X}},\mathbf{P}_{\mathbf{y}}}{\operatorname{argmin}}\;\alpha\,\lVert(\mathbf{X}-\mathbf{XWP}_{\mathbf{X}}^{\prime})\rVert^{2}+(1-\alpha)\,\lVert(y-\mathbf{XWP}_{\mathbf{y}}^{\prime})\rVert^{2}$
(13)

subject to the constraint
$\mathbf{W}^{\prime}\mathbf{X}^{\prime}\mathbf{XW}=\mathbf{T}^{\prime}\mathbf{T}=\mathbf{I}$.

The $\alpha$ parameter defines which reconstruction error is being
prioritized. When $\alpha=1$, the emphasis is exclusively placed on
reconstructing $\mathbf{X}$, casting PCR as a special case of PCovR. When
$\alpha=0.5$, the importance of $\mathbf{X}$ and $\mathbf{y}$ is equally
weighted, a case that resembles partial least square regression (PLSR), which
we discuss in subsection 3.1.3. In practice, the value of $\alpha$ can be
found by cross-validation or according to a sequential procedure based on
maximum likelihood principles (Vervloet ., 2013). In particular,

$\alpha_{ML}=\frac{\lVert\mathbf{X}\lVert^{2}}{\lVert\mathbf{X}\lVert^{2}+\lVert
y\lVert^{2}\frac{\hat{\sigma}_{\mathbf{E}_{\mathbf{X}}}^{2}}{\hat{\sigma}_{e_{\mathbf{y}}}^{2}}}$
(14)

where $\hat{\sigma}_{\mathbf{E}_{\mathbf{X}}}^{2}$ can be obtained as the
unexplained variance by components computed according to classical PCA on
$\mathbf{X}$ and $\hat{\sigma}_{e_{\mathbf{y}}}^{2}$ can be estimated as the
unexplained variance by the linear model regressing $\mathbf{y}$ on
$\mathbf{X}$. Note that, for the same data set, the more components are
retained, the smaller $\hat{\sigma}_{\mathbf{E}_{\mathbf{X}}}$ is, the higher
$\alpha_{ML}$ is, and the closer PCovR becomes to PCR. Retaining the same
number of components as the number of variables in $\mathbf{X}$ results in
$\hat{\sigma}_{\mathbf{E}_{\mathbf{X}}}=0$ and $\alpha=1$, which casts PCR as
a special case of PCovR.

Compared to PCR, PCovR allows estimating PCs that not only represent well the
predictor variables but also predict well the dependent variable. Compared to
SPCR, the PCs computed with PCovR are always linear combinations of _all_
variables in $\mathbf{X}$. PCovR can downweigh irrelevant variables for the
prediction of $\mathbf{y}$, but it will never exclude them entirely. We refer
to the PCovR approach to supervision as _continuous_ , as opposed to the
_discrete_ supervision of SPCR.

When applied to the MICE algorithm, we can use PCovR as a univariate
imputation model in a similar way to how we can use PCR and SPCR. At every
iteration of the MICE algorithm, the steps described in Algorithm 3 are
followed to impute $\mathbf{z}_{j}$. We refer to this use of PCovR within MICE
as MI-PCovR.

Algorithm 3 Imputation under the PCovR model with bootstrap

1:Define $Q$, the number of PCs to compute.

2:For a given $\mathbf{z}_{j}$ variable under imputation, draw a bootstrap
version $\mathbf{z}^{*}_{j,\text{obs}}$ with replacement from the observed
cases $\mathbf{z}_{j,\text{obs}}$, and store as
$\mathbf{Z}^{*}_{-j,\text{obs}}$ the corresponding values on
$\mathbf{Z}^{(i)}_{-j}$.

3:Center and scale $\mathbf{Z}^{*}_{-j,\text{obs}}$ and store the result as
$\tilde{\mathbf{Z}}^{*}_{-j,\text{obs}}$.

4:Center and scale $\mathbf{Z}^{(i)}_{-j,\text{mis}}$ based on based on the
means and standard deviations of $\mathbf{Z}^{*}_{-j,\text{obs}}$ and store
the result as $\tilde{\mathbf{Z}}_{-j,\text{mis}}$.

5:Center $\mathbf{z}^{*}_{j,\text{obs}}$ on its mean value
$\bar{\mathbf{z}}^{*}_{j,\text{obs}}$ and store it in
$\tilde{\mathbf{z}}^{*}_{j,\text{obs}}$.

6:Compute the value of $\alpha$ based on 14.

7:Compute $Q$ PCs by estimating the PCovR $\mathbf{W}$, $\mathbf{P}_{-j}$, and
$\mathbf{P}_{j}$ based on $\tilde{\mathbf{z}}^{*}_{j,\text{obs}}$ and
$\tilde{\mathbf{Z}}^{*}_{-j,\text{obs}}$. Note that $\mathbf{P}_{-j}$, and
$\mathbf{P}_{j}$ correspond to $\mathbf{P}_{\mathbf{X}}$, and
$\mathbf{P}_{\mathbf{y}}$ of equation (13).

8:Estimate the residual error variance $\sigma^{2}$ as the ratio between the
residual sum of square ($RSS$) and the degrees of freedom ($df$):
$\displaystyle\sigma^{2}$ $\displaystyle=RSS/df$
$\displaystyle=\frac{\Sigma\,(\tilde{\mathbf{z}}^{*}_{j,\text{obs}}-\tilde{\mathbf{Z}}^{*}_{-j,\text{obs}}\mathbf{W}\mathbf{P}_{j}^{\prime})^{2}}{N-Q}$

9:Obtain the predicted values for $\mathbf{z}_{j,miss}$ by
$\hat{z}_{j,miss}=\tilde{\mathbf{Z}}_{-j,\text{mis}}\mathbf{W}\mathbf{P}_{-j}^{\prime}$

10:Obtain imputations by adding normally distributed errors scaled by
$\sigma^{2}$ to these predictions and by adding
$\bar{\mathbf{z}}^{*}_{j,\text{obs}}$ to center them appropriately.

####  3.1.3 Partial least square regression

Partial least square regression (PLS, Wold, 1975) is a dimensionality
reduction technique that seeks linear combinations (or PLS components) that
account for a large proportion of the variance in the predictors and correlate
strongly with the dependent variable. Like PCR, PLSR finds independent linear
combinations of the predictors in $\mathbf{X}$ that summarize the data well
and uses these linear combinations to predict the dependent variable. Like
PCovR, the weights defining the linear combinations of the predictors are
computed using all the predictor variables and the outcome, resulting in a
_continuous_ supervision: irrelevant variables will be downweighted in the
linear combinations, but they will not be completely ignored. Unlike PCR,
SPCR, and PCovR, PLSR computes one linear combination at a time and stops at
the required number of PLS components $Q$.

PLSR estimates $\mathbf{t}_{1}$, the first PLS component, by:

  1. 1.

Computing the vector of weights $\mathbf{w}_{1}$ with elements
$w_{1,j}=\mathbf{x}_{j}^{\prime}\mathbf{y}$ for $j\in\\{1,\dots,P\\}$, the
inner products between each predictor $\mathbf{x}_{j}$ and the dependent
variable $\mathbf{y}$.

  2. 2.

Deriving the constructed variable
$\mathbf{t}_{1}=\sum_{j=1}^{P}\mathbf{x}_{j}w_{1,j}=\mathbf{X}\mathbf{w}_{1}$.

  3. 3.

Orthogonalizing $\mathbf{x}_{1}$ to $\mathbf{x}_{P}$ with respect to
$\mathbf{t}_{1}$.

The second linear combination ($\mathbf{t}_{2}=\mathbf{X}\mathbf{w}_{2}$) is
then derived by repeating the same procedure but replacing each
$\mathbf{x}_{j}$ with their versions orthogonalized with respect to
$\mathbf{t}_{1}$. In PLSR, the $q$th weight vector ($\mathbf{w}_{q}$)
maximizes the following optimization criterion (Stone Brooks, 1990; Frank
Friedman, 1993):

$\underset{\mathbf{w}_{q}}{\operatorname{argmax}}\;\text{Corr}^{2}(\mathbf{y},\mathbf{X}\mathbf{w}_{q})\text{Var}(\mathbf{X}\mathbf{w}_{q})$
(15)

where $\text{Corr}^{2}(.)$ is the squared correlation of the vectors between
brackets. As with all other methods, the linear combinations derived by the
PLS algorithm are constrained to be mutually orthogonal.

As with SPCR and PCovR, at every iteration of the MICE algorithm, we can use
PLSR to obtain imputations. Algorithm 4 describes the univariate imputation
method based on PLSR222A similar version of PLS as a univariate imputation
method has also been implemented in the R package ‘miceadds’ (Robitzsch Grund,
2022). that we used to impute $\mathbf{z}_{j}$. We refer to this use of PLSR
within MICE as MI-PLSR. In Table 1, we summarize the differences between the
univariate imputation methods used by MI-PCR and the SDR-based approaches we
described.

Algorithm 4 Imputation under the PLSR model with bootstrap

1:Define $Q$, the number of components to compute.

2:For a given $\mathbf{z}_{j}$ variable under imputation, draw a bootstrap
version $\mathbf{z}^{*}_{j,\text{obs}}$ with replacement from the observed
cases $\mathbf{z}_{j,\text{obs}}$, and store as
$\mathbf{Z}^{*}_{-j,\text{obs}}$ the corresponding values on
$\mathbf{Z}^{(i)}_{-j}$.

3:Estimate PLSR with $Q$ components by regressing
$\mathbf{z}^{*}_{j,\text{obs}}$ onto $\mathbf{Z}^{*}_{-j,\text{obs}}$.

4:Estimate the residual error variance $\sigma^{2}$ as the ratio between the
residual sum of square ($RSS$) and the degrees of freedom ($df$).

5:Obtain the predicted values for $\mathbf{z}_{j,miss}$ based on the trained
PLSR model.

6:Obtain imputations by adding noise scaled by $\sigma^{2}$ to these
predictions.

Method |  | Supervision
---
type
| Optimization
---
criterion
| Estimated
---
parameters
| Tuning
---
parameters
MI-PCR | none | $\displaystyle\underset{\mathbf{W},\mathbf{P}}{\operatorname{argmin}}\;\lVert\tilde{\mathbf{Z}}^{*}_{-j,\text{obs}}-\tilde{\mathbf{Z}}^{*}_{-j,\text{obs}}\mathbf{WP}^{\prime}\rVert^{2}$ | $\displaystyle\mathbf{W},\mathbf{P}$ | -
MI-SPCR | discrete | $\displaystyle\underset{\mathbf{W},\mathbf{P}}{\operatorname{argmin}}\;\lVert\tilde{\mathbf{Z}}^{*}_{s,obs}-\tilde{\mathbf{Z}}^{*}_{s,obs}\mathbf{WP}^{\prime}\rVert^{2}$ | $\displaystyle\mathbf{W},\mathbf{P}$ | $\displaystyle\rho_{s}$
MI-PCovR | continuous | $\begin{aligned} \underset{\mathbf{W},\mathbf{P}_{\mathbf{X}},\mathbf{P}_{\mathbf{y}}}{\operatorname{argmin}}&\;\alpha\,\lVert(\tilde{\mathbf{Z}}^{*}_{-j,\text{obs}}-\tilde{\mathbf{Z}}^{*}_{-j,\text{obs}}\mathbf{WP}_{\mathbf{X}}^{\prime})\rVert^{2}+\\\ &+(1-\alpha)\,\lVert(\tilde{\mathbf{z}}^{*}_{j,\text{obs}}-\tilde{\mathbf{Z}}^{*}_{-j,\text{obs}}\mathbf{WP}_{\mathbf{y}}^{\prime})\rVert^{2}\end{aligned}$ | $\displaystyle\mathbf{W},\mathbf{P}_{\mathbf{X}},\mathbf{P}_{\mathbf{y}}$ | $\displaystyle\alpha$
MI-PLSR | continuous | $\displaystyle\underset{\mathbf{w}_{q}}{\operatorname{argmax}}\;\text{Corr}^{2}(\tilde{\mathbf{z}}^{*}_{j,\text{obs}},\tilde{\mathbf{Z}}^{*}_{-j,\text{obs}}\mathbf{w}_{q})\text{Var}(\tilde{\mathbf{Z}}^{*}_{-j,\text{obs}}\mathbf{w}_{q})$ | $\displaystyle\mathbf{w}_{q}$ | -
|  |  |  |

Table 1:  Summary of the differences between the univariate imputation models
used within the MICE algorithm.

##  4 Simulation study

We investigated the relative performance of unsupervised PCR and the
supervised alternatives described above with a Monte Carlo simulation study.
In particular, we investigated the estimation bias, confidence interval width,
and confidence interval coverage of several analysis model parameters obtained
after imputation. We varied the dimensionality of the latent structure in the
data-generating model, the proportion of missing values, the missing data
mechanism, and the number of components used as predictors in the imputation
models. The proportion of missing values and the missing data mechanism both
influence the statistical properties of any missing data treatment, while the
number of latent variables in the data-generating model and the number of
components used in the imputation model affect the extent to which supervision
can improve upon the limitations of MI-PCR. In table 2, we summarize the
experimental factors we varied and the levels we used with each factor.

Experimental factor | Label |  Levels
---|---|---
Number of latent variables | $L$ |  2, 10, 50
Missing data mechanism used | $mech$ |  MCAR, MAR
Proportion of missing values | $pm$ |  0.1, 0.25, 0 50
Missing data treatment | $method$ |  MI-PCR, MI-SPCR, MI-PcoVR, MI-PLS, MI-QP, MI-AM, MI-ALL, CC, FO
Number of components | $nc$ |  0, 1 to 12, 20, 29, 30, 40, 48, 49, 50, 51, 52, 60, 149
|  |
Table 2:  Summary of experimental factors for the simulation study.

###  4.1 Procedure

The simulation study involved four steps:

  1. 1.

Data generation: We generated $R=240$ data sets from a confirmatory factor
analysis model, following the procedure described in Section 4.1.1.

  2. 2.

Missing data generation: We generated missing values on three target items in
each data set, following the procedure described in Section 4.1.2.

  3. 3.

Imputation: We generated $d=5$ multiply imputed versions of each generated
data set using different imputation methods, as described in Section 4.1.3.

  4. 4.

Analysis: We estimated the means, variances, covariances, and correlations of
the three items with missing values on the $d$ imputed data sets, and we
pooled the estimates according to Rubin’s rules (Rubin, 1986, p. 76). We then
assessed each imputation method by computing the bias of different parameter
estimates, and their confidence interval widths and coverages as described in
Section 4.1.4.

####  4.1.1 Data generation

For each of the $R$ replications, we generated a $1000\times P$ data matrix
$\mathbf{Z}$. The sample size should be large enough to generate data sets
that have statistical properties similar to large social science data sets.
Each data set was generated based on the following model:

$\mathbf{Z}=\mathbf{F}\mathbf{\Lambda}^{\prime}+\mathbf{E},$ (16)

where $\mathbf{F}$ is a $1000\times L$ matrix of latent variables scores, $L$
is the number of latent variables, $\mathbf{\Lambda}$ is a $3\times L$ matrix
of factor loadings, where $3$ is the number of items measuring each latent
variable, and $\mathbf{E}$ is a $1000\times P$ matrix of measurement errors,
where $P=3*L$. The dimensionality of the data resembles that of the many large
social surveys that use short scales to measure respondents’ attitudes such as
_political engagement_ and _anti-immigrant attitudes_. For example, consider
the European Values Study (EVS, 2020) which measures a variety of attitudes
with 3, 4, or 5 items.

The factor loading matrix $\mathbf{\Lambda}$ in Equation 16 described a simple
structure (Bollen, 1989, p. 234) where each item loads on exactly one of the
$L$ latent variables. In real data applications, this factor structure is
uncommon but not implausible. For example, a relatively clear simple structure
can be found when analyzing personality inventories (e.g., NEO-PR-I, Costa Jr
., 1991) with the _Neuroticism-Extroversion-Openness_ three-factor model
(McCrae Costa Jr, 1983). A simple structure is also often assumed when
performing exploratory factor analysis because it provides the most
parsimonious explanation (Costa McCrae, 2008, pp. 183-184).

We sampled $\mathbf{F}$ from a multivariate normal distribution with mean
$\mathbf{0}$ and covariance matrix $\mathbf{\Psi}$. The correlation between
the first and second latent variables was fixed at $0.8$, while the
correlation between all other latent variables and the first two was fixed at
$0.1$. Together with factor loadings fixed at $\lambda=0.85$, these choices
resulted in correlations of approximately $0.72$, $0.58$, and $0.07$ between
items measuring the same latent variable, items measuring the first and second
latent variable, and items measuring the first latent variable and the others,
respectively. These values represent plausible, but reasonably high, item-
scale associations and they should mitigate the impact of measurement error on
our findings without resorting to implausibly precise data.

The matrix of measurement errors $\mathbf{E}$ was sampled from a multivariate
normal distribution with mean vector $\mathbf{0}$ and covariance matrix
$\mathbf{\Theta}$. The off-diagonal elements of $\mathbf{\Theta}$ were set to
0 to reflect uncorrelated errors, while the diagonal elements were specified
as $1-\lambda^{2}$ to give the simulated items unit variances. After sampling,
the columns of $\mathbf{Z}$ were rescaled to have approximately a mean of 5
and a variance of 6.5, which are common values for Likert items in social
surveys measured on a 10-point scale (for example in EVS, 2020).

In this data-generating procedure, we considered the number of latent
variables used $(L)$ as an experimental factor with levels $2,10,50$,
resulting in data sets containing 6, 30, and 150 total items. Costantini .
(2023) showed that the number of components used in MI-PCR needs to be at
least as high as the number of latent variables in the data-generating model.
So, we expected MI-PCR to require at least $L$ components to achieve
satisfactory performance. One of this study’s main objectives was to
understand how well supervision can overcome this limitation of MI-PCR. We
generated data according to a confirmatory factor analysis model instead of
generating the data directly based on true principal components to avoid the
misleading results that can occur when using a single model for both data
generation and imputation (see Oberman Vink, 2023, p. 4).

####  4.1.2 Missing data generation

We generated missing data on the three items measuring the first latent
variable ($\mathbf{z}_{1}$, $\mathbf{z}_{2}$, $\mathbf{z}_{3}$). The
proportion of missing values per variable ($pm$) was defined as an
experimental factor taking three levels $pm\in\\{0.1,0.25,0.5\\}$. The missing
data mechanism ($mech$) was a factor with two levels:

  * •

Missing completely at random (MCAR): To test how the methods performed in the
simplest possible missing data mechanism, we generated missing values on each
item based on a missing data indicator ($\delta$) sampled from a binomial
distribution with success probability $pm$. If $\delta=1$, the item score was
set to missing. If $\delta=0$, the item score was set to observed. As a
result, every variable had a proportion of missing values approximately equal
to $pm$.

  * •

Missing at random (MAR): To test how the methods performed in a more realistic
situation, we generated missing values based on a MAR mechanism with the three
items measuring the second latent variable ($\mathbf{z}_{4}$,
$\mathbf{z}_{5}$, $\mathbf{z}_{6}$) used as predictors of missingness. We
sampled $\delta$ from Bernoulli distributions with probabilities defined based
on the following logit model:

$logit(\delta=1)=\beta_{0}+\mathbf{Z}_{(4,5,6)}\mathbf{\beta},$ (17)

where $\beta_{0}$ is an intercept parameter, and $\mathbf{\beta}$ is a vector
of slope parameters. All slopes in $\mathbf{\beta}$ were fixed to 1, while the
value of $\beta_{0}$ was chosen with an optimization algorithm that minimized
the difference between the actual and desired proportion of missing values on
the variable. The pseudo R-squared for the logistic regression of the missing
value indicator on the predictors of missingness was approximately 14%. The
AUC for the logistic regression was approximately 0.74. To create realistic
missing data patterns, the location of missing data was fixed to right for
$\mathbf{z}_{1}$, left for $\mathbf{z}_{2}$, and tails for $\mathbf{z}_{3}$.

####  4.1.3 Imputation

We imputed the missing values using the four dimension reduction-based methods
described above (MI-PCR, MI-SPCR, MI-PCovR, MI-PLSR) as well as three
traditional approaches:

  * •

MI with all the available variables used as predictors in the imputation
models (MI-ALL), which represents the most naive way to define the imputation
model.

  * •

MI with a correlation-based threshold strategy to select the subset of
important predictors (MI-QP). As a pragmatic point of comparison, this method
used the _quickpred_ function from the R package mice (Van Buuren Groothuis-
Oudshoorn, 2011) to select the predictors for the univariate imputation models
via the correlation-based threshold strategy described by Van Buuren . (1999,
pp. 687-–688). To implement this approach, we selected only those predictors
that correlated with the imputation targets (or their associated missingness
indicators) higher than $0.1$.

  * •

MI with the analysis model variables used as sole predictors in the imputation
models (MI-AM). This method produces the simplest possible congenial
imputation model and is interesting due to its popularity in the social
scientific literature (Costantini ., 2022).

As reference points, we also treated the missing values with complete case
analysis (CC) and estimated the analysis model from the original, fully
observed data.

Every imputation-based method used simple random draws from the observed data
as starting values and was run to obtain 5 imputed data sets. Convergence was
achieved after 20 iterations for all methods333Convergence plots are reported
in the interactive results dashboard that we developed to accompany this
article. See Section 4.2 for more details.. All our dimensionality reduction-
based algorithms need the user to define $Q$, the optimal number of
components. In practice, researchers will choose a single value of $Q$ with a
cross-validation procedure, but in this study, we are interested in evaluating
the performance of the four dimension-reduction imputation approaches while
varying the number of retained components. Therefore, we defined $Q$ as one of
our main experimental factors (the number of components, $nc$) taking values
$\\{1,\ldots,12,20,29,30,40,48,49,50,51,52,60,149\\}$. These values were
chosen to both cover the range of possible choices (i.e., $1,\dots,[P-1]$) and
to provide more granularity around the true number of latent variables (i.e.,
3, 10, 50).

Finally, for MI-SPCR, we used the cross-validation procedure described in
Section 3.1.1 to select $\rho_{s}$ from the vector of values
$\\{0.05,0.1,0.15,\ldots,0.95\\}$. For a given number of components, some
threshold values can exclude enough predictors to preclude computing the
required number of components. To avoid this possibility, the cross-validation
algorithm only considered values of $\rho$ that retained enough variables to
compute the required number of components. The weighting parameter ($\alpha$)
for MI-PCovR was selected using the sequential MLE-based estimation procedure
described in Vervloet . (2013). The degrees of freedom for the PLSR imputation
model were computed based on the naive approach described by Krämer Sugiyama
(2011)

####  4.1.4 Analysis and comparison criteria

For a given parameter $\phi$ (e.g., the mean of $\mathbf{z}_{1}$, the
correlation between $\mathbf{z}_{1}$ and $\mathbf{z}_{2}$), we used the
absolute percent relative bias (PRB) to quantify the estimation bias
introduced by the imputation procedure:

$\text{PRB}=\displaystyle\left\lvert\
\frac{\bar{\hat{\phi}}-\phi}{\phi}\right\rvert\times 100$ (18)

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

To measure the statistical efficiency of the imputation methods we computed
the average width of the confidence intervals (CIW).

$\text{CIW}=\frac{\sum_{r=1}^{R}(\widehat{\text{CI}}^{upper}_{r}-\widehat{\text{CI}}^{lower}_{r})}{R},$
(19)

with $\widehat{\text{CI}}^{upper}_{r}$ and $\widehat{\text{CI}}^{lower}_{r}$
being the upper and lower bounds of the estimated confidence interval for the
$r$th replication. Narrower CIWs indicate higher efficiency. However, narrower
CIWs are not preferred if they come at the expense of good confidence interval
coverage (CIC) of the parameter values. CIC is the proportion of confidence
intervals that contain the true value of the parameter, across the $R$ data
samples:

$\text{CIC}=\frac{\sum_{r=1}^{R}I(\phi\in\widehat{\text{CI}}_{r})}{R},$ (20)

where $\widehat{\text{CI}}_{r}$ is the confidence interval of the parameter
estimate $\hat{\phi}_{r}$ in the $r$th replication, and $I(.)$ is the
indicator function that returns 1 if the argument is true and 0 otherwise. CIC
depends on both the bias and the variability of the CIW for a parameter
estimate. In particular, for a given level of bias, a narrower CIW leads to
lower CIC, and, for a given CIW, a larger bias leads to lower CIC. An
imputation method with good coverage should result in CICs greater than or
equal to the nominal rate. For 95% CIs, CIC below 0.9 is usually considered
problematic (e.g., Van Buuren, 2018, p. 52; Collins ., 2001, p. 340) as it
implies inflated Type I error rates. High CIC (e.g., 0.99) implies inflated
Type II error rates.

###  4.2 Results

We report only the results for the correlation between $\mathbf{z}_{1}$ and
$\mathbf{z}_{2}$ for the conditions with _mech_ = MAR and $pm=0.5$ because the
type of parameter and the different levels of these two factors did not impact
the relative performances of the imputation methods. We focused on the
correlation between two items with missing values because this parameter
differentiated the performances of the methods the most. The full set of
results is available via the interactive results dashboard that we developed
to accompany this article (Costantini, 2022). The dashboard can be downloaded
and installed as an R package, and it can be used as an R Shiny app.

In Figures 1, 2, and 3, we report the PRB, CIW, and CIC for the correlation
coefficient between $\mathbf{z}_{1}$ and $\mathbf{z}_{2}$ for different
numbers of latent variables in the data-generating model ($L$) and numbers of
components retained by the methods ($nc$). Across all values of $L$, MI-PCR
resulted in a smaller bias and coverage closer to nominal the more components
were retained. However, MI-PCR required the number of components to be greater
or equal to the number of latent variables used in the data-generating model
to return acceptable bias and coverages. In particular, for $L=2$ and $L=10$,
MI-PCR resulted in acceptable bias ($\text{PRB}<10$) and close to nominal
coverage ($\text{CIC}>0.9$) only when using $nc\geq 2$ and $nc\geq 10$,
respectively. Contrary to expectation, this trend did not persist for all
higher values of $L$ and $nc$. For $L=50$, MI-PCR resulted in high bias
($\text{PRB}>20$) even for $nc=50$. Furthermore, for $L=2$ and $L=10$, MI-PCR
resulted in large deviations from nominal coverage ($\text{CIC}<0.9$) for
$nc=5$ and for $nc\in\\{11,12\\}$, respectively.

Compared to MI-PCR, MI-SPCR performed much better, especially when using just
a few components. MI-SPCR resulted in the lowest bias, smallest confidence
interval width, and closest to nominal coverage when using between 2 and 5
components. For all values of $L$, using 2 components instead of 1, led to a
large reduction in PRB and improvement in CIC. Using 6 or more components had
only a minor negative impact on the performance of the method in the condition
with $L=10$, resulting in a negligible increase in bias. However, for $L=50$,
using 6 or more components did lead to high bias and low coverage. Finally,
the maximum number of components led to algorithmic failures, so there were no
results to report for $nc=29$ and $nc=149$ in the $L=10$ and $L=50$
conditions, respectively.

MI-PCovR resulted in acceptable bias for all values of $L$, for all $nc$
values reported. However, its bias performance was less stable than that of
MI-SPCR across the values of $L$ and $nc$. For $L=10$, MI-PCovR led to smaller
bias when using $nc=2$ instead of $nc=1$, but the bias increased when using
$nc\in\\{3,\dots,9\\}$, only to decrease again for $nc\geq 10$. In the $L=50$
condition, MI-PCovR resulted in decreasing bias for the range
$nc\in\\{1,\dots,9\\}$, and the lowest bias was achieved with $nc=L$, just
after a small increase. Furthermore, the CIC resulted in larger deviations
from nominal coverage than MI-SPCR, resulting in acceptable CIC only with a
few of the many $nc$ values considered.

Compared to MI-SPCR and MI-PCovR, the reduction in bias obtained by MI-PLSR
for higher values of $nc$ was more gradual. For $L=10$ and $50$, using 2
components instead of 1, led to a reduction in bias, but 3 components were
necessary to achieve $\text{PRB}<10$, while both MI-SPCR, and MI-PCovR only
needed 2 to achieve the same result. However, the PRB remained small for
$nc\in\\{6,7,8,9,10\\}$, even when that of MI-SPCR and MI-PCovR increased.
Despite this good bias performance, the CIW and CIC of MI-PLSR fluctuated
between acceptable and not, with only a few values of $nc$ resulting in close-
to-nominal coverage.

To put these results into perspective, we reported the same performance
metrics for three traditional MI approaches, and complete case analysis in
Figure 4. MI-QP and MI-ALL resulted in acceptable bias ($PRB<10$) for all
values of $L$, although larger values of $L$ did result in increased bias and
decreasing coverage for both methods, but especially for MI-ALL. MI-AM and CC
did not result in a higher bias or lower coverages for larger values of $L$,
but both returned relatively high bias and low coverage across all conditions.

Figure 1:  The PRB for the estimated correlation coefficient between the first
two items imputed is reported (Y-axis) as a function of the number of
components ($nc$) used by the PCA-based imputation methods (X-axis). The plot
is divided into a grid where the rows distinguish the results obtained after
imputing the data with the four PCA-based methods and the columns distinguish
the number of latent variables used to generate the data ($L$). All results
plotted in this figure were obtained on data generated with _mech_ = MAR and
$pm=0.5$.  Figure 2:  The CIW for the estimated correlation coefficient
between the first two items imputed is reported (Y-axis) as a function of the
number of components ($nc$) used by the PCA-based imputation methods (X-axis).
The plot is divided into a grid where the rows distinguish the results
obtained after imputing the data with the four PCA-based methods and the
columns distinguish the number of latent variables used to generate the data
($L$). All results plotted in this figure were obtained on data generated with
_mech_ = MAR and $pm=0.5$. The light gray color indicates the parameter
estimate for which the CIW is reported had both acceptable bias
($\text{PRB}<10$) and coverage ($\text{CIC}>0.9$). The dark gray color
indicates the parameter estimate for which the CIW is reported had large bias
($\text{PRB}>10$) or low coverage ($\text{CIC}<0.9$), or both. The black
horizontal lines represent the average CIW obtained on the original fully
observed data.  Figure 3:  The CIC for the estimated correlation coefficient
between the first two items imputed is reported (Y-axis) as a function of the
number of components ($nc$) used by the PCA-based imputation methods (X-axis).
The plot is divided into a grid where the rows distinguish the results
obtained after imputing the data with the four PCA-based methods and the
columns distinguish the number of latent variables used to generate the data
($L$). All results plotted in this figure were obtained on data generated with
_mech_ = MAR and $pm=0.5$. For CIC values below 0.8, we reported the precise
value within the corresponding bar.  Figure 4:  PRB, CIC, and CIW for the
estimated correlation coefficient between the first two items imputed are
reported for the traditional missing data handling methods considered. The
black horizontal lines represent the average CIW obtained for the parameter of
interest when analyzing the original fully observed data. For CIC values below
0.8, we reported the precise value within the corresponding bar.

##  5 Discussion

###  5.1 Supervised dimensionality reduction

Our simulation study outlined some clear advantages of using supervised
dimensionality reduction techniques over standard PCA with MICE. We found that
MI-PCR requires the use of at least as many components as the number of latent
variables in the data-generating model, which is in line with the results
presented by Costantini . (2023). The simulation study presented here also
showed that meeting this requirement is not sufficient to obtain good
imputations for a large number of latent variables. We found that the
performance of MI-PCR does not only depend on knowing the number of latent
variables in the data-generating model, but also on the number itself. On the
contrary, the SDR-based methods retaining just a few components resulted in
small bias and good confidence interval coverage, _independent_ of the number
of latent variables in the data-generating model. Furthermore, when using any
given number of components, except the maximum, using the SDR-based methods
resulted in smaller bias, narrower confidence intervals, and closer to nominal
coverage than MI-PCR. Considering these results, SDR-based MICE seems to be
more appropriate than PCR-based MICE for automatic imputation model
specification.

Among the SDR-based methods, MI-SPCR had the best statistical properties. MI-
SPCR returned smaller bias and better coverage for a wider range of retained
components compared to MI-PCovR. MI-SPCR also achieved a smaller bias than MI-
PLSR when retaining fewer components, and it resulted in consistently closer-
to-nominal coverages. Based on our results it seems that, at least in the
context of the imputation of data with a latent structure, the _discrete_ type
of supervision employed by MI-SPCR should be preferred to the _continuous_
supervision employed by MI-PCovR and MI-PLSR.

###  5.2 Supervision and the number of principal components

Based on our simulation study results, irrespective of which type of
supervised dimensionality reduction is used, the implementation of SDR-based
methods in MICE should aim for computing a small number of components. Despite
this general trend, MI-SPCR and MI-PCovR showed different performances in
relation to the different numbers of components retained.

As described in the results section, the bias obtained by MI-SPCR was the
smallest when using between 2 and 5 components, independently of the number of
latent variables in the data-generating model, and it led to algorithmic
failures for large numbers of components. These results can be explained by
considering how the number of PC retained influences MI-SPCR. In MI-SPCR,
supervision is introduced by pre-screening the columns of possible predictors
set to exclude any predictors that are not correlated strongly enough with the
variable under imputation. By reducing the number of columns in the predictor
set, this supervision reduces the maximum number of components that can be
estimated. The inverse constraint also holds. Fixing the number of components
puts an upper-bound on the number of variables that can be excluded during the
screening process. So, for example, when using 50 components, MI-SPCR must
retain, at least, 50 predictors during the screening step, regardless of how
weakly some of these variables may associate with any given variable under
imputation. Retaining more components forces the threshold values to be
smaller and results in keeping more predictors that are less strongly related
to the dependent variable, and, by doing so, it limits the advantage of using
supervision in the PCA. In conclusion, considering our results and the
relationship between the number of components and supervision, we recommend
retaining between 2 and 5 components when using MI-SPCR.

In the results section, we noted that in the condition with 10 latent
variables, MI-PCovR led to smaller bias when using 2 components instead of 1,
but its bias increased when retaining 3 to 9 components, only to decrease
again when retaining 10 or more components. A similar but less extreme trend
was also detected in the condition with 50 latent variables, a result which we
explore further in the appendix. This fluctuating bias performance can be
explained by considering the relationship between the number of components
retained and the way we computed $\alpha$ in our simulation study. As
described in equation 14, for the same data set, $\alpha_{ML}$ is bigger when
the unexplained variance by the components retained is smaller. The more
components we retained in our simulation study, the closer the value of
$\alpha$ was to 1, and the closer MI-PCovR became to MI-PCR. As a result, MI-
PCovR resulted in smaller bias than MI-PCR when retaining the first
components, as supervision helped to compute leading components that were
important predictors for the imputation task. However, as more components were
retained, the effect of supervision started to diminish, which drove the bias
closer to that of MI-PCR. After this initial increase, the bias achieved by
MI-PCovR dropped when retaining as many components as the number of latent
variables, mirroring the drop in bias presented by MI-PCR for the same number
of components. As a result, the best performances for MI-PCovR could be
achieved by retaining the first components or by retaining a number of
components just above the number of latent variables. Because of this
fluctuating performance, the optimal range of components to consider when
using MI-PCovR is not as clear as for MI-SPCR.

##  6 Limitations and future directions

An important aspect to consider in deciding which version of supervised
dimensionality reduction to use with MICE is how flexible these approaches are
to deviations from normality of the data. This topic was not covered by our
simulation study. MI-SPCR can easily be adapted to impute binary and
categorical variables, and the only complication would be in defining a
suitable threshold parameters. One option would be to estimate the
associations via fit measures derived from simple (multinomial) logistic
regression models. Much research has been dedicated to extending PLSR to
categorical outcomes (e.g., Ding Gentleman, 2005; Chung Keles, 2010), and
these approaches could be used in a similar way to the standard PLS
implementation we used in this study. The development of PCovR for
classification tasks has not received much attention (Park ., , is the only
example of which we are aware), but the same approaches used to fit PLS in the
generalized linear framework should also apply to PCovR. However, the maximum
likelihood estimation of the $\alpha$ parameter can only be done for
continuous dependent variables. To impute categorical variables, PCovR would
require cross-validation to estimate the value of $\alpha$, adding to the
computational intensity of the procedure.

The set of predictors used to compute the components can also include
categorical variables. There are different ways of accommodating these
categorical variables when estimating the components, including the naive
application of traditional PCA Filmer Pritchett (2001) and the PCAMIX
algorithm (Kiers, 1991; Chavent ., 2012, 2014) specifically designed for this
purpose. Which of these approaches is more appropriate for the predictive task
involved in MICE is yet to be tested.

Finally, all of the PCA-based methods considered here (both supervised and
unsupervised) entail a high computational load. For every variable and every
iteration of the MICE algorithm, complex matrix operations need to be
performed to estimate the components. When cross-validating the tuning
parameters, the supervised approaches described in this article can increase.
Future research should explore possible computational shortcuts to perform
supervised dimensionality reduction faster (e.g., Abraham Inouye, 2014; Halko
., 2011).

##  7 Conclusions

Based on the simulation study presented here, it can be concluded that adding
a supervision element to the classical use of PCR as a univariate imputation
method can improve significantly the performance of MI-PCR, especially when
the data contain hundreds of variables. Although there is room to assess the
performance of these imputation methods in more complex data scenarios, MI-
SPCR was particularly effective for the imputation of missing values and seems
to be preferable to MI-PCovR and MI-PLSR.

##  8 Code and data availability

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

## References

  * Abraham Inouye (2014) abrahamInouye:2014Abraham, G. Inouye, M.  2014\.  Fast principal component analysis of large-scale genome-wide data Fast principal component analysis of large-scale genome-wide data.  PloS one94e93766.
  * Bair . (2006) bairEtAl:2006Bair, E., Hastie, T., Paul, D. Tibshirani, R.  2006\.  Prediction by supervised principal components Prediction by supervised principal components.  Journal of the American Statistical Association101473119–137.
  * Bollen (1989) bollen:1989Bollen, KA.  1989\.  Structural Equations with Latent Variables Structural equations with latent variables ( 210).  John Wiley & Sons.
  * Breiman (2001) breiman:2001Breiman, L.  2001\.  Random forests Random forests.  Machine learning4515–32.
  * 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
  * Chavent . (2014) chaventEtAl:2017Chavent, M., Kuentz-Simonet, V., Labenne, A. Saracco, J.  2014\.  Multivariate analysis of mixed data: The R Package PCAmixdata Multivariate analysis of mixed data: The r package pcamixdata.  arXiv preprint arXiv:1411.4911.
  * Chavent . (2012) chaventEtAl:2012Chavent, M., Kuentz-Simonet, V. Saracco, J.  2012\.  Orthogonal rotation in PCAMIX Orthogonal rotation in pcamix.  Advances in Data Analysis and Classification62131–146.
  * Chung Keles (2010) chungKeles:2010Chung, D. Keles, S.  2010\.  Sparse partial least squares classification for high dimensional data Sparse partial least squares classification for high dimensional data.  Statistical applications in genetics and molecular biology91.
  * 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
  * Costa McCrae (2008) costaMcCrae:2008Costa, PTJ. McCrae, RR.  2008\.  The SAGE Handbook of Personality Theory and Assessment: Volume 2 - Personality Measurement and Testing The sage handbook of personality theory and assessment: Volume 2 - personality measurement and testing.  ( 2, 179–198).  Sage Publications, Inc.  https://psycnet.apa.org/doi/10.4135/9781849200479.n9
  * Costa Jr . (1991) costaEtAl:1991Costa Jr, PT., McCrae, RR. Dye, DA.  1991\.  Facet scales for agreeableness and conscientiousness: A revision of the NEO Personality Inventory Facet scales for agreeableness and conscientiousness: A revision of the neo personality inventory.  Personality and individual Differences129887–898.
  * Costantini (2022) costantini:2022bCostantini, E.  202212\.  EdoardoCostantini/plotmispcr: glitter-surprise. Edoardocostantini/plotmispcr: glitter-surprise.  Zenodo.  https://doi.org/10.5281/zenodo.7451802 10.5281/zenodo.7451802
  * Costantini . (2022) costantiniEtAl:2022aCostantini, E., Lang, KM., Reeskens, T. Sijtsma, K.  2022\.  High-dimensional imputation for the social sciences: a comparison of state-of-the-art methods High-dimensional imputation for the social sciences: a comparison of state-of-the-art methods.  arXiv preprint arXiv:2208.13656.
  * Costantini . (2023) costantiniEtAl:2023Costantini, E., Lang, KM., Sijtsma, K. Reeskens, T.  2023\.  Solving the many-variables problem in MICE with principal component regression. Solving the many-variables problem in mice with principal component regression.  Behavior Research Methods.  doi.org/10.3758/s13428-023-02117-1
  * De Jong Kiers (1992) deJongKiers:1992De Jong, S. Kiers, HA.  1992\.  Principal covariates regression: part I. Theory Principal covariates regression: part i. theory.  Chemometrics and Intelligent Laboratory Systems141-3155–164.
  * 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
  * Ding Gentleman (2005) dingGentleman:2005Ding, B. Gentleman, R.  2005\.  Classification using generalized partial least squares Classification using generalized partial least squares.  Journal of Computational and Graphical Statistics142280–298.
  * 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.
  * Enders (2010) enders:2010Enders, CK.  2010\.  Applied missing data analysis Applied missing data analysis.  New York, NYThe Guilford Press.
  * EVS (2020) EVS:2017EVS.  2020\.  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
  * Filmer Pritchett (2001) filmerPritchett:2001Filmer, D. Pritchett, LH.  2001\.  Estimating wealth effects without expenditure data—or tears: an application to educational enrollments in states of India Estimating wealth effects without expenditure data—or tears: an application to educational enrollments in states of india.  Demography381115–132.
  * Frank Friedman (1993) frankFriedman:1993Frank, LE. Friedman, JH.  1993\.  A statistical view of some chemometrics regression tools A statistical view of some chemometrics regression tools.  Technometrics352109–135.
  * Halko . (2011) halkoTropp:2011Halko, N., Martinsson, PG. Tropp, JA.  2011\.  Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions.  SIAM review532217–288.
  * 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.
  * Hoerl Kennard (1970) hoerlKennard:1970Hoerl, AE. Kennard, RW.  1970\.  Ridge regression: Biased estimation for nonorthogonal problems Ridge regression: Biased estimation for nonorthogonal problems.  Technometrics12155–67.
  * 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
  * Jolliffe (2002) jolliffe:2002Jolliffe, IT.  2002\.  Principal component analysis Principal component analysis.  Springer.
  * Kiers (1991) kiers:1991Kiers, HA.  1991\.  Simple structure in component analysis techniques for mixtures of qualitative and quantitative variables Simple structure in component analysis techniques for mixtures of qualitative and quantitative variables.  Psychometrika562197–212.
  * Krämer Sugiyama (2011) kramerMasashi:2011Krämer, N. Sugiyama, M.  2011\.  The degrees of freedom of partial least squares regression The degrees of freedom of partial least squares regression.  Journal of the American Statistical Association106494697–705.
  * Little Rubin (2002) littleRubin:2002Little, RJA. Rubin, DB.  2002\.  Statistical analysis with missing data Statistical analysis with missing data (2 ).  Hoboken, NJWiley-Interscience.
  * McCrae Costa Jr (1983) mcCraeCosta:1983McCrae, RR. Costa Jr, PT.  1983\.  Joint factors in self-reports and ratings: Neuroticism, extraversion and openness to experience Joint factors in self-reports and ratings: Neuroticism, extraversion and openness to experience.  Personality and individual differences43245–255.
  * Meng (1994) meng:1994Meng, XL.  1994\.  Multiple-imputation inferences with uncongenial sources of input Multiple-imputation inferences with uncongenial sources of input.  Statistical Science94538–558.
  * 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.
  * Oberman Vink (2023) obermanVink:2023Oberman, HI. Vink, G.  2023\.  Towards a standardized evaluation of imputation methodology Towards a standardized evaluation of imputation methodology.  Biometrical Journal2200107.  10.1002/bimj.202200107
  * Park . () parkEtAl:0000Park, S., Ceulemans, E. Van Deun, K.  .  Logistic Regression with Sparse Common and Distinctive Logistic regression with sparse common and distinctive.
  * Robitzsch Grund (2022) miceaddsRobitzsch, A. Grund, S.  2022\.  miceadds: Some Additional Multiple Imputation Functions, Especially for ’mice’ miceadds: Some additional multiple imputation functions, especially for ’mice’ [].  https://CRAN.R-project.org/package=miceadds R package version 3.15-21
  * 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.
  * Schafer Graham (2002) schaferGraham:2002Schafer, JL. Graham, JW.  2002\.  Missing Data: Our View of State of the Art Missing data: Our view of state of the art.  Psychological Methods72147–177.  10.1037//1082-989X.7.2.147
  * 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.
  * Stone Brooks (1990) stoneBrooks:1990Stone, M. Brooks, RJ.  1990\.  Continuum regression: cross-validated sequentially constructed prediction embracing ordinary least squares, partial least squares and principal components regression Continuum regression: cross-validated sequentially constructed prediction embracing ordinary least squares, partial least squares and principal components regression.  Journal of the Royal Statistical Society: Series B (Methodological)522237–258.
  * 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 (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.
  * Van Buuren Groothuis-Oudshoorn (2011) miceVan Buuren, S. Groothuis-Oudshoorn, K.  2011\.  mice: Multivariate Imputation by Chained Equations in R mice: Multivariate imputation by chained equations in R.  Journal of Statistical Software4531–67.
  * Van Buuren Oudshoorn (2000) vanBuurenOudshoorn:2000Van Buuren, S. Oudshoorn, C.  2000\.  Multivariate imputation by chained equations: MICE V1.0 user’s manual Multivariate imputation by chained equations: Mice v1.0 user’s manual PG/VGZ/00.038.  TNO Prevention and Health.
  * Vervloet . (2016) vervloetEtAl:2016Vervloet, M., Van Deun, K., Van den Noortgate, W. Ceulemans, E.  2016\.  Model selection in principal covariates regression Model selection in principal covariates regression.  Chemometrics and Intelligent Laboratory Systems15126-33.  https://www.sciencedirect.com/science/article/pii/S0169743915003093 https://doi.org/10.1016/j.chemolab.2015.12.004
  * Vervloet . (2013) vervloetEtAl:2013Vervloet, M., Van Deun, K., Van den Noortgate, W. Ceulemans, E.  2013\.  On the selection of the weighting parameter value in principal covariates regression On the selection of the weighting parameter value in principal covariates regression.  Chemometrics and Intelligent Laboratory Systems12336–43.
  * Wold (1975) wold:1975Wold, H.  1975\.  Path models with latent variables: The NIPALS approach Path models with latent variables: The nipals approach.  Quantitative sociology Quantitative sociology ( 307–357).  Elsevier.
  * 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) zouEtAL:2006Zou, H., Hastie, T. Tibshirani, R.  2006\.  Sparse principal component analysis Sparse principal component analysis.  Journal of computational and Graphical statistics152265–286.