Camilo Aguilar

Mathias Ortner

Josiane Zerubia

Small moving target MOT tracking

Keywords: Object Detection, Satellite Object Tracking, MOT, CNN, Deep Learning, GM-PHD

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Multi-object tracking (MOT) constitutes a key task in numerous biological, surveillance, and remote sensing applications. Particularly, remote sensing applications have experienced an increased demand for MOT: novel Earth observation satellites are able to capture large-scale images at sub-meter resolutions that allow to detect and track small vehicles and pedestrians. For instance, the Chinese Jilin-1 satellite constellation yields ground videos at resolutions up to 0.72m, and Planet's Terra Bella Skysat-1 satellite constellation records videos up to 0.8m. Large-scale ground data will contribute to improve applications that currently rely on local ground/aerial sensors such as traffic estimation [START_REF] Kaack | Truck traffic monitoring with satellite images[END_REF], traffic pattern research [START_REF] Chen | Spatial temporal analysis of traffic patterns during the covid-19 epidemic by vehicle detection using Planet remote-sensing satellite images[END_REF], or rare object tracking [START_REF] Michael | A satellite image set for the evolution of image transforms for defense applications[END_REF]. detecting and tracking small targets in satellite images. Vehicles are often represented as blobs with an average size of 6 × 6 pixels, such as in Fig. 1(b). Also, texture-rich backgrounds contain blob-like objects such as air-conditioning units and electrical transformers, which generate numerous false positives for object detection. Finally, training data is sparse and often imperfect due to the difficulty of labeling small objects across numerous frames. We approach the challenging task of small object detection by using a patch-based CNN to detect tiny moving vehicles in satellite images. We achieve this task by combining the three frame difference algorithm to drive the region selection of a CNN object detector. This approach facilitates the training of the CNN when training data is sparse because the network learns directly from regions with objects rather than whole images and it also reduces the amount of false positives as it focuses on patches with moving targets.

In addition, we employ the Probability Hypothesis Filter (PHD) to perform object tracking. The PHD filter offers a sophisticated Bayesian data-association framework [START_REF] Mahler | Multitarget Bayes filtering via firstorder multitarget moments[END_REF] that allows estimating the states of a varying number of targets given past measurements. This framework is robust to misdetections, false alarms, and appearing and disappearing objects. It also presents a significant advantage over commonly employed MOT algorithms such as probabilistic data association (PDA) [START_REF] Bar-Shalom | Tracking and data association / Yaakov Bar-Shalom[END_REF] and multiple hypothesis tracking (MHT) [START_REF] Reid | An algorithm for tracking multiple targets[END_REF] thanks to its relatively low computational complexity.

We propose an improved PHD filter that receives the measurements from the CNN object detector, discriminates surviving objects from appearing objects, and propagates the labels in time. The contribution of this paper is namely a trackby-detection approach by using 1) a novel method to detect small moving vehicles in satellite images and 2) the integration with an improved PHD filter to convert the measurements into tracks with a MOT statistical model.

The paper is organized as follows: Section 2 presents a literature review of state-of-the-art methods, Section 3 describes the proposed approach, first discussing the patch selection mechanism and then developing on the improvements of the PHD filter. Section 4 discusses the experiment section and Section 5 develops the conclusions and insights.

RELATED WORK

Common approaches for satellite object tracking involve using motion information or correlation filters. Du et al. utilized a combination of frame-difference and correlation filter to track a single object [START_REF] Du | Object tracking in satellite videos by fusing the kernel correlation filter and the three-frame-difference algorithm[END_REF], Xuan et al. employed correlation filters, together with Kalman filters and linear motion to make a tracker robust to occlusions [START_REF] Shiyu Xuan | Object tracking in satellite videos by improved correlation filters with motion estimations[END_REF]. These approaches show promising results for satellite object tracking but are limited to single object tracking, require an initialization, and they generally focus on larger targets (planes/trains vs cars/pedestrians). Jiao et al. performed a survey [START_REF] Jiao | New generation deep learning for video object detection: A survey[END_REF] with new generation deep learning tracking approaches where numerous state-of-the-art deep learning approaches consist in learning deep features and matching them in future/past frames. However, these approaches would experience severe drawbacks when tracking small objects in satellite videos: satellite videos do not have large training datasets, targets do not exhibit distinct features, and the objects of interest are tiny relative to the satellite's field of view.

Common baselines for online multi-target trackers include Tracktor++ [START_REF] Bergmann | Tracking without bells and whistles[END_REF] and DeepSORT [START_REF] Wojke | Simple online and realtime tracking with a deep association metric[END_REF]. These methods obtain remarkable results in surveillance camera videos by using CNNs for detection and Kalman filters for tracking. However, their CNNs are based on Faster-RCNN and would have difficulties learning directly from the whole image.

Wei et al. presented the closest state-of-the-art approach comparable to our problem in [START_REF] Ao | Needles in a haystack: Tracking city-scale moving vehicles from continuously moving satellite[END_REF]. The authors used motion information, several post-processing steps, and multiple Kalman filters to detect and track small objects in satellite videos. This approach proved promising but would underperform if cluttered objects, slow moving objects, and the data association were done using as a post processing step with the Hungarian algorithm [START_REF] Kuhn | The Hungarian method for the assignment problem[END_REF]. In our approach, we incorporate deep learning for object detection, and we opt for a direct-MOT tracking with the PHD filter rather than using several individual Kalman filters.

PROPOSED TRACKING ALGORITHM

We use the Random Finite Set (RFS) framework to model the collection of state vectors as

X k = {x 1 k , x 2 k , . . . x N k k }
, and the frame measurements as

Z k = {z 1 k , z 2 k , . . . , z M k k }
, where N k and M k denote the number of objects and measurements at time k. We define a single target state vector as

x j k = [p j k , v j k , a j k ]
T where j and k are the object and time indices respectively, and p j k , v j k , a j k ∈ R 2 denote the target's position, velocity, and acceleration vectors. Similarly, we define a measurement vector z i k = [p i k ] T , where pi k ∈ R 2 denotes the coordinates obtained by an object detector for the i th measurement.

During each frame, we use the CNN-based object detector to obtain new measurements, i.e. Z k , and then we use the PHD filter to approximate the posterior multi-target RFS distribution: p k|k (X k |Z 1:k ). The summary of the proposed approach is depicted in Fig. 2.

Patch-based CNN Object Detector

We propose a patch-driven approach to segment moving vehicles. This approach filters overwhelming image information and performs inference on relevant areas of interest. We first apply a motion-based detector to obtain rough object locations. In our experiments, the frame difference algorithm proved to be an effective and computationally lightweight operation and it is defined as:

∆I k (i, j) = I k (i, j) -I k-1 (i, j) (1)
3F D k (i, j) = 1 if |∆I k (i, j)| + |∆I k+1 (i, j)| > T k 0 otherwise (2)
Where I k (i, j) denotes the image intensity at coordinates i, j during time k and T k is an adaptive threshold to binarize The three-frame difference segmentation outputs rough object locations but it does not regularize objects by area or shape and it misses slow-moving objects. For example, Fig. 3 (c) shows nearby object merging, measurement duplication for a single target, and missed objects toward the bottom of the patch. Therefore we upsample the patch and apply Faster-RCNN [START_REF] Shaoqing Ren | Faster R-CNN: Towards real-time object detection with region proposal networks[END_REF] to obtain precise object locations. Fig. 3 (d) shows that Faster-RCNN detects and discriminates nearby objects, and corrects for the missed detections. In fact, with visual inspection, the proposed object detection provides a better fit than the hand-annotated labels.

Finally, we convert the patch coordinates back to global coordinates and perform global non-maximum suppression (NMS) to discard repeated objects. This approach allows us to exploit low-availability of labeled data by focusing the network training on patches rather than whole images/videos. The summarized pipeline for this approach is depicted in Fig. 4.

Gaussian Mixture Probability Hypothesis Filter

We use the PHD filter to approximate the target's posterior RFS p k|k (X k |Z 1:k ) by recursively propagating its first order moment, namely its intensity function, D k (x). The GM-PHD, proposed by Vo and Ma [START_REF] Ba | The Gaussian mixture probability hypothesis density filter[END_REF] approximates a closed form solution to the PHD recursion with a weighted Gaussian mixture. It assumes that the intensity function from the previous frame has the form:

D k-1 (x) = J k-1 j=1 w j k-1 N (x; m j k-1 , P j k-1 ) (3)
Where ω j k-1 , m j k-1 , and P j k-1 are the weight, mean, and covariance for the j th component and J k-1 is the number of components. The closed form solution for the GM-PHD prediction step is given by the equation:

D k|k-1 (x) = λ(x) + p s J k-1 j=1 ω k-1 N (x; F m j k-1 , Q + F P j k-1 F T )
Where F and Q are the transition and motion covariance matrices in the same format as in [START_REF] Pace | Stochastic models and methods for multiobject tracking[END_REF], λ(x) is the birth RFS intensity and p s is a hyper-parameter to denotes the survival probability. Once a new set of measurements Z k arrives from the object detector, we update the GM-PHD posterior following the equation:

D k (x) = (1 -p D ) D k|k-1 (x)+ z∈Z k J k|k-1 j=1 ω j k (z) N (x; m j k|k (z), P j k|k )
Where D k|k-1 (x) denotes the predicted GM components and p D is a hyper-parameter to denote the detection probability. The terms m j k|k (z), P j k|k are the target-measurement association mean and covariance, which are calculated using the Kalman equations, and ω j k (z) denotes the updated weight for the data association and is defined as:

ω j k (z) = p D ω j k|k-1 l j k (z) κ k (z) + p D J k|k-1 i=1 w i k|k-1 l i k (z) (4)
Where κ k (z) denotes the clutter process intensity and l i k (z) denotes the target-measurement association likelihood.

PHD Filter Improvements

The PHD filter requires an accurate birth intensity λ(x) modeling in order to deal with a varying number of targets. We estimate the birth intensity λ(x) by using information from the CNN object detector. We use the decomposed measurement set similar to [START_REF] Fu | GM-PHD filter based online multiple human tracking using deep discriminative correlation matching[END_REF], defined as

Z k = Z s k ∪ Z b k , where Z s k
denotes the surviving targets and Z b k denotes the appearing (birth) targets. We classify a new measurement as surviving measurement (i.e. z i k ∈ Z s k ) if it is the closest measurement to the location of a predicted GM component and if this distance is less than a threshold T birth , otherwise we classify the measurement as a birth measurement (i.e. z i k ∈ Z b k ). We approximate λ(x) by initializing a new mixture component centered around the spatial coordinates from z i k ∈ Z b k and its velocity and acceleration components set to zero.

In addition, the PHD filter does not keep tracks of labels by default, therefore we added the label tracking scheme proposed by [START_REF] Panta | Data association and track management for the Gaussian mixture probability hypothesis density filter[END_REF]. This approach propagates labels in time in a tree-like structure without affecting the filter's performance.

EXPERIMENTS

Evaluation Metrics

We evaluate our approach using object detection metrics, tracking quality measures, and the ClearMOT [START_REF] Bernardin | Evaluating Multiple Object Tracking Performance: The CLEAR MOT Metrics[END_REF] framework. We use the Hungarian algorithm [START_REF] Kuhn | The Hungarian method for the assignment problem[END_REF] at each frame to match predicted hypotheses to ground truth labels and we call a hypothesis T P if it is located within a distance d T hrs from the ground truth label. We choose d T hrs = 20 pixels to account for the noise introduced by the image stabilization procedure. Similarly, we call F P any unmatched hypothesis and F N any unmatched ground truth label.

We report object detection metrics by using precision (P), recall (R), F 1 score (F1), and the Jacard similarity index (J). These scores provide an assesment to the quality of object detection and were computed by counting TPs, FPs, and FNs over all the frames. In addition, we report tracking quality measures where we report the percentage of mostly tracked (MT) trajectories and mostly lost (ML) trajectories. We classify tracks as MT if at least 80% of the ground truth trajectory is recovered, and ML if less than 20% of the ground truth track is recovered.

Finally, we follow the ClearMOT [START_REF] Bernardin | Evaluating Multiple Object Tracking Performance: The CLEAR MOT Metrics[END_REF] as it has become a popular and robust metric for tracking algorithms. We report the multiple object tracking accuracy (MOTA) which considers FPs, FNs, and identity switches (IDs), and we also report the multiple object tracking precision (MOTP), which considers the average distance error between the detected objects and the ground truth objects.

Satellite Dataset

The experimental satellite video can be obtained by CGSTL (available at https://mall.charmingglobe.com). This video contains 3071×4096 pixels and it represents the city of Valencia, Spain. The video contains 580 frames, lasts 29 seconds, and was imaged at 20 fps at a resolution of 1.0m by the Jilin-1 Satellite on March 7, 2017. The labels were created by the authors of [START_REF] Ao | Needles in a haystack: Tracking city-scale moving vehicles from continuously moving satellite[END_REF] and they contain labels for one every 10 frames. Specifically, three regions of the image were actively labeled for evaluation purposes. These regions span 500 × 500 pixels and are shown in [START_REF] Ao | Needles in a haystack: Tracking city-scale moving vehicles from continuously moving satellite[END_REF]. The approximate locations are: [520, 1616] for area 1, [450, 2810] for area 2, and [1074, 1895] for area 3. Due to the low availability of ground truth annotations, we tested our method on area 1 and area 2, and trained the network using 2D patches from area 3 and areas outside area 1 and area 2.

The network was trained using an NVIDIA Quadro T2000 with 4 GBs of RAM and the training time was 6 hours. In addition, we performed a pre-processing step of image stabilization by using the first frame's strong edges as reference. The stabilization step was fundamental for the three-frame difference algorithm as it contributed to significantly reduce noise introduced by the satellite motion. Experimentally, we set the three-frame difference threshold parameter c = 0.15, and the GM-PHD filter parameters p D = 0.90 and p S = 0.95.

Table 1 denotes numerical comparisons for our proposed approach vs competing methods, namely Needles in a Haystack (NIAH) [START_REF] Ao | Needles in a haystack: Tracking city-scale moving vehicles from continuously moving satellite[END_REF], ViBe [START_REF] Yang | An improved ViBe for video moving object detection based on evidential reasoning[END_REF] ture Models (GMM) [START_REF] Nurhadiyatna | Background subtraction using Gaussian mixture model enhanced by hole filling algorithm (GMMHF)[END_REF]. We implemented these methods following the directives in [START_REF] Ao | Needles in a haystack: Tracking city-scale moving vehicles from continuously moving satellite[END_REF] and used SORT [START_REF] Bewley | Simple online and realtime tracking[END_REF] as the data-association step. We stabilized the images and its respective ground truth annotations in order to reduce noise across all the methods; hence our results are slightly different than the ones reported in [START_REF] Ao | Needles in a haystack: Tracking city-scale moving vehicles from continuously moving satellite[END_REF].

Our approach outperforms competing methods in comprehensive object detection and object tracking metrics. For instance, NIAH [START_REF] Ao | Needles in a haystack: Tracking city-scale moving vehicles from continuously moving satellite[END_REF] and GMM [START_REF] Nurhadiyatna | Background subtraction using Gaussian mixture model enhanced by hole filling algorithm (GMMHF)[END_REF] obtain better precision scores in Area 1, but their recall values in Area 1 are significantly lower. The addition of Faster RCNN allows our approach to discriminate nearby objects that are merged or not detected in NIAH [START_REF] Ao | Needles in a haystack: Tracking city-scale moving vehicles from continuously moving satellite[END_REF] or GMM [START_REF] Nurhadiyatna | Background subtraction using Gaussian mixture model enhanced by hole filling algorithm (GMMHF)[END_REF] such as in Fig 3 .  Area 2 presents numerous tiny objects such as motorcycles or small cars that are undetected by all object detectors; hence our approach obtains similar recall values. However, the combination of our object detector and GM-PHD filter allows our method to reduce FP, hence obtaining better precision scores. The trade-offs and advantages of the proposed method are denoted in Table 1 where our approach obtains the highest F 1 scores and J scores for both areas.

The tracking performance is also outlined in the number of mostly tracked (MT) trajectories, the MOTA score, and MOTP score. Our approach obtains a slightly worse MOTP score (1 pixel difference) due to the modeling of uncertainty in the GM-PHD-filter; however, our method obtains more MT trajectories, lower identity switches (IDs), and better MOTA thanks to the employment of the GM-PHD filter. Particularly, the GM-PHD reduces clutter detections and keeps track of undetected objects for several frames.

CONCLUSION

In this paper, we presented a track-by-detection approach combining a patch-based CNN object detector and the PHD filter. The patch selection mechanism contributes to filter unnecessary information and to train the neural network in a dataset with sparse labels. In addition, we use the GM-PHD algorithm to track multiple targets while reducing the amount of clutter. Our approach presents an increased computational burden due to the patch selection and CNN inference, but its results outperform competing methods in both object detection and tracking scores for the challenging task of detecting small satellite objects.

ACKNOWLEDGMENT

The authors are grateful to the OPAL infrastructure from Université Côte d'Azur for providing resources and support. Additionally, the authors would like to thank BPI France for the financial support under the LiChiE contract.

Fig. 1 (Fig. 1 .

 11 Fig. 1. Sample Satellite Image and Tracks.

Fig. 2 .

 2 Fig. 2. Proposed Approach.

Fig. 3 .

 3 Fig. 3. Object Detector Comparison.

Fig. 4 .

 4 Fig. 4. Patch Based Object Detection.

Fig. 5 .

 5 Fig. 5. Area 1 Tracking Results.

Table 1 .

 1 and Gaussian Mix-Tracking scores. The arrow's direction represents better scores.

	Area Method	FP ↓	FN ↓ IDs↓ Precision↑ Recall↑	F1↑	J↑	MT↑	ML↓	MOTA↑ MOTP↓
		Proposed	1560 3755 264	89.12	77.29	76.01 70.63 65.31% 10.20%	66.3	3.64
	1	NIAH [12] ViBe [20]	935 53993 3908 310 7877 901	90.25 18.96	52.37 76.37	56.64 49.56 40.99 17.91 51.02 % 4.08% 6.12% 6.12 %	41.30 -252.0	3.66 3.64
		GMM [21]	556	5687 210	95.13	65.61	74.32 63.48 44.90 % 18.37%	61.0	2.50
		Proposed	108	441	32	87.61	63.40	65.09 58.19 68.96% 17.24%	51.80	3.23
	2	NIAH [12] ViBe [20]	197 38218 2046 2953 121 82	93.77 9.20	50.08 65.42	60.55 48.47 20.69% 34.48% 18.69 8.76 55.17% 20.69%	44.71 -582.0	2.35 3.69
		GMM [21] 1696 1915	43	70.00	67.39	40.07 52.29 58.62% 27.58%	37.80	2.25


Franco Quezada
email: franco.quezada@usach.cl

Céline Gicquel
email: celine.gicquel@lri.fr

Safia Kedad-Sidhoum
email: safia.kedadsidhoum@cnam.fr

New valid inequalities for a multi-echelon multi-item lot-sizing problem with returns and lost sales

Keywords: Production planning, Lot-sizing, Remanufacturing, Mixedinteger linear programming, Valid inequalities

This work studies a multi-echelon multi-item lot-sizing problem with remanufacturing and lost sales. The problem is formulated as a mixed-integer linear program. A new family of valid inequalities taking advantage of the problem structure is introduced and used in a customized branch-and-cut algorithm. The provided numerical results show that the proposed algorithm outperforms both the generic branch-andcut algorithm embedded in a standard-alone mathematical solver and a previously published customized branch-and-cut algorithm.

Introduction

Industrial companies face an increasing pressure from customers and governments to become more environmentally responsible and mitigate the environmental impact of their products. One way of achieving this objective is to remanufacture the products once they have reached their end of life. Remanufacturing is defined as a set of processes transforming used products into like-new products, mainly by rehabilitating damaged components. By reusing the materials and components embedded in used products, remanufacturing both contributes in reducing pollution emissions and natural resource consumption, making production processes more environment-friendly.

The present work considers a remanufacturing system involving three production echelons: disassembly of used products brought back by customers, refurbishing of the recovered parts and reassembly into like-new finished products. We aim at optimizing the production planning for the corresponding three-echelon system over a multi-period horizon. Within a remanufacturing context, production planning includes making decisions on the used products returned by customers, such as how much and when used products should be disassembled, refurbished or reassembled in order to build new or like-new products. The main objective is to meet customers' demand for the remanufactured products in the most cost-effective way.

We consider the case where the production on a machine requires setup operations such as machine calibration and incurs fixed setup costs. As a naive perception, to reduce these setup costs, production should be run using large lot sizes. However, this generates desynchronized patterns between the customers' demand and the production plan leading to costly high levels of inventory. Lotsizing models thus aim at reaching the best possible trade-off between minimizing the setup costs and minimizing the inventory holding costs, taking into account both the customers' demand satisfaction and the practical limitations of the system. In the present work, we investigate the problem of reaching the best possible trade-off between setup and inventory holding costs within a remanufacturing environment and introduce an additional lost sales cost to be paid when the customers' demand is not satisfied on time. We thus study a 3-echelon lot-sizing problem with returns and lost sales.

Only a few works have addressed such multi-echelon production systems through exact solution approaches. A first attempt at tackling this difficulty can be found in [START_REF] Quezada | A multi-stage stochastic integer programming approach for a multi-echelon lot-sizing problem with returns and lost sales[END_REF]. Quezada et al. [START_REF] Quezada | A multi-stage stochastic integer programming approach for a multi-echelon lot-sizing problem with returns and lost sales[END_REF] considered the problem in a stochastic setting, taking into account uncertainties on the problem input parameters. They proposed a multi-stage stochastic approach based on the use of scenario trees. The problem was formulated as a MILP and solved through a new customized branch-and-cut algorithm. This algorithm relied on valid inequalities focused on strengthening the formulation of the single-echelon uncapacitated lot-sizing sub-problems embedded in the main problem. Although this approach was successful at providing near optimal solutions for small to medium size instances, some numerical difficulties were encountered to solve the larger instances. Intuitively, this difficulty might be partly due to the fact that the valid inequalities used to strengthen the formulation considered uncapacitated single-echelon subproblems. They did not take into account the fact that, even if the production resources are assumed uncapacitated, the amount of products that can be processed on a resource at a given time period is limited among others by the amount of available used products returned up to this time period and by the yield of the disassembly process, i.e. by the proportion of disassembled parts that are in a sufficiently good state to be refurbished and reused in a remanufactured product. Hence, using valid inequalities taking into account this aspect of the problem might contribute in further strengthening its MILP formulation and decrease the computational effort needed to solve large-size instances.

To the best of our knowledge, the formulation of valid inequalities that explicitly take into account the impact of a limited returns quantity on the production plan has not yet been studied for a multi-echelon remanufacturing system.

The present work aims at partially closing this gap by proposing new valid inequalities for this problem. However, in view of the theoretical and numerical difficulties encountered when developing new valid inequalities, we focus on the deterministic variant of the problem. Our contributions are thus twofold. First, we propose a new family of valid inequalities for the problem under study. These valid inequalities can be seen as an extension of valid inequalities previously known for the uncapacitated single-echelon lot-sizing problem with lost sales (the (k, U ) inequalities first introduced in [START_REF] Loparic | The uncapacitated lot-sizing problem with sales and safety stocks[END_REF]) to take into account, at each echelon of the studied multi-echelon production system, the constraints on the production plan coming from the limited availability of the returns. We prove that these new valid inequalities are at least as strong as the previously known inequalities. Second, we develop a branch-and-cut algorithm based on the newly proposed valid inequalities and seek to assess its computational performance by comparing it with the one of a stand-alone mathematical programming solver and the one of a branch-and-cut algorithm based on previously known valid inequalities. The numerical results show the usefulness of the proposed inequalities at solving the problem under study.

The remainder of this paper is organized as follows. We first provide a brief overview of the related literature in Section 2. The problem description, together with its MILP formulation, are provided in Section 3. We then present the proposed new family of valid inequalities in Section 4. Computational results are summarized in Section 5. Finally, Section 6 gives a conclusion and some research perspectives.

Related works

Throughout the last decade, several works sought to strengthen the MILP formulation of single-echelon lot-sizing problems involving remanufacturing, either through extended reformulations or through valid inequalities.

Helmrich et al. [START_REF] Retel Helmrich | Economic lot-sizing with remanufacturing: complexity and efficient formulations[END_REF] discussed several MILP formulations of the uncapacitated single-item single-echelon lot-sizing problem with remanufacturing and introduced new valid inequalities by adapting the previously known (l, S, W W ) proposed by [START_REF] Pochet | Polyhedra for lot-sizing with wagner-whitin costs[END_REF] to their problem. The inequalities developed in [START_REF] Retel Helmrich | Economic lot-sizing with remanufacturing: complexity and efficient formulations[END_REF] are based on the assumption that all returned products are either processed or kept in stock and do not consider the possibility that some of the returns may be discarded in case of an unbalance between returns quantity and demand. They can therefore not be directly used for the problem under study here. Similarly, [START_REF] Cunha | On multi-item economic lot-sizing with remanufacturing and uncapacitated production[END_REF] proposed a multi-commodity reformulation and a new set of valid inequalities for this problem. In particular, they further strengthened the (l, S, W W ) inequalities presented in [START_REF] Retel Helmrich | Economic lot-sizing with remanufacturing: complexity and efficient formulations[END_REF] by considering that the amount of finished products remanufactured in a given period t is limited by the cumulative quantity of returned products brought back up to t. Ali et al. [START_REF] Ali | A theoretical and computational analysis of lot-sizing in remanufacturing with separate setups[END_REF] enriched the previous works by highlighting a theoretical property with regards to the equivalence of the shortest path and facility location reformulations. They also carried out a polyhedral analysis of a related sub-problem based on the single node fixed-charge network problem, proving the validity of several flow cover inequalities and their facetdefining conditions as well. Akartunali and Arulselvan [START_REF] Akartunalı | Economic lot-sizing problem with remanufacturing option: complexity and algorithms[END_REF] studied both the uncapacitated and capacitated variants of this single-item single-echelon lot-sizing problem. They showed that the uncapacitated problem cannot have a fully polynomial time approximation scheme (FPTAS) and provided a pseudo-polynomial algorithm to solve the problem. They also provided valid inequalities based on the flow-cover inequalities for the problem with a limited production capacity. Finally, we refer the reader to [START_REF] Brahimi | Single-item dynamic lot-sizing problems: An updated survey[END_REF] for a recent survey on single-item lot-sizing problems with remanufacturing.

We note that all the above mentioned works focus on single-item singleechelon problems and do not consider the fact that remanufacturing may involve several processing steps, i.e. several production echelons, in order to transform the returned used products into like-new products. Moreover, these works consider hybrid manufacturing/remanufacturing systems and assume that, thanks to the presence of an uncapacitated manufacturing system, it will always be possible to satisfy the demand for finished products on time. In contrast, we investigate a pure remanufacturing system and consider that the demand will be lost in case the quantity and/or the quality of the returned products are insufficient to meet it on time. This means among others that the inequalities introduced in [START_REF] Retel Helmrich | Economic lot-sizing with remanufacturing: complexity and efficient formulations[END_REF], [START_REF] Cunha | On multi-item economic lot-sizing with remanufacturing and uncapacitated production[END_REF] and [START_REF] Ali | A theoretical and computational analysis of lot-sizing in remanufacturing with separate setups[END_REF] are not valid for our problem and that new valid inequalities taking into account its specific features are needed.

3 Problem description and modeling

Production system

We consider a remanufacturing system comprising three main production echelons: disassembly, refurbishing and reassembly. We seek to plan the production activities in this system over a horizon comprising a discrete set T = {1, .., T } of periods. The system involves a set I of items. Among these ones, item i = 0 represents the used products returned by customers in limited quantities at each period. A used product is composed of I parts. Let α i be the number of parts i embedded in a used product. The returned products are first disassembled to obtain a set I r = {1, ..., I} of recoverable parts. Due to the usage state of the used products, some of the parts obtained during disassembly have to be discarded. In order to reflect the variations in the quality of the used products, i.e. the yield of the disassembly process, we let π t i denote the proportion of parts which will be recoverable at each time period t for each item i = {1, . . . , I}. The recoverable parts are then refurbished on dedicated refurbishing processes. The set of I s = {I + 1, ..., 2I} of serviceable parts obtained after refurbishing are reassembled into remanufactured products which have the same bill-of-material as the used products. These remanufactured products, indexed by i = 2I + 1, are used to satisfy the dynamic demand of customers.

The system comprises a set P = {0, ..., I + 1} of production processes: p = 0 corresponds to the disassembly process, p ∈ {1, ..., I} corresponds to the process refurbishing the recoverable part indexed by p into the serviceable part indexed by p + I and p = I + 1 corresponds to the reassembly process. All these processes are assumed to be uncapacitated. However, the system might not be able to satisfy the customer demand on time due to part shortages if there are not enough used products returned by customers or if their quality is low. In this situation, the corresponding demand is lost incurring a high penalty cost to account for the loss of customer goodwill. Moreover, some used products and recoverable parts are allowed to be discarded. This option might be useful in case more used products are returned than what is needed to satisfy the demand for remanufactured products and in case there is a strong unbalance between the part-dependent disassembly yields leading to an unnecessary accumulation in inventory of the easy-to-recover parts.

All input parameters of the problem are time-dependent: r t denotes the quantity of collected used products, d t the customers' demand and π t i the proportion of recoverable parts i ∈ I r obtained by disassembling one unit of returned product at period t. As for the costs, for each period t, we have the setup cost f t p for process p ∈ P, the unit inventory cost h t i for part i ∈ I, the unit lost-sales penalty cost l t , the unit cost q t i for discarding item i ∈ I r ∪ {0} and the unit cost g t for discarding the unrecoverable parts obtained while disassembling one unit of returned product.

Natural formulation

In order to build a mathematical model for the problem, we introduce the following decision variables at time period t ∈ T : X t p the quantity of parts processed on process p ∈ P, Y t p ∈ {0, 1} the setup variable for process p ∈ P, S t i the inventory level of part i ∈ I, Q t i the quantity of part i ∈ I r ∪ {0} discarded and L t the lost sales of remanufactured products. This leads to the following MILP model.

min t∈T p∈J f t p Y t p + i∈I h t i S t i + l t L t + i∈Ir∪{0} q t i Q t i + g t X t 0 (1) X t p ≤ M t p Y t p ∀p ∈ J , ∀t ∈ T (2)
S t 0 = S t-1 0 + r t -X t 0 -Q t 0 ∀t ∈ T (3)
S t i = S t-1 i + π t i α i X t 0 -X t i -Q t i ∀i ∈ I r , ∀t ∈ T (4)
S t i = S t-1 i + X t i-I -α i X t I+1 ∀i ∈ I s , ∀t ∈ T (5)
S t 2I+1 = S t-1 2I+1 + X t I+1 -d t + L t ∀t ∈ T (6)
S 0 i = 0 ∀i ∈ I (7)
S t i ≥ 0 ∀i ∈ I, ∀t ∈ T (8) X t p ≥ 0, Y t p ∈ {0, 1} ∀p ∈ J , ∀t ∈ T (9)
The objective function (1) aims at minimizing the total remanufacturing cost over the whole planning horizon, i.e., the sum of the setup, inventory holding, lost sales and disposal costs. Constraints (2) link the production quantity variables to the setup variables. Constraints (3)-( 6) are the inventory balance constraints. More specifically, Constraints (3) ensure that any returned product is either disassembled, discarded, or kept in stock. Constraints (4) guarantee that any item obtained from the disassembly process is either refurbished, discarded, or kept in stock. Constraints [START_REF] Cunha | On multi-item economic lot-sizing with remanufacturing and uncapacitated production[END_REF] ensure that any refurbished item is either used in the reassembly process or kept in stock. Constraints [START_REF] Guan | Cutting planes for multistage stochastic integer programs[END_REF] ensure that any remanufactured/finished product is either used to satisfy the demand or kept in stock and that, if there is not enough remanufactured products to satisfy the demand, the unsatisfied demand is lost. Without loss of generality, we assume that the initial inventory, S 0 i , is set to 0 for each item i ∈ I ( see Constraints [START_REF] Leung | Facets and algorithms for capacitated lot sizing[END_REF]). Finally, Constraints (8)-( 9) provide the domain of the decision variables.

Note that the value of each constant M t p can be set by using an upper bound on the quantity that can be processed on process p at each time period t. This quantity is limited by two elements: the availability of the used products already returned by customers and the future demand for remanufactured products. We thus have, for each period t:

-M t 0 = min 1≤κ≤t r κ , t≤κ≤T d κ min i=1,...,I π t i -M t p = α p min 1≤κ≤t r κ πκ,t p , 1≤κ≤T d κ , for p ∈ I r . -M t I+1 = min min p∈Ir 1≤κ≤t r κ πκ,t p , t≤κ≤T d κ
where πκ,t p = argmax{π θ p , θ = κ, . . . , t} denotes the maximum disassembly yield that can be obtained for recoverable item p over the time interval [κ, t].

Echelon stock reformulation

We now provide a reformulation of the problem using the echelon stock concept [START_REF] Pochet | Production planning by mixed integer programming[END_REF]. The echelon stock of a product in a multi-echelon production system corresponds to the total quantity of the product held in inventory, either as such or as a component within its successors in the bill-of-material. We thus denote by E t i the echelon stock level of item i ∈ I \ {0} at the end of period t. Replacing variables S t i by variables E t i in Problem ( 1)-( 9) leads to the following MILP formulation:

min t∈T p∈J f t p Y t p + h t i S t 0 + i∈I\{0} eh t i E t i + i∈Ir∪{0} q t i Q t i + g t 0 X t 0 ( 10
)
X t p ≤ M t p Y t p ∀p ∈ J , ∀t ∈ T ( 11
)
S t 0 = S t-1 0 + r t -X t 0 -Q t 0 ∀t ∈ T ( 12
)
E t i = E t-1 i + π t i α i X t 0 -α i (d t -L t ) -Q t i ∀i ∈ I r , ∀t ∈ T ( 13
)
E t i = E t-1 i + X t i-I -α i (d t -L t ) ∀i ∈ I s , ∀t ∈ T ( 14
)
E t 2I+1 = E t-1 2I+1 + X t I+1 -d t + L t ∀t ∈ T (15)
S 0 0 = 0 (16)

E 0 i = 0 ∀i ∈ I \ {0} (17)
E t i -E t I+i ≥ 0 ∀i ∈ I r , ∀n ∈ T (18)
E t i -α i E t 2I+1 ≥ 0 ∀i ∈ I s , ∀n ∈ T (19)
E t i ≥ 0 ∀i ∈ I, ∀t ∈ T (20)
S t 0 , L t ≥ 0 ∀t ∈ T (21)
X t p ≥ 0, Y t p ∈ {0, 1} ∀p ∈ J , ∀t ∈ T (22)
The objective function [START_REF] Pochet | Polyhedra for lot-sizing with wagner-whitin costs[END_REF] aims at minimizing the total cost over the whole planning horizon. Constraints [START_REF] Pochet | Production planning by mixed integer programming[END_REF] link the production quantity variables to the setup variables. Constraints ( 12)-(15) are the inventory balance constraints. Constraints [START_REF] Quezada | A multi-stage stochastic integer programming approach for a multi-echelon lot-sizing problem with returns and lost sales[END_REF] use the classical inventory variables, whereas Constraints ( 13)-(15) make use of the echelon inventory variables. Constraints ( 16)-(17) translate the fact that the initial inventory of each item is assumed to be equal to 0. Constraints (18)-(19) ensure consistency between the echelon inventory at the different echelons of the bill-of-material and guarantee that the physical inventory of each product, S t i , remains non-negative for all i ∈ I . Finally, Constraints (20)-( 22) define the domain of the decision variables.

The use of the echelon stock reformulation ( 11)-( 22) enables us to decompose the initial problem into a series of single-echelon sub-problems by relaxing the linking constraints (18)-( 19). Each of these sub-problems is an uncapacitated single-echelon single-item lot-sizing problem with lost sales, whose formulation can be strengthened by the (k, U ) valid inequalities proposed in [START_REF] Loparic | The uncapacitated lot-sizing problem with sales and safety stocks[END_REF]. We refer the reader to [START_REF] Quezada | A multi-stage stochastic integer programming approach for a multi-echelon lot-sizing problem with returns and lost sales[END_REF] for a detailed description of each subproblem and the single-echelon (k, U ) inequalities applied to each subproblem. Nonetheless, this decomposition into single-echelon uncapacitated sub-problems overlooks the fact that the production on each process at a given period is limited by the amount of used products returned up to this period. Therefore, in what follows, we investigate a class of valid inequalities in which these aspects of the problem are explicitly considered.

Single-echelon ( , k, U ) inequalities

We now seek to strengthen the single-echelon (k, U ) inequalities investigated in [START_REF] Loparic | The uncapacitated lot-sizing problem with sales and safety stocks[END_REF] and [START_REF] Quezada | A multi-stage stochastic integer programming approach for a multi-echelon lot-sizing problem with returns and lost sales[END_REF] by considering the limited quantity of returned products available at each time period in the system. The ( , k, U ) inequalities are defined as follows:

Proposition 1. Let 0 ≤ ≤ k ≤ T be two periods of the planning horizon. Let U ⊆ {k + 1, ..., T } be a subset of periods and t * = max{τ : τ ∈ U } be the last time period belonging to U . The following inequalities are valid for Problem (11)-( 22):

S 0 π ,t * i + α -1 i E k i + k<t≤t * φ t i Y t 0 ≥ t∈U (d t -L t ) ∀i ∈ Ir (23) S 0 π ,t * i + α -1 i (E i -E i+I ) + α -1 i E k i+I + k<t≤t * φ t i Y t i ≥ t∈U (d t -L t ) ∀i ∈ Ir (24) S 0 π ,t * i + (α -1 i E i -E 2I+I ) + E k 2I+1 + k<t≤t * φ t i Y t I+1 ≥ t∈U (d t -L t )
(d t -L t ).
We consider two cases.

-Case 1: there is no production on p = 0 over the interval [k + 1; τ -1]

In this case, Y t 0 = 0 and X t 0 = 0 for all periods t in [k+1, τ -1]. As no disassembly occurs over [k + 1, τ -1], all the recoverable items needed to satisfy the demand over this time interval, and in particular needed to satisfy t∈U ;t≤τ -1

(d t -L t ),
should already have been disassembled previously and be in stock at the end of period k. This gives α -1 i E k i ≥ t∈U :t≤τ -1 (d t -L t ). We thus have:

S 0 π ,t * i + α -1 i E k i + k<t≤t * φ t i Y t 0 ≥ α -1 i E k i + φ τ i Y τ 0 ≥ t∈U (d t -L t )
Inequality ( 23) is thus valid in this case.

-Case 2: there is at least one production period on p = 0 over interval [k+1; τ -1]

Let θ be the last period of production on p = 0 over the interval [k + 1; τ -1]. By definition of θ, we have:

φ θ i = <ν≤θ (r ν πν,θ i )
. By summing up the inventory balance constraints (13) over periods k + 1,. . . ,τ -1 and using the fact that variables E k i and Q t i , ∀t = k + 1, . . . , τ -1, are non-negative, we have:

α -1 i E k i + τ -1 t=k+1 π t i X t 0 ≥ τ -1 t=k+1 (d t -L t ) (26)
By definition of τ , θ and , we have:

τ -1 t=k+1 π t i X t 0 = θ t=k+1 π t i X t 0 ≤ θ t= +1 π t i X t 0 ( 27
)
This gives:

α -1 i E k i + θ t= +1 π t i X t 0 α -1 i ≥ E k i + τ -1 t=k+1 π t i X t 0 (28) ≥ τ -1 t=k+1 (d t -L t ) (29) ≥ t∈U :t≤τ -1 (d t -L t ) (30)
We now compute an upper bound of θ t= +1 π t i X t 0 . This one is obtained by first computing the linear combination [START_REF] Quezada | A multi-stage stochastic integer programming approach for a multi-echelon lot-sizing problem with returns and lost sales[END_REF] t . This gives:

θ t= +1 πt,θ i × (
θ t= +1 πt,θ i S t 0 = θ t= +1 πt,θ i S t-1 0 + r t -X t 0 -Q t 0 ( 31
)
By the non-negativity of variables Q t 0 and S t 0 and the fact that πt,θ i ≥ πt+1,θ i , we have:

θ t= +1 π t i X t 0 ≤ θ t= πt,θ i X t 0 (32) ≤ θ t= +1 πt,θ i S t-1 0 - θ t= +1 πt,θ i S t 0 + θ t= +1 πt,θ i r t (33) ≤ π ,θ i S 0 + θ t= +1 πt,θ i r t (34) ≤ π ,t * i S 0 + φ θ i Y θ 0 ( 35
)
Replacing θ t= π t i X t 0 in Inequalities (30) by its upper bound provided by (35), we have:

α -1 i E k i + π ,t * i S 0 + φ θ i Y θ 0 ≥ t∈U :t≤τ -1 (d t -L t ) (36)
Finally, we have:

S 0 π ,t * i + α -1 i E k i + k<t≤t * φ t i Y t 0 ≥ S 0 π ,t * i + α -1 i E k i + φ θ i Y θ 0 + φ τ i Y τ 0 ≥ t∈U :t≤τ -1 (d t -L t ) + t∈U :t≥τ (d t -L t ) ≥ t∈U (d t -L t )
This concludes the proof of validity for Inequality (23). The same arguments can be used to prove the validity of Inequalities ( 24) and (25).

It is worth mentioning that the (k, U ) inequalities used in [START_REF] Quezada | A multi-stage stochastic integer programming approach for a multi-echelon lot-sizing problem with returns and lost sales[END_REF] to strengthen the formulation ( 10)-( 22) can be seen as a particular case of the more general family of ( , k, U ) inequalities ( 23 25) is at least as tight as the linear relaxation strengthened by the (k, U ) valid inequalities used in [START_REF] Quezada | A multi-stage stochastic integer programming approach for a multi-echelon lot-sizing problem with returns and lost sales[END_REF].

Proof. Let P LR be the linear relaxation of polyhedron given by inequalities ( 11)-( 22), ( 23)-( 25) and PLR be the linear relaxation of polyhedron given by inequalities ( 11)-( 22) and the (k, U ) inequalities. As any (k, U ) inequality is a valid inequality ( 23)-( 25) with φ t = ν∈U :t≤ν d ν and = 0, we have P LR ⊆ PLR .

The main implication of Proposition 2 is that the lower bound obtained by strengthening the formulation ( 10)-( 22) with the ( , k, U ) inequalities is at least as tight as the lower bound obtained while using the single-echelon (k, U ) inequalities.

We now briefly discuss the resolution of the separation problem for the ( , k, U ) valid inequalities. Recall that this problem consists in finding an inequality ( 23)-(25) violated by a given solution ( X, Ỹ , S, Ẽ, L, Q) of the linear relaxation of Problem ( 11)-( 22) or prove that no such inequality exists. In the present case, in order to find the most violated inequality among e.g. inequalities (23), we should find, for each period k = 1, ..., T , the set U of time periods and the period that maximize the difference between the right-hand and the left-hand side of the inequality. This is not trivial, in particular because the value of each coefficient φ t i simultaneously depends on U and . We thus consider a heuristic separation algorithm in our computational experiments. This one can be summarized as follows. For a given process p and time period k:

1. For each period t = k + 1, ..., T , add t to U if d t (1 - t τ =k+1 Ỹ τ p ) -Lt > 0. 2.
For each period = 0, ..., k, -compute the value of each coefficient φ t i = min <ν≤t r ν πν,t i , ν∈U :t≤ν d ν -compute the left-hand side of the inequality (23) (resp. ( 24) and ( 25)). 3. Set to the period index which minimizes this left-hand side value. This algorithm has a time complexity of O(T 2 ) as the computation of set U in step 1 and of coefficients φ i in step 2 both require O(T 2 ) operations.

Computational Experiments

In this section, we focus on assessing the performance of the proposed valid inequalities when used within a customized branch-and-cut algorithm. We compare the performance of this algorithm with the one of the generic branch-and-cut algorithm embedded in a mathematical programming solver and the one of a branch-and-cut algorithm using single-echelon (k, U ) inequalities.

Instance Generation

We considered two sets of instances: Set 1 instances involve T = 25 periods and I = 10 parts whereas Set 2 instances involve T = 35 periods and I = 10 parts. Within each set, the instances were randomly generated by adapting the procedure presented in [START_REF] Guan | Cutting planes for multistage stochastic integer programs[END_REF]. More precisely, we considered four values of the setupholding cost ratio f /h ∈ {600, 1200, 1800, 2400}, two values for the productionholding cost ratio g/h ∈ {2, 4} and three values of the returns-demand quantity ratio r/d ∈ {1, 2, 3}. For each set and each possible combination of f /h, g/h, r/d, ten random instances were generated, resulting in a total of 480 instances.

For each instance, the value of each problem parameter was set as follows. -Demand d t was uniformly distributed in the interval [0, 100] and the returns quantity r t was uniformly distributed in the interval -The bill-of-materials coefficients α i = α i+I , i ∈ I r , were randomly generated following a discrete uniform distribution over [1; 6] and we set α 0 = α 2I+1 = 1.

-The holding cost h t 0 for the returned product i = 0 was fixed to 1. The holding cost h t i for each recoverable item i ∈ I r was randomly generated following a discrete uniform distribution over interval [START_REF] Ali | A theoretical and computational analysis of lot-sizing in remanufacturing with separate setups[END_REF][START_REF] Leung | Facets and algorithms for capacitated lot sizing[END_REF]. Similarly, the holding cost h t i for each serviceable item i ∈ I s was randomly generated following a discrete uniform distribution over interval [START_REF] Leung | Facets and algorithms for capacitated lot sizing[END_REF][START_REF] Quezada | A multi-stage stochastic integer programming approach for a multi-echelon lot-sizing problem with returns and lost sales[END_REF]. Finally, in order to ensure non negative echelon costs, we set the value of the inventory holding cost for the remanufactured product, h t 2I+1 , to The unit penalty cost for lost sales, l n , was fixed to 10000 per

Results

We carried out extensive numerical experiments in order to assess the computational performance of the proposed valid inequalities. This was achieved by solving each considered instance using three alternatives branch-and-cut algorithms: -CP X: the generic branch-and-cut algorithm embedded in CPLEX 12.8 using the echelon-stock formulation ( 11)-( 22). -(k, U ): a customized branch-and-cut algorithm using the (k, U ) inequalities to strengthen the echelon-stock formulation ( 11)-( 22) similarly to what was done in [START_REF] Quezada | A multi-stage stochastic integer programming approach for a multi-echelon lot-sizing problem with returns and lost sales[END_REF]. -( , k, U ): a customized branch-and-cut algorithm using the newly introduced ( , k, U ) inequalities to strengthen the echelon-stock formulation (11)-( 22). This algorithm is based on the solver CPLEX 12.8. It generates inequalities of type (23)-(25) through a cutting-plane generation algorithm at the root node and at intermediate nodes of the branch-and-bound search tree using the UserConstraints callbacks provided by the solver.

All related linear programs and mixed-integer linear programs were solved using CPLEX 12.8 with the solver default settings. The algorithms were implemented in C++ using the Concert Technology environment. All tests were run on the computing infrastructure of the Laboratoire d'Informatique de Paris VI (LIP6), which consists of a cluster of Intel Xeon Processors X5690. We set the cluster to use two 3.46GHz cores and 12GB RAM to solve each instance. We imposed a time limit of 3600 seconds.

The corresponding results are displayed in Tables 1 for Set 1 instances and 2 for Set 2 instances. Column Method indicates the branch-and-cut algorithm used to solve the instances. Column R.LP gap reports the gap between the value of the linear relaxation strengthened by the corresponding valid inequalities and the best feasible solution found through the branch-and-bound search. For the CPX method, it reports the gap between the value of the initial linear relaxation and the best feasible solution found through the branch-and-bound search. Column R.MIP gap reports the gap between the lower bound at the root node (after the In general, we observe that the customized branch-and-cut algorithms based either on the (k, U ) or on the ( , k, U ) inequalities outperform method CPX, providing solutions of better quality within shorter computation times. Specifically, the total computation time is reduced on average by 20% when using the branch-and-cut algorithm based on the (k, U ) and by 26% when using the branch-and-cut algorithm based on the ( , k, U ) inequalities.

Regarding the relative performance of the (k, U ) and ( , k, U ) inequalities, we note that the branch-and-cut algorithm based on the ( , k, U ) outperforms the algorithm based on (k, U ) when the value of the demand-returns ratio is small, ie. when r/d ∈ {1, 2}. Thus, over the 160 instances corresponding to a value of r/d equal to 1, the total average computation time is reduced from 1732s when using (k, U ) inequalities to 1482s when using ( , k, U ) inequalities. Similarly, over the 160 instances corresponding to a value of r/d equal to 2, the average MIP gap is reduced from 2.72% when using (k, U ) inequalities to 2.38% when using ( , k, U ) inequalities. We note however that the relative performance of the proposed ( , k, U ) inequalities deteriorates for the instances corresponding to the largest considered value of the demand-returns ratio. Namely, when r/d is set to 3, the branch-and-cut algorithm based on the (k, U ) inequalities provides a smaller MIP gap and/or smaller computation times. This might be explained by the fact that the corresponding instances involve a large amount of returned products so that the quantity processed on a resource at a given period is not (or at least to a lesser extent) limited by the availability of the returned products. This means that the proposed refinements in the expression of the valid inequalities are less relevant in this case.

Conclusion and perspectives

We considered a lot-sizing problem aiming at planning production for a multiechelon remanufacturing system. This problem can be formulated as a mixedinteger linear program. We focused on strengthening this formulation in order to be able to provide optimal or near-optimal solutions of this problem. Our main contribution is the development of a new set of valid inequalities which take into account, at each production echelon, the limitations on the produced quantities coming from the limited availability of the returned products. The results of our computational experiments show that a branch-and-cut algorithm based on these new valid inequalities performs well as compared to the generic branch-andcut algorithm of CPLEX solver and to a previously published branch-and-cut algorithm based on less general valid inequalities. A first possible research direction could be to develop an exact separation algorithm for the ( , k, U ) inequalities in the presented branch-and-cut algorithm as this may further improve their performance when used in a branch-and-cut algorithm. Moreover, it would also be worth studying whether valid inequalities previously proposed for capacitated lot-sizing problems (see e.g. [START_REF] Bansal | Facets for single module and multi-module capacitated lot-sizing problems without backlogging[END_REF], [START_REF] Leung | Facets and algorithms for capacitated lot sizing[END_REF], [START_REF] Miller | On the capacitated lot-sizing and continuous 0-1 knapsack polyhedra[END_REF]) might be useful to help solving our problem. Additional computational experiments are also needed to assess the size of the largest instances that may be solved with the proposed exact solution approach.

On a longer perspective, we could seek to extend the proposed valid inequalities to lot-sizing problems with returns involving complicating features such as a limited capacity, backlogging, safety stocks or minimum production levels. Finally, extending the proposed valid inequalities and the cutting-plane generation to solve the stochastic version of the problem studied in [START_REF] Quezada | A multi-stage stochastic integer programming approach for a multi-echelon lot-sizing problem with returns and lost sales[END_REF] is also worth investigating.

Fig. 1 :

 1 Fig. 1: Studied remanufacturing system

  )-(25) proposed in this work. Namely, by setting to 0 and by computing the value of φ t without taking the returns into account (i.e. by setting φ t to ν∈U :t≤ν d ν ), each ( , k, U ) inequality (23)-(25) becomes a (k, U ) inequality. Proposition 2. The linear relaxation of formulation (10)-(22) strengthened by valid inequalities (23)-(

  [0.8(r/d) d, 1.2(r/d) d], where d = d t T is the average demand per period. -The proportion of recoverable parts π t i was uniformly distributed in the interval [0.4, 0.6].

Table 1 :

 1 Performance of CPLEX and branch-and-cut methods over instance in Set 1. generation of CPLEX generic cutting planes) and the best feasible solution found through the branch-and-bound search. Column MIP gap reports the gap between the best lower and the best feasible solution found through the branch-andbound search. The average CPU time for the cutting-plane generation of each method is reported in column C.Time, the CPU time spent at the root node in Column R.Time and the average total CPU time in Column T.Time. Note that each line corresponds to the average value of the corresponding 40 instances.

	r/d g/h Method R.LPgap R.MIPgap MIPgap C.Time R.Time T.Time
	1 2 CP X	8.26	4.10	0.06	0.04	1.16 898.00
	(k, U )	6.08	3.66	0.05	0.92	2.93 1,011.79
	( , k, U ) 4.17	3.55	0.03	0.97	2.35 792.64
	4 CP X	5.38	2.49	0.06	0.03	0.97 1,360.48
	(k, U )	3.95	2.33	0.05	0.50	1.47 1,088.78
	( , k, U ) 2.84	2.22	0.04	0.85	1.98 1,056.29
	2 2 CP X	40.86	11.28	0.97	0.03	1.81 2,570.62
	(k, U )	19.36	9.17	0.46	0.35	2.74 1,883.71
	( , k, U ) 11.88	8.68	0.22	0.74	2.67 1,558.09
	4 CP X	37.79	10.71	2.29	0.03	1.59 3,280.23
	(k, U )	18.39	9.25	1.06	0.19	1.47 2,856.09
	( , k, U ) 11.89	8.48	0.94	0.62	2.27 2,883.50
	3 2 CP X	41.86	13.47	0.04	0.03	1.38 545.15
	(k, U )	20.60	9.68	0.02	0.29	1.96 296.65
	( , k, U ) 14.16	9.11	0.02	0.48	1.91 250.52
	4 CP X	38.29	12.44	0.09	0.02	0.98 817.93
	(k, U )	20.81	9.71	0.02	0.15	1.00 458.41
	( , k, U ) 15.36	8.99	0.06	0.41	1.57 522.70

Table 2 :

 2 Performance of CPLEX and branch-and-cut methods over instance in Set 2.

	r/d g/h Method R.LPgap R.MIPgap MIPgap C.Time R.Time T.Time
	1 2 CP X	7.48	3.65	0.25	0.06	2.26 2,544.42
	(k, U )	5.09	3.08	0.25	8.08	13.32 2,300.22
	( , k, U ) 3.49	2.95	0.11	7.59	11.07 1,764.24
	4 CP X	4.96	2.33	0.19	0.06	1.91 2,930.72
	(k, U )	3.31	2.00	0.11	3.66	6.58 2,529.28
	( , k, U ) 2.40	1.91	0.10	5.09	7.95 2,319.77
	2 2 CP X	44.73	12.89	6.76	0.06	4.43 3,599.07
	(k, U )	20.10	10.00	4.34	0.94	5.71 3,508.81
	( , k, U ) 12.85	9.54	3.77	2.09	6.08 3,362.24
	4 CP X	43.15	12.54	7.29	0.06	3.57 3,599.07
	(k, U )	20.02	10.40	5.03	0.53	3.08 3,599.54
	( , k, U ) 13.14	9.49	4.61	2.42	6.01 3,599.44
	3 2 CP X	43.25	14.11	2.67	0.04	2.93 3,044.18
	(k, U )	18.97	9.48	0.55	0.87	4.26 1,855.28
	( , k, U ) 13.38	9.11	0.64	1.40	4.30 1,614.35
	4 CP X	37.56	11.83	1.60	0.04	2.55 2,802.69
	(k, U )	16.82	8.30	0.19	0.49	2.30 1,035.84
	( , k, U ) 12.67	7.86	0.20	1.21	3.67 1,029.84

This work was partially funded by the National Agency for Research and Development (ANID) / Scholarship Program / DOCTORADO BECAS CHILE/2018 -72190160.


Marie-Cécile Fournier

Yohann Foucher

Paul Blanche

Fanny Buron

Magali Giral

Etienne Dantan
email: etienne.dantan@univ-nantes.fr

A joint model for longitudinal and time-to-event data to better assess the specific role of donor and recipient factors on long-term kidney transplantation outcomes

Keywords: joint modeling, time-to-event data, repeated measurements, serum creatinine, graft failure, kidney transplantation

Background In renal transplantation, serum creatinine (SCr) is the main biomarker routinely measured to assess patient's health, with chronic increases being strongly associated with longterm graft failure risk (death with a functioning graft or return to dialysis). Joint modeling may be useful to identify the specific role of risk factors on chronic evolution of kidney transplant recipients: some can be related to the SCr evolution, finally leading to graft failure, whereas others can be associated with graft failure without any modification of SCr.

Methods Sample data for 2749 patients transplanted between 2000 and 2013 with a functioning kidney at 1-year post-transplantation were obtained from the DIVAT cohort. A shared random effect joint model for longitudinal SCr values and time to graft failure was performed.

Results

We show that graft failure risk depended on both the current value and slope of the SCr.

Deceased donor graft patient seemed to have a higher SCr increase, similar to patient with diabetes history, while no significant association of these two features with graft failure risk was found. Patient with a second graft was at higher risk of graft failure, independent of changes in SCr values. Anti-HLA immunization was associated with both processes simultaneously.

Conclusion Joint models for repeated and time-to-event data bring new opportunities to improve the epidemiological knowledge of chronic diseases. For instance in renal transplantation, several features should receive additional attention as we demonstrated their correlation with graft failure risk was independent of the SCr evolution.

INTRODUCTION

Recently, Asar et al. [START_REF] Asar | Joint modelling of repeated measurement and time-to-event data: an introductory tutorial[END_REF] pointed out that joint models for repeated and time-to-event data should be used to correctly consider longitudinal and survival processes and their relationship when they are strongly dependent [START_REF] Rizopoulos | Joint Models for Longitudinal and Time-to-Event Data: With Applications in R[END_REF][START_REF] Wulfsohn | A joint model for survival and longitudinal data measured with error[END_REF]. Joint models could allow identification of specific associations for each risk factor: 1) association with an event risk, 2) association with an event risk throughout a previous marker evolution modification, and 3) associations on both longitudinal and survival processes. With this type of suitable model, the precise description of specific associations could provide new insights in the knowledge of a disease pathway mechanism. Indeed, from a methodological standpoint, longitudinal measurements and time-to-event data are typically analyzed separately, using a mixed model and survival model respectively, without considering their possible relationship [START_REF] Rizopoulos | Joint Models for Longitudinal and Time-to-Event Data: With Applications in R[END_REF]. However, it is known that these two processes are often mutually dependent in a chronic disease context. Inferences from mixed models may be biased in cases of an informative censoring process [START_REF] Leffondre | Analysis of risk factors associated with renal function trajectory over time: a comparison of different statistical approaches[END_REF][START_REF] Ibrahim | Missing data methods in longitudinal studies: a review[END_REF][START_REF] Tsiatis | joint modeling of longitudinal and time-to-event data: an overview[END_REF]. Similarly, the time-dependent Cox model fails to correctly handle a time-dependent endogenous variable, which is a variable generated by the patient themselves (e.g. creatinine) in contrast to an exogenous variable which is not (e.g. air pollution level) [START_REF] Rizopoulos | Tools & techniques--statistics: Dealing with time-varying covariates in survival analysis--joint models versus Cox models[END_REF][START_REF] Andrinopoulou | An introduction to mixed models and joint modeling: analysis of valve function over time[END_REF][START_REF] Kalbfleisch | The Statistical Analysis of Failure Time Data[END_REF]. It also often fails to correctly account for measurement error.

In many chronic diseases, the occurrence of major events and the assessment of corresponding risk factors guide physicians in implementing the most beneficial care for patients. In renal transplantation, serum creatinine (SCr) is a well-known longitudinal marker used to assess the health of kidney transplant recipients especially after the first year post-transplantation, i.e. in the chronic phase of the disease evolution [START_REF] Levey | National Kidney Foundation practice guidelines for chronic kidney disease: evaluation, classification, and stratification[END_REF]. In the chronic phase, graft failure, defined as return to dialysis or death with a functioning graft, is often preceded by a continuous deterioration in renal function and is associated with an irreversible increase in SCr levels. In contrast, during the first year post-transplantation, patients are submitted to a risk of early graft failure due to acute clinical events such as delayed graft function, acute rejection episode, and infections or complications. Therefore the disease evolution mechanism is very different between the acute and chronic phases [START_REF] Galichon | Clinical and histological predictors of long-term kidney graft survival[END_REF][START_REF] Foucher | A clinical scoring system highly predictive of long-term kidney graft survival[END_REF]. Few authors have studied the association of risk factors on renal function evolution [START_REF] Ferro | Bayesian analysis of glomerular filtration rate trajectories in kidney transplant recipients: a pilot study[END_REF][START_REF] Marcén | Long-term graft function changes in kidney transplant recipients[END_REF] while the risk factors associated with long-term graft failure have been well described [START_REF] Foucher | A clinical scoring system highly predictive of long-term kidney graft survival[END_REF][START_REF] Pascual | Strategies to improve long-term outcomes after renal transplantation[END_REF][START_REF] Debout | Each additional hour of cold ischemia time significantly increases the risk of graft failure and mortality following renal transplantation[END_REF]. Renal function has been shown to be one of the most important risk factors [START_REF] De Bruijne | Predicting kidney graft failure using time-dependent renal function covariates[END_REF][START_REF] Kasiske | A thirty percent chronic decline in inverse serum creatinine is an excellent predictor of late renal allograft failure[END_REF]. However, the precise mechanism is not well known: the majority of risk factors leading to graft failure may be associated with chronic SCr changes, finally leading to graft failure, but one can also hypothesize that some features may be related to graft failure risk independently from their association on SCr evolution.

Whilst several authors have used joint models for longitudinal and time-to-event data in renal transplantation [START_REF] Abdi | Impact of longitudinal exposure to mycophenolic acid on acute rejection in renal-transplant recipients using a joint modeling approach[END_REF][START_REF] Moranne | Rate of renal graft function decline after one year is a strong predictor of all-cause mortality[END_REF][START_REF] Rizopoulos | A Bayesian semiparametric multivariate joint model for multiple longitudinal outcomes and a time-to-event[END_REF][START_REF] Garre | A joint latent class changepoint model to improve the prediction of time to graft failure[END_REF], none have precisely studied the specific association of each baseline explicative variable. The precise study of such associations would be of primary importance for physicians in order to improve their appraisal of kidney transplant recipients' health. Therefore, we present for the first time a shared random effect multivariable joint model to study the baseline characteristics that could be related to long-term kidney graft outcomes. By studying SCr evolution and graft failure from one year post-transplantation, it brings an epidemiological approach to understand the risk factors associated with the disease evolution in its chronic phase.

MATERIALS AND METHODS

Study population

Data were extracted from the French observational and prospective DIVAT cohort

Available data

Most classical risk factors susceptible to influence SCr evolution and/or graft failure risk were extracted from the database. Donor features included: age, gender, last SCr level, deceased (from cardiovascular cause vs other) / living donation. Recipient characteristics were: age, gender, Body Mass Index (BMI), history of comorbidities (diabetes, hypertension, dyslipidemia, neoplasia, cardiovascular), duration of dialysis before transplantation, preemptive graft, hemodialysis or peritoneal dialysis, pre-transplant anti-class I or class II Human Leucocyte Antigen (HLA) immunization and cause of initial renal disease (recurrent nephropathy or not).

Transplantation parameters were: cold ischemia time and number of HLA-A-B-DR incompatibilities. The following variables were collected within the first year posttransplantation: occurrence of delayed graft function (defined as the need for dialysis after transplantation), occurrence of acute rejection episodes and SCr levels at 3 and 6 months posttransplantation.

Outcomes

The baseline was the 1-year post-transplantation anniversary. The best marker of renal function should be the measured GFR (mGFR) [START_REF] Stevens | Assessing kidney function--measured and estimated glomerular filtration rate[END_REF]. Unfortunately, this measurement is costly and is not performed in practice for routine patient follow-up. Different equations to estimate GFR (eGFR) have been proposed [START_REF] Levey | A new equation to estimate glomerular filtration rate[END_REF][START_REF] Levey | A more accurate method to estimate glomerular filtration rate from serum creatinine: a new prediction equation. Modification of Diet in Renal Disease Study Group[END_REF][START_REF] Nankivell | Predicting glomerular filtration rate after kidney transplantation[END_REF][START_REF] Cockcroft | Prediction of creatinine clearance from serum creatinine[END_REF], and are principally based on SCr values adjusted on recipient age, gender and ethnicity. A major limitation lies in the fact that these eGFR formulae have been developed from general population data or from patients with chronic kidney disease and thus may not be applicable to kidney transplant patients [START_REF] Buron | Estimating glomerular filtration rate in kidney transplant recipients: performance over time of four creatinine-based formulas[END_REF][START_REF] White | Estimating Glomerular Filtration Rate in Kidney Transplantation: Is the New Chronic Kidney Disease Epidemiology Collaboration Equation Any Better?[END_REF][START_REF] Gaspari | Performance of different prediction equations for estimating renal function in kidney transplantation[END_REF]. Despite this, SCr or eGFR are equivalent in terms of relative evolution, the second resulting from a transformation of the first parameter. We decided to study the SCr (µmol/L), which was the longitudinal marker routinely recorded yearly until patient death with a functioning graft or return to dialysis. The time-toevent (graft failure) was defined as the delay between 1-year post-transplantation and the first event between return to dialysis or death with a functioning graft.

Statistical analysis

We used a shared random effect model. It combines a mixed model for the longitudinal process and a parametric survival model for the time-to-event process, for which underlying hypotheses were checked in an independent manner. For the longitudinal assessment, we used a logarithmic transformation of SCr values in order to respect both assumptions related to residual's homoscedasticity and linear relationship over time. Two random effects were considered for the baseline value and the slope. For survival, hazard proportionality and log-linearity were assessed.

In a joint shared random effect framework, longitudinal and survival processes are linked through common random effects. Rizopoulos has previously described the possible parameterizations to model this dependence [START_REF] Rizopoulos | Joint Models for Longitudinal and Time-to-Event Data: With Applications in R[END_REF]. For instance, the survival process can be modeled as dependent on the current level of the marker, on the intensity of marker evolution (i.e. the slope), on both current level and slope, on cumulative effects or on lagged effects.

In the first step of model building, we defined a baseline risk function and the dependence between the two processes from a joint model without baseline explicative variables. We graphically retained a Weibull distribution for the baseline risk function. According to the Bayesian Information Criteria (BIC), the dependence between the two processes was characterized by the instantaneous hazard of graft failure depending on both the level and the slope of the longitudinal marker at the current time. As recommended by Rizopoulos to solve optimization difficulties, all quantitative variables were standardized in order to scale the coefficients (2).

In the second step, we performed the selection of baseline explicative variables. Univariable models were composed using three effects of each variable: on baseline value, on the slope (interaction with time) and on the graft failure risk. Among these parameters, those which were not significant (p>0.05) were removed in a hierarchical manner: if the association on the slope was significant, the corresponding association on baseline value was also considered. Finally, a multivariable joint model was generated by including effects retained in the univariable models, and a forward stepwise selection was performed (always using a 5% type-I error rate).

In order to study the relevance of the joint modeling, we also performed separate analyses: i) a linear mixed model to study the SCr evolution and ii) a time-dependent Cox model to study the graft failure risk. We used the same variables selection procedure.

As sensitivity analyses, we performed two joint models in a cause specific approach: i) time-toreturn to dialysis by censoring death, and in contrast ii) time-to-death with a functioning graft by censoring return to dialysis.

Joint model parameters were estimated by likelihood maximization. The complete mathematical formulation of the joint model is shown in appendix 1. Due to the logarithmic transformation of SCr, coefficients for the longitudinal process have an interpretation as relative change rather than absolute change. Details related to interpretations are presented in appendix 2. Confidence intervals for relative change were obtained using parametric simulations (5 000 iterations) [START_REF] Aalen | A Markov model for HIV disease progression including the effect of HIV diagnosis and treatment: application to AIDS prediction in England and Wales[END_REF].

All analyses were performed using the 3.0.1 version of the R software (32) with the 1.3-0 version of the JM package [START_REF] Rizopoulos | An R Package for the Joint Modelling of Longitudinal and Time-to-Event Data[END_REF].

RESULTS

Baseline characteristics

Baseline characteristics are presented in Table 1. Sixty percent of the recipients were male, with a mean age of 49.7 ± 13.6 years. Histories of cardiovascular disease or dyslipidemia were observed in one third of recipients, 11.6% had history of diabetes, 82.6% had hypertension, and 8.3% had a cancer before the transplantation. Second transplantations were realized in 17.2% of studied patients. Immunologic characteristics included: 12.8% of patients presented more than 4 HLA-A-B-DR incompatibilities, and around one third were immunized against class I or class II HLA prior to transplantation. Donors were mainly deceased (84.8%) with a mean age of 50.7 ± 15.5 years and 56.4% were male. Delayed graft function occurred for 714 patients (26.1%). SCr at 3 and 6 months were on average 138.3 ± 53.4 and 136.6 ± 53.2 µmol/L respectively. Finally, 21.5% of recipients presented at least one episode of acute rejection before the first anniversary of the graft.

Follow-up description

During follow-up, 278 patients returned to dialysis and 203 died with a functioning graft. The median event-free follow-up time was 3.99 years. The patient-graft survival curve and its corresponding 95% confidence interval (95%CI) are presented in Figure 1. Patient-graft survival rates at 10 years after the first anniversary of the graft was 58.4% [95%CI: 54.8% ; 62.3%].

Additionally, 12 843 SCr measurements were collected, with a median of 4 measurements per patient (ranging from 1 to 14). The median time between two measurements was 11.7 months (interquartile range: 9.2, 12.5).

Joint modeling

Table 2 presents the estimations related to the final multivariable joint model.

Dependence between SCr dynamic and graft failure risk

For any time one year after transplantation (t>1), the graft failure risk depended on both the current value and the current slope of the SCr. If a patient had a 25% higher SCr, graft failure risk was twice as high (HR=1.92, 95%CI: [1.75 ; 2.11]). Moreover, for a given SCr value, where a patient had a steeper increase in SCr graft failure risk was significantly worse (HR=1.89, 95%CI: [1.17 ; 3.06] for an increase of 25% in SCr value in one year).

Factors associated with 1-year post-transplantation SCr

Several factors appeared significantly correlated with a higher 1-year SCr without significant association with the SCr evolution or with graft failure risk. An increase of 50 µmol/L in the 6month SCr level was associated with a 1-year SCr increase of 17.99% (95%CI: [16.62% ; 19.34%]). Patient with a graft provided from a donor ten years older compared to other donor had a 3.68% higher SCr at 1 year (95%CI: [2.97% ; 4.39%]).

Factors associated with SCr evolution during follow-up

A history of diabetes was associated with a higher SCr increase. After 5 years, the presence of this comorbidity for a patient was associated with a SCr 14.45% higher (95%CI: [7.76% ; 21.46%]) compared to its absence, while we did not observe any significant difference at 1 year. For a given patient, a graft from a deceased donor due to a cerebrovascular cause was associated with a relative increase of 12.52% in expected 5-year SCr (95%CI: [6.50% ; 18.89%]), compared to a graft from a living donor, but only a 7.14% increase when compared to a deceased donor due to another cause (95%CI: [1.48% ; 13.00%]).

Factors associated with both 1-year post-transplantation SCr and SCr evolution during the follow-up

Male recipient had a 7.68% higher 1-year SCr (95%CI: [5.83% ; 9.51%]). Nevertheless, he was likely to have a slower increase: after 5 years follow-up the gap reduced to 3.98% (95%CI: [0.18% ; 7.81%]).

Factors associated with the graft failure risk

Two factors were significantly associated with graft failure risk. Without significant correlation with the SCr, a patient transplanted for the second time had a higher graft failure risk compared to his first transplantation (HR=1.32, 95%CI: [1.02 ; 1.73]). On the contrary, for a given patient, graft from male tended to be associated with a smaller risk of graft failure compared to graft from female donor (HR=0.83, 95%CI: [0.69 ; 1.01]).

Factors associated with both SCr and the graft failure processes

A 10-year older patient was associated with a 2.04% lower 1-year SCr (95%CI: [1.31% ; 2.77%]). Moreover, this difference increased during the follow-up, ie 5.57% lower at 5years post-transplantation (95%CI: [4.20% ; 6.95%]). This could be explained by lower creatinine production in the elderly population due to a smaller muscle mass. However, for a given SCr evolution during the follow-up, the situation where the recipient is older was associated with a higher graft failure risk (HR=1.35 for a 10-year older patient, 95%CI: [1.25 ; 1.46]). For a given patient, higher 3-month SCr was associated with a significantly higher 1-year SCr level. For instance, for a 50 µmol/L difference at 3-months, the 1-year SCr level increased by 8.08% (95%CI: [6.83 ; 9.32]). Additionally, for a given evolution of SCr from 1-year posttransplantation, a higher 3-month SCr was associated with a lower graft failure risk (HR=0.85 for an increase of 50 µmol/L in 3-month SCr level, 95%CI: [0.75 ; 0.95]).

Major risk factors included acute rejection in the first year post-transplantation, immunization, and cardiovascular history. When a patient had an acute rejection episode during the first year post-transplantation, a 5.65% higher 1-year SCr was observed compared to cases where no acute rejection occurred (95%CI: [3.65 ; 7.71]). Nevertheless, independently of the current value and the slope of SCr, the situation where acute rejection has occurred appeared with a higher risk of graft failure (HR=1.46, 95%CI: [1.17 ; 1.83]). A patient with cardiovascular history was more likely to have an increased SCr compared to the same patient without this history, and a higher graft failure risk independently from this increase (HR=1.39, 95%CI: [1.14 ; 1.69]). Similarly, a significant SCr increase was demonstrated for pre-transplant immunized patient, with an additional graft failure risk not related to this SCr increase, compared to the same patient nonimmunized.

Separate models

The linear mixed model estimations and those of the time-dependent Cox model are presented in tables 3 and 4 respectively. One can note differences in the retained variables. An acute rejection episode in the first year post-transplantation was significantly associated with the SCr evolution by using a linear mixed model, in opposition to the results obtained by using joint models. No relationship between cardiovascular history or donor type with SCr evolution were retained by using the linear mixed model, while we concluded from the final joint model that patients with cardiovascular history may have a significantly higher SCr increase during the follow-up (p= 0.0371) and SCr evolution could be different given the donor type status (p=0.0022). Slight underestimations of hazard ratios were obtained from the time-dependent Cox model compared to the joint model. For example, the hazard ratio related to the recipient age was 1.25 (95%CI [1.19 -1.38]) for 10 years older by using the time-dependent Cox model against 1.35 (95%CI [1.25 ; 1.46]) by using the joint model. Additionally, diabetes was retained as a risk factor for graft failure by using the Cox model, while it does not by using the joint model. In contrast, acute rejection episode was not retained as a risk factor for graft failure by using the time-dependent Cox model, while it was by using the joint model.

Cause specific approach

Using a cause specific approach (tables S1 and S2 in supplementary materials), we observed that current SCr level was more importantly associated with return to dialysis (HR=2.51, 95%CI [2.22 ; 2.84]) compared to death with a functioning graft (HR=1.47, 95%CI [1.24 ; 1.74]). As expected, this higher association was also observed for acute rejection episode (HR=1.63, 95%CI [1.20 ; 2.20]

DISCUSSION

Our results show that during the chronic phase of renal transplantation, elevated SCr levels as well as the magnitude of SCr increases are associated with a higher risk of graft failure.

Accordingly, physicians routinely supervise both the current SCr level and its increase. The large majority of baseline explicative variables are firstly associated with the baseline SCr level or its evolution, finally leading to graft failure. Interestingly, we demonstrated that besides the association of cardiovascular history with increased SCr, this risk factor was additionally associated with an increase in the risk of graft failure. Therefore, at a given time for a given SCr level and slope, the presence of cardiovascular history should be considered as a risk factor for graft failure. Similarly, patient transplanted for a second time seemed at higher risk of graft failure, regardless of the SCr level or its slope, compared to its first graft. Other factors independent of SCr leading to increased graft failure risk may result from stronger immunosuppression or undetected immunization against donor specific antigens. In addition to retransplantation or the presence of cardiovascular history, particular attention should also be paid to patients with a high 3-month SCr level, transplantation in older patients or when an acute rejection episode during the first year has occurred. These patients may be more susceptible to graft failure without having previously displayed aberrant SCr levels. On the contrary, if a patient received a deceased donor graft but had a SCr evolution analogous to those which would be observed if the graft had come from a living donor, the monitoring of this patient should be the same regardless of the donor status.

In renal transplantation, numerous studies have focused on only one or two measurements of renal function to study their association with graft failure [START_REF] Kasiske | A simple tool to predict outcomes after kidney transplant[END_REF][START_REF] Lenihan | MDRD-estimated GFR at one year post-renal transplant is a predictor of long-term graft function[END_REF][START_REF] Hariharan | Posttransplant renal function in the first year predicts long-term kidney transplant survival[END_REF]. However, the joint modeling approach allows the whole trajectory of longitudinal SCr measurements to be taken into account.

In this paper, we used for the first time a shared random effect joint modeling to more precisely specify the association between chronic SCr evolution and graft failure risk. Different types of dependence can be considered such as the current marker level, the evolution intensity during the follow-up, cumulative effects or lagged effects [START_REF] Rizopoulos | Joint Models for Longitudinal and Time-to-Event Data: With Applications in R[END_REF].

More generally, our approach illustrates that joint modeling constitutes a powerful approach for time-to-event analysis with endogenous time-dependent variable [START_REF] Asar | Joint modelling of repeated measurement and time-to-event data: an introductory tutorial[END_REF], which supports a real mechanistic evolution for many chronic diseases. However, their use in observational studies is still uncommon. As previously acknowledged by Asar et al. [START_REF] Asar | Joint modelling of repeated measurement and time-to-event data: an introductory tutorial[END_REF], differing results and interpretations between the joint modeling and the separate approaches reinforce the necessity to use joint modeling in the presence of endogenous variable. We also highlighted the differences in our application in kidney transplantation. In other diseases, the informative censoring or the endogenous nature for the longitudinal variable can result in even higher differences.

One limitation in our study may be the graft failure definition: the first event between the return to dialysis and death with a functioning graft. Because it is very difficult to distinguish the cause of death related or not to the disease, we performed a sensitivity analyses. The results illustrated the overall robustness of the results but with a probable underestimation of the association between the SCr and the acute rejection episode on the risk of graft failure.

In conclusion, our results illustrate the importance of joint models and their potential usefulness in improving chronic disease research. It brings a more complete epidemiological view of the risk factors and the related natural disease history mechanisms. The use of this novel statistical model on a large cohort of kidney transplant recipients highlights that several risk factors were associated with SCr evolution while others were associated with graft failure risk independently of the initial SCr value or its subsequent evolution. These included elderly or immunized recipients, second transplantations, grafts coming from female donors, patients experiencing an acute rejection episode in the first year post-transplantation, patients with cardiovascular history or with a high gap between 3 and 12 month SCr measurements, features that should receive additional attention. = exp (β 2 Δ/ sd W1 )exp(β 3 t ij Δ/ sd W1 ))

Tables and Figures

2.2/ Hazard ratio for the longitudinal marker

As we have seen in appendix 1, the instantaneous risk function is written as follows: hi(t) = h0(t) exp(γ T X3i + α1mi(t) + α2mi'(t))

= h0(t) exp(γ T X3i + α1 mi(t) + α2 As we use a log transformation of SCr measurement (Y(t) = log(SCr(t))), the hazard ratio which quantifies the association between the longitudinal marker and the risk of event was expressed for a clinically relevant difference.

-For the current level of the marker, we can rewrite the HR for a difference of 25% in SCr values at the same time for the same patient and the same slope: HR 1.25SCr(t)vs SCr(t) = h 0 (t)exp (γ T X 3i + α 1 log(1.25SCr(t)) + α 2 m i '(t)) h 0 (t)exp (γ T X 3i + α 1 log(SCr(t)) + α 2 m i '(t))

= exp (α 1 (log(1.25SCr(t)) -log(SCr(t))))

= 1.25 α 1

(

  www.divat.fr) of kidney transplant recipients from 6 University hospitals (French Research Ministry: RC12_0452, last agreement No 13 334, No CNIL for the cohort: 891735, No CNIL for the study: 914226). According to the following inclusion criteria, 2749 patients were studied: adult recipients who received a first or second renal transplant between January 2000 and August 2013 from a living or heart beating deceased donor, alive with a functioning graft at 1year post-transplantation and maintained under Tacrolimus and Mycofenolic acid. All study participants gave informed consent.

Fig. 1 :

 1 Fig. 1: Patient and graft survival according to the time since the first anniversary of the transplantation (n = 2749) from Kaplan-Meier estimator and their corresponding 95% confidence interval

Table 1 : Description of recipients, donors, and transplantation characteristics of the studied population (n = 2749)

 1

	sd: standard deviation ; BMI: Body Mass Index; SCr: Serum Creatinine; HLA Human Leukocyte Antigen

ACKNOWLEDGEMENTS

We wish to thank the DIVAT scientific board (C Legendre, H Kreis, L Rostaing, N Kamar, E Morelon, G Mourad, V Garrigue, M Kessler and M Ladrière) as well as members of the clinical research assistant team (S Le Floch, A Petit, J Posson, C Scellier, V Eschbach, K Zurbonsen, C Dagot, F M'Raiagh, V Godel, X Longy and P Przednowed). The DIVAT cohort is partially supported by Roche Laboratory since 1994.

Conflict of Interest: The authors declare that they have no conflict of interest.

Table 2: Multivariable joint model for longitudinal evolution of logarithmic transformation of serum creatinine (SCr) and risk of graft failure (return to dialysis or death with a functioning graft) in kidney transplant patients (n=2584 patients, 165 patients excluded due to missing data)

Coef: coefficient; HR: Hazard Ratio; CI: confidence interval. Referential value for 1-year SCr was 4.024, 95%CI: [3.982 ; 4.065]. Referential value for SCr evolution was 0.034 95%CI: [0.018 ; 0.050]. This model is adjusted on a time effect with a threshold at 2008 (before 2008 vs after): coefficient for the relation to the SCr at 1-year: 0.018 95%CI: [0.002 ; 0.034] and to the SCr evolution: 0.013 95%CI: [0.005 ; 0.020] and HR = 0.73 [0.57 ; 0.94]. Parameters of the Weibull baseline risk function were: intercept -20.247 ± 0.982 ; log(shape): 0.337 ± 0.046. α= 2.93 ; α2 = 3. [START_REF] White | Estimating Glomerular Filtration Rate in Kidney Transplantation: Is the New Chronic Kidney Disease Epidemiology Collaboration Equation Any Better?[END_REF] Longitudinal process Association with the Association with the Survival process log (

Appendix 1: Mathematical formulation of the shared random effect joint model

Let Y be the longitudinal marker and tij the time of measurement of the j th (j = 1, …, ni) measure for the patient i (i = 1, …, N). Let h(.) denotes the instantaneous risk function of graft failure. The joint model combines a linear mixed model (equation 1) with a parametric regression model (equation 2). They share the random effects

with (b0i ; b1i) T ~MVN(0,B), B an unstructured variance-covariance matrix, X1i a vector of baseline covariates influencing the baseline value of longitudinal marker, X2i another vector of baseline covariates that may change marker evolution over time and β0, β1 two scalars defining the referential value of the baseline level and the slope of the longitudinal biomarker Y (.) respectively, and β2, β3 two p-vectors of the same dimension as X1 and X2 respectively. The evolution of the measurements Yij(tij) are defined by the sum of a subject specific trend mi(tij) plus an error term εij~N(0,σε 2 ). For the instantaneous risk function of graft failure, h0(t) denotes the baseline risk function, and X3i is a vector of baseline covariates that could influence the graft failure risk, with a corresponding vector of fixed regression coefficients γ. g is a function of the true level of the marker mi, which specifies the type of dependence between the longitudinal and the survival processes. Classically, it may be the current level of the marker (g(mi(t)) = αmi(t)), the intensity of marker deterioration during the follow-up i.e. the slope (g(mi(t)) = α2mi'(t)), or both (g(mi(t)) = α1mi(t) + α2mi'(t)) (2). This latter is the retained association of the model presented in table 2. Let Z be a qualitative variable associated with:

-the 1-year SCr only (Z ⊆ X1 ; Z ⊆ X2). The excess of SCr for a patient with Z = 1 as compared to the case where Z = 0 for the same patient is:

This gap of SCr is constant beyond 1-year post-transplantation.

-both the 1-year SCr and the SCr increase during the follow-up (Z ⊆ X1 ; Z ⊆ X2)

This gap of SCr value is increasing or decreasing during the follow-up according to the sign of β3. For clinical purposes, in the interpretations, we used the time t = 5 to quantify a relative change at 5 years after the first year post-transplantation.

2.1.2/ Quantitative variables:

Let W1 be a quantitative variable with sdW1 its standard deviation, w a value of W1 and Δ a relevant clinical increase.

-For the intensity of the marker, the HR which compares the situation in which )

This leads to: HR = exp(α 2 (log(1 + x/100) -log(1 + y/100))) which is the HR which compares an increase of x% between t-1 and t to an increase of y %. In our paper, we choose to compare an increase of 25% compare to the mean evolution (a growth of 3% each year).

2.3/ Hazard ratio for the quantitative variables

Because the quantitative variables have been standardized, the HR for these factors were expressed for an increase of one standard deviation. In order to calculate them for an increase of relevant threshold in the variable unit, we can proceed as follows:

Let X1 be the standardization of W1 with sd1 its standard deviation. HRX is the HR obtained for the standardized variable and HRW is the one for an increase of Δ unit of W1.

𝑊 = HR X ( Δ sd 1 )