C-IMHI: High-Dimensional Imputation

• C-IMHI is a Bayesian high-dimensional imputation method that addresses missing neuroimaging data by utilizing robust posterior-mode covariance estimation.
• It significantly reduces computational cost from O(p³) to O(p) per iteration while maintaining stability when the number of variables far exceeds the sample size.
• Empirical results show that C-IMHI outperforms traditional methods such as MICE, achieving lower error rates and rapid convergence in neuroimaging studies.

High-Dimensional Multiple Imputation (HIMA) is a Bayesian methodology for addressing missing data in multivariate neuroimaging datasets characterized by high dimensionality and complex dependence structures. Traditional multiple imputation approaches are computationally prohibitive and unstable in such settings due to the size and properties of the covariance matrices involved. HIMA introduces a computationally efficient and numerically stable strategy for multiple imputation by leveraging robust posterior-mode-based covariance estimation, specifically tailored for scenarios where the number of variables far exceeds the sample size (Lu et al., 2023).

1. Bayesian Model Specification

Let $Y_i = (y_{i1},...,y_{ip})'$ denote the $p$-dimensional neuroimaging measurement vector for subject $i = 1,...,n$. Each vector $Y_i$ is partitioned into observed ($Y_i^{obs} \in \mathbb{R}^{p_o}$) and missing ($Y_i^{mis} \in \mathbb{R}^{p_m}$) blocks, such that $p_o + p_m = p$. A joint multivariate normal model is assumed under a Missing At Random (MAR) mechanism: $$(Y_i^{obs}, Y_i^{mis})' \sim N_p( \mu = [\mu_o; \mu_m],\, \Sigma )$$ with

$$\Sigma = \begin{pmatrix} \Sigma_{oo} & \Sigma_{om} \ \Sigma_{mo} & \Sigma_{mm} \end{pmatrix}$$

The observed-data likelihood is given by: $$L(\mu, \Sigma; Y^{obs}) = \prod_{i=1}^n (2\pi)^{-p_o/2} |\Sigma_{oo}|^{-1/2} \exp\left\{-\frac{1}{2} (Y_i^{obs} - \mu_o)' \Sigma_{oo}^{-1} (Y_i^{obs} - \mu_o)\right\}$$ Conjugate priors are adopted: $p$0 yielding a joint posterior: $p$1 with $p$2 as the sample covariance and $p$3 as the mean of the completed data. The full conditional distributions are: $p$4

2. Computational Strategy for High-Dimensional Covariance Estimation

Direct posterior sampling of $p$5 from an Inverse-Wishart distribution is computationally expensive—the cost is $p$6 per draw due to the requisite matrix inversion and Cholesky decomposition—and numerically unstable when $p$7 due to near-singularity. HIMA circumvents this by substituting full posterior sampling with repeated computation of the posterior-mode (MAP) estimate $p$8.

The posterior mode under an Inv–Wishart prior is: $p$9 Reparameterizing the prior as $i = 1,...,n$0, where $i = 1,...,n$1 is the prior precision and $i = 1,...,n$2 the prior mean, and setting $i = 1,...,n$3 yields: $i = 1,...,n$4 Hyperparameters $i = 1,...,n$5 are estimated by empirical Bayes via minimization of the Kullback–Leibler divergence: $i = 1,...,n$6 Closed-form updates for this procedure are available. Theoretical results establish consistency of $i = 1,...,n$7 in operator norm even in extreme high-dimensional regimes ($i = 1,...,n$8). Replacing posterior draws with mode estimation reduces per-iteration computational cost to $i = 1,...,n$9 and enhances numerical stability by avoiding degeneracy or near-singularity of the covariance estimate (Lu et al., 2023).

3. HIMA Algorithmic Workflow

HIMA generates $Y_i$0 multiply imputed datasets as follows:

1. Initialization: For each imputation $Y_i$1, initialize $Y_i$2 with marginal column means. Set $Y_i$3, $Y_i$4.
2. Iterative Steps (for $Y_i$5):
   • Imputation Step (I-step): For each $Y_i$6, sample missing values $Y_i$7 from the conditional normal distribution:

   $Y_i$8

   where

   $Y_i$9

   $Y_i^{obs} \in \mathbb{R}^{p_o}$0

- Parameter Step (P-step): Compute sample covariance $Y_i^{obs} \in \mathbb{R}^{p_o}$1 and mean $Y_i^{obs} \in \mathbb{R}^{p_o}$2 of the completed data. Update covariance:

$Y_i^{obs} \in \mathbb{R}^{p_o}$3

Sample new mean:

$Y_i^{obs} \in \mathbb{R}^{p_o}$4

1. Output: After $Y_i^{obs} \in \mathbb{R}^{p_o}$5 iterations (typically $Y_i^{obs} \in \mathbb{R}^{p_o}$6–$Y_i^{obs} \in \mathbb{R}^{p_o}$7 suffices), $Y_i^{obs} \in \mathbb{R}^{p_o}$8 is recorded as the $Y_i^{obs} \in \mathbb{R}^{p_o}$9-th imputed dataset.

Convergence is assessed via trace plots of selected entries in $Y_i^{mis} \in \mathbb{R}^{p_m}$0, $Y_i^{mis} \in \mathbb{R}^{p_m}$1, or representative imputed values. The average number of iterations required for stabilization is approximately $Y_i^{mis} \in \mathbb{R}^{p_m}$2 in empirical settings (Lu et al., 2023).

4. Empirical Performance and Evaluation

4.1 Simulation Studies

Simulation studies use real MRI data: $Y_i^{mis} \in \mathbb{R}^{p_m}$3 subjects, $Y_i^{mis} \in \mathbb{R}^{p_m}$4--$Y_i^{mis} \in \mathbb{R}^{p_m}$5 voxels per region. Entries are randomly deleted to achieve $Y_i^{mis} \in \mathbb{R}^{p_m}$620% missingness per voxel, with variables exceeding 40% missingness removed. The following metrics are reported, standardized per voxel variance:

• Weighted Mean Absolute Error (wMAE)
• Weighted Mean Squared Error (wMSE)
• Weighted Mean Bias Error (wMBE)

HIMA is benchmarked against mean-substitution and MICE (Multivariate Imputation by Chained Equations, $Y_i^{mis} \in \mathbb{R}^{p_m}$7, $Y_i^{mis} \in \mathbb{R}^{p_m}$8). Results for the right hippocampus region ($Y_i^{mis} \in \mathbb{R}^{p_m}$9) are summarized below.

Method	Mean wMAE (sd)	Mean wMSE (sd)	Mean wMBE (sd)	Time (hrs)
Mean-imp	0.452 (0.02)	0.320 (0.03)	0.041 (0.01)	<0.1
MICE	0.129 (0.01)	0.046 (0.01)	0.008 (0.004)	2300
HIMA	0.072 (0.005)	0.018 (0.002)	0.004 (0.002)	1.1

HIMA reduces error by approximately 40% relative to MICE and gives more than 2000-fold speedup in computational time (Lu et al., 2023). Trace plots demonstrate that HIMA achieves smaller posterior variance and rapid convergence in about 10 iterations.

4.2 Application to Real Schizophrenia MRI Data

Applied to 99 ICBM brain regions (missing-rate 14.6%, $p_o + p_m = p$0--$p_o + p_m = p$1), HIMA is run with $p_o + p_m = p$2, $p_o + p_m = p$3. Imputation quality is verified via:

• Scatterplots of observed vs imputed values for randomly sampled voxels, which cluster near the diagonal.
• Overlapping kernel density estimates for observed and imputed marginal distributions.
• Whole-brain heatmaps showing restoration of coverage in data-sparse anatomical regions, especially frontal and temporal lobes.

The imputed datasets enable more stable group-level analysis of water-exchange function in schizophrenia by restoring statistical power in key brain regions (Lu et al., 2023).

5. Limitations and Extensions

HIMA relies on a multivariate normality assumption, which may be violated for non-Gaussian contrasts (e.g., task-fMRI). The approach loses stability when more than 40% data are missing per variable. It does not explicitly model spatial or graphical structure in the covariance matrix beyond what is captured in $p_o + p_m = p$4; for example, there is no use of Gaussian-Markov random fields.

Potential extensions include introducing structured (spatial/graphical) priors on $p_o + p_m = p$5—such as banded or sparse-precision matrices—to better represent spatial dependencies; and dimensionality reduction (e.g., PCA) for very high variable counts ($p_o + p_m = p$6).

For successful application, the following recommendations are provided:

• Variables with $p_o + p_m = p$7 missingness should be dropped before applying HIMA.
• Marginal densities for each subject may be smoothed to mitigate inter-subject distributional variation.
• Convergence diagnostics should be run using trace plots for selected variables.
• For typical downstream usage, $p_o + p_m = p$8–$p_o + p_m = p$9 and $$(Y_i^{obs}, Y_i^{mis})' \sim N_p( \mu = [\mu_o; \mu_m],\, \Sigma )$$0–$$(Y_i^{obs}, Y_i^{mis})' \sim N_p( \mu = [\mu_o; \mu_m],\, \Sigma )$$1 suffice, though $$(Y_i^{obs}, Y_i^{mis})' \sim N_p( \mu = [\mu_o; \mu_m],\, \Sigma )$$2 may be increased for improved uncertainty quantification.

6. Significance and Domain Impact

HIMA delivers a numerically stable, $$(Y_i^{obs}, Y_i^{mis})' \sim N_p( \mu = [\mu_o; \mu_m],\, \Sigma )$$3 per-iteration computational solution for high-dimensional covariance estimation and multiple imputation. It achieves substantial speed advantages over traditional Gibbs/MICE approaches while also improving the accuracy and stability of imputations in high-dimensional neuroimaging data (Lu et al., 2023). Restoring missing data in neuroimaging studies of disorders such as schizophrenia enables more comprehensive and reliable group-level statistical analyses, particularly in anatomically large or data-sparse brain regions. The method is broadly adaptable to multivariate neuroimaging data with similar missingness and dimensionality profiles, with possible extensions to incorporate spatial structure and alternative distributional assumptions.