Thibaut Neveux
email: thibaut.neveux@edf.fr

Maxime Bretaud

Norinda Chhim

Khalil Shakourzadeh

Sophie Rapenne

Pilot plant experiments and modeling of CaCO 3 growth inhibition by the use of antiscalant polymers in recirculating cooling circuits

Keywords: Cooling water, Precipitation, Antiscalant, Inhibition, Adsorption, Polymer

come

Pilot plant experiments and modeling of CaCO3 growth inhibition by the use of antiscalant polymers in recirculating cooling circuits

Introduction

Precipitation of CaCO 3 is a constraining issue in industrial operations, particularly in energy-intensive industries such as cooling circuits of thermal and nuclear power plants (Liu et al., 2012a;[START_REF] Rahmani | Evaluation of inhibitors and biocides on the corrosion, scaling and biofouling control of carbon steel and copper-nickel alloys in a power plant cooling water system[END_REF], as well as separation processes such as distillation [START_REF] Al-Rawajfeh | Simultaneous desorption-crystallization of CO2-CaCO3 in multi stage flash (MSF) distillers[END_REF][START_REF] Ghani | Impacts of different antiscalant dosing rates and their thermal performance in Multi Stage Flash (MSF) distiller in Kuwait[END_REF], membrane desalination and reverse osmosis [START_REF] Phuntsho | Membrane scaling and flux decline during fertiliser-drawn forward osmosis desalination of brackish groundwater[END_REF][START_REF] Greenlee | The effect of antiscalant addition on calcium carbonate precipitation for a simplified synthetic brackish water reverse osmosis concentrate[END_REF]. In these processes, the precipitation of calcium carbonate occurs by the increase in water temperature and/or salt concentration, attributed to the removal of water from the liquid phase by evaporation or permeation. Hence, CaCO 3 precipitates on surfaces, and thermal or separation efficiency decreases, resulting in increased operative costs for the prevention or control of scaling. In addition, structures are weakened, requiring considerable maintenance, thereby decreasing the equipment lifespan. Among several possible treatments, the addition of soluble polymers (such as antiscalant polymer or ASP) in water is an economic way for controlling the formation of scale as they significantly decrease the precipitation of CaCO 3 at low polymer concentration, typically less than 10 ppm of active matter. A wide range of chemical compounds can be utilized, such as polymers based on acrylates, amines, phosphonates, sulfonates, and aspartates [START_REF] Lin | Inhibition of calcite crystal growth by polyphosphates[END_REF][START_REF] Hasson | State of the art of friendly "green" scale control inhibitors: a review article[END_REF][START_REF] Macadam | Calcium carbonate scale formation and control[END_REF]. As some chemicals alter the environment, e.g., phosphorous compounds, recent studies have reported the development of environment-friendly inhibitors [START_REF] Hasson | State of the art of friendly "green" scale control inhibitors: a review article[END_REF][START_REF] Chaussemier | State of art of natural inhibitors of calcium carbonate scaling. A review article[END_REF].

Several studies have been reported either on the efficiency comparison of antiscalants between several polymers or identification of mechanisms [START_REF] Macadam | Calcium carbonate scale formation and control[END_REF][START_REF] Liu | Research on performance evaluation and anti-scaling mechanism of green scale inhibitors by static and dynamic methods[END_REF][START_REF] Tang | Investigation of CaCO3 scale inhibition by PAA, ATMP and PAPEMP[END_REF][START_REF] Amjad | Scale inhibition in desalination applications-an overview[END_REF][START_REF] Li | Mineral precipitation and deposition in cooling systems using impaired waters: mechanisms, kinetics, and inhibition[END_REF]; however, relatively few studies have been reported on the modeling of inhibition efficiency. Two general inhibition models have been reported in literature according to the inhibition mechanism: nucleation inhibition and growth inhibition. Nucleation inhibition is used so that the induction time (period before the appearance of the first crystalline germ) becomes higher than the process residence time so as to ensure that no scaling occurs in the process. Most previous studies have proposed dose calculations based on this mechanism, i.e., the minimum dose required for preventing germ formation. For example, experience-based orders of magnitude of inhibitor dosages have been suggested by Puckorius (reported by [START_REF] Cavano | Understanding scaling indices and calculating inhibitor dosages[END_REF] with respect to the Langelier, Ryznar, and Puckorius saturation indexes for three inhibitors. [START_REF] He | Inhibition of calcium carbonate precipitation in NaCl brines from 25 to 90 •[END_REF][START_REF] He | Inhibition of mineral scale precipitation by polymers[END_REF] have proposed a model based on nucleation inhibition by measuring the deviations with respect to induction time, linking the variation of surface energies to the inhibitor concentration; they have utilized their model for calculating the inhibitor dose required for injection with respect to the saturation index, temperature, and desired elongated induction time. [START_REF] Fergusson | Modeling scale formation and optimizing scale inhibitor dosages in Membrane Systems[END_REF] has proposed a similar model for calculating the induction time with respect to inhibitor dose for several crystals such as calcite, gypsum, or barite. [START_REF] Vanderpool | Calculating minimum threshold inhibitor dosage[END_REF] has also proposed a procedure for calculating the minimum dose required for several inhibitors, including a dose for the inhibition of homogeneous nucleation, a dose for dispersed matter, and a dose for maintaining equilibrium. On the other hand, some studies have proposed models for describing the kinetic growth inhibitor by polymers, for systems where nucleation has already occurred; hence, polymers are used for decreasing the precipitation kinetics of CaCO 3 , instead of preventing nucleation. Such models typically consider the adsorption of polymers on crystallization growth sites, with the use of adsorption isotherms such as the Langmuir, Freundlich, or Langmuir-Freundlich [START_REF] Reddy | Calcite crystal growth inhibition by phosphonates[END_REF][START_REF] Kubota | A kinetic model for crystal growth from aqueous solution in the presence of impurity[END_REF][START_REF] Kan | Adsorption and precipitation of an aminoalkylphosphonate onto calcite[END_REF][START_REF] Lisitsin | Modeling the effect of anti-scalant on CaCO3 precipitation in continuous flow[END_REF]. The use of either a nucleation inhibition or a crystal growth inhibition model will depend on the system under study. For processes with a low residence time, such as one-through continuous flow (e.g., reverse osmosis and open cooling circuit), the inhibition of nucleation is efficient for preventing scaling. For this purpose, small quantities of polymers are sufficient to make the induction time higher than the residence time. However, for processes with higher residence times, such as in recirculating wet cooling circuits, the pragmatic viewpoint is that nucleation will occur at some point. Hence, the formation of scale is not limited by nucleation kinetics, but by crystal growth kinetics; hence, a crystal growth inhibition model is required.

In this study, emphasis is placed on recirculating cooling water systems, as shown in Figure 1, although the results obtained herein can be applied to other systems with similar operating conditions, especially with respect to water quality, temperature range, and polymer concentration. In wet recirculating cooling circuits, water is heated in an exchanger called condenser, which is then cooled in a cooling tower predominantly by evaporation, and a large amount of water is then recycled (around 10-50 times the make-up flow rate) to the condenser after cooling. Typically, for a power plant, up to several hundreds of metric tons of CaCO 3 deposits are observed annually, corresponding to a significant economic impact [START_REF] Walker | Economic impact of condenser fouling in existing thermoelectric power plants[END_REF]. Owing to the large volume of water in the circuit (implying a water mean residence time of several hours) and a high recirculation ratio, a large amount of chemicals are required to avoid nucleation, such as polymers to extend the induction time and/or acids to decrease the pH of water. However, this would lead to unacceptably high treatment costs and environmental impact, and hence, it cannot be performed. As a result, nucleation eventually occurs during process operation, even for initially cleaned systems; hence, polymers are added for limiting the growth kinetic rate. On the other hand, the models proposed for the inhibition of growth by polymers, such as the one reported by [START_REF] Lisitsin | Modeling the effect of anti-scalant on CaCO3 precipitation in continuous flow[END_REF], are unfortunately not directly compatible with such industrial applications as they have been developed for a limited number of inhibitors, including highly legislated ones (e.g., phosphate- based inhibitors). Such models are applicable to laboratory-scale setup using synthetic waters to mimic natural raw water; and these modeling approaches are not yet validated, to the best of our knowledge, on a larger scale with natural raw water. Notably, experimental measurements have been performed in recirculating cooling pilot plants [START_REF] Chien | Pilot-scale cooling tower to evaluate corrosion, scaling, and biofouling control strategies for cooling system makeup water[END_REF] using polymaleic acid (PMA), but the makeup water utilized originates from secondary-treated [START_REF] Li | Control of mineral scale deposition in cooling systems using secondary-treated municipal wastewater[END_REF] and tertiary-treated (Liu et al., 2012b) municipal wastewater and not river water, and no growth inhibition model has been proposed.

Hence, this study aims to fill the void with respect to the lack of models developed using natural water as the cooling fluid, validated for pilot-plant-scale systems. First, experiments are conducted in industrial pilot plants located in a nuclear power plant using raw river water for evaluating the efficiency of model polymers. Second, a model is proposed on the basis of the state-of-the-art models and experiments from previous studies. Finally, system-scale efficiency and model limitations, as well as perspectives for further development, are discussed.

Model development will assist researchers and engineers in designing water treatments and defining scale control strategies, focusing on the evaluation of ASP efficiency on a system scale for processes with high residence time, such as wet recirculating cooling circuits.

Experimental methods

Industrial pilot plant description

In this study, experimental setup is composed of four identical pilot plants, as shown in Figure 2, located in the Nogent-sur-Seine nuclear power plant (in France) owned and operated by the EDF ( Électricité de France) company. Each pilot plant is mainly composed of a low-pressure steam condenser (with heat duty up to 24 kW th ), a counter-current cooling tower, an ambient air fan equipped with humidity and temperature control, a washing tower, hot and cold basins, a recirculation pump, and injections pumps for various additives (e.g., acids, polymers, and biocides).

The pilot plant is designed to be representative of hydrodynamics and the water chemistry as makeup water is drawn off from the power plant. The residence time in each compartment is equal to the power plant values, and the pilot plant can run utilizing the same operating conditions as those for the power plant cooling circuit. The setup is also representative of materials in the main active components as the same metals are used in the condenser tube (e.g., stainless steel, brass, or titanium) and for the cooling tower packing (e.g., structured fills); ducts, and basins are nonetheless made of plastic instead of concrete. These pilot-plant units have been operated since 2001 in several EDF facilities for studies on fouling and biological growth, as well as for the development of associated treatments.

Pilot plant monitoring

Most of the operating parameters like the cooling water flowrate, the outlet condenser temperature, the cooling tower temperature approach, the makeup water flowrate, the concentration and flowrates of each additives are fixed by the operator. Optionally, the inlet air humidity and temperature can be fixed as well if necessary. The makeup water originating from the river, the water quality is not controlled and vary over time. In order to achieve the operating parameters targets, the air flow rate is controlled as well as the condenser heat duty and the blowdown water flowrate.

Since the purpose of this study is to evaluate the antiscale effect of polymers, the emphasis is put on how much calcium carbonate precipitates over time. In order to evaluate the precipitation rate, one could for example measure the deposit mass on specimens but specimen only represent scaling where they are installed. A more systemic indicator is the "water calcium loss" (∆Ca), which is basically a simple mass balance in steady-state conditions between the calcium in the makeup water and that in the blowdown water. It can be expressed as:

∆Ca [mg/L as CaCO3] = Q m [Ca 2+ ] m -Q b [Ca 2+ ] b Q b = CoC[Ca 2+ ] m -[Ca 2+ ] b (1)
Here, Q m and Q b are the makeup and blowdown water flow rates (in m 3 /s), CoC is the cycles of concentration (ratio of makeup and blowdown water flowrates) and [Ca 2+ ] is the calcium ion concentration expressed in terms of calcium hardness (in mg/L as CaCO 3 ). The numerator represents the amount of calcium precipitated in the whole circuit (in g/s of CaCO 3 ), it is divided by the blowdown water flow rate so as to have a common indicator compatible with both pilot-plant and power-plant scales.

Test campaign

Moreover, the overall campaign was conducted by maintaining one pilot plant as a reference (module No.4), where no polymer was added, and the other three modules were treated using polymers (modules No.1 to 3, respectively). The test campaign lasted from October 2014 to April 2015, implying a variable water calcium loss (∆Ca) during the whole campaign attributed to changes in river water quality. For compensating these variations and comparing the effects of polymers regardless of water quality, the overall efficiency of ASP η ASP is defined as follows:

η ASP = 1 - ∆Ca ASP ∆Ca ref mod (2)
with ∆Ca ASP representing the water calcium loss in the module operating with polymer, and ∆Ca ref mod representing the reference module operating under the same conditions but without polymers. Two sequential series of tests were performed using the four pilot plants. For comparing the efficiencies of the ASPs, three commercial polymers were used, one per module, under the same operating conditions (temperature, flow rate, and concentrations) for 9 weeks. The main available characteristics of the three ASP polymers are as follows: an acrylic acid copolymer of 4000 g/mol (ASP1), an acrylic acid homopolymer of 2000 g/mol (ASP2), and a carboxylate copolymer of 2000 g/mol (ASP3). During the second phase (10 weeks), ASP1 was used in all three modules for model calibration. During the whole campaign, cycles of concentration were maintained constant at 1.6, and the difference of temperatures in the cooling tower was 10°C. The investigated influence parameters are the polymer concentration in ppm of active matter, hot temperature (at the condenser outlet), and the amount of acid (for potential acid/polymer interaction). The main experimental conditions are reported in Table 1 in the Appendix, with make-up water parameters and operating conditions of the pilot plants, calcite saturation ratio Ω (i.e., the ratio of ion activity product and calcite equilibrium constant K sp ) reported in Table 1 have been calculated using PHREEQC software.

As the objective of the campaign is not to obtain zero scaling but rather to observe the inhibition of growth, scaling has to be observed during the campaign. This leads to polymer efficiencies of less than 100%. However, as the campaign lasted for several weeks, we did not want to induce massive scaling in the reference module; hence, sulfuric acid (0.18 mol/L) is always injected in all pilot plant units together with polymers so as to limit the observed scaling. It was required to adjust the acid flow rate to the variations of river chemistry; its addition is aimed at controlling the water chemistry in the cooling circuit and testing the potential dependency between the polymer efficiency and water alkalinity.

Modeling growth inhibition

Precipitation kinetic modeling in absence of inhibitor

For modeling the inhibition of growth by polymers, a model for the precipitation of CaCO 3 without polymers is required. The in-house "CooliSS" software [START_REF] Alos-Ramos | CooliSS (Cooling circuit Simulation Software): A computer program to simulate the water chemistry in recirculating cooling water circuits[END_REF] is employed herein, coded as a Scilab toolbox [START_REF] Rapenne | A Scilab toolbox for scaling prediction in condenser cooling circuits with cooling towers and its coupling with OpenTURNS[END_REF], with the possible use of a graphical interface. An overview of the CooliSS software is provided online as Supplementary Material, with the contributions of [START_REF] Alos-Ramos | CooliSS (Cooling circuit Simulation Software): A computer program to simulate the water chemistry in recirculating cooling water circuits[END_REF]; [START_REF] Marconnet | CooliSS: a modeling tool for scaling prediction in condenser cooling circuits with cooling towers[END_REF] and [START_REF] Mabrouk | Modelling scale formation in wet cooling system and validation with pilot plant data[END_REF]. In brief, the software utilizes a canonical chemical engineering approach for representing physicochemical phenomena in circuits. Chemical equilibria are modeled using the PHREEQC software [START_REF] Parkhurst | User's guide to PHREEQC (version 2) : a computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations[END_REF]. The model uses a rate-based formulation, considering the kinetic limitation for the precipitation or dissolution of CaCO 3 , as well as gas-liquid mass-transfer in the cooling tower for stripping as well as for the evaporation of water. The kinetics of calcium carbonate precipitation is taken from [START_REF] Appelo | A hydrogeochemical transport model for an oxidation experiment with pyrite/calcite/exchangers/organic matter containing sand[END_REF], which is based on an earlier model of [START_REF] Plummer | The kinetics of calcite dissolution in CO2-water systems at 5 to 60 C and 0.0 to 1.0 atm CO2[END_REF] This kinetic law has been derived from laboratory measurements using ideal solutions and has been adapted by the introduction of a critical saturation index [START_REF] Mabrouk | Modelling scale formation in wet cooling system and validation with pilot plant data[END_REF] for accurately predicting the saturation threshold observed in recirculating cooling circuits; the modified equation is expressed as follows:

R CaCO3 = r f 1 -10 2 3 (SI-CSI)
(3)

r f = k 1 a H + + k 2 a CO2 + k 3 a H2O (4)
Here, R CaCO3 represents the surface kinetic rate (mol.m -2 .s -1 , positive for dissolution and negative for precipitation), k 1 , k 2 , k 3 represent the temperature-dependent kinetic constants (m 3 .m -2 .s -1 , see [START_REF] Plummer | The kinetics of calcite dissolution in CO2-water systems at 5 to 60 C and 0.0 to 1.0 atm CO2[END_REF], a represents the activities of the species (mol.m -3 ), SI = log a Ca 2+ a CO 2- 3 /K sp represents the calcite saturation index and CSI the critical saturation index (substituted as 0.5 based on previous tests). The activities and calcite thermodynamic precipitation constant K sp are calculated using the PHREEQC software.

The mass-transfer between water and air in the cooling tower is modeled using the two-film theory (similarly as in [START_REF] Safari | Effect of CO2 stripping on pH in open-recirculating cooling water systems[END_REF] by considering the mass-transfer limitation and thermodynamic driving force:

J = K OL a gl (C L -H/C G ) (5)
With J represents the mass-transfer flux from the liquid to the gas phase (mol.m -3 .s -1 ), K OL represents the overall mass-transfer coefficient (m.s -1 ), a gl the gas-liquid interfacial area (m 2 .m -3 ), C L and C G represent the concentration in the liquid and the gas phase, respectively, and H represents the Henry's constant (-).

The pH affects the precipitation rate through the proton activity in eq. ( 4), its value is calculated using the PHREEQC software considering the water chemical equilibria, especially the calco-carbonic equilibrium modified by acid injection, and the release of dissolved gases to the air in the cooling tower using the eq. ( 5).

The reaction and stripping rates, equations (3) and ( 5), respectively, are valid only locally in the circuit and need to be implemented in mass balances in each equipment. Hence, a plug-flow regime is assumed in both the condenser and cooling tower, as well as perfect stirring behavior in basins and injections points. As a result, the solutions of transport, transfer, and reaction equations represent the concentrations and fluxes in the whole system.

Solving the transport, reaction and transfer equations requires the knowledge of water temperature path along the circuit as well as the water evaporation rate in the cooling tower. Experimental values can be used if these information are available to the user, thermal calculations are then performed by assuming linear temperature profiles in condenser and cooling tower and performing heat balances in each unit operation. Alternatively, a dedicated tool called TEFERI [START_REF] Bourillot | TEFERI: numerical model for calculating the performance of an evaporative cooling tower[END_REF]) can be used for thermal calculations. The program TEFERI calculates the cooling performances of cooling towers as a function of inlet water (temperature, flow rate), ambient conditions (e.g. air temperature and humidity) and the characteristics of the evaporative cooling towers provided by the packing manufacturer or by in-house measurements (e.g. Height of a Transfer Unit, HTU vs. the ratio of water-air mass flow rates). The software is based on the Poppe and Merkel methods for the calculation of evaporative cooling. For more details on cooling tower calculations, see for example the chapter dedicated to cooling towers written by [START_REF] Erens | VDI Heat Atlas[END_REF] in the VDI Heat Atlas.

As an example of typical temperature path, the outlet condenser temperature is around 40 °C (generally between 35 and 45 °C), the outlet cooling tower temperature is around 30 °C (the cooling water temperature difference could be between 5 and 15 °C depending on cooling tower design). The water evaporative rate is around 12 kg/h for a pilot plant with a make-up flow rate of 34 kg/h, and cycles of concentration (CoC) of 1.6. These values are order of magnitude and are highly dependent on operating and ambient conditions. More details on the counter-current cooling tower design procedure and order of magnitudes of various quantities can be found in [START_REF] Chien | Pilot-scale cooling tower to evaluate corrosion, scaling, and biofouling control strategies for cooling system makeup water[END_REF].

On the other hand, the software has been calibrated with pilot-plant measurements and validated using industrial data from EDF nuclear power plants. Currently, this software is employed by engineers for evaluating the impact of treatments on the cooling circuit scaling and water quality, for example, to define dynamic acid injection control strategies [START_REF] Marconnet | CooliSS: a modeling tool for scaling prediction in condenser cooling circuits with cooling towers[END_REF]. Before using the CooliSS software for developing a new polymer inhibition model, it is essential to verify that these simulations accurately predict the precipitation of CaCO 3 in the absence of polymers. Hence, simulations are performed using the makeup water characteristics (concentration, temperature, and flow rate) and pilot-plant operating conditions (heat duty, temperature, cycle of concentration, acid flow rate). Experimental water calcium losses are shown in Figure 1, with corresponding simulations for the reference module without the polymer. Along the test campaigns, the agreement between simulations and experiments is satisfactory as variations (attributed to the evolution of makeup water and/or modifications of operating conditions) are accurately predicted with an average deviation of 13 mg/L as CaCO 3 (around 9% relative deviation). Given the experimental uncertainties of the water calcium loss estimated at around 15 mg/L as CaCO 3 , the model is thus considered to be sufficiently robust for supporting the complementary inhibition model.

Inhibition model development

As explained above, the adsorption of polymers on growth sites is the main mode of action for the inhibition of scale. Hence, the kinetic rate in the presence of inhibitor R CaCO3, ASP can be expressed as the kinetic rate in the absence of inhibitor R CaCO3 multiplied by the fraction of available growth site (1 -θ), where θ represents the fractional coverage by the adsorbed polymers:

R CaCO3, ASP = (1 -θ)R CaCO3 (6)
The kinetic rate R CaCO3 has already been implemented in CooliSS software, and θ can be described by a dedicated model detailed below. Within this framework, the effects of water quality and cooling circuit operation (e.g., concentration factor and injected acid) are considered by the precipitation rate R CaCO3 , mitigated by the calculation of the fraction coverage θ, dependent of temperature and inhibitor concentration. The fractional coverage phenomenon can be either the consequence of the adsorption of polymers on growth sites or the accumulation of inhibitors at the diffusion boundary in the narrow crystal growth sites [START_REF] Gill | Mechanism of Scale Inhibition by Phosphonates[END_REF]. In both cases, polymers prevent crystal precursors from reaching the growth sites, attributed to steric hindrance and repulsive force caused by static electricity.

Laboratory-scale measurements by [START_REF] Lisitsin | Modeling the effect of anti-scalant on CaCO3 precipitation in continuous flow[END_REF] as well as results obtained from the pilot-plant test campaign (see Results section),suggest a type I adsorption isotherm, typically well represented by the Langmuir adsorption model. In this paper, the mathematical model of [START_REF] Kubota | A kinetic model for crystal growth from aqueous solution in the presence of impurity[END_REF] is utilized, which modifies the original Langmuir adsorption model by introducing an impurity effectiveness factor α to account for the effectiveness of a given impurity for growth suppression (for the physical meaning of α, please see the original paper of [START_REF] Kubota | A kinetic model for crystal growth from aqueous solution in the presence of impurity[END_REF]. Hence, the mathematical model is expressed as follows:

θ(T, [ASP]) = α K a (T )[ASP] 1 + K a (T )[ASP] (7)
© 2016. This manuscript version is made available under the CC-BY-NC-ND 4.0 license

K a (T ) = K 0 a e -∆H 0 a /(RT ) (8)
Here, θ represents the fractional coverage of Eq. 6 (dimensionless), α represents the effectiveness factor (dimensionless), [ASP] represents the ASP concentration (in mg/L of active matter), K a represents the Langmuir equilibrium constant (here in L/mg) whose dependency on temperature is given by Eq. 8, K 0 a is the pre-exponential factor (in L/mg), ∆H 0 a is the standard enthalpy change (J/mol), R represents the ideal gas constant (J/mol/K) and T represents the absolute temperature (K). The three parameters α, K 0 a and ∆H 0 a are identified by experimental measurements obtained during the test campaign.

Results and discussions

Experimental results

In the following sections, the reported polymer concentration is the overall concentration injected in the circuit (in ppm of active matter), which should not be confused with its residual concentration in water. Figure 4 shows the water calcium loss ∆Ca during the test campaign (see the Appendix for details on the operating conditions and river water quality during the whole campaign).

Comparison of polymer efficiency. The first phase aims to compare the ASP efficiency of three commercial polymers: ASP1, ASP2, and ASP3; the efficiency can be observed by comparing the water calcium loss in the circuits treated with polymers (pilots 1-3) and that in pilot 4 without the polymer. As shown in Figure 4 (top), ASP1, ASP2, and ASP3 exhibit comparable efficiencies and similar responses to changes in operating conditions such as polymer concentration, temperature, and water chemistry. Despite experimental uncertainties (estimated around 15 mg/L as CaCO 3 ), ranking is observed, with ASP1 being the most efficient; this ranking can be explained by the differences in the polymer characteristics such as chemical functions and molar mass (as illustrated by [START_REF] Jada | Effect of sodium polyacrylate molecular weight on the crystallogenesis of calcium carbonate[END_REF][START_REF] Al-Hamzah | Inhibition of homogenous formation of calcium carbonate by poly (acrylic acid). the effect of molar mass and end-group functionality[END_REF]; however, it is difficult to proceed with the analysis with the marginal information given by polymer suppliers.

Nucleation inhibition and growth inhibition. Irrespective of the polymer chosen, a sharp increase for the water calcium loss is observed between weeks 3 and 4, shifting from a null calcium balance (i.e., 100% polymer efficiency) to positive values; this shift matches the increase of the calcite saturation ratio of the makeup water during the week end between those weeks (see Appendix), without the adaptation of the acid flow rate. Furthermore, notably, there is no return to the original state as similar operating conditions before and after this event lead to different water calcium losses. Conditions are similar during week 2 (outlet condenser temperature of 40 °C, 3 ppm of ASP) and weeks 8 and 9 (outlet condenser temperature of 40 °C, ASP concentrations of 2.5 and 4 ppm, respectively), and the water calcium losses during week 2 are almost equal to zero in the three pilot plants with ASP, and between 50 and 120 mg/L during weeks 8 and 9. This observation is interpreted as a change in the mode of action of the polymer. As indicated in the introduction, the polymer can prevent the appearance of crystalline germ, i.e., inhibits nucleation or limits crystalline growth, that is, inhibits growth kinetics. Initially, with the cleaned pilot plants, the addition of polymers was probably sufficient for preventing nucleation for three weeks. However, with a higher saturation index, the nucleation energy barrier decreases; hence, it is necessary to increase the polymer concentration, and it is difficult to reach the concentration threshold. Consequently, it is assumed that the polymer acts as a nucleation inhibitor during the three first week, while the circuit was clean, then acts as a crystal growth inhibitor once nucleation occurred, attributed to a high residence time of the water in the recirculating circuit and a high calcite saturation ratio.

Effect of operating parameters on polymer efficiency. The individual effect of operating parameters is observed by gathering data obtained under identical conditions. The overall efficiency of ASP polymer η ASP is plotted in Figure 5 with respect to the polymer concentration (Figure 5, top) and condenser outlet temperature (Figure 5, bottom), with other parameters maintained constant while screening the investigated parameter. As shown in Figure 5 (top), the ASP efficiency rapidly increases with polymer concentration and reaches a constant value after 4 ppm. However, the efficiency decreases with increasing temperature, as shown in Figure 5 (bottom). These two observations, already observed previously [START_REF] Al-Hamzah | Inhibition of homogenous formation of calcium carbonate by poly (acrylic acid). the effect of molar mass and end-group functionality[END_REF][START_REF] Sparks | Adsorption of poly acrylic acid onto the surface of calcite: an experimental and simulation study[END_REF], reinforce the choice of an adsorption model as adsorption is favored at low temperature, and Figure 5 ΔCa, water calcium loss

Modeling results

As ASP1 exhibits higher efficiency, it was used in the second series of experiments in the three pilot plants. The goal was to acquire more data for the parametric identification of the model using the ASP1 polymer. Three parameters, the effectiveness factor α, pre-exponential factor of the Langmuir equilibrium constant K 0 a and the standard enthalpy change ∆H 0 a (see equations 7 and 8) are to be identified. The effectiveness factor α is assumed to be constant; however, it could represent a function of the physicochemical properties in a further version of the model. In total, 37 points are selected among the stabilized data for covering the concentration range from 0.25 to 9 ppm and the temperature range in the circuit from 20 to 45°C, and the three parameters regress with a least-square method on water calcium loss. The optimal set of parameters for ASP1 is as follows:

   α = 0.97 ln K 0 a = -33.26 (K 0 a in L.mg -1 ASP ) ∆H 0 a = -91.8 kJ.mol -1 (9)
The average absolute deviation between regressed and experimental data is 10 mg/L as CaCO 3 for the 37 chosen points; hence, a satisfactory fit is obtained as the deviation is below the experimental uncertainties. For validating the model, 43 additional points from the campaign are considered, as well as 25 points from a former campaign in 2012. The points from the former 2012 tests are obtained under more constraining conditions; the cycles of concentration were high, with less or no acid at all, and more alkaline river water compared to the 2014-2015 tests. Figure 6 shows the parity plot representing calcium losses for the simulated and pilot plant water for all 105 points. The points used for regression are well distributed along the first bisector in the regression range between 30 and 180 mg/L, as well as the validation points in the same range. In extrapolation, the points are also well distributed for a high water calcium loss, corresponding to the 2012 tests under conditions prone to the precipitation of calcite. However, the model tends to overpredict water calcium losses in the low-scaling region as it simulates scaling around 50 mg/L as CaCO 3 , whereas no scaling is practically observed on pilot plants. For all 105 points, the average absolute deviation between measurements and simulations is 21 mg/L as CaCO 3 .

Figure 7 shows the comparison between the water calcium loss obtained with ASP1 on pilot plant 1 (same values as on Figure 4) and the corresponding simulations with CooliSS software, modified with equations ( 6) to (8) and parameters of equation ( 9). As the circuit was clean during the three first weeks, simulations are not performed for the period as the model is not adapted; hence, a nucleation model could be used instead. The model accurately reacts to changes in river water quality as well as modifications of operating parameters, except for two weeks (13 th and 14 th ) where the software predicts a positive water calcium loss, but no scaling is observed on pilot plants. These weeks correspond to the points of higher deviations in Figure 6 in the lowscaling region at high acid flow rate, suggesting a missing interaction in the model construction (see discussions on model limitations). Weeks -pilot plant n°1

Simulations with growth inhibition model Pilot plant

Clean circuit

Scaled circuit

[mg/L as CaCO 3 ]
ΔCa, water calcium loss 4.3. Comments on recirculating cooling circuits Effect of polymer treatment. As mentioned earlier, two of the main characteristics of recirculating cooling circuits are the salt concentration by water evaporation (represented by the cycles of concentration, CoC) and the important water residence time, attributed to the large water volume in the recirculation loop. The higher the CoC, the higher the calcite saturation ratio in the circuit. If Ω river represents the calcite saturation ratio in the river, the saturation ratio in the circuit would then be CoC 2 × Ω river in the absence of scaling. This increase of the saturation ratio lowers both nucleation and crystal-growth thermodynamic energy barriers. In addition, the mean residence time is typically around a few hours, which is quite superior to the calcite induction period.

In the pilot plants, the mean residence time was set to 5 h, implying that it takes at least a day to evacuate 99% of inlet salts, considering a first-order response, i.e., assuming that the circuit is a stirred tank without dead volumes (which is obviously not the case). In comparison, the induction period is estimated between 3 and 15 min during the test campaign using the formula from [START_REF] He | Inhibition of calcium carbonate precipitation in NaCl brines from 25 to 90 •[END_REF].

Owing to both thermodynamic and kinetic aspects, maintaining a clean recirculating cooling circuit would be possible for low CoC and/or low saturated waters. The minimal polymer dose could be estimated by the models reported by [START_REF] Vanderpool | Calculating minimum threshold inhibitor dosage[END_REF] or [START_REF] He | Inhibition of calcium carbonate precipitation in NaCl brines from 25 to 90 •[END_REF]. In other cases, the polymer acts as a growth inhibitor with efficiencies described herein. According to equation ( 7), the fractional growth sites coverage θ approaches α (very close to 1) at high polymer concentrations; θ represents the relative amount of inaccessible growth sites, caused by the presence of the polymer, either by the direct adsorption of the polymer or by electrostatic repulsion. However, the relationship between θ and the inhibition of calcite precipitation in the circuit η ASP is not linear, and a fractional coverage of θ = 0.5 would not lead to a 50% inhibition efficiency. For illustrating this point, Figure 8 shows the antiscalant efficiency, systemically measured for pilot plants, versus the mean fractional coverage calculated by the model for all tests. The water calcium loss exhibits a moderate decrease with fractional coverage (i.e., overall efficiency increase), followed by a stiff increase above a fractional coverage of 0.8. This nonlinear curve is caused by the plug-flow regime in the condenser and cooling tower, and the reaction rate is integrated along the flow path; hence, the concentration exiting such equipment is not proportional to the reaction rate, but it follows a logarithmic evolution. This result explains the nature of the dependency between fraction coverage, which directly modifies the kinetic rate, and the observed efficiency in the circuit. The complete curve shown in Figure 8 is obtained with CooliSS by manually modifying the value of θ from 0 (no coverage) to 1 (full coverage).

Complementary effect with acid injection. The injection of acid, as a stoichiometric treatment of water alkalinity, would lead to excessive use of chemicals for cooling circuits with high-scaling issues. Hence, the addition of a few parts per million of polymers would result in an adverse decrease in the precipitation of calcium carbonate. For instance, Figure 9 shows the (normalized) amount of acid to inject in the circuit for maintaining scaling below an arbitrary water calcium loss of 30 mg/L as CaCO 3 , for three river water hardness values. For soft water (dotted curve), a small amount of acid is required to limit scaling without the polymer (a fifth of the amount for hard water), and the addition of polymer decreases this amount to zero after 2 ppm. Moreover, for moderately hard water (dashed curve) and hard water (plain curve), the addition of polymer also decreases the acid flow rate, but acid is still required even for high polymer concentration, attributed to the adsorption plateau. For hard water, a combination of both treatments would probably be the most effective and economic for limiting scaling without the consumption of large amounts of chemicals and polymer injection complemented by acid injection as 100% efficiency with only use of polymers is difficult to achieve because of the adsorption plateau.

Model limitations and further development

The developed model allows for the prediction of the efficiency of the chosen ASP with respect to the cooling circuit design, make-up water quality, and operating conditions. Although the model is validated for an extended range of water qualities and operating conditions, experimental results are representative for pilot plants with specific surface materials and suspended solids in the Seine river natural water. Further development is required for enhancing model predictability, including the effect of surfaces (e.g., concrete, heat exchanger metallic surfaces, plastic packing, and suspended matter) and additional water properties on model parameters for ensuring better transposability from one circuit to another. For example, discrepancies between simulations and experiments (see Figure 7) suggest the dependency of the pH on the polymer efficiency. For instance, [START_REF] Wisniewska | Comparison of adsorption affinity of polyacrylic acid for surfaces of mixed silica-alumina[END_REF] have reported the effects of solution pH on the adsorption isotherms of polyacrylic acid (PAA) on mixed silica alumina, caused by the modification of surface charge density; [START_REF] Sparks | Adsorption of poly acrylic acid onto the surface of calcite: an experimental and simulation study[END_REF] © 2016. This manuscript version is made available under the CC-BY-NC-ND 4.0 license have reported the effect of pH on the energies of adsorption of PAA on calcite. Hence, the adsorption model could include such dependencies, in addition to temperature and polymer concentration.

In addition, it would be interesting to investigate the effect of the polymer characteristics (e.g., repeat units, structure, and molar mass distribution) on model parameters. This would let us limit the need of a long test campaign on pilot plants for new polymers. For instance, [START_REF] Kubota | A kinetic model for crystal growth from aqueous solution in the presence of impurity[END_REF] have investigated the crystal growth of KBr in the presence of aliphatic carboxylic acids and found that the effectiveness factor α increases with the carbon number of the acid, while the adsorption equilibrium constant K remains nearly constant. Such quantification between polymer properties and antiscalant efficiency would allow for the estimation of the model parameters with laboratory-scale measurements; nevertheless, the chemical nature of the antiscalants would need to be known.

Conclusions

This study provides comprehensive knowledge on the antiscalant efficiency of polymers in a recirculating cooling circuit by using industrial pilot plants and natural river water as the coolant. Three commercial polymers were tested, and they exhibited comparable antiscalant efficiencies and similar responses to changes in operating conditions. This study demonstrates the characteristics of a recirculating cooling circuit, where the inhibition of calcite occurs because of crystal-growth inhibition rather than nucleation inhibition, attributed to the high saturation ratio and long residence time in the system.

The results obtained from the industrial pilot plant are in agreement with laboratory-scale observations, particularly the identification of crystal-growth inhibition by an adsorption mechanism. In addition, interesting insights on scale mitigation are obtained for pilot plants, such as the non-linear relation between the fraction growth-site coverage and the scale inhibition in the circuit. The use of an industrial pilot plant is critical for accurately simulating a full-scale cooling circuit and scale inhibition.

The developed adsorption model related the raw river water quality and process operating conditions (e.g., cycles of concentration, temperatures, and amount of treatments) to the ASP efficiency. This model could be used to design treatments and define control strategies. The model could be improved to consider the dependency of model parameters to additional operating parameters such as solution pH and to the polymer nature.

Table 1: Experimental conditions (weekly averaged) during the test campaign; three different ASP were used during the first phase (efficiency comparison, weeks 1 to 9), and only ASP1 was used during the second phase (parametric screening, weeks 10 to 19) Week

Figure 1 :

 1 Figure 1: Wet recirculating cooling water system in a power plant

Figure 2 :

 2 Figure 2: Schematic of a pilot plant (left) and photograph of the cooling tower skids in Nogent-sur-Seine (right)

Figure 3 :

 3 Figure 3: Comparison between experimental water calcium loss for the reference module (without the polymer) and CooliSS simulations (no fitting performed)

  (top) is similar to a type I adsorption isotherm. Manuscript accepted in Desalination. DOI: 10.1016/j.desal.2016.06.018 © 2016. This manuscript version is made available under the CC-BY-NC-ND 4

Figure 5 :Figure 6 :

 56 Figure 5: Effect of polymer concentration and, condenser outlet temperature on the overall efficiency of an antiscalant polymer

Figure 7 :

 7 Figure 7: Experimental water calcium balance and CooliSS (modified with ASP inhibition model) simulations for all the data of pilot 1

Figure 8 :

 8 Figure 8: Link between the value of θ and the precipitation inhibition

Figure 9 :

 9 Figure9: Examples of the trade-off between polymer and acid injection to maintain scaling below a water calcium loss of 30 mg/L as CaCO 3 , simulated for several river water hardness

© 2016. This manuscript version is made available under the CC-BY-NC-ND 4.0 license

Supplementary data

Supplementary material related to this article (overview of the original CooliSS software, i.e. without polymer addition) can be found online at 10.1016/j.desal.2016.06.018.

Appendix: experimental conditions during tests


Xavier Bry

Théo Simac

Salah Eddine El Ghachi

Philippe Antoine

Bridging data-exploration and modeling in event-history analysis: the Supervised-Component Cox Regression method

Keywords: Supervised components, Cox regression, PLS Cox regression, regularization, SCGLR, survival analysis

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Statistical models based on stochastic counting processes are one of the most useful mathematical tools for population study [START_REF] Andersen | Statistical Models Based on Counting Processes[END_REF].

When they include the effects of covariates in a generalized linear formulation of their intensity, they provide a strong and handy explanatory way to link response-variables to what supposedly influences them, and quantify this statistical link. Among such explanatory models, Cox's proportional hazard regression model is most popular for having long been successfully used in epidemiology for survival analysis [START_REF] Cox | Regression Models and Life-Tables[END_REF][START_REF] Cox | Analysis of Survival Data[END_REF] ; [START_REF] Collett | Modelling Survival Data in Medical Research[END_REF] and in demography for event-history analysis [START_REF] Courgeau | Analyse démographique des biographies[END_REF]. These generalized linear counting-process models, as powerful as they are, have the same Achilles' heel as all classical regression models: they cannot handle high-dimensional data, i.e. data where explanatory variables outnumber observations. Nor can they handle highly correlated explanatory variables. High dimension and high correlation make such models unidentified if one does not add some constraints or penalty to the goodness-of-fit criterion they currently maximize. Theory-based crisp identification constraints, when they do exist, are usually too rare to make the model identified. Moreover, in practice, they do not exist at all, for the excess of variables over observations and the presence of highly redundant variables arise from an "excess of measurement" : as the "true" explanatory variables can usually not be observed directly, they are replaced with a number of proxies supposedly revolving "about" some of them, which is why they are so correlated. This lack of constraints makes the penalty-based regularized models, such as ridge-or least absolute shrinkage and selec-tion operator-penalised regression models [START_REF] Tibshirani | The Lasso method for variable selection in the Cox model[END_REF][START_REF] Perperoglou | Reduced-rank hazard regression for modelling non-proportional hazards[END_REF], an all-terrain and very much appreciated way of breaking down the curse of dimensionality. But these regularized models still miss one point: if they do provide a linear predictor of the response, this linear predictor is very hard to interpret, for the variables in it are too many and redundant. What the analyst needs is a small set of explanatory, i.e. predictive and interpretable dimensions which capture the statistical link between the response and the set of its regressors. So, true explanatory modeling in a high-dimensional framework requires, on top of regularization, analytical dimension reduction.

On the other hand, classical exploratory analysis methods as principal component analysis (PCA) and correspondence analysis are justly reputed as powerful dimension-reduction techniques, but have no explanatory power of their own. Now, in-between regression and principal component analysis lies partial least square (PLS) regression. The basic idea of partial least squares regression [START_REF] Wold | The collinearity problem in linear regression. the partial least squares (PLS) approach to gener-59 alized inverses[END_REF] is to bridge principal component analysis and ordinary least squares (OLS) regression by maximizing a trade-off criterion between ordinary least squares regression's goodness-of-fit criterion (square correlation of dependent variable with component) and the compo-nent's "strength", as measured by its variance under a unit-norm constraint on the coefficient-vector. But partial least squares regression is not at all straightforward to extend to generalized linear models. [START_REF] Bry | / Exploring explanatory models ; an event history application[END_REF], [START_REF] Bastien | PLS generalised linear regression[END_REF] and [START_REF] Bry | Extension de l'Analyse en Composantes Thématiques univariée au modèle linéaire généralisé[END_REF] have proposed different ways to do it.

In our view, all had to be improved in terms of generality, flexibility, and mathematical consistency. More recently, [START_REF] Bry | THEME: THEmatic Model Exploration through Multiple Co-Structure maximization[END_REF] have proposed to extend the measure of a component's strength beyond its variance, by introducing a flexible measure of structural relevance. This measure can act as a bonus given to components close to certain variable-structures, e.g. bundles of highly correlated variables. This is an important asset in "explanatory exploration" , for it produces components all the easier to interpret as they are close to enough explanatory variables. Below, we develop in detail what was briefly introduced by [START_REF] Bry | Supervised-Component based Cox Regression[END_REF]. We propose to flexibly combine the structural relevance with the partial likelihood of the proportional hazard model into a new criterion. Maximizing this criterion under a unit-norm constraint on the coefficient vector allows to calculate a first structurally relevant explanatory component. A rank-h component can then be calculated by maximizing the same criterion under the additional constraint that it be orthogonal to all lower-rank components. As the calculated components are supervised ones, they are not exogenous, which entails that standard significance-tests are invalid. As a consequence, once a given (and presumably too large) number of components is obtained, a cross-validation procedure is required to only keep the really predictive ones. Eventually, one gets a regularized linear predictor of the hazard function, decomposed on a small number of dimensions spanning an interpretable explanatory subspace. By varying the parameter tuning the combination of the structural relevance with the partial likelihood, one can make the method go smoothly from standard non-regularized Cox regression to Cox Regression on such exogenous (non-supervised) "strong" components as principal components.

These two extremities are usually not very interesting: on the one hand, in a high-dimensional setting, standard Cox regression yields highly unstable and hardly interpretable coefficients, if any, for want of prior dimension reduction; on the other hand, exogenous components are not optimized to predict the survival time, so do not usually provide the best explanatory view of it. In practice, the most interesting components will be given by some intermediate values of the tuning parameter.

The paper is organized as follows. In section 2, after a short reminder on Cox's proportional hazard model, we present the corresponding componentmodel and its estimation technique in a formal way. In section 3, a simulation study is performed with intent to show how the tuning-parameters condition the estimation results, and how to interpret them. Finally, in section 4, the hints given in section 3 are used to analyse life-history data of men in Dakar, in order to find which variables are useful to model their change from monogamy to polygamy.

Model and technique

We consider a survival time Y , depending on a set of covariates X, plus a set of extra-covariates Z. X and Z can be time-dependent. There may be noninformative right-censoring on Y . Variables in Z are few and exhibit low or no correlation. By contrast, variables in X are many and possibly redundant, so that the proportional hazard model demands regularization with respect to X.

Notations:

• y i is the survival-time or censoring-time of unit i ∈ {1, ..., n}.

• x i,t and z i,t are the values of vectors of covariates x and z respectively for unit i at time t.

• Event-indicator δ is defined through: ∀i, δ i = 1 if the event occurs for i at time y i , and δ i = 0 if i is censored at time y i .

• R(t) denotes the set of all individuals at risk at time t.

• Π A , where A is a matrix, denotes the orthogonal projector on the space spanned by the column-vectors of A, with respect to a given metric.

• A being a matrix, A denotes the transpose of A.

The proportional hazard model

Cox's model is based on the following formulation of the hazard function of unit i at time t:

h(t; x i,t , z i,t ) = h 0 (t)e β x i,t +γ z i,t , (1)
where h 0 (t) is the baseline hazard function. The survival function having hazard function h 0 (t) is the baseline survival function.

The partial likelihood defined by [START_REF] Cox | Regression Models and Life-Tables[END_REF][START_REF] Cox | Partial Likelihood[END_REF] is, when there are no simultaneous events:

L p (β , γ; X, Z) = n ∏ i=1 p i , where p i = e β x i,y i +γ z i,y i ∑ j∈R(y i ) e β x j,y i +γ z i,y i δ i (2)
The partial likelihood can, under some assumptions, be interpreted as a marginal likelihood of events' ranks [START_REF] Cox | Partial Likelihood[END_REF]. The partial likelihood is rid of h 0 , involving only (β , γ). When the model is identified (i.e. matrix [X, Z] is full column rank), maximizing the partial likelihood with respect to (β , γ) through a Newton-Raphson algorithm yields estimates ( β , γ) based on which [START_REF] Kalbfleisch | Marginal likelihoods based on Cox's regression and life model[END_REF] and [START_REF] Breslow | A large-sample study of the life table and product limit estimates under random censorship[END_REF], among others, proposed an estimation of the baseline survival function.

The problem with our data is the large amount of collinearity within X, which causes the X-part of the linear predictor to be unstable, if identified.

We shall make it identified and stable by prompting it to lean on components that have some structural relevance. The idea is to replace regressor-block X with a block F = XU of H orthogonal components, in the model. The regressors being time-dependent, so will the components. Let X be the matrix the columns of which are the X-regressors and the N rows of which are the individuals-at-risk-at-time-points: (i,t). Orthogonality of the f 's will be taken with respect to matrix W = N -1 I N . The hazard function of a unit i at time t will thus be:

h(t; x i,t , z i,t ) = h 0 (t)e α f i,t +γ z i,t = h 0 (t)e α f i,t +γ z i,t (3)
Components will be calculated hierarchically, starting with one, and each extra-component being constrained to be orthogonal to the former ones.

Now, what we want is that the components be structurally relevant, in that they should be as close as possible to directions of a pre-defined type, such as explanatory variables, or other relevant subspaces.

Here is how we propose to measure that.

Structural relevance of a component

To X we associate a p × p symmetric positive definite (s.p.d.) matrix A, such that principal component analysis of X with metric matrix A and weightmatrix W is relevant.

Let τ ∈ [0; 1], and M = (τA -1 + (1τ)X WX) -1 . The purpose of coefficient τ is to tune the regularization of the model as follows. Matrix M is s.p.d., and component f = Xu will be constrained by : u 2 M -1 = 1, so that:

• τ = 0 ⇒ M = (X WX) -1 , so that u 2 M -1 = 1 ⇔ f 2 W = 1 . All
directions in X are equivalent, be them close to subsets of variables or not: X is then seen as a pure vector-space.

• τ = 1 ⇒ M = A . Recall that the program of the principal component analysis of (X, A,W ) is:

max w Aw=1 XAw 2 W ⇔ max u A -1 u=1 Xu 2 W with u = Aw (4) So, s.t. u 2 A -1 = 1 , f 2 W is the inertia of observations (rows) along w = A -1 u ∈ R p .
Components with a larger inertia will hence be favoured over those with a smaller one.

Structural relevance was introduced by [START_REF] Bry | THEME: THEmatic Model Exploration through Multiple Co-Structure maximization[END_REF] as a possible extension of the component's variance to measure the ability of a component to capture information in its variable-block. Given a set of "reference" s.p.d. matrices N = {N j ; j = 1, ..., J} encoding types of structures of interest (target-spaces, e.g. variables in X), a weight system Ω = {ω j ; j = 1, ..., J}, and a scalar l ≥ 1, the associated structural relevance measure is defined as a generalized average of quadratic forms of u:

φ N,Ω,l (u) := J ∑ j=1 ω j (u N j u) l 1 l (5)
In Eq. ( 5), the value of l tunes the locality of the bundles of structures coded in N. The larger the value of l, the more local the bundle. An example is given below.

Particular instances of structural relevance measures

• Component Variance:

φ (u) = V (Xu) = Xu 2 W = u (X WX)u
Under constraint u M -1 u = 1, this is the inertia of units along direction u , and is maximized by the first (direct) eigenvector in the principal component analysis of (X, M,W ) .

In practice, explanatory variables are often a mixture of numeric and nominal variables. Assume that X = x 1 , ..., x K , X 1 , ..., X L , where: x 1 , ..., x K are column-vectors coding the numeric regressors, and X 1 , ..., X L are blocks of centred indicator variables, each block coding a nominal regressor ( X l has q l -1 columns if the corresponding variable has q l levels, the removed level being taken as "reference level"). We should then consider the following A , which bridges ordinary principal component analysis of numeric variables with multiple correspondence analysis of nominal variables:

A := diag (x 1 W x 1 ) -1 , ..., (x K W x K ) -1 , (X 1 W X 1 ) -1 , ..., (X L W X L ) -1 (6)
• Variable powered inertia :

For a block X consisting of p standardised numeric variables x j , the variable powered inertia is defined as:

φ (u) = p ∑ j=1 ω j Xu|x j 2l 1 l = p ∑ j=1 ω j (u X W x j x j W Xu) l 1 l (7) ⇔ φ (u) = p ∑ j=1 ω j ρ 2l ( f , x j ) 1 l f 2 W (8)
Assuming we take τ = 0, we have f 2 W = 1 and:

<u> ϕ 1 (u) ϕ 2 (u) ϕ 4 (u)
In the elementary case of 4 coplanar variables x with ∀ j, ω j = 1 , fig. 1 graphs φ l X (u) in polar coordinates (using the complex notation, where vector u is identified with complex number e iθ : z(θ ) = φ l X (e iθ )e iθ ;θ ∈ [0, 2π) ) for various values of l . Note that φ l X (u) was graphed instead of φ X (u) so that curves would not overlap. One can see how the value of l tunes the locality of bundles considered: the greater the l , the more local the bundle.

For a block X consisting of p categorical variables X j , each of which is coded through the set of its centred indicator variables less one, the variable powered inertia is:

φ (u) = p ∑ j=1 ω j cos 2l (Xu, X j ) 1 l = p ∑ j=1 ω j Xu|Π X j Xu l W 1 l (10)
where:

Π X j = X j (X j W X j ) -1 X j W

Combining the partial likelihood with the structural relevance

We propose to combine the partial likelihood with the structural relevance using a geometric average:

c(u, α, γ;Z;s) = [L p (u, α, γ)] 1-s [φ (u)] s with 0 ≤ s ≤ 1 (11)
The scales of L p and φ are not comparable, and the geometric average has the obvious advantage to make the compound criterion insensitive to these scales in that, at the optimum, the relative variations of L p and φ compensate with a fixed rate:

dc = 0 ⇔ d ln c = 0 ⇔ (1 -s)d ln L p + sd ln φ = 0 ⇔ dL p L p = - s 1 -s dφ φ (12)

The Supervized Component Cox Regression algorithm

Maximizing the criterion c in Eq. ( 11) subject to constraint u M -1 u = 1 yields vector u 1 , hence the first component :

f 1 = Xu 1 .
When looking for rank-h component f h = Xu h , we take former compo-

nents F h-1 = [ f 1 , ..., f h-1 ]
as known, and impose the additional orthogonality constraint:

F h-1 W f h = 0 ⇔ D h u h = 0 with D h = X W F h-1 (13)
Besides, F h-1 is taken as a block of extra-covariates, thus appended to

Z. Let Z h = [Z, F h-1 ].
Vector u h is obtained as the solution of the following program:

P : max u, α, γ s.t.: u M -1 u = 1; D h u = 0 c(u, α, γ; Z h ; s) (14)
Program P is solved by iteratively maximizing the criterion with respect to u and (α, γ), in turn:

• Maximization with respect to u : given (α, γ) , the criterion is maximized on u s.t. u M -1 u = 1; D h u = 0 through the projected iterated normed gradient algorithm given in appendix B. The formula for calculating gradient ∇c may be found in appendix A.

• Maximization with respect to (α, γ) : given u , the criterion is maximized on (α, γ) exactly as in the classical Cox regression on covariates

{Xu, Z h } , with Z h = Z ∪ F h-1 .
Overall, we get the following algorithm:

-A number H of components to be calculated is chosen. Let D 1 = nullmatrix, F 0 = / 0 and Z 1 = Z. Then:

For h = 1 to H :

-Calculate u h as the solution of P with Z h as extra-covariates.

-Set

f h = Xu h , F h = [F h-1 ; f h ] , D h+1 = X W F h and Z h+1 = [Z; F h ] ;
End for.

-Perform Cox regression on regressors [F H , Z] . This yields linear predictor:

η = F R α + Zγ = XU R α + Zγ = Xβ + Zγ, whereβ = U R α (15)
If X denotes the original uncentered-variable matrix, then:

X = X -1 N W X , so: (16) η = β 0 + X β + Zγ , with β = U R α and β 0 = -1 N W X β . (17)

Model-assessment

As said earlier, when s < 1, the components are supervised, hence nonexogenous, which requires that their assessment be made through crossvalidation.

Cross-validation quality coefficient

To each triplet (s, l, H) we associate a corresponding model M .

CV k (M ) = l(θ -k (M )) -l -k (θ -k (M )) (18)
where θ = (β 0 , β , γ) , l -k is the log-partial likelihood excluding part k of the data, and θ -k (M ) is the θ (M ) estimated on the non-left out data. The overall CV (M ) is the average of {CV k (M ); k = 1, ..., K} .

Obtaining good values for the tuning-parameters

We shall, from now on, use the variable powered inertia to measure structural relevance. The tuning parameters are many (s, l, τ, H), so that using cross-validation to compare all combinations of them on a cross-product grid is out of the question in practice. Therefore, heuristics are needed. It is important to note that, even if these parameters have different purposes, which can to some extent be served sequentially, they are not completely independent. We propose to momentarily separate the exploratory goal from the predictive one, these two goals being served in this very order. Indeed, contrary to prediction, exploration does not require regularization, which simplifies tuning, and should lead to interpretable explanatory dimensions used then for reliable prediction.

Exploratory phase: The primary key-parameter to be tuned in this stage is s (importance of structural relevance with respect to partial likelihood), which entails the approximate number H of components to be retained ultimately. Choosing s = 0 (pure goodness-of-fit maximization) should lead to a single component: the predictor of the classical Cox regression. To the other end, s = 1 leads to exogenous components that do not take the goodness of fit into account, which means that many such components may be necessary to correctly fit the model. The higher s is, the higher H should be. In practice, we propose the following heuristic way to identify the explanatory dimensions. First, take a low value of τ, strictly positive to ensure identification but producing only minimal regularization, e.g. τ = 10 -2 .

Then:

(a) Perform principal component analysis on X (e.g. with s = 0 and (e) Select the value of τ leading to the best CV * (τ) , and try to fine-tune τ so as to yet improve it.

l = 1),
(f) In the process, see that the interpretation of components is not altered, which may be the case with large values of τ. Should this happen, either decrease τ, or re-interpret components.

Simulation study

The simulation scheme is intended to show how the tuning of parameters s, l and τ influence the estimation results, which gives some insight on their roles and how to use them in practice.

Simulation scheme

We defined the time-span as [0,30] , divided in 30 unit-length elementary intervals. over [0,30], a baseline hazard function was simulated as:

h 0 (t) = a + b(t -t m ) 2 with t m = 12, a = .2, b = 10 -3 (19)
Time-dependent covariates describing n = 75 subjects were then simulated with bundle-structures, as follows:

Firstly, three independent latent variables taking time-independent values at subject-level were simulated :

ψ j i ∼ N(0; 1), j ∈ {1, 2, 3}, i ∈ {1, ..., 75}.
Secondly, three independent latent variables were simulated on the 75 × 30 (subject, elementary interval) pairs: φ j it ∼ N(0; 1), j ∈ {1, 2, 3}, i ∈ {1, ..., 75},t ∈ {1, ..., 30}. Thirdly, these variables were combined into three latent vari-ables directing the covariate-bundles:

∀(i,t, j) : ξ j it = ψ j i + φ j it ( 20
)
Revolving around latent variables ξ 's,"observed" covariates were simulated as follows:

∀ j ∈ {1, 2, 3} : bundle B j = {x jk , k = 1, ..., K j } with x jk it = ξ j it + ε k it , ε k it ∼ N(0; σ 2 ) , with σ = .3. We took K 1 = 4, K 2 = 6, K 3 = 10.
Finally, 20 independent pure noise covariates were simulated as N(0;1).

Covariates were renumbered as:

B 1 = {x 1 , ..., x 4 } ; B 2 = {x 5 , ..., x 10 } ; B 3 = {x 11 , ..., x 20 } ; Noise = {x 21 , ..., x 40 }
The linear predictor was calculated as the following function of latent variables f j :

∀(i,t, j) : η it = .25 + ξ 1 it -.5ξ 2 it (21)
Hence the hazard function of each subject i: h i (t) = h 0 (t)e η it .

Note that : 1) Within every variable-bundle B j , correlation is very high (0.96), which precludes standard Cox regression.

2) The bundle playing the major role is B 1 , which is only the third in size, and the latent variable structuring it, ξ 1 , has a positive coefficient.

3) The bundle playing the secondary role is B 2 , which is the second in size, and the associated latent variable, ξ 2 , has a negative coefficient. 4) ξ 3 does not play any part in the hazard.

Associated bundle B 3 is thus a decoy-structure, and it is the first in size, weighing the same as B 1 and B 2 together. 5) Pure noise outnumbers every bundle, weighing the same as all bundles together.

The challenge is to identify B 1 as the primary explanatory structure, with positive effect, and B 2 as the second, with negative effect.

Independent right-censoring was simulated according to an exponential distribution (constant hazard rate) so that 21 subjects were censored before undergoing the event, and 54 underwent the event.

Estimation results

The parameter-varying scheme we follow is the one proposed in section 2.6:

• We start with s = 1 (taking no account of the goodness-of-fit), and see how gradually taking the goodness-of-fit into account by decreasing s changes the results.

• Parameter l is regarded as secondary tuning-parameter which may improve interpretation of components when raised, once a reasonably good value of s has been found.

• Taking a low value (e.g. 0 < τ < .05) is enough for pure exploratory purposes, but stable prediction requires higher values. Starting with τ = 10 -2 merely to ensure identification in presence of high collinearities, we vary τ only in the prediction-stage.

Let us now see the evolution of the results when we follow this parametervarying scheme. Decreasing s to .95 has a dramatic effect: ξ 1 captures important parts of B 1 and B 2 with opposite effects (4). But B 3 still has a non-negligible projection in plane ( f 1 , f 2 ). Raising l to 4 allows the plane to focus on B 1 and B 2 exclusively (fig. 6). Now, fig. 7 shows that f 3 aligns B 3 perfectly, and that f 4 aligns noise-variable x 39 , which is obvious over-fitting, yet given away by a high p-value. We witness on it the opposite roles of B 1 and B 2 , but the component is not so close to the direction of either, and is likely to contain over-fitting. Component f 2 is not really determined. Starting to take some structural relevance into account by raising s to .1 brings back B 3 along component f 2 but does not change much f 1 , essentially determined by the maximum likelihood estimation.

The regularization step

Once we retain a suitable (s, l) pair, e.g. here (.95, 4), we can tune τ so as to regularize the coefficient-vector, with intent to lessen fitting on noise, so distribute coefficients more evenly among highly correlated variables. Using the estimated coefficients of the first two components calculated every time, we get the estimated linear predictor η, and calculate its correlation with the actual simulated η. The results are given in table 1 and show regularization Table 1: Coefficients of variables in components f 1 , f 2 and correlation of estimated predictor η with actual one η for increasing τ (s = .95, l = 4). 4 Application to real data

τ = 0 τ = .1 τ = .3 τ = .5 τ = .7 f 1 f 2 f 1 f 2 f 1 f 2 f 1 f 2 f 1 f 2 B 1 x 1 -.
The supervized component Cox regression method will now be used to study the advent of polygamy among men living in Dakar.

Poverty in Senegal is high. Dakar gathers about one fourth of the senegalese population and attracts an important migration. The society in Dakar is moslem and women are dominated by men to a notable extent. The 2002 census revealed that among men, polygamy is on the increase mostly from the age of 45. Among married men living in Dakar aged 45 to 49, the percentage of the polygamous is over 20%. Among those aged 60 and over, this percentage is rather stable, ranging between 30% and 40%.

The mechanisms making polygamy possible have been studied and are well-documented [START_REF] Pilon | Contribution à l'analyse de la polygamie[END_REF]. The existence of polygamy rests upon an important age-gap at marriage between men and women: women enter marriage much younger than men. The frequent and quick remarrying of divorced women and widows also favors this practice. Among other things, polygamy allows the man to maximize his posterity [START_REF] Chojnacka | Early marriage and polygyny : feature characteristics of nuptiality in Africa, Genus[END_REF].

Having many descendants means benefitting from an important manpower.

It also means hoping for a better assistance in the old age: in case of health problems, it raises the chances to be looked after by at least one wife and her children [START_REF] Møller | Polygamy, economic security and well-being of retired Zulu migrant workers[END_REF][START_REF] Gning | Polygamie et personnes âgées au Sénégal[END_REF]. In that, polygamy may serve as a safety net. Finally, polygamy is, for a man, a way to reconcile the preferences of the group (social choice) when his parents impose a spouse on him, with his own (individual choice).

We are going to confront these elements with the results of a first and somewhat naive application of supervized component Cox regression. Our intention here is not to complete a thorough analysis of the data, but to Cross-validated supervised components can be used as such a criterion.

To illustrate the gain provided by supervised-component Cox regression,

we first performed Cox regression on principal components, which is one end of the continuum spanned by SCCoxR, associated with s = 1. Then, to the other end of the continuum (s arbitrarily close to 0), we performed

Cox regression on the whole explanatory subspace, in order to see whether the linear predictor is close to an interpretable dimension. In view of the findings provided by these two attempts, we finally tune the parameters of SCCoxR so as to track down interpretable explanatory dimensions.

Cox regression on the first 8 principal components

Cox regression on principal components is obtained setting s = 1 and l = 1.

Cox regression on the full explanatory space

To the other end of the continuum is value s = 0. Setting s = 0 is forbidden here because of the exact linear dependence of some variables, e. once by the first one. Since we took a small s, component 2 is identified, but almost random, and plays no role in the prediction. By taking τ = 1, we ensured regularisation of the linear predictor, i.e. that confusion be minimized between correlated variables. So, the signs of coefficients reflect better the effects of the associated variables, provided their absolute value is not too low. These coefficients are given in line β of table 3.

They should a priori not be considered completely reliable, insofar as they are not based on strong structural dimensions of the explanatory variables.

Neither of these extreme choices has provided substantial insight as to the main dimensions influencing the risk. Let us track these by fine-tuning the parameters, as in the simulations: first s, then l, and finally τ.

Fine-tuning supervized component Cox regression

As in the simulations, we kept l = 1 and started decreasing s from s = 1 to s = 0.95, then 0.9, 0.85 and 0.8. The best results were obtained for s = 0.9, with p-values (viewed as descriptive indicators) of components 2 to 5 lower than 1e -5. Then, in order to improve interpretability of the supervised components by dragging them towards close variable-structures, we raised l from 1 to 8.

Cox regression on the obtained first 10 supervized components shows that components 2 to 5 have a p-value lower than 1e-4, and component 8 has a p-value lower than 1e-2. In order to assess the real predictive power of the components, we perform cross-validation on models built with a number of We may now take a look at some correlation scatterplots associated with 6 components and try to interpret the space they span in terms of explanatory sets out the secondary education, to which it is highly positively correlated.

It is negatively, but poorly correlated with low education levels. Figure 15 shows an interesting plane, spanned by components 3 and 4. Its interpretation is wholly connected to places of birth and infancy, with a triangle opposing the three factor-levels: Dakar ("Dak"), other city ("urb") , rural zone ("rur"). As in most situations where the variable-vectors illustrating the components sum up to zero (this is the case of the three birth-place levels, and of the three infancy-place ones), variables don't align components really well, but the plane they span is clearly related to the risk. Finally, component 5 aligns the age-gap, revealing its role in polygamy, and component 6 is loosely related to the kinship of the first spouse to ego.

At this stage, we may use the components to try and select promising subsets of explanatory variables for a standard Cox regression, bearing in mind that doing so, we may lose valuable predictive information, since predictive components combine many more variables to be efficient. The estimated effects of some variable-subsets, with corresponding p-values, are given in table 2.

Supervised components do seem to help variable-selection. Yet, one can see that some confusion remains: the p-values of effects sometimes change dramatically when a variable (e.g. "Infancy in Dakar") is replaced with a Components are combinations of variables, and the proper way to use them is as follows: select the best number of supervized components through cross-validation; then, as they are regularized combinations of variables, use them to calculate the corresponding β , and interpret the signs and magnitudes of the largest effects. In the table displayed on fig. 3, line β 5 displays the coefficients calculated on the basis of the first 5 components. These are the ones to be used for interpretation of variable-effects, because the corresponding linear predictor is based on the best predictive component-model.

The linear predictor has been regularised so as to minimise confusion between effects of correlated variables. One consequence of this is that the signs of the betas can overall be trusted, at least if the absolute value is not negligible.

Bootstrapping can of course provide confidence intervals associated with the coefficients, but it is computationally very costly.

From coefficients in table given on fig. 3, it seems possible to calculate relative risk ratios. For instance, the risk-ratio of a man having spent his infancy in the country relative to a man having the same characteristics except having spent his infancy in Dakar would be: But one should never forget that this "having the same other characteristics" is sheer fiction: in the situation being dealt with, variables may be linked up to perfection. Picture a variable denoted a, linked to no other, and having coefficient β a in the estimated model. Now, picture duplicating a into a 1 and a 2 in a new model, with all other variables unchanged. Then, regularization will entail that the former effect of a be distributed evenly on a 1 and a 2 , so that their coefficients in the estimated new model will both be equal to β a 2 . So, one must keep in mind the statistical links between explanatory variables before trying to interpret the magnitude of these relative risk ratios.

For a binary variable without its complement, or for a quantitative variable, the rule would be the same as in Cox regression. For instance: the risk-ratio of a man having given his consent for his first marriage relative to a man not having, but sharing all other characteristics, would be:

r(Consent / no Consent) = exp(β 5 Consent ) = 0.93 (23)
And the effect of an extra-year in age-gap, all other characteristics equal, would be to multiply the risk by: r(AgeGap) = exp(β 5 AgeGap ) = 0.66 (24)

We shall now interpret the signs and sometimes the magnitudes of the coefficients of variables linked to the first dimensions found by supervisedcomponent Cox regression. We chose to include component 6 because it is illustrated by a kinship variable, and even if on average it tends to produce overfitting, there are many cross-validations samples with which it contributes to improve the CV coefficient. As a general rule, we advise to possibly consider a few extra-components after the optimum, provided the fall in CV keeps moderate and such components are reasonably correlated to some variables.

Linked to component 1 are the offspring variables. On the whole, the effect of the number of children on the hazard-rate is negative. Having no son appears to have a special effect, increasing the risk with respect to all other numbers of sons. This is not the case with daughters.

Linked to component 2 are the ego-education variables, and more specifically secondary education. Indeed, the coefficients show that its effect relative to all other levels is to decrease the risk. For example, the relative risk ratio of a man with secondary education with respect to a man having no Table 3: Coefficients of variables in the regularized linear predictor obtained for: s = 1e -3, l = 1, τ = 1, 1 component (β ), and s = .9, l = 8, τ = 1, 5 components (β 5)

.21e-1 -1.33e-1 1.76e-1 1.02e-1 -2.01e-1 -1.54e-1 -2.04e-2 -4.81e-2 β 5 -9.08e-2 4.29e-2 1.41e-1 -8.54e-2 2.21e-1 1.34e-1 -1.76e-1 -3.00e-1 -3.79e-3 -3.
-1.74e-2 1.74e-2 1.21e-1 -5.27e-2 1.47e-1 -2.21e-1 -1.38e-1 β 5 -3.47e-2 3.47e-2 1.96e-1 2.50e-2 1.37e-1 -3.81e-1 -1.55e-1
education, and all other characteristics identical, would be:

r(Edu Sec/Edu Non) = exp(β 5 Edu Sec -β 5 Edu Non ) = exp(-0.143 -0.0635) = 0.81 (25)
Men having had a secondary education or above tend to enter marriage later and with a partner they chose. Their conception of marriage is less traditional. This may account for a lesser risk of becoming polygamous.

Linked to components 3 and 4 are the places of birth and infancy. Of course, these variables being strongly linked, their effect on the risk will tend to be shared and distributed evenly across very correlated ones. As a consequence, one cannot for instance speak of the effect of "infancy in Dakar" with respect to "infancy in a rural zone" with all other characteristics identical, because infancy is much linked to place of birth, so that their effect is common to a large extent. Where standard Cox regression, as any non-regularized regression, tries to separate the effects of explanatory variables -thus having severe problems when these are not separable -supervised-component Cox regression does all the less so as variables are correlated. The effects of uncorrelated explanatory variables will be their "proper" effect indeed, just as in plain Cox regression, but the effects of highly correlated explanatory variables will be an even distribution across them of their "common" effect. How do we make use of that here? By considering the bundle { Inf rur , BP rur } as a whole and opposing it to { Inf Dak , BP Dak }, for instance. The relative risk ratio of a man born and raised in a rural zone with respect to one born and raised in Dakar, with all other characteristics identical, would then be: = exp(0.052 + 0.058 -0.011 + 0.12) = 1.29 (26)

Likewise, the relative risk ratio of a man born and raised in a city other than Dakar with respect to one born and raised in Dakar, with all other characteristics identical, would be: 1.17.

Of course, if the variables are not strictly linked through a logical constraint, it is always conceivable to vary the situation. For instance, the relative risk ratio of a man born and raised in a rural zone with respect to one born in a rural zone but raised in Dakar, with all other characteristics identical, would then be: down the risk by one third. This may be linked to the "safety net" approach to polygamy, namely the man's wish to increase his offspring when his first wife, having an age close to his, is getting old or unwilling to, or his wish to found a younger family to support him in his old age.

Component 6 is related to the kinship of ego to his first wife. The relative risk ratio of a man having a first wife unrelated to him with respect to a man having one related through his mother is: exp(β 5 ParSp1 Noβ 5 ParSp1 Mat ) = 0.62 , and the ratio is 0.66 when the first spouse is related to ego's father.

A kinship between spouses is very often the sign of a prescribed marriage, which may easily lead to polygamy because, if it is difficult to part from a spouse more or less imposed by the family, ego can make for it by marrying a second wife more freely chosen, most often unrelated.

It must be emphasized that we have given these ratio calculations for the mere sake of illustration. The reader will naturally be aware that such calculations should not be taken at face value. They merely indicate tendencies obtained here on a rather small sample, and always under the costly hypothesis of "all other characteristics kept identical".

This toy-application on real data allowed us to see what supervized component Cox regression brings to the analyst in terms of explanatory dimen-sion reduction, but we only demonstrated the first step of the analytical process. On a second step, the selected variables should be taken as extracovariates Z, and SCCoxR relaunched so as to track down more complementary predictive dimensions, if any. And so on until the predictive dimensionality is exhausted.

Conclusion

The Yet, staying "too close" to the current starting point on the arc (v [t] , γ [t] ) may make the algorithm too slow to reach the maximum. To avoid that, we propose to use a Gauss-Newton unidimensional maximization to find the maximum of g(v) on the arc (v [t] , γ [t] ), and take it as v [t+1] .

The fixed point of the resulting algorithm is a critical point of (1), hence a local maximum of g s.t. C v = 0 .

Figure 1 :

 1 Figure 1: Polar representation of the variable powered inertia according to the value of l

  For a given model M , the cross-validation quality coefficient (CV) is calculated according to the technique proposed for the proportional hazard model by van Houwelingen et al. (2006): the sample is split into K parts, and for each part k, we calculate:

Figures 2

 2 Figures 2 to 7 show the correlation-scatterplots of variables with com-

Figure 2 :Figure 3 :Figure 4 :Figure 5 :Figure 6 :

 23456 Figure 2: Correlation scatterplot a : s = 1, l = 1 f 1 : coefft = -.03 ; p=.83 ; f 2 : coefft = -.42; p=.004

Figure 7 :

 7 Figure 7: Correlation scatterplot f : s = .95, l = 4 f 3 :coefft = .06 ; p=.70 ; f 4 : coefft = .18; p=.20

Figure 9 :Figure 10 :

 910 Figure 9: Correlation scatterplot g : s = 0 f 1 : coefft = -1.85 ; p<2.e-16 ; f 2 : coefft = -.12; p=.35

  provide a first illustration of how the supervised-component Cox regression method works, how it should be used, and the kind of help it may give to the data analyst. Supervised-component Cox regression is applied to life-histories of 222 married men born before 1967 and living in Dakar. The data comes from the 2001 retrospective survey conducted by Philippe Antoine (IRD) and Abdou Salam Fall: (IFAN-UCAD) "Crise, passage à l'âge adulte et devenir de la famille dans les classes moyennes et pauvres à Dakar"[START_REF] Antoine | Crise, passage à l'âge adulte et devenir de la famille dans les classes moyennes et pauvres à Dakar[END_REF], partially funded by the CODESRIA. Our aim is to find out the main explanatory dimensions that influence the risk of becoming polygamic by marrying a second wife. The survival time is the duration of the first union in the monogamic state. The terminal event is marrying a second wife. Death of the first wife and divorce were considered independent right-censoring, just as was the date of the survey. Of course, divorce is likely not to be independent of the event under study, so should be treated as a competing event in future work. The subjects are described with 107 time-varying numerical variables, most of which are indicators of nominal variable values.All are put in matrix X, matrix Z being null here. The total number of rows of the data matrix (number of (subject, date) pairs) is 2293, and 55 events were observed. The high number of variables and their collinearities make direct unregularized Cox regression impossible. Selection of variables by means of mere intuition is tricky and may be opposed on the grounds that the many confusion effects change the results according to the variables selected, so, at least, variable-selection should lean on some solid criterion.

  Figure 11: Correlation scatterplot with principal components (4,5)

Figure 13 :

 13 Figure 13: CV coefficient against number of components for s = 0.9, l = 8, τ = 1

Figure 14 :Figure 16 :

 1416 Figure 14: Correlation scatterplot with supervised components (1,2) for s = 0.9, l = 8, τ = 1

  r(Inf rur/Inf Dak) = exp(β 5 Inf rurβ 5 Inf Dak ) = exp(.059 + .12) = 1.20 (22)

  r(Rural/Dakar) = exp(β 5 BP rur + β 5 Inf rur -(β 5 BP Dak + β 5 Inf Dak ))

r(

  Inf rural/Inf Dakar) = exp(β 5 Inf rurβ 5 Inf Dak ) = 1.20 (27) Component 5 is very strongly correlated to the age-gap (age of ego -age of first wife). The corresponding relative risk has been given above: every extra-year in age-gap, all other characteristics kept equal, would roughly cut

  Figure 17: PING vectors

  For the retained value of s, start raising l gently (integer values are sufficient): l ∈ {1, 2, 3, 4, ..., 10, ...}. Rather soon, one notices on the scatterplots that if the interpretation of components was enhanced by an early raise

	The secondary key-parameter to be tuned in this stage is l (locality of
	variable-bundles components should align on).
	(c)
	of the potential modeling value of the component. Keep a value of s giving
	interpretable components with reasonably low p-values.

and choose H so that principal components with rank > H may be considered noise. Perform Cox's regression on the H principal components.

These being exogenous, the standard significance tests are valid. Downsize H to the highest rank of a significant component, which is its maximum possibly interesting value.

(b) Perform supervised-component Cox regression with a decreasing sequence of s values starting slightly below s = 1, e.g. 0.95 and using a reasonable step (e.g. 0.05). This is gradually letting the goodness-of-fit into the criterion. For each value of s, use the p-values of Cox's regression on the components as a mere descriptive indicator of their roles. These p-values are not taken here as rigorous inference values, which they aren't because components are not exogenous when s < 1 , but as "descriptive indicators" of l , raising it further only changes them negligibly. For the few values retained of l, examine the p-values of components and see that an easier interpretation is not paid for by some non-negligible rise (it may even decrease them).

Predictive phase:

We have now a good idea of the values to give to s, l, H. These mays still change slightly with the value of τ (tuning regularization), but not dramatically.

(d) Start raising τ gently (e.g.: τ ∈ {.1, .2, .3, .4, ..., 1} ). For each value, calculate the CV coefficient for models ranging from h = 1 to H components and select the number of components H * (τ) giving the best CV, then denoted CV * (τ).

Table 2 :

 2 Results of Cox regression on three variable-subsets selected through correlation with supervised components

	Variable	β	exp(β )	p
	Son 0	0.6590 1.9329 0.0553
	Edu Sec	-0.6567 0.5186 0.0315
	Inf Dak	-0.5452 0.5797 0.0697
	AgeGap	-0.0970 0.9075 0.0012
	ParSp1 No -0.9242 0.3968 0.0016
	Variable	β	exp(β )	p
	Son 0	0.7193 2.0531 0.0362
	Edu Sec	-0.6468 0.5237 0.0343
	BP Dak	-0.3060 0.7364 0.3058
	AgeGap	-0.0975 0.9071 0.0014
	ParSp1 No -0.9188 0.3990 0.0016
	Variable	β	exp(β )	p
	NumChild -0.1309 0.8773 0.0884
	Edu Sec	-0.6868 0.5032 0.0250
	Inf Dak	-0.6182 0.5389 0.0371
	AgeGap	-0.0964 0.9081 0.0012
	ParSp1 No -0.9459 0.3883 0.0012

 et al. , 2016), and application to counting processes is currently considered.

Appendix A : The gradient formulas

On every time-point t i when an event occurs for some subject i, δ i = 1 , we have the following partial likelihood term:

e αx it i u+z it i γ ∑ j∈R t i e

The overall partial likelihood is :

with:

61

Putting ∀t i , ∀ j ∈ R t i : ω j,t i := e αx jt i u+z jt i γ and Ω t i := ∑ j∈R t i ω j,t i , we get:

And finally:

x it i -1

• Variable powered inertia:

62

Letting Ω(u) = diag(ϖ l-1 j ) , we have:

Appendix B : The Projected Iterated Normed Gradient (PING) algorithm

Notation: the current value of any quantity a on iteration t is denoted:

Putting v = M -1/2 u , g(x) = h(M 1/2 x) and C = M 1/2 D , this is strictly equivalent to:

The corresponding lagrangian is:

Put back into (4), this yields:

In the particular case where C = 0 , we shall take:

Finally, ( 5) and (1) imply:

This gives the basic iteration of the PING algorithm:

65

Let us show that this iteration follows a direction of ascent. Because, by construction: ∀s: v [s] ⊥ C, we have:

which has the sign of:

Picking a point on a direction of ascent does not guarantee that g actually increases, for we may "go too far" in this direction. Let γ

. Staying "close enough" to the current starting point on the arc (v [t] , γ [t] )

guarantees that g increases. Indeed, let ϖ be the plane tangent to the sphere on v [t] and let w denote the vector tangent to arc (v [t] , γ [t] ) on v [t] (cf. fig. 17).

Then:

∃τ > 0, w = τΠ ϖ γ [t] ⇒ w|γ [t] = τ Π ϖ γ [t] |γ [t] = τ cos 2 (γ [t] , ϖ) > 0 (55)