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 "North" "East" "North" "East" ...
#> $ employed : logi FALSE TRUE TRUE TRUE FALSE
#> $ survey_date: Date, format: "2024-10-28" "2024-10-28" ...
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.
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