Hicham Mhamdi

Emmanuelle Stoetzel

Oussama Bouarakia

Hassan Aouraghe

Mohamed Souhir

Morphological, biometric and taphonomic analyses of Holocene gerbils from Guenfouda (Eastern Morocco) and first mention of Gerbillus henleyi in a fossil record

Keywords: Rodents, Gerbillinae, Holocene, eastern Morocco, palaeontology, taphonomy, palaeoenvironments

The Guenfouda cave, located in eastern Morocco, has yielded an abundant macro-and microvertebrate fauna associated with a rich lithic industry attributed to several cultures (Middle Stone Age, Later Stone Age and Neolithic). Among the microfauna, on the basis of new morphological and biometric data of the upper and lower first molars, we identified two gerbillines of the genus Gerbillus collected from the upper levels of an Holocene age. Gerbillus henleyi is here identified for the first time in the fossil record of North Africa, while G. campestris has been cited previously at this site and in many other Pleistocene and Holocene sites in the Maghreb. The taphonomic study concluded that a medium-sized nocturnal predator (possibly Tyto alba) was likely to be the source of the accumulation of small mammals in the filling of the Guenfouda cave. The palaeoecological reconstruction based on the whole faunal spectrum indicates a semiopen steppic landscape, with patches of more wooded areas and water ponds or streams, which is also characteristic for this region at present.

Introduction

As in other countries of North Africa, the microvertebrates of Morocco are still poorly documented and often underrepresented in the archaeological context [START_REF] Stoetzel | Adaptations and Dispersals of Anatomically Modern Humans in the Changing Environments of North Africa: the Contribution of Microvertebrates[END_REF][START_REF] Stoetzel | Quaternary terrestrial microvertebrates from Mediterranean northwestern Africa: state-of-the-art focused on recent multidisciplinary studies[END_REF]. The information obtained from the discovery of abundant microvertebrate remains in the Holocene levels of the Guenfouda cave in eastern Morocco (Aouraghe et al. 2010;[START_REF] López-García | The small mammals from the Holocene site of Guenfouda (Jerada, Eastern Morocco): chronological and paleoecological implications[END_REF] shed new lights on the paleoenvironments and paleoclimates in which the prehistoric human populations of the region evolved. This discovery also filled some gaps in the fossil record of the concerned taxa in eastern Morocco and in the Maghreb in general.

By their abundance, and adaptations to well-defined habitats and precise climatic conditions, micromammals are good ecological markers and allow a fine reconstruction of palaeoenvironments (e.g., [START_REF] Chaline | La signification des Rongeurs dans les dépôts Quaternaires[END_REF]. However, the accumulation of the remains of small vertebrates at an archaeological or palaeontological site is generally (but not exclusively) the result of the intervention of one or more predators, and also of various post-deposition processes (e.g., [START_REF] Denys | Nouveaux critères de reconnaissance des concentrations de microvertébrés d'après l'étude des pelotes de chouettes du Botswana (Afrique australe). Paris (4èm)[END_REF][START_REF] Andrews | Owls, caves & fossils[END_REF]Fernandez-Jalvo and Andrews 1992;[START_REF] Desclaux | Origine(s) des accumulations de micromammifères dans quelques sites préhistoriques du Pléistocène moyen et supérieur d'Europe méridionale[END_REF]. A taphonomic analysis is therefore necessary in order to determine the main accumulating agent as well as the factors that may have affected the assemblage secondarily. This, in turn, ensures the reliability of the palaeoenvironmental reconstructions deduced from the faunal lists and the application of palaeoecological indices.

Previous studies in Guenfouda focused on the fauna in general (Aouraghe et al. 2010), and small mammals in particular [START_REF] López-García | The small mammals from the Holocene site of Guenfouda (Jerada, Eastern Morocco): chronological and paleoecological implications[END_REF]). In the present work, we focused on the gerbillines of the genus Gerbillus, because during the review of the material it became clear that at least two species of very distinct size were present, something unheard of in a Moroccan Holocene site. The genus Gerbillus appeared in the Maghreb during the Pliocene [START_REF] Tong | Origine et évolution des Gerbillidae (Mammalia, Rodentia) en Afrique du Nord[END_REF][START_REF] Geraads | Rongeurs et Insectivores (Mammalia) du Pliocène final de Ahl Al Oughlam (Casablanca, Maroc)[END_REF], and underwent significant diversification during the Pleistocene, a period characterised by strong climatic and environmental variations [START_REF] Stoetzel | Late cenozoic micromammal biochronology of northwestern Africa[END_REF][START_REF] Ndiaye | Evolutionary systematics and biogeography of the arid habitat-adapted rodent genus Gerbillus (Rodentia, Muridae): a mostly Plio-Pleistocene African history[END_REF]. It is probably the establishment of increasingly long and marked arid periods in the Lower and Middle Pleistocene that would have particularly favoured the diversification of Gerbillinae in North Africa. But at the boundary between the Middle and Upper Pleistocene, we observe a significant decrease in diversity within micromammals, with the disappearance of many fossil species (and even entire lineages) including within the subfamily Gerbillinae [START_REF] Geraads | The rodents from the late middle Pleistocene hominid-bearing site of J'bel Irhoud, Morocco, and their chronological and paleoenvironmental implications[END_REF][START_REF] Stoetzel | Late cenozoic micromammal biochronology of northwestern Africa[END_REF]. At the Upper Pleistocene, within the genus Gerbillus, only G. campestris, and very occasionally G. maghrebi, are recorded in the archaeological and palaeontological sites of the Maghreb.

The systematics, phylogeny and phylogeography of the current Gerbillinae have undergone a significant development in recent years (e.g., [START_REF] Chevret | Systematics and evolution of the subfamily Gerbillinae (Mammalia, Rodentia, Muridae)[END_REF][START_REF] Pavlinov | A Review of Phylogeny and Classification of Gerbillinae (Mammalia: rodentia)[END_REF][START_REF] Abiadh | Molecular phylogenetics of the genus Gerbillus (Rodentia, Gerbillinae): implications for systematics, taxonomy and chromosomal evolution[END_REF][START_REF] Ito | Molecular phylogeny of the subfamily Gerbillinae (Muridae, Rodentia) with emphasis on species living in the Xinjiang-Uygur autonomous region of China and based on the mitochondrial cytochrome b and cytochrome c oxidase subunit II genes[END_REF][START_REF] Ndiaye | Evolutionary systematics and biogeography of endemic gerbils (Rodentia, Muridae) from Morocco: an integrative approach[END_REF][START_REF] Ndiaye | Evolutionary systematics and biogeography of the arid habitat-adapted rodent genus Gerbillus (Rodentia, Muridae): a mostly Plio-Pleistocene African history[END_REF][START_REF] Aghova | Multilocus phylogeny of East African gerbils (Rodentia, Gerbilliscus) illuminates the history of the Somali-Masai savanna[END_REF][START_REF] Alhajeri | Craniomandibular Variation in the Taxonomically Problematic Gerbil Genus Gerbillus (Gerbillinae, Rodentia): assessing the Influence of Climate, Geography, Phylogeny, and Size[END_REF][START_REF] Nanova | Phylogeography, phylogeny, and taxonomical revision of the Midday jird (Meriones meridianus) species complex from Dzungaria[END_REF][START_REF] Bryja | Rodents of the Afar Triangle (Ethiopia): geographical isolation causes high level of endemism[END_REF]. However, there is still no consensus on the position of some taxa such as Dipodillus, which is, according to the authors, considered as a synonym of Gerbillus (e.g., [START_REF] Chevret | Systematics and evolution of the subfamily Gerbillinae (Mammalia, Rodentia, Muridae)[END_REF], a separate genus (e.g., [START_REF] Musser | Superfamily Muroidea[END_REF][START_REF] Pavlinov | A Review of Phylogeny and Classification of Gerbillinae (Mammalia: rodentia)[END_REF]) or a sub-genus of Gerbillus (e.g., [START_REF] Abiadh | Molecular phylogenetics of the genus Gerbillus (Rodentia, Gerbillinae): implications for systematics, taxonomy and chromosomal evolution[END_REF]Ndiaye et al;2012;[START_REF] Nicolas | Phylogeography of the North African dipodil (Rodentia: muridae) based on cytochrome-b sequences[END_REF][START_REF] Alhajeri | Molecular systematics of gerbils and deomyines (Rodentia: Gerbillinae, Deomyinae) and a test of desert adaptation in the tympanic bulla[END_REF]. This debate is not part of the objectives of the present study, and consequently we chose to use the name Gerbillus in a broad sense.

The aim of this study is twofold: first, to identify the gerbilline species present in the Guenfouda cave using morphological and biometric criteria, and second, to study the taphonomic context of the assemblage in order to ensure the age of these remains and the reliability of the palaeoenvironmental reconstructions.

The site

The Guenfouda cave (also known as Ghar Zebouj) is located 6 km from the village of the same name and 30 km SW of the city of Oujda in eastern Morocco (Figure 1a).

The cave cut through a cliff of dolomite and oolitic limestone (Figure 1b) of Aaleno-Bajocian age (about 175 million years old), a formation deposited in a marine platform environment. Orographically, it is located on the eastern part of the Oujda Mountains at 930 m above sea level, and it occupies a panoramic position overlooking the plain of Beni Yala (Figure 1c). The arched porch, which opens to the east, rises 4.35 m from the ground at the entrance. The cave is mostly unichamber from the entrance until a central stalagmitic pillar located at the back. From that pillar, the cave starts to divide into two cavities (Figure 1d): a chamber on the right which is about 6 m long and 3.5 m wide, and a smaller room on the left of 3.5 m length by 2.5 m width.

Its large area, good aeration, easy access and its location in a place which overlooks the surrounding areas (Figure 1c) make the cave a suitable place for a human habitat. Indeed, regularly excavated since 2004 by Professor H. Aouraghe and his team from Mohamed First University, the Guenfouda cave yielded abundant and varied archaeological material (lithic and bone industries, ceramics, faunal and human remains) attributed to several cultures (mainly from the Neolithic, but also to the Later Stone Age, and to the Middle Stone Age) [START_REF] Aouraghe | La grotte préhistorique de Guenfouda, Maroc Oriental : les premiers résultats (fouilles 2004-2007)[END_REF].

The filling of the excavated area extends over two metres in depth and comprises the succession of several stratigraphic levels of different sedimentological natures (Figure 2). These levels have recorded several human occupations ranging from the Middle Stone Age (MSA) to the Neolithic. The cave's stratigraphy studied so far is divided into five layers (Figure 2), which are, from top to bottom:

 Layer C0: this surface layer consists of powdery silts and limestone blocks. It yielded bone remains, lithic and bone tools and ceramic shards attributed to the Neolithic.  Layer C1a: characterised by an alternation of beds rich in ashes (whitish levels), charcoals (dark levels) or clays (brown levels). The little disturbed appearance of these layers suggests the good conservation of the underlying strata. A polished axe and a handstone characteristic of the Neolithic have been discovered in this layer. This layer is also particularly rich in terrestrial gastropod shells.  Layer C1b: made up of light brown powdery silts, without any visible structure, encasing limestone blocks mainly fallen from the walls of the cave and various faunal and lithic remains. This layer yielded ceramic furniture of the Cardial type, corresponding to an ancient phase of the Neolithic  Layer C2: made up of dark brown clay, with an aggregate structure, thicker than layer C1b. The transition with the latter is irregular, but clear. Charcoals are relatively numerous and allowed a 14C dating which resulted in an age of 11,150 ± 50 BP (11,904 cal BC) at the top of Layer C2. This layer has yielded archaeological material attributed to Neolithic cultures (lithic and bone industries). Other levels have been reached during the recent excavations, but were not yet studied in detail and are not represented on Figure 2:  Layer C3: divided into three sublevels (C3a, C3b, C3c). The basal part (C3c) of this layer is dominated by powdery ashy clays. This is then topped by an irregular stony level (C3b) formed by centimetric, or even decimetric, angular limestone elements. An irregular and finely laminated limestone crust forms the third sublevel (C3a) which caps the layer. Sub-layer C3c has provided lithic tools of the ‗slats with knocked back' type (Later Stone Age).  Layer C4: provisionally divided into two parts. The basal part (C4b) is relatively friable, quite ashy, and rich in charcoal debris. The top part (C4a) is characterised by its orange colour, a relative induration and a clear lamination. Also, the charcoal remains are very abundant, and scattered centimetric limestone fragments are present. In 2021, this layer provided lithic tools of the Levallois type (Middle Stone Age).  Both layers (C3 and C4) yielded abundant faunal remains (which are still under study). Archaeological excavations have reached a depth of up to two metres in some places, but they have not yet reached the base of the filling. To specify its extent, an electric tomography campaign was carried out by Professor D. Khattach [START_REF] Aouraghe | Actes de la 6ème rencontre des quaternaristes marocains (RQM6). Rapport de fouille: Mission archéologique Guenfouda[END_REF] and assessed the depth of the filling at 5 m at the edges next to the cave walls.

A Neolithic burial was discovered in Layer C2 as well as isolated human remains including teeth, post cranial bones, and some bones of hands and feet [START_REF] Aouraghe | Mission Guenfouda (Maroc oriental) Rapport sur la campagne de Fouille[END_REF]Aouraghe et al., 2005).

The Guenfouda cave also yielded abundant faunal remains (Table 1) represented by at least 15 taxa of larger fauna [START_REF] Aouraghe | La grotte préhistorique de Guenfouda, Maroc Oriental : les premiers résultats (fouilles 2004-2007)[END_REF](Aouraghe et al. , 2010;;[START_REF] Aouraghe | Mission Guenfouda (Maroc oriental) Rapport sur la campagne de Fouille[END_REF][START_REF] Bougariane | Les Vertébrés du Pléistocène terminal-Holocène de quelques sites marocains : Paléontologie, Taphonomie et Archéozoologie[END_REF], and 34 taxa of small vertebrates. Micromammals are composed of two soricomorphs (Crocidura russula, C. whitakeri), one macroscelid (Petrosaltator rozeti), two bats (Myotis blythii, Rhinolophus hipposideros) and seven rodents (Mus spretus, Meriones shawii, Apodemus sylvaticus, Gerbillus campestris, Eliomys sp., Atlantoxerus getulus, Ctenodactylus sp.) [START_REF] López-García | The small mammals from the Holocene site of Guenfouda (Jerada, Eastern Morocco): chronological and paleoecological implications[END_REF].

However, recent reviews of the material have shown that some of the previously identified taxa need reappraisal. For instance, [START_REF] Stoetzel | Adaptations and Dispersals of Anatomically Modern Humans in the Changing Environments of North Africa: the Contribution of Microvertebrates[END_REF] showed that it is rather Meriones grandis, and not M. shawii, that is represented at Guenfouda. We also revealed that the remains previously attributed to Ctenodactylus sp. belong actually to Jaculus sp. (Mhamdi and Stoetzel, pers. obs.), and that a second species of gerbil is present (present work).

Materials and methods

Sampling

The material studied in this work originates from the Layers C0 to C2, excavated in 2016, 2017 and 2019. These excavations were carried out by independent subdivisions of 5 cm thick layers, while respecting the limits of sedimentary layers and archaeological soils. Each bucket of sediment removed during the excavations bore the indications of its origin (number of the square, depth [Z] level measured by laser theodolite). The sediments were then water-sieved on two superimposed sieves of 3 mm and 1 mm mesh in order to recover the smaller archaeological materials (bones of small animals, ostrich eggshell beads, micro-flint shards, etc). The sieved material came mainly from squares I39, I40, J39, J40, K39 and L39 (Figure 1d) which had been the subject of extensive excavations. After drying in the sun, a first selective sorting using tweezers allowed us to recover the remains of microvertebrates. These were then packed in micro-mount boxes bearing the indications of the sample's origin (number of the square, layer, Z level and date).

Because of their proximity to the stalagmitic column, the majority of the fossils from square J40 were covered by concretions, which made their taxonomic and taphonomic study very difficult. To clean them, the samples were submerged in a 10% acetic acid solution for 10 minutes before being rinsed with clean water. These were considered separately for the taphonomic analysis.

Species identification

Gerbillus species were identified mainly according to the anatomical identification criteria established by [START_REF] Tong | Origine et évolution des Gerbillidae (Mammalia, Rodentia) en Afrique du Nord[END_REF], [START_REF] Aulagnier | Dipodillus maghrebi schlitter et setzer, 1972 et Gerbillus campestris levaillant, 1857 (Rodentia, Gerbillidae) dans le nord du Maroc : morphologie et biométrie crâniennes, éléments de répartition[END_REF], [START_REF] Stoetzel | Les Microvertébrés du site d'occupation humaine d'El Harhoura 2[END_REF][START_REF] López-García | The small mammals from the Holocene site of Guenfouda (Jerada, Eastern Morocco): chronological and paleoecological implications[END_REF] and [START_REF] Aulagnier | Mammifères sauvages du Maroc: peuplement, répartition, écologie[END_REF]. The dental nomenclature followed here is that used by [START_REF] Tong | Origine et évolution des Gerbillidae (Mammalia, Rodentia) en Afrique du Nord[END_REF] (Figure 3).

Direct observation of the morphological characters was made thanks to a stereomicroscope (OPTIKA SZM GEM1, magnification x7 to x40). For comparisons, we consulted current genotyped specimens from the collection of O. Bouarakia deposited at the ‗Biodiversity, Ecology and Genome Laboratory', Faculty of Sciences, Mohammed V University, Rabat, as well as photos of specimens from the MNHN in Paris and BMHN in London taken by E. Stoetzel. Materials recovered from rejected pellets and dead individuals collected by H. Mhamdi in the vicinity of the site, made available at the Oujda University Museum, were also used.

The dimensions of the upper (M1) and lower (m1) first molars can be very useful at differentiating certain groups of rodents [START_REF] Geraads | Rongeurs et Lagomorphes du Pleistocène moyen de la Grotte des Rhinocéros, Carrière Oulad Hamida 1 à Casablanca, Maroc[END_REF]. The list of current and fossil material used in this work is presented in Table 2. The Guenfouda fossil assemblage studied includes 107 molars (40 M1: 13 right M1 and 27 left M1; 67 m1: 32 right m1 and 35 left m1), corresponding to a minimum of 35 individuals. The molars were photographed in occlusal view with a camera (OPTIKAM 2350) coupled to a stereomicroscope (OPTIKA SZM GEM1, magnification x40) at the Faculty of Sciences of Oujda. Measurements of different molar distances, indicated in mm, were then taken using the Digimizer software (trial period). The measurements used in this work follow those of [START_REF] Tong | The Gerbillinae (Rodentia) from Tighennif (Pleistocene of Algeria) and their significance[END_REF] and [START_REF] Stoetzel | Les Microvertébrés du site d'occupation humaine d'El Harhoura 2[END_REF]: the length and width of the molars correspond to the greatest anteroposterior distance and the width of the median lobe, respectively.

After entering all data in an Excel table, the usual statistical calculations (sum, mean, standard deviation, and variance) as well as plotting the graphs were carried out using the PAleontological STatistics (PAST) software (version 3.22).

High-resolution photos used for illustrations and some taphonomic observations were taken with the Scanning Electron Microscope (SEM) of the Faculty of Sciences of Oujda.

Taphonomic observations

The taphonomic study of a fossil assemblage requires observing several criteria such as faunal and skeletal representation, fragmentation, and surface alterations such as traces of digestion (e.g., [START_REF] Denys | Nouveaux critères de reconnaissance des concentrations de microvertébrés d'après l'étude des pelotes de chouettes du Botswana (Afrique australe). Paris (4èm)[END_REF][START_REF] Denys | Taphonomie des sites à petits vertébrés, importance de la prédation et des altérations post-prédation[END_REF][START_REF] Denys | Taphonomie des Petits Vertébrés: Référentiels et Transferts aux Fossiles: Actes de la Table Ronde du RTP Taphonomie[END_REF][START_REF] Andrews | Owls, caves & fossils[END_REF]Fernandez-Jalvo and Andrews 1992). A criterion used alone is rarely sufficient to interpret a fossil assemblage; this is why it is often necessary to combine several approaches to obtain a global vision of the taphonomic processes that have affected an assemblage. In the case of an accumulation by predation, the faunal spectrum, the percentage of digested elements and the intensity of digestion, however, remain the most reliable criteria to identify the predator. The anatomical representation and the fragmentation are often subject to post-deposition factors which may mask the signal of the predator's digestion [START_REF] Stoetzel | Implication de l'étude taphonomique des microvetébrés de Témara[END_REF]. Thus, in the context of this preliminary taphonomic study, we only focused on the observation of bone surfaces.

Results

Species identification

In this work, we identified, based on morphological and biometric criteria, two gerbil species from the Guenfouda material: Gerbillus campestris and Gerbillus henleyi. Below we present a description of the material and discuss the geographical range and ecology of the current and fossil populations.

Gerbillus campestris (Loche, 1867)

List of the studied material 6 maxilla (with M1, M2 and M3 in situ); 14 maxilla (with M1 and M2 in situ); 10 maxilla (with M1 in situ); 5 isolated M1; 2 mandibles (with m1, m2 and m3 in situ); 28 mandibles (with m1 and m2 in situ); 15 mandibles (with m1 in situ) and 9 isolated m1.

Previous mentions in the Maghreb fossil record

Forms related to G. campestris are present in the Maghreb fossil record from the Lower Pleistocene, in particular at Jebel Ressas in Tunisia or at Oued Melah in Algeria [START_REF] Ameur-Chehbeur | Biochronologie des formations continentales du Néogène et du Quaternaire de l'Algérie[END_REF][START_REF] Tong | Origine et évolution des Gerbillidae (Mammalia, Rodentia) en Afrique du Nord[END_REF][START_REF] Mein | Gisements karstiques pléistocènes au Djebel Ressas, Tunisie[END_REF]. Later, the species became more frequently recorded during the Middle Pleistocene in Morocco (Hominid cave and Rhinoceros Cave, [START_REF] Geraads | Mammalian faunas from the Pliocene and Pleistocene of Casablanca (Morocco)[END_REF]Salé, Tong 1989;Jebel Irhoud;[START_REF] Geraads | The rodents from the late middle Pleistocene hominid-bearing site of J'bel Irhoud, Morocco, and their chronological and paleoenvironmental implications[END_REF]Sidi Abderrahmane, Tong 1989), as well as in Algeria (Aïn Mefta; Tong 1989) or in Tunisia (Jebel Ressas; [START_REF] Mein | Gisements karstiques pléistocènes au Djebel Ressas, Tunisie[END_REF]. In the Upper Pleistocene and the Holocene, with the exception of a mention of G. maghrebi at Chrafate and Ez-Zarka in northern Morocco [START_REF] Ouahbi | Recent Quaternary fossil mammals of Chrafate and Ez Zarka. The origin of modern fauna in the Northern Rif (NW Morocco, Northern Africa)[END_REF], G. campestris became the only gerbil species cited throughout different sites in Morocco (El Harhoura 2, [START_REF] Stoetzel | Pleistocene and Holocene Small Vertebrates of El Harhoura 2 cave (Rabat-Témara, Morocco) : an annotated preliminary taxonomic list[END_REF]Contrebandiers Cave, Reed and Barr 2010;Felids Cave, Geraads et al. 2010;Chrafate and Ez Zarka;[START_REF] Ouahbi | Recent Quaternary fossil mammals of Chrafate and Ez Zarka. The origin of modern fauna in the Northern Rif (NW Morocco, Northern Africa)[END_REF]Oued Assaka, Wengler et al. 2002;Ifri-el-Baroud;[START_REF] Potì | Human occupation and environmental change in the western Maghreb during the Last Glacial Maximum (LGM) and the Late Glacial. New evidence from the Iberomaurusian site Ifri El Baroud (northeast Morocco)[END_REF]Guenfouda, M et al. 2013;Ifri N ‗Ammar, Mouhsine 2003; Kehf-el-Baroud; [START_REF] Barton | The late upper palaeolithic occupation of the Moroccan Northwest Maghreb during the last glacial maximum[END_REF]Kahf-That-El-Ghar, Ouchaou 2000;Gazelles Cave, Geraads et al. 2010;Hassi Ouenzga, Linstädter 2004) and in Algeria (Gueldaman, [START_REF] Saidani | Nouvelles données sur le contexte taphonomique et environnemental des occupations néolithiques de la grotte de Gueldaman GLD1 (Algérie) d'après l'étude des microvertébrés[END_REF].

Description and comparison

The upper first molars (M1) (Figure 3a) present a protocone that is offset mesially from the paracone. These two tubercles are not individualised and are linked on little-worn teeth by a short ridge. The longitudinal ridges (anterolophule and neolophe), are relatively well developed (high) and alternate [START_REF] Aulagnier | Dipodillus maghrebi schlitter et setzer, 1972 et Gerbillus campestris levaillant, 1857 (Rodentia, Gerbillidae) dans le nord du Maroc : morphologie et biométrie crâniennes, éléments de répartition[END_REF][START_REF] Stoetzel | Les Microvertébrés du site d'occupation humaine d'El Harhoura 2[END_REF]). On the lower first molars (m1) (Figure 3b), the anteroconid is asymmetrical. The anterolophulid connects the latter to the protoconid with a labial position (taking sometimes a more median position). The metaconid and the protoconid are not individualised and are in the opposite position. The longitudinal ridges are well developed, reaching almost the same height as the tubercles (as is the case on M1). Finally, anterolophulid and neolophulid are aligned [START_REF] Tong | Origine et évolution des Gerbillidae (Mammalia, Rodentia) en Afrique du Nord[END_REF][START_REF] Stoetzel | Les Microvertébrés du site d'occupation humaine d'El Harhoura 2[END_REF].

In addition, a m1 (G115-L39) presents an anterior cingulum, as well as a posterior cingulum (Figure 3c). This characteristic can be found frequently in certain fossil gerbil species (G. cingulatus, G. grandis), or in old populations of G. campestris1989; [START_REF] Stoetzel | Les Microvertébrés du site d'occupation humaine d'El Harhoura 2[END_REF]). Thus, our material is similar to the extant populations of G. campestris, and the presence of an anterior and posterior cingulum on some m1 would reflect a primitive character displayed by some individuals.

Measurement

The upper, and especially the lower, molars of G. campestris from Guenfouda appear on average somewhat smaller than other fossil (Table 3, Figure 5) or current (Table 2, Figure 6) populations of this species. But this difference, which does not seem to be significant, could just as well be due to a difference in the measurement method between observers. The values remain within the overall size variability of the species (Table 2,Figures 5,6).

Current distribution and ecology

The North African gerbil G. campestris is widespread in the Mediterranean region of North Africa, from Morocco to Egypt, west of the Nile. It is also present, but more scattered, in the Sahara Desert and the Sahel region (in Mali, Niger and Sudan) [START_REF] Happold | Mammals of Africa Volume III. Rodents, Hares And Rabbits[END_REF][START_REF] Granjon | Gerbillus campestris. The IUCN Red List of Threatened Species 2016[END_REF]. In Morocco, this species is widely distributed from the north of the country to the Drâa Valley in the south, with some isolated populations in the Sahara. This gerbil lives in a variety of habitats, ranging from sub-humid to arid regions, with the exception of high mountains, forests and sandy desert areas. It occupies steppes, arable land, rocky habitats and oases [START_REF] Nicolas | Phylogeography of the North African dipodil (Rodentia: muridae) based on cytochrome-b sequences[END_REF][START_REF] Aulagnier | Mammifères sauvages du Maroc: peuplement, répartition, écologie[END_REF]Bouarakia et al. 2019a).

Gerbillus henleyi (Winton, 1903)

List of the studied material 1 maxilla (with M1, M2 and M3 in situ); 3 maxilla (with M1 and M2 in situ); 1 maxilla (with M1 in situ); 6 mandibles (with m1 and m2 in situ); 6 mandibles (with m1 in situ) and 1 isolated m1.

Previous mentions in the Maghreb fossil record

To our knowledge, this is the first mention of this species in a fossil or sub-fossil state in the Maghreb, or even in the whole of North Africa. Note that a mention of another small gerbil, G. nanus, was made at Ti-n-Torha in Libya, also dated to the beginning of the Holocene [START_REF] Gautier | The fauna of Kadero and the arrival of pastoralism in the Nile Valley of Central Sudan[END_REF]Van Neer, 2011, in Masseti 2010). However, African populations of G. nanus were reattributed to G. amoenus according to recent molecular studies [START_REF] Ndiaye | Molecular variation and chromosomal stability within Gerbillus nanus (Rodentia, Gerbillinae): taxonomic and biogeographic implications[END_REF]Bouarakia et al. 2019b).

Description and comparison

On the M1 (Figure 4a) the anterocone is squeezed on its lingual and labial flanks. The protocone and the paracone, which are slightly offset mesially, are individualised on little worn teeth. The ridges joining the tubercles are lower than these. With wear, save rare exceptions, tubercles progressively fuse in the order anteroconus-protocone-paracone-metacone. The hypocone and metacone are not very distinct on young teeth; with wear, they form a single transverse lophe.

On m1 (Figure 4b), the anteroconid joins the metaconid by the anterolophulid, and later with wear it joins the protoconid as well. The metaconid and the protoconid are individualised and are in opposite position along a lingo-labial axis. The posterior lophe (hypoconid-entoconid), is individualised and directed inwards and a little backwards. The mesosinusid (between metaconid and entoconid) is directed forward in its labial part, almost isolating the metaconid when the teeth start to wear. The longitudinal ridges are less developed, and are slightly offset in young teeth [START_REF] Baker | A morphometric and taxonomic revision of the genus Gerbillus (Mammalia, Rodentia, Gerbillidae) in Jordan with notes on its current distribution[END_REF][START_REF] Siahsarvie | New Records of Naked-footed Gerbil Gerbillus nanus and Pygmy Gerbil Gerbillus cfr. henleyi (Rodentia, Muridae) from Iran[END_REF][START_REF] Stoetzel | Taphonomy and palaeoecology of the late Pleistocene to middle Holocene small mammal succession of El Harhoura 2 cave (Rabat-Témara, Morocco)[END_REF][START_REF] Aouraghe | Mission Guenfouda (Maroc oriental) Rapport sur la campagne de Fouille[END_REF][START_REF] Geraads | The rodents from the late middle Pleistocene hominid-bearing site of J'bel Irhoud, Morocco, and their chronological and paleoenvironmental implications[END_REF][START_REF] Amr | Systematics, distribution and ecological analysis of rodents in Jordan[END_REF]).

Measurements

G. henleyi from Guenfouda exhibit slightly smaller molars than the comparative material, but they remain within the overall size variability of current populations (Table 3,Figures 5,6). However, these can be considered as preliminary observations, pending on a larger sample from current populations.

Current distribution and ecology

The pygmy gerbil G. henleyi lives in the Sahel region and in the southern and northern limits of the Sahara Desert, with an extension into the Middle East and Arabia [START_REF] Musser | Superfamily Muroidea[END_REF][START_REF] Granjon | Gerbillus henleyi Henley's Gerbil[END_REF][START_REF] Monadjem | Rodents of sub-Saharan Africa. A Biogeographic and Taxonomic Synthesis[END_REF]. In Morocco, the distribution area appears sparse with a few capture points or traces of presence reported in the south and southeast of the country [START_REF] Aulagnier | Mammifères sauvages du Maroc: peuplement, répartition, écologie[END_REF]. Recent captures and molecular identifications have confirmed the extension of its distribution range in northeastern Morocco in the Hauts Plateaux [START_REF] Bouarakia | Notes on the distribution and phylogeography of two rare small Gerbillinae (Rodentia, Muridae) in Morocco: Gerbillus simoni and Gerbillus henleyi[END_REF]. This gerbil occupies desert to semi-desert regions and prefers coastal sand dunes covered with halophytes, salt marshes, cultivated fields, wadi beds, stony plains and hamadas [START_REF] Granjon | Gerbillus henleyi Henley's Gerbil[END_REF][START_REF] Aulagnier | Mammifères sauvages du Maroc: peuplement, répartition, écologie[END_REF].

Discussion

The morphology as well as the size of the G. campestris Guenfouda material is similar to specimens described from the neighbouring and contemporary Ifri N'Ammar deposit (Afsou, Nador) by [START_REF] Mouhsine | Contribution à la connaissance de la faune du Pléistocène supérieur du Rif Oriental : la grotte Ifri N'Ammar[END_REF]. However, some differences were observed compared to the current populations of the species. In particular, we noted a slightly smaller size of the molars (Table 3, Figures 5,6) and the persistence of primitive characters (anterior and posterior cingulum) in certain individuals. But these differences remain within the overall variability of the species.

The morphological and biometric study enabled us to identify with certainty G. henleyi at Guenfouda. These specimens did not differ from their present-day counterparts except by their slightly smaller size (Table 3, Figures 5,6). G. henleyi is now part of the current fauna of eastern Morocco, with an expansion of 370 km to the north-east of the country according to the most recent record of this species from south of Tafilalt from rejection pellets [START_REF] Aulagnier | Mammifères sauvages du Maroc: peuplement, répartition, écologie[END_REF][START_REF] Bouarakia | Notes on the distribution and phylogeography of two rare small Gerbillinae (Rodentia, Muridae) in Morocco: Gerbillus simoni and Gerbillus henleyi[END_REF]. This expansion could be interpreted as the result of the ongoing desertification process in North Africa and the increasing aridification advancing northwards to regions such as eastern Morocco [START_REF] Jalut | Holocene circum-Mediterranean vegetation changes: climate forcing and human impact[END_REF][START_REF] Tabel | Environmental changes over the past 25,000 years in the southern Middle Atlas, Morocco[END_REF][START_REF] Bouarakia | Notes on the distribution and phylogeography of two rare small Gerbillinae (Rodentia, Muridae) in Morocco: Gerbillus simoni and Gerbillus henleyi[END_REF][START_REF] Bouarakia | First genetic identification of Sundevall's jird Meriones crassus (Rodentia, Muridae) in Morocco and the hypothesis of range extension[END_REF]. Similar to the case of G. henleyi, Sundevall's Jird Meriones crassus Sundevall, 1842, a gerbilline species widespread in the desert regions of North Africa and Asia, was also recently discovered in eastern Morocco [START_REF] Bouarakia | First genetic identification of Sundevall's jird Meriones crassus (Rodentia, Muridae) in Morocco and the hypothesis of range extension[END_REF].

Taphonomic observations

For a study of digestion, a total of 357 Gerbil fossil remains were considered (114 femoral heads, 136 incisors and 107 M1 and m1 molars). When considering the entire material, the analysis of the bone and dental surface under a stereomicroscope and with a scanning electron microscope (SEM) shows a low percentage of digested elements (13.2% of the femoral heads, 16.2% of the incisors and 7.5% of the molars) (Table 4), with mostly low to moderate digestion grades (Figures 7 and8).

While the exact identification of the raptor species responsible for the accumulation of the assemblage is always a difficult exercise, it is nevertheless possible to determine a predation category comprising a number of potential predators with similarities in their taphonomic signals. In Guenfouda, the alteration produced by digestion highlights the probable role of category 1 or 2 nocturnal predator(s) according to the classification by [START_REF] Andrews | Owls, caves & fossils[END_REF][START_REF] Fernández-Jalvo | Taphonomy for taxonomists: implications of predation in small mammal studies[END_REF]. These categories include, for example, the barn owl Tyto alba, an accumulator of small vertebrates which is also present in the region. Note that the presence of a few remains with stronger alterations may come from the occasional contribution of other predators, or young barn owls which tend to digest more strongly than adults [START_REF] Bruderer | Inventaire taxonomique et taphonomique d'un assemblage de pelotes d'un site de nidification de Tyto alba de Mauritanie[END_REF].

Concerning the post-depositional processes, we observed traces of plant roots which are relatively few, and traces of corrosion caused by soil acidity. We also noted some rodent burrows especially at the level of layers C1a and C1b, which testify a very local disturbance of the deposits, but these were well identified during the excavations and the archeo-stratigraphic context seems to be well preserved.

The slight alteration of the tooth and bone elements suggests that their accumulation is probably linked to a category 1 or 2 nocturnal raptor, possibly Tyto alba. Such predators are known to be opportunistic hunters and to give a good representation of the diversity of the preyspecies in the environment, in the limit of their size, nocturnal habits and hunting territory (e.g., [START_REF] Andrews | Owls, caves & fossils[END_REF]).

Paleoenvironments

The presence of gerbils such as G. campestris, and especially G. henleyi, reflects open and arid environments in the vicinity of the cave during the period in which they were accumulated. But the use of only gerbils to reconstruct palaeoenvironments is too reductive, and we need to enlarge the spectrum of the considered taxa.

In the previous study of [START_REF] López-García | The small mammals from the Holocene site of Guenfouda (Jerada, Eastern Morocco): chronological and paleoecological implications[END_REF] on Guenfouda small mammals, the authors used the method of Taxonomic Habitat Index (THI). We adjusted these results in replacing Ctenodactylus sp. by Jaculus orientalis, Meriones shawii by Meriones grandis, and in adding Gerbillus henleyi to the faunal list. The taxonomic name of some other species has been updated as well. We used the habitat categories established by [START_REF] Stoetzel | Implication de l'étude taphonomique des microvetébrés de Témara[END_REF] and [START_REF] Saidani | Nouvelles données sur le contexte taphonomique et environnemental des occupations néolithiques de la grotte de Gueldaman GLD1 (Algérie) d'après l'étude des microvertébrés[END_REF] in including the whole microvertebrate faunal spectrum, and not only micromammals (Table 5). However, Gekkonidae sp., Lacertidae sp. and Chalcides sp. were not included, due to the large differences in the ecology of the species belonging to these groups. We also used the Climatogram method, by weighting the habitats frequented by the species by their abundance (reflected by the MNI = Minimum Number of Individuals) and using the ecological groups defined by [START_REF] Saidani | Nouvelles données sur le contexte taphonomique et environnemental des occupations néolithiques de la grotte de Gueldaman GLD1 (Algérie) d'après l'étude des microvertébrés[END_REF] (Table 6). This method, however, assigns each species to a single habitat, and does not take into account their ecological plasticity. Moreover, it may be biased by the preferred prey and hunting territory of the predator responsible for the accumulation. Nevertheless, this method provides an interesting and complementary insight to the THI.

The results obtained by both methods (Figure 9) indicate that during the early-middle Holocene the environment of the Guenfouda cave was dominated by a semi-open steppic landscape, with patches of more wooded areas and water ponds. Considering the climatogram method, the open habitats (meadows, steppes) are much more dominant than in the THI, mainly due the abundance of Mus spretus in the assemblage. But this species can also frequent ‗intermediate' habitats and open forests. All the identified species still occur today in the region.

Although the presence of Natrix maura and amphibian species indicate the punctual presence of water points, it should be mentioned that only Pelophylax saharicus strictly needs permanent water; Sclerophrys mauritanicus, Bufotes boulengeri and even Discoglossus pictus (to a lesser extent) are more terrestrial and can burrow themselves during the dry season [START_REF] Bons | Amphibiens et reptiles du Maroc (Sahara occidental compris): atlas biogéographique[END_REF][START_REF] Escoriza | Amphibians of North Africa[END_REF]. Rocky soils were dominant, but soft and sandy soils were also represented, notably through the presence of burrowing species, and especially G. henleyi. Our results are globally in line with those previously obtained by [START_REF] López-García | The small mammals from the Holocene site of Guenfouda (Jerada, Eastern Morocco): chronological and paleoecological implications[END_REF], with however a lower representation of open-dry and rocky habitats with the THI method, mainly due to the actual absence of Ctenodactylus in the assemblage. Comparative analyses with other contemporaneous (Holocene) sites are difficult, because small vertebrate analyses are scarce in eastern Morocco, and there exist important regional differences in the Maghreb, in particular in Morocco [START_REF] Stoetzel | Quaternary terrestrial microvertebrates from Mediterranean northwestern Africa: state-of-the-art focused on recent multidisciplinary studies[END_REF]. The results obtained for the Holocene levels of Guenfouda indicate an environment close to what is known today in the region, i.e., an arid Mediterranean climate, and a landscape dominated by steppes and scrubs (especially with the climatogram method), together with patches of more wooded areas and water ponds or streams.

Conclusions

This work contributes to the knowledge of the diversity, geographical and chronological occurrence of gerbil species in North Africa. A total of 405 fossil remains were analysed and measured (114 femoral heads, 136 incisors and 107 molars). The morphological and biometric analysis, as well as the comparison with current and fossil specimens, allowed us to attribute the fossil material of the Guenfouda gerbils to G. campestris and G. henleyi. The latter is cited in a fossil state for the first time in North Africa.

The taphonomic analysis shows that a nocturnal raptor (probably Tyto alba) was the source of the assemblage. Some burrows observed during excavations argue for few local disturbances of the deposits, without major mixing of sediments or archaeological material.

The abundance of gerbillines in the Guenfouda assemblage indicates an open and dry environment, which is also characteristic for this region at present.

Figure 1 .

 1 Figure 1. a) Location of the Guenfouda Cave and representation of the main archaeological sites of Eastern Morocco; b) Close-up view of the cave, carved out of the Jurassic limestone (Photo Aouraghe. H); c) Panoramic view from the interior of the cave and the strategic position of the place (Photo Aouraghe. H).; d) Plan of the excavation of Guenfouda cave (Aouraghe et al. 2014).

Figure 2 .

 2 Figure 2. Stratigraphic section of the Guenfouda cave (square J40) (picture and infography Hamid Haddoumi).

Figure 3 .

 3 Figure 3. Examples of first molars of G. campestris from Genfouda. a = M1L (G054-I39); b = m1R (G123-K39); c = m1R (G115-L39) displaying anterior and posterior cingulum; d = upper (left) and lower (right) dental rows of a modern G. campestris from Fritissa (Morocco), specimen (MAR15-FRT17) trapped and genotyped by O. Bouarakia in 2015 (drawing H. Mhamdi). Scale = 1 mm. Dental nomenclature of the Gerbillinae (following Tong 1989): ac = anterocone; acd = anteroconid; al = anterolophule; ald = anterolophulid; ecd = entoconid; hc = hypocone; hcd = hypoconid; mc = metacone; mcd = metaconid; nl = neolophe; nld = neolophulid; pa = paracone; prd = protoconid; pr = protocone.

Figure 4 .

 4 Figure 4. Examples of first molars of G. henleyi from Guenfouda. a = M1R (G055-I39); b = m1L (G125-K39); c = upper (left) and lower (right) dental rows of a modern G. henleyi from Ain Beni Mathar (Morocco), specimen (MAR14-BMT1) trapped and genotyped by O. Bouarakia in 2014 (drawing H. Mhamdi). Scale = 1 mm.

Figure 5 .

 5 Figure 5. Comparison of mean size of upper and lower first molars of different fossil and modern Gerbillus species, standard deviation in brackets. n = number; MNHN = Muséum National d'Histoire Naturelle, Paris, France; FSR = Faculty of Sciences of Rabat, Morocco; BMNH = British Museum of Natural History, London, UK.

Figure 6 .

 6 Figure 6. Bivariate diagram of measurements (length x width in mm) of upper (M1) and lower (m1) first molars of fossil and modern Gerbillus species. n = number; MNHN = Muséum National d'Histoire Naturelle, Paris, France; FSR = Faculty of Sciences of Rabat, Morocco; BMNH = British Museum of Natural History, London, UK.

Figure 7 .

 7 Figure 7. Examples of rodent digested incisors from Guenfouda. a: non digested gerbil in situ lower incisor; b: lightly digested gerbil lower incisor (enamel surface pitting on the tip); c-d: moderately digested rodent upper incisor (enamel and dentine dissolved on the tip); e-f: heavily digested rodent upper incisor (enamel removed from most of the tooth, dentine strongly altered).

Figure 8 .

 8 Figure 8. Examples of digested Gerbillus femoral heads from Guenfouda. a: Femur head without any digestion; b: light digestion (slight penetration of the head); c: moderate digestion (advanced penetration of the whole head with osteocyte lacunae well visible). Scale = 1 mm.

Figure 9 .

 9 Figure 9. Representation of the THI (a) and climatogram (b) results based on the microvertebrate species from the Holocene levels of Guenfouda.

Table 4 .

 4 Percentages and grades of digestion of Gerbillus species incisors, molars, and femoral heads from Guenfouda. PD = Percentage of Digestion %.

			Bush -	Open
		Forest	open	(meadow,	Wetland	Rocky	Sandy
			forest	steppe)
	Chiroptera
	Myotis punicus	0.00	0.60	0.30	0.00	0.10	0.00
	Rhinolophus hipposideros	0.40	0.40	0.10	0.00	0.10	0.00
	Macroscelidea
	Petrosaltator rozeti	0.00	0.00	0.30	0.00	0.50	0.20
	Soricomorpha
	Crocidura russula	0.20	0.40	0.30	0.05	0.05	0.00
	Crocidura whitakeri	0.10	0.30	0.40	0.05	0.10	0.05
	Rodentia
	Eliomys munbyanus	0.30	0.30	0.10	0.00	0.30	0.00
	Apodemus sylvaticus	0.60	0.30	0.05	0.05	0.00	0.00
	Mus spretus	0.10	0.30	0.40	0.10	0.05	0.05
	Gerbillus campestris	0.05	0.20	0.55	0.00	0.10	0.10
	Gerbillus henleyi	0.00	0.00	0.60	0.00	0.10	0.30
	Meriones grandis	0.00	0.20	0.60	0.00	0.10	0.10
	Atlantoxerus getulus	0.00	0.10	0.20	0.00	0.70	0.00
	Jaculus orientalis	0.00	0.00	0.60	0.00	0.30	0.10
	Amphibia
	Sclerophrys mauritanica	0.05	0.15	0.30	0.30	0.15	0.05
	Bufotes boulengeri	0.20	0.20	0.20	0.30	0.05	0.05
	Pelophylax saharicus	0.20	0.20	0.20	0.40	0.00	0.00
	Discoglossus pictus	0.20	0.20	0.20	0.30	0.05	0.05
	Squamata
	Trogonophis wiegmanni	0.15	0.35	0.20	0.00	0.10	0.20
	Agama bibronii	0.00	0.00	0.30	0.00	0.70	0.00
	Chamaeleo chamaeleon	0.10	0.80	0.10	0.00	0.00	0.00
	Eumeces algeriensis	0.05	0.25	0.40	0.00	0.15	0.15
	Malpolon monspessulanus	0.05	0.20	0.50	0.05	0.15	0.05
	Coronella girondica	0.10	0.50	0.30	0.05	0.05	0.00
	Natrix maura	0.10	0.20	0.20	0.50	0.00	0.00
	Total	2.95	6.15	7.40	2.15	3.90	1.45

Table 5 .

 5 Values of the Taxonomic Habitat Index (THI) of the microvertebrate species from the Holocene levels of Guenfouda. Adapted and completed from López-García et al. (2013),[START_REF] Stoetzel | Implication de l'étude taphonomique des microvetébrés de Témara[END_REF] and[START_REF] Saidani | Nouvelles données sur le contexte taphonomique et environnemental des occupations néolithiques de la grotte de Gueldaman GLD1 (Algérie) d'après l'étude des microvertébrés[END_REF].

Acknowledgments

The Faculty of Sciences of Mohamed First University at Oujda provided us with the logistical means to carry out excavations in the Guenfouda cave; our thanks go to the excavation team. 14C dating was performed on charcoal remains at the Centre de Datation par le Radiocarbon, Lyon, France (Lyon-15389-SacA-53713). We also thank the technical team of the Scanning Electron Microscope (SEM) of the Faculty of Sciences of Oujda for helping us to take highresolution pictures, and the -Biodiversity, Ecology and Genome‖ Laboratory of the Faculty of Sciences of Rabat, for allowing us to use its collections of current specimens. We warmly thank Salvador Bailon (MNHN, Paris) for providing us information about the identification and counts of amphibians and squamates of Guenfouda. We thank also Pr Hamid Hadoumi for his contribution to the description and illustration of the stratigraphic sequence, and Mohamed Amezian for the English revision. Finally, we would like to thank the two anonymous reviewers for their valuable corrections and comments that allowed to improve the manuscript.

Data availability statement:

The authors confirm that the data supporting the findings of this study are available within the article.

Disclosure statement:

No potential conflict of interest was reported by the author(s). 1. Taxonomic list of the large and small vertebrates of Guenfouda cave [START_REF] Aouraghe | La grotte préhistorique de Guenfouda, Maroc Oriental : les premiers résultats (fouilles 2004-2007)[END_REF][START_REF] Aouraghe | Mission Guenfouda (Maroc oriental) Rapport sur la campagne de Fouille[END_REF][START_REF] Bougariane | Les Vertébrés du Pléistocène terminal-Holocène de quelques sites marocains : Paléontologie, Taphonomie et Archéozoologie[END_REF][START_REF] López-García | The small mammals from the Holocene site of Guenfouda (Jerada, Eastern Morocco): chronological and paleoecological implications[END_REF]Stoetzel et al. 2017, present work). The names of some amphibian species have been updated according to recent works in taxonomy and systematics [START_REF] Escoriza | Amphibians of North Africa[END_REF]. + = presence, -= absence, HS = not in stratigraphy.

Groups


R Adatorwovor
email: radatorwovor@uky.edu

A Latouche

A Parametric Approach to Relaxing the Independence Assumption in Relative Survival Analysis

Keywords: Competing risks, Copula, Dependence Modelling, Net survival, Relative survival 1 Introduction

With known cause of death (CoD), competing risk survival methods are applicable in estimating disease-specific survival. Relative survival analysis may be used to estimate disease-specific survival when cause of death is either unknown or subject to misspecification and not reliable for practical usage. This method is popular for population-based cancer survival studies using registry data and does not require CoD information. The standard estimator is the ratio of all-cause survival in the cancer cohort group to the known expected survival from a general reference population. Diseasespecific death competes with other causes of mortality, potentially creating dependence among the CoD. The standard ratio estimate is only valid when death from disease and death from other causes are independent. To relax the independence assumption, we formulate dependence using a copulabased model. Likelihood-based parametric method is used to fit the distribution of disease-specific death without CoD information, where the copula is assumed known and the distribution of other cause of mortality is derived from the reference population. We propose a sensitivity analysis, where the analysis is conducted across a range of assumed dependence structures. We demonstrate the utility of our method through simulation studies and an application to French breast cancer data.

K = 2 and ε = 1 implies death from cancer and ε = 2 implies death from other competing causes. Standard methods for independently right censored survival data without competing risks cannot generally be used to make inference about disease-specific survival. Under dependent competing risks, where T1 and T2 are dependent, the Kaplan-Meier [START_REF] Kaplan | Nonparametric estimation from incomplete observations[END_REF] curve estimates a function of the cause-specific hazard function, defined in Section 2. The logrank test [START_REF] Bland | The logrank test[END_REF]) assesses group differences between the cause-specific hazard function, while the standard proportional hazards model [START_REF] Cox | Regression models and life tables (with discussion)[END_REF] formulates the effects of covariates on the cause-specific hazard function. The cumulative incidence function, defined in Section 2, gives diseasespecific survival in the presence of competing events. This quantity has been widely adopted in applications, with the Aalen-Johanson estimator [START_REF] Aalen | An empirical transition matrix for non-homogeneous Markov chains based on censored observations[END_REF]), Gray's test [START_REF] Gray | A class of K-sample tests for comparing the cumulative incidence of a competing risk[END_REF], and the Fine-Gray model [START_REF] Fine | A proportional hazards model for the subdistribution of a competing risk[END_REF], providing analogs to the Kaplan-Meier curve, the logrank test, and the proportional hazards model for the cumulative incidence function. Without cause of death information, these methods are not applicable.

To address disease-specific survival without reliable cause of death information, relative survival methods have been proposed. Relative survival, SR(t), is the ratio of the observed survival rate in a group of cancer patients, during a specified period, to the expected survival rate in a general reference population, [START_REF] Ederer | The relative survival rate: a statistical methodology[END_REF]. Mathematically,

SR(t) = SO(t) SP (t) (1)
where at time t, SO(t) is the observed survival probability from the registry and SP (t) is the expected survival from mortality tables. Existing literature has focused exclusively on the estimation of SR(t) under the independence assumption, T1 ⊥ T2. Under independence, SO(t) = ST 1 (t) • ST 2 (t), with SP (t) = ST 2 (t) which implies SR(t) = ST 1 (t) where ST 1 (t) and ST 2 (t) are the survival probabilities corresponding to T1 and T2 respectively. The relationship (1) can be rewritten in terms of hazard functions as λO(t) = λE(t) + λP (t) [START_REF] Cronin | Cumulative cause-specific mortality for cancer patients in the presence of other causes: a crude analogue of relative survival[END_REF], where λO(t) is the hazard in the disease registry, λE(t) is the so called excess hazard among the cancer cohort, and λP (t) is the hazard from mortality tables. Under independence, λE(t) = λT 1 (t) and λP (t) = λT 2 (t), where λT j (t) = -dlogS T j (t) dt , j = 1, 2, are the net hazard functions for cancer and other cause mortality, respectively. The net survival probability (under independence assumption) ST 1 (t) is the target of relative survival analysis and corresponds to a hypothetical population in which death from competing causes does not exist. It differs from the cumulative incidence function which is commonly used to quantify disease-specific survival in analyses with known cause of death information. Under dependence, SR(t) in ( 1) has an excess hazard [START_REF] Suissa | Relative excess risk: an alternative measure of comparative risk[END_REF] interpretation since the estimator eliminates the effect of background risk, and is no longer the survival probability ST 1 (t).

Relative survival method was pioneered by [START_REF] Berkson | Calculation of survival rates for cancer[END_REF], and [START_REF] Ederer | The relative survival rate: a statistical methodology[END_REF] for nonparametric estimation of ST 1 (t). A variant of this method was proposed by [START_REF] Hakulinen | Cancer survival corrected for heterogeneity in patient withdrawal[END_REF] to address the bias due to heterogeneity of patient withdrawal within subgroups. Pohar [START_REF] Perme | On estimation in relative survival[END_REF] demonstrated that these classical methods may be biased under certain censoring patterns. For example, in population comparisons, such bias may arise from unmeasured covariates affecting the cancer cohort group and the reference population from which rates of expected mortality are drawn. Rebolj Kodre and Pohar Perme, (2013) studied biases associated with censoring and age distribution (at the time of cancer diagnosis) and proposed weighting corrections. [START_REF] Nixon | Relationship of patient age to pathologic features of the tumor and prognosis for patients with stage I or II breast cancer[END_REF] documented that event times and censoring times are dependent on the age of the patients in a cancer study. Stratified methods [START_REF] Sasieni | On standardized relative survival[END_REF] based on age standardization of relative survival ratios may reduce such biases. [START_REF] Hakulinen | Choosing the relative survival method for cancer survival estimation[END_REF] developed alternative estimators valid under weaker assumption.

Several authors [START_REF] Bolard | Assessing time-by-covariate interactions in relative survival models using restrictive cubic spline functions[END_REF][START_REF] Giorgi | A relative survival regression model using B-spline functions to model non-proportional hazards[END_REF][START_REF] Nelson | Flexible parametric models for relative survival, with application in coronary heart disease[END_REF][START_REF] Mahboubi | Flexible modeling of the effects of continuous prognostic factors in relative survival[END_REF] proposed variations of spline functions to address some of limitations of earlier methods. These models include mainly flexible parametric functions that examine the effects of covariates that potentially influence the estimation of the excess survival. [START_REF] Charvat | A multilevel excess hazard model to estimate net survival on hierarchical data allowing for non-linear and non-proportional effects of covariates[END_REF] introduced spline models for grouped data, where they fitted flexible excess hazard models to cluster data, while [START_REF] Rubio | On a general structure for hazard-based regression models: an application to population-based cancer research[END_REF] proposed extensions (beyond proportional hazard method) to the general parametrization of hazard functions and implemented flexible parametric distribution in modelling the event time. However, the above estimation methods for SR(t) do not model the dependence event times concurrently.

To relax the independence assumption, we formulate the dependence between the latent failure times distributions for death from disease and death from competing causes using copula models [START_REF] Deheuvels | Caractérisation complète des lois extrêmes multivariées et de la convergence des types extrêmes[END_REF]. A copula function generates a joint distribution for the two event times, taking as input their marginal distributions. Copulas allow a broad range of dependence structures and have been employed widely in survival analysis, including bivariate event times [START_REF] Oakes | A model for association in bivariate survival data[END_REF], competing risks with known cause of failure [START_REF] Heckman | The identifiability of the competing risks model[END_REF], and semi-competing risks where one event time censors the other but not vice versa [START_REF] Fine | On semi-competing risks data[END_REF]. We employ such models with competing risks data from disease registries where cause of death information is either not reliable or not available. Because the joint distribution of the latent failure times is nonparametrically nonidentifiable [START_REF] Tsiatis | A nonidentifiability aspect of the problem of competing risks[END_REF], we treat the copula function as known. The marginal distribution of the time to disease-specific death is modelled parametrically with the distribution of death from other causes drawn from the reference population. Likelihood-based inference is proposed. Because the joint distribution is unidentifiable nonparametrically and unverifiable from the observed registry data, a sensitivity analysis is suggested in which disease-specific survival is estimated across a range of rich dependence structures, specified via the copula func-tion. To our knowledge, this is the first attempt in modelling dependence between T1 and T2 in relative survival analysis.

The main purpose of this method is to provide an alternative estimator for net survival (survival in a hypothetical world where other competing causes of death do not exist) under dependence between cancer mortality and other cause mortality. This is accomplished by explicitly modelling the dependence between cancer mortality and other cause mortality. The rest of this paper proceeds as follows. In section 2, we present the data and copula model formulation for competing risks data. Section 3 describes the likelihood estimation and inference procedure without cause of death information, as well as the proposed sensitivity analysis. In section 4, we present the numerical illustrations including simulation results and application to French breast cancer data. Section 5 discusses and concludes the paper.

Data and Model

We begin by defining traditional endpoints for competing risk data with known cause of death. The cause-specific hazard, λ k (t) is the instantaneous failure rate for occurrence of event ε = k at time t [START_REF] Prentice | The analysis of failure times in the presence of competing risks[END_REF],

λ k (t) = lim δt→0 P (t ≤ T < t + δt, K = k|T > t) δt (2)
and the cumulative incidence function C k (t) is the proportion of patients who died from cause k by time t in the presence of patients who might die from other causes. The cumulative incidence function can be expressed as

C k (t) = P (T ≤ t, ε = k) = t 0 λ k (s) • S(s)
ds where S(t) = P (T > t) is the overall survival probability. Standard competing risks methods with known cause of failure focus on estimation of λ k (t) and C k (t).

Without cause of death information, the registry data is simply time to death from any cause, T, which may be right censored by lost to follow up. Let C be the time to right censoring, with the common assumption being that T and C are independent. The observed data consist of Xi = min(Ti, Ci) and δi = I(Ti ≤ Ci), where Ti and Ci are the failure and censoring times on individual i = 1, 2, 3, • • • , n. Relative survival methods employing such data do not focus on the traditional competing risks endpoints λ k (t) and C k (t) but rather on the latent failure time distributions with the corresponding survival functions ST 1 (t) and ST 2 (t).

To capture the dependence between T1 and T2, we employ copula models, which completely describe the dependence structure and provide scale invariant measures of association [START_REF] Venter | Tails of copulas[END_REF][START_REF] Müller | Orderings of risks: A comparative study via stop-loss transforms[END_REF][START_REF] Bäuerle | Modeling and comparing dependencies in multivariate risk portfolios[END_REF][START_REF] Denuit | On s-convex stochastic extrema for arithmetic risks[END_REF]. Suppose ψ is a copula generator function defined such that ψ : [0, 1] → [0, +∞], then the j-dimensional copula function is given by C(u1, • • • , uj) = ψ ψ -1 (u1), • • • , ψ -1 (uj) with marginal distributions, uj = P (Tj ≤ tj) = FT j (tj) = 1 -ST j (tj), ∀j ∈ N. When j = 2, then, the copula model for the joint distributions of T1 and T2 is:

C(u1, u2) = P (T1 ≤ t1, T2 ≤ t2) = ψ ψ -1 (u1) + ψ -1 (u2) = FT 1 ,T 2 (t1, t2)
where ψ -1 is the inverse of ψ and ψ satisfies the Laplace-Stiltjes transform and [START_REF] Bernstein | Sur les fonctions absolument monotones[END_REF] theorem. [START_REF] Mcneil | Multivariate Archimedean copulas, d-monotone functions and l-norm symmetric distributions[END_REF] showed that the generator function ψ is completely monotone for non-negative random variables with ψ(0) = 1, ψ (•) < 0 and ψ (•) < 0.

The most widely used scale invariant measures of association to characterize dependence are Spearman's rho (ρS) and Kendall's tau (τ ken ) correlation coefficients. The connection between the latter and the copula generator function was shown by Genest and MacKay, (1986) as:

τ ken = 1 + 4 1 0 ψ -1 (u) ψ -1 (u) du = 1 -4 ∞ 0 u(ψ(u)) 2 du
with ψ -1 being the derivative of ψ -1 . While in theory, any copula may be used to link the marginal distributions of T1 and T2, in this paper, we focus on two popular Archimedean copulas, indexed by a single dependence parameter θ having simple interpretations. The Gumbel copula is:

C(u1, u2) = exp -{(-log(u1)) θ + (-log(u2)) θ } 1 θ (3)
with θ ∈ (1, +∞) and the Clayton copula is:

C(u1, u2) = (u -θ 1 + u -θ 2 -1) -1 θ (4)
with θ ∈ (0, +∞). A special case of product copula model: C(u1, u2) = u1 • u2 is obtained when θ = 1 and when θ → 0 for Gumbel and Clayton copulas respectively. The product copula model gives independence of T1 and T2. When θ > 0, the Clayton copula is bounded by: C(u1, u2) ≤ θ(1 -u1 -u2) + (1 + θ)u1u2. As dependence increases, that is θ → +∞, the Clayton copula approximates the Frećhet-Hoeffding upper bound, [START_REF] Fréchet | Sur les tableaux de corrélation dont les marges sont données[END_REF][START_REF] Hoeffding | Masstabinvariante korrelationstheorie[END_REF] giving perfect positive dependence.

We first formulate our model without covariates for the potentially dependent latent failure times T1 and T2. In the sequel, the distribution of T1 involves an unknown parameter η, while that of T2 is assumed from a reference population. The distribution of T implicitly involves η via T = min(T1, T2). The survival function for all-cause mortality at time t given η is:

ST (t|η) = ST 1 (t|η) + ST 2 (t) -1 + FT 1 ,T 2 (t, t|η) = 1 -FT 1 (t|η) -FT 2 (t) + FT 1 ,T 2 (t, t|η) (5)
with the corresponding density function using

FT (t|η) = FT 1 (t|η) + FT 2 (t) -FT 1 ,T 2 (t, t|η)) fT (t|η) = fT 1 (t|η) + fT 2 (t) -fT 1 ,T 2 (t, t|η) (6)
where

fT 1 (t|η) = dF T j (t|η) dt , fT 1 ,T 2 (t, t|η) = dF T 1 ,T 2 (t,t|η) dt
and fT 2 (t) is derived from the reference population. If censoring of T by C is noninformative, then the likelihood contribution for individual i is:

Li = fX i ,∆ i (Xi, δi|η) = [fT (Xi|η)] δ i [ST (Xi|η)] 1-δ i (7)
From equation ( 7), the full log-likelihood function based on n independent observations is:

l(X, ∆|η) = n i=1 (δi * log fT (Xi|η) + (1 -δi) * log ST (Xi|η)) = n i=1 δi * log [fT 1 (Xi|η) + fT 2 (Xi) -fT 1 ,T 2 (Xi, Xi|η)] + n i=1 (1 -δi) * log [ST 1 (Xi|η) + ST 2 (Xi) -1 + FT 1 ,T 2 (Xi, Xi|η)] (8)
where

(X, ∆) = (Xi, ∆i, i = 1, 2, 3, • • • , n).
We specify a parametric model for FT 1 (t), with finite dimensional parameter of interest η. The general form of the probability density function for T1 at time t is fT 1 (t|η) with survival probability

ST 1 (t|η) = 1 - FT 1 (t|η) = ∞ t fT 1 (s|η)ds.
The distribution of T2 is assumed known and extracted from the reference population with the usual assumption that disease-specific death is negligible in this reference population [START_REF] Ederer | The relative survival rate: a statistical methodology[END_REF]. This is illustrated in the French breast cancer data analysis in section 4.2. The copula distribution linking FT 1 (t|η) and FT 2 (t) may be specified using simple parametric copula models such as the Archemedean copulas. The parameters in the copula model may be chosen for a pre-specified dependence between T1 and T2, for example, Kendall's tau (τ ken ). In the numerical illustrations, T1 was assumed to follow an exponentiated Weibull distribution with parameter η = (λ, κ, α) and probability density function

f (t : λ, κ, α) = α κ λ • t λ κ-1 • exp -t λ κ 1 -exp -t λ κ
α-1 because of its versatility to accommodate a wide range of hazard shapes.

We consider the Gumbel and Clayton copulas for the joint distribution of T1 and T2 as both copulas exhibit tail behaviours that mimic the mortality trend observed in the cancer registry data. The bivariate joint distribution and density functions under the Gumbel copula are given below.

FT 1 ,T 2 (t|η, t) = exp -(-log (FT 1 (t|η))) θ + (-log(FT 2 (t))) θ 1 θ fT 1 ,T 2 (t|η, t) = FT 1 ,T 2 (t|η, t) • -log (FT 1 (t|η)) θ + -log (FT 2 (t)) θ 1 θ -1 × -log(FT 1 (t|η)) θ-1 • fT 1 (t|η) FT 1 (t|η) + -log (FT 2 (t|η)) θ-1 • fT 2 (t) FT 2 (t) (9)
Under the Clayton copula, the bivariate joint distribution and density functions are:

FT 1 ,T 2 (t|η, t) = FT 1 (t|η) -θ + FT 2 (t) -θ -1 -1 θ fT 1 ,T 2 (t|η, t) = FT 1 ,T 2 (t|η, t) (FT 1 (t|η) -θ + FT 2 (t) -θ -1) • fT 1 (t|η) FT 1 (t|η) θ+1 + fT 2 (t) FT 2 (t) θ+1 (10)
The maximum likelihood estimator (MLE) of η can be obtained by maximizing the log-likelihood function in (8) using any optimization algorithm such as Nelder-Mead [START_REF] Nelder | A simplex method for function minimization[END_REF] or better. In the simulation in section 4.1, and because the model is highly nonlinear for hazard shapes and the estimation of dependent parameter θ, we suggest using multiple starting values and taking the MLE to be the maximizer giving the largest value of the log likelihood across all starting values. Under a fixed and correctly specified copula, the usual regularity conditions for the MLE holds and the estimator converges in probability, that is η P -→ η and is asymptotically normal, η ∼ N η, IO(η) -1 with variance estimated using the inverse of the observed information matrix

IO(η) = - ∂ 2 l(X, ∆|η) ∂η∂η T | η=η = -δi • n i=1      [fT (Xi|η)] • ∂ 2 ∂η 2 fT (Xi|η) -∂ ∂η fT (Xi|η) T ∂ ∂η fT (Xi|η) [fT (Xi|η)] 2      +(1 -δi) • n i=1      [ST (Xi|η)] • ∂ 2 ∂η 2 ST (Xi|η) -∂ ∂η ST (Xi|η) T ∂ ∂η ST (Xi|η) [ST (Xi|η)] 2      | η=η (11)
Since the dependence structure for time to disease-specific mortality (T1) and time to other competing mortality (T2) is nonidentifiable and unverifiable from the observed registry data, we propose a sensitivity analysis, where the analysis is conducted across a range of assumed dependence structures. The levels of dependence represent the varying levels of dependent competing mortality possible in the observed registry data. For each copula dependence structure with known θ, we estimate η with η and compute FT 1 (t|η) to estimate relative survival. The corresponding standard errors are obtained as the square root of the Delta method variance:

V ar( ST 1 (X)) = g( ST 1 (X)) • IO(η) -1 • g T ( ST 1 (X)) where g(η)
is the derivative of ST 1 (t|η) with respect to η. Due to the complex nature of the likelihood, numerical approximation is used to estimate the information matrix in the numerical illustrations in Section 4.

In the presence of informative censorship where T and C are dependent, we propose conditioning on additional covariates Z in FT 2 , (Sasieni and Brentnall, 2017; and Pohar [START_REF] Perme | On estimation in relative survival[END_REF], where FT 2 (t|Z) is the conditional distribution of T2 given Z. Such covariates might include age, sex, period, as well as other relevant demographic variables. Let Zi be the covariate observed on individual i = 1, • • • , n. The log-likelihood function ( 8) is easily modified, where the likelihood contribution for individual i (= 1, • • • , n) is ( 7) with FT 2 (t|Zi) replacing FT 2 (t) in fT (Xi|η) and ST (Xi|η). Here, we estimate η in FT 1 (t|η) conditionally on Z to mitigate against the bias (caused by informative censoring) associated with these covariates (Pohar Perme et al, 2012; and Sasieni and Brentnall, 2017, [START_REF] Schaffar | Estimation of net survival for cancer patients: relative survival setting more robust to some assumption violations than cause-specific setting, a sensitivity analysis on empirical data[END_REF]. The usual likelihood regularity conditions continue to hold, with the resulting estimator η being consistent and asymptotically normal with variance which may be estimated using the inverse of the observed information matrix evaluated at η.

Numerical Illustrations

Simulation Studies

To evaluate the performance of our proposed method, first, we simulated a general survival time data Tj ∼ ExpW eibull(λj, κj, αj), and second, to mimic the French breast cancer data set for sample sizes; 1000, 2500 and 5000 with 500 replications. The latent failure times for Tj with probability density function defined above in Section 3 and parameters as in figure 1. The exponentiated Weibull distribution degenerates to a 2-parameter Weibull distribution when α = 1. The parameters for the reparametrized Weibull distribution for T1 were λ1 = 0.182 and α1 = 1.609, while those for T2 were λ2 = 0.742 and α2 = 0.693. In the estimation of λ1, α1 for T1, λ2, α2 are assumed known for T2 and vice versa for the estimation of λ2 and α2. Noninformative censoring times were generated from a uniform distribution (0, γ), where γ was chosen for 10, 30 (omitted from table) and 50% censoring. We consider the Gumbel copula with Kendall's tau, τ ken = 1 -1 θ = 0, 0.25, 0.50, and 0.75. Initial parameter values were randomly chosen from uniform distributions, with multiple starting values (wherever possible) as described in Section 3. We also simulated data from the Clayton copula. The results are similar to those for the Gumbel copula and are described in the appendix. Tables 1, 5 and 6 show the results for estimation of the model for T1 treating T2 as a competing event and for T2 treating T1 (for brevity, tables 5 and 6 show results for selected sample sizes) as a competing event. The bias is small decreasing to zero as the sample size increases for each of the censoring levels. The empirical variance and the model based variance tend to agree and the coverage is close to the nominal 0.95 level, particularly at larger sample sizes. The empirical variance decreases as the sample size increases at roughly the expected root n rate. Table 4 shows the robust survival estimates for the mis-specified model (data simulation from exponentiated Weibull but Weibull survival estimated instead) for T1 treating T2 as competing risk. We presented the results for the lower, mean and upper quantiles for each of the sample sizes and dependence levels. This shows decreasing bias for increasing sample sizes for each of the quantiles. 2 and3 show the estimates of ST 1 (t) for cancer mortality both overall and stratified by age. The parametric estimates under independence are similar to those from the Pohar-Perme method. This suggests that the Weibull assumption is a reasonable fit to the data. One observes that as dependence increases, cancer survival generally decreases. For a fixed dependence level, younger women tend to have higher cancer survival rates than do older women, with marked reductions for the 65-74 and 75-99 age groups. There is some instability in survival estimates at 15 years, especially for the older age groups, as evidenced by the large standard errors. This may be due to small numbers of patients at risk at longer follow-up times.

The net survival function for cancer corresponds to a hypothetical world where the only cause of death is breast cancer. This quantity can only be estimated under unverifiable dependence assumptions between T1 and T2 using disease registry data. To account for uncertainty in dependence, we recommend reporting a range of probabilities corresponding to differing levels of dependence. For example, using results from table 2, the overall 5 year net breast cancer survival from 1980 -2011 is estimated to be between 84.0-87.4% under dependence ranging from Kendall's tau equal to 0 (independence) to 0.75 (strong dependence). These cancer survival probabilities may be meaningfully compared with those in other populations having different background mortality rates and different dependence levels between T1 and T2.

The sensitivity analysis was conducted across different levels of dependence representing different competing mortality potentially observable in the registry data. Figures 2 and3 show the 2, 5, 10 and 15-yr overall net breast survival plots across a spectrum of dependence structures for women between the ages of 18 and 96-yr living in France during 2008 and 2011. As the dependence level increases, the net breast cancer survival decreases dramatically. Perhaps, this might be due to increase hazard for the patients compounding the effect of competing mortality hereby decreasing the chances of survival.

Discussion and Conclusion

Our model formulation for competing risk data without cause of failure information is general, permitting arbitrary but known copula functions. The distribution of other cause mortality is obtained from external reference data (Sarfati et al., 2010; Pohar Perme et al., 2012; Sasieni and Brentnall, 2017). We have undertaken preliminary investigations of simultaneous estimation of the dependence parameter and the parameter in the disease-specific survival distribution. There is evidence of instability in the estimation process especially at the boundary values, with care needed in the model specification to aid estimability of the model parameters. This is expected (see [START_REF] Zahl | Frailty modelling for the excess hazard[END_REF] for challenges), as there are similar issues even when the cause of failure is known. The proposed sensitivity analysis is a practical solution to this issue, providing a range of estimates across different dependence levels not requiring simultaneous estimation of the dependence parameter. The parametric model for diseasespecific mortality is restrictive but may be flexible enough for applications where the hazard is smooth over time, which is the case in cancer registry data. Estimating expected hazard or the distribution of T2 from life tables is limiting as mismatches in covariates and other stratifying variables may also induce biases, Rubio et al (2021). However, it is important to note that our estimator is robust and arguably invariant to the parameter estimates under different specifications of the copula models as shown in the parameter estimates for both Gumbel and Clayton copulas (tables 5 and6). To relax the parametric assumption, nonparametric techniques are currently being developed for use in more complex failure patterns.

The focus of relative or net survival analysis is the distribution of the latent event time for death from disease. This endpoint has been advocated by many practitioners [START_REF] Slud | Dependent competing risks and the latent-failure model[END_REF][START_REF] Reason | The contribution of latent human failures to the breakdown of complex systems[END_REF][START_REF] Louzada | The log-Weibull-negative-binomial regression model under latent failure causes and presence of randomized activation schemes[END_REF], as it removes the impact of other cause mortality on the risk of disease-specific mortality, permitting comparisons across populations with different background mortality. As an alternative, other work has considered estimation of the crude disease-specific survival, C k (t), using the relative survival estimates and the known reference hazard for other cause mortality [START_REF] Cronin | Cumulative cause-specific mortality for cancer patients in the presence of other causes: a crude analogue of relative survival[END_REF]. An analogous procedure could be implemented using our copula based estimate of the distribution of T1 and would provide an assessment of the sensitivity of the estimator of C k under independence of T1 and T2. Such procedure would be of interest to individuals who prefer crude disease-specific mortality to net disease-specific mortality. This and time-dependent models are topics for future research.

Our proposed estimator performed well overall and by subgroup analysis. We observed that in cases of elderly patients, long-term survival decreases dramatically as expected. Perhaps this might be due to elderly patients experiencing higher expected mortality rates than younger patients particularly in terms of higher risk of death from other competing causes leading to loss of patients at longer follow-up times. This loss of information induces higher variability in net survival estimates for this elderly populations, thereby inducing a higher variability as observed in the variance. Net survival under dependence competing risk assumption is observable and does not require additive model as in the case for net survival (in hypothetical world) under the independence assumption. Our estimator modelled both the dependence between times to disease-specific event and competing risk event with covariates like age, sex, and period (date of diagnosis and date of event: death or censored). Since it's been known (Nanieli, et. al., 2012) that such covariates affect both excess and expected hazards estimates. We conditioned on these covariates (Pohar-Perme (2012) by matching cancer cohort data with registry data to alleviate biases associated with informative censoring as done in multivariable modelling techniques [START_REF] Bolard | Assessing time-by-covariate interactions in relative survival models using restrictive cubic spline functions[END_REF], Giorgi, et. al., (2003), [START_REF] Lambert | Additive and multiplicative covariate regression models for relative survival incorporating fractional polynomials for time-dependent effects[END_REF], [START_REF] Remontet | An overall strategy based on regression models to estimate relative survival and model the effects of prognostic factors in cancer survival studies[END_REF]) where estimators are adjusted for life-table covariates. Researchers are encouraged to use a patients' medical history in determining the levels of competing risks and use a corresponding dependence survival estimate (independence, low, moderate or high) as a measure for disease-specific prognosis.

In conclusion, our proposed methodology provides estimates for net survival under both independent and dependent competing mortality. On the contrary, Pohar-Perme et al., (2012) estimator is only valid under the independent competing risk assumption. Additionally, Pohar-Perme et al., (2012) estimator may exceed 1 in the left tail. [START_REF] Schaffar | Estimation of net survival for cancer patients: relative survival setting more robust to some assumption violations than cause-specific setting, a sensitivity analysis on empirical data[END_REF] showed that these erratic results may occur with longer follow-up times. Our estimator provides comparable relative or net survival estimates under both independent and dependent competing risk assumptions without the need for cause of disease-specific event in the competing risks setting.

Table 1 :

 1 Estimated parameters for the exponentiated Weibull model for T 1 across samples sizes (N), dependence levels (τ ken ) and for 15% censoring treating T 2 as a competing event. η: estimated parameters: λ, scale,α, and κ are shape parameters, ModB: model-based variance, EMP: empirical variance, CP: 95% coverage probability.

	τ ken	N	η Mean	Bias	ModB EMP	CP
	0.00 1000 λ 2.936 -0.064	0.069	0.056 0.960
			κ 4.198	0.198	8.663	3.873 0.848
			α 0.120	0.020	0.004	0.003 0.942
		2500 λ 2.965 -0.035	0.025	0.023 0.949
			κ 4.180	0.180	2.839	2.305 0.885
			α 0.109	0.009	0.001	0.002 0.894
		5000 λ 2.994 -0.006	0.012	0.012 0.947
			κ 4.264	0.264	1.446	1.250 0.947
			α 0.100 1.6e-4	0.001	0.001 0.931
	0.25 1000 λ 3.001	0.001	0.051	0.062 0.924
			κ 4.912	0.912	3.042 14.274 0.922
			α 0.101	0.001	0.002	0.002 0.896
		2500 λ 3.005	0.005	0.019	0.018 0.960
			κ 4.213	0.213	0.765	0.829 0.958
			α 0.100	0.000	0.000	0.001 0.928
		5000 λ 2.998 -0.002	0.009	0.008 0.960
			κ 4.082	0.082	0.300	0.297 0.958
			α 0.100 1.0e-5 2.0e-4 2.5e-4 0.958
	0.50 1000 λ 3.003	0.003	0.036	0.038 0.954
			κ 4.376	0.376	1.645	1.898 0.952
			α 0.100 3.3e-4	0.001	0.001 0.920
		2500 λ 3.006	0.006	0.014	0.014 0.950
			κ 4.147	0.147	0.439	0.523 0.938
			α 0.100 1.7e-5 2.6e-4	0.001 0.934
		5000 λ 3.002	0.002	0.007	0.007 0.964
			κ 4.072	0.072	0.185	0.187 0.948
			α 0.099 -5.7e-4 1.2e-4 1.4e-4 0.940
	0.75 1000 λ 3.003	0.003	0.027	0.042 0.942
			κ 4.266	0.266	1.101	1.146 0.950
			α 0.102	0.002	0.001	0.007 0.926
		2500 λ 3.006	0.006	0.011	0.010 0.942
			κ 4.132	0.132	0.323	0.330 0.954
			α 0.099 -0.001 2.0e-4 2.0e-4 0.938
		5000 λ 3.004	0.004	0.005	0.005 0.952
			κ 4.066	0.066	0.145	0.145 0.956
			α 0.099 -0.001 9.7e-5 1.1e-4 0.934

Table 2 :

 2 The 2, 5, 10 and 15-yr overall net survival for French women diagnosed with breast cancer between 1980 and 2011. a : ×10 -2 , b : ×10 -3 , τ ken : dependence, PP: Pohar-Perme, S T1 (t): survival estimate at year t, SE: standard error for the relative survival estimate.In this section we analyze data from women between the ages of 18 and 96 years surviving breast cancer in France from 1980 to 2011. The data were obtained from the Institut Curie breast cancer database. This database contains records from 24, 458 nonmetastatic breast cancer patients treated at the Institut Curie. Out of the 24, 458 breast cancer patients, 9, 885 (40.4%) died while 14, 573 were alive and administratively censored on December 31 st 2011. Five age group categories were considered for the estimation of relative survival. 3, 970 were between the ages of15 -44, 6, 895 between the ages of45 -54, 6, 420 between the ages of 55 -64, 4, 675 between the ages of 65 -74 and 2, 498 were in the 75 -99 age group category. We individually matched the observed death or censoring time in the disease cohort group with a corresponding time in the general reference population on age, sex, and year (date of diagnosis and the date of death or censored) for each participant and for each follow-up period. The background mortality data from the Human Mortality Database (https://www.mortality.org) was last modified on June 28, 2018. Within each follow-up year, we assumed that λP (t) is piecewise constant (Dickman et al.

	Independent Competing Risks		Dependent Competing Risks
	τ ken	0.00		0.25		0.50		0.75
	t P P a S T1 (t) a	SE b S T1 (t) a	SE b S T1 (t) a	SE b S T1 (t) a	SE b
	2 95.6	96.0	6.99	95.8	6.96	95.4	7.23	94.7	7.74
	5 84.8	87.4	9.01	86.6	9.10	85.5	9.31	84.0	9.53
	10 71.0	72.8	11.01	71.4 10.99	69.8 10.91	68.0 10.67
	15 59.5	59.5	12.22	57.9 12.08	56.3 11.74	54.9 11.19
	4.2 Application to French Breast Cancer Data

, 2004) for each period up to time X. The cumulative hazard for each period based on λP (t) is calculated from the background survival function at the beginning and end of the period. The cumulative hazard is then used to obtain λP (t) under the piecewise constant assumption. The goal of matching in determining λT 2 = λP is to mitigate the impact of age and calendar year on potentially dependent censoring by C (Pohar Perme et al., 2012). Thus, our approach is stratified by age, sex and year. We estimate 2, 5, 10, and 15-year relative survival assuming a Weibull distribution for T1 and a Gumbel copula model with differing levels of dependence to specify the joint distribution for T1 and T2. We compared our estimates with estimates from Pohar Perme et al. (2012), which require independence of T1 and T2 with ST 2 (t) derived from the background reference population.

Tables

Table 3 :

 3 The 2, 5, 10 and 15-yr age group specific net breast cancer survival among French women diagnosed between 1980 and 2011. a : ×10 -2 , b : ×10 -3 , τ ken : dependence, PP: Pohar-Perme, S T1 (t): relative survival estimate at year t, SE: standard error for the relative survival estimate.

		Independent Competing Risks		Dependent Competing Risks
	τ ken		0.00		0.25		0.50	0.75
	t	Agegp P P a S T1 (t) a	SE b S T1 (t) a	SE b S T1 (t) a	SE b S T1 (t) a	SE b
	2	15-44 95.8	94.9	20.90	94.9 20.73	94.8 20.73	94.8 20.68
		45-54 97.1	96.6	16.44	96.5 16.13	96.3 16.27	96.2 16.40
		55-64 95.7	96.1	13.72	96.0 13.49	95.7 13.70	95.3 14.12
		65-74 95.1	97.0	08.50	96.8 08.54	96.2 09.61	95.1 11.60
		75-99 91.5	96.5	07.94	95.6 08.93	93.4 12.44	89.9 17.16
	5	15-44 85.1	86.9	23.70	86.8 23.64	86.7 23.62	86.7 23.35
		45-54 88.6	90.4	19.39	90.1 19.36	89.8 19.45	89.7 19.28
		55-64 85.8	88.1	17.72	87.6 17.71	86.9 17.87	86.6 17.66
		65-74 84.1	86.9	16.71	85.8 17.01	84.2 17.71	82.5 18.09
		75-99 72.3	77.1	24.21	72.7 24.85	67.1 25.08	61.7 24.00
	10	15-44 71.9	74.4	26.88	74.2 26.84	74.0 26.75	74.1 26.62
		45-54 78.3	80.1	22.83	79.6 22.80	79.2 22.73	79.2 22.34
		55-64 73.4	74.5	22.03	73.5 21.97	72.7 21.74	72.7 21.08
		65-74 68.4	67.2	25.38	65.0 25.32	63.0 24.72	62.3 23.20
		75-99 44.6	43.1	34.83	37.0 32.55	33.0 28.61	31.1 24.35
	15	15-44 62.5	63.2	29.03	63.0 28.96	62.9 28.83	63.0 28.72
		45-54 70.8	70.5	25.31	69.8 25.24	69.4 25.00	69.6 24.51
		55-64 63.5	61.9	24.81	60.7 24.62	59.9 24.11	60.3 23.20
		65-74 50.3	48.7	30.06	46.2 29.46	44.7 27.92	45.3 25.47

Table 4 :

 4 Robust survival probability for a misspecified exponentiated Weibull model for T 1 across samples sizes (N), dependence levels (τ ken ) and for 16% censoring treating T 2 as a competing event. Ŝw : estimated Weibull survival probability, and associated bias at the 0.25, 0.50, and 0.75 quantiles.

			Survival Quantiles	Bias Quantiles
	τ ken	N	sw.25 sw.50 sw.75 Bias.25 Bias.50 Bias.75
	0.00 1000 0.373 0.532 0.751 -19.5e-4 25.4e-3 11.7e-3
		2500 0.372 0.532 0.750 -16.8e-4 26.0e-3 12.3e-3
		5000 0.372 0.532 0.751 -15,1e-4 26.1e-3 12.3e-3
	0.25 1000 0.363 0.533 0.755 -14.7e-3 18.5e-3 68.6e-4
		2500 0.362 0.532 0.754 -14.2e-3 19.1e-3 74.9e-4
		5000 0.362 0.532 0.754 -13.6e-3 19.6e-3 78.0e-4
	0.50 1000 0.345 0.527 0.758 -17.1e-3 17.8e-3 31.2e-4
		2500 0.345 0.527 0.758 -16.7e-3 18.4e-3 36.2e-4
		5000 0.345 0.527 0.758 -16.3e-3 18.7e-3 38.7e-4
	0.75 1000 0.326 0.521 0.761 -74.6e-4 22.5e-3	-5.1e-4
		2500 0.325 0.519 0.761 -64.8e-4 23.4e-3 -1.9e-05
		5000 0.325 0.519 0.761 -61.6e-4 23.8e-3 2.1e-05

Table 5 :

 5 Gumbel Model: Estimated parameters of the Weibull model for T 1 across samples sizes (N), dependence levels (τ ken ) and levels of censoring (C) treating T 2 as a competing event and vice versa.

	η: estimated parameters, ModB: model-based variance, EMP: empirical variance, CP: 95% coverage
	probability. a : ×10 -3 .
		C				0.10			0.50
	τ ken	N	η	Mean Bias a ModB a EMP a	CP	Mean Bias a ModB a EMP a	CP
	0.00 1000 λ1 0.182 -0.080	0.090	0.090	0.940 0.182 -0.290	0.150	0.170 0.928
			α1 1.610	0.790	1.420	1.560	0.938 1.611	1.980	2.620	2.720 0.950
		5000 λ1 0.182 -0.050	0.020	0.020	0.948 0.182	0.050	0.030	0.030 0.958
			α1 1.610	0.520	0.280	0.280	0.954 1.610	0.290	0.520	0.530 0.956
		1000 λ2 0.748	5.980	9.940	10.270 0.936 0.746	3.650	14.410 15.320 0.922
			α2 0.694	0.840	6.790	7.020	0.944 0.697	4.070	8.490	0.010 0.948
		5000 λ2 0.743	0.630	1.870	1.640	0.962 0.743	1.000	2.700	2.460 0.968
			α2 0.693	0.100	1.340	1.210	0.962 0.693	0.280	1.680	1.530 0.958
	0.25 1000 λ1 0.182 -0.480	0.080	0.080	0.948 0.182 -0.450	0.140	0.130 0.954
			α1 1.610	1.490	1.310	1.340	0.952 1.613	3.960	2.480	2.840 0.934
		5000 λ1 0.182 -0.050	0.020	0.020	0.956 0.182 -0.130	0.030	0.030 0.940
			α1 1.609 -0.250	0.260	0.240	0.956 1.610	0.590	0.500	0.460 0.950
		1000 λ2 0.753 10.920 14.700	14.430 0.938 0.760 18.430 20.460 20.810 0.946
			α2 0.690 -3.170	8.240	7.930	0.954 0.687 -5.930 10.080	9.900 0.944
		5000 λ2 0.741 -0.950	2.690	2.530	0.950 0.742	0.260	3.620	3.580 0.964
			α2 0.695	1.940	1.610	1.480	0.958 0.695	2.260	1.960	1.940 0.948
	0.50 1000 λ1 0.182 -0.030	0.070	0.070	0.956 0.182 -0.260	0.120	0.120 0.954
			α1 1.611	2.270	1.270	1.300	0.948 1.613	3.840	2.380	2.710 0.928
		5000 λ1 0.182	0.010	0.010	0.020	0.946 1.824	0.040	0.020	0.030 0.956
			α1 1.609 -0.340	0.250	0.240	0.954 1.610	0.450	0.480	0.510 0.932
		1000 λ2 0.759 17.440 19.080	20.140 0.932 0.767 25.540 25.050 25.330 0.932
			α2 0.688 -5.180	9.440	9.910	0.944 0.684 -9.170 11.210 11.510 0.940
		5000 λ2 0.740 -1.580	3.360	3.380	0.944 0.744	1.620	4.340	4.660 0.946
			α2 0.695	1.720	1.820	1.870	0.936 0.692 -0.750	2.150	2.270 0.944
	0.75 1000 λ1 0.182 -0.200	0.060	0.070	0.956 0.182 -0.260	0.100	0.100 0.948
			α1 1.610	0.490	1.060	1.520	0.936 1.612	2.660	2.090	2.370 0.942
		5000 λ1 0.182	0.050	0.010	0.010	0.948 0.182 -0.020	0.020	0.020 0.952
			α1 1.609 -0.010	0.210	0.210	0.944 1.609 -0.120	0.420	0.450 0.948
		1000 λ2 0.760 17.780 20.190	20.360 0.938 0.766 24.200 26.040 25.790 0.952
			α2 0.689 -4.600	9.870	10.440 0.946 0.685 -7.510 11.580 11.590 0.956
		5000 λ2 0.742 -0.190	3.540	3.770	0.946 0.743	0.830	4.440	4.540 0.950
			α2 0.694	0.550	1.900	1.990	0.930 0.693	0.350	2.210	2.280 0.944

Table 6 :

 6 Clayton Model: Estimated parameters of the Weibull model for T 1 across samples sizes (N), dependence levels (τ ken ) and levels of censoring (C) treating T 2 as a competing event and vice versa. η: estimated parameters, ModB: model-based variance, EMP: empirical variance, CP: 95% coverage probability. a : ×10 -3 .

		C			0.10			0.50
	τ ken	N	η	Mean Bias a ModB a EMP a	CP	Mean Bias a ModB a EMP a	CP
	0.00 1000 λ1 0.182 0.000	0.080	0.090	0.948 0.182 -0.340	0.140	0.160 0.930
			α1 1.610 0.710	1.390	1.520	0.942 1.611 1.820	2.490	2.540 0.948
		5000 λ1 0.182 -0.020	0.020	0.020	0.942 0.182 0.060	0.030	0.030 0.946
			α1 1.610 0.520	0.280	0.280	0.956 1.610 0.710	0.500	0.510 0.952
		1000 λ1 0.747 5.910	9.750	9.910	0.940 0.746 3.570	14.010 14.150 0.922
			α1 0.694 0.780	6.720	6.950	0.948 0.696 3.840	8.360	8.590 0.956
		5000 λ1 0.742 0.690	1.840	1.610	0.962 0.742 0.880	2.640	2.380 0.964
			α1 0.693 0.050	1.330	1.200	0.962 0.693 0.320	1.650	1.490 0.956
	0.25 1000 λ1 0.182 -0.010	0.070	0.080	0.952 0.182 -0.670	0.130	0.110 0.964
			α1 1.610 0.560	1.280	1.410	0.930 1.611 1.870	2.390	0.220 0.952
		5000 λ1 0.182 -0.180	0.010	0.010	0.946 0.182 -0.120	0.020	0.020 0.940
			α1 1.609 -0.310	0.230	0.230	0.950 1.609 -0.030	0.430	0.450 0.958
		1000 λ1 0.751 9.200	16.630	17.220 0.924 0.749 7.170	21.000	1.610 0.914
			α2 0.693 0.200	8.520	8.840	0.940 0.696 2.890	10.090 10.680 0.942
		5000 λ2 0.742 0.460	2.980	2.690	0.952 0.744 2.080	3.840	3.690 0.940
			α2 0.693 0.550	1.660	1.54	0.950 0.693 0.140	1.980	1.930 0.942
	0.50 1000 λ1 0.182 -0.090	0.060	0.050	0.956 0.182 -0.440	0.100	0.100 0.964
			α1 1.608 -0.980	1.030	1.090	0.952 1.612 2.580	2.020	2.060 0.968
		5000 λ1 0.182 -0.180	0.010	0.010	0.936 0.182 -0.180	0.020	0.020 0.940
			α1 1.609 -0.190	0.200	0.210	0.952 1.610 0.240	0.400	0.410 0.948
		1000 λ2 0.750 9.020	17.470	18.370 0.924 0.748 6.530	21.500 22.090 0.912
			α2 0.693 0.390	8.850	9.390	0.928 0.696 3.120	10.350 11.080 0.930
		5000 λ 2 0.742 0.070	3.170	2.870	0.948 0.744 2.400	3.970	3.760 0.944
			α 2 0.693 0.460	1.730	1.620	0.950 0.693 -0.080	2.040	1.960 0.948
	0.75 1000 λ1 0.182 0.130	0.040	0.040	0.944 0.182 -0.050	0.080	0.080 0.937
			α1 1.609 0.030	7e-04	0.860	0.924 1.610 1.120	1.410	1.450 0.947
		5000 λ1 0.182 -0.120	0.010	0.010	0.940 0.182 -0.040	0.020	0.020 0.948
			α1 1.609 0.040	0.140	0.160	0.936 1.610 0.780	0.270	0.320 0.926
		1000 λ1 0.747 4.900	13.000	14.790 0.924 0.744 1.780	16.030 17.270 0.928
			α1 0.694 1.460	8.370	9.630	0.924 0.697 4.230	9.660	10.980 0.932
		5000 λ1 0.743 1.160	2.440	2.120	0.966 0.744 1.920	3.040	2.590 0.966
			α1 0.693 -0.090	1.650	1.490	0.970 0.693 -0.470	1.910	1,700 0.962

Conflict of Interest None declared

Appendix Simulation Results for Gumbel and Clayton Copula Models

We simulated data to mimic the French breast cancer data set for sample sizes; 1000, and 5000 with 500 replications. The latent failure times for Tj ∼ W eibull(αj, λj) with probability density function defined in section 3. The parameters for the Weibull distribution for T1 were λ1 = 0.182 and α1 = 1.609, while those for T2 were λ2 = 0.742 and α2 = 0.693. In the estimation of λ1, α1 for T1, λ2, α2 are assumed known for T2, and vice versa for estimation of λ2 ,α2. Noninformative censoring times were generated from a uniform distribution (0, γ), where γ was chosen for 10, 30 and 50% censoring. We consider the Clayton copula with Kendall's tau, τ ken = θ θ+2 = 0, 0.25, 0.50, 0.75. Initial parameter values were randomly chosen from uniform distributions, with multiple starting values as described in section 3. The simulation results based on the Clayton copula are presented in the table 4 below.