Gregory Guernec
email: gregory.guernec@inserm.fr

Valerie Gares

Jeremy Omer

Philippe Saint-Pierre

Nicolas Savy

OTrecod: An R Package for Data Fusion using Optimal Transportation Theory

The advances of information technologies often confront users with a large amount of data which is essential to integrate easily. In this context, creating a single database from multiple separate data sources can appear as an attractive but complex issue when same information of interest is stored in at least two distinct encodings. In this situation, merging the data sources consists in finding a common recoding scale to fill the incomplete information in a synthetic database. The OTrecod package provides R-users two functions dedicated to solve this recoding problem using optimal transportation theory. Specific arguments of these functions enrich the algorithms by relaxing distributional constraints or adding a regularization term to make the data fusion more flexible. The OTrecod package also provides a set of support functions dedicated to the harmonization of separate data sources, the handling of incomplete information and the selection of matching variables. This paper gives all the keys to quickly understand and master the original algorithms implemented in the OTrecod package, assisting step by step the user in its data fusion project.

Introduction

The large amount of data produced by information technology requires flexible tools to facilitate its handling. Among them, the field of data fusion [START_REF] Hall | An introduction to multisensor data fusion[END_REF][START_REF] Klein | Sensor and data fusion: a tool for information assessment and decision making[END_REF][START_REF] Castanedo | A review of data fusion techniques[END_REF] also known as statistical matching [START_REF] Adamek | Fusion: combining data from separate sources[END_REF][START_REF] D'orazio | Statistical matching: Theory and practice[END_REF][START_REF] Vantaggi | Statistical matching of multiple sources: A look through coherence[END_REF] aims to integrate the overall information from multiple data sources for a better understanding of the phenomena that interact in the population. Assuming that two heterogeneous databases A and B share a set of common variables X while an information of interest is encoded in two distinct scales respectively: Y in A and Z in B. If Y and Z are never jointly observed, a basic data fusion objective consists in the recoding of Y in the same scale of Z (or conversely), to allow the fusion between the databases as illustrated in Table 1.

Providing a solution to this recoding problem is often very attractive because it aims at giving access to more accurate and consistent information with no additional costs in a unique and bigger database. Despite this, if we exclude all R data integration packages applied in the context of genetic area like the MultitDataset package [START_REF] Hernandez-Ferrer | Multidataset: an r package for encapsulating multiple data sets with application to omic data integration[END_REF], the OMICsPCA package [START_REF] Das | OMICsPCA: An R package for quantitative integration and analysis of multiple omics assays from heterogeneous samples[END_REF], or the MixOmics package [START_REF] Rohart | mixomics: An r package for 'omics feature selection and multiple data integration[END_REF] which are often only effective for quantitative data integration. To our knowledge, the StatMatch package (D' Orazio and D'Orazio, 2019) is actually the only one that provide a concrete solution to the problem using hot deck imputation procedures. The main reason for this relative deficiency is that this problem is, in fact, often assimilated and solved like a missing data imputation problem. According to this idea, a very large amount of works and reference books now exist about the handling of missing data [START_REF] Zhu | High-dimensional principal component analysis with heterogeneous missingness[END_REF][START_REF] Little | Statistical analysis with missing data[END_REF]. Moreover, we can enumerate several R packages that we can sort by types of imputation methods [START_REF] Mayer | R-miss-tastic: a unified platform for missing values methods and workflows[END_REF]: mice (Van Buuren and Groothuis-Oudshoorn, 2011) and missforest [START_REF] Stekhoven | Missforest-non-parametric missing value imputation for mixedtype data[END_REF] which use conditional models, softImpute [START_REF] Hastie | Matrix completion and low-rank svd via fast alternating least squares[END_REF] and missMDA [START_REF] Josse | missmda: A package for handling missing values in multivariate data analysis[END_REF] which apply low-rank based models. For all these packages, imputation performances can sometimes fluctuate a lot according to the structure and the proportion of non-response encountered. Regressions and non parametric imputation approaches (like hot-deck methods from donor based family) seem to use partially, or not at all, the available information of the two databases to provide the individual predictions. Contrary to our approach, all these methods only use the set of shared variables X for the prediction of Z in A (or conversely Y in B) without really taking into account Y and its interrelations with X in their process of predictions. The purpose of this paper is to present a new package to the R community (R Core [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF] called OTrecod which provides a simple and intuitive access to two original algorithms [START_REF] Gares | On the use of optimal transportation theory to recode variables and application to database merging[END_REF][START_REF] Gares | Regularized optimal transport of covariates and outcomes in data recoding[END_REF] dedicated to solve these recoding problems in the data fusion context by considering them as applications of Optimal Transportation theory (OT). In fact, this theory was already applied in many areas: for example, the transport package [START_REF] Schuhmacher | transport: Computation of Optimal Transport Plans and Wasserstein Distances[END_REF] solves optimal transport problems in the field of image processing. A specific package called POT: Python Optimal Transport also exists in the Python software (https:/pythonot.github.io) to solve optimization problems using optimal transportation theory in the fields of signal theory, image processing and domain adaptation. Nevertheless, all these available tools were still not really adapted to our recoding problem and the performances established by OT-based approaches to predict missing information compared to more standard processes [START_REF] Gares | On the use of optimal transportation theory to recode variables and application to database merging[END_REF][START_REF] Gares | Regularized optimal transport of covariates and outcomes in data recoding[END_REF][START_REF] Muzellec | Missing data imputation using optimal transport[END_REF] have finished to confirm our decision.

In OTrecod the first provided algorithm, called outcome and integrated in the OT_outcome function consists in finding a map that pushes the distribution of Y forward to the distribution of Z [START_REF] Gares | On the use of optimal transportation theory to recode variables and application to database merging[END_REF] while the second one, called joint and integrated in the OT_joint function, pushes the distribution of (Y, X) forward to the distribution of (Z, X). Consequently, by building, these two algorithms take advantage of all the potential relationships between Y , Z, and X for the prediction of the incomplete information of Y and/or Z in A and B. Enrichments related to these algorithms and described in [START_REF] Gares | On the use of optimal transportation theory to recode variables and application to database merging[END_REF]Cuturi, 2013) are also available via the optional arguments of these two functions. In its current version, these algorithms are accompanied by original preparation (merge_dbs, select_pred) and validation (verif_OT) functions which are key steps of any standard statistical matching project and can be used independently of the outcome and joint algorithms.

2 Solving recoding problems using optimal transportation theory

The optimal transportation problem

The optimal transportation (OT) problem was originally stated by [START_REF] Monge | Mémoire sur la Théorie des Déblais et des Remblais[END_REF] and consists in finding the cheapest way to transport a pile of sand to fill a hole. Formally the problem writes as follows: Consider two (Radon) spaces X and Y, µ X a probability measure on X, and µ Y a probability measure on Y and c a Borel-measurable function from X×Y to [0, ∞]. The Kantorovich's formulation of the optimal transportation problem [START_REF] Kantorovich | On the transfer of masses[END_REF] consists in finding a measure γ ∈ Γ(µ X , µ Y ) that realizes the infimum:

inf X×Y c(x, y) dγ(x, y) γ ∈ Γ(µ X , µ Y ) , (1)
where Γ(µ X , µ Y ) is the set of measures on X × Y with marginals µ X on X and µ Y on Y. This theory is applied here to solve a recoding problem of missing distributions in a data fusion area. To do so, we make use of Kantovorich's formulation adapted to the discrete case, known as Hitchcock's problem [START_REF] Hitchcock | The distribution of a product from several sources to numerous localities[END_REF]. Therefore, by construction, the proposed algorithms are usable for specific target variables only: categorical variables, ordinal or nominal, and discrete variables with finite number of values.

Optimal transportation of outcomes applied to data recoding

Let A and B be two databases corresponding to two independent sets of subjects. We assume without loss of generality that the two databases have equal sizes, so that they can be written as A = {i1, . . . , in} and B = {j1, . . . , jn}.

Let (Xi, Yi, Zi) i∈A and (Xj, Yj, Zj) j∈B be two sequences of i.i.d. discrete random variables with values in X × Y × Z, where X is a finite subset of R P , and Y and Z are finite subsets of R.

Variables (Xi, Yi, Zi), i ∈ A, are i.i.d copies of (X A , Y A , Z A ) and (Xj, Yj, Zj), j ∈ B, are i.i.d copies of (X B , Y B , Z B ). Moreover assume that (Xi, Yi, Zi), i ∈ A are independent of (Xj, Yj, Zj), j ∈ B . The first version using the optimal transportation algorithm approach, described in [START_REF] Gares | On the use of optimal transportation theory to recode variables and application to database merging[END_REF], assumes that:

Assumption 1. Y A and Z A respectively follow the same distribution as Y B and Z B .

Assumption 2. For all x ∈ X the probability distributions of Y A and Z A given that X A = x are respectively equal to those of Y B and Z B given that X B = x.

In this setting, the aim is to solve the recoding problem given by equation ( 1) that pushes

µ Y A forward to µ Z A .
The variable γ of ( 1) is a discrete measure with marginals µ Y A and µ Z A , represented by a |Y|×|Z| matrix. The cost function denoted as c is a |Y|×|Z| matrix, (cy,z)y∈Y,z∈Z . The goal is in the identification of:

γ * ∈ argmin γ∈R |Y|×|Z| + ⟨γ|c⟩ : γ1 |Z| = µ Y A , γ T 1 |Y| = µ Z A , (2)
where ⟨•|•⟩ is the dot product, 1 is a vector of ones with appropriate dimension and M T is the transpose of matrix M . The cost function considered by [START_REF] Gares | On the use of optimal transportation theory to recode variables and application to database merging[END_REF], cy,z, measures the average distance between the profiles of shared variables of A satisfying Y = y and subjects of B satisfying Z = z, that is:

cy,z = E d(X A , X B ) | Y A = y, Z B = z , (3)
where d is a given distance function to choose on X × X .

In fact, the above situation cannot be solved in reality, since the distributions of X A , X B , Y A and Z A are never jointly observed. As a consequence, the following unbiased empirical estimators are used: μX A n of µ X A and μX B n of µ X B . Because Y and Z are only available in A and B respectively, two distinct empirical estimators have to be defined:

μY A n,y = 1 n i∈A 1 {Y i =y} , ∀y ∈ Y μZ A n,z = 1 n j∈B 1 {Z j =z} , ∀z ∈ Z, (4)
where

1 {Y =y} = 1 if Y = y and 0 otherwise.
The assumption 1 gives:

µ Z A = µ Z B from which we can conclude that μZ B n,z = μZ A n,z . Finally, denoting: κn,y,z ≡ i∈A j∈B
1 {Yi=y,Zj=z} the number of pairs (i, j) ∈ A × B such that Yi = y and Zj = z, the cost matrix c is estimated by: ĉn,y,z =

1 κn,y,z i∈A j∈B 1 {Yi=y,Zj=z} × d(Xi, Xj), ∀y ∈ Y, z ∈ Z : κn,y,z ̸ = 0, 0, ∀y ∈ Y, z ∈ Z : κn,y,z = 0. (5)
Plugging the values observed for these estimators in (2) yields to a linear programming model denoted:

P0 n :                    min < ĉn, γ > s.t. z∈Z γy,z = µ Y A n,y , ∀y ∈ Y y∈Y γy,z = µ Z A n,z , ∀z ∈ Z γy,z ≥ 0, ∀y ∈ Y, ∀z ∈ Z (6)
The solution γn can then be interpreted as an estimator μ(Y A ,Z A ) n of the joint distribution of Y A and Z A , µ (Y A ,Z A ) . If this estimate is necessary, it is nevertheless insufficient here to provide the individual predictions on Z in A. These predictions are done in a second step using a nearest neighbor algorithm from which we deduce an estimation of µ Z A |X A =x,Y A =y (see [START_REF] Gares | Regularized optimal transport of covariates and outcomes in data recoding[END_REF] for details). In the remainder, the overall algorithm described in this section is referred to as outcome. To improve the few drawbacks of this algorithm described in [START_REF] Gares | On the use of optimal transportation theory to recode variables and application to database merging[END_REF], derived algorithms from outcome have been developed [START_REF] Gares | Regularized optimal transport of covariates and outcomes in data recoding[END_REF] and described in the following part.

Optimal transportation of outcomes and covariates

Using the same notations, [START_REF] Gares | On the use of optimal transportation theory to recode variables and application to database merging[END_REF] propose to search for an optimal transportation map between the two joint distributions of (X A , Y A ) and (X A , Z A ) with marginals µ (X A ,Y A ) and µ (X A ,Z A ) respectively. Under Kantorovich's formulation in a discrete setting, they search for:

γ * ∈ argmin γ∈D < c, γ >,
where c is a given cost matrix and D is the set of joint distributions with marginals µ (X A ,Y A ) and µ (X A ,Z A ) . It is natural to see any element γ ∈ D as the vector of joint probabilities P((X A = x, Y A = y), (X A = x ′ , Z A = z)) for any x, x ′ ∈ X 2 , y ∈ Y and z ∈ Z. Since this probability nullifies for all x ̸ = x ′ , γ ∈ D is defined as a vector of R |X |×|Y|×|Z| , where γx,y,z stands for an estimation of the joint probability P(X A = x, Y A = y, Z A = z). These notations lead to the more detailed model:

P :                    min < c, γ > s.t. z∈Z γx,y,z = µ (X A ,Y A ) x,y , ∀x ∈ X , ∀y ∈ Y y∈Y γx,y,z = µ (X A ,Z A ) x,z , ∀x ∈ X , ∀z ∈ Z γx,y,z ≥ 0, ∀x ∈ X , ∀y ∈ Y, ∀z ∈ Z (7)
The above algorithm can be solved only if the marginals µ (X A ,Y A ) and µ (X A ,Z A ) are known, but, based on assumption 2, unbiased estimators μX A ,Y A n and μX A ,Z A n can be built according to the previous subsection. For the first one it gives:

μX A ,Y A n = 1 n i∈A 1 {Y i =y,X i =x} , ∀x ∈ X , ∀y ∈ Y, (8)
The cost matrix introduced in the outcome algorithm is used (3) and estimated by ( 5). Formally we can write: cx,y,z = cy,z, ∀x ∈ X , ∀y ∈ Y, ∀z ∈ Z.

(9)

which does not depend on the value of x. Plugging the values observed for these estimators in (7) yield a linear programming model denoted as Pn. In contrast to outcome, the algorithm that consists in solving Pn to solve the recoding problem is referred to as joint in what follows.

An estimation of the distribution of Z A given the values of X A and Y A is then given by:

μZ A |X A =x,Y A =y n,z =      γn,x,y,z μ(X A ,Y A ) n,x,y , ∀x ∈ X , y ∈ Y, z ∈ Z : μ(X A ,Y A ) n,x,y ̸ = 0, 0, ∀x ∈ X , y ∈ Y, z ∈ Z : μ(X A ,Y A ) n,x,y = 0. (10)
and an individual prediction of Z A is then deduced using the maximum a posterior rule:

z A i = argmax z∈Z μZ A |X A =x i ,Y A =y i n,z
Due to potential errors in the estimations of P, the constraints of Pn may derive from the true values of the marginals of µ (X A ,Y A ,Z A ) . To deal with this situation, small violations of the constraints of Pn are allowed by enriching the initial algorithm as described in [START_REF] Gares | Regularized optimal transport of covariates and outcomes in data recoding[END_REF].

The equality constraints of Pn are then relaxed as follows:

z∈Z γx,y,z = μ(X A ,Y A ) n,x,y + e X,Y x,y , ∀x ∈ X , ∀y ∈ Y (11) y∈Y γx,y,z = μ(X A ,Z A ) n,x,z + e X,Z x,z , ∀x ∈ X , ∀z ∈ Z (12) x∈X ,y∈Y e X,Y x,y = 0, x∈X ,z∈Z e X,Z x,z = 0 (13) -e X,Y,+ x,y ≤ e X,Y x,y ≤ e X,Y,+ x,y , ∀x ∈ X , ∀y ∈ Y (14) -e X,Z,+ x,z ≤ e X,Z x,z ≤ e X,Z,+ x,z , ∀x ∈ X , ∀z ∈ Z (15) x∈X ,y∈Y e X,Y,+ x,y ≤ αn, x∈X ,z∈Z e X,Z,+ x,z ≤ αn ( 16
)
This relaxation is possible by introducing extra-variables e X,Y,+ and e X,Z,+ as additional constraints ( 14)-( 15). [START_REF] Gares | Regularized optimal transport of covariates and outcomes in data recoding[END_REF] suggests to consider αn := α √ n from ( 16), with a parameter α to calibrate numerically but proposes also a default value fixed to 0.4.

A regularization term λ given by ( πx,y,z

μX A n,x
)x∈X,y∈Y,z∈Z can also be added to improve regularity in the variations of the conditional distribution µ Y A ,Z A |X A =x with respect to x. The corresponding regularized algorithm is:

PR n :                              min < ĉn, γ > +λ (x i ,x j )∈E X wi,j y∈Y,z∈Z r + i,j,y,z s.t. constraints (11)-(16) γx i ,y,z μX A n,x i - γx j ,y,z μX A n,x j ≤ r + i,j,y,z , ∀{xi, xj} ∈ EX , y ∈ Y, z ∈ Z γx i ,y,z μX A n,x i - γx j ,y,z μX A n,x j ≥ -r + i,j,y,z , ∀{xi, xj} ∈ EX , y ∈ Y, z ∈ Z γx,y,z ≥ 0, ∀x ∈ X , ∀y ∈ Y, ∀z ∈ Z (17)
The constant λ ∈ R + is a regularization parameter to be calibrated numerically (0.1 can be considered as default value) and EX ⊂ X 2 includes the pairs of elements of X defined as neighbors: {xi, xj} ∈ EX if xj is among the k nearest neighbors of xi for some parameter k ≥ 1. The method that computes a solution to the recoding problem with regularization and relaxation is called rjoint.

In a same way, a relaxation of the assumption 1 is also proposed and added to the outcome algorithm: This resulting method is denoted r-outcome and is related to the following program: Note that algorithms joint and r-joint do not require assumption 1. The relaxation in routcome alleviates its dependence to the satisfaction of this assumption. However, algorithms joint and r-joint require that X is a set of discrete variables (factors ordered or not are obviously allowed) while the absence of X in the linear algorithms outcome and r-outcome allow X to be a set of discrete and/or continuous variables. In this case, the nature of the variables X need to be considered when choosing the distance d.

P0-R n :                                                    min < ĉn, γ > z∈Z γy,z = μY A n,y + e Y y , ∀y ∈ Y y∈Y γy,z = μZ A n,z + e Z z ,
3 Package installation and description

Installation

The OTrecod package can be installed from the Comprehensive R Archive Network (CRAN) by running the install.packages("OTrecod") command. The development version of OTrecod (https: //github.com/otrecoding/OTrecod) is also available and can be directly installed from GitHub by loading the devtools package and following the command: R> devtools::install_github("otrecoding/OTrecod")

Main functions

The two types of optimal transportation (OT) algorithms previously introduced (outcome and joint) and their respective enrichments (r-outcome and r-joint) are available in the OTrecod package via two core functions denoted OT_outcome and OT_joint. Details about their implementations in R are described in the following section. In this package, these algorithms of recoding are seen as fundamental steps of data fusion projects that also require often adapted preparation and validation functions. In this way, the table 2 introduces the main functions proposed by OTrecod to handle a data fusion project.

R Function Description Pre-process functions merge_dbs

Harmonization of the data sources select_pred Selection of matching variables Functions of data fusion OT_outcome Data fusion with OT theory using the outcome or r-outcome algorithms.

OT_joint

Data fusion with OT theory using the joint or rjoint algorithms.

Post-process function verif_OT

Quality assessment of the data fusion avg_dist_closest, indiv_grp_closest, indiv_grp_optimal), and their related documentations, are all included and usable separately in the package. They have been kept available for users to ensure a great flexibility as other interesting functions like power_set that returns the power set of a set. This function did not exist on R until now and could be of interest for specialists of algebra. These functions are not described here but detailed in the related pdf manual of the package.

Functionalities overview

Expected structure of the input databases as arguments

The functions of recoding OT_outcome and OT_joint require a specific structure of data.frame as input arguments. Described in Table 3, it must be the result of two overlayed databases made of at least four variables:

• A first variable, discrete or categorical, corresponding to the database identifier, stored in factor or not, but with only two classes or levels (for example: A and B, 1 and 2 or otherwise).

• The target variable of the first database (or top database) denoted Y for example, whose values related to the second database are missing. This variable can be discrete or categorical stored in factor, ordered factor or not.

• In the same way, the target variable of the second database (or below database) denoted Z for example, whose values related to the first database are missing.

• At least one shared variable (defined as a variable with the same label and the same encoding in the two distinct data sources). The type of shared variables can be continuous, categorical stored in factor or not, complete or not. Nevertheless, few constraints must be noticed related to this question. First, in a critical situation where only one shared variable exists, this latter cannot be incomplete. Second, continuous shared variables are actually not allowed in the current version of the function OT_joint, therefore, these variables must be transformed beforehand.

DB

Y Z X 1 X 2 X 3 1 (600-800] <NA> M Yes 50 1 (600-800] <NA> M No 32 1 [200-600] <NA> W No 31 2 <NA> G1 M No 47 2 <NA> G3 W Yes 43 2 <NA> G2 W No 23 2 <NA> G4 M Yes 22 2 <NA> G2 W Yes 47
Table 3: Example of expected structure for two databases 1 and 2 with three shared variables X 1 , X 2 , X 3

As additional examples, users can also refer to the databases simu_data and tab_test provided in the package with expected structures. Note that class objects are not expected here as input arguments of these functions to allow users to freely work with or without the use of the pre-process functions provided in the package.

Choice of solver

The package OTrecod uses the ROI optimization infrastructure [START_REF] Theußl | Roi: The r optimization infrastructure package[END_REF] to solve the optimization problems related to the outcome and joint algorithms. The solver GLPK (The GNU Linear Programming Kit Makhorin ( 2011)) is the default solver actually integrated in the OT_outcome and OT_joint functions for handling linear problems with linear constraints. The ROI infrastructure makes easy for users to switch solvers for comparisons. In many situations, some of them can noticeably reduce the running time of the functions.

For example, the solver Clp [START_REF] Forrest | Clp user guide[END_REF] for COINT-OR Linear Programming, known to be particularly convenient in linear and quadratic situations, can be easily installed by users via the related plug-in available in ROI (searchable with the instruction ROI_available_solvers()) and following the few instructions detailed in [START_REF] Theußl | Roi: An extensible r optimization infrastructure[END_REF] or via the dedicated website.

An illustrative example

In California (United States), from 1999 to 2018, the Public Schools Accountability Act (PSAA) imposed on its California Department of Education (CDE) to provide annually the results of an Academic Performance Index (API) which established a ranking of the best public schools of the state.

This numeric score, indicator of school's performance levels, could vary from 200 to 1000 and the performance objective to reach for each school was 800. Information related to the 418 schools (identified by cds) of Nevada (County 29) and to the 362 schools of San Benito (County 35), was respectively collected in two databases, api29 and api35, available in the package. The distributions of all the variables in the two databases are provided by following the R commands:

R> library(OTrecod) R> data(api29); data(api35) R> summary(api29) #- ------------------------------------------------ ------------------------------------------------ The two databases seem to share a majority of variables (same labels, same encodings) of different types, therefore inferential statistics could be ideally considered by combining all the information of the two databases to study the effects of social factors on the results of the API score in 2000.

Nevertheless, while this target variable called apicl_2000 is correctly stored in the database api29 and encoded as a three levels ordered factors clearly defined: [200-600], (600-800] and (800-1000]1 , the only information related to the API score in the database api35 is the variable apicl_1999 for the API score collected in 1999, encoded in four unknown balanced classes (G1, G2, G3 and G4). As no school is common to the two counties, we easily deduce that these two variables have never been jointly observed.

By choosing these two variables as outcomes (called also target variables), the objective of the following examples consists in creating a synthetic database where the missing information related to the API score in 2000 is fully completed in api35 by illustrating the use of the main functions of the package.

Harmonization of two datasources using merge_dbs

The function merge_dbs is an optional pre-process function dedicated to data fusion projects that merges two raw databases by detecting possible discrepancies among the respective sets of variables from one database to another. The current discrepancy situations detected by the function follow specific rules described below:

• any variable (other than a target one) whose label (or name) is not common to the two data sources is automatically excluded. By default, the remaining variables are denoted shared variables.

• among the subset of shared variables, any variable stored in a different format (or type) from one data source to another will be automatically discarded from the subset previously generated and its label saved in an output object called REMOVE1.

• among the remaining subset of shared variables, a factor variable (ordered or not) stored with different levels or number of levels from one data source to another will be automatically discarded from the subset and its label saved in an output object called REMOVE2.

These situations sometimes require reconciliation actions which are not supported by the actual version of the merge_dbs function. Therefore, when reconciliation actions are required by the user, they will have to be treated a posteriori and outside the function.

Applied to our example and according to the introduced rules, the first step of our data fusion project consists in studying the coherence between the variables of the databases api29 and api35 via the following R code: As entry, the raw databases must be declared separately in the DB1 and DB2 arguments and the name of the related target variables of each database must be specified via the NAME_Y and NAME_Z for DB1 and DB2 respectively. In presence of row identifiers, the respective column indexes of each database must be set in the argument row_ID1 and row_ID2. The arguments ordinal_DB1 and ordinal_DB2 list the related column indexes of all the variables defined as ordinal in the two databases (including also the indexes of the target variables if necessary). Here, apicl_2000 is clearly an ordinal variable, and, by default, we suppose that the unknown encoding related to apicl_1999 is also ordinal: the corresponding indexes (2 and 3) are so added in these two arguments.

After running, the function informs users that no row was dropped from the databases during the merging because each target variable is fully completed in the two databases. Nine potential predictors are kept while only one variable is discarded because of discrepancies between the databases: its identity is consequently stored in output and informs user about the nature of the problem: the mobility factor has different levels from one database to the other.

R> summary(step1) R> step1$REMOVE1

# List of removed variables because of type's problem NULL R> step1$REMOVE2 # List of removed factors because of levels' problem "mobility" R> levels(api29$mobility) # Verification "[0-20]" "(20-100]" R> levels(api29$mobility); levels(api35$mobility) "[0-20]" "(20-100]" # levels in api29 "1" "2" # levels in api35 R> step1$REMAINING_VAR "acs.core" "acs.k3.20" "api.stu" "awards" "ell" "full" "grad.sch" "meals" "stype"

The function returns a list object and notably DB_READY, a data.frame whose structure corresponds to the expected structure introduced in the previous subsection: a unique database, here result of the superimposition of api29 on api35 where the first column (DB) corresponds to the database identifier (1 for api29 and 2 for api35), the second and third columns (Y, Z respectively) corresponds to the target variables of the two databases. Missing values are automatically assigned by the function to the unknown part of each target variable: in Y when the identifier equals to 2, in Z when the identifier equals to 1. The subset of shared variables whose information is homogeneous between the databases are now stored from the fourth column to the last one. Their identities are available in the output object called REMAINING_VAR.

The merge_dbs function can handle incomplete shared variables by using the impute argument. This option allows to keep the missing information unchanged (the choice by default, and also the case in this example), to work with complete cases only, to do fast multivariate imputations by chained equations approach (the function integrates the main functionalities of the mice function of the mice package), or to impute data using the dimensionality reduction method introduced in the missMDA package (see the pdf manual for details).

Selection of matching variables using select_pred

In data fusion projects, a selection of shared variables (also called matching variables) appears as an essential step for two main reasons. First, the proportion of shared variables X between the two databases can be important (much higher than three variables) and keeping a pointless part of variability between variables could strongly reduce the quality of the fusion. Second, this selection greatly influences the quality of the predictions regardless of the matching technique which is chosen a posteriori [START_REF] Adamek | Fusion: combining data from separate sources[END_REF]. The specific context of data fusion is subject to the following constraints:

1. Y , Z and X are never jointly observed in database A or B, so the relationships between Y and X, and between Z and X must be investigated separately.

2. matching variables need to be good predictors of both target variables (outcomes) Y and Z [START_REF] Cohen | Statistical matching and microsimulation models, improving information for social policy decisions, the use of microsimulation modeling[END_REF], in the sense that they explain relevant parts of variability of the targets.

These particularities suppose that the investigations have to be done separately but in the same manner in both databases. In this way, a predictor which appears at the same time as highly correlated to Y and Z will be automatically selected for the fusion. On the contrary, predictors whose effects on Y and Z seem not obvious will be discarded. Additional recommended rules also emerge from literature:

• concerning the subset of variables that would predict only one of the two targets Y or Z, D 'Orazio et al. (2006) suggests a moderate parsimonious selection remaining too many predictors could complicate the fusion procedure.

• [START_REF] Cibella | How to choose the matching variables[END_REF] and [START_REF] Scanu | Recommendations on statistical matching[END_REF] suggests to select quality predictors with no error and just a small amount of missing data.

The function of the package dedicated to this task is select_pred. From the DB_READY database generated in the previous subsection, studying the outputs related to each database and produced by the following R commands assist users in selecting the best predictors: The quanti, nominal, and ordinal arguments requires vectors of column indexes related to the type of each variable of the input database. The ID argument specifies the column index of the database identifier. Y and Z expected the respective names of the target variables while OUT provides the target variable to predict (Y or Z).

To detect the subset of best predictors, select_pred studies the standard pairwise associations between each shared variable and the outcomes Y and Z, taken separately (by only varying the argument OUT), according to its type: for numeric and/or ordered factor variables, select_pred calculates the associations using rank correlation coefficients (Spearman) while the Cramer's V criterion [START_REF] Bergsma | A bias-correction for cramér's v and tschuprow's t[END_REF] is used for categorical variables and not ordered factors. The related ranking tables of top scoring predictors are available in two distinct output objects: cor_OUTC_num and cor_OUTC_cat. In our example, the corresponding results, for each database, are: The two first tables related to api29 highlights the strong associations between Y and the variables meals, ell, full and grad.sch in this order of importance, while the summary tables related to api35 highlights the strong associations between Z and the variables meals, grad.sch, ell and full.

### ASSOCATIONS
It is often not unusual to observe that one or more predictors are in fact linear combinations of others. In supervised learning areas, these risks of collinearity increase with the number of predictors, and must be detected beforehand to keep only the most parsimonious subset of predictors for fusion. To avoid collinearity situations, the result of a Farrar and Glauber (FG) test is provided [START_REF] Farrar | Multicollinearity in regression analysis: the problem revisited[END_REF]. This test is based on the determinant of the rank correlation matrix of the shared variables D (Spearman) and the corresponding test statistic is given by:

SF G = -n -1 - 1 6 (2k + 5) × ln(det(D))
where n is the number of rows and k is the number of covariates. Th null hypothesis supposes that SF G follows a chi square distribution with k(k -1)/2 dof, and its acceptation indicates an absence of collinearity. In presence of a large number of numeric covariates and/or ordered factors, the approximate Farrar-Glauber test, based on the normal approximation of the null distribution [START_REF] Kotz | Continuous multivariate distributions[END_REF]) can be more adapted and the statistic test becomes:

√ 2SF G - √ 2k -1
Users will note that these two tests can often appear highly sensitive in the sense that they tend to easily conclude in favor of multicollinearity. Thus, it is suggested to consider these results as indicators of collinearity between predictors rather than an essential condition of acceptability. The results related to this test are stored in the FG_test object of the select_pred output.

In presence of collinearity, select_pred tests the pairwise associations between all the potential predictors according to their types (Spearman or Cramér's V ). The thres_num argument fixed the threshold beyond which two ranked predictors are considered as redundant while thres_cat fixed the threshold of acceptability for the Cramér's V criterion in the subgroup of factors. In output, the results stored in the collinear_PB object permit to identify the corresponding variables. In our example, we observe: The FG test warns the user against the risks of collinearity between predictors, and the function notably detects strong collinearities between the variables meals, ell, full in the api29 (less strong trends in api35): an information that have to be taken into account during the selection. The part of predictors finally kept for the data fusion must be small to improve its quality. When the initial number of shared variables is not too important as in this example, choosing the best candidates between groups of redundant predictors can be made manually by selecting highest ranked predictors in the summary tables previously described. In this way, the variable meals could be preferred to ell, and full, while the variable stype could be dropped. Consequently, a possible final list of predictors could be only composed of the variables meals and grad.sch.

### DETECTION OF REDUNDANT PREDICTORS ## Results for the api29 dataset ----- ##
When the number of predictors is important, or when users prefer that an automatic process performs the variable selection, a random forest procedure can also be used via the select_pred function. Random forest approaches [START_REF] Breiman | Random forests[END_REF] are here particularly convenient [START_REF] Grajski | Classification of eeg spatial patterns with a tree-structured methodology: Cart[END_REF] for multiple reasons: it works fine when the number of variables exceeds the number of rows, whatever the types of covariates, it allows to deal with non linearity, to consider correlated predictors, ordered or not ordered outcomes and to rank good candidate predictors through an inbuilt criterion: the variable importance measure.

In few words, random forest processes aggregates many CART models [START_REF] Breiman | Classification and regression trees[END_REF] with RF_ntree bootstrap samples from the raw data source and averaging accuracies from each model permits to reduce the related variances and also the errors of prediction. A standard random forest process provides two distinct measures of importance of each variable for the prediction of the outcome, the Gini importance criterion and the permutation importance criterion, which depends on the appearance frequency of the predictor but also on its place taken up in each tree. For more details about random forest theory, user can consult [START_REF] Breiman | Random forests[END_REF] and/or the pdf manual of the randomForest package. Strobl et al. (2009) suggests that the permutation importance criterion, which works with permuted samples (subsampling without replacements) instead of bootstrap ones, is particularly conve-nient with uncorrelated predictors, but must be replaced by a conditional permutation measurement in presence of strong correlations. select_pred provides these assessments by integrating the main functionalities of the cforest and varimp functions of the package party (Hothorn et al., 2006b;[START_REF] Zeileis | Model-based recursive partitioning[END_REF]Hothorn et al., 2006a;[START_REF] Strobl | Bias in random forest variable importance measures: Illustrations, sources and a solution[END_REF][START_REF] Strobl | Conditional variable importance for random forests[END_REF]. However, these measurements must be used with caution, by accounting the following constraints:

• the Gini importance criterion can produce bias in favor of continuous variables and variables with many categories. This criterion is thus not available in the function.

• the permutation importance criterion can overestimate the importance of highly correlated predictors and therefore redundant predictors will be discarded beforehand using the first steps of the process integrated in the function.

The select_pred function allows to proceed with different scenarios according to the type of predictors (Table 4 can help users to choose). The first one consists in boiling down to a set of categorical variables (ordered or not) by categorizing all the continuous predictors using the dedicated argument (convert_num and convert_clss) and to work with the conditional importance assessments that directly provide unbiased estimations (by setting the RF_condi argument to TRUE). Users can consult (Hothorn et al., 2006b) for more details about the approach and consult the pdf manual of the package for details about the related arguments of the select_pred function. This approach does not take into account incomplete information, so that the method will be applied to complete data only (incomplete rows will be temporarily removed from the study). It is nevertheless possible to impute missing data beforehand by using dedicated pre-existing R packages like mice [START_REF] Van Buuren | Multiple imputation of missing blood pressure covariates in survival analysis[END_REF] or by using the imput_cov function provided in this package.

The second possible scenario (always usable in presence of mixed type predictors), consists in the execution of a standard random forest procedure after taking care to rule out collinearity issues by first selecting unique candidates between potential redundant predictors (in this case, the discarded predictors are stored in the drop_var output object). This is the default approach used by select_pred as soon as the RF_condi argument is set to FALSE while RF is set to TRUE. This scenario can work in presence of incomplete predictors. By constructing, note that results from random forest procedures stay dependent on successive random draws carried out for the constitution of trees, and it is so suggested to check this stability by testing different random seeds (RF_SEED argument) before concluding.

The R command previously written provides automatically the results related to the second approach as additional results. The results from the two datasets show here the permutation importance estimates of each variable ranked in order of importance and expressed as percentage, after resolving collinearity problems: The results confirm that the variable meals appeared as the best predictor of the target variables in api29 and api35 respectively. The variable ell is not present in the first list (see RF_PRED from step2a) because the variables meals and ell has been detected as strongly collinear (according to the initial chosen threshold) and so select_pred keep the best choice between the two predictors: meals (the reasoning is the same for full which disappeared from the list).

R> step2a$RF_PRED #
The ell variable has been kept in the second list (not found as collinear enough to be removed here) and appears moreover as a good predictor of apicl_1999 in api35. Nevertheless its potential collinearity problem with meals encourages us not to keep it for the rest of the study. According to this discussion, we finally keep meals, stype and grad.sch which combines the advantages of being good predictors for the two target variables while not presenting major problems of collinearity between them.

The following synthetic database (called here bdd_ex) is now ready for the prediction of the missing API scores in api29, api35 or both, using the function OT_outcome or OT_joint: R> bdd_ex = step1$DB_READY [, c(1:3,10:12)

Data fusion using OT_outcome

The OT_outcome function provides individual predictions of Z in A (and/or Y in B) by considering the recoding problem involving optimal transportation of outcomes. In a first step, the aim of the function is so to determine γ from (2) while this estimate is used in a second step to provide the predictions. The main input arguments of the function and the output values are respectively described in Tables 5 and6. The arguments datab, index_DB_Y_Z, quanti, nominal, ordinal, and logic are not optional and must be carefully completed before each execution. A unique synthetic data.frame made of two overlayed databases (called A and B for example) (see Table 3) is expected as datab argument. If this data.frame corresponds to the output objects DB_USED or DB_READY of the select_pred or merge_dbs functions respectively, then the expected object has the required structure. Otherwise users must be sure that their data.frames are made up of two overlayed databases with at least 4 variables as described in the subsection Optimal transportation of outcomes applied to data recoding. The order of the variables have no importance in the raw database but will have to be specified a posteriori in the index_DB_Y_Z argument if necessary.

The subset of remaining predictors used for fusion may requires prior transformations according to the distance function chosen in input by the user. This process is fixed by setting the dist.choice argument. The distances actually implemented in the function are: the standard Euclidean ("E") and Manhattan ("M") distances, the Hamming distance ("H", for binary covariates), and the Gower distance ("G" sometimes preferred with mixed variables). Automatic transformations prior to the use of each distance function are summarized in Table 7.

The first version of the OT algorithm described in [START_REF] Gares | On the use of optimal transportation theory to recode variables and application to database merging[END_REF] was tested on numeric coordinates extracted from a factor analysis of mixed data (FAMD) fitted on mixed raw covariates [START_REF] Pages | Analyse factorielle multiple appliquée aux variables qualitatives et aux données mixtes[END_REF]. This transformation is here available by setting the FAMD.coord argument to "YES". The minimal percentage of information explained by the FAMD is also fixable using the FAMD.perc argument. The OT_outcome functions proposes four types of models for the prediction of Y in B (and/or) Z in A, according to the values of the method and maxrelax arguments:

• When maxrelax = 0 and indiv.method = "sequential" (default options), the fitted model corresponds to the outcome algorithm described by P0 n in equation ( 6). Assuming that Y and Z in A follow the same distribution as Y and Z in B respectively (assumption 1, page 2), this related algorithm derives the joint distribution of Y and Z in A (respectively in B) in a first step, and uses in a second step, a nearest neighbor procedure to predict missing values of Z in A (resp. Y in B). This algorithm calculates averaged distances between each subject from A (resp. B) and subgroups of subjects from B (resp. A) having same levels of Z in B (resp. Y in A). These calculations can be done using all subjects of each subgroups (by default, percent.knn = 1) or only restricted parts of them (percent.knn< 1).

• When maxrelax > 0 and indiv.method = "sequential", assumption 1 is alleviated by relaxing the constraints on marginal distributions and an r-outcome algorithm (with relaxation) is fitted. In this case, the maxrelax argument corresponds to the αn parameter in (18).

• When maxrelax = 0 and indiv.method = "optimal", the second step of the original algorithm (nearest neighbor procedure) is replaced by a linear optimization problem: searching for the individual predictions of Z in A (resp. Y in B) by minimizing the total of the distances between each individual of A and individuals of each levels of Z in B.

• When maxrelax > 0 and indiv.method = "optimal", the constraints on marginal distributions of the previous model are also relaxed and the function fits an r-outcome algorithm with relaxation. For these three last situations, the corresponding r-outcome algorithm is so described by P0-R n (18).

When maxrelax > 0, it is recommended to calibrate the maxrelax parameter by testing different values according to the stability of individual predictions. In output, the gamma_A and gamma_B matrices correspond to the estimates of the joint distributions γ of (Y A , Z A ) and (Y B , Z B ) re- spectively (equation ( 2)). Moreover, a posteriori estimates of conditional distributions probabilities (Z A |Y A , X A ) and (Y B |Z B , X B ) (see equation 10) are also provided in two lists called estimatorZA, and estimatorYB respectively and the completed databases A and B are stored in the DATA1_OT and DATA2_OT objects. In particular, the individual predictions are stored in the OTpred column of these data.frames. Moreover, the profile object is a data.frame that stores the profile of covariates encountered in the two data sources while the res_prox object stored all the distance computations that will be used in the validation step (users can consult details of the proxim_dist function of the pdf manual). The computation of conditional probabilities implies to define the part of individuals considered as neighbors of each encountered profile of covariates. The argument prox.dist fixes this threshold for each profile, following the decision rule for A (and respectively for B): A subject i of A (or a statistical unit corresponding to a row of A) will be considered as neighbor of a given profile of shared variables C as soon as:

dist(subjecti, C) < prox.dist × max k=1,...,n A dist(subject k , C)
When the number of shared variables is small and all of them are categorical (optimal situation), it is suggested to set the prox.dist argument to 0. Finally, using the which.DB argument, users can choose to estimate the individual predictions of Y and Z in the two databases (which.DB = "BOTH") or only one of the both (which.DB = "A" for the predictions of Z in A or which.DB = "B" for the predictions Y in B).

From the data.frame bdd_ex built previously, we use the OT_outcome function to illustrate the prediction of the target variable apicl_2000 in the api35 dataset via an outcome algorithm and using an Euclidean distance function: --------------------------------------In bdd_ex, all variables except the first one (DB: the database identifier) are or can be considered as ordinal factors, and the quanti and ordinal arguments are filled in accordingly. For the prediction of apicl_2000 in the api35 dataset (the steps would be the same for the prediction of apicl_1999 in api29 by setting which.DB = "A" or "BOTH"), the optimal transportation theory determines a map γ that pushes, in the distribution of apicl_1999 forward to the distribution of apicl_2000 in the database api35. In this case, γ B is an estimator of the joint distribution of apicl_2000 and apicl_1999 in api35. The estimation of γ B is available in the gamma_B output object while all the profiles of predictors met in the two databases are listed in the profile object: 0.22248804 0.0000000 0.00000000 0.0000000 #(600-800] 0.02889318 0.2486188 0.15311005 0.0000000 #(800-1000] 0.00000000 0.0000000 0.09550874 0.2513812

R> outc1$gamma_B # G1 G2 G3 G4 #

R> outc1$profile[1,]

# the first profile ID grad.sch meals stype P_1 3 1 3

The output object estimatorYB is a list that corresponds to the estimations of the conditional probabilities of apicl_2000 in the api35 dataset for a given profile of predictors. For example, the conditional probabilities related to the first profile of predictors are:

R> outc1$estimatorYB[1,,] # conditional probabilities (1st profile) [,1] [,2]
[,3] G1 0.3333333 0.3333333 0.3333333 G2 0.3333333 0.3333333 0.3333333 G3 0.0000000 0.0000000 1.0000000 G4 0.0000000 0.0000000 1.0000000

According to these results, we can conclude that: for a subject from api35 with a related profile of predictors P_1 and whose levels of apicl_1999 is 'G1', the probability for apicl_2000 to be predicted '[200,600]' is about 0.63. Finally, the individual predictions of apicl_2000 in api35 are available in the DATA2_OT object corresponding to the OTpred column: The outcome algorithm uses here a nearest neighbor procedure to assign the individual predictions from the estimation of γ which can be a drawback as described in [START_REF] Gares | On the use of optimal transportation theory to recode variables and application to database merging[END_REF]. The r-outcome algorithm proposes an enrichment that directly assigns the individual predictions of apicl_2000 from the estimation of γ, without using the nearest neighbor approach. In this case, a linear optimization problem is solved to determine the individual predictions that minimize the sum of the individual distances in api35 having the modalities of the target apicl_2000 in api29. The corresponding R command is: R> R_outc3 = OT_outcome(bdd_ex, quanti = 1, ordinal = 2:6, dist.choice = "E" , indiv.method = "optimal", which.DB = "B")

Moreover, the outcome algorithm assumes that the distribution of apicl_2000 is identically distributed in api29 and api35 (assumption 1 from the subsection 2.2) which can appear as a strong hypothesis to hold in many situations. To overcome this constraint, the R_OUTCOME algorithm also allows to relax the assumption 1 by adding a relaxation parameter (18) in the optimization system. This relaxation can be done by varying the maxrelax arguments of the function:

### R-OUTCOME algorithm: optimal assignments + relaxation parameter = 0.4 R> R_outc4 = OT_outcome(bdd_ex, quanti = 1, ordinal = 2:6, dist.choice = "E", indiv.method = "optimal", maxrelax = 0.4, which.DB = "B")

The running times of these two models can take few minutes. The quality of these models will be compared in Tables 8, 9 and 10 ( subsection Validation of the data fusion using verif_OT).

Data fusion using OT_joint

The OT_joint function provides individual predictions of Z in A (and/or Y in B) by considering the recoding problem as an optimal transportation problem of covariates and outcomes, that pushes the conditional distribution of (Y A |X A ) forward to the conditional distribution of (Z A |X A ) (and conversely (Z B |X B ) forward to the conditional distribution of (Y |X)). The aim is to determine γ from (7). Because joint distributions of outcomes and covariates are now mapped together (it was not the case with the outcome family of algorithms), it is not required for target variables to be equally distributed (assumption 1).

The call of OT_joint was thought to propose globally the same structure as those of the function OT_outcome as described in Tables 5 and6 with few differences described below. joint and r-joint algorithms are usable via OT_joint for solving the recoding problem depending on the values related to the maxrelax and lambda_reg arguments:

• When maxrelax = 0 and lambda.reg = 0 (default values), the fitted model corresponds to the joint algorithm described by Pn in equation ( 7).

• When at least one of these two arguments differs from 0, the r-joint algorithm is called. Using r-joint, it is so possible to relax constraints on marginal distributions (maxrelax > 0) and/or add an eventual more or less strong L1 regularisation term among the constraints of the algorithm by filling the lambda_reg argument. maxrelax parameter correspond to parameter αn in ( 16) and lambda.reg parameter correspond to parameter λ in (17).

When maxrelax > 0 and/or lambda.reg > 0, it is recommended to calibrate these two parameters by testing many different values and so studying the stability of the related individual predictions. Nevertheless, by default or as a starting point, [START_REF] Gares | On the use of optimal transportation theory to recode variables and application to database merging[END_REF] suggests the use of default values determined from simulated databases: 0.1 for the regularization parameter (lambda.reg) and 0.4 for the relaxation one (maxrelax). Finally, the output objects are similar to those previously described in the OT_outcome function.

The implementation of the function is very similar to that of OT_outcome, and applied to our example, the R code related to a joint algorithm is: R> outj1 = OT_joint(bdd_ex, nominal = 1, ordinal = c(2:6), dist.choice = "E" , which.DB = "B") distribution estimated from two variables which are never jointly observed [START_REF] Kadane | Some statistical problems in merging data files[END_REF][START_REF] Rodgers | An evaluation of statistical matching[END_REF]. When the study aims to identify estimators of interests that improve the understanding of the relationships between variables in the different databases (called macro approach in D 'Orazio et al. (2006)) without going until the creation of a complete and unique dataset, uncertainty analyses are usually proposed to estimate the sensitivity of the assessments [START_REF] Rässler | Statistical matching: A frequentist theory, practical applications, and alternative Bayesian approaches[END_REF]. On the contrary, if the objective is to create a synthetic dataset where the information is available for every unit, the use of auxiliary information or proxy variables must be privileged to assess the quality of the results [START_REF] Paass | Statistical match: evaluation of existing procedures and improvements by using additional information[END_REF]. In the absence of complementary information, [START_REF] Rubin | Statistical matching using file concatenation with adjusted weights and multiple imputations[END_REF] suggests to study the preservation of the relationships at best between the distributions of the target variables. In this way, the verif_OT function proposes tools to assess quality of the individual predictions provided by the algorithms previously introduced.

The first expected input argument of verif_OT is an 'otres' object from the OT_outcome or OT_joint functions. In our example, this function is firstly applied to the out_c1 model by following the R command: In output, the nb.profil value gives the number of profiles of predictors globally detected in the two databases (27). In the example, a profile will be characterized by a combination of values related to the three predictors kept: grad.sch, meals and stype.

The first step of the function is dedicated to the study of the proximity between the two target variables. Therefore, standard criteria (Cramer's V and Spearman's rank correlation coefficient) are used to evaluate the pairwise association between Y and Z, globally or in one of the two databases only, according to the predictions provided by the output of OT_outcome or OT_joint. Only the database api35 was completed in the example (because which.DB = "B" as argument of OT_outcome) and stored in the DATA2_OT object, therefore, the criteria compare here the proximity distribution between the predicted values of apicl_2000 (Y) and the observed value of apicl_1999 (Z) in the database api35 (B). Regarding independence or a small correlation between Y A and Z A (or Y B and Z B ) must question about the reliability of the predictions especially when Y and Z are supposed to summarize a same information. In the example, whatever the criteria used, we notice via the res.prox object, a strong association between the two target variables in api35 which reassures about the consistency of the predictions. The related confusion matrix between the predicted values of apicl_2000 (Y) and the observed value of apicl_1999 (Z) is stored in the conf.mat object.

Second, the function proposes an optional tool (by setting group.clss= TRUE) which evaluates the impact of grouping levels of one factor on the association of Y and Z. When users have initially no information about one of the two encodings of interest, this approach can be particularly useful to detect ordinal from nominal ones and its principle is as follow. Assuming that Y ∈ Y, and Z ∈ Z (Y and Z are the respective levels) and that |Y| ≥ |Z|. From Y , successive new variables Y ′ ∈ Y ′ are built, as soon as Y ′ is a partition of Y such as |Y ′ | = |Z| (and inversely for Z if the levels of Z is the greatest). The related associations between Z and Y ′ (with now equal number of levels) are then studied using: Cramer's V, rank correlation, Kappa coefficient and confusion matrix and the results are stored in a It appears from these results that grouping the levels G2 and G3 of apicl_1999 (Z) strongly improves the association of this variable with apicl_2000 (Y) (the error rate of the confusion matrix varies from 43.4 to 13.3). Moreover the structure of conf.mat confirms that the encoding of apicl_1999 seems to be ordinal.

The third step of the quality process integrated in the function corresponds to the study of the Hellinger distance [START_REF] Liese | Statistical decision theory[END_REF]). This distance function is used as a measure of the discrepancies between the observed and predicted distributions of Y (L(Y A ) versus L( Y B )) and/or (L( Z A ) versus L(Z B )). For Y and Z, the definition of the distance is respectively:

dist hell (Y A , Y B ) = 1 2 y∈Y µ Y A n,y -µ Y B n,y 2 ) and dist hell (Z B , Z A ) = 1 2 z∈Z µ Z B n,z -µ Z A n,z 2
)
where µ Y B n,y and µ Z A n,z corresponds to the empirical estimators of µ Y B and µ Z A respectively. The Hellinger distance varies between 0 (identical) and 1 (strong dissimilarities) while 0.05 can be used as an acceptable threshold below which two distributions can be considered as similar. When outcome models are applied on datasets, this criterion shows that predictions respect assumption 1. It can also be used to determine the best relaxation parameter of an r-outcome model. On the contrary, there is no need to interpret this criterion when an r-joint model is applied, because in this case, the assumption 1 is not required. Results related to this criterion are stored in the hell object: R> verif_outc1$hell YA_YB ZA_ZB Hellinger dist. 0.008 NA With a p-value equals to 0.08, the assumption 1 is difficult to hold here and must challenge the user on the need to switch the algorithm or to improve the current one by adding relaxation and/or regularization parameters.

Finally, the verif_OT function uses the mean and standard deviance of the conditional probabilities P(Z = ẑi|Y = yi, X = xi) estimated by each model, as indicators of the stability of the individual predictions (provided that stab.prob = TRUE). It is nevertheless possible that conditional probabilities are computed from too few individuals (according to the frequency of each profile of shared variables met), to be considered as a reliable estimate of the reality. To avoid this problem, trimmed means and standard deviances are suggested by removing these specific probabilities from the computation, using the min.neigb parameter. In output, the results related to this last study are stored in the res.stab object:

R>

Conclusion and perspectives

To our knowledge, OTrecod is the first R package that takes advantage of the optimal transportation theory (Monge ( 1781)) in the context of data fusion and the comparative study of methods described in [START_REF] Gares | Regularized optimal transport of covariates and outcomes in data recoding[END_REF] underlines the promising performances of these approaches.

For more details and examples about the functions, users are invited to consult the "ARTICLES" section of the related website https://otrecoding.github.io/OTrecod/.

Drawbacks

The functions of OTrecod only apply to a specific data fusion context where there is no overlapping part between the two raw data sources. This package does not deal with record linkage which is already the focus of extensive works [START_REF] Sariyar | The recordlinkage package: Detecting errors in data[END_REF]. This package is not adapted when the target variables (outcomes) of the two databases are continuous with an infinite number of values (like weights in grams with decimals for example). Moreover, if the data fusion requires only one shared variable to work, the quality of the fusion depends on the presence of a subset of good shared predictors in the raw databases. Finally, if the function OT_outcome allows all types of predictors, the current version of the function OT_joint imposes categorical matching variables only (scale variables are allowed): this restriction should be relaxed in a next version.

Perspectives

A number of more advanced research still need further investigation:

1. the possibility of extending the OT algorithm for recoding variables to multidimensional frameworks.

2. the stability of the algorithm when the matching variables are incomplete and the non-response processes are missing at random or not.

3. the contribution of calibration techniques in the quality process assessment [START_REF] Deming | On a least squares adjustment of a sampled frequency table when the expected marginal totals are known[END_REF] 4. the creation of a tuning function which defines a grid search to find the optimal combination of parameters related to the OT algorithms, could be added in the future versions of the package.

  R> step1 = merge_dbs(DB1 = api29, DB2 = api35, NAME_Y = "apicl_2000", NAME_Z = "apicl_1999", row_ID1 = 1, row_ID2 = 1, ordinal_DB1 = c(2:3, 8:12), ordinal_DB2 = c(2:3, 8:12))DBS MERGING in progress. Please wait ...DBS MERGING OK -----------------------SUMMARY OF DBS MERGING:Nb of removed subjects due to NA on targets: 0(0%) Nb of removed covariates due to differences between the 2 bases: 1 Nb of remained covariates: 9 Imputation on incomplete covariates: NO

#

  For the dataset api29 --------------R> step2a = select_pred(step1$DB_READY, Y = "Y", Z = "Z", ID = 1, OUT = "Y", quanti = c(4,6), nominal = c(1,5,7), ordinal = c(2,3,8:12), thresh_cat = 0.50, thresh_num = 0.70, RF_SEED = 3011) # For the dataset api35 --------------R> step2b = select_pred(step1$DB_READY, Y = "Y", Z = "Z", ID = 1, OUT = "Z", quanti = c(4,6), nominal = c(1,5,7), ordinal = c(2,3,8:12), thresh_cat = 0.50, thresh_num = 0.70, RF_SEED = 3011)

  R> outc1 = OT_outcome(bdd_ex, quanti = 1, ordinal = 2:6, dist.choice = "E", indiv.method = "sequential", which.DB = "B")---------------------------------------OT PROCEDURE in progress ... ------------------------------------

#

  ## Quality criteria for outc1 (OUTCOME model) R> verif_outc1 = verif_OT(outc1, group.class = TRUE, ordinal = FALSE, stab.prob = TRUE, min.neigb = 5) ### First results related to outc1: R> verif_outc1$nb.

Table 1 :

 1 Formulation of a recoding problem for 2 distinct databases A and B

	1

Table 2 :

 2

A brief description of the main functions of OTrecod

All the intermediate functions integrated in the OT_outcome and OT_joint functions (proxim_dist,

Table 4 :

 4 Completing the RF_condi argument according to predictors

			For the api29 dataset
	meals	stype grad.sch	awards acs.core
	71.5529 11.6282 11.1181	5.4674	0.2334
	R> step2b$RF_PRED	# For the api35 dataset
	Possible scenarios	Correlation	State of the	Incomplete
				between	RF_condi	information
				predictors	argument
	Same type predictors	NO	FALSE	Allowed
	Same type predictors	YES	TRUE	Not allowed
	Different type predictors	NO	FALSE	Allowed
	Different type predictors	YES	TRUE	Not Allowed

Table 5 :

 5 Main arguments of the OT_outcome and OT_joint functions. (*: NULL as default value)

	Argument	OT_outcome OT_joint Description (default value)
	datab	✓	✓	Data.frame in the expected structure
	index_DB_Y_Z	✓	✓	Indexes of the ID,Y, and Z columns
				(1,2,3)
	nominal	✓	✓	Column indexes of nominal variables
				(*)
	ordinal	✓	✓	Col. indexes of ordinal variables (*)
	logic	✓	✓	Col. indexes of boolean variables (*)
	quanti	✓		Col. indexes of quantitative variables
				(*)
	convert.num	✓	✓	Col. indexes of the quantitative vari-
				ables to convert (* or =quanti in
				OT_joint)
	convert.clss	✓	✓	Corresponding numbers of desired
				classes for conversion (*)
	which.DB	✓	✓	Specify the target variables to com-
				plete: Both or only one (BOTH)
	solvR	✓	✓	Choice of the solver to solve the opti-
				mization problem (glpk)
	dist.choice	✓	✓	Distance function (Euclidean). See Ta-
				ble 7
	percent.knn	✓	✓	Ratio of closest neighbors involved int
				the computations (1)
	indiv.method	✓		Type of individual predictions (sequen-
				tial) for outcome and r-outcome al-
				gorithms
	maxrelax	✓	✓	Adding of a relaxation parameter (0)
	lambda.reg		✓	Adding of regularization parameter (0)

Table 6 :

 6 Values of the OT_outcome and OT_joint functions. (*: NULL if not required)

	Value	OT_outcome OT_joint Description
	time_exe	✓	✓	Running time of the algorithm
	gamma_A	✓	✓	Estimation of the joint distribution of
				(Y, Z) for the prediction of Z in A (*)
	gamma_B	✓	✓	Estimation of the joint distribution of
				(Y, Z) for the prediction of Y in B (*)
	profile	✓	✓	The list of detected profiles of covari-
				ates
	res.prox	✓	✓	A list that provides all the information
				related to the estimated proximities be-
				tween profiles and groups of profiles
	estimator_ZA	✓	✓	Estimates of the probability distribu-
				tion of Z conditional to X and Y in
				database A (*)
	estimator_YB	✓	✓	Estimates of the probability distribu-
				tion of Y conditional to X and Z in
				database B (*)
	DATA1_OT	✓	✓	The database A fully completed (if re-
				quired in input by the which.DB argu-
				ment)
	DATA2_OT	✓	✓	The database B fully completed (if re-
				quired in input by the which.DB argu-
				ment)

  table called res.grp. The corresponding outputs of the example are:

	predY	G1 G2 G3 G4 Sum
	[200-600]	81	1	1	1 84
	(600-800]	10 89 55	1 155
	(800-1000]	0	0 34 89 123
	Sum	91 90 90 91 362
	R> verif_outc1$res.grp
	grp levels Z to Y error_rate Kappa Vcramer RankCor
	G1/G2 G3/G4		13.3 0.794	0.83	0.877
	G1/G2/G3 G4		19.1 0.714	0.78	0.839
	G1 G3/G2/G4		28.2 0.593	0.68	0.624
	G1 G2/G3/G4		37.6 0.457	0.64	0.813
	G1 G4/G2/G3		43.4 0.374	0.59	0.115
	G1/G2 G4/G3		43.4 0.326	0.64	0.574
	R> verif_outc1$conf.mat
		Z

Table 12 :

 12 Confusion matrices in the api35 dataset for the models (a) R_outc4 and (b) R_outj1

	verif_outc1$eff.neig
	Nb.neighbor Nb.Prob
	1	14
	2	18
	3	18
	4	28
	5	20
	6	6
	7	14
	8	16
	9	27
	10	20
	R> verif_outc1$res.stab
	N min.N mean	sd
	2nd DB 284	5 0.968 0.122

a parenthesis removed

Acknowledgments

The authors would especially thank Pierre Navaro (CNRS UMR 6625) for their advices during the implementation of the OTrecod package. This research has received the help from Region Occitanie Grant RBIO-2015-14054319 and Mastodons-CNRS Grant.

---------------------------------------OT JOINT PROCEDURE in progress ... --------------------------------------- For relaxing the constraints stated on marginal distributions, we use the maxrelax argument that corresponds to varying the α parameter of a r-joint algorithm (see 17) and a parameter of regularization λ can be simply added to the previous model as follows:

### R-JOINT algorithm (relaxation parameter = 0.4) R> R_outj1 = OT_joint(bdd_ex, nominal = 1, ordinal = c(2:6), dist.choice = "E", maxrelax = 0.4, which.DB = "B")

dist.choice = "E", maxrelax = 0.4, lambda.reg = 0.1, which.DB = "B")

As previously, these models will be compared in the next subsection.

Validation of the data fusion using verif_OT

Assessing the quality of individual predictions obtained through statistical matching techniques can be a complex task notably because it is legitimate to wonder about the true reliability of a joint

The first result shows that 14 individual predictions among 362 (the total number of rows in api35) have been assigned to subjects that exhibits a unique combination of predictors and outcome. From these subjects, it would be obviously overkill to draw conclusions about the reliability of the predictions provided by the model. We therefore decide to fix here a threshold of 5 below which it is difficult to extract any information related to the prediction. From the remaining ones, we simply perform the average (0.968) which could be interpreted as follows: When the fitted model is confronted with a same profile of predictors and a same level of apicl_1999, more than 96 times of a hundred, it will return the same individual prediction.

We run a total of 11 models to determine the optimal one for the prediction of apicl_2000 in api35. Every arguments from the syntax of the OT_outcome and OT_joint functions described on pages 15, 16 and 17, stay unchanged with the exception of indiv.method, maxrelax, and lambda.reg which vary. The values linked to each criterion of the verif_OT function are summarized in Tables 8, 9 and 10. The syntax related to verif_OT stay unchanged for each model (notably same min.neigb arguments). Note that comparisons of stability predictions between outcome and joint models impose that the prox.dist argument of the OT_outcome function is fixed to 0.

Model

Summary table

From these results, we can conclude that:

• whatever the algorithm used (outcome or joint), adding a relaxation parameter improves here the association between the target variables (see Table 8).

• according to the results of Table 9, the outc1 and R_outc3 models are not optimal because the observed values from the Hellinger criterion reflects a violation of assumption 1 (see Table 9). Therefore, for these two models, a relaxation parameter must be added or increased respectively.

• conversely, adding a relaxation parameter that is too high seems to affect the quality of the association. In this example, 0.4 seems to be a good value for this parameter whatever the algorithm used. Among the remaining models (R_outc1, R_outc4, R_outj1, and R_outj4), R_outj1 and R_outj4 seem to be those with the most stable predictive potential (Table 10).

Moreover, R_outc4 appears here as the best model from the tested outcome algorithms.

• adding a regularization parameter to the R_outj1 model caused a decrease in the stability of the predictions (see the value of R_outj4) and we thus conclude in favor of R_outj1 as best model among those tested in the joint family of algorithms.

• Table 11 shows that the three remaining models counted between 86 an 91% of common predictions and these last result also reassures the user about the quality of the provided predictions. The confusion matrices related to R_outc4 and R_outj1 are described in Table 12 and seems to confirm that the encoding of apicl_2000 in three groups, could simply correspond to the grouping of levels G2 and G3 of the api_cl1999 variable.

Model

Finally, note that the running time of each model took less than 15 seconds with R version 4.0.3 for Windows (10 Pro-64bits/ Process Intel 2.66 GHz). = verif_outj1$res.prox , R_outj1 = verif_R_outj1$res.prox, R_outj2 = verif_R_outj2$res.prox , R_outj3 = verif_R_outj3$res.prox, R_outj4 = verif_R_outj4$res.prox ) R> res.prx = data.frame(Model = c("outc1","R_outc2","R_outc3","R_outc4", "R_outc5","R_outc6","outj1", "R_outj1", "R_outj2","R_outj3","R_outj4"), Type = c("OUTCOME",rep("R-OUTCOME",5),"JOINT", rep("R-JOINT",4)), Relax = c(0,0.4,0.6,0,0.4,0.6,0,0.4,0.6,0.8,0.4), Regul = c(rep(0,10),0.1), res.prx) R> row.names(res.prx) = NULL; head(res.prx,3) R_outc2 = verif_R_outc2$hell, R_outc3 = verif_R_outc3$hell, R_outc4 = verif_R_outc4$hell, R_outc5 = verif_R_outc5$hell, outj1 = verif_outj1$hell , R_outj1 = verif_R_outj1$hell, R_outj2 = verif_R_outj2$hell, R_outj3 = verif_R_outj3$hell, R_outj4 = verif_R_outj4$hell ) res.helld = data.frame(res.prx[,1:4], res.helld) row.names(res.helld) = NULL; res.helld #-----## Stability of the prediction: Summary Table 10 # same R code as for Summary Table 8, changing res.prox by res.stab #----## Ratio of common predictions: Table 11 R> stoc = list(outc1$DATA2_OT$OTpred , R_outc1$DATA2_OT$OTpred, R_outc4$DATA2_OT$OTpred, outj1$DATA2_OT$OTpred ,