Romain Reocreux

Minh Huynh

Carine Michel
email: carine.michel@ens-lyon.fr

Philippe Sautet
email: philippe.sautet@ens-lyon.fr

Romain Réocreux

Controlling the Adsorption of Aromatic Compounds on Pt(111) with Oxygenate Substituents: From DFT to Simple Molecular Descriptors

Keywords: aromatic, Pt(111), adsorption, DFT, HSAB, dispersion, deformation

  L'archive ouverte pluridisciplinaire

The adsorption of substituted aromatics at metallic surfaces is of great interest in a wide range of domains, from corrosion protection [START_REF] Gattinoni | Understanding corrosion inhibition with van der Waals DFT methods: the case of benzotriazole[END_REF][START_REF] Özcan | Insights into surfaceadsorbate interactions in corrosion inhibition processes at the molecular level[END_REF][START_REF] Obot | Density functional theory (DFT) as a powerful tool for designing new organic corrosion inhibitors. Part 1: An overview[END_REF] to the degradation of pollutants. Among them, oxygenates are of particular interest. Besides some examples of phenol synthesis from the heterogeneous oxidation of benzene, [START_REF] Orita | Simulation of phenol formation from benzene with a Pd membrane reactor: Ab initio periodic density functional study[END_REF][START_REF] Niwa | A One-Step Conversion of Benzene to Phenol with a Palladium Membrane[END_REF] most of the studies reported in the literature have focused on their decomposition. Phenolic compounds are indeed very toxic polluants [START_REF] Nousir | Deactivation phenomena during catalytic wet air oxidation (CWAO) of phenol over platinum catalysts supported on ceria and ceriazirconia mixed oxides[END_REF] (e.g. dyes and endocrine disruptors) and, at the same time, commonly used as models to understand the decomposition and deoxygenation of lignin, a biomass polymer that might become a sustainable source of aromatics in the future. [START_REF] Zakzeski | The Catalytic Valorization of Lignin for the Production of Renewable Chemicals[END_REF][START_REF] Nimmanwudipong | Catalytic Conversion of Guaiacol Catalyzed by Platinum Supported on Alumina: Reaction Network Including Hydrodeoxygenation Reactions[END_REF][START_REF] Hensley | Adsorption of guaiacol on Fe (110) and Pd (111) from first principles[END_REF][START_REF] Gao | Chemo-and Regioselective Hydrogenolysis of Diaryl Ether CO Bonds by a Robust Heterogeneous Ni/C Catalyst: Applications to the Cleavage of Complex Lignin-Related Fragments[END_REF] Those transformations can be performed in reducing or oxidative media using a metallic surface, the catalytic activity of which is related to the adsorption energies of the intermediates as the Sabatier principle suggests. This has been developed quantitatively by the mean of volcano plots and Brønsted-Evans-Polanyi linear relationships. [START_REF] Nørskov | Towards the computational design of solid catalysts[END_REF] Using such tools, one can rapidly predict the activity of a catalyst knowing how it interacts with the reactant or the intermediates. This is all the more important that today there is a need in rapid screening of catalysts to optimize industrial processes.

Theoretical chemistry and surface science are important methods to understand how those particular chemicals adsorb on metallic surfaces. Vibrational spectroscopies [START_REF] Lehwald | Vibration spectroscopy of benzene adsorbed on Pt(111) and Ni(111)[END_REF][START_REF] Tan | Catalyst-induced changes in a substituted aromatic: A combined approach via experiment and theory[END_REF] and Scanning Tunneling Microscopy [START_REF] Weiss | Site Dependence of the Apparent Shape of a Molecule in Scanning Tunneling Microscope Images: Benzene on Pt(111)[END_REF][START_REF] Sautet | A theoretical analysis of the site dependence of the shape of a molecule in STM images[END_REF] provide structural features while Temperature Program Desorption (TPD) [START_REF] Ihm | Stepwise dissociation of thermally activated phenol on Pt(111)[END_REF] and microcalorimetry experiments [START_REF] Ihm | Calorimetric Measurement of the Heat of Adsorption of Benzene on Pt(111)[END_REF] allow us to quantify the affinity of the aromatic adsorbate to the metallic surface by the mean of heats of adsorption. On the theoretical side, Density Functional Theory (DFT) is now widely used to characterize adsorption geometries with their computed adsorption energies. The early studies by Morin et al. [START_REF] Morin | Chemisorption of Benzene on Pt(111), Pd(111), and Rh(111) Metal Surfaces : A Structural and Vibrational Comparison from First Principles[END_REF] determined the predominant structures of benzene on Pt(111) even though the GGA (generalized gradient approximation) calculations markedly underestimate the experimental adsorption energies. Many other aromatics on different metal surfaces have been reported with such method. [START_REF] Ghiringhelli | Phenol near Ni(111), Ni(110), and Ni(221) surfaces in a vertical ring geometry: A density functional study of the oxygensurface bonding and O-H cleavage[END_REF][START_REF] Delle Site | Adsorption energies and geometries of phenol on the (111) surface of nickel: An ab initio study[END_REF][START_REF] Bonalumi | Theoretical and Spectroscopic Study of the Effect of Ring Substitution on the Adsorption of Anisole on Platinum[END_REF][START_REF] Honkela | Computational study of the adsorption and dissociation of phenol on Pt and Rh surfaces[END_REF] Since then, more advanced DFT methods have been developed to account accurately for dispersion interactions. [START_REF] Steinmann | Comprehensive benchmarking of a densitydependent dispersion correction[END_REF][START_REF] Klimeš | Perspective: Advances and challenges in treating van der Waals dispersion forces in density functional theory[END_REF] Some of those functionals provide a much better estimation of the aromatic adsorption energies. [START_REF] Yildirim | Trends in Adsorption Characteristics of Benzene on Transition Metal Surfaces: Role of Surface Chemistry and van der Waals Interactions[END_REF][START_REF] Gautier | Molecular adsorption at Pt(111). How accurate are DFT functionals?[END_REF] Other DFT methods, such as the Random Phase Approximation, are known to be among the most accurate calculations for adsorption on metals but remain computationally extremely expensive. [START_REF] Schimka | Accurate surface and adsorption energies from many-body perturbation theory[END_REF] Despite that remarkable evolution toward experimental and theoretical accuracy over the years, all those techniques are very timedemanding and they rarely provide chemical trends that can easily be exploited by both experimentalists and theoreticians. [START_REF] Bonalumi | Theoretical and Spectroscopic Study of the Effect of Ring Substitution on the Adsorption of Anisole on Platinum[END_REF][START_REF] Peköz | Tuning the Adsorption of Aromatic Molecules on Platinum via Halogenation[END_REF][START_REF] Hensley | Adsorption of phenol on Fe (110) and Pd (111) from first principles[END_REF] Notwithstanding, there are many examples of simple aromatic descriptors widely used to account for various chemical properties. Those can be geometry descriptors (Holleman's product rule) or electronic descriptors (Hammett parameters or frontier orbital and Pearson's descriptors). They do not aim at describing the properties very accurately but they definitely are very useful and popular rules of thumb.

In this work, we investigate the adsorption of a range of alkyl and alkoxy substituted aromatics (see Figure 1) on a Pt(111) surface with periodic PBE-dDsC calculations (see SI for a description of methods). Beyond the DFT numbers, we propose simple trends that enable the understanding and the prediction of the influence of substituents on aromatic adsorption. In the set, there are different kinds of aromatic compounds, mostly biomassderivated oxygenates, the substituents of which are mainly linked to the aromatic ring via an oxygen or a carbon atom. Therefore we will distinguish between O-substituents and Csubstituents. To that set, we added toluene and benzene: the former to isolate the effect of Csubstituents, and the latter to compare with the simplest aromatic chemical.

The very studied benzene on Pt(111), albeit not substituted, is a clear starting point. It predominantly adsorbs at bri30 and hcp0 sites (see Figure 2), the later being the most stable. [START_REF] Morin | Chemisorption of Benzene on Pt(111), Pd(111), and Rh(111) Metal Surfaces : A Structural and Vibrational Comparison from First Principles[END_REF][START_REF] Yildirim | Trends in Adsorption Characteristics of Benzene on Transition Metal Surfaces: Role of Surface Chemistry and van der Waals Interactions[END_REF][START_REF] Gautier | Molecular adsorption at Pt(111). How accurate are DFT functionals?[END_REF] In that nomenclature the position of the centre of the flat-lying aromatic adsorbate is first given (hcp, bri, fcc or top) followed by the ide- alized value of the angle between the aromatic C-C direction and the Pt(111) close-packed orientation (see Figure 2). Upon adsorption the ring gets slightly deformed and the hydrogens are bent away from the surface. To account for the latter effect, we introduce the out-ofplane deformation angle Θ (see Scheme 1 for the graphical definition). At the hcp0 site, the aromatic carbons are strictly equivalent: they are bound -in pairs -to the surface in a way that is very similar to π-bonded ethylene (with Θ = 20 • ). They will therefore be labelled π. At the bri30 site, we can find two pairs of π carbons as well (with Θ = 18 • ) but also two other carbons -in para position to each other -that show a bigger deformation angle (Θ = 34 • ) and a much more tetrahedral geometry as σ-bonded carbonated species would. They will therefore be labelled σ. What about the substituted aromatics? Where shall we put the substituent(s) to get the most stable structure? On π or σ carbons ? Do(es) the substituent(s) influence the position of the aromatic ring on the surface? To answer these questions, we limited the DFT inspection to those most stable bri30 and hcp0 sites since the other adsorption modes are at least about 0.3 eV less stable. Let's start with two simple monosubstituted aromatic compounds: an O-substituted one, namely anisole, and a C-substituted one, namely toluene. We have to consider two geometries at the bri30 site, the substituent being either on a π or an σ carbon, and one at the hcp0 site. Anisole and toluene show the same preference for the bri30 adsorption site over hcp0 with an energy difference of about 0.3 eV (see Table 1). This energy difference does not seem to be related to any increased deformation (similar deformation angle as in the most stable bri30 geometry) but rather a commun feature to many aromatic compounds. The same energy difference was indeed computed for benzene by Morin et al. [START_REF] Morin | Chemisorption of Benzene on Pt(111), Pd(111), and Rh(111) Metal Surfaces : A Structural and Vibrational Comparison from First Principles[END_REF] and the bri30 structure was also reported for anisole [START_REF] Bonalumi | Competition at Chiral Metal Surfaces: Fundamental Aspects of the Inversion of Enantioselectivity in Hydrogenations on Platinum[END_REF] and other aromatic compounds such as phenol, [START_REF] Honkela | Computational study of the adsorption and dissociation of phenol on Pt and Rh surfaces[END_REF] guaiacol [START_REF] Lee | Guaiacol Hydrodeoxygenation Mechanism on Pt(111): Insights from Density Functional Theory and Linear Free Energy Relations[END_REF][START_REF] Lu | Theoretical Investigation of the Reaction Mechanism of the Guaiacol Hydrogenation over a Pt(111) Catalyst[END_REF] and halogenobenzenes. [START_REF] Peköz | Tuning the Adsorption of Aromatic Molecules on Platinum via Halogenation[END_REF] Regarding substitution at the bri30 adsorption site, toluene does not really show any preference for a π or σ position even if the deformation angle is almost doubled between the first and the second geometry. Unlike toluene, anisole is much more sensitive to the substitution position: with a very similar variation of the deformation angle, the σ position is less stable than the π position by 0.11 eV. From those observations, one can go further and search for simple molecular descriptors giving a picture of the adsorption at the most stabilizing bri30 site. Since deformation seems to play a role, we will decompose the adsorption energy into three terms as follows:

E ads = E def + E φ + E χ (1)
where E def , E φ and E χ stand for the deformation energy, the physical (dispersion interaction) and the chemical part of the total interaction energy respectively (see the associated Hess cycle in Scheme 1). What molecular descriptors do those energetic contributions depend on?

Dispersion interaction. The dispersion interaction is known to be proportional to the polarizability of the molecules. According to the Clausius-Mossotti relation, 33 the polarizability is proportional to the molecular weight. From semi-classical physics (see details in SI), we can estimate the dispersion interaction energy of a compound (molecular weight M at a height of Z) with a metallic surface as follows:

E disp -k M Z 3 (2)
with k 0.3 eV Å3 mol g -1 for typical organic compounds. For the aromatic compounds we considered (molecular weight M ), we have therefore the contributions of both the aromatic ring (molecular weight M 0 , height Z 0 3 Å) that lies closer to the surface than the substituent(s) (molecular weight M , height Z 3.5 Å). Adding the two contributions and considering that M = M -M 0 , we get:

E disp = -k M 0 Z 3 0 + M Z 3 = - k Z 3 × M -k M 0 Z 3 0 - M 0 Z 3 (3)
We obtain from equation (3) the relation E disp = aM + b with a = -7 • 10 -3 eV mol/g and b = -0.3 eV. That agrees very well with the DFT PBE-dDsC results that show a very similar linear relationship (E disp in eV and M in g/mol, Figure 3a) with a very small mean absolute deviation (MAD):

E disp = -6.7 • 10 -3 M -0.25 MAD = 0.04 eV
Although dispersion interactions are important to estimate properly the adsorption energies of flat-flying conjugated molecules, [START_REF] Gautier | Molecular adsorption at Pt(111). How accurate are DFT functionals?[END_REF] it does interestingly not depend on the π/σ distinction and probably explains why GGA calculations are still relevant for energy differences but not for absolute energies. Deformation energies. The deformation energy should depend on the total number of substituents n since there are more functional groups to be bent from the ring plane. That is why we plotted the deformation energy as a function of the total number of substituents (Figure 3b). In order to account for the sensitivity of O-substituents to the π/σ positions, the whole set was divided into two sub-sets: the first sub-set consisting of adsorbates without any O-substituents at the σ positions (no constraints on the C-substituents though) and the second sub-set consisting of adsorbates with one, and one only, O-substituent at a σ position. For the former -referred to as O 0 σ -we can see that the deformation energy (in eV) increases linearly with the total number of substituents n:

E def (O 0 σ ) = 0.19n + 1.88 MAD = 0.04 eV
for the later -referred to as O 1 σ -the same behaviour is observed :

E def (O 1 σ ) = 0.21n + 2.21 MAD = 0.06 eV
It means there is the same deformation energy cost of ca. 0.2 eV per substituent and it is ad-That was already pointed out by Peköz et al. [START_REF] Peköz | Tuning the Adsorption of Aromatic Molecules on Platinum via Halogenation[END_REF] who studied different polyhalogenobenzenes on Pt(111): the deformation being about 0.3 eV per substituent. However they did not distinguish between π and σ positions, the influence of which is a δe π/σ = 0.33 eV off-set for O-substituents.

Chemical interaction energies. To rationalise chemical interactions, there are several approaches. One of them is Pearson's Hard and Soft Acids and Bases (HSAB) theory. [START_REF] Parr | Absolute hardness: companion parameter to absolute electronegativity[END_REF][START_REF] Pearson | Hard and soft acids and bases, HSAB, part 1: Fundamental principles[END_REF] This theory describes the interaction between two species A and B in terms of their chemical hardness η A and η B and their electron chemical potential µ A and µ B that can be seen as the opposite of their electronegativity. The interaction induces an electronic flow (from high to lower chemical potentials) associated with the following HSAB stabilisation energy E:

E = - 1 4 (µ A -µ B ) 2 η A + η B < 0 (4)
Within the finite difference approximation and the frontier orbital approximation, Pearson proposed to estimate the chemical hardness η i and the electron chemical potential µ i as follows:

µ i 1 2 (E LU M O + E HOM O ) for molecules ε F
for solids (5)

η i 1 2 (E LU M O -E HOM O ) for molecules 1 2 BG
for solids (6) where E LU M O , E HOM O , ε F and BG stand for the energies of the lowest unoccupied molecular orbital energy, the highest occupied molecular orbital energy, the Fermi level and the band gap respectively. For a metal, the band gap is zero. Therefore the HSAB interaction energy simply writes:

E = - (µ -ε F ) 2 4η < 0 (7)
where µ and η stand for the adsorbate parameters. Literature shows examples where this interaction energy E is used quantitatively. [START_REF] Romero | The Local HSAB Principle and Bond Dissociation Energy of p -Substituted Phenol[END_REF] Such a quantitative treatment is however accurate only if the geometries of the two interacting parts do not change significantly upon interaction. Otherwise HSAB theory gets much more complex to account for the external potential variations but still allows -in its simpler form -to understand chemical interaction at metal surfaces. [START_REF] Kokalj | On the HSAB based estimate of charge transfer between adsorbates and metal surfaces[END_REF] For the purpose of simplicity (and predictability) we will correlate the actual DFT interaction energy with its HSAB counterpart evaluated from the totally relaxed gas phase structures (µ 0 and η 0 ) referred to as E 0 :

E 0 = - (µ 0 -ε F ) 2 4η 0 < 0 ( 8
)
where both (µ 0 -ε F ) and η 0 are positive quantities. Plotting the DFT chemical interaction energies as a function of the HSAB interaction energies (Figure 3c), we get for the first sub-set O 0 σ :

E χ (O 0 σ ) = 0.45 • E 0 -3.14 MAD = 0.05 eV
and for the second sub-set O 1 σ , which consists of more deformed adsorbates:

E χ (O 1 σ ) = 1.20 • E 0 -2.92 MAD = 0.13 eV
Because of the approximations we did, we expected the DFT chemical interaction energies E χ to be different from the HSAB idealized interaction energies E 0 . Despite of that, we interestingly get a fair correlation with reasonable MAD and a positive slope (the bigger E 0 , the bigger E χ ) that is still of the order of 1. As for the intercept, it corresponds to the interaction energy between an adsorbate and a surface of identical electronegativities (or electron chemical potentials) µ 0 = ε F . In other words, the intercept is the covalent contribution to the overall chemical interaction and in the two subsets they are quite comparable (about 3 eV, hence 0.5 eV per aromatic carbon). Moreover, the chemical interaction in the more deformed O 1 σ sub-set is stronger than the less deformed O 0 σ sub-set. It can be interpreted in terms of chemical hardness. The Principle of Maximum Hardness 38 suggests that the gas phase structures maximize the chemical hardness: thus the molecules should get softer upon deformation. Since the metallic surface is extremely soft (η P t = 0), the deformation process indeed increases the intensity of the soft-soft interactions. Nonetheless, we showed that the deformation had its own energetic cost. Altogether the chemical part of the interaction energy favors the presence of an O-substituent at the σ position, and therefore competes with the effect of the deformation energy as discussed in the next section.

Prediction of adsorption energies and structures. By summing the fitted linear expressions obtained in the previous parts, one can give an estimation of the adsorption energies E ads of the two sub-sets O 0 σ and O 1 σ :

E ads (O 0 σ ) = αM + 0.45E 0 + 0.19n -1.51 (9)
E ads (O 1 σ ) = αM + 1.20E 0 + 0.21n -0.95 (10) with α = -6.7 • 10 -3 eV mol/g. Taking into account the molecular weight M , the total number of substituents n and the electronic effects through E 0 , our approach remarkably well ap- proximates the actual adsorption energies computed with DFT (see the parity plot given in Figure 4). The overall mean absolute deviation is only 0.08 eV with a maximum absolute deviation of 0.25 eV, so our approach appears quite reliable. Thus, those relations are a fast screening tool for the adsorption of substituted aromatic oxygenates, predicting the energy and the structure, in line with the widely used scaling relations [START_REF] Lee | Guaiacol Hydrodeoxygenation Mechanism on Pt(111): Insights from Density Functional Theory and Linear Free Energy Relations[END_REF][START_REF] Ferrin | Modeling ethanol decomposition on transition metals: a combined application of scaling and Brønsted-Evans-Polanyi relations[END_REF][START_REF] Liu | Decomposition pathways of glycerol via C-H, O-H, and C-C bond scission on Pt(111): A density functional theory study[END_REF][START_REF] Sutton | A Theoretical and Computational Analysis of Linear Free Energy Relations for the Estimation of Activation Energies[END_REF][START_REF] Calle-Vallejo | Finding optimal surface sites on heterogeneous catalysts by counting nearest neighbors[END_REF][START_REF] Calle-Vallejo | Introducing structural sensitivity into adsorptionenergy scaling relations by means of coordination numbers[END_REF] More interestingly, our model also seems to depict -at least partly -the balance between the deformation energy (which tends to favor adsorption configuration without O substituent in the σ position) and the chemical interaction energy (which tends to favor adsorption modes with a O substituent at the σ position). The energy difference ∆ 0/1 σ writes: ∆ 0/1 σ = E ads (O 1 σ ) -E ads (O 0 σ ) = 0.75E 0 + 0.02n + 0.56 (11) The HSAB parameter E 0 being negative, two regimes are possible. At small E 0 , the most favorable adsorption mode is the O 0 σ , with no O susbstituent at the σ position when possible. For large enough absolute values of E 0 , one should in principle see an inversion of this trend. This is indeed the case for MG that has the largest HSAB idealized interaction energy (E 0 = -0.94 eV): the O 1 σ geometry is slightly more stable than the O 0 σ geometry with an energy difference of 0.01 eV (see SI for details).

Our correlations do not only allow to predict the adsorption energies and geometries from very simple molecular descriptors but they also help to understand the nature of the bond between aromatic oxygenates and Pt(111), a very typical surface for heterogeneous catalysis applications. It opens new strategies to rationalize chemical reactivities or the design of metallic catalysts since we showed that besides the Fermi level of the metal (or the related dband centre), other parameters control adsorption such as the electronegativity (or electron chemical potential), the chemical hardness, the molecular weight and the adsorption geometry.

Figure 2 :

 2 Figure 2: The two most stable adsorption sites of benzene: hcp0 on the right (all the substitution positions being equivalent) and bri30 on the left (with the two substitution positions: σ and π).

Scheme 1 :

 1 Scheme 1: Hess cycle for anisole adsorption involving the deformations of the molecule (quantified by the out-of-plane deformation angle Θ) and the slab.

Figure 3 :

 3 Figure 3: Energy decomposition over the set of molecules defined in Figure 1: (a) physical interaction (b) deformation energy (c) chemical interaction energy. In the last two diagrams: × adsorption mode with all hydroxy/alkoxy groups at a π position (sub-set O 0 σ ); + adsorption mode with one hydroxy/alkoxy group at a σ position (sub-set O 1 σ ).

Figure 4 :

 4 Figure 4: Descriptor-based adsorption energies as a function of the actual DFT computed adsorption energies. The full line corresponds to the identity.

Table 1 :

 1 Energy and geometry features of anisole and toluene structures on Pt(111).

	Adsorption site	bri30	hcp0
	Substituent position π	σ	π
	Anisole
	Deformation angle Θ 24 • 44 •	24 •
	Relative stability (eV) 0	0.11 0.29
	Toluene
	Deformation angle Θ 26 • 40 •	27 •
	Relative stability (eV) 0 0.002 0.25
	Benzene
	Deformation angle Θ 18 • 34 •	20 •
	Relative stability (eV)	0	0.26

Acknowledgements

The authors are grateful to the Centre Blaise Pascal (CBP), the Pôle Scientifique de Modélisation Numérique (PSMN) -both at the École Normale Supérieure de Lyon -and the GENCI (CINES), project x2015080609 that provided assistance and calculation resources. They also thank the LIA FUNCAT collaboration for fundings.

Supporting

Information Available: Computational details, pictures and full energy analysis of all the adsorbates mentioned, and derivation of equation ( 2). This material is available free of charge via the Internet at http://pubs.acs.org/.


Perrine Soret

Marta Fernandez Avalos

Linda Wittkop

Daniel Commenges

Rodolphe Thiébaut

Lasso regularization for left-censored Gaussian outcome and high-dimensional predictors

Keywords: Limit of detection, Buckley-James least squares procedure, HIV viral load, Drug resistance, HIV genotypic mutations, Cross-sectional studies

Background

Left-censoring due to the lower detection limit of an assay is a common problem in many fields including biology, chemistry, and the environmental sciences. One example is the quantification of the human immunodeficiency virus (HIV) viral load in plasma. The sensitivity of assays has improved and the detection threshold has decreased from 10,000 copies/mL to 20 or fewer copies/mL today. Several statistical methods have been proposed to account for left-censoring of such quantitative variables in crosssectional (with one measure per subject) and longitudinal (with several measures per subject) studies. Standard methods include multiple imputation [START_REF] Paxton | Longitudinal analysis of quantitative virologic measures in human immunodeficiency virus-infected subjects with > or = 400 CD4 lymphocytes: implications for applying measurements to individual patients. National Institute of Allergy and Infectious Diseases AIDS Vaccine Evaluation Group[END_REF][START_REF] Helsel | More than obvious: Better methods for interpreting nondetect data[END_REF][START_REF] Lee | Multiple imputation for left-censored biomarker data based on Gibbs sampling method[END_REF][START_REF] Greco | Bayesian analysis of censored response data in family-based genetic association studies[END_REF], reverse survival analysis methods [START_REF] Helsel | More than obvious: Better methods for interpreting nondetect data[END_REF][START_REF] Marschner | Clinical trials using HIV-1 RNA-based primary endpoints: Statistical analysis and potential biase[END_REF][START_REF] Gillespie | Estimating population distributions when some data are below a limit of detection by using a reverse Kaplan-Meier estimator[END_REF][START_REF] Dinse | Accomodating measurements below a limit of detection: A novel application of Cox regression[END_REF], quantile regression [START_REF] Wang | Quantile regression in partially linear varying coefficient models[END_REF][START_REF] Eilers | Quantile regression for the statistical analysis of immunological data with many non-detects[END_REF] and censored quantile regression [START_REF] Powell | Least absolute deviations estimation for the censored regression model[END_REF][START_REF] Powell | Censored regression quantiles[END_REF]. Furthermore, the Tobit model with censored outcome which is supposed to be normally distributed can be estimated by maximum likelihood [START_REF] Tobin | Estimation of relationships for limited dependent variables[END_REF][START_REF] Hughes | Mixed effects models with censored data with application to HIV RNA levels[END_REF][START_REF] Jacqmin-Gadda | Analysis of left-censored longitudinal data with application to viral load in HIV infection[END_REF][START_REF] Lynn | Maximum likelihood inference for left-censored HIV RNA data[END_REF][START_REF] Nie | Linear regression with an independent variable subject to a detection limit[END_REF][START_REF] Fu | A comparative investigation of methods for longitudinal data with limits of detection through a case study[END_REF][START_REF] Wiegand | Comparison of models for analyzing two-group, cross-sectional data with a gaussian outcome subject to a detection limit[END_REF] or by the Buckley-James estimator [START_REF] Wiegand | Comparison of models for analyzing two-group, cross-sectional data with a gaussian outcome subject to a detection limit[END_REF][START_REF] Buckley | Linear regression with censored data[END_REF]. Indeed, HIV viral load appears to have an underlying Gaussian distribution truncated by the detection limit that justifies the normality hypothesis [13-15, 17, 18]. As expected, approaches accounting for left-censoring outperform simple imputation of a constant [2, 4, 13-16, 18, 20-22].

Another issue may arise when the number of predictors (p) is high compared to the number of statistical units (n), without excluding the possibility that n < p. This is known as high dimensionality. In the context of HIV infection, this can be illustrated by analyzing the association between the presence of HIV mutations and the response to antiretroviral therapy which is measured by HIV viral load. HIV strains circulating in a given individual can present mutations associated with antiretroviral treatment failure (detectable HIV viral load), also called HIV drug resistance mutations. Thus, genotypic tests allowing the detection of HIV drug resistance mutations are commonly performed in patients starting a new antiretroviral regimen or even in newly HIV-infected patients because of the transmission of resistant strains [START_REF] Hirsch | Antiretroviral drug resistance testing in adult HIV-1 infection: 2008 recommendations of an International AIDS Society-USA panel[END_REF][START_REF] Wittkop | Effects of transmitted drug resistance on virological and immunological response to initial combination antiretroviral therapy for HIV (euro-coord-chain joint project): a european multicohort study[END_REF][START_REF] Hofstra | Transmission of HIV drug resistance and the predicted effect on current first-line regimens in europe[END_REF][START_REF] Wensing | update of the drug resistance mutations in HIV-1[END_REF]. Lasso linear [START_REF] Rabinowitz | Accurate prediction of HIV-1 drug response from the reverse transcriptase and protease amino acid sequences using sparse models created by convex optimization[END_REF][START_REF] Beerenwinkel | The individualized genetic barrier predicts treatment response in a large cohort of HIV-1 infected patients[END_REF] and logistic regressions [START_REF] Cozzi-Lepri | Can linear regression modeling help clinicians in the interpretation of genotypic resistance data? an application to derive a lopinavir-score[END_REF], principal component and partial least square logistic regressions [START_REF] Wittkop | Alternative methods to analyse the impact of HIV mutations on virological response to antiviral therapy[END_REF], and multiple testing correction [START_REF] Assoumou | Standardization and clinical relevance of HIV drug resistance testing project from the forum for collaborative HIV research. Relative contributions of baseline patient characteristics and the choice of statistical methods to the variability of genotypic resistance scores: the example of didanosine[END_REF] have been used to deal with more than 100 predictors and fewer than a few hundred of patients, a common situation in this context [START_REF] Rhee | Genotypic predictos of human immunodeficiency cirus type 1 drug resistance[END_REF].

These studies use a dichotomized outcome or simple imputation by a constant to circumvent the problem of censoring. One limitation of dichotomizing a continuous outcome is the loss of information and hence power. In addition, success is usually defined as achieving an undetectable HIV viral load. However, the detection limit, although not random, depends on several factors that differ from one study to another. Thus, there is no reason, except convenience, for the detection limit to correspond to the threshold for dichotomization.

We hypothesize that approaches accounting for leftcensoring will exhibit better results compared to simple imputation strategies in a high-dimensional setting similar to what has been found in low-dimensional settings.

Some works have simultaneously addressed both censoring and high-dimensional problems using the Lasso [START_REF] Tibshirani | The Lasso method for variable selection in the Cox model[END_REF][START_REF] Huang | Regularized estimation in the accelerated failure time model with high-dimensional covariates[END_REF][START_REF] Datta | Predicting patient survival from microarray data by accelerated failure time modeling using partial least squares and Lasso[END_REF][START_REF] Johnson | Variable selection in semiparametric linear regression with censored data[END_REF][START_REF] Wang | Doubly penalized Buckley-James method for survival data with high-dimensional covariates[END_REF][START_REF] Cai | Regularized estimation for the accelerated failure time model[END_REF][START_REF] Ueki | A note on automatic variable selection using smooth-threshold estimating equations[END_REF][START_REF] Wang | Buckley-james boosting for survival analysis with high-dimensional biomarker data[END_REF][START_REF] Shows | Sparse estimation and inference for censored median regression[END_REF][START_REF] Wang | Variable selection for censored quantile regression[END_REF][START_REF] Chung | A tutorial on rank-based coefficient estimation for censored data in small-and large-scale problems[END_REF], partial least squares [START_REF] Huang | Modeling the relationship between LVAD support time and gene expression changes in the human heart by penalized partial least squares[END_REF], random forests [START_REF] Ishwaran | Random survival forests[END_REF], support vector machines [START_REF] Wang | Support vector hazards machine: A counting process framework for learning risk scores for censored outcomes[END_REF], and deep learning [START_REF] Van Der Burgh | Deep learning predictions of survival based on MRI in amyotrophic lateral sclerosis[END_REF]. These examples were developed for right-censored survival data. A main approach to left-censored data analysis is based on methods typically used with rightcensored survival data such as the Buckley-James estimator. Left-censored data are then previously reversed to right-censored data. While from a statistical point of view, the nature of the outcome (time-to-event or quantitative measurement below a limit of detection) is secondary, this can impact the choice of adequate probability distribution functions and other practical issues.

We propose a Lasso-regularized Buckley-James least squares method with both, non-parametric imputation using Kaplan-Meier and parametric imputation based on the Gaussian distribution. The non-parametric Buckley-James estimator, which simply replaces censored residuals by their conditional expectations in an iterative way, has been previously applied to left-censored HIV viral load data in a cross-sectional study [START_REF] Wiegand | Comparison of models for analyzing two-group, cross-sectional data with a gaussian outcome subject to a detection limit[END_REF]. On the other hand, the Lasso extension of the non-parametric Buckley-James method has been proposed for right-censored data [START_REF] Johnson | Variable selection in semiparametric linear regression with censored data[END_REF][START_REF] Cai | Regularized estimation for the accelerated failure time model[END_REF][START_REF] Wang | Buckley-james boosting for survival analysis with high-dimensional biomarker data[END_REF][START_REF] Johnson | On lasso for censored data[END_REF]. Our contribution consists in using the latter method for left-censored outcomes and highdimensional predictors. Furthermore, we propose an original parametric version of the Buckley-James method, which is adapted to the typical assumption of a Gaussian distribution of HIV viral load. We demonstrate the value of these approaches by comparing them to Lasso linear regression with simple imputation [START_REF] Beerenwinkel | The individualized genetic barrier predicts treatment response in a large cohort of HIV-1 infected patients[END_REF] for predicting the response to antiretroviral therapy by HIV genotypic mutations.

Our primary objective is to predict as accurately as possible responses in future patients who will switch to a similar regimen. Thus, comparisons are based on mean square prediction error. The prediction performances of the different methods were assessed on simulated data that reproduced the observed data. Then, methods were applied to data obtained in a collaborative study from clinical trials and cohorts provided by the Standardization in Clinical Relevance of HIV Drug Resistance Testing Project from the Forum for Collaborative HIV Research [START_REF] Cozzi-Lepri | Initiatives for developing and comparing genotype interpretation systems: external validation of existing rule-based interpretation systems for abacavir against virological response[END_REF]. The actual data presented a moderate censoring rate of 26 %, i.e. a realistic magnitude [START_REF] Wiegand | Comparison of models for analyzing two-group, cross-sectional data with a gaussian outcome subject to a detection limit[END_REF][START_REF] Marks | The spectrum of engagement in HIV care: do more than 19% of HIV-infected persons in the US have undetectable viral load?[END_REF]. However, high (around 50 %) or even severe (around 70 %) censoring rates could be observed in older studies with a high limit of detection (LOD) or particular populations with low treatment failure rate, e.g. HIV controllers [START_REF] Wiegand | Comparison of models for analyzing two-group, cross-sectional data with a gaussian outcome subject to a detection limit[END_REF][START_REF] Dao | Study to Understand the Natural History of HIV and AIDS in the Era of Effective Therapy (SUN) Investigators. Low vitamin D among HIV-infected adults: prevalence of and risk factors for low vitamin D levels in a cohort of HIV-infected adults and comparison to prevalence among adults in the US general population[END_REF][START_REF] Leon | Rate and predictors of progression in elite and viremic HIV-1 controllers[END_REF]. Thus, we also explored the impact of high and severe censoring rates on performance.

We detail how to use publicly available R packages to compute Lasso estimates with left-censored data.

Finally, we discuss possible extensions and applications of our work.

Methods

Methods to analyze left-censored outcome

In this section, we review the simplest models and estimation methods used to deal with left-censoring in cross-sectional studies. For a more extensive and comprehensive review of these methods, see [START_REF] Helsel | More than obvious: Better methods for interpreting nondetect data[END_REF][START_REF] Wiegand | Comparison of models for analyzing two-group, cross-sectional data with a gaussian outcome subject to a detection limit[END_REF]. Thereafter, we consider the Lasso extension of those methods that support simple implementations.

The linear model

First, consider the general linear regression model

Y i = X i β + ε i , i = 1, • • • , n (1)
where X i is a p-vector of fixed predictors, Y i is the uncensored continuous random outcome variable, β = β 1 , • • • , β p is a p-vector of unknown regression parameters and ε i are independent and identically normally distributed random variables with mean 0 and constant variance σ 2 . Let X be the n × p matrix X = (X 1 , . . . , X n ) and Y the n × 1 vector Y = (Y 1 , . . . , Y n ) . The intercept is omitted in the model for simplicity, and all predictor variables are assumed to be standardized (i.e. zero mean and unit variance).

Lasso on complete data

The Lasso (Least Absolute Shrinkage and Selection Operator) [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF] is one of the most popular methods in high-dimensional data analyses. It allows for simultaneous estimation and variable selection and has efficient algorithms available. It is considered here as the Gold Standard for our simulation studies. The Lasso estimator of parameters in model ( 1) is:

β(λ) GoldS = argmin β Y -Xβ 2 2 + λ β 1 (2)
where 2 is the quadratic loss, β 1 = p j=1 β j is the Lasso penalty on the parameter size, and λ > 0 controls the amount of regularization. When λ is large enough (which depends on data), all coefficients are forced to be exactly zero. Inversely, λ = 0 corresponds to the unpenalized ordinary leastsquares estimate.

Y -Xβ 2 2 = n i=1 (Y i -X i β)
This model on complete data (no left-censored measures) is considered as a reference when comparing the other methods applied to incomplete datasets that include left-censored values.

The Tobit model

Because of the detection limit, Y i can be left-censored. Let LOD i be the (fixed and known) censoring threshold of subject i. To simplify, we consider LOD i = LOD. Z i is the observed response. The so-called Tobit model [START_REF] Tobin | Estimation of relationships for limited dependent variables[END_REF] can be defined as:

Z i = Y i if Y i > LOD LOD if Y i ≤ LOD (3)
where Y i is the response variable defined in model [START_REF] Paxton | Longitudinal analysis of quantitative virologic measures in human immunodeficiency virus-infected subjects with > or = 400 CD4 lymphocytes: implications for applying measurements to individual patients. National Institute of Allergy and Infectious Diseases AIDS Vaccine Evaluation Group[END_REF]. We can equivalently write:

Z i = max (Y i , LOD) , or Z i = δ i Y i + (1 -δ i ) LOD, (4)
where

δ i = I (Y i >LOD)
is a censoring indicator. The idea behind the Tobit regression model is to deal with the left-censored variable Z as the outcome of a normally distributed latent variable Y.

Simple imputation

Simple imputation is a substitution method that replaces left-censored values with a single value, LOD. LOD/2 is another common choice. Let βLOD be the ordinary least squares estimate of model:

Z i = X i β + ε i , i = 1, • • • , n (5)
Simple imputation is widely used for its simplicity. However, replacing any censored observation by a single value may lead to biased parameter estimates. Beerenwinkel et al. [START_REF] Beerenwinkel | The individualized genetic barrier predicts treatment response in a large cohort of HIV-1 infected patients[END_REF] applied Lasso-regularized linear regression with the naïve approach of replacing the unobserved undetectable value with the limit of detection of the assay. Then the Lasso estimator of parameters in model ( 5) is:

β(λ) LOD = argmin β Z -Xβ 2 2 + λ β 1 . ( 6
)

Maximum likelihood estimation

In the Tobit model ( 1)-(3), one can assume that when Z = LOD, the density function of Z is equal to the probability of observing Y ≤ LOD and for Z > LOD the density function of Z is the same as the density of Y. The likelihood function takes the form:

L β, σ 2 = n i=1 P(Y i |Y i >LOD, X i ) δ i P(Y i |Y i ≤ LOD, X i ) 1-δ i
When the Gaussian distribution for the outcome is assumed, the log-likelihood function can be written as:

lnL β, σ 2 = n i=1 δ i ln f G Y i , X i , β, σ 2 + (1 -δ i ) ln F G LOD, X i , β, σ 2 (7) with f G u, X i , β, σ 2 = e (u-X i β) 2 -2σ 2 √
2πσ 2 the Gaussian probability density function of Y i with mean X i β and constant variance σ 2 evaluated at u and

F G v, X i , β, σ 2 = v -∞ f G u, X i , β, σ 2 du
is the corresponding Gaussian cumulative distribution function evaluated at v. Let βMLE be the maximum likelihood estimation obtained by maximizing [START_REF] Dinse | Accomodating measurements below a limit of detection: A novel application of Cox regression[END_REF]. Extensions to other distributions have also been explored [START_REF] Sigrist | Using the censored gamma distribution for modeling fractional response variables with an application to loss given default[END_REF]. Several works have shown the superiority of this method [START_REF] Wiegand | Comparison of models for analyzing two-group, cross-sectional data with a gaussian outcome subject to a detection limit[END_REF][START_REF] Hewett | A comparison of several methods for analyzing censored data[END_REF][START_REF] Uh | Evaluation of regression methods when immunological measurements are constrained by detection limits[END_REF][START_REF] Kafatos | Regression models for censored serological data[END_REF]. However, when the parametric model is misspecified, the sample size is small or the percent censoring is high, the maximum-likelihood estimation method has been shown to perform poorly [START_REF] Helsel | More than obvious: Better methods for interpreting nondetect data[END_REF].

The Lasso penalty applied to some likelihood function has become an established and relatively standard technique. However, when the likelihood function is a more complex function of the model parameter, such as the likelihood function for the Tobit model ( 7), adding a nondifferentiable penalty leads to a computational challenging optimization.

Quantile regression and censored quantile regression

Quantile regression, particularly least absolute deviations (LAD) regression, has been applied to left-censored data [START_REF] Eilers | Quantile regression for the statistical analysis of immunological data with many non-detects[END_REF]:

βLAD = argmin β 1 n n i=1 |Z i -X i β|
Median regression is a natural alternative to the usual mean regression in the presence of heteroscedasticity or when the normality assumption is violated. Simple imputation using robust regression may be less sensitive to the influence of censored observations. Lasso-regularized least absolute deviations regression has been investigated in the literature (e.g. [START_REF] Belloni | L1-penalized quantile regression in high-dimensional sparse models[END_REF]).

Powell [START_REF] Powell | Least absolute deviations estimation for the censored regression model[END_REF] proposed the LAD estimate specifically for censored data:

βCLAD = argmin β 1 n n i=1 |Z i -max {LOD, X i β}|
Later, the approach was extended to more general quantiles [START_REF] Powell | Censored regression quantiles[END_REF].

The Lasso extension of censored quantile regression (basically, censored LAD) has been analyzed for rightcensored survival data [START_REF] Shows | Sparse estimation and inference for censored median regression[END_REF][START_REF] Wang | Variable selection for censored quantile regression[END_REF][START_REF] Xue | A censored quantile regression approach for the analysis of time to event data[END_REF] and specifically for left-censored data [START_REF] Zhanfeng | A lasso-type approach to variable selection and estimation for censored regression model[END_REF][START_REF] Yue | Bayesian tobit quantile regression model for medical expenditure panel survey data[END_REF][START_REF] Liu | Group variable selection and estimation in the tobit censored response model[END_REF][START_REF] Zhou | LAD-lasso variable selection for doubly censored median regression models[END_REF][START_REF] Alhamzawi | Bayesian elastic net tobit quantile regression[END_REF][START_REF] Müller | Censored linear model in high dimensions[END_REF]. Yu et al. [START_REF] Yue | Bayesian tobit quantile regression model for medical expenditure panel survey data[END_REF] and Alhamzawi et al. [START_REF] Alhamzawi | Bayesian elastic net tobit quantile regression[END_REF] proposed Bayesian approaches using different hyperparameters priors. These methods rely on computationally intensive algorithms. In practice, applications are limited to the n >> p case. Others [START_REF] Zhanfeng | A lasso-type approach to variable selection and estimation for censored regression model[END_REF][START_REF] Liu | Group variable selection and estimation in the tobit censored response model[END_REF][START_REF] Zhou | LAD-lasso variable selection for doubly censored median regression models[END_REF][START_REF] Müller | Censored linear model in high dimensions[END_REF] derived theoretical properties of the Lasso-regularized censored least absolute deviations regression, but the algorithmic development was not a priority in these works and the practical use was limited to n >> p or not addressed. To our knowledge, there are no publicly available software tools that implement the Lasso extension of Powell's approach and no simple implementation relying on existing packages seems straightforward.

Non-parametric Buckley-James

Left-censored outcome data can be analyzed using methods designed for right-censored survival data by reversing the outcome scale. For instance, Gillespie et al. [START_REF] Gillespie | Estimating population distributions when some data are below a limit of detection by using a reverse Kaplan-Meier estimator[END_REF] proposed the reverse Kaplan-Meier and Dinse et al. [START_REF] Dinse | Accomodating measurements below a limit of detection: A novel application of Cox regression[END_REF] reversed the Cox method (though in the case of leftcensored exposures and uncensored outcome). After the Cox model, the accelerated failure time model is the most frequently used regression model for right-censored survival data. It directly links the expected response to predictors, analogously to the classical linear regression approach. A popular method for fitting the accelerated failure time model is the Buckley-James estimator [START_REF] Buckley | Linear regression with censored data[END_REF], an extension of the least squares principle. The idea is to impute the censored values by their estimated conditional mean to provide censoring and predictor values:

Z * i = δ i Y i + (1 -δ i ) E (Y i |Y i ≤ LOD, X i ) (8)
with

E (Y i |Y i ≤ LOD, X i ) = LOD -∞ uf (u, X i , β) du F (LOD, X i , β)
where f (u, X i , β) is the (unknown) probability density function of Y i with mean X i β evaluated at u and F (u, X i , β) is the corresponding cumulative distribution function. By "flipping" the data (turning it from leftcensored to right-censored), the application of algorithms previously developed for right-censoring is direct and has been performed in other contexts [START_REF] Marschner | Clinical trials using HIV-1 RNA-based primary endpoints: Statistical analysis and potential biase[END_REF]. We consider M an arbitrary constant that equals or exceeds the largest observation. Then subtract all uncensored and left-censored outcomes from M. The left-censored at LOD variable Z is then replaced by M -Z which is right-censored at

M -LOD. Let (M -Z i ) * be imputed as δ i (M -Y i ) + (1 -δ i ) E (M -Y i |M -Y i ≥ M -LOD, X i ).
Then, we can calculate the conditional expectation by

E (M -Y i |M -Y i ≥ M -LOD, X i ) = ∞ M-LOD uf (u, X i , β) du 1 -F (M -LOD, X i , β) (9)
where F (u, X i , β) is now the (unknown) cumulative distribution function of M -Y i with mean M -X i β evaluated at u, which can be estimated, for example, by Kaplan-Meier. The Buckley-James estimate, βNonParBJ , can be computed using a semiparametric iterative algorithm that alternates between imputation of censored values according to [START_REF] Eilers | Quantile regression for the statistical analysis of immunological data with many non-detects[END_REF] and least-squares estimation.

The main drawback of this method is that convergence of the algorithm is not guaranteed. Due to the discontinuous nature of the estimating function (formulation (8) makes βNonParBJ to be a piecewise linear function in β), the iterative procedure may oscillate between different parameter values. The problem is of practical importance in situations where the effect of predictors is small or in small samples [START_REF] Wu | A large-scale monte carlo study of the Buckley-James estimator with censored data[END_REF] (which could be worse in highdimensional settings). To circumvent this problem, a onestep algorithm that stops at the first iteration is used in some works [START_REF] Johnson | Variable selection in semiparametric linear regression with censored data[END_REF][START_REF] Wang | The Buckley-James estimator and induced smoothing[END_REF]. This approach is close to a substitution method in which values below the detection limit are replaced by expected values of the missing measurements, provided they are less than the detection limit [START_REF] Gleit | Estimation for small normal data sets with detection limits[END_REF].

Several authors have proposed combining the iterative Buckley-James imputation and methods handling highdimensional predictors: Johnson et al. [START_REF] Johnson | Variable selection in semiparametric linear regression with censored data[END_REF][START_REF] Johnson | On lasso for censored data[END_REF] and Cai et al. [START_REF] Cai | Regularized estimation for the accelerated failure time model[END_REF] used the Lasso, Wang et al. [START_REF] Wang | Doubly penalized Buckley-James method for survival data with high-dimensional covariates[END_REF], used boosting, Wang et al. [START_REF] Wang | Buckley-james boosting for survival analysis with high-dimensional biomarker data[END_REF] used ElasticNet, Johnson et al. [START_REF] Johnson | On path restoration for censored outcomes[END_REF] and Li et al. [START_REF] Zhao | The Dantzig selector for censored linear regression models[END_REF] used the Dantzig selector, and Dirienzo et al. [START_REF] Dirienzo | Parsimonious covariate selection with censored outcomes[END_REF] used parsimonious covariate selection. The Buckley-James estimate can be computed using an iterative algorithm that alternates between imputation of censored values according to [START_REF] Eilers | Quantile regression for the statistical analysis of immunological data with many non-detects[END_REF] and the Lasso:

β(λ) NonParBJ = argmin β (M-Z) * -(M-Xβ) 2 2 + λ β 1 (10)

Gaussian Buckley-James

Alternatively, the Buckley-James imputation (8), assuming the logarithm of HIV viral load follows a Gaussian distribution, can be calculated with the conditional expectation:

E (Y i |Y i ≤ LOD, X i ) = LOD -∞ uf G u, X i , β, σ 2 du F G LOD, X i , β, σ 2 (11)
where f G and F G are the Gaussian density and cumulative distribution functions defined in [START_REF] Dinse | Accomodating measurements below a limit of detection: A novel application of Cox regression[END_REF]. Again, the solution can be computed by iteratively alternating between imputation based on (11) and parameter estimation using ordinary least squares, βGaussBJ , or the Lasso:

β(λ) GaussBJ = argmin β Z * -Xβ 2 2 + λ β 1 (12)

Graphical illustration in a low-dimensional setting

To illustrate the difference between estimation methods we generated data from the simple linear model (p = 1):

Y i = X i β+ε i with i = 1, • • • , n, X ∼ N(0, 1), ε ∼ N(0, σ 2 ).
β was set to 10 and σ 2 was chosen such that the signalto-noise ratio was 4:3. A limit of detection was then fixed to obtain the desired censoring rate: moderate, 20 %, high, 50 % or severe, 70 %. In Fig. 1, predicted regression lines are obtained using different methods: the true model that generated the data, the gold standard (ordinary least squares with uncensored data), maximum likelihood estimation (MLE), which is identical to the Gaussian Buckley-James estimation (BJ) when p = 1, non-parametric Buckley-James (BJ), least absolute deviations (LAD) and censored LAD regressions, simple imputation by the limit of detection (LOD) and by LOD/2.

Notice that simple imputation by LOD and LOD/2 are the most distant regression lines from the true and gold standard lines, in an opposite way: simple imputation by LOD tends to overestimate the response values while simple imputation by LOD/2 tends to underestimate them.

Maximum likelihood estimation shows one of the best behaviors, but the computational complexity dramatically increases with p (results not showed). Mean and median regressions with simple imputation by LOD are quite close and are the closest when the estimation situation is easy (high n, low censoring rate). The censored LAD shows better results than censored mean regression (MLE and Buckley-James) for small sample size while the inverse is observed when n = 100. Gaussian Buckley-James and MLE are identical, but their differences increase when p increases (results not showed). In this i.i.d. generated from a Gaussian distribution example, the Gaussian Buckley-James estimate shows better behavior than nonparametric Buckley-James, the difference being higher when n is small.

Tuning parameter selection

K-fold cross-validation is routinely applied to select the optimal regularization parameter when the main goal of the study is prediction. Data D is randomly chunked into K disjoint blocks of approximately equal size. To avoid a potentially unbalanced partition, we consider stratified K-fold cross-validation, i.e. each fold contains roughly the same proportion of censoring as in the whole sample. D k is the learning data, used to estimate coefficients. D k is the test data, not used in the estimation process and then used to evaluate the loss function L. This K-fold cross-validation can be written as:

CV(λ) = 1 K K k=1 L β(λ) D k , D k ( 13
)
Fig. 1 Observed and predicted simulated censored data. Regression lines are obtained using: the true model that generated the data, the gold standard, maximum likelihood estimation (MLE), non-parametric Buckley-James (BJ), least absolute deviations (LAD) and censored LAD regressions, simple imputation by the limit of detection (LOD) and by LOD/2

CV is evaluated on a grid of λ-values. The highest value, λ max , corresponds to the smallest value of λ for which all coefficients are zero. The lowest value, λ min , corresponds to the unpenalized solution (when feasible). We choose the λ value that minimizes the CV function.

Squared error loss is one of the most widely used loss functions:

L( β(λ) D k , D k ) = 1 n k i∈D k Y i -X i β(λ) D k 2 , (14)
where n k is the sample size of D k . However, Y is a latent variable not fully observed due to the detection limit. This loss function could be used only for the gold standard in (2), with simulated data. Again, the simplest imputation strategy consists in replacing Y with Z, in D k . Alternatively, Buckley-James strategies could replace censored Y i values in the test data D k by their conditional expectation estimated using the learning data [START_REF] Johnson | On lasso for censored data[END_REF].

On the other hand, a loss function differentiating the contribution of uncensored and censored data would be useful. Assuming the Gaussian distribution of the HIV viral load [START_REF] Dinse | Accomodating measurements below a limit of detection: A novel application of Cox regression[END_REF], the following loss function could be derived:

L G β(λ) D k , D k = 1 n unc k i∈D k i uncensored Y i -X i β(λ) D k 2 + 2 σ 2 D k n unc k i∈D k i censored -ln F G LOD, X i , β(λ) D k , σ 2 D k (15)
where n unc k is the number of uncensored observations in D k . The loss function L G in ( 15) is proportional and equivalent to the negative Gaussian log-likelihood loss function, but allows for comparison with the squared loss in [START_REF] Jacqmin-Gadda | Analysis of left-censored longitudinal data with application to viral load in HIV infection[END_REF].

Implementation issues

All statistical analyses, comparisons and implementations were performed using the computing environment R (R Development Core Team, 2017) [START_REF] Core | R: A language and environment for statistical computing[END_REF]. We used the function cv.glmnet from package glmnet [START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF] to choose the optimal λ value of Lasso linear regression on complete data (GoldS) and Lasso linear regression with simple substitution of left-censored values by the detection limit (LOD). We implemented the Lasso nonparametric Buckley-James (NonParBJ) using the bujar package [START_REF] Wang | bujar: Buckley-James regression for survival data with high-dimensional covariates[END_REF]. We modified the function to support stratified K-fold cross-validation and conserve the same proportion of censoring in all the folds. Lasso Gaussian Buckley-James (Gaussian BJ) was implemented in a new function cvGaussBJ. Algorithm 1 specifies how to solve the problem. The stopping criterion is based on the difference between current and previous regression coefficient estimates, variance estimates, and imputed data. Because of the tendency to oscillate between different parameter values of the iterative procedure, the algorithm is also stopped if the number of oscillations is high [START_REF] Wang | Buckley-james boosting for survival analysis with high-dimensional biomarker data[END_REF]. Alternatively, we also considered the one-step algorithm, which stops at the first iteration [START_REF] Johnson | Variable selection in semiparametric linear regression with censored data[END_REF][START_REF] Wang | The Buckley-James estimator and induced smoothing[END_REF]. Crossvalidation based on both imputation and loss function accounting for censored and uncensored contributions is considered. The Lasso estimation step depends on package glmnet.

All these implementations and an artificial example are available at: https://github.com/psBiostat/left-censored-Lasso.

Prediction of HIV viral load from HIV genotypic mutations: real and simulated data

HIV is highly replicative and thus presents high mutation and recombination rates which could lead to the development of HIV drug resistance and consequently reduce the efficacy of antiretroviral treatment. To optimize the control of the evolution of HIV drug resistance, HIV viral load is routinely monitored to identify treatment failure, and HIV genotypic tests are commonly performed before a switch to a new treatment regimen in patients already treated or at the initiation of the first treatment in naive HIV-infected patients [START_REF] Iyidogan | Current perspectives on HIV-1 antiretroviral drug resistance[END_REF].

Algorithm 1 Lasso-regularized Gaussian Buckley-James Initialization:

β(0) ← argmin β n i=1 (Z i -X i β) 2 + λ β 1 σ 2(0) ← 1 n n i=1 Z i -X i β(0) 2
while the stopping criterion is not satisfied do Imputation step

Z * (k) i ← δ i Y i + (1 -δ i ) LOD -∞ uf G (u, X i , β(k-1) , σ 2(k-1) ) du F G (LOD, X i , β(k-1) , σ 2(k-1) )
Lasso estimation step

β(k) ← argmin β n i=1 (Z * (k) i -X i β (k) ) 2 + λ β (k) 1 σ 2(k) ← 1 n n i=1 Z * (k) i -X i β(k) 2

end while

Our objective is to compare methods that handle leftcensoring by conditional imputation with methods that handle left-censoring by imputing a single constant value (that is, the Lasso-regularized linear regression with simple imputation by LOD and LOD/2) to predict of HIV viral load by HIV genotypic mutations. The methods accounting for left-censoring by imputing the estimated conditional mean given censoring and predictor values are the Lasso-regularized Buckley-James least square algorithms (with/without Gaussian assumption, with complete convergence/1-step, using cross-validation based on imputation/loss function accounting for censored and uncensored contributions).

Real data

The database used in this study was provided by the Standardization and Clinical Relevance of HIV Drug Resistance Testing Project for the Forum for Collaborative HIV Research [START_REF] Cozzi-Lepri | Initiatives for developing and comparing genotype interpretation systems: external validation of existing rule-based interpretation systems for abacavir against virological response[END_REF]. Patients included in this study were all treatment-experienced and switched to an abacavircontaining regimen. The investigated drug, abacavir, is a nucleoside reverse transcriptase inhibitor (NRTI) that blocks HIV reverse transcriptase.

The sample size n = 99 was slightly smaller than the number of predictors p = 121. 54 of the 121 predictors correspond to the presence or absence of specific mutations in the reverse transcriptase gene (RTG), which were reported to be probably associated with resistance to abacavir, multi-NRTI, NRTI (other than abacavir) or nonnucleoside reverse transcriptase inhibitors (NNRTI) at the time of the study [START_REF] Shafer | HIV-1 drug resistance mutations: an updated framework for the second decade of HAART[END_REF][START_REF] Johnson | Update of the drug resistance mutations in HIV-1: December 2009[END_REF]. The number of mutations reported to be probably associated with resistance to abacavir or multi-NRTI is low (14%). The other 67 predictors correspond to the presence or absence of specific mutations in the protease gene (PG) reported to be probably associated with resistance to one or several protease inhibitors (PI) at the time of the study [START_REF] Shafer | HIV-1 drug resistance mutations: an updated framework for the second decade of HAART[END_REF][START_REF] Johnson | Update of the drug resistance mutations in HIV-1: December 2009[END_REF]. The number of molecules, including abacavir, ranged from 1 to 6 (with the median number of molecules being 3 and interquartile range 2). In particular, a PI was prescribed in 59% of the patients and 43% received an NNRTI. The response variable is the log-HIV viral load measured at t 8 (8 weeks after treatment initiation at t 0 ). LOD was fixed at 100 copies/mL and the censoring rate was moderate (26%).

Generation of simulated data

HIV viral load appears to have an underlying Gaussian distribution when log-transformed. Therefore, our outcome, Y [START_REF] Wang | Quantile regression in partially linear varying coefficient models[END_REF] i is generated from a Gaussian distribution. We simulated 200 data sets of size n = 100 and p = 100 predictors from the model:

Y (8) i = β 0 + β (0) 1 Y (0) i + X i β + ε i i = 1, • • • , n where • Y (0)
i , the HIV viral load at t 0 is generated by a normal distribution with mean 12 (log 10 copies/mL) and variance 1.

• β (0)

1 represents the change of the slope between the HIV viral load on the day of treatment, t 0 , and 8 weeks later, at t 8 , when no mutations are present and for 1 log 10 /mL higher concentration of viral load at t 0 . We fix β 0 and β (0) 1 to obtain the desired censoring rates: 20% (moderate), 50% (high), and 70% (severe).

• X (n×p) , representing the presence or absence of HIV mutations, is generated by a multinomial distribution with mean 0.15 (the fixed prevalence for all the 100 mutations) and covariance matrix where ij = 0.4 |i-j| (the closer the mutations, the more positively they are correlated).

• β = β 1 , • • • , β p . Among p = 100 candidate mutations only 10% are relevant with effects β j = 1, if j = 1, . . . , 10 and 0 if j > 10. A 1-unit increase in HIV viral load is expected per occurrence of these relevant mutations for a given baseline HIV viral load. • ε is generated from a normal distribution with mean 0 and variance σ 2 chosen such that the signal-to-noise ratio is fixed at 3 : 1.

Our primary goal was to compare competing methods in terms of prediction accuracy. Consequently, we simulated training and test datasets. The former were used to estimate, the latter were used to evaluate the prediction performance. We ensured that training and test datasets contained roughly the same proportions of censoring. We computed the mean squared error on test data as:

MSE = 1 n test n test i=1 Y test i -X test i β( λ)
2 with β( λ) estimated on training data by stratified 5fold cross-validation using a given regression method (Gold standard, Lasso-regularized Non-parametric and Gaussian Buckley-James -with complete convergence and 1-step-and Lasso-regularized linear regression with simple imputation -by LOD and LOD/2-) and the corresponding loss function in the cross-validation criterion. The gold standard uses [START_REF] Jacqmin-Gadda | Analysis of left-censored longitudinal data with application to viral load in HIV infection[END_REF], the others replace the censoring values according to their imputation strategy, and the Gaussian Buckley-James also uses the loss function in [START_REF] Lynn | Maximum likelihood inference for left-censored HIV RNA data[END_REF].

Results

Simulation results

Figure 2 shows the mean prediction error results. "OLS on True Model" corresponds to ordinary least squares for the linear model with uncensored data and only relevant predictors (the so-called oracle estimator). Gold standard (GoldS) corresponds to the Lasso estimation for the linear The imputation by LOD/2 led to poorer results than the imputation by LOD. Thus, for the simple imputation, only LOD imputation (LOD) results are shown. For the Gaussian Buckley-James algorithms (Gaussian BJ), the error is calculated by using both cross-validation with imputation and cross-validation with the loss function indicated in [START_REF] Lynn | Maximum likelihood inference for left-censored HIV RNA data[END_REF], but only the best results are shown.

The Gaussian Buckley-James method presented an oscillating behavior in 9.5% of the generated samples when the convergence rate was 20%. This percentage rose to 82.5% and 95.0% when the convergence rates were 50% and 70%, respectively. For the Gaussian Buckley-James using the 1-step algorithm (Gaussian BJ 1-Step), results using the two cross-validation approaches were almost identical. Nevertheless, for the Gaussian Buckley-James with complete convergence, a notable improvement was obtained when applying [START_REF] Lynn | Maximum likelihood inference for left-censored HIV RNA data[END_REF].

The higher the rate of censoring, the less information is available to train the models and, unsurprisingly, the higher is the prediction error. For a moderate rate of censoring (20%), all methods show a good performance close to that of the gold standard GoldS. When the rate of censoring is 50%, Gaussian BJ shows the lowest prediction error, followed by Gaussian BJ 1-step, NonParBJ and finally simple imputation, which shows more errors.

The same patterns but more pronounced were observed with a severe rate of censoring (70%). Taking the knowledge about the distribution into account appears to have only a slight impact. Simple imputation yields the poorest results. In addition, it showed high variability with some extreme errors.

Application to real data

The Lasso-regularized Buckley-James least square algorithm that showed the best behavior in the simulation study (with complete convergence and cross-validation based on the loss function L G in [START_REF] Lynn | Maximum likelihood inference for left-censored HIV RNA data[END_REF]) was applied to real data, as well as the Lasso-regularized non-parametric Buckley-James method and simple imputation (by LOD and LOD/2).

Regularization parameters were estimated by stratified cross-validation in order to ensure that each fold had the same proportion of censoring as in the corresponding data set (26%). In addition, because some mutations were relatively infrequent, we used 20-fold cross-validation. Indeed, the higher the number of folds, the lower the probability of randomly obtaining test sets with no subject exposed to infrequent HIV drug resistance mutations.

Figure 3 shows two examples of the observed HIV viral load at t 0 and the observed and estimated HIV viral load at t 8 . As in the low-dimensional case shown in Fig. 1, simple imputation by LOD predicted the highest values of HIV viral load. Inversely, simple imputation by LOD/2 estimates the lowest values of HIV viral load. The difference between the two estimates at 8 weeks is >0.5 log 10 copies/mL, which is clinically relevant. Lasso-regularized Buckley-James least square algorithms (with/without Gaussian assumption), often gave a prediction in between. This tendency was increased when the censoring rate was high or when the sensitivity of the assay was low.

When applying the Gaussian Buckley-James method to real data, no oscillating behavior was observed.

Table 1 indicates the number of HIV genotypic mutations selected from the list of mutations that may contribute to a reduced virologic response known at the time of the study [START_REF] Shafer | HIV-1 drug resistance mutations: an updated framework for the second decade of HAART[END_REF][START_REF] Johnson | Update of the drug resistance mutations in HIV-1: December 2009[END_REF], according to each method applied. The Lasso-regularized Gaussian Buckley-James selected several HIV genotypic mutations suspected of being associated with abacavir or multi-NRTI resistance. Furthermore, it selected a high number of HIV genotypic mutations probably associated with PI resistance and a few probably associated with NNRTI resistance. This selection of a large number of candidate predictors seems to be relevant because all patients received an abacavir-containing regimen and a high percentage of patients received regimens including a PI-and/or NNRTI. The Lasso-regularized non-parametric Buckley-James selected fewer mutations, and especially fewer mutations in PG, probably due to PI resistance. Simple imputation of the LOD or LOD/2 selected few mutations. In particular, only 1 of the 5 mutations in RTG probably associated with abacavir resistance was retained.

Discussion

Simple imputation of the detection limit or of half of this limit is an ad hoc approach to address left-censored outcome data. However, in standard (low-dimensional) settings, it leads to biased estimates of parameters and standard errors. In our high-dimensional simulation study, simple imputation using Lasso-regularized leastsquares showed poor performance. As in low-dimensional Fig. 3 Two examples of individual predicted values using Lasso regularization techniques (Gaussian Buckley-James, non-parametric Buckley-James, simple imputation by LOD and LOD/2) and observed values. One example is censored (at left) and the other is uncensored (at right) settings, approaches accounting for left-censoring outperformed simple imputation.

In this work, we propose a Lasso-regularized Gaussian Buckley-James algorithm, according to the usual Gaussian assumption of log-transformed HIV viral load. Because of the well-known convergence problems of the iterative Buckley-James procedure, we implemented two algorithms, the first algorithm running until convergence and the second one being stopped after one step [START_REF] Johnson | Variable selection in semiparametric linear regression with censored data[END_REF][START_REF] Wang | The Buckley-James estimator and induced smoothing[END_REF]. This one-step algorithm showed similar results in the simulation study. Other solutions have been proposed to deal with convergence problems in low-dimensional settings [START_REF] Ueki | A note on automatic variable selection using smooth-threshold estimating equations[END_REF][START_REF] Wang | The Buckley-James estimator and induced smoothing[END_REF] and could be investigated in future research.

As in other works [START_REF] Johnson | On lasso for censored data[END_REF], we implemented a crossvalidation criterion for the tuning parameter based on imputing values to Y i in the test set from conditional expectations estimated using the learning set. We also proposed a cross-validation criterion based on a loss function that accounts for the different contribution of censored and uncensored values. Almost identical results were obtained when applying the two cross-validation criteria to the one-step algorithm. However, when running the algorithm until convergence, better results were obtained with the cross-validation criterion based on a loss function that accounts for censored and uncensored contributions.

On the other hand, we reversed the Lasso-regularized non-parametric Buckley-James method previously applied to right-censored survival data [START_REF] Johnson | Variable selection in semiparametric linear regression with censored data[END_REF][START_REF] Cai | Regularized estimation for the accelerated failure time model[END_REF][START_REF] Wang | Buckley-james boosting for survival analysis with high-dimensional biomarker data[END_REF][START_REF] Johnson | On lasso for censored data[END_REF] in order to apply to left-censoring due to detection limits. Foreseeably, in our homoscedastic Gaussian outcome data scenario, the Gaussian Buckley-James showed better behavior than the non-parametric algorithm. However, accounting for the knowledge about the distribution seems to have had a slight influence. When the Gaussian assumption is violated, non-parametric imputation using Kaplan Meier is perhaps the best option.

We provide a publicly available R code to compute the methods introduced in this work (https://github. com/psBiostat/left-censored-Lasso). It would be interesting to compare the Lasso-regularized Buckley-James least squares method to Lasso-regularized censored LAD method. The Lasso extension of censored LAD has been proposed in different works [START_REF] Shows | Sparse estimation and inference for censored median regression[END_REF][START_REF] Wang | Variable selection for censored quantile regression[END_REF][START_REF] Xue | A censored quantile regression approach for the analysis of time to event data[END_REF][START_REF] Zhanfeng | A lasso-type approach to variable selection and estimation for censored regression model[END_REF][START_REF] Yue | Bayesian tobit quantile regression model for medical expenditure panel survey data[END_REF][START_REF] Liu | Group variable selection and estimation in the tobit censored response model[END_REF][START_REF] Zhou | LAD-lasso variable selection for doubly censored median regression models[END_REF][START_REF] Alhamzawi | Bayesian elastic net tobit quantile regression[END_REF][START_REF] Müller | Censored linear model in high dimensions[END_REF]. However, to our knowledge, there is no publicly available implementation, and no simple implementation relying on existing packages seem straightforward. Moreover, several works have shown the superiority of maximum likelihood estimation in low-dimensional settings when the Gaussian assumption is valid [START_REF] Wiegand | Comparison of models for analyzing two-group, cross-sectional data with a gaussian outcome subject to a detection limit[END_REF][START_REF] Hewett | A comparison of several methods for analyzing censored data[END_REF][START_REF] Uh | Evaluation of regression methods when immunological measurements are constrained by detection limits[END_REF][START_REF] Kafatos | Regression models for censored serological data[END_REF]. Nevertheless, optimization strategies for complex likelihood functions (such as that in Eq. 7) including penalties that are not smooth are not obvious.

To illustrate the application of the methods on real data, we consider a data set from the Standardization and Clinical Relevance of HIV Drug Resistance Testing Project for the Forum for Collaborative HIV Research. The data set used to illustrate the initial data set is characterized by a sample size-to-predictors ratio of around 1. There is no gold standard to measure and compare predictive performance of the different methods when using censored outcome data. All patients were being treated with abacavir, an NRTI, so we expected our methods to select a high number of HIV genotypic mutations known to contribute to abacavir and NRTI resistance. Furthermore, a high number of patients were on PI-and/or NNRTI-containing regimens, and a selection of several HIV genotypic mutations reported to be probably associated with resistance to any of these molecules was also expected [START_REF] Shafer | HIV-1 drug resistance mutations: an updated framework for the second decade of HAART[END_REF][START_REF] Johnson | Update of the drug resistance mutations in HIV-1: December 2009[END_REF]. In that sense, the Gaussian and non-parametric Buckley-James methods showed more coherent results with the literature compared to simple imputation.

Otherwise, the data presented a moderate censoring rate of 26%, which is a realistic magnitude [START_REF] Wiegand | Comparison of models for analyzing two-group, cross-sectional data with a gaussian outcome subject to a detection limit[END_REF][START_REF] Marks | The spectrum of engagement in HIV care: do more than 19% of HIV-infected persons in the US have undetectable viral load?[END_REF] in studies measuring HIV viral load. However, high or even severe censoring rates were found in older studies with a high limit of detection (LOD) or particular populations with a low treatment failure rate [START_REF] Wiegand | Comparison of models for analyzing two-group, cross-sectional data with a gaussian outcome subject to a detection limit[END_REF][START_REF] Dao | Study to Understand the Natural History of HIV and AIDS in the Era of Effective Therapy (SUN) Investigators. Low vitamin D among HIV-infected adults: prevalence of and risk factors for low vitamin D levels in a cohort of HIV-infected adults and comparison to prevalence among adults in the US general population[END_REF][START_REF] Leon | Rate and predictors of progression in elite and viremic HIV-1 controllers[END_REF]. Furthermore, left-censoring due to the lower detection limit of an assay is a problem in many fields such as biology, immunology, chemistry, and the environmental sciences in which high censoring rates may be frequent. Our simulation study shows that the difference in performance between Lasso-regularized Buckley-James methods and Lasso-regularized simple imputation methods increased with the censoring rate.

In our simulations and real application, the detection threshold was the same for all subjects. The detection threshold may vary among subjects, for example, in multicentric studies. Our R code also supports multiple lower limits of quantification. However, the findings should be interpreted with caution: differences in technological equipment could be a confounding factor that might help explain the differences in patient response to HIV treatment (in addition to HIV mutations). Adjusting or stratifying for the hospital would then be necessary.

In this study we focused on the prediction performance of Lasso-regularized methods. In clinical applications, even when prediction accuracy is the main objective, researchers aim to identify which predictors are more strongly associated with outcome. Our proposal could be easily extended or adapted to support other Lassotype penalties. When the primary goal is to infer the set of truly relevant variables, the adaptive Lasso and the bootstrap-enhanced Lasso could thus be considered.

Fig. 2

 2 Fig.2Mean Square Error (MSE) calculated on test data of size n test = 100 and 200 replications using gold standard method (GoldS), Gaussian Buckley-James (Gaussian BJ) and 1-Step version (Gaussian BJ 1-Step), non-Parametric Buckley-James (NonParBJ) and simple imputation by LOD (LOD). Data were generated to present different censoring rates: 0% (uncensored), 20% (moderate), 50% (high), and 70% (severe), which are represented by different gray degrees

Table 1

 1 Distribution of 121 HIV genotypic mutations included in real data study according to knowledge at study time and reported in[START_REF] Shafer | HIV-1 drug resistance mutations: an updated framework for the second decade of HAART[END_REF][START_REF] Johnson | Update of the drug resistance mutations in HIV-1: December 2009[END_REF] and number of HIV genotypic mutations selected by Lasso regularized methods

	Number of HIV genotypic mutations	Gaussian BJ	NonPar BJ	LOD	LOD/2
	present in real data study
	5 in RTG probably associated	3 (60%)	3 (60%)	1 (20%)	1 (20%)
	with abacavir resistance
	13 in RTG probably associated	7 (54%)	4 (31%)	4 (31%)	4 (31%)
	with multi-NRTI resistance
	6 in RTG probably associated	3 (50%)	2 (33%)	1 (17%)	1 (17%)
	with NRTI resistance (other than abacavir)
	30 in RTG probably associated	12 (40%)	11 (37%)	8 (27%)	7 (23%)
	with NNRTI resistance
	67 in PG probably associated	40 (60%)	22 (33%)	15 (22%)	14 (21%)
	with PI resistance
	121 Total	65 (54%)	42 (35%)	29 (24%)	27 (22%)

Acknowledgements

We acknowledge members of the Standardization and Clinical Relevance of HIV Drug Resistance Testing Project for the Forum for Collaborative HIV Research. We thank "Sidaction, Ensemble contre le Sida", France, for their continuous support. We would like to thank Binbin Xu, postdoctoral researcher at SISTM research team from Inserm BPH U1219 & Inria BSO, for his help on testing the R code.

Funding

This work was partially supported by the "Investissements d'Avenir" program managed by the ANR under reference ANR-10-LABX-77

Availability of data and materials R code and artificial example are available at https://github.com/psBiostat/ left-censored-Lasso.

Authors' contributions

PS developed the algorithms and corresponding R code, carried out the statistical analysis and helped to draft the manuscript. MA developed the algorithms, revised the R code and drafted the manuscript. RT designed and supervised the applied research. RT and LW interpreted the results of the analysis. DC revised the methodology. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.