---
title: "Getting started with SynthesizerPlus"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Getting started with SynthesizerPlus}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7,
                      fig.height = 4.2, dpi = 72, out.width = "100%")
```

**SynthesizerPlus** generates synthetic data that mimic a real data set: the
synthetic records follow the same marginal distributions and the same
dependence structure, but they are new records. This vignette introduces the
modelling approach and the main functions. Other vignettes cover
[evaluation](evaluating-synthetic-data.html),
[time series](time-series.html),
[multivariate distributions](multivariate-distributions.html) and
[file formats](file-formats.html).

```{r}
library(SynthesizerPlus)
```

## 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

1. 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;
2. 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`);
3. 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

```{r}
fit <- fit_synthesizer(iris)
fit
syn <- generate(fit, n = 500, seed = 1)
head(syn)
```

`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:

```{r}
syn2 <- synthesize(mtcars, n = 100, seed = 2)
str(syn2[, 1:5])
```

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:

```{r}
out_csv <- file.path(tempdir(), "mtcars_synthetic.csv")
synthesize(mtcars, n = 100, seed = 2, file = out_csv)
file.exists(out_csv)
readLines(out_csv, n = 3)          # the first lines of the file
```

`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:

```{r}
syn_back <- read_data(out_csv, template = mtcars)
str(syn_back[, 1:5])
all.equal(syn_back, syn2, check.attributes = FALSE)   # the file holds the same data
```

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

```{r}
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)
str(generate(fit_p, n = 5, seed = 3))
```

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:

```{r}
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))
```

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:

```{r}
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)
```

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:

```{r}
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)
```

```{r, fig.height = 5}
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:

```{r}
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")
})
cor(mtcars$mpg, mtcars$wt, method = "spearman")
```

## 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}
r <- as.data.frame(diff(log(EuStockMarkets)))
fit_t <- fit_synthesizer(r, copula = "t")
fit_t
```

## Identifiers, model size and privacy

* Columns listed in `id_cols` are not modelled; synthetic records get fresh
  identifiers.
* By default the continuous marginals store all observed values.
  `knots = 200`, say, compresses each quantile function to 200 quantiles:
  the model gets smaller and no longer contains the original values.

```{r}
d <- data.frame(id = sprintf("R%04d", 1:1000), x = rexp(1000))
fit_k <- fit_synthesizer(d, id_cols = "id", knots = 200)
head(generate(fit_k, n = 3, seed = 1))
format(object.size(fit_k), units = "Kb")
```

## Multiple synthetic data sets

To reflect the uncertainty of the synthesis (for inference with multiple
synthetic data sets), draw several independent copies with `simulate()`:

```{r}
sims <- simulate(fit_mt, nsim = 5, seed = 1)
sapply(sims, function(s) coef(lm(mpg ~ wt, data = s))[["wt"]])
```

## Where next

* `compare_synthetic()` and the plotting functions assess the quality of the
  synthetic data — see `vignette("evaluating-synthetic-data")`.
* `generate_to_file()` writes very large synthetic data sets in chunks to
  CSV, Parquet, HDF5 and other formats — see `vignette("file-formats")`.