BBGC: Bayesian Bootstrap Gaussian Copula Model

• BBGC is a fully Bayesian, nonparametric Gaussian copula framework that integrates Bayesian bootstrap priors for robust imputation of mixed continuous, ordinal, and binary data.
• It employs an ordinal-threshold construction and random marginal distributions to preserve joint dependence and accurately propagate uncertainty under MCAR and MAR conditions.
• Simulation and real-data experiments demonstrate that BBGC often achieves lower NRMSE and improved uncertainty quantification compared to traditional imputation methods at moderate to high missingness rates.

The Bayesian Bootstrap-based Gaussian Copula model (BBGC) is a fully Bayesian, nonparametric Gaussian copula framework for imputing missing values in multivariate mixed data, with continuous and ordinal variables and binary variables treated as a special case of ordinal. It combines a Gaussian copula for joint dependence, Bayesian bootstrap priors for each marginal cumulative distribution function, and an ordinal-threshold construction for mixed-variable integration. BBGC is formulated for missingness mechanisms that are MCAR or MAR, and is motivated in particular by settings with high missingness rates, where treating marginal distributions as fixed can induce biased marginal estimates and overconfident imputations (Kim et al., 9 Jul 2025). In a broader methodological sense, it is a concrete realization of Bayesian copula modelling, which separates marginal behaviour from joint dependence and treats dependence parameters as objects of Bayesian inference (Smith, 2011).

1. Conceptual setting and motivation

Copula-based modelling starts from Sklar’s theorem, which decomposes a multivariate distribution into marginal CDFs and a copula. In the Gaussian copula case,

$$\left(\Phi^{-1}(F_1(X_1)), \ldots, \Phi^{-1}(F_p(X_p))\right)^T \sim N(0,R),$$

where $F_j$ are marginal CDFs and $R$ is a correlation matrix. This decomposition is central to Bayesian copula analysis because it permits the marginal components and the dependence structure to be modelled separately, while preserving a coherent joint distribution (Smith, 2011).

BBGC is motivated by the observation that many standard imputation procedures either distort dependence or fail to propagate uncertainty. Deletion methods can be acceptable when the missing rate is below $15\%$, but become biased and inefficient at higher rates. Mean imputation, regression imputation, and EM-based single imputation tend to underestimate variability and distort dependencies. MICE is flexible but relies on specifying conditional models variable by variable and may struggle with complex multivariate dependence, nonlinearity, and mixed data. KNN and missForest can capture nonlinearity but do not explicitly account for joint uncertainty, and their performance degrades with very high missingness. Matrix factorization and autoencoder approaches are sensitive to model tuning and are not designed primarily to quantify uncertainty. Dirichlet Process Mixture Models are Bayesian nonparametric but computationally heavy in high-dimensional mixed-data settings (Kim et al., 9 Jul 2025).

Within Gaussian copula imputation, the principal limitation identified by BBGC is the treatment of the marginals. Existing approaches often plug in empirical CDFs or use rank-likelihood methods that avoid estimating $F_j$ explicitly. In either case, uncertainty in the marginal CDFs is not propagated into inference on $R$ or into posterior predictive imputation. Under high missingness, the observed sample for each variable is small, empirical marginals become unstable, and treating them as fixed yields overconfident and potentially biased imputations (Kim et al., 9 Jul 2025).

2. Formal specification of the BBGC model

Let $X=(x_{ij})$ be an $n\times p$ data matrix, with variables $j=1,\dots,p$ and units $i=1,\dots,n$. Continuous variables satisfy $F_j$0, ordinal variables satisfy $F_j$1, and binary variables are treated as ordinal with $F_j$2. The observed and missing entries are partitioned into $F_j$3, with

$F_j$4

BBGC is formulated under MAR, with MCAR treated as a special case; MNAR is not treated (Kim et al., 9 Jul 2025).

The copula component is defined through a latent Gaussian vector. Let

$F_j$5

where $F_j$6 is a correlation matrix. The joint CDF then takes the Gaussian copula form

$F_j$7

The distinctive feature of BBGC is its treatment of the marginals $F_j$8. For a single variable $F_j$9 with observed values $R$0, the Bayesian bootstrap defines

$R$1

This random CDF is supported on the observed values, with random Dirichlet weights. When $R$2, it coincides with the empirical CDF.

BBGC does not use $R$3 directly. It introduces the adjusted transform

$R$4

If

$R$5

then $R$6 is sign-invariant. The chosen adjustment is

$R$7

This serves two stated purposes. First, it avoids boundary divergence because $R$8, so $R$9 remains finite. Second, it restores sign invariance in expectation (Kim et al., 9 Jul 2025).

A finite-sample guarantee is given for the adjusted Bayesian bootstrap CDF. If $15\%$0 with $15\%$1 continuous, then for all $15\%$2 and $15\%$3,

$15\%$4

This establishes a non-asymptotic bound on the sup-norm distance between $15\%$5 and the true CDF (Kim et al., 9 Jul 2025).

The hierarchical BBGC model can be summarized as

$15\%$6

$15\%$7

$15\%$8

3. Representation of continuous and ordinal variables

For continuous variables, BBGC uses a monotone Gaussianization. For each observed continuous entry $15\%$9,

$F_j$0

There is therefore a one-to-one mapping between an observed continuous value and its latent Gaussian score through the random marginal CDF (Kim et al., 9 Jul 2025).

For ordinal variables, including binary variables, BBGC uses a threshold construction. If $F_j$1, define cutoffs

$F_j$2

with

$F_j$3

Then

$F_j$4

For observed ordinal entries, the latent Gaussian variable is constrained to the corresponding interval, and for imputation or prediction the mapping back to the data scale is defined by the same interval rule (Kim et al., 9 Jul 2025).

The observed-data mechanism enters through the MAR posterior decomposition

$F_j$5

BBGC approximates $F_j$6 by drawing each marginal independently from the Bayesian bootstrap posterior. Conditional on $F_j$7, inference proceeds on $F_j$8 and the latent Gaussian variables $F_j$9, and the final posterior is obtained by Bayesian model averaging across multiple $R$0-draws (Kim et al., 9 Jul 2025).

This construction directly targets the main deficiency of fixed-marginal Gaussian copula methods. The dependence model remains Gaussian-copular, but the map $R$1 is itself random because the $R$2 are random. A plausible implication is that posterior uncertainty in the marginals is transmitted to posterior uncertainty in both $R$3 and the imputed values, rather than being suppressed by a deterministic empirical transform.

4. Posterior inference, Gibbs sampling, and imputation

Conditional on a draw $R$4 of the marginal CDFs, the joint posterior has the form

$R$5

Inference is carried out by Gibbs sampling with three blocks: latent $R$6, correlation $R$7, and missing $R$8 (Kim et al., 9 Jul 2025).

For each observation $R$9 and variable $X=(x_{ij})$0, the full conditional of $X=(x_{ij})$1 given $X=(x_{ij})$2 and $X=(x_{ij})$3 is

$X=(x_{ij})$4

with

$X=(x_{ij})$5

This Gaussian conditional is then modified according to data type and missingness pattern. For observed continuous entries, $X=(x_{ij})$6 is deterministic: $X=(x_{ij})$7 For missing continuous entries,

$X=(x_{ij})$8

without truncation. For observed ordinal entries in category $X=(x_{ij})$9,

$n\times p$0

and for missing ordinal entries the conditional is again untruncated Gaussian (Kim et al., 9 Jul 2025).

Given the latent vectors $n\times p$1, the dependence matrix is updated by first sampling a covariance matrix

$n\times p$2

and then normalizing: $n\times p$3 This guarantees that $n\times p$4 is symmetric positive-definite with unit diagonal (Kim et al., 9 Jul 2025).

Missing values are imputed from the posterior predictive distribution conditional on $n\times p$5 and $n\times p$6. For continuous variables,

$n\times p$7

Because $n\times p$8 is discrete, the implementation uses the empirical inverse CDF of $n\times p$9. For ordinal variables,

$j=1,\dots,p$0

Marginal uncertainty is integrated by repeating this Gibbs procedure over multiple independent draws of $j=1,\dots,p$1. The conceptual workflow is: sample Dirichlet weights for each variable, construct $j=1,\dots,p$2, run the Gibbs sampler conditional on $j=1,\dots,p$3, and aggregate imputations or parameter summaries across runs. The paper does not provide pseudocode, but this is the stated inferential structure (Kim et al., 9 Jul 2025).

The computational bottleneck is matrix algebra in dimension $j=1,\dots,p$4. Each Gibbs iteration requires conditional-Gaussian calculations involving $j=1,\dots,p$5, and sampling $j=1,\dots,p$6 from an inverse-Wishart is also $j=1,\dots,p$7. The authors explicitly note that posterior sampling is computationally expensive, particularly because inversion of the correlation matrix poses challenges in high-dimensional settings. Truncated normal sampling for many ordinal variables adds further overhead. Code is provided at https://github.com/zlatjdals/BBGC (Kim et al., 9 Jul 2025).

5. Simulation studies and empirical performance

The simulation study uses $j=1,\dots,p$8 with correlation structure

$j=1,\dots,p$9

Latent normals are generated from $i=1,\dots,n$0. Variables $i=1,\dots,n$1 through $i=1,\dots,n$2 are ordinal or discrete normals,

$i=1,\dots,n$3

variables $i=1,\dots,n$4 through $i=1,\dots,n$5 are continuous uniforms,

$i=1,\dots,n$6

and variables $i=1,\dots,n$7 through $i=1,\dots,n$8 are continuous exponentials,

$i=1,\dots,n$9

Artificial missingness is imposed at rates $F_j$00 under both MCAR and MAR, with $F_j$01 replications per scenario. Performance is measured by normalized RMSE,

$F_j$02

Lower values are better (Kim et al., 9 Jul 2025).

In those simulations, BBGC achieves the best NRMSE among all methods at missing rates $F_j$03, $F_j$04, and $F_j$05, under both MCAR and MAR. At the very high missing rate of $F_j$06, mean imputation yields the lowest NRMSE, while BBGC is close but slightly worse, and missForest and other methods degrade more strongly. The reported interpretation is that, at extremely high missingness, modelling complex dependence may become unstable and a very simple imputation can be more robust in the NRMSE sense (Kim et al., 9 Jul 2025).

The paper also evaluates uncertainty quantification for the adjusted Bayesian bootstrap marginals under $F_j$07 MCAR missingness. For three representative variables—$F_j$08, $F_j$09, and $F_j$10—many draws of $F_j$11 are generated and pointwise $F_j$12 credible bands are formed using the $F_j$13th and $F_j$14th percentiles. The true CDF lies within the $F_j$15 credible band over most of the range, and coverage of the true CDF at observed points is reported as $F_j$16 for $F_j$17, $F_j$18 for $F_j$19, and $F_j$20 for $F_j$21 (Kim et al., 9 Jul 2025).

Three real-data experiments further characterize empirical behaviour.

Dataset	Design	Reported outcome
Wine Quality (Red Wine)	$F_j$22, $F_j$23; 3 ordinal variables; $F_j$24–$F_j$25 MCAR/MAR	missForest best at $F_j$26 and $F_j$27; BBGC best at $F_j$28 and $F_j$29
Breast Cancer Wisconsin (Diagnostic)	$F_j$30, $F_j$31; 1 binary and 30 continuous variables	BBGC lowest NRMSE across all missingness rates and mechanisms
Semiconductor manufacturing process data	$F_j$32, $F_j$33; intrinsic missing rate $F_j$34	BBGC lowest NRMSE for additional 50, 100, or 150 masked entries

For the Wine Quality dataset, DPMM is reported to be close to BBGC across all missing rates, especially at high missingness. For the Breast Cancer dataset, MICE is the second-best method at $F_j$35 missingness, whereas DPMM becomes the second-best at $F_j$36, but both remain worse than BBGC. For the semiconductor dataset, DPMM is again second-best, and mean imputation nearly matches missForest in some scenarios, underscoring the difficulty of the problem. BBGC nonetheless maintains the lowest NRMSE despite the overall missingness of approximately $F_j$37 (Kim et al., 9 Jul 2025).

6. Relation to Bayesian copula methodology, limitations, and extensions

BBGC sits naturally within the larger Bayesian copula framework in which the joint model is decomposed into marginal distributions and a copula, and Bayesian inference targets the dependence parameters through likelihood-based posterior sampling (Smith, 2011). In that broader framework, Gaussian copulas are parameterized by a correlation matrix, Bayesian methods are used when likelihood optimization is difficult or model selection is required, and data augmentation generalizes multivariate probit and related latent-Gaussian models to broader copula constructions for discrete and ordinal data (Smith, 2011).

Earlier Bayesian copula work explicitly allows margins to be modelled nonparametrically and notes that Bayesian nonparametric methods can be used for $F_j$38, while retaining Gaussian-copula or vine-based dependence modelling (Smith, 2011). BBGC can therefore be viewed as a specific nonparametric marginal construction inside that program: the copula remains Gaussian, but the marginals are endowed with Bayesian bootstrap priors rather than being fixed at empirical CDFs. This suggests that BBGC is less a departure from Bayesian copula modelling than a particular answer to the marginal-uncertainty problem in high-missing mixed data.

The limitations reported for BBGC are fourfold. First, computational cost is substantial because inversion and sampling of a $F_j$39 covariance matrix is $F_j$40, and truncated-normal updates for many ordinal variables add further expense. Second, the dependence structure is constrained to a Gaussian copula, so performance may degrade when the true copula exhibits behaviour such as heavy tail dependence. Third, scalability is limited because no built-in sparsity or factor structure is imposed on $F_j$41. Fourth, the missingness assumption is MAR; MNAR is not handled (Kim et al., 9 Jul 2025).

The authors suggest incorporating structural assumptions on $F_j$42, including sparsity through graphical models or shrinkage priors, and low-rank factorization through Gaussian copula factor models such as Murray et al. (2013), in order to improve scalability (Kim et al., 9 Jul 2025). Other natural extensions, identified as consistent with the framework rather than fully developed within it, include replacing the Gaussian copula with alternatives such as a $F_j$43-copula, extending to time series or longitudinal settings via dynamic copulas or state-space models, and explicitly modelling MNAR mechanisms by coupling the missingness model with the copula-based data model.