Boyam Fabrice Yameogo

Pierre-Olivier Vandanjon

Pascal Gastineau

Pierre Hankach

Generating a two-layered synthetic population for French municipalities: Results and evaluation of four synthetic reconstruction methods

Keywords: Synthetic Population Generation, Multi-level, Microsimulation, Simultaneous Control

This article describes the generation of a detailed two-layered synthetic population of households and individuals for French municipalities. Using French census data, four synthetic reconstruction methods associated with two probabilistic integerization methods are applied. The paper offers an in-depth description of each method through a common framework. A comparison of these methods is then carried out on the basis of various criteria. Results show that the tested algorithms produce realistic synthetic populations with the most efficient synthetic reconstruction methods assessed being the Hierarchical Iterative Proportional Fitting and the relative entropy minimization algorithms. Combined with the Truncation Replication Sampling allocation method for performing integerization, these algorithms generate household-level and individual-level data whose values lie closest to those of the actual population.

Introduction

Agent-Based Models (ABMs) have grown in popularity since the 1990's and are now applied in a range of sectors: healthcare [START_REF] Tomintz | The geography of smoking in leeds: estimating individual smoking rates and the implications for the location of stop smoking services[END_REF][START_REF] Edwards | Simobesity: combinatorial optimisation (deterministic) model[END_REF], economic policy evaluation [START_REF] Avram | The distributional effects of fiscal consolidation in nine eu countries. Tech. rep., Euromod working[END_REF]Sutherland & Figari 2013), geography (O'Sullivan 2008), and transport [START_REF] Kickhöfer | Pricing local emission exposure of road traffic: An agent-based approach[END_REF][START_REF] Hörl | A first look at bridging discrete choice modeling and agent-based microsimulation in matsim[END_REF]. These models require comprehensive data on the demographic and socioeconomic characteristics of individuals and households. However, for privacy reasons, there is no complete dataset can be compiled on the socio-demographic characteristics of individuals at a small geographic scale. To perform a microsimulation, one necessary step consists of generating a "synthetic population" that is representative of the actual population. During this process, the characteristics of (all) the individuals within in a given study area are normally inferred from the characteristics of individuals in the sample from that area, as well as from the marginal distributions (aggregate data). The resulting synthetic population is a simplified microscopic representation of the actual population because only the variables of interest are to be reproduced [START_REF] Chapuis | A brief review of synthetic population generation practices in agent-based social simulation[END_REF]. Most approaches developed to generate a synthetic population have focused on deriving either individual-centered or household-centered populations. For example, Iterative Proportional Fitting (IPF), which by far is the most widely used algorithm for generating a synthetic population, does not yield populations linking households and individuals (and thus controlling at both levels); use of this algorithm, outputs a population yet does not link individual characteristics to household information. In many cases, this absence of a link clearly constitutes a shortcoming since an individual's decision depends on both his/her characteristics and family situation, which highlights the need to generate synthetic populations that take into account not only the individual level but also the household-level information. This article evaluates and tests, for the French case, the most appropriate methods to generate a two-layered population capable of satisfying the following conditions: maintaining the hierarchical structure of the data by associating individual and household variables in the most optimal manner; reflecting the heterogeneity of the distribution of households and individuals across geographic areas [START_REF] Münnich | On the simulation of complex universes in the case of applying the german microcensus[END_REF]; reproducing the interdependencies among agents in the same household [START_REF] Sun | A hierarchical mixture modeling framework for population synthesis[END_REF]; possessing the ability to fit with aggregate data.

Many methods serve to generate a synthetic population of individuals and households; they differ depending on the assumptions made or according to the total or partial use of the sample and aggregate data. Along the lines of [START_REF] Sun | A hierarchical mixture modeling framework for population synthesis[END_REF], we have classified the methods into three categories: synthetic reconstruction (SR), combinatorial optimization (CO), and statistical learning (SL). The SR approach combines information from the sample and the aggregate data and moreover computes weights that reflect the representativeness of each household in the the sample within a given zone. The CO methods also use the sample and aggregate data in order to select an appropriate combination of households that best fits the marginals. The third and last methods (SL) merely consider the sample and focus on the joint distribution of all attributes by estimating a probability for each combination.

The choice of appropriate method closely depends on the amount, type and quality (representativeness and comprehensiveness) of available data [START_REF] Rich | Large-scale spatial population synthesis for denmark[END_REF]. The overwhelming majority of statistical institutes make available to the public two kinds of data. The first is a disaggregated dataset, consisting of data for a sample of the population. Such a sample is typically referred to as a Public Use Micro Sample (PUMS); it is commonly compiled from census data and provides information on the socio-demographic characteristics of individuals or households (gender, profession, household size, household income, etc.) for a specific zone. The second source consists of aggregate data that provide the marginal distributions of socio-demographic variables covering a specific zone. These variables and distributions are referred to as the marginals or control variables [START_REF] Templ | Simulation of synthetic complex data: The r package simpop[END_REF], and their aspects differ from one country to the next. France differs from many countries in two regards: the sample made available is quite large (30% of the population as opposed to often less than 5% elsewhere); and all data stem from the same source (French census), which ensures data consistency.

Based on a review of the methods available to generate a synthetic population that jointly controls household and individual attributes [START_REF] Yaméogo | Comparing methods for generating a two-layered synthetic population[END_REF], it can be concluded that, given the characteristics of the French population data, SR methods are the most appropriate. We introduce herein four different algorithms from the SR family, namely Hierarchical Iterative Proportional Fitting (HIPF), Iterative Proportional Update (IPU), Generalized Raking (GR), and relative entropy minimization (ent), within a common framework so as to harmonize notations. We then test and compare the four algorithms by generating a two-layered population for each municipality within the Nantes Urban Area (western France). These methods however produce fractions of households and individuals, a problem that can be solved by converting the fractions into integers through an integerization process. To achieve this step, we apply two probabilistic integerization methods: the proportional probabilities approach and the truncation replication and sampling (TRS) method.

The objective of this paper therefore is to introduce and assess these various methods. To the best of our knowledge, no published research has quantitatively compared these specific approaches on the basis of a common conceptual framework (featuring a harmonization of notations, detailed description of each method, use of a case study, and application of quantitative performance metrics proposed in the scholarly literature).

The remainder of the paper is organized as follows. The second section reviews the existing population synthesis methods. The third then formally introduces the algorithms used for population generation within a common framework (in harmonizing notations). The fourth section is devoted to presenting the data and case study. The fifth section provides and discusses the results of our analyses followed by a conclusion offering perspectives on this paper.

Literature review

The methods utilized to generate a synthetic population can be grouped into three main categories: Synthetic Reconstruction (SR), Combinatorial Optimization (CO), and Statistical Learning (SL) [START_REF] Sun | A hierarchical mixture modeling framework for population synthesis[END_REF]. These methods will be described and compared hereinafter.

Synthetic Reconstruction

This category of methods is the most widely used to generate synthetic populations. A synthetic population is produced according to a two-step procedure: fitting and allocation. The fitting step involves assigning positive weights to the individuals and households contained in the sample with the resulting weights typically being non-integers. During the allocation step, these non-integer weights are converted into integer weights in order to replicate individuals and households.

SR methods are deterministic methods, meaning that depending on the sample studied, the weights obtained during the fitting step never vary. The prerequisite to applying SR methods is to possess both a sample and aggregate data. The underlying assumptions here are twofold: the sample represents the true correlation structure among the attributes [START_REF] Farooq | Simulation based population synthesis[END_REF]; and the interactions present in the sample are, to a great extent, preserved for the synthetic agents (Müller & Axhausen 2010). The sample therefore needs to be consistent, representative and composed of at least one observation for each type of individual in the actual population.

One of the commonly used SR techniques is Iterative Proportional Fitting (IPF) [START_REF] Beckman | Creating synthetic baseline populations[END_REF][START_REF] Pritchard | Advances in population synthesis: fitting many attributes per agent and fitting to household and person margins simultaneously[END_REF], which adjusts a contingency table constructed from the sample so as to match marginal distributions.

In its original formulation, IPF cannot simultaneously estimate both household and individuallevel attributes. Some IPF-based algorithms have attempted to address both household and individual attributes [START_REF] Arentze | Creating synthetic household populations: problems and approach[END_REF][START_REF] Guo | Population synthesis for microsimulating travel behavior[END_REF][START_REF] Auld | Efficient methodology for generating synthetic populations with multiple control levels[END_REF][START_REF] Zhu | Synthetic population generation at disaggregated spatial scales for land use and transportation microsimulation[END_REF][START_REF] Pritchard | Advances in population synthesis: fitting many attributes per agent and fitting to household and person margins simultaneously[END_REF]. However, in all these studies, the joint distribution of household and individual-level attributes is fitted either separately or sequentially which fails to guarantee the consistency between these two levels. Another approach consists of a fitting stage using IPF and a simulation stage where individuals are grouped into households with a household allocation procedure using the concept of "spouse matching" and "kids matching" [START_REF] Rich | Large-scale spatial population synthesis for denmark[END_REF].

In order to generate a two-layered synthetic population, four main algorithms have been proposed: Iterative Proportional Update (IPU) [START_REF] Ye | A methodology to match distributions of both household and person attributes in the generation of synthetic populations[END_REF], Hierarchical Iterative Proportional Fitting (HIPF) [START_REF] Müller | Multi-level fitting algorithms for population synthesis[END_REF][START_REF] Müller | A generalized approach to population synthesis[END_REF], relative entropy minimization (ent) [START_REF] Lee | Cross-entropy optimization model for population synthesis in activitybased microsimulation models[END_REF], and Generalized Raking (GR) [START_REF] Deville | Generalized raking procedures in survey sampling[END_REF]. In effect, these techniques generate populations of individuals grouped into households by computing household-level weights that satisfy the marginals at both the household and individual levels. Such algorithms will prove to be the most appropriate for the case study presented below, in considering the available input data, and will be presented in greater detail in the third section.

Combinatorial Optimization

The second category of approaches falls under to Combinatorial Optimization (CO) techniques. CO based techniques are two-layered since they can directly generate a list of households and individuals [START_REF] Ma | Synthetic population generation with multilevel controls: A fitnessbased synthesis approach and validations[END_REF].

Similar to SR methods, CO requires information on both the sample and marginal level, with the synthetic population being obtained by replicating individuals (without explicitly determining the joint distribution across all controlled attributes). But unlike SR methods, Combinatorial Optimization is a stochastic process. The data requirements for CO methods are less restrictive than those for SR methods [START_REF] Templ | Simulation of synthetic complex data: The r package simpop[END_REF], though on the other hand they do suffer from computational complexity when the population size is large [START_REF] Lee | Cross-entropy optimization model for population synthesis in activitybased microsimulation models[END_REF]. A description of this method has been given by [START_REF] Voas | An evaluation of the combinatorial optimisation approach to the creation of synthetic microdata[END_REF] and [START_REF] Templ | Simulation of synthetic complex data: The r package simpop[END_REF].

Statistical Learning

The third approach available to to generate a two-layered synthetic population is Statistical Learning (SL), also known as the simulation-based approach. SL focuses on the joint distribution of all attributes in the sample by directly estimating a probability for each combination, including those not observed in the sample [START_REF] Sun | A hierarchical mixture modeling framework for population synthesis[END_REF].

SL methods offer greater flexibility in terms of data requirements and data sources; in general, they display good performance both in treating the lack of heterogeneity problem encountered in SR and CO [START_REF] Sun | A hierarchical mixture modeling framework for population synthesis[END_REF]) and with small samples [START_REF] Borysov | How to generate micro-agents? a deep generative modeling approach to population synthesis[END_REF][START_REF] Sun | A hierarchical mixture modeling framework for population synthesis[END_REF]). However, a major drawback of SL methods is their inability to satisfy the conditional distributions while satisfying the marginal distributions of all variables simultaneously. During the population synthesis process, when marginals are available, it is indeed necessary to precisely match the observed marginal distributions with the population generated at the zonal level. Some of the two-layered SL-based algorithms derived for synthetic population generation include: the hierarchical Chain Monte Carlo method (hMCMC) [START_REF] Farooq | Simulation based population synthesis[END_REF], the Bayesian Networks-based method [START_REF] Sun | A bayesian network approach for population synthesis[END_REF][START_REF] Zhang | Connected population synthesis for transportation simulation[END_REF], hierarchical mixture modeling (HM) [START_REF] Sun | A hierarchical mixture modeling framework for population synthesis[END_REF], and deep generative modeling based on a Variational Autoencoder (VAE) [START_REF] Borysov | How to generate micro-agents? a deep generative modeling approach to population synthesis[END_REF]. [START_REF] Ye | On iterative proportional updating: Limitations and improvements for general population synthesis[END_REF] proposed a tensor decomposition method to guarantee the consistency between three levels of constraints: individual, household and enterprise.

A comparison of methods

Two-layered SR methods (IPU, HIPF, ent and GR) can generate high-quality two-layered synthetic populations that closely represent the actual population. Nonetheless, such techniques require a major preprocessing effort and are very stringent in terms of data needs. In fact, they require a representative sample and aggregate statistics at both the individual and household levels (Chapuis & Taillandier 2019). CO methods are less restrictive on data quality than SR methods for generating a two-layered synthetic population, yet they cannot always guarantee the optimal solution with respect to matching marginals and moreover require too much computing time. This CO category is better suited for generating small synthetic populations. SL methods are able to produce consistent results even for small sample sizes and generate a synthetic population from sample data only when necessary. On the other hand, SL methods make it impossible to satisfy marginal distributions of variables, which constitutes a major drawback when these marginals are available. In some configurations, combining SR and SL methods could be the most relevant option to have an accurate synthetic population satisfying marginal distributions. For example, combining a Variational Autoencoder model with IPF and quota-based random sampling [START_REF] Borysov | How to generate micro-agents? a deep generative modeling approach to population synthesis[END_REF] or Bayesian Networks with Generalized Raking techniques [START_REF] Sun | A bayesian network approach for population synthesis[END_REF].

The aim of this paper is to generate two-layered synthetic populations using French census data. The particularity of this dataset is the availability of a representative sample at the municipality level. The sample size is roughly 30% of the total municipal population; furthermore, aggregate statistics for both individual and household attributes are available which ensures data consistency. The data requirements for using SR methods in order to generate a synthetic population of individuals and households are therefore being met. Hence, SR methods are best suited since neither CO nor SL methods will not provide any advantage over SR. CO methods will in fact limit the population size potentially generated while SL methods prevent fitting to the marginals. The SR methods adopted to generate the synthetic population will be detailed in the next section.

Synthetic population generation methodology

The synthetic population is generated using a two-step procedure: 1-fitting, and 2-allocation (see Figure 1). In the first subsection, the four two-layered SR methods (Iterative Proportional Update (IPU), Hierarchical Iterative Proportional Fitting (HIPF), relative entropy minimization (ent) and Generalized Raking (GR)) available for use during the first step are presented. The second subsection then describes the two methods (proportional probabilities approach and TRS method) that convert non-integer weights resulting from the fitting step into integer weights in order to replicate individuals and households.

The fitting step

The objective of this step is to find the vector of household weights: W = (w h ), where h = 1 . . . n s h . n s h is the number of households in the sample and w h is a positive real measuring the importance of the corresponding household. This weight will be used in the allocation step to repeat or draw its corresponding household. Marginals are modeled as constraints on the weights vector.

We propose formulating this problem within the framework of the regularization of ill-posed inverse problems 1 in order to clarify the comparison among the various algorithms. From this point of view, the objective here is to find W that satisfies the marginal constraints, i.e.:

O H • W = M H O I • W = M I W ≥ 0 (1)
M H (resp. M I ) is the n mh × 1 (resp. n mi × 1) vector of marginals for the households (resp. individuals) (with n mh constraints on the households and n mi on the individuals). O H (resp. O I ) is the n mh × n s h (resp. n mi × n s h ) occurrence matrix that codes the sample according to the marginals on the households (resp. on the individuals). The next section on IPU will detail these equations.

This problem is ill-posed inasmuch as there are more variables than constraints (n m = n mh + n mi < n s h ). The constraints are often inconsistent with one another. The solution therefore must be regularized. An intuitive regularization consists of seeking a solution that is not too far removed from the sample, i.e. the vector solution Ŵ is not too distant from the vector of the prior weights, W prior , which models the sample. The following optimization problem serves to translate this idea.

Ŵ = arg min

W ≥0 O H •W M H O I •W M I dist(W, W prior ) (2)
with dist being a measurement of the distance between the vector W and the prior weight2 . Without consistent information on sampling, all the components of the prior weight vector have the same value: 1 (W prior = 1). The proposed methods tolerate some small deviations to the marginal constraints, which is why the constraints are no longer strict equalities. In the following, O is the concatenation of the occurrence matrices

O H and O I , O = O H t O I t t . Using the same notation, M = M H t M I t t .
The Statistical Reconstruction (SR) methods described in this paper, i.e. Iterative Proportional Update (IPU), Hierarchical Iterative Proportional Fitting (HIPF), Relative Entropy Minimization (ent) and Generalized Raking (GR), can all be interpreted within this common framework: these methods are in fact different views of the regularized Problem 2 of ill-posed Problem 1. For ent and GR, the minimization is explicit though the distance measurement differs. For IPU and HIPF, the minimization is implicit.

Iterative Proportional Update offers a geometric point of view of Problem 2. The IPU method starts from the sample, with initial weights being uniform. This vector is projected onto the hyperplane corresponding to household constraints before being projected onto a second hyperplane corresponding to the constraints on individuals. The process is iterative in support of the purpose of this algorithm to find a solution which is not too far removed from the initial sample and consistent with the constraints. It can be interpreted as a heuristic solution of Problem 2. IPU has been proven to have some limitations when generating a synthetic population at both individual and household levels. In particular, [START_REF] Ye | On iterative proportional updating: Limitations and improvements for general population synthesis[END_REF] have shown that theoretically, IPU is unable to converge to an optimal population distribution that simultaneously satisfies the constraints from individual and household levels. The authors have proposed an extension of IPU in order to address IPU failures. However, in our use case, IPU generates suitable solutions because the sample is large.

Hierarchical Iterative Proportional Fitting presents a dual view of Problem 2: it minimizes the distance to both constraint types on households and individuals, starting from uniform initial weights. The weights are modified as little as possible while optimizing the distance to the constraints.

Relative Entropy Minimization conveys a probabilistic point of view of Problem 2. The objective is to determine a probability, p h , associated with each household that can be interpreted as the weight, w h , divided by the number of households in the target population, n h . The solution must satisfy the marginal constraints and minimize the relative entropy to a prior, nearly uniform probability. By dividing in Problem 2, the weights vectors by the size of the target population and then by replacing the distance operator dist(p h , p prior h ) by p h log( p h p prior h ), the entropy formulation can be derived.

Generalized Raking provides an optimization point of view of Problem 2. It proposes solving this problem by setting up the Lagrangian.

After this more comprehensive introduction of the framework for treating inverse problems, the various methods will now be presented in greater depth.

Iterative Proportional Update

The Iterative Proportional Update (IPU), developed by [START_REF] Ye | A methodology to match distributions of both household and person attributes in the generation of synthetic populations[END_REF], is an iterative heuristic algorithm that simultaneously controls individual and household-level marginals during the fitting procedure. The corresponding mathematical optimization problem can be formulated with the following objective function [START_REF] Ye | A methodology to match distributions of both household and person attributes in the generation of synthetic populations[END_REF]:

min w h j h o j,h w h -m j /m j 2 (3)
Subject to w h 0, where: h denotes a household (h=1, 2,...,n s h ); j denotes the constraint or population characteristic of interest (j =1, 2,...,n mh ); and o j,h represents the frequency of the constraint j in household h (i.e. the occurrence), as one element of the matrix of occurrence, O. Moreover, w h is the weight attributed to the h th household and m j the value of constraint j.

The objective function measures the inconsistency between the weighted sample and the given constraints. At the first iteration, all households have a weight of one. IPU typically starts by adjusting weights to satisfy household constraints first, then updating them to satisfy individual constraints. At each iteration, a statistical measurement δ provides a goodness-of-fit result; it is the average of the absolute value of the relative difference between the weighted sum and the constraints, i.e.:

δ = j [|( h o j,h w h -m j )| /m j ] n m (4)
with n m = n mh + n mi being the number of marginals.

The gain in fit between two consecutive iterations is then calculated (∆ = |δ a -δ b |). The entire process is continued until the gain in fit is negligible or below a preset tolerance level. This tolerance level serves as the convergence criterion for terminating the algorithm [START_REF] Ye | A methodology to match distributions of both household and person attributes in the generation of synthetic populations[END_REF]).

Hierarchical Iterative Proportional Fitting

The HIPF algorithm [START_REF] Müller | Hierarchical ipf: Generating a synthetic population for switzerland[END_REF][START_REF] Müller | A generalized approach to population synthesis[END_REF] converts the householdlevel weights into individual-level weights and vice versa. It also proceeds in iterations and the procedure can be defined as follows [START_REF] Müller | Hierarchical ipf: Generating a synthetic population for switzerland[END_REF]:

k ← 0 w 0 h ← 1 for all h ∈ S repeat w (k+1) h ← FIT(w (k) h , m h a , m h b , ...) for all h ∈ S w (k+2) hi ← w (k+1) h for all i ∈ S for all i ∈ I(h) w (k+3) hi ← FIT(w (k+2) hi , m i α , m i β , ...) for all h ∈ S for all i ∈ I(h) w (k+4) h ← 1 nmi(h) i∈I(h) w (k+3) hi for all h ∈ S estimate w (k+5) h from w (k+4) h
by adjusting the individuals-per-household ratio using the relative entropy minimizing.

k ← k + 5 until convergence return w (k) h
where: h denotes a household, i an individual, k the iteration number, w

h the weight attributed to the h th household, w At the first iteration, all households have a weight of one. For all households, weights are computed to fit household-level control totals and converted to individual-level weights. These weights are then used as initial values to estimate new individual-level weights to fit the individuallevel control totals. The next step (Step 5) is to convert these new individual-level weights to household-level weights by considering that the weight of each household equals the average of the sum of the weights of the individuals in that household.

The sixth step consists of recomputing new household weights (w

(k+5) h
) by minimizing the relative entropy from weights obtained in Step 5 (w

(k+4) h
) to these latest weights, as defined below:

D w (k+5) h ||w (k+4) h = h w (k+5) h ln w (k+5) h w (k+4) h (5)
subject to the following constraints:

n s h h=1 w (k+5) h = n h (6) n s h h=1 nmi(h) i=1 w (k+5) hi = n (7)
where: n h represents households totals, n individuals totals, and n mi (h) the number of individuals in household h. This process is then repeated until convergence.

Entropy minimization

A mathematical formulation, using the relative entropy minimization function as the objective function, to generate synthetic data was proposed by Bar-Gera et al. (2009) and [START_REF] Lee | Cross-entropy optimization model for population synthesis in activitybased microsimulation models[END_REF]. According to this approach, both household and individual-level characteristics are contained in the constraints. The entropy optimization (ent) method described in this section closely follows that of [START_REF] Lee | Cross-entropy optimization model for population synthesis in activitybased microsimulation models[END_REF].

Let's consider the following notations: n h and n are respectively the total number of households and total population in the research area; n v and n u respectively the number of household-level and individual-level characteristics (factors); α and β are two subsets of respectively {1, 2, . . . , n v } and {1, 2, . . . , n u } (α and β will be used to model the marginals); x h v represents one household-level characteristic and x i u represents one individual-level characteristic; x h α represents the householdlevel characteristics associated with subset α (x h α = (x h v ) v∈α ), while x i β represents individual-level characteristics associated with β (x i β = (x i u ) u∈β ). Furthermore, let's consider: nx i u is the number of people in household h with person-level characteristic x i u where nx i u ∈ N and where u = 1, 2, . . . , n u ; nx i , is the vector of possible number of people in a household with a given individual-level characteristic and nx i equals {nx i 1 , . . . , nx i u , . . . , nx i nu }; h v denotes one possible value of x h v , where h v ∈ Ω v , and Ω v is a finite domain of values of x h v , where v is equal to 1, 2, . . . , n v ; h α denotes one possible value of x h α , where h α ∈ v∈α Ω v ; and h is the vector of all possible values of x h . nx i β is the number of people in household h with person-level characteristics x i β , where nx i β ∈ u∈β N.

Using the above notations, p α (h α ), E β (nx i β ), and p [h,nx i ] are defined as follows:

• p α (h α ) = joint distribution across household-level characteristics x α , where p α (h α ) = m h (h α )/n h and where m h (h α ) is the aggregate summary count across household-level characteristic x h α ;

• E β (nx i β ) = expected number of people in one household across person-level characteristics

x i β , where E β (nx i β )) = m i (x i β )/n h ; and m i (x i β ) is the count of person-level characteristics x i β , with β m i (x i β ) = n; • p [h,nx i ] =
min p [h,nx i ] D(p [h,nx i ] ||p prior [h,nx i ] ) = h,nx i p [h,nx i ] ln p [h,nx i ] p prior [h,nx i ] (8) subject to {h,nx i |v / ∈α} p [h,nx i ] = p α (h α ) ∀h v ∈ Ω v , v = 1, 2, ..., n v , h α ∈ v∈α Ω v (9) nx i β nx i β   {h,nx i |u / ∈β} p [h,nx i ]   = E β (nx i β ) ∀nx i u ∈ N, u = 1, 2, ..., n u , nx i β ∈ u∈β N (10) p [h,nx i ] ≥ 0 h,nx i p [h,nx i ] = 1 (11)
This formulation is an implementation of Problem 2, in considering probability p [h,nx i ] instead of weight w h , by inputting in Equation 9 and 10 the constraints on households and on individuals and by instantiating the distance measurement dist p h , p prior h with p h log p h p prior h .

Generalized Raking

The Generalized Raking (GR) techniques were developed by [START_REF] Deville | Generalized raking procedures in survey sampling[END_REF] to generate a synthetic population of both individuals and households. These techniques allow adjusting sampling weights in order to match known population totals. The problem formulation aligns with that of [START_REF] Deville | Generalized raking procedures in survey sampling[END_REF] and [START_REF] Müller | A generalized approach to population synthesis[END_REF].

Let's now consider a finite population U = {1, . . . , h, . . . , n h } with a response variable y h ∈ R. A sample s of size n s h is drawn from U with a given sampling design such that the inclusion probabilities p inclusion h,x = P (x ∈ s) are known. Let y h be the value of a variable of interest y , for the h th population element. The objective then is to estimate the finite population total t y = h∈U y h . An unbiased commonly used estimator of y is the Horvitz-Thompson estimator:

ŷ = h∈s 1 p inclusion h,x y h = h∈s w prior h y h ( 12
)
where w prior h are the prior weights (inverse of the prior probabilities). However, we have no information about y and only have n m auxiliary variables X i = (x i1 , . . . , x ij , . . . , x inm ) ∈ R for each member of the population. Also, the vector-valued population total M := h∈U X h is known accurately (i.e. the auxiliary variables and vector-valued population total correspond respectively to the occurrence matrix and marginals vector in Problem 2). In order to estimate y, we must seek new weights denoted w h , by modifying the prior weights w prior h in light of the auxiliary information while remaining close to the original weights. Let's consider a distance function G to minimize the gap between w h and w prior h subject to the constraints h∈s w h X h = h∈U X h = M . G must be positive and strictly convex, with G(1) = G (1) = 0 and G (1) = 1. In the context of synthetic population generation, we hold a sample to match the aggregate data, and the auxiliary variables are the marginals. The objective then is to minimize the difference existing between initial weights and final weights in order to fit the constraints at both the individual and household levels. The objective function is given by the following formula:

min w h h∈s w prior h G(w h /w prior h ) (13)
subject to:

h∈s w h X h = M (14)
This problem can now be solved by introducing a vector of Lagrange multipliers as demonstrated in [START_REF] Deville | Generalized raking procedures in survey sampling[END_REF].

This formulation is an implementation of Problem 2, by instantiating the distance measurement

dist w h , w prior h with w prior h G( w h w prior h
) and by identifying X h with the occurrence matrix O.

Generation process comparison across the four methods

Iterative Proportional Update (IPU), Hierarchical Iterative Proportional Fitting (HIPF) and relative entropy minimization (ent) all generate populations of individuals grouped into households by computing household-level weights that satisfy the marginals at both the household and individual levels. The HIPF algorithm constantly switches between household and individual domains, in employing an entropy-optimizing adjustment step [START_REF] Müller | Hierarchical ipf: Generating a synthetic population for switzerland[END_REF]. With the IPU and ent algorithms, weights are adjusted to satisfy household-level constraints first and then updated to satisfy individual-level constraints. The difference between IPU and ent lies in the procedure applied to adjust weights for a given individual-level control: if a household contains two or more individuals of the same category, ent reweights this household more heavily than a household with just one individual from this category, while IPU makes no distinction [START_REF] Müller | A generalized approach to population synthesis[END_REF]. The Generalized Raking method directly adjusts weights to satisfy both individual and household-level constraints.

The allocation step

All the methods described above generate fractional weights of households and individuals, making the results difficult to analyze. To construct the final population, we thus need to integerize these weights. The integerization process refers to converting these fractional weights into integer weights. To achieve this, two probabilistic methods are used: the proportional probabilities approach, and the truncate replicate sample (TRS) method. According to [START_REF] Lovelace | Evaluating the performance of iterative proportional fitting for spatial microsimulation: new tests for an established technique[END_REF], both of them outperform deterministic methods (simple rounding, threshold approach) in terms of final population counts and accuracy.

The proportional probabilities approach

The proportional probabilities (PP) approach considers fractional weights as probabilities [START_REF] Lovelace | Evaluating the performance of iterative proportional fitting for spatial microsimulation: new tests for an established technique[END_REF][START_REF] Joubert | Synthetic populations of south african urban areas[END_REF]. For example, the probability p h of a given household lies in the final synthetic population is thus given by: p h = w h / w h . The higher the fractional weight, the more likely an individual/household lies in the final population. As a result, an individual with a very high weight may be replicated several times, while one with a very low weight might not be included in the final synthetic population.

The TRS approach

The TRS approach (Lovelace & Ballas 2013) combines deterministic and probabilistic sampling in order to generate integer weights according to a three-step process: truncation, replication, and sampling.

1. The truncation step yields integer values by removing all information to the right of the decimal point. The decimal remainders (between 0 and 1) are then kept. As an illustration, a household with a weight of 4.65 will have a truncated value of 4. Its decimal remainder is 0.65.

2. During the second step, individuals/households are replicated depending on their integer weights obtained during the truncation step. Only truncated weights greater than 0 are replicated. For example, the household with a weight of 4.65 will be replicated 4 times.

Another household with a weight of 0.99 would not be replicated in this step (its truncated value is 0). When performing truncation and replication, no chance of oversampling exists (i.e. the sum of all integer weights will always be less than the population size).

3. During the last step, only the decimal weight remainders are included in applying a weighted random sampling without replacement. The rest of the individuals/households are selected from the entire sample, with selection probabilities set equal to the decimal weight remainders.

In our example, the household with the starting weight of 4.65 will have a 0.65 probability of being chosen again, while the other household will have a 0.99 probability.

This section has demonstrated how four SR methods function in abstract terms; a test scenario is now needed to conduct a practical comparison. The next section will describe the case study and data implemented.

Case study

The performance of the various methods described above will now be assessed using data drawn from the French census. This dataset has been collected by the French National Institute of Statistics and Economic Studies (INSEE). Since 2004, this census has covered all municipalities and is valid over a five-year period. By compiling successive five-year surveys, an array of population statistics could be obtained. To build a more robust database, the collected data were then adjusted to a single reference date, thus ensuring that all municipalities were being treated equally. This reference date was set on January 1 st of the median five-year survey period. 4The data provided by INSEE are available in two distinct forms: a sample of individuals and households, and control variables, both at the level of an IRIS (acronym for "aggregated units for statistical information"), which represents the basic unit for dissemination of intra-municipal data. Municipalities with over 10,000 inhabitants, and a large proportion of those with 5,000 to 10,000 population, are divided into several IRIS units and, by extension, all municipalities not divided into IRIS units constitute IRIS units in themselves.

We are specifically using census data from the Nantes Urban Area5 (NUA) from 2015 (these data were collected from 2013 to 2017). The total population of the NUA was approximately 949,000 individuals, residing in 418,000 households within 307 IRIS or equivalent units. The sample included 287,000 individuals and 136,000 households. Each observation in the sample represents a unique individual with his or her personal characteristics, as well as the household and main residence characteristics. Table 1 describes the attributes used in the generation process.

Descriptive statistics

Table 1 describes the variables collected for all 307 IRIS included in the sample. For our analysis, we considered 5 variables at the individual level and 4 variables at the household level. An IRIS contains on average of 1,363.3 households (±631.3) and 3,092.2 individuals (±1, 334.6); the samples contain on average 32.30% of households (±11.36) and 30.87% of individuals (±9.96) from the actual population.

Figure 2 and 3 display the distributions of the shares of the various categories of individual-level and household-level variables within the 307 IRIS. For most of these distributions, a fairly large variability can be observed.

Validation

The next section will compare the four previously described approaches to generating a synthetic population of households and individuals: IPU, HIPF, ent and Generalized Raking. 6 For each of these, we have used the proportional probabilities (PP) and truncate replicate sample (TRS) methods to integerize the weights. We have thus evaluated not only the performance of the four generation approaches but also that of the two integerization techniques. Two main aspects can be considered regarding an evaluation of the accuracy of a synthetic population: internal validation and external validation. Internal validation consists of comparing the variables of the synthetic population with the marginals in order to test the reliability of the generated data (e.g. does the estimated distribution of family composition correspond to distribution given by the census data?). In other words, an internal validation tests the ability of the population to fit with aggregate data. A validation is external if the estimated variables of the synthetic population are compared with another data source not used in the estimation process. Our case study does not feature a data source external to the French census at the IRIS level. Hence, we have solely focused on the internal validation.

According to the literature, internal validation can be carried out on either variables (marginals are compared with corresponding ones in the synthetic population), cells or the entire synthetic population. Many quantitative methods are available for internal validation [START_REF] Lovelace | Evaluating the performance of iterative proportional fitting for spatial microsimulation: new tests for an established technique[END_REF][START_REF] Timmins | Validation of spatial microsimulation models: a proposal to adopt the bland-altman method[END_REF]. The following performance metrics have been considered herein:

The coefficient of determination R 2 is the square of the Pearson correlation; it is a quantitative indicator that varies between 0 and 1 and moreover reveals how closely the simulated values fit the census data. An R 2 value of 1 denotes a perfect fit, while an R 2 value close to zero suggests no correspondence between constraints and simulated values [START_REF] Lovelace | Evaluating the performance of iterative proportional fitting for spatial microsimulation: new tests for an established technique[END_REF].

Total absolute error (TAE) and the standardized absolute error (SAE). TAE is the sum of the difference between simulated values and the marginals and SAE is TAE divided by the total population.

Standardized Root Mean Squared Error (SRMSE). This indicator focuses on error dispersion and is used to evaluate the goodness of fit between the estimated synthetic population and the marginals; it is the one of the most common indicators used [START_REF] Lee | Cross-entropy optimization model for population synthesis in activitybased microsimulation models[END_REF][START_REF] Lovelace | Evaluating the performance of iterative proportional fitting for spatial microsimulation: new tests for an established technique[END_REF][START_REF] Sun | A bayesian network approach for population synthesis[END_REF][START_REF] Saadi | Hidden markov model-based population synthesis[END_REF]. A zero value indicates a perfect match between census data and synthetic population, while a high SRMSE value suggests a poor fit.

The Bland-Altman method. Widely employed in healthcare studies to compare two measurements of the same variable, this graphical method can also be used to complement the other indicators [START_REF] Timmins | Validation of spatial microsimulation models: a proposal to adopt the bland-altman method[END_REF]. The Bland-Altman method consists of plotting of the difference between simulated and census counts versus the averages of the two counts.

Results and Discussion

This section presents the results of the internal validation procedure. The four SR algorithms have been implemented in the open-source MultiLevelIPF7 extension to the R statistical software package.8 A synthetic population has been generated for each IRIS of the NUA.

Internal validation with R 2 , TAE, SAE and SRMSE

The validation results show that all the proposed methods produce synthetic populations that are representative of the actual population, yet some methods prove to be more efficient. The first indicator, R 2 , revealed that even though all methods tested performed well, the TRS integerization method yielded better results than the proportional probabilities method. Moreover, the results of the Generalized Raking method results were less accurate compared than the other three generation methods. A more detailed description of the R 2 results follows: HIPF or IPU combined with TRS (HIPF+TRS or IPU+TRS) yields coefficients above 0.99 for all individual and household-level variables; HIPF or IPU combined with the proportional probabilities method (HIPF+PP or IPU+PP) and entropy minimization combined with either the TRS or proportional probabilities method (ent+TRS or ent+PP) yield coefficients greater than or equal to 0.98 for all individual and household-level variables;

GR combined with either TRS or proportional probabilities method (GR+TRS or GR+PP) yield coefficients greater than or equal to 0.91 for all individual and household-level variables.

The R 2 validation method merely provides an indication of fit and is influenced by outliers. A further analysis based on three other indicators (TAE, SAE and SRMSE), is therefore displayed in Table 2. These results confirm that all methods are globally efficient, but entropy minimization and HIPF do outperform the others. Based on Table 2, the method can be ranked in the following order from most to least accurate: entropy minimization, HIPF, IPU and GR for the individual level; HIPF, entropy mimimization, IPU and GR for the household level;

TRS and proportional probabilities.

According to all the validation indicators considered (R 2 , TAE, SAE and SRMSE), it can be concluded as regards the generation methods, slight differences exist between entropy minimization and HIPF. Moreover, these two methods outperform IPU and GR. For the integerization methods, the TRS approach outperforms the proportional probabilities approach. HIPF and entropy minimization combined with TRS therefore provide the best possible approximation of the actual population.

IRIS-level analysis

In addition to the global analysis given above, a local analysis, IRIS by IRIS, has been conducted in order to identify the zones with the highest errors (i.e. IRIS with the highest SAE values). For each method tested, whether at the individual or household level, three IRIS always stood out. Table 3 presents the values for the two best methods.

A qualitative analysis of the constraints from these three IRIS underscores their particular characteristics. IRIS 136 and 237 are activity zones with a small number of households and individuals. The population of IRIS 136 (resp. 237) is 359 (resp. 327) households and 964 (resp. 810) individuals. In these two IRIS, 67% (resp. 51%) of households have just one member; also, most of the inhabitants of these IRIS are men (73% (resp. 62%)) and belong to the 18-54 age group. IRIS 162 is a residential area with 1,163 households and 2,133 individuals. However, a significant portion of the territory is occupied by a psychiatric hospital. In this IRIS, 65% of the households are single-member and 65% of the individuals are between 15 and 64 years old. In conclusion, the simulation runs prove to be accurate for all IRIS except a few due to the particular population breakdown of these IRIS.

Bland-Altman approach

A Bland-Altman plot analysis of the data has been performed for comparing the census and simulated values of each IRIS for a given variable. This graphical method studies the mean difference and constructs limits of agreement [START_REF] Bland | Measuring agreement in method comparison studies[END_REF]. The X-axis corresponds to the mean of the two values, and the Y-axis is the difference between these two values. The limits of agreement are defined by ± 1.96 × the standard deviation of the mean difference. Analysis of the plot can help to identify some anomalies such as systematic overestimation or underestimation of census values by a synthetic reconstruction approach [START_REF] Kalra | Decoding the bland-altman plot: basic review[END_REF].

Our analysis has been applied to the 400 possible cases (50 categories of variables × 8 synthetic reconstruction approaches). The average of the differences (in both real and absolute terms) between simulated and census values by category for each of method has been computed. Figure 4 shows the mean and standard deviation of these mean values for each fitting method. Let's note that the average of the differences between simulated values and census values, expressed in real terms lies close to zero for the HIPF and entropy methods with a rather low standard deviation. The two measurements listed in Figure 4 would seem to confirm the HIPF and entropy methods outperform Generalized Raking and IPU. For purposes of illustration, Figure 5 plots Bland-Altman values for the five categories of the household-level variable "Family composition", generated with the HIPF method (in association with TRS allocation approach). The Y-axis shows the difference between the two populations (synthetic and actual), while the X-axis depicts the average of the two values. Depending on the IRIS, the simulated values are in some cases higher and in other cases lower than the census values. The average of the differences (middle line) is close to zero. 95% of the data points lie within the 'limits of agreement'(top and bottom lines) which indicates that there is agreement between census and simulated values. Depending on the IRIS, the simulated values are in some cases higher and in other cases lower than the census values, then it suggests there is no consistent bias.

Conclusion

This paper has provided a synopsis of the synthetic methods aimed at generating a population of individuals and households. We offered a detailed description of four synthetic reconstruction methods for the fitting step through use of a common framework. These methods are Hierarchical Iterative Proportional Fitting (HIPF), Iterative Proportional Update (IPU), Generalized Raking (GR), and relative entropy minimization (ent). Two integerization methods were also discussed, namely proportional probabilities and truncation replication and sampling (TRS). Next, an evaluation was performed of the most relevant method for generating a two-layered synthetic population. These methods were implemented using the R language. A case study involving the synthesis of agents (418,000 households, 949,000 individuals) from the Nantes Urban Area (western France) was considered, beginning with a sample of 136,000 households, including 287,000 individuals and 15,350 marginals. The synthetic population was generated with four household-level attributes and six individual-level attributes.

Results were evaluated using four indicators: R 2 , TAE, SAE, and SRMSE. The validation findings indicate that all methods considered yield good results, i.e. a two-layered synthetic population whose aggregate characteristics lie close to the census marginals. However, some methods output better results than others. For the fitting step, entropy minimization (ent) and Hierarchical Iterative Proportional Fitting (HIPF) prove to be the most efficient methods. For the allocation step, the truncation, replication and sampling (TRS) approach outperforms proportional probabilities.

We believe that the comparison of different statistical reconstruction (SR) algorithms performed in this paper with common notations and a common theoretical framework will facilitate a better dissemination of existing algorithms. This common framework will stimulate the development of new algorithms and position them with respect to existing methods.

The next step is to spatially allocate households or add to the demographic characteristics of households and individuals other socio-economic variables such as income by using other databases, e.g. fiscal database. The mathematical model used in this article inspires our current research to propose data fusion algorithms that enrich the synthetic population.

2. For each control variable, we must have a R dataframe. In our case study, we have 10 variables, so we must have 10 R dataframes.

Conditions: Ensure data consistency. Each row of the sample dataframe must represent an individual with his or her personal and households characteristics, a unique personnal ID number, a household ID and IRIS ID.

Each row of a control variable dataframe must represent a category with the number of people/households in this category and iris ID.

Step 3 : Fitting step Use of the four two-layered SR methods (Iterative Proportional Update (IPU), Hierarchical Iterative Proportional Fitting (HIPF), relative entropy minimization (ent) and Generalized Raking (GR)).

Conditions:

-Install MultiLevelIPF package under the following link: devtools::install_github("krlmlr/MultiLevelIPF") and then call library(MultiLevelIPF).

Outputs: households weights for each algorithm.

Step 4 : Allocation step Use of TRS method.

Outputs: synthetic population of households to merge with individuals by household ID.

Step 5 : Validation step Use of classical performance metrics (R 2 , TAE, SRMSE and Bland-Altman) to compare algorithms.

Figure 1 :

 1 Figure 1: Simplified flowchart of Synthetic Reconstruction methods to generate a two-layered synthetic population (households and individuals)

  hi the weight attributed to the i th individual in household h, m h a and m h b are household-level control totals, and m i α and m i β individual-level control totals. Moreover, S = {1 . . . n s h } with n s h being the number of households in the sample. P (h) = {1 . . . n mi (h)}, whereby n mi (h) is the number of individuals in household h.

Figure 2 :

 2 Figure 2: Distributions of the shares of the various categories of individual-level variables within the 307 IRIS

Figure 4 :

 4 Figure 4: Bland-Altman approach: mean of difference (in both real and absolute terms) between simulated and census values for each fitting method.

Figure 5 :

 5 Figure5: Bland-Altman approach: differences between simulated and census values for the five categories of the household-level variable "Family composition" generated with the HIPF fitting method (associated with the TRS allocation approach)

  multiway proportion of households in the research area with household-level characteristics h = {h 1 , . . . , h v , . . . , h nv } and number of individuals with person-level characteristicsx i , nx i = {nx i 1 , . . . , nx i u , . . . , nx i nu } and u = 1, 2, . . . , n u . • p prior [h,nx i ] = prior p [h,nx i ] ,easily computed from the disaggregated sample.The objective is to minimize the relative entropy between p [h,nx i ] and p prior [h,nx i ] (i.e. the estimation of p [h,nx i ] must be discriminated from p prior [h,nx i ] with a minimum difference). 3 This objective function can be written as follows:

Table 1 :

 1 Individual and Household-level attributes

	Level	Variable	Definition	Categories
			[number of categories]
	Household	Fam	Family composition [5]	Single member; The nuclear family is a couple without chil-
				dren; The nuclear family is a couple with children; The nuclear
				family is a single-parent family; Other composition
		ProfRP	Profession of the reference	Farmers, tradespeople; Executive; Intermediate occupations;
			person [7]	Clerical support workers; Lower-skilled technical occupations;
				Retiree; Unemployed
		Size	Household size [2]	One person; Two persons or more
		Cars	Number of cars [3]	No car; One; Two or more
	Individual	Age	Age [12]	0-2; 3-5; 6-10; 11-14; 15-17; 18-24; 25-29; 30-39; 40-54; 55-64;
				65-79; 80/+
		Sex	Gender [2]	Female; Male
		Relate	Relationship to the house-	Household reference person; Other household member
			hold reference person [2]
		Prof	Profession [7]	Farmers, tradespeople; Executive; Intermediate occupations;
				Clerical support workers; Lower-skilled technical occupations;
				Retiree; Unemployed
		Wstat	Work status [7]	In fixed-term employment; Permanent employment; Self-
				employed; Unpaid apprenticeships for those 15 or older; Un-
				employed; Under 15 years old; Other non-active persons
		Wtime	Working time [3]	Full-time worker; Part-time worker; Not applicable

Table 2 :

 2 Summary results of TAE, SAE and SRMSE between simulated and constrained data

			Individual-level			Household-level
	Method	TAE	SAE (%)	SRMSE	TAE	SAE (%)	SRMSE
	IPU+TRS	87,046	1.53	0.0024	17,082	1.02	0.0012
	IPU+PP	188,191	3.30	0.0032	56,529	3.37	0.0027
	HIPF+TRS	53,436	0.94	0.0013	14,134	0.84	0.0007
	HIPF+PP	176,612	3.10	0.0027	54,564	3.26	0.0025
	ent+TRS	50,412	0.88	0.0008	16,621	0.99	0.0009
	ent+PP	168,567	2.96	0.0024	55,830	3.33	0.0026
	GR+TRS	252,621	4.43	0.0090	71,778	4.28	0.0098
	GR+PP	337,630	5.93	0.0093	108,368	6.47	0.0128

Note: IPU: iterative proportional update; HIPF: hierarchical iterative proportional fitting; GR: generalized raking; ent: entropy minimization; TRS: truncation, replication, sampling; PP: proportional probabilities.

Table 3 :

 3 IRIS with the highest SAE values

			Individual-level	Household-level
	Method	Iris Id	SAE (%)	SAE (%)
		136	35.25	6.63
	HIPF+TRS	237	13.71	2.82
		162	10.08	3.09
		136	14.22	36.20
	ent+TRS	237	3.71	10.37
		162	2.40	7.51
	Note: HIPF: hierarchical iterative proportional fitting; ent: entropy minimiza-
	tion; TRS: truncation, replication, sampling.

An abundance of literature exists of this subject ever since the seminal work by Tikhonov[START_REF] Tikhonov | Solutions of ill posed problems[END_REF].

In this paper, distance is not intended in its strict mathematical definition.

In the literature on synthetic population generation, this method is often called cross-entropy minimization; from a strictly mathematical point of view, it is not valid. We have chosen to replace the term cross-entropy by relative entropy, also known as Kullback-Leibler divergence, which is correct and consistent with the notation of D for the measurement of this relative entropy.

https://www.insee.fr/fr/information/2383265, Consulted on 22 April 2020

According to the INSEE Institute, an urban area is a group of adjoining municipalities, without pockets of clear land, encompassing an urban centre (urban unit) providing at least 10 000 jobs, and whose neighboring rural districts or suburban units (urban periphery) account for at least 40% of the employed residents working in the center or in the municipalities attracted by this center.

For the Generalized Raking approach, four distance functions G can be used: the linear method, the raking ratio method, the logit method, and the truncated linear method. We tested all these functions, but convergence was only achieved for the logit method.

https://github.com/krlmlr/MultiLevelIPF, Consulted on 24 April 2020

We used a computer of 2 x 2.60GHz CPU cores and 16 GB RAM.

Appendix A: Model Documentation

In this appendix, we describe in five steps the approach used in the paper in a more detailed fashion with a R script. The first step presents the databases used and how to download them. The other steps describe the statistical analyses performed with a toy model. Interested readers can directly contact the authors to get the complete codes.

Step 1: Databases access The databases used are available under the following links (consulted on 6 November 2020) :

1 Step 2: Data processing Inputs : downloaded databases Outputs :

1. A R dataframe for the sample.