Anicet Beauvais
email: beauvais@cerege.fr

Nicolas J Bonnet

Dominique Chardon

Nicolas Arnaud

Mudlappa Jayananda

Very long-term stability of passive margin escarpment constrained by 40 Ar-39 Ar dating of K-Mn oxides

The post-rift denudation history of high-elevation divergent continental margins is central to decipher "source-to-sink" systems across such margins and to unravel the topographic evolution of their escarpment. We perform 40 Ar-39 Ar dating of supergene cryptomelane (K-Mn oxide) from supergene manganese ore deposits to constrain the age of in-situ formed laterites of the low-relief lowland and highland separated by the Western Ghats Escarpment (WGE) of Peninsular India. Documentation of laterites as old as 53 Ma on the highland and 47 Ma at the foot of the WGE shows that the escarpment stabilized before 47 Ma (possibly 60 Ma ago). The geomorphic setting of the dated weathering mantles further allows constraining post-Early Eocene denudation of the lowland and the highland to maximums of 5 and 15 m/My, respectively. These results allow refining apatite fission track thermochronology and cosmogenic radionuclides studies that overestimate 4 to 20 times denudation rates, particularly after escarpment stabilization.

INTRODUCTION

Constraining the denudation history of high elevation (i.e., escarpment bearing) divergent continental margins documents their post-rift topographic evolution that has implications to quantify their vertical motion, drainage pattern and sedimentary supply over geological time scales. Current evolution models of high-elevation passive margin topography predict contrasted denudation histories. Models based on apatite fission track thermochronology (AFTT) predict a short period (≤ 10-15 Ma) of high denudation focused on a coastal strip resulting in escarpment formation shortly after the onset of rifting (Brown et al., 2002;Braun and van der Beek, 2004; see also Matmon et al., 2002). An alternative model integrating AFTT, the elevation of paleosurfaces and indirect geological constraints argues for successive burials, uplifts and topographic rejuvenations after continental break-up implying that escarpments result from the last uplift episode (Japsen et al., 2012). However, possible thermal instabilities in the shallow crust (< 2 km, i.e., T°< 60°C) and the failure of AFTT models to quantify erosion of crustal sections < 1 km (Brown and Summerfield, 1997) result in large uncertainties in the timing and magnitude of denudation predicted by these models (see also Kohn et al., 2005). Whether and how escarpments are rejuvenated after their formation (e.g., Moore and Blenkinsop, 2006) is also debated because the mode and velocity of their retreat have crucial implications on the evolution of the drainage networks of continental margins and their sedimentary budgets (e.g., Tucker and Slingerland, 1994). Studies based on cosmogenic radionuclides (CRN) predict slow denudation and retreat of escarpments, i.e., 3 to 16 m/My [START_REF] Bierman | Slow rates of rock surface erosion and sediment production across the Namib desert and escarpment, southern Africa[END_REF], but are only valid over millennial time scales.

Once slow denudation regimes (< 20 m/My) install on continental margins, remnant landscapes may be commonly preserved as a result of limited to moderate relief inversion [START_REF] Beauvais | Modes, tempo, and spatial variability of Cenozoic cratonic denudation: The West African example: Geochem[END_REF]. In the tropical belt, such paleolandscape remnants are mantled by laterites that may host K-Mn oxydes (cryptomelane) datable by 40 Ar/ 39 Ar geochronology that already proved beneficial (e.g., Vasconcelos, 1999;Vasconcelos and Conroy, 2003;[START_REF] Beauvais | Chemical and physical erosion rhythms of the West African Cenozoic morphogenesis: The 39 Ar-40 Ar dating of supergene K-Mn oxides[END_REF]. If properly mapped and dated, lateritic paleolandscape remnants are very useful to calibrate erosion over geological time scales [START_REF] Beauvais | Modes, tempo, and spatial variability of Cenozoic cratonic denudation: The West African example: Geochem[END_REF]. Here we combine geomorphology with the dating of cryptomelane from in-situ formed lateritic Mn ore deposits preserved at the foot and above the Western Ghats Escarpment (WGE) of Peninsular India as a test of post break-up denudation scenarios deduced from AFTT and CRN studies (Gunnell et al., 2003;Mandal et al., 2015). Our results lead to very low denudation rates (< 5-15 m/My) since stabilization of the escarpment at least 50 Ma ago and argue for a great stability of both the topography and relief of the margin since then.

GEOMORPHOLOGICAL AND GEOCHRONOLOGICAL OUTLINES

The WGE is up to 1000 m high and was carved into the Deccan traps and their Archean basement (Fig. 1a). The escarpment separates the coastal lowland from the highland plateau and coincides with the continental divide along most of its trace (Fig. 1a). The lowland consists of a weakly dissected lateritic pediment (Figs. 1b and1c) making the seaward piedmont of the WGE (Widdowson and Gunnell, 1999), on which laterites have been loosely dated to Mid-and Late Tertiary by paleomagnetism (Schmidt et al., 1983). The poorly defined highland laterites would have Late Cretaceous-Eocene paleomagnetic ages (Schmidt et al., 1983). Those laterites are mostly bauxites and mantle a low-relief relict landscape best preserved as 1000-1100 m high mesas near the WGE (e.g., Krishna Rao et al. 1989). Within the highland plateau, relicts of that paleolandscape are sparser (Fig. 1b). At Sandur, 40 Ar-39 Ar dating of supergene cryptomelane formed in Mn ore deposits carried by such an highland relict (Figs. 1a and1b) documents a period of intense lateritic weathering at c. 36-26 Ma (Bonnet et al., 2014), after the formation of Eocene bauxites (Krishna Rao et al., 1989) and before Neogene landscape dissection (Radhakrishna, 1993).

Inversion of apatite fission track data predicts higher denudation in the lowland (up to 120 m/My) than in the highland (< 20 m/My) between 80 and 50 Ma, and low (< 20 m/My) denudation on either side of the escarpment after 40 Ma (Gunnell et al., 2003). This suggests that the WGE formed at c. 50 Ma, after the extrusion of Deccan traps (65 ± 2 Ma) and largely after rifting between India and Madagascar (88 ± 3 Ma).

FIELD RELATIONSHIPS, MATERIAL AND METHOD

The lowland pediment consists of a c. 40 km wide concave surface ranging in elevation from 30 to 300 m, and incised by 25 to 100 m deep valleys (Figs. 1b-c). Field observations indicate that the pediment has truncated a bauxitic profile before being in turn weathered and cemented by a ferricrete, which caps a 25 to 60 m thick lateritic weathering profile (Fig. 1c). The studied lowland Caurem and Naveli Mn ore deposits are lentoid pockets within the pediment weathering profiles developed from Archean manganiferous schists and phyllites. The two open cast pits are located at the foot of the WGE at ~ 100 m and 140 m elevation, respectively (Fig. 1c), and four samples were collected at altitudes of 68 m and 45 m in Caurem pit and 113 m and 85 m in Naveli pit (see supplemental material). The highland Sandur Mn ore deposit formed upon Archean manganiferous phyllites and is exposed on a relict lateritic paleosurface capped by a ferricrete at 1012-1015 m elevation (Fig. 1b). Two samples were collected on benches at altitudes of 890 to 975 m (see supplemental material). All samples are massive with botryoidal or cavity filling structures enabling cryptomelane crystallization.

Thorough optical microscopy and X-ray micro-fluorescence analyses of 200-300 m thin sections allowed separation of eight cryptomelane grains from 300-500 m symmetrical sections slabs using a binocular magnifier (see Bonnet et al., 2014). The separated grains were characterized using X-ray diffraction, electron microprobe analyses and scanning electron microscopy before irradiation. Cryptomelane (KxMn IV 8-xMn III xO16) crystallizes into a monoclinic prismatic system with a typical tunnel-type crystal lattice framed by a double chain of MnO6 octahedra and K + cations in the large tunnel to insure the electronic neutrality of the lattice (Turner and Buseck, 1979). High retentiveness of potassium (content up to 5.5 wt.% K) and radiogenic argon ( 40 Ar*) in the intra-crystalline tunnels warranties the suitability of cryptomelane for 40 Ar-39 Ar dating (see Vasconcelos, 1999).

Gas was extracted from the irradiated cryptomelane grains either from a step-wise heating procedure in a double staudacher-type furnace, or from step incremental heating of the grains with a CO2 laser probe power. The gas fractions were then cleaned and analyzed using a VG3600 or Argus IV mass spectrometer. The 40 Ar/ 39 Ar ages are calculated from plateaus encompassing at least three consecutive 39 Ar release steps comprising up to 50% of total 39 ArK released, and from best-fit inverse isochrones in 36 Ar/ 40 Ar vs. 39 Ar/ 40 Ar correlation diagrams (see also supplemental material).

RESULTS AND INTERPRETATIONS

The 39 Ar release spectra of irradiated cryptomelane grains are stacked in the Figures 2a and2b. Each age spectrum allows calculation of a precise age from a plateau encompassing heating steps overlapping at 2 confidence level while degassing at least 5% of the total 39 Ar released from the highly retentive intra-crystalline sites of cryptomelane structure. The calculated plateau ages also agree well with the isochron ages (Table 1; see also data repositories DR1 to DR2 in supplemental material). The 40 Ar-39 Ar ages from the highland Sandur deposit document a weathering period from c.

53 to c. 50 Ma (Fig. 2a and2c). (1989).

Lateritic profiles of the lowland pediment record weathering periods at c. 47-45 Ma, c.

24-19 Ma, and an episode at c. 9 Ma (Fig. 2b and2c). The 47-45 Ma ages date the early (bauxitic) weathering of the lowland, whereas the ages between 24 and 19 Ma document renewed weathering and formation of the pediment's capping ferricrete, the 9

Ma age reflecting a discrete weathering pulse.

DISCUSSION

The preservation of in-situ formed laterites as old as 53-50 Ma above and 47-45

Ma below the WGE indicates that the current topographic envelope of the SW Indian margin corresponds to a bauxitic paleosurface dating from the Early Eocene, which already included the WGE with its present amplitude (Fig. 3). The 47-45 Ma old laterites at 55 m depth underneath the lowland pediment belong to a bauxitic profile that did not exceed 120 m in thickness (e.g., Bardossy and Aleva, 1990). Even if the lowland bauxitic profile had been totally eroded (i.e., 120 m of stripping), the maximum denudation rate of the piedmont between early (bauxitic) weathering (47 Ma) and abandonment of the pediment (19 Ma) would be lower than 5 m/My. Given that incision of the pediment does not exceed 100 m, the incision rate of the WGE' piedmont is less than 6 m/My over the last 19 Ma. Anyhow, the remarkable preservation of the pediment surface argues for negligible net denudation of the piedmont's envelope after c. 20 Ma (Fig. 1c).

Different ages recorded at a same depth in the lowland weathering profiles indicate that preserved old lateritic mantles of the escarpment piedmont have undergone several weathering episodes under a slow mechanical denudation regime. Preservation of 47 Ma old bauxitic weathering mantles at least 60 m thick thus attests that the current escarpment established at least 47 Ma ago, and did not retreat since then. Therefore, escarpment formation or retreat driven by Neogene rejuvenation of the lowland (e.g., Radhakrishna, 1993;Widdowson and Gunnell, 1999) is excluded.

Topographic inversion of the low-relief bauxitic paleolandscape throughout the highland plateau rarely attains 450 m (e.g., Figs. 1b and3). Therefore, the ages (53-50

Ma) obtained on laterites mantling this relict landscape imply a post-Early Eocene incision rate of less that 9 m/My over the last 50 Ma. Given the very low preservation rate of the old paleolandscape (Fig. 3), this incision rate approximates the denudation rate [START_REF] Beauvais | Modes, tempo, and spatial variability of Cenozoic cratonic denudation: The West African example: Geochem[END_REF] Formation of 47 Ma-old lateritic bauxites at the foot of the WGE would be in agreement with the last denudation pulse predicted at c. 50 Ma in the lowland (Fig. 2c) from AFTT inversion model (Gunnell et al., 2003). But this also suggests that the escarpment is even older as bauxitic profiles form slowly at 3 to 10 m/My in contexts of tectonic quiescence and limited denudation (Boulangé et al., 1997). Therefore, if the lowland bauxitic profile comprising 47 Ma-old laterite was 120 m thick, it would have taken at least 12 My to develop under an escarpment that should have been stabilized at least 60 Ma ago.

Our results imply denudation rates lower than 5 m/My in the lowland for the last 47 Ma, as opposed to 5-20 m/My (between 40 and 0 Ma) and even 20-70 m/My (between 50 and 40 Ma) from the AFTT model (blue curve in Fig. 2c). Even erosion rates of c. 10 -2 to 10 -1 m/ky (10-100 m/My) derived from CRN studies in riverine sediments (Mandal et al., 2015) are also questionable given the remarkable preservation of thick weathering profiles as old as 47 Ma on the WGE' piedmont. Likewise, the preservation of c. 100 m thick highland Eocene lateritic profiles formed upon the latest basaltic flow (63 Ma) of the Deccan traps (Widdowson and Gunnell, 1999) would be unlikely under such erosion rates or those (15-25 m/My) derived from AFTT between 53 and 45 Ma (red curve in Fig. 2c). Furthermore, AFTT fails to precisely measure or detect less than 1 km of denudation, which corresponds to the height of most passive margin escarpments. In other words, our results suggest that AFTT-or CRN derived denudation rates are largely overestimated and not realistic once the escarpment there.

Given the early installation of very slow denudation regimes across high-elevation margins such as that of Peninsular India, we suggest that the combination of radiometric dating ( 40 Ar/ 39 Ar of K-Mn oxides, or 4 He/ 3 He of Fe-oxides) of lateritic landscape remnants with apatite (U-Th)/He thermochronological models be a powerful tool to refine the denudation history of divergent margins of the tropical belt.

3-22.

Chaubey, A.K., Gopala Rao, D., Srinivas, K., Ramprasad, T., Ramana, M.V., and Subrahmanyam, V., 2002, Analyses of multichannel seismic reflection, gravity and magnetic data along a regional profile across the central-western continental margin of India : Marine Geol., v. 182, p. 303-323. Gunnell, Y., Gallagher, K., Carter, A., Widdowson, M., and Hurford, A.J., 2003 Topographic cross-section from Sandur to the escarpment edge in the highland and Goa in the lowland (section trace on Fig. 1a). (c) Synthetic Cross-section of the lowland piedmont passing through the two dated Mn ore pits. Cryptomelane grains were separated from 300-500 m thick slabs by hand picking.

The separated grains were ultrasonically cleaned in absolute ethanol and conditioned in aluminium foil packets, to be irradiated for 50 hours in the TRIGA Mark-II reactor of Pavia University (Italia). The Factor J was determined from the analysis of the standard Taylor Creek Rhyolite sanidine-2 (TCRs-2) monitor, with an age of 28 ± 0.08 Ma (Baksi et al., 1996). The standard was analyzed after every ten unknown samples. After a two-month "cooling" period, the irradiated cryptomelane grains were loaded in a double vacuum Staudacher-type furnace for step heating Ar isotopes measurements. The furnace temperature was calibrated by means of a classical thermocouple, and the gas purification was accomplished using a cold trap with liquid air and Al-Zr AP10 getters (one hot, one cold) for 8 minutes before the introduction into the VG3600 mass spectrometer. One minute was allowed for equilibration before analysis. 40 Ar and 39 Ar were measured on a Faraday cup with a resistor of 10 11 ohm, while 39 Ar, 38 Ar, 37 Ar, and 36 Ar were analyzed using a scintillator and photomultiplier after interaction on a Daly plate. The analytical data are reported in data repositories (Figures DR1 and Table DR2), and the errors are quoted at the 1σ level. Plateau ages are calculated from at least three consecutive 39 Ar release steps comprising up to 50% of total 39 ArK released and overlapping at the 2σ confidence level (Fleck et al., 1977). Isochrone ages are accepted when mean square weighted deviation (MSWD) is less than 2.5 and the 40 Ar/ 36 Ar intercept within 2 from the ( 40 Ar/ 36 Ar)atm value of 298.56 ± 0.31 (Lee et al., 2006 ;Renne et al., 2009).

References cited

Baksi

TABLE DR2

Relative Abundances

36Ar [V] 39Ar [V] %1 %1 38Ar [V] %1 37Ar [V] Relative Abundances 36Ar [V] 39Ar [V] 40Ar [V] %1 38Ar [V] %1 37Ar [V] %1  (Ma) %1 40(r)/39(k)  K/Ca %1  (Ma) %1 40(r)/39(k)  K/Ca 40Ar [V] %1 CAU-1a

CONCLUSION

  40 Ar-39 Ar geochronology of supergene cryptomelane formed in situ in supergene manganese ore deposits on either side of the Western Ghats escarpment indicate the development and preservation of an Eocene weathering paleolandsurface mantling the escarpment, attesting of its stability since at least 47-50 Ma, and possibly since 60 Ma.The ages obtained also constrain slow denudation rates of the lowland (< 5 m/My) and the highland (< 15 m/My) of the escarpment after its stabilization. Our results allow refining erosion rates estimated from AFTT and show that absolute dating of lateritic paleolandscape markers may accurately constrain post break-up denudation history of passive margins, with important implications on their sediment delivery histories on geological time scales. Bonnet, N.J.,Beauvais, A., Arnaud, N., Chardon, D., and Jayananda, M., 2014, First 40 Ar/ 39 Ar dating of intense Late Palaeogene lateritic weathering in Peninsular India: Earth Planet.Sci. Lett., v. 386,. Boulangé, B., Ambrosi, J-P., and Nahon, D., 1997, Laterites and bauxites: In Soils & Sediments: Mineralogy and Geochemistry, H. Paquet & N. Clauer, Eds., Springer Berlin, Heidelberg, pp. 49-65. Braun, J., and van der Beek, P., 2004, Evolution of passive margin escarpments: what can we learn from low-temperature thermochronology?: J. Geophys. Res., v. 109, F04009. Brown, R.W., Summerfield, M.A., 1997, Some uncertainties in the derivation of rates of denudation from thermochronological data: Earth Surf. Proc. Landf., v. 22, 239-248. Brown, R.W., Summerfield, M.A., and Gleadow, A.J.W., 2002, Denudational history along a transect across the Drakensberg escarpment of southern Africa derived from apatite fission track thermochronology: J. Geophys. Research, v. 107, B12, 2350. Campanile, D., Nambiar, C. G., Bishop, P., Widdowson, M., and Brown, R., 2008, Sedimentation record in the Konkan-Kerala basin: implications fro the evolution of the Western Ghats and the Western Indian passive margin: Basin research, v. 20, p.

  Figures caption

Figure 2 .

 2 Figure 2. Stacked 40 Ar-39 Ar age spectra of cryptomelane from (a) the highland and (b)

Figure 3 .

 3 Figure 3. Cross-section of the Southwestern Indian divergent margin from offshore

Figure. DR1 .

 DR1 Figure. DR1.39 Ar releasing spectra showing well defined plateau ages with K/Ca (grey) and

,

  Denudation history of the continental margin of western peninsular India since the early Mesozoic -Reconciling apatite fission-track data with geomorphology: Earth Planet.Sci. Lett., v. 215, Japsen, P., Chalmers, J.A., Green, P.F., and Bonow, J.W., 2012, Elevated, passive continental margins: Not rift shoulders, but expressions of episodic, post-rift burial and exhumation: Global Planet.Change, v. 90-91, p. 73-86. Schmidt, P.W., Prasad, V., and Ramam, P.K., 1983, Magnetic ages of some Indian laterites: Palaeogeogr.Palaeoclimatol. Palaeoecol., v. 44,

	Tucker, G. E., and Slingerland, R. L., 1994, Erosional dynamics, flexural isostasy, and
	long-lived escarpments: A numerical modeling study: J. Geophys. Res. 99, B6,
	12229-12243.
	Turner, S., and Buseck, P.R., 1979, Manganese Oxide Tunnel Structures and Their
	Intergrowths: Science, v. 203, p. 456-458.
	Letters, v. 425, p. 154-167.
	Matmon, A., Bierman, P., and Enzel, Y., 2002, Pattern and tempo of great escarpment
	erosion: Geology, v. 30, p. 1135-1138.
	Radhakrishna, B.P., 1993, Neogene uplift and geomorphic rejunevation of the Indian

Kohn, B.P., Gleadow, A.J.W., Brown, R.W., Gallagher, K., Lorencak, M., and Noble, W.P., 2005, Visualizing thermotectonic and denudation histories using apatite fission track thermochronology: Rev.

Miner. Geochem., v. 58, p. 527-565. Krishna Rao, B., Satish, P.N., and Sethumadhav, M.S., 1989, Syngenetic and epigenetic

features and genesis of the bauxite-bearing laterite of Boknur-Navge, Belgaum district, Karnataka: J. Geol. Soc.

India, v. 34,

Lee J-Y, Marti, K., Severinghaus, J.P., Kawamura, K., Yoo, H-S., Lee, J.B., and Kim, J.S., 2006, A redetermination of the isotopic abundances of atmospheric Ar: Geochim. Cosmochim.

Acta, v. 70,[4507][4508][4509][4510][4511][4512]

Mandal, S.K., Lupker, M., Burg, J-P., Valla, P.G., Haghipour, N., and Christl M., 2015, Spatial variability of 10Be-derived erosion rates across the southern peninsular India escarpement: A key to landscape evolution across passive margins: Earth Planet. Sci. Peninsula: Curr. Sci., v. 64, no. 11/12, p. 787-793. Vasconcelos, P.M., 1999, K-Ar and 40 Ar/ 39 Ar geochronology of weathering processes:

Annu. Rev. Earth Planet. Sci., v. 27, p. 183-229. Vasconcelos, P.M., and Conroy, M., 2003

, Geochronology of weathering and landscape evolution, Dugald River valley, NW Queensland, Australia: Geochim. Cosmochim.

Acta, v. 67, no. 16, p. 2913-2930. Widdowson, M., and Gunnell, Y., 1999

, Tertiary palaeosurfaces and lateritization of the coastal lowlands of western peninsula India: in Palaeoweathering, Palaeosurfaces and Related Continental Deposits, Thiry, M., and Simon-Coinçon, R., Eds., International Association of

Sedimentologists, Sp. Pub., v. 27, p. 245-274. Whiting, B.M., Karner, G.D., and Driscoll, N.W., 1994, Flexural and stratigraphic

development of the West Indian continental margin: J. Geophys.

Research, v. 99, p. 13791-13811.

Table caption Table 1 .

 caption1 Plateau and isochrone ages obtained for the cryptomelane separated from

samples collected in highland Sandur and lowland Caurem (CAU) and Naveli (NAV) Mn ore pits. Depth indicates the sampling level below the surface. MSWD = Mean

TABLE 1 Table

 1

Beauvais et al., 2015, Very long term stability of passive margin escarpment constrained by 40 Ar-39 Ar dating of K-Mn oxides Supplementary material

	Manganiferous ore samples have been sampled in upland Sandur ore deposits (KPA-8:
	15° 0' 2.27''N/76° 32' 41.6''E; KMK-2: 14° 59' 48.34''N/76° 34' 40.34''E) and in lowland ore
	deposits, Naveli (NAV-3: 15° 07' 55.10''N/74° 09' 49.72''E; NAV-4: 15° 08' 07.12''N/74° 09'
	32.26''E), and Caurem (CAU-1: 15° 07' 2.19''N/74° 08' 39.66''E; CAU-3: 15° 07' 3.47''N/74°
	08' 38.97''E).

  , A.K.,Archibald, D.A., and Farrar, E., 1996, Intercalibration of 40 Ar 39 Ar dating standards: Chem. Geol., v. 129, p. 307-324. Fleck, R.J., Sutter, J.F., and Elliot, D.H., 1977, Interpretation of discordant 40 Ar/ 39 Ar agespectra of mesozoic tholeiites from Antarctica: Geochim. Cosmochim. Acta, v. 41, p. 15-32.

Lee J-Y, Marti, K., Severinghaus, J.P., Kawamura, K., Yoo, H.S.,

Lee, J.B., and Kim, J.S., 2006

, A redetermination of the isotopic abundances of atmospheric Ar: Geochim.

Cosmochim. Acta, v. 70, p. 4507-4512. Renne, P.R., Cassata, W.S., and Morgan, L.E., 2009, The isotopic composition of atmospheric argon and 40Ar/39Ar geochronology: time for a change:

Quat. Geochron., v. 4,

  The concentrations of 36 Ar,37 Ar,38 Ar, 39 Ar and 40 Ar with their respective 1 error are provided for each step heating. The amount of 40 Ar* (%), of 39 Ar released (% 39 Ar) and the K/Ca ratio (derived from 39 Ar/ 37 Ar) are also given. Finally, this table show ratios 40 Ar*/ 39 Ark used to determine the corresponding apparent ages, which are presented in the last column with their associated 2error. The different J-Factor values are also provided for each irradiated grains.

	Supplemental file_Figure.DR1-A Supplemental file_Figure.DR1-B Supplemental file_Figure.DR1-C
	KMK-2		J = 0.00402950 ± 0.00000363
															Age 	40Ar(r)	39Ar(k)
																(%)	(%)
	48_KMK-2_	°C	4	1.1839652	0.511	231.8581	72.426	0.3215167	6.386	4.95915	0.668	371.0028	0.252	8.52452 ± 6.07881	61.07 ± 42.82	11.01	1.46	0.0089 ± 0.0129
	49_KMK-2_	°C	4	0.5680986	0.691	134.9113	172.430	0.2579779	9.837	6.52910	0.550	201.6829	0.468	7.03545 ± 6.19151	50.55 ± 43.87	22.43	1.95	0.0205 ± 0.0707
	51_KMK-2_	°C	4	0.2762329	1.149	265.9626	60.978	0.2666950	9.407	13.01527	0.628	148.9358	0.627	7.00618 ± 2.17753	50.34 ± 15.43	60.32	3.89	0.0207 ± 0.0253
	52_KMK-2_	°C	4	0.0984086	3.023	303.8510	62.665	0.2356036	8.707	15.12151	0.301	107.7680	0.866	7.00719 ± 2.19461	50.35 ± 15.55	96.89	4.53	0.0211 ± 0.0264
	53_KMK-2_	°C	4	0.0963079	3.187	216.0023	93.560	0.2716634	7.769	14.61645	0.193	114.3836	0.819	7.05201 ± 2.34370	50.67 ± 16.61	91.08	4.49	0.0288 ± 0.0539
	55_KMK-2_	°C	4	0.0999899	2.979	29.7181	649.669	0.2428391	8.840	14.53952	0.185	121.8410	0.766	6.52232 ± 2.27418	46.91 ± 16.15	77.72	4.41	0.2101 ± 2.7294
	56_KMK-2_	°C	4	0.1380666	2.228	-169.60	90.823	0.6075632	4.609	44.83913	0.130	369.8846	0.260	6.83980 ± 0.57628	49.16 ± 4.09	84.95	13.95	-0.1140 ± 0.2071
	57_KMK-2_	°C	4	0.1912850	1.639	143.1068	85.817	1.3379733	1.752	100.95234	0.115	758.6071	0.123	6.92242 ± 0.20665	49.75 ± 1.46	94.03	31.30	0.3030 ± 0.5201
	4_KMK-2_	°C	4	0.0770691	11.247	294.5770	98.924	0.4479906	5.973	33.31179	0.128	219.7423	1.161	6.55654 ± 1.48708	47.15 ± 10.56	100.90	10.27	0.0483 ± 0.0956
	5_KMK-2_	°C	4	0.0584562	14.822	403.3393	72.429	0.2896831	11.745	21.71829	0.152	117.2863	2.176	6.26068 ± 2.35533	45.05 ± 16.74	114.37	6.51	0.0228 ± 0.0331
	6_KMK-2_	°C	4	0.0724425	12.081	617.6113	48.098	0.3092851	10.266	22.86000	0.140	114.7255	2.224	6.49998 ± 2.29767	46.75 ± 16.31	126.97	6.81	0.0156 ± 0.0150
	8_KMK-2_	10 °C	4	0.0912916	9.525	704.0786	40.752	0.3388404	7.878	21.13750	0.219	103.3896	2.469	6.60114 ± 2.41926	47.47 ± 17.17	131.69	6.27	0.0126 ± 0.0103
	9_KMK-2_	12 °C	4	0.0525406	16.497	573.1992	66.571	0.2170933	14.152	14.13494	0.357	76.3872	3.342	7.98376 ± 4.79143	57.26 ± 33.83	143.38	4.17	0.0103 ± 0.0137
			S	3.0041547	0.734	3748.6166	23.075	5.1447248	1.833	327.73499	0.060	2825.6367	0.223
	KPA-8		J = 0.0039438 ± 0.00000355
															Age 	40Ar(r)	39Ar(k)
																(%)	(%)
	31_KPA-8B_	°C		2.8089159	0.536	234.32	130.296	0.6593686	3.542	3.20068	0.820	843.8298	0.072	11.13103 ± 17.59558	77.69 ± 120.21	4.00	1.52	0.0056 ± 0.0145
	32_KPA-8B_	°C		2.6898971	0.550	195.28	150.750	0.6711098	4.116	7.48457	0.411	829.3414	0.068	7.10266 ± 7.13672	49.96 ± 49.51	6.15	3.61	0.0162 ± 0.0488
	33_KPA-8B_	°C	4	1.7738239	0.533	261.41	118.841	0.5120498	5.930	10.48252	0.495	574.4086	0.096	7.18749 ± 5.30643	50.55 ± 36.80	12.60	5.06	0.0169 ± 0.0402
	35_KPA-8B_	°C	4	0.4962284	0.761	97.74	380.046	0.2574942	9.157	13.08041	0.307	212.9395	0.259	5.85080 ± 4.97662	41.25 ± 34.69	34.98	6.40	0.0572 ± 0.4350
	36_KPA-8B_	°C	4	0.2897844	0.903	213.48	232.238	0.2457546	6.535	11.22189	0.282	142.6239	0.386	6.78306 ± 7.64231	47.74 ± 53.08	52.63	5.56	0.0223 ± 0.1035
	37_KPA-8B_	°C	4	0.3112543	0.847	130.12	250.994	0.2640977	12.342	12.42831	0.238	162.4586	0.337	6.60376 ± 4.51844	46.49 ± 31.41	50.14	6.20	0.0408 ± 0.2046
	39_KPA-8B_	°C	4	0.1868086	1.810	23.81	913.916	0.9692660	2.109	71.99564	0.134	593.5181	0.148	7.49733 ± 0.51736	52.69 ± 3.58	90.92	36.17	1.2999 ± 23.7605
	40_KPA-8B_	°C	4	0.0683805	3.439	-745.35	46.385	0.3885421	7.195	32.44277	0.200	364.3965	0.454	8.68685 ± 1.84452	60.91 ± 12.72	76.96	16.22	-0.0190 ± 0.0177
	41_KPA-8B_	°C		0.0430701	2.523	-813.85	69.654	0.1578973	15.361	12.94441	0.371	248.4135	0.063	12.55212 ± 7.34728	87.37 ± 49.93	66.94	6.66	-0.0072 ± 0.0100
	43_KPA-8B_	°C		0.0425698	2.407	-572.74	96.309	0.1136482	26.458	8.10058	0.637	181.0269	0.085	14.35041 ± 11.36444	99.55 ± 76.71	66.08	4.19	-0.0064 ± 0.0123
	44_KPA-8B_	°C		0.0562862	1.242	-483.32	65.192	0.1757371	11.013	9.27267	0.519	174.7342	0.091	12.14450 ± 5.62442	84.60 ± 38.28	66.89	4.84	-0.0086 ± 0.0112
	45_KPA-8B_	°C		0.0601309	1.322	-240.19	109.146	0.1056815	20.501	7.09495	0.579	172.0936	0.089	18.80735 ± 6.36101	129.39 ± 42.23	77.76	3.58	-0.0130 ± 0.0284
			S	8.8271500	0.273	-1699.29	77.343	4.5206470	1.934	199.74940	0.087	4499.7847	0.052
	Tables. DR2-DR3. Analytical results obtained for highland and lowland cryptomelane grains,

either from Laser energy (spectrometer Argus IV) or double vacuum Staudacher-type furnace temperature, T °C, (spectrometer VG 3600) for each irradiated crypromelane grain.

TABLE DR3 (continued)

 DR3

	Relative Abundances		36Ar [V]	%1s	37Ar [V]	%1s	38Ar [V]	%1s	39Ar [V]	%1s	40Ar [V]	%1s	40*/39(k) ± 2s	(Ma)	K/Ca	± 2s
	Step	Temperature 40Ar/39Ar 38Ar/39Ar 37Ar/39Ar 36Ar/39Ar	39Ar	%39Ar	%40Ar*	40Ar*/39K	Age ± s.d.
	NAV-4f	°C				(E-3)	(E-14)
	Spectro : VG					moles
	J= 0.002727
	1		233329.16	1.04	66.39745	1039.77	0	0	99.87	0		0 ± 0
	2		477.918	0.409	0.01935	1647.584	0.015	3.91	-1.88	-8.96		0 ± 0
	3		8.185	0.118	0.0005	22.722	0.016	8.2	17.72	1.45	7.121 ± 3.503
	4		1644181.3	1.04	66.39745 5898367.2	0	8.2	-6.01	-103237.8		0 ± 0
	5		8.153	0	0.00269	77.269	0.003	9.01	-180.3	-14.7		0 ± 0
	6		6.563	0.021	0.00012	7.719	0.067	26.42	64.94	4.26	20.848 ± 0.988
	7		6.456	0.008	0.00016	7.304	0.052	40.09	66.26	4.28	20.923 ± 1.073
	8		10.924	0.015	0.00011	22.598	0.073	59.13	38.69	4.23	20.676 ± 1.731
	9		6.353	0.014	0.00011	0.002	0.075	78.77	99.68	6.33	30.89 ± 0.089
	10		8.322	0.012	0.00021	0.003	0.039	89.07	99.75	8.3	40.385 ± 0.211
	11		26.545	0.032	0.00157	24.211	0.005	90.45	72.97	19.37	92.865 ± 18.47
	12	1000	31.473	0	0.00246	18.469	0.003	91.33	82.6	26	123.561 ± 8.969
	13	1200	6.321	0.002	0.00038	8.767	0.022	97.05	58.7	3.71	18.162 ± 2.383
	Relative Abundances 14 1450	8.936	36Ar [V] 0	%1s 0.00073	37Ar [V] 25.99	%1s 0.011	38Ar [V] 100	%1s 13.83	39Ar [V] 1.24	%1s 6.07 ± 6.384 40Ar [V]	%1s	40(r)/39(k) ± 2s	(Ma)	K/Ca	± 2s
	Step	Temperature 40Ar/39Ar 38Ar/39Ar 37Ar/39Ar 36Ar/39Ar	39Ar	%39Ar	%40Ar*	40Ar*/39K	Age ± s.d.
	CAU-3f	°C				(E-3)	(E-14)
	Spectro : VG					moles
	J= 0.002727
	1		5.946	0.211	0.00042	27.703	0.019	4.68	-38	-2.26		0 ± 0
	2		4.964	0.071	0.00056	0.009	0.014	8.2	99.54	4.94	24.146 ± 0.216
	3		6.494	0.026	0.00027	3.348	0.029	15.57	84.46	5.49	26.786 ± 1.404
	4		6.145	0.021	0.00017	6.785	0.047	27.41	67.05	4.12	20.159 ± 2.29
	5		5.065	0.02	0.00008	3.36	0.097	52.01	80	4.05	19.826 ± 1.234
	6		9.738	0.021	0.00009	19.76	0.088	74.3	39.84	3.88	18.985 ± 0.934
	7		5.171	0.02	0.00013	2.684	0.061	89.72	84.27	4.36	21.315 ± 0.494
	8		11.343	0.048	0.00039	9.709	0.02	94.84	74.53	8.45	41.119 ± 1.734
	9	1000	28.501	0	0.00398	61.908	0.002	95.34	35.75	10.19	49.439 ± 11.841
	10	1200	13.136	0.009	0.53894	10.137	0.015	99.15	77.29	10.16	49.29 ± 1.805
	11	1450	52.209	0	30714.419	59.822	0.003	100	99.96	-99.93		0 ± 0
	Relative Abundances		36Ar [V]	%1s	37Ar [V]	%1s	38Ar [V]	%1s	39Ar [V]	%1s	40Ar [V]	%1s	40(r)/39(k) ± 2s	(Ma)	K/Ca	± 2s

TABLE DR3 (continued)

 DR3

	released	Ma
	released	Ma

Acknowledgements

-This work was funded by IFCPAR project 5007-1, the IRD (UR 161) and the CNRS (INSU 2011-CT2). The French Ministry of Research granted N.J.B. with a three years PhD scholarship (ED251, AMU, OSU Pytheas). Two referees are thanked for their comments on an earlier version of the manuscript.


Thibaud Cochereau

Lucie Bailly

Laurent Orgéas

Nathalie Henrich Bernardoni

Philippe Chaffanjon

A micro-mechanical model of the vocal-fold upper layers

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The vocal folds are soft multi-layered laryngeal tissues, owning remarkable vibro-mechanical performances. Composed of collagen and elastin microfibrils' networks, the upper layers play a major role in the vocal-fold vibrations. However, the impact of these tissues' histological features on their mechanical behavior is still poorly known. This is mainly ascribed to their challenging experimental characterization: vocal folds together with their fibrous architectures are not easily observable in vivo; ex vivo mechanical tests are rare and complex to interpret [1]. Consequently, most of the vocal-fold mechanical models developed so far rely on phenomenological macroscopic approaches, roughly assuming a homogeneous vocal tissue with linear-elastic properties. Since 2010, a few authors have started to investigate and model the vocaltissue collagenous fibrous microstructure, opening a new insight into voice biomechanics [2]. Theoretical formulations at the fiber scale still need to be developed. Thereby, this study aims at: (i) proposing an idealized but relevant model of the fibrous architecture of the vocal-fold upper layers; (ii) building a mechanical model able to predict the layers' multiscale properties from the above idealized architectures; (iii) assessing its relevance by comparison with a reference tensile database.

Methods

Idealization of the layer fibrous architectures -

The vocal-fold upper layers can be seen as a 3D composite structure made of a soft isotropic hyperelastic matrix, reinforced by a periodic network of bundles of wavy collagen fibrils displaying preferred orientations. This fibrous architecture is idealized by the Representative Elementary Cell (REC) shown in Fig. 1(a). Most of the REC parameters can be extracted from existing histological data [1,2,4]: fibril diameter (𝑑 0 ≈ 0.1 µm) and length 𝑙 𝑓 , initial tortuosity (𝜉 0 ≈ 0.84; 0.92) and lattice length (𝑙 0 ≈ 𝜉 0 𝑙 𝑓 ≈ 20; 30 µm), collagen fibrils content (𝜙 ≈ 0.15), initial orientation of fibril bundles (𝜃 0 ≈ 67.5°) with respect to vocal-fold main direction. Multiscale model -To determine the macroscale mechanical response of the as-idealized upper layers under simple uniaxial tensile loading, a non-linear and anisotropic multiscale model is proposed. This theoretical model is written using the homogenization method for periodic discrete structures, as already applied to the context of vascular biomechanics [3]. For that purpose, attention is paid to build a relevant micromechanical model describing the unfolding and the tensile deformation of wavy collagen fibrils, accounting for their elastic properties [4,5].

Results and discussion

The prediction of collagen fibrils' tensile behavior (𝜎 𝑓 , 𝜆 𝑓 ) is plotted in Fig. 1(a) for two typical fibril tortuosities: 𝜉 0 ≈ 0.84; 0.92 [2]. Unfolding the fibrils conducts to typical J-shaped stress-elongation curves depending on 𝜉 0 . Once the fibril unfolded, curves show a linear evolution with a tangent modulus of 1 GPa [4].

Figure 1 :

 1 Figure 1: Stress-elongation curves for tensile tests along e y : (a) model predictions for the fibril (b) model predictions and experimental data for the tissue.The macroscale predictions of the multiscale model are given by the stress-elongation curves plotted in Fig.1(b) for the same 𝜉 0 values. The mechanical behavior of the vocal ligament under tension is also illustrated in Fig.1(b), as an experimental reference recently reported using DIC measurements[1]. As observed at the microscale in Fig.1(a), both macroscale curves exhibit a J-shaped non-linear mechanical behavior, characteristic of collagenous soft tissues. Fig.1(b) also demonstrates a fairly good quantitative agreement with the reference mechanical data, keeping in mind that all model parameters derived from histological evidences. Further used to simulate in vivo loadings, this model will help to understand the role of the fibrous networks in vocal-fold vibratory properties. It will also provide proper guidelines to design new composites mimicking vocal-fold biomechanics.

Acknowledgements

This work is supported by the Labex Tec21 (Investissements d'Avenir -grant agreement ANR-11-LABX-0030).


Matthieu Marbac

Christophe Biernacki

Vincent Vandewalle

Model-based clustering of Gaussian copulas for mixed data

Keywords: Clustering, Gaussian copula, Metropolis-within-Gibbs algorithm, Mixed data, Mixture models, Visualization

Clustering of mixed data is important yet challenging due to a shortage of conventional distributions for such data. In this paper, we propose a mixture model of Gaussian copulas for clustering mixed data. Indeed copulas, and Gaussian copulas in particular, are powerful tools for easily modeling the distribution of multivariate variables. This model clusters data sets with continuous, integer and ordinal variables (all having a cumulative distribution function) by considering the intra-component dependencies in a similar way to the Gaussian mixture. Indeed, each component of the Gaussian copula mixture produces a correlation coefficient for each pair of variables and its univariate margins follow standard distributions (Gaussian, Poisson and ordered multinomial) depending on the nature of the variable (continuous, integer or ordinal). As an interesting by-product, this model generalizes many well-known approaches and provides tools for visualization based on its parameters. The Bayesian inference is achieved with a Metropolis-within-Gibbs sampler. The numerical experiments, on simulated and real data, illustrate the benefits of the proposed model: flexible and meaningful parametrization combined with visualization features.

Introduction

In a probabilistic framework, clustering is often managed by modeling the distribution of the observed variables using finite mixture models of parametric distributions [START_REF] Mclachlan | Finite mixture models[END_REF]. A class is defined as the subset of the individuals arising from the same mixture component. The literature covering homogeneous data (composed of variables of the same type) is extensive and presents Gaussian mixture models [START_REF] Banfield | Model-based Gaussian and non-Gaussian clustering[END_REF], multinomial mixture models [START_REF] Goodman | Exploratory latent structure analysis using both identifiable and unidentifiable models[END_REF] and Poisson mixture models [START_REF] Karlis | Exact Bayesian modeling for bivariate Poisson data and extensions[END_REF] as the standard models used to cluster such data sets. The use of conventional distributions for mixture components explains the success of these models, since the components can be easily interpreted. Although many data sets contain mixed data (variables of different types), few mixture models can manage these data [START_REF] Hunt | Clustering mixed data[END_REF] due to the shortage of multivariate distributions.

The locally independent mixture model [START_REF] Moustaki | Latent class models for mixed variables with applications in archaeometry[END_REF][START_REF] Lewis | Naive (Bayes) at forty: The independence assumption in information retrieval[END_REF][START_REF] Hand | Idiot's Bayes -Not So Stupid after All?[END_REF]) is a convenient approach for clustering mixed data since it assumes independence within-component between variables. Thus, each component is defined by a product of standard univariate distributions that facilitate their interpretation. However, this model can lead to severe bias when its main assumption is violated [START_REF] Van Hattum | Market Segmentation Using Brand Strategy Research: Bayesian Inference with Respect to Mixtures of Log-Linear Models[END_REF]. Therefore, two models have been introduced to relax this assumption.

The location mixture model [START_REF] Krzanowski | The location model for mixtures of categorical and continuous variables[END_REF][START_REF] Willse | Identifiable finite mixtures of location models for clustering mixed-mode data[END_REF] has been proposed for clustering a data set with continuous and categorical variables. It assumes that, for each component, the categorical variables follow a multinomial distribution and the continuous variables follow a multivariate Gaussian distribution conditionally on the categorical variables. Therefore, the intra-component dependencies are taken into account. However, the model requires too many parameters. [START_REF] Hunt | Theory & Methods: Mixture model clustering using the MULTIMIX program[END_REF] extended this approach by splitting the variables into within-component independent blocks. Each block contains no more than one categorical variable and follows a location model. The interpretation of this model can be complex since, for a given component, the univariate marginal of a continuous variable follows a Gaussian mixture model. Moreover, the estimation of the repartition of the variables into blocks is a difficult problem that the authors achieve with an ascending method that is sub-optimal.

The underlying variables mixture model [START_REF] Everitt | A finite mixture model for the clustering of mixed-mode data[END_REF] has been proposed for clustering a data set with continuous and ordinal variables. It assumes that each ordinal variable arises from a latent continuous variable and that all continuous variables (observed and unobserved) follow a Gaussian mixture model. The distribution of observed variables is obtained by integrating each Gaussian component into the subset of latent variables. However, in practice, this computation is not feasible when there are more than two ordinal variables. In an effort to study data sets with numerous binary variables, [START_REF] Morlini | A latent variables approach for clustering mixed binary and continuous variables within a Gaussian mixture model[END_REF] expanded this model by estimating the scores of latent variables from those of binary variables. However, the interpretation of the mixture components can be complex since it is based on the score-related parameters (not those related to observed variables).

Previous models illustrate the difficulty for clustering mixed data with a model for which interpretation and inference are easy. Moreover, they do not take account of cases where some variables are integer. The main difficulty is due to a shortage of conventional distributions for mixed data. However, copulas are standard tools for defining multivariate distributions in a systematic way, and they therefore have good potentiality for providing a sensible answer.

Copulas [START_REF] Joe | Multivariate models and dependence concepts[END_REF][START_REF] Nelsen | An introduction to copulas[END_REF] can be used to build a multivariate model by defining, on the one hand, the univariate marginal distributions, and, on the other, the dependency model. Recently, [START_REF] Smith | Estimation of copula models with discrete margins via Bayesian data augmentation[END_REF] and [START_REF] Murray | Bayesian Gaussian copula factor models for mixed data[END_REF] modeled the distribution of mixed variables using one Gaussian copula. As pointed out by [START_REF] Pitt | Efficient Bayesian inference for Gaussian copula regression models[END_REF], the maximum likelihood inference is very difficult for a Gaussian copula with discrete margins. Therefore, it is often replaced by the Inference Function for Margins method performing the inference in two steps [START_REF] Joe | Multivariate models and dependence concepts[END_REF][START_REF] Joe | Asymptotic efficiency of the two-stage estimation method for copula-based models[END_REF]. When all the variables are continuous, the fixed-point-based algorithm proposed by [START_REF] Song | Maximization by parts in likelihood inference[END_REF] achieves the maximum likelihood estimation, but this approach is not doable for mixed data. Therefore, as shown by [START_REF] Smith | Estimation of copula models with discrete margins via Bayesian data augmentation[END_REF], it is more convenient to work in a Bayesian framework since this simplifies the inference by using the latent structure of the model. This paper is organized as follows. Section 2 introduces the Gaussian copula mixture model and explains its links to well-known models. Section 3 presents the Metropolis-within-Gibbs algorithm used to perform Bayesian inference. Section 4 illustrates the algorithm behavior and the model robustness through numerical experiments. Section 5 presents two applications of the new model by clustering real data sets. Section 6 concludes the study.

2 Mixture model of Gaussian copulas 2.1 Finite mixture model An observation x = (x 1 , . . . , x e ) ∈ R c ×X is composed of e = c+d variables. Its first c elements, denoted by x c , correspond to the subset of the continuous variables defined on the space R c . Its last d elements, denoted by x d , correspond the subset of the discrete variables (integer, ordinal or binary) defined on the space X . Note that if x j is an ordinal variable with m j levels, then it uses a numeric coding {1, . . . , m j }. An observation is assumed to arise from the mixture model of g parametric distributions whose probability distribution function (pdf) is defined by

p(x|θ) = g k=1 π k p(x|α k ), (1)
where θ = (π, α) denotes the whole parameters. Vector π = (π 1 , . . . , π g ) is defined on the simplex of size g and groups the component proportions, where π k is the proportion of component k. Vector α = (α 1 , . . . , α g ) groups the component parameters, where α k denotes the parameters of component k.

Component modeled by a Gaussian copula

The Gaussian copula mixture model considers that each component follows a Gaussian copula.

Component k is also parametrized by α k = (Γ k , β k ) where β k = (β k1 , . . . , β ke ) groups the parameters of the univariate margin, β kj being the parameters of the j-th univariate margin, and where Γ k is the correlation matrix of size e × e. The cumulative distribution function (cdf) of component k is written as

P (x|α k ) = Φ e (Φ -1 1 (u 1 k ), . . . , Φ -1 1 (u e k )|0, Γ k ), (2)
where u j k = P (x j |β kj ) is the value of the cdf of the univariate marginal distribution of variable j for component k evaluated at x j , where Φ e (.|0, Γ k ) is the cdf of the e-variate centred Gaussian distribution with correlation matrix Γ k and where Φ -1 1 (.) is the inverse cumulative distribution function of the standard univariate Gaussian N 1 (0, 1).

The model implies the categorical variable z ∈ {1, . . . , g} which follows the multinomial distribution M g (π 1 , . . . , π g ) and which indicates the individual's component membership. In cluster analysis, the realization z is not observed while x is observed. Hence, mixture models are often interpreted as the marginal distribution of x based on the distribution of the variable pair (x, z). The Gaussian copula mixture model involves a second latent variable y = (y 1 , . . . , y e ) ∈ R e , such that y|z = k follows an e-variate centered Gaussian distribution N e (0, Γ k ). Thus, the Gaussian copula mixture can be interpreted as the marginal distribution of x based on the distribution of the variable triplet (x, y, z). Conditionally on (y, z = k), each element of x is defined by

x j = P -1 (Φ 1 (y j )|β kj ), ∀j = 1, . . . , e. (3)
Thus, the generative model of the Gaussian copula mixture is written as

• Class membership sampling: z ∼ M g (π 1 , . . . , π g ),
• Gaussian copula sampling: y|z = k ∼ N e (0, Γ k ),

• Observed data deterministic computation: x is obtained from (3).

Specific distributions for mixed-type variables

The cdf of component k defined by (2) implies the cdf of the univariate marginal distributions. Hence, it requires the definition of the distributions of the univariate margins (i.e. distribution of x j |z = k). We use conventional parametric distributions to facilitate the component interpretation. The parameters of margin j for component k are denoted by β kj . Hence,

• if x j is continuous then x j |z = k follows a Gaussian distribution with mean µ kj and variance σ 2 kj and β kj = (µ kj , σ kj ),

• if x j is integer then x j |z = k follows a Poisson distribution with parameter

β kj ∈ R + * ,
• if x j is ordinal then x j |z = k follows an ordered multinomial distribution with parameter β kj defined on the simplex of size m j . Note that the order between the levels is crucial since it permits the definition of the cdf.

Given that the first c variables of x (x c ) are continuous while the last d variables (x d ) are discrete, the pdf of component k can be decomposed as

p(x|α k ) = p(x c |α k ) × p(x d |x c , α k ). (4)
We use the decomposition into sub-matrices Γ k = Γ kcc Γ kcd Γ kdc Γ kdd , for instance Γ kcc is the submatrix of Γ k composed by the rows and the columns related to the observed continuous variables. Under component k, the knowledge of the continuous variable x j implies that y j =

x j -µ kj σ kj . Denoting y c = x j -µ kj σ kj ; j = 1, . . . , c , p(x c |α k ) = φc(y c |0,Γ kcc ) c j=1 σ kj where φ c (.|0, Γ kcc ) denotes the pdf of c-variate Gaussian distribution with mean 0 correlation matrix Γ kcc . If the variable j is discrete, any value y j in the interval S j k (x j ) =]b k (x j ), b ⊕ k (x j )] produces the same observation x j under component k, where b k (x j ) = Φ -1 1 (P (x j -1|β kj )) and b ⊕ k (x j ) = Φ -1 1 (P (x j |β kj )). Under component k, the distribution of the continuous latent variable y d = y j ; j = c + 1, . . . , e conditionally on y c is a Gaussian distribution with mean

µ d k = Γ kdc Γ -1 kcc Ψ(x c ; α k ) and covariance matrix Σ d k = Γ kdd -Γ kdc Γ -1 kcc Γ kcd .
Thus, the pdf of component k is written as

p(x|α k ) = φ c (y c |0, Γ kcc ) c j=1 σ kj × S k (x d ) φ d (u|µ d k , Σ d k )du, (5)
where

S k (x d ) = S c+1 k (x c+1 ) × . . . × S e k (x e
). Remark 2.1 (Model identifiability). The Gaussian copula mixture model is identifiable if at least one variable is continuous or integer (see Appendix A).

Strengths of the Gaussian copula mixture model

Related models

The Gaussian copula mixture model generalizes many conventional mixture models, including the four cases mentioned below.

• If the correlation matrices are diagonal (i.e. Γ k = I, ∀k = 1, . . . , g), then the model is equivalent to the locally independent mixture model.

• If all the variables are continuous (i.e. c = e and d = 0), then the model is equivalent to the Gaussian mixture model without constraints among parameters [START_REF] Banfield | Model-based Gaussian and non-Gaussian clustering[END_REF]. In the spirit of the homescedastic Gaussian mixture, we also propose a parsimonious version of the Gaussian copula mixture model by assuming equality between the correlation matrices over component. This model is named homoscedastic since the covariance matrices of the latent Gaussian variables are equal between components (i.e. Γ 1 = . . . = Γ g ). The free correlation model will be now called the heteroscedastic model).

• The model is linked to the binned Gaussian mixture model. For example, when variables are ordinal, it is equivalent to the mixture model presented by [START_REF] Gouget | Utilisation des modèles de mélange pour la classification automatique de données ordinales[END_REF]. In such cases, the model is stable through fusion of modalities.

• If the variables are both continuous and ordinal, then the model is a new parametrization of the model proposed by [START_REF] Everitt | A finite mixture model for the clustering of mixed-mode data[END_REF]. It should be noted that Everitt directly estimates the space S k (x d ) containing the antecedents of x d . Moreover, he uses a simplex algorithm to perform maximum likelihood inference, but this method dramatically limits the number of ordinal variables. The new parametrization of the proposed mixture allows the univariate marginal parameters β kj of each component to directly estimate (see details in Section 3), whereas Everitt's parametrization implies a difficult estimation of the bounds of integration. Thus, the parameter inference is easier.

Standardized coefficient of correlation per class

The Gaussian copula provides a user-friendly correlation coefficient for each pair of variables. Indeed, when both variables are continuous, it is equal to the upper boundary of the correlation coefficients obtained by monotonic transformation of the variables [START_REF] Klaassen | Efficient estimation in the bivariate normal copula model: normal margins are least favourable[END_REF]. Furthermore, when both variables are discrete, it is equal to the polychoric correlation coefficient [START_REF] Olsson | Maximum likelihood estimation of the polychoric correlation coefficient[END_REF].

Data visualization per component: a by-product of Gaussian copulas

By using the latent vectors of the Gaussian copulas y|z, a PCA-type method allows visualization of the individuals per component which permits the identification of main within-component dependencies. The visualization of component k is performed by computing the coordinates E[y|x, z = k; α k ] and then projecting them onto the PCA region associated with the Gaussian copula of component k. This space is obtained directly through spectral decomposition of Γ k . The individuals arising from component k follow a centered Gaussian distribution on this factorial map. Those arising from another component have an expectation not equal to zero. Therefore, an individual located far away from the origin arises from a distribution significantly different from the distribution of component k. Finally, the correlation circle summarizes the within-component correlations and avoids the direct interpretation of the correlation matrix Γ k , which can be tedious if e is large. The following example illustrates these properties.

Example 2.2. Let three variables-one continuous, one integer and one binary-arise, in this order, from the bi-component Gaussian copula mixture model parametrized by π = (0.5, 0.5), β 11 = (-2, 1), β 12 = 5, β 13 = β 23 = (0.5, 0.5), β 21 = (2, 1),

β 22 = 15, Γ 1 =   1 -0.4 0.4 -0.4 1 0.4 0.4 0.4 1   and Γ 2 =   1 0.8 0.1 0.8 1 0.1 0.1 0.1 1   .
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -4 -2 0 2 4 5 10 15 20 25

x 1

x 2 (a) Individuals described by three variables: one continuous (abscissa), one integer (ordinate) and one binary (symbol). Colors indicate the true class memberships q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -10 -8 a centred one (red) and a second one (black) located on the left side. More precisely, the first axis (explained by the continuous and the integer variables) is strongly discriminative while the second axis (explained exclusively by the binary variable) is not discriminative. Figure 1(c) shows the correlation circle of the first PCA-map of the red component. It allows a strong correlation to be identified, for the red component, between the continuous and the integer variables.

3 Bayesian inference

Sampling layout on data and parameters

We observe the sample x = (x 1 , . . . , x n ) composed of n independent realizations x i ∈ R c × X assumed to arise from a Gaussian copula mixture model. As pointed out by [START_REF] Smith | Estimation of copula models with discrete margins via Bayesian data augmentation[END_REF], the Bayesian framework simplifies the inference considerably since it uses the latent structure of the model (y, z). Without prior information about the data, we assume independence between the prior distributions. The proportions and the parameters of the univariate marginal distributions of each component β kj follow the classical conjugate prior distributions [START_REF] Robert | The Bayesian choice: from decision-theoretic foundations to computational implementation[END_REF]. Finally, the conjugate prior of the covariance matrices is derived from an Inverse Wishart distribution as proposed by [START_REF] Hoff | Extending the rank likelihood for semiparametric copula estimation[END_REF]. Details on the prior distributions are given in Appendix B.

Gibbs and Metropolis-within-Gibbs samplers

The Bayesian estimation is managed by a Gibbs sampler (described in Algorithm 3.1) which is the most popular approach for inferring mixture models since it uses the latent structure of the data. Its stationary distribution is p(θ, y, z|x) where z = (z 1 , . . . , z n ) denotes the class memberships of x and where y = (y 1 , . . . , y n ) denotes the Gaussian vector related to x. Note that the Gaussian variable y is twice sampled during one iteration of the algorithm to manage the strong dependencies between y and z, and between y j

[rk] = {y j i : z (r) i

= k} and β kj . Obviously, the stationary distribution stays unchanged. Thus, the sequence of parameters is sampled from the marginal posterior distribution p(θ|x), and a consistent estimate of θ can be obtained by taking the mean of the sampled parameters.

Algorithm 3.1 (The Gibbs sampler). Starting from an initial value θ (0) , its iteration (r) consists in the following four steps (k ∈ {1, . . . , g}, j ∈ {1, . . . , e})

z (r) , y (r-1/2) ∼ z, y|x, θ (r-1) (6) β (r) kj , y j(r) [rk] ∼ β kj , y j [rk] |x, y (r) [rk] , z (r) , β (r) k , Γ (r-1) k (7) π (r) ∼ π|z (r) (8) Γ (r) k ∼ Γ k |y (r) , z (r) (9)
where

y [rk] = y {i:z (r) i =k} , y (r) i = (y 1(r) i , . . . , y j-1(r) i , y j+1(r-1/2) i , . . . , y e(r-1/2) i ) and β (r) k = (β (r) k1 , . . . , β (r) kj-1 , β (r-1) kj+1 , . . . , β (r-1) ke ).
The samplings according to (8) and ( 9) are classical but the samplings from ( 6) and ( 7) are not easy. They are therefore replaced by one iteration of a Metropolis-Hastings algorithm that does not change the stationary distribution. The resulting algorithm is a Metropolis-within-Gibbs sampler [START_REF] Robert | Monte Carlo statistical methods[END_REF] whose four steps are detailed in Appendix C.

Label switching problem

The label switching problem is generally solved by specific procedures [START_REF] Stephens | Dealing with label switching in mixture models[END_REF]. However, based on the argument of [START_REF] Jacques | Model-based clustering for multivariate partial ranking data[END_REF], these techniques are mainly effective when g is known.

However, when the model is used for clustering, the number of classes is unknown, and the model selection is performed using the BIC criterion [START_REF] Schwarz | Estimating the dimension of a model[END_REF] which simultaneously avoids the label switching phenomenon. Indeed, on the one hand, this criterion selects quite separate classes when the sample size is small. Hence, label switching is not present (with high probability) in practice because of this class separation. On the other hand, even though it can select more classes when the sample size increases, the label switching problem does not occur since this phenomenon vanishes asymptotically.

Obviously, when the number of classes is fixed and the sample size is small, the label switching problem can occur. In such cases, we obviously advise using the procedures of [START_REF] Stephens | Dealing with label switching in mixture models[END_REF].

Simulations

In this section, two simulations are used to illustrate the new model. The first simulation shows the relevance of the estimation procedure. The second simulation illustrates the robustness of the proposed model by analyzing data sampled from a mixture of Poisson distributions.

Experiment conditions For each situation, 100 samples are generated. Parameters are estimated by taking the mean of the parameters sampled by 10 3 iterations of Algorithm 3.1 after a burn-in period of 10 2 iterations. Algorithm 3.1 is initialized with the maximum likelihood estimator of the locally independent model (particularly relevant when within-class dependencies are small). The Kullback-Leibler divergence [START_REF] Kullback | On information and sufficiency[END_REF] is used to compare the estimated distribution and the distribution used to sample the data. This divergence is approximated via 10 4 iterations of a Monte Carlo simulation.

Estimation efficiency

Data sets are composed of one continuous variable, one integer variable and one binary variable. They are sampled from the distribution described in the example in Section 2.4. The results are presented in Figure 2. According to Figure 2(a), the estimated distribution converges to the true distribution when the sample size increases. Indeed, the Kullback-Leibler divergence of the estimated model from the true model decreases as a function of sample size and converges to zero. Moreover, as shown by Figure 2(b), the misclassification rate converges to the theoretical misclassification rate (equal to 0.005) when n increases. This simulation illustrates the convergence of the estimator computed by averaging the parameters sampled using the algorithm. Finally, the estimation procedure is not particularly time consuming since samples of size n = 100 and n = 1600 take 15 and 64 seconds respectively, on an Intel Core i5-3320M processor.

Robustness

Samples are generated from the bivariate Poisson mixture model [START_REF] Karlis | Exact Bayesian modeling for bivariate Poisson data and extensions[END_REF] with π = (1/3, 2/3), whose univariate margin parameters α k = (λ k1 , λ k2 , λ k3 ) take on the following values: λ 1h = h and λ 2h = 3 + h, for h = 1, 2, 3 (see notation detailed in [START_REF] Karlis | Exact Bayesian modeling for bivariate Poisson data and extensions[END_REF]). Figure 3 presents the results and shows the robustness of the Gaussian copula mixture model since it efficiently manages such data sets, as detailed below.

As shown by Figure 3(a), the resulting misclassification rate converges to the theoretical misclassification rate (equal to 0.0967). Moreover, Figure 3 q q q q q q q q q q q 100 200 400 800 1600 0.00 4). In the same way, Figure 3(d) shows that the univariate margin parameter of variable 1, for component 1, converges to its theoretical value (equal to λ 13 / (λ 11 + λ 13 )(λ 12 + λ 13 ) = 3/ √ 20 0.67). Finally, the estimation procedure is not excessively time consuming since it takes 12 and 54 seconds to analyze samples of size n = 100 and n = 1600 respectively, on an Intel Core i5-3320M processor.

Applications

In this section, we analyze two real data sets with the proposed Gaussian copula mixture model. For each number of components, 10 runs of Gibbs sampler are performed with 1000 iterations after a burn-in period of 100 iterations.

Information criteria (BIC criterion [START_REF] Schwarz | Estimating the dimension of a model[END_REF], ICL criterion [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]...) can be used to perform the model selection. These asymptotic criteria are computed with the estimate provided by the Gibbs sampler and not with the maximum likelihood estimate. However, this procedure remains valid since the BIC criterion can be computed with the estimate of the maximum a posteriori [START_REF] Lebarbier | Une introduction au critère BIC : fondements théoriques et interprétation[END_REF]. Moreover, the ICL criterion can be computed by penalizing the BIC criterion with a term of entropy.

These criteria require the computation of the number of parameters. The locally independent model involves ν Loc = (g -1) + g e j=1 ν j parameters, where ν j is the number of parameters involved by the univariate margin distribution of one component (i.e. ν j = 2 if x j is continuous, ν j = 1 if x j is integer and ν j = m j -1 if x j is ordinal with m j modalities). The heteroscedastic Gaussian copula mixture model involves ν He = ν Loc +g e(e-1) 2 parameters and the homoscedastic Gaussian copula mixture model requires ν Ho = ν Loc + e(e-1) 2 parameters. (a) Misclassification error rate q q q q q q q q q q q q q 100 200 400 800 1600 0.00 (c) Correlation coefficient between both variables for component 1 q q q q q q q q q 100 200 400 800 1600

South African Hearth data set

The data

The data are available at http://sci2s.ugr.es/keel/dataset.php?cod=184. This data set is a retrospective sample of males in a heart-disease high-risk region of the Western Cape, South Africa. Many of the coronary heart disease (CHD) positive men have undergone blood pressure reduction treatment and other programs to reduce their risk factors after their CHD event.

In some cases the measurements were made after these treatments. The class label indicates if the person has coronary heart disease (negative or positive) and is hidden for our analysis. Individuals are described by the following nine variables. The continuous variables are systolic blood pressure (sbp), cumulative tobacco (tobacco), low density lipoprotein cholesterol (ldl), adiposity, obesity and current alcohol consumption (alcohol). The integer variables are type-A behavior (typea) and age at onset (age). Finally, a binary variable indicates the presence or not of heart disease in the family history (famhist).

Model selection

Three mixture models (locally independent, heteroscedastic and homoscedastic mixture of Gaussian copulas) are fitted for various numbers of components. Table 1 presents the values of information criteria used to select the homoscedastic tri-component Gaussian copula mixture model. This model obtains the best results for fitting the data distribution (BIC) and for fitting the best partition (ICL). Moreover, this model detects less components than the locally independent model, thus its interpretation is easier. Note that the BIC criterion selects five components for the locally independent model while the ICL criterion selects four components.

Partition study

The competing approaches overestimate the number of components since the true number of classes is two. However, the Gaussian copula mixture models (homo. and hetero.) select less classes than the locally independent mixture, which needs five components. Note that the partition provided by all of these models is strongly different to the true one since these models obtain an ARI of 0.02. We note that the equality constraint of the covariance matrix increases the value of the information criteria obtained by the homoscedastic model, but it also affects its resulting partition according to the partition resulting from the heteroscedastic model (see Figure 4 shows the PCA visualization for component 3. The tri-component homoscedastic Gaussian copula mixture model provides three well-separated classes as shown by the factorial representation of Figure 4(a). We can see that class 2 (red triangles) is an "intermediate" class, while class 1 (red dots) is strongly different to class 3 (black dots). We now detail the model interpretation. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q -4 -2

Interpretation of best-fit model

The PCA visualization (Figure 4) shows that the individuals of class 3 have riskier behaviour than the others (high tobacco and alcohol consumption, older population, high level of obesity).

A three-level interpretation (proportions, univariate marginal distributions and within-class dependencies) is done using the model parameters. The main characteristics of the variables are summarized in Figure 5: • Class 1 (weak-risk behaviours): this is the smallest class (π 1 = 0.07). It contains the young individuals with low alcohol and tobacco consumption. This class groups 57 individuals where only one has a coronary heart disease.

• Class 2 (moderate-risk behaviours): this class is moderated in size (π 2 = 0.24). This class is composed of individuals which high alcohol consumption. The other variables take intermediate values. Among the 106 individuals belonging to this class, 20 have coronary heart disease.

• Class 3 (high-risk behaviours): this is the biggest class (π 3 = 0.69). It contains the individuals with the highest risk behaviour. Of the 299 individuals in this class, 139 have coronary heart disease.

Finally, for all the components, age and the consumption of alcohol and tobacco are strongly linked (see Figure 4(b)). Moreover, adiposity and obesity are also strongly linked.

Conclusion

For such data, the Gaussian copula mixture model reduces the drawbacks of the locally independent model. By decreasing the number of components, it yields a more interpretable model that better fits the data (BIC criterion) and provides a pertinent partition (ICL criterion). Finally, the estimation of main within-class dependencies, based on PCA outputs per component, is an efficient tool for refining the interpretation.

Forest fire data set

The data

The data are composed of 517 forest fires that have occurred in the northeast region of Portugal [START_REF] Cortez | A data mining approach to predict forest fires using meteorological data[END_REF]. These forest fires are described by the following meteorological variables: seven continuous variables (four fire weather index (FWI) system variables, i.e. fine fuel moisture code (FFMC), duff moisture code (DMC), drought code (DC), initial season index (ISI), and three meteorological variables, i.e. temperature (Temp), relative humidity (RH) and wind) and three binary variables indicating the presence of rain, the season (summer or other) and the day of the week (weekend or other).

Model selection

Table 3 presents the values of information criteria used to distinctly select the heteroscedastic tri-component Gaussian copula mixture model. Note that the locally independent model degenerates with five components. The heteroscedastic model with three components obtains better values for the information criteria than the locally independent model since it models the within-component dependencies. Moreover, as shown by Table 4, these within-component dependencies influence the resulting partition since four individuals are affiliated in different classes by both models.

Interpretation of best-fit model

The three-step interpretation of the heteroscedastic tri-component Gaussian copula mixture model is presented using the parameters summarized in Figure 6 Table 4: Confusion matrix for the partitions provided by the locally independent model with three components (rows) and the heteroscedastic Gaussian copula mixture model with three components (columns).

• Class 1 (unpredictable fires): this is the smallest class (π 1 = 0.09). It contains fires that are difficult to predict since they occur with small values of the four FWI system variables (especially ISI). These fires occur during all the year but only when the weather is dry.

• Class 2 (predictable summer fires): this is the biggest class (π 2 = 0.78). This class is composed of fires occurring mainly in summer. They appear with high values of the four FWI system variables and a high temperature.

• Class 3 (winter fires): this is a class of moderate size (π 3 = 0.13). It contains fires occurring in winter, so with low temperature. They occur in dry weather and very small values of the four FWI systems (except ISI).

The correlation matrices highlight dependencies between the summer period and high temperatures, and between the FFMC and DMC values (see Figure 7(b) for component 3). Finally, it should be noted that the variable indicating the day of the week roughly follows the same distribution for all three classes.

The results of the PCA done according to component 3 is shown in Figure 7. Obviously, Figure 7(a) shows that the individuals belonging to component 3 are strongly different to the other ones. Indeed, few individuals belonging to component 2 are visible in this map. The other individuals are too far away from the origin. Thus, the distribution of component 3 is strongly different the other distributions.

Conclusion

The cluster analysis obtained with the Gaussian copula mixture model is more accurate than the one obtained with the locally independent model. Indeed, it provides a meaningful model which needs less components. This model sheds light on three kinds of fires: fires predictable with the FWI system (class 2), fires occurring in winter (class 3) and unpredictable fires (class 1). q q q q q q q q q -20 -10 0 10 20

Conclusion and future extensions

A Gaussian copula mixture model has been introduced and used to cluster mixed data. Using Gaussian copulas, the univariate marginal distributions of each component follow conventional distributions, and within-class dependencies are effectively modelled. Thus, the model results can be easily interpreted. Using the continuous latent variables of Gaussian copulas, a PCAtype method allows for component-based visualization of individuals. Moreover, this approach provides a summary of within-component dependencies, which avoids tedious interpretation of correlation matrices.

In the description of numerical experiments and applications, we pointed out that this model is sufficiently robust to fit data obtained from another model. Furthermore, it can reduce the bias produced by the locally independent model (e.g. reduction of the number of components).

The number of parameters increases with the number of components and number of variables, particularly due to the correlation matrices of the Gaussian copulas. In order to overcome this drawback, we have proposed a homoscedastic version of the model which assumes equality between correlation matrices. However, the number of parameters of this model is still a quadratic function of the number of variables. Therefore, more parsimonious correlation matrices could be proposed in future studies for clustering high-dimensional mixed data.

Since the distribution of all the variables is modeled, this model could be used to manage data sets with missing values. By assuming that values missing at random, the Gibbs sampler could also be adapted, but the underlying principle would remain roughly the same.

Finally, the proposed model cannot cluster non-ordinal categorical variables having more than two modalities. In such cases, the cumulative distribution function is not defined. An artificial order between modalities could be added to define a cumulative distribution function, but this method presents three potential difficulties that require attention: it assumes regular dependencies between the modalities of two variables, its estimation would slow down the estimation algorithm, and its stability would have to be verified.

MixCluster (https://r-forge.r-project.org/R/?group_id=1939) is an R package which performs the cluster analysis method described in the article. It also contains the data sets used in this paper.

A Proof of model identifiability

Model identifiability is proved by two propositions. The first proposition proves model identifiability when the variables are continuous and/or integer. This proposition presents the reasoning in a simple case since it does not consider the ordinal variables. The second proposition proves that the model requires at least one continuous or integer variable to be identifiable.

Proposition A.1 (Identifiability with continuous and integer variables). The Gaussian copula mixture model is weakly identifiable [START_REF] Teicher | Identifiability of Finite Mixtures[END_REF] if the variables are continuous and integer ones ( i.e. the univariate marginal distributions of the components are Gaussian or Poisson distributions). Thus,

∀x ∈ R c × N d , g k=1 π k p(x|α k ) = g k=1 π k p(x|α k ) (10) ⇒ g = g , π = π , α = α . ( 11
)
Proof. The identifiability of the multivariate Gaussian mixture models and of the univariate Poisson mixture model [START_REF] Teicher | Identifiability of Finite Mixtures[END_REF][START_REF] Yakowitz | On the identifiability of finite mixtures[END_REF]) means that (10) implies g = g , π = π , β kj = β kj and Γ kcc = Γ kcc .

We now show that Γ kcd = Γ kcd and Γ kdd = Γ kdd . Let j ∈ {1, . . . , c} and h ∈ {c + 1, . . . , e}. We denote by

ρ k = Γ k (j, h), ρ k = Γ k (j, h), v k = Φ -1 1 (P (x j ; β kj )), ε k (x j ) = π k φ 1 (v k ) σ kj , a k = b ⊕ k (x j )-ρ k v k √ 1-ρ 2 k and a k = b ⊕ k (x j )-ρ k v k √ 1-ρ 2 k . Without
loss of generality, we order the components such that σ kj > σ k+1j and if σ kj = σ k+1j then µ kj > µ k+1j , then (10) implies that

1 + g k=2 (ε k (x j )Φ(a k ))/(ε 1 (x j )Φ(a 1 )) = g k=1 ε k (x j )Φ(a k )/(ε 1 (x j )Φ(a 1 )). Let γ t = {(x j , x h ) ∈ R × N : a 1 = t}. Then, letting x h → ∞ such that (x j , x h ) ∈ γ t , ∀t, a 1 t φ(u)du Φ(t) = 0. ( 13
)
Thus a 1 = a 1 , so ρ 1 = ρ 1 . Repeating this argument for k = 2, . . . , g and for all the couples (j, h), we conclude that Γ kcd = Γ kcd . When both variables are integer, we use the same argument with γ

(t,ξ) = {(x j , x h ) ∈ N × N : a 1 ∈ B(t, ξ)}. Note that if ρ 1 = ρ 1 then ∃n 0 such that ∀x j > n 0 a 1 > t+ξ. Letting x h → ∞ such that (x j , x h ) ∈ γ (t,ξ) , we obtain the following contradiction a 1 t+ξ φ(u)du Φ(t-ξ) = 0 and a 1 t+ξ φ(u)du Φ(t-ξ) > 0. So, a 1 = a 1 then ρ 1 = ρ 1 .
Repeating this argument for k = 2, . . . , g and for all the couples (j, h), we conclude that Γ kdd = Γ kdd .

Proposition A.2 (Identifiability of the Gaussian copula mixture model). The Gaussian copula mixture model is weakly identifiable [START_REF] Teicher | Identifiability of Finite Mixtures[END_REF] if at least one variable is continuous or integer.

Proof. In this proof, we consider only one continuous variable and two binary variables. Obviously, the same reasoning can be extended to the other cases. We now show that Γ kcd = Γ kcd and Γ kdd = Γ kdd .

Let j = 1 and let h ∈ {2, 3}. We note

ρ k = Γ k (j, h), ρ k = Γ k (j, h), v k = Φ -1 1 (P (x j ; β kj )), ε k (x j ) = π k φ(v k ;0,1) σ kj , a k = b ⊕ k (x j )-ρ k v k √ 1-ρ 2 k and a k = b ⊕ k (x j )-ρ k v k √ 1-ρ 2 k
. Without loss of generality, we order the components such that σ kj > σ [k+1]j and if

σ kj = σ [k+1]j then µ kj > µ [k+1]j . Note that (10) implies that 1 + g k=2 (ε k (x j )Φ(a k ))/(ε 1 (x j )Φ(a 1 )) = g k=1 ε k (x j )Φ(a k )/(ε 1 (x j )Φ(a 1 )).
Letting x 1 → ∞ and assuming that ρ k > 0 then

Φ(a k ) Φ(a k ) = 1. So, sign(ρ k ) = sign(ρ k ). By denoting κ = lim a→∞ φ(a) Φ(a) and letting x 1 → ∞, κ 1 κ φ(a k ) φ(a k ) = 1. Thus a 1 = a 1 , so ρ 1 = ρ 1 and b ⊕ k (x j ) = b ⊕ k (x j ) so β kh = β kh .
Note that the same result can be obtained by tending x 1 to -∞ if ρ k < 0. Repeating this argument for k = 2, . . . , g and for all the couples (j, h), we conclude that Γ kcd = Γ kcd then Γ kdd = Γ kdd .

B Prior distributions

We assume independence between the parameters as follows

p(θ) = p(π) g k=1   p(Γ k ) e j=1 p(β kj )   . (14)
The classical conjugate prior distribution of the proportion vector is the Jeffreys non informative one, which is the following Dirichlet distribution

π ∼ D g 1 2 , . . . , 1 2 . ( 15
)
If x j is continuous, then β kj denotes the parameters of a univariate Gaussian distribution so p(β kj ) = p(µ kj |σ 2 kj )p(σ 2 kj ) with

σ 2 kj ∼ G -1 (c 0 , C 0 ) and µ kj |σ 2 kj ∼ N 1 (b 0 , σ 2 kj /N 0 ), (16)
where G -1 (., .) denotes the inverse gamma distribution. With an empirical Bayesian approach, the hyper-parameters (c 0 , C 0 , b 0 , N 0 ) are fixed as proposed by [START_REF] Raftery | Hypothesis testing and model selection[END_REF], so c 0 = 1.28, C 0 = 0.36Var(x j ), b 0 = 1 n n i=1 x j i and N 0 = 2.6 argmax x j -argmin x j . If x j is integer, β kj denotes the parameter of a Poisson distribution and

β kj ∼ G(a 0 , A 0 ). (17)
According to [START_REF] Frühwirth-Schnatter | Finite mixture and Markov switching models[END_REF], the values of hyper-parameters a 0 and A 0 are empirically fixed to a 0 = 1 and A 0 = a 0 n/ n i=1 x j i . If x j is ordinal, β kj denotes the parameter of a multinomial distribution and its Jeffreys non informative conjugate prior means that

β kj ∼ D m j 1 2 , . . . , 1 2 . ( 18
)
The conjugate prior of a covariance matrix is the Inverse Wishart distribution denoted by W -1 (., .). Therefore, it is natural to define the prior of the correlation matrix Γ k from the prior of the correlation matrix Λ k . Indeed, Γ k |Λ k is deterministic [START_REF] Hoff | Extending the rank likelihood for semiparametric copula estimation[END_REF]. So,

Λ k ∼ W -1 (s 0 , S 0 ) and ∀1 ≤ h, ≤ e, Γ k [h, ] = Λ k [h, ] Λ k [h, h]Λ k [ , ] , (19)
where (s 0 , S 0 ) are two hyper-parameters. However, an empirical Bayesian approach cannot be fitted to these parameters since y is not observed. Uniform distribution on ] -1, 1[ is also obtained for the margin distributions of each correlation coefficient by setting s 0 = e + 1 and S 0 equal to the identity matrix [START_REF] Barnard | Modeling covariance matrices in terms of standard deviations and correlations, with application to shrinkage[END_REF].

C Metropolis-within-Gibbs sampler

This section explains the four samplings used in Algorithm 3.1. The first two samplings are difficult to perform directly, so they are done by one iteration of two Metropolis-Hastings algorithms.

For both Metropolis-Hastings algorithms, the instrumental distributions assume conditional independence between parameters. So the smaller the within-class dependencies are, the closer the instrumental distributions of both algorithms are to the stationary distributions. Finally, the last two samplings used in Algorithm 3.1 are classical.

C.1 Class membership and Gaussian vector sampling

The aim is to sample from (6) but the sampling from z, y|x, θ (r-1) cannot be performed directly. So it is achieved by one iteration of the following Metropolis-Hastings algorithms. The sampling from ( 6) is performed in two steps using independence between the individuals, which implies that p(z, y|x, θ (r-1) ) = n i=1 p(z i |x i , θ (r-1) )p(y i |x i , z i , θ (r-1) ). ( 20)

Firstly, each z (r) i is independently sampled from the multinomial distribution z i |x i , θ (r-1) ∼ M g (t i1 (θ (r-1) ), . . . , t ig (θ (r-1) )),

where t ik (θ (r-1) ) = . Note that t ik (θ (r-1) ) is the posterior probability that x i arises from component k with the parameters θ (r-1) .

Secondly, each y (r-1/2) i is independently sampled given (x i , z

i , θ (r-1) ) ∝ φ d (y d i ; µ

d(r-1) z (r) i , Σ d(r-1) z (r) i )1 {y d i ∈S z (r) i (x d i )} , (22)
where µ d(r-1) z .

The computation of t ik (θ (r-1) ) involves the calculation of the integral defined in (5) which can be time consuming if d is large (d > 6). In such cases, the sampling from ( 6) is replaced by one iteration of the Metropolis-Hastings algorithm, which independently samples each couple (z i , y i ). Its stationary distribution is p(z i , y i |x i , θ (r-1) ) ∝ π z i p(x i , y i |z i , θ (r-1) ).

(23) Note that p(x i , y i |z i , θ (r-1) ) = φ e (y i ; 0, Γ (r-1) z i

)1 {y c i =Ψ(x c i ;α (r-1) z i

)} 1 {y d i ∈Sz i (x d i )} . The Metropolis-Hastings algorithm samples a candidate (z i , y i ) by the instrumental distribution q 1 (.|x i , θ (r-1) ), which uniformly samples z i , then samples y i |z i as follows. Its first c elements, denoted by y c i , are equal to y c i = Ψ(x c i ; α (r-1) z i

). Its last d elements, denoted by y d i , follow a multivariate independent Gaussian distribution truncated on S z i (x d i ). Thus, q 1 (z i , y i |x i , θ (r-1) ) = 1 g φ d (y d i |0, I)1 {y d i ∈Sz i (x d i )} e j=c+1 p(x j i ; β (r-1)

z i j ) 1 {y c i =Ψ(x c i |α (r-1) z i )} . ( 24
)
The candidate is accepted with the probability

ρ (r) 1i = min      q 1 (z (r-1) i
, y (r-1) i |x i ) q 1 (z i , y i |x i ) π z i φ e (y i ; 0, Γ (25) Thus, at iteration (r) of Algorithm 3.1, the sampling according to ( 6) is performed via one iteration of the following Metropolis-Hastings algorithm having p(z i , y i |x i , θ (r-1) ) as its stationary distribution.

Figure 1 :

 1 Figure 1: Example of data visualization.

Figure 1

 1 Figure1provides an example of data visualization. Figure1(a)shows the scatterplot of the individuals in their native space. Figure1(b) presents the scatterplot of the individuals in the first PCA-map of the second component (red). It allows two classes to be easily distinguished: a centred one (red) and a second one (black) located on the left side. More precisely, the first axis (explained by the continuous and the integer variables) is strongly discriminative while the second axis (explained exclusively by the binary variable) is not discriminative. Figure1(c)shows the correlation circle of the first PCA-map of the red component. It allows a strong correlation to be identified, for the red component, between the continuous and the integer variables.

  (b) shows that the Kullback-Leibler divergence almost vanishes when the sample size increases, thus demonstrating the flexibility of the Gaussian copula mixture model. Furthermore, the resulting parameters reflect the main properties of the true distribution. Indeed, Figure3(c) shows that the correlation coefficient between both variables, for component 1, converges to its theoretical value (equal to λ 11 + λ 13 = Kullback-Leibler divergence of the estimated model from the true model.

Figure 2 :

 2 Figure 2: Boxplots of indicators of the good behavior of the estimation procedure for different sample sizes.

Figure 3 :

 3 Figure 3: Boxplots of the indicators related to the estimated model. Values obtained with the true Poisson mixture model are indicated by the horizontal black lines.

Figure 4 :

 4 Figure 4: Visualization using parameters of component 3 for the tri-component homoscedastic Gaussian copula mixture model.

Figure 5 :

 5 Figure 5: One dimensional marginal distributions for the South African Hearth data: whole distribution (gray), component 1 (black), component 2 (red) and component 3 (green) of the tri-component homoscedastic Gaussian copula mixture model.

Figure 6 :

 6 Figure 6: One dimensional marginal distributions for the forest fire data: whole distribution (gray), component 1 (black), component 2 (red) and component 3 (green) of the tri-component homoscedastic Gaussian copula mixture model.

Figure 7 :

 7 Figure 7: Visualization based on the parameters of component 3 for the tri-component heteroscedastic Gaussian copula mixture model.

  i |θ (r-1) )

  last d elements, denoted by y d(r-1/2) i , are sampled from the d-variate Gaussian distribution N d (0, Γ

Table 1 :

 1 Values of the BIC and ICL criteria obtained on the South African Hearth data set (best values are in bold).

	g	1	2	3	4	5
	BIC loc. indpt. -14127.26 -13131.88 -12813.92 -12829.68 -12738.66
	homo.	-14724.98 -13016.09 -12739.94 -12774.15 -12927.45
	hetero.	-14724.98 -13076.93 -12971.72 -13071.92 -13253.06
	ICL loc. indpt. -14127.26 -13144.21 -12832.12 -12887.19 -12805.68
	homo.	-14724.98 -13028.07 -12762.79 -12816.44 -12979.06
	hetero.	-14724.98 -13085.52 -12989.06 -13103.61 -13299.16

Table 2

 2

	).
			hetero.
		Class 1 Class 2 Class 3
	Class 1	46	11	0
	homo. Class 2	0	92	14
	Class 3	0	3	296

Table 2 :

 2 Confusion matrices between the tri-component homoscedatic model (row) and the tri-component homoscedastic model (column).

Table 3 :

 3 Values of the BIC and ICL criteria obtained on the forest fire data set (best values are in bold).

	:

Algorithm C.1.

i , y

(r-1/2) i ) = (z i , y i ) with probability ρ (r) 1i

(z (r-1) i

, y (r-1) i

) with probability 1 -ρ (r) 1i .

(27)

C.2 Margin parameter and Gaussian vector sampling

The aim is to sample from (7) but the sampling from β kj |x, y (r)

[rk] , z (r) , β

k , Γ k cannot be performed directly. So, it is done by means by one iteration of the following Metropolis-Hastings algorithms. The sampling from ( 7) is performed using the following decomposition

). ( 28)

kj is sampled first. The full conditional distribution of β kj is defined up to a normalizing constant such that p(β kj |x, y

The distribution of x j i |y

where the real μi = Γ

is the full conditional mean of y j i , Γ k [j, ] being the row j of Γ k deprived of element j and Γ k [, ] being the matrix Γ k deprived of the row and the column j, and where σ2 i is the full conditional variance of y j i defined by σ2

As the normalizing constant of ( 29) is unknown, β (r) kj cannot be directly sampled. This problem is avoided by one iteration of the Metropolis-Hastings algorithm. The instrumental distribution of this Metropolis-Hastings algorithm q 2 (.|x, z) samples a candidate β kj according to the posterior distribution of β kj under a conditional independence assumption (this distribution is explicit since the conjugate prior distributions are used). So, q 2 (.|x, z) = p(β kj |x, z, Γ k = I). Thus, according to (29), the candidate β kj is accepted with the probability

.

Thus, at iteration (r) of Algorithm 3.1, step ( 7) is performed via one iteration of the following Metropolis-Hastings algorithm whose stationary distribution is p( (32)

Vector y j(r)

[rk] is easily sampled after β (r) kj . Indeed, independence between the individuals defines the full conditional distribution of y j

[rk] by

). ( 33) If x j is a continuous variable (i.e. 1 ≤ j ≤ c), when z (r) i

= k, the full conditional distribution of y j i is a Dirac distribution at

. If x j is a discrete variable (i.e. c + 1 ≤ j ≤ e), when z (r) i = k, the full conditional distribution of y j i is a truncated Gaussian distribution such as,

So, step ( 7) is performed in two stages. First, β

kj is sampled via one iteration of the Metropolis-Hastings algorithm of which the stationary distribution is p(β kj |x, y

[rk] is sampled from (34).

C.3 Proportion vector sampling

The aim is to sample from (8). The sampling from ( 8) is classical. Indeed, the conjugate Jeffreys non informative prior means that

where n (r)

C.4 Correlation matrix sampling

The aim is to sample from (9). To sample from (9), we use the approach proposed by [START_REF] Hoff | Extending the rank likelihood for semiparametric copula estimation[END_REF] in the case of semiparameteric Gaussian copulas. First, a covariance matrix is generated by its explicit posterior distribution, and second, the correlation matrix is deduced by normalizing the covariance matrix. As (y, z) are known in this step, we are in the well-known case of a multivariate Gaussian mixture model with known means. Thus, the sampling according to Γ k |y (r) , z (r) is performed by the following two steps

where

. As the homoscedastic model assumes equality between the correlation matrices, in this case we only sample one Λ so (36) is replaced by

and we put Λ k = Λ for k = 1, . . . , g.