--- title: "Hybrid copula and conditional CART/GAM synthesis for nonlinear relationships" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Hybrid copula and conditional CART/GAM synthesis for nonlinear relationships} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") ``` ## 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 ```{r, eval=requireNamespace("rpart", quietly=TRUE)} 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 - Supports one double-precision continuous target per model (possibly with missing observations), with all other modelled columns as predictors. - Missing predictors are encoded by median-imputation plus explicit indicators (numeric) or an extra category (categorical). - The copula still models predictor distributions and missingness; the conditional response modification may change the final dependence. - The selected response is calibrated to its **observed** marginal. Under nonignorable missingness this need not estimate the population marginal. - Non-default `closeness`, `dependence`, and `perturb` controls are rejected. - The stored tree and residual pools are **not privacy protected**. Do not publish the fitted model or treat this option as differential privacy. - Tree smoothing can fail on high-dimensional or highly extrapolative data. Evaluate utility on an independent test set and compare against the original copula and established sequential synthesis tools. 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`): ```{r hybrid-ablation-demo, eval=FALSE} 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). ```{r gam-hybrid, eval=requireNamespace("mgcv", quietly=TRUE)} 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. ```{r auto-engine-example, eval=FALSE} 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.