Original Article [1]Equally contributing authors. Till Koebe, Department of
Economics, Freie Universität Berlin, Garystrasse 21, 14195 Berlin, Germany
<EMAIL_ADDRESS>The authors did not receive any specific funding for
this work.

# Releasing survey microdata with exact cluster locations and additional
privacy safeguards

Till Koebe  Alejandra Arias-Salazar

###### Abstract

Household survey programs around the world publish fine-granular georeferenced
microdata to support research on the interdependence of human livelihoods and
their surrounding environment. To safeguard the respondents’ privacy, micro-
level survey data is usually (pseudo)-anonymized through deletion or
perturbation procedures such as obfuscating the true location of data
collection. This, however, poses a challenge to emerging approaches that
augment survey data with auxiliary information on a local level. Here, we
propose an alternative microdata dissemination strategy that leverages the
utility of the original microdata with additional privacy safeguards through
synthetically generated data using generative models. We back our proposal
with experiments using data from the 2011 Costa Rican census and satellite-
derived auxiliary information. Our strategy reduces the respondents’ re-
identification risk for any number of disclosed attributes by 60-80% even
under re-identification attempts.

###### keywords:

Generative models, Statistical disclosure control, Geomasking, Copula,
Official statistics, Satellite imagery

##  1 Introduction

Since almost hundred years, sample surveys are dominating knowledge generation
in empirical research. The advantages of survey sampling are obvious: with an
appropriate sampling design representative results for a population can be
collected by surveying only a fraction of it. With computer assistance, the
time from collecting data to publishing results can be sped up significantly
[1]. Two trends, however, increasingly challenge the way data is collected via
surveys. On the one hand, the growing demand for fast and granular information
drives up sample size and thus costs. As a response, recent years have seen a
large amount of academic research on augmenting surveys with secondary data
from non-traditional data sources such as social networks, mobile phones or
remote sensing in order to overcome shortcomings in coverage, frequency and
granularity with applications in fields as diverse as population dynamics [2,
3], socio-demographic analysis [4, 5, 6, 7, 8], policy targeting [9, 10],
environmental mapping [11] and health research [12, 13]. This augmentation is
usually done via geographic matching, i.e. combining area-level averages [14].
Since the number of matched areas corresponds to the sample size for
subsequent supervised learning tasks, finding the smallest common geographical
denominator is essential to avoid running into small sample problems. However,
this is not always trivial as sample surveys usually provide data only for a
fraction of small geographic areas. On the other hand, digital transformations
across various sectors such as health care have led to an explosion of digital
personal data. It is the abundance of secondary data that amplifies re-
identification risks in published surveys as some of the information could be
used to link pseudoanonymized survey responses back to the actual respondents
[15, 16, 17]. Together with new privacy regulations such as the European
General Data Protection Regulation (GDPR) this calls for additional
precautionary measures to safeguard the individual’s privacy. For aggregated
data releases, the introduction of differential privacy has provided a solid
mathematical framework to manage re-identification risks independent of a
potential attacker’s capabilities or prior knowledge [18]. With regard to
microdata dissemination strategies, a common de-identification practice today
is a combination of deletion and perturbation procedures, which include
removing (unique) identifiers such as first and last name and replacing the
individual’s true location with aggregated (i.e. area-level) and randomized
information (see e.g. [19, 20, 21]).

For example, in the Demographic and Health Survey (DHS), a major global
household survey program, urban survey clusters are re-located within a 2km-
radius and rural clusters within a 5km-, sometimes even 10km-radius [22]. This
location privacy procedure has two main advantages: it does not affect the
quality of the remaining (non-spatial) survey information and it reduces the
need for other privacy safeguards, e.g. deleting or perturbing sensitive
information. However, it does not provide a similar rigorous measure for
privacy protection as already small sets of attributes can quickly increase
the chances of re-identification, even in incomplete, pseudonymous datasets
[23]. In addition, it obviously affects the utility of the published data when
it comes to matching with auxiliary data as this type of analysis relies on
the congruence of its geographic links [24, 25, 26, 27].

In that regard, advances in synthetic data generation have introduced new ways
to narrow the void between information loss and privacy protection. These
methods allow for the generation of synthetic records that resemble the real
data by reproducing relationships learned from the latter. While all
approaches have in common that they try to capture the joint distribution in
the original data, the ways to do so vastly differ. For example, [28] and [29]
use imputation processes to decompose the multidimensional joint distribution
into conditional univariate distributions. [30] and [31] use parametric models
in combination with conditional re-sampling to synthesize hierarchical
relationships. As an alternative to these fully parametric approaches, [32]
and [33] make use of classification and regression trees (CART), while more
recently, [34, 35, 23, 36, 37, 38] and others have used Bayesian networks,
Generative Adversarial Networks or copulas to capture the underlying linear
and non-linear relationships between the attributes.

The challenge for data producers is to define adequate microdata dissemination
strategies that allow users to satisfy their needs, i.e. release survey
microdata that can be used for statistical analysis and that are compatible
with other sources of information allowing to answer new and more detailed
research questions and – at the same time – it must be ensured that the
identities of the respondents are protected. In that regard, the Spatial Data
Repository of the DHS program [39] is a good example for facilitating new
types of research by combining survey microdata with geospatial covariates and
gridded interpolation surfaces. However, also those products are based on
perturbed cluster locations, thus incurring a certain information loss.

In this paper, we propose an alternative microdata dissemination strategy that
leverages the utility of the original microdata with additional privacy
safeguards through copula-generated synthetic data. Specifically, we propose
to adopt a strategy of publishing two sets of micro-level survey data: first,
the original microdata stripped of geographic identifiers below the strata-
level. Second, synthetic microdata including the true cluster locations. We
show in an experiment using Costa Rican census data from 2011 and satellite-
derived auxiliary information from WorldPop [40] that we can reduce the re-
identification risk vis-à-vis common spatial perturbation procedures, while
maintaining data utility for non-spatial analysis and improving data utility
for spatial analysis.

From the plethora of options, we choose copulas as our synthetic data
generation approach. Copulas facilitate fine-tuning as they allow us to model
the marginal distributions separately from the joint distribution. Dating back
to 1959 [41] with diverse applications since, their theoretical properties are
well understood. In comparison with alternatives like GANs, copula-based
synthetic data generation has lower computational cost [36] and it is easier
to interpret [42]. Furthermore, the procedure is in general less cumbersome,
in comparison with the steps followed by [30] to generate the synthetic
population data AAT-SILC (Artificial Austrian Statistics on Income and Living
Conditions [43]. Finally, copulas are also attractive for data producers such
as National Statistical Offices as only new nationally representative margins
are required to update the synthetic microdata file (cf. [44]). In addition,
well-documented open-source tools such as the Synthetic Data Vault [45] are
available to users with important features such as data transformation and
constraints specification.

##  2 Results

###  2.1 Geomasking to obfuscate true survey locations

We consider a survey $D$ as a random sample with sample size $n_{D}$ from a
given population $C$. Our unit of observations are individuals $i$ living
together in a household $h$. Each individual is described by a set of
attributes denoted as $X_{1},X_{2},\dots X_{d}$. While different sampling
designs are possible, we assume a commonly used complex design for larger
household surveys such as the DHS: a stratified two-stage cluster design. In
the first stage, the primary sampling units (PSUs) - usually enumeration areas
from the latest census denoted as $j$ \- are selected for each stratum $s$
with a probability proportional to (population) size $\pi_{j}$. In the second
stage, households within each selected PSU are sampled with a fixed
probability $\pi_{h|j}$. Consequently, the sampling weights defined as the
inverses of the household-level inclusion probabilities are given for each
stratum separately by:

$w_{hj}=\frac{1}{\pi_{hj}},\hskip
28.45274pt\pi_{hj}=\pi_{h|j}*\pi_{j}\quad\text{with}\quad\pi_{j}=\frac{n_{s}}{N_{s}}.$
(1)

PSUs, called clusters in the following, are geo-located as point locations
$r_{j}$ via their geographic centroids. In the following, we describe the
original survey attributes together with the original geo-locations of the
clusters as our true survey. The true survey builds our starting point for
further anonymization approaches, notably the geomasking approach and the
copula-based synthetic data generation approach. We follow the geomasking
methodology outlined in [22] by perturbing the centroids of the selected
clusters within a given larger administrative area $l$ using a rejection
sampling procedure described in Algorithm 1:

for  _$j\in D$_ do

while  _$r_{j}^{\text{masked}}\notin l_{r_{j}}$_ do

$\text{angle}\leftarrow U_{[0,360]}*\frac{\pi}{180}$ ;

/* Random displacement angle */

if  _$j$ is Urban_ then

$\text{dist}\leftarrow U_{[0,2000]}$ ;

/* Random displacement distance (in meters) for urban clusters */

end if

if  _$j$ is Rural_ then

if  _$j$ is selected as 1% of rural clusters_ then

$\text{dist}\leftarrow U_{[0,10000]}$ ;

/* Random displacement distance for 1% of rural clusters */

else

$\text{dist}\leftarrow U_{[0,5000]}$

end if

end if

$r_{x,j}^{\text{masked}}\leftarrow r_{x,j}+\text{dist}*\cos(\text{angle})$ ;

/* Displaced x-coordinate */

$r_{y,j}^{\text{masked}}\leftarrow r_{y,j}+\text{dist}*\sin(\text{angle})$ ;

/* Displaced y-coordinate */

end while

end for

Algorithm 1 Geomasked survey: DHS cluster displacement algorithm

We denote the masked point locations of the selected clusters with the
superscript masked. As the overall inclusion probability for a household is
not affected by geomasking, direct estimates and corresponding variances for
area-level aggregates $l$ and above remain the same. However, this does not
hold for area-level aggregates smaller than $l$. In the following, we describe
the original survey attributes together with the masked locations of the
clusters as our geomasked survey.

For our experiment using Costa Rican census data from 2011, point locations
for the corresponding enumeration areas are not available. Therefore, we
randomly sample them from the smallest available area denoted with $k$ \- in
our experiment the districts (at the same time the zip codes) in Costa Rica.
The zip code therefore corresponds to the smallest geographic identifier
available in the survey. Through the displacement procedure, roughly 30% of
the clusters are assigned to a new zip code, representing approx. 30% of the
sampled individuals in each simulation round.

###  2.2 Copula-based synthetic data generation

As an alternative to geomasking, we use synthetically generated survey
attributes for protecting the respondents’ privacy while keeping the true
point locations of the selected clusters. To do so, we fit a Gaussian copula
model on the original survey attributes $X_{d}$ and sample from the learned
joint distribution for each cluster individually with the originally sample
size $n_{j}$. A copula allows to describe the dependence structure - also
called association structure \- independently from the marginal distributions
(also called allocation structure). Several copula families are available. We
focus on the Gaussian copula that allows us to represent the association
structure of random variables irrespective of their true distribution through
a multivariate standard normal distribution [46]. Since we also assume the
marginals to be normally distributed, which may certainly constitute a mis-
specification for some of the variables, we regard the results rather as a
lower bound in terms of goodness-of-fit. Further, a copula is uniquely defined
only for continuous variables [47], meaning that in principle, copulas cannot
model non-continuous variables. Since socio-economic surveys are largely made
up of categorical variables, data transformation, e.g. via one-hot encoding,
is needed. In addition, we impose constraints on the marginals to account for
censoring (e.g. to avoid negative synthetic age records) or between-variable
dependencies (e.g. female and male household members need to add up to the
total household size) via rejection sampling.

Thus, the process to generate synthetic data from a survey dataset $\tilde{D}$
with transformed categorical attributes (details of the data transformation
are described in Algorithms 3 and 4 in Section 4) using a Gaussian copula
model is summarized in Algorithm 2:

Input $\tilde{D}=(\tilde{X}_{1},\dots,\tilde{X}_{d})$

Output $\tilde{S}=(\tilde{Y}_{1},\dots,\tilde{Y}_{d})$

for  _$s\in\tilde{D}$_ do

$\Psi\leftarrow$ Estimated marginal distributions of $\tilde{X}$ with
$\psi_{d}\sim\mathcal{N}(\mu,\sigma^{2})$

$\Sigma\leftarrow$ Estimated covariance matrix of $\Psi$

$U\leftarrow F(\Psi)$ ;

/* Probability integral transforms */

$C^{G}_{\Sigma}(u_{1},\dots,u_{d})\leftarrow\phi_{\Sigma}\big{(}\phi_{1}^{-1}(u_{1}),\dots,\phi_{d}^{-1}(u_{d})\big{)}$
;

/* d-dimensional Gaussian copula */

for  _$j\in\tilde{D}_{s}$_ do

for  _$i\leftarrow 1$ to $n_{D_{j}}$ _ do

while  _$\textbf{y}^{\\{i\\}}$ not meets constraints_ do

$\textbf{v}\sim\mathcal{N}(\mu,\Sigma)$

$\tilde{\textbf{y}}^{\\{i\\}}\leftarrow
F^{-1}\big{(}\phi_{1}(v_{1}),\dots,\phi_{d}(v_{d})\big{)}$ ;

/* Convert back to original space */

end while

end for

end for

end for

Algorithm 2 Synthetic survey: Copula-based synthetic data generation algorithm

$\phi_{\Sigma}$ is the cumulative distribution function (cdf) of a
multivariate normal distribution with $\mathcal{N}(\mu,\Sigma)$ and $\phi_{d}$
the cdf of a standard normal distribution. By fitting our model to the true
survey, it learns the parameters of both the allocation and association
structure, i.e. of the marginal distributions $\Psi$ and the multivariate
Gaussian copula $C_{\Sigma}^{G}(u_{1},\dots,u_{d})$. Based on these learned
relationships, new synthetic records $\tilde{\textbf{y}}^{\\{i\\}}$ are
sampled from the multivariate probability function
$c_{\Sigma}^{G}(\textbf{u})$ using the inverse probability integral transform
for each component $F_{d}^{-1}(u_{d})$ (cf. [48]). Since we sample in our
experiment for each cluster individually to ensure a synthetic cluster-level
sample size of exactly $n_{j}$, we use the parameters of a conditional
multivariate normal distribution. In case no conditions are applied, the
scenario is simplified to drawing from a multivariate standard normal
distribution. In the following, we call the synthetic attributes together with
the true cluster locations our synthetic survey. Further details about the
copula-based synthetic data generation procedure can be found in the Section 4
and in [49].

Figure 1 provides a first impression on the overall goodness-of-fit of the
three different survey datasets. Specifically, Figures 1(a) \- 1(c) show the
normalized KL divergence $Z_{KL}$ of the survey attributes y from the true
census attributes $\textbf{x}\in C$ defined as

$Z_{KL}(F_{d,k}(X_{d,k})||F_{d,k}(Y_{d,k}))=\frac{1}{1+D_{KL}(F_{d,k}(X_{d,k})||F_{d,k}(Y_{d,k}))}$
(2)

averaged across simulation runs for each attribute $d$ and zip code $k$,
respectively.

(a) True survey

(b) Geomasked survey

(c) Synthetic survey

Figure 1: Normalized Kullback-Leibler divergence (in bits) from the true
census distribution for each attribute and zip code, averaged across 100
simulation rounds.
The attributes on the y-axis are ordered by their respective number of
classes, the zip codes on the x-axis are ordered by their average sample size
across simulation rounds. Values close to one (yellow) represent little
divergence from the true census distribution and therefore indicate a high
goodness-of-fit. The number of attribute classes range from two to 111. Across
attributes and zip codes, the true survey scores best with $Z_{KL}=0.76$ in
total, followed by the synthetic survey with $Z_{KL}=0.74$ and the geomasked
survey at $Z_{KL}=0.73$.

Clearly visible is a gradient from the top left to the bottom right indicating
that the overall goodness-of-fit of the sample distributions improve the
larger the underlying sample sizes and the lower the number of classes per
categorical attribute. In addition, as expected, attributes with high levels
of non-response (visible through the white spots across the horizontal axis)
are stronger affected by sampling and anonymization compared to attributes
with little or no non-response.

###  2.3 Population uniqueness of survey respondents

To approach the utility-risk trade-off in (pseudo)-anonymized microdata, we
define two risk-related measures: the population uniqueness of the survey
respondents and the re-identification risk of a sensitive attribute in the
original data using the perturbed data. We define population uniqueness
$\Xi_{\textbf{x}}$ as the share of survey respondents being unique in the
population $C$ for a given set of attributes $\textbf{x}=(X_{1},\dots,X_{d})$:

$\Xi_{\textbf{x}}=\frac{1}{n_{D}}\sum_{i}^{n_{D}}\mathbb{1}_{\textbf{x}^{\\{i\\}}}\quad\text{with}\quad\mathbb{1}_{\textbf{x}^{\\{i\\}}}=\begin{cases}1,&\text{if}\
\textbf{x}^{\\{i\\}}\text{ unique in }C(\textbf{x})\\\
0,&\text{otherwise}\end{cases}$ (3)

Figure 2 shows how $\Xi_{\textbf{x}}$ changes with the increasing number of
attributes across 100 simulation runs. Naturally, the share constantly
increases for the true survey. For the geomasked survey, the population
uniqueness increases to a level of roughly 70%. Recalling that the only
difference between the geomasked survey and the true survey is the perturbed
zip code, the remaining 30% corresponds to the average number of survey
respondents assigned to a new zip code due to the spatial anonymization
process. Thus, not considering the zip code in the re-identification effort
would let the population uniqueness of the geomasked survey also converge
towards 1 similar to the true survey, even though requiring further knowledge
on additional attributes. For the synthetic survey, the curve remains almost
flat. The initial bump can largely be explained by the probability of a random
combination of attributes representing an actual population unique in a small
(area) sample size setting. Besides this theoretical argument, synthetic data
always provides plausible deniability to the survey respondents. Similarly to
our definition, [23] use a Gaussian copula model to estimate the empirical
likelihood of population uniqueness in incomplete datasets such as $D$ by
assuming
$\Xi_{\textbf{x}}\sim\text{Binomial}(\mathbb{1}_{\textbf{x}^{\\{i\\}}},n_{D})$
with $\textbf{x}^{\\{i\\}}i.i.d.$. While this approach is an excellent
alternative to measure the re-identification risk in micro-level survey data
when no validation data (in our experiment the 2011 Costa Rican census) is
available, it assumes that the individual records are independent and
identically distributed, which may be contestable in the presence of
hierarchical dependencies and complex sampling designs.

Figure 2: Population uniqueness across survey types.
Share of population-unique survey respondents for 100 simulation runs. The
thick lines represent the average population uniqueness across the 100
simulation runs, the thin lines individual simulation runs. In the true
survey, no attribute is perturbed. In the geomasked survey, the cluster
identifier is perturbed. In the synthetic survey, all variables but the
cluster identifier are perturbed. In the re-identified survey, the synthetic
survey is used to predict the "private" attribute – i.e. the cluster – in the
original dataset along the lines of the proposed microdata dissemination
strategy. The re-identified original survey is then used to calculate
population uniqueness. Both the re-identified and the synthetic survey provide
significant privacy gains vis-à-vis the other survey types.

Therefore, Figure 2 gives a strong indication that geomasking provides little
additional safeguards for the respondents’ privacy compared to the true survey
in the presence of third-party information on a subset of the contained
attributes. Hence, we consider an alternative microdata dissemination
strategy: instead of publishing original microdata with perturbed cluster
locations, we investigate the option of publishing two datasets - original
microdata stripped of geographic identifiers below the level of strata and
synthetic microdata with the original cluster locations. The choice is
motivated by adopting a user-centric perspective: official household survey
publications predominantly report on results up to the strata-level as results
below are usually considered not representative. Analysis that benefits from
below strata-level data often investigates proximity-related questions such as
distances to certain locations and surrounding habitat. For the former,
cluster locations are of minor importance, for the latter, however, the
perturbation procedure introduces significant levels of uncertainty to the
analysis [25]. The alternative microdata dissemination strategy obviously
conserves data utility for analysis on the representative level via the first
dataset, while the second dataset allows for the accurate capture of
proximity-related information. However, two potential shortcomings need to be
considered: first, can we use the synthetic dataset to predict the ‘private’
attribute in the original dataset, i.e. the small area identifier, thus
bypassing the privacy protection measures? Second, is the uncertainty we
introduce by synthesizing the non-spatial attributes for spatial analysis
smaller than the uncertainty from perturbing the cluster locations?

###  2.4 Risk of re-identifying private geocodes

To investigate the first shortcoming, we train a random forest model on the
small area identifier - the zip code - in the anonymized surveys for each
stratum separately. We use the trained models on the original data to predict
the zip code for each record. Finally, we evaluate our predicted label against
the original label. In addition, we compare the outcomes to randomly guessing
the correct label in order to account for the number of small areas within
each stratum. Figure 3 shows the median accuracy of the approaches across 100
simulation runs. While we are able to successfully re-construct the original
zip code in most cases for the geomasked survey, it does not work much better
for the synthetic data than for the random guess.

Figure 3: Re-identification of the zip code as private attribute in the true
survey for each stratum across 100 simulation runs.
Accuracy is measured by the share of successfully re-identified zip code
labels in the true survey. A random forest model is trained on perturbed data,
i.e. the geomasked and the synthetic survey, respectively. We evaluate the
results against the true zip code labels in the true survey and compare them
against random guesses of the private attribute.

In our experiment, only one stratum consequently hosts more than ten small
areas across all simulation runs, with one stratum hosting only two small
areas in some simulation runs, giving the random guess also a good chance to
predict correctly. Recalling from Figure 2 that roughly 70% of the displaced
clusters stay within the same zip code in the geomasked survey, even
predicting the sensitive attribute for strata hosting as little as two small
areas, average population uniqueness in the synthetic data would not exceed
much the 50/50-chance of the random guess%, thus providing better privacy
protection in the re-identified original survey than the geomasked
alternative.

###  2.5 Utility for survey augmentation

To give an indication about the utility of the different anonymization
approaches, we use a setup common in recent academic literature (cf. [4, 3,
5]): we augment the surveys with auxiliary information from geospatial (big)
data. Specifically, we construct zip code-level aggregates from gridded
satellite-derived features available from the WorldPop repository [40] and
combine them with zip code-level survey aggregates to provide predictions,
especially for areas not sampled in the survey. As our target variable, we
select the Unsatisfied Basic Needs index (Necesidades Básicas Insatisfechas
(NBI)) - a composite indicator similar to the multidimensional poverty index
(MPI) [50, 51] used as a key statistical indicator in Costa Rica. Details on
the index can be found in the Supplementary Information. We evaluate our
predictions against the census in terms of adjusted $R^{2}$, bias and the Mean
Squared Error (MSE). Figures 4(a) \- 4(c) show the performance along these
three evaluation criteria across 100 simulation runs.

(a) Adjusted R2

(b) Relative Bias

(c) MSE

(d) NBI densities

Figure 4: Performance metrics of survey-based NBI estimates on the zip code-
level.
(a) Adjusted $R^{2}$ is based on the in-sample zip codes. (b) and (c) are
based on the full sample and predictions are evaluated against the census
across 100 simulation runs. (d) compares zip code-level NBI averages for a
single simulation run.

Surprisingly, the synthetic approach not only outperforms the geomasked
survey, it also provides predictions more in line with the census results than
the true survey. A possible explanation could be that the copula approach
reduces the impact of outliers on the zip code-specific NBI sample averages.
This explanation is supported by Figure 4(d) that shows the distribution of
zip code-level NBI averages grouped into quartiles for one simulation run as
both the synthetic survey and the census showcase smaller tails in their
distributions, respectively. We run additional experiments to compare the
directly synthesized NBI and its underlying indicators with their counterparts
computed from synthetic survey variables (see Supplementary Table 2 and
Supplementary Fig. 4 in the Supplementary Information).

##  3 Discussion

In this paper, we proposed and evaluated an alternative data dissemination
strategy for micro-level survey data that improves the trade-off between
privacy risk and data utility. Specifically, we showed that by publishing two
datasets, namely the original survey data with limited geographic identifiers
and a synthetically-generated survey dataset with the true cluster locations,
re-identification risks can be reduced significantly vis-à-vis popular
geomasking approaches without incurring additional losses in terms of data
utility for survey augmentation. This could help mapping initiatives such as
WorldPop or GRID3 to improve their products as more accurate spatial data is
available. In addition, by separating the marginals from the dependence
structure, it provides data producers such as National Statistical Offices
also with a useful tool to update the respective synthetic microdata files for
the following years by updating the margins with nationally representative new
data as sub-nationally representative surveys may only be conducted every few
years. In the Supplementary Information, we further investigate the stability
of our results by alternating the experiment design.

First, while we chose the strata for the main analysis as they provide ’large-
enough’ sample sizes at the same time explicitly accounting for at least high-
level regional variation, we study in further experiments whether fitting on
smaller or larger geographic levels may better capture local variation at the
expense of running into the risk of small sample problems or vice versa.
Supplementary Fig. 1 summarizes the results for our copula model being fitted
on the whole survey, the twelve strata and the zip code-level, respectively.
It shows that by selecting the strata as our fitting level, we strike a
balance between the underlying sample size (usually the larger the better) and
capturing regional variation (usually the more disaggregated the better) both
in terms of utility and risk. In addition, by using subsets of the full
microdata for model fitting, the approach becomes computationally tractable
also for larger surveys.

Second, since generative models allow us to sample an arbitrary number of
synthetic observations, we look at the impact of the synthetic sample size on
the outcomes of the survey augmentation experiment, notably the adjusted
$R^{2}$ and a measure of confidence in the direct survey estimates of the Fay-
Herriot model (cf. Section 4.2) - the shrinkage factor $\gamma$. Supplementary
Fig. 2 shows that with an increasing sample size, $\gamma$ increases as well,
thus shifting more weight to the direct estimate. Even though intuitive as the
sampling variance naturally decreases in $n_{D}$, at some point it may become
misleading with potentially negative effects on the model performance as the
synthetic data generating process still relies on the same information
conveyed in the true survey with sample size $n_{D}$. However, in our
experiment the adjusted $R^{2}$ does not exhibit a bump, but increases
monotonically, thus hinting at little additional explanatory power of our
satellite-derived covariates vis-à-vis the area-level direct survey estimate
for the in-sample areas.

Third, since our target variable NBI is a composite indicator, we compare the
different composition levels of the synthetic NBI with the NBI constructed
from synthetic data. While the divergence measure shows an overall good fit
for the underlying indicators (see Supplementary Table 2 and Supplementary
Fig. 4), correlations are low, especially for higher-level compositions as the
dimensions or the NBI itself.

Lastly, we test alternative encoding schemes for the transformation of
categorical data. Also, we relax our assumption of the normally distributed
margins by opening up to a wider group of parametric copulas (such as beta,
gamma or uniform distributions) selected for each margin individually based on
the two-sample Kolmogorov-Smirnov (KS) statistic to study the effect of the
specification choice on the normalized KL divergence. Supplementary Fig. 3
shows that neither the encoding scheme nor the specification of the marginal
distributions have large effects on the quality of the synthetically generated
data.

Nevertheless, our approach is not without limitations. The copula-based
approach towards synthetic data generation largely fails to correctly capture
lower-level hierarchical relationships such as individuals - line numbers -
households - houses from the original data. As said before, since we see our
analysis using a naïve Gaussian copula model as providing somewhat a lower
bound for improving the utility-risk trade-off by adopting the proposed
microdata dissemination strategy vis-à-vis common geomasking approaches, there
is much room for improvement. To name a few, latent copula designs can be
considered to avoid data transformations, marginal distributions can be
modelled non-parametrically, hierarchical structures can be accounted for more
rigorously by either modelling the hierarchies separately as suggested by [20]
or by modelling the relationships explicitly. In addition, synthetic data may
- under some circumstances - leak private information, e.g. through the
generated value ranges. As a response, differentially-private implementations
of existing generative models have been proposed such as PrivBayes [35],
PrivSyn [52] and PATE-GAN [53].

##  4 Methods

###  4.1 Fitting Gaussian copulas to survey attributes

As an alternative to geomasking, we use synthetically generated survey
attributes for protecting the respondents’ privacy while keeping the true
point locations of the selected clusters. To do so, we fit a Gaussian copula
model on the original survey attributes x and sample from the learned joint
distribution for each cluster individually with the originally sample size
$n_{j}$. Therefore, consider our survey $D$, where $X_{1}$ represents a random
variable with a continuous marginal cumulative distribution function (cdf)
denoted by $F_{1}(x_{1})=P(X_{1}\leq x_{1})$. For the multivariate case, the
joint cdf can be generalized to $F_{1,\dots,d}(x_{1},\dots,x_{d})=P(X_{1}\leq
x_{1},\dots,X_{d}\leq x_{d})$.

A copula, firstly introduced in the work of [41], is a cumulative density
function with uniform marginals between [0,1]. Thus - based on Sklar´s theorem
[41] \- when all variables are continuous, the d-dimensional random vector
$X_{1},\dots,X_{d}$ can be defined in a uniform space $[0,1]^{d}$, creating a
random vector $U_{1},\dots,U_{d}$ via the probability integral transform
$u_{d}=F_{d}(x_{d})$. In this case, a unique d-dimensional copula $C(u_{d})$
exists:

$C(u_{1},\dots,u_{d})=F\big{(}F_{1}^{-1}(u_{1}),\dots,F_{d}^{-1}(u_{d})\big{)}$
(4)

As motivated in Section 2.2, we account for the fact that household surveys
largely consist of categorical variables by applying data transformation.
Among the plethora of possible encoding schemes, the most common encoding
scheme is one-hot encoding, where for each class of a categorical variable a
binary dummy variable is created [54]. A disadvantage of this option is that
it may become computationally challenging and prone to multicollinearity in
the presence of variables with a high cardinality, i.e. with a large number of
classes, since each possible class creates a new variable [55]. Interestingly,
there is – to the best of our knowledge – little comprehensive, comparative
and conclusive scientific evidence on the properties and performance of
different categorical encoding schemes. Therefore, we explore two other well-
known alternatives with more favourable computation times: ordinal and
frequency encoding. Ordinal encoding uses integers to represent each classes,
e.g. from 0 to $v$, the number of classes in a categorical variable. Assigning
an unreal order to nominal variables is the main pitfall of this alternative
[56]. Frequency encoding – as used in medical imaging [57] and similar to the
concept of term frequency in Natural Language Processing [58] – assigns an
interval in [0,1] to each class based on and ordered by its proportion of
occurrence. Then, it uses the middle point of each interval as float
representative of the respective class. Back-transformation is done by
assigning a new point to a class via the respective interval it falls into. In
this sense, this alternative conveys information of the importance of each
class [59]. Based on the results of the different encoding schemes shown in
Supplementary Fig. 3, we opt for the frequency encoding scheme in the
following. Algorithms 3 and 4 provide details on the chosen scheme.

Input $D=(X_{1},\dots,X_{d})$

Output $\tilde{D}=(\tilde{X}_{1},\dots,\tilde{X}_{d})$

for  _$X_{d}\in D$_ do

if  _$X_{d}$ is Continuous_ then

$\tilde{X}_{p}\leftarrow X_{d}$

end if

if  _$X_{d}$ is Non-continuous_ then

$\tilde{X}_{q}\leftarrow\mathrm{T}(X_{d})$

end if

end for

$\tilde{D}\leftarrow(\tilde{X}_{p},\tilde{X}_{q})\quad\forall\quad p\in
P\text{ and }q\in Q$

Algorithm 3 Transform categorical variables

Input $\tilde{S}=(\tilde{Y}_{1},\dots,\tilde{Y}_{d})$

Output $S=(Y_{1},\dots,Y_{d})$

for  _$\tilde{Y}_{d}\in\tilde{S}$_ do

if  _$\tilde{Y}_{d}$ is not indexed as variable in $Q$_ then

$Y_{p}\leftarrow\tilde{Y}_{d}$

end if

if  _$\tilde{Y}_{d}$ is indexed as variable in $Q$_ then

$Y_{q}\leftarrow T^{-1}(\tilde{Y}_{d})$

end if

end for

$S\leftarrow(Y_{p},Y_{q})\quad\forall\quad p\in P\text{ and }q\in Q$

Algorithm 4 Back-transform frequency encoded variables

Thus, the d-dimensional Gaussian copula $C_{\Sigma}^{G}(\textbf{u})$ is
defined as the cdf of a multivariate normal distribution
$\mathcal{N}(\mu,\Sigma)$ with $\Sigma\in\mathbb{R}^{dxd}$ represented on the
unit cube $[0,1]^{d}$:

$C_{\Sigma}^{G}(u_{1},\dots,u_{d})=\phi_{\Sigma}\big{(}\phi^{-1}(u_{1}),\dots,\phi^{-1}(u_{d})\big{)}$
(5)

The density of a Gaussian copula is then defined as:

$c_{\Sigma}^{G}(\textbf{u})=\frac{1}{\sqrt{\text{det
}\Sigma}}\text{exp}\Big{(}-\frac{1}{2}\phi^{-1}(\textbf{u})^{T}\cdot(\Sigma^{-1}-I)\cdot\phi^{-1}(\textbf{u})\Big{)}$
(6)

with $\textbf{u}\in[0,1]^{d}$, $I\in\mathbb{R}^{dxd}$ being the identity
matrix, and $\phi^{-1}$ being the inverse cumulative distribution function of
a standard normal distribution. $\Sigma$ is a positive semi-definite
covariance matrix that we estimate based on Pearson’s correlation coefficient
$\rho$ [34].

As noted in Section 2.2, we sample for each cluster individually with a sample
size of $n_{j}$. While rejection sampling could be an option for ensuring only
synthetic rows with the respective cluster identifier are selected, it proves
computationally inefficient. With copulas being multivariate cdfs, we
introduce conditions instead. Hence, we sample from a multivariate normal
distribution conditional on cluster $j$. Thus, our transformed dataset
$\tilde{D}$ with one conditional variable becomes
$\tilde{D}=(\tilde{\textbf{X}}_{a}|\tilde{\textbf{X}}_{b})$ with
$\tilde{\textbf{X}}_{a}:=\tilde{X}_{1},\dots,\tilde{X}_{d-1}$ being the
transformed attributes to be synthesized and
$\tilde{\textbf{X}}_{b}:=\tilde{X}_{d}$ being the transformed cluster
identifier. The parameters of the respective multivariate normal distributions
are thus partitioned into:

$\tilde{D}=\begin{bmatrix}\tilde{\textbf{X}}_{a}\\\
\tilde{\textbf{X}}_{b}\end{bmatrix}\text{,}{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}x}\mu=\begin{bmatrix}\mu_{a}\\\
\mu_{b}\end{bmatrix}\text{and}{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}x}\Sigma=\begin{bmatrix}\Sigma_{aa}&\Sigma_{ab}\\\
\Sigma_{ba}&\Sigma_{bb}\end{bmatrix}$ (7)

with $\mu_{a}\in\mathbb{R}^{d-1}$ and $\mu_{b}\in\mathbb{R}^{1}$ and
$\Sigma_{aa}\in\mathbb{R}^{(d-1)\times(d-1)}$,
$\Sigma_{ab}\in\mathbb{R}^{(d-1)\times 1}$,
$\Sigma_{ba}\in\mathbb{R}^{1\times(d-1)}$, and
$\Sigma_{bb}\in\mathbb{R}^{1\times 1}$ being the means and positive semi-
definite covariance matrices, respectively. Following Algorithm 2, the
parameters of our estimated marginal distributions $\Psi$ and of the copula
$C_{\Sigma}^{G}(\textbf{u})$ need to be adapted to mirror the conditionality
such that $\Psi_{a|b}(\tilde{\textbf{X}}_{a}|\tilde{\textbf{X}}_{b})$ and
$C_{\Sigma}^{G}(\textbf{u}_{a}|\textbf{u}_{b})$. Consequently, we sample from
$\sim\mathcal{N}(\bar{\mu},\bar{\Sigma}$) with:

$\bar{\mu}=\mu_{a}+\Sigma_{ab}\Sigma^{-1}_{bb}(X_{b}-\mu_{b})\in\mathbb{R}^{d-1}$
(8)

and

$\bar{\Sigma}=\Sigma_{aa}-\Sigma_{ab}\Sigma^{-1}_{bb}\Sigma_{ba}\in\mathbb{R}^{(d-1)\times(d-1)}.$
(9)

We iterate the copula-based fitting and sampling procedure for every stratum
separately as it allows to better capture sub-national variation using
representative sub-samples and as it proves computationally more tractable.
For sampling designs with varying household- or individual-level inclusion
probabilities (e.g. in the DHS, women - in comparison to men - are usually
oversampled), [20] suggests to sample a synthetic population and re-iterate
the sampling procedure to produce valid synthetic sampling weights. As in our
design sampling weights are identical across households for a given PSU due to
the systematic sampling approach in the second stage, the original sampling
weights remain valid. The virtue in our model choice is the relative
simplicity, little requirements in terms of ex-ante knowledge about the
individual distributions $X_{d}$ and its computational efficiency. For further
experiments on the robustness and sensitivity of our modelling choices, we
refer to the Supplementary Information.

###  4.2 Area-level survey augmentation methods

Survey data can be augmented with the use of area-level models, e.g. the Fay-
Herriot model [60] by linking direct estimators gathered from survey data to
relevant auxiliary information. Both, direct estimators, and auxiliary data
are aggregated on $k=1,\dots,D$ areas. Traditionally, these auxiliary
covariates $\textbf{x}_{k}$ are obtained from recent censuses, administrative
records or other geospatial (big) data sources. In this paper, we make use of
satellite imagery features as area-level covariates. The Fay-Herriot is a two-
level model, the first part is composed by the sampling model:

$\hat{\theta}^{\text{Dir}}_{k}=\theta_{k}+e_{k},\hskip 14.22636pte_{k}\sim
N(0,\sigma^{2}_{e_{k}}),$ (10)

where the sampling error is represented by $e_{k}$ and
$\hat{\theta}^{\text{Dir}}_{k}$ is the direct estimator of $\theta_{k}$ (e.g.
sample mean). The linking model provides the second part, where relevant area-
level covariates are considered:

$\theta_{k}=\textbf{x}^{\prime}_{k}\hat{\beta}+u_{k}.$ (11)

Here, the random area effects $u_{k}$ are assumed to be independent with mean
0 and variance $\sigma^{2}_{u}$. The empirical best linear unbiased predictor
(EBLUP) estimator is given by:

$\hat{\theta}^{\text{FH}}_{k}=\gamma_{k}\hat{\theta}^{Dir}_{k}+(1-\gamma_{k})\textbf{x}^{\prime}_{k}\hat{\beta}=\textbf{x}^{\prime}_{k}\hat{\beta}+\hat{u}_{k},{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color<EMAIL_ADDRESS>(12)

with
$\gamma_{k}=\frac{\hat{\sigma}^{2}_{u}}{\hat{\sigma}^{2}_{u}+\hat{\sigma}^{2}_{e_{k}}}$
denoting the shrinkage factor for each area $k$. The parameter estimates of
this model can be obtained via maximum likelihood (ML) or restricted ML
(REML). Note that the shrinkage factor allows to weight in favor of the direct
estimator when sampling variances are small; on the contrary the synthetic
estimator $\textbf{x}^{\prime}_{k}\hat{\beta}$ receives more weight when the
sampling variance is larger. Results on an experiment studying the sensitivity
of the shrinkage factor and adjusted $R^{2}$ for varying synthetic sample
sizes are shown in Supplementary Fig. 2. Further details on the Fay-Herriot
model can be found in [61].

## References

  * [1] Granello, D. H. & Wheaton, J. E.  Online data collection: Strategies for research.  J Couns Dev. 82, 387–393 (2004).  URL https://onlinelibrary.wiley.com/doi/10.1002/j.1556-6678.2004.tb00325.x.
  * [2] Stevens, F. R., Gaughan, A. E., Linard, C. & Tatem, A. J.  Disaggregating census data for population mapping using Random forests with remotely-sensed and ancillary data.  PLoS ONE 10, e0107042 (2015).  URL https://doi.org/10.1371/journal.pone.0107042.
  * [3] Leasure, D. R., Jochem, W. C., Weber, E. M., Seaman, V. & Tatem, A. J.  National population mapping from sparse survey data: A hierarchical Bayesian modeling framework to account for uncertainty.  Proc. Natl. Acad. Sci. U.S.A. 117, 24173–24179 (2020).  URL https://www.pnas.org/doi/abs/10.1073/pnas.1913050117.
  * [4] Pokhriyal, N. & Jacques, D. C.  Combining disparate data sources for improved poverty prediction and mapping.  Proc. Natl. Acad. Sci. U.S.A. 114, E9783–E9792 (2017).  URL https://www.pnas.org/doi/full/10.1073/pnas.1700319114.
  * [5] Schmid, T., Bruckschen, F., Salvati, N. & Zbiranski, T.  Constructing sociodemographic indicators for national statistical institutes by using mobile phone data: estimating literacy rates in Senegal.  J. R. Stat. Soc. Ser. A Stat. Soc. 180, 1163–1190 (2017).  URL https://rss.onlinelibrary.wiley.com/doi/abs/10.1111/rssa.12305.
  * [6] Subash, S. P., Kumar, R. R. & Aditya, K. S.  Satellite data and machine learning tools for predicting poverty in rural India.  Agric. Econ. Res. Rev. 31, 231–240 (2018).  URL https://ageconsearch.umn.edu/record/284254.
  * [7] Fatehkia, M., Coles, B., Ofli, F. & Weber, I.  The relative value of facebook advertising data for poverty mapping.  Proc. Int. AAAI Conf. Web Soc. Media 14, 934–938 (2020).  URL https://ojs.aaai.org/index.php/ICWSM/article/view/7361.
  * [8] Chi, G., Fang, H., Chatterjee, S. & Blumenstock, J. E.  Microestimates of wealth for all low- and middle-income countries.  Proc. Natl. Acad. Sci. U.S.A. 119, e2113658119 (2022).  URL https://www.pnas.org/doi/abs/10.1073/pnas.2113658119.  2104.07761.
  * [9] Blumenstock, J. E.  Estimating economic characteristics with phone data.  AEA Pap. Proc. 108, 72–76 (2018).  URL https://www.aeaweb.org/articles?id=10.1257/pandp.20181033.
  * [10] Aiken, E., Bellue, S., Karlan, D., Udry, C. & Blumenstock, J. E.  Machine learning and phone data can improve targeting of humanitarian aid.  Nature 2022 603, 864–870 (2022).  URL https://www.nature.com/articles/s41586-022-04484-9.
  * [11] Grace, K. et al. Integrating Environmental Context into DHS Analysis While Protecting Participant Confidentiality: A New Remote Sensing Method.  Popul. Dev. Rev. 45, 197–218 (2019).  URL https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6446718/.
  * [12] Brown, M. E., Grace, K., Shively, G., Johnson, K. B. & Carroll, M.  Using satellite remote sensing and household survey data to assess human health and nutrition response to environmental change.  Popul. Environ. 36, 48–72 (2014).  URL https://doi.org/10.1007/s11111-013-0201-0.
  * [13] Arambepola, R. et al. Spatiotemporal mapping of malaria prevalence in madagascar using routine surveillance and health survey data.  Sci. Rep. 10, 18129 (2020).  URL https://doi.org/10.1038/s41598-020-75189-0.
  * [14] Koebe, T.  Better coverage, better outcomes? Mapping mobile network data to official statistics using satellite imagery and radio propagation modelling.  PLoS ONE 15, e0241981 (2020).  URL https://dx.plos.org/10.1371/journal.pone.0241981.
  * [15] Armstrong, M. P., Rushton, G. & Zimmerman, D. L.  Geographically masking health data to preserve confidentiality.  Stat. Med. 18, 497–525 (1999).  URL {https://onlinelibrary.wiley.com/doi/10.1002/%28SICI%291097-0258%2819990315%2918%3A5%3C497%3A%3AAID-SIM45%3E3.0.CO%3B2-%23}.
  * [16] Kroll, M. & Schnell, R.  Anonymisation of geographical distance matrices via Lipschitz embedding.  Int. J. Health Geogr. 15, 1–14 (2016).  URL https://doi.org/10.1186/s12942-015-0031-7.
  * [17] West, B. T. et al. Establishing Infrastructure for the Use of Big Data to Understand Total Survey Error, chap. 21, 457–485 (John Wiley & Sons, Ltd, 2017).  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/9781119041702.ch21.
  * [18] Dwork, C.  Differential privacy: A survey of results.  In Theory and Applications of Models of Computation. TAMC 2008. Lecture Notes in Computer Science, vol. 4978, 1–19 (Springer, Berlin, Heidelberg, 2008).  URL https://doi.org/10.1007/978-3-540-79228-4_1.
  * [19] Andrés, M. E., Bordenabe, N. E., Chatzikokolakis, K. & Palamidessi, C.  Geo-indistinguishability: Differential privacy for location-based systems.  In Proc. 2013 ACM SIGSAC Conf. Comput. Commun. Secur., 901–914 (2013).  URL https://doi.org/10.1145/2508859.2516735.
  * [20] Templ, M.  Statistical disclosure control for microdata (Springer, 2017).
  * [21] de Jonge, E. & de Wolf, P.-P.  sdcSpatial: Statistical Disclosure Control for Spatial Data (2019).  URL https://CRAN.R-project.org/package=sdcSpatial.  R package version 0.1.1.
  * [22] Burgert, C. R., Colston, J., Roy, T. & Zachary, B.  Geographic Displacement Procedure and Georeferenced Data Release Health Surveys. DHS Spatial Analysis Reports.  Tech. Rep. 7, ICF International, USAID, Calverton, Maryland, USA (2013).  URL https://dhsprogram.com/pubs/pdf/SAR7/SAR7.pdf.
  * [23] Rocher, L., Hendrickx, J. M. & de Montjoye, Y. A.  Estimating the success of re-identifications in incomplete datasets using generative models.  Nat. Commun. 10, 3069 (2019).  URL https://doi.org/10.1038/s41467-019-10933-3.
  * [24] Elkies, N., Fink, G. & Bärnighausen, T.  “Scrambling” geo-referenced data to protect privacy induces bias in distance estimation.  Popul. Environ. 37, 83–98 (2015).  URL https://doi.org/10.1007/s11111-014-0225-0.
  * [25] Warren, J. L., Perez-Heydrich, C., Burgert, C. R. & Emch, M. E.  Influence of demographic and health survey point displacements on distance-based analyses.  Spatial Demography 4, 155–173 (2016).  URL https://doi.org/10.1007/s40980-015-0014-0.
  * [26] Blankespoor, B., Croft, T., Dontamsetti, T., Mayala, B. & Murray, S.  Spatial anonymization: Guidance note prepared for the Inter-Secretariat working group on household surveys.  Tech. Rep., UN Inter-secretariat Working Group on Household Surveys Task Force on Spatial Anonymization in Public-Use Household Survey Datasets (2021).  URL https://unstats.un.org/iswghs/task-forces/documents/Spatial_Anonymization_Report_submit01272021_ISWGHS.pdf.
  * [27] Hunter, L. M. et al. Working toward effective anonymization for surveillance data: innovation at South Africa’s Agincourt Health and Socio-Demographic Surveillance Site.  Popul. Environ. 42, 445–476 (2021).  URL https://doi.org/10.1007/s11111-020-00372-4.
  * [28] 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.  Adv. Stat. Anal. 92, 439–458 (2008).  URL https://doi.org/10.1007/s10182-008-0090-1.
  * [29] Heldal, J. & Iancu, D.-C.  Synthetic data generation for anonymization purposes. Application on the Norwegian Survey on living conditions/EHIS.  In Joint UNECE/Eurostat Work Session on Statistical Data Confidentiality (2019).  URL https://unece.org/fileadmin/DAM/stats/documents/ece/ces/ge.46/2019/mtg1/SDC2019_S1_Norway_Heldal_Iancu_AD.pdf.
  * [30] Alfons, A., Kraft, S., Templ, M. & Filzmoser, P.  Simulation of close-to-reality population data for household surveys with application to EU-SILC.  Stat. Methods Appt. 20, 383–407 (2011).  URL https://doi.org/10.1007/s10260-011-0163-2.
  * [31] Templ, M., Meindl, B., Kowarik, A. & Dupriez, O.  Simulation of synthetic complex data: The R package simPop.  J. Stat. Softw. 79, 1–38 (2017).  URL https://www.jstatsoft.org/index.php/jss/article/view/v079i10.
  * [32] Reiter, J. P.  Using CART to generate partially synthetic public use microdata.  J. Off. Stat. 21, 441–462 (2005).  URL https://www.scb.se/contentassets/ca21efb41fee47d293bbee5bf7be7fb3/using-cart-to-generate-partially-synthetic-public-use-microdata.pdf.
  * [33] Wang, H. & Reiter, J. P.  Multiple imputation for sharing precise geographies in public use data.  Ann. Appl. Stat. 6, 229–252 (2012).  URL https://doi.org/10.1214/11-AOAS506.
  * [34] Li, H., Xiong, L. & Jiang, X.  Differentially private synthesization of multi-dimensional data using copula functions.  In Advances in database technology: proceedings. International conference on extending database technology, 475–486 (2014).  URL https://doi.org/10.5441/002/edbt.2014.43.
  * [35] Zhang, J., Cormode, G., Procopiuc, C. M., Srivastava, D. & Xiao, X.  Privbayes: Private data release via bayesian networks.  ACM Trans. Database Syst. 42, 1–41 (2017).  URL https://doi.org/10.1145/3134428.
  * [36] Sun, Y., Cuesta-Infante, A. & Veeramachaneni, K.  Learning vine copula models for synthetic data generation.  Proc. AAAI Conf. Artif. Intell. 33, 5049–5057 (2019).  URL https://doi.org/10.1609/aaai.v33i01.33015049.
  * [37] Torkzadehmahani, R., Kairouz, P. & Paten, B.  Dp-cgan: Differentially private synthetic data and label generation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (2019).  URL https://openaccess.thecvf.com/content_CVPRW_2019/html/CV-COPS/Torkzadehmahani_DP-CGAN_Differentially_Private_Synthetic_Data_and_Label_Generation_CVPRW_2019_paper.html.
  * [38] Xu, L., Skoularidou, M., Cuesta-Infante, A. & Veeramachaneni, K.  Modeling tabular data using conditional GAN.  Adv. Neural Inf. Process. Syst. 32, 1–11 (2019).  URL https://proceedings.neurips.cc/paper/2019/hash/254ed7d2de3b23ab10936522dd547b78-Abstract.html.
  * [39] ICF.  The DHS Program Spatial Data Repository. (2022).  URL https://spatialdata.dhsprogram.com/home/.
  * [40] WorldPop.  Global High Resolution Population Denominators Project (2018).  URL www.worldpop.org.
  * [41] Sklar, A.  Fonctions de répartition à n dimensions et leurs marges.  Publ. Inst. Statist. Univ. Paris 8, 229–231 (1959).
  * [42] Kamthe, S., Assefa, S. & Deisenroth, M.  Copula flows for synthetic data generation.  arXiv preprint arXiv:2101.00598 (2021).
  * [43] Alfons, A. et al. Synthetic Data Generation of SILC Data. Research Project Report WP6, D6.2.  Tech. Rep., The AMELI Project (2011).  URL https://www.uni-trier.de/fileadmin/fb4/projekte/SurveyStatisticsNet/Ameli_Delivrables/AMELI-WP6-D6.2-240611.pdf.
  * [44] Koebe, T., Arias-Salazar, A., Rojas-Perilla, N. & Schmid, T.  Intercensal updating using structure-preserving methods and satellite imagery.  J. R. Stat. Soc. Ser. A Stat. Soc. (2022).  URL https://onlinelibrary.wiley.com/doi/full/10.1111/rssa.12802.
  * [45] MIT Data To AI Lab.  The synthetic data vault (SDV).  https://sdv.dev/ (2022).
  * [46] Patki, N., Wedge, R. & Veeramachaneni, K.  The synthetic data vault.  In 2016 IEEE International Conference on Data Science and Advanced Analytics (DSAA), 399–410 (2016).
  * [47] Jeong, B., Lee, W., Kim, D.-S. & Shin, H.  Copula-based approach to synthetic population generation.  PloS ONE 11, e0159496 (2016).  URL https://doi.org/10.1371/journal.pone.0159496.
  * [48] Janke, T., Ghanmi, M. & Steinke, F.  Implicit generative copulas.  In Ranzato, M., Beygelzimer, A., Dauphin, Y., Liang, P. & Vaughan, J. W. (eds.) Adv. Neural Inf. Process. Syst., vol. 34, 26028–26039 (Curran Associates, Inc., 2021).  URL https://proceedings.neurips.cc/paper/2021/file/dac4a67bdc4a800113b0f1ad67ed696f-Paper.pdf.
  * [49] Nelsen, R. B.  An introduction to copulas (Springer Science & Business Media, 2007).
  * [50] Méndez, F. & Bravo, O.  Costa Rica Mapas de Pobreza 2011.  Tech. Rep., INEC Costa Rica, San José, Costa Rica (2011).  URL https://www.inec.cr/sites/default/files/documentos/pobreza_y_presupuesto_de_hogares/pobreza/metodologias/documentos_metodologicos/mepobrezacenso2011-01.pdf.pdf.
  * [51] Alkire, S., Kanagaratnam, U. & Suppa, N.  The Global Multidimensional Poverty Index (MPI) 2019\. OPHI MPI Methodological Note 47.  Tech. Rep., Oxford Poverty and Human Development Initiative, University of Oxford (2019).  URL https://www.ophi.org.uk/wp-content/uploads/OPHI_MPI_MN_47_2019_vs2.pdf.
  * [52] Zhang, Z. et al. PrivSyn: Differentially private data synthesis.  In Proceedings of the 30th USENIX Security Symposium, 929–946 (2021).  URL https://www.usenix.org/system/files/sec21fall-zhang-zhikun.pdf.  2012.15128.
  * [53] Jordon, J., Yoon, J. & Van Der Schaar, M.  PATE-GaN: Generating synthetic data with differential privacy guarantees.  In International Conference on Learning Representations (2019).  URL https://openreview.net/forum?id=S1zk9iRqF7.
  * [54] Benali, F., Bodénès, D., Labroche, N. & de Runz, C.  MTCopula: Synthetic complex data generation using copula.  In 23rd International Workshop on Design, Optimization, Languages and Analytical Processing of Big Data (DOLAP), 51–60 (2021).  URL https://hal.archives-ouvertes.fr/hal-03188317.
  * [55] Bourou, S., El Saer, A., Velivassaki, T. H., Voulkidis, A. & Zahariadis, T.  A Review of Tabular Data Synthesis Using GANs on an IDS Dataset.  Information 12, 375–389 (2021).  URL https://www.mdpi.com/2078-2489/12/9/375.
  * [56] Jiang, D., Lin, W. & Raghavan, N.  A novel framework for semiconductor manufacturing final test yield classification using machine learning techniques.  IEEE Access 8, 197885–197895 (2020).
  * [57] Mansfield, P. & Maudsley, A. A.  Medical imaging by nmr.  Br. J. Radiol. 50, 188–194 (1977).
  * [58] Aizawa, A.  An information-theoretic perspective of tf–idf measures.  Inf. Process. Manag. 39, 45–65 (2003).
  * [59] Sabharwal, N. & Agrawal, A.  Introduction to Word Embeddings, 41—63 (Apress, Berkeley, CA, 2021).  URL https://doi.org/10.1007/978-1-4842-6664-9_3.
  * [60] Fay, R. E. & Herriot, R. A.  Estimates of income for small places: an application of James-Stein procedures to census data.  J. Am. Stat. Assoc. 74, 269–277 (1979).  URL https://doi.org/10.1080/01621459.1979.10482505.
  * [61] Rao, J. N. K. & Molina, I.  Small Area Estimation (Wiley, New York, 2015), 2nd edn.
  * [62] Instituto Nacional de Estadistica y Censos.  X Censo Nacional de Población y VI de Vivienda. Catálogo central de datos (2022).  URL http://sistemas.inec.cr/pad5/index.php/catalog/113.

##  5 Supplementary information

###  5.1 Data description

As our reference dataset in this project, we use data from Costa Rica –
notably the Xth Population and VIth Housing Census of Costa Rica, 2011 (Censo
Nacional de población y Viviendas de Costa Rica 2011) – to produce three
different data file types: First, we draw survey samples from a census
population using a stratified two-stage cluster sample design without applying
any statistical disclosure control mechanisms. We use these survey samples
(called true surveys in the study) as starting point for creating file types
two and three: By re-assigning clusters to new zip codes based on the
displacement algorithm described in Algorithm 1 , we perturb the zip code
identifier in the true surveys, thereby creating the geomasked surveys. Again
based on the true surveys, we apply the copula-based synthetic data generation
algorithm described in Algorithm 2 to generate synthetic data for each
attribute except the zip code, which keeps it original structure. In addition,
in order to test the robustness of our specifications, we create additional
datasets with alternating data generating process designs. The censuses are
carried out every ten years by the national statistic office of Costa Rica
(INEC) and collect information of people, households, and dwellings on topics
such as access to education, employment, social security, technology necessary
for the planning, execution, and evaluation of public policies [50].

Administratively, Costa Rica had in 2011 four disaggregation levels: two
zones, six planning regions, 81 cantons and 473 districts (municipalities).
The sampling design used for the main National Household Survey (Encuesta
Nacional de Hogares, ENAHO) specifies twelve strata - each planning region
divided by urban and rural areas. In this case, the strata coincide with the
study domains. For our experiment, we use a 10% random sample of the original
2011 census, which can be obtained from [62] as a pseudo-population. The
smallest geographical information available in this dataset are the 473
districts. In the first stage, we select districts as our PSUs for each
stratum separately with a selection probability proportional to population
size. In the second stage, we select a minimum of 10 households in each PSU by
using simple random sampling without replacement. PSUs with less than 10
households are discarded from this procedure, affecting roughly 4% of all
PSUs.

$N_{C}$ | $n_{D}$ |  # of PSUs in $C$ |  # of PSUs in $D$ | # of attributes
---|---|---|---|---
427830 | [7638; 11914] | 767 | 123 | 106
Table 1: Descriptive statistics on the census-derived data across 100
simulation runs

As auxiliary information, we use covariates derived from satellite imagery.
Specifically, we use features derived from satellite imagery provided by [40]
in our survey augmentation setup. The advantages of using satellite imagery
here are five-fold: Data with virtually global coverage at high spatial
resolutions for frequent time intervals on human-made impact provided in a
structured format enables us to extract covariates for all administrative
areas in Costa Rica at the time of the census. Therefore, we can use area-
level survey augmentation (cf. Methods Section) to provide estimates,
especially for areas not covered by the respective survey. WorldPop data are
provided in the tagged image file format (TIFF) with a pixel representing
roughly a 100m $\times$ 100m grid square in an open data repository under
CC4.0 licence ([40]). Pixel values are aggregated to the administrative areas
of Costa Rica via their centroids. Specifically, we generate area-level
averages for the distances to different types of natural areas (e.g.
cultivated, woody-tree, and shrub areas, coastlines etc.) and to
infrastructure such as roads and waterways, the intensity of night-time
lights, topographic information and information on the presence of human
settlements.

###  5.2 Sensitivity of copula vis-à-vis geographic fitting level

In order to study the effect of the geographic level on the copula modelling
performed for synthetic data generation, we run Algorithm 2 on the whole
survey (‘Country), the twelve strata (‘Strata’) and the roughly 110 zip code
areas (‘Zip Code’), respectively. Results are provided in Supplementary Fig.
5. It appears that fitting the copula model on the whole survey limits the
ability of the approach to capture regional variations. On the other hand,
model fitting on the zip code-level does neither increase the re-
identification risk of the zip code identifier as a private attribute and nor
affect the overall prediction performance of the outcome variable, hinting at
overfitting not being a problem on that level. Striking a balance between
underlying sample size and a certain level of disaggregation shows better
results. Also, it allows to scale computations to settings with larger samples
and more attributes.

(a) Adjusted $R^{2}$

(b) Relative Bias

(c) MSE

(d) Re-identification Risk

(e) Population uniqueness across attributes

Figure 5: Evaluation metrics for different geographical copula fitting levels.
(a) - (c) The copula model is fitted on the whole survey (‘Country’), for each
of the twelve strata (‘Strata’) and for each of the roughly 110 zip codes
(‘Zip Code’) separately. As a reference, the metrics for the geomasked and the
true survey are provided as well. (d) The accuracy to successfully re-identify
the zip code as a private attribute in the original data using a random forest
model trained on synthetic data across fitting levels remains similar. (e) The
share of population-unique survey respondents is virtually not affected by the
copula fitting level.

###  5.3 Effects of synthetic sample size on prediction outcomes

Generative models can be used to create synthetic samples of an arbitrary size
regardless the amount of underlying data. While the advantages of that are
similar to those of other resampling procedures such as bootstrapping (i.e. to
estimate the precision of the sample statistics or to perform cross-
validation), it can also mislead modelling approaches that ‘borrow strength’
from auxiliary data by overestimating the strength of the synthetic direct
estimates eventually resulting in losses of explanatory power of the model. In
our survey augmentation setup, the shrinkage factor $\gamma$ indicates whether
final estimates rather rely on the direct estimates from the synthetic survey
or on the satellite-derived covariates for the in-sample predictions depending
on the sampling variance. Supplementary Figure 6 shows that larger sample
sizes lead to increasing gamma values (via decreasing sampling variances of
the direct estimator), however, not incurring losses in the goodness-of-fit of
our estimation model. This hints at the fact that the contribution of the
auxiliary information to the explanatory power of the model for the in-sample
predictions is negligible.

Figure 6: Sensitivity of model performance on changes in synthetic sample
size.
Samples are drawn from a synthetic population. The synthetic population is
generated using the copula-based approach described in the Results Section.
Sample sizes are determined by the sampling rate (shown on the x-axis).
Results are evaluated against the true census population. The shrinkage factor
$\gamma$ is averaged across zip codes. The thick lines represent the metric
averages across the 100 simulation runs, the thin lines individual simulation
runs.

###  5.4 Choosing marginal distributions & encoding schemes

As already mentioned in the Section 2.2, assuming normally-distributed margins
may represent a misspecification of the true univariate distribution of
$X_{d}$. In addition, computationally tractable alternatives to one-hot
encoding exist. We compare two different ways to model the marginal
distributions together with two different encoding schemes. The results are
presented in Supplementary Fig. 7. Measured by the normalized KL divergence
averaged across 100 simulation runs, frequency encoding produces slightly
better goodness-of-fit of the synthetic data ($Z_{KL}=0.74$ for frequency
encoding versus $Z_{KL}=0.72$ for ordinal encoding with gaussian marginals).
Surprisingly, the naïve assumption of normally distributed marginals
outperforms the KS-based parametric marginals with $Z_{KL}=0.74$ and
$Z_{KL}=0.70$, respectively.

Figure 7: Effect of encoding schemes and marginal distribution choice on the
overall goodness-of-fit of the synthetic data measured by the normalized KL
divergence $Z_{KL}$ (in bits).
The attributes on the y-axis are ordered by their respective number of
classes, the zip codes on the x-axis are ordered by their average sample size
across simulation rounds. Values close to one (yellow) represent little
divergence from the true census distribution and therefore indicate a high
goodness-of-fit.

###  5.5 Detailed analysis of the NBI as composite indicator

The NBI is a composite indicator computed from approx. 20 underlying survey
variables grouped into four dimensions (i.e. access to decent housing (Acceso
albergue digno), access to a healthy life (Acceso a vida saludable), access to
knowledge (Acceso al conocimiento) and access to other goods and services
(Acceso a otros bienes y servicios)) using 19 indicators in total. All
indicators and dimensions are binary (yes/no). An identified need in one of
the indicators leads to a positive needs status in higher dimensions. The
sensitivity for false positives is thus assumed to be high for the NBI as a
small change (e.g. one year age difference) in one of the 19 underlying
variables can turn a NBI-negative to a NBI-positive survey respondent.

Generally, two strategies for computed indicators exist to create synthetic
counterparts: a) directly synthesize the computed indicators or b) re-
construct the indicator based on synthetic survey variables. While the former
is more likely to reflect the original distribution, it may not be
consistently decomposable into its underlying indicators; vice-versa holds for
the latter. The strength of these effects are largely determined by the
complexity and sensitivity of the composite indicator and the overall
goodness-of-fit of the synthetic data. Thus, if both approaches produce
similar compositions, it can be regarded as a strong indication that the
underlying synthetic data also successfully captures relationships across
multiple variables in the dataset, not only the composite index. Supplementary
Table 2 shows that this not fully holds for the NBI.

Indicators | # of indicators | Pearson’s $\rho$ | $Z_{KL}$ | Incidence
---|---|---|---|---
1.x | 5 | 0.42 | 0.99 | 100
Dimension 1 |  | 0.24 | 0.98 | 647
2.x | 5 | 0.22 | 0.98 | 85
Dimension 2 |  | 0.19 | 0.98 | 455
3.x | 2 | 0.02 | 0.89 | 507
Dimension 3 |  | 0.02 | 0.84 | 1845
4.x | 7 | 0.02 | 0.99 | 60
Dimension 4 |  | 0.03 | 1.00 | 622
Composite NBI | 19 | 0.07 | 0.97 | 3253
Table 2: Relationship between synthetic and computed NBI indicators across 100
simulation runs.
Indicator-level results (e.g. 1.x) are averaged across indicators. The
incidence describes the average number of respondents across 100 simulated
surveys with unsatisfied needs in the respective indicator/dimension.

Although the overall number of survey respondents with unsatisfied needs are
captured with a high accuracy as measured by the normalized KL divergence
$Z_{KL}$ for binary data, the NBI status on the individual level strongly
diverges following Pearson’s $\rho$ (cf. Supplementary Table 2). Supplementary
Figure 8 shows that the lack of linear correlation is mainly due to improperly
captured relationships in the underlying variables than in the synthetic NBI
as the former is outperformed by the latter for survey augmentation expressed
in terms of adjusted $R^{2}$, bias and MSE. However, it remains on par with
the geomasked survey at lower privacy risks.

(a) Adjusted R2

(b) Relative Bias

(c) MSE

(d) NBI densities

Figure 8: Performance of the synthetic vs. computed composite NBI.
(a) - (c) show of the different survey types in our survey augmentation
experiment across 100 simulation runs. (d) shows the densities of the
composite NBI by quartiles for one simulation run.


# TIP: A Trust Inference and Propagation Model in Multi-Human Multi-Robot
Teams

Yaohui Guo  0000-0002-7064-9249 University of MichiganAnn ArborMichiganUSA
<EMAIL_ADDRESS>,  X. Jessie Yang  0000-0001-6071-0387 University of
MichiganAnn ArborMichiganUSA<EMAIL_ADDRESS>and  Cong Shi
0000-0003-3564-3391 University of MichiganAnn ArborMichiganUSA
<EMAIL_ADDRESS>
(2023)

###### Abstract.

Trust has been identified as a central factor for effective human-robot
teaming. Existing literature on trust modeling predominantly focuses on dyadic
human-autonomy teams where one human agent interacts with one robot. There is
little, if not no, research on trust modeling in teams consisting of multiple
human agents and multiple robotic agents. To fill this research gap, we
present the trust inference and propagation (TIP) model for trust modeling in
multi-human multi-robot teams. We assert that in a multi-human multi-robot
team, there exist two types of experiences that any human agent has with any
robot: direct and indirect experiences. The TIP model presents a novel
mathematical framework that explicitly accounts for both types of experiences.
To evaluate the model, we conducted a human-subject experiment with 15 pairs
of participants ($N=30$). Each pair performed a search and detection task with
two drones. Results show that our TIP model successfully captured the
underlying trust dynamics and significantly outperformed a baseline model. To
the best of our knowledge, the TIP model is the first mathematical framework
for computational trust modeling in multi-human multi-robot teams.

team of teams, multi-operator multi-autonomy (MOMA)

††journalyear: 2023††copyright: acmlicensed††conference: Companion of the 2023
ACM/IEEE International Conference on Human-Robot Interaction; March 13–16,
2023; Stockholm, Sweden††booktitle: Companion of the 2023 ACM/IEEE
International Conference on Human-Robot Interaction (HRI ’23 Companion), March
13–16, 2023, Stockholm, Sweden††price: 15.00††doi:
10.1145/3568294.3580164††isbn: 978-1-4503-9970-8/23/03

##  1\. Introduction

Trust has been identified as one central factor for effective human-robot
teaming (Sheridan, 2016; Yang et al., 2021, 2023). Despite research efforts
over the past thirty years, existing literature predominantly focuses on
dyadic human-robot teams where one human agent interacts with one robot
(National Academies of Sciences, Engineering, and Medicine, 2022). There is
little, if not no, research on trust modeling in teams consisting of multiple
human agents and multiple robots.

Consider a scenario where two human agents, $x$ and $y$, and two robots, $A$
and $B$, are to perform a task. The four agents are allowed to form sub-teams
to enhance task performance (e.g., maximizing throughput and/or minimizing
task completion time). For instance, they could initially form two dyadic
human-robot teams to complete the first part of the task, merge to complete
the second part, and split again with a different configuration to complete
the third part of the task, and so on (see Fig 1).

Figure 1. Four agents can form sub-teams. In Part 1, human $x$ and robot $A$
form a dyad, and human $y$ and robot $B$ form a dyad. In part 2, two dyads
merge. In part 3, human $x$ and robot $B$ form a dyad, and human $y$ and robot
$A$ form a dyad.

In this scenario, we assert that there exist two types of experiences that a
human agent has with a robot: _direct and indirect experiences_. Direct
experience, by its name, means that a human agent’s interaction with a robot
is by him-/her-self; indirect experience means that a human agent’s
interaction with a robot is mediated by another party. Fig. 2 illustrates Part
_3_ of the task shown in Fig. 1. Human $x$ works directly with robot $B$
(i.e., direct experience). Even though there is no direct interaction between
$x$ and $A$ in part 3, we postulate that $x$ could still update his or her
trust in $A$ by learning $y$’s experience with $A$, i.e., $y$’s direct
experience with $A$ becomes $x$’s indirect experience. $y$’s trust in $A$
_propagates_ to $x$.

Figure 2. An arrow points from a trustor to a trustee, representing the trust
$t^{\text{trustor},\text{trustee}}$. Human $x$ updates her trust in robot $B$
via direct experience. Even though $x$ does not have direct interaction with
$A$, $x$ could still update his or her trust toward $A$ through a third party,
$y$.

Under the direct and indirect experience framework, prior work on trust
modeling in dyadic human-robot teams can be regarded as examining how _direct_
experience influences a person’s trust in a robot. In multi-human multi-robot
teams, we postulate that both direct and indirect experiences drive a human
agent’s trust in a robot.

In this study, we develop the Trust Inference and Propagation (TIP) model for
multi-human multi-robot teams, explicitly accounting for the direct and
indirect experiences a human agent may have with a robot. We examine trust
dynamics under the TIP framework and prove theoretically that trust converges
after repeated (direct and indirect) interactions. To evaluate the proposed
TIP model, we conducted a human-subject experiment with 15 pairs of
participants ($N=30$). Each pair worked with two drones to perform a threat
detection task for 15 sessions. We compared the TIP model (i.e., accounts for
both the direct and indirect experiences) and a direct-experience-only model
(i.e., only accounts for the direct experience a human agent has with a
robot). Results show that the TIP model successfully captured people’s trust
dynamics with a significantly smaller root-mean-square error (RMSE) compared
to the direct-experience-only model. To the best of our knowledge, the
proposed TIP model is the first mathematical framework for computational trust
modeling in multi-human multi-robot teams.

##  2\. Related Work

Several computational trust models in dyadic human-robot teams exist (Chen et
al., 2018; Xu and Dudek, 2015; Guo and Yang, 2021; Pippin and Christensen,
2014). Notably, Xu and Dudek (Xu and Dudek, 2015) proposed the online
probabilistic trust inference model (OPTIMo) utilizing Bayesian networks to
estimate human trust from automation’s performance and human behavior. Hu et
al. (Hu et al., 2016) proposed to classify trust or distrust based on
electroencephalography (EEG) and galvanic skin response (GSR) signals. Soh et
al. (Soh et al., 2020) modeled trust as a latent dynamic function for
predicting human trust in robots across different tasks. Guo and Yang (Guo and
Yang, 2021) and Bhat et al. (Bhat et al., 2022) proposed to model trust as a
Beta random variable and predicted trust value in a Bayesian framework. For a
detailed review, refer to (Kok and Soh, 2020).

Even though the amount of research on trust modeling in multi-human multi-
robot teams is extremely limited, some inspiration can be drawn from research
on reputation/trust management. Central to the reputation management system is
a propagation mechanism that allows a buyer to obtain the
reputation/trustworthiness values of a seller, especially when the buyer had
no prior transactions with the seller (Hendrikx et al., 2015). Examples
include the Beta reputation system (Josang and Ismail, 2002) and the FIRE
trust management model (Huynh et al., 2006). Variants of such propagation
mechanisms can be found in several multi-agent systems, such as e-commerce
(Cen et al., 2020) and social networks (Fan et al., 2019).

##  3\. Mathematical Model

We present the detailed TIP model. Our objective is to develop a fully
computational trust propagation model that works in general human-robot
interaction settings.

###  3.1. Assumptions

We make two major assumptions. First, we assume that each human agent
communicates trust as a single-dimensional value (Guo and Yang, 2021; Guo et
al., 2021). Second, we assume that human agents are cooperative, i.e., they
report their trust truthfully to their human teammates.

###  3.2. Proposed Model

Trust as a Beta random variable. We define a trust value as a real number in
$[0,1]$, where 0 stands for “[do] not trust at all” and 1 stands for “trust
completely”. We take a probabilistic view to model trust as in (Guo and Yang,
2021). At time $k$, the trust $t_{k}^{a,b}$ that a human agent $a$ has towards
another agent $b$ follows a Beta distribution, i.e.,
$t_{k}^{a,b}\sim\operatorname{Beta}\left(\alpha_{k}^{a,b},\beta_{k}^{a,b}\right)$,
where $\alpha_{k}^{a,b}$ and $\beta_{k}^{a,b}$ are the positive and negative
experiences $a$ had about $b$ up to time $k$, respectively, $k=0,1,2,\dots$.
When $k=0$, $\alpha_{0}^{a,b}$ and $\beta_{0}^{a,b}$ are the prior experiences
that $a$ had before any interaction with $b$. The expected trust is given by
$\mu_{k}^{a,b}=\alpha_{k}^{a,b}/\left(\alpha_{k}^{a,b}+\beta_{k}^{a,b}\right)$.
Here we note that $t_{k}^{a,b}$ is the self-reported trust given by the agent
$a$, which has some randomness due to subjectivity, while $\mu_{k}^{a,b}$ is
the expected trust determined by the experiences.

Trust update through direct experience. Similar to (Guo and Yang, 2021), we
update the direct trust experience at time ${k}$ by setting

(1) $\begin{aligned} \alpha_{k}^{a,b}&=\alpha_{k-1}^{a,b}+s^{a,b}\cdot
p_{k}^{b}\\\
\beta_{k}^{a,b}&=\beta_{k-1}^{a,b}+f^{a,b}\cdot\overline{p}_{k}^{b}\end{aligned}\text{.}$

Here $p_{k}^{b}$ and $\overline{p}_{k}^{b}$ are the measurements of $b$’s
success and failure during time $k$, respectively; $s^{a,b}$ and $f^{a,b}$ are
$a$’s unit experience gains with respect to success or failure of $b$,
respectively. We require $s^{a,b}$ and $f^{a,b}$ to be positive to ensure that
cumulative experiences are non-decreasing. The updated trust $t_{k}^{a,b}$
follows the distribution
$\operatorname{Beta}(\alpha_{k}^{a,b},\beta_{k}^{a,b})$.

Trust update through indirect experience propagation. Let $x$ and $y$ denote
two human agents and let $A$ denote a robot agent, as illustrated in Fig. 2.
At time $k$, $y$ communicates his or her trust $t_{k}^{y,A}$ on $A$ with $x$,
and then $x$ updates his or her indirect trust experience by

(2) $\begin{aligned} \alpha_{k}^{x,A}&=\alpha_{k-1}^{x,A}+\hat{s}^{x,A}\cdot
t_{k}^{x,y}\cdot\left[t_{k}^{y,A}-t_{k-1}^{x,A}\right]^{+}\\\
\beta_{k}^{x,A}&=\beta_{k-1}^{x,A}+\hat{f}^{x,A}\cdot
t_{k}^{x,y}\cdot\left[t_{k-1}^{x,A}-t_{k}^{y,A}\right]^{+}\end{aligned}\text{,}$

where the superscript ‘$+$’ means taking the positive part of the
corresponding number, i.e., $t^{+}=\max\\{0,t\\}$ for a real number $t$, and
$t_{k}^{x,A}\sim\operatorname{Beta}(\alpha_{k}^{x,A},\beta_{k}^{x,A})$.

The intuition behind this model is that $x$ needs to reason upon
$t_{k}^{y,A}$, i.e., $y$’s trust towards $A$. First, $x$ compares $y$’s trust
$t_{k}^{y,A}$ with his or her previous trust $t_{k-1}^{x,A}$. Let $\Delta
t:=t_{k}^{y,A}-t_{k-1}^{x,A}$ be the difference. If $\Delta t\geq 0$, $x$
gains positive indirect experience about $A$, which amounts to the product of
the trust difference $\Delta t$, a coefficient $\hat{s}^{x,A}$, and a
discounting factor $t_{k}^{x,y}$, i.e., $x$’s trust on $y$; if $\Delta t<0$,
then $x$ gains negative indirect experience about $A$, which is defined
similarly.

###  3.3. Parameter Inference

The proposed model characterizes a human agent’s trust on a robot by six
parameters. For instance, the parameter of $x$ on robot $A$, which is defined
as
$\theta^{x,A}=\left(\alpha_{0}^{x,A},\beta_{0}^{x,A},s^{x,A},f^{x,A},\hat{s}^{x,A},\hat{f}^{x,A}\right)$,
including $x$’s prior experiences $\alpha_{0}^{x,A}$ and $\beta_{0}^{x,A}$,
the unit direct experience gains $s_{A}^{x}$ and $f_{A}^{x}$, and the unit
indirect experience gains $\hat{s}_{A}^{x}$ and $\hat{f}_{A}^{x}$. We denote
the indices of $x$’s direct and indirect trust updates with $A$ up to time $k$
as $D_{k}$ and $\overline{D}_{k}$, respectively. Then, we can compute
$\alpha_{k}^{x,A}$ and $\beta_{k}^{x,A}$, according to Eqs. (1) and (2), as

(3) $\begin{aligned} \alpha_{k}^{x,A}=&\alpha_{0}^{x,A}+s^{x,A}\sum_{j\in
D_{k}}p_{j}^{A}+\hat{s}\sum_{j\in\overline{D}_{k}}t_{j}^{x,y}\left[t_{j}^{y,A}-t_{j-1}^{x,A}\right]^{+}\\\
\beta_{k}^{x,A}=&\beta_{0}^{x,A}+f^{x,A}\sum_{j\in
D_{k}}\overline{p}_{j}^{A}+\hat{f}\sum_{j\in\overline{D}_{k}}t_{j}^{x,y}\left[t_{j-1}^{x,A}-t_{j}^{y,A}\right]^{+}\end{aligned}.$

The optimal parameter $\theta_{*}^{x,A}$ maximizes the log likelihood function

(4)
$H\left(\theta^{x,A}\right):=\sum_{k=0}^{K}\log\operatorname{Beta}\left(t_{k}^{x,A}\middle|\alpha_{k}^{x,A},\beta_{k}^{x,A}\right),$

where $\alpha_{k}^{x,A}$ and $\beta_{k}^{x,A}$ are defined in Eq. (3).

We note that
$\log\operatorname{Beta}(t_{k}^{x,A}|\alpha_{k}^{x,A},\beta_{k}^{x,A})$ is
concave in $\theta^{x,A}$ by the composite rule that the function is concave
in $(\alpha_{k}^{x,A},\beta_{k}^{x,A})$ and $\alpha_{k}^{x,A}$ and
$\beta_{k}^{x,A}$ are non-decreasing linear functions of $\theta^{x,A}$.
Consequently, $H(\theta^{x,A})$ is concave in $\theta^{x,A}$ because it is a
summation of several concave functions. Therefore, we can run the gradient
descent method to compute the optimal parameters.

##  4\. Human-Subject Study

We conducted a human-subject experiment to evaluate the proposed model. The
experiment, inspired by (Yang et al., 2017), simulated a threat detection
task, where two human agents work with two smart drones to search for threats
at multiple sites.

###  4.1. Participants

A total of $N=30$ participants (average age = 25.3 years, SD = 4.3 ages, 16
females, 14 males) with normal or corrected-to-normal vision formed 15 teams
and participated in the experiment. Each participant received a base payment
of $15 and a bonus of up to $10 depending on their team performance. To
promote cooperation between a pair of players, team performance instead of
individual performance was used to calculate the bonus.

###  4.2. Experimental Task and Design

In the experiment, a pair of participants performed a simulated threat
detection task with two assistant drones for $K=15$ sessions on two separate
desktop computers. At each session, each participant was assigned one drone
and worked on the detection tasks. After the session, they were asked to
report their trust in each drone and their trust in their human teammate. For
clarity, we named the two drones $A$ and $B$ and colored them in red and blue,
respectively; and we denoted the participants $x$ and $y$. A trust rating is
denoted as $t^{a,b}_{k}$, where the superscript $a\in\\{x,y\\}$ stands for the
trustor, the superscript $b\in\\{x,y,A,B\\}$ stands for the trustee, and the
subscript $k$ is the session index. For example, $t^{x,A}_{2}$ is person $x$’s
trust in drone $A$ after the 2nd session. The range of a trust rating is
$[0,1]$, where 0 stands for “(do) not trust at all” and 1 stands for “trust
completely”. The flow of the experimental task is illustrated in Fig. 3(a).

((a)) Flow of experimental task

((b)) Task interface. The drone will highlight the potential threat in bright
red. The participant is asked to click the ‘Danger’ button if a threat is
present and to click the ‘Clear’ button otherwise.

Figure 3. Experimental task and design

Initial trust rating: At the start, each participant gave their initial trust
in the two drones based on their prior experience with automation/robots.
Additionally, they gave their initial trust in each other. These trust ratings
were indexed by 0, e.g., $x$’s initial trust rating on $A$ was denoted as
$t^{x,A}_{0}$.

Robot assignment: At each session, each participant was randomly assigned one
drone as his or her assistant robot.

Detection task: Each session consisted of 10 locations to detect. As shown in
Fig. 3(b), four views were present at each location. If a threat, which
appeared like a combatant, was in any of the views, the participant should
click the ‘Danger’ button; otherwise, they should click the ‘Clear’ button.
Meanwhile, his or her drone would assist and highlight a view if the drone
detected a threat there. In addition, a 3-second timer was set for each
location. If a participant did not click either button before the timer
counted down to zero, the testbed would move to the next location
automatically. After all the 10 locations, an end-of-session screen was shown,
displaying how many correct choices the participant and the drone had made in
the current session. Correct choices mean correctly identifying threats or
declaring ‘Clear’ within 3 seconds.

Figure 4. Fitting results. Red curves are for drone $A$ while blue curves are
for drone $B$. The solid lines are the participants’ self-reported trust,
while the dashed lines are the expected trust value predicted by the model.
The shaded areas indicate the 90% probability interval of the Beta
distribution at each session. The index $i$-$j$ stands for the $j$th
participant in the $i$the group.

Trust rating: After each session, participants reported three trust values.
First, each participant updated his or her trust in the drone s/he just worked
with, i.e., through direct experience. Next, each participant submitted and
communicated their trust score to their human teammate. After that, each
participant updated his or her trust in the drone the human teammate just
worked with (i.e., the other drone) and his or her trust in the human
teammate. After participants completed all 15 sessions, the experiment ended.

###  4.3. Experimental Procedure

Before the experiment, each participant signed a consent form and filled out a
demographic survey. Two practice sessions were provided, wherein a practice
drone was used to assist the participants. The participants were told that the
practice drone differed from the two drones used in the real experiment. After
the experiment started, the assignment of drones was randomized in each group.
Specifically, we assigned drone $A$ with equal change to either participant
and then assigned drone $B$ to the other participant. The threat detection
accuracy of the practice drone, drone $A$, and drone $B$ were set to 80%, 90%,
and 60%, respectively.

##  5\. Results and Discussion

We use the gradient descent method in Sec. 3.3 to compute the optimal
parameters $\theta^{p_{i},A}_{*}$ and $\theta^{p_{i},B}_{*}$ for each
participant $p_{i}$. The fitting results are shown in Fig. 4. We calculate the
performance measurements of drone $A$ at session $k$ as $p_{k}^{A}=A_{k}/10$
and $\overline{p}_{k}^{A}=1-p_{k}^{A}$, where $A_{k}$ is the number of correct
choices drone $A$ made in the $k$th session; and we define $p_{k}^{B}$ and
$\overline{p}_{k}^{B}$ similarly. To measure the performance of the model, we
calculate the fitting error at each session for each participant as
$e_{k}^{p_{i},R}=|\mu_{k}^{p_{i},R}-t_{k}^{p_{i},R}|$, $R\in\\{A,B\\}$, where
$t_{k}^{p_{i},R}$ is the participant’s self-reported trust while
$\mu_{k}^{p_{i},R}$ is the expected trust defined in section 3.2 and computed
based on $\theta_{*}^{p_{i},R}$; and, we calculate the root-mean-square error
(RMSE) between the ground truth and the expected trust value as

${\textstyle\text{RMSE}^{R}=\left[\frac{1}{N}\sum_{i=1}^{N}\frac{1}{K+1}\sum_{k=0}^{K}\left(e_{k}^{p_{i},R}\right)^{2}\right]^{1/2},}$

for $R\in\\{A,B\\}$. The RMSE results for the TIP model are
$\text{RMSE}^{A}=0.057$ and $\text{RMSE}^{B}=0.082$.

Fig. 4 shows the fitting results of the TIP model. The shaded regions indicate
the 90% confidence interval of the Beta distribution at each session. We
observe that for most participants, such as 7-2 and 10-2, the proposed TIP
model can accurately fit the trust curve with a narrow confidence interval;
but for some other participants, such as 5-2 and 8-1, the model cannot fit the
trust curve due to trust oscillation. However, in the latter case, the fitted
curve has a similar trend with the ground truth and can cover most data points
with the 90% confidence interval.

For comparison, we consider a direct-update-only model that only accounts for
the direct experience a human agent has with a robot. The direct-update-only
model is equivalent to the TIP model with zero unit indirect experience gains,
i.e., $\hat{s}^{x,A}=\hat{f}^{x,A}=0$. We recompute the model parameters for
the direct-update-only model, and the corresponding RMSE errors are
$\text{RMSE}^{\prime A}=0.085$ and $\text{RMSE}^{\prime B}=0.107$.
Furthermore, we compare each participant’s fitting error
$\bar{e}^{p_{i},R}:={1}/({K+1})\sum_{k=0}^{K}e_{k}^{p_{i},R}$ of the TIP model
($A$: $0.044\pm 0.037$; $B$: $0.069\pm 0.045$) and that of the direct-update-
only model ($A$: $0.075\pm 0.041$; $B$: $0.095\pm 0.051$) using a paired-
sample t-test. Results show that the former is significantly smaller than the
latter, with $t(29)=-6.18,p<.001$ for drone $A$, and $t(29)=-7.31$, $p<.001$
for drone $B$.

##  6\. Acknowledgement

This work is supported by the National Science Foundation under Grant No.
2045009 and the Air Force Office of Scientific Research under Grant No.
FA9550-20-1-0406.

## References

  * (1)
  * Bhat et al. (2022) Shreyas Bhat, Joseph B. Lyons, Cong Shi, and X. Jessie Yang. 2022\.  Clustering Trust Dynamics in a Human-Robot Sequential Decision-Making Task.  _IEEE Robotics and Automation Letters_ 7, 4 (2022), 8815–8822.  https://doi.org/10.1109/LRA.2022.3188902
  * Cen et al. (2020) Yukuo Cen, Jing Zhang, Gaofei Wang, Yujie Qian, Chuizheng Meng, Zonghong Dai, Hongxia Yang, and Jie Tang. 2020\.  Trust Relationship Prediction in Alibaba E-Commerce Platform.  _IEEE Transactions on Knowledge and Data Engineering_ 32, 5 (2020), 1024–1035.  https://doi.org/10.1109/TKDE.2019.2893939
  * Chen et al. (2018) Min Chen, Stefanos Nikolaidis, Harold Soh, David Hsu, and Siddhartha Srinivasa. 2018.  Planning with Trust for Human-Robot Collaboration. In _Proceedings of the 2018 ACM/IEEE International Conference on Human-Robot Interaction_ (Chicago, IL, USA) _(HRI ’18)_. Association for Computing Machinery, New York, NY, USA, 307–315.  https://doi.org/10.1145/3171221.3171264
  * Fan et al. (2019) Wenqi Fan, Yao Ma, Qing Li, Yuan He, Eric Zhao, Jiliang Tang, and Dawei Yin. 2019\.  Graph Neural Networks for Social Recommendation. In _The World Wide Web Conference on - WWW ’19_. ACM Press, San Francisco, CA, USA, 417–426.  https://doi.org/10.1145/3308558.3313488
  * Guo et al. (2021) Yaohui Guo, Cong Shi, and X. Jessie Yang. 2021.  Reverse psychology in trust-aware human-robot interaction.  _IEEE Robotics and Automation Letters_ 6, 3 (2021), 4851–4858.  https://doi.org/10.1109/LRA.2021.3067626
  * Guo and Yang (2021) Yaohui Guo and X. Jessie Yang. 2021.  Modeling and Predicting Trust Dynamics in Human–Robot Teaming: A Bayesian Inference Approach.  _International Journal of Social Robotics_ 13 (2021), 1899–1909.  https://doi.org/10.1007/s12369-020-00703-3
  * Hendrikx et al. (2015) Ferry Hendrikx, Kris Bubendorfer, and Ryan Chard. 2015\.  Reputation systems: A survey and taxonomy.  _J. Parallel Distrib. Comput._ 75 (2015), 184–197.  https://doi.org/10.1016/j.jpdc.2014.08.004
  * Hu et al. (2016) Wan Lin Hu, Kumar Akash, Neera Jain, and Tahira Reid. 2016\.  Real-Time Sensing of Trust in Human-Machine Interactions.  _IFAC-PapersOnLine_ 49, 32 (Jan. 2016), 48–53.  https://doi.org/10.1016/j.ifacol.2016.12.188
  * Huynh et al. (2006) Trung Dong Huynh, Nicholas R. Jennings, and Nigel R. Shadbolt. 2006.  An integrated trust and reputation model for open multi-agent systems.  _Autonomous Agents and Multi-Agent Systems_ 13, 2 (2006), 119–154.  https://doi.org/10.1007/s10458-005-6825-4
  * Josang and Ismail (2002) Audun Josang and Roslan Ismail. 2002.  The beta reputation system. In _Proceedings of the 15th Bled Electronic Commerce Conference_ , Vol. 5. 2502–2511.
  * Kok and Soh (2020) Bing Cai Kok and Harold Soh. 2020.  Trust in robots: Challenges and opportunities.  _Current Robotics Reports_ 1, 4 (2020), 297–309.  https://doi.org/10.1007/S43154-020-00029-Y
  * National Academies of Sciences, Engineering, and Medicine (2022) National Academies of Sciences, Engineering, and Medicine. 2022.  _Human-AI Teaming: State-of-the-Art and Research Needs_.  The National Academies Press, Washington, DC.  https://doi.org/10.17226/26355
  * Pippin and Christensen (2014) Charles Pippin and Henrik Christensen. 2014.  Trust modeling in multi-robot patrolling. In _2014 IEEE International Conference on Robotics and Automation (ICRA)_. 59–66.  https://doi.org/10.1109/ICRA.2014.6906590
  * Sheridan (2016) Thomas B. Sheridan. 2016\.  Human–Robot Interaction: Status and Challenges.  _Human Factors_ 58, 4 (2016), 525–532.  https://doi.org/10.1177/0018720816644364
  * Soh et al. (2020) Harold Soh, Yaqi Xie, Min Chen, and David Hsu. 2020\.  Multi-task trust transfer for human–robot interaction.  _The International Journal of Robotics Research_ 39, 2-3 (2020), 233–249.  https://doi.org/10.1177/0278364919866905
  * Xu and Dudek (2015) Anqi Xu and Gregory Dudek. 2015.  OPTIMo. In _Proceedings of the Tenth Annual ACM/IEEE International Conference on Human-Robot Interaction - HRI ’15_. ACM Press, New York, New York, USA, 221–228.  https://doi.org/10.1145/2696454.2696492
  * Yang et al. (2023) X. Jessie Yang, Yaohui Guo, and Christoper Schemanske. 2023\.  From Trust to Trust Dynamics: Combining Empirical and Computational Approaches to Model and Predict Trust Dynamics In Human-Autonomy Interaction.  In _Human-Automation Interaction: Transportation_ , Vincent G. Duffy, Steven J. Landry, John D. Lee, and Neville Stanton (Eds.). Springer International Publishing, Cham, 253–265.  https://doi.org/10.1007/978-3-031-10784-9_15
  * Yang et al. (2021) X. Jessie Yang, Christopher Schemanske, and Christine Searle. 2021.  Toward Quantifying Trust Dynamics: How People Adjust Their Trust After Moment-to-Moment Interaction With Automation.  _Human Factors_ (2021), 00187208211034716\.  https://doi.org/10.1177/00187208211034716
  * Yang et al. (2017) X. Jessie Yang, Vaibhav V. Unhelkar, Kevin Li, and Julie A. Shah. 2017.  Evaluating Effects of User Experience and System Transparency on Trust in Automation. In _Proceedings of the 2017 ACM/IEEE International Conference on Human-Robot Interaction - HRI ’17_. ACM, New York, NY, USA, 408–416.  https://doi.org/10.1145/2909824.3020230