2 Multiple imputation
The clusterMI package offers various multiple imputation methods dedicated to clustered individuals. They can be divided into two categories: joint modelling (JM) imputation and fully conditional specification (FCS). The first assumes a joint distribution for all variables and imputation is performed using this model. The second proceeds variable per variable in a sequential manner by regression. FCS methods are more time consuming, but also more flexible, allowing a better fit of the imputation model. Multiple imputation is performed with the imputedata function. We start by presenting JM approaches in Section (2.1), while FCS approaches will be presented in Section (2.2).
2.1 Joint modelling imputation
The package proposes two JM methods: JM-GL and JM-DP. Both are based on a multivariate gaussian mixture model.
JM-GLis implemented in themixpackage (Joseph L. Schafer. (2022)). Initially, this method is dedicated to the imputation of mixed data, but it can be used by considering the partition variable as fully incomplete categorical variable. The method assumes constant variance in each cluster (for continuous data).JM-DPis a joint modelling method implemented in the R packagesDPImputeCont(Kim (2020)),NPBayesImputeCat(Wang et al. (2022)),MixedDataImpute(Murray and Reiter (2015)) for continuous, categorical or mixed data respectively. Such a method has the advantage to automatically determine a number of clusters (but this number needs to by bounded bynb.clust). Furthermore, it allows various covariance matrices according to the clusters (for continuous data).
JM-GL is the default imputation method used in imputedata. To perform multiple imputation with this default method, we proceed as follows:
m <- 20 # Number of imputed datasets
res.imp.JM <- imputedata(data.na = wine.na,
nb.clust = nb.clust,
m = m)
#> JM-GL
#> 1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...18...19...20...done!and we specify method = "JM-DP" for imputation using the other JM method:
Note that the imputation requires a pre-specified number of clusters (nb.clust). Here it is tuned to 3 (corresponding to the number of varieties). We explain in Section 4.3 how it can be tuned. Furthermore, the number of imputed data sets is tuned to m = 20 which is generally enough (this choice will be discussed in Section 4.2).
2.2 Fully conditional specification
Fully conditional specification methods consist in a variable per variable imputation. The two fully conditional imputation methods proposed are FCS-homo and FCS-hetero (Audigier, Niang, and Resche-Rigon (2021)). They essentially differ by the assumption about the covariance in each cluster (constant or not respectively).
To perform multiple imputation, we proceed as follows:
maxit <- 20 # Number of iterations for FCS imputation, should be larger in practice
res.imp.FCS <- imputedata(data.na = wine.na,
method = "FCS-homo",
nb.clust = nb.clust,
maxit = maxit,
m = m,
verbose = FALSE)With FCS methods, the imputedata function alternates cluster analysis and imputation given the partition of individuals. When the cluster analysis is performed, the imputedata function calls the mice function from the mice R package (van Buuren and Groothuis-Oudshoorn (2011)). The mice package proposes various methods for imputation. By default, imputedata uses the default method used in mice (predictive mean matching for continuous data), but others can be specified by tuning the method.mice argument. For instance, to imputation under the normal model, use
imputedata(data.na = wine.na,
method = "FCS-homo",
nb.clust = nb.clust,
maxit = maxit,
m = m,
method.mice = "norm")FCS-hetero allows imputation of continuous variables according to linear mixed models (various methods are available in the micemd R package Audigier and Resche-Rigon (2023)). Furthermore, contrary to FCS-homo, FCS-hetero updates the partition without assuming constant variance in each cluster.
2.2.1 Convergence
FCS imputation consists in imputing each variable sequentially several times. Many iterations can be required (maxit argument). For checking convergence, the within and between inertia of each imputed variable can be plotted at each iteration, as proposed by the choosemaxit function
Note that by default, only the five first imputed datasets are plotted (corresponding to the number of curves plotted for each variable). The plotm argument can be tuned to modify which curves should be drawn.
In this case, the number of iterations could be potentially increased. For achieving this goal, the imputation should be rerun by increasing the maxit argument as follows:
2.2.2 Specifying imputation models
Fully conditional imputation methods are quickly limited when the number of variables is large since imputation models become overfit. To address this issue, we can use penalised regression as proposed in mice by specifying method.mice = lasso.norm for instance. Another way consists in specifying conditional imputation models by tuning the predictmat argument. This argument is a binary matrix where each row indicates which explanatory variables (in column) should be used for imputation.
2.2.2.1 The varselbest procedure
To tune this matrix in an automatic way, the varselbest function proposes to perform variable selection following Bar-Hen and Audigier (2022). Briefly, varselbest performs variable selection on random subsets of variables and then, combines them to recover which explanatory variables are related to the response. More precisely, the outlines of the algorithm are as follows: let consider a random subset of sizeblock among p variables. Then, any selection variable scheme can be applied (lasso, stepwise and knockoff are proposed by tuning the method.select argument). By resampling (B times) a sample of size sizeblock among the p variables, we may count how many times a variable is considered as significantly related to the response and how many times it is not. We need to define a threshold (r) to conclude if a given variable is significantly related to the response (by default, r = 0.3). The main advantage of this function is that it handles both missing values and high-dimensional data.
By default, the varselbest function performs variable selection by knockoff (Barber and Candès (2015)) based on B = 200 bootstrap subsets. The threshold r is tuned at 0.3 allowing to omit only variables very poorly predictive. The choices of B and r are discussed in next sections.
Since the method is time consuming, the function allows parallel computing by tuning the nnodes argument. In the next example, the imputation model for the variable alco is obtained using the algorithm previously described.
nnodes <- 2
# Number of CPU cores used for parallel computation.
# Use parallel::detectCores() to choose an appropriate number
# variable selection to impute the "alco" variable
B <- 50 # number of bootstrap subsets, should be increased in practice
res.varsel <- varselbest(res.imputedata = res.imp.FCS, B = B, listvar = "alco",
nnodes = nnodes, graph = FALSE)
res.varsel$predictormatrix["alco", ]#> alco malic ash alca mg phe fla nfla pro col hue ratio prol
#> 0 0 1 1 1 0 1 0 1 1 1 1 1The function suggests considering the variables ash, alca, mg, fla, pro, col, hue, ratio, prol to impute the alco variable. Then, imputation can be rerun by specifying the predictmat argument returned by the varselbest function
# multiple imputation with the new model
res.imp.select <- imputedata(data.na = wine.na,
method = "FCS-homo",
nb.clust = nb.clust,
maxit = maxit,
m = m,
predictmat = res.varsel$predictormatrix,
verbose = FALSE)Note that for specifying all conditional imputation models you should use
2.2.2.2 Convergence
The number of iterations B should be large so that the proportion of times a variable is selected becomes stable. The chooseB function plots the proportion according to the number of iterations.
B=50 iterations seems not enough.
2.2.2.3 Tuning the threshold r
To tune the threshold r, the vector of proportions can be graphically investigated.
# check the variable importance
round(res.varsel$proportion["alco",], 2)
#> alco malic ash alca mg phe fla nfla pro col hue ratio prol
#> 0.00 0.26 0.36 0.45 0.33 0.10 0.33 0.14 0.45 0.95 0.71 0.60 0.36
barplot(sort(res.varsel$proportion["alco",], decreasing=TRUE),
ylab = "proportion",
main = "alco",
ylim = c(0, 1),
las = 2,
cex.names = .5)
r <- 0.2 # a new threshold value (r = 0.3 by default)
abline(h = r, col = 2, lty = 2)
Then, the predictor matrix can be easily updated at hand as follows
predictormatrix <- res.varsel$predictormatrix
predictormatrix[res.varsel$proportion>r] <- 1
predictormatrix[res.varsel$proportion<=r] <- 0
predictormatrix["alco", ]
#> alco malic ash alca mg phe fla nfla pro col hue ratio prol
#> 0 1 1 1 1 0 1 0 1 1 1 1 1By decreasing the threshold, more variables are used in the imputation model.
A more automatic way to tune r consists in using the function chooser. chooser computes an optimal threshold using K-fold cross-validation. The call to chooser is highly time consuming, but the optimal value of r can be follows during the process through graphical outputs. By this way, the user can stop the process early. This can be achieved as follows:
On the wine dataset, the number of clusters suggested is clearly 3.
