Gabriel Capson-Tojo

Chloé Girard

Maxime Rouez

Marion Crest

Steyer

Nicolas Bernet

Jean-Philippe Delgenès

Renaud Escudié
email: renaud.escudie@inra.fr

Biochar and industrial FeCl 3 as additives to favor the consumption of volatile fatty acids 1 during anaerobic digestion of food waste

Keywords: Biogas, Anaerobic Digestion, Char, Design of Experiment

BACKGROUND 15

Although anaerobic digestion is a promising alternative for the valorization of complex 16 substrates such as food waste, this process is yet to be optimized and options for its 17 stabilization must be developed. The goal of this study was to assess for the first time if the 18 addition of biochar, together with trace elements in the form of an industrial FeCl 3 solution, 19 could serve as method for stabilizing anaerobic digestion of food waste. 20

RESULTS

21

Results from batch reactors demonstrated that the supplementation of both biochar and 22 industrial FeCl 3 favored the digestion kinetics. Their addition improved the maximum 23 methane production rates (from 886 up to 1498 ml•g VS -1 •d -1 ) and the average daily methane 24 production rates (from 280 up to 376 ml•d -1 ), related to acetate and propionate consumption, 25

Introduction 35

The increasing production of food waste (FW) and novel international regulations call for the 36 development of novel sustainable technologies for its treatment and valorization [START_REF] Capson-Tojo | Food waste 396 valorization via anaerobic processes: a review[END_REF]. Among 37 the existing options, anaerobic digestion (AD) is an environmental-friendly process that 38 provides an efficient waste treatment, producing at the same time biogas and digestate. [START_REF] Cooney | Start up performance of biochar 504 packed bed anaerobic digesters[END_REF] Nevertheless, FW AD is a complex biological process that often leads to inefficient results or 40 reactor acidification. The first complication occurring during FW AD is related to the fast 41 biodegradability of this substrate, which is mainly composed of easily degradable 42 carbohydrates. Therefore, the reactors can be easily overloaded, especially in batch systems 43 [START_REF] Capson-Tojo | Dry 399 anaerobic digestion of food waste and cardboard at different substrate loads[END_REF], causing an unbalance between the acidogenesis/acetogenesis and the methanogenesis 44 steps, which results in an initial accumulation of volatile fatty acids (VFAs) and a consequent solution may be industrial FeCl 3 , which consists on an acid solution, highly concentrated in 76 different metals (mainly Fe) that is commonly applied in wastewater treatment plants and 77 anaerobic digesters for pH control worldwide. Concerning GAC, an affordable substitute may 78 be biochar, which is also a carbon-conductive material with the capacity of accepting and 79 donating electrons [START_REF] Prévoteau | The electron donating capacity of 481 biochar is dramatically underestimated[END_REF]. In addition, biochar may also improve the AD process by providing 80 buffering capacity and by sorption of inhibitors [START_REF] Fagbohungbe | The 483 challenges of anaerobic digestion and the role of biochar in optimizing anaerobic digestion[END_REF]. It also favors nutrient retention in the 81 digestates, facilitating nutrient uptake if the digestate is spread on land for plant cultivation 82 [START_REF] Fagbohungbe | The 483 challenges of anaerobic digestion and the role of biochar in optimizing anaerobic digestion[END_REF][START_REF] Wang | Improving anaerobic 486 digestion of easy-acidification substrates by promoting buffering capacity using biochar 487 derived from vermicompost[END_REF]. Moreover, as biochar can be produced from digestate [START_REF] Monlau | A new concept for 489 enhancing energy recovery from agricultural residues by coupling anaerobic digestion and 490 pyrolysis process[END_REF] or directly from FW [START_REF] Rago | Assessing the potential of biofuel (biochar) 492 production from food wastes through thermal treatment[END_REF], 83 its addition for AD improvement clearly fits within the approach of environmental 84 biorefinery. 85

Biochar has been previously used as amendment for AD: (i) for biogas purification (removal 86 of H 2 S or CO 2 ) [START_REF] Linville | In-situ 495 biogas upgrading during anaerobic digestion of food waste amended with walnut shell biochar 496 at bench scale[END_REF][START_REF] Shen | Towards a sustainable paradigm of waste-to-energy process: Enhanced anaerobic digestion of 499 sludge with woody biochar[END_REF][START_REF] Shen | Producing 501 pipeline-quality biomethane via anaerobic digestion of sludge amended with corn stover 502 biochar with in-situ CO2 removal[END_REF], (ii) as reactor packing for biofilm support [START_REF] Cooney | Start up performance of biochar 504 packed bed anaerobic digesters[END_REF], (iii) as AD substrate 87 [START_REF] Mumme | Use of biochars in anaerobic digestion[END_REF], (iv) for nutrient supplementation [START_REF] Shen | Yearlong semi-continuous 508 operation of thermophilic two-stage anaerobic digesters amended with biochar for enhanced 509 biomethane production[END_REF], (v) as matrix for sorption of inhibitors [START_REF] Fagbohungbe | Impact of 511 biochar on the anaerobic digestion of citrus peel waste[END_REF][START_REF] Torri | Biochar enables anaerobic digestion of aqueous phase from 514 intermediate pyrolysis of biomass[END_REF] or 88 (vi) as means to increase the buffer capacity of the system [START_REF] Wang | Improving anaerobic 486 digestion of easy-acidification substrates by promoting buffering capacity using biochar 487 derived from vermicompost[END_REF]. However, only three studies 89 have been carried out dealing with low-dosage of biochar for improving the AD kinetics. Luo 90 et al. [START_REF] Luo | Application of eco-compatible biochar in anaerobic 516 digestion to relieve acid stress and promote the selective colonization of functional microbes[END_REF] studied the addition of biochar for improving the kinetics of methane production 91 and VFA degradation (HAc and HBu) with glucose as substrate, concluding that a positive 92 effect existed. Zhao et al. [START_REF] Zhao | Potential 467 enhancement of direct interspecies electron transfer for syntrophic metabolism of propionate 468 and butyrate with biochar in up-flow anaerobic sludge blanket reactors[END_REF] investigated the effect of biochar and ethanol supplementation 93 in reactors degrading HPr and HBu, observing that the degradation of both VFAs was 94 improved by adding biochar and concluding that bacteria known to participate in DIET (i.e. 95

Geobacter species) were attached onto the biochar surface. Finally, Sunyoto et al. [START_REF] Sunyoto | Effect of biochar addition on hydrogen 519 and methane production in two-phase anaerobic digestion of aqueous carbohydrates food 520 waste[END_REF] 96 assessed the influence of biochar addition in a two-phase AD reactor treating aqueous 97 carbohydrates FW, concluding that biochar supplementation increased the hydrogen and 98 methane production rates. They attributed this improvement to a promotion of the biofilm 99 formation. 100

F o r P e e r R e v i e

w

To the knowledge of the authors, no study has been carried out so far to assess if addition of 101 biochar can favor VFA consumption during AD of complex substrates such as FW. Moreover, 102 the stabilizing effect of jointly adding industrial FeCl 3 (to provide TEs) and biochar on the 103 performance of FW AD has been investigated for the first time. Therefore, the goal of this 104 study was to: (i) optimize the concentrations of both biochar and industrial FeCl 3 using a 105 batch experimental design to improve the kinetics of FW AD and (ii) evaluate the effect of 106 industrial FeCl 3 and biochar supplementation on the AD performance of semi-continuous 107 pilot reactors treating commercial FW. 108 109 Experimental 110

Inoculum, substrate and AD additives 111

The reactors were inoculated with digestate from an industrial plant treating different organic 112 streams at high TAN/FAN concentrations (5.04 g TAN•l -1 ; 615 mg FAN•l -1 ). Thus, it was 113 assumed that the microorganisms would be adapted to high the high TAN/FAN levels 114 associated with FW AD [START_REF] Capson-Tojo | 417 Accumulation of propionic acid during consecutive batch anaerobic digestion of commercial 418 food waste[END_REF]. The TS content of the inoculum was 5.81 ± 0.02 %, with a 115 proportion of volatile solids (VS) of 59.13 ± 0.08 %. The activity of the inoculum was 116 verified prior its utilization using ethanol as substrate. The commercial FW was collected 117 from five mayor producers from the region of the Grand Narbonne, in the south of France. A 118 proportional mixture (wet weight) was used as substrate. The main characteristics of the FWs, 119 as well as those of the mixture used as substrate and the inoculum are presented in Table 1. 120

The characteristics of the FW mixture were in agreement with typical values presented in the 121 literature [START_REF] Capson-Tojo | Food waste 396 valorization via anaerobic processes: a review[END_REF], which indicated that it could be considered as a representative sample of a 122 general FW. It had a TS contents of 21 % (90.3 % VS), it was mainly composed of 123 carbohydrates (618 g•kg TS -1 ) and it had relatively low C/N ratios (16.1). No sulfur was 124 detected. A more extensive characterization and a deeper discussion of the results can be [START_REF] Capson-Tojo | 417 Accumulation of propionic acid during consecutive batch anaerobic digestion of commercial 418 food waste[END_REF]. 126

The industrial FeCl 3 was provided by an industrial AD plant, where it is still used for pH 127 controlling purposes. The characteristics and composition of this industrial solution, rich in 128 TEs, are shown in Table 2. 129

The biochar was natural slow-pyrolyzed wood charcoal, commonly found in the market for 130 several applications. Before utilization, the biochar was grinded and sieved at 600 µm. 131

Batch experimental design for dosage optimization 132

A multilevel factorial design was used to optimize the dosage of biochar and industrial FeCl 3 133 and to evaluate their individual effect on the methane production and the VFA production-134 consumption kinetics. Digestate from continuous reactors digesting FW was used as inoculum 135 (after consumption of the remaining VFAs). This sludge had a TS content of 5.17 %, with 136 60.2 % corresponding to VS and had a TAN content of 7.27 g•l -1 . Sixty g of FW were added 137 as substrate in all the reactors. A substrate to inoculum (S/X) ratio of 1 g VS•g VS -1 was 138 applied, leading to initial FW concentrations of approximately 27 g VS FW•l -1 . The reactors 139 were incubated at 37 °C. The working volumes ranged from 487 ml to 529 ml. Two different 140 concentrations of the FeCl 3 solution (0.1 and 0.2 g Fe•l -1 ) and three of biochar (10, 55 and 100 141 g•l -1 ) were tested. As a result, an experimental design with 12 conditions was defined (6 142 conditions in duplicate; Table 3). Two consecutive batch feeding were carried out. The first 143 feeding served for inoculum adaptation (results not presented), and the results of the second 144 one were used for modelling purposes. It must be mentioned that, based on previous results, it 145 was decided not to include into the experimental design reactors containing only biochar and 146 was also used to define the levels of industrial FeCl 3 and biochar to be tested (30). 151

This experimental design was chosen because it allows analyzing the effects of the selected 152 factors (biochar and FeCl 3 concentrations) on a chosen output through the entire experimental 153 region covered. This implied that it was possible to distinguish between the effect of biochar 154 and FeCl 3 without the need of single additive treatments. It allows so by predicting the 155 responses using a quadratic model (Eq. 1): 156 157

ݕ = ܽ + ∑ ܽ ିଵ • ݔ + ∑ ܽ ିଵ • ݔ ଶ + ∑ ܽ ழ • ݔ • ݔ Equation 1 158 159
Where y is the response to be predicted, x i are the studied factors and a i , a ii and a ij are the 160 parameters corresponding to each factor. These parameters represent the linear effects, the 161 quadratic effects and the interactional effects, respectively. The first coefficient a 0 is required 162 for fitting the mathematical model. The p-values from F-tests (95 % confidence) and the 163 coefficient of determination R 2 were used to evaluate the fitness of the model. The experiment 164 was designed and evaluated using the software STATGRAPHICS Centurion XVI Version 165 16.1.03 (©StatPoint Technologies Inc.). 166

The reactors were incubated in an Automated Methane Potential Testing System (AMPTSII) 167 (Bioprocess Control, Sweden). Twelve reactors from the AMPTSII system, with a total 168 volume of 600 ml, were used. According to the manufacturer instructions, they were 169 connected to CO 2 traps (NaOH solutions) and to gas flow meters to determine continuously 170 the methane flow rates. These reactors also allowed the follow-up of the dynamics of VFA 171 production-consumption. A hole present in each vessel was used as sampling port, facilitating 172 an easy sampling of the digestate. The reactors were agitated during one minute every 10 173 minutes at 40 rpm. 174

The criterion followed to decide when the batch had finished was the total consumption of the

Continuous pilot scale reactors 178

In parallel to the batch reactors, two different pilot scale reactors were run: a control system 179 simply digesting FW and a doped reactor supplemented with biochar and the FeCl 3 solution. 180 Both reactors were incubated at 37 °C and had a working volume of 12 l. The reactors were 181 fed once per day, initially with an OLR of 1.4 g VS•l -1 •d -1 , corresponding to a hydraulic 182 retention time (HRT) of 110 days. This value was increased to 2.8 g VS•l -1 •d -1 (HRT of 55 183 days) after 77 days of operation. An equivalent amount of digestate was withdrawn to keep 184 the volume of the reactors constant. The FeCl 3 solution was diluted with water (x20 vol:vol) 185 and dosed into the supplemented reactor to keep a constant concentration of 100 mg Fe•l -1 186 (value calculated from optimal results reported in the literature [START_REF] Zhang | Enhanced anaerobic digestion of food waste by trace metal 420 elements supplementation and reduced metals dosage by green chelating agent [S, S]-EDDS 421 via improving metals bioavailability[END_REF][START_REF] Banks | Trace element requirements for stable food 423 waste digestion at elevated ammonia concentrations[END_REF](12)). This corresponded to 187 concentrations of 0.38 mg Mn•l -1 , 0.19 mg Zn•l -1 , 0.03 mg Ni•l -1 , 0.01 mg Co•l -1 and 0.01 mg 188 Mg•l -1 in the reactor. The initial concentration of biochar in the supplemented reactor was 10 189 g•l -1 . As it will be further explained, this concentration was increased up to 50 g•l -1 to favor the 190 consumption of the accumulated VFAs, according to the results obtained from the batch 191 experimental design. 192

The pilot reactors consisted of jacketed cylindrical vessels made of stainless steel that had 193 inner stirring blades to provide continuous agitation. A more detailed description of the 194 reactors can be found elsewhere [START_REF] Capson-Tojo | 417 Accumulation of propionic acid during consecutive batch anaerobic digestion of commercial 418 food waste[END_REF]. The experiments lasted for 196 days. 195

Analytical methods 196

Physico-chemical characterization of the FW 197

The TS and VS contents were measured as described in the standard methods of the American 198 Public Health Association [START_REF] Apha | Standard Methods for the Examination of Water and Wastewater[END_REF]. The concentrations of carbohydrates and lipids were 199 determined using the Dubois method [START_REF] Dubois | Colorimetric Method for 524 Determination of Sugars and Related Substances[END_REF] and by a gravimetric method based on accelerated 200 [START_REF] Apha | Standard Methods for the Examination of Water and Wastewater[END_REF], respectively. The protein content was calculated using the total 201 Kjeldahl nitrogen (TKN) contents (6.25 g protein•g N -1 (48)). The TKN and NH 4 + 202 concentrations were determined using an AutoKjehdahl Unit K-370, BUCHI. The contents of 203 organic (TOC) and inorganic carbon (IC) were measured with a Shimadzu (Kyoto, Japan) 204 TOC-V CSN Total Organic Carbon Analyzer coupled to a Shimadzu ASI-V tube rack. The C/N 205 ratio was calculated as TOC/TKN. A WTW (London, United Kingdom) pHmeter series 206 inoLab pH720 was used for pH measurement. Finally, the biochemical methane potentials 207 (BMPs) of the substrates were measured according to Motte et al. [START_REF] Motte | Morphological structures of wheat straw strongly impacts its anaerobic digestion[END_REF]. 208

F o r P e e r R e v i e w solvent extraction

Gas quantification and analysis 209

The amount of methane produced was automatically measured in the AMPTSII system and 210 the volume of biogas produced in the pilot reactors was continuously measured using Ritter 211 MilliGascounters MGC-1 V3.0. The composition of the biogas was determined by gas 212 chromatography as described in Cazier et al. [START_REF] Cazier | Biomass hydrolysis inhibition at high 532 hydrogen partial pressure in solid-state anaerobic digestion[END_REF]. 213

Analysis of metabolites and final products of the digestion 214

A sample of digestate from the pilot reactors was taken once per week for measurement of the 215 concentrations of VFAs and ionic species in the reactors. Concerning the batch study, a 216 plastic tube submerged in the reaction media served for digestate sampling when required. 217 Before sampling, a clip was used for blocking the gas output and the equivalent volume of 218 digestate to be removed was injected as nitrogen gas, avoiding an overestimation of the 219 produced gas. The concentrations of VFAs and ionic species in the digestates were measured 220 as described in Motte et al. [START_REF] Motte | Total solids 535 content: a key parameter of metabolic pathways in dry anaerobic digestion[END_REF], by gas and ion chromatography, respectively. Starting with the methane yields, the modeled results indicated that both linear coefficients 251 from the biochar and the FeCl 3 (i.e. a biochar and a FeCl3 ) had a significant effect on this variable 252 (p-values linear coefficients < 0.05). However, the quadratic parameter a biochar-biochar was 253 statistically different than zero with a 99 % of certainty (p-value < 0.01) and the high value of 254 this parameter indicated that the experimental design was not properly centered for predicting 255 the methane yield. The relatively low R 2 value (86.7 %) also suggests the lack of fit of the 256 model. In addition, the Durbin-Wilson test (with a p-value < 0.05) could not exclude that 257 correlations existed due to the order in which the data were used as input. All these results 258 suggest that the model was not able to predict any direct effect of the biochar or the FeCl 3 on 259 the final methane yields obtained, indicating that no significant differences existed between 260 the obtained results. Considering that the biodegradable matter content of these two additives 261 should be negligible, this is a logical outcome. 262 Nevertheless, the model was able to reproduce precisely the maximum methane production 263 rates (R 2 of 94.2 % and p-value Durbin-Wilson 0.20) and the average daily methane 264 production rates (R 2 of 100 %), indicating that highly significant differences existed between 265 the reactors. The model responses for these two variables are presented in Figure 1. 266

The maximum methane production rate was mainly affected by the biochar concentration, 267 with a negligible influence of the industrial FeCl 3 concentration. This can be verified by the p-268 values of the parameters in Table 5, where it can be observed that the only parameter with a p-269 value < 0.5 was a biochar . This indicates that addition of biochar clearly improved the 270 degradation of the HAc that accumulated at the beginning of the AD process, when the 271 maximum methane production rates were registered (Figure S1 and Figure S2). This is in 272 agreement with previous results, where GAC addition enhanced the HAc consumption during 273 batch FW AD (30). Regarding the average daily methane production rates (estimated by 274 dividing the total volume of methane produced by the time to completely degrade the HPr), the value of this variable also increased greatly at higher biochar concentrations. However, in 276 this case the effect of the FeCl 3 was significant, with both a biochar and a FeCl3 showing p-values 277 lower than 0.5. This variable was mainly affected by the batch duration, which was 278 determined by the time required for HPr consumption. Thus, it can be concluded that both 279 biochar and industrial FeCl 3 addition favored the degradation of the HPr accumulated during 280 the first days of the batch process. This can be related to favored syntrophic interactions and 281 to the occurrence of DIET in the case of biochar and to the synthesis of enzymes for the 282 industrial FeCl 3 (i.e., TEs contained in this solution). Again, this is in agreement with 283 previous results presented in the literature. Adding GAC and TEs into batch reactors, a 284 reduction of the time required to completely degraded the accumulated HPr of around 40 % 285 was observed during batch FW AD (30). 286

Although optimal values could be extrapolated (162 g•l -1 of biochar at 0.1 g Fe•l -1 for 287 maximizing the methane production rate and 111 g•l -1 of biochar at 0.2 g Fe•l -1 for 288 maximizing the average methane production rates), it is clear that the experimental design 289 was no properly centered and therefore further experiments must be performed to precisely 290 optimize the dosing of these reagents. 291

Performance of the continuous pilot reactors 292

Two continuous pilot reactors were run in parallel to the batch reactors during 196 days: a 293 control reactor fed only with FW and a reactor supplemented with biochar and industrial 294 FeCl 3. The main operational parameters and the obtained results are presented in Figure 2. methane yields, lower concentrations of HAc and HPr (i.e. peaks of HAc of 13 and 16 g•l -1 , 301 respectively) and higher pH values. This suggests that the added AD enhancers favored the 302 consumption of the VFAs accumulated during the start-up period, even at a biochar 303 concentration of 10 g•l -1 . This is in agreement with the batch results, confirming the positive 304 effect of adding these reagents on the initial HAc consumption (leading to a lower extent of 305 HPr accumulation) and on the HPr degradation. 306

Moving forwards, after 60 days both reactors achieved efficient methane productions at the 307 first OLR applied, with yields around 400 ml CH 4 •g VS -1 (95 % of the BMP), with relatively 308 low concentrations of VFAs in the reactors (2 g•l -1 of HAc and 0.3 g•l -1 for HPr) and with high 309 pH values (around 8.1). It must be mentioned that high TAN concentrations were already 310 present in the reactors at this point, with values around 8 g•l -1 in both reactors at day 68. It 311 must also be commented that the average TAN concentrations in the reactors were 8513 ± 362 312 mg•l -1 in control reactor and 8432 ± 777 mg•l -1 in the biochar-supplemented reactor 313 (corresponding to 398 and 550 mg FAN•l -1 , respectively). This indicates that no significant 314 differences existed and that the observed improvement of the AD performance was no caused 315 by a reduction in the TAN/FAN concentrations related to its adsorption on the biochar 316 particles. 317 However, due to the low OLR applied, relatively low methane production rates (around 10 l 318 CH 4 •d -1 ) were obtained. Therefore, the OLR was doubled in day 77 to reach 2.8 g VS•l -1 •d -1 . 319 This caused a sudden drop in the methane yields to 160 ml CH 4 •g VS -1 , which was associated 320 with an increase in the HAc and HPr concentrations in the reactors. In agreement with the 321 previous results, the levels of both VFAs were always lower in the supplemented reactor. 322

In an attempt to reduce the intensity of VFA accumulation in the reactor containing biochar 323 and FeCl 3 , the biochar concentrations were increased to 20 g•l -1 on day 105 and to 50 g•l -1 on 324 day 146 (based on the results from the batch optimization experiment described above). While the first increase did not have significant effects, the second one (to 50 g•l -1 ) caused a drop in 326 the HAc concentration from 13 to 7 g•l -1 , raising the methane yields up to 350 ml CH 4 •g VS -1 . 327 A consequent decrease in the HPr concentration was observed (from 3.1 to 1.8 g•l -1 ). Sadly, 328 on day 167 a problem occurred with the heating of the biochar-supplemented reactor. This 329 caused a temperature drop, which led to a sudden decrease of the methane yields obtained 330 and, again, an accumulation of HAc at the end of the operational period (days 167 to 196). 331

Besides the aforementioned complication, Figure 2 clearly shows that addition of biochar and 332 industrial FeCl 3 decreased the HPr concentrations in reactors, with considerable differences 333 between the supplemented reactor and the control. This discrepancy was particularly 334 important after increasing the biochar concentrations, with HPr levels of 7.2 g•l -1 in the 335 control reactor and of 1.8 g•l -1 in the reactor containing biochar and industrial FeCl 3 at the end 336 of the operational period. 337

Although the obtained results further suggest that addition of biochar and industrial FeCl 3 can 338 improve the AD kinetics and favor VFA consumption during FW AD, it is also clear that 339 further continuous experiments must be carried out, allowing longer operational periods. 340

These experiments should aim to reach an operational steady-state, with results that can be 341 extrapolated to a potential industrial-scale facility. 342

Biochar as feasible option for enhancing FW AD 343

The results obtained both in batch and continuous reactors are in agreement with different 344 studies carried out to study the effect of biochar on AD. The kinetics of consumption of HAc Using GAC as carbon-based AD enhancer, different studies have suggested that HAc, HBu 351 and HPr can be directly metabolized through DIET, improving the kinetics of consumption of 352 these VFAs [START_REF] Dang | Enhancing 453 anaerobic digestion of complex organic waste with carbon-based conductive materials[END_REF][START_REF] Zhao | Communities stimulated with 461 ethanol to perform direct interspecies electron transfer for syntrophic metabolism of 462 propionate and butyrate[END_REF][START_REF] Viggi | Magnetite particles 464 triggering a faster and more robust syntrophic pathway of methanogenic propionate 465 degradation[END_REF][START_REF] Zhao | Potential 467 enhancement of direct interspecies electron transfer for syntrophic metabolism of propionate 468 and butyrate with biochar in up-flow anaerobic sludge blanket reactors[END_REF]. In addition, biochar has also been found to promote the growth of 353 bacteria known to participate in DIET (i.e. Geobacter species) onto its surface [START_REF] Zhao | Potential 467 enhancement of direct interspecies electron transfer for syntrophic metabolism of propionate 468 and butyrate with biochar in up-flow anaerobic sludge blanket reactors[END_REF]. 354 Therefore, it can be hypothesized that the improved VFA degradation kinetics could be 355 related to an enhancement of the syntrophic interactions between microorganisms via biofilm 356 formation and to the occurrence of DIET. The degradation of HAc through DIET has already 357 been proposed in the literature and, although being more limited thermodynamically [START_REF] Capson-Tojo | 417 Accumulation of propionic acid during consecutive batch anaerobic digestion of commercial 418 food waste[END_REF], 358 DIET may have also played an important role in the oxidation of HPr. Besides, even if direct 359 DIET of HPr might not have occurred extensively, its degradation would be favored anyway 360 due to lower HAc and hydrogen/formate concentrations. Further studies analyzing the 361 microbial communities attached on the biochar, as well as the properties of the biochar used, 362 should be performed to verify this hypothesis. Concerning the industrial FeCl 3 addition, this 363 additive favored the HPr degradation due to the supplementation of TEs, critical for enzyme 364 synthesis [START_REF] Banks | Trace element requirements for stable food 423 waste digestion at elevated ammonia concentrations[END_REF]. 365

These experiments proved that a regular biochar (natural slow-pyrolyzed wood charcoal) 366 could also improve greatly the AD performance using complex FW as substrate. It must also 367 be considered that, other than the concentration applied, many parameters and variables which 368 have not been considered in this study have a huge potential for optimization when 369 considering biochar as AD enhancer. It is clear that the textural characteristics of the biochar 370 (e.g. specific surface, pore volume, pore size or pore distribution) as well as its surface 371 chemistry (e.g. hydrophobicity) or its particle size play a major role on biofilm formation. In 372 addition, also its resistivity (conductivity) might have a huge impact on its capability for 373 favoring DIET. All these characteristics are dependent on different variables that clearly 374 deserve further study, such as the raw material used for biochar production [START_REF] Shen | Yearlong semi-continuous 508 operation of thermophilic two-stage anaerobic digesters amended with biochar for enhanced 509 biomethane production[END_REF]

HCl (g•l -1 ) 2.2 Mn (mg•l -1 ) 780
Zn (mg•l -1 ) 390
Pb (mg•l -1 ) 220
Ni (mg•l -1 ) 67
Co (mg•l -1 ) 28
Cu (mg•l -1 ) 65
Cr (mg•l -1 ) 45
Ca (mg•l -1 ) 540
Na (mg•l -1 ) 110
Al (mg•l -1 ) 100

  -Tojo et al.

FeCl 3 .

 3 As it was demonstrated in Capson-Tojo et al. (30), although the separate addition of 147 GAC and TEs also enhanced the VFAs degradation kinetics during FW AD, the best 148 performances were achieved when both of these reagents were added simultaneously. With 149 these results as starting point, the current experimental designed was defined. The same study 150

  endogenous respiration of the inoculum. 177

  mass of VS of substrates. They were corrected to take into account the digestate 226 removed when sampling. The concentrations of FAN in the reactors were determined 227 according to Rajagopal et al.[START_REF] Rajagopal | A critical review on inhibition of anaerobic 415 digestion process by excess ammonia[END_REF], as a function of temperature, pH, ionic strength and TAN 228 concentration. The concentrations of the main ions present in the reactors (Cl -, Na + , NH 4 + , K + , 229 Mg 2+ , H + , VFAs and Ca 2+ ) were taken into account in this calculation. 230 231Results and discussion 232Optimization of biochar and FeCl 3 dosing in batch reactors 233As aforementioned, the first objective of this study was to evaluate the effect of different 234 concentrations of both biochar and industrial FeCl 3 , aiming at optimizing their dosage. The 235 results of the experimental design (using the second consecutive batch) are shown in Table4. 236All the reactors produced methane efficiently, with yields ranging from 456 to 505 ml CH 4 •g 237 VS -1 and final pH values between 8.15 and 8.42. These high pH values were caused by the 238 high TAN concentrations, up to 9.75 g•l -1 . As it can be observed, the main differences were 239 related to the AD kinetics. The maximum methane production rates (consequence of the 240 initial consumption of HAc) varied widely, from 907 to 1498 ml•g VS -1 •d -1 . The total time for 241 consumption of HPr (most recalcitrant VFA to be degraded) ranged from 18.9 to 24.4 days. 242The average daily methane production was calculated by dividing the total volume of methane 243 produced by the batch duration (i.e. time to completely degrade HPr), according to Capson-244 Tojo et al.(30). It varied from 302 to 376 ml CH 4 •d -1 . The values of the methane yields, the 245 maximum methane rates and the average daily methane production rates were used as inputs 246 for the quadratic model (Equation1), obtaining the results presented in Table5. The raw 247 kinetic curves of the methane yields and production rates, the pH values and the 248 concentrations of the different VFAs are presented in FigureS1and FigureS2(supplementary 249 material). 250

FigureFigure 1 .Figure 2 .

 12 Figure and Table captions 551

Figure S1 .Figure S2 .

 S1S2 Figure S1. Evolution of (A) the methane yields, (B) the methane production rates and (C) the pH

  the obtained results suggest for the first time that biochar and industrial FeCl 3 382 can be a feasible alternative for stabilizing AD of FW, favoring the consumption of VFAs and 383 improving the methane productivities. When compared to substrate dilution (1:1 vol:vol), 384 which is the most commonly applied stabilization method in industrial facilities (increasing 385 greatly the HRTs and the volumes of digestate produced and requiring an input of clean 386 water), the proposed technology is clearly a more environmental-friendly option that can be 387 coupled with other waste treatment processes (i.e. green waste pyrolysis). Niu Q, Chi Y, Li Y. Trace metals requirements for continuous thermophilic 451 methane fermentation of high-solid food waste. Chem Eng J 222: 330-336 (2013).

	o r P e e r R e v i e temperature and pressure applied during pyrolysis (i.e. slow or fast pyrolysis) or the grinding) (32). Although deep techno-economic analyses must be carried out before considering its application at industrial scale (particularly at the high concentrations of biochar applied, which would require recirculation of the solid fraction of the digestate for biochar reutilization), 388 378 379 380 381 F o r P e e r R e v i e and co-digestion proportions. Bioresour Technol 233: 166-75 (2017). F o r P e e r R e v i e 12. F o r P e e r R e v i e 21. Qiang H, F o r P e e r R e v i e F o r P e e r R e v i e F o r P e e r R e v i e F o r P e e r R e v i e , the F w w w w w w w w

Table 1 .

 1 Characteristics of the food waste mixture and the inoculum (9) 585

		Parameter	Food waste mixture	Inoculum
		TS (%)	21.0	5.81
		VS/TS (%)	90.3	59.1
		Carbohydrates (g•kg TS -1 )	618	n.m. 1
		Proteins (g•kg TS -1 )	187	n.m. 1
		Lipids (g•kg TS -1 )	121	n.m. 1
		BMPs (ml CH 4 •g VS -1 )	420	n.m. 1
		pH	5.02	8.10
		TAN (g•kg TS -1 )	0.90	5.04
		TKN (g•kg TS -1 )	30.0	93.0
		C/N	16.1	3.04
	586 587 588	1. n.m. stands for "not measured" F o r
	589 590	P e
	591 592		e r
	593 594		R e
	595		v
	596 597 598		i e w
	599
	600
	601
	602
	603
	604
	605
	606

Table 2 .

 2 Characteristics and composition of the industrial FeCl 3 solution 607

	Parameter	Unit	Value
	Density at 20 °C	(g•cm -3 )	1.45
	FeCl 3	(%)	41.1
	Cl total	(g•l -1 )	397
	Fe total	(g•l -1 )	206
	F o r
	P e
	e r
		R e
		v
		i e
		w

Table 3 .

 3 Experimental design of the batch experiment. All the reactors were fed with 60 g of

	625	FW at an S/X ratio of 1 g VS•g VS -1 and incubated at 37 °C
		Reactor	Normalized level Biochar FeCl 3	Real concentrations (g•l -1 ) Biochar FeCl 3
		1	0	-1	55	0.1
		2	-1	1	10	0.2
		3	1	1	100	0.2
		4	-1	-1	10	0.1
		5	1	-1	100	0.1
		6	0	1	55	0.2
		7	0	-1	55	0.1
		8	-1	-1	10	0.1
		9	0	1	55	0.2
	626	10 11 12	1 F -1 o 1 r	1 1 -1	100 10 100	0.2 0.2 0.1
	627		P e
	628 629			e r
	630			R e
	631			v
	632			i e
	633			w
	634
	635
	636
	637
	638
	639

Table 4 .

 4 Results of the batch experimental design 640

	Reactor	Methane yield (ml•g
		F o r
		P e
		e r
		R e
		v
		i e
		w

VS -1 ) Maximum methane rate (ml•g VS -1 •d -1 ) Time for HPr consumption (d) Average daily methane production rates (ml•d -1 ) 1

  Calculated as final methane yield divided by the time required for HPr consumption (batch duration)

						Final
						pH
	1	483	1327	18.9	363	8.15
	2	505	948	23.0	302	8.22
	3	484	1498	18.9	376	8.19
	4	509	886	24.4	280	8.21
	5	489	1436	18.9	374	8.18
	6	459	1249	18.9	356	8.29
	7	461	1281	20.1	329	8.29
	8	501	907	21.5	316	8.17
	9	456	1142	20.1	323	8.29
	10	466	1489	18.9	361	8.27
	11	496	913	21.5	326	8.42
	12	496	1457	20.1	360	8.19
	1. 641
	642
	643
	644
	645
	646
	647
	648
	649
	650
	651
	652
	653
	654
	655
	656

Table 5 .

 5 Coefficients of the quadratic model for the main responses of the experimental

	658	design
		Parameter/coefficient	Methane yield (ml•g
		F o r
			P e
			e r
			R e
			v
			i e
			w

VS -1 ) Maximum methane rate (ml•g VS -1 •d -1 ) Average daily methane production rates (ml•d -1 ) 1

  Calculated as final methane yield divided by the time required for HPr consumption (batch duration)

	a 0	465	1250	343
	a biochar	-9.54 *	278.2 **	37.3 ***
	a FeCl3	-9.04 *	-4.71	4.02 ***
	a biochar-biochar	28.3 **	-57.9	-12.2 ***
	a biochar-FeCl3	-8.03	3.25	-6.94 ***
	R 2	86.7 %	94.2 %	100 %
	p-value Durbin-Wilson	0.025	0.20	-
	1.

pH drop[START_REF] Capson-Tojo | Food waste 396 valorization via anaerobic processes: a review[END_REF][START_REF] Capson-Tojo | Kinetic study 402 of dry anaerobic co-digestion of food waste and cardboard for methane production[END_REF][START_REF] Stabnikova | The use of hybrid anaerobic 405 solid-liquid (HASL) system for the treatment of lipid-containing food waste[END_REF][START_REF] Kim | Effect of enzymatic pretreatment 408 on acid fermentation of food waste[END_REF][START_REF] Wang | Lactic acid fermentation of food waste using 410 integrated glucoamylase production[END_REF]. The second issue to be dealt with is related to the high protein content of this

substrate and its low water content (~20 % total solids; TS). During AD, organic nitrogen

(usually in the form of proteins) is reduced into ammonia nitrogen (total ammonia nitrogen;

TAN), which, in its free form (free ammonia nitrogen; FAN) is toxic to microorganisms,

especially to methanogenic archaea, eventually resulting also in an accumulation of VFAs.

http://mc.manuscriptcentral.com/jctb-wiley Journal of Chemical Technology & Biotechnology

After starting the reactors at an OLR of 1.4 g VS•l -1 •d -1 , it was required to wait around 30-40

and HBu were reported to be faster when adding biochar using glucose as substrate[START_REF] Luo | Application of eco-compatible biochar in anaerobic 516 digestion to relieve acid stress and promote the selective colonization of functional microbes[END_REF]. Also

pretreatment applied to the biochar before its addition into the AD reactor (i.e. mechanical

elements. Waste Manag 32 (8): 1509-1515 (2012).


J C Boulet

M Brandolini-Bunlon

G Chaix

M Lesnoff

A Mallet

S Mas

M Metz

J M Roger

V Rossard

R Servien

ChemHouse: A research and development centre for chemometrics

The aim of this article is to present the different activities of ChemHouse and to give an overview of the first years of operation. Detailed information can be found on the ChemHouse website (http://chemproject.org/chemhouse). A few years ago, some researchers from Montpellier, France, created ChemHouse, a multiinstitute research cluster (INRAE, CIRAD, Irstea, University of Montpellier). This creation was guided by the development of three tools, mainly applied to near infrared spectrometry (NIRS): CheMoocs, a MOOC dedicated to chemometrics for NIRS; ChemFlow, a free and open software tool, allowing anyone to implement the techniques learned in CheMoocs and ChemData, an open database. ChemHouse aims at ensuring an open and shared scientific animation to encourage national and international collaborations in chemometrics, in particular in the form of hosting researchers, and to allow the collaborative development of own research. ChemHouse also hosts the forges of the three tools: CheMoocs, ChemFlow and ChemData. Today, ChemHouse has 48 members (Cf https://chemproject.org/chemHouse/team). Every fortnight, ChemHouse members are invited to meet to discuss the operational and research issues of the cluster, without any restrictions. At each session, a member (or an outsider, if invited by a member) leads a scientific seminar around a presentation on a topic of their choice. More than 40 scientific presentations have been held in ChemHouse over the two years: 2019 and 2020. Some specific sessions are organised in the form of collective work on data and processing methodology, with for example participation in scientific conference shootouts. A list of these presentations and their content is available on the ChemHouse website at http://chemproject.org/chemhouse/ressources. Many ChemHouse seminars have been devoted to topical research issues:

• The calibration of NIRS models is certainly the research question that has attracted the most contributions. The use of NIRS in areas of massive data acquisition, such as precision agriculture, poses problems of model calibration. The adaptation of multivariate calibration to large databases is therefore a current topic, which has been the subject of seven ChemHouse seminars. A first solution consists of performing linear models (PLS) on a neighbourhood of the point to be predicted. This approach, known as Local-PLS, makes it possible to manage the non-linearities inherent in large databases. Three seminars were dedicated to local PLS: weighted PLS, ParSketch-PLS and RoBoost-PLS. [START_REF] Lesnoff | Comparison of locally weighted PLS strategies for regression and discrimination on agronomic NIR data[END_REF][START_REF] Metz | Un algorithme big-data pour la local-PLS[END_REF][START_REF] Metz | RoBoost-PLSR: a novel robust PLS regression method inspired from boosting principles[END_REF] Some seminars were held on other non-linear methods, such as artificial neural networks.

• Another important topic in ChemHouse is the preprocessing of spectra. Several contributions have been made on orthogonal projections, which can be used to avoid unwanted spectral variations. [START_REF] Herrero-Langreo | Orthogonal projection as a spectral pre-treatment method to reduce the interference of polystyrene signal in NIR imaging of agar on petri-dishes[END_REF][START_REF] Boulet | Further investigations on the relationship between the OPLS preprocessing and the NAS[END_REF][START_REF] Ryckewaert | Reduction of repeatability error for analysis of variance-Simultaneous Component Analysis (REP-ASCA): application to NIR spectroscopy on coffee sample[END_REF][START_REF] Mishra | Parallel preprocessing through orthogonalization (PORTO) and its application to near-infrared spectroscopy[END_REF][START_REF] Roger | Sequential preprocessing through ORThogonalization (SPORT) and its application to near-infrared spectroscopy[END_REF] However, completely new methods have also been developed. The Variable Sorting for Normalization (VSN) [START_REF] Rabatel | VSN: variable sorting for normalization[END_REF] offers an interesting alternative to the widely used but also much criticised Standard Normal Variate (SNV) method. The combination of several preprocessings via multiblock methods is also a new idea. Two new methods, Sequential preprocessing through ORThogonalization (SPORT) [START_REF] Roger | Sequential preprocessing through ORThogonalization (SPORT) and its application to near-infrared spectroscopy[END_REF] and Parallel pre-processing through orthogonalization (PORTO), [START_REF] Mishra | Parallel preprocessing through orthogonalization (PORTO) and its application to near-infrared spectroscopy[END_REF] have recently been developed in ChemHouse, in collaboration with external partners. • Another topic of interest for ChemHouse researchers is variable selection. Thus, Successive Orthogonalized Covariate Selection (SO-CovSel), [START_REF] Biancolillo | SO-CovSel: a novel method for variable selection in a multiblock framework[END_REF] a method for selecting variables in multi-block data sets, was developed. A seminar was also devoted to the selection of variables in functional data.

11

• In addition, many other seminars were devoted to open questions that did not lead to publication, such as the selection of the dimension of a PLS model, the estimation of prediction uncertainties, the estimation of the Net Analyse Signal (NAS), robust principal component analysis, the statistical properties of standard prediction errors, cross-validation, oblique projections or the boosting/bagging/stacking of models.

• Furthermore, ChemHouse regularly organises workshops on key research issues. These workshops, lasting from half a day to a full day, bring together ChemHouse researchers and external partners. Among others, the topics of local PLS, discrimination, hyperspectral imaging, soil NIRS, Python development, prediction uncertainties, software versioning and Git, use of the RNIRS package (https://github.com/mlesnoff/rnirs), etc. have been studied and discussed.

ChemHouse has produced papers and seminars dedicated to NIRS applications in the fields of phenotyping, wood quality, [START_REF] Chaix | Impact of drought on eucalyptus wood chemistry by nearinfrared hyperspectral imaging[END_REF] waste characterisation, [START_REF] Mallet | Unveiling nonlinear water effects in near-infrared spectroscopy: a study on organic wastes during drying using chemometrics[END_REF] process monitoring [START_REF] Rey-Bayle | Multiangle near-infrared spectroscopy associated with common components and specific weights analysis for in-line monitoring[END_REF] or calibration transfer.

15

ChemHouse is also a place for hosting researchers. In total, 57 weeks of hosting have been carried out since the cluster was created. These stays are of very diverse natures. Doctoral students are hosted for short periods to assist them in processing their data. More experienced researchers are hosted for longer periods to carry out methodological developments, such as SO-CovSel, 10 SPORT 8 or VSN. [START_REF] Rabatel | VSN: variable sorting for normalization[END_REF] International partnerships are established with the University of Rome "La Sapienza" (IT), the University of Aquilla (IT), the University of Barcelona (ES), the University of Wageningen (NL), the University of Dublin (IRL), the University of Bilbao (ES), the University of Tarragona (ES) and the Polytechnic University of Madrid (ES). In addition to these direct hosts, ChemHouse has a wider impact due to the privileged interactions between its members and their scientific partners. This is the case, for example, of CIRAD's ChemHouse members who work in partnership in tropical countries and transfer the chemometrics skills and know-how produced and exchanged in ChemHouse, in Africa (Madagascar, Côte d'Ivoire) and South America (Brazil, Argentina). ChemHouse also maintains a private partnership. The company Ondalys is actively involved in the scientific activities and in the organisation of events. The company Pellenc Selective Technologies participates as a donor. The CheMoocs MOOC (chemproject.org/chemoocs) was developed by a broad collective of French-speaking chemometricians, many of whom have joined ChemHouse to continue its development. Since 2016, this MOOC offers a complete teaching in chemometrics, requiring very few prerequisites. Each year, about 1500 people register for CheMoocs and about a hundred people successfully complete the course. ChemHouse is also the forge of ChemFlow [START_REF] Rossard | ChemFlow, chemometrics using Galaxy[END_REF] (chemproject.org/chemflow), a free and open source chemometrics software. Based on a Galaxy-project platform, this software allows any user, after registration, to access a large number of chemometrics tools such as spectral preprocessing, unsupervised analysis, multivariate regression, discriminant analysis, multiblock analysis, variable selection methods. ChemFlow can host code written in different languages (R, Python, Scilab, Octave), which allows the collection of scripts from the ChemHouse community and its partners. ChemFlow can generate processing workflows, making process automation and collaboration between users possible. ChemFlow is also used as an educational tool, even internationally (Cf https://chemproject.org/ressources/trainings). It is also the support of CheMoocs, for the realization of exercises and challenges. ChemFlow training courses are for an audience of 10 to 50 people. These training sessions aim to make ChemFlow users autonomous and to teach good practices. ChemHouse is part of an open science approach. With this in mind, ChemFlow and ChemData (https://chemproject.org/ChemData) have been developed. This project proposes to share the data of ChemHouse members in a three-sided scheme. The data are described in a data paper, [START_REF] Zgouz | Dataset of visible-near infrared handheld and micro-spectrometerscomparison of the prediction accuracy of sugarcane properties[END_REF] they are the purpose of CheMoocs exercises and they are made available to the community in a dataverse.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: ChemHouse is supported by the sponsorship of the company PELLENC Selective Technologies (www.pellencst.com), via a SupAgro Foundation project and grants from the INRAe MathNum department. (https://www.inrae.fr/departements/mathnum).

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.


J Jacques
email: julien.jacques@polytech-lille.fr

C Biernacki

Extension of model-based classification for binary data when training and test populations differ

Standard discriminant analysis supposes that both the training sample and the test sample are derived from the same population. When these samples arise from populations differing from their descriptive parameters, a generalization of discriminant analysis consists in adapting the classification rule related to the training population to another rule related to the test population, by estimating a link map between both populations. This paper extends an existing work in the multinormal context to the case of binary data. In order to solve the problem of defining a link map between the two binary populations, it is assumed that the binary data result from the discretization of latent Gaussian data. An estimation method and a robustness study are presented, and two applications in a biological context illustrate this work.

Introduction

Consider a first (training) sample of individuals described by explanatory variables, for which a partition in groups or classes is known. Consider also a new sample of individuals (test sample), drawn from the same population of the training sample. Discriminant analysis consists in estimating an allocation rule from the training sample in order to classify individuals of the test sample (see [START_REF] Mclachlan | Discriminant analysis and statistical pattern recognition[END_REF] for a survey). Since works of [START_REF] Fisher | he use of multiple measurements in taxonomic problems[END_REF], who introduced a linear discriminant rule between two groups, numerous evolutions have been proposed. All of them concern the nature of the discriminant rule: Parametric quadratic rule (see for example [START_REF] Tomassone | A generalized discriminant for sexing fulmarine petrels from external measurements[END_REF]), semi parametric rule (as logistic discrimination, [START_REF] Anderson | Separate sample logistic discrimination[END_REF]) or non parametric rule [START_REF] Fix | Discriminatory analysis -non parametric discrimination: Consistency properties[END_REF], [START_REF] Friedman | Projection pursuit regression[END_REF], [START_REF] Hand | of Electronic & Electrical Engineering Research Studies: Pattern Recognition & Image Processing Series[END_REF], [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]). An alternative approach, introduced by Van Franeker and Ter Brack [1993] and developed further by [START_REF] Biernacki | A generalized discriminant rule when training population and test population differ on their descriptive parameters[END_REF], considers the case in which the training sample does not necessarily arise from the same population as the one of the test sample. Biernacki et al. define several models of generalized discriminant analysis in a multinormal context, and conduct experiments for biological data consisting of birds from the same species, but with different geographical origins. In many domains (insurance, medicine, biology, etc.) a large number of applications deals with binary data as well. The goal of the present paper is to extend the generalized discriminant analysis of [START_REF] Biernacki | A generalized discriminant rule when training population and test population differ on their descriptive parameters[END_REF] to the case of binary data. The paper is organized as follows. The next section presents the data and the latent class model for both training and test populations. Section 3 makes the assumption that these binary data are discretized latent continuous variables in which the order information is lost. This hypothesis is the key to establish a general stochastic link map between the two populations, from which many pertinent parsimonious sub-models can be obtained. In Section 4, estimation is performed based on the maximum likelihood principle using the EM algorithm. Then, a robustness study for the Gaussian assumption is carried out in Section 5 involving both theoretical and experimental (simulated data) arguments. In Section 6, two applications in a biological context illustrate realistic situations where the proposed generalized discriminant analysis outperforms standard discriminant analysis and clustering. Finally, the last section concludes this paper by discussing possible extensions of the present work.

The data and the latent class model

The data consist in two samples: The first sample S, labelled and drawn from the training population P , and the second sample S, unlabelled and drawn from the test population P . A fundamental assumption of the present work is that populations P and P may be different.

The training sample S is composed of n pairs (x 1 , z 1 ), . . . , (x n , z n ), where x i is the binary explanatory vector for the ith object with x i = (x i1 , . . . , x id ) T ∈ {0, 1} d and where z i = (z i1 , . . . , z iK ) T is the group membership with z ik being equal to 1 if the ith object belongs to the kth group and being equal to 0 otherwise (i = 1, . . . , n, k = 1, . . . , K). The number of binary explanatory variables and the number of groups are respectively denoted by d and K. Each pair (x i , z i ) is assumed to be an independent realization of the random vector (X 1 , Z 1 ) with distribution: X 1j |Z 1k =1 ∼ B(α kj ) for all j = 1, . . . , d and Z 1 ∼ M(1, p 1 , . . . , p K ),

where B(α kj ) is the Bernoulli distribution of parameter α kj (0 < α kj < 1), and M(1, p 1 , . . . , p K ) defines the one order multinomial distribution of parameters p 1 , . . . , p K (0 < p k < 1, K k=1 p k = 1). Moreover, using the latent class model assumption of conditional independence of the explanatory variables [START_REF] Celeux | Clustering criteria for discrete data and latent class models[END_REF], [START_REF] With Discussion | An introduction to latent variable models[END_REF]), the probability function of X 1 conditionally to the group membership is:

f k (x 11 , . . . , x 1d ) = d j=1 α kj x 1j (1 -α kj ) 1-x 1j .
(2) This is one of the most popular generative method for discriminating categorical data since it is straightforward to implement and it is often efficient (see for instance experimental comparisons with alternative methods, not necessarily generative, in ?, ? Chap. 9 or also ?).

Similarly, the test sample S is composed of ñ pairs (x 1 , z1 ), . . . , (x ñ, zñ ), where the d explanatory variables are the same as in the training sample, but where the zi are unknown. These pairs are assumed to be independent realizations of ( X1 , Z1 ) with distribution:

X1j | Z1k =1 ∼ B( αkj )
for all j = 1, . . . , d and Z1 ∼ M(1, p1 , . . . , pK ).

(3)

The explanatory variables X1j , for j = 1, . . . , d, are also assumed to be conditionally independent.

Basically, the distribution of ( X1 , Z1 ) differs from this one of (X 1 , Z 1 ) only by the values of the parameters α kj and p k .

Our goal is to estimate the unknown labels z1 , . . . , zñ by using information from both training and test samples. The challenge resides in finding a link map between the populations P and P . Remarks

• In fact, since both labelled and unlabelled data are used together in the inference process, our problem is related to the so-called semi-supervised learning purpose. Obviously, the originality of our work is that the data sets do not necessarily arise from the same population.

• The use of the terminology "test" for the sample S (and the population P ) is abusive because this sample is used to determine the discrimination rule. Nevertheless, this terminology is adopted in order to facilitate the link with the standard discrimination methods. Moreover, it appears to be a usual notation in the semi-supervised classification community (see for instance ?, Chap. 1).

Relationship between test and training populations

Formalizing the link between populations

In a multinormal context, a linear stochastic relationship between P and P is not only justified (under very few assumptions that will be recalled later) but also intuitive [START_REF] Biernacki | A generalized discriminant rule when training population and test population differ on their descriptive parameters[END_REF]). In the binary context, since such an intuitive relationship seems more difficult to exhibit, an additional assumption is stated: The binary variables are supposed to result from the discretization of some latent Gaussian variables. For instance, if a binary variable is a product purchased by a customer, it is assumed that the customer gives a score to the product, and buys it only if this score is greater than a given threshold. This assumption is not new in statistics: See for example [START_REF] Thurstone | A law of comparative judgement[END_REF], who used this idea in his comparative judgment model to choose between two stimuli. ? also modelled multivariate ordered categorical variables as a latent multinormal distribution involving a possibly full correlation matrix (see ? for a more recent reference). Moreover, [START_REF] Everitt | A finite mixture model for the clustering of mixed-mode data[END_REF] proposed a classification algorithm for binary, categorical and continuous data. Thus, the explanatory variables X 1j |Z 1k =1 having a Bernoulli distribution B(α kj ) are assumed to arise from the discretization of latent continuous variables Y 1j |Z 1k =1 in the following manner:

X 1j |Z 1k =1 = 0 if λ j Y 1j |Z 1k =1 < λ j s j 1 if λ j Y 1j |Z 1k =1 ≥ λ j s j for j = 1, . . . , d, (4)
where s j ∈ R is the discretization threshold, and λ j ∈ {-1, 1} is introduced to avoid choosing to which value of X ) T is assumed to be multivariate normal, with marginal distributions N (µ kj , σ 2 kj ), such the obtained discretized variables X 1j |Z 1k =1 (j = 1, . . . , d) are independent in order to retrieve the conditional independence assumption on the binary variables. It should be noted that conditional independence of the latent Gaussian variables Y 1j |Z 1k =1 (j = 1, . . . , d) is a sufficient assumption for conditional independence of the binary variables. However, this assumption may not be necessary. From (1) and (4) the following relationship between α kj , λ j , µ kj and σ kj can be derived:

α kj = p(X 1j |Z 1k =1 = 1) = 1 -Φ λ j s j -µ kj σ kj ( 5
)
where Φ is the N (0, 1) cumulative density function.

As for the variable X 1j , the binary variable X1j is also assumed to arise from the discretization of a latent Gaussian variable Ỹ1j with distribution N (μ kj , σ2 kj ). The equations are the same as (4) and ( 5), by changing α kj into αkj , µ kj into μkj and σ kj into σkj . The thresholds sj are naturally supposed to be equal to s j (s j = s j ) since, with the previous example, it is equivalent to assume that customers of both populations buy the product if the score is greater than the same threshold. In the same spirit λj is supposed to be equal to λ j , so the rule of purchase -lower or higher than the threshold -is the same for both populations. In a Gaussian setting, [START_REF] Biernacki | A generalized discriminant rule when training population and test population differ on their descriptive parameters[END_REF] showed that the only possible link map between the latent continuous variable Y 1 |Z 1k =1 of P and Ỹ1 | Z1k =1 of P is linear when the two following reasonable hypotheses are satisfied: (i) The transformation between P and P is C 1 and (ii) the jth component

Ỹ1j | Z1k =1 of Ỹ1 | Z1k =1 is only related to the jth component Y 1j |Z 1k =1 of Y 1 |Z 1k =1
. More precisely, this relationship is expressed by

Ỹ1 | Z1k =1 ∼ A k Y 1 |Z 1k =1 + b k , (6)
where A k is a diagonal matrix of R d×d containing the elements a kj and b k is a vector of R d containing the elements b kj (1

≤ k ≤ K, 1 ≤ j ≤ d).
By using ( 6) and ( 5) the following relationship between the parameters αkj and α kj can be obtained (see details in Appendix A):

αkj = Φ δ kj Φ -1 (α kj ) + λ j γ kj , (7)
where δ kj ∈ R + \ {0}, λ j ∈ {-1, 1} and γ kj ∈ R. Note that this relationship corresponds to a linear link between the probit functions of both α kj and αkj . Conditionally to the fact that α kj are known (they will be estimated in practice), estimation of the Kd continuous parameters αkj is thus obtained from the estimated parameters of the link between P and P (plug-in method): δ kj , γ kj and λ j . Note that the choice of the discretization thresholds s j is not important. However, estimating the number of parameters for the link map is 2Kd and one thus obtain that the model is overparameterized. This fact should not be surprising since the underlying Gaussian model is by far more complex (in terms of the number of parameters) than the Bernoulli model. Hence we need to reduce the number of free continuous parameters in (7), and one way to do this is to propose sub-models defined via imposing natural additional constraints on the transformation between both populations P and P .

Models of constraints on the stochastic link map

The parameters δ kj (1 ≤ k ≤ K and 1 ≤ j ≤ d) will be successively constrained to be equal to 1 (denoted by 1), to be class-and dimension-independent (δ), to be only class-dependent (δ k ) or only dimension-dependent (δ j ). In the same way, γ kj can be constrained to be equal to 0, γ (constant w.r.t. k and j), γ k (constant w.r.t. j) or γ j (constant w.r.t. k). Thus, 16 models can be defined and indexed using the following ad-hoc notation: [1 0] means δ kj = 1 and γ kj = 0 (it corresponds to the usual discriminant analysis model),

[δ k γ j ] means δ kj = δ k and γ kj = γ j (1 ≤ k ≤ K and 1 ≤ j ≤ d), etc.
For these 16 models, an additional assumption on the group proportions is taken into account: Either the proportions of P are constrained to be equal to those of P , or they have to be estimated. In the following, [p k 1 0] denotes the model [1 0] with equal proportions whereas [p k 1 0] denotes this model with free proportions. The number of constrained models is thus growing to 32. Table 1 gives the number of (continuous) parameters to be estimated for these models. If the mixing proportions are different from P to P , K -1 must be added to these values.

Table 1: Number of continuous parameters (param.) to be estimated for the constrained models. model

[p k 1 0] [p k 1 γ] [p k 1 γ k ] [p k 1 γ j ] [p k δ 0] [p k δ γ] [p k δ γ k ] [p k δ γ j ] param. 0 1 K d 1 2 K + 1 d + 1 model [p k δ k 0] [p k δ k γ] [p k δ k γ k ] [p k δ k γ j ] [p k δ j 0] [p k δ j γ] [p k δ j γ k ] [p k δ j γ j ] param. K K + 1 2K K + d d d + 1 K + d 2d
Finally, we chose to use the BIC criterion (Bayesian Information Criterion, [START_REF] Schwarz | Estimating the dimension of a model[END_REF]) for automatic selection among the 32 generalized discriminant models. However, other criteria such AIC (An Information Criterion, ?) could be used as well. BIC is defined by: BIC = -2l( θ) + ν log(ñ), where θ = (p k , δ kj , λ j , γ kj ) for 1 ≤ k ≤ K and 1 ≤ j ≤ d, l( θ) is the maximum log-likelihood corresponding to the estimation θ of θ, and ν is the number of free continuous parameters associated to the given model. The model leading to the smallest BIC value is then selected. Before estimating the parameter θ by the maximum likelihood method, we need to discuss identifiability of each parametrization.

Model identifiability

Some of the constrained models previously defined can be non-identifiable. It is necessary to clarify these identifiability problems, which arise at two levels: Identifiability of the model parameters when α kj is transformed into αkj , and identifiability of the transformation which ensures that α kj can be only transformed into αkj and not into αk ′ j (with k ′ = k). We call respectively them intragroup and inter-group identifiability. The reader can find theoretical and experimental discussion about these two kinds of identifiability in Appendix B. The conclusion of this discussion is identifiability will occur in practical situations.

Parameter estimation

In this section, only the situation where proportions are unknown is presented, otherwise it is straightforward.

The three estimation steps

Generalized discriminant analysis needs three estimation steps. The first step consists in estimating parameters p k and α kj (1 ≤ k ≤ K and 1 ≤ j ≤ d) from population P based on training sample S. Since S is a labelled sample, the maximum likelihood estimate is simply given by [START_REF] With Discussion | An introduction to latent variable models[END_REF], [START_REF] Celeux | Clustering criteria for discrete data and latent class models[END_REF]):

pk = 1 n n i=1 z ik and αkj = 1 n n i=1
x ij z ik .

The second step consists in estimating parameters pk and αkj (1 ≤ k ≤ K and 1 ≤ j ≤ d) of population P by using pk , αkj (1 ≤ k ≤ K and 1 ≤ j ≤ d) and θ. Thus, for estimating αkj , the parameters δ kj , γ kj and λ j of the link between P and P have to be estimated, and then an estimate of αkj is deduced by plug-in inside Equation ( 7). This step is described below. Finally, the third step consists in estimating group membership of individuals from the test sample S, by maximum a posteriori.

Estimation of the link parameters

For the second step above, maximum likelihood estimation can be efficiently based on the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]). The likelihood is given by:

L(θ) = ñ i=1 K k=1 pk d j=1 αx ij kj (1 -αkj ) 1-x ij ,
and the completed log-likelihood is:

l c (θ; z1 , . . . , zñ ) = ñ i=1 K k=1 zik log pk d j=1 αx ij kj (1 -αkj ) (1-x ij )
.

The E step From a current value θ (q) of the parameter θ, the E step of the EM algorithm consists in computing the conditional expectation of the completed log-likelihood:

Q(θ; θ (q) ) = E θ (q) [l c (θ; Z1 , . . . , Zñ )|x 1 , . . . , xñ ] = ñ i=1 K k=1 t (q) ik log(p k ) + d j=1 log αkj xij (1 -αkj ) 1-x ij where t (q) ik = p( Zik = 1|x 1 , . . . , xñ ; θ (q) ) = p(q) k d j=1 ( α(q) kj ) xij (1 - α(q) kj ) (1-x ij ) K κ=1 p(q) κ d j=1 ( α(q) κj ) xij (1 - α(q) κj ) (1-x ij )
is the conditional probability for the individual i to belong to the group k.

The M step The M step of the EM algorithm consists in choosing the value θ (q+1) which maximizes the conditional expectation Q computed at the E step:

θ (q+1) = argmax θ∈Θ Q(θ; θ (q) ) (8)
where Θ is a parameter space depending on the model at hand. The M step is now described for each component of θ = {p k , δ kj , λ j , γ kj }.

For proportions, the maximum is:

p(q+1) k = 1 ñ ñ i=1 t (q)
ik .

The parameters δ kj and γ kj are never considered because the full generalized discriminant analysis model is overparameterized. Thus, only the constrained models are to be estimated. In this context, it is proved in Appendix C that Q is a strictly concave function of δ k , δ j , δ, γ k , γ j and γ. Thus, the maximisation of Q(θ; θ (q) ) is computed by an alternated iterative algorithm which consists in a succession, componentwise, of simplex algorithms if the optimization is unconstrained (γ k , γ j ∈ R). If the optimization is constrained (δ k , δ j > 0), the same algorithm is used but if the optimization leads to a negative result, the value 0 is retained (in this case the EM algorithm is used on its generalized form: GEM, [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]). The starting point of the alternating algorithm is θ (q) , and this one of the EM algorithm, θ (0) , is the point which corresponds to P = P . For the discrete parameters λ j , if the dimension d is sufficiently low, the maximization is carried out by computing Q(θ; θ (q) ) for all 2 d possible values for these discrete parameters. Since computational limits are quickly reached, a relaxation method can be used, which consists in assuming that λ j is not a binary parameter in {-1, 1} but a continuous one in [-1, 1], named λ * j (see [START_REF] Wolsey | Integer programming[END_REF] for instance). Optimization is thus performed on this continuous parameter, with the previous alternated algorithm since Q is a strictly concave function of λ j (Appendix C), and the solution λ * (q+1) j is then discretized to obtain a binary solution λ (q+1) j as follows: λ

(q+1) j = sgn(λ * (q+1) j
), where sgn denotes the sign function. This relaxed approach is not used in the experiments of this paper but see [START_REF] Jacques | Contributions à l'analyse de sensibilité et à l'analyse discriminante generalisée[END_REF] for some examples of use. Remark: Here, the estimation of (p k , α kj ) and θ is performed in a sequential fashion. This procedure enjoys the advantage to be algorithmically straightforward unlike the procedure involving the full likelihood function of all parameters. Moreover, experiments below indicate good behaviour of this strategy. However, since full likelihood estimates are expected to have less bias, this new way could be explored in future works (see Section 7).

Robustness study to the Gaussian assumption

In this section, we first prove that the Gaussian hypothesis can be weakened into a new assumption and, then, some experiments illustrate and evaluate the robustness of the methodology when this weakened assumption is itself violated.

Theoretical robustness study

Under the hypothesis that the link between the latent variables of both populations P and P is linear, the assumption of normality of the latent variables Y 1j |Z 1k =1 and Ỹ1j | Z1k =1 is in fact not necessary, and sufficient conditions are the following:

F Y 1j |Z 1k =1 (y) = Φ( y -µ kj σ kj ) for y = s j and y = s j -b kj a kj , (9)
where F Y 1j |Z 1k =1 denotes the conditional cumulative density function of variable Y 1j |Z 1k =1 and µ kj , σ kj its mean and standard deviation respectively. Note the difference between initial assumption in Subsection 3.1: This equality was stated for all y values in R, whereas it is now stated for only two specific values of y.

However it is straightforward to verify that Equations ( 5) and ( 11) remain valid in this context. Since they are the key for establishing Relationship (7), this latter is still right.

Practical robustness study

To illustrate the robustness against Condition (9) itself described in the previous subsection, consider the following example: For all 1 ≤ j ≤ 5 and

1 ≤ k ≤ 2, X 1j |Z 1k =1 (respectively X1j | Z1k =1
) is the binary discretization (with the threshold ) are given in (Figure 1).

s j = s) of Y 1j |Z 1k =1 (resp. Ỹ1j | Z1k =1 ),
The values of the discretization threshold, of the transformation parameters between P and P (a kj and b kj ) and of the first and second moments of Y 1j |Z 1k =1 (µ kj , σ 2 kj ) are arbitrary chosen: s = 13; a kj = 1.2 and b kj = 1 for 1 ≤ k ≤ 2 and 1 ≤ j ≤ 5; (µ 1j , µ 2j ) = (10, 13) and (σ 1j , σ 2j ) = (1.8, 1.7) for 1 ≤ j ≤ 5. Now, let choose (by numerical optimization) the parameters of the mixture density of Y 1j |Z k =1 with respect to the following constraints, for 0 ≤ ǫ ≤ 1 :

F Y 1j |Z 1k =1 (s j ) = Φ( s j -µ kj σ kj ) + (1 -Φ( s j -µ kj σ kj )) × ǫ F Y 1j |Z 1k =1 ( s j -b kj a kj ) = Φ( s j -a kj µ kj -b kj a kj σ kj ) + (1 -Φ( s j -a kj µ kj -b kj a kj σ kj )) × ǫ. (10)
When ǫ = 0, these constraints correspond to (9) and then satisfy the assumptions of Generalized Discriminant Analysis (GDA), and the greater is ǫ, the less (9) is respected. (Figure 2) illustrates the cumulative density function of the latent variables for different values of ǫ. The sample size is set to 50, the experiments are repeated 20 times and the mean error rate, estimated on "out-of-sample" data, is presented on (Figure 3) for different strategies: GDA with model [p k δ k γ k ], standard discriminant analysis and clustering. Moreover, the optimal error rate is also given. GDA outperforms usual methods when ǫ is moderate (ǫ < 0.22), and for higher ǫ, the assumed model of GDA is too incorrect and then clustering becomes better.

Comparison of methods on biological data 6.1 Discretized continuous data

The first motivations for which GDA was developed are biological applications [START_REF] Biernacki | A generalized discriminant rule when training population and test population differ on their descriptive parameters[END_REF], [START_REF] Tomassone | A generalized discriminant for sexing fulmarine petrels from external measurements[END_REF]), in which the aim was to predict sex of birds from biometrical variables. Very powerful results have been obtained with multinormal assumptions. The species of birds considered in the present application is Cory's Shearwater Calanectris diomedea [START_REF] Thibault | Cory's shearwater calonectris diomedia[END_REF]). Two subspecies can be identified: borealis which lives in the Atlantic islands (the Azores, Canaries, etc.) and diomedea which lives in the Mediterranean islands (Balearics, Corsica, etc.). A sample of Borealis (n = 206, 45% females) was measured using skins in several National Museums. Five morphological variables were measured: Culmen (bill length), tarsus, wings and tail lengths, and culmen depth. Similarly, a sample of diomedea (n = 38, 58% females) was measured using the same set of variables. In this example, two groups are present, males and females, and all the birds are of known sex (from dissection). (Figure 4) illustrates differences between the two subspecies borealis and diomedea, for two biometrical variables.

To provide an application of the present work, the continuous biometrical variables are discretized into binary data. As it can be shown on (Figure 4), discretization must be carried out carefully, especially concerning the choice of the discretization threshold. Indeed, if this threshold equals the mean of the biometrical variables for one subspecies for instance, then all the values for the other subspecies will be on the same side of this discretization threshold. Consequently, conditionally to each biometrical variables, the threshold is chosen such that there is (roughly) a maximum of individuals of each population on each side of the discretization threshold. The subspecies borealis is selected as the training population and the subspecies diomedea as the test population. The 32 GDA models, including standard discriminant analysis [p k 1 0], and clustering are applied on these data. The classification error rate and the value of the BIC criterion are given in Table 2.

[p k 1 0] [p k 1 γ] [p k 1 γ k ] [p k 1 γ j ] [p k δ 0] [p k δ γ] [p k δ γ k ] [p k δ γ j ] error
clustering. Indeed, by adapting the classification rule derived from the training population to the test population, GDA gives lower classification error rates than by applying directly the rule derived from the training population (standard discriminant analysis), or by omitting the training population and applying directly clustering on the test population.

It is worth pointing out that the assumption for binary data to be derived from the discretization of Gaussian variables (biometrical variables) is relatively realistic in this application. Nevertheless, there exists a strong correlation between the five biometrical variables, which violates the assumption that discretized variables are independent. Remark Although the "test sample" is used for estimating the discriminant rule through the unlabelled data of S and is used also for estimating the error rate but this once through the labelled data of S, this estimated error rate is usually not an optimistic measure of the classification method's performance. This fact is well-known in the semi-supervised setting (?, Subsection 21.1.2). In addition, we have also verified this claim with our biological data by performing a hold-out procedure (results not reported here).

Real binary data

The second application considers sea birds from the species puffins (?). Two groups of subspecies are considered: The first one is composed of subspecies living in Pacific Islands -subalaris (Galapagos Island), polynesial, dichrous (Enderbury and Palau Islands) and gunax -and the second one is composed of subspecies living in Atlantic Islands -boydi (Cap Verde Islands). Here, the difference between populations is the geographical range (Pacific vs. Atlantic Islands). A sample of Pacific birds (n = 171) was measured using skins in several National Museums. Four variables are measured on these birds: Coller, stripe and piping (absence or presence for these three variables) and under-caudal (self couloured or not). Similarly, a sample of Atlantic birds (n = 19) was measured using the same set of variables. Like in the previous example, two groups are present (males and females) and the sex of all the birds is known. Pacific birds are chosen as the training population and Atlantic ones as the test population. Choosing Atlantic birds as the test population corresponds to a realistic situation because it could be hazardous to perform a clustering process on a sample of such a small size. This is a typical situation where our methodology could be expected to provide a parsimonious and meaningful alternative. According to the biologist Vincent Bretagnolle, the morphological variables which are used in this application are not very discriminative, and then one can not expect that the error rate will be better than 40 -45%.

The 32 GDA models, among which standard discriminant analysis [p k 1 0], are applied on these data and the results are presented in Table 3. Clustering is also applied, and the obtained error rate is 49.05%.

As in the previous study, GDA is more efficient than standard discriminant analysis (50.94%) and clustering (49.05%) to classify birds according to their sex. Moreover the BIC criterion leads to choose the model with the smallest error rate. The relatively poor classification results (the minimal error rate is 43%) confirm the assumption of the biologist.

Conclusion

Generalized discriminant analysis extends standard discriminant analysis by allowing training and test samples to arise from different but stochastically linked populations. Our contribution consists in extending previous original work, derived in a multinormal context, to the case of binary data. Applications to a biological problem illustrate the power of our methodology. A classification of birds according to their sex is provided by using generalized discriminant analysis, and this classification is better than those obtained by standard discriminant analysis or by clustering. Perspectives for this work are numerous. Firstly, the parameters of both populations P and P are estimated successively: α kj and p k are estimated in a first step and then αkj and pk are deduced from these estimations and from those of all parameters. It should be very (computationally) useful to consider a joint estimation of this four parameters. In particular, ? extended the earlier work of [START_REF] Biernacki | A generalized discriminant rule when training population and test population differ on their descriptive parameters[END_REF] to the full likelihood estimation of parameters in the multinormal situation. It appears that error rates obtained by sequential estimate and by joint estimate are quite similar when the learning sample is large (situation of the current paper). But, when the learning sample has a small sample size, joint estimation could significantly improve the error rate. We could expect the same behaviour for our current model but this assumption needs to be confirmed by a future specific study. Secondly, the link between both populations was defined by using Gaussian cumulative density function. Although it seemed initially difficult to find this link, a simple link involving the probit function was obtained. It was not easy to imagine it, but it is meaningful afterwards. It would be interesting to try other types of cumulative density functions; Obviously theoretical reasons will have to be developed and practical tests will have to be carried out. Thirdly, with this contribution generalized discriminant analysis is now developed for continuous data and for binary data. To allow to analyse a large number of practical cases, it is important to study the case of categorical variables (i.e. more than two modalities), and thereafter the case of mixed variables (binary, categorical and continuous together). Everitt's works [START_REF] Everitt | A finite mixture model for the clustering of mixed-mode data[END_REF]), which defined a classification algorithm for mixed variables, can be helpful for this topics.

Finally, it would be also interesting to extend other classical discriminant method like non-parametric discrimination or semi-parametric discrimination. See [START_REF] Biernacki | Apprentissage sur une sous population et prédiction sur une autre : une extension à (et de) la discrimination logistique[END_REF] for logistic regression.

A Proof of the relation between test and training populations

From (6) it follows that Ỹ1j | Z1k =1 is Gaussian with mean μkj = a kj µ kj + b kj and with standard deviation σkj = |a kj |σ kj . However, this transformation is clearly non-identifiable: There exists more than one couple (A k , b k ) which satisfies Relationship (6). To solve this problem, all the a kj are assumed to be non negative (a kj ≥ 0). It is possible to derived from Equation ( 5):

λ j s j σ kj = -Φ -1 (α kj ) + λ j µ kj σ kj ,
and then, using expressions of μkj and σkj just given above:

αkj = 1 -Φ λ j s j -μkj σkj = 1 -Φ 1 a kj λ j s j σ kj -λ j a kj µ kj + b kj a kj σ kj = Φ Φ -1 (α kj ) a kj + λ j a kj µ kj + b kj -µ kj a kj σ kj . (11)
Alternatively, it is also possible to write from (5) again:

λ j µ kj σ kj = Φ -1 (α kj ) + λ j s j σ kj ,
and thus, following the same process as in ( 11),

αkj = 1 -Φ λ j s j -μkj σkj = 1 -Φ -λ j µ kj σ kj + λ j s j -b kj a kj σ kj = Φ Φ -1 (α kj ) + λ j a kj s j + b kj -s j a kj σ kj . (12)
Of course, Equations ( 11) and ( 12) are equivalent but the first one will be retained since Parametrization (11) will be more convenient for proposing later numerous parsimonious models of constraints on the link between P and P . Consequently, the following relationship between α kj and αkj is obtained:

αkj = Φ δ kj Φ -1 (α kj ) + λ j γ kj ,
where δ kj ∈ R + \ {0}, λ j ∈ {-1, 1} and γ kj ∈ R.

B Model identifiability B.1 Intra-group identifiability

Firstly, the identifiability of the couple (λ j , γ kj ) in different constrained model situations is reaching as follows:

• Models involving γ kj = 0 or γ kj = γ j : By setting for instance λ j = +1 (j = 1, . . . , d),

• Models involving γ kj = γ or γ kj = γ k : By setting for instance λ 1 = +1.

By this way, the product λ j γ kj is allways identifiable and obviously all these constraints on λ j have no impact on the estimation of the product λ j γ kj , as the reader can easily convince himself.

Secondly, the identifiability of the couple of parameters (δ kj , γ kj ) conditionally to λ j is discussed. Equation ( 7) leads to

Φ -1 ( αkj ) = δ kj Φ -1 (α kj ) + λ j γ kj ( 13
)
which can be expressed as the following linear system

Φ = Φ λ × u δ,γ
where Φ = (Φ -1 ( α11 ), . . . , Φ -1 ( αkj ), . . . , Φ -1 ( αKd )) T ∈ R Kd and Φ λ and u δ,γ are respectively a matrix and a vector, with dimension depending on the model at hand, representing the values of Φ -1 (α kj ) and λ j for Φ λ , and the values of δ kj and γ kj for u δ,γ .

Identifiability is obtained if and only if the matrix Φ λ is of full rank. It is easily noticed (see the example below) that Φ λ is not of full rank only for very particular values of α kj (typically α kj = 1 2 for some k, j for instance). Moreover these theoretical non-identifiable situations lead seldom to practical non-identifiable situations since an estimator αkj is used instead of the unknown true value α kj . To be definitively convinced of this fact, the reader can take a look at the robustness study (Section 5) where simulations involve a theoretically non-identifiable model (all α 1j = 1 2 ) but where practical identifiability is observed.

Example with model

Φ λ =                      Φ -1 (α 11 ) 0 . . . 0 λ 1 0 . . . 0 . . . . . . Φ -1 (α 1d ) 0 . . . 0
. . . 0 Φ -1 (α K1 ) 0 . . . 0 λ 1 . . . . . . 0 . . . 0 Φ -1 (α Kd ) 0 . . . 0 λ d                      .
Φ λ is not of full rank (min(Kd, 2K)) if and only if

• Possibility 1: There exists k ∈ {1, . . . , K} such that α kj = 1 2 for all j ∈ {1, . . . , d},

• Possibility 2: There exists k 1 , . . . , k τ ∈ {1, . . . , K}, k 1 , . . . , k τ being all different, with τ ∈ {1, . . . , K} satisfying d(Kτ ) < 2K, such that for all j, j ′ ∈ {1, . . . , d}, for all k ∈ {k 1 , . . . , k τ }, λ j ′ Φ -1 (α kj ) = λ j Φ -1 (α kj ′ ).

B.2 Inter-group identifiability

This non-identifiability problem means that one group of the population P can be transformed into more than one group of the population P . It cannot happen if the group proportions p k in P are all different (it is likely in practice) and if simultaneously all the constrained models are with fixed proportions between populations (models [p k . . ]). But if these proportions conditions are not verified, the non-identifiability problem is the following: For fixed k 1 and k 2 (k 1 , k 2 ∈ {1, ..., K}, k 1 = k 2 ) and for fixed j (j ∈ {1, . . . , d}) there exists two sets of parameters (δ k 1 j , λ j , γ k 1 j ) = (δ ′ k 1 j , λ ′ j , γ ′ k 1 j ) which transform respectively α k 1 j into αk 1 j and into αk 2 j . In fact, if the group k 1 of P is transformed into the group k 2 of P (instead of the group k 1 of P ) then necessarily the group k 2 of P is not transformed into the group k 2 of P , but into a group k 3 = k 2 of P ; The simplest solution

is k 3 = k 1 but it is not certain if K > 2.
Thus, the identifiability problem can be rewritten equivalently: There exists k 2 = k 1 and k 3 = k 2 such that, for any j and (δ k 1 j , λ j , γ k 1 j ), (δ k 2 j , λ j , γ k 2 j ) such that

Φ -1 ( αk 1 j ) = δ k 1 j Φ -1 (α k 1 j ) + λ j γ k 1 j (14) Φ -1 ( αk 2 j ) = δ k 2 j Φ -1 (α k 2 j ) + λ j γ k 2 j , (15)
there exists

(δ ′ k 1 j , λ ′ j , γ ′ k 1 j ) = (δ k 1 j , λ j , γ k 1 j ) and (δ ′ k 2 j , λ ′ j , γ ′ k 2 j ) = (δ k 2 j , λ j , γ k 2 j ) such that Φ -1 ( αk 2 j ) = δ ′ k 1 j Φ -1 (α k 1 j ) + λ ′ j γ ′ k 1 j (16) Φ -1 ( αk 3 j ) = δ ′ k 2 j Φ -1 (α k 2 j ) + λ ′ j γ ′ k 2 j . ( 17
)
It follows from ( 14) and ( 16)

Φ -1 ( αk 2 j ) = δ ′ k 1 j δ k 1 j Φ -1 ( αk 1 j ) -λ j γ k 1 j δ k 1 j δ ′ k 1 j + λ ′ j γ ′ k 1 j (18)
and similarly from ( 15) and ( 17)

Φ -1 ( αk 3 j ) = δ ′ k 2 j δ k 2 j Φ -1 ( αk 2 j ) -λ j γ k 2 j δ k 2 j δ ′ k 2 j + λ ′ j γ ′ k 2 j . (19)
Equations ( 18) and ( 19) lead to the following linear system

Φ(k 2 ,k 3 ) = Φ (k 1 ,k 2 ) δ,γ,λ,λ ′ × u (k 1 ,k 2 ) δ ′ ,γ ′ (20)
where Φ(k 2 ,k 3 ) = (Φ -1 ( αk 2 1 ), . . . , Φ -1 ( αk 2 d ), Φ -1 ( αk 3 1 ), . . . , Φ -1 ( αk

3 d )) ′ ∈ R 2d and Φ (k 1 ,k 2 ) δ,γ,λ,λ ′ and u (k 1 ,k 2 ) δ ′ ,γ ′
are respectively a matrix and a vector, with dimension depending on the model at hand, representing values of Φ -1 ( αk 1 j ), Φ -1 ( αk 2 j ),

δ k 1 j , δ k 2 j , γ k 1 j , γ k 2 j , λ j , λ ′ j for Φ (k 1 ,k 2 ) δ,γ,λ,λ ′ and values of δ ′ k 1 j , δ ′ k 2 j , γ ′ k 1 j , γ ′ k 2 j for u (k 1 ,k 2 ) δ ′ ,γ ′ .
Conditionally to the values of λ ′ j , the problem is identifiable if no solution exists to System (20) or, in other words, if the number of linearly independent lines of Φ

(k 1 ,k 2 ) δ,γ,λ,λ ′ is no less than the number of free parameters in u (k 1 ,k 2 ) δ ′ ,γ ′ . The number of lines of Φ (k 1 ,k 2 )
δ,γ,λ,λ ′ is equal to 2d and in most cases of interest these lines are independent for any values λ ′ j . See for instance the following example with the model [p k δ k γ k ] for a discussion on this subject. The number of free parameters of u

(k 1 ,k 2 ) δ ′ ,γ ′
corresponds to the values in Table 1 by artificially fixing K = 2 in this table. Except for model [p k δ j γ j ] the maximum number of free parameters is equal to max(4, d + 2). Thus, a sufficient condition of inter-group identifiability is 2d > max(4, d + 2). Since d > 2 usually for binary data, all models are identifiable except the model [p k δ j γ j ] since the number of free parameters is 2d. Although this model is non-identifiable, it can nevertheless be used by bearing in mind that a label switching may occur from P to P . Remark Only a single switch between 2 classes ((k 1 , k 2 ) of P becomes (k 2 , k 1 ) of P ) may occur with the model [p k δ j γ j ]. Indeed in case of a switch between more than 2 classes, for instances (k 1 , k 2 , k 3 ) of P becomes (k 2 , k 3 , k 1 ) of P or (k 1 , k 2 , k 3 , k 4 ) of P becomes (k 2 , k 1 , k 4 , k 3 ) of P , the number of equations reaches at least 3d for 2d free parameters: The model [p k δ j γ j ] is then identifiable.

Example with model

[p k δ k γ k ] for inter-group identifiability. In this situation Φ (k 1 ,k 2 ) δ,γ,λ,λ ′ = [M|N]
is a (2d × 4)-matrix formed by two matrix blocks M = (M ls ) 1≤s≤2

Φ (k 1 ,k 2 ) δ,γ,λ,λ ′ =            1 δ k 1 1 Φ -1 ( αk 1 1 ) -λ 1 γ k 1 1 δ k 1 1 0 λ ′ 1 0 . . . . . . 1 δ k 1 d Φ -1 ( αk 1 d ) -λ d γ k 1 d δ k 1 d 0 λ ′ d 0 0 1 δ k 2 1 Φ -1 ( αk 2 1 ) -λ 1 γ k 2 1 δ k 2 1 0 λ ′ 1 . . . . . . 0 1 δ k 2 d Φ -1 ( αk 2 d ) -λ d γ k 2 d δ k 2 d 0 λ ′ d           
The reader can easily convince himself that the 2d lines are generally linearly independent in practice. For instance, the simulation involved in the robustness study (Section 5) considers this model and none identifiability problem is encountered (the classification error rates are always strictly lower than 50% and thus none label switching has been obtained).

C Proof of the concavity of the function Q(θ, θ (q) )

The aim of this appendix is to prove that Q(θ, θ (q) ) is a strictly concave function of γ k , γ j and γ, and also of δ k , δ j , δ and λ j . The fundamental key of the proof is to study first concavity of Q(θ, θ (q) ) for γ kj . Let Q(γ kj ) be the function Q(θ, θ (q) ) with parameters pk , δ kj , λ j and θ (q) fixed, and let prove that Q(γ kj ) is strictly concave: (21)

Q(γ kj ) =
To prove that Q is strictly concave, it is sufficient to prove that both functions g 1 and g 2 are strictly positive:

• For all x ∈ R : g 1 (x) = xΦ(x) + φ(x) > 0, because lim x→-∞ g 1 (x) = 0 and g 1 is strictly increasing since g ′ 1 (x) = Φ(x) + xφ(x)xφ(x) = Φ(x) > 0,

• For all x ∈ R : g 2 (x) = xΦ(x)x + φ(x) > 0, because lim x→+∞ g 2 (x) = 0 and g 2 is strictly decreasing since g ′ 2 (x) = Φ(x) -1 + xφ(x)xφ(x) = Φ(x) -1 < 0.

Thus

∂ 2 Q(γ kj ) ∂γ 2 kj < 0 and Q(γ kj ) is strictly concave. If Q is no longer function of γ kj but now of γ k (respectively of γ j , of γ) the expression of the second derivative is the same as ( 21) by removing the sum on j (resp. on k, on (k, j)), and thus Q is still strictly concave. Consider now the function Q(δ kj ) with the above convention (p k , γ kj , λ j and θ (q) fixed). Same type of calculus leads to: and thus Q(δ kj ) is strictly concave. By using the above arguments, it arises immediately that Q(δ k ), Q(δ j ), Q(δ) are also strictly concave. Consider finally the function Q(λ j ). The second derivative is:

∂ 2 Q(δ kj
∂ 2 Q(λ j ) ∂λ 2 j = - ñ i=1 K k=1 t ik d j=1 γ 2 kj φ(ζ kj ) xij [Φ(ζ kj )] 2 g 1 (ζ kj ) + (1 -xij ) [1 -Φ(ζ kj )] 2 g 2 (ζ kj ) ,
and Q(λ j ) is strictly concave.

Figure 1 :Figure 2 :

 12 Figure 1: Density distribution of Y 1j |Z 1k =1 and Ỹ1j | Z1k =1 for the robustness study.

Figure 3 :

 3 Figure3: Classification error rate for usual discriminant analysis, generalized discriminant analysis (GDA) and clustering compared to the optimal error rate for different values of ǫ.

Figure 4 :

 4 Figure 4: Wings and tarsus lengths for diomedea and borealis.

  [p k δ k γ k ] for intra-group identifiability In this situation, Φ λ = [M|N] is a (Kd × 2K)-matrix formed by two block matrices M = (M lk ) 1≤k≤K 1≤l≤Kd and N = (N lk ) 1≤k≤K 1≤l≤Kd defined by

0

 0

  xij ) log 1 -Φ(ζ kj ) with ζ kj = δ kj Φ -1 (α kj ) + λ j γ kj . The derivative of the function Q(γ kj ) is: xij )λ j -φ(ζ kj ) 1 -Φ(ζ kj ) ,where φ is the probability density function of N (0, 1).Using that ∂φ(bx+a) ∂x = -b(bx + a)φ(bx + a) with a, b ∈ R, the second derivative is: -λ j ζ kj φ(ζ kj )Φ(ζ kj )λ j φ(ζ kj ) 2 [Φ(ζ kj )] 2 +(1xij ) λ j ζ kj φ(ζ kj )(1 -Φ(ζ kj ))λ j φ(ζ kj ) 2 [1 -Φ(ζ kj )] 2 , ζ kj )] 2 ( g 1 (ζ kj ) ζ kj Φ(ζ kj ) + φ(ζ kj )) + (1xij ) [1 -Φ(ζ kj )] 2 ( g 2 (ζ kj )ζ kj (φ(ζ kj ) -1) + φ(ζ kj )) .

  ζ kj )] 2 g 1 (ζ kj ) + (1xij ) [1 -Φ(ζ kj )] 2 g 2 (ζ kj ) ,

  1j |Z 1k =1 , 0 or 1, corresponds a positive value of Y 1j |Z 1k =1 , and then to avoid binary variables to inherit from the natural order induced by continuous variables. Moreover, the joint distribution of Y 1 |Z 1k =1 = (Y 11 |Z 1k =1 , . . . , Y 1d |Z 1k =1

  whose distribution is a mixture of two Gaussians. The mixture density distribution of Y 1j |Z 1k =1 and Ỹ1j | Z1k =1 (for fixed j and k), and the Gaussian corresponding density (with the same moments as Y 1j |Z 1k =1 and Ỹ1j | Z1k =1

Table 2 :

 2 Classification error rates (%) and value of the BIC criterion for test population diomedea with training population borealis. model

Acknowledgements

Authors thank Vincent Bretagnolle of CEBC-CNRS for providing the biological data.