Martin Bromberger
email: mbromber@mpi-inf.mpg.de

A Reduction from Unbounded Linear Mixed Arithmetic Problems into Bounded Problems

We present a combination of the Mixed-Echelon-Hermite transformation and the Double-Bounded Reduction for systems of linear mixed arithmetic that preserve satisfiability and can be computed in polynomial time. Together, the two transformations turn any system of linear mixed constraints into a bounded system, i.e., a system for which termination can be achieved easily. Existing approaches for linear mixed arithmetic, e.g., branch-and-bound and cuts from proofs, only explore a finite search space after application of our two transformations. Instead of generating a priori bounds for the variables, e.g., as suggested by Papadimitriou, unbounded variables are eliminated through the two transformations. The transformations orient themselves on the structure of an input system instead of computing a priori (over-)approximations out of the available constants. Experiments provide further evidence to the efficiency of the transformations in practice. We also present a polynomial method for converting certificates of (un)satisfiability from the transformed to the original system.

Introduction

Efficient linear arithmetic decision procedures are important for various independent research lines, e.g., optimization, system modeling, and verification. We are interested in feasibility of linear arithmetic problems in the context of the combination of theories, as they occur, e.g., in SMT solving or theorem proving.

The SMT and theorem proving communities have presented several interesting and efficient approaches for pure linear rational arithmetic [START_REF] Dutertre | A fast linear-arithmetic solver for DPLL(T)[END_REF] as well as linear integer arithmetic [START_REF] Bobot | A simplex-based extension of fouriermotzkin for solving linear integer arithmetic[END_REF][START_REF] Bromberger | Linear integer arithmetic revisited[END_REF][START_REF] Dillig | Cuts from proofs: A complete and practical technique for solving linear inequalities over integers[END_REF][START_REF] Griggio | A practical approach to satisfiability modulo linear integer arithmetic[END_REF]. SMT research also starts to extend into linear mixed arithmetic [START_REF] Christ | Cutting the mix[END_REF][START_REF] Dutertre | A fast linear-arithmetic solver for DPLL(T)[END_REF] because some applications require both rational and integer variables, e.g., planning/scheduling problems and verification of timed automata and hybrid systems.

We are interest in decision procedures for mixed arithmetic because of a possible combination with superposition [START_REF] Althaus | Superposition modulo linear arithmetic SUP(LA)[END_REF][START_REF] Baumgartner | Hierarchic Superposition with Weak Abstraction[END_REF][START_REF] Fietzke | Superposition as a decision procedure for timed automata[END_REF]. In the superposition context, arithmetic constraints are part of the first-order clauses. The problems are typically unbounded due to transformations that turn the input formula into a superposition specific input format. Since these problems are unbounded, the search space becomes infinite, which is the case where termination becomes difficult for most linear arithmetic approaches. Unbounded problems appear also in other areas of automated reasoning. Either because of bad encodings, necessary but complicating transformations, e.g., slacking (see Section 5), or the sheer complexity of the verification goal. Hence, efficient techniques for handling unbounded problems are necessary for a generally reliable combined procedure.

It is theoretically very easy to achieve termination for linear integer and mixed arithmetic because of so called a priori bounds. For example, the a priori bounds presented by Papadimitriou [START_REF] Papadimitriou | On the complexity of integer programming[END_REF] guarantee that a problem has a mixed solution if and only if the problem extended by the bounds |x i | ≤ 2n(ma) 2m+1 for every variable x i has a mixed solution. In these a priori bounds, n is the number of variables, m the number of inequalities, and a the largest absolute value of any integer coefficient or constant in the problem. By extending a problem with those a priori bounds, we reduce the search space for a branch-and-bound solver (and many other mixed arithmetic decision procedures) to a finite search space. So branch-andbound is guaranteed to terminate.

However, these bounds are so large that the resulting search space cannot be explored in reasonable time for many practical problems. One reason for the impracticability of a priori bounds is that they only take parameter sizes into account and not actually the structure of each problem. A priori bounds are not integrated in any state-of-the-art SMT solvers [START_REF] Barrett | CVC4. In CAV[END_REF][START_REF] Christ | SMTInterpol: an interpolating SMT solver[END_REF][START_REF] Cimatti | The Math-SAT5 SMT Solver[END_REF][START_REF] De Moura | Z3: An efficient SMT solver[END_REF][START_REF] Dutertre | Yices 2.2[END_REF] since they are no help in practice. As far as we know, none of the state-of-the-art SMT solvers use any method that guarantees termination for linear integer or mixed arithmetic.

In this paper, we present satisfiability preserving transformations that reduce unbounded problems into bounded problems. On these bounded problems, most linear mixed decision procedures become terminating, which we show on the example of branch-and-bound. Our reduction works by eliminating unbounded variables. First, we use the Double-Bounded reduction (Section 4) to eliminate all unbounded inequalities from our constraint system. Then we use the Mixed-Echelon-Hermite transformation (Section 3) to shift the variables of our system to ones that are either bounded or do not appear in the new inequalities and are, therefore, eliminated. With Corollary 2 & Lemma 13 we explain how to efficiently convert certificates of (un)satisfiability between the transformed and the original system. Our method is efficient because it is fully guided by the structure of the problem. This is confirmed by experiments (Section 5). We also show how to efficiently determine when a problem is unbounded (Lemma 10). This prevents our solver from applying our transformations on bounded problems.

An extended version of this paper is available on arXiv [START_REF] Bromberger | A reduction from unbounded linear mixed arithmetic problems into bounded problems[END_REF]. It contains an appendix, where we explain how to implement the presented procedures in an incrementally efficient way. This is relevant for the implementation of an efficient SMT theory solver. The extended version also contains several new examples as well as additional implementation tricks.

Preliminaries

While the difference between matrices, vectors, and their components is always clear in context, we generally use upper case letters for matrices (e.g., A), lower case letters for vectors (e.g., x), and lower case letters with an index i or j (e.g., b i , x j ) as components of the associated vector at position i or j, respectively. The only exceptions are the row vectors a T i = (a i1 , . . . , a in ) of a matrix A = (a 1 , . . . , a m ) T , which already contain an index i that indicates the row's position inside A. We also abbreviate the n-dimensional origin (0, . . . , 0) T as 0 n . Moreover, we denote by piv(A, j) the row index of the pivot of a column j, i.e., the smallest row index i with a non-zero entry a ij or m + j if there are no non-zero entries in column j.

A system of constraints Ax ≤ b is just a set of non-strict inequalities

1 {a T 1 x ≤ b 1 , . . . , a T m x ≤ b m }
and the rational solutions of this system are exactly those points x ∈ Q n that satisfy all inequalities in this set. The row coefficients are given by A = (a 1 , . . . , a m ) T ∈ Q m×n , the variables are given by x = (x 1 , . . . , x n ) T , and the inequality bounds are given by b = (b 1 , . . . , b m ) T ∈ Q m . Moreover, we assume that any constant rows a i = 0 n were eliminated from our system during an implicit preprocessing step. This is a trivial task and eliminates some unnecessarily complicated corner cases.

In this paper, we consider mixed constraint systems, i.e., variables are assigned a type: either rational or integer. Due to convenience, we assume that the first n 1 variables (x 1 , . . . , x n 1 ) are rational and the remaining n 2 variables (x n 1 +1 , . . . , x n ) are integer, where

n = n 1 + n 2 . A mixed solution is a point x ∈ (Q n 1 × Z n 2 ) that satisfy all inequalities in Ax ≤ b and we denote by M(Ax ≤ b) = {x ∈ (Q n 1 × Z n 2 ) : Ax ≤ b}
the set of mixed solutions to the system of inequalities Ax ≤ b. We sometimes need to relax the variables to be completely rational. Therefore, we denote by Q(Ax ≤ b) = {x ∈ Q n : Ax ≤ b} the set of rational solutions to the system of inequalities Ax ≤ b.

Since Ax ≤ b and A x ≤ b are just sets, we can write their combination as (Ax ≤ b) ∪ (A x ≤ b ). A special system of inequalities is a system of equations Dx = c, which is equivalent to the combined system of inequalities (Dx ≤ c) ∪ (-Dx ≤ -c). We say that a constraint system implies an inequality h T x ≤ g, where h ∈ Q n , h = 0 n , and g ∈ Q, if h T x ≤ g holds for all x ∈ Q(Ax ≤ b). In the same manner, a constraint system implies an equality h T x = g, where h ∈ Q n , h = 0 n , and

g ∈ Q, if h T x = g holds for all x ∈ Q(Ax ≤ b). A constraint implied by Ax ≤ b is explicit if it does appear in Ax ≤ b. Otherwise, it is called implicit.
Most deductions on linear inequalities are based on Farkas' Lemma:

Lemma 1 (Farkas' Lemma [START_REF] Boyd | Convex Optimization[END_REF]). Q(Ax ≤ b) = ∅ iff there exists a y ∈ Q m with y ≥ 0 m and y T A = 0 n so that y T b < 0, i.e., there exists a non-negative linear combination of inequalities in Ax ≤ b that results in an inequality y T Ax ≤ y T b that is constant and unsatisfiable. If such a y exists, then we call it a certificate of unsatisfiability.

We also frequently use the following lemma, which is just a reformulation of Farkas' Lemma:

Lemma 2 (Linear Implication Lemma). Let Q(Ax ≤ b) = ∅, h ∈ Q n \ {0 n }, and g ∈ Q. Then, Ax ≤ b implies h T x ≤ g iff
there exists a y ∈ Q m with y ≥ 0 m and y T A = h T so that y T b ≤ g, i.e., there exists a non-negative linear combination of inequalities in Ax ≤ b that results in the inequality h T x ≤ g.

As we mentioned in the introduction, this paper describes equisatisfiable transformations for constraint systems. We transform the systems in such a way that most linear mixed decision procedures become terminating and still retain their general efficiency. We even show this on the example of branchand-bound. Although we do not have the time to discuss all facets of branchand-bound [START_REF] Schrijver | Theory of Linear and Integer Programming[END_REF], we still want to give a short summary of the algorithm. Branch-and-bound is a recursive algorithm that computes mixed solutions for constraint systems. In each call of the algorithm, it first computes a rational solution s to a constraint system Ax ≤ b2 . If there are none, then we know that Ax ≤ b has no mixed solution. We are also done in the case that s is a mixed solution. Otherwise, we select one of the integer variables x i assigned to a fractional value s i ∈ Z and call branch-and-bound recursively on (Ax ≤ b) ∪ (x i ≥ s i ) and (Ax ≤ b) ∪ (x i ≤ s i ). If none of the recursive calls returns a mixed solution, then Ax ≤ b also does not have a mixed solution. Likewise, if one of them returns a mixed solution s, then it also is a mixed solution to Ax ≤ b.

Branch-and-bound alone is already complete on bounded constraint systems, i.e., systems where all directions are bounded:

Definition 1 (Bounded Direction). A direction/vector h ∈ Q n \ {0 n } is bounded in the constraint system Ax ≤ b if there exist l, u ∈ Q such that Ax ≤ b implies h T x ≤ u and -h T x ≤ -l. Otherwise, it is called unbounded. Definition 2 (Bounded System). A constraint system Ax ≤ b is bounded if all directions h ∈ Q n \{0 n } are bounded. Otherwise, it is called unbounded.
For bounded systems, branch-and-bound is one of the most popular and efficient algorithms. It may, however, diverge if the system has unbounded directions. Even so, not all unbounded systems are equally difficult. For instance, a system where all directions are unbounded has always a mixed solution:

Lemma 3 (Absolutely Unbounded [START_REF] Bromberger | Fast cube tests for LIA constraint solving[END_REF]). If all directions are unbounded in a constraint system Ax ≤ b, then the constraint system has an integer solution.

In a previous article, we described two cube tests that detect and solve constraint systems with infinite lattice width (another name for absolutely unbounded systems) in polynomial time [START_REF] Bromberger | Fast cube tests for LIA constraint solving[END_REF]. The case of absolutely unbounded systems is, therefore, trivial and branch-and-bound can be easily extended so it also becomes complete for absolutely unbounded systems.

The actual difficult case is when some directions are bounded and others unbounded. We call these systems partially unbounded. Here, branch-andbound and most other algorithms diverge or become inefficient in practice. The transformations, which we present, are designed to efficiently handle this subclass of problems.

Mixed-Echelon-Hermite Transformation

Our overall goal is to present an equisatisfiable transformation that turns any constraint system into a system that is bounded, i.e., a system on which branch-and-bound and many other arithmetic decision procedures terminate. In this section, we only present such a transformation for a subset of constraint systems, which we call double-bounded constraint systems. We then show in the next section that each constraint system can be reduced to an equisatisfiable double-bounded system. We also show how to efficiently transform a mixed solution from the double-bounded reduction to a mixed solution for the original system.

Definition 3 (Double-Bounded Constraint System). A constraint system Dx ≤ u is double-bounded if Dx ≤ u implies Dx ≥ l for l ∈ Q m .
For such a double-bounded system, we call the bounds u the upper bounds of Dx and the bounds l the lower bounds of Dx. Moreover, we typically write l ≤ Dx ≤ u instead of Dx ≤ u although the lower bounds l are only implicit.

Note that only the inequalities in a double-bounded constraint system are guaranteed to be bounded. Variables might still be unbounded. For instance, in the constraint system 1 ≤ 3x 1 -3x 2 ≤ 2 both inequalities are bounded but the variables x 1 and x 2 are not. Moreover, the above constraint system is also an example where branch-and-bound diverges. This means that even bounding all inequalities does not yet guarantee termination. So for our purposes, a double-bounded constraint system is still too complex.

This changes, however, if we also require that the coefficient matrix D of our constraint system is a lower triangular matrix with gaps: Definition 4 (Lower Triangular Matrix with Gaps). A matrix A ∈ Q m×n is lower triangular with gaps if it holds for each column j that piv(A, j) > m or that piv(A, j) < piv(A, k) for all columns k with j < k ≤ n, i.e., column j either has only zero entries or all pivoting entries right of j have a higher row index.

A matrix is lower triangular if and only if the row indices of its pivots are strictly increasing, i.e., piv(A, 1) < . . . < piv(A, n). If we also allow it to have gaps, only the row indices of pivots with non-zero columns have to be strictly increasing. Now we get termination for free because of our restrictions:

Lemma 4 (Lower Triangular Double-Bounded Systems). Let D ∈ Q m×n be a lower triangular matrix with gaps and l ≤ Dx ≤ u be a double-bounded constraint system. Then each variable x j is either bounded, i.e., l ≤ Dx ≤ u implies that l j ≤ x j ≤ u j or its column in D has only zero entries.

Proof. Proof by induction. Assume that the above property already holds for all variables x k with k < j. Let p = piv(D, j). If p > m, then the column j of D is zero and we are done. If p ≤ m, then the pivoting entry d pj of column j is non-zero. Because of Definition 4 and our induction hypothesis, this also means that each column k with k < j has either a zero entry in row p or the variable x k is bounded by our induction hypothesis, i.e., l ≤ Dx ≤ u implies l k ≤ x k ≤ u k . Since Definition 4 also implies that row p has only zero entries to the right of d pj , the row p has only one unbounded variable with a non-zero entry, viz., x j . This means we can transform the row l p ≤ d T p x ≤ u p into the following two inequalities:

l p -j-1 k=1 d pk x k ≤ d pj x j and u p -j-1 k=1 d pk x k ≥ d pj x j ,
where the variables x k on the left sides are either bounded or d pk = 0. Hence, we can derive an upper and lower bound for x j via bound propagation/refinement [START_REF] Jovanović | Cutting to the chase[END_REF].

Corollary 1 (BnB-LTDB-Termination). Branch-and-bound terminates on every double-bounded system l ≤ Dx ≤ u where D is lower triangular with gaps.

Our next goal is to efficiently transform every double-bounded system l ≤ Dx ≤ u into an equisatisfiable system that also has a lower triangular coefficient matrix with gaps. We start by defining a class of transformations that do not only preserve mixed equisatisfiability, but are also very expressive. [START_REF] Christ | Cutting the mix[END_REF]). Given a mixed constraint system. A matrix V ∈ Q n×n is a mixed column transformation matrix if it is invertible and consists of an invertible matrix

Definition 5 (Mixed Column Transformation Matrix

V (Q) ∈ Q n 1 ×n 1 , a unimodular matrix V (Z) ∈ Z n 2 ×n 2 , and a matrix V (M ) ∈ Q n 1 ×n 2 such that V = V (Q) V (M ) 0 n 2 ×n 1 V (Z)
.

The inverse of a mixed column transformation matrix V is also a mixed column transformation matrix and can be used to undo the transformation V :

Lemma 5 (Mixed Column Transformation Inversion [START_REF] Christ | Cutting the mix[END_REF]). Given a mixed constraint system. Let V ∈ Q n×n be a mixed column transformation matrix. Then V -1 is also a mixed column transformation matrix.

This means that each mixed column transformation matrix defines a bijection from (Q

n 1 × Z n 2 ) to (Q n 1 × Z n 2 )
. Hence, they guarantee mixed equisatisfiability:

Lemma 6 (Mixed Column Transformation Equisatisfiability [START_REF] Christ | Cutting the mix[END_REF]). Let Ax ≤ b be a mixed constraint system. Let V ∈ Q n×n be a mixed column transformation matrix. Then every solution y ∈ M((AV )y ≤ b)) can be converted into a solution V y = x ∈ M(Ax ≤ b) and vice versa.

Moreover, the mixed column transformation matrix V also establishes a direct relationship between the linear combinations of the original constraint system and the transformed one:

Lemma 7 (Mixed Column Transformation Implications). Let Ax ≤ b be a constraint system. Let V ∈ Q n×n be a mixed column transformation matrix. Let Ax ≤ b imply h T x ≤ g. Then AV z ≤ b implies h T V z ≤ g.
Proof. By Lemma 2, Ax ≤ b implies h T x ≤ g iff there exists a non-negative linear combination y ∈ Q n such that y ≥ 0, y T A = h T and y T b ≤ g. Multiplying y T A = h T with V results in y T AV = h T V and thus y is also the non-negative linear combination of inequalities AV z ≤ b that results in h T V z ≤ g.

Corollary 2 (Mixed Column Transformation Certificates). Let Ax ≤ b be a constraint system. Let V ∈ Q n×n be a mixed column transformation matrix. Then y is a certificate of unsatisfiability for Ax ≤ b iff it is one for AV z ≤ b. Now we only need a mixed column transformation matrix V for every coefficient matrix A such that H = AV is lower triangular with gaps. One such matrix V is the one that transforms A into Mixed-Echelon-Hermite normal form:

Definition 6 (Mixed-Echelon-Hermite Normal Form [12]). A matrix H ∈ Q m×n is in Mixed-Echelon-Hermite normal form if H = E 0 r×(n 1 -r) 0 r×n 2 E 0 (m-r)×(n 1 -r) H ,
where E is an r × r identity matrix (with r ≤ n 1 ), E ∈ Q (m-r)×r , and H ∈ Q (m-r)×n 2 is a matrix in hermite normal form, i.e., a lower triangular matrix without gaps, where each entry h piv(H ,j)k in the row piv(H , j) is non-negative and smaller than h piv(H ,j)j ).

The following proof for the existence of the Mixed-Echelon-Hermite normal form is constructive and presents the Mixed-Echelon-Hermite transformation.

Lemma 8 (Mixed-Echelon-Hermite Transformation). Let A ∈ Q m×n be a matrix, where the upper left r × n 1 submatrix has the same rank r as the complete left m × n 1 submatrix. Then there exists a mixed transformation matrix V ∈ Q n×n such that H = AV is in Mixed-Echelon-Hermite normal form.

Proof. Proof from [START_REF] Christ | Cutting the mix[END_REF] with slight modifications so it also works for singular matrices. We subdivide A into

A = A 11 A 12 A 21 A 22 such that A 11 ∈ Q r×n 1 , A 12 ∈ Q r×n 2 , A 21 ∈ Q m-r×n 1 , and A 21 ∈ Q m-r×n 2 .
Then we bring A 11 with an invertible matrix

V 11 ∈ Q n 1 ×n 1 into reduced echelon column form H 11 = (E 0 r×(n 1 -r) ) = A 11 V 11
, where E is an r × r identity matrix. We get V 11 and H 11 by using Bareiss algorithm instead of the better known Gaussian elimination as it is polynomial in time [START_REF] Bareiss | Sylvester's identity and multistep integer-preserving gaussian elimination[END_REF]. 3 Note that the last n 1 -r columns of H 21 = (H 21 0 (m-r)×(n 1 -r) ) = A 21 V 11 are also zero because all rows in A 21 are linear dependent of A 11 (due to the rank). Next we notice that

A 12 -A 11 V 11 A 12 0 (n 1 -r)×n 2 = A 12 -(E 0 r×(n 1 -r) ) A 12 0 (n 1 -r)×n 2
= 0 r×n 2 so we can reduce the upper right submatrix A 12 to zero by adding multiples of the n 1 columns with rational variables to the n 2 columns with integer variables. However, this also transforms the lower right submatrix A 22 into

H 22 = A 22 -A 21 V 11 A 12 0 (n 1 -r)×n 2 .
Finally, we transform this new submatrix H 22 into hermite normal form H 22 via the algorithm of Kannan and Bachem (or a similar polynomial time algorithm). 3 This algorithm also returns a unimodular matrix V 22 ∈ Z n 2 ×n 2 such that H 22 = H 22 V 22 . To summarize: our total mixed transformation matrix is

V =   V 11 -V 11 • A 12 0 (n 1 -r)×n 2 • V 22 0 n 2 ×n 1 V 22   and H = AV = H 11 0 r×n 2 H 21 H 22 .
It is not possible to transform every matrix A ∈ Q m×n into Mixed-Echelon-Hermite normal form. We have to restrict ourselves to matrices, where the upper left r × n 1 submatrix has the same rank r as the complete left m × n 1 submatrix. However, this is very easy to accomplish for a system of linear mixed arithmetic constraints l ≤ Ax ≤ u. The reason is that the order of inequalities does not change the set of satisfiable solutions. Hence, we can swap the inequalities and, thereby, the rows of A until its upper left r × n 1 submatrix has the desired form. This also means that there are usually multiple possible inequality orderings that each have their own Mixed-Echelon-Hermite normal form H.

To conclude this section: whenever we have a double-bounded constraint system l ≤ Dx ≤ u, we can transform it (after some row swapping) into an equisatisfiable system l ≤ Hy ≤ u where H = DV is in Mixed-Echelon-Hermite normal form and V y = x. Since H is also a lower triangular matrix with gaps, branch-and-bound terminates on l ≤ Hy ≤ u with a mixed solution t or it will return unsatisfiable (Corollary 1). Moreover, we can convert any mixed solution t for l ≤ Hy ≤ u into a mixed solution s for l ≤ Dx ≤ u by setting s := V t. Hence, we have a complete algorithm for double-bounded constraint systems.

Double-Bounded Reduction

In the previous Section, we have shown how to solve a double-bounded constraint system. Now we show how to reduce any constraint system A x ≤ b to an equisatisfiable double-bounded system l ≤ Dx ≤ u. Moreover, we explain how to take any solution of l ≤ Dx ≤ u and turn it into a solution for A x ≤ b .

As the first step of our reduction, we reformulate the constraint system into a so called split system:

Definition 7 (Split System). (Ax ≤ b) ∪ (l ≤ Dx ≤ u) is a split system if: (i) all directions are unbounded in Ax ≤ b; (ii) all row vectors a i from A are also unbounded in (Ax ≤ b) ∪ (l ≤ Dx ≤ u)
. Moreover, we call Ax ≤ b the unbounded part and l ≤ Dx ≤ u the bounded part of the split system.

A split system consists of an unbounded part Ax ≤ b that is guaranteed to have (infinitely many) integer solutions (see Lemma 3) and a doublebounded part l ≤ Dx ≤ u. Any constraint system can be brought into the above form. We just have to move all unbounded inequalities into the unbounded part and all bounded inequalities into the bounded part.

Lemma 9 (Split Equivalence). Let A x ≤ b be a constraint system with A ∈ Q m×n . Then there exists an equivalent split system (Ax ≤ b) ∪ (l ≤ Dx ≤ u) where: (i) A ∈ Q m 1 ×n and D ∈ Q m 2 ×n so that m 1 + m 2 = m; (ii) all rows d T i of D and a T k of A appear as rows in A ; and (iii) Dx ≤ u implies l ≤ Dx.

Proof. For (i), (ii), and the equivalence, it is enough to move all bounded inequalities a T i x ≤ b i of A x ≤ b into Dx ≤ u and all unbounded inequalities into Ax ≤ b. For (iii), we assume for a contradiction that Dx ≤ u does not imply l i ≤ d T i x but (Dx ≤ u)∪(Ax ≤ b) does. By Lemma 2, this means that there exists a y ∈ Q m 2 with y ≥ 0 m 2 and a z ∈ Q m 1 with z ≥ 0 m 1 so that y T D + z T A = -d T i and y T u + z T b ≤ -l i . We also know that there exists a z k > 0 because Dx ≤ u alone does not imply l i ≤ d T i x. We use this fact to reformulate y

T D+z T A = -d T i into -a T k = 1 z k y T D + d T i + m 1
j=1,j =k z j a T j , and use the bounds of the inequalities in Dx ≤ u and Ax ≤ b to derive a lower bound for a T k x:

-a T k x ≤ 1 z k y T u + u i + m 1 j=1,j =k z j b j .
Hence, a T k is bounded in A x ≤ b and we have our contradiction.

The above Lemma also shows that the bounded part of a constraint system is self-contained, i.e., a constraint system implies that a direction is bounded if and only if its bounded part does. The actual difficulty of reformulating a system into a split system is not the transformation per se, but finding out which inequalities are bounded or not. There are many ways to detect whether an inequality is bounded by a constraint system. Most work even in polynomial time. For instance, solving the linear rational optimization problem "minimize a T i x such that Ax ≤ b" returns -∞ if a i is unbounded, ∞ if Ax ≤ b has no rational solution, and the optimal lower bound l i for a T i x otherwise. However, it still requires us to solve m linear optimization problems.

A, in our opinion, more efficient alternative is based on our previously presented algorithm for finding equality bases [START_REF] Bromberger | Computing a complete basis for equalities implied by a system of LRA constraints[END_REF]. This is due to the following relationship between bounded directions and equalities: Lemma 10 (Bounds and Equalities).

Let Q(Ax ≤ b) = ∅. Then h is bounded in Ax ≤ b iff Ax ≤ 0 m implies that h T x = 0.
Proof. By Definition 1, h is bounded in Ax ≤ b means that there exists l, u ∈ Q such that Ax ≤ b implies h T x ≤ u and -h T x ≤ -l. By Lemma 2, this is equivalent to: there exist l, u ∈ Q, y, z ∈ Q m with y, z ≥ 0 m , and y T A = h T = -z T A so that y T b ≤ u and z T b ≤ -l. Symmetrically, Ax ≤ 0 implies that h T x = 0 is equivalent to: there exist a y, z ∈ Q m with y, z ≥ 0 m and y T A = h T = -z T A so that y T 0 m ≤ 0 and z T 0 m ≤ 0. Since u and l only have to exists, we can trivially choose them as u := y T b and l := -z T b. This means that y T b ≤ u, z T b ≤ -l, y T 0 m ≤ 0, and z T 0 m ≤ 0 are all trivially satisfied by any pair of linear combinations y, z ∈ Q m with y, z ≥ 0 m such that y T A = h T = -z T A. Hence, the two definitions are equivalent and our lemma holds.

It is easy and efficient to compute an equality basis for Ax ≤ 0 m and to determine with it the inequalities in Ax ≤ b that are bounded [START_REF] Bromberger | Computing a complete basis for equalities implied by a system of LRA constraints[END_REF]. The only disadvantage towards the optimization approach is that we do not derive an optimal lower bound l for the inequalities. This is no problem because only the existence of lower bounds is relevant and not the actual bound values.

In a split system (Ax ≤ b)∪(l ≤ Dx ≤ u), the unbounded part is actually inconsequential to the rational/mixed satisfiability of the system. It may reduce the number of rational/mixed solutions, but it never removes them all. Hence, (Ax

≤ b) ∪ (l ≤ Dx ≤ u) is equisatisfiable to just l ≤ Dx ≤ u.
We first show this equisatisfiability for the rational case:

Lemma 11 (Rational Extension). Let (Ax ≤ b) ∪ (l ≤ Dx ≤ u) be a split system. Let s ∈ Q n be a rational solution to the bounded part l ≤ Dx ≤ u such that Ds = g, where g ∈ Q m 2 . Then (Ax ≤ b) ∪ (Dx = g) has a solution s .
Proof. Assume for a contradiction that (Ax ≤ b) ∪ (Dx = g) has no solution. By Lemma 1, this means that there exist a y ∈ Q m 1 with y ≥ 0 m 1 and z, z ∈ Q m 2 with z, z ≥ 0 m 2 such that y T A + z T D -z T D = 0 n and y T b + z T g -z T g < 0. Since Dx = g is satisfiable by itself, there must exist a y i > 0. Now we use this fact to reformulate the equation

y T A + z T D -z T D = 0 n into -a T i = 1 y i m 1 j=1j =i y j a T j + z T D -z T D ,
from which we deduce a lower bound for

a T i x in (Ax ≤ b) ∪ (l ≤ Dx ≤ u): -a T i x ≤ 1 y i m 1 j=1j =i y j b j + z T u -z T l . Therefore, a i is bounded in (Ax ≤ b) ∪ (l ≤ Dx ≤ u), which is a contradic- tion.
Note that the bounded part l ≤ Dx ≤ u of a split system can still have unbounded directions (not inequalities). Some of these unbounded directions in l ≤ Dx ≤ u are the orthogonal directions to the row vectors d i , i.e., vectors v j ∈ Z n such that d T i v j = 0 for all i ∈ {1, . . . , m 2 }. This also means that the existence of one mixed solution s ∈ (Q n 1 × Z n 2 ) and one unbounded direction proves the existence of infinitely many mixed solutions. We just need to follow the orthogonal directions, i.e., for all λ ∈ Z, s = λ • v j + s is also a mixed solution because

d T i s = λ • d T i v j + d T i s = d T i s.
In the next two steps, we prove that Ax ≤ b cannot cut off all of these orthogonal solutions because it is completely unbounded. The first step proves that Ax ≤ b remains absolutely unbounded even if we settle on one set of orthogonal solutions, i.e., enforce Dx = Ds for some solution s.

Lemma 12 (Persistence of Unboundedness). Let (Ax ≤ b) ∪ (l ≤ Dx ≤ u) be a split system. Let s ∈ Q n be a rational solution for l ≤ Dx ≤ u such that Ds = g (with g ∈ Q m 2 ). Then all row vectors a i from A are still unbounded in (Ax ≤ b) ∪ (Dx = g).

Proof. By Lemma 11, (Ax ≤ b) ∪ (Dx = g) has at least a rational solution s * . Moreover, (Ax ≤ 0) ∪ (Dx = 0) does not imply a T i x = 0 because of Lemma 10 and the assumption that the row vectors a i from A are unbounded in (Ax ≤ b) ∪ (l ≤ Dx ≤ u). In reverse, (Ax ≤ b) ∪ (Dx = g) having a real solution, (Ax ≤ 0) ∪ (Dx = 0) does not imply a T i x = 0, and Lemma 10 prove together that the row vectors a i from A are also unbounded in (Ax ≤ b) ∪ (Dx = g).

The next step proves how to extend the mixed solution from the bounded part to the complete system with the help of the Mixed-Echelon-Hermite normal form and the absolute unboundedness of Ax ≤ b.

Lemma 13 (Mixed Extension

). Let (Ax ≤ b) ∪ (l ≤ Dx ≤ u) be a split system. Let s ∈ (Q n 1 × Z n 2 ) be a mixed solution for l ≤ Dx ≤ u. Then (Ax ≤ b) ∪ (l ≤ Dx ≤ u) has a mixed solution s .
Proof. Let g = Ds. Without loss of generality we assume that the upper left r × n 1 submatrix of D has the same rank r as the complete left m 1 × n 1 submatrix of D. (Otherwise, we just reorder the rows accordingly.) Therefore, there exists a mixed column transformation matrix V such that H = DV is in mixed-echelon-hermite normal form (see Lemma 8). By Lemma 6, there exists a mixed vector t ∈ (Q n 1 × Z n 2 ) such that s = V t and t is a mixedsolution to l ≤ Hy ≤ u as well as Hy = g. Let U be the set of indices with 0 columns in H and B the column indices with bounded variables. Then the equation system (Hy = g) fixes each variable y j with j ∈ B to the value t j because H is lower triangular with gaps. Hence, ((AV )y ≤ b) ∪ (Hy = g) is equivalent to

A    j∈U    v 1j . . . v nj    • y j    ≤ b -A    j∈B    v 1j . . . v nj    • t j    . (1)
Due to Lemma 12 and 7, all directions are unbounded in [START_REF] Althaus | Superposition modulo linear arithmetic SUP(LA)[END_REF]. This means (1) has an integer solution (Lemma 3) assigning each variable y j with j ∈ U to a t j ∈ Z. (Can be computed via the unit cube test [START_REF] Bromberger | New Techniques for Linear Arithmetic: Cubes and Equalities[END_REF]). We extend this solution to all variables y by setting t j := t j for j ∈ B and we have a mixed solution t ∈ (Q n 1 × Z n 1 ) for ((AV )y ≤ b) ∪ (l ≤ Hy ≤ u). Hence, we have via Lemma 6 a mixed solution s ∈ (Q

n 1 × Z n 2 ) for (Ax ≤ b) ∪ (l ≤ Dx ≤ u) with s = V t .
Corollary 3 (Double-Bounded Reduction). The split system (Ax ≤ b)∪(l ≤ Dx ≤ u) is mixed equisatisfiable to (l ≤ Dx ≤ u).

Experiments

We integrated the Double-Bounded reduction and the Mixed-Echelon-Hermite transformation into our own theory solver SPASS-IQ v0.2 4 and ran it on four families of newly constructed benchmarks 4 . Once with the transformations turned on (SPASS-IQ) and once with the transformations turned off (SPASS-IQ-Off ). If SPASS-IQ encounters a system Ax ≤ b that is not explicitly bounded, i.e., where not all variables have an explicit upper and lower bound, then it computes an equality basis for Ax ≤ 0 m . This basis is used to determine whether the system is implicitly bounded, absolutely unbounded or partially bounded, as well as which of the inequalities are bounded. Our solver only applies our two transformations if the problem is partially unbounded. The resulting equisatisfiable but bounded problem is then solved via branch-and-bound. The other two cases, absolutely unbounded and implicitly bounded, are solved respectively via the unit cube test [START_REF] Bromberger | New Techniques for Linear Arithmetic: Cubes and Equalities[END_REF] and branch-and-bound on the original system. Our solver also converts any mixed solutions from the transformed system into mixed solutions for the original system following the proof of Lemma 13. Rational conflicts are converted between the two systems by using Corollary 2. We tried to restrict our benchmarks to partially unbounded problems since we only apply our transformations on those problems. We even found some partially unbounded problems in the SMT-LIB benchmarks for QF LIA (quantifier free linear arithmetic). However, there are not many such bench- SlackedQFLIA: are linear integer benchmarks based on the SMT-LIB classes CAV-2009 [START_REF] Dillig | Cuts from proofs: A complete and practical technique for solving linear inequalities over integers[END_REF], cut lemmas [START_REF] Griggio | A practical approach to satisfiability modulo linear integer arithmetic[END_REF], and dillig [START_REF] Dillig | Cuts from proofs: A complete and practical technique for solving linear inequalities over integers[END_REF]. We simply took all of the unsatisfiable benchmarks and replaced in them all variables x with x + -x -where x + and x -are two new variables such that x + , x -≥ 0. This transformation, called slacking, is equisatisfiable and the slacked version of the dillig-benchmarks, called slacked [START_REF] Jovanović | Cutting to the chase[END_REF], is already in the SMT-LIB. Slacking turns any unsatisfiable problem into a partially unbounded one. Hence, all problems in SlackedQFLIA are partially unbounded. Slacking is commonly used to integrate absolute values into linear systems or for solvers that require non-negative variables [START_REF] Schrijver | Theory of Linear and Integer Programming[END_REF]. RandomUnbd : are linear integer benchmarks that are all partially unbounded and satisfiable with 10, 25, 50, 75, and 100 variables. All problems are randomly created via a sagemath script 4 .

FlippedQFLIA and FlippedRandomUnbd : are linear mixed benchmarks that are all partially unbounded. They are based on SlackedQFLIA and RandomUnbd. We constructed them by first copying ten versions of the integer benchmarks and then randomly flipping the type of some of the variables to rational (probability of 20%). Some of the flipped instances of SlackedQFLIA became satisfiable.

We compared our solver with some of the state-of-the-art SMT solvers currently available for linear mixed arithmetic: cvc4-1.5 [START_REF] Barrett | CVC4. In CAV[END_REF], mathsat5-5.1 [START_REF] Cimatti | The Math-SAT5 SMT Solver[END_REF], SMTInterpol 2.1-335-g4c543a5 [13], yices2.5.4 [START_REF] Dutertre | Yices 2.2[END_REF], and z3-4.6.0 [START_REF] De Moura | Z3: An efficient SMT solver[END_REF]. Most of these solvers employ a branch-and-bound approach with an underlying dual simplex solver [START_REF] Dutertre | A fast linear-arithmetic solver for DPLL(T)[END_REF], which is also the basis for our own solver. As far as we are aware, none of them employ any techniques that guarantee termination.

SMTInterpol extends branch-and-bound via the cuts from proofs approach, which uses the Mixed-Echelon-Hermite transformation to find more versatile branches and cuts [START_REF] Christ | Cutting the mix[END_REF]. Although the procedure is not complete, the similarities to our own approach make an interesting comparison. Actually, the Double-Bounded reduction alone would be sufficient to make SMTInterpol terminating since it already builds branches via a Mixed-Echelon-Hermite transformation.

We also compared our solver with the ctrl-ergo solver [START_REF] Bobot | A simplex-based extension of fouriermotzkin for solving linear integer arithmetic[END_REF] although it is restricted to pure integer arithmetic. Ctrl-ergo is complete over linear integer arithmetic and uses the most similar approach to our transformations that we found in the literature. It dynamically eliminates one linear independent bounded direction at a time via transformation. The disadvantages of the dynamic approach are that it is very restrictive and does not leave enough freedom to change strategies or to add complementing techniques. Moreover, ctrl-ergo uses this transformation approach for all problems and not only the partially unbounded ones, which sometimes leads to a massive overhead on bounded problems.

For the experiments, we used a Debian Linux cluster and allotted to each problem and solver combination 2 cores of an Intel Xeon E5620 (2.4 GHz) processor, 4 GB RAM, and 40 minutes. The only solver benefiting from multiple cores is SMTInterpol. The plots in Figures 1 and2 depict the results of the different solvers. In the legends of the plots, the numbers behind the solver names are the number of solved instances. For FlippedQFLIA, there are two numbers to indicate the number of satisfiable/unsatisfiable instances solved. This is only necessary for FlippedQFLIA because it is the only tested benchmark family with satisfiable and unsatisfiable instances. (We verified that the results match if two solvers solved the same problem.)

Although our solver could not solve all problems (due to time and memory limits) it was still able to solve more problems than the other solvers. It was also faster on most instances than the other solvers. In some of the unsatisfiable, partially unbounded benchmarks ctrl-ergo is better than SPASS-IQ. This is due to its conflict focused, dynamic approach. For the same reason, ctrl-ergo is slower on the satisfiable, partially unbounded benchmarks. Only SPASS-IQ, ctrl-ergo, and yices solved all of the ten original SMT-LIB benchmarks that are partially unbounded, though the complete methods were still a lot faster (SPASS-IQ took 23s, ctrl-ergo took 42s, and yices took 1273s). On one of these benchmarks, 20-14.slacks.smt2 from slacks, all other solvers seem to diverge. Another interesting result of our experiments is that relaxing some integer variables to rational variables seems to make the problems harder instead of easier. We expected this for our transformations because the resulting systems become more complex and less sparse, but it is also true for the other solvers. The reason might be that bound refinement, a technique used in most branch-and-bound implementations, is less effective on mixed problems.

The time SPASS-IQ needs to detect the bounded inequalities and to apply our transformations is negligible. This is even true for the implicitly bounded problems we tested. As mentioned before, we do not have to apply our transformations to terminate on bounded problems. This is also the only advantage we gain from detecting that a problem is implicitly bounded. Since there is no noticeable difference in the run time, we do not further elaborate the results on bounded problems, e.g. with graphs.

An actual disadvantage of our approach is that the Mixed-Echelon-Hermite transformation increases the density of the coefficient matrix as well as the absolute size of the coefficients. Both are important factors for the efficiency of the underlying simplex solver. Moreover, SPASS-IQ reaches more often the memory limit than the time limit because it needs a (too) large number of branches and bound refinements before terminating.

Conclusion

We have presented the Mixed-Echelon-Hermite transformation (Lemma 8) and the Double-Bounded reduction (Lemma 9 & Corollary 3). We have shown that both transformations together turn any constraint system into an equisatisfiable system that is also bounded (Lemma 4). This is sufficient to make branch-and-bound, and many other linear mixed decision procedures, complete and terminating. We have also shown how to convert certificates of (un)satisfiability efficiently between the transformed and original systems (Corollary 2 & Lemma 13). Moreover, experimental results on partially unbounded benchmarks show that our approach is also efficient in practice.

Our approach can be nicely combined with the extensive branch-andbound framework and its many extensions, where other complete techniques cannot be used in a modular way [START_REF] Bobot | A simplex-based extension of fouriermotzkin for solving linear integer arithmetic[END_REF][START_REF] Bromberger | Linear integer arithmetic revisited[END_REF]. For future research, we plan to test our transformations in combination with other algorithms, e.g., cuts from proofs, or as a dynamic version similar to the approach used by ctrl-ergo [START_REF] Bobot | A simplex-based extension of fouriermotzkin for solving linear integer arithmetic[END_REF]. We also want to test whether our transformations are useful preprocessing

Figure 1 :

 1 Figure 1: horizontal axis: # of solved instances; vertical axis: time (seconds)

Figure 2 :

 2 Figure 2: horizontal axis: # of solved instances; vertical axis: time (seconds)

All techniques discussed in this paper can be extended to strict inequalities with the help of δ-rationals[START_REF] Dutertre | A fast linear-arithmetic solver for DPLL(T)[END_REF]. We will omit the strict inequalities and focus only on non-strict inequalities due to lack of space.

A rational solution can be computed in polynomial time[START_REF] Schrijver | Theory of Linear and Integer Programming[END_REF].

We do actually use less efficient, Gaussian-elimination-based transformations in our own implementation[START_REF] Bromberger | A reduction from unbounded linear mixed arithmetic problems into bounded problems[END_REF]. The reason is that these transformations are incrementally efficient. Our experiments show that the transformation cost still remains negligible in practice.

Available on http://www.spass-prover.org/spass-iq

steps for select constraint system classes that are bounded.


Boris P Hejblum
email: bhejblum@hsph.harvard.edu

Griffin M Weber

Katherine P Liao

Nathan P Palmer

Susanne Churchill

Nancy A Shadick

Peter Szolovits

Shawn N Murphy

Isaac S Kohane

Tianxi Cai

Probabilistic record linkage of de-identified research datasets with discrepancies using diagnosis codes Authors

We develop an algorithm for probabilistic linkage of de-identified research datasets at the patient level, when only diagnosis codes with discrepancies and no personal health identifiers such as name or date of birth are available. It relies on Bayesian modelling of binarized diagnosis codes, and provides a posterior probability of matching for each patient pair, while considering all the data at once. Both in our simulation study (using an administrative claims dataset for data generation) and in two real use-cases linking patient electronic health records from a large tertiary care network, our method exhibits good performance and compares favourably to the standard baseline Fellegi-Sunter algorithm. We propose a scalable, fast and efficient open-source implementation in the ludic R package available on CRAN, which also includes the anonymized diagnosis code data from our real use-case. This work suggests it is possible to link de-identified research databases stripped of any personal health identifiers using only diagnosis codes, provided sufficient information is shared between the data sources.

Introduction

An increasing quantity and variety of health data, including administrative claims data, electronic health records (EHR) data, and data generated from biomedical research studies, are becoming available for discovery research. For a subset of patients who contribute information to multiple datasets, synthesizing information from these multiple datasets may provide a more complete picture of the patients' health condition and enable more comprehensive population studies. For example, if research database A contains genomic data and a small set of phenotypic data and claims database B contains a complete list of diagnostic codes associated with patients' medical histories, linking databases A and B would enable researchers to investigate the role of genomics in a wide range of phenotypic conditions. More broadly, analyses of linked biomedical data could potentially lead to better understanding of disease risks and treatment effects 1,2 . Repurposing existing datasets through record linkage can reduce the cost of research through better exploitation of costly data already collected or curated. Concern for patients' privacy often makes such research datasets available only in deidentified forms, thus preventing linkage by use of a universal identifier or another nearlyunique personal characteristic of each patient. Our goal here is to investigate the use of informative characteristics of the de-identified data themselves that permit probabilistic linkage among different datasets featuring discrepancies (a discrepancy designates a disparity between two data sources concerning the same information).

Existing record linkage methods can be classified into two broad categories: i) deterministic approaches, where a fixed set of rules determine which records are matching and which are not; and ii) probabilistic methods, where discrepancies can be accounted for and a linkage probability is estimated for each record pair. Many factors can contribute to data discrepancy, such as missing information when a patient seeks care at multiple institutions and not all records are stored in one EHR, or administrative recording errors. Deterministic approaches are highly susceptible to these discrepancies, which has pushed the development of probabilistic record linkage methods [3][4][START_REF] Schmidlin | Privacy Preserving Probabilistic Record Linkage (P3RL): a novel method for linking existing health-related data and maintaining participant confidentiality[END_REF][START_REF] Sayers | Probabilistic record linkage[END_REF] , for instance to link cohort studies with registry data [START_REF] Moore | Poor record linkage sensitivity biased outcomes in a linked cohort analysis[END_REF][START_REF] Sengayi | Record linkage to correct under-ascertainment of cancers in HIV cohorts: The Sinikithemba HIV clinic linkage project[END_REF] . A statistical framework for probabilistic record linkage was originally proposed by Fellegi & Sunter [START_REF] Fellegi | A Theory for Record Linkage[END_REF] , formalizing Newcombe's ideas [START_REF] Newcombe | Automatic Linkage of Vital Records[END_REF] . Since then, several improvements have been proposed [START_REF] Winkler | Using the EM Algorithm for Weight Computation in the Fellegi-Sunter Model of Record Linkage[END_REF][START_REF] Winkler | Frequency-Based Matching in the Fellegi-Sunter Model of Record Linkage[END_REF][START_REF] Winkler | Improved Decision Rules In The Fellegi-Sunter Model Of Record Linkage[END_REF][START_REF] Larsen | Iterative Automated Record Linkage Using Models Mixture[END_REF][START_REF] Grannis | Analysis of a probabilistic record linkage technique without human review[END_REF][START_REF] Ravikumar | A Hierarchical Graphical Model for Record Linkage[END_REF][START_REF] Bhattacharya | A Latent Dirichlet Model for Unsupervised Entity Resolution[END_REF][START_REF] Murray | Probabilistic Record Linkage and Deduplication after Indexing, Blocking, and Filtering[END_REF] , relying on EM-based latent-class algorithms. Other approaches use machine learning methods, but they systematically require representative training data for which the true matching status is known, and this "gold standard" is typically not available [START_REF]Methodological Developments in Data Linkage[END_REF] . Recently, privacy-preserving record linkage (PPRL) methods relying on encryption of identifiable patient data have also been developed [START_REF] Trepetin | Privacy-Preserving String Comparisons in Record Linkage Systems: A Review[END_REF][START_REF] Kum | Privacy preserving interactive record linkage (PPIRL)[END_REF][START_REF] Kho | Design and implementation of a privacy preserving electronic health record linkage tool in Chicago[END_REF] .

Regardless of the algorithm, most approaches to patient record linkage rely on individually identifiable health information, such as name, date of birth, and social security number, which the Health Insurance Portability and Accountability Act (HIPAA) classifies as protected health information (PHI). In contrast, datasets generated for research purposes are often deidentified, where all PHI has intentionally been removed to protect patient privacy. The objective of our study is to determine whether it is possible to link different de-identified clinical research datasets using informative characteristics of the clinical research data such as patients' diagnoses (although linking records in two de-identified datasets does not actually identify any individual patients, caution must be taken since it potentially increases the risk that larger patterns unique to a specific individual may help an ill-intentioned adversary correlate these with other identifying information).

To our knowledge, Bennett et al. 4 published the only attempt to date at doing this. However, they assume that there are no discrepancies between the different data sources, and they never address the choice of priors in their modelling (e.g., on the number of matches, or any other priors). In addition, they require a subsample for which true matches are known, in order to tune certain parameters of their model. Finally, the details of their modelling and of their algorithm are not described in their publication, and the software they use is only available through purchase.

We address these limitations in this study, specifically testing whether it is possible to link different de-identified clinical research datasets using only a patient's diagnosis codes when discrepancies are present. Many types of data other than diagnoses can also exist in deidentified datasets, such as medications and procedures, and these might even be better than diagnoses for record linkage. However, our focus here is just on diagnoses because of how commonly they appear in research datasets, and because of the widespread use of standard vocabularies (e.g., ICD-9 and ICD-10). Furthermore, diagnosis codes have been shown to uniquely identify most patients within a single dataset [START_REF] Loukides | The disclosure of diagnosis codes can breach research participants' privacy[END_REF] . We ignore dates because they are considered PHI, and a de-identified dataset would have no more than the year of an event.

In this study, we evaluate the conditions under which diagnoses contain enough concordant information across datasets to enable linkage. This is important because contrary to PHI and other demographics, which are expected to be the same across datasets (assuming there are no recording errors), diagnosis codes often have valid differences. For example, a patient might be receiving ongoing care for diabetes at one hospital and visit a separate urgent care clinic closer to home for a common cold; or, a healthcare facility might use some refined diagnosis codes internally that are not necessarily sent out for billing, and thus will not always be present in the associated claims data. This means it is important to account for and to model such discrepancies, and that in practice our approach will only be successful in specific situations where the overlapping diagnosis information between the datasets is sufficient. Such situations include linking registry data -previously linked with billing data from a health system source using PHI -with further EHR data, or linking EHR data with a claims database.

Here we propose a probabilistic record linkage method that only uses diagnosis codes to estimate a posterior probability of matching between records. We rely on a likelihood framework using a Bayesian approach and a generative model to formally estimate the posterior probability of each patient record pair being a match.

Results

Numerical study

We first performed simulation studies based on the 2010 ICD9 code data from 3,153 patients in an un-identifiable claims database from a nationwide US health insurance plan. A total of 2,850 unique ICD9 codes were recorded for these patients with an average of 4.7 codes per patient. To create a second dataset we perturbed the original data with correlated multivariate normal noise as follows: ! * = $ %⋅'(),+,).()/%)⋅'(/0,+,) where Σ is the empirical correlation matrix of the 2,850 codes, X the original binary data, X * the second dataset, 4 the level of perturbation (noise) introduced and $ 5 the indicator function that is 1 if 7 > 0 and 0 otherwise. Then a subset of 2,500 patients from the original data was matched against a subset of 1,253 patients from the perturbed data, with only partial overlap (not all patients have a true match). The proportion of patient overlap varied between different simulation settings. Because the datasets are derived from the same source, the true matches were known.

We compare our method to the Fellegi-Sunter approach, estimated via an E-M algorithm from Grannis et al. [START_REF] Grannis | Analysis of a probabilistic record linkage technique without human review[END_REF] , itself based on Winkler's method [START_REF] Winkler | Using the EM Algorithm for Weight Computation in the Fellegi-Sunter Model of Record Linkage[END_REF] . Such an approach is still widely used and considered as a competitive unsupervised probabilistic approach for record linkage [START_REF]Methodological Developments in Data Linkage[END_REF] , even though it is usually used for matching records on a small set of PHI and has never been evaluated for matching clinical data based only on a large number of diagnosis codes. The Fellegi-Sunter method does not support 1-1 matching (the constraint that a given patient record should be matched at most once), yielding exceedingly low positive predictive value (PPV) due to a large number of potential matches for each patient. To reduce this number of false matches in the Fellegi-Sunter method, we also implemented a variant where only the largest matching score is considered as a potential match for each patient (in case of ties at the maximum matching score values, both matches were kept). This enforcement of the 1-1 matching improves statistical performance compared to the original Fellegi-Sunter method. Fig. 1 displays both the Positive Predictive Value (PPV) and the True Positive Rate (TPR) under various simulation settings. As expected, the 0.9 cutoff has higher PPV, while the 0.5 cutoff has higher TPR under every circumstance. For both cutoffs, our method exhibits good PPV for most settings and a TPR > 0.5 when (a) the noise perturbation level is less than 1.0 and (b) at least 1000 codes with low prevalence were used for matching (Fig. 1a and1b). On the contrary, the Fellegi-Sunter approach, even after the 1-1 matching modification, performs much worse with a PPV always below 20% and a TPR not much higher. This poor performance of the Fellegi-Sunter method is partly due to the fact that it does not distinguish between a (0,0) concordance and a (1,1) concordance, resulting in a high matching score for pairs with a high number of matching 0's. In addition, the Fellegi-Sunter method also does not account for discrepancies between the data sources. From Fig. 1b we also see that the more codes that are included in the matching procedure, the better the performance. In addition, although performance is acceptable with only a small overlap proportion between the two datasets, there is a noticeable improvement of the matching as the overlap proportion increases (Fig. 1c). Value (PPV). F-S refers to the Fellegi-Sunter method while ludic denotes our proposed Bayesian approach.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • 0.00 0.
• ludic -0.5 cutoff ludic -0.9 cutoff F-S F-S 1-1 • • • • • • • • • • • • • • • • • • • • • • • • 0.00 0.
• • 0.5 cutoff 0.9 cutoff F-S F-S 1-1 • • • • • • • • • • • • • • • • 0.00 0.25 0.50 0.75 1.00 0.1 0.5 1 2
Noise perturbation level introduced in the simulation Criterion TPR PPV Matching method In our simulations, we observed that the matching results are quite insensitive to the user input value for the hyper-parameter : . (representing the vector of marginal discrepancy rates between records, see Methods) provided that the simulated datasets shared enough diagnosis information, as shown in Fig. 2. A larger value of : . can slightly improve the performance if the true discrepancy is high, although a small fixed value of : . tends to work well. Concerning : / we have seen that it is generally better to provide a low value -this is expected in a context where diagnoses have low prevalence, making it unlikely to observe a discrepancy for a diagnosis given that it has been recorded in one of the datasets.

• • 0.5 cutoff 0.9 cutoff F-S F-S 1-1 • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0.
Note that our linkage algorithm requires that the two datasets contain overlapping diagnosis information (otherwise linkage is not feasible). For this reason, our approach will likely fail to link most de-identified datasets (which is the point of de-identifying a dataset). Yet it should be able to succeed in specific situations where diagnosis information overlap is sufficient (e.g., linking EHR data and claims data, or linking data collected over overlapping time span). In the next section, we apply our linkage algorithm to two sets of de-identified datasets on patients with rheumatoid arthritis (RA) with some matching patients, containing diagnosis codes recorded over overlapping time periods. . which models the probability of discrepancy between the two datasets. The top panel displays results when true value of ; . is 3%, while the bottom panel is for a true value of 20%.

Real world use-cases

We first applied our linkage algorithm to two de-identified datasets containing electronic medical records (EMR) from Partners HealthCare (a large hospital and physician network) in Boston, Massachusetts. Both EMR datasets were generated as part of prior research studies to look at patients with RA. The RA1 dataset contained 26,681 patients who had at least one Ideally, the performance of the algorithm would be evaluated compared to a gold standard where the match between subjects in RA1 and RA2 are known. However, the true matching status for subjects in the two cohorts was not known. Fortunately, the full RA1 and RA2 datasets have laboratory test results and dates, which we used to create a silver standard as a proxy for true matching status in order to benchmark our method (the full RA1 and RA2 datasets lack demographic information, but they are not strictly de-identified since they contain dates). Subjects in RA1 with the same laboratory test results and dates in the RA2 dataset were assumed to be a true match, because it would be very unlikely that two different patients would have the exact same test result at the exact same date. This silver standard identified 3,831 subjects who are likely the same person in RA1 and RA2 according to this silver standard. Since not all patients have laboratory test results and the databases were created during different years, the silver standard is likely missing some true matches. Nevertheless, these silver standard labels should be sufficiently accurate to form a basis for evaluating the performance of our algorithm compared to the Fellegi-Sunter approach. We then removed the laboratory test results and dates from RA1 and RA2 to generate truly de-identified datasets where only diagnosis codes are available for linkage.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • true ε + =0
The two datasets show very few discrepancies (i.e., when for a true match a diagnosis is recorded in one database and not in the other) when only 6 years of records are considered, and exhibit similar average statistics in Table 1. This is expected, because most patients in the RA2 set are expected to be a subset of the RA1 set, although we still anticipate some discrepancy between the two datasets (even with the same 6-year time span) since the two datasets were frozen at different times and records got updated over time due to administrative reasons. When considering the whole 11 years of records available, the extra 5 years allowed the second dataset to record additional diagnosis codes not present in the first one, thus introducing more discrepancies between the two. Such discrepancies can be seen in the average number of diagnoses per patient among the silver standard matches, for instance (see Table 1). Our method exhibits very good performance with both TPR and PPV always above 80%, even when considering the 11-year time span, as can be seen from Table 2. Due to computational limitations of the F-S method, we had to implement a blocking strategy using age to compute the F-S results. As in the numerical study, the Fellegi-Sunter method shows extremely poor performance without the 1-1 matching enforcement modification. And even with it, the PPV does not exceed 26%, even dropping to a low of 17% when considering the 11-year time span. Similarly, and to further illustrate the relevance of our proposed method, we have also linked a subsample of the RA1 dataset with data from the BRASS registry cohort [START_REF] Iannaccone | Using genetic and clinical data to understand response to disease-modifying anti-rheumatic drug therapy: data from the Brigham and Women's Hospital Rheumatoid Arthritis Sequential Study[END_REF] . The RA1 dataset contains electronic medical record data, collected as part of routine care, as well as a Genetic Risk Score (GRS) for RA on 1,454 subjects, calculated as part of a research study [START_REF] Chibnik | Genetic Risk Score Predicting Risk of Rheumatoid Arthritis Phenotypes and Age of Symptom Onset[END_REF] . The BRASS cohort is a prospective observational study that started in 2003 to enroll patients diagnosed with RA at the Robert Breck Brigham Arthritis Clinic in Boston, MA (USA). Notably it recorded longitudinal measurements of RA Disease Activity Score (DAS) for 1,130 patients, a composite score consisting of patient and physician reported information. DAS, a validated measure of RA disease activity used in clinical trials, is not routinely measured as part of clinical care. Linking those two sub-cohorts would enable studying the association between this GRS and the DAS in RA patients, an analysis that could not be performed in either of the datasets alone. There were 4,126 diagnosis codes recorded at least once in either database between 2002 and 2008 that were used for the linkage. Setting both ; . and ; / to 0.01, and with a threshold at 0.5 (respectively 0.9) on the matching posterior probability, this linkage exhibited a PPV of 0.92 (resp. 0.98) and a TPR of 0.76 (resp. 0.71). Here, performance indicators were again computed against 442 silver standard matches as a proxy (created from date-time stamps on the diagnosis and lab codes).

Discussion

In this work, we developed a method that can effectively use diagnosis codes to link two deidentified research datasets. Our approach presents three innovations compared to available state-of-the-art probabilistic linkage algorithms: i) our method uses only diagnosis codes, and no unique or direct identifiers are required; ii) our method takes into account the expected discrepancies across datasets in modeling; iii) our method does not need any training subsample for tuning, nor does it need any human intervention. In addition, our method is scalable to the use of thousands of diagnosis codes at once. Traditional approaches for probabilistic record linkage, based on Fellegi-Sunter seminal method were designed for using only a few informative PHI variables to identify matching patients across datasets. However, such methods are not equipped to deal with the high dimensionality of diagnosis codes and their sparse information content (see Supplementary Information for details).

Our method takes advantage of the likelihood framework to give more weight either to the presence or absence of a diagnosis. This is particularly useful for diagnosis codes as the presence of a diagnosis is usually much more informative than its absence. However, in case a particular code would be quite common and its absence would actually be more informative, the use of empirical Bayes priors for the marginal prevalence of codes takes this into account and adaptively tunes those weights according to the actual prevalence of each code in the data. Thus, our method naturally accounts for each diagnosis frequency. Because extremely rare codes might dominate other information too much, it can be useful to use a pre-filtering to discard such extremely rare codes. In addition, our method uses the constraint that a given patient record should be matched at most once, leveraging the information from all the available patient records when estimating the posterior probability of a match for a given pair.

In both our simulation study and our real use case, the Fellegi-Sunter method performs quite poorly. This can be explained because it was not designed to deal with either high-dimensional or sparse-binary variables, but rather with a small number of highly informative PHI variables. One way to improve the results from Fellegi-Sunter could be to compute agreement between two records differently, for instance treating a concurring diagnosis absence as noninformative, even though this could turn out to be harmful for highly prevalent diagnoses. On the contrary, our method tackles this issue by automatically tuning itself to the prevalence of each diagnosis.

Record linkage always requires some overlapping of information among the true matches between the datasets. If two datasets contain completely different information for the true matches, linkage will fail -regardless of the method used. Concerning diagnosis codes, there could be various reasons why information across databases would not overlap, such as disjoint time-periods for the record collection, or records from different healthcare systems. With less overlapping information, matches are all the less likely to be recovered by our approach. For instance, in the extreme case of no overlapping years, a few true matches might still be recovered for patients who have rare chronic conditions, but the sensitivity of the linkage would be quite poor. In practice, such degradation of shared information, across time or space, would give mediocre linkage performance indicated by low posterior matching probabilities.

EHRs can sometimes contain duplicate patients' records. In this work, we assume that each record is unique in one dataset, which we take advantage of by enforcing 1-1 matching. In practical cases, one might be required to perform a deduplication step on each dataset before linking them. Deduplication itself can be viewed as a matching problem and our approach can be easily adapted to solve it, for instance by linking a dataset to itself. Another straightforward extension of the proposed method is the inclusion of demographic blocking if available and needed. This can be done by adding a likelihood term comparing those demographics and assigning a very high prior weight for disagreement. Also, one limitation of the method proposed here is that it uses only binarized diagnosis codes. It ignores the actual number of times clinicians have assigned the same diagnosis code to a patient through multiple visits. A natural extension of our method would be to derive the likelihood for actual count-data instead of binary distributed codes, in order to use more information. This could be done for instance by using bivariate zero-inflated distributions 28 in the likelihood. A final limitation is the potential risk to patients from re-identifying a de-identified dataset by linking it to a dataset containing PHI. Often the data use agreements for research datasets forbid any attempts at doing record linkage for this reason, but the use of formal privacy models could help mitigate that risk [START_REF] Gkoulalas-Divanis | Publishing data from electronic health records while preserving privacy: A survey of algorithms[END_REF][START_REF] Poulis | Anonymizing datasets with demographics and diagnosis codes in the presence of utility constraints[END_REF] .

Although we expect the proposed patient linkage algorithm for de-identified research data using diagnosis codes to be broadly applicable, the good performance of the algorithm in real applications thus far has only been demonstrated in linking RA cohorts. Additional empirical assessment of its linkage performance in a broader set of real applications is warranted to ascertain its general applicability.

Various factors can influence the linkage performance, such as the record time-periods, the record health care systems, the prevalence of diagnoses, quality and accuracy of the records, health care utilization, or the proportion of patient overlap. However, based on both realistic simulations and two real world use-cases, we show that it is possible to achieve good linkage performance using only diagnosis codes as long as i) there is sufficient diagnosis information shared between the two datasets, and ii) there are enough rare diagnoses. We provide an open-source implementation of the proposed method, with complete documentation, in the R package ludic (Linkage Using DIagnosis Codes) available on CRAN. In addition, the ludic package also includes two anonymized diagnosis code datasets from our real use-case (along with a silver standard of true matches) for benchmarking future record linkage methods.

Methods

Data sources

We used three data sources for this study. The first was a de-identified administrative claims dataset from a U.S. based nationwide health insurer, which included all billing diagnosis codes recorded in 2010 for women between 50 and 69 years old with zip codes 021XX. Procedure codes are not included here, but their inclusion would be straightforward in the proposed approach. This dataset provided us a basis for generating diagnosis codes with realistic frequency distributions in our simulation studies. The second data source is from the EHR of patients with diagnostic codes of rheumatoid arthritis (RA) at Partners Healthcare in Boston, Massachusetts. The third data source is from the BRASS RA registry cohort which enrolls patients diagnosed with RA at the Robert Breck Brigham Arthritis Clinic in Boston, MA (USA).

We used the second and third data sources for validating the performance of our algorithm in real world settings. We first tested the algorithm's performance using two different deidentified EHR datasets of patients with diagnostic codes of RA. The two datasets contain overlapping but not identical subjects with RA identified via an EHR phenotyping algorithm [START_REF] Liao | Electronic medical records for discovery research in rheumatoid arthritis[END_REF][START_REF] Liao | Methods to Develop an Electronic Medical Record Phenotype Algorithm to Compare the Risk of Coronary Artery Disease across 3 Chronic Disease Cohorts[END_REF] in 2009 for the first dataset (RA1), and then again in 2014 for the second dataset (RA2). For each cohort, the structured data, including all diagnosis codes, were previously extracted from the EHR and stored as a de-identified dataset. The original datasets were actually limited datasets, containing not only diagnoses, but also the exact dates and results of laboratory tests. We first used these additional variables to deterministically link the datasets, based on the assumption that the dates and test results would be unique enough in these relatively small cohorts to perform accurate deterministic linkage. Because we cannot guarantee that this linkage is perfectly correct, we call this a "silver standard" rather than gold standard mapping. We then removed the dates and test results to create two de-identified RA datasets, which we could use to test our probabilistic record linkage algorithm using only diagnosis codes and compare our results to this silver standard mapping. We further validated our algorithm by linking the RA1 dataset with the BRASS cohort which contains longitudinal disease activity scores, inflammation markers as well as EHR data on diagnostic codes.

The RA studies are examples of real-life use cases demonstrating the need for a probabilistic algorithm. The RA1 cohort contains structured EHR data, single nucleotide polymorphism (SNP) data, auto-antibody marker measurements, as well as abstracted information from clinical notes obtained using chart review. In the RA2 cohort, additional data were abstracted from clinical notes, and biomarkers as well as genetic methylation levels were measured.

While the genetic and biomarker data are non-overlapping between the two cohorts, there is partial overlap for the structured EHR data, namely diagnosis codes among subjects in both cohorts. Being able to link the datasets would enable combined analysis of the two complementary sets of genetic and biomarker data between the two cohorts. Linking RA1 with the BRASS cohort would enable association studies on the effect of auto-antibody biomarkers on the longitudinal disease activity or inflammation markers in RA patients.

The Institutional Review Board at Harvard University approved the use of the administrative claims data. The Institutional Review Board at Partner's Healthcare Systems approved the use of the two RA datasets for probabilistic record linkage.

Linkage algorithm

Our goal is to link dataset < consisting of = > records to dataset ? consisting of = @ records. Assume that A common features, such as diagnosis codes, are recorded in both datasets. For B = 1, … , = > and D = 1, … , = @ , let E FG be the binary indicator of whether these two records are a true match, < (F) = (H ) (F) , … , H I (F) ) represent the A features for patient B from <, and

? (G) = (J ) (G) , … , J I (G)
) represent the A features for patient D from ?. All K features are assumed to be binary, e.g., presence or absence of a diagnosis code. A key component in our linkage algorithm is the posterior probability of being a match for each possible pair {B, D} of patients from < and ?: {M FG = ℙOE FG = 1P<, ?Q: B = 1, … , = > , D = 1, … , = @ }. We next detail key components of the algorithm. See Fig. 4 for the complete workflow of our algorithm.

Bayesian model

To compute the posterior matching probability M FG , we first compute a similarity measure for < (F) and ? (G) through a likelihood formulation. To this end, we assume that the respective distributions of < (F) and ? (G) do not depend on whether they can be matched with M > T = U(H V (F) = 1) and M @ T = U(J V (G) = 1). To allow discrepancy between the records in the two datasets, let

; V/ = ℙ(J V (G) = 0 | H V (F) = 1, E FG = 1) and ; V. = ℙ(J V (G) = 1 | H V (F) = 0, E FG =
1) as the respective marginal discrepancy rates for the X YZ feature between records. These hyper-parameters can be related to the [ V and \ V parameters from the Fellegi-Sunter approach by the following equations:

[ V = M > T (1 -; V/ ) + (1 -M > T )(1 -; V. ) and \ V = M > T M @ T + (1 -M > T )(1 -M @ T ).
The log-likelihood ratio for observing H V (F) and J V (G) is ℒ V `HV (F) , J V (G) a, where

ℒ V (b, c) = log g ℙ(H V (F) = b, J V (G) = c | E (FG) = 1) ℙ(H V (F) = b, J V (G) = c | E (FG) = 0) h = $ ij),kj) log l ; V/ M @ T m +$ ijn,kjn log l ; V. M @ T m + $ ij),kjn log l ; V/ M @ T m + $ ijn,kj) log l ; V. M @ T m
where for brevity we write ; = 1 -; for any ;. Then we define the similarity measure between < (F) and ? (G) as ℒ (FG) = ∑ Vj) I ℒ V `HV (F) , J V (G) a which essentially is the log-likelihood under a naive Bayes assumption, that is commonly used in the context of record linkage [START_REF] Winkler | Using the EM Algorithm for Weight Computation in the Fellegi-Sunter Model of Record Linkage[END_REF][START_REF] Grannis | Analysis of a probabilistic record linkage technique without human review[END_REF][START_REF] Murray | Probabilistic Record Linkage and Deduplication after Indexing, Blocking, and Filtering[END_REF] and that maintains good classification performance even when the independence assumption does not hold in practice [START_REF] Zhang | Exploring conditions for the optimality of naïve bayes[END_REF][START_REF] Zhang | Naive Bayes for optimal ranking[END_REF][START_REF] Manning | Introduction to Information Retrieval[END_REF] .

Figure 4: Workflow of the linkage algorithm.

To calculate the posterior probability M FG , we further impose an assumption that the probabilities of a record in <, < (F) , being matched to either one record in ? or none sum up to 1. This essentially would assign this probability to 0.5 if < (F) is matched to two identical records in ?. For patient B from <, we define the random variable p F ∈ {0,1,2, … , D, … , s @ } indexing the patient in ? who is a match, with p F = 0 indicating that no match is found. The posterior probability of p F = D can be shown to be (see Supplementary Information for details):

M FG (>→@) = ℙ(p F = D | < (F) , ?) = exp `ℒ(FG) + logit(π n )a 1 + ∑ exp { | ℓj) (ℒ (Fℓ) + logit(π n ))

Empirical prior

Diagnosis codes prevalence

User input: matching threshold J where M n is the prior probability of matching for a random pair from < and ?. Similarly, we let the random variable ~G be the matching index of patient D from ? in < and its posterior probability given the data is:

M FG (@→>) = ℙ(~G = B | <, ? (G) ) = exp `ℒ(FG) + logit(π n )a 1 + ∑ exp { ℓj) (ℒ (ℓG) + logit(π n ))

Hyper-parameters

The posterior probabilities involve several unknown hyper-parameters including M @ T , ; V/ and ; V. for each X in {1, … , A} and M n . It is straightforward to estimate M @ T empirically using sample fractions in the data directly, as this represents the prevalence of the feature X in the dataset ?. To estimate the prior matching probability M n , we note that the observed similarity measures ℒ ⃗ = {ℒ (FG) , B = 1, … , = > , D = 1, … , = @ } follow a mixture distribution with density M n Å + (1 -M n )Ç, where Å and Ç are the respective density functions of ℒ (FG) | E FG = 1 and ℒ (FG) | E FG = 0. Since a vast majority of the pairs are not matches and the similarity measures among matched pairs are expected to follow a distribution very different from Ç, we estimate M n , as the fraction of pairs with ℒ (FG) exceeding a threshold value Ñn chosen as the furthest inflexion point of Ç on the right with Ç estimated by fitting a skew-t distribution to ℒ ⃗ . See Supplementary Information for details. Finally, the discrepancy rates {(; V/ , ; V. ), X = 1, … , A} are generally not directly estimable in datasets without gold standard labels on which pairs are matched. Again, they can instead be roughly estimated from the data using only pairs with high ℒ (FG) (for instance those exceeding Ñn), or one can provide an educated guess according to prior knowledge.

Final matching rule

For the pair < (F) and ? (G) , we estimate the posterior probability of being a match as M Ü FG ={M Ü FG (@→>) + M Ü FG (>→@) }/2, where M Ü FG (@→>) and M Ü FG (>→@) are the respective estimates of M FG (@→>)

and M FG (>→@) obtained by plugging in the estimated hyper-parameters. Finally, a pair is identified as a match if M Ü FG ≥ a, where the threshold a requires user input and its choice depends on the goal of the study. A straightforward sensible cutoff value is 0.5 (more chances of the pair being a match than not), while a more stringent threshold of 0.9 can be used if the tolerance for false matches in the study is very low.

Anonymized data usage notes

The R package ludic also includes an anonymized version of the binarized diagnosis code data from the RA1 and RA2 datasets, for both 6-year and 11-year time span. In accordance with the IRB, the ICD-9 diagnosis codes have also been masked and randomly reordered, replaced by meaningless names. Finally, the silver standard matching pairs are also provided to allow the benchmarking of methods for probabilistic record linkage using diagnosis codes.

B Fitting a skew-t distribution to observed L

There are several way to parametrize skew t distribution. We use the following:

f ST (x, m, s, ⌫, ⇠) = 2 s(⇠ + 1/⇠) r ⌫ ⌫ 2 f T ✓r ⌫ ⌫ 2 1 ⇠ ✓ x m s + µ ◆ , ⌫ ◆ 1 {( x m s +µ) 0} (4)
2

s(⇠ + 1/⇠) r ⌫ ⌫ 2 f T ✓r ⌫ ⌫ 2 ⇠ ✓ x m s + µ ◆ , ⌫ ◆ 1 {( x m s +µ)<0} (5)
where

• µ = m 1 (⇠ 1/⇠) • m 1 = 2 p ⌫ 2 (⌫ 1) • = ✓ 1 2 ◆ ⇣ ⌫ 2 ⌘ ✓ ⌫ + 1 2 ◆ • = p (1 m 2 1 )(⇠ 2 + 1/⇠ 2 ) + 2m 2 1 1)
• f T is the density function of a Student's t distribution:

f T (y, ⌫) = ✓ ⌫ + 1 2 ◆ p ⇡⌫ ⇣ ⌫ 2 ⌘ ✓ 1 + y 2 ⌫ ◆ ⌫ + 1 2
In order to estimate ⇡ 0 , we compute the first and second x-derivative of f ST to determine the inflexion point the further left of density curve. Therefore, we focus on the case where The dataset we used here has two more years of follow-up compared to the one described in Liao et al. [1]. The 2 additional years of follow-up allowed us to identify more patients when the algorithm used in Liao et al. was thus ran on a larger set of patients. Also, there was a change in the Electronic Health Record system that occurred in late 2001, so we only used records starting in 2002 onwards to ensure data coherence.

D On the failure of the Fellegi-Sunter method for matching using diagnosis codes

The Fellegi-Sunter method for record linkage [2] fails in the context of diagnosis codes because it does not di↵erentiate agreement (between both datasets) for the presence or for the absence of a diagnosis code. However, agreement for the presence of a diagnosis code is often much more informative for matching than agreement for its absence. Our proposed approach has the advantage of not only di↵erentiating the agreement weights between presence and absence, but also to automatically tune them according to the prevalence of a diagnosis (considering rarer codes as more informative). For this reason, Fellegi-Sunter will fail when using diagnosis codes, even when reaching matching weight higher than 16.6 in the context presented of Figure 1 (the lower bound given by Cook et al. [3] approach supposed to ensure a successful linkage with a probability of selecting true matches above 95%).

E EHR data for the BRASS cohort patients

As described in Iannaccone et al. [4], BRASS patients were recruited prospectively through the clinic and data are collected through patient interviews and questionnaires. When BRASS patients are enrolled, they are assigned a BRASS study ID which is linked to their medical record number (MRN).

For the study presented here, IRB approval was granted to extract the EHR data for the BRASS patients using their MRNs. Patients in the BRASS study have blood samples drawn and analyzed at a separate laboratory, that are not available in the EHR. The labs used to construct the silver standard in the study presented here are part of their routine follow-up in the Arthritis Center and are extracted from the EHR using the method described above.

Figure 1 :

 1 Figure 1: Matching accuracy on simulated data under various settings. (a) Impact of the discordance between the two datasets. (b) impact of rare codes. (c) impact of the proportion of overlapping patient between the two datasets. The figure shows the performance of the matching according to various simulation scenarios in terms of True Positive Rate (TPR) and Positive Predictive Value (PPV). F-S refers to the Fellegi-Sunter method while ludic denotes our proposed Bayesian approach.

Figure 2 :

 2 Figure 2: Matching accuracy according to : . . The figure shows the performance of the matching according to various value of the hyper-parameter ;. which models the probability of discrepancy between the two datasets. The top panel displays results when true value of ; . is 3%, while the bottom panel is for a true value of 20%.

Figure 3 :

 3 Figure 3: Histogram of diagnosis code prevalence in the RA datasets.

  prevalence: number of patients diagnosed (logarithmic scale) Number of diagnosis codes

K

  diagnosis-codes K diagnosis-codes De-identified database A De-identified database B : unaccessible direct and indirect identifiers Computation of the similarity score ℒ (&') for each possible pair {*, ,} of patients from A & B respectively, using only the K binarized diagnosis codes Final matching rule An index of the matching pairs *, , for which 0 1 &' ≥ J Number of expected matches

1 ⇠ 2 (dx 2 f 1 ⇠ 2 (

 12212 f ST (x, m, s, ⌫, ⇠) = c 1 + is a constant for x. Then d dx f ST (x, m, s, ⌫, ⇠) = d ST (x, m, s, ⌫, ⇠) =d s 1 + (x, m, s, ⌫, ⇠) < " and d 2 dx 2 f ST (x, m, s, ⌫, ⇠) < "C EHR data from Partners RA patients

Table 1 : Characteristics of the real datasets.

 1 ICD9 code of RA or had been tested for anti-cyclic citrullinated peptide. This dataset contained diagnosis codes recorded between 1/1/2002 and 12/31/2007, and these 26,681 patients together had 7,868 different unique diagnosis codes. The RA2 dataset contained 6,394 patients identified as having RA according to a phenotyping algorithm[START_REF] Liao | Electronic medical records for discovery research in rheumatoid arthritis[END_REF][START_REF] Liao | Methods to Develop an Electronic Medical Record Phenotype Algorithm to Compare the Risk of Coronary Artery Disease across 3 Chronic Disease Cohorts[END_REF] , and it contained codes recorded between 1/1/2002 and 12/31/2012 (note the same start date but a later end date in RA2). The two datasets likely have overlapping patients, and we wanted to know if we can link the data from these two studies. The 26,681 patients in the RA1 dataset together had 7,868 unique diagnosis codes from 2002 through 2007. Over the same time span as the RA1 dataset (2002 through 2007), 5,707 of the 6,394 patients in the RA2 dataset had at least one diagnosis recorded. These patients have 4,981 unique diagnosis codes through 2007, of which 4,936 also appear in the RA1 dataset. All patients in the RA2 dataset had at least one diagnosis when considering the full date range, from 2002 through 2012, with 6,086 unique diagnosis codes, of which 5,593 appear in the RA1 dataset. Patients in RA1 have a mean of 30.2 unique diagnosis codes; while patients in RA2 have a mean of 29.0 unique diagnosis codes recorded between 2002 and 2007 and 44.2 between 2002 and 2012. Table 1 summarizes these characteristics and Fig. 3 displays the histogram of diagnosis codes prevalence.

	•	ludic -0.5 cutoff
	•	ludic -0.9 cutoff
		F-S
		F-S 1-1

*The 6-year time span includes codes from 1/1/2002 through 12/31/2007, while the 11-year time span includes codes from 1/1/2002 through 12/31/2012.

Table 2 . Performance matching 2 real use case datasets. Data time span Matching method

 2

			Number of	TPR*	PPV*	Computing time
			matches
	6 years	0.5 cutoff	4,369	0.93	0.81	96 sec §
	6 years	0.9 cutoff	4,179	0.91	0.84	96 sec §
	6 years	F-S blocked	2,594,443	0.81	< 0.01	49 min §
	6 years	F-S blocked 1-1	5,696	0.38	0.26	49 min §
	6 years	F-S	-	-	-	> 4 days †
	6 years	F-S 1-1	-	-	-	> 4 days †
	11 years	0.5 cutoff	4,043	0.84	0.80	96 sec §
	11 years	0.9 cutoff	3,625	0.80	0.84	96 sec §
	11 years	F-S blocked	2,898,367	0.80	< 0.01	62 min §
	11 years	F-S blocked 1-1	6,356	0.29	0.17	62 min §
	11 years	F-S	-	-	-	> 4 days †
	11 years	F-S 1-1	-	-	-	> 4 days †

*based on the 3,831 silver standard true matches. § using a 3.5 GHz Intel Core i7 processor with 32 GB of memory available. † using a 3.6 GHz Intel Xeon 5600 series processor with 96 GB of memory available.

Acknowledgements

This work was supported by the National Institutes of Health [grant number U54-HG007963]. BPH had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.

Code availability

Open-source software implementing the proposed method is available together with complete documentation in the R package ludic (Linkage Using DIagnosis Codes) from the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=ludic.

Data Availability

The "RA1" and "RA2" data that support the findings of this study are available from the Comprehensive R Archive Network (CRAN) as part of the R package ludic, at https://CRAN.Rproject.org/package=ludic. See the Anonymized data usage notes for more information.

Author contributions

BPH developed and implemented the linkage method, analysed the data, drafted the manuscript and approved its final version. GMW contributed to the method development, to the analysis of the data, drafted the manuscript and approved its final version. KPL contributed to the analysis of the data, drafted the manuscript and approved its final version. NPP provided the claims data for the simulation study and approved the final version of the manuscript. SC contributed to drafting the manuscript and approved its final version. NAS provided the BRASS registry cohort data and approved the final version of the manuscript. PS contributed to the method development, to drafting the manuscript and approved its final version. SNM contributed to the method development, to drafting the manuscript and approved its final version. ISK contributed to the method development, to the drafting the manuscript and approved its final version. TC developed the linkage method, contributed to the analysis of the data, drafted the manuscript and approved its final version.

Competing interests

The authors declare no competing interests. ⇤ bhejblum@hsph.harvard.edu

1

A Posterior probabilities calculation details

Applying the Bayes theorem to the log-ratio, one finds:

where

k . It then follows that:

, B (j 0 ) , j 0 = 1, 2, . . . , N B )

+ log

Note that:

, B (j 0 ) , j 0 = 1, 2, . . . , N B ) P(X i = 0 | A (i)

, B (j 0 ) , j 0 = 1, 2, . . . , N B )

P(X i = `| A (i) , B (j 0 ) , j 0 = 1, 2, . . . , N B ) P(X i = 0 | A (i)

, B (j 0 ) , j 0 = 1, 2, . . . , N

exp L (i`) + logit(⇡ 0 )