fit_synthesizer() learns a generative model from a data set. The model
combines
a marginal model per column: an interpolated empirical quantile function for continuous variables (linear or monotone-spline interpolation between order statistics), and an empirical probability table for discrete, categorical, ordinal and logical variables. Interpolation means that synthetic values of a continuous variable lie between, rather than on, the observed values, always within the observed range. With small samples this slightly reduces the spread in the tails compared with the data;
a copula (Gaussian, Student-t or independence) that couples the columns. Discrete and categorical variables take part in the dependence structure: their pseudo-observations use the distributional transform, and their latent correlations are estimated by maximum likelihood (polyserial / polychoric), see Estimation details;
an optional missing-data model: for each column with missing values a missingness indicator is added to the copula, so that the rate of missingness and its (monotone) association with other variables are reproduced;
optional stratification (by), fitting separate models within
groups, with a pooled fallback for small groups.
Use generate() to draw synthetic data from the fitted model, or
synthesize() to fit and generate in one step.
fit_synthesizer(data, ...)# S3 method for data.frame
fit_synthesizer(
data,
copula = c("gaussian", "t", "independence"),
by = NULL,
min_stratum_size = 20L,
missing = c("joint", "independent", "drop"),
interpolation = c("linear", "spline"),
discrete_threshold = 20L,
knots = NULL,
id_cols = NULL,
pd_method = c("auto", "higham", "eigen", "none"),
nonlinear_target = NULL,
nominal_target = NULL,
nominal_decay = 1e-04,
nominal_maxit = 1000L,
nominal_spline_df = 3L,
t_profile = FALSE,
nonlinear_engine = c("cart", "gam", "auto"),
nonlinear_cv_folds = 5L,
nonlinear_cv_rule = c("one_se", "min"),
nonlinear_cv_score = c("crps", "mse"),
nonlinear_gam_k = 10L,
nonlinear_gam_interactions = FALSE,
nonlinear_local_k = 30L,
nonlinear_minbucket = 20L,
nonlinear_cp = 0.001,
nonlinear_maxdepth = 10L,
nonlinear_residual = c("leaf", "global", "none"),
nonlinear_calibration = c("quantile", "none"),
seed = NULL,
...
)
# S3 method for matrix
fit_synthesizer(data, ...)
# S3 method for default
fit_synthesizer(data, ...)
An object of class sp_synthesizer: a list whose main
components are pooled (the model fitted to all rows: marginals,
copula and, when requested, nonlinear or nominal), strata
(with by: stratum keys, values, counts, probs, per-stratum
models, NULL for strata smaller than min_stratum_size, and
conditional_fallback, a data frame of the strata whose own model has
no conditional model for the target), and settings.
With nonlinear_engine = "auto", pooled$nonlinear$selection (and the
same element of each stratum model) holds the selected candidate
(selected, engine, interactions), the reason ("cv", or a
fallback reason), rule, score ("crps" or "mse"), folds, a
scores table with per-candidate mean_score, se_score, mean_mse,
se_mse and any issue, and the per-fold matrices fold_scores and
fold_mse.
A data frame, matrix, atomic vector (numeric, integer,
logical, character, factor, Date, POSIXct) or ts object. Time
series are handled by fit_synthesizer.ts().
Further arguments passed to methods.
Dependence model: "gaussian" (default), "t" (heavier
joint tails) or "independence". See Estimation details.
Optional character vector of column names used to stratify the model. Each stratum gets its own marginals and copula; strata are sampled in proportion to their training size.
Strata with fewer rows than this use a pooled model fitted on all rows.
How missing values are treated: "joint" (default; model
missingness inside the copula), "independent" (reproduce the missing
rate of each column independently) or "drop" (do not reproduce
missingness: partially observed columns are generated without missing
values; a column with no observed training value has no distribution
to draw from and necessarily stays entirely missing). Under "joint",
missingness may depend on the other variables (missing at random); it
cannot depend on the unobserved value itself, which the data do not
identify (see Estimation details).
Interpolation of the empirical quantile function for
continuous variables: "linear" (default) or "spline" (monotone Hyman
spline, smoother densities).
Numeric columns with at most this many distinct values are treated as discrete (sampled from their observed support), provided their values are whole numbers or actually repeat (fewer distinct values than half the observations). A small sample of a continuous variable is therefore still modelled as continuous.
Optional maximum number of quantile knots stored per
continuous column. NULL (default) keeps every observation; a smaller
value shrinks the model and avoids storing the original values.
Optional character vector of identifier columns. They are not modelled; synthetic data get fresh sequential identifiers.
How to repair a correlation matrix that is not positive
definite: "auto"/"higham" (nearest correlation matrix), "eigen"
(eigenvalue flooring) or "none" (fail if not positive definite).
Optional name of ONE continuous numeric response to
generate using a hybrid copula--conditional model. NULL (default)
leaves the existing copula synthesizer unchanged. The copula generates
all predictors; CART or a penalised GAM predicts the response, local
residuals can be resampled, and the response is quantile-calibrated to its
empirical marginal. Requires rpart for CART or mgcv for GAM. Not supported
for a ts model or a response with discrete/categorical marginals.
The fitted tree and residual pools contain training-data information;
do not release the fitted model as if it were privacy-safe.
Optional name of one unordered categorical response
(factor, character or logical). When supplied, predictors are generated
using the copula and this response is drawn from a fitted multinomial
conditional model (nnet::multinom). Nominal categories are not
thresholded along a one-dimensional latent variable in the final draw.
Mutually exclusive with nonlinear_target. Experimental, maximum of 20
observed categories and at least 40 observed responses; strata without a
converged model fall back to copula synthesis.
Nonnegative weight-decay penalty passed to
nnet::multinom() (default 0.0001).
Maximum number of optimisation iterations for the nominal model (default 1000). A model that has not converged within this limit is not used.
Degrees of freedom of the natural cubic spline
basis used for each numeric predictor (with at least five distinct
values) in the nominal multinomial model, so category probabilities can
depend nonlinearly on numeric predictors. Default 3. The spline model
and the linear model are both fitted and the one with the smaller AIC
is used (recorded in nominal$spline_df and nominal$selection); 0
fits the linear model only.
If TRUE with copula = "t", add exploratory elapsed-time
stage measurements under fit$pooled$copula$profile and per-stratum
copulas. Profile fields are absent by default.
Conditional response model for nonlinear_target:
"cart" (default, unchanged), "gam" (penalised additive model
fitted by mgcv), or "auto" (experimental training-only K-fold
selection among CART, additive GAM and GAM with tensor interactions).
When "auto", the fitted model records the selection scores in
fit$pooled$nonlinear$selection (and separately per fitted stratum).
The option is ignored when nonlinear_target = NULL.
Integer number of folds for automatic engine selection, default 5 (2--10). Only non-missing target observations are used for validation; predictors and target are held out together.
Either "one_se" (default: prefer CART, then
additive GAM, then interaction GAM when their average fold score falls
within one standard error of the best) or "min" (select the minimum
average fold score). All choices are made using original training
data only; no benchmark holdout observations are consulted.
Criterion for automatic engine selection:
"crps" (default) scores the whole conditional predictive distribution
(fitted mean plus residual resampling, 50 draws per held-out row) by the
continuous ranked probability score, so a candidate that misrepresents
the conditional spread is penalised; "mse" scores only the conditional
mean by its squared prediction error. Both are reported in the
selection table.
Maximum basis dimension per GAM smooth (default 10, permitted 4--50), truncated when too few distinct covariate values occur.
Whether to include pairwise tensor
interactions among up to three eligible numeric smooth predictors in the
GAM (default FALSE). An additive GAM may miss interactions.
Number of nearest training fitted values for
local GAM residual resampling (default 30). For nonlinear_engine = "gam",
the existing residual mode "leaf" means this local fitted-mean
neighbourhood, since a GAM has no CART terminal leaves. This does not
necessarily reproduce variance changes orthogonal to the fitted mean.
Minimum observations per CART leaf (default 20).
Complexity penalty for CART splitting (default 0.001).
Maximum CART depth (default 10).
Experimental residual resampling mode for a
hybrid nonlinear_target: "leaf" (default: CART leaf groups, or
GAM local neighbourhoods defined by fitted conditional means), "global" (resample from all fitted residuals),
or "none" (use only fitted conditional means). With "none",
identical predictions are calibrated to the same value, so no artificial
within-leaf randomness is added.
Experimental marginal recalibration of the
hybrid response: "quantile" (default) maps each simulated value
through the empirical CDF of a reference sample of conditional draws,
stored at fitting time, and then through the observed response's
empirical quantile function. The result does not depend on how many
rows are generated per call, chunk or stratum. "none" uses raw
conditional draws, so its synthetic marginal need not equal the observed
response marginal and values can exceed its observed range.
Optional random seed for the fit. Fitting is stochastic
(random tie-breaking and the randomised distributional transform used
for the pseudo-observations, row subsampling for the t copula), so a
seed makes the fitted model reproducible when fitting and generation are
separate steps. The global random number generator state is restored
afterwards. (synthesize() has its own seed, covering both steps.)
Continuous response (nonlinear_target). The copula generates all
columns; the response is then replaced by a draw from a conditional model
given the generated predictors. Predictors are encoded from all training
rows (numeric columns median-imputed plus a missingness indicator;
categories, at most 30, plus a MISSING level); a category never observed
with the response is passed to the tree as unknown and routed by surrogate
splits. CART (rpart) gives a piecewise-constant conditional mean
\(\hat m(x)\) and, with residual mode "leaf", residuals are resampled
from the centred residuals of the same leaf. The GAM (mgcv::gam, REML,
cubic regression splines for numeric predictors with at least four
distinct values, select = TRUE so terms can be penalised to zero,
optional pairwise ti() interactions among at most three of them)
resamples centred residuals from the nonlinear_local_k training
observations whose fitted means are nearest. Calibration "quantile"
maps a raw draw \(y^*\) to \(F_Y^{-1}(\tilde G(y^*))\), where
\(F_Y^{-1}\) is the empirical quantile function of the observed response
and \(\tilde G\) the mid-point empirical CDF of 20000 raw conditional
draws for resampled training predictor rows, drawn once at fitting time
from a fixed random stream. The synthetic response marginal therefore
matches the observed one up to Monte Carlo error, independently of batch
size.
Automatic engine selection. CART, an additive GAM and a GAM with tensor
interactions are compared by K-fold cross-validation on the training rows
with an observed response. With nonlinear_cv_score = "crps" each
held-out response \(y\) is scored by the continuous ranked probability
score of \(M = 50\) conditional draws \(X_1,\dots,X_M\) from the
fold model, \(\frac1M\sum_i|X_i-y| - \frac1{2M^2}\sum_{i,j}|X_i-X_j|\);
with "mse" by \((y-\hat m(x))^2\). A candidate failing in any fold is
discarded. The "one_se" rule picks the simplest candidate (CART, then
additive GAM, then interaction GAM) whose mean score is within one
standard error of the best; the chosen engine is then refitted on all
rows.
Unordered categorical response (nominal_target). The response is
removed from the copula (its missingness indicator stays) and drawn from
a multinomial logit model (nnet::multinom, weight decay
nominal_decay), \(P(Y=k\mid x) \propto \exp(\beta_{k0} + \sum_j
f_{kj}(x_j))\), with \(f_{kj}\) a natural cubic spline with
nominal_spline_df degrees of freedom (knots at quantiles) for each
numeric predictor with at least five distinct values, and linear or
dummy terms otherwise. The model with linear terms only is also fitted,
and of the two converged models the one with the smaller AIC is used:
splines are kept when the data support the extra parameters, while for
(nearly) separated classes the better-determined linear logit is used.
Category probabilities have no interactions between predictors.
In a stratified fit, a stratum whose own model cannot fit the requested
conditional model (too few observed responses, fewer than two categories,
or no convergence) generates the target from its stratum copula; this is
reported by a warning and in strata$conditional_fallback. The fitted
conditional models contain training-data information and are not
disclosure-safe artefacts.
Latent correlations. For a pair of continuous variables the Gaussian copula uses the normal scores \(\Phi^{-1}(U)\) and the t copula inverts Kendall's tau, \(\rho = \sin(\pi\tau/2)\). As Kendall's tau is quadratic in the sample size, it is computed on a random subsample of 2000 rows; a pair that the subsample leaves with fewer than \(\min(n_{ij}, 500)\) jointly observed rows is re-estimated from up to 2000 of its own \(n_{ij}\) joint rows, so every pair with \(n_{ij} \ge 3\) is estimated whatever the subsample.
Discrete and categorical variables. The pseudo-observations of a discrete margin fill each category's probability interval with independent uniform noise (the distributional transform). Correlating these jittered scores would bias every dependence involving a discrete variable towards zero, most strongly for binary variables. Pairs with at least one discrete margin are therefore estimated by maximum likelihood under the Gaussian copula, with the thresholds fixed at the normal quantiles of the marginal cumulative probabilities: the polyserial correlation for a discrete-continuous pair (Olsson, Drasgow and Dorans, 1982) and the polychoric correlation for a discrete-discrete pair (Olsson, 1979), using the bivariate normal distribution function of Genz (2004). The t copula uses the same estimates for these pairs.
Order of unordered categories. A categorical variable occupies one
latent coordinate, so its categories need an order on it. Ordered factors
keep their level order and logical variables FALSE < TRUE. For unordered
factors and character columns the order is estimated from the data (a
form of optimal scaling): each numeric, ordered or logical column is
converted to normal scores \(\Phi^{-1}(r/(n+1))\) from its ranks
\(r\); the matrix \(M\) of category means of these scores (one row
per category) is centred with the category frequencies \(w\) as
weights, and the categories are sorted by their scores on the first
principal component of \(\mathrm{diag}(\sqrt{w})M\), the sign fixed so
that the first non-zero loading is positive, with ties broken by
decreasing frequency and then by label. Categories that are adjacent in
their association with the other columns are thus adjacent on the latent
coordinate, and renaming or reordering the categories does not change the
fitted model. Without informative reference columns, the declared level
order (factors) or sorted labels (character) are used. The same order is
used when generating, so the latent correlations are the maximum
likelihood estimates of the model actually fitted. One coordinate cannot
represent every association pattern of a nominal variable; see
nominal_target and stratification (by).
Missing data. All estimates use pairwise complete observations. The
rows in which two variables are both observed may be a selected subset
(e.g. y is observed mainly for low x), with a restricted range. For
continuous pairs the estimate is therefore the pairwise covariance of
the normal scores, scaled by each variable's standard deviation over all
of its observed values (for the t copula, the tau-based correlation
rescaled in the same way). This reproduces the relation among the
observed rows under such selection; with complete data it equals the
ordinary correlation. With missing = "joint", the missingness
indicators are binary variables of the copula, so missingness can depend
on the other variables (missing at random). The latent correlation
between a variable and its own indicator is fixed at 0: the data do not
identify it (the value is never seen when missing), and this choice keeps
the distribution of the observed values of each variable as in the
training data. Pairwise estimates based on different subsets of rows
need not form a coherent joint estimate, and the resulting matrix need not
be positive definite; it is then projected onto the set of valid
correlation matrices according to pd_method (Higham's nearest
correlation matrix by default). The projection guarantees a valid model,
not an optimal joint estimator.
Unidentified pairs. If two variables are observed together in fewer than
3 rows (e.g. questions from different survey modules), their correlation
cannot be estimated (with 2 rows it is trivially \(\pm 1\)); the same
applies when the pairwise estimate is not finite. Such
entries are not set to zero, which would assert independence and, through
the projection, distort the identified entries. They are filled by the
maximum-determinant positive-definite completion of the identified
entries (Dempster, 1972), which keeps every identified correlation and
makes the pair conditionally independent given the other variables. If
the identified entries cannot be completed, the unknown entries fall back
to 0. In either case the matrix may then still be modified by the
positive-definiteness repair. A warning reports how many pairs were
affected and how they were filled; copula_correlation(fit, what = "n_pair") gives the pairwise counts, and the fitted copula records the
method (completion: "none", "maxdet" or "zero_fallback") and the
largest change made by the repair (pd_adjustment), both shown by
summary(). Pair counts and identification always refer to the full
data.
t-copula degrees of freedom. With the correlation matrix held fixed,
df is chosen by profile likelihood over the grid
2, 3, 4, 5, 6, 8, 10, 15, 20, 30, 50, 100 using the complete rows (at most
5000), so it is a discrete grid estimate rather than a continuous maximum
likelihood estimate. If fewer than 10 complete rows are available, df
is not estimated and is fixed at 10; a warning is issued and the
fitted copula stores df_estimated = FALSE (see summary()).
Sklar, A. (1959). Fonctions de répartition à n dimensions et leurs marges. Publications de l'Institut de Statistique de l'Université de Paris, 8, 229--231.
Dempster, A. P. (1972). Covariance selection. Biometrics, 28(1), 157--175. tools:::Rd_expr_doi("10.2307/2528966")
Genz, A. (2004). Numerical computation of rectangular bivariate and trivariate normal and t probabilities. Statistics and Computing, 14(3), 251--260. tools:::Rd_expr_doi("10.1023/B:STCO.0000035304.20635.31")
Olsson, U. (1979). Maximum likelihood estimation of the polychoric correlation coefficient. Psychometrika, 44(4), 443--460. tools:::Rd_expr_doi("10.1007/BF02296207")
Olsson, U., Drasgow, F. and Dorans, N. J. (1982). The polyserial correlation coefficient. Psychometrika, 47(3), 337--347. tools:::Rd_expr_doi("10.1007/BF02294164")
Rüschendorf, L. (2009). On the distributional transform, Sklar's theorem, and the empirical copula process. Journal of Statistical Planning and Inference, 139(11), 3921--3927. tools:::Rd_expr_doi("10.1016/j.jspi.2009.05.030")
Higham, N. J. (2002). Computing the nearest correlation matrix -- a problem from finance. IMA Journal of Numerical Analysis, 22(3), 329--343. tools:::Rd_expr_doi("10.1093/imanum/22.3.329")
generate(), synthesize(), compare_synthetic()
fit <- fit_synthesizer(iris)
fit
syn <- generate(fit, n = 300, seed = 1)
head(syn)
# Stratified model and missing-data model
fit_aq <- fit_synthesizer(airquality, by = "Month")
colSums(is.na(generate(fit_aq, seed = 2)))
Run the code above in your browser using DataLab