S Daggumati

I De Baere

W Van Paepegem

J Degrieck

J Xu

S V Lomov

I Verpoest

Local damage in a 5 -harness satin weave composite under static tension: part i -experimental analysis

Keywords: Textile composite, Transverse cracking, Acoustic emission, Multiscale modeling, Weft yarn damage

Local damage in a 5 -harness satin weave composite under static tension: part i -experimental analysis

Introduction

Damage accumulation in textile composites is a complicated process, and development begins on the micro scale with the fibre matrix debonding, the matrix cracking and the fibre failure (the micro scale defines the arrangement of fibres in an impregnated yarn or fibrous ply). On the meso scale, damage develops by intra-yarn cracking and delaminations (the meso scale defines the internal structure of the reinforcement, variation of the fibre direction and the fibre volume fraction inside the yarns and the fibrous plies). Finally, the macro failure of the composite is characterized by dense cracking, intersection of several small cracks (crack conjunction) and fibre rupture [START_REF] Ivanov | Failure Analysis of Triaxial Braided Composite[END_REF] (the macro scale defines the three-dimensional geometry of the composite part and the distribution of local reinforcement properties).

Starting from the micro scale, Figure 1 schematically depicts the weft yarn damage that is affected by the restriction of local ply deformation caused by the surrounding layers. The influence of local constraints is manifested by a change in the weft yarn crack location (edge / centre) based on the location of the chosen ply in the laminate (inner / surface).

In order to obtain comprehensive knowledge of the damage phenomena occurring in textile composites, it is necessary to understand the factors that contribute to the damage at different length scales. Apart from the natural variability in yarn spacing and its dimensions, the unit cells of textile reinforcement are theoretically the same. However, parameters such as the fibre orientation, localized fibre spacing and packing often exhibit a wide statistical variation when evaluated on the micro scale in a processed composite.

Therefore, some localized micro volumes are stressed more than others. The stress inhomogeneity is further enhanced by the inhomogeneity of the elastic properties of the composite constituents. The inhomogeneity of the stress field, coupled with the inhomogeneity of the strength properties of the reinforcing elements, the matrix and the interface, lead to the gradual damage development in composites [START_REF] Dzenis | Analysis of microdamage evolution histories in composites[END_REF]. Moreover, the micro scale stress-strain state defines the number of transverse cracks in a particular yarn and how they are placed over its entire cross-section [START_REF] Koissin | Fibre distribution inside yarns of textile composite: geometrical and FE modelling[END_REF].

In addition to the micro scale parameters mentioned above, at the ply level (meso scale), geometrical parameters such as the yarn crimp and the variations in intra-yarn volume fraction can contribute to the stochastic nature of the stress concentration in the micro volumes [START_REF] Kurashiki | A numerical simulation of damage development for laminated woven composites ECCM-13 Proceedings[END_REF][START_REF] John | Longitudinal and transverse damage taxonomy in woven composite components[END_REF]. Finally on the scale of laminate, in the textile composite, different layers have different yarn nesting patterns. The nesting conditions of the particular laminate are defined by the (random) shifting of its layers during the manufacturing process. Hence, the stress-strain state inside the ply (as in the micro volumes) depends on the placement of ply inside the laminate [START_REF] Anzelotti | Mesomechanic strain analysis of twillweave composite lamina under unidirectional in-plane tension[END_REF][START_REF] Ivanov | Damage analysis in textile composites[END_REF][START_REF] Ivanov | Stress distribution in outer and inner plies of textile laminates and novel boundary conditions for unit cell analysis[END_REF]. The above mentioned parameters of the textile composites at various hierarchical levels can lead to different damage initiation conditions, different patterns of progressive damage in the chosen ply of a laminate. In general, micro scale damage initiation in a textile composite is not only affected by the stochastic micro level geometrical variables, but also by the meso and laminate level geometrical parameters.

A combination of different experimental techniques is used in the current study to understand the detailed damage phenomena and the damage history. Damage initiation is detected using the acoustic emission (AE) technique during the tensile test. AE provides the damage initiation threshold and the critical stresses of the composite. In order to detect the location of the damage, microscopic analysis is used. The above mentioned methodology has been proposed by many researchers [START_REF] Mattsson | Damage in NCF composites under tension: Effect of layer stacking sequence[END_REF][START_REF] Morscher | Modal acoustic emission of damage accumulation in a woven SiC/SiC composite[END_REF][START_REF] Giordano | An acoustic-emission characterization of the failure modes in polymer-composite materials[END_REF] to study the damage initiation and propagation in textile composites. From the microscopic images captured during the tensile test, the major damage observed in the satin weave composite is weft yarn cracking before the final failure of the laminate. Based on the microscopic images, apart from the detection of damage initiation stress and locations, a novel procedure is introduced to interpret the layer-wise sequential damage at different locations in the laminate. The aim of using this technique is to provide the information regarding the history of the weft yarn transverse damage in different layers of the laminate. This then also provides information on the influence of adjacent layers (internal yarn nesting) on the meso damage behavior of a local ply. Finally, the correlation between AE event data and the microscopic crack density provided confidence in the interpretation of the complex damage behavior.

The mechanical testing procedure and the material used

Tensile tests on the composite specimens were performed using a standard Instron machine (Instron 4505, test speed 1 mm/min). Composite samples used for the tensile tests were prepared according to the ASTM standard (width 25mm, gauge length 170mm, thickness 2.5 mm), with the end tabs of the same material (Carbon-PPS). The local strain during the tension test is measured using the digital image correlation technique (LIMESS). At the end of the loading process, the average strain is computed across the entire window of the speckle pattern ( mm 37 23×

) used for the full-field strain registration [START_REF] Ivanov | Strain mapping analysis of textile composites[END_REF]. Along with the LIMESS, the AE measurement and data acquisition system (VALLEN) is used to detect the damage initiation stress, followed by the microscopic analysis of the composite specimens for the inspection of damage locations.

Acoustic emission

It is important to determine the damage initiation and the critical strains over the stressstrain curve obtained from the tension test on the textile composite specimen. Based on the damage initiation and critical strains observed, the composite fatigue and durability life can be determined [START_REF] Lomov | A comparative study of tensile properties of non-crimp 3d orthogonal weave and multi-layer plain weave e-glass composites[END_REF][START_REF] Ivanov | A comparative study of tensile properties of non-crimp 3D orthogonal weave and multi-layer plain weave E-glass composites[END_REF]. To this end, AE has proven to be a valuable technique for online damage detection as well as for tracking the damage progression. In general, the characteristic failure mechanisms in fibre reinforced composites are initiated at the micro level and result in a spontaneous release of elastic strain energy, which is dissipated as a wave that propagates from the failure source through the medium [START_REF] Giordano | An acoustic-emission characterization of the failure modes in polymer-composite materials[END_REF]. These emitted acoustic events are detected by one of the sensors attached to the composite specimen.

Once the composite test specimen is mounted on the tensile machine, two sensors are attached at the ends of the specimen at a distance of 150mm from each other. Special vacuum grease is used to mount the sensors onto the composite, which will transfer the acoustic events from the composite material to the AE sensors. The complete AE VALLEN system used in the current experiment can be obtained from the VALLEN data manual [START_REF] Vallen | AE Testing (AT) Fundamentals -Equipment -Data Analysis[END_REF]. Details of the AE equipment used for the event registration are shown in Table 1. At the end of the tensile test, AE signals occurred outside the area between sensors, and the external noise was removed from the entire AE data. The obtained AE output combined with the strain data, results are plotted as a logarithm of cumulative acoustic energy versus tensile strain (Figure 3a). By analyzing the different strain levels in Figure 3a, the critical stress at which damage occurs can be detected [START_REF] Ivanov | A comparative study of tensile properties of non-crimp 3D orthogonal weave and multi-layer plain weave E-glass composites[END_REF].

Microscopic analysis

The AE technique described in the above section provides the damage initiation stress quantitatively, but it can not give any indication of the damage locations. In order to visualize the damage locations, microscopic analysis of the composite specimen is a possible technique [START_REF] Searles | An image analysis technique for evaluating internal damage in graphite-fabric/polyimide composites[END_REF][START_REF] Gao | Damage accumulation in woven-fabric CFRP laminates under tensile loading: Part 1. Observations of damage accumulation[END_REF]. The technique presented in this paper (in later sections) is based on capturing and performing the qualitative as well as the quantitative analysis of the microscopic images. Quasi -static tensile tests were performed on the polished composite specimen to allow for the microscopic inspection of damage initiation and propagation. In order to visualize the damage, an OLYMPUS AX-70 microscope is used along with the 'cell-D' (www.microscopy.olympus.eu) software for the analysis of the microscopic images that were captured.

Material

The material under study is a thermo-plastic 5-harness satin weave composite (CD0286 supplied by 'Ten Cate'), which has T300 JB carbon fibres as a reinforcement and PPS as a matrix. The composite laminate is manufactured using the press forming technique, which has eight layers of satin weave fabric with 50:50 weight ratios in the warp and weft directions. The details of the 5-harness satin weave fabric are listed in Table 2.

Due to the sensitivity of the micro scale geometrical parameters on the damage behaviour, there is a special focus on investigating the variation of the local yarn crosssection (which directly controls the fibre packing fraction) at various locations of the laminate (Fig 2a -2d). The yarn dimensions, which were measured in different locations of the laminate, revealed that the 5-harness satin weave weft yarn maintained its elliptical shape, with a major diameter of 1.32 17 . 0 ± mm, and the minor diameter varying at around 02 . 0 15 . 0 ± mm. The above measured yarn dimensions with 3000 filaments in a bundle showed a constant packing fraction that is independent of the location of the chosen weft yarn.

Results and discussion

Stress -strain diagram, analysis of AE data

In order to obtain critical strains for the damage analysis, the stress versus strain curve is plotted together with the strain versus acoustic energy (Figure 3a). To show consistency in the event pattern for different samples, the event distribution (in a linear scale) is plotted with the tensile stress (Figure 3b). These curves were plotted using the acoustic events obtained from the AE VALLEN, the stress data obtained from the INSTRON loadcell and the LIMESS strain data.

The first phase of the AE data (Figure 3a) is associated with the occurrence of low energy events starting at the global tensile strain around 0.14% in the warp yarn direction (corresponding to approximately 80 MPa stress), which is designated as the damage initiation strain min ε [START_REF] Lomov | A comparative study of tensile properties of non-crimp 3d orthogonal weave and multi-layer plain weave e-glass composites[END_REF]. In the second phase ( min ε to 1 ε -0.14 to 0.3%), the rise in the slope of the AE curve is very small, and this phase is related to the occurrence of sporadic cracks which were associated with a low energy release. In the next phase of the AE Moreover, based on the literature on thermo-plastic composites [START_REF] Parlevliet | Residual stresses in thermoplastic composites--A study of the literature--Part I: Formation of residual stresses[END_REF][START_REF] Parlevliet | Residual stresses in thermoplastic composites -a study of the literature. Part III: Effects of thermal residual stresses[END_REF][START_REF] Parlevliet | Residual stresses in thermoplastic composites--A study of the literature--Part II: Experimental techniques[END_REF], it is observed that the amount of thermal residual stresses accumulated during the manufacturing process can not be neglected. For the carbon-PPS composite plates, thermal stresses were accumulated due to the shrinkage of PPS during the cooling process (micro-mechanical residual stress). In addition, at the micro-level, fibre waviness can influence the residual stresses and hence the failure behavior [START_REF] Parlevliet | Residual stresses in thermoplastic composites--A study of the literature--Part II: Experimental techniques[END_REF]. At the macro level, unsymmetrical weaving of the 5-harness satin weave fabric causes the warping of composite plates. In general, due to the differences in processing conditions such as layer stacking (internal yarn nesting), different composite plates produced with the same materials tend to exhibit differences in the local warping, which is manifested by the partial relieving of thermal stresses.

In the context of analyzing the variations in thermal stresses in different composite plates, and their influence on the critical strains ( min ε , 1 ε , 2 ε ), further investigations were carried out by testing the other batch of carbon-PPS samples made from a different plate. AE tests on the other batch of carbon-PPS samples predicted the damage initiation at an early strain of 0.09% [40 MPa] (Figure 3c). From the AE tests on batch 1 and 2 (Figure 3b,c), it is evident that there is almost a 50% difference in the damage initiation stress due to the variation in thermal residual stresses. Figure 3 (b &c) provides evidence that thermal stresses cause the early damage events with low energy release, which is termed as the damage caused by manufacture induced inhomogeneities. However, in both the cases that are illustrated in Figure 3 (b &c), the rise in event count and the energy release is strongly increased approximately around 200 MPa stress or 0.3-0.4% of average strain. This is a typical value for the intra-yarn damage onset observed experimentally for many textile composites [START_REF] Ivanov | Stress distribution in outer and inner plies of textile laminates and novel boundary conditions for unit cell analysis[END_REF]. Table 3 provides the summary of mechanical properties for the 5harness satin weave carbon-PPS composite under tension.

Microscopic analysis of the damage

Microscopic analysis of the composite specimens is a possible technique for evaluating the damage initiation and propagation in the textile laminate. Along with the AE, microscopic analysis allows for the study of the comprehensive damage behaviour of the composite [START_REF] John | Longitudinal and transverse damage taxonomy in woven composite components[END_REF][START_REF] Lomov | Experimental methodology of study of damage initiation and development in textile composites in uniaxial tensile test[END_REF]. Quasi-static tensile tests were performed on the polished specimens of the composite laminate, which allowed for the microscopic inspection of the damage initiation and propagation along the thickness direction as well as on the traction free surface of the composite. From the captured microscopic images, apart from the traditional crack counting (crack density) on the polished edges, a new technique is introduced to study the effect of layer stacking on the local damage behavior of the weft yarn. With the help of a satin weave laminate cross section (approximately 8mm) (Figure 4a), details of the damage phenomena starting from the crack initiation to the final failure of the composite is explained in this section.

Due to the sensitivity of the initial damage, up to 100 MPa stress, the load on the INSTRON machine is increased in two steps. After the initial damage, the load is increased with a step size of 100 MPa until the final failure. Damage developed in the 5harness satin weave under tensile load takes the form of weft yarn cracking, which is perpendicular to the loading direction (Figure 4a). Moreover, damage initiation occurs at the yarn crimp location and then propagates into the matrix. This type of weft yarn cracking began early in the loading process and almost all weft yarns possess this kind of damage before the final failure of the composite.

To study the effect of local constraints that are posed by the surrounding composite constituents on the weft yarn damage, detailed sequential damage analysis was conducted on the composite specimen. The above mentioned four locations in the laminate (Figure 2a to d) were chosen for this purpose, which can represent different cases where the weft yarn local warping (bending) is allowed or constrained by the adjacent layers [START_REF] Ito | Elastic moduli and stress field of plain-weave composites under tensile loading[END_REF]. Also, the above mentioned four locations represent the various local fabric geometries that are available in the laminate. By constantly monitoring the weft yarn damage at the above mentioned locations during the loading process, the effects of local yarn constraints on the damage initiation stress and location can be realized.

During the loading process of the composite, initial damage is detected at around the 100MPa stress level on the polished edges of the composite. At these low stress levels, damage occurs at the edges of the weft yarn inside the laminate, in the region where there is no contact between the load carrying warp yarn and the perpendicular weft yarn (Figure 4a). By increasing the stress level from 100 to 400 MPa in three steps, the damage is observed at the nested yarn configurations at various locations in the laminate (Figure 4b). Moreover, the number of cracks between these stress levels can be compared to the event count obtained from AE. By increasing the stress from 400 to 500 MPa, the next new damage location is observed on the surface of the laminate (Figure 4c).

Moreover, 80% of the microscopic images on the surface (top and bottom) of the laminate show the occurrence of intra-yarn damage at the centre of the weft yarn. At this stress level (400-500 MPa), along with the damage on the surface layers, the crack density in the nested yarns inside the laminate is increased at various locations. Finally, between the stress levels of 500-680MPa, the next damage is observed at the edges of the weft yarn, where it is tightly packed by load carrying warp yarns (Figure 4d). Any further loading causes the catastrophic failure of the composite specimen. From the above discussed sequential damage analysis, it can be observed that the local fabric geometry whose local bending is restricted tends to show the initial damage. In contrast, the yarn crimps on the surface layers are relatively more free to deform in the out-of-plane direction, which causes the damage on the surface weft yarns at later stages of the loading.

Besides the qualitative analysis, microscopic images were used for the quantitative analysis as well. The accumulation of transverse matrix cracking in the thickness direction is quantified on the polished edges of the composite. The crack density is measured as a function of the applied stress for the layers of the laminate, by dividing the number of cracks measured in each layer with the laminate gauge length [START_REF] Gao | Damage accumulation in woven-fabric CFRP laminates under tensile loading: Part 1. Observations of damage accumulation[END_REF]. Figure 5 shows the crack density as a function of applied tensile stress for three composite samples. The microscopic damage analysis described above is in correlation with the AE event pattern. At around the 250 MPa stress level, the event count obtained from AE as well as the crack density tends to increase.

Conclusion

A comprehensive damage analysis is conducted on the carbon-PPS 5-harness satin weave composite under uni-axial static tensile load. The damage initiation strain is detected using acoustic emission technique, while the damage locations were detected using microscopic analysis on the polished edges of composite specimens in an interrupted tensile testing process. The analysis of the experimental data leads to the following conclusions.

• Damage initiation in the carbon-PPS satin weave composite varies from 40 to 100MPa stress. The variation of the experimental damage initiation stress in different composite plates can be attributed to the variation in micro level (fibre spacing, intra yarn volume fraction etc) as well as the laminate level (nesting) parameters combined with the thermal residual stresses in the composite.

• From the research output of Le Page et al. [START_REF] Page | Finite element simulation of woven fabric composites[END_REF] and the AE event data obtained from the satin weave composite under tensile load, it can be concluded that the energy release rate and the associated local bending of the satin weave weft yarn increases approximately around 200 MPa global tensile stress.

• From the extensive microscopic analysis of damage in several satin weave laminates, the damage initiation in different layers of the laminate depends on the position of the ply in the laminate. The earliest damage occurs in the inner layers, followed by the transverse crack on the surface weft yarns. This behaviour is related to the constraints posed by the neighbouring layers on the local ply.

• Initial transverse cracks tend to be located near the edges of the weft yarns crosssection rather than the centre.

• The above conclusions are entirely supported by the work of Ivanov. [START_REF] Ivanov | Damage analysis in textile composites[END_REF], in which the author reported that the meso stress distribution, and crack density is different for inner and outer layers of the laminate. Moreover, meso stress distribution depends upon the stacking sequence, ply shift and the number of plies in the laminate [START_REF] Ivanov | Stress distribution in outer and inner plies of textile laminates and novel boundary conditions for unit cell analysis[END_REF].

The damage phenomena discussed above focus mainly on the qualitative damage behaviour of the satin weave composite. In part II of this series, deterministic numerical simulations provide overall insight into the mechanics of damage and the local stress behaviour of the yarns. This cannot be determined using the conventional experimental techniques.

curve ( 1 ε to 2 ε

 12 -0.3 to 0.4%), the slope increases rapidly with the occurrence of low and middle energy acoustic events. The final part of the AE curve (after 2 ε ) is associated with the continuous events with the high energy release, indicating the critical crack propagation and eventually catastrophic failure of the composite. Once the stress level reached 500 MPa, the AE sensors were removed from the composite specimen to avoid damage to the sensors due to the catastrophic failure of the composite laminate.

Figure 1 .

 1 Figure 1. Schematic representation of the damage at various levels in a satin weave carbon-pps thermoplastic composite.

Figure 3 .

 3 Figure 3. a) Stress vs. strain and strain vs. AE registration curve (test1); b) AE event pattern (batch 1); c) AE event pattern (batch 2).

  location 3-1 (Surface yarn) top view Damage location 3-1(Surface yarn) side view d) Damage at location 4-1

Figure 4 .

 4 Figure 4. Microscopic sequential damage analysis : a) laminate at100MPa tensile stress; b) tensile stress between 100 -400 MPa; c) tensile stress between 400 -500 MPa; d) tensile stress between 500 -680 MPa (Size bar is 1mm).

Figure 5 .

 5 Figure 5.Crack density as a function of applied tensile stress.

Table 1 . Parameters of the AE (Acoustic Emission) equipment.

 1

	Software	Vallen AMSY-5
	Amplifiers	Vallen AEP4
	Amplification, dB	34
	Discrimination time, ms	0.4
	Rearm time, ms	3.2
	Range, MHz	0.025..1.6
	Sample rate, MHz	5
	Sensors	Digital wave 1025
	Sensor diameter, cm	0.93
	Threshold, dB	35

Table 2 . Satin weave carbon/PPS composite information.

 2

	Carbon fibre type	T300JB
	TEX, (g/km)	198
	End/Pick count, yarns per 10 cm	70
	Yarn filament count	3000
	Filament diameter, (mm)	0.007
	Carbon fibre density, (g/cm 3 )	1.75
	Fabric areal density, (g/m 2 )	285
	Overall fibre volume fraction	50 ± 3%

Table 3 . Experimental mechanical properties of the carbon-PPS satin weave composite.

 3

	Warp direction	Weft direction

ACKNOWLEDGEMENTS

The authors would like to acknowledge the FWO -Vlaanderen for the financial support provided through the project G.0233.06H. Assistance from the laboratory staff at the department of MTM -Kris van der Staey and Johan Vanhuist is gratefully acknowledged. The authors also express their gratitude to 'Ten Cate' for supplying the composite plates.

Weft yarn damage at the laminate level

Ply level-Weft yarn damage on the surface layer (at the weft yarn centre) Ply level -Weft yarn damage at nested yarn configuration (at the edge of the weft yarn) Micro level-fibre matrix debonding in the weft yarn


Marie Walschaerts
email: walschaerts.m@chu-toulouse.fr

Eve Leconte
email: leconte@cict.fr

Philippe Besse
email: philippe.besse@math.univ-toulouse.fr

Stable variable selection for right censored data: comparison of methods

Keywords: censored data, variable selection, survival trees, survival random forests, Lasso, Cox model, bootstrap

The instability in the selection of models is a major concern with data sets containing a large number of covariates. This paper deals with variable selection methodology in the case of high-dimensional problems where the response variable can be right censored. We focuse on new stable variable selection methods based on bootstrap for two different methodologies commonly used in survival analysis: the Cox proportional hazard model and survival trees. As far as the Cox model is concerned, we investigate the bootstrapping applied to two variable selection techniques: the stepwise algorithm based on the AIC criterion and the L 1 -penalization of Lasso. Regarding survival trees, we review two methodologies: the bootstrap node-level stabilization and random survival forests. We apply these different approaches to two real data sets, a classical breast cancer data set and an original infertility data set. We compare the methods on two criteria: the prediction error rate based on the Harrell concordance index and the relevance of the interpretation of the corresponding selected models, focusing on the original infertility data set. The aim is to find a compromise between a good prediction performance and ease to interpretation for clinicians. Results suggest that in the case of a small number of individuals, a bootstrapping adapted to L 1 -penalization in the Cox model or a bootstrap node-level stabilization in survival trees give a good alternative to the random survival forest methodology, known to give the smallest prediction error rate but difficult to interprete by non-statisticians. In a clinical perspective, the complementarity between the methods based on the Cox model and those based on survival trees would permit to built reliable models easy to interprete by the clinician.

Introduction

Problems of variable selection arouse a growing interest in the processing of data sets containing more and more variables. In the last twenty years, many methods of variable selection have been proposed to handle these high-dimensional problems, especially when the number of covariates p exceeds the number of observations n. To avoid a wrong estimation due to collinearity problems and to improve interpretation, the scientific community has developed tools to select the most relevant variables. A large literature concentrates on the case of the linear regression. A classical well-known method is the stepwise algorithm based on the Akaike Information Criterion (AIC). Recently, another field of research has focused on optimization problems, such as L 1 -penalty approaches. On the other way, tree-based algorithms provide an interesting alternative to handle non-parametrically a large number of covariates.

We consider the special case where the response variable is right censored. In this context, the Cox proportional hazards model [START_REF] Cox | Regression models and life tables (with discussion)[END_REF] has become the gold-standard tool for the statistical analysis, especially in the medical field. However, in the case of a large amount of covariates, it may be very unstable, even when stepwise selection or L 1 -penalization of Lasso are added to the classical procedure. The instability of these selection approaches have also been encountered in the context of the linear regression [START_REF] Harrell | Regression modelling strategies for improved prognostic prediction[END_REF][START_REF] Meinshausen | Stability selection[END_REF]. To remedy this problem, some authors have proposed to use bootstrapping to investigate the reliability of the choice of the variables in the final model. Bach [START_REF] Bach | Model-consistent sparse estimation through the bootstrap[END_REF] introduced stability in the selection by using bootstrapping in a Lasso algorithm, a method called Bolasso, and Meinshausen and Bühlmann [START_REF] Meinshausen | Stability selection[END_REF] proposed to improve the L 1 -penalization by randomizing the selection process of the covariates. These bootstrapped Lasso methods have only been considered in the framework of linear regression: we propose to extend them to the stabilization of the selection of covariates in a Cox model.

Alternatively to the Cox model, survival trees procedures permit to take into account non linear relationships between the censored variable and the covariates and yield easily interpretable classification rules to the clinicians, but they tend to overfit the data and suffer also from instability [5 -8]. A stabilization technique known to improve the prediction performance of a single tree consists in aggregating a family of tree models using boostrapped samples and a random selection of the covariates at each node of the trees. This procedure, called "random forests" [START_REF] Breiman | Heuristics of instability and stabilization in model selection[END_REF] has been adapted to the survival framework by Ishwaran et al. [START_REF] Ishwaran | Random survival forests[END_REF]. This Random Survival Forest (RSF) methodology is considered as the best modelization in terms of prediction performance but it it is not easy to interprete as it does not provide a single tree. On the other way, Dannegger [START_REF] Dannegger | Tree stability diagnostics and some remedies for instability[END_REF] introduced tools to control the stability of the selection of the covariates at each node of a tree and propose resampling techniques at node-level to stabilize the choice of the split. This methodology has the advantage to keep the simplicity of interpretation of a single final tree.

We review and compare these different stable variable selection methods derived on the Cox model or on survival trees on real data sets. We will base the comparison of the proposed methods on statistical criteria as the prediction accuracy but also on the capacity of the method to supply a final model whose clinical interpretation is easy and pertinent. As a matter of fact, the objective of the statistician is to provide to clinicians a model which is robust enough to permit a proper understanding of the relationships between the variable of interest and the covariates in a predictive and explanatory purpose.

We consider two real data sets. The first data set concerns the survival of breast cancer patients in relation with their gene-expression signature [START_REF] Van De Vijver | A gene-expression signature as a predictor of survival in breast cancer[END_REF] and is a very usual case encountered in survival analysis. This classical example will validate the different approaches discussed earlier. The second data set is much more original as it concerns couples consulting for infertility and very few studies were led in this domain. The censored variable of interest is the duration between the first visit of the couple and the birth of a alive child, which is difficult to model. The study takes into account several covariates of the couple or man and woman, but also the associated medical treatment according to the diagnosed causes of infertility. For example, a man with severe oligospermia will receive ICSI (IntraCytoplasmic Sperm Injection), while a woman with a dysovulation or with no tubal factor will enter an IVF (In Vitro Fertilization) program. In many couples, causes of infertility are not identified and physicians can not clearly answer the question of on-going pregnancy in couples.

In section 2, a presentation of the variable selection methods based on the Cox model is given, followed, in section 3, by a review of the selection methods based on trees algorithms. Section 4 compare the different approaches on the two data sets, giving the prediction erro rate and the selected variables in the final model for each approach. Concluding remarks and perspectives are presented in the last section.

Stable variable selection methods based on the Cox model Notations

We denote by T the variable of interest, which is a time of failure. We suppose that T may be right censored at a non-informative censoring time C such that C is independent of T conditionally on Z, a p-vector of covariates. The observations are therefore the possibly censored time X = min(T, C) and the censoring indicator ∆ = 1{T ≤ C}, where 1{•} is the indicator function. The observed sample is thus (X i , ∆ i , Z i ) i=1,...,n .

In survival analysis, the Cox model [START_REF] Cox | Regression models and life tables (with discussion)[END_REF] has become the most popular method to modelize the relationship between a survival time and one or more predictors (i.e. covariates). This model has the advantage to be semi-parametric in the sense that it does not require assumptions on the survival time distribution. Morevoer, it is of easy interpretation for clinicians in providing estimates of the effect of the covariates on survival time after adjustment on the other covariates.

In the Cox regression model in its simplest form, the hazard function for the failure time of an individual takes the form

λ(t|Z) = λ 0 (t) exp(β Z),
where β is a p-vector of unknown regression parameters and λ 0 (t) is an unknown baseline hazard function. Denote t (1) , . . . , t (k) the k ordered uncensored survival times. The parameter of interest β is estimated by maximising the logarithm of the partial likelihood function

L(β) = k i=1   β Z i -log   j∈R i e β Z j     (1)
where R i is the set of indices of the individuals at risk at time t (i) . Notice that the Cox model is a log-linear model based on the proportional hazard hypothesis, which stipulates that the ratio of the hazard functions of two individuals is constant over time.

In the following, we review three stable methods to select the relevant covariates using the Cox model: the bootstrap stepwise selection, the bootstrap Lasso selection and the bootstrap randomized Lasso selection.

Bootstrap stepwise selection (BSS)

In case of a large number of covariates, the selection of the predictors in the Cox model is usually made by a stepwise algorithm minimizing the Akaike Information Criteron (AIC). However, the stability of this method is questionable. Using a data splitting approach, Harrell et al. [START_REF] Harrell | Regression modelling strategies for improved prognostic prediction[END_REF] showed variation in the selected predictors. This result was confirmed by Chen and George [START_REF] Chen | The bootstrap and identification of prognostic factors via cox's proportional hazards regression model[END_REF] who applied the stepwise procedure on 100 bootstrapped samples from a study of acute lymphocytic leukaemia. The original model which was built from the entire set by a Cox regression model coincide in only 2% of the cases with the models selected from the bootstrapped samples. Sauerbrei and Schumacher [START_REF] Sauerbrei | A bootstrap resampling procedure for model building: application to the cox regression model[END_REF] developped a bootstrap selection procedure which combined the bootstrap method with stepwise selection in Cox regression. They examined the inclusion frequencies of the variables selected by the stepwise algorithm into the models derived from the bootstrapped samples and keep in the final model the variables for which the inclusion frequency exceeds a given cut-off value κ in (0, 1). The choice of κ is arbitrary. The authors applied the methodology to a data set on brain tumours including 447 patients and 12 covariates, using two values of κ (0.3 and 0.6). They showed that the choice of κ = 0.6 gives the same final model whatever the number of bootstrap samples (from 100 to 1000).

Bootstrap Lasso and bootstrap randomized Lasso selection

Lasso selection

An alternative to the stepwise selection procedure is the Lasso selection. Adapted by Tibshirani [START_REF] Tibshirani | The lasso method for variable selection in the cox model[END_REF] to the Cox model, the method estimates the β parameter via maximising the log partial likelihood function [START_REF] Cox | Regression models and life tables (with discussion)[END_REF] with the constraint p j=1 |β j | ≤ λ where λ is a regularization parameter. The Lasso constraint selects variables by shrinking estimated coefficients towards 0. This leads to coefficients exactly equal to zero and allows a parcimonious and interpretable model. The choice of a proper λ is sensitive and leads to variation in model selection as demonstrated by Meinshausen and Bühlmann [START_REF] Meinshausen | Stability selection[END_REF]. In a linear regression framework, they studied the stability of the selection in a gene expression data set. They introduced noise in variables in permuting all but 6 of the 4088 gene expressions. They showed that, when λ increases, the final model retained as well the 6 unpermuted variables as the irrelevant noised variables. Although cross-validation seems a natural solution to choose λ, it is not a good alternative for high-dimensional problems. As a matter of fact, the authors showed in their example that 14 permuted variables are still retained in the final model obtained in choosing λ by cross-validation. Moreover, Meinshausen and Bühlmann [START_REF] Meinshausen | High dimensional graphs and variable selection with the lasso[END_REF] prooved that if the number of covariates tends to infinity then the probability of selecting a wrong variable converges to 1.

Bootstrap Lasso selection (BLS)

In order to obtain a consistent model selection, Bach [START_REF] Bach | Model-consistent sparse estimation through the bootstrap[END_REF] proposed for the linear model a bootstrapped version of the Lasso, referred as the Bolasso. He defined the Bolasso model estimate as the final model composed by only the variables which are selected in all bootstrapped samples; in other words, using the terminology of the BSS, it corresponds to a cut-off value κ equal to 1. However, Meinshausen and Bühlmann [START_REF] Meinshausen | Stability selection[END_REF] showed with the permuted gene data set that even with random subsampling (which is a procedure close to the bootstrap [START_REF] Freedman | A remark on the difference between sampling with and without replacement[END_REF]), the Lasso algorithm could select irrelevant variables when λ is too large.

Bootstrap randomized Lasso selection (BRLS)

To deal with the choice of λ, Meinshausen and Bühlmann [START_REF] Meinshausen | Stability selection[END_REF] proposed a generalisation of the bootstrap Lasso procedure called bootstrap randomized Lasso where the covariates are penalized by different values randomly chosen in the range [λ, λ/α] with α in (0, 1). This turns out to estimate the β parameter with the constraint p j=1 β j W j ≤ λ. In practice, the set of covariates {Z j : j = 1, . . . , p} are weighted by the set {W j : j = 1, . . . , p} randomly generated where P (W j = α) = p w and P (W j = 1) = 1 -p w with p w in (0, 1).

The authors give no indication on the choice of p w and in the following we take p w = 0.5.

This procedure is very simple to implement and the authors showed that choosing α in the range of (0.2, 0.8) gives a consistent variable selection whatever the choice of the penalty λ. The authors explained that a low value of α will decrease the selection probabilities of irrelevant variables even if the penalty λ is large. As the BLS and BRLS methods are based on a penalization of the parameter β, it is necessary to first normalize the covariates such as

||Z j || 2 = ( n i=1 Z 2 ij ) 1/2 =
1 for all j in {1, . . . , p}.

Stable variable selection methods based on survival trees

Although they are not so popular than the Cox model, tree-based methods in survival analysis (the so-called survival trees) have known a great development in the last decades. They provide a good alternative to the Cox regression model in identifying covariates which play a role on the survival outcome and in predicting the individual risk of failure. In addition to be easy to interpret in a large frame of applications, survival trees methods can incorporate non linear effects, and also take into account interactions between covariates. First developed for basic classification trees, the Classification and Regression Trees (CART) algorithm of Breiman [START_REF] Breiman | Classification and Regression Trees[END_REF] is based on binary recursive partitioning. This is an iterative process which splits the data into two subgroups (daughter nodes) according to the value of one of the predictors. The splitting rule maximises the difference between nodes. Let denote {Z j : j = 1, . . . , p} the set of covariates. Formally, a split is induced by a question of the form "Is Z j ≤ b ?" for a continuous covariate where b is a cut-off value in the set of realizations of Z j . For nominal covariates with possible values in the set B = {b 1 , . . . , b r }, the question is of the form "Is Z j in S?" where S ⊂ B. To determine the best split s among all covariates in the current node h, a measure of improvement G(s, h) is evaluated, producing the most homogeneous daughter nodes. The best split s best is the optimal split among all possible splits in the set S h such as G(s best , h) = max s∈S h G(s, h). The process goes on until each node reaches a user-specified minimum node size and becomes a terminal node, or is homogeneous. To control the size of the tree, a stopping rule is used to prune the "large" tree containing "pure" nodes. The method employed by CART is called cost-complexity pruning. The complexity of a tree Θ is :

R cp (Θ) = R(Θ) + cp Θ ,
where R(Θ) is the raw error of measure (sum of the error of measure over the terminal nodes), Θ is the number of terminal nodes, and cp is an arbitrary penalty weight between 0 and ∞ called the complexity parameter. A sequence of nested trees is built. The final tree is the smallest tree for which R cp (Θ) is minimized. To deal with the choice of cp, cross-validation techniques can be used to determine the optimal size of the tree.

Various splitting and pruning approaches have been proposed to adapt regression trees to survival data. Gordon and Olshen [START_REF] Gordon | Tree-structured survival analysis[END_REF] used the Wasserstein metrics to measure the distance between two Kaplan-Meier estimates of the survival distribution. The split criterion chooses the predictor (and if necessary the cut-off) which maximises this distance for the left and right daugther nodes. To prune the survival tree, they generalized the cost-complexity measure for censored data. A more usual way to measure the difference between survival curves is the logrank test statistic. This was done by Ciampi et al. [START_REF] Ciampi | Stratification by stepwise regression, correspondence analysis and recursive partition: A comparison of three methods of analysis for survival data with covariates[END_REF], Segal [START_REF] Segal | Regression trees for censored data[END_REF] and LeBlanc and Crowley [START_REF] Leblanc | Survival trees by goodness-of-split[END_REF] who suggested to use the logrank test statistic as a between-node heterogeneity measure. As an alternative to the standard cost-complexity pruning approach using proportional hazards martingale residuals as error measures, LeBlanc and Crowley [START_REF] Leblanc | Survival trees by goodness-of-split[END_REF] proposed a "goodness-of-split" complexity based on the sum of the standardized splitting logrank test statistics of the internal nodes. Others authors suggested to use split criteria based on the likelihood function. Davis and Anderson [START_REF] Davis | Exponential survival trees[END_REF] assume that the survival function within each node is an exponential function with a constant hazard, while LeBlanc and Crowley [START_REF] Leblanc | Relative risk trees for censored survival data[END_REF] only assume proportionality for the hazard functions of two daughter nodes. These latter used for estimation the full or partial likelihood function in the Cox proportional hazards model.

Instability in the selection of covariates by regression trees has been observed and demonstrated by many authors [5 -8] This instability may be due to an overfitting of data. The variance observed may also come from arbitrary cutpoints defined by the dichotomization of continuous covariates. We review in the following two remedies to instability of survival trees: the bootstrap node-level stabilization proposed by Dannegger [START_REF] Dannegger | Tree stability diagnostics and some remedies for instability[END_REF] and the random survival forests of Ishwaran et al. [START_REF] Ishwaran | Random survival forests[END_REF].

Bootstrap node-level stabilization (BNLS)

Dannegger [START_REF] Dannegger | Tree stability diagnostics and some remedies for instability[END_REF] proposed a bootstrap node-level stabilization procedure for survival trees. The algorithm consists, at node h, in drawing bootstrapped samples from the original set, and for each of them, in finding the best split. The split which appears the most of the time at the node h is selected. For a continuous variable, the cut-off value b chosen in the set of realizations of the split variable is the median of all the b-values proposed at each bootstrap. As Dannegger [START_REF] Dannegger | Tree stability diagnostics and some remedies for instability[END_REF] did not propose a choice of the cut-off for a categorical variable, we decide to affect to the daughter node the level of the categorical variable which was mainly chosen by the boostrapped samples. To find the optimal complexity parameter cp to prune the tree, we use a ten-fold cross-validation procedure as suggested by Dannegger [START_REF] Dannegger | Tree stability diagnostics and some remedies for instability[END_REF].

Using simulated data, Dannegger [START_REF] Dannegger | Tree stability diagnostics and some remedies for instability[END_REF] compared his methodology (for 100 bootstraps at each node) to the CART algorithm and to the bagging method and found that the BNLS reduces the predictor error rate from 26% (the value obtained with CART) to 6%. The bagging method has the best performance with a prediction error of 3%. However the model obtained with this latter method is not easy to use and interpret by nonstatisticians, unlike BNLS.

Random Survival Forests (RSF)

In order to stabilize the trees obtained by the CART procedure, Breiman [START_REF] Breiman | Bagging predictors[END_REF] proposed the bagging method (so called for bootstrap aggregating) based on a family of random trees: multiple versions of predictor models obtained from bootstrapped samples are aggregated in order to construct a robust model estimator. Bagging was then adapted to the survival framework by Efron [START_REF] Efron | Censored data and the bootstrap[END_REF] and Akritas [START_REF] Akritas | Bootstrapping the kaplan-meier estimator[END_REF]. Another way to improve stability of trees is the boosting, developped by Freund and Schapire [START_REF] Freund | A short introduction to boosting[END_REF]. As for the bagging, the boosting consists in aggregating a family of models. Each tree is built iteratively from a weighted sample (an individual who is misclassified gains weight and an individual who is classified correctly looses weight) and then evaluated according to its ability to classify data. Considered as better than the bagging, the boosting is yet limited when data are too noised: the algorithm gives a heavy weight to noised data which leads to a bad overfit [START_REF] Dietterich | Ensemble methods in machine learning[END_REF].

Breiman [START_REF] Breiman | Bagging predictors[END_REF] proposed a random selection approach which combines the bagging method with a random selection of the covariates at each node of the tree. This methodology, called random forests, is more stable than the two previous ones [START_REF] Breiman | Bagging predictors[END_REF][START_REF] Freund | A short introduction to boosting[END_REF]. The method was adapted to the survival framework in an approach called "random survival forests" by Ishwaran [START_REF] Ishwaran | Random survival forests[END_REF]. Bootstrapped samples were drawn from the original data set. Notice that each bootstrapped sample excludes on average 37% of the data, a set called out-of-bag (OOB) data. For each bootstrapped sample b in (1, . . . , ntree), a survival tree is built: at each node, a subset of covariates is randomly selected among all the covariates. The splitting process continues under the constraint that each terminal node contains at least a fixed number of unique event times. The RSF algorithm computes an ensemble estimate for the cumulative hazard function (CHF), which is used as a predictor. For each terminal node h, let t (1,h) < . . . < t (N h ,h) be the distinct ordered uncensored event times, and define d l,h and Y l,h as the number of events and individuals at risk at time t (l,h) , for l = 1, ..., N h respectively. The CHF estimator for the node h is the Nelson-Aalen estimator

Ĥh (t) = t (l,h) ≤t d l,h Y l,h
•

Each tree provides such a CHF estimate at each terminal node. Let Ĥb (t|z) denote the cumulative hazard estimate for tree b conditionnaly on the covariate z. To determine the CHF estimate for individual i with covariate vector z i obtained from the tree b, drop z i down the tree. It will fall in a unique terminal node h. So we have Ĥb

(t|z i ) = Ĥh (t) if z i ∈ h. Let I i,b = 1 if i is an OOB point for b and 0 otherwise. The OOB ensemble CHF estimator for i is: Ĥ * e (t|z i ) = ntree b=1 I i,b Ĥb (t|z i ) ntree b=1 I i,b • (2)
Note that the above estimator is obtained by avering only over bootstrap samples in which i is excluded.

Comparison of the methods

In order to compare the different variable selection approaches, we apply them to two real data sets, the first concerning breast cancer and the second male fertility. The five following procedures have been compared: the bootstrap Cox stepwise procedure (BSS), the bootstrap Cox Lasso procedure (BLS), the bootstrap Cox randomized Lasso procedure (BRLS) with three different values of α (0.2, 0.4 and 0.6), the bootstrap nodelevel stabilization procedure (BNLS) and the random survival forest method (RSF).

Software

We use the methods implemented in R software (CRAN). We adapted the bootstrap to the stepwise Cox algorithm and to the R package penalized. We also modified the R package rpart in order to introduce bootstrap at node-level in the building of the tree (BNLS procedure). Instead of the default exp method which maximises an exponential likelihood, we used as a splitting criterion in rpart the logrank test statistic, following the recommandations of Radespiel-Troger et al. [START_REF] Radespiel-Troger | Association between split selection instability and predictive error in survival trees[END_REF]. They showed that the prediction error of the tree is related to the instability of the covariate selection at each node.

Compared to other splitting criteria on a real survival data set on gallbladder stones, the logrank test statistic yields the lowest prediction error (evaluated by the Brier score). The R package randomSurvivalForests provides four different splitting rules: logrank, conservation of events, logrank score, and approximate logrank (for more details, see Ishwaran and Kogalur [START_REF] Ishwaran | Random survival forests for R[END_REF]). We used the default splitting criterion i.e. the logrank test statistic. Concerning the number of bootstrapped samples, we take N = 100 for the Cox model based procedures, the default value N = 1000 for RSF, and N = 1000 for the BNLS procedure.

Comparison criteria

We compare the five procedures on the basis on their prediction performance but also on the usefulness of the model for non-statisticians. The prediction performance is measured by the prediction error rate. A model will be considered as useful if the statistical method involved is the most adapted to the problematic of the clinician and if the results are easy to interprete.

Comparison of the prediction error rates

We compare the stability of the five procedures using a statistical criterion: the prediction error rate, which allows to quantify the prediction performance of the final model selected.

To this aim, the original data set is divided into two subsets, the training set and the test set. The procedures, computed on the training set, give a final model which is then applied on the test set to calculate the prediction error rate. We obtain the prediction error rates in three steps:

1. We graphically chose sensible values of the cut-off κ for the methods based on the Cox model and also values of λ for BLS and BLRS. For that, considering the whole set, we plot the values of κ for each selected covariable for BSS and we plot the values of κ with respect to values of λ for BLS and BLRS. The aim is to identify the most relevant variables whose frequencies of inclusion κ are larger than those of the other covariates whatever the value of λ. Then we chose values of λ and κ which maximize the gap between the two subgroups of covariates. For BNLS, the choice of the optimal complexity parameter cp is achieved by cross-validation techniques.

2. We then apply the procedures on the training set with the values of λ, κ and cp chosen at step 1. As far as the methods based on the Cox model are concerned, for each method, we adjust on the training set a simple Cox model including the subset of covariates for which the inclusion frequencies are greater than κ. This permits to obtain the estimations of the β parameter in order to calculate the predicted outcomes.

3. We finally calculate the prediction error rates on the test set for each method.

The most commonly used prediction error rate for survival models is based on the Harrell's concordance index C [START_REF] Harrell | A new distribution-free quantile estimator[END_REF]. To compute the C-index, the observed and predicted outcomes are compared, as follows:

• Form all possible pairs of observations.

• Eliminate those pairs where the shorter survival time is censored, and also pairs where the survival times and the censoring indicators are equal.

• For each permissible pair, count 1 if the shorter survival time has the worse predicted outcome and count 0.5 if the predicted outcomes are tied. Sum over all permissible pairs.

• The C-index is the ratio of the latter sum over the number of permissible pairs.

The prediction error rate is defined as 1 -C and is in {0, 1}. Note that a value of 0.5 is not better than random guessing. For the Cox model based procedures, the prediction for an individual i with covariate vector Z i is based on the linear predictor β Z i . Individual i has worse predicted outcome than j if β Z i > β Z j . The predictions of the survival trees based procedures are derived from the CHF estimator. For the RSF procedure, they are derived from the OOB ensemble CHF estimator of formula (2). Let t 1 , . . . , t m be the unique times in the test set. Individual i is said to have a worse predicted outcome than j if

m k=1 Ĥ * e (t k |Z i ) > m k=1
Ĥ * e (t k |Z j ).

A total of 30 training sets and test sets were drawn to obtain a sample of error rates. Boxplots of the error rates are presented in order to compare the variability of the prediction performances.

Comparison of the selected covariates

We compare the set of covariates selected by the five different approaches applied on the entire set. For BSS, BLS and BRL, we take the values of the cut-off κ and of the penalties λ determined above. The BNLS algorithm provides the final tree which has been pruned with a fixed complexity parameter cp. For the RSF procedure, the impact of the variables is measured by their "variable importance" (VIMP), which is the difference between the prediction error obtained with the original RSF procedure and the prediction error obtained using randomized assignments whenever a split for the considered variable is encountered [START_REF] Ishwaran | Random survival forests[END_REF].

Results for the breast cancer data set 4.3.1 Description

The breast cancer data set contains the metastasis-free survival times from the study of Vijver et al. [START_REF] Van De Vijver | A gene-expression signature as a predictor of survival in breast cancer[END_REF] who classified a serie of 295 patients with primary breast carcinomas as having a gene-expression signature associated with either a poor or a good prognosis. We restricted the study to the 144 patients who had lymph node positive disease. Vijver et al. [START_REF] Van De Vijver | A gene-expression signature as a predictor of survival in breast cancer[END_REF] evaluated the predictive value of the gene-expression profile of patients for the 70 genes previously determined by Veer et al. [START_REF] Van 't Veer | Gene expression profiling predicts clinical outcome of breast cancer[END_REF] based on a supervised learning method. The data set can be found in the R package penalized. Five clinical risk factors and 70 gene expression measurements found to be prognostic for metastasis-free survival have been recorded. The censoring rate is 66%. The variables in the data set are:

• time: metastasis-free follow-up time,

• event: censoring indicator (1 = metastasis or death; 0 = censored),

• diam: diameter of the tumor (two levels),

• N: number of affected lymph nodes (two levels),

• ER: estrogen receptor status (two levels), • grade: grade of the tumor (three ordered levels),

• age: age of the patient at diagnosis,

• TSPYL5...C20orf46: gene expression measurements of 70 prognostic genes.

Prediction error rates

The Cox stepwise algorithm did not converge for the BSS procedure, probably due to too many covariates and not enough events in the boostrapped samples. As suggested by Meinshausen and Bühlmann [START_REF] Meinshausen | Stability selection[END_REF], we search graphically the penalty λ which leads to the best split between covariates in the BLS and BRLS procedures. However, the graphical determination of λ is not so easy as in Meinshausen and Bühlmann [START_REF] Meinshausen | Stability selection[END_REF] which presents results on an ad-hoc example. Figure 1 suggests that selecting between four and five covariates seems to be the most relevant whatever the value of λ and κ. So we decide to choose a value of λ = 0.4 and a value of κ = 0.2 for BLS and BRLS whatever the value of α. The complexity parameter cp of the BNLS procedure is obtained by a ten-fold cross-validation procedure: we chose graphically cp = 0.002. Figure 2 presents the boxplots for the different approaches based on 30 iterations and the numerical summaries are in table 1. As expected, RSF gives the lowest mean, median and standard deviation of the prediction error rates. We verify that the single survival tree has the worse prediction performance, but the BNLS procedure is hardly better. Whatever the value of α, the BRLS procedure shows a large dispersion of error rates compared to the BLS procedure which has a standard deviation of 0.06, similar to this of RSF. However, a good compromise between dispersion and mean of error rates is observed for BRLS with α = 0.2 or α = 0.4.

Variables selected in the final model

Figure 3 shows the covariates selected by the Cox based procedures BLS and BRLS according to the choice of the inclusion frequency κ. Only 23 covariates of the 75 initial covariates appear in the selection done by the procedures applied to the boostrapped samples (i.e. with κ = 0) and no clinical covariate is selected. For κ = 0.2, six covariates are selected by BLS whereas only four are selected by the BRLS procedures whatever the value of α (these four covariates are a subset of the previous six). We notice a clear gap between the four first selected variables and the other ones in the BRLS procedures. We can observe that the BLS and BRLS procedures find the same four most relevant covariates but not in the same order: PRC1, ZNF533, QSCN6L1 and IGFBP5.1.

For the tree-based procedures RSF and BNLS, figure 4 and5 show that the same most discriminant variable is ZNF533. BNLS selects three other covariates among the first most important variables selected by RSF: COL4A2, PRC1 and N. The selection with a single survival tree procedure exhibits the same first discriminant variable ZNF533 followed by PRC1, RP5.860F19.3, HRASLS, IGFBP5, and SCUBE2. But the additional selected covariates by the CART algorithm never appear in the stable selection procedures. We can also notice that no clinical factor is included in the survival tree built by CART procedure unlike the BNLS and RSF procedures.

As far as the comparison of the Cox based methods and survival trees methods is

BBC3 RP5.860F19.3 EXT1 Contig35251_RC AYTL2 GSTM3 PECI.1 PITRM1 C20orf46 AA555029_RC SCUBE2 ECT2 CDC42BPA UCHL5 PECI IGFBP5 Contig20217_RC SERF1A MTDH Contig40831_RC RFC4 CENPA ALDH4A1 GPR180 NUSAP1 FLT1 OXCT1 RUNDC1 NM_004702 GNAZ STK32B ESM1 C9orf30 LGP2 EGLN1 ORC6L RTN4RL1 MS4A7 Contig32125_RC Contig63649_RC COL4A2 ZNF533 IGFBP5.1 QSCN6L1 PRC1
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8

κ FGF18 Contig32125_RC FLT1 RUNDC1 GPR180 RAB6B UCHL5 MTDH TGFB3 PECI.1 EGLN1 ALDH4A1 ECT2 GMPS PECI Contig40831_RC STK32B GNAZ OXCT1 CDC42BPA MELK DTL C9orf30 SCUBE2 RFC4 ESM1 IGFBP5 Contig63649_RC LGP2 CENPA RTN4RL1 NM_004702 COL4A2 ORC6L MS4A7 NUSAP1 IGFBP5.1 QSCN6L1 ZNF533 PRC1
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.1 0.2 0.3 0.4

κ 1. BLS, λ = 0.4 2. BRLS, α = 0.2, λ = 0.4 DIAPH3 FGF18 WISP1 GSTM3 UCHL5 TGFB3 FBXO31 SLC2A3 AP2B1 C20orf46 RP5.860F19.3 EXT1 ECT2 CDC42BPA PECI MELK GPR126 PECI.1 PITRM1 ALDH4A1 Contig40831_RC EGLN1 SCUBE2 GNAZ OXCT1 GMPS MTDH DTL RUNDC1 Contig35251_RC ESM1 GPR180 RFC4 C9orf30 NM_004702 STK32B MS4A7 CENPA Contig32125_RC Contig63649_RC LGP2 NUSAP1 RTN4RL1 COL4A2 IGFBP5 ORC6L QSCN6L1 IGFBP5.1 ZNF533 PRC1
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.1 0.2 0.3 0.4

κ SCUBE2 EXT1 PECI MTDH TGFB3 MELK DTL PECI.1 C20orf46 GNAZ RUNDC1 Contig35251_RC ECT2 GPR180 C9orf30 EGLN1 ESM1 Contig32125_RC FLT1 STK32B ALDH4A1 GMPS Contig40831_RC IGFBP5 LGP2 NM_004702 ORC6L MS4A7 COL4A2 CENPA Contig63649_RC RFC4 NUSAP1 RTN4RL1 ZNF533 IGFBP5.1 QSCN6L1 PRC1
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.1 0.2 0.3 0.4 0.5 concerned, we notice similarity between the selected covariates: for example, ZNF533 is found in all models. Regarding clinical factors, none was selected in Cox models contrary to the survival trees based procedures. As they can take into account interactions between covariates, tree based procedures show that N, the number of affected lymph nodes and the age of the patient at diagnostic are relevant predictors in metastasis-free survival.

Notice that these two factors are also known to be clinically relevant.

Results for the fertility data set 4.4.1 Description

The fertility data were obtained from 2138 couples consulting for male infertility during the period from 2000 to 2004 at Toulouse Male Sterility Center (TMSC) located in Universitary Public Hospital (France). Patients were followed from entry and during treatment by an andrologist specialist until either discontinuation of treatment or delivery of an alive infant. The maximum follow-up duration is 9 years. The outcome assessment was based on the delivery of an alive infant obtained at TMSC (pregnancies after medical treatment -medecine and/or surgical treatment, or assisted reproductive technologies (ART) -as well as spontaneous pregnancies). The event considered here is the birth of an alive infant and right-censored events correspond to miscarriage or loss to follow-up. The "survival" time is the delay in months from the first visit of the couple to the birth of its alive infant. We will work on the subset of the 1773 couples with covariates without missing values. 40% of the couples succeed in their parental project, leading to a censoring rate of 60 %.

FGF18 RAB6B Contig32125_RC BBC3 FLT1 OXCT1 ORC6L NMU Contig20217_RC RUNDC1 AP2B1 LGP2 DIAPH3.1 SLC2A3 ESM1 C20orf46 DIAPH3 RP5.860F19.3 TGFB3 ER NUSAP1 GMPS WISP1 MELK IGFBP5 PECI.1 MMP9 ALDH4A1 GPR126 Grade QSCN6L1 C16orf61 HRASLS GPR180 CDCA7 RFC4 NM_004702 UCHL5 MCM6 Contig63649_RC MS4A7 N IGFBP5.1 COL4A2 Age PRC1 ZNF533
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.000 0.005 0.010 0.015 0.020 0.025 0.030 In agreement with clinicians, we decided to keep 32 covariates, among which:

• for the man: age, clinical investigation including medical histories as: histories of orchitis, epididymitis, urogenital infections (UI), sexually transmitted infection (STI), inguinal hernia, testis trauma, cryptorchidism (testicular maldescent), cancer, as well as testis, epididymis or vas deferens surgery, and clinical examination as: aspects of scrotum, testicular migration (presence of the testes into the scrotum and their position), epididymides and vasa deferentia evaluation, presence of a varicocele or hydrocele, testicular volumes, and if he received a non-ART treatment (both medical and surgical, including pharmacological and hormonal treatment, varicocelectomy, and vas deferens/epididymis surgery).

• for the woman: age, tubal factor, ovarian factor, cervical factor, ovulation, and if she received a non-ART treatment (hormonal treatment).

• for the couple: fecundity type (primary or secondary), infertility duration, type of ART including IUI (intra-uterine insemination), IVF (in vitro fertilization), ICSI (intracytoplasmic sperm injection) using male sperm cells, and ART with donor sperm.

Prediction error rates

As for the breast cancer data set, we chose λ and κ graphically. A first sight of figure 6 reveals that the contribution of the BRLS procedure with regard to the BLS procedure for the choice of the cut-off values is not obvious. For BLS and BRLS, only two variables seem to be the most relevant factors whatever the value of λ. But for BLS and BLRS with α = 0.6, five other covariates follow the two first when λ = 0.4, for a value of κ = 0.3 whereas for BRLS with α = 0.2 and α = 0.4, no group of variables seems relevant. So we decided to choose a fixed penalty λ = 0.4 for a value of κ = 0.3 for BLS and BRLS.

For the BSS procedure, we decide to choose a value of κ = 0.5 which seems to split the set of covariates in two distinct subgroups as suggested by figure 8. Thus, the first eight variables seem to be the most relevant to predict the duration to the birth. For the BNLS procedure, we observe by a 10-fold cross-validation procedure that the optimal value of the complexity parameter is cp = 0.0035.

The boxplots of the prediction error rates for the five procedures are presented in figure 7 and table 2 shows the summary statistics. Notice that the mean and median values obtained are higher (between 0.41 and 0.45) than those obtained from the data set on breast cancer which reflects the difficulty to predict the delay to the birth. However, the variabilities of the error rates are lower. As seen previously, the RSF method gives the best predictive model. However, it appears that RSF is similar in dispersion with other procedures. Moreover, we can notice that the BNLS procedure is not better than a single survival tree and the same remark can be done between the BSS procedure compared to a single stepwise Cox regression. These results can be explained by the size of the sample, which is sufficient to produce low error rates without bootstrapping. Regarding BLS and BRLS, these procedures show less variation in error rates but their means and medians are close to 0.5 which suggests that these models do no better prediction than random guessing.

Selected variables in the final model

If we compare the covariates selected by the different approaches, we can see in figure 8 that the first four selected covariates do not differ for the BLS and BRLS procedures for a value of κ = 0.3: we find tubal factor, IUI, sperm donor and infertility duration. The additional covariates included by BLS are epididymis, varicocelis, inguinal hernia, fecundity type, testicular trauma and testicular volume. For the BSS procedure, for a value of κ = 0.5, we find in the selected covariates the first four covariates selected by BLS but also female age, testicular volume, male treatment, varicocelis, testicular trauma, scrotum, female treatment and epididymis. We observe also that the BSS procedure includes more variables than the BLS and BRLS procedures, which leads to a lower prediction error rate. We can notice that the cut-off value κ = 0.3 appears very clearly in the graphs of figure 8 for the BRLS procedures with α = 0.4 and α = 0.6, suggesting that the BRLS procedure is a good stable variable selection procedure. However, only four covariates are selected, which leads to poor prediction performances.

If we adjust a single Cox stepwise regression, we find the same covariates selected by the BSS procedure except for the variable scrotum. However, epididymis, testicular trauma and male treatment are not statistically significant at the 5% level in the Cox stepwise model.

Regarding the tree-based RSF and BNLS procedures, we observe in figures 8 and 9 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.2 0. q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0. q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0. q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0. q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.000 0.005 0.010 0.015 0.020 0.025 is the same for the two types of procedures but the following splits are different. This may be explained by the fact that the tree based procedures take into account interactions contrary to the Cox model based procedures. Moreover, the covariates selected for the BNLS procedure are found in the most important variables in RSF. The variables selected by a single survival tree are almost the same than those selected by BNLS, which can be explained by the large sample size of the fertility data set.

Discussion

Computational requirements

On a 2.4GHz processor, for the fertility data set with 1773 individuals and 32 covariates, running the RSF procedure is 200 times longer than running a single survival tree from the CART procedure (which takes about 1 second), whereas the stepwise algorithm takes about 10 seconds. Compared to the RSF method, the BLS and BRLS procedures have similar computing times. The BNLS procedure is much longer than RSF (4 times) as a bootstrap is realized at each node. Finally, the BSS procedure, for which a bootstrap sample is used at each stepwise, is the most expensive procedure in time (8 times longer than RSF). Thus, the running times seem reasonable for each procedure taken separately. However, a systematic comparison of errors becomes heavy, as well as the running time to exhibit the optimal complexity parameter cp in the cross-validation procedure for BNLS, and also the penalty λ for the BLS and BRLS methods.

Prediction error rates

Compared to the data set on breast cancer, the prediction error rates obtained from the fertility data set are higher, suggesting that data on fertility are more complex and difficult to study and that finding a good prediction model is a difficult task. Moreover, the sizes of the two data sets are very different, which could lead to observe discrepancy in the results. However, we find from the two data sets that the RSF procedure is the best procedure for prediction and that it allows to exhibit the most relevant variables to explain the survival durations. The Cox model and RSF have been previously compared by Omurlu et al. [START_REF] Omurlu | The comparisons of random survival forests and cox regression analysis with simulation and an application related to breast cancer[END_REF] using the Harrell'concordance index. On the basis on Monte Carlo simulations, they show that the Cox model has the best predictive performance whatever the size of the sample (n = 50, 100, 250, 500). However, on a real data set on breast cancer, they find that RSF has the lowest prediction error rate. These contradictory results are not surprising if we consider the fact that the simulated data have been generated from a Cox model. Moreover, the RSF procedure is easy to use and does not require the choice of tuning values as do the BSS, BLS, BRLS and BNLS procedures. However, even if the selected covariates are identified and sorted by their importance, as no final tree is provided, the RSF results stay a black-box not easy to interpret and use for clinicians. We find that the BRLS procedure, whatever the value of α, does not seem to improve the BLS procedure. Although the gap between the selected variables appears more clearly for BRLS, the final model does explain these results. Moreover, the good results obtained with this method by Meinshausen and Bühlmann [START_REF] Meinshausen | Stability selection[END_REF] may come from the fact that they presented an ad-hoc example. As noticed by these authors, we find on our data sets that choosing a value of κ equal to 1, as suggested by Bach [START_REF] Bach | Model-consistent sparse estimation through the bootstrap[END_REF], is too restrictive.

Even if some authors [START_REF] Chen | The bootstrap and identification of prognostic factors via cox's proportional hazards regression model[END_REF][START_REF] Sauerbrei | A bootstrap resampling procedure for model building: application to the cox regression model[END_REF] showed that the bootstrap method adapted to the stepwise algorithm improves the stability of the variable selection, our results suggest no improvement with the BSS procedure. As a matter of fact, the BSS algorithm does not converge on the breast cancer data set containing too few events with respect to the number of the covariates. On the other side, for the fertility data set whose size is sufficient to assure the convergence of the algorithm, it does not give better results than a single stepwise Cox model.

As far as the BNLS procedure is concerned, it seems not perform better than a single survival tree on the two data sets, contrary to the results found by Dannegger [START_REF] Dannegger | Tree stability diagnostics and some remedies for instability[END_REF]. It can be explained by the fact that it is difficult to tune a sensitive value for the complexity parameter by cross-validation. It may also result from the sufficient size of the fertility data set.

Selected variables in the final model

The breast cancer data set was originally analyzed by Vijver et al. [START_REF] Van De Vijver | A gene-expression signature as a predictor of survival in breast cancer[END_REF] who included in a multivariate Cox regression model the clinical factors and a prognosis variable based on the gene profile. Indeed, the 70 genes expressions were not included as independent covariates in the Cox model, as we do in our analysis. We notice that the procedures based on the Cox model selected no clinical factors in the breast cancer data set. It can result from the fact that the genes which were introduced in the model are those whose expression is the most correlated positively or negatively with the breast cancer survival duration among the multitude of genes studied. Thus, these genes have a more important weight in the multivariate analysis than clinical factors. In the BNLS and RSF procedures, which take naturally into account interactions between covariates, two clinical factors appears in the selected set of covariates. As it is not possible to include interactions in the Cox model in a high-dimensional context, an alternative would be to force the inclusion of some clinical factors clearly identified by the clinicians in the procedures based on the Cox model.

Regarding the fertility data set, we obtain different final models with the procedures based on the Cox model and those based on survival trees. However, all the selected covariates are relevant clinical factors in reproductive issues. Indeed, the BSS, BNLS and RSF procedures showed that female age is an important factor, with a cut-off value of 34.5 years, which is similar to results observed in other studies [34 -36]. Moreover, the infertility duration is selected by all the procedures, with a cut-off of 33.5 months in the BNLS method. It is not surprising as it is well-known that duration of infertility is predictive in the occurrence of pregnancy [START_REF] Ducot | Male factors and the likelihood of pregnancy in infertile couples. II. study of clinical characteristics -practical consequences[END_REF]. It is also interesting to notice that the BLS, BSS and BNLS procedures found a relationship between the live-birth and male clinical factors whereas BRLS included very few covariates which are essentially limited to female aspect and ART treatment. In most of the publications about reproductive issue in infertile couples, authors analysed the survival event of live-birth according only to female aspect and type of ART treatment [38 -41]. However, we observe in our study that varicocelis is included in most of the procedures, which can be explained by the fact that varicocelis has an impact on an impaired spermatogenesis [START_REF] Redmon | Varicocele: the most common cause of male factor infertility?[END_REF][START_REF] Health | The influence of varicocele on parameters of fertilityin a large group of men presenting to infertility clinics[END_REF]. The same observation could be done for the testicular volume which seems to be an important variable in the BNLS procedure, but with less impact in the other final models.

Concluding remarks

As far as the breast cancer data set is concerned, the objective is reasonably well achieved. The method based on randomized Lasso selection combined with bootstrapping provides prediction error rates very similar to those obtained by the RSF method, in association with an easily interpretable final model for the medical field. Regarding the procedure that combines bootstrap and survival trees, this method especially allows to identify interactions between the relevant covariates.

On the contrary, the data set on fertility seems to be much more difficult to analyze. The bootstrapping adapted to the various procedures based on the Cox model or on survival trees seems not improve the standard procedures (single Cox stepwise selection or single survival tree). Moreover, the wide variations observed in the prediction error rates may be due to the weak predictive performance of the selected covariates.

Finally, these results suggest that the Cox and tree based procedures should be performed in a complementary way to identify the most relevant covariates and provide to clinicians a stable and reliable model. Each procedure shows indeed a particular interest, either in terms of its prediction performance, either in the selection of the relevant covariates.

6 Figure 1 :

 61 Figure 1: Inclusion frequencies κ with respect to λ for Bootstrap Lasso Selection and Bootstrap Randomized Lasso Selection.

Figure 2 :

 2 Figure 2: Boxplots of the prediction error rates for seven methods : Bootstrap Lasso selection, Bootstrap Randomized Lasso selection (three values of α), Bootstrap node-level selection, Random Survival Forest and a simple survival tree.

κ 3 . 2 Figure 3 :

 323 Figure 3: Selected covariates with respect to the inclusion frequency κ for Bootstrap Lasso Selection and Bootstrap Randomized Lasso Selection

Figure 4 :

 4 Figure 4: The final tree obtained by the Bootstrap Node-Level Selection.

Figure 5 :

 5 Figure 5: Importance of variables for Random Survival Forest.

6 Figure 6 :

 66 Figure 6: Inclusion frequencies κ with respect to λ for Bootstrap Lasso Selection and Bootstrap Randomized Lasso Selection.

Figure 7 :

 7 Figure 7: Boxplots of the prediction error rates for nine methods : Bootstrap Stepwise selection, Booststrap Lasso selection, Bootstrap Randomized Lasso selection (three values of α), Bootstrap node-level selection, Random Survival Forest, Cox stepwise selection and a simple survival tree.

Table 1 :

 1 Mean, standard deviation and median of error rates for the different procedures: Bootstrap Lasso selection, Bootstrap Randomized Lasso selection (three values of α), Bootstrap node-level selection, Random Survival Forest and a single survival tree.

		Mean Standard deviation Median
	BLS	0.319	0.066	0.310
	BRLS α = 0.2 0.309	0.075	0.292
	BRLS α = 0.4 0.309	0.081	0.287
	BRLS α = 0.6 0.311	0.086	0.302
	BNLS	0.376	0.073	0.394
	RSF	0.279	0.062	0.286
	TREE	0.389	0.0741	0.382

Table 2 :

 2 Mean, standard deviation and median of error rates for the different procedures: Bootstrap Stepwise selection, Booststrap Lasso selection, Bootstrap Randomized Lasso selection (three values of α), Bootstrap node-level selection, Random Survival Forest, Cox stepwise selection and a single survival tree.

		Mean Standard deviation Median
	BSS	0.453	0.017	0.454
	BLS	0.471	0.017	0.470
	BRLS α = 0.2 0.483	0.013	0.482
	BRLS α = 0.4 0.484	0.013	0.482
	BRLS α = 0.6 0.482	0.015	0.480
	BNLS	0.447	0.021	0.446
	RSF	0.413	0.017	0.414
	STEP	0.451	0.015	0.450
	TREE	0.445	0.018	0.443

  FigureThe final tree obtained from the Bootstrap Node-Level selection that the selected covariates are substantially different from the most relevant covariates selected by the Cox based procedures. The first split, which corresponds to female age,

	5. BRLS α = 0.6, λ = 0.2		6. RSF
	Figure 8: Second part: Selected variables for Bootstrap Stepwise Selection, Bootstrap
	Lasso Selection, Bootstrap Randomized Lasso Selection and importance of variables for
	Random Survival Forest.
		female_age>=34.5 | 1773/717
	observations/events
			infertility_duration>=33.5
	RR=1.17		1372/622
	(401/95)
	female_treatment< 0.5		testicular_volume>=0.5
				varicocelis< 0.5	905/443
			male_age< 33.5	347/160	RR=1.52
				RR=1.63	(558/283)
			ICSI< 0.5	(97/55)
	RR=1.13	RR=0.50	RR=0.77
	(290/138)	(177/41)	(178/65)
			RR=1.29 (62/32)	RR=4.68 (10/8)


Venceslas Biri
email: biri@univ-mlv.fr

Single Scattering Effects for Computer Games

Keywords: Single scattering, Real Time, Hardware rendering

This paper addresses the rendering of single scattering effects such as glows and shafts of light, along with volumetric shadows induced by shadow casters in the participating media in real-time. Our method is easy to integrate in a video game graphics engine using the shadow volume technique since it requires only a little additional texture memory and is implemented with simple shaders. Realistic images can be produced in real-time for usual graphic scenes and at a high level framerate for complex scenes, allowing changes in the properties of participating medium, animations of objects and even light sources movements.

I. GOALS AND CONTRIBUTIONS

In this paper we show how to render a complete single scattering illumination model. The originality of the model is the introduction of the indirect single scattering. Indeed, we will call direct single scattering the effect of the first scattering of light along view rays. The indirect single scattering is the same effect but along illumination rays of any point of scene.

Our contributions in this subject are :

• Define an analytic and comprehensive formulation of light scattering along view rays and illumination rays.

Based on a angular formulation of the radiative transfer equation, we present precomputed 2D tables to compute directly these contributions. • Integration of the volumetric shadows. We build a method based on the shadow volume technique and using spatial coherence strategy, allowing the rendering of volumetric shadows in the scene, especially discernable around light sources. • Hardware implementation. Except for the determination of object's silhouette, all the work is done with GPU. We store our precomputed 2D tables in textures and use simple shaders to render the illumination of objects and participating media. The rendering of volumetric shadows also involves an intensive use of the stencil buffer. Our method is not restricted to isotropic participating media which can be lit by one or several, static or moving, point light sources since no precomputation are done involving either lights or camera. Our technique produces high resolution images and takes into account volumetric shadows, cast by occluders contained in the media. With very few texture memory cost, but using intensively graphics hardware, our method can render images at a high frame rate and is realtime for classical graphics scene. Our method is also easy to implement in traditional graphics engines since it follows the same strategy than the shadow volume algorithm, and use only shaders and textures. Therefore, it is straightforward with our method to obtain animations where objects or even light sources can move.

II. THE SINGLE SCATTERING ILLUMINATION MODEL

In this section, we will present shortly the analytic formulation of the single scattering illumination model. The precise theory can be found in [START_REF] Biri | Real Time Single Scattering Effects[END_REF].

A. Single scattering along a ray

Considering a view ray immersed in a participating medium, the radiance L observed along u P from a point P can be written (see figure 2) :

L u P = L dss ( u P ) + e -ktOP L(P ) + L a e -ktOP + L oiss ( u P )
(1) where L dss is the direct single scattering and L(P ) the radiance of P attenuated by both absorption and out-scattering. The third term represents both emission and multiple scattering inducing a constant radiance L a along the ray. The last term is the object indirect single scattering.

The radiance of P that can be written :

L(P ) = L iss (P ) + e -ktSP L d (P ) + L a e -ktSP (2)
L iss is the indirect single scattering received on point P and L d (P ) an usual direct illumination, like Phong model, that is attenuated by absorption and out-scattering. The third term is once again the contribution of emission and multiple scattering on the illumination ray. Only L dss and L iss remain unknown. They both are the contribution of the first light scattering along, respectively, a view ray and an illumination ray. To resume we give directly the expressions of these two terms. Developments can be found in [START_REF] Biri | Real Time Single Scattering Effects[END_REF].

If we denote :

Γ(λ, γ) = γ 0 p(α+ π 2 )e -ktλ sin(α)+1 cos(α) dα (3)
the direct single scattering could be written :

L dss (d, φ) = k 2 t ΩI S 4π e -k t t kth Γ(k t h, atan( d -t h )) -Γ(k t h, atan( -t h )) (4)
Note that the 2D function Γ is purely numerical and only depends on the shape of the phase function. It can therefore be precomputed and stored in a 2D table.

The new single scattering Phong model can be with this equation :

L iss (P ) = k 2 t ΩI S 4π k s k t l Γ P (k t l) (5)
with

Γ P (φ ′ , λ)= 2π 0 π 0 π 2 φ ′ -π 2 p(α+ π 2 )e -λ cos(α-φ ′ )+1 cos(α) cos n (φ)dα (6)
We want to integrate now the effect of occlusions along any view rays. Equation (4) describes the particular case where the view ray remains totally lit. Using the laplace formula, it is straightforward to see that equation (4) becomes :

L dss (d, φ) = k 2 t ΩI S 4π e -ktt [Γ(k t h, γ P ) -Γ(k t h, γ B ) +Γ(k t h, γ A ) -Γ(k t h, γ O )] (7)

III. HARDWARE IMPLEMENTATION

A. Overview of our method Our algorithm is easy to implement. We present here the 5 steps of this method and we will precise for each step if the computation is done by CPU or by GPU.

1) (CPU) The silhouettes of every moving shadow caster are computed. If the light source is moving, every silhouette needs to be recomputed.

2) (GPU) Scene is rendered using the conventional polygonal rendering method to obtain the direct illumination and the indirect single scattering. Surface shadows can be obtained using shadow planes algorithms [START_REF] Heidman | Real Shadows Real Time[END_REF], [START_REF] Everitt | Practical and Robust Shadow Volumes[END_REF].

The stencil buffer now contains lit areas of the scene. An ambient fog is added to take into account both absorption and multiple scattering. 3) (GPU) Scene is rendered once more time and direct single scattering is computed for each vertex of the scene. Depth test is set to equality. Only lit parts of the scene are rendered thanks to the stencil buffer. 4) (CPU) Shadow planes determined by the object's silhouettes are sorted in a back to front order. 5) (GPU) Shadow planes are rendered in that precise order.

The depth test function accepts only planes that are closer to the camera. Front facing planes add their contribution when back facing planes subtract them. Stencil function is set to allow fragments if the stencil is equal to 1 for front facing planes and 0 for back facing ones. Front facing planes always decrement the stencil buffer and back facing ones always increment it. All stages have to be done for each light source. As in [START_REF] Everitt | Practical and Robust Shadow Volumes[END_REF], a initialization stage is done to obtain the ambiant lighting and the first depth map. Each stage is detailed in [START_REF] Biri | Real Time Single Scattering Effects[END_REF].

Fig. 2 .

 2 Fig. 2. Notations for our model

Fig. 1 .

 1 Fig. 1. The same scene lit a. (left) classically, b. (center) with single scattering and c. (right) with single scattering and volumetric shadows. Who could say precisely where are the lights in the left picture ? This is obvious in the right picture