# Probabilistic Forecast Reconciliation with Kullback-Leibler Divergence
Regularization

Guanyu Zhang  AI Lab
Lenovo Research
Beijing, China
<EMAIL_ADDRESS>Feng Li2  School of Statistics and Mathematics
Central University of Finance and Economics
Beijing, China (Corresponding author)
<EMAIL_ADDRESS>Yanfei Kang Presented in The 2023 IEEE International
Conference on Data Mining (ICDM): Artificial Intelligence for Time Series
Analysis Workhsop (AI4TS). School of Economics and Management
Beihang University
Beijing, China
<EMAIL_ADDRESS>
###### Abstract

As the popularity of hierarchical point forecast reconciliation methods
increases, there is a growing interest in probabilistic forecast
reconciliation. Many studies have utilized machine learning or deep learning
techniques to implement probabilistic forecasting reconciliation and have made
notable progress. However, these methods treat the reconciliation step as a
fixed and hard post-processing step, leading to a trade-off between accuracy
and coherency. In this paper, we propose a new approach for probabilistic
forecast reconciliation. Unlike existing approaches, our proposed approach
fuses the prediction step and reconciliation step into a deep learning
framework, making the reconciliation step more flexible and soft by
introducing the Kullback-Leibler divergence regularization term into the loss
function. The approach is evaluated using three hierarchical time series
datasets, which shows the advantages of our approach over other probabilistic
forecast reconciliation methods.

###### Index Terms:

Forecasting, Hierarchical time series, Probabilistic forecast reconciliation

##  I Introduction

Multivariate time series in many practical applications are often arranged in
hierarchies, i.e., time series at upper levels of the hierarchy are aggregates
of those at lower levels, which can be represented by a tree diagram. When
tackling the hierarchical time series forecasting problem, for
interpretability and better helping downstream tasks to make decisions, we
require that the forecasts meet the coherency, i.e. forecasts of aggregated
time series are the sum of forecasts of the corresponding disaggregated time
series. To achieve coherency, a rich literature has emerged on forecast
reconciliation, such as bottom-up (BU), top-down (TD), and MinT [1].

Currently, most reconciliation methods are based on point forecasts, and there
is little research on probabilistic forecast reconciliation. However,
probabilistic forecast is very important in practical applications and is
helpful for risk management and decision-making. Especially in the power,
wind, and other new energy industries, data is very unstable, and the
probabilistic forecast is more practical than point forecast, which helps
staff allocate energy rationally and maximize the efficiency of energy use
[2]. Probabilistic forecast reconciliation is very challenging because it does
not just reconcile the mean prediction values but reconciles multiple related
predictive distributions. Most of the probabilistic forecast reconciliation
methods proposed in existing research are realized by sampling, reordering,
specifying hierarchical dependencies, and choosing reconciliation weights by
optimizing a scoring rule [3]. These approaches have two drawbacks: First, the
whole process contains too many artificial assumptions and is not totally
data-driven. Second, the prediction and reconciliation steps are regarded as
two separate independent processes, which is not conducive to parameter tuning
and model updating under data mode conversion, and the reconciliation step
cannot be adjusted accordingly for different prediction methods. Although
Rangapuram _et al._ proposes an end-to-end framework to solve this problem,
the reconciliation step in this framework is actually a hard mapping step and
not flexible and soft enough [4].

Recently, deep learning has shown great potential in the field of time series
forecasting. Deep neural networks can learn complex data representations,
which reduces the need for manual feature engineering and model design.
Outputs and loss functions are flexible, allowing for the customization of
different aims, such as point forecasting, probabilistic forecasting, and
other special requirements. Some novel deep learning architectures developed
for time series forecasting include DeepAR [5], N-BEATS [6], and Temporal
Fusion Transformer [7]. However, the application of deep learning in the field
of probabilistic forecast reconciliation is still lacking.

Inspired by the end-to-end method proposed by Rangapuram _et al._ [4], we
propose a new probabilistic forecast reconciliation method based on deep
learning to solve the above problems. First, we use a deep learning model to
predict multivariate time series and introduce the Kullback–Leibler divergence
(KLD) regularization term into the loss function, which measures the
“distance” between multiple distributions, i.e. the distance between the
parent node distribution and the sum of the child node distributions. In this
way, the prediction and reconciliation steps are fused into a deep learning
framework, instead of treating reconciliation as a post-processing step, the
KLD regularization term will help the model achieve approximate coherency of
the predictive distribution at training time. We set a coefficient for the KLD
regularization term to control the degree of reconciliation, which can be
tuned as a hyperparameter based on different objectives. For example, we can
select the optimal coefficient based on the accuracy of a validation set to
avoid sacrificing accuracy for the sake of coherency. Second, the multivariate
predictive distribution given by the trained model does not meet the complete
coherency, so the idea of bottom-up can be used to transform the incompletely
coherent multivariate distribution into a strictly coherent multivariate
distribution.

Compared with the existing probabilistic forecast reconciliation methods, our
proposed method has the following advantages: First, the prediction step and
reconciliation step are integrated into a deep learning framework, which is
conducive to direct optimization of the target, avoiding error accumulation
and facilitating model updating. Second, MinT and other traditional
reconciliation methods have limited consideration of the correlation between
multiple time series. However, the deep learning model is more able to capture
relevant information. Third, the coefficient in the KLD regularization term is
adjustable, which is conducive to adjusting the degree of reconciliation.
Sometimes, there is a trade-off between coherency and accuracy, and we can set
a reasonable degree of reconciliation to meet different research goals.

The rest of the paper is organized as follows. In Section II, we introduce the
notation and review the existing common forecast reconciliation methods. We
then demonstrate our new probabilistic forecast reconciliation method in
Section III. In Section IV, we conduct experiments on three data sets to
compare our approach with existing common probabilistic forecast
reconciliation methods. We conclude the paper and discuss future work in
Section V. The code for reproducing the results is available at
https://github.com/guanyu0316/Probabilistic-Forecast-Reconciliation-with-DL.

##  II Preliminaries

###  II-A Key Terms and Definitions

---
Figure 1:  A hierarchy example.

Let $\bm{y}_{\bm{t}}$ denote an n-vector containing observations at time $t$
of all nodes in a hierarchical structure. Let $\bm{b}_{\bm{t}}$ be an m-vector
with the observations at time $t$ of bottom nodes. In the field of
hierarchical forecasting, an important matrix is the summing matrix $\bm{S}$,
which is able to map the bottom time series to $\bm{y}_{\bm{t}}$, i.e.
$\bm{y}_{\bm{t}}=\bm{S}\bm{b}_{\bm{t}}$. Once the hierarchical structure is
determined, $\bm{S}$ is determined. Taking the hierarchical structure shown in
Figure 1 as an example, a company can organize its five SKUs into two
categories, then

$\left[\begin{array}[]{c}y_{h}\\\ y_{A,h}\\\ y_{B,h}\\\ y_{AA,h}\\\
y_{AB,h}\\\ y_{AC,h}\\\ y_{BA,h}\\\
y_{BB,h}\end{array}\right]=\left[\begin{array}[]{lllll}1&1&1&1&1\\\
1&1&1&0&0\\\ 0&0&0&1&1\\\ 1&0&0&0&0\\\ 0&1&0&0&0\\\ 0&0&1&0&0\\\ 0&0&0&1&0\\\
0&0&0&0&1\end{array}\right]\left[\begin{array}[]{l}y_{AA,h}\\\ y_{AB,h}\\\
y_{AC,h}\\\ y_{BA,h}\\\ y_{BB,h}\end{array}\right].$

The goal of point forecast is to estimate the $h$ period-ahead conditional
expectation

${E}\left[\bm{y}_{T+h}\mid\bm{y}_{1},\ldots,\bm{y}_{T}\right].$

We can compute forecasts for all series at all levels, which we call base
forecasts $\hat{\bm{y}}_{h}$. However, these forecasts may not meet coherency,
i.e. $\hat{\bm{y}}_{h}=\bm{S}\hat{\bm{b}}_{h}$, so we need to transform the
base forecasts to meet the coherency, and this process is called
reconciliation. $\tilde{\bm{y}}_{h}$ represents forecast after being
reconciled, which must meet coherency.

While the goal of probabilistic forecast is to estimate, for each node in the
hierarchy, the $h$ period-ahead conditional cumulative distribution function

$F_{i,T+h}\left(y\mid\bm{y}_{1},\ldots,\bm{y}_{T}\right)={P}\left(y_{i,T+h}\leq
y\mid\bm{y}_{1},\ldots,\bm{y}_{T}\right).$

In the field of probabilistic forecast, coherency is defined as follows: Let
$X_{i}\sim\hat{F}_{i}\ for\ i=1,...,n$ and let $i_{1},\ldots,i_{n_{k}}$ denote
the $n_{k}$ children of series $i$. The forecasts $\hat{F}_{i}$ are
probabilistic coherent if $X_{i}\stackrel{{\scriptstyle
d}}{{=}}X_{i_{1}}+\cdots+X_{i_{n_{k}}}\ for\ i=1,...,r$, where
$\stackrel{{\scriptstyle d}}{{=}}$ denoted equality in distribution [8]. In
other words, the forecast is coherent if the predictive distribution of each
aggregate series is equal to the distribution of the sum of the children
series. The process of generating coherent probabilistic forecasts using a
specific method is referred to as probabilistic forecast reconciliation, which
is the focus of this study.

###  II-B Benchmark Models

In this subsection, we will introduce the benchmark models used for
comparison. For illustrative purposes, we take the hierarchical structure
shown in Figure 1 as an example. In order to clearly introduce the existing
probabilistic forecast reconciliation methods, let’s first introduce point
forecast reconciliation methods. Common point forecast reconciliation methods
include bottom-up, top-down, and MinT. The essence of these reconciliation
methods is the linear transformation

$\widetilde{\bm{y}}_{h}=\bm{S}\bm{P}\widehat{\bm{y}}_{h}.$

The $\bm{P}$ matrix can map the base forecast $\widehat{\bm{y}}_{h}$ to the
forecast on the bottom level, then $\bm{S}$ can restore this bottom level
forecast to the forecast vector of all nodes. Once the hierarchical structure
is determined, $\bm{S}$ is determined, while $P$ depends on different
reconciliation methods. For example, in the bottom-up approach

$\bm{P}=\left[\begin{array}[]{llllllll}0&0&0&1&0&0&0&0\\\ 0&0&0&0&1&0&0&0\\\
0&0&0&0&0&1&0&0\\\ 0&0&0&0&0&0&1&0\\\ 0&0&0&0&0&0&0&1\end{array}\right],$

i.e., the reconciled forecast is obtained by aggregating the bottom-level
forecast. In the top-down method

$\bm{P}=\left[\begin{array}[]{llllllll}p_{1}&0&0&0&0&0&0&0\\\
p_{2}&0&0&0&0&0&0&0\\\ p_{3}&0&0&0&0&0&0&0\\\ p_{4}&0&0&0&0&0&0&0\\\
p_{5}&0&0&0&0&0&0&0\end{array}\right],$

i.e. the forecast of the top node is disassembled to the bottom level
according to a certain proportion vector, and the bottom level forecast will
then be aggregated into the forecast of all nodes using $\bm{S}$.
Wickramasuriya _et al._ proved that when $\bm{P}$ takes the next formula, the
trace of the variance-covariance matrix $\bm{V}_{h}$ of the reconciled
prediction error is the smallest [9]. Therefore, the optimal reconciliation
method is also called MinT (Minimize Trace).

$\bm{P}=\left(\bm{S}^{\prime}\bm{W}_{h}^{-1}\bm{S}\right)^{-1}\bm{S}^{\prime}\bm{W}_{h}^{-1},\\\
$
$\bm{V}_{h}=\operatorname{Var}\left[\bm{y}_{T+h}-\widetilde{\bm{y}}_{h}\right]=\bm{S}\bm{P}\bm{W}_{h}\bm{P}^{\prime}\bm{S}^{\prime},$

where
$\bm{W}_{h}=\operatorname{Var}\left[\left(\bm{y}_{T+h}-\widehat{\bm{y}}_{h}\right)\right]$
is the variance-covariance matrix of base forecast error, which is unknown and
difficult to be estimated accurately. In practice, it is common to make
simplifying assumptions about $\bm{W}_{h}$. A commonly used assumption is
$\bm{W}_{h}=k_{h}\bm{\Lambda},\bm{\Lambda}=\operatorname{diag}(\mathbf{S1})$,
i.e. the base forecast error variance of a node is proportional to the number
of child nodes it contains, which is called STRUCT assumption.

Compared with the research on point forecast reconciliation, few studies on
probabilistic forecast reconciliation exist. Jeon _et al._ proposed a
relatively complete probabilistic forecast reconciliation approach [3].
Specifying the hierarchical sample arrangement form, the point forecast
reconciliation methods can be applied to probabilistic forecast
reconciliation. Suppose we can obtain a sample from the predictive joint
distribution:

$\widehat{\bm{Y}}=\left(\widehat{\bm{y}}_{1}^{t+h\mid
t},\ldots,\widehat{\bm{y}}_{N}^{t+h\mid t}\right).$

Every column of $\widehat{\bm{Y}}$ is a sample vector from the predictive
joint distribution. We can apply the formula of point forecast reconciliation
to probabilistic forecast as follows. Point forecast reconciliation is to
reconcile a forecast vector, while probabilistic forecast reconciliation is to
reconcile a matrix:

$\bm{\widetilde{Y}=SP\widehat{Y}s}.$

Jeon _et al._ proposed three ways to construct joint distribution samples,
which are abbreviated as stack, rank and random [3]. The stack method refers
to that we sample from the predictive distribution of all nodes separately,
and then the samples are concatenated directly:

$\widehat{\bm{Y}}^{S}=\left[\begin{array}[]{c}\widehat{\mathbf{Z}}_{1}\\\
\widehat{\mathbf{Z}}_{2}\\\ \vdots\\\ \widehat{\bm{Z}}_{L}\end{array}\right],$

where $\widehat{\bm{Z}}_{l}$ is vector of $N$ samples from predictive
distribution of node $l$. The essence of this approach is to assume that the
time series of all nodes are independent. The rank method means that the
values of column $i$ of the constructed joint distribution sample
$\widehat{\bm{Y}}^{S}$ are all the $(i/N)$ quantiles of the predictive
distributions of all nodes. The rank method assumes that the series of all
nodes or all levels are highly correlated. The random method refers to the
random rearrangement of the samples drawn from each node and then
concatenating them together like the stack method. The implicit assumption of
the random method is that the series of all nodes are weakly correlated.

After obtaining the joint distribution samples through different permutation
methods, we can use MinT or BU methods to reconcile the predictive joint
distribution samples just like point forecast reconciliation, which can be
understood as treating each column of $\widehat{\bm{Y}}$ as a point forecast
sequence and reconciliating them separately.

Currently, both the point forecast reconciliation method and the probabilistic
forecast reconciliation method can be regarded as a two-step process, i.e.,
the first step is to obtain the base forecast, and the second is to reconcile
the base forecast. This idea is not conducive to tuning parameters and model
adjustment under data mode transformation. In addition, the above
probabilistic forecasting reconciliation method implies too many artificial
assumptions of hierarchical dependence and is not flexible enough.

##  III Probabilistic Forecast Reconciliation with Kullback-Leibler Divergence
Regularization

In this section, we present a probabilistic forecast reconciliation approach
based on deep learning with Kullback-Leibler divergence regularization, with
the objective of overcoming some of the aforementioned limitations of the
existing probabilistic forecast reconciliation method. Our approach can
realize end-to-end training and ensure the coherency of probabilistic forecast
results. In our approach, the deep learning method is not specified, all
global probabilistic forecasting models for multivariate time series based on
deep learning can be used in our approach.

###  III-A Proposed approach

Taking the hierarchical structure in Figure 1 as an example, there are eight
time series in total, and the value of the time series $i$ at time $t$ is
defined as $y_{i,t}$, $[1,t_{0}-1]$ is the history period or conditioning
range, $[t_{0},T]$ is the prediction range. Suppose we use a deep learning
model for producing probabilistic forecasts, which can model multiple time
series simultaneously. In this deep learning model, the input to the network
is history information and available covariates, and the output is the
parameter of the predictive probabilistic distribution. Given a data set of
time series and associated covariates, the parameters of this deep learning
model can be learned by minimizing the loss:

$L=\sum_{i=1}^{N}\sum_{t=t_{0}}^{T}l\left(y_{i,t},\theta_{i,t}\right).$

Where $\theta_{i,t}$ is the parameter of the predictive probabilistic
distribution, and N is 8 in the example of Figure 1. $l$ can be the negative
log-likelihood function. In order to meet the coherency of multiple series
probabilistic forecasts, we add a KLD regularization term to the loss
function:

$L=\sum_{i=1}^{N}\sum_{t=t_{0}}^{T}l\left(y_{i,t},\theta_{i,t}\right)+KL_{reg}.$

Taking the hierarchical structure in Figure 1 as an example, $KL_{reg}$ can be
represented by the following equation:

$\begin{gathered}KL_{reg}=\lambda_{1}D_{KL}\left(F_{Y_{.t}}\|F_{\bm{Y}_{A.t}+\bm{Y}_{B.t}}\right)\\\
+\lambda_{2}D_{KL}\left(F_{\bm{Y}_{A.t}}\|F_{\bm{Y}_{AA.t}+\bm{Y}_{AB.t}+\bm{Y}_{AC.t}}\right)\\\
+\lambda_{3}D_{KL}\left(F_{\bm{Y}_{B.t}}\|F_{\bm{Y}_{BA.t}+\bm{Y}_{BB.t}}\right),\end{gathered}$

where $\bm{Y}_{i,t}\sim F_{\bm{Y}_{i,t}}$. Multiple KLD terms are added, and
each term is the KLD between the distribution of a parent node and the
distribution of the sum of its child nodes. In mathematical statistics, the
KLD denoted $D_{KL}(P\|Q)$, measuring how one probabilistic distribution $P$
differs from a second, reference probabilistic distribution $Q$. The specific
definition of KLD is: for discrete probabilistic distributions $P$ and $Q$
defined on the same sample space, the KLD from $Q$ to $P$ is defined to be

$D_{KL}(P\|Q)=\sum_{x\in\mathcal{X}}P(x)\log\left(\frac{P(x)}{Q(x)}\right).$

For distributions $P$ and $Q$ of a continuous random variable, KLD is defined
to be the integral:

$D_{KL}(P\|Q)=\int_{-\infty}^{\infty}p(x)\log\left(\frac{p(x)}{q(x)}\right)dx.$

So KLD is not symmetric in the two distributions. When calculating KLD, we
swap positions between the two distributions and calculate the mean of the two
KLDs. Taking the first item in $KL_{reg}$ as an example, let’s amend it to

$\begin{gathered}D_{symmetric}=\\\
[D_{KL}\left(F_{Y,t}\|F_{Y_{A,t}+Y_{B,t}}\right)+D_{KL}\left(F_{Y_{A,t}+Y_{B,t}}\|F_{Y_{,t}}\right)]/2.\end{gathered}$

The KLD regularization term will help the model achieve approximate coherency
of the predictive probabilistic distribution at training time. The point is
how to guarantee the differentiability of the regularization term in the
neural network. KLD is usually calculated using samples from the two
distributions. In order to ensure the differentiability of the sampling step,
we use the reparameterization technique to make the sampling step independent
of the model, so that the sample from the two distributions is only related to
the distribution parameters, i.e. only related to parameters in the whole
neural network, which ensures the differentiability of this step. For example,
suppose the predictive distribution is Gaussian distribution. Before training,
we draw samples from the standard multivariate normal distribution, i.e.
$\bm{z\sim N(0,I)}$. When we want to sample from the predictive distribution,
we do the reparameterization:

$\bm{y_{t}=\mu_{t}+\Sigma_{t}^{1/2}z},$

then the sampling step is outside of the network, which ensures the
differentiability. When calculating the distribution of the sum of child
nodes, we assume that the time series of each child node is independent of
each other, but this does not mean that our method does not consider the
correlation information between each time series, because the deep learning
method itself has taken into account the correlation, which is the advantage
of deep learning.

The $\lambda$ coefficient in the KLD regularization term is used to control
the degree of reconciliation. In practice, coherency usually improves accuracy
[9], and subsequent experiments in this paper confirm this. However, for deep
learning models, adding regularization terms to the loss function may affect
the accuracy, and there is a trade-off between accuracy and coherency, so we
set $\lambda$ as a hyperparameter. $\lambda$ close to zero does not enforce
coherency, while larger $\lambda$ make predictions more coherent. Depending on
the practical application, we can choose a soft penalty version or a hard
constraint version. In our experiment, we set the same $\lambda$ for all items
in $KL_{reg}$. However, different $\lambda$ can be set if having sufficient
GPU computing power.

In the actual training process, it is necessary to specify the batch size as
an integer multiple of the total number of nodes in the hierarchical
structure, i.e. each training batch must exactly contain the time series of
all nodes so that the KLD regularization term can be calculated. We need to be
careful to generate multiple training instances and mark the node
corresponding to each instance, which facilitates the construction of batches
and the calculation of loss. The forecast given by the above model after
training is not entirely coherent, because the penalty term only helps the
model give an approximately coherent forecast. One more step of bottom-up
aggregation is required to give coherent forecasts.

The following subsection presents one deep learning model, DeepAR, in this
study for reconciliation. The basic idea of DeepAR, training mode, and details
of implementation will be given.

###  III-B DeepAR

DeepAR is a time series probabilistic forecast method based on deep learning
proposed by Flunkert _et al._ [5]. It has been applied to various forecasting
problems with excellent results. DeepAR is a nonlinear form of the classical
autoregressive model, which can be applied to single time series or extended
to multiple time series. DeepAR is based on Recurrent Neural Network (RNN).
The objective of DeepAR is to obtain the joint conditional probabilistic
distribution
$P\left(\bm{y}_{i,t_{0}:T}\mid\bm{y}_{i,1:t_{0}-1},\bm{x}_{i,1:T}\right)$,
i.e. according to the existing data $\bm{y}_{i,1:t_{0}-1}$ and covariates
$\bm{x}_{i,1:T}$, modeling the future sequence $\bm{y}_{i,t_{0}:T}$. The above
distribution can be written in the following likelihood form:

$\begin{gathered}Q_{\Theta}\left(\bm{y}_{i,t_{0}:T}\mid\bm{y}_{i,1:t_{0}-1},\bm{x}_{i,1:T}\right)=\prod_{t=t_{0}}^{T}Q_{\Theta}\left(y_{i,t}\mid\bm{y}_{i,1:t_{0}-1},\bm{x}_{i,1:T}\right)\\\
=\prod_{t=t_{0}}^{T}\ell\left(y_{i,t}\mid\theta\left(\mathbf{h}_{i,t},\Theta\right)\right),\end{gathered}$

where
$\mathbf{h}_{i,t}=h\left(\mathbf{h}_{i,t-1},y_{i,t-1},\bm{x}_{i,t},\Theta\right)$
is the output of RNN. We input the hidden state $\mathbf{h}_{i,t-1}$ at the
last moment and data $y_{i,t-1}$, as well as covariate $\bm{x}_{i,t}$, then
RNN will output $\mathbf{h}_{i,t}$, which is then transformed into the
parameters of a given distribution by a specific transformation
$\theta(\cdot)$, and then the likelihood function is obtained. In the
experiment of this paper, we use Gaussian distribution to calculate the
likelihood function. The mean value $\mu$ of Gaussian distribution is given by
the affine function of network output, and the standard deviation $\sigma$ is
obtained by applying affine transformation and softplus activation function to
ensure $\sigma>0$, i.e.

$\mu\left(\mathbf{h}_{i,t}\right)=\mathbf{w}_{\mu}^{T}\mathbf{h}_{i,t}+b_{\mu},\sigma\left(\mathbf{h}_{i,t}\right)=\log\left(1+\exp\left(\mathbf{w}_{\sigma}^{T}\mathbf{h}_{i,t}+b_{\sigma}\right)\right).$

Then the negative log-likelihood loss function is as follows:

$L=-\sum_{i=1}^{N}\sum_{t=t_{0}}^{T}\log\ell\left(y_{i,t}\mid\theta\left(\mathbf{h}_{i,t}\right)\right).$

We use multivariate DeepAR, so the likelihood function is multiplied by the
likelihood function of all nodes. When applying DeepAR to our proposed
reconciliation method, we add KLD regularization term to this loss function.

---
Figure 2:  DeepAR training and prediction process.

The training and prediction process of DeepAR can be represented by Figure 2.
There is a difference between the training process and the prediction process,
the label values are known at the time of training, so each network input is
the time lag term of the real value. For the prediction process, the input of
each round of the network is derived from the mean of the distribution of the
previous round of network output.

In subsequent experiments, we apply DeepAR to our proposed reconciliation
method for probabilistic forecasting and name it DeepAR-Hier.

##  IV Experiments

###  IV-A Data

In order to verify the effectiveness of the proposed approach, we conduct
experiments on three classical data sets in the field of hierarchical
forecast, i.e. the Australian infant mortality data set, the Australian
tourism data set, and the Wikipedia page view data set. The background and
features of these data sets are detailed below.

The Australian infant mortality data set, abbreviated as Infant, is from the
Australian Social Science Data Archive and includes infant mortality rates
across eight Australian states: New South Wales (NSW), Victoria (VIC),
Queensland (QLD), South Australia (SA), Western Australia (WA), Northern
Territory (NT), Australian Capital Territory (ACTOT), and Tasmania (TAS). This
data set is grouped by gender and state. The time range of available data is
from 1901 to 2003. The data set is publicly available in R’s addb package
[10]. Due to missing values, we use annual data from 1933 to 2003 in our
analysis, and the forecast horizon is 4.

The Australian tourism data set, abbreviated as Tourism, contains overnight
trips in Australia from the first quarter of 1998 to the fourth quarter of
2016 [11]. The data are quarterly and grouped by three variables: tourism
region, eight Australian states, and purpose of travel. For this data set, the
forecast horizon is set to 4.

Wikipedia page views data set, abbreviated as Wiki, gives the daily page views
of 145,000 different Wikipedia articles from 2015-07-01 to 2016-12-31, from
which we sample 32 bottom time series and aggregate them to get the upper time
series. Specifically, the data are grouped by three variables: agent type,
access type, and country code. In the experiment, the forecast horizon is set
to 15 days.

###  IV-B Experimental set-up

We choose DeepAR as the deep learning model in our proposed approach. We
conducted a grid search for some important hyperparameters of DeepAR-Hier,
including the hidden recurrent size between (10,30), the number of RNN layers
between (2,5), dropout rate between (0.1,0.2), the number of epochs (20,60)
and the coefficient of the KLD regularization term $\lambda$ between(0,1).
According to the time length of different data sets, we use cross validation
with a specific number of folds to search parameters.

TABLE I: All methods available for comparison. DeepAR-Hier | Arima-Random-Bu-None
---|---
Pure-DeepAR | Ets-Stack-MinT-Struct
Arima-Stack-MinT-Struct | Ets-Stack-Bu-None
Arima-Stack-Bu-None | Ets-Rank-MinT-Struct
Arima-Rank-MinT-Struct | Ets-Rank-Bu-None
Arima-Rank-Bu-None | Ets-Random-MinT-Struct
Arima-Random-MinT-Struct | Ets-Random-Bu-None

We compare the proposed approach with the existing probabilistic forecast
reconciliation method mentioned in section II. In a word, there are 14 models
participating in the comparison summarised in Table I. For example, Arima-
Stack-MinT-Struct means that the base forecast is given by Arima, the joint
distribution sample is constructed by stack, and the reconciliation method is
MinT in which error covariance is constructed by STRUCT. The BU method does
not need to assume the error covariance, which is represented by None. Pure-
DeepAR refers to the original DeepAR method without changing the loss
function. We use PyTorch to implement DeepAR-Hier and Pure-DeepAR, use
Python’s Sktime library to implement Arima and ETS, and automatically tune the
model parameters of Arima and ETS.

We evaluate the forecasting performance of these methods in terms of
Continuous Ranked Probability Score (CRPS). CRPS is a commonly used evaluation
metric in the field of probabilistic forecasting. The definition of CRPS is as
follows:

$\operatorname{CRPS}\left(F^{f},F^{0}\right)=\int_{-\infty}^{+\infty}\left[F^{f}(x)-F^{0}(x)\right]^{2}\mathrm{~{}d}x,$

where $F^{f}$ is predictive CDF and $F^{0}$ is real CDF. CRPS measures the
differences between predictive distribution and real distribution. The lower
the CRPS, the more accurate the probabilistic forecast results. In most cases,
the real distribution is unknown, and we only know the observation value. If
we have real observations of prediction steps
$\xi_{1},\xi_{2},\ldots\ldots,\xi_{n}$, and the corresponding predictive
distributions $F_{1},F_{2},\ldots\ldots,F_{n}$, CRPS can be estimated by the
following equation:

$\displaystyle\operatorname{CRPS}(F,\xi)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\operatorname{crps}\left(F_{i},\xi_{i}\right)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\int_{-\infty}^{+\infty}\left[F_{i}(x)-\varepsilon\left(x-\xi_{i}\right)\right]^{2}\mathrm{~{}d}x,$

where $\varepsilon(t)$ is unit step function. There is no fixed threshold for
how small the CRPS is to be a more accurate probabilistic forecast and it
depends on different practical applications. The CRPS values under the same
situation or data set are comparable.

###  IV-C Result

TABLE II: Average CRPS. Method Dataset | Tourism | Infant | Wiki
---|---|---|---
DeepAR-Hier | 23.854 | 8.468 | 3187.048
Pure-DeepAR | 28.909 | 12.155 | 3841.136
Arima-Stack-MinT-Struct | 25.036 | 11.532 | 18758.610
Arima-Stack-Bu-None | 27.140 | 12.136 | 28156.813
Arima-Rank-MinT-Struct | 23.364 | 15.140 | 17829.531
Arima-Rank-Bu-None | 24.206 | 16.349 | 27377.863
Arima-Random-MinT-Struct | 25.329 | 11.527 | 18587.533
Arima-Random-Bu-None | 27.138 | 12.161 | 27913.228
Ets-Stack-MinT-Struct | 21.312 | 11.609 | 16312.947
Ets-Stack-Bu-None | 26.589 | 12.996 | 23884.263
Ets-Rank-MinT-Struct | 20.724 | 15.563 | 15663.066
Ets-Rank-Bu-None | 23.998 | 16.545 | 23162.240
Ets-Random-MinT-Struct | 21.072 | 11.602 | 16386.476
Ets-Random-Bu-None | 26.673 | 12.960 | 23055.768

Using the average CRPS of all series in one data set as the evaluation metric,
the performance of each method is shown in Table II, where the minimum average
CRPS value of each data set has been bolded. For Infant and Wiki data sets,
DeepAR-Hier is the best approach compared with other probabilistic forecast
reconciliation methods. For the Tourism data set, the performance of DeepAR-
Hier is worse than Ets-Rank-MinT-Struct.

TABLE III: Average CRPS at each level. Dataset | Level | $\lambda$ | DeepAR-Hier
---|---|---|---
Tourism | 1 | 0.006 | 1014.807
2 | 241.295
3 | 35.892
4 | 12.685
Infant | 1 | 0.12 | 32.280
2 | 9.694
3 | 6.366
Wiki | 1 | 0.24 | 31991.549
2 | 6273.763
3 | 2469.683
4 | 2001.627

We calculated the average CRPS at each level for the three data sets, i.e. the
mean of CRPS values of all time series at a certain level. The results are
shown in Table III, where the $\lambda$ column gives the optimal $\lambda$
value determined on the validation set of each data set. The CRPS value in
bold indicates that DeepAR-Hier is the best method at this level. It can be
seen that for the Tourism data set, although DeepAR-Hier is not as good as
other methods in terms of the average CRPS of all nodes, at the bottom level,
the average CRPS of the DeepAR-Hier is the smallest. For the Wiki and Infant
data set, DeepAR-Hier outperforms the other methods at all levels. In general,
DeepAR-Hier is better than other probabilistic forecast reconciliation methods
in most situations, but specifically, DeepAR-Hier can improve the
probabilistic forecast accuracy of the bottom level, and performance at upper
levels is not as good as that at lower levels.

---
Figure 3:  MCB result of Tourism dataset.
---
Figure 4:  MCB result of Infant dataset.
---
Figure 5:  MCB result of Wiki dataset.

In order to evaluate the performance of each method more appropriately, we use
multiple comparisons with the best (MCB) proposed by Koning _et al._ to test
the statistical significance of the performance differences of multiple
methods. MCB focuses on the average rank, i.e. the average rank of the
evaluation metrics of all nodes [12]. MCB calculates the average rank and a
confidence interval for each method to be compared. In Figure 3, 4 and 5, we
present the average rank of each method and its confidence interval, and the
results for each data set are shown in different panels. The best method for
each data set is shown at the bottom of the panel, and the worst-performing
method is shown at the top of each panel. If the intervals do not overlap, it
indicates that the difference between the two methods in forecasting
performance is statistically significant. It can be seen that for the Infant
and Wiki datasets, the DeepAR-Hier method performs best. Although the average
CRPS of DeepAR-Hier for the Tourism dataset is not the smallest, the forecast
performance of DeepAR-Hier does not differ significantly from that of other
methods because of interval overlapping.

##  V Conclusions and discussion

In this paper, we propose a probabilistic forecast reconciliation method based
on deep learning with Kullback-Leibler divergence regularization. We use a
deep learning model to predict multivariate time series and introduce the
Kullback-Leibler divergence regularization term into the loss function. This
approach can fuse the prediction step and reconciliation step into a deep
learning framework, instead of treating the reconciliation step as a fixed and
hard post-processing step, making the reconciliation step more flexible and
soft.

In our approach, the deep learning method is not specified, all global
probabilistic forecasting models for multivariate time series based on deep
learning can be used in our approach. We apply DeepAR to our proposed
reconciliation method for probabilistic forecasting and conduct experiments on
three hierarchical time series data sets. The experiment result shows that our
approach is better than other probabilistic forecast reconciliation methods in
most situations, but specifically, our approach can improve the probabilistic
forecast accuracy of the bottom level, and performance at upper levels is not
as good as that at lower levels.

There are several valuable directions worthy of further investigation. First,
the approach we propose is aimed at cross-sectional data, in the future, we
can try to apply the idea to temporal hierarchy and cross-temporal hierarchy.
Second, DeepAR-Hier performs best on the Wikipedia page views data set, which
is related to the fact that deep learning methods are suitable for large-
volume data. In the future, deep learning frameworks suitable for small sample
size data should be studied, and the minimum length of training time series
should also be explored. It can be considered to conduct an experiment to test
the trend of accuracy change relative to the training scale to determine the
optimal training length. Finally, in our experiments, we assume that the
predictive distribution is Gaussian distribution and we can try to replace the
predictive distribution in the future.

## Acknowledgments

We thank the reviewers for helpful comments that improve the contents of this
paper. This research was supported by the National Social Science Fund of
China (22BTJ028).

## References

  * [1] R. J. Hyndman, A. J. Lee, and E. Wang, “Fast computation of reconciled forecasts for hierarchical and grouped time series,” _Computational Statistics & Data Analysis_, vol. 97, pp. 16–32, 2016.
  * [2] Y. Zhang, J. Wang, and X. Wang, “Review on probabilistic forecasting of wind power generation,” _Renewable and Sustainable Energy Reviews_ , vol. 32, pp. 255–270, 2014.
  * [3] J. Jeon, A. Panagiotelis, and F. Petropoulos, “Probabilistic forecast reconciliation with applications to wind power and electric load,” _European Journal of Operational Research_ , vol. 279, no. 2, pp. 364–379, 2019.
  * [4] S. S. Rangapuram, L. Werner, K. Benidis, P. Mercado, J. Gasthaus, and T. Januschowski, “End-to-end learning of coherent probabilistic forecasts for hierarchical time series,” in _International Conference on Machine Learning_ , 2021.
  * [5] D. Salinas, V. Flunkert, J. Gasthaus, and T. Januschowski, “Deepar: Probabilistic forecasting with autoregressive recurrent networks,” _International Journal of Forecasting_ , vol. 36, no. 3, pp. 1181–1191, 2020\.
  * [6] B. N. Oreshkin, D. Carpov, N. Chapados, and Y. Bengio, “N-BEATS: neural basis expansion analysis for interpretable time series forecasting,” _CoRR_ , vol. abs/1905.10437, 2019. [Online]. Available: http://arxiv.org/abs/1905.10437
  * [7] B. Lim, S. . Arık, N. Loeff, and T. Pfister, “Temporal fusion transformers for interpretable multi-horizon time series forecasting,” _International Journal of Forecasting_ , vol. 37, no. 4, pp. 1748–1764, 2021\.
  * [8] S. B. Taieb, J. W. Taylor, and R. J. Hyndman, “Hierarchical probabilistic forecasting of electricity demand with smart meter data,” _Journal of the American Statistical Association_ , vol. 116, no. 533, pp. 27–43, 2021.
  * [9] S. L. Wickramasuriya, G. Athanasopoulos, and R. J. Hyndman, “Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization,” _Journal of the American Statistical Association_ , vol. 114, no. 526, pp. 804–819, 2019.
  * [10] R. Hyndman, _addb: Australian Demographic Data Bank_ , 2010.
  * [11] R. Hyndman and G. Athanasopoulos, _Forecasting: Principles and Practice_ , 3rd ed. Australia: OTexts, 2021.
  * [12] A. J. Koning, P. H. Franses, M. Hibon, and H. Stekler, “The m3 competition: Statistical tests of the results,” _International Journal of Forecasting_ , vol. 21, no. 3, pp. 397–409, 2005.


11institutetext: KTH Royal Institute of Technology, Stockholm, Sweden
11email: {msnkhan ${}^{\textrm{({\char 0\relax})}}$, nreje<EMAIL_ADDRESS>
# Utility Assessment of Synthetic Data Generation Methods

Md Sakib Nizam Khan ${}^{\textrm{({\char 0\relax})}}$ 11 Niklas Reje 11 Sonja
Buchegger 11

###### Abstract

Big data analysis poses the dual problem of privacy preservation and utility,
i.e., how accurate data analyses remain after transforming original data in
order to protect the privacy of the individuals that the data is about - and
whether they are accurate enough to be meaningful. In this paper, we thus
investigate across several datasets whether different methods of generating
fully synthetic data vary in their utility a priori (when the specific
analyses to be performed on the data are not known yet), how closely their
results conform to analyses on original data a posteriori, and whether these
two effects are correlated. We find some methods (decision-tree based) to
perform better than others across the board, sizeable effects of some choices
of imputation parameters (notably number of released datasets), no correlation
between broad utility metrics and analysis accuracy, and varying correlations
for narrow metrics. We did get promising findings for classification tasks
when using synthetic data for training machine learning models, which we
consider worth exploring further also in terms of mitigating privacy attacks
against ML models such as membership inference and model inversion.

###### Keywords:

Synthetic Data Utility Metrics Analysis Correlation.

##  1 Introduction

In this era of big data and artificial intelligence, technologies are becoming
increasingly dependent on data processing and analysis. Since a fair share of
the data used by these technologies is the privacy-sensitive data of
individuals, there is a growing need for methods that facilitate privacy-
preserving data analysis and sharing. Differential Privacy (DP) [13] and
k-anonymity-related [36] mechanisms focus on providing provable privacy
guarantees, at the cost of utility [1] due to the noise additions required by
the mechanisms. The mechanisms either transform data/models/results or
restrict possible queries from a database to achieve the guarantees. For
instance, DP was originally designed for interactive statistical queries to a
database and later on extended to non-interactive settings including microdata
release. Even though for interactive settings it can provide mathematically
quantifiable guarantees, it has been shown that for microdata releases it is
not able to provide similar specific confidentiality guarantees [22].

Statistical disclosure control (SDC) is another research domain that has been
very active in addressing the issue of protecting the confidentiality of
individuals while maximizing the utility provided by the released data. An
alternative to traditional SDC approaches is publication of synthetic data, an
idea that is generally considered as part of SDC literature [34]. Synthetic
data is generated by drawing from a model that is fitted on original data.
Firstly, with synthetic data, there is no given direct linkability from
records to individuals as the records are not real 111In [34] Ruiz et al.
however argue that, based on a scenario where the attacker has access to the
original dataset in its entirety, individuals’ representations in synthetic
and original datasets remain linked by the information they convey.. Secondly,
given the large variety of ways to generate synthetic data, there is potential
for better utility to arrive at an acceptable level of accuracy for a given
analysis.

The main idea of statistical methods for synthetic data generation is
imputation [33] which is originally developed for replacing missing data in
survey results. There are different choices of parameters for the synthesis
process such as which other variables to choose to synthesize one variable,
the order of variables chosen for synthesis, etc. which can impact the utility
of the generated data. Similarly, depending on the dataset the choice of the
synthesizer can have a large impact on the utility of the synthetic data since
some synthesizers perform better for numerical variables and some perform
better with categorical variables. There has been a lot of research on
synthetic data utility in recent years however, more research is needed
concerning the impact of different choices during the synthetic data
generation process on the utility. Further investigation is also required to
find out if there is any correlation between the general utility metrics
commonly used to evaluate synthetic data and how well the synthetic data
performs for a given analysis.

In this work, we thus investigate the utility provided by different synthetic
data generation techniques and synthesis parameters, on separate but similar
datasets (Adult and Polish, both census-type data) and on a dataset with
different characteristics (Avila). First, we determine the effects on utility
as measured by various metrics for the similarity between the original and the
synthetic data and, second, by comparing the similarity of results from
analyses performed on the original versus on the synthetic data. Third, we
investigate to what extent the first can predict the second when the analysis
is not known beforehand. In summary, our main contributions are:

  * •

A comprehensive evaluation of synthetic data utility using three publicly
available datasets.

  * •

Identification of the individual impact of the choice of variables during
synthesis, synthesis order, proper and non-proper synthesis, and number of
datasets on the utility of synthetic data.

  * •

A comprehensive study on the correlation between different utility metrics and
how well the synthetic data performs for any given analysis.

Organization. The rest of the paper is organized as follows. In Section 2, we
discuss the related work, followed by an overview on the synthetic data
generation process and the utility metrics in Section 3. We then present our
experimental results in Section 4, followed by our discussion and conclusion
in Section 5.

##  2 Related Work

There have been various research works concerning synthetic data generation
and the utility provided by them. The related works that we have distilled
from the literature mostly focus on the comparison between imputation methods
(single and multiple imputation [37], hierarchical Bayes and conventional
generalized linear imputation models [15]), evaluation of different mechanisms
and tools for synthetic data generation (fully and partially synthetic
data[3], tools for regression [17] and classification [16] tasks), evaluation
of the usefulness of synthetic data [23, 14, 38, 10, 21], and what utility
metrics to use for comparing original and synthetic data [5, 35]. In this
section, we review some of these works in detail.

Among the works on imputation methods, Taub et al. [37] evaluated the real
word analyses replicability of singly and multiply imputed CART generated
synthetic data using Purdam and Elliot’s [27] methodology. To evaluate the
impact of disclosure control on the analysis outcomes, Purdam and Elliot [27]
replicated the published analysis on datasets using the disclosure controlled
versions of the same datasets. Based on this methodology, Taub et al. [37]
replicated 9 different sets of analyses involving 28 different tests and
models. According to the findings of the authors, multiply imputation
performed better for some analyses, but not for all, and depending on the
complicacy of the analysis, single imputation can be useful in some scenarios.
The authors also investigated the relationship between the utility metrics and
how well a synthetic dataset performs for a given analysis and found that
there is no clear relationship.

To investigate the usefulness of synthetic data, Nowok [23] evaluated the
performance of non-parametric tree-based synthetic data generation methods
(Classification and Regression Trees (CART), bagging, and random forests)
using synthpop. Besides general utility measures, the authors used some
hypothetical analyses to evaluate analysis specific utility of the synthetic
data. From the empirical evaluation, the authors conclude that it is possible
to produce useful completely synthetic data using automated methods.
Similarly, Drechsler et al. [10] replicated an already published analysis of a
dataset using synthetic data and found that the regression coefficients of the
synthetic and the original were almost identical and concluded that the
authors of the published analysis would have drawn the same conclusion using
the synthetic dataset. Lee et al. [21] used univariate, bivariate, and linear
regression-based exploratory data analysis to evaluate the utility of
synthetic data generated using CART models and found that for univariate
analysis synthetic data results match the original data, whereas for bivariate
and linear regression it was not true. According to the authors, the reason
for this could be that the CART model underestimates the strong correlation
between the variables in the dataset.

Among the works on utility metrics, Dankar et al. [5], first classified the
available utility metrics for synthetic data comparison into different
categories based on the measure they attempt to preserve. Then the authors
chose one metric from each category depending on popularity and consistency
and used them to compare the utility of four data synthesizers (i.e Data-
Synthesizer (DS), Synthetic Data Vault (SDV), Synthpop Parametric, Synthpop
Non-parametric). According to the authors, their experimental results show
that the Synthpop Non-parametric is the best performing synthesizer overall
and provides the best average values across all metrics as well as the best
(overall) stability and consistency. Snoke et al. [35] investigated how
general utility compares to the specific utility for synthetic data and also
presented two contrasting example evaluations. The authors conclude that
general utility measures can be used to tailor the methods that are used to
synthesize datasets for specific purposes whereas specific utility measures
can be used for reassurance after performing any standard exploratory analysis
on the synthetic data that the data was not misleading.

Even though there have been a lot of studies on different aspects of synthetic
data utility, more research is needed to better understand how choice and
order of variables during synthesis, proper and non-proper synthesis, number
of datasets, etc. affect the utility of generated synthetic data. Similarly,
there is a lack of concrete work on the correlation between commonly used
utility metrics for synthetic data and the performance of synthetic data in
any given analysis.

##  3 Synthetic Data Generation

Figure 1 provides an overview of synthetic generation process and its two main
components, Data Synthesis and Evaluation.

Figure 1: Synthetic Data Generation Process

###  3.1 Data Synthesis

####  3.1.1 Choice of Input Parameters.

The goal of the data synthesis component is to generate synthetic datasets
based on the original dataset. The synthesis process requires a set of input
parameters alongside the original data.

In statistics, imputation is commonly used for solving the problem of missing
data in surveys and the idea is to replace them with substitute values drawn
from similar records. Rubin [32] developed an approach called multiple
imputation and later on [33] proposed to use it for generating fully synthetic
datasets. The basic idea proposed by him is to treat the population not
selected in a sample as missing data and then use multiple imputation to
create synthetic datasets. For multiple imputations, two general approaches
have emerged over the years: the joint modeling approach and the fully
conditional specification (FCS) approach [7]. With the FCS approach variables
are replaced in a specific order whereas with the joint modeling approach
multiple variables are replaced together. More details about these approaches
can be found in [7]. Here, we present the three main parameter choices for
synthesis.

Synthesis Order. For the FCS approach, a particular order of variables needs
to be selected for the synthesis. The order that seems most logical is usually
chosen. There are also other approaches [30] for choosing the order.
Nonetheless, the order needs to be chosen based on the relations between the
variables and also the computational complexity. Since the order impacts which
relationships are kept between the variables, it can also impact the utility
of synthetic data

Simple & Selective Synthesis. For the FCS approach, it is possible to select
only a subset of variables for the synthesis of another variable (we term this
approach Selective) instead of using all other variables (which we call
Simple). The Selective approach can be taken in an effort to reduce
computation time or to preserve the relationships between variables that are
most essential. However, the downside is that the correlations that are not
chosen can get lost, and knowing which relationships to preserve is also not
straightforward.

Proper & Non-proper Synthesis. Depending on whether the new values of the
imputation parameters are drawn from their posterior predictive distribution
(Proper) or not (Non-proper), synthesis can be of two types.

####  3.1.2 Choice of Synthesizers.

There are several methods developed over the years for generating synthetic
data using multiple imputation. We use three of the most widely used methods:
Parametric and Decision Trees based on the FCS approach and Saturated model
based on the joint modeling approach. The parametric method uses different
kinds of regression models depending on what type of data values are meant to
be synthesized with linear regression used for numerical variables and
logistic regression used for categorical variables. Decision Trees are a non-
parametric alternative to the parametric methodology [25]. The Saturated Model
approach works by fitting a model that perfectly reproduces the data, to a
cross-tabulated table of all the categorical variables [20, 26]. Since this
method only works with categorical variables, a parametric method or decision
tree is needed for the numerical variables. More details on each of these
methods can be found in [31].

###  3.2 Evaluation

For evaluating the utility of synthetic datasets a common approach found in
the literature is to use different utility metrics. Besides this, some studies
have also used the accuracy of analysis to measure how well synthetic data
imitate the original data for any given analysis. In this section, we briefly
discuss these two comparison approaches.

####  3.2.1 General Utility Metrics.

The existing works on synthetic data use different metrics to estimate the
utility provided by such data. However, there is no such metric that can
individually capture all the different aspects of utility concerning synthetic
data. Hence, a combination of different metrics is needed to compare the
utility provided by synthetic data and the original data. There are two
categories [18, 11] of utility measures which are commonly termed broad
measures and narrow measures of utility. The broad measures mostly compare the
utility between the entire distributions of synthetic and original data using
some statistical distance measures such as Kullback-Leibler (KL) divergence
[18] score. Narrow measures on the other hand tend to compare specific models
(e.g., regression, point estimation, etc.) between the synthetic and original
dataset and are more widely used in the literature for the utility comparison
of synthetic data. Confidence Interval Overlap (CIO) developed by Karr et al.
[18] is a widely used [9, 11, 37] narrow measure of utility for synthetic
data. In this work, we use the both CIO and KL divergence as the general
utility metrics.

Confidence Interval Overlap (CIO). To measure the utility of synthetic data
based on confidence interval overlap (CIO), we use a combination of mean point
and regression fit coefficient estimation. We calculate the 95% Confidence
Interval (CI) for the mean points and the regression fit coefficients from
both the synthetic and the original data and then use the formula proposed by
Karr et al in [18] to measure the CIO. We calculate the average CIO for all
the coefficients of each regression fit and point estimation. We then take the
average over all the tests for each $m$ and $k$ where $m$ is number of
datasets to release and $k$ is the number tests performed (i.e., the number of
different sets of $m$ datasets). For the point estimates and regression fits,
we need to estimate the variance. There are different equations for variance
estimation based on combining rules for imputation proposed by Raghunathan et
al. [29]. For more details about the imputation combining rules and variance
estimation see [8, 37]. In this work, we use the equation proposed by Raab et
al. [28] for the variance estimation. The reason for choosing the equation is
that it is not dependent on the number of datasets and thus valid even if
there is only one synthetic dataset.

Kullback-Leibler Divergence. Kullback-Leibler Divergence is used to measure
how close two probability distributions are, based on information entropy and
is a global measure of data utility [39]. We measure the KL divergence score
of each variable individually and then take the average over all the variables
for the comparison between synthesizers.

Average Percentage Overlap (APO) above 90%. A CIO of 90% or higher has been
regarded as good. With the APO measure, we can have a measure of confidence in
the utility beyond the average CIO by calculating the ratio of the CIOs over
90% to all CIOs.

####  3.2.2 Analysis Accuracy.

The problem with broad and narrow measures is that none of them is
individually sufficient for measuring the utility of synthetic data. They both
have certain weaknesses. For instance, narrow measures perform really well for
certain analyses but are unable to provide any insights on the overall utility
of the dataset whereas broad measures perform well for many analyses but
really well for none [18]. To overcome these issues, multiple studies use
either already published analysis [27, 37] on original data or ad hoc analysis
[23, 20] to measure the utility provided by synthetic data. Since this
approach uses multiple analyses for the utility comparison, it avoids the
possibility of reporting unreasonably high or low utility measures based on a
single measure which is the case for broad and narrow measures. In this work,
for measuring the utility, we use machine learning-based classification
accuracy along with the accuracy of multiple ad hoc analyses which compare the
multivariate relationship and univariate distributions between synthetic and
original data.

##  4 Experimental Results

###  4.1 Datasets and Experimental Environment

In this work, for synthetic data generation, we use the R package synthpop
developed by Nowok et al. [24]. For the datasets, we use three publicly
available ones, i.e., Polish quality of life dataset (SD2011) [2] from
synthpop example datasets, Adult dataset [19], and Avila dataset [6] from UCI
Machine Learning Repository [12]. While both the Polish and the Adult datasets
are standard census surveys, Avila is a completely different type that is not
commonly used in existing works on synthetic data generation. We choose Avila
because we want to see how well the methods which are mostly developed for
census or similar datasets perform for a completely different dataset. More
details about the datasets and the modifications that we performed on them can
be found in Appendix.

Table 1: Naming Convention

Synthesizer | Symbol |  | Name | Suffix |  |  |
---|---|---|---|---|---|---|---
Standard Parametric (SAP) | P |  | Original Ordering | No Suffix |  | Name | Suffix
| Decision Tree / Classification &
---
Regression Tree (CART)
D |  | Opposite Ordering | O |  |  |
| Saturated Model (Catall) with
---
Predictive Mean Matching (CAP)
CP |  | Own Ordering | V |  | Non-Proper | No Suffix
| Saturated Model (Catall) with
---
CART (CAC)
CC |  |  | Largest Categorical
---
Variables First
H |  |  |
Sample | S |  Imputation Order |  | Largest Categorical
---
Variables Last
L |  Synthesis Method | Proper | T

For the experiments, we generated multiple synthetic datasets for Standard
Parametric (P), Cart (D), Catall (i.e., saturated model) with PMM (CP), and
Catall with Cart (CC) 222We do not preprocess the data since there is no clear
benefit of doing so [4].. For P and D, we generated synthetic datasets with 3
different orders, a total of 6 different datasets. Next, for CP, CC, and the 6
datasets for P and D we generated Proper and Non-Proper combinations which sum
up to 16 datasets. Additionally, for baseline comparison, we created a fifth
synthesizer, Sample, which just samples values from the original dataset for
each variable without attempting to retain any relationships between the
variables. Thus, for each of the 3 original datasets, we generate 17 dataset
collections. In the case of synthesis order, just for the sake of comparison,
we generated our own ordering where we chose an order that we thought made the
most sense concerning preserving relationships between variables. Table 1
shows how we symbolize each method and parameter in this work. We use the
naming convention shown in the table for the rest of this paper. Our codes for
all the experiments are available on
GitHub333https://github.com/sakib570/synthetic-data-utility.

###  4.2 General Utility Metrics

####  4.2.1 Confidence Interval Overlap (CIO).

To evaluate the utility of synthetic datasets using confidence interval
overlap (CIO), we perform mean point estimation and regression fitting. The
mean point estimations are performed on the numerical variables of each
dataset. In the case of regression fitting, we perform 20 regression fits on
the Polish dataset and 24 regression fits each on the Adult and Avila
datasets. The details of regression fits for each dataset can be found in
[31].

##### Mean Point Estimation.

The mean point estimation results for all three datasets show that at least
$m=3$ or higher is required for all variables to reach APO above $90\%$. In
terms of synthesizer, though CP performed well for Adult and Polish, it
struggled for Avila whereas D and CC performed well for all three datasets. P
struggled for both Adult and Avila and was unable to reach APO above $90\%$
regardless of the number of datasets released. The results for mean point
estimation can be found in Table 3 under Appendix.

##### Regression Fit.

For regression fit, we calculate both the average CIO and the APO above $90$%
over all the regression fits for each dataset. Figure 2 shows the average CIO
(left column) and APO above 90% (right column) for the regression fits for all
three datasets.

The average CIO shows that CC is the best synthesizer by being the only one to
reach above 90% for Polish and achieving the highest scores for both Adult and
Avila. However, for Adult and Avila, none of the synthesizers can reach above
90% and in the case of Avila, the best one is only able to reach 60%. CP and D
are the next best ones for all three datasets. The Sample method performs far
worse than the others which confirms that the other synthesizers have been
able to retain relationships between the variables. We also see that there is
a steady increase in overlap as $m$ increases, though there appears to be a
somewhat diminishing return with the largest increases being below 10 to 20
datasets and the overlaps above only raising the results a few percentage
points at most. A key difference between the results of Polish and Adult/Avila
is that for the latter the average CIO drops for many of the synthesizers as
$m$ increases, something which only happens for Sample in the case of Polish.
Non-proper methods perform better overall than proper ones. In terms of
synthesis order, for Polish the original order (O), for Adult in the case of D
the L ordering and for P the H ordering, and for Avila own ordering (V)
perform better than other orderings.

(a) CIO -Polish Dataset

(b) APO above 90% - Polish Dataset

(c) CIO - Adult Dataset

(d) APO above 90% - Adult Dataset

(e) CIO - Avila Dataset

(f) APO above 90% - Avila Dataset

Figure 2: CIO and APO above 90% for Regression Fits

The APO above 90% also shows similar results as average CIO in terms of
synthesizer performance where CC, CP, and D perform well for all three
datasets. The major difference between average CIO and APO above 90% is the
results for lower $m$ where APO 90% shows that lower $m$ (i.e., $m<\sim 3$)
perform very poorly utility wise compare to higher $m$ which is not the case
for average CIO. We can also see that the APO above $90\%$ had a steady
increase as the number of datasets to release increases. The increase of APO
above $90$% as $m$ increases is also a lot more than for average CIO and is
more significant for higher $m$ as well. For all three datasets for lower $m$
(e.g. $m=\\{1,2,3\\}$) proper methods perform much worse than Non-Proper
methods regardless of synthesizer combinations. In terms of synthesis order,
we see similar results as average CIO.

####  4.2.2 KL-Divergence.

For the KL-Divergence score (Appendix Table 4), we calculate the score for all
variables and normalize them over the score of Sample and then take the
average over all variables for each synthesizer for easier comparison. The KL
divergence results conform to the CIO results that CC and D perform well for
all three datasets and also CP is not far off. In terms of synthesis order,
the original ordering performs better with few exceptions for all three
datasets. Similarly, non-proper performs much better than proper in most
cases.

###  4.3 Analysis Accuracy

Figure 3: Classification Accuracy

####  4.3.1 Classification Accuracy.

Figure 3 shows the classification accuracy results for all three datasets.
From the results we find that CC, CP and D perform well and P is consistently
poor for all three datasets. We tested the synthetic model accuracy using the
original dataset and found that the accuracy score is similar to the original
model. We also look at whether the original model and the synthetic model
correctly and incorrectly classify same records or they differ. The
investigation reveals that for all three datasets the synthetic model and the
original model classify the same records correctly and incorrectly and only
differ depending on the difference in accuracy score. The details of the
classification tasks and model generation process can be found under Appendix.

(a) Polish: Male Education Vocational/Grammar

(b) Adult: US Male Income $\leq 50K$

(c) Avila: Exploit Variable Distribution

Figure 4: Ad-hoc Analysis on the Datasets

####  4.3.2 Ad hoc Analysis.

To determine the effectiveness of synthetic data in different analyses, we
look at the multivariate relationship and univariate distributions of the
variables present in the datasets. Figure 4 shows the results of ad hoc
analysis on all three datasets. For Polish, Adult, and Avila, we perform
bivariate, multivariate (i.e., 3 variables), and univariate analysis
respectively. In the case of Polish and Adult synthesizers CC, CP, and D
perform well. In terms of synthesis order, some orders have significant
performance improvements over the others for both Polish and Adult. Finally,
for Avila, except CPT all synthesizers have similar distributions as the
original dataset. In terms of synthesis method, in all three scenarios, we see
that non-proper synthesis performs better overall than proper.

Table 2: Correlation between Utility Metrics and Analysis Accuracy

Polish Dataset | Adult Dataset | Avila Dataset
---|---|---
Relation | P Value |  | Correlation
---
Co-efficient
Relation | P Value |  | Correlation
---
Co-efficient
Relation | P Value |  | Correlation
---
Co-efficient
Correlation between Classification Accuracy and Utility Metrics
| Accuracy of Income Classification
---
& CIO
1.35E-05 | -0.85 |  | Accuracy of Income Classification
---
& CIO
7.48E-05 | -0.81 |  | Accuracy of Author Classification
---
& CIO
5.47E-07 | -0.91
| Accuracy of Income Classification
---
& KL Divergence
0.5338 | -0.16 |  | Accuracy of Income Classification
---
& KL Divergence
0.6213 | -0.13 |  | Accuracy of Author Classification
---
& KL Divergence
2.43E-02 | 0.54
| Accuracy of Income Classification
---
& APO 90%
1.54E-03 | -0.71 |  | Accuracy of Income Classification
---
& APO 90%
2.15E-02 | -0.55 |  | Accuracy of Author Classification
---
& APO 90%
3.54E-06 | -0.88
Correlation between Different Utility Metrics
CIO & APO above 90% | 4.14E-07 | 0.91 | CIO & APO 90% | 5.54E-06 | 0.87 | CIO & APO 90% | 3.04E-11 | 0.98
CIO & KL Divergence | 3.97E-01 | 0.22 | CIO & KL Divergence | 1.69E-01 | -0.35 | CIO & KL Divergence | 9.80E-02 | -0.41
APO 90% & KL Divergence | 8.43E-01 | -0.05 | APO 90% & KL Divergence | 2.28E-01 | -0.31 | APO 90% & KL Divergence | 4.66E-02 | -0.49
Correlation between Ad hoc Analysis and Utility Metrics
APO 90% & Female Inc 1K to 2K | 0.06245 | -0.46 | APO 90% & Male <=50K US | 8.42E-03 | -0.62 | APO 90% & Units Class A | 0.2718 | -0.28
APO 90% & Male Edu Voc/Grm | 4.80E-04 | -0.75 | APO 90% & Units Income <=50K | 0.4505 | -0.20 | CIO & Units Class A | 0.3071 | -0.26
APO 90% & Units No Income | 0.5931 | -0.14 | APO 90% & White below 50K | 0.0776 | -0.44 |  |  |
CIO & Female Inc 1K to 2K | 0.04712 | -0.49 | CIO & Male <=50K US | 1.12E-04 | -0.80 |  |  |
CIO & Male Edu Voc/Grm | 1.14E-06 | -0.90 | CIO & Units Income <=50K | 0.4568 | -0.19 |  |  |
CIO & Married alcabuse Yes | 0.8505 | -0.05 | CIO & White above 50K | 0.08458 | -0.43 |  |  |
CIO & Units No Income | 0.7308 | -0.09 | CIO & White below 50K | 0.02481 | -0.54 |  |  |

###  4.4 Correlation between Utility Metrics and Analysis Accuracy

To determine whether there is any relationship between the general utility
metrics and analysis accuracy, we conduct several pearson correlation tests
between them. For the tests, we measure the correlation of different analysis
accuracy results with CIO and APO individually and in some cases with KL
divergence as well. We also measure the correlation between the utility
metrics CIO, APO, and KL divergence. For the classification accuracy and ad
hoc analysis, the relationship is tested between the deviation from the
original for each score and the corresponding utility metrics. Table 2 shows
the results of the correlation tests.

For the classification accuracy, we see that there exists a strong correlation
between accuracy score and CIO for all three datasets. The correlation
coefficient for all of them is above 0.8 in the negative direction. There is
also a moderate correlation between classification accuracy and APO above 90%
with the lowest coefficient being Adult with a value of -0.55. However, when
we look at the correlation between classification accuracy and KL divergence
score, there is no such consistent correlation.

In the case of ad hoc analysis, we see a mixed batch of results. For
univariate distributions and CIO relationship (i.e., Polish: CIO Raab & Units
No Income, Adult: CIO Raab & Units Income $\leq 50K$, and Avila: CIO Raab &
Units Class A), we see that there is no strong correlation for all three
datasets. However, for multivariate relationships where two or more variables
are involved, we see mixed results. For instance, Male Edu Voc/Grm and Male
$\leq 50K$ US for Polish and Adult respectively show strong correlation with
CIO with a correlation coefficient above -0.80. We also see moderate (e.g.,
Polish: CIO Raab & Female Inc 1K to 2K, Adult: CIO Raab & Female Inc 1K to 2K)
and no correlation (e.g., Polish: CIO Raab & Married alcabuse Yes) between the
multivariate analysis and CIO. Similar tests with APO above 90% show that APO
either has a similar or lower correlation coefficient than CIO.

The correlation test between the general utility metrics shows that there is a
strong correlation between CIO and APO above 90%. Nonetheless, there seems to
be no such correlation between CIO or APO above 90% with KL divergence. From
the correlation test, we find that CIO has a better correlation with analysis
accuracy than APO 90% and KL divergence. However, CIO also has limitations and
it is not possible to guarantee that a higher CIO score will ensure higher
accuracy in all types of analysis performed on the synthetic data.

###  4.5 Generation Time

For generation time, we look at APO above $90$% with the time it took to
generate $100$ datasets for each synthesizer combination for both the Simple
and the Selective methodology (Appendix Figure 5). The result of generation
time reveals that CC and D are faster overall and achieve better accuracy.
However, CP is much slower (i.e., 40 times slower in some cases) than CC and D
though they all achieve similar accuracy. For numerical variables, P performs
much faster than other synthesizers, however, the accuracy it achieves is much
worse than best performing synthesizers. While Selective is faster for most,
it is not significantly faster than Simple overall. The choice of synthesis
order can have an huge impact on the generation time in some scenarios (e.g.,
for P in Adult the slowest order takes 10 times more than the fastest). We do
not see any huge difference between proper and non-proper in terms of
generation time.

##  5 Discussion & Conclusion

In this work, we look at the impact of different parameters and synthesizers
on the utility of synthetic data. We also perform a thorough investigation of
the correlation between the utility metrics and the analysis accuracy of
synthetic data. Our investigation reveals that the choice of synthetic data
generation method, the number of datasets to release, and sometimes the
synthesis order can impact the utility of synthetic data. We find that
Kulback-Leibler divergence, a broad utility metric, does not correlate with
analysis accuracy and some narrow metrics such as confidence-interval overlaps
(CIO) show varying correlations for specific univariate and multivariate
analyses. Simply put, these metrics compare data similarity but cannot
guarantee that for any given analysis the results are similar enough.
Nevertheless, the synthesizers with the best CIOs also performed best in terms
of analysis accuracy. We found another more promising effect when comparing
the results of classification tasks by machine learning models trained on the
synthetic versus on the original data, both tested on original data. There,
CIO correlates with how similarly the machine learning models perform the
classification task, both in terms of accuracy and which records they
classified correctly or incorrectly, respectively. Machine learning models
have been shown to be susceptible to privacy threats such as membership
inference and model inversion attacks. In future work, we will therefore
investigate whether our results generalize in terms of utility and to what
extent training on synthetic data can mitigate privacy threats to machine
learning models.

## Acknowledgment

This work was partially supported by the Wallenberg AI, Autonomous Systems &
Software Program (WASP) funded by Knut & Alice Wallenberg Foundation.

## References

  * [1] Alvim, M., Chatzikokolakis, K., Palamidessi, C., Pazii, A.: Local Differential Privacy on Metric Spaces: Optimizing the Trade-off with Utility. In: 2018 IEEE 31st Computer Security Foundations Symposium (CSF). pp. 262–267. IEEE (2018)
  * [2] Czapiński, J., Panek, T.: Social Diagnosis 2011. Objective and Subjective Quality of Life in Poland. Contemporary Economics 5(3) (2011)
  * [3] Dandekar, A., Zen, R.A.M., Bressan, S.: A Comparative Study of Synthetic Dataset Generation Techniques. In: Hartmann, S., Ma, H., Hameurlain, A., Pernul, G., Wagner, R.R. (eds.) Database and Expert Systems Applications. pp. 387–395. Springer International Publishing, Cham (2018)
  * [4] Dankar, F.K., Ibrahim, M.: Fake It till You Make It: Guidelines for Effective Synthetic Data Generation. Applied Sciences 11(5), 2158 (2021)
  * [5] Dankar, F.K., Ibrahim, M.K., Ismail, L.: A Multi-Dimensional Evaluation of Synthetic Data Generators. IEEE Access (2022)
  * [6] De Stefano, C., Maniaci, M., Fontanella, F., di Freca, A.S.: Reliable Writer Identification in Medieval Manuscripts through Page Layout Features: The “Avila” Bible Case. Engineering Applications of Artificial Intelligence 72, 99–110 (2018)
  * [7] Drechsler, J.: Background on Multiply Imputed Synthetic Datasets, pp. 7–11. Springer New York, New York, NY (2011). https://doi.org/10.1007/978-1-4614-0326-5_2, https://doi.org/10.1007/978-1-4614-0326-5_2
  * [8] Drechsler, J.: Some Clarifications Regarding Fully Synthetic Data. In: Domingo-Ferrer, J., Montes, F. (eds.) Privacy in Statistical Databases. pp. 109–121. Springer International Publishing, Cham (2018)
  * [9] Drechsler, J., Bender, S., Rässler, S.: Comparing Fully and Partially Synthetic Datasets for Statistical Disclosure Control in the German IAB Establishment Panel. Trans. Data Privacy 1(3), 105–130 (dec 2008)
  * [10] Drechsler, J., Dundler, A., Bender, S., Rässler, S., Zwick, T.: A New Approach for Disclosure Control in the IAB Establishment Panel — Multiple Imputation for a Better Data Access. AStA Advances in Statistical Analysis 92(4), 439–458 (2008)
  * [11] Drechsler, J., Reiter, J.: Disclosure Risk and Data Utility for Partially Synthetic Data: An Empirical Study Using the German IAB Establishment Survey. Journal of Official Statistics 25(4), 589 (2009)
  * [12] Dua, D., Graff, C.: UCI Machine Learning Repository (2017), http://archive.ics.uci.edu/ml
  * [13] Dwork, C.: Differential Privacy: A Survey of Results. In: International conference on theory and applications of models of computation. pp. 1–19. Springer (2008)
  * [14] El Emam, K., Mosquera, L., Jonker, E., Sood, H.: Evaluating the Utility of Synthetic COVID-19 Case Data. JAMIA open 4(1), ooab012 (2021)
  * [15] Graham, P., Young, J., Penny, R.: Multiply Imputed Synthetic Data: Evaluation of Hierarchical Bayesian Imputation Models. Journal of Official Statistics 25(2), 245 (2009)
  * [16] Hittmeir, M., Ekelhart, A., Mayer, R.: On the Utility of Synthetic Data: An Empirical Evaluation on Machine Learning Tasks. In: Proceedings of the 14th International Conference on Availability, Reliability and Security. pp. 1–6 (2019)
  * [17] Hittmeir, M., Ekelhart, A., Mayer, R.: Utility and Privacy Assessments of Synthetic Data for Regression Tasks. In: 2019 IEEE International Conference on Big Data (Big Data). pp. 5763–5772. IEEE (2019)
  * [18] Karr, A.F., Kohnen, C.N., Oganian, A., Reiter, J.P., Sanil, A.P.: A Framework for Evaluating the Utility of Data Altered to Protect Confidentiality. The American Statistician 60(3), 224–232 (2006). https://doi.org/10.1198/000313006X124640, https://doi.org/10.1198/000313006X124640
  * [19] Kohavi, R., Becker, B.: Adult Dataset. UCI machine learning repository 5, 2093 (1996)
  * [20] Lee, A.: Generating Synthetic Microdata from Published Marginal Tables and Confidentialised Files. Statistics New Zealand (2009)
  * [21] Lee, J.H., Kim, I.Y., O’Keefe, C.M.: On Regression-tree-based Synthetic Data Methods for Business Data. Journal of Privacy and Confidentiality 5(1) (2013)
  * [22] Muralidhar, K., Domingo-Ferrer, J., Martínez, S.: $\epsilon$ -Differential Privacy for Microdata Releases Does Not Guarantee Confidentiality (Let Alone Utility). In: International Conference on Privacy in Statistical Databases. pp. 21–31. Springer (2020)
  * [23] Nowok, B.: Utility of Synthetic Microdata Generated Using Tree-based Methods. UNECE Statistical Data Confidentiality Work Session (2015)
  * [24] Nowok, B., Raab, G., Dibben, C.: synthpop: Bespoke Creation of Synthetic Data in R. Journal of Statistical Software, Articles 74(11), 1–26 (2016). https://doi.org/10.18637/jss.v074.i11, https://www.jstatsoft.org/v074/i11
  * [25] Nowok, B., Raab, G.M., Dibben, C.: Providing Bespoke Synthetic Data for the UK Longitudinal Studies and Other Sensitive Data with the synthpop Package for R. Statistical Journal of the IAOS 33(3), 785–796 (2017)
  * [26] Nowok, B., Raab, G.M., Dibben, C.: synthpop: Catall (2019), https://cran.r-project.org/web/packages/synthpop/synthpop.pdf#nameddest=syn.catall
  * [27] Purdam, K., Elliot, M.: A Case Study of the Impact of Statistical Disclosure Control on Data Quality in the Individual UK Samples of Anonymised Records. Environment and Planning A 39(5), 1101–1118 (2007)
  * [28] Raab, G.M., Nowok, B., Dibben, C.: Practical Data Synthesis for Large Samples. Journal of Privacy and Confidentiality 7(3), 67–97 (2016)
  * [29] Raghunathan, T.E., Reiter, J.P., Rubin, D.B.: Multiple Imputation for Statistical Disclosure Limitation. Journal of official statistics 19(1), 1 (2003)
  * [30] Reiter, J.P.: Using CART to Generate Partially Synthetic Public Use Microdata. Journal of Official Statistics 21, 441–462 (2005)
  * [31] Reje, N.: Synthetic Data Generation for Anonymization. Master’s thesis, KTH, School of Electrical Engineering and Computer Science (EECS) (2020)
  * [32] Rubin, D.B.: Multiple Imputation for Survey Nonresponse (1987)
  * [33] Rubin, D.B.: Statistical Disclosure Limitation. Journal of official Statistics 9(2), 461–468 (1993)
  * [34] Ruiz, N., Muralidhar, K., Domingo-Ferrer, J.: On the Privacy Guarantees of Synthetic Data: A Reassessment from the Maximum-knowledge Attacker Perspective. In: International Conference on Privacy in Statistical Databases. pp. 59–74. Springer (2018)
  * [35] Snoke, J., Raab, G.M., Nowok, B., Dibben, C., Slavkovic, A.: General and Specific Utility Measures for Synthetic Data. Journal of the Royal Statistical Society: Series A (Statistics in Society) 181(3), 663–688 (2018)
  * [36] Sweeney, L.: k-anonymity: A Model for Protecting Privacy. International Journal of Uncertainty, Fuzziness and Knowledge-Based Systems 10(05), 557–570 (2002)
  * [37] Taub, J., Elliot, M., Sakshaug, J.W.: The Impact of Synthetic Data Generation on Data Utility with Application to the 1991 UK Samples of Anonymised Records. Transactions on Data Privacy 13(1), 1–23 (2020)
  * [38] Wang, Z., Myles, P., Tucker, A.: Generating and Evaluating Synthetic UK Primary Care Data: Preserving Data Utility & Patient Privacy. In: 2019 IEEE 32nd International Symposium on Computer-Based Medical Systems (CBMS). pp. 126–131. IEEE (2019)
  * [39] Woo, M.J., Reiter, J., Oganian, A., Karr, A.: Global Measures of Data Utility for Microdata Masked for Disclosure Limitation. Journal of Privacy and Confidentiality 1, 111–124 (2009)

##  Appendix 0.A Appendix

Details of the Datasets

  * •

Polish Dataset: The Polish dataset is from synthpop example datasets which is
a Polish census on the quality of life in Poland [2]. We use 14 out of the 35
variables to simplify the computation. The dataset contains records with
missing values which are replaced with similar values using random sampling.
The total number of records is 5000 for the polish dataset.

  * •

Adult Dataset: The Adult dataset [19] is from UCI Repository intended for
predicting whether income exceeds 50K/yr based on census data. The dataset
contains 15 variables where we removed numerical education level which is a
redundant version of the categorical education variable. We also modify the
levels of the native country variable from 44 countries to 7 coarser ones, a
change that was required for most of the synthesizers to work. To simplify
computation, we use the first 10000 records of the dataset.

  * •

Avila Dataset: The Avila dataset [6] is also from UCI Repository intended for
predicting the copyist based on the patterns of segments of a Latin bible. The
dataset has no missing values and all 11 variables and all 10430 records are
used.

Table 3: APO above $90$% for Mean Point Estimations

m | S | P | D | CP | CC
---|---|---|---|---|---
1 | 0.00 | 0.08 | 0.00 | 0.00 | 0.00
2 | 0.90 | 0.99 | 0.85 | 0.86 | 0.86
3 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
5 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
10 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
20 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
50 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
100 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00

(a) Polish Dataset

m | S | P | D | CP | CC
---|---|---|---|---|---
1 | 0.01 | 0.00 | 0.00 | 0.01 | 0.00
2 | 0.90 | 0.50 | 0.86 | 0.83 | 0.82
3 | 1.00 | 0.50 | 0.99 | 1.00 | 0.98
4 | 1.00 | 0.50 | 1.00 | 1.00 | 1.00
5 | 1.00 | 0.50 | 1.00 | 1.00 | 1.00
6 | 1.00 | 0.50 | 1.00 | 1.00 | 1.00
7 | 1.00 | 0.50 | 1.00 | 1.00 | 1.00
8 | 1.00 | - | 1.00 | 1.00 | 1.00

(b) Adult Dataset

m | S | P | D | CP | CC
---|---|---|---|---|---
1 | 0.05 | 0.08 | 0.06 | 0.14 | 0.05
2 | 0.80 | 0.51 | 0.80 | 0.92 | 0.78
3 | 0.90 | 0.68 | 0.91 | 0.93 | 0.90
5 | 0.95 | 0.75 | 0.96 | 0.92 | 0.95
10 | 0.98 | 0.82 | 0.99 | 0.91 | 0.98
20 | 1.00 | 0.86 | 1.00 | 0.90 | 1.00
50 | 1.00 | 0.89 | 1.00 | 0.89 | 1.00
100 | 1.00 | 0.90 | 1.00 | 0.90 | 1.00

(c) Avila Dataset

Table 4: KL-Divergence Score

Dataset | metric | S | P | D | PT | DT | PO | DO | POT | DOT | PV | DV | PVT | DVT | CP | CPT | CC | CCT
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Polish | Average | 1.00 | 1.21 | 1.06 | 1.96 | 2.28 | 1.13 | 1.29 | 1.86 | 2.33 | 1.23 | 1.07 | 1.96 | 2.29 | 1.07 | 2.37 | 1.03 | 2.1
|  | S | P | D | PT | DT | PH | DH | PHT | DHT | PL | DL | PLT | DLT | CP | CPT | CC | CCT
Adult | Average | 1.00 | 9.5 | 1.04 | 10.17 | 2.17 | 4.72 | 1.05 | 5.62 | 2.21 | 16.92 | 1.05 | 17.21 | 2.12 | 1.06 | 2.48 | 1.00 | 2.03
|  | S | P | D | PT | DT | PO | DO | POT | DOT | PV | DV | PVT | DVT | CP | CPT | CC | CCT
Avila | Average | 1.00 | 6.96 | 1.08 | 7.32 | 2.31 | 8.74 | 1.6 | 8.96 | 2.78 | 5.7 | 1.08 | 5.97 | 2.28 | 1.13 | 25.36 | 1.08 | 2.31

Details of the Classification Tasks Adult and Avila dataset are primarily used
for classification tasks using machine learning. The classification task for
Adult is to classify whether the income is above or below 50K. Similarly, for
Avila it is to classify the author based on the patterns. For Polish dataset
we did not find any such common analysis. Since it is a census dataset similar
to Adult and has a income variable, we perform a similar income classification
on the dataset. For each original dataset, we train the 10 machine learning
models using 10 different synthetic datasets and take the average accuracy
score. We test the accuracy of each synthetic model using the original
dataset.

(a) Polish Dataset - Simple

(b) Adult Dataset - Simple

(c) Avila Dataset - Simple

(d) Polish Dataset - Selective

(e) Adult Dataset - Selective

(f) Avila Dataset - Selective

Figure 5: APO above 90% and Generation Time for 100 Datasets

Comparison between Simple and Selective Synthesis. For the comparison between
Simple and Selective synthesis, we looked at the APO 90% using both the
techniques with Non-Proper synthesis and $m=10$ (Appendix Table 5). The
results show that for majority of the synthesizers in all of the three
datasets simple works better than selective. For Polish, the simple method was
better regardless of the synthesizers. However, for Adult and Avila, we find
that few synthesizers performed better with selective synthesis. For Adult, CP
and CC with Selective performed slightly better than simple. Nonetheless, P
with synthesis order largest categorical variable last (PL) achieved a
significant improvement with selective. For Avila, except for P with opposite
synthesis order (PO), simple performed better than selective for all the other
synthesizers. Thus, in terms of accuracy simple seems to be a better choice.

Table 5: APO above $90$% Comparison between Simple and Selective with $m=10$

Dataset | metric | P | PSe | D | DSe | PO | POSe | DO | DOSe | PV | PVSe | DV | DVSe | CP | CPSe | CC | CCSe
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
| Average | 0.38 | 0.22 | 0.50 | 0.30 | 0.33 | 0.16 | 0.21 | 0.16 | 0.39 | 0.21 | 0.45 | 0.29 | 0.49 | 0.40 | 0.58 | 0.48
Polish | Ratio | 0.78 | 0.22 | 0.85 | 0.15 | 0.83 | 0.17 | 0.63 | 0.37 | 0.76 | 0.24 | 0.90 | 0.10 | 0.65 | 0.35 | 0.75 | 0.25
|  | P | PSe | D | DSe | PH | PHSe | DH | DHSe | PL | PLSe | DL | DLSe | CP | CPSe | CC | CCSe
| Average | 0.29 | 0.20 | 0.27 | 0.20 | 0.31 | 0.20 | 0.24 | 0.21 | 0.26 | 0.39 | 0.31 | 0.23 | 0.51 | 0.55 | 0.51 | 0.55
Adult | Ratio | 0.81 | 0.19 | 0.79 | 0.21 | 0.88 | 0.12 | 0.70 | 0.30 | 0.29 | 0.71 | 0.75 | 0.25 | 0.30 | 0.70 | 0.46 | 0.54
|  | P | PSe | D | DSe | PO | POSe | DO | DOSe | PV | PVSe | DV | DVSe | CP | CPSe | CC | CCSe
| Average | 0.10 | 0.08 | 0.14 | 0.04 | 0.06 | 0.08 | 0.18 | 0.14 | 0.11 | 0.10 | 0.24 | 0.07 | 0.23 | 0.10 | 0.24 | 0.16
Avila | Ratio | 0.68 | 0.32 | 0.91 | 0.09 | 0.37 | 0.63 | 0.61 | 0.39 | 0.58 | 0.42 | 0.96 | 0.04 | 0.88 | 0.12 | 0.74 | 0.26