The model in one paragraph
By Sklar’s theorem every multivariate distribution can be written as
a copula (the dependence structure) applied to the
marginal distributions. SynthesizerPlus estimates
- one marginal model per column — an
interpolated empirical quantile function for continuous
variables (linear or monotone-spline interpolation between order
statistics, so synthetic values fall between observed ones and within
the observed range), and the empirical probability table for discrete,
categorical, ordinal and logical variables;
- a copula — Gaussian (default), Student-t
or independence — whose correlations are estimated from the normal
scores of the data. Discrete and categorical columns are part of the
dependence structure instead of being sampled independently; their
correlations are estimated by maximum likelihood (polyserial and
polychoric correlations, see
?fit_synthesizer);
- a missing-data model — each column with missing
values gets a missingness indicator that is also part of the copula, so
both the missingness rate and its (monotone) relation to other variables
are reproduced.
Generation reverses the steps: draw from the copula, push each
uniform through the corresponding quantile function, and restore the
original column classes (factor levels, Date,
POSIXct with time zone, integer, …).
Fit once, generate many times
fit <- fit_synthesizer(iris)
fit
#> <sp_synthesizer>
#> trained on : 150 rows, 5 modelled column(s)
#> marginals : 1 categorical, 4 continuous
#> copula : gaussian
#> missing : joint
syn <- generate(fit, n = 500, seed = 1)
head(syn)
#> Sepal.Length Sepal.Width Petal.Length Petal.Width Species
#> 1 5.1 3.056437 4.208853 1.500000 versicolor
#> 2 6.0 2.920090 5.100000 1.400000 versicolor
#> 3 5.0 2.600000 1.500000 1.001553 versicolor
#> 4 7.2 3.000000 5.961895 1.900000 virginica
#> 5 6.1 3.400000 4.500000 1.506789 versicolor
#> 6 5.0 3.800000 1.300000 0.200000 setosa
generate() can be called as often as needed;
n may be larger or smaller than the training data. The
seed argument makes the output reproducible
without changing the random-number state of your session.
Fitting is itself slightly random (ties among the pseudo-observations
are broken at random, and discrete variables use a randomised
transform), so for a fully reproducible workflow give the fit a seed
too: fit_synthesizer(mtcars, seed = 1).
For a quick one-off, synthesize() fits and generates in
one call:
syn2 <- synthesize(mtcars, n = 100, seed = 2)
str(syn2[, 1:5])
#> 'data.frame': 100 obs. of 5 variables:
#> $ mpg : num 17.9 15.1 23.8 12.2 14.9 ...
#> $ cyl : num 8 8 4 8 8 8 6 8 6 4 ...
#> $ disp: num 269 456 120 455 393 ...
#> $ hp : num 149 175 63 224 230 245 110 208 170 96 ...
#> $ drat: num 3.83 3.04 3.7 3.12 2.76 ...
With a file argument, the synthetic data are also
written to disk (and returned invisibly). The format follows the file
extension; here a CSV file in a temporary folder:
out_csv <- file.path(tempdir(), "mtcars_synthetic.csv")
synthesize(mtcars, n = 100, seed = 2, file = out_csv)
file.exists(out_csv)
#> [1] TRUE
readLines(out_csv, n = 3) # the first lines of the file
#> [1] "\"mpg\",\"cyl\",\"disp\",\"hp\",\"drat\",\"wt\",\"qsec\",\"vs\",\"am\",\"gear\",\"carb\""
#> [2] "17.8873149337211,8,269.003259089779,149,3.82749154182601,3.20574434430353,15.4601901371236,0,1,4,4"
#> [3] "15.0848824562055,8,456.467090919786,175,3.03626534302339,4.95925000836009,17.9902637328502,0,0,3,2"
read_data() reads the file back. With
template = the column classes of the original data are
restored, which matters for formats such as CSV that do not store
them:
syn_back <- read_data(out_csv, template = mtcars)
str(syn_back[, 1:5])
#> 'data.frame': 100 obs. of 5 variables:
#> $ mpg : num 17.9 15.1 23.8 12.2 14.9 ...
#> $ cyl : num 8 8 4 8 8 8 6 8 6 4 ...
#> $ disp: num 269 456 120 455 393 ...
#> $ hp : num 149 175 63 224 230 245 110 208 170 96 ...
#> $ drat: num 3.83 3.04 3.7 3.12 2.76 ...
all.equal(syn_back, syn2, check.attributes = FALSE) # the file holds the same data
#> [1] TRUE
The same works for any supported format,
e.g. file = "synthetic.parquet",
"synthetic.xlsx" or "synthetic.h5" (see
supported_formats() and
vignette("file-formats")). A fitted model can also write a
new file whenever needed, with
write_data(generate(fit, n = 1000), "synthetic.rds").
Supported column types
set.seed(10)
n <- 200
people <- data.frame(
age = as.integer(round(runif(n, 18, 80))),
income = round(rlnorm(n, 10, 0.6)),
sector = factor(sample(c("agri", "industry", "services"), n, TRUE,
prob = c(0.1, 0.3, 0.6))),
education = factor(sample(c("primary", "secondary", "tertiary"), n, TRUE),
levels = c("primary", "secondary", "tertiary"), ordered = TRUE),
region = sample(c("North", "South", "East", "West"), n, TRUE),
employed = runif(n) < 0.7,
survey_date = as.Date("2024-01-01") + sample(0:365, n, TRUE),
stringsAsFactors = FALSE
)
people$income[sample(n, 20)] <- NA
fit_p <- fit_synthesizer(people)
summary(fit_p)
#> Synthesizer trained on 200 rows
#> Copula: gaussian | missing data: joint
#>
#> column class model support missing
#> age integer continuous [19, 79] 0.0%
#> income numeric continuous [5105, 94650] 10.0%
#> sector factor categorical 3 values 0.0%
#> education ordered categorical 3 values 0.0%
#> region character categorical 4 values 0.0%
#> employed logical categorical 2 values 0.0%
#> survey_date Date continuous [19720, 20090] 0.0%
str(generate(fit_p, n = 5, seed = 3))
#> 'data.frame': 5 obs. of 7 variables:
#> $ age : int 29 40 52 25 51
#> $ income : num 23010 23894 41967 10228 44065
#> $ sector : Factor w/ 3 levels "agri","industry",..: 2 2 2 3 3
#> $ education : Ord.factor w/ 3 levels "primary"<"secondary"<..: 1 1 1 3 2
#> $ region : chr "South" "West" "South" "West" ...
#> $ employed : logi TRUE TRUE TRUE TRUE TRUE
#> $ survey_date: Date, format: "2024-11-03" "2024-11-11" ...
Numeric columns with at most discrete_threshold (default
20) distinct values are modelled as discrete, so a variable such as the
number of cylinders only takes observed values. Continuous values are
always within the observed range.
Order of unordered categories
In the copula each categorical column occupies one latent coordinate,
on which its categories are arranged in an order. Ordered factors keep
their level order. For unordered factors and character columns the order
is estimated from the data: categories are sorted by
the first principal component of their mean normal scores on the
numeric, ordered and logical columns, so categories that behave alike
with respect to the other columns are neighbours on the latent
coordinate (details in ?fit_synthesizer). As a consequence
the fitted model does not depend on the names or the level order of the
categories:
fit_iris <- fit_synthesizer(iris, seed = 1)
fit_iris$pooled$marginals$Species$values # estimated order
#> [1] "setosa" "versicolor" "virginica"
iris2 <- iris
iris2$Species <- factor(iris2$Species,
levels = c("virginica", "setosa", "versicolor"))
fit_iris2 <- fit_synthesizer(iris2, seed = 1)
all.equal(fit_iris$pooled$copula$R, fit_iris2$pooled$copula$R)
#> [1] TRUE
(Versions before 0.4.0 used the order of the labels, so synthetic
data involving unordered categorical columns differ from those of
earlier versions.) One latent coordinate cannot represent every
association of a nominal variable; when it matters, stratify on it
(by, below) or, for one nominal response, use
nominal_target
(vignette("nominal-synthesis")).
Missing values
The default missing = "joint" reproduces the share of
missing values and how missingness depends on the other variables, as
long as that dependence is monotone (missing more often for
larger, or for smaller, values of another variable). Here income is
missing more often for older people:
set.seed(11)
age <- round(runif(2000, 18, 80))
income <- exp(7 + 0.01 * age + rnorm(2000, sd = 0.5))
income[runif(2000) < plogis(-4 + 0.05 * age)] <- NA
svy <- data.frame(age, income)
syn_svy <- generate(fit_synthesizer(svy, seed = 1), n = 20000, seed = 1)
bands <- function(a) cut(a, c(17, 30, 45, 60, 80))
rbind(real = tapply(is.na(svy$income), bands(svy$age), mean),
synthetic = tapply(is.na(syn_svy$income), bands(syn_svy$age), mean))
#> (17,30] (30,45] (45,60] (60,80]
#> real 0.05808081 0.09832636 0.2205882 0.3938462
#> synthetic 0.05137940 0.13036338 0.2030924 0.3854819
A single copula cannot represent a non-monotone pattern. In
airquality, ozone is missing far more often in June than in
the months before and after it; the pooled model can only spread that
excess into a trend over the months. Stratifying by month reproduces the
pattern:
aq_pooled <- generate(fit_synthesizer(airquality, seed = 1), n = 5000, seed = 4)
aq_strat <- generate(fit_synthesizer(airquality, by = "Month", seed = 1), n = 5000, seed = 4)
miss_by_month <- function(d) tapply(is.na(d$Ozone), d$Month, mean)
round(rbind(real = miss_by_month(airquality), pooled = miss_by_month(aq_pooled),
`by = "Month"` = miss_by_month(aq_strat)), 2)
#> 5 6 7 8 9
#> real 0.16 0.70 0.16 0.16 0.03
#> pooled 0.42 0.31 0.24 0.17 0.10
#> by = "Month" 0.16 0.71 0.15 0.16 0.02
Missingness that depends on the unobserved value itself (e.g. high
incomes are not reported because they are high) cannot be
learned from the data, and is not reproduced.
Use missing = "independent" to reproduce only the rates,
or "drop" to generate complete data.
Stratified models
A single Gaussian copula cannot represent every kind of dependence.
Strong, non-monotone relations with a categorical variable — such as the
species clusters in iris — are attenuated. Fitting separate
models within strata solves this:
syn_pooled <- generate(fit_synthesizer(iris), seed = 1)
syn_strat <- generate(fit_synthesizer(iris, by = "Species"), seed = 1)
quality <- rbind(
pooled = summary(compare_synthetic(iris, syn_pooled, seed = 1))$value[c(3, 5)],
stratified = summary(compare_synthetic(iris, syn_strat, seed = 1))$value[c(3, 5)]
)
colnames(quality) <- c("mean |assoc. diff|", "discriminator AUC")
round(quality, 3)
#> mean |assoc. diff| discriminator AUC
#> pooled 0.077 0.652
#> stratified 0.032 0.409
plot_pairs(iris, syn_strat, vars = c("Sepal.Length", "Petal.Length", "Petal.Width"))
Strata smaller than min_stratum_size (default 20) fall
back to a pooled model, which avoids over-fitting (and reproducing) tiny
groups. Stratum sizes in the synthetic data are drawn from a multinomial
distribution, or set proportionally with
generate(..., stratum_sizes = "proportional").
Controlling the dependence
dependence shrinks the copula correlation matrix towards
the identity without changing any marginal distribution — useful for
sensitivity analyses or to reduce disclosure risk:
fit_mt <- fit_synthesizer(mtcars)
sapply(c(0, 0.5, 1), function(l) {
s <- generate(fit_mt, n = 2000, dependence = l, seed = 5)
cor(s$mpg, s$wt, method = "spearman")
})
#> [1] 0.02856503 -0.46944277 -0.91193002
cor(mtcars$mpg, mtcars$wt, method = "spearman")
#> [1] -0.886422
Heavy tails: the t copula
For data with joint extremes (e.g. financial returns) the Gaussian
copula underestimates the probability of simultaneous large values. The
t copula estimates its degrees of freedom by profile
likelihood:
r <- as.data.frame(diff(log(EuStockMarkets)))
fit_t <- fit_synthesizer(r, copula = "t")
fit_t
#> <sp_synthesizer>
#> trained on : 1859 rows, 4 modelled column(s)
#> marginals : 4 continuous
#> copula : t (df = 8)
#> missing : joint
The degrees of freedom are chosen on a grid by maximising the copula
likelihood with the correlation matrix held fixed. To see where the
fitting time goes on large data, t_profile = TRUE records
the elapsed time of each stage (correlation estimation, latent
correlations of discrete columns, matrix repair and the
degrees-of-freedom search):
fit_tp <- fit_synthesizer(r, copula = "t", t_profile = TRUE)
fit_tp$pooled$copula$profile
#> kendall_seconds discrete_latent_seconds
#> 0.295 0.000
#> correlation_repair_seconds df_grid_seconds
#> 0.001 0.033
#> other_seconds total_seconds
#> 0.000 0.329