Hybrid copula and conditional CART/GAM synthesis for nonlinear relationships

Scope

Elliptical copulas are useful joint models, but they do not necessarily reproduce nonlinear conditional relationships, even if their marginal quantiles and pairwise rank correlations are accurate. SynthesizerPlus offers an optional, experimental semi-parametric hybrid to address the narrower problem of synthesizing one continuous response conditionally on a set of copula-generated predictors. The default remains the unmodified copula.

Let \(X\) denote the predictors and \(Y\) the target. A regression tree estimates \(\widehat m(X)\); within a group of training observations with the same fitted terminal prediction, we retain centered empirical residuals. Given a new copula-generated predictor vector \(X^*\), the provisional target is

\[Y^*_{\mathrm{raw}} = \widehat m(X^*) + \varepsilon^*, \]

where \(\varepsilon^*\) is sampled from the empirical residual pool corresponding to \(\widehat m(X^*)\). A final quantile calibration maps \(Y^*_{\mathrm{raw}}\) to the observed empirical marginal distribution of \(Y\):

\[Y^* = \widehat F_Y^{-1}\bigl(\widetilde G(Y^*_{\mathrm{raw}})\bigr),\]

where \(\widehat F_Y^{-1}\) is the empirical quantile function of the observed response and \(\widetilde G\) is the mid-point empirical distribution function of a reference sample of 20000 raw conditional draws, generated once at fitting time for resampled training predictor rows. Because each value is calibrated against this fixed reference, and not against the other values generated in the same call, the result does not depend on the batch size: generating one row at a time, in chunks, or within small strata gives the same distribution. The calibration is monotone, so it preserves the ordering of the conditional draws, but it can alter conditional means. The hybrid is not equivalent to sampling an elliptical copula.

Predictor categories are encoded from all training rows. A category that occurs only in rows with a missing response can therefore be generated by the copula; the tree treats it as unknown and routes it by surrogate splits.

Reproducible use case

library(SynthesizerPlus)
set.seed(31)
x <- runif(500, -2, 2)
d <- data.frame(x = x, y = sin(3 * x) + 0.3 * rnorm(500))

base <- fit_synthesizer(d, seed = 1)
hybrid <- fit_synthesizer(d, nonlinear_target = "y",
                          nonlinear_minbucket = 15, seed = 1)
original_draw <- generate(base, n = 500, seed = 2)
hybrid_draw <- generate(hybrid, n = 500, seed = 2)

# Compare nonlinear curves and marginal error without treating a plot as
# proof of general superiority.
plot(d$x, d$y, col = grDevices::adjustcolor("black", 0.25),
     xlab = "x", ylab = "y", main = "Nonlinear response")
points(original_draw$x, original_draw$y, col = "red", pch = 3)
points(hybrid_draw$x, hybrid_draw$y, col = "blue", pch = 1)
legend("topleft", c("original", "copula", "hybrid"),
       col = c("black", "red", "blue"), pch = c(1, 3, 1), bty = "n")

Limitations and responsible use

See system.file("examples", "hybrid_nonlinear_benchmark.R", package = "SynthesizerPlus") for a Monte Carlo comparison template.

Controlled ablation: what drives the improvement?

The hybrid method contains three distinguishable components: a fitted CART conditional mean, a stochastic residual generator, and optional empirical quantile calibration of the response. The default hybrid uses nonlinear_residual = "leaf" and nonlinear_calibration = "quantile". To separate the effects, select a residual generator (leaf, global, or none) and a calibration strategy (quantile or none):

fit <- fit_synthesizer(d, nonlinear_target = "y", seed = 101,
                       nonlinear_residual = "global",
                       nonlinear_calibration = "none")
synthetic <- generate(fit, n = nrow(d), seed = 202)

With none residuals, the generated response takes only the fitted tree’s terminal predictions before calibration. With global residuals, the procedure ignores differences in conditional noise variance across leaves. quantile calibration reproduces the unconditional marginal up to Monte Carlo error but may alter conditional means, variances and interactions. none calibration retains raw conditional draws, which can extend outside the training response range. These settings are experimental analysis controls, not general claims of improved statistical fidelity or disclosure control. The script system.file("examples", "hybrid_nonlinear_benchmark.R", package = "SynthesizerPlus") is a template for paired comparisons on held-out data. Fitted trees and residual pools are not safe to publish with confidential input data.

GAM conditional engine

The original optional CART hybrid is a piecewise-constant estimate of the conditional mean. For smooth nonlinear relationships, a penalised GAM may provide a better approximation without changing the copula model used for the predictors. Fit it by selecting nonlinear_engine = "gam"; the default remains "cart".

For a continuous target \(Y\), the additive variant models

\[ E(Y\mid X) = \beta_0 + \sum_j f_j(X_j) + \sum_k \gamma_k I(X_{\mathrm{cat}}=k), \]

where the numeric terms use penalised cubic regression splines fitted by mgcv::gam(..., method = "REML", select = TRUE). The optional nonlinear_gam_interactions = TRUE adds tensor-product interactions between eligible smooth numerical predictors (at most three).

fit_gam <- fit_synthesizer(d, nonlinear_target = "y",
                           nonlinear_engine = "gam",
                           nonlinear_gam_k = 10,
                           nonlinear_residual = "leaf",
                           nonlinear_calibration = "quantile", seed = 1)
syn_gam <- generate(fit_gam, n = 500, seed = 2)
stopifnot(nrow(syn_gam) == 500)

A GAM has no terminal tree leaves. With the existing default nonlinear_residual = "leaf", the GAM engine instead samples residuals from the nonlinear_local_k training predictions nearest to the generated fitted mean; the selected pool is centred. This provides a local conditional approximation, not a fully nonparametric conditional distribution estimator. It cannot generally capture noise variance varying in directions orthogonal to the fitted mean. "global" and "none" remain available as controlled ablations. "quantile" calibration is retained by default.

The saved GAM model contains training-derived smooth coefficients and residual information, and must not be treated as disclosure-safe. Comparisons against CART, original copula, synthesizer and synthpop should be performed on independent holdouts and multiple data-generating processes before drawing conclusions about relative utility.

Automatic conditional-engine selection (experimental)

nonlinear_engine = "auto" selects between CART, additive GAM, and GAM with pairwise smooth interactions by training-only K-fold cross-validation. By default each candidate is scored by the continuous ranked probability score (CRPS) of its conditional predictive distribution, i.e. the fitted conditional mean together with the residual resampling used for generation, evaluated on the held-out rows. The CRPS penalises a wrong conditional spread as well as a wrong conditional mean; nonlinear_cv_score = "mse" scores the conditional mean only. The default one-standard-error rule selects the simplest model whose mean CV loss is within one fold-level standard error of the best candidate. A model is not considered if any CV fold fails. The final selected engine is fitted anew on the full training data before generating synthetic observations.

fit_auto <- fit_synthesizer(
  d, nonlinear_target = "y", nonlinear_engine = "auto",
  nonlinear_cv_folds = 5L, nonlinear_cv_rule = "one_se",
  seed = 101
)
fit_auto$pooled$nonlinear$selection$scores
fit_auto$pooled$nonlinear$selection$selected
syn_auto <- generate(fit_auto, n = nrow(d), seed = 102)

For a held-out response \(y\) with \(M = 50\) conditional draws \(X_1, \dots, X_M\) from the fold model, the score is

\[\mathrm{CRPS} = \frac{1}{M}\sum_{i}|X_i - y| - \frac{1}{2M^2}\sum_{i,j}|X_i - X_j|,\]

averaged over the held-out rows. The table selection$scores reports mean_score and se_score for the criterion used, together with mean_mse and se_mse; the per-fold values are in selection$fold_scores and selection$fold_mse.

In a stratified fit each stratum model with enough observed responses has its own conditional model (and, with "auto", its own selection). Strata without one generate the response by their copula; fit_synthesizer() warns and lists them in fit$strata$conditional_fallback.

The selector uses only the original training rows; it does not inspect any external evaluation or benchmark holdout. CV scores describe the conditional distribution of the response, not the joint quality of the synthetic data, confidence intervals or privacy. The chosen engine can be suboptimal for those other purposes; when the conditional mean improvement of a GAM outweighs its coarser residual model, the CRPS can prefer the GAM even if CART reproduces heteroskedasticity better. The auto path needs rpart and uses mgcv when installed. A sparse stratum that permits fitting the full CART model but lacks enough rows for CV falls back to CART and records its reason.

Automatic selection is expected to be slower than a fixed conditional engine, particularly if GAM interactions are considered. The default nonlinear_engine = "cart" remains unchanged for existing applications.