• Breast

• 94% of cases are local or regional at diagnosis but 30% will relapse Pollard, N Eng J Med, 2016

• Estimate the amount of residual distant disease at diagnosis in order to personalize the adjuvant (chemo)-therapy

• Avoid unnecessary, heavy toxicities

• Lung

• 57% of cases are metastatic

• Decide whether to use whole brain radiation therapy or just (stereotactic) surgery

• Avoid cognitive impairment of the patient Steeg, Nat. Rev. Cancer, 2016

• Metastases are the main cause of death (>90%) from solid cancers Lambert and Weinberg, Cell, 2017 Institut Bergonié, Bordeaux

Metastasis (µ!"# = beyond, $"#$%& = place)

!"#$%&'()'# *+,"-.+/.*")%()%)'# ,'(("*'&%)'+& %&,.01+2)$

Growth rates of primary and secondary tumors gp and g Dissemination rate

d(Vp) = !Vp
Size distribution of the metastases "(t,v)

!"#$%#$#&' ()*+", -#.#$/ ,)0("* .1 0"#$%#$#&' '"//%2 -3.04"*#52 Model

_ [ Y _ ] _ [ ∂ t ρ(t, v)+∂ v (g(v)ρ(t, v))=0 g(V 0 )ρ(t, V 0 )=d(V p (t)) ρ(0,v)=ρ 0 M (t)= ⁄ +∞ V 0 vρ(t, v)dv = ⁄ t 0 d (V p (t -s)) V (s)ds 67$#$ "#8$/9:8;8<=".* >&./:8?@@@ dV p dt =(α p -β p ln(V p )) V p g(v) =(α -β ln(v)) v Iwata, Kawasaki
θ i = θ pop + η i , η i ∼ N (0, ω 2 ) 4&#
dVp dt =(αp -βp ln(Vp)) Vp      ∂tρ(t, v)+∂v ((α -β ln(v)) vρ(t, v))=0 g(V0)ρ(t, V0)=µVp(t)γ ρ(0,v)=0
TTR i = inf Ó t>0; N vis 1 t, θ i 2 > 1 Ô x i
• Implement clinical variables and biomarkers as biologically meaningful covariates

Model

Diagnosis

dVp dt =(αp -βp ln(Vp)) Vp      ∂tρ(t, v)+∂v ((α -β ln(v)) vρ(t, v))=0 g(V0)ρ(t, V0)=µVp(t)γ ρ(0,v)=0
Clinical database biomarkers at diagnosis ln

1 µ i 2 = ln (µ pop )+β T µ x i µ + η i µ , η i µ ∼ N (0, ω 2 µ ) ln 1 α i 2 = ln (α pop )+β T α x i α + η i α , η i α ∼ N (0, ω 2 α )
344%'5*$(*+/!"%2)2'26"%'1$(*5'")($/0$%0"'%%+% Semi-mechanistic mathematical model

-=0 = E = 1+ - ) = -) ) =0 =0 G = ) -G G =0 =0 H = G -H H =0 =0 =+ ) + G + H E 0

  Median model primary tumor 10th and 90th percentiles model primary tumor Data metastatic burden Median model metastatic burden 10th and 90th percentiles model metastatic burden !"#$%#&'( *$+,%-,..,/&0 )*'*%)*%+'$ ,"-&$./"(.%#*&(0 '#%,'$.1'(%'2%$%*3

  5&6(3 &*.'$7.892"):;.<'#+&(.=&);.>?@A ! 0"*,)12'3&4 /"(.BC.'#-.,&*)D.! " #$!% & " #$& 92") $'2 =")E&$$.B'(6.<'#+&(.F#)*%*G*& Benzekry et al. (Ebos), Cancer Res, 2016 Current work: implementing the mechanistic model into a statistical learning model • Use the model for predictions of time-to-relapse (TTR)

Features

  Growth law: g p (V p )=V p (! p -" p ln(V p ))Growth law: g(V)=V(!-"ln(V))Dissemination law: d(V p )= #(V p )

  What is the optimal time gap between administration of bevacizumab and cytotoxic chemotherapy? How to capture inter-individual variability for designing personalized therapies?• Possible explanation: transient normalization of the otherwise abnormal (leaky, tortuous) vascular architecture

Acknowledgements

Thanks for listening!

We are looking for postdocs! sebastien.benzekry@inria.fr


Aur Élien Grolet
email: aurelien.grolet@ecl2009.ec-lyon.fr

Fabrice Thouverez

VIBRATION ANALYSIS OF A NONLINEAR SYSTEM WITH CYCLIC SYMMETRY

Keywords: Geometrical nonlinearity, nonlinear normal modes, bifurcation, localization Density. ω, Ω Free frequency, Excitation frequency

This work is devoted to the study of non linear dynamics of structures with cyclic symmetry under geometrical nonlinearity using the harmonic balance method (HBM). In order to study the influence of the non-linearity due to large deflection of blades a simplified model has been developed. It leads to nonlinear differential equations of the second order, linearly coupled, in which the nonlinearity appears by cubic terms. Periodic solutions in both free and forced cases are sought by the HBM coupled with an arc length continuation and stability analysis.

In this study, a specific attention has been paid to the evaluation of nonlinear modes and to the influence of excitation on dynamic responses. Indeed, several cases of excitation have been analyzed: punctual one and tuned or detuned low engine order. The paper shows that for a localized, or sufficiently detuned, excitation, several solutions can coexist, some of them being represented by closed curves in the Frequency-Amplitude domain. Those different kinds of solution meet up when increasing the force amplitude, leading to forced nonlinear localization. As the closed curves are not tied with the basic nonlinear solution they are easily missed. They were calculated using a sequential continuation with the force amplitude as a parameter.

Introduction

The study considers both free and forced vibrations of structures with cyclic symmetry under geometrical nonlinearity. This class of systems appears in the model of bladed disc or space antenna [START_REF] Georgiades | Modal analysis of a nonlinear periodic structure with cyclic symmetry[END_REF], [START_REF] Vakakis | A multiple-scales analysis of nonlinear, localized modes in a cyclic periodic system[END_REF]. It leads to nonlinear differential equations of the second order, linearly coupled, in which the nonlinearity appears by cubic terms.

In the linear case, most of the natural frequencies appear in pair due to the perfect symmetry of the problem. These natural frequencies are related to deformed shapes with nodal diameters [START_REF] Samaranayake | Subharmonic oscillations in harmonically excited mechanical systems with cyclic symmetry[END_REF]. For weakly coupled and weakly mistuned systems, lo-calization can take place leading to motions that are confined on only a few substructures.

In the nonlinear case, the study of free vibrations relies on the definitions of nonlinear normal modes (NNMs) which have leads to many scientific papers [START_REF] Vakakis | Nonlinear normal mode and their application in vibration theory: An overview[END_REF] [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF]. Unlike linear system, the number of NNMs can exceed the number of dofs. The other NNMs arise from bifurcations. In both free and forced cases, if the ratio between the coupling of the substructure and the nonlinearity is small, Vakakis [START_REF] Vakakis | Normal mode and localiation in nonlinear systems[END_REF] showed that nonlinear localization can take place in a perfectly symmetric system. When this ratio increases, bifurcations occur and nonlinear localization disappear. Traveling waves motions have also been detected in such systems [START_REF] Vakakis | Dynamics of a nonlinear periodic structure with cyclic symmetry[END_REF], [START_REF] Georgiades | Modal analysis of a nonlinear periodic structure with cyclic symmetry[END_REF].

In recent works, Peeters [START_REF] Peeters | Nonlinear normal modes, part ii: toward a practical computation using continuation technique[END_REF] [1] used a shooting method coupled with a continuation algorithm to study the NNMs of a system with cyclic symmetry. Similar and non similar NNMs have been found. Moreover, he studied the modal interactions between modes. He showed that these interactions can occur even if the natural frequencies of the modes are not commensurable, and he detected a countable infinity of such interactions.

The aim of this paper is to study both free and forced nonlinear vibrations of a bladed disc, using the HBM coupled with an arc length continuation. This study emphases the numerous bifurcation that can happen in this kind of system. Attention has been paid to the localization phenomenon, and particularly to the link between non symmetric loading and localization. The effect of the force amplitude on the solutions is also studied.

Simplified model

Consider a structure with cyclic symmetry made of n identical substructures called sectors. As a result of its dimension, its materials, and because of the external effort, such sectors can experiment large deflections leading to geometrical nonlinearity.

Many papers have studied the effect of geometrical nonlinearity on thin structures. Among others, Benamar [START_REF] Benamar | The effect of large vibration amplitudes on the mode shapes and natural frequencies of thin elastic structures, part ii: fully clamped rectangular isotropic plates[END_REF] and Amabili [START_REF] Amabili | Theory and experiments for largeamplitude vibrations of rectangular plates with geometric imperfections[END_REF], focused on the case of thin rectangular plates with large deflection. The case of circular shells have been studied by Touzé [START_REF] Touze | Hardening/softening behaviour in non linear oscillation of structural systems using non linear normal modes[END_REF], and the case of beams by Lewandowski [START_REF] Lewandowski | Computational formulation for periodic vibration of geometrically nonlinear structures-part 1:theoretical background[END_REF] [START_REF] Lewandowski | Computational formulation for periodic vibration of geometrically nonlinear structurespart 2:numerical strategy and examples[END_REF].

Here the computation of our simplified mathematical model for a structure with cyclic symmetry is describe. In this paper, each sector is modeled by a thin rectangular plate clamped on one edge to a fixed frame. The coupling between the substructure is realized by a linear stiffness (fig. 1 a). First, only one sector will be considered in order to obtain the expression of the different kinds of energy. These preliminary results will then be applied to the system with cyclic symmetry.

Preliminary considerations for one sector

This section focuses on a single sector. The retained hypothesis will be presented and expressions for the different kinds of Consider a plate P with dimensions L x , L y , and thickness h in a Cartesian system of coordinates (O, x, y, z), clamped on the edge (x = 0). The displacement of a point with coordinates (x, y, z) in the direction (Ox) (respectively (Oy), (Oz)) is denoted by u (respectively v, w) (fig. 1 b).

Elastic strain energy of a plate

The Love-Kirchhoff hypothesis for the displacement are used so that the displacement u u u is expressed as:

u u u(x, y, z) =   u(x, y, z) v(x, y, z) w(x, y, z)   =   -z ∂w ∂x -z ∂w ∂y w(x, y)   (1)
In addition, the use of the Von Karman nonlinear straindisplacement relationships and the standard bidimensional Hooke's law [START_REF] Ribeiro | Non-linear free vibration of isotropic plates with internal resonance[END_REF], leads to the following expression for the elastic strain energy U:

U = 1 2 Eh (1 -ν 2 ) L x 0 L y /2 -L y /2 ( 1 2 ( ∂w ∂x ) 2 ) 2 + ( 1 2 ( ∂w ∂y ) 2 ) 2 + 2ν( 1 2 ( ∂w ∂x ) 2 1 2 ( ∂w ∂y ) 2 ) + 1 -ν 2 ( ∂w ∂x ∂w ∂y ) 2 dxdy + 1 2 Eh 3 12(1 -ν 2 ) L x 0 L y /2 -L y /2 ( ∂ 2 w ∂x 2 ) 2 + ( ∂ 2 w ∂y 2 ) 2 + 2ν ∂ 2 w ∂x 2 ∂ 2 w ∂y 2 + 1 -ν 2 (2 ∂ 2 w ∂x∂y ) 2 dxdy ( 2
)
where E is the Young's modulus and ν is the Poisson's ratio.

Kinetic energy and external load

By neglecting rotary inertia, the kinetic energy T of a rectangular plate is given by :

T = 1 2 ρh L x x=0 L y /2 y=-L y /2 ( u2 + v2 + ẇ2 )dxdy (3)
In this paper the terms u2 + v2 will be omitted since it is in O(h 2 ) which is small compared to ẇ2 which is in O [START_REF] Georgiades | Modal analysis of a nonlinear periodic structure with cyclic symmetry[END_REF]. The kinetic energy is therefore given by :

T = 1 2 ρh L x x=0 L y /2 y=-L y /2
ẇ2 dxdy (4)

The work W due to external forces is written as:

W = L x x=0 L y /2 y=-L y /2 ( f x u + f y v + f z w)dxdy ( 5
)
where f x , f y , f z represent the distributed forces per unit area acting respectively in (Ox), (Oy), (Oz) directions. In this study only a localized harmonic force orthogonal to the plate is considered, thus f x = f y = 0. The external load f z of a localized force is given by :

f z = A f cos(Ωt)δ(x -x f )δ(y -y f ) ( 6
)
where Ω is the excitation frequency, t is the time, δ is the Dirac function, A f is the force amplitude and (x f , y f ) = (L x , 0). In this case, equation ( 5) can be rewritten as:

W = w(L x , 0)A f cos(Ωt) (7)
2.2 Equation of motion for the structure with cyclic symmetry Let's come back to the structure with cyclic symmetry. For 1 ≤ j ≤ n, w j denotes the transverse displacement of the plate number j, U j its strain energy, T j its kinetic energy, and W j the work due to external forces on this plate. U j , T j and W j are obtain by substituting w by w j in equations ( 2),( 4), [START_REF] Vakakis | Dynamics of a nonlinear periodic structure with cyclic symmetry[END_REF]. The coupling between substructures is modeled by massless linear stiffness of value k. These stiffness are positioned in a point (x r , y r ) = (L x /4, 0) for all plates. The energy V j of such stiffness between the plates j and j + 1 is given by:

V j = 1 2 k(w j (x r , y r ) -w j+1 (x r , y r )) 2 for 1 ≤ j ≤ n with convention j + 1 = 1 if j = n (8)
The total energies U t , T t , V t and W t are then given by the sum over the number of plates n of the different local energies U j , T j , V j and W j . The discretization of the transverse displacement w j in now introduced. This discretization is done by a Rayleigh-Ritz method [START_REF] Liew | pb2-rayleigh-ritz method for general plate analysis[END_REF], and it allows us to have a simplified model of the system. Transverse displacements w j are interpolated by the following expression:

w j (x, y,t) = N ∑ i=1 λ j i (t)Φ i (x, y) (9)
where (Φ i ) are kinematically admissible shape functions, λ j i is the contribution of the shape function Φ i in the displacement w j , and N is the number of shape functions used for the interpolation. This method leads to a discretized problem with nN unknowns (λ j i ) 1≤i≤N, 1≤ j≤n . Equations of motion are then given by Lagrange's equations (Appendix).

In this paper, we propose to study a system with six identical substructures (n = 6), of which displacements are interpolated by a single Ritz shape function (N = 1). This will leads to a problem with six degrees of freedom which will be used for numerical computation.

The retained shape function for the example is Φ = ( x L x ) 2 , it satisfies the boundary condition clamped at x = 0. By using Lagrange's equations we finally obtained the following equations of motion: Ẍ X X + KX X X + βX X X 3 = F F F(t) [START_REF] Amabili | Theory and experiments for largeamplitude vibrations of rectangular plates with geometric imperfections[END_REF] with X X X = (λ j ) 1≤ j≤6 , the convention X X X 3 = (x 3 j ) 1≤ j≤6 , and K a matrix defined by:

K =         α + 2c -c 0 0 0 -c -c α + 2c -c 0 0 0 0 -c α + 2c -c 0 0 0 0 -c α + 2c -c 0 0 0 0 -c α + 2c -c -c 0 0 0 -c α + 2c         (11)
where α, c, β, and F F F are given in appendix 3 Resolution method: Harmonic Balance Method

Principle of the HBM

The HBM is a widely spread method for solving non linear differential equations. One of the advantages of this method is that it can treat weakly or strongly non linear system in the same way. It can be used for problem with friction [START_REF] Laxalde | Qualitative analysis of forced responce of blisks with friction ring dampers[END_REF] [17], or geometrical nonlinearity [START_REF] Lewandowski | Computational formulation for periodic vibration of geometrically nonlinear structures-part 1:theoretical background[END_REF], [START_REF] Ribiero | Nonlinear vibration of plates by the hierarchical finite element and contination method[END_REF]. Consider the following system of non linear differential equations

M Ẍ X X + KX X X + K K K nl ( ( (X X X) ) ) = F F F(t) (12)
where X X X is a vector of unknowns of dimension n, M and K the mass and linear stiffness matrix of the system, K K K nl the non linear forces, and F F F the external forces. The HBM consists in finding of periodic solutions X X X of the form:

X X X(t) = A A A 0 + N h ∑ k=1 A A A k cos(kωt) + B B B k sin(kωt) ( 13
)
where N h is the number of retained harmonics. Substituting equation ( 13) in equation ( 12), and projecting the result on the trigonometric base (1, (cos(kωt), sin(kωt)) 1≤k≤N h ), one obtains a set of n(2N h + 1) algebraic nonlinear equations with n(2N h + 1) + 1 unknowns

(A A A k ) 0≤k≤N h ,(B B B k ) 1≤k≤N h and ω.
The most important parameter in this method is the number N h of retained harmonics. This number is not known a priori, so convergence studies are needed to ensure a good estimation of the solution. In one hand the higher N h is the better the solution. But in another hand, when N h is too high, computations can require many time and memory. However, in many cases, few harmonics are needed to ensure a good convergence leading to reasonable system size.

Numerical procedure: Arc length continuation

In a general way, the algebraic system obtains by the HBM can be reformulated as following:

G G G(Y Y Y , ω) = 0 ( 14
)
where Y Y Y is a vector of dimension n(2N + 1) containing the unknown

(A A A k ) 0≤k≤N and (B B B k ) 1≤k≤N ,
and G G G is a nonlinear function taking its values in a space of dimension n(2N + 1). This system is going to be solved by an arc length continuation method. Continuation methods consist in finding numerically a series of points (Y Y Y i , ω i ) 0≤i≤M which verifies the following criterion :

G G G(Y Y Y i , ω i ) ≤ ε ( 15
)
where ε is a small parameter which determines the accuracy of the solution. In an arc-length continuation, solutions are parameterized by the arc-length s, so that

(Y Y Y i , ω i ) = (Y Y Y (s i ), ω(s i )).
The method takes place in two steps : a predictor step which gives an estimation of the solution and a corrector step which correct the solution until the criterion ( 15) is achieve. In the study a tangent predictor and the Newton-Raphson corrector have been chosen. Further explanations can be found in [START_REF] Nayfey | Applied nonlinear dynamics[END_REF]. A sum up of the continuation algorithm is given in Fig. 2. This algorithm has been implemented in the Matlab software.

Stability of solutions and bifurcations

The stability of periodic solutions is determined by Floquet theory. It needs the computation of the monodromy matrix and of its eigenvalues. A detailed presentation of the method is given in [START_REF] Nayfey | Applied nonlinear dynamics[END_REF].

Here, the focus is on the bifurcation that the system can withstand. There is two way of detecting a bifurcation. First, by monitoring the eigenvalues of the monodromy matrix: if some eigenvalues go outside the unit circle then there is a bifurcation, the type of the bifurcation is determined by the way that the eigenvalues leave the unit circle. The second method to detect bifurcations is by monitoring the determinant |J y | of the jacobian matrix

[J y ] = ∂G G G
∂Y Y Y of system ( 14): when this determinant is zero at some point, then the matrix [J y ] is singular and this point is a bifurcation point. The type of bifurcation is determined by the range of the matrix [J y J ω ] (where J ω = ∂G G G ∂ω ). For a turning point, the range of [J y J ω ] is n, and there is no particular treatment to apply since the arc length continuation can handle turning points. For a branching point, the range of [J y J ω ] is at most n -1, so there is at least two tangent vectors T T T = [T T T y y y T ω ] t such that [J y J ω ]T T T = 0. There is two way of computing the bifurcated branches. The first method is to add a small perturbation to the function G G G(Y Y Y , ω) in the vicinity of the bifurcation point. Since branching points are structurally unstable, the perturbation breaks the bifurcation. Once another branch have been detected with the continuation algorithm, one remove the perturbation and carry on the continuation [START_REF] Nayfey | Applied nonlinear dynamics[END_REF]. The other method is to compute the eigenvectors Ψ i of the matrix [J y ] associated with the zero eigenvalue at the bifurcation point. The eigenvectors Ψ i indicate the directions to be followed [START_REF] Ribeiro | Non-linear free vibration of isotropic plates with internal resonance[END_REF].

Free vibrations of a system with six dofs

The first objective of this paper is to estimate the free solutions of system [START_REF] Amabili | Theory and experiments for largeamplitude vibrations of rectangular plates with geometric imperfections[END_REF]. The natural approach is to use the nonlinear normal modes (NNMs) which have lead to numerous publications. A survey about NNMs is given in [START_REF] Vakakis | Nonlinear normal mode and their application in vibration theory: An overview[END_REF] and [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF]. The most used definition of NNM is the one given by Rosenberg in [START_REF] Rosenberg | On nonlinear vibration of systems with many degrees of freedom[END_REF], in which he defines the NNMs by the vibrations in unison of a free and undamped system. Other approaches to define NNMs exist, they are based on geometrical arguments like the invariant manifold definition proposed by Shaw and Pierre [START_REF] Shaw | Normal modes for non-linear vibratory systems[END_REF], or based on the normal form theory [START_REF] Touze | Hardening/softening behaviour in non linear oscillation of structural systems using non linear normal modes[END_REF].

In this section, the underlying linear system of problem (10) will first be studied. This will give starting points for the study of the nonlinear system. Then, the search for the NNMs of system (10) by the HBM coupled with an arc length continuation will be carried out.

Numerical values for the physical parameters are introduced :

L x = 1, 5m, L y = 0, 3m, h = 0, 03m, E = 210GPa, ν = 0.3, k = 8.10 5 N.m -1 (16)
which correspond to the following values for the parameters of equation ( 10):

α = 8, 7662.10 3 s -2 , c = 148, 36s -2 , β = 4, 6752.10 7 m -2 .s -2 (17)
Modal analysis of the underlying linear system The underlying linear system of problem ( 10) is given by: Ẍ

X X + KX X X = F F F(t) ( 18
)
Linear normal modes (LNMs) are given by the resolution of the eigenvalue problem KΦ = ω 2 Φ. Matrix K have single and double eigenvalues, but double eigenvalues represent the majority [START_REF] Samaranayake | Subharmonic oscillations in harmonically excited mechanical systems with cyclic symmetry[END_REF]. With each double eigenvalue are associated two distinct deformed shapes which are not uniquely defined. The mode shapes are characterized by their number of nodal diameters p.

The eigenvalues and eigenvectors of matrix K are given by (the linear frequencies are given in rad.s -1 ):

ω 0 = √ α = 93.63 Φ 0 = [ 1, 1, 1, 1, 1, 1 ]
ω 1 = √ α + c = 94.42 Φ c 1 = [ 1, 1, 0, -1, -1, 0 ] Φ s 1 = [ 1, 1 2 , -1 2 , -1, -1 2 , 1 2 ] ω 2 = √ α + 3c = 95.98 Φ c 2 = [ 1, -1, 0, 1, -1, 0 ] Φ s 2 = [ 1, -1 2 , -1 2 , 1, -1 2 , -1 2 ] ω 3 = √ α + 4c = 96.75 Φ 3 = [ 1, -1, 1, -1, 1, -1 ] (19)
The deformed shapes of equation [START_REF] Nayfey | Applied nonlinear dynamics[END_REF] are depicted in Fig. 3 Nonlinear normal modes In this section, the nonlinear normal modes of system [START_REF] Amabili | Theory and experiments for largeamplitude vibrations of rectangular plates with geometric imperfections[END_REF] are computed with the algorithm described earlier. NNMs branches are computed by starting from the corresponding linear normal mode at low amplitude. These branches are termed the backbone curves, and they are represented in an Energy-Frequency plot as in [START_REF] Georgiades | Modal analysis of a nonlinear periodic structure with cyclic symmetry[END_REF]. The energy considered here is defined by:

E(X X X) = A 0 2 + ∑ k ( A k 2 + B k 2 ) ( 20
)
where A A A k and B B B k are the HBM coefficients defined in relation [START_REF] Lewandowski | Computational formulation for periodic vibration of geometrically nonlinear structurespart 2:numerical strategy and examples[END_REF]. Since the non linearity is odd and because there is no damping in system [START_REF] Amabili | Theory and experiments for largeamplitude vibrations of rectangular plates with geometric imperfections[END_REF], only odd harmonics and cosine terms are retained in the development of X X X (the last condition corresponds to a phase condition in which all initial velocities are set to zero). A convergence study with the number of harmonic on the zero nodal diameter mode showed that a good approximation can be obtained when retaining only the first and the third harmonic. This approximation holds for amplitudes of vibration up to the plate thickness h. We suppose that this remark holds for all NNMs. Therefore the NNMs of system [START_REF] Amabili | Theory and experiments for largeamplitude vibrations of rectangular plates with geometric imperfections[END_REF] are sought of the form:

X X X(t) = A A A 1 cos(ωt) +

A A A 3 cos(3ωt) (21)

Natural nonlinear normal modes

We call a natural NNM, a NNM which is directly born from a LNM at low magnitude of vibration. Since there are six LMNs, there will be six natural NNMs. They are represented in Fig. 4 and Fig. 5. The first noticeable feature in these figures is the Frequency-Energy (or amplitude) dependence of the NNMs. In addition to the dependence of their oscillation frequency, the NNMs may also have their modal shapes varying with the amplitude of vibration: these are nonsimilar NNMs. NNMs born from the LNMs shapes Φ c i (i = 0, 1, 2, or 3) defined in equation ( 19) are similar NNMs, and NNMs born from the LNM shapes Φ s i (i = 1, 2) are nonsimilar NNMs.

For double modes, the two backbones curves associated with the p nodal diameter mode (p = 1 or 2) are different (fig. 5). This is because the number of sector is of the form = 2m with m odd.

As we shall see, other NNM branches bifurcate from these backbone curves.

Bifurcated NNMs

This section focus on the NNM with 3 nodal diameters (NNM3D). By using the procedure described in section 3.3, several branching point bifurcations have been detected and bifurcated branches have been computed. They are It means that only one (strongly localized) or several (weakly localized) sectors have a non negligible amplitude of motion, the others remaining virtually motionless [START_REF] Vakakis | Normal mode and localiation in nonlinear systems[END_REF]. Such localization phenomenon was already observed in a linear mistuned system [START_REF] Yana | Vibration mechanism of a mistuned bladed-disk[END_REF], but here, and in general nonlinear system with cyclic symmetry, it occurs without structural disorder. The spatial confinement of the energy causes the responses of some sectors to be high and might lead to premature failure of the blades.

Traveling waves All the NNMs computed so far correspond to standing wave motions in the sense that the dofs vibrate in a synchronous way as in Rosenberg's definition of NNMs. Some papers [START_REF] Vakakis | Dynamics of a nonlinear periodic structure with cyclic symmetry[END_REF], [START_REF] Georgiades | Modal analysis of a nonlinear periodic structure with cyclic symmetry[END_REF] refer to traveling wave motions in which vibrations are not synchronous anymore. For Vakakis [START_REF] Vakakis | Dynamics of a nonlinear periodic structure with cyclic symmetry[END_REF], these traveling waves arise from a 1:1 internal resonance between the two modes associated with the same number of nodal diameter.

Because of the phase difference between the dofs, a different phase condition is considered for traveling wave NNMs computation. Sine terms are re-introduce in the HBM development of solution X X X in order to take into account the phase difference. Therefore X X X is now sought of the following form:

X X X(t) = A A A 1 cos(ωt) + B B B 1 sin(ωt) + A A A 2 cos(3ωt) + B B B 2 sin(3ωt) (22)
The starting point, at low amplitude of vibration, for the continuation algorithm is given by:

A A A 1 = ε cos(φ φ φ p ) A A A 2 = 0 0 0 B B B 1 = ε sin(-φ φ φ p ) B B B 2 = 0 0 0 ( 23
)
where p correspond to the number of nodal diameter (p = 1, or 2), ε is the amplitude of vibration (considered small) and φ φ φ p is a vector which coordinates ((φ φ φ p ) i ) 1≤i≤6 are defined by:

(φ φ φ p ) i = (i -1) 2pπ 6 (24)
These starting points enable the computation of traveling waves propagating in the anticlockwise direction. Their companions, propagating in the clockwise direction, can be obtain by replacing φ φ φ p by -φ φ φ p in equation 23.

The backbone curves of traveling waves are given in Fig. 8 (curves a and b). By monitoring the jacobian matrix as indicated in section 3.3, a branch that bifurcate from the traveling wave with 2 nodal diameters has been found (fig. 8 curve c). This bifurcated branch is quite noticeable because it can be seen as a weakly localized traveling wave since only three of the six sectors vibrate with non negligible amplitude. To illustrate the phase difference between coordinates, we plotted in fig. 9 the time series of the six blades for the 2 nodal diameters localized traveling wave motions.

Forced vibrations

Results from the previous section (free vibrations) are in agreement with the literature regarding structure with cyclic symmetry under geometric (cubic) nonlinearity [START_REF] Georgiades | Modal analysis of a nonlinear periodic structure with cyclic symmetry[END_REF], [START_REF] Vakakis | Dynamics of a nonlinear periodic structure with cyclic symmetry[END_REF] [START_REF] Vakakis | A multiple-scales analysis of nonlinear, localized modes in a cyclic periodic system[END_REF]. They highlight the complexity of the free dynamic of such structures. Obviously, this complexity is going to come back in forced vibrations. We here propose to study the impact of the different external loading on the forced response. Simulations are carried out for the structure with six dofs which correspond to system [START_REF] Amabili | Theory and experiments for largeamplitude vibrations of rectangular plates with geometric imperfections[END_REF]. Numerical parameters are the same as in the previous section. In order to have finite amplitude response, a damping term is added to system (10) which becomes:

Ẍ X X + δ Ẋ X X + KX X X + βX X X 3 = F F F(t) (25)
with δ a damping coefficient defined by δ = ω 0 /200. The right hand side of equation ( 25) is considered to be of the following form:

F F F(t) = A A A F cos(Ωt) ( 26
)
where A A A F stands for the shape and magnitude of the external force and Ω stands for the excitation frequency. Solutions are computed by the HBM. Stability of solutions is investigated using Floquet theory, and bifurcations are detected as described in section 3.3. Since the nonlinearity is odd, only odd harmonics are retained in the HBM approximation. A convergence study of the stability with the number of retained harmonics showed that keeping the first and the third harmonic is sufficient enough to ensure valid results on stability. Therefore, solutions of equation ( 25) are sought of the form:

X X X(t) = A A A 1 cos(Ωt) + B B B 1 sin(Ωt) + A A A 2 cos(3Ωt) + B B B 2 cos(3Ωt)
(27) Three kinds of excitationwill be considered. First, low engine order excitation, then localized excitation, and finally, detuned excitation. In all cases, a parametric study with the force amplitude as parameter has been carried out. Because of symmetry, only one (or four) of the six sectors are going to be represented.

Low engine order excitation with zero nodal diameters

In this section, the forced response of system ( 25) is studied for an excitation F(t) which takes the shape of the zero nodal diameters LNM. This means that F(t) is of the form:

F F F(t) = A A A F cos(Ωt) with A A A F = a f Φ Φ Φ (28)
where Φ Φ Φ is equal to Φ Φ Φ 0 0 0 of equation ( 19)

The effect of the force amplitude a f on the form and on the stability of solutions is studied. For low force amplitude (Fig. 10 a), the response of the nonlinear system is close to the one of the underlying linear system, and the solution is stable for all frequencies Ω. By increasing the force amplitude (Fig. 10 b), an unstable zone is generated between two turning points. When the force amplitude is increased again (Fig. 11), a second zone of instability is generated after the second turning point. This second zone takes birth through a branching point bifurcation. By using the branch switching method described in section 3.3, a bifurcated branch of solution has been computed. The bifurcated branch is stable for a few points at the beginning of the curve and then an instability is generated through a Hopf bifurcation. The area near the peak of the bifurcated solution is stable and correspond to a weakly localized motion (framed areas of Fig. 11)

Localized excitation

In this section, the external loading acts on only one sector (the first one), therefore the force is given by:

F F F(t) = a f         1 0 0 0 0 0         cos(Ωt) (29)
The forced response has been computed for several values of the amplitude a f . Results are depicted from Fig. 12 to Fig. 15. Because of symmetry only the response for four of the six sectors has been represented.

For values of amplitude from a f = 0 to a f = a s f ≈ 0.8 4 , the nonlinear response is topologically equivalent to the response of the underlying linear system (fig. 12).

When a f is larger than a loc f ≈ 0.9125 4 , stable forced localization occurs (fig. 15): for a well chosen excitation frequency (between 98 and 100 rad.s -1 ), only the excited sector vibrate with a non negligible amplitude, and this kind of motion is stable. Moreover the part of the curve corresponding to forced localization is positioned around the backbone curve of the strongly localized NNM which have been computed in the previous section (free vibration). These results are in agreement with those of Vakakis presented in [START_REF] King | A very copmplicated structure of resonances in a nonlinear system with cyclic symmetry: non forced localization[END_REF].

When a s f ≤ a f ≤ a loc f (Fig. 14 and Fig. 13), it is possible to find another type of solution. These new solutions are represented by closed curve in the Amplitude-Frequency diagram. They are not tied up with the basic nonlinear response, but they are positioned around the backbone curve of the strongly localized NNM. Stability analysis shows that some parts of the closed curves can be stable for points that have very larger amplitude than the basic nonlinear solution.

Since the closed curves are not tied up with the basic nonlinear response, they are unreachable with a unique continuation scheme. To compute the closed curves, a sequential continuation with fixed frequency and with the force amplitude a f as a parameter has been used. The starting point was taken from the part of the curve that correspond to forced localization (for a f ≥ a loc f ), then the force amplitude was decreased by continuation until a f ≤ a loc f . The results obtained this way have been used as starting points for the arc length continuation. This method is illustrated in Fig. 16.

The new kind of solution makes us think that forced localization is the result of the fusion of the closed curve with the basic nonlinear solution when a f = a loc f . Whereas the basic nonlinear solution seems to be weakly nonlinear since it is close to the linear solution, there exist stable solutions (the closed curve) which have a large vibration amplitude and which could be attained depending on initial conditions. This remark is important because the closed curves could be easily missed by a classical continuation scheme, thus leading to a wrong design of the structure. Frequency

Detuned low engine order with zero nodal diameter excitation

It has been seen that in the case of low engine order excitation the system response is simpler than in the case of localized excitation. However we have to keep in mind that realistic excitations are not perfectly symmetric. So now, the external forces are considered detuned, that is a perturbation is added to the low Sequential continuation with the force amplitude as a parameter Figure 16. Illustration of the method used to computes closed curve solutions engine order excitation. The perturbation is assumed to act only on the force amplitude. The external forces are then given by:

F F F(t) = (Φ + ε) (30)
where Φ is a linear deformed shape and ε is a vector characterizing the perturbation Because of the projection in the HBM, adding a perturbation to the force is equivalent to adding a perturbation to the function G G G(Y Y Y , ω) defined in equation [START_REF] Ribeiro | Non-linear free vibration of isotropic plates with internal resonance[END_REF]. Then the branching point bifurcations which are structurally unstable will disappear [START_REF] Nayfey | Applied nonlinear dynamics[END_REF]. Thus, bifurcated branches will be part of the nonlinear solution and will be obtained by the continuation method.

Here, only a low engine order with zero nodal diameters is considered, but similar results could be obtained with the other vibration modes. Finally, the perturbation is assumed to act only on the first coordinate of the force. Therefore we have:

F F F(t) = a f         1 + ε 1 1 1 1 1 1         cos(Ωt) (31)
where ε 1 correspond to the detuning percentage of the force with zero nodal diameter. In the remainder of the paper the force amplitude a f will be set to a f = 1 4 . Simulations start with ε = 1%. Results are represented in Fig. 17. As expected the bifurcated branch of the case ε = 0% (Fig. 11) is now a part of the solution.

When ε 1 becomes larger than 26% it is possible to compute a secondary solution which correspond to a closed curve. This kind of solution is represented in Fig. 18 for ε 1 = 40%. This closed curve have stable parts, and it merges to the basic nonlinear solution for ε 1 = 70% (Fig. 19).

In this example the closed curve have been detected when the detuning is larger than 26% which is a quite significant detuning level. However, the threshold for the detection of closed curve is highly dependent on the damping of the system: for a smaller damping the closed curves would have appeared sooner.

Conclusion

In this study, the free and forced response of a system with cyclic symmetry under geometric nonlinearity has been computed. The discretized equation of has been obtained by a Rayleigh-Ritz method. Solutions were computed by the harmonic balance method, stability was investigated using Floquet theory and bifurcations were computed through a branch switching algorithm.

In the free case, in addition to natural nonlinear normal modes, other vibration modes were detected by studying the branching point bifurcation of the system. Some of the bifurcated solutions correspond to localized nonlinear modes. Traveling wave motions have also been detected in this study.

In the forced case, when the excitation is localized to one of the substructure or sufficiently detuned, in addition to the basic nonlinear response, another kind of solution has been detected. This secondary solution is represented by closed curve in the Amplitude Frequency diagram, and it merge with the basic nonlinear solution when increasing the force amplitude, leading to forced nonlinear localization. For low engine excitations, we computed several bifurcated branches which also correspond to stable localized motions.

This study is based on a simple model, therefore the results remains at a theoretical level. However, this paper highlights the complex dynamic of system with cyclic symmetry under geomet- rical nonlinearity. Of importance is the detection of closed curve solutions which can be stable, and which may play an important role during the design of such a mechanical system.

Appendix: Expression of parameters for the system with six dofs

Here, we give the definition of parameters of equation [START_REF] Amabili | Theory and experiments for largeamplitude vibrations of rectangular plates with geometric imperfections[END_REF].

With the definition of the Ritz shape function Φ given in section 2.2, the transverse displacement is given by: j (x, y,t) = λ j ( x L x ) 2 for 1 ≤ j ≤ 6 (32)

Lagrange's equations are given by:

d dt ∂T t ∂ λ j i + ∂U t ∂λ j i + ∂V t ∂λ j i = ∂W t ∂λ j i for 1 ≤ i ≤ N, 1 ≤ j ≤ n (33)
Then the different terms of Lagrange's equations (33) become: F j (t) = Φ(x j f , y j f ) f j (t) m = f j (t) m (35)

  (a) Simplified model of a bladed disc (b) Rectangular plate and system of coordinates

Figure 1 .

 1 Figure 1. Retained model and system of coordinates for a rectangular plate

Figure 2 .

 2 Figure 2. Algorithm for the computation of solutions by Arc-Length continuation

Figure 3 .

 3 Figure 3. Mode shapes of the linear system[START_REF] Ribiero | Nonlinear vibration of plates by the hierarchical finite element and contination method[END_REF]

Figure 4 .Figure 5 .

 45 Figure 4. Backbone curves and deformed shapes for natural NNMs with 0 (a) or 3 (b) nodal diameters

Figure 6 .Figure 7 .

 67 Figure 6. Backbone curve of NNM3D (-) and its localized bifurcations (-• -)(a: 1 sector, b: 2 sectors, c: 3 sectors, d: 4 sectors)

Figure 8 .

 8 Figure 8. Backbone curves of traveling wave motions with one nodal diameter (a), two nodal diameters (b), two nodal diameters localized on three sectors (c)

Figure 9 .

 9 Figure 9. Time series of the two nodal diameters localized traveling wave

a f = 1/ 8 Figure 10 .Figure 11 .

 81011 Figure 10. Nonlinear forced response for a low engine order with zero nodal diameter (-: stable, • • • : unstable, -• -: linear response)

Figure 12 .

 12 Figure 12. Nonlinear forced response for a localized excitation on sector 1 with a f = 2,5 16 (-: stable, • • • : unstable, -• -: linear response)

Figure 13 .Figure 14 .

 1314 Figure 13. Nonlinear forced response for a localized excitation on sector 1 with a f = 3,222 16 (-: stable, • • • : unstable, -• -: linear response)

Figure 15 .

 15 Figure 15. Nonlinear forced response for a localized excitation on sector 1 with a f = 1 4 (-: stable, • • • : unstable, -• -: linear response)

3 Figure 17 .Figure 18 .

 31718 Figure 17. Nonlinear forced response for a detuned excitation with ε 1 = 1% (-: stable, • • • : unstable, -• -: linear response)

Figure 19 .

 19 Figure 19. Nonlinear forced response for a detuned excitation with ε 1 = 70% (-: stable, • • • : unstable, -• -: linear response)

  d dt ∂T ∂ λ j = m λ j ∂U ∂λ j = mαλ j + mβ(λ j ) 3 ∂V ∂λ j = mc(2λ jλ j+1 -λ j-1 ) ∂W ∂λ j = mF j (t)


Marie Morvan

Emilie Devijver

Madison Giacofci

Valérie Monbet

PREDICTION OF THE NASH THROUGH PENALIZED MIXTURE OF LOGISTIC REGRESSION MODELS

Keywords: MSC 2010 subject classifications: 62H30, 62J02, 62P10 Mixture regression model, Prediction, Variable selection, Heterogeneous data, Spectrometry data

In this paper, a suitable and interpretable diagnosis statistical model is proposed to predict the Non-Alcoholic Steatohepatitis (NASH) from near infrared spectrometry data. In this disease, unknown patients profiles are expected to lead to different diagnosis. The model has then to take into account the heterogeneity of the data and the dimension of the spectrometric data. To this end, we propose to fit a mixture on the joint distribution of the diagnosis binary variable and the covariates selected in the spectra. Because of the high dimension of the data, a penalized maximum likelihood estimator is considered. In practice, a twofold penalty on both regression coefficients and covariance parameters is imposed. Automatic selection criteria such as the AIC and BIC are used to select the amount of shrinkage and the number of clusters. Performance of the overall procedure is evaluated through a simulation study and its application on the NASH data set is analysed. The model leads to higher prediction performance than competitive methods and provides highly interpretable results.

Introduction

Non-Alcoholic Fatty Liver Diseases (NAFLD) is nowadays one of the leading cause of liver disease in Western countries. NAFLD is characterized by a build-up of fat in the liver. Due to the worlwide increase of obesity and type 2 diabetes, its prevalence, currently estimated at 24%, is expected to further grow in the future (Younossi et al., 2018a). Combined with liver cell injuries and inflammation, subgroups of NAFLD patients may derive to Non-Alcoholic Steatohepatitis (NASH), a more serious form of NAFLD, which represents the second cause of liver transplants in the USA. While being crucial for patients and healthcare systems, diagnosis of NASH is still an issue as the disease is essentially asymptomatic with low specific syndromes. Properly detecting and staging the severity of NASH patients requires a liver biopsy with the disadvantages of being invasive, costly, source of potential surgical complications and subject to sampling and inter-observer variability. In addition, liver biopsy can not be envisaged at a large scale, nor be repeated in time. At the time being, although diagnosis of NASH based on noninvasive modalities is an active field of research, no prediction procedure achieved consensus in the medical community (see e.g. Younossi et al., 2018b). Mid-infrared spectroscopy provides a molecular fingerprint of biological fluid, such as blood sera or urine, and may represent a promising approach to predict and understand the physiological consequences of the disease. In this project we study a data set containing blood serum spectra of 395 morbidly obese patients [START_REF] Anty | A new composite model including metabolic syndrome, alanine aminotransferase and cytokeratin-18 for the diagnosis of non-alcoholic steatohepatitis in morbidly obese patients[END_REF], including 66 patients diagnosed as NASH. Diagnosis has been based on liver biopsy conducted within the Hepatology unit of the Nice University Hospital. From a statistical perspective, our goal is to propose a statistical learning model to assign a score to a patient spectrum. Basically, mid-infrared spectra represent the absorbance of a biological sample discretely measured on a given range of wavelength. Such data are complex to analyze as they raise several statistical issues.

1. Dimension. Each individual spectrum comprises hundreds of measurement points and more specifically, more variables than individuals in the sample. Hence the problem lies in a high-dimensional framework. 2. Inner nature of the data. Mid-infrared spectrum contains all the molecular composition of a biological sample. However, not every molecular group is expected to be associated with the disease. It is then required to use appropriate methods to select the range of wavelengths that are associated with the disease progression. 3. Inter-individual variability. The studied population is heterogeneous and yields high individual fluctuations. This may be attributed to external or metabolic factors that are directly linked to the pathology development.

Taking account of subject-specific variability is a central point in our approach. It is now commonly accepted that individual metabolisms may strongly differ depending on lifestyle, feeding or global medical path. Compelling a cohort of patients to fit a rigid model appears then to be naive. A popular approach is to decompose the cohort onto few reference profiles summarizing as much as possible metabolic behaviors. Usually referred to as disease trajectories in the literature [START_REF] Ross | Nonparametric mixture of gaussian processes with constraints[END_REF], such reference profiles provide experts and practitioners interpretable knowledge on patient conditions and valuable information to predict the diagnosis.

In supervised learning, Generalized Linear Mixed-effects Models (GLMM) are a flexible setting to structure variability across individuals through a grouping structure [START_REF] Breslow | Approximate inference in generalized linear mixed models[END_REF]. However, it requires the grouping structure to be known in advance which is a difficult task as the NASH disease causes remain poorly known. Mixture models become then the appropriate tool by modelling the unknown sub-population through a discrete latent variable. Such approaches have been widely studied in the regression context through the mixture of regression models paradigm [START_REF] Grün | Fitting finite mixtures of generalized linear regressions in R[END_REF][START_REF] Khalili | Variable selection in finite mixture of regression models[END_REF][START_REF] Städler | 1-penalization for mixture regression models[END_REF]. In dedicated research works the central idea is to model the conditional distribution of the response variable given the predictors as a mixture distribution.

The covariates are then treated as non-random variables and the heterogeneity is assume to be fully contained in the conditional distribution. Such a modelling is often unrealistic for observational data as predictors can display significantly different behaviours depending on cluster. Most importantly, those models are not appropriate in a prediction framework [START_REF] Hoshikawa | Mixture regression for observational data, with application to functional regression models[END_REF]: posterior cluster membership probabilities indeed depends on the unobserved response variable and cannot be computed. As a matter of fact, [START_REF] Grün | Fitting finite mixtures of generalized linear regressions in R[END_REF] examine model fitting strategies for finite mixture of generalized linear regressions but do not discuss the prediction issue. Several works examine alternative strategies to allow computation of posterior probabilities as in [START_REF] Misiti | Mixture of linear regression models for short term pm10 forecasting in haute normandie (france)[END_REF], [START_REF] Ahonen | Prediction with a flexible finite mixture-of-regressions[END_REF] or [START_REF] Bougeard | Clusterwise analysis for multiblock component methods[END_REF] but require to use a separate learning algorithm to do so. Concur-rently, the machine learning community developed mixture of experts (MoE) models. Such models are dedicated to prediction when dealing with heterogeneous regions in the input space. Based on the principle of divide-and-conqueer, MOE aim at estimating the distribution of the response variable conditionally on the covariates as for mixture of regressions but model the prior cluster size as functional mixing weights depending on the covariates. Classical mixing weights include exponential weights with the softmax gating network or more general ones [START_REF] Yuksel | Twenty years of mixture of experts[END_REF]. As pointed by [START_REF] Ahonen | Prediction with a flexible finite mixture-of-regressions[END_REF] MOE by strictly partitioning the covariates space can not be considered as proper mixture models. Hence, although appealing from a prediction perspective, estimated model parameters have no clear biological interpretation.

In this paper, we consider an in-between approach where the joint distribution of the predictors and the response is defined as a mixture. Thus, our model has the advantage to exploit grouping structure information carried out in the conditional distribution as well as in the predictors. In addition, prediction of the response variable can be straightforwardly computed as the posterior cluster membership probabilities do not depend on the unobserved response for a new observation. For the NASH disease data, we are dealing with binary response. Hence the model considered is a mixture of logistic regression models, where the mixture is defined for the joint distribution, and the covariates are assumed to be Gaussian. Inference is performed through the Expectation-Maximization (EM) algorithm which is adapted to the latent variable setting.

As previously mentioned, spectrometry data contain all the molecular information contained in a biological sample whereas only few wavelengths are expected to be informative regarding to the prediction of the disease. The regression coefficient vector is thus expected to be sparse and variable selection should be considered to achieve accurate parameter estimation of our model. We choose to consider an 1 -penalized likelihood approach to simultaneously achieve variable selection and parameter estimation. Such approach can be straightforwardly slot into the EM mechanism and benefit from theoretical guarantees in the mixture of regression framework [START_REF] Khalili | Variable selection in finite mixture of regression models[END_REF][START_REF] Städler | 1-penalization for mixture regression models[END_REF]. Similarly, GLasso estimator is considered for the precision matrix of the predictors in the clusters. Besides, this second penalization highlights dependence between covariates and reduces the dimension.

Our method leads to the estimation of patient profiles and the estimated parameters allow the interpretation of the molecular variables involved in the disease. Their interactions can be represented with graphical models using the estimated precision matrices and has brought valuable insights about the NASH disease. First, the mixture of logistic regression model is presented in Section 2.1, and the prediction step is described thereby. Then, parameter estimation by maximization of the likelihood is described in Section 2.2 and in Section 2.3, regularization is introduced to reduce the dimension and to allow interpretation. Model selection is discussed in Section 2.4 to select the number of clusters. Section 3 investigates the numerical performance in a simulation study. An R Markdown file available at http://rpubs.com/morvan_ma/PMLR allows to replay some of the analysis on simulated data. Finally, Section 4 illustrates our result on the data set concerning the NASH disease.

2. Penalized mixture of logistic regressions model with random design

Random design model specification and prediction

Let (X, Y ) ∈ X × {0, 1} be a pair of random variables with X ⊂ R p , for p ∈ N. The variable Y stands for a binary dependent variable whereas X is a random vector of p quantitative covariates. In a learning context, one seek to describe the conditional relationship f (y|X = x) between the response Y and the predictors X. We assume this relationship depends on a latent (unobserved) class variable Z = (Z 1 , . . . , Z k ), K ∈ N, which follows a multinomial distribution such that

Z ∼ M(1, π 1 , • • • , π K ),
where π k = P (Z k = 1) for k = 1, . . . , K with π k > 0 and K k=1 π k = 1. Conditionally on Z, we define a logistic regression model with random design hierarchically as

X|{Z k = 1} ∼ N p (µ k , Σ k ), (1)
Y |{X = x, Z k = 1} ∼ B p (k) (x) .
The notations N p , B denote respectively the Gaussian and binomial distributions. Parameters µ k ∈ R p and Σ k ∈ R p×p are respectively the mean and the covariance matrix of the random covariates X conditionally on {Z k = 1}. Moreover, given {Z k = 1}, the covariates X are related to the response variable Y through the logistic link function such as

logit p (k) (X) = X T β k ,
where logit : x → log(x/(1 -x)) and β k = (β k,1 , . . . , β k,p ) are the regression coefficients of the generalized linear model in cluster k. The marginal joint distribution of (X, Y ) (unconditional on Z) is then defined as a mixture of logistic regression with a random design. In particular, the joint density can be written as

(2) f Y,X (y, x) = K k=1 π k f Y |{X,Z} (y; β k )f X|Z (x; µ k , Σ k ).
In the sequel, the set of parameters for the mixture density (2) with K clusters will be denoted by Φ K = (π 1 , ..., π K , φ 1 , ..., φ K ), where φ k = (µ k , Σ k , β k ), k = 1, . . . , K, is the vector of parameters of the k-th cluster.

The originality of our model compared to traditional finite mixture of regressions model usually explored in the literature (see e.g [START_REF] Grün | Fitting finite mixtures of generalized linear regressions in R[END_REF][START_REF] Khalili | Variable selection in finite mixture of regression models[END_REF] is the random design framework. It has two important consequences. First, as stated in Model (1), the random covariates are distributed themselves as a mixture and hence are also impacted by the unknown grouping structure. As pointed out by [START_REF] Hoshikawa | Mixture regression for observational data, with application to functional regression models[END_REF], it is then possible to infer cluster membership probabilities conditionally on the values observed for the covariates. Besides, the random design framework allows to consider prediction of the response for out-of-sample observations. Regression based on a fixed design indeed suffer from an optimism excess as emphasized in [START_REF] Rosset | From fixed-x to random-x regression: Bias-variance decompositions, covariance penalties, and prediction error estimation[END_REF] and hence should be used for prediction only with in-sample observations.

Formally, for a given out-of-sample observation x 0 drawn from Model (1), the prediction rule is given by,

(3) E(Y 0 |X 0 = x 0 ) = K k=1 P(Y 0 = 1|Z 0k = 1, X 0 = x 0 )P(Z 0k = 1|X 0 = x 0 ),
with Z 0 = (Z 01 , ..., Z 0K ) being the latent random variable associated to the new observation X 0 .

Replacing both quantities with their expressions in the joint mixture of logistic regressions framework, one obtains a predicted value

Y 0 = E(Y 0 |X 0 = x 0 ) = K k=1 τ 0,k exp(x t 0 β k ) 1 + exp(x t 0 β k )
, where (4)

τ 0,k = P(Z 0k = 1|X 0 = x 0 ) = π k f N (x 0 ; µ k , Σ k ) K =1 π f N (x 0 ; µ , Σ ) , with f N ( • ; µ, Σ)
being the multivariate Gaussian density with mean µ and covariance Σ.

Model (1) can be compared to mixture of experts models (MoE), initially introduced in Jacobs et al. (1991) and reviewed in [START_REF] Yuksel | Twenty years of mixture of experts[END_REF]. Such class of models are also dedicated for prediction in the mixture context. The mixing proportions depend on the covariates through a link function which is chosen according to the application. The outcome variable distribution depends then on both the covariates X and on the latent cluster membership variable Z. In our setting, the cluster variable Z is assumed to follow a multinomial distribution as prior distribution, inducing posterior probabilities as weights depending on the covariates, and in that sense, can be viewed as some particular case of mixture of experts using Gaussian parametric forms in the gate [START_REF] Yuksel | Twenty years of mixture of experts[END_REF].

Estimation using penalized likelihood

Given a sample {(x i , y i ) i=1,...,n } of n independent realizations of the random variables (X, Y ), the unknown parameters Φ K = (π 1 , ..., π K , φ 1 , ..., φ K ) are estimated by maximum likelihood. Our work takes place in a moderate dimension context for the covariates X. However, even under this assumption, some variables may not be relevant for predicting the response variable Y and finding a consistent estimation of a full unstructured covariance matrix remains challenging. A penalized likelihood approach is therefore adopted with a double penalty. The first penalty is a Lasso penalty [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] dedicated to perform clusterwise selection of features relevant for the logistic regressions in the model. The second part of the penalization controls the estimation of the covariates covariance matrices involved in the clustering part of the model. The maximum likelihood estimator of covariance matrices is indeed known to poorly behave even in a moderate dimension. A Graphical Lasso penalty [START_REF] Friedman | Sparse Inverse Covariance Estimation with the Graphical Lasso[END_REF] is thus used to estimate covariance matrices, under the assumption that their inverses are sparse. Besides, this assumption seems reasonable for the real data studied in the sequel.

The penalized log-likelihood is finally defined, for λ k ≥ 0, ρ k ≥ 0, for all k = 1, . . . , K, as

(5) ln L(y 1 , . . . , y n , x 1 , . . . ,

x n ; Φ K ) - K k=1 λ k β k 1 - K k=1 ρ k Θ k 1 ,
where

β k 1 = p j=1 |β k,j |, Θ k = Σ -1
k is the precision matrix in the kth cluster and Θ k 1 denotes the sum of the absolute values of the elements of Θ k . The regularization constants λ k and ρ k drive the amount of shrinkage on the parameters β k and Θ k for every cluster k, k = 1, . . . , K. Their selection is described in Section 2.3.

Penalized maximum likelihood estimation consists in maximizing the convex function ( 5) with respect to parameters Φ K . In a latent variable framework, this objective is usually achieved through an Expectation-Maximisation (EM) algorithm [START_REF] Dempster | Maximum likelihood from Incomplete Data via the EM Algorithm[END_REF], which alternates between two steps, E-step and M-step. Convergence properties of the EM algorithm for penalized likelihood are studied in [START_REF] Green | On use of the em algorithm for penalized likelihood estimation[END_REF]. At iteration [h], the E-step of the algorithm resumes to compute, for all individuals i = 1, . . . , n, the cluster membership posterior probabilities τ

[h] ik given the current parameters value Φ [h-1] K
. The posterior probabilities are computed thanks to ( 6)

τ [h] ik = P(Z ik = 1|X i = x i , Y i = y i ) = π k f B (y i ; p β k (x i ))f N (x i ; µ k , Σ k ) K =1 π f B (y i ; p β (x i ))f N (x 0 ; µ , Σ )
,

where f B (•; p β k (x i )) denotes the density of the Bernoulli distribution with parameters p β k (x i ) = exp(x t i β k )/(1 + exp(x t i β k )).
Then, in the M-step, the conditional expectation of the penalized likelihood given the current parameters value is maximized. It leads to the following updates of every parameter at iteration [h], for all k = 1, . . . , K:

π [h] k = 1 n n i=1 τ [h] ik , µ [h] k = n i=1 τ [h] ik -1 n i=1 τ [h] ik x i , Θ[h] k = arg max Θ k n i=1 τ [h] ik log det Θ k - 1 2 x i -µ [h] k T Θ k x i -µ [h] k -ρ k Θ k 1 , β[h] k = arg max β k n i=1 τ [h] ik y i x T i β k -ln 1 + exp(x T i β k ) -λ k β k 1 .
Few iterations of a Newton-Raphson algorithm are run to compute

Θ[h] k and β[h]
k . Note that it is not necessary to let the optimization run until convergence since an improvement of the penalized likelihood at each iteration is enough to guarantee the convergence of the EM algorithm to a local maxima [START_REF] Meng | Maximum likelihood estimation via the ECM algorithm: A general framework[END_REF].

The convergence of the EM algorithm to the global maxima is known to highly depend on the initial parameters. Here, we adopt a Search/Run/Select (S/R/S) strategy as developed in [START_REF] Biernacki | Choosing starting values for the EM algorithm for getting the highest likelihood in multivariate Gaussian mixture models[END_REF] and we perform the initialization in 3 steps:

1. Find t initial positions of parameters: obtain a partition of the set of observations of the explanatory variables (x i ) i=1,...,n into K clusters with a k-means algorithm [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF]. According to this clustering, compute the logistic regression estimators in each cluster, for each of the t trials.

2. Run a small fixed number of iterations of the EM algorithm at the t initial positions previously found.

3. Among these t possible starting values, select the one which maximizes the loglikelihood to start the EM algorithm.

Tuning the regularization parameters

The tuning parameters λ = (λ 1 , ..., λ K ) and ρ = (ρ 1 , ..., ρ K ) are related to the amount of regularization, and their selection is a critical issue in a penalized likelihood approach. It is usually based on a trade-off between bias and variance: large values of tuning parameters tend to select a simple model whose parameters estimates have smaller variance, whereas small values of the tuning parameters lead to complex models, with smaller bias. Generalized crossvalidation is a popular method for the tuning parameters selection in penalized likelihood approaches (see for example [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF][START_REF] Khalili | Variable selection in finite mixture of regression models[END_REF]). However, this method is computationally heavy and [START_REF] Wang | Tuning parameter selectors for the smoothly clipped absolute deviation method[END_REF] show that Generalized Cross-Validation may lead to the selection of some irrelevant variables. Besides, different studies (for example [START_REF] Wang | Tuning parameter selectors for the smoothly clipped absolute deviation method[END_REF][START_REF] Khalili | Regularization in finite mixture of regression models with diverging number of parameters[END_REF][START_REF] Jiang | Model selection for the localized mixture of experts models[END_REF][START_REF] Lloyd-Jones | A globally convergent algorithm for lasso-penalized mixture of linear regression models[END_REF] suggest the use of the Bayesian Information Criterion (BIC) to that purpose. The BIC [START_REF] Schwarz | Estimating the dimension of a model[END_REF] is constructed to realise a trade-off between the number of free parameters and the quality of the model measured by the likelihood. For a given number of clusters K, the BIC is defined as

BIC (λ,ρ) = -2 ln L y 1 , . . . , y n , x 1 , . . . , x n ; Φ(λ,ρ) K + ν (λ,ρ) ln(n), with Φ(λ,ρ)
K the arguments of the maxima of the penalized log-likelihood function with tuning parameters λ and ρ. The quantity ν (λ,ρ) counts the number of free parameters of the model, corresponding to the number of non-zero coefficients of the model. As tuning parameters λ k and ρ k have to be chosen for each cluster k = 1, . . . , K, an automatic procedure is proposed. First a classification by MAP rule is derived after few iterations of the EM algorithm, and grids for λ and ρ are deduced for each cluster k. Parameters are then estimated for each combination of the possible values of the two dimensional grid and finally, the values of parameters minimizing the BIC are retained.

Selection of the number of clusters

The number of clusters K is a sensible parameter because it describes the heterogeneity of the population. In an unsupervised setting, K is unknown and thus has to be estimated as well. Besides variable selection, the BIC is also commonly used to determine the number of clusters K in a mixture models framework [START_REF] Keribin | Consistent estimation of the order of mixture models[END_REF]. For a given number of clusters K, the BIC is defined as

BIC K = -2 ln L(y 1 , . . . , y n , x 1 , . . . , x n ; Φ(λ,ρ) K ) + ν K ln(n), with Φ(λ,ρ) K
the maximum likelihood estimator restricted to the relevant variables and ν K the number of free parameters of the model estimated with K clusters. However, the BIC was designed for non-structured data, and [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF] show that in some misspecified situations, it can lead to a wrong choice of K. For that reason, the Integrated Classification Likelihood (ICL) criterion, better adapted to the model-based clustering framework, was developed (see [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF][START_REF] Mclachlan | Finite Mixture Models[END_REF]. An entropy term taking into account the concentration shape of the clusters is added to the usual BIC. For a model estimated with K clusters, the ICL is defined as

ICL K = BIC K -2 n i=1 K k=1
τ ik lnτ ik with τ ik being the posterior probabilities. In the case of well-separated clusters, the entropy term is close to zero, and the ICL criterion value is close to the BIC value. In case of non separated clusters, the entropy term is highly negative and the value of ICL K increases. In other words, the ICL criterion favors models with well separated clusters. Finally, for predictive purpose, the Akaike Information Criterion (AIC) is known to be more suitable [START_REF] Shmueli | To explain or to predict[END_REF]. With the previous notations, the AIC is defined as

AIC K = -2 ln L(y 1 , . . . , y n , x 1 , . . . , x n ; Φ(λ,ρ) K ) + 2ν K .
Those three criteria will be compared in the simulation study.

Simulation study

Performance of our procedure are assessed through a simulation study. In particular, the study focuses on the model selection properties and the ability of the model to correctly predict a diagnostic, i.e the response variable.

Synthetic data generation

Simulated data are randomly generated according to a mixture of K = 3 clusters with proportions π = (0.3, 0.3, 0.4). The dimension of the predictors is set to p = 40 variables with only 8 variables being relevant in each cluster. Four different cases described in details in Appendix 7 are investigated. As a brief summary, the first two cases (Case 1 and 2) consider the same support for the regression coefficients in every cluster, i.e. the non-zero coefficients are associated with the same variables, but the concentration in each cluster differs and the parameters of the regression models lead to unequal class balances with respect to the cluster. The clustering task in Case 1 can be considered to be easier than the one in Case 2 because the ratio of means difference over the within variance (also called signal to noise ratio) is larger. The two last cases (Case 3 and 4) consider different supports for the regression coefficients, i.e. the non-zero coefficients are associated with different variables depending on the cluster. In Case 4, an extra difficulty is further added: the cluster means µ 1 and µ 3 are assumed to be equal so that information about clusters is carried by the conditional distributions Y |X, Z, in other words, by the values of the regression parameters. For each model, the impact of the sample size is measured, n being respectively equal to 250 and 500. These values are chosen because they are plausible for the applications in mind. Each experiment is repeated 30 times and majority votes or confidence intervals are computed.

Model selection

The proposed penalized mixture of logistic regression model, referred to as PMLR, is compared to a finite mixture of logistic regression model introduced in Grün and Leisch ( 2007) and denoted MLR. Note that in MLR, the clustering structure is only driven by the conditional distribution of Y |X so that the posterior distribution of the latent class variable can not be defined and the model can not be used for any prediction task.

Number of clusters

Let us first consider the selection of the number of clusters. As stated earlier, several information criteria are compared. Table 1 reports the majority vote of each of them for the PMLR and the MLR models. For PMLR, the AIC leads, whatever of the size of the dataset, to the selection of 3 clusters which is the right number. For Case 4, where two of the three clusters only differ through the coefficients of the regression model, the criteria fails and selects models with only 2 clusters. The penalty terms of BIC and ICL criteria are stronger than the one of the AIC and they tend to under estimate the number of clusters. These trend is more severe for the MLR than for PLMR. Finally, we decide to use AIC for the selection of K.

3 3 1 1 2 2 1 1 2 2 1 1 Difficult case (2) 3 3 1 1 3 3 1 1 3 3 1 1 Different supports Easy case (3) 3 3 1 2 3 2 1 1 3 2 1 1 Difficult case (4) 2 2 1 2 1 2 1 1 1 1 1 1 Table 1
Selected number of clusters. We compare our method PMLR with the mixture of logistic regression MLR. For each method, we provide the number of clusters leading to the smallest value of the following criterion: AIC, BIC and ICL, as a majority vote over the 30 repetitions for each of the 4 simulation cases. The simulations are done for two sample sizes (n=250 and n=500). The true number of clusters, 3, is in bold.

The quality of the clustering can be further evaluated using the Ajusted Rank Index (see [START_REF] Vinh | Information theoretic measures for clusterings comparison: Variants, properties, normalization and correction for chance[END_REF]) which is abbreviated ARI and which measures the similarity between the simulated and the estimated clusterings. This index is bounded by 1 and the greater the better. The mean of ARI values computed for the models selected by the AIC are reported in Table 2. They are close to 1 for the PLRM models except for the difficult Case 4, where two clusters (among the three) are close. For the MLR model all the values are close to 0 and it is mainly due to the bad selection of the number of clusters. But it is also explained by the lack of information brought by the model about the cluster. Indeed, in MLR model, the clusters are only defined by the coefficients of the regression models while in PMLR the covariates are also involved in the clustering structure. 2 Performance in clustering via ARI. We compare PMLR with the mixture of logistic regressions MLR. For each method, for a number of clusters fixed to K = 3, we compute the Adjusted Rand Index: closest to 1, better the clustering. It is done over the 30 repetitions for each of the 4 simulation cases. The simulation are also done for two sample sizes, n=250 and n=500.

Method

Variable selection In order to evaluate the ability of the penalization method to select the relevant predictors, relevant variable detection (RVD) and irrelevant variable elimination (IVE) criteria1 are used. The RVD associated to an estimate βk = ( βk,1 , ..., βk,p ) is defined as RV D k = T P k /(T P k + F N k ) with T P k the number of coefficients (β k,j ) j=1,...,p correctly predicted as non zero and F N k the number of coefficients (β k,j ) j=1,...,p predicted zero while being non zero. The IVE is defined as

IV E k = T N k /(T N k + F P k ) with T N k the number of coefficients (β k,j ) j=1,.
..,p correctly predicted as zero and F P k the number of coefficients (β k,j ) j=1,...,p predicted non zero while being equal to zero. The RDV is equal to 1 if all the relevant variables are retained in the model, whereas the IVE is equal to 1 if all the irrelevant variables are eliminated from the model. Then, a selection is good is both criteria are close to 1. However, in practice, a trade-off between both is often observed. The distributions of RDV and IVE estimated on the 30 repetitions of the simulation-estimation tasks are represented by boxplots in Figure 1. The regularization procedure leads to high values of the RVD and low values of the IVE for the PMLR and the opposite occurs for the MLR. It means that, in the PMLR, the effect of the penalization is conservative in the sense that it tends to select models keeping to much predictors. It is usually a good strategy for the prediction, where models slightly too rich generally perform better. At the opposite, in the MLR, the regularization eliminates too much variables including some of the important ones. Distributions of Figure 1 also illustrate that RVD significantly decreases when the clusters are difficult to identify (as in Case 4).

Prediction performances

To evaluate the quality of the prediction of the binary response, we consider the PMLR model for a fixed number of clusters K = 3. For a sake of comparison, results are also reported for the one cluster model referred to as PLR for Penalized Logistic Regression and the logistic regression (LR) model which corresponds to the model with one cluster and no Figure 1: Performance in variable selection for β. We provide boxplots for the Relevant Variable Detection (RVD) on the left and boxplots for the Irrelevant Variable Elimination (IVE) on the right. We compare our method (PMLR) with the Penalized Logistic Regression (PLR), which selects relevant variables in the regression matrix with an 1 -penalty. Regularization parameters are selected with the BIC. Every case introduced in the simulation setting is described in abscissa.

penalization. For each of them, the Area Under the Receiver Operating Characteristic curve (AUROC) is computed from out of sample observations. AUROC is a common summary statistic for the goodness of a predictor in a binary classification task. It is equal to the probability that a predictor will rank a randomly chosen positive instance higher than a randomly chosen negative one. The distributions of the AUROC are plotted in Figure 2. In Case 1, the classes of the binary variable are easy to discriminate so that all the models behave well although the PMLR model performs slightly better. In simulation cases 2 and 3 the population is highly heterogeneous and the regression models differ depending on the clusters. In theses cases, the PMLR model, which takes into account the heterogeneity, is significantly better than the two other models without any mixture. In order to go into more details, Received Operating Characteristic curves are plotted for Cases 2 and 3. The curves corresponding to the 30 estimations-predictions of the PMLR model are almost all above the ones of the PLR and LR models highlighting better predictions. It is interesting to note that the curves of the PLR and LR models are not very different indicating that over-learning in the LR model does not deteriorate significantly the prediction. The identification of the heterogeneity is clearly more important than the selection of the relevant predictor for the prediction task. Finally in simulation Case 4, where two clusters have the same mean, all the models have similar poor performances. In general, the prediction becomes better according to the AUROC criteria if the model is easier to calibrate. This is illustrated by the decreasing of the median AUROC from simulation Case 1 to simulation Case 4.

The conclusions of the simulation study are that the proposed model selection procedures are efficient if the clusters are distinct enough. Furthermore, the introduction of the latent class variable is critical for the prediction of the response variable. Finally a sample size of 250 seems to be large enough to get reasonable estimations. Bias and variance of estimators of μ, Σ and β have been computed (not shown here). The bias is close to zero and the variances decreased with the sample size, as expected for likelihood estimates, and are already small enough for n = 250.

Application to the NASH data set

On biological applications, we usually face to individual effects changing the prediction rule. We assume here that there exist homogeneous clusters of observations, relying on biological or genetic similarities. Those similarities might be independent of the severity of the disease that we want to diagnose. A better prediction of the disease is achieved considering this cluster structure and more importantly, a better understanding of the disease evolution process is given by our modeling. First, we describe in more details the data set. Then, we describe the results get by our procedure PMLR, from estimation, interpretation and biological points of view.

The NASH data set

Non Alcoholic Steatohepatitis (NASH) is a disease affecting the liver, and characterized by a fat deposit in the liver cells. In the long term, this disease results in important perturbations of the patient's metabolism. Currently, the diagnosis is obtained after a liver biopsy and an histological study of the sample, which is an invasive method. The method we propose is based on spectrum measured on blood serum, thus non invasive, and leads to a prediction of the disease. Moreover, as the method is model-based, interpretation of the prediction results is possible, leading to a better understanding from experts of the disease and its evolution. Experts suggest that different unknown patient typologies exist, and for each typology the molecular signature to establish the diagnosis is different. The proposed model, based on the existence of a discrete latent class (cluster), takes into account this feature of our data.

As we deal with high-dimensional data (large number of wavelengths in each spectrum), experts suspect that among the information available, some variables are irrelevant for the NASH diagnosis. One of the aim is to select the relevant variables to improve the diagnosis and to allow a better interpretation of the data.

The data set we consider is the following. We observe 395 patients, including 66 NASH patients (∼ 17%), coming from Nice hospital, in France. Clinical variables and spectrometric measures on sera samples are available. The spectrometric curves represent a molecular fingerprint of the sample and reflect the metabolic profile of each patient, affected by the liver condition. Portions of the spectrometric curves are selected by experts for there ability to describe metabolism variations that could be linked to the liver condition of the patients. Spectral variables are used to construct the prediction model, whereas biological and clinical variables only help for the interpretation.

Analyses and results

Model selection

The data set is randomly split in a calibration set containing 4/5 of the individuals (316 individuals including 53 NASH patients) and a validation set containing the individuals left (79 individuals including 13 NASH patients). These sets are randomly chosen but contain the same proportion of NASH patients and no significant differences are observed between clinical variables within the two sets. The model is estimated on the calibration set, for K =1 to 3 clusters. The model selection criteria are evaluated for each model and represented in Table 3. The lowest AIC and BIC values are obtained for the model estimated with two clusters. The lowest ICL value is obtained for the model estimated with one cluster, but with a slight difference with the ICL value corresponding to the model with two clusters. Following the conclusions detailed in Section 3.2 we select the model with 2 clusters according to the AIC value. Remark that the model PMLR with one cluster is equivalent to PLR, so performance is summarized by PLR in Table 4.

K = 1 K = 2 K = 3 AIC -
In Table 4 Estimators and models Graphical models obtained from the sparsely estimated precision matrices for each cluster are represented in Figure 5. We observe different relationships between variables according to the cluster. Considering only the relationship between the variables, we can see that for the first cluster, a group of variables from X2 to X11 is densely connected. For the second cluster, we observe two groups of strongly linked variables: the first pool is concerned with the variables X2, X3, X4, X6, X7, X8, and the second with the variables X1, X5, X12, X14, X17, X19. The links between variables strongly differ according to the cluster. The node color represents the coefficient value of the variable for the model applied to the considered cluster. In the first cluster, the estimated model is very sparse, a large proportion of regression coefficients are close or equal to zero. In the second cluster, regression coefficients take more extreme values. The difference observed in the graphical structures and regression models for each cluster suggest different metabolic mechanisms of the patients, depending on the cluster.

Biological interpretation of the constructed model We characterize the clusters obtained with the selected model through the clinical variables available in Table 5. Val-

Figure 5: Graphical models. The selected model has 2 clusters, the network related to the first cluster is on the left and the network related to the second cluster is on the right. Precision matrices are sparsely estimated, so the resulting networks are sparse. Arrow colors correspond to the sign of the partial correlation (green for positive correlation, red for negative correlation) and the edge widths correspond to the value of the correlation (stronger is the color, stronger is the absolute correlation). Node color is associated with the value of the regression coefficient for each variable and cluster: strongest the color, larger is the value (green for positive values and red for negatives ones). No color indicates a value of zero for the regression coefficient (irrelevant variable). 5 Characterization of the clusters obtained with the clinical variables. The model selected has two clusters. Values correspond to the mean over individuals for each cluster. We compare the mean in each cluster with a t-test, reporting the p-value and the significance. For the nominal variable Sex, the women's rate is indicated, and a Fisher test is used for the comparison.

ues correspond to the mean over individuals for each cluster. We compare the mean in each cluster with a t-test, reporting the p-value and the significance. We observe significant difference between the two clusters for the variables ALT (corresponding to the alanine transaminase), AST.ALT (ratio aspartate aminotransferase-alanine transaminase), GGT (Gamma-glutamyltransferase), HBA1C (glycated hemoglobin), chol (cholesterol) and TG (triglyceride). In the first cluster, the variables associated with diabete (Gluc and HBA1C) have higher values as well as the variables indicating liver problems (ALT, GGT). More generally, patients from the first cluster seem to have more severe liver complications than patients from the second cluster according to the seric indicators. Moreover, there is no significant difference between the two clusters for the morphological variables (weight, height, BMI), so that model allows to recognize the severity of the liver injury even when patients are not different for morphological variables.

We also represent the distribution of the predicted scores according to the different stage of histological variables. We can see in Figure 6 that the predicted score is a good indicator of the histological caracteristics of the patient. Indeed, the score increases with the steatosis, ballooning and inflammation stage, indicators used to establish the diagnosis. Fibrosis does not enter in the NASH definition and thus is not predicted by our model.

Conclusion

In this paper we have presented a predictive method that allows to build a model on data structured in unknown clusters, including non-relevant variables. This method provides interpretable tools to help for a better understanding of the data, with similar or higher prediction performance than competitive predictive models. This work was suggested by a real data problem concerning the use of the spectrometric technology to develop a noninvasive diagnosis tool for the prediction of the NASH disease. We obtained very encouraging results both in terms of prediction performance and in terms of interpretation, with clusters characterized by clinical variables and a prediction score linked to histological variables. Moreover, it is common in medical problems to face structured data with unknown patients profiles. Thus, our method could be much broader used.

In this reported work, the analysis focused on a preselection of wavenumbers performed by experts, but an interesting direction for further research would be to consider the whole spectra. Thus, we would like to adapt this method to handle functional data and perform automatic selection of specific areas of the spectra. This would highlight the type of molecules involved in the disease and offer the possibility to link the different areas of the spectra. Moreover, the discriminant information allowing the best prediction could be held at the same time by the intensity values of the spectra at particular wavenumbers and by the shape of the spectra at specific areas of the spectra. A different projection scale could be considered according to the area of the spectrum selected.

Acknowledgements

Part of this work was supported by the CNRS and AMIES institutions through the exploratory projects PEPS I3A AppSpec and PEPS I3A STATOPO. The authors are grateful to Rodolphe Anty (Centre Hospitalier Universitaire de Nice) and the Hepatology unit for

Figure 2 :

 2 Figure 2: Prediction performance. We provide boxplots for Area Under the Receiver Operating Characteristic (AUROC) obtained for the 4 simulation cases for each method: Logistic Regression (LR), Penalized Logistic Regression (PLR) and Penalized Mixture of Logistic Regression with 3 clusters (PMLR-3). Simulations are repeated 30 times.

Figure 3 :

 3 Figure 3: Prediction performance. We provide Receiver Operating Characteristic curves (ROC) for each method: Penalized Mixture of Logistic Regression (PMLR) with 3 clusters, Penalized Logistic Regression (PLR), Logistic Regression (LR). Plot on the left corresponds to the case 2, whereas plot on the right corresponds to case 3. Simulations are replicated over 30 data sets, and we plot a curve per data set to highlight the variance.

Figure 4 :

 4 Figure 4: Performance in prediction. We provide boxplots of the scores with respect to the true class, for the Penalized Mixture of Logistic Regression with 3 clusters (PMLR). The red dashed line corresponds to the threshold automatically learned.

Figure 6 :

 6 Figure 6: Boxplots of the repartition of the score predicted with the selected model according to the stage of the histological variable considered. Top left: steatosis, top rigth: ballooning, bottom left: inflammation, bottom right: fibrosis.

  Model selection. Comparison of the model selection criteria AIC, BIC and ICL for models estimated by PMLR on the calibration set of the NASH data set for 1 to 3 clusters. The bold values indicate the best values of the criteria obtained.Statistical interpretation of the constructed model The proportions of each cluster are 0.66 and 0.34. The proportion of diseased patients changes according to the cluster: 19 % in cluster 1 and 12 % in cluster 2. Comparison of the prediction performance obtained with different methods: our method with 2 to 3 clusters (PMLR-2, PMLR-3), penalized logistic regression (PLR) and logistic regression (LR). The chosen model is PMLR-2. The bold values indicate the best values. The quantities we use for the comparison are the following: Area Under the Receiver Operating Characteristic (AUROC), sensibility (Se), specificity (Sp), negative predictive value (NPV), positive predictive value (PPV), classification rate (CR).

		50180	-50459 -30957
	BIC	-49936	-49970 -30627
	ICL	-49936	-49901	-30365
		Table 3
		PMLR-2 PMLR-3 PLR	LR
	AUROC	0.75	0.68	0.64 0.67
	Se	0.77	0.85	0.62 0.69
	Sp	0.76	0.5	0.62	0.7
	NPV	0.94	0.94	0.89 0.92
	PPV	0.38	0.25	0.24 0.31
	CR	0.76	0.56	0.62	0.7
		Table 4

  , we observe the highest AUROC values and good classification rate values for the model estimated with two clusters, that shows the lowest AIC and BIC values. The model selection with the AIC and BIC is consistent with the cross validation performance obtained. Compared to competing methods, the chosen model has the best performance with the highest AUROC (0.75) and good classification rate values (0.76), and a high negative predictive value indicating a good screening test. Distribution of the predicted scores according to the real class of the individuals from the validation set is represented by Figure 4.

These criteria correspond to the sensibility and specificity criteria in binary classifier evaluation, renamed here to avoid confusion.

Appendix: parameters definition

We summarize the parameters chosen to generate synthetic datasets in Table 6. Table 6 Summary of the parameters for the 4 cases. To define a model, β, µ and Σ have to be defined. As the mixture model has 3 clusters, there are 3 different parameters in each case (the cluster is represented with k). In Case 1 and Case 2, the relevant variables are the same within clusters. The differences are on the concentration of the clusters and the balance between the two classes in each cluster. In Case 3 and Case 4, the relevant variables are different within clusters. In Case 4, the clustering relies on Y |X whereas µ 1 and µ 3 are the same, where in Case 3, the clustering relies on Y |X and X. For the covariance matrices Σ, we have used diagonal matrices and banded matrices.

Address of the First, Third and Fourth author Institut de recherche mathématique de Rennes 263 avenue du Général Leclerc 35000 Rennes, France E-mail: marie.morvan@univ-rennes1.fr E-mail: joyce.giacofci@univ-rennes2.fr E-mail: valerie.monbet@univ-rennes1.fr

Address of the Second author Laboratoire d'Informatique de Grenoble -LIG Université Grenoble Alpes 700 Avenue centrale 38 400 Saint Martin d'Hères, France E-mail: emilie.devijver@univ-grenoble-alpes.fr