Multiple imputations have become one of the standard methods in drawing inferences in many in- complete data applications. Applications of multiple imputations in relatively more complex settings with clustered high-dimensional data require specialized methods to overcome the computational burden. Using mixed-effects models, we develop methods that can be applied to continuous, binary, or categorical incomplete data. These methods are computationally feasible as they rely on variational Bayesian inference for sampling the posterior predictive distribution missing data. These methods specifically target high-dimensional covariates and work with spike-and-slab priors. The individual regression computation is then incorporated in a sequential hierarchical regression imputation. To further improve the computational efficiency for the categorical data, we apply the calibration-based imputation methods using the approximating variational distribution. We present a simulation study to validate the usefulness of the proposed approach.